c program recfunk c 9/26/00 JJP c multitaper program to invert for receiver functions c SINGLE STATION WITH SPECTRAL WEIGHTING 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., in press 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 and time intervals c the data must be binary SAC format c horizontals must be rotated to radial and transverse 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, 599 events max c c code does NOT combine data with different sample rates c data files limited to 15K 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 c c these files are overwritten everytime you run the program c so rename them to save them c c to compile recfunk: c c xf77 -o recfunk recfunk.f plotlib c c The "-fast -native" option will generate a faster executable. Note that c c alias xf77 'f77 -I${OPENWINHOME}/include -L${OPENWINHOME}/lib -lxview -lolgx -lX11' c c plotlib is PLOTIT plotting package, used to plot results thruout the code c PLOTIT resides at love.geology.yale.edu FTP site at pub/park/Plotxy c c EISPACK routines TINVIT and TRIDIB compute Slepian tapers for spectrum analysis c implicit real*4(a-h,o-z) real*8 ar,ai,el(12),aaa complex*8 zc,zero complex*8 afft,rf,rfs character*80 name,subtext,string character*28 title(2) character*10 cmps(3) character comp(3) common/nnaamm/iwflag common/npiprol/anpi common/stap2/tap(4100,16) common/taperz/ta(4100,8,3),tai(4100,8,3) common/stap0/ar(4100),ai(4100),bb(4100),time(2000) common/specs/afft(4100,8,3),rf(4100,2),sss(4100,6),crf(4100,2) common/saveit/rft(4100,599,2),crft(4100,2),rfs(2050,599,2), x drfs(2050,599,2),bazs(599),epic(599) common/datastuff/a(15000,3),pdat(4100),tdat(4100),drf(2050,2) common/header/ahead(158) common/distrib/xi(72,36),yi(72,36),sumij(72,36) 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. zero=cmplx(0.,0.) fmax=5. type *,'enter fmax (Hz)' read(5,*) fmax c generate tapers for npad=2048 point records anpi=2.5 nwin=3 itt=4096 npad=4096 nnf=npad/2 do i=1,4000 crft(i,1)=0. crft(i,2)=0. end do irec=0 open(10,file='in_recfunk',form='formatted') 10 type *,'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 type *,name type *,'no filename?' stop 98 continue do kcmp=1,3 name(namlen:namlen)=comp(kcmp) type *,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.15000) then type *,'Careful! data length greater than 15000:', nscan pause endif end do c naturally, you are in deep poop if dt is not the same for all records npre=10./dt npost=30./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 irec=irec+1 baz=ahead(53) bazs(irec)=ahead(53) epic(irec)=ahead(54) c do i=0,79,10 c type *,(ahead(ii+i),ii=1,10) c end do type *,baz,'= back azimuth ',epic(irec),'= epicentral dist' 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 fr(1)=0. fr(2)=ddf fr(3)=ddf C pre and post-window durations are the same c program pads with zeroes if start time is too early type *,'P window start time, duration (sec)' read(10,*) pstart,postwin npts=postwin/dt nst=pstart/dt call taper(npts,nwin,el) type *,'eigenvalues ',(el(k),k=1,nwin) c loop over components type *,'npad,nf,nst,npts',npad,nf,nst,npts type *,'dt,nscan,ddf',dt,nscan,ddf c pause npts0=min0(nst+1,npts) if(npts0.lt.npts) then call taper(npts0,nwin,el) type *,'eigenvalues ',(el(k),k=1,nwin) endif do kcmp=1,3 c first: compute spectrum of pre-event noise do i=1,nf sss(i,kcmp)=0. end do do k=1,nwin do i=1,npts0 c OK to reverse order of data here, as we only use mod-sq spectrum ar(i)=a(nst+1-i,kcmp)*tap(i,k) ai(i)=0.d0 end do do i=npts0+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 type *,'we use the complex FFT routine fft2' call fft2(ar,ai,npad) do i=1,nf sss(i,kcmp)=sss(i,kcmp)+(ar(i)**2+ai(i)**2) end do end do c if pre-event noise duration is truncated, boost the noise spectrum by c relative length factor if(npts0.lt.npts) then do i=1,nf sss(i,kcmp)=sss(i,kcmp)*float(npts)/float(npts0) end do endif end do if(npts0.lt.npts) then call taper(npts,nwin,el) type *,'eigenvalues ',(el(k),k=1,nwin) endif do kcmp=1,3 do i=1,nf sss(i,kcmp+3)=0. end do c first: compute spectrum of P coda do k=1,nwin do i=1,npts ar(i)=a(i+nst,kcmp)*tap(i,k) bb(i)=a(i+nst,kcmp)*tap(i,k) ai(i)=0.d0 end do c call plotit(tt,bb,dum,npts,'tapered data', c x 'time',' ',1,0,0,0,1) 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 type *,'we use the complex FFT routine fft2' call fft2(ar,ai,npad) do i=1,nf afft(i,k,kcmp)=dcmplx(ar(i),ai(i)) sss(i,kcmp+3)=sss(i,kcmp+3)+(ar(i)**2+ai(i)**2) end do c call plotit(fr,bb,dum,nf,'fft', c x 'frequency (Hz)','H(f)',1,0,0,0,1) end do end do c do kcmp=1,3 c call plotit_axes(0.,0.,0.,0.) c call plotit(fr,sss(1,kcmp),dum,nf,' ',' ',' ',1,2,0.1,0,0) c call plotit(fr,sss(1,kcmp+3),dum,nf,cmps(kcmp),'Freq(Hz)', c x ' ',1,2,0,0,60+kcmp) c end do c call plotit_axes(0.,0.,0.,0.) c do kcmp=1,3 c call plotit(fr,sss(1,kcmp+3),dum,nf,'P spectrum','Freq(Hz)', c x ' ',1,0,(3-kcmp)*0.05,0,64*(kcmp/3)) c end do c simple correlation coefficient with vertical (kcmp=3) do l=1,2 do n=1,nf zc=zero do k=1,nwin zc=zc+conjg(afft(n,k,3))*afft(n,k,l) end do rf(n,l)=zc/(sss(n,6)+sss(n,3)) crf(n,l)=cabs(zc)**2/(sss(n,6)*sss(n,3+l)) drf(n,l)=sqrt((1.0-crf(n,l))/(crf(n,l)*(nwin-1))) x *cabs(rf(n,l)) crft(n,l)=crft(n,l)+crf(n,l) c type *,n,l,zc,aaa,rf(n,l) rfs(n,irec,l)=rf(n,l) drfs(n,irec,l)=drf(n,l) end do c call plotit_axes(0.,0.,0.,0.) c call plotit(fr,crf(1,l),dum,nf,'Coherence: '//cmps(l), c x ' ',' ',1,0,0,0,64+l) end do sq2=sqrt(2.) do l=1,2 do n=1,nf fac=cos(pih*(n-1)/nf)**2 c fac=1.0 ar(n)=real(rf(n,l))*fac ar(npad+2-n)=real(rf(n,l))*fac ai(n)=-imag(rf(n,l))*fac ai(npad+2-n)=imag(rf(n,l))*fac pdat(n)=drf(n,l)*fac/sq2 tdat(n)=fr(1)+(n-1)*fr(2) end do do n=nf+1,npad/2+2 ar(n)=0.d0 ar(npad+2-n)=0.d0 ai(n)=0.d0 ai(npad+2-n)=0.d0 end do dtdat=(tdat(nf)-tdat(1))/30. tdata=tdat(1)-dtdat tdatb=tdat(nf)+dtdat bbmx=0. do n=1,npad bb(n)=dsqrt(ar(n)**2+ai(n)**2) bbmx=amax1(bbmx,bb(n)) 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) c call plotit_axes(0.,0.,0.,0.) 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)=datan2d(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.) c invert the Fourier transform 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=float(nnf)/(float(npad)*float(nf)) type *,'ddf,npad,nnf,nf,fac',ddf,npad,nnf,nf,fac 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 do i=1,ntot rft(i,irec,l)=baz+50.*bb(i) end do c call plotit(ttt,bb,dum,ntot/2,'magnified version', c x 'time(sec)','H(t)',1,0,0,0,59+3*l) c call plotit(ttt,bb,dum,ntot,title(l), c x 'time(sec)','H(t)',1,0,0,0,60+3*l) end do go to 10 111 continue close(10) open(12,file='outr_baz.grid',form='formatted') open(13,file='outt_baz.grid',form='formatted') do l=1,2 do i=1,nf crft(i,l)=crft(i,l)/irec end do end do type *,irec c call plotit_axes(0.,0.,0.,0.) c call plotit(fr,crft(1,2),dum,nf,'Average Transverse Coherence', c x 'Freq(Hz)','Avg C\\sup{2}(f)',1,0,0,0,21) c call plotit_axes(0.,0.,0.,0.) c call plotit(fr,crft(1,1),dum,nf,'Average Radial Coherence', c x 'Freq(Hz)','Avg C\\sup{2}(f)',1,0,0,0,22) call plotit_axes(tdata,tdatb,0.0,1.0) bb(1)=0. bb(2)=fmax bb(3)=1./3. bb(4)=bb(3) call plotit(bb,bb(3),dum,2,' ',' ',' ',2,0,0,0,0) call plotit(fr,crft(1,2),dum,nf,' ', x ' ',' ',1,0,0.05,0,0) call plotit(fr,crft(1,1),dum,nf,'Average P-Coda Coherence', x 'Freq(Hz)','Avg C\\sup{2}(f) R-solid T-dashed',1,0,0,0,1) c one more time in case you wanted to save the figure c call plotit(fr,crft(1,1),dum,nf,'Aggregate Radial Coherence', c x 'Freq(Hz)','Coherence',1,0,0,0,21) c call plotit(fr,crft(1,2),dum,nf,'Aggregate Transv Coherence', c x 'Freq(Hz)','Coherence',1,0,0,0,22) c call plotit_axes(0.,0.,0.,0.) c c there are ntot points available in the raw c can only plot 500 total traces in plotit iirec=min(irec,250) nss=7./dt npss=18./dt call plotit_axes(0.,0.,0.,0.) do n=1,iirec call plotit(rft(nss,n,1),time(nss),dum,npss, x 'Radial RF Back-Azimuth Section', x 'Back-Azimuth (degrees CW from N)', c x 'time(sec)',2,0,0,0,21*(n/iirec)) x 'time(sec)',2,0,0,0,(n/iirec)) end do call plotit_axes(0.,0.,0.,0.) do n=1,iirec call plotit(rft(nss,n,2),time(nss),dum,npss, x 'Transverse RF Back-Azimuth Section', x 'Back-Azimuth (degrees CW from N)', c x 'time(sec)',2,0,0,0,22*(n/iirec)) x 'time(sec)',2,0,0,0,(n/iirec)) end do c plot distribution of back-azimuth and epicentral distance do j=1,36 do i=1,72 sumij(i,j)=0.0 xi(i,j)=(i-1)*5. yi(i,j)=(j-1)*5. end do end do do baz=0,355,5 ii=(baz+0.00000001)/5+1 do i=1,irec test=abs(baz-bazs(i)) if(test.lt.5..or.test.gt.355.) then jj=epic(i)/5+1 sumij(ii,jj)=sumij(ii,jj)+1.0 sumij(ii,jj+1)=sumij(ii,jj+1)+1.0 endif end do end do nij=72*36 call plotit(xi,yi,sumij,nij,'Data Distribution','Back Azimuth', x 'Epicentral Distance',5,0,0,3,1) do j=1,36 do i=1,72 if(sumij(i,j).lt.0.05) then sumij(i,j)=0. else sumij(i,j)=0.5+alog(sumij(i,j)) endif end do end do call plotit(xi,yi,sumij,nij,'Data Distribution','Back Azimuth', x 'Epicentral Distance',5,0,0,3,1) type *,irec type *,'enter gcarc halfwidth for summing (gc-limit,gc+limit)' type *,'enter an integer! 0 or 180 --> no limit' read(5,*) limit if(limit.eq.0) limit=180 kaz=0 gc1=0. gc2=180. do baz=355,0,-5 c do baz=0,355,5 c do baz=150,265,5 c code to limit the Epicentral distance range for this sum c first find the bin with most seismograms iii=1+(baz+0.0000001)/5. amost=0. jmost=1 do j=1,36 if(sumij(iii,j).gt.amost) then amost=sumij(iii,j) jmost=j endif end do gc1=jmost*5.-limit gc2=jmost*5.+limit c end code to limit the Epicentral distance range for this sum do l=1,2 do n=1,nf rf(n,l)=zero drf(n,l)=0. end do jrec=0 do i=1,irec test=abs(baz-bazs(i)) if(test.lt.5..or.test.gt.355.) then c branch to limit the Epicentral distance range for this sum if(epic(i).ge.gc1.and.epic(i).le.gc2) then jrec=jrec+1 do n=1,nf sig=drfs(n,i,l) rf(n,l)=rf(n,l)+rfs(n,i,l)/sig drf(n,l)=drf(n,l)+1.0/sig end do endif endif end do print *,baz,jrec if(jrec.gt.0) then 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 fac=cos(pih*(n-1)/nf)**2 rf(n,l)=rf(n,l)/drf(n,l) 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*sqrt(float(jrec))/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 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)=datan2d(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)) 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 write(iunit,1022) t3,baz,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,60+3*l) write(12,101) '>' write(13,101) '>' endif end do call plotit_axes(0.,0.,0.,0.) if(jrec.gt.0) then do i=1,ntot bb(i)=(rft(i,kz,1)-baz)/50. end do call plotit(ttt,bb,dum,npost,string, x 'time(sec)','H(t)',1,0,0.0,0,0) do i=1,ntot bb(i)=(rft(i,kz,2)-baz)/100. end do call plotit(ttt,bb,dum,npost,string, x 'time(sec)','H(t)',1,0,0.1,0,1) endif end do kaz=kaz/2 type *,kaz,' traces' do iagain=1,2 call plotit_axes(-10.,370.,3.,-16.) c call plotit_axes(140.,275.,3.,-16.) 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,21*(n/kaz)) end do call plotit_axes(-10.,370.,3.,-16.) c call plotit_axes(140.,275.,3.,-16.) 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,22*(n/kaz)) end do end do call plotit_axes(-10.,370.,3.,-26.) c call plotit_axes(140.,275.,3.,-26.) 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,(n/kaz)) end do call plotit_axes(-10.,370.,3.,-26.) c call plotit_axes(140.,275.,3.,-26.) 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,(n/kaz)) end do close(12) close(13) cc NOW PLOT TRACES THAT ARE BINNED WITH EPICENTRAL DISTANCE type *,'compute RFs binned with epicentral distance' type *,'enter back-azimuth limits ib1,ib2 (integers!)' type *,' 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. type *,'epicentral bins from',ep1,' to',ep2 type *,'with halfbin width',dep type *,' do you want to change this spacing? (1=yes)' read(5,*) ickk if(ickk.eq.1) then type *,'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 type *,'HEY! There are ',nep,' traces!' endif c end code for changing the epicentral bin width do epi=ep2,ep1,-dep do l=1,2 do n=1,nf rf(n,l)=zero drf(n,l)=0. end do jrec=0 do i=1,irec test=abs(epi-epic(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 if(ick1.eq.1) then jrec=jrec+1 do n=1,nf sig=drfs(n,i,l) rf(n,l)=rf(n,l)+rfs(n,i,l)/sig drf(n,l)=drf(n,l)+1.0/sig end do endif endif end do if(jrec.gt.0) then type *,epi,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 fac=cos(pih*(n-1)/nf)**2 rf(n,l)=rf(n,l)/drf(n,l) 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*sqrt(float(jrec))/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 call plotit_axes(tdata,tdatb,-0.3*bbmx,1.03*bbmx) call plotit(tdat,bb,pdat,nf,'Frequency Domain', 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)=datan2d(ai(n),ar(n)) end do call plotit_axes(tdata,tdatb,-180.,180.) call plotit(tdat,bb,pdat,nf,'Frequency Domain', 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)) 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 write(iunit,1022) t3,epi,bb(i) end do call plotit_axes(0.,0.,0.,0.) call plotit(ttt,bb,dum,npost,string, x 'time(sec)','H(t)',1,0,0,0,60+3*l) endif end do call plotit_axes(0.,0.,0.,0.) if(jrec.gt.0) then do i=1,ntot bb(i)=(rft(i,kz,1)-epi)/50. end do call plotit(ttt,bb,dum,npost,string, x 'time(sec)','H(t)',1,0,0.0,0,0) do i=1,ntot bb(i)=(rft(i,kz,2)-epi)/100. end do call plotit(ttt,bb,dum,npost,string, x 'time(sec)','H(t)',1,0,0.1,0,1) write(12,101) '>' write(13,101) '>' endif end do 1022 format(f7.3,f6.1,f7.3) close(12) close(13) kaz=kaz/2 type *,kaz,' traces' do iagain=1,2 call plotit_axes(ep1-5.,ep2+5.,3.,-16.) 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,21*(n/kaz)) end do call plotit_axes(ep1-5.,ep2+5.,3.,-16.) 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,22*(n/kaz)) end do end do call plotit_axes(ep1-5.,ep2+5.,3.,-26.) 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,(n/kaz)) end do call plotit_axes(ep1-5.,ep2+5.,3.,-26.) 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,(n/kaz)) end do 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(65536) ! we use this common block for scratch space common/stap2/ta(4100,16) dimension a(4100,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 type *,'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 type *,ierr,rlb,rlu type *,'eigenvalues for the eigentapers' c type *,(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 type *,'eigenvalues for the eigentapers (small k)' c type *,(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 type *,'eigenvalues for the eigentapers (k near 2n)' c type *,(ell(i),i=1,nwin) do i=1,nwin el(i)=dmax1(ell(i),el(i)) end do c type *,'composite asymptotics for taper eigenvalues' c type *,(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 type *,'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 TINVIT(NM,N,D,E,E2,M,W,IND,Z, X IERR,RV1,RV2,RV3,RV4,RV6) C INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP REAL*8 D(N),E(N),E2(N),W(M),Z(NM,M), X RV1(N),RV2(N),RV3(N),RV4(N),RV6(N) REAL*8 U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER,MACHEP REAL*8 DSQRT,DABS,DFLOAT INTEGER IND(M) C C THIS SUBROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH- C NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL C SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, C USING INVERSE ITERATION. C C ON INPUT: C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX; C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX; C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY; C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E, C WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E. C E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN C THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM C OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST CONTAIN C 0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0D0 C IF THE EIGENVALUES ARE IN DESCENDING ORDER. IF BISECT, C TRIDIB, OR IMTQLV HAS BEEN USED TO FIND THE EIGENVALUES, C THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE; C C M IS THE NUMBER OF SPECIFIED EIGENVALUES; C C W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER; C C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC. C C ON OUTPUT: C C ALL INPUT ARRAYS ARE UNALTERED; C C Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS. C ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO; C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH C EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS; C C RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C :::::::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C MACHEP = 16.0D0**(-13) FOR LONG FORM ARITHMETIC C ON S360 :::::::::: C FOR F_FLOATING DEC FORTRAN C DATA MACHEP/1.1D-16/ C FOR G_FLOATING DEC FORTRAN DATA MACHEP/1.25D-15/ C IERR = 0 IF (M .EQ. 0) GO TO 1001 TAG = 0 ORDER = 1.0D0 - E2(1) Q = 0 C :::::::::: ESTABLISH AND PROCESS NEXT SUBMATRIX :::::::::: 100 P = Q + 1 C DO 120 Q = P, N IF (Q .EQ. N) GO TO 140 IF (E2(Q+1) .EQ. 0.0D0) GO TO 140 120 CONTINUE C :::::::::: FIND VECTORS BY INVERSE ITERATION :::::::::: 140 TAG = TAG + 1 S = 0 C DO 920 R = 1, M IF (IND(R) .NE. TAG) GO TO 920 ITS = 1 X1 = W(R) IF (S .NE. 0) GO TO 510 C :::::::::: CHECK FOR ISOLATED ROOT :::::::::: XU = 1.0D0 IF (P .NE. Q) GO TO 490 RV6(P) = 1.0D0 GO TO 870 490 NORM = DABS(D(P)) IP = P + 1 C DO 500 I = IP, Q 500 NORM = NORM + DABS(D(I)) + DABS(E(I)) C :::::::::: EPS2 IS THE CRITERION FOR GROUPING, C EPS3 REPLACES ZERO PIVOTS AND EQUAL C ROOTS ARE MODIFIED BY EPS3, C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW :::::::::: EPS2 = 1.0D-3 * NORM EPS3 = MACHEP * NORM UK = DBLE(Q-P+1) EPS4 = UK * EPS3 UK = EPS4 / DSQRT(UK) S = P 505 GROUP = 0 GO TO 520 C :::::::::: LOOK FOR CLOSE OR COINCIDENT ROOTS :::::::::: 510 IF (DABS(X1-X0) .GE. EPS2) GO TO 505 GROUP = GROUP + 1 IF (ORDER * (X1 - X0) .LE. 0.0D0) X1 = X0 + ORDER * EPS3 C :::::::::: ELIMINATION WITH INTERCHANGES AND C INITIALIZATION OF VECTOR :::::::::: 520 V = 0.0D0 C DO 580 I = P, Q RV6(I) = UK IF (I .EQ. P) GO TO 560 IF (DABS(E(I)) .LT. DABS(U)) GO TO 540 C :::::::::: WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF C E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY :::::::::: XU = U / E(I) RV4(I) = XU RV1(I-1) = E(I) RV2(I-1) = D(I) - X1 RV3(I-1) = 0.0D0 IF (I .NE. Q) RV3(I-1) = E(I+1) U = V - XU * RV2(I-1) V = -XU * RV3(I-1) GO TO 580 540 XU = E(I) / U RV4(I) = XU RV1(I-1) = U RV2(I-1) = V RV3(I-1) = 0.0D0 560 U = D(I) - X1 - XU * V IF (I .NE. Q) V = E(I+1) 580 CONTINUE C IF (U .EQ. 0.0D0) U = EPS3 RV1(Q) = U RV2(Q) = 0.0D0 RV3(Q) = 0.0D0 C :::::::::: BACK SUBSTITUTION C FOR I=Q STEP -1 UNTIL P DO -- :::::::::: 600 DO 620 II = P, Q I = P + Q - II RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I) V = U U = RV6(I) 620 CONTINUE C :::::::::: ORTHOGONALIZE WITH RESPECT TO PREVIOUS C MEMBERS OF GROUP :::::::::: IF (GROUP .EQ. 0) GO TO 700 J = R C DO 680 JJ = 1, GROUP 630 J = J - 1 IF (IND(J) .NE. TAG) GO TO 630 XU = 0.0D0 C DO 640 I = P, Q 640 XU = XU + RV6(I) * Z(I,J) C DO 660 I = P, Q 660 RV6(I) = RV6(I) - XU * Z(I,J) C 680 CONTINUE C 700 NORM = 0.0D0 C DO 720 I = P, Q 720 NORM = NORM + DABS(RV6(I)) C IF (NORM .GE. 1.0D0) GO TO 840 C :::::::::: FORWARD SUBSTITUTION :::::::::: IF (ITS .EQ. 5) GO TO 830 IF (NORM .NE. 0.0D0) GO TO 740 RV6(S) = EPS4 S = S + 1 IF (S .GT. Q) S = P GO TO 780 740 XU = EPS4 / NORM C DO 760 I = P, Q 760 RV6(I) = RV6(I) * XU C :::::::::: ELIMINATION OPERATIONS ON NEXT VECTOR C ITERATE :::::::::: 780 DO 820 I = IP, Q U = RV6(I) C :::::::::: IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE C WAS PERFORMED EARLIER IN THE C TRIANGULARIZATION PROCESS :::::::::: IF (RV1(I-1) .NE. E(I)) GO TO 800 U = RV6(I-1) RV6(I-1) = RV6(I) 800 RV6(I) = U - RV4(I) * RV6(I-1) 820 CONTINUE C ITS = ITS + 1 GO TO 600 C :::::::::: SET ERROR -- NON-CONVERGED EIGENVECTOR :::::::::: 830 IERR = -R XU = 0.0D0 GO TO 870 C :::::::::: NORMALIZE SO THAT SUM OF SQUARES IS C 1 AND EXPAND TO FULL ORDER :::::::::: 840 U = 0.0D0 C DO 860 I = P, Q 860 U = U + RV6(I)**2 C XU = 1.0D0 / DSQRT(U) C 870 DO 880 I = 1, N 880 Z(I,R) = 0.0D0 C DO 900 I = P, Q 900 Z(I,R) = RV6(I) * XU C X0 = X1 920 CONTINUE C IF (Q .LT. N) GO TO 100 1001 RETURN C :::::::::: LAST CARD OF TINVIT :::::::::: END SUBROUTINE TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,M,W,IND,IERR,RV4,RV5) C INTEGER I,J,K,L,M,N,P,Q,R,S,II,M1,M2,M11,M22,TAG,IERR,ISTURM REAL*8 D(N),E(N),E2(N),W(M),RV4(N),RV5(N) REAL*8 U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,MACHEP REAL*8 DABS,DMAX1,DMIN1,DFLOAT INTEGER IND(M) C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BISECT, C NUM. MATH. 9, 386-393(1967) BY BARTH, MARTIN, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 249-256(1971). C C THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL C SYMMETRIC MATRIX BETWEEN SPECIFIED BOUNDARY INDICES, C USING BISECTION. C C ON INPUT: C C N IS THE ORDER OF THE MATRIX; C C EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED C EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, C IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, C NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE C PRECISION AND THE 1-NORM OF THE SUBMATRIX; C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX; C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY; C C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2(1) IS ARBITRARY; C C M11 SPECIFIES THE LOWER BOUNDARY INDEX FOR THE DESIRED C EIGENVALUES; C C M SPECIFIES THE NUMBER OF EIGENVALUES DESIRED. THE UPPER C BOUNDARY INDEX M22 IS THEN OBTAINED AS M22=M11+M-1. C C ON OUTPUT: C C EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS C (LAST) DEFAULT VALUE; C C D AND E ARE UNALTERED; C C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. C E2(1) IS ALSO SET TO ZERO; C C LB AND UB DEFINE AN INTERVAL CONTAINING EXACTLY THE DESIRED C EIGENVALUES; C C W CONTAINS, IN ITS FIRST M POSITIONS, THE EIGENVALUES C BETWEEN INDICES M11 AND M22 IN ASCENDING ORDER; C C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.; C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C 3*N+1 IF MULTIPLE EIGENVALUES AT INDEX M11 MAKE C UNIQUE SELECTION IMPOSSIBLE, C 3*N+2 IF MULTIPLE EIGENVALUES AT INDEX M22 MAKE C UNIQUE SELECTION IMPOSSIBLE; C C RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. C C NOTE THAT SUBROUTINE TQL1, IMTQL1, OR TQLRAT IS GENERALLY FASTER C THAN TRIDIB, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C :::::::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C MACHEP = 16.0D0**(-13) FOR LONG FORM ARITHMETIC C ON S360 :::::::::: C FOR F_FLOATING DEC FORTRAN C DATA MACHEP/1.1D-16/ C FOR G_FLOATING DEC FORTRAN DATA MACHEP/1.25D-15/ C IERR = 0 TAG = 0 XU = D(1) X0 = D(1) U = 0.0D0 C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ENTRIES AND DETERMINE AN C INTERVAL CONTAINING ALL THE EIGENVALUES :::::::::: DO 40 I = 1, N X1 = U U = 0.0D0 IF (I .NE. N) U = DABS(E(I+1)) XU = DMIN1(D(I)-(X1+U),XU) X0 = DMAX1(D(I)+(X1+U),X0) IF (I .EQ. 1) GO TO 20 IF (DABS(E(I)) .GT. MACHEP * (DABS(D(I)) + DABS(D(I-1)))) X GO TO 40 20 E2(I) = 0.0D0 40 CONTINUE C X1 = DMAX1(DABS(XU),DABS(X0)) * MACHEP * DBLE(N) XU = XU - X1 T1 = XU X0 = X0 + X1 T2 = X0 C :::::::::: DETERMINE AN INTERVAL CONTAINING EXACTLY C THE DESIRED EIGENVALUES :::::::::: P = 1 Q = N M1 = M11 - 1 IF (M1 .EQ. 0) GO TO 75 ISTURM = 1 50 V = X1 X1 = XU + (X0 - XU) * 0.5D0 IF (X1 .EQ. V) GO TO 980 GO TO 320 60 IF (S - M1) 65, 73, 70 65 XU = X1 GO TO 50 70 X0 = X1 GO TO 50 73 XU = X1 T1 = X1 75 M22 = M1 + M IF (M22 .EQ. N) GO TO 90 X0 = T2 ISTURM = 2 GO TO 50 80 IF (S - M22) 65, 85, 70 85 T2 = X1 90 Q = 0 R = 0 C :::::::::: ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING C INTERVAL BY THE GERSCHGORIN BOUNDS :::::::::: 100 IF (R .EQ. M) GO TO 1001 TAG = TAG + 1 P = Q + 1 XU = D(P) X0 = D(P) U = 0.0D0 C DO 120 Q = P, N X1 = U U = 0.0D0 V = 0.0D0 IF (Q .EQ. N) GO TO 110 U = DABS(E(Q+1)) V = E2(Q+1) 110 XU = DMIN1(D(Q)-(X1+U),XU) X0 = DMAX1(D(Q)+(X1+U),X0) IF (V .EQ. 0.0D0) GO TO 140 120 CONTINUE C 140 X1 = DMAX1(DABS(XU),DABS(X0)) * MACHEP IF (EPS1 .LE. 0.0D0) EPS1 = -X1 IF (P .NE. Q) GO TO 180 C :::::::::: CHECK FOR ISOLATED ROOT WITHIN INTERVAL :::::::::: IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940 M1 = P M2 = P RV5(P) = D(P) GO TO 900 180 X1 = X1 * DBLE(Q-P+1) LB = DMAX1(T1,XU-X1) UB = DMIN1(T2,X0+X1) X1 = LB ISTURM = 3 GO TO 320 200 M1 = S + 1 X1 = UB ISTURM = 4 GO TO 320 220 M2 = S IF (M1 .GT. M2) GO TO 940 C :::::::::: FIND ROOTS BY BISECTION :::::::::: X0 = UB ISTURM = 5 C DO 240 I = M1, M2 RV5(I) = UB RV4(I) = LB 240 CONTINUE C :::::::::: LOOP FOR K-TH EIGENVALUE C FOR K=M2 STEP -1 UNTIL M1 DO -- C (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) :::::::::: K = M2 250 XU = LB C :::::::::: FOR I=K STEP -1 UNTIL M1 DO -- :::::::::: DO 260 II = M1, K I = M1 + K - II IF (XU .GE. RV4(I)) GO TO 260 XU = RV4(I) GO TO 280 260 CONTINUE C 280 IF (X0 .GT. RV5(K)) X0 = RV5(K) C :::::::::: NEXT BISECTION STEP :::::::::: 300 X1 = (XU + X0) * 0.5D0 IF ((X0 - XU) .LE. (2.0D0 * MACHEP * X (DABS(XU) + DABS(X0)) + DABS(EPS1))) GO TO 420 C :::::::::: IN-LINE PROCEDURE FOR STURM SEQUENCE :::::::::: 320 S = P - 1 U = 1.0D0 C DO 340 I = P, Q IF (U .NE. 0.0D0) GO TO 325 V = DABS(E(I)) / MACHEP IF (E2(I) .EQ. 0.0D0) V = 0.0D0 GO TO 330 325 V = E2(I) / U 330 U = D(I) - X1 - V IF (U .LT. 0.0D0) S = S + 1 340 CONTINUE C GO TO (60,80,200,220,360), ISTURM C :::::::::: REFINE INTERVALS :::::::::: 360 IF (S .GE. K) GO TO 400 XU = X1 IF (S .GE. M1) GO TO 380 RV4(M1) = X1 GO TO 300 380 RV4(S+1) = X1 IF (RV5(S) .GT. X1) RV5(S) = X1 GO TO 300 400 X0 = X1 GO TO 300 C :::::::::: K-TH EIGENVALUE FOUND :::::::::: 420 RV5(K) = X1 K = K - 1 IF (K .GE. M1) GO TO 250 C :::::::::: ORDER EIGENVALUES TAGGED WITH THEIR C SUBMATRIX ASSOCIATIONS :::::::::: 900 S = R R = R + M2 - M1 + 1 J = 1 K = M1 C DO 920 L = 1, R IF (J .GT. S) GO TO 910 IF (K .GT. M2) GO TO 940 IF (RV5(K) .GE. W(L)) GO TO 915 C DO 905 II = J, S I = L + S - II W(I+1) = W(I) IND(I+1) = IND(I) 905 CONTINUE C 910 W(L) = RV5(K) IND(L) = TAG K = K + 1 GO TO 920 915 J = J + 1 920 CONTINUE C 940 IF (Q .LT. N) GO TO 100 GO TO 1001 C :::::::::: SET ERROR -- INTERVAL CANNOT BE FOUND CONTAINING C EXACTLY THE DESIRED EIGENVALUES :::::::::: 980 IERR = 3 * N + ISTURM 1001 LB = T1 UB = T2 RETURN C :::::::::: LAST CARD OF TRIDIB :::::::::: END subroutine sacread(name,a,ierr) c read a SAC-format file real*4 ahead character*80 name character*4 chead(158) common/header/ahead(158) dimension a(1) dimension iah(158) equivalence (ahead,iah),(ahead,chead) ierr=0 open(8,file=name,access='direct',recl=512,err=999) c read the 158-word header, can use its info read(8,rec=1,err=998)(ahead(i),i=1,128) type *,(ahead(i),i=50,55) type *,(iah(70+i),i=1,10) nscan=iah(80) irec=2 read(8,rec=2,err=998)(ahead(i+128),i=1,30),(a(j),j=1,98) c type *,'rec=2' if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then nword=nscan+158 nrec=(nword-1)/128+1 last=nword-128*(nrec-1) else type *,'reading greens functions' nword=6*nscan+158 nrec=(nword-1)/128+1 last=nword-128*(nrec-1) endif if(nrec.gt.3) then do irec=3,nrec-1 c type *,'rec=',irec read(8,rec=irec,err=998)(a((irec-3)*128+98+i),i=1,128) end do endif c type *,'rec=',nrec read(8,rec=nrec,err=998)(a((nrec-3)*128+98+i),i=1,last) close(8) return 999 type *,'open error' ierr=1 return 998 type *,'read error: reading',irec,' out of',nrec ierr=1 print *,irec close(8) return end