[cig-commits] r13377 - in seismo/3D/ADJOINT_TOMO/iterate_adj: . UTILS UTILS/vtk matlab misfit_plot misfit_plot/OLD_version
carltape at geodynamics.org
carltape at geodynamics.org
Sun Nov 23 16:05:12 PST 2008
Author: carltape
Date: 2008-11-23 16:05:12 -0800 (Sun, 23 Nov 2008)
New Revision: 13377
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/README
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/9983429_T006_T030_LDF_CI_m13_m00_seis.pdf
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/DATA/
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/README
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/compare_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt.m
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt_run.m
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/plot_misfit.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/seis_norm.m
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/setup_misfit_dir.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/README
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/SYN/
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis_all.pl
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/info/
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/MISFIT_in_development/
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/OUTPUT_MISFIT/
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs_pmax.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_local.tcl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_smooth_local.tcl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model_local.tcl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update_local.tcl
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
Log:
Updated run scripts for subspace method. Added scripts for comparing seismograms and window picks generated from two different models.
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/README (from rev 13376, seismo/3D/ADJOINT_TOMO/iterate_adj/info/README)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/README 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,19 @@
+Carl Tape, 23-Nov-2008
+
+iterate_adj is a collection of scripts needed to do adjoint tomography.
+It is the third package in the sequence:
+
+flexwin
+measure_adj
+iterate_adj
+
+Several operations are done on pangu, like computing the Hessian matrix for the subspace method. And several operations are done in Matlab, like choosing the damping parameter.
+
+Thus, the two main directories for scripts are
+
+pangu
+matlab
+
+The scripts in matlab are executed NOT on a cluster, while the scripts in pangu are.
+
+---------------------
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/README
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs.pl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -7,12 +7,12 @@
# This script makes the event kernel plots by calling vtk scripts.
#
# EXAMPLE (unsmoothed kernels):
-# ~/UTILS/tomo_make_figs.pl 1 1 all m10 0 0
-# ~/UTILS/tomo_make_figs.pl 1 0 all m10 0 0 (make vtu files only)
+# ~/UTILS/tomo_make_figs.pl 1 1 all m13 0 0
+# ~/UTILS/tomo_make_figs.pl 1 0 all m13 0 0 (make vtu files only)
#
# EXAMPLE (smoothed kernels):
-# ~/UTILS/tomo_make_figs.pl 1 1 all m10 6 1
-# ~/UTILS/tomo_make_figs.pl 1 0 all m10 6 1 (make vtu files only)
+# ~/UTILS/tomo_make_figs.pl 1 1 all m13 6 1
+# ~/UTILS/tomo_make_figs.pl 1 0 all m13 6 1 (make vtu files only)
#
#-----------------------------------
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs_pmax.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs_pmax.pl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/tomo_make_figs_pmax.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -8,7 +8,7 @@
# by calling vtk scripts.
#
# EXAMPLE:
-# ~/UTILS/tomo_make_figs_pmax.pl 1 0 dm10 0.80 0.0 mu_kernel_smooth beta
+# ~/UTILS/tomo_make_figs_pmax.pl 1 0 dm13 0.80 0.0 mu_kernel_smooth beta_window
#
# ~/UTILS/tomo_make_figs_pmax.pl 1 0 dm08 0.80 0.0 mu_kernel_smooth beta
# ~/UTILS/tomo_make_figs_pmax.pl 1 0 dm08 0.40 0.0 kappa_kernel_smooth bulk
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel.pl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -73,7 +73,7 @@
`cp $cmax_file .`;
open(IN,"cmax_kernel"); $cmax = <IN>; close(IN);
} else {
- $cmax = 5e-11; # default
+ $cmax = 5e-11; # defalt
}
#$cmax = 1e-11;
@@ -144,7 +144,8 @@
@dlayers = (-0.25,0.001,2,4,6,8,10,15,20,25,30,35,40);
$Nz = @dlayers;
-#$imin = 1; $imax = $Nz; # default
+$imin = 1; $imax = $Nz; # default
+$imin = 1; $imax = 10; # testing
$imin = 4; $imax = $imin; # testing
for ($i = $imin; $i <= $imax; $i++) {
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_local.tcl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_local.tcl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_local.tcl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -40,13 +40,13 @@
# === load the unstructured grid data ===
vtkXMLUnstructuredGridReader kReader1
- kReader1 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/14079184/m11/MESH_all/mu_kernel_low_1.vtu
+ kReader1 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/9716853/m13/MESH_all/mu_kernel_low_1.vtu
vtkXMLUnstructuredGridReader kReader2
- kReader2 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/14079184/m11/MESH_all/mu_kernel_low_2.vtu
+ kReader2 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/9716853/m13/MESH_all/mu_kernel_low_2.vtu
vtkXMLUnstructuredGridReader kReader3
- kReader3 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/14079184/m11/MESH_all/mu_kernel_low_3.vtu
+ kReader3 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/9716853/m13/MESH_all/mu_kernel_low_3.vtu
## When the files are merged, everything crashes.
#
@@ -267,7 +267,7 @@
#=== text title ===
vtkTextActor titleActor
titleActor SetDisplayPosition 280 60
- titleActor SetInput "14079184 : m11 kernel for SHEAR MODULUS -- 105 stations -- Cut at z = 4.0 km"
+ titleActor SetInput "9716853 : m13 kernel for SHEAR MODULUS -- 85 stations -- Cut at z = 4.0 km"
set tprop [titleActor GetTextProperty]
$tprop SetJustificationToCentered
$tprop SetColor 0 0 1
@@ -315,7 +315,7 @@
w2i SetInput renWin
vtkPostScriptWriter writer
writer SetInputConnection [w2i GetOutputPort]
- writer SetFileName "117_14079184_mu_all_kernel_04.ps"
+ writer SetFileName "132_9716853_mu_all_kernel_04.ps"
writer Write
wm withdraw .
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_smooth_local.tcl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_smooth_local.tcl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_kernel_smooth_local.tcl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -40,7 +40,7 @@
# === load the unstructured grid data ===
vtkXMLUnstructuredGridReader kReader1
- kReader1 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/SMOOTH_EVENT_KERNELS/m11/14079184/mu_kernel_smooth_h006km_v001km.vtu
+ kReader1 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/SMOOTH_EVENT_KERNELS/m13/9716853/mu_kernel_smooth_h006km_v001km.vtu
# vtkXMLUnstructuredGridReader kReader2
# kReader2 SetFileName vtu_files/lin_model/vp_low_2.vtu
@@ -267,7 +267,7 @@
#=== text title ===
vtkTextActor titleActor
titleActor SetDisplayPosition 280 60
- titleActor SetInput "14079184 : m11 kernel for SHEAR MODULUS -- 105 stations -- Cut at z = 4.0 km"
+ titleActor SetInput "9716853 : m13 kernel for SHEAR MODULUS -- 85 stations -- Cut at z = 4.0 km"
set tprop [titleActor GetTextProperty]
$tprop SetJustificationToCentered
$tprop SetColor 0 0 1
@@ -315,7 +315,7 @@
w2i SetInput renWin
vtkPostScriptWriter writer
writer SetInputConnection [w2i GetOutputPort]
- writer SetFileName "117_14079184_mu_all_kernel_h006km_v001km_04.ps"
+ writer SetFileName "132_9716853_mu_all_kernel_h006km_v001km_04.ps"
writer Write
wm withdraw .
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model.pl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -9,17 +9,21 @@
#
# CALLED BY: tomo_make_figs.pl
#
-# EXAMPLE:
-# view_model.pl 5 vs_m12
-# view_model.pl 5 vs
-# view_model.pl 7 vs_m12_m00
-# view_model.pl 7 vs_m12_m11
-# view_model.pl 5 vs_low
-# view_model.pl 7 Beta
+# EXAMPLES:
#
-# view_model.pl 4 vp_m12
-# view_model.pl 4 vp
-# view_model.pl 6 vp_m12_m00
+# view_model.pl 4 vp_m14
+# view_model.pl 4 vp_m00
+# view_model.pl 5 vs_m14
+# view_model.pl 5 vs_m00
+# view_model.pl 6 vb_m14
+# view_model.pl 6 vb_m00
+# view_model.pl 7 poisson_m14
+# view_model.pl 7 poisson_m00
+# view_model.pl 8 vp_m14_m13
+# view_model.pl 8 vp_m14_m00
+# view_model.pl 9 vs_m14_m13
+# view_model.pl 9 vs_m14_m00
+# view_model.pl 10 vb_m14_m00
#
#---------------------------------
@@ -43,10 +47,12 @@
#`cp /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/m09/*vtu .`; $title_tag = "Model m09 for";
#`cp /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/m10/*vtu .`; $title_tag = "Model m10 for";
#`cp /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/m11/*vtu .`; $title_tag = "Model m11 for";
-`cp /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/m12/*vtu .`; $title_tag = "Model m12 for";
+#`cp /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/m12/*vtu .`; $title_tag = "Model m12 for";
+#`cp /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/m13/*vtu .`; $title_tag = "Model m13 for";
+`cp /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/m14/*vtu .`; $title_tag = "Model m14 for";
-$title_tag = "ln(m12\\/m11) for"; $cmax = 0.05;
-#$title_tag = "ln(m12\\/m00) for"; $cmax = 0.30;
+$title_tag = "LN(m14\\/m13) for"; $cmax = 0.05;
+$title_tag = "LN(m14\\/m00) for"; $cmax = 0.20;
$tcl_tag = "view_model";
$tcl_file = "${tcl_tag}.tcl";
@@ -62,8 +68,8 @@
#$dir_ker_lab = "\\/net\\/sierra\\/raid1\\/carltape\\/socal\\/socal_3D\\/RUNS\\/${eid}\\/${smodel}\\/MESH_${ftag}";
# kernel options
- at klabs = ("kappa","mu","rho","vp","vs","vp","vs");
- at ktitles = ("BULK MODULUS","SHEAR MODULUS","DENSITY","P WAVE-SPEED","S WAVE-SPEED","P WAVE-SPEED","S WAVE-SPEED");
+ at klabs = ("kappa","mu","rho","vp","vs","vb","nu","vp","vs","vb");
+ at ktitles = ("BULK MODULUS","SHEAR MODULUS","DENSITY","P WAVE-SPEED","S WAVE-SPEED","BULK WAVE-SPEED","POISSON RATIO","P WAVE-SPEED","S WAVE-SPEED","BULK WAVE-SPEED");
$klab = $klabs[$iker-1];
$ktitle = $ktitles[$iker-1];
#$file_tag = "${klab}_low"; $ftag = "${klab}_low";
@@ -96,6 +102,11 @@
# percent perturbation for each layer
@players = (15,15,15,15,15,15,10,10,10,10,10,10,5);
+# mean values for each depth slice
+ at cmean_vp = (5.20,5.20,5.50,5.9,6.11,6.11,6.22,6.46,6.50,6.55,6.7,7.30,7.80);
+ at cmean_vs = (2.95,2.95,3.1,3.25,3.4,3.45,3.55,3.65,3.7,3.7,3.75,4.2,4.5);
+#@cmean_vb = (3.92,3.92,4.17,4.55,4.68,4.63,4.67,4.89,4.89,4.96,5.11,5.45,5.81);
+
$pert = 10;
if($iker == 1) {
$cmin = 0e10; $cmax = 10e10; $cunits = "Pa";
@@ -106,17 +117,34 @@
} elsif($iker == 4) {
$cmin = 1000; $cmax = 6000; $cunits = "m\\/s";
- @cmean = (5.20,5.20,5.50,5.9,6.11,6.11,6.22,6.46,6.50,6.55,6.7,7.30,7.80);
+ @cmean = @cmean_vp;
} elsif($iker == 5) {
$cmin = 1000; $cmax = 5000; $cunits = "m\\/s";
- @cmean = (2.95,2.95,3.1,3.25,3.4,3.45,3.55,3.65,3.7,3.7,3.75,4.2,4.5);
+ @cmean = @cmean_vs;
} elsif($iker == 6) {
- $cmin = -$cmax; $cunits = " ";
+ $cmin = 1000; $cmax = 5000; $cunits = "m\\/s";
+ #@cmean = @cmean_vb;
+ for ($i = 1; $i <= $Nz; $i++) {
+ $vp = $cmean_vp[$i-1];
+ $vs = $cmean_vs[$i-1];
+ $cmean[$i-1] = sqrt( $vp**2 - (4/3)*$vs**2 );
+ }
+ #print "\n @cmean_vp \n @cmean_vs \n @cmean \n"; die("TESTING");
} elsif($iker == 7) {
+ $cmin = 0.1; $cmax = 0.4; $cunits = " ";
+
+} elsif($iker == 8) {
$cmin = -$cmax; $cunits = " ";
+
+} elsif($iker == 9) {
+ $cmin = -$cmax; $cunits = " ";
+
+} elsif($iker == 10) {
+ $cmin = -$cmax; $cunits = " ";
+
}
# corners of the UTM mesh
@@ -171,7 +199,7 @@
$Nz = @dlayers;
$imin = 1; $imax = $Nz; # default
#$imin = 9; $imax = 13; # testing
-#$imin = 3; $imax = $imin; # testing
+#$imin = 2; $imax = $imin; # testing
for ($i = $imin; $i <= $imax; $i++) {
@@ -185,11 +213,16 @@
# title for plot
$title = sprintf("%s -- Cut at z = %.1f km",$tlab,-$zcen_km);
- if($iker <= 5) {
+ if($iker <= 6) {
$clabel = sprintf("%.0f $cunits +- %.0f percent",$cmid*1000,$pert);
- $scalar_low = $cmid* 1000 * (1 - $pert/100);
+ $scalar_low = $cmid * 1000 * (1 - $pert/100);
$scalar_high = $cmid * 1000 * (1 + $pert/100);
+ } elsif($iker == 7) {
+ $clabel = " ";
+ $scalar_low = sprintf("%.4e",$cmin);
+ $scalar_high = sprintf("%.4e",$cmax);
+
} else {
$clabel = "Percent Change";
$scalar_low = sprintf("%.4e",$cmin);
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model_local.tcl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model_local.tcl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_model_local.tcl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -40,7 +40,7 @@
# === load the unstructured grid data ===
vtkXMLUnstructuredGridReader kReader1
- kReader1 SetFileName vs_m12_m11.vtu
+ kReader1 SetFileName vs_m13_m00.vtu
vtkXMLUnstructuredGridReader kReader2
kReader2 SetFileName vtu_files/lin_model/vp_low_2.vtu
@@ -67,7 +67,7 @@
# data mappers
vtkDataSetMapper kMapper1
kMapper1 SetInput [kReader1 GetOutput]
- kMapper1 SetScalarRange -5.0000e-02 5.0000e-02
+ kMapper1 SetScalarRange -2.0000e-01 2.0000e-01
# data actors
vtkActor kActor1
@@ -83,7 +83,7 @@
vtkDataSetMapper hrMapper1
hrMapper1 SetInput [hrCut1 GetOutput]
hrMapper1 InterpolateScalarsBeforeMappingOn
- hrMapper1 SetScalarRange -5.0000e-02 5.0000e-02
+ hrMapper1 SetScalarRange -2.0000e-01 2.0000e-01
hrMapper1 SetLookupTable lut
vtkActor hrActor1
hrActor1 SetMapper hrMapper1
@@ -176,7 +176,7 @@
#=== text title ===
vtkTextActor titleActor
titleActor SetDisplayPosition 280 60
- titleActor SetInput "ln(m12/m11) for S WAVE-SPEED -- Cut at z = 40.0 km"
+ titleActor SetInput "LN(m13/m00) for S WAVE-SPEED -- Cut at z = 40.0 km"
set tprop [titleActor GetTextProperty]
$tprop SetJustificationToCentered
$tprop SetColor 0 0 1
@@ -220,7 +220,7 @@
w2i SetInput renWin
vtkPostScriptWriter writer
writer SetInputConnection [w2i GetOutputPort]
- writer SetFileName "vs_m12_m11_13.ps"
+ writer SetFileName "vs_m13_m00_13.ps"
writer Write
wm withdraw .
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update.pl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -44,6 +44,10 @@
# view_update.pl 2 dm07 44 1.00 0 dm_mu_kernel_smooth_p044 beta
# view_update.pl 2 dm10 28 1.60 0 dm_mu_kernel_smooth_p028 beta
# view_update.pl 2 dm11 40 1.60 0 dm_mu_kernel_smooth_p040 beta
+# view_update.pl 2 dm12 84 1.20 0 dm_mu_kernel_smooth_p084 beta
+# view_update.pl 2 dm13 82 1.00 0 dm_mu_kernel_smooth_p082 beta_window
+# view_update.pl 2 dm13 74 2.80 0 dm_mu_kernel_smooth_p074 beta_window
+# view_update.pl 2 dm13 78 2.40 0 dm_mu_kernel_smooth_p074 beta_window
#
# view_update.pl 2 dm00 1 0.05 6 mu_kernel_smooth_dm beta_cg_smooth_06km
# view_update.pl 1 dm00 1 0.05 6 kappa_kernel_smooth_dm bulk_cg_smooth_06km
@@ -148,8 +152,8 @@
@dlayers = (-0.25,0.001,2,4,6,8,10,15,20,25,30,35,40);
$Nz = @dlayers;
-$imin = 1; $imax = $Nz; # default
-#$imin = 4; $imax = $imin; # testing
+#$imin = 1; $imax = $Nz; # default
+$imin = 2; $imax = $imin; # testing
for ($i = $imin; $i <= $imax; $i++) {
#$zcen = $utm_zmax - ($i-1)*$dzinc;
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update_local.tcl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update_local.tcl 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/UTILS/vtk/view_update_local.tcl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -40,7 +40,7 @@
# === load the unstructured grid data ===
vtkXMLUnstructuredGridReader kReader1
- kReader1 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/dm11/beta/dm_mu_kernel_smooth_p040.vtu
+ kReader1 SetFileName /net/sierra/raid1/carltape/socal/socal_3D/RUNS/MODELS/dm13/beta_window/dm_mu_kernel_smooth_p074.vtu
# vtkXMLUnstructuredGridReader kReader2
# kReader2 SetFileName vtu_files/lin_model/vp_low_2.vtu
@@ -67,7 +67,7 @@
# data mappers
vtkDataSetMapper kMapper1
kMapper1 SetInput [kReader1 GetOutput]
- kMapper1 SetScalarRange -1.6000e+00 1.6000e+00
+ kMapper1 SetScalarRange -2.4000e+00 2.4000e+00
# vtkDataSetMapper kMapper2
# kMapper2 SetInput [kReader2 GetOutput]
# kMapper2 SetScalarRange 4687.47 5729.13
@@ -85,7 +85,7 @@
# === generate plane cut ====
vtkPlane hrPlane1
- hrPlane1 SetOrigin 385719.8 3823329.2 -40000.0
+ hrPlane1 SetOrigin 385719.8 3823329.2 -1.0
hrPlane1 SetNormal 0 0 1
vtkCutter hrCut1
hrCut1 SetInput [kReader1 GetOutput]
@@ -93,7 +93,7 @@
vtkDataSetMapper hrMapper1
hrMapper1 SetInput [hrCut1 GetOutput]
hrMapper1 InterpolateScalarsBeforeMappingOn
- hrMapper1 SetScalarRange -1.6000e+00 1.6000e+00
+ hrMapper1 SetScalarRange -2.4000e+00 2.4000e+00
hrMapper1 SetLookupTable lut
vtkActor hrActor1
hrActor1 SetMapper hrMapper1
@@ -267,7 +267,7 @@
#=== text title ===
vtkTextActor titleActor
titleActor SetDisplayPosition 280 60
- titleActor SetInput "Subspace update for SHEAR-WAVE-SPEED -- model 40 -- Cut at z = 40.0 km"
+ titleActor SetInput "Subspace update for SHEAR-WAVE-SPEED -- model 78 -- Cut at z = 0.0 km"
set tprop [titleActor GetTextProperty]
$tprop SetJustificationToCentered
$tprop SetColor 0 0 1
@@ -284,8 +284,8 @@
iren SetRenderWindow renWin
vtkCamera cam1
- cam1 SetPosition 385719.8 3823329.2 10000.0
- cam1 SetFocalPoint 385719.8 3823329.2 -40000.0
+ cam1 SetPosition 385719.8 3823329.2 49999.0
+ cam1 SetFocalPoint 385719.8 3823329.2 -1.0
cam1 SetViewUp 0 1 0
#=== scene ===============
@@ -315,7 +315,7 @@
w2i SetInput renWin
vtkPostScriptWriter writer
writer SetInputConnection [w2i GetOutputPort]
- writer SetFileName "dm_mu_kernel_smooth_p040_13.ps"
+ writer SetFileName "dm_mu_kernel_smooth_p074_02.ps"
writer Write
wm withdraw .
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2008-11-24 00:05:12 UTC (rev 13377)
@@ -26,8 +26,8 @@
dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
% min and max periods for the different bandpassed datasets
-Tminvec = [2 6]; Tmaxvec = [30 30];
-%Tminvec = [2 3 6]; Tmaxvec = [30 30 30];
+%Tminvec = [2 6]; Tmaxvec = [30 30];
+Tminvec = [2 3 6]; Tmaxvec = [30 30 30];
comps = {'BHZ','BHR','BHT'};
% regardless of the component label for the DATA, the measurement code
@@ -38,7 +38,8 @@
% strings for labeling
sTbp = repmat(cellstr(' '),1,nper);
for tt = 1:nper
- sTbp{tt} = sprintf('[%is,%is]',Tminvec(tt),Tmaxvec(tt));
+ %sTbp{tt} = sprintf('[%is,%is]',Tminvec(tt),Tmaxvec(tt));
+ sTbp{tt} = sprintf('%i-%is',Tminvec(tt),Tmaxvec(tt));
end
%-------------------------
@@ -63,11 +64,11 @@
iwrite = input(' Enter 1 to write files (0 otherwise) : ');
% files: event IDs, CMTSOLUTIONS, stations
-dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_10/';
+dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_11/';
eid_file = [dir_source 'EIDs_only_loc'];
%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_' stmod];
eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
-cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v10'];
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v11'];
stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
%-------------------------
@@ -170,6 +171,7 @@
Tmin = Tminvec(tt);
Tmax = Tmaxvec(tt);
Ttag = sprintf('T%3.3i_T%3.3i',Tmin,Tmax);
+ Ttag2 = [Ttag '_' stmod]
%Ttag = sprintf('T%3.3i',Tmin); % old version
disp('-----------------------------------------');
@@ -177,7 +179,7 @@
disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_' Ttag '/'];
- mfile = [eids{ii} '_' Ttag '_window_chi'];
+ mfile = [eids{ii} '_' Ttag2 '_window_chi'];
filename = [dir1 mfile];
if ~exist(filename,'file')
@@ -323,8 +325,8 @@
end
if iwrite == 1
- sTfmti = repmat('%10i',1,nper);
- sTfmts = repmat('%10s',1,nper);
+ sTfmti = repmat('%8i',1,nper);
+ sTfmts = repmat('%8s',1,nper);
sTdash = repmat(cellstr('--'),1,nper);
stks = {'event','network','receiver','component'};
@@ -412,9 +414,9 @@
if iwrite == 1
% display sorted from greatest to least norm
- data_norms = [eids_num Ns nrec_all_event dnorm_sq*nevent dnorm_sq dnorm ws];
+ data_norms = [eids_num [1:nevent]' Ns nrec_all_event dnorm_sq*nevent dnorm_sq dnorm ws];
- % list only records that fit criteria
+ % This is useful if you want to only list events matching a certain set of criteria.
if 0==1
eid_reject = load('/net/sierra/raid1/carltape/results/KERNELS/kernels_exclude');
X1 = 20;
@@ -427,17 +429,17 @@
end
klabs = {'dnorm','nrec'};
- kind = -([4 3]+nper); % columns to sort
+ kind = -([5 4]+nper); % columns to sort
for kk = 1:2
dsort = sortrows(data_norms,kind(kk));
fid = fopen([odir ftag '_data_norms_sort_by_' klabs{kk} '.txt'],'w');
- fprintf(fid,['%12s%8s' sTfmts '%10s%10s%10s%10s%10s\n'],'eid','Nwin',...
+ fprintf(fid,['%12s%4s%8s' sTfmts '%8s%10s%10s%10s%10s\n'],'eid','ind','Nwin',...
sTbp{:},'Nrec','dnorm2*E','dnorm2','dnorm','weight');
- fprintf(fid,['%12s%8i' sTfmts '%10s%10.2f%10.4f%10s%10s\n'],...
- 'TOTAL',sum(Ns),sTdash{:},'--',sum(dnorm_sq*nevent),sum(dnorm_sq),'--','--');
+ fprintf(fid,['%12s%4i%8i' sTfmts '%8s%10.2f%10.4f%10s%10s\n'],...
+ 'TOTAL',' ',sum(Ns),sTdash{:},'--',sum(dnorm_sq*nevent),sum(dnorm_sq),'--','--');
for ii = 1:nevent
if any(dsort(ii,1)==eids_save)
- fprintf(fid,['%12i%8i' sTfmti '%10i%10.4f%10.4f%10.4f%10.4f\n'],dsort(ii,:));
+ fprintf(fid,['%12i%4i%8i' sTfmti '%8i%10.4f%10.4f%10.4f%10.4f\n'],dsort(ii,:));
end
end
%fprintf(fid,'%12s%8i%6s%6s%6s%12.4f%12s%12s\n','TOTAL',sum(Ns),'--','--','--',sum(dnorm_sq),'--','--');
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2008-11-23 22:24:02 UTC (rev 13376)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2008-11-24 00:05:12 UTC (rev 13377)
@@ -276,7 +276,8 @@
%lamvec = 10.^lampwr * sqrt(sum(diag(H'*H))); % based on trace of H'*H
else
%minlampwr = -2; maxlampwr = 5; % m03
- minlampwr = -3; maxlampwr = 3; % m04
+ minlampwr = -3; maxlampwr = 3; % m04 - m11
+ minlampwr = -5; maxlampwr = 0; % m12
lampwr = linspace(minlampwr,maxlampwr,numlam);
lamvec = 10.^lampwr;
end
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/9983429_T006_T030_LDF_CI_m13_m00_seis.pdf
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/9983429_T006_T030_LDF_CI_m13_m00_seis.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/README 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,7 @@
+Carl Tape, 20-Oct-2008
+
+/svn/cig/seismo/3D/ADJOINT_TOMO/iterate_adj_work/matlab/MISFIT_in_development/README
+
+These scripts are in development. The objective is to have a way to combine synthetics and measurements from many different models into a series of plots.
+
+---------
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/compare_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/compare_misfit.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/compare_misfit.m 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,302 @@
+%
+% compare_misfit.m
+% CARL TAPE, 23-July-2008
+% printed xxx
+%
+% This script reads in a set of window_chi output measurement files from
+% mt_measure_adj.f90, and it tabulates and plots misfit function values.
+%
+% /cig/seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/
+%
+% Copied from compute_misfit.m on 23-July-2008
+%
+% calls read_window_chi.m, plot_histo.m, readCMT.m, read_station_SPECFEM.m
+% called by xxx
+%
+
+%-------------------------
+
+clear
+close all
+format compact
+
+dir0 = '/net/sierra/raid1/carltape/socal/socal_3D/RUNS/';
+
+Tvec = [6 2];
+ % LOWEST period used for each measurement band-pass (upper is 40 s)
+comps = {'BHZ','BHR','BHT'};
+ % regardless of the component label for the DATA, the measurement code
+ % defaults so that the first two letters are always BH
+ncomp = length(comps);
+nper = length(Tvec);
+
+%DT = 0.011; % really this is only needed because it was left out from
+ % integrating the waveforms in mt-measure_adj.f90
+
+%-------------------------
+% USER INPUT
+
+iwrite = input(' Enter 1 to write files (0 otherwise) : ');
+
+% files: event IDs, CMTSOLUTIONS, stations
+dir_source = '/net/sierra/raid1/carltape/results/SOURCES/socal_09/';
+%eid_file = [dir_source 'EIDs_only_loc'];
+%eid_file = ['/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_' stmod];
+eid_file = '/net/sierra/raid1/carltape/results/EID_LISTS/kernels_use_m05';
+cmt_file_all = [dir_source 'SOCAL_FINAL_CMT_v09'];
+stations_file = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
+
+%-------------------------
+% read in list of event IDs and sources
+
+% load the event IDs corresponding to the kernels
+% load these as strings, since event IDs could have letters
+eids = textread(eid_file,'%s');
+eids_num = textread(eid_file,'%f'); % numeric version
+%eids_num = str2num(char(eids));
+nevent = length(eids);
+
+% read in CMT solutions
+[date,tshift,hdur,slat,slon,dep,M,eid,elabel] = readCMT(cmt_file_all, 13, 1);
+
+%-----------------------
+% read in stations file (used for synthetics)
+
+[rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(stations_file);
+nrec = length(stnm);
+for ii = 1:nrec, strec{ii} = [stnm{ii} '.' netwk{ii}]; end
+strec = strec(:);
+
+% make list of networks
+stnet = unique(netwk);
+nnet = length(stnet);
+
+%-----------------------
+
+iread = input(' Enter 1 to read in measurement files (0 otherwise) : ');
+
+imodel_min = input(' Enter minimum model index (0-5) : ');
+imodel_max = input(' Enter maximum model index (1-5) : ');
+nmodel = imodel_max - imodel_min;
+
+wdiff_array = zeros(nevent,nrec,ncomp,nper,nmodel);
+
+if iread == 1
+
+ % range of events
+ imin = 1; imax = nevent; % default
+ %imin = 4; imax = 7;
+ %imin = 4; imax = imin;
+
+ for imod = imodel_min:imodel_max
+
+ stmod = sprintf('m%2.2i',imod);
+ disp('====================================');
+ disp(['MODEL ' stmod]);
+
+ k2 = 0;
+ meas_array = [];
+ %for tt = 1:1
+ for tt = 1:nper
+ Tper = Tvec(tt);
+ sTper = sprintf('%3.3i',Tper);
+ if imod <= 2, sTper2 = sprintf('%2.2i',Tper); else sTper2 = sTper; end % CHT
+ disp('-----------------------------------------');
+
+ for ii = imin:imax % loop over events
+
+ disp(sprintf('Event %i out of %i -- %s',ii,nevent,eids{ii}));
+ dir1 = [dir0 eids{ii} '/' stmod '/MEASURE_T' sTper '/'];
+ mfile = [eids{ii} '_T' sTper2 '_window_chi'];
+ filename = [dir1 mfile];
+
+ if ~exist(filename,'file')
+ disp('--> file does not exist (or nwin = 0)');
+ else
+ % read the window_chi file
+ [stnet0,strec0,comp,iwin,iker,t1,t2,...
+ chiMT_dT,chiMT_dA,chiCC_dT,chiCC_dA,...
+ measMT_dT,measMT_dA,measCC_dT,measCC_dA,...
+ sigmaMT_dT,sigmaMT_dA,sigmaCC_dT,sigmaCC_dA,...
+ wind2,wins2,win_chi0,windur, recd2,recs2,rec_chi0,recdur,...
+ tr_chi,am_chi] = read_window_chi(filename);
+
+ % waveform differences normalized by data-squared
+ win_chi = win_chi0 ./ wind2;
+ rec_chi = rec_chi0 ./ recd2;
+
+ nwin = length(strec0);
+ disp(['--> nwin = ' num2str(nwin)]);
+
+ % assign the string network an integer index
+ index_net = NaN * ones(nwin,1);
+ for inet = 1:nnet
+ ind_net = find( strcmp(stnet(inet), stnet0) == 1 );
+ index_net(ind_net) = inet;
+ end
+
+ % assign the string receiver an integer index
+ index_rec = NaN * ones(nwin,1);
+ for irec = 1:nrec
+ ind_rec = find( strcmp(strec(irec), strec0) == 1 );
+ index_rec(ind_rec) = irec;
+ end
+
+ % assign the string component an integer index
+ index_comp = NaN * ones(nwin,1);
+ for icomp = 1:ncomp
+ ind_comp = find( strcmp(comps(icomp), comp) == 1 );
+ index_comp(ind_comp) = icomp;
+ end
+
+ % measurement indices
+ k1 = k2+1;
+ k2 = k1+nwin-1;
+ kinds = [k1:k2]';
+
+ meas_array(kinds,:) = [kinds tt*ones(nwin,1) ii*ones(nwin,1) index_net index_rec index_comp iwin ...
+ iker measCC_dT sigmaCC_dT measCC_dA sigmaCC_dA measMT_dT sigmaMT_dT win_chi rec_chi tr_chi];
+ % 1 kinds
+ % 2 index_T
+ % 3 index_event
+ % 4 index_network
+ % 5 index_rec
+ % 6 index_comp
+ % 7 iwin
+ % 8 iker
+ % 9 measCC_dT
+ % 10 sigmaCC_dT
+ % 11 measCC_dA
+ % 12 sigmaCC_dA
+ % 13 measMT_dT
+ % 14 sigmaMT_dT
+ % 15 win_chi
+ % 16 rec_chi
+ % 17 tr_chi
+ end
+
+ end
+ end
+
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+ for kk = 1:ncomp
+ imatch = find( and( and( meas_array(:,2)==tt, meas_array(:,3)==ii), ...
+ and( meas_array(:,5)==jj, meas_array(:,6)==kk) ));
+ if ~isempty(imatch)
+ % take the first only
+ wdiff_array(ii,jj,kk,tt,imod+1) = meas_array(imatch(1),16);
+ end
+ end
+ end
+ end
+ end
+
+ end
+
+ save('wdiff_array.mat','wdiff_array');
+
+else
+ load('wdiff_array.mat');
+end
+
+whos wdiff_array
+
+% total number of windows
+N = length(meas_array);
+
+%----------------------------------------------
+% CODE IN PROGRESS
+% for each record for each model (m0, m1) that has at least one window picked,
+% we save the integrated waveform difference as the purest measure of misfit
+% between the data and synthetics
+
+[m1,m2,m3,m4,m5] = size(wdiff_array)
+
+if 0==1
+
+ ratio_array = []; ll = 0;
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+ for kk = 1:ncomp
+ % check that there is a measurement for all models
+ if prod(wdiff_array(ii,jj,kk,tt,:)) > 0
+ ll = ll+1;
+ rat = wdiff_array(ii,jj,kk,tt,imodel_max+1) / wdiff_array(ii,jj,kk,tt,imodel_min+1);
+ ratio_array(ll,:) = [tt ii jj kk rat];
+ end
+ end
+ end
+ end
+ end
+
+ %ratio_sort = sortrows(ratio_array,5);
+ ratio_sort = sortrows(ratio_array,[-1 5]);
+
+ hmax = length(ratio_sort);
+ hmax = 100;
+ for hh = 1:hmax
+ tt = ratio_sort(hh,1);
+ ii = ratio_sort(hh,2);
+ jj = ratio_sort(hh,3);
+ kk = ratio_sort(hh,4);
+ disp(sprintf('%10s --> %6.1f %7s %5s %6.3f',char(eids{ii}),...
+ Tvec(tt),char(strec{jj}),char(comps{kk}),ratio_sort(hh,5) ));
+ end
+
+ %------------------------
+
+ rec_array = []; ll = 0;
+ for tt = 1:nper
+ for ii = imin:imax
+ disp(sprintf('Event %i out of %i',ii,nevent));
+ for jj = 1:nrec
+
+ % check that there is a measurement for all components for all models
+ if prod( wdiff_array(ii,jj,:,tt,:) ) > 0
+ ll = ll+1;
+ rat = prod( wdiff_array(ii,jj,:,tt,imodel_max+1) ./ wdiff_array(ii,jj,:,tt,imodel_min+1) );
+ rec_array(ll,:) = [tt ii jj rat];
+ end
+ end
+
+ end
+ end
+ rec_sort = sortrows(rec_array,[-1 4]);
+
+ for hh = 1:length(rec_sort)
+ tt = rec_sort(hh,1);
+ ii = rec_sort(hh,2);
+ jj = rec_sort(hh,3);
+ disp(sprintf('%10s --> %6.1f %7s %5s %10.6f',char(eids{ii}),...
+ Tvec(tt),char(strec{jj}),rec_sort(hh,4) ));
+ end
+
+ %------------------------
+
+% % loop over all stations
+% tarray = []; ll = 0;
+% for jj = 1:nrec
+% imatch = find( ratio_array(:,3)==jj);
+% if length(imatch)==6
+% ll = ll+1;
+% tarray(ll,:) = [jj prod(ratio_array(imatch,5)) ];
+% end
+% end
+% tarray_sort = sortrows(tarray,2);
+% for hh = 1:20
+% jj = tarray_sort(hh,1);
+% disp(sprintf('%7s %6.3f',char(strec{jj}),tarray_sort(hh,2) ));
+% end
+
+ figure; N = plot_histo(ratio_array(:,end),[0:0.2:5]);
+ figure; N = plot_histo(log(ratio_array(:,end)),[-4:0.2:4]);
+
+ break
+end
+
+%=======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt.m 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,186 @@
+%
+% misfit_gmt.m
+% CARL TAPE, 03-Sept-2008
+% printed xxx
+%
+%
+% calls xxx
+% called by xxx
+%
+
+%-------------------------
+
+function norms = misfit_gmt(odir,Tmin,Tmax,smodel1,smodel2,eid,stnm,netwk)
+
+disp('running misfit_gmt.m...');
+
+
+%-------------------------
+
+ifigure = 0;
+
+Ttag = sprintf('T%3.3i_T%3.3i',Tmin,Tmax);
+ssuffix1 = ['semd.sac.d.' Ttag '.' smodel1];
+ssuffix2 = ['semd.sac.d.' Ttag '.' smodel2];
+dsuffix = ['sac.d.' Ttag];
+
+comps = {'Z','R','T'};
+%chans = {'HL','BL','HN','BN','HH','BH'};
+chans = {'HH','BH'};
+ncomp = length(comps);
+nchan = length(chans);
+
+%dir0 = pwd;
+%odir = [dir0 '/OUTPUT_MISFIT/EVENTS/' eid '/' stnm '.' netwk '/'];
+
+%-------------------------
+
+fhits = zeros(ncomp,3);
+ffiles = repmat(cellstr(''),ncomp,3);
+for icomp = 1:ncomp
+ for ichan = 1:nchan
+
+ comp = [chans{ichan} comps{icomp}];
+ dfile = [odir eid '.' netwk '.' stnm '.' comp '.' dsuffix];
+ sfile1 = [odir stnm '.' netwk '.' comp '.' ssuffix1];
+ sfile2 = [odir stnm '.' netwk '.' comp '.' ssuffix2];
+
+ if exist(dfile,'file')
+ fhits(icomp,1) = fhits(icomp,1)+1;
+ ffiles{icomp,1} = dfile;
+ end
+ if exist(sfile1,'file')
+ fhits(icomp,2) = fhits(icomp,2)+1;
+ ffiles{icomp,2} = sfile1;
+ end
+ if exist(sfile2,'file')
+ fhits(icomp,3) = fhits(icomp,3)+1;
+ ffiles{icomp,3} = sfile2;
+ end
+
+ end
+end
+fhits_sum = sum(fhits,2);
+
+if ifigure==1
+ figure; nr=5; nc=3;
+ st1 = ['syn(' smodel1 ')'];
+ st2 = ['syn(' smodel2 ')'];
+end
+
+norms = zeros(9,ncomp);
+
+for icomp = 1:ncomp
+%for icomp = 3
+
+ if fhits_sum(icomp) == 3
+ dfile = ffiles{icomp,1};
+ sfile1 = ffiles{icomp,2};
+ sfile2 = ffiles{icomp,3};
+
+ [seisD,HdrData,tnu,pobj,timeD] = readsac(dfile,0,'l');
+ [seisS1,HdrData,tnu,pobj,timeS1] = readsac(sfile1,0,'l');
+ [seisS2,HdrData,tnu,pobj,timeS2] = readsac(sfile2,0,'l');
+
+ ti = timeD(:);
+ n = length(ti);
+
+ % interpolate synthetics exactly onto the data time vector
+ % (In theory, they should already be within two time-steps at this point.)
+ s1 = interp1(timeS1,seisS1,ti,'linear','extrap');
+ s2 = interp1(timeS2,seisS2,ti,'linear','extrap');
+
+ % residuals for the two models; also perturbation to seismogram
+ res1 = s1 - seisD;
+ res2 = s2 - seisD;
+ spert = s2 - s1;
+
+ ymax = 1.05*max(abs([seisD' s1' s2' res1' res2']));
+ tmin = min(ti); tmax = max(ti);
+ ax0 = [tmin tmax -ymax ymax];
+
+ % compute nine different norms -- IS THIS EXACTLY WHAT WE WANT?
+ norms(1,icomp) = norm(seisD);
+ norms(2,icomp) = norm(s1);
+ norms(3,icomp) = norm(s2);
+ norms(4,icomp) = norm(res1);
+ norms(5,icomp) = norm(res2);
+ norms(6,icomp) = norm(s2-s1);
+ norms(7,icomp) = norm(res1)^2 / norm(seisD)^2; % chi^2
+ norms(8,icomp) = norm(res2)^2 / norm(seisD)^2; % chi^2
+ norms(9,icomp) = 100*(1 - (norm(res2)^2 / norm(res1)^2)); % VR
+
+ if ifigure==1
+ subplot(nr,nc,icomp);
+ plot(ti,s1,'r',ti,seisD,'b'); legend(st1,'data'); axis(ax0);
+ subplot(nr,nc,nc*1+icomp);
+ plot(ti,res1,'k'); legend([st1 ' - data']); axis(ax0);
+ subplot(nr,nc,nc*2+icomp);
+ plot(ti,s2,'r',ti,seisD,'b'); legend(st2,'data'); axis(ax0);
+ subplot(nr,nc,nc*3+icomp);
+ plot(ti,res2,'k'); legend([st2 ' - data']); axis(ax0);
+ subplot(nr,nc,nc*4+icomp);
+ plot(ti,spert,'k'); legend([st2 ' - ' st1]); axis(ax0);
+ end
+
+ % WRITE FILES FOR GMT PLOTTING
+ filename = [odir 'GMT_time_series_' Ttag '_' comps{icomp} '.dat'];
+ fid = fopen(filename,'w');
+ for ii = 1:n
+ fprintf(fid,[repmat('%12.4e',1,7) '\n'],...
+ ti(ii),seisD(ii),s1(ii),s2(ii),res1(ii),res2(ii),spert(ii) );
+ end
+ fclose(fid);
+
+ filename = [odir 'GMT_axes_' Ttag '_' comps{icomp} '.dat'];
+ fid = fopen(filename,'w');
+ fprintf(fid,[repmat('%12.4e',1,4) '\n'],...
+ ax0(1),ax0(2),ax0(3),ax0(4) );
+ fclose(fid);
+
+ filename = [odir 'GMT_norms_' Ttag '_' comps{icomp} '.dat'];
+ fid = fopen(filename,'w');
+ for kk = 1:9, fprintf(fid,'%12.4e\n',norms(kk,icomp) ); end
+ fclose(fid);
+
+ end
+
+end
+
+if ifigure == 1
+ orient tall, wysiwyg, fontsize(6)
+end
+
+return
+
+%-----------------------------------
+% if all three components exist, then compute additional time series for plots
+
+if sum(fhits_sum) == 9
+
+ % compute the time vector to use for all components
+ mint = zeros(ncomp,3);
+ maxt = zeros(ncomp,3);
+ for icomp = 1:ncomp
+ dfile = ffiles{icomp,1};
+ sfile1 = ffiles{icomp,2};
+ sfile2 = ffiles{icomp,3};
+ [seisD,HdrData,tnu,pobj,timeD] = readsac(dfile,0,'l');
+ [seisS1,HdrData,tnu,pobj,timeS1] = readsac(sfile1,0,'l');
+ [seisS2,HdrData,tnu,pobj,timeS2] = readsac(sfile2,0,'l');
+
+ mint(icomp,1) = min(timeD);
+ mint(icomp,2) = min(timeS1);
+ mint(icomp,3) = min(timeS2);
+ maxt(icomp,1) = max(timeD);
+ maxt(icomp,2) = max(timeS1);
+ maxt(icomp,3) = max(timeS2);
+ end
+
+ tmin = max(mint(:));
+ tmax = min(maxt(:));
+ ti = [tmin:0.05:tmax]';
+
+end
+
+%=======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt_run.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt_run.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/misfit_gmt_run.m 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,72 @@
+%
+% misfit_gmt_run.m
+% CARL TAPE, 03-Sept-2008
+% printed xxx
+%
+%
+% calls xxx
+% called by xxx
+%
+
+%-------------------------
+
+clear
+close all
+format compact
+
+Tmin = 6;
+Tmax = 40;
+smodel1 = 'm00';
+smodel2 = 'm07';
+%stnm = 'USC'; netwk = 'CI';
+
+ncomp = 3;
+nper = 1;
+nmodel = 2;
+
+dir0 = pwd;
+dir1 = [dir0 '/OUTPUT_MISFIT/EVENTS/'];
+
+stations_file = [dir1 'STATIONS'];
+
+%-------------------------
+% read in list of event IDs and sources
+
+nevent = 1;
+eid = '14383980';
+dir2 = [dir1 eid '/'];
+
+%-----------------------
+% read in stations file (used for synthetics)
+
+[rlon0,rlat0,relev0,rburial0,stnm0,netwk0] = read_station_SPECFEM(stations_file);
+nrec = length(stnm0);
+%for ii = 1:nrec, strec{ii} = [stnm{ii} '.' netwk{ii}]; end
+%strec = strec(:);
+
+for irec = 1:nrec
+ disp(sprintf('%i out of %i -- %s.%s',irec,nrec,stnm0{irec},netwk0{irec}));
+end
+
+%-------------------------
+
+nevent = 1; ievent = 1;
+VR_array = zeros(nevent,nrec,ncomp,nper);
+
+for irec = 1:nrec
+
+ stnm = stnm0{irec};
+ netwk = netwk0{irec};
+ disp(sprintf('%i out of %i -- %s.%s',irec,nrec,stnm,netwk));
+
+ odir = [dir2 stnm '.' netwk '/'];
+
+ norms = misfit_gmt(odir,Tmin,Tmax,smodel1,smodel2,eid,stnm,netwk);
+
+ VR_array(ievent,irec,1,1) = norms(end,1); % Z
+ VR_array(ievent,irec,2,1) = norms(end,2); % R
+ VR_array(ievent,irec,3,1) = norms(end,3); % T
+
+end
+
+%=======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/plot_misfit.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/plot_misfit.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/plot_misfit.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,349 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 02-Sept-2008
+# plot_misfit.pl
+#
+#
+# EXAMPLE:
+# plot_misfit.pl 9818433 USC/CI 6/40
+#
+#-----------------------------------
+
+#use Getopt::Std;
+#use POSIX;
+
+if (@ARGV < 3) {die("Usage: plot_misfit.pl eid sta/net Tmin/Tmax\n");}
+($eid,$stanet,$Ts) = @ARGV;
+
+#($sta,$net) = split(/\//,$stanet);
+($sta,$net) = split("/",$stanet);
+
+# bandpass range
+($Tmin,$Tmax) = split("/",$Ts);
+$sTmin = sprintf("T%3.3i",$Tmin);
+$sTmax = sprintf("T%3.3i",$Tmax);
+$Ttag = "${sTmin}_${sTmax}";
+
+#---------------------------
+# USER INPUT
+
+# models to compare
+$minm = 0;
+$maxm = 7;
+$smodel_min = sprintf("m%2.2i",$minm);
+$smodel_max = sprintf("m%2.2i",$maxm);
+
+# base directory (contains STATIONS and eid_list)
+#$dir0 = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/iterate_adj_work/matlab/OUTPUT_MISFIT";
+
+# CMTSOLUTION file
+$cmtfile = "../CMTSOLUTION";
+ at cmtvals = &get_cmt($cmtfile); # see subroutine below
+
+# STATIONS file
+$stafile = "../STATIONS_ADJOINT";
+
+# plot faults (southern California)
+$ifault = 1;
+
+# file to use to extract the SAC headers
+$Zsyn = "${sta}.${net}.BHZ.semd.sac.d.${Ttag}.${smodel_max}";
+
+$saclst = "/opt/seismo-util/bin/saclst";
+
+# axes positions
+$x0_title = 0.5; $y0_title = 10.5;
+$dY_under_title = 0.4;
+$dY_under_subtitle = 0.4;
+$dX_between_col = 0.2;
+$dY_under_seis = 0.3;
+
+# dimensions of seismogram subplots
+$xwid = 2.4;
+$ywid = 1.0;
+$J = "-JX${xwid}/${ywid}";
+
+# subtitles
+ at subs = ("(a) Vertical component, Z","(b) Radial component, R","(c) Transverse component, T","(d) Variance reduction for each component","(e) Source-receiver geometry");
+
+ at slabs = ("d and s($smodel_min)","d and s($smodel_max)","s($smodel_min) - d","s($smodel_max) - d","s($smodel_max) - s($smodel_min)");
+$nrows = 5; # <= 5
+
+$x0_map = 5.0; $y0_map = 2.3 + (5-$nrows)*$ywid;
+$x0_misfit = 1.5; $y0_misfit = 2.9 + (5-$nrows)*$ywid;
+
+#---------------------------
+
+# create psmeca file from CMTSOLUTION file
+ at M[0..5] = @cmtvals[12..17];
+$year = $cmtvals[0]; $month = $cmtvals[1]; $day = $cmtvals[2];
+$elat = $cmtvals[9]; $elon = $cmtvals[10]; $edep = $cmtvals[11];
+$eid = $cmtvals[18];
+$cmtpsmeca = "temp.txt";
+open(MIJ,">$cmtpsmeca");
+#printf MIJ "$elon $elat $edep $M[0] $M[1] $M[2] $M[3] $M[4] $M[5] 1 X Y $eid";
+printf MIJ "$elon $elat $edep $M[0] $M[1] $M[2] $M[3] $M[4] $M[5] 1";
+close(MIJ);
+
+# compute moment using Dahlen and Tromp, Eq. (5.91)
+# compute moment magnitude using Kanamori (1977) -- see Carl's Latex notes cmt.pdf
+$M0 = 1/sqrt(2) * sqrt($M[0]**2 + $M[1]**2 + $M[2]**2 + 2*($M[3]**2 + $M[4]**2 + $M[5]**2 ) );
+$k = -11.2 + (2/3)*log(5)/log(10);
+$Mw = (2/3)*log($M0)/log(10) + $k;
+$stM0 = sprintf("M0 = %.3e dyne-cm",$M0);
+$stMw = sprintf("Mw = %.2f",$Mw);
+print "\n moment : $stM0 \n moment magnitude : $stMw\n";
+
+# get the receiver location, azimuth, distance (from the $Zsyn file)
+# We want the azimuth for the file name, so that all the PDFs will be sorted by azimuth.
+(undef,$tmin0,$tmax0,$slon,$slat,$dist,$az,$Tchan) = split(" ",`$saclst b e stlo stla dist az kcmpnm f $Zsyn`);
+$edate = sprintf("%4.4i . %2.2i . %2.2i",$year,$month,$day);
+$stedep = sprintf("depth = %.1f km",$edep);
+$strec = "$sta . $net";
+$stdist = sprintf("dist = %.1f km",$dist);
+$staz = sprintf("az = %.1f deg",$az);
+
+#---------------------------
+
+$fontno = 1; # font number to use for text
+
+# set GMT defaults
+`gmtset BASEMAP_TYPE plain PAPER_MEDIA letter ANNOT_FONT_SIZE = 10p LABEL_FONT_SIZE = 10p PAGE_ORIENTATION = portrait PLOT_DEGREE_FORMAT D TICK_LENGTH 4p`;
+
+# name of the output files
+$azlabel = sprintf("%3.3i",$az);
+$name = "${azlabel}_${sta}_${net}_${Ttag}";
+$psfile = "$name.ps";
+$jpg_file = "$name.jpg";
+
+# overall title
+$title = "Event $eid, station ${sta}.${net}, bandpass ${Tmin}-${Tmax} s, models ${smodel_min} and ${smodel_max}";
+`pstext -R0/1/0/1 -JX1 -N -Xa${x0_title} -Ya${y0_title} -K -P > $psfile <<EOF\n 0 0 14 0 1 LM $title \nEOF\n`; # START (no -O)
+
+# subtitles
+$x0_sub = $x0_title;
+$y0_sub = $y0_title - $dY_under_title;
+
+for ($k = 1; $k <= 3; $k = $k+1) {
+ $x_sub = $x0_sub + ($k-1)*($xwid+$dX_between_col);
+ $origin_sub = "-Xa${x_sub} -Ya${y0_sub}";
+ `pstext -R0/1/0/1 -JX1 -N $origin_sub -K -O >> $psfile <<EOF\n 0 0 10 0 1 LM $subs[$k-1] \nEOF\n`;
+}
+
+# subtitle for misfit plot
+$x0_sub2 = $x0_misfit - 0.8;
+$y0_sub2 = $y0_title - $dY_under_title - $dY_under_subtitle - $nrows*$ywid - $dY_under_seis;
+$origin_sub = "-Xa${x0_sub2} -Ya${y0_sub2}";
+`pstext -R0/1/0/1 -JX1 -N $origin_sub -K -O >> $psfile <<EOF\n 0 0 10 0 1 LM $subs[3] \nEOF\n`;
+
+# subtitle for map
+$x0_sub3 = $x0_map - 0.8;
+$y0_sub3 = $y0_sub2;
+$origin_sub = "-Xa${x0_sub3} -Ya${y0_sub3}";
+`pstext -R0/1/0/1 -JX1 -N $origin_sub -K -O >> $psfile <<EOF\n 0 0 10 0 1 LM $subs[4] \nEOF\n`;
+
+#---------------------------------------------------------------------------
+# PLOT SEISMOGRAMS AND RESIDUALS
+
+$x0_seis = $x0_title;
+$y0_seis = $y0_title - $dY_under_title - $dY_under_subtitle - $ywid;
+
+$B0 = "-Ba50f10g50:\" \":/a1f1::wesN";
+$B = "-Ba50f10g50:\" \":/a1f1::wesn";
+
+#$slabinfo = "-C2p -W0/255/255o,0/255/255 -G0/0/0 -N";
+$slabinfo = "-C3p -W255/255/255o,0/0/0 -G0/0/0 -N";
+
+ at comps = ("Z","R","T");
+for ($k = 1; $k <= 3; $k = $k+1) {
+ $comp = $comps[$k-1];
+
+ $Pfile = "GMT_time_series_${Ttag}_${comp}.dat";
+ $bounds = "GMT_axes_${Ttag}_${comp}.dat";
+ $Nfile = "GMT_norms_${Ttag}_${comp}.dat";
+
+ # if the seismograms exist, then plot them
+ if ( (-f $Pfile) && (-f $bounds) ) {
+
+ open(IN,"$bounds"); @lines = <IN>;
+ ($tmin,$tmax,$ymin,$ymax) = split(" ",$lines[0]);
+ $tmin = -10;
+ if($Tmin==2) {$tmax = $tmax/2;} # USER: reduce the records for 2s
+ $R = "-R$tmin/$tmax/$ymin/$ymax";
+
+ # plot seismograms
+ $x_seis = $x0_seis + ($k-1)*($xwid+$dX_between_col);
+ $x_seis_lab = $x_seis + $xwid - 0.05;
+
+ for ($j = 1; $j <= $nrows; $j = $j+1) {
+ $y_seis = $y0_seis - ($j-1)*$ywid;
+ $originP = "-Xa${x_seis} -Ya${y_seis}";
+
+ $y_seis_lab = $y_seis + $ywid - 0.05;
+ $originPlab = "-Xa${x_seis_lab} -Ya${y_seis_lab}";
+
+ if($j==1) {`psbasemap $J $R $B0 -O -K $originP >> $psfile`}
+ else {`psbasemap $J $R $B -O -K $originP >> $psfile`}
+ if ($j==1) {
+ `awk \'{print \$1,\$2}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `awk \'{print \$1,\$3}\' $Pfile | psxy $J $R -W0.5p,255/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[0]\nEOF\n`;
+ } elsif ($j==2) {
+ `awk \'{print \$1,\$2}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `awk \'{print \$1,\$4}\' $Pfile | psxy $J $R -W0.5p,255/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[1]\nEOF\n`;
+ } elsif ($j==3) {
+ `awk \'{print \$1,\$5}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[2]\nEOF\n`;
+ } elsif ($j==4) {
+ `awk \'{print \$1,\$6}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[3]\nEOF\n`;
+ } elsif ($j==5) {
+ `awk \'{print \$1,\$7}\' $Pfile | psxy $J $R -W0.5p,0/0/0 -O -K $originP >> $psfile`;
+ `pstext -R0/1/0/1 -JX1 $slabinfo -K -O $originPlab >> $psfile <<EOF\n0 0 8 0 $fontno RT $slabs[4]\nEOF\n`;
+ }
+ }
+
+ } else {
+ print "seismogram files $Pfile and $bounds do not both exist\n";
+ }
+
+}
+
+#---------------------------------------------------------------------------
+# PLOT THE MISFIT
+
+$originN = "-Xa${x0_misfit} -Ya${y0_misfit}";
+
+$xwid = 1.5;
+$ywid = 1.3;
+$J = "-JX${xwid}/${ywid}";
+$B = "-Ba1f1:\"Model index\":/a50f10:\"Variance reduction\":WeSn";
+$kmax = $maxm + 3;
+$R = "-R-0.5/$kmax/-50/100";
+
+$Zinfo = "-Sc12p -W0.5p,0/0/0 -G255 -N";
+
+print "J option --> $J\n";
+print "R option --> $R\n";
+print "B option --> $B\n";
+
+ at comps = ("Z","R","T");
+for ($k = 1; $k <= 3; $k = $k+1) {
+ $comp = $comps[$k-1];
+
+ $Nfile = "GMT_norms_${Ttag}_${comp}.dat";
+
+ # if the norms file exists
+ if ( -f $Nfile ) {
+ # reduction in the misfit
+ open(IN,"$Nfile"); @Nlines = <IN>;
+ #$yred = 100*(1.0 - $Nlines[4]/$Nlines[3]);
+ $yred = $Nlines[8]; chomp($yred);
+ print "yred --> $yred\n";
+
+ # plot the misfit
+ if($k==1) {`psbasemap $J $R $B -O -K $originN >> $psfile`}
+ `psxy $R $J -W1p -K -O $originN >> $psfile <<EOF\n$minm 0\n$maxm $yred\nEOF\n`;
+ `psxy $R $J $Zinfo -K -O $originN >> $psfile <<EOF\n$minm 0\nEOF\n`;
+ `psxy $R $J $Zinfo -K -O $originN >> $psfile <<EOF\n$maxm $yred\nEOF\n`;
+ `pstext $R $J -N -K -O $originN >> $psfile <<EOF\n$maxm $yred 8 0 $fontno CM $comp\nEOF\n`;
+
+ } else {
+ print "norm file $Nfile does not exist\n";
+ }
+
+}
+
+#---------------------------------------------------------------------------
+# PLOT BASE MAP with stations, faults, CMT
+
+$originM = "-Xa${x0_map} -Ya${y0_map}";
+$proj = "-JM2.5";
+$textinfo = "-G0 -S1p,255/255/255";
+
+$xmin = -122; $xmax = -114;
+$ymin = 32; $ymax = 37;
+$bounds="-R$xmin/$xmax/$ymin/$ymax ";
+$xtick = 1; $ytick = 1;
+$tick="-B${xtick}/${ytick}::wesn";
+$cinfo="-W0.5 -Dh -A100";
+
+# plot the base map
+`psbasemap $bounds $proj $tick $originM -K -O >> $psfile`;
+`pscoast $bounds $proj $cinfo -S150/200/255 -W2 -G200/255/150 -N1/1p -N2/0.5p -O -K $originM >> $psfile`;
+
+# plot southern California faults
+if($ifault==1) {
+ $dir1 = "/home/carltape/gmt";
+ $fault_file = "$dir1/faults/jennings.xy";
+ $kcf_file = "$dir1/faults/kcf.xy";
+ $breck_file = "$dir1/faults/breck.xy";
+ $fault_infoK = "-M -W0.5p,0/0/0";
+ `psxy $fault_file $bounds $proj $fault_infoK -O -K $originM >> $psfile`;
+ `psxy $kcf_file $bounds $proj $fault_infoK -O -K $originM >> $psfile`;
+ `psxy $breck_file $bounds $proj $fault_infoK -O -K $originM >> $psfile`;
+}
+
+# plot the stations
+`awk \'{print \$4,\$3,0,4,B,C}\' $stafile | psxy $proj $bounds -M -St0.10 -W0.25p -G200/200/200 -N -O -K $originM >> $psfile`;
+
+#------------
+
+# plot labels next to the map
+ at mapstrings = ($edate,$eid,$stMw,$stedep,$strec,$stdist,$staz);
+$nstrings = @mapstrings;
+$x1 = $xmin - 4*$xtick;
+for ($i=0; $i < $nstrings; $i++) {
+ $y1 = $ymax - 0.5 - $i*$ytick*0.7;
+ `pstext $bounds $proj $tick -N -K -O $originM >> $psfile <<EOF\n $x1 $y1 10 0 $fontno LM $mapstrings[$i] \nEOF\n`;
+}
+
+# plot ray path, station, and CMT source
+`psxy $bounds $proj -W1p -K -O $originM >> $psfile <<EOF\n$elon $elat\n$slon $slat\nEOF\n`; # ray path
+#$cmtfont = 10; $cmt_scale = 0.3; $cmtinfo = "-Sm${cmt_scale}/$cmtfont -L0.5p/0/0/0 -G255/0/0";
+#`psmeca $cmtpsmeca $bounds $proj $cmtinfo -O -K $originM >> $psfile`; # CMT source
+$cmt_scale = 0.3; $cmtinfo = "-Sm${cmt_scale} -L0.5p/0/0/0 -G255/0/0";
+`psmeca $cmtpsmeca $bounds $proj $cmtinfo -O -K $originM >> $psfile`; # CMT source
+$elat_shift = $elat + 0.6;
+`pstext $bounds $proj $textinfo -K -O $originM >> $psfile <<EOF\n $elon $elat_shift 9 0 $fontno CM $eid\nEOF\n`; # station label
+`rm $cmtpsmeca`;
+
+`psxy $bounds $proj -St0.15 -W2 -G255/0/0 -K -O $originM >> $psfile <<EOF\n$slon $slat\nEOF\n`; # red station
+$slat_shift = $slat - 0.3;
+`pstext $bounds $proj $textinfo -K -O $originM >> $psfile <<EOF\n $slon $slat_shift 8 0 $fontno CM $sta \nEOF\n`; # station label
+
+#---------------------------------------------------------------------------
+
+`pstext -R0/1/0/1 -JX1 -N -Xa${x0_title} -Ya${y0_title} -O >> $psfile <<EOF\n 10 10 14 0 1 LM $title\nEOF\n`; # FINISH (no -K)
+
+#`convert $psfile $jpg_file`; `xv $jpg_file &`;
+`ps2pdf $psfile`;
+#`rm $psfile`; # remove PS file
+#system("xv $jpg_file &");
+
+#================================================
+sub get_cmt {
+ my ($cmtfile)=@_;
+ open(CMT,"$cmtfile") or die("Error opening cmtfile $cmtfile\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,$eid)=split(" ",$cmt[1]);
+ 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,$eid);
+}
+
+#================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/plot_misfit.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/seis_norm.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/seis_norm.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/seis_norm.m 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,24 @@
+%
+% function x = seis_norm(t,y)
+% Carl Tape, 03-Sept-2008
+%
+% This function computes the L2 norm of a seismogram.
+%
+% t2 2 N 2 N 2
+% INT s(t) dt dt SUM (s_i) SUM (s_i)
+% 2 t1 i=1 i=1
+% x = ------------------- = ---------------- = ---------
+% t2 - t1 dt*N N
+%
+% Note that the number will not depend on the length of the time series.
+%
+% calls xxx
+% called by xxx
+%
+
+function x = seis_norm(t,y)
+
+x = sqrt( sum(y.^2) / length(y) );
+
+%============================================================
+
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/setup_misfit_dir.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/setup_misfit_dir.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/setup_misfit_dir.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,304 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 03-Sept-2008
+# setup_misfit_dir.pl
+#
+# This script xxx
+#
+# EXAMPLE: setup_misfit_dir.pl m00 m07 6/40
+#
+#-----------------------------------
+
+if (@ARGV < 3) {die("Usage: setup_misfit_dir.pl smodel_min smodel_max Tmin/Tmax\n")}
+($minm,$maxm,$Ts) = @ARGV;
+
+$pwd = $ENV{PWD};
+print "\n$pwd\n";
+
+# labels for bandpass-filtered data and synthetics
+($Tmin,$Tmax) = split("/",$Ts);
+$sTmin = sprintf("T%3.3i",$Tmin);
+$sTmax = sprintf("T%3.3i",$Tmax);
+$Ttag = "${sTmin}_${sTmax}";
+
+#-------------------------------------------
+# USER INPUT
+
+# directory containing all CMTSOLUTION files
+$dir_src = "/net/sierra/raid1/carltape/results/SOURCES/socal_09";
+$dir_cmt = "${dir_src}/CMT_files_pre_inverted";
+if (not -e $dir_src) {die("check if dir_src $dir_src exist or not\n")}
+if (not -e $dir_cmt) {die("check if dir_cmt $dir_cmt exist or not\n")}
+
+# list of event IDs to use
+$file_eid0 = "/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_m07";
+if (not -f ${file_eid0}) {die("check if file_eid ${file_eid0} exist or not\n")}
+
+# FULL stations file
+$file_stations0 = "/net/denali/home2/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
+if (not -f ${file_stations0}) {die("check if file_stations ${file_stations0} exist or not\n")}
+
+# data and synthetics directories
+$dir_data0 = "/net/sierra/raid1/carltape/socal/socal_3D/DATA/FINAL";
+$dir_syn01 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${minm}";
+$dir_syn02 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${maxm}";
+if (not -e ${dir_data0}) {die("check if dir_data ${dir_data0} exist or not\n")}
+if (not -e ${dir_syn01}) {die("check if dir_syn ${dir_syn01} exist or not\n")}
+if (not -e ${dir_syn02}) {die("check if dir_syn ${dir_syn02} exist or not\n")}
+
+# directory containing the output for windows, measurements, adjoint sources, etc
+$dir_run = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS";
+if (not -e ${dir_run}) {die("check if dir_run ${dir_run} exist or not\n")}
+
+#-------------------------------------------
+
+# make the base output directory if it does not exist
+$dir0 = "OUTPUT_MISFIT";
+`mkdir -p $dir0`;
+$dir1 = "$dir0/EVENTS";
+`mkdir -p $dir1`;
+
+# link various files
+$file_eid = "$dir0/eid_list";
+`ln -s ${file_eid0} ${file_eid}`;
+$file_stations = "$dir0/STATIONS";
+`ln -s ${file_stations0} ${file_stations}`;
+
+# open EID list
+open(IN,"${file_eid}"); @eids = <IN>; close(IN);
+$nevent = @eids;
+print "\n $nevent events in the list";
+
+# open STATIONS file
+open(IN,"${file_stations}"); @slines = <IN>; close(IN);
+$nrec = @slines - 1;
+print "\n $nrec stations in the list";
+
+#-------------------------------------------
+
+if (0==1) {
+ for ($ievent = 1; $ievent <= $nevent; $ievent++) {
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ }
+ die("testing");
+}
+
+#=============================================
+# MAKE DIRECTORIES
+
+#$imin = 1; $imax = $nevent;
+#$imin = 1; $imax = 100;
+$imin = 142; $imax = $imin;
+
+$isetup1 = 0;
+if ($isetup1 == 1) {
+
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+
+ # make base directory for each event
+ $dir2 = "$dir1/$eid";
+ `mkdir -p $dir2`;
+
+ # link CMTSOLUTION file
+ $cmtfile0 = "${dir_cmt}/CMTSOLUTION_${eid}";
+ $cmtfile = "$dir2/CMTSOLUTION";
+ `ln -s $cmtfile0 $cmtfile`;
+
+ # link STATIONS file from the windowing/measurement codes (different for each period range)
+ # --> use the max model to show the stations used
+ $stafile0 = "${dir_run}/$eid/$maxm/MEASURE_${sTmin}/ADJOINT_${sTmin}/STATIONS_ADJOINT";
+ if (not -f $stafile0) {die("check if stafile0 $stafile0 exist or not\n")}
+ $stafile = "$dir2/STATIONS_ADJOINT";
+ `ln -s $stafile0 $stafile`;
+
+ # make a directory for each station in the STATIONS list
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $dir_rec = "$dir2/$sta.$net";
+ print "\n $irec out of $nrec -- $dir_rec";
+ `mkdir -p $dir_rec`;
+ }
+ }
+
+}
+
+#=============================================
+# COPY IN DATA AND SYNTHETICS
+# NOTE: For now we are copying in all available records,
+# but really we only want BHZ,BHR,BHT (or HHZ,HHR,HHT).
+
+$isetup2 = 0;
+
+if ($isetup2 == 1) {
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+
+ # data and synthetics directories for each event
+ $dir_data = "${dir_data0}/$eid/PROCESSED/PROCESSED_${Ttag}";
+ $dir_syn1 = "${dir_syn01}/$eid/PROCESSED/PROCESSED_${Ttag}";
+ $dir_syn2 = "${dir_syn02}/$eid/PROCESSED/PROCESSED_${Ttag}";
+ if (not -e ${dir_data}) {die("check if dir_data ${dir_data} exist or not\n");}
+ if (not -e ${dir_syn1}) {die("check if dir_syn ${dir_syn1} exist or not\n");}
+ if (not -e ${dir_syn2}) {die("check if dir_syn ${dir_syn2} exist or not\n");}
+
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ #for ($irec = 128; $irec <= 128; $irec++) {
+
+ cd_directory($pwd);
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ # link plotting script
+ $plotfile0 = "${pwd}/plot_misfit.pl";
+ $plotfile = "${dir_rec}/plot_misfit.pl";
+ `rm $plotfile`;
+ `ln -s $plotfile0 $plotfile`;
+
+ $dtags = "${dir_data}/*$net.$sta*"; @dfiles = glob($dtags); `cp @dfiles ${dir_rec}`;
+
+ $stags1 = "${dir_syn1}/$sta.$net*"; @sfiles1 = glob($stags1); #`cp @sfiles1 ${dir_rec}`;
+ cd_directory($dir_rec);
+ `cp @sfiles1 .`;
+ `mv_files.pl -x "*semd.sac.d.${Ttag}" "*semd.sac.d.${Ttag}.$minm"`;
+
+ $stags2 = "${dir_syn2}/$sta.$net*"; @sfiles2 = glob($stags2); #`cp @sfiles2 ${dir_rec}`;
+ `cp @sfiles2 .`;
+ `mv_files.pl -x "*semd.sac.d.${Ttag}" "*semd.sac.d.${Ttag}.$maxm"`;
+ }
+ }
+
+}
+
+#=============================================
+# NOW RUN MATLAB CODE misfit_gmt_run.m
+
+#=============================================
+# GENERATE GMT FIGURES
+
+$isetup3 = 0;
+
+if ($isetup3 == 1) {
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ #for ($irec = 104; $irec <= 104; $irec++) {
+
+ cd_directory($pwd);
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ cd_directory($dir_rec);
+ @datfiles = glob("*.dat");
+ $ndat = @datfiles;
+
+ if($ndat > 0) {
+ `plot_misfit.pl $eid $sta/$net $Ts`;
+ }
+ }
+ }
+
+}
+
+#=============================================
+# CONCATENATE PDF FIGURES
+
+$isetup4 = 0;
+
+if ($isetup4 == 1) {
+
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ @pdffiles = glob("${dir_rec}/*${Ttag}.pdf");
+ $nfile = @pdffiles;
+ #$pdffile = "${dir_rec}/${sta}_${net}_${Ttag}.pdf";
+ if($nfile == 1) {
+ $pdffile = $pdffiles[0]; chomp($pdffile);
+ `cp $pdffile $dir2`;
+ }
+ }
+
+ # concatenate into a single pdf
+ $outfile = "$dir1/${eid}_${Ttag}_waveform_az.pdf";
+ `rm $outfile`;
+ @pdfall = glob("${dir2}/*${Ttag}.pdf");
+ print "\n /home/carltape/bin/pdcat -r @pdfall $outfile\n";
+ `/home/carltape/bin/pdcat -r @pdfall $outfile`;
+ `rm @pdfall`;
+ }
+}
+
+#=============================================
+# DELETE FILES
+
+$iclean = 0;
+
+if ($iclean == 1) {
+
+ for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ $dir2 = "$dir1/$eid";
+
+ for ($irec = 1; $irec <= $nrec; $irec++) {
+ #for ($irec = 128; $irec <= 128; $irec++) {
+
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ $dir_rec = "${dir2}/${sta}.${net}";
+ print "\n $irec out of $nrec -- ${dir_rec}";
+
+ `rm ${dir_rec}/*pdf`; `rm ${dir_rec}/*ps`;
+ }
+ }
+}
+
+#================================================
+print "\n done with setup_misfit_dir.pl\n\n";
+#================================================
+
+sub cd_directory {
+ my($new) = @_;
+ my $old = `pwd`;
+ chomp($old);
+ check_directory($new);
+ #print "$prog: cd $new\n";
+ chdir $new;
+ return($old);
+}
+
+sub check_directory {
+ if(! -e $_[0] ) {
+ print "Directory not found: $_[0]\n";
+ exit(-1);
+ }
+}
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/OLD_version/setup_misfit_dir.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/README 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,7 @@
+Carl Tape, 20-Nov-2008
+
+/svn/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/README
+
+These scripts are in development. The objective is to have a way to combine synthetics and measurements from many different models into a series of plots.
+
+---------
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,604 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 17-Nov-2008
+# plot_seis.pl
+#
+# Text here.
+#
+#
+# EXAMPLE:
+#
+# plot_seis.pl -m CMTSOLUTION_3321590 -n PHL/CI -l 10/190 -a 3321590_T006_T030_m13_STATIONS_ADJOINT -d DATA -s SYN -w MEASUREMENT_WINDOWS_3321590_T006_T030_m13/MEASUREMENT_WINDOWS_3321590_T006_T030_m00 -x 3321590_T006_T030_m13_window_chi/3321590_T006_T030_m00_window_chi -i 3321590/m13/m00/6/30
+#
+# plot_seis.pl -m CMTSOLUTION_3321590 -n PHL/CI -l 10/190 -a 3321590_T006_T030_m13_STATIONS_ADJOINT -d DATA -s SYN -i 3321590/m13/m00/6/30 MEASUREMENT_WINDOWS_3321590_T006_T030_m13 MEASUREMENT_WINDOWS_3321590_T006_T030_m00 3321590_T006_T030_m13_window_chi 3321590_T006_T030_m00_window_chi
+#
+#-----------------------------------
+
+use Getopt::Std;
+use POSIX;
+
+sub Usage{
+ print STDERR <<END;
+
+ plot_seis.pl -m CMTFILE -n station/network -r -a STATION_FILE -d data_dir -s syn_dir -w winfile
+ with
+
+ -m -- CMTSOLUTION file to plot event focal mechanism (beach ball) on the map
+ -l -- Start/End of trace to be cut
+ -n -- station/network
+ -r -- plot 3-compoment seismograms and adjoint sources in relatively scaled amplitude
+ -a -- station file to plot all the stations on the map
+ -d -- data directory
+ -s -- synthetics directory
+ -w -- window file for model 1 / window file for model 2
+ -x -- measurement file for model 1 / measurement file for model 2
+ -i -- event ID / model label 1 / model label 2 / Tmin / Tmax
+END
+ exit(1);
+}
+
+if (@ARGV == 0) { Usage(); }
+if (!getopts('m:l:n:ra:d:s:i:')) {die('Check input arguments\n');}
+print "STATIONS file is $opt_a\n";
+if ($opt_m) {@cmtvals = &get_cmt($opt_m);}
+if ($opt_l) {($tmin,$tmax) = split(/\//,$opt_l);$trange = $tmax-$tmin;} # search for opt_l
+#if ($opt_w) {($win_file1,$win_file2) = split(/\//,$opt_w);}
+#if ($opt_x) {($meas_file1,$meas_file2) = split(/\//,$opt_x);}
+if ($opt_n) {($sta,$net)=split(/\//,$opt_n);}
+if (!$opt_a) {$opt_a = "STATIONS_ADJOINT"}
+if (!$opt_d) {$opt_d = "DATA"}
+if (!$opt_s) {$opt_s = "SYN"}
+if ($opt_i) {($eid,$smodel1,$smodel2,$Tmin,$Tmax) = split(/\//,$opt_i);}
+
+# final four files are: window file 1, window file 2, measurement file 1, measurement file 2
+if(@ARGV != 4) {die("EXIT: must have two window files and two measurement files");}
+ at fourfiles = @ARGV;
+$win_file1 = $fourfiles[0];
+$win_file2 = $fourfiles[1];
+$meas_file1 = $fourfiles[2];
+$meas_file2 = $fourfiles[3];
+
+$sTmin = sprintf("T%3.3i",$Tmin);
+$sTmax = sprintf("T%3.3i",$Tmax);
+$Ttag = "${sTmin}_${sTmax}";
+$Ttag_title = sprintf("bandpass %i s to %i s",$Tmin,$Tmax);
+$Ttag_plot = sprintf("bp [%i s, %i s]",$Tmin,$Tmax);
+
+# name of file containing measurements (see mt_sub.f90)
+#$meas_file1 = "${eid}_${Ttag}_${smodel1}_window_chi";
+#$meas_file2 = "${eid}_${Ttag}_${smodel2}_window_chi";
+
+if (not -f $opt_a) {die("\n check if STATIONS file $opt_a exists\n")}
+if (not -f $meas_file1) {die("\n check if measurment file $meas_file1 exists\n")}
+if (not -f $meas_file2) {die("\n check if measurment file $meas_file2 exists\n")}
+if (not -f $win_file1) {die("\n check if window file $win_file1 exists\n")}
+if (not -f $win_file2) {die("\n check if window file $win_file2 exists\n")}
+
+print "\nCHECKING USER INPUT:\n";
+print "Event ID: ${eid}\n";
+print "Station $sta network $net\n";
+print "Models $smodel1 and $smodel2\n";
+print "Period range $Ttag\n";
+print "STATIONS file is ${opt_a}\n";
+print "Plotting interval is $trange seconds: $tmin to $tmax\n";
+print "Window file 1 is $win_file1\n";
+print "Window file 2 is $win_file2\n";
+print "Measurement file 1 is $meas_file1\n";
+print "Measurement file 2 is $meas_file2\n";
+print "CMTSOLUTION: @cmtvals\n";
+#die("CHECKING USER INPUT");
+
+#---------------------------
+# ADDITIONAL USER INPUT
+
+# plot faults (southern California)
+$ifault = 1;
+
+# width of seismograms depends on whether you plot adjoint sources
+$iplot_win = 1; # plot windows
+$iplot_CClabel = 1; # plot CC measurement labels for each window
+if($iplot_CClabel==1) {$iplot_win = 1;}
+$sdX = 4.5; $sdY = 1.4;
+
+#---------------------------
+
+# extract values from the array cmtvals
+ at M[0..5] = @cmtvals[12..17];
+$year = $cmtvals[0]; $month = $cmtvals[1]; $day = $cmtvals[2];
+$elat = $cmtvals[9]; $elon = $cmtvals[10]; $edep = $cmtvals[11];
+$eid = $cmtvals[18];
+$cmtpsmeca = "temp.txt";
+open(MIJ,">$cmtpsmeca");
+#printf MIJ "$elon $elat $edep $M[0] $M[1] $M[2] $M[3] $M[4] $M[5] 1 X Y $eid";
+printf MIJ "$elon $elat $edep $M[0] $M[1] $M[2] $M[3] $M[4] $M[5] 1";
+close(MIJ);
+
+# compute moment using Dahlen and Tromp, Eq. (5.91)
+# compute moment magnitude using Kanamori (1977) -- see Carl's Latex notes cmt.pdf
+$M0 = 1/sqrt(2) * sqrt($M[0]**2 + $M[1]**2 + $M[2]**2 + 2*($M[3]**2 + $M[4]**2 + $M[5]**2 ) );
+$k = -11.2 + (2/3)*log(5)/log(10);
+$Mw = (2/3)*log($M0)/log(10) + $k;
+$stM0 = sprintf("M0 = %.3e dyne-cm",$M0);
+$stMw = sprintf("Mw = %.2f",$Mw);
+print "\n moment : $stM0 \n moment magnitude : $stMw\n";
+
+$fontno = 1; # font number to use for text
+
+# pstext info for measurement labels
+$textinfo1 = "-C2p -W0/255/255o,0/255/255 -G0/0/0";
+$textinfo2 = "-C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+$textinfo3 = "-C4p -W255/255/255o,1.0p,255/255/255,solid -G0/0/0";
+$textinfo4 = "-C2p -W255/255/255o,1.0p,255/255/255,solid -G0/0/0";
+
+#---------------------------
+
+# set GMT defaults
+`gmtset BASEMAP_TYPE plain PAPER_MEDIA letter ANNOT_FONT_SIZE = 9p LABEL_FONT_SIZE = 9p PAGE_ORIENTATION = landscape PLOT_DEGREE_FORMAT D`;
+
+#$name = "${sta}_${net}_win_adj_${klab}";
+$name = "${eid}_${Ttag}_${sta}_${net}_${smodel1}_${smodel2}_seis";
+$ps_file = "$name.ps";
+$jpg_file = "$name.jpg";
+
+$saclst="/opt/seismo-util/bin/saclst";
+
+#------------
+$xmin = -122; $xmax = -114;
+$ymin = 32; $ymax = 37;
+$bounds="-R$xmin/$xmax/$ymin/$ymax ";
+#$bounds="-R-120/-115/32/36.5";
+$xtick = 1; $ytick = 1;
+$tick="-Ba2f1/a2f1::wEsN";
+$cinfo="-W0.5 -Dh -A100";
+$proj="-JM2.5";
+
+# plot the base map
+$origin = "-X2.25 -Y5.5";
+`psbasemap $bounds $proj $tick $origin -K > $ps_file`; # START
+`pscoast $bounds $proj $cinfo -S150/200/255 -W2 -G200/255/150 -N1/1p -N2/0.5p -O -K >> $ps_file`;
+
+# plot southern California faults
+if($ifault==1) {
+ $dir0 = "/home/carltape/gmt";
+ $fault_file = "$dir0/faults/jennings.xy";
+ $kcf_file = "$dir0/faults/kcf.xy";
+ $breck_file = "$dir0/faults/breck.xy";
+ $fault_infoK = "-M -W0.5p,0/0/0";
+ `psxy $fault_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
+ `psxy $kcf_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
+ `psxy $breck_file $bounds $proj $fault_infoK -O -K >> $ps_file`;
+}
+
+# plot the stations
+`awk \'{print \$4,\$3,0,4,B,C}\' $opt_a | psxy $proj $bounds -M -St10p -W0.5p -G200/200/200 -N -O -K >>$ps_file`;
+
+#------------
+# data, synthetics
+
+# DATA
+$d_tag = "$eid.$net.$sta";
+$Tdat = `ls ${opt_d}/${d_tag}.*T.sac*`; chomp($Tdat);
+$Rdat = `ls ${opt_d}/${d_tag}.*R.sac*`; chomp($Rdat);
+$Zdat = `ls ${opt_d}/${d_tag}.*Z.sac*`; chomp($Zdat);
+if(-f $Tdat) {$Tflag = 1} else {$Tflag = 0}
+if(-f $Rdat) {$Rflag = 1} else {$Rflag = 0}
+if(-f $Zdat) {$Zflag = 1} else {$Zflag = 0}
+#print "\n Tflag : $Tflag, Rflag : $Rflag, Zflag : $Zflag \n";
+
+$smodel = $smodel1;
+
+# SYNTHETICS -- component label always has the form BH_,
+# even if the data file is something else (HH_, LH_).
+# We assume that ALL componenets exist.
+$s_tag = "semd.sac.${smodel}.${Ttag}";
+$Tlab = "$sta.$net.BHT";
+$Rlab = "$sta.$net.BHR";
+$Zlab = "$sta.$net.BHZ";
+$Tsyn = `ls ${opt_s}/${Tlab}*${s_tag}`; chomp($Tsyn);
+$Rsyn = `ls ${opt_s}/${Rlab}*${s_tag}`; chomp($Rsyn);
+$Zsyn = `ls ${opt_s}/${Zlab}*${s_tag}`; chomp($Zsyn);
+
+#------------
+
+#print "\n checking files: \n$Tdat\n$Rdat\n$Zdat\n$Tsyn\n$Rsyn\n$Zsyn\n"; die("TESTING FILES");
+
+# get the receiver location, azimuth, distance (from the $Zsyn file)
+
+(undef,$tmin0,$tmax0,$slon,$slat,$dist,$az,$Tchan) = split(" ",`$saclst b e stlo stla dist az kcmpnm f $Zsyn`);
+$edate = sprintf("%4.4i . %2.2i . %2.2i",$year,$month,$day);
+$stedep = sprintf("depth = %.2f km",$edep);
+$strec = "$sta . $net";
+$stdist = sprintf("dist = %.1f km",$dist);
+$staz = sprintf("az = %.1f deg",$az);
+$stmodel = "model $smodel";
+
+# if the time cut is not specified, then use the whole record
+if (!$opt_l) {$tmin = $tmin0; $tmax = $tmax0;}
+$trange = $tmax - $tmin;
+
+# plot labels next to the map
+ at mapstrings = ($edate,$eid,$stMw,$stedep,$strec,$stdist,$staz,"--",${Ttag_plot},$stmodel);
+$nstrings = @mapstrings;
+$x1 = $xmin - 4*$xtick;
+for ($i=0; $i < $nstrings; $i++) {
+ $y1 = $ymax - $i*$ytick/2;
+ `pstext $bounds $proj $tick -N -K -O >> $ps_file <<EOF\n $x1 $y1 9 0 $fontno LM $mapstrings[$i] \nEOF\n`;
+}
+
+#`pstext $bounds $proj $tick -N -K -O >> $ps_file <<EOF\n $x1 $y0 12 0 $fontno LM $eid \nEOF\n`;
+#`pstext $bounds $proj $tick -N -K -O >> $ps_file <<EOF\n $x1 $y1 12 0 $fontno LM $stMw \nEOF\n`;
+#`pstext $bounds $proj $tick -N -K -O >> $ps_file <<EOF\n $x1 $y2 12 0 $fontno LM $stedep \nEOF\n`;
+#`pstext $bounds $proj $tick -N -K -O >> $ps_file <<EOF\n $x1 $y3 12 0 $fontno LM $strec \nEOF\n`;
+#`pstext $bounds $proj $tick -N -K -O >> $ps_file <<EOF\n $x1 $y4 12 0 $fontno LM $stdist \nEOF\n`;
+#`pstext $bounds $proj $tick -N -K -O >> $ps_file <<EOF\n $x1 $y5 12 0 $fontno LM $staz \nEOF\n`;
+
+# plot ray path, station, and CMT source
+`psxy $bounds $proj -W1.5p -K -O >> $ps_file <<EOF\n$elon $elat\n$slon $slat\nEOF\n`; # ray path
+$cmtfont = 10; $cmt_scale = 0.35;
+$cmtinfo = "-Sm${cmt_scale}/$cmtfont -L0.5p/0/0/0 -G255/0/0";
+`psmeca $cmtpsmeca $bounds $proj $cmtinfo -O -K >> $ps_file`; # CMT source
+`rm $cmtpsmeca`;
+`psxy $bounds $proj -St0.15 -W2 -G255/0/0 -K -O>> $ps_file <<EOF\n$slon $slat\nEOF\n`; # red station
+$slat_shift=$slat-0.35;
+`pstext $bounds $proj $textinfo4 -K -O>> $ps_file <<EOF\n $slon $slat_shift 8 0 $fontno CM $sta \nEOF\n`; # station label
+
+#-----------------------------------
+# limits for seismogram axes
+$Smin = 0;
+$Smax = 2;
+$Srange = $Smax - $Smin;
+
+$xtextpos1 = $tmin - 0.00*$trange; # xmax label
+$xtextpos2 = $tmin + 0.02*$trange; # first window measurement label
+$xtextpos3 = $tmin + 0.50*$trange; # title position
+$xtextpos4 = $tmin + 0.98*$trange; # final window measurement label
+
+$ytextpos5 = 1.35*$Srange; # titles
+#$ytextpos4 = 1.03*$Srange; # dT label
+#$ytextpos3 = 0.90*$Srange; # dlnA label
+$ytextpos4 = 0.90*$Srange; # dT label
+$ytextpos3 = 0.80*$Srange; # dlnA label
+$ytextpos2 = 0.50*$Srange; # Z, R, T label
+$ytextpos1 = 0.15*$Srange; # ymax label
+
+# dimension of seismograms
+$proj = "-JX${sdX}/${sdY}";
+$dYseis = $sdY + 0.00;
+
+$Bseis0 = "-Ba50f10:\" \":/a1f1";
+#$tick = "-Ba50f10:\"Time (s)\":/a1f1S";
+#$tick2 = "-Ba50f10/a1f1S";
+
+$pen = "0.5p"; # line thickness for seismograms
+$red_pen = "-W${pen},250/0/0";
+$black_pen = "-W${pen},0/0/0";
+$recon_pen = "-W${pen},0/255/255";
+
+#============================================================================
+
+for ($kk = 1; $kk <= 2; $kk++) {
+
+ my(@Twinb); my(@Twine);
+ my(@Rwinb); my(@Rwine);
+ my(@Zwinb); my(@Zwine);
+ my($undef);
+
+ if ($kk == 1) {
+ $smodel = $smodel1;
+ $measfile = $meas_file1;
+ $winfile = $win_file1;
+ #$shift = "-X-2 -Y4.5";
+ $shift = "-X-1.5 -Y-4.5";
+
+ } elsif($kk == 2) {
+ $smodel = $smodel2;
+ $measfile = $meas_file2;
+ $winfile = $win_file2;
+ $dX = $sdX + 0.5;
+ $dY = -2*$dYseis;
+ $shift = "-X$dX -Y$dY";
+ }
+
+
+# SYNTHETICS -- component label always has the form BH_,
+# even if the data file is something else (HH_, LH_).
+# We assume that ALL componenets exist.
+$s_tag = "semd.sac.${smodel}.${Ttag}";
+$Tlab = "$sta.$net.BHT";
+$Rlab = "$sta.$net.BHR";
+$Zlab = "$sta.$net.BHZ";
+$Tsyn = `ls ${opt_s}/${Tlab}*${s_tag}`; chomp($Tsyn);
+$Rsyn = `ls ${opt_s}/${Rlab}*${s_tag}`; chomp($Rsyn);
+$Zsyn = `ls ${opt_s}/${Zlab}*${s_tag}`; chomp($Zsyn);
+
+# windows and measurements: transverse component
+$ntline=`grep -n \"$Tsyn\" $winfile | awk -F: '{ print \$1 }'`;
+chomp($ntline);
+if ($ntline) {
+ $nline=$ntline;
+ $nline=$nline+1;
+ $nT=`sed -n ${nline}p $winfile`;
+ @mlines = `grep $Tlab $measfile`; # measurements
+ for ($i=0; $i < $nT; $i++) {
+ $nline=$nline+1;
+ $win=`sed -n ${nline}p $winfile`; chomp($win);
+ (undef,$winb,$wine)=split(/ +/,$win);
+ push @Twinb, "$winb";
+ push @Twine, "$wine";
+
+ (undef,undef,undef,undef,undef,$kplotT,undef,undef,$chi2T,$chi3T,$chi4T,$chi5T,$meas2T,$meas3T,$meas4T,$meas5T,$sigma2T,$sigma3T,$sigma4T,$sigma5T) = split(" ",$mlines[$i]);
+ $stlabT_dT[$i] = sprintf("dT = %.2f \\261 %.2f",$meas4T,$sigma4T);
+ $stlabT_dA[$i] = sprintf("dA = %.2f \\261 %.2f",$meas5T,$sigma5T);
+ }
+}
+#print "Transverse component: ntline = $ntline, nT = $nT\n";
+
+# windows and measurements: radial component
+$nrline=`grep -n \"$Rsyn\" $winfile | awk -F: '{ print \$1 }'`;
+chomp($nrline);
+#print "$nrline\n";
+if ($nrline) {
+ $nline=$nrline;
+ $nline=$nline+1;
+ $nR=`sed -n ${nline}p $winfile`;
+ @mlines = `grep $Rlab $measfile`; # measurements
+ for ($i=0; $i < $nR; $i++) {
+ $nline=$nline+1;
+ $win=`sed -n ${nline}p $winfile`;
+ chomp($win);
+ # print "$win\n";
+ (undef,$winb,$wine)=split(/ +/,$win);
+ # print "$winb\n";
+ # print "$wine\n";
+ push @Rwinb, "$winb";
+ push @Rwine, "$wine";
+
+ (undef,undef,undef,undef,undef,$kplotR,undef,undef,$chi2R,$chi3R,$chi4R,$chi5R,$meas2R,$meas3R,$meas4R,$meas5R,$sigma2R,$sigma3R,$sigma4R,$sigma5R) = split(" ",$mlines[$i]);
+ $stlabR_dT[$i] = sprintf("dT = %.2f \\261 %.2f",$meas4R,$sigma4R);
+ $stlabR_dA[$i] = sprintf("dA = %.2f \\261 %.2f",$meas5R,$sigma5R);
+ }
+}
+#print "Radial component: nrline = $nrline, nR = $nR\n";
+
+# windows and measurements: vertical component
+$nzline=`grep -n \"$Zsyn\" $winfile | awk -F: '{ print \$1 }'`;
+chomp($nzline);
+#print "$nzline\n";
+if ($nzline) {
+ $nline=$nzline;
+ $nline=$nline+1;
+ $nZ=`sed -n ${nline}p $winfile`;
+ @mlines = `grep $Zlab $measfile`; # measurements
+ for ($i=0; $i < $nZ; $i++) {
+ $nline=$nline+1;
+ $win=`sed -n ${nline}p $winfile`;
+ chomp($win);
+ # print "$win\n";
+ (undef,$winb,$wine)=split(/ +/,$win);
+ # print "$winb\n";
+ # print "$wine\n";
+ push @Zwinb, "$winb";
+ push @Zwine, "$wine";
+
+ (undef,undef,undef,undef,undef,$kplotZ,undef,undef,$chi2Z,$chi3Z,$chi4Z,$chi5Z,$meas2Z,$meas3Z,$meas4Z,$meas5Z,$sigma2Z,$sigma3Z,$sigma4Z,$sigma5Z) = split(" ",$mlines[$i]);
+ $stlabZ_dT[$i] = sprintf("dT = %.2f \\261 %.2f",$meas4Z,$sigma4Z);
+ $stlabZ_dA[$i] = sprintf("dA = %.2f \\261 %.2f",$meas5Z,$sigma5Z);
+ }
+}
+#print "Vertical component: nrline = $nrline, nZ = $nZ\n";
+
+#==========================================
+
+# pssac2 scaling
+# -M vertical scaling in sacfile_unit/MEASURE_UNIT = size<required>
+# size: each trace will normalized to size (in MEASURE_UNIT)
+# scale = (1/size) * [data(max) - data(min)]
+
+# limits for synthetic records
+(undef,$minT2,$maxT2) = split(" ",`$saclst depmin depmax f $Tsyn`);
+$maxT2 = max(abs($minT2),abs($maxT2));
+(undef,$minR2,$maxR2) = split(" ",`$saclst depmin depmax f $Rsyn`);
+$maxR2 = max(abs($minR2),abs($maxR2));
+(undef,$minZ2,$maxZ2) = split(" ",`$saclst depmin depmax f $Zsyn`);
+$maxZ2 = max(abs($minZ2),abs($maxZ2));
+
+# use the same vertical scale for m00 and m13 Z, the same vertical scale for m00 and m13 R, etc
+if ($kk == 1) {
+
+ # limits for data records
+ $maxT = 0; $maxR = 0; $maxZ = 0;
+ if ($Tflag==1) {(undef,$minT,$maxT)=split(" ",`$saclst depmin depmax f $Tdat`); $maxT=max(abs($minT),abs($maxT));}
+ if ($Rflag==1) {(undef,$minR,$maxR)=split(" ",`$saclst depmin depmax f $Rdat`); $maxR=max(abs($minR),abs($maxR));}
+ if ($Zflag==1) {(undef,$minZ,$maxZ)=split(" ",`$saclst depmin depmax f $Zdat`); $maxZ=max(abs($minZ),abs($maxZ));}
+
+ # overall max values
+ $maxTDS = max($maxT,$maxT2);
+ $maxRDS = max($maxR,$maxR2);
+ $maxZDS = max($maxZ,$maxZ2);
+ $maxD = max($maxT,max($maxR,$maxZ));
+ $maxS = max($maxT2,max($maxR2,$maxZ2));
+ $max = max($maxD,$maxS);
+ #print "\n Overall max values: $maxD (data), $maxS (syn) $max (overall)\n";
+
+ $bounds = "-R$tmin/$tmax/$Smin/$Smax";
+ $scale = 1.2; # KEY: scaling for plotting seismograms
+ $size = 1/$scale;
+
+}
+
+if ($opt_r) {
+ # normalize by the same value (max)
+ $sizeT = $size*$maxT/$max;
+ $sizeT2 = $size*$maxT2/$max;
+ $sizeR = $size*$maxR/$max;
+ $sizeR2 = $size*$maxR2/$max;
+ $sizeZ = $size*$maxZ/$max;
+ $sizeZ2 = $size*$maxZ2/$max;
+
+} else {
+ # normalize by different values
+ $sizeT = $size*$maxT/$maxTDS;
+ $sizeT2 = $size*$maxT2/$maxTDS;
+ $sizeR = $size*$maxR/$maxRDS;
+ $sizeR2 = $size*$maxR2/$maxRDS;
+ $sizeZ = $size*$maxZ/$maxZDS;
+ $sizeZ2 = $size*$maxZ2/$maxZDS;
+}
+
+# strings for plotting -- show the yaxis limits for each trace (based on synthetics)
+$strT = sprintf('%.2e',$maxT2/$sizeT2);
+$strR = sprintf('%.2e',$maxR2/$sizeR2);
+$strZ = sprintf('%.2e',$maxZ2/$sizeZ2);
+
+#==========================================
+
+# TRANSVERSE component: data, synthetics, and windows
+$Bseis = "${Bseis0}::weSn";
+`pssac2 $shift $proj $Tsyn $bounds -Ent-3 -M${sizeT2} $red_pen -N -K -O $Bseis >> $ps_file`; # synthetics
+if ($ntline) {
+ 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);
+ }
+ }
+ 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 $textinfo1 $proj $bounds -N -O -K >> $ps_file`;
+ `echo \"$xtext $ytextpos3 7 0 1 $just $stlabT_dA[$i]\" | pstext $textinfo1 $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 && ($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 $textinfo2 -N -O -K >> $ps_file`;
+`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strT\" | pstext $proj $bounds $textinfo3 -N -O -K >> $ps_file`;
+
+# RADIAL component: data, synthetics, and windows
+$Bseis = "${Bseis0}::wesn";
+`pssac2 -Y$dYseis $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_pen -N -K -O $Bseis >> $ps_file`;
+if ($nrline) {
+ 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);
+ }
+ }
+ 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 $textinfo1 $proj $bounds -N -O -K >> $ps_file`;
+ `echo \"$xtext $ytextpos3 7 0 1 $just $stlabR_dA[$i]\" | pstext $textinfo1 $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`;}
+`pssac2 $proj $Rsyn $bounds -Ent-3 -M${sizeR2} $red_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 $textinfo2 -N -O -K >> $ps_file`;
+`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strR\" | pstext $proj $bounds $textinfo3 -N -O -K >> $ps_file`;
+
+# VERTICAL component: data, synthetics, and windows
+$Bseis = "${Bseis0}::wesn";
+`pssac2 -Y$dYseis $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -K -O $Bseis >> $ps_file`;
+if ($nzline) {
+ 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);
+ }
+ }
+ 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 $textinfo1 $proj $bounds -N -O -K >> $ps_file`;
+ `echo \"$xtext $ytextpos3 7 0 1 $just $stlabZ_dA[$i]\" | pstext $textinfo1 $proj $bounds -N -O -K >> $ps_file`;
+ }
+ }
+}
+if ($Zflag==1) {`pssac2 $proj $Zdat $bounds -Ent-3 -M${sizeZ} $black_pen -N -K -O $Bseis >> $ps_file`;}
+`pssac2 $proj $Zsyn $bounds -Ent-3 -M${sizeZ2} $red_pen -N -O -K >> $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 $textinfo2 -N -K -O >> $ps_file`;
+`echo \"$xtextpos1 $ytextpos1 9 0 1 CM $strZ\" | pstext $proj $bounds $textinfo3 -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`;
+
+
+} # for loop
+
+#----------------------------------------
+
+# 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 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 ; xv $jpg_file &`;
+`ps2pdf $ps_file`;
+#`rm $ps_file`; # remove PS file
+system("ghostview $ps_file &");
+
+print "done with ${name}.pdf\n";
+
+#================================================
+
+# subroutine name: max
+# Input: number1, number2
+# returns greater of 2 numbers
+sub max {
+ if ($_[0]<$_[1]) {return $_[1]} else {return $_[0]};
+}
+
+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,$eid)=split(" ",$cmt[1]);
+ 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,$eid);
+}
+
+#================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis_all.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis_all.pl 2008-11-24 00:05:12 UTC (rev 13377)
@@ -0,0 +1,137 @@
+#!/usr/bin/perl -w
+
+#-----------------------------------
+# Carl Tape, 03-Sept-2008
+# plot_seis_all.pl
+#
+# This script xxx
+#
+# EXAMPLE: plot_seis_all.pl m13 m00 6/30
+#
+#-----------------------------------
+
+if (@ARGV < 3) {die("Usage: plot_seis_all.pl smodel1 smodel2 Tmin/Tmax\n")}
+($smodel1,$smodel2,$Ts) = @ARGV;
+
+$pwd = $ENV{PWD};
+print "\n$pwd\n";
+
+# labels for bandpass-filtered data and synthetics
+($Tmin,$Tmax) = split("/",$Ts);
+$sTmin = sprintf("T%3.3i",$Tmin);
+$sTmax = sprintf("T%3.3i",$Tmax);
+$Ttag = "${sTmin}_${sTmax}";
+
+#-------------------------------------------
+# USER INPUT
+
+# directory containing all CMTSOLUTION files
+$dir_src = "/net/sierra/raid1/carltape/results/SOURCES/socal_11";
+$dir_cmt = "${dir_src}/CMT_files_pre_inverted";
+if (not -e $dir_src) {die("check if dir_src $dir_src exist or not\n")}
+if (not -e $dir_cmt) {die("check if dir_cmt $dir_cmt exist or not\n")}
+
+# list of event IDs to use
+$file_eid0 = "/net/sierra/raid1/carltape/results/EID_LISTS/syn_run_m12";
+if (not -f ${file_eid0}) {die("check if file_eid ${file_eid0} exist or not\n")}
+
+# FULL stations file
+$file_stations0 = "/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem";
+if (not -f ${file_stations0}) {die("check if file_stations ${file_stations0} exist or not\n")}
+
+# data and synthetics directories
+$dir_data0 = "/net/sierra/raid1/carltape/socal/socal_3D/DATA/FINAL";
+$dir_syn01 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${smodel1}";
+$dir_syn02 = "/net/sierra/raid1/carltape/socal/socal_3D/SYN/model_${smodel2}";
+if (not -e ${dir_data0}) {die("check if dir_data ${dir_data0} exist or not\n")}
+if (not -e ${dir_syn01}) {die("check if dir_syn ${dir_syn01} exist or not\n")}
+if (not -e ${dir_syn02}) {die("check if dir_syn ${dir_syn02} exist or not\n")}
+
+# directory containing the output for windows, measurements, adjoint sources, etc
+$dir_run = "/net/sierra/raid1/carltape/socal/socal_3D/RUNS";
+if (not -e ${dir_run}) {die("check if dir_run ${dir_run} exist or not\n")}
+
+#-------------------------------------------
+
+# open EID list
+open(IN,"${file_eid0}"); @eids = <IN>; close(IN);
+$nevent = @eids;
+print "\n $nevent events in the list";
+
+# open STATIONS file
+open(IN,"${file_stations0}"); @slines = <IN>; close(IN);
+$nrec = @slines - 1;
+print "\n $nrec stations in the list";
+
+#-------------------------------------------
+
+if (0==1) {
+ for ($ievent = 1; $ievent <= $nevent; $ievent++) {
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+ }
+ die("testing");
+}
+
+#=============================================
+# GENERATE GMT FIGURES
+
+#$imin = 1; $imax = $nevent;
+#$imin = 1; $imax = 100;
+$imin = 7; $imax = $imin;
+
+for ($ievent = $imin; $ievent <= $imax; $ievent++) {
+
+ $eid = $eids[$ievent-1]; chomp($eid);
+ print "\n $ievent : $eid";
+
+ $cmtfile = "${dir_cmt}/CMTSOLUTION_${eid}";
+ $stafile = "${dir_run}/${eid}/${smodel1}/MEASURE_${Ttag}/ADJOINT_${Ttag}/STATIONS_ADJOINT";
+ $winfile1 = "${dir_run}/${eid}/${smodel1}/WINDOW_${Ttag}/MEASUREMENT_WINDOWS_${eid}_${Ttag}_${smodel1}";
+ $winfile2 = "${dir_run}/${eid}/${smodel2}/WINDOW_${Ttag}/MEASUREMENT_WINDOWS_${eid}_${Ttag}_${smodel2}";
+ $measfile1 = "${dir_run}/${eid}/${smodel1}/MEASURE_${Ttag}/${eid}_${Ttag}_${smodel1}_window_chi";
+ $measfile2 = "${dir_run}/${eid}/${smodel2}/MEASURE_${Ttag}/${eid}_${Ttag}_${smodel2}_window_chi";
+
+ if (not -f $cmtfile) {die("check if cmtfile $cmtfile exist or not\n")}
+ if (not -f $stafile) {die("check if stafile $stafile exist or not\n")}
+ if (not -f $winfile1) {die("check if winfile1 $winfile1 exist or not\n")}
+ if (not -f $winfile2) {die("check if winfile2 $winfile2 exist or not\n")}
+ if (not -f $measfile1) {die("check if measfile1 $measfile1 exist or not\n")}
+ if (not -f $measfile2) {die("check if measfile2 $measfile2 exist or not\n")}
+
+ # link the data and synthetics into the temporary directory
+ $datadir = "${dir_data0}/${eid}/PROCESSED/PROCESSED_${Ttag}";
+ $syndir1 = "${dir_syn01}/${eid}/PROCESSED/PROCESSED_${Ttag}";
+ $syndir2 = "${dir_syn02}/${eid}/PROCESSED/PROCESSED_${Ttag}";
+ if (not -e $datadir) {die("check if datadir $datadir exist or not\n")}
+ if (not -e $syndir1) {die("check if syndir1 $syndir1 exist or not\n")}
+ if (not -e $syndir2) {die("check if syndir2 $syndir2 exist or not\n")}
+ `rm DATA/* SYN/*`; # remove previous links
+ `ln -s $datadir/*BHZ* DATA`; `ln -s $datadir/*BHR* DATA`; `ln -s $datadir/*BHT* DATA`;
+ `ln -s $datadir/*HHZ* DATA`; `ln -s $datadir/*HHR* DATA`; `ln -s $datadir/*HHT* DATA`;
+ `ln -s $syndir1/*BHZ* SYN`; `ln -s $syndir1/*BHR* SYN`; `ln -s $syndir1/*BHT* SYN`;
+ `ln -s $syndir2/*BHZ* SYN`; `ln -s $syndir2/*BHR* SYN`; `ln -s $syndir2/*BHT* SYN`;
+
+ #$rmin = 1; $rmax = $nrec;
+ $rmin = 272; $rmax = $rmin;
+
+ for ($irec = $rmin; $irec <= $rmax; $irec++) {
+
+ ($sta,$net,$slat,$slon,undef,undef) = split(" ",$slines[$irec]);
+ $station = "$sta.$net";
+ print "\n $irec out of $nrec -- $station";
+
+ #$lcut = "-20/100";
+ $lcut = "10/180";
+
+ print "\n plot_seis.pl -m $cmtfile -n $sta/$net -l $lcut -a $stafile -d DATA -s SYN -i ${eid}/${smodel1}/${smodel2}/${Ts} $winfile1 $winfile2 $measfile1 $measfile2\n";
+ `plot_seis.pl -m $cmtfile -n $sta/$net -l $lcut -a $stafile -d DATA -s SYN -i ${eid}/${smodel1}/${smodel2}/${Ts} $winfile1 $winfile2 $measfile1 $measfile2`;
+
+ }
+}
+
+#================================================
+print "\n done with plot_seis_all.pl\n\n";
+#================================================
+
+#=================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot/plot_seis_all.pl
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list