ogun

marchenko3D/demo/test/test.scr

+#!/bin/bash
+
+./../../test file_shot=shot_gy.su
+./../../test file_shot=shot_sxgy.su
+./../../test file_shot=shot_sxsygy.su

marchenko3D/getFileInfo3D.c

+#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;
+}

marchenko3D/readShotData3D.c

+#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;
+}

marchenko3D/readTinvData3D.c

+#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;
+}
+

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 */
+
+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;
+}