program autosplit c######################################################################## c c program autosplit c c author - M Evans, (using some code by E Sandvol and N Teanby) c date - July 2003 c c Purpose : c Written to compute shear wave splitting parameters using the c method of Silver & Chan (1991), and the bootstrapping error c analysis of Sandvol and Hearn (1994), minimising either c transverse energy of the smaller eigenvalue of the covariance c matrix. c c Windows length is selected by finding the center point of the c longest period of stable polarization angle and lag time c measurements, (independently for each method). c c Parameters : max - Maximum size of input array (integer) c nshift - Size of gridsearch in time [units of 'dt'] (integer) c nwin - Maximum number of different time windows (integer) c c input : f1 - Filename of north component (character*60) c f2 - Filename of east component (character*60) c nbag - Number of bootstrapping iterations (integer) c method - Measurement method, 1=eigenvalue, 2=transverse, 3=both (integer) c c output : All relevant variables printed to standard output c c internal : count - Array element identifier in main loop (integer) c cov - Covariance matrix (real(2,2)) c doeig - Flag set if eigenvalue method to be performed (integer) c dotrans - Flag set if transverse energy method to be performed (integer) c eigangerr - Bootstrapping angle error, eigenvalue method (real) c eigangwin - Array containing eigenvalue fast polarization direction results (real(nwin)) c eigcur - Identfies array elements containing result of selected window (integer) c eigdelerr - Bootstrapping lag time error, eigenvalue method (real) c eigdelwin - Array containing eigenvalue lag time results (real(nwin)) c eigpolwin - Polarization of first eigenvector (real(nwin)) c i - General loop variable (integer) c lam1 - First eigenvalue (real) c lam2 - Second eigenvalue (real) c nptsmax - Maximum size of waveform (integer) c nptswin - Array identifying window length of each result (integer(nwin)) c traangerr - Bootstrapping angle error, transverse energy method (real) c traangwin - Array containing transverse energy fast polarization direction results (real(nwin)) c tracur - Identfies array elements containing result of selected window (integer) c tradelerr - Bootstrapping lag time error, eigenvalue method (real) c tradelwin(nwin) - Array containing transverse energy lag time results (real(nwin)) c vec1 - First eigenvector (real(2)) c vec2 - Second eigenvector (real(2)) c c common : /eigen/ c eigang - Eigenvalue fast polorization angle result (real) c eigdelay - Eigenvalue lag time result (real) c eigjmin - Number of 'dt' time steps corresponding to 'eigdelay' (integer) c eigminlam1 - First eigenvalue (real) c eigminlam2 - Second eigenvalue (real) c eigpol - Polorization of first eigenvector (real) c minvec1 - First eigenvector (real(2)) c minvec2 - Second eigenvector (real(2)) c c /trans/ c traang - Transverse energy fast polorization angle result (real) c tradelay - Transverse energy lag time result (real) c trajmin - Number of 'dt' time steps corresponding to 'tradelay' (integer) c traminlam1 - First eigenvalue (real) c traminlam2 - Second eigenvector (real(2)) c c /waveforms/ c baz - Back azimuth {station-event azimuth} (real) c caz1 - Component 1 azimuth {North please!} (real) c caz2 - Component 2 azimuth {East please!} (real) c cmp1 - Component 1 time series (real(max)) c cmp2 - Component 2 time series (real(max)) c dt - Increment between samples in 'caz1' and 'caz2' (real) c plt - Flag, set to 1 if output for GMT files is required c c calls subroutines : boot Performs bootstrapping error analysis c covariance Calculates covariance matrix of two time series c eigen2x2 Calculates eigenvalues and eigenvectors of covariance matrix c grid_sch Performs grid search to find best results for both methods c sacin Reads in SAC files c window_select Identifies which window bootstrapping should be performed on c c other subroutines : convolve Convolves two time series c copy Copy one array into another c energy Calculates energy on time series c part_motion Generates output for GMT particle motion plots c pm_polar Calculate direction associated with largest eigenvector c rotate Rotate two seismograms through set angle c rotbaz Rotate seismograms into radial and transverse component c shift Add or remove time delay to seismogram c zero Clear contents of array c c c MODIFICATIONS MADE FOR S&C RESULTS AND REORDERING LOOPS c c################################################################### implicit none integer max, nrot, nshift, nwin parameter(max=600,nrot=181,nshift=121,nwin=600) integer npts character f1*60, f2*60 integer method, nbag c23456789012345678901234567890123456789012345678901234567890123456789012 real cov(2,2), eigangerr, eigangwin(nwin), eigdelwin(nwin) real eigpolwin(nwin), eigdelerr, lam1, lam2, traangerr, tradelerr real traangwin(nwin), tradelwin(nwin), vec1(2), vec2(2) integer count, doeig, dotrans, eigcur, i, nptsmax, nptswin(nwin) integer tracur c *** COMMON BLOCKS *** real eigang, eigdelay, eigminlam1, eigminlam2, eigpol, minvec1 real minvec2 integer eigjmin real eig_grid, tra_grid real traang, tradelay, traminlam1, traminlam2 integer trajmin real data real baz, caz1, caz2, cmp1, cmp2, dt integer plt common /eigen/ eigang, eigdelay, eigjmin, eigminlam1, eigminlam2, & eigpol, minvec1(2), minvec2(2) common /grid/ eig_grid(nrot,nshift), tra_grid(nrot,nshift) common /trans/ traang, tradelay, trajmin, traminlam1, traminlam2 common /values/ data(nrot,nshift,nwin,4) common /waveforms/ baz, caz1, caz2, cmp1(max), cmp2(max), dt, & plt c c ... READ ORIGINAL HORIZONTAL COMPONENTS IN c write(6,*)' ENTER FILENAMES (N FIRST)?' read(5,'(A)')f1 read(5,'(A)')f2 c write(6,*)' ENTER NUMBER OF BOOTSTRAPPING ITERATIONS' read(5,*)nbag c write(6,*)' ENTER METHOD' write(6,*)' (1=EIGENVALUE, 2=TRANSVERSE ENERGY, 3=BOTH)' read(5,*)method c c23456789012345678901234567890123456789012345678901234567890123456789012 write(6,*)' DO YOU WANT GMT OUTPUT' write(6,*)'(If you dont understand the question the answer is NO)' write(6,*)' (1=YES, 0=NO)' read(5,*)plt c call sacin(f1,f2,npts) c c Do measurement c if (method.eq.1.) then doeig = 1. end if if (method.eq.2.) then dotrans = 1. end if if (method.eq.3.) then doeig = 1. dotrans = 1. end if print*,'caz1=',caz1 print*,'caz2=',caz2 print*,'baz=',baz print*,'dt=',dt print*,'npts=',npts nptsmax = npts - (nshift-1) count = 1. call main_calc(method,npts) do 10 i=121,nptsmax npts = i call main_loop(method,npts) if (plt.eq.1.) then write(*,99)i,eigang,eigdelay,traang,tradelay endif if (i.eq.323.) then c write(0,*)'EIG', eigang, eigdelay, eigjmin, eigminlam1, c & eigminlam2, eigpol endif c nptswin(count) = npts traangwin(count) = traang tradelwin(count) = tradelay eigangwin(count) = eigang eigdelwin(count) = eigdelay c lagwin(count) = lag count = count + 1. 10 continue count = count - 1. if (doeig.eq.1.) then call window_select(eigangwin,eigdelwin,count,eigcur) print*,'EIGENVALUE METHOD', nptswin(eigcur) call covariance(cmp1,cmp2,nptswin(eigcur),cov) call eigen2x2(cov,lam1,lam2,vec1,vec2) print*,'orginal lam1',lam1 print*,'orginal lam2',lam2 print*,'eigenvalue original ratio',lam2/lam1 print*,'Window chosen: ',nptswin(eigcur),eigangwin(eigcur), & eigdelwin(eigcur),eigpolwin(eigcur) print*,'1',nptswin(eigcur),nbag,eigangwin(eigcur),eigjmin, & eigpol,eigangerr,eigdelerr c write(0,*)'EIG', eigang, eigdelay, eigjmin, eigminlam1, c & eigminlam2, eigpol c write(0,*)'EIGCUR', eigcur, nptswin(eigcur) call grid_sch(1,nptswin(eigcur)) c write(0,*)'EIG', eigang, eigdelay, eigjmin, eigminlam1, c & eigminlam2, eigpol print*,'POL. DIRECTION',eigpol print*,'fast direction=',eigang print*,'delay=',eigdelay c write(0,*)'CHECK',nptswin(eigcur),nbag,eigangwin(eigcur),eigjmin, c & eigpol,eigangerr,eigdelerr print*,'eigminlam1',eigminlam1 print*,'eigminlam2',eigminlam2 print*,'eigenvalue corrected ratio',eigminlam2/eigminlam1 print*,'minvec1',minvec1 print*,'minvec2',minvec2 print*,'Doing Eigenvalue Bootstrap Error Analysis ...' print*,'Eigenvalue window',float(nptswin(eigcur)*dt) call boot(1,nptswin(eigcur),nbag,eigangwin(eigcur),eigjmin, & eigpol,eigangerr,eigdelerr) end if if (dotrans.eq.1.) then c tracur = 165. c nptswin(tracur) = 285. c traangwin(tracur) = 44. c tradelwin(tracur) = 0.8 call window_select(traangwin,tradelwin,count,tracur) print*,'TRANSVERSE METHOD', nptswin(tracur) call covariance(cmp1,cmp2,nptswin(tracur),cov) call eigen2x2(cov,lam1,lam2,vec1,vec2) print*,'orginal lam1',lam1 print*,'orginal lam2',lam2 print*,'transverse original ratio',lam2/lam1 print*,'Window chosen: ',nptswin(tracur),traangwin(tracur), & tradelwin(tracur) call grid_sch(2,nptswin(tracur)) print*,'fast direction=',traang print*,'delay=',tradelay print*,'traminlam1',traminlam1 print*,'traminlam2',traminlam2 print*,'transverse corrected ratio',traminlam2/traminlam1 print*,'Doing Transverse Bootstrap Error Analysis ...' print*,'Transverse window',float(nptswin(tracur)*dt) call boot(2,nptswin(tracur),nbag,traangwin(tracur),trajmin, & 0,traangerr,tradelerr) end if c23456789012345678901234567890123456789012345678901234567890123456789012 c print*,'imin=',imin write(*,66) write(*,66) write(*,78)eigangwin(eigcur),eigangerr,eigdelwin(eigcur),eigdelerr write(*,79)traangwin(tracur),traangerr,tradelwin(tracur),tradelerr write(*,66) 66 format(/) 78 format(4x,'EIGENVALUE: ANGLE = ',f6.2,' +/-',f6.2,3x, +'DELAY = ',f6.3,' +/-',f6.2,1x,'sec.') 79 format(4x,'TRANSVERSE: ANGLE = ',f6.2,' +/-',f6.2,3x, +'DELAY = ',f6.3,' +/-',f6.2,1x,'sec.') 99 format(2x,'Window length',2x,i4,4x,f6.2,2x,f6.3,4x,f6.2,2x,f6.3) stop end c c######################################################################## c c subroutine boot (npts,ang,lag,angerr,delerr,pdi) c c author - E Sandvol, (simplified by M Evans, July 2003) c date - unknown c c Purpose : c Perform bootstrap error analysis c c Parameters : max - Maximum size of input array (integer) c maxbag - Maximum size of array used in bootstrapping (integer) c nshift - Number of data points searched through in lag time (integer) c c input : method - Measurement method, 1=eigenvalue, 2=transverse (integer) c npts - Number of points in time series (integer) c nbag - Number of bootstrap iterations (integer) c ang - Fast polorization direction (real) c lag - Delay time [in timesteps NOT seconds] (integer) c pdi - Polarization direction [degrees] (real) c c output : angerr - Angle error, calculated from 2 standard deviations (real) c delerr - Time lag error, calculated from 2 standard deviations (real) c c internal : ampmax1 - RMS amplitude of original noise sequence (real) c ampmax2 - RMS amplitude of generated noise sequence (real) c angbag - Array of fast polarization measurements (real(maxbag)) c angdif - Difference between 'angmn' and each element of 'angbag' (real) c angmn - Mean polarization angle [degrees] (real) c az1 - Angle of first time series used in rotate commands (real) c az2 - Angle of second time series used in rotate commands (real) c cmp1_cop - Copy of 'cmp1' in common block '/waveforms/' (real(max)) c cmp2_cop - Copy of 'cmp2' in common block '/waveforms/' (real(max)) c cov - Covariance martix of particle motion (real(2,2)) c delbag - Array of lag time measurements (real(maxbag)) c delmn - Mean lag time [seconds] (real) c doeig - Flag set to 1 if eigenvalue measurement to be made (integer) c dotrans - Flag set to 1 if transverse energy measurement to be made (integer) c drot - Unit of gridsearch for angle (integer) c dum2(2*max) - Temp copy of slwdum (real(2*max)) c dumf(2*max) - Fast component of each bootstrap iteration (real(2*max)) c dums(2*max) - Slow component of each bootstrap iteration (real(2*max)) c dumnv - Stores 'slow' values between each iteration (real(-2*nshift:2*max)) c emin - Current minimum transverse energy (real) c en - Transverse energy on current iteration (real) c fast - Time series of fast polarized wavelet (real(2*max)) c i - General loop variable (integer) c ii - Counter variable in one loop (integer) c imin - Stores 'i' value of found minimum (integer) c j - Used to loop over all values of 'nshift' (integer) c jmin - Value of 'j' at current minimum (integer) c lam1 - First eigenvalue (real) c lam2 - Second eigenvalue (real) c nn - Loop over each bootstrap iteration (integer) c nrot - Number of rotations to be performed (integer) c nro2 - Used to calculate 'rota' (real) c nstart - Used to randomise offset of noise sequence (integer) c phi - Azimuth of noise sequences (real) c rand_dump - Output from inbuilt 'rand' function when seeded with current time (real) c rv - Array of Gaussian white noise (real(2*max)) c rota - Angle seismograms rotated by over each iteration of 'i' (real) c slow - Time series of slow polarized wavelet, offset by 'lag' (real(-2*nshift:2*max)) c slwdum - Version of 'slow' same size as 'fast' (real(2*max)) c vec1 - First eigenvector (real(2)) c vec2 - Second eigenvector (real(2)) c xaz - Azimuth of 'xx' (real) c x - Lagged time series used in calculations (real(2*max)) c xfast - Copy of 'fast', fast polarized wavelet (real(2*max)) c xx - Rotated time series calculated from 'dumf' and 'dums' (real(2*max)) c yaz - Azimuth of 'yy' (real) c y - Lagged time series used in calculations (real(2*max)) c yy - Rotated time series calculated from 'dumf' and 'dums' (real(2*max)) c c common : /waveforms/ c c calls subroutines : convolve c copy c eigen2x2 c energy c part_motion c rand (inbuilt) c rotate c rotbaz c shift c time (inbuilt) c zero c c################################################################### subroutine boot(method,npts,nbag,ang,lag,pdi,angerr,delerr) implicit none integer max, maxbag, nrot, nshift parameter (max=600,maxbag=200,nrot=181,nshift=121) real ang, pdi integer lag, method, nbag, npts real angerr, delerr real ampmax1, ampmax2, angbag(maxbag), angdif, angmn, az1, az2 real cmp1_cop(max), cmp2_cop(max), cov(2,2), delbag(maxbag) real delmn, dum2(2*max), dumf(2*max), dums(2*max) real dumnv(-2*nshift:2*max), emin, en, fast(2*max), imin, lam1 real lam2, nro2, phi, rand_dump, rv(2*max), rota real slow(-2*nshift:2*max), slwdum(2*max), vec1(2), vec2(2), xaz real x(2*max), xfast(2*max), xx(2*max), yaz, y(2*max), yy(2*max) real lowang,uppang,lowdel,uppdel c TEST STUFF integer ndf real eig_grid, tra_grid real grid_min,idfast,idtlag, dtlag, dfast c END TEST STUFF integer doeig, dotrans, drot, i, ii, j, jmin, nn, nstart real baz, caz1, caz2, cmp1, cmp2, dt integer plt common /grid/ eig_grid(nrot,nshift), tra_grid(nrot,nshift) common /waveforms/ baz, caz1, caz2, cmp1(max), cmp2(max), dt, & plt real rand integer time c c ... Rotate oringial components for fast and slow components c cc Assume phase is radially polarized: ccc phi = baz for trans, phi = pdi for eigen doeig = 0. dotrans = 0. if (method.eq.1.) doeig = 1. if (method.eq.2.) dotrans = 1. if (method.eq.1.) then phi = pdi else if (method.eq.2.) then phi = baz end if print*,'Phi= ',phi if (plt.eq.1.) then call part_motion(method,cmp1,cmp2,npts,ang,lag) endif do 10 i=1,npts cmp1_cop(i) = cmp1(i) cmp2_cop(i) = cmp2(i) 10 continue cc C Assume E-W and N-S data input az1=caz1 az2=caz2 c call rotate(cmp1_cop,cmp2_cop,ang,npts,az1,az2,1) c c shift array by lag time, keeping shifted noise c ii=0 do 20 i=-lag+1,0 ii=ii+1 slow(i)=cmp2_cop(ii) 20 continue do 30 i=1,npts-lag slow(i)=cmp2_cop(i+lag) fast(i)=cmp1_cop(i) 30 continue do 40 i=npts-lag+1,npts slow(i)=0. fast(i)=cmp1_cop(i) 40 continue c c Rotate to phi so now slow is the noise vector: c (slwdum is a kluge designed to avoid problems with negative array indexes) c c Need to put 'slow' into 'slwdum' as it has same array dimensions as 'fast' c do 50 i=1,npts slwdum(i)=slow(i) 50 continue call rotbaz(fast,slwdum,phi,npts,az1,az2) c do 55 i=1,2*max c print*,fast(i),slwdum(i) c 55 continue c ******************************************************************************* c ## Calculate Silver & Chan Errors here!!! call zndf(slwdum,npts,npts,ndf) print*,'NDF',ndf if (method.eq.1.) then call zerror_min(eig_grid,nrot,nshift,grid_min) call zerror95(eig_grid,nrot,nshift,ndf,grid_min,idfast,idtlag) else call zerror_min(tra_grid,nrot,nshift,grid_min) call zerror95(tra_grid,nrot,nshift,ndf,grid_min,idfast,idtlag) endif c ** Formulae in Nicks code: c ** dtlag = delta * idtlag * itlag_step / real(f) c ** dfast = 180. * idfast / real(nshift-1) c ** Simplify to (without interpolation): dtlag = dt * idtlag dfast = idfast c23456789012345678901234567890123456789012345678901234567890123456789012 if (method.eq.1.) then print*,'Eigenvalue S&C results: Angle error = ', & dfast,' Time error ',dtlag else print*,'Transverse S&C results: Angle error = ', & dfast,' Time error ',dtlag endif c ******************************************************************************* do 60 i=1,npts slow(i)=slwdum(i) 60 continue do 70 i=-lag+1,npts-lag dumnv(i)=slow(i) 70 continue c c Create nbag number of different hopefully random inversions: c c Seed random number generator using current time (in seconds) c rand_dump = rand(time()) do 80 nn=1,nbag c c Generate Random noise vector c c Create Gaussian White noise sequence: c do 90 i=-lag+1,npts-lag rv(i+lag)=2*(rand(0)-.5) slwdum(i+lag)=dumnv(i) 90 continue c ***Extend White noise vector to avoid edge effects do 100 i=npts+1,2*npts rv(i)=2*(rand(0)-.5) 100 continue c c Calculate RMS amplitude of original noise sequence: c ampmax1=0. do 110 i=1,npts ampmax1=ampmax1+(slwdum(i)*slwdum(i)) 110 continue ampmax1=sqrt(ampmax1) c***Cyclic stuff here nstart=ifix(float(npts)*(rand(0))) c print*,'nstart=',nstart ii=0. do 120 i=1,npts 120 dum2(i)=slwdum(i) do 130 i=1,npts if (i+nstart.le.npts) then slwdum(i)=-dum2(i+nstart) else slwdum(i)=-dum2(ii+1) ii=ii+1 endif 130 continue call convolve(slwdum,npts,rv,2*npts) c****Normalize noise vector ampmax2=0. do 140 i=1,npts ampmax2=ampmax2+(slwdum(i)*slwdum(i)) 140 continue ampmax2=sqrt(ampmax2) do 150 i=1,npts slwdum(i)=(ampmax1/ampmax2)*slwdum(i) 150 continue c c Uncorrect the two components for the solved fast direction and lag time c call copy(npts,fast,xfast) az1=-ang az2=-ang+90 call rotbaz(xfast,slwdum,-phi,npts,az1,az2) do 170 i=1,npts slow(i)=slwdum(i) 170 continue do 180 i=1,npts dums(i)=slow(i-lag) dumf(i)=xfast(i) 180 continue az1=0. az2=90. call rotate(dumf,dums,-ang,npts,az1,az2,1) c c Resolve Shear Wave splitting Parameters for Rerotated orthogonal pair : c c if(drot.eq.1.00) nrot = 180 c if(drot.eq.2.00) nrot = 90 c if(drot.eq.5.00) nrot = 36 c if(drot.eq.10.0) nrot = 18 c drot = 1. c nrot = 180 nro2 = ((nrot-1) / 2) + 1 emin=1.0e30 c c ... START ROTATING ORIGINAL SEISMOGRAMS c c nrot = nrot + 1 c do 190 i = 1,nrot c call zero (x,max) call zero (y,max) call zero (xx,max) call zero (yy,max) c call copy(npts,dumf,xx) call copy(npts,dums,yy) c xaz=0.0 yaz=90.0 rota = drot * float(nro2-i) c call rotate(xx,yy,rota,npts,xaz,yaz,1) c c ... START SHIFTING ROTATED SEISMOGRAMS c do 200 j = 0,nshift-1 c call copy(npts,xx,x) call copy(npts,yy,y) c c call shift(y,npts,j) c if(doeig.eq.1.0) then c call covariance(x,y,npts-j,cov) call eigen2x2(cov,lam1,lam2,vec1,vec2) en = lam2 endif if(dotrans.eq.1.0) then call rotbaz(x,y,baz,npts,xaz,yaz) call energy(npts,y,en) en=en/float(npts-j) endif c if(en.ge.emin) goto 200 emin=en imin=i jmin=j 200 continue 190 continue c c ... FINDING BOOTSTRAPPED FAST POLARIZATION DIRECTION AND TIME DELAY c angbag(nn) = drot * float(nro2-imin) delbag(nn) = float(jmin)*dt 80 continue c c Do boot strap statistics with results of randomized inversions above c delmn=0. angmn=0. c c first, switch parameter spaces from azimuth space c -90 to 90 to ang-45 to ang+45: c do 220 i=1,nbag if (angbag(i).lt.ang-90) angbag(i)=angbag(i)+180. if (angbag(i).gt.ang+90) angbag(i)=angbag(i)-180. if (angbag(i).gt.45.+ang) then angbag(i)=angbag(i)-90 delbag(i)=-delbag(i) endif if (angbag(i).lt.ang-45.) then delbag(i)=-delbag(i) angbag(i)=angbag(i)+90 endif c c Output to screen now corrected.... c if (plt.eq.1.) then print*,'Bootstrap values',method,angbag(i),delbag(i) endif 220 continue call percentile(nbag,angbag,16,lowang,84,uppang) if (method.eq.1.) then print *,'EIGENVALUE ANGLE PERCENTILES',lowang,uppang else print *,'TRANSVERSE ANGLE PERCENTILES',lowang,uppang endif call percentile(nbag,delbag,16,lowdel,84,uppdel) if (method.eq.1.) then print *,'EIGENVALUE LAG PERCENTILES',lowdel,uppdel else print *,'TRANSVERSE LAG PERCENTILES',lowdel,uppdel endif c c c Compute standard error c do 230 i=1,nbag delmn=delmn+delbag(i) angmn=angmn+angbag(i) 230 continue angmn=angmn/float(nbag) delmn=delmn/float(nbag) angerr=0. delerr=0. do 240 i=1,nbag angdif=abs(angmn-angbag(i)) if (angdif.gt.90) then angdif=90-mod(angdif,90) endif angerr=angerr+(angdif*angdif) delerr=delerr+(delmn-delbag(i))**2 240 continue angerr=2*(sqrt((angerr)/float(nbag-1))) delerr=2*(sqrt((delerr)/float(nbag-1))) return end c######################################################################## c c subroutine convolve (x,nx,y,ny) c c author - E Sandvol c date - unknown c c Purpose : c Convolve arrays x and y c c Parameters : max - Maximum size of input array (integer) c c input : nx - Size of array 'x' (integer) c ny - Size of array 'y' (integer) c x - Input noise vector (real(max)) c y - Input noise vector (real(max)) c c output : x - Output two vectors convolved together (real(max)) c c internal : i - Loop variable (integer) c j - Loop variable (integer) c y1 - Temporary array created in maths (real(-max:max)) c y2 - Temporary array created in maths (real(-max:max)) c y3 - Temporary array created in maths (real(max)) c z - Temporary value (real) c c common : none c c calls subroutines : none c c################################################################### subroutine convolve(x,nx,y,ny) implicit none integer max parameter(max=600) real x(max), y(max) integer nx, ny dimension y1(-max:max), y2(-max:max) real y1, y2, y3(max), z integer i, j j = 0 do i = 1, ny y2(i-1) = y(i) end do do i = 1, ny y2(-i+1) = y(i) end do do i = 1, nx y1(i-1)=x(i) end do c**** do convolution do 10 i = 0 , nx + ny - 2 if ( i.lt.nx-1) then z = 0. do 20 j = 0,i z = z + y2(i-j)*y1(j) 20 continue elseif(i.ge.(nx-1).and.i.le.(ny-1)) then z = 0. do 30 j = 0,nx-1 z = z + y2(i-j)*y1(j) 30 continue else z = 0. do 40 j = i-ny-1, nx-1 z = z + y2(i-j)*y1(j) 40 continue endif y3(i)=z 10 continue do 50 i=nx,2*nx x(i-nx+1)=y3(i) 50 continue return end c######################################################################## c c subroutine copy (np,x1,x2) c c author - E Sandvol c date - unknown c c Parameters : none c c Purpose : c Copies array x1 into x2 c c input : np - array dimension (integer) c x1 - array to be copied (real(np)) c c output : x2 - copy of x1 (real(np)) c c internal : i - loop variable (integer) c c common : none c c calls subroutines : none c c################################################################### subroutine copy(np,x1,x2) implicit none real x1(np) integer np real x2(np) integer i do 10 i =1,np x2(i)=x1(i) 10 continue return end c######################################################################## c c subroutine covariance (x,y,np,cov) c c author - N Teanby c date - July 2002 c c Purpose : c Calculates covariance matrix of two c vectors, (x and y) c c Parameters : none c c input : x - Time series (real(np)) c y - Time series (real(np)) c np - Number ot points (integer) c c output : cov - Covariance matrix (real(2,2)) c c internal : i - Loop variable (integer) c c common : none c c calls subroutines : none c c################################################################### subroutine covariance(x,y,np,cov) implicit none real x(np), y(np) integer np real cov(2,2) integer i c ** calc normalised cross correlation coefficient ** cov(1,1) = 0. cov(2,2) = 0. cov(1,2) = 0. cov(2,1) = 0. do 10 i=1,np cov(1,1) = cov(1,1) + x(i)**2 cov(2,2) = cov(2,2) + y(i)**2 cov(1,2) = cov(1,2) + x(i)*y(i) 10 continue cov(2,1) = cov(1,2) return end c######################################################################## c c subroutine eigen2x2 (matrix,lambda1,lambda2,vec1,vec2) c c author - N Teanby c date - July 2002 c c Purpose : c Calculates eigenvalues & eigenvectors of a 2x2 matrix c c Parameters : none c c input : matrix - 2x2 matrix (real(2,2)) c c output : lambda1 - Larger eigenvalue (real) c lambda2 - 2nd eigenvalue (real) c vec1 - Eigenvector corresponding to larger eigenvalue (real(2)) c vec2 - 2nd eigenvector (real(2)) c c internal : a \ c b > Variables used in quadratic equation (real) c c / c temp - Used to flip lambda1 & lambda2 [if required] (real) c norm - Used to normalise eigenvalues (real) c c common : none c c calls subroutines : none c c################################################################### subroutine eigen2x2(matrix,lambda1,lambda2,vec1,vec2) implicit none real matrix(2,2) real lambda1, lambda2, vec1(2), vec2(2) real a, b, c, temp, norm c ** EIGENVALUES ** c ** eigenvalue are the solution of a quadratic eqn with c.f.s ** a = 1. b = - matrix(1,1) - matrix(2,2) c = matrix(1,1)*matrix(2,2) - matrix(1,2)*matrix(2,1) lambda1 = 0.5*( -b + sqrt(b**2 - 4*a*c)) lambda2 = 0.5*( -b - sqrt(b**2 - 4*a*c)) c ** order the eigenvalues so that lambda2 is the smallest ** if (lambda2.gt.lambda1) then temp = lambda1 lambda1 = lambda2 lambda2 = temp endif C ** EIGENVECTORS (normalised) ** vec1(1)=1. vec1(2)=(lambda1-matrix(1,1))/matrix(1,2) norm = sqrt(vec1(1)**2 + vec1(2)**2) vec1(1)=vec1(1)/norm vec1(2)=vec1(2)/norm vec2(1)=1. vec2(2)=(lambda2-matrix(1,1))/matrix(1,2) norm = sqrt(vec2(1)**2 + vec2(2)**2) vec2(1)=vec2(1)/norm vec2(2)=vec2(2)/norm return end c######################################################################## c c subroutine energy (np,cc,sum) c c author - E Sandvol c date - July 2003 c c Purpose : c Calculates energy on time series c c Parameters - none c c input : np - Array dimension (integer) c cc - Time series (real(np)) c c output : sum - Energy (real) c c internal : i - Loop variable (integer) c tsq - Energy at each component (real) c c common : none c c calls subroutines : none c c################################################################### subroutine energy(np,cc,sum) implicit none real cc(np) integer np real sum real tsq integer i sum = 0.0 do 10 i=1,np tsq = cc(i)*cc(i) sum = sum + tsq 10 continue return end c######################################################################## c c subroutine grid_sch(method,npts) c c author - M Evans, (based on E.Sandvol's code) c date - July 2003 c c Purpose : c Grid searches for best splitting parameters c c Parameters : max - Maximum array size (integer) c nshift - Number of data points searched through in lag time (integer) c c input : method - Measurement method, 1=eigenvalue, 2=transverse, 3=both (integer) c npts - Number of points in arrays (integer) c c output : none - all output to common blocks c c internal : cov - Covariance martix of particle motion (real(2,2)) c doeig - Flag set to 1 if eigenvalue measurement to be made (integer) c dotrans - Flag set to 1 if transverse energy measurement to be made (integer) c drot - Unit of gridsearch for angle (integer) c eigimin - Stores 'i' value of eigenvalue minimum (integer) c emin - Current minimum energy (real) c en - Energy on current corrected trace (real) c i - Used to loop over all values of 'nrot' (integer) c j - Used to loop over all values of 'nshift' (integer) c lam1 - First eigenvalue (real) c lam2 - Second eigenvalue (real) c nrot - Number of rotations to be performed (integer) c nro2 - Used to calculate 'rota' (real) c rota - Angle cmp1 and cmp2 are rotated through for each i value c traimin - Stores 'i' value of transverse minimum (integer) c vec1 - First eigenvector (real(2)) c vec2 - Second eigenvector (real(2)) c x - Lagged time series used in calculations (real(max)) c xaz - Azimuth of time series in array 'xx' (real) c xx - Rotated time series calculated from 'cmp1' and 'cmp2' (real(max)) c y - Lagged time series used in calculations (real(max)) c yaz - Azimuth of time series in array 'yy' (real) c yy - Rotated time series calculated from 'cmp1' and 'cmp2' (real(max)) c c common : /eigen/ c /trans/ c /waveforms/ c c calls subroutines : copy c covariance c eigen2x2 c energy c pm_polar c rotate c rotbaz c shift c zero c c################################################################### subroutine grid_sch(method,npts) implicit none integer max, nrot, nshift parameter(max=600, nrot=181, nshift=121) integer method, npts real cov(2,2), emin, en, lam1, lam2, nro2, rota, vec1(2), vec2(2) real x(max), xaz, xx(max),y(max), yaz, yy(max) integer doeig, dotrans, drot, eigimin, i, j, traimin real eigang, eigdelay, eigminlam1, eigminlam2, eigpol, minvec1 real minvec2 real eig_grid, tra_grid integer eigjmin real traang, tradelay, traminlam1, traminlam2 integer trajmin real baz, caz1, caz2, cmp1, cmp2, dt integer plt common /eigen/ eigang, eigdelay, eigjmin, eigminlam1, eigminlam2, & eigpol, minvec1(2), minvec2(2) common /grid/ eig_grid(nrot,nshift), tra_grid(nrot,nshift) common /trans/ traang, tradelay, trajmin, traminlam1, traminlam2 common /waveforms/ baz, caz1, caz2, cmp1(max), cmp2(max), dt, & plt doeig = 0. dotrans = 0. if (method.eq.1.) then print *,'EIGMINLAM', eigminlam2 doeig = 1. eigminlam2=1.0e30 print *,'EIGMINLAM', eigminlam2 end if if (method.eq.2.) then dotrans = 1. emin=1.0e30 end if if (method.eq.3.) then doeig = 1. eigminlam2=1.0e30 dotrans = 1. emin=1.0e30 end if c c if(drot.eq.1.00) nrot = 180 c if(drot.eq.2.00) nrot = 90 c if(drot.eq.5.00) nrot = 36 c if(drot.eq.10.0) nrot = 18 c drot=1. c nrot = 180 nro2 = ((nrot-1) / 2) + 1 c c ... START ROTATING ORIGINAL SEISMOGRAMS c c nrot = nrot + 1 c do 10 i = 1,nrot call zero (x,max) call zero (y,max) call zero (xx,max) call zero (yy,max) call copy(max,cmp1,xx) call copy(max,cmp2,yy) c xaz=caz1 yaz=caz2 rota = drot * float(nro2-i) call rotate(xx,yy,rota,max,xaz,yaz,1) c c ... START SHIFTING ROTATED SEISMOGRAMS c do 20 j = 0,nshift-1 call copy(max,xx,x) call copy(max,yy,y) c c call shift(y,max,j) if(doeig.eq.1.0) then c call covariance(x,y,npts,cov) call eigen2x2(cov,lam1,lam2,vec1,vec2) eig_grid(i,j+1) = lam2 c if (i.eq.45.) then c if (j.eq.9.) then c write(0,*)'X= ',x c write(0,*)'Y= ',y c write(0,*)'COV= ',COV c endif c endif if(lam2.lt.eigminlam2) then eigimin=i eigjmin=j c Added v2.0 eigminlam1=lam1 eigminlam2=lam2 minvec1(1)=vec1(1) minvec1(2)=vec1(2) minvec2(1)=vec2(1) minvec2(2)=vec2(2) end if endif if(dotrans.eq.1.0) then call rotbaz(x,y,baz,npts,xaz,yaz) call energy(npts,y,en) en=en/float(npts-j) tra_grid(i,j+1) = en if(en.lt.emin) then emin = en call covariance(x,y,npts,cov) call eigen2x2(cov,traminlam1,traminlam2,vec1,vec2) traimin=i trajmin=j end if endif c c 20 continue 10 continue c c ... FINDING FAST POLARIZATION DIRECTION AND TIME DELAY c if(doeig.eq.1.0) then c eigang = drot * float(nro2-eigimin) eigdelay = float(eigjmin)*dt call pm_polar(eigang,minvec1,eigpol) endif if(dotrans.eq.1.0) then traang = drot * float(nro2-traimin) tradelay = float(trajmin)*dt endif print *,'EIGMINLAM', eigminlam2 c c c print*,'EIGENVALUE METHOD' c print*,'fast direction=',eigang c print*,'delay=',eigdelay c c print*,'eigminlam1',eigminlam1 c print*,'eigminlam2',eigminlam2 c c print*,'minvec1',minvec1 c print*,'minvec2',minvec2 c c print*,'TRANSVERSE METHOD' c print*,'fast direction=',traang c print*,'delay=',tradelay c c write(*,99)npts,eigang,eigdelay,traang,tradelay c c 99 format(2x,i4,4x,f6.2,2x,f6.3,4x,f6.2,2x,f6.3) return end c######################################################################## c c subroutine part_motion(method,cmp1,cmp2,npts,ang,lag) c c author - M Evans c date - August 2003 c c Purpose : c Outputs particle motion for GMT (if required) c c Parameters : max - Maximum size of input array (integer) c c input : ang - Fast polorization direction (real) c cmp1 - North component (real(max)) c cmp2 - East component (real(max)) c lag - Delay time [in timesteps NOT seconds] (integer) c method - Measurement method, 1=eigenvalue, 2=transverse (integer) c npts - Number of points in seismogram (integer) c c output : none c c internal : az1 - Azimuth of 'xx' (real) c az2 - Azimuth of 'yy' (real) c cmax - Maximum displacement from all seismograms (real) c i - Loop variable (integer) c xx - Rotated time series (real(max)) c yy - Rotated time series (real(max)) c c common : none c c calls subroutines : copy c rotate c shift c c################################################################### subroutine part_motion(method,cmp1,cmp2,npts,ang,lag) implicit none integer max parameter(max=600) c INPUT real ang, cmp1(max), cmp2(max) integer lag, method, npts c INTERNAL real az1, az2, cmax, xx(max), yy(max) integer i az1 = 0. az2 = 90. cmax = 0. call copy(npts,cmp1,xx) call copy(npts,cmp2,yy) c write(0,*)'DATA IN ',method,npts,ang,lag c do 4 i=1,npts c write(0,*)'PRE-JUNK ',xx(i),yy(i) c 4 continue call rotate(xx,yy,ang,npts,az1,az2,1) call shift(yy,npts,lag) c do 5 i=1,npts c write(0,*)'JUNK ',xx(i),yy(i) c 5 continue call rotate(xx,yy,-ang,npts,az1,az2) do 10 i=1,npts if (abs(cmp1(i)).gt.cmax) cmax = abs(cmp1(i)) if (abs(cmp2(i)).gt.cmax) cmax = abs(cmp2(i)) if (abs(xx(i)).gt.cmax) cmax = abs(xx(i)) if (abs(yy(i)).gt.cmax) cmax = abs(yy(i)) 10 continue do 20 i=1,npts print*,'PM PLOTS',method,cmp1(i)/cmax,cmp2(i)/cmax, & xx(i)/cmax,yy(i)/cmax 20 continue return end c######################################################################## c c subroutine pm_polar (fast, vec1, pol) c c author - M Evans c date - July 2003 c c Purpose : c Calculates polarization direction of particle motion c from larger eigenvector and fast direction c c Parameters: pi - Mathematical constant c c input : fast - Fast polarization direction (real) c vec1 - Eigenvector (real(2)) c c output : pol - Polarization of particle motion (real) c c internal : none c c common : none c c calls subroutines : none c c################################################################### subroutine pm_polar(fast,vec1,pol) implicit none real pi parameter (pi=3.141592654) real fast,vec1(2) real pol c ** polarization in azimuth degrees(clockwise from north) c note, eigenvectors are in rotated coord frame and must add c fast to the angle. ** pol = fast + atan2(vec1(2),vec1(1))*180./pi return end c######################################################################## c c subroutine rotate (x,y,ang,np,azx,azy,kp) c c author - E Sandvol c date - unknown c c Purpose : c Rotates two orthogonal time series c through requested angle & polarity c c Parameters : pi - Mathematical constant c c input : x - Time series (real(np)) c y - Time series (real(np)) c ang - Angle to be rotated through (real) c np - Array dimension of x & y (integer) c azx - Azimuth of x (real) c azy - Azimuth of y (real) c kp - Polarity flag [set to -1 for flip] (integer) c c output : x - rotated component (real(np)) c y - transverse component [positive clockwise from radial] (real(np)) c c internal : angr - Angle 'ang' converted to radians (real) c c - Cosine of 'angr' (real) c i - Loop variable (integer) c s - Sine of 'angr' (real) c xc - New value of 'x(i)' (real) c yc - New value of 'y(i)' (real) c c common : none c c calls subroutines : none c c################################################################### subroutine rotate(x,y,ang,np,azx,azy,kp) implicit none real pi parameter (pi=3.141592654) real ang, azx, azy, x(np), y(np) integer kp, np real angr, c, s, xc, yc integer i azx = azx + ang azy = azy + ang angr = ang * pi / 180. s = sin(angr) c = cos(angr) c do 10 i = 1,np xc = x(i)*c + y(i)*s yc =-x(i)*s + y(i)*c x(i) = xc y(i) = yc 10 continue c if(kp.eq.-1) then do 11 i = 1,np 11 y(i)=-y(i) azy=180-azy else endif return end c######################################################################## c c subroutine rotbaz (x1,x2,baz,np,az1,az2) c c author - E Sandvol c date - unknown c c Purpose : c Rotates two horizontal time series into c radial and transverse components c c Parameters : pi - Mathematical constant c c input : x1 - Time series (real(np)) c x2 - Time series (real(np)) c baz - Back azimuth (real) c np - Array dimension of x1 & x2 (int) c az1 - Azimuth of x1 (real) c az2 - Azimuth of x2 (real) c c output : x1 - Radial component [positive away from source] (real(np)) c x2 - Transverse component [positive clockwise from radial] (real(np)) c c internal : aa - Used to calculate sine of 'radbaz' - 'radaz1' (real) c azck - Angle perpendicular to az1 (real) c bb - Used to calculate cosine of 'radbaz' - 'radaz1' (real) c diff - Angle between 'az2' and 'azck' (real) c i - Loop variable (integer) c radaz1 - Angle 'az1' converted to radians (real) c radbaz - Angle 'baz' converted to radians (real) c radial - New value of 'x1(i)' (real) c transv - New value of 'x2(i)' (real) c c common : none c c calls subroutines : none c c################################################################### subroutine rotbaz(x1,x2,baz,np,az1,az2) implicit none real pi parameter (pi=3.141592654) real az1, az2, baz, x1(np), x2(np) integer np integer i real aa, azck, bb, diff, radaz1, radbaz, radial, transv c azck = az1 + 90.0 diff = abs(azck-az2) if(diff.lt.0.01) goto 10 if(diff.lt.179.)then write(6,*)' PROBLEM IN ROTATION ' return else endif 10 radbaz = baz * pi / 180. radaz1 = az1 * pi / 180. aa=sin(radbaz-radaz1) bb=cos(radbaz-radaz1) do 20 i = 1,np radial = - x1(i)*bb - x2(i)*aa transv = x1(i)*aa - x2(i)*bb x1(i) = radial x2(i) = transv 20 continue return end c c######################################################################## c c subroutine sacin(f1,f2,npts) c c author - E Sandvol, (simplified by M Evans, July 2003) c date - unknown c c Purpose : c Read in both SAC data files c c Parameters : max - Maximum trace size (integer) c c input : f1 - Filename of north component (character) c f2 - Filename of east component (character) c npts - Number of points (integer) c c output : none (all output to common /waveforms/ block) c c internal : az1 - Azimuth of 'xx' (real(max)) c az2 - Azimuth of 'yy' (real(max)) c beg - Begining value of dependant variable (see SAC manual) c cmax - Maximum amplitude on radial traces (real) c i - Loop variable (integer) c nerr - SAC error flag (see SAC manual) c xx - Copy of 'cmp1' in common block '/waveforms/' (real(max)) c yy - Copy of 'cmp2' in common block '/waveforms/' (real(max)) c c common : /waveforms/ c c calls subroutines : copy c getfhv (SAC library) c getnhv (SAC library) c rotbaz c rsac1 (SAC library) c c################################################################### subroutine sacin(f1,f2,npts) c implicit none integer max parameter(max=600) integer npts character*60 f1, f2 real az1, az2, beg, cmax, nerr, xx(max), yy(max) integer i real baz, caz1, caz2, cmp1, cmp2, dt integer plt common /waveforms/ baz, caz1, caz2, cmp1(max), cmp2(max), dt, & plt c print*,' ABOUT TO READ SAC FILES' call rsac1(f1,cmp1,npts,beg,dt,max,nerr) call getfhv('cmpaz',caz1,nerr) call rsac1(f2,cmp2,npts,beg,dt,max,nerr) call getfhv('cmpaz',caz2,nerr) print*,' SAC STUFF DONE OK' c if(nerr.eq.0) goto 10 c write(6,*)' ATTENTION MAX DIMENSION IS EXCEDEED,/, * TERMINATION .... ' return c 10 call getfhv('baz',baz,nerr) call getnhv('npts',npts,nerr) c c IF required, plot seismograms c if (plt.eq.1.) then call copy(npts,cmp1,xx) call copy(npts,cmp2,yy) az1 = caz1 az2 = caz2 call rotbaz(xx,yy,baz,npts,az1,az2) cmax = 0. do 20 i=1,npts if (abs(xx(i)).gt.cmax) cmax = abs(xx(i)) if (abs(yy(i)).gt.cmax) cmax = abs(yy(i)) 20 continue do 30 i=1,npts print*,'Seismograms',xx(i)/cmax,yy(i)/cmax 30 continue endif return end c######################################################################## c c subroutine shift(x,np,ish) c c author - E Sandvol, (simplified by M Evans, July 2003) c date - unknown c c Purpose : c moves data in array x by ish positions c open holes are filled with zeroes c c Parameters : none c c input : x - Time series (real(np)) c np - Size of array (integer) c ish - Number of positions moved (integer) c c output : x - shifted time series (real(np)) c c internal : i - Loop variable (integer) c istp - Number of iterations of first loop (integer) c istrt - Number of iterations of second loop (integer) c c common : none c c calls subroutines : none c c################################################################### subroutine shift(x,np,ish) implicit none real x(np) integer ish, np integer i, istp, istrt c c ... data will be shifted up the array : advance c istp = np - ish do 30 i=1,istp x(i) = x(i+ish) 30 continue istrt = istp + 1 do 40 i=istrt,np x(i) = 0.0 40 continue return end c######################################################################## c c subroutine subroutine window_select(angwin,delwin,np,win_end) c c author - M Evans, (based on E Sandvols' code) c date - July 2003 c c Purpose : c Selects appropriate window using largest plateau method c c Parameters : none c c input : angwin - Array of polarization angle measurements, (real(np)) c delwin - Array of lag time measurements, (real(np)) c np - Number of points in arrays (integer) c c output : win_end - Array element of selected window (integer) c c internal : cur_angle - Polarization angle of current plateau (real) c cur_start - Array element of first value in current plateau (integer) c cur_time - Lag time of current plateau (real) c cur_value - Number of iterations in current plateau (integer) c i - Used to check both angle and lag times match (integer) c j - Loop variable (integer) c max_angle - Polarization angle of best plateau (real) c max_start - Array element of first value in best plateau (integer) c max_time - Lag time of best plateau (real) c max_value - Number of iterations in best plateau (integer) c c common : none c c calls subroutines : none c c################################################################### subroutine window_select(angwin,delwin,np,win_end) implicit none real angwin(np), delwin(np) integer np integer win_end real cur_angle, cur_time, max_angle, max_time integer cur_start, cur_value, i, j, max_start, max_value cur_angle = -100. cur_time = -1. cur_value = 1. max_angle = -100. max_time = -1. max_value = 1. do 10 j=1,np i = 0. if (angwin(j).eq.cur_angle) i=i+1. if (delwin(j).eq.cur_time) i=i+1. if (i.eq.2.) then cur_value = cur_value+1. if (cur_value.gt.max_value) then max_angle = cur_angle max_time = cur_time max_value = cur_value max_start = cur_start end if else cur_angle = angwin(j) cur_time = delwin(j) cur_value = 1. cur_start = j end if 10 continue c ** Following line only works due to F77 integer division** c win_end = max_start + (max_value / 2.) return end c######################################################################## c c subroutine zero(x,np) c c author - E Sandvol c date - unknown c c Purpose : c Zeroes out array c c Parameters : none c c input : np - Array dimension (integer) c x - Array to be zeroed out (real(np)) c c output : x - Blank array (real(np)) c c internal : i - Loop variable (integer) c c common : none c c calls subroutines : none c c################################################################### subroutine zero(x,np) implicit none real x(np) integer i, np do 10 i=1,np x(i) = 0.0 10 continue return end c################################################################### subroutine percentile(np,data,lower,lowper,upper,uppper) implicit none real data(np) integer np, lower, upper real lowper, uppper real x integer i, j, e do 10 i = 1, np do 20 j=i, np if (data(i).gt.data(j)) then x = data(i) data(i) = data(j) data(j) = x endif 20 continue 10 continue c** NOTE: Rounds off to lowest integer due to integer division e = (np * lower) / 100 lowper = data(e) print*,'RESULT',lower,e,lowper e = (np * upper) / 100 uppper = data(e) print*,'RESULT',upper,e,uppper return end c----------------------------------------------------------------------- subroutine zndf(y,n,norig,ndf) c----------------------------------------------------------------------- c c subroutine to calculate number of degrees of freedom (ndf) of a time c series y, with length n, and physical dimension max. If y has been c interpolated then norig is the original number of points in the shear c wave window c c calculated according to appendix in Silver and Chan 1991 c c For a gaussian white noise ndf should be equal to n. c The true ndf should be less than n c c variables c in: c y(max) real time series c n int number of points c norig int original (uninterpolated) number of points c in the shear wave window c out: c ndf int number of degrees of freedom c local: c max int array dimension c (read from SIZE_max.h at compile time) c n2 int next power of 2 up from n c yint(max) real y interpolated onto n2 points c yfft(max*2) real fft of interpolated (to n2) y in NR format c yfft_amp(max) real amplitude spectum c nf int number of freq points for n data points c yfft_amp real interpolated amp spec on nf points c f2,f4 real f2 and f4 from Silver and Chan 1991 c used to calc ndf c c----------------------------------------------------------------------- c N. Teanby 5-7-02 Original code, based on ndf_fun2.f c and ndf_spect.f in DTM c----------------------------------------------------------------------- implicit none integer ndf integer n,max,n2,nf,norig parameter (max=600) real y(max),yint(max),yfft(max*2),yfft_amp(max) real f2,f4 integer i c -- calc the fft of y -- c ** calc n2, the next power of 2 from n (or n if n=power of 2) ** n2 = 2**(int(log(real(n)-0.1)/log(2.0)) + 1) c ** interpolate to have n2 points ** call zlinint(y,n,max,n2,yint) c ** perform fft (1 means do forward fft)** call zfft(yint,n2,max,1,yfft) c ** determine number of frequency points for original time series. nf=(norig+2)/2 c ** calculate amplitude spectrum ** do 2 i=1,nf yfft_amp(i) = sqrt(yfft(2*i-1)**2 + yfft(2*i)**2) 2 continue c -- calc f2 and f4 from Silver and Chan 1991 -- f2=0. f4=0. do 7 i=1,nf f2 = f2 + yfft_amp(i)**2 f4 = f4 + yfft_amp(i)**4 c -- for zero frequency and for Nyquist: if(i.eq.1.or.i.eq.nf)then f2 = f2 - 0.5*yfft_amp(i)**2 f4 = f4 - 0.5*yfft_amp(i)**4 endif 7 continue ndf = nint( 2.0 * (2.0*f2**2/f4 - 1.) ) print*,ndf,f2,f4 return end c----------------------------------------------------------------------- subroutine zlinint(y,n,np,ninterp,yint) c----------------------------------------------------------------------- c c linear interpolate the series y from n points to ninterp points c c variables c input: c y(np) real series to interpolate c n int number of points c np int array dimension c ninterp int number of points to interpolate to c c yint(np) real interpolated series c c----------------------------------------------------------------------- c N. Teanby 31-7-02 Original code c----------------------------------------------------------------------- implicit none integer i,n,np,ninterp,index1,index2 real y(np),yint(np) real x,dx c ** assume spacing of orignal series (y) = 1 ** c ** spacing of interpolated series is then dx, where ** dx = (n-1.)/(ninterp-1.) c ** interpolation ** do 1 i=1,ninterp c ** calc index of the two points used to do the interpolation ** x = ( i - 1.) * dx index1 = int( x ) index2 = index1 + 1 c ** interpolate ** if (index1.eq.index2) then yint(i)=y(index1) else yint(i)=y(index1) + (y(index2)-y(index1))*(x-aint(x)) endif 1 continue return end c----------------------------------------------------------------------- SUBROUTINE zfft(y,nn,np,isign,yfft) c----------------------------------------------------------------------- c c routine to calc fft of of a real time series c c INPUT VARIABLES c y(np) real time seris c nn int number of data points c np int size of array c isign int set = 1 for fft c set = -1 for inverse fft c OUTPUT VARIABLES c yfft(2*np) real fft of y c -odd and even indices are real and imag parts, respectively c -contains nn pairs of points c c NOTES: nn must be a power of 2 (not checked for in routine) c c----------------------------------------------------------------------- c n. teanby 5-7-02 modified from Num Rec four1.f c----------------------------------------------------------------------- INTEGER isign,nn,np REAL y(np),yfft(2*np) INTEGER i,istep,j,m,mmax,n REAL tempi,tempr DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp c -- this bit added -- c ** copy y to yfft, because NR routine replaces y with it's fft ** c ** yfft: odd and even indices are real and imag parts, respectively ** do 100 i=1,nn c ** real part ** yfft(2*i-1) = y(i) c ** imaginary part (=0 because real data) ** yfft(2*i) = 0.0 100 continue c -- this bit unchanged from NR routine -- n=2*nn j=1 do 11 i=1,n,2 if(j.gt.i)then tempr=yfft(j) tempi=yfft(j+1) yfft(j)=yfft(i) yfft(j+1)=yfft(i+1) yfft(i)=tempr yfft(i+1)=tempi endif m=n/2 1 if ((m.ge.2).and.(j.gt.m)) then j=j-m m=m/2 goto 1 endif j=j+m 11 continue mmax=2 2 if (n.gt.mmax) then istep=2*mmax theta=6.28318530717959d0/(isign*mmax) wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0 wi=0.d0 do 13 m=1,mmax,2 do 12 i=m,n,istep j=i+mmax tempr=sngl(wr)*yfft(j)-sngl(wi)*yfft(j+1) tempi=sngl(wr)*yfft(j+1)+sngl(wi)*yfft(j) yfft(j)=yfft(i)-tempr yfft(j+1)=yfft(i+1)-tempi yfft(i)=yfft(i)+tempr yfft(i+1)=yfft(i+1)+tempi 12 continue wtemp=wr wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi 13 continue mmax=istep goto 2 endif return END C (C) Copr. 1986-92 Numerical Recipes Software *%&&,1{. c----------------------------------------------------------------------- subroutine zerror95(error,np1,np2,ndf,lambda2_min,ierror,jerror) c----------------------------------------------------------------------- c c subroutine to calculate normalise contours so that 1=95% confidence c contour. Also calc the errorbars on tlag and fast direction c c calculated confidence interval according to appendix in c Silver and Chan 1991 using the formula: c c lambda2 <= 1 + k f_(k,ndf-k) (1-alpha) c ------- ------- c lambda2_min ndf - k c c where: c k = number of parameters = 2 c alpha = confindence level = 0.05 for 95% confidence c f is the inverse of the F distribution c ndf = number of degrees of freedom c lambda2 = second eigen value of covaraiance matrix (= elements of error) c lambda2_min = minimum lambda2 c c calc errorbars on tlag and fast by finding the bounding rectangle of the c 95% confidence contour. c c variables c in: c error(np1,np2) real array of lambda2 (altered by routine) c np1/2 int array dimension c ndf int number of degrees of freedom c lambda2_min real minimum lambda2 (corresponds to solution) c out: c error(np1,np2) real returned normalised array of lambda2 c ierror real extent (half width) of contour in i dirn c jerror real extent (half width) of contour in j dirn c other: c lambda2_max real lambda2 corresponding to 95% confidence c c----------------------------------------------------------------------- c N. Teanby 1-8-02 Original code c----------------------------------------------------------------------- implicit none integer ndf,np1,np2,npc integer i,j,jmin,jmax,irange,jrange parameter (npc=500) integer line(npc) real ierror,jerror integer k,k1,irange_min,irange_max,istart,line_test(npc) real error(np1,np2),lambda2_min,lambda_max c external fftable real ffvalue c ** check that npc is beg enough ** if (npc.lt.np1) then pause 'ERROR: zerror95: npc < np1' endif c ** calc value of lambda at 95% confidence limit from tabulated values ** if (ndf.ge.3) then call fftable(ndf-2,ffvalue) lambda_max = lambda2_min*( 1. + 2.*ffvalue/real(ndf-2)) else print*,'WARNING: zerror95: ndf <=2, set to 3)' call fftable(1,ffvalue) lambda_max = lambda2_min*( 1. + 2.*ffvalue) endif c ** normalise errors by lambda_max, so that and error of 1 = 95% confidence ** do i=1,np1 do j=1,np2 error(i,j)=error(i,j)/lambda_max enddo enddo c ** find i/j error (half width of 95% confidence contour) ** c ** find min and max j, simply search the array ** jmin=np2 jmax=1 do i=1,np1 do j=1,np2 if (error(i,j).le.1.0) then jmin = min0(jmin,j) jmax = max0(jmax,j) endif enddo enddo jrange=jmax-jmin c ** finding min max i is more difficult because of cyclicity of angles ** c ** sweep a line over all j, set point on line equal to 1 if it falls within the 95% convidence c contour for any j. The height of the bounding rectangle is defined by the shortest line c which includes all points with a value of line(i)=1. This line is found by searching all c line lengths from the minimum = sum_i linr(i) to maximum = np1** do i=1,np1 line(i)=0 enddo do j=1,np2 do i=1,np1 if (error(i,j).le.1.0) then line(i)=1 endif enddo enddo c ** min line length ** irange_min = 0 do i=1,np1 irange_min = irange_min + line(i) enddo c ** max line length ** irange_max = np1 c ** search all line length and starting points to find irange ** do i = irange_min , irange_max do istart = 1,np1 do k = 1,np1 line_test(k)=0 enddo do k = istart,istart+i if (k.gt.np1) then k1 = k - np1 else k1 = k endif line_test(k1) = 1 enddo do k = 1,np1 if ((line(k).eq.1).and.(line_test(k).ne.1)) then goto 1 endif enddo irange = i goto 11 1 continue enddo enddo 11 continue c ** one standard deviation = 0.5*95% error half width c (so x full width of 95% contour by 0.25 to get 1s.d.)** ierror = 0.25*real(irange) jerror = 0.25*real(jrange) return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c----------------------------------------------------------------------- subroutine fftable(ndf,ffvalue) c----------------------------------------------------------------------- c c function to return value of f statistic c c specialised for k=2 and alpha=0.05 (2 params and 95% confidence) c c table downloaded from: c http://www.itl.nist.gov/div898/handbook/eda/section3/eda3673.htm c c ndf int number of degrees of freedom c c----------------------------------------------------------------------- c N. Teanby 1-8-02 Original code c----------------------------------------------------------------------- implicit none real ffvalue integer ndf real ftable_data(100) data ftable_data / 199.500, 19.000, 9.552, 6.944, 5.786, >5.143, 4.737, 4.459, 4.256, 4.103, >3.982, 3.885, 3.806, 3.739, 3.682, >3.634, 3.592, 3.555, 3.522, 3.493, >3.467, 3.443, 3.422, 3.403, 3.385, >3.369, 3.354, 3.340, 3.328, 3.316, >3.305, 3.295, 3.285, 3.276, 3.267, >3.259, 3.252, 3.245, 3.238, 3.232, >3.226, 3.220, 3.214, 3.209, 3.204, >3.200, 3.195, 3.191, 3.187, 3.183, >3.179, 3.175, 3.172, 3.168, 3.165, >3.162, 3.159, 3.156, 3.153, 3.150, >3.148, 3.145, 3.143, 3.140, 3.138, >3.136, 3.134, 3.132, 3.130, 3.128, >3.126, 3.124, 3.122, 3.120, 3.119, >3.117, 3.115, 3.114, 3.112, 3.111, >3.109, 3.108, 3.107, 3.105, 3.104, >3.103, 3.101, 3.100, 3.099, 3.098, >3.097, 3.095, 3.094, 3.093, 3.092, >3.091, 3.090, 3.089, 3.088, 3.087 / if (ndf.le.0) then c ** if ndf is below range 1-100 error ** pause 'ERROR: fftable: ndf.le.0' else if (ndf.le.100) then c ** if ndf is in range 1-100 take directly from tabulated values ** ffvalue=ftable_data(ndf) else if (ndf.le.999) then c ** if ndf is in range 100-999 use linear interpolation ** c ** f=3.087 at ndf = 100, and f=3.000 at ndf=999 ** ffvalue = 3.087 - 0.087*real(ndf-100)/899. else if (ndf.gt.999) then c ** if ndf is over range return value at 999 ** ffvalue = 3.0 print*,'WARNING: ndf.gt.999, f statistic at ndf=999 returned' endif return end c----------------------------------------------------------------------- subroutine zerror_min(error,np1,np2,lambda2_min) c----------------------------------------------------------------------- c c grid search for the minimum value of lambda2 on the interpolated c error surface c c variables c in: c error(np1,np2) real error surface (i.e. lambda2) c np1/2 int array dimensions c out: c lambda2_min real minimum value of lambda2 c itlag int index of lag corresponding to lambda2_min c ifast int index of fast dirn corresponding to lambda2_min c c----------------------------------------------------------------------- c N. Teanby 4-8-02 Original code c----------------------------------------------------------------------- implicit none integer np1,np2 real error(np1,np2),lambda2_min integer i,j c ** find the minimum lambda2 position ** lambda2_min=error(1,1) do 1 i=1,np1 do 2 j = 1,np2 if (error(i,j).lt.lambda2_min) then lambda2_min = error(i,j) endif 2 continue 1 continue return end c######################################################################## c c subroutine main_calc(method,npts) c c author - M Evans, (based on E.Sandvol's code) c date - July 2003 c c Purpose : c Grid searches for best splitting parameters c c Parameters : max - Maximum array size (integer) c nshift - Number of data points searched through in lag time (integer) c c input : method - Measurement method, 1=eigenvalue, 2=transverse, 3=both (integer) c npts - Number of points in arrays (integer) c c output : none - all output to common blocks c c internal : cov - Covariance martix of particle motion (real(2,2)) c doeig - Flag set to 1 if eigenvalue measurement to be made (integer) c dotrans - Flag set to 1 if transverse energy measurement to be made (integer) c drot - Unit of gridsearch for angle (integer) c eigimin - Stores 'i' value of eigenvalue minimum (integer) c emin - Current minimum energy (real) c en - Energy on current corrected trace (real) c i - Used to loop over all values of 'nrot' (integer) c j - Used to loop over all values of 'nshift' (integer) c lam1 - First eigenvalue (real) c lam2 - Second eigenvalue (real) c nrot - Number of rotations to be performed (integer) c nro2 - Used to calculate 'rota' (real) c rota - Angle cmp1 and cmp2 are rotated through for each i value c traimin - Stores 'i' value of transverse minimum (integer) c vec1 - First eigenvector (real(2)) c vec2 - Second eigenvector (real(2)) c x - Lagged time series used in calculations (real(max)) c xaz - Azimuth of time series in array 'xx' (real) c xx - Rotated time series calculated from 'cmp1' and 'cmp2' (real(max)) c y - Lagged time series used in calculations (real(max)) c yaz - Azimuth of time series in array 'yy' (real) c yy - Rotated time series calculated from 'cmp1' and 'cmp2' (real(max)) c c common : /eigen/ c /trans/ c /waveforms/ c c calls subroutines : copy c covariance c eigen2x2 c energy c pm_polar c rotate c rotbaz c shift c zero c c################################################################### subroutine main_calc(method,npts) implicit none integer max, nrot, nshift, nwin parameter(max=600, nrot=181, nshift=121, nwin=600) integer method, npts real cov(2,2), emin, en, lam1, lam2, nro2, rota, vec1(2), vec2(2) real x(max), xaz, xx(max),y(max), yaz, yy(max) integer doeig, dotrans, drot, eigimin, i, j, traimin integer k real xsq, ysq, xy, ccsq real eigang, eigdelay, eigminlam1, eigminlam2, eigpol, minvec1 real minvec2 real eig_grid, tra_grid integer eigjmin real traang, tradelay, traminlam1, traminlam2 integer trajmin real data real baz, caz1, caz2, cmp1, cmp2, dt integer plt common /eigen/ eigang, eigdelay, eigjmin, eigminlam1, eigminlam2, & eigpol, minvec1(2), minvec2(2) common /grid/ eig_grid(nrot,nshift), tra_grid(nrot,nshift) common /trans/ traang, tradelay, trajmin, traminlam1, traminlam2 common /values/ data(nrot,nshift,nwin,4) common /waveforms/ baz, caz1, caz2, cmp1(max), cmp2(max), dt, & plt c drot=1. c nrot = 180 nro2 = ((nrot-1) / 2) + 1 c c ... START ROTATING ORIGINAL SEISMOGRAMS c c do 10 i = 1,nrot call zero (x,max) call zero (y,max) call zero (xx,max) call zero (yy,max) call copy(max,cmp1,xx) call copy(max,cmp2,yy) c xaz=caz1 yaz=caz2 rota = drot * float(nro2-i) call rotate(xx,yy,rota,max,xaz,yaz,1) c c ... START SHIFTING ROTATED SEISMOGRAMS c do 20 j = 0,nshift-1 call copy(max,xx,x) call copy(max,yy,y) c c call shift(y,npts,j) xsq = 0. ysq = 0. xy = 0. ccsq = 0. do 30 k = 1, npts - j xsq = xsq + x(k)**2 ysq = ysq + y(k)**2 xy = xy + x(k)*y(k) data(i,j+1,k,1) = xsq data(i,j+1,k,2) = ysq data(i,j+1,k,3) = xy c if (k.eq.281.) then c if (i.eq.45.) then c if (j.eq.9.) then c write(0,*)'X= ',x c write(0,*)'Y= ',y c write(0,*)'COV= ',COV c endif c endif c endif 30 continue call rotbaz(x,y,baz,npts,xaz,yaz) do 40 k = 1, npts - j ccsq = ccsq + y(k)**2 data(i,j+1,k,4) = ccsq 40 continue c c 20 continue 10 continue c c ... FINDING FAST POLARIZATION DIRECTION AND TIME DELAY c c c c print*,'EIGENVALUE METHOD' c print*,'fast direction=',eigang c print*,'delay=',eigdelay c c print*,'eigminlam1',eigminlam1 c print*,'eigminlam2',eigminlam2 c c print*,'minvec1',minvec1 c print*,'minvec2',minvec2 c c print*,'TRANSVERSE METHOD' c print*,'fast direction=',traang c print*,'delay=',tradelay c c write(*,99)npts,eigang,eigdelay,traang,tradelay c c 99 format(2x,i4,4x,f6.2,2x,f6.3,4x,f6.2,2x,f6.3) return end c######################################################################## c c subroutine main_loop(method,npts) c c author - M Evans, (based on E.Sandvol's code) c date - July 2003 c c Purpose : c Grid searches for best splitting parameters c c Parameters : max - Maximum array size (integer) c nshift - Number of data points searched through in lag time (integer) c c input : method - Measurement method, 1=eigenvalue, 2=transverse, 3=both (integer) c npts - Number of points in arrays (integer) c c output : none - all output to common blocks c c internal : cov - Covariance martix of particle motion (real(2,2)) c doeig - Flag set to 1 if eigenvalue measurement to be made (integer) c dotrans - Flag set to 1 if transverse energy measurement to be made (integer) c drot - Unit of gridsearch for angle (integer) c eigimin - Stores 'i' value of eigenvalue minimum (integer) c emin - Current minimum energy (real) c en - Energy on current corrected trace (real) c i - Used to loop over all values of 'nrot' (integer) c j - Used to loop over all values of 'nshift' (integer) c lam1 - First eigenvalue (real) c lam2 - Second eigenvalue (real) c nrot - Number of rotations to be performed (integer) c nro2 - Used to calculate 'rota' (real) c rota - Angle cmp1 and cmp2 are rotated through for each i value c traimin - Stores 'i' value of transverse minimum (integer) c vec1 - First eigenvector (real(2)) c vec2 - Second eigenvector (real(2)) c x - Lagged time series used in calculations (real(max)) c xaz - Azimuth of time series in array 'xx' (real) c xx - Rotated time series calculated from 'cmp1' and 'cmp2' (real(max)) c y - Lagged time series used in calculations (real(max)) c yaz - Azimuth of time series in array 'yy' (real) c yy - Rotated time series calculated from 'cmp1' and 'cmp2' (real(max)) c c common : /eigen/ c /trans/ c /waveforms/ c c calls subroutines : copy c covariance c eigen2x2 c energy c pm_polar c rotate c rotbaz c shift c zero c c################################################################### subroutine main_loop(method,npts) implicit none integer max, nrot, nshift, nwin parameter(max=600, nrot=181, nshift=121, nwin=600) integer method, npts real cov(2,2), emin, en, lam1, lam2, nro2, rota, vec1(2), vec2(2) real x(max), xaz, xx(max),y(max), yaz, yy(max) integer doeig, dotrans, drot, eigimin, i, j, traimin real eigang, eigdelay, eigminlam1, eigminlam2, eigpol, minvec1 real minvec2 real eig_grid, tra_grid integer eigjmin real traang, tradelay, traminlam1, traminlam2 integer trajmin real data real baz, caz1, caz2, cmp1, cmp2, dt integer plt common /eigen/ eigang, eigdelay, eigjmin, eigminlam1, eigminlam2, & eigpol, minvec1(2), minvec2(2) common /grid/ eig_grid(nrot,nshift), tra_grid(nrot,nshift) common /trans/ traang, tradelay, trajmin, traminlam1, traminlam2 common /values/ data(nrot,nshift,nwin,4) common /waveforms/ baz, caz1, caz2, cmp1(max), cmp2(max), dt, & plt doeig = 0. dotrans = 0. if (method.eq.1.) then doeig = 1. eigminlam2=1.0e30 end if if (method.eq.2.) then dotrans = 1. emin=1.0e30 end if if (method.eq.3.) then doeig = 1. eigminlam2=1.0e30 dotrans = 1. emin=1.0e30 end if c c if(drot.eq.1.00) nrot = 180 c if(drot.eq.2.00) nrot = 90 c if(drot.eq.5.00) nrot = 36 c if(drot.eq.10.0) nrot = 18 c drot=1. c nrot = 180 nro2 = ((nrot-1) / 2) + 1 c c ... START ROTATING ORIGINAL SEISMOGRAMS c c nrot = nrot + 1 c do 10 i = 1,nrot xaz=caz1 yaz=caz2 rota = drot * float(nro2-i) call rotate(xx,yy,rota,npts,xaz,yaz,1) c c ... START SHIFTING ROTATED SEISMOGRAMS c do 20 j = 0,nshift-1 if(doeig.eq.1.0) then c cov(1,1) = data(i,j+1,npts,1) cov(2,2) = data(i,j+1,npts,2) cov(1,2) = data(i,j+1,npts,3) cov(2,1) = cov(1,2) c if (npts.eq.281.) then c if (i.eq.45.) then c if (j.eq.9.) then c write(0,*)'X= ',x c write(0,*)'Y= ',y c write(0,*)'COV= ',COV c endif c endif c endif call eigen2x2(cov,lam1,lam2,vec1,vec2) eig_grid(i,j+1) = lam2 if(lam2.lt.eigminlam2) then eigimin=i eigjmin=j c Added v2.0 eigminlam1=lam1 eigminlam2=lam2 minvec1(1)=vec1(1) minvec1(2)=vec1(2) minvec2(1)=vec2(1) minvec2(2)=vec2(2) end if endif if(dotrans.eq.1.0) then c call rotbaz(x,y,baz,npts,xaz,yaz) c call energy(npts,y,en) en = data(i,j+1,npts,4) en=en/float(npts-j) tra_grid(i,j+1) = en if(en.lt.emin) then emin = en c call covariance(x,y,npts,cov) c call eigen2x2(cov,traminlam1,traminlam2,vec1,vec2) traimin=i trajmin=j c end if endif c c 20 continue 10 continue c c ... FINDING FAST POLARIZATION DIRECTION AND TIME DELAY c if(doeig.eq.1.0) then c eigang = drot * float(nro2-eigimin) eigdelay = float(eigjmin)*dt c call pm_polar(eigang,minvec1,eigpol) endif if(dotrans.eq.1.0) then traang = drot * float(nro2-traimin) tradelay = float(trajmin)*dt endif c print*,eigang,eigdelay,eigminlam2,traang,tradelay,emin return end