c program rfmigrate_estack c 5/10/02 c adapted from recfunk_estack c 5/11/01 JJP --> REV2 5/16/01 c updated 08/17/07 JJP c c xf77 -o /park/backup/bin/rfmigrate_estack rfmigrate_estack.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a c xf77 -o /Users/jjpark/bin/rfmigrate_estack rfmigrate_estack.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a c c will compute RFs using single eigenspectrum estimate per record c instead of cross-correlating eigenspectra in one event c --> cross-correlate across events c steps: c LOOP1 c read the filenames to analyse, compute 1-pi prolate spectra c and coherence estimates c save baz/epicenter info c LOOP2 c choose criteria for estack c for each group of events: loop over freq: c assemble the spectral vectors, cross-correlate c with coherence weighting c compute the H(f), display c c c BASED ON c program recfunk_noplot c 5/25/99 JJP c multitaper program to invert for receiver functions c SINGLE STATION WITH SPECTRAL WEIGHTING c c "Multiple Taper Correlation" (MTC) Receiver functions are described in c c Park, J., and V. Levin, Receiver functions from multiple-taper spectral c correlation estimates, Bull. Seismo. Soc Amer., v90, pp1507-1520, 2000. c c spectral estimates of pre-event noise c weight the inversion in the frequency domain, thereby avoiding the c usual problem with microseism noise in 0.1-0.5 Hz c c the main advantage of MTC is not the noise damping, however -- the c cross-correlation of the horizontal and vertical components is used c in place of spectral division when estimating the RF. c this approach appears to be more stable and gives RF-uncertainties c in the frequency domain that are useful for stacking. c c this version of the RF code reads a file of data filenames c you have two choices: either read the time intervals in the filename file c or read them in the sac header c the data must be binary SAC format c horizontals must be rotated to radial and transverse c c for start times in the file: c the file is "in_recfunk" and has lines of the form: c c 1997.070.19.33.bh? <-- code replaces ? with r,t,z c 57 52 <-- start time of analysis window, window duration (sec) c 1997.076.08.15.bh? c 62 62 c ... c ... c ... c stop <-- tells code that data is finished, 799 events max c c c for start times in the SAC header c reads seismic record start times from the sac header c will search the SAC header for specific markers of P phases c T1 - P, Pdiff ahead(12) c T2 - PKP,PKIKP ahead(13) c T3 - PP ahead(14) c T1=T2=T3=0 ==> use original A-marker ahead(9) c c code does NOT combine data with different sample rates c data files limited to 99K pnts. To increase, see common block /datastuff/ c c many intermediate quantities are plotted with PLOTIT as the code proceeds. c other intermediate quantities can be plotted by uncommenting calls to PLOTIT c c the program computes RFs in the Freq domain for multiple records c then it computes composite RFs that are dependent on either c Back-azimuth or epicentral distance. The standard bin width is 10 degrees c with stacked RFs computed every 5 degrees. You can change this in the code. c RF are not migrated in this process, so the user is prompted for parameters c that narrow the focus of the stacks c e.g. a BAZ-sector over which to compute an epicentral RF profile. c c the PLOTIT graphs of RF sweeps can be improved upon, natch c the code writes the BAZ- and EPICEN-dependent RFs to files c in a format easily digested by GMT (traces are separated by ">" lines) c c filenames: out[rt]_baz.grid out[rt]_epi.grid oum[rt]_baz.grid c c these files are overwritten everytime you run the program c so rename them to save them c implicit real*4(a-h,o-z) real*8 ar,ai,el(12),aaa complex*8 zc,zero complex*8 afft,rf,rff character*80 name,subtext,string,name_stack character*28 title(2) character*10 cmps(3) character comp(3) common/nnaamm/iwflag common/npiprol/anpi common/stap2/tap(16400,18) common/taperz/ta(16400,2,3),tai(16400,2,3) common/stap0/ar(16400),ai(16400),bb(16400),time(16400) common/specs3/sss(4100,10),crf(4100,2) common/specs/afft(4100,799,3),cft(4100,799,3) common/specs2/rft(4100,100,2),crft(4100,2),rf(4100,2),rff(4100,2) common/saveit/bazs(799),gcarc(799),iflag(799),epic(799),bbaz(799) c NOTE THAT MODEL PARAMETERS ARE REAL*4 IN THIS CODE, NOT REAL*8 AS IN ANIREC c ONLY 20 LAYERS OVER HALFSPACE ARE ALLOWED common/model/z(21),dz(21),rho(21),vp(21),vs(21),vp2(21), x vp4(21),vs2(21),vss(21) common/model2/xmu(21),xla(21),xmu2(21),xla2(21),xla4(21) common/modstk/tau(21),itau(21),psd(21),gam(21),xik(21) common/migstk/arm(4100,2,21),aim(4100,2,21),dam(4100,2,21), x asumz(4100,2,21),asumh(4100,2,21) common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2) common/header/ahead(158) common/distrib/xi(72,36),yi(72,36),sumij(72,36) common/tt_p/d_p(200),d_pkp(200),s_p(200),s_pkp(200),n_p,n_pkp dimension iah(158),chead(158),mmonth(12),fr(3),tt(3),ttt(3) equivalence (ahead,iah),(ahead,chead) data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/ data comp/'r','t','z'/ data cmps/'Radial ','Transverse','Vertical '/ data title/'Radial Receiver Function ', x 'Transverse Receiver Function'/ con=180./3.14159265358979 pi=3.14159265358979 pih=pi/2. sq2=sqrt(2.) zero=cmplx(0.,0.) fmax=5. print *,'enter velocity model for stacking', x ' =stack_model' read(5,101) name_stack if(name_stack(1:1).eq.' ') then open(7,file='stack_model',form='formatted') else open(7,file=name_stack,form='formatted') endif read(7,101) subtext print *,subtext read(7,*) nl if(nl.gt.20) then print *,'only twenty layers in stack model are allowed!' stop endif c read in theta,phi in degrees - polar coords of fast axis nlp=nl+1 nlm=nl-1 do i=1,nlp read(7,*) theta,phi c read depth to ith interface, vp (m/sec), pk-to-pk cos(2th) relative P pert c pk-to-pk cos(4th) relative P pert, v_s, pk-to-pk cos(2th) relative S pert c density (kg/m**3) read(7,*) z(i),vp(i),vp2(i),vp4(i),vs(i),vs2(i),rho(i) c recall that we interpret fractional values of b,c,e c as peak-to-peak relative velocity perts. c therefore, e=0.02 is 2% pert to mu from slowest to fastest c DO NOT NORMALIZE THE MODEL xmu(i)=rho(i)*vs(i)**2 xmu2(i)=vs2(i)*xmu(i) xla(i)=rho(i)*vp(i)**2 xla2(i)=vp2(i)*xla(i) xla4(i)=vp4(i)*xla(i) vs(i)=vs(i) vp(i)=vp(i) rho(i)=rho(i) z(i)=z(i) end do close(7) do i=2,nl dz(i)=z(i)-z(i-1) end do dz(1)=z(1) call read_ttimes 104 format(2f10.1) print *,'enter fmax (Hz)' read(5,*) fmax c generate tapers for npts-point records, zero-pad to npad=4096 npad=16384 nnf=npad/2 do i=1,4100 crft(i,1)=0. crft(i,2)=0. end do irec=0 1010 print *,'time intervals in file (0) or in SAC-header (1)' read(5,*) isegment if(isegment.eq.0) then print *, x'file of seismic records to read? (space-ret: in_recfunk_estack)' read(5,101) name if(name(1:1).eq.' ') then open(10,file='in_recfunk_estack',form='formatted') else open(10,file=name,form='formatted') endif elseif(isegment.eq.1) then print *, x'file of seismic records to read? (space-ret: in_recpick_estack)' read(5,101) name if(name(1:1).eq.' ') then open(10,file='in_recpick_estack',form='formatted') else open(10,file=name,form='formatted') endif print *,'duration of data windows for analysis' read(5,*) postwin else go to 1010 endif print *,'minimum number of events in a binsum (eg 1 or 2)' read(5,*) ibskip print *,'rotate to LQT coordinates to isolate P and SV? (0=no)' read(5,*) lqt 10 print *,'enter filename e.g. 1998.361.bhz' read(10,101) name 101 format(a) if(name(1:4).eq.'stop') go to 111 c we assume that the last character is either 'r','t',or 'z' do i=80,1,-1 if(name(i:i).ne.' ') then namlen=i go to 98 endif end do print *,name print *,'no filename?' stop 98 continue do kcmp=1,3 name(namlen:namlen)=comp(kcmp) print *,name call sacread(name,a(1,kcmp),ierr) dt=ahead(1) tt(1)=0. tt(2)=dt tt(3)=dt nscan=iah(80) if(nscan.gt.99000) then print *,'Careful! data length greater than 99000:', nscan pause endif end do if(isegment.eq.0) then print *,'P window start time, duration (sec)' read(10,*) pstart,postwin imark=0 else c READING TIME INTERVAL FROM SAC HEADER c ahead(12) is SAC's T1 parameter, for arrival time. Set by pressing T & 1 in plotpk if(ahead(12).gt.1) then tmark=ahead(12) imark=1 c ahead(13) is SAC's T2 parameter, for arrival time. Set by pressing T & 2 in plotpk elseif(ahead(13).gt.1) then tmark=ahead(13) imark=2 c ahead(14) is SAC's T3 parameter, for arrival time. Set by pressing T & 3 in plotpk elseif(ahead(14).gt.1) then tmark=ahead(14) imark=3 c ahead(9) is SAC's A parameter, for arrival time. Set by pressing A in plotpk elseif(ahead(9).gt.1) then tmark=ahead(9) imark=0 else print *,' hey! no P-wave mark for file ',name stop endif pstart=tmark-ahead(6)-3. print *,'sach params b and tmark, pstart',ahead(6),tmark,pstart endif c naturally, you are in deep poop if dt is not the same for all records npre=10./dt npost=85./dt ntot=npre+npost ttt(1)=-npre*dt ttt(2)=dt ttt(3)=dt do i=1,ntot time(i)=-ttt(1)-(i-1)*ttt(2) end do c here we set a trap for deep earthquakes - if too deep, skip the record if(ahead(39).gt.50000) go to 10 irec=irec+1 if(irec.gt.799) then print *,'too many records! 799 MAX' stop endif baz=ahead(53) bazs(irec)=ahead(53) gcarc(irec)=ahead(54) epic(irec)=ahead(54) c do i=0,79,10 c print *,(ahead(ii+i),ii=1,10) c end do print *,baz,'= back azimuth ',gcarc(irec),'= epicentral dist' if(ahead(53).lt.-0.1.or.ahead(53).gt.360.1) then print *,'HEY! THE BACK AZIMUTH OF THIS RECORD IS OUT OF RANGE!' stop elseif(ahead(54).lt.-0.1.or.ahead(54).gt.180.1) then print *,'HEY! THE GCARC OF THIS RECORD IS OUT OF RANGE!' stop endif c kluge: wire the horizontal comps as linear combinations of z-comp c do i=100,nscan ccc a(i,1)=-0.1*a(i-10,3) c a(i,1)=-0.1*a(i,3)+0.1*a(i-20,3)+0.1*a(i-100,3) c a(i,2)=0.3*a(i,3)+0.1*a(i-80,3)-0.1*a(i-100,3) c end do c end kluge c call plotit_axes(0.,0.,0.,0.) c call manyplot(15000,1,nscan,3,tt,a,1.1,name,'time(sec)',' ') c ddf is cycles/day ddf=1./(dt*npad) nf=fmax/ddf+1 if(nf.gt.4100) then print *,'nf is too big', nf stop endif fr(1)=0. fr(2)=ddf fr(3)=ddf npts=postwin/dt nst=pstart/dt anpi=1. nwin=1 call taper(npts,nwin,el) print *,'eigenvalues ',(el(k),k=1,nwin) c loop over components print *,'npad,nf,nst,npts',npad,nf,nst,npts print *,'dt,nscan,ddf',dt,nscan,ddf epipi=epic(irec) c preferred option -- use parameterized relation between GCARC and zrrot c compute approximate angle between Z and R components for P wave c if the phase is PP, its moveout is the same as a P wave at gcarc/2 iflag(irec)=1 if(imark.eq.3) then epipi=epipi/2. gcarc(irec)=epipi endif if(imark.eq.2) then zrrot=9.0 iflag(irec)=2 else if(epipi.lt.80.) then zrrot=42.+(20.-epipi)*(20./60.) elseif(epipi.lt.118.) then zrrot=16.+(6./900.)*(110.-epipi)**2 elseif(imark.eq.0) then zrrot=9.0 iflag(irec)=2 else zrrot=16.+(6./900.)*(110.-epipi)**2 endif endif c here is the option to rotate the Radial and vertical Components to c a P and SV coordinate system if(lqt.ne.0) then zrrot=zrrot/con cs=cos(zrrot) sn=sin(zrrot) c rotate from Z-R coords to L-Q, similar to P-SV do i=1,nst+npts plog=cs*a(i,3)+sn*a(i,1) a(i,1)=cs*a(i,1)-sn*a(i,3) a(i,3)=plog end do endif do kcmp=1,3 c first: compute spectrum of P coda do k=1,nwin do i=1,npts ar(i)=a(i+nst,kcmp)*tap(i,k) ai(i)=0.d0 end do do i=npts+1,npad ar(i)=0.d0 ai(i)=0.d0 end do c we use the complex FFT routine fft2 -- uses real*8 input arrays c print *,'we use the complex FFT routine fft2' call fft2(ar,ai,npad) do n=1,nf afft(n,irec,kcmp)=dcmplx(ar(n),ai(n)) end do end do end do c compute coherency estimates for regression weighting c to do this without fudge factors, recompute the spectra for 3 npi=2.5 tapers anpi=2.5 nwin=3 call taper(npts,nwin,el) print *,'eigenvalues ',(el(k),k=1,nwin) c loop over components print *,'npad,nf,nst,npts',npad,nf,nst,npts print *,'dt,nscan,ddf',dt,nscan,ddf c first part of calculation: cft accumulates hi-res spectra of 3 comps c rf accumulates unnormalized correlation c 2nd part of calculation: cft(*,*,[12]) has [RT] coherence c cft(*,*,3) has total spectrum for weighting do j=1,3 do n=1,nf cft(n,irec,j)=0. end do end do do j=1,2 do n=1,nf rf(n,j)=zero end do end do c first part: compute spectrum of P coda do k=1,nwin do kcmp=1,3 do i=1,npts ar(i)=a(i+nst,kcmp)*tap(i,k) ai(i)=0.d0 end do do i=npts+1,npad ar(i)=0.d0 ai(i)=0.d0 end do c we use the complex FFT routine fft2 -- uses real*8 input arrays c print *,'we use the complex FFT routine fft2' call fft2(ar,ai,npad) if(kcmp.le.2) then do n=1,nf rff(n,kcmp)=dcmplx(ar(n),ai(n)) cft(n,irec,kcmp)=cft(n,irec,kcmp)+(ar(n)**2+ai(n)**2) end do elseif(kcmp.eq.3) then do n=1,nf zc=dcmplx(ar(n),ai(n)) rf(n,1)=rf(n,1)+conjg(zc)*rff(n,1) rf(n,2)=rf(n,2)+conjg(zc)*rff(n,2) cft(n,irec,kcmp)=cft(n,irec,kcmp)+(ar(n)**2+ai(n)**2) end do endif end do end do c second part: assemble quantities for later regression weighting do n=1,nf cohrad=cabs(rf(n,1))**2/(cft(n,irec,1)*cft(n,irec,3)) cohtrn=cabs(rf(n,2))**2/(cft(n,irec,2)*cft(n,irec,3)) sum=0. do kcmp=1,3 sum=sum+cft(n,irec,kcmp) end do cft(n,irec,1)=cohrad cft(n,irec,2)=cohtrn cft(n,irec,3)=sum end do go to 10 111 continue close(10) c plot the inferred angle between radial and vertical against epicentral dist c do i=1,irec c if(epic(i).lt.80.) then c bb(i)=42.+(20.-epic(i))*(20./60.) c elseif(epic(i).lt.118.) then c bb(i)=16.+(6./900.)*(110.-epic(i))**2 c else c bb(i)=9.0 c endif c end do c call plotit(epic,bb,dum,irec,' ',' ',' ',2,0,0.0,0,0) c call plotit(epic,bbaz,dum,irec,'P-SV Rotation Check', c x 'Epicentral Distance','Inferred R/Z Angle',2,0,0.1,3,1) open(12,file='outr_baz.grid',form='formatted') open(13,file='outt_baz.grid',form='formatted') print *,irec,' events' 1009 print *,'enter gcarc interval for summing (GCMIN,GCMAX)' print *,'enter integers between 0 and 180' read(5,*) igcmin,igcmax if(igcmin.ge.igcmax) go to 1009 print *,'bin interval in degrees (e.g. 10)' read(5,*) dbaz kaz=0 gc1=igcmin gc2=igcmax dbazh=dbaz/2. do baz=360.-dbazh,0.,-dbaz do l=1,2 ! loop over R,T do n=1,nf ! loop over freq zc=zero sumz=0. sumh=0. jrec=0 ! count records in the binstack do i=1,irec ! loop over seismic events test=abs(baz-bazs(i)) if(test.lt.dbaz.or.test.gt.360.-dbaz) then c branch to limit the Epicentral distance range for this sum if(gcarc(i).ge.gc1.and.gcarc(i).le.gc2) then jrec=jrec+1 dcoh=sqrt((1.0-cft(n,i,l))/(2.*cft(n,i,l)))*cft(n,i,3) zc=zc+conjg(afft(n,i,3))*afft(n,i,l)/dcoh sumz=sumz+cabs(afft(n,i,3))**2/dcoh sumh=sumh+cabs(afft(n,i,l))**2/dcoh endif endif end do rf(n,l)=zc/sumz crft(n,l)=cabs(zc)**2/(sumz*sumh) if(jrec.gt.1) then drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*(jrec-1))) x *cabs(rf(n,l)) else ! KLUGE TO AVOID ZERODIVIDE drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*0.5)) x *cabs(rf(n,l)) endif end do if(jrec.ge.ibskip) then print *,baz,jrec kaz=kaz+1 c do the expectation of variance for the weighted sum yourself c if you doubt the formula here -- basically the weighted terms in rf-sum c have unit variance, so that variance of the total variance is jrec/drf**2 do n=1,nf c if(n.ge.nf/2) then c fac=cos(pi*(n-nf/2)/nf)**2 c else c fac=1.0 c endif fac=cos(pi*(n-1)/(2*nf))**2 ar(n)=fac*real(rf(n,l)) ai(n)=-fac*imag(rf(n,l)) ar(npad+2-n)=fac*real(rf(n,l)) ai(npad+2-n)=fac*imag(rf(n,l)) drf(n,l)=fac*drf(n,l) end do do n=nf+1,npad-nf+1 ar(n)=0.d0 ai(n)=0.d0 end do write(string,109) baz bbmx=0. do n=1,nf bb(n)=dsqrt(ar(n)**2+ai(n)**2) bbmx=amax1(bbmx,bb(n)) pdat(n)=drf(n,l)/sq2 end do do n=1,npad if(bb(n).gt.1e-4) then pdat(n)=con*(pdat(n)/bb(n)) else pdat(n)=360. endif bb(n)=con*datan2(ai(n),ar(n)) end do do n=1,nf bb(n)=ar(n) end do call fft2(ar,ai,npad) c normalization factor: c divide by npad for fft2 routine normalization c mult by nnf/nf to compensate for losing high freqs c mult by 2 if cosine**2 taper is applied to spectral RF fac=2.*float(nnf)/(float(npad)*float(nf)) c fac=(4./3.)*float(nnf)/(float(npad)*float(nf)) do i=1,npost bb(npre+i)=ar(i)*fac end do do i=1,npre bb(npre+1-i)=ar(npad+1-i)*fac end do kz=(kaz-l)/2+1 do i=1,ntot rft(i,kz,l)=baz+50.*bb(i)*l t3=-time(i) iunit=11+l if(2*((i-1)/2).eq.i-1) write(iunit,1022) t3,baz,bb(i) c write(iunit,1022) t3,baz,bb(i) end do endif end do write(12,101) '>' write(13,101) '>' call plotit_axes(0.,0.,0.,0.) end do kaz=kaz/2 print *,kaz,' traces' nss=7./dt npss=nss+81./dt do iagain=1,2 call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,2),time(nss),dum,npss, x 'Composite Transverse RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,21*(n/kaz)) end do call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,1),time(nss),dum,npss, x 'Composite Radial RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,22*(n/kaz)) end do end do call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,2),time(nss),dum,npost, x 'Composite Transverse RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,(n/kaz)) end do call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,1),time(nss),dum,npost, x 'Composite Radial RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,(n/kaz)) end do close(12) close(13) cc NOW PLOT TRACES THAT ARE BINNED WITH EPICENTRAL DISTANCE print *,'compute RFs binned with epicentral distance' print *,'enter back-azimuth limits ib1,ib2 (integers!)' print *,' ib1=ib2 -> 0,360 and 360-wraparound if ib1 > ib2 ' read(5,*) ib1,ib2 if(ib1.eq.ib2) then ick=1 baz1=0. baz2=360. elseif(ib1.lt.ib2) then ick=1 baz1=ib1 baz2=ib2 else ick=2 baz1=ib1 baz2=ib2 endif kaz=0 open(12,file='outr_epi.grid',form='formatted') open(13,file='outt_epi.grid',form='formatted') c code for changing the epicentral bin width ep1=0. ep2=175. dep=5. print *,'epicentral bins from',ep1,' to',ep2 print *,'with halfbin width',dep print *,' do you want to change this spacing? (1=yes)' read(5,*) ickk if(ickk.eq.1) then print *,'enter delta range and spacing (ep1,ep2,dep)' read(5,*) ep1,ep2,dep nep=(ep2-ep1)/dep+1 if(nep.gt.500.or.nep.lt.1) x print *,'HEY! There are ',nep,' traces!' endif c end code for changing the epicentral bin width do n=1,nf tdat(n)=fr(1)+(n-1)*fr(2) end do do epi=ep2,ep1,-dep do l=1,2 ! loop over R,T do n=1,nf ! loop over freq zc=zero sumz=0. sumh=0. jrec=0 ! count records in the binstack do i=1,irec ! loop over seismic events test=abs(epi-gcarc(i)) if(test.lt.dep) then c branch to limit the back-azimuth range for this sum ick1=0 if(ick.eq.1) then if(bazs(i).ge.baz1.and.bazs(i).le.baz2) ick1=1 else if(bazs(i).ge.baz1.or.bazs(i).le.baz2) ick1=1 endif c normalize the records to enhance highly coherent records c and normalize to equalize the signal size, so large events do not dominate if(ick1.eq.1) then jrec=jrec+1 dcoh=sqrt((1.0-cft(n,i,l))/(2.*cft(n,i,l)))*cft(n,i,3) zc=zc+conjg(afft(n,i,3))*afft(n,i,l)/dcoh sumz=sumz+cabs(afft(n,i,3))**2/dcoh sumh=sumh+cabs(afft(n,i,l))**2/dcoh endif endif end do rf(n,l)=zc/sumz crft(n,l)=cabs(zc)**2/(sumz*sumh) if(jrec.gt.1) then drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*(jrec-1))) x *cabs(rf(n,l)) else ! KLUGE TO AVOID ZERODIVIDE drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*0.5)) x *cabs(rf(n,l)) endif end do if(jrec.ge.ibskip) then print *,epi,jrec kaz=kaz+1 c no weighted sum yet -- just a crude sum with simple error estimate c if you doubt the formula here -- basically the weighted terms in rf-sum c have unit variance, so that variance of the total variance is jrec/drf**2 do n=1,nf c if(n.ge.nf/2) then c fac=cos(pi*(n-nf/2)/nf)**2 c else c fac=1.0 c endif fac=cos(pi*(n-1)/(2*nf))**2 ar(n)=fac*real(rf(n,l)) ai(n)=-fac*imag(rf(n,l)) ar(npad+2-n)=fac*real(rf(n,l)) ai(npad+2-n)=fac*imag(rf(n,l)) drf(n,l)=fac*drf(n,l) end do do n=nf+1,npad-nf+1 ar(n)=0.d0 ai(n)=0.d0 end do write(string,1099) epi bbmx=0. do n=1,nf bb(n)=dsqrt(ar(n)**2+ai(n)**2) bbmx=amax1(bbmx,bb(n)) pdat(n)=drf(n,l)/sq2 end do c call plotit_axes(tdata,tdatb,-0.3*bbmx,1.03*bbmx) c call plotit(tdat,bb,pdat,nf,'Frequency Domain', c x 'frequency (Hz)','|H(f)|',3,0,0.05,0,58+l*3) do n=1,npad if(bb(n).gt.1e-4) then pdat(n)=con*(pdat(n)/bb(n)) else pdat(n)=360. endif bb(n)=con*datan2(ai(n),ar(n)) end do c call plotit_axes(tdata,tdatb,-180.,180.) c call plotit(tdat,bb,pdat,nf,'Frequency Domain', c x 'frequency (Hz)','H(f) phase',3,0,0.05,0,59+l*3) call plotit_axes(0.,0.,0.,0.) do n=1,nf bb(n)=ar(n) end do call plotit_axes(0.,0.,0.,0.) call fft2(ar,ai,npad) c normalization factor: c divide by npad for fft2 routine normalization c mult by nnf/nf to compensate for losing high freqs c mult by 2 if cosine**2 taper is applied to spectral RF fac=2.*float(nnf)/(float(npad)*float(nf)) c fac=(4./3.)*float(nnf)/(float(npad)*float(nf)) do i=1,npost bb(npre+i)=ar(i)*fac end do do i=1,npre bb(npre+1-i)=ar(npad+1-i)*fac end do kz=(kaz-l)/2+1 do i=1,ntot rft(i,kz,l)=epi+50.*bb(i)*l t3=-time(i) iunit=11+l if(2*((i-1)/2).eq.i-1) write(iunit,1022) t3,epi,bb(i) c write(iunit,1022) t3,epi,bb(i) end do c call plotit_axes(0.,0.,0.,0.) c call plotit(ttt,bb,dum,npost,string, c x 'time(sec)','H(t)',1,0,0.0,0,60+3*l) endif end do call plotit_axes(0.,0.,0.,0.) write(12,101) '>' write(13,101) '>' end do 1022 format(f7.3,f8.1,f9.3) close(12) close(13) kaz=kaz/2 print *,kaz,' traces' nss=7./dt npss=nss+81./dt do iagain=1,2 call plotit_axes(ep1-dep,ep2+dep,0.,0.) do n=1,kaz call plotit(rft(nss,n,2),time(nss),dum,npss, x 'Composite Transverse RF Section', x 'Epicentral Distance (degrees)', x 'time(sec)',2,0,0.0,0,21*(n/kaz)) end do call plotit_axes(ep1-dep,ep2+dep,0.,0.) do n=1,kaz call plotit(rft(nss,n,1),time(nss),dum,npss, x 'Composite Radial RF Section', x 'Epicentral Distance (degrees)', x 'time(sec)',2,0,0.0,0,22*(n/kaz)) end do end do call plotit_axes(ep1-dep,ep2+dep,0.,0.) do n=1,kaz call plotit(rft(nss,n,2),time(nss),dum,npost, x 'Composite Transverse RF Section', x 'Epicentral Distance (degrees)', x 'time(sec)',2,0,0.0,0,(n/kaz)) end do call plotit_axes(ep1-dep,ep2+dep,0.,0.) do n=1,kaz call plotit(rft(nss,n,1),time(nss),dum,npost, x 'Composite Radial RF Section', x 'Epicentral Distance (degrees)', x 'time(sec)',2,0,0.0,0,(n/kaz)) end do c NOW TO MIGRATE THE EVENT STACK open(12,file='oumr_baz.grid',form='formatted') open(13,file='oumt_baz.grid',form='formatted') c compute some parameters for the migration scheme c these specify the tie points between migration layers c the lowest layer (nl+1) is taken as evanescent region. sdelay=0. pdelay=0. do i=1,nl sdelay=sdelay+(dz(i)/vs(i)) pdelay=pdelay+(dz(i)/vp(i)) tau(i)=sdelay-pdelay itau(i)=tau(i)/dt print *,itau(i),z(i),tau(i) ztot=ztot+dz(i) end do itau(nlp)=4100 kaz=0 do baz=360.-dbazh,0.,-dbaz c do baz=355.,0.,-5. do l=1,2 ! loop over R,T do k=1,nlp do n=1,nf arm(n,l,k)=0. aim(n,l,k)=0. dam(n,l,k)=0. asumz(n,l,k)=0. asumh(n,l,k)=0. end do end do jrec=0 ! count records in the binstack do i=1,irec ! loop over seismic events test=abs(baz-bazs(i)) if(test.lt.dbaz.or.test.gt.360.-dbaz) then c branch to limit the Epicentral distance range for this sum if(gcarc(i).ge.gc1.and.gcarc(i).le.gc2) then jrec=jrec+1 c TO DO THIS PROPERLY, ONE MUST SPLINE INTERPOLATE THE FREQUENCY AXIS OF C THE EIGENSPECTRA THEN STACK igcarc=gcarc(i) fac=gcarc(i)-igcarc print *,'igcarc,i,gcarc(i),fac',igcarc,i,gcarc(i),fac if(iflag(i).eq.1) then ! P-phase if(igcarc.ge.180.or.igcarc.le.0) then print *,'WHOA! igcarc,gcarc = ',igcarc,gcarc(i) pause endif c allow extreme P_diff phases (moveout is constant for P_diff) if(igcarc.gt.129) igcarc=129 slow=(1.-fac)*s_p(igcarc)+fac*s_p(igcarc+1) elseif(iflag(i).eq.2) then ! PKP-phase if(igcarc.ge.180.or.igcarc.lt.90) then print *,'WHOA! igcarc,gcarc = ',igcarc,gcarc(i) pause endif ig=igcarc-89 slow=(1.-fac)*s_pkp(ig)+fac*s_pkp(ig+1) else print *,'WHOA! iflag(i),i = ',iflag(i),i pause endif c the stack-model velocities are in SI units, s_p is in sec/km slw=(slow/1000.) print *,slow,slw,' SLOWNESS' c compute the slowness-dependent time delays at the interfaces of stack-model do k=1,nlp coss=sqrt(1.-vs(k)*vs(k)*slw*slw) cosp=sqrt(1.-vp(k)*vp(k)*slw*slw) cc trapdoor to avoid incoming slownesses that are evanescent (e.g. head waves) cc basically, one stacks with the last acceptable interval velocity if(vp(k)*slw.lt.1.) then gam(k)=(vp(k)*coss-vs(k)*cosp)/(vp(k)-vs(k)) else gam(k)=gam(k-1) endif end do psd(1)=gam(1)*tau(1) xik(1)=0. do k=2,nlp psd(k)=psd(k-1)+gam(k)*(tau(k)-tau(k-1)) xik(k)=psd(k-1)/gam(k)-tau(k-1) print *,gam(k),psd(k),xik(k) end do c spline-interpolate the cross-correlation in freq domain eps=1.e-10 do n=1,nf dcoh=sqrt((1.0-cft(n,i,l))/(2.*cft(n,i,l)))*cft(n,i,3) zc=conjg(afft(n,i,3))*afft(n,i,l)/dcoh sss(n,1)=real(zc) sss(n,2)=imag(zc) sss(n,3)=dcoh sss(n,4)=alog(eps+cabs(afft(n,i,3))**2/dcoh) sss(n,5)=alog(eps+cabs(afft(n,i,l))**2/dcoh) end do c print *,(itau(k),k=1,nl) c print *,(tau(k),k=1,nl) c print *,(gam(k),k=1,nlp) c print *,(psd(k),k=1,nlp) c print *,nf,slow call splneq(nf,sss(1,1),sss(1,6)) call splneq(nf,sss(1,2),sss(1,7)) call splneq(nf,sss(1,3),sss(1,8)) call splneq(nf,sss(1,4),sss(1,9)) call splneq(nf,sss(1,5),sss(1,10)) c sum into nlp RF stacks, each for a different layer of the stacking model do k=1,nlp dfm=ddf/gam(k) dxm=1.0/gam(k) dcs=cos(2.*pi*dfm*xik(k)) dsn=sin(2.*pi*dfm*xik(k)) c we start at f=0., so that the migration phase factor = 1 c csm and snm are real/imag parts of the migration phase factor c and the position in the spline is x=1 (fm is "migration frequency") csm=1. snm=0. xm=1. do n=1,nf dre=evaleq(xm,nf,sss(1,1),sss(1,6),0,1.) dim=evaleq(xm,nf,sss(1,2),sss(1,7),0,1.) sig=evaleq(xm,nf,sss(1,3),sss(1,8),0,1.) dsumz=evaleq(xm,nf,sss(1,4),sss(1,9),0,1.) dsumh=evaleq(xm,nf,sss(1,5),sss(1,10),0,1.) dsumz=exp(dsumz) dsumh=exp(dsumh) c note -- we do not need to divide by SIG here, unlike in rfmigrate c BECAUSE we have divided by DCOH in defining sss(.,[12]) c DITTO for the summing variables ASUMZ and ASUMH dre=dre/(gam(k)) dim=dim/(gam(k)) arm(n,l,k)=arm(n,l,k)+dre*csm+dim*snm aim(n,l,k)=aim(n,l,k)-dre*snm+dim*csm dam(n,l,k)=dam(n,l,k)+1.0/sig asumz(n,l,k)=asumz(n,l,k)+dsumz asumh(n,l,k)=asumh(n,l,k)+dsumh csm=csm*dcs-snm*dsn snm=snm*dcs+csm*dsn xm=xm+dxm end do c lines are carried over from recfunk_estack, do not apply here c zc=zc+conjg(afft(n,i,3))*afft(n,i,l)/dcoh c sumz=sumz+cabs(afft(n,i,3))**2/dcoh c sumh=sumh+cabs(afft(n,i,l))**2/dcoh end do ! k=1,nlp endif ! epicentral limits endif ! baz binning limits end do ! end the event loop c do k=1,6 c call plotit_axes(0.,0.,0.,0.) c call plotit(tdat,asumz(1,l,k),dum,nf,'summed vert RF spec', c x 'freq (Hz)','amp',1,0,0.05,1,60+k) c end do c do k=1,6 c call plotit_axes(0.,0.,0.,0.) c call plotit(tdat,asumh(1,l,k),dum,nf,'summed horiz RF spec', c x 'freq (Hz)','amp',1,0,0.05,1,60+k) c end do c do k=1,6 c call plotit_axes(0.,0.,0.,0.) c call plotit(tdat,arm(1,l,k),dum,nf,'Real RF spectra', c x 'freq (Hz)','amp',1,0,0.05,1,60+k) c end do c do k=1,6 c call plotit_axes(0.,0.,0.,0.) c call plotit(tdat,aim(1,l,k),dum,nf,'Imag RF spectra', c x 'freq (Hz)','amp',1,0,0.05,1,60+k) c end do end do ! l=1,2 print *,'end the event loop' c pause c we have computed the proper cross-correlations by looping thru the events c now we normalize to obtain the nlp different RFs if(jrec.ge.ibskip) then print *,baz,jrec kaz=kaz+1 do k=1,nlp do l=1,2 c rf(n,l)=zc/sumz c crft(n,l)=cabs(zc)**2/(sumz*sumh) c if(jrec.gt.1) then c drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*(jrec-1))) c x *cabs(rf(n,l)) c else ! KLUGE TO AVOID ZERODIVIDE c drf(n,l)=sqrt((1.0-crft(n,l))/(crft(n,l)*0.5)) c x *cabs(rf(n,l)) c endif c end do c do the expectation of variance for the weighted sum yourself c if you doubt the formula here -- basically the weighted terms in rf-sum c have unit variance, so that variance of the total variance is jrec/drf**2 do n=1,nf c if(n.ge.nf/2) then c fac=cos(pi*(n-nf/2)/nf)**2 c else c fac=1.0 c endif fac=cos(pi*(n-1)/(2*nf))**2 ar(n)=fac*arm(n,l,k)/asumz(n,l,k) ai(n)=-fac*aim(n,l,k)/asumz(n,l,k) ar(npad+2-n)=fac*arm(n,l,k)/asumz(n,l,k) ai(npad+2-n)=fac*aim(n,l,k)/asumz(n,l,k) drf(n,l)=fac*dam(n,l,k) end do do n=nf+1,npad-nf+1 ar(n)=0.d0 ai(n)=0.d0 end do write(string,109) baz bbmx=0. do n=1,nf bb(n)=dsqrt(ar(n)**2+ai(n)**2) bbmx=amax1(bbmx,bb(n)) pdat(n)=drf(n,l)/sq2 end do c call plotit_axes(tdata,tdatb,-0.3*bbmx,1.03*bbmx) c call plotit(tdat,bb,pdat,nf,'Frequency Domain', c x 'frequency (Hz)','|H(f)|',3,0,0.05,0,58+l*3) do n=1,npad if(bb(n).gt.1e-4) then pdat(n)=con*(pdat(n)/bb(n)) else pdat(n)=360. endif bb(n)=con*datan2(ai(n),ar(n)) end do c call plotit_axes(tdata,tdatb,-180.,180.) c call plotit(tdat,bb,pdat,nf,'Frequency Domain', c x 'frequency (Hz)','H(f) phase',3,0,0.05,0,59+l*3) do n=1,nf bb(n)=ar(n) end do call fft2(ar,ai,npad) c normalization factor: c divide by npad for fft2 routine normalization c mult by nnf/nf to compensate for losing high freqs c mult by 2 if cosine**2 taper is applied to spectral RF fac=2.*float(nnf)/(float(npad)*float(nf)) c fac=(4./3.)*float(nnf)/(float(npad)*float(nf)) do i=1,npost bb(npre+i)=ar(i)*fac end do do i=1,npre bb(npre+1-i)=ar(npad+1-i)*fac end do kz=kaz c we store nlp migrated receiver functions in rft(i,kz+k,l), k=1,nlp do i=1,ntot rft(i,kz+k,l)=bb(i) end do c call plotit_axes(0.,0.,0.,0.) c call plotit(ttt,bb,dum,npost,string, c x 'time(sec)','H(t)',1,0,0.0,0,60+3*l) c if(k.le.6) then c call plotit(ttt,bb,dum,npost,string, c x 'time(sec)','H(t)',1,0,(l-1)*0.05,0,(60+k)*(l/2)) c endif end do ! l=1,2 end do ! k=1,nlp print *,'ntot=',ntot do l=1,2 do i=1,ntot rft(i,kz,l)=0. end do i1=1 do k=1,nlp i2=min(ntot,npre+itau(k)) print *,i1,i2 do i=i1,i2 rft(i,kz,l)=rft(i,kz,l)+rft(i,kz+k,l) end do i1=npre+1+itau(k) end do end do do l=1,2 c call plotit_axes(0.,0.,0.,0.) c call manyplot(4100,1,ntot,nlp+1,ttt,rft(1,kz,l),1.1,string, c x 'time(sec)','Layer-Migrated Receiver Functions') do i=1,ntot t3=-time(i) iunit=11+l if(2*((i-1)/2).eq.i-1) x write(iunit,1022) t3,baz,rft(i,kz,l) c write(iunit,1022) t3,baz,rft(i,kz,l) rft(i,kz,l)=baz+50.*rft(i,kz,l)*l end do end do endif write(12,101) '>' write(13,101) '>' end do print *,kaz,' traces' nss=7./dt npss=nss+81./dt do iagain=1,2 call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,2),time(nss),dum,npss, x 'Composite Transverse RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,21*(n/kaz)) end do call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,1),time(nss),dum,npss, x 'Composite Radial RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,22*(n/kaz)) end do end do call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,2),time(nss),dum,npost, x 'Composite Transverse RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,(n/kaz)) end do call plotit_axes(-10.,370.,0.,0.) do n=1,kaz call plotit(rft(nss,n,1),time(nss),dum,npost, x 'Composite Radial RF Section', x 'Back-Azimuth (degrees CW from N)', x 'time(sec)',2,0,0.0,0,(n/kaz)) end do close(12) close(13) 109 format('Back Azimuth Centered on ',f4.0,' degrees') 1099 format('Epicentral Distance Centered on ',f4.0,' degrees') stop end subroutine taper(n,nwin,el) c c to generate slepian tapers c ta is a real*4 array c c j. park c implicit real*4(a-h,o-z) real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb real*8 dfac,drat,gamma,bh,ell common/nnaamm/iwflag common/npiprol/anpi common/tapsum/tapsum(20),ell(20) common/work/ip(16400) common/taperzz/z(270000) ! we use this common block for scratch space common/stap2/ta(16400,18) dimension a(16400,8),el(10) equivalence (a(1,1),ta(1,1)) data iwflag/0/,pi/3.14159265358979d0/ an=dfloat(n) ww=dble(anpi)/an cs=dcos(2.d0*pi*ww) c initialize matrix for eispack subroutine c print *,'ww,cs,an',ww,cs,an do i=0,n-1 ai=dfloat(i) a(i+1,1)=-cs*((an-1.d0)/2.d0-ai)**2 a(i+1,2)=-ai*(an-ai)/2.d0 a(i+1,3)=a(i+1,2)**2 ! required by eispack routine end do eps=1.e-13 m11=1 call tridib(n,eps,a(1,1),a(1,2),a(1,3),rlb,rlu,m11,nwin,el,ip, x ierr,a(1,4),a(1,5)) c print *,ierr,rlb,rlu print *,'eigenvalues for the eigentapers' c print *,(el(i),i=1,nwin) call tinvit(n,n,a(1,1),a(1,2),a(1,3),nwin,el,ip,z,ierr, x a(1,4),a(1,5),a(1,6),a(1,7),a(1,8)) c we calculate the eigenvalues of the dirichlet-kernel problem c i.e. the bandwidth retention factors c from slepian 1978 asymptotic formula, gotten from thomson 1982 eq 2.5 c supplemented by the asymptotic formula for k near 2n from slepian 1978 eq 61 dfac=an*pi*ww drat=8.d0*dfac dfac=4.d0*dsqrt(pi*dfac)*dexp(-2.d0*dfac) do k=1,nwin el(k)=1.d0-dfac dfac=dfac*drat/k ! is this correct formula? yes,but fails as k -> 2n end do c print *,'eigenvalues for the eigentapers (small k)' c print *,(el(i),i=1,nwin) gamma=dlog(8.d0*an*dsin(2.d0*pi*ww))+0.5772156649d0 do k=1,nwin bh=-2.d0*pi*(an*ww-dfloat(k-1)/2.d0-.25d0)/gamma ell(k)=1.d0/(1.d0+dexp(pi*bh)) end do c print *,'eigenvalues for the eigentapers (k near 2n)' c print *,(ell(i),i=1,nwin) do i=1,nwin el(i)=dmax1(ell(i),el(i)) end do c print *,'composite asymptotics for taper eigenvalues' c print *,(el(i),i=1,nwin) c normalize the eigentapers to preserve power for a white process c i.e. they have rms value unity do k=1,nwin kk=(k-1)*n tapsum(k)=0. tapsq=0. do i=1,n aa=z(kk+i) ta(i,k)=aa tapsum(k)=tapsum(k)+aa tapsq=tapsq+aa*aa end do aa=sqrt(tapsq/n) tapsum(k)=tapsum(k)/aa do i=1,n ta(i,k)=ta(i,k)/aa end do end do c print *,'tapsum',(tapsum(i),i=1,nwin) 101 format(80a) c refft preserves amplitudes with zeropadding c zum beispiel: a(i)=1.,i=1,100 will transform at zero freq to b(f=0.)=100 c no matter how much zero padding is done c therefore we need not doctor the taper normalization, c but wait until the fft to force the desired units iwflag=1 return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c routine c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine fft2(ar,ai,n) c fft routine with 2 real input arrays rather than complex - see ffttwo c for comments. c fft2 subroutine mults by exp(i\omega t) c OUTPUT: f=0 in ar(1), f=f_N in ar(n/2+1) c ar(i)=ar(n+2-i), ai(i)=-ai(n+2-i) c fft2 is NOT a unitary transform, mults the series by sqrt(n) c the inverse FT can be effected by running fft2 on the conjugate of c the FFT-expansion, then taking the the conjugate of the output, c and dividing thru by N. to wit: c c assume Xr, Xi is in freq domain, xr, xi in the time domain c c (Xr,Xi)=fft2(xr,xi,N) c (xr,-xi)=fft2(Xr,-Xi,N)/N c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) dimension ar(1),ai(1) mex=dlog(dble(float(n)))/.693147d0 nv2=n/2 nm1=n-1 j=1 do 7 i=1,nm1 if(i .ge. j) go to 5 tr=ar(j) ti=ai(j) ar(j)=ar(i) ai(j)=ai(i) ar(i)=tr ai(i)=ti 5 k=nv2 6 if(k .ge. j) go to 7 j=j-k k=k/2 go to 6 7 j=j+k pi=3.14159265358979d0 do 20 l=1,mex le=2**l le1=le/2 ur=1.0 ui=0. wr=dcos(pi/le1 ) wi=dsin (pi/le1) do 20 j=1,le1 do 10 i=j,n,le ip=i+le1 tr=ar(ip)*ur - ai(ip)*ui ti=ai(ip)*ur + ar(ip)*ui ar(ip)=ar(i)-tr ai(ip)=ai(i) - ti ar(i)=ar(i)+tr ai(i)=ai(i)+ti 10 continue utemp=ur ur=ur*wr - ui*wi ui=ui*wr + utemp*wi 20 continue return end subroutine read_ttimes implicit real*4(a-h,o-z) implicit integer*4(i-n) common/tt_p/d_p(200),d_pkp(200),s_p(200),s_pkp(200),n_p,n_pkp c read in the P- and PKP-slowness values from ttimes c the code interpolates across the upper-mantle triplication for DELTA=19-27 c and chooses P or PKP based on the start time of the data-analysis window c for DELTA.ge.105 do i=1,200 d_p(i)=0. end do c fill in the head-wave slowness for DELTA.lt.14 c 13.2 sec/degree is 8.4+eps km/sec do i=1,13 d_p(i)=float(i) s_p(i)=13.2 end do c open(7,file='/park/backup/Vary/Layers/ttimes_P.dat', c x form='formatted') open(7,file='/Users/jjpark/Research/Layers/ttimes_P.dat', x form='formatted') do i=14,200 read(7,*,end=110) d_p(i),s_p(i) end do 110 close(7) print *,i,' = counter at end of ttimes file' do i=1,200 if(d_p(i).eq.0.) go to 120 end do 120 n_p=i-1 c interpolates the triplication c linearly interpolate slowness between Delta=17 and p=11.06 sec/degree c and Delta=28 and p=8.92 sec/degree do i=18,27 d_p(i)=float(i) s_p(i)=11.06-(i-17)*(11.06-8.92)/11. end do do i=45,n_p d_p(i-18)=d_p(i) s_p(i-18)=s_p(i) end do n_p=n_p-18 do i=1,200 d_pkp(i)=0. end do c open(7,file='/park/backup/Vary/Layers/ttimes_PKP.dat', c x form='formatted') open(7,file='/Users/jjpark/Research/Layers/ttimes_PKP.dat', x form='formatted') do i=1,200 read(7,*,end=130) d_pkp(i),s_pkp(i) end do 130 close(7) print *,i,' = counter at end of ttimes file' do i=1,200 if(d_pkp(i).eq.0.) go to 140 end do 140 n_pkp=i-1 factor=180./(3.14159265358979*6371.) do i=1,n_p s_p(i)=s_p(i)*factor end do do i=1,n_pkp s_pkp(i)=s_pkp(i)*factor end do print *,'n_p,n_pkp',n_p,n_pkp c call plotit_axes(0.,0.,0.,0.) c call plotit(d_p,s_p,dum,n_p,' ',' ',' ',2,0,0.0,0,0) c call plotit(d_pkp,s_pkp,dum,n_pkp,' ', c x 'Epicentral Distance (degrees)','Slowness (sec/km)',2,0,0.0,0,1) return end