diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 8367a6d18441993c9e2e55ed9c2764dbecad3178..303417a0bd011ab60240be798e91db9ea1732c27 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -22,7 +22,7 @@ SRCT = marchenko3D.c \
imaging3D.c \
homogeneousg3D.c \
readSnapData3D.c \
- writeDataIter.c \
+ writeDataIter3D.c \
wallclock_time.c \
name_ext.c \
verbosepkg.c \
diff --git a/marchenko3D/homogeneousg3D.c b/marchenko3D/homogeneousg3D.c
index 8a09ecb483f3ce1e1a0bbd0802970c93bcd616f9..e5b2cfd497565f04da3f5d1c7175198febd9385f 100644
--- a/marchenko3D/homogeneousg3D.c
+++ b/marchenko3D/homogeneousg3D.c
@@ -32,7 +32,8 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt);
void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt);
-void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose)
+void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx, long ny, long nt, float dx, float dy,
+ float dt, long Nfoc, long *sx, long *sy, long *sz, long verbose)
{
char *file_inp;
long i, j, k, l, ret, scheme, count=0, n_source;
@@ -73,6 +74,16 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx,
readSnapData3D(file_inp, input, hdr_inp, n_source, nx, ny, nt, 0, nx, 0, ny, 0, nt);
zsrc = labs(hdr_inp[0].selev);
+ for (l = 0; l < nx*ny; l++) {
+ sx[l] = hdr_inp[0].sx;
+ sy[l] = hdr_inp[0].sy;
+ sz[l] = hdr_inp[0].sdepth;
+ }
+ for (l = 0; l < n_source; l++) {
+ sx[l] = hdr_inp[l*nx*ny].sx;
+ sy[l] = hdr_inp[l*nx*ny].sy;
+ sz[l] = hdr_inp[l*nx*ny].sdepth;
+ }
if (scheme==1) { // Scale the Green's functions if the classical scheme is used
if (verbose) vmess("Classical Homogeneous Green's function retrieval");
diff --git a/marchenko3D/makeWindow3D.c b/marchenko3D/makeWindow3D.c
index 5930c5032ae96c25d81a3f74167d0f4f83e4954d..2809e630aa62d37dd65ed02db0dbcca3f85a35a0 100644
--- a/marchenko3D/makeWindow3D.c
+++ b/marchenko3D/makeWindow3D.c
@@ -73,46 +73,49 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
/* Defining mute window using raytimes */
vmess("Using raytime for mutewindow");
- hdrs_mute = (segy *) calloc(Nfoc,sizeof(segy));
+ hdrs_mute = (segy *) calloc(Nfoc*ny,sizeof(segy));
fp = fopen( file_ray, "r" );
if ( fp == NULL ) {
perror("Error opening file containing ray");
}
fclose(fp);
- readSnapData3D(file_ray, timeval, hdrs_mute, Nfoc, 1, 1, nx, 0, 1, 0, 1, 0, nx);
- // for (is=0; is<Nfoc; is++) {
- // readData3D(fp, &timeval[is*nx*ny], &hdrs_mute[is], nx);
- // }
+ readSnapData3D(file_ray, timeval, hdrs_mute, Nfoc, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+
+ /* Set the scaling for gy */
+ if (hdrs_mute[0].scalel < 0) {
+ scl = 1.0/((float)(labs(hdrs_mute[0].scalel)));
+ }
+ else if (hdrs_mute[0].scalel > 0) {
+ scl = ((float)hdrs_mute[0].scalel);
+ }
+ else {
+ scl = 1.0;
+ }
/*Check whether the amplitude is also used*/
if (file_amp != NULL) {
vmess("Using ray-amplitudes");
- hdrs_amp = (segy *) calloc(Nfoc,sizeof(segy));
+ hdrs_amp = (segy *) calloc(Nfoc*ny,sizeof(segy));
fp = fopen( file_amp, "r" );
if ( fp == NULL ) {
perror("Error opening file containing ray-amplitude");
}
fclose(fp);
- readSnapData3D(file_amp, amp, hdrs_amp, Nfoc, 1, 1, nx, 0, 1, 0, 1, 0, nx);
- // for (is=0; is<Nfoc; is++) {
- // readData3D(fp, &[is*nx*ny], &hdrs_amp[is], nx);
- // }
+ readSnapData3D(file_amp, amp, hdrs_amp, Nfoc, 1, ny, nx, 0, 1, 0, ny, 0, nx);
}
/*Define source and receiver locations from the raytime*/
for (is=0; is<Nfoc; is++) {
for (iy=0; iy<ny; iy++) {
for (ix=0; ix<nx; ix++) {
- xrcv[is*nx*ny+iy*nx+ix] = (hdrs_mute[is].f1 + hdrs_mute[is].d1*((float)ix));
- yrcv[is*nx*ny+iy*nx+ix] = hdrs_mute[is].gy;
+ xrcv[is*nx*ny+iy*nx+ix] = (hdrs_mute[is*ny].f1 + hdrs_mute[is*ny].d1*((float)ix));
+ yrcv[is*nx*ny+iy*nx+ix] = ((float)hdrs_mute[is*ny].gy)*scl;
}
}
- xnx[is]=hdrs_mute[is].ns;
- if (hdrs_mute[is].scalco < 0) scl=-1.0/hdrs_mute[is].scalco;
- else scl=hdrs_mute[is].scalco;
- xsrc[is] = hdrs_mute[is].sx*scl;
- ysrc[is] = hdrs_mute[is].sy*scl;
- zsrc[is] = hdrs_mute[is].sdepth*scl;
+ xnx[is]=hdrs_mute[is*ny].ns;
+ xsrc[is] = ((float)hdrs_mute[is*ny].sx)*scl;
+ ysrc[is] = ((float)hdrs_mute[is*ny].sy)*scl;
+ zsrc[is] = ((float)hdrs_mute[is*ny].sdepth)*scl;
}
@@ -125,8 +128,8 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
if (file_wav!=NULL) { /*Apply the wavelet to create a first arrival*/
if (file_amp != NULL) {
for (ig=0; ig<nfreq; ig++) {
- cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]*amp[j*ny*nx+l*nx+i]);
- cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]*amp[j*ny*nx+l*nx+i]);
+ cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]);
+ cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]);
}
}
else { /*Use the raytime only to determine the mutewindow*/
@@ -146,7 +149,8 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
}
}
}
-
+ free(timeval); free(hdrs_mute);
+ if (file_amp!=NULL) free(amp); free(hdrs_amp);
return;
}
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 5f892f7f125097761ff727947debb6cafcc2fde4..f9b55d1d953f248556298eeac0ed715176143baf 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -72,7 +72,8 @@ long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2,
void imaging3D(float *Image, float *Gmin, float *f1plus, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
-void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
+void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx, long ny, long nt, float dx, float dy,
+ float dt, long Nfoc, long *sx, long *sy, long *sz, long verbose);
long linearsearch(long *array, size_t N, long value);
@@ -155,7 +156,7 @@ int main (int argc, char **argv)
FILE *fp_out, *fp_f1plus, *fp_f1min, *fp_imag, *fp_homg;
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 size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad, *sx, *sy, *sz;
long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
long reci, countmin, mode, n2out, n3out, verbose, ntfft;
long iter, niter, tracf, *muteW, *tsynW, ampest, plane_wave;
@@ -249,7 +250,6 @@ int main (int argc, char **argv)
dys = d3;
fxsb = f2;
fysb = f3;
- vmess("nxs:%li dxs:%.3f",nxs,dxs);
}
else {
ret = getFileInfo3D(file_tinv, &n1, &n2, &n3, &ngath, &d1, &d2, &d3, &f1, &f2, &f3, &scl, &ntraces);
@@ -534,6 +534,7 @@ int main (int argc, char **argv)
if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
}
+
/*================ initializations ================*/
if (reci) {
@@ -751,7 +752,6 @@ int main (int argc, char **argv)
free(Ni);
free(Nig);
- if (ampest < 1) free(G_d);
/* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
for (l = 0; l < Nfoc; l++) {
@@ -945,9 +945,14 @@ int main (int argc, char **argv)
/* Determine homogeneous Green's function*/
if (file_homg!=NULL) {
+ // Allocate the headers for the source info
+ sx = (long *)calloc(nxs*nys,sizeof(long));
+ sy = (long *)calloc(nxs*nys,sizeof(long));
+ sz = (long *)calloc(nxs*nys,sizeof(long));
+
// Determine Image
HomG = (float *)calloc(Nfoc*ntfft,sizeof(float));
- homogeneousg3D(HomG, green, f2p, zsyn, nxs, nys, ntfft, dxs, dys, dt, Nfoc, verbose);
+ homogeneousg3D(HomG, green, f2p, zsyn, nxs, nys, ntfft, dxs, dys, dt, Nfoc, sx, sy, sz, verbose);
fp_homg = fopen(file_homg, "w+");
if (fp_homg==NULL) verr("error on creating output file %s", file_homg);
@@ -964,11 +969,12 @@ int main (int argc, char **argv)
hdrs_Nfoc[l*nxim+j].trid = 2;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
- hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sx = sx[l*nxim+j];
+ hdrs_Nfoc[l*nxim+j].sy = sy[l*nxim+j];
hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sdepth = sz[l*nxim+j];
+ hdrs_Nfoc[l*nxim+j].selev = -sz[l*nxim+j];
hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
hdrs_Nfoc[l*nxim+j].d1 = dzim;
@@ -1050,8 +1056,8 @@ int main (int argc, char **argv)
hdrs_out[k*n2+i].sy = NINT(ysyn[l]*1000);
hdrs_out[k*n2+i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
hdrs_out[k*n2+i].tracf = tracf++;
- hdrs_out[k*n2+i].selev = NINT(zsyn[l]*1000);
- hdrs_out[k*n2+i].sdepth = NINT(-zsyn[l]*1000);
+ hdrs_out[k*n2+i].selev = NINT(-zsyn[l]*1000);
+ hdrs_out[k*n2+i].sdepth = NINT(zsyn[l]*1000);
hdrs_out[k*n2+i].f1 = f1;
}
}
diff --git a/marchenko3D/writeDataIter3D.c b/marchenko3D/writeDataIter3D.c
index 09a7cd83680f07ada4e43cce99cb2f103206c389..8d63f7dd0109fd69661239f7e7a58cceaef464fb 100644
--- a/marchenko3D/writeDataIter3D.c
+++ b/marchenko3D/writeDataIter3D.c
@@ -20,7 +20,7 @@ long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2,
{
FILE *fp_iter;
size_t nwrite;
- int i, l, j, ip, ret, tracf, size, ix, iy;
+ int i, l, j, ret, tracf, size, ix, iy;
char number[16], filename[1024];
float *trace;
@@ -31,10 +31,11 @@ long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2,
fp_iter = fopen(filename, "w+");
if (fp_iter==NULL) verr("error on creating output file %s", filename);
size=n1*n2*n3;
+ tracf=1;
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
- ix = ixpos[ip]; /* select proper position */
- iy = iypos[ip];
+ ix = ixpos[i]; /* select proper position */
+ iy = iypos[i];
hdrs[i].fldr = l+1;
hdrs[i].sx = NINT(xsyn[l]*1000);
hdrs[i].sy = NINT(ysyn[l]*1000);
@@ -61,7 +62,6 @@ long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2,
assert(nwrite == TRCBYTES);
nwrite = fwrite(trace, sizeof(float), n1, fp_iter);
assert (nwrite == n1);
- ip++;
}
}
ret = fclose(fp_iter);
diff --git a/raytime3d/getParameters3d.c b/raytime3d/getParameters3d.c
index 13fc3737bbdf5393724e44e31f759dccc41fc6ad..24c70664c7b9ed01245e4cf9f2444a431c146397 100644
--- a/raytime3d/getParameters3d.c
+++ b/raytime3d/getParameters3d.c
@@ -41,7 +41,7 @@ long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPa
float rsrc, oxsrc, ozsrc, oysrc, dphisrc, ncsrc;
size_t nsamp;
long nxsrc, nzsrc, nysrc;
- long is;
+ long is, shotnxyz, rcvn;
char *src_positions, *file_cp=NULL;
if (!getparlong("verbose",&verbose)) verbose=0;
@@ -322,6 +322,26 @@ long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPa
recvPar3D(rec, sub_x0, sub_y0, sub_z0, dx, dy, dz, nx, ny, nz);
+ /* Calculate the maximum radius */
+ shotnxyz = shot->nx*shot->ny*shot->nz-1;
+ rcvn = rec->n-1;
+ shot->maxrad = labs(rec->x[0]-shot->x[0]);
+ if (shot->maxrad < labs(rec->x[rcvn]-shot->x[0])) shot->maxrad = labs(rec->x[rcvn]-shot->x[0]);
+ if (shot->maxrad < labs(rec->x[0]-shot->x[shotnxyz])) shot->maxrad = labs(rec->x[0]-shot->x[shotnxyz]);
+ if (shot->maxrad < labs(rec->x[rcvn]-shot->x[shotnxyz])) shot->maxrad = labs(rec->x[rcvn]-shot->x[shotnxyz]);
+
+ if (shot->maxrad < labs(rec->y[0]-shot->y[0])) shot->maxrad = labs(rec->y[0]-shot->y[0]);
+ if (shot->maxrad < labs(rec->y[rcvn]-shot->y[0])) shot->maxrad = labs(rec->y[rcvn]-shot->y[0]);
+ if (shot->maxrad < labs(rec->y[0]-shot->y[shotnxyz])) shot->maxrad = labs(rec->y[0]-shot->y[shotnxyz]);
+ if (shot->maxrad < labs(rec->y[rcvn]-shot->y[shotnxyz])) shot->maxrad = labs(rec->y[rcvn]-shot->y[shotnxyz]);
+
+ if (shot->maxrad < labs(rec->z[0]-shot->z[0])) shot->maxrad = labs(rec->z[0]-shot->z[0]);
+ if (shot->maxrad < labs(rec->z[rcvn]-shot->z[0])) shot->maxrad = labs(rec->z[rcvn]-shot->z[0]);
+ if (shot->maxrad < labs(rec->z[0]-shot->z[shotnxyz])) shot->maxrad = labs(rec->z[0]-shot->z[shotnxyz]);
+ if (shot->maxrad < labs(rec->z[rcvn]-shot->z[shotnxyz])) shot->maxrad = labs(rec->z[rcvn]-shot->z[shotnxyz]);
+
+ shot->maxrad++;
+
if (verbose) {
if (rec->n) {
dyrcv = rec->yr[rec->nx]-rec->yr[0];
@@ -341,6 +361,7 @@ long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPa
vmess("izmin = %li izmax = %li ", rec->z[0], rec->z[rec->n-1]);
vmess("ixmin = %li ixmax = %li ", rec->x[0], rec->x[rec->n-1]);
vmess("iymin = %li iymax = %li ", rec->y[0], rec->y[rec->n-1]);
+ vmess("maxrad = %li",shot->maxrad);
fprintf(stderr,"\n");
}
else {
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
index 7d5e65a0ead722ade4f4634ffe5bdfe991006eca..f93bb963090d3480524bcd46d24524ae2b5f6f7e 100644
--- a/raytime3d/raytime3d.c
+++ b/raytime3d/raytime3d.c
@@ -189,9 +189,6 @@ void main(int argc, char *argv[])
if (slow0 == NULL) verr("Out of memory for slow0 array!");
readModel3d(mod.file_cp, slow0, nz, nx, ny, h, verbose);
-
- if (verbose) vmess("h = %.2f nx = %li nz = %li ny = %li", h, nx, nz, ny);
-
}
else {
if(!getparfloat("c",&c)) verr("c not defined");
@@ -287,8 +284,8 @@ private (is,time0,ampl,nrx,nry,nrz,nr,cp_average,i,j,k,ix,iy,iz,hdrs,tmpdata,nwr
#pragma omp flush(writer)
}
- if (verbose>2) vmess("Writing src %li of %li sources",is+1,shot.n);
- if (verbose) vmess("xsrc[%li]=%f ysrc[%li]=%f zsrc[%li]=%f",shot.x[is],shot.xs[is],shot.y[is],shot.ys[is],shot.z[is],shot.zs[is]);
+ if (verbose) vmess("Writing src %li of %li sources",is+1,shot.n);
+ if (verbose>1) vmess("xsrc[%li]=%f ysrc[%li]=%f zsrc[%li]=%f",shot.x[is],shot.xs[is],shot.y[is],shot.ys[is],shot.z[is],shot.zs[is]);
hdrs = (segy *) calloc(1,sizeof(segy));
tmpdata = (float *)malloc((rec.nx)*sizeof(float));
@@ -303,7 +300,7 @@ private (is,time0,ampl,nrx,nry,nrz,nr,cp_average,i,j,k,ix,iy,iz,hdrs,tmpdata,nwr
hdrs[0].sy = (long)(f3+(shot.y[is]-1)*d3)*1000;
hdrs[0].sdepth = (long)(f1+(shot.z[is]-1)*d1)*1000;
hdrs[0].selev = -(long)(f1+(shot.z[is]-1)*d1)*1000;
- hdrs[0].gy = mod.y0+rec.yr[i*rec.nx];
+ hdrs[0].gy = (long)(mod.y0+rec.yr[i*rec.nx])*1000;
hdrs[0].ns = rec.nx;
hdrs[0].ntr = rec.ny*shot.n;
hdrs[0].trwf = rec.ny*shot.n;
diff --git a/raytime3d/raytime3d.h b/raytime3d/raytime3d.h
index dbf479ac40df596daebffefb34759c8922b19538..14867d9687f70841a2e460043d533ae31df8757f 100644
--- a/raytime3d/raytime3d.h
+++ b/raytime3d/raytime3d.h
@@ -118,6 +118,7 @@ typedef struct _shotPar { /* Shot Parameters */
float *zs;
float *xs;
float *ys;
+ long maxrad;
} shotPar;
typedef struct _raypar { /* ray-tracing parameters */
diff --git a/raytime3d/src3D.c b/raytime3d/src3D.c
index 369cd993233172b3a0adb14fc2bd525394eef535..efdd2f44562368277bd78fd7b79ef08b91418fa8 100644
--- a/raytime3d/src3D.c
+++ b/raytime3d/src3D.c
@@ -78,13 +78,13 @@ void src3D(float *time0, float *slow0, long nz, long nx, long ny, float h, float
dvz = (1/s0(xx,yy,zz)-1/s0(xs,ys,zs))/(zz-zs);
dv = fabs(dvz);
if (dv == 0.) {
- t0(xx,yy,zz) = s0(xs,ys,zs)*DIST(fxs,fys,fzs,xx,yy,zz);
+ t0(xx,yy,zz) = s0(xs,ys,zs)*DIST(xs,ys,zs,xx,yy,zz);
continue;
}
rzc = -v0/dv;
- rx = h*(xx - fxs);
- ry = h*(yy - fys);
- rz = h*(zz - fzs);
+ rx = h*(xx - xs);
+ ry = h*(yy - ys);
+ rz = h*(zz - zs);
rz1 = rz*dvz/dv;
rxy1 = sqrt(rx*rx+ry*ry+rz*rz-rz1*rz1);
if (rxy1<=h/1.e6)
diff --git a/raytime3d/vidale3d.c b/raytime3d/vidale3d.c
index d7e560baa6810a5c645388facb1b4a2b0724c5b5..c7107bd0cdda9783463a0e9e2104b8a09a06d7ca 100644
--- a/raytime3d/vidale3d.c
+++ b/raytime3d/vidale3d.c
@@ -1264,42 +1264,42 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
else {
head=0;
if (headw[1]>0) {
- if(verbose) vmess("Head waves found on left: %li",headw[1]);
+ if(verbose>3) vmess("Head waves found on left: %li",headw[1]);
if (headw[1]>head) {
head = headw[1];
srcwall = 1;
}
}
if (headw[2]>0) {
- if(verbose) vmess("Head waves found on right: %li",headw[2]);
+ if(verbose>3) vmess("Head waves found on right: %li",headw[2]);
if (headw[2]>head) {
head = headw[2];
srcwall = 2;
}
}
if (headw[3]>0) {
- if(verbose) vmess("Head waves found on front: %li",headw[3]);
+ if(verbose>3) vmess("Head waves found on front: %li",headw[3]);
if (headw[3]>head) {
head = headw[3];
srcwall = 3;
}
}
if (headw[4]>0) {
- if(verbose) vmess("Head waves found on back: %li",headw[4]);
+ if(verbose>3) vmess("Head waves found on back: %li",headw[4]);
if (headw[4]>head) {
head = headw[4];
srcwall = 4;
}
}
if (headw[5]>0) {
- if(verbose) vmess("Head waves found on top: %li",headw[5]);
+ if(verbose>3) vmess("Head waves found on top: %li",headw[5]);
if (headw[5]>head) {
head = headw[5];
srcwall = 5;
}
}
if (headw[6]>0) {
- if(verbose) vmess("Head waves found on bottom: %li",headw[6]);
+ if(verbose>3) vmess("Head waves found on bottom: %li",headw[6]);
if (headw[6]>head) {
head = headw[6];
srcwall = 6;
@@ -1332,7 +1332,7 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
vmess("RESTART at bottom side of model"); }
else { z1= -1; dz1=0; }
if (reverse == 0)
- vwarn("RESTART CANCELLED by choice of input parameter `reverse`");
+ if (verbose>3) vwarn("RESTART CANCELLED by choice of input parameter `reverse`");
}
reverse--;
diff --git a/utils/combine.c b/utils/combine.c
index 45fb07d4d48698c796390148ed3d064c5da6c6b8..c08792ec069405fce004eaf74471396895c411c6 100755
--- a/utils/combine.c
+++ b/utils/combine.c
@@ -121,7 +121,12 @@ int main (int argc, char **argv)
sprintf(fins,"%li",numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
nt = 1;
- getFileInfo3D(fin2, &nz, &nx, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
+ if (transpose==0) {
+ getFileInfo3D(fin2, &nz, &nx, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
+ }
+ else {
+ getFileInfo3D(fin2, &nx, &nz, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
+ }
nxyz = nx*ny*nzs*nz;
@@ -145,20 +150,34 @@ int main (int argc, char **argv)
* Combine the separate files
*----------------------------------------------------------------------------*/
for (is = 0; is < nzs; is++) {
- if (verbose) vmess("Combining file %li out of %li",iz+1,nzs);
- sprintf(fins,"%li",iz*dnumb+numb);
+ if (verbose) vmess("Combining file %li out of %li",is+1,nzs);
+ sprintf(fins,"%li",is*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin2, "r");
if (fp_in == NULL) {
verr("Error opening file");
}
fclose(fp_in);
- readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
- for (it = 0; it < nt; it++) {
- for (iy = 0; iy < ny; iy++) {
- for (ix = 0; ix < nx; ix++) {
- for (iz = 0; iz < nz; iz++) {
- outdata[it*ny*nx*nz*nzs+iy*nx*nz*nzs+ix*nz*nzs+iz+is] = indata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
+ if (transpose==0) {
+ readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
+ for (it = 0; it < nt; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*ny*nx*nz*nzs+iy*nx*nz*nzs+ix*nz*nzs+iz+is] = indata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
+ }
+ }
+ }
+ }
+ }
+ else {
+ readSnapData3D(fin2, indata, hdr_in, nt, nz, ny, nx, 0, nz, 0, ny, 0, nx);
+ for (it = 0; it < nt; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*ny*nx*nz*nzs+iy*nx*nz*nzs+ix*nz*nzs+iz+is] = indata[it*ny*nz*nx+iy*nz*nx+iz*nx+ix];
+ }
}
}
}
diff --git a/utils/combine_induced.c b/utils/combine_induced.c
index 4b766d3f2c64b2e732e9268443d0e3354633d45b..bdb37148f5f4fe00b96cff1597113798253685bc 100755
--- a/utils/combine_induced.c
+++ b/utils/combine_induced.c
@@ -51,10 +51,10 @@ int main (int argc, char **argv)
{
FILE *fp_in, *fp_out;
char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
- float *indata, *outdata, fz, fy, fx, shift, dtshift, dt_time;
- float dt, dy, dx, t0, x0, y0, sclsxgx, dt2, dy2, dx2, t02, x02, y02, sclsxgx2, dxrcv, dyrcv, dzrcv;
- long nshots, nt, ny, nx, ntraces, nshots2, nt2, ny2, nx2, ntraces2, ix, iy, it, is, iz, pos, file_det, nxs, nys, nzs;
- long numb, dnumb, ret, nzmax, verbose, nshot_out, ishift, nshift;
+ float *indata, *outdata, shift, dtshift, dt_time;
+ float dt, dz, dy, dx, t0, x0, y0, z0, scl, dxrcv, dyrcv, dzrcv;
+ long nt, nz, ny, nx, nxyz, ntr, ix, iy, it, is, iz, pos, file_det, nxs, nys, nzs;
+ long numb, dnumb, ret, nzmax, verbose, nt_out, ishift, nshift, *sx, *sy, *sz;
segy *hdr_in, *hdr_bin, *hdr_out;
initargs(argc, argv);
@@ -123,72 +123,77 @@ int main (int argc, char **argv)
*----------------------------------------------------------------------------*/
sprintf(fins,"%li",numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
- nshots = 0;
- getFileInfo3D(fin2, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &sclsxgx, &ntraces);
+ nt = 0;
+ getFileInfo3D(fin2, &nz, &nx, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
- sprintf(fins,"%li",numb+dnumb);
- sprintf(fin2,"%s%s%s",fbegin,fins,fend);
- nshots = 0;
- getFileInfo3D(fin2, &nt2, &nx2, &ny2, &nshots2, &dt2, &dx2, &dy2, &t02, &x02, &y02, &sclsxgx2, &ntraces2);
-
- dxrcv=dx*1000;
- dyrcv=dy*1000;
- dzrcv=t02-t0;
+ nxyz = nx*ny*nz;
- if (nshots==0) nshots=1;
- nxs = ntraces;
+ if (verbose) {
+ vmess("number of time samples: %li", nt);
+ vmess("Number of virtual receivers: %li, x: %li, y: %li, z: %li",nxyz,nx,ny,nz);
+ vmess("Starting distance for x: %.3f, y: %.3f, z: %.3f",x0,y0,z0);
+ vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
+ vmess("Number of virtual sources: %li",nzs);
+ }
/*----------------------------------------------------------------------------*
* Read in a single file to determine if the header values match
* and allocate the data
*----------------------------------------------------------------------------*/
- hdr_bin = (segy *)calloc(nxs,sizeof(segy));
- indata = (float *)calloc(nxs*nt,sizeof(float));
+ hdr_in = (segy *)calloc(nx*ny,sizeof(segy));
+ indata = (float *)calloc(nxyz*nt,sizeof(float));
+
+ readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
- readSnapData3D(fin2, &indata[0], &hdr_bin[0], nshots, nxs, ny, nt, 0, nxs, 0, ny, 0, nt);
- nshots = hdr_bin[nxs-1].fldr;
- nxs = hdr_bin[nxs-1].tracf;
+ dt = ((float)hdr_in[0].dt)/1E6;
+ sx = (long *)calloc(nx*ny,sizeof(float));
+ sy = (long *)calloc(nx*ny,sizeof(float));
+ sz = (long *)calloc(nx*ny,sizeof(float));
- nshot_out = nshots/2;
- free(indata);
+ for (ix = 0; ix < nx*ny; ix++) {
+ sx[ix] = hdr_in[0].sx;
+ sy[ix] = hdr_in[0].sy;
+ sz[ix] = hdr_in[0].sdepth;
+ }
- hdr_out = (segy *)calloc(nshot_out*nxs,sizeof(segy));
- outdata = (float *)calloc(nshot_out*nxs*nt,sizeof(float));
+ nt_out = nt/2+1;
+ free(indata); free(hdr_in);
- /*----------------------------------------------------------------------------*
- * Write out the file info in case of verbose
- *----------------------------------------------------------------------------*/
- if (verbose) vmess("Number of virtual receivers: %li, nx=%li, ny=%li, nz=%li",nx*ny*nt,nx,ny,nt);
+ hdr_out = (segy *)calloc(nx*ny,sizeof(segy));
+ outdata = (float *)calloc(nt_out*nxyz,sizeof(float));
/*----------------------------------------------------------------------------*
* Parallel loop for reading in and combining the various shots
*----------------------------------------------------------------------------*/
#pragma omp parallel default(shared) \
- private(indata,iz,hdr_in,fins,fin2,fp_in,is,ix,iy,it,ishift)
+ private(indata,hdr_in,fins,fin2,fp_in,is,ix,iy,iz,it,ishift)
{
- indata = (float *)calloc(ntraces*nt,sizeof(float));
- hdr_in = (segy *)calloc(ntraces,sizeof(segy));
+ indata = (float *)calloc(nxyz*nt,sizeof(float));
+ hdr_in = (segy *)calloc(nx*ny*nt,sizeof(segy));
#pragma omp for
- for (iz = 0; iz < nzs; iz++) {
- if (verbose) vmess("Depth:%li out of %li",iz+1,nzs);
- sprintf(fins,"%li",iz*dnumb+numb);
+ for (is = 0; is < nzs; is++) {
+ if (verbose) vmess("Depth:%li out of %li",is+1,nzs);
+ sprintf(fins,"%li",is*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin2, "r");
if (fp_in == NULL) {
verr("Error opening file");
}
fclose(fp_in);
- readSnapData3D(fin2, &indata[0], &hdr_in[0], nshots, nxs, ny, nt, 0, nxs, 0, ny, 0, nt);
- if (iz==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
- if (iz==1) dzrcv=hdr_in[0].f1-fz;
+ readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
+ sx[is] = hdr_in[0].sx;
+ sy[is] = hdr_in[0].sy;
+ sz[is] = hdr_in[0].sdepth;
- ishift = nshift*iz;
+ ishift = nshift*is;
if (verbose) vmess("Shifting %li timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)iz)));
- for (is = ishift; is < nshot_out; is++) {
- for (ix = 0; ix < nxs; ix++) {
- for (it = 0; it < nt; it++) {
- outdata[is*nxs*nt+ix*nt+it] += indata[(is-ishift+(nshots/2))*nxs*nt+ix*nt+it];
+ for (it = ishift; it < nt_out; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*nxyz+iy*nx*nz+ix*nz+iz] += indata[(it-ishift+(nt/2))*nxyz+iy*nx*nz+ix*nz+iz];
+ }
}
}
}
@@ -201,26 +206,31 @@ int main (int argc, char **argv)
*----------------------------------------------------------------------------*/
fp_out = fopen(fout, "w+");
- for (is = 0; is < nshot_out; is++) {
- for (ix = 0; ix < nxs; ix++) {
- hdr_out[ix].fldr = is+1;
- hdr_out[ix].tracl = is*nxs+ix+1;
- hdr_out[ix].tracf = ix+1;
- hdr_out[ix].scalco = -1000;
- hdr_out[ix].scalel = -1000;
- hdr_out[ix].sdepth = hdr_bin[0].sdepth;
- hdr_out[ix].trid = 1;
- hdr_out[ix].ns = nt;
- hdr_out[ix].trwf = nxs;
- hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
- hdr_out[ix].f1 = fz;
- hdr_out[ix].f2 = fx;
- hdr_out[ix].dt = dt_time*(1E6);
- hdr_out[ix].d1 = dzrcv;
- hdr_out[ix].d2 = dxrcv;
- hdr_out[ix].sx = (int)roundf(fx + (ix*hdr_out[ix].d2));
- hdr_out[ix].gx = (int)roundf(fx + (ix*hdr_out[ix].d2));
- hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ for (it = 0; it < nt_out; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ hdr_out[iy*nx+ix].fldr = it+1;
+ hdr_out[iy*nx+ix].tracl = it*ny*nx+iy*nx+ix+1;
+ hdr_out[iy*nx+ix].tracf = iy*nx+ix+1;
+ hdr_out[iy*nx+ix].scalco = -1000;
+ hdr_out[iy*nx+ix].scalel = -1000;
+ hdr_out[iy*nx+ix].sdepth = sz[iy*nx+ix];
+ hdr_out[iy*nx+ix].selev = -sz[iy*nx+ix];
+ hdr_out[iy*nx+ix].trid = 2;
+ hdr_out[iy*nx+ix].ns = nz;
+ hdr_out[iy*nx+ix].trwf = nx*ny;
+ hdr_out[iy*nx+ix].ntr = hdr_out[iy*nx+ix].fldr*hdr_out[iy*nx+ix].trwf;
+ hdr_out[iy*nx+ix].f1 = z0;
+ hdr_out[iy*nx+ix].f2 = x0;
+ hdr_out[iy*nx+ix].dt = dt*1E6;
+ hdr_out[iy*nx+ix].d1 = dz;
+ hdr_out[iy*nx+ix].d2 = dx;
+ hdr_out[iy*nx+ix].sx = sx[iy*nx+ix];
+ hdr_out[iy*nx+ix].gx = (int)roundf(x0 + (ix*dx))*1000;
+ hdr_out[iy*nx+ix].sy = sy[iy*nx+ix];
+ hdr_out[iy*nx+ix].gy = (int)roundf(y0 + (iy*dy))*1000;
+ hdr_out[iy*nx+ix].offset = (hdr_out[iy*nx+ix].gx - hdr_out[iy*nx+ix].sx)/1000.0;
+ }
}
ret = writeData3D(fp_out, &outdata[is*nxs*nt], hdr_out, nt, nxs);
if (ret < 0 ) verr("error on writing output file.");