3D

marchenko3D/ampest3D_backup.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+#include "par.h"
+#include <genfft.h>
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+    float r,i;
+} complex;
+#endif/* complex */
+
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+long loptncr(long n);
+long maxest3D(float *data, long nt);
+long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
+void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
+void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift);
+void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
+
+void AmpEst3D(float *f1d, float *Gd, float *ampest, long Nfoc, long nxs, long nys, long ntfft, long *ixpos, long npos,
+    char *file_wav, float dx, float dy, float dt)
+{
+	
+	long 	l, i, ix, iw, nfreq;
+	float 	scl, sclt, *wavelet, *scaled, *conv, *f1dsamp;
+	float   dtm, dxm, cpm, rom, *trace;
+	FILE 	*fp_wav;
+	segy 	*hdrs_wav;
+
+	scl = dx*dy;
+    sclt = 1.0*dt/((float)ntfft);
+
+	conv	= (float *)calloc(nys*nxs*ntfft,sizeof(float));
+	wavelet	= (float *)calloc(ntfft,sizeof(float));
+	scaled	= (float *)calloc(ntfft,sizeof(float));
+	f1dsamp	= (float *)calloc(nys*nxs*ntfft,sizeof(float));
+
+	if (file_wav!=NULL) {
+		trace	= (float *)calloc(ntfft,sizeof(float));
+		hdrs_wav = (segy *)calloc(1, sizeof(segy));
+    	fp_wav = fopen(file_wav, "r");
+        if (fp_wav==NULL) verr("error on opening wavelet file %s", file_wav);
+    	readData3D(fp_wav, trace, hdrs_wav, 0);
+    	fclose(fp_wav);
+		corr(trace, trace, wavelet,  1, ntfft, dt, 0);
+		free(hdrs_wav); free(trace);
+		/* For a monopole source the scaling is (2.0*dt*cp*cp*rho)/(dx*dx) */
+		for (iw=0; iw<ntfft; iw++){
+			wavelet[iw] *= dt;
+		}
+	}
+
+	for (l=0; l<Nfoc; l++) {
+		for (i=0; i<npos; i++) {
+			ix = ixpos[i];
+			for (iw=0; iw<ntfft; iw++) {
+				f1dsamp[i*ntfft+iw] = f1d[l*nxs*nys*ntfft+ix*ntfft+iw];
+			}
+		}
+		if (file_wav==NULL){
+			corr(f1dsamp, f1dsamp, conv,  nxs*nys, ntfft, dt, 0);
+			for (i=0; i<nxs*nys; i++) {
+				for (iw=0; iw<ntfft; iw++) {
+					wavelet[iw] += dt*scl*conv[i*ntfft+iw];
+				}
+			}
+		}
+		memset(&conv[0],0.0, sizeof(float)*ntfft*nxs*nys);
+		convol(f1dsamp, &Gd[l*nxs*nys*ntfft], conv, nxs*nys, ntfft, dt, 0);
+		for (i=0; i<nxs*nys; i++) {
+			for (iw=0; iw<ntfft; iw++) {
+				scaled[iw] += dt*scl*conv[i*ntfft+iw];
+			}
+		}
+		ampest[l] = sqrtf(wavelet[0]/scaled[0]);
+		memset(&conv[0],0.0,    sizeof(float)*ntfft*nxs*nys);
+		memset(&scaled[0],0.0,  sizeof(float)*ntfft);
+	}
+	free(wavelet);free(scaled);free(conv);free(f1dsamp);
+
+	return;
+}
+
+long maxest3D(float *data, long nt)
+{
+	float maxt;
+	long it;
+
+	maxt = data[0];
+	for (it = 0; it < nt; it++) {
+		if (fabs(data[it]) > fabs(maxt)) maxt=data[it];
+	}
+
+	return maxt;
+}
+
+/**
+* Calculates the time convolution of two arrays by 
+* transforming the arrayis to frequency domain,
+* multiplies the arrays and transform back to time.
+*
+**/
+
+void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift)
+{
+	long 	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 = loptncr(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], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)sign);
+	rcmfft(&rdata2[0], &cdata2[0], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)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 = 1.0/(2.0*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], (int)optn, (int)nrec, (int)nfreq, (int)optn, (int)sign);
+	scl_data(rdata1,optn,nrec,scl,con,nsam);
+
+	free(ccon);
+	free(rdata1);
+	free(rdata2);
+	return;
+}
+
+/**
+* Calculates the time correlation of two arrays by 
+* transforming the arrayis to frequency domain,
+* multiply the arrays and transform back to time.
+*
+**/
+
+
+void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift)
+{
+	long 	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 = loptncr(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], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)sign);
+	rcmfft(&rdata2[0], &cdata2[0], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)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], (int)optn, (int)nrec, (int)nfreq, (int)optn, (int)sign);
+	scl_data(rdata1,optn,nrec,scl,cov,nsam);
+
+	free(ccov);
+	free(rdata1);
+	free(rdata2);
+	return;
+}
+
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout)
+{
+	long 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, long nsam, long nrec, float scl, float *datout, long nsamout)
+{
+	long it,ix;
+	for (ix = 0; ix < nrec; ix++) {
+		for (it = 0 ; it < nsamout ; it++)
+			datout[ix*nsamout+it] = scl*data[ix*nsam+it];
+	}
+}
\ No newline at end of file

marchenko3D/imaging3D.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>
+
+double wallclock_time(void);
+
+void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
+void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
+
+
+void imaging3D(float *Image, float *Gmin, float *f1plus, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose)
+{
+    long     i, l, count=0;
+	float   *conv;
+    double  t0, t2;
+
+    t0   = wallclock_time();
+
+#pragma omp parallel default(shared) \
+  private(i,conv) 
+{	
+	conv	= (float *)calloc(nx*ny*nt,sizeof(float));
+
+#pragma omp for
+	for (l = 0; l < Nfoc; l++) {
+		count+=1;
+		if (verbose > 2) vmess("Imaging location %d out of %d",count,Nfoc);
+
+		convol(&Gmin[l*nx*ny*nt], &f1plus[l*nx*ny*nt], conv, nx*ny, nt, dt, 0);
+		for (i=0; i<nx*ny; i++) {
+        	Image[l] += conv[i*nt]*dx*dy*dt;
+		}
+	}
+    free(conv);
+}
+
+    t2 = wallclock_time();
+    if (verbose) {
+        vmess("Total Imaging time = %.3f", t2-t0);
+    }
+
+    return;
+}
\ No newline at end of file

marchenko3D/makeWindow3D.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+    float r,i;
+} complex;
+#endif/* 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) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+void findShotInMute(float *xrcvMute, float xrcvShot, long nxs, long *imute);
+long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
+long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps, long nx, long ny, long nz, long sx, long ex, long sy, long ey, long sz, long ez);
+
+
+
+void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, 
+    long *xnx, long Nfoc, long nx, long ny, long ntfft, long *maxval, float *tinv, long verbose)
+{
+	FILE *fp;
+	segy hdr, *hdrs_mute, *hdrs_amp;
+	size_t nread;
+	long ig, is, i, iy, ix, j, l, nfreq, ntwav;
+	float *wavelet, *wavtmp, scl, *timeval, dw, *amp;
+	complex *cmute, *cwav;
+
+	/*Define parameters*/
+	nfreq   = ntfft/2+1;	
+	wavelet = (float *)calloc(ntfft,sizeof(float));
+	cwav	= (complex *)malloc(nfreq*sizeof(complex));
+	cmute	= (complex *)malloc(nfreq*sizeof(complex));
+	dw  	= 2*M_PI/(ntfft*dt);
+
+	/*Create wavelet using parameters or read in wavelet*/
+    if (file_wav != NULL) {
+        //Determine the amount of sample
+        if (verbose>0) vmess("Reading in wavelet");
+        fp = fopen( file_wav, "r" );
+        if ( fp == NULL ) {
+            perror("Error opening file containing wavelet");
+        }
+        nread = fread( &hdr, 1, TRCBYTES, fp );
+        ntwav = hdr.ns;
+        fclose(fp);
+        //Read in the wavelet
+        fp = fopen( file_wav, "r" );
+	    wavtmp = (float *)calloc(ntwav,sizeof(float));
+        readData3D(fp, wavtmp, &hdr, ntwav);
+        //Fit the wavelet into the same time-axis as the Marchenko scheme
+        for (i=0; i<ntwav; i++) {
+            wavelet[i] = wavtmp[i];
+            wavelet[ntfft-1-i] = wavtmp[ntwav-1-i];
+        }
+        rc1fft(wavelet,cwav,ntfft,-1);
+        free(wavtmp);
+        free(wavelet);
+    }
+
+
+	timeval = (float *)calloc(Nfoc*nx*ny,sizeof(float));
+	if (file_amp!=NULL) amp = (float *)calloc(Nfoc*nx*ny,sizeof(float));
+
+
+    /* Defining mute window using raytimes */
+    vmess("Using raytime for mutewindow");
+    hdrs_mute = (segy *) calloc(Nfoc,sizeof(segy));
+    fp = fopen( file_ray, "r" );
+    if ( fp == NULL ) {
+        perror("Error opening file containing ray");
+    }
+    fclose(fp);
+    readSnapData3D(file_ray, timeval, hdrs_mute, Nfoc, 1, 1, nx, 0, 1, 0, 1, 0, nx);
+    // for (is=0; is<Nfoc; is++) {
+    //     readData3D(fp, &timeval[is*nx*ny], &hdrs_mute[is], nx);
+    // }
+
+    /*Check whether the amplitude is also used*/
+    if (file_amp != NULL) {
+        vmess("Using ray-amplitudes");
+        hdrs_amp = (segy *) calloc(Nfoc,sizeof(segy));
+        fp = fopen( file_amp, "r" );
+        if ( fp == NULL ) {
+            perror("Error opening file containing ray-amplitude");
+        }
+        fclose(fp);
+        readSnapData3D(file_amp, amp, hdrs_amp, Nfoc, 1, 1, nx, 0, 1, 0, 1, 0, nx);
+        // for (is=0; is<Nfoc; is++) {
+        //     readData3D(fp, &amp[is*nx*ny], &hdrs_amp[is], nx);
+        // }
+    }
+
+    /*Define source and receiver locations from the raytime*/
+    for (is=0; is<Nfoc; is++) {
+        for (iy=0; iy<ny; iy++) {
+            for (ix=0; ix<nx; ix++) {
+                xrcv[is*nx*ny+iy*nx+ix] = (hdrs_mute[is].f1 + hdrs_mute[is].d1*((float)ix));
+                yrcv[is*nx*ny+iy*nx+ix] = hdrs_mute[is].gy;
+            }
+        }
+        xnx[is]=hdrs_mute[is].ns;
+        if (hdrs_mute[is].scalco < 0) scl=-1.0/hdrs_mute[is].scalco;
+        else scl=hdrs_mute[is].scalco;
+        xsrc[is] = hdrs_mute[is].sx*scl;
+        ysrc[is] = hdrs_mute[is].sy*scl;
+        zsrc[is] = hdrs_mute[is].sdepth*scl;
+    }
+
+
+	/*Determine the mutewindow*/
+	for (j=0; j<Nfoc; j++) {
+        for (l=0; l<ny; l++) {
+            for (i=0; i<nx; i++) {
+                maxval[j*ny*nx+l*nx+i] = (long)roundf(timeval[j*ny*nx+l*nx+i]/dt);
+                if (maxval[j*ny*nx+l*nx+i] > ntfft-1) maxval[j*ny*nx+l*nx+i] = ntfft-1;
+                if (file_wav!=NULL) { /*Apply the wavelet to create a first arrival*/
+                    if (file_amp != NULL) {
+                        for (ig=0; ig<nfreq; ig++) {
+                            cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]*amp[j*ny*nx+l*nx+i]);
+                            cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]*amp[j*ny*nx+l*nx+i]);
+                        }
+                    }
+                    else { /*Use the raytime only to determine the mutewindow*/
+                        for (ig=0; ig<nfreq; ig++) {
+                            cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0));
+                            cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0));
+                        }
+                    }
+                }
+                else {
+                    for (ig=0; ig<nfreq; ig++) {
+                        cmute[ig].r = (1.0/sqrtf((float)ntfft))*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0);
+                        cmute[ig].i = (1.0/sqrtf((float)ntfft))*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0);
+                    }
+                }
+                cr1fft(cmute,&tinv[j*ny*nx*ntfft+l*nx*ntfft+i*ntfft],ntfft,1);
+            }
+        }
+    }
+
+
+	return;
+}
+

marchenko3D/readSnapData3D.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;
+
+long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps, long nx, long ny, long nz, long sx, long ex, long sy, long ey, long sz, long ez)
+{
+	FILE *fp;
+	segy hdr;
+	size_t nread;
+	long nt, it, ix, iy, iz, dx, dy, dz;
+	float *tmpdata;
+
+	tmpdata = (float *)malloc(nsnaps*nx*ny*nz*sizeof(float));
+	/* 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;
+	}
+    //nread = fread(&hdr, 1, TRCBYTES, fp);
+    for (it = 0; it < nsnaps*nx*ny; it++) {
+		nread = fread(&hdr, 1, TRCBYTES, fp);
+		if (nread != TRCBYTES) {
+			break;
+		}
+		assert(nread == TRCBYTES);
+        nread = fread(&tmpdata[it*nz], sizeof(float), nz, fp);
+        assert (nread == nz);
+		memcpy(&hdrs[it], &hdr, TRCBYTES);
+    }
+	dx = ex-sx;
+	dy = ey-sy;
+	dz = ez-sz;
+	for (iz = sz; iz < ez; iz++) {
+        for (iy = sy; iy < ey; iy++) {
+            for (ix = sx; ix < ex; ix++) {
+                for (it = 0; it < nsnaps; it++) {
+                    data[it*dy*dx*dz+(iy-sy)*dx*dz+(ix-sx)*dz+iz-sz]=tmpdata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
+                }
+            }
+        }
+    }
+	fclose(fp);
+	free(tmpdata);
+	return 0;
+}

utils/ampdet.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 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 readData(FILE *fp, float *data, segy *hdrs, int n1);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+void complex_sqrt(complex *z);
+void deconv_small(complex *c1, complex *c2, complex *c3, float nkx, float nfreq, float reps, float eps);
+void conv_small(complex *c1, complex *c2, complex *c3, float nkx, float nfreq);
+void corr_small(complex *c1, complex *c2, complex *c3, float nkx, float nfreq);
+void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout);
+void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout);
+void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout);
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" ampdet - Determine amplitude",
+" ",
+" author  : Jan Thorbecke : 19-04-1995 (janth@xs4all.nl)",
+" product : Originates from DELPHI software",
+"                         : revision 2010",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+    FILE    *fp;
+	char    *file_gp, *file_fp, *file_wav;
+    int     nx, nt, ngath, ntraces, ret, size, nxwav;
+    int     ntfft, nfreq, nxfft, nkx, i, j, n;
+    float   dx, dt, fx, ft, xmin, xmax, scl;
+    float   df, dw, dkx, eps, reps;
+    float   *Gpd, *f1pd, *G_pad, *f_pad, *wav, *wav_pad;
+    complex *G_w, *f_w, *Gf, *amp, *wav_w, *S, *ZS, *SS;
+    segy    *hdr_gp, *hdr_fp, *hdr_wav;
+
+	initargs(argc, argv);
+	requestdoc(1);
+
+	if(!getparstring("file_gp", &file_gp)) file_gp=NULL;
+    if (file_gp==NULL) verr("file %s does not exist",file_gp);
+    if(!getparstring("file_gp", &file_fp)) file_fp=NULL;
+    if (file_fp==NULL) verr("file %s does not exist",file_fp);
+    if(!getparstring("file_wav", &file_wav)) file_wav=NULL;
+    if (file_wav==NULL) verr("file %s does not exist",file_wav);
+	if(!getparfloat("eps", &eps)) eps=0.00;
+	if(!getparfloat("reps", &reps)) reps=0.01;
+
+    ngath = 1;
+    ret = getFileInfo(file_gp, &nt, &nx, &ngath, &dt, &dx, &ft, &fx, &xmin, &xmax, &scl, &ntraces);
+
+    size    = nt*nx;
+
+	Gpd     = (float *)malloc(size*sizeof(float));
+	hdr_gp  = (segy *) calloc(nx,sizeof(segy));
+    fp      = fopen(file_gp, "r");
+	if (fp == NULL) verr("error on opening input file_in1=%s", file_gp);
+    nx      = readData(fp, Gpd, hdr_gp, nt);
+    fclose(fp);
+
+	f1pd    = (float *)malloc(size*sizeof(float));
+	hdr_fp  = (segy *) calloc(nx,sizeof(segy));
+    fp      = fopen(file_fp, "r");
+	if (fp == NULL) verr("error on opening input file_in1=%s", file_fp);
+    nx      = readData(fp, f1pd, hdr_fp, nt);
+    fclose(fp);
+
+    wav     = (float *)malloc(nt*sizeof(float));
+	hdr_wav = (segy *) calloc(1,sizeof(segy));
+    fp      = fopen(file_wav, "r");
+	if (fp == NULL) verr("error on opening input file_in1=%s", file_fp);
+    nxwav   = readData(fp, wav, hdr_wav, nt);
+    fclose(fp);
+    vmess("test:%d",nxwav);
+
+    /* Start the scaling */
+    ntfft   = optncr(nt);
+	nfreq   = ntfft/2+1;
+	df      = 1.0/(ntfft*dt);
+    dw      = 2.0*PI*df;
+	nkx     = optncc(nx);
+	dkx     = 2.0*PI/(nkx*dx);
+
+    vmess("ntfft:%d, nfreq:%d, nkx:%d",ntfft,nfreq,nkx);
+
+    /* Allocate the arrays */
+    G_pad = (float *)malloc(ntfft*nkx*sizeof(float));
+	if (G_pad == NULL) verr("memory allocation error for G_pad");
+    f_pad = (float *)malloc(ntfft*nkx*sizeof(float));
+	if (f_pad == NULL) verr("memory allocation error for f_pad");
+    wav_pad = (float *)malloc(ntfft*sizeof(float));
+	if (wav_pad == NULL) verr("memory allocation error for wav_pad");
+    G_w   = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (G_w == NULL) verr("memory allocation error for G_w");
+    f_w   = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (f_w == NULL) verr("memory allocation error for f_w");
+    Gf    = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (Gf == NULL) verr("memory allocation error for Gf");
+    wav_w = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (wav_w == NULL) verr("memory allocation error for wav_w");
+    amp   = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (amp == NULL) verr("memory allocation error for amp");
+    S   = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (S == NULL) verr("memory allocation error for S");
+    ZS   = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (ZS == NULL) verr("memory allocation error for ZS");
+    SS   = (complex *)malloc(nfreq*nkx*sizeof(complex));
+	if (SS == NULL) verr("memory allocation error for SS");
+
+    /* pad zeroes in 2 directions to reach FFT lengths */
+	pad2d_data(Gpd, nt,nx,ntfft,nkx,G_pad);
+	pad2d_data(f1pd,nt,nx,ntfft,nkx,f_pad);
+    pad_data(wav, nt, 1, ntfft, wav_pad);
+
+    /* double forward FFT */
+	xt2wkx(&G_pad[0], &G_w[0], ntfft, nkx, ntfft, nkx, 0);
+	xt2wkx(&f_pad[0], &f_w[0], ntfft, nkx, ntfft, nkx, 0);
+    rcmfft(&wav_pad[0], &wav_w[0], ntfft, 1, ntfft, nfreq, -1);
+
+    for (i=1; i<nkx; i++) {
+        for (j=0; j<nfreq; j++) {
+            wav_w[i*nfreq+j] = wav_w[j];
+        }
+    }
+
+    /* Create Z*(|S|*)/(|S|*(|S|*)) */
+    conv_small(  G_w,   f_w,   Gf,  nkx, nfreq); // Z
+    corr_small(  wav_w, wav_w, S,   nkx, nfreq); //|S|
+    corr_small(  Gf,    G_w,   ZS,  nkx, nfreq); // Z *(|S|*)
+    corr_small(  G_w,   G_w,   SS,  nkx, nfreq); //|S|*(|S|*)
+    deconv_small(ZS,    SS,    amp, nkx, nfreq, reps, eps); // amp
+
+    for (i=0; i<nkx*nfreq; i++) {
+        complex_sqrt(&amp[i]);
+    }
+    
+    conv_small(G_w, amp, Gf, nkx, nfreq); // Scaled data
+
+    /* inverse double FFT */
+	wkx2xt(&Gf[0], &G_pad[0], ntfft, nkx, nkx, ntfft, 0);
+	/* select original samples and traces */
+	scl = 1.0;
+	scl_data(G_pad,ntfft,nx,scl,Gpd ,nt);
+
+    fp      = fopen("out.su", "w+");
+    ret = writeData(fp, Gpd, hdr_gp, nt, nx);
+	if (ret < 0 ) verr("error on writing output file.");
+    fclose(fp);
+
+    free(f1pd);free(Gpd);free(hdr_gp);free(hdr_fp);
+
+	return 0;
+}
+
+void conv_small(complex *c1, complex *c2, complex *c3, float nkx, float nfreq)
+{
+
+    float   *qr, *qi, *p1r, *p1i, *p2r, *p2i;
+    int     n, j;
+
+    /* apply convolution */
+	p1r = (float *) &c1[0];
+	p2r = (float *) &c2[0];
+	qr = (float *) &c3[0].r;
+	p1i = p1r + 1;
+	p2i = p2r + 1;
+	qi = qr + 1;
+	n = nkx*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;
+	}
+}
+
+void corr_small(complex *c1, complex *c2, complex *c3, float nkx, float nfreq)
+{
+
+    float   *qr, *qi, *p1r, *p1i, *p2r, *p2i;
+    int     n, j;
+
+    /* apply convolution */
+	p1r = (float *) &c1[0];
+	p2r = (float *) &c2[0];
+	qr = (float *) &c3[0].r;
+	p1i = p1r + 1;
+	p2i = p2r + 1;
+	qi = qr + 1;
+	n = nkx*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;
+	}
+}
+
+void deconv_small(complex *c1, complex *c2, complex *c3, float nkx, float nfreq, float reps, float eps)
+{
+
+    float   *qr, *qi, *p1r, *p1i, *p2r, *p2i, maxden, *den, leps;
+    int     n, j;
+
+    den = (float *)malloc(nfreq*nkx*sizeof(float));
+	if (den == NULL) verr("memory allocation error for den");
+
+    /* apply deconvolution */
+	p1r = (float *) &c1[0];
+	p2r = (float *) &c2[0];
+	p1i = p1r + 1;
+	p2i = p2r + 1;
+	n = nkx*nfreq;
+	maxden=0.0;
+	for (j = 0; j < n; j++) {
+		den[j] = *p2r**p2r + *p2i**p2i;
+		maxden = MAX(den[j], maxden);
+		p2r += 2;
+		p2i += 2;
+	}
+	p1r = (float *) &c1[0];
+	p2r = (float *) &c2[0];
+	qr = (float *) &c3[0].r;
+	p1i = p1r + 1;
+	p2i = p2r + 1;
+    qi = qr + 1;
+	leps = reps*maxden+eps;
+	for (j = 0; j < n; j++) {
+
+		if (fabs(*p2r)>=fabs(*p2i)) {
+			*qr = (*p2r**p1r+*p2i**p1i)/(den[j]+leps);
+			*qi = (*p2r**p1i-*p2i**p1r)/(den[j]+leps);
+		} else {
+			*qr = (*p1r**p2r+*p1i**p2i)/(den[j]+leps);
+			*qi = (*p1i**p2r-*p1r**p2i)/(den[j]+leps);
+		}
+		qr += 2;
+		qi += 2;
+		p1r += 2;
+		p1i += 2;
+		p2r += 2;
+		p2i += 2;
+	}
+}
+
+void complex_sqrt(complex *z)
+{
+    float zmod, zmodr, zzmr, zzmi, zzm;
+
+    zmod  = sqrtf(z[0].r*z[0].r+z[0].i*z[0].i);
+    zmodr = sqrtf(zmod);
+    zzmr  = z[0].r + zmod;
+    zzmi  = z[0].i;
+    zzm   = sqrtf(zzmr*zzmr+zzmi*zzmi);
+
+    z[0].r = (zmodr*zzmr)/zzm;
+    z[0].i = (zmodr*zzmi)/zzm;
+}
+
+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 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 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];
+	}
+}
\ No newline at end of file