diff --git a/marchenko3D/demo/test/test.scr b/marchenko3D/demo/test/test.scr
new file mode 100755
index 0000000000000000000000000000000000000000..af4ff97379d801e03f38fd63eefcf28235ff5e3a
--- /dev/null
+++ b/marchenko3D/demo/test/test.scr
@@ -0,0 +1,5 @@
+#!/bin/bash
+
+./../../test file_shot=shot_gy.su
+./../../test file_shot=shot_sxgy.su
+./../../test file_shot=shot_sxsygy.su
diff --git a/marchenko3D/getFileInfo3D.c b/marchenko3D/getFileInfo3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..814d816d67fe2375d1dd7c7bebf1a400df85be34
--- /dev/null
+++ b/marchenko3D/getFileInfo3D.c
@@ -0,0 +1,216 @@
+#define _FILE_OFFSET_BITS 64
+#define _LARGEFILE_SOURCE
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <math.h>
+#include "segy.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* gets sizes, sampling and min/max values of a SU file
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+void vmess(char *fmt, ...);
+void verr(char *fmt, ...);
+int optncr(int n);
+
+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)
+{
+ FILE *fp;
+ size_t nread, data_sz;
+ off_t bytes, ret, trace_sz, ntraces;
+ int sx_shot, sy_shot, gx_start, gx_end, gy_start, gy_end, itrace, one_shot, igath, end_of_file, fldr_shot;
+ int verbose=1, igy, nsx, nsy;
+ float scl, *trace, dxsrc, dxrcv, dysrc, dyrcv;
+ segy hdr;
+
+ if (filename == NULL) { /* read from input pipe */
+ *n1=0;
+ *n2=0;
+ *n3=0;
+ return -1; /* Input pipe */
+ }
+ else fp = fopen( filename, "r" );
+ if (fp == NULL) verr("File %s does not exist or cannot be opened", filename);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ ret = fseeko( fp, 0, SEEK_END );
+ if (ret<0) perror("fseeko");
+ bytes = ftello( fp );
+
+ *n1 = hdr.ns;
+ if ( (hdr.trid == 1) && (hdr.dt != 0) ) {
+ *d1 = ((float) hdr.dt)*1.e-6;
+ *f1 = ((float) hdr.delrt)/1000.;
+ }
+ else {
+ *d1 = hdr.d1;
+ *f1 = hdr.f1;
+ }
+ *f2 = hdr.f2;
+
+ data_sz = sizeof(float)*(*n1);
+ trace_sz = sizeof(float)*(*n1)+TRCBYTES;
+ ntraces = (int) (bytes/trace_sz);
+
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+
+ *sclsxgxsygy = scl;
+ /* check to find out number of traces in shot gather */
+
+ one_shot = 1;
+ itrace = 0;
+ igy = 1;
+ fldr_shot = hdr.fldr;
+ sx_shot = hdr.sx;
+ sy_shot = hdr.sy;
+ gx_start = hdr.gx;
+ gy_start = hdr.gy;
+ gy_end = gy_start;
+ trace = (float *)malloc(hdr.ns*sizeof(float));
+ fseeko( fp, TRCBYTES, SEEK_SET );
+
+ while (one_shot) {
+ nread = fread( trace, sizeof(float), hdr.ns, fp );
+ assert (nread == hdr.ns);
+ itrace++;
+ if (hdr.gy != gy_end) {
+ gy_end = hdr.gy;
+ igy++;
+ }
+ gx_end = hdr.gx;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread==0) break;
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr) ) break;
+ }
+
+ if (itrace>1) {
+ *n2 = itrace/igy;
+ *n3 = igy;
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ *d2 = fabs(dxrcv)*scl;
+ *d3 = fabs(dyrcv)*scl;
+ if (NINT(dxrcv*1e3) != NINT(fabs(hdr.d2)*1e3)) {
+ if (dxrcv != 0) *d2 = fabs(dxrcv)*scl;
+ else *d2 = hdr.d2;
+ }
+ }
+ else {
+ *n2 = MAX(hdr.trwf, 1);
+ *n3 = 1;
+ *d2 = hdr.d2;
+ *d3 = 0.0;
+ dxrcv = hdr.d2;
+ dyrcv = 0.0;
+ }
+
+/* check if the total number of traces (ntraces) is correct */
+
+/* expensive way to find out how many gathers there are */
+
+// fprintf(stderr, "ngath = %d dxrcv=%f d2=%f scl=%f \n", *ngath, dxrcv, *d2, scl);
+ if (*ngath == 0) {
+ *n2 = 0;
+ *n3 = 0;
+
+ end_of_file = 0;
+ one_shot = 1;
+ igath = 0;
+ fseeko( fp, 0, SEEK_SET );
+ dxrcv = *d2;
+ dyrcv = *d3;
+
+ while (!end_of_file) {
+ itrace = 0;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { break; }
+ fldr_shot = hdr.fldr;
+ sx_shot = hdr.sx;
+ gx_start = hdr.gx;
+ gx_end = hdr.gx;
+ sy_shot = hdr.sy;
+ gy_start = hdr.gy;
+ gy_end = hdr.gy;
+
+ itrace = 0;
+ igy = 1;
+ while (one_shot) {
+ fseeko( fp, data_sz, SEEK_CUR );
+ itrace++;
+ if (hdr.gx != gx_end) dxrcv = MIN(dxrcv,abs(hdr.gx-gx_end));
+ if (hdr.gy != gy_end) {
+ igy++;
+ gy_end = hdr.gy;
+ dyrcv = MIN(dyrcv,abs(hdr.gy-gy_end));
+ }
+ gx_end = hdr.gx;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
+ }
+ if (itrace>1) {
+ *n2 = itrace/igy;
+ *n3 = igy;
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ dxsrc = (float)(hdr.sx - sx_shot)*scl;
+ dysrc = (float)(hdr.sy - sy_shot)*scl;
+ }
+ if (verbose>1) {
+ fprintf(stderr," . Scanning shot %d (%d) with %d traces dxrcv=%.2f dxsrc=%.2f %d %d dyrcv=%.2f dysrc=%.2f %d %d\n",sx_shot,igath,itrace,dxrcv*scl,dxsrc,gx_end,gx_start,dyrcv*scl,dysrc,gy_end,gy_start);
+ }
+ if (itrace != 0) { /* end of shot record */
+ fseeko( fp, -TRCBYTES, SEEK_CUR );
+ igath++;
+ }
+ else {
+ end_of_file = 1;
+ }
+ }
+ *ngath = igath;
+ *d2 = dxrcv*scl;
+ *d3 = dyrcv*scl;
+ }
+ else {
+ /* read last trace header */
+
+ fseeko( fp, -trace_sz, SEEK_END );
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ *ngath = ntraces/((*n2)*(*n3));
+ }
+// *nxm = NINT((*xmax-*xmin)/dxrcv)+1;
+ *nxm = (int)ntraces;
+
+ fclose( fp );
+ free(trace);
+
+ return 0;
+}
+
+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)
+{
+ vmess("file %s contains", file);
+ vmess("*** n1 = %d n2 = %d n3 = %d ntftt=%d", n1, n2, n3, optncr(n1));
+ vmess("*** d1 = %.5f d2 = %.5f d3 = %.5f", d1, d2, d3);
+ vmess("*** f1 = %.5f f2 = %.5f f3 = %.5f", f1, f2, f3);
+ vmess("*** fldr = %d sx = %d sy = %d", hdrs[0].fldr, hdrs[0].sx, hdrs[0].sy);
+
+ return 0;
+}
diff --git a/marchenko3D/readShotData3D.c b/marchenko3D/readShotData3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..db5ab8b10a07833d454febe754c29b252db5c00a
--- /dev/null
+++ b/marchenko3D/readShotData3D.c
@@ -0,0 +1,138 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+int optncr(int n);
+void cc1fft(complex *data, int n, int sign);
+void rc1fft(float *rdata, complex *cdata, int n, int sign);
+
+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)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ int fldr_shot, sx_shot, sy_shot, itrace, one_shot, igath, iw;
+ int end_of_file, nt, nxy;
+ int *isx, *igx, *isy, *igy, k, l, m, j;
+ int samercv, samesrc, nxrk, nxrm, maxtraces, ixsrc;
+ float scl, scel, *trace, dt;
+ complex *ctrace;
+
+ nxy = nx*ny;
+
+ /* Reading first header */
+
+ if (filename == NULL) fp = stdin;
+ else fp = fopen( filename, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", filename);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+ if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
+ else if (hdr.scalel == 0) scel = 1.0;
+ else scel = hdr.scalel;
+
+ fseek(fp, 0, SEEK_SET);
+
+ nt = hdr.ns;
+ dt = hdr.dt/(1E6);
+
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+ isx = (int *)malloc((nxy*nshots)*sizeof(int));
+ igx = (int *)malloc((nxy*nshots)*sizeof(int));
+ isy = (int *)malloc((nxy*nshots)*sizeof(int));
+ igy = (int *)malloc((nxy*nshots)*sizeof(int));
+
+
+ end_of_file = 0;
+ one_shot = 1;
+ igath = 0;
+
+ /* Read shots in file */
+
+ while (!end_of_file) {
+
+ /* start reading data (shot records) */
+ itrace = 0;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { /* no more data in file */
+ break;
+ }
+
+ sx_shot = hdr.sx;
+ sy_shot = hdr.sy;
+ fldr_shot = hdr.fldr;
+ isx[igath] = sx_shot;
+ isy[igath] = sy_shot;
+ xsrc[igath] = sx_shot*scl;
+ ysrc[igath] = sy_shot*scl;
+ zsrc[igath] = hdr.selev*scel;
+ xnx[igath]=0;
+ while (one_shot) {
+ igx[igath*nxy+itrace] = hdr.gx;
+ igy[igath*nxy+itrace] = hdr.gy;
+ xrcv[igath*nxy+itrace] = hdr.gx*scl;
+ yrcv[igath*nxy+itrace] = hdr.gy*scl;
+
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* transform to frequency domain */
+ if (ntfft > hdr.ns)
+ memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
+
+ rc1fft(trace,ctrace,ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ cdata[igath*nxy*nw+iw*nxy+itrace].r = scale*ctrace[nw_low+iw].r;
+ cdata[igath*nxy*nw+iw*nxy+itrace].i = scale*mode*ctrace[nw_low+iw].i;
+ }
+ itrace++;
+ xnx[igath]+=1;
+
+ /* read next hdr of next trace */
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
+ }
+ if (verbose>2) {
+ vmess("finished reading shot x=%d y=%d (%d) with %d traces",sx_shot,sy_shot,igath,itrace);
+ }
+
+ if (itrace != 0) { /* end of shot record */
+ fseek( fp, -TRCBYTES, SEEK_CUR );
+ igath++;
+ }
+ else {
+ end_of_file = 1;
+ }
+ }
+
+ free(ctrace);
+ free(trace);
+
+ return 0;
+}
diff --git a/marchenko3D/readTinvData3D.c b/marchenko3D/readTinvData3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..9c616218ecfedc88e08d57642e98fe48c4106758
--- /dev/null
+++ b/marchenko3D/readTinvData3D.c
@@ -0,0 +1,267 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+#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 findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute);
+
+int readTinvData(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)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ int fldr_shot, sx_shot, sy_shot, itrace, one_shot, ig, isyn, i, j;
+ int end_of_file, nt, gx0, gx1, gy0, gy1;
+ int nx1, ny1, jmax, imax, tstart, tend, nxy;
+ float xmax, tmax, lmax, ixmax, iymax;
+ float scl, scel, *trace, dxrcv;
+ complex *ctrace;
+
+ nxy = nx*ny;
+
+ /* Reading first header */
+
+ if (filename == NULL) fp = stdin;
+ else fp = fopen( filename, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", filename);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+ if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
+ else if (hdr.scalel == 0) scel = 1.0;
+ else scel = hdr.scalel;
+ fseek(fp, 0, SEEK_SET);
+
+ nt = hdr.ns;
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+
+ end_of_file = 0;
+ one_shot = 1;
+ isyn = 0;
+
+ /* Read shots in file */
+
+ while (!end_of_file) {
+
+ /* start reading data (shot records) */
+ itrace = 0;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { /* no more data in file */
+ break;
+ }
+
+ sx_shot = hdr.sx;
+ sy_shot = hdr.sy;
+ fldr_shot = hdr.fldr;
+ gx0 = hdr.gx;
+ gy0 = hdr.gy;
+ gy1 = gy0;
+ xsrc[isyn] = sx_shot*scl;
+ ysrc[isyn] = sy_shot*scl;
+ zsrc[isyn] = hdr.selev*scel;
+ xnx[isyn] = 0;
+ ig = isyn*nxy*ntfft;
+ ny1 = 1;
+ while (one_shot) {
+ xrcv[isyn*nxy+itrace] = hdr.gx*scl;
+ yrcv[isyn*nxy+itrace] = hdr.gy*scl;
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* copy trace to data array */
+ memcpy( &tinv[ig+itrace*ntfft], trace, nt*sizeof(float));
+
+ gx1 = hdr.gx;
+ if (gy1 != hdr.gy) {
+ gy1 = hdr.gy;
+ ny1++;
+ }
+ itrace++;
+
+ /* read next hdr of next trace */
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
+ }
+ if (verbose>2) {
+ fprintf(stderr,"finished reading shot x=%d y=%d (%d) with %d traces\n",sx_shot,isyn,itrace);
+ }
+
+ /* look for maximum in shot record to define mute window */
+ /* find consistent (one event) maximum related to maximum value */
+ nx1 = itrace/ny1;
+ xnx[isyn]=itrace;
+
+ /* alternative find maximum at source position */
+ dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
+ dyrcv = (gy1 - gy0)*scl/(float)(ny1-1);
+ //imax = NINT(((sx_shot-gx0)*scl)/dxrcv);
+ ixmax = NINT(((sx_shot-gx0)*scl)/dxrcv);
+ iymax = NINT(((sy_shot-gy0)*scl)/dyrcv);
+ tmax=0.0;
+ jmax = 0;
+ for (j = 0; j < nt; j++) {
+ lmax = fabs(tinv[ig+iymax*nx*ntfft+ixmax*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ xmax=lmax;
+ }
+ }
+ }
+ maxval[isyn*nxy+iymax*nx+ixmax] = jmax;
+ if (verbose >= 3) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, iymax, maxval[isyn*nxy+iymax*nx+ixmax]);
+
+ /* search forward in x-trace direction from maximum in file */
+ for (i = ixmax+1; i < nx1; i++) {
+ tstart = MAX(0, (maxval[isyn*nxy+iymax*nx+(i-1)]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nxy+iymax*nx+(i-1)]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+iymax*nx*ntfft+i*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nxy+iymax*nx+i] = jmax;
+ }
+ /* search backward in x-trace direction from maximum in file */
+ for (i = ixmax-1; i >=0; i--) {
+ tstart = MAX(0, (maxval[isyn*nxy+iymax*nx+i+1]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nxy+iymax*nx+i+1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+iymax*nx*ntfft+i*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nxy+iymax*nx+i] = jmax;
+ }
+
+ /* search forward in y-trace direction from maximum in file */
+ for (i = iymax+1; i < ny1; i++) {
+ tstart = MAX(0, (maxval[isyn*nxy+(i-1)*nx+ixmax]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nxy+(i-1)*nx+ixmax]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+i*nx*ntfft+ixmax*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nxy+i*nx+ixmax] = jmax;
+ }
+
+ /* search backward in y-trace direction from maximum in file */
+ for (i = iymax-1; i >= 0; i--) {
+ tstart = MAX(0, (maxval[isyn*nxy+(i+1)*nx+ixmax]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nxy+(i+1)*nx+ixmax]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+i*nx*ntfft+ixmax*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nxy+i*nx+ixmax] = jmax;
+ }
+
+ if (itrace != 0) { /* end of shot record, but not end-of-file */
+ fseek( fp, -TRCBYTES, SEEK_CUR );
+ isyn++;
+ }
+ else {
+ end_of_file = 1;
+ }
+
+ /* copy trace to data array for mode=-1 */
+ /* time reverse trace */
+ if (mode==-1) {
+ for (i = 0; i < nx1; i++) {
+ memcpy( trace, &tinv[ig+i*ntfft], ntfft*sizeof(float));
+ j=0;
+ tinv[ig+i*ntfft+j] = trace[j];
+ for (j=1; j<ntfft; j++) tinv[ig+i*ntfft+ntfft-j] = trace[j];
+ }
+ }
+ }
+
+ free(ctrace);
+ free(trace);
+
+ return 0;
+}
+
+
+/* simple sort algorithm */
+void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute)
+{
+ int i, sign;
+ float diff1, diff2;
+
+ *imute=0;
+
+ if (xrcvMute[0] < xrcvMute[1]) sign = 1;
+ else sign = -1;
+
+ if (sign == 1) {
+ i = 0;
+ while (xrcvMute[i] < xrcvShot && i < nxs) {
+ i++;
+ }
+ /* i is now position larger than xrcvShot */
+ }
+ else {
+ i = 0;
+ while (xrcvMute[i] > xrcvShot && i < nxs) {
+ i++;
+ }
+ /* i is now position smaller than xrcvShot */
+ }
+
+ diff1 = fabsf(xrcvMute[i]-xrcvShot);
+ diff2 = fabsf(xrcvMute[i-1]-xrcvShot);
+ if (diff1 < diff2) *imute = i;
+ else *imute = i-1;
+
+ return;
+}
+
diff --git a/marchenko3D/test.c b/marchenko3D/test.c
new file mode 100755
index 0000000000000000000000000000000000000000..157f20d8946be16d021401d68376672f4583d265
--- /dev/null
+++ b/marchenko3D/test.c
@@ -0,0 +1,84 @@
+#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 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);
+
+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;
+ float dt, dx, dy, ft, fx, fy, scl, fmin, fmax;
+ float *xrcv, *yrcv, *xsrc, *ysrc, *zsrc;
+ complex *Refl;
+ int nt, nx, ny, nshots, ntraces, ret, *xnx;
+ int ntfft, nfreq, nw_low, nw_high, nw, mode, verbose;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
+ if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (!getparfloat("fmax", &fmax)) fmax = 70.0;
+ if (file_shot == NULL) verr("Give me a file name for the love of God");
+
+ nshots=0;
+ mode=1;
+ verbose=8;
+
+ ret = getFileInfo3D(file_shot, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+ vmess("n1=%d n2=%d n3=%d ngath=%d ntraces=%d",nt,nx,ny,nshots,ntraces);
+ vmess("d1=%.5f d2=%.5f d3=%.5f f1=%.5f f2=%.5f f3=%.5f scl=%.5f",dt,dx,dy,ft,fx,fy,scl);
+
+ ntfft = optncr(nt);
+ nfreq = ntfft/2+1;
+ nw_low = (int)MIN((fmin*ntfft*dt), nfreq-1);
+ nw_low = MAX(nw_low, 1);
+ nw_high = MIN((int)(fmax*ntfft*dt), nfreq-1);
+ nw = nw_high - nw_low + 1;
+ scl = 1.0/((float)ntfft);
+
+ Refl = (complex *)malloc(nw*nx*ny*nshots*sizeof(complex)); // reflection response in frequency domain
+ xrcv = (float *)calloc(nshots*nx*ny,sizeof(float)); // x-rcv postions of shots
+ xsrc = (float *)calloc(nshots,sizeof(float)); // x-src position of shots
+ yrcv = (float *)calloc(nshots*nx*ny,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, nx, ny, ntfft, mode, scl, verbose);
+
+ free(Refl);free(xrcv);free(yrcv);free(xsrc);free(ysrc);free(zsrc);free(xnx);
+
+ return 0;
+}