diff --git a/fdelmodc3D/defineSource3D.c b/fdelmodc3D/defineSource3D.c
index 01cc29b40617b99b454f8f3d9a53bfb85244f4fa..3f5d89ab290db022a87a2a151c170270b0898f77 100644
--- a/fdelmodc3D/defineSource3D.c
+++ b/fdelmodc3D/defineSource3D.c
@@ -28,17 +28,17 @@ typedef struct _dcomplexStruct { /* complex number */
} dcomplex;
#endif/* complex */
-int optncr(int n);
+long loptncr(long n);
void rc1fft(float *rdata, complex *cdata, int n, int sign);
void cr1fft(complex *cdata, float *rdata, int n, int sign);
-int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2, float d1, float d2);
-int writesufilesrcnwav(char *filename, float **src_nwav, wavPar wav, int n1, int n2, float f1, float f2, float d1, float d2);
+long writesufile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2);
+long writesufilesrcnwav3D(char *filename, float **src_nwav, wavPar wav, long n1, long n2, float f1, float f2, float d1, float d2);
float gaussGen();
float normal(double x,double mu,double sigma);
-int comp (const float *a, const float *b);
+long comp (const float *a, const float *b);
void spline3(float x1, float x2, float z1, float z2, float dzdx1, float dzdx2, float *a, float *b, float *c, float *d);
-int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend, int verbose);
+long randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend, long verbose);
/* random number generators */
double dcmwc4096();
@@ -50,14 +50,14 @@ void seedCMWC4096(void);
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
-int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwav, int reverse, int verbose)
+long defineSource3D(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwav, long reverse, long verbose)
{
FILE *fp;
size_t nread;
- int optn, nfreq, i, j, k, iwmax, tracesToDo;
- int iw, n1, namp, optnscale, nfreqscale;
+ long optn, nfreq, i, j, k, iwmax, tracesToDo;
+ long iw, n1, namp, optnscale, nfreqscale;
float scl, d1, df, deltom, om, tshift;
float amp1, amp2, amp3;
float *trace, maxampl, scale;
@@ -72,7 +72,6 @@ int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwa
for (i=0; i<8192; i++) {
amp1 = dcmwc4096();
}
-
}
else {
@@ -100,12 +99,12 @@ int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwa
fclose(fp);
}
- optn = optncr(n1);
+ optn = loptncr(n1);
nfreq = optn/2 + 1;
if (wav.nt != wav.ns) {
vmess("Sampling in wavelet is %e while for modeling is set to %e", wav.ds, mod.dt);
vmess("Wavelet sampling will be FFT-interpolated to sampling of modeling");
- vmess("file_src Nt=%d sampling after interpolation=%d", wav.ns, wav.nt);
+ vmess("file_src Nt=%li sampling after interpolation=%li", wav.ns, wav.nt);
optnscale = wav.nt;
nfreqscale = optnscale/2 + 1;
}
@@ -113,7 +112,7 @@ int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwa
optnscale = optn;
nfreqscale = optnscale/2 + 1;
}
-// fprintf(stderr,"define S optn=%d ns=%d %e nt=%d %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
+// fprintf(stderr,"define S optn=%li ns=%li %e nt=%li %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
ctrace = (complex *)calloc(nfreqscale,sizeof(complex));
trace = (float *)calloc(optnscale,sizeof(float));
@@ -192,7 +191,7 @@ int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwa
maxampl = MAX(maxampl,fabs(src_nwav[i][j]));
}
}
- if (verbose > 3) vmess("Wavelet sampling (FFT-interpolated) done for trace %d", i);
+ if (verbose > 3) vmess("Wavelet sampling (FFT-interpolated) done for trace %li", i);
}
}
/* set values smaller than 1e-5 maxampl to zero */
@@ -212,12 +211,12 @@ int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwa
for (i=0; i<wav.nx; i++) {
if (src.distribution) {
scl = gaussGen()*src.amplitude;
- k = (int)MAX(MIN(namp*(scl+5*src.amplitude)/(10*src.amplitude),namp-1),0);
+ k = (long)MAX(MIN(namp*(scl+5*src.amplitude)/(10*src.amplitude),namp-1),0);
d1 = 10.0*src.amplitude/(namp-1);
}
else {
scl = (float)(drand48()-0.5)*src.amplitude;
- k = (int)MAX(MIN(namp*(scl+1*src.amplitude)/(2*src.amplitude),namp-1),0);
+ k = (long)MAX(MIN(namp*(scl+1*src.amplitude)/(2*src.amplitude),namp-1),0);
d1 = 2.0*src.amplitude/(namp-1);
}
@@ -229,41 +228,21 @@ int defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwa
src_nwav[i][j] *= scl;
}
}
- if (verbose>2) writesufile("src_ampl.su", trace, namp, 1, -5*src.amplitude, 0.0, d1, 1);
-/*
- qsort(trace,wav.nx,sizeof(float), comp);
- for (i=0; i<wav.nx; i++) {
- scl = trace[i];
- trace[i] = normal(scl, 0.0, src.amplitude);
- }
- if (verbose>2) writesufile("src_ampl.su", trace, wav.nx, 1, -5*src.amplitude, 0.0, d1, 1);
-*/
+ if (verbose>2) writesufile3D("src_ampl.su", trace, namp, 1, -5*src.amplitude, 0.0, d1, 1);
free(trace);
}
- if (verbose>3) writesufilesrcnwav("src_nwav.su", src_nwav, wav, wav.nt, wav.nx, 0.0, 0.0, wav.dt, 1);
-
-/* set maximum amplitude in source file to 1.0 */
-/*
- assert(maxampl != 0.0);
- scl = wav.dt/(maxampl);
- scl = 1.0/(maxampl);
- for (i=0; i<wav.nx; i++) {
- for (j=0; j<n1; j++) {
- src_nwav[i*n1+j] *= scl;
- }
- }
-*/
+ if (verbose>3) writesufilesrcnwav3D("src_nwav.su", src_nwav, wav, wav.nt, wav.nx, 0.0, 0.0, wav.dt, 1);
return 0;
}
-int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend, int verbose)
+long randomWavelet3D(wavPar wav, srcPar src, float *trace, float tbeg, float tend, long verbose)
{
- int optn, nfreq, j, iwmax;
- int iw, n1, itbeg, itmax, nsmth;
+ long optn, nfreq, j, iwmax;
+ long iw, n1, itbeg, itmax, nsmth;
float df, amp1;
float *rtrace;
float x1, x2, z1, z2, dzdx1, dzdx2, a, b, c, d, t;
@@ -271,7 +250,7 @@ int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend,
n1 = wav.nt; /* this is set to the maximum length (tlength/dt) */
- optn = optncr(2*n1);
+ optn = loptncr(2*n1);
nfreq = optn/2 + 1;
ctrace = (complex *)calloc(nfreq,sizeof(complex));
rtrace = (float *)calloc(optn,sizeof(float));
@@ -322,7 +301,6 @@ int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend,
dzdx1 = 0.0;
x2 = nsmth;
z2 = rtrace[itbeg+nsmth];
-// dzdx2 = (rtrace[itbeg+(nsmth+1)]-rtrace[itbeg+(nsmth-1)])/(2.0);
dzdx2 = (rtrace[itbeg+nsmth-2]-8.0*rtrace[itbeg+nsmth-1]+
8.0*rtrace[itbeg+nsmth+1]-rtrace[itbeg+nsmth+2])/(12.0);
spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
@@ -333,19 +311,16 @@ int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend,
x1 = 0.0;
z1 = rtrace[itmax-nsmth];
-// dzdx1 = (rtrace[itmax-(nsmth-1)]-rtrace[itmax-(nsmth+1)])/(2.0);
dzdx1 = (rtrace[itmax-nsmth-2]-8.0*rtrace[itmax-nsmth-1]+
8.0*rtrace[itmax-nsmth+1]-rtrace[itmax-nsmth+2])/(12.0);
x2 = nsmth;
z2 = 0.0;
dzdx2 = 0.0;
-// fprintf(stderr,"x1=%f z1=%f d=%f x2=%f, z2=%f d=%f\n",x1,z1,dzdx1,x2,z2,dzdx2);
spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
for (j=0; j<nsmth; j++) {
t = j;
rtrace[itmax-nsmth+j] = a*t*t*t+b*t*t+c*t+d;
-// fprintf(stderr,"a=%f b=%f c=%f d=%f rtrace%d=%f \n",a,b,c,d,j,rtrace[itmax-nsmth+j]);
}
for (j=itbeg; j<itmax; j++) trace[j-itbeg] = rtrace[j];
@@ -361,7 +336,7 @@ float normal(double x,double mu,double sigma)
return (float)(1.0/(2.0*M_PI*sigma*sigma))*exp(-1.0*(((x-mu)*(x-mu))/(2.0*sigma*sigma)) );
}
-int comp (const float *a, const float *b)
+long comp (const float *a, const float *b)
{
if (*a==*b)
return 0;
diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c
index 559849883a79b92d43b0294ca8a122e831e5db11..a79ede1973ad96f941f5eeb096f8486954c7a5d2 100644
--- a/fdelmodc3D/fdelmodc3D.c
+++ b/fdelmodc3D/fdelmodc3D.c
@@ -21,7 +21,7 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, float *l2m, float *lam, float *muu, float *tss, float *tes, float *tep);
-long defineSource(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwav, long reverse, long verbose);
+long defineSource3D(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwav, long reverse, long verbose);
long writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
@@ -111,11 +111,6 @@ char *sdoc[] = {
" bottom=4 .......... type of boundary on bottom edge of model",
" front=4 ........... type of boundary on front edge of model",
" back=4 ............ type of boundary on back edge of model",
-//" tapleft=0 ......... =1: taper left edge of model",
-//" tapright=0 ........ =1: taper right edge of model",
-//" taptop=0 .......... =1: taper top edge of model",
-//" tapbottom=0 ....... =1: taper bottom edge of model",
-//" cfree=0 ........... 1=free surface",
" grid_dir=0 ........ direction of time modeling (1=reverse time)",
" Qp=15 ............. global Q-value for P-waves in visco-elastic (ischeme=2,4)",
" file_qp= .......... model file Qp values as function of depth",
@@ -135,8 +130,6 @@ char *sdoc[] = {
" - 1=monopole",
" - 2=dipole +/- vertical oriented",
" - 3=dipole - + horizontal oriented",
-//" - 4=dipole +/0/-",
-//" - 5=dipole + -",
" dip=0 ............. dip for double-couple source",
" strike=0 .......... strike for double-couple source",
" xsrc=middle ....... x-position of (first) shot ",
@@ -221,7 +214,6 @@ char *sdoc[] = {
" zrcv2=zrcv1 ....... last z-position of linear receiver array(s)",
" dzrcv=0.0 ......... z-position increment of receivers in linear array(s)",
" dtrcv=.004 ........ desired sampling in receiver data (seconds)",
-//" max_nrec=15000 .... maximum number of receivers", not needed anymore
" xrcva= ............ defines receiver array x-positions",
" yrcva= ............ defines receiver array y-positions",
" zrcva= ............ defines receiver array z-positions",
@@ -231,11 +223,8 @@ char *sdoc[] = {
" ozrcv=0.0 ......... z-center position of circle",
" dphi=2 ............ angle between receivers on circle ",
" rcv_txt=........... text file with receiver coordinates. Col 1: x, Col. 2: z",
-//" largeSUfile=0 ..... writing large SU file (nt > 64000)",
" rec_ntsam=nt ...... maximum number of time samples in file_rcv files",
" rec_delay=0 ....... time in seconds to start recording: recorded time = tmod - rec_delay",
-//" dxspread=0 ........ if nshot > 1: x-shift of rcv spread",
-//" dzspread=0 ........ if nshot > 1: z-shift of rcv spread",
" rec_type_p=1 ...... p registration _rp",
" rec_type_vz=1 ..... Vz registration _rvz",
" rec_type_vy=0 ..... Vy registration _rvy",
@@ -303,7 +292,7 @@ int main(int argc, char **argv)
float sinkvel, npeshot;
double t0, t1, t2, t3, tt, tinit;
size_t size, sizem, nsamp;
- long n1, ix, iy, iz, ir, ishot, i;
+ long n1, n2, ix, iy, iz, ir, ishot, i;
long ioPx, ioPy, ioPz;
long it0, it1, its, it, fileno, isam;
long ixsrc, iysrc, izsrc, is0, is1;
@@ -331,6 +320,7 @@ int main(int argc, char **argv)
/* allocate arrays for model parameters: the different schemes use different arrays */
n1 = mod.naz;
+ n2 = mod.nax;
sizem=mod.nax*mod.naz*mod.nay;
rox = (float *)calloc(sizem,sizeof(float));
@@ -386,44 +376,56 @@ int main(int argc, char **argv)
}
}
- defineSource(wav, src, mod, rec, src_nwav, mod.grid_dir, verbose);
+ defineSource3D(wav, src, mod, rec, src_nwav, mod.grid_dir, verbose);
/* allocate arrays for wavefield and receiver arrays */
vx = (float *)calloc(sizem,sizeof(float));
+ vy = (float *)calloc(sizem,sizeof(float));
vz = (float *)calloc(sizem,sizeof(float));
tzz = (float *)calloc(sizem,sizeof(float)); /* =P field for acoustic */
if (mod.ischeme>2) {
txz = (float *)calloc(sizem,sizeof(float));
+ txy = (float *)calloc(sizem,sizeof(float));
+ tyz = (float *)calloc(sizem,sizeof(float));
txx = (float *)calloc(sizem,sizeof(float));
+ tyy = (float *)calloc(sizem,sizeof(float));
}
size = rec.n*rec.nt;
if (rec.type.vz) rec_vz = (float *)calloc(size,sizeof(float));
+ if (rec.type.vy) rec_vy = (float *)calloc(size,sizeof(float));
if (rec.type.vx) rec_vx = (float *)calloc(size,sizeof(float));
if (rec.type.p) rec_p = (float *)calloc(size,sizeof(float));
if (rec.type.txx) rec_txx = (float *)calloc(size,sizeof(float));
+ if (rec.type.tyy) rec_tyy = (float *)calloc(size,sizeof(float));
if (rec.type.tzz) rec_tzz = (float *)calloc(size,sizeof(float));
if (rec.type.txz) rec_txz = (float *)calloc(size,sizeof(float));
+ if (rec.type.txy) rec_txy = (float *)calloc(size,sizeof(float));
+ if (rec.type.tyz) rec_tyz = (float *)calloc(size,sizeof(float));
if (rec.type.pp) rec_pp = (float *)calloc(size,sizeof(float));
if (rec.type.ss) rec_ss = (float *)calloc(size,sizeof(float));
if (rec.type.ud) {
rec_udvz = (float *)calloc(mod.nax*rec.nt,sizeof(float));
rec_udp = (float *)calloc(mod.nax*rec.nt,sizeof(float));
}
- /* get velcity and density at first receiver location */
- ir = mod.ioZz + rec.z[0]+(rec.x[0]+mod.ioZx)*n1;
+ /* get velocity and density at first receiver location */
+ ir = mod.ioZz + rec.z[0]+(rec.x[0]+mod.ioZx)*n1+(rec.y[0]+mod.ioZy)*n1*n2;
rec.rho = mod.dt/(mod.dx*roz[ir]);
rec.cp = sqrt(l2m[ir]*(roz[ir]))*mod.dx/mod.dt;
if(sna.beam) {
size = sna.nz*sna.nx;
if (sna.type.vz) beam_vz = (float *)calloc(size,sizeof(float));
+ if (sna.type.vy) beam_vy = (float *)calloc(size,sizeof(float));
if (sna.type.vx) beam_vx = (float *)calloc(size,sizeof(float));
if (sna.type.p) beam_p = (float *)calloc(size,sizeof(float));
if (sna.type.txx) beam_txx = (float *)calloc(size,sizeof(float));
+ if (sna.type.tyy) beam_tyy = (float *)calloc(size,sizeof(float));
if (sna.type.tzz) beam_tzz = (float *)calloc(size,sizeof(float));
if (sna.type.txz) beam_txz = (float *)calloc(size,sizeof(float));
+ if (sna.type.txy) beam_txy = (float *)calloc(size,sizeof(float));
+ if (sna.type.tyz) beam_tyz = (float *)calloc(size,sizeof(float));
if (sna.type.pp) beam_pp = (float *)calloc(size,sizeof(float));
if (sna.type.ss) beam_ss = (float *)calloc(size,sizeof(float));
}
@@ -446,52 +448,67 @@ int main(int argc, char **argv)
of the first receiver to sink through to the next layer. */
ioPx=mod.ioPx;
+ ioPy=mod.ioPy;
ioPz=mod.ioPz;
if (bnd.lef==4 || bnd.lef==2) ioPx += bnd.ntap;
+ if (bnd.fro==4 || bnd.fro==2) ioPy += bnd.ntap;
if (bnd.top==4 || bnd.top==2) ioPz += bnd.ntap;
- if (rec.sinkvel) sinkvel=l2m[(rec.x[0]+ioPx)*n1+rec.z[0]+ioPz];
+ if (rec.sinkvel) sinkvel=l2m[(rec.y[0]+ioPy)*n1*n2+(rec.x[0]+ioPx)*n1+rec.z[0]+ioPz];
else sinkvel = 0.0;
/* sink receivers to value different than sinkvel */
for (ir=0; ir<rec.n; ir++) {
iz = rec.z[ir];
+ iy = rec.y[ir];
ix = rec.x[ir];
- while(l2m[(ix+ioPx)*n1+iz+ioPz] == sinkvel) iz++;
+ while(l2m[(iy+ioPy)*n1*n2+(ix+ioPx)*n1+iz+ioPz] == sinkvel) iz++;
rec.z[ir]=iz+rec.sinkdepth;
rec.zr[ir]=rec.zr[ir]+(rec.z[ir]-iz)*mod.dz;
-// rec.zr[ir]=rec.z[ir]*mod.dz;
- if (verbose>3) vmess("receiver position %d at grid[ix=%d, iz=%d] = (x=%f z=%f)", ir, ix+ioPx, rec.z[ir]+ioPz, rec.xr[ir]+mod.x0, rec.zr[ir]+mod.z0);
+ if (verbose>3) vmess("receiver position %li at grid[ix=%li, iy=%li iz=%li] = (x=%f y=%f z=%f)", ir, ix+ioPx, iy+ioPy, rec.z[ir]+ioPz, rec.xr[ir]+mod.x0, rec.yr[ir]+mod.y0, rec.zr[ir]+mod.z0);
}
/* sink sources to value different than zero */
for (ishot=0; ishot<shot.n; ishot++) {
iz = shot.z[ishot];
+ iy = shot.y[ishot];
ix = shot.x[ishot];
- while(l2m[(ix+ioPx)*n1+iz+ioPz] == 0.0) iz++;
+ while(l2m[(iy+ioPy)*n1*n2+(ix+ioPx)*n1+iz+ioPz] == 0.0) iz++;
shot.z[ishot]=iz+src.sinkdepth;
}
/* scan for free surface boundary in case it has a topography */
- for (ix=0; ix<mod.nx; ix++) {
- iz = ioPz;
- while(l2m[(ix+ioPx)*n1+iz] == 0.0) iz++;
- bnd.surface[ix+ioPx] = iz;
- if ((verbose>3) && (iz != ioPz)) vmess("Topgraphy surface x=%.2f z=%.2f", mod.x0+mod.dx*ix, mod.z0+mod.dz*(iz-ioPz));
+ for (iy=0; iy<mod.ny; iy++) {
+ for (ix=0; ix<mod.nx; ix++) {
+ iz = ioPz;
+ while(l2m[(iy+ioPy)*n1*n2+(ix+ioPx)*n1+iz] == 0.0) iz++;
+ bnd.surface[(iy+ioPy)*n2+ix+ioPx] = iz;
+ if ((verbose>3) && (iz != ioPz)) vmess("Topgraphy surface x=%.2f y=%.2f z=%.2f", mod.x0+mod.dx*ix, mod.y0+mod.dy*iy, mod.z0+mod.dz*(iz-ioPz));
+ }
}
- for (ix=0; ix<ioPx; ix++) {
- bnd.surface[ix] = bnd.surface[ioPx];
+ for (iy=0; iy<ioPy; iy++) {
+ for (ix=0; ix<ioPx; ix++) {
+ bnd.surface[iy*n2+ix] = bnd.surface[ioPy*n2+ioPx];
+ }
+ for (ix=ioPx+mod.nx; ix<mod.iePx; ix++) {
+ bnd.surface[iy*n2+ix] = bnd.surface[ioPy*n2+mod.iePx-1];
+ }
}
- for (ix=ioPx+mod.nx; ix<mod.iePx; ix++) {
- bnd.surface[ix] = bnd.surface[mod.iePx-1];
+ for (iy=ioPy+mod.ny; iy<mod.iePy; iy++) {
+ for (ix=0; ix<ioPx; ix++) {
+ bnd.surface[iy*n2+ix] = bnd.surface[(mod.iePy-1)*n2+ioPx];
+ }
+ for (ix=ioPx+mod.nx; ix<mod.iePx; ix++) {
+ bnd.surface[iy*n2+ix] = bnd.surface[(mod.iePy-1)*n2+mod.iePx-1];
+ }
}
- if (verbose>3) writeSrcRecPos(&mod, &rec, &src, &shot);
+ if (verbose>3) writeSrcRecPos3D(&mod, &rec, &src, &shot);
/* Outer loop over number of shots */
#ifdef MPI
npeshot = MAX((((float)shot.n)/((float)npes)), 1.0);
is0=ceil(pe*npeshot);
is1=MIN(ceil((pe+1)*npeshot), shot.n);
- if (verbose>1) vmess("MPI: pe=%d does shots is0 %d - is1 %d\n", pe, is0, is1);
+ if (verbose>1) vmess("MPI: pe=%li does shots is0 %li - is1 %li\n", pe, is0, is1);
#else
is0=0;
is1=shot.n;
@@ -500,10 +517,12 @@ int main(int argc, char **argv)
for (ishot=is0; ishot<is1; ishot++) {
izsrc = shot.z[ishot];
+ iysrc = shot.y[ishot];
ixsrc = shot.x[ishot];
fileno= 0;
memset(vx,0,sizem*sizeof(float));
+ memset(vy,0,sizem*sizeof(float));
memset(vz,0,sizem*sizeof(float));
memset(tzz,0,sizem*sizeof(float));
if (mod.ischeme==2) {
@@ -511,7 +530,10 @@ int main(int argc, char **argv)
}
if (mod.ischeme>2) {
memset(txz,0,sizem*sizeof(float));
+ memset(txy,0,sizem*sizeof(float));
+ memset(tyz,0,sizem*sizeof(float));
memset(txx,0,sizem*sizeof(float));
+ memset(tyy,0,sizem*sizeof(float));
}
if (mod.ischeme==4) {
memset(r,0,sizem*sizeof(float));
@@ -520,12 +542,14 @@ int main(int argc, char **argv)
}
if (verbose) {
if (!src.random) {
- vmess("Modeling source %d at gridpoints ix=%d iz=%d", ishot, shot.x[ishot], shot.z[ishot]);
- vmess(" which are actual positions x=%.2f z=%.2f", mod.x0+mod.dx*shot.x[ishot], mod.z0+mod.dz*shot.z[ishot]);
+ vmess("Modeling source %li at gridpoints ix=%li iy=%li iz=%li", ishot, shot.x[ishot], shot.y[ishot], shot.z[ishot]);
+ vmess(" which are actual positions x=%.2f y=%.2f z=%.2f", mod.x0+mod.dx*shot.x[ishot], mod.y0+mod.dy*shot.y[ishot], mod.z0+mod.dz*shot.z[ishot]);
}
- vmess("Receivers at gridpoint x-range ix=%d - %d", rec.x[0], rec.x[rec.n-1]);
+ vmess("Receivers at gridpoint x-range ix=%li - %li", rec.x[0], rec.x[rec.n-1]);
vmess(" which are actual positions x=%.2f - %.2f", mod.x0+rec.xr[0], mod.x0+rec.xr[rec.n-1]);
- vmess("Receivers at gridpoint z-range iz=%d - %d", rec.z[0], rec.z[rec.n-1]);
+ vmess("Receivers at gridpoint y-range iy=%li - %li", rec.y[0], rec.y[rec.n-1]);
+ vmess(" which are actual positions y=%.2f - %.2f", mod.y0+rec.yr[0], mod.y0+rec.yr[rec.n-1]);
+ vmess("Receivers at gridpoint z-range iz=%li - %li", rec.z[0], rec.z[rec.n-1]);
vmess(" which are actual positions z=%.2f - %.2f", mod.z0+rec.zr[0], mod.z0+rec.zr[rec.n-1]);
}
@@ -548,13 +572,13 @@ int main(int argc, char **argv)
for (it=it0; it<it1; it++) {
#pragma omp parallel default (shared) \
-shared (rox, roz, l2m, lam, mul, txx, txz, tzz, vx, vz) \
+shared (rox, roy, roz, l2m, lam, mul, txx, tyy, tzz, txz, tyz, txy, vx, vy, vz) \
shared (tss, tep, tes, r, q, p) \
shared (tinit, it0, it1, its) \
-shared(beam_vx, beam_vz, beam_txx, beam_tzz, beam_txz, beam_p, beam_pp, beam_ss) \
-shared(rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz, rec_p, rec_pp, rec_ss) \
+shared(beam_vx, beam_vy, beam_vz, beam_txx, beam_tyy, beam_tzz, beam_txz, beam_tyz, beam_txy, beam_p, beam_pp, beam_ss) \
+shared(rec_vx, rec_vy, rec_vz, rec_txx, rec_tyy, rec_tzz, rec_txz, rec_tyz, rec_txy, rec_p, rec_pp, rec_ss) \
shared (tt, t2, t3) \
-shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
+shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbose)
{
if (it==it0) {
threadAffinity();
@@ -618,7 +642,7 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
#pragma omp master
{
if ( (((it-rec.delay) % rec.skipdt)==0) && (it >= rec.delay) ) {
- int writeToFile, itwritten;
+ long writeToFile, itwritten;
writeToFile = ! ( (((it-rec.delay)/rec.skipdt)+1)%rec.nt );
itwritten = fileno*(rec.nt)*rec.skipdt;
@@ -673,14 +697,14 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
if (rec.scale==1) { /* scale receiver with distance src-rcv */
float xsrc, zsrc, Rrec, rdx, rdz;
- int irec;
+ long irec;
xsrc=mod.x0+mod.dx*ixsrc;
zsrc=mod.z0+mod.dz*izsrc;
for (irec=0; irec<rec.n; irec++) {
rdx=mod.x0+rec.xr[irec]-xsrc;
rdz=mod.z0+rec.zr[irec]-zsrc;
Rrec = sqrt(rdx*rdx+rdz*rdz);
- fprintf(stderr,"Rec %d is scaled with distance %f R=%.2f,%.2f S=%.2f,%.2f\n", irec, Rrec,rdx,rdz,xsrc,zsrc);
+ fprintf(stderr,"Rec %li is scaled with distance %f R=%.2f,%.2f S=%.2f,%.2f\n", irec, Rrec,rdx,rdz,xsrc,zsrc);
for (it=0; it<rec.nt; it++) {
rec_p[irec*rec.nt+it] *= sqrt(Rrec);
}
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 415254f5306f8a774072bb3d0eaf79cc84e867a5..691bfc1ff3a2456edec0b0b0480d989e8853b4ea 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -17,7 +17,10 @@ SRCT = marchenko3D.c \
synthesis3D.c \
applyMute3D.c \
writeData3D.c \
+ makeWindow3D.c \
ampest3D.c \
+ imaging3D.c \
+ readSnapData3D.c \
writeDataIter.c \
wallclock_time.c \
name_ext.c \
diff --git a/marchenko3D/ampest3D.c b/marchenko3D/ampest3D.c
index 820c720ae374f2f49f33eeeecb16445917487cf3..d6c6984e8af602880f8234bd2ca60282c0d7a344 100644
--- a/marchenko3D/ampest3D.c
+++ b/marchenko3D/ampest3D.c
@@ -22,94 +22,43 @@ void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long
void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift);
void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
+void deconv(float *data1, float *data2, float *decon, long nrec, long nsam,
+ float dt, float eps, float reps, long shift);
void AmpEst3D(float *f1d, float *Gd, float *ampest, long Nfoc, long nxs, long nys, long ntfft, long *ixpos, long npos,
char *file_wav, float dx, float dy, float dt)
{
long l, i, ix, iw, nfreq;
- float scl, sclt, *wavelet, *scaled, *conv, *f1dsamp;
- float dtm, dxm, cpm, rom, *trace;
+ float scl, sclt, *f1dsamp;
+ float dtm, dxm, cpm, rom, *trace, eps, reps;
FILE *fp_wav;
segy *hdrs_wav;
+ if(!getparfloat("eps", &eps)) eps=0.01;
+ if(!getparfloat("reps", &reps)) reps=0.0;
+
scl = dx*dy;
sclt = 1.0*dt/((float)ntfft);
- conv = (float *)calloc(nys*nxs*ntfft,sizeof(float));
- wavelet = (float *)calloc(ntfft,sizeof(float));
- scaled = (float *)calloc(ntfft,sizeof(float));
f1dsamp = (float *)calloc(nys*nxs*ntfft,sizeof(float));
- for (i=0; i<npos; i++) {
- ix = ixpos[i];
- for (iw=0; iw<ntfft; iw++) {
- f1dsamp[i*ntfft+iw] = f1d[l*nxs*nys*ntfft+ix*ntfft+iw];
- }
- }
-
- if (file_wav==NULL){
- corr(f1dsamp, f1dsamp, conv, nxs*nys, ntfft, dt, 0);
- for (i=0; i<nxs*nys; i++) {
- for (iw=0; iw<ntfft; iw++) {
- wavelet[iw] += scl*conv[i*ntfft+iw];
- }
- }
- }
- else {
- trace = (float *)calloc(ntfft,sizeof(float));
- hdrs_wav = (segy *)calloc(1, sizeof(segy));
- fp_wav = fopen(file_wav, "r");
- readData3D(fp_wav, trace, hdrs_wav, 0);
- fclose(fp_wav);
- corr(trace, trace, wavelet, 1, ntfft, dt, 0);
- free(hdrs_wav); free(trace);
-
- if (!getparfloat("dtm", &dtm)) dtm = 0.004;
- if (!getparfloat("dxm", &dxm)) dxm = 1.0;
- if (!getparfloat("cpm", &cpm)) cpm = 1000.0;
- if (!getparfloat("rom", &rom)) rom = 1000.0;
-
- vmess("dtm:%f dxm:%f cpm:%f rom:%f",dtm,dxm,cpm,rom);
-
- /* For a monopole source the scaling is (2.0*dt*cp*cp*rho)/(dx*dx) */
- for (iw=0; iw<ntfft; iw++){
- wavelet[iw] *= dt*(2.0*dtm*cpm*cpm*rom)/(dx*dx);
- }
- }
-
for (l=0; l<Nfoc; l++) {
- memset(&conv[0],0.0, sizeof(float)*ntfft*nxs*nys);
- convol(f1dsamp, &Gd[l*nxs*nys*ntfft], conv, nxs*nys, ntfft, dt, 0);
- for (i=0; i<nxs*nys; i++) {
- for (iw=0; iw<ntfft; iw++) {
- scaled[iw] += dt*scl*conv[i*ntfft+iw];
+ for (i=0; i<npos; i++) {
+ ix = ixpos[i];
+ iw = 0;
+ f1dsamp[i*ntfft+iw] = f1d[l*nxs*nys*ntfft+ix*ntfft+iw];
+ for (iw=1; iw<ntfft; iw++) {
+ f1dsamp[i*ntfft+iw] = f1d[l*nxs*nys*ntfft+ix*ntfft+ntfft-iw];
}
}
- ampest[l] = (wavelet[0]/scaled[0]);
- vmess("wavelet[0]=%f scaled[0]=%f",wavelet[0],scaled[0]);
- memset(&conv[0],0.0, sizeof(float)*ntfft*nxs*nys);
- memset(&wavelet[0],0.0, sizeof(float)*ntfft);
- memset(&scaled[0],0.0, sizeof(float)*ntfft);
+ deconv(&f1dsamp[0], &Gd[l*nxs*nys*ntfft], &est[l*nxs*nys*ntfft], nxs*nys, ntfft, dt, eps, reps, 0);
}
- free(wavelet);free(scaled);free(conv);free(f1dsamp);
+ free(f1dsamp);
return;
}
-long maxest3D(float *data, long nt)
-{
- float maxt;
- long it;
-
- maxt = data[0];
- for (it = 0; it < nt; it++) {
- if (fabs(data[it]) > fabs(maxt)) maxt=data[it];
- }
-
- return maxt;
-}
-
/**
* Calculates the time convolution of two arrays by
* transforming the arrayis to frequency domain,
@@ -305,4 +254,112 @@ void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long
for (it = 0 ; it < nsamout ; it++)
datout[ix*nsamout+it] = scl*data[ix*nsam+it];
}
+}
+
+/**
+* Calculates the time deconvolution of two arrays by
+* transforming the arrayis to frequency domain,
+* divides the arrays and transform back to time.
+*
+**/
+
+void deconv(float *data1, float *data2, float *decon, long nrec, long nsam,
+ float dt, float eps, float reps, long shift)
+{
+ long i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, *den, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2, maxden, leps;
+ complex *cdata1, *cdata2, *cdec, tmp;
+
+ optn = loptncr(nsam);
+ nfreq = optn/2+1;
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ cdec = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdec == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+ den = (float *)malloc(nfreq*nrec*sizeof(float));
+ if (den == NULL) verr("memory allocation error for rdata1");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
+ rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
+
+ /* apply deconvolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ n = nrec*nfreq;
+ maxden=0.0;
+ for (j = 0; j < n; j++) {
+ den[j] = *p2r**p2r + *p2i**p2i;
+ maxden = MAX(den[j], maxden);
+ p2r += 2;
+ p2i += 2;
+ }
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &cdec[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ leps = reps*maxden+eps;
+ for (j = 0; j < n; j++) {
+
+ if (fabs(*p2r)>=fabs(*p2i)) {
+ *qr = (*p2r**p1r+*p2i**p1i)/(den[j]+leps);
+ *qi = (*p2r**p1i-*p2i**p1r)/(den[j]+leps);
+ } else {
+ *qr = (*p1r**p2r+*p1i**p2i)/(den[j]+leps);
+ *qi = (*p1i**p2r-*p1r**p2i)/(den[j]+leps);
+ }
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+ free(den);
+
+ if (shift) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+ tau = dt*(nsam/2);
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = cdec[j*nfreq+i].r*cos(om*tau) + cdec[j*nfreq+i].i*sin(om*tau);
+ tmp.i = cdec[j*nfreq+i].i*cos(om*tau) - cdec[j*nfreq+i].r*sin(om*tau);
+ cdec[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdec[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,decon,nsam);
+
+ free(cdec);
+ free(rdata1);
+ free(rdata2);
+ return;
}
\ No newline at end of file
diff --git a/marchenko3D/getFileInfo3D.c b/marchenko3D/getFileInfo3D.c
index 13553faf2f1c4d53f0cd995e6e766e238f3d0209..cdf000077ef7ab7bb5fac24fe6b72c8ea0e28104 100644
--- a/marchenko3D/getFileInfo3D.c
+++ b/marchenko3D/getFileInfo3D.c
@@ -72,7 +72,7 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
/* check to find out number of traces in shot gather */
one_shot = 1;
- itrace = 0;
+ itrace = 1;
igy = 1;
fldr_shot = hdr.fldr;
sx_shot = hdr.sx;
@@ -86,7 +86,6 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
while (one_shot) {
nread = fread( trace, sizeof(float), hdr.ns, fp );
assert (nread == hdr.ns);
- itrace++;
if (hdr.gy != gy_end) {
gy_end = hdr.gy;
igy++;
@@ -95,13 +94,24 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread==0) break;
if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr) ) break;
+ itrace++;
}
if (itrace>1) {
*n2 = itrace/igy;
*n3 = igy;
- dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
- dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ if (*n2>1) {
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ }
+ else {
+ dxrcv = 1.0/scl;
+ }
+ if (*n3>1) {
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ }
+ else {
+ dyrcv = 1.0/scl;
+ }
*d2 = fabs(dxrcv)*scl;
*d3 = fabs(dyrcv)*scl;
if (NINT(dxrcv*1e3) != NINT(fabs(hdr.d2)*1e3)) {
@@ -113,7 +123,7 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
*n2 = MAX(hdr.trwf, 1);
*n3 = 1;
*d2 = hdr.d2;
- *d3 = 0.0;
+ *d3 = 1.0;
dxrcv = hdr.d2;
dyrcv = 0.0;
}
@@ -135,7 +145,6 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
dyrcv = *d3;
while (!end_of_file) {
- itrace = 0;
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread != TRCBYTES) { break; }
fldr_shot = hdr.fldr;
@@ -146,11 +155,10 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
gy_start = hdr.gy;
gy_end = hdr.gy;
- itrace = 0;
+ itrace = 1;
igy = 1;
while (one_shot) {
fseeko( fp, data_sz, SEEK_CUR );
- itrace++;
if (hdr.gx != gx_end) dxrcv = MIN(dxrcv,abs(hdr.gx-gx_end));
if (hdr.gy != gy_end) {
igy++;
@@ -165,15 +173,34 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
break;
}
if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
+ itrace++;
}
if (itrace>1) {
*n2 = itrace/igy;
*n3 = igy;
- dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
- dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ if (*n2>1) {
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ }
+ else {
+ dxrcv = 1.0/scl;
+ }
+ if (*n3>1) {
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ }
+ else {
+ dyrcv = 1.0/scl;
+ }
dxsrc = (float)(hdr.sx - sx_shot)*scl;
dysrc = (float)(hdr.sy - sy_shot)*scl;
}
+ else {
+ *n2 = MAX(hdr.trwf, 1);
+ *n3 = 1;
+ *d2 = hdr.d2;
+ *d3 = 1.0;
+ dxrcv = hdr.d2/scl;
+ dyrcv = 1.0/scl;
+ }
if (verbose>1) {
fprintf(stderr," . Scanning shot %li (%li) with %li traces dxrcv=%.2f dxsrc=%.2f %li %li dyrcv=%.2f dysrc=%.2f %li %li\n",sx_shot,igath,itrace,dxrcv*scl,dxsrc,gx_end,gx_start,dyrcv*scl,dysrc,gy_end,gy_start);
}
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index ace2b3179f7110a2d0884615c06fefa5ac7096bc..f1b615a9cbdc5d1429c53364ffe77e6761c0af9a 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -38,12 +38,12 @@ long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
long nw, long nw_low, long nshots, long nx, long ny, long ntfft, long mode, float scale, long verbose);
long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
long *xnx, long Nfoc, long nx, long ny, long ntfft, long mode, long *maxval, float *tinv, long hw, long verbose);
-// int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn,
-// float *zsyn, int *ixpos, int npos, int iter);
long unique_elements(float *arr, long len);
void name_ext(char *filename, char *extension);
+void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
+
void applyMute3D(float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nt, long *xrcvsyn, long npos, long shift);
long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
@@ -56,6 +56,9 @@ double wallclock_time(void);
void AmpEst3D(float *f1d, float *Gd, float *ampest, long Nfoc, long nxs, long nys, long ntfft, long *ixpos, long npos,
char *file_wav, float dx, float dy, float dt);
+void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
+ long *xnx, long Nfoc, long nx, long ny, long ntfft, long *maxval, float *tinv, long verbose);
+
void synthesisPositions3D(long nx, long ny, long nxs, long nys, long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc,
long *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys, long nshots, long nxsrc, long nysrc,
long *ixpos, long *npos, long reci, long verbose);
@@ -65,6 +68,8 @@ void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, long nx, l
float dx, float dy, long ntfft, long nw, long nw_low, long nw_high, long mode, long reci, long nshots, long nxsrc, long nysrc,
long *ixpos, long npos, double *tfft, long *isxcount, long *reci_xsrc, long *reci_xrcv, float *ixmask, long verbose);
+void imaging3D(float *Image, float *Gmin, float *f1plus, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
+
long linearsearch(long *array, size_t N, long value);
/*********************** self documentation **********************/
@@ -119,26 +124,27 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_out, *fp_f1plus, *fp_f1min;
- FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin;
+ FILE *fp_out, *fp_f1plus, *fp_f1min, *fp_imag;
+ FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin, *fp_amp;
long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
long reci, countmin, mode, n2out, n3out, verbose, ntfft;
long iter, niter, tracf, *muteW, ampest;
- long hw, smooth, above, shift, *ixpos, npos, ix;
+ long hw, smooth, above, shift, *ixpos, npos, ix, nzim, nxim, nyim;
long nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
float d1, d2, d3, f1, f2, f3, fxsb, fxse, fysb, fyse, ft, fx, fy, *xsyn, *ysyn, dxsrc, dysrc;
float *green, *f2p, *pmin, *G_d, dt, dx, dy, dxs, dys, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
- float xmin, xmax, ymin, ymax, scale, tsq, Q, f0;
- float *ixmask, *iymask, *ampscl, *Gd;
+ float xmin, xmax, ymin, ymax, scale, tsq, Q, f0, *tmpdata;
+ float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim, dyim, dxim;
complex *Refl, *Fop;
- char *file_tinv, *file_shot, *file_green, *file_iter, *file_wav;
+ char *file_tinv, *file_shot, *file_green, *file_iter, *file_imag, *file_ampscl;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
- segy *hdrs_out;
+ char *file_ray, *file_amp, *file_wav;
+ segy *hdrs_out, *hdrs_Nfoc;
initargs(argc, argv);
requestdoc(1);
@@ -148,6 +154,8 @@ int main (int argc, char **argv)
if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
+ if (!getparstring("file_ray", &file_ray)) file_ray = NULL;
+ if (!getparstring("file_amp", &file_amp)) file_amp = NULL;
if (!getparstring("file_f1plus", &file_f1plus)) file_f1plus = NULL;
if (!getparstring("file_f1min", &file_f1min)) file_f1min = NULL;
if (!getparstring("file_gplus", &file_gplus)) file_gplus = NULL;
@@ -157,6 +165,8 @@ int main (int argc, char **argv)
if (!getparstring("file_pmin", &file_pmin)) file_pmin = NULL;
if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
if (!getparstring("file_wav", &file_wav)) file_wav = NULL;
+ if (!getparstring("file_imag", &file_imag)) file_imag = NULL;
+ if (!getparstring("file_ampscl", &file_ampscl)) file_ampscl = NULL;
if (!getparlong("verbose", &verbose)) verbose = 0;
if (file_tinv == NULL && file_shot == NULL)
verr("file_tinv and file_shot cannot be both input pipe");
@@ -187,15 +197,28 @@ int main (int argc, char **argv)
/*================ Reading info about shot and initial operator sizes ================*/
ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
- ret = getFileInfo3D(file_tinv, &n1, &n2, &n3, &ngath, &d1, &d2, &d3, &f1, &f2, &f3, &scl, &ntraces);
- Nfoc = ngath;
- nxs = n2;
- nys = n3;
- nts = n1;
- dxs = d2;
- dys = d3;
- fxsb = f2;
- fysb = f3;
+ if (file_ray!=NULL) {
+ ret = getFileInfo3D(file_ray, &n2, &n1, &n3, &ngath, &d2, &d1, &d3, &f2, &f1, &f3, &scl, &ntraces);
+ Nfoc = ngath;
+ nxs = n2;
+ nys = n3;
+ nts = n1;
+ dxs = d2;
+ dys = d3;
+ fxsb = f2;
+ fysb = f3;
+ }
+ else {
+ ret = getFileInfo3D(file_tinv, &n1, &n2, &n3, &ngath, &d1, &d2, &d3, &f1, &f2, &f3, &scl, &ntraces);
+ Nfoc = ngath;
+ nxs = n2;
+ nys = n3;
+ nts = n1;
+ dxs = d2;
+ dys = d3;
+ fxsb = f2;
+ fysb = f3;
+ }
ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
ret = getFileInfo3D(file_shot, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
@@ -253,12 +276,25 @@ int main (int argc, char **argv)
/*================ Read and define mute window based on focusing operator(s) ================*/
/* G_d = p_0^+ = G_d (-t) ~ Tinv */
-
- mode=-1; /* apply complex conjugate to read in data */
- readTinvData3D(file_tinv, xrcvsyn, yrcvsyn, xsyn, ysyn, zsyn, xnxsyn, Nfoc,
- nxs, nys, ntfft, mode, muteW, G_d, hw, verbose);
+ if (file_ray!=NULL) {
+ makeWindow3D(file_ray, file_amp, file_wav, dt, xrcvsyn, yrcvsyn, xsyn, ysyn, zsyn,
+ xnxsyn, Nfoc, nx, ny, ntfft, muteW, G_d, verbose);
+ }
+ else {
+ mode=-1; /* apply complex conjugate to read in data */
+ readTinvData3D(file_tinv, xrcvsyn, yrcvsyn, xsyn, ysyn, zsyn, xnxsyn, Nfoc,
+ nxs, nys, ntfft, mode, muteW, G_d, hw, verbose);
+ }
/* reading data added zero's to the number of time samples to be the same as ntfft */
nts = ntfft;
+
+ /*Determine the shape of the focal positions*/
+ nzim = unique_elements(zsyn,Nfoc);
+ if (nzim>1) dzim = zsyn[1]-zsyn[0];
+ else dzim = 1.0;
+ nyim = unique_elements(ysyn,Nfoc);
+ nxim = unique_elements(xsyn,Nfoc);
+
/* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
@@ -619,7 +655,7 @@ int main (int argc, char **argv)
applyMute3D(green, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
/* compute upgoing Green's function G^+,- */
- if (file_gmin != NULL) {
+ if (file_gmin != NULL || file_imag!= NULL) {
Gmin = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
/* use f1+ as operator on R in frequency domain */
@@ -675,27 +711,92 @@ int main (int argc, char **argv)
/* Estimate the amplitude of the Marchenko Redatuming */
if (ampest>0) {
if (verbose>0) vmess("Estimating amplitude scaling");
+ // Create the first arrival data
Gd = (float *)calloc(Nfoc*nxs*nys*ntfft,sizeof(float));
memcpy(Gd,Gplus,sizeof(float)*Nfoc*nxs*nys*ntfft);
applyMute3D(Gd, muteW, smooth, 2, Nfoc, nxs*nys, nts, ixpos, npos, shift);
- ampscl = (float *)calloc(Nfoc,sizeof(float));
- AmpEst3D(G_d,Gd,ampscl,Nfoc,nxs,nys,ntfft,ixpos,npos,file_wav,dxs,dys,dt);
+ // Determine the amplitude
+ ampscl = (float *)calloc(nxs*nys*ntfft,sizeof(float));
+ tmpdata = (float *)calloc(nxs*nys*ntfft,sizeof(float));
+ // Scale all wavefields
for (l=0; l<Nfoc; l++) {
- for (j=0; j<nxs*nys*nts; j++) {
- green[l*nxs*nts+j] *= ampscl[l];
- if (file_gplus != NULL) Gplus[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_gmin != NULL) Gmin[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f2 != NULL) f2p[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_pmin != NULL) pmin[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f1plus != NULL) f1plus[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f1min != NULL) f1min[l*nxs*nys*nts+j] *= ampscl[l];
- }
- if (verbose>1) vmess("Amplitude of focal position %li is equal to %.3e",l,ampscl[l]);
+ AmpEst3D(&G_d[l*nxs*nys*ntfft],&Gd[l*nxs*nys*ntfft],ampscl,1,nxs,nys,ntfft,ixpos,npos,file_wav,dxs,dys,dt);
+ for (j=0; j<nxs*nys*nts; j++) {
+ tmpdata[j] = green[l*nxs*nys*nts+j];
+ }
+ convol(tmpdata, ampscl, &green[l*nxs*nys*nts], nxs*nys, ntfft, dt, 0);
+ if (file_gplus != NULL) {
+ for (j=0; j<nxs*nys*nts; j++) {
+ tmpdata[j] = Gplus[l*nxs*nys*nts+j];
+ }
+ convol(tmpdata, ampscl, &Gplus[l*nxs*nys*nts], nxs*nys, ntfft, dt, 0);
+ }
+ if (file_gmin != NULL || file_imag!=NULL) {
+ for (j=0; j<nxs*nys*nts; j++) {
+ tmpdata[j] = Gmin[l*nxs*nys*nts+j];
+ }
+ convol(tmpdata, ampscl, &Gmin[l*nxs*nys*nts], nxs*nys, ntfft, dt, 0);
+ }
+ //if (verbose>4) vmess("Amplitude of focal position %li is equal to %.3e",l,ampscl[l]);
}
+ free(tmpdata);
+ // if (file_ampscl!=NULL) { //Write the estimation of the amplitude to file
+ // hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
+ // for (l=0; l<nyim; l++){
+ // for (j=0; j<nxim; j++){
+ // hdrs_Nfoc[l*nxim+j].ns = nzim;
+ // hdrs_Nfoc[l*nxim+j].sx = xsyn[j];
+ // hdrs_Nfoc[l*nxim+j].sy = ysyn[l];
+ // hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l];
+ // hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
+ // hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
+ // hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
+ // hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
+ // }
+ // }
+ // // Write the data
+ // fp_amp = fopen(file_ampscl, "w+");
+ // if (fp_amp==NULL) verr("error on creating output file %s", file_ampscl);
+ // ret = writeData3D(fp_amp, (float *)&scl[0], hdrs_Nfoc, nzim, nxim*nyim);
+ // if (ret < 0 ) verr("error on writing output file.");
+ // fclose(fp_amp);
+ // free(hdrs_Nfoc);
+ // free(ampscl);
+ // }
free(Gd);
if (file_gplus == NULL) free(Gplus);
}
+ /* Apply imaging*/
+ if (file_imag!=NULL) {
+ // Determine Image
+ Image = (float *)calloc(Nfoc,sizeof(float));
+ imaging3D(Image, Gmin, f1plus, nxs, nys, ntfft, dxs, dys, dt, Nfoc, verbose);
+ if (file_gmin==NULL) free(Gmin);
+ // Set headers
+ hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
+ for (l=0; l<nyim; l++){
+ for (j=0; j<nxim; j++){
+ hdrs_Nfoc[l*nxim+j].ns = nzim;
+ hdrs_Nfoc[l*nxim+j].sx = xsyn[j];
+ hdrs_Nfoc[l*nxim+j].sy = ysyn[l];
+ hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l];
+ hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
+ hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
+ hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
+ hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
+ }
+ }
+ // Write out image
+ fp_imag = fopen(file_imag, "w+");
+ if (fp_imag==NULL) verr("error on creating output file %s", file_imag);
+ ret = writeData3D(fp_imag, (float *)&Image[0], hdrs_Nfoc, nzim, nxim*nyim);
+ if (ret < 0 ) verr("error on writing output file.");
+ fclose(fp_imag);
+ free(hdrs_Nfoc);
+ free(Image);
+ }
+
t2 = wallclock_time();
if (verbose) {
vmess("Total CPU-time marchenko = %.3f", t2-t0);
diff --git a/raytime/raytime.c b/raytime/raytime.c
index 532f9c1e692c1a9bf99594caa47d1622580c409a..414673e43627ba094f6f70bd2b34f64c08907ec1 100644
--- a/raytime/raytime.c
+++ b/raytime/raytime.c
@@ -262,7 +262,7 @@ int main(int argc, char **argv)
memset(&hdr,0,sizeof(hdr));
hdr.scalco = -1000;
hdr.scalel = -1000;
- hdr.trid = 1;
+ hdr.trid = 0;
t1=wallclock_time();
tinit = t1-t0;
diff --git a/utils/green3D b/utils/green3D
index 3513a0edb3c0096f4a53c9e3fa73a5474de53bf5..c2f71642dbc9a737f12146aaac155419fe32ad9c 100755
Binary files a/utils/green3D and b/utils/green3D differ