Just an informational item regarding Snow Leopard's treatment of Fortran. The first comment is that programmers need to install XCode from the Snow Leapard install DVD. I bulled ahead without checking this and ran into a roadblock in trying to compile a Fortran program -- the computer could not find the unix command "as" which assembles object files. Once I loaded the XCodes package from the Install disk, this problem went away, and others cropped up.
The innovation of Snow Leopard (and Leopard, I think, Im not exactly sure when the transition occurred relative to Tiger (MacOS 10.4)) is the full utilization of 64-bit chip architecture. The Snow Leopard gcc compiler creates object files that are 64bit, but the Fortran compiler g77 does not. If you try to assemble Fortran and C routines into one archive or a single program, the computer balks. You can tell the nature of the object file with the command "file *.o". I sampled the internets for hints on solving this problem. The simplest solution, and the one Im adopting for now, is to force the gcc compiler to generate a 32-bit object file, rather than a 64bit version, with the "-m32" runtime option. Less powerful than 64bit, and Im still getting warnings from the compiler when I create large codes, but the executables seem to work OK.
The "gfortran" option, rather than "g77" for the Fortran compiler, might work with the 64bit option, but the internets are not unanimous on its success. There are compilers available for download from the HPC site (http://hpc.sourceforge.net/) that claim to vectorize gcc (version 4.4) in Snow Leopard. I read comments from programmers that full-optimization of the compiler has bugs, so buyer beware. I might try it later, though.
NOTE: In 2012-2013 I switched to compile codes with gfortran compiler.
Wednesday, October 28, 2009
Tuesday, October 20, 2009
grn_sph.f for computing greens function for synthetic seismograms using summed free oscillations for spherical earth models
click title to see Google Drive directory JParkCodes
grn_sph.f
The task of computing synthetic seismograms from summed free oscillations for an aspherical earth model is typically accomplished in two or three steps. The first step computes the hybrid free-oscillations for a chosen earth model, comprised of the "simple" spherical-earth modes whose particle motions follow single vector spherical harmonics. The second step takes those hybrid modes and computes greens functions for particle motion at a given source location (allowing six different time series, one for each component of the seismic moment tensor) and station location (for a given sensor and three time series for three components: vertical/horizontal). The last step is the computation of a synthetic seismogram by summing the six time series for each motion component by weighting with the components of the moment tensor.
Historically, I have used the greens function routine as the location for adding the instrument response, but I have reconsidered this choice. Computers are fast enough for the FFTs an instrument response to be computed on the fly in the last step. This also avoids the problem of updating the instrument response codes (notoriously unstable, owing to the various formats for instrument responses) in multiple greensfunction codes. My final code (sac_syndat.f) is used for all, and should be the place where instrument response either is (or is not) applied. The output of this routine has units of displacement in meters.
So, for making these codes more public, lets start with the simple code for synthetic seismograms on a spherical earth model.
The code is grn_sph.f, which has existed in one form or another for over 20 years (this explains the dead code inside it).
print *,'typical input for this program:'
print *,'0 20 --> freq interval to sum (mHz)'
print *,'/home/jjpark/Mineos/Models/premiss --> 1-D earth model'
print *,'/home/jjpark/Mineos/premissnl --> spheroidal info'
print *,'/home/jjpark/Mineos/premiss_eig --> spheroidal modes'
print *,'/home/jjpark/Mineos/premistnl --> toroidal info'
print *,'/home/jjpark/Mineos/premist_eig --> toroidal modes'
print *,'1000 --> # of overtone branches (here: all of them)'
print *,'0 --> # of data points to synthesize (0-> #=data: NPTS)'
print *,'/park/nobackup/Sumatra/VH/PET.VHZ --> data file'
print *,'.....'
print *,'. . . more data files'
print *,'.....'
print *,'stop'
print *,' '
print *,'The format of the greens fct files is impure SAC-format'
print *,'readable only by sac_syndat and related JPark codes'
print *,'It has a SAC header with one ASCII variable set to GRN'
print *,'and six*NPTS data values that follow the header'
print *,'NPTS for each component of the source moment tensor'
print *,'The greens-function file names are created from your input file'
print *,'by prepending gsph_ to the input filename'
print *,'and replacing the last character with z,r,t'
This code reads files full of free-oscillation displacement functions for a standard earth model, computes the strain of the mode at the desired source depth, and sums 6 traces for surface displacement at a specified distance from the source. Because the earth model is a spherical-earth model, only the source depth and source-receiver distance affect the waveforms that represent oscillatory displacement. The code, however, takes a data file to specify the source and receiver locations and creates synthetics for the geometry. A subsequent processing step (syndat.f) will take the greens function file, sum with the components of a moment tensor, and apply the instrument responses appropriate to the station's sensors and channel.
The start time for the earthquake determines where the decaying sinusoids begin. Older versions of green-fct codes used an explicit time that the user had to supply. This version assumes that the information is in the omarker time in the SAC header. There is also a bmarker time that is used in case the start time of the record does not correspond to the reference time in the data header.
c the reference time of the file could be before or after the first data point
c the ymdhmsec parameters define the reference time
c the origin quake of the quake and start of record are wrt reference time
origin=ahead(8) ! event start is the omarker in sac header
recstart=ahead(6) ! record start is the bmarker in sac header
grn_sph.f
The task of computing synthetic seismograms from summed free oscillations for an aspherical earth model is typically accomplished in two or three steps. The first step computes the hybrid free-oscillations for a chosen earth model, comprised of the "simple" spherical-earth modes whose particle motions follow single vector spherical harmonics. The second step takes those hybrid modes and computes greens functions for particle motion at a given source location (allowing six different time series, one for each component of the seismic moment tensor) and station location (for a given sensor and three time series for three components: vertical/horizontal). The last step is the computation of a synthetic seismogram by summing the six time series for each motion component by weighting with the components of the moment tensor.
Historically, I have used the greens function routine as the location for adding the instrument response, but I have reconsidered this choice. Computers are fast enough for the FFTs an instrument response to be computed on the fly in the last step. This also avoids the problem of updating the instrument response codes (notoriously unstable, owing to the various formats for instrument responses) in multiple greensfunction codes. My final code (sac_syndat.f) is used for all, and should be the place where instrument response either is (or is not) applied. The output of this routine has units of displacement in meters.
So, for making these codes more public, lets start with the simple code for synthetic seismograms on a spherical earth model.
The code is grn_sph.f, which has existed in one form or another for over 20 years (this explains the dead code inside it).
print *,'typical input for this program:'
print *,'0 20 --> freq interval to sum (mHz)'
print *,'/home/jjpark/Mineos/Models/premiss --> 1-D earth model'
print *,'/home/jjpark/Mineos/premissnl --> spheroidal info'
print *,'/home/jjpark/Mineos/premiss_eig --> spheroidal modes'
print *,'/home/jjpark/Mineos/premistnl --> toroidal info'
print *,'/home/jjpark/Mineos/premist_eig --> toroidal modes'
print *,'1000 --> # of overtone branches (here: all of them)'
print *,'0 --> # of data points to synthesize (0-> #=data: NPTS)'
print *,'/park/nobackup/Sumatra/VH/PET.VHZ --> data file'
print *,'.....'
print *,'. . . more data files'
print *,'.....'
print *,'stop'
print *,' '
print *,'The format of the greens fct files is impure SAC-format'
print *,'readable only by sac_syndat and related JPark codes'
print *,'It has a SAC header with one ASCII variable set to GRN'
print *,'and six*NPTS data values that follow the header'
print *,'NPTS for each component of the source moment tensor'
print *,'The greens-function file names are created from your input file'
print *,'by prepending gsph_ to the input filename'
print *,'and replacing the last character with z,r,t'
This code reads files full of free-oscillation displacement functions for a standard earth model, computes the strain of the mode at the desired source depth, and sums 6 traces for surface displacement at a specified distance from the source. Because the earth model is a spherical-earth model, only the source depth and source-receiver distance affect the waveforms that represent oscillatory displacement. The code, however, takes a data file to specify the source and receiver locations and creates synthetics for the geometry. A subsequent processing step (syndat.f) will take the greens function file, sum with the components of a moment tensor, and apply the instrument responses appropriate to the station's sensors and channel.
The start time for the earthquake determines where the decaying sinusoids begin. Older versions of green-fct codes used an explicit time that the user had to supply. This version assumes that the information is in the omarker time in the SAC header. There is also a bmarker time that is used in case the start time of the record does not correspond to the reference time in the data header.
c the reference time of the file could be before or after the first data point
c the ymdhmsec parameters define the reference time
c the origin quake of the quake and start of record are wrt reference time
origin=ahead(8) ! event start is the omarker in sac header
recstart=ahead(6) ! record start is the bmarker in sac header
Friday, August 1, 2008
Rayleigh Wave Dispersion Code chenray.f
click title to see Google Drive directory JParkCodes
This a a computer code that performs the 1D layered-media reflectivity method to compute Rayleigh-wave dispersion. Was an exercise for me when adapting the computational scheme of the GJI paper by Xiaofei Chen (1993) for layered anisotropic media. Uses model files in the same format used by the anisotropic surface-wave code aniprop.f and the anisotropic body-wave codes anirec.f and related. The read statements read the anisotropic parameters into dummy variables and ignore them. The code only plots dispersion curves for phase velocity CVEL, but does not compute GVEL, the group velocity. The values of CVEL are *not* written out for later use. You can modify the code yourself to write out cvel values.
c chenray - to calculate propagating Rayleigh modes of layered structure
c revised to calculate dispersion curves 7/7/95, debugged 7/27/08
c 9/2/94 - cannibalized from aniprop to test its structure
c xf77 -o /home/jjpark/bin/chenray chenray.f /home/jjpark/Plotxy/plotlib.a /home/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/chenray chenray.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/jlib.a
c
c will calculate Rayleigh wave dispersion & eigenfunctions using formulas
c of Chen (1993)
c max is 100 layers over a halfspace
c model is read from file 'model'
c reads nl=#_of_layers, then depth_i,rho_i,alpha_i,beta_i, i=1,nl+1
c ee, e1 and e2 are 4x4 matrices that get partitioned into 2x2 submatricies to
c separate out the reflection and transmission of waves at an interface.
c
c rt(j) is generalized reflection coef \tilde{R}_du at jth interface
c tt(j) is generalized transmission coef \tilde{T}_d at jth interface
c rt0 is gnrlzd refl coef \tilde{R}_ud^(0) at free surface (0th interface)
c rtm(*,*,j) is matrix of modified R/T coefs at the jth interface
c rtm0(2,2) is matrix of modified R/T coefs at free surface - 0th interface
c
c OUTPUT: makes plots of dispersion, but does not write phase velocity (cvel) to a file
c or compute the group velocity (gvel). The code has not worked in 1994, and I chose to
c push on to the general anisotropic surface wave code aniprop.f, rather than to debug this
c Rayleigh-wave code. This code was revived and debugged in 2008(!) by JPark
c
This a a computer code that performs the 1D layered-media reflectivity method to compute Rayleigh-wave dispersion. Was an exercise for me when adapting the computational scheme of the GJI paper by Xiaofei Chen (1993) for layered anisotropic media. Uses model files in the same format used by the anisotropic surface-wave code aniprop.f and the anisotropic body-wave codes anirec.f and related. The read statements read the anisotropic parameters into dummy variables and ignore them. The code only plots dispersion curves for phase velocity CVEL, but does not compute GVEL, the group velocity. The values of CVEL are *not* written out for later use. You can modify the code yourself to write out cvel values.
c chenray - to calculate propagating Rayleigh modes of layered structure
c revised to calculate dispersion curves 7/7/95, debugged 7/27/08
c 9/2/94 - cannibalized from aniprop to test its structure
c xf77 -o /home/jjpark/bin/chenray chenray.f /home/jjpark/Plotxy/plotlib.a /home/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/chenray chenray.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/jlib.a
c
c will calculate Rayleigh wave dispersion & eigenfunctions using formulas
c of Chen (1993)
c max is 100 layers over a halfspace
c model is read from file 'model'
c reads nl=#_of_layers, then depth_i,rho_i,alpha_i,beta_i, i=1,nl+1
c ee, e1 and e2 are 4x4 matrices that get partitioned into 2x2 submatricies to
c separate out the reflection and transmission of waves at an interface.
c
c rt(j) is generalized reflection coef \tilde{R}_du at jth interface
c tt(j) is generalized transmission coef \tilde{T}_d at jth interface
c rt0 is gnrlzd refl coef \tilde{R}_ud^(0) at free surface (0th interface)
c rtm(*,*,j) is matrix of modified R/T coefs at the jth interface
c rtm0(2,2) is matrix of modified R/T coefs at free surface - 0th interface
c
c OUTPUT: makes plots of dispersion, but does not write phase velocity (cvel) to a file
c or compute the group velocity (gvel). The code has not worked in 1994, and I chose to
c push on to the general anisotropic surface wave code aniprop.f, rather than to debug this
c Rayleigh-wave code. This code was revived and debugged in 2008(!) by JPark
c
Love-Wave Dispersion Code: chenlove.f
click title to see Google Drive directory JParkCodes
chenlove.f
This a a computer code that performs the 1D layered-media reflectivity method to compute Love-wave dispersion. Was an exercise for me when adapting the computational scheme of the GJI paper by Xiaofei Chen (1993) for layered anisotropic media. Uses model files in the same format used by the anisotropic surface-wave code aniprop.f and the anisotropic body-wave codes anirec.f and related. The read statements read the anisotropic parameters into dummy variables and ignore them.
c will calculate love wave dispersion & eigenfunctions using formulas
c of Chen (1993)
c max is 100 layers over a halfspace
c model is read from file 'model'
c reads nl=#_of_layers, then depth_i,rho_i,alpha_i,beta_i, i=1,nl+1
c the variables:
c c(j,1) are C_d; c(j,2) are the C_u from Chen (1993) for the jth layer
c c(nlayer+1,*) are coefs in the halfspace.
c rt(j) is generalized reflection coef \tilde{R}_du at jth interface
c tt(j) is generalized transmission coef \tilde{T}_d at jth interface
c rt0 is gnrlzd refl coef \tilde{R}_ud^(0) at free surface (0th interface)
c rtm(*,*,j) is matrix of modified R/T coefs at the jth interface
c rtm0(2,2) is matrix of modified R/T coefs at free surface - 0th interface
c
c OUTPUT: makes plots of dispersion, and also writes phase velocity (cvel)
c group velocity (gvel) and three integration scalars that are used in the variational
c computation of group velocity from phase velocity
c do the sum
c ...
c ci1=ci1+rho(il)*(dci1+dci2)
c ci2=ci2+xmu(il)*(dci1+dci2)
c ci4=ci4+xmu(il)*(dci1-dci2)*xnu(il)**2
c end do
c gvel=ci2/(cvel*ci1)
c
c
chenlove.f
This a a computer code that performs the 1D layered-media reflectivity method to compute Love-wave dispersion. Was an exercise for me when adapting the computational scheme of the GJI paper by Xiaofei Chen (1993) for layered anisotropic media. Uses model files in the same format used by the anisotropic surface-wave code aniprop.f and the anisotropic body-wave codes anirec.f and related. The read statements read the anisotropic parameters into dummy variables and ignore them.
c will calculate love wave dispersion & eigenfunctions using formulas
c of Chen (1993)
c max is 100 layers over a halfspace
c model is read from file 'model'
c reads nl=#_of_layers, then depth_i,rho_i,alpha_i,beta_i, i=1,nl+1
c the variables:
c c(j,1) are C_d; c(j,2) are the C_u from Chen (1993) for the jth layer
c c(nlayer+1,*) are coefs in the halfspace.
c rt(j) is generalized reflection coef \tilde{R}_du at jth interface
c tt(j) is generalized transmission coef \tilde{T}_d at jth interface
c rt0 is gnrlzd refl coef \tilde{R}_ud^(0) at free surface (0th interface)
c rtm(*,*,j) is matrix of modified R/T coefs at the jth interface
c rtm0(2,2) is matrix of modified R/T coefs at free surface - 0th interface
c
c OUTPUT: makes plots of dispersion, and also writes phase velocity (cvel)
c group velocity (gvel) and three integration scalars that are used in the variational
c computation of group velocity from phase velocity
c do the sum
c ...
c ci1=ci1+rho(il)*(dci1+dci2)
c ci2=ci2+xmu(il)*(dci1+dci2)
c ci4=ci4+xmu(il)*(dci1-dci2)*xnu(il)**2
c end do
c gvel=ci2/(cvel*ci1)
c
c
Tuesday, July 1, 2008
agu.f Viewing the output of the sacw_corr_many code
click title to see Google Drive directory JParkCodes
The code agu.f (guess where I wrote it!) reads the cvel_YYYY.DDD.HH files and generates some plots of the phase-velocity as fct of period using a collection of algorithms: weighted sum, trimmed weighted sum, median. Version posted can process up to 100 cvel files. You can expand the arrays for more.
input:
cvel_2006.175.02 --> input file (one per event , just "ls cvel_20* >! c_agu"
cvel_2006.225.10
cvel_2006.225.19
cvel_2006.240.22
stop --------------------> end of input files
3 --------------------> make a plot
4 --------------------> make a plot and write to PS file
cvel_eastwest.ps ---------> PS plotfile
The code agu.f (guess where I wrote it!) reads the cvel_YYYY.DDD.HH files and generates some plots of the phase-velocity as fct of period using a collection of algorithms: weighted sum, trimmed weighted sum, median. Version posted can process up to 100 cvel files. You can expand the arrays for more.
input:
cvel_2006.175.02 --> input file (one per event , just "ls cvel_20* >! c_agu"
cvel_2006.225.10
cvel_2006.225.19
cvel_2006.240.22
stop --------------------> end of input files
3 --------------------> make a plot
4 --------------------> make a plot and write to PS file
cvel_eastwest.ps ---------> PS plotfile
sacw_corr_many.f -- Using the wavelet-based surface wave code
click title to see Google Drive directory JParkCodes
sacw_corr_many.f does a 2-station estimation of phase velocity.
sacw --> sac-format data, analyzed with wavelets
corr --> cross-correlation is used to estimated phase shits and phase velocity
many --> code is set up to run a series of data-record pairs.
The posted code is designed for 10 sps data, and computes Slepian wavelets in 16 banks logarithmically spaced between 8 and 56 second center-period. Three complex-valued Slepian wavelets are used, in order to estimate phased cross-correlations with 2 and 4 degrees of freedom. The phase of the cross-correlation is scaled with the center-period of the wavelet to compute a phase-delay. This phase delay is computed at evenly spaced points in the seismogram. The station locations and the back-azimuth of the earthquake are used to estimate a phase velocity from the phase delay, but only at the time point where the maximum wavelet spectrum is found. The wavelet spectrum is computed only from the first Slepian wavelet, because experience demonstrates that this affords a more stable determination of the max-amplitude arrival time.
The input is at the terminal, so the command
sacw_corr_many < c_files_process
will be standard analysis method. This command file has the form (for stations PIIR and VOLR)
2003.325.04.12.54.PIIR.DHZ
2003.325.04.12.56.VOLR.DHZ
3 ---> plot to screen, could be "0" to skip plotting
4 ---> plot to screen and write output to PS file (also could skip plotting with "0")
2003.325.ps ---> name of the PS file for phase-velocity graph
2004.024.13.05.37.PIIR.DHZ
2004.024.13.05.40.VOLR.DHZ
3
4
2004.024.ps
stop ----> end of data input
The data files must have a common start time (despite the differing numbers in the filenames in this example) In practice, one can use SAC commands to achieve this easily. It is convenient to drop the minutes and seconds entries in the filenames, saving the year/julianday/hour filename format, because these parameters will match for different stations, and the minutes/seconds wont.. For instance, if you have a file "infiles" of the filename pairs to analyze, then
ls 20*VOLR*DHZ | sed 's/\./ /g' | awk '{print "read " $1 "." $2 "." $3 "*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off"}' >! c_synch
will generate the SAC macro to synchronize and cut to a common start point:
read 2003.325.04*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.009.19*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.024.13*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.038.21*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
and so on . . . .
The code writes a small ASCII file of phase-vel information, cvel_YYYY.DDD.HH, which looks like this
cat cvel_2004.339.10
56.00 1.547 1.858 13.76 236.6 341.4 7.340
49.18 4.596 12.74 13.76 236.6 342.0 2.470
43.20 2.893 1.626 13.76 236.6 319.2 3.925
37.93 4.331 1.237 13.76 236.6 321.0 2.621
33.32 3.411 0.3311 13.76 236.6 322.8 3.328
29.27 3.258 0.9897E-01 13.76 236.6 334.8 3.484
25.70 2.934 0.2715 13.76 236.6 334.8 3.869
22.58 2.948 0.1586 13.76 236.6 347.4 3.850
19.83 2.939 0.2432 13.76 236.6 358.2 3.863
17.42 2.759 1.053 13.76 236.6 372.0 4.114
15.30 3.739 2.711 13.76 236.6 370.2 3.036
13.43 5.761 2.135 13.76 236.6 417.6 1.971
11.80 4.417 0.4976 13.76 236.6 456.0 2.570
10.37 1.556 0.1375 13.76 236.6 559.2 7.295
9.100 2.039 0.4498 13.76 236.6 725.4 5.568
8.000 -27.40 99.12 13.76 236.6 995.4 -0.4143
period cvel d(cvel) gcarc baz timeOcorrelation delay-time
sacw_corr_many.f does a 2-station estimation of phase velocity.
sacw --> sac-format data, analyzed with wavelets
corr --> cross-correlation is used to estimated phase shits and phase velocity
many --> code is set up to run a series of data-record pairs.
The posted code is designed for 10 sps data, and computes Slepian wavelets in 16 banks logarithmically spaced between 8 and 56 second center-period. Three complex-valued Slepian wavelets are used, in order to estimate phased cross-correlations with 2 and 4 degrees of freedom. The phase of the cross-correlation is scaled with the center-period of the wavelet to compute a phase-delay. This phase delay is computed at evenly spaced points in the seismogram. The station locations and the back-azimuth of the earthquake are used to estimate a phase velocity from the phase delay, but only at the time point where the maximum wavelet spectrum is found. The wavelet spectrum is computed only from the first Slepian wavelet, because experience demonstrates that this affords a more stable determination of the max-amplitude arrival time.
The input is at the terminal, so the command
sacw_corr_many < c_files_process
will be standard analysis method. This command file has the form (for stations PIIR and VOLR)
2003.325.04.12.54.PIIR.DHZ
2003.325.04.12.56.VOLR.DHZ
3 ---> plot to screen, could be "0" to skip plotting
4 ---> plot to screen and write output to PS file (also could skip plotting with "0")
2003.325.ps ---> name of the PS file for phase-velocity graph
2004.024.13.05.37.PIIR.DHZ
2004.024.13.05.40.VOLR.DHZ
3
4
2004.024.ps
stop ----> end of data input
The data files must have a common start time (despite the differing numbers in the filenames in this example) In practice, one can use SAC commands to achieve this easily. It is convenient to drop the minutes and seconds entries in the filenames, saving the year/julianday/hour filename format, because these parameters will match for different stations, and the minutes/seconds wont.. For instance, if you have a file "infiles" of the filename pairs to analyze, then
ls 20*VOLR*DHZ | sed 's/\./ /g' | awk '{print "read " $1 "." $2 "." $3 "*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off"}' >! c_synch
will generate the SAC macro to synchronize and cut to a common start point:
read 2003.325.04*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.009.19*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.024.13*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
read 2004.038.21*HZ;synch;write over;cuterr fillz;cut 0 e; read; write over; cut off
and so on . . . .
The code writes a small ASCII file of phase-vel information, cvel_YYYY.DDD.HH, which looks like this
cat cvel_2004.339.10
56.00 1.547 1.858 13.76 236.6 341.4 7.340
49.18 4.596 12.74 13.76 236.6 342.0 2.470
43.20 2.893 1.626 13.76 236.6 319.2 3.925
37.93 4.331 1.237 13.76 236.6 321.0 2.621
33.32 3.411 0.3311 13.76 236.6 322.8 3.328
29.27 3.258 0.9897E-01 13.76 236.6 334.8 3.484
25.70 2.934 0.2715 13.76 236.6 334.8 3.869
22.58 2.948 0.1586 13.76 236.6 347.4 3.850
19.83 2.939 0.2432 13.76 236.6 358.2 3.863
17.42 2.759 1.053 13.76 236.6 372.0 4.114
15.30 3.739 2.711 13.76 236.6 370.2 3.036
13.43 5.761 2.135 13.76 236.6 417.6 1.971
11.80 4.417 0.4976 13.76 236.6 456.0 2.570
10.37 1.556 0.1375 13.76 236.6 559.2 7.295
9.100 2.039 0.4498 13.76 236.6 725.4 5.568
8.000 -27.40 99.12 13.76 236.6 995.4 -0.4143
period cvel d(cvel) gcarc baz timeOcorrelation delay-time
Code for surface-wave phase-velocity estimation sacw_corr_many.f
click title to see Google Drive directory JParkCodes
c program sacw_corr_many
c JJP 04/13/07 updated 7/24/07 updated 12/11/07, again 6/25/08
c compute phase delays as function of frequency using the multiwavelet transform
c input two fil,es, typically BHZ components at different stations
c cross-correlate the wavelet transforms
c output is a colormap of time and frequency dependent time delay
c plotted with the wavelet transform of the two time series
c
c f77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c f77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr.f /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr_many.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c to read SAC-format data, filename.bh? <-- br="" e="" for="" n="" routine="" subs="" z="">c
c generates wavelets with eigenvector decomposition + interpolation
c calculate padded FFT of data for convolution with wavelet FFTs
c general use - calculates set of wavelets for a specified number of points
c interpolates and fft other lengths on the fly
c max=12 real-wavelets max = 6 complex-valued wavelets
c
c code is based on the algorithms described in
c
c Lilly, J., and J. Park, Multiwavelet spectral and polarization analysis
c of seismic records, Geophysical J. International, v122, 1001-1021, 1995.
c
c Park, J. and M. E. Mann, Interannual temperature events and shifts in
c global temperature: A multiwavelet correlation approach, Earth Interactions,
c v. 4, online paper 4-001, 2000.
c
C SEE ALSO:
c Bear, L. K., and G. L. Pavlis, Estimation of slowness vectors and
c their uncertainties using multi-wavelet seismic array processing,
c Bull. Seismol. Soc. Am., v87, 755--769, 1997.
c
c this code generates 2-D arrays of spectra, polarization and statistical stuff
c it writes GMT scripts to display these arrays in dazzling color
c
c for an input filename sacfile.bh? the GMT scripts are written to
c GMD_sacfile.bh
c
c
c the ASCII files to plot are named
c cpt_nilcut.sh spw_nilcut.sh1 spw_nilcut.sh2
c dat_nilcut.sh
c cpd_nilcut.sh del_nilcut.sh
c
c dat_* is data to plot
c cpt_* and cohere.cpt are colorfiles for GMT
c del_* dep_* and dem_* are the time delay arrays
c spw_* wavelet spectra
c
c you can fill a lot of disk space if you are not tidy with old files!
c
c user must make scripts executable ("chmod +x GM*_sacfile.bh"),
c execute the scripts and display the postscript files
c program sacw_corr_many
c JJP 04/13/07 updated 7/24/07 updated 12/11/07, again 6/25/08
c compute phase delays as function of frequency using the multiwavelet transform
c input two fil,es, typically BHZ components at different stations
c cross-correlate the wavelet transforms
c output is a colormap of time and frequency dependent time delay
c plotted with the wavelet transform of the two time series
c
c f77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /park/backup/bin/sacw_corr_many sacw_corr.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c f77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr.f /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/sacw_corr_many sacw_corr_many.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c to read SAC-format data, filename.bh? <-- br="" e="" for="" n="" routine="" subs="" z="">c
c generates wavelets with eigenvector decomposition + interpolation
c calculate padded FFT of data for convolution with wavelet FFTs
c general use - calculates set of wavelets for a specified number of points
c interpolates and fft other lengths on the fly
c max=12 real-wavelets max = 6 complex-valued wavelets
c
c code is based on the algorithms described in
c
c Lilly, J., and J. Park, Multiwavelet spectral and polarization analysis
c of seismic records, Geophysical J. International, v122, 1001-1021, 1995.
c
c Park, J. and M. E. Mann, Interannual temperature events and shifts in
c global temperature: A multiwavelet correlation approach, Earth Interactions,
c v. 4, online paper 4-001, 2000.
c
C SEE ALSO:
c Bear, L. K., and G. L. Pavlis, Estimation of slowness vectors and
c their uncertainties using multi-wavelet seismic array processing,
c Bull. Seismol. Soc. Am., v87, 755--769, 1997.
c
c this code generates 2-D arrays of spectra, polarization and statistical stuff
c it writes GMT scripts to display these arrays in dazzling color
c
c for an input filename sacfile.bh? the GMT scripts are written to
c GMD_sacfile.bh
c
c
c the ASCII files to plot are named
c cpt_nilcut.sh spw_nilcut.sh1 spw_nilcut.sh2
c dat_nilcut.sh
c cpd_nilcut.sh del_nilcut.sh
c
c dat_* is data to plot
c cpt_* and cohere.cpt are colorfiles for GMT
c del_* dep_* and dem_* are the time delay arrays
c spw_* wavelet spectra
c
c you can fill a lot of disk space if you are not tidy with old files!
c
c user must make scripts executable ("chmod +x GM*_sacfile.bh"),
c execute the scripts and display the postscript files
Subscribe to:
Posts (Atom)