[cig-commits] r16250 - in seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab: . matlab_scripts
carltape at geodynamics.org
carltape at geodynamics.org
Tue Feb 9 16:23:05 PST 2010
Author: carltape
Date: 2010-02-09 16:23:04 -0800 (Tue, 09 Feb 2010)
New Revision: 16250
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/README
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/boxpts.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/colors.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/cpt2cmap.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/getsubset.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/merge_arrays.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/raypath.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_delimited_file.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_station_SPECFEM.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m
Removed:
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/utm2ll.m
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m
Log:
Added additional matlab scripts for generating plotting grids for tomographic cross sections. Moved matlab scripts into sub-folder.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/README 2010-02-10 00:23:04 UTC (rev 16250)
@@ -0,0 +1,12 @@
+Carl Tape
+22-May-2009
+
+/ADJOINT_TOMO/iterate_adj/model_plot/matlab/
+
+Matlab scripts for generating plotting grids that are used to extract values from SEM meshes. Then, back in Matlab, we create new data files that are read into GMT.
+
+The main script is socal_tomo_plots.m
+
+This calls several non-built-in Matlab scripts that could probably be cleaned up or, in some cases, replaced.
+
+------------------
\ No newline at end of file
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m 2010-02-10 00:21:05 UTC (rev 16249)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -1,35 +0,0 @@
-%
-% 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
-
-%==========================================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m 2010-02-10 00:21:05 UTC (rev 16249)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -1,38 +0,0 @@
-%
-% 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) ]';
-
-%===================================================
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m (from rev 15075, seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/axes_utm2ll.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -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
+
+%==========================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/axes_utm2ll.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/boxpts.m (from rev 15075, seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/boxpts.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/boxpts.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/boxpts.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -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) ]';
+
+%===================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/boxpts.m
___________________________________________________________________
Name: svn:mergeinfo
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/colors.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/colors.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/colors.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -0,0 +1,53 @@
+%
+% colors.m
+% Carl Tape 01-Dec-2002
+% printed xxx
+%
+% This gets the color-mapping matrices.
+%
+% calls cpt2cmap.m
+% called by trigridN.m, plotcmapsN.m, raysAGU.m, xxx
+%
+
+% blue = (0,0,1)
+% red = (1,0,0)
+% green = (0,1,0)
+% yellow = (0,1,1)
+% white = (1,1,1)
+% white = (0,0,0)
+
+numc = 65;
+white = ones(numc, 3); % white (no good for contour plots)
+bw = gray(numc); % black-to-white
+wb = -gray(numc) + 1; % white-to-black
+nw = ceil(numc/2);
+inc = 1/(nw-1);
+br = [[ones(nw,1) [0:inc:1]' [0:inc:1]']; [[1-inc:-inc:0]' [1-inc:-inc:0]' ones(nw-1,1)]];
+ % blue-to-red (top-to-bottom), with white in the middle
+br = [[ones(nw,1) [0:inc:1]' [0:inc:1]']; [[1-inc:-inc:0]' [1-inc:-inc:0]' ones(nw-1,1)]];
+ % blue-to-red (top-to-bottom), with white in the middle
+for uy = 1:numc,
+ rb(uy,:) = br(numc+1-uy, :);
+end
+
+% blue = (1,0,0)
+% red = (0,0,1)
+% green = (0,1,0)
+% yellow = (0,1,1)
+% white = (1,1,1)
+
+% values from GMT 'seis' color palette
+seis = [170 0 0; 206 0 0; 243 0 0; 255 24 0; 255 60 0; 255 97 0;
+255 133 0; 255 170 0; 255 206 0; 255 243 0; 255 255 0; 255 255 0;
+231 255 4; 161 255 17; 90 255 30; 51 249 64; 13 242 99; 0 194 152;
+0 125 214; 0 68 248; 0 34 226];
+
+% values for color palette used in TW phase velocity maps
+ana = [202 80 0; 255 90 0; 255 110 0; 255 130 0; 255 150 0; 255 170 0; 255 190 0; ...
+ 255 205 0; 190 190 240; 170 170 240; 150 150 240; 130 130 240; 100 100 240; ...
+ 70 70 240; 30 30 220; 30 30 140];
+
+br = cpt2cmap(ana);
+seis = cpt2cmap(seis);
+
+%=========================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/cpt2cmap.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/cpt2cmap.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/cpt2cmap.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -0,0 +1,28 @@
+%
+% function outmat = cpt2cmap(inmat)
+% Carl Tape, 14-March-2005
+% printed xxx
+%
+% This function inputs a 3 column matrix of colors (RGB)
+% and outputs a matrix suitable for Matlab plots.
+%
+% calls xxx
+% called by xxx
+%
+
+function outmat = cpt2cmap(inmat)
+
+numc = 65;
+len = length(inmat);
+temp = inmat/max(max(inmat)); % normalize to 0-1 range
+rv = temp(:,1);
+gv = temp(:,2);
+bv = temp(:,3);
+
+% I don't know a better way to do this -- this interpolates to numc points.
+numc = 65;
+outmat = [ interp1(linspace(1,numc,len), rv', [1:numc])' ...
+ interp1(linspace(1,numc,len), gv', [1:numc])' ...
+ interp1(linspace(1,numc,len), bv', [1:numc])'];
+
+%==================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/getsubset.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/getsubset.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/getsubset.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -0,0 +1,45 @@
+%
+% function [inds,xval,yval] = getsubset(x,y,ax0)
+% CARL TAPE, 21-Sept-2006
+% printed xxx
+%
+%
+%
+% calls xxx
+% called vby test_getsubset.m
+%
+
+function [inds,xval,yval] = getsubset(x,y,ax0)
+
+inds = find( and( and( x >= ax0(1) , x <= ax0(2)), and(y >= ax0(3), y <= ax0(4))) );
+xval = x(inds);
+yval = y(inds);
+
+disp(sprintf('%i points in the subset out of %i',length(inds),length(x)));
+
+% % correct if xmin and xmax are reversed
+% if xmin > xmax
+% xtemp = xmax;
+% xmax = xmin;
+% xmin = xtemp;
+% end
+%
+% if xmin <= min(x)
+% iq1 = 1;
+% else
+% itemp = find(abs(x - xmin) < trsh);
+% iqt1 = round(length(itemp)/2);
+% iq1 = itemp(iqt1);
+% end
+%
+% if xmax >= max(x)
+% iq2 = length(x);
+% else
+% itemp = find(abs(x - xmax) < trsh); % find nearest x-values
+% iqt2 = round(length(itemp)/2); % take the middle index
+% iq2 = itemp(iqt2);
+% end
+%
+% iran = [iq1 : iq2];
+
+%======================================================
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/merge_arrays.m (from rev 15075, seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/merge_arrays.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/merge_arrays.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -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
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/merge_arrays.m
___________________________________________________________________
Name: svn:mergeinfo
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/raypath.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/raypath.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/raypath.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -0,0 +1,71 @@
+%
+% function [rdist,lon,lat] = raypath(lat1,lon1,preseg,postseg,dx)
+% CARL TAPE, 05-Dec-2008
+% printed xxx
+%
+% This function returns the values of points on a great-circle "ray path",
+% given a starting point and a finishing point. The ray path is
+% constructed from two segments, each extending from the starting point.
+% The spacing between points will not necessarily be the same for the two
+% segments, and the finishing point will not necessarily be exactly one of
+% the digitized points on the ray path.
+%
+% calls xxx
+% called by xxx
+%
+
+function [lon,lat,rdist,azis] = raypath(lon1,lat1,lon2,lat2,preseg,postseg,dx)
+
+% distance and azimuth from point 1 to point 2
+[rng, az] = distance(lat1,lon1,lat2,lon2);
+
+% endpoints of the ray path, including the extensions
+%[latA,lonA] = reckon(lat1,lon1,km2deg(preseg),wrapTo360(az+180) );
+%[latB,lonB] = reckon(lat1,lon1,rng + km2deg(postseg),wrapTo360(az+180) );
+
+% lengths of two segments and number of legs per segment
+dist1 = preseg;
+dist2 = postseg+deg2km(rng);
+n1 = ceil(dist1/dx)+1;
+n2 = ceil(dist2/dx)+1;
+
+% two segments
+[lats1,lons1] = track1(lat1,lon1,wrapTo360(az+180),km2deg(dist1),[],[],n1);
+[lats2,lons2] = track1(lat1,lon1,az,km2deg(dist2),[],[],n2);
+npt1 = length(lats1);
+npt2 = length(lats2);
+npt = npt1 + npt2;
+
+% connect everything
+lon = zeros(npt,1);
+lat = zeros(npt,1);
+lon = [flipud(lons1) ; lons2(2:end) ];
+lat = [flipud(lats1) ; lats2(2:end) ];
+n = length(lon);
+
+[ddist,azis] = distance(lat1*ones(n,1),lon1*ones(n,1),lat,lon);
+rdist = deg2km(ddist);
+rdist(1:npt1) = -rdist(1:npt1);
+
+%------------------------------------
+
+% EXAMPLE
+if 0==1
+ lon1 = -120; lat1 = 36;
+ lon2 = -117; lat2 = 33;
+ preseg = 50; postseg = 50;
+ dx = 10;
+
+ [lon,lat,rdist,azis] = raypath(lon1,lat1,lon2,lat2,preseg,postseg,dx);
+ figure; hold on;
+ plot(lon,lat,'.');
+ plot(lon1,lat1,'rp',lon2,lat2,'ko','markersize',16);
+ for ii=1:length(lon), text(lon(ii),lat(ii),num2str(ii)); end
+
+ for ii=1:length(lon)
+ disp(sprintf('lon, lat, dist, az: %10.3f %10.3f %10.3f %10.3f',...
+ lon(ii),lat(ii),rdist(ii),azis(ii)));
+ end
+end
+
+%======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_delimited_file.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_delimited_file.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_delimited_file.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -0,0 +1,45 @@
+%
+% [data,ilabs,labels] = read_delimited_file(filename,ncol);
+% CARL TAPE, 22-August-2007
+% printed xxx
+%
+% This function reads in a delimited file. It was designed to read files
+% that were output from GMT using the pscoast command.
+%
+% NOTE: THERE MUST BE A DELIMITER TO START AND END THE FILE.
+%
+% calls xxx
+% called by xxx
+%
+
+function [data,ilabs,labels,inds] = read_delimited_file(filename,ncol)
+
+disp([' reading ' filename '...']);
+
+% read all lines
+lines = textread(filename,'%s','delimiter','\n','whitespace','');
+nlines = length(lines);
+
+% read in the LUR file and string delimiters
+jj = 0; kk = 0;
+data = NaN * ones(nlines,ncol);
+for ii = 1:nlines
+ tline = lines(ii);
+ temp = str2num(char(tline));
+
+ if ~isempty(temp)
+ jj = jj + 1;
+ data(ii,:) = temp;
+ else
+ kk = kk + 1;
+ labels(kk) = tline;
+ ilabs(kk) = ii;
+ end
+end
+ilabs = ilabs(:);
+labels = labels(:);
+
+% indices for each segment
+inds = [ilabs(1:end-1)+1 ilabs(2:end)-1];
+
+%=========================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_station_SPECFEM.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_station_SPECFEM.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/read_station_SPECFEM.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -0,0 +1,23 @@
+%
+% function [rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(filename)
+% CARL TAPE, 15-March-2007
+% printed xxx
+%
+% This function inputs a STATIONS file formatted for SPECFEM3D.f90.
+% See also write_station_SPECFEM.m
+%
+% In Perl, this can be done as follows:
+% open(IN,${stations_seis}); @temp = <IN>; $nrec = $temp[0]; chomp($nrec);
+% print CSH "tail -n $nrec ${stations_seis} > stemp0 \n";
+% print CSH "awk '{print \$4,\$3}' stemp0 > temp \n";
+%
+% calls xxx
+% called by xxx
+%
+
+function [rlon,rlat,relev,rburial,stnm,netwk] = read_station_SPECFEM(filename)
+
+% read in the STATIONS file
+[stnm,netwk,rlat,rlon,relev,rburial] = textread(filename,'%s%s%f%f%f%f','headerlines',1);
+
+%=======================================================================
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m (from rev 15075, seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -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
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2ll.m
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m (from rev 15075, seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m)
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -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);
+
+%==========================================================================
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/matlab_scripts/utm2llvec.m
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m 2010-02-10 00:21:05 UTC (rev 16249)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/merge_arrays.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -1,32 +0,0 @@
-%
-% 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
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 2010-02-10 00:21:05 UTC (rev 16249)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/socal_tomo_plots.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -1,7 +1,6 @@
%
% socal_tomo_plots.m
-% CARL TAPE, 27-March-2009
-% printed xxx
+% Carl Tape, 27-March-2009
%
% This file generates xyz points for which we want to extract an
% interpolation value from a SPECFEM volumetric field. Once we have the
@@ -16,21 +15,30 @@
close all
format short, format compact
+%-----------------------------------
+
+% add path to additional matlab scripts
+path(path,[pwd '/matlab_scripts']);
+
+%-----------------------------------
+
% USER OUTPUT DIRECTORIES
dir0 = '/home/carltape/ADJOINT_TOMO/ADJOINT_TOMO_OUTPUT/';
pdir = [dir0 'model_plot_matlab_OUTPUT/'];
mdir = [dir0 'OUTPUT_subspace_save_m16/'];
% KEY: vertical or horizontal cross-sections
-ivert = 1;
+ivert = 0;
iplot_vert = 0;
ihorz = 0;
iplot_horz = 0;
iplot_horz_coverage = 0;
+iplot_horz_ker = 1;
-iwrite = 1;
+iwrite = 0;
ifigure = 0;
+imask = 1; % write out masked versions with NaN (=1) or not (=0)
%-----------------------------------
@@ -93,7 +101,7 @@
if iwrite==1
% list of depths for each index
- fid = fopen([dir1 'horz_xc_elevation_cuts'],'w');
+ fid = fopen([dir1 'horz_xc_cuts'],'w');
for ii=1:ndep
fprintf(fid,'%10i%16i\n',ii,depvec(ii));
end
@@ -109,7 +117,7 @@
% xyz files
for ii=1:ndep
zdep = depvec(ii)
- filename = sprintf('horz_xc_elevation_%3.3i',ii);
+ filename = sprintf('horz_xc_%3.3i',ii);
fid = fopen([dir1 filename],'w');
for kk=1:N
fprintf(fid,'%16.6e%16.6e%16.6e\n',Xvec(kk),Yvec(kk),zdep);
@@ -138,7 +146,8 @@
% 4 LON
% 5 LAT
irun = 5; % KEY: which set of cross-sections to make
- dir1 = [pdir 'vert_' sprintf('%2.2i',irun) '/']; % KEY COMMAND
+ stirun = sprintf('%2.2i',irun);
+ dir1 = [pdir 'vert_' stirun '/']; % KEY COMMAND
switch irun
@@ -308,7 +317,7 @@
sdeps = zeros(nrays,1);
end
- fid = fopen([dir1 'ALL_rays'],'w');
+ fid = fopen([dir1 'ALL_rays_' stirun],'w');
aziall = zeros(nrays,2);
for ii = 1:nrays
azi1 = azimuth(rlats(ii),rlons(ii),slats(ii),slons(ii));
@@ -352,7 +361,7 @@
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');
+ fid0 = fopen([dir1 'ALL_rays_' stirun '_fault_positions'],'w');
imin = 1; imax = nrays;
%imin = 80; imax = nrays;
@@ -474,7 +483,7 @@
if iwrite == 1
% xyz files
- filename = sprintf('vert_xc_%3.3i_input',ii);
+ filename = sprintf('vert_%2.2i_xc_%3.3i_input',irun,ii);
fid = fopen([dir1 filename],'w');
for kk=1:hnum % loop over horizontal points in the ray path
for jj=1:znum
@@ -485,7 +494,7 @@
fclose(fid);
% ray path
- filename = sprintf('vert_xc_%3.3i_ray_path',ii);
+ filename = sprintf('vert_%2.2i_xc_%3.3i_ray_path',irun,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',...
@@ -494,7 +503,7 @@
fclose(fid);
% lon-lat points of source and receiver
- filename = sprintf('vert_xc_%3.3i_A_B',ii);
+ filename = sprintf('vert_%2.2i_xc_%3.3i_A_B',irun,ii);
fid = fopen([dir1 filename],'w');
[dq,azq] = distance(slats(ii),slons(ii),rlats(ii),rlons(ii)); % distance to receiver
iflip = 0;
@@ -518,18 +527,21 @@
%irun = 1; pmax = 86;
%irun = 2; pmax = 32;
- irun = 3; pmax = 10;
+ %irun = 3; pmax = 10;
+ %irun = 4; pmax = 13;
+ irun = 5; pmax = 9;
- dir1 = [pdir 'vert_' sprintf('%2.2i',irun) '/'];
+ stirun = sprintf('%2.2i',irun);
+ dir1 = [pdir 'vert_' stirun '/'];
smodel1 = 'm00';
smodel2 = 'm16';
modlab1 = 'vs';
modlab2 = 'vb';
smodeltag = [smodel2 '_' smodel1];
- pmin = 1;
+ pmin = 1; % default
%pmin = 80; pmax = 60;
- %pmin = pmax; pmax = pmin;
+ %pmin = 5; pmax = pmin;
for ip = pmin:pmax
stip = sprintf('%3.3i',ip);
@@ -542,13 +554,14 @@
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'];
+ fz0 = [dir1 'vert_' stirun '_xc_'];
+ file0 = [fz0 stip '_A_B'];
+ file1 = [fz0 stip '_ray_path'];
+ file2 = [fz0 stip '_input'];
+ file3 = [fz0 tag3 '_' stip '.gmt'];
+ file4 = [fz0 tag4 '_' stip '.gmt'];
+ file5 = [fz0 tag5 '_' stip '.gmt'];
+ file6 = [fz0 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');
@@ -628,7 +641,7 @@
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');
+% fid = fopen([fz0 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));
@@ -636,7 +649,7 @@
% 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');
+ fid = fopen([fz0 '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),...
@@ -652,16 +665,17 @@
end
%----------------------------
-% PLOTTING CROSS-SECTIONS
+% PLOTTING MODEL CROSS-SECTIONS
if iplot_horz==1
colors;
irun = 1;
- dir1 = [pdir 'horz_' sprintf('%2.2i',irun) '/'];
+ stirun = sprintf('%2.2i',irun);
+ dir1 = [pdir 'horz_' stirun '/'];
smodel1 = 'm00';
smodel2 = 'm01'; % m16, m01, etc
- modlab = 'vb';
+ modlab = 'vs';
smodeltag = [smodel2 '_' smodel1];
% KEY: TWO DIFFERENT KINDS OF SUMMED KERNELS
@@ -675,9 +689,16 @@
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');
+ [inds,depvec] = textread([dir1 'horz_' stirun '_xc_cuts'],'%f%f');
+ [xlon,ylat,xutm,yutm] = textread([dir1 'horz_' stirun '_xc_lonlat'],'%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;
+
% replace values far west in regions of no coverage with NaN
%xvpts = [-119 -122 -122 -119];
%yvpts = [32 36 32 32];
@@ -685,45 +706,35 @@
%figure; hold on;
%plot(xlon,ylat,'b.');
- %plot(xlon(imask),ylat(imask),'ro');
+ %plot(xlon(indmask),ylat(indmask),'ro');
- pmin = 1; pmax = length(depvec);
- %pmin = 21; pmax = pmin;
- %pmin = 3; pmax = length(depvec);
+ %pmin = 1; pmax = length(depvec);
+ pmin = 3; pmax = pmin;
+ %pmin = 1; pmax = 12;
for ip = pmin:pmax
stip = sprintf('%3.3i',ip);
- disp(sprintf('%i out of the set %i - %i',ip,pmin,pmax));
+ disp(sprintf('%i out of the set %i to %i',ip,pmin,pmax));
% load output from cluster code
tag1 = [modlab '_' smodel1];
tag2 = [modlab '_' smodel2];
- file1 = [dir1 'horz_xc_elevation_' tag1 '_' stip '.gmt'];
+ file1 = [dir1 tag1 '/horz_xc_' tag1 '_' stip '.gmt'];
[junk1,junk2,junk3,cval1,cdist] = textread(file1,'%f%f%f%f%f');
- file2 = [dir1 'horz_xc_elevation_' tag2 '_' stip '.gmt'];
+ file2 = [dir1 tag2 '/horz_xc_' 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;
+ if imask==1
+ % load coverage mask
+ tagm = [modlab '_m16_' covertag];
+ %stip0=stip;
+ %if ip==3, stip0 = '002'; end
+ filem = [dir1 tagm '/horz_xc_' tagm '_' stip '.gmt'];
+ [junk1,junk2,junk3,cvalm,cdist] = textread(filem,'%f%f%f%f%f');
+ end
- %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
+ % remove points above the topography -- these are set as negative
+ % numbers in the interpolation algorithm
ineg = find(cval1 < 0); cval1(ineg) = NaN;
ineg = find(cval2 < 0); cval2(ineg) = NaN;
@@ -734,15 +745,17 @@
%----------------------------------------------------
% COVERAGE MASK
- % take log10 of the field
- cvalmlog = log10(cvalm / Nfac);
+ if imask==1
+ % take log10 of the field
+ cvalmlog = log10(cvalm / Nfac);
+
+ % replace points outside coverage with a mask
+ indmask = find(cvalmlog <= cfac);
+ cval1(indmask) = NaN;
+ cval2(indmask) = NaN;
+ cval3(indmask) = NaN;
+ end
- % 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 ')']};
@@ -755,10 +768,6 @@
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);
@@ -783,7 +792,7 @@
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');
+ fid = fopen([dir1 'horz_xc_' modlab '_' smodeltag '_' stip '_mask' num2str(imask) '.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));
@@ -800,9 +809,10 @@
if iplot_horz_coverage==1
colors;
- irun = 1;
- dir1 = [pdir 'horz_' sprintf('%2.2i',irun) '/'];
- modlab = 'vb';
+ irun = 2;
+ stirun = sprintf('%2.2i',irun);
+ dir1 = [pdir 'horz_' stirun '/'];
+ modlab = 'vb_m16';
%smodeltag = [smodel2 '_' smodel1];
% load number of window picks per event
@@ -812,25 +822,26 @@
% KEY: TWO DIFFERENT KINDS OF SUMMED KERNELS
cfac = log10(1e-16);
- smodeltag = 'coverage_sum_abs'; Nfac = nwin;
- %smodeltag = 'coverage_sum_abs_nwin'; Nfac = nevent;
+ %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');
+ [inds,depvec] = textread([dir1 'horz_' stirun '_xc_cuts'],'%f%f');
+ [xlon,ylat,xutm,yutm] = textread([dir1 'horz_' stirun '_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);
+ %indmask = inpolygon(xlon,ylat,xvpts,yvpts);
%figure; hold on;
%plot(xlon,ylat,'b.');
- %plot(xlon(imask),ylat(imask),'ro');
+ %plot(xlon(indmask),ylat(indmask),'ro');
pmin = 1; pmax = length(depvec);
- %pmin = 21; pmax = pmin;
+ %pmin = 1; pmax = pmin;
+ pmin = 13; pmax = length(depvec);
for ip = pmin:pmax
stip = sprintf('%3.3i',ip);
@@ -838,7 +849,7 @@
% load output from cluster code
tag1 = [modlab '_' smodeltag];
- file1 = [dir1 'horz_xc_elevation_' tag1 '_' stip '.gmt'];
+ file1 = [dir1 tag1 '/horz_xc_' 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,
@@ -858,9 +869,9 @@
cval2 = log10(cval1);
% replace points outside coverage with a mask
- imask = find(cval2 <= cfac);
- cval3 = cval2; cval3(imask) = NaN;
- cval1mask = cval1; cval1mask(imask) = NaN;
+ indmask = find(cval2 <= cfac);
+ cval3 = cval2; cval3(indmask) = NaN;
+ cval1mask = cval1; cval1mask(indmask) = NaN;
% convex hull -- convhull
@@ -919,4 +930,100 @@
end
+%----------------------------
+% PLOTTING KERNEL CROSS-SECTIONS
+
+if iplot_horz_ker == 1
+ colors;
+
+ irun = 2;
+ stirun = sprintf('%2.2i',irun) ;
+ dir1 = [pdir 'horz_' stirun '/'];
+
+ %fsmooth = 'h003km_v001km';
+ fsmooth = 'h002km_v001km';
+ %fsmooth = 'h000km_v000km';
+ ftags = {'Kmu_rayleigh','Kmu_reflection_one','Kmu_reflection_two'};
+
+ % load files needed for horizontal cross-sections
+ [inds,depvec] = textread([dir1 'horz_' stirun '_xc_cuts'],'%f%f');
+ [xlon,ylat,xutm,yutm] = textread([dir1 'horz_' stirun '_xc_lonlat'],'%f%f%f%f');
+
+ kvec = [1 5 9 13 17 21 25 29];
+
+ %pmin = 1; pmax = length(depvec);
+ %pmin = 13; pmax = pmin;
+ %pmin = 13; pmax = length(depvec);
+
+ for ik = 1:length(kvec)
+ ip = kvec(ik)
+ stip = sprintf('%3.3i',ip);
+ %disp(sprintf('%i out of the set %i to %i',ip,pmin,pmax));
+
+ % load output from cluster code
+ file1 = [dir1 'kernels/horz_xc_' ftags{1} '_' fsmooth '_' stip '.gmt'];
+ [junk1,junk2,junk3,K1,junk4] = textread(file1,'%f%f%f%f%f');
+ file2 = [dir1 'kernels/horz_xc_' ftags{2} '_' fsmooth '_' stip '.gmt'];
+ [junk1,junk2,junk3,K2,junk4] = textread(file2,'%f%f%f%f%f');
+ file3 = [dir1 'kernels/horz_xc_' ftags{3} '_' fsmooth '_' stip '.gmt'];
+ [junk1,junk2,junk3,K3,junk4] = textread(file3,'%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;
+
+ %----------------------------------------------------
+
+ clabs = {'rayleigh-1','rayleigh-2','rayleigh-3'};
+ for kk = 1:3
+ %for kk = 2:2
+ switch kk
+ case 1, cval = K1;
+ case 2, cval = K2;
+ case 3, cval = K3;
+ 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, K1, K2, K3
+ fid = fopen([dir1 'horz_xc_ker_' fsmooth '_' 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),K1(kk),K2(kk),K3(kk));
+ end
+ fclose(fid);
+ end
+ end
+
+end
+
%====================================================================
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m 2010-02-10 00:21:05 UTC (rev 16249)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2ll.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -1,171 +0,0 @@
-%
-% 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
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m 2010-02-10 00:21:05 UTC (rev 16249)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/model_plot/matlab/utm2llvec.m 2010-02-10 00:23:04 UTC (rev 16250)
@@ -1,29 +0,0 @@
-%
-% 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