3D

marchenko3D/Makefile

@@ -7,7 +7,7 @@ LIBS    += -L$L -lgenfft -lm $(LIBSM)
 
 #ALL: fmute marchenko marchenko2
 
-ALL: fmute marchenko test fmute3D
+ALL: fmute marchenko marchenko3D fmute3D
 
 SRCJ	= fmute.c \
 		getFileInfo.c  \
@@ -69,8 +69,8 @@ fmute:	$(OBJJ)
 
 OBJT	= $(SRCT:%.c=%.o)
 
-test:	$(OBJT) 
-	$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o test $(OBJT) $(LIBS)
+marchenko3D:	$(OBJT) 
+	$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko3D $(OBJT) $(LIBS)
 
 OBJH	= $(SRCH:%.c=%.o)
 
@@ -90,13 +90,13 @@ fmute3D:	$(OBJJ3)
 install: fmute marchenko test fmute3D
 	cp fmute $B
 	cp marchenko $B
-	cp test $B
+	cp marchenko3D $B
 	cp fmute3D $B
 
 #	cp marchenko2 $B
 
 clean:
-		rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2) test $(OBJT) fmute3D $(OBJJ3)
+		rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2) marchenko3D $(OBJT) fmute3D $(OBJJ3)
 
 realclean: clean
-		rm -f $B/fmute $B/marchenko $B/marchenko2 $B/test $B/fmute3D
+		rm -f $B/fmute $B/marchenko $B/marchenko2 $B/marchenko3D $B/fmute3D

marchenko3D/applyMute3D.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <assert.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))
+
+void applyMute3D( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift)
+{
+    int ix, iy, i, j, l, isyn;
+    float *costaper, scl;
+    int imute, tmute;
+
+    if (smooth) {
+        costaper = (float *)malloc(smooth*sizeof(float));
+        scl = M_PI/((float)smooth);
+        for (ix=0; ix<smooth; ix++) {
+            costaper[ix] = 0.5*(1.0+cos((ix+1)*scl));
+        }
+    }
+
+    for (isyn = 0; isyn < Nfoc; isyn++) {
+        if (above==1) {
+            for (i = 0; i < npos; i++) {
+                imute = ixpos[i];
+                tmute = mute[isyn*nxs+imute];
+                for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+                    data[isyn*nxs*nt+i*nt+j] = 0.0;
+                }
+                for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+                }
+            }
+        }
+        else if (above==0){
+            for (i = 0; i < npos; i++) {
+                imute = ixpos[i];
+                tmute = mute[isyn*nxs+imute];
+                if (tmute >= nt/2) {
+                    memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
+                    continue;
+                }
+                for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+                }
+                for (j = MAX(0,tmute-shift+smooth)+1; j < MIN(nt,nt+1-tmute+shift-smooth); j++) {
+                    data[isyn*nxs*nt+i*nt+j] = 0.0;
+                }
+                for (j = MIN(nt,nt-tmute+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+                }
+            }
+        }
+        else if (above==-1){
+            for (i = 0; i < npos; i++) {
+                imute = ixpos[i];
+                tmute = mute[isyn*nxs+imute];
+                for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+                }
+                for (j = MAX(0,tmute-shift+smooth); j < nt; j++) {
+                    data[isyn*nxs*nt+i*nt+j] = 0.0;
+                }
+            }
+        }
+        else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
+            for (i = 0; i < npos; i++) {
+                imute = ixpos[i];
+                tmute = mute[isyn*nxs+imute];
+                for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+                }
+                for (j = 0; j < MAX(0,tmute-shift-smooth-1); j++) {
+                    data[isyn*nxs*nt+i*nt+j] = 0.0;
+                }
+                for (j = MIN(nt,nt+1-tmute+shift+smooth); j < nt; j++) {
+                    data[isyn*nxs*nt+i*nt+j] = 0.0;
+                }
+                for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+                }
+            }
+        }
+        else if (above==2){//Separates the direct part of the wavefield from the coda
+            for (i = 0; i < npos; i++) {
+                imute = ixpos[i];
+                tmute = mute[isyn*nxs+imute];
+                for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+                    data[isyn*nxs*nt+i*nt+j] = 0.0;
+                }
+                for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+                }
+                for (j = MAX(0,tmute+shift),l=0; j < MAX(0,tmute+shift+smooth); j++,l++) {
+                    data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+                }
+                for (j = MAX(0,tmute+shift+smooth); j < nt; j++) {
+                    data[isyn*nxs*nt+i*nt+j] = 0.0;
+                }
+            }
+        }
+    }
+
+    if (smooth) free(costaper);
+
+    return;
+}
+

marchenko3D/fmute3D.c

+#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 getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, float *sclsxgxsygy, int *nxm);
+int readData3D(FILE *fp, float *data, segy *hdrs, int n1);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int disp_fileinfo3D(char *file, int n1, int n2, int n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
+void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift);
+double wallclock_time(void);
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" fmute3D - mute in time domain file_shot along 3D curves of maximum amplitude in file_mute ",
+" ",
+" fmute3D file_shot= {file_mute=} [optional parameters]",
+" ",
+" Required parameters: ",
+" ",
+"   file_mute= ................ input file with event that defines the mute line",
+"   file_shot= ................ input data that is muted",
+" ",
+" Optional parameters: ",
+" ",
+"   file_out= ................ output file",
+"   above=0 .................. mute after(0), before(1) or around(2) the maximum times of file_mute",
+"   .......................... options 4 is the inverse of 0 and -1 the inverse of 1",
+"   shift=0 .................. number of points above(positive) / below(negative) maximum time for mute",
+"   check=0 .................. plots muting window on top of file_mute: output file check.su",
+"   scale=0 .................. scale data by dividing through maximum",
+"   hw=15 .................... number of time samples to look up and down in next trace for maximum",
+"   smooth=0 ................. number of points to smooth mute with cosine window",
+//"   nxmax=512 ................ maximum number of traces in input file",
+//"   ntmax=1024 ............... maximum number of samples/trace in input file",
+"   verbose=0 ................ silent option; >0 display info",
+" ",
+" author  :  Jan Thorbecke     : 2012 (janth@xs4all.nl)",
+" author 3D: Joeri Brackenhoff : 2019 (j.a.brackenhoff@tudelft.nl)"
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+    FILE    *fp_in1, *fp_in2, *fp_out, *fp_chk, *fp_psline1, *fp_psline2;
+    int        verbose, shift, k, nx1, ny1, nt1, nx2, ny2, nt2, nxy;
+    int     ntmax, nxmax, nymax, ret, i, j, l, jmax, ixmax, iymax, above, check;
+    int     size, ntraces, ngath, *maxval, hw, smooth;
+    int     tstart, tend, scale, *xrcv;
+    float   dt, dt1, dx1, dy1, ft1, fx1, fy1, t0, t1, dt2, dx2, dy2, ft2, fx2, fy2;
+    float    w1, w2, dxrcv, dyrcv;
+    float     *tmpdata, *tmpdata2, *costaper;
+    char     *file_mute, *file_shot, *file_out;
+    float   scl, sclsxgx, sclshot, xmin, xmax, tmax, lmax;
+    segy    *hdrs_in1, *hdrs_in2;
+
+    t0 = wallclock_time();
+    initargs(argc, argv);
+    requestdoc(1);
+
+    if(!getparstring("file_mute", &file_mute)) file_mute=NULL;
+    if(!getparstring("file_shot", &file_shot)) file_shot=NULL;
+    if(!getparstring("file_out", &file_out)) file_out=NULL;
+    if(!getparint("ntmax", &ntmax)) ntmax = 1024;
+    if(!getparint("nxmax", &nxmax)) nxmax = 512;
+    if(!getparint("above", &above)) above = 0;
+    if(!getparint("check", &check)) check = 0;
+    if(!getparint("scale", &scale)) scale = 0;
+    if(!getparint("hw", &hw)) hw = 15;
+    if(!getparint("smooth", &smooth)) smooth = 0;
+    if(!getparfloat("w1", &w1)) w1=1.0;
+    if(!getparfloat("w2", &w2)) w2=1.0;
+    if(!getparint("shift", &shift)) shift=0;
+    if(!getparint("verbose", &verbose)) verbose=0;
+
+/* Reading input data for file_mute */
+
+    if (file_mute != NULL) {
+        ngath = 1;
+        ret = getFileInfo3D(file_mute, &nt1, &nx1, &ny1, &ngath, &dt1, &dx1, &dy1, &ft1, &fx1, &fy1, &sclsxgx, &ntraces);
+        
+        if (!getparint("ntmax", &ntmax)) ntmax = nt1;
+        if (!getparint("nxmax", &nxmax)) nxmax = nx1;
+        if (!getparint("nymax", &nymax)) nymax = ny1;
+        if (verbose>=2 && (ntmax!=nt1 || nxmax!=nx1 || nymax!= ny1))
+            vmess("dimensions overruled: %d x %d y %d",ntmax,nxmax,nymax);
+        if(!getparfloat("dt", &dt)) dt=dt1;
+
+        fp_in1 = fopen(file_mute, "r");
+        if (fp_in1 == NULL) verr("error on opening input file_mute=%s", file_mute);
+
+        size = ntmax * nxmax *nymax;
+        tmpdata = (float *)malloc(size*sizeof(float));
+        hdrs_in1 = (segy *)calloc(nxmax*nymax,sizeof(segy));
+        
+        nx1,ny1 = readData3D(fp_in1, tmpdata, hdrs_in1, nt1);
+        nxy = nx1*ny1;
+        if (nxy == 0) {
+            fclose(fp_in1);
+            if (verbose) vmess("end of file_mute data reached");
+        }
+
+        if (verbose) {
+            disp_fileinfo3D(file_mute, nt1, nx1, ny1, ft1, fx1, fy1, dt, dx1, dy1, hdrs_in1);
+        }
+    }
+
+/* Reading input data for file_shot */
+
+    ngath = 1;
+    ret = getFileInfo3D(file_shot, &nt2, &nx2, &ny2, &ngath, &dt2, &dx2, &dy2, &ft2, &fx2, &fy2, &sclshot, &ntraces);
+    
+    if (!getparint("ntmax", &ntmax)) ntmax = nt2;
+    if (!getparint("nxmax", &nxmax)) nxmax = nx2;
+    if (!getparint("nymax", &nymax)) nymax = ny2;
+
+    size = ntmax * nxmax * nymax;
+    tmpdata2 = (float *)malloc(size*sizeof(float));
+    hdrs_in2 = (segy *)calloc(nxmax*nymax,sizeof(segy));
+
+    if (file_shot != NULL) fp_in2 = fopen(file_shot, "r");
+    else fp_in2=stdin;
+    if (fp_in2 == NULL) verr("error on opening input file_shot=%s", file_shot);
+    
+    nx2,ny2 = readData3D(fp_in2, tmpdata2, hdrs_in2, nt2);
+    nxy = nx2*ny2;
+    if (nxy == 0) {
+        fclose(fp_in2);
+        if (verbose) vmess("end of file_shot data reached");
+    }
+    nt2 = hdrs_in2[0].ns;
+    ft2 = hdrs_in2[0].f1;
+    fx2 = hdrs_in2[0].f2;
+    dt2 = (float)hdrs_in2[0].dt*1e-6;
+        
+    if (verbose) {
+        disp_fileinfo3D(file_shot, nt2, nx2, ny2, ft2, fx2, fy2, dt2, dx2, dy2, hdrs_in2);
+    }
+    
+    /* file_shot will be used as well to define the mute window */
+    if (file_mute == NULL) {
+        nx1=nx2;
+        nt1=nt2;
+        ny1=ny2;
+        dt=dt2;
+        ft1=ft2;
+        fx1=fx2;
+        fy1=fy2;
+        tmpdata = tmpdata2;
+        hdrs_in1 = hdrs_in2;
+        sclsxgx = sclshot;
+    }
+
+    if (verbose) vmess("sampling file_mute=%d, file_shot=%d", nt1, nt2);
+
+/*================ initializations ================*/
+
+    nxy    = nx1*ny1;
+    maxval = (int *)calloc(nxy,sizeof(int));
+    xrcv   = (int *)calloc(nxy,sizeof(int));
+    
+    if (file_out==NULL) fp_out = stdout;
+    else {
+        fp_out = fopen(file_out, "w+");
+        if (fp_out==NULL) verr("error on ceating output file");
+    }
+    if (check!=0){
+        fp_chk = fopen("check.su", "w+");
+        if (fp_chk==NULL) verr("error on ceating output file");
+        fp_psline1 = fopen("pslinepos.asci", "w+");
+        if (fp_psline1==NULL) verr("error on ceating output file");
+        fp_psline2 = fopen("pslineneg.asci", "w+");
+        if (fp_psline2==NULL) verr("error on ceating output file");
+        
+    }
+    if (smooth) {
+        costaper = (float *)malloc(smooth*sizeof(float));
+        scl = M_PI/((float)smooth);
+        for (i=0; i<smooth; i++) {
+            costaper[i] = 0.5*(1.0+cos((i+1)*scl));
+/*            fprintf(stderr,"costaper[%d]=%f\n",i,costaper[i]);*/
+        }
+    }
+
+/*================ loop over all shot records ================*/
+
+    k=1;
+    while (nxy > 0) {
+        if (verbose) vmess("processing input gather %d", k);
+
+/*================ loop over all shot records ================*/
+
+        /* find consistent (one event) maximum related to maximum value */
+        
+        /* find global maximum 
+        xmax=0.0;
+        for (i = 0; i < nx1; i++) {
+            tmax=0.0;
+            jmax = 0;
+            for (j = 0; j < nt1; j++) {
+                lmax = fabs(tmpdata[i*nt1+j]);
+                if (lmax > tmax) {
+                    jmax = j;
+                    tmax = lmax;
+                    if (lmax > xmax) {
+                        imax = i;
+                        xmax=lmax;
+                    }
+                }
+            }
+            maxval[i] = jmax;
+        }
+        */
+
+        /* alternative find maximum at source position */
+        dxrcv = (hdrs_in1[nxy-1].gx - hdrs_in1[0].gx)*sclsxgx/(float)(nx1-1);
+        dyrcv = (hdrs_in1[nxy-1].gy - hdrs_in1[0].gy)*sclsxgx/(float)(ny1-1);
+        ixmax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
+        iymax = NINT(((hdrs_in1[0].sy-hdrs_in1[0].gy)*sclsxgx)/dxrcv);
+        if (iymax > ny1-1) {
+            vmess("source of y is past array, snapping to nearest y");
+            iymax = ny1-1;
+        }
+        if (iymax < 0) {
+            vmess("source of y is before array, snapping to nearest y");
+            iymax = 0;
+        }
+        if (ixmax > nx1-1) {
+            vmess("source of x is past array, snapping to nearest x");
+            ixmax = nx1-1;
+        }
+        if (ixmax < 0) {
+            vmess("source of x is before array, snapping to nearest x");
+            ixmax = 0;
+        }
+        tmax=0.0;
+        jmax = 0;
+        xmax=0.0;
+        for (j = 0; j < nt1; j++) {
+            lmax = fabs(tmpdata[iymax*nx1*nt1+ixmax*nt1+j]);
+            if (lmax > tmax) {
+                jmax = j;
+                tmax = lmax;
+                   if (lmax > xmax) {
+                       xmax=lmax;
+                   }
+            }
+        }
+        maxval[iymax*nx1+ixmax] = jmax;
+        if (verbose >= 3) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, iymax, maxval[iymax*nx1+ixmax]);
+
+        /* search forward in x-trace direction from maximum in file */
+        for (i = ixmax+1; i < nx1; i++) {
+            tstart = MAX(0, (maxval[iymax*nx1+i-1]-hw));
+            tend   = MIN(nt1-1, (maxval[iymax*nx1+i-1]+hw));
+            jmax=tstart;
+            tmax=0.0;
+            for(j = tstart; j <= tend; j++) {
+                lmax = fabs(tmpdata[iymax*nx1*nt1+i*nt1+j]);
+                if (lmax > tmax) {
+                    jmax = j;
+                    tmax = lmax;
+                }
+            }
+            maxval[iymax*nx1+i] = jmax;
+        }
+        /* search backward in x-trace direction from maximum in file */
+        for (i = ixmax-1; i >=0; i--) {
+            tstart = MAX(0, (maxval[iymax*nx1+i+1]-hw));
+            tend   = MIN(nt1-1, (maxval[iymax*nx1+i+1]+hw));
+            jmax=tstart;
+            tmax=0.0;
+            for(j = tstart; j <= tend; j++) {
+                lmax = fabs(tmpdata[iymax*nx1*nt1+i*nt1+j]);
+                if (lmax > tmax) {
+                    jmax = j;
+                    tmax = lmax;
+                }
+            }
+            maxval[iymax*nx1+i] = jmax;
+        }
+
+        /* search forward in y-trace direction from maximum in file */
+        for (i = iymax+1; i < ny1; i++) {
+            tmax=0.0;
+            jmax = 0;
+            for (j = 0; j < nt1; j++) {
+                lmax = fabs(tmpdata[i*nx1*nt1+ixmax*nt1+j]);
+                if (lmax > tmax) {
+                    jmax = j;
+                    tmax = lmax;
+                }
+            }
+            maxval[i*nx1+ixmax] = jmax;
+            if (verbose >= 8) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, i, maxval[i*nx1+ixmax]);
+            /* search forward in x-trace direction from maximum in file */
+            for (l = ixmax+1; l < nx1; l++) {
+                tstart = MAX(0, (maxval[i*nx1+l-1]-hw));
+                tend   = MIN(nt1-1, (maxval[i*nx1+l-1]+hw));
+                jmax=tstart;
+                tmax=0.0;
+                for(j = tstart; j <= tend; j++) {
+                    lmax = fabs(tmpdata[i*nx1*nt1+l*nt1+j]);
+                    if (lmax > tmax) {
+                        jmax = j;
+                        tmax = lmax;
+                    }
+                }
+                maxval[i*nx1+l] = jmax;
+            }
+            /* search backward in x-trace direction from maximum in file */
+            for (l = ixmax-1; l >=0; l--) {
+                tstart = MAX(0, (maxval[i*nx1+l+1]-hw));
+                tend   = MIN(nt1-1, (maxval[i*nx1+l+1]+hw));
+                jmax=tstart;
+                tmax=0.0;
+                for(j = tstart; j <= tend; j++) {
+                    lmax = fabs(tmpdata[i*nx1*nt1+l*nt1+j]);
+                    if (lmax > tmax) {
+                        jmax = j;
+                        tmax = lmax;
+                    }
+                }
+                maxval[i*nx1+l] = jmax;
+            }
+        }
+
+        /* search backward in y-trace direction from maximum in file */
+        for (i = iymax-1; i >= 0; i--) {
+            tmax=0.0;
+            jmax = 0;
+            for (j = 0; j < nt1; j++) {
+                lmax = fabs(tmpdata[i*nx1*nt1+ixmax*nt1+j]);
+                if (lmax > tmax) {
+                    jmax = j;
+                    tmax = lmax;
+                }
+            }
+            maxval[i*nx1+ixmax] = jmax;
+            if (verbose >= 8) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, i, maxval[i*nx1+ixmax]);
+            /* search forward in x-trace direction from maximum in file */
+            for (l = ixmax+1; l < nx1; l++) {
+                tstart = MAX(0, (maxval[i*nx1+l-1]-hw));
+                tend   = MIN(nt1-1, (maxval[i*nx1+l-1]+hw));
+                jmax=tstart;
+                tmax=0.0;
+                for(j = tstart; j <= tend; j++) {
+                    lmax = fabs(tmpdata[i*nx1*nt1+l*nt1+j]);
+                    if (lmax > tmax) {
+                        jmax = j;
+                        tmax = lmax;
+                    }
+                }
+                maxval[i*nx1+l] = jmax;
+            }
+            /* search backward in x-trace direction from maximum in file */
+            for (l = ixmax-1; l >=0; l--) {
+                tstart = MAX(0, (maxval[i*nx1+l+1]-hw));
+                tend   = MIN(nt1-1, (maxval[i*nx1+l+1]+hw));
+                jmax=tstart;
+                tmax=0.0;
+                for(j = tstart; j <= tend; j++) {
+                    lmax = fabs(tmpdata[i*nx1*nt1+l*nt1+j]);
+                    if (lmax > tmax) {
+                        jmax = j;
+                        tmax = lmax;
+                    }
+                }
+                maxval[i*nx1+l] = jmax;
+            }
+        }
+
+/* scale with maximum ampltiude */
+
+        if (scale==1) {
+            for (l = 0; l < ny2; l++) {
+                for (i = 0; i < nx2; i++) {
+                    lmax = fabs(tmpdata2[l*nx2*nt2+i*nt2+maxval[i]]);
+                    for (j = 0; j < nt2; j++) {
+                        tmpdata2[l*nx2*nt2+i*nt2+j] = tmpdata2[l*nx2*nt2+i*nt2+j]/lmax;
+                    }
+                }
+            }
+        }
+
+        for (l = 0; l < ny2; l++) {
+            for (i = 0; i < nx2; i++) {
+                xrcv[l*nx2+i] = l*nx2+i;
+            }
+        }
+
+/*================ apply mute window ================*/
+        
+        applyMute(tmpdata2, maxval, smooth, above, 1, nx2*ny2, nt2, xrcv, nx2*ny2, shift);
+
+/*================ write result to output file ================*/
+
+        ret = writeData(fp_out, tmpdata2, hdrs_in2, nt2, nx2*ny2);
+        if (ret < 0 ) verr("error on writing output file.");
+
+        /* put mute window in file to check correctness of mute */
+        if (check !=0) {
+            for (l=0; l<ny1; l++) {
+                for (i = 0; i < nx1; i++) {
+                    jmax = maxval[l*nx1+i]-shift;
+                    tmpdata[l*nx1*nt1+i*nt1+jmax] = 2*xmax;
+                }
+            }
+            if (above==0){
+                for (l=0; l<ny1; l++) {
+                    for (i = 0; i < nx1; i++) {
+                        jmax = nt2-maxval[l*nx1+i]+shift;
+                        tmpdata[l*nx1*nt1+i*nt1+jmax] = 2*xmax;
+                    }
+                }
+            }
+            ret = writeData(fp_chk, tmpdata, hdrs_in1, nt1, nx1*ny1);
+            if (ret < 0 ) verr("error on writing check file.");
+            for (l=0; l<ny1; l++) {
+                for (i=0; i<nx1; i++) {
+                    jmax = maxval[l*nx1+i]-shift;
+                    ret = fprintf(fp_psline1, "%.5f %.5f \n",jmax*dt,hdrs_in1[l*nx1+i].gx*sclshot,hdrs_in1[l*nx1+i].gy*sclshot);
+                    jmax =-maxval[l*nx1+i]+shift;
+                    ret = fprintf(fp_psline2, "%.5f %.5f \n",jmax*dt,hdrs_in1[l*nx1+i].gx*sclshot,hdrs_in1[l*nx1+i].gy*sclshot);
+                }
+            }
+        }
+
+/*================ Read next record for muting ================*/
+
+        if (file_mute != NULL) {    
+            nx1, ny1 = readData3D(fp_in1, tmpdata, hdrs_in1, nt1);
+            nxy = nx1*ny1;
+            if (nxy == 0) {
+                fclose(fp_in1);
+                if (verbose) vmess("end of file_mute data reached");
+                fclose(fp_in2);
+                if (fp_out!=stdout) fclose(fp_out);
+                if (check!=0) fclose(fp_chk);
+                if (check!=0) {
+                    fclose(fp_psline1);
+                    fclose(fp_psline2);
+                }
+                break;
+            }
+            nt1 = (int)hdrs_in1[0].ns;
+            if (nt1 > ntmax) verr("n_samples (%d) greater than ntmax", nt1);
+            if (nx1 > nxmax) verr("n_traces  (%d) greater than nxmax", nx1);
+            if (ny1 > nymax) verr("n_traces  (%d) greater than nymax", ny1);
+            if (verbose) {
+                disp_fileinfo3D(file_mute, nt1, nx1, ny1, ft1, fx1, fy1, dt, dx1, dy1, hdrs_in1);
+            }
+        }
+
+/*================ Read next shot record(s) ================*/
+
+        nx2,ny2 = readData3D(fp_in2, tmpdata2, hdrs_in2, nt2);
+        nxy = nx2*ny2;
+        if (nxy == 0) {
+            if (verbose) vmess("end of file_shot data reached");
+            fclose(fp_in2);
+            break;
+        }
+        nt2 = (int)hdrs_in2[0].ns;
+        if (nt2 > ntmax) verr("n_samples (%d) greater than ntmax", nt2);
+        if (nx2 > nxmax) verr("n_traces  (%d) greater than nxmax", nx2);
+        if (ny2 > nymax) verr("n_traces  (%d) greater than nymax", ny2);
+        if (verbose) {
+            disp_fileinfo3D(file_shot, nt2, nx2, ny2, ft2, fx2, fy2, dt, dx2, dy2, hdrs_in2);
+        }
+
+        if (file_mute == NULL) {
+            nx1=nx2;
+            ny1=ny2;
+            nt1=nt2;
+            hdrs_in1 = hdrs_in2;
+            tmpdata = tmpdata2;
+        }
+
+        k++;
+    }
+
+    t1 = wallclock_time();
+    if (verbose) vmess("Total CPU-time = %f",t1-t0);
+    
+
+    return 0;
+}
\ No newline at end of file

marchenko3D/readData3D.c

+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "segy.h"
+
+/**
+* reads SU file and returns header and 2D array
+*
+*   AUTHOR:
+*           Jan Thorbecke (janth@xs4all.nl)
+*           The Netherlands 
+**/
+
+
+int readData3D(FILE *fp, float *data, segy *hdrs, int n1)
+{
+	size_t nread;
+	int oneshot, itrace, sx, sy, fldr, gy, nx, ny;
+	segy hdr;
+
+	oneshot = 1;
+	itrace  = 0;
+    ny      = 1;
+
+	nread = fread(&hdrs[0], 1, TRCBYTES, fp);
+	if (nread == 0) return 0; /* end of file reached */
+
+	if (n1==0) n1 = hdrs[0].ns;
+	sx   = hdrs[0].sx;
+    sy   = hdrs[0].sy;
+	fldr = hdrs[0].fldr;
+    gy   = hdrs[0].gy;
+	while (oneshot) {
+		nread = fread(&data[itrace*n1], sizeof(float), n1, fp);
+		assert (nread == n1);
+		itrace++;
+		nread = fread(&hdr, 1, TRCBYTES, fp);
+		if (nread == 0) break;
+		assert(nread == TRCBYTES);
+		if ( (sx != hdr.sx) || (sy != hdr.sy) || (fldr != hdr.fldr)) { /* end of shot record */
+			fseek( fp, -TRCBYTES, SEEK_CUR );
+			break;
+		}
+		if ( (gy != hdr.gy)) {
+            ny++;
+            gy = hdr.gy;
+        }
+		memcpy(&hdrs[itrace], &hdr, TRCBYTES);
+	}
+    nx = itrace/ny;
+
+	return nx, ny;
+}

marchenko3D/synthesis3D.c

+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+int omp_get_max_threads(void);
+int omp_get_num_threads(void);
+void omp_set_num_threads(int num_threads);
+
+#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))
+int compareInt(const void *a, const void *b) 
+{ return (*(int *)a-*(int *)b); }
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+    float r,i;
+} complex;
+#endif/* complex */
+
+void synthesisPositions3D(int nx, int ny, int nxs, int nys, int Nfoc, float *xrcv, float *yrcv,
+float *xsrc, float *ysrc, int *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys,
+int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose)
+{
+    int     j, l, ixsrc, iysrc, isrc, k, *count, nxy;
+    float   fxb, fxe, fyb, fye;
+
+    if (fxsb < 0) fxb = 1.001*fxsb;
+    else          fxb = 0.999*fxsb;
+    if (fysb < 0) fyb = 1.001*fysb;
+    else          fyb = 0.999*fysb;
+    if (fxse > 0) fxe = 1.001*fxse;
+    else          fxe = 0.999*fxse;
+    if (fyse > 0) fye = 1.001*fyse;
+    else          fye = 0.999*fyse;
+
+    nxy = nx*ny;
+
+    count   = (int *)calloc(nxs*nys,sizeof(int)); // number of traces that contribute to the integration over x
+
+/*================ SYNTHESIS ================*/
+
+    for (l = 0; l < 1; l++) { /* assuming all focal operators cover the same lateral area */
+//    for (l = 0; l < Nfoc; l++) {
+        *npos=0;
+
+        if (reci == 0 || reci == 1) {
+            for (k=0; k<nshots; k++) {
+
+                ixsrc = NINT((xsrc[k] - fxsb)/dxs);
+                iysrc = NINT((ysrc[k] - fysb)/dys);
+                isrc  = iysrc*nxs + ixsrc;
+                if (verbose>=3) {
+                    vmess("source position:         x=%.2f y=%.2f in operator x=%d y=%d pos=%d", xsrc[k], ysrc[k], ixsrc, iysrc, isrc);
+                    vmess("receiver positions:      x:%.2f <--> %.2f y:%.2f <--> %.2f", xrcv[k*nxy+0], xrcv[k*nxy+nxy-1], yrcv[k*nxy+0], yrcv[k*nxy+nxy-1]);
+                    vmess("focal point positions:   x:%.2f <--> %.2f y:%.2f <--> %.2f", fxsb, fxse, fysb, fyse);
+                }
+        
+                if ((NINT(xsrc[k]-fxse) > 0)           || (NINT(xrcv[k*nxy+nxy-1]-fxse) > 0) ||
+                    (NINT(xrcv[k*nxy+nxy-1]-fxsb) < 0) || (NINT(xsrc[k]-fxsb) < 0)           || 
+                    (NINT(xrcv[k*nxy+0]-fxsb) < 0)     || (NINT(xrcv[k*nxy+0]-fxse) > 0)     || 
+                    (NINT(ysrc[k]-fyse) > 0)           || (NINT(yrcv[k*nxy+nxy-1]-fyse) > 0) ||
+                    (NINT(yrcv[k*nxy+nxy-1]-fysb) < 0) || (NINT(ysrc[k]-fysb) < 0)           || 
+                    (NINT(yrcv[k*nxy+0]-fysb) < 0)     || (NINT(yrcv[k*nxy+0]-fyse) > 0)       ) {
+                    vwarn("source/receiver positions are outside synthesis aperture");
+                    vmess("xsrc = %.2f xrcv_1 = %.2f xrvc_N = %.2f", xsrc[k], xrcv[k*nxy+0], xrcv[k*nxy+nxy-1]);
+                    vmess("ysrc = %.2f yrcv_1 = %.2f yrvc_N = %.2f", ysrc[k], yrcv[k*nxy+0], yrcv[k*nxy+nxy-1]);
+                    vmess("source position x:       %.2f in operator %d", xsrc[k], ixsrc);
+                    vmess("source position y:       %.2f in operator %d", ysrc[k], iysrc);
+                    vmess("receiver positions x:    %.2f <--> %.2f", xrcv[k*nxy+0], xrcv[k*nxy+nxy-1]);
+                    vmess("receiver positions y:    %.2f <--> %.2f", yrcv[k*nxy+0], yrcv[k*nxy+nxy-1]);
+                    vmess("focal point positions x: %.2f <--> %.2f", fxsb, fxse);
+                    vmess("focal point positions y: %.2f <--> %.2f", fysb, fyse);
+                }
+                
+                if ( (xsrc[k] >= fxb) && (xsrc[k] <= fxe) &&
+                     (ysrc[k] >= fyb) && (ysrc[k] <= fye) ) {
+                    
+				    j = linearsearch(ixpos, *npos, isrc);
+				    if (j < *npos) { /* the position (at j) is already included */
+					    count[j] += xnx[k];
+				    }
+				    else { /* add new postion */
+            		    ixpos[*npos] =  isrc;
+					    count[*npos] += xnx[k];
+                   	    *npos += 1;
+				    }
+    //                vmess("source position %d is inside synthesis model %f *npos=%d count=%d", k, xsrc[k], *npos, count[*npos]);
+			    }
+    
+    	    } /* end of nshots (k) loop */
+   	    } /* end of reci branch */
+    } /* end of Nfoc loop */
+
+    if (verbose>=4) {
+	    for (j=0; j < *npos; j++) { 
+            vmess("ixpos[%d] = %d count=%d", j, ixpos[j], count[j]);
+		}
+    }
+    free(count);
+
+/* sort ixpos into increasing values */
+    qsort(ixpos, *npos, sizeof(int), compareInt);
+
+
+    return;
+}
+
+int linearsearch(int *array, size_t N, int value)
+{
+	int j;
+/* Check is position is already in array */
+    j = 0;
+    while (j < N && value != array[j]) {
+        j++;
+    }
+	return j;
+}
+
+/*================ Convolution and Integration ================*/
+
+void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int ny, int nt, int nxs, int nys, int nts, float dt, float *xsyn, float *ysyn, 
+int Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, int *xnx, float fxse, float fxsb, float fyse, float fysb, float dxs, float dys, float dxsrc, 
+float dysrc, float dx, float dy, int ntfft, int nw, int nw_low, int nw_high,  int mode, int reci, int nshots, int nxsrc, int nysrc, 
+int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc,  int *reci_xrcv, float *ixmask, int verbose)
+{
+    int     nfreq, size, inx;
+    float   scl;
+    int     i, j, l, m, iw, ix, k, isrc, il, ik, nxy, nxys;
+    float   *rtrace, idxs, idys;
+    complex *sum, *ctrace;
+    int     npe;
+    static int first=1, *ircv;
+    static double t0, t1, t;
+
+    nxy     = nx*ny;
+    nxys    = nxs*nys;
+
+    size  = nxys*nts;
+    nfreq = ntfft/2+1;
+    /* scale factor 1/N for backward FFT,
+     * scale dt for correlation/convolution along time, 
+     * scale dx (or dxsrc) for integration over receiver (or shot) coordinates */
+    scl   = 1.0*dt/((float)ntfft);
+
+#ifdef _OPENMP
+    npe   = omp_get_max_threads();
+    /* parallelisation is over number of virtual source positions (Nfoc) */
+    if (npe > Nfoc) {
+        vmess("Number of OpenMP threads set to %d (was %d)", Nfoc, npe);
+        omp_set_num_threads(Nfoc);
+    }
+#endif
+
+    t0 = wallclock_time();
+
+    /* reset output data to zero */
+    memset(&iRN[0], 0, Nfoc*nxys*nts*sizeof(float));
+    ctrace = (complex *)calloc(ntfft,sizeof(complex));
+
+    if (!first) {
+    /* transform muted Ni (Top) to frequency domain, input for next iteration  */
+        for (l = 0; l < Nfoc; l++) {
+            /* set Fop to zero, so new operator can be defined within ixpos points */
+            memset(&Fop[l*nxys*nw].r, 0, nxys*nw*2*sizeof(float));
+            for (i = 0; i < npos; i++) {
+                   rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
+                   ix = ixpos[i];
+                   for (iw=0; iw<nw; iw++) {
+                       Fop[l*nxys*nw+iw*nxys+ix].r = ctrace[nw_low+iw].r;
+                       Fop[l*nxys*nw+iw*nxys+ix].i = mode*ctrace[nw_low+iw].i;
+                   }
+            }
+        }
+    }
+    else { /* only for first call to synthesis using all nxs traces in G_d */
+    /* transform G_d to frequency domain, over all nxs traces */
+        first=0;
+        for (l = 0; l < Nfoc; l++) {
+            /* set Fop to zero, so new operator can be defined within all ix points */
+            memset(&Fop[l*nxys*nw].r, 0, nxys*nw*2*sizeof(float));
+            for (i = 0; i < nxys; i++) {
+                   rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
+                   for (iw=0; iw<nw; iw++) {
+                       Fop[l*nxys*nw+iw*nxys+i].r = ctrace[nw_low+iw].r;
+                       Fop[l*nxys*nw+iw*nxys+i].i = mode*ctrace[nw_low+iw].i;
+                   }
+            }
+        }
+        idxs = 1.0/dxs;
+        idys = 1.0/dys;
+        ircv = (int *)malloc(nshots*nxy*sizeof(int));
+        for (k=0; k<nshots; k++) {
+            for (i = 0; i < nxy; i++) {
+                ircv[k*nxy+i] = NINT((yrcv[k*nxy+i]-fysb)*idys)*nx+NINT((xrcv[k*nxy+i]-fxsb)*idxs);
+            }
+        }
+    }
+    free(ctrace);
+    t1 = wallclock_time();
+    *tfft += t1 - t0;
+
+/* Loop over total number of shots */
+    if (reci == 0 || reci == 1) {
+        for (k=0; k<nshots; k++) {
+            if ((xsrc[k] < 0.999*fxsb) || (xsrc[k] > 1.001*fxse) || (ysrc[k] < 0.999*fysb) || (ysrc[k] > 1.001*fyse)) continue;
+            isrc = NINT((ysrc[k] - fysb)/dys)*nxs+NINT((xsrc[k] - fxsb)/dxs);
+            inx = xnx[k]; /* number of traces per shot */
+
+/*================ SYNTHESIS ================*/
+
+#pragma omp parallel default(none) \
+ shared(iRN, dx, dy, npe, nw, verbose) \
+ shared(Refl, Nfoc, reci, xrcv, xsrc, yrcv, ysrc, xsyn, ysyn) \
+ shared(fxsb, fxse, fysb, fyse, nxs, nys, nxys, dxs, dys) \
+ shared(nx, ny, nxy, dysrc, dxsrc, inx, k, nfreq, nw_low, nw_high) \
+ shared(Fop, size, nts, ntfft, scl, ircv, isrc) \
+ private(l, ix, j, m, i, sum, rtrace)
+{ /* start of parallel region */
+            sum   = (complex *)malloc(nfreq*sizeof(complex));
+            rtrace = (float *)calloc(ntfft,sizeof(float));
+
+#pragma omp for schedule(guided,1)
+            for (l = 0; l < Nfoc; l++) {
+		        /* compute integral over receiver positions */
+                /* multiply R with Fop and sum over nx */
+                memset(&sum[0].r,0,nfreq*2*sizeof(float));
+                for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
+                    for (i = 0; i < inx; i++) {
+                        ix = ircv[k*nxy+i];
+                        sum[j].r += Refl[k*nw*nxy+m*nxy+i].r*Fop[l*nw*nxys+m*nxys+ix].r -
+                                    Refl[k*nw*nxy+m*nxy+i].i*Fop[l*nw*nxys+m*nxys+ix].i;
+                        sum[j].i += Refl[k*nw*nxy+m*nxy+i].i*Fop[l*nw*nxys+m*nxys+ix].r +
+                                    Refl[k*nw*nxy+m*nxy+i].r*Fop[l*nw*nxys+m*nxys+ix].i;
+                    }
+                }
+
+                /* transfrom result back to time domain */
+                cr1fft(sum, rtrace, ntfft, 1);
+
+                /* place result at source position ixsrc; dx = receiver distance */
+                for (j = 0; j < nts; j++) 
+                    iRN[l*size+isrc*nts+j] += rtrace[j]*scl*dx*dy;
+            
+            } /* end of parallel Nfoc loop */
+            free(sum);
+            free(rtrace);
+
+#ifdef _OPENMP
+#pragma omp single 
+            npe   = omp_get_num_threads();
+#endif
+} /* end of parallel region */
+
+        if (verbose>4) vmess("*** Shot gather %d processed ***", k);
+
+        } /* end of nshots (k) loop */
+    }     /* end of if reci */
+
+    t = wallclock_time() - t0;
+    if (verbose) {
+        vmess("OMP: parallel region = %f seconds (%d threads)", t, npe);
+    }
+
+    return;
+}
\ No newline at end of file

marchenko3D/test.c

-#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 */
-
-void cr1fft(complex *cdata, REAL *rdata, int n, int sign);
-int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, float *sclsxgxsygy, int *nxm);
-int disp_fileinfo3D(char *file, int n1, int n2, int n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
-int readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots, int nx, int ny, int ntfft, int mode, float scale, int verbose);
-int unique_elements(float *arr, int len);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, int *xnx, int Nfoc, int nx, int ny, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
-void synthesisPosistions3D(int nx, int ny, int nxs, int nys, int Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, int *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys, int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose);
-void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int ny, int nt, int nxs, int nys, int nts, float dt, float *xsyn, float *ysyn, 
-int Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, int *xnx, float fxse, float fxsb, float fyse, float fysb, float dxs, float dys, float dxsrc, 
-float dysrc, float dx, float dy, int ntfft, int nw, int nw_low, int nw_high,  int mode, int reci, int nshots, int nxsrc, int nysrc, 
-int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc,  int *reci_xrcv, float *ixmask, int verbose)
-
-char *sdoc[] = {
-" ",
-" test 3D - Test 3D functions ",
-" ",
-"fin=.......................... File name of input",
-" ",
-" authors  : Joeri Brackenhoff 	(J.A.Brackenhoff@tudelft.nl)",
-"		   : Jan Thorbecke		(janth@xs4all.nl)",
-NULL};
-
-int main (int argc, char **argv)
-{
-	char    *file_shot, *file_out, *file_tinv, *file_shot_out;
-    FILE    *fp_out;
-	float   dts, dxs, dys, fts, fxs, fys, scl, fmin, fmax;
-    float   dtd, dxd, dyd, ftd, fxd, fyd, scld;
-    float   *xrcv, *yrcv, *xsrc, *ysrc, *zsrc;
-    float   *xrcvsyn, *yrcvsyn, *xsyn, *ysyn, *zsyn, *f1d, *iRN;
-    float   fxsb, fysb, fxse, fyse;
-    complex *Refl;
-	int     nts, nxs, nys, nxys, nshots, ntrs, ret, *xnx, *xnxsyn;
-    int     ntd, nxd, nyd, nxyd, Nfoc, ntrd, *muteW;
-    int     ntfft, nfreq, nw_low, nw_high, nw, mode, verbose;
-    int     i, j, l, it, nxsrc, nysrc, hw, shift, smooth, tracf;
-    int     *ixpos, npos;
-    segy    *hdrs_out, *hdrs_shot;
-
-	initargs(argc, argv);
-	requestdoc(1);
-
-	if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
-    if (!getparstring("file_out", &file_out)) file_out = "out.su";
-    if (!getparstring("file_shot_out", &file_shot_out)) file_shot_out = "shot_out.su";
-    if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
-    if (!getparfloat("fmin", &fmin)) fmin = 0.0;
-    if (!getparfloat("fmax", &fmax)) fmax = 70.0;
-    if (!getparint("hw", &hw)) hw = 15;
-    if (!getparint("smooth", &smooth)) smooth = 5;
-    if (!getparint("shift", &shift)) shift = 5;
-	if (file_shot == NULL) verr("Give me a shot name for the love of God");
-    if (file_tinv == NULL) verr("Give me a tinv name for the love of God");
-
-    nshots=0;
-    mode=1;
-    verbose=10;
-        
-    ret = getFileInfo3D(file_shot, &nts, &nxs, &nys, &nshots, &dts, &dxs, &dys, &fts, &fxs, &fys, &scl, &ntrs);
-    vmess("nts=%d nxs=%d nys=%d nshots=%d ntrs=%d",nts,nxs,nys,nshots,ntrs);
-    vmess("dts=%.5f dxs=%.5f dys=%.5f fts=%.5f fxs=%.5f fys=%.5f scl=%.5f",dts,dxs,dys,fts,fxs,fys,scl);
-    
-    ntfft = optncr(nts); 
-    nfreq = ntfft/2+1;
-    nw_low = (int)MIN((fmin*ntfft*dts), nfreq-1);
-    nw_low = MAX(nw_low, 1);
-    nw_high = MIN((int)(fmax*ntfft*dts), nfreq-1);
-    nw  = nw_high - nw_low + 1;
-    scl   = 1.0/((float)ntfft);
-
-    Refl    = (complex *)malloc(nw*nxs*nys*nshots*sizeof(complex));   // reflection response in frequency domain
-    xrcv    = (float *)calloc(nshots*nxs*nys,sizeof(float));          // x-rcv postions of shots
-    xsrc    = (float *)calloc(nshots,sizeof(float));                // x-src position of shots
-    yrcv    = (float *)calloc(nshots*nxs*nys,sizeof(float));          // y-rcv postions of shots
-    ysrc    = (float *)calloc(nshots,sizeof(float));                // y-src position of shots
-    zsrc    = (float *)calloc(nshots,sizeof(float));                // z-src position of shots
-    xnx     = (int *)calloc(nshots,sizeof(int));                    // number of traces per shot
-
-    readShotData3D(file_shot, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw, nw_low, nshots, nxs, nys, ntfft, mode, scl, verbose);
-
-    for (i=0;i<nshots;i++) {
-        vmess("xsrc=%.3f ysrc=%.3f zsrc=%.3f xrcv1=%.3f xrcv2=%.3f yrcv1=%.3f yrcv2=%.3f",xsrc[i],ysrc[i],zsrc[i],xrcv[i*nxs*nys],xrcv[(i+1)*nxs*nys-1],yrcv[i*nxs*nys],yrcv[(i+1)*nxs*nys-1]);
-    }
-
-    nxsrc = unique_elements(xsrc,nshots);
-    nysrc = unique_elements(ysrc,nshots);
-
-    vmess("nxsrc=%d nysrc=%d",nxsrc,nysrc);
-
-    ret = getFileInfo3D(file_tinv, &ntd, &nxd, &nyd, &Nfoc, &dtd, &dxd, &dyd, &ftd, &fxd, &fyd, &scld, &ntrd);
-    vmess("ntd=%d nxd=%d nyd=%d Nfoc=%d ntrd=%d",ntd,nxd,nyd,Nfoc,ntrd);
-    vmess("dtd=%.5f dxd=%.5f dyd=%.5f ftd=%.5f fxd=%.5f fyd=%.5f scld=%.5f",dtd,dxd,dyd,ftd,fxd,fyd,scld);
-
-    f1d     = (float *)calloc(Nfoc*nxd*nyd*ntfft,sizeof(float));
-    muteW   = (int *)calloc(Nfoc*nxd*nyd,sizeof(int));
-    xrcvsyn = (float *)calloc(Nfoc*nxd*nyd,sizeof(float)); // x-rcv postions of focal points
-    yrcvsyn = (float *)calloc(Nfoc*nxd*nyd,sizeof(float)); // x-rcv postions of focal points
-    xsyn    = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points
-    ysyn    = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points
-    zsyn    = (float *)malloc(Nfoc*sizeof(float)); // z-src position of focal points
-    xnxsyn  = (int *)calloc(Nfoc,sizeof(int)); // number of traces per focal point
-    ixpos   = (int *)calloc(nxd,sizeof(int)); // x-position of source of shot in G_d domain (nxd with dxd)
-
-    mode=-1;
-    readTinvData3D(file_tinv, xrcvsyn, yrcvsyn, xsyn, ysyn, zsyn, xnxsyn, Nfoc, nxd, nyd, ntfft, mode, muteW, f1d, hw, verbose);
-
-    fxsb = fxd;
-    fysb = fyd;
-    if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxsb = xrcvsyn[0];
-    fxse = fxsb + (float)(nxd-1)*dxd;
-    if (yrcvsyn[0] != 0 || yrcvsyn[1] != 0 ) fysb = yrcvsyn[0];
-    fyse = fysb + (float)(nyd-1)*dyd;
-
-    vmess("fxsb=%.3f fysb=%.3f fxse=%.3f fyse=%.3f",fxsb,fysb,fxse,fyse);
-
-    synthesisPosistions3D(nxs, nys, nxd, nyd, Nfoc, xrcv, yrcv,
-        xsrc, ysrc, xnx, fxse, fyse, fxsb, fysb, dxs, dys,
-        nshots, nxsrc, nysrc, ixpos, &npos, 0, verbose);
-
-    vmess("npos:%d",npos);
-
-    iRN     = (float *)calloc(Nfoc*nxd*nyd*ntfft,sizeof(float));
-    nxys    = nxs*nys;
-    nxyd    = nxd*nyd;
-
-    synthesis3D(Refl, Fop, f1d, iRN, nxs, nys, nts, nxd, nyd, ntd, dts, xsyn, ysyn,
-            Nfoc, xrcv, yrcv, xsrc, ysrc, xnx, fxse, fxsb, fyse, fysb, dxs, dys, dxsrc, dysrc, dx, dy, ntfft, nw, nw_low, nw_high, mode,
-            0, nshots, nxsrc, nysrc, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
-
-
-    hdrs_out     = (segy *)calloc(Nfoc*nxd*nyd,sizeof(segy));
-
-    tracf = 1;
-
-    for (l = 0; l < Nfoc; l++) {
-        for (j = 0; j < nyd; j++) {
-            for (i = 0; i < nxd; i++) {
-                hdrs_out[l*nyd*nxd+j*nxd+i].ns      = ntfft;
-                hdrs_out[l*nyd*nxd+j*nxd+i].fldr    = l+1;
-                hdrs_out[l*nyd*nxd+j*nxd+i].trid    = 1;
-                hdrs_out[l*nyd*nxd+j*nxd+i].dt      = dtd*1000000;
-                hdrs_out[l*nyd*nxd+j*nxd+i].f1      = ftd;
-                hdrs_out[l*nyd*nxd+j*nxd+i].f2      = fxd;
-                hdrs_out[l*nyd*nxd+j*nxd+i].d1      = dtd;
-                hdrs_out[l*nyd*nxd+j*nxd+i].d2      = dxd;
-                hdrs_out[l*nyd*nxd+j*nxd+i].trwf    = nxd*nyd;
-                hdrs_out[l*nyd*nxd+j*nxd+i].scalco  = -1000;
-                hdrs_out[l*nyd*nxd+j*nxd+i].sx      = NINT(1000*(xsyn[l]));
-                hdrs_out[l*nyd*nxd+j*nxd+i].sy      = NINT(1000*(ysyn[l]));
-                hdrs_out[l*nyd*nxd+j*nxd+i].gx      = NINT(1000*(xrcvsyn[l*nyd*nxd+j*nxd+i]));
-                hdrs_out[l*nyd*nxd+j*nxd+i].gy      = NINT(1000*(yrcvsyn[l*nyd*nxd+j*nxd+i]));
-                hdrs_out[l*nyd*nxd+j*nxd+i].scalel  = -1000;
-                hdrs_out[l*nyd*nxd+j*nxd+i].tracl   = i+1;
-                hdrs_out[l*nyd*nxd+j*nxd+i].offset  = (long)NINT(xrcvsyn[l*nyd*nxd+j*nxd+i] - xsyn[l]);
-                hdrs_out[l*nyd*nxd+j*nxd+i].tracf   = tracf++;
-                hdrs_out[l*nyd*nxd+j*nxd+i].selev   = NINT(zsyn[l]*1000);
-                hdrs_out[l*nyd*nxd+j*nxd+i].sdepth  = NINT(-zsyn[l]*1000);
-            }
-        }
-    }
-
-    fp_out = fopen(file_out,"w+");
-    ret = writeData(fp_out, iRN, hdrs_out, ntfft, Nfoc*nxd*nyd);
-    if (ret < 0 ) verr("error on writing output file.");
-    fclose(fp_out);
-
-    free(Refl);free(xrcv);free(yrcv);free(xsrc);free(ysrc);free(zsrc);free(xnx);
-    free(f1d);free(xrcvsyn);free(yrcvsyn);free(xsyn);free(ysyn);free(zsyn);free(xnxsyn);free(muteW);
-
-	return;
-}
-
-
-int unique_elements(float *arr, int len)
-{
-     if (len <= 0) return 0;
-     int unique = 1;
-     int outer, inner, is_unique;
-
-     for (outer = 1; outer < len; ++outer)
-     {
-        is_unique = 1;
-        for (inner = 0; is_unique && inner < outer; ++inner)
-        {  
-             if ((arr[inner] >= arr[outer]-1.0) && (arr[inner] <= arr[outer]+1.0)) is_unique = 0;
-        }
-        if (is_unique) ++unique;
-     }
-     return unique;
-}
\ No newline at end of file