diff --git a/marchenko_applications/HomG_backup21jun2018.c b/marchenko_applications/HomG_backup21jun2018.c
new file mode 100755
index 0000000000000000000000000000000000000000..5ea5136fa41dc7875521c6f443daabc78f6d30f0
--- /dev/null
+++ b/marchenko_applications/HomG_backup21jun2018.c
@@ -0,0 +1,428 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+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);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+int topdet(float *data, int nt);
+
+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);
+void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift);
+
+char *sdoc[] = {
+" ",
+" HomG - Calculate a Homogeneous Green's function ",
+" ",
+" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. First file of the array of receivers",
+" file_shot= ............... File containing the shot location",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first snapshot file",
+" dnumb= ................... integer number of increment in snapshot files",
+" zmax= .................... Integer number of maximum depth level",
+" inx= ..................... Number of sources per depth level",
+" zrcv= .................... z-coordinate of first receiver location",
+" xrcv= .................... x-coordinate of first receiver location",
+" dzrcv= ................... z-spacing of receivers",
+" dxrcv= ................... x-spacing of receivers",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_shot, *fp_out;
+ char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
+ float *indata, *Ghom, *shotdata, *rtrace, *costaper, scl, rho;
+ float dt, dx, t0, x0, xmin, xmax1, sclsxgx, f1, f2, dxrcv, dzrcv, dxpos, offset, dw, *taper;
+ int nshots, nt, nw, nx, ntraces, ret, ix, it, is, ir, pos, ifile, file_det, nxs, nzs, sxmin, sxmax;
+ int pos1, xcount, zcount, npos, zmax, file_cl, ht, inx, numb, dnumb, indrcv, shift;
+ int rmt, smooth, *tol, tolside, tolset, mode, ntap, maxoffset, offset_tmp, count;
+ complex *chom, *cshot, *ctrace;
+ segy *hdr_in, *hdr_out, *hdr_shot;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("fin", &fin)) fin = NULL;
+ if (!getparstring("fshot", &fshot)) fshot = NULL;
+ if (!getparstring("fout", &fout)) fout = "out.su";
+ if (!getparint("zmax", &zmax)) zmax = 0;
+ if (!getparint("inx", &inx)) inx = 0;
+ if (!getparfloat("zrcv", &f1)) f1 = 0;
+ if (!getparfloat("xrcv", &f2)) f2 = 0;
+ if (!getparfloat("dzrcv", &dzrcv)) dzrcv = -1;
+ if (!getparfloat("dxrcv", &dxrcv)) dxrcv = -1;
+ if (!getparfloat("rho", &rho)) rho=1000.0;
+ if (!getparint("numb", &numb)) numb=0;
+ if (!getparint("dnumb", &dnumb)) dnumb=1;
+ if (!getparint("tolset", &tolset)) tolset=10;
+ if (!getparint("mode", &mode)) mode=0;
+ if (!getparint("ntap", &ntap)) ntap=0;
+ if (!getparint("maxoffset", &maxoffset)) maxoffset=0;
+ if (fin == NULL) verr("Incorrect f2 input");
+ if (fshot == NULL) verr("Incorrect Green input");
+
+ if (dnumb == 0) dnumb = 1;
+
+ sprintf(fins,"z%d",numb);
+
+ ptr = strstr(fin,fins);
+ pos1 = ptr - fin + 1;
+
+ sprintf(fbegin,"%*.*s", pos1-1, pos1-1, fin);
+ sprintf(fend,"%s", fin+pos1+1);
+
+ file_det = 1;
+ zcount=0;
+ nzs=0;
+
+ while (file_det) {
+ sprintf(fins,"z%d",nzs*dnumb+numb);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ if (nzs == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzs == 1) {
+ vmess("1 file detected");
+ file_det = 0;
+ break;
+ }
+ else {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzs++;
+ }
+
+ if (inx < 1) {
+ inx = 1;
+ }
+
+ if (zmax < 1) zmax=1;
+ if (zmax < nzs) nzs=zmax;
+
+ nxs = inx;
+ count=0;
+ npos = nxs*nzs;
+
+ vmess("nxs: %d, nzs: %d",nxs,nzs);
+
+ nshots = 0;
+ getFileInfo(fshot, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax1, &sclsxgx, &ntraces);
+
+ if (dxrcv < 0) dxrcv=dx;
+ if (dzrcv < 0) dzrcv=dx;
+
+ // ngath zijn het aantal schoten
+ shotdata = (float *)malloc(nt*nx*nshots*sizeof(float));
+ hdr_shot = (segy *)calloc(nx*nshots,sizeof(segy));
+
+ fp_shot = fopen(fshot,"r");
+ if (fp_shot == NULL) {
+ verr("Could not open file");
+ }
+ vmess("nt: %d nx: %d nshots: %d",nt,nx,nshots);
+ //nx = readData(fp_shot, shotdata, hdr_shot, nt);
+ fclose(fp_shot);
+ readSnapData(fshot, &shotdata[0], &hdr_shot[0], 1, nx, nt, 0, nx, 0, nt);
+
+ hdr_out = (segy *)calloc(nxs,sizeof(segy));
+ Ghom = (float *)malloc(nt*npos*sizeof(float));
+ ht = (int)ceil(nt/2);
+ nw = ht+1;
+ dw = 2.0*(M_PI)/(dt*nt);
+ cshot = (complex *)malloc(nw*nx*sizeof(complex));
+ tol = (int *)malloc(nxs*sizeof(float));
+ taper = (float *)malloc(nx*sizeof(float));
+
+ /*for (ix=0; ix<nx; ix++) {
+ taper[ix] = 1.0;
+ }
+ if (ntap > 0) {//Create taper
+ vmess("Tapering of %d points applied at edges",ntap);
+ for (ix=0; ix<ntap; ix++) {
+ taper[ix] = (cos((M_PI)*(ix-ntap)/ntap)+1)/2.0;
+ taper[nx-1-ix] = (cos((M_PI)*(ix-ntap)/ntap)+1)/2.0;
+ }
+ }*/
+
+ for (ix = 0; ix < nx; ix++) {
+ /*for (it=0; it<nt; it++) {
+ shotdata[ix*nt+it] *= taper[ix];
+ }*/
+ if (hdr_shot[ix].scalco < 0) {
+ offset_tmp = (hdr_shot[ix].sx-hdr_shot[ix].gx)/-hdr_shot[ix].scalco;
+ }
+ else if (hdr_shot[ix].scalco == 0) {
+ offset_tmp = hdr_shot[ix].sx-hdr_shot[ix].gx;
+ }
+ else {
+ offset_tmp = (hdr_shot[ix].sx-hdr_shot[ix].gx)*hdr_shot[ix].scalco;
+ }
+ if (maxoffset > 0 ) {
+ if (abs(offset_tmp) > maxoffset) {
+ for (it=0;it<nt;it++) {
+ shotdata[ix*nt+it] = 0.0;
+ }
+ vmess("Removed offset:%d",offset_tmp);
+ }
+ }
+ rc1fft(&shotdata[ix*nt],&cshot[ix*nw],nt,-1);
+ }
+
+#pragma omp parallel default(shared) \
+ private(offset,ctrace,rtrace,chom,indrcv,rmt,ix,it,is) \
+ private(indata, hdr_in,fins,fin2,fp_in,offset_tmp)
+{
+ chom = (complex *)calloc(nw,sizeof(complex));
+ ctrace = (complex *)malloc(nw*sizeof(complex));
+ rtrace = (float *)malloc(nt*sizeof(float));
+ indata = (float *)malloc(nt*nx*nxs*sizeof(float));
+ hdr_in = (segy *)calloc(nx*nxs,sizeof(segy));
+#pragma omp for
+ for (ir = 0; ir < nzs; ir++) {
+ sprintf(fins,"z%d",ir*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin2, "r");
+ if (fp_in == NULL) {
+ verr("Danger Will Robinson");
+ }
+ fclose(fp_in);
+ readSnapData(fin2, &indata[0], &hdr_in[0], nxs, nx, nt, 0, nx, 0, nt);
+ for (is = 0; is < nxs; is++) {
+ for (ix = 0; ix < nx; ix++) {
+ /*for (it=0; it<nt; it++) {
+ indata[is*nt*nx+ix*nt+it] *= taper[ix];
+ }*/
+ if (hdr_in[is*nx+ix].scalco < 0) {
+ offset_tmp = (hdr_in[is*nx+ix].sx-hdr_in[is*nx+ix].gx)/-hdr_in[is*nx+ix].scalco;
+ }
+ else if (hdr_in[is*nx+ix].scalco == 0) {
+ offset_tmp = hdr_in[is*nx+ix].sx-hdr_in[is*nx+ix].gx;
+ }
+ else {
+ offset_tmp = (hdr_in[is*nx+ix].sx-hdr_in[is*nx+ix].gx)*hdr_in[is*nx+ix].scalco;
+ }
+ if (maxoffset > 0 ) {
+ if (abs(offset_tmp) > maxoffset) {
+ for (it=0;it<nt;it++) {
+ indata[is*nt*nx+ix*nt+it] = 0.0;
+ }
+ //vmess("Removed offset:%d",offset_tmp);
+ }
+ }
+ rc1fft(&indata[is*nt*nx+ix*nt],ctrace,nt,-1);
+ if (mode==0) { //Single source
+ for (it = 1; it < nw; it++) {
+ //chom[it].r += (4/(rho*dw*it))*(ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r);
+ //chom[it].r += (4/(rho*dw*it))*2*(ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);
+ chom[it].r += 2*(ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);
+ //chom[it].r += ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i;
+ //chom[it].r += ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r;
+ }
+ }
+ else { //Multiple sources
+ for (it = 0; it < nw; it++) {
+ /*chom[it].r -= (2/(rho*dw*it))*(ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r);
+ chom[it].i += (2/(rho*dw*it))*(ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);*/
+ chom[it].r -= (ctrace[it].r*cshot[ix*nw+it].i + ctrace[it].i*cshot[ix*nw+it].r);
+ chom[it].i += (ctrace[it].r*cshot[ix*nw+it].r - ctrace[it].i*cshot[ix*nw+it].i);
+ }
+ }
+ }
+ cr1fft(&chom[0],rtrace,nt,1);
+ indrcv = 0;
+ rmt = MIN(nt-indrcv,indrcv)-shift;
+ for (it = 0; it < ht; it++) {
+ if (it > ht-rmt) {
+ Ghom[it*nxs*nzs+is*nzs+ir] = 0.0;
+ }
+ else {
+ Ghom[it*nxs*nzs+is*nzs+ir] = rtrace[ht+it];
+ }
+ }
+ for (it = ht; it < nt; it++) {
+ if (it < ht+rmt) {
+ Ghom[it*nxs*nzs+is*nzs+ir] = 0.0;
+ }
+ else {
+ Ghom[it*nxs*nzs+is*nzs+ir] = rtrace[it-ht];
+ }
+ }
+ memset(&chom[0].r, 0, nw*2*sizeof(float));
+ }
+ //vmess("Creating Homogeneous Green's function at depth %d from %d depths",ir+1,nzs);
+ count+=1;
+ vmess("Creating Homogeneous Green's function at depth %d from %d depths",count,nzs);
+ }
+ free(chom);free(ctrace);free(rtrace);
+ free(indata);free(hdr_in);
+}
+ free(shotdata);
+
+ vmess("nxs: %d nxz: %d f1: %.7f",nxs,nzs,f1);
+
+ fp_out = fopen(fout, "w+");
+
+ for (ir = 0; ir < nt; ir++) {
+ for (ix = 0; ix < nxs; ix++) {
+ hdr_out[ix].fldr = ir+1;
+ hdr_out[ix].tracl = ir*nxs+ix+1;
+ hdr_out[ix].tracf = ix+1;
+ hdr_out[ix].scalco = hdr_shot[0].scalco;
+ hdr_out[ix].scalel = hdr_shot[0].scalel;
+ hdr_out[ix].sdepth = hdr_shot[0].sdepth;
+ hdr_out[ix].trid = 1;
+ hdr_out[ix].ns = nzs;
+ hdr_out[ix].trwf = nxs;
+ hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
+ hdr_out[ix].f1 = f1;
+ hdr_out[ix].f2 = f2/1000;
+ hdr_out[ix].dt = dt*(1E6);
+ hdr_out[ix].d1 = dzrcv;
+ hdr_out[ix].d2 = dxrcv;
+ hdr_out[ix].sx = hdr_shot[0].sx;
+ hdr_out[ix].gx = (int)roundf(f2 + (ix*hdr_out[ix].d2)*1000.0);
+ hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &Ghom[ir*nxs*nzs], hdr_out, nzs, nxs);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ return 0;
+}
+
+void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift)
+{
+ int i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccon, tmp;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccon = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (ccon == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
+ rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
+
+ /* apply convolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccon[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ n = nrec*nfreq;
+ for (j = 0; j < n; j++) {
+ *qr = (*p2r**p1r-*p2i**p1i);
+ *qi = (*p2r**p1i+*p2i**p1r);
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ if (shift) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+ tau = dt*(nsam/2);
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = ccon[j*nfreq+i].r*cos(om*tau) + ccon[j*nfreq+i].i*sin(om*tau);
+ tmp.i = ccon[j*nfreq+i].i*cos(om*tau) - ccon[j*nfreq+i].r*sin(om*tau);
+ ccon[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/((float)(optn));
+ crmfft(&ccon[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,con,nsam);
+
+ free(ccon);
+ free(rdata1);
+ free(rdata2);
+ return;
+}
+
+void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout)
+{
+ int it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsamout+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsamout+it]=0.0;
+ }
+}
+
+void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout)
+{
+ int it,ix;
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsamout ; it++)
+ datout[ix*nsamout+it] = scl*data[ix*nsam+it];
+ }
+}
diff --git a/marchenko_applications/combine_backup4dec2018.c b/marchenko_applications/combine_backup4dec2018.c
new file mode 100755
index 0000000000000000000000000000000000000000..30d410680817c84562edd3809a36648cbbc74fcc
--- /dev/null
+++ b/marchenko_applications/combine_backup4dec2018.c
@@ -0,0 +1,185 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+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);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+
+char *sdoc[] = {
+" ",
+" combine - Combine results into a single result ",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. File containing the first data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first file",
+" dnumb= ................... integer number of increment in files",
+" nzmax= ................... Maximum number of files read",
+NULL};
+
+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, *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;
+ segy *hdr_in, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparint("numb", &numb)) numb=0;
+ if (!getparint("dnumb", &dnumb)) dnumb=0;
+ if (!getparint("nzmax", &nzmax)) nzmax=0;
+ if (!getparint("verbose", &verbose)) verbose=0;
+ if (fin == NULL) verr("Incorrect downgoing input");
+
+ if (dnumb < 1) dnumb = 1;
+
+ sprintf(numb1,"%d",numb);
+
+ ptr = strstr(fin,numb1);
+ pos = ptr - fin + 1;
+
+ sprintf(fbegin,"%*.*s", pos-1, pos-1, fin);
+ sprintf(fend,"%s", fin+pos);
+
+ file_det = 1;
+ nzs=0;
+
+ while (file_det) {
+ sprintf(fins,"%d",nzs*dnumb+numb);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ if (nzs == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzs == 1) {
+ vmess("1 file detected");
+ }
+ else {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzs++;
+ if (nzmax!=0 && nzs == nzmax) {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+
+ sprintf(fins,"%d",numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+
+ sprintf(fins,"%d",numb+dnumb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt2, &nx2, &nshots2, &dt2, &dx2, &t02, &x02, &xmin2, &xmax2, &sclsxgx2, &ntraces2);
+
+ dxrcv=dx*1000;
+ dzrcv=t02-t0;
+
+ 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));
+
+ 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;
+
+ 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");
+ 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];
+ }
+ }
+ }
+ free(indata);
+
+ fp_out = fopen(fout, "w+");
+
+ for (is = 0; is < nshots; 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_in[0].sdepth;
+ hdr_out[ix].trid = 1;
+ hdr_out[ix].ns = nzs;
+ 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*(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;
+ }
+ ret = writeData(fp_out, &outdata[is*nxs*nzs], hdr_out, nzs, nxs);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ return 0;
+}
+
diff --git a/marchenko_applications/combine_induced.c b/marchenko_applications/combine_induced.c
new file mode 100755
index 0000000000000000000000000000000000000000..0b7601e48a0c994bdaad6108e90dca1f6ba71492
--- /dev/null
+++ b/marchenko_applications/combine_induced.c
@@ -0,0 +1,208 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+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);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+
+char *sdoc[] = {
+" ",
+" combine_induced - Combine induced seismicity results together ",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. File containing the first data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first file",
+" dnumb= ................... integer number of increment in files",
+" nzmax= ................... Maximum number of files read",
+NULL};
+
+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, *rtrace, fz, fx, shift, dtshift, dt_time;
+ 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, iz, pos, ifile, file_det, nxs, nzs;
+ int xcount, numb, dnumb, ret, nzmax, verbose, nshot_out, ishift, nshift;
+ segy *hdr_in, *hdr_bin, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparint("numb", &numb)) numb=0;
+ if (!getparint("dnumb", &dnumb)) dnumb=0;
+ if (!getparint("nzmax", &nzmax)) nzmax=0;
+ if (!getparint("verbose", &verbose)) verbose=0;
+ if (!getparint("nshift", &nshift)) nshift=0;
+ if (!getparfloat("shift", &shift)) shift=0.0;
+ if (!getparfloat("dtshift", &dtshift)) dtshift=0.0;
+ if (!getparfloat("dt_time", &dt_time)) dt_time=0.004;
+ if (fin == NULL) verr("Incorrect downgoing input");
+
+ if (dnumb < 1) dnumb = 1;
+
+ sprintf(numb1,"%d",numb);
+
+ ptr = strstr(fin,numb1);
+ pos = ptr - fin + 1;
+
+ sprintf(fbegin,"%*.*s", pos-1, pos-1, fin);
+ sprintf(fend,"%s", fin+pos);
+
+ file_det = 1;
+ nzs=0;
+
+ while (file_det) {
+ sprintf(fins,"%d",nzs*dnumb+numb);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ if (nzs == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzs == 1) {
+ vmess("1 file detected");
+ }
+ else {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzs++;
+ if (nzmax!=0 && nzs == nzmax) {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+
+ sprintf(fins,"%d",numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+
+ sprintf(fins,"%d",numb+dnumb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt2, &nx2, &nshots2, &dt2, &dx2, &t02, &x02, &xmin2, &xmax2, &sclsxgx2, &ntraces2);
+
+ dxrcv=dx*1000;
+ dzrcv=t02-t0;
+
+ if (nshots==0) nshots=1;
+ nxs = ntraces;
+
+ // ngath zijn het aantal schoten
+ hdr_bin = (segy *)calloc(nxs,sizeof(segy));
+ indata = (float *)calloc(nxs*nt,sizeof(float));
+
+ readSnapData(fin2, &indata[0], &hdr_bin[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ nshots = hdr_bin[nxs-1].fldr;
+ nxs = hdr_bin[nxs-1].tracf;
+
+ nshot_out = nshots/2;
+ free(indata);
+
+ hdr_out = (segy *)calloc(nshot_out*nxs,sizeof(segy));
+ outdata = (float *)calloc(nshot_out*nxs*nt,sizeof(float));
+
+ vmess("nshot_out:%d nxs=%d nt:%d",nshot_out,nxs,nt);
+
+#pragma omp parallel default(shared) \
+ private(indata,iz,hdr_in,fins,fin2,fp_in,is,ix,it,ishift)
+{
+ indata = (float *)calloc(ntraces*nt,sizeof(float));
+ hdr_in = (segy *)calloc(ntraces,sizeof(segy));
+
+#pragma omp for
+ 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 (iz==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
+ if (iz==1) dzrcv=hdr_in[0].f1-fz;
+ //ishift = (int)((shift+(dtshift*((float)iz)))/dt_time);
+ ishift = nshift*iz;
+ if (verbose) vmess("Shifting %d 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];
+ }
+ }
+ }
+ }
+ free(indata);free(hdr_in);
+}
+
+ 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;
+ }
+ ret = writeData(fp_out, &outdata[is*nxs*nt], hdr_out, nt, nxs);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ return 0;
+}
+
diff --git a/marchenko_applications/homogeneousg_backup26nov2018.c b/marchenko_applications/homogeneousg_backup26nov2018.c
new file mode 100644
index 0000000000000000000000000000000000000000..0fbd59d6a4e39e24b863dcc838f33abc27b323f4
--- /dev/null
+++ b/marchenko_applications/homogeneousg_backup26nov2018.c
@@ -0,0 +1,550 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+#include "marchenko.h"
+#include "raytime.h"
+
+int omp_get_max_threads(void);
+int omp_get_num_threads(void);
+void omp_set_num_threads(int num_threads);
+
+void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift, int pad, int nt0);
+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 conjugate(float *data, int nsam, int nrec, float dt);
+
+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);
+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 synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, 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, double *tfft, int first, int verbose);
+
+void kxwfilter(float *data, int nt, int nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc);
+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);
+
+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, is, kxwfilt, shift_num;
+ float scl, *conv, *greenf2, fmin, fmax, alpha, cp, rho, perc;
+ 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;
+ if (!getparint("niter",&niter)) niter=10;
+ 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 (!getparfloat("source_shift", &source_shift)) source_shift = 0.1;
+
+ tfft = 0.0;
+ ret = 0;
+ t0 = wallclock_time();
+
+ if (scheme==1) {
+ if (verbose) vmess("Classical Homogeneous Green's function retrieval");
+ greenf2 = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ greenf2jkz = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ greenjkz = (float *)calloc(nxs*ntfft,sizeof(float));
+
+ if (niter<1) {
+ vmess("Single event");
+ }
+ else {
+ vmess("Multiple events");
+ }
+
+ for (l = 0; l < Nsyn; l++) {
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ /* set green to zero if mute-window exceeds nt/2 */
+ if (muteW[l*nxs+ixpossyn[i]] >= nts/2) {
+ memset(&greenf2[l*nxs*nts+i*nts],0, sizeof(float)*nt);
+ memset(&greenf2jkz[l*nxs*nts+i*nts],0, sizeof(float)*nt);
+ continue;
+ }
+ greenf2[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
+ greenf2jkz[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ greenf2[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
+ greenf2jkz[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
+ }
+ }
+ depthDiff(&greenf2jkz[l*nxs*nts], ntfft, nshots, dt, dx, fmin, fmax, cp, 1);
+ }
+
+ for (i = 0; i < npossyn; i++) {
+ for (j = 0; j < nts; j++) {
+ greenjkz[i*nts+j] = green[i*nts+j];
+ }
+ }
+ conjugate(greenjkz, ntfft, nshots, dt);
+ 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,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++) {
+
+ 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);
+ 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;
+ 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);
+ for (i = 0; i < npossyn; i++) {
+ for (j = 0; j < nts; j++) {
+ conv[i*nts+j] = tmp1[i*nts+j]+tmp2[i*nts+j];
+ }
+ }
+ 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;
+ 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, -2);
+ timeDiff(conv, ntfft, nshots, dt, fmin, fmax, -2);
+ 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;
+ }
+ }
+ }
+ }
+ }
+ free(conv);
+ if (scheme==1) {
+ free(tmp1);
+ free(tmp2);
+ }
+ if (scheme==3) free(tmp1);
+}
+ if (scheme==1) {
+ free(greenf2);
+ free(greenf2jkz);
+ }
+
+ t2 = wallclock_time();
+ if (verbose) {
+ vmess("Total Homogeneous G time = %.3f", t2-t0);
+ }
+
+ return;
+}
+
+void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt, int shift)
+{
+ int i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccov, tmp;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccov = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (ccov == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
+ rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
+
+ /* apply correlation */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccov[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ n = nrec*nfreq;
+ for (j = 0; j < n; j++) {
+ *qr = (*p1r * *p2r + *p1i * *p2i);
+ *qi = (*p1i * *p2r - *p1r * *p2i);
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ /* shift t=0 to middle of time window (nsam/2)*/
+ if (shift) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+ tau = dt*(nsam/2);
+
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = ccov[j*nfreq+i].r*cos(om*tau) + ccov[j*nfreq+i].i*sin(om*tau);
+ tmp.i = ccov[j*nfreq+i].i*cos(om*tau) - ccov[j*nfreq+i].r*sin(om*tau);
+ ccov[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&ccov[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,cov,nsam);
+
+ free(ccov);
+ free(rdata1);
+ free(rdata2);
+ return;
+}
+
+void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt)
+{
+ int optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, 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));
+ iomin = MAX(iomin, 1);
+ 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;
+ }
+ if (opt == 1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ cdatascl[ix*nfreq+iom].r = -om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt == -1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = -om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt == -2) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 4.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].r;
+ cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].i;
+ }
+ }
+ }
+ 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;
+}
+
+void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt)
+{
+ int optn, iom, iomin, iomax, nfreq, ix, ikx, nkx, ikxmax;
+ float omin, omax, deltom, df, dkx, *rdata, kx, scl;
+ float kx2, kz2, kp2, kp;
+ complex *cdata, *cdatascl, kz, kzinv;
+
+ optn = optncr(nsam);
+ nfreq = optncr(nsam)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nrec);
+ dkx = 2.0*PI/(nkx*dx);
+ cdata = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad2d_data(data,nsam,nrec,optn,nkx,rdata);
+
+ /* double forward FFT */
+ xt2wkx(&rdata[0], &cdata[0], optn, nkx, optn, nkx, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (int)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((int)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (iom = 0; iom < iomin; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+ }
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkx2xt(&cdatascl[0], &rdata[0], optn, nkx, nkx, optn, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout)
+{
+ int it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsam+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+ for (ix=nrec;ix<nrecout;ix++) {
+ for (it=0;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+}
+void conjugate(float *data, int nsam, int nrec, float dt)
+{
+ int optn, nfreq, j, ix, it, sign, ntdiff;
+ float *rdata, scl;
+ complex *cdata;
+
+ optn = optncr(nsam);
+ ntdiff = optn-nsam;
+ nfreq = optn/2+1;
+
+ 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);
+
+ /* take complex conjugate */
+ for(ix = 0; ix < nrec; ix++) {
+ for(j = 0; j < nfreq; j++) cdata[ix*nfreq+j].i = -cdata[ix*nfreq+j].i;
+ }
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdata[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsam ; it++)
+ data[ix*nsam+it] = scl*rdata[ix*optn+it+ntdiff];
+ }
+ //scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdata);
+ free(rdata);
+
+ return;
+}
+
diff --git a/marchenko_applications/imaging_backup26nov2018.c b/marchenko_applications/imaging_backup26nov2018.c
new file mode 100644
index 0000000000000000000000000000000000000000..16185f45fb89bf96001517c405de0b64b5bda088
--- /dev/null
+++ b/marchenko_applications/imaging_backup26nov2018.c
@@ -0,0 +1,272 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+#include "marchenko.h"
+#include "raytime.h"
+
+int omp_get_max_threads(void);
+int omp_get_num_threads(void);
+void omp_set_num_threads(int num_threads);
+
+void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift, int pad, int nt0);
+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);
+void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift);
+
+void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, 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, double *tfft, int *first, 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)
+{
+ FILE *fp, *fp_out;
+ 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;
+
+ tfft = 0.0;
+ 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));
+ iRN = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ Gmin = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+
+ if (WP.wav) {
+ wavelet = (float *)calloc(ntfft,sizeof(float));
+ 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,
+ xrcv, xsrc, fxs2, fxs, dxs, dxsrc, dx, ixa, ixb, ntfft, nw, nw_low, nw_high, mode,
+ reci, nshots, ixpossyn, npossyn, &tfft, &first, verbose);
+
+ /* compute upgoing Green's G^-,+ */
+ for (l = 0; l < Nsyn; l++) {
+ for (i = 0; i < npossyn; i++) {
+ j=0;
+ Gmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j] - f1min[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Gmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j] - f1min[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+ free(Fop);free(iRN);
+
+ /* Apply mute with window for Gmin */
+ applyMute(Gmin, muteW, smooth, 4, Nsyn, nxs, nts, ixpossyn, npossyn, shift, pad, nt0);
+
+
+ /*fp_out = fopen("Gm.dat","w+");
+ for (i=0; i<nts; i++) {
+ for (j=0; j<npossyn; j++) {
+ fprintf(fp_out,"%.8f\t",Gmin[j*nts+i]);
+ }
+ fprintf(fp_out,"\n");
+ }
+ fclose(fp_out);*/
+
+#pragma omp parallel default(shared) \
+ private(i,conv,wav)
+{
+ conv = (float *)calloc(nxs*ntfft,sizeof(float));
+ if (WP.wav) wav = (float *)calloc(nxs*ntfft,sizeof(float));
+
+#pragma omp for
+ for (l = 0; l < Nsyn; l++) {
+
+ count+=1;
+
+ if (verbose > 2) vmess("Imaging location %d out of %d",count,Nsyn);
+
+ if (WP.wav) {
+ for (i=0; i<npossyn; i++) {
+ convol(&Gmin[l*nxs*nts+i*nts], wavelet, &wav[i*nts], 1, nts, dt, 0);
+ }
+ convol(wav, &f1plus[l*nxs*nts], conv, nxs, nts, dt, 0);
+ }
+ 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]] += costaper[i]*conv[i*nts]/((float)(nxs*ntfft));
+ }
+ }
+ free(conv);
+ if (WP.wav) free(wav);
+}
+
+ if (WP.wav) free(wavelet);
+ free(Gmin);free(costaper);
+
+ t2 = wallclock_time();
+ if (verbose) {
+ vmess("Total Imaging time = %.3f", t2-t0);
+ }
+
+ return;
+}
+
+void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift)
+{
+ int i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccon, tmp;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccon = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (ccon == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
+ rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
+
+ /* apply convolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccon[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ n = nrec*nfreq;
+ for (j = 0; j < n; j++) {
+ *qr = (*p2r**p1r-*p2i**p1i);
+ *qi = (*p2r**p1i+*p2i**p1r);
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ if (shift==1) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+ tau = dt*(nsam/2);
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = ccon[j*nfreq+i].r*cos(om*tau) + ccon[j*nfreq+i].i*sin(om*tau);
+ tmp.i = ccon[j*nfreq+i].i*cos(om*tau) - ccon[j*nfreq+i].r*sin(om*tau);
+ ccon[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+ if (shift==-2) {
+ for (j = 0; j < nrec; j++) {
+ for (i = 0; i < nfreq; i++) {
+ ccon[j*nfreq+i].r = ccon[j*nfreq+i].i;
+ ccon[j*nfreq+i].i = 0.0;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/((float)(optn));
+ crmfft(&ccon[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,con,nsam);
+
+ free(ccon);
+ free(rdata1);
+ free(rdata2);
+ return;
+}
+
+void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout)
+{
+ int it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsamout+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsamout+it]=0.0;
+ }
+}
+
+void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout)
+{
+ int it,ix;
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsamout ; it++)
+ datout[ix*nsamout+it] = scl*data[ix*nsam+it];
+ }
+}