[cig-commits] r15994 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: . PLOTTING matlab
carltape at geodynamics.org
carltape at geodynamics.org
Wed Nov 18 08:36:37 PST 2009
Author: carltape
Date: 2009-11-18 08:36:36 -0800 (Wed, 18 Nov 2009)
New Revision: 15994
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_geometry.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m
Log:
Minor changes for plotting first default simulation results.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2009-11-18 14:11:22 UTC (rev 15993)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2009-11-18 16:36:36 UTC (rev 15994)
@@ -4,7 +4,7 @@
#
# joint_inversion.pl
# Carl Tape
-# 21-Nov-2008
+# 18-Nov-2009
#
# This script plots figures for the joint source-structure inversions
# produced by the 2D SEM wave propagation code.
@@ -26,13 +26,24 @@
($colors,$iker,$irun0,$ipoly,$qmax,$ihessian) = @ARGV;
$iter = 0;
-# USER change
-$dirq = "/home/carltape/wave2d/SEM2D_iterate";
+#-----------------------------------
+# USER INPUT
-# directory with the "reference" conjugate gradient simulation is
-#$irun0_cg = $irun0;
-$irun0_cg = 600;
+# base directory
+$dirq = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work";
+# OPTIONAL: directory with the "reference" conjugate gradient simulation is
+$irun0_cg = $irun0;
+#$irun0_cg = 600;
+
+# resolution of color plots
+#$interp = "-I0.5m/0.5m -S4m";
+$interp = "-I2m/2m -S4m";
+
+$edir = "event_005"; # event for first event (sr.txt)
+
+#-----------------------------------
+
$sth = sprintf("%1i",$ihessian);
# directories
@@ -42,14 +53,11 @@
$odir = "${dirq}/OUTPUT/run_";
if (not -e $ddir) {die("Check if ddir $ddir exist or not\n");}
-$edir = "event_005"; # event for first event (sr.txt)
-
# misfit function variable
$misfitvar = "\@~\143\@~";
$misfitvar = "S";
-# resolution of color plots
-$interp = "-I0.5m/0.5m -S4m"; # key information
+# grd file for interpolation
$grdfile = "temp.grd";
$mfile_dat_str = "structure_dat.dat";
@@ -641,7 +649,7 @@
print CSH "pstext -N $J_title $R_title -O -V >>$psfile<<EOF\n 10 10 $fsize_title 0 $fontno CM $title \nEOF\n"; # FINISH
print CSH "convert $psfile $jpgfile\n";
print CSH "echo done with $psfile\n";
-print CSH "ghostview $psfile &\n";
+print CSH "gv $psfile &\n";
close (CSH);
system("csh -f $cshfile");
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_geometry.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_geometry.pl 2009-11-18 14:11:22 UTC (rev 15993)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_geometry.pl 2009-11-18 16:36:36 UTC (rev 15994)
@@ -4,13 +4,13 @@
#
# plot_geometry.pl
# Carl Tape
-# 19-Nov-2008
+# 18-Nov-2009
#
# This plots the source-receiver geometry and source time function (STF) for two possible cases
# 1. PSV wavefield, view from side, three components of STF
# 2. SH wavefield, top view from surface, one component of STF
#
-# EXAMPLES (from /home/carltape/wave2d/2d_adjoint/PLOTTING/FIGURES):
+# EXAMPLES (from /ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/FIGURES):
# ../plot_geometry.pl 0 100 events_dat_lonlat.dat recs_lonlat.dat stf_00001 structure_syn.dat 10.0
# ../plot_geometry.pl 0 200 events_dat_lonlat.dat recs_lonlat.dat stf_00001 structure_syn.dat 10.0
# ../plot_geometry.pl 0 300 events_dat_lonlat.dat recs_lonlat.dat stf_00001 structure_syn.dat 10.0
@@ -22,8 +22,10 @@
if (@ARGV < 7) {die("Usage: plot_geometry.pl xxx \n");}
($ibody,$irun,$efile_syn,$rfile_syn,$stf_tag,$mfile_syn,$hdur) = @ARGV;
+# USER INPUT
+$rundir = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work";
+
# directories
-$rundir = "/home/carltape/wave2d/SEM2D_iterate"; # USER change
$plotdir = "${rundir}/PLOTTING";
$odir = "${rundir}/OUTPUT";
$figdir = "${plotdir}/FIGURES";
@@ -221,6 +223,6 @@
close (CSH);
system("csh -f $cshfile");
-system("ghostview $psfile &")
+system("gv $psfile &")
#=================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2009-11-18 14:11:22 UTC (rev 15993)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/plot_surf_model.pl 2009-11-18 16:36:36 UTC (rev 15994)
@@ -4,7 +4,7 @@
#
# plot_surf_model.pl
# Carl Tape
-# 23-Nov-2008
+# 18-Nov-2009
#
# Plot the source time function, the velocity model for the synthetics,
# and the velocity model for the data (optional).
@@ -30,6 +30,9 @@
($iopt,$irun0,$imodel,$ievents,$pbar,$stf_file1,$ifinite) = @ARGV;
$comp = 1;
+# USER INPUT
+$rundir = "/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work";
+
($pmax,$ptick) = split("/",$pbar);
($ievent_all,$ievent_one,$ievent) = split("/",$ievents);
$rec_file1 = "sr.txt";
@@ -38,7 +41,6 @@
$mfile_dat = "structure_dat.dat";
# directories
-$rundir = "/home/carltape/wave2d/SEM2D_iterate"; # USER change
$plotdir = "${rundir}/PLOTTING";
$odir = "${rundir}/OUTPUT";
$idir = "${rundir}/INPUT";
@@ -312,7 +314,7 @@
close (CSH);
system("csh -f $cshfile");
- system("ghostview $psfile &")
+ system("gv $psfile &")
#if($iopt <= 2) {system("xv $jpgfile &")}
#=================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2009-11-18 14:11:22 UTC (rev 15993)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2009-11-18 16:36:36 UTC (rev 15994)
@@ -22,6 +22,38 @@
SEM2D_iterate/gji_paper/codes_18_06Aug2006
--------------------------------
+DIRECTORIES AND FILES
+
+README
+INPUT
+OUTPUT
+Makefile
+Makefile_ifort
+src
+mod
+obj
+PLOTTING
+matlab
+gji_paper
+
+The default assumes that you have an OUTPUT directory linked from
+ ../SEM2D_iterate_OUTPUT
+
+--------------------------------
+COMPILATION
+
+Last compiled on 18-Nov-2009 using (gfortran --version)
+ GNU Fortran (Ubuntu 4.3.3-5ubuntu4) 4.3.3
+on a 64-bit Linux with Ubuntu 9.04 (jaunty).
+
+make clean
+make wave2d
+
+Lots of warning with the -fbounds-check flag.
+Need to modify the simple loops:
+ Warning: Deleted feature: Step expression in DO loop at (1) must be integer
+
+--------------------------------
EXAMPLES TO RUN
The user is encouraged to first try to run the set of seven examples.
@@ -46,20 +78,6 @@
/PLOTTING/FIGURES/Example7_run_0700.pdf
--------------------------------
-COMPILATION
-
-Last compiled on 12-Nov-2009 using (gfortran --version)
- GNU Fortran (Ubuntu 4.3.3-5ubuntu4) 4.3.3
-This was on a 64-bit Linux with Ubuntu 9.04 (jaunty).
-
-make clean
-make wave2d
-
-Lots of warning with the -fbounds-check flag.
-Need to modify the simple loops:
- Warning: Deleted feature: Step expression in DO loop at (1) must be integer
-
---------------------------------
DETAILS FOR RUNNING EACH EXAMPLE
EXAMPLE 0 (run_0000)
@@ -71,24 +89,29 @@
OUTPUT/run_0001/chi.dat 75.2856269798
OUTPUT/run_0002/chi.dat 31.8072588982
+------------------
EXAMPLE 1 (run_0100)
Make the following modifications to src/wave2d_constants.f90
IRUNZ = 0 --> IRUNZ = 100
NITERATION = 1 --> NITERATION = 16
Then compile and execute:
- make clean ; make wave2d ; wave2d
-After the simulation finishes
- cd PLOTTING/FIGURES
+ make clean ; make wave2d ; wave2d
+
Note that the output files are in the directories beginning with OUTPUT/run_0100.
-Then execute the example lines in plot_geometry.pl and plot_surf_model.pl.
-This will require first changing the variables for the main directory (search USER).
-Then open matlab and the file wave2d_cg_poly.m and wave2d_cg_figs.m.
-Run wave2d_cg_poly.m and then wave2d_cg_figs.m, which generate data files that
- will be used in joint_inversion.pl.
-Next execute the example line in joint_inversion.pl.
+ cd PLOTTING
+Change the user base directory
+in plot_geometry.pl, plot_surf_model.pl, joint_inversion.pl
+ cd FIGURES
+Execute the example lines --> geometry_0100.ps, model_0100.ps
+ cd ../../matlab/
+Open matlab and the files wave2d_cg_poly.m and wave2d_cg_figs.m.
+Run wave2d_cg_poly.m and then wave2d_cg_figs.m; this generates
+data files that will be used in joint_inversion.pl.
+ cd ../PLOTTING/FIGURES
+Execute the example line in joint_inversion.pl --> joint_inversion_0100_h0.ps
+------------------
EXAMPLE 2 (run_0200)
-
Make the following modifications to src/wave2d_constants.f90
IRUNZ = 100 --> IRUNZ = 200
PERT_STRUCT = 1 --> PERT_STRUCT = 0
@@ -101,6 +124,7 @@
make clean ; make wave2d ; wave2d
Follow steps from Example 1 to make figures.
+------------------
EXAMPLE 3 (run_0300)
Make the following modifications to src/wave2d_constants.f90
IRUNZ = 200 --> IRUNZ = 300
@@ -110,6 +134,7 @@
make clean ; make wave2d ; wave2d
Follow steps from Example 1 to make figures.
+------------------
EXAMPLE 4 (run_0400)
Make the following modification to src/wave2d_constants.f90
IRUNZ = 300 --> IRUNZ = 400
@@ -120,6 +145,7 @@
make clean ; make wave2d ; wave2d
Follow steps from Example 1 to make figures.
+------------------
EXAMPLE 5 (run_0500)
Make the following modification to src/wave2d_constants.f90
IRUNZ = 400 --> IRUNZ = 500
@@ -130,6 +156,7 @@
make clean ; make wave2d ; wave2d
Follow steps from Example 1 to make figures.
+------------------
EXAMPLE 6 (run_0600)
Make the following modification to src/wave2d_constants.f90
IRUNZ = 500 --> IRUNZ = 600
@@ -143,6 +170,7 @@
make clean ; make wave2d ; wave2d
Follow steps from Example 1 to make figures.
+------------------
EXAMPLE 7 (run_0700)
Make the following modification to src/wave2d_constants.f90
IRUNZ = 600 --> IRUNZ = 700
@@ -180,8 +208,6 @@
hmod = 1 --> hmod = 2
This tells the code where to look for the new model.
-As you can see, it is not ideal to go back-and-forth between Matlab and Fortran.
-However, many aspects are still in development, and Matlab provides superior
-quick graphics, which are crucial for debugging, checking, and testing.
+As you can see, it is not ideal to go back-and-forth between Matlab and Fortran. However, many aspects are still in development, and Matlab provides quick graphics that are crucial for debugging, checking, and testing.
==================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m 2009-11-18 14:11:22 UTC (rev 15993)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_figs.m 2009-11-18 16:36:36 UTC (rev 15994)
@@ -24,9 +24,9 @@
%n2 = 4; ax_chi = [n1-1 n2 10^-2 10^3];
-% USER: change this
-dir0 = '/home/carltape/wave2d/SEM2D_iterate/PLOTTING/DATA_FILES/';
-if ~exist(dir0), error([dir0 ' does not exist']); end
+% dir0 from wave2d_cg_poly.m
+dir1 = [dir0 '/PLOTTING/DATA_FILES/'];
+if ~exist(dir1), error([dir1 ' does not exist']); end
%==================================================================
@@ -47,7 +47,7 @@
if 1==1
% write curves to file
ww = ['poly_curve_' num2str(sprintf(stfm, irun0)) '.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
for ii=1:length(xpts)
fprintf(fid,'%16.7e%16.7e%16.7e%16.7e%16.7e%16.7e\n',...
xpts(ii),g1_line(ii),g2_line(ii),g1_quad_test(ii),g1_cube_fit(ii),g1_quad_fit(ii) );
@@ -56,7 +56,7 @@
% write points to file
ww = ['poly_points_' num2str(sprintf(stfm, irun0)) '.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
xptsQ = [x1 x2 x2 xmin xmin xmin 0];
yptsQ = [y1 0 y2 0 ymin chi chi];
for ii=1:length(xptsQ)
@@ -66,7 +66,7 @@
% write chi fit to file
ww = ['chi_' num2str(sprintf(stfm, irun0)) '.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
for ii=1:length(chis)
fprintf(fid,'%16.7e %16.7e \n',its(ii),chis(ii) );
end
@@ -74,7 +74,7 @@
% write chi fit to file
ww = ['chi_curve_' num2str(sprintf(stfm, irun0)) '.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
for ii=1:length(chifit_smooth)
fprintf(fid,'%16.7e %16.7e \n',its_smooth(ii),chifit_smooth(ii) );
end
@@ -82,7 +82,7 @@
% write var fit to file
ww = ['var_fit_' num2str(sprintf(stfm, irun0)) '.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
for ii=1:length(var_red1_fit)
fprintf(fid,'%16.7e %16.7e \n',x_var_fit(ii),var_red1_fit(ii) );
end
@@ -90,7 +90,7 @@
% write var to file
ww = ['var_' num2str(sprintf(stfm, irun0)) '.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
for ii=1:length(var_red1)
fprintf(fid,'%16.7e %16.7e \n',x_var(ii),var_red1(ii) );
end
@@ -98,7 +98,7 @@
% write axes to file
ww = ['axes_' num2str(sprintf(stfm, irun0)) '.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
fprintf(fid,'%16.7e %16.7e %16.7e %16.7e \n',axpoly(1),axpoly(2),axpoly(3),axpoly(4));
fprintf(fid,'%16.7e %16.7e %16.7e %16.7e \n',ax_chi(1),ax_chi(2),ax_chi(3),ax_chi(4));
fprintf(fid,'%16.7e %16.7e %16.7e %16.7e \n',ax_var(1),ax_var(2),ax_var(3),ax_var(4));
@@ -218,11 +218,11 @@
% load the chi values for this irun0 sequence
stfm = '%3.3i';
- dir0 = '/home/carltape/wave2d/2d_adjoint/OUTPUT/';
+ dir1 = '/home/carltape/wave2d/2d_adjoint/OUTPUT/';
chi_vec = zeros(nrun,1);
for ii=1:nrun
irun = irun0_vec(ii) + k*2;
- stf = [dir0 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [dir1 'run_' num2str(sprintf(stfm,irun)) '/'];
ww = 'summed_chi_all';
load([stf ww '.dat']);
chi = eval(ww);
@@ -271,18 +271,18 @@
stc = {'r','y','g','b','m','k'};
-dir0 = '/home/carltape/wave2d/2d_adjoint/gji_paper/figures/';
+dir1 = '/home/carltape/wave2d/2d_adjoint/gji_paper/figures/';
for ii=1:nrun
irun0 = irun0_vec(ii);
% fitting curve
ww = ['chi_curve_' num2str(sprintf(stfm, irun0)) ];
- load([dir0 ww '.dat']); temp = eval(ww);
+ load([dir1 ww '.dat']); temp = eval(ww);
its_smooth = temp(:,1);
chi_smooth = temp(:,2);
%ww = ['chi_' num2str(sprintf(stfm, irun0)) ];
- %load([dir0 ww '.dat']); temp = eval(ww);
+ %load([dir1 ww '.dat']); temp = eval(ww);
%its = temp(:,1);
%chis = temp(:,2);
@@ -295,7 +295,7 @@
for ii=1:nrun
irun0 = irun0_vec(ii);
ww = ['chi_' num2str(sprintf(stfm, irun0)) ];
- load([dir0 ww '.dat']); temp = eval(ww);
+ load([dir1 ww '.dat']); temp = eval(ww);
its = temp(:,1);
chis = temp(:,2);
plot(its,log10(chis),[stc{ii} '.'],'markersize',16);
@@ -303,17 +303,17 @@
if 0==1
- dir0 = '/home/store2/carltape/OUTPUT/run_001/';
- dir1 = [dir0 'event_001/'];
+ dir1 = '/home/store2/carltape/OUTPUT/run_001/';
+ dir2 = [dir1 'event_001/'];
% load source time function
- ww = 'stffor_00001_1'; load([dir1 ww]); temp = eval(ww);
+ ww = 'stffor_00001_1'; load([dir2 ww]); temp = eval(ww);
ti = temp(:,1); f = temp(:,2);
% load data, synthetics, adjoint source
- ww = 'dat_00001_1'; load([dir1 ww]); temp = eval(ww); s_dat = temp(:,2);
- ww = 'syn_00001_1'; load([dir1 ww]); temp = eval(ww); s_syn = temp(:,2);
- ww = 'stfadj_00001_1'; load([dir1 ww]); temp = eval(ww); f_adj = temp(:,2);
+ ww = 'dat_00001_1'; load([dir2 ww]); temp = eval(ww); s_dat = temp(:,2);
+ ww = 'syn_00001_1'; load([dir2 ww]); temp = eval(ww); s_syn = temp(:,2);
+ ww = 'stfadj_00001_1'; load([dir2 ww]); temp = eval(ww); f_adj = temp(:,2);
% compute velocity
sv_dat = gradient(s_dat,ti);
@@ -346,7 +346,7 @@
% write the data to a file
ww = ['time_series.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
for ii=1:length(ti)
fprintf(fid,'%16.7e%16.7e%16.7e%16.7e%16.7e%16.7e%16.7e\n',...
tall(ii,1),tall(ii,2),tall(ii,3),tall(ii,4),tall(ii,5),tall(ii,6),tall(ii,7));
@@ -355,7 +355,7 @@
% write the axes info to file
ww = ['time_series_axes.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
fprintf(fid,'%16i%16i%16i%16i\n',pwr_vec(1),pwr_vec(2),pwr_vec(3),pwr_vec(4));
fprintf(fid,'%16i%16i%16i%16i\n',cmx_vec(1),cmx_vec(2),cmx_vec(3),cmx_vec(4));
for ii=1:length(ax_mat)
@@ -366,13 +366,13 @@
end
if 0==1
- dir0 = '/home/store2/carltape/OUTPUT/run_002/'; % note run 002
- dir1 = [dir0 'event_001/'];
+ dir1 = '/home/store2/carltape/OUTPUT/run_002/'; % note run 002
+ dir2 = [dir1 'event_001/'];
% load source time function, synthetics, adjoint source
- ww = 'stffor_00001_1'; load([dir1 ww]); temp = eval(ww); ti = temp(:,1); f = temp(:,2);
- ww = 'syn_00001_1'; load([dir1 ww]); temp = eval(ww); s_syn = temp(:,2);
- ww = 'stfadj_00001_1'; load([dir1 ww]); temp = eval(ww); f_adj = temp(:,2);
+ ww = 'stffor_00001_1'; load([dir2 ww]); temp = eval(ww); ti = temp(:,1); f = temp(:,2);
+ ww = 'syn_00001_1'; load([dir2 ww]); temp = eval(ww); s_syn = temp(:,2);
+ ww = 'stfadj_00001_1'; load([dir2 ww]); temp = eval(ww); f_adj = temp(:,2);
% compute velocity
sv_syn = gradient(s_syn,ti);
@@ -409,7 +409,7 @@
% write the data to a file
ww = ['time_series.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
for ii=1:length(ti)
fprintf(fid,'%16.7e%16.7e%16.7e%16.7e%16.7e%16.7e%16.7e\n',...
tall(ii,1),tall(ii,2),tall(ii,3),tall(ii,4),tall(ii,5),tall(ii,6),tall(ii,7));
@@ -418,7 +418,7 @@
% write the axes info to file
ww = ['time_series_axes.dat'];
- fid = fopen([dir0 ww],'w');
+ fid = fopen([dir1 ww],'w');
fprintf(fid,'%16i%16i%16i%16i\n',pwr_vec(1),pwr_vec(2),pwr_vec(3),pwr_vec(4));
fprintf(fid,'%16i%16i%16i%16i\n',cmx_vec(1),cmx_vec(2),cmx_vec(3),cmx_vec(4));
for ii=1:length(ax_mat)
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2009-11-18 14:11:22 UTC (rev 15993)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2009-11-18 16:36:36 UTC (rev 15994)
@@ -24,8 +24,9 @@
close all
clear
-stf0 = '/home/carltape/wave2d/SEM2D_iterate/';
-if ~exist(stf0), error([stf0 ' does not exist']); end
+% USER: change this
+dir0 = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work/';
+if ~exist(dir0), error([dir0 ' does not exist']); end
colors;
npts = 100;
@@ -60,7 +61,7 @@
k = 0;
while k < 50
irun = irun0 + k;
- chifile = [stf0 odir 'run_' sprintf(stfm,irun) '/chi.dat'];
+ chifile = [dir0 odir 'run_' sprintf(stfm,irun) '/chi.dat'];
iexist = exist(chifile);
if iexist == 2
chi = load(chifile);
@@ -93,7 +94,7 @@
% chits = zeros(niter,1);
%
% % load initial chi value
-% stf = [stf0 odir 'run_' sprintf(stfm,irun0) '/'];
+% stf = [dir0 odir 'run_' sprintf(stfm,irun0) '/'];
% ww = 'summed_chi_all';
% load([stf ww '.dat']); chi = eval(ww); chis(1) = chi;
%
@@ -101,8 +102,8 @@
% for ii=1:niter
% irun = irun_vec(ii);
%
-% stft = [stf0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
-% stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
+% stft = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+% stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
% chit = load([stft 'summed_chi_all.dat']);
% chi = load([stf 'summed_chi_all.dat']);
% chis(ii+1) = chi;
@@ -133,7 +134,7 @@
%if or(ifinal_test==1, ii < niter)
if icubic == 1
ww = 'cubic_poly';
- stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
@@ -141,7 +142,7 @@
g1 = temp(5); g2 = temp(6);
else
ww = 'quad_poly';
- stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
@@ -167,7 +168,7 @@
% load actual chi-value computed from the next run
% wave2d: xx1,xx2,yy1,yy2,g1,g2,Pa,Pb,Pc,Pd,xmin
%ww = 'summed_chi_all';
- %stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
+ %stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
%load([stf ww '.dat']);
%chi = eval(ww);
%chis(ii+1) = chi;
@@ -206,7 +207,7 @@
%-------------------------
% write information to file for GMT plotting
- dir = [stf0 'gji_paper/figures/'];
+ dir = [dir0 'gji_paper/figures/'];
ww = 'chi_cubic_quad.dat';
fid = fopen([dir ww],'w');
for ii=1:nchi
@@ -410,7 +411,7 @@
% specify number of iterations to plot; load initial chi
niter = 1;
stirun = [' irun0 = ' num2str(irun0) '; niter = ' num2str(niter) ];
-stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun0)) '/'];
+stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun0)) '/'];
ww = 'chi';
load([stf ww '.dat']); chi = eval(ww); chis(1) = chi;
@@ -422,14 +423,14 @@
if icubic == 1
ww = 'cubic_poly';
- stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
y1 = temp(3); y2 = temp(4);
g1 = temp(5); g2 = temp(6);
else
ww = 'quad_poly';
- stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
+ stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun)) '/'];
load([stf ww '.dat']); temp = eval(ww)';
x1 = temp(1); x2 = temp(2);
y1 = temp(3); y2 = temp(4);
@@ -548,7 +549,7 @@
% load actual chi-value computed from the next run
% wave2d: xx1,xx2,yy1,yy2,g1,g2,Pa,Pb,Pc,Pd,xmin
ww = 'chi';
- stf = [stf0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
+ stf = [dir0 odir 'run_' num2str(sprintf(stfm,irun+1)) '/'];
load([stf ww '.dat']);
chi = eval(ww);
More information about the CIG-COMMITS
mailing list