-
Joeri Brackenhoff authoredJoeri Brackenhoff authored
getParameters3D.c 52.15 KiB
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<assert.h>
#include"par.h"
#include"fdelmodc3D.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
/**
*
* The routine getParameters reads in all parameters to set up a FD modeling.
* Model and source parameters are used to calculate stability and dispersion relations
* Source and receiver positions are calculated and checked if they fit into the model.
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
float gaussGen();
long loptncr(long n);
long getModelInfo3D(char *file_name, long *n1, long *n2, long *n3,
float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
float *min, float *max, long *axis, long zeroch, long verbose);
long getWaveletInfo3D(char *file_src, long *n1, long *n2, float *d1, float *d2,
float *f1, float *f2, float *fmax, long *nxm, long verbose);
long getWaveletHeaders3D(char *file_src, long n1, long n2, float *gx, float *sx,
float *gy, float *sy, float *gelev, float *selev, long verbose);
long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
float dx, float dy, float dz, long nx, long ny, long nz);
long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *src,
shotPar *shot, bndPar *bnd, long verbose)
{
long isnapmax1, isnapmax2, isnapmax, sna_nrsna;
long n1, n2, n3, nx, ny, nz, nsrc, ix, axis, ioPz, is0, optn;
long npxsrc, npysrc, isx0, isy0, isx, isy;
long idzshot, idxshot, idyshot, nsrctext;
long src_ix0, src_iy0, src_iz0, src_ix1, src_iy1, src_iz1;
long disable_check;
float cp_min, cp_max, cs_min, cs_max, ro_min, ro_max;
float stabfactor,dispfactor, cmin, dt, fmax, scl, wfct, tapfact;
float zstart, xstart, ystart, d1, d2, d3, f1, f2, f3, sub_x0, sub_y0, sub_z0;
float srcendx, srcendy, srcendz, dx, dy, dz;
float xsrc, ysrc, zsrc, dxshot, dyshot, dzshot, dtshot;
float dxrcv, dyrcv, dzrcv, dxspread, dyspread, dzspread;
float tsnap1, tsnap2, dtsnap, dxsnap, dysnap, dzsnap, dtrcv;
float xsnap1, xsnap2, ysnap1, ysnap2, zsnap1, zsnap2, xmax, ymax, zmax;
float xsrc1, xsrc2, ysrc1, ysrc2, zsrc1, zsrc2, tsrc1, tsrc2, tlength, tactive;
float src_anglex, src_angley, src_velox, src_veloy, px, py, grad2rad, rdelay, scaledt;
float *xsrca, *ysrca, *zsrca, rrcv;
float Mxx, Myy, Mzz, Mxy, Myz, Mxz;
float rsrc, oxsrc, oysrc, ozsrc, dphisrc, ncsrc;
size_t nsamp;
long i, j, l;
long cfree;
long tapleft,tapright,taptop,tapbottom,tapfront, tapback;
long nxsrc, nysrc, nzsrc;
long largeSUfile;
long is,ntraces,length_random;
float rand;
char *src_positions, tmpname[1024]; char* src_txt;
FILE *fp;
if (!getparlong("verbose",&verbose)) verbose=0;
if (!getparlong("disable_check",&disable_check)) disable_check=0;
if (!getparlong("iorder",&mod->iorder)) mod->iorder=4;
if (!getparlong("ischeme",&mod->ischeme)) mod->ischeme=1;
if (!getparlong("sh",&mod->sh)) mod->sh=0;
if (!getparstring("file_cp",&mod->file_cp)) {
verr("parameter file_cp required!");
}
if (!getparstring("file_den",&mod->file_ro)) {
verr("parameter file_den required!");
}
if (mod->ischeme>2 && mod->ischeme!=5) {
if (!getparstring("file_cs",&mod->file_cs)) {
verr("parameter file_cs required!");
}
}
if (!getparstring("file_src",&wav->file_src)) wav->file_src=NULL;
if (!getparstring("file_snap",&sna->file_snap)) sna->file_snap="snap.su";
if (!getparstring("file_beam",&sna->file_beam)) sna->file_beam="beam.su";
if (!getparstring("file_rcv",&rec->file_rcv)) rec->file_rcv="recv.su";
if (!getparlong("grid_dir",&mod->grid_dir)) mod->grid_dir=0;
if (!getparlong("src_at_rcv",&src->src_at_rcv)) src->src_at_rcv=1;
/* read model parameters, which are used to set up source and receivers and check stability */
getModelInfo3D(mod->file_cp, &nz, &nx, &ny, &dz, &dx, &dy, &sub_z0, &sub_x0, &sub_y0, &cp_min, &cp_max, &axis, 1, verbose);
getModelInfo3D(mod->file_ro, &n1, &n2, &n3, &d1, &d2, &d3, &zstart, &xstart, &ystart, &ro_min, &ro_max, &axis, 0, verbose);
mod->cp_max = cp_max;
mod->cp_min = cp_min;
mod->ro_max = ro_max;
mod->ro_min = ro_min;
assert( (ro_min != 0.0) );
if (NINT(100*(dx/d2)) != 100)
vwarn("dx differs for file_cp and file_den!");
if (NINT(100*(dy/d3)) != 100)
vwarn("dy differs for file_cp and file_den!");
if (NINT(100*(dz/d1)) != 100)
vwarn("dz differs for file_cp and file_den!");
if (nx != n2)
vwarn("nx differs for file_cp and file_den!");
if (ny != n3)
vwarn("nx differs for file_cp and file_den!");
if (nz != n1)
vwarn("nz differs for file_cp and file_den!");
if (mod->ischeme>2 && mod->ischeme!=5) {
getModelInfo3D(mod->file_cs, &n1, &n2, &n3, &d1, &d2, &d3, &zstart, &xstart, &ystart, &cs_min, &cs_max, &axis, 1, verbose);
mod->cs_max = cs_max;
mod->cs_min = cs_min;
if (NINT(100*(dx/d2)) != 100)
vwarn("dx differs for file_cp and file_cs!");
if (NINT(100*(dy/d3)) != 100)
vwarn("dy differs for file_cp and file_cs!");
if (NINT(100*(dz/d1)) != 100)
vwarn("dz differs for file_cp and file_cs!");
if (nx != n2)
vwarn("nx differs for file_cp and file_cs!");
if (ny != n3)
vwarn("ny differs for file_cp and file_cs!");
if (nz != n1)
vwarn("nz differs for file_cp and file_cs!");
}
if (mod->ischeme==5) {
cs_max=0.0; cs_min=0.0;
mod->cs_max = cs_max; mod->cs_min = cs_min;
}
mod->nfz = nz;
mod->nfx = nx;
mod->nfy = ny;
/* check if 1D, 2D or full 3D gridded model is given as input model */
if (nx==1 && ny==1 ) { // 1D model
if (!getparlong("nx",&nx)) nx=nz;
if (!getparlong("ny",&ny)) ny=nx;
dx=dz;
dy=dx;
sub_x0=-0.5*(nx-1)*(dx);
sub_y0=-0.5*(ny-1)*(dy);
}
else if (ny==1) { // 2D model
if (!getparlong("ny",&ny)) ny=nx;
dy=dx;
sub_y0=-0.5*(ny-1)*(dy);
fprintf(stderr,"get param ny=%ld dy=%f y0=%f\n", ny, dy, sub_y0);
}
mod->dz = dz;
mod->dx = dx;
mod->dy = dy;
mod->nz = nz;
mod->nx = nx;
mod->ny = ny;
// end of part1 ################################################################################################
/* define wavelet(s), modeling time and wavelet maximum frequency */
if (wav->file_src!=NULL) {
getWaveletInfo3D(wav->file_src, &wav->ns, &wav->nx, &wav->ds, &d2, &f1, &f2, &fmax, &ntraces, verbose);
if (wav->ds <= 0.0) {
vwarn("dt in wavelet (file_src) equal to 0.0 or negative.");
vwarn("Use parameter dt= to overule dt from file_src.");
}
wav->nt = wav->ns;
wav->dt = wav->ds;
if(!getparfloat("tmod",&mod->tmod)) mod->tmod = (wav->nt-1)*wav->dt;
if(!getparfloat("dt",&mod->dt)) mod->dt=wav->dt;
if (NINT(wav->ds*1000000) != NINT(mod->dt*1000000)) {
if (wav->dt > mod->dt) {
scaledt = wav->dt/mod->dt;
scaledt = floorf(wav->dt/mod->dt);
optn = loptncr(wav->ns);
wav->nt = floorf(scaledt*optn);
vmess("file_src dt-scalefactor=%f : wav.dt=%e ==interpolated==> mod.dt=%e", scaledt, wav->dt, mod->dt);
wav->dt = mod->dt;
}
else {
wav->dt = mod->dt; /* in case if wav.dt is smaller than 1e-7 and can not be read by SU-getpar */
}
}
if(!getparfloat("fmax",&wav->fmax)) wav->fmax=fmax;
}
else {
fmax = 50;
if(!getparfloat("dt",&mod->dt)) verr("dt must be given or use file_src=");
if(!getparfloat("tmod",&mod->tmod)) verr("tmod must be given");
if(!getparfloat("fmax",&wav->fmax)) wav->fmax=fmax;
fmax = wav->fmax;
wav->dt=mod->dt;
}
assert(mod->dt!=0.0);
/* check if receiver delays is defined; option inactive: add delay time to total modeling time */
if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
rec->delay=NINT(rdelay/mod->dt); mod->nt = NINT(mod->tmod/mod->dt);
dt = mod->dt;
if (!getparlong("src_type",&src->type)) src->type=1;
if (!getparfloat("src_Mxx",&src->Mxx)) src->Mxx=0.0;
if (!getparlong("src_orient",&src->orient)) {
src->orient=1;
if (getparlong("dipsrc",&src->orient)) src->orient=2; // for compatability with DELPHI's fdacmod
}
if (mod->ischeme<=2) {
if (src->type>1 && src->type<6)
verr("Invalid src_type for acoustic scheme!");
}
if (mod->ischeme==2 || mod->ischeme==4) {
if (!getparstring("file_qp",&mod->file_qp)) mod->file_qp=NULL;
if (!getparstring("file_qs",&mod->file_qs)) mod->file_qs=NULL;
if (!getparfloat("Qp",&mod->Qp)) mod->Qp=1;
if (!getparfloat("Qs",&mod->Qs)) mod->Qs=mod->Qp;
if (!getparfloat("fw",&mod->fw)) mod->fw=0.5*wav->fmax;
}
/* dissipative medium option for Evert */
if (mod->ischeme==-1) {
if (!getparfloat("qr",&mod->qr)) mod->qr=0.1;
}
assert(src->type > 0);
/* dispersion factor to 10 points per wavelength (2nd order)
or 5 points per wavelength (4th order) */
if (mod->iorder == 2) {
dispfactor=10;
stabfactor = 1.0/sqrt(2.0);
}
else {
dispfactor = 6;
//stabfactor = 0.606; /* courant number */
stabfactor = 0.495; //Joeri check: number changes in 3D case. Sei&Simes, 1995.
/*However as we will see in the sequel, the stability condition
has only an indicative role. We will have to choose p = c. A t/Ax much less
than the maximum allowed by the stability condition to fulfill the precision
criteria we have imposed.*/
}
// end of part2 ################################################################################################
/* origin of model in real (non-grid) coordinates */
mod->x0 = sub_x0;
mod->y0 = sub_y0;
mod->z0 = sub_z0;
xmax = sub_x0+(nx-1)*dx;
ymax = sub_y0+(ny-1)*dy;
zmax = sub_z0+(nz-1)*dz;
if (verbose) {
vmess("*******************************************");
vmess("************** general info ***************");
vmess("*******************************************");
vmess("tmod = %f",mod->tmod);
vmess("ntsam = %li dt = %f(%e)",mod->nt, mod->dt, mod->dt);
if (mod->ischeme == 1) vmess("Acoustic staggered grid, pressure/velocity");
if (mod->ischeme == 2) vmess("Visco-Acoustic staggered grid, pressure/velocity");
if (mod->ischeme == 3) vmess("Elastic staggered grid, stress/velocity");
if (mod->ischeme == 4) vmess("Visco-Elastic staggered grid, stress/velocity");
if (mod->ischeme == 5) vmess("Acoustic staggered grid, Txx/Tzz/velocity");
if (mod->grid_dir) vmess("Time reversed modelling");
else vmess("Forward modelling");
vmess("*******************************************");
vmess("*************** model info ****************");
vmess("*******************************************"); if (mod->nfx == 1) {
vmess(" 1-dimensional model is read from file");
}
else if (mod->nfy == 1) {
vmess(" 2-dimensional model is read from file");
}
else {
vmess(" 3-dimensional model is read from file");
}
vmess("nz = %li ny = %li nx = %li", nz, ny, nx);
vmess("dz = %8.4f dy = %8.4f dx = %8.4f", dz, dy, dx);
vmess("zmin = %8.4f zmax = %8.4f", sub_z0, zmax);
vmess("ymin = %8.4f ymax = %8.4f", sub_y0, ymax);
vmess("xmin = %8.4f xmax = %8.4f", sub_x0, xmax);
vmess("min(cp) = %9.3f max(cp) = %9.3f", cp_min, cp_max);
if (mod->ischeme>2 && mod->ischeme!=5) vmess("min(cs) = %9.3f max(cs) = %9.3f", cs_min, cs_max);
vmess("min(ro) = %9.3f max(ro) = %9.3f", ro_min, ro_max);
if (mod->ischeme==2 || mod->ischeme==4) {
if (mod->file_qp!=NULL) vmess("Qp from file %s ", mod->file_qp);
else vmess("Qp = %9.3f ", mod->Qp);
vmess("at freq = %5.3f", mod->fw);
}
if (mod->ischeme==4) {
if (mod->file_qs!=NULL) vmess("Qs from file %s ", mod->file_qs);
else vmess("Qs = %9.3f ", mod->Qs);
vmess("at freq = %5.3f", mod->fw);
}
}
if (mod->ischeme <= 2) {
cmin = cp_min;
}
else {
cmin = cs_min;
if ( (cmin<1e-20) || (cp_min<cs_min) ) cmin=cp_min;
}
if (verbose) {
vmess("*******************************************");
vmess("******** dispersion and stability *********");
vmess("*******************************************");
vmess("Dispersion criterion is %3d points per wavelength: ", NINT(dispfactor));
vmess(" ====> wavelength > %f m [dx*disp]", dx*dispfactor);
vmess("The maximum frequency in source wavelet must be:");
vmess(" ====> frequency < %f Hz. [Cmin/dx*disp]", cmin/(dx*dispfactor));
vmess("Stability criterion for current settings: ");
vmess(" ====> Cp < %f m/s [dx*disp/dt]", dx*stabfactor/dt);
if (wav->file_src != NULL) vmess(" For wavelet(s) in file_src fmax = %f", fmax);
vmess("Optimal discretisation for current model:");
vmess(" With maximum velocity = %f dt <= %e", cp_max,dx*stabfactor/cp_max);
vmess(" With maximum frequency = %f dx <= %e", wav->fmax, cmin/(wav->fmax*dispfactor));
}
/* Check stability and dispersion setting */
if (cp_max > dx*stabfactor/dt) {
vwarn("************ ! Stability ! ****************");
vwarn("From the input file maximum P-wave velocity");
vwarn("in the current model is %f !!", cp_max);
vwarn("Hence, adjust dx >= %.4f,",cp_max*dt/stabfactor);
vwarn(" or adjust dt <= %f,",dx*stabfactor/cp_max);
vwarn(" or lower the maximum velocity below %.3f m/s.",dx*stabfactor/dt);
vwarn("***************** !!! *********************");
if (!disable_check) verr("********* leaving program *********");
}
if (wav->fmax > cmin/(dx*dispfactor)) {
vwarn("*********** ! Dispersion ! ****************");
vwarn("The maximum frequency in the source wavelet is");
vwarn("%.3f for stable modeling fmax < %.3f ", wav->fmax, cmin/(dx*dispfactor));
vwarn("Hence, adjust dx <= %.4f",cmin/(wav->fmax*dispfactor)); vwarn(" or adjust fmax <= %f (overruled with parameter fmax=),",cmin/(dx*dispfactor));
vwarn(" or increase the minimum velocity above %.3f m/s.",dx*dispfactor*wav->fmax);
vwarn("***************** !!! *********************");
if (!disable_check) verr("********* leaving program *********");
}
/* to support old parameter interface */
if (!getparlong("cfree",&cfree)) taptop=1;
if (!getparlong("tapleft",&tapleft)) tapleft=0;
if (!getparlong("tapright",&tapright)) tapright=0;
if (!getparlong("taptop",&taptop)) taptop=0;
if (!getparlong("tapbottom",&tapbottom)) tapbottom=0;
if (!getparlong("tapfront",&tapfront)) tapfront=0;
if (!getparlong("tapback",&tapback)) tapback=0;
if (tapleft) bnd->lef=4;
else bnd->lef=1;
if (tapright) bnd->rig=4;
else bnd->rig=1;
if (taptop) bnd->top=4;
else bnd->top=1;
if (tapbottom) bnd->bot=4;
else bnd->bot=1;
if (tapfront) bnd->fro=4;
else bnd->fro=1;
if (tapback) bnd->bac=4;
else bnd->bac=1;
/* define the type of boundaries */
/* 1=free 2=pml 3=rigid 4=taper */
if (!getparlong("left",&bnd->lef) && !tapleft) bnd->lef=4;
if (!getparlong("right",&bnd->rig)&& !tapright) bnd->rig=4;
if (!getparlong("top",&bnd->top) && !taptop) bnd->top=1;
if (!getparlong("bottom",&bnd->bot) && !tapbottom) bnd->bot=4;
if (!getparlong("front",&bnd->fro) && !tapfront) bnd->fro=4;
if (!getparlong("back",&bnd->bac) && !tapback) bnd->bac=4;
/* calculate default taper length to be three wavelenghts */
if (!getparlong("ntaper",&bnd->ntap)) bnd->ntap=0; // bnd->ntap=3*NINT((cp_max/wav->fmax)/dx);
if (!bnd->ntap) if (!getparlong("npml",&bnd->ntap)) bnd->ntap=3*NINT((cp_max/wav->fmax)/dx);
if (!getparfloat("R",&bnd->R)) bnd->R=1e-5;
if (!getparfloat("m",&bnd->m)) bnd->m=2.0;
bnd->npml=bnd->ntap;
if (bnd->ntap) {
bnd->tapx = (float *)malloc(bnd->ntap*sizeof(float));
bnd->tapy = (float *)malloc(bnd->ntap*sizeof(float));
bnd->tapz = (float *)malloc(bnd->ntap*sizeof(float));
bnd->tapxz = (float *)malloc(bnd->ntap*bnd->ntap*sizeof(float));
bnd->tapxyz = (float *)malloc(bnd->ntap*bnd->ntap*bnd->ntap*sizeof(float));
if(!getparfloat("tapfact",&tapfact)) tapfact=0.30;
scl = tapfact/((float)bnd->ntap);
for (i=0; i<bnd->ntap; i++) {
wfct = (scl*i);
bnd->tapx[i] = exp(-(wfct*wfct));
bnd->tapy[i] = exp(-(wfct*wfct));
wfct = (scl*(i+0.5));
bnd->tapz[i] = exp(-(wfct*wfct));
}
//corner
for (j=0; j<bnd->ntap; j++) {
for (i=0; i<bnd->ntap; i++) {
wfct = (scl*sqrt(i*i+j*j));
bnd->tapxz[j*bnd->ntap+i] = exp(-(wfct*wfct));
}
}
//double corner
for (j=0; j<bnd->ntap; j++) {
for (i=0; i<bnd->ntap; i++) {
for (l=0; l<bnd->ntap; l++) {
wfct = (scl*sqrt(i*i+j*j+l*l));
bnd->tapxyz[l*(bnd->ntap)*(bnd->ntap)+j*(bnd->ntap)+i] = exp(-(wfct*wfct));
}
}
}
}
/* Vx: rox */
mod->ioXx=mod->iorder/2;
mod->ioXy=mod->iorder/2-1;
mod->ioXz=mod->iorder/2-1;
/* Vy: roy */
mod->ioYx=mod->iorder/2-1;
mod->ioYy=mod->iorder/2;
mod->ioYz=mod->iorder/2-1;
/* Vz: roz */
mod->ioZx=mod->iorder/2-1;
mod->ioZy=mod->iorder/2-1;
mod->ioZz=mod->iorder/2;
/* P, Txx, Tzz: lam, l2m */
mod->ioPx=mod->iorder/2-1;
mod->ioPy=mod->ioPx;
mod->ioPz=mod->ioPx;
/* Txz: mul */
mod->ioTx=mod->iorder/2;
mod->ioTy=mod->ioTx;
mod->ioTz=mod->ioTx;
/* end loop iteration in FD kernels */
/* Vx: rox */
mod->ieXx=nx+mod->ioXx;
mod->ieXy=ny+mod->ioXy;
mod->ieXz=nz+mod->ioXz;
/* Vy: roy */
mod->ieYx=nx+mod->ioYx;
mod->ieYy=ny+mod->ioYy;
mod->ieYz=nz+mod->ioYz;
/* Vz: roz */
mod->ieZx=nx+mod->ioZx;
mod->ieZy=ny+mod->ioZy;
mod->ieZz=nz+mod->ioZz;
/* P, Txx, Tzz: lam, l2m */
mod->iePx=nx+mod->ioPx;
mod->iePy=ny+mod->ioPy;
mod->iePz=nz+mod->ioPz;
/* Txz: muu */
mod->ieTx=nx+mod->ioTx;
mod->ieTy=ny+mod->ioTy;
mod->ieTz=nz+mod->ioTz;
mod->naz = mod->nz+mod->iorder;
mod->nay = mod->ny+mod->iorder;
mod->nax = mod->nx+mod->iorder;
/* for tapered and PML extra points are needed at the boundaries of the model */
if (bnd->top==4 || bnd->top==2) {
mod->naz += bnd->ntap;
mod->ioXz += bnd->ntap;
mod->ioYz += bnd->ntap;
mod->ioZz += bnd->ntap;
mod->ieXz += bnd->ntap;
mod->ieYz += bnd->ntap;
mod->ieZz += bnd->ntap;
/* For P/Tzz, Txx and Txz fields the tapered boundaries are calculated in the main kernels */
mod->iePz += bnd->ntap;
mod->ieTz += bnd->ntap; }
if (bnd->bot==4 || bnd->bot==2) {
mod->naz += bnd->ntap;
/* For P/Tzz, Txx and Txz fields the tapered boundaries are calculated in the main kernels */
mod->iePz += bnd->ntap;
mod->ieTz += bnd->ntap;
}
if (bnd->lef==4 || bnd->lef==2) {
mod->nax += bnd->ntap;
mod->ioXx += bnd->ntap;
mod->ioYx += bnd->ntap;
mod->ioZx += bnd->ntap;
mod->ieXx += bnd->ntap;
mod->ieYx += bnd->ntap;
mod->ieZx += bnd->ntap;
/* For Tzz, Txx and Txz fields the tapered boundaries are calculated in the main kernels */
mod->iePx += bnd->ntap;
mod->ieTx += bnd->ntap;
}
if (bnd->rig==4 || bnd->rig==2) {
mod->nax += bnd->ntap;
/* For P/Tzz, Txx and Txz fields the tapered boundaries are calculated in the main kernels */
mod->iePx += bnd->ntap;
mod->ieTx += bnd->ntap;
}
if (bnd->fro==4 || bnd->fro==2) {
mod->nay += bnd->ntap;
mod->ioXy += bnd->ntap;
mod->ioYy += bnd->ntap;
mod->ioZy += bnd->ntap;
mod->ieXy += bnd->ntap;
mod->ieYy += bnd->ntap;
mod->ieZy += bnd->ntap;
/* For Tzz, Txx and Txz fields the tapered boundaries are calculated in the main kernels */
mod->iePy += bnd->ntap;
mod->ieTy += bnd->ntap;
}
if (bnd->bac==4 || bnd->bac==2) {
mod->nay += bnd->ntap;
/* For P/Tzz, Txx and Txz fields the tapered boundaries are calculated in the main kernels */
mod->iePy += bnd->ntap;
mod->ieTy += bnd->ntap;
}
/* Intialize the array which contains the topography surface */
if (bnd->top==4 || bnd->top==2) ioPz=mod->ioPz - bnd->ntap;
else ioPz=mod->ioPz;
ioPz=mod->ioPz;
bnd->surface = (long *)malloc((mod->nax*mod->nay+mod->naz)*sizeof(long));
for (ix=0; ix<mod->nax*mod->nay+mod->naz; ix++) {
bnd->surface[ix] = ioPz;
}
if (verbose) {
vmess("*******************************************");
vmess("************* boundary info ***************");
vmess("*******************************************");
vmess("*** 1=free 2=pml 3=rigid 4=tapered ***");
vmess("Top boundary : %li",bnd->top);
vmess("Left boundary : %li",bnd->lef);
vmess("Right boundary : %li",bnd->rig);
vmess("Bottom boundary : %li",bnd->bot);
vmess("Front boundary : %li",bnd->fro);
vmess("Back boundary : %li",bnd->bac);
vmess("taper length = %li points",bnd->ntap);
}
/* define the number and type of shots to model */
/* each shot can have multiple sources arranged in different ways */
if (!getparfloat("xsrc",&xsrc)) xsrc=sub_x0+((nx-1)*dx)/2.0;
if (!getparfloat("ysrc",&ysrc)) ysrc=sub_y0+((ny-1)*dy)/2.0;
if (!getparfloat("zsrc",&zsrc)) zsrc=sub_z0;
if (!getparlong("nshot",&shot->n)) shot->n=1;
if (!getparfloat("dxshot",&dxshot)) dxshot=dx;
if (!getparfloat("dyshot",&dyshot)) dyshot=dy;
if (!getparfloat("dzshot",&dzshot)) dzshot=0.0;
if (!getparfloat("dip",&src->dip)) src->dip=0.0;
if (!getparfloat("strike",&src->strike)) src->strike=0.0;
if (!getparfloat("rake",&src->rake)) src->rake=0.0;
src->dip = M_PI*(src->dip/180.0);
src->strike = M_PI*(src->strike/180.0);
src->rake = M_PI*(src->rake/180.0);
if (shot->n>1) {
idxshot=MAX(0,NINT(dxshot/dx));
idyshot=MAX(0,NINT(dyshot/dy));
idzshot=MAX(0,NINT(dzshot/dz));
}
else {
idxshot=0.0;
idyshot=0.0;
idzshot=0.0;
}
/* calculate the shot positions */
src_ix0=MAX(0,NINT((xsrc-sub_x0)/dx));
src_ix0=MIN(src_ix0,nx);
src_iy0=MAX(0,NINT((ysrc-sub_y0)/dy));
src_iy0=MIN(src_iy0,ny);
src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
src_iz0=MIN(src_iz0,nz);
srcendx=(shot->n-1)*dxshot+xsrc;
srcendy=(shot->n-1)*dyshot+ysrc;
srcendz=(shot->n-1)*dzshot+zsrc;
src_ix1=MAX(0,NINT((srcendx-sub_x0)/dx));
src_ix1=MIN(src_ix1,nx);
src_iy1=MAX(0,NINT((srcendy-sub_y0)/dy));
src_iy1=MIN(src_iy1,ny);
src_iz1=MAX(0,NINT((srcendz-sub_z0)/dz));
src_iz1=MIN(src_iz1,nz);
shot->x = (long *)calloc(shot->n,sizeof(long));
shot->y = (long *)calloc(shot->n,sizeof(long));
shot->z = (long *)calloc(shot->n,sizeof(long));
for (is=0; is<shot->n; is++) {
shot->x[is] = src_ix0+is*idxshot;
shot->y[is] = src_iy0+is*idyshot;
shot->z[is] = src_iz0+is*idzshot;
if (shot->x[is] > nx-1) shot->n = is-1;
if (shot->y[is] > ny-1) shot->n = is-1;
if (shot->z[is] > nz-1) shot->n = is-1;
}
/* check if source array is defined */
nxsrc = countparval("xsrca");
nysrc = countparval("ysrca");
nzsrc = countparval("zsrca");
if (nxsrc != nzsrc) {
verr("Number of sources in array xsrca (%li), ysrca (%li), zsrca (%li) are not equal",nxsrc, nysrc, nzsrc);
}
/* source positions defined through txt file */
if (!getparstring("src_txt",&src_txt)) src_txt=NULL;
/* check if sources on a circle are defined */
if (getparfloat("rsrc", &rsrc)) {
if (!getparfloat("dphisrc",&dphisrc)) dphisrc=2.0;
if (!getparfloat("oxsrc",&oxsrc)) oxsrc=0.0; if (!getparfloat("oysrc",&oysrc)) oysrc=0.0;
if (!getparfloat("ozsrc",&ozsrc)) ozsrc=0.0;
ncsrc = NINT(360.0/dphisrc);
src->n = nsrc;
src->x = (long *)malloc(ncsrc*sizeof(long));
src->y = (long *)malloc(ncsrc*sizeof(long));
src->z = (long *)malloc(ncsrc*sizeof(long));
for (ix=0; ix<ncsrc; ix++) {
src->x[ix] = NINT((oxsrc-sub_x0+rsrc*cos(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dx);
src->y[ix] = NINT((oysrc-sub_y0+rsrc*sin(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dy);
src->z[ix] = NINT((ozsrc-sub_z0+rsrc*sin(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dz);
if (verbose>4) fprintf(stderr,"Source on Circle: xsrc[%li]=%li ysrc=%li zsrc=%li\n", ix, src->x[ix], src->y[ix], src->z[ix]);
}
}
/* TO DO propagate src_positions parameter and structure through code */
if (!getparstring("src_positions",&src_positions)) src_positions="single";
wav->random=0;
src->random=0;
src->plane=0;
src->array=0;
src->single=0;
if (strstr(src_positions, "single")) src->single=1;
else if (strstr(src_positions, "array")) src->array=1;
else if (strstr(src_positions, "random")) src->random=1;
else if (strstr(src_positions, "plane")) src->plane=1;
else src->single=1;
/* to maintain functionality of older parameters usage */
if (!getparlong("src_random",&src->random)) src->random=0;
if (!getparlong("plane_wave",&src->plane)) src->plane=0;
if (src->random) {
if (!getparlong("wav_random",&wav->random)) wav->random=1;
src->plane=0;
src->array=0;
src->single=0;
}
else {
if (!getparlong("wav_random",&wav->random)) wav->random=0;
}
if (src->plane) {
src->random=0;
src->array=0;
src->single=0;
}
if (!wav->random) assert (wav->file_src != NULL);
if (wav->random) {
wav->nt=mod->nt;
wav->dt=mod->dt;
wav->nx=1;
}
/* number of sources per shot modeling */
if (!getparlong("src_nxwindow",&src->nxwindow)) src->nxwindow=0;
if (!getparlong("src_nywindow",&src->nywindow)) src->nywindow=0;
src->window = 0;
if (!getparfloat("src_anglex",&src_anglex)) src_anglex=0.;
if (!getparfloat("src_angley",&src_angley)) src_angley=0.;
if (!getparfloat("src_velox",&src_velox)) src_velox=1500.;
if (!getparfloat("src_veloy",&src_veloy)) src_veloy=1500.;
if (!getparlong("distribution",&src->distribution)) src->distribution=0;
if (!getparlong("src_multiwav",&src->multiwav)) src->multiwav=0;
if (!getparfloat("amplitude", &src->amplitude)) src->amplitude=0.0;
if (!getparfloat("tlength", &tlength)) tlength=mod->dt*(mod->nt-1); if (!getparlong("src_injectionrate", &src->injectionrate)) src->injectionrate=0;
if (src->random && nxsrc==0) {
if (!getparlong("nsrc",&nsrc)) nsrc=1;
if (!getparlong("seed",&wav->seed)) wav->seed=10;
if (!getparfloat("xsrc1", &xsrc1)) xsrc1=sub_x0;
if (!getparfloat("xsrc2", &xsrc2)) xsrc2=xmax;
if (!getparfloat("ysrc1", &ysrc1)) ysrc1=sub_y0;
if (!getparfloat("ysrc2", &ysrc2)) ysrc2=ymax;
if (!getparfloat("zsrc1", &zsrc1)) zsrc1=sub_z0;
if (!getparfloat("zsrc2", &zsrc2)) zsrc2=zmax;
if (!getparfloat("tsrc1", &tsrc1)) tsrc1=0.0;
if (!getparfloat("tsrc2", &tsrc2)) tsrc2=mod->tmod;
if (!getparfloat("tactive", &tactive)) tactive=tsrc2;
tsrc2 = MIN(tsrc2, mod->tmod);
if (!getparfloat("tlength", &tlength)) tlength=tsrc2-tsrc1;
if (!getparlong("length_random", &length_random)) length_random=1;
dxshot = xsrc2-xsrc1;
dyshot = ysrc2-ysrc1;
dzshot = zsrc2-zsrc1;
dtshot = tsrc2-tsrc1;
if (wav->random) {
if (!getparlong("src_multiwav",&src->multiwav)) src->multiwav=1;
if (src->multiwav) wav->nx = nsrc;
else wav->nx = 1;
}
if (wav->random) wav->nt = NINT(tlength/mod->dt)+1;
src->tbeg = (float *)malloc(nsrc*sizeof(float));
src->tend = (float *)malloc(nsrc*sizeof(float));
wav->nsamp = (size_t *)malloc((nsrc+1)*sizeof(size_t));
src->x = (long *)malloc(nsrc*sizeof(long));
src->y = (long *)malloc(nsrc*sizeof(long));
src->z = (long *)malloc(nsrc*sizeof(long));
nsamp = 0;
srand48(wav->seed);
for (is=0; is<nsrc; is++) {
rand = (float)drand48();
src->x[is] = NINT((xsrc1+rand*dxshot-sub_x0)/dx);
rand = (float)drand48();
src->y[is] = NINT((ysrc1+rand*dyshot-sub_y0)/dy);
rand = (float)drand48();
src->z[is] = NINT((zsrc1+rand*dzshot-sub_z0)/dz);
if (length_random) rand = (float)drand48();
else rand = 0.0;
src->tbeg[is] = tsrc1+rand*(dtshot);
if (wav->random) {
if (src->distribution) rand = fabsf(tlength+gaussGen()*tlength);
else rand = (float)drand48()*tlength;
if (length_random!=1) rand = tlength;
src->tend[is] = MIN(src->tbeg[is]+rand, tactive);
wav->nsamp[is] = (size_t)(NINT((src->tend[is]-src->tbeg[is])/mod->dt)+1);
}
else {
src->tend[is] = MIN(src->tbeg[is]+(wav->nt-1)*wav->dt,mod->tmod);
wav->nsamp[is] = wav->nt;
}
nsamp += wav->nsamp[is];
if (verbose>3) {
vmess("Random xsrc=%f ysrc=%f zsrc=%f src_tbeg=%f src_tend=%f nsamp=%ld",src->x[is]*dx, src->y[is]*dy, src->z[is]*dz, src->tbeg[is], src->tend[is], wav->nsamp[is]);
}
}
wav->nsamp[nsrc] = nsamp; /* put total number of samples in last position */
wav->nst = nsamp; /* put total number of samples in nst part */
/* write time and length of source signals */
if (verbose>3) {
float *dum;
dum = (float *)calloc(mod->nt, sizeof(float));
for (is=0; is<nsrc; is++) {
dum[(long)floor(src->tbeg[is]/mod->dt)] = src->tend[is]-src->tbeg[is];
} FILE *fp;
sprintf(tmpname,"srcTimeLengthN=%li.bin",mod->nt);
fp = fopen(tmpname, "w+");
fwrite(dum, sizeof(float), mod->nt, fp);
fclose(fp);
free(dum);
}
}
else if ( (nxsrc != 0) || (src_txt != NULL) ) {
/* source array is defined */
if (src_txt!=NULL) {
/* Sources from a Text File */
/* Open text file */
nsrctext=0;
fp=fopen(src_txt,"r");
assert(fp!=NULL);
/* Get number of lines */
while (!feof(fp)) if (fgetc(fp)=='\n') nsrctext++;
fseek(fp,-1,SEEK_CUR);
if (fgetc(fp)!='\n') nsrctext++; /* Checks if last line terminated by /n */
if (verbose) vmess("Number of sources in src_txt file: %li",nsrctext);
rewind(fp);
nsrc=nsrctext;
}
else {
nsrc=nxsrc;
}
/* Allocate arrays */
src->x = (long *)malloc(nsrc*sizeof(long));
src->y = (long *)malloc(nsrc*sizeof(long));
src->z = (long *)malloc(nsrc*sizeof(long));
src->tbeg = (float *)malloc(nsrc*sizeof(float));
src->tend = (float *)malloc(nsrc*sizeof(float));
xsrca = (float *)malloc(nsrc*sizeof(float));
ysrca = (float *)malloc(nsrc*sizeof(float));
zsrca = (float *)malloc(nsrc*sizeof(float));
if (src_txt!=NULL) {
/* Read in source coordinates */
for (i=0;i<nsrc;i++) {
if (fscanf(fp,"%e %e %e\n",&xsrca[i],&ysrca[i],&zsrca[i])!=3) vmess("Source Text File: Can not parse coordinates on line %li.",i);
}
/* Close file */
fclose(fp);
}
else {
getparfloat("xsrca", xsrca);
getparfloat("ysrca", ysrca);
getparfloat("zsrca", zsrca);
}
/* Process coordinates */
for (is=0; is<nsrc; is++) {
src->x[is] = NINT((xsrca[is]-sub_x0)/dx);
src->y[is] = NINT((ysrca[is]-sub_y0)/dy);
src->z[is] = NINT((zsrca[is]-sub_z0)/dz);
src->tbeg[is] = 0.0;
src->tend[is] = (wav->nt-1)*wav->dt;
if (verbose>3) fprintf(stderr,"Source Array: xsrc[%li]=%f ysrc=%f zsrc=%f\n", is, xsrca[is], ysrca[is], zsrca[is]);
}
src->random = 1;
wav->nsamp = (size_t *)malloc((nsrc+1)*sizeof(size_t));
if (wav->random) {
if (!getparlong("src_multiwav",&src->multiwav)) src->multiwav=1;
if (src->multiwav) wav->nx = nsrc;
else wav->nx = 1;
wav->nt = NINT(tlength/mod->dt)+1;
nsamp=0;
for (is=0; is<nsrc; is++) {
rand = (float)drand48()*tlength; src->tend[is] = MIN(src->tbeg[is]+rand, mod->tmod);
wav->nsamp[is] = (size_t)(NINT((src->tend[is]-src->tbeg[is])/mod->dt)+1);
nsamp += wav->nsamp[is];
}
wav->nsamp[nsrc] = nsamp; /* put total number of samples in last position */
wav->nst = nsamp; /* put total number of samples in nst part */
}
else {
nsamp=0;
for (is=0; is<nsrc; is++) {
wav->nsamp[is] = wav->nt;
nsamp += wav->nsamp[is];
}
wav->nsamp[nsrc] = nsamp; /* put total number of samples in last position */
wav->nst = nsamp; /* put total number of samples in nst part */
}
free(xsrca);
free(ysrca);
free(zsrca);
}
else if (wav->nx > 1) {
/* read file_src for number of sources and receiver positions */
if (!getparlong("src_multiwav",&src->multiwav)) src->multiwav=1;
float *gx, *sx, *gy, *sy, *gelev, *selev;
gx = (float *)malloc(wav->nx*sizeof(float));
sx = (float *)malloc(wav->nx*sizeof(float));
gy = (float *)malloc(wav->nx*sizeof(float));
sy = (float *)malloc(wav->nx*sizeof(float));
gelev = (float *)malloc(wav->nx*sizeof(float));
selev = (float *)malloc(wav->nx*sizeof(float));
getWaveletHeaders3D(wav->file_src, wav->ns, wav->nx, gx, sx, gy, sy, gelev, selev, verbose);
nsrc = wav->nx;
src->x = (long *)malloc(nsrc*sizeof(long));
src->y = (long *)malloc(nsrc*sizeof(long));
src->z = (long *)malloc(nsrc*sizeof(long));
src->tbeg = (float *)malloc(nsrc*sizeof(float));
src->tend = (float *)malloc(nsrc*sizeof(float));
wav->nsamp = (size_t *)malloc((nsrc+1)*sizeof(size_t));
nsamp=0;
for (is=0; is<nsrc; is++) {
if (src->src_at_rcv>0){
src->x[is] = NINT((gx[is]-sub_x0)/dx);
src->y[is] = NINT((gy[is]-sub_y0)/dy);
src->z[is] = NINT((gelev[is]-sub_z0)/dz);
if (verbose>3) fprintf(stderr,"Source Array: xsrc[%li]=%f %li ysrc=%f %li zsrc=%f %li\n", is, gx[is], src->x[is], gy[is], src->y[is], gelev[is], src->z[is]);
}
else {
src->x[is]=NINT((sx[is]-sub_x0)/dx);
src->y[is]=NINT((sy[is]-sub_y0)/dy);
src->z[is]=NINT((selev[is]-sub_z0)/dz);
if (verbose>3) fprintf(stderr,"Source Array: xsrc[%li]=%f %li ysrc=%f %li zsrc=%f %li\n", is, sx[is], src->x[is], sy[is], src->y[is], selev[is], src->z[is]);
}
src->tbeg[is] = 0.0;
src->tend[is] = (wav->nt-1)*wav->dt;
wav->nsamp[is] = (size_t)(NINT((src->tend[is]-src->tbeg[is])/mod->dt)+1);
nsamp += wav->nsamp[is];
}
wav->nsamp[nsrc] = nsamp; /* put total number of samples in last position */
free(gx);
free(sx);
free(gy);
free(sy);
free(gelev);
free(selev);
}
else {
if (src->plane) {
if (!getparlong("npxsrc",&npxsrc)) npxsrc=1;
if (!getparlong("npysrc",&npysrc)) npysrc=1;
} else {
npxsrc=1;
npysrc=1;
}
if (npxsrc > nx) {
vwarn("Number of sources used in plane wave (%li) is larger than ",npxsrc);
vwarn("number of gridpoints in X (%li). Plane wave will be clipped to the edges of the model",nx);
npxsrc = mod->nx;
}
if (npysrc > ny) {
vwarn("Number of sources used in plane wave (%li) is larger than ",npysrc);
vwarn("number of gridpoints in Y (%li). Plane wave will be clipped to the edges of the model",ny);
npysrc = mod->ny;
}
nsrc = npxsrc*npysrc;
/* for a source defined on mutliple gridpoint calculate p delay factor */
src->x = (long *)malloc(nsrc*sizeof(long));
src->y = (long *)malloc(nsrc*sizeof(long));
src->z = (long *)malloc(nsrc*sizeof(long));
src->tbeg = (float *)malloc(nsrc*sizeof(float));
src->tend = (float *)malloc(nsrc*sizeof(float));
grad2rad = 17.453292e-3;
px = sin(src_anglex*grad2rad)/src_velox;
py = sin(src_angley*grad2rad)/src_veloy;
if (py < 0.0) {
for (isy=0; isy<npysrc; isy++) {
if (px < 0.0) {
for (isx=0; isx<npxsrc; isx++) {
src->tbeg[isy*npxsrc+isx] = fabsf((npysrc-isy-1)*dy*py) + fabsf((npxsrc-isx-1)*dx*px);
}
}
else {
for (isx=0; isx<npxsrc; isx++) {
src->tbeg[isy*npxsrc+isx] = fabsf((npysrc-isy-1)*dy*py) + isx*dx*px;
}
}
}
}
else {
for (isy=0; isy<npysrc; isy++) {
if (px < 0.0) {
for (isx=0; isx<npxsrc; isx++) {
src->tbeg[isy*npxsrc+isx] = isy*dy*py + fabsf((npxsrc-isx-1)*dx*px);
}
}
else {
for (isx=0; isx<npxsrc; isx++) {
src->tbeg[isy*npxsrc+isx] = isy*dy*py + isx*dx*px;
}
}
}
}
for (is=0; is<nsrc; is++) {
src->tend[is] = src->tbeg[is] + (wav->nt-1)*wav->dt;
}
isx0 = -1*floor((npxsrc-1)/2);
isy0 = -1*floor((npysrc-1)/2);
for (isy=0; isy<npysrc; isy++) {
for (isx=0; isx<npxsrc; isx++) {
src->x[isy*npxsrc+isx] = isx0 + isx;
src->y[isy*npxsrc+isx] = isy0 + isy;
src->z[isy*npxsrc+isx] = 0;
}
}
if (wav->random) {
if (!getparlong("src_multiwav",&src->multiwav)) src->multiwav=1; if (src->multiwav) wav->nx = nsrc;
else wav->nx = 1;
wav->nt = NINT(tlength/mod->dt)+1;
wav->nsamp = (size_t *)malloc((wav->nx+1)*sizeof(size_t));
nsamp=0;
for (is=0; is<wav->nx; is++) {
rand = (float)drand48()*tlength;
src->tend[is] = MIN(src->tbeg[is]+rand, mod->tmod);
wav->nsamp[is] = (size_t)(NINT((src->tend[is]-src->tbeg[is])/mod->dt)+1);
nsamp += wav->nsamp[is];
}
wav->nsamp[nsrc] = nsamp; /* put total number of samples in last position */
wav->nst = nsamp; /* put total number of samples in nst part */
}
else {
wav->nsamp = (size_t *)malloc((nsrc+1)*sizeof(size_t));
nsamp=0;
for (is=0; is<nsrc; is++) {
wav->nsamp[is] = wav->nt;
nsamp += wav->nsamp[is];
}
wav->nsamp[nsrc] = nsamp; /* put total number of samples in last position */
wav->nst = nsamp; /* put total number of samples in nst part */
}
src->nx = npxsrc;
src->ny = npysrc;
}
if (src->multiwav) {
if (wav->nx != nsrc) {
vwarn("src_multiwav has been defined but number of traces in");
vwarn("file_src = %li is not equal to nsrc = %li", wav->nx, nsrc);
vwarn("last trace in file_src will be repeated.");
}
else {
if (wav->file_src != NULL) vmess("Using all traces in file_src for a real shot");
}
}
src->n=nsrc;
if (verbose) {
vmess("*******************************************");
vmess("************* wavelet info ****************");
vmess("*******************************************");
vmess("wav_nt = %6li wav_nx = %li", wav->ns, wav->nx);
vmess("src_type = %6li src_orient = %li", src->type, src->orient);
vmess("number of sources = %li", shot->n);
vmess("fmax = %8.2f", fmax);
fprintf(stderr," %s: Source type : ",xargv[0]);
switch ( src->type ) {
case 1 : fprintf(stderr,"P "); break;
case 2 : fprintf(stderr,"Txz "); break;
case 3 : fprintf(stderr,"Tzz "); break;
case 4 : fprintf(stderr,"Txx "); break;
case 5 : fprintf(stderr,"S-potential"); break;
case 6 : fprintf(stderr,"Fx "); break;
case 7 : fprintf(stderr,"Fz "); break;
case 8 : fprintf(stderr,"P-potential"); break;
case 9 : fprintf(stderr,"Double-couple"); break;
case 10 : fprintf(stderr,"Moment tensor"); break;
}
fprintf(stderr,"\n");
if (src->type==9) {
Mxx = -1.0*(sin(src->dip)*cos(src->rake)*sin(2.0*src->strike)+sin(src->dip*2.0)*sin(src->rake)*sin(src->strike)*sin(src->strike));
Myy = sin(src->dip)*cos(src->rake)*sin(2.0*src->strike)-sin(src->dip*2.0)*sin(src->rake)*cos(src->strike)*cos(src->strike);
Mzz = sin(src->dip*2.0)*sin(src->rake);
Mxz = -1.0*(cos(src->dip)*cos(src->rake)*cos(src->strike)+cos(src->dip*2.0)*sin(src->rake)*sin(src->strike));
Mxy = sin(src->dip)*cos(src->rake)*cos(src->strike*2.0)+0.5*(sin(src->dip*2.0)*sin(src->rake)*sin(src->strike*2.0));
Myz = -1.0*(cos(src->dip)*cos(src->rake)*sin(src->strike)-cos(src->dip*2.0)*sin(src->rake)*cos(src->strike)); vmess("Strike %.3f (%.2f degrees) Rake %.3f (%.2f degrees) Dip %.3f (%.2f degrees)",src->strike,180.0*src->strike/M_PI,src->rake,180.0*src->rake/M_PI,src->dip,180.0*src->dip/M_PI);
vmess("Mxx %.2f Myy %.2f Mzz %.2f",Mxx,Myy,Mzz);
vmess("Mxy %.2f Mxz %.2f Myz %.2f",Mxy,Mxz,Myz);
}
if (wav->random) vmess("Wavelet has a random signature with fmax=%.2f", wav->fmax);
if (src->n>1) {
vmess("*******************************************");
vmess("*********** source array info *************");
vmess("*******************************************");
vmess("Areal source array is defined with %li sources (x=%li, y=%li).",nsrc,npxsrc,npysrc);
vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
if (src->plane) vmess("Computed px-value = %f. and py-value = %f",px,py);
}
if (src->random) {
vmess("Sources are placed at random locations in domain: ");
vmess(" x[%.2f : %.2f] y[%.2f : %.2f] z[%.2f : %.2f] ", xsrc1, xsrc2, ysrc1, ysrc2, zsrc1, zsrc2);
vmess(" and all start in time window t[%.3f : %.3f].", tsrc1, tsrc2);
vmess(" after time %.3f the sources will not be active anymore.", tactive);
}
}
/* define snapshots and beams */
if (!getparfloat("tsnap1", &tsnap1)) tsnap1=0.1;
if (!getparfloat("tsnap2", &tsnap2)) tsnap2=0.0;
if (!getparfloat("dtsnap", &dtsnap)) dtsnap=0.1;
if (!getparfloat("dxsnap", &dxsnap)) dxsnap=dx;
if (!getparfloat("dysnap", &dysnap)) dysnap=dy;
if (!getparfloat("dzsnap", &dzsnap)) dzsnap=dz;
if (!getparfloat("xsnap1", &xsnap1)) xsnap1=sub_x0;
if (!getparfloat("xsnap2", &xsnap2)) xsnap2=xmax;
if (!getparfloat("ysnap1", &ysnap1)) ysnap1=sub_y0;
if (!getparfloat("ysnap2", &ysnap2)) ysnap2=ymax;
if (!getparfloat("zsnap1", &zsnap1)) zsnap1=sub_z0;
if (!getparfloat("zsnap2", &zsnap2)) zsnap2=zmax;
if (!getparlong("sna_vxvztime", &sna->vxvztime)) sna->vxvztime=0;
if (!getparlong("beam", &sna->beam)) sna->beam=0;
if (!getparlong("snapwithbnd", &sna->withbnd)) sna->withbnd=0;
if (!getparlong("sna_type_vz", &sna->type.vz)) sna->type.vz=1;
if (!getparlong("sna_type_vy", &sna->type.vy)) sna->type.vy=0;
if (!getparlong("sna_type_vx", &sna->type.vx)) sna->type.vx=0;
if (mod->ischeme>2) {
sna->type.p=0;
if (!getparlong("sna_type_txx", &sna->type.txx)) sna->type.txx=0;
if (!getparlong("sna_type_tyy", &sna->type.tyy)) sna->type.tyy=0;
if (!getparlong("sna_type_tzz", &sna->type.tzz)) sna->type.tzz=0;
if (!getparlong("sna_type_txz", &sna->type.txz)) sna->type.txz=0;
if (!getparlong("sna_type_txy", &sna->type.txy)) sna->type.txy=0;
if (!getparlong("sna_type_tyz", &sna->type.tyz)) sna->type.tyz=0;
if (!getparlong("sna_type_pp", &sna->type.pp)) sna->type.pp=0;
if (!getparlong("sna_type_ss", &sna->type.ss)) sna->type.ss=0;
}
else {
if (!getparlong("sna_type_p", &sna->type.p)) sna->type.p=1;
sna->type.txx=0;
sna->type.tyy=0;
sna->type.tzz=0;
sna->type.txz=0;
sna->type.txy=0;
sna->type.tyz=0;
sna->type.pp=0;
sna->type.ss=0;
}
sna->nsnap = 0;
if (tsnap2 >= tsnap1) {
sna_nrsna = 1+NINT((tsnap2-tsnap1)/dtsnap);
sna->skipdt = MAX(1,NINT(dtsnap/dt));
sna->skipdx = MAX(1,NINT(dxsnap/dx)); sna->skipdy = MAX(1,NINT(dysnap/dy));
sna->skipdz = MAX(1,NINT(dzsnap/dz));
sna->delay = NINT(tsnap1/dt);
isnapmax1 = (sna_nrsna-1)*sna->skipdt;
isnapmax2 = floor( (mod->nt-(sna->delay + 1))/sna->skipdt) * sna->skipdt;
isnapmax = (sna->delay + 1) + MIN(isnapmax1,isnapmax2);
sna->nsnap = floor((isnapmax-(sna->delay + 1))/sna->skipdt) + 1;
sna->x1=NINT((MIN(MAX(sub_x0,xsnap1),xmax)-sub_x0)/dx);
sna->x2=NINT((MIN(MAX(sub_x0,xsnap2),xmax)-sub_x0)/dx);
sna->y1=NINT((MIN(MAX(sub_y0,ysnap1),ymax)-sub_y0)/dy);
sna->y2=NINT((MIN(MAX(sub_y0,ysnap2),ymax)-sub_y0)/dy);
sna->z1=NINT((MIN(MAX(sub_z0,zsnap1),zmax)-sub_z0)/dz);
sna->z2=NINT((MIN(MAX(sub_z0,zsnap2),zmax)-sub_z0)/dz);
dxsnap=dx*sna->skipdx;
dysnap=dy*sna->skipdy;
dzsnap=dz*sna->skipdz;
sna->nx=1+(((sna->x2-sna->x1))/sna->skipdx);
sna->ny=1+(((sna->y2-sna->y1))/sna->skipdy);
sna->nz=1+(((sna->z2-sna->z1))/sna->skipdz);
if (verbose) {
vmess("*******************************************");
vmess("************* snap shot info **************");
vmess("*******************************************");
vmess("tsnap1 = %f tsnap2 = %f ", tsnap1, tsnap2);
vmess("dtsnap = %f Nsnap = %li ", dtsnap, sna->nsnap);
vmess("nzsnap = %li nysnap = %li nxsnap = %li ", sna->nz, sna->ny, sna->nx);
vmess("dzsnap = %f dysnap = %f dxsnap = %f ", dzsnap, dysnap, dxsnap);
vmess("zmin = %f zmax = %f ", sub_z0+dz*sna->z1, sub_z0+dz*sna->z2);
vmess("ymin = %f ymax = %f ", sub_y0+dy*sna->y1, sub_y0+dy*sna->y2);
vmess("xmin = %f xmax = %f ", sub_x0+dx*sna->x1, sub_x0+dx*sna->x2);
if (sna->vxvztime) vmess("vx/vy/vz snapshot time : t+0.5*dt ");
else vmess("vx/vy/vz snapshot time : t-0.5*dt ");
fprintf(stderr," %s: Snapshot types : ",xargv[0]);
if (sna->type.vz) fprintf(stderr,"Vz ");
if (sna->type.vy) fprintf(stderr,"Vy ");
if (sna->type.vx) fprintf(stderr,"Vx ");
if (sna->type.p) fprintf(stderr,"p ");
if (mod->ischeme>2) {
if (sna->type.txx) fprintf(stderr,"Txx ");
if (sna->type.tyy) fprintf(stderr,"Tyy ");
if (sna->type.tzz) fprintf(stderr,"Tzz ");
if (sna->type.txz) fprintf(stderr,"Txz ");
if (sna->type.txy) fprintf(stderr,"Txy ");
if (sna->type.tyz) fprintf(stderr,"Tyz ");
if (sna->type.pp) fprintf(stderr,"P ");
if (sna->type.ss) fprintf(stderr,"S ");
}
fprintf(stderr,"\n");
}
}
else {
sna->nsnap = 0;
if (verbose) vmess("*************** no snapshots **************");
}
if (sna->beam) {
sna->skipdx = MAX(1,NINT(dxsnap/dx));
sna->skipdy = MAX(1,NINT(dysnap/dy));
sna->skipdz = MAX(1,NINT(dzsnap/dz));
sna->x1=NINT((MIN(MAX(sub_x0,xsnap1),xmax)-sub_x0)/dx);
sna->x2=NINT((MIN(MAX(sub_x0,xsnap2),xmax)-sub_x0)/dx);
sna->y1=NINT((MIN(MAX(sub_y0,ysnap1),ymax)-sub_y0)/dy);
sna->y2=NINT((MIN(MAX(sub_y0,ysnap2),ymax)-sub_y0)/dy);
sna->z1=NINT((MIN(MAX(sub_z0,zsnap1),zmax)-sub_z0)/dz);
sna->z2=NINT((MIN(MAX(sub_z0,zsnap2),zmax)-sub_z0)/dz);
dxsnap=dx*sna->skipdx;
dysnap=dy*sna->skipdy;
dzsnap=dz*sna->skipdz;
sna->nx=1+(((sna->x2-sna->x1))/sna->skipdx); sna->ny=1+(((sna->y2-sna->y1))/sna->skipdy);
sna->nz=1+(((sna->z2-sna->z1))/sna->skipdz);
if (verbose) {
vmess("*******************************************");
vmess("**************** beam info ****************");
vmess("*******************************************");
vmess("nzsnap = %li nysnap =%li nxsnap = %li ", sna->nz, sna->ny, sna->nx);
vmess("dzsnap = %f dysnap = %f dxsnap = %f ", dzsnap, dysnap, dxsnap);
vmess("zmin = %f zmax = %f ", sub_z0+dz*sna->z1, sub_z0+dz*sna->z2);
vmess("ymin = %f ymax = %f ", sub_y0+dy*sna->y1, sub_y0+dy*sna->y2);
vmess("xmin = %f xmax = %f ", sub_x0+dx*sna->x1, sub_x0+dx*sna->x2);
fprintf(stderr," %s: Beam types : ",xargv[0]);
if (sna->type.vz) fprintf(stderr,"Vz ");
if (sna->type.vy) fprintf(stderr,"Vy ");
if (sna->type.vx) fprintf(stderr,"Vx ");
if (sna->type.p) fprintf(stderr,"p ");
if (mod->ischeme>2) {
if (sna->type.txx) fprintf(stderr,"Txx ");
if (sna->type.tyy) fprintf(stderr,"Tyy ");
if (sna->type.tzz) fprintf(stderr,"Tzz ");
if (sna->type.txz) fprintf(stderr,"Txz ");
if (sna->type.txy) fprintf(stderr,"Txy ");
if (sna->type.tyz) fprintf(stderr,"Tyz ");
if (sna->type.pp) fprintf(stderr,"P ");
if (sna->type.ss) fprintf(stderr,"S ");
}
fprintf(stderr,"\n");
}
}
else {
if (verbose) vmess("**************** no beams *****************");
}
/* define receivers */
if (!getparlong("largeSUfile",&largeSUfile)) largeSUfile=0;
if (!getparlong("sinkdepth",&rec->sinkdepth)) rec->sinkdepth=0;
if (!getparlong("sinkdepth_src",&src->sinkdepth)) src->sinkdepth=0;
if (!getparlong("sinkvel",&rec->sinkvel)) rec->sinkvel=0;
if (!getparfloat("dtrcv",&dtrcv)) dtrcv=0.004;
/* TODO check if dtrcv is integer multiple of dt */
rec->skipdt=NINT(dtrcv/dt);
dtrcv = mod->dt*rec->skipdt;
if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
if (!getparlong("rec_ntsam",&rec->nt)) rec->nt=NINT((mod->tmod-rdelay)/dtrcv)+1;
if (!getparlong("rec_int_p",&rec->int_p)) rec->int_p=0;
if (!getparlong("rec_int_vx",&rec->int_vx)) rec->int_vx=0;
if (!getparlong("rec_int_vy",&rec->int_vy)) rec->int_vy=0;
if (!getparlong("rec_int_vz",&rec->int_vz)) rec->int_vz=0;
if (!getparlong("max_nrec",&rec->max_nrec)) rec->max_nrec=15000;
if (!getparlong("scale",&rec->scale)) rec->scale=0;
if (!getparfloat("dxspread",&dxspread)) dxspread=0;
if (!getparfloat("dyspread",&dyspread)) dyspread=0;
if (!getparfloat("dzspread",&dzspread)) dzspread=0;
rec->nt=MIN(rec->nt, NINT((mod->tmod-rdelay)/dtrcv)+1);
/* allocation of receiver arrays is done in recvPar */
/* calculates the receiver coordinates */
recvPar3D(rec, sub_x0, sub_y0, sub_z0, dx, dy, dz, nx, ny, nz);
if (!getparlong("rec_type_vz", &rec->type.vz)) rec->type.vz=1;
if (!getparlong("rec_type_vy", &rec->type.vy)) rec->type.vy=0;
if (!getparlong("rec_type_vx", &rec->type.vx)) rec->type.vx=0;
if (!getparlong("rec_type_ud", &rec->type.ud)) rec->type.ud=0; if (mod->ischeme!=1 && rec->type.ud==1) {
warn("Receiver decomposition only implemented for acoustis scheme (1)");
}
if (mod->ischeme>2) {
rec->type.p=0;
if (!getparlong("rec_type_txx", &rec->type.txx)) rec->type.txx=0;
if (!getparlong("rec_type_tyy", &rec->type.tyy)) rec->type.tyy=0;
if (!getparlong("rec_type_tzz", &rec->type.tzz)) rec->type.tzz=0;
if (!getparlong("rec_type_txz", &rec->type.txz)) rec->type.txz=0;
if (!getparlong("rec_type_txy", &rec->type.txy)) rec->type.txy=0;
if (!getparlong("rec_type_tyz", &rec->type.tyz)) rec->type.tyz=0;
if (!getparlong("rec_type_pp", &rec->type.pp)) rec->type.pp=0;
if (!getparlong("rec_type_ss", &rec->type.ss)) rec->type.ss=0;
/* for up and downgoing waves store all x-positons for Vz, Vx, Txz, Tzz into an array */
}
else {
if (!getparlong("rec_type_p", &rec->type.p)) rec->type.p=1;
rec->type.txx=0;
rec->type.tyy=0;
rec->type.tzz=0;
rec->type.txz=0;
rec->type.txy=0;
rec->type.tyz=0;
rec->type.pp=0;
rec->type.ss=0;
/* for up and downgoing waves store all x-positons for P and Vz into an array */
}
/* receivers are on a circle, use default interpolation to real (not on a grid-point) receiver position */
if (getparfloat("rrcv", &rrcv)) {
if (!getparlong("rec_int_p",&rec->int_p)) rec->int_p=3;
if (!getparlong("rec_int_vx",&rec->int_vx)) rec->int_vx=3;
if (!getparlong("rec_int_vy",&rec->int_vy)) rec->int_vy=3;
if (!getparlong("rec_int_vz",&rec->int_vz)) rec->int_vz=3;
}
if (rec->int_p==3) {
rec->int_vx=3;
rec->int_vy=3;
rec->int_vz=3;
}
if (verbose) {
if (rec->n) {
dxrcv = rec->xr[MIN(1,rec->n-1)]-rec->xr[0];
dyrcv = rec->yr[MIN(1,rec->n-1)]-rec->yr[0];
dzrcv = rec->zr[MIN(1,rec->n-1)]-rec->zr[0];
vmess("*******************************************");
vmess("************* receiver info ***************");
vmess("*******************************************");
vmess("ntrcv = %li nrcv = %li ", rec->nt, rec->n);
vmess("dtrcv = %f ", dtrcv );
vmess("dzrcv = %f dyrcv = %f dxrcv = %f ", dzrcv, dyrcv, dxrcv);
vmess("time-delay = %f = points = %li", rdelay, rec->delay);
if ( fmax > (1.0/(2.0*dtrcv)) ) {
vwarn("Receiver time sampling (dtrcv) is aliased.");
vwarn("time sampling should be < %.6f", 1.0/(2.0*fmax) );
}
vmess("Receiver sampling can be => %.6e", 1.0/(2.0*fmax));
vmess("Receiver array at coordinates: ");
vmess("zmin = %f zmax = %f ", rec->zr[0]+sub_z0, rec->zr[rec->n-1]+sub_z0);
vmess("ymin = %f ymax = %f ", rec->yr[0]+sub_y0, rec->yr[rec->n-1]+sub_y0);
vmess("xmin = %f xmax = %f ", rec->xr[0]+sub_x0, rec->xr[rec->n-1]+sub_x0);
vmess("which are gridpoints: ");
vmess("izmin = %li izmax = %li ", rec->z[0], rec->z[rec->n-1]);
vmess("iymin = %li iymax = %li ", rec->y[0], rec->y[rec->n-1]);
vmess("ixmin = %li ixmax = %li ", rec->x[0], rec->x[rec->n-1]);
if (rec->type.p) {
fprintf(stderr," %s: Receiver interpolation for P: ",xargv[0]);
if(rec->int_p==0) fprintf(stderr,"p->p\n");
if(rec->int_p==1) fprintf(stderr,"p->vz\n"); if(rec->int_p==2) fprintf(stderr,"p->vx\n");
if(rec->int_p==3) fprintf(stderr,"interpolate to actual (no-grid) position of receiver\n");
}
if (rec->type.vx) {
fprintf(stderr," %s: Receiver interpolation for Vx: ",xargv[0]);
if(rec->int_vx==0) fprintf(stderr,"vx->vx\n");
if(rec->int_vx==1) fprintf(stderr,"vx->vz\n");
if(rec->int_vx==2) fprintf(stderr,"vx->txx/tzz\n");
if(rec->int_vx==3) fprintf(stderr,"interpolate to real(no-grid) position of receiver\n");
}
if (rec->type.vy) {
fprintf(stderr," %s: Receiver interpolation for Vx: ",xargv[0]);
if(rec->int_vy==0) fprintf(stderr,"vy->vy\n");
if(rec->int_vy==1) fprintf(stderr,"vy->vz\n");
if(rec->int_vy==2) fprintf(stderr,"vy->tyy/tzz\n");
if(rec->int_vy==3) fprintf(stderr,"interpolate to real(no-grid) position of receiver\n");
}
if (rec->type.vz) {
fprintf(stderr," %s: Receiver interpolation for Vz: ",xargv[0]);
if(rec->int_vz==0) fprintf(stderr,"vz->vz\n");
if(rec->int_vz==1) fprintf(stderr,"vz->vx\n");
if(rec->int_vz==2) fprintf(stderr,"vz->txx/tzz(P)\n");
if(rec->int_vz==3) fprintf(stderr,"interpolate to real(no-grid) position of receiver\n");
}
fprintf(stderr," %s: Receiver types : ",xargv[0]);
if (rec->type.vz) fprintf(stderr,"Vz ");
if (rec->type.vy) fprintf(stderr,"Vy ");
if (rec->type.vx) fprintf(stderr,"Vx ");
if (rec->type.p) fprintf(stderr,"p ");
if (rec->type.ud) fprintf(stderr,"P+ P- ");
if (mod->ischeme>2) {
if (rec->type.txx) fprintf(stderr,"Txx ");
if (rec->type.tyy) fprintf(stderr,"Tyy ");
if (rec->type.tzz) fprintf(stderr,"Tzz ");
if (rec->type.txz) fprintf(stderr,"Txz ");
if (rec->type.txy) fprintf(stderr,"Txy ");
if (rec->type.tyz) fprintf(stderr,"Tyz ");
if (rec->type.pp) fprintf(stderr,"P ");
if (rec->type.ss) fprintf(stderr,"S ");
}
fprintf(stderr,"\n");
if ( ( ((mod->nt*mod->dt-rec->delay)/rec->skipdt)+1) > 16384) {
vwarn("Number of samples in receiver file is larger that SU can handle ");
vwarn("use the paramater rec_ntsam=nt (with nt < 16384) to avoid this");
}
if ((mod->nt-rec->delay)*mod->dt > rec->nt*dtrcv) {
long nfiles = ceil((mod->nt*mod->dt)/(rec->nt*dtrcv));
long lastn = floor((mod->nt)%(rec->nt*rec->skipdt)/rec->skipdt)+1;
vmess("Receiver recordings will be written to %li files",nfiles);
vmess("Last file will contain %li samples",lastn);
}
}
else {
vmess("*************** no receivers **************");
}
}
return 0;
}