[cig-commits] r14507 - in seismo/3D/ADJOINT_TOMO/iterate_adj: . model_plot model_plot/gmt model_plot/matlab
carltape at geodynamics.org
carltape at geodynamics.org
Fri Mar 27 18:07:41 PDT 2009
Author: carltape
Date: 2009-03-27 18:07:39 -0700 (Fri, 27 Mar 2009)
New Revision: 14507
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/README
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/SAVE/
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/PPRmpa.sty
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/plot_vert_models_basemap.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_xc_seis_all.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/ps2eps.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/xc_and_seis.tex
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m
Log:
Matlab and GMT scripts for plotting cross sections that have been generated from binary-stored, large meshes from SPECFEM3D simulations.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/README 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,17 @@
+Carl Tape
+27-March-2009
+
+A preliminary set of scripts for generating xyz files for cross-sections to be extracted from volumetric .bin files (for example, a Vp model or a sensitivity kernel). These have been used for the southern California crustal tomography study. Significant work is needed to generalize the codes for any user.
+
+Basic order of operations:
+1. matlab -- Generate target cross section (horizontal or vertical)
+2. CLUSTER -- extract model value of the nearest GLL point to each target point
+3. gmt -- plot cross sections
+
+Choices to be made include:
+ (1) Vertical or horizontral cross section
+ (2) Sampling density in the horizontal dimension
+ (3) Sampling density in the vertical dimension (for vertical xc only)
+ (4) Area of the cross section
+
+------------------
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/PPRmpa.sty
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/PPRmpa.sty (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/PPRmpa.sty 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,58 @@
+%==============================================================================
+% Prosper -- (PPRmpa.sty) Style file
+% A LaTeX class for creating slides
+% Author: Georg Drenkhahn (georg at mpa-garching.mpg.de)
+% 2001-10-18
+%==============================================================================
+\NeedsTeXFormat{LaTeX2e}[1995/12/01]
+\ProvidesPackage{PPRmpa}[2001/10/18]
+\typeout{`MPA' style for prosper}
+
+%% already included by prosper:
+% seminar, graphicx, hyperref, ifthen
+
+%% everybody needs color :-)
+%% not required for MPA style file itself
+\RequirePackage{color}
+
+%% redefine some standard fonts: Seminar-Helvetica
+%% \RequirePackage{semhelv}
+
+%% Type 1 CM fonts
+\RequirePackage{ae,aecompl}
+
+\newrgbcolor{mpagreena}{0.188 0.384 0.371}
+\newrgbcolor{mpagreenb}{0.235 0.486 0.447}
+\newrgbcolor{mpagreenc}{0.263 0.537 0.506}
+
+\FontTitle{%
+ \sffamily\bfseries\fontsize{20.74pt}{18pt}\selectfont\mpagreenb}{%
+ \sffamily\bfseries\fontsize{20.74pt}{18pt}\selectfont\mpagreenb}
+
+\FontText{%
+ \sffamily\mdseries\fontsize{12.4}{12}\selectfont\black}{%
+ \sffamily\mdseries\fontsize{12.4}{12}\selectfont\black}
+
+% itemisation bullets
+\myitem{1}{{\scriptsize\mpagreena\raisebox{2pt}{\ensuremath{\bullet}}}}
+\myitem{2}{{\scriptsize\mpagreenb\raisebox{2pt}{\ensuremath{\bullet}}}}
+\myitem{3}{{\scriptsize\mpagreenc\raisebox{2pt}{\ensuremath{\bullet}}}}
+
+% Positionning of the title of a slide.
+% A4 format: 29.7x21cm 11.7x8.3in
+\newcommand{\slidetitle}[1]{%
+ \rput[t](5.45,4.5){\fontTitle{#1}}%
+}
+
+\newcommand{\theMPAStyle}[1]{%
+ \PutLogo%
+ {#1}}
+
+\LogoPosition{-1.1,-1.2}
+
+%% coordinate zero point is on the left edge at middle height
+%% the top center of the box is set at this point.
+\NewSlideStyle{t}{5.45,3.5}{theMPAStyle}
+%% does it have any effect?
+\PDFCroppingBox{10 40 594 820}
+\endinput
Added: 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 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,385 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_horz_coverage.pl
+# Carl Tape
+# 07-Dec-2008
+#
+# horz_xc_vs_coverage_003.dat
+# This script reads in a horizontal cross-section data file with five columns:
+# lon, lat, log10(Ksum), log10(Ksum)+mask
+# 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 21/21 0/4 0/1/3
+# plot_horz_coverage.pl 1 21/21 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");}
+($imask,$pinds,$cvals,$parms) = @ARGV;
+
+$icolor = 1;
+($pmin,$pmax) = split(/\//,$pinds);
+($cmean,$cpert) = split(/\//,$cvals);
+($itoplab,$isidelab,$iB) = split(/\//,$parms);
+
+# KEY COMMAND: HOW DID YOU WEIGHT THE SUMMED KERNELS
+$smodeltag = "coverage_sum_abs_nwin";
+$smodeltag = "coverage_sum_abs";
+
+# detailed flags for publication figures
+#$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
+$iinsetlab = 0; # label the model inan inset
+$iletter = 0; # label for publication figures
+$letter = "A";
+
+$cgray = "-G220/220/220";
+
+if($ivs==1) {$modlab = "vs"; $mtit = "K_Vs";} else {$modlab = "vb"; $mtit = "K_Vb";}
+
+#---------------------------------------------------------
+# USER INPUT
+
+# directory containing the data files
+$dirdat = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/horz_full";
+
+# KEY: file showing the cuts and the plotting range for the values
+$fcuts = "${dirdat}/horz_xc_elevation_cuts_mod";
+if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
+
+# load values and check that the files exist
+open(IN,$fcuts); @lines = <IN>; $nump = @lines; close(IN);
+for ($p = 1; $p <= $nump; $p ++ ) {
+ $stip = sprintf("%3.3i",$p);
+ (undef,$zcut,$vs_norm,$vb_norm,$vpert) = split(" ",$lines[$p-1]);
+ $vperts[$p-1] = $vpert;
+ $vsnorms[$p-1] = $vs_norm;
+ $vbnorms[$p-1] = $vb_norm;
+ $zcuts[$p-1] = -$zcut;
+ $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}.dat";
+ #if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
+}
+print "\nNumber of cuts is $nump\n";
+if($pmax > $nump) {$pmax = $nump;}
+if($pmin < 1) {$pmin = 1;}
+
+# subtitles
+ at titles = ("$mtit coverage");
+ at labs = ("$mtit");
+
+# parameters for color scales
+$cpert1 = 0.2;
+#$cpert2 = 0.2;
+$scale_color = 65.0;
+$colorbar = "seis";
+
+# resolution of color plots
+#$interp = "-I2m/2m -S4m"; # key information
+$interp = "-I2m";
+$grdfile = "temp.grd";
+
+# position of titles and labels
+$xtx1 = 0.5; $ytx1 = 1.00;
+$xtx2 = -0.15; $ytx2 = 0.4;
+$xtx3 = 0.5; $ytx3 = -0.1;
+$xtx4 = $xtx3; $ytx4 = -0.3;
+$xtx5 = -0.15; $ytx5 = 0.85; # A, B, C, etc
+$xtx6 = 0.0; $ytx6 = 0.05; # m00, m16, ln(m16/m00)
+
+#print "$dtitle\n"; die("TESTING");
+
+#---------------------------------------------------------
+# DATA FILES
+
+$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_tomo = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO";
+$dir_sources = "/net/sierra/raid1/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";
+
+# check that files exist
+#if (not -f $xxx) { die("Check if $xxx exist or not\n") }
+
+# bounds and dimensions
+# 0 = LA
+# 1 = southern California
+$iregion = 1;
+
+if ($iregion==1) {
+ # southern California
+ $wid = 2.8; # 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; $tick1 = 1; $tick2 = 0.5;
+ $xmin = -122; $xmax = -114; $ymin = 32; $ymax = 37;
+ $xtick1 = 2; $xtick2 = 0.5; $ytick1 = 1; $ytick2 = 0.5;
+
+ $cmax = 5000; $cmin = -$cmax;
+ $y_title = 0.98; $y_title2 = 0.92; $x_title = 0.5; $x_title2 = 0.5;
+ $iportrait = 0;
+ $stag = "Southern California";
+ $topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
+ $fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
+ $name = "plot_horz_coverage";
+ $origin = "-X1i -Y5i";
+
+ # inset map
+ $iinset = 1;
+ $origin_inset = "-Xa7.5 -Ya4.5";
+ $Jinset = "-JM1.5";
+ $Rinset = "-R-132/-110/25/50";
+ $Binset = "-B100wesn";
+ $coast_res = "-Df -A0/0/4";
+}
+
+$dX = $wid + 0.5;
+$hwid = $wid/2;
+#$R = "-Rg";
+$R = "-R$xmin/$xmax/$ymin/$ymax";
+#$J = "-JK180/${wid}i";
+$J = "-JM${wid}";
+#$J = "-Jm0.70";
+
+# scale for plots
+ at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+
+$B0 = "-Ba${xtick1}f${xtick2}d:\" \":/a${ytick1}f${ytick2}d:\" \":";
+#$B0 = "-Ba${xtick1}f${xtick2}d:.\" \":";
+#$B0 = "-Ba${tick1}f${tick1}g${tick1}:.\" \":"; # with gridlines
+
+$B = "$B0".$Bopts[0];
+
+if($iportrait == 1){$orient = "-P"; $rotangle = 0}
+else {$orient = " "; $rotangle = 90}
+
+# color bar
+$Dlen = $wid*0.4;
+$Dx = $Dlen*0.5;
+$Dy = -0.2;
+$Dscale = "-D$Dx/$Dy/$Dlen/0.15h";
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "12";
+$fsize2 = "10";
+$fsize3 = "6";
+$fontno = "1"; # 1 or 4
+$tick = "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_info = "${coast_res} -W0.75p -Na/0.75p";
+$fault_info = "-M -W0.75p,255/0/0";
+$fault_infoK = "-M -W0.75p,0/0/0";
+
+$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
+
+$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 = 1;
+#$ixv = 1; $ipdf = 0;
+
+# plot title
+$J_title = "-JM${wid}";
+$R_title = "-R0/1/0/1";
+$fsize_title = $fsize0;
+
+#==================================================
+
+$cshfile = "plot_horz_coverage.csh";
+open(CSH,">$cshfile");
+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";
+
+ at shifts = ($origin,"-X$dX","-X$dX");
+#if($ilonlab==1) {@iBs = (2,14,2);} else {@iBs = (15,9,15);}
+
+for ($p = $pmin; $p <= $pmax; $p ++ ) {
+
+ $stip = sprintf("%3.3i",$p);
+ $fname = "horz_xc_${modlab}_${smodeltag}_${stip}_${imask}";
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+ $vsnorm = $vsnorms[$p-1]; # reference velocity for m00 and m16
+ $vbnorm = $vbnorms[$p-1]; # reference velocity for m00 and m16
+ $vpert = $vperts[$p-1]; # perturbation from reference velocity to plot
+ $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";
+ 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;
+ #print "-- $cnorm --\n"; die("TESTING:");
+
+ $dtitle = sprintf("Depth = %.1f km",$zdep/1000);
+
+ # datafile -- COLUMNS: lon, lat, m00, m16
+ $dfile = $dfiles[$p-1];
+ if (not -f $dfile) {die("Check if dfile $dfile exist or not\n");}
+ print "Data file is $dfile\n";
+
+ # colorpoint file for m00 and m16 -- VARIES with depth -- in ln(m00/cnorm)
+ $cptfile = "color.cpt";
+# if($cmean >= 0) {
+# $cmin = $cmean * exp(-$cpert);
+# $cmax = $cmean * exp($cpert);
+# } else {
+# $cmin = $cmean * exp($cpert);
+# $cmax = $cmean * exp(-$cpert);
+# }
+ $cmin = -$cpert;
+ $cmax = $cpert;
+ $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 > $cptfile\n";
+ print CSH "makecpt -Chot $T -D > $cptfile\n";
+ #print "\n colorbar is -- $T -- $cmin, $cmean, $cmax \n"; die("TESTING");
+
+ $ctick = $cmax;
+ $Bscale = sprintf("-Ba%2.2ef%2.2e:\" \": -E10p",$ctick,$ctick/2);
+
+ for ($k = 1; $k <= 1; $k ++ ) {
+
+ $shift = $shifts[$k-1];
+ #$iB = $iBs[$k-1];
+ $B = "$B0".$Bopts[$iB];
+ $title = $titles[$k-1];
+ $kcol = 3 + $imask;
+
+ if ($k==1) {
+ print CSH "psbasemap $J $R $B -K -V $orient $origin > $psfile\n"; # START
+ } else {
+ 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 ($icolor==1) {
+ #print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "awk '{print \$1,\$2,log(\$${kcol}/$cnorm)}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ 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);
+ 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";
+ if($ishelf==1) {print CSH "awk '{print \$2,\$1}' ${shelf_file} |psxy $J $R $Wshelf -K -O -P -V >> $psfile\n";}
+
+ #--------------------------------
+
+ #print CSH "psxy ${plate_file} $J $R $plate_info -K -V -O >> $psfile \n";
+ print CSH "psxy ${fault_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${kcf_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${breck_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${fault_file2} $J $R $fault_infoK -K -V -O >> $psfile \n";
+
+ if ($imask==1) {
+ # boundaries of simulation
+ print CSH "psxy ${outer_boundary} $J $R -W2p,0/0/0 -K -O -V >>$psfile\n";
+ #print CSH "psxy ${inner_boundary} $J $R -W1.5p,0/0/0,-- -K -O -V >>$psfile\n";
+
+ $ibox = 0;
+ if ($ibox==1) {
+ $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";
+ print CSH "psxy ${lin_boundary} $J $R $boxinfo -K -O -V >>$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";
+
+ #----------------------------
+
+ # 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 && $k==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 && $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$xtx5 $ytx5 18 0 $fontno TL $letter\nEOF\n";
+ }
+
+ } # 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 "ghostview $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");
+
+# print "convert $psfile $jpgfile \n";
+# system("convert $psfile -rotate $rotangle $jpgfile");
+# if ($ipdf==1) {system("ps2pdf $psfile")}
+# if ($ixv==1) {system("ghostview $psfile &");}
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_coverage.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: 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 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,437 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_horz_models.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_horz_models.pl 1/50 m00/m16 1 0.2 # Vs
+# plot_horz_models.pl 1/50 m00/m16 0 0.3 # Vb
+# plot_horz_models.pl 1/6 m00/m01 1 0.2 # Vs
+#
+# SCIENCE PAPER:
+# 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 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 21/21 m00/m16 0 0/0/1/0.3 # Vb
+#
+#==========================================================
+
+if (@ARGV < 4) {die("Usage: plot_horz_models.pl xxx\n");}
+($pinds,$smodels,$ivs,$parms) = @ARGV;
+
+$icolor = 1;
+($pmin,$pmax) = split(/\//,$pinds);
+($smodel1,$smodel2) = split(/\//,$smodels);
+$smodeltag = "${smodel2}_${smodel1}";
+if($ivs==1) {$modlab = "vs"; $mtit = "Vs";} else {$modlab = "vb"; $mtit = "Vb";}
+($ilonlab,$itoplab,$isidelab,$ctick1b_km) = split(/\//,$parms);
+
+# 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 = 1; # label for publication figures
+$letter = "C";
+
+$cgray = "-G220/220/220";
+
+#---------------------------------------------------------
+# USER INPUT
+
+# directory containing the data files
+$dirdat = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/horz_full";
+
+# KEY: file showing the cuts and the plotting range for the values
+$fcuts = "${dirdat}/horz_xc_elevation_cuts_mod";
+if (not -f $fcuts) {die("Check if fcuts $fcuts exist or not\n")}
+
+# load values and check that the files exist
+open(IN,$fcuts); @lines = <IN>; $nump = @lines;
+for ($p = 1; $p <= $nump; $p ++ ) {
+ $stip = sprintf("%3.3i",$p);
+ (undef,$zcut,$vs_norm,$vb_norm,$vpert) = split(" ",$lines[$p-1]);
+ $vperts[$p-1] = $vpert;
+ $vsnorms[$p-1] = $vs_norm;
+ $vbnorms[$p-1] = $vb_norm;
+ $zcuts[$p-1] = -$zcut;
+ $dfiles[$p-1] = "${dirdat}/horz_xc_${modlab}_${smodeltag}_${stip}.dat";
+ #if (not -f $dfile) {die("Check if dfile $dfile exist or not\n")}
+}
+print "\nNumber of cuts is $nump\n";
+if($pmax > $nump) {$pmax = $nump;}
+if($pmin < 1) {$pmin = 1;}
+
+# subtitles
+ at titles = ("$mtit model $smodel1","$mtit model $smodel2","ln ( ${smodel2} / ${smodel1} )");
+ at labs = ("$smodel1","$smodel2","ln(${smodel2}/${smodel1})");
+
+# parameters for color scales
+$cpert1 = 0.2;
+#$cpert2 = 0.2;
+$scale_color = 65.0;
+$colorbar = "seis";
+
+# resolution of color plots
+#$interp = "-I2m/2m -S4m"; # key information
+$interp = "-I2m";
+$grdfile = "temp.grd";
+
+# 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;
+$xtx5 = -0.15; $ytx5 = 0.85; # A, B, C, etc
+$xtx6 = 0.0; $ytx6 = 0.05; # m00, m16, ln(m16/m00)
+
+#print "$dtitle\n"; die("TESTING");
+
+#---------------------------------------------------------
+# DATA FILES
+
+$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_tomo = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO";
+$dir_sources = "/net/sierra/raid1/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";
+
+# check that files exist
+#if (not -f $xxx) { die("Check if $xxx exist or not\n") }
+
+# bounds and dimensions
+# 0 = LA
+# 1 = southern California
+$iregion = 1;
+
+if ($iregion==1) {
+ # southern California
+ $wid = 2.8; # 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;
+ $cmax = 5000; $cmin = -$cmax;
+ $y_title = 0.98; $y_title2 = 0.92; $x_title = 0.5; $x_title2 = 0.5;
+ $iportrait = 0;
+ $stag = "Southern California";
+ $topo_labels = "$dir0/socal_2005/socal_topo_labs.xyz";
+ $fault_labels = "$dir0/socal_2005/socal_fault_labs.xyz";
+ $name = "plot_horz_models";
+ $origin = "-X1i -Y5i";
+
+ # inset map
+ $iinset = 1;
+ $origin_inset = "-Xa7.5 -Ya4.5";
+ $Jinset = "-JM1.5";
+ $Rinset = "-R-132/-110/25/50";
+ $Binset = "-B100wesn";
+ $coast_res = "-Df -A0/0/4";
+}
+
+$dX = $wid + 0.5;
+$hwid = $wid/2;
+#$R = "-Rg";
+$R = "-R$xmin/$xmax/$ymin/$ymax";
+#$J = "-JK180/${wid}i";
+$J = "-JM${wid}";
+#$J = "-Jm0.70";
+
+# scale for plots
+ at Bopts = ("WESN","Wesn","wesN","wEsn","weSn","WesN","wEsN","wESn","WeSn","WEsn","weSN","wESN","WESn","WeSN","WEsN","wesn");
+$B0 = "-Ba${tick1}f${tick2}d:.\" \":";
+#$B0 = "-Ba${tick1}f${tick1}g${tick1}:.\" \":"; # with gridlines
+
+$B = "$B0".$Bopts[0];
+
+if($iportrait == 1){$orient = "-P"; $rotangle = 0}
+else {$orient = " "; $rotangle = 90}
+
+# color bar
+$Dlen = $wid*0.6;
+$Dx = $Dlen*0.5;
+$Dy = -0.5;
+$Dscale = "-D$Dx/$Dy/$Dlen/0.15h";
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "12";
+$fsize2 = "10";
+$fsize3 = "6";
+$fontno = "1"; # 1 or 4
+$tick = "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_info = "${coast_res} -W0.75p -Na/0.75p";
+$fault_info = "-M -W0.75p,255/0/0";
+$fault_infoK = "-M -W0.75p,0/0/0";
+
+$city_info = "-Sc10p -W1.0p/0/0/0 -G255/255/0";
+
+$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 = 1;
+#$ixv = 1; $ipdf = 0;
+
+# plot title
+$J_title = "-JM${wid}";
+$R_title = "-R0/1/0/1";
+$fsize_title = $fsize0;
+
+#==================================================
+
+$cshfile = "plot_horz_models.csh";
+open(CSH,">$cshfile");
+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";
+
+ at shifts = ($origin,"-X$dX","-X$dX");
+if($ilonlab==1) {@iBs = (2,14,2);} else {@iBs = (15,9,15);}
+
+for ($p = $pmin; $p <= $pmax; $p ++ ) {
+
+ $stip = sprintf("%3.3i",$p);
+ $fname = "horz_xc_${modlab}_${smodeltag}_${stip}";
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+ $vsnorm = $vsnorms[$p-1]; # reference velocity for m00 and m16
+ $vbnorm = $vbnorms[$p-1]; # reference velocity for m00 and m16
+ $vpert = $vperts[$p-1]; # perturbation from reference velocity to plot
+ $zdep = $zcuts[$p-1]; # depth of the cross-section (m)
+
+ # make less of a range for Vb compared with Vs -- factor
+ if($ivs == 0) {$vpert = 1.0*$vpert;}
+
+ $dtitle = sprintf("Depth = %.1f km",$zdep/1000);
+
+ # datafile -- COLUMNS: lon, lat, m00, m16
+ $dfile = $dfiles[$p-1];
+ if (not -f $dfile) {die("Check if dfile $dfile exist or not\n");}
+ print "Data file is $dfile\n";
+
+ if($ivs == 1) {$cnorm = $vsnorm;} else {$cnorm = $vbnorm;}
+ $cpert1 = $vpert;
+ $cpert2 = $vpert; # NOTE: over-rule the input value
+
+ # colorpoint file for m00 and m16 -- VARIES with depth -- in ln(m00/cnorm)
+ $cptfile1a = "color1a.cpt";
+ $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
+ $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";
+
+ # colorpoint file for perturbation ln(m16/m00) -- FIXED with depth
+ $cptfile2 = "color2.cpt";
+ $cmin = -$cpert2*1.01; $cmax = $cpert2*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 > $cptfile2\n";
+
+ #
+ $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);
+
+# 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");}
+
+ for ($k = 1; $k <= 3; $k ++ ) {
+ #for ($k = 1; $k <= 1; $k ++ ) {
+
+ $shift = $shifts[$k-1];
+ $iB = $iBs[$k-1];
+ $B = "$B0".$Bopts[$iB];
+ $title = $titles[$k-1];
+
+ if ($k==1) {
+ print CSH "psbasemap $J $R $B -K -V $orient $origin > $psfile\n"; # START
+ } else {
+ 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 ($icolor==1) {
+ if ($k == 1) {
+ #print CSH "awk '{print \$1,\$2,log(\$3/$cnorm)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "awk '{print \$1,\$2,log(\$3/$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} $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";
+
+ } elsif ($k == 2) {
+ #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} $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";
+
+ } else {
+ #print CSH "awk '{print \$1,\$2,log(\$4/\$3)}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ #print CSH "awk '{print \$1,\$2,\$5}' $dfile | nearneighbor -G$grdfile $R $interp\n";
+ print CSH "awk '{print \$1,\$2,\$5}' $dfile | xyz2grd -G$grdfile $R $interp\n";
+ print CSH "grdimage $grdfile -C$cptfile2 $J -K -O -V -Q >> $psfile\n";
+
+ print CSH "psscale -C${cptfile2} $Dscale $Bscale2 -K -O -V >> $psfile\n";
+ $slab1 = "ln($smodel2 / $smodel1)";
+ 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";
+ if($ishelf==1) {print CSH "awk '{print \$2,\$1}' ${shelf_file} |psxy $J $R $Wshelf -K -O -P -V >> $psfile\n";}
+
+ #--------------------------------
+
+ #print CSH "psxy ${plate_file} $J $R $plate_info -K -V -O >> $psfile \n";
+ print CSH "psxy ${fault_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${kcf_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${breck_file} $J $R $fault_infoK -K -V -O >> $psfile \n";
+ #print CSH "psxy ${fault_file2} $J $R $fault_infoK -K -V -O >> $psfile \n";
+
+ if (0==1) {
+ # boundaries of simulation
+ print CSH "psxy ${outer_boundary} $J $R -W2p,0/0/0 -K -O -V >>$psfile\n";
+ print CSH "psxy ${inner_boundary} $J $R -W1.5p,0/0/0,-- -K -O -V >>$psfile\n";
+
+ $ibox = 0;
+ if ($ibox==1) {
+ $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";
+ print CSH "psxy ${lin_boundary} $J $R $boxinfo -K -O -V >>$psfile\n";
+ }
+ }
+
+ #--------------------------------
+
+ $finfo1 = "-W3p,0/0/0";
+ $finfo2 = "-W1.5p,0/255/255";
+
+ # plot a line at -119 longitude
+ if($p==11 && $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) {
+ @irs = (2,5); $nray = @irs ;
+ for ($ik = 1; $ik <= $nray; $ik ++ ) {
+ $ir = $irs[$ik-1];
+ $stir = sprintf("%3.3i",$ir);
+ $pfile = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/vert_01/vert_xc_${stir}_ray_path";
+ if (not -f $pfile) {die("Check if pfile $pfile exist or not\n");}
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo1 -K -O -V >> $psfile\n";
+ print CSH "awk '{print \$1,\$2}' ${pfile} | psxy $J $R $finfo2 -K -O -V >> $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";
+
+ #----------------------------
+
+ # 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 && $k==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 && $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$xtx5 $ytx5 18 0 $fontno TL $letter\nEOF\n";
+ }
+
+ } # 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 "ghostview $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");
+
+# print "convert $psfile $jpgfile \n";
+# system("convert $psfile -rotate $rotangle $jpgfile");
+# if ($ipdf==1) {system("ps2pdf $psfile")}
+# if ($ixv==1) {system("ghostview $psfile &");}
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_horz_models.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: 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 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,816 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_vert_models.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.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 1/10 m00/m16 -40/5 1 0.10 3.0 3 1.5 # Vs -- Garlock-perp profiles
+#
+# 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
+# plot_vert_models.pl 1/8 m00/m16 -40/5 3 0.10 3.0 1 1.5 # Vp
+#
+# plot_vert_models.pl 1/3 m00/m16 -40/5 1 0.15 3.0 2 1.5 # Vs
+#
+#==========================================================
+
+if (@ARGV < 8) {die("Usage: plot_vert_models.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 = 0; # plot a 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 = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/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_tomo = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO";
+$dir_sources = "/net/sierra/raid1/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 = "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.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";
+$srcinfo_map = "-Sa20p $colinfo_map -N";
+$recinfo_map = "-Si20p $colinfo_map -N";
+$srcinfo_mapinset = "-Sa10p -W0.5p,0/0/0 -G0/255/255 -N";
+$recinfo_mapinset = "-Si10p -W0.5p,0/0/0 -G0/255/255 -N";
+$colinfo_xc = "-W2p,0/0/0";
+$srcinfo_xc = "-Sa20p $colinfo_xc -N";
+$recinfo_xc = "-Si20p $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;
+
+# 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_rays_fault_positions";
+ $ffile = "$dirdat/ALL_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.csh";
+open(CSH,">$cshfile");
+print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter 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_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";
+ 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 = 5.8; $y0_xc = 5; # 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.15; # 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
+
+
+ #==================================================
+
+ $fname = "vert_${stirun}_xc_${modlab}_${smodeltag}_${stip}";
+ if($irun==1) {$fname = "vert_${stirun}_xc_${sta}_${net}_${sazi1}_${eid}_${modlab}_${smodeltag}_${stip}";}
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+ #-------------------
+ # 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 - $heightxc; $oxc2 = "-Xa${x2} -Ya${y2}";
+ $x3 = $x0_xc - $wid - $xgap_xc; $y3 = $y0_xc - $heightxc; $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);
+
+ # overall label
+ $x0_label = $x3;
+ $y0_label = $y0_xc + $heightxc;
+ $olab = "-Xa${x0_label} -Ya${y0_label}";
+ $x0_label2 = $x3;
+ $y0_label2 = $y0_xc + 0.8*$heightxc;
+ $olab2 = "-Xa${x0_label2} -Ya${y0_label2}";
+
+ # map
+ $xgap_lab_map = 0.4;
+ $ygap_xc_map = 0.2;
+ #$x0_map = $x3 + $xgap_lab_map;
+ $x0_map = $x3 + $wid/2 - $widmap/2;
+ $y0_map = $y3 + $heightxc + $ygap_xc_map;
+ $omap = "-Xa${x0_map} -Ya${y0_map}";
+
+ # 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
+
+ #------------------------------------
+
+ print CSH "gmtset TICK_LENGTH $tlen_map\n"; # ticks for map
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -V $orient $omap > $psfile\n"; # START
+
+ print CSH "pscoast $Jmap $Rmap $coast_info -K -O -V $omap >> $psfile \n";
+ if ($ishelf==1) {
+ print CSH "awk '{print \$2,\$1}' ${shelf_file} | psxy $Jmap $Rmap $Wshelf -K -O -V $omap >> $psfile\n";
+ }
+
+ #print CSH "psxy ${plate_file} $Jmap $Rmap $plate_info -K -V -O $omap >> $psfile \n";
+ print CSH "psxy ${fault_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
+ #print CSH "psxy ${kcf_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
+ #print CSH "psxy ${breck_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
+ #print CSH "psxy ${fault_file2} $Jmap $Rmap $fault_infoK -K -V -O $omap >> $psfile \n";
+
+ if ($isimbox==1) {
+ $outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
+ if (not -f ${outer_boundary}) {
+ die("Check if outer_boundary ${outer_boundary} exist or not\n");
+ }
+ print CSH "psxy ${outer_boundary} $Jmap $Rmap -W1p,0/0/0 -K -O -V $omap >> $psfile\n";
+ }
+
+ # plot ray path
+ $pfile = "${dirdat}/vert_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";
+
+ # plot source and receiver lon-lat
+ $rfile = "${dirdat}/vert_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";
+ print CSH "awk '{print \$1,\$2}' ${rfile} | psxy $Jmap $Rmap $srcinfo_map -K -O -V $omap >> $psfile\n";
+ }
+
+ # plot cities
+ ##print CSH "psxy $Jmap $Rmap $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} $Jmap $Rmap -K -V -O >> $psfile \n";
+ #print CSH "awk '{print \$2,\$1,12,0,1,\"LB\",\$3}' $cities | pstext $textinfo3 -D0.05/0.05 $Jmap $Rmap -K -V -O >> $psfile \n";
+
+ #print CSH "pstext ${topo_labels} $Jmap $Rmap $textinfo -N -K -V -O >> $psfile\n";
+ #print CSH "pstext ${fault_labels} $textinfo2 $Jmap $Rmap -K -V -O >> $psfile \n";
+
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -O -V $omap >> $psfile\n";
+ print CSH "gmtset TICK_LENGTH $tlen\n"; # restore tick
+
+ #-----------------------------------------------
+
+ for ($k = 1; $k <= 3; $k ++ ) {
+ #for ($k = 1; $k <= 1; $k ++ ) {
+
+ #$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 -O -V $oxc >> $psfile\n";
+ print CSH "psbasemap $J $R $B -G${sky_color} -K -O -V $oxc >> $psfile\n";
+
+ #if ($k==1) {
+ # print CSH "psbasemap $J $R $B -K -V $orient $origin > $psfile\n"; # START
+ #} else {
+ # print CSH "psbasemap $J $R $B -K -O -V $shift >> $psfile\n";
+ #}
+
+ if ($k == 1) {
+ if ($icolor==1) {
+ #print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | pscontour -C$cptfile1 -A- -I $J $R -K -O -V $oxc >> $psfile\n";
+ #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";
+ }
+
+ 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";
+
+ # vertical exaggeration label
+ if ($ivertlab==1) {
+ $vtext = sprintf("Vertical Exaggeration = %.1f",$vexag);
+ print CSH "pstext -N $R_title $J_title -K -O -V $ov >>$psfile<<EOF\n0 0 $fsize2 0 $fontno LM $vtext\nEOF\n";
+ }
+
+
+ } elsif ($k == 2) {
+ if ($icolor==1) {
+ #print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | pscontour -C$cptfile1 -A- -I $J $R -K -O -V $oxc >> $psfile\n";
+ #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 "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";
+
+ } elsif ($k == 3) {
+ if ($icolor==1) {
+ #print CSH "awk '{print \$1,\$2,\$($icol+$k-1)}' $dfile | pscontour -C$cptfile2 -A- -I $J $R -K -O -V $oxc >> $psfile\n";
+ #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$cptfile2 $J -Q -K -O -V $oxc >> $psfile\n";
+ }
+
+ print CSH "psscale -C${cptfile2} $Dscale $Bscale2 -K -O -V $ocb >> $psfile\n";
+ #$slab1 = "ln($smodel2 / $smodel1)";
+ #print CSH "pstext -N $R_title $J_title -K -O -V $ocbarlab >>$psfile<<EOF\n0 0 $fsize2 0 $fontno LM $slab1\nEOF\n";
+ }
+
+ print CSH "psbasemap $J $R $B -K -O -V $oxc >> $psfile\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
+ if(1==1) {
+ $zminf0 = -5;
+ $finfo = "-W1.5p,0/0/0";
+ 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";
+
+ if($irun==2) {print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$saf $zmax0\n$saf $zmin0\nEOF\n";}
+ if($irun==3) {print CSH "psxy $J $R $finfo -K -O -V $oxc >>$psfile<<EOF\n$gf $zmax0\n$gf $zmin0\nEOF\n";}
+
+ }
+ #if ($p==6) {
+ # $xlarse = $saf+60;
+ # print CSH "psxy $J $R -W2p,0/0/0 -K -O -V $oxc >>$psfile<<EOF\n$saf -19\n$xlarse -21\nEOF\n"; # LARSE reflector
+ #}
+ #}
+
+ # 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 && $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";
+ 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 horizontal title above each column
+ #$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 RB $title\nEOF\n";
+
+ # # plot vertical title left of the row
+ # if ($k == 1) {
+ # print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n$xtx2 $ytx2 $fsize0 90 $fontno CM $dtitle\nEOF\n";
+ # }
+
+ } # loop over k
+
+ # 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 10 0 $fontno TC $eid\nEOF\n";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab2 >>$psfile<<EOF\n0 0 10 0 $fontno TC $stanet\nEOF\n";
+ }
+ }
+
+ # plot label on each cross section
+ for ($k = 1; $k <= 3; $k ++ ) {
+ $title = $titles[$k-1];
+ #$x3 = $x0_xc + $wid - $xgap_lab;
+ #$y3 = $y0_xc - ($k-1)*$heightxc + $ygap_lab;
+ #$talign = "RB";
+ #$x3 = $x0_xc + $xgap_lab;
+ #$y3 = $y0_xc - ($k-1)*$heightxc + $ygap_lab;
+ $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 12 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 "ghostview $psfile &\n";}
+ if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+
+} # loop over p
+
+#==================================================
+# figure with map with m16 beneath
+
+if ($imfinal == 1) {
+
+ #print CSH "gmtset COLOR_NAN $sky_color 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";
+
+ @iBs = (5,8,4);
+
+ for ($p = $pmin; $p <= $pmax; $p ++ ) {
+
+ $stip = sprintf("%3.3i",$p);
+ $fname = "${fname}_single";
+
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+ # datafile -- COLUMNS: lon, lat, m00, m16
+ $dfile = "${dirdat}/vert_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";
+ 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 = $x2;
+ $y0_label = $y0_xc + $heightxc;
+ $olab = "-Xa${x0_label} -Ya${y0_label}";
+
+ # map
+ $xgap_lab_map = 0.4;
+ $ygap_xc_map = 0.2;
+ $x0_map = $x2 + $wid/2 - $widmap/2;
+ $y0_map = $y2 + $heightxc + $ygap_xc_map;
+ $omap = "-Xa${x0_map} -Ya${y0_map}";
+
+ #---------------------
+
+ print CSH "gmtset TICK_LENGTH $tlen_map\n"; # ticks for map
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -V $orient $omap > $psfile\n"; # START
+
+ print CSH "pscoast $Jmap $Rmap $coast_info -K -O -V $omap >> $psfile \n";
+ if ($ishelf==1) {print CSH "awk '{print \$2,\$1}' ${shelf_file} | psxy $Jmap $Rmap $Wshelf -K -O -V $omap >> $psfile\n";}
+
+ print CSH "psxy ${fault_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
+ #print CSH "psxy ${kcf_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
+
+ if ($isimbox==1) {
+ $outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
+ if (not -f ${outer_boundary}) {die("Check if outer_boundary ${outer_boundary} exist or not\n");}
+ print CSH "psxy ${outer_boundary} $Jmap $Rmap -W1p,0/0/0 -K -O -V $omap >> $psfile\n";
+ }
+
+ # plot ray path
+ $pfile = "${dirdat}/vert_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";
+
+ # plot source and receiver lon-lat
+ $rfile = "${dirdat}/vert_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 $omap >> $psfile\n";
+ print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmap $Rmap $recinfo_map -K -O -V $omap >> $psfile\n";
+
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -O -V $omap >> $psfile\n";
+
+ # plot additional ray paths
+ if ($p==10) {
+ $x0_mapinset = $x0_map + 2;
+ $y0_mapinset = $y0_map + 0;
+ $omapinset = "-Xa${x0_mapinset} -Ya${y0_mapinset}";
+
+ $Rmapinset = "-R-118.6/-117/33.2/34.5";
+ $Jmapinset = "-JM1i";
+ $xlon1 = -118; $xlat1 = 33.9;
+ $xlon2 = -118.3; $xlat2 = 33.8;
+ $rayinfo = "-W1.0p,0/0/0";
+ print CSH "psbasemap $Jmapinset $Rmapinset $Bmap -K -O -V $omapinset >> $psfile\n";
+ print CSH "pscoast $Jmapinset $Rmapinset $coast_info2 -K -O -V $omapinset >> $psfile \n";
+
+ print CSH "psxy $Jmapinset $Rmapinset $rayinfo -K -O -V $omapinset >>$psfile<<EOF\n$slon $slat\n\n$rlon $rlat\nEOF\n";
+ print CSH "psxy $Jmapinset $Rmapinset $rayinfo -K -O -V $omapinset >>$psfile<<EOF\n$slon $slat\n$xlon1 $xlat1\n$rlon $rlat\nEOF\n";
+ print CSH "psxy $Jmapinset $Rmapinset $rayinfo -K -O -V $omapinset >>$psfile<<EOF\n$slon $slat\n$xlon2 $xlat2\n$rlon $rlat\nEOF\n";
+
+ print CSH "psxy ${fault_file} $Jmapinset $Rmapinset $fault_info -K -V -O $omapinset >> $psfile \n";
+ #print CSH "psxy ${kcf_file} $Jmapinset $Rmapinset $fault_info -K -V -O $omapinset >> $psfile \n";
+ print CSH "awk '{print \$1,\$2}' ${rfile} | psxy $Jmapinset $Rmapinset $srcinfo_mapinset -K -O -V $omapinset >> $psfile\n";
+ print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmapinset $Rmapinset $recinfo_mapinset -K -O -V $omapinset >> $psfile\n";
+ print CSH "psbasemap $Jmapinset $Rmapinset $Bmap -K -O -V $omapinset >> $psfile\n";
+ }
+
+ print CSH "gmtset TICK_LENGTH $tlen\n"; # restore tick
+
+ #-----------------------------------------------
+
+ for ($k = 2; $k <= 2; $k ++ ) {
+
+ #$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 -O -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 "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";
+ #if ($p==6) {
+ # $xlarse = $saf+60;
+ # print CSH "psxy $J $R -W2p,0/0/0 -K -O -V $oxc >>$psfile<<EOF\n$saf -19\n$xlarse -21\nEOF\n"; # LARSE reflector
+ #}
+ #}
+
+ # plot source and receiver
+ 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 = 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 $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";}
+ #}
+
+ } # loop over k
+
+ # 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 {
+ $textinfo = "-N";
+ print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 12 0 $fontno TR $eid to $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 12 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 "ghostview $psfile &\n";
+ }
+ if ($ipdf==1) {
+ print CSH "ps2pdf $psfile\n";
+ }
+
+ } # loop over p
+
+} # imfinal == 1
+
+#==================================================
+
+#------------------------------------
+
+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("ghostview $psfile &");}
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,284 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_vert_models_basemap.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_basemap.pl 1
+#
+#==========================================================
+
+if (@ARGV < 1) {die("Usage: plot_vert_models_basemap.pl xxx\n");}
+ at pinds = @ARGV;
+$nump = @pinds;
+$stp = sprintf("%3.3i",$nump);
+
+$irun = 1;
+$icolor = 1;
+$ivertlab = 0; # vertical exaggeration label
+$isimbox = 1; # box around simulation region
+$iletter = 0; # label for publication figures
+$letter = "A";
+$imfinal = 0; # plot a 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 = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/vert_${stirun}";
+
+
+#---------------------------------------------------------
+
+# subtitles
+#@titles = ("$mtit $smodel1","$mtit $smodel2","ln(${smodel2} / ${smodel1})");
+
+# parameters for color scales
+#$cpert1 = 0.2;
+#$cpert2 = 0.2;
+$scale_color = 65.0;
+$colorbar = "seis";
+
+# KEY: resolution of color plots
+$interp = "-I0.5/0.5 -S1.5"; # 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/socal_faults_xc.lonlat";
+#$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_tomo = "/net/denali/scratch1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO";
+$dir_sources = "/net/sierra/raid1/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.2;
+#$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";
+
+$iportrait = 1;
+if($iportrait == 1){$orient = "-P"; $rotangle = 0}
+else {$orient = " "; $rotangle = 90}
+
+# plotting specifications
+$fsize0 = "14";
+$fsize1 = "12";
+$fsize2 = "11";
+$fsize3 = "10";
+$fontno = "1"; # 1 or 4
+$tlen = "6p";
+$fpen = "1.5p";
+$tpen = "1.5p";
+
+# 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";
+$fault_info2 = "-M -W2p,255/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";
+$srcinfo_map = "-Sa20p $colinfo_map -N";
+$recinfo_map = "-Si20p $colinfo_map -N";
+$srcinfo_mapinset = "-Sa10p -W0.5p,0/0/0 -G0/255/255 -N";
+$recinfo_mapinset = "-Si10p -W0.5p,0/0/0 -G0/255/255 -N";
+$colinfo_xc = "-W2p,0/0/0";
+$srcinfo_xc = "-Sa20p $colinfo_xc -N";
+$recinfo_xc = "-Si20p $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;
+
+# plot title
+$J_title = "-JM1";
+$R_title = "-R0/1/0/1";
+$fsize_title = $fsize0;
+
+#---------------------------------------------------------
+
+# load values and check that the files exist
+($nump0,undef,undef) = split(" ",`ls -1 ${dirdat}/*path | wc`);
+if($nump0 == 0) {die("No ray paths available")}
+#if($pmax > $nump) {$pmax = $nump;}
+#if($pmin < 1) {$pmin = 1;}
+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";
+ if (not -f $ffile) {die("Check if ffile $ffile exist or not\n");}
+ open(IN,$ffile); @faults = <IN>; close(IN);
+
+ $widmap = 6.5; # width of map
+ $Jmap = "-JM${widmap}";
+ $B0map = "-Ba${tick1}f${tick2}d:.\" \":";
+ #$Bmap = "$B0map".$Bopts[7];
+ $Bmap = "$B0map".$Bopts[0];
+
+ # parameters controlling the placement of subplots and labels
+ $omap = "-Xa1 -Ya3";
+
+ $fname = "vert_${stirun}_basemap_${stp}";
+ $psfile = "$fname.ps"; $jpgfile = "$fname.jpg";
+
+#==============================================
+
+$cshfile = "plot_vert_models_basemap.csh";
+open(CSH,">$cshfile");
+print CSH "gmtset COLOR_NAN $sky_color PAPER_MEDIA letter 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";
+
+ print CSH "psbasemap $Jmap $Rmap $Bmap -K -V $orient $omap > $psfile\n"; # START
+
+ print CSH "pscoast $Jmap $Rmap $coast_info -K -O -V $omap >> $psfile \n";
+ if ($ishelf==1) {
+ print CSH "awk '{print \$2,\$1}' ${shelf_file} | psxy $Jmap $Rmap $Wshelf -K -O -V $omap >> $psfile\n";
+ }
+
+ print CSH "psxy ${fault_file} $Jmap $Rmap $fault_info -K -V -O $omap >> $psfile \n";
+ print CSH "psxy ${fault_file2} $Jmap $Rmap $fault_info2 -K -V -O $omap >> $psfile \n";
+
+ if ($isimbox==1) {
+ $outer_boundary = "${dir_sources}/SPECFEM_outer_points.dat";
+ if (not -f ${outer_boundary}) {
+ die("Check if outer_boundary ${outer_boundary} exist or not\n");
+ }
+ print CSH "psxy ${outer_boundary} $Jmap $Rmap -W2p,0/0/0 -K -O -V $omap >> $psfile\n";
+ }
+
+#==============================================
+
+# modify the height of each cross-section based on the length
+#for ($p = $pmin; $p <= $pmax; $p ++ ) {
+for ($ip = 1; $ip <= $nump; $ip ++ ) {
+ $p = $pinds[$ip-1];
+ $stip = sprintf("%3.3i",$p);
+ print "-- $ip -- $nump -- $p --\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";
+ 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]);
+
+ # positions of six faults along each profile
+ # (1) SAF, (2) GF, (3) SGF, (4) MCF, (5) SYF, (6) CRF
+ (undef,$eid,$stanet,$azi1,$azi2,$saf,$gf,$sgf,$mcf,$syf,$crf) = split(" ",$faults[$p-1]);
+ ($sta,$net) = split("\\.",$stanet);
+ $sazi1 = sprintf("%3.3i",$azi1);
+ $sazi2 = sprintf("%3.3i",$azi2);
+ print "-- $stanet -- $sta -- $net -- $eid -- $sazi1 -- $sazi2 --\n";
+
+ #---------------------------------------------------------
+
+ # plot ray path
+ $pfile = "${dirdat}/vert_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";
+
+ # plot source and receiver lon-lat
+ $rfile = "${dirdat}/vert_xc_${stip}_A_B";
+ if (not -f $rfile) {
+ die("Check if rfile $rfile exist or not\n");
+ }
+
+ # plot source
+ print CSH "awk '{print \$1,\$2}' ${rfile} | psxy $Jmap $Rmap $srcinfo_map -K -O -V $omap >> $psfile\n";
+
+ # plot receiver
+ print CSH "awk '{print \$5,\$6}' ${rfile} | psxy $Jmap $Rmap $recinfo_map -K -O -V $omap >> $psfile\n";
+
+ # plot cities
+ ##print CSH "psxy $Jmap $Rmap $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} $Jmap $Rmap -K -V -O >> $psfile \n";
+ #print CSH "awk '{print \$2,\$1,12,0,1,\"LB\",\$3}' $cities | pstext $textinfo3 -D0.05/0.05 $Jmap $Rmap -K -V -O >> $psfile \n";
+
+ #print CSH "pstext ${topo_labels} $Jmap $Rmap $textinfo -N -K -V -O >> $psfile\n";
+ #print CSH "pstext ${fault_labels} $textinfo2 $Jmap $Rmap -K -V -O >> $psfile \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 {
+# $textinfo = "-N";
+# print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab >>$psfile<<EOF\n0 0 10 0 $fontno TC $eid\nEOF\n";
+# print CSH "pstext $R_title $J_title $textinfo -K -O -V $olab2 >>$psfile<<EOF\n0 0 10 0 $fontno TC $stanet\nEOF\n";
+# }
+
+} # loop over p
+
+ # fault labels
+ $fault_labs = "$dir0/faults/socal_fault_labs.lonlat";
+ $textinfo = "-N -C4p -W255/255/255o,1.0p,255/0/0,solid -G0/0/0";
+ print CSH "awk '{print \$1,\$2,10,0,1,\"CM\",\$3}' ${fault_labs} | pstext $textinfo $Jmap $Rmap -K -V -O $omap >> $psfile\n";
+
+ print CSH "pstext -N $R_title $J_title -K -O -V >>$psfile<<EOF\n 10 10 16 0 $fontno CM \nEOF\n";
+ print CSH "psbasemap $Jmap $Rmap $Bmap -O -V $omap >> $psfile\n"; # FINISH (no K)
+
+ #if($ixv==1) {print CSH "convert $psfile -rotate $rotangle $jpgfile\n";}
+ if ($ixv==1) {print CSH "ghostview $psfile &\n";}
+ if ($ipdf==1) {print CSH "ps2pdf $psfile\n";}
+
+#------------------------------------
+
+close (CSH);
+system("csh -f $cshfile");
+
+system("ghostview $psfile &")
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_vert_models_basemap.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_xc_seis_all.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_xc_seis_all.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_xc_seis_all.pl 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,138 @@
+#!/usr/bin/perl -w
+
+#==========================================================
+#
+# plot_xc_seis_all.pl
+# Carl Tape
+# 07-Dec-2008
+#
+# EXAMPLES:
+# plot_xc_seis_all.pl T006_T030 1 1
+# plot_xc_seis_all.pl T006_T030 0 1
+# plot_xc_seis_all.pl T003_T030 1 1
+# plot_xc_seis_all.pl T003_T030 0 1
+# plot_xc_seis_all.pl T002_T030 1 1
+# plot_xc_seis_all.pl T002_T030 0 1
+#
+#==========================================================
+
+if (@ARGV < 3) {die("Usage: plot_xc_seis_all.pl xxx\n");}
+($Ttag,$iwin,$i1D) = @ARGV;
+
+# make output directory
+$otag = "${Ttag}_win${iwin}_1D${i1D}";
+$odir = "./$otag";
+#`rm -rf $odir`;
+`mkdir -p $odir`;
+
+# directory containing the data files
+$irun = 1;
+$stirun = sprintf("%2.2i",$irun);
+
+$dirdat = "/net/sierra/raid1/carltape/results/PLOTTING/SAVE/vert_${stirun}";
+$seisdir = "/net/denali/raid1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot_work/${otag}";
+$xcdir = "/net/sierra/raid1/carltape/results/PLOTTING/gmt";
+if (not -e $dirdat) {die("Check if dirdat $dirdat exist or not\n");}
+if (not -e $xcdir) {die("Check if xcdir $xcdir exist or not\n");}
+if (not -e $seisdir) {die("Check if seisdir $seisdir exist or not\n");}
+
+#---------------------------------------------------------
+
+# load fault positions along each profile
+$ffile = "$dirdat/ALL_rays_fault_positions_mod";
+if (not -f $ffile) {die("Check if ffile $ffile exist or not\n");}
+open(IN,$ffile); @faults = <IN>; close(IN);
+$nump = @faults;
+
+for ($p = 1; $p <= $nump; $p ++ ) {
+ (undef,$eid,$stanet,$azi1,$azi2,$saf,$gf,$sgf,$mcf,$syf,$crf,$ef,$snf,$kcf) = split(" ",$faults[$p-1]);
+ print "-- $p -- $stanet -- $eid --\n";
+}
+#die("TESTING");
+
+$pmin = 1; $pmax = $nump;
+#$pmin = 1; $pmax = 10;
+#$pmin = 77; $pmax = $pmin;
+
+for ($p = $pmin; $p <= $pmax; $p ++ ) {
+ $stp = sprintf("%3.3i",$p);
+
+ # 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 "-- $p -- $stanet -- $sta -- $net -- $eid -- $sazi1 -- $sazi2 --\n";
+
+ # length of cross section
+ $rayfile = "$dirdat/vert_xc_${stp}_ray_path";
+ if (not -f $rayfile) {die("Check if rayfile $rayfile exist or not\n");}
+ (undef,undef,undef,undef,undef,undef,undef,undef,$dmin,$dmax) = split(" ",`minmax -C -I2 $rayfile`);
+ $dran = ($dmax - $dmin)/1000;
+ $heightxc = $heightxc0;
+ if ($dran < 350) {$heightxc = "width=4.5cm";}
+ else {$heightxc = "height=12cm";}
+ print "length of cross section is $dran km\n";
+
+ # find cross section file
+ @xcfiles = glob("${xcdir}/*$sta*$net*$eid*eps");
+ $nxcfiles = @xcfiles;
+ if($nxcfiles != 1) {
+ print "$ftag\n @xcfiles\n $nxcfiles\n";
+ die("Must have only one cross section file");
+ }
+ $xcfile = $xcfiles[0];
+ #$xcfile = "${xcdir}/vert_01_xc_LDF_CI_269_9983429_vs_m16_m00_001.eps";
+ if (not -f $xcfile) {die("Check if xcfile $xcfile exist or not\n");}
+
+ # find seismofile
+ $stag = "seis_2";
+ if($i1D==1) {$stag = "seis_3";}
+ $ftag = "${seisdir}/${eid}_${Ttag}_${sta}_${net}*${stag}.eps";
+ @seisfiles = glob($ftag);
+ $nseisfiles = @seisfiles;
+ if($nseisfiles != 1) {
+ print "$ftag\n @seisfiles\n $nseisfiles\n";
+ die("Must have only one seismogram file");
+ }
+ $seisfile = $seisfiles[0];
+ #$seisfile = "${seisdir}/9983429_T006_T030_LDF_CI_m16_seis_2.eps";
+ if (not -f $seisfile) {die("Check if seisfile $seisfile exist or not\n");}
+
+ # write latex file for combining figures
+ $name = "xc_and_seis";
+ $texfile = "${name}.tex";
+ open(TEX,">$texfile");
+print TEX "\\documentclass[pdf,mpa]{prosper}\n";
+print TEX "\\usepackage{portland}\n";
+print TEX "\\begin{document}\n";
+print TEX "\\begin{slide}{}\n";
+print TEX "\\vspace{-1.7cm}\n";
+print TEX "\\fcolorbox{white}{white}{\n";
+print TEX "\\hspace{-1cm} \n";
+print TEX "\\begin{tabular}{c}\n";
+#print TEX "\\includegraphics[width=3.8cm,angle=-90]{$xcfile}\n";
+print TEX "\\includegraphics[${heightxc},angle=-90]{$xcfile}\n";
+print TEX "\\\\\n";
+print TEX "\\includegraphics[width=3.8cm,angle=-90]{$seisfile}\n";
+print TEX "\\end{tabular}\n";
+print TEX "}\n";
+print TEX "\\end{slide}\n";
+print TEX "\\end{document}\n";
+close(TEX);
+
+ # execute latex file and convert to pdf
+ $fname = "xc_and_seis_${eid}_${sta}_${net}_${Ttag}";
+ `latex ${name}`;
+ `dvips -o ${name}.ps ${name}`;
+ `mv ${name}.ps ${fname}.ps`;
+ `ps2pdf ${fname}.ps`;
+
+}
+
+# move all pdf and ps files to output directory
+`mv xc_and_seis*pdf xc_and_seis*ps $odir`;
+`/home/carltape/bin/pdcat -r $odir/xc_and_seis*pdf $odir/ALL_xc_and_seis_${otag}.pdf`;
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/plot_xc_seis_all.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/ps2eps.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/ps2eps.pl (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/ps2eps.pl 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,25 @@
+#!/usr/bin/perl -w
+
+$irun = 3;
+$stip = sprintf("%2.2i",$irun);
+
+ at epsfiles = glob("vert_${stip}_xc*.ps");
+$neps = @epsfiles;
+for ($k = 1; $k <= $neps; $k = $k+1) {
+ $efile = $epsfiles[$k-1];
+ `ps2eps -f $efile`;
+}
+
+if ($irun==1) {
+ `ps2ps vert_01_xc_HEC_CI_282_9968977_vs_m16_m00_003.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_HEC_CI_282_9968977_vs_m16_m00_003.eps`;
+ `ps2ps vert_01_xc_LCG_CI_122_10215753_vs_m16_m00_021.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_LCG_CI_122_10215753_vs_m16_m00_021.eps`;
+ `ps2ps vert_01_xc_LGU_CI_353_10097009_vs_m16_m00_067.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_LGU_CI_353_10097009_vs_m16_m00_067.eps`;
+ `ps2ps vert_01_xc_BZN_AZ_306_14138080_vs_m16_m00_013.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_BZN_AZ_306_14138080_vs_m16_m00_013.eps`;
+ `ps2ps vert_01_xc_EDW2_CI_134_14236768_vs_m16_m00_011.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_EDW2_CI_134_14236768_vs_m16_m00_011.eps`;
+ `ps2ps vert_01_xc_ISA_CI_160_14383980_vs_m16_m00_049.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_ISA_CI_160_14383980_vs_m16_m00_049.eps`;
+ `ps2ps vert_01_xc_HEC_CI_296_14418600_vs_m16_m00_007.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_HEC_CI_296_14418600_vs_m16_m00_007.eps`;
+ `ps2ps vert_01_xc_ISA_CI_113_14418600_vs_m16_m00_044.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_ISA_CI_113_14418600_vs_m16_m00_044.eps`;
+ `ps2ps vert_01_xc_DPP_CI_314_9753485_vs_m16_m00_079.ps temp.ps ; ps2eps temp.ps ; mv temp.eps vert_01_xc_DPP_CI_314_9753485_vs_m16_m00_079.eps`;
+}
+
+#==================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/ps2eps.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/xc_and_seis.tex
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/xc_and_seis.tex (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/gmt/xc_and_seis.tex 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,15 @@
+\documentclass[pdf,mpa]{prosper}
+\usepackage{portland}
+\begin{document}
+\begin{slide}{}
+\vspace{-1.7cm}
+\fcolorbox{white}{white}{
+\hspace{-1cm}
+\begin{tabular}{c}
+\includegraphics[width=4.5cm,angle=-90]{/net/sierra/raid1/carltape/results/PLOTTING/gmt/vert_01_xc_DSC_CI_214_14155260_vs_m16_m00_086.eps}
+\\
+\includegraphics[width=3.8cm,angle=-90]{/net/denali/raid1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/misfit_plot_work/T002_T030_win0_1D0/14155260_T002_T030_DSC_CI_m16_seis_2.eps}
+\end{tabular}
+}
+\end{slide}
+\end{document}
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,35 @@
+%
+% function
+% Carl Tape, 01-Aug-2007
+% printed xxx
+%
+% This function converts between a bounding region in lat-lon and in UTM.
+%
+% EXAMPLE:
+
+% ax_utm = axes_utm2ll([-121.6 -114.7 32.2 36.8],11,0)
+%
+%
+% calls utm2ll.m
+% called by xxx
+%
+
+function ax_out = axes_utm2ll(ax_in,i_zone,i_type)
+
+xmin0 = ax_in(1);
+xmax0 = ax_in(2);
+ymin0 = ax_in(3);
+ymax0 = ax_in(4);
+
+[xmin,ymin] = utm2ll(xmin0,ymin0,i_zone,i_type);
+[xmax,ymax] = utm2ll(xmax0,ymax0,i_zone,i_type);
+
+ax_out = [xmin xmax ymin ymax];
+
+if 0==1
+ ax_ll = [-121.6 -114.7 32.2 36.8]
+ ax_utm = axes_utm2ll(ax_ll,11,0)
+ axes_utm2ll(ax_utm,11,1)
+end
+
+%==========================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,38 @@
+%
+% function [xbox,ybox] = boxpts(ax0,n)
+% CARL TAPE, 06-Sept-2005
+% printed xxx
+%
+% This function inputs an axes box and outputs a set of points describing
+% the boundary of the box.
+%
+% EXAMPLE:
+% ax0 = [-121.6 -114.7 32.2 36.8]; [xbox,ybox] = boxpts(ax0,100); figure, plot(xbox,ybox,'.');
+%
+% calls xxx
+% called by xxx
+%
+
+function [xbox,ybox] = boxpts(ax0,n)
+
+xmin = ax0(1);
+xmax = ax0(2);
+ymin = ax0(3);
+ymax = ax0(4);
+
+% total length of box
+xran = xmax - xmin;
+yran = ymax - ymin;
+len = 2*xran + 2*yran;
+
+% spacing increment
+dx = len/n;
+nx = round(xran/dx);
+ny = round(yran/dx);
+xpts = linspace(xmin,xmax,nx);
+ypts = linspace(ymin,ymax,ny);
+
+xbox = [xpts ones(1,ny)*xmax fliplr(xpts) ones(1,ny)*xmin ]';
+ybox = [ones(1,nx)*ymin ypts ones(1,nx)*ymax fliplr(ypts) ]';
+
+%===================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,32 @@
+%
+% function out_array = merge_arrays(array1,array2,stconnect)
+% CARL TAPE, 08-Aug-2007
+% printed xxx
+%
+% This file inputs two arrays of strings and merges them with a connecting
+% character.
+%
+% See also split_array.m, which is the opposite.
+%
+% calls xxx
+% called by xxx
+%
+
+function out_array = merge_arrays(array1,array2,stconnect)
+
+num = length(array1);
+if length(array1) ~= length(array2)
+ error(' array1 and array2 must have the same length');
+end
+
+% allocate array; initialize array
+% THIS SPEEDS THINGS UP BY A FACTOR OF ABOUT 100.
+out_array = cellstr(repmat(' ',num,1));
+
+for ii = 1:num
+ %if mod(ii,500)==0, disp([num2str(ii) ' out of ' num2str(num)]); end
+ out_array{ii} = [array1{ii} stconnect array2{ii}];
+end
+out_array = out_array(:);
+
+%=======================================================================
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,857 @@
+%
+% socal_tomo_plots.m
+% CARL TAPE, 27-March-2009
+% printed xxx
+%
+% This file generates xyz points for which we want to extract an
+% interpolation value from a SPECFEM volumetric field. Once we have the
+% function value, we can plot the cross-section in GMT.
+%
+% calls merge_arrays.m, xxx
+% called by xxx
+%
+
+clc
+clear
+close all
+format short, format compact
+
+dir0 = '/net/denali/raid1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/';
+pdir = [dir0 'iterate_adj/model_plot'];
+mdir = [dir0 'OUTPUT_SUBSPACE_m16/'];
+%mdir = [dir0 'iterate_adj/matlab_work/OUTPUT_SUBSPACE/'];
+
+% KEY: do you want vertical and/or horizontal cross-sections
+ivert = 0;
+iplot_vert = 1;
+
+ihorz = 0;
+iplot_horz = 0;
+iplot_horz_coverage = 0;
+
+iwrite = 1;
+ifigure = 0;
+
+%-----------------------------------
+
+% get a set of plotting gridpoints to extract a function value at a
+% given x,y,z from a set of .bin files
+
+% bounds for SPECFEM3D simulations
+ax0 = [-121.6 -114.7 32.2 36.8];
+ax0_utm = axes_utm2ll(ax0,11,0);
+%ax0_utm = [ 0.3271 0.4698 3.6792 3.7914] * 1e6; % zoom
+xmin = ax0_utm(1); xmax = ax0_utm(2);
+ymin = ax0_utm(3); ymax = ax0_utm(4);
+
+ % outline of SPECFEM region
+[xbox_utm,ybox_utm] = boxpts(ax0_utm,500);
+[xbox_ll,ybox_ll] = utm2llvec(xbox_utm,ybox_utm,11,1);
+
+figure; plot(xbox_utm, ybox_utm, '.');
+
+%-----------------------------------
+% HORIZONTAL CROSS-SECTIONS
+
+if ihorz == 1
+ % horizontal grid
+ numx = 300; % key command for resolution
+ xvec = linspace(xmin,xmax,numx);
+ dx = xvec(2) - xvec(1);
+ yvec = [ymin:dx:ymax];
+ numy = length(yvec);
+ [X,Y] = meshgrid(xvec,yvec);
+ Xvec = X(:);
+ Yvec = Y(:);
+ N = length(Xvec);
+
+ % convert back to lat-lon
+ disp('Entering utm2llvec.m...');
+ [Xvec_ll,Yvec_ll] = utm2llvec(Xvec,Yvec,11,1);
+ figure; plot(Xvec_ll, Yvec_ll, '.','markersize',2);
+
+ disp(sprintf('horizontal mesh is %i (x) by %i (y) for %i points total',numx,numy,N));
+ disp(sprintf('horizontal spacing of plotting points is %.2f meters',dx));
+
+ % set of elevation values, in meters
+ depvec = 1e3 * [ 0:-1:-40];
+ %depvec = 1e3 * [-1];
+ ndep = length(depvec);
+
+ %----------------------------
+
+ if iwrite==1
+ % list of depths for each index
+ fid = fopen([pdir 'horz_xc_elevation_cuts'],'w');
+ for ii=1:ndep
+ fprintf(fid,'%10i%16i\n',ii,depvec(ii));
+ end
+ fclose(fid);
+
+ % reference horizontal surface
+ fid = fopen([pdir 'horz_xc_lonlat'],'w');
+ for kk=1:N
+ fprintf(fid,'%16.6e%16.6e%16.6e%16.6e\n',Xvec_ll(kk),Yvec_ll(kk),Xvec(kk),Yvec(kk));
+ end
+ fclose(fid);
+
+ % xyz files
+ for ii=1:ndep
+ zdep = depvec(ii)
+ filename = sprintf('horz_xc_elevation_%3.3i',ii);
+ fid = fopen([pdir filename],'w');
+ for kk=1:N
+ fprintf(fid,'%16.6e%16.6e%16.6e\n',Xvec(kk),Yvec(kk),zdep);
+ end
+ fclose(fid);
+ end
+ end
+end
+
+%----------------------------
+% VERTICAL CROSS-SECTIONS
+
+if ivert ==1
+
+ % read in sources and receivers
+ recfile = '/home/carltape/gmt/stations/seismic/Matlab_output/STATIONS_CALIFORNIA_TOMO_INNER_specfem';
+ [rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(recfile);
+ stanet = merge_arrays(stnm,netwk,'.');
+ srcfile = '/net/sierra/raid1/carltape/results/SOURCES/socal_16/EIDs_lonlat_otime';
+ [junk1,eid,junk2,junk3,sdep,junk5,junk6,slon,slat] = textread(srcfile,'%f%f%s%s%f%f%f%f%f');
+
+ irun = 3; % KEY: which set of cross-sections to make
+ dir1 = [pdir 'SAVE/vert_' sprintf('%2.2i',irun) '/']; % KEY COMMAND
+
+ switch irun
+
+ case 1 % 1: PICK SOURCE-RECEIVER PAIRS
+ % specify the src-rec pairs and find their lon,lat,dep
+ eids = [9983429 14096196 9968977 14179736 14383980 9818433 14418600 9967901 ...
+ 14383980 9703873 14236768 10992159 14138080 14186612 14186612 9688709 ...
+ 14077668 14095540 9627721 14236768 10215753 10215753 14096196 14096196 ...
+ 9983429 10006857 9173365 14095540 14077668 14095628 9674049 3321590 ...
+ 3320736 14179292 9674049 14418600 9173365 9155518 9718013 14179736 ...
+ 9154092 14169456 9045109 14418600 3320736 13935988 14179292 14236768 ...
+ 14383980 9695397 9173365 9882325 14077668 10006857 9967901 9695397 ...
+ 14155260 14077668 14096196 9818433 10006857 10215753 14151344 9173365 ...
+ 13935988 14155260 10097009 9882329 14169456 9674213 9753485 9096972 ...
+ 14000376 10370141 14383980 9093975 12659440 14077668 9753485 14155260 ...
+ 14155260 14155260 14155260 14155260 14155260 14155260];
+ recs = {'LDF.CI','WGR.CI','HEC.CI','LAF.CI','STC.CI','CLC.CI','HEC.CI','OSI.CI',...
+ 'SMS.CI','RVR.CI','EDW2.CI','SMM.CI','BZN.AZ','PSR.CI','FMP.CI','BRE.CI',...
+ 'LDF.CI','LCG.CI','FMP.CI','WSS.CI','LCG.CI','HAST.TA','SMM.CI','SCI2.CI',...
+ 'DAN.CI','LGB.CI','FRD.AZ','SNCC.CI','SNCC.CI','SNCC.CI','SNCC.CI','SNCC.CI',...
+ 'SNCC.CI','SNCC.CI','LDF.CI','LDF.CI','LDF.CI','LDF.CI','LDF.CI','LDF.CI',...
+ 'LDF.CI','ISA.CI','ISA.CI','ISA.CI','ISA.CI','ISA.CI','ISA.CI','ISA.CI',...
+ 'ISA.CI','ISA.CI','ISA.CI','ISA.CI','ISA.CI','ISA.CI','ISA.CI','CWC.CI',...
+ 'SCZ2.CI','SCZ2.CI','SCZ2.CI','SCZ.CI','LGU.CI','LGU.CI','LGU.CI','LGU.CI',...
+ 'LGU.CI','LGU.CI','LGU.CI','LGU.CI','LGU.CI','LGU.CI','LGU.CI','LGU.CI',...
+ 'LGU.CI','LGU.CI','LGU.CI','LGU.CI','LGU.CI','LGU.CI','DPP.CI','TEH.CI',...
+ 'LDR.CI','LRL.CI','CLC.CI','SLA.CI','LMR2.CI','DSC.CI'};
+
+ slons = zeros(length(eids),1);
+ slats = zeros(length(eids),1);
+ sdeps = zeros(length(eids),1);
+ nrays = length(eids);
+ for ii=1:nrays
+ imatch = find(eids(ii) == eid);
+ slons(ii) = slon(imatch);
+ slats(ii) = slat(imatch);
+ sdeps(ii) = sdep(imatch);
+ end
+ rlons = zeros(nrays,1); rlats = zeros(nrays,1); rdeps = zeros(nrays,1);
+ for ii=1:nrays
+ imatch = find(strcmp(recs(ii),stanet)==1);
+ rlons(ii) = rlon(imatch);
+ rlats(ii) = rlat(imatch);
+ end
+
+ % specify parameters for sampling
+ preseg_vec = 100*ones(nrays,1);
+ postseg_vec = 100*ones(nrays,1);
+
+ case 2 % 2: PICK SOURCE-RECEIVER PAIRS
+
+ % select profiles perpendicular to the SAF -- saf.m
+ saflonlat_cuts = [
+ -120.9420 36.3617 -118.4337 38.1771
+ -120.7410 36.1754 -118.2194 37.9734
+ -120.5640 36.0081 -117.9244 37.6869
+ -120.3880 35.8542 -117.9540 37.7197
+ -120.2380 35.6755 -117.8501 37.5748
+ -119.9610 35.4046 -117.3811 37.1230
+ -119.7050 35.1649 -117.3195 37.0532
+ -119.5490 35.0317 -117.3772 37.0805
+ -119.4150 34.9367 -117.4435 37.1133
+ -119.2290 34.8725 -118.1995 37.4393
+ -119.0590 34.8386 -117.5635 37.2505
+ -118.9380 34.8158 -118.0720 37.4215
+ -118.7480 34.7690 -117.7609 37.3463
+ -118.5930 34.7257 -117.2598 37.1985
+ -118.4080 34.6632 -117.0520 37.1275
+ -118.1260 34.5614 -116.3847 36.8594
+ -117.8940 34.4664 -116.2315 36.8018
+ -117.6070 34.3428 -115.8980 36.6545
+ -117.4670 34.2764 -115.7606 36.5887
+ -117.2620 34.1669 -115.5181 36.4592
+ -116.9810 34.0651 -115.5469 36.4959
+ -116.7950 33.9496 -116.7731 36.6475
+ -116.5660 33.9238 -115.3337 36.4268
+ -116.3040 33.8099 -114.4192 36.0215
+ -116.0630 33.6362 -113.9734 35.7163
+ -115.8780 33.4662 -113.7552 35.5201
+ -115.6280 33.2626 -113.3162 35.1667
+ -115.5940 33.0397 -113.1488 34.8184
+ -115.4120 32.7096 -112.9413 34.4540 ];
+ nsafcuts = length(saflonlat_cuts);
+
+ slons = [-118.1 -118.1 -118.1 saflonlat_cuts(:,1)'];
+ slats = [33.9 33.9 33.9 saflonlat_cuts(:,2)'];
+ rlons = [-117 -119.4 -120 saflonlat_cuts(:,3)'];
+ rlats = [35 35.4 35 saflonlat_cuts(:,4)'];
+ nrays = length(slons);
+ preseg_vec = [100 150 150 200*ones(1,nsafcuts)];
+ postseg_vec = [100 150 150 0*ones(1,nsafcuts)];
+ sdeps = zeros(nrays,1);
+
+ case 3 % 2: PICK SOURCE-RECEIVER PAIRS
+
+ % select profiles perpendicular to the Garlock -- saf.m
+ garlonlat_cuts = [
+ -118.7309 34.8803 -117.4704 33.4150
+ -118.5238 34.9737 -117.6324 33.3333
+ -118.1984 35.1372 -116.8821 33.7032
+ -118.0406 35.2404 -116.5687 33.9113
+ -117.8138 35.3766 -116.5167 33.9286
+ -117.4834 35.4934 -116.6783 33.8211
+ -117.1826 35.5674 -116.6713 33.8186
+ -116.9706 35.6063 -116.9243 33.8081
+ -116.6402 35.5966 -116.6391 33.7980
+ -116.4528 35.5966 -116.2956 33.8026 ];
+ ngarcuts = length(garlonlat_cuts);
+
+ slons = garlonlat_cuts(:,1)';
+ slats = garlonlat_cuts(:,2)';
+ rlons = garlonlat_cuts(:,3)';
+ rlats = garlonlat_cuts(:,4)';
+ nrays = length(slons);
+ preseg_vec = 200*ones(1,ngarcuts);
+ postseg_vec = 0*ones(1,ngarcuts);
+ sdeps = zeros(nrays,1);
+ end
+
+ fid = fopen([dir1 'ALL_rays'],'w');
+ aziall = zeros(nrays,2);
+ for ii = 1:nrays
+ azi1 = azimuth(rlats(ii),rlons(ii),slats(ii),slons(ii));
+ azi2 = azimuth(slats(ii),slons(ii),rlats(ii),rlons(ii));
+ aziall(ii,:) = [azi1 azi2];
+ if irun==1
+ fprintf(fid,'%3i%12i%10s%10.2f%10.2f%10.2f%10.2f%10.2f%10.2f\n',...
+ ii,eids(ii),recs{ii},slons(ii),slats(ii),rlons(ii),rlats(ii),aziall(ii,:));
+ else
+ fprintf(fid,'%3i%12i%10s%10.2f%10.2f%10.2f%10.2f%10.2f%10.2f\n',...
+ ii,NaN,'NAN',slons(ii),slats(ii),rlons(ii),rlats(ii),aziall(ii,:));
+ end
+ end
+ fclose(fid);
+
+ dx = 1; % horizontal sampling interval, km
+ zmax = 5*1e3; zmin = -45*1e3; znum = 100;
+
+% % update the modified file if you add any columns (i.e., faults)
+% if 0==1
+% dir1 = '/net/sierra/raid1/carltape/results/PLOTTING/SAVE/vert_01/';
+% [i1,i2,i3,i4,i5,f1,f2,f3,f4,f5,f6] = ...
+% textread([dir1 'ALL_rays_fault_positions_mod_old'],'%f%f%s%f%f%f%f%f%f%f%f');
+% [i1b,i2b,i3b,i4b,i5b,f1b,f2b,f3b,f4b,f5b,f6b,f7b,f8b,f9b] = ...
+% textread([dir1 'ALL_rays_fault_positions'],'%f%f%s%f%f%f%f%f%f%f%f%f%f%f');
+% fid = fopen([dir1 'ALL_rays_fault_positions_new'],'w');
+% stfmt = repmat('%10.0f',1,9);
+% for ii = 1:length(i1b)
+% fprintf(fid,['%3i%12i%10s%10.1f%10.1f' stfmt '\n'],...
+% i1b(ii),i2b(ii),i3b{ii},i4b(ii),i5b(ii),...
+% f1(ii),f2(ii),f3(ii),f4(ii),f5(ii),f6(ii),...
+% f7b(ii),f8b(ii),f9b(ii));
+% end
+% fclose(fid);
+% end
+
+ % load six faults
+ ffile = '/home/carltape/gmt/faults/socal_faults_xc.lonlat';
+ [data,ilabs,labels,finds] = read_delimited_file(ffile,2);
+ nfault = length(finds);
+ DMIN_KM = 10; % min distance must be DMIN_KM from the fault to mark the label
+
+ % compute distances to six selected faults
+ fid0 = fopen([dir1 'ALL_rays_fault_positions'],'w');
+
+ imin = 1; imax = nrays;
+ %imin = 80; imax = nrays;
+ %imin = nrays; imax = imin;
+
+ for ii = imin:imax
+
+ % source and receiver points
+ lon1 = slons(ii); lat1 = slats(ii);
+ lon2 = rlons(ii); lat2 = rlats(ii);
+
+ if irun==1
+ stit = sprintf('event %i to %s',eids(ii),recs{ii});
+ else
+ stit = 'model cut';
+ end
+ disp(sprintf('%s -- %i out of %i to %i',stit,ii,imin,imax));
+
+% switch ii
+% case 1
+% % this example has input points going outside the mesh:
+% % the bottom, top, and south side
+% znum = 100;
+% hnum = 200;
+% zmax = 5*1e3; zmin = -70*1e3;
+% lon1 = -119.6; lat1 = 35.6; lon2 = -115.9; lat2 = 31.2;
+% %rng1 = 600; az1 = 144;
+% case 2
+% znum = 100;
+% hnum = 100;
+% zmax = 5*1e3; zmin = -45*1e3;
+% lon1 = -119.6; lat1 = 34.6; lon2 = -117.4; lat2 = 33.4;
+% %rng1 = 220; az1 = 110;
+% end
+
+ % depth sampling
+ depvec = linspace(zmin,zmax,znum)';
+
+ % get points on the ray path
+ preseg = preseg_vec(ii);
+ postseg = postseg_vec(ii);
+ [lon_gc1,lat_gc1,rdist,azis] = raypath(lon1,lat1,lon2,lat2,preseg,postseg,dx); % KEY COMMAND
+ [lon_gc1_utm,lat_gc1_utm] = utm2llvec(lon_gc1,lat_gc1,11,0);
+
+ % exclude the points on the ray that are outside the mesh,
+ % then convert back to lat-lon
+ inds = getsubset(lon_gc1_utm,lat_gc1_utm,ax0_utm);
+ rdist = rdist(inds);
+ azis = azis(inds);
+ lon_gc1_utm = lon_gc1_utm(inds);
+ lat_gc1_utm = lat_gc1_utm(inds);
+ [lon_gc1,lat_gc1] = utm2llvec(lon_gc1_utm,lat_gc1_utm,11,1);
+ hnum = length(lon_gc1);
+
+ % compute distance to each of the six faults
+ dmins = 1e9*ones(nfault,1);
+ rinds = zeros(nfault,1);
+ for kk = 1:nfault
+ indf = [finds(kk,1):finds(kk,2)];
+ flon = data(indf,1); flat = data(indf,2); nf = length(flon);
+
+ dminall = 1e6*ones(nf,1);
+ iminall = zeros(nf,1);
+ for ij = 1:nf
+ [dtemp,itemp] = min(deg2km(distance(...
+ flon(ij)*ones(hnum,1),flat(ij)*ones(hnum,1),lon_gc1,lat_gc1)));
+ dminall(ij) = dtemp;
+ iminall(ij) = itemp;
+ end
+ [gmin,gind] = min(dminall);
+ rinds(kk) = iminall(gind);
+ dmins(kk) = dminall(gind);
+ end
+
+ rcuts = 9999*ones(1,nfault);
+ igood = find(dmins < DMIN_KM);
+ rcuts(igood) = rdist(rinds(igood));
+
+ % write the fault cuts to file
+ stfmt = repmat('%10.0f',1,nfault);
+ if irun==1
+ fprintf(fid0,['%3i%12i%10s%10.1f%10.1f' stfmt '\n'],...
+ ii,eids(ii),recs{ii},aziall(ii,:),rcuts);
+ else
+ fprintf(fid0,['%3i%12i%10s%10.1f%10.1f' stfmt '\n'],...
+ ii,NaN,'NaN',aziall(ii,:),rcuts);
+
+ end
+
+% % get points on the ray path
+% %[lat_gc1,lon_gc1] = track1(lat1,lon1,az1,km2deg(rng1),[],'degrees',hnum);
+% [lat_gc1,lon_gc1] = gcwaypts(lat1,lon1,lat2,lon2,hnum-1);
+% [lon_gc1_utm,lat_gc1_utm] = utm2llvec(lon_gc1,lat_gc1,11,0);
+%
+% % compute the distance in km from one point to the end --
+% % this is used in the 2D plots
+% [rng,az] = distance(lat1*ones(hnum,1),lon1*ones(hnum,1),lat_gc1,lon_gc1);
+% rdist = deg2km(rng);
+
+ figure; hold on;
+ plot(xbox_ll, ybox_ll, '.');
+ plot(lon_gc1,lat_gc1,'b.');
+ plot(lon_gc1(rinds),lat_gc1(rinds),'ko','markersize',10);
+ for kk=1:nfault, indf = [finds(kk,1):finds(kk,2)]; plot(data(indf,1),data(indf,2)); end
+ text(lon_gc1(1),lat_gc1(1),num2str(1),'fontsize',16);
+ text(lon_gc1(hnum),lat_gc1(hnum),num2str(hnum),'fontsize',16);
+ plot(lon1,lat1,'rp',lon2,lat2,'ko','markersize',20);
+ title(stit);
+
+% % reference points at the surface
+% % COLUMNS: dist from start (km), utm-x, utm-y, lon, lat
+% filename = sprintf('socal_xc_%3.3i_ref',ii);
+% fid = fopen([pdir filename],'w');
+% for kk=1:hnum % loop over horizontal points in the ray path
+% fprintf(fid,'%16.6e%16.6e%16.6e%16.6e%16.6e\n',rdist(kk),...
+% lon_gc1_utm(kk),lat_gc1_utm(kk),lon_gc1(kk),lat_gc1(kk));
+% end
+% fclose(fid);
+
+ if iwrite == 1
+ % xyz files
+ filename = sprintf('vert_xc_%3.3i_input',ii);
+ fid = fopen([dir1 filename],'w');
+ for kk=1:hnum % loop over horizontal points in the ray path
+ for jj=1:znum
+ fprintf(fid,'%16.6e%16.6e%16.6e%16.6e\n',...
+ lon_gc1_utm(kk),lat_gc1_utm(kk),depvec(jj),1000*rdist(kk));
+ end
+ end
+ fclose(fid);
+
+ % ray path
+ filename = sprintf('vert_xc_%3.3i_ray_path',ii);
+ fid = fopen([dir1 filename],'w');
+ for kk=1:hnum % loop over horizontal points in the ray path
+ fprintf(fid,'%16.6e%16.6e%16.6e%16.6e%16.6e\n',...
+ lon_gc1(kk),lat_gc1(kk),lon_gc1_utm(kk),lat_gc1_utm(kk),1000*rdist(kk));
+ end
+ fclose(fid);
+
+ % lon-lat points of source and receiver
+ filename = sprintf('vert_xc_%3.3i_A_B',ii);
+ fid = fopen([dir1 filename],'w');
+ [dq,azq] = distance(slats(ii),slons(ii),rlats(ii),rlons(ii)); % distance to receiver
+ iflip = 0;
+ %if azq > 180, iflip = 1; end
+ if rlons(ii) < slons(ii), iflip = 1; end
+ fprintf(fid,[repmat('%16.6e',1,9) '%4i\n'],...
+ slons(ii),slats(ii),0,sdeps(ii),rlons(ii),rlats(ii),deg2km(dq),0,azq,iflip);
+ fclose(fid);
+ end
+
+ end % for ii=1:nrays
+ fclose(fid0); % fault cut positions
+
+end
+
+%----------------------------
+% PLOTTING CROSS-SECTIONS
+
+if iplot_vert==1
+ colors;
+
+ %irun = 1; pmax = 86;
+ %irun = 2; pmax = 32;
+ irun = 3; pmax = 10;
+
+ dir1 = [pdir 'SAVE/vert_' sprintf('%2.2i',irun) '/'];
+ smodel1 = 'm00';
+ smodel2 = 'm16';
+ modlab1 = 'vs';
+ modlab2 = 'vb';
+ smodeltag = [smodel2 '_' smodel1];
+
+ pmin = 1;
+ %pmin = 80; pmax = 60;
+ %pmin = pmax; pmax = pmin;
+
+ for ip = pmin:pmax
+ stip = sprintf('%3.3i',ip);
+ disp(sprintf('model %i out of %i',ip,pmax));
+
+ % four possible combinations
+ tag3 = [modlab1 '_' smodel1];
+ tag4 = [modlab1 '_' smodel2];
+ tag5 = [modlab2 '_' smodel1];
+ tag6 = [modlab2 '_' smodel2];
+
+ % load files needed for depth cross-sections
+ file0 = [dir1 'vert_xc_' stip '_A_B'];
+ file1 = [dir1 'vert_xc_' stip '_ray_path'];
+ file2 = [dir1 'vert_xc_' stip '_input'];
+ file3 = [dir1 'vert_xc_' tag3 '_' stip '.gmt'];
+ file4 = [dir1 'vert_xc_' tag4 '_' stip '.gmt'];
+ file5 = [dir1 'vert_xc_' tag5 '_' stip '.gmt'];
+ file6 = [dir1 'vert_xc_' tag6 '_' stip '.gmt'];
+
+ [lon1,lat1,d1,z1,lon2,lat2,d2,z2,az,iflip] = textread(file0,'%f%f%f%f%f%f%f%f%f%f');
+ [lon_gc1,lat_gc1,lon_gc1_utm,lat_gc1_utm,rdist] = textread(file1,'%f%f%f%f%f');
+ [utm_x,utm_y,z,d] = textread(file2,'%f%f%f%f');
+ [junk1,junk2,junk3,vs1,cdist] = textread(file3,'%f%f%f%f%f');
+ [junk1,junk2,junk3,vs2,cdist] = textread(file4,'%f%f%f%f%f');
+ [junk1,junk2,junk3,vb1,cdist] = textread(file5,'%f%f%f%f%f');
+ [junk1,junk2,junk3,vb2,cdist] = textread(file6,'%f%f%f%f%f');
+ n = length(vs1);
+
+ %lon1 = lon_gc1(1); lat1 = lat_gc1(1);
+ %lon2 = lon_gc1(end); lat2 = lat_gc1(end);
+
+ % turn negative points into NaN values
+ ineg = find(vs1 < 0); vs1(ineg) = NaN;
+ ineg = find(vs2 < 0); vs2(ineg) = NaN;
+ ineg = find(vb1 < 0); vb1(ineg) = NaN;
+ ineg = find(vb2 < 0); vb2(ineg) = NaN;
+
+ % compute Vp
+ vp1 = zeros(n,1); vp2 = zeros(n,1);
+ vp1 = sqrt( 4/3 * vs1.^2 + vb1.^2 );
+ vp2 = sqrt( 4/3 * vs2.^2 + vb2.^2 );
+
+ % compute the difference -- GMT does not seem to like operating on
+ % NaN vectors, so we have to do it here and then write it out
+ vs3 = log(vs2 ./ vs1);
+ vb3 = log(vb2 ./ vb1);
+ vp3 = log(vp2 ./ vp1);
+
+ if ifigure == 1
+ for kk = 1:3
+ switch kk
+ case 1, cval = vs1;
+ case 2, cval = vs2;
+ case 3, cval = vs3;
+ end
+
+ % construct mesh for Matlab plotting
+ xvec = d*1e-3;
+ yvec = z*1e-3;
+ zvec = cval;
+ xlin = linspace(min(xvec), max(xvec), 200);
+ ylin = linspace(min(yvec), max(yvec), 100);
+ [X,Y] = meshgrid(xlin,ylin);
+
+ % determine interpolated function using xvec,yvec input
+ %whos xvec yvec zvec
+ Z = griddata(xvec, yvec, zvec, X, Y, 'nearest');
+
+ figure; nr=3; nc=1;
+ %clims = [2000 8500];
+
+ subplot(nr,nc,1); hold on;
+ plot(xbox_ll, ybox_ll, '.');
+ plot(lon_gc1,lat_gc1,'b.');
+ plot(lon1,lat1,'rp','markersize',20);
+ axis equal, axis tight
+
+ subplot(nr,nc,2);
+ pcolor(X,Y,Z); shading flat;
+ xlabel('Distance along cross-section, km');
+ ylabel('Depth, km');
+ %caxis(clims);
+ colorbar;
+
+ subplot(nr,nc,3);
+ pcolor(X,Y,Z); shading flat;
+ xlabel('Distance along cross-section, km');
+ ylabel('Depth, km');
+ %caxis(clims);
+ colorbar; axis equal, axis tight
+
+ colormap(seis);
+ end
+ end
+
+ if iwrite == 1
+% % write out a file with five columns: xdist, zdep, m00 (m/s), m16 (m/s), ln(m16/m00)
+% fid = fopen([dir1 'vert_xc_' modlab '_' smodeltag '_' stip '.dat'],'w');
+% for kk=1:length(d)
+% fprintf(fid,'%16.6e%16.6e%16.6e%16.6e%16.6e\n',...
+% d(kk)*1e-3,z(kk)*1e-3,cval1(kk)*1e-3,cval2(kk)*1e-3,cval3(kk));
+% end
+% fclose(fid);
+
+ % write out a file with 11 columns: xdist, zdep, m00 (m/s), m16 (m/s), ln(m16/m00)
+ fid = fopen([dir1 'vert_xc_vsvbvp_' smodeltag '_' stip '.dat'],'w');
+ for kk=1:length(d)
+ fprintf(fid,[repmat('%16.6e',1,11) '\n'],...
+ d(kk)*1e-3,z(kk)*1e-3,vs1(kk)*1e-3,vs2(kk)*1e-3,vs3(kk),...
+ vb1(kk)*1e-3,vb2(kk)*1e-3,vb3(kk),...
+ vp1(kk)*1e-3,vp2(kk)*1e-3,vp3(kk) );
+ end
+ fclose(fid);
+ %whos vs1 vs2 vs3 vb1 vb2 vb3 vp1 vp2 vp3, length(d)
+ end
+
+ end
+
+end
+
+%----------------------------
+% PLOTTING CROSS-SECTIONS
+
+if iplot_horz==1
+ colors;
+
+ irun = 1;
+ dir1 = [pdir 'SAVE/horz_' sprintf('%2.2i',irun) '/'];
+ %dir1 = [pdir 'SAVE/horz_full/'];
+ smodel1 = 'm00';
+ smodel2 = 'm16'; % m16, m01, etc
+ modlab = 'vs';
+ smodeltag = [smodel2 '_' smodel1];
+
+ % KEY: TWO DIFFERENT KINDS OF SUMMED KERNELS
+ % load number of window picks per event
+ nwins_all = load([mdir 'm16/window_tomo/m16_window_picks_nwin']);
+ nevent = length(nwins_all);
+ nwin = sum(nwins_all);
+ cfac = log10(1e-16);
+ %covertag = 'coverage_sum_abs'; Nfac = nwin;
+ covertag = 'coverage_sum_abs_nwin'; Nfac = nevent;
+ cnorm = 10^cfac;
+
+ % load files needed for horizontal cross-sections
+ [inds,depvec] = textread([dir1 'horz_xc_elevation_cuts'],'%f%f');
+ [xlon,ylat,xutm,yutm] = textread([dir1 'horz_xc_lonlat'],'%f%f%f%f');
+
+ % replace values far west in regions of no coverage with NaN
+ %xvpts = [-119 -122 -122 -119];
+ %yvpts = [32 36 32 32];
+ %imask = inpolygon(xlon,ylat,xvpts,yvpts);
+
+ %figure; hold on;
+ %plot(xlon,ylat,'b.');
+ %plot(xlon(imask),ylat(imask),'ro');
+
+ pmin = 1; pmax = length(depvec);
+ %pmin = 21; pmax = pmin;
+ %pmin = 3; pmax = length(depvec);
+
+ for ip = pmin:pmax
+ stip = sprintf('%3.3i',ip);
+ disp(sprintf('%i out of the set %i - %i',ip,pmin,pmax));
+
+ % load output from cluster code
+ tag1 = [modlab '_' smodel1];
+ tag2 = [modlab '_' smodel2];
+ file1 = [dir1 'horz_xc_elevation_' tag1 '_' stip '.gmt'];
+ [junk1,junk2,junk3,cval1,cdist] = textread(file1,'%f%f%f%f%f');
+ file2 = [dir1 'horz_xc_elevation_' tag2 '_' stip '.gmt'];
+ [junk1,junk2,junk3,cval2,cdist] = textread(file2,'%f%f%f%f%f');
+
+ % load coverage mask
+ tagm = [modlab '_' covertag];
+ %stip0=stip;
+ %if ip==3, stip0 = '002'; end
+ filem = [dir1 'horz_xc_elevation_' tagm '_' stip '.gmt'];
+ [junk1,junk2,junk3,cvalm,cdist] = textread(filem,'%f%f%f%f%f');
+
+ % for some reason matlab interpolation needs to work with smaller numbers,
+ % so we subtract a reference value
+ xutm0 = min(xutm);
+ yutm0 = min(yutm);
+ xutm = xutm - xutm0;
+ yutm = yutm - yutm0;
+
+ %cnorm = 3000; % plot values w.r.t this one
+
+ % replace values far west in regions of no coverage with NaN
+ %cval1(imask) = NaN;
+ %cval2(imask) = NaN;
+
+ % remove points above the topography
+ ineg = find(cval1 < 0); cval1(ineg) = NaN;
+ ineg = find(cval2 < 0); cval2(ineg) = NaN;
+
+ % compute the difference -- GMT does not seem to like operating on
+ % NaN vectors, so we have to do it here and then write it out
+ cval3 = log(cval2 ./ cval1);
+
+ %----------------------------------------------------
+ % COVERAGE MASK
+
+ % take log10 of the field
+ cvalmlog = log10(cvalm / Nfac);
+
+ % replace points outside coverage with a mask
+ imask = find(cvalmlog <= cfac);
+ cval1(imask) = NaN;
+ cval2(imask) = NaN;
+ cval3(imask) = NaN;
+
+ %----------------------------------------------------
+
+ clabs = {smodel1,smodel2,['ln(' smodel2 '/' smodel1 ')']};
+ for kk = 1:3
+ %for kk = 2:2
+ switch kk
+ case 1, cval = cval1;
+ case 2, cval = cval2;
+ case 3, cval = cval3;
+ end
+
+ if ifigure ==1
+ % replace negative values with NaN
+ %ineg = find(cval < 0);
+ %cval(ineg) = NaN;
+
+ % construct mesh for plotting
+ xlin = linspace(min(xutm), max(xutm), 200);
+ ylin = linspace(min(yutm), max(yutm), 200);
+ [X,Y] = meshgrid(xlin,ylin);
+
+ % determine interpolated function using xvec,yvec input
+ Z = griddata(xutm, yutm, cval, X, Y, 'nearest');
+
+ figure;
+ clims = [-1 1]*0.2;
+ pcolor((X+xutm0)*1e-3,(Y+yutm0)*1e-3,Z);
+ shading flat;
+ %caxis(clims);
+ xlabel('X distance, km');
+ ylabel('Y distance, km');
+ colorbar; axis equal, axis tight
+ title(sprintf('Depth = %.1f km -- %s',depvec(ip)*1e-3,clabs{kk}));
+ %title(sprintf('Depth = %.1f km (Vref = %.1f km/s)',depvec(ip)*1e-3,cnorm*1e-3));
+ colormap(seis);
+ end
+ end
+
+ if iwrite == 1
+ % write out a file with five columns: lon, lat, m00 (m/s), m16 (m/s), ln(m16/m00)
+ fid = fopen([dir1 '/horz_xc_' modlab '_' smodeltag '_' stip '.dat'],'w');
+ for kk=1:length(xlon)
+ fprintf(fid,'%16.6e%16.6e%16.6e%16.6e%16.6e\n',...
+ xlon(kk),ylat(kk),cval1(kk),cval2(kk),cval3(kk));
+ end
+ fclose(fid);
+ end
+ end
+
+end
+
+%----------------------------
+% PLOTTING COVERAGE CROSS-SECTIONS
+
+if iplot_horz_coverage==1
+ colors;
+
+ irun = 1;
+ dir1 = [pdir 'SAVE/horz_' sprintf('%2.2i',irun) '/'];
+ %dir1 = [pdir 'SAVE/horz_full/'];
+ modlab = 'vb';
+ %smodeltag = [smodel2 '_' smodel1];
+
+ % load number of window picks per event
+ nwins_all = load([mdir 'm16/window_tomo/m16_window_picks_nwin']);
+ nevent = length(nwins_all);
+ nwin = sum(nwins_all);
+
+ % KEY: TWO DIFFERENT KINDS OF SUMMED KERNELS
+ cfac = log10(1e-16);
+ smodeltag = 'coverage_sum_abs'; Nfac = nwin;
+ %smodeltag = 'coverage_sum_abs_nwin'; Nfac = nevent;
+ cnorm = 10^cfac;
+
+ % load files needed for depth cross-sections
+ [inds,depvec] = textread([dir1 'horz_xc_elevation_cuts'],'%f%f');
+ [xlon,ylat,xutm,yutm] = textread([dir1 'horz_xc_lonlat'],'%f%f%f%f');
+
+ % replace values far west in regions of no coverage with NaN
+ %xvpts = [-119 -122 -122 -119];
+ %yvpts = [32 36 32 32];
+ %imask = inpolygon(xlon,ylat,xvpts,yvpts);
+
+ %figure; hold on;
+ %plot(xlon,ylat,'b.');
+ %plot(xlon(imask),ylat(imask),'ro');
+
+ pmin = 1; pmax = length(depvec);
+ %pmin = 21; pmax = pmin;
+
+ for ip = pmin:pmax
+ stip = sprintf('%3.3i',ip);
+ disp(sprintf('%i out of %i - %i',ip,pmin,pmax));
+
+ % load output from cluster code
+ tag1 = [modlab '_' smodeltag];
+ file1 = [dir1 'horz_xc_elevation_' tag1 '_' stip '.gmt'];
+ [junk1,junk2,junk3,cval1,cdist] = textread(file1,'%f%f%f%f%f');
+
+ % for some reason matlab interpolation needs to work with smaller numbers,
+ % so we subtract a reference value
+ xutm0 = min(xutm);
+ yutm0 = min(yutm);
+ xutm = xutm - xutm0;
+ yutm = yutm - yutm0;
+
+ % remove points above the topography
+ ineg = find(cval1 < 0); cval1(ineg) = NaN;
+
+ % normalize by number of windows or number of events
+ cval1 = cval1 / Nfac;
+
+ % take -log10 of the field
+ cval2 = log10(cval1);
+
+ % replace points outside coverage with a mask
+ imask = find(cval2 <= cfac);
+ cval3 = cval2; cval3(imask) = NaN;
+ cval1mask = cval1; cval1mask(imask) = NaN;
+
+ % convex hull -- convhull
+
+ clabs = {tag1,['log10(' tag1 ')'],['log10(' tag1 ')']};
+ for kk = 1:3
+ %for kk = 2:2
+ switch kk
+ case 1, cval = cval1;
+ case 2, cval = cval2;
+ case 3, cval = cval3;
+ end
+
+ if ifigure ==1
+ % replace negative values with NaN
+ %ineg = find(cval < 0);
+ %cval(ineg) = NaN;
+
+ % construct mesh for plotting
+ xlin = linspace(min(xutm), max(xutm), 200);
+ ylin = linspace(min(yutm), max(yutm), 200);
+ [X,Y] = meshgrid(xlin,ylin);
+
+ % determine interpolated function using xvec,yvec input
+ Z = griddata(xutm, yutm, cval, X, Y, 'nearest');
+
+ figure;
+ clims = [-1 1]*0.2;
+ pcolor((X+xutm0)*1e-3,(Y+yutm0)*1e-3,Z);
+ shading flat;
+ %caxis(clims);
+ xlabel('X distance, km');
+ ylabel('Y distance, km');
+ colorbar; axis equal, axis tight
+ title(sprintf('Depth = %.1f km -- %s',depvec(ip)*1e-3,clabs{kk}));
+ %title(sprintf('Depth = %.1f km (Vref = %.1f km/s)',depvec(ip)*1e-3,cnorm*1e-3));
+ colormap(seis);
+ end
+ end
+
+ if iwrite == 1
+ % write out a file with five columns: lon, lat, log10(Ksum), log10(Ksum)+mask
+ filetag = ['horz_xc_' modlab '_' smodeltag '_' stip];
+ fid = fopen([dir1 filetag '.dat'],'w');
+ for kk=1:length(xlon)
+ fprintf(fid,'%16.6e%16.6e%16.6e%16.6e\n',...
+ xlon(kk),ylat(kk),cval1(kk),cval1mask(kk));
+ end
+ fclose(fid);
+
+ % threshold value for mask
+ fid = fopen([dir1 filetag '_mask_value.dat'],'w');
+ fprintf(fid,'%.6e\n',cnorm);
+ fclose(fid);
+ end
+ end
+
+end
+
+%====================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,171 @@
+%
+% function [xo,yo,i_zone] = utm2ll(xi,yi,i_zone,i_type)
+% Carl Tape, 22-April-2004
+% printed xxx
+%
+% This function converts between lat-lon and UTM coordinates.
+% Note: input is NOT in vector form (one value at a time).
+%
+% FUNCTION INPUT:
+% xi longitude or x-UTM
+% yi latitude or y-UTM
+% i_type lat-lon --> UTM (=0) or UTM --> lat-lon (=1)
+% i_zone UTM zone (if i_type=1)
+%
+%
+%****************************************************************
+%**
+%** FILE NAME: utmtoll
+%**
+%** DATE WRITTEN: 7/22/93
+%**
+%** PROGRAMMER: Scott Hensley
+%**
+%** FUNCTIONAL DESCRIPTION: This routine converts between lat
+%** lon and utm coordinates for a datum determined from the input
+%** a and e2.
+%**
+%** ROUTINES CALLED:none
+%**
+%** NOTES: none
+%**
+%** UPDATE LOG:
+%** MATLAB version, Yuri Fialko 11/25/99
+%****************************************************************
+%
+% INPUT VARIABLES:
+% real*8 xi,yi !input coordinates
+% integer i_zone !UTM zone (needed if i_type=2)
+% integer i_type ! [1 = lat,lon to utm], [2 = utm to lat,lon]
+% old stuff:
+% real*8 r_a !ellispoid semi-major axis
+% real*8 r_e2 !ellipsoid eccentricity squared
+% real*8 r_v(2) !geocentric vector (meters)
+% real*8 r_lat !latitude (deg -90 to 90)
+% real*8 r_lon !longitude (deg -180 to 180)
+% character*1 a_grid !UTM North-South grid
+%
+% OUTPUT VARIABLES:
+% real*8 xo,yo !output coordinates
+%
+% LOCAL VARIABLES:
+% integer i_ft,i_gi
+% real*8 pi,r_dtor
+% real*8 r_ep2,r_k0,r_k
+% real*8 r_fe,r_fn(2)
+% real*8 r_e4,r_e6,r_n,r_t,r_t2,r_c,r_c2,r_ba
+% real*8 r_a2,r_a3,r_a4,r_a5,r_a6
+% real*8 r_d,r_d2,r_d3,r_d4,r_d5,r_d6
+% real*8 r_lon0,r_lat1,r_m,r_m0,r_mu,r_lat0
+% real*8 r_et,r_e1,r_e12,r_e13,r_e14,r_r,r_vu(2)
+% character*1 a_griddes(20)
+%
+%****************************************************************
+
+function [xo,yo,i_zone] = utm2ll(xi,yi,i_zone,i_type)
+
+if and(i_type ~= 1, i_type ~= 0), error('i_type must be 1 or 0'); end
+
+% DATA STATEMENTS
+
+r_dtor=1.74532925199d-2;
+i_ft=0;
+a_griddes=['C','D','E','F','G','H','J',...
+'K','L','M','N','P','Q','R','S','T','U',...
+'V','W','X'];
+r_a=6378206.4d0;
+r_e2=0.00676865799761d0;
+r_k0=0.9996d0; %scale at center
+r_lat0=0.0d0;
+r_fe=5e5;
+r_fn(1)=0;
+r_fn(2)=1e7; % N-S hemispheres??
+
+% PROCESSING STEPS
+
+r_ep2 = r_e2/(1.d0 - r_e2);
+r_e4 = r_e2^2;
+r_e6 = r_e2^3;
+r_dtor = pi/180;
+
+if i_type == 0 % convert lat,lon to UTM
+
+ xi=xi*r_dtor;
+ yi=yi*r_dtor;
+ if (i_zone==0)
+ i_zone = fix(mod(xi+3.d0*pi,2.d0*pi)/(r_dtor*6.d0)) + 1;
+ i_zone = max(min(i_zone,60),1);
+ end
+ r_lon0 = -pi + 6.d0*r_dtor*(i_zone-1) + 3.d0*r_dtor; % central meridian
+
+ r_n = r_a/sqrt(1.d0 - r_e2*sin(yi)^2);
+ r_t = tan(yi)^2;
+ r_t2 = r_t^2;
+ r_c = r_ep2*cos(yi)^2;
+ r_ba = (xi - r_lon0)*cos(yi);
+ r_a2 = r_ba^2;
+ r_a3 = r_ba*r_a2;
+ r_a4 = r_ba*r_a3;
+ r_a5 = r_ba*r_a4;
+ r_a6 = r_ba*r_a5;
+ r_m = r_a*((1.0d0-r_e2/4 - 3.0d0*r_e4/64.0d0-5.d0*r_e6/256.d0)*yi - (3.d0*r_e2/8.d0 + 3.d0*r_e4/32.d0 +45.d0*r_e6/1024.d0)*sin(2.d0*yi) + (15.d0*r_e4/256.d0 +45.d0*r_e6/1024.d0)*sin(4.d0*yi) - (35.d0*r_e6/3072.d0)*sin(6.d0*yi));
+ r_m0 = r_a*((1.d0-r_e2/4 - 3.d0*r_e4/64.d0 -5.d0*r_e6/256.d0)*r_lat0 - (3.d0*r_e2/8.d0 + 3.d0*r_e4/32.d0 +45.d0*r_e6/1024.d0)*sin(2.d0*r_lat0) + (15.d0*r_e4/256.d0 +45.d0*r_e6/1024.d0)*sin(4.d0*r_lat0) - (35.d0*r_e6/3072.d0)*sin(6.d0*r_lat0));
+
+ r_v(1) = r_k0*r_n*(r_ba+(1.d0-r_t+r_c)*r_a3/6.d0 +(5.d0-18.d0*r_t+r_t2+72.d0*r_c-58.d0*r_ep2)*r_a5/120.d0);
+ r_v(1) = r_v(1) + r_fe;
+
+ r_v(2) = r_k0*(r_m - r_m0 + r_n*tan(yi)*( r_a2/2.d0 + (5.d0-r_t+9.d0*r_c+4.d0*r_c^2)*(r_a4/24.d0) + (61.d0-58.d0*r_t+r_t2+600.d0*r_c-330.d0*r_ep2)*(r_a6/720.d0) ));
+ if yi >= 0
+ r_v(2) = r_v(2) + r_fn(1);
+ else
+ r_v(2) = r_v(2) + r_fn(2);
+ end
+
+ r_k = r_k0*(1.d0+(1.d0+r_ep2*cos(yi)^2)*(r_v(1)-r_fe)^2/(2.d0*(r_k0^2)*r_n^2));
+
+ i_gi = fix((yi/r_dtor+80.d0)/8.d0) + 1;
+ i_gi = max(min(i_gi,20),1);
+ a_grid = a_griddes(i_gi);
+
+ xo=r_v(1);
+ yo=r_v(2);
+
+elseif i_type == 1 % convert UTM to lat,lon
+
+ r_vu(1) = xi - r_fe;
+ r_vu(2) = yi;
+ if r_vu(2) >= r_fn(2) % LOOK at this if in the S hemisphere!!
+ r_vu(2) = yi - r_fn(2);
+ end
+ r_lon0 = -pi + 6.d0*r_dtor*(i_zone-1) + 3.d0*r_dtor;
+
+ r_et = sqrt(1.d0-r_e2);
+ r_e1 = (1.d0-r_et)/(1.d0+r_et);
+ r_e12 = r_e1^2;
+ r_e13 = r_e1*r_e12;
+ r_e14 = r_e1*r_e13;
+ r_m = r_vu(2)/r_k0;
+ r_mu = r_m/(r_a*(1.d0-r_e2/4.d0-3.d0*r_e4/64.d0-5.d0*r_e6/256.d0));
+ r_lat1 = r_mu + (3.d0*r_e1/2.d0-27.d0*r_e13/32.d0)*sin(2.d0*r_mu)+(21.d0*r_e12/16.d0-55.d0*r_e14/32.d0)*sin(4.d0*r_mu) +(51.d0*r_e13/96.d0)*sin(6.d0*r_mu) +(1097.d0*r_e14/512.d0)*sin(8.d0*r_mu);
+
+ r_n = r_a/sqrt(1.d0 - r_e2*sin(r_lat1)^2);
+ r_r = (r_a*(1.d0-r_e2))/sqrt(1.d0 - r_e2*sin(r_lat1)^2)^3;
+ r_t = tan(r_lat1)^2;
+ r_t2 = r_t^2;
+ r_c = r_ep2*cos(r_lat1)^2;
+ r_c2 = r_c^2;
+ r_d = r_vu(1)/(r_n*r_k0);
+ r_d2 = r_d^2;
+ r_d3 = r_d2*r_d;
+ r_d4 = r_d3*r_d;
+ r_d5 = r_d4*r_d;
+ r_d6 = r_d5*r_d;
+
+ yo = r_lat1 - (r_n*tan(r_lat1)/r_r)*(r_d2/2.d0+(5.d0+3.d0*r_t+10.d0*r_c-4.d0*r_c2-9.d0*r_ep2)*r_d4/24.d0 +(61.d0+90*r_t+298.d0*r_c+45.d0*r_t2-252.d0*r_ep2-3.d0*r_c2)*(r_d6/720.d0));
+ xo = r_lon0 + (r_d - (1.d0+2.d0*r_t+r_c)*r_d3/6.d0 +(5.d0-2.d0*r_c+28.d0*r_t-3.d0*r_c2+8.d0*r_ep2+24.d0*r_t2)*(r_d5/120.d0))/cos(r_lat1);
+ xo=xo/r_dtor;
+ yo=yo/r_dtor;
+
+end
+
+%==========================================================================
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m 2009-03-28 01:07:39 UTC (rev 14507)
@@ -0,0 +1,29 @@
+%
+% function [x,y] = utm2llvec(xi,yi,i_zone,i_type)
+% Carl Tape, 10-Nov-2004
+% printed xxx
+%
+%
+% calls utm2ll.m
+% called by omsplotsC.m
+%
+
+function [x,y] = utm2llvec(xi,yi,i_zone,i_type)
+
+[a,b] = size(xi);
+num = a*b;
+xvec0 = reshape(xi,num,1);
+yvec0 = reshape(yi,num,1);
+xvec = zeros(num,1);
+yvec = zeros(num,1);
+
+for ii=1:length(xi)
+ %if mod(length(xi),1000)==0
+ % disp(sprintf('UTM to LATLON : %i out of %i',ii,length(xi)));
+ %end
+ [xvec(ii),yvec(ii)] = utm2ll(xvec0(ii),yvec0(ii),i_zone,i_type);
+end
+x = reshape(xvec,a,b);
+y = reshape(yvec,a,b);
+
+%==========================================================================
More information about the CIG-COMMITS
mailing list