[cig-commits] r14599 - seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab
carltape at geodynamics.org
carltape at geodynamics.org
Sun Apr 5 16:09:44 PDT 2009
Author: carltape
Date: 2009-04-05 16:09:44 -0700 (Sun, 05 Apr 2009)
New Revision: 14599
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m
Log:
Updated Matlab script for tomography model cross sections.
Modified: 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 2009-04-05 23:08:26 UTC (rev 14598)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m 2009-04-05 23:09:44 UTC (rev 14599)
@@ -17,9 +17,9 @@
format short, format compact
% USER OUTPUT DIRECTORIES
-dir0 = '/net/denali/raid1/carltape/svn/cig/seismo/3D/ADJOINT_TOMO/ADJOINT_TOMO_OUTPUT/';
+dir0 = '/home/carltape/ADJOINT_TOMO/ADJOINT_TOMO_OUTPUT/';
pdir = [dir0 'model_plot_matlab_OUTPUT/'];
-mdir = [dir0 'OUTPUT_SUBSPACE_m16/'];
+mdir = [dir0 'OUTPUT_subspace_save_m16/'];
% KEY: vertical or horizontal cross-sections
ivert = 1;
@@ -29,7 +29,7 @@
iplot_horz = 0;
iplot_horz_coverage = 0;
-iwrite = 0;
+iwrite = 1;
ifigure = 0;
%-----------------------------------
@@ -55,20 +55,32 @@
if ihorz == 1
- irun = 1;
+ irun = 2; % 1 for 0 km, 1 km, 2 km, etc
+ % 2 for 0.5 km, 1.5 km, etc
dir1 = [pdir 'horz_' sprintf('%2.2i',irun) '/'];
-
- % 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);
+
+ switch irun
+ case 1
+ depvec = 1e3 * [ 0:-1:-40]; % KEY
+ numx = 300; % key command for resolution
+ xvec = linspace(xmin,xmax,numx);
+ dx = xvec(2) - xvec(1);
+ yvec = [ymin:dx:ymax];
+ numy = length(yvec);
+ case 2
+ depvec = 1e3 * [ 0:-0.5:-40]; % KEY
+ dx = 2000; % key command for resolution
+ xvec = [xmin:dx:xmax];
+ yvec = [ymin:dx:ymax];
+ numx = length(xvec);
+ numy = length(yvec);
+ end
+ ndep = length(depvec);
[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);
@@ -77,23 +89,18 @@
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');
+ fid = fopen([dir1 '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');
+ fid = fopen([dir1 '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
@@ -103,7 +110,7 @@
for ii=1:ndep
zdep = depvec(ii)
filename = sprintf('horz_xc_elevation_%3.3i',ii);
- fid = fopen([pdir filename],'w');
+ fid = fopen([dir1 filename],'w');
for kk=1:N
fprintf(fid,'%16.6e%16.6e%16.6e\n',Xvec(kk),Yvec(kk),zdep);
end
@@ -121,10 +128,16 @@
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';
+ srcfile = '/home/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
+ % irun
+ % 1 src-rec
+ % 2 SAF + 3
+ % 3 Garlock
+ % 4 LON
+ % 5 LAT
+ irun = 5; % KEY: which set of cross-sections to make
dir1 = [pdir 'vert_' sprintf('%2.2i',irun) '/']; % KEY COMMAND
switch irun
@@ -175,7 +188,7 @@
preseg_vec = 100*ones(nrays,1);
postseg_vec = 100*ones(nrays,1);
- case 2 % 2: PICK SOURCE-RECEIVER PAIRS
+ case 2 % 2: San Andreas
% select profiles perpendicular to the SAF -- saf.m
saflonlat_cuts = [
@@ -219,7 +232,7 @@
postseg_vec = [100 150 150 0*ones(1,nsafcuts)];
sdeps = zeros(nrays,1);
- case 3 % 2: PICK SOURCE-RECEIVER PAIRS
+ case 3 % 3: Garlock
% select profiles perpendicular to the Garlock -- saf.m
garlonlat_cuts = [
@@ -243,6 +256,56 @@
preseg_vec = 200*ones(1,ngarcuts);
postseg_vec = 0*ones(1,ngarcuts);
sdeps = zeros(nrays,1);
+
+ case 4
+ % cuts of constant longitude
+ lon_cuts = [
+ -115.0860 32.2367 -115.0860 37.1830
+ -115.5080 32.9593 -115.5080 37.9056
+ -115.9960 33.5937 -115.9960 38.5400
+ -116.4850 33.8932 -116.4850 38.8395
+ -116.9810 34.0651 -116.9810 39.0114
+ -117.5030 34.2957 -117.5030 39.2420
+ -118.0110 34.5118 -118.0110 39.4581
+ -118.5090 34.6976 -118.5090 39.6439
+ -118.9980 34.8209 -118.9980 39.7672
+ -119.4850 34.9868 -119.4850 39.9331
+ -119.9980 35.4483 -119.9980 40.3946
+ -120.4920 35.9436 -120.4920 40.8899
+ -120.9890 36.4062 -120.9890 41.3525 ];
+ nloncuts = length(lon_cuts);
+
+ slons = lon_cuts(:,1)';
+ slats = lon_cuts(:,2)';
+ rlons = lon_cuts(:,3)';
+ rlats = lon_cuts(:,4)';
+ nrays = length(slons);
+ preseg_vec = 250*ones(1,nloncuts);
+ postseg_vec = 0*ones(1,nloncuts);
+ sdeps = zeros(nrays,1);
+
+ case 5
+ % cuts of constant latitude
+ lat_cuts = [
+ -115.2100 32.5132 -114.7000 32.5132
+ -115.5940 33.0397 -114.7000 33.0397
+ -115.9090 33.4981 -114.7000 33.4981
+ -116.8660 34.0004 -114.7000 34.0004
+ -117.9800 34.5007 -114.7000 34.5007
+ -119.4850 34.9868 -114.7000 34.9868
+ -120.0410 35.4953 -114.7000 35.4953
+ -120.5640 36.0081 -114.7000 36.0081
+ -121.0880 36.5124 -114.7000 36.5124 ];
+ nlatcuts = length(lat_cuts);
+
+ slons = lat_cuts(:,1)';
+ slats = lat_cuts(:,2)';
+ rlons = lat_cuts(:,3)';
+ rlats = lat_cuts(:,4)';
+ nrays = length(slons);
+ preseg_vec = 300*ones(1,nlatcuts);
+ postseg_vec = 0*ones(1,nlatcuts);
+ sdeps = zeros(nrays,1);
end
fid = fopen([dir1 'ALL_rays'],'w');
@@ -266,7 +329,7 @@
% % 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/';
+% dir1 = '/home/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] = ...
@@ -293,7 +356,7 @@
imin = 1; imax = nrays;
%imin = 80; imax = nrays;
- imin = nrays; imax = imin;
+ %imin = nrays; imax = imin;
for ii = imin:imax
@@ -597,8 +660,8 @@
irun = 1;
dir1 = [pdir 'horz_' sprintf('%2.2i',irun) '/'];
smodel1 = 'm00';
- smodel2 = 'm16'; % m16, m01, etc
- modlab = 'vs';
+ smodel2 = 'm01'; % m16, m01, etc
+ modlab = 'vb';
smodeltag = [smodel2 '_' smodel1];
% KEY: TWO DIFFERENT KINDS OF SUMMED KERNELS
More information about the CIG-COMMITS
mailing list