dise.f 14.6 KB
Newer Older
1 2 3 4 5 6 7
c this is <dise.f>
c------------------------------------------------------------------------------
c
c 19/04/99 by Thomas Forbriger (IfG Stuttgart)
c
c DIfferential SEismograms
c
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
c ----
c This program is free software; you can redistribute it and/or modify
c it under the terms of the GNU General Public License as published by
c the Free Software Foundation; either version 2 of the License, or
c (at your option) any later version. 
c 
c This program is distributed in the hope that it will be useful,
c but WITHOUT ANY WARRANTY; without even the implied warranty of
c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c GNU General Public License for more details.
c 
c You should have received a copy of the GNU General Public License
c along with this program; if not, write to the Free Software
c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
c ----
c
24 25
c REVISIONS and CHANGES
c    19/04/99   V1.0   Thomas Forbriger
26
c    20/09/00   V1.1   set matching source time
thomas.forbriger's avatar
thomas.forbriger committed
27
c               V1.2   added normalized datasets
28 29 30 31 32 33
c
c==============================================================================
c
      program dise
c
      character*79 version
thomas.forbriger's avatar
thomas.forbriger committed
34
      parameter(version='DISE   V1.2   DIfferential SEismograms')
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
c
c commandline
      integer maxopt, lastarg, iargc
      character*80 argument
      parameter(maxopt=2)
      character*2 optid(maxopt)
      character*40 optarg(maxopt)
      logical optset(maxopt), opthasarg(maxopt)
c 
      character*80 infile1,infile2,outfile,intrace1,intrace2,ssecbeg,ssecend
      integer ifirst,ilast,itrace1,itrace2
      real secbeg,secend
      character*132 wid2line
c 
      integer maxsamp
      parameter(maxsamp=100000)
      real data1(maxsamp), data2(maxsamp), data3(maxsamp)
      integer idata1(maxsamp), idata2(maxsamp), idata3(maxsamp)
      real dataout1(maxsamp), dataout2(maxsamp),dataout3(maxsamp)
thomas.forbriger's avatar
thomas.forbriger committed
54
      real dataout4(maxsamp), dataout5(maxsamp),dataout6(maxsamp)
55
      integer idataout1(maxsamp), idataout2(maxsamp),idataout3(maxsamp)
thomas.forbriger's avatar
thomas.forbriger committed
56
      integer idataout4(maxsamp), idataout5(maxsamp),idataout6(maxsamp)
57 58 59 60 61 62
      equivalence(data1,idata1)
      equivalence(data2,idata2)
      equivalence(data3,idata3)
      equivalence(dataout1,idataout1)
      equivalence(dataout2,idataout2)
      equivalence(dataout3,idataout3)
thomas.forbriger's avatar
thomas.forbriger committed
63 64 65
      equivalence(dataout4,idataout4)
      equivalence(dataout5,idataout5)
      equivalence(dataout6,idataout6)
66 67 68
c 
      integer nsamp1, nsamp2, nout
      real dt1, dt2, tanf1, tanf2
thomas.forbriger's avatar
thomas.forbriger committed
69
      double precision rms1, rms2, rms3, rms4
70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
c 
      integer i, lu, ierr
      parameter(lu=10)
      logical last
c debugging
      logical debug, verbose
c here are the keys to our commandline options
      data optid/2h-d, 2h-v/
      data opthasarg/2*.FALSE./
      data optarg/2*1h-/
c
c------------------------------------------------------------------------------
c basic information
c
c
      argument=' '
      if (iargc().eq.1) call getarg(1, argument)
      if ((argument(1:5).eq.'-help').or.(iargc().ne.7)) then
        print *,version
        print *,'Usage: dise infile1 n1 infile2 n2 outfile tb te'
        print *,'   or: dise -help'
        if (argument(1:5).ne.'-help') stop 'ERROR: wrong number of arguments'
        print *,' '
        print *,'DIfferential SEismograms'
        print *,' '
        print *,'infile1      SFF input file for trace A'
        print *,'n1           trace number of trace A in file ''infile1'' '
        print *,'infile2      SFF input file for trace B'
        print *,'n2           trace number of trace B in file ''infile2'' '
        print *,'outfile      results will be written to this file'
        print *,'tb           time of first sample (in sec)'
        print *,'te           time of last sample (in sec)'
        print *,' '
        print *,'about DISE'
        print *,'=========='
        print *,' '
        print *,'  This program compares two time series. These may be'
        print *,'  real data and synthetics or synthetics calculated with'
        print *,'  two different methods. The two time series are named A'
        print *,'  and B and are read from the input file. In order to'
        print *,'  calculate a difference trace (A-B), trace B will be'
        print *,'  resampled to match the sampling interval and time range'
        print *,'  of trace A. Root mean square (rms) values are calculated'
        print *,'  from all three traces within the specified time range'
        print *,'  (as defined by ''tb'' and ''te''). Ratios of the rms values'
        print *,'  are reported on the terminal.'
thomas.forbriger's avatar
thomas.forbriger committed
116 117
        print *,' '
        print *,'  nA and nB are datasets normalized to an rms of 1.'
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
        stop
      endif
c
c------------------------------------------------------------------------------
c read command line arguments
c
      call tf_cmdline(1, lastarg, maxopt, optid,
     &                optarg, optset, opthasarg)
      debug=optset(1)
      verbose=optset(2)
      call getarg(1,infile1)
      call getarg(2,intrace1)
      call getarg(3,infile2)
      call getarg(4,intrace2)
      call getarg(5,outfile)
      call getarg(6,ssecbeg)
      call getarg(7,ssecend)
c
c------------------------------------------------------------------------------
c go
      read(intrace1, *)itrace1
      read(intrace2, *)itrace2
      read(ssecbeg, *)secbeg
      read(ssecend, *)secend
      call getfile(infile1, itrace1, data1, idata1, maxsamp, nsamp1,
     &  tanf1, dt1)
      call getfile(infile2, itrace2, data2, idata2, maxsamp, nsamp2,
     &  tanf2, dt2)
c 
c do only as much as necessary
      ifirst=1+int((secbeg-tanf1)/dt1)
      ilast=1+int((secend-tanf1)/dt1)
      if ((ifirst.gt.ilast).or.(ifirst.lt.1).or.(ilast.gt.nsamp1))
thomas.forbriger's avatar
thomas.forbriger committed
151
     &  stop 'ERROR: check time range'
152 153 154 155 156 157 158 159
      nout=ilast-ifirst+1
      call resamp(data2, data3, nsamp1, nsamp2, tanf1, dt1, tanf2, dt2,
     &  ifirst, ilast)
c 
      print *,'calculating difference trace (A-B)...'
      rms1=0.
      rms2=0.
      rms3=0.
thomas.forbriger's avatar
thomas.forbriger committed
160
      rms4=0.
161 162 163 164 165 166 167 168 169 170 171
      do i=1,nout
        dataout2(i)=data1(i-1+ifirst)
        dataout3(i)=data3(i)
        dataout1(i)=dataout2(i)-dataout3(i)
        rms1=rms1+dataout1(i)*dataout1(i)
        rms2=rms2+dataout2(i)*dataout2(i)
        rms3=rms3+dataout3(i)*dataout3(i)
      enddo
      rms1=sqrt(rms1/nout)
      rms2=sqrt(rms2/nout)
      rms3=sqrt(rms3/nout)
thomas.forbriger's avatar
thomas.forbriger committed
172 173 174 175 176 177 178 179 180 181 182 183
c
c notice related values:
c   data1   dataout2    rms2
c   data3   dataout3    rms3  (resampled version of data2)
c 
      do i=1,nout
        dataout5(i)=data1(i-1+ifirst)/rms2
        dataout6(i)=data3(i)/rms3
        dataout4(i)=dataout5(i)-dataout6(i)
        rms4=rms4+dataout4(i)*dataout4(i)
      enddo
      rms4=sqrt(rms4/nout)
184 185 186 187
c 
      last=.false.
      print *,'going to write A-B, A and B to ',
     &  outfile(1:index(outfile, ' ')),'...'
188 189
      call sff_WOpenS(lu, outfile, 'nil                           ',
     &                'C', 0., 0., 0., '990420','000000.000', ierr)
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
      if (ierr.ne.0) stop 'ERROR: opening output file'
c 
      call sff_PrepWID2(nout, 1./dt1, 'NSP   ', 1999, 4, 20, 0, 0,
     &  'A-B', 'NSP   ', 'NSP   ', 0., 1., 1., -1., -1., wid2line, ierr)
      if (ierr.ne.0) stop 'ERROR: preparing WID2 line for trace one'
      call sff_WTrace(lu, wid2line, nout, dataout1, idataout1, last, ierr)
      if (ierr.ne.0) stop 'ERROR: writing first trace to output file'
c 
      call sff_PrepWID2(nout, 1./dt1, 'NSP   ', 1999, 4, 20, 0, 0,
     &  'A', 'NSP   ', 'NSP   ', 0., 1., 1., -1., -1., wid2line, ierr)
      if (ierr.ne.0) stop 'ERROR: preparing WID2 line for trace two'
      call sff_WTrace(lu, wid2line, nout, dataout2, idataout2, last, ierr)
      if (ierr.ne.0) stop 'ERROR: writing second trace to output file'
c 
      call sff_PrepWID2(nout, 1./dt1, 'NSP   ', 1999, 4, 20, 0, 0,
     &  'B', 'NSP   ', 'NSP   ', 0., 1., 1., -1., -1., wid2line, ierr)
      if (ierr.ne.0) stop 'ERROR: preparing WID2 line for trace three'
      call sff_WTrace(lu, wid2line, nout, dataout3, idataout3, last, ierr)
      if (ierr.ne.0) stop 'ERROR: writing third trace to output file'
thomas.forbriger's avatar
thomas.forbriger committed
209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
c 
      call sff_PrepWID2(nout, 1./dt1, 'NSP   ', 1999, 4, 20, 0, 0,
     &  'nA-nB', 'NSP   ', 'NSP   ', 0., 1., 1., -1., -1., wid2line, ierr)
      if (ierr.ne.0) stop 'ERROR: preparing WID2 line for trace one'
      call sff_WTrace(lu, wid2line, nout, dataout4, idataout4, last, ierr)
      if (ierr.ne.0) stop 'ERROR: writing fourth trace to output file'
c 
      call sff_PrepWID2(nout, 1./dt1, 'NSP   ', 1999, 4, 20, 0, 0,
     &  'nA', 'NSP   ', 'NSP   ', 0., 1., 1., -1., -1., wid2line, ierr)
      if (ierr.ne.0) stop 'ERROR: preparing WID2 line for trace two'
      call sff_WTrace(lu, wid2line, nout, dataout5, idataout5, last, ierr)
      if (ierr.ne.0) stop 'ERROR: writing fifth trace to output file'
c 
      call sff_PrepWID2(nout, 1./dt1, 'NSP   ', 1999, 4, 20, 0, 0,
     &  'nB', 'NSP   ', 'NSP   ', 0., 1., 1., -1., -1., wid2line, ierr)
      if (ierr.ne.0) stop 'ERROR: preparing WID2 line for trace three'
      last=.true.
      call sff_WTrace(lu, wid2line, nout, dataout6, idataout6, last, ierr)
      if (ierr.ne.0) stop 'ERROR: writing sixth trace to output file'
228 229 230 231 232
c 
      print *,'results are:'
      print *,'  Arms/Brms:     ',rms2,'/',rms3,'=',rms2/rms3
      print *,'  (A-B)rms/Arms: ',rms1,'/',rms2,'=',rms1/rms2
      print *,'  (A-B)rms/Brms: ',rms1,'/',rms3,'=',rms1/rms3
thomas.forbriger's avatar
thomas.forbriger committed
233
      print *,'     (nA-nB)rms: ',rms4
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408
c
      stop
      end
c
c======================================================================
c
      subroutine getfile(filename, itrace, data, idata, maxsamp,
     &                   nsamp, tanf, dt)
c 
      character filename*(*)
      integer itrace, nsamp, maxsamp
      real data(maxsamp)
      integer idata(maxsamp)
      real tanf, dt
c
      character cs*1, code*20, time*20, date*20, tist*20, type*30, rs*1
      character wid2line*132
      real inv
      real c1,c2,c3,r1,r2,r3,sffu_tfirst
      integer ierr, lu, nstack, i
      parameter(lu=10)
      logical last
c 
      print *,'going to read ',filename(1:index(filename, ' ')),'...'
c 
      call sff_ROpenS(lu, filename, inv, tist, code, type, cs, c1, c2, c3,
     &  date,time,ierr)
      if (ierr.ne.0) stop 'ERROR opening input file'
      nsamp=maxsamp
      last=.false.
      do i=1,itrace
        if (last) stop 'ERROR: too few traces in input file'
        call sff_RTraceI(lu, tanf, dt, wid2line, nsamp, data, idata, code,
     &                   last, rs, r1, r2, r3, nstack, ierr)
        if (ierr.ne.0) stop 'ERROR: reading input file'
      enddo
      print *,'  read trace no.',itrace
      if (.not.(last)) close(lu)
      tanf=sffu_tfirst(wid2line, time, date)  
      print *,'  trace begins at ',tanf,'sec'
      print *,'  sampling interval is ',dt,'sec'
      return
      end
c 
c----------------------------------------------------------------------
c
      subroutine resamp(data2, data3, nsamp1, nsamp2, tanf1, dt1, tanf2, dt2,
     &  ifirst, ilast)
c 
      real data2(nsamp2), data3(nsamp1)
      integer nsamp1, nsamp2, ifirst, ilast
      real tanf1, tanf2, dt1, dt2
c 
      logical verbose
      parameter(verbose=.true.)
c
      integer maxspec
      parameter(maxspec=100000)
      complex spect(maxspec)
      integer nspec
c 
      real pi, t, sum,dper,selbeg,selend,oselbeg,oselend
      integer ispec,i,oifirst,oilast,onspec,nout
      integer limit,percount,nextper
      complex ime, expfac, expfacbase
      parameter(pi=3.1415926535897931159979,ime=(0.,1.))
c 
      if (verbose) then
        print *,' '
        print *,'do some resampling...'
      endif
c 
      selbeg=(ifirst-1)*dt1+tanf1
      selend=(ilast-1)*dt1+tanf1
      oifirst=int((selbeg-tanf2)/dt2)
      oilast=2+int((selend-tanf2)/dt2)
      if ((oifirst.gt.oilast).or.(oifirst.lt.1).or.(oilast.gt.nsamp2))
     &  stop 'ERROR: check times'
      onspec=oilast-oifirst+1
      nout=ilast-ifirst+1
      oselbeg=(oifirst-1)*dt2+tanf2
      oselend=(oilast-1)*dt2+tanf2
c 
      print *,'  original trace:'
      print *,'    first sample: ',tanf2,'sec; last sample: ',
     &  (nsamp2-1)*dt2+tanf2,'sec'
      print *,'    sampling interval: ',dt2,'sec'
      print *,'    selected first sample: ',oifirst
      print *,'    selected last sample:  ',oilast
      print *,'    selected start time:   ',oselbeg,'sec'
      print *,'    selected end time:     ',oselend,'sec'
      print *,'  master trace:'
      print *,'    first sample: ',tanf1,'sec; last sample: ',
     &  (nsamp1-1)*dt1+tanf1,'sec'
      print *,'    sampling interval: ',dt1,'sec'
      print *,'    selected first sample: ',ifirst
      print *,'    selected last sample:  ',ilast
      print *,'    selected start time:   ',selbeg,'sec'
      print *,'    selected end time:     ',selend,'sec'
c 
      nspec=2
      do while(nspec.lt.onspec)
        nspec=nspec*2
      enddo
      if (nspec.gt.maxspec) stop 'ERROR: maxspec too small!'
c
      print *,'preparing spectrum with ',nspec,' samples'
      do i=1,nspec
        spect(i)=(0.,0.)
      enddo
c
      print *,'   extracting'
      print *,'     original trace ',oifirst,'-',oifirst-1+onspec
      print *,'     to spectrum trace ',1,'-',onspec
      do i=1,onspec
        spect(i)=cmplx(data2(oifirst+i-1))
      enddo
c 
      call tf_fork(nspec,spect,-1.)
c 
      if (verbose) then
        print *,'  we hold the spectrum now'
        print *,'  going to perform resampling DFT (which will be slow)...'
      endif
c 
c
c dt*sqrt(nsamples) is the factor we have to multiply in order to get
c   spectral coefficients corresponding to a foruier integral transform
c
c dt*(nsamples-1) is the time length of the original trace, which is the 
c   recurrence time of the time limited fourier signal
c 
c 1./(dt*nsamples) is the frequency step size
c
      if (dt1.gt.dt2) then
        print *,'      old sampling interval (',dt2,'sec)'
        print *,'      is less than new one (',dt1,'sec) !'
        print *,
     &  '      it''s up to you to perform a correct low pass filtering'
      endif
c 
c 
      expfacbase=2.*pi*ime/(dt2*float(nspec))
      limit=nspec/2
      dper=10.
      percount=1
      nextper=nout*dper/100.
c 
      do i=1,nout
        t=(i-1)*dt1+selbeg
c        print *,'calc sample ',i,' at ',t,'sec'
        expfac=expfacbase*(t-oselbeg)
        if ((t.lt.oselbeg).or.(t.gt.oselend)) then
          data3(i)=0.
          print *,'trapped'
        else
           sum=0.5*real(spect(1))
           do ispec=2,limit
             sum=sum+real(spect(ispec)*cexp(expfac*float(ispec-1)))
           enddo
           sum=sum+0.5*real(spect(limit+1)*cexp(expfac*(limit)))
           data3(i)=2.*sum/sqrt(float(nspec))
        endif
        if (i.gt.nextper) then
          print *,'    finished ',percount*dper,'%'
          percount=percount+1
          nextper=percount*nout*dper/100.
        endif
      enddo
      print *,'    finished 100.%'
c 
      return
      end
c 
c ----- END OF dise.f -----