diff --git a/fdelmodc/Makefile b/fdelmodc/Makefile
index 63107542ecbf6f0cee54aff375244967c7d2a1ec..37f88801d8c6e93aee11b2f6825a91da971293a4 100644
--- a/fdelmodc/Makefile
+++ b/fdelmodc/Makefile
@@ -8,7 +8,7 @@ ALLINC = -I.
#LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC = -g -Wall -fsignaling-nans -O0 -fopenmp
#OPTC += -fopenmp -Waddress
-OPTC += -qopt-report=5
+#OPTC += -qopt-report=5 -g
#OPTC := $(subst -O3 -ffast-math, -O1 -g ,$(OPTC))
#PGI options for compiler feedback
#OPTC += -Mprof=lines
diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index 3aeb972c1404aaa91a8e9b3edf17685edea4e40a..934f16bc5a727526e2c2c1265445567df55a0b85 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -156,6 +156,10 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
/* stable implementation changes amplitude and more work is needed */
//vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
//vz[ix*n1+iz] = 0.25*(vz[ix*n1+iz-2]+vz[ix*n1+iz-1]+vz[ix*n1+iz]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
+ }
+ else if (src.type == 10) {
+ tzz[ix*n1+iz-1] -= src_ampl;
+ tzz[ix*n1+iz+1] += src_ampl;
} /* src.type */
diff --git a/fdelmodc/demo/modelOilGas.scr b/fdelmodc/demo/modelOilGas.scr
index 5fa26904388a3d67e07cc72a6bc7861d1b5bec16..6227c3739cffebeb4815b8ba8315ac5dd01980d2 100755
--- a/fdelmodc/demo/modelOilGas.scr
+++ b/fdelmodc/demo/modelOilGas.scr
@@ -16,19 +16,19 @@ makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-export OMP_NUM_THREADS=1
+export OMP_NUM_THREADS=4
zsrc=1100
zsrc=0
which fdelmodc
-fdelmodc \
+../fdelmodc \
file_cp=syncl_cp.su ischeme=1 iorder=4 \
file_den=syncl_ro.su \
file_src=wave.su \
file_rcv=shot_fd.su \
- src_type=7 \
+ src_type=1 \
src_orient=1 \
src_injectionrate=0 \
rec_type_vz=1 \
@@ -37,13 +37,13 @@ fdelmodc \
dtrcv=0.004 \
rec_delay=0.1 \
verbose=2 \
- tmod=2.01 \
+ tmod=0.01 \
dxrcv=10.0 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=$zsrc \
file_snap=snapF_$zsrc \
- tsnap1=0.1 tsnap2=2.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
+ tsnap1=0.1 tsnap2=0.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
ntaper=101 \
snapwithbnd=1 \
left=2 right=2 top=2 bottom=2
diff --git a/fdelmodc/demo/modelOilGas.scr.ok b/fdelmodc/demo/modelOilGas.scr.ok
deleted file mode 100755
index 89907c0d5c8891a821c3b06d27cdbd1f25b94fdf..0000000000000000000000000000000000000000
--- a/fdelmodc/demo/modelOilGas.scr.ok
+++ /dev/null
@@ -1,95 +0,0 @@
-#!/bin/bash
-
-#export PATH=$HOME/bin:$HOME/src/OpenSource/utils:$PATH:
-
-cp=2000
-rho=2500
-dx=2.5
-dt=0.0005
-
-
-makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
- orig=-2500,0 file_base=syncl.su \
- intt=def x=-2500,0,2500 z=250,250,250 poly=0 cp=2300 ro=2000 \
- intt=def x=-2500,-2000,-1000,-800,0,800,2500 z=650,650,700,750,900,750,600 poly=2 cp=2600 ro=2500 \
- intt=def x=-2500,0,2500 z=1390,1390,1390 poly=0 cp=2000 ro=2000
-
-makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-
-export OMP_NUM_THREADS=8
-
-zsrc=1100
-zsrc=0
-
-../fdelmodc \
- file_cp=syncl_cp.su ischeme=1 iorder=4 \
- file_den=syncl_ro.su \
- file_src=wave.su \
- file_rcv=shot_fd.su \
- src_type=7 \
- src_orient=1 \
- src_injectionrate=0 \
- rec_type_vz=1 \
- rec_type_p=1 \
- rec_int_vz=2 \
- dtrcv=0.004 \
- rec_delay=0.1 \
- verbose=2 \
- tmod=2.01 \
- dxrcv=10.0 \
- xrcv1=-2250 xrcv2=2250 \
- zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=$zsrc \
- file_snap=snapF_$zsrc \
- tsnap1=0.1 tsnap2=2.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
- ntaper=101 \
- snapwithbnd=1 \
- left=2 right=2 top=2 bottom=2
-
-suxmovie < snapF_${zsrc}_svz.su loop=1 clip=1e-13
-
-exit
-
-makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
- orig=-2500,0 file_base=hom.su
-
-fdelmodc \
- file_cp=hom_cp.su ischeme=1 iorder=4 \
- file_den=hom_ro.su \
- file_src=wave.su \
- file_rcv=shot_hom_fd.su \
- src_type=7 \
- src_orient=1 \
- src_injectionrate=1 \
- rec_type_vz=1 \
- rec_type_p=1 \
- rec_int_vz=2 \
- dtrcv=0.004 \
- rec_delay=0.1 \
- verbose=2 \
- tmod=4.195 \
- dxrcv=10.0 \
- xrcv1=-2250 xrcv2=2250 \
- zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=0 \
- ntaper=250 \
- left=4 right=4 top=1 bottom=4
-
-cp=2000
-rho=1000
-dx=10
-
-makemod sizex=5000 sizez=2500 dx=$dx dz=2.5 cp0=$cp ro0=$rho \
- orig=-2500,0 file_base=syncl_migr.su \
- intt=def x=-2500,0,2500 z=250,250,250 poly=0 cp=2300 ro=5000 \
- intt=def x=-2500,-2000,-1000,-800,0,800,2500 z=650,650,700,750,900,750,600 poly=2 cp=2600 ro=1000 \
- intt=def x=-2500,0,2500 z=1390,1390,1390 poly=0 cp=2000 ro=5000
-
-sudiff shot_fd_rvz.su shot_hom_fd_rvz.su > diff_rvz.su
-
-makewave fp=20 dt=0.004 file_out=wavedt.su nt=1024 t0=0.0
-
-migr file_shot=diff_rvz.su file_src=wavedt.su file_vel=syncl_migr_cp.su nshots=1 \
- file_image=migr0.su verbose=3 imc=0
-
-
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index f780b78fadff0c9c9a1937dc45663e2df6a2739a..6f6242a26899474661f8ce1ddb8249b0bea4b364 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -127,7 +127,7 @@ char *sdoc[] = {
" verbose=0 ......... silent mode; =1: display info",
" ",
" SHOT AND GENERAL SOURCE DEFINITION:",
-" src_type=1 ........ 1=P 2=Txz 3=Tzz 4=Txx 5=S-pot 6=Fx 7=Fz 8=P-pot 9=double-couple",
+" src_type=1 ........ 1=P 2=Txz 3=Tzz 4=Txx 5=S-pot 6=Fx 7=Fz 8=P-pot 9=double-couple 10=Fz by P",
" src_orient=1 ...... orientation of the source",
" - 1=monopole",
" - 2=dipole +/- vertical oriented",
@@ -749,6 +749,23 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
free(p);
free(q);
}
+ if (bnd.ntap) {
+ free(bnd.tapx);
+ free(bnd.tapz);
+ free(bnd.tapxz);
+ }
+ free(bnd.surface);
+ free(shot.x);
+ free(shot.z);
+ free(src.x);
+ free(src.z);
+ free(src.tbeg);
+ free(src.tend);
+ free(rec.x);
+ free(rec.z);
+ free(rec.xr);
+ free(rec.zr);
+ if(wav.nsamp!=NULL) free(wav.nsamp);
#ifdef MPI
MPI_Finalize();
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index d882ba97c05a56c85c4314bc10ba87480ddcdb25..a316a1f53635357c73247199b83486bd47b12810 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -173,7 +173,6 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
/* 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->tmod += rdelay;
mod->nt = NINT(mod->tmod/mod->dt);
dt = mod->dt;
@@ -925,6 +924,9 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
wav->nst = nsamp; /* put total number of samples in nst part */
}
}
+ if (src->type==7) { /* set also src_injectionrate=1 */
+ src->injectionrate=1;
+ }
if (src->multiwav) {
if (wav->nx != nsrc) {
@@ -956,6 +958,8 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
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,"Fz on P grid with +/-"); break;
}
fprintf(stderr,"\n");
if (wav->random) vmess("Wavelet has a random signature with fmax=%.2f", wav->fmax);
@@ -1110,7 +1114,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
rec->skipdt=NINT(dtrcv/dt);
dtrcv = mod->dt*rec->skipdt;
if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
- if (!getparint("rec_ntsam",&rec->nt)) rec->nt=NINT((mod->tmod-rdelay)/dtrcv)+1;
+ if (!getparint("rec_ntsam",&rec->nt)) rec->nt=(int)round((mod->tmod-rdelay+0.01*mod->dt)/dtrcv)+1;
if (!getparint("rec_int_p",&rec->int_p)) rec->int_p=0;
if (!getparint("rec_int_vx",&rec->int_vx)) rec->int_vx=0;
if (!getparint("rec_int_vz",&rec->int_vz)) rec->int_vz=0;
@@ -1118,7 +1122,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparint("scale",&rec->scale)) rec->scale=0;
if (!getparfloat("dxspread",&dxspread)) dxspread=0;
if (!getparfloat("dzspread",&dzspread)) dzspread=0;
- rec->nt=MIN(rec->nt, NINT((mod->tmod-rdelay)/dtrcv)+1);
+ rec->nt=MIN(rec->nt, NINT((mod->tmod-rdelay+0.01*mod->dt)/dtrcv)+1);
/* allocation of receiver arrays is done in recvPar */
/*
diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr
index 4cff40073965795e47d7403239829ad0349c7754..94b37e1ca0081efa4b5b1581d24b6e444478d768 100755
--- a/marchenko/demo/oneD/primaries.scr
+++ b/marchenko/demo/oneD/primaries.scr
@@ -1,8 +1,9 @@
#!/bin/bash -x
#SBATCH -J marchenko_primaries
-#SBATCH --cpus-per-task=8
+#SBATCH --cpus-per-task=20
#SBATCH --ntasks=1
-#SBATCH --time=0:15:00
+##SBATCH --time=1:29:00
+#SBATCH -p max2h
# Generate the model and Reflection data
#modelpm.scr
@@ -10,7 +11,7 @@
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
-export OMP_NUM_THREADS=8
+export OMP_NUM_THREADS=20
#shot record to remove internal multiples
R=shotsdx5_rp.su
@@ -20,7 +21,7 @@ makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
- niter=31 smooth=10 niterec=2 niterskip=50 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
+ niter=31 shift=20 smooth=10 niterec=2 niterskip=50 file_rr=pred_rr.su T=0 file_update=update.su
exit;
diff --git a/marchenko/demo/oneD/primariesPlane.scr b/marchenko/demo/oneD/primariesPlane.scr
index 5146bb4c3b2266c918cdb027781223d24045cc29..94fe87958e2e20a687b2b4c6c6b2af0e50a74f85 100755
--- a/marchenko/demo/oneD/primariesPlane.scr
+++ b/marchenko/demo/oneD/primariesPlane.scr
@@ -2,7 +2,7 @@
#SBATCH -J marchenko_primaries
#SBATCH --cpus-per-task=10
#SBATCH --ntasks=1
-#SBATCH --time=2:40:00
+#SBATCH --time=0:30:00
#cd $SLURM_SUBMIT_DIR
@@ -10,12 +10,9 @@
export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=10
-#shot record to remove internal multiples
-select=451
-
makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
-../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=10.0 \
- nshots=901 verbose=2 istart=40 iend=1000 fmax=90 pad=1024 \
- niter=20 smooth=10 niterskip=60 niterec=2 shift=20 file_rr=plane2_10_rr.su T=0
+../../marchenko_primaries1 file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=10.0 \
+ nshots=901 verbose=4 istart=40 iend=500 fmax=90 pad=1024 \
+ niter=30 smooth=10 niterskip=60 niterec=2 shift=20 file_rr=plane1_10_rr.su T=0
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index 16c7dc0fc54a7fd05168b931f98324ad2b77e29c..6c572a4111c8bfe0bbe6a0c88afcfd4febd2fcb3 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -18,6 +18,7 @@ void omp_set_num_threads(int num_threads);
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define ISODD(n) ((n) & 01)
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -37,9 +38,7 @@ int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
-void synthesisp(complex *Refl, complex *Fop, float *Top, float *RNi, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
-Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int
-nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
+void synthesisp(complex *Refl, complex *Fop, float *Top, float *RNi, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
void synthesisPositions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose);
@@ -109,7 +108,7 @@ NULL};
int main (int argc, char **argv)
{
FILE *fp_out, *fp_rr, *fp_w, *fp_up;
- size_t nread, sizea, size;
+ size_t nread, size;
int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath, nacq;
int n1, n2, ntap, tap, di, ntraces, tr;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
@@ -119,7 +118,7 @@ int main (int argc, char **argv)
int smooth, *ixpos, *ixp, npos, ix, ixrcv, m, pad, T, isms, isme, perc;
int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift, plane_wave;
float fmin, fmax, tom, deltom, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
- double t0, t1, t2, t3, t4, ttime, t5, tsyn, tread, tfft, tcopy, tii;
+ double t0, t1, t2, t3, t4, t5, ttime, tsyn, tread, tfft, tcopy, tii;
double energyMi, *energyM0;
float tt0, d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *M0, *DD, *RR, *SRC, dt, dx, dxs, scl, mem, scltap;
@@ -171,8 +170,9 @@ int main (int argc, char **argv)
if(!getparint("ishot", &ishot)) ishot=300;
if(!getparint("plane_wave", &plane_wave)) plane_wave=0;
if(!getparfloat("src_angle", &src_angle)) src_angle = 0.0;
- if (!getparfloat("src_velo",&src_velo)) src_velo=1500.;
- if (!getparfloat("t0",&tt0)) tt0=0.1;
+ if(!getparfloat("src_velo",&src_velo)) src_velo=1500.;
+ if(!getparfloat("t0",&tt0)) tt0=0.1;
+ if( (niterskip>1) && ISODD(niter) ) niter++;
if (T>0) {
T=-1;
@@ -421,15 +421,13 @@ int main (int argc, char **argv)
/* use ishot from Refl, complex-conjugate(time reverse), scale with -1 and convolve with wavelet */
if (file_tinv == NULL) {
if (verbose) vmess("Selecting M0 from Refl of %s", file_shot);
- nts = ntfft;
- //nxs = xnx[ishot];
-
- scl = 1.0/((float)2.0*ntfft);
+ nts = ntfft;
+ scl = 1.0/((float)2.0*ntfft);
rtrace = (float *)calloc(ntfft,sizeof(float));
ctrace = (complex *)calloc(nfreq+1,sizeof(complex));
for (i = 0; i < xnx[ishot]; i++) {
- ixrcv = NINT((xrcv[ishot*nx+i]-fxsb)/dxs);
+ ixrcv = NINT((xrcv[ishot*nx+i]-fxsb)/dxs);
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
ctrace[j].r = Refl[ishot*nw*nx+m*nx+i].r*cwave[j].r + Refl[ishot*nw*nx+m*nx+i].i*cwave[j].i;
ctrace[j].i = -Refl[ishot*nw*nx+m*nx+i].i*cwave[j].r + Refl[ishot*nw*nx+m*nx+i].r*cwave[j].i;;
@@ -437,13 +435,8 @@ int main (int argc, char **argv)
/* transfrom result back to time domain */
cr1fft(ctrace, rtrace, ntfft, 1);
for (j = 0; j < nts; j++) {
- DD[0*nxs*nts+ixrcv*nts+j] = -1.0*scl*rtrace[j];
+ DD[0*nxs*nts+ixrcv*nts+j] = scl*rtrace[j];
}
- /* remove taper at edges for selected shot record */
- //if (tap == 2 || tap == 3) scltap= scl/tapersh[i];
- //for (j = 0; j < nts; j++) {
- // DD[0*nxs*nts+i*nts+j] = -1.0*scltap*rtrace[j];
- //}
}
/* set timereversal for searching First Break */
/* for experimenting with non flat truncation windows */
@@ -456,14 +449,14 @@ int main (int argc, char **argv)
/* construct plane wave (time reversed and multiplied with -1) from all shot records */
if (plane_wave==1) {
- nxs = nshots;
- for (l=0; l<nxs; l++) {
+ for (l=0; l<nshots; l++) {
ix = ixpos[l];
memset(ctrace, 0, sizeof(complex)*(nfreq+1));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- tom = j*deltom*twplane[ix];
csum.r=0.0; csum.i=0.0;
for (i = 0; i < xnx[l]; i++) {
+ // ToDo for general acquisitons use ix at receiver position */
+ tom = j*deltom*twplane[i];
csum.r += Refl[l*nw*nx+m*nx+i].r*cos(-tom) - Refl[l*nw*nx+m*nx+i].i*sin(-tom);
csum.i += Refl[l*nw*nx+m*nx+i].i*cos(-tom) + Refl[l*nw*nx+m*nx+i].r*sin(-tom);
}
@@ -474,11 +467,11 @@ int main (int argc, char **argv)
/* transfrom result back to time domain */
cr1fft(ctrace, rtrace, ntfft, 1);
for (j = 0; j < nts; j++) {
- DD[0*nxs*nts+ix*nts+j] = -1.0*scl*rtrace[j];
+ DD[0*nxs*nts+ix*nts+j] = 1.0*scl*rtrace[j];
}
/* compute Source wavelet for plane-wave imaging */
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- tom = j*deltom*twplane[ix];
+ tom = j*deltom*twplane[l];
csum.r = cos(-tom);
csum.i = sin(-tom);
/* Optional add wavelet to SRC-field */
@@ -499,17 +492,14 @@ int main (int argc, char **argv)
xsyn[0] = xsrc[ishot];
zsyn[0] = zsrc[ishot];
xnxsyn[0] = xnx[ishot];
- /* find minimum and maximum postions in header values */
-/*
fxse = fxsb = xrcv[0];
+ /* check consistency of header values */
for (l=0; l<nshots; l++) {
for (i = 0; i < nx; i++) {
fxsb = MIN(xrcv[l*nx+i],fxsb);
fxse = MAX(xrcv[l*nx+i],fxse);
}
}
-*/
-/* not needed xmin and xmax are the same */
dxf = dx;
dxs = dx;
}
@@ -630,22 +620,22 @@ int main (int argc, char **argv)
perc=(iend-istart)/100;if(!perc)perc=1;
if (plane_wave) {
- writeDataIter("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, NINT(src_angle));
+ writeDataIter("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, NINT(src_angle));
}
- /* make DD causal again and undo the -1 multiplication */
+ /* make DD causal again */
for (l=0; l<npos; l++) {
ix = ixpos[l];
j=0;
- SRC[l*nts+j] = -1.0*DD[ix*nts+j];
+ SRC[l*nts+j] = DD[ix*nts+j];
for (j = 1; j < nts; j++) {
- SRC[l*nts+j] = -1.0*DD[ix*nts+nts-j];
+ SRC[l*nts+j] = DD[ix*nts+nts-j];
}
}
if (plane_wave) {
- writeDataIter("DDplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, NINT(src_angle));
- }
- else {
- writeDataIter("DDshot.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, ishot);
+ writeDataIter("DDplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, NINT(src_angle));
+ }
+ else {
+ writeDataIter("DDshot.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, ishot+1);
}
free(SRC);
@@ -666,31 +656,27 @@ int main (int argc, char **argv)
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
ix = ixpos[i];
- iw = NINT((ii*dt+twplane[ix])/dt);
+ iw = NINT((ii*dt+twplane[i])/dt);
for (j = 0; j < nts; j++) {
- M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts+j];
+ M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j];
}
- /* apply mute window for samples above nts-ii */
+ /* apply mute window for samples above nts-ii */
for (j = 0; j < MIN(nts, nts-iw+isms); j++) {
M0[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-iw+isms, k=1; j < MIN(nts, nts-iw+isme); j++, k++) {
- //fprintf(stderr,"j=%d smooth-k=%d taper=%f\n",j,smooth-k,costaper[smooth-k]);
M0[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
}
for (i = 0; i < npos; i++) {
ix = ixpos[i];
j = 0;
- k1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j];
+ k1min[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
- k1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j];
+ k1min[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts+nts-j];
}
}
}
- if (file_iter != NULL) {
- writeDataIter("M0.su", M0, hdrs_out, ntfft, nacq, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii);
- }
}
else { /* use k1min from previous iteration as starting point and do niterec iterations */
niterrun=niterec;
@@ -698,14 +684,13 @@ int main (int argc, char **argv)
if (verbose>1) vmess("Doing %d iterations using previous result at time-sample %d",niterrun,ii);
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
- j=0;
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
- //iw = ii;
- M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - k1min[l*nxs*nts+i*nts+j];
+ M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts] - k1min[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - k1min[l*nxs*nts+i*nts+nts-j];
+ M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts+nts-j] - k1min[l*nxs*nts+i*nts+nts-j];
}
+
/* apply mute window for samples above nts-ii */
for (j = 0; j < MIN(nts,nts-iw+isms); j++) {
M0[l*nxs*nts+i*nts+j] = 0.0;
@@ -716,6 +701,10 @@ int main (int argc, char **argv)
}
}
}
+ if (file_iter != NULL) {
+ writeDataIter("M0.su", M0, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii);
+ }
+
/*================ initialization ================*/
memcpy(Mi, M0, Nfoc*nxs*ntfft*sizeof(float));
@@ -739,15 +728,25 @@ int main (int argc, char **argv)
if (verbose >=2) {
for (l = 0; l < Nfoc; l++) {
- energyMi = 0.0;
- for (i = 0; i < npos; i++) {
- ix = ixpos[i];
- for (j = 0; j < nts; j++) {
- energyMi += RMi[l*nacq*nts+ix*nts+j]*RMi[l*nacq*nts+ix*nts+j];
- }
- }
- if (iter % 2 == 0) {
- if (iter==0) energyM0[l] = energyMi;
+ if (iter % 2 == 0) {
+ energyMi = 0.0;
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ if (recur==0) {
+ for (j = 0; j < nts; j++) {
+ energyMi += RMi[l*nacq*nts+ix*nts+j]*RMi[l*nacq*nts+ix*nts+j];
+ }
+ }
+ else {
+ mem = RMi[l*nacq*nts+ix*nts+nts-1]+DD[l*nacq*nts+ix*nts];
+ energyMi += mem*mem;
+ for (j = 1; j < nts; j++) {
+ mem = RMi[l*nacq*nts+ix*nts+nts-j]+DD[l*nacq*nts+ix*nts+j];
+ energyMi += mem*mem;
+ }
+ }
+ }
+ if ( (iter==0) && (recur==0) ) energyM0[l] = energyMi;
vmess(" - ii %d: Mi at iteration %d has energy %e; relative to M0 %e", ii, iter, sqrt(energyMi), sqrt(energyMi/energyM0[l]));
}
}
@@ -761,9 +760,9 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) {
j = 0;
ix = ixpos[i];
- Mi[l*nxs*nts+i*nts+j] = -RMi[l*nacq*nts+ix*nts+j];
+ Mi[l*nxs*nts+i*nts+j] = RMi[l*nacq*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
- Mi[l*nxs*nts+i*nts+j] = -RMi[l*nacq*nts+ix*nts+nts-j];
+ Mi[l*nxs*nts+i*nts+j] = RMi[l*nacq*nts+ix*nts+nts-j];
}
}
}
@@ -774,7 +773,6 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) {
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
- //iw = ii;
/* apply mute window for samples after ii */
for (j = MAX(0,iw-isme); j < nts; j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0;
@@ -795,9 +793,9 @@ int main (int argc, char **argv)
}
}
}
- if (file_iter != NULL) {
- writeDataIter("v1plus.su", v1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
- }
+ if (file_iter != NULL) {
+ writeDataIter("v1plus.su", v1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
+ }
}
else {/* odd iterations, convolution => k_1^-(t) */
for (l = 0; l < Nfoc; l++) {
@@ -805,19 +803,24 @@ int main (int argc, char **argv)
ix = ixpos[i];
if (recur==1) { /* use k1min from previous iteration */
for (j = 0; j < nts; j++) {
- Mi[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j];
+ Mi[l*nxs*nts+i*nts+j] += -DD[l*nxs*nts+ix*nts+j];
}
j = 0;
k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j];
}
+ //j = 0;
+ //k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
+ //for (j = 1; j < nts; j++) {
+ // k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
+ //}
if (file_update != NULL) {
j=0;
- Mup[l*nxs*nts+i*nts+j] += k1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+ix*nts+j];
- for (j = 1; j < nts; j++) {
- Mup[l*nxs*nts+i*nts+j] += k1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+ix*nts+nts-j];
- }
+ Mup[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j] - k1min[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Mup[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+nts-j] - k1min[l*nxs*nts+i*nts+j];
+ }
}
}
else {
@@ -828,24 +831,14 @@ int main (int argc, char **argv)
}
if (file_update != NULL) {
j=0;
- Mup[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
+ Mup[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Mup[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
+ Mup[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+nts-j];
}
}
}
- /* apply mute window for delta function at t=0
- iw = NINT((twplane[ix])/dt);
- for (j = nts-shift+smooth+iw; j < nts; j++) {
- Mi[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = nts-shift+iw, k=0; j < MIN(nts, nts-shift+smooth+iw); j++, k++) {
- Mi[l*nxs*nts+i*nts+j] *= costaper[k];
- }
- */
/* apply mute window for samples above nts-ii */
iw = NINT((ii*dt+twplane[ix])/dt);
- //iw = ii;
for (j = 0; j < MIN(nts,nts-iw+isms); j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0;
}
@@ -858,6 +851,7 @@ int main (int argc, char **argv)
writeDataIter("k1min.su", k1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
}
} /* end else (iter) branch */
+
if (file_iter != NULL) {
writeDataIter("Mi.su", Mi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
}
@@ -962,5 +956,3 @@ int main (int argc, char **argv)
exit(0);
}
-
-