diff --git a/marchenko_applications/HomG.c b/marchenko_applications/HomG.c
index 3d855cb25ea373d857ee80e36d93b3d5d6a30263..570bf7ff84a6bfca26244d72a77376285242ca93 100755
--- a/marchenko_applications/HomG.c
+++ b/marchenko_applications/HomG.c
@@ -33,6 +33,7 @@ void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsa
void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout);
void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift);
void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt, int shift);
+void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax);
void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt);
void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt);
void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout);
@@ -60,6 +61,9 @@ char *sdoc[] = {
" xrcv= .................... x-coordinate of first receiver location",
" dzrcv= ................... z-spacing of receivers",
" dxrcv= ................... x-spacing of receivers",
+" shift=0.0 ................ shift per shot",
+" scheme=0 ................. Scheme used for retrieval. 0=Marchenko,",
+" 1=Marchenko with multiple sources, 2=classical",
NULL};
int main (int argc, char **argv)
@@ -68,9 +72,9 @@ int main (int argc, char **argv)
char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
float *indata, *Ghom, *shotdata, *shotdata_jkz, rho, fmin, fmax;
float dt, dx, t0, x0, xmin, xmax1, sclsxgx, f1, f2, dxrcv, dzrcv;
- float *conv, *tmp1, *tmp2, cp;
- int nshots, nt, nx, ntraces, ret, ix, it, is, ir, file_det, nxs, nzs, verbose;
- int pos1, npos, zmax, inx, numb, dnumb, count, scheme, ntmax, ntshift;
+ float *conv, *conv2, *tmp1, *tmp2, cp, shift;
+ int nshots, nt, nx, ntraces, ret, ix, it, is, ir, ig, file_det, nxs, nzs, verbose;
+ int pos1, npos, zmax, inx, numb, dnumb, count, scheme, ntmax, ntshift, shift_num;
segy *hdr_in, *hdr_out, *hdr_shot;
initargs(argc, argv);
@@ -89,6 +93,7 @@ int main (int argc, char **argv)
if (!getparfloat("cp", &cp)) cp = 1500.0;
if (!getparfloat("fmin", &fmin)) fmin=0.0;
if (!getparfloat("fmax", &fmax)) fmax=100.0;
+ if (!getparfloat("shift", &shift)) shift=0.0;
if (!getparint("numb", &numb)) numb=0;
if (!getparint("dnumb", &dnumb)) dnumb=1;
if (!getparint("scheme", &scheme)) scheme = 0;
@@ -168,7 +173,8 @@ int main (int argc, char **argv)
}
vmess("nt: %d nx: %d nshots: %d",nt,nx,nshots);
fclose(fp_shot);
- readSnapData(fshot, &shotdata[0], &hdr_shot[0], 1, nx, nt, 0, nx, 0, nt);
+ readSnapData(fshot, &shotdata[0], &hdr_shot[0], nshots, nx, nt, 0, nx, 0, nt);
+
hdr_out = (segy *)calloc(nxs,sizeof(segy));
Ghom = (float *)malloc(nt*npos*sizeof(float));
@@ -182,6 +188,7 @@ int main (int argc, char **argv)
}
}
conjugate(shotdata_jkz, nt, nx, dt);
+ conjugate(shotdata, nt, nx, dt);
depthDiff(shotdata_jkz, nt, nx, dt, dx, fmin, fmax, cp, 1);
if (verbose) vmess("Applied jkz to source data");
}
@@ -190,14 +197,18 @@ int main (int argc, char **argv)
}
else if (scheme==1) {
vmess("Marchenko representation with multiple sources");
- }
+ }
+ else if (scheme==3) {
+ vmess("Marchenko representation with multiple shot gathers");
+ }
#pragma omp parallel default(shared) \
- private(ix,it,is,indata, hdr_in,fins,fin2,fp_in,conv,tmp1,tmp2)
+ private(ix,it,is,indata, hdr_in,fins,fin2,fp_in,conv,ig,conv2,tmp1,tmp2)
{
indata = (float *)malloc(nt*nx*nxs*sizeof(float));
hdr_in = (segy *)calloc(nx*nxs,sizeof(segy));
conv = (float *)calloc(nx*nt,sizeof(float));
+ conv2 = (float *)calloc(nx*nt,sizeof(float));
if (scheme==2) {
tmp1 = (float *)calloc(nx*nt,sizeof(float));
tmp2 = (float *)calloc(nx*nt,sizeof(float));
@@ -218,9 +229,6 @@ int main (int argc, char **argv)
convol(shotdata, &indata[is*nx*nt], conv, nx, nt, dt, -2);
timeDiff(conv, nt, nx, dt, fmin, fmax, -2);
for (ix=0; ix<nx; ix++) {
- //it=0;
- //Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] = -conv[ix*nt+it];
- //Ghom[it*nxs*nzs+is*nzs+ir] = -conv[ix*nt+(it+nt/2)];
for (it=0; it<nt/2; it++) {
Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv[ix*nt+it]/rho;
Ghom[it*nxs*nzs+is*nzs+ir] += conv[ix*nt+(it+nt/2)]/rho;
@@ -232,8 +240,6 @@ int main (int argc, char **argv)
convol(shotdata, &indata[is*nx*nt], conv, nx, nt, dt, 0);
timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
for (ix=0; ix<nx; ix++) {
- //it=0;
- //Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] -= 2*conv[ix*nt+it];
for (it=0; it<nt/2; it++) {
Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += 2*conv[ix*nt+it]/rho;
Ghom[it*nxs*nzs+is*nzs+ir] += 2*conv[ix*nt+(it+nt/2)]/rho;
@@ -243,28 +249,53 @@ int main (int argc, char **argv)
else if (scheme==2) { //classical representation
convol(&indata[is*nx*nt], shotdata_jkz, tmp1, nx, nt, dt, 0);
depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- corr(&indata[is*nx*nt], shotdata, tmp2, nx, nt, dt, 0);
+ convol(&indata[is*nx*nt], shotdata, tmp2, nx, nt, dt, 0);
+ //corr(&indata[is*nx*nt], shotdata, tmp2, nx, nt, dt, 0);
for (ix = 0; ix < nx; ix++) {
for (it = 0; it < nt; it++) {
- conv[ix*nt+it] = (tmp2[ix*nt+it]+tmp1[ix*nt+it]);
+ conv[ix*nt+it] = tmp2[ix*nt+it]+tmp1[ix*nt+it];
}
}
timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
for (ix=0; ix<nx; ix++) {
- //it=0;
- //Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv[ix*nt+it];
for (it=0; it<nt/2; it++) {
Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv[ix*nt+it]/rho;
Ghom[it*nxs*nzs+is*nzs+ir] += conv[ix*nt+(it+nt/2)]/rho;
}
}
+ }
+ if (scheme==3) { //Marchenko representation with multiple shot gathers
+ depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ for (ig=0; ig<nshots; ig++) {
+ convol(&shotdata[ig*nx*nt], &indata[is*nx*nt], conv, nx, nt, dt, -2);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -2);
+ shift_num = ig*((int)(shift/dt));
+ for (ix = 0; ix < nx; ix++) {
+ for (it = nt/2+1; it < nt; it++) {
+ conv[ix*nt+it] = 0.0;
+ }
+ for (it = shift_num; it < nt; it++) {
+ conv2[ix*nt+it] = conv[ix*nt+it-shift_num];
+ }
+ for (it = 0; it < shift_num; it++) {
+ conv2[ix*nt+it] = conv[ix*nt+nt-shift_num+it];
+ }
+ }
+ for (ix=0; ix<nx; ix++) {
+ Ghom[(-1+nt/2)*nxs*nzs+is*nzs+ir] += conv2[ix*nt+nt-1]/rho;
+ for (it=0; it<nt/2; it++) {
+ Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv2[ix*nt+it]/rho;
+ //Ghom[it*nxs*nzs+is*nzs+ir] += conv2[ix*nt+(it+nt/2)]/rho;
+ }
+ }
+ }
}
}
count+=1;
if (verbose) vmess("Creating Homogeneous Green's function at depth %d from %d depths",count,nzs);
}
- free(conv); free(indata); free(hdr_in);
+ free(conv); free(indata); free(hdr_in); free(conv2);
if (scheme==2) {
free(tmp1);free(tmp2);
}
@@ -763,3 +794,65 @@ void conjugate(float *data, int nsam, int nrec, float dt)
return;
}
+
+void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax)
+{
+ int optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, tom, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (int)MIN((omin/deltom), (nfreq));
+ iomax = MIN((int)(omax/deltom), (nfreq));
+
+ cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ tom = om*shift;
+ cdatascl[ix*nfreq+iom].r = cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom);
+ cdatascl[ix*nfreq+iom].i = cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom);
+ }
+ }
+ free(cdata);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
diff --git a/marchenko_applications/Makefile b/marchenko_applications/Makefile
index ff32af7089be1e8599948d503c18cf21bf059e93..ce724ae80214c18591a1e329cac2ea1c39803aa5 100644
--- a/marchenko_applications/Makefile
+++ b/marchenko_applications/Makefile
@@ -7,7 +7,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#ALL: fmute marchenko marchenko2
-ALL: fmute marchenko_app combine gmshift MuteSnap HomG iba reshape_su
+ALL: fmute marchenko_app combine gmshift MuteSnap HomG iba reshape_su combine_induced
SRCJ = fmute.c \
getFileInfo.c \
@@ -66,6 +66,17 @@ SRCC = combine.c \
docpkge.c \
readSnapData.c
+
+SRCCI = combine_induced.c \
+ getFileInfo.c \
+ writeData.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
+
SRCRS = reshape_su.c \
getFileInfo.c \
writeData.c \
@@ -132,6 +143,11 @@ OBJC = $(SRCC:%.c=%.o)
combine: $(OBJC)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine $(OBJC) $(LIBS)
+OBJCI = $(SRCCI:%.c=%.o)
+
+combine_induced: $(OBJCI)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine_induced $(OBJCI) $(LIBS)
+
OBJRS = $(SRCRS:%.c=%.o)
reshape_su: $(OBJRS)
@@ -157,7 +173,7 @@ OBJIBA = $(SRCIBA:%.c=%.o)
iba: $(OBJIBA)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o iba $(OBJIBA) $(LIBS)
-install: fmute marchenko_app combine gmshift MuteSnap HomG iba reshape_su
+install: fmute marchenko_app combine gmshift MuteSnap HomG iba reshape_su combine_induced
cp fmute $B
cp marchenko_app $B
cp combine $B
@@ -166,14 +182,15 @@ install: fmute marchenko_app combine gmshift MuteSnap HomG iba reshape_su
cp HomG $B
cp iba $B
cp reshape_su $B
+ cp combine_induced $B
# cp marchenko2 $B
clean:
- rm -f core fmute $(OBJJ) marchenko_app $(OBJH) combine $(OBJC) gmshift $(OBJG) MuteSnap $(OBJMS) HomG $(OBJHG) iba $(OBJIBA) reshape_su $(OBJRS)
+ rm -f core fmute $(OBJJ) marchenko_app $(OBJH) combine $(OBJC) gmshift $(OBJG) MuteSnap $(OBJMS) HomG $(OBJHG) iba $(OBJIBA) reshape_su $(OBJRS) combine_induced $(OBJCI)
realclean: clean
- rm -f $B/fmute $B/marchenko_app $B/combine $B/gmshift $B/MuteSnap $B/HomG $B/iba $B/reshape_su
+ rm -f $B/fmute $B/marchenko_app $B/combine $B/gmshift $B/MuteSnap $B/HomG $B/iba $B/reshape_su $B/combine_induced
diff --git a/marchenko_applications/MuteSnap.c b/marchenko_applications/MuteSnap.c
index ee2bba5517623b74f1103f46acf15012b0c8e5c9..b90681ec6965d7a2ff110c80c1f771ee1aee29fb 100755
--- a/marchenko_applications/MuteSnap.c
+++ b/marchenko_applications/MuteSnap.c
@@ -132,7 +132,12 @@ int main (int argc, char **argv)
vmess("First arrival determined through tolerance (=%.4f)",tol);
}
else if (mode == 2) {
- vmess("First arrival determined through raytimes:%d",nz*nx);
+ vmess("First arrival determined through raytimes");
+ fp_snap = fopen(fray,"r");
+ if (fp_snap == NULL) {
+ verr("Could not open file");
+ }
+ fclose(fp_snap);
hdrs_mute = (segy *) calloc(nz,sizeof(segy));
timeval = (float *)calloc(nz*nx,sizeof(float));
readSnapData(fray, timeval, hdrs_mute, nz, 1, nx, 0, 1, 0, nx);
diff --git a/marchenko_applications/combine.c b/marchenko_applications/combine.c
index 30d410680817c84562edd3809a36648cbbc74fcc..a93a8814fcd21d41cee3106d861efeb018ef3959 100755
--- a/marchenko_applications/combine.c
+++ b/marchenko_applications/combine.c
@@ -52,8 +52,8 @@ int main (int argc, char **argv)
char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
float *indata, *outdata, *rtrace, fz, fx;
float dt, dx, t0, x0, xmin, xmax, sclsxgx, dt2, dx2, t02, x02, xmin2, xmax2, sclsxgx2, dxrcv, dzrcv;
- int nshots, nt, nx, ntraces, nshots2, nt2, nx2, ntraces2, ix, it, is, ir, pos, ifile, file_det, nxs, nzs;
- int xcount, numb, dnumb, ret, nzmax, verbose;
+ int nshots, nt, nx, ntraces, nshots2, nt2, nx2, ntraces2, ix, it, is, iz, pos, ifile, file_det, nxs, nzs;
+ int xcount, numb, dnumb, ret, nzmax, transpose, verbose;
segy *hdr_in, *hdr_out;
initargs(argc, argv);
@@ -65,6 +65,7 @@ int main (int argc, char **argv)
if (!getparint("dnumb", &dnumb)) dnumb=0;
if (!getparint("nzmax", &nzmax)) nzmax=0;
if (!getparint("verbose", &verbose)) verbose=0;
+ if (!getparint("transpose", &transpose)) transpose=0;
if (fin == NULL) verr("Incorrect downgoing input");
if (dnumb < 1) dnumb = 1;
@@ -122,31 +123,43 @@ int main (int argc, char **argv)
if (nshots==0) nshots=1;
nxs = ntraces;
+
+
// ngath zijn het aantal schoten
- hdr_out = (segy *)calloc(nxs,sizeof(segy));
- outdata = (float *)calloc(nxs*nzs,sizeof(float));
- hdr_in = (segy *)calloc(nxs,sizeof(segy));
- indata = (float *)calloc(nxs,sizeof(float));
+ hdr_out = (segy *)calloc(nxs*nt,sizeof(segy));
+ outdata = (float *)calloc(nxs*nzs*nt,sizeof(float));
+ hdr_in = (segy *)calloc(nxs*nt,sizeof(segy));
+ indata = (float *)calloc(nxs*nt,sizeof(float));
readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
nshots = hdr_in[nxs-1].fldr;
- nxs = hdr_in[nxs-1].tracf;
+ if (transpose==0) {
+ nxs = hdr_in[nxs-1].tracf;
+ }
+ else {
+ nxs = hdr_in[nxs-1].ns;
+ }
- for (ir = 0; ir < nzs; ir++) {
- if (verbose) vmess("Depth:%d out of %d",ir+1,nzs);
- sprintf(fins,"%d",ir*dnumb+numb);
- sprintf(fin2,"%s%s%s",fbegin,fins,fend);
- fp_in = fopen(fin2, "r");
+ for (iz = 0; iz < nzs; iz++) {
+ if (verbose) vmess("Depth:%d out of %d",iz+1,nzs);
+ sprintf(fins,"%d",iz*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);
- readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
- if (ir==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
- if (ir==1) dzrcv=hdr_in[0].f1-fz;
- for (is = 0; is < nxs; is++) {
- for (it = 0; it < nshots; it++) {
- outdata[it*nxs*nzs+is*nzs+ir] = indata[it*nxs+is];
+ if (transpose==0) {
+ readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ }
+ else {
+ readSnapData(fin2, &indata[0], &hdr_in[0], nshots, 1, nxs, 0, 1, 0, nxs);
+ }
+ if (iz==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
+ if (iz==1) dzrcv=hdr_in[0].f1-fz;
+ for (ix = 0; ix < nxs; ix++) {
+ for (is = 0; is < nshots; is++) {
+ outdata[is*nxs*nzs+ix*nzs+iz] = indata[is*nxs+ix];
}
}
}
diff --git a/marchenko_applications/homogeneousg.c b/marchenko_applications/homogeneousg.c
index bf5f7e75843ca6c7ca6aa28bb59b1e9940cf0498..bdda95beaea96c6779ff0fd30de13007be54329a 100644
--- a/marchenko_applications/homogeneousg.c
+++ b/marchenko_applications/homogeneousg.c
@@ -31,14 +31,13 @@ void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax,
void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt);
void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout);
-void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose)
+void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int n_source, int verbose)
{
int i, j, l, ret, scheme,count=0;
- int iter, niter, ix, kxwfilt;
+ int iter, niter, ix, is, kxwfilt, shift_num;
float scl, *conv, *greenf2, fmin, fmax, alpha, cp, rho, perc;
- float *greenjkz, *greenf2jkz, *tmp1, *tmp2;
+ float *greenjkz, *greenf2jkz, *tmp1, *tmp2, source_shift;
double t0, t2, tfft;
- FILE *fp, *fp1, *fp2, *fp3;
if (!getparint("scheme", &scheme)) scheme = 0;
if (!getparint("kxwfilt", &kxwfilt)) kxwfilt = 0;
@@ -49,6 +48,7 @@ void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int
if (!getparfloat("cp", &cp)) cp = 1500.0;
if (!getparfloat("rho", &rho)) rho = 1000.0;
if (!getparfloat("perc", &perc)) perc = 0.15;
+ if (!getparfloat("source_shift", &source_shift)) source_shift = 0.1;
tfft = 0.0;
ret = 0;
@@ -95,18 +95,22 @@ void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int
conjugate(green, ntfft, nshots, dt);
depthDiff(greenjkz, ntfft, nshots, dt, dx, fmin, fmax, cp, 1);
}
+ else if (scheme==2) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with multiple sources");
+ }
else {
if (verbose) vmess("Marchenko Homogeneous Green's function retrieval");
}
#pragma omp parallel default(shared) \
- private(i,j,conv,tmp1,tmp2)
+ private(i,j,is,conv,tmp1,tmp2)
{
conv = (float *)calloc(nxs*ntfft,sizeof(float));
if (scheme==1) {
tmp1 = (float *)calloc(nxs*ntfft,sizeof(float));
tmp2 = (float *)calloc(nxs*ntfft,sizeof(float));
}
+ if (scheme==3) tmp1 = (float *)calloc(nxs*ntfft,sizeof(float));
#pragma omp for
for (l = 0; l < Nsyn; l++) {
@@ -114,23 +118,25 @@ void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int
count+=1;
if (verbose > 2) vmess("Creating Homogeneous G at location %d out of %d",count,Nsyn);
+ if (scheme==3) vmess("Looping over %d source positions",n_source);
+
if (scheme==0) { //Marchenko representation
depthDiff(&f2p[l*nxs*nts], ntfft, nshots, dt, dx, fmin, fmax, cp, 1);
- convol(green, &f2p[l*nxs*nts], conv, nxs, nts, dt, -2);
- timeDiff(conv, ntfft, nshots, dt, fmin, fmax, -2);
+ convol(green, &f2p[l*nxs*nts], conv, nxs, nts, dt, 0);
+ timeDiff(conv, ntfft, nshots, dt, fmin, fmax, -3);
if (kxwfilt) {
kxwfilter(conv, ntfft, nshots, dt, dx, fmin, fmax, alpha, cp, perc);
}
for (i=0; i<npossyn; i++) {
- for (j=0; j<nts/2; j++) {
- HomG[(j+nts/2)*Nsyn+synpos[l]] += conv[i*nts+j]/rho;
+ for (j=0; j<nts/2; j++) {
+ HomG[(j+nts/2)*Nsyn+synpos[l]] += conv[i*nts+j]/rho;
HomG[j*Nsyn+synpos[l]] += conv[i*nts+(j+nts/2)]/rho;
- }
- }
+ }
+ }
}
else if (scheme==1) { //classical representation
- convol(&greenf2jkz[l*nxs*nts], green, tmp1, nxs, nts, dt, 0);
- convol(&greenf2[l*nxs*nts], greenjkz, tmp2, nxs, nts, dt, 0);
+ convol(&greenf2jkz[l*nxs*nts], green, tmp1, nxs, nts, dt, 0);
+ convol(&greenf2[l*nxs*nts], greenjkz, tmp2, nxs, nts, dt, 0);
for (i = 0; i < npossyn; i++) {
for (j = 0; j < nts; j++) {
conv[i*nts+j] = tmp1[i*nts+j]+tmp2[i*nts+j];
@@ -138,18 +144,58 @@ void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int
}
timeDiff(conv, ntfft, nshots, dt, fmin, fmax, -1);
for (i=0; i<npossyn; i++) {
- for (j=0; j<nts/2; j++) {
- HomG[(j+nts/2)*Nsyn+synpos[l]] += conv[i*nts+j]/rho;
+ for (j=0; j<nts/2; j++) {
+ HomG[(j+nts/2)*Nsyn+synpos[l]] += conv[i*nts+j]/rho;
HomG[j*Nsyn+synpos[l]] += conv[i*nts+(j+nts/2)]/rho;
- }
- }
+ }
+ }
+ }
+ else if (scheme==2) { //Marchenko representation with multiple sources
+ depthDiff(&f2p[l*nxs*nts], ntfft, nshots, dt, dx, fmin, fmax, cp, 1);
+ convol(green, &f2p[l*nxs*nts], conv, nxs, nts, dt, 0);
+ timeDiff(conv, ntfft, nshots, dt, fmin, fmax, -1);
+ for (i=0; i<npossyn; i++) {
+ for (j=0; j<nts/2; j++) {
+ HomG[(j+nts/2)*Nsyn+synpos[l]] += 2*conv[i*nts+j]/rho;
+ HomG[j*Nsyn+synpos[l]] += 2*conv[i*nts+(j+nts/2)]/rho;
+ }
+ }
}
+ if (scheme==3) { //Marchenko representation with multiple shot gathers
+ depthDiff(&f2p[l*nxs*nts], ntfft, nshots, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<n_source; is++) {
+ convol(&green[is*nxs*nts], &f2p[l*nxs*nts], conv, nxs, nts, dt, 0);
+ timeDiff(conv, ntfft, nshots, dt, fmin, fmax, -3);
+ //shift_num = is*((int)(source_shift/dt));
+ if (kxwfilt) {
+ kxwfilter(conv, ntfft, nshots, dt, dx, fmin, fmax, alpha, cp, perc);
+ }
+ for (i=0; i<npossyn; i++) {
+ /*for (j = nts/2+1; j < nts; j++) {
+ tmp1[i*nts+j] = 0.0;
+ }
+ for (j = shift_num; j < nts; j++) {
+ tmp1[i*nts+j] = conv[i*nts+j-shift_num];;
+ }
+ for (j = shift_num; j < nts; j++) {
+ tmp1[i*nts+j] = conv[i*nts+nts-shift_num+j];;
+ }*/
+ //HomG[(nts/2-1)*Nsyn+synpos[l]] += tmp1[i*nts+nts-1]/rho;
+ for (j=0; j<nts/2; j++) {
+ //HomG[(j+nts/2)*Nsyn+synpos[l]] += tmp1[i*nts+j]/rho;
+ HomG[is*nts*Nsyn+(j+nts/2)*Nsyn+synpos[l]] += conv[i*nts+j]/rho;
+ HomG[is*nts*Nsyn+j*Nsyn+synpos[l]] += conv[i*nts+(j+nts/2)]/rho;
+ }
+ }
+ }
+ }
}
free(conv);
if (scheme==1) {
free(tmp1);
free(tmp2);
}
+ if (scheme==3) free(tmp1);
}
if (scheme==1) {
free(greenf2);
@@ -307,6 +353,13 @@ void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax,
cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].r;
cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].i;
}
+ }
+ else if (opt == -3) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = 2*om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
}
}
free(cdata);
diff --git a/marchenko_applications/imaging.c b/marchenko_applications/imaging.c
index d6fa5c4250eca77d3e0cd31ad5f4a4749b869ea5..16185f45fb89bf96001517c405de0b64b5bda088 100644
--- a/marchenko_applications/imaging.c
+++ b/marchenko_applications/imaging.c
@@ -17,6 +17,7 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+void kxwfilter(float *data, int nt, int nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc);
void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout);
void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout);
@@ -27,9 +28,9 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
void imaging(float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose)
{
FILE *fp, *fp_out;
- int i, j, l, ret, count=0;
- int iter, ix, nfreq, first=0;
- float *iRN, *conv, *Gmin, *wavelet, *wav;
+ int i, j, l, ret, count=0, kxwfilt;
+ int iter, ix, nfreq, first=0, im_shift, im_smooth;
+ float *iRN, *conv, *Gmin, *wavelet, *wav, fmin, fmax, alpha, cp, rho, perc, *costaper;
complex *Fop;
double t0, t2, tfft;
segy *hdrs;
@@ -38,6 +39,41 @@ void imaging(float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, i
ret = 0;
t0 = wallclock_time();
nfreq = ntfft/2+1;
+
+ if (!getparint("kxwfilt", &kxwfilt)) kxwfilt = 0;
+ if (!getparint("im_shift", &im_shift)) im_shift = 0;
+ if (!getparint("im_smooth", &im_smooth)) im_smooth = 0;
+ if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (!getparfloat("fmax", &fmax)) fmax = 100.0;
+ if (!getparfloat("alpha", &alpha)) alpha = 65.0;
+ if (!getparfloat("cp", &cp)) cp = 1500.0;
+ if (!getparfloat("rho", &rho)) rho = 1000.0;
+ if (!getparfloat("perc", &perc)) perc = 0.15;
+
+ if (im_shift<0) im_shift=0;
+ if (im_smooth<0) im_smooth=0;
+
+ costaper = (float *)calloc(nxs,sizeof(float));
+
+ if (im_shift>0 && im_smooth>0) {
+ vmess("Applying shift of %d samples and taper of %d samples",im_shift,im_smooth);
+ }
+
+ for (j = 0; j < im_shift; j++) {
+ costaper[j] = 0.0;
+ }
+ for (j = im_shift; j < im_shift+im_smooth; j++) {
+ costaper[j] = (cos(PI*(j-im_shift+im_smooth)/im_smooth)+1)/2.0;
+ }
+ for (j = im_shift+im_smooth; j < nxs-(im_shift+im_smooth); j++) {
+ costaper[j] = 1.0;
+ }
+ for (j = nxs-(im_shift+im_smooth); j < nxs-im_shift; j++) {
+ costaper[j] = (cos(PI*(j-(nxs-(im_shift+im_smooth)))/im_smooth)+1)/2.0;
+ }
+ for (j = nxs-(im_shift); j < nxs; j++) {
+ costaper[j] = 0.0;
+ }
//Image = (float *)malloc(Nsyn*sizeof(float));
Fop = (complex *)calloc(nxs*nw*Nsyn,sizeof(complex));
@@ -49,8 +85,8 @@ void imaging(float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, i
if (verbose>3) vmess("Modeling wavelet for Image");
freqwave(wavelet, WP.nt, WP.dt, WP.fp, WP.fmin, WP.flef, WP.frig, WP.fmax,
WP.t0, WP.db, WP.shift, WP.cm, WP.cn, WP.w, WP.scale, WP.scfft, WP.inv, WP.eps, verbose);
- }
-
+ }
+
/* use f1+ as operator on R in frequency domain */
mode=1;
synthesis(Refl, Fop, f1plus, iRN, nx, nt, nxs, nts, dt, xsyn, Nsyn,
@@ -104,9 +140,12 @@ void imaging(float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, i
else{
convol(&Gmin[l*nxs*nts], &f1plus[l*nxs*nts], conv, nxs, nts, dt, 0);
}
-
+ if (kxwfilt) {
+ if (verbose>8) vmess("Applying kxwfilter");
+ kxwfilter(conv, ntfft, nshots, dt, dx, fmin, fmax, alpha, cp, perc);
+ }
for (i=0; i<nxs; i++) {
- Image[synpos[l]] += conv[i*nts]/((float)(nxs*ntfft));
+ Image[synpos[l]] += costaper[i]*conv[i*nts]/((float)(nxs*ntfft));
}
}
free(conv);
@@ -114,7 +153,7 @@ void imaging(float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, i
}
if (WP.wav) free(wavelet);
- free(Gmin);
+ free(Gmin);free(costaper);
t2 = wallclock_time();
if (verbose) {
diff --git a/marchenko_applications/makeWindow.c b/marchenko_applications/makeWindow.c
index e4ec0498f2f5fa1763fb0ca10e5786302b1dab23..db562b615ddacbeb94d8ea38626c6074c6cd7612 100644
--- a/marchenko_applications/makeWindow.c
+++ b/marchenko_applications/makeWindow.c
@@ -128,8 +128,8 @@ void makeWindow(WavePar WP, char *file_ray, char *file_amp, float dt, float *xrc
if (WP.wav) { /*Apply the wavelet to create a first arrival*/
if (file_amp != NULL || geosp==1) {
for (ig=0; ig<nfreq; ig++) {
- cmute[ig].r = (cwav[ig].r*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
- cmute[ig].i = (cwav[ig].i*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
+ cmute[ig].r = (cwav[ig].r*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/(amp[j*nx+i]*amp[j*nx+i]*ntfft*sqrtf(timeval[j*nx+i]));
+ cmute[ig].i = (cwav[ig].i*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/(amp[j*nx+i]*amp[j*nx+i]*ntfft*sqrtf(timeval[j*nx+i]));
}
}
else { /*Use the raytime only to determine the mutewindow*/
@@ -139,7 +139,6 @@ void makeWindow(WavePar WP, char *file_ray, char *file_amp, float dt, float *xrc
}
}
cr1fft(cmute,&tinv[j*nx*ntfft+i*ntfft],ntfft,1);
- tinv[j*nx*ntfft+i*ntfft] /= ntfft;
}
}
}
diff --git a/marchenko_applications/marchenko.c b/marchenko_applications/marchenko.c
index 7cd6708985fc40a7f92c982e57e2847e201d051d..1b834fd7fa90cb31ed777bcc29a49e13b06d9abb 100644
--- a/marchenko_applications/marchenko.c
+++ b/marchenko_applications/marchenko.c
@@ -51,7 +51,7 @@ double wallclock_time(void);
void makeWindow(WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
void iterations (complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *first, int niter, int verbose);
void imaging (float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose);
-void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose);
+void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int n_source, int verbose);
void AmpEst(float *amp, WavePar WP, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose);
@@ -214,7 +214,7 @@ int main (int argc, char **argv)
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn, *synpos;
int reci, mode, ixa, ixb, n2out, verbose, ntfft;
int iter, niter, niterh, tracf, *muteW, pad, nt0, *hmuteW, *hxnxsyn;
- int hw, smooth, above, shift, *ixpossyn, npossyn, ix, first=1;
+ int hw, smooth, above, shift, *ixpossyn, npossyn, ix, first=1, shot_transpose;
float fmin, fmax, *tapersh, *tapersy, fxf, dxf, fxs2, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
float *hzsyn, *hxsyn, *hxrcvsyn, *hG_d, xloc, zloc, *HomG;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, *J;
@@ -265,10 +265,6 @@ int main (int argc, char **argv)
verr("file_tinv and file_shot cannot be both input pipe");
if (file_img == NULL && file_homg == NULL)
verr("file_img and file_homg cannot both be empty");
- if (!getparstring("file_green", &file_green)) {
- if (verbose) vwarn("parameter file_green not found, assume pipe");
- file_green = NULL;
- }
if (!getparfloat("fmin", &fmin)) fmin = 0.0;
if (!getparfloat("fmax", &fmax)) fmax = 70.0;
if (!getparint("ixa", &ixa)) ixa = 0;
@@ -289,6 +285,7 @@ int main (int argc, char **argv)
if(!getparint("above", &above)) above = 0;
if(!getparint("shift", &shift)) shift=12;
if(!getparint("nb", &nb)) nb=0;
+ if(!getparint("shot_transpose", &shot_transpose)) shot_transpose=0;
if (!getparfloat("bstart", &bstart)) bstart = 1.0;
if (!getparfloat("bend", &bend)) bend = 1.0;
@@ -352,8 +349,14 @@ int main (int argc, char **argv)
WP.xloc = -123456.0;
WP.zloc = -123456.0;
synpos = (int *)calloc(ngath,sizeof(int));
- shot.nz = 1;
- shot.nx = ngath;
+ if (shot_transpose==1) {
+ shot.nx = 1;
+ shot.nz = ngath;
+ }
+ else {
+ shot.nx = ngath;
+ shot.nz = 1;
+ }
shot.n = shot.nx*shot.nz;
for (l=0; l<shot.nz; l++) {
for (j=0; j<shot.nx; j++) {
@@ -420,7 +423,6 @@ int main (int argc, char **argv)
/*================ Allocating all data arrays ================*/
- green = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
f2p = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
pmin = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
f1plus = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
@@ -570,7 +572,6 @@ int main (int argc, char **argv)
vmess("time sampling = %e ", dt);
if (ampest > 0) vmess("Amplitude correction estimation is switched on");
if (nb > 0) vmess("Scaling estimation in %d step(s) from %.3f to %.3f (db=%.3f)",nb,bstart,bend,db);
- 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);
if (file_pmin != NULL) vmess("Pmin output file = %s ", file_pmin);
@@ -603,25 +604,35 @@ int main (int argc, char **argv)
fclose(fp_src);
ret = getFileInfo(file_src, &nt_src, &nx_src, &n_src, &dt_src, &dx_src, &fzsrc, &fxsrc, &srcmin, &srcmax, &srcscl, &ntr_src);
- hdr_src = (segy *) calloc(nx_src,sizeof(segy));
- src_data = (float *)calloc(nx_src*nt_src,sizeof(float));
+ hdr_src = (segy *) calloc(n_src*nx_src,sizeof(segy));
+ src_data = (float *)calloc(n_src*nx_src*nt_src,sizeof(float));
+
+ vmess("n_src:%d, nx_src:%d, nt_src:%d",n_src,nx_src,nt_src);
+
readSnapData(file_src, &src_data[0], hdr_src, n_src, nx_src, nt_src, 0, nx_src, 0, nt_src);
- //pad_data(src_data, nt_src, nx_src, ntfft, green);
+
+ vmess("n_src:%d",n_src);
+
+ green = (float *)calloc(n_src*nxs*ntfft,sizeof(float));
if (nt_src < ntfft) {
- for (l=0;l<nxs;l++) {
- for (j=0;j<nt_src;j++) {
- green[l*ntfft+j] = src_data[l*nt_src+j];
- }
- for (j=nt_src;j<ntfft;j++) {
- green[l*ntfft+j] = 0.0;
+ for (i=0;i<n_src;i++) {
+ for (l=0;l<nxs;l++) {
+ for (j=0;j<nt_src;j++) {
+ green[i*ntfft*nxs+l*ntfft+j] = src_data[i*ntfft*nx_src+l*nt_src+j];
+ }
+ for (j=nt_src;j<ntfft;j++) {
+ green[i*ntfft*nxs+l*ntfft+j] = 0.0;
+ }
}
}
}
- else if (nt_src >= ntfft) {
- for (l=0;l<nxs;l++) {
- for (j=0;j<ntfft;j++) {
- green[l*ntfft+j] = src_data[l*nt_src+j];
+ else if (nt_src >= ntfft) {
+ for (i=0;i<n_src;i++) {
+ for (l=0;l<nxs;l++) {
+ for (j=0;j<ntfft;j++) {
+ green[i*ntfft*nxs+l*ntfft+j] = src_data[i*ntfft*nx_src+l*nt_src+j];
+ }
}
}
}
@@ -678,6 +689,8 @@ int main (int argc, char **argv)
ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
f2p,hG_d,hmuteW,smooth,shift,above,pad,nt0,&first,niterh,verbose);
+ green = (float *)calloc(nxs*ntfft,sizeof(float));
+
/* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -817,40 +830,43 @@ int main (int argc, char **argv)
hdrs_homg = (segy *) calloc(shot.nx,sizeof(segy));
if (hdrs_homg == NULL) verr("allocation for hdrs_out");
- HomG = (float *)calloc(Nsyn*ntfft,sizeof(float));
+ HomG = (float *)calloc(n_src*Nsyn*ntfft,sizeof(float));
homogeneousg(HomG,green,Refl,nx,nt,nxs,nts,dt,xsyn,Nsyn,xrcv,xsrc,fxs2,fxs,dxs,dxsrc,dx,ixa,ixb,
ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
- f2p,G_d,muteW,smooth,shift,above,pad,nt0,synpos,verbose);
+ f2p,G_d,muteW,smooth,shift,above,pad,nt0,synpos,n_src,verbose);
/*============= write output files ================*/
fp_out = fopen(file_homg, "w+");
if (fp_out==NULL) verr("Homogeneous Green's function file %s could not be found",file_homg);
- for (j = 0; j < ntfft; j++) {
- for (i = 0; i < shot.nx; i++) {
- hdrs_homg[i].fldr = j+1;
- hdrs_homg[i].tracl = j*shot.nx+i+1;
- hdrs_homg[i].tracf = i+1;
- hdrs_homg[i].scalco = -1000;
- hdrs_homg[i].scalel = -1000;
- hdrs_homg[i].sdepth = (int)(zloc*1000.0);
- hdrs_homg[i].trid = 1;
- hdrs_homg[i].ns = shot.nz;
- hdrs_homg[i].trwf = shot.nx;
- hdrs_homg[i].ntr = hdrs_homg[i].fldr*hdrs_homg[i].trwf;
- hdrs_homg[i].f1 = zsyn[0];
- hdrs_homg[i].f2 = xsyn[0];
- hdrs_homg[i].dt = dt*(1E6);
- hdrs_homg[i].d1 = (float)zsyn[shot.nx]-zsyn[0];
- hdrs_homg[i].d2 = (float)xsyn[1]-xsyn[0];
- hdrs_homg[i].sx = (int)roundf(xsyn[0] + (i*hdrs_homg[i].d2));
- hdrs_homg[i].gx = (int)roundf(xsyn[0] + (i*hdrs_homg[i].d2));
- hdrs_homg[i].offset = (hdrs_homg[i].gx - hdrs_homg[i].sx)/1000.0;
- }
- ret = writeData(fp_out, &HomG[j*shot.n], hdrs_homg, shot.nz, shot.nx);
- if (ret < 0 ) verr("error on writing output file.");
+ for (l = 0; l < n_src; l++) {
+ for (j = 0; j < ntfft; j++) {
+ for (i = 0; i < shot.nx; i++) {
+ hdrs_homg[i].fldr = j+1;
+ hdrs_homg[i].tracl = j*shot.nx+i+1;
+ hdrs_homg[i].tracf = i+1;
+ hdrs_homg[i].scalco = -1000;
+ hdrs_homg[i].scalel = -1000;
+ hdrs_homg[i].sdepth = (int)(zloc*1000.0);
+ hdrs_homg[i].trid = 1;
+ hdrs_homg[i].ns = shot.nz;
+ hdrs_homg[i].trwf = shot.nx;
+ hdrs_homg[i].ntr = hdrs_homg[i].fldr*hdrs_homg[i].trwf;
+ hdrs_homg[i].f1 = zsyn[0];
+ hdrs_homg[i].f2 = xsyn[0];
+ hdrs_homg[i].dt = dt*(1E6);
+ hdrs_homg[i].d1 = (float)zsyn[shot.nx]-zsyn[0];
+ hdrs_homg[i].d2 = (float)xsyn[1]-xsyn[0];
+ hdrs_homg[i].sx = (int)roundf(xsyn[0] + (i*hdrs_homg[i].d2));
+ hdrs_homg[i].sy = l+1;
+ hdrs_homg[i].gx = (int)roundf(xsyn[0] + (i*hdrs_homg[i].d2));
+ hdrs_homg[i].offset = (hdrs_homg[i].gx - hdrs_homg[i].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &HomG[l*ntfft*shot.n+j*shot.n], hdrs_homg, shot.nz, shot.nx);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
}
fclose(fp_out);
diff --git a/movies/field/Real_line_source.mp4 b/movies/field/Real_line_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..b569725705d7dad61a99abde2994cef4b6e19065
Binary files /dev/null and b/movies/field/Real_line_source.mp4 differ
diff --git a/movies/field/Virtual_line_source.mp4 b/movies/field/Virtual_line_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..9520c1638e4800a12942b272b3281de233dd9b13
Binary files /dev/null and b/movies/field/Virtual_line_source.mp4 differ
diff --git a/movies/field/Virtual_point_source.mp4 b/movies/field/Virtual_point_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..c926b42ff9fa2e7c316efd4f9d07fb6306c17937
Binary files /dev/null and b/movies/field/Virtual_point_source.mp4 differ
diff --git a/movies/line/Modeling.mp4 b/movies/line/Modeling.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..3a5506aa2c4b4db4dfa424bbffa2a08a52b7311d
Binary files /dev/null and b/movies/line/Modeling.mp4 differ
diff --git a/movies/line/Random_real_source.mp4 b/movies/line/Random_real_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..90e513cd954a0fcff029483d104a2c465a9b3c17
Binary files /dev/null and b/movies/line/Random_real_source.mp4 differ
diff --git a/movies/line/Uniform_virtual_source.mp4 b/movies/line/Uniform_virtual_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..415936e9715085edf7f59fad642526919d7adf9c
Binary files /dev/null and b/movies/line/Uniform_virtual_source.mp4 differ
diff --git a/movies/line/random_virtual_source.mp4 b/movies/line/random_virtual_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..719a70d0ef45d0c31112803f4ddfe2b6f58c75b7
Binary files /dev/null and b/movies/line/random_virtual_source.mp4 differ
diff --git a/movies/point/Modeling.mp4 b/movies/point/Modeling.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..135e3a3e1588cb02bf568c66f68cbdc8890cbf0d
Binary files /dev/null and b/movies/point/Modeling.mp4 differ
diff --git a/movies/point/Monopole_source.mp4 b/movies/point/Monopole_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..c47f03bd1d039f0cfef58e4fc825f5fbcf0314b5
Binary files /dev/null and b/movies/point/Monopole_source.mp4 differ
diff --git a/movies/point/Monopole_source_artefacts.mp4 b/movies/point/Monopole_source_artefacts.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..999700cd883877f989c53f08c935a4f772d5d85c
Binary files /dev/null and b/movies/point/Monopole_source_artefacts.mp4 differ
diff --git a/movies/point/double_couple_source.mp4 b/movies/point/double_couple_source.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..08d353369b171abeafa97d06a992c73412e28bc8
Binary files /dev/null and b/movies/point/double_couple_source.mp4 differ