Fortran Coder
标题: 读取SAC单道记录 [打印本页]
作者: zhangzhipeng 时间: 2016-2-28 23:58
标题: 读取SAC单道记录
[Fortran] 纯文本查看 复制代码
Module SACParameter
Integer, parameter :: FileUnit = 100
Character( Len = * ), Parameter :: FileName = "..\WAV\WAVFormYiliang\yiliang\AST.GZ.2012251031817.00.BHE"
Type :: EventPara
Integer( kind = 4 ) :: StationNumber = 1
Integer( kind = 4 ) :: FileNumber
Character( Len = 4 ) :: StationName
Character( Len = 4 ) :: FileName
Real(kind = 4), allocatable :: StationWaveformLong(:)
Real(kind = 4), allocatable :: StationMaximumAmplitude(:)
Real(kind = 4), allocatable :: StationDELTA(:)
End Type EventPara
Type :: SACHeader
Real( Kind = 4 ), Allocatable :: DELTA, DEPMIN, DEPMAX, SCALE, ODELTA, B, E, &
O, A, INTERNAL1, T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, F, &
RESP0,RESP1,RESP2,RESP3,RESP4,RESP5,RESP6,RESP7,RESP8,RESP9, &
STLA,STLO,STEL,STDP,EVLA,EVLO,EVEL,EVDP,MAG, &
USER0, USER1, USER2, USER3, USER4, USER5, USER6, USER7, USER8, USER9, &
DIST, AZ, BAZ, GCARC, INTERNAL2, INTERNAL3, &
DEPMEN, CMPAZ, CMPINC, XMINIMUM, XMAXIMUM, YMINIMUM, YMAXIMUM, &
ADJTM, UNUSED1, UNUSED2, UNUSED3, UNUSED4, UNUSED5, UNUSED6
Integer( Kind = 4 ), Allocatable :: NZYEAR, NZJDAY, NZHOUR, NZMIN, NZSEC, NZMSEC, NVHDR, &
NORID, NEVID, NPTS, NSPTS, NWFID, NXSIZE, NYSIZE, UNUSED7, &
IFTYPE, IDEP, IZTYPE, UNUSED8, IINST, ISTREG, IEVREG, IEVTYP, IQUAL, ISYNTH, IMAGTYP, IMAGSRC, &
UNUSED9, UNUSED10, UNUSED11, UNUSED12, UNUSED13, UNUSED14, UNUSED15, UNUSED16, &
LEVEN, LPSPOL, LOVROK, LCALDA, UNUSED17
Character( Len = 16 ), Allocatable :: KEVNM
Character( Len = 8 ), Allocatable :: KSTNM, KHOLE, KO, KA, KT0, KT1, KT2, KT3, KT4, KT5, KT6, KT7, KT8, KT9, &
KF, KUSER0, KUSER1, KUSER2, KCMPNM, KNETWK, KWaveformRD, KINST
Integer(kind = 4),Dimension(70) :: Header1
Real(kind = 4),Dimension(40) :: Header2
Character(Len = 8), Dimension(22) :: Header3
real( kind = 4 ), Allocatable :: Waveform(:)
Integer(kind = 4) :: WaveformLong
End Type SACHeader
End module SACParameter
Program TestReadSAC
USE SACParameter
Implicit none
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
integer :: i, j, k, DELTAdt
real( kind = 4 ), Allocatable :: WaveformOne(:)
Open( Unit = FileUnit, File = 'FileName', Form = 'Unformatted', Status = 'OLD', Action = 'READ' )
Call GetFileN( FileUnit , FileName )
close( FileUnit )
Call RecordSACHeader( SACFILE1%Header1, SACFILE1%Header2, SACFILE1%Header3, SACFILE1%WaveformLong)
Do i = 1, EVENT1%StationNumber
Do j = 1, SACFILE1%WaveformLong
WaveformOne(j) = SACFILE1%Waveform(j)
End do
End do
Write(*,*)"please input the resample delta T!"
read(*,*)DELTAdt
Call Resample(WaveformOne, SACFILE1%WaveformLong, SACFILE1%DELTA, DELTAdt)
Call Toonewaveform( EVENT1%stationnumber, SACFILE1%waveformlong, SACFILE1%waveform, EVENT1%StationWaveformLong, EVENT1%StationMaximumAmplitude)
Call WavePlots(EVENT1%stationnumber, SACFILE1%waveformlong, SACFILE1%waveform, SACFILE1%DIST, SACFILE1%BAZ, SACFILE1%A, dtss, ntss, SACFILE1%KSTNM, nstr_sta, tlong, nnk)
End Program TestReadSAC
!---------------------------------------------
Subroutine GetFileN( iFileUnit )
USE SACParameter
USE, INTRINSIC :: ISO_FORTRAN_ENV
Implicit None
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Logical, Parameter :: FileIsNotEnd = .True.
Integer, Intent( INOUT ) :: iFileUnit
Integer :: i = 1, j = 1, k = 1
!Character( len = 1 ) ::cDummy
!GetFileN = 0
Rewind( iFileUnit )
Do i = 1,70
Read( Unit = iFileUnit, iostat = k, FMT = "(70G15.7)" ) SACFILE1%Header1(i)
!GetFileN = GetFileN + 1
End do
Do i = 1,40
Read( Unit = iFileUnit, iostat = k, FMT = "(40I10)" ) SACFILE1%Header2(i)
!GetFileN = GetFileN + 1
End do
Read( Unit = iFileUnit, iostat = k, FMT = "(A8)" ) SACFILE1%Header3(1)
!GetFileN = GetFileN + 1
Read( Unit = iFileUnit, iostat = k, FMT = "(A16)" ) SACFILE1%KEVNM
!GetFileN = GetFileN + 1
Do i = 2,22
Read( Unit = iFileUnit, iostat = k, FMT = "(21A8)" )SACFILE1%Header3(i)
End do
SACFILE1%WaveformLong = 0
Do while( FileIsNotEnd )
Read( Unit = iFileUnit, iostat = k) SACFILE1%Waveform( SACFILE1%WaveformLong )
!GetFileN = GetFileN + 1
SACFILE1%WaveformLong = SACFILE1%WaveformLong + 1
If( k == IOSTAT_END ) Exit
End do
Rewind( iFileUnit )
End Subroutine GetFileN
!-------------------------------------------
Subroutine RecordSACHeader(Header1, Header2, Header3, WaveformLong)
USE SACPArameter
Implicit None
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Integer(kind = 4), Intent( IN ), Dimension(70) :: Header1
Real(kind = 4), Intent( IN ), Dimension(40) :: Header2
Character(Len = 8), Intent( IN ), Dimension(22) :: Header3
Integer(kind = 4), Intent( IN ) :: WaveformLong
!此处为头段数据的记录，见下一篇帖子
End Subroutine RecordSACHeader
!----------------------------------
Subroutine toonewaveform(stationnumber, waveformlong, waveform, StationWaveformLong, StationMaximumAmplitude)
USE SACParameter
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Real(kind = 4), allocatable :: waveform(:,:)
Real(kind = 4), allocatable :: StationWaveformLong(:)
Real(kind = 4), allocatable :: StationMaximumAmplitude(:)
Integer(kind = 4), Intent( IN ) :: WaveformLong
Integer(kind = 4), Intent( IN ) :: stationnumber
allocate(waveform(waveformlong, stationnumber))
allocate(StationWaveformLong(stationnumber))
allocate(StationMaximumAmplitude(stationnumber))
StationWaveformLong(stationnumber) = WaveformLong
Do i = 1, stationnumber
StationMaximumAmplitude(i) = 1.E-15
Do j = 1, ntss(i)
If (StationMaximumAmplitude(i)<abs(waveform(j,i))) StationMaximumAmplitude(i) = abs(waveform(j,i))
End Do
Do j = 1, StationWaveformLong(i)
waveform(j, i) = waveform(j, i) / StationMaximumAmplitude(i)
End Do
End Do
Return
End Subroutine toonewaveform
!------------------------------------
Subroutine Resample(Waveform, WaveformLong, Delta, DeltaTdt)
!resampling for one dimension signal
USE SACParameter
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Integer(kind = 4), Intent( IN ) :: WaveformLong
real( kind = 4 ), Allocatable :: Waveform(:), WaveformTemp(:)
Integer(kind = 4) :: tLong
!Allocate(WaveformTemp(SACFILE1%WaveformLong))
!Allocate(Waveform(SACFILE1%WaveformLong))
tLong = ( SACFILE1%WaveformLong - 1 ) * SACFILE1%Delta
WaveformLongTemp = int( tLong / Deltadt + 0.5) + 1
Do i = 1, WaveformLongTemp
t = ( i - 1 ) * DELTAdt
Do j = 1, WaveformLong
t1 = ( j - 1 ) * Deltadt
t2 = j * Deltadt
If (t >= t1 .And. t < t2) Then
WaveformTemp( j ) = SACFILE1%Waveform( j ) + ( SACFILE1%Waveform( j + 1 ) - SACFILE1%Waveform( j )) / SACFILE1%Delta * ( t - t1 )
Goto 301
End If !j
End Do
301 Continue
End Do
Do i = 1, WaveformLongTemp
SACFILE1%Waveform(i) = WaveformTemp(i)
End Do
SACFILE1%WaveformLong = WaveformLongTemp
SACFILE1%Delta = DeltaTdt
Return
End Subroutine Resample
!----------------------------------------------
Subroutine BubbleSort(x, ny, n)
Dimension x(n), ny(n)
nsig = 1
m = n - 1
Do While ((m>0) .And. (nsig==1))
nsig = 2.
Do j = 1, m
If (x(j)>x(j+1)) Goto 100
nsig = 1
temp = x(j)
x(j) = x(j+1)
x(j+1) = temp
temp = ny(j)
ny(j) = ny(j+1)
ny(j+1) = temp
100 End Do
m = m - 1
End Do
Return
End Subroutine BubbleSort
作者: zhangzhipeng 时间: 2016-2-29 00:00
标题: SAC头段记录与说明
本帖最后由 zhangzhipeng 于 2016-2-29 09:15 编辑
SAC头段读取并存储。其中字段的文含义随后补上，如有问题欢迎回复[Fortran] 纯文本查看 复制代码
Subroutine RecordSACHeader(Header1, Header2, Header3, WaveformLong)
USE SACPArameter
Implicit None
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Integer(kind = 4), Intent( IN ), Dimension(70) :: Header1
Real(kind = 4), Intent( IN ), Dimension(40) :: Header2
Character(Len = 8), Intent( IN ), Dimension(22) :: Header3
Integer(kind = 4), Intent( IN ) :: WaveformLong
!
SACFILE1%DELTA = SACFILE1%Header1(1) ! 采样间隔
SACFILE1%DEPMIN = SACFILE1%Header1(2) ! 因变量最小值
SACFILE1%DEPMAX = SACFILE1%Header1(3) ! 因变量最大值
SACFILE1%SCALE = SACFILE1%Header1(4) ! 因变量放大系数
SACFILE1%ODELTA = SACFILE1%Header1(5) ! 观测仪器的采样间隔
SACFILE1%B = SACFILE1%Header1(6) ! 自变量的起始值
SACFILE1%E = SACFILE1%Header1(7) ! 自变量的结束值
SACFILE1%O = SACFILE1%Header1(8) ! 发震时刻
SACFILE1%A = SACFILE1%Header1(9) ! 初动到时
SACFILE1%INTERNAL1 = SACFILE1%Header1(10)! 内部变量
SACFILE1%T0 = SACFILE1%Header1(11)! 用户拾取的震相的到时
SACFILE1%T1 = SACFILE1%Header1(12)! 用户拾取的震相的到时
SACFILE1%T2 = SACFILE1%Header1(13)! 用户拾取的震相的到时
SACFILE1%T3 = SACFILE1%Header1(14)! 用户拾取的震相的到时
SACFILE1%T4 = SACFILE1%Header1(15)! 用户拾取的震相的到时
SACFILE1%T5 = SACFILE1%Header1(16)! 用户拾取的震相的到时
SACFILE1%T6 = SACFILE1%Header1(17)! 用户拾取的震相的到时
SACFILE1%T7 = SACFILE1%Header1(18)! 用户拾取的震相的到时
SACFILE1%T8 = SACFILE1%Header1(19)! 用户拾取的震相的到时
SACFILE1%T9 = SACFILE1%Header1(20)! 用户拾取的震相的到时
SACFILE1%F = SACFILE1%Header1(21)! 事件结束时间(相对基准时间的秒数)
SACFILE1%RESP0 = SACFILE1%Header1(22)!
SACFILE1%RESP1 = SACFILE1%Header1(23)!
SACFILE1%RESP2 = SACFILE1%Header1(24)!
SACFILE1%RESP3 = SACFILE1%Header1(25)!
SACFILE1%RESP4 = SACFILE1%Header1(26)!
SACFILE1%RESP5 = SACFILE1%Header1(27)!
SACFILE1%RESP6 = SACFILE1%Header1(28)!
SACFILE1%RESP7 = SACFILE1%Header1(29)!
SACFILE1%RESP8 = SACFILE1%Header1(30)!
SACFILE1%RESP9 = SACFILE1%Header1(31)!
SACFILE1%STLA = SACFILE1%Header1(32)! 台站维度
SACFILE1%STLO = SACFILE1%Header1(33)! 台站经度
SACFILE1%STEL = SACFILE1%Header1(34)! 台站高程
SACFILE1%STDP = SACFILE1%Header1(35)! 台站地表深度
SACFILE1%EVLA = SACFILE1%Header1(36)! 事件维度
SACFILE1%EVLO = SACFILE1%Header1(37)! 事件经度
SACFILE1%EVEL = SACFILE1%Header1(38)! 事件高程
SACFILE1%EVDP = SACFILE1%Header1(39)! 事件地表深度
SACFILE1%MAG = SACFILE1%Header1(40)! 震级
SACFILE1%USER0 = SACFILE1%Header1(41)!
SACFILE1%USER1 = SACFILE1%Header1(42)!
SACFILE1%USER2 = SACFILE1%Header1(43)!
SACFILE1%USER3 = SACFILE1%Header1(44)!
SACFILE1%USER4 = SACFILE1%Header1(45)!
SACFILE1%USER5 = SACFILE1%Header1(46)!
SACFILE1%USER6 = SACFILE1%Header1(47)!
SACFILE1%USER7 = SACFILE1%Header1(48)!
SACFILE1%USER8 = SACFILE1%Header1(49)!
SACFILE1%USER9 = SACFILE1%Header1(50)!
SACFILE1%DIST = SACFILE1%Header1(51)! 震中距
SACFILE1%AZ = SACFILE1%Header1(52)! 方位角
SACFILE1%BAZ = SACFILE1%Header1(53)! 反方位角
SACFILE1%GCARC = SACFILE1%Header1(54)! 台站到事件的大圆弧长度
SACFILE1%INTERNAL2 = SACFILE1%Header1(55)! 内部变量
SACFILE1%INTERNAL3 = SACFILE1%Header1(56)! 内部变量
SACFILE1%DEPMEN = SACFILE1%Header1(57)! 因变量平均值
SACFILE1%CMPAZ = SACFILE1%Header1(58)! 分量方位角
SACFILE1%CMPINC = SACFILE1%Header1(59)! 分量入射角
SACFILE1%XMINIMUM = SACFILE1%Header1(60)!
SACFILE1%XMAXIMUM = SACFILE1%Header1(61)!
SACFILE1%YMINIMUM = SACFILE1%Header1(62)!
SACFILE1%YMAXIMUM = SACFILE1%Header1(63)!
SACFILE1%ADJTM = SACFILE1%Header1(64)!
SACFILE1%UNUSED1 = SACFILE1%Header1(65)! 未使用
SACFILE1%UNUSED2 = SACFILE1%Header1(66)! 未使用
SACFILE1%UNUSED3 = SACFILE1%Header1(67)! 未使用
SACFILE1%UNUSED4 = SACFILE1%Header1(68)! 未使用
SACFILE1%UNUSED5 = SACFILE1%Header1(69)! 未使用
SACFILE1%UNUSED6 = SACFILE1%Header1(70)! 未使用
!
SACFILE1%NZYEAR = SACFILE1%Header2(1) ! GMT年
SACFILE1%NZJDAY = SACFILE1%Header2(2) ! GMT天
SACFILE1%NZHOUR = SACFILE1%Header2(3) ! GMT时
SACFILE1%NZMIN = SACFILE1%Header2(4) ! GMT分
SACFILE1%NZSEC = SACFILE1%Header2(5) ! GMT秒
SACFILE1%NZMSEC = SACFILE1%Header2(6) ! GMT毫秒
SACFILE1%NVHDR = SACFILE1%Header2(7) ! 头段版本号，现在是6
SACFILE1%NORID = SACFILE1%Header2(8) !
SACFILE1%NEVID = SACFILE1%Header2(9) !
SACFILE1%NPTS = SACFILE1%WaveformLong! 记录的长度
SACFILE1%NSPTS = SACFILE1%Header2(11)!
SACFILE1%NWFID = SACFILE1%Header2(12)!
SACFILE1%NXSIZE = SACFILE1%Header2(13)!
SACFILE1%NYSIZE = SACFILE1%Header2(14)!
SACFILE1%UNUSED7 = SACFILE1%Header2(15)!
SACFILE1%IFTYPE = SACFILE1%Header2(16)!
SACFILE1%IDEP = SACFILE1%Header2(17)!
SACFILE1%IZTYPE = SACFILE1%Header2(18)!
SACFILE1%UNUSED8 = SACFILE1%Header2(19)!
SACFILE1%IINST = SACFILE1%Header2(20)!
SACFILE1%ISTREG = SACFILE1%Header2(21)!
SACFILE1%IEVREG = SACFILE1%Header2(22)!
SACFILE1%IEVTYP = SACFILE1%Header2(23)!
SACFILE1%IQUAL = SACFILE1%Header2(24)!
SACFILE1%ISYNTH = SACFILE1%Header2(25)!
SACFILE1%IMAGTYP = SACFILE1%Header2(26)!
SACFILE1%IMAGSRC = SACFILE1%Header2(27)!
SACFILE1%UNUSED9 = SACFILE1%Header2(28)!
SACFILE1%UNUSED10 = SACFILE1%Header2(29)!
SACFILE1%UNUSED11 = SACFILE1%Header2(30)!
SACFILE1%UNUSED12 = SACFILE1%Header2(31)!
SACFILE1%UNUSED13 = SACFILE1%Header2(32)!
SACFILE1%UNUSED14 = SACFILE1%Header2(33)!
SACFILE1%UNUSED15 = SACFILE1%Header2(34)!
SACFILE1%UNUSED16 = SACFILE1%Header2(35)!
SACFILE1%LEVEN = SACFILE1%Header2(36)!
SACFILE1%LPSPOL = SACFILE1%Header2(37)!
SACFILE1%LOVROK = SACFILE1%Header2(38)!
SACFILE1%LCALDA = SACFILE1%Header2(39)!
SACFILE1%UNUSED17 = SACFILE1%Header2(40)!
!
SACFILE1%KSTNM = SACFILE1%Header3(1)
SACFILE1%KHOLE = SACFILE1%Header3(2)
SACFILE1%KO = SACFILE1%Header3(3)
SACFILE1%KA = SACFILE1%Header3(4)
SACFILE1%KT0 = SACFILE1%Header3(5)
SACFILE1%KT1 = SACFILE1%Header3(6)
SACFILE1%KT2 = SACFILE1%Header3(7)
SACFILE1%KT3 = SACFILE1%Header3(8)
SACFILE1%KT4 = SACFILE1%Header3(9)
SACFILE1%KT5 = SACFILE1%Header3(10)
SACFILE1%KT6 = SACFILE1%Header3(11)
SACFILE1%KT7 = SACFILE1%Header3(12)
SACFILE1%KT8 = SACFILE1%Header3(13)
SACFILE1%KT9 = SACFILE1%Header3(14)
SACFILE1%KF = SACFILE1%Header3(15)
SACFILE1%KUSER0 = SACFILE1%Header3(16)
SACFILE1%KUSER1 = SACFILE1%Header3(17)
SACFILE1%KUSER2 = SACFILE1%Header3(18)
SACFILE1%KCMPNM = SACFILE1%Header3(19)
SACFILE1%KNETWK = SACFILE1%Header3(20)
SACFILE1%KWaveformRD = SACFILE1%Header3(21)
SACFILE1%KINST = SACFILE1%Header3(22)
!
End Subroutine RecordSACHeader
作者: zhangzhipeng 时间: 2016-2-29 09:06
标题: 准备写个读写SAC的小工具—这是第一部分
读单道SAC文件，其中头段的记录在下一篇帖子。未完待续。如果运行有问题请回复，不胜感激[Fortran] 纯文本查看 复制代码
Module SACParameter
Integer, parameter :: FileUnit = 100
Character( Len = * ), Parameter :: FileName = "..\WAV\WAVFormYiliang\yiliang\AST.GZ.2012251031817.00.BHE"
Type :: EventPara
Integer( kind = 4 ) :: StationNumber = 1
Integer( kind = 4 ) :: FileNumber
Character( Len = 4 ) :: StationName
Character( Len = 4 ) :: FileName
Real(kind = 4), allocatable :: StationWaveformLong(:)
Real(kind = 4), allocatable :: StationMaximumAmplitude(:)
Real(kind = 4), allocatable :: StationDELTA(:)
End Type EventPara
Type :: SACHeader
Real( Kind = 4 ), Allocatable :: DELTA, DEPMIN, DEPMAX, SCALE, ODELTA, B, E, &
O, A, INTERNAL1, T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, F, &
RESP0,RESP1,RESP2,RESP3,RESP4,RESP5,RESP6,RESP7,RESP8,RESP9, &
STLA,STLO,STEL,STDP,EVLA,EVLO,EVEL,EVDP,MAG, &
USER0, USER1, USER2, USER3, USER4, USER5, USER6, USER7, USER8, USER9, &
DIST, AZ, BAZ, GCARC, INTERNAL2, INTERNAL3, &
DEPMEN, CMPAZ, CMPINC, XMINIMUM, XMAXIMUM, YMINIMUM, YMAXIMUM, &
ADJTM, UNUSED1, UNUSED2, UNUSED3, UNUSED4, UNUSED5, UNUSED6
Integer( Kind = 4 ), Allocatable :: NZYEAR, NZJDAY, NZHOUR, NZMIN, NZSEC, NZMSEC, NVHDR, &
NORID, NEVID, NPTS, NSPTS, NWFID, NXSIZE, NYSIZE, UNUSED7, &
IFTYPE, IDEP, IZTYPE, UNUSED8, IINST, ISTREG, IEVREG, IEVTYP, IQUAL, ISYNTH, IMAGTYP, IMAGSRC, &
UNUSED9, UNUSED10, UNUSED11, UNUSED12, UNUSED13, UNUSED14, UNUSED15, UNUSED16, &
LEVEN, LPSPOL, LOVROK, LCALDA, UNUSED17
Character( Len = 16 ), Allocatable :: KEVNM
Character( Len = 8 ), Allocatable :: KSTNM, KHOLE, KO, KA, KT0, KT1, KT2, KT3, KT4, KT5, KT6, KT7, KT8, KT9, &
KF, KUSER0, KUSER1, KUSER2, KCMPNM, KNETWK, KWaveformRD, KINST
Integer(kind = 4),Dimension(70) :: Header1
Real(kind = 4),Dimension(40) :: Header2
Character(Len = 8), Dimension(22) :: Header3
real( kind = 4 ), Allocatable :: Waveform(:)
Integer(kind = 4) :: WaveformLong
End Type SACHeader
End module SACParameter
Program TestReadSAC
USE SACParameter
Implicit none
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
integer :: i, j, k, DELTAdt
real( kind = 4 ), Allocatable :: WaveformOne(:)
Open( Unit = FileUnit, File = 'FileName', Form = 'Unformatted', Status = 'OLD', Action = 'READ' )
Call GetFileN( FileUnit , FileName )
close( FileUnit )
Call RecordSACHeader( SACFILE1%Header1, SACFILE1%Header2, SACFILE1%Header3, SACFILE1%WaveformLong)
Do i = 1, EVENT1%StationNumber
Do j = 1, SACFILE1%WaveformLong
WaveformOne(j) = SACFILE1%Waveform(j)
End do
End do
Write(*,*)"please input the resample delta T!"
read(*,*)DELTAdt
Call Resample(WaveformOne, SACFILE1%WaveformLong, SACFILE1%DELTA, DELTAdt)
Call Toonewaveform( EVENT1%stationnumber, SACFILE1%waveformlong, SACFILE1%waveform, EVENT1%StationWaveformLong, EVENT1%StationMaximumAmplitude)
Call WavePlots(EVENT1%stationnumber, SACFILE1%waveformlong, SACFILE1%waveform, SACFILE1%DIST, SACFILE1%BAZ, SACFILE1%A, dtss, ntss, SACFILE1%KSTNM, nstr_sta, tlong, nnk)
End Program TestReadSAC
!---------------------------------------------
Subroutine GetFileN( iFileUnit )
USE SACParameter
USE, INTRINSIC :: ISO_FORTRAN_ENV
Implicit None
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Logical, Parameter :: FileIsNotEnd = .True.
Integer, Intent( INOUT ) :: iFileUnit
Integer :: i = 1, j = 1, k = 1
!Character( len = 1 ) ::cDummy
!GetFileN = 0
Rewind( iFileUnit )
Do i = 1,70
Read( Unit = iFileUnit, iostat = k, FMT = "(70G15.7)" ) SACFILE1%Header1(i)
!GetFileN = GetFileN + 1
End do
Do i = 1,40
Read( Unit = iFileUnit, iostat = k, FMT = "(40I10)" ) SACFILE1%Header2(i)
!GetFileN = GetFileN + 1
End do
Read( Unit = iFileUnit, iostat = k, FMT = "(A8)" ) SACFILE1%Header3(1)
!GetFileN = GetFileN + 1
Read( Unit = iFileUnit, iostat = k, FMT = "(A16)" ) SACFILE1%KEVNM
!GetFileN = GetFileN + 1
Do i = 2,22
Read( Unit = iFileUnit, iostat = k, FMT = "(21A8)" )SACFILE1%Header3(i)
End do
SACFILE1%WaveformLong = 0
Do while( FileIsNotEnd )
Read( Unit = iFileUnit, iostat = k) SACFILE1%Waveform( SACFILE1%WaveformLong )
!GetFileN = GetFileN + 1
SACFILE1%WaveformLong = SACFILE1%WaveformLong + 1
If( k == IOSTAT_END ) Exit
End do
Rewind( iFileUnit )
End Subroutine GetFileN
!-------------------------------------------
Subroutine RecordSACHeader(Header1, Header2, Header3, WaveformLong)
USE SACPArameter
Implicit None
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Integer(kind = 4), Intent( IN ), Dimension(70) :: Header1
Real(kind = 4), Intent( IN ), Dimension(40) :: Header2
Character(Len = 8), Intent( IN ), Dimension(22) :: Header3
Integer(kind = 4), Intent( IN ) :: WaveformLong
!此处为头段数据的记录，见下一篇帖子
End Subroutine RecordSACHeader
!----------------------------------
Subroutine toonewaveform(stationnumber, waveformlong, waveform, StationWaveformLong, StationMaximumAmplitude)
USE SACParameter
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Real(kind = 4), allocatable :: waveform(:,:)
Real(kind = 4), allocatable :: StationWaveformLong(:)
Real(kind = 4), allocatable :: StationMaximumAmplitude(:)
Integer(kind = 4), Intent( IN ) :: WaveformLong
Integer(kind = 4), Intent( IN ) :: stationnumber
allocate(waveform(waveformlong, stationnumber))
allocate(StationWaveformLong(stationnumber))
allocate(StationMaximumAmplitude(stationnumber))
StationWaveformLong(stationnumber) = WaveformLong
Do i = 1, stationnumber
StationMaximumAmplitude(i) = 1.E-15
Do j = 1, ntss(i)
If (StationMaximumAmplitude(i)<abs(waveform(j,i))) StationMaximumAmplitude(i) = abs(waveform(j,i))
End Do
Do j = 1, StationWaveformLong(i)
waveform(j, i) = waveform(j, i) / StationMaximumAmplitude(i)
End Do
End Do
Return
End Subroutine toonewaveform
!------------------------------------
Subroutine Resample(Waveform, WaveformLong, Delta, DeltaTdt)
!resampling for one dimension signal
USE SACParameter
type(EventPara) :: EVENT1
type(SACHeader) :: SACFILE1
Integer(kind = 4), Intent( IN ) :: WaveformLong
real( kind = 4 ), Allocatable :: Waveform(:), WaveformTemp(:)
Integer(kind = 4) :: tLong
!Allocate(WaveformTemp(SACFILE1%WaveformLong))
!Allocate(Waveform(SACFILE1%WaveformLong))
tLong = ( SACFILE1%WaveformLong - 1 ) * SACFILE1%Delta
WaveformLongTemp = int( tLong / Deltadt + 0.5) + 1
Do i = 1, WaveformLongTemp
t = ( i - 1 ) * DELTAdt
Do j = 1, WaveformLong
t1 = ( j - 1 ) * Deltadt
t2 = j * Deltadt
If (t >= t1 .And. t < t2) Then
WaveformTemp( j ) = SACFILE1%Waveform( j ) + ( SACFILE1%Waveform( j + 1 ) - SACFILE1%Waveform( j )) / SACFILE1%Delta * ( t - t1 )
Goto 301
End If !j
End Do
301 Continue
End Do
Do i = 1, WaveformLongTemp
SACFILE1%Waveform(i) = WaveformTemp(i)
End Do
SACFILE1%WaveformLong = WaveformLongTemp
SACFILE1%Delta = DeltaTdt
Return
End Subroutine Resample
!----------------------------------------------
Subroutine BubbleSort(x, ny, n)
Dimension x(n), ny(n)
nsig = 1
m = n - 1
Do While ((m>0) .And. (nsig==1))
nsig = 2.
Do j = 1, m
If (x(j)>x(j+1)) Goto 100
nsig = 1
temp = x(j)
x(j) = x(j+1)
x(j+1) = temp
temp = ny(j)
ny(j) = ny(j+1)
ny(j+1) = temp
100 End Do
m = m - 1
End Do
Return
End Subroutine BubbleSort
作者: fcode 时间: 2016-2-29 16:39
SACFILE1%KSTNM = SACFILE1%Header3(1)
SACFILE1%KHOLE = SACFILE1%Header3(2)
SACFILE1%KO = SACFILE1%Header3(3)
SACFILE1%KA = SACFILE1%Header3(4)
SACFILE1%KT0 = SACFILE1%Header3(5)
SACFILE1%KT1 = SACFILE1%Header3(6)
SACFILE1%KT2 = SACFILE1%Header3(7)
SACFILE1%KT3 = SACFILE1%Header3(8)
SACFILE1%KT4 = SACFILE1%Header3(9)
这一堆真是费劲，楼主真有耐心。为何不直接读 SACFILE 结构体？
| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) |
Powered by Discuz! X3.2 |