[cig-commits] r14894 - seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process
leif at geodynamics.org
leif at geodynamics.org
Wed May 6 12:10:14 PDT 2009
Author: leif
Date: 2009-05-06 12:10:14 -0700 (Wed, 06 May 2009)
New Revision: 14894
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl.in
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl.in
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl.in
Removed:
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/Makefile.am
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/configure.ac
Log:
Removed most hard-coded paths from Perl scripts.
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/Makefile.am
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/Makefile.am 2009-05-06 18:35:14 UTC (rev 14893)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/Makefile.am 2009-05-06 19:10:14 UTC (rev 14894)
@@ -1,5 +1,7 @@
## Process this file with automake to produce Makefile.in
+SUFFIXES = .pl.in .pl
+
bin_PROGRAMS = asc2sac convolve_stf
bin_SCRIPTS = process_data.pl process_syn.pl rotate.pl
@@ -8,3 +10,12 @@
asc2sac_SOURCES = asc2sac.c
convolve_stf_SOURCES = convolve_stf.c
+
+
+do_build = sed -e s:[@]sacaux[@]:$(SACAUX):g -e s:[@]prefix[@]:$(prefix):g
+
+.pl.in.pl:
+ $(do_build) $< > $@ || (rm -f $@ && exit 1)
+
+CLEANFILES = $(bin_SCRIPTS)
+EXTRA_DIST = process_data.pl.in process_syn.pl.in rotate.pl.in
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/configure.ac
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/configure.ac 2009-05-06 18:35:14 UTC (rev 14893)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/configure.ac 2009-05-06 19:10:14 UTC (rev 14894)
@@ -11,6 +11,9 @@
if test "x$SACAUX" = x; then
AC_MSG_ERROR([environment variable SACAUX is not set])
fi
+if test ! -d "$SACAUX"; then
+ AC_MSG_ERROR([SACAUX (=$SACAUX) does not name a directory])
+fi
# Checks for programs.
AC_PROG_CC
Deleted: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl 2009-05-06 18:35:14 UTC (rev 14893)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl 2009-05-06 19:10:14 UTC (rev 14894)
@@ -1,216 +0,0 @@
-#!/usr/bin/perl
-
-use Time::Local;
-use Getopt::Std;
-use POSIX;
-
-sub Usage{
- print STDERR <<END;
-
- process_data_new.pl -m CMTFILE -l Start/End -t Ts/Tl -i -p -x Ext -d OutDir
- -s Sps -P n/p -T Taper_width -c -y t1/v1/...
- data_sac_files
- with
-
- -m -- use CMTSOLUTION file to set event info
- -l Start/End -- start and end of record from o
- -t Ts/Tl -- shortest and longest period
- -i -- transfer record to displacement
- -p -- pick event p and s first arrival into sac headers (t1 and t2)
- -x Ext -- extension of new file name
- -d Dir -- directory to put the output files (default .)
- -s Sps -- sample per second (default 1.0)
- -P n/p -- number of poles and passes for butterworth filter (default 4/2)
- -T taper_width -- taper width to use in 'taper' command of SAC
- -c check sac output on the screen, otherwise print out summary
- -y t0/vel/t0/vel... --- pick times into headers t3/t4/... (rayleigh = 0/3.8; love = 0/4.0)
-
- data_sac_files --- names of the data sac files to be processed
-
- Notice:
- 0. Make sure you have sac, saclst, phtimes in PATH before execution
- 1. We require that polezero files are in the same directory as the sac
- data files which is generally satisfied. All the needed info are
- 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
-
- Qinya Liu, May 2007, Caltech
-
-END
- exit(1);
-}
-
-if (@ARGV == 0) { Usage(); }
-
-if (!getopts('m:l:t:ipx:d:s:P:T:y:c')) {die('Check input arguments\n');}
-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);}
-}
-if ($opt_x) {$ext='.'.$opt_x;}else {$ext="";}
-if (!$opt_s) {$dt=1.0;} else {$dt = 1.0/$opt_s;}
-if (!$opt_P) {$poles=4;$pass=2;}
-else{($poles,$pass)=split(/\//,$opt_P);
- if(not defined $pass or $pass < 1) {$pass=2;}
-}
-if (!$opt_T) {$opt_T = 0.05;}
-if ($opt_l) {($lmin,$lmax) = split(/\//,$opt_l);} else {$lmin = 0; $lmax = 3600;}
-
-$saclst="/opt/seismo-util/source/saclst/saclst";
-$phtimes="/opt/seismo-util/bin/phtimes";
-if (! -e $saclst) {die(" No $saclst file\n");}
-if (! -e $phtimes) {die("No $phtimes file\n");}
-
-$eps=1e-5; $undef=-12345.0; $cundef="-12345";
-
-foreach $file (@ARGV) {
- if (! -f $file) {die("No such file to be processed!!\n"); }
- print "Processing data file $file \n";
- if (not $opt_d) {$outfile = $file.$ext;}
- else { ($filename) = split(" ",`basename $file`);
- $outfile = "$out_dir/${filename}${ext}";}
- ($filedir) = split(" ",`dirname $file`);
- if ($ext or $opt_d) {system("\\cp -f $file $outfile");}
-
- # process data
- if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
- else {open(SAC,"|sac > /dev/null") || die("Can't open sac\n");}
- print SAC "echo on\n";
- print SAC "r $outfile\n";
-
- # add event information from CMTSOLUTION FILE
- if ($opt_m) {
- print " Adding event info from $opt_m\n";
- ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,undef,
- $elat,$elon,$edep,undef) = get_cmt($opt_m);
- ($oyear1,$ojday1,$ohr1,$omin1,$osec1,$omsec1)
- =tdiff($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift);
- print SAC "ch o gmt $oyear1 $ojday1 $ohr1 $omin1 $osec1 $omsec1\n";
- print SAC "evaluate to tmp 0 - &1,o\n";
- print SAC "ch allt %tmp% iztype io\n";
- print SAC "ch evla $elat evlo $elon evdp $edep \n";
- print SAC "w $outfile\nquit\n";
- close(SAC);
- if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
- else {open(SAC,"|sac > /dev/null") || die("Can't open sac\n");}
- print SAC "echo on\n";
- print SAC "r $outfile\n";}
-
- 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");}
- print SAC "setbb begin ( max &1,b ( &1,o + $lmin ) ) \n";
- print SAC "setbb end ( min &1,e ( &1,o + $lmax ) ) \n";
- print SAC "cut %begin% %end% \n";
- print SAC "r $outfile\n cut off\n w over \nquit\n";
- close (SAC);
- if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
- else {open(SAC,"|sac > /dev/null") || die("Can't open sac\n");}
- print SAC "echo on\n";
- print SAC "r $outfile\n";}
-
- if ($opt_s) {
- print SAC "interp delta $dt\n";}
-
- if ($opt_t) {# 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) { # transfer instrument response
- if (! $opt_t) {die(" Specify bandpass range by -t\n");}
- $f0=$f1*0.8;
- $f3=$f2*1.2;
- (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");}
- $pzfile=`ls -1 $filedir/SAC_PZs_${network}_${sta}_${comp}_${khole}* | head -1`;
- chomp($pzfile);
- if (!-f $pzfile) {die(" No pzfile $pzfile \n");}
- print " Transfer instrument response using pz file $pzfile\n";
- printf SAC ("trans from polezero s %20s to none freq %12.6f%12.6f%12.6f%12.6f \n",
- $pzfile,$f0,$f1,$f2,$f3);}
-
- 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`);
- if (abs($gcarc-$undef) < $eps or abs($edep-$undef) < $eps) {die("No gcarc and depth info\n");}
- ($Pph,$Ptime)=split(" ",`$phtimes $edep $gcarc P`);
- ($Sph,$Stime)=split(" ",`$phtimes $edep $gcarc S`);
- if (not $opt_m) {$tshift = 0;}
- print SAC "evaluate to tmp1 $Ptime - $tshift\n";
- print SAC "evaluate to tmp2 $Stime - $tshift\n";
- print SAC "ch t1 %tmp1% t2 %tmp2%\n";
- print SAC "ch kt1 $Pph kt2 $Sph\n";}
-
- if ($opt_y) { # add arrival times to headers
- @numbers=split(/\//,$opt_y); $npairs = floor(@numbers/2);
- if (@numbers != $npairs*2) {die("Enter -y t0/Vel/t0/Vel/...\n");}
- for ($i=0;$i<$npairs;$i++) {
- $int = $numbers[$i*2]; $slope = $numbers[$i*2+1];
- print " Add arrival time for waves with group velocity: $slope\n";
- (undef,$dist,$begin,$end)=split(" ",`$saclst dist b e f $outfile`);
- if (abs($dist-$undef) < $eps) {die("Not defined dist\n");}
- $h1 = 3+$i; $h1 = "t$h1"; $k1 = "k$h1"; $v1 = $int + $dist/$slope;
- (undef,$v1,undef) = sort {$a <=> $b} ($begin, $v1 ,$end);
- 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);
-}
-print " Done! \n";
-
-#**********************************************************
-sub mday2jday {
- my($oyear,$omonth,$oday)=@_;
- $omonth = $omonth-1; #months range from 0..11
- $time_sec = timegm(3,3,3,$oday,$omonth,$oyear);
- @t = gmtime($time_sec);
- my($jday) = $t[7];
- $jday += 1;
- return ($jday);
-}
-
-sub tdiff{
- my ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tadd)=@_;
- $time = timegm($osec, $omin, $ohr, $oday , $omonth-1, $oyear);
- $time += ($tadd +$omsec/1000); #event_time in machine format
- $msec = sprintf("%03.0f",($time - (floor($time)))*1000);
- $time = floor($time);
-#event-time:
- my($sec, $min, $hr, $day , $month, $year,$weekday,$jday) = gmtime($time);
- $month += 1;
- $year += 1900;
- $jday +=1;
- return ($year,$jday,$hr,$min,$sec,$msec);
-}
-
-sub get_cmt {
- my ($cmt_file)=@_;
- open(CMT, "$cmt_file") or die("Error opening $cmt_file\n");
- @cmt = <CMT>;
- my($pde,$oyear,$omonth,$oday,$ohr,$omin,$osec1)=split(" ",$cmt[0]);
- my($osec,$omsec)=split(/\./,$osec1); $omsec=$omsec*10;
- my(undef,undef,$tshift)=split(" ",$cmt[2]);
- my(undef,undef,$hdur)=split(" ",$cmt[3]);
- my(undef,$elat)=split(" ",$cmt[4]);
- my(undef,$elon)=split(" ",$cmt[5]);
- my(undef,$edep)=split(" ",$cmt[6]);
- my(undef,$Mrr)=split(" ",$cmt[7]);
- my(undef,$Mtt)=split(" ",$cmt[8]);
- my(undef,$Mpp)=split(" ",$cmt[9]);
- my(undef,$Mrt)=split(" ",$cmt[10]);
- my(undef,$Mrp)=split(" ",$cmt[11]);
- my(undef,$Mtp)=split(" ",$cmt[12]);
- close(CMT);
- return ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep,$Mrr,$Mtt,$Mpp,$Mrt,$Mrp,$Mtp);
-}
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl.in (from rev 14886, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl.in (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl.in 2009-05-06 19:10:14 UTC (rev 14894)
@@ -0,0 +1,219 @@
+#!/usr/bin/perl
+
+use Time::Local;
+use Getopt::Std;
+use POSIX;
+
+sub Usage{
+ print STDERR <<END;
+
+ process_data.pl -m CMTFILE -l Start/End -t Ts/Tl -i -p -x Ext -d OutDir
+ -s Sps -P n/p -T Taper_width -c -y t1/v1/...
+ data_sac_files
+ with
+
+ -m -- use CMTSOLUTION file to set event info
+ -l Start/End -- start and end of record from o
+ -t Ts/Tl -- shortest and longest period
+ -i -- transfer record to displacement
+ -p -- pick event p and s first arrival into sac headers (t1 and t2)
+ -x Ext -- extension of new file name
+ -d Dir -- directory to put the output files (default .)
+ -s Sps -- sample per second (default 1.0)
+ -P n/p -- number of poles and passes for butterworth filter (default 4/2)
+ -T taper_width -- taper width to use in 'taper' command of SAC
+ -c check sac output on the screen, otherwise print out summary
+ -y t0/vel/t0/vel... --- pick times into headers t3/t4/... (rayleigh = 0/3.8; love = 0/4.0)
+
+ data_sac_files --- names of the data sac files to be processed
+
+ Notice:
+ 0. Make sure you have sac, saclst, phtimes in PATH before execution
+ 1. We require that polezero files are in the same directory as the sac
+ data files which is generally satisfied. All the needed info are
+ 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
+
+ Qinya Liu, May 2007, Caltech
+
+END
+ exit(1);
+}
+
+$sacaux = $ENV{SACAUX} or "@sacaux@";
+unless (-d $sacaux) {die("missing SAC aux directory: '$sacaux'\n");}
+
+if (@ARGV == 0) { Usage(); }
+
+if (!getopts('m:l:t:ipx:d:s:P:T:y:c')) {die('Check input arguments\n');}
+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);}
+}
+if ($opt_x) {$ext='.'.$opt_x;}else {$ext="";}
+if (!$opt_s) {$dt=1.0;} else {$dt = 1.0/$opt_s;}
+if (!$opt_P) {$poles=4;$pass=2;}
+else{($poles,$pass)=split(/\//,$opt_P);
+ if(not defined $pass or $pass < 1) {$pass=2;}
+}
+if (!$opt_T) {$opt_T = 0.05;}
+if ($opt_l) {($lmin,$lmax) = split(/\//,$opt_l);} else {$lmin = 0; $lmax = 3600;}
+
+$saclst="$sacaux/../bin/saclst";
+$phtimes="phtimes";
+if (! -e $saclst) {die(" No $saclst file\n");}
+#if (! -e $phtimes) {die("No $phtimes file\n");}
+
+$eps=1e-5; $undef=-12345.0; $cundef="-12345";
+
+foreach $file (@ARGV) {
+ if (! -f $file) {die("No such file to be processed!!\n"); }
+ print "Processing data file $file \n";
+ if (not $opt_d) {$outfile = $file.$ext;}
+ else { ($filename) = split(" ",`basename $file`);
+ $outfile = "$out_dir/${filename}${ext}";}
+ ($filedir) = split(" ",`dirname $file`);
+ if ($ext or $opt_d) {system("\\cp -f $file $outfile");}
+
+ # process data
+ if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
+ else {open(SAC,"|sac > /dev/null") || die("Can't open sac\n");}
+ print SAC "echo on\n";
+ print SAC "r $outfile\n";
+
+ # add event information from CMTSOLUTION FILE
+ if ($opt_m) {
+ print " Adding event info from $opt_m\n";
+ ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,undef,
+ $elat,$elon,$edep,undef) = get_cmt($opt_m);
+ ($oyear1,$ojday1,$ohr1,$omin1,$osec1,$omsec1)
+ =tdiff($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift);
+ print SAC "ch o gmt $oyear1 $ojday1 $ohr1 $omin1 $osec1 $omsec1\n";
+ print SAC "evaluate to tmp 0 - &1,o\n";
+ print SAC "ch allt %tmp% iztype io\n";
+ print SAC "ch evla $elat evlo $elon evdp $edep \n";
+ print SAC "w $outfile\nquit\n";
+ close(SAC);
+ if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
+ else {open(SAC,"|sac > /dev/null") || die("Can't open sac\n");}
+ print SAC "echo on\n";
+ print SAC "r $outfile\n";}
+
+ 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");}
+ print SAC "setbb begin ( max &1,b ( &1,o + $lmin ) ) \n";
+ print SAC "setbb end ( min &1,e ( &1,o + $lmax ) ) \n";
+ print SAC "cut %begin% %end% \n";
+ print SAC "r $outfile\n cut off\n w over \nquit\n";
+ close (SAC);
+ if ($opt_c) {open(SAC,"|sac ") || die("Can't open sac\n");}
+ else {open(SAC,"|sac > /dev/null") || die("Can't open sac\n");}
+ print SAC "echo on\n";
+ print SAC "r $outfile\n";}
+
+ if ($opt_s) {
+ print SAC "interp delta $dt\n";}
+
+ if ($opt_t) {# 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) { # transfer instrument response
+ if (! $opt_t) {die(" Specify bandpass range by -t\n");}
+ $f0=$f1*0.8;
+ $f3=$f2*1.2;
+ (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");}
+ $pzfile=`ls -1 $filedir/SAC_PZs_${network}_${sta}_${comp}_${khole}* | head -1`;
+ chomp($pzfile);
+ if (!-f $pzfile) {die(" No pzfile $pzfile \n");}
+ print " Transfer instrument response using pz file $pzfile\n";
+ printf SAC ("trans from polezero s %20s to none freq %12.6f%12.6f%12.6f%12.6f \n",
+ $pzfile,$f0,$f1,$f2,$f3);}
+
+ 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`);
+ if (abs($gcarc-$undef) < $eps or abs($edep-$undef) < $eps) {die("No gcarc and depth info\n");}
+ ($Pph,$Ptime)=split(" ",`$phtimes $edep $gcarc P`);
+ ($Sph,$Stime)=split(" ",`$phtimes $edep $gcarc S`);
+ if (not $opt_m) {$tshift = 0;}
+ print SAC "evaluate to tmp1 $Ptime - $tshift\n";
+ print SAC "evaluate to tmp2 $Stime - $tshift\n";
+ print SAC "ch t1 %tmp1% t2 %tmp2%\n";
+ print SAC "ch kt1 $Pph kt2 $Sph\n";}
+
+ if ($opt_y) { # add arrival times to headers
+ @numbers=split(/\//,$opt_y); $npairs = floor(@numbers/2);
+ if (@numbers != $npairs*2) {die("Enter -y t0/Vel/t0/Vel/...\n");}
+ for ($i=0;$i<$npairs;$i++) {
+ $int = $numbers[$i*2]; $slope = $numbers[$i*2+1];
+ print " Add arrival time for waves with group velocity: $slope\n";
+ (undef,$dist,$begin,$end)=split(" ",`$saclst dist b e f $outfile`);
+ if (abs($dist-$undef) < $eps) {die("Not defined dist\n");}
+ $h1 = 3+$i; $h1 = "t$h1"; $k1 = "k$h1"; $v1 = $int + $dist/$slope;
+ (undef,$v1,undef) = sort {$a <=> $b} ($begin, $v1 ,$end);
+ 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);
+}
+print " Done! \n";
+
+#**********************************************************
+sub mday2jday {
+ my($oyear,$omonth,$oday)=@_;
+ $omonth = $omonth-1; #months range from 0..11
+ $time_sec = timegm(3,3,3,$oday,$omonth,$oyear);
+ @t = gmtime($time_sec);
+ my($jday) = $t[7];
+ $jday += 1;
+ return ($jday);
+}
+
+sub tdiff{
+ my ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tadd)=@_;
+ $time = timegm($osec, $omin, $ohr, $oday , $omonth-1, $oyear);
+ $time += ($tadd +$omsec/1000); #event_time in machine format
+ $msec = sprintf("%03.0f",($time - (floor($time)))*1000);
+ $time = floor($time);
+#event-time:
+ my($sec, $min, $hr, $day , $month, $year,$weekday,$jday) = gmtime($time);
+ $month += 1;
+ $year += 1900;
+ $jday +=1;
+ return ($year,$jday,$hr,$min,$sec,$msec);
+}
+
+sub get_cmt {
+ my ($cmt_file)=@_;
+ open(CMT, "$cmt_file") or die("Error opening $cmt_file\n");
+ @cmt = <CMT>;
+ my($pde,$oyear,$omonth,$oday,$ohr,$omin,$osec1)=split(" ",$cmt[0]);
+ my($osec,$omsec)=split(/\./,$osec1); $omsec=$omsec*10;
+ my(undef,undef,$tshift)=split(" ",$cmt[2]);
+ my(undef,undef,$hdur)=split(" ",$cmt[3]);
+ my(undef,$elat)=split(" ",$cmt[4]);
+ my(undef,$elon)=split(" ",$cmt[5]);
+ my(undef,$edep)=split(" ",$cmt[6]);
+ my(undef,$Mrr)=split(" ",$cmt[7]);
+ my(undef,$Mtt)=split(" ",$cmt[8]);
+ my(undef,$Mpp)=split(" ",$cmt[9]);
+ my(undef,$Mrt)=split(" ",$cmt[10]);
+ my(undef,$Mrp)=split(" ",$cmt[11]);
+ my(undef,$Mtp)=split(" ",$cmt[12]);
+ close(CMT);
+ return ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep,$Mrr,$Mtt,$Mpp,$Mrt,$Mrp,$Mtp);
+}
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_data.pl.in
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mergeinfo
+
Deleted: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl 2009-05-06 18:35:14 UTC (rev 14893)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl 2009-05-06 19:10:14 UTC (rev 14894)
@@ -1,250 +0,0 @@
-#!/usr/bin/perl
-
-use Time::Local;
-use Getopt::Std;
-use POSIX;
-
-sub Usage{
-print STDERR <<END;
-
-Usage: process_syn_new.pl
- -S -m CMTFILE -h -o offset -lStart/End -tTmin/Tmax -P n/p
- -i Dir -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
- -m CMTSOLUTION file for event origin, location and half duration
- -h convolve a triangle stf with half duration from -m
- -o moves the synthetics back by 'offset' from centroid time.
- -l lmin/lmax start/end cut of trace
- -t Tmin/Tmax -- shortest/longest period of bandpass filter
- -P -- number of poles and passes in butterworth filter(default 4/2)
- -i Dir -- convolve synthetics with instrument response 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"
- to get information from Vala''s station file
- -s sps -- resampling rate for synthetics (default is 1.0)
- -y t0/V0/t1/V1 -- add to t3/t4 the arrival times of t + dist/V
- -c -- show sac output
- -d OutDir -- directory to dump output (default: same as input files)
- -x Ext -- add extension Ext
-
- names of files -- name of syn files to be processed
-
- Make sure you have sac, saclst, phtimes, asc2sac, ascii2sac.csh
- in PATH before execution
-
- Qinya Liu, Caltech, May 2007
-
-END
-exit(1);
-}
-
- at ARGV > 1 or Usage();
-
-if (!getopts('Sm:ho:l:t:i:pP:a:cd:x:vs:y:A:')) {die(" check input arguments\n");}
-
-if ($opt_t) {($tmin, $tmax) = split(/\//,$opt_t);
- $f1 = 1./$tmax; $f2=1./$tmin;}
-if ($opt_x) {$ext=".".$opt_x;} else {$ext="";}
-if (!$opt_s) {$dt=1.0;} else {$dt = 1.0/$opt_s;}
-if (!$opt_P) {$poles=4;$pass=2;}
-else{($poles,$pass)=split(/\//,$opt_P);
- 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-util/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_A) {if ($opt_A !~ /^\d/) {die("-A option should be numbers\n");}}
-
-$undef = -12345;
-$eps = 0.1;
-
-$saclst="/opt/seismo-util/source/saclst/saclst";
-$phtimes="/opt/seismo-util/bin/phtimes";
-$asc2sac="/opt/seismo-util/bin/ascii2sac.csh";
-if (! -e $saclst) {die(" No $saclst file\n");}
-if (! -e $phtimes) {die("No $phtimes file\n");}
-if (! -e $asc2sac) {die("No $asc2sac file\n");}
-
-$min_hdur=1.0;
-
-if ($opt_m) { # extract event information from CMTSOLUTION
- if (!-f $opt_m) {die(" No cmtsolution file $opt_m\n");}
- ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep) = get_cmt($opt_m);
- ($oyear1,$ojday1,$ohr1,$omin1,$osec1,$omsec1) = tdiff($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift);}
-
-$sta_text="";
-
-foreach $file (@ARGV) {
-
- print "\nProcessing file $file\n";
- if (! -f $file) {die(" No such file : $file\n");}
-
- # transfer ascii file to sac file according to name of files
- if (not $opt_S) {# ascii file
- print "Transfer ascii file to sac file\n";
- system("$asc2sac $file >/dev/null");
- $file = $file.".sac";}
-
- (undef,$begin_time) = split(" ",`$saclst b f $file`);
- if (not defined $begin_time) {die("Check if the file is SAC format\n")};
-
- ($filename) = split(" ",`basename $file`);
- ($sta,$net,$comp)=split(/\./,$filename);
- if (not $opt_d) {$outfile = $file.$ext;}
- else {$outfile = "$opt_d/${filename}${ext}";}
- if ($ext or $opt_d) {system("\\cp -f $file $outfile");}
- print "Output to file $outfile\n";
-
- if ($opt_c) { open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
- print SAC "echo on\n";
-
- print SAC "r $outfile\n";
- if ($opt_m) { # add event location and original time
- print "Add event information\n";
- if($opt_o){$tmp_o = -$opt_o; print SAC "ch allt $tmp_o\n";}
- print SAC "ch nzyear $oyear1 nzjday $ojday1 nzhour $ohr1 nzmin $omin1 nzsec $osec1 nzmsec $omsec1\n";
- print SAC "ch o 0\n";
- print SAC "ch evla $elat evlo $elon evdp $edep \n wh\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_h and $hdur > $min_hdur) { # convolve source time function
- system("/opt/seismo-util/bin/convolve_stf g $hdur $outfile");
- system("mv ${outfile}.conv ${outfile}");}
-
-
- if ($opt_a) { # add station information
- print "Add station information\n";
- print SAC "ch kstnm $sta knetwk $net kcmpnm $comp\n ";
- 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 (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`);
- if (not defined $sta_name) {
- print " No such station as $sta in the $opt_a file\n";
- $sta_text .="$sta ";}}
- print SAC "ch stla $sta_lat stlo $sta_lon stel $sta_ele stdp $sta_bur\nwh\n";}
-
- if ($opt_l) {
- print SAC "setbb begin ( max &1,b ( &1,o + $lmin ) ) \n";
- print SAC "setbb end ( min &1,e ( &1,o + $lmax ) ) \n";
- print SAC "cut %begin% %end% \n";
- print SAC "r $outfile\n";
- print SAC "cut off\n";}
-
- if ($opt_A) {print SAC "mul $opt_A\n";}
-
- if ($opt_s) {print SAC "interp delta $dt\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";
- printf SAC ("bp n $poles p $pass co %10.5f %10.5f\n",$f1,$f2);
- print SAC "rtrend\n rmean\n taper\n";}
-
- if ($opt_i) {# convolve with instrument response
- print "Convolving instrument response...\n";
- $pzfile="SAC_PZs_${net}_${sta}_${comp}_";
- @nfiles=`ls -l $opt_i/${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");}
- print SAC "transfer from none to polezero s ${pz}\n";
- print SAC "w over\n";}
-
- if ($opt_p) { # pick arrival
- print "Picking arrivals...\n";
- (undef,$gcarc)=split(" ",`$saclst gcarc f $outfile`);
- if (abs($gcarc-$undef)< $eps) {die("No gcarc is defined\n");}
- ($Pph,$Ptime)=split(" ",`$phtimes $edep $gcarc P`);
- ($Sph,$Stime)=split(" ",`$phtimes $edep $gcarc S`);
- print SAC "evaluate to tmp1 $Ptime - $tshift\n";
- print SAC "evaluate to tmp2 $Stime - $tshift\n";
- print SAC "ch t1 %tmp1% t2 %tmp2%\n";
- print SAC "ch kt1 $Pph kt2 $Sph\n wh\n";}
-
- if ($opt_y) { # add arrival time for surface waves
- @numbers=split(/\//,$opt_y); $npairs = floor(@numbers/2);
- if (@numbers != $npairs*2) {die("Enter -y in1/slope1/in2/slope2/...\n");}
- for ($i=0;$i<$npairs;$i++) {
- $int = $numbers[$i*2]; $slope = $numbers[$i*2+1];
- print "Add arrival time for waves with group velocity: $slope\n";
- (undef,$dist,$begin,$end)=split(" ",`$saclst dist b e f $outfile`);
- if (abs($dist-$undef) < $eps) {die("Not defined dist\n");}
- $h1 = 3+$i; $h1 = "t$h1"; $k1 = "k$h1"; $v1 = $int + $dist/$slope;
- (undef,$v1,undef) = sort {$a <=> $b} ($begin, $v1 ,$end);
- printf SAC ("ch $h1 %12.2f\n",$v1);
- print SAC "ch $k1 $h1\n";}}
-
- print SAC "w over\nquit\n";
- close(SAC);
-
-}
-if ($sta_text) {
- print "\n\nNotice the following stations are not found in the station file\n $sta_text\n\n";}
-print "DONE\n";
-
-
-#****************************************************************
-sub mday2jday {
- my($oyear,$omonth,$oday)=@_;
- $omonth = $omonth-1; #months range from 0..11
- $time_sec = timegm(3,3,3,$oday,$omonth,$oyear);
- @t = gmtime($time_sec);
- $jday = @t[7];
- $jday += 1;
- return ($jday);
-}
-
-sub tdiff {
- my ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tadd)=@_;
- # die("$oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tadd\n");
- $time = timegm($osec, $omin, $ohr, $oday , $omonth-1, $oyear);
- $time += ($tadd +$omsec/1000); #event_time in machine format
- $msec = sprintf("%03.0f",($time - (floor($time)))*1000);
- $time = floor($time);
- #event-time:
- ($sec, $min, $hr, $day , $month, $year,$weekday,$jday) = gmtime($time);
- $month += 1;
- $year += 1900;
- $jday +=1;
- return ($year,$jday,$hr,$min,$sec,$msec);
-}
-
-sub get_cmt {
- local ($cmt_file)=@_;
- open(CMT, "$cmt_file") or die("Error opening $cmt_file\n");
- @cmt = <CMT>;
- local($pde,$oyear,$omonth,$oday,$ohr,$omin,$osec1)=split(" ",$cmt[0]);
- local($osec,$omsec)=split(/\./,$osec1); $omsec=$omsec*10;
- local(undef,undef,$tshift)=split(" ",$cmt[2]);
- local(undef,undef,$hdur)=split(" ",$cmt[3]);
- local(undef,$elat)=split(" ",$cmt[4]);
- local(undef,$elon)=split(" ",$cmt[5]);
- local(undef,$edep)=split(" ",$cmt[6]);
- local(undef,$Mrr)=split(" ",$cmt[7]);
- local(undef,$Mtt)=split(" ",$cmt[8]);
- local(undef,$Mpp)=split(" ",$cmt[9]);
- local(undef,$Mrt)=split(" ",$cmt[10]);
- local(undef,$Mrp)=split(" ",$cmt[11]);
- local(undef,$Mtp)=split(" ",$cmt[12]);
- close(CMT);
- return ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep,$Mrr,$Mtt,$Mpp,$Mrt,$Mrp,$Mtp);
-}
-
-
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl.in (from rev 14886, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl.in (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl.in 2009-05-06 19:10:14 UTC (rev 14894)
@@ -0,0 +1,253 @@
+#!/usr/bin/perl
+
+use Time::Local;
+use Getopt::Std;
+use POSIX;
+
+sub Usage{
+print STDERR <<END;
+
+Usage: process_syn_new.pl
+ -S -m CMTFILE -h -o offset -lStart/End -tTmin/Tmax -P n/p
+ -i Dir -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
+ -m CMTSOLUTION file for event origin, location and half duration
+ -h convolve a triangle stf with half duration from -m
+ -o moves the synthetics back by 'offset' from centroid time.
+ -l lmin/lmax start/end cut of trace
+ -t Tmin/Tmax -- shortest/longest period of bandpass filter
+ -P -- number of poles and passes in butterworth filter(default 4/2)
+ -i Dir -- convolve synthetics with instrument response 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"
+ to get information from Vala''s station file
+ -s sps -- resampling rate for synthetics (default is 1.0)
+ -y t0/V0/t1/V1 -- add to t3/t4 the arrival times of t + dist/V
+ -c -- show sac output
+ -d OutDir -- directory to dump output (default: same as input files)
+ -x Ext -- add extension Ext
+
+ names of files -- name of syn files to be processed
+
+ Make sure you have sac, saclst, phtimes, asc2sac, ascii2sac.csh
+ in PATH before execution
+
+ Qinya Liu, Caltech, May 2007
+
+END
+exit(1);
+}
+
+$sacaux = $ENV{SACAUX} or "@sacaux@";
+unless (-d $sacaux) {die("missing SAC aux directory: '$sacaux'\n");}
+
+ at ARGV > 1 or Usage();
+
+if (!getopts('Sm:ho:l:t:i:pP:a:cd:x:vs:y:A:')) {die(" check input arguments\n");}
+
+if ($opt_t) {($tmin, $tmax) = split(/\//,$opt_t);
+ $f1 = 1./$tmax; $f2=1./$tmin;}
+if ($opt_x) {$ext=".".$opt_x;} else {$ext="";}
+if (!$opt_s) {$dt=1.0;} else {$dt = 1.0/$opt_s;}
+if (!$opt_P) {$poles=4;$pass=2;}
+else{($poles,$pass)=split(/\//,$opt_P);
+ 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-util/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_A) {if ($opt_A !~ /^\d/) {die("-A option should be numbers\n");}}
+
+$undef = -12345;
+$eps = 0.1;
+
+$saclst="$sacaux/../bin/saclst";
+$phtimes="phtimes";
+$asc2sac="@prefix@/bin/ascii2sac.csh";
+if (! -e $saclst) {die(" No $saclst file\n");}
+#if (! -e $phtimes) {die("No $phtimes file\n");}
+if (! -e $asc2sac) {die("No $asc2sac file\n");}
+
+$min_hdur=1.0;
+
+if ($opt_m) { # extract event information from CMTSOLUTION
+ if (!-f $opt_m) {die(" No cmtsolution file $opt_m\n");}
+ ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep) = get_cmt($opt_m);
+ ($oyear1,$ojday1,$ohr1,$omin1,$osec1,$omsec1) = tdiff($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift);}
+
+$sta_text="";
+
+foreach $file (@ARGV) {
+
+ print "\nProcessing file $file\n";
+ if (! -f $file) {die(" No such file : $file\n");}
+
+ # transfer ascii file to sac file according to name of files
+ if (not $opt_S) {# ascii file
+ print "Transfer ascii file to sac file\n";
+ system("$asc2sac $file >/dev/null");
+ $file = $file.".sac";}
+
+ (undef,$begin_time) = split(" ",`$saclst b f $file`);
+ if (not defined $begin_time) {die("Check if the file is SAC format\n")};
+
+ ($filename) = split(" ",`basename $file`);
+ ($sta,$net,$comp)=split(/\./,$filename);
+ if (not $opt_d) {$outfile = $file.$ext;}
+ else {$outfile = "$opt_d/${filename}${ext}";}
+ if ($ext or $opt_d) {system("\\cp -f $file $outfile");}
+ print "Output to file $outfile\n";
+
+ if ($opt_c) { open(SAC,"|sac");}else {open(SAC,"|sac > /dev/null");}
+ print SAC "echo on\n";
+
+ print SAC "r $outfile\n";
+ if ($opt_m) { # add event location and original time
+ print "Add event information\n";
+ if($opt_o){$tmp_o = -$opt_o; print SAC "ch allt $tmp_o\n";}
+ print SAC "ch nzyear $oyear1 nzjday $ojday1 nzhour $ohr1 nzmin $omin1 nzsec $osec1 nzmsec $omsec1\n";
+ print SAC "ch o 0\n";
+ print SAC "ch evla $elat evlo $elon evdp $edep \n wh\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_h and $hdur > $min_hdur) { # convolve source time function
+ system("@prefix@/bin/convolve_stf g $hdur $outfile");
+ system("mv ${outfile}.conv ${outfile}");}
+
+
+ if ($opt_a) { # add station information
+ print "Add station information\n";
+ print SAC "ch kstnm $sta knetwk $net kcmpnm $comp\n ";
+ 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 (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`);
+ if (not defined $sta_name) {
+ print " No such station as $sta in the $opt_a file\n";
+ $sta_text .="$sta ";}}
+ print SAC "ch stla $sta_lat stlo $sta_lon stel $sta_ele stdp $sta_bur\nwh\n";}
+
+ if ($opt_l) {
+ print SAC "setbb begin ( max &1,b ( &1,o + $lmin ) ) \n";
+ print SAC "setbb end ( min &1,e ( &1,o + $lmax ) ) \n";
+ print SAC "cut %begin% %end% \n";
+ print SAC "r $outfile\n";
+ print SAC "cut off\n";}
+
+ if ($opt_A) {print SAC "mul $opt_A\n";}
+
+ if ($opt_s) {print SAC "interp delta $dt\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";
+ printf SAC ("bp n $poles p $pass co %10.5f %10.5f\n",$f1,$f2);
+ print SAC "rtrend\n rmean\n taper\n";}
+
+ if ($opt_i) {# convolve with instrument response
+ print "Convolving instrument response...\n";
+ $pzfile="SAC_PZs_${net}_${sta}_${comp}_";
+ @nfiles=`ls -l $opt_i/${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");}
+ print SAC "transfer from none to polezero s ${pz}\n";
+ print SAC "w over\n";}
+
+ if ($opt_p) { # pick arrival
+ print "Picking arrivals...\n";
+ (undef,$gcarc)=split(" ",`$saclst gcarc f $outfile`);
+ if (abs($gcarc-$undef)< $eps) {die("No gcarc is defined\n");}
+ ($Pph,$Ptime)=split(" ",`$phtimes $edep $gcarc P`);
+ ($Sph,$Stime)=split(" ",`$phtimes $edep $gcarc S`);
+ print SAC "evaluate to tmp1 $Ptime - $tshift\n";
+ print SAC "evaluate to tmp2 $Stime - $tshift\n";
+ print SAC "ch t1 %tmp1% t2 %tmp2%\n";
+ print SAC "ch kt1 $Pph kt2 $Sph\n wh\n";}
+
+ if ($opt_y) { # add arrival time for surface waves
+ @numbers=split(/\//,$opt_y); $npairs = floor(@numbers/2);
+ if (@numbers != $npairs*2) {die("Enter -y in1/slope1/in2/slope2/...\n");}
+ for ($i=0;$i<$npairs;$i++) {
+ $int = $numbers[$i*2]; $slope = $numbers[$i*2+1];
+ print "Add arrival time for waves with group velocity: $slope\n";
+ (undef,$dist,$begin,$end)=split(" ",`$saclst dist b e f $outfile`);
+ if (abs($dist-$undef) < $eps) {die("Not defined dist\n");}
+ $h1 = 3+$i; $h1 = "t$h1"; $k1 = "k$h1"; $v1 = $int + $dist/$slope;
+ (undef,$v1,undef) = sort {$a <=> $b} ($begin, $v1 ,$end);
+ printf SAC ("ch $h1 %12.2f\n",$v1);
+ print SAC "ch $k1 $h1\n";}}
+
+ print SAC "w over\nquit\n";
+ close(SAC);
+
+}
+if ($sta_text) {
+ print "\n\nNotice the following stations are not found in the station file\n $sta_text\n\n";}
+print "DONE\n";
+
+
+#****************************************************************
+sub mday2jday {
+ my($oyear,$omonth,$oday)=@_;
+ $omonth = $omonth-1; #months range from 0..11
+ $time_sec = timegm(3,3,3,$oday,$omonth,$oyear);
+ @t = gmtime($time_sec);
+ $jday = @t[7];
+ $jday += 1;
+ return ($jday);
+}
+
+sub tdiff {
+ my ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tadd)=@_;
+ # die("$oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tadd\n");
+ $time = timegm($osec, $omin, $ohr, $oday , $omonth-1, $oyear);
+ $time += ($tadd +$omsec/1000); #event_time in machine format
+ $msec = sprintf("%03.0f",($time - (floor($time)))*1000);
+ $time = floor($time);
+ #event-time:
+ ($sec, $min, $hr, $day , $month, $year,$weekday,$jday) = gmtime($time);
+ $month += 1;
+ $year += 1900;
+ $jday +=1;
+ return ($year,$jday,$hr,$min,$sec,$msec);
+}
+
+sub get_cmt {
+ local ($cmt_file)=@_;
+ open(CMT, "$cmt_file") or die("Error opening $cmt_file\n");
+ @cmt = <CMT>;
+ local($pde,$oyear,$omonth,$oday,$ohr,$omin,$osec1)=split(" ",$cmt[0]);
+ local($osec,$omsec)=split(/\./,$osec1); $omsec=$omsec*10;
+ local(undef,undef,$tshift)=split(" ",$cmt[2]);
+ local(undef,undef,$hdur)=split(" ",$cmt[3]);
+ local(undef,$elat)=split(" ",$cmt[4]);
+ local(undef,$elon)=split(" ",$cmt[5]);
+ local(undef,$edep)=split(" ",$cmt[6]);
+ local(undef,$Mrr)=split(" ",$cmt[7]);
+ local(undef,$Mtt)=split(" ",$cmt[8]);
+ local(undef,$Mpp)=split(" ",$cmt[9]);
+ local(undef,$Mrt)=split(" ",$cmt[10]);
+ local(undef,$Mrp)=split(" ",$cmt[11]);
+ local(undef,$Mtp)=split(" ",$cmt[12]);
+ close(CMT);
+ return ($oyear,$omonth,$oday,$ohr,$omin,$osec,$omsec,$tshift,$hdur,$elat,$elon,$edep,$Mrr,$Mtt,$Mpp,$Mrt,$Mrp,$Mtp);
+}
+
+
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/process_syn.pl.in
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mergeinfo
+
Deleted: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl 2009-05-06 18:35:14 UTC (rev 14893)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl 2009-05-06 19:10:14 UTC (rev 14894)
@@ -1,105 +0,0 @@
-#!/usr/bin/perl
-
-use Getopt::Std;
-use POSIX;
-
-sub Usage {
- print STDERR <<EOF;
-
- Usage: rotate.pl -lStart -LEnd -c -d [East|1]_Component_file_name
-
- rotate.pl rotates iris data(-d), trinet data and synthetics (with
- out -d option)
- -l -L specifies the start and end point of the traces to rotate
- -c check SAC output on screen
- ex. rotate.pl -d *.LHE.SAC for iris data
- rotate.pl PAS.*.LHE.sac for synthetics
- 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)
-}
-
-
-if (!getopts('l:L:cdts')) {die('Check input arguments\n');}
- at ARGV > 0 or Usage();
-if (!$opt_l) {$opt_l = 0;}
-$saclst = "saclst";
-if (not -f $saclst) {die("No such file as $saclst\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");}
- (undef,$comp)=split(" ",`$saclst kcmpnm f $file`);
- if ($comp eq "-12345") {die("No component name defined in the file\n");}
- if (not ($comp=~/E/ or $comp=~/1/)) {die("Please input only E/1 comp\n");}
- $ecomp = $comp; $ncomp = $ecomp; $rcomp = $ecomp; $tcomp = $ecomp;
- $ncomp=~s/E/N/; $ncomp =~s/1/2/;
- $tcomp=~s/E/T/; $tcomp =~s/1/T/;
- $rcomp=~s/E/R/; $rcomp =~s/1/R/;
- ($dir) = split(" ",`dirname $file`);$east = $file;
- if ($opt_d) {
- (undef,undef,undef,undef,undef,undef,$network,undef)=split(/\./,`basename $file`);
- (undef,$east1) = split(/\.$network\./,$file);
- $east1 = "$network.$east1";
- $north1 = $east1; $north1=~s/$ecomp/$ncomp/;
- ($north) = split(" ",`ls -1 $dir/*$north1`);
- $tang1=$east1;$radial1=$east1;
- $tang1=~s/$ecomp/$tcomp/; $radial1=~s/$ecomp/$rcomp/;
- $tang = "$dir/$tang1"; $radial = "$dir/$radial1";
- } else {
- $north = $east; $north =~s/$ecomp/$ncomp/;
- $tang = $east; $tang =~s/$ecomp/$tcomp/;
- $radial = $east; $radial =~s/$ecomp/$rcomp/;
-}
-
-# print " rotate $north $east \n to $tang $radial\n";
- if (!-f $north or !-f $east ) {
- print " no such files: $north and $east, check if -d option is missing\n";
- next;}
-
-# check some header variables
- @tmp=`$saclst b kcmpnm npts gcarc cmpinc cmpaz nzyear nzjday nzhour nzmin nzsec nzmsec f $north $east`;
- (undef,$nb,$ncomp,$npts_n,$gcarc_n,$ninc,$naz,$ny,$nj,$nh,$nm,$ns,$nms)=split(" ",$tmp[0]);
- (undef,$eb,$ecomp,$npts_e,$gcarc_e,$einc,$eaz,$ey,$ej,$eh,$em,$es,$ems)=split(" ",$tmp[1]);
-
- 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");}
-
- 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,"| 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
- if(!$opt_L){
- print SAC "setbb begin ( max &1,b &2,b ) \n";
- print SAC "setbb end ( min &1,e &2,e ) \n";}
- else{
- print SAC "setbb begin ( max &1,b &2,b ( &1,o + $opt_l ) ( &2,o + $opt_l ) ) \n";
- print SAC "setbb end ( min &1,e &2,e ( &1,o + $opt_L ) ( &2,o + $opt_L ) ) \n"; }
- print SAC "cut %begin% %end%\n r $north $east\n cut off\n";
-
- }
- 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";
- print SAC "w $radial $tang\nquit\n";
- close(SAC);
-}
-
-print " Done !\n";
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl.in (from rev 14886, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl.in (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl.in 2009-05-06 19:10:14 UTC (rev 14894)
@@ -0,0 +1,107 @@
+#!/usr/bin/perl
+
+use Getopt::Std;
+use POSIX;
+
+sub Usage {
+ print STDERR <<EOF;
+
+ Usage: rotate.pl -lStart -LEnd -c -d [East|1]_Component_file_name
+
+ rotate.pl rotates iris data(-d), trinet data and synthetics (with
+ out -d option)
+ -l -L specifies the start and end point of the traces to rotate
+ -c check SAC output on screen
+ ex. rotate.pl -d *.LHE.SAC for iris data
+ rotate.pl PAS.*.LHE.sac for synthetics
+ 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)
+}
+
+$sacaux = $ENV{SACAUX} or "@sacaux@";
+unless (-d $sacaux) {die("missing SAC aux directory: '$sacaux'\n");}
+
+if (!getopts('l:L:cdts')) {die('Check input arguments\n');}
+ at ARGV > 0 or Usage();
+if (!$opt_l) {$opt_l = 0;}
+$saclst = "$sacaux/../bin/saclst";
+if (not -f $saclst) {die("No such file as $saclst\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");}
+ (undef,$comp)=split(" ",`$saclst kcmpnm f $file`);
+ if ($comp eq "-12345") {die("No component name defined in the file\n");}
+ if (not ($comp=~/E/ or $comp=~/1/)) {die("Please input only E/1 comp\n");}
+ $ecomp = $comp; $ncomp = $ecomp; $rcomp = $ecomp; $tcomp = $ecomp;
+ $ncomp=~s/E/N/; $ncomp =~s/1/2/;
+ $tcomp=~s/E/T/; $tcomp =~s/1/T/;
+ $rcomp=~s/E/R/; $rcomp =~s/1/R/;
+ ($dir) = split(" ",`dirname $file`);$east = $file;
+ if ($opt_d) {
+ (undef,undef,undef,undef,undef,undef,$network,undef)=split(/\./,`basename $file`);
+ (undef,$east1) = split(/\.$network\./,$file);
+ $east1 = "$network.$east1";
+ $north1 = $east1; $north1=~s/$ecomp/$ncomp/;
+ ($north) = split(" ",`ls -1 $dir/*$north1`);
+ $tang1=$east1;$radial1=$east1;
+ $tang1=~s/$ecomp/$tcomp/; $radial1=~s/$ecomp/$rcomp/;
+ $tang = "$dir/$tang1"; $radial = "$dir/$radial1";
+ } else {
+ $north = $east; $north =~s/$ecomp/$ncomp/;
+ $tang = $east; $tang =~s/$ecomp/$tcomp/;
+ $radial = $east; $radial =~s/$ecomp/$rcomp/;
+}
+
+# print " rotate $north $east \n to $tang $radial\n";
+ if (!-f $north or !-f $east ) {
+ print " no such files: $north and $east, check if -d option is missing\n";
+ next;}
+
+# check some header variables
+ @tmp=`$saclst b kcmpnm npts gcarc cmpinc cmpaz nzyear nzjday nzhour nzmin nzsec nzmsec f $north $east`;
+ (undef,$nb,$ncomp,$npts_n,$gcarc_n,$ninc,$naz,$ny,$nj,$nh,$nm,$ns,$nms)=split(" ",$tmp[0]);
+ (undef,$eb,$ecomp,$npts_e,$gcarc_e,$einc,$eaz,$ey,$ej,$eh,$em,$es,$ems)=split(" ",$tmp[1]);
+
+ 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");}
+
+ 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,"|$sacaux/../bin/sac");}
+ else {open(SAC,"|$sacaux/../bin/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
+ if(!$opt_L){
+ print SAC "setbb begin ( max &1,b &2,b ) \n";
+ print SAC "setbb end ( min &1,e &2,e ) \n";}
+ else{
+ print SAC "setbb begin ( max &1,b &2,b ( &1,o + $opt_l ) ( &2,o + $opt_l ) ) \n";
+ print SAC "setbb end ( min &1,e &2,e ( &1,o + $opt_L ) ( &2,o + $opt_L ) ) \n"; }
+ print SAC "cut %begin% %end%\n r $north $east\n cut off\n";
+
+ }
+ 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";
+ print SAC "w $radial $tang\nquit\n";
+ close(SAC);
+}
+
+print " Done !\n";
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/UTILS/seis_process/rotate.pl.in
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mergeinfo
+
More information about the CIG-COMMITS
mailing list