[cig-commits] r14491 - in seismo/3D/ADJOINT_TOMO/measure_adj: . PLOTS scripts_meas scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Fri Mar 27 13:23:29 PDT 2009
Author: carltape
Date: 2009-03-27 13:23:28 -0700 (Fri, 27 Mar 2009)
New Revision: 14491
Modified:
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_station.pl
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats.pl
seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl
seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.m
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_3_adj_src_all.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl
seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
Log:
Updated run scripts based on socal tomography runs.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_station.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_station.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/make_pdf_by_station.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -10,7 +10,7 @@
# 1 distance
# 2 azimuth
#
-# EXAMPLE (from /net/sierra/raid1/carltape/results/MEASUREMENTS/m12/PDF_ALL_T006_T030):
+# EXAMPLE (from /net/sierra/raid1/carltape/results/MEASUREMENTS/m16/PDF_ALL_T006_T030):
# make_pdf_by_station.pl
#
#==========================================================
@@ -22,8 +22,8 @@
# USER INPUT
$emin = 1; # number of events recorded by a station required to make a composite PDF
-$smodel = "m12";
-$Ttag = "T006_T030";
+$smodel = "m16";
+$Ttag = "T003_T030";
$isort = 2; # sorted stations by distance (1) or azimuth (2)
# full stations file
@@ -35,6 +35,11 @@
# (See /ADJOINT_TOMO/iterate_adj/UTILS/station_lists/ )
$edir = "/net/sierra/raid1/carltape/results/SOURCES/EID_STATION_LISTS";
+# subset list of events that you want to use
+#$file_eid_sub = "/net/sierra/raid1/carltape/results/SOURCES/socal_16/EIDs_only_loc";
+$file_eid_sub = "/net/sierra/raid1/carltape/results/EID_LISTS/eids_simulation";
+if (not -f ${file_eid_sub}) {die("check if ${file_eid_sub} exist or not\n")}
+
# list of all possible event IDs
#$file_eid = "/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_${smodel}";
#if (not -f ${file_eid}) {die("check if ${file_eid} exist or not\n")}
@@ -59,7 +64,7 @@
#--------------------------------------
$imin = 1; $imax = $nrec;
-#$imin = 308; $imax = $imin;
+#$imin = 120; $imax = $imin;
# loop over all the stations in the sorted list
for ($ik = $imin; $ik <= $imax; $ik = $ik+1) {
@@ -71,7 +76,7 @@
($stalon,$stalat,$station,$network,undef,undef) = split(" ",$stalines[$ik-1]);
print "$ik out of $nrec : station $station.$network\n";
- # get list of events sorted by azimuth or distance
+ # get list of all possible events sorted by azimuth or distance
# EXAMPLE: EIDS_by_az_from_LT2.NP
$event_sort = "${edir}/EIDS_by_${sortlab}_from_${station}.${network}";
if (not -f "${event_sort}") {die("check if event_sort ${event_sort} exist or not\n")}
@@ -91,22 +96,27 @@
# get the event ID in the sorted list
($eid,$elon,$elat,undef,undef,undef,undef,undef) = split(" ",$elines[$j-1]);
- # EXAMPLE: 9828889_T006_T030_GSC_CI_m12_cc_win_adj.pdf
- @files = glob("${eid}_${Ttag}_${station}_${network}_${smodel}*pdf");
- $numf = @files;
- if ($numf == 0) {
- #print "--> no pdf file exists\n";
+ # check if event is in the subset list
+ ($nmatch,undef,undef) = split(" ",`grep $eid ${file_eid_sub} | wc`);
- } elsif ($numf == 1) {
- #print "--> pdf file exists\n";
- $pdffile = $files[0]; chomp($pdffile);
- $pdcat[$k] = $pdffile; $k = $k+1;
+ if ($nmatch == 1) {
+ # EXAMPLE: 9828889_T006_T030_GSC_CI_m16_cc_win_adj.pdf
+ @files = glob("${eid}_${Ttag}_${station}_${network}_${smodel}*pdf");
+ $numf = @files;
+ if ($numf == 0) {
+ #print "$j -- $eid --> no pdf file exists\n";
- } else {
- print "$eid\n";
- die("more than one pdf file exists\n");
+ } elsif ($numf == 1) {
+ #print "$j -- $eid --> pdf file exists\n";
+ $pdffile = $files[0]; chomp($pdffile);
+ $pdcat[$k] = $pdffile; $k = $k+1;
+
+ } else {
+ print "$eid\n";
+ die("more than one pdf file exists\n");
+ }
}
- }
+ } # for
# if there is at least one file, then make the composite PDF
if ($k > $emin+1) {
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -71,7 +71,8 @@
if($iker == 0) {$klab = "ik"; $ktitle = "waveform";}
elsif($iker == 1) {$klab = "mtm"; $ktitle = "multitaper TT";}
elsif($iker == 2) {$klab = "cc"; $ktitle = "cross-correlation TT";}
-else {die("ERROR: iker ($iker) must be 0,1, or 2\n");}
+elsif($iker == 3) {$klab = "bdcc"; $ktitle = "banana-doughnut CC-TT";}
+else {die("ERROR: iker ($iker) must be 0,1, 2, or 3\n");}
$sTmin = sprintf("T%3.3i",$Tmin);
$sTmax = sprintf("T%3.3i",$Tmax);
@@ -101,7 +102,8 @@
$iplot_CClabel = 1; # plot CC measurement labels for each window
$iplot_recon = 1; # plot reconstructed synthetics
if($iplot_CClabel==1) {$iplot_win = 1;}
-if($iplot_adj==1) {$swid = 3.5;} else {$swid = 7.4};
+if($iplot_adj==1) {$swid = 3.5; $fsizetitle = 10;}
+else {$swid = 7.4; $fsizetitle = 12;};
#---------------------------
@@ -545,9 +547,9 @@
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strZ\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# plot titles
-`echo \"$xtextpos3 $ytextpos5 12 0 $fontno BC Data (Black) and Synthetics (Red) -- ${Ttag_title}\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+`echo \"$xtextpos3 $ytextpos5 $fsizetitle 0 $fontno BC Data (Black) and Synthetics (Red) -- ${Ttag_title}\" | pstext $proj $bounds -N -K -O >> $ps_file`;
if ($iplot_adj==1) {
- `echo \"$xtextpos3 $ytextpos5 12 0 $fontno BC Adjoint Source ($ktitle)\" | pstext $proj $bounds -N -X$dX -K -O >> $ps_file`;
+ `echo \"$xtextpos3 $ytextpos5 $fsizetitle 0 $fontno BC Adjoint Source ($ktitle)\" | pstext $proj $bounds -N -X$dX -K -O >> $ps_file`;
}
#----------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -60,14 +60,14 @@
# plot southern California faults
$ifault = 1;
if($ifault==1) {
- $dir_fault = "/home/carltape/gmt";
- $fault_file = "${dir_fault}/faults/jennings.xy";
- $kcf_file = "${dir_fault}/faults/kcf.xy";
- $breck_file = "${dir_fault}/faults/breck.xy";
+ $dir_fault = "/home/carltape/gmt/faults";
+ $fault_file = "${dir_fault}/jennings_more.xy";
+ #$kcf_file = "${dir_fault}/kcf.xy";
+ #$breck_file = "${dir_fault}/breck.xy";
$fault_infoK = "-M -W0.5p,0/0/0";
print CSH "psxy $fault_file $R $J $fault_infoK -O -K -V >> $psfile\n";
- print CSH "psxy $kcf_file $R $J $fault_infoK -O -K -V >> $psfile\n";
- print CSH "psxy $breck_file $R $J $fault_infoK -O -K -V >> $psfile\n";
+ #print CSH "psxy $kcf_file $R $J $fault_infoK -O -K -V >> $psfile\n";
+ #print CSH "psxy $breck_file $R $J $fault_infoK -O -K -V >> $psfile\n";
}
# plot the stations
@@ -115,6 +115,7 @@
# create a new measurement file for plotting the histograms
$file0 = "${eid}_cc_dT_dA_dTsigma_dAsigma"; # see mt_measure_adj.f90
+# KEY: extract columns from window_chi file
#`awk '{print \$15,\$16,\$20,\$21,\$15/\$20,\$16/\$21}' $meas_file > $file0`;
`awk '{print \$15,\$16,\$19,\$20,\$15/\$19,\$16/\$20}' $meas_file > $file0`;
@@ -137,7 +138,8 @@
$max_sigma_DA = $max_sigma_DT; $min_sigma_DA = $min_sigma_DT ;
$max_DTerr = $max_DT; $min_DTerr = -$max_DTerr;
$max_DAerr = $max_DA; $min_DAerr = -$max_DAerr;
-if($Tmin == 2) {
+# change limits for shorter-period data
+if($Tmin == 2 || $Tmin == 3) {
$max_DT = 5; $min_DT = -$max_DT;
$max_DA = 1.5; $min_DA = -$max_DA;
$max_sigma_DT = 2.0; $min_sigma_DT = 0;
@@ -184,7 +186,7 @@
$Bhist5 = "-Ba4f1:\"dT / sigma-dT ($st5)\":/a${count_tick}:\"$sty\":WeSn";
$Bhist6 = "-Ba1f0.25:\"dlnA / sigma-dlnA ($st6)\":/a${count_tick}:\"$sty\":WeSn";
- $hist_info = "-G255 -W5 -G0/255/255 -L0.5p,0/0/0 -Z1";
+ $hist_info = "-W5 -G0/255/255 -L0.5p,0/0/0 -Z1";
# plot histograms
print CSH "pshistogram $file0 $hist_info $Jhist $Bhist1 $Rhist1 -T0 -W${DT_bin} -K -O -V $origin_hist1 >> $psfile\n";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -10,10 +10,8 @@
# The set of final histograms are copied to an output directory.
#
# EXAMPLE:
-# plot_win_stats_all.pl 6/30 m12 0
-# plot_win_stats_all.pl 6/30 m12 1 (plot histogram)
-# plot_win_stats_all.pl 2/30 m12 0
-# plot_win_stats_all.pl 2/30 m12 1 (plot histogram)
+# plot_win_stats_all.pl 6/30 m16 0
+# plot_win_stats_all.pl 6/30 m16 1 (plot histogram)
#
#-----------------------------------
@@ -27,11 +25,16 @@
$sTmax = sprintf("T%3.3i",$Tmax);
$Ttag = "${sTmin}_${sTmax}";
$sTper = sprintf("T = %i - %i s",$Tmin,$Tmax);
+$ftag = "${Ttag}_${smodel}";
-$dir_source = "/home/carltape/results/SOURCES/socal_12";
-$dir_source_text = "${dir_source}/v12_files_text";
-$file_eids = "${dir_source}/SOCAL_FINAL_CMT_v12_eid";
-#$file_eids = "/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_m12";
+$mid = 16; # index for directory with earthquake sources
+$dir_source = "/home/carltape/results/SOURCES/socal_${mid}";
+$dir_source_text = "${dir_source}/v${mid}_files_text";
+
+#$file_eids = "${dir_source}/SOCAL_FINAL_CMT_v${mid}_eid";
+#$file_eids = "/net/sierra/raid1/carltape/results/EID_LISTS/eids_simulation";
+$file_eids = "/net/sierra/raid1/carltape/results/EID_LISTS/eids_tomo";
+#$file_eids = "/net/sierra/raid1/carltape/results/EID_LISTS/eids_extra";
if (not -f $file_eids) {die("\n check if $file_eids exists\n")}
open(IN,$file_eids); @eids = <IN>; $nevent = @eids;
@@ -43,13 +46,13 @@
if (not -e $dir_run) {die("check if $dir_run exist or not\n")}
# name of file containing ALL source focal mechanisms in PSMECA format (see socal_tomo_sources.m)
-$cmtall_psmeca = "${dir_source}/SOCAL_FINAL_CMT_v11_psmeca_eid";
+$cmtall_psmeca = "${dir_source}/SOCAL_FINAL_CMT_v${mid}_psmeca_eid";
if (not -f $cmtall_psmeca) {die("\n check if $cmtall_psmeca exists\n")}
#$nevent = 5;
# make a directory to dump all the individual measurement files into
-$dir_meas_collect = "meas_files";
+$dir_meas_collect = "meas_files_${ftag}";
`rm -rf ${dir_meas_collect} ; mkdir ${dir_meas_collect}`;
# remove plot-related files from local directory
@@ -57,16 +60,24 @@
print "\n\n Running plot_win_stats_all.pl for $nevent possible events...\n";
+if (0==1) {
+ for ($i = 1; $i <= $nevent; $i = $i+1) {
+ $eid = $eids[$i-1]; chomp($eid);
+ print "-- $i, $eid --\n";
+ }
+ die("LISTING EIDS ONLY");
+}
+
# loop over all events
$k = 0;
$imin = 1; $imax = $nevent; # default
#$imin = 1; $imax = 20;
-#$imin = 204; $imax = $imin;
+#$imin = 81; $imax = $imin;
for ($i = $imin; $i <= $imax; $i = $i+1) {
$eid = $eids[$i-1]; chomp($eid);
- $ftag = "${eid}_${Ttag}";
+ $ftag = "${eid}_${Ttag}_${smodel}";
print "\n------------------------------";
print "\n $i, $eid \n";
@@ -178,7 +189,7 @@
$max_sigma_DA = $max_sigma_DT; $min_sigma_DA = $min_sigma_DT ;
$max_DTerr = $max_DT; $min_DTerr = -$max_DTerr;
$max_DAerr = $max_DA; $min_DAerr = -$max_DAerr;
-if($Tmin == 2) {
+if($Tmin == 2 || $Tmin == 3) {
$max_DT = 4; $min_DT = -$max_DT;
$max_DA = 1.5; $min_DA = -$max_DA;
$max_sigma_DT = 1.5; $min_sigma_DT = 0;
@@ -226,7 +237,7 @@
$Bhist5 = "-Ba4f1:\"dT / sigma-dT ($st5)\":/a${count_tick}:\"$sty\":WeSn";
$Bhist6 = "-Ba1f0.25:\"dlnA / sigma-dlnA ($st6)\":/a${count_tick}:\"$sty\":WeSn";
- $hist_info = "-G255 -W5 -G0/255/255 -L0.5p,0/0/0 -Z1";
+ $hist_info = "-W5 -G0/255/255 -L0.5p,0/0/0 -Z1";
print CSH "pshistogram $file0 $hist_info $Jhist $Bhist1 $Rhist1 -T0 -W${DT_bin} -K -V $origin_hist1 > $psfile\n"; # START
print CSH "pshistogram $file0 $hist_info $Jhist $Bhist2 $Rhist2 -T1 -W${DA_bin} -K -O -V $origin_hist2 >> $psfile\n";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_measure_adj.f90 2009-03-27 20:23:28 UTC (rev 14491)
@@ -27,10 +27,10 @@
double precision :: BEFORE_DLNA ! temporary: this should be an input parameter
double precision :: tshift_xc, sigma_dt_cc, dlnA, sigma_dlnA_cc, sigma_dt, sigma_dlnA
- double precision :: tr_chi, am_chi, tshift_f1f2, cc_max_f1f2, cc_max
+ double precision :: tr_chi, am_chi, tshift_f1f2, cc_max_f1f2, cc_max, T_pmax_dat, T_pmax_syn
double precision, dimension(NPT) :: dtau_w, dlnA_w, err_dt, err_dlnA, dzr0, dzr20
complex*16, dimension(NPT) :: trans_mtm
- integer :: nlen, nlen0, i_left, i_pmax, i_right, i_right0, istart,sta_len,net_len,chan_len, ipair, npairs, nwin, itmax
+ integer :: nlen, nlen0, i_left, i_pmax_dat, i_pmax_syn, i_right, i_right0, istart,sta_len,net_len,chan_len, ipair, npairs, nwin, itmax
logical :: use_trace, flag_short_window
!double precision :: trbdndw, a
!integer :: iord, passes
@@ -41,10 +41,10 @@
! input file MEASUREMENT.PAR
! is_mtm -- measurement type: 1 multi-taper; 2 cosine; 3 boxcar
! iker -- adjoint source type:
+ ! 0 waveform adjoint sources
! 1 multi-taper adjoint sources
! 2 xcorr adjoint sources
! 3 banana doughnut adjoint sources
- ! 0 waveform adjoint sources.
! TLONG, TSHORT -- starting and ending period of the signals
! tt, dtt, nn -- adj source interpolation scheme
! out_dir -- seismograms output directory
@@ -193,9 +193,11 @@
adj_file_prefix = trim(file_prefix2) // '.mtm'
else if (iker == 2) then
adj_file_prefix = trim(file_prefix2) // '.cc'
+ else if (iker == 3) then
+ adj_file_prefix = trim(file_prefix2) // '.bdcc'
else
print *, 'iker = ', iker
- stop 'iker must be 0, 1, or 2'
+ stop 'iker must be 0, 1, 2, or 3'
endif
print *
print *, trim(file_prefix2), ' --- '
@@ -232,8 +234,9 @@
window_chi(:) = 0.
! compute integrated waveform difference, normalized by duration of the record
- ! NOTE: (1) for comparison with waveform_chi, we include the 0.5 factor
- ! (2) we might want to include dt as an integration factor (also for waveform_chi),
+ ! NOTE: (1) this is for the FULL record, not the windowed record
+ ! (2) for comparison with waveform_chi, we include the 0.5 factor
+ ! (3) we might want to include dt as an integration factor (also for waveform_chi),
! but the ratio (d-s)^2 / d^2 avoids the need for dt, nstep, or length of record
window_chi(17) = 0.5 * sum( data**2 )
window_chi(18) = 0.5 * sum( syn**2 )
@@ -248,18 +251,22 @@
! CHT: also compute reconstructed synthetics for CC (and MT, if specified) measurements
call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,dzr20,dzr0,nlen,&
tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,tseis_recon_cc,&
- i_pmax,i_right0,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
+ i_pmax_dat,i_pmax_syn,i_right0,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
i_right = i_right0
i_left = 1
+ ! period of the max power of the synthetic record
+ T_pmax_dat = (dt*NPT) / dble(i_pmax_dat)
+ T_pmax_syn = (dt*NPT) / dble(i_pmax_syn)
+
! adjust measurements for adjoint source
if (iker == 1) then ! CHT: is_mtm --> iker
fstart = fstart0 ; fend = fend0
- call mt_measure_select(nlen,tshift_xc,i_pmax,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
+ call mt_measure_select(nlen,tshift_xc,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
print *, 'fstart0/fend0 :', fstart0, fend0
- print *, ' fstart/fend :', fstart, fend
- ! CHT: tshift_xc --> abs(tshift_xc) and tshift_f1f2 --> abs(tshift_f1f2)
+ print *, 'fstart/fend :', fstart, fend
+ print *, 'Tpmax :', T_pmax_dat, T_pmax_syn
if (use_trace .and. (cc_max < BEFORE_QUALITY .or. cc_max_f1f2 < AFTER_QUALITY .or. abs(tshift_xc) > BEFORE_TSHIFT .or. abs(tshift_f1f2) > AFTER_TSHIFT)) then
use_trace = .false.
if (DISPLAY_DETAILS) then
@@ -273,15 +280,15 @@
print *, ' Status: ', trim(measure_file_prefix), use_trace
! CHT: If MT measurement window is rejected by mt_measure_select, then use a CC measurement.
- ! --> But what about quality checks for the cross-correlation measurements?
+ ! --> May also want to consider quality checks for CC as well.
if(.not. use_trace) then
- !stop 'why was this MT measurement rejected?'
+ !stop 'Ceck why this MT measurement was rejected'
print *, 'Reverting from multitaper measurement to cross-correlation measurement'
iker = 2
is_mtm = 3
call mt_measure(is_mtm,iker,datafile,measure_file_prefix,data,syn,t0,dt,npts,tstart,tend,istart,dzr20,dzr0,nlen,&
tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,tseis_recon_cc,&
- i_pmax,i_right,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA)
+ i_pmax_dat,i_pmax_syn,i_right,trans_mtm,dtau_w,dlnA_w,sigma_dt,sigma_dlnA)
use_trace = .true.
endif
@@ -289,9 +296,9 @@
use_trace = .true.
endif
- ! CHT: If the DT is for some reason larger than the max allowable,
+ ! CHT: If the CC timeshift is for some reason larger than the allowable max,
! then effectively eliminate the window by zeroing the
- ! cross-correlation travel-time and amplitude measurements.
+ ! cross-correlation traveltime and amplitude measurements.
! See BEFORE_TSHIFT in subroutine compute_time_shift in mt_sub.f90.
! NOTE: BEFORE_DLNA should be included in the PARFILE
if ( abs(tshift_xc) > BEFORE_TSHIFT ) then
@@ -324,13 +331,15 @@
endif
! KEY: write misfit function values to file (two for each window)
+ ! Here are the 20 columns of the vector window_chi
! 1: MT-TT chi, 2: MT-dlnA chi, 3: XC-TT chi, 4: XC-dlnA chi
! 5: MT-TT meas, 6: MT-dlnA meas, 7: XC-TT meas, 8: XC-dlnA meas
! 9: MT-TT error, 10: MT-dlnA error, 11: XC-TT error, 12: XC-dlnA error
! WINDOW : 13: data power, 14: syn power, 15: (data-syn) power, 16: window duration
! FULL RECORD: 17: data power, 18: syn power, 19: (data-syn) power, 20: record duration
- write(13,'(a14,a8,a3,a5,i4,i4,2e14.6,20e14.6,2e14.6)') file_prefix0,sta,net,chan,j,iker,&
- tstart,tend,window_chi(:),tr_chi,am_chi
+ ! Example of a reduced file: awk '{print $2,$3,$4,$5,$6,$31,$32}' window_chi > window_chi_sub
+ write(13,'(a14,a8,a3,a5,i4,i4,2e14.6,20e14.6,2e14.6,2f14.6)') file_prefix0,sta,net,chan,j,iker,&
+ tstart,tend,window_chi(:),tr_chi,am_chi,T_pmax_dat,T_pmax_syn
print *, ' tr_chi = ', sngl(tr_chi), ' am_chi = ', sngl(am_chi)
! combine adjoint sources from different measurement windows
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/mt_sub.f90 2009-03-27 20:23:28 UTC (rev 14491)
@@ -41,7 +41,7 @@
subroutine mt_measure(is_mtm,iker,datafile,file_prefix,data,syn,t0,dt,npts,tstart,tend, &
istart,dzr20,dzr0,nlen,tshift_xc,sigma_dt_cc,dlnA,sigma_dlnA_cc,cc_max,tshift_f1f2,cc_max_f1f2,tseis_recon_cc, &
- i_pmax,i_right,trans_w,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
+ i_pmax_dat,i_pmax_syn,i_right,trans_w,dtau_w,dlnA_w,sigma_dt,sigma_dlnA,err_dt,err_dlnA)
integer, intent(in) :: is_mtm, iker, npts
double precision, dimension(:), intent(in) :: data, syn
@@ -51,7 +51,7 @@
double precision, intent(out) :: tshift_xc, sigma_dt_cc, dlnA, sigma_dlnA_cc, cc_max, tshift_f1f2,cc_max_f1f2, sigma_dt,sigma_dlnA
complex*16, dimension(:), intent(out) :: trans_w
double precision, dimension(:), intent(out) :: dtau_w, dlnA_w, dzr0, dzr20
- integer, intent(out) :: nlen, i_right, istart, i_pmax
+ integer, intent(out) :: nlen, i_right, istart, i_pmax_dat, i_pmax_syn
double precision, dimension(:), intent(out), optional :: err_dt, err_dlnA
!double precision, intent(out), optional :: sigma_dt,sigma_dlnA
@@ -82,7 +82,6 @@
!call dwrite_ascfile_f(trim(file_prefix)//'_syn.txt',t0,dt,npts,syn)
endif
-
!--------------------------------------------------------------------------
! window and interpolate data and synthetics
!--------------------------------------------------------------------------
@@ -178,7 +177,9 @@
! write cross-correlation measurement to file
call write_average_meas(file_prefix, IKER_CC, tshift_xc,dlnA,sigma_dt_cc,sigma_dlnA_cc)
- ! CHT: if you wave a simple waveform difference, then return
+ !========================================
+
+ ! CHT: if you want a simple waveform difference, then return
if (is_mtm == 0) return
!-----------------------------------------------------------------------------
@@ -214,6 +215,16 @@
call fft(LNPT,wseis_syn,FORWARD_FFT,dt)
call fft(LNPT,wseis_dat,FORWARD_FFT,dt)
+ ! index of the freq of the max power in the windowed data
+ ampmax_unw = 0.
+ i_pmax_dat = 1
+ do i = 1, fnum ! loop over frequencies
+ if( abs(wseis_dat(i)) > ampmax_unw) then
+ ampmax_unw = abs(wseis_dat(i))
+ i_pmax_dat = i
+ endif
+ enddo
+
! water level based untapered synthetics
! used to determine the i_right values (maximum frequency for measurement)
ampmax_unw = 0.
@@ -226,16 +237,16 @@
wtr_use_unw = cmplx(ampmax_unw * WTR, 0.)
! index of the freq of the max power in the windowed synthetics
- i_pmax = i_amp_max_unw
+ i_pmax_syn = i_amp_max_unw
i_right = fnum
i_right_stop = 0
do i = 1,fnum
- if (abs(wseis_syn(i))<=abs(wtr_use_unw) .and. i_right_stop==0 .and. i>i_amp_max_unw ) then
+ if( abs(wseis_syn(i)) <= abs(wtr_use_unw) .and. i_right_stop==0 .and. i > i_amp_max_unw ) then
i_right_stop = 1
i_right = i
endif
- if (abs(wseis_syn(i))>=10.*abs(wtr_use_unw) .and. i_right_stop==1 .and. i>i_amp_max_unw) then
+ if( abs(wseis_syn(i)) >= 10.*abs(wtr_use_unw) .and. i_right_stop==1 .and. i > i_amp_max_unw) then
i_right_stop = 0
i_right = i
endif
@@ -243,10 +254,15 @@
if (DISPLAY_DETAILS) then
print *, 'Frequency of max power in windowed synthetic (Hz):'
- print *, ' i_pmax = ', i_pmax, ', f_pmax = ', sngl(i_pmax * df), ', T_pmax = ', sngl(1./(i_pmax*df))
+ print *, ' i_pmax_syn = ', i_pmax_syn, ', f_pmax = ', sngl(i_pmax_syn * df), ', T_pmax = ', sngl(1./(i_pmax_syn*df))
print *, 'FFT freq spacing df = ', sngl(df)
print *, 'measurement spacing df_new = ', sngl(df_new)
print *, ' i_right = ', i_right, ', stopping freq = ', sngl(i_right * df)
+
+ ! write out power for each signal
+ call dwrite_ascfile_f(trim(file_prefix)//'.obs.power',df,df,i_right,abs(wseis_dat(1:i_right)) )
+ call dwrite_ascfile_f(trim(file_prefix)//'.syn.power',df,df,i_right,abs(wseis_syn(1:i_right)) )
+
endif
!-------------------------------------------------------------------------------
@@ -317,8 +333,8 @@
enddo
! for cosine or boxcar tapers only -- SEE COMMENTS BELOW for the multitaper case
- ! NOTE: here we are using trans_w, not trans_mtm
- ! BEWARE: The single-taper transfer function should give you a PERFECT FIT,
+ ! NOTE 1: here we are using trans_w, not trans_mtm
+ ! NOTE 2: The single-taper transfer function should give you a PERFECT FIT,
! but it is not relevant from the perspective of obtaining a measurement.
if (is_mtm /= 1) then
! phase, abs(trans), travel-time and amplitude as a func of freq for single-tapered measurements
@@ -601,6 +617,7 @@
! ----------------------
! frequency-domain tapers
+ ! THIS CHOICE CAN HAVE A SUBSTANTIAL EFFECT ON THE ADJOINT SOURCES
ipwr_w = 10
w_taper(:) = 0.
do i = i_left, i_right ! CHT: 1 --> i_left
@@ -629,11 +646,11 @@
wp_taper(i) = w_taper(i) / ffac
wq_taper(i) = w_taper(i) / ffac
- elseif (INCLUDE_ERROR == 1) then ! error estimate similar to CC error estimate
+ elseif (INCLUDE_ERROR == 1) then ! MT error estimate is assigned the CC error estimate
wp_taper(i) = w_taper(i) / ffac / (sigma_dt ** 2)
wq_taper(i) = w_taper(i) / ffac / (sigma_dlnA ** 2)
- elseif (INCLUDE_ERROR == 2) then ! jack-knife error estimate
+ elseif (INCLUDE_ERROR == 2) then ! MT jack-knife error estimate
err_t = err_dtau(i)
if (err_dtau(i) < dtau_wtr) err_t = err_t + dtau_wtr
err_A = err_dlnA(i)
@@ -814,7 +831,7 @@
tr_adj_src(i1) = -(tshift_xc / sigma_dt_cc**2 ) * ft_bar_t(i) * time_window(i)
am_adj_src(i1) = -(dlnA / sigma_dlnA_cc**2 ) * fa_bar_t(i) * time_window(i)
- ! banana-doughnut kernel adjoint source (no measurement) -- not yet an option!
+ ! banana-doughnut kernel adjoint source (no measurement)
elseif(iker==3) then
tr_adj_src(i1) = ft_bar_t(i) * time_window(i)
am_adj_src(i1) = fa_bar_t(i) * time_window(i)
@@ -867,7 +884,7 @@
tr_chi = window_chi(1)
am_chi = window_chi(2)
- elseif(iker==2) then ! cross_correlation
+ elseif(iker==2 .or. iker==3) then ! cross_correlation
tr_chi = window_chi(3)
am_chi = window_chi(4)
@@ -880,9 +897,9 @@
!-----------------------------------------------------------------------------
- subroutine mt_measure_select(nlen,tshift_xc,i_pmax,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
+ subroutine mt_measure_select(nlen,tshift_xc,i_pmax_syn,dtau_w,dlnA_w,err_dt,err_dlnA,dt,i_left,i_right,fstart,fend,use_trace)
- integer, intent(in) :: nlen, i_pmax
+ integer, intent(in) :: nlen, i_pmax_syn
double precision, intent(in) :: tshift_xc, dt
double precision, dimension(:), intent(inout) :: dtau_w, dlnA_w, err_dt, err_dlnA
double precision, intent(inout) :: fstart, fend
@@ -896,7 +913,7 @@
use_trace = .true.
df = 1./(dt*NPT)
- f_pmax = df * i_pmax
+ f_pmax = df * i_pmax_syn
T_pmax = 1./ f_pmax
Wlen = dt*nlen
@@ -914,7 +931,7 @@
! DECREASE the frequency range of the measurement (and adjoint source)
! --> note NCYCLE_IN_WINDOW and window length
- ! We arbitrarily state that we want at least 10 frequency points for the multitaper measurement.
+ ! We subjectively state that we want at least 10 frequency points for the multitaper measurement.
fstart = max(fstart, NCYCLE_IN_WINDOW/Wlen)
fend = min(fend, 1./(2.0*dt))
if (fstart >= fend - 10.0*df ) then
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.m 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.m 2009-03-27 20:23:28 UTC (rev 14491)
@@ -1,7 +1,7 @@
-function plot_mtm(data,syn,new,new_dt,dlnA,dt,dlnA_full,dt_full,dlnA_ave,dt_ave,...
+function plot_mtm(data,syn,dataf,synf,new,new_dt,dlnA,dt,dlnA_full,dt_full,dlnA_ave,dt_ave,...
new_cc,new_cc_dt,dlnA_cc,dt_cc,flims,iall)
-fsize = 11;
+fsize = 9;
lsize = 1.0;
msize = 12;
@@ -17,27 +17,35 @@
stdT_mtave = sprintf(' dT-mt-ave (blue) = %.3f +/- %.3f',dt_ave(1),dt_ave(2));
stdA_mtave = sprintf(' dlnA-mt-ave (blue) = %.3f +/- %.3f',dlnA_ave(1),dlnA_ave(2));
-subplot(2,1,1), hold on
+subplot(3,1,1), hold on
%title('Am','fontsize',fsize)
-plot(data(:,1),data(:,2),'-','linewidth',lsize,'Color','k')
-plot(syn(:,1),syn(:,2),'--','linewidth',lsize,'Color','r')
+plot(data(:,1),data(:,2),'-','linewidth',2*lsize,'Color','k')
+plot(syn(:,1),syn(:,2),'-','linewidth',2*lsize,'Color','r')
-% HOW CAN WE GET THE LEGEND TO PLOT WHERE WE WANT IT TO?
-% MY QUICK FIX WAS TO PUT IT OUTSIDE, BUT I HAD TO SHORTEN THE LABELS.
if 0==1
- plot(new_dt(:,1),new_dt(:,2),'--','linewidth',lsize,'Color','b')
- plot(new(:,1),new(:,2),'-','linewidth',lsize,'Color','b')
+ plot(new_dt(:,1),new_dt(:,2),'--','linewidth',lsize,'Color','cyan')
+ plot(new(:,1),new(:,2),'--','linewidth',lsize,'Color','cyan')
plot(new_cc_dt(:,1),new_cc_dt(:,2),'--','linewidth',lsize,'Color',cgreen)
- plot(new_cc(:,1),new_cc(:,2),'-','linewidth',lsize,'Color',cgreen)
- legend('Data','Syn','Syn-MT-1','Syn-MT-2','Syn-CC-1','Syn-CC-2','location','NorthWestOutside')
+ plot(new_cc(:,1),new_cc(:,2),'--','linewidth',lsize,'Color',cgreen)
+ legend('Data','Syn','Syn-MT-1','Syn-MT-2','Syn-CC-1','Syn-CC-2','location','BestOutside')
else
- plot(new(:,1),new(:,2),'-','linewidth',lsize,'Color','b')
- plot(new_cc(:,1),new_cc(:,2),'-','linewidth',lsize,'Color',cgreen)
- legend('Data','Syn','Syn-MT','Syn-CC','location','NorthWestOutside')
+ plot(new(:,1),new(:,2),'--','linewidth',lsize,'Color','cyan')
+ plot(new_cc(:,1),new_cc(:,2),'--','linewidth',lsize,'Color',cgreen)
+ legend('Data','Syn','Syn-MT','Syn-CC','location','BestOutside')
end
xlabel('Time [s]','fontsize',fsize)
-subplot(2,2,3), hold on
+% plot power of data and synthetics
+subplot(3,3,4), hold on
+plot(dataf(:,1),dataf(:,2),'-','linewidth',2*lsize,'Color','k')
+plot(synf(:,1),synf(:,2),'-','linewidth',2*lsize,'Color','r')
+xlim([f1 f2])
+ax2 = axis; [xmat, ymat] = vertlines(flims,ax2(3),ax2(4)); plot(xmat(:,1:4),ymat(:,1:4),'k',xmat(:,5:6),ymat(:,5:6),'r');
+ylabel('Power','fontsize',fsize)
+xlabel('Frequency [Hz]','fontsize',fsize)
+title('Power of data and synthetics','fontsize',fsize)
+
+subplot(3,3,5), hold on
% plot multitaper measurements
if iall == 1
errorbar(dlnA(:,1),dlnA(:,2),dlnA(:,3),'.','markersize',msize,'linewidth',lsize)
@@ -54,7 +62,7 @@
xlabel('Frequency [Hz]','fontsize',fsize)
title({'Amplitude measurement',stdA_cc,stdA_mtave},'fontsize',fsize)
-subplot(2,2,4), hold on
+subplot(3,3,6), hold on
% plot multitaper measurements
if iall==1
errorbar(dt(:,1),dt(:,2),dt(:,3),'.','markersize',msize,'linewidth',lsize)
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_meas/plot_mtm.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -10,7 +10,8 @@
print STDERR <<END;
plot_mtm.pl mtmfiles.dlnA
- Example (from OUTPUT_FILES): ../scripts/plot_mtm.pl *.dlnA
+ EXAMPLE (execute from OUTPUT_FILES):
+ ../scripts_meas/plot_mtm.pl *.dlnA
original coding by Vala, adapted by Qinya for the MT_MEASURE_ADJ package
@@ -29,7 +30,7 @@
# USER MUST CHANGE THIS PATH
#$progdir = "/home/lqy/SVN/mt_measure_adj/scripts";
#$progdir = "/opt/seismo-util/source/multitaper";
-$progdir = "/home/carltape/svn/specfem/mt_measure_adj_work/scripts";
+$progdir = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/measure_adj_work/scripts_meas";
my(@inputdirs);
@@ -55,6 +56,8 @@
print MATLAB "disp(['$filebase'])\n";
print MATLAB "data=dlmread('$inputdir/$filebase.obs');\n";
print MATLAB "syn=dlmread('$inputdir/$filebase.syn');\n";
+ print MATLAB "dataf=dlmread('$inputdir/$filebase.obs.power');\n";
+ print MATLAB "synf=dlmread('$inputdir/$filebase.syn.power');\n";
print MATLAB "new=dlmread('$inputdir/$filebase.recon_syn');\n";
print MATLAB "new_dt=dlmread('$inputdir/$filebase.recon_syn_dt');\n";
if ($iall == 1) {
@@ -89,7 +92,7 @@
print MATLAB "figure;\n";
# KEY COMMAND (plot_mtm.m)
- print MATLAB "plot_mtm(data,syn,new,new_dt,dlnA,dt,dlnA_full,dt_full,dlnA_ave,dt_ave,new_cc,new_cc_dt,dlnA_cc,dt_cc,flims,$iall);\n";
+ print MATLAB "plot_mtm(data,syn,dataf,synf,new,new_dt,dlnA,dt,dlnA_full,dt_full,dlnA_ave,dt_ave,new_cc,new_cc_dt,dlnA_cc,dt_cc,flims,$iall);\n";
print MATLAB "orient tall; suptitle(title);\n";
#print MATLAB "print ('-depsc2','$inputdir/${filebase}_measure.eps');\n";
print MATLAB "print ('-depsc2','$epsfile');\n";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_3_adj_src_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_3_adj_src_all.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/combine_3_adj_src_all.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -4,11 +4,13 @@
# Carl Tape, 03-Oct-2008
# combine_3_adj_src_all.pl
#
-# This script calls combine_adj_src.pl to plot make composite adjoint sources.
+# This script calls combine_adj_src.pl to plot make composite adjoint sources
# for a set of events.
#
# EXAMPLE:
-# combine_3_adj_src_all.pl 2/30 3/30 6/30 cc cc cc m13
+# combine_3_adj_src_all.pl 2/30 3/30 6/30 mtm mtm mtm m15
+# combine_3_adj_src_all.pl 2/30 3/30 6/30 cc cc cc m14
+# combine_3_adj_src_all.pl 2/30 3/30 6/30 bdcc bdcc bdcc m16
#
#-----------------------------------
@@ -26,7 +28,6 @@
# list of event IDs
$file_eids = "/net/sierra/raid1/carltape/results/EID_LISTS/kernels_run_${smodel}";
-#$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/socal_11/EIDs_only_eid";
if (not -f $file_eids) {die("\n check if $file_eids exists\n")}
open(IN,$file_eids); @eids = <IN>; $nevent = @eids;
@@ -54,8 +55,8 @@
# loop over all events
$k = 0;
$imin = 1; $imax = $nevent; # default
-#$imin = 1; $imax = 45;
-#$imin = 54; $imax = $imin;
+#$imin = 1; $imax = 128;
+#$imin = 130; $imax = $imin;
for ($i = $imin; $i <= $imax; $i = $i+1) {
@@ -93,6 +94,9 @@
} else {
print " --> all three STATIONS_ADJOINT files do NOT exist\n";
+ if (not -f $sta_file1) {print "sta_file1 --> $sta_file1\n";}
+ if (not -f $sta_file2) {print "sta_file2 --> $sta_file2\n";}
+ if (not -f $sta_file3) {print "sta_file3 --> $sta_file3\n";}
#die("STOPPING: all three STATIONS_ADJOINT files do NOT exist\n");
}
}
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/prepare_mt_measure_adj.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -53,10 +53,11 @@
# directories for data, synthetics, and window files : MUST BE CHANGED FOR EACH USER
if ($iwindow==1) {
# windowing code run directory
- $dir_win_run = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/flexwin_run";
+ $dir_win_run = "/net/denali/raid1/carltape/svn/cig/seismo/3D/flexwin_run";
$dir_win_run_syn = "${dir_win_run}/${dir_syn}";
$dir_win_run_data = "${dir_win_run}/${dir_data}";
$win_in = "${dir_win_run}/${dir_meas}/MEASUREMENT_WINDOWS";
+ $par_in = "${dir_win_run}/${dir_meas}/PAR_FILE";
#$dir_win_run_meas = "${dir_win_run}/${dir_meas}";
} else {
@@ -66,6 +67,7 @@
$dir_win_run_syn = "$dir1/SYN";
$dir_win_run_data = "$dir1/DATA";
$win_in = "$dir1/MEASUREMENT_WINDOWS_${eid}_${Ttag}_${smodel}";
+ $par_in = "$dir1/PAR_FILE";
#$dir0 = "/net/sierra/raid1/carltape/socal/socal_3D";
#$dir_win_run_syn = "$dir0/SYN/model_m0/${eid}/PROCESSED";
@@ -76,6 +78,7 @@
if (not -e ${dir_win_run_syn}) {die("check if ${dir_win_run_syn} exist or not\n");}
if (not -e ${dir_win_run_data}) {die("check if ${dir_win_run_data} exist or not\n");}
if (not -e ${win_in}) {die("check if ${win_in} exist or not\n");}
+if (not -e ${par_in}) {die("check if ${par_in} exist or not\n");}
# directories for measurement code
$dir_meas_syn = ${dir_syn};
@@ -89,8 +92,8 @@
open(CSH,">$cshfile");
# remove folders in the measurement directory
-print CSH "\\rm -rf ${dir_meas_syn} ${dir_meas_data} \n";
-#print CSH "mkdir ${dir_meas_syn} ${dir_meas_data} \n";
+print CSH "\\rm -rf ${dir_meas_syn} ${dir_meas_data}\n";
+#print CSH "mkdir ${dir_meas_syn} ${dir_meas_data}\n";
# copy ALL data and synthetics into measurement directory (BHZ,BHR,BHT)
# NOTE: copy the whole directory, because the list might be too long for using cp
@@ -102,8 +105,7 @@
#print CSH "cp ${dir_win_run_data}/* ${dir_meas_data}\n";
# CMTSOLUTION file
-#$dir_cmt = "/net/sierra/raid1/carltape/results/SOURCES/socal_6/CMT_files_post_inverted";
-$dir_cmt = "/net/sierra/raid1/carltape/results/SOURCES/socal_11/CMT_files_pre_inverted";
+$dir_cmt = "/net/sierra/raid1/carltape/results/SOURCES/socal_16/v16_files";
$cmtfile = "${dir_cmt}/CMTSOLUTION_${eid}";
if (not -f $cmtfile) {die("check if cmtfile $cmtfile exist or not\n")}
print CSH "\\rm CMTSOLUTION*\n";
@@ -158,9 +160,23 @@
# print CSH "\\cp @files ${win_out}\n";
#}
-# copy windowing file, which must have the prefix MEASUREMENT
+# check that the number of records matches the first line of the file
+# NOTE: this is important when manually removing windows
+($nseis,undef,undef) = split(" ",`grep DATA ${win_in} | wc`); # number of listed data records
+open(IN,${win_in}); @lines = <IN>; close(IN); $nlist = $lines[0]; chomp($nlist); # header number
+print "$eid -- nlist = $nlist, nseis = $nseis\n";
+if($nlist != $nseis) {
+ print "${win_in}\n";
+ die("check window files\n");
+}
+
+# copy windowing file and PAR file
$win_out = "./MEASUREMENT.WINDOWS";
+$par_out = "./PAR_FILE";
+print CSH "\\rm ${win_out}\n";
print CSH "\\cp ${win_in} ${win_out}\n";
+print CSH "\\rm ${par_out}\n";
+print CSH "\\cp ${par_in} ${par_out}\n";
# copy the windows file to PLOTS
if ($iplot == 1) {print CSH "cp $win_out PLOTS\n"}
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_cc_plot.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -111,13 +111,14 @@
print CSH "\\mv STATIONS_ADJOINT PLOTS\n";
print CSH "cd PLOTS\n";
-print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b 1 -k 1 -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS -i $smodel\n";
+print CSH "plot_win_adj_all.pl -l $lcut -m ../$cmtfile -b 1 -k 1 -a $stafile2 -d $dir_data -s $dir_syn -c $dir_recon -w MEASUREMENT.WINDOWS -i $smodel -j $Ts\n";
+
print CSH "\\rm test*.pdf\n";
# make a single PDF file
-$isort = 3;
-$otag = "mt_cc_all";
-print CSH "make_pdf_by_event.pl $isort $otag\n";
+#$isort = 3;
+#$otag = "mt_cc_all";
+#print CSH "make_pdf_by_event.pl $isort $otag\n";
#print CSH "/home/carltape/bin/pdcat -r *.pdf mt_cc_all.pdf\n"; # replace with make_pdf_by_event.pl
print CSH "cd -\n";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_mt_measure_adj.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -73,6 +73,35 @@
# empty the output files directory and compile the measurement code
print CSH "\\rm -rf OUTPUT_FILES ; mkdir OUTPUT_FILES ; make clean\n make\n";
+#----------------------------
+# choose ipar2 based on the values used in the windowing code PAR_FILE (Feb 2009)
+# This is a more sensible option, since you guarantee consistency.
+if(1==1) {
+ # Allow for some "loosening" of the parameters, since the exact nature
+ # of the measurement in the measurement code may exclude some input windows.
+ $dt_fac = 0.5; # in seconds
+ $cc_fac = 0.02; # from 0.0 to 1.0
+
+ (undef,undef,$tshift_base) = split(" ",`grep TSHIFT_BASE PAR_FILE`);
+ (undef,undef,$tshift_ref) = split(" ",`grep TSHIFT_REFERENCE PAR_FILE`);
+ (undef,undef,$dlna_base) = split(" ",`grep DLNA_BASE PAR_FILE`);
+ (undef,undef,$dlna_ref) = split(" ",`grep DLNA_REFERENCE PAR_FILE`);
+ (undef,undef,$cc_base) = split(" ",`grep CC_BASE PAR_FILE`);
+ print "Values from PAR_FILE:\n";
+ print "-- $tshift_base -- $tshift_ref -- $dlna_base -- $dlna_ref -- $cc_base --\n";
+ $BEFORE_QUALITY = $cc_base - $cc_fac;
+ $AFTER_QUALITY = $cc_base - $cc_fac;
+ $BEFORE_TSHIFT = abs($tshift_ref) + $tshift_base + $dt_fac;
+ $BEFORE_DLNA = abs($dlna_ref) + $dlna_base;
+
+ $par2_old = $par2;
+ ($BEFORE_QUALITY_OLD,$AFTER_QUALITY_OLD,$BEFORE_TSHIFT_OLD,$AFTER_TSHIFT_OLD) = split("/",$par2);
+ $AFTER_TSHIFT = $AFTER_TSHIFT_OLD;
+ $par2 = "${BEFORE_QUALITY}/${AFTER_QUALITY}/${BEFORE_TSHIFT}/${AFTER_TSHIFT}";
+ print "Updating par2:\n old -- $par2_old --\n new -- $par2 --\n";
+}
+#----------------------------
+
# create MEASUREMENT.PAR file
$wtr = 0.02; $npi = 2.5; # "default" values
#$itaper = 3;
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-03-27 20:05:55 UTC (rev 14490)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2009-03-27 20:23:28 UTC (rev 14491)
@@ -4,35 +4,35 @@
#
# run_tomo.pl
# Carl Tape
-# 19-Oct-2008
+# 27-March-2009
#
# This is the master program for running the measurement code on all the
# events for which the windowing code has been run.
#
# EXAMPLE (from working measure_adj directory) :
-# run_tomo.pl m14 0 6/30 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/5.0/1.0 1.0/0.5 1.0
-# run_tomo.pl m14 0 3/30 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/3.0/1.0 1.0/0.5 1.0
-# run_tomo.pl m14 0 2/30 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/2.0/1.0 1.0/0.5 1.0
+# run_tomo.pl m16 0 6/30 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/5.0/1.0 1.0/0.5 1.0
+# run_tomo.pl m16 0 3/30 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/3.0/1.0 1.0/0.5 1.0
+# run_tomo.pl m16 0 2/30 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/2.0/1.0 1.0/0.5 1.0
#
# socal m00 - m04
# run_tomo.pl m04 0 6/40 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.2 0.8
# run_tomo.pl m04 0 2/40 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/8.0/1.0 0.2/0.2 0.5
#
-# New m00 run:
-# run_tomo.pl m00 0 6/30 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/12.0/1.0 1.0/0.5 1.0
-# run_tomo.pl m00 0 3/30 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/6.0/1.0 1.0/0.5 1.0
-# run_tomo.pl m00 0 2/30 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/4.0/1.0 1.0/0.5 1.0
+# Multitaper (m00, m01, m02, m03, m15, m16):
+# run_tomo.pl m16 0 6/30 0.011 1/1/1 2.0/2.5/3.5/1.5 0.7/0.7/5.0/1.0 1.0/0.5 1.0 # 0.8
+# run_tomo.pl m16 0 3/30 0.011 1/1/1 0.25/0.5/3.5/1.5 0.7/0.7/3.0/1.0 1.0/0.5 1.0 # 0.6
+# run_tomo.pl m16 0 2/30 0.011 1/1/1 0.25/0.5/3.5/1.5 0.7/0.7/2.0/1.0 1.0/0.5 1.0 # 0.4
#
# socal m00 - m04, 0.2/0.2
# socal m05 - m08, 0.4/0.4
# socal m09, MAX: 4.0 pm 0.6 (2-30s), 4.0 pm 0.8 (3-30s), 5.0 pm 0.8 (6-30s)
# socal m10, MAX: 2.0 pm 0.6 (2-30s), 3.0 pm 0.8 (3-30s), 5.0 pm 0.8 (6-30s)
# socal m12, MAX: 2.0 pm 0.6 (2-30s), 3.0 pm 0.8 (3-30s), 5.0 pm 1.0 (6-30s)
-# socal m13 --> , MAX: 2.0 pm 1.0 (2-30s), 3.0 pm 1.0 (3-30s), 5.0 pm 1.0 (6-30s)
+# socal m13 --> m16, MAX: 2.0 pm 1.0 (2-30s), 3.0 pm 1.0 (3-30s), 5.0 pm 1.0 (6-30s)
#
#==========================================================
-if (@ARGV < 8) {die("Usage: run_tomo.pl smodel Tmin/Tmax dt itest iparbools par1 par2 par3\n")}
+if (@ARGV < 8) {die("Usage: run_tomo.pl smodel ibandpass Tmin/Tmax dt iparbools par1 par2 par3 lcut_frac\n")}
($smodel,$ibp,$Ts,$dt,$iparbools,$par1,$par2,$par3,$lcut_frac) = @ARGV;
($Tmin,$Tmax) = split("/",$Ts);
@@ -41,7 +41,7 @@
$Ttag = "${sTmin}_${sTmax}";
# cut times for plotting purposes only
-($lcut_min,$lcut_max) = split("/",$lcut);
+#($lcut_min,$lcut_max) = split("/",$lcut);
#================================================================
# USER INPUT
@@ -55,14 +55,14 @@
if (not -e $dir_pdf_all) {die("check if dir_pdf_all $dir_pdf_all exist or not\n")}
# file containing event IDs, half-durations for source, and number of time-steps for simulations
-$file_eid = "$dir0/results/SOURCES/socal_12/EIDs_simulation_duration";
+$file_eid = "$dir0/results/SOURCES/socal_16/EIDs_simulation_duration";
if (not -f ${file_eid}) {die("check if ${file_eid} exist or not\n")}
open(IN,"${file_eid}"); @elines = <IN>; close(IN);
$nevent = @elines;
print "\n $nevent events in the list \n";
# full stations file (copy into PLOTS directory)
-$file_stations = "/net/denali/home1/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
+$file_stations = "/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
if (not -f ${file_stations}) {die("check if file_stations ${file_stations} exist or not\n")}
`cp $file_stations PLOTS/STATIONS_TOMO`;
@@ -79,13 +79,14 @@
#die("testing");
$imin = 1; $imax = $nevent;
-#$imin = 1; $imax = 91;
-#$imin = 133; $imax = $imin;
+#$imin = 100; $imax = $nevent;
+#$imin = 2; $imax = $imin;
# EXAMPLE lines from EIDs_simulation_duration
-# 5 10222753 300.0 27300 0.26 332.0 4.010
-# 6 9700217 300.0 27300 0.16 331.7 3.600
-# 7 10222697 300.0 27300 0.21 331.0 3.800
+# 2 9967025 200.000 18200 0.29 4.084
+# 3 9967137 200.000 18200 0.33 4.211
+# 4 9967249 200.000 18200 0.32 4.191
+# 5 9967541 200.000 18200 0.20 3.790
#=============================================
# write the C-shell script to file
@@ -96,7 +97,7 @@
for ($ievent = $imin; $ievent <= $imax; $ievent++) {
# read in a row of numbers
- ($eindex,$eid,$Ldur_round,$nstep,$hdur,$Ldur,$mw) = split(" ",$elines[$ievent-1]);
+ ($eindex,$eid,$Ldur_round,$nstep,$hdur,$mw) = split(" ",$elines[$ievent-1]);
# get sorted stations file for ordering the plots
$file_stations_sort = "$dir0/results/SOURCES/EID_STATION_LISTS/STATIONS_by_${sortlab}_from_${eid}";
@@ -136,7 +137,7 @@
} else {
if (not -f ${window_file}) {
- print "--> DO NOT RUN: no window file exists ($window_file)\n";
+ print "--> DO NOT RUN: no window file exists\n${window_file}\n";
} else {
@@ -160,9 +161,19 @@
# print CSH "run_mt_cc_plot.pl $smodel $ibp $Ts $lcut $tvec 0 $iparbools $par1 $par2 $par3\n";
# }
- # SOCAL: as of model m04, we are no longer making MT measurements
- print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+# # SOCAL: as of model m04, we are no longer making MT measurements
+# print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+
+ # SOCAL: CC only
+ #print CSH "run_mt_measure_adj.pl $smodel 2 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+
+ # SOCAL: multitaper is back for m15 and m16
+ #print CSH "run_mt_measure_adj.pl $smodel 1 1 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+
+ # SOCAL: coverage kernels for m16 (CC banana-doughnut adjoint sources)
+ print CSH "run_mt_measure_adj.pl $smodel 3 3 $ibp $Ts $lcut $tvec 0 1 0 $iparbools $par1 $par2 $par3\n";
+
# copy run files into RUN directory
$infile1 = "window_chi"; $outfile1 = "${dir_meas}/${ftag}_${infile1}";
$infile2 = "run_file"; $outfile2 = "${dir_meas}/${ftag}_${infile2}";
More information about the CIG-COMMITS
mailing list