[cig-commits] r12548 - in seismo/3D/ADJOINT_TOMO/measure_adj: PLOTS scripts_tomo
carltape at geodynamics.org
carltape at geodynamics.org
Wed Aug 6 09:12:54 PDT 2008
Author: carltape
Date: 2008-08-06 09:12:54 -0700 (Wed, 06 Aug 2008)
New Revision: 12548
Modified:
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/scripts_tomo/run_tomo.pl
Log:
Updated plotting scripts to give the option of NOT plotting the adjoint sources. This makes it possible to better analyze the data and synthetics, and to pick out bad measurement windows.
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2008-08-06 15:43:37 UTC (rev 12547)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_adj.pl 2008-08-06 16:12:54 UTC (rev 12548)
@@ -2,6 +2,7 @@
#-----------------------------------
# Min Chen, 09-July-2007
+# Carl Tape, 02-Aug-2008
# plot_win_adj.pl
#
# This script plots a single figure for a single source-receiver pair
@@ -64,16 +65,20 @@
if (!$opt_i) {$opt_i = "m00"}
$smodel = $opt_i;
-# name of file containing measurements (see mt_sub.f90)
-$meas_file = "window_chi";
-if (not -f $meas_file) {die("\n check if $meas_file exists\n")}
-
# suffix for adjoint sources
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");}
+#---------------------------
+# ADDITIONAL USER INPUT
+
+# name of file containing measurements (see mt_sub.f90)
+$meas_file = "window_chi";
+if (not -f $meas_file) {die("\n check if $meas_file exists\n")}
+
+# plot faults (southern California)
$ifault = 1;
# with iboth = 1, you need to have made both the CC and MT adjoint sources
@@ -82,6 +87,14 @@
$ktitle = "TT -- MT = black ; CC = red";
}
+# width of seismograms depends on whether you plot adjoint sources
+$iplot_adj = 0; # plot adjoint sources
+$iplot_win = 1; # plot windows
+$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};
+
#---------------------------
@M[0..5] = @cmtvals[12..17];
@@ -123,7 +136,7 @@
$proj="-JM4";
# plot the base map
-$origin = "-X3.0";
+$origin = "-X2.75";
`psbasemap $bounds $proj $tick $origin -K -P > $ps_file`; # START
`pscoast $bounds $proj $cinfo -S150/200/255 -W2 -G200/255/150 -N1/1p -N2/0.5p -O -K >> $ps_file`;
@@ -143,7 +156,7 @@
`awk \'{print \$4,\$3,0,4,B,C}\' $opt_a | psxy $proj $bounds -M -St0.15 -W2 -G200/200/200 -N -O -K >>$ps_file`;
#------------
-# data and synthetics and reconstructed synthetics
+# data, synthetics, and (optional) reconstructed synthetics
$suffix = ".d"; # USER MUST CHANGE
@@ -232,7 +245,7 @@
$Smax = 2;
$Srange = $Smax - $Smin;
-$xtextpos1 = -0.10*$tmax; # ymax label
+$xtextpos1 = -0.07*$tmax; # ymax label
$xtextpos2 = 0.05*$tmax; # first window measurement label
$xtextpos3 = 0.50*$tmax; # title position
$xtextpos4 = 0.95*$tmax; # final window measurement label
@@ -246,8 +259,9 @@
# pstext info for measurement labels
$textinfo = "-C2p -W0/255/255o,0/255/255 -G0/0/0";
-$proj = "-JX3/1";
-$dX = 4;
+# dimension of seismograms
+$proj = "-JX${swid}/1";
+$dX = $swid + 0.5;
$dY = 1.5;
$tick = "-Ba50f10:\"Time (s)\":/a1f1S";
@@ -426,202 +440,260 @@
#==========================================
# TRANSVERSE component: data, synthetics, and windows
-`pssac2 -X-2 -Y4.5 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $tick >> $ps_file`; # synthetics
+`pssac2 -X-2 -Y4.5 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $tick >> $ps_file`; # synthetics
if ($ntline) {
- for ($i=0; $i<$nT; $i++) {
- open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
- print GMT "$Twinb[$i] 0\n $Twine[$i] 0\n $Twine[$i] 2\n $Twinb[$i] 2\n";
- #print GMT "$Twinb[$i] 0\n";
- #print GMT "$Twine[$i] 0\n";
- #print GMT "$Twine[$i] 2\n";
- #print GMT "$Twinb[$i] 2\n";
- close(GMT);
+ if ($iplot_win==1) {
+ for ($i=0; $i<$nT; $i++) {
+ open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
+ print GMT "$Twinb[$i] 0\n $Twine[$i] 0\n $Twine[$i] 2\n $Twinb[$i] 2\n";
+ #print GMT "$Twinb[$i] 0\n";
+ #print GMT "$Twine[$i] 0\n";
+ #print GMT "$Twine[$i] 2\n";
+ #print GMT "$Twinb[$i] 2\n";
+ close(GMT);
+ }
}
- for ($i=0; $i<$nT; $i++) {
- $xtext = $Twinb[$i]; $just = "LB";
- if($i == 0) {$xtext = $xtextpos2; $just = "LB";}
- if($i == $nT-1) {$xtext = $xtextpos4; $just = "RB";}
- `echo \"$xtext $ytextpos4 7 0 1 $just $stlabT_dT[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
- `echo \"$xtext $ytextpos3 7 0 1 $just $stlabT_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ if ($iplot_CClabel) {
+ for ($i=0; $i<$nT; $i++) {
+ $xtext = $Twinb[$i]; $just = "LB";
+ if ($i == 0) {
+ $xtext = $xtextpos2; $just = "LB";
+ }
+ if ($i == $nT-1) {
+ $xtext = $xtextpos4; $just = "RB";
+ }
+ `echo \"$xtext $ytextpos4 7 0 1 $just $stlabT_dT[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ `echo \"$xtext $ytextpos3 7 0 1 $just $stlabT_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ }
}
}
-if($Tflag==1) {`pssac2 $proj $Tdat $bounds -Ent-3 -M${sizeT} $black_pen -N -K -O >> $ps_file`;} # data
-`pssac2 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O >> $ps_file`; # synthetics
-if($ntline) {`pssac2 $proj $Trecon $bounds -Ent-3 -M${sizeT3} $recon_pen -N -K -O >> $ps_file`;} # reconstructed
+if ($Tflag==1) {`pssac2 $proj $Tdat $bounds -Ent-3 -M${sizeT} $black_pen -N -K -O >> $ps_file`;} # data
+`pssac2 $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O >> $ps_file`; # synthetics
+if ($ntline && ($iplot_recon==1)) {`pssac2 $proj $Trecon $bounds -Ent-3 -M${sizeT3} $recon_pen -N -K -O >> $ps_file`;} # reconstructed
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC T\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strT\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# RADIAL component: data, synthetics, and windows
`pssac2 -Y$dY $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O $tick2 >> $ps_file`;
if ($nrline) {
- for ($i=0; $i<$nR; $i++) {
- open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
- print GMT "$Rwinb[$i] 0\n $Rwine[$i] 0\n $Rwine[$i] 2\n $Rwinb[$i] 2\n";
- close(GMT);
+ if ($iplot_win==1) {
+ for ($i=0; $i<$nR; $i++) {
+ open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
+ print GMT "$Rwinb[$i] 0\n $Rwine[$i] 0\n $Rwine[$i] 2\n $Rwinb[$i] 2\n";
+ close(GMT);
+ }
}
- for ($i=0; $i<$nR; $i++) {
- $xtext = $Rwinb[$i]; $just = "LB";
- if($i == 0) {$xtext = $xtextpos2; $just = "LB";}
- if($i == $nR-1) {$xtext = $xtextpos4; $just = "RB";}
- `echo \"$xtext $ytextpos4 7 0 1 $just $stlabR_dT[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
- `echo \"$xtext $ytextpos3 7 0 1 $just $stlabR_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ if ($iplot_CClabel==1) {
+ for ($i=0; $i<$nR; $i++) {
+ $xtext = $Rwinb[$i]; $just = "LB";
+ if ($i == 0) {
+ $xtext = $xtextpos2; $just = "LB";
+ }
+ if ($i == $nR-1) {
+ $xtext = $xtextpos4; $just = "RB";
+ }
+ `echo \"$xtext $ytextpos4 7 0 1 $just $stlabR_dT[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ `echo \"$xtext $ytextpos3 7 0 1 $just $stlabR_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ }
}
}
-if($Rflag==1) {`pssac2 $proj $Rdat $bounds -Ent-3 -M${sizeR} $black_pen -N -K -O >> $ps_file`;}
+if ($Rflag==1) {`pssac2 $proj $Rdat $bounds -Ent-3 -M${sizeR} $black_pen -N -K -O >> $ps_file`;}
`pssac2 $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O >> $ps_file`;
-if($nrline) {`pssac2 $proj $Rrecon $bounds -Ent-3 -M${sizeR3} $recon_pen -N -K -O >> $ps_file`;}
+if ($nrline && ($iplot_recon==1)) {`pssac2 $proj $Rrecon $bounds -Ent-3 -M${sizeR3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC R\" | pstext $proj $bounds -N -O -K >> $ps_file`;
`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strR\" | pstext $proj $bounds -N -O -K >> $ps_file`;
# VERTICAL component: data, synthetics, and windows
`pssac2 -Y$dY $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -K -O $tick2 >> $ps_file`;
if ($nzline) {
- for ($i=0; $i<$nZ; $i++) {
- open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
- print GMT "$Zwinb[$i] 0\n $Zwine[$i] 0\n $Zwine[$i] 2\n $Zwinb[$i] 2\n";
- close(GMT);
+ if ($iplot_win==1) {
+ for ($i=0; $i<$nZ; $i++) {
+ open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
+ print GMT "$Zwinb[$i] 0\n $Zwine[$i] 0\n $Zwine[$i] 2\n $Zwinb[$i] 2\n";
+ close(GMT);
+ }
}
- for ($i=0; $i<$nZ; $i++) {
- $xtext = $Zwinb[$i]; $just = "LB";
- if($i == 0) {$xtext = $xtextpos2; $just = "LB";}
- if($i == $nZ-1) {$xtext = $xtextpos4; $just = "RB";}
- `echo \"$xtext $ytextpos4 7 0 1 $just $stlabZ_dT[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
- `echo \"$xtext $ytextpos3 7 0 1 $just $stlabZ_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ if ($iplot_CClabel==1) {
+ for ($i=0; $i<$nZ; $i++) {
+ $xtext = $Zwinb[$i]; $just = "LB";
+ if ($i == 0) {
+ $xtext = $xtextpos2; $just = "LB";
+ }
+ if ($i == $nZ-1) {
+ $xtext = $xtextpos4; $just = "RB";
+ }
+ `echo \"$xtext $ytextpos4 7 0 1 $just $stlabZ_dT[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ `echo \"$xtext $ytextpos3 7 0 1 $just $stlabZ_dA[$i]\" | pstext $textinfo $proj $bounds -N -O -K >> $ps_file`;
+ }
}
}
-if($Zflag==1) {`pssac2 $proj $Zdat $bounds -Ent-3 -M${sizeZ} $black_pen -N -K -O $tick2 >> $ps_file`;}
+if ($Zflag==1) {`pssac2 $proj $Zdat $bounds -Ent-3 -M${sizeZ} $black_pen -N -K -O $tick2 >> $ps_file`;}
`pssac2 $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -O -K >> $ps_file`;
-if($nzline) {`pssac2 $proj $Zrecon $bounds -Ent-3 -M${sizeZ3} $recon_pen -N -K -O >> $ps_file`;}
+if ($nzline && ($iplot_recon==1)) {`pssac2 $proj $Zrecon $bounds -Ent-3 -M${sizeZ3} $recon_pen -N -K -O >> $ps_file`;}
`echo \"$xtextpos1 $ytextpos2 14 0 1 BC Z\" | pstext $proj $bounds -N -K -O >> $ps_file`;
`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 vs Synthetics\" | pstext $proj $bounds -N -K -O >> $ps_file`;
-`echo \"$xtextpos3 $ytextpos5 12 0 $fontno BC Adjoint Source ($ktitle)\" | pstext $proj $bounds -N -X$dX -K -O >> $ps_file`;
+`echo \"$xtextpos3 $ytextpos5 12 0 $fontno BC Data (Black) and Synthetics (Red)\" | 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`;
+}
#----------------------------------------
+# PLOTTING ADJOINT SOURCES
-$tick="-Ba50f10:\"Time (s)\":/a1f1S";
-$tick2="-Ba50f10/a1f1S";
-$proj="-JX3/1";
-#$red_pen = "-W1p,250/0/0";
-#$black_pen = "-W1p,0/0/0";
+if($iplot_adj == 1) {
-# adjoint sources (always have the BH_ channel, even if the data do not)
-$Zadj="ADJOINT_SOURCES/$sta.$net.BHZ.${klab}.adj";
-$Radj="ADJOINT_SOURCES/$sta.$net.BHR.${klab}.adj";
-$Tadj="ADJOINT_SOURCES/$sta.$net.BHT.${klab}.adj";
+ $tick="-Ba50f10:\"Time (s)\":/a1f1S";
+ $tick2="-Ba50f10/a1f1S";
+ $proj="-JX${swid}/1";
+ #$red_pen = "-W1p,250/0/0";
+ #$black_pen = "-W1p,0/0/0";
-#$maxZ=1.0; $maxR=1.0; $maxT=1.0;
+ # adjoint sources (always have the BH_ channel, even if the data do not)
+ $Zadj="ADJOINT_SOURCES/$sta.$net.BHZ.${klab}.adj";
+ $Radj="ADJOINT_SOURCES/$sta.$net.BHR.${klab}.adj";
+ $Tadj="ADJOINT_SOURCES/$sta.$net.BHT.${klab}.adj";
-$scale=1.2; # adjust for plotting (this times max gives the axis limits)
-if (-f $Zadj) {
- ($minZ,$maxZ)=split(/\//,`minmax $Zadj | awk '{print \$6}'`);
- (undef,$minZ)=split(/</,$minZ);($maxZ,undef)=split(/>/,$maxZ);
- $maxZ=$scale*max(abs($maxZ),abs($minZ)); }
-if (-f $Radj) {
- ($minR,$maxR)=split(/\//,`minmax $Radj | awk '{print \$6}'`);
- (undef,$minR)=split(/</,$minR);($maxR,undef)=split(/>/,$maxR);
- $maxR=$scale*max(abs($maxR),abs($minR)); }
-if (-f $Tadj) {
- ($minT,$maxT)=split(/\//,`minmax $Tadj | awk '{print \$6}'`);
- (undef,$minT)=split(/</,$minT);($maxT,undef)=split(/>/,$maxT);
- $maxT=$scale*max(abs($maxT),abs($minT)); }
+ #$maxZ=1.0; $maxR=1.0; $maxT=1.0;
-# overall max value
-$maxA = max($maxT,max($maxR,$maxZ));
-#print "\n Overall adjoint source max value: $maxA\n";
-if ($opt_r) {$maxZ=$maxA;$maxR=$maxA;$maxT=$maxA;}
+ $scale=1.2; # adjust for plotting (this times max gives the axis limits)
+ if (-f $Zadj) {
+ ($minZ,$maxZ)=split(/\//,`minmax $Zadj | awk '{print \$6}'`);
+ (undef,$minZ)=split(/</,$minZ);($maxZ,undef)=split(/>/,$maxZ);
+ $maxZ=$scale*max(abs($maxZ),abs($minZ));
+ }
+ if (-f $Radj) {
+ ($minR,$maxR)=split(/\//,`minmax $Radj | awk '{print \$6}'`);
+ (undef,$minR)=split(/</,$minR);($maxR,undef)=split(/>/,$maxR);
+ $maxR=$scale*max(abs($maxR),abs($minR));
+ }
+ if (-f $Tadj) {
+ ($minT,$maxT)=split(/\//,`minmax $Tadj | awk '{print \$6}'`);
+ (undef,$minT)=split(/</,$minT);($maxT,undef)=split(/>/,$maxT);
+ $maxT=$scale*max(abs($maxT),abs($minT));
+ }
-# CHT: 10e-3 tolerance too large for waveform difference
-$tol = 10e-8;
-if ($maxZ<$tol) {$maxZ=1.0;}
-if ($maxR<$tol) {$maxR=1.0;}
-if ($maxT<$tol) {$maxT=1.0;}
+ # overall max value
+ $maxA = max($maxT,max($maxR,$maxZ));
+ #print "\n Overall adjoint source max value: $maxA\n";
+ if ($opt_r) {
+ $maxZ=$maxA;$maxR=$maxA;$maxT=$maxA;
+ }
-#$maxZ=1.0; $maxR=1.0; $maxT=1.0; # testing
+ # CHT: 10e-3 tolerance too large for waveform difference
+ $tol = 10e-8;
+ if ($maxZ<$tol) {
+ $maxZ=1.0;
+ }
+ if ($maxR<$tol) {
+ $maxR=1.0;
+ }
+ if ($maxT<$tol) {
+ $maxT=1.0;
+ }
-# strings for plotting
-$strT = sprintf('%.2e',$maxT);
-$strR = sprintf('%.2e',$maxR);
-$strZ = sprintf('%.2e',$maxZ);
+ #$maxZ=1.0; $maxR=1.0; $maxT=1.0; # testing
-# adjoint sources for cross-correlation plot in red
-if($iboth==1) {
- $Zadj2="ADJOINT_SOURCES/$sta.$net.BHZ.cc.adj";
- $Radj2="ADJOINT_SOURCES/$sta.$net.BHR.cc.adj";
- $Tadj2="ADJOINT_SOURCES/$sta.$net.BHT.cc.adj";
-}
+ # strings for plotting
+ $strT = sprintf('%.2e',$maxT);
+ $strR = sprintf('%.2e',$maxR);
+ $strZ = sprintf('%.2e',$maxZ);
-# VERTICAL component: adjoint source and windows
-$bounds="-R$tmin/$tmax/-$maxZ/$maxZ";
-`psbasemap $proj $bounds $tick2 -K -O >> $ps_file`;
-if ($nzline) {
-for($i =0; $i<$nZ; $i++){
- open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
- print GMT "$Zwinb[$i] -$maxZ\n";
- print GMT "$Zwine[$i] -$maxZ\n";
- print GMT "$Zwine[$i] $maxZ\n";
- print GMT "$Zwinb[$i] $maxZ\n";
- close(GMT);
- }}
-if (-f $Zadj) {
- $ytextpos2 = 0; $ytextpos1 = -$maxZ + 0.15*(2*$maxZ);
- `psxy $proj $Zadj $bounds $black_pen -N -K -O >> $ps_file`;
- if($iboth==1) {`psxy $proj $Zadj2 $bounds $red_pen -N -K -O >> $ps_file`;}
- `echo \"$xtextpos1 $ytextpos2 14 0 1 BC Z\" | pstext $proj $bounds -N -K -O >> $ps_file`;
- `echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strZ\" | pstext $proj $bounds -N -K -O >> $ps_file`;
-}
+ # adjoint sources for cross-correlation plot in red
+ if ($iboth==1) {
+ $Zadj2="ADJOINT_SOURCES/$sta.$net.BHZ.cc.adj";
+ $Radj2="ADJOINT_SOURCES/$sta.$net.BHR.cc.adj";
+ $Tadj2="ADJOINT_SOURCES/$sta.$net.BHT.cc.adj";
+ }
-# RADIAL component: adjoint source and windows
-$bounds="-R$tmin/$tmax/-$maxR/$maxR";
-`psbasemap -Y-$dY $proj $bounds $tick2 -K -O >> $ps_file`;
-if ($nrline) {
-for($i =0; $i<$nR; $i++){
- open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
- print GMT "$Rwinb[$i] -$maxR\n";
- print GMT "$Rwine[$i] -$maxR\n";
- print GMT "$Rwine[$i] $maxR\n";
- print GMT "$Rwinb[$i] $maxR\n";
- close(GMT);
- }}
-if (-f $Radj) {
- $ytextpos2 = 0; $ytextpos1 = -$maxR + 0.15*(2*$maxR);
- `psxy $proj $Radj $bounds $black_pen -N -K -O >> $ps_file`;
- `echo \"$xtextpos1 $ytextpos2 14 0 1 BC R\" | pstext $proj $bounds -N -K -O >> $ps_file`;
- if($iboth==1) {`psxy $proj $Radj2 $bounds $red_pen -N -K -O >> $ps_file`;}
- `echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strR\" | pstext $proj $bounds -N -K -O >> $ps_file`;
-}
+ # VERTICAL component: adjoint source and windows
+ $bounds="-R$tmin/$tmax/-$maxZ/$maxZ";
+ `psbasemap $proj $bounds $tick2 -K -O >> $ps_file`;
+ if ($iplot_win==1) {
+ if ($nzline) {
+ for ($i =0; $i<$nZ; $i++) {
+ open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
+ print GMT "$Zwinb[$i] -$maxZ\n";
+ print GMT "$Zwine[$i] -$maxZ\n";
+ print GMT "$Zwine[$i] $maxZ\n";
+ print GMT "$Zwinb[$i] $maxZ\n";
+ close(GMT);
+ }
+ }
+ }
+ if (-f $Zadj) {
+ $ytextpos2 = 0; $ytextpos1 = -$maxZ + 0.15*(2*$maxZ);
+ `psxy $proj $Zadj $bounds $black_pen -N -K -O >> $ps_file`;
+ if ($iboth==1) {
+ `psxy $proj $Zadj2 $bounds $red_pen -N -K -O >> $ps_file`;
+ }
+ `echo \"$xtextpos1 $ytextpos2 14 0 1 BC Z\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ `echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strZ\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ }
-# TRANSVESE component: adjoint source and windows
-$bounds="-R$tmin/$tmax/-$maxT/$maxT";
-`psbasemap -Y-$dY $proj $bounds $tick -K -O >> $ps_file`;
-if ($ntline) {
-for($i =0; $i<$nT; $i++){
- open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
- print GMT "$Twinb[$i] -$maxT\n";
- print GMT "$Twine[$i] -$maxT\n";
- print GMT "$Twine[$i] $maxT\n";
- print GMT "$Twinb[$i] $maxT\n";
- close(GMT);
- }}
-if (-f $Tadj) {
- $ytextpos2 = 0; $ytextpos1 = -$maxT + 0.15*(2*$maxT);
- `psxy $proj $Tadj $bounds $black_pen -N -K -O >> $ps_file`;
- if($iboth==1) {`psxy $proj $Tadj2 $bounds $red_pen -N -K -O >> $ps_file`;}
- `echo \"$xtextpos1 $ytextpos2 14 0 1 BC T\" | pstext $proj $bounds -N -K -O >> $ps_file`;
- `echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strT\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ # RADIAL component: adjoint source and windows
+ $bounds="-R$tmin/$tmax/-$maxR/$maxR";
+ `psbasemap -Y-$dY $proj $bounds $tick2 -K -O >> $ps_file`;
+ if ($iplot_win==1) {
+ if ($nrline) {
+ for ($i =0; $i<$nR; $i++) {
+ open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
+ print GMT "$Rwinb[$i] -$maxR\n";
+ print GMT "$Rwine[$i] -$maxR\n";
+ print GMT "$Rwine[$i] $maxR\n";
+ print GMT "$Rwinb[$i] $maxR\n";
+ close(GMT);
+ }
+ }
+ }
+ if (-f $Radj) {
+ $ytextpos2 = 0; $ytextpos1 = -$maxR + 0.15*(2*$maxR);
+ `psxy $proj $Radj $bounds $black_pen -N -K -O >> $ps_file`;
+ `echo \"$xtextpos1 $ytextpos2 14 0 1 BC R\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ if ($iboth==1) {
+ `psxy $proj $Radj2 $bounds $red_pen -N -K -O >> $ps_file`;
+ }
+ `echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strR\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ }
+
+ # TRANSVESE component: adjoint source and windows
+ $bounds="-R$tmin/$tmax/-$maxT/$maxT";
+ `psbasemap -Y-$dY $proj $bounds $tick -K -O >> $ps_file`;
+ if ($iplot_win==1) {
+ if ($ntline) {
+ for ($i =0; $i<$nT; $i++) {
+ open(GMT,"|psxy $proj $bounds -W2 -G220/220/220 -O -K >> $ps_file");
+ print GMT "$Twinb[$i] -$maxT\n";
+ print GMT "$Twine[$i] -$maxT\n";
+ print GMT "$Twine[$i] $maxT\n";
+ print GMT "$Twinb[$i] $maxT\n";
+ close(GMT);
+ }
+ }
+ }
+ if (-f $Tadj) {
+ $ytextpos2 = 0; $ytextpos1 = -$maxT + 0.15*(2*$maxT);
+ `psxy $proj $Tadj $bounds $black_pen -N -K -O >> $ps_file`;
+ if ($iboth==1) {
+ `psxy $proj $Tadj2 $bounds $red_pen -N -K -O >> $ps_file`;
+ }
+ `echo \"$xtextpos1 $ytextpos2 14 0 1 BC T\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ `echo \"$xtextpos1 $ytextpos1 9 0 $fontno CM $strT\" | pstext $proj $bounds -N -K -O >> $ps_file`;
+ }
+
}
+#-----------------------
+
# plot titles
#`echo \"$xtextpos3 $maxT 12 0 $fontno BC Adjoint Source ($ktitle)\" | pstext $proj $bounds -N -Y3.2 -K -O >> $ps_file`;
-#`echo \"$xtextpos3 $maxT 12 0 $fontno BC Data vs Synthetics\" | pstext $proj $bounds -N -O -X-4 >> $ps_file`;
-`echo \"0 0\" | psxy $proj $bounds -N -O -X15 >> $ps_file`;
+#`echo \"$xtextpos3 $maxT 12 0 $fontno BC Data and Synthetics\" | pstext $proj $bounds -N -O -X-4 >> $ps_file`;
+`echo \"0 0\" | psxy $proj $bounds -N -O -X15 >> $ps_file`; # FINISH
-#`convert $ps_file $jpg_file`;
+#`convert $ps_file $jpg_file ; xv $jpg_file &`;
`ps2pdf $ps_file`;
`rm $ps_file`; # remove PS file
-#`xv $jpg_file &`;
#system("xv $jpg_file &");
#================================================
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats.pl 2008-08-06 15:43:37 UTC (rev 12547)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats.pl 2008-08-06 16:12:54 UTC (rev 12548)
@@ -88,7 +88,8 @@
$cmt_scale = 0.4;
$cmtinfo = "-N -Sm${cmt_scale}/$cmtfont -L0.5p/0/0/0 -G255/0/0";
$cmtinfo_dc = "-N -Sd${cmt_scale}/$cmtfont -L0.5p/0/0/0 -G255/0/0";
- at psmeca_line = `grep $eid $cmtall_psmeca`; chomp(@psmeca_line);
+ at psmeca_line = `grep $eid ${cmtall_psmeca}`; chomp(@psmeca_line);
+print "\n -- @psmeca_line -- \n";
print CSH "psmeca $R $J $cmtinfo -O -K -V >>$psfile<<EOF\n @psmeca_line \nEOF\n";
#-----------------------------------------
@@ -130,8 +131,8 @@
# KEY: commands for histograms axes
$Nbin = 20; # number of bins
-$max_DT = 12; $min_DT = -$max_DT;
-$max_DA = 3.5; $min_DA = -$max_DA;
+$max_DT = 8; $min_DT = -$max_DT;
+$max_DA = 2.5; $min_DA = -$max_DA;
$max_sigma_DT = 2.5; $min_sigma_DT = 0;
$max_sigma_DA = $max_sigma_DT; $min_sigma_DA = $min_sigma_DT ;
$max_DTerr = $max_DT; $min_DTerr = -$max_DTerr;
@@ -164,7 +165,7 @@
$origin_hist5 = "-Xa$ohistX[4] -Ya$ohistY[4]";
$origin_hist6 = "-Xa$ohistX[5] -Ya$ohistY[5]";
- $max_count = 30; $count_tick = 5;
+ $max_count = 50; $count_tick = 10;
$Rhist1 = "-R$min_DT/$max_DT/0/$max_count";
$Rhist2 = "-R$min_DA/$max_DA/0/$max_count";
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 2008-08-06 15:43:37 UTC (rev 12547)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/PLOTS/plot_win_stats_all.pl 2008-08-06 16:12:54 UTC (rev 12548)
@@ -10,10 +10,10 @@
# The set of final histograms are copied to an output directory.
#
# EXAMPLE:
-# plot_win_stats_all.pl 6/40 m0 0
-# plot_win_stats_all.pl 6/40 m0 1 (plot histogram)
-# plot_win_stats_all.pl 2/40 m0 0
-# plot_win_stats_all.pl 2/40 m0 1 (plot histogram)
+# plot_win_stats_all.pl 6/40 m07 0
+# plot_win_stats_all.pl 6/40 m07 1 (plot histogram)
+# plot_win_stats_all.pl 2/40 m07 0
+# plot_win_stats_all.pl 2/40 m07 1 (plot histogram)
#
#-----------------------------------
@@ -23,11 +23,13 @@
$pwd = $ENV{PWD};
($Tmin,$Tmax) = split("/",$Ts);
-$sTmin = sprintf("T%2.2i",$Tmin);
-$sTmax = sprintf("T%2.2i",$Tmax);
+$sTmin = sprintf("T%3.3i",$Tmin);
+$sTmax = sprintf("T%3.3i",$Tmax);
$sTper = sprintf("T = [%i s, %i s]",$Tmin,$Tmax);
-$file_eids = "/net/sierra/raid1/carltape/results/SOURCES/socal_5/SOCAL_FINAL_CMT_v5_eid";
+$dir_source = "/home/carltape/results/SOURCES/socal_9";
+$dir_source_text = "${dir_source}/v9_files_text";
+$file_eids = "${dir_source}/SOCAL_FINAL_CMT_v9_eid";
if (not -f $file_eids) {die("\n check if $file_eids exists\n")}
open(IN,$file_eids); @eids = <IN>; $nevent = @eids;
@@ -39,8 +41,7 @@
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)
-$dir_source = "/home/carltape/results/SOURCES/socal_5";
-$cmtall_psmeca = "${dir_source}/CMT_files_post_inverted_mod/CMT_ALL_post_inverted_mod_psmeca_eid";
+$cmtall_psmeca = "${dir_source}/SOCAL_FINAL_CMT_v9_psmeca_eid";
if (not -f $cmtall_psmeca) {die("\n check if $cmtall_psmeca exists\n")}
#$nevent = 5;
@@ -57,8 +58,9 @@
# loop over all events
$k = 0;
$imin = 1; $imax = $nevent; # default
-#$imin = 1; $imax = 3;
-#$imin = 1; $imax = $imin;
+#$imin = 1; $imax = 10;
+$imin = 204; $imax = $imin;
+
for ($i = $imin; $i <= $imax; $i = $i+1) {
$eid = $eids[$i-1]; chomp($eid);
@@ -70,25 +72,33 @@
# name of file containing measurements (see mt_sub.f90)
$dir_meas = "${dir_run}/${eid}/${smodel}/MEASURE_${sTmin}";
$meas_file = "${dir_meas}/${ftag}_window_chi";
+ if (-f $meas_file) {
- # name of stations file
- $sta_file = "${dir_meas}/ADJOINT_${sTmin}/STATIONS_ADJOINT";
+ # name of stations file
+ $sta_file = "${dir_meas}/ADJOINT_${sTmin}/STATIONS_ADJOINT";
- # file containing text for plotting here
- $dir_source_text = "/home/carltape/results/SOURCES/socal_5/CMT_files_post_inverted_text";
- $eid_text = "${dir_source_text}/${eid}_text";
+ # file containing text for plotting here
+ $eid_text = "${dir_source_text}/${eid}_text";
+ if (not -f $eid_text) {
+ die("\n check if eid_text $eid_text exists\n");
+ }
- # plot figure
- $name = sprintf("%3.3i_${ftag}_win_stats",$i);
- `plot_win_stats.pl $Ts $eid $name $meas_file $sta_file $eid_text $cmtall_psmeca`;
- #print "\n plot_win_stats.pl $Ts $eid $name $meas_file $sta_file $eid_text $cmtall_psmeca";
+ # plot figure
+ $name = sprintf("%3.3i_${ftag}_win_stats",$i);
+ print "\n plot_win_stats.pl $Ts $eid $name $meas_file $sta_file $eid_text ${cmtall_psmeca}\n\n";
+ #die("TESTING");
+ `plot_win_stats.pl $Ts $eid $name $meas_file $sta_file $eid_text ${cmtall_psmeca}`;
- # check if the figure was made
- if (-f "${name}.pdf") {
- $k = $k+1;
- print "\n $i, $eid, ${name}.pdf was made";
+ # check if the figure was made
+ if (-f "${name}.pdf") {
+ $k = $k+1;
+ print "\n $i, $eid, ${name}.pdf was made";
+ } else {
+ print "\n $i, $eid, ${name}.pdf was NOT made";
+ }
+
} else {
- print "\n $i, $eid, ${name}.pdf was NOT made";
+ print "meas_file $meas_file does not exist\n";
}
}
@@ -160,16 +170,16 @@
# KEY: commands for histograms axes
$Nbin = 20;
-$max_DT = 12; $min_DT = -$max_DT;
-$max_DA = 3.5; $min_DA = -$max_DA;
+$max_DT = 8.0; $min_DT = -$max_DT;
+$max_DA = 3.0; $min_DA = -$max_DA;
$max_sigma_DT = 2.5; $min_sigma_DT = 0;
$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) {
- $max_DT = 5; $min_DT = -$max_DT;
+ $max_DT = 4; $min_DT = -$max_DT;
$max_DA = 1.5; $min_DA = -$max_DA;
- $max_sigma_DT = 2.0; $min_sigma_DT = 0;
+ $max_sigma_DT = 1.5; $min_sigma_DT = 0;
$max_sigma_DA = $max_sigma_DT; $min_sigma_DA = $min_sigma_DT ;
$max_DTerr = 10; $min_DTerr = -$max_DTerr;
$max_DAerr = 4; $min_DAerr = -$max_DAerr;
@@ -195,7 +205,7 @@
$origin_hist6 = "-Xa$ohistX[5] -Ya$ohistY[5]";
# plot histograms of DA and DT
- $max_count = 30; $count_tick = 5;
+ $max_count = 40; $count_tick = 5;
$Rhist1 = "-R$min_DT/$max_DT/0/$max_count";
$Rhist2 = "-R$min_DA/$max_DA/0/$max_count";
Modified: seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2008-08-06 15:43:37 UTC (rev 12547)
+++ seismo/3D/ADJOINT_TOMO/measure_adj/scripts_tomo/run_tomo.pl 2008-08-06 16:12:54 UTC (rev 12548)
@@ -12,7 +12,8 @@
# EXAMPLE (from working measure_adj directory) :
# socal m06 -->
# run_tomo.pl m07 0 6/40 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/5.0/1.0 0.4/0.4 1.0
-# run_tomo.pl m07 0 2/40 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/4.0/1.0 0.2/0.2 0.5
+# run_tomo.pl m07 0 4/40 0.011 1/1/1 2.0/2.5/3.5/3.0 0.7/0.7/5.0/1.0 0.4/0.4 1.0
+# run_tomo.pl m07 0 2/40 0.011 1/1/1 0.25/0.5/3.5/3.0 0.7/0.7/4.0/1.0 0.4/0.4 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
@@ -53,8 +54,8 @@
#die("testing");
$imin = 1; $imax = $nevent;
-#$imin = 1; $imax = 160;
-$imin = 204; $imax = $imin;
+$imin = 1; $imax = 100;
+#$imin = 204; $imax = $imin;
# EXAMPLE lines from EIDs_simulation_duration
# 5 10222753 300.0 27300 0.26 332.0 4.010
More information about the cig-commits
mailing list