[cig-commits] r17968 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . UTILS UTILS/SeismiQuery UTILS/lib UTILS/seis_process
liuqy at geodynamics.org
liuqy at geodynamics.org
Wed Feb 23 19:40:55 PST 2011
Author: liuqy
Date: 2011-02-23 19:40:54 -0800 (Wed, 23 Feb 2011)
New Revision: 17968
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/dhi2mseed.py
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/iris.txt
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/config.h
Removed:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/sph2utm_Qinya_Liu.tar.gz
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/convolve_stf.c
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/rotate.pl
seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90
Log:
Update processing scripts.
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/dhi2mseed.py
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/dhi2mseed.py (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/dhi2mseed.py 2011-02-24 03:40:54 UTC (rev 17968)
@@ -0,0 +1,98 @@
+#!/usr/bin/python
+
+# requesting data for events and convert them to sac format
+# written by Qinya Liu, Nov 2009, for global shakemovie website
+
+import os, sys, glob, getopt, datetime
+
+usage='\nUsage: dhi2mseed.py -m CMT -d duration (min) [-j DHI2mseed.jar_file -p iris.txt -s station_file -o out_datadir]\n'
+
+try:
+ opts,arg=getopt.getopt(sys.argv[1:],'j:m:d:p:s:o:')
+except getopt.GetoptError:
+ sys.exit(usage)
+if (len(opts) == 0):
+ sys.exit(usage)
+
+# dir where _GSN.dataless or II/IU.dataless is placed
+dataless_dir='/data2/Datalib/stations/'
+jar_file='/data2/Datalib/iris_request/DHI2mseed.jar'
+iris_file='/data2/Datalib/iris_request/iris.txt'
+sta_file='/data2/Datalib/stations/STATIONS_LH_II_IU'
+out_dir='data'
+
+for o,a in opts:
+ if o == '-j':
+ jar_file=a
+ if o == '-m':
+ cmt_file=a
+ if o == '-d':
+ dur=float(a)
+ if o == '-p':
+ iris_file=a
+ if o == '-s':
+ sta_file=a
+ if o == '-o':
+ out_dir=a
+
+if (not os.path.isfile(jar_file)):
+ sys.exit(jar_file+' does not exist')
+if (not os.path.isfile(cmt_file)):
+ sys.exit(cmt_file+' does not exist')
+if (not os.path.isfile(iris_file)):
+ sys.exit(iris_file+' does not exist')
+if (not os.path.isfile(sta_file)):
+ sys.exit(sta_file+' does not exist')
+if (not os.path.isdir(out_dir)):
+ os.mkdir(out_dir)
+
+# write station request file
+f=open(out_dir+'/scnl.list','w')
+lines = file(sta_file, 'r').readlines()
+for i in range(0,len(lines)):
+ sp = lines[i].split()
+# print "%s LH* %s *" % (sp[0], sp[1])
+ f.write("%s LH* %s *\n" % (sp[0], sp[1]))
+f.close()
+
+# calculate request start time
+[year,mon,day,hr,min,sec]=os.popen("awk 'NR == 1 {print $0}' "+cmt_file).readline().split()[1:7]
+
+starttime=(datetime.datetime(int(year), int(mon), int(day), int(hr), int(min), int(float(sec)))+datetime.timedelta(minutes=-3)).isoformat(' ')
+
+command='java -jar '+jar_file+' -starttime "'+starttime+'" -duration '+str(dur)+' -p '+iris_file+' -f '+out_dir+'/scnl.list -o '+out_dir
+
+print '\n**** requesting data from IRIS ***'
+print command+'\n'
+if (os.system(command+' |grep writingi 1>&2') != 0):
+ sys.exit('Error requesting data')
+
+print '\n**** convert miniSEED files to sac files ****'
+if (not os.path.isdir(dataless_dir)):
+ sys.exit(dataless_dir+' does not exist')
+
+# use different dataless seed file for files from different networks
+net_input={};
+for file in glob.glob(out_dir+'/*.mseed'):
+ bfile=os.path.basename(file)
+ net=bfile.split('.')[2]
+ if (net in net_input.keys()):
+ net_input[net]=net_input[net]+bfile+'\n\n\nd\n\n\n\n\n\n\n\n\n3\n\n\n\n\n\n'
+ else:
+ net_input[net]=bfile+'\n\n\nd\n\n\n\n\n\n\n\n\n3\n\n\n\n\n\n'
+
+for net,ni in net_input.iteritems():
+ dataless_file=dataless_dir+'/'+net+'.dataless'
+# dataless_file=dataless_dir+'/_GSN.dataless'
+ if (not os.path.isfile(dataless_file)):
+ print dataless_file+' does not exist'; exit()
+ print 'Using dataless file ' + dataless_file+' for network '+ net
+
+ os.environ['ALT_RESPONSE_FILE']=dataless_file
+
+# print 'cd '+out_dir+'; rdseed <<EOF\n'+ni+'Quit\nEOF\n' # debug
+ if (os.system('cd '+out_dir+'; rdseed <<EOF 2> /dev/null | grep Writing 1>&2 \n'+ni+'Quit\nEOF\n') != 0):
+ sys.exit('Error converting miniseed to sac for '+net)
+
+# move original mseed files to a sub-directory
+os.system('rm -f scnl.list; cd '+out_dir+'; mkdir -p mseed_files; mv *.mseed mseed_files/')
Property changes on: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/dhi2mseed.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/iris.txt
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/iris.txt (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/SeismiQuery/iris.txt 2011-02-24 03:40:54 UTC (rev 17968)
@@ -0,0 +1,5 @@
+NameService=corbaloc:iiop:dmc.iris.washington.edu:6371/NameService
+NetworkDNS=edu/iris/dmc
+NetworkDC=IRIS_NetworkDC
+SeismogramDNS=edu/iris/dmc
+SeismogramDC=IRIS_DataCenter
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/config.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/config.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/config.h 2011-02-24 03:40:54 UTC (rev 17968)
@@ -0,0 +1,4 @@
+// 64-bit
+typedef int sac_int_t;
+// 32-bit
+// typedef long int sac_int_t;
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/convolve_stf.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/convolve_stf.c 2011-02-24 03:35:46 UTC (rev 17967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/lib/convolve_stf.c 2011-02-24 03:40:54 UTC (rev 17968)
@@ -193,11 +193,11 @@
/* creat source time function time series */
if(min_nhdur * dt / 2 > hdur) {
- fprintf(stderr,"The half duration %f is too small to convolve\n", hdur);
+ fprintf(stderr,"The half duration %f is too small \n", hdur);
return 1;
}
- nstf = (int)ceil(2 * hdur/dt)+1;
+ nstf = (int)ceil(6 * hdur/dt)+1;
hstf = (nstf-1)*dt/2;
if((stf = (float *) malloc(nstf*sizeof(float))) == NULL) {
@@ -206,21 +206,22 @@
}
if(cstf == 't') {
/* triangular */
- for (i=0; i<nstf/2; i++)
- stf[i] = i * dt / (hstf * hstf);
- for (i=nstf/2; i<nstf; i++)
- stf[i] = (2 * hstf - i * dt)/ (hstf*hstf);
- } else {
- /* gaussian */
const float decay_rate = 1.628;
- float alpha = decay_rate / hdur;
- float divisor = sqrt(4*atan(1.0) /*pi*/);
- for (i=0; i<nstf; i++) {
- float tao_i = fabs(i*dt - hstf);
- stf[i] = alpha * exp(- alpha*alpha* tao_i*tao_i) / divisor;
- }
+ hdur = hdur/decay_rate;
+
+// for (i=0; i<nstf/2; i++)
+// stf[i] = i * dt / (hstf * hstf);
+// for (i=nstf/2; i<nstf; i++)
+// stf[i] = (2 * hstf - i * dt)/ (hstf*hstf);
}
+ printf("convolved by a gaussian function with hdur %f\n",hdur);
+ float alpha = 1. / hdur;
+ float divisor = sqrt(pi);
+ for (i=0; i<nstf; i++) {
+ float tao_i = fabs(i*dt - hstf);
+ stf[i] = alpha * exp(- alpha*alpha* tao_i*tao_i) / divisor;
+ }
/* creat convolution time series */
convolve(&conv,&nconv,data,npts,stf,nstf);
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl 2011-02-24 03:35:46 UTC (rev 17967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_data.pl 2011-02-24 03:40:54 UTC (rev 17968)
@@ -7,7 +7,8 @@
sub Usage{
print STDERR <<END;
- process_data_new.pl -m CMTFILE -l Start/End -t Ts/Tl -f -i -R resp_dir
+ process_data_new.pl -m CMTFILE -l Start/End
+ -t Ts/Tl -h hdur -f -i -R resp_dir
-p -s Sps -P n/p -T Taper_width -c -y t1/v1/...
-x Ext -d OutDir
data_sac_files
@@ -17,9 +18,10 @@
-a STAFILE -- add station information from STAFILE, use '-a none' to get
information from Vala''s station file.
-l Start/End -- start and end of record from o
- -t Ts/Tl -- shortest and longest period
- -f -- do not apply butterworth filter before deconvolution
- -i -- transfer record to displacement
+ -t Ts/Tl -- shortest and longest period (used in filters and ins transfer)
+ -i -- transfer record to displacement based on instrument response
+ -f -- transfer record only, no extra filtering
+ (possible scenarios: -i/-R -t; -i/-R -t -h; -i/-R -t -f;
-R resp_dir -- transfer record using RESP file in resp_dir
-p -- pick event p and s first arrival into sac headers (t1 and t2)
-x Ext -- extension of new file name
@@ -35,32 +37,35 @@
Notice:
1. We require that polezero files in the same directory as the sac
data files which is generally satisfied. We require that resp dir
- be specified even if it is current(.). All needed info is
+ be specified even if it is current(.). All the needed info is
taken from sac headers
2. Origin time is set to PDE + time_shift (given by the CMTSOLUTION)
3. The displacement outputs after -i option are in the unit of meters
4. For BH? components, please set -s 20, otherwise interpolation to
1 sample/second will be performed
-
- NOTE: Please make sure that SAC, saclst and IASP91 packages are installed properly on
- your system, and that all related env variables are set properly before
- running the script.
- Qinya Liu, originally written in Oct 2002; updated in Jan 2011
+ Qinya Liu, Originally written in Oct, 2002; updated in Nov 2009
END
exit(1);
}
+$saclst="saclst";
+$phtimes="./phtimes.csh";
+
if (@ARGV == 0) { Usage(); }
-if (!getopts('m:a:l:t:fiR:px:d:s:P:T:y:c')) {die('Check input arguments\n');}
+if (!getopts('m:a:l:t:iR:px:d:s:P:T:y:ch:f')) {die('Check input arguments\n');}
if ($opt_m and not -f $opt_m) {die("Check if file $opt_m exists\n");}
-if ($opt_a and not -f $opt_a) {$opt_a="/opt/seismo-util/data/STATIONS_new";}
+if ($opt_a and not -f $opt_a) {$opt_a="/opt/seismo/datalib/stations/STATIONS_LH";}
if ($opt_R and not -d $opt_R) {die("Check if dir $opt_R exists\n");}
+# check the use of -h to with -I
+if ($opt_h and !($opt_I or $opt_R)) {die("-h and -I/-R should appear in pairs\n");}
+if ($opt_h) {$hdur_input=$opt_h;}
if ($opt_t) {($tmin, $tmax) = split(/\//,$opt_t);
$f1 = 1./$tmax; $f2=1./$tmin;}
+
if ($opt_d) {
$out_dir=$opt_d; if (not -d $out_dir) {mkdir($out_dir,0777);}
}
@@ -72,11 +77,8 @@
}
if (!$opt_T) {$opt_T = 0.05;}
if ($opt_l) {($lmin,$lmax) = split(/\//,$opt_l);} else {$lmin = 0; $lmax = 3600;}
-
if ($opt_f and not ($opt_R or $opt_i)) {die("-f option goes together with -i or -R\n");}
-$saclst="saclst";
-$phtimes="./phtimes.csh";
$eps=1e-5; $undef=-12345.0; $cundef="-12345";
@@ -132,7 +134,7 @@
print SAC "echo on\n";
print SAC "r $outfile\n";}
- if ($opt_l){ # cut record
+ if ($opt_l){ # cut record and rtrend and rmean
print " Cut the record from o+$lmin to o+$lmax\n";
(undef,$tmp_o)=split(" ",`$saclst o f $outfile`);
if (abs($tmp_o-$undef) < $eps) {die("O has not been set to cut record\n");}
@@ -146,18 +148,10 @@
print SAC "echo on\n";
print SAC "r $outfile\n";}
-
- if ($opt_t and not $opt_f) {# filter records
- print " Filter record at periods $tmin and $tmax\n";
- print SAC "rtrend\n rmean\n taper width $opt_T\n";
- printf SAC ("bp n %4d p $pass cor %12.6f %12.6f\n",$poles,$f1,$f2);
- printf SAC " rtrend\n rmean\n taper width $opt_T\n";}
-
if ($opt_i or $opt_R) { # transfer instrument response
if (! $opt_t) {die(" Specify bandpass range by -t\n");}
$f0=$f1*0.8;
$f3=$f2*1.2;
- print SAC "rtrend\n rmean\n taper width $opt_T\n";
(undef,$network,$sta,$comp,$khole)=split(" ",`$saclst knetwk kstnm kcmpnm khole f $outfile`);
if ($network eq $cundef or $sta eq $cundef or $comp eq $cundef or $khole eq $cundef) {
die("No network station name or component name or khole defined in $outfile\n");}
@@ -177,12 +171,28 @@
printf SAC ("trans from evalresp fname $respfile to none freq %12.6f%12.6f%12.6f%12.6f \n",
$f0,$f1,$f2,$f3);
printf SAC "mul 1e-9\n"; }}
- printf SAC " rtrend\n rmean\n taper width $opt_T\n";
}
- if ($opt_s) {print SAC "interp delta $dt\n";
- print SAC "w over\n";}
+ if ($opt_t and not $opt_h and not $opt_f) {# butterworth filter
+ print " Filter record at periods $tmin and $tmax\n";
+ print SAC "rtrend\n rmean\n taper width $opt_T\n";
+ printf SAC ("bp n %4d p $pass cor %12.6f %12.6f\n",$poles,$f1,$f2);
+ printf SAC " rtrend\n rmean\n taper width $opt_T\n";}
+ elsif ($opt_h and not $opt_f) { # gaussian filter
+ print SAC "rtrend\n rmean\n taper width $opt_T\n";
+ print SAC "w $outfile \nquit\n";
+ close(SAC);
+ system("convolve_stf t $hdur_input 1 $outfile");
+ if ($? != 0) {die("Check if |convolve_stf t $hdur_input 1 $outfile| was successfully run\n");}
+ system("mv -f ${outfile}.conv $outfile");
+ if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
+ else {open(SAC,"|sac >sac_out.log ") || die("Can't open sac\n");}
+ print SAC "echo on\n r $outfile\n";
+ }
+
+ if ($opt_s) {print SAC "interp delta $dt\n";}
+
if ($opt_p) { # add p and s arrival info
print " Add first P and S arrival information\n";
(undef,$gcarc,$edep)=split(" ",`$saclst gcarc evdp f $outfile`);
@@ -208,7 +218,6 @@
printf SAC ("ch $h1 %12.2f\n",$v1);
print SAC "ch $k1 $h1\n";}}
-# print " write file $outfile\n";
print SAC "w $outfile\n";
print SAC "echo off\nquit\n";
close(SAC);
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl 2011-02-24 03:35:46 UTC (rev 17967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/process_syn.pl 2011-02-24 03:40:54 UTC (rev 17968)
@@ -9,7 +9,7 @@
Usage: process_syn_new.pl
-S -m CMTFILE -h -o offset -lStart/End -tTmin/Tmax -f -P n/p
- -i Dir -A amp_factor -p -a STAFILE -s sps -y t0/V0/t1/V1...
+ -i Dir/data-comp -A amp_factor -p -a STAFILE -s sps -y t0/V0/t1/V1...
-c -d OutDir -x Ext synthetics_files
where
-S already sac file, will not apply ascii2sac.csh
@@ -20,7 +20,8 @@
-t Tmin/Tmax -- shortest/longest period of bandpass filter
-f -- apply transfer filter instead of butterworth filter for -t
-P -- number of poles and passes in butterworth filter(default 4/2)
- -i Dir -- convolve synthetics with instrument response in Dir
+ -i Dir/data-comp -- convolve synthetics with instrument response for
+ given data channel in Dir
-A amp_factor -- multiple all files with an amplification factor
-p -- pick P and S arrivals according to iasp91 model to t1/t2 headers
-a STAFILE -- add station information from STAFILE, use "-a none"
@@ -57,10 +58,14 @@
if(not defined $pass or $pass<1){$pass=2;}}
if ($opt_l) {($lmin,$lmax) = split(/\//,$opt_l);}
else {$lmin = 0; $lmax = 3600;}
-if ($opt_a and not -f $opt_a) {$opt_a="/opt/seismo/data/STATIONS_new";}
+if ($opt_a) {$stafile=$opt_a;
+ if (not -f $stafile) {$sta_file="/opt/seismo/data/STATIONS_new";}}
if ($opt_o and not $opt_m) {die("Specify centroid time first\n");}
if ($opt_d and not -d $opt_d) {die("No such directory as $opt_d\n");}
-if ($opt_i and not -d $opt_i) {die("No such directory as $opt_i\n");}
+if ($opt_i) {
+ ($ins_dir,$cmpd) = split(/\//,$opt_i);
+ if (not -d $ins_dir) {die("No such resp directory as $ins_dir\n");}
+ if (not defined $cmpd) {$cmpd="LH";} @cmpd=split(//,$cmpd);}
if ($opt_A) {if ($opt_A !~ /^\d/) {die("-A option should be numbers\n");}}
if ($opt_f) {
if (not $opt_t) {die("-t tmin/tmax has to be specified for -f\n");}
@@ -70,8 +75,9 @@
$eps = 0.1;
$saclst="saclst";
-$phtimes="./phtimes.csh";
+$phtimes="phtimes.csh";
$asc2sac="ascii2sac.csh";
+$convolve_stf="convolve_stf";
#if (! -e $saclst) {die(" No $saclst file\n");}
#if (! -e $phtimes) {die("No $phtimes file\n");}
@@ -125,7 +131,7 @@
if ($opt_h and $hdur > $min_hdur) { # convolve source time function
print SAC "quit\n";
close(SAC);
- system("/opt/seismo-utils/bin/convolve_stf g $hdur $outfile");
+ system("$convolve_stf t $hdur $outfile");
system("mv ${outfile}.conv ${outfile}");
if ($opt_c) {open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
print SAC "echo on\n r $outfile\n";}
@@ -137,21 +143,16 @@
if ($comp=~/E/) {print SAC "ch cmpaz 90 cmpinc 90\n";}
elsif ($comp=~/N/) {print SAC "ch cmpaz 0 cmpinc 90\n";}
elsif ($comp=~/Z/) {print SAC "ch cmpaz 0 cmpinc 0 \n";}
- if (! -f $opt_a ) {die(" No such files: $opt_a");}
- ($sta_name,$sta_net,$sta_lat,$sta_lon,$sta_ele,$sta_bur)=split(" ",`egrep '$sta +$net' $opt_a`);
+ if (! -f $stafile ) {die(" No such files: $stafile");}
+ ($sta_name,$sta_net,$sta_lat,$sta_lon,$sta_ele,$sta_bur)=split(" ",`egrep '$sta +$net' $stafile`);
if (not defined $sta_name) {
- print("No such station $sta+$net in $opt_a file\n");
- ($sta_name,$sta_net,$sta_lat,$sta_lon,$sta_ele,$sta_bur)=split(" ",`egrep '$sta' $opt_a`);
+ print("No such station $sta+$net in $stafile file\n");
+ ($sta_name,$sta_net,$sta_lat,$sta_lon,$sta_ele,$sta_bur)=split(" ",`egrep '$sta' $stafile`);
if (not defined $sta_name) {
- print " No such station as $sta in the $opt_a file\n";
+ print " No such station as $sta in the $stafile file\n";
$sta_text .="$sta ";}}
print SAC "ch stla $sta_lat stlo $sta_lon stel $sta_ele stdp $sta_bur\nwh\n";}
- print SAC "w over\nquit\n";
- close(SAC);
- if ($opt_c) { open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
- print SAC "echo on\n r $outfile\n ";
-
if ($opt_l) { #cut record
print SAC "setbb begin ( max &1,b ( &1,o + $lmin ) ) \n";
print SAC "setbb end ( min &1,e ( &1,o + $lmax ) ) \n";
@@ -161,6 +162,11 @@
if ($opt_A) {print SAC "mul $opt_A\n";}
+ print SAC "w over\nquit\n";
+ close(SAC);
+ if ($opt_c) { open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
+ print SAC "echo on\n r $outfile\n ";
+
if ($opt_t) { # filter record
print "Filtering records...\n";
print SAC "rtrend\n rmean\n taper\n";
@@ -172,8 +178,10 @@
if ($opt_i) {# convolve with instrument response
print "Convolving instrument response...\n";
+ @cmpc=split(//,$comp);
+ $comp=$cmpd[0].$cmpd[1].$cmpc[2];
$pzfile="SAC_PZs_${net}_${sta}_${comp}_";
- @nfiles=`ls -l $opt_i/${pzfile}* | awk '{print \$9}'`;
+ @nfiles=`ls -l $ins_dir/${pzfile}* | awk '{print \$9}'`;
if (@nfiles != 1) {die("Pzfile is incorrect: \n at nfiles files\n");}
$pz=$nfiles[0]; chomp($pz);
if (! -f $pz) {die("Not a pz file $pz\n");}
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/rotate.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/rotate.pl 2011-02-24 03:35:46 UTC (rev 17967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/rotate.pl 2011-02-24 03:40:54 UTC (rev 17968)
@@ -17,10 +17,6 @@
For iris data (with -d), the timing part of the name will be ignored
and only the component part of the name will be changed correspondingly
- Make sure you have sac, saclst in the PATH before execution
-
- Qinya Liu, Caltech, May 2007
-
EOF
exit(1)
}
@@ -29,12 +25,11 @@
if (!getopts('l:L:cdts')) {die('Check input arguments\n');}
@ARGV > 0 or Usage();
if (!$opt_l) {$opt_l = 0;}
-$saclst = "/opt/seismo-util/bin/saclst";
-if (not -f $saclst) {die("No such file as $saclst\n");}
+$saclst = "saclst";
+if (system("which $saclst >/dev/null") != 0) {die(" No $saclst file\n");}
$undef=-12345.0;
$eps=0.1;
-
foreach $file (@ARGV) {
print "processing $file\n";
if (! -f $file) {die (" check to see if $file exists or not\n");}
@@ -73,16 +68,22 @@
if (abs($ninc-$undef)<$eps or abs($einc-$undef)<$eps or abs($naz-$undef)<$eps or abs($eaz-$undef)<$eps ) {die("Check header cmpinc and cmpaz for $north and $east\n");}
- if ($ny != $ey or $nj != $ej or $nh != $eh or $nm != $em or $ns != $es or $nms != $ems)
- {die(" Not same reference time for $north and $east\n");}
+ # allow n and e component .2 degrees from vertical
+ $daz = abs($naz-$eaz);
+ if ($daz > 180 and abs(abs($daz-270) >= 0.18) or $daz < 180 and abs(abs($daz-90) >=0.18) ) {
+ print "cmpaz not properly set for $north and $east: $naz,$eaz, skip\n";
+ next;}
+ if ($ny != $ey or $nj != $ej or $nh != $eh or $nm != $em or $ns != $es or abs($nms-$ems)>1)
+ {die(" Not the same reference time for $north and $east\n");}
+
if (abs($gcarc_e-$undef)<$eps or abs($gcarc_n-$undef)<$eps )
{die(" Check to see if GCARC is defined in the header\n");}
# check if the reference time is the same or not!
- if ($opt_c) {open(SAC,"|sac2000");}
- else {open(SAC,"|sac2000 > /dev/null");}
+ if ($opt_c) {open(SAC,"|sac");}
+ else {open(SAC,"|sac > /dev/null");}
print SAC "echo on\n";
print SAC "r $north $east\n ";
if ($opt_L or $npts_n != $npts_e or abs($nb-$eb) > $eps) { #cut properly
@@ -95,6 +96,7 @@
print SAC "cut %begin% %end%\n r $north $east\n cut off\n";
}
+ print SAC "ch cmpinc 90.0\n"; # make sure both are horizontal comp.
print SAC "rot\n"; # this is the default rotate to normal
print SAC "ch file 1 kcmpnm $rcomp\n";
print SAC "ch file 2 kcmpnm $tcomp\n";
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/sph2utm_Qinya_Liu.tar.gz
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90 2011-02-24 03:35:46 UTC (rev 17967)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90 2011-02-24 03:40:54 UTC (rev 17968)
@@ -741,6 +741,7 @@
! Mrt = -Mzn
! Mrp = Mze
! Mtp = -Mne
+ ! for consistency, location derivatives are in the order of [Xr,Xt,Xp]
! minus sign for sloc_der(3,irec_local) to get derivative for depth instead of radius
write(27,'(g16.5)') moment_der(3,3,irec_local) * 1e-7
More information about the CIG-COMMITS
mailing list