diff --git a/utils/pwshift.c b/utils/pwshift.c
new file mode 100644
index 0000000000000000000000000000000000000000..54999a06efd7135edefe7b638bf14c03d58de2d5
--- /dev/null
+++ b/utils/pwshift.c
@@ -0,0 +1,374 @@
+#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) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3,
+ float *f1, float *f2, float *f3, float *sclsxgxsygy, long *nxm);
+double wallclock_time(void);
+long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
+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 timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax);
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
+void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
+
+char *sdoc[] = {
+" ",
+" combine - Combine results into a single result ",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_gmin= ................. File containing the G- data",
+" file_ray= .................. File containing the ray time data",
+" file_amp= .................. File containing the ray amplitude data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out=out.su .......... Filename of the output",
+" verbose=1 ................ Give detailed information of process",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_gmin, *fp_ray, *fp_amp, *fp_out;
+ char *file_gmin, *file_ray, *file_amp, *file_out, fbr[100], fer[100], fba[100], fea[100], fins[100], fin2[100], numb1[100], *ptr;
+ float *gmin, *conv, *time, *amp, *image, fmin, fmax, *delay;
+ float dt, dy, dx, dz, t0, y0, x0, scl;
+ float dt_ray, dy_ray, dx_ray, t0_ray, y0_ray, x0_ray, scl_ray, px, py, src_velox, src_veloy, src_anglex, src_angley, grad2rad;
+ long nshots, nt, ny, nx, ntr;
+ long nray, nt_ray, ny_ray, nx_ray, ntr_ray;
+ long verbose, ix, iy, it, iz, is, *gx, *gy, *gz, numb, dnumb, pos, nzs, file_det;
+ size_t ret;
+ segy *hdr_gmin, *hdr_time, *hdr_amp, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
+ if (!getparstring("file_gmin", &file_gmin)) file_gmin = NULL;
+ if (!getparstring("file_ray", &file_ray)) file_ray = NULL;
+ if (!getparstring("file_amp", &file_amp)) file_amp = NULL;
+ if (!getparstring("file_out", &file_out)) file_out = "out.su";
+ if (!getparlong("verbose", &verbose)) verbose=1;
+ if (!getparfloat("fmin", &fmin)) fmin=0.0;
+ if (!getparfloat("fmax", &fmax)) fmax=70.0;
+ if (!getparfloat("src_velox", &src_velox)) src_velox=1500.0;
+ if (!getparfloat("src_veloy", &src_veloy)) src_veloy=1500.0;
+ if (!getparfloat("src_anglex", &src_anglex)) src_anglex=0.0;
+ if (!getparfloat("src_angley", &src_angley)) src_angley=0.0;
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=1;
+ if (file_gmin == NULL) verr("Incorrect gmin input");
+ if (file_ray == NULL) verr("Incorrect ray time input");
+ if (file_amp == NULL) verr("Incorrect ray amplitude input");
+
+ /*----------------------------------------------------------------------------*
+ * Determine the position of the number in the string
+ * and split the file into beginning, middle and end
+ *----------------------------------------------------------------------------*/
+ if (dnumb < 1) dnumb = 1;
+ sprintf(numb1,"z%li",numb);
+
+ ptr = strstr(file_ray,numb1);
+ pos = ptr - file_ray + 1;
+ sprintf(fbr,"%*.*s", pos-1, pos-1, file_ray);
+ sprintf(fer,"%s", file_ray+pos+1);
+
+ ptr = strstr(file_amp,numb1);
+ pos = ptr - file_amp + 1;
+ sprintf(fba,"%*.*s", pos-1, pos-1, file_ray);
+ sprintf(fea,"%s", file_ray+pos+1);
+
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files that are present
+ *----------------------------------------------------------------------------*/
+ file_det = 1;
+ nzs=0;
+ while (file_det) { // Check for a file with the filename
+ sprintf(fins,"z%li",nzs*dnumb+numb);
+ sprintf(file_ray,"%s%s%s",fbr,fins,fer);
+ fp_ray = fopen(file_ray, "r");
+ if (fp_ray == NULL) { // If the filename does not exist
+ if (nzs == 0) { // The filename is wrong to begin with
+ verr("error on opening basefile=%s", file_ray);
+ }
+ else if (nzs == 1) { // There is only a single file
+ vmess("1 file detected");
+ file_det = 0;
+ break;
+ }
+ else { // Stop after the final file has been detected
+ vmess("%li files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_ray);
+ nzs++;
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Read in the first two files and determine the header values
+ * of the output
+ *----------------------------------------------------------------------------*/
+ getFileInfo3D(file_gmin, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &scl, &ntr);
+
+ if (verbose) {
+ vmess("************************ Gmin info ************************");
+ vmess("Number of depth levels : %li",nshots);
+ vmess("Number of samples x: %li, y: %li, t: %li",nx,ny,nt);
+ vmess("Starting distance for x: %.3f, y: %.3f, t: %.3f",x0,y0,t0);
+ vmess("Sampling distance for x: %.3f, y: %.3f, t: %.3f",dx,dy,dt);
+ vmess("***********************************************************");
+ }
+
+ sprintf(fins,"z%li",0);
+ sprintf(file_ray,"%s%s%s",fbr,fins,fer);
+ getFileInfo3D(file_ray, &nx_ray, &nt_ray, &ny_ray, &nray, &dx_ray, &dt_ray, &dy_ray, &x0_ray, &t0_ray, &y0_ray, &scl_ray, &ntr_ray);
+
+ if (verbose) {
+ vmess("************************ ray info *************************");
+ vmess("Number of depth levels : %li",nzs);
+ vmess("Number of focal points : %li",nray);
+ vmess("Number of samples x: %li, y: %li, t: %li",nx_ray,ny_ray,nt_ray);
+ vmess("Starting distance for x: %.3f, y: %.3f, t: %.3f",x0_ray,y0_ray,t0_ray);
+ vmess("Sampling distance for x: %.3f, y: %.3f, t: %.3f",dx_ray,dy_ray,dt_ray);
+ vmess("***********************************************************");
+ }
+
+ if (nshots!=nzs) verr("The depth levels of the rays (%li) do not match those of the plane waves (%li)",nzs,nshots);
+
+ /*----------------------------------------------------------------------------*
+ * Read in a single file to determine if the header values match
+ * and allocate the data
+ *----------------------------------------------------------------------------*/
+ hdr_gmin = (segy *)calloc(nshots*nx*ny,sizeof(segy));
+ gmin = (float *)calloc(nshots*nx*ny*nt,sizeof(float));
+ hdr_out = (segy *)calloc(nray*nshots,sizeof(segy));
+ delay = (float *)calloc(nx*ny,sizeof(float));
+
+ image = (float *)calloc(nray*nshots,sizeof(float));
+ gx = (long *)calloc(nray*nshots,sizeof(long));
+ gy = (long *)calloc(nray*nshots,sizeof(long));
+ gz = (long *)calloc(nray*nshots,sizeof(long));
+
+ readSnapData3D(file_gmin, gmin, hdr_gmin, nshots, nx, ny, nt, 0, nx, 0, ny, 0, nt);
+
+ /*----------------------------------------------------------------------------*
+ * Add the delay in case the plane wave is at an angle
+ *----------------------------------------------------------------------------*/
+
+ grad2rad = 17.453292e-3;
+ px = sin(src_anglex*grad2rad)/src_velox;
+ py = sin(src_angley*grad2rad)/src_veloy;
+
+ if (py < 0.0) {
+ for (iy=0; iy<ny; iy++) {
+ if (px < 0.0) {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + fabsf((ny-1-iy)*dy*py);
+ }
+ }
+ else {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = ix*dx*px + fabsf((ny-1-iy)*dy*py);
+ }
+ }
+ }
+ }
+ else {
+ for (iy=0; iy<ny; iy++) {
+ if (px < 0.0) {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + iy*dy*py;
+ }
+ }
+ else {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = ix*dx*px + iy*dy*py;
+ }
+ }
+ }
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Apply the imaging condition
+ *----------------------------------------------------------------------------*/
+#pragma omp parallel for schedule(static,1) default(shared) \
+ private(is,it,ix,fins,fin2,time,amp,hdr_time,hdr_amp,conv)
+ for (iy = 0; iy < nshots; iy++) {
+
+ hdr_time = (segy *)calloc(ny*nray,sizeof(segy));
+ time = (float *)calloc(nx*ny*nray,sizeof(float));
+ hdr_amp = (segy *)calloc(ny*nray,sizeof(segy));
+ amp = (float *)calloc(nx*ny*nray,sizeof(float));
+ conv = (float *)calloc(nx*ny*nt,sizeof(float));
+
+ vmess("Depth level %li out of %li",iy+1,nshots);
+ sprintf(fins,"z%li",iy*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbr,fins,fer);
+ readSnapData3D(fin2, time, hdr_time, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+ sprintf(fin2,"%s%s%s",fba,fins,fea);
+ readSnapData3D(file_amp, amp, hdr_amp, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+ for (is = 0; is < nray; is++) {
+ gx[is*nshots+iy] = hdr_time[is*ny].sx;
+ gy[is*nshots+iy] = hdr_time[is*ny].sy;
+ gz[is*nshots+iy] = hdr_time[is*ny].sdepth;
+ for (it = 0; it < nx*ny*nt; it++) {
+ conv[it] = gmin[iy*nx*ny*nt+it];
+ }
+ timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],delay,fmin,fmax);
+ for (ix = 0; ix < ny*nx; ix++) {
+ image[is*nshots+iy] += conv[ix*nt]*dx*dy*dt;
+ }
+ }
+
+ free(hdr_time); free(time); free(hdr_amp); free(amp); free(conv);
+ }
+ free(gmin); free(hdr_gmin);
+
+ /*----------------------------------------------------------------------------*
+ * Write out the data
+ *----------------------------------------------------------------------------*/
+ fp_out = fopen(file_out, "w+");
+
+ if (nshots>1) dz = ((float)(gz[1] - gz[0]))/1000.0;
+ else dz = 1.0;
+
+ for (it = 0; it < nray; it++) {
+ hdr_out[it].fldr = 1;
+ hdr_out[it].tracl = it+1;
+ hdr_out[it].tracf = it+1;
+ hdr_out[it].scalco = -1000;
+ hdr_out[it].scalel = -1000;
+ hdr_out[it].trid = 1;
+ hdr_out[it].ns = nshots;
+ hdr_out[it].trwf = nray;
+ hdr_out[it].ntr = nray;
+ hdr_out[it].f1 = (((float)gz[0])/1000.0);
+ hdr_out[it].f2 = (((float)gx[0])/1000.0);
+ hdr_out[it].dt = ((int)(dt*1E6));
+ hdr_out[it].d1 = roundf(dz*1000.0)/1000.0;
+ hdr_out[it].d2 = roundf(dx*1000.0)/1000.0;
+ hdr_out[it].gx = gx[it*nshots];
+ hdr_out[it].gy = gy[it*nshots];
+ }
+
+ ret = writeData3D(fp_out, &image[0], hdr_out, nshots, nray);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ fclose(fp_out);
+
+ free(image); free(hdr_out);
+
+ vmess("Wrote data");
+
+ return 0;
+}
+
+void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, tom, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (long)MIN((omin/deltom), (nfreq));
+ iomax = MIN((long)(omax/deltom), (nfreq));
+
+ cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ tom = om*-1.0*(time[ix]+delay[ix]);
+ cdatascl[ix*nfreq+iom].r = (cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom))/(amp[ix]*amp[ix]);
+ cdatascl[ix*nfreq+iom].i = (cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom))/(amp[ix]*amp[ix]);
+ }
+ }
+ free(cdata);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void 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
diff --git a/utils/zfpmar.h b/utils/zfpmar.h
new file mode 100644
index 0000000000000000000000000000000000000000..6442357821c1e54a233886c8f02eeb0d377f77f9
--- /dev/null
+++ b/utils/zfpmar.h
@@ -0,0 +1,65 @@
+#ifndef ZFPMAR_H
+#define ZFPMAR_H
+#define MARBYTES 48
+#define TOPBYTES 80
+
+/* Size of the different types
+long - 8 bytes
+int - 4 bytes
+char - 1 byte
+float - 4 bytes
+*/
+
+/* TYPEDEFS */
+typedef struct { /* zfpmar - headers for the compression and decompression of data for marchenko */
+
+ long nx; /* Number of samples in x-direction */
+
+ long ny; /* Number of samples in y-direction */
+
+ int sx; /* Source coordinate in x-direction */
+
+ int sy; /* Source coordinate in y-direction */
+
+ int sz; /* Source coordinate in z-direction */
+
+ int gx; /* Receiver coordinate in x-direction */
+
+ long gy; /* Receiver coordinate in y-direction */
+
+ long compsize; /* Size of the compressed data */
+
+} zfpmar;
+
+typedef struct { /* zfptop - headers for the compression and decompression of data for marchenko */
+
+ long ndim; /* Number of dimension is between 1 and 4 */
+
+ long nz; /* Number of samples in z-direction */
+
+ long ns; /* Number of shots */
+
+ long nt; /* Number of time samples */
+
+ float dx; /* Sampling distance in x-direction */
+
+ float dy; /* Sampling distance in x-direction */
+
+ float dz; /* Sampling distance in x-direction */
+
+ float scale; /* Scaling of the coordinates and the sampling */
+
+ double tolerance; /* Set the tolerance of the zfp (de)compression */
+
+ float fmin; /* Minimum frequency of the signal */
+
+ float fmax; /* Maximum frequency of the singal */
+
+ float fx; /* First location of receiver in x-direction */
+
+ float fy; /* First location of receiver in y-direction */
+
+ float fz; /* First location of receiver in y-direction */
+
+} zfptop;
+#endif
\ No newline at end of file