From f528715e9ec471f45cd23391ad7280afe3aeb12d Mon Sep 17 00:00:00 2001
From: "joeri.brackenhoff" <joeri.brackenhoff@login0.ogun.local>
Date: Wed, 30 Jan 2019 15:50:29 -0300
Subject: [PATCH] 3D
---
marchenko3D/Makefile | 12 +-
marchenko3D/applyMute3D.c | 115 ++++++
marchenko3D/fmute3D.c | 508 ++++++++++++++++++++++++
marchenko3D/marchenko3D.c | 812 ++++++++++++++++++++++++++++++++++++++
marchenko3D/readData3D.c | 54 +++
marchenko3D/synthesis3D.c | 275 +++++++++++++
marchenko3D/test.c | 214 ----------
7 files changed, 1770 insertions(+), 220 deletions(-)
create mode 100644 marchenko3D/applyMute3D.c
create mode 100644 marchenko3D/fmute3D.c
create mode 100644 marchenko3D/marchenko3D.c
create mode 100644 marchenko3D/readData3D.c
create mode 100644 marchenko3D/synthesis3D.c
delete mode 100755 marchenko3D/test.c
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index dd2432e..491a1c0 100644
--- a/marchenko3D/Makefile
+++ b/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
diff --git a/marchenko3D/applyMute3D.c b/marchenko3D/applyMute3D.c
new file mode 100644
index 0000000..cc6f86d
--- /dev/null
+++ b/marchenko3D/applyMute3D.c
@@ -0,0 +1,115 @@
+#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;
+}
+
diff --git a/marchenko3D/fmute3D.c b/marchenko3D/fmute3D.c
new file mode 100644
index 0000000..533c7eb
--- /dev/null
+++ b/marchenko3D/fmute3D.c
@@ -0,0 +1,508 @@
+#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
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
new file mode 100644
index 0000000..7323516
--- /dev/null
+++ b/marchenko3D/marchenko3D.c
@@ -0,0 +1,812 @@
+/*
+ * Copyright (c) 2017 by the Society of Exploration Geophysicists.
+ * For more information, go to http://software.seg.org/2017/00XX .
+ * You must read and accept usage terms at:
+ * http://software.seg.org/disclaimer.txt before use.
+ */
+
+#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))
+
+
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+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 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);
+int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn,
+float *zsyn, int *ixpos, int npos, int iter);
+int unique_elements(float *arr, int len);
+
+void name_ext(char *filename, char *extension);
+
+void applyMute(float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npos, int shift);
+
+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 readData(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);
+double wallclock_time(void);
+
+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);
+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 linearsearch(int *array, size_t N, int value);
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" MARCHENKO - Iterative Green's function and focusing functions retrieval",
+" ",
+" marchenko file_tinv= file_shot= [optional parameters]",
+" ",
+" Required parameters: ",
+" ",
+" file_tinv= ............... direct arrival from focal point: G_d",
+" file_shot= ............... Reflection response: R",
+" ",
+" Optional parameters: ",
+" ",
+" INTEGRATION ",
+" tap=0 .................... lateral taper focusing(1), shot(2) or both(3)",
+" ntap=0 ................... number of taper points at boundaries",
+" fmin=0 ................... minimum frequency in the Fourier transform",
+" fmax=70 .................. maximum frequency in the Fourier transform",
+" MARCHENKO ITERATIONS ",
+" niter=10 ................. number of iterations",
+" MUTE-WINDOW ",
+" above=0 .................. mute above(1), around(0) or below(-1) the first travel times of file_tinv",
+" shift=12 ................. number of points above(positive) / below(negative) travel time for mute",
+" hw=8 ..................... window in time samples to look for maximum in next trace",
+" smooth=5 ................. number of points to smooth mute with cosine window",
+" REFLECTION RESPONSE CORRECTION ",
+" scale=2 .................. scale factor of R for summation of Ni with G_d",
+" pad=0 .................... amount of samples to pad the reflection series",
+" reci=0 ................... 1; add receivers as shots 2; only use receivers as shot positions",
+" countmin=0 ............... 0.3*nxrcv; minumum number of reciprocal traces for a contribution",
+" OUTPUT DEFINITION ",
+" file_green= .............. output file with full Green function(s)",
+" file_gplus= .............. output file with G+ ",
+" file_gmin= ............... output file with G- ",
+" file_f1plus= ............. output file with f1+ ",
+" file_f1min= .............. output file with f1- ",
+" file_f2= ................. output file with f2 (=p+) ",
+" file_pplus= .............. output file with p+ ",
+" file_pmin= ............... output file with p- ",
+" file_iter= ............... output file with -Ni(-t) for each iteration",
+" verbose=0 ................ silent option; >0 displays info",
+" ",
+" ",
+" author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+ FILE *fp_out, *fp_f1plus, *fp_f1min;
+ FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin;
+ int i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
+ int size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
+ int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
+ int reci, countmin, mode, n2out, n3out, verbose, ntfft;
+ int iter, niter, tracf, *muteW;
+ int hw, smooth, above, shift, *ixpos, npos, ix;
+ int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
+ float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
+ double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
+ float d1, d2, d3, f1, f2, f3, fxsb, fxse, fysb, fyse, ft, fx, fy, *xsyn, *ysyn, dxsrc, dysrc;
+ float *green, *f2p, *pmin, *G_d, dt, dx, dy, dxs, dys, scl, mem;
+ float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
+ float xmin, xmax, ymin, ymax, scale, tsq, Q, f0;
+ float *ixmask, *iymask;
+ complex *Refl, *Fop;
+ char *file_tinv, *file_shot, *file_green, *file_iter;
+ char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
+ segy *hdrs_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ tsyn = tread = tfft = tcopy = 0.0;
+ t0 = wallclock_time();
+
+ if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
+ if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
+ if (!getparstring("file_f1plus", &file_f1plus)) file_f1plus = NULL;
+ if (!getparstring("file_f1min", &file_f1min)) file_f1min = NULL;
+ if (!getparstring("file_gplus", &file_gplus)) file_gplus = NULL;
+ if (!getparstring("file_gmin", &file_gmin)) file_gmin = NULL;
+ if (!getparstring("file_pplus", &file_f2)) file_f2 = NULL;
+ if (!getparstring("file_f2", &file_f2)) file_f2 = NULL;
+ if (!getparstring("file_pmin", &file_pmin)) file_pmin = NULL;
+ if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
+ if (!getparint("verbose", &verbose)) verbose = 0;
+ if (file_tinv == NULL && file_shot == NULL)
+ verr("file_tinv and file_shot cannot be both input pipe");
+ if (!getparstring("file_green", &file_green)) {
+ if (verbose) vwarn("parameter file_green not found, assume pipe");
+ file_green = NULL;
+ }
+ if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (!getparfloat("fmax", &fmax)) fmax = 70.0;
+ if (!getparint("reci", &reci)) reci = 0;
+ if (!getparfloat("scale", &scale)) scale = 2.0;
+ if (!getparfloat("tsq", &tsq)) tsq = 0.0;
+ if (!getparfloat("Q", &Q)) Q = 0.0;
+ if (!getparfloat("f0", &f0)) f0 = 0.0;
+ if (!getparint("tap", &tap)) tap = 0;
+ if (!getparint("ntap", &ntap)) ntap = 0;
+ if (!getparint("pad", &pad)) pad = 0;
+
+ if(!getparint("niter", &niter)) niter = 10;
+ if(!getparint("hw", &hw)) hw = 15;
+ if(!getparint("smooth", &smooth)) smooth = 5;
+ if(!getparint("above", &above)) above = 0;
+ if(!getparint("shift", &shift)) shift=12;
+
+ if (reci && ntap) vwarn("tapering influences the reciprocal result");
+
+/*================ Reading info about shot and initial operator sizes ================*/
+
+ ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
+ ret = getFileInfo3D(file_tinv, &n1, &n2, &n3, &ngath, &d1, &d2, &d3, &f1, &f2, &f3, &scl, &ntraces);
+ Nfoc = ngath;
+ nxs = n2;
+ nys = n3;
+ nts = n1;
+ dxs = d2;
+ dys = d3;
+ fxsb = f2;
+ fysb = f3;
+
+ ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
+ ret = getFileInfo3D(file_shot, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+ nshots = ngath;
+ assert (nxs*nys >= nshots);
+
+ if (!getparfloat("dt", &dt)) dt = d1;
+
+ ntfft = optncr(MAX(nt+pad, nts+pad));
+ 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);
+ if (!getparint("countmin", &countmin)) countmin = 0.3*nx*ny;
+
+/*================ Allocating all data arrays ================*/
+
+ Fop = (complex *)calloc(nys*nxs*nw*Nfoc,sizeof(complex));
+ green = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ f2p = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ pmin = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ f1plus = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ f1min = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ iRN = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ Ni = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ G_d = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ muteW = (int *)calloc(Nfoc*nys*nxs,sizeof(int));
+ trace = (float *)malloc(ntfft*sizeof(float));
+ tapersy = (float *)malloc(nxs*sizeof(float));
+ xrcvsyn = (float *)calloc(Nfoc*nys*nxs,sizeof(float)); // x-rcv postions of focal points
+ yrcvsyn = (float *)calloc(Nfoc*nys*nxs,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(nys*nxs,sizeof(int)); // x-position of source of shot in G_d domain (nxs*nys with dxs, dys)
+
+ Refl = (complex *)malloc(nw*ny*nx*nshots*sizeof(complex));
+ tapersh = (float *)malloc(nx*sizeof(float));
+ xrcv = (float *)calloc(nshots*ny*nx,sizeof(float)); // x-rcv postions of shots
+ yrcv = (float *)calloc(nshots*ny*nx,sizeof(float)); // x-rcv postions of shots
+ xsrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
+ ysrc = (float *)calloc(nshots,sizeof(float)); //x-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
+
+ if (reci!=0) {
+ reci_xsrc = (int *)malloc((nxs*nxs*nys*nys)*sizeof(int));
+ reci_xrcv = (int *)malloc((nxs*nxs*nys*nys)*sizeof(int));
+ isxcount = (int *)calloc(nxs*nys,sizeof(int));
+ ixmask = (float *)calloc(nxs*nys,sizeof(float));
+ }
+
+/*================ Read and define mute window based on focusing operator(s) ================*/
+/* G_d = p_0^+ = G_d (-t) ~ Tinv */
+
+ mode=-1; /* apply complex conjugate to read in data */
+ readTinvData3D(file_tinv, xrcvsyn, yrcvsyn, xsyn, ysyn, zsyn, xnxsyn, Nfoc,
+ nxs, nys, ntfft, mode, muteW, G_d, hw, verbose);
+ /* reading data added zero's to the number of time samples to be the same as ntfft */
+ nts = ntfft;
+
+ /* define tapers to taper edges of acquisition */
+ if (tap == 1 || tap == 3) {
+ for (j = 0; j < ntap; j++)
+ tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
+ for (j = ntap; j < nxs-ntap; j++)
+ tapersy[j] = 1.0;
+ for (j = nxs-ntap; j < nxs; j++)
+ tapersy[j] =(cos(PI*(j-(nxs-ntap))/ntap)+1)/2.0;
+ }
+ else {
+ for (j = 0; j < nxs; j++) tapersy[j] = 1.0;
+ }
+ if (tap == 1 || tap == 3) {
+ if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
+ for (l = 0; l < Nfoc; l++) {
+ for (k = 0; k < nys; k++) {
+ for (i = 0; i < nxs; i++) {
+ for (j = 0; j < nts; j++) {
+ G_d[l*nys*nxs*nts+k*nxs*nts+i*nts+j] *= tapersy[i];
+ }
+ }
+ }
+ }
+ }
+
+ /* check consistency of header values */
+ if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxsb = xrcvsyn[0];
+ if (yrcvsyn[0] != 0 || yrcvsyn[nys*nxs-1] != 0 ) fysb = yrcvsyn[0];
+ if (nxs>1) {
+ fxse = fxsb + (float)(nxs-1)*dxs;
+ dxf = (xrcvsyn[nys*nxs-1] - xrcvsyn[0])/(float)(nxs-1);
+ }
+ else {
+ fxse = fxsb;
+ dxs = 1.0;
+ dx = 1.0;
+ d2 = 1.0;
+ dxf = 1.0;
+ }
+ if (nys>1) {
+ fyse = fysb + (float)(nys-1)*dys;
+ dyf = (yrcvsyn[nys*nxs-1] - yrcvsyn[0])/(float)(nys-1);
+ }
+ else {
+ fyse = fysb;
+ dys = 1.0;
+ d3 = 1.0;
+ dy = 1.0;
+ dyf = 1.0;
+ }
+ if (NINT(dxs*1e3) != NINT(fabs(dxf)*1e3)) {
+ vmess("dx in hdr.d2 (%.3f) and hdr.gx (%.3f) not equal",d2, dxf);
+ if (dxf != 0) dxs = fabs(dxf);
+ vmess("dx in operator => %f", dxs);
+ }
+ if (NINT(dys*1e3) != NINT(fabs(dyf)*1e3)) {
+ vmess("dy in hdr.d3 (%.3f) and hdr.gy (%.3f) not equal",d3, dyf);
+ if (dyf != 0) dys = fabs(dyf);
+ vmess("dy in operator => %f", dys);
+ }
+
+/*================ Reading shot records ================*/
+
+ mode=1;
+ readShotData3D(file_shot, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw,
+ nw_low, nshots, nx, ny, ntfft, mode, scale, verbose);
+
+ tapersh = (float *)malloc(nx*sizeof(float));
+ if (tap == 2 || tap == 3) {
+ for (j = 0; j < ntap; j++)
+ tapersh[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
+ for (j = ntap; j < nx-ntap; j++)
+ tapersh[j] = 1.0;
+ for (j = nx-ntap; j < nx; j++)
+ tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
+ }
+ else {
+ for (j = 0; j < nx; j++) tapersh[j] = 1.0;
+ }
+ if (tap == 2 || tap == 3) {
+ if (verbose) vmess("Taper for shots applied ntap=%d", ntap);
+ for (l = 0; l < nshots; l++) {
+ for (j = 1; j < nw; j++) {
+ for (k = 0; k < ny; k++) {
+ for (i = 0; i < nx; i++) {
+ Refl[l*nx*ny*nw+j*nx*ny+k*nx+i].r *= tapersh[i];
+ Refl[l*nx*ny*nw+j*nx*ny+k*nx+i].i *= tapersh[i];
+ }
+ }
+ }
+ }
+ }
+ free(tapersh);
+
+ /* check consistency of header values */
+ nxshot = unique_elements(xsrc,nshots);
+ nyshot = nshots/nxshot;
+
+ fxf = xsrc[0];
+ if (nx > 1) dxf = (xrcv[ny*nx-1] - xrcv[0])/(float)(nx-1);
+ else dxf = d2;
+ if (NINT(dx*1e3) != NINT(fabs(dxf)*1e3)) {
+ vmess("dx in hdr.d2 (%.3f) and hdr.gx (%.3f) not equal",dx, dxf);
+ if (dxf != 0) dx = fabs(dxf);
+ else verr("gx hdrs not set");
+ vmess("dx used => %f", dx);
+ }
+ fyf = ysrc[0];
+ if (ny > 1) dyf = (yrcv[ny*nx-1] - yrcv[0])/(float)(ny-1);
+ else dyf = d3;
+ if (NINT(dy*1e3) != NINT(fabs(dyf)*1e3)) {
+ vmess("dy in hdr.d3 (%.3f) and hdr.gy (%.3f) not equal",dy, dyf);
+ if (dyf != 0) dy = fabs(dyf);
+ else verr("gy hdrs not set");
+ vmess("dy used => %f", dy);
+ }
+
+ dxsrc = (float)xsrc[1] - xsrc[0];
+ if (dxsrc == 0) {
+ vwarn("sx hdrs are not filled in!!");
+ dxsrc = dx;
+ }
+ dysrc = (float)ysrc[nxshot-1] - ysrc[0];
+ if (dysrc == 0) {
+ vwarn("sy hdrs are not filled in!!");
+ dysrc = dy;
+ }
+
+/*================ Check the size of the files ================*/
+
+ if (NINT(dxsrc/dx)*dx != NINT(dxsrc)) {
+ vwarn("x: source (%.2f) and receiver step (%.2f) don't match",dxsrc,dx);
+ if (reci == 2) vwarn("x: step used from operator (%.2f) ",dxs);
+ }
+ if (NINT(dysrc/dy)*dy != NINT(dysrc)) {
+ vwarn("y: source (%.2f) and receiver step (%.2f) don't match",dysrc,dy);
+ if (reci == 2) vwarn("y: step used from operator (%.2f) ",dys);
+ }
+ dxi = NINT(dxf/dxs);
+ if ((NINT(dxi*dxs) != NINT(dxf)) && verbose)
+ vwarn("dx in receiver (%.2f) and operator (%.2f) don't match",dx,dxs);
+ dyi = NINT(dyf/dys);
+ if ((NINT(dyi*dys) != NINT(dyf)) && verbose)
+ vwarn("dy in receiver (%.2f) and operator (%.2f) don't match",dy,dys);
+ if (nt != nts)
+ vmess("Time samples in shot (%d) and focusing operator (%d) are not equal",nt, nts);
+ if (verbose) {
+ vmess("Number of focusing operators = %d", Nfoc);
+ vmess("Number of receivers in focusop = x:%d y:%d total:%d", nxs, nys, nxs*nys);
+ vmess("number of shots = %d", nshots);
+ vmess("number of receiver/shot = x:%d y:%d total:%d", nx, ny, nx*ny);
+ vmess("first model position = x:%.2f y:%.2f", fxsb, fysb);
+ vmess("last model position = x:%.2f y:%.2f", fxse, fyse);
+ vmess("first source position = x:%.2f y:%.2f", fxf, fyf);
+ vmess("source distance = x:%.2f y:%.2f", dxsrc, dysrc);
+ vmess("last source position = x:%.2f y:%.2f", fxf+(nxshot-1)*dxsrc, fyf+(nyshot-1)*dysrc);
+ vmess("receiver distance = x:%.2f y:%.2f", dxf, dyf);
+ vmess("direction of increasing traces = x:%d y:%d", dxi, dyi);
+ vmess("number of time samples (nt,nts) = %d (%d,%d)", ntfft, nt, nts);
+ vmess("time sampling = %e ", dt);
+ if (file_green != NULL) vmess("Green output file = %s ", file_green);
+ if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
+ if (file_gplus != NULL) vmess("Gplus output file = %s ", file_gplus);
+ if (file_pmin != NULL) vmess("Pmin output file = %s ", file_pmin);
+ if (file_f2 != NULL) vmess("f2 (=pplus) output file = %s ", file_f2);
+ if (file_f1min != NULL) vmess("f1min output file = %s ", file_f1min);
+ if (file_f1plus != NULL)vmess("f1plus output file = %s ", file_f1plus);
+ if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
+ }
+
+/*================ initializations ================*/
+
+ if (reci) {
+ n2out = nxs;
+ n3out = nys;
+ }
+ else {
+ n2out = nxshot;
+ n3out = nyshot;
+ }
+ mem = Nfoc*n2out*n3out*ntfft*sizeof(float)/1048576.0;
+ if (verbose) {
+ vmess("number of output traces = x:%d y:%d total:%d", n2out, n3out, n2out*n3out);
+ vmess("number of output samples = %d", ntfft);
+ vmess("Size of output data/file = %.1f MB", mem);
+ }
+
+
+ /* dry-run of synthesis to get all x-positions calcalated by the integration */
+ synthesisPositions3D(nx, ny, nxs, nys, Nfoc, xrcv, yrcv, xsrc, ysrc, xnx,
+ fxse, fyse, fxsb, fysb, dxs, dys, nshots, nxshot, nyshot, ixpos, &npos, reci, verbose);
+ if (verbose) {
+ vmess("synthesisPosistions: nxshot=%d nyshot=%d nshots=%d npos=%d", nxshot, nyshot, nshots, npos);
+ }
+
+/*================ set variables for output data ================*/
+
+ n1 = nts; n2 = n2out; n3 = n3out;
+ f1 = ft; f2 = xrcvsyn[ixpos[0]]; f3 = yrcvsyn[ixpos[0]];
+ d1 = dt;
+ if (reci == 0) { // distance between sources in R
+ d2 = dxsrc;
+ d3 = dysrc;
+ }
+ else if (reci == 1) { // distance between traces in G_d
+ d2 = dxs;
+ d3 = dys;
+ }
+ else if (reci == 2) { // distance between receivers in R
+ d2 = dx;
+ d3 = dy;
+ }
+
+ hdrs_out = (segy *) calloc(n2*n3,sizeof(segy));
+ if (hdrs_out == NULL) verr("allocation for hdrs_out");
+ size = nys*nxs*nts;
+
+ for (k = 0; k < n3; k++) {
+ for (i = 0; i < n2; i++) {
+ hdrs_out[k*n2+i].ns = n1;
+ hdrs_out[k*n2+i].trid = 1;
+ hdrs_out[k*n2+i].dt = dt*1000000;
+ hdrs_out[k*n2+i].f1 = f1;
+ hdrs_out[k*n2+i].f2 = f2;
+ hdrs_out[k*n2+i].d1 = d1;
+ hdrs_out[k*n2+i].d2 = d2;
+ hdrs_out[k*n2+i].trwf = n2out*n3out;
+ hdrs_out[k*n2+i].scalco = -1000;
+ hdrs_out[k*n2+i].gx = NINT(1000*(f2+i*d2));
+ hdrs_out[k*n2+i].gy = NINT(1000*(f3+k*d3));
+ hdrs_out[k*n2+i].scalel = -1000;
+ hdrs_out[k*n2+i].tracl = k*n2+i+1;
+ }
+ }
+ t1 = wallclock_time();
+ tread = t1-t0;
+
+/*================ initialization ================*/
+
+ memcpy(Ni, G_d, Nfoc*nys*nxs*ntfft*sizeof(float));
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ ix = ixpos[i]; /* select the traces that have an output trace after integration */
+ f2p[l*nys*nxs*nts+i*nts+j] = G_d[l*nys*nxs*nts+ix*nts+j];
+ f1plus[l*nys*nxs*nts+i*nts+j] = G_d[l*nys*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ f2p[l*nys*nxs*nts+i*nts+j] = G_d[l*nys*nxs*nts+ix*nts+j];
+ f1plus[l*nys*nxs*nts+i*nts+j] = G_d[l*nys*nxs*nts+ix*nts+j];
+ }
+ }
+ }
+
+/*================ start Marchenko iterations ================*/
+
+ for (iter=0; iter<niter; iter++) {
+
+ t2 = wallclock_time();
+
+/*================ construction of Ni(-t) = - \int R(x,t) Ni(t) ================*/
+
+ synthesis3D(Refl, Fop, Ni, iRN, nx, ny, nt, nxs, nys, nts, dt, 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, reci, nshots,
+ nxshot, nyshot, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
+ ixmask, verbose);
+
+ t3 = wallclock_time();
+ tsyn += t3 - t2;
+
+ if (file_iter != NULL) {
+ writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs*nys, d2, f2, n2out*n3out, Nfoc, xsyn, zsyn, ixpos, npos, iter);
+ }
+ /* N_k(x,t) = -N_(k-1)(x,-t) */
+ /* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
+ for (l = 0; l < Nfoc; l++) {
+ energyNi = 0.0;
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ ix = ixpos[i];
+ Ni[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+ix*nts+j];
+ pmin[l*nys*nxs*nts+i*nts+j] += iRN[l*nys*nxs*nts+ix*nts+j];
+ energyNi += iRN[l*nys*nxs*nts+ix*nts+j]*iRN[l*nys*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ Ni[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+ix*nts+nts-j];
+ pmin[l*nys*nxs*nts+i*nts+j] += iRN[l*nys*nxs*nts+ix*nts+j];
+ energyNi += iRN[l*nys*nxs*nts+ix*nts+j]*iRN[l*nys*nxs*nts+ix*nts+j];
+ }
+ }
+ if (iter==0) energyN0 = energyNi;
+ if (verbose >=2) vmess(" - iSyn %d: Ni at iteration %d has energy %e; relative to N0 %e",
+ l, iter, sqrt(energyNi), sqrt(energyNi/energyN0));
+ }
+
+ /* apply mute window based on times of direct arrival (in muteW) */
+ applyMute(Ni, muteW, smooth, above, Nfoc, nxs*nys, nts, ixpos, npos, shift);
+
+ /* update f2 */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f2p[l*nys*nxs*nts+i*nts+j] += Ni[l*nys*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f2p[l*nys*nxs*nts+i*nts+j] += Ni[l*nys*nxs*nts+i*nts+j];
+ }
+ }
+ }
+
+ if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f1min[l*nys*nxs*nts+i*nts+j] -= Ni[l*nys*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nys*nxs*nts+i*nts+j] -= Ni[l*nys*nxs*nts+i*nts+nts-j];
+ }
+ }
+ }
+ }
+ else {/* odd iterations update: => f_1^+(t) */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f1plus[l*nys*nxs*nts+i*nts+j] += Ni[l*nys*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1plus[l*nys*nxs*nts+i*nts+j] += Ni[l*nys*nxs*nts+i*nts+j];
+ }
+ }
+ }
+ }
+
+ t2 = wallclock_time();
+ tcopy += t2 - t3;
+
+ if (verbose) vmess("*** Iteration %d finished ***", iter);
+
+ } /* end of iterations */
+
+ free(Ni);
+ free(G_d);
+
+ /* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ /* set green to zero if mute-window exceeds nt/2 */
+ if (muteW[l*nys*nxs+ixpos[i]] >= nts/2) {
+ memset(&green[l*nys*nxs*nts+i*nts],0, sizeof(float)*nt);
+ continue;
+ }
+ green[l*nys*nxs*nts+i*nts+j] = f2p[l*nys*nxs*nts+i*nts+j] + pmin[l*nys*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ green[l*nys*nxs*nts+i*nts+j] = f2p[l*nys*nxs*nts+i*nts+nts-j] + pmin[l*nys*nxs*nts+i*nts+j];
+ }
+ }
+ }
+
+ /* compute upgoing Green's function G^+,- */
+ if (file_gmin != NULL) {
+ Gmin = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+
+ /* use f1+ as operator on R in frequency domain */
+ mode=1;
+ synthesis3D(Refl, Fop, f1plus, iRN, nx, ny, nt, nxs, nys, nts, dt, 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, reci, nshots,
+ nxshot, nyshot, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
+ ixmask, verbose);
+
+ /* compute upgoing Green's G^-,+ */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j=0;
+ ix = ixpos[i];
+ Gmin[l*nys*nxs*nts+i*nts+j] = iRN[l*nys*nxs*nts+ix*nts+j] - f1min[l*nys*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Gmin[l*nys*nxs*nts+i*nts+j] = iRN[l*nys*nxs*nts+ix*nts+j] - f1min[l*nys*nxs*nts+i*nts+j];
+ }
+ }
+ }
+ /* Apply mute with window for Gmin */
+ applyMute(Gmin, muteW, smooth, 1, Nfoc, nxs*nys, nts, ixpos, npos, shift);
+ } /* end if Gmin */
+
+ /* compute downgoing Green's function G^+,+ */
+ if (file_gplus != NULL) {
+ Gplus = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+
+ /* use f1-(*) as operator on R in frequency domain */
+ mode=-1;
+ synthesis3D(Refl, Fop, f1min, iRN, nx, ny, nt, nxs, nys, nts, dt, 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, reci, nshots,
+ nxshot, nyshot, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
+ ixmask, verbose);
+
+ /* compute downgoing Green's G^+,+ */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j=0;
+ ix = ixpos[i];
+ Gplus[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+ix*nts+j] + f1plus[l*nys*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Gplus[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+ix*nts+j] + f1plus[l*nys*nxs*nts+i*nts+nts-j];
+ }
+ }
+ }
+ } /* end if Gplus */
+
+ t2 = wallclock_time();
+ if (verbose) {
+ vmess("Total CPU-time marchenko = %.3f", t2-t0);
+ vmess("with CPU-time synthesis = %.3f", tsyn);
+ vmess("with CPU-time copy array = %.3f", tcopy);
+ vmess(" CPU-time fft data = %.3f", tfft);
+ vmess("and CPU-time read data = %.3f", tread);
+ }
+
+/*================ write output files ================*/
+
+
+ fp_out = fopen(file_green, "w+");
+ if (fp_out==NULL) verr("error on creating output file %s", file_green);
+ if (file_gmin != NULL) {
+ fp_gmin = fopen(file_gmin, "w+");
+ if (fp_gmin==NULL) verr("error on creating output file %s", file_gmin);
+ }
+ if (file_gplus != NULL) {
+ fp_gplus = fopen(file_gplus, "w+");
+ if (fp_gplus==NULL) verr("error on creating output file %s", file_gplus);
+ }
+ if (file_f2 != NULL) {
+ fp_f2 = fopen(file_f2, "w+");
+ if (fp_f2==NULL) verr("error on creating output file %s", file_f2);
+ }
+ if (file_pmin != NULL) {
+ fp_pmin = fopen(file_pmin, "w+");
+ if (fp_pmin==NULL) verr("error on creating output file %s", file_pmin);
+ }
+ if (file_f1plus != NULL) {
+ fp_f1plus = fopen(file_f1plus, "w+");
+ if (fp_f1plus==NULL) verr("error on creating output file %s", file_f1plus);
+ }
+ if (file_f1min != NULL) {
+ fp_f1min = fopen(file_f1min, "w+");
+ if (fp_f1min==NULL) verr("error on creating output file %s", file_f1min);
+ }
+
+
+ tracf = 1;
+ for (l = 0; l < Nfoc; l++) {
+ if (reci) {
+ f2 = fxsb;
+ f3 = fysb;
+ }
+ else {
+ f2 = fxf;
+ f3 = fyf;
+ }
+
+ for (k = 0; k < n3; k++) {
+ for (i = 0; i < n2; i++) {
+ hdrs_out[k*n2+i].fldr = l+1;
+ hdrs_out[k*n2+i].sx = NINT(xsyn[l]*1000);
+ hdrs_out[k*n2+i].sy = NINT(ysyn[l]*1000);
+ hdrs_out[k*n2+i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
+ hdrs_out[k*n2+i].tracf = tracf++;
+ hdrs_out[k*n2+i].selev = NINT(zsyn[l]*1000);
+ hdrs_out[k*n2+i].sdepth = NINT(-zsyn[l]*1000);
+ hdrs_out[k*n2+i].f1 = f1;
+ }
+ }
+
+ ret = writeData(fp_out, (float *)&green[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ if (file_gmin != NULL) {
+ ret = writeData(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_gplus != NULL) {
+ ret = writeData(fp_gplus, (float *)&Gplus[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_f2 != NULL) {
+ ret = writeData(fp_f2, (float *)&f2p[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_pmin != NULL) {
+ ret = writeData(fp_pmin, (float *)&pmin[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_f1plus != NULL) {
+ /* rotate to get t=0 in the middle */
+ for (i = 0; i < n2*n3; i++) {
+ hdrs_out[i].f1 = -n1*0.5*dt;
+ memcpy(&trace[0],&f1plus[l*size+i*nts],nts*sizeof(float));
+ for (j = 0; j < n1/2; j++) {
+ f1plus[l*size+i*nts+n1/2+j] = trace[j];
+ }
+ for (j = n1/2; j < n1; j++) {
+ f1plus[l*size+i*nts+j-n1/2] = trace[j];
+ }
+ }
+ ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_f1min != NULL) {
+ /* rotate to get t=0 in the middle */
+ for (i = 0; i < n2*n3; i++) {
+ hdrs_out[i].f1 = -n1*0.5*dt;
+ memcpy(&trace[0],&f1min[l*size+i*nts],nts*sizeof(float));
+ for (j = 0; j < n1/2; j++) {
+ f1min[l*size+i*nts+n1/2+j] = trace[j];
+ }
+ for (j = n1/2; j < n1; j++) {
+ f1min[l*size+i*nts+j-n1/2] = trace[j];
+ }
+ }
+ ret = writeData(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ }
+ ret = fclose(fp_out);
+ if (file_gplus != NULL) {ret += fclose(fp_gplus);}
+ if (file_gmin != NULL) {ret += fclose(fp_gmin);}
+ if (file_f2 != NULL) {ret += fclose(fp_f2);}
+ if (file_pmin != NULL) {ret += fclose(fp_pmin);}
+ if (file_f1plus != NULL) {ret += fclose(fp_f1plus);}
+ if (file_f1min != NULL) {ret += fclose(fp_f1min);}
+ if (ret < 0) verr("err %d on closing output file",ret);
+
+ if (verbose) {
+ t1 = wallclock_time();
+ vmess("and CPU-time write data = %.3f", t1-t2);
+ }
+
+/*================ free memory ================*/
+
+ free(hdrs_out);
+ free(tapersy);
+
+ exit(0);
+}
+
+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
diff --git a/marchenko3D/readData3D.c b/marchenko3D/readData3D.c
new file mode 100644
index 0000000..09b121a
--- /dev/null
+++ b/marchenko3D/readData3D.c
@@ -0,0 +1,54 @@
+#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;
+}
diff --git a/marchenko3D/synthesis3D.c b/marchenko3D/synthesis3D.c
new file mode 100644
index 0000000..8ab12a9
--- /dev/null
+++ b/marchenko3D/synthesis3D.c
@@ -0,0 +1,275 @@
+#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
diff --git a/marchenko3D/test.c b/marchenko3D/test.c
deleted file mode 100755
index 613719b..0000000
--- a/marchenko3D/test.c
+++ /dev/null
@@ -1,214 +0,0 @@
-#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
--
GitLab