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
Friday, August 1, 2008
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
Tuesday, May 27, 2008
changes in the inverse FFT in recfunk codes
click title to see Google Drive directory JParkCodes
After some comparison runs, Vadim Levin convinced me that the "improvement" I made to the inverse FFT in the recfunk*, rfmig* codes was a bad innovation. The original RF-estimation codes applied a cos^2 factor to the frequency-domain RFs H_R(f) and H_T(f) with weighting unity at the f=0 point, decreasing to zero at f_max while following a cos-squared dependence. I changed all codes to have a sharper bandpass (I was unhappy that the half-amplitude point occurred at f_max/2!) with unit factor from f=0 to f=0.5*f_max, then a cos-squared rolloff within (0.5*f_max , f_max). This sharper bandpass led to small negative sidelobes at the flanks of sharp Ps pulses. Although one can live with this, the temptation to interpret a wiggle could be strong, so it is safer to use the original weighting of H(f) before inverse FFT to obtain the time-domain RF. The sidelobes do not disappear completely, but are far smaller.
After some comparison runs, Vadim Levin convinced me that the "improvement" I made to the inverse FFT in the recfunk*, rfmig* codes was a bad innovation. The original RF-estimation codes applied a cos^2 factor to the frequency-domain RFs H_R(f) and H_T(f) with weighting unity at the f=0 point, decreasing to zero at f_max while following a cos-squared dependence. I changed all codes to have a sharper bandpass (I was unhappy that the half-amplitude point occurred at f_max/2!) with unit factor from f=0 to f=0.5*f_max, then a cos-squared rolloff within (0.5*f_max , f_max). This sharper bandpass led to small negative sidelobes at the flanks of sharp Ps pulses. Although one can live with this, the temptation to interpret a wiggle could be strong, so it is safer to use the original weighting of H(f) before inverse FFT to obtain the time-domain RF. The sidelobes do not disappear completely, but are far smaller.
recfunk_pick.f
recfunk_estack.f
rfmigrate_estack.f
rfmigrate.f
recfunk_svd.f
recfunk_mwm.f
recfunk_mwmj.f
recfunk_mwms.f
rfmig_boot.f
rfmig_mboot.f
rfmig_cboot.f
rfmig_mcboot.f
Update of anirec_synth.f and anirec_synth_circle.f
click title to see Google Drive directory JParkCodes
A small amount of random noise was added to synthetics to prevent zero-divide errors (in cases where the recfunk codes tries to compute an RF for an interval of x=0 data). In the update, this bit of noise was made smaller, and explicitly zero-mean. There were some low-frequency artifacts cropping up in synthetic RFs. The updated versions seem to avoid this problem:
anirec_synth.f
anirec_synth_circle.f
A small amount of random noise was added to synthetics to prevent zero-divide errors (in cases where the recfunk codes tries to compute an RF for an interval of x=0 data). In the update, this bit of noise was made smaller, and explicitly zero-mean. There were some low-frequency artifacts cropping up in synthetic RFs. The updated versions seem to avoid this problem:
anirec_synth.f
anirec_synth_circle.f
for info see link
Monday, March 24, 2008
The erstwhile version of aniprop.f has had a bugfix
click title to see Google Drive directory JParkCodes
The erstwhile version of aniprop.f that was (is?) available on my website has been superceded by an all-in-one version that computes dispersion curves and synthetic seismograms all in one program. Adapted and debugged on an Intel Mac in gnu fortran.
Google Drive, directory JParkCodes/aniprop_072307.f
after a bugfix is now
Google Drive, directory JParkCodes/aniprop_032408.f
The code computes surface wave modes at evenly-spaced frequency points between 0 and 0.5 Hz. You can change this by tinkering with the code. It computes synthetics for a source at a specified downrange distance and depth, and plots them using PLOTIT. The dispersion values are output to files out_cvel and out_gvel (phase and group velocities, natch), and the synthetics are written out to SAC-format files.
The bugfix concerned a LU-decomposition with partial pivoting that could fail in cases where the radial and transverse components of motion in a layered structure were totally uncoupled e.g., isotropic media, or media in which the axis of symmetry for anisotropy had strike angle phi=0.
The erstwhile version of aniprop.f that was (is?) available on my website has been superceded by an all-in-one version that computes dispersion curves and synthetic seismograms all in one program. Adapted and debugged on an Intel Mac in gnu fortran.
Google Drive, directory JParkCodes/aniprop_072307.f
after a bugfix is now
Google Drive, directory JParkCodes/aniprop_032408.f
The code computes surface wave modes at evenly-spaced frequency points between 0 and 0.5 Hz. You can change this by tinkering with the code. It computes synthetics for a source at a specified downrange distance and depth, and plots them using PLOTIT. The dispersion values are output to files out_cvel and out_gvel (phase and group velocities, natch), and the synthetics are written out to SAC-format files.
The bugfix concerned a LU-decomposition with partial pivoting that could fail in cases where the radial and transverse components of motion in a layered structure were totally uncoupled e.g., isotropic media, or media in which the axis of symmetry for anisotropy had strike angle phi=0.
Subscribe to:
Posts (Atom)