[cig-commits] r16244 - in seismo/3D/ADJOINT_TOMO/iterate_adj/matlab: . matlab_scripts
carltape at geodynamics.org
carltape at geodynamics.org
Tue Feb 9 15:52:49 PST 2010
Author: carltape
Date: 2010-02-09 15:52:49 -0800 (Tue, 09 Feb 2010)
New Revision: 16244
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/writeCMT_psmeca.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_station_SPECFEM.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_window_chi_all.m
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m
seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
Log:
Updated matlab scripts to reflect new ordering in window_chi measurement vector that is output from the adjoint source code measure_adj.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m 2010-02-09 23:44:05 UTC (rev 16243)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compare_misfit.m 2010-02-09 23:52:49 UTC (rev 16244)
@@ -18,7 +18,7 @@
format compact
% add path to additional matlab scripts
-path(path,'./matlab_scripts')
+path(path,[pwd '/matlab_scripts']);
dir0 = '/home/carltape/RUNS/';
@@ -28,7 +28,7 @@
iwrite = 0;
isciencehist = 0; % 0 or 1
-odirsci = '/home/carltape/manuscripts/2009/tomo_gji/latex/figures/scripts/hist_data/';
+odirpub = '/home/carltape/manuscripts/2009/tomo_gji/latex/figures/scripts/hist_data/';
DTMIN = 1.0; % see description below (1.0 or 0.0 for socal)
@@ -494,7 +494,7 @@
% fpre = 'hist_all_wdiff';
% %fpre = 'hist_T006_T030_DT';
%
-% fname = [odirsci fpre '_' DTMIN_tag '_' dtag '_' wtag '_' klabs{kk} '.dat'];
+% fname = [odirpub fpre '_' DTMIN_tag '_' dtag '_' wtag '_' klabs{kk} '.dat'];
% fid = fopen(fname,'w');
% pp = 0;
% for ii = 1:length(dsort)
@@ -524,7 +524,7 @@
klabs = {'DT2','DT_chi2','seis_chi1','seis_chi2','win_chi1','win_chi2','VRseis','VRwin'};
kind = [-26 30 31 32 33 34 -35 -36]; % columns to sort, corresponding to klabs
else
- odir = odirsci;
+ odir = odirpub;
klabs = {'seis_chi1','win_chi1'};
kind = [31 33];
end
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2010-02-09 23:44:05 UTC (rev 16243)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/compute_misfit.m 2010-02-09 23:52:49 UTC (rev 16244)
@@ -1,6 +1,6 @@
%
% compute_misfit.m
-% CARL TAPE, 06-Aug-2008
+% CARL TAPE, 14-Oct-2009
% printed xxx
%
% This script reads in a set of window_chi output measurement files from
@@ -26,12 +26,13 @@
dir0 = '/home/carltape/RUNS/';
% add path to additional matlab scripts
-path(path,'./matlab_scripts')
+path(path,[pwd '/matlab_scripts']);
% min and max periods for the different bandpassed datasets
-%Tminvec = [2 6]; Tmaxvec = [30 30];
-%Tminvec = [2 3 6]; Tmaxvec = [30 30 30];
-Tminvec = [6]; Tmaxvec = [30];
+%Trange = [ 2 30 ; 3 30 ; 6 30 ];
+Trange = [ 6 30 ];
+Tminvec = Trange(:,1);
+Tmaxvec = Trange(:,2);
comps = {'BHZ','BHR','BHT'};
% regardless of the component label for the DATA, the measurement code
@@ -45,6 +46,12 @@
%sTbp{tt} = sprintf('[%is,%is]',Tminvec(tt),Tmaxvec(tt));
sTbp{tt} = sprintf('%i-%is',Tminvec(tt),Tmaxvec(tt));
end
+
+% file label
+stBtag = 'BP';
+for ii=1:nper
+ stBtag = [stBtag num2str(Tminvec(ii))];
+end
%-------------------------
% USER INPUT
@@ -79,9 +86,8 @@
%eid_file = ['/home/carltape/results/EID_LISTS/kernels_run_' stmod];
%eid_file = ['/home/carltape/results/EID_LISTS/kernels_use_' stmod];
-eid_file = '/home/carltape/results/EID_LISTS/eids_simulation';
-%eid_file = '/home/carltape/results/EID_LISTS/eids_tomo';
-%eid_file = '/home/carltape/results/EID_LISTS/eids_extra';
+etag = 'extra'; % simulation, tomo, extra
+eid_file = ['/home/carltape/results/EID_LISTS/eids_' etag];
%-------------------------
% read in list of event IDs and sources
@@ -174,7 +180,7 @@
if iread == 1
% range of events
imin = 1; imax = nevent; % default
- %imin = 139; imax = imin;
+ %imin = 139; imax = imin; % 9818433 for simulation is index 139
meas_array = read_window_chi_all(imin,imax,stmod,Tminvec,Tmaxvec,dir0,eids,strec,stnet,comps);
save('meas_array.mat','meas_array');
@@ -184,14 +190,43 @@
whos meas_array
+% comparing CC measurments for DA and DT
+if 0==1
+ figure; nr=2; nc=1;
+ subplot(nr,nc,1);
+ x = abs(meas_array(:,21));
+ edges = [0:0.1:2.0]; ylims = [0 0.2];
+ Nbin = plot_histo(x,edges);
+ xlabel('|DA - CC for each window|'); ylim(ylims); grid on;
+
+ subplot(nr,nc,2);
+ x = abs(meas_array(:,19));
+ edges = [0:0.25:6.0]; ylims = [0 0.4];
+ Nbin = plot_histo(x,edges);
+ xlabel('|DT - CC for each window|'); ylim(ylims); grid on;
+
+ odirpub = '/home/carltape/manuscripts/2009/tomo_gji/latex/figures/scripts/hist_data/';
+ ofile = [odirpub 'window_chi_' stmod '_' stBtag '_' etag];
+ write_window_chi_all(meas_array,ofile,stmod,Tminvec,Tmaxvec,eids,strec,stnet,comps);
+
+ break
+end
+
% comparing multitaper measurments with cross-correlation measurements
if 0==1
iMT = find( meas_array(:,25) ~= 0 );
-
- figure; plot( meas_array(iMT,21),meas_array(iMT,25),'.');
+ amx = max([meas_array(iMT,21) ; meas_array(iMT,25)]);
+ length(iMT)
+
+ %inds = find( and( (abs( meas_array(:,25) - -0.31 ) <= 0.01) == 1, (abs( meas_array(:,21) - 0.19 ) <= 0.01) == 1 ) )
+ %display_meas(meas_array(inds,:),Tminvec,eids,strec,comps)
+
+ figure; hold on;
+ plot( meas_array(iMT,21),meas_array(iMT,25),'.');
+ plot([-1 1]*amx,[-1 1]*amx,'r--');
xlabel('DA - CC'); ylabel('DA - MT');
- axis equal; grid on
-
+ axis equal, axis([-1 1 -1 1]*amx); grid on
+
figure; nr=2; nc=1;
subplot(nr,nc,1);
x = meas_array(iMT,21);
@@ -206,9 +241,13 @@
%----------------------------------------
- figure; plot( meas_array(iMT,19),meas_array(iMT,23),'.');
+ amx = max([meas_array(iMT,19) ; meas_array(iMT,23)]);
+
+ figure; hold on;
+ plot( meas_array(iMT,19),meas_array(iMT,23),'.');
+ plot([-1 1]*amx,[-1 1]*amx,'r--');
xlabel('DT - CC'); ylabel('DT - MT');
- axis equal; grid on
+ axis equal, axis([-1 1 -1 1]*amx); grid on
figure; nr=2; nc=1;
subplot(nr,nc,1);
@@ -222,6 +261,16 @@
Nbin = plot_histo(x,edges);
xlabel('DT - MT for each window'); ylim(ylims); grid on;
+ % check differences between MT and CC measurements
+ %itest = find( abs(meas_array(iMT,19) - meas_array(iMT,23)) > 1.5 );
+ %itest = find( abs(meas_array(iMT,21) - meas_array(iMT,25)) > 1.0 );
+ %display_meas(meas_array(iMT(itest),:),Tminvec,eids,strec,comps);
+
+ itest = find( abs(meas_array(:,19) ) > 6.0 );
+ %itest = find( meas_array(:,19) < -3.0 );
+ %itest = find( meas_array(:,21) < -1.3 );
+ display_meas(meas_array(itest,:),Tminvec,eids,strec,comps);
+
break
end
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m 2010-02-09 23:44:05 UTC (rev 16243)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/display_meas.m 2010-02-09 23:52:49 UTC (rev 16244)
@@ -1,6 +1,6 @@
%
% display_meas.m
-% CARL TAPE, 18-Mar-2008
+% CARL TAPE, 13-Oct-2009
% printed xxx
%
% /cig/seismo/3D/mt_measure_adj_work/scripts_tomo/matlab/
@@ -12,8 +12,8 @@
function display_meas(meas_array,Tvec,eids,strec,comps)
X = meas_array;
-[m,n] = size(X);
-if n ~= 25, error(' should be 25 columns in the input matrix'); end
+[n,m] = size(X);
+if m ~= 28, error(' should be 28 columns in the input matrix'); end
index_win = X(:,1);
index_per = X(:,2);
@@ -42,15 +42,21 @@
sigmaCC_dA = X(:,22);
measMT_dT = X(:,23);
sigmaMT_dT = X(:,24);
-tr_chi = X(:,25);
+measMT_dA = X(:,25);
+sigmaMT_dA = X(:,26);
+tr_chi = X(:,27);
+am_chi = X(:,28);
disp('----------------');
disp('INDEX (window, period, event, receiver, component)');
-for ii = 1:m
+for ii = 1:n
disp(sprintf('%10s --> %6.1f %7s %6s %3i DT = %6.2f +- %6.2f DA = %6.2f chi = %6.2f -- %i',...
char(eids{index_event(ii)}),Tvec(index_per(ii)),char(strec{index_rec(ii)}),...
char(comps{index_comp(ii)}),isub_win(ii),...
measCC_dT(ii),sigmaCC_dT(ii),measCC_dA(ii),tr_chi(ii),index_win(ii) ));
+
+ % optional: list MT measurment values
+ %disp(sprintf('%40s DT = %6.2f +- %6.2f DA = %6.2f',' ',measMT_dT(ii),sigmaMT_dT(ii),measMT_dA(ii)));
end
disp('INDEX (window, period, event, receiver, component)');
disp('====================');
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/writeCMT_psmeca.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/writeCMT_psmeca.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/writeCMT_psmeca.m 2010-02-09 23:52:49 UTC (rev 16244)
@@ -0,0 +1,91 @@
+%
+% writeCMT_psmeca.m
+% CARL TAPE, 02-Feb-2007
+% printed xxx
+%
+% This inputs a set of moment tensors in units of N-m and outputs four
+% different files for plotting in GMT using psmeca, which assumes units of
+% dyne-cm. The four files differ only in what label to use plotting above
+% the beach balls.
+%
+% moment tensor M = [Mrr Mtt Mpp Mrt Mrp Mtp]
+% From psmeca (-Sm) : mrr, mtt, mff, mrt, mrf, mtf in 10*exponent dynes-cm
+%
+% calls CMT2m0.m
+% called by slab_readCMT.m, socal_quakes_2007.m
+%
+
+%function writeCMT_psmeca(filename,date,lat,lon,dep,M,eid,opts)
+function writeCMT_psmeca(filename,date,lat,lon,dep,M,eid,isource)
+
+ncmt = length(date);
+
+% make sure M is N by 6
+[a,b] = size(M); if b ~= 6, M = M'; end
+[a,b] = size(M); if b ~= 6, error('dimension of M must be N by 6'); end
+
+% if there is no isource input, then just default to all ones
+if nargin < 8, isource = ones(ncmt,1); end
+
+% convert moment tensor from N-m to dyne-cm
+M = 1e7 * M;
+
+% exponent for computing magnitude in psmeca
+M0 = CMT2m0(1,M);
+iexp_all = floor(log10(M0));
+
+% for labeling the moment magnitude
+Mw = m02mw(1,CMT2m0(1,M*1e-7)); % M0 must be in N-m
+
+% options
+%ilab = opts{1}; % type of label for beach ball
+
+% write 6 different versions of the file, each having different labels
+for ilab = 0:5
+ if ilab==0, ext = ''; end
+ if ilab==1, ext = '_eid'; end
+ if ilab==2, ext = '_Mw'; end
+ if ilab==3, ext = '_year'; end
+ if ilab==4, ext = '_date'; end
+ if ilab==5, ext = '_all'; end
+ disp(['writeCMT_psmeca.m : extension is ' ext]);
+
+ % write to file for GMT plotting
+ file1 = [filename '_psmeca' ext];
+ fid = fopen(file1,'w');
+ for ii = 1:ncmt
+
+ % title for beach ball
+ switch ilab
+ case 0, cmtlabel = ' ';
+ case 1, cmtlabel = char(eid{ii});
+ case 2, cmtlabel = sprintf('%.2f',Mw(ii));
+ case 3, cmtlabel = datestr(date(ii),10);
+ case 4, cmtlabel = datestr(date(ii),29);
+ case 5, cmtlabel = sprintf('%4i-M%.1f-%i',year(date(ii)),Mw(ii),isource(ii));
+ end
+
+ fac = 10^-iexp_all(ii);
+
+ fprintf(fid,'%14.6f%14.6f%14.6f%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%4i%14.6e%14.6e%16s\n',...
+ lon(ii), lat(ii), dep(ii),...
+ M(ii,1)*fac, M(ii,2)*fac, M(ii,3)*fac,...
+ M(ii,4)*fac, M(ii,5)*fac, M(ii,6)*fac,...
+ iexp_all(ii),...
+ lon(ii), lat(ii), cmtlabel);
+ end
+ fclose(fid);
+end
+
+% write a list of event IDs
+file2 = [filename '_eid'];
+fid = fopen(file2,'w');
+for ii=1:ncmt, fprintf(fid,'%s\n',char(eid{ii})); end
+fclose(fid);
+
+disp(' writing psmeca file for GMT plotting...');
+disp([' output file : ' file1]);
+disp([' output file : ' file2]);
+disp([' number of CMT solutions : ' num2str(ncmt)]);
+
+%======================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_station_SPECFEM.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_station_SPECFEM.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_station_SPECFEM.m 2010-02-09 23:52:49 UTC (rev 16244)
@@ -0,0 +1,47 @@
+%
+% function write_station_SPECFEM(filename,rlon,rlat,relev,rburial,stnm,netwk)
+% CARL TAPE, 15-March-2007
+% printed xxx
+%
+% This function inputs a set of values and names for stations, and it
+% outputs a stations file suitable for:
+% (1) plotting in GMT
+% (2) running SPECFEM3D
+%
+% 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 > stemp \n";
+% print CSH "psxy stemp ${station_info2} $J $R -K -V -O >> $psfile \n";
+%
+% calls xxx
+% called by xxx
+%
+
+function write_station_SPECFEM(filename,rlon,rlat,relev,rburial,stnm,netwk)
+
+% number of receivers
+nrec = length(rlat);
+
+% GMT format
+ofile = [filename '_gmt'];
+disp([' Writing the stations file ' ofile ]);
+fid = fopen(ofile,'w');
+for ii=1:nrec
+ fprintf(fid,'%14.6f%14.6f%7s%4s%12.2f%12.2f\n',...
+ rlon(ii),rlat(ii),stnm{ii},netwk{ii},relev(ii),rburial(ii));
+end
+fclose(fid);
+
+% SPECFEM format
+ofile = [filename '_specfem'];
+disp([' Writing the stations file ' ofile ]);
+fid = fopen(ofile,'w');
+fprintf(fid,'%10i\n',nrec);
+for ii=1:nrec
+ fprintf(fid,'%7s%4s%14.6f%14.6f%12.2f%12.2f\n',...
+ stnm{ii},netwk{ii},rlat(ii),rlon(ii),relev(ii),rburial(ii));
+end
+fclose(fid);
+
+%=======================================================================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_window_chi_all.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_window_chi_all.m (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/matlab_scripts/write_window_chi_all.m 2010-02-09 23:52:49 UTC (rev 16244)
@@ -0,0 +1,93 @@
+%
+% write_window_chi_all.m
+% CARL TAPE, 30-July-2009
+% printed xxx
+%
+% This reads in an array of window_chi misfit files and outputs a text file
+% that can be used for additional analysis and plots (e.g., GMT histograms).
+%
+% calls xxx
+% called by xxx
+%
+
+function write_window_chi_all(X,ofile,stmod,Tminvec,Tmaxvec,eids,strec,stnet,comps)
+
+% total number of measurements
+n = length(X);
+
+index_win = X(:,1);
+index_per = X(:,2);
+index_event = X(:,3);
+index_net = X(:,4);
+index_rec = X(:,5);
+index_comp = X(:,6);
+isub_win = X(:,7);
+iker = X(:,8);
+%---------------------------
+seisdur = X(:,9);
+windur = X(:,10);
+T_pmax_dat = X(:,11);
+T_pmax_syn = X(:,12);
+%---------------------------
+seis_d2 = X(:,13);
+win_d2 = X(:,14);
+seis_s2 = X(:,15);
+win_s2 = X(:,16);
+seis_diff2 = X(:,17);
+win_diff2 = X(:,18);
+%---------------------------
+measCC_dT = X(:,19);
+sigmaCC_dT = X(:,20);
+measCC_dA = X(:,21);
+sigmaCC_dA = X(:,22);
+measMT_dT = X(:,23);
+sigmaMT_dT = X(:,24);
+measMT_dA = X(:,25);
+sigmaMT_dA = X(:,26);
+tr_chi = X(:,27);
+am_chi = X(:,28);
+
+%n = 1000; % testing
+
+fid = fopen(ofile,'w');
+
+fprintf(fid,['%4s%10s%7s%7s%9s%5s%3s%7s%7s%6s%6s' ...
+ repmat('%10s',1,6) repmat('%7s',1,10) '\n'],...
+ 'mod','eid','Tmin','Tmax','rec','comp','iw',...
+ 'sdur','wdur','Tdat','Tsyn','sd2','wd2','ss2','ws2','sd2','wd2',...
+ 'DTCC','sigma','DACC','sigma','DTMT','sigma','DAMT','sigma',...
+ 'tr_chi','am_chi');
+
+for ii = 1:n
+ fprintf(fid,['%4s%10s%7.1f%7.1f%9s%5s%3i%7.1f%7.1f%6.1f%6.1f' ...
+ repmat('%10.2e',1,6) repmat('%7.3f',1,10) '\n'],...
+ stmod,char(eids{index_event(ii)}),Tminvec(index_per(ii)),Tmaxvec(index_per(ii)),...
+ char(strec{index_rec(ii)}),char(comps{index_comp(ii)}),isub_win(ii),...
+ X(ii,9:12),X(ii,13:18),X(ii,19:28) );
+end
+fclose(fid);
+
+%------------
+
+% find MT measurements
+nMT = 0; % temporary
+
+% values to go with the histograms
+% 1 - nwinCC
+% 2 - m16 DT-CC mean
+% 3 - m16 DT-CC std
+% 4 - m16 DT-MT mean
+% 5 - m16 DT-MT std
+% 6 - nwinMT
+% 7 - m16 DlnA mean
+% 8 - m16 DlnA std
+% 9 - m16 DlnA-MT mean
+% 10 - m16 DlnA-MT std
+
+disp('writing out statistic values...');
+fid = fopen([ofile '_stats'],'w');
+fprintf(fid,'%10i%12.4e%12.4e%12.4e%12.4e%10i%12.4e%12.4e%12.4e%12.4e',...
+ n,mean(measCC_dT),std(measCC_dT),0,0,nMT,mean(measCC_dA),std(measCC_dA),0,0);
+fclose(fid);
+
+%======================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2010-02-09 23:44:05 UTC (rev 16243)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/matlab/subspace_specfem.m 2010-02-09 23:52:49 UTC (rev 16244)
@@ -18,7 +18,7 @@
format compact
% add path to additional matlab scripts
-path(path,'./matlab_scripts')
+path(path,[pwd '/matlab_scripts']);
ax1 = [-121 -114 31 37];
stfm = '%4.4i';
More information about the CIG-COMMITS
mailing list