[cig-commits] [commit] devel, master: More detailed instructions on how to use initial input plane waves (initialfield=true) an example is checked into EXAMPLES/init_plane prepare_initial_field is modified to allow plane wave incidence from the right (angleforce=negative vals) output snapshots are in jpg files which are easier to convert to avi movie (1c258d0)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 18 15:19:07 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branches: devel,master
Link : https://github.com/geodynamics/specfem2d/compare/fc67e6fd7ad890705b2b72b4b3c509accb22249e...e9ca46c40131588d89d7b0883250bc6584ce6b4c
>---------------------------------------------------------------
commit 1c258d0082ebcf7e2919625c12825a9c88bcde3b
Author: Qinya Liu <liuqy at physics.utoronto.ca>
Date: Wed Jul 20 21:36:27 2011 +0000
More detailed instructions on how to use initial input plane waves (initialfield=true)
an example is checked into EXAMPLES/init_plane
prepare_initial_field is modified to allow plane wave incidence from the right (angleforce=negative vals)
output snapshots are in jpg files which are easier to convert to avi movie
>---------------------------------------------------------------
1c258d0082ebcf7e2919625c12825a9c88bcde3b
canyon/SOURCE_canyon | 2 +-
.../Par_file_Slave | 40 +++++-----
.../Par_file_Slave_for | 40 +++++-----
.../Par_file_Slave_kernel | 40 +++++-----
.../SOURCE_Slave | 12 +--
M2_UPPA/process.sh => init_plane/adj.sh | 26 +------
init_plane/array_geo.m | 44 +++++++++++
.../interfaces_Slave.dat | 9 ++-
init_plane/isochron.m | 18 +++++
init_plane/plot_field.m | 14 ++++
init_plane/plot_kernel.m | 90 ++++++++++++++++++++++
{Tromp2005_kernel => init_plane}/process.sh | 6 +-
init_plane/readme | 12 +++
init_plane/run.bash | 12 +++
init_plane/run_for_ker.sh | 69 +++++++++++++++++
init_plane/select_adj.m | 74 ++++++++++++++++++
init_plane/set_source.bash | 29 +++++++
init_plane/stack_kernel.m | 58 ++++++++++++++
init_plane/ttplane.m | 50 ++++++++++++
19 files changed, 550 insertions(+), 95 deletions(-)
diff --git a/canyon/SOURCE_canyon b/canyon/SOURCE_canyon
index f3256d2..787006e 100644
--- a/canyon/SOURCE_canyon
+++ b/canyon/SOURCE_canyon
@@ -2,7 +2,7 @@
source_surf = .false. # source inside the medium or at the surface
xs = 2. # source location x in meters
zs = 3. # source location z in meters
-source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
+source_type = 1 # 1 for plane P waves, 2 for plane SV waves, 3 for Rayleigh wave
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 1.0 # dominant source frequency (Hz) if not Dirac or Heaviside
tshift = 0.0 # time shift when multi sources (if one source, must be zero)
diff --git a/Tromp2005/Par_file_Tromp2005_s100 b/init_plane/Par_file_Slave
similarity index 82%
copy from Tromp2005/Par_file_Tromp2005_s100
copy to init_plane/Par_file_Slave
index 67035ff..c5c6cf7 100644
--- a/Tromp2005/Par_file_Tromp2005_s100
+++ b/init_plane/Par_file_Slave
@@ -1,5 +1,5 @@
# title of job
-title = Tromp-Tape-Liu (GJI 2005)
+title = Slave Craton
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
@@ -10,7 +10,7 @@ nproc = 1 # number of processes
partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
ngnod = 4 # number of control nodes per element (4 or 9)
-initialfield = .false. # use a plane wave as source or not
+initialfield = .true. # use a plane wave as source or not
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
@@ -21,8 +21,8 @@ freq0 = 10 # frequency for viscous attenua
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 3000 # total number of time steps
-deltat = 2.0d-4 # duration of a time step
+nt = 6000 # total number of time steps
+deltat = 0.01 # duration of a time step
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
@@ -32,7 +32,7 @@ force_normal_to_surface = .false. # angleforce normal to surface
N_SLS = 2 # number of standard linear solids for attenuation
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
-# receiver line parameters for seismograms
+# eceiver line parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
@@ -40,19 +40,19 @@ anglerec = 0.d0 # angle to rotate components at
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line
-nrec = 1 # number of receivers
-xdeb = 1500. # first receiver x in meters
-zdeb = 400. # first receiver z in meters
-xfin = 700. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 0. # last receiver z in meters (ignored if onlyone receiver)
+nrec = 30
+xdeb = 500000.
+zdeb = 250000. # first receiver z in meters
+xfin = 100000.
+zfin = 250000. # last receiver z in meters (ignored if onlyone receiver)
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
-NTSTEP_BETWEEN_OUTPUT_INFO = 400 # display frequency in time steps
+NTSTEP_BETWEEN_OUTPUT_INFO = 200 # display frequency in time steps
output_postscript_snapshot = .true. # output Postscript snapshot of the results
output_color_image = .true. # output color image of the results
imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
-cutsnaps = 1. # minimum amplitude in % for snapshots
+cutsnaps = 5. # minimum amplitude in % for snapshots
meshvect = .true. # display mesh on vector plots or not
modelvect = .false. # display velocity model on vector plots
boundvect = .true. # display boundary conditions on plots
@@ -66,14 +66,15 @@ output_energy = .false. # compute and output acoustic a
output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
# velocity and density models
-nbmodels = 1 # nb of different models
+nbmodels = 2 # nb of different models
# define models as
# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 1 2600.d0 5800.d0 3198.6d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+1 1 2800.d0 6500.d0 3700d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+2 1 3200.d0 8300.d0 4700d0 0 0 10.d0 10.d0 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -97,12 +98,12 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = ../interfaces_Tromp2005.dat
+interfacesfile = ../interfaces_Slave.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
-xmax = 2000.d0 # abscissa of right side of the model
-nx = 80 # number of elements along X
+xmax = 600000.d0 # abscissa of right side of the model
+nx = 120 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
@@ -111,5 +112,6 @@ absorbtop = .false.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 1 # nb of regions and model number for each
-1 80 1 32 1
+nbregions = 2 # nb of regions and model number for each
+1 120 43 50 1
+1 120 1 42 2
diff --git a/Tromp2005_kernel/Par_file_Tromp2005 b/init_plane/Par_file_Slave_for
similarity index 82%
copy from Tromp2005_kernel/Par_file_Tromp2005
copy to init_plane/Par_file_Slave_for
index 83461c2..c5cc35f 100644
--- a/Tromp2005_kernel/Par_file_Tromp2005
+++ b/init_plane/Par_file_Slave_for
@@ -1,5 +1,5 @@
# title of job
-title = Tromp-Tape-Liu (GJI 2005)
+title = Slave Craton
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
@@ -10,7 +10,7 @@ nproc = 1 # number of processes
partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
ngnod = 4 # number of control nodes per element (4 or 9)
-initialfield = .false. # use a plane wave as source or not
+initialfield = .true. # use a plane wave as source or not
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
@@ -21,8 +21,8 @@ freq0 = 10 # frequency for viscous attenua
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 3000 # total number of time steps
-deltat = 2.0d-2 # duration of a time step
+nt = 6000 # total number of time steps
+deltat = 0.01 # duration of a time step
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
@@ -32,7 +32,7 @@ force_normal_to_surface = .false. # angleforce normal to surface
N_SLS = 2 # number of standard linear solids for attenuation
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
-# receiver line parameters for seismograms
+# eceiver line parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
@@ -40,19 +40,19 @@ anglerec = 0.d0 # angle to rotate components at
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line
-nrec = 1 # number of receivers
-xdeb = 150000. # first receiver x in meters
-zdeb = 40000. # first receiver z in meters
-xfin = 70000. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 0. # last receiver z in meters (ignored if onlyone receiver)
+nrec = 30
+xdeb = 500000.
+zdeb = 250000. # first receiver z in meters
+xfin = 100000.
+zfin = 250000. # last receiver z in meters (ignored if onlyone receiver)
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
-NTSTEP_BETWEEN_OUTPUT_INFO = 400 # display frequency in time steps
+NTSTEP_BETWEEN_OUTPUT_INFO = 200 # display frequency in time steps
output_postscript_snapshot = .true. # output Postscript snapshot of the results
output_color_image = .true. # output color image of the results
imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
-cutsnaps = 1. # minimum amplitude in % for snapshots
+cutsnaps = 5. # minimum amplitude in % for snapshots
meshvect = .true. # display mesh on vector plots or not
modelvect = .false. # display velocity model on vector plots
boundvect = .true. # display boundary conditions on plots
@@ -66,14 +66,15 @@ output_energy = .false. # compute and output acoustic a
output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
# velocity and density models
-nbmodels = 1 # nb of different models
+nbmodels = 2 # nb of different models
# define models as
# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 1 2600.d0 5800.d0 3198.6d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+1 1 2800.d0 6500.d0 3700d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+2 1 3200.d0 8300.d0 4700d0 0 0 10.d0 10.d0 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -97,12 +98,12 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = ../interfaces_Tromp2005.dat
+interfacesfile = ../interfaces_Slave.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
-xmax = 200000.d0 # abscissa of right side of the model
-nx = 80 # number of elements along X
+xmax = 600000.d0 # abscissa of right side of the model
+nx = 120 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
@@ -111,5 +112,6 @@ absorbtop = .false.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 1 # nb of regions and model number for each
-1 80 1 32 1
+nbregions = 2 # nb of regions and model number for each
+1 120 43 50 1
+1 120 1 42 2
diff --git a/Tromp2005/Par_file_Tromp2005_s100 b/init_plane/Par_file_Slave_kernel
similarity index 82%
copy from Tromp2005/Par_file_Tromp2005_s100
copy to init_plane/Par_file_Slave_kernel
index 67035ff..d60dd56 100644
--- a/Tromp2005/Par_file_Tromp2005_s100
+++ b/init_plane/Par_file_Slave_kernel
@@ -1,8 +1,8 @@
# title of job
-title = Tromp-Tape-Liu (GJI 2005)
+title = Slave Craton
# forward or adjoint simulation
-SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+SIMULATION_TYPE = 2 # 1 = forward, 2 = adjoint + kernels
SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
# parameters concerning partitioning
@@ -21,8 +21,8 @@ freq0 = 10 # frequency for viscous attenua
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 3000 # total number of time steps
-deltat = 2.0d-4 # duration of a time step
+nt = 6000 # total number of time steps
+deltat = 0.01 # duration of a time step
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
@@ -32,7 +32,7 @@ force_normal_to_surface = .false. # angleforce normal to surface
N_SLS = 2 # number of standard linear solids for attenuation
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
-# receiver line parameters for seismograms
+# eceiver line parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
@@ -40,19 +40,19 @@ anglerec = 0.d0 # angle to rotate components at
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line
-nrec = 1 # number of receivers
-xdeb = 1500. # first receiver x in meters
-zdeb = 400. # first receiver z in meters
-xfin = 700. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 0. # last receiver z in meters (ignored if onlyone receiver)
+nrec = 30
+xdeb = 500000.
+zdeb = 250000. # first receiver z in meters
+xfin = 100000.
+zfin = 250000. # last receiver z in meters (ignored if onlyone receiver)
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
-NTSTEP_BETWEEN_OUTPUT_INFO = 400 # display frequency in time steps
+NTSTEP_BETWEEN_OUTPUT_INFO = 200 # display frequency in time steps
output_postscript_snapshot = .true. # output Postscript snapshot of the results
output_color_image = .true. # output color image of the results
imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
-cutsnaps = 1. # minimum amplitude in % for snapshots
+cutsnaps = 5. # minimum amplitude in % for snapshots
meshvect = .true. # display mesh on vector plots or not
modelvect = .false. # display velocity model on vector plots
boundvect = .true. # display boundary conditions on plots
@@ -66,14 +66,15 @@ output_energy = .false. # compute and output acoustic a
output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
# velocity and density models
-nbmodels = 1 # nb of different models
+nbmodels = 2 # nb of different models
# define models as
# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 1 2600.d0 5800.d0 3198.6d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+1 1 2800.d0 6500.d0 3700d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+2 1 3200.d0 8300.d0 4700d0 0 0 10.d0 10.d0 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -97,12 +98,12 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = ../interfaces_Tromp2005.dat
+interfacesfile = ../interfaces_Slave.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
-xmax = 2000.d0 # abscissa of right side of the model
-nx = 80 # number of elements along X
+xmax = 600000.d0 # abscissa of right side of the model
+nx = 120 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
@@ -111,5 +112,6 @@ absorbtop = .false.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 1 # nb of regions and model number for each
-1 80 1 32 1
+nbregions = 2 # nb of regions and model number for each
+1 120 43 50 1
+1 120 1 42 2
diff --git a/Tromp2005/SOURCE_Tromp2005_s100 b/init_plane/SOURCE_Slave
similarity index 59%
copy from Tromp2005/SOURCE_Tromp2005_s100
copy to init_plane/SOURCE_Slave
index d6d6aae..12f7c3a 100644
--- a/Tromp2005/SOURCE_Tromp2005_s100
+++ b/init_plane/SOURCE_Slave
@@ -1,12 +1,12 @@
#source 1. The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
source_surf = .false. # source inside the medium or at the surface
-xs = 500. # source location x in meters
-zs = 400. # source location z in meters
-source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
-time_function_type = 2 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
-f0 = 42. # dominant source frequency (Hz) if not Dirac or Heaviside
+xs = -150000
+zs = 250000. # source location z in meters
+source_type = 1 # 1 for plane P waves, 2 for plane SV waves, 3 for Rayleigh wave
+time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0 = 0.4 # dominant source frequency (Hz) if not Dirac or Heaviside
tshift = 0.0 # time shift when multi sources (if one source, must be zero)
-angleforce = 270.0 # angle of the source (for a force only)
+angleforce = 23
Mxx = 1. # Mxx component (for a moment tensor source only)
Mzz = -1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
diff --git a/M2_UPPA/process.sh b/init_plane/adj.sh
similarity index 61%
copy from M2_UPPA/process.sh
copy to init_plane/adj.sh
index 4f94f18..07f1810 100755
--- a/M2_UPPA/process.sh
+++ b/init_plane/adj.sh
@@ -7,27 +7,6 @@
echo "running example: `date`"
currentdir=`pwd`
-echo
-echo "(will take about 2 minutes)"
-echo
-
-# sets up directory structure in current example directoy
-echo
-echo " setting up example..."
-echo
-
-mkdir -p OUTPUT_FILES
-mkdir -p DATA
-
-# sets up local DATA/ directory
-cd DATA/
-cp ../Par_file_M2_UPPA Par_file
-cp ../SOURCE_M2_UPPA SOURCE
-cd ../
-
-# cleans output files
-rm -rf OUTPUT_FILES/*
-
# compiles executables in root directory
cd ../../
make > tmp.log
@@ -39,8 +18,7 @@ ln -s ../../bin/xmeshfem2D
ln -s ../../bin/xspecfem2D
# stores setup
-cp DATA/Par_file OUTPUT_FILES/
-cp DATA/SOURCE OUTPUT_FILES/
+cp DATA/Par_file OUTPUT_FILES/Par_file_kernel
# runs database generation
echo
@@ -52,7 +30,7 @@ echo
echo
echo " running solver..."
echo
-./xspecfem2D > OUTPUT_FILES/output_solver.txt
+./xspecfem2D > OUTPUT_FILES/output_kernel.txt
# stores output
cp DATA/SOURCE_xz.dat OUTPUT_FILES/
diff --git a/init_plane/array_geo.m b/init_plane/array_geo.m
new file mode 100644
index 0000000..c37410e
--- /dev/null
+++ b/init_plane/array_geo.m
@@ -0,0 +1,44 @@
+% given the length of the station array L_s and the average Moho
+% depth H, figure out the appropriate geometry for a plane wave
+% with incident angle th (in degrees) that will propagate (as
+% initial condition to the free surface and cover the entire array
+
+close all; clear all;
+L_s=400; %km
+H=40; %km
+L_1=1./4*L_s; % km
+angle_inc=23;
+th=angle_inc*pi/180;
+
+% plane waves start below the Moho
+L_1=max(L_1,H*tan(th));
+
+x0=0.; z0=0.;
+x1=x0+L_1;
+x2=x1+L_s;
+xm=x2+x1;
+
+%
+fac=1.2;
+xp=x0-H*fac/tan(th);
+x4=x1-(x1-xp)*(sin(th).^2);
+z4=(x1-xp)*sin(th)*cos(th);
+x5=x2-(x2-xp)*(sin(th).^2);
+z5=(x2-xp)*sin(th)*cos(th);
+zm=z5*1.1;
+
+% draw the setup geometry
+set(0,'Defaultlinelinewidth',2);
+rectangle('position',[x0,z0,xm,zm],'EdgeColor','g'); hold on;
+% extending -x axis
+plot([x0,xp],[zm,zm],'g:');
+% rays going to x1 and x2 receivers
+plot([x1,x4],[zm,zm-z4]);
+plot([x2,x5],[zm,zm-z5]);
+% wavefront connecting (x4,z4), (x5,z5)
+plot([xp,x4,x5],[zm,zm-z4,zm-z5],'r:');
+% plot Moho
+plot([x0,xm],[zm-H,zm-H],'m')
+axis equal;
+
+xp,xm,zm
\ No newline at end of file
diff --git a/Tromp2005/interfaces_Tromp2005.dat b/init_plane/interfaces_Slave.dat
similarity index 89%
copy from Tromp2005/interfaces_Tromp2005.dat
copy to init_plane/interfaces_Slave.dat
index b8be387..81571fb 100644
--- a/Tromp2005/interfaces_Tromp2005.dat
+++ b/init_plane/interfaces_Slave.dat
@@ -6,13 +6,14 @@
# interface number 1 (bottom of the mesh)
2
0 0
- 200000 0
+ 600000 0
# interface number 5 (topography, top of the mesh)
2
- 0 80000
- 200000 80000
+ 0 250000
+ 600000 250000
#
# for each layer, we give the number of spectral elements in the vertical direction
#
# layer number 1
- 32
\ No newline at end of file
+50
+
diff --git a/init_plane/isochron.m b/init_plane/isochron.m
new file mode 100644
index 0000000..35bc022
--- /dev/null
+++ b/init_plane/isochron.m
@@ -0,0 +1,18 @@
+% read in forward seismograms
+
+s=load('OUTPUT_FILES/S0005.AA.BXZ.semd');
+t=s(:,1); disp=s(:,2);
+tmin=12.0; tmax=21.0;
+tshift=5.0; factor=0.2;
+dt=t(2)-t(1);
+ind=find(t>tmin & t<tmax);
+ind_shift=floor(tshift/dt);
+
+adj=zeros(size(disp));
+adj(ind+ind_shift)=disp(ind)*factor;
+
+plot(t,adj)
+
+fid = fopen('S0005.AA.BXZ.adj','w');
+fprintf(fid,'%6.2f %15.4g\n',[t';adj']); % the variable list has to be two rows and n columns
+fclose(fid);
\ No newline at end of file
diff --git a/init_plane/plot_field.m b/init_plane/plot_field.m
new file mode 100644
index 0000000..0eff305
--- /dev/null
+++ b/init_plane/plot_field.m
@@ -0,0 +1,14 @@
+close all; clear all
+sframe=200; eframe=5000; dframe=sframe
+i=0;
+for iframe = [5, sframe:dframe:eframe]
+ i = i+1
+ file=sprintf('OUTPUT_FILES/wavefield%07d_02_000.txt',iframe);
+ s=load(file);
+ if (i == 1)
+ x=s(:,1); z=s(:,2);
+ end
+ sx=s(:,3); sz=s(:,5);
+ figure(i)
+ scatter(x,z,30,sz); % or sx
+end
\ No newline at end of file
diff --git a/init_plane/plot_kernel.m b/init_plane/plot_kernel.m
new file mode 100644
index 0000000..eeb5d79
--- /dev/null
+++ b/init_plane/plot_kernel.m
@@ -0,0 +1,90 @@
+clear all; close all;
+LENGTH=250000; L1=LENGTH/1000; H=40;
+angle_inc=18;
+th=angle_inc*pi/180;
+dirr='OUTPUT_FILES_FORWARD/';
+dir = 'OUTPUT_FILES_KERNEL/';
+% receiver
+file=[dirr 'receiver.dat'];
+rec=load(file);
+rec=rec/LENGTH; nrec=length(rec);
+
+% kernel file
+file=[dir 'proc000000_rhop_alpha_beta_kernel.dat'];
+s=load(file);
+xx=s(:,1)/LENGTH;zz=s(:,2)/LENGTH;
+minx=min(xx); maxx=max(xx); minz=min(zz); maxz=max(zz);
+nx=400; nz=400;
+% define mesh
+xa=linspace(minx,maxx,nx);za=linspace(minz,maxz,nz);
+[X,Z] = meshgrid(xa,za);
+rhop=s(:,3); alpha=s(:,4); beta=s(:,5);
+
+subplot(3,1,1)
+%scatter(xx,zz,30,rhop,'filled'); hold on;
+FV=griddata(xx,zz,rhop,X,Z);
+imagesc(xa,za,FV); set(gca,'YDir','normal'); hold on;
+plot(rec(:,1),rec(:,2),'kx')
+text(0.2,0.2,'K_{\rho_p} ')
+quiver(0.2,1-H/L1, 0.1,0,'MaxHeadsize', 30)
+quiver(0,0,0.4*sin(th),0.4*cos(th),'MaxHeadsize', 30)
+title(['inc angle=', num2str(angle_inc), '; nrec=', num2str(nrec)]);
+jpeg_header='right-kernel-';
+jpeg_file=[jpeg_header num2str(angle_inc) '-' num2str(nrec) '.jpg'];
+
+subplot(3,1,2)
+%scatter(xx,zz,30,beta,'filled'); hold on;
+FV=griddata(xx,zz,beta,X,Z);
+imagesc(xa,za,FV); set(gca,'YDir','normal'); hold on;
+plot(rec(:,1),rec(:,2),'kx')
+plot(rec(:,1),rec(:,2),'kx')
+text(0.2,0.2,'K_\beta')
+quiver(0.2,1-H/L1, 0.1,0,'MaxHeadsize', 30)
+quiver(0,0,0.4*sin(th),0.4*cos(th),'MaxHeadsize', 30)
+
+
+subplot(3,1,3)
+%scatter(xx,zz,30,alpha,'filled'); hold on;
+FV=griddata(xx,zz,alpha,X,Z);
+imagesc(xa,za,FV); set(gca,'YDir','normal'); hold on;
+plot(rec(:,1),rec(:,2),'kx')
+text(0.2,0.2,'K_\alpha')
+quiver(0.2,1-H/L1, 0.1,0,'MaxHeadsize', 30)
+quiver(0,0,0.4*sin(th),0.4*cos(th),'MaxHeadsize', 30)
+
+print('-djpeg',jpeg_file)
+%file='angle_23/OUTPUT_FILES_KERNEL/proc000000_rhop_alpha_beta_kernel.dat';
+%s=load(file);
+%xx=s(:,1)/LENGTH;zz=s(:,2)/LENGTH;
+%rhop=s(:,3); alpha=s(:,4); beta=s(:,5);
+%field2=rhop;
+% define surface grid
+%maxv=max(abs(field))/3; ind=find(abs(field)>maxv);
+%field(ind)=maxv;
+%xmax=max(xx); zmax=max(zz);
+%xg = linspace(0,xmax,900);
+%zg = linspace(0,zmax,300);
+%[XI,YI] = meshgrid(xg,zg);
+%FI = griddata(xx,zz,field,XI,YI);
+%mesh(XI,YI,FI); view(2);
+
+%hold on;
+% receivers
+
+
+%subplot(3,1,2);
+%scatter(xx,zz,30,field2+field1,'filled');hold on;
+
+% plot assumed Moho at 40 km depth
+%plot([0,3],[160,160]/L1,'k:')
+
+% isochrons
+% $$$ th=23*pi/180;
+% $$$ beta=3.7; alpha=6.5;
+% $$$ phi=linspace(0,pi);
+% $$$ t0=5.0;
+% $$$ r=t0./(1/beta-1/alpha*cos(th+pi/2-phi));
+% $$$ r2=t0./(1/alpha-1/alpha*cos(th+pi/2-phi));
+% $$$ plot(r.*cos(phi)/L1+rec(1),rec(2)-r.*sin(phi)/L1); hold on;
+% $$$ plot(r2.*cos(phi)/L1+rec(1),rec(2)-r2.*sin(phi)/L1);
+% $$$ axis([0,3,0,1])
diff --git a/Tromp2005_kernel/process.sh b/init_plane/process.sh
similarity index 92%
copy from Tromp2005_kernel/process.sh
copy to init_plane/process.sh
index 08d1763..c1da134 100755
--- a/Tromp2005_kernel/process.sh
+++ b/init_plane/process.sh
@@ -18,12 +18,12 @@ echo
mkdir -p OUTPUT_FILES
mkdir -p DATA
-mkdir -p SEM
# sets up local DATA/ directory
cd DATA/
-ln -s ../Par_file_Tromp2005 Par_file
-ln -s ../SOURCE_Tromp2005 SOURCE
+rm -f Par_file SOURCE
+ln -s ../Par_file_Slave Par_file
+ln -s ../SOURCE_Slave SOURCE
cd ../
# cleans output files
diff --git a/init_plane/readme b/init_plane/readme
new file mode 100644
index 0000000..15baf5e
--- /dev/null
+++ b/init_plane/readme
@@ -0,0 +1,12 @@
+This directory gives an example of propgation of incoming plane wave into a layer over half space model.
+
+sequence of script:
+1) set_source .bash angle nrec
+
+2) process.sh 1 0 # forward
+
+3) matlab/select_adj.m, choose t1 and t2 (only once now)
+
+4) process.sh 0 1 # kernel
+
+5) matlab/plot_kernel.m (rho, vp and vs kernels)
diff --git a/init_plane/run.bash b/init_plane/run.bash
new file mode 100755
index 0000000..825bbbb
--- /dev/null
+++ b/init_plane/run.bash
@@ -0,0 +1,12 @@
+#!/bin/bash
+
+for dir in + -; do
+ for ang in 23 20 18; do
+ for nrec in 10 20 30; do
+ set_source.bash $dir $ang $nrec
+ run_for_ker.sh 1 1
+ mkdir -p ang-$ang; mkdir -p ang-$ang/nrec-$nrec-dir-$dir
+ mv OUTPUT_FILES_FORWARD OUTPUT_FILES_KERNEL ang-$ang/nrec-$nrec-dir-$dir
+ done
+ done
+done
diff --git a/init_plane/run_for_ker.sh b/init_plane/run_for_ker.sh
new file mode 100755
index 0000000..ab3d95e
--- /dev/null
+++ b/init_plane/run_for_ker.sh
@@ -0,0 +1,69 @@
+#!/bin/bash
+#
+# script runs mesher and solver (in serial)
+# using this example setup
+#
+
+if [ $# != 2 ]; then
+ echo 'process.sh forward kernel'; exit
+fi
+forward=$1
+kernel=$2
+echo ' '
+echo `date`
+if [ $forward -eq 1 ] ; then
+currentdir=`pwd`
+mkdir -p OUTPUT_FILES
+mkdir -p DATA
+rm -rf OUTPUT_FILES/* OUTPUT_FILES_FORWARD OUTPUT_FILES_KERNEL
+
+cd ../../
+make > tmp.log
+cd $currentdir
+rm -f xmeshfem2D xspecfem2D
+ln -s ../../bin/xmeshfem2D
+ln -s ../../bin/xspecfem2D
+
+# sets up local DATA/ directory
+cd DATA/
+rm -f Par_file SOURCE
+ln -s ../Par_file_Slave_for Par_file
+ln -s ../SOURCE_Slave SOURCE
+cd ../
+
+echo 'Forward simulation ...'
+echo " running mesher..."
+./xmeshfem2D > OUTPUT_FILES/output_mesher.txt
+echo " running solver..."
+./xspecfem2D > OUTPUT_FILES/output_solver.txt
+cp DATA/SOURCE_xz.dat DATA/STATIONS_* DATA/Par_file DATA/SOURCE OUTPUT_FILES/
+grep Receiver OUTPUT_FILES/output_mesher.txt | awk '{print $4,$5}' > OUTPUT_FILES/receiver.dat
+mv OUTPUT_FILES OUTPUT_FILES_FORWARD
+echo "see results in directory: OUTPUT_FILES_FORWARD/"
+fi
+
+
+if [ $kernel -eq 1 ]; then
+# run select_adj in matlab
+matlab.sh "select_adj"
+
+echo 'Kernel simulation ...'
+mkdir -p OUTPUT_FILES/
+cd SEM; rm -f *; cd ..
+cd OUTPUT_FILES_FORWARD; cp lastframe* absorb* ../OUTPUT_FILES; cp *.adj ../SEM; cd ..
+cd SEM; rename 's/semd.adj/adj/' *.adj; cd ..
+cd DATA; rm -f Par_file; ln -s ../Par_file_Slave_kernel Par_file; cd ..
+echo " running mesher..."
+./xmeshfem2D > OUTPUT_FILES/output_mesher.txt
+echo " running solver..."
+./xspecfem2D > OUTPUT_FILES/output_solver.txt
+cp DATA/Par_file OUTPUT_FILES/
+mv OUTPUT_FILES OUTPUT_FILES_KERNEL
+echo "see results in directory: OUTPUT_FILES_KERNEL/"
+
+matlab.sh "plot_kernel"
+fi
+
+
+echo "done"
+echo `date`
diff --git a/init_plane/select_adj.m b/init_plane/select_adj.m
new file mode 100644
index 0000000..dd53624
--- /dev/null
+++ b/init_plane/select_adj.m
@@ -0,0 +1,74 @@
+close all; clear all;
+% input parameters
+H=40; Vp1=6.5; Vp2=8.3; Vs1=3.7; Vs2=4.7; X0=-150;
+dir='OUTPUT_FILES_FORWARD/';
+
+% receivers
+Xs=100;Xe=500;
+nsta=30;
+if nsta == 1
+ xarray=[Xs];
+else
+ xarray=linspace(Xs,Xe,nsta);
+end
+
+% incident waves
+angle_inc=18;
+th=angle_inc*pi/180; p=sin(th)/Vp2;
+
+tmin=zeros(nsta,1); tmax=zeros(nsta,1);
+[tP, tPs, tPpp, tPps, tPss]=ttplane(X0,xarray,p,H,Vp2,Vs2,Vp1,Vs1);
+
+% reduced time
+rt=tP;
+% plot distance between traces
+dy=0.5;
+% start of seismograms
+
+%%% select t1 and t2 of adjoint source
+t1=2; t2=7; rtaper=0; % select Ps phase
+t1=2; t2=30;
+
+comp = ['BXZ'; 'BXX'];
+% plot seismograms
+for i = 1 : nsta
+ X=(Xe-Xs)/(nsta-1)*(i-1)+Xs;
+
+ for j = 1 : 2
+ figure(j);
+ file=sprintf([dir 'S00%02d.AA.' comp(j,:) '.semd'],i);
+ s=load(file);
+ ym = (i-1)*dy;
+ plot(s(:,1)-rt(i),s(:,2)+ym); hold on; grid on;
+
+ % pick start and end of adjoint seismograms (after aligning P)
+ tmin(i)=t1;
+ tmax(i)=min(t2,s(end,1)-rt(i));
+ plot(tmin(i),ym,'rx',tmax(i),ym,'mx');
+
+ % generate adjint source
+ ind=find(s(:,1)-rt(i)>tmin(i) & s(:,1)-rt(i)<tmax(i));
+ sadj=zeros(size(s(:,2)));
+ taper=tukeywin(length(ind),rtaper);
+ sadj(ind)=s(ind,2).*taper;
+ plot(s(ind,1)-rt(i),sadj(ind)+ym,'g');
+
+ % write adjoint source
+ fid = fopen(strcat(file,'.adj'),'w');
+ fprintf(fid,'%6.2f %12.4g\n',[s(:,1)';sadj']);
+ fclose(fid);
+ end
+
+ file=sprintf([dir 'S00%02d.AA.BXY.semd'],i);
+ fid = fopen(strcat(file,'.adj'),'w');
+ fprintf(fid,'%6.2f %12.4g\n',[s(:,1)';zeros(size(s(:,1)))']);
+ fclose(fid);
+end
+
+for j = 1 : 2
+ figure(j)
+ plot(tPs-tP,(0:nsta-1)*dy); text(tPs(1)-tP(1),-0.5,'tPs');
+ plot(tPss-tP,(0:nsta-1)*dy); text(tPss(1)-tP(1),-0.5,'tPss');
+ plot(tPps-tP,(0:nsta-1)*dy); text(tPps(1)-tP(1),-0.5,'tPps');
+ plot(tPpp-tP,(0:nsta-1)*dy); text(tPpp(1)-tP(1),-0.5,'tPpp');
+end
diff --git a/init_plane/set_source.bash b/init_plane/set_source.bash
new file mode 100755
index 0000000..fbd4c57
--- /dev/null
+++ b/init_plane/set_source.bash
@@ -0,0 +1,29 @@
+#!/bin/bash
+if [ $# != 3 ]; then
+ echo 'set_source +/- angle nrec'; exit
+fi
+sign=$1
+angle=$2
+nrec=$3
+
+sed -i "s/^angleforce .*$/angleforce = $sign$angle/" SOURCE_Slave
+if [ $sign = "+" ]; then
+ sed -i "s/^xs .*$/xs = -150000/" SOURCE_Slave
+ sed -i "s/^xdeb .*$/xdeb = 100000./" Par_file_Slave_for Par_file_Slave_kernel
+ sed -i "s/^xfin .*$/xfin = 500000./" Par_file_Slave_for Par_file_Slave_kernel
+ sed -i "s/^jpeg_header.*$/jpeg_header='left-kernel-';/" plot_kernel.m
+elif [ $sign = "-" ]; then
+ sed -i "s/^xs .*$/xs = 750000/" SOURCE_Slave
+ sed -i "s/^xdeb .*$/xdeb = 500000./" Par_file_Slave_for Par_file_Slave_kernel
+ sed -i "s/^xfin .*$/xfin = 100000./" Par_file_Slave_for Par_file_Slave_kernel
+ sed -i "s/^jpeg_header.*$/jpeg_header='right-kernel-';/" plot_kernel.m
+fi
+
+sed -i "s/^nrec .*$/nrec = $nrec/" Par_file_Slave_for
+sed -i "s/^nrec .*$/nrec = $nrec/" Par_file_Slave_kernel
+
+# set select_adj.m and plot_kernel.m
+# X0 is hard-wired into selected_adj.m
+sed -i "s/^angle_inc=.*$/angle_inc=$angle;/" select_adj.m
+sed -i "s/^nsta=.*$/nsta=$nrec;/" select_adj.m
+sed -i "s/^angle_inc=.*$/angle_inc=$angle;/" plot_kernel.m
diff --git a/init_plane/stack_kernel.m b/init_plane/stack_kernel.m
new file mode 100644
index 0000000..f790efa
--- /dev/null
+++ b/init_plane/stack_kernel.m
@@ -0,0 +1,58 @@
+close all; clear all;
+% choose incident angles and directions to stack
+LENGTH=250000; L1=LENGTH/1000; H=40;
+nrec=30;
+% modify directions of incident waves for stacked kernels
+dir='+-'; ndir=length(dir);
+angle=[18,20,23];
+
+for idir = 1 : ndir
+ cdir=dir(idir);
+ for j = 1 : length(angle)
+ cang=angle(j);
+
+ kdir=sprintf('ang-%2.2i/nrec-%2.2i-dir-%s/',cang,nrec,cdir)
+ file=[kdir 'OUTPUT_FILES_KERNEL/proc000000_rhop_alpha_beta_kernel.dat'];
+ s=load(file);
+ if idir == 1 & j == 1
+ rhop=s(:,3); alpha=s(:,4); beta=s(:,5);
+ file=[kdir 'OUTPUT_FILES_FORWARD/receiver.dat']; rec=load(file);
+ rec=rec/LENGTH; nrec=length(rec);
+ else
+ rhop=rhop+s(:,3); alpha=alpha+s(:,4); beta=beta+s(:,5);
+ end
+ end
+end
+
+xx=s(:,1)/LENGTH;zz=s(:,2)/LENGTH;
+minx=min(xx); maxx=max(xx); minz=min(zz); maxz=max(zz);
+nx=400; nz=400;
+% define mesh
+xa=linspace(minx,maxx,nx);za=linspace(minz,maxz,nz);
+[X,Z] = meshgrid(xa,za);
+
+subplot(3,1,1)
+FV=griddata(xx,zz,rhop,X,Z);
+imagesc(xa,za,FV); set(gca,'YDir','normal'); hold on;
+plot(rec(:,1),rec(:,2),'kx')
+text(0.2,0.2,'K_{\rho_p} ')
+quiver(0.2,1-H/L1, 0.1,0,'MaxHeadsize', 30)
+title(['stacked by inc angle & dir=',num2str(ndir), ', nrec=', num2str(nrec)]);
+
+subplot(3,1,2)
+FV=griddata(xx,zz,beta,X,Z);
+imagesc(xa,za,FV); set(gca,'YDir','normal'); hold on;
+plot(rec(:,1),rec(:,2),'kx')
+plot(rec(:,1),rec(:,2),'kx')
+text(0.2,0.2,'K_\beta')
+quiver(0.2,1-H/L1, 0.1,0,'MaxHeadsize', 30)
+
+subplot(3,1,3)
+FV=griddata(xx,zz,alpha,X,Z);
+imagesc(xa,za,FV); set(gca,'YDir','normal'); hold on;
+plot(rec(:,1),rec(:,2),'kx')
+text(0.2,0.2,'K_\alpha')
+quiver(0.2,1-H/L1, 0.1,0,'MaxHeadsize', 30)
+
+jpeg=['stacked-kernel-' num2str(nrec) '-' num2str(ndir) '.jpg']
+print('-djpeg',jpeg)
\ No newline at end of file
diff --git a/init_plane/ttplane.m b/init_plane/ttplane.m
new file mode 100644
index 0000000..f9d5aa6
--- /dev/null
+++ b/init_plane/ttplane.m
@@ -0,0 +1,50 @@
+function [tP, tPs, tPpp, tPps, tPss] = ttplane(x0, xarray, p, H, Vp2, Vs2, Vp1, Vs1)
+% this function computes the travel time of converted phases at
+% Moho in a layer (Vp1, Vs1) over half space (Vp2, Vs2) model.
+% [tP, tPs, tPpp, tPps, tPss] = ttplane(x0, xarray, p, H, Vp2, Vs2, Vp1, Vs1)
+% the incoming plane wave is defined as f((x-x_0)*p+z*eta-t), where
+% z is vertical pointing downwards and x is the horizontal direction.
+% x can be a vector; H is in km, V is in km/s and p in s/km
+%
+nx=length(xarray);
+thp1=asin(p*Vp1); thp2=asin(p*Vp2);
+ths1=asin(p*Vs1); ths2=asin(p*Vs2);
+etap1=cos(thp1)/Vp1; etap2=cos(thp2)/Vp2;
+etas1=cos(ths1)/Vs1; etas2=cos(ths2)/Vs2;
+% initialization
+tP=zeros(size(xarray)); tPs=tP; tPpp=tP; tPps=tP; tPss=tP;
+for i = 1: nx
+ x=xarray(i);
+ if (x0+H*cot(thp2) > x)
+ error('Can not handle incidence from the right')
+ end
+ tt=p*(x-x0-H*cot(thp2));
+ tP(i)=tt+etap1*H;
+ tPs(i)=tt+etas1*H;
+ tPpp(i)=tt+etap1*3*H;
+ tPss(i)=tt+etas1*2*H+etap1*H;
+ tPps(i)=tt+etap1*2*H+etas1*H;
+% $$$ tP(i)=p*(x-xs)+etap1*H+etap2*zs;
+% $$$ %Ps
+% $$$ xx=xx1-H*tan(ths1);
+% $$$ zs=xx*cos(thp2)*sin(thp2);
+% $$$ xs=x2+xx*cos(thp2)^2;
+% $$$ tPs(i)=p*(x-xs)+etas1*H+etap2*zs;
+% $$$ %Ppp
+% $$$ xx=xx1-2*H*tan(thp1);
+% $$$ zs=xx*cos(thp2)*sin(thp2);
+% $$$ xs=x2+xx*cos(thp2)^2;
+% $$$ tPpp(i)=p*(x-xs)+etap1*2*H+etap2*zs;
+% $$$ %Pps+Psp
+% $$$ xx=xx1-H*tan(thp1)-H*tan(ths1);
+% $$$ zs=xx*cos(thp2)*sin(thp2);
+% $$$ xs=x2+xx*cos(thp2)^2;
+% $$$ tPps(i)=p*(x-xs)+etap1*H+etas1*H+etap2*zs;
+% $$$ %Pss
+% $$$ xx=xx1-2*H*tan(ths1);
+% $$$ zs=xx*cos(thp2)*sin(thp2);
+% $$$ xs=x2+xx*cos(thp2)^2;
+% $$$ tPss(i)=p*(x-xs)+etas1*2*H+etap2*zs;
+end
+
+
\ No newline at end of file
More information about the CIG-COMMITS
mailing list