[cig-commits] r16252 - in seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt: . scripts
carltape at geodynamics.org
carltape at geodynamics.org
Tue Feb 9 16:31:10 PST 2010
Author: carltape
Date: 2010-02-09 16:31:09 -0800 (Tue, 09 Feb 2010)
New Revision: 16252
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/GJI_EDW2_ker.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_one.pl
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl
Log:
Updated Perl-GMT plotting scripts for tomographic cross sections.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl 2010-02-10 00:30:55 UTC (rev 16251)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -4,7 +4,7 @@
#
# plot_horz_coverage.pl
# Carl Tape
-# 07-Dec-2008
+# 16-Nov-2009
#
# horz_xc_vs_coverage_003.dat
# This script reads in a horizontal cross-section data file with five columns:
@@ -12,16 +12,13 @@
# and plots the three cross-sections.
#
# EXAMPLE (SCIENCE PAPER):
-# plot_horz_coverage.pl 0 3/3 0/4 0/1/6 # no mask
-# plot_horz_coverage.pl 1 3/3 0/4 0/0/6 # with mask
-# plot_horz_coverage.pl 0 11/11 0/4 0/1/3
-# plot_horz_coverage.pl 1 11/11 0/4 0/0/3
+# plot_horz_coverage.pl 0 5/5 0/4 0/1/6 # no mask
+# plot_horz_coverage.pl 1 5/5 0/4 0/0/6 # with mask
# plot_horz_coverage.pl 0 21/21 0/4 0/1/3
# plot_horz_coverage.pl 1 21/21 0/4 0/0/3
+# plot_horz_coverage.pl 0 41/41 0/4 0/1/3
+# plot_horz_coverage.pl 1 41/41 0/4 0/0/3
#
-# plot_horz_coverage.pl 0 1/1 0/4 0/1/6 # no mask
-# plot_horz_coverage.pl 1 1/1 0/4 0/0/6 # with mask
-#
#==========================================================
if (@ARGV < 4) {die("Usage: plot_horz_coverage.pl xxx\n");}
@@ -40,10 +37,10 @@
#$ilonlab = 0; # display longitude tick labels at the top
#$itoplab = 1; # label the model above the plot
#$isidelab = 1; # label of the depth on the left of the plot
-$ivs = 0; # KEY COMMAND
+$ivs = 1; # KEY COMMAND
$iinsetlab = 0; # label the model inan inset
$iletter = 0; # label for publication figures
-$letter = "A";
+$letter = "C";
$cgray = "-G220/220/220";
@@ -53,10 +50,11 @@
# USER INPUT
# directory containing the data files
-$dirdat = "INPUT/horz_01";
+$stirun = "02";
+$dirdat = "INPUT/horz_${stirun}";
# KEY: file showing the cuts and the plotting range for the values
-$fcuts = "${dirdat}/horz_xc_elevation_cuts_mod";
+$fcuts = "${dirdat}/horz_${stirun}_xc_cuts_mod";
if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
# load values and check that the files exist
@@ -68,7 +66,7 @@
$vsnorms[$p-1] = $vs_norm;
$vbnorms[$p-1] = $vb_norm;
$zcuts[$p-1] = -$zcut;
- $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}.dat";
+ $dfiles[$p-1] = "${dirdat}/horz_${stirun}_xc_${modlab}_m16_${smodeltag}_${stip}.dat";
#if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
}
print "\nNumber of cuts is $nump\n";
@@ -233,7 +231,7 @@
for ($p = $pmin; $p <= $pmax; $p ++ ) {
$stip = sprintf("%3.3i",$p);
- $fname = "horz_xc_${modlab}_${smodeltag}_${stip}_${imask}";
+ $fname = "horz_${stirun}_xc_${modlab}_m16_${smodeltag}_${stip}_${imask}";
$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
$vsnorm = $vsnorms[$p-1]; # reference velocity for m00 and m16
@@ -242,7 +240,7 @@
$zdep = $zcuts[$p-1]; # depth of the cross-section (m)
# zero-level for mask normalization
- $cfile = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}_mask_value.dat";
+ $cfile = "${dirdat}/horz_${stirun}_xc_${modlab}_m16_${smodeltag}_${stip}_mask_value.dat";
if (not -f $cfile) {die("Check if cfile $cfile exist or not\n")}
open(IN,$cfile); $cnorm = <IN>; close(IN); chomp($cnorm);
#$cnorms[$p-1] = $cnorm;
@@ -295,7 +293,9 @@
print CSH "grdimage $grdfile -C${cptfile} $J -Q -K -O -V >> $psfile\n";
print CSH "psscale -C${cptfile} $Dscale $Bscale -K -O -V >> $psfile\n";
- $slab1 = sprintf("ln ( %s / %.1e )",$mtit,$cnorm);
+ #$slab1 = sprintf("ln ( %s / %.1e )",$mtit,$cnorm);
+ $slab1 = sprintf("ln ( %s / 10\@+%2.2i\@+ )",$mtit,log($cnorm)/log(10));
+
print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
}
print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile \n";
@@ -381,4 +381,6 @@
# if ($ipdf==1) {system("ps2pdf $psfile")}
# if ($ixv==1) {system("ghostview $psfile &");}
+system("gv $psfile &")
+
#==================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2010-02-10 00:30:55 UTC (rev 16251)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -11,16 +11,16 @@
# and plots the three cross-sections.
#
# EXAMPLE:
-# plot_horz_models.pl 1/50 m00/m16 1 1/1/1/0.2 # Vs
-# plot_horz_models.pl 1/50 m00/m16 0 1/1/1/0.3 # Vb
+# plot_horz_models.pl 1/81 m00/m16 1 1/1/1/0.2 # Vs
+# plot_horz_models.pl 1/81 m00/m16 0 1/1/1/0.3 # Vb
#
# FIGURES:
-# plot_horz_models.pl 3/3 m00/m16 1 1/0/1/0.2 # Vs
-# plot_horz_models.pl 11/11 m00/m16 1 0/0/1/0.2 # Vs
+# plot_horz_models.pl 5/5 m00/m16 1 1/0/1/0.2 # Vs
# plot_horz_models.pl 21/21 m00/m16 1 0/0/1/0.2 # Vs
-# plot_horz_models.pl 3/3 m00/m16 0 1/0/1/0.3 # Vb
-# plot_horz_models.pl 11/11 m00/m16 0 0/0/1/0.3 # Vb
+# plot_horz_models.pl 41/41 m00/m16 1 0/0/1/0.2 # Vs
+# plot_horz_models.pl 5/5 m00/m16 0 1/0/1/0.3 # Vb
# plot_horz_models.pl 21/21 m00/m16 0 0/0/1/0.3 # Vb
+# plot_horz_models.pl 41/41 m00/m16 0 0/0/1/0.3 # Vb
#
# plot_horz_models.pl 3/3 m00/m01 1 1/0/1/0.2 # Vs
#
@@ -39,21 +39,26 @@
# detailed flags for publication figures
#$ilonlab = 0; # display longitude tick labels at the top
#$itoplab = 0; # label the model above the plot
-$iinsetlab = 1; # label the model inan inset
-$iletter = 0; # label for publication figures
-$letter = "c";
-$iextra = 0; # extra information on figures
+$iinsetlab = 1; # label the model inan inset
+$iletter = 0; # label for publication figures
+$letter = "C";
+$iextra = 0; # extra information on figures
+$irun = 1; # default 2; 1 for m01/m00
+$imask = 1;
+$isingle = 0; # single figure showing only final model
$cgray = "-G220/220/220";
+$stirun = sprintf("%2.2i",$irun);
+
#---------------------------------------------------------
# USER INPUT
# directory containing the data files
-$dirdat = "INPUT/horz_01";
+$dirdat = "INPUT/horz_${stirun}";
# KEY: file showing the cuts and the plotting range for the values
-$fcuts = "${dirdat}/horz_xc_elevation_cuts_mod";
+$fcuts = "${dirdat}/horz_${stirun}_xc_cuts_mod";
if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
# load values and check that the files exist
@@ -65,7 +70,7 @@
$vsnorms[$p-1] = $vs_norm;
$vbnorms[$p-1] = $vb_norm;
$zcuts[$p-1] = -$zcut;
- $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}.dat";
+ $dfiles[$p-1] = "${dirdat}/horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}_mask${imask}.dat";
#if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
}
print "\nNumber of cuts is $nump\n";
@@ -90,7 +95,7 @@
# position of titles and labels
$xtx1 = 0.5; $ytx1 = 1.00;
$xtx2 = -0.15; $ytx2 = 0.4;
-$xtx3 = 0.7; $ytx3 = -0.2; # cbar 1
+$xtx3 = 0.7; $ytx3 = -0.19; # cbar 1
$xtx4 = $xtx3; $ytx4 = -0.1; # cbar 2
$xtx5 = -0.15; $ytx5 = 0.85; # A, B, C, etc
$xtx6 = 0.0; $ytx6 = 0.05; # m00, m16, ln(m16/m00)
@@ -209,7 +214,7 @@
$textinfo3 = "-G0/0/0 -S2p,255/255/0";
# BOOLEAN: WHICH FIGURES TO PLOT
-$ixv = 0; $ipdf = 1;
+$ixv = 1; $ipdf = 1;
#$ixv = 1; $ipdf = 0;
# plot title
@@ -229,7 +234,8 @@
for ($p = $pmin; $p <= $pmax; $p ++ ) {
$stip = sprintf("%3.3i",$p);
- $fname = "horz_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}";
+ $fname0 = "horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}";
+ $fname = $fname0;
$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
$vsnorm = $vsnorms[$p-1]; # reference velocity for m00 and m16
@@ -435,14 +441,88 @@
print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n $x_title $y_title 16 0 $fontno CM \nEOF\n"; # FINISH
#if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
+ #if($ixv==1) {print CSH "gv $psfile &\n";}
+ if($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+
+ #-------------------------------------------
+
+ if ($isingle == 1) {
+
+ $fname = "${fname0}_single";
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+ $k = 2;
+
+ $shift = $shifts[$k-1];
+ $B = "$B0"."wENs";
+ $title = $titles[$k-1];
+
+ print CSH "psbasemap $J $R $B -K -V $orient $origin > $psfile\n"; # START
+
+ print CSH "psbasemap $J $R $B $cgray -K -O -V >> $psfile\n";
+ if ($icolor==1) {
+ print CSH "awk '{print \$1,\$2,\$4/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1b $J -Q -K -O -V >> $psfile\n";
+ print CSH "psscale -C${cptfile1b} $Dscale2 $Bscale1b -K -O -V >> $psfile\n";
+ $slab1 = sprintf("(%.2f \261 %.0f \045)",$cnorm/1000,$cpert1*100);
+ $slab2 = "$mtit km/s";
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 $fsize2 0 $fontno LM $slab2\nEOF\n";
+ }
+ print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile \n";
+ if($ishelf==1) {print CSH "awk '{print \$2,\$1}' ${shelf_file} |psxy $J $R $Wshelf -K -O -P -V >> $psfile\n";}
+ print CSH "psxy ${fault_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+
+ print CSH "psbasemap $J $R $B -K -V -O >> $psfile\n";
+
+ # plot cities
+ ##print CSH "psxy $J $R $city_info -K -O -V >>$psfile<<EOF\n -118.1717 34.1483 \nEOF\n"; # Pasadena
+ #print CSH "awk '{print \$2,\$1}' $cities | psxy ${city_info} $J $R -K -V -O >> $psfile \n";
+ #print CSH "awk '{print \$2,\$1,12,0,1,\"LB\",\$3}' $cities | pstext $textinfo3 -D0.05/0.05 $J $R -K -V -O >> $psfile \n";
+
+ #print CSH "pstext ${topo_labels} $J $R $textinfo -N -K -V -O >> $psfile\n";
+ #print CSH "pstext ${fault_labels} $textinfo2 $J $R -K -V -O >> $psfile \n";
+
+ # fault labels
+ $fault_labs = "$dir0/faults/socal_fault_labs_red.lonlat";
+ $textinfo = "-N -C2p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfo $J $R -K -V -O >> $psfile\n";
+
+ #----------------------------
+
+ # plot horizontal title above each column
+ if ($itoplab==1) {
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx1 $ytx1 $fsize0 0 $fontno CM $title\nEOF\n";
+ }
+
+ # inset label for each plot
+ if ($iinsetlab==1) {
+ $lab = $labs[$k-1];
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V >>$psfile<<EOF\n$xtx6 $ytx6 11 0 $fontno LB $lab\nEOF\n";
+ }
+
+ # plot vertical title left of the row
+ if ($isidelab==1) {
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx2 $ytx2 $fsize0 90 $fontno CM $dtitle\nEOF\n";
+ }
+
+ # plot overall label (for publication)
+ if ($iletter==1) {
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V >>$psfile<<EOF\n$xtx5 $ytx5 18 0 $fontno TL $letter\nEOF\n";
+ }
+
+ print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n $x_title $y_title 16 0 $fontno CM \nEOF\n"; # FINISH
+ if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
if($ixv==1) {print CSH "gv $psfile &\n";}
if($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+ }
+
} # loop over p
#------------------------------------
-# print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n $x_title $y_title 16 0 $fontno CM \nEOF\n"; # FINISH
-
close (CSH);
system("csh -f $cshfile");
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl 2010-02-10 00:30:55 UTC (rev 16251)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -12,8 +12,10 @@
#
# EXAMPLE:
# plot_vert_models.pl 1/86 m00/m16 -40/5 1 0.10 3.0 1 1.5 # Vs -- src-rec paths
-# plot_vert_models.pl 1/32 m00/m16 -40/5 1 0.10 3.0 2 1.5 # Vs -- SAF-perp profiles and three more
+# plot_vert_models.pl 4/32 m00/m16 -40/5 1 0.10 3.0 2 1.5 # Vs -- SAF-perp profiles and three more
# plot_vert_models.pl 1/10 m00/m16 -40/5 1 0.10 3.0 3 1.5 # Vs -- Garlock-perp profiles
+# plot_vert_models.pl 1/13 m00/m16 -40/5 1 0.10 3.0 4 1.5 # Vs -- N-S (constant longitude)
+# plot_vert_models.pl 1/9 m00/m16 -40/5 1 0.10 3.0 5 1.5 # Vs -- E-W (constant latitude)
#
# plot_vert_models.pl 1/8 m00/m16 -40/5 1 0.10 3.0 1 1.5 # Vs
# plot_vert_models.pl 1/8 m00/m16 -40/5 2 0.10 3.0 1 1.5 # Vb
@@ -22,9 +24,17 @@
# plot_vert_models.pl 1/3 m00/m16 -40/5 1 0.15 3.0 2 1.5 # Vs
#
#--------------------
+# Science figure: 25,5
# GJI figures: 15,25,2,5,10,9,6,79,19,23,7,49,11
# plot_vert_models.pl 11/11 m00/m16 -40/5 1 0.10 3.0 1 1.5
#
+# plot_vert_models.pl 25/25 m00/m16 -40/5 1 0.10 3.0 1 1.5
+#
+# plot_vert_models.pl 6/6 m00/m16 -40/5 1 0.10 3.0 1 1.5 # GJI - LARSE Vs
+# plot_vert_models.pl 6/6 m00/m16 -40/5 2 0.10 3.0 1 1.5 # GJI - LARSE Vb
+# plot_vert_models.pl 6/6 m00/m16 -40/5 3 0.10 3.0 1 1.5 # GJI - LARSE Vp
+# plot_vert_models.pl 9/9 m00/m16 -40/5 1 0.10 3.0 4 1.5 # GJI - 119 long
+#
#==========================================================
if (@ARGV < 8) {die("Usage: plot_vert_models.pl xxx\n");}
@@ -35,7 +45,7 @@
$isimbox = 1; # box around simulation region
$iletter = 0; # label for publication figures
$letter = "a";
-$imfinal = 1; # plot an additional figure with only m16 and the map
+$imfinal = 0; # plot an additional figure with only m16 and the map
($pmin,$pmax) = split(/\//,$pinds);
($smodel1,$smodel2) = split(/\//,$smodels);
@@ -149,7 +159,7 @@
# BOOLEAN: WHICH FIGURES TO PLOT
#$ixv = 1; $ipdf = 0;
$ixv = 0; $ipdf = 1;
-$ijpg = 0;
+$ijpg = 1;
# plot title
$J_title = "-JM1";
@@ -166,8 +176,7 @@
print "\nNumber of cuts is $nump\n";
# load fault positions along each profile
- #$ffile = "$dirdat/ALL_rays_fault_positions";
- $ffile = "$dirdat/ALL_rays_fault_positions_mod";
+ $ffile = "$dirdat/ALL_${stirun}_rays_fault_positions_mod";
if (not -f $ffile) {die("Check if ffile $ffile exist or not\n");}
open(IN,$ffile); @faults = <IN>; close(IN);
@@ -182,13 +191,13 @@
$stip = sprintf("%3.3i",$p);
# datafile -- COLUMNS: lon, lat, m00, m16
- $dfile = "${dirdat}/vert_xc_vsvbvp_${smodeltag}_${stip}.dat";
+ $dfile = "${dirdat}/vert_${stirun}_xc_vsvbvp_${smodeltag}_${stip}.dat";
if (not -f $dfile) {die("Check if dfile $dfile exist or not\n");}
print "Data file is $dfile\n";
# source and receiver points -- NOT necessarily the endpoints of the ray path
# COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
- $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
open(IN,$rfile); @line = <IN>; close(IN);
($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
@@ -372,7 +381,7 @@
}
# plot ray path
- $pfile = "${dirdat}/vert_xc_${stip}_ray_path";
+ $pfile = "${dirdat}/vert_${stirun}_xc_${stip}_ray_path";
if (not -f $pfile) {
die("Check if pfile $pfile exist or not\n");
}
@@ -380,7 +389,7 @@
print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W3p,0/255/255 -K -O -V $omap >> $psfile\n";
# plot source and receiver lon-lat
- $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
if($irun==1) {
print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmap $Rmap $recinfo_map -K -O -V $omap >> $psfile\n";
@@ -433,6 +442,14 @@
#print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J -Q -T -K -O -V $oxc >> $psfile\n";
+
+ # -A+r0.005 OR -A-
+ # -C0.25 -L2.0/4.0 -Q100
+ $anot = "-A+r0.005+a0 -Gd4.5c:20p";
+ if($ivel==1) {$coninfo = "-C0.5 -L2.0/4.0 -Q100 -W0.75p $anot";}
+ if($ivel==2) {$coninfo = "-C0.50 -L2.0/5.0 -Q100 -W0.75p $anot";}
+ if($ivel==3) {$coninfo = "-C0.50 -L3.0/7.0 -Q100 -W0.75p $anot";}
+ print CSH "grdcontour $grdfile $coninfo $J -K -O -V $oxc >> $psfile\n";
}
print CSH "psscale -C${cptfile1} $Dscale $Bscale1 -K -O -V $ocb >> $psfile\n";
@@ -452,6 +469,8 @@
#print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J -Q -K -O -V $oxc >> $psfile\n";
+
+ print CSH "grdcontour $grdfile $coninfo $J -K -O -V $oxc >> $psfile\n";
}
print CSH "psscale -C${cptfile1} $Dscale $Bscale1 -K -O -V $ocb >> $psfile\n";
@@ -574,9 +593,9 @@
print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n"; # FINISH
#if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
- if ($ixv==1) {print CSH "ghostview $psfile &\n";}
+ if ($ixv==1) {print CSH "gv $psfile &\n";}
if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
- if ($ijpg==1) {print CSH "convert $psfile $jpgfile\n";}
+ if ($ijpg==1) {print CSH "convert -rotate $rotangle $psfile $jpgfile\n";}
#==================================================
# figure with map with m16 beneath
@@ -593,7 +612,7 @@
$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
# datafile -- COLUMNS: lon, lat, m00, m16
- $dfile = "${dirdat}/vert_xc_vsvbvp_${smodeltag}_${stip}.dat";
+ $dfile = "${dirdat}/vert_${stirun}_xc_vsvbvp_${smodeltag}_${stip}.dat";
if (not -f $dfile) {
die("Check if dfile $dfile exist or not\n");
}
@@ -601,7 +620,7 @@
# source and receiver points -- NOT necessarily the endpoints of the ray path
# COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
- $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
if (not -f $rfile) {
die("Check if rfile $rfile exist or not\n");
}
@@ -664,6 +683,8 @@
#---------------------
+ $imap = 0;
+ if($imap==1) {
print CSH "gmtset TICK_LENGTH $tlen_map\n"; # ticks for map
print CSH "psbasemap $Jmap $Rmap $Bmap -K -V $orient $omap > $psfile\n"; # START
@@ -680,7 +701,7 @@
}
# plot ray path
- $pfile = "${dirdat}/vert_xc_${stip}_ray_path";
+ $pfile = "${dirdat}/vert_${stirun}_xc_${stip}_ray_path";
if (not -f $pfile) {die("Check if pfile $pfile exist or not\n");}
print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W6p,0/0/0 -K -O -V $omap >> $psfile\n";
print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W3p,0/255/255 -K -O -V $omap >> $psfile\n";
@@ -721,6 +742,7 @@
}
print CSH "gmtset TICK_LENGTH $tlen\n"; # restore tick
+ } # imap
#-----------------------------------------------
@@ -741,13 +763,19 @@
#------------------
- print CSH "psbasemap $J $R $B -K -O -V $oxc >> $psfile\n";
+ if ($imap==1) {
+ print CSH "psbasemap $J $R $B -K -O -V $oxc > $psfile\n";
+ } else {
+ print CSH "psbasemap $J $R $B -K -V $oxc > $psfile\n";
+ }
print CSH "psbasemap $J $R $B -G${sky_color} -K -O -V $oxc >> $psfile\n";
if ($icolor==1) {
print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
#print CSH "awk '{print \$1,\$2,\$4}' $dfile | xyz2grd -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1 $J -Q -K -O -V $oxc >> $psfile\n";
+
+ print CSH "grdcontour $grdfile -Cconfile2_${modlab} -A- -W0.75p $J -K -O -V $oxc >> $psfile\n";
}
print CSH "psscale -C${cptfile1} $Dscale $Bscale1 -K -O -V $ocb >> $psfile\n";
@@ -777,6 +805,8 @@
#}
#}
+ print CSH "psxy $J $R -W1p,0/0/0,-- -K -O -V $oxc >>$psfile<<EOF\n$dmin -10\n$dmax -10\nEOF\n";
+
# plot source and receiver
if ($irun == 1) {
print CSH "awk '{print \$3,-\$4}' $rfile | psxy $J $R $srcinfo_xc -K -O -V $oxc >> $psfile\n";
@@ -828,9 +858,9 @@
print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n"; # FINISH
#if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
- if ($ixv==1) {print CSH "ghostview $psfile &\n";}
+ if ($ixv==1) {print CSH "gv $psfile &\n";}
if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
- if ($ijpg==1) {print CSH "convert $psfile $jpgfile\n";}
+ if ($ijpg==1) {print CSH "convert -rotate $rotangle $psfile $jpgfile\n";}
} # imfinal == 1
} # loop over p
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/GJI_EDW2_ker.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/GJI_EDW2_ker.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/GJI_EDW2_ker.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -0,0 +1,633 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# GJI_EDW2_ker.pl
+# Carl Tape
+# 26-May-2009
+#
+# This script reads in a horizontal cross-section data file with five columns:
+# lon, lat, model-1, model-2, ln(model-2 / model-1)
+# and plots the three cross-sections.
+#
+# GJI_EDW2_ker.pl -40/5 1 0.10 1.5 1 0.75
+#
+#==========================================================
+
+if (@ARGV < 6) {
+ die("Usage: GJI_EDW2.pl xxx\n");
+}
+($zran,$ivel,$cpert2,$vexag,$irun,$heightxc0) = @ARGV;
+
+$icolor = 1;
+$ivertlab = 0; # vertical exaggeration label
+$isimbox = 1; # box around simulation region
+$iletter = 1; # label for publication figures
+$letter = "a";
+
+$pmin = 11; $pmax = 11;
+$smodel1 = "m00"; $smodel2 = "m16";
+
+#($pmin,$pmax) = split(/\//,$pinds);
+#($smodel1,$smodel2) = split(/\//,$smodels);
+$smodeltag = "${smodel2}_${smodel1}";
+if ($ivel==1) {
+ $modlab = "vs"; $mtit = "Vs"; $vmin = 2.0; $vmax = 4.21; $icolv = 3; $vtick = 1;
+} elsif ($ivel==2) {
+ $modlab = "vb"; $mtit = "Vb"; $vmin = 2.0; $vmax = 5.5; $icolv = 6; $vtick = 1;
+} elsif ($ivel==3) {
+ $modlab = "vp"; $mtit = "Vp"; $vmin = 2.5; $vmax = 7.5; $icolv = 9; $vtick = 1;
+}
+($zmin0,$zmax0) = split(/\//,$zran);
+
+# directory containing the data files
+$stirun = sprintf("%2.2i",$irun);
+$dirdat = "INPUT/vert_${stirun}";
+
+#---------------------------------------------------------
+# absolute origins
+
+$widmap = 2.5; # with of horizontal cross sections
+$height = 2.1; # approximate (depends on projection)
+$dxgap = 0.5;
+$dygap = 0.5;
+$x1 = 1.5; $y1 = 7.5; # Vs m16 horizontal cross section
+$or1 = "-Xa${x1} -Ya${y1}";
+
+$x2 = $x1 + $widmap + $dxgap; $y2 = $y1; $or2 = "-Xa${x2} -Ya${y2}";
+$x3 = $x1; $y3 = $y2 - $ygap - $height; $or3 = "-Xa${x3} -Ya${y3}";
+$x4 = $x2; $y4 = $y3; $or4 = "-Xa${x4} -Ya${y4}";
+$x5 = $x1; $y5 = $y3 - $ygap - $height; $or5 = "-Xa${x5} -Ya${y5}";
+$x6 = $x2; $y6 = $y5; $or6 = "-Xa${x6} -Ya${y6}";
+ at omaps = ($or1,$or2,$or3,$or4,$or5,$or6);
+
+$Dthick = 0.1; # thickness of colorbar
+$xgap = 0.2;
+$ygap = 0.5;
+$xo1 = $x1 - $xgap - $Dthick; $yo1 = $y1 + $ygap; $ocb1 = "-Xa${xo1} -Ya${yo1}";
+$xo2 = $x2 + $widmap + $xgap; $yo2 = $yo1; $ocb2 = "-Xa${xo2} -Ya${yo2}";
+$xo3 = $xo1; $yo3 = $y3 + $ygap; $ocb3 = "-Xa${xo3} -Ya${yo3}";
+$xo4 = $xo2; $yo4 = $yo3; $ocb4 = "-Xa${xo4} -Ya${yo4}";
+$xo5 = $xo1; $yo5 = $y5 + $ygap; $ocb5 = "-Xa${xo5} -Ya${yo5}";
+$xo6 = $xo2; $yo6 = $yo5; $ocb6 = "-Xa${xo6} -Ya${yo6}";
+ at ocbs = ($ocb1,$ocb2,$ocb3,$ocb4,$ocb5,$ocb6);
+
+# vertical cross section
+$xvert = $x1 + 0.0;
+$yvert = $y1 + 1.9;
+$overt = "-Xa${xvert} -Ya${yvert}";
+
+#---------------------------------------------------------
+
+# subtitles
+ at titles = ("$mtit $smodel1","$mtit $smodel2","ln(${smodel2} / ${smodel1})");
+
+# parameters for color scales
+#$cpert1 = 0.2;
+#$cpert2 = 0.2;
+$scale_color = 65.0; # 65.0
+$colorbar = "seis";
+
+# KEY: resolution of color plots
+$interp = "-I0.5/0.5 -S1.5"; # nearneighbor
+#$interp = "-I1/1 -S2"; # nearneighbor
+#$interp = "-I1.0"; # xyz2grd
+$grdfile = "temp.grd";
+
+ at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+
+$dir0 = "/home/carltape/gmt";
+$dir = "socal_2005";
+
+$fault_file = "$dir0/faults/jennings_more.xy";
+$fault_file2 = "$dir0/faults/scitex.lin";
+#$kcf_file = "$dir0/faults/kcf.xy";
+#$breck_file = "$dir0/faults/breck.xy";
+
+$plate_file = "$dir0/plates/plate_boundaries/bird_boundaries";
+$plate_labels = "$dir0/socal_2005/socal_plate_labs.xyz";
+
+# continental shelf
+$ishelf = 0;
+$shelf_file = "/home/carltape/gmt/coastlines/socal_shelf";
+$Wshelf = "-W0.75p,0/0/0,--";
+
+# CMT sources
+$dir_sources = "/home/carltape/results/SOURCES/socal_16";
+$outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
+$inner_boundary = "${dir_sources}/SPECFEM_inner_points.dat";
+
+# cities
+$cities = "$dir0/cities/socal_tomo_cities";
+
+# socal map
+#$tick1 = 100; $tick2 = 100; # no ticks for map
+$tick1 = 1; $tick2 = 0.5;
+$tlen_map = "3p";
+#$Rmap = "-R-122/-114/32/37";
+$Rmap = "-R-119/-115/33/35.2";
+
+$topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
+$fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
+$coast_res = "-Df -A0/0/4";
+
+$B = "$B0".$Bopts[0];
+
+$iportrait = 1;
+if ($iportrait == 1) {
+ $orient = "-P"; $rotangle = 0;
+} else {
+ $orient = " "; $rotangle = 90;
+}
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "12";
+$fsize2 = "9";
+$fsize3 = "6";
+$fontno = "1"; # 1 or 4
+$tlen = "6p";
+$fpen = "1p";
+$tpen = "1p";
+
+# plotting specificiations
+# NOTE: pen attribute (-W) does not work for psvelo
+#$coast_info = "-A5000 -Dl -W1.0p -S220/255/255 -C220/255/255 -G200"; # A : smallest feature plotted, in km^2; D : resolution
+$coast_res = "-Df -A0/0/4";
+$coast_info = "${coast_res} -W0.75p -Na/0.75p";
+#$coast_info = "${coast_res} -W1p -Na/1.0p -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
+$coast_info2 = "${coast_res} -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
+$fault_info = "-M -W0.75p,0/0/0";
+
+$sky_color = "100/100/100";
+$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
+
+#$colinfo = "-W1p,0/0/0 -G255/255/255";
+$colinfo_map = "-W1.0p,0/0/0 -G0/255/255";
+$srcinfo_map = "-Sa15p $colinfo_map -N";
+$recinfo_map = "-Si15p $colinfo_map -N";
+$colinfo_xc = "-W1.0p,0/0/0 -G0/255/255";
+$srcinfo_xc = "-Sa15p $colinfo_xc -N";
+$recinfo_xc = "-Si15p $colinfo_map -N";
+
+$textinfo = "-G255 -S1.5p";
+$textinfo2 = "-G0/0/0 -S2p,255/255/255";
+$textinfo3 = "-G0/0/0 -S2p,255/255/0";
+
+# BOOLEAN: WHICH FIGURES TO PLOT
+#$ixv = 1; $ipdf = 0;
+$ixv = 0; $ipdf = 1;
+$ijpg = 0;
+
+# plot title
+$J_title = "-JM1";
+$R_title = "-R0/1/0/1";
+$fsize_title = $fsize0;
+
+#---------------------------------------------------------
+
+# load values and check that the files exist
+($nump,undef,undef) = split(" ",`ls -1 ${dirdat}/*path | wc`);
+if ($nump == 0) {
+ die("No ray paths available");
+}
+if ($pmax > $nump) {
+ $pmax = $nump;
+}
+if ($pmin < 1) {
+ $pmin = 1;
+}
+print "\nNumber of cuts is $nump\n";
+
+# load fault positions along each profile
+#$ffile = "$dirdat/ALL_${stirun}_rays_fault_positions";
+$ffile = "$dirdat/ALL_${stirun}_rays_fault_positions_mod";
+if (not -f $ffile) {
+ die("Check if ffile $ffile exist or not\n");
+}
+open(IN,$ffile); @faults = <IN>; close(IN);
+
+$cshfile = "GJI_EDW2.csh";
+open(CSH,">$cshfile");
+print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter MEASURE_UNIT inch BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tlen LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1 FRAME_PEN $fpen TICK_PEN $tpen \n";
+
+#@iBs = (9,3,14,6,14,6); # what sides of each cross section to show tick marks
+ at iBs = (7,4,7,4,7,4);
+
+# modify the height of each cross-section based on the length
+$p = 11;
+$stip = sprintf("%3.3i",$p);
+
+# datafile -- COLUMNS: lon, lat, m00, m16
+$dfile = "${dirdat}/vert_${stirun}_xc_vsvbvp_${smodeltag}_${stip}.dat";
+if (not -f $dfile) {
+ die("Check if dfile $dfile exist or not\n");
+}
+print "Data file is $dfile\n";
+
+# source and receiver points -- NOT necessarily the endpoints of the ray path
+# COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
+$rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
+if (not -f $rfile) {
+ die("Check if rfile $rfile exist or not\n");
+}
+open(IN,$rfile); @line = <IN>; close(IN);
+($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
+#print "-- $slon -- $slat -- $sdist -- $sdep --\n";
+#print "$rlon -- $rlat -- $rdist -- $rdep -- $raz -- $iflip--\n";
+
+# positions of six faults along each profile
+# (1) SAF, (2) GF, (3) SGF, (4) MCF, (5) SYF, (6) CRF, (7) EF, (8) SN, (9) KC
+(undef,$eid,$stanet,$azi1,$azi2,$saf,$gf,$sgf,$mcf,$syf,$crf,$ef,$snf,$kcf) = split(" ",$faults[$p-1]);
+($sta,$net) = split("\\.",$stanet); # only for irun = 1
+$sazi1 = sprintf("%3.3i",$azi1);
+$sazi2 = sprintf("%3.3i",$azi2);
+print "-- $stanet -- $sta -- $net -- $eid -- $sazi1 -- $sazi2 --\n";
+
+# axes limits for cross sections
+($dmin,$dmax,$zmin,$zmax) = split(" ",`minmax -C -I2 $dfile`);
+$dran = $dmax - $dmin; $zran = $zmax - $zmin;
+$heightxc = $heightxc0;
+
+# compute the height, including vertical exaggeration
+$wid = $heightxc*$dran/$zran / $vexag;
+$dfile = $dfile;
+$R = "-R$dmin/$dmax/$zmin0/$zmax0";
+$J = "-JX$wid/$heightxc";
+if ($iflip==1) {
+ $J = "-JX-$wid/$heightxc";
+}
+
+#$widmap = $heightxc*(1.9/1.5); # width of map -- tied to height of xc
+$Jmap = "-JM${widmap}";
+$B0map = "-Ba${tick1}f${tick2}d:.\" \":";
+#$Bmap = "$B0map".$Bopts[7];
+$Bmap = "$B0map".$Bopts[15];
+
+#---------------------------------------------------------
+# USER INPUT
+
+# parameters controlling the placement of subplots and labels
+$x0_map = 1;
+$y0_map = 8;
+$x0_xc = $x0_map + 3.5; $y0_xc = $y0_map;
+$ygap_lab = 0.1; $xgap_lab = 0.6*$ygap_lab; # labels on the xc
+$xgap_cbar = 0.2; # gap between plot and cbar
+
+# factors (times $heightxc) controlling y position of colorbar and labels
+$fygap_bot = 0.0;
+$f_cbar = 0.9;
+$fygap_top = 0.0;
+$f_lab = $fygap_bot + $f_cbar + $fygap_top;
+$Dlen = $heightxc*$f_cbar; # length of colorbar
+$Dscale0 = "-D0/0/${Dlen}/${Dthick}";
+
+#---------------------------------------------------------
+
+# scale for plots
+$tick1x = 50; $tick2x = 10;
+$tick1z = 20; $tick2z = 10;
+ at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+$B0 = sprintf("-Ba%3.3ff%3.3fd:\" \":/a%3.3ff%3.3fd:\" \"::.\" \":WesN",$tick1x,$tick2x,$tick1z,$tick2z);
+#$B0 = sprintf("-Ba%3.3ff%3.3fg%3.3f:\" \":/a%3.3ff%3.3fg%3.3f:\" \"::.\" \":WesN",$tick1x,$tick2x,$tick1x,$tick1z,$tick2z,$tick1z); # with gridlines
+
+#==================================================
+
+#-------------------
+# absolute origins -- these depend on the WIDTH of the cross section
+
+# overall label
+$xgap = 0.1; $ygap = -0.1;
+$x0_label = $xvert + $xgap;
+$y0_label = $yvert + $heightxc0 + $ygap;
+$olab = "-Xa${x0_label} -Ya${y0_label}";
+
+# colorbar
+$fac = $heightxc*$fygap_bot + $Dlen/2;
+$xcb1 = $xvert + $wid + $xgap_cbar; $ycb1 = $yvert + $fac;
+$ocb = "-Xa${xcb1} -Ya${ycb1}";
+
+# label for colorbar
+$xcblab1 = $xcb1 + 0.6; $ycblab1 = $ycb1 + 0.5;
+$ocblab1 = "-Xa${xcblab1} -Ya${ycblab1}";
+$xcblab2 = $xcblab1; $ycblab2 = $ycblab1 - 0.2;
+$ocblab2 = "-Xa${xcblab2} -Ya${ycblab2}";
+
+#-------------------
+
+#$dtitle = sprintf("Depth = %.1f km",$zdep/1000);
+
+$cpert1 = $vpert;
+
+# colorpoint file for m00 and m16 -- VARIES with depth -- in KILOMETERS
+$cptfile1 = "color1.cpt";
+$dc = ($vmax-$vmin)/${scale_color};
+$T = sprintf("-T%3.3e/%3.3e/%3.3e",$vmin,$vmax,$dc);
+print CSH "makecpt -C$colorbar $T -D -Z > color.cpt\n"; # -Z for continuous (not for pscontour)
+print CSH "sed 's/^F.*/F 200 200 200/' color.cpt > $cptfile1\n";
+
+# colorbars for cross sections
+$Bscale1 = sprintf("-B%2.2ef0.5:\" \": -E8p",$vtick);
+
+#==================================================
+
+$stip = sprintf("%3.3i",$p);
+$fname = "GJI_EDW2_ker";
+$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+#---------------------
+# 1: VERTICAL CROSS SECTION OF MODEL
+
+# datafile -- COLUMNS: lon, lat, m00, m16
+$dfile = "${dirdat}/vert_${stirun}_xc_vsvbvp_${smodeltag}_${stip}.dat";
+if (not -f $dfile) {
+ die("Check if dfile $dfile exist or not\n");
+}
+print "Data file is $dfile\n";
+
+# source and receiver points -- NOT necessarily the endpoints of the ray path
+# COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
+$rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
+if (not -f $rfile) {
+ die("Check if rfile $rfile exist or not\n");
+}
+open(IN,$rfile); @line = <IN>; close(IN);
+($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
+
+# axes limits for cross sections
+($dmin,$dmax,$zmin,$zmax) = split(" ",`minmax -C -I2 $dfile`);
+$dran = $dmax - $dmin; $zran = $zmax - $zmin;
+
+# compute the height, including vertical exaggeration
+#$heightxc = $zran*$wid/$dran * $vexag;
+$wid = $heightxc*$dran/$zran / $vexag;
+$R = "-R$dmin/$dmax/$zmin0/$zmax0";
+$J = "-JX$wid/$heightxc";
+if ($iflip==1) {
+ $J = "-JX-$wid/$heightxc";
+}
+
+#-----------------------------------------------
+
+$k = 2;
+
+#$shift = $shifts[$k-1];
+$B = "$B0".$Bopts[5];
+$title = $titles[$k-1];
+
+#------------------
+# get absolute origins for various pieces
+
+$Dscale = $Dscale0;
+
+#------------------
+
+print CSH "psbasemap $J $R $B -K -V $overt $orient > $psfile\n"; # START
+print CSH "psbasemap $J $R $B -G${sky_color} -K -O -V $overt >> $psfile\n";
+
+if ($icolor==1) {
+ print CSH "awk '{print \$1,\$2,\$($icolv+$k-1)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$4}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1 $J -Q -K -O -V $overt >> $psfile\n";
+}
+
+print CSH "psscale -C${cptfile1} $Dscale $Bscale1 -K -O -V $ocb >> $psfile\n";
+$slab1 = "Vs m16";
+print CSH "pstext -N $R_title $J_title -K -O -V $ocblab1 >>$psfile<<EOF\n0 0 $fsize2 0 $fontno CM $slab1\nEOF\n";
+$slab2 = "km/s";
+print CSH "pstext -N $R_title $J_title -K -O -V $ocblab2 >>$psfile<<EOF\n0 0 $fsize2 0 $fontno CM $slab2\nEOF\n";
+
+# horizontal line
+$finfo = "-W1.0p,0/0/0,--";
+$zdep = -4; $x1 = -500; $x2 = 500;
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$x1 $zdep\n$x2 $zdep\nEOF\n";
+
+# vertical lines
+$zminf0 = -5;
+$finfo = "-W1.5p,0/0/0";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$saf $zmax0\n$saf $zminf0\nEOF\n";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$gf $zmax0\n$gf $zminf0\nEOF\n";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$sgf $zmax0\n$sgf $zminf0\nEOF\n";
+#print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$mcf $zmax0\n$mcf $zminf0\nEOF\n";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$syf $zmax0\n$syf $zminf0\nEOF\n";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$crf $zmax0\n$crf $zminf0\nEOF\n";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$ef $zmax0\n$ef $zminf0\nEOF\n";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$snf $zmax0\n$snf $zminf0\nEOF\n";
+print CSH "psxy $J $R $finfo -K -O -V $overt >>$psfile<<EOF\n$kcf $zmax0\n$kcf $zminf0\nEOF\n";
+
+# plot source and receiver
+if ($irun == 1) {
+ print CSH "awk '{print \$3,-\$4}' $rfile | psxy $J $R $srcinfo_xc -K -O -V $overt >> $psfile\n";
+ print CSH "awk '{print \$7,$zmax0}' $rfile | psxy $J $R $recinfo_xc -K -O -V $overt >> $psfile\n";
+}
+
+# details for specific cross-sections
+$talign = "CB";
+$flsize = 9;
+$textinfo = "-N -C3p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+if ($saf != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$saf $zmax0 $flsize 0 $fontno $talign SA\nEOF\n";
+}
+if ($gf != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$gf $zmax0 $flsize 0 $fontno $talign G\nEOF\n";
+}
+if ($sgf != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$sgf $zmax0 $flsize 0 $fontno $talign SG\nEOF\n";
+}
+#if($mcf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$mcf $zmax0 $flsize 0 $fontno $talign MC\nEOF\n";}
+if ($syf != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$syf $zmax0 $flsize 0 $fontno $talign SY\nEOF\n";
+}
+if ($crf != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$crf $zmax0 $flsize 0 $fontno $talign CR\nEOF\n";
+}
+if ($ef != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$ef $zmax0 $flsize 0 $fontno $talign E\nEOF\n";
+}
+if ($snf != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$snf $zmax0 $flsize 0 $fontno $talign SN\nEOF\n";
+}
+if ($kcf != 9999) {
+ print CSH "pstext $J $R $textinfo -K -O -V $overt >>$psfile<<EOF\n$kcf $zmax0 $flsize 0 $fontno $talign KC\nEOF\n";
+}
+
+# plot overall label (for publication)
+if ($iletter==1) {
+ $textinfo = "-N -C3p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 14 0 $fontno TL $letter\nEOF\n";
+}
+
+# # plot label on each cross section
+# for ($k = 2; $k <= 2; $k ++ ) {
+# $title = $titles[$k-1];
+# $talign = "LB";
+# #$otitle = "-Xa${x3} -Ya${y3}";
+# $otitle = $overtlabs[$k-1];
+# $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+# print CSH "pstext $R_title $J_title $textinfo -K -O -V $otitle >>$psfile<<EOF\n0 0 12 0 $fontno $talign $title\nEOF\n";
+# }
+
+#---------------------
+# 2: HORIZONTAL CROSS SECTIONS OF MODEL
+
+$icut = 13; # depth to plot
+$sticut = sprintf("%3.3i",$icut);
+$irunh = 2; $stirunh = sprintf("%2.2i",$irunh);
+$dirdath = "INPUT/horz_${stirunh}";
+
+# KEY: file showing the cuts and the plotting range for the values
+$fcuts = "${dirdath}/horz_${stirunh}_xc_cuts_mod";
+if (not -f $fcuts) {
+ die("Check if fcuts $fcuts exist or not\n");
+}
+open(IN,$fcuts); @lines = <IN>; $nump = @lines;
+
+# factors (times $heightxc) controlling y position of colorbar and labels
+$Dthick = 0.1; # thickness of colorbar
+$Dlen = 1.0; # length of colorbar
+$Dscale = "-D0/0/${Dlen}/${Dthick}";
+
+#---------------------
+# 3: HORIZONTAL CROSS SECTIONS OF KERNELS
+
+# parameters controlling which kernel slices to plot
+ at icols = (4,4, 3,3,4,5); # column in the datafile
+ at icuts = (1,13, 1,13,13,13); # index into depth slice
+ at letters = ("b","c","d","e","f","g");
+
+# KEY: data files
+ at dtags = ("${modlab}_${smodeltag}","${modlab}_${smodeltag}","ker_h000km_v000km","ker_h000km_v000km","ker_h003km_v001km","ker_h003km_v001km");
+
+# color scales
+ at norms = (1e-3,1e-3, 1e12,1e12,1e12,1e12);
+ at kvals = (99,99, 1.21,0.81,1.01,0.81); # max value for color scale
+ at labels = ("Vs m16, depth 0 km","Vs m16, depth 4 km","K1, depth 0 km","K1, depth 4 km","K2, depth 4 km","K3, depth 4 km");
+ at clabels1 = (" "," ","10\@+-12\@+","10\@+-12\@+","10\@+-12\@+","10\@+-12\@+",);
+$uker = "m\@+-3\@+ s"; # BDK units
+ at clabels2 = ("km/s","km/s",$uker,$uker,$uker,$uker);
+ at tick1s = (0.4,0.2,0.8,0.8,0.8,0.8);
+ at tick2s = (0.2,0.2,0.2,0.2,0.2,0.2);
+ at Aflip = ("-A"," ","-A"," ","-A"," ",);
+ at xtxs = (-0.1,0.2,-0.1,0.2,-0.1,0.2);
+
+#for ($ik = 1; $ik <= 5; $ik ++ ) {
+for ($ik = 1; $ik <= 6; $ik ++ ) {
+
+ $origin = $omaps[$ik-1];
+ $iB = $iBs[$ik-1];
+ $Bmap = "$B0map".$Bopts[$iB];
+ $Bscale1 = sprintf("-B%2.2ef%2.2e:\" \": -E8p %s",$tick1s[$ik-1],$tick2s[$ik-1],$Aflip[$ik-1]);
+
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -O -V $origin >> $psfile\n";
+
+ $dtag = $dtags[$ik-1];
+ $norm = $norms[$ik-1];
+ $icol = $icols[$ik-1];
+ $cmax = $kvals[$ik-1];
+ $icut = $icuts[$ik-1]; $sticut = sprintf("%3.3i",$icut);
+
+ # minmax for color scale
+ if ($ik <= 2) {
+ # load values and check that the files exist
+ (undef,$zcut,$vs_norm,$vb_norm,$vpert) = split(" ",$lines[$icut-1]);
+ $cnorm = $vs_norm;
+ $cpert = $vpert;
+ if ($ik==2) {
+ $cpert = 0.10;
+ }
+ $cmin = $cnorm*exp(-$cpert)*$norm;
+ $cmax = $cnorm*exp($cpert)*$norm;
+
+ } else {
+ $cmin = -$cmax;
+ }
+
+ # data file
+ if($ik <= 2) {
+ $dfileh = "${dirdath}/horz_${stirunh}_xc_${dtag}_${sticut}_mask1.dat";
+ } else {
+ $dfileh = "${dirdath}/horz_${stirunh}_xc_${dtag}_${sticut}.dat";
+ }
+ if (not -f ${dfileh}) {
+ die("Check if dfileh ${dfileh} exist or not\n");
+ }
+
+ # colorpoint file
+ $cptfileh = "colorK.cpt";
+ $dc = ($cmax-$cmin)/${scale_color};
+ $T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
+ print CSH "makecpt -C$colorbar $T -D -Z > $cptfileh\n";
+
+ if ($icolor==1) {
+ $interph = "-I1m -S2m";
+ $grdfileh = "temph.grd";
+ #print CSH "awk '{print \$1,\$2,\$${icol}*1e12}' $dfileh | xyz2grd -G$grdfile $Rmap $interph\n";
+ print CSH "awk '{print \$1,\$2,\$${icol}*$norm}' $dfileh | nearneighbor -G$grdfileh $Rmap $interph\n";
+ print CSH "grdimage $grdfileh -C$cptfileh $Jmap -Q -K -O -V $origin >> $psfile\n";
+ }
+
+ print CSH "pscoast $Jmap $Rmap $coast_info -K -O -V $origin >> $psfile \n";
+ if ($ishelf==1) {
+ print CSH "awk '{print \$2,\$1}' ${shelf_file} | psxy $Jmap $Rmap $Wshelf -K -O -V $origin >> $psfile\n";
+ }
+
+
+# colorbar
+$ocb = $ocbs[$ik-1];
+
+# color scale
+print CSH "psscale -C${cptfileh} $Dscale $Bscale1 -K -O -V $ocb >> $psfile\n";
+$slab1 = $clabels1[$ik-1];
+print CSH "pstext -N $R_title $J_title -K -O -V $ocb >>$psfile<<EOF\n$xtxs[$ik-1] 1.0 $fsize2 0 $fontno CM $slab1\nEOF\n";
+$slab2 = $clabels2[$ik-1];
+print CSH "pstext -N $R_title $J_title -K -O -V $ocb >>$psfile<<EOF\n$xtxs[$ik-1] 0.8 $fsize2 0 $fontno CM $slab2\nEOF\n";
+
+ print CSH "psxy ${fault_file} $Jmap $Rmap $fault_info -K -V -O $origin >> $psfile \n";
+ #print CSH "psxy ${kcf_file} $Jmap $Rmap $fault_info -K -V -O $origin >> $psfile \n";
+
+ # plot ray path
+ $pfile = "${dirdat}/vert_${stirun}_xc_${stip}_ray_path";
+ if (not -f $pfile) {
+ die("Check if pfile $pfile exist or not\n");
+ }
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W1.0p,0/0/0,-- -K -O -V $origin >> $psfile\n";
+ #print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W6p,0/0/0 -K -O -V $origin >> $psfile\n";
+ #print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W3p,0/255/255 -K -O -V $origin >> $psfile\n";
+
+ # plot source and receiver lon-lat
+ if ($irun==1) {
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
+ if (not -f $rfile) {
+ die("Check if rfile $rfile exist or not\n");
+ }
+ print CSH "awk '{print \$1,\$2}' ${rfile} | psxy $Jmap $Rmap $srcinfo_map -K -O -V $origin >> $psfile\n";
+ print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmap $Rmap $recinfo_map -K -O -V $origin >> $psfile\n";
+ }
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -O -V $origin >> $psfile\n";
+
+ # plot overall label (for publication)
+ if ($iletter==1) {
+ $letter = @letters[$ik-1];
+ #$textinfo = "-N -C3p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ #print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 14 0 $fontno TL $letter\nEOF\n";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $origin >>$psfile<<EOF\n0.1 1.6 14 0 $fontno TL $letter\nEOF\n";
+ $label = @labels[$ik-1];
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $origin >>$psfile<<EOF\n2.4 1.6 10 0 $fontno TR $label\nEOF\n";
+ }
+
+}
+
+#------------------------------------
+
+print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n"; # FINISH
+#if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
+if ($ixv==1) {print CSH "gv $psfile &\n";}
+if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+if ($ijpg==1) {print CSH "convert $psfile $jpgfile\n";}
+
+print CSH "gv $psfile &\n";
+
+#==================================================
+close (CSH);
+system("csh -f $cshfile");
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/GJI_EDW2_ker.pl
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl 2010-02-10 00:30:55 UTC (rev 16251)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_one.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -11,7 +11,7 @@
# and plots the m16 cross-sections only.
#
# FIGURES:
-# plot_horz_models_one.pl 11/11 m00/m16 1 0/0/1/0.2 # Vs
+# plot_horz_models_one.pl 21/21 m00/m16 1 0/0/1/0.2 # Vs
#
#==========================================================
@@ -33,16 +33,21 @@
$letter = "b";
$iextra = 1; # extra information on figures
+$irelief = 1; # shaded relief
+
$cgray = "-G220/220/220";
#---------------------------------------------------------
# USER INPUT
+$irun = 2;
+$stirun = sprintf("%2.2i",$irun);
+
# directory containing the data files
-$dirdat = "INPUT/horz_01";
+$dirdat = "INPUT/horz_${stirun}";
# KEY: file showing the cuts and the plotting range for the values
-$fcuts = "${dirdat}/horz_xc_elevation_cuts_mod";
+$fcuts = "${dirdat}/horz_${stirun}_xc_cuts_mod";
if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
# load values and check that the files exist
@@ -54,7 +59,7 @@
$vsnorms[$p-1] = $vs_norm;
$vbnorms[$p-1] = $vb_norm;
$zcuts[$p-1] = -$zcut;
- $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}.dat";
+ $dfiles[$p-1] = "${dirdat}/horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}_mask1.dat";
#if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
}
print "\nNumber of cuts is $nump\n";
@@ -73,14 +78,13 @@
# resolution of color plots
#$interp = "-I2m/2m -S4m"; # key information
-$interp = "-I2m";
$interp = "-I1.5m";
$grdfile = "temp.grd";
# position of titles and labels
$xtx1 = 0.5; $ytx1 = 1.00;
$xtx2 = -0.2; $ytx2 = 0.4;
-#$xtx3 = 1.0; $ytx3 = 0.6;
+$xtx3 = 1.04; $ytx3 = 0.7;
$xtx4 = 1.04; $ytx4 = 0.6;
$xtx5 = -0.15; $ytx5 = 0.85; # A, B, C, etc
$xtx6 = 0.0; $ytx6 = 0.05; # m00, m16, ln(m16/m00)
@@ -124,7 +128,7 @@
if ($iregion==1) {
# southern California
- $wid = 4; # width of figure (inches)
+ $wid = 5.5; # width of figure (inches)
$xmin = -122; $xmax = -114; $ymin = 31.5; $ymax = 37;
$xmin = -121.2; $xmax = -114.8; $ymin = 32.3; $ymax = 36.7; # SPECFEM
$tick1 = 1; $tick2 = 0.5;
@@ -135,7 +139,7 @@
$topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
$fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
$name = "plot_horz_models_one";
- $origin = "-X2i -Y4i";
+ $origin = "-X1.5i -Y4i";
# inset map
$iinset = 1;
@@ -165,7 +169,7 @@
else {$orient = " "; $rotangle = 90}
# color bar
-$Dwid = 0.15;
+$Dwid = 0.20;
$Dlen = $wid*0.5;
$Dx = $wid + 0.25;
$Dy = $Dlen/2;
@@ -175,22 +179,22 @@
$Dscale2 = "-D$Dx/$Dy2/$Dlen/${Dwid}";
# plotting specifications
-$fsize0 = "16";
-$fsize1 = "12";
-$fsize2 = "11";
-$fsize3 = "8";
+$fsize0 = "18";
+$fsize1 = "16";
+$fsize2 = "14";
+$fsize3 = "10";
$fontno = "1"; # 1 or 4
$tick = "6p";
-$fpen = "1p";
-$tpen = "1p";
+$fpen = "2p";
+$tpen = "2p";
# plotting specificiations
# NOTE: pen attribute (-W) does not work for psvelo
#$coast_info = "-A5000 -Dl -W1.0p -S220/255/255 -C220/255/255 -G200"; # A : smallest feature plotted, in km^2; D : resolution
-$coast_info = "${coast_res} -W0.75p -Na/0.75p";
+$coast_info = "${coast_res} -W1.0p -Na/1.0p";
$fault_info = "-M -W0.75p,255/0/0";
-$fault_infoK = "-M -W0.75p,0/0/0";
+$fault_infoK = "-M -W1.5p,0/0/0";
$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
@@ -199,8 +203,10 @@
$textinfo3 = "-G0/0/0 -S2p,255/255/0";
# BOOLEAN: WHICH FIGURES TO PLOT
-$ixv = 0; $ipdf = 1;
-#$ixv = 1; $ipdf = 0;
+$ixv = 0;
+$ipdf = 0;
+$ieps = 1;
+$ijpg = 1;
# plot title
$J_title = "-JM${wid}";
@@ -219,8 +225,10 @@
for ($p = $pmin; $p <= $pmax; $p ++ ) {
$stip = sprintf("%3.3i",$p);
- $fname = "horz_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}_single";
- $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+ $fname = "horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}_extra${iextra}_single";
+ $psfile = "$fname.ps";
+ $epsfile = "$fname.eps";
+ $jpgfile = "$fname.jpg";
$vsnorm = $vsnorms[$p-1]; # reference velocity for m00 and m16
$vbnorm = $vbnorms[$p-1]; # reference velocity for m00 and m16
@@ -287,6 +295,17 @@
print CSH "psbasemap $J $R $B -K -O -V $shift >> $psfile\n";
}
print CSH "psbasemap $J $R $B $cgray -K -O -V >> $psfile\n";
+
+if ($irelief==1) {
+ $grdfile0 = "/home/carltape/datalib/Topography/GLOBAL/ETOPO1/ETOPO1_Ice_g.grd";
+ $grddat = "ETOPO.xyz";
+ $grdnew = "ETOPO.grd";
+ $gradnew = "ETOPO.grad";
+ print CSH "grd2xyz $grdfile0 $R > $grddat\n";
+ print CSH "xyz2grd $grddat $R $interp -G$grdnew\n";
+ print CSH "grdgradient $grdnew -A110 -G$gradnew -Nt0.8 -V\n";
+}
+
if ($icolor==1) {
#print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
#print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
@@ -294,14 +313,27 @@
#print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
print CSH "awk '{print \$1,\$2,\$4/1000}' $dfile | xyz2grd -G$grdfile $R $interp\n";
- print CSH "grdimage $grdfile -C$cptfile1b $J -Q -K -O -V >> $psfile\n";
+if ($irelief==1) {
+ print CSH "grdimage $grdfile -C$cptfile1b -I$gradnew $J -Q -K -O -V >> $psfile\n";
+} else {
+ print CSH "grdimage $grdfile -C$cptfile1b $J -Q -K -O -V >> $psfile\n";
+}
+
#print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
print CSH "psscale -C${cptfile1b} $Dscale2 $Bscale1b -K -O -V >> $psfile\n";
- $slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000); $slab2 = "$mtit km/s";
- #print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
- print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 14 0 $fontno LM $slab2\nEOF\n";
+# $slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000);
+# $slab2 = "$mtit km/s";
+# print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
+# print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 14 0 $fontno LM $slab2\nEOF\n";
+ $slab1 = "$mtit km/s";
+ $slab2 = sprintf("(%.2f \261 %.0f \045)",$cnorm/1000,$cpert1*100);
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize1 0 $fontno LM $slab1\nEOF\n";
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 $fsize2 0 $fontno LM $slab2\nEOF\n";
+
+
+
}
print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile \n";
if($ishelf==1) {print CSH "awk '{print \$2,\$1}' ${shelf_file} |psxy $J $R $Wshelf -K -O -P -V >> $psfile\n";}
@@ -336,19 +368,19 @@
$finfo2 = "-W1.5p,0/255/255";
# plot a line at -119 longitude
- if ($p==11 && $ivs==1) {
+ if ($p==21 && $ivs==1) {
$xmark = -119;
print CSH "psxy $J $R $finfo1 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
print CSH "psxy $J $R $finfo2 -K -O -V >>$psfile<<EOF\n$xmark $ymin\n$xmark $ymax\nEOF\n";
}
# plot two cross section ray paths
- if ($p==21 && $ivs==1) {
+ if ($p==41 && $ivs==1) {
@irs = (2,5); $nray = @irs ;
for ($ik = 1; $ik <= $nray; $ik ++ ) {
$ir = $irs[$ik-1];
$stir = sprintf("%3.3i",$ir);
- $pfile = "./INPUT/vert_01/vert_xc_${stir}_ray_path";
+ $pfile = "./INPUT/vert_01/vert_01_xc_${stir}_ray_path";
if (not -f $pfile) {
die("Check if pfile $pfile exist or not\n");
}
@@ -358,10 +390,16 @@
}
# fault labels
- $fault_labs = "$dir0/faults/socal_fault_labs_red.lonlat";
- $textinfo = "-N -C2p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
- print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfo $J $R -K -V -O >> $psfile\n";
+ #$fault_labs = "$dir0/faults/socal_fault_labs_red.lonlat";
+ $textinfof = "-N -C2p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ #print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfof $J $R -K -V -O >> $psfile\n";
+ print CSH " pstext $textinfof $J $R -K -V -O >> $psfile <<EOF
+-120.7000 36.2000 10 -45 1 CM SA
+-116.2000 33.8000 10 -45 1 CM SA
+-118.1678 35.1562 10 30 1 CM G
+EOF\n";
+
}
#--------------------------------
@@ -387,7 +425,7 @@
if ($iinsetlab==1) {
$lab = $labs[$k-1];
$textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
- print CSH "pstext $R_title $J_title $textinfo -K -O -V >>$psfile<<EOF\n$xtx6 $ytx6 14 0 $fontno LB $lab\nEOF\n";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V >>$psfile<<EOF\n$xtx6 $ytx6 $fsize0 0 $fontno LB $lab\nEOF\n";
}
# plot vertical title left of the row
@@ -404,9 +442,10 @@
} # loop over k
print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n $x_title $y_title 16 0 $fontno CM \nEOF\n"; # FINISH
- #if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
if($ixv==1) {print CSH "gv $psfile &\n";}
if($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+ if ($ieps==1) {print CSH "ps2eps -f -l $psfile\n";}
+ if ($ijpg==1) {print CSH "convert $epsfile $jpgfile\n";}
} # loop over p
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl 2010-02-10 00:30:55 UTC (rev 16251)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_horz_models_three.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -15,7 +15,8 @@
# plot_horz_models_three.pl 1/50 m00/m16 0 1/1/1/0.3 # Vb
# plot_horz_models_three.pl 1/6 m00/m01 1 1/1/1/0.2 # Vs
#
-# plot_horz_models_three.pl 3/11/21 m00/m16 1 1/0/0/0.2 # Vs
+# plot_horz_models_three.pl 3/11/21 m00/m16 1 1/0/0/0.2 # Vs -- irun = 1
+# plot_horz_models_three.pl 5/21/41 m00/m16 1 1/0/0/0.2 # Vs -- irun = 2
#
#==========================================================
@@ -23,7 +24,7 @@
($pinds,$smodels,$ivs,$parms) = @ARGV;
$icolor = 1;
-$icuts = 1;
+$icuts = 0;
($pmin,$pmid,$pmax) = split(/\//,$pinds);
@pinds = ($pmin,$pmid,$pmax);
@@ -48,10 +49,12 @@
# USER INPUT
# directory containing the data files
-$dirdat = "INPUT/horz_01";
+$stirun = "02";
+#$dirdat = "INPUT/horz_${stirun}";
+$dirdat = "INPUT/horz_${stirun}";
# KEY: file showing the cuts and the plotting range for the values
-$fcuts = "${dirdat}/horz_xc_elevation_cuts_mod";
+$fcuts = "${dirdat}/horz_${stirun}_xc_cuts_mod";
if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
# load values and check that the files exist
@@ -63,7 +66,7 @@
$vsnorms[$p-1] = $vs_norm;
$vbnorms[$p-1] = $vb_norm;
$zcuts[$p-1] = -$zcut;
- $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}.dat";
+ $dfiles[$p-1] = "${dirdat}/horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}_mask1.dat";
#if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
}
print "\nNumber of cuts is $nump0\n";
@@ -88,8 +91,8 @@
# position of titles and labels
$xtx1 = 0.5; $ytx1 = 1.00;
$xtx2 = -0.15; $ytx2 = 0.4;
-$xtx3 = 0.7; $ytx3 = -0.1;
-$xtx4 = $xtx3; $ytx4 = -0.3;
+$xtx3 = 0.7; $ytx3 = -0.19; # cbar 1
+$xtx4 = $xtx3; $ytx4 = -0.1; # cbar 2
$xtx5 = -0.15; $ytx5 = 0.85; # A, B, C, etc
$xtx6 = 0.0; $ytx6 = 0.05; # m00, m16, ln(m16/m00)
@@ -154,7 +157,7 @@
$coast_res = "-Df -A0/0/4";
}
-$dX = $wid + 0.5;
+$dX = $wid + 0.25;
$hwid = $wid/2;
#$R = "-Rg";
$R = "-R$xmin/$xmax/$ymin/$ymax";
@@ -173,10 +176,14 @@
else {$orient = " "; $rotangle = 90}
# color bar
-$Dlen = $wid*0.6;
-$Dx = $Dlen*0.5;
-$Dy = -0.5;
-$Dscale = "-D$Dx/$Dy/$Dlen/0.15h";
+$Dwid = 0.15;
+$Dlen = $wid*0.5;
+$Dx = 0.2 + $Dlen*0.5;
+$Dy = -0.2;
+$Dy2 = $Dy-$Dwid;
+$Dy2 = $Dy;
+$Dscale = "-D$Dx/$Dy/$Dlen/${Dwid}h";
+$Dscale2 = "-D$Dx/$Dy2/$Dlen/${Dwid}h";
# plotting specifications
$fsize0 = "14";
@@ -203,8 +210,8 @@
$textinfo3 = "-G0/0/0 -S2p,255/255/0";
# BOOLEAN: WHICH FIGURES TO PLOT
-$ixv = 0; $ipdf = 1;
-#$ixv = 1; $ipdf = 0;
+#$ixv = 0; $ipdf = 1;
+$ixv = 1; $ipdf = 0;
# plot title
$J_title = "-JM${wid}";
@@ -218,14 +225,15 @@
print CSH "gmtset PAPER_MEDIA letter BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tick LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1 FRAME_PEN $fpen TICK_PEN $tpen \n";
@shifts = ($origin,"-X$dX","-X$dX");
-if($ilonlab==1) {@iBs = (2,14,2);} else {@iBs = (15,9,15);}
+#if($ilonlab==1) {@iBs = (2,14,2);} else {@iBs = (15,9,15);}
+ at iBs = (5,2,2);
for ($w = 1; $w <= $nump; $w ++ ) {
$p = $pinds[$w-1];
$stip = sprintf("%3.3i",$p);
- #$fname = "horz_xc_${modlab}_${smodeltag}_${stip}";
+ #$fname = "horz_${stirun}_xc_${modlab}_${smodeltag}_${stip}";
$fname = "horz_three_${modlab}_${smodel2}_${stips}_cuts${icuts}";
$psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
@@ -253,16 +261,17 @@
$cmin = -$cpert1*1.01; $cmax = $cpert1*1.01;
$dc = ($cmax-$cmin)/${scale_color};
$T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
- #print CSH "makecpt -C$colorbar $T -D > $cptfile1a\n";
print CSH "makecpt -C$colorbar $T -D > $cptfile1a\n";
# colorpoint file for m00 and m16 -- VARIES with depth -- in KILOMETERS
+ # THIS IS FOR THE COLORBAR ONLY
$cptfile1b = "color1b.cpt";
$cmin = $cnorm*exp(-$cpert1)*1e-3;
$cmax = $cnorm*exp($cpert1)*1e-3;
$dc = ($cmax-$cmin)/${scale_color};
$T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
print CSH "makecpt -C$colorbar $T -D > $cptfile1b\n";
+ $Bscale1b = sprintf("-B%2.2ef0.1:\" \": -E10p",$ctick1b_km);
# colorpoint file for perturbation ln(m16/m00) -- FIXED with depth
$cptfile2 = "color2.cpt";
@@ -271,16 +280,14 @@
$T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin,$cmax,$dc);
print CSH "makecpt -C$colorbar $T -D > $cptfile2\n";
- #
+ # ticks for colorbars
$Bscale1a = sprintf("-Ba%2.2ef0.05:\" \": -E10p -A",$cpert1);
- $Bscale1b = sprintf("-B%2.2ef0.1:\" \": -E10p",$ctick1b_km);
- $Bscale2 = sprintf("-Ba%2.2ef0.05:\" \": -E10p -A",$cpert2);
+ $Bscale2 = sprintf("-Ba%2.2ef0.05:\" \": -E10p",$cpert2);
# attempting to get shaded relief superimposed
#$gradfile0 = "/home/datalib/Topography/GLOBAL/ETOPO1/ETOPO1_Ice_g.grad";
#if (not -f $gradfile0) {die("Check if gradfile $gradfile0 exist or not\n");}
-
$shift = $shifts[$w-1];
$iB = $iBs[$w-1];
$B = "$B0".$Bopts[$iB];
@@ -296,12 +303,20 @@
#print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
print CSH "awk '{print \$1,\$2,log(\$4/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
print CSH "grdimage $grdfile -C$cptfile1a $J -Q -K -O -V >> $psfile\n";
+ #print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
+ print CSH "psscale -C${cptfile1b} $Dscale2 $Bscale1b -K -O -V >> $psfile\n";
- print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
- print CSH "psscale -C${cptfile1b} $Dscale $Bscale1b -K -O -V >> $psfile\n";
- $slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000); $slab2 = "$mtit km/s";
+ #$slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000);
+ $slab1 = sprintf("(%.2f \261 %.0f \045)",$cnorm/1000,$cpert1*100);
+ $slab2 = "$mtit km/s";
print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 $fsize2 0 $fontno LM $slab2\nEOF\n";
+
+# print CSH "psscale -C${cptfile1a} $Dscale $Bscale1a -K -O -V >> $psfile\n";
+# print CSH "psscale -C${cptfile1b} $Dscale $Bscale1b -K -O -V >> $psfile\n";
+# $slab1 = sprintf("ln(%s / %.2f)",$mtit,$cnorm/1000); $slab2 = "$mtit km/s";
+# print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx3 $ytx3 $fsize2 0 $fontno LM $slab1\nEOF\n";
+# print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx4 $ytx4 $fsize2 0 $fontno LM $slab2\nEOF\n";
}
print CSH "pscoast $J $R $coast_info -K -O -V >> $psfile \n";
if($ishelf==1) {print CSH "awk '{print \$2,\$1}' ${shelf_file} |psxy $J $R $Wshelf -K -O -P -V >> $psfile\n";}
@@ -324,7 +339,7 @@
$boxinfo = "-W1.5p,0/0/255"; # -A : suppress drawing line segments as great circle arcs
# Lin model (2007) box
- $lin_boundary = "/net/denali/home2/carltape/gmt/tomography/lin_2007/lin_boundary_points.dat";
+ $lin_boundary = "/home/carltape/gmt/tomography/lin_2007/lin_boundary_points.dat";
print CSH "psxy ${lin_boundary} $J $R $boxinfo -K -O -V >>$psfile\n";
}
}
@@ -360,6 +375,19 @@
#--------------------------------
+ # fault labels
+ $fault_labs = "/home/carltape/gmt/faults/socal_fault_labs.lonlat";
+ $boxborder = "0/0/0";
+ $textinfof = "-N -C2p -W255/255/255o,0.5p,${boxborder},solid -G0/0/0";
+ #print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfof $J $R -K -V -O >> $psfile\n";
+ print CSH " pstext $textinfof $J $R -K -V -O >> $psfile <<EOF
+-120.7000 36.2000 10 -45 1 CM SA
+-116.2000 33.8000 10 -45 1 CM SA
+-118.1678 35.1562 10 30 1 CM G
+EOF\n";
+
+ #--------------------------------
+
print CSH "psbasemap $J $R $B -K -V -O >> $psfile\n";
# plot cities
@@ -399,7 +427,7 @@
print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n $x_title $y_title 16 0 $fontno CM \nEOF\n"; # FINISH
#if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
- if($ixv==1) {print CSH "ghostview $psfile &\n";}
+ if($ixv==1) {print CSH "gv $psfile &\n";}
if($ipdf==1) {print CSH "ps2pdf $psfile\n";}
#------------------------------------
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl 2010-02-10 00:30:55 UTC (rev 16251)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_basemap.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -26,7 +26,7 @@
$nump = @pinds;
$stp = sprintf("%3.3i",$nump);
-$irun = 2; # KEY: 1, 2 (SAF), 3 (Garlock)
+$irun = 5; # KEY: 1, 2 (SAF), 3 (Garlock)
$icolor = 1;
$ivertlab = 0; # vertical exaggeration label
$isimbox = 1; # box around simulation region
@@ -160,8 +160,7 @@
print "\nNumber of total cuts is $nump0\n";
# load ALL fault positions along each profile
- #$ffile = "$dirdat/ALL_rays_fault_positions";
- $ffile = "$dirdat/ALL_rays_fault_positions_mod";
+ $ffile = "$dirdat/ALL_${stirun}_rays_fault_positions_mod";
if (not -f $ffile) {die("Check if ffile $ffile exist or not\n");}
open(IN,$ffile); @faults = <IN>; close(IN);
@@ -212,7 +211,7 @@
# source and receiver points -- NOT necessarily the endpoints of the ray path
# COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
- $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
open(IN,$rfile); @line = <IN>; close(IN);
($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
@@ -228,7 +227,7 @@
#---------------------------------------------------------
# plot ray path
- $pfile = "${dirdat}/vert_xc_${stip}_ray_path";
+ $pfile = "${dirdat}/vert_${stirun}_xc_${stip}_ray_path";
if (not -f $pfile) {
die("Check if pfile $pfile exist or not\n");
}
@@ -236,7 +235,7 @@
print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $Jmap $Rmap -W3p,0/255/255 -K -O -V $omap >> $psfile\n";
# plot source and receiver lon-lat
- $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
if (not -f $rfile) {
die("Check if rfile $rfile exist or not\n");
}
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_one.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_one.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_one.pl 2010-02-10 00:31:09 UTC (rev 16252)
@@ -0,0 +1,531 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_vert_models_one.pl
+# Carl Tape
+# 07-Dec-2008
+#
+# This script reads in a horizontal cross-section data file with five columns:
+# lon, lat, model-1, model-2, ln(model-2 / model-1)
+# and plots the three cross-sections.
+#
+# EXAMPLE:
+# plot_vert_models_one.pl 1/86 m00/m16 -40/5 1 0.10 3.0 1 2.2 # Vs -- src-rec paths
+# plot_vert_models_one.pl 4/32 m00/m16 -40/5 1 0.10 3.0 2 2.2 # Vs -- SAF-perp profiles and three more (cut 14,25)
+# plot_vert_models_one.pl 1/10 m00/m16 -40/5 1 0.10 3.0 3 2.2 # Vs -- Garlock-perp profiles
+# plot_vert_models_one.pl 1/13 m00/m16 -40/5 1 0.10 3.0 4 2.2 # Vs -- N-S (constant longitude)
+# plot_vert_models_one.pl 1/9 m00/m16 -40/5 1 0.10 3.0 5 2.2 # Vs -- E-W (constant latitude)
+#
+#
+#==========================================================
+
+if (@ARGV < 8) {die("Usage: plot_vert_models_one.pl xxx\n");}
+($pinds,$smodels,$zran,$ivel,$cpert2,$vexag,$irun,$heightxc0) = @ARGV;
+
+$icolor = 1;
+$ivertlab = 0; # vertical exaggeration label
+$isimbox = 1; # box around simulation region
+$iletter = 0; # label for publication figures
+$letter = "a";
+$imfinal = 1; # plot an additional figure with only m16 and the map
+
+($pmin,$pmax) = split(/\//,$pinds);
+($smodel1,$smodel2) = split(/\//,$smodels);
+$smodeltag = "${smodel2}_${smodel1}";
+if($ivel==1) {$modlab = "vs"; $mtit = "Vs"; $vmin = 2.0; $vmax = 4.21; $icol = 3; $vtick = 1;}
+elsif($ivel==2) {$modlab = "vb"; $mtit = "Vb"; $vmin = 2.0; $vmax = 5.5; $icol = 6; $vtick = 1;}
+elsif($ivel==3) {$modlab = "vp"; $mtit = "Vp"; $vmin = 2.5; $vmax = 7.5; $icol = 9; $vtick = 1;}
+($zmin0,$zmax0) = split(/\//,$zran);
+
+# directory containing the data files
+$stirun = sprintf("%2.2i",$irun);
+$dirdat = "INPUT/vert_${stirun}";
+
+#---------------------------------------------------------
+
+# subtitles
+ at titles = ("$mtit $smodel1","$mtit $smodel2","ln(${smodel2} / ${smodel1})");
+
+# parameters for color scales
+#$cpert1 = 0.2;
+#$cpert2 = 0.2;
+$scale_color = 65.0; # 65.0
+$colorbar = "seis";
+
+# KEY: resolution of color plots
+$interp = "-I0.5/0.5 -S1.5"; # nearneighbor
+#$interp = "-I1/1 -S2"; # nearneighbor
+#$interp = "-I1.0"; # xyz2grd
+$grdfile = "temp.grd";
+
+ at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+
+$dir0 = "/home/carltape/gmt";
+$dir = "socal_2005";
+
+$fault_file = "$dir0/faults/jennings_more.xy";
+$fault_file2 = "$dir0/faults/scitex.lin";
+#$kcf_file = "$dir0/faults/kcf.xy";
+#$breck_file = "$dir0/faults/breck.xy";
+
+$plate_file = "$dir0/plates/plate_boundaries/bird_boundaries";
+$plate_labels = "$dir0/socal_2005/socal_plate_labs.xyz";
+
+# continental shelf
+$ishelf = 0;
+$shelf_file = "/home/carltape/gmt/coastlines/socal_shelf";
+$Wshelf = "-W0.75p,0/0/0,--";
+
+# CMT sources
+$dir_sources = "/home/carltape/results/SOURCES/socal_16";
+$outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
+$inner_boundary = "${dir_sources}/SPECFEM_inner_points.dat";
+
+# cities
+$cities = "$dir0/cities/socal_tomo_cities";
+
+# socal map
+#$tick1 = 100; $tick2 = 100; # no ticks for map
+$tick1 = 1; $tick2 = 1;
+$tlen_map = "3p";
+#$Rmap = "-R-121.2/-114.8/32.3/36.7";
+$Rmap = "-R-122/-114/32/37";
+
+$topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
+$fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
+$coast_res = "-Df -A0/0/4";
+
+$B = "$B0".$Bopts[0];
+
+$iportrait = 0;
+if($iportrait == 1){$orient = "-P"; $rotangle = 0}
+else {$orient = " "; $rotangle = 90}
+
+# plotting specifications
+$fsize0 = "18";
+$fsize1 = "16";
+$fsize2 = "14";
+$fsize3 = "12";
+$fontno = "1"; # 1 or 4
+$tlen = "6p";
+$fpen = "1p";
+$tpen = "1p";
+
+# plotting specificiations
+# NOTE: pen attribute (-W) does not work for psvelo
+#$coast_info = "-A5000 -Dl -W1.0p -S220/255/255 -C220/255/255 -G200"; # A : smallest feature plotted, in km^2; D : resolution
+$coast_res = "-Df -A0/0/4";
+#$coast_info = "${coast_res} -W0.75p -Na/0.75p";
+$coast_info = "${coast_res} -W1p -Na/1.0p -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
+$coast_info2 = "${coast_res} -S220/255/255 -C220/255/255 -G200"; # -I1/1.0p for rivers
+$fault_info = "-M -W0.5p,255/0/0";
+$fault_infoK = "-M -W0.5p,0/0/0";
+
+$sky_color = "100/100/100";
+$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
+
+#$colinfo = "-W1p,0/0/0 -G255/255/255";
+$colinfo_map = "-W1.5p,0/0/0 -G0/255/255";
+$colinfo_xc = "-W2p,0/0/0";
+$srcinfo_xc = "-Sa30p $colinfo_xc -N";
+$recinfo_xc = "-Si30p $colinfo_map -N";
+
+$textinfo = "-G255 -S1.5p";
+$textinfo2 = "-G0/0/0 -S2p,255/255/255";
+$textinfo3 = "-G0/0/0 -S2p,255/255/0";
+
+# BOOLEAN: WHICH FIGURES TO PLOT
+$ixv = 0;
+$ipdf = 0;
+$ieps = 1;
+$ijpg = 1;
+
+# plot title
+$J_title = "-JM1";
+$R_title = "-R0/1/0/1";
+$fsize_title = $fsize0;
+
+#---------------------------------------------------------
+
+# load values and check that the files exist
+($nump,undef,undef) = split(" ",`ls -1 ${dirdat}/*path | wc`);
+if($nump == 0) {die("No ray paths available")}
+if($pmax > $nump) {$pmax = $nump;}
+if($pmin < 1) {$pmin = 1;}
+print "\nNumber of cuts is $nump\n";
+
+# load fault positions along each profile
+ $ffile = "$dirdat/ALL_${stirun}_rays_fault_positions_mod";
+ if (not -f $ffile) {die("Check if ffile $ffile exist or not\n");}
+ open(IN,$ffile); @faults = <IN>; close(IN);
+
+$cshfile = "plot_vert_models_one.csh";
+open(CSH,">$cshfile");
+print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter MEASURE_UNIT inch BASEMAP_TYPE plain PLOT_DEGREE_FORMAT D TICK_LENGTH $tlen LABEL_FONT_SIZE $fsize2 ANOT_FONT_SIZE $fsize2 HEADER_FONT $fontno ANOT_FONT $fontno LABEL_FONT $fontno HEADER_FONT_SIZE $fsize1 FRAME_PEN $fpen TICK_PEN $tpen \n";
+
+ at iBs = (5,8,4); # what sides of each cross section to show tick marks
+
+# modify the height of each cross-section based on the length
+for ($p = $pmin; $p <= $pmax; $p ++ ) {
+ $stip = sprintf("%3.3i",$p);
+
+ # datafile -- COLUMNS: lon, lat, m00, m16
+ $dfile = "${dirdat}/vert_${stirun}_xc_vsvbvp_${smodeltag}_${stip}.dat";
+ if (not -f $dfile) {die("Check if dfile $dfile exist or not\n");}
+ print "Data file is $dfile\n";
+
+ # source and receiver points -- NOT necessarily the endpoints of the ray path
+ # COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
+ if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
+ open(IN,$rfile); @line = <IN>; close(IN);
+ ($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
+ #print "-- $slon -- $slat -- $sdist -- $sdep --\n";
+ #print "$rlon -- $rlat -- $rdist -- $rdep -- $raz -- $iflip--\n";
+
+ # positions of six faults along each profile
+ # (1) SAF, (2) GF, (3) SGF, (4) MCF, (5) SYF, (6) CRF, (7) EF, (8) SN, (9) KC
+ (undef,$eid,$stanet,$azi1,$azi2,$saf,$gf,$sgf,$mcf,$syf,$crf,$ef,$snf,$kcf) = split(" ",$faults[$p-1]);
+ ($sta,$net) = split("\\.",$stanet); # only for irun = 1
+ $sazi1 = sprintf("%3.3i",$azi1);
+ $sazi2 = sprintf("%3.3i",$azi2);
+ print "-- $stanet -- $sta -- $net -- $eid -- $sazi1 -- $sazi2 --\n";
+
+ # axes limits for cross sections
+ ($dmin,$dmax,$zmin,$zmax) = split(" ",`minmax -C -I2 $dfile`);
+ $dran = $dmax - $dmin; $zran = $zmax - $zmin;
+ $heightxc = $heightxc0;
+ #if ($dran > 450) {
+ # $heightxc = 1/($dran/700 + 0.03); # empirical relationship
+ #}
+
+ # compute the height, including vertical exaggeration
+ $wid = $heightxc*$dran/$zran / $vexag;
+ $dfile = $dfile;
+ $R = "-R$dmin/$dmax/$zmin0/$zmax0";
+ $J = "-JX$wid/$heightxc";
+ if ($iflip==1) {$J = "-JX-$wid/$heightxc";}
+
+ $widmap = $heightxc*(1.9/1.5); # width of map -- tied to height of xc
+ $Jmap = "-JM${widmap}";
+ $B0map = "-Ba${tick1}f${tick2}d:.\" \":";
+ #$Bmap = "$B0map".$Bopts[7];
+ $Bmap = "$B0map".$Bopts[15];
+
+ #---------------------------------------------------------
+ # USER INPUT
+
+ # parameters controlling the placement of subplots and labels
+ $x0_xc = 0.75; $y0_xc = 3; # lower left of m00 cross section
+ #$x0_map = 0.6;
+ #$y0_map = 5;
+ #$xgap_map_xc = 0.8;
+ #$x0_xc = $x0_map + $widmap + $xgap_map_xc;
+ #$y0_xc = $y0_map;
+ $ygap_lab = 0.1; $xgap_lab = 0.6*$ygap_lab; # labels on the xc
+ $xgap_cbar = 0.2; # gap between plot and cbar
+ #$omap = "-Xa${x0_map} -Ya${y0_map}";
+ #$oxc = "-Xa${x0_xc} -Ya${y0_xc}";
+
+ # factors (times $heightxc) controlling y position of colorbar and labels
+ $fygap_bot = 0.1;
+ $f_cbar = 0.6;
+ $fygap_top = 0.20;
+ $f_lab = $fygap_bot + $f_cbar + $fygap_top;
+ $Dthick = 0.20; # thickness of colorbar
+ $Dlen = $heightxc*$f_cbar; # length of colorbar
+ $Dscale0 = "-D0/0/${Dlen}/${Dthick}";
+ #@Dscales = ($Dscale0,$Dscale0,"${Dscale0}h");
+ @Dscales = ($Dscale0,$Dscale0,$Dscale0);
+
+ # factors controlling position of the vertical exaggeration label
+ $xv = $x0_xc;
+ $yv = $y0_xc + $heightxc + 0.5;
+ $ov = "-Xa${xv} -Ya${yv}";
+
+ #---------------------------------------------------------
+
+ # scale for plots
+ $tick1x = 50; $tick2x = 10;
+ $tick1z = 10; $tick2z = 5;
+ if ($dran > 550) {
+ $tick1x = 100; $tick1z = 20;
+ }
+ @Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+ $B0 = sprintf("-Ba%3.3ff%3.3fd:\" \":/a%3.3ff%3.3fd:\" \"::.\" \":WesN",$tick1x,$tick2x,$tick1z,$tick2z);
+ #$B0 = sprintf("-Ba%3.3ff%3.3fg%3.3f:\" \":/a%3.3ff%3.3fg%3.3f:\" \"::.\" \":WesN",$tick1x,$tick2x,$tick1x,$tick1z,$tick2z,$tick1z); # with gridlines
+
+
+ #==================================================
+
+ #-------------------
+ # absolute origins -- these depend on the WIDTH of the cross section
+
+ # cross sections
+ $xgap_xc = 0.5; # hspace between cross sections
+ #if($wid < 3.5) {$xgap_xc = 1.0;}
+ $x1 = $x0_xc; $y1 = $y0_xc; $oxc1 = "-Xa${x1} -Ya${y1}";
+ $x2 = $x0_xc; $y2 = $y0_xc; $oxc2 = "-Xa${x2} -Ya${y2}";
+ $x3 = $x0_xc; $y3 = $y0_xc; $oxc3 = "-Xa${x3} -Ya${y3}";
+ @oxcs = ($oxc1,$oxc2,$oxc3);
+
+ # cross-section label
+ $x1_label = $x1 + $xgap_lab; $y1_label = $y1 + $ygap_lab; $oxclab1 = "-Xa${x1_label} -Ya${y1_label}";
+ $x2_label = $x2 + $xgap_lab; $y2_label = $y2 + $ygap_lab; $oxclab2 = "-Xa${x2_label} -Ya${y2_label}";
+ $x3_label = $x3 + $xgap_lab; $y3_label = $y3 + $ygap_lab; $oxclab3 = "-Xa${x3_label} -Ya${y3_label}";
+ @oxclabs = ($oxclab1,$oxclab2,$oxclab3);
+
+ # colorbar
+ $fac = $heightxc*$fygap_bot + $Dlen/2;
+ $xcb1 = $x1 + $wid + $xgap_cbar; $ycb1 = $y1 + $fac; $ocb1 = "-Xa${xcb1} -Ya${ycb1}";
+ $xcb2 = $x2 + $wid + $xgap_cbar; $ycb2 = $y2 + $fac; $ocb2 = "-Xa${xcb2} -Ya${ycb2}";
+ $xcb3 = $x3 - $xgap_cbar - $Dthick; $ycb3 = $y3 + $fac; $ocb3 = "-Xa${xcb3} -Ya${ycb3}";
+ #$xcb3 = $x0_map + $widmap + $Dlen/2 + 0.4; $ycb3 = $y3 + $heightxc + 1.0; $ocb3 = "-Xa${xcb3} -Ya${ycb3}";
+ #$xcb3 = $x3 + $wid - $Dlen/2; $ycb3 = $y3 + $heightxc + 1.0; $ocb3 = "-Xa${xcb3} -Ya${ycb3}";
+ @ocbs = ($ocb1,$ocb2,$ocb3);
+
+ # label for colorbar
+ $xcblab1 = $xcb1 + $xgap_cbar; $ycblab1 = $y1 + $heightxc*$f_lab; $ocblab1 = "-Xa${xcblab1} -Ya${ycblab1}";
+ $xcblab2 = $xcb2 + $xgap_cbar; $ycblab2 = $y2 + $heightxc*$f_lab; $ocblab2 = "-Xa${xcblab2} -Ya${ycblab2}";
+ $ocblab3 = "-Xa0 -Ya0";
+ @ocbarlabs = ($ocblab1,$ocblab2,$ocblab3);
+
+ #-------------------
+
+ #$dtitle = sprintf("Depth = %.1f km",$zdep/1000);
+
+ $cpert1 = $vpert;
+
+ # colorpoint file for m00 and m16 -- VARIES with depth -- in KILOMETERS
+ $cptfile1 = "color1.cpt";
+ $dc = ($vmax-$vmin)/${scale_color};
+ $T = sprintf("-T%3.3e/%3.3e/%3.3e",$vmin,$vmax,$dc);
+ print CSH "makecpt -C$colorbar $T -D -Z > color.cpt\n"; # -Z for continuous (not for pscontour)
+ print CSH "sed 's/^F.*/F 200 200 200/' color.cpt > $cptfile1\n";
+
+ # colorpoint file for perturbation ln(m16/m00) -- FIXED with depth
+ $cptfile2 = "color2.cpt";
+ $cmin = -$cpert2; $cmax = $cpert2;
+ $dc = ($cmax-$cmin)/${scale_color};
+ $T = sprintf("-T%3.3e/%3.3e/%3.3e",$cmin*1.01,$cmax*1.01,$dc);
+ print CSH "makecpt -C$colorbar $T -D > $cptfile2\n";
+
+ # colorbars for cross sections
+ $Bscale1 = sprintf("-B%2.2ef0.5:\" \": -E10p",$vtick);
+ #$Bscale2 = sprintf("-Ba%2.2ef0.05:\"$titles[2]\": -E10p",$cpert2);
+ $Bscale2 = sprintf("-Ba%2.2ef0.05:\" \": -E10p -A",$cpert2); # -A for ticks on other side
+
+ #------------------------------------
+
+ # plot ray path
+ $pfile = "${dirdat}/vert_${stirun}_xc_${stip}_ray_path";
+ if (not -f $pfile) {
+ die("Check if pfile $pfile exist or not\n");
+ }
+
+ # plot source and receiver lon-lat
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
+ if (not -f $rfile) {die("Check if rfile $rfile exist or not\n");}
+
+ #-----------------------------------------------
+
+ $k = 2;
+
+ #$shift = $shifts[$k-1];
+ $iB = $iBs[$k-1];
+ $B = "$B0".$Bopts[$iB];
+ $title = $titles[$k-1];
+
+ #------------------
+
+ # vertical lines
+ $zminf0 = -5;
+ $finfo = "-W1.5p,0/0/0";
+
+ # details for specific cross-sections
+ #if ($irun==1 && $k == 1) {
+ (undef,undef,undef,undef,undef,$saf,$gf,$sgf,$mcf,$syf) = split(" ",$faults[$p-1]);
+ $talign = "CB";
+ $flsize = 9;
+ $textinfo = "-N -C3p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+
+#==================================================
+# figure with map with m16 beneath
+
+if ($imfinal == 1) {
+
+ @iBs = (5,8,4);
+
+ #for ($p = $pmin; $p <= $pmax; $p ++ ) {
+
+ $stip = sprintf("%3.3i",$p);
+
+ $fname = "vert_${stirun}_xc_${modlab}_${smodeltag}_${stip}_single";
+ if($irun==1) {$fname = "vert_${stirun}_xc_${sta}_${net}_${sazi1}_${eid}_${modlab}_${smodeltag}_${stip}";}
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg"; $epsfile = "$fname.eps";
+
+ # datafile -- COLUMNS: lon, lat, m00, m16
+ $dfile = "${dirdat}/vert_${stirun}_xc_vsvbvp_${smodeltag}_${stip}.dat";
+ if (not -f $dfile) {
+ die("Check if dfile $dfile exist or not\n");
+ }
+ print "Data file is $dfile\n";
+
+ # source and receiver points -- NOT necessarily the endpoints of the ray path
+ # COLUMNS: slons(ii),slats(ii),0,rlons(ii),rlats(ii),deg2km(dq),azq,iflip
+ $rfile = "${dirdat}/vert_${stirun}_xc_${stip}_A_B";
+ if (not -f $rfile) {
+ die("Check if rfile $rfile exist or not\n");
+ }
+ open(IN,$rfile); @line = <IN>; close(IN);
+ ($slon,$slat,$sdist,$sdep,$rlon,$rlat,$rdist,$rdep,$raz,$iflip) = split(" ",$line[0]);
+
+ # axes limits for cross sections
+ ($dmin,$dmax,$zmin,$zmax) = split(" ",`minmax -C -I2 $dfile`);
+ $dran = $dmax - $dmin; $zran = $zmax - $zmin;
+
+ # compute the height, including vertical exaggeration
+ #$heightxc = $zran*$wid/$dran * $vexag;
+ $wid = $heightxc*$dran/$zran / $vexag;
+ $R = "-R$dmin/$dmax/$zmin0/$zmax0";
+ $J = "-JX$wid/$heightxc";
+ if ($iflip==1) {
+ $J = "-JX-$wid/$heightxc";
+ }
+
+ #-------------------
+ # absolute origins -- these depend on the WIDTH of the cross section
+
+ # overall label
+ $x0_label = $x0_xc;
+ $y0_label = $y0_xc + $heightxc + 1;
+ $olab = "-Xa${x0_label} -Ya${y0_label}";
+ $x0_label2 = $x0_label;
+ $y0_label2 = $y0_xc + 0.8*$heightxc + 1;
+ $olab2 = "-Xa${x0_label2} -Ya${y0_label2}";
+
+ #-----------------------------------------------
+
+ #$shift = $shifts[$k-1];
+ $iB = $iBs[$k-1];
+ $B = "$B0".$Bopts[$iB];
+ $title = $titles[$k-1];
+
+ #------------------
+ # get absolute origins for various pieces
+
+ $oxc = $oxcs[$k-1];
+ $ocb = $ocbs[$k-1];
+ $ocbarlab = $ocbarlabs[$k-1];
+ $Dscale = $Dscales[$k-1];
+
+ #------------------
+
+ print CSH "psbasemap $J $R $B -K -V $oxc > $psfile\n";
+ print CSH "psbasemap $J $R $B -G${sky_color} -K -O -V $oxc >> $psfile\n";
+
+ if ($icolor==1) {
+ print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$4}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile1 $J -Q -K -O -V $oxc >> $psfile\n";
+
+ print CSH "grdcontour $grdfile -Cconfile1 -A- -W1.5p -Q100 $J -K -O -V $oxc >> $psfile\n";
+ }
+
+ print CSH "psscale -C${cptfile1} $Dscale $Bscale1 -K -O -V $ocb >> $psfile\n";
+ $slab1 = "km/s";
+ print CSH "pstext -N $R_title $J_title -K -O -V $ocbarlab >>$psfile<<EOF\n0 0 $fsize2 0 $fontno CM $slab1\nEOF\n";
+
+ # details for specific cross-sections
+ #if ($irun==1) {
+ # positions of five faults along each profile
+ # (1) SAF, (2) GF, (3) SGF, (4) MCF, (5) SYF, (6) CRF
+ #(undef,undef,undef,$saf,$gf,$sgf,$mcf,$syf,$crf) = split(" ",$faults[$p-1]);
+ #print "$p -- SAF $saf -- GF -- $gf -- SGF $sgf -- MCF $mcf -- SYF $syf\n"; die("TESTING");
+
+ # vertical lines
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$saf $zmax0\n$saf $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$gf $zmax0\n$gf $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$sgf $zmax0\n$sgf $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$mcf $zmax0\n$mcf $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$syf $zmax0\n$syf $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$crf $zmax0\n$crf $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$ef $zmax0\n$ef $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$snf $zmax0\n$snf $zminf0\nEOF\n";
+ print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$kcf $zmax0\n$kcf $zminf0\nEOF\n";
+
+ #print CSH "psxy $J $R -W1p,0/0/0,-- -K -O -V $oxc >>$psfile<<EOF\n$dmin -10\n$dmax -10\nEOF\n";
+
+ # plot source and receiver
+ if ($irun == 1) {
+ print CSH "awk '{print \$3,-\$4}' $rfile | psxy $J $R $srcinfo_xc -K -O -V $oxc >> $psfile\n";
+ print CSH "awk '{print \$7,$zmax0}' $rfile | psxy $J $R $recinfo_xc -K -O -V $oxc >> $psfile\n";
+ }
+
+ # details for specific cross-sections
+ #if ($irun==1) {
+ #(undef,undef,undef,$saf,$gf,$sgf,$mcf,$syf) = split(" ",$faults[$p-1]);
+ $talign = "CB";
+ $flsize = $fsize1;
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ if($saf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$saf $zmax0 $flsize 0 $fontno $talign SA\nEOF\n";}
+ if($gf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$gf $zmax0 $flsize 0 $fontno $talign G\nEOF\n";}
+ if($sgf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$sgf $zmax0 $flsize 0 $fontno $talign SG\nEOF\n";}
+ if($mcf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$mcf $zmax0 $flsize 0 $fontno $talign MC\nEOF\n";}
+ if($syf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$syf $zmax0 $flsize 0 $fontno $talign SY\nEOF\n";}
+ if($crf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$crf $zmax0 $flsize 0 $fontno $talign CR\nEOF\n";}
+ if($ef != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$ef $zmax0 $flsize 0 $fontno $talign E\nEOF\n";}
+ if($snf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$snf $zmax0 $flsize 0 $fontno $talign SN\nEOF\n";}
+ if($kcf != 9999) {print CSH "pstext $J $R $textinfo -K -O -V $oxc >>$psfile<<EOF\n$kcf $zmax0 $flsize 0 $fontno $talign KC\nEOF\n";}
+ #}
+
+ # plot overall label (for publication)
+ if ($iletter==1) {
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 18 0 $fontno TL $letter\nEOF\n";
+ } else {
+ if ($irun==1) {
+ $textinfo = "-N";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 $fsize1 0 $fontno TC $eid\nEOF\n";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab2 >>$psfile<<EOF\n0 0 $fsize1 0 $fontno TC $stanet\nEOF\n";
+ }
+ }
+
+ # plot label on each cross section
+ for ($k = 2; $k <= 2; $k ++ ) {
+ $title = $titles[$k-1];
+ $talign = "LB";
+ #$otitle = "-Xa${x3} -Ya${y3}";
+ $otitle = $oxclabs[$k-1];
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,0/0/0,solid -G0/0/0";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $otitle >>$psfile<<EOF\n0 0 $fsize0 0 $fontno $talign $title\nEOF\n";
+ }
+
+ #------------------------------------
+
+ print CSH "pstext -N $R_title $J_title -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n"; # FINISH
+ #if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
+ if ($ixv==1) {print CSH "gv $psfile &\n";}
+ #if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+ if ($ieps==1) {print CSH "ps2eps -f -R + -l $psfile\n";}
+ #if ($ijpg==1) {print CSH "convert -rotate $rotangle $psfile $jpgfile\n";}
+ if ($ijpg==1) {print CSH "convert $epsfile $jpgfile\n";}
+
+ } # imfinal == 1
+} # loop over p
+
+#==================================================
+
+#------------------------------------
+
+close (CSH);
+system("csh -f $cshfile");
+
+# print "convert $psfile $jpgfile \n";
+# system("convert $psfile -rotate $rotangle $jpgfile");
+# if ($ipdf==1) {system("ps2pdf $psfile")}
+# if ($ixv==1) {system("gv $psfile &");}
+system("gv $psfile &");
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/scripts/plot_vert_models_one.pl
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list