diff --git a/marchenko3D/applyMute3D.c b/marchenko3D/applyMute3D.c
index b1e0abe60e02d4d6265e4a968df3d296f71e6e4f..7356aa6baa4b554fa227a8639ff92e78db4578bd 100644
--- a/marchenko3D/applyMute3D.c
+++ b/marchenko3D/applyMute3D.c
@@ -200,7 +200,7 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
Nig[j] *= costaper[l];
}
- for (j = MAX(0,tmute-shift+smooth+1); j < MIN(nt,nt+1-tmute+2*ts+shift-smooth); j++) {
+ for (j = MAX(0,tmute-shift+smooth); j < MIN(nt,nt-tmute+2*ts+shift-smooth); j++) {
Nig[j] = 0.0;
}
for (j = MIN(nt-1,nt-tmute+2*ts+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
@@ -218,6 +218,24 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
Nig[j] = 0.0;
}
}
+ else if (above==-2){ //New Theta window keeps Gd
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
+ if (tmute >= nt/2) {
+ memset(&Nig[0],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = MAX(0,-2*ts+tmute+shift+smooth); j < MIN(nt,nt-tmute-shift-smooth); j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ }
else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
@@ -228,7 +246,7 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth-1); j++) {
Nig[j] = 0.0;
}
- for (j = nt+1-tmute+shift+smooth; j < nt; j++) {
+ for (j = nt-tmute+shift+smooth; j < nt; j++) {
Nig[j] = 0.0;
}
for (j = nt-tmute+shift,l=0; j < nt-tmute+shift+smooth; j++,l++) {
diff --git a/marchenko3D/applyMute3D_tshift.c b/marchenko3D/applyMute3D_tshift.c
index 1111c29b164a7958b516dead018a258c6fcf4a51..e58381eae411a49846af9dc08f516792bbfe9a43 100644
--- a/marchenko3D/applyMute3D_tshift.c
+++ b/marchenko3D/applyMute3D_tshift.c
@@ -90,14 +90,32 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
Nig[j] = 0.0;
}
}
+ else if (above==-2){ //New Theta window keeps Gd
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
+ if (tmute >= nt/2) {
+ memset(&Nig[0],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
+ Nig[isyn*nxys*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,-2*ts+tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
+ Nig[isyn*nxys*nt+i*nt+j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ Nig[isyn*nxys*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ }
else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
ts = tsynW[isyn*nxys+imute];
- for (j = -2*ts+tmute-shift-smooth,l=0; j < -2*ts+tmute-shift; j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
Nig[j] *= costaper[smooth-l-1];
}
- for (j = 0; j < -2*ts+tmute-shift-smooth-1; j++) {
+ for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth-1); j++) {
Nig[j] = 0.0;
}
for (j = nt+1-tmute+shift+smooth; j < nt; j++) {
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 3eb355ef66a09b17e5748c5a99b64fde15225bb6..f8fd501224a60f4f22d7b5e9b76a3c19a423ee4d 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -397,7 +397,7 @@ int main (int argc, char **argv)
px = sin(src_anglex*grad2rad)/src_velox;
py = sin(src_angley*grad2rad)/src_veloy;
- tshift = fabs((nys-1)*dys*py) + fabs((nxs-1)*dxs*px);
+ tshift = fabs((nys-1)*dys*py) + fabs((nxs-1)*dxs*px) - dt;
for (j = 0; j < Nfoc; j++) {
itmin[j] = nt;
@@ -602,7 +602,13 @@ int main (int argc, char **argv)
vmess("number of time samples (nt,nts) = %li (%li,%li)", ntfft, nt, nts);
vmess("frequency cutoffs = min:%.3f max:%.3f",fmin,fmax);
vmess("time sampling = %e ", dt);
- if (plane_wave) vmess("Plane wave focusing is applied");
+ if (plane_wave) {
+ vmess("Plane wave focusing is applied");
+ vmess("Plane wave angle = x:%.2f y%.2f",src_anglex,src_angley);
+ vmess("Plane wave velocity = x:%.2f y%.2f",src_velox,src_veloy);
+ vmess("Plane wave p value = x:%.6f y%.6f",px,py);
+ vmess("tshift = x:%.3f",tshift);
+ }
if (file_green != NULL) vmess("Green output file = %s ", file_green);
if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
if (file_gplus != NULL) vmess("Gplus output file = %s ", file_gplus);
@@ -781,7 +787,14 @@ int main (int argc, char **argv)
}
}
}
- if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ if (above==-2) {
+ if (plane_wave==1) {
+ applyMute3D_tshift(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, iter, tsynW);
+ }
+ else {
+ applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ }
+ }
}
else {/* odd iterations update: => f_1^+(t) */
for (l = 0; l < Nfoc; l++) {
diff --git a/utils/pwshift.c b/utils/pwshift.c
index 74ccaa1ba8f5590eeed6059a59862436de08251a..e72dd4a17ae2b59f30e800c9c092235ca8bf88f6 100644
--- a/utils/pwshift.c
+++ b/utils/pwshift.c
@@ -28,7 +28,7 @@ long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps,
long nx, long ny, long nz, long sx, long ex, long sy, long ey, long sz, long ez);
-void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax);
+void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float delay, float fmin, float fmax);
void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
@@ -64,9 +64,9 @@ int main (int argc, char **argv)
FILE *fp_gmin, *fp_ray, *fp_amp, *fp_out;
char *file_gmin, *file_ray, *file_amp, *file_out, fbr[100], fer[100], fba[100], fea[100], fins[100], fin2[100], numb1[100], *ptr;
float *gmin, *conv, *time, *amp, *image, fmin, fmax;
- float dt, dy, dx, dz, t0, y0, x0, scl, *shift, *block;
+ float dt, dy, dx, dz, t0, y0, x0, te, ye, xe, scl, *shift, *block, sy, sx;
float dt_ray, dy_ray, dx_ray, t0_ray, y0_ray, x0_ray, scl_ray, px, py, src_velox, src_veloy, src_anglex, src_angley, grad2rad;
- long nshots, nt, ny, nx, ntr, delay;
+ long nshots, nt, ny, nx, ntr, *delay, nx1, ny1, tmp;
long nray, nt_ray, ny_ray, nx_ray, ntr_ray;
long verbose, ix, iy, it, iz, is, *gx, *gy, *gz, numb, dnumb, pos, nzs, file_det;
size_t ret;
@@ -146,11 +146,16 @@ int main (int argc, char **argv)
*----------------------------------------------------------------------------*/
getFileInfo3D(file_gmin, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &scl, &ntr);
+ xe = x0 + ((float)(nx-1))*dx;
+ ye = y0 + ((float)(ny-1))*dy;
+ te = t0 + ((float)(nt-1))*dt;
+
if (verbose) {
vmess("************************ Gmin info ************************");
vmess("Number of depth levels : %li",nshots);
vmess("Number of samples x: %li, y: %li, t: %li",nx,ny,nt);
vmess("Starting distance for x: %.3f, y: %.3f, t: %.3f",x0,y0,t0);
+ vmess("Ending distance for x: %.3f, y: %.3f, t: %.3f",xe,ye,te);
vmess("Sampling distance for x: %.3f, y: %.3f, t: %.3f",dx,dy,dt);
vmess("***********************************************************");
}
@@ -186,6 +191,7 @@ int main (int argc, char **argv)
gz = (long *)calloc(nray*nshots,sizeof(long));
block = (float *)calloc(nray*nshots*nt,sizeof(float));
+ delay = (long *)calloc(nray,sizeof(long));
readSnapData3D(file_gmin, gmin, hdr_gmin, nshots, nx, ny, nt, 0, nx, 0, ny, 0, nt);
if (verbose) vmess("Read in Gmin data");
@@ -193,41 +199,67 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Add the delay in case the plane wave is at an angle
*----------------------------------------------------------------------------*/
+
+ hdr_time = (segy *)calloc(ny*nray,sizeof(segy));
+ time = (float *)calloc(nx*ny*nray,sizeof(float));
+
+ sprintf(fins,"z%li",numb);
+ sprintf(fin2,"%s%s%s",fbr,fins,fer);
+ readSnapData3D(fin2, time, hdr_time, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+
+ ny1 = 1;
+ tmp = hdr_time[0].sy;
+ for (it=1; it<nray; it++) {
+ if (tmp!=hdr_time[it*ny].sy)
+ {
+ ny1++;
+ tmp = hdr_time[it*ny].sy;
+ }
+ }
+ nx1=nray/ny1;
+ vmess("ny1=%li nx1=%li nray=%li",ny1,nx1,nray);
- shift = (float *)calloc(nx*ny,sizeof(float));
+ shift = (float *)calloc(nray,sizeof(float));
grad2rad = 17.453292e-3;
- px = sin(src_anglex*grad2rad)/src_velox;
- py = sin(src_angley*grad2rad)/src_veloy;
- if (verbose) vmess("px value is %f and py value is %f",px,py);
-
- // if (py < 0.0) {
- // for (iy=0; iy<ny; iy++) {
- // if (px < 0.0) {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + fabsf((ny-1-iy)*dy*py);
- // }
- // }
- // else {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = ix*dx*px + fabsf((ny-1-iy)*dy*py);
- // }
- // }
- // }
- // }
- // else {
- // for (iy=0; iy<ny; iy++) {
- // if (px < 0.0) {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + iy*dy*py;
- // }
- // }
- // else {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = ix*dx*px + iy*dy*py;
- // }
- // }
- // }
- // }
+ if (src_anglex==0.0) px=0.0;
+ else px = sin(-1.0*src_anglex*grad2rad)/src_velox;
+ if (src_angley==0.0) py=0.0;
+ else py = sin(-1.0*src_angley*grad2rad)/src_veloy;
+ if (verbose) {
+ vmess("************************ plane wave info *************************");
+ vmess("Plane wave angle x: %.3f, y: %.3f",src_anglex,src_angley);
+ vmess("Plane wave velocity x: %.3f, y: %.3f",src_velox,src_veloy);
+ vmess("Plane wave slowness x: %.3f, y: %.3f",px,py);
+ vmess("***********************************************************");
+ }
+
+ for (iy=0; iy<ny1; iy++) {
+ for (ix=0; ix<nx1; ix++) {
+ sy = ((float)(hdr_time[(iy*nx1+ix)*ny].sy))/1000.0;
+ sx = ((float)(hdr_time[(iy*nx1+ix)*ny].sx))/1000.0;
+ //vmess("sy=%.3f sx=%.3f",sy,sx);
+ if (py < 0.0) {
+ if (px < 0.0) {
+ shift[iy*nx1+ix] = fabsf((xe-sx)*px) + fabsf((ye-sy)*py);
+ }
+ else {
+ shift[iy*nx1+ix] = fabsf((x0-sx)*px) + fabsf((ye-sy)*py);
+ }
+ }
+ else {
+ if (px < 0.0) {
+ shift[iy*nx1+ix] = fabsf((xe-sx)*px) + fabsf((y0-sy)*py);
+ }
+ else {
+ shift[iy*nx1+ix] = fabsf((x0-sx)*px) + fabsf((y0-sy)*py);
+ }
+ }
+ }
+ }
+
+
+ free(hdr_time); free(time);
+
/*----------------------------------------------------------------------------*
* Apply the imaging condition
@@ -247,7 +279,7 @@ int main (int argc, char **argv)
sprintf(fin2,"%s%s%s",fbr,fins,fer);
readSnapData3D(fin2, time, hdr_time, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
sprintf(fin2,"%s%s%s",fba,fins,fea);
- readSnapData3D(file_amp, amp, hdr_amp, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+ readSnapData3D(file_amp, amp, hdr_amp, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
for (is = 0; is < nray; is++) {
gx[is*nshots+iy] = hdr_time[is*ny].sx;
@@ -257,34 +289,35 @@ int main (int argc, char **argv)
x0_ray = ((float)gx[is*nshots+iy])/1000.0;
y0_ray = ((float)gy[is*nshots+iy])/1000.0;
- if (py < 0.0) {
- if (px < 0.0) {
- delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
- }
- else {
- delay = NINT((x0_ray-x0)*px/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
- }
- }
- else {
- if (px < 0.0) {
- delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT((y0_ray-y0)*py/dt);
- }
- else {
- delay = NINT((x0_ray-x0)*px/dt) + NINT((y0_ray-y0)*py/dt);
- }
- }
-
- if (delay > nt-1) delay = nt-1;
+ // if (py < 0.0) {
+ // if (px < 0.0) {
+ // delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
+ // }
+ // else {
+ // delay = NINT((x0_ray-x0)*px/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
+ // }
+ // }
+ // else {
+ // if (px < 0.0) {
+ // delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT((y0_ray-y0)*py/dt);
+ // }
+ // else {
+ // delay = NINT((x0_ray-x0)*px/dt) + NINT((y0_ray-y0)*py/dt);
+ // }
+ // }
+
+ // if (delay > nt-1) delay = nt-1;
for (it = 0; it < nx*ny*nt; it++) {
conv[it] = gmin[iy*nx*ny*nt+it];
}
- timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],shift,fmin,fmax);
+ timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],shift[is],fmin,fmax);
for (ix = 0; ix < ny*nx; ix++) {
- image[is*nshots+iy] += conv[ix*nt+delay+2]*dx*dy*dt;
+ image[is*nshots+iy] += conv[ix*nt]*dx*dy*dt;
for (it=0; it<nt; it++) {
block[is*nshots*nt+iy*nt+it] += conv[ix*nt+it];
}
+ block[is*nshots*nt+iy*nt] = 1e20;
}
}
@@ -324,37 +357,37 @@ int main (int argc, char **argv)
fclose(fp_out);
- fp_out = fopen("block.su", "w+");
-
- if (nshots>1) dz = ((float)(gz[1] - gz[0]))/1000.0;
- else dz = 1.0;
-
- for (is = 0; is < nshots; is++) {
- for (it = 0; it < nray; it++) {
- hdr_out[it*nshots+is].fldr = it+1;
- hdr_out[it*nshots+is].tracl = it+1;
- hdr_out[it*nshots+is].tracf = it+1;
- hdr_out[it*nshots+is].scalco = -1000;
- hdr_out[it*nshots+is].scalel = -1000;
- hdr_out[it*nshots+is].trid = 1;
- hdr_out[it*nshots+is].ns = nt;
- hdr_out[it*nshots+is].trwf = nray;
- hdr_out[it*nshots+is].ntr = nray;
- hdr_out[it*nshots+is].f1 = (((float)gz[0])/1000.0);
- hdr_out[it*nshots+is].f2 = (((float)gx[0])/1000.0);
- hdr_out[it*nshots+is].dt = ((int)(dt*1E6));
- hdr_out[it*nshots+is].d1 = roundf(dz*1000.0)/1000.0;
- hdr_out[it*nshots+is].d2 = roundf(dx*1000.0)/1000.0;
- hdr_out[it*nshots+is].gx = gx[it*nshots+is];
- hdr_out[it*nshots+is].gy = gy[it*nshots+is];
- hdr_out[it*nshots+is].sdepth = gz[it*nshots+is];
- }
- }
+ // fp_out = fopen("block.su", "w+");
+
+ // if (nshots>1) dz = ((float)(gz[1] - gz[0]))/1000.0;
+ // else dz = 1.0;
+
+ // for (is = 0; is < nshots; is++) {
+ // for (it = 0; it < nray; it++) {
+ // hdr_out[it*nshots+is].fldr = it+1;
+ // hdr_out[it*nshots+is].tracl = it+1;
+ // hdr_out[it*nshots+is].tracf = it+1;
+ // hdr_out[it*nshots+is].scalco = -1000;
+ // hdr_out[it*nshots+is].scalel = -1000;
+ // hdr_out[it*nshots+is].trid = 1;
+ // hdr_out[it*nshots+is].ns = nt;
+ // hdr_out[it*nshots+is].trwf = nray;
+ // hdr_out[it*nshots+is].ntr = nray;
+ // hdr_out[it*nshots+is].f1 = (((float)gz[0])/1000.0);
+ // hdr_out[it*nshots+is].f2 = (((float)gx[0])/1000.0);
+ // hdr_out[it*nshots+is].dt = ((int)(dt*1E6));
+ // hdr_out[it*nshots+is].d1 = dt;
+ // hdr_out[it*nshots+is].d2 = roundf(dx*1000.0)/1000.0;
+ // hdr_out[it*nshots+is].gx = gx[it*nshots+is];
+ // hdr_out[it*nshots+is].gy = gy[it*nshots+is];
+ // hdr_out[it*nshots+is].sdepth = gz[it*nshots+is];
+ // }
+ // }
- ret = writeData3D(fp_out, &block[0], hdr_out, nt, nray*nshots);
- if (ret < 0 ) verr("error on writing output file.");
+ // ret = writeData3D(fp_out, &block[0], hdr_out, nt, nray*nshots);
+ // if (ret < 0 ) verr("error on writing output file.");
- fclose(fp_out);
+ // fclose(fp_out);
free(image); free(hdr_out); free(block);
@@ -363,7 +396,7 @@ int main (int argc, char **argv)
return 0;
}
-void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax)
+void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float delay, float fmin, float fmax)
{
long optn, iom, iomin, iomax, nfreq, ix, sign;
float omin, omax, deltom, om, tom, df, *rdata, scl;
@@ -406,7 +439,7 @@ void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *
}
for (iom = iomin ; iom < iomax ; iom++) {
om = deltom*iom;
- tom = om*-1.0*(time[ix]-delay[ix]);
+ tom = om*-1.0*(time[ix]+delay);
cdatascl[ix*nfreq+iom].r = (cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom))/(amp[ix]*amp[ix]);
cdatascl[ix*nfreq+iom].i = (cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom))/(amp[ix]*amp[ix]);
}
@@ -443,4 +476,4 @@ 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];
}
-}
\ No newline at end of file
+}