JParkCodes

GMT script to generate Figure 4 of Park and Royer (2011)

2011-06-14T09:58:00.001-04:00

#!/bin/csh

gmtset ANNOT_FONT_SIZE 12p LABEL_FONT_SIZE 15p HEADER_FONT_SIZE 15p

gmtset ANNOT_FONT Times-Roman LABEL_FONT Times-Roman HEADER_FONT Times-Roman

awk '{print $1, $2, $3/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.025/0.135 -I0.4/0.01 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.02:'ACT'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X1.0 -Y1.0 -K >! aaplot

awk '{print $1, $5, $6/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.15/0.85 -I0.4/0.05 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.1:'FERT'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X0.0 -Y3.5 -K -O >> aaplot

awk '{print $1, $8, $9/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.05/0.55 -I0.4/0.025 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.1:'LIFE'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X4.0 -Y0.0 -K -O >> aaplot

awk '{print $1, $11, $12/1000}' out_ppdf_control.dat >! blob.dat

surface blob.dat -R0.2/10.6/0.45/1.25 -I0.4/0.05 -Gblobg -T0.25

grdimage blobg -JX3/2 -Cppdf.cpt -R -Ba2f1:'CO@-2@- Sensitivity (degC)':/0.2:'GYM'::'.Fraction of ModelRuns within 3-@~s@~':/SWen -X0.0 -Y-3.5 -K -O >> aaplot

psscale -D3.5/1/2/0.25c -Cppdf.cpt -Ba0.2f0.1:'.Fraction <3-@~s@~': -O >> aaplot

gv aaplot &

CO2 proxy data set for Park and Royer (2011)

2011-06-14T09:35:00.003-04:00

http://earth.geology.yale.edu/~jjpark/Code/co2proxy_2010.dat

carbon dioxide proxy data set used by

Park, J., and D. L. Royer, Geologic constraints on the glacial amplification of Phanerozoic climate sensitivity, American J. Science, 311, 1-26, doi:10.2475/01.2011.01, 2011.

The sources of data are enumerated in the reference list. The dataset is an updated version of the data compilation reported by

Royer DL. 2006. CO2-forced climate thresholds during the Phanerozoic. Geochimica et Cosmochimica Acta, 70: 5665-5675.

CO2 Proxy data set from Royer et al 2007

2011-06-14T08:29:00.003-04:00

Carbon dioxide proxy data set for the Phanerozoic used by

Royer, D. L., R. A. Berner, and J. Park, Climate sensitivity constrained by carbon dioxide concentrations over the last 420 million years, Nature, v446, 530-532, 2007.

data compilation described in

Royer DL. 2006. CO2-forced climate thresholds during the Phanerozoic. Geochimica et Cosmochimica Acta, 70: 5665-5675.

post-processing of GEOCARBSULFvolc simulation results

2011-06-14T07:57:00.000-04:00

c program gcsv_ppdf.f

c 5/30/10 JJP

c f77 -o /Users/jjpark/bin/gcsv_ppdf gcsv_ppdf.f

c program to read variance reduction matrices for gcsfdata_export/search

c reads relative variance residual for CO2-proxy datafit to ndat*10000

c geocarb model runs

GEOCARBSULFvolc carbon-cycle simulation code

2011-06-14T07:50:00.003-04:00

http://earth.geology.yale.edu/~jjpark/Code/gcsv10_export.f

c program gcsv10_export.f

c fortran version of bob berner's GEOCARBSULF BASIC code

c gcsv10 stands for GeoCarbSulfVolc algorithm for 2010

c THIS CODE CAN READ EITHER THE 2006/7 proxy-CO2 data file

c OR THE 2010 DATA FILE.

c 8/17/2010 -- confirmed the VNV parameter as 4, though Berner asked for 5,

c as motivated by Aaron Taylors thesis -- but other Taylor source is vnv=3

c this code developed in MacOS 10.4, using gcc compiler,

c takes 90minutes to run on a 2006-vintage Macbook Pro with 2 GHz Intel Core Duo processor

c takes 45minutes to run on a 2009-vintage Macbook Pro with 2.53 GHz Intel Core 2 Duo processor

c -- looped over range of input parameters, set up to fit proxy data for CO2 ppm

c 11/22/06 JJP

c proxy data has uncertainties that are sometimes large. To make the data-fitting robust,

c compute chi-squared using an uncertainty in the log-domain.

c xf77 -o /Users/jjpark/bin/gcsv10_export gcsv10_export.f /Users/jjpark/Plotxy/plotlib.a

c code uses plotit, a standalone fortran/C plotting subroutine in Xwindows

c this code can be found on JPark's website http://earth.geology.yale.edu/~jjpark/Code/plotit.tar

c instructions for PLOTIT at http://earth.geology.yale.edu/~jjpark/Code/plotit.html

c code updated from the code used for the paper

c Royer, D. L., Berner, R. A., and Park, J., 2007, Climate sensitivity constrained by carbon

c dioxide concentrations over the last 420 million years: Nature, v.446, p.530-532.

c and used for this paper

c Park, J., and D. L. Royer, Geologic constraints on the glacial amplification of Phanerozoic

c climate sensitivity, American J. Science, V.311, 2011, P. 1-26, DOI 10.2475/01.2011.01

c there are many test computations in Park and Royer (2011) that were effected

c with modified versions of this gcsv10_export.f code

c Output of this code is written to two files

c outoutout.dat is a simple table of CO2proxy data values averaged in 10My intervals

c out_gcsv10.dat is a large ASCII file (>50 Mb) that contains the GEOCARB carbon-cycle simulations

c for all parameter choices. This allows the user to examine the statistical behavior of the

c carbon-cycle simulations offline this code (e.g. in gcsv_ppdf.f)

c code converted from BASIC to run grid search over parameters

c Newton-Raphson damping applied here, to avoid NaN.

c divergent transition over 380-350My is interpolated better,

c correcting earlier convergence error

c the code reads a file proxydata.txt that contains proxy CO2 PPM data

c the GEOCARBSULF carbon-cycle model is computed at time steps of 1-My from 570Ma to present

c the proxy data ranges only from 420Ma to the present

c Five nested loops comprise the heart of the code, from outermost to innermost

c DeltaT = temp increase with CO2 doubling

c ACT = dimensionless activation-energy for silicate weathering

c FERT = plant CO2-fertilization coefficient

c LIFE = liverwort-based weathering factor -- less than vascular plants

c GYM = gymnosperm-based weathering factor (1=angiosperm weathering parameter)

c gymnosperms (e.g. conifers) might be better at weathering than angiosperms

c for each choice of parameters, we compute the carbon flux balance between terms

c that depend on silicate weathering (fBBS) and carbon burial/degassing (fB),

c each normalized with common factors derived by Bob Berner many years ago

c The carbon-cycle balance depends on CO2 level,

c because the weathering of silicates accelerates with warmer temps induced by greenhouse gases

c the code solves for the CO2 level that achieves the balance of carbon fluxes

c The CO2-dependent parameter fB is calculated at each time step from carbon and

c carbon isotope mass balance and values of all other parameters, based on

c geological and biological data, that affect the carbon cycle. Then the value of

c RCO2 (CO2 concentration normalized to the mean value for the past 1 million

c years = 250 ppm) is calculated by inversion from a complex expression for fB based on the

c greenhouse effect, CO2 fertilization of plant weathering, solar evolution, and

c changes in land temperature).

c for detailed exposition of the theoretical GEOCARB models, see

c Berner, R. A., The Phanerozoic Carbon Cycle: CO2 and O2, 150pp. Oxford University Press, 2004.

c and a collection of GEOCARB and GEOCARBSULF papers published over the past decade by

c Berner and Colleagues.

c Berner, R. A., 2006a, GEOCARBSULF: A combined model for Phanerozoic atmospheric O2 and CO2:

c Geochimica et Cosmochimica Acta, v.70, p.5653-5664.

c Berner, R. A., 2006b, Inclusion of the weathering of volcanic rocks in the GEOCARBSULF model:

c American Journal of Science, v.306, p.295-302, doi:10.2475/052006.01.

c Berner, R. A., 2008, Addendum to Inclusion of the weathering of volcanic rocks in the

c GEOCARBSULF model (Berner, R. A., 2006, v.306, p.295-302) American Journal of Science,

c v.308, p.100-103, doi:10.2475/01.2008.04.

c Berner, R. A., 2009, Phanerozoic atmospheric oxygen: New results using the GEOCARBSULF Model:

c American Journal of Science, v.309, p.603-606, doi:10.2475/07.2009.03.

c The code computes the aggregate data fit of GEOCARBSULF for each choice

c of the combined parameters, expressed as the fractional variance residual when expressed

c in the log domain without uncertainty weighting.

c description of the proxy data can be found in

c Royer et al, CO2 as a primary driver of Phanerozoic climate, GSA Today, March 2004.

Reminder for the Fortran plotting subroutine PLOTIT

2011-06-14T07:36:00.001-04:00

http://earth.geology.yale.edu/~jjpark/Code/plotit.tar

information README at

http://earth.geology.yale.edu/~jjpark/Code/plotit.html

Disclaimer and Free Software Statement

2009-11-11T21:22:00.002-05:00

The computer codes that I post to the JParkCodes blog are provided as free software.

You can redistribute the codes and/or modify them under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,

but WITHOUT ANY WARRANTY; without even the implied warranty of

MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

GNU General Public License for more details.

You may obtain a copy of the GNU General Public License

along with these program at the link above. If not active, see http://www.gnu.org/licenses/.

grn_propzon.f -- code to compute greensfunctions for zonally-symmetric earth models.

2009-11-05T11:32:00.003-05:00

c program grn_propzon -- three-component gfcts
c for zonally-symmetric earth models
c f77 -o /Users/jjpark/bin/grn_propzon grn_propzon.f /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/grn_propzon grn_propzon.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/jlib.a
c the code reads a special compressed hybrid eigenvector file
c produced by program propzon.
c This routine does a coupling calculation in frequency intervals, storing all
c hybrid singlets in a central interval of the coupling bandwidth.
c the trick is to read the coupling partners for each block, calculate
c receiver and source scalars, then read the eigenvectors one by one.

As with the spherical-earth greensfunction code, this code creates waveforms for a source-receiver geometry stored in the SAC header of a common seismic data file. Unlike the sph-earth code, the output is tied to geographic coordinates (east/north), not radial/transverse. The code syndat.f takes the greensfunctions and sums with a source mechanism to form synthetic seismograms.

sample input:

0 --> # of point to compute, 0=same number as in the data file
ffc.lhz --> data file for source-receiver-station info
gzB0.05_ffc --> filename root for greenfct files, code will append ".lhz" ".lhe" ".lhn"
../Vary/ymodel10_B0.05.zon20 --> input file of hybrid singlets, computed with zonchn.f
/Users/jjpark/Mineos/Models/premiss --> earth model (isotropic PREM in this case)
/Users/jjpark/Mineos/sphnl --> N,L file of spheroidal mode parameters
/Users/jjpark/Mineos/sphe --> direct-access Fortran file of spheroidal modal eigenfunctions
/Users/jjpark/Mineos/tornl --> N,L file of toroidal mode parameters
/Users/jjpark/Mineos/tore --> direct-access Fortran file of toroidal modal eigenfunctions
stop

zonsubs.f -- the subroutines needed to make zonchn.f work

2009-11-05T11:22:00.001-05:00

This set of subroutines serves zonchn.f, see the earlier post. The first routine yread reads the parameters of the zonally symmetric model example file here

1 1 ! useless parameters for option I never implemented
0.0 0.05 0.0 0.0 0.0 A,B,C,D,E relative perts to the Backus anisotropy parameters for azim anisotropy with an axis of symmetry.
1 152 181 ! lithosphere in model premiss
10 ! k=10 wavenumber on the sphere

Some relevant passages of code are the following:

c angint is original theta integral for cos{(l+1/2)\theta} structure
c angint2 uses 3-j tricks to speed up the calculation
c angint2X uses 3-j tricks, but for either spherical cap or equatorial band
c angint2XX is s=20 expansion of equatorial band model, no discontinuites
c call angint2(l1,l2,ls(kii),angi)
call angint2X(l1,l2,ls(kii),angi)
c call angint2XX(l1,l2,ls(kii),angi)

c combine isotropic and anisotropic parts, use b,c,e & a,d
c **************************** note that code can incorporate only
c ***************************** one aspherical wavenumber, despite
c ***************************** the do loops
c NOTE: some lines have alternate versions for phi and theta as axis of sym
if(ick.eq.3) then ! tor-tor
do j=1,nom
coup(j)=coup(j)
x +angi(3,j)*(c*ysc(1,1)+e*ysc(2,1)) ! anisotropic
c x -angi(5,j)*e*ysc(3,1)-angi(6,j)*c*ysc(1,1) ! anisotropic-phi
x +angi(5,j)*e*ysc(3,1)-angi(6,j)*c*ysc(1,1) ! anisotropic-theta
x +angi(2,j)*d*ysc(3,1)+angi(3,j)*d*ysc(2,1) ! isotropic
end do

comment or uncomment the phi and theta terms as you wish, but dont comment both or uncomment both.

zonchn.f -- code to couple free oscillations for a zonally symmetric earth model

2009-11-05T10:13:00.002-05:00

zonchn.f is a code to couple together free oscillations of Earth for a set of highly restricted aspherical earth models. The models are zonally symmetric, dependent only on co-latitude. The code has subroutines that integrate the modes against two types of models (1) a model that depends on a single wavenumber k, and (2) a model that is an equatorial band or a symmetric spherical cap. The first model type was used to examine the wavenumber behavior of surface-wave coupling, and the second type simulates the interaction of surface waves with a sharp boundary. The paper that describes the code's workings is

Park, J., The sensitivity of seismic free oscillations to upper mantle anisotropy I. Zonal symmetry, J. Geophys. Res., v98, 19933-19949, 1993.

Zonal symmetry in the model allows the coupling to vanish for all free-oscillation singlets with differing azimuthal numbers m, m'. In this manner, if there are two modes nSl and n'Sl' with a total of 2l'+1+2l+1 = 2(l+l'+1) singlets, the coupling problem splits apart into (for l'>l) 2l+1 2x2 coupling matrices, and l'-l uncoupled singlets. As a result of this symmetry, zonchn can readily accept a large number of modes in a narrow frequency band without breaking the computer memory or the CPU. The general scheme is to march up in frequency in steps of df and Df: All modes in a half-bandwidth Df about some fiducial frequency fo are coupled, but only those hybrid modes in a narrower bandwidth df about fo are saved. The code increases fo by 2*df and iterates up to some cutoff frequency. The user specifies how many overtones n are included. The uploaded version of zonchn.f has been tested with 0.2.le.f.le.20 mHz and df=0.1 mHz, Df=0.5 mHz and n-max=5.

The code depends on a set of subroutines in the file zonsubs.f, which are posted next in this blog.

The code reads from an input file "phzonchn" that contains

print *,'open(15,file=phzonchn,form=formatted)'
open(15,file='phzonchn',form='formatted')
call yread
print *,'output chain file '
read(15,102)namchn
print *,namchn
102 format(80a)
print *,'initial open of output file in append mode (0=no) '
read(15,*) iappend
print *,'iappend= ',iappend
print *,'enter inner freq halfwidth, outer freq halfwidth (mHz)'
c outer band is modes that couple, inner band contains modes that are saved
read(15,*) hfin,hfout
print *,hfin,hfout
print *,'enter freq band for hybrid modes (mHz)'
read(15,*) fb1,fb2
print *,fb1,fb2
fb1=dmax1(fb1,0.25d0)
print *,'enter max overtone branch to consider'
read(15,*) nbrmax

The option of writing the output file in "append" mode is an artifact of slow computers of the 1990s, in which one might need to stop and restart a long computation. I doubt that the user will need this option in 2009, though my test calculation required the greater part of an hour to complete.

example for phzonchn:

ymodel10_B0.05 --> aspherical model file
ymodel10_B0.05.zon20 --> output file of hybrid modes (linear comb of sph-earth modes, new frequencies)
0 --> dont open output file in append mode (usual case)
0.1 0.5 --> inner band contains modes that are saved, outer band is modes that couple (mHz)
0.2 20 --> frequency band to compute (mHz)
5 --> max number of overtones
/Users/jjpark/Mineos/Models/premiss --> earth model (isotropic PREM in this case)
/Users/jjpark/Mineos/sphnl --> N,L file of spheroidal mode parameters
/Users/jjpark/Mineos/sphe --> direct-access Fortran file of spheroidal modal eigenfunctions
/Users/jjpark/Mineos/tornl --> N,L file of toroidal mode parameters
/Users/jjpark/Mineos/tore --> direct-access Fortran file of toroidal modal eigenfunctions

The model file is based on the A,B,C,D,E Backus parameterization of hexagonally-symmetric anisotropy. A and D are isotropic P and S wavespeed perturbations (peak2peak variation: A=0.05 implies 5%). B is the 2\phi azimuthal Vp term, C is the 4\phi azimuthal Vp term, and E is the 2\phi azimuthal Vs term. The input, as explained in the post for zonsubs.f iare the values of the ABCDE and the wavenumber k of a corrugation in colatitude. The code does not allow the user to choose the equatorial band model on input. You must change a subroutine call and recompile the code. I may fix this soon, as it is cumbersome. Similarly, there is a slight alteration in the coupling formulas that switches the cases for an anisotropic symmetry axis w-hat aligned with theta-hat, or with phi-hat, that is, with colatitude and longitude, respectively. Due to symmetry restrictions, only w-hat parallel to theta-hat or phi-hat is allowed by this code. No tilted axes of symmetry, no general azimuthal angle.

syndat.f -- sum greens functions with moment tensor and centroid time correction

2009-10-30T10:11:00.002-04:00

This code is a version of the old sac_syndat.f that incorporates the IRIS-suppied EVRESP subroutine to compute and apply instrument response correction to convert displacement to instrument counts.

c f77 -o /Users/jjpark/bin/syndat syndat.f /Users/jjpark/Ritz/libevresp.a /Users/jjpark/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/syndat syndat.f /Users/jjpark/Ritz/libevresp.a /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/jlib.a
c program to make synthetics using green fcts
c the greenfcts are raw displacement in meters. This code applies the instrument
c response in the frequency domain with a FFT --> multiply by H(f) --> IFFT process.

program asks for the centroid time delay, which is applied with a phase-ramp in the frequency domain
The centroid time delay occurs because the start time of a large event is typically not the best time to place the stepfunction onset of the long-period course. If moment release occurs over T seconds, then a better fit to data is found by delaying the synthetic onset by T/2 seconds. (Assuming a symmetric source time function, natch).

print *,'enter time delay e.g. centroid-time correction'
read(5,*) tconst
40 print *,'do you want to use strike,dip and rake (0=no) '
read(5,*) ick
if(ick.ne.0) then
call fault(sol)
else
print *,'enter f(x6) : '
read(5,*) (sol(i),i=1,6)
endif
print *,'routine expects the gfct filename to start with g'
print *,'the code substitutes s for g in the output filename'
30 print *,'green fns file name(=stop) : '

The program will read a succession of greensfct files, one per station-component, until it reads 'stop' e.g.

0 --> zero centroid time delay
1 --> read strike,dip,slip
329 8 110 --> strike dip slip
400 --> seismic moment (units of 10**27 dyne-cm/10**20 nt-m)
/Users/jjpark/Synth/gsph_ffc.lhz --> greensfct file
/Users/jjpark/Synth/gsph_ffc.lhr
/Users/jjpark/Synth/gsph_ffc.lht
stop

There are PLOTIT calls within the code (commented out for now) that allow the user to view his/her handiwork, testing the causality of the instrument response, etc.

The EVRESP routine is found in the libevresp.a library that you can create from the IRIS EVALRESP software distribution. For reasons only the DMC programmers know, the latest versions of the EVALRESP distribution hide the libevresp.a library in a subdirectory ".lib" that wont show up in a "ls" command unless you are looking for it. Just another convenience feature, I guess.

Evalresp follies

2009-10-30T09:33:00.002-04:00

The IRIS-supplied instrument response software was not straightforward to incorporate into our codes. The key C-subroutine EVRESP is built for incorporating into Fortran codes, but seems to have developed a bug over its past revisions, and no one may have complained of it. The subroutine generates the filename of a RESP file, nominally output by rdseed. In previous versions of the subroutine I have used (circa 2005) the routine opens the file, reads the contents and generates a frequency domain instrument response for convolving synthetics or deconvolving data. The 2008-9 versions of EVRESP cant seem to find the file whose name they create (seemingly correctly) and the program fails. Diddling with the SEEDRESP environmental variable does not seem to help. The fix that I found was to feed EVRESP the RESP.* filename explicitly, by generating the filename myself from the same group of character-string variables that EVRESP uses. To wit:

file='RESP.'//net//'.'//sta//'.'//locid//'.'//cha
print *,file
c remove blank spaces
do i=30,7,-1
if(file(i:i).eq.' ') file=file(1:i-1)//file(i+1:80)//' '
end do
print *,file
print *,'entering evresp'
print *,sta,cha,net,locid,datime,unts,nscan,npad,nha
iflag=evresp(sta,cha,net,locid,datime,unts,file,freq,nha,resp,
x rtyp,vbs,start_stage,stop_stage,stdio_flag)
print *,'done evresp'

The character string "file" is passed to the routine and the code looks for it successfully, if the file is there!

sacread.f -- subroutines to read and write SAC-format files -- adjusted to write 6-component greens functions

2009-10-29T12:19:00.005-04:00

sacread.f

subroutines to read and write SAC-format files -- adjusted to write 6-component greens functions. The routines provided by SAC are fine, of course, but I use a hydrid file format for greens-function files. These files contain six displacement traces rather than one, one for each independent component of the seismic moment tensor. There is a flag in the SAC header that I took over for indicating that six traces need to be read instead of one.

if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
nword=nscan+158
nrec=(nword-1)/128+1
last=nword-128*(nrec-1)
else
print *,'reading greens functions'
nword=6*nscan+158
nrec=(nword-1)/128+1
last=nword-128*(nrec-1)
endif

The subroutines are sacread, sacout, sacread_swap and sacout_swap. The last two routines swap big-endian and little-endian data. You will be happier if you never need to swap bytes, but they are there for when you are painted into that corner.

csolve.f -- solve a system of equations with complex coefficients

2009-10-29T12:14:00.002-04:00

csolve.f is a subroutine for doing Gaussian elimination with partial pivoting to solve a system of linear equations A.x=y with complex matrix A and vector y. Useful in several routines.

Snow Leopard transition: Fortran

2009-10-28T14:04:00.001-04:00

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.

Code for computing greens function for synthetic seismograms using summed free oscillations for spherical earth models

2009-10-20T11:49:00.008-04:00

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

Rayleigh Wave Dispersion Code

2008-08-01T13:04:00.005-04:00

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

2008-08-01T12:34:00.003-04:00

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

Viewing the output of the sacw_corr_many code

2008-07-01T15:31:00.003-04:00

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

Using the wavelet-based surface wave code

2008-07-01T14:34:00.006-04:00

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

2008-07-01T14:27:00.002-04:00

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? <-- routine subs e,n,z for "?"
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

changes in the inverse FFT in recfunk codes

2008-05-27T22:09:00.003-04:00

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.

These codes have been updated on the website:

Update of anirec_synth.f and anirec_synth_circle.f

2008-05-27T00:38:00.004-04:00

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

The erstwhile version of aniprop.f has had a bugfix

2008-03-24T18:46:00.003-04:00

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.

http://earth.geology.yale.edu/~jjpark/Code/aniprop_072307.f

after a bugfix is now

http://earth.geology.yale.edu/~jjpark/Code/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.

plotting script for rfmig_cboot.f and rfmig_mcboot.f output WITH UNCERTAINTIES

2007-12-06T12:49:00.001-05:00

# ./Plotj_cexp outr_cexp.grid outr1_cexp.grid outr2_cexp.grid outt_cexp.grid outt1_cexp.grid outt2_cexp.grid 0.15 KEG

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_cboot.f or rfmig_mcboot.f
# Input - $1 - mean R+T $2 lower R+T $3 upper R+T $4 mean R-T $5 lower R-T $6 upper R-T $7 scale $8 TITLE

plotting script for rfmig_boot.f and rfmig_mboot.f output WITH UNCERTAINTIES

2007-12-06T12:47:00.000-05:00

# ../Plotj_bexp outr_bexp.grid outr1_bexp.grid outr2_bexp.grid outt_bexp.grid outt1_bexp.grid outt2_bexp.grid 0.3 title

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_cboot.f or rfmig_mcboot.f
# Input - $1 - mean R $2 lower R $3 upper R $4 mean T $5 lower T $6 upper T $7 scale $8 TITLE

plotting script for rfmig_cboot.f and rfmig_mcboot.f output

2007-12-06T12:46:00.000-05:00

#../Plot_cexp outr_cexp.grid outt_cexp.grid 0.15 KEG

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_cboot.f or rfmig_mcboot.f
#Input - $1 - file name $2 file name $3 - scale $4 TITLE

plotting script for rfmig_boot.f and rfmig_mboot.f output

2007-12-06T12:43:00.000-05:00

#../Plot_bexp outr_bexp.grid outt_bexp.grid 0.15 KEG

#shell to plot RF wiggles from the "grid" files generated by
# rfmig_boot.f or rfmig_mboot.f
#Input - $1 - file name $2 file name $3 - scale $4 TITLE

rfmig_mcboot.f -- code to compute complex-valued RF harmonic expansion in back-azimuth, moving-window moveout correction

2007-12-04T23:57:00.000-05:00

c computes moving-window moveout correction for MTC receiver functions
c applied in the frequency domain.
c requires a stacking model in the anirec format
c such a model may have anisotropy parameters,
c but migration code only uses the isotropic velocities.
c
c code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. The constant term should be zero for the transverse RF.
c The 2-lobed terms govern dipping interface effects and tilted symmetry-axis ansotropy.
c The 4-lobed term is anisotropy with a horizontal axis
c The code regresses for the harmonic expansion, using combined radial/transverse stack
c bootstrap-resamples the data to estimate the
c uncertainty of the harmonic terms.
c
c output files are out[rt]_cexp.grid -- harmonic-expansions of the RFs, in time domain
c out[rt]1_cexp.grid -- harmonic-expansion RFs plus bootstrap uncertainty
c out[rt]2_cexp.grid -- harmonic-expansion RFs minus bootstrap uncertainty
c out[rt]_bbaz.grid -- harmonic-expansion RFs computed for ordered baz values
c
c has kluge to cheat the pre-event noise for synthetic records 3/12/00 JJP
c check to see if the kluge is commented out
c
c this version of the RF code reads a file of data filenames
c you have two choices: either read the time intervals in the filename file
c or read them in the sac header
c the data must be binary SAC format
c horizontals must be rotated to radial and transverse

nboot=0: compute simple regression, no bootstrap computation of variance.

rfmig_mboot.f -- code to compute RF harmonic expansion in back-azimuth, moving-window moveout correction

2007-12-04T23:54:00.000-05:00

c program rfmig_mboot
c 10/12/04 JJP -- adapted from rfmigrate
c
c xf77 -o /park/backup/bin/rfmig_mboot rfmig_mboot.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/rfmig_mboot rfmig_mboot.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c computes moving-window moveout correction for MTC receiver functions
c applied in the frequency domain.
c requires a stacking model in the anirec format
c such a model may have anisotropy parameters,
c but migration code only uses the isotropic velocities.
c
c nboot=0 -- only compute a single regression for RF harmonic expansion --> no bootstrap uncertainty estimate
c
c code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. The constant term should be zero for the transverse RF.
c The 2-lobed terms govern dipping interface effects and tilted symmetry-axis ansotropy.
c The 4-lobed term is anisotropy with a horizontal axis
c The code regresses for the harmonic expansion, and bootstrap-resamples the data to estimate the
c uncertainty of the harmonic terms.
c
c output files are out[rt]_bexp.grid -- harmonic-expansions of the RFs, in time domain
c out[rt]1_bexp.grid -- harmonic-expansion RFs plus bootstrap uncertainty
c out[rt]2_bexp.grid -- harmonic-expansion RFs minus bootstrap uncertainty
c out[rt]_bbaz.grid -- harmonic-expansion RFs computed for ordered baz values
c
c has kluge to cheat the pre-event noise for synthetic records 3/12/00 JJP
c check to see if the kluge is commented out

rfmig_cboot.f -- code to compute harmonic expansion of RFs with back-azimuth

2007-12-02T22:31:00.000-05:00

rfmig_boot.f -- code to estimate harmonic expansion of RFs in back azimuth

2007-12-01T21:14:00.000-05:00

c code computes frequency-domain stacks of receiver functions that follow a harmonic expansion in baz
c for both radial and transverse RFs there are constant terms and sin/cos terms for 2- and 4-lobed
c amplitude dependence. The constant term should be zero for the transverse RF.
c The 2-lobed terms govern dipping interface effects and tilted symmetry-axis ansotropy.
c The 4-lobed term is anisotropy with a horizontal axis
c The code regresses for the harmonic expansion, and bootstrap-resamples the data to estimate the
c uncertainty of the harmonic terms. The posted version only
c
c output files are out[rt]_bexp.grid -- harmonic-expansions of the RFs, in time domain
c out[rt]1_bexp.grid -- harmonic-expansion RFs plus bootstrap uncertainty
c out[rt]2_bexp.grid -- harmonic-expansion RFs minus bootstrap uncertainty
c out[rt]_bbaz.grid -- harmonic-expansion RFs computed for ordered baz values
c
c migrates MTC receiver functions in the frequency domain.
c requires a stacking model in the anirec format (see previous posts on codes that apply
c a frequency-domain moveout correction for the variation of Ps delay time with epicentral distance.)
c such a model may have anisotropy parameters,
c but migration code only uses the isotropic velocities.
c
c has kluge to cheat the pre-event noise for synthetic records 3/12/00 JJP
c check to see if the kluge is commented out
c
c this version of the RF code reads a file of data filenames
c you have two choices: either read the time intervals in the filename file
c or read them in the sac header
c the data must be binary SAC format
c horizontals must be rotated to radial and transverse
c
c for start times in the file:
c the file is "in_recfunk" and has lines of the form:
c
c 1997.070.19.33.bh? <-- code replaces ? with r,t,z
c 57 52 <-- start time of analysis window, window duration (sec)
c 1997.076.08.15.bh?
c 62 62
c ...
c ...
c ...
c stop <-- tells code that data is finished, 799 events max
c
c
c for start times in the SAC header
c reads seismic record start times from the sac header
c will search the SAC header for specific markers of P phases
c T1 - P, Pdiff ahead(12)
c T2 - PKP,PKIKP ahead(13)
c T3 - PP ahead(14)
c T1=T2=T3=0 ==> use original A-marker ahead(9)
c
c code does NOT combine data with different sample rates
c data files limited to 99K pnts. To increase, see common block /datastuff/
c
c many intermediate quantities are plotted with PLOTIT as the code proceeds.
c other intermediate quantities can be plotted by uncommenting calls to PLOTIT

anirec_synth_circle.f -- make P-coda synthetics at evenly-spaced back azimuth and constant epicentral distance

2007-12-01T21:11:00.000-05:00

anirec_synth.f --- a code to compute synthetic P coda from a collection of SAC-format data

2007-12-01T21:03:00.000-05:00

c anirec_synth - anirec to compute synthetics for a specified model and
c a specified collection of seismic records (read from in_recpick-format file of filenames)
c modified to accept SACfiles with start time info in the header, in A-marker (ahead(9))
c because the timing information is not used, the value of A-marker is not important.
c However, in the synthetics, the A-marker and the T1-marker (or T2) variables in the header are set to 5.0 seconds
c for easy analysis by recfunk codes that use the header information for timing.
c code divides between T1 (P) or T2 (PKP) markers based on epicentral distance, with 120 degrees the divide
c the synthetic files are s_filename.bh[zrt], where filename is the name of the datafile.
c the code writes in_recpick_synth, a list of the filenames for recfunk_pick and other recfunk codes that
c expect the start time in the SAC header
c

The goal of this code is to generate P-coda form a specified model to replicate the earthquake distribution of a real data set. This allows the analyst to test whether the oddities in an RF sweep are likely caused by an imperfect data distribution, or by a shortcoming of his/her model.

recfunk_mwms.f -- moving-window migration that splices together a long RF from a series of target depths

2007-08-22T11:28:00.000-04:00

This code can be accessed by clicking the title link

recfunk_mwms.f generates a succession of RF estimates for a succession of migration target depths. It splices together the results into a long RF. I have used this code to reconstruct BAZ-sweeps of RFs from 0 to 75-s delay, encompassing the crust to the 660-km discontinuity.

This code will warm your laptop. For a splice from 0 to 680 depth, targeted every 40 km with 524 events, my Macbook Pro required 15 minutes of run time.

INPUT
'velocity model for stacking (blank = stack_model)'
'depth range and inc for the RFs (km) (dep1,dep2,ddep)'
'enter fmax (Hz)'
'time intervals in file (0) or in SAC-header (1)'
'file of seismic records to read? (space-ret: in_recpick)'
'duration of data windows for analysis'
'minimum number of events in a binsum (eg 1 or 2)'
'rotate to LQT coordinates to isolate P and SV? (0=no)'

--->> then the first pass through the seismic data -- typically for target depth = 0 km.

'change baz-interval or baz-spacing? (1=yes)'
'enter back-azimuth range: bazmax, bazmin, bazinc'

--->> compute the BAZ-interval binsums

'compute RFs binned with epicentral distance: enter back-azimuth limits ib1,ib2 (integers!)'
' do you want to change 'epicentral spacing? (1=yes)'
'enter delta range and spacing (ep1,ep2,dep)'

Then Computer code takes off, remembering the binsum parameters for successive target depths.

OUTPUT

grid files for short interval of time around each target depth NNN-km

out[rt]_baz_NNN.grid, out[rt]_epi_NNN.grid

grid file for spliced baz- and epicen-sweeps

ous[rt]_baz.grid, ous[rt]_epi.grid

You can plot these files with the posted GMT scripts, after adjusting the BOX variable to change plot limits.

recfunk_mwmj.f -- moving-window moveout RF estimation with jackknife uncertainty on the bin-sums

2007-08-20T15:54:00.000-04:00

c program recfunk_mwmj
c 6/3/04 JJP --- UPDATED 08/20/07 JJP
c MIGRATION/MOVEOUT variant to compute cross-correlation in a moving time window
c This version performs a jackknife uncertainty estimate of bin-summed RFs
c
c you tell the program the depth at which you want to focus the cross-correlation
c and it computes the moving-window delay from an input model
c using the epicentral distance to look up the horizontal slowness

c the code writes the BAZ- and EPICEN-dependent RFs to files
c in a format easily digested by GMT (traces are separated by ">" lines)
c
c filenames: out[rt]_baz_NNN.grid out[rt]_epi_NNN.grid
c where NNN corresponds to the target depth in kilometers
c jackknife uncertainty envelope (lower) out[rt]1_baz_NNN.grid out[rt]1_epi_NNN.grid
c jackknife uncertainty envelope (upper) out[rt]2_baz_NNN.grid out[rt]2_epi_NNN.grid
c
c these files are overwritten everytime you run the program
c so rename them to save them
c
c use the GMT scripts in the previous few blog posts for plotting these files.

GMT shellscript syntax that is not showing up in the blog posts CAVEAT EMPTOR

2007-08-20T15:38:00.000-04:00

GMT script for the moving-window with jackknife uncertainty -- epicentral bin-sums

2007-08-20T14:59:00.001-04:00

NOTE -- there is a bit of shellscript syntax that is not showing up in the blog post. There is a shell command

process leftangle-leftangle END rightangle-rightangle plot.ps

that instructs "process" to read the succeeding lines of data until it reaches a line that contains just "END"
In my blog posts the phrase "leftangle-leftangle END rightangle-rightangle" is printing as "leftangle-rightangle"

So you must correct these lines for the shellscripts to work

#!/bin/csh -f

#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue

echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4

/bin/rm plot.ps

set scale = $7

set title = $8

pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps

pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 180
END

pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps

pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 180
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 180
END

echo gv is JPark's alias for ghostscript
gv plot.ps

GMT script for the moving-window with jackknife uncertainty -- back-azimuth bin-sums

2007-08-20T14:50:00.000-04:00

#!/bin/csh -f

#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1 $2 $3 - radial RF filenames $4 $5 $6 transverse RF filenames $7 - scale $8 title
# filenames are for mean RF estimate, and lower and upper jackknife uncertainty bounds on the RF estimate
# the script uses PS wiggle to draw a green envelope that represents the uncertainty range
# Then plots the mean RF estimate within the green uncertainty
# two calls to pswiggle for each jackknife file allow the "certain" RF excursion to be red or blue

echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-12/12/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4

/bin/rm plot.ps

set scale = $7

set title = $8
gmtset ANOT_FONT Times-Roman HEADER_FONT Times-Roman LABEL_FONT Times-Roman

pswiggle $2 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SWen -P -K > plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps
pswiggle $3 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SWen -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\Swen -P -K -O >> plot.ps

pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $7 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 380
END

pswiggle $5 $FRAME $BOX -Z$scale -M -G200/255/200 -Ba4f1/f30\SwEn -P -O -K -X3.5 >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G200/255/200 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps
pswiggle $6 $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/f30\SwEn -P -O -K >> plot.ps
pswiggle $5 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/f30\SwEn -P -K -O >> plot.ps

pswiggle $4 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $4 scale $7 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $8 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 380
END

echo gv is JPark's alias for ghostview
gv plot.ps $

stacking model example

2007-08-20T13:41:00.000-04:00

stacking model example (used for GSN Station RAYN) 4 layers over a halfspace
4
50 145
5000 5500 -0.06 0 3070 -0.06 2700
0 0
8000 6350 0 0 3494 0 2700
0 0
21000 6420 0.02 0 3510 0.02 2700
0 0
40000 6930 0.03 0 3804 0.03 2700
0 0
220000 7900 -0.15 0 4300 -0.15 3350

GMT script for plotting epicentral RF sweep from moving-window-moveout grid files

2007-08-20T13:26:00.000-04:00

#!/bin/csh -f

#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by EPI
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
# The first bit of business is to re-order the traces

/bin/rm wigtab? tmp? plot.ps

echo standard usage
echo "./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-10/10/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4

awk '{print $2}' $1 | sort -u -n > bazlist
foreach i ( `cat bazlist` )
echo working on $i
awk '$2 == A && $1 >= -10' A=$i $1 >> wigtabR
awk '$2 == A && $1 >= -10' A=$i $2 >> wigtabT

#in case you decide to demean...
#awk '$2 == A && $1 >= -2' A=$i $1 > tmp1
#set mean = `awk '{a+=$3} END {print a/NR}' tmp1 `
#echo baz $i mean $mean
#awk '{print $1,$2,$3-mn}' mn=$mean tmp1 >> wigtab

echo ">" >> wigtabR
echo ">" >> wigtabT

end

set scale = $3
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SWen -P -K > plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SWen -P -K -O >> plot.ps

pswiggle wigtabR $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END
# psxy dipping_6.5_3.6_43_15_20.dat -: $FRAME $BOX -W1.5p/255/0/0 -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
5 -10
5 180
END

pswiggle wigtabT $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SwEn -P -O -K -X3.5 >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SwEn -P -K -O >> plot.ps

pswiggle wigtabT $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $2 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END

psxy $FRAME $BOX -W1.5p/255/0/0 -O <> plot.ps
5 -10
5 180
END

gv plot.ps &
echo gv is aliased by JPark to ghostview

GMT script for plotting backazimuthal RF sweep from moving-window-moveout grid files

2007-08-20T13:22:00.000-04:00

#!/bin/csh -f

#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_baz outr_baz.grid outt_baz.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles

echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
# set BOX = -R25/65/-5/365
set BOX = -R-10/10/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4

/bin/rm plot.ps

set scale = $3

set title = $4

pswiggle $1 $FRAME $BOX -Z$scale -M -G0/0/255 -P -K > plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SWen -P -K -O >> plot.ps

pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 380
END

pswiggle $2 $FRAME $BOX -Z$scale -M -G0/0/255 -P -O -K -X3.5 >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SwEn -P -K -O >> plot.ps

pswiggle $2 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 CT $2 scale $3 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -1.25 2.7 24 0 5 CT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 380
END

gv plot.ps &
echo gv is aliased by JPark to ghostview

recfunk_mwm.f -- moving-window moveout

2007-08-20T13:03:00.000-04:00

click the title of the post to access the code.

This code applies a moveout correction by shifting the time window of the horizontal-component time series before DFT and multiple-taper spectral correlation. The code prompts the user to specify a stacking model in the anirec format, and a target depth. The code uses the phase type (P,PP or PKP) and epicentral distance to compute a predicted time that a Ps phase from the target depth should arrive. The horizontal time series is shifted by this time, so that the RF at t=0 is the Ps pulse at the target depth. The zero-time pulse on the radial RF associated with the direct P wave is shifted to negative times, and will often lose focus in back-azimuthal stacks. The GMT script for plotting is similar to that used for other recfunk codes, but has a shifted X-axis to focus the eye on the signals at Ps delay-times both before and after the target Ps arrival.

The code writes output files that flag the target depth:

outr_baz_040.grid,outt_baz_040.grid,outr_epi_040.grid,outt_epi_040.grid

where the integer 040 stands for a target depth of 40 km in the stacking model. GMT scripts for plotting are executed this way:

./Plot_baz outr_baz_040.grid outt_baz_040.grid 1. RAYN_mwm

./Plot_epi outr_epi_040.grid outt_epi_040.grid 1. RAYN_mwm

rfmigrate.f -- recfunk version with migration

2007-08-19T21:11:00.000-04:00

Click the title of the post to access the code.

Version of MTC RF estimation that computes single-station RFs with 3 2.5-pi Slepian eigentapers, similar to recfunk_pick.f. This code uses a migration scheme to align Ps pulses to the time expected for Ps from a vertical-incidence wave. See comments for the rfmigrate_estack code, posted yesterday (August 18, 2007).

recfunk_svd.f version of recfunk with SVD in the inner loop

2007-08-19T09:07:00.000-04:00

click the title of post to obtain the code

From unpublished manuscript on RF estimation:

Receiver Functions from Multiple-Taper Spectral Correlation:
Statistics, Single-Station Migration, and Jackknife Uncertainty
J. Park and V. Levin, submitted to GJI

RF-computation suggests that the coherence between P and SH-converted energy in the P coda is comparable to the coherence between P and SV-converted energy. Because aggregate coherence C(f)< 0.5, one can question whether the SV and SH converted phases are mutually correlated. The alternative hypothesis is that SV and SH scattering arise from statistically distinct portions of the incoming P wave. Poor correlation between SV and SH scattering can be tested by simultaneous cross-correlation of all three particle-motion components, computed via the singular-value decomposition (SVD). Define a K ×3 matrix A(f) of eigenspectral estimates in the P-SV -SH coordinate system (in our application K = 3):

For noise-free, perfectly correlated data, the matrix A(f) = u (f) v(f), a single dyad in which the 3-component vector v(f) contains the correlation coefficients between components and the K-component vector u(f) describes the spectral signature of the correlated signal in a narrow frequency band about f. The receiver functions HSV (f),HSH(f) correspond to the component ratios vSV /vP and vSH/vP , where v(f) = [vP , vSV , vSH]. For real noisy data, we associate the SVD receiver function with the singular vectors of A that correspond to the largest singular value, similar to the ”pure state” filtering algorithm (Samson, 1983). We SVD decompose A(f). The uncertainty of the SVD receiver function can be estimated by analogy to uncertainties of RF computed from complex cross-correlation and coherence.

In tests with real data, the added sophistication of the SVD receiver function algorithm leads only to small changes in the estimated RFs. We infer that P-SV conversion and P-SH conversion in the P coda are co-dependent, in a statistical sense. It should therefore be feasible to estimate HSV (f) and HSH(f) jointly. On the other hand, our experiments suggest that the ordinary pairwise MTC RF estimator is adequate in practice. The SVD-based RF estimator may have utility in complex situations, such as the estimation of S receiver functions for S-to-P converted energy.

rfmigrate_estack.f -- code to migrate the RFs in the frequency domain

2007-08-18T10:14:00.000-04:00

Click the link to access the code.

Long-time window event-stacking of RFs. This code takes three swipes at the data. First it does a normal event-stack without migration in back-azimuth bins (output files out[rt]_baz.grid). Then a normal event-stack without migration in epicentral bins (output files out[rt]_epi.grid). Lastly it does an event-stack with migration in back-azimuth bins (output files oum[rt]_baz.grid).
For GMT scripts to plot these output files, see blog posts from July.

ttimes_PKP.dat -- accessory file for RF migration and moveout corrections

2007-08-18T10:09:00.001-04:00

This file is used by RF migration and moveout codes (soon to be posted). It uses IASPEI traveltimes to scale the expected Ps traveltime delay as a function of epicentral distance for PKP/PKIKP. Numbers massaged from ttimes, the public code available from IRIS. Click the title to access the file.

ttimes_P.dat -- accessory file for receiver function moveout correction

2007-08-18T09:56:00.000-04:00

This file is used by RF migration and moveout codes (soon to be posted). It uses IASPEI traveltimes to scale the expected Ps traveltime delay as a function of epicentral distance for P and Pdiff. Numbers massaged from ttimes, the public code available from IRIS.

recfunk_estack.f -- an event-stacking algorithm for RF with long delay times

2007-08-17T14:19:00.000-04:00

Click the title to access the Fortran code.

recfunk_estack is a version of the multiple taper correlation receiver-function estimator that allows the time-domain RF to extend to times longer than (roughly) T/10, where T is the duration of the analysis window. Recfunk_estack does NOT apply a moveout or migration correction, so the ability to resolve Ps from deeper interfaces is limited. Other codes to be posted soon will have these properties.

The major shortcoming of the MTC RF estimator of Park & Levin (2000) is its difficulty retrieving Ps conversions at long delay times. Receiver function studies of upper mantle structure at the base of the lithosphere and the transition zone involve Ps converted phases with delay times between 20 and 70 seconds (e.g., Bostock, 1998; Sheehan et al 2000), but the frequency averaging properties of the Slepian data tapers with time-bandwidth product p\ge 2 tends to limit typical RFs to delay times \tau<20 s.

Frequency averaging affords the MTC algorithm desirable statistical properties, especially for single-event RF estimates. We can achieve a narrower frequency average for the retrieval of Ps from interfaces deeper than 100 km if we sum single spectral estimates over K seismic events, rather than K eigenspectra from one event, in the RF cross-correlation formulas (2) and(3). Using only one spectrum estimate per component for each P wave, we estimate the spectra with a Slepian taper of time-bandwidth p=1, that is, the "one-\pi-prolate" taper. This taper averages local frequency information over a central bandwidth similar to the simple boxcar taper, but is optimized for spectral leakage resistance.

Modified recfunk_pick_through.f

2007-08-16T23:07:00.000-04:00

More options for specifying the start time of the data record are recognised.

c program recfunk_pick_through
c version of recfunk_pick that reads seismic record start times from
c the sac header, and writes to in_recpick the filenames that have been picked
c using EITHER the A-marker or the T1, T2, or T3 markers in the SAC header, placed by running plotpk
c within SAC. Set by pressing "T" and 1,2 or 3 in plotpk sequentially
c 10/13/05 JJP
c input file: in_recpick_raw
c output file: in_recpick
c
c xf77 -o /park/backup/bin/recfunk_pick_through recfunk_pick_through.f /park/backup/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/recfunk_pick_through recfunk_pick_through.f /Users/jjpark/Ritz/jlib.a
c
......implicit real*4(a-h,o-z)
......real*8 ar,ai,el(12),aaa
......complex*8 zc,zero
......complex*8 afft,rf,rfs
......character*80 name,subtext,string
......character*28 title(2)
......character*10 cmps(3)
......character comp(3)
......common/nnaamm/iwflag
......common/npiprol/anpi
......common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
......common/header/ahead(158)
......dimension iah(158),chead(158),fr(3),tt(3),ttt(3)
......equivalence (ahead,iah),(ahead,chead)
......data comp/'r','t','z'/
......print *,'This code reads start of interval from sac header'
......print *,'it writes filenames of those files that have been picked'
......open(10,file='in_recpick_raw',form='formatted')
......open(11,file='in_recpick',form='formatted')
10 print *,'enter filename e.g. 1998.361.bhz'
......read(10,101) name
101 format(a)
......if(name(1:4).eq.'stop') go to 111
......subtext=name
c we assume that the last character is either 'r','t',or 'z'
......do i=80,1,-1
......if(name(i:i).ne.' ') then
......namlen=i
......go to 98
......endif
......end do
......print *,name
......print *,'no filename?'
......stop
98 continue
......do kcmp=3,3
......name(namlen:namlen)=comp(kcmp)
......print *,name
......call sacread(name,a(1,kcmp),ierr)
......if(ierr.eq.1) stop
......dt=ahead(1)
......tt(1)=0.
......tt(2)=dt
......tt(3)=dt
......nscan=iah(80)
......if(nscan.gt.99000) then
......print *,'Careful! data length greater than 99000:', nscan
......pause
......endif
......end do
c end kluge
c if(baz.lt.70) then
c call plotit_axes(0.,0.,0.,0.)
c call manyplot(99000,1,nscan,3,tt,a,1.1,name,'time(sec)',' ')
c endif
c plot kluge
c KLUGE: CHANGES FOR READING TIME INTERVAL FROM SAC HEADER
c ahead(9) is SAC's A parameter, for arrival time. Set by pressing A in plotpk
c ahead(12,13,14) are SAC's T1,T2,T3 parameters, for arrival time. Set by pressing "T" and "N" in plotpk
......if(ahead(9).gt.1.0.or.ahead(12).gt.1.0.or.ahead(13).gt.1.0.
......x or.ahead(14).gt.1.0) write(11,101) subtext
......pstart=ahead(9)-ahead(6)-3.
......print *,'sach params b and a, pstart',ahead(6),ahead(9),pstart,
......x (ahead(n),n=12,14)
......go to 10
111 write(11,101) name
......close(10)
......close(11)
......stop
......end

test case

2007-08-16T22:29:00.001-04:00

test case

recfunk_pick.f -- update to allow the identification of P phase in the SAC header

2007-08-16T22:03:00.000-04:00

Update suggested by Vadim Levin

http://earth.geology.yale.edu/~jjpark/Code/recfunk_pick.f

c program recfunk_pick -- UPDATED
c reads seismic record start times from the sac header
c will search the SAC header for specific markers of P phases
c T1 - P, Pdiff ahead(12)
c T2 - PKP,PKIKP ahead(13)
c T3 - PP ahead(14)
c T1=T2=T3=0 ==> use original A-marker ahead(9)
c
c 08/16/07 JJP

This code also asks the user whether he/she wishes to apply the LQT coordinates (not everyone does!)

When you prepare P-wave files using the previously blogged procedures, or with your own algorithm, you can select "T" & "1" at the start of a P wave or Pdiff, "T" & "2" at the start of a PKP/PKIKP window, and "T" & "3" at the start of a PP window. The code change for recfunk_pick creates an epicentral distance variable to determine the LQ-rotation angle and (in other codes) the moveout/migration parameters. GCARC --> GCARC/2 for the PP phase, because PP has two equal P-phase legs.

Vadim Levin links to RF codes

2007-08-03T14:27:00.000-04:00

From my good friend and colleague Vadim Levin:

hey, cool blog :)
why don't you put a plug for my version of the RF software there,
it's at

http://www.ldeo.columbia.edu/~vadim/RF/RF-manual.html

aniprop_downrange.f -- version of aniprop to emulate a linear array of seismometers

2007-07-23T14:52:00.000-04:00

c aniprop_downrange - program to calculate propagating modes of anisotropic layer
c writes a file of dispersion curves at evenly-spaced freq points in files out_cvel out_gvel
c also writes a set of SAC_Format data files SYN_nn.BH[ZNE] for downrange distances of
c 75km 150 225 300 ... 1200 1275 1350 1424 1500km downrange
c E-component is the x component (radial) and N-component is the y component (transverse)

http://earth.geology.yale.edu/~jjpark/Code/aniprop_downrange.f

aniprop.f -- a standalone program to generate surface waves in layered anisotropic structures

2007-07-23T14:05:00.001-04:00

Sportsfans,

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.

http://earth.geology.yale.edu/~jjpark/Code/aniprop_072307.f

(after a 24 March 2008 bugfix is now http://earth.geology.yale.edu/~jjpark/Code/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.

c aniprop - program to calculate propagating modes of anisotropic layer
c writes a file of dispersion curves at evenly-space freq points
c
c reads a layered model like the following file (ignore the "c "s)
c
c K&H SoCal model, deep crust horizontal anisotropy TITLE
c 3 # OF LAYERS OVER HSPACE
c 45 45 THETA,PHI ORIENTATION ANGLES
c 4000 5500 0.06 0.00 3175 0.03 2600 FOR SYMMETRY AXIS
c 0 0
c 27400 6300 0.00 0.00 3637 0.00 2800 DEPTH (M), VP (M/S), "B", "C",
c 90 45 VS (M/S), "E"
c 32400 6800 0.04 0.00 3925 0.02 2900 B,C,E ARE ANISOTROPIC PARAMETERS
c 0 0
c 60000 7800 0.00 0.00 4500 0.00 3200 NOTE: HSPACE MUST BE ISOTROPIC
c

One quirk of this code is misbehavior for isotropic models. The code solves for Rayleigh and Love dispersion simultaneously, because the two polarizations are hybridized in anisotropic models. Paradox: for a perfectly isotropic model, the rootfinder sometimes misses roots, typically at occasions where the Rayleigh and Love dispersion curves cross. The bookkeeping for solutions is not very rigorous, so loss of a mode can happen. With anisotropy, the hybridization of Rayleigh/Love motion seems to repel the roots and keep some minimal separation. As far as I can tell, the total effect on actual synthetic seismograms is negligible.

Wave propagation in a truly 1-D lossless model leads to some interesting reverbarations that one rarely sees in nature. With all the overtones included, you can see a succession of bouncing S waves that precede the main surface wave pulse. In Earth's crust, such phases would scatter away their high frequency energy from surface topography and sedimentary layers. There is a loop for mode-summing that can be altered to eliminate the overtones, leaving the fundamental Love and Rayleigh waves: (lines 681-3)

c loop over overtones at a certain frequency
do iov=1,maxbr
c do iov=1,2

Notes on the angle conventions for w-hat, the axis of symmetry:

In the anisotropic reflectivity code, subroutine zzget *assumes* a coordinate system in which z is down (anti-vertical), x is the radial direction, and y is anti-transverse. Therefore, the position angles theta,phi for w-hat are tilt relative to down, and azimuth defined as a rotation from x towards y. This rotation is CCW if viewed from below, and CW if viewed from above. Since w-hat and -(w-hat) define the same axis of symmetry, the position angles *also* can be defined as theta=(tilt from vertical) and phi=(rotation from anti-x (anti-radial) towards anti-y (transverse)). Viewed from above, this phi rotation is CW, and defines the strike of w-hat relative to the arrival azimuth of the wave.

In order to compute seismograms for a variety of back-azimuths, the synthetic code genrecd_az.f accepts a layered model with w-hat position angles defined as theta=(tilt from vertical) phi=(strike CW from N). For an event at back-azimuth psi (CW from N), the code rotates w-hat from geographic coordinates to ray-based coordinates before passing it to subroutine zzget. If a wave arrives at back-azimuth psi, the strike of the axis of symmetry w-hat relative to its arrival azimuth is phi'=phi-psi. The code performs this rotation with this code in subroutine matget, for w-hat azimuth "az":

caz=dcosd(az)

saz=-dsind(az) ! sin(-az)

do n=1,nlp

ww(3)=w(3,n)

ww(1)=w(1,n)*caz-w(2,n)*saz

ww(2)=w(1,n)*saz+w(2,n)*caz

...

In this manner, the axes of symmetry of the model, saved in array w(.,.), are never modified.

Journal Reference:

Park, J., Surface waves in layered anisotropic structures, Geophys. J. Int. v126, 173-184, 1996.

GMT script for plotting epicentral RF sweep from recfunk *.grid files

2007-07-21T16:40:00.000-04:00

GMT script for plotting epicentral RF sweep from recfunk *.grid files

Anyone planning to use this post should know what GMT commands are. If not, google "Generic Mapping Tools"
This shell script is based on an Ur-script written by Vadim Levin

#!/bin/csh -f

#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by EPI
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles
# The first bit of business is to re-order the traces

/bin/rm wigtab? tmp? plot.ps

echo standard usage
echo "./Plot_epi outr_epi.grid outt_epi.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-2/15/0/180
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4

awk '{print $2}' $1 | sort -u -n > bazlist
foreach i ( `cat bazlist` )

echo working on $i
awk '$2 == A && $1 >= -2' A=$i $1 >> wigtabR
awk '$2 == A && $1 >= -2' A=$i $2 >> wigtabT
echo ">" >> wigtabR
echo ">" >> wigtabT

end

set scale = $3
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SWen -P -K > plot.ps
pswiggle wigtabR $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SWen -P -K -O >> plot.ps

pswiggle wigtabR $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END
# psxy dipping_6.5_3.6_43_15_20.dat -: $FRAME $BOX -W1.5p/255/0/0 -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
5 -10
5 180
END

pswiggle wigtabT $FRAME $BOX -Z$scale -M -G0/0/255 -Ba4f1/a30/SwEn -P -O -K -X3.5 >> plot.ps
pswiggle wigtabT $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/a30f15/SwEn -P -K -O >> plot.ps

pswiggle wigtabT $FRAME $BOX -Z$scale -W1p/0 -M -P -K -O >> plot.ps
echo 0 2.4 14 0 5 6 $2 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
echo -2.5 2.7 24 0 5 LT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/255/0/0 -O -K <> plot.ps
0 -10
0 180
END

psxy $FRAME $BOX -W1.5p/255/0/0 -O <> plot.ps
5 -10
5 180
END

gv plot.ps &
echo gv is aliased by JPark to ghostview

GMT script for plotting backazimuthal RF sweep from recfunk *.grid files

2007-07-21T16:27:00.000-04:00

Anyone planning to use this post should know what GMT commands are. If not, google "Generic Mapping Tools"
This shell script is based on an Ur-script written by Vadim Levin

#!/bin/csh -f
#shell to plot RF wiggles from the "grid" files generated by
#spectral coherence code of JP, arranged by BAZ
#Input - $1: radialRF filename; $2: tranvRF filename; $3: scale; $4: TITLE
# e.g. ./Plotbaz outr_baz.grid outt_baz.grid 0.5 RAYN_2.0Hz
# the scale parameter has inverse dependence -- larger scale = smaller wiggles

echo sample usage
echo "./Plotbaz outr_baz.grid outt_baz.grid 0.5 RAYN_2.0Hz"
echo "EDIT PLOTTING PARAMETERS TO SUIT YOUR DATA"
set BOX = -R-2/15/-5/365
set FRAME = -JX3/7.5
set FRAME_BIG = -JX3/9.4

set scale = $3
set title = $4

pswiggle $1 $FRAME $BOX -Z$scale -M -G0/0/255 -P -K >! plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SWen -P -K -O >> plot.ps
pswiggle $1 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps

echo 0 2.4 14 0 5 CT $1 scale $3 | pstext $FRAME -R-1/1/-2/2 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
7.5 -10
7.5 380
END

pswiggle $2 $FRAME $BOX -Z$scale -M -G0/0/255 -P -O -K -X3.5 >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -M -G255/0/0 -N -Ba4f1/30/SwEn -P -K -O >> plot.ps
pswiggle $2 $FRAME $BOX -Z$scale -W1/0 -M -P -K -O >> plot.ps

echo 0 2.4 14 0 5 CT $2 scale $3 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
echo -1.25 2.7 24 0 5 CT $4 | pstext $FRAME_BIG -R-1/1/-2/3 -N -O -K >> plot.ps
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
0 -10
0 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
2.5 -10
2.5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O -K <> plot.ps
5 -10
5 380
END
psxy $FRAME $BOX -W1.5p/0/150/0 -O <> plot.ps
7.5 -10
7.5 380
END

gv plot.ps &
echo gv is aliased by JPark to ghostview

recfunk_pick.f -- version of multiple-taper RF estimation

2007-07-21T15:29:00.000-04:00

This code is a simple variant of the multiple-taper RF estimator recfunk.f, which has been public for years on my Yale website. This version reads the start time of the P-wave analysis window from the "AMARKER" entry of the SAC header. It also rotates the vertical and radial seismogram components into LQT coordinates using a formula based on epicentral distance. This formula for LQ rotation angle is very crude, but simple measures of the actual rotation angle scatter quite a bit, and every station has unique behavior. This variability suggests that careful derivation of a precise LQ rotation angle is a fool's game.

http://earth.geology.yale.edu/~jjpark/Code/recfunk_pick_072107.f

(you didnt think that I would post a 1000-line Fortran code!)

The code needs the EISPACK library for matrix eigenvalue subroutines. Also Plotit and Sacread, which are the subject of prior posts.

Code is superceded on 8/16/07 to diversify the P-phase identification in the SAC header -- see later post

recfunk_pick_through.f -- read list of data files and retain only those which have a P-wave window marked for RF estimation

2007-07-21T13:53:00.001-04:00

c program recfunk_pick_through
c version of recfunk_pick that reads seismic record start times from
c the sac header, and writes to in_recpick the filenames that have been picked
c 10/13/05 JJP
c input file: in_recpick_raw
c output file: outfile (to be renamed in_recpick)
c
c xf77 -o /park/backup/bin/recfunk_pick_through recfunk_pick_through.f /park/backup/Plotxy/plotlib.a /park/backup/Ritz/eislib.a /park/backup/Ritz/jlib.a
c xf77 -o /Users/jjpark/bin/recfunk_pick_through recfunk_pick_through.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
......implicit real*4(a-h,o-z)
......real*8 ar,ai,el(12),aaa
......complex*8 zc,zero
......complex*8 afft,rf,rfs
......character*80 name,subtext,string
......character*28 title(2)
......character*10 cmps(3)
......character comp(3)
......common/nnaamm/iwflag
......common/npiprol/anpi
......common/stap2/tap(16400,16)
......common/taperz/ta(16400,8,3),tai(16400,8,3)
......common/stap0/ar(16400),ai(16400),bb(16400),time(99000)
......common/specs/afft(4100,8,3),rf(4100,2),sss(4100,6),crf(4100,2)
......common/saveit/rft(4100,799,2),crft(4100,2),rfs(2050,799,2),
x drfs(2050,799,2),bazs(799),epic(799),s2n(799),ip(799)
......common/saveit2/crfts(4100,799,2)
......common/datastuff/a(99000,3),pdat(16400),tdat(16400),drf(4100,2)
......common/header/ahead(158)
......common/distrib/xi(72,36),yi(72,36),sumij(72,36)
......common/chisq/chi(160000),enchi(160000),bbaz(160000),
x chim(200,2),enchim(200,2)
......dimension iah(158),chead(158),mmonth(12),fr(3),tt(3),ttt(3)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data comp/'r','t','z'/
......data cmps/'Radial ','Transverse','Vertical '/
......data title/'Radial Receiver Function ',
x 'Transverse Receiver Function'/
......con=180./3.14159265358979
......pi=3.14159265358979
......pih=pi/2.
......zero=cmplx(0.,0.)
......fmax=5.
......print *,'This code reads start of interval from sac header'
......print *,'it writes filenames of those files that have been picked'
c some default values (often superceded)
......fmax0=fmax
......fmin=0
......anpi=2.5
......nwin=3
......npad=16384
......nnf=npad/2
......do i=1,1000
...... crft(i,1)=0.
...... crft(i,2)=0.
......end do
......irec=0
......open(10,file='in_recpick_raw',form='formatted')
......open(11,file='in_recpick',form='formatted')
10 print *,'enter filename e.g. 1998.361.bhz'
......read(10,101) name
101 format(a)
......if(name(1:4).eq.'stop') go to 111
......subtext=name
c we assume that the last character is either 'r','t',or 'z'
......do i=80,1,-1
...... if(name(i:i).ne.' ') then
...... namlen=i
...... go to 98
...... endif
......end do
......print *,name
......print *,'no filename?'
......stop
98 continue
......do kcmp=3,3
...... name(namlen:namlen)=comp(kcmp)
...... print *,name
...... call sacread(name,a(1,kcmp),ierr)
...... if(ierr.eq.1) stop
...... dt=ahead(1)
...... tt(1)=0.
...... tt(2)=dt
...... tt(3)=dt
...... nscan=iah(80)
...... if(nscan.gt.99000) then
...... print *,'Careful! data length greater than 99000:', nscan
...... pause
...... endif
......end do
c naturally, you are in deep poop if dt is not the same for all records
......npre=10./dt
......npost=30./dt
......ntot=npre+npost
......ttt(1)=-npre*dt
......ttt(2)=dt
......ttt(3)=dt
......do i=1,ntot
...... time(i)=-ttt(1)-(i-1)*ttt(2)
......end do
......irec=irec+1
......if(irec.gt.799) then
...... print *,'max records exceeded'
...... stop
......endif
......baz=ahead(53)
......bazs(irec)=ahead(53)
......epic(irec)=ahead(54)
......print *,baz,'= back azimuth'
c ddf is cycles/sec
......ddf=1./(dt*npad)
......nf1=1
......nf2=fmax/ddf+1
......nf=nf2-nf1+1
......if(nf.le.0) stop
......fr(1)=0.
......fr(2)=ddf
......fr(3)=ddf
c FOR READING TIME INTERVAL FROM SAC HEADER
c ahead(9) is SAC's A parameter, for arrival time. Set by pressing A in plotpk
c......print *,'P window start time, duration (sec)'
c......read(10,*) pstart,postwin
......pstart=ahead(9)-ahead(6)-3.
......print *,'sach params b and a, pstart',ahead(6),ahead(9),pstart
......if(ahead(9).gt.0.1) write(11,101) subtext
......go to 10
111 write(11,101) name
......close(10)
......close(11)
......stop
......end

Fortran subroutine for X-window plotting of simple arrays

2007-07-21T13:45:00.000-04:00

see the following link

http://earth.geology.yale.edu/~jjpark/Code/plotit.html

examine_sac.f -- write a SAC macro to visually assess your body-wave data before receiver function estimation

2007-07-21T13:43:00.000-04:00

c program examine_sac
c 7/8/99 JJP - adapted for WEED output, after running event_prune & rotate_sac
c 4/1/02 altered to order the input filenames by epicentral distance
c...... to facilitate quality control (P and S wave pattern varies smoothly)
c
c xf77 -o /Users/jjpark/bin/examine_sac examine_sac.f /Users/jjpark/Ritz/jlib.a /Users/jjpark/Plotxy/plotlib.a
c to read filenames of data files (bhr,t,z) for a single event
c
......character*80 name,nami,namo*15,dumf(3)*2,bigbanner*85
c slat, slon are in radians, sdep is in km, slat=colatitude
c SAC headers are degrees and meters, respectively, and lat is lat
c ---> so must convert
......common/sour/slat,slon,sdep,sol(6),am0,strike,dip,slip
......common/ptitle/bigbanner,name,namo,dumf
c we read in data in scratch space
......common/sing2/a(240000)
......common/allfiles/nami(3,6000),gcarc(6000),gcord(6000),ip(6000)
......common/header/ahead(158)
......dimension iah(158),chead(158),mmonth(12),xx(3)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data dumf/'ae','an','az'/
......con=180./3.14159265358979
......open(7,file='infiles2',form='formatted') ! sacfiles 19*.bh[rtz] 20*.bh[rtz]
......open(18,file='outfile1',form='formatted') ! to be in_recfunk
......open(9,file='outfile',form='formatted') ! to be sac_examine
......open(13,file='outfile3',form='formatted') ! to be infiles3 (for pruning SKS data sets)
......write(18,103) '0.4'
......j=1
10 print *,'starting again'
......do k=1,3
...... read(7,101,end=998) nami(k,j)
......end do
......name=nami(1,j)
......print *,name
......call sacread(name,a,ierr)
......gcarc(j)=ahead(54)
......j=j+1
......go to 10
998 nevents=j-1
......close(7)
......xx(1)=1.
......xx(2)=1.
c order the events by epicentral distance
......call plotit(xx,gcarc,dum,nevents,'Events','event number',
x 'epicentral distance',1,0,0.05,2,21)
......do l=1,nevents
...... blob=200.
...... do m=1,nevents
...... if(blob.gt.gcarc(m)) then
............blob=gcarc(m)
............mm=m
...... endif
...... end do
...... ip(l)=mm
...... gcord(l)=gcarc(mm)
...... gcarc(mm)=200.
......end do
......call plotit(xx,gcord,dum,nevents,'Events','event number',
x 'epicentral distance',1,0,0.05,2,22)
......do j=1,nevents
...... name=nami(1,ip(j))
...... write(18,103) name(1:17),'?'
...... write(18,103)'117 75 75'
......end do
......write(18,103) 'stop'
......close(18)
......write(9,101) 'qdp 2000'
......do j=1,nevents
...... name=nami(1,ip(j))
........write(13,103) nami(1,ip(j))
........write(13,103) nami(2,ip(j))
........write(13,103) nami(2,ip(j))
...... write(9,103) 'read ',name(1:17),'?; chnhdr a 115 '
...... write(9,103)'rtrend; lp c 0.15 np 8; write lowr lowt lowz'
...... write(9,103) 'read ',name(1:17),'?'
...... write(9,103)'rtrend; hp c .5 np 8; write hir hit hiz'
...... write(9,103) 'readhdr lowz;listhdr gcarc baz evla evlo evdp ',
x 'kzdate'
c to pick arrival times for pulses, change hi?;plot1;pause
c to hi?;plotpk markall on; write over
...... write(9,103) 'read ',name(1:17),'? low? hi?;plot1;pause'
......end do
999 close(9)
......close(13)
101 format(a)
103 format(7a)
102 format(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f4.1)
105 format(22x,f8.2,1x,f9.2,1x,f4.0)
104 format(a,f10.3,a,f8.3,a,f10.2)
......print *,'to pick arrival times, change '
......print *,'to '
......stop
......end

subroutines sacread and sacout -- SAC-binary I/O

2007-07-21T13:29:00.000-04:00

These subroutines lack the neat feature of the new SAC release, which can detect and compensate for byte-swapped data. Beware the swapped byte! There are variants of the routines included that swap bytes, dogswork to be sure.

c......subroutine sacread_noswap(name,a,ierr)
......subroutine sacread(name,a,ierr)
c read a SAC-format file, dont swap the bytes
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......ierr=0
......open(8,file=name,access='direct',recl=512,err=999)
c read the 158-word header, can use its info
......read(8,rec=1,err=998)(ahead(i),i=1,128)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......irec=2
......read(8,rec=2,err=998)(ahead(i+128),i=1,30),(a(j),j=1,98)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......else......
...... print *,'reading greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......endif
......if(nrec.gt.3) then
...... do irec=3,nrec-1
...... read(8,rec=irec,err=998)(a((irec-3)*128+98+i),i=1,128)
...... end do
......endif
c......print *,'rec=',nrec,' Im trying to read the last partial record'
......read(8,rec=nrec,err=998)(a((nrec-3)*128+98+i),i=1,last)
......close(8)
......return
999 print *,'open error'
......ierr=1
......return
998 print *,'read error: reading',irec,' out of',nrec
......ierr=1
......print *,irec
......close(8)
......return
......end
......subroutine sacread_swap(name,a,ierr)
c......subroutine sacread(name,a,ierr)
c read a SAC-format file -- swap bytes for Intel mac for all numbers, not characters
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......ierr=0
......open(8,file=name,access='direct',recl=512,err=999)
c read the 158-word header, can use its info
......read(8,rec=1,err=998)(ahead(i),i=1,128)
c swap bytes
......do i=1,120
...... chd=chead(i)
........do j=1,4
........ chhead(j,i)=chd(5-j:5-j)
........end do
......end do
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......irec=2
......read(8,rec=2,err=998)(ahead(i+128),i=1,30),(a(j),j=1,98)
c swap bytes
......do i=1,98
...... ahd=a(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........a(i)=ahhd
......end do
c......print *,'rec=2'
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......else......
...... print *,'reading greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
...... last=nword-128*(nrec-1)
......endif
......if(nrec.gt.3) then
...... do irec=3,nrec-1
c...... print *,'rec=',irec
...... read(8,rec=irec,err=998)(a((irec-3)*128+98+i),i=1,128)
c swap bytes
...... do i=1,128
............ahd=a((irec-3)*128+98+i)
...... .....do j=1,4
..............chhd(j:j)=chd(5-j:5-j)
........ end do
........ a((irec-3)*128+98+i)=ahhd
...... end do
...... end do
......endif
......print *,'byteswapping read'
c this here is a big fat kluge --- I cant read the last partial record
c so I skip the last fragment and adjust the header to discard the points.
c......print *,'NOT READING ',last,' points in last record -- SACBUGFIX'
......read(8,rec=nrec,err=998)(a((nrec-3)*128+98+i),i=1,last)
c swap bytes
......do i=1,last
...... ahd=a((nrec-3)*128+98+i)
...... .do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........a((nrec-3)*128+98+i)=ahhd
......end do
c......nscan=nscan-last
c......iah(80)=nscan
......close(8)
......return
999 print *,'open error'
......ierr=1
......return
998 print *,'read error: reading',irec,' out of',nrec
......ierr=1
......print *,irec
......close(8)
......return
......end
c......subroutine sacout_noswap(name,a)
......subroutine sacout(name,a)
c write a SAC-format file
......real*4 ahead
......character*80 name
......character*4 chead(158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158)
......equivalence (ahead,iah),(ahead,chead)
......open(8,file=name,access='direct',recl=512)
c write the 158-word header
......write(8,rec=1)(ahead(i),i=1,128)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
......nscan=iah(80)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
......else......
...... print *,'writing greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
......endif
......write(8,rec=2)(ahead(i+128),i=1,30),(a(j),j=1,98)
......do irec=3,nrec
...... write(8,rec=irec)(a((irec-3)*128+98+i),i=1,128)
......end do
......print *,nrec,' records written'
......close(8)
......return
......end
c......subroutine sacout(name,a)
......subroutine sacout_swap(name,a)
c write a SAC-format file
......real*4 ahead,ahd,ahhd
......character*80 name
......character*4 chead(158),chd,chhd
......character*1 chhead(4,158)
......common/header/ahead(158)
......dimension a(1)
......dimension iah(158),output(158)
......equivalence (ahead,iah),(ahead,chead),(ahead,chhead(1,1))
......equivalence (chd,ahd),(chhd,ahhd)
......open(8,file=name,access='direct',recl=512)
......print *,(ahead(i),i=50,55)
......print *,(iah(70+i),i=1,10)
c write the 158-word header
c swap bytes
......do i=1,120
...... chd=chead(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........output(i)=ahhd
......end do
......do i=121,128
...... output(i)=ahead(i)
......end do
......write(8,rec=1)(output(i),i=1,128)
......nscan=iah(80)
......if(chead(145).ne.' GRN'.and.chead(145).ne.' grn') then
...... nword=nscan+158
...... nrec=(nword-1)/128+1
......else......
...... print *,'writing greens functions'
...... nword=6*nscan+158
...... nrec=(nword-1)/128+1
......endif
c swap bytes
......do i=1,30
...... output(i)=ahead(128+i)
......end do
......do i=1,98
...... ahd=a(i)
........do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........end do
........output(30+i)=ahhd
......end do
......write(8,rec=2)(output(i),i=1,128)
......do irec=3,nrec
...... iii=(irec-3)*128+98
...... do i=1,128
...... ahd=a(iii+i)
........ do j=1,4
........ chhd(j:j)=chd(5-j:5-j)
........ end do
........ output(i)=ahhd
...... end do
...... write(8,rec=irec)(output(i),i=1,128)
......end do
......print *,nrec,'byteswapped records written'
......close(8)
......return
......end

rotate_sac.f -- Preprocessing of SAC data for receiver function analysis

2007-07-21T13:20:00.000-04:00

c program rotate_sac
c 7/8/99 JJP -- adapted for JWEED output
c f77 -o /Users/jjpark/bin/rotate_sac rotate_sac.f /Users/jjpark/Ritz/jlib.a
c to read infiles1=filenames of data files (E,N,Z) for a single event
c cross-ref with file *events* of source locs and times
c writes a SAC macro to read three components of data, synchronize times, cut to equal length,
c put the event location in the header, and rotate the horizontals to radial and transverse.
c
c there must be three lines in infiles1 for every line in events, and the events must match the data
c take out the dots and replace with spaces before compiling
c
......character*80 name,nami,namo*15,dumf(3)*2,bigbanner*85,chh*51
c slat, slon are in radians, sdep is in km, slat=colatitude
c SAC headers are degrees and meters, respectively, and lat is lat
c ---> so must convert
......common/sour/slat,slon,sdep,sol(6),am0,strike,dip,slip
......common/ptitle/bigbanner,nami(3),namo,name,dumf
c we read in data in scratch space
......common/sing2/a(240000)
......common/header/ahead(158)
......dimension iah(158),chead(158),mmonth(12)
......equivalence (ahead,iah),(ahead,chead)
......data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
......data dumf/'ae','an','az'/
......con=180./3.14159265358979
......print *,'open(7,file=infiles1) ! sacfilenames'
......print *,'open(18,file=events) ! source locations and times'
......open(7,file='infiles1',form='formatted') ! sacfilenames
......open(18,file='events',form='formatted') ! source locations and times
......open(9,file='outfile',form='formatted')
......write(9,101) 'qdp 2000'
10 print *,'starting again'
c......if(jjc.eq.1) go to 112
......read(18,101,end=110) name
......print *,name
c......go to 10
c 112 continue
......read(name,102) jy,jm,jd,jh,jmi,ssec
......print *, jy,jm,jd,jh,jmi,ssec
c......nami(1:38)=name(43:80)
c......read(name,105) slat,slon,sdep
......read(name(44:80),*) slat,slon,sdep
......print *, slat,slon,sdep
101 format(a)
103 format(7a)
102 format(19x,i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,f6.3)
c here reads to column 45 -- if evdp can be read OK here, doesnt matter
c if the spacing of columns is perfect
105 format(22x,f8.2,f9.2,f4.0)
c 105 format(25x,3f8.2)
104 format(a,f10.3,a,f8.3,a,f10.2)
......leap=1-(mod(jy,4)+3)/4
......jd_julian=mmonth(jm)+jd
......if(jm.gt.2) jd_julian=jd_julian+leap
......write(9,103) 'cuterr fillz'
......sdep=sdep*1000.
......print *,slat,slon,sdep,' is depth in meters?'
......do k=1,3
...... read(7,101) nami(k)
......end do
......write(9,103) '/bin/rm ae an az'
......do k=1,3
...... name=nami(k)
...... ii=80
...... do while (name(ii:ii).eq.' ')
...... ii=ii-1
...... end do
c...... ii=ii-4 ! earlier SEED-output name convention
...... ii=ii-6 ! newer SEED-output name convention
...... print *,ii
...... if(name(ii:ii).eq.'Z') then
...... j=3
...... elseif(name(ii:ii).eq.'E'.or.name(ii:ii).eq.'2') then
...... j=1
...... elseif(name(ii:ii).eq.'N'.or.name(ii:ii).eq.'1') then
...... j=2
...... else
...... print *,'bad character: ',name(ii:ii)
...... endif
...... write(9,103) 'sc ln ',name(1:ii+7),dumf(j)
......end do
c LOVROK is the flag for allowing the user to overwrite a data file
......chh='synch;rmean;rtrend;p1;chnhdr LOVROK TRUE;write over'
......write(9,103)'read ae an az;',chh
......write(9,103) 'cut 0 600'
......write(9,103)'read;write over;cut off'
......do k=1,3
...... name=nami(k)
c read in the file to obtain its reference time
...... call sacread(name,a,ierr)
...... if(ierr.ne.0) go to 111 ! in case of mistyped filename
...... iy=iah(71)
...... id=iah(72)
...... ih=iah(73)
...... im=iah(74)
...... iss=iah(75)
...... ims=iah(76)
c compute origin time - year, month are assumed equal
...... torg=ssec-iss-float(ims)/1000.+
x......60.*(jmi-im+60*(jh-ih+24*(jd_julian-id)))
...... write(9,103) 'readhdr ',nami(k)
...... write(9,104) 'chnhdr evla ',slat,' evlo ',slon,' evdp ',sdep
...... write(9,104) 'chnhdr o ',torg
...... write(6,104) 'chnhdr o ',torg
...... write(9,103) 'writehdr'
......end do
......write(9,103) 'read ae an az'
......write(9,103) 'plot1'
......write(9,103) 'pause period 0.25'
......write(9,103) ' '
......write(9,103) 'read ae an'
......write(9,103) 'rotate to gcp reversed'
......write(9,103) 'read more az'
......write(9,103) 'plot1'
......write(9,103) 'pause period 0.25'
......write(9,103) ' '
......namo=name(1:15)
......write(9,103) 'write ',namo,'bhr ',namo,'bht ',namo,'bhz'
......go to 10
110 close(7)
......close(18)
......close(9)
......stop
111 print *,'error! ',name
......stop
......end

Synthetic Seismogram normalizations

2007-07-21T12:56:00.001-04:00

I have received queries about constants in the synthetic seismogram code anirec.f

In order to avoid round-off error and numerical overflow in numerical computations, it is conventional to nondimensionalize physical units in synthetic seismogram calculations. My codes use the convetion that I inherited from Freeman Gilbert, my dissertation advisor, and the UCSD Free Oscillation crowd. The rules for normalizing Mass, Length and Time are derived from three constraints.

density is normalized by rho-bar = average Earth density = 5515 kg/m^3

length is normalized by the mean Earth radius R=6.371*10^6 m

The last constraint is that pi*G=1, where pi=3.14159265358979d0 and G=6.6723*10^{-11} m^3/(kg-s^2) is the Gravitational constant.

From the public free-oscillation code MINEOS

c the model has been normalised such that a density of 5515 mg/m**3 = 1 ;
c pi*g=1 where g is the grav constant (6.6723e-11) and the radius of the
c earth (rn=6371000 m) is 1. these normalizations result in
c acceleration normalisation = pi*g*rhobar*rn = 7.365081 m/s^2
c velocity normalisation = rn*sqrt(pi*g*rhobar)= vn = 6850.0396 m/s
c frequency normalization = vn/rn = sqrt(pi*g*rhobar) = 1.075190645e-3 s^1

This last constraint leads to a normalization constant for frequency of REN=1.075190645d-3. I have found a slightly different value for G quoted on the web, which is understandable considering that the seismological constants were set when G was not known as precisely as now. There can be danger in mistyping the numbers. One colleague of mine lost days and much sleep when the rho-bar constant was mistyped 5517 kg/m^3 and a 10^{-4} error crept into calculations that should have been good to double precision. Believe me, no one was lauging when the error was discovered!

Converting JWeed event files to julian days

2007-07-21T12:40:00.000-04:00

JWeed event files identify earthquakes in month-day format, but seismic data filenames follow Julian-day format. Unless you have a magic patch of brain cells that converts the calendar to julian days, this is a big hassle when matching files with events. I adapted a little fortran code to convert, and some awk lines that convert a JWeed event file into a format that is used by other codes to pre-process body wave data. Here are the awk lines:

sed 's/,/ /g' events*.events | sed 's[/[ [g' >! blob.dat
awk '{print "echo " $3 " " $4 " " $2 " >! blobb; echo " $0 " >>blobb ; Julian1 < blobb >> ev_julian " }' blob.dat >! c_events
sed 's/xxxx yyyyy//g' c_events >! c_ev

the JWeed event filenaming convention (for me) is events_YYYY.events, with variants. The first step strips the syntax separating numbers in the events file, so that awk can parse the numbers. The second line creates a command file c_events that uses Julian1 to read the output of the awk step. The last step takes out some superfluous matter from c_events, then you "source c_ev" to make ev_julian. I rename this file events.julian as a template for the file "events" that each station needs in preprocessing (see previous post).

Oh yes, there is a Fortran code

c program Julian1
c f77 -o /Users/jjpark/bin/Julian1 Julian1.f
character string*80
dimension mmonth(12)
data mmonth/0,31,59,90,120,151,181,212,243,273,304,334/
c print *,'enter month, day, year'
read(5,*) jmo,jd,jy
leap=1-(mod(jy,4)+3)/4
jd_julian=mmonth(jmo)+jd
if(jmo.gt.2) jd_julian=jd_julian+leap
c print *,'julian day = ',jd_julian
read(5,101) string
print 102,jy,jd_julian,string
101 format(a)
102 format(2i5,4x,a)
stop
end

Installing JWeed on Linux

2007-07-20T16:04:00.000-04:00

There is a nasty glitch in the Linux installation for JWeed. The IRIS website downloads install.bin, but it has two bad lines. FIRST, run the installer thusly to fix one bad line without changing the checksum

cat install.bin | sed "s/export LD_ASSUME_KERNEL/#xport LD_ASSUME_KERNEL/" > temp
sh temp

Then comment out lines 1335 and 1375 of JWeedv2.9.8, which is the executable script that embodies JWeed. These are the lines

1335 export LD_ASSUME_KERNEL=2.2.5

1375 LD_ASSUME_KERNEL=2.2.5

After these changes, there is still a little complaint about an extra backslash, which I tracked down to line 1343:
split ($3, ver, "[\._-]")

I took the backslash out and the code still worked.

Thanks to Alex Nikulin and Garrett Leahy for the tip!

Preparing RETREAT body waves for Receiver-Function Estimation

2007-07-20T14:11:00.000-04:00

Work performed on an IntelMac laptop -- computer codes and scripts will be posted in successive blog posts. This post documents the preparation of body wave data from the RETREAT experiment. You can use the concepts to guide the preparation of data for other stations.

The directory Teleseismic has subdirectories PIIR, VLC, VOLR, etc, one for each station name. The station subdirectories have data files in SAC format, but without the source location in the header. There are three steps in preparing a station's data for processing. First, we associate events with data-file triplets. The extraction by RDSEED from the SEED files should have synchronized the data files for each component of motion. This first step can be tedious, however, because there will be more events than datafile triplets whenever the station was down, and more data files than events when the data files are fragmented. The goal is to prune the files "events" to only the events that have data, and the file "infiles1" to be the triplets of datafile names -- infiles1 will have three times as many lines and events, natch.

cd PIIR (from the Teleseismic directory)
ls *.SAC > infiles_raw
cp infiles_raw infiles1
cp ../events.julian events
te events &
te infiles1 &

events_julian is a master file of events used to request data from the IRIS DMC for the experiment, based on a JWeed event file:

2003 356 MHDF 2003 12 22 19:15:56.000 35.706 -121.102 7.6 MB 6.1 MS 6.4
2003 359 MHDF 2003 12 25 07:11:11.590 8.416 -82.824 33 MB 6.0 MS 6.4
2003 359 MHDF 2003 12 25 20:42:33.720 -22.252 169.488 10 MB 6.3 MS 6.3
2003 360 MHDF 2003 12 26 01:56:52.440 28.995 58.311 10 MB 6.0 MS 6.8
2003 360 MHDF 2003 12 26 21:26:04.100 -22.273 169.314 10 MB 6.1 MS 6.8
2003 361 MHDF 2003 12 27 16:00:59.450 -22.015 169.766 10 MB 6.1 MS 7.1
2003 361 MHDF 2003 12 27 22:38:01.880 -21.672 169.835 10 MB 5.8 MS 6.7

Note that the month/day/year format of the JWeed event file is replaced by year julianday, which matches the output filename format of rdseed.

The user now opens infiles1 and events in a text editor and synchronizes their lines. It is smart to keep an Xterminal open to run SAC or pql. Sometimes there will be more than three datafiles associated with an event start time, and you must display their contents in order to choose the three files to use, or else to merge some files.

pql yyyy.ddd.* &

is a good command for this. Or

echo "read yyyy.ddd.*; qdp 2000; p1" | sac

because you can retrieve previous commands in a Unix shell and replace the yyyy.ddd for different cases.

The second step is to run the code "rotate_sac." this code reads the info in infiles1 and events. It generates a SAC macro file that will read the data triplets, put the earthquake hypocentre into the SAC headers, cut the files to 600-s length and rotate the horizontals to radial and transverse direction. You run

rotate_sac
/bin/mv outfile sac_rotate
echo exit >> sac_rotate
sac < sac_rotate

and the files *yyyy.ddd.hh.*BH[ENZ]*.SAC are converted to truncated and rotated data files "yyyy.ddd.hh.mm.bh[zrt]" Now comes the third step of processing. The code examine_sac reads the data files and composes a SAC macro to visualize the data in raw and bandpassed traces, in order to assess the signal to noise ratio of the P wave.

ls *bh? >! infiles2
echo 3 >! c_c_c
examine_sac < c_c_c
/bin/mv outfile1 in_recfunk
/bin/mv outfile sac_examine
sed -e 's/plot1;pause/plotpk markall on; writehdr over/' sac_examine >! sac_examine_pick
sac < sac_examine_pick

The SAC-macro sac_examine displays data from each event in raw, lowpassed and highpassed traces. Oftentimes a Pwave is evident in some bandpasses and not others -- I typically keep events where the low or high bandpass has SNR>2, but SNR.le.1 elsewhere. The last "sed" command converts this display macro into a phase-picking macro. "sac_examine_pick" asks you to click "A" on the desired start of a RF data window (the duration of the window is user-input in the RF codes, but the start time of analysis is most convenient if stored in the SAC header). The SAC macro marks "A" in the lowpassed traces at the point 120s from the start of the record, the typical timing of the P wave. A PKP or PP wave will arrive later than this, except for events that are so far away that JWEED neglected to use the Pdiff as a reference. GCARC~150 is the transition, but Im not certain where exactly this transition occurs. For P and Pdiff selection, you can leave the cursor close to 120s after trace onset, and adjust to capture the P wave with a 5-10 seconds to spare. For the PKP phase, you must move the cursor to later in the record.

The data-evaluation and marking exercise is done with the A and Q keys.
If the P or PKP wave is determined "analyzable" click A to mark the start time of the data window, then Q to move to the next record. SAC marks all data traces with the A-time. If you make a mistake in selecting A, either too early or too late, just mark it again.
If the P or PKP wave is weak or bad (e.g. the PP overlaps the PKP too much), then simply click Q and SAC will leave the datafile unmarked. This flags the event and datafiles to be discarded. If you click Q by accident before clicking A before a good P phase, you must return to the data later by running the relevant lines of sac_examine_pick. Write down the event time (yyyy.ddd.hh) and go back to it later.
examine_sac takes the data files and orders them in increasing epicentral distance. This way the user will assess the data first from nearby earthquakes with high-frequency P, moving on to low-frequency Pdiff, to an interval where only the largest quakes have visible Pdiff, and then to GCARC~120 where the PKP phase precedes the PP phase by enough time to be useful for RF analysis. Finally the time-reference for JWeed shifts and the data interval has PKP/PKIKP starting at 120s.

After the desired P waves have been identified with A-markers in their SAC headers, it is time to prune the list of events down to the chosen few. To prune the in_recfunk file,

grep bh in_recfunk >! in_recpick_raw
echo stop >> in_recpick_raw
recfunk_pick_through

--> reads in_recpick_raw and writes in_recpick, a list of wildcarded filenames that includes only those whose headers define a data window to analyze. Then this in_recpick can be used by recfunk_pick to compute simple RF sweeps (no migration, only epicentral and back-azimuths bin-sums).

recfunk_pick

There are GMT scripts in the Teleseismic directory, some files from recfunk_pick can be plotted to see what the station RFs are like:

../Plotss outr_baz.grid outt_baz.grid 0.75 plot_baz
cp plot.ps plot_baz.ps

../Plot_epi3 outr_epi.grid outt_epi.grid 0.75 plot_epi
cp plot.ps plot_epi.ps

Using JWeed to read seismic data from the RETREAT Project (Post 1)

2007-07-20T13:57:00.000-04:00

Notes on getting RETREAT data from the IRIS DMC (6/13/07) -- Experience working on an Intel Mac laptop

FIRST: obtain JWeed from the IRIS website. Install on your computer.

SECOND: obtain ccrypt from sourceforge on the web. There are versions that have executables for many processors and operating systems, but it might be safest to get the source code and compile it in your own environment. First "./configure" then "make" then "sudo make install" on a Mac.

Running JWeed

Three items to remember: 1) the RETREAT network code is YI. There are some other stations in Antarctica that are included in network YI. You might be able to avoid these in JWeed by restricting the stations to a box around Italy. 2) the definition of YI is time-variable. If you ask for the station info when JWeed is set for a time interval in 2007, you wont get the RETREAT stations. 3) the RETREAT and GAIA stations record on the BH? channels, and the INGV network stations record on the HH? channels. You must request both channel types to get data from both sets of stations.

Dont try to move back and forth in JWeed while you are retrieving something from the internet. In particular, the download of stations from the database is tedious, and there is temptation to wander in JWeed. However, I have found that JWeed will halt when you move to a different tab in the GUI. It is smarter to use the Network DataCenter once for each time interval of interest, then use the JWeed-written station file for later requests.

To get the INGV HH? stations, Ive wasted a day trying. HH? channel is not a channel normally found by JWeed. Also the IRIS database cant always find the INGV stations in JWeed. It is a pisser. The key is that JWeed wont find a station that does not have channel that is being requested. So if you ask for BH? channels, only the PASSCAL and GAIA stations appear in YI, not the INGV stations. I find that the request works if you ask for BOTH B* and HH* channels.

The other JWeed "feature" is that you cant create multiple breqfast files without exiting the software. Either the YI network decreases down to a single station on the Antarctic ice cap, the breqfast file hangs, or something else surprises you. Give up! Exit JWeed and start over!

If you have event files and station files, it is possible to run JWeed without accessing the outside DHI servers. However, I had to turn off the internet access for the computer to avoid problems. Also, you cant run JWeed multiple times to create multiple files. Ive had to close JWeed and start over every time.

Here is what seems to work:

1) set IRIS netDC for Jweed in first tab.

2) Set the whole world in the second tab -- it is the default box. You can also specify an epicentral range around a map point. Use the crosshairs icon to the left

3) in third tab, specify magnitude range, depth range, MB MS ML as the magnitudes to use for body-wave records, load the catalogs if they are not there already, choose WHDF (maybe) and specify the time period for gathering events (VERY important!)

4) Save events if you like

5) Map the events, an operation that pops open the second tab again -- in the map specify a small rectangle over Italy for choosing stations

6) in the 4th tab, access the DHI to bring down a list of networks and stations from the server. This can take a long time. When "finished" open the YI folder to verify that the stations desired are there. If only the Antarctic station is there, search the DHI again. If only BH? stations are there and not the INGV stations, start over if you want the INGV stations.

7) Select the YI folder and map the stations. save them to a file if you wish.

8) Go to the generate-request tab. click the "listing" button, rather than the "map" button. wait until the summary file appears in the text window. The program prompts for a file name and write out the summary file automatically. JWeed write the file to the JWEED_HOME directory. Make sure that you have set this environmental variable.

9) Go to the mail-request tab. Fill in the info, choose the Breqfast-IRIS tab to specify the IRIS DMC. Do NOT email the request directly, save it to a file. This allows you to fix any syntax issues easily.

Running CCRYPT

the routine ccrypt has better encryption protection than the old Unix version of crypt. This is great, but IRIS encodes with the old version. Therefore you must use the "u" option:

ccrypt -u < file.seed.crypted > file.seed

and the routine prompts you for the crypting key. I put alias crypt "ccrypt -u" into my .cshrc file

Using Breq_fast

JWeed will generate a breq_fast request. It is smart to save the breqfast file on your home computer in order to fix syntax problems. Some syntax issues crop up intermittently. For instance, JWeed will let you specify multiple channels in the request
e.g. BHZ HHZ, but in the last few days the two items in the breqfast file are not parsed by the DMC software, only the first one is. Parsing the multiple-channel request by changing the breqfast file to "BHZ HHZ" with quote marks WONT WORK. For the YI case you can try to specify ?H? as the channel and the request will obtain both HH? and BH?. BUT, watch the last number before the channel identifier. It should be the number 1, else there will be a syntax error. However, the ?H? trick works in some cases and not others. Therefore I have requested BH? and HH? separately. One can use the same breqfast file twice, since no station has both BH and HH channels.

A number of tests failed to work. I hypothesize that Breqfast wildcarding wont allow '?', only '*' In any case, the last time I sent a breakfast request with the '*' wiildcard and two channel items, the DMC responded as documented. AAAAAUUUUUGGGGHHH!!!!!

Running rdseed on the SEED files

This seems to be OK. The use of summary files is important, because these will synchronize start times for the three components of motion, if possible. RDSEED is not transferring event informationj from the summary file, so there are no epicentral distances and azimuth values (or correct ones!) in the SAC header. It is possible to run rdseed in command-line mode:

rdseed -d -f filename.seed -x filename.summary

Welcome to this Blog

2007-07-20T13:50:00.000-04:00

JParkCodes is an expedient for me to disseminate information about computer codes used by the Yale University Seismology Group. It will take me some time to post useful code for folks in Yale and colleagues worldwide. Send email with queries, and add comments to posted programs and shell-scripts to flag bugs and suggest improvements.