diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index a2aad66ea839d90fe6338fef730c48d1d3d8ee45..8367a6d18441993c9e2e55ed9c2764dbecad3178 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -51,6 +51,8 @@ OBJJ3 = $(SRCJ3:%.c=%.o)
fmute3D: $(OBJJ3)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute3D $(OBJJ3) $(LIBS)
+OBJMS = $(SRCMS:%.c=%.o)
+
install: marchenko3D fmute3D
cp marchenko3D $B
cp fmute3D $B
diff --git a/marchenko3D/ampest3D2.c b/marchenko3D/ampest3D2.c
deleted file mode 100644
index db7dbac8bb8bef709c135f56adc273a9de250aee..0000000000000000000000000000000000000000
--- a/marchenko3D/ampest3D2.c
+++ /dev/null
@@ -1,94 +0,0 @@
-#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 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 Wmax, Amax, *wavelet, *At, scl, sclt;
- FILE *fp_wav;
- complex *Gdf, *f1df, *Af, *cwav, tmp;
- segy *hdrs_wav;
-
- nfreq = ntfft/2+1;
- scl = dx*dy;
- sclt = 1.0*dt/((float)ntfft);
-
- Gdf = (complex *)malloc(nfreq*sizeof(complex));
- f1df = (complex *)malloc(nfreq*sizeof(complex));
- Af = (complex *)calloc(nfreq,sizeof(complex));
- At = (float *)malloc(ntfft*sizeof(complex));
- wavelet = (float *)calloc(ntfft,sizeof(complex));
-
- if (file_wav == NULL) {
- Wmax = 1.0;
- }
- else {
- hdrs_wav = (segy *)calloc(1, sizeof(segy));
- fp_wav = fopen(file_wav, "r");
- readData3D(fp_wav, wavelet, hdrs_wav, 0);
- fclose(fp_wav);
- cwav = (complex *)calloc(nfreq,sizeof(complex));
- rc1fft(wavelet,cwav,(int)ntfft,-1);
- for (i=0; i<nfreq; i++) {
- tmp.r = cwav[i].r*cwav[i].r - cwav[i].i*cwav[i].i;
- tmp.i = 2*cwav[i].r*cwav[i].i;
- cwav[i].r = tmp.r*sclt;
- cwav[i].i = tmp.i*sclt;
- }
- cr1fft(cwav,wavelet,(int)ntfft,1);
- Wmax = maxest3D(wavelet,ntfft);
- vmess("Wmax: %.3e",Wmax);
- }
-
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- ix = ixpos[i];
- rc1fft(&Gd[l*nxs*nys*ntfft+i*ntfft],Gdf,(int)ntfft,-1);
- rc1fft(&f1d[l*nxs*nys*ntfft+ix*ntfft],f1df,(int)ntfft,-1);
- for (iw=0; iw<nfreq; iw++) {
- Af[iw].r += scl*sclt*(f1df[iw].r*Gdf[iw].r-f1df[iw].i*Gdf[iw].i);
- Af[iw].i += scl*sclt*(f1df[iw].r*Gdf[iw].i+f1df[iw].i*Gdf[iw].r);
- }
- }
- cr1fft(&Af[0],At,(int)ntfft,1);
- Amax = maxest3D(At,ntfft);
- ampest[l] = sqrtf(Wmax/Amax);
- memset(&Af[0],0.0, sizeof(float)*2*nfreq);
- }
- free(Gdf);free(f1df);free(Af);free(At);free(cwav); free(wavelet);
-
- 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;
-}
diff --git a/marchenko3D/ampest3D_backup.c b/marchenko3D/ampest3D_backup.c
deleted file mode 100644
index 08b480b642ec861b4b0d6222e5ed238c6b7a1a17..0000000000000000000000000000000000000000
--- a/marchenko3D/ampest3D_backup.c
+++ /dev/null
@@ -1,298 +0,0 @@
-#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
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index dc0ad063b7e38f87d92204eaab2e86779b3ef94e..b54ae64e65fc9f9dcd8e8a570a2bbcbe4e65fe90 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -786,7 +786,7 @@ int main (int argc, char **argv)
hdrs_Nfoc[l*nxim+j].fldr = 1;
hdrs_Nfoc[l*nxim+j].tracl = 1;
hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
- hdrs_Nfoc[l*nxim+j].trid = 1;
+ hdrs_Nfoc[l*nxim+j].trid = 2;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
@@ -832,7 +832,7 @@ int main (int argc, char **argv)
hdrs_Nfoc[l*nxim+j].fldr = 1;
hdrs_Nfoc[l*nxim+j].tracl = 1;
hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
- hdrs_Nfoc[l*nxim+j].trid = 1;
+ hdrs_Nfoc[l*nxim+j].trid = 2;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
@@ -879,7 +879,7 @@ int main (int argc, char **argv)
hdrs_Nfoc[l*nxim+j].fldr = i+1;
hdrs_Nfoc[l*nxim+j].tracl = 1;
hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
- hdrs_Nfoc[l*nxim+j].trid = 1;
+ hdrs_Nfoc[l*nxim+j].trid = 2;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
diff --git a/marchenko3D/marchenko3D_backup.c b/marchenko3D/marchenko3D_backup.c
deleted file mode 100644
index 42ca94123ee53f6b75fa017fed0b9e8a8796ed73..0000000000000000000000000000000000000000
--- a/marchenko3D/marchenko3D_backup.c
+++ /dev/null
@@ -1,979 +0,0 @@
-/*
- * 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) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
-
-
-
-#ifndef COMPLEX
-typedef struct _complexStruct { /* complex number */
- float r,i;
-} complex;
-#endif/* complex */
-
-long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, long *xnx, complex *cdata,
- long nw, long nw_low, long nshots, long nx, long ny, long ntfft, long mode, float scale, long verbose);
-long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
- long *xnx, long Nfoc, long nx, long ny, long ntfft, long mode, long *maxval, float *tinv, long hw, long verbose);
-long unique_elements(float *arr, long len);
-
-void name_ext(char *filename, char *extension);
-
-void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
-
-void applyMute3D(float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nt, long *xrcvsyn, long npos, long shift);
-
-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);
-long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
-long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
-long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
-double wallclock_time(void);
-
-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);
-
-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);
-
-void synthesisPositions3D(long nx, long ny, long nxs, long nys, long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc,
- long *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys, long nshots, long nxsrc, long nysrc,
- long *ixpos, long *npos, long reci, long verbose);
-void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, long nx, long ny, long nt, long nxs, long nys, long nts, float dt,
- float *xsyn, float *ysyn, long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, long *xnx,
- float fxse, float fxsb, float fyse, float fysb, float dxs, float dys, float dxsrc, float dysrc,
- float dx, float dy, long ntfft, long nw, long nw_low, long nw_high, long mode, long reci, long nshots, long nxsrc, long nysrc,
- long *ixpos, long npos, double *tfft, long *isxcount, long *reci_xsrc, long *reci_xrcv, float *ixmask, long verbose);
-
-void imaging3D(float *Image, float *Gmin, float *f1plus, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
-
-void homogeneousg3D(float *HomG, float *green, float *f2, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
-
-long linearsearch(long *array, size_t N, long value);
-
-/*********************** self documentation **********************/
-char *sdoc[] = {
-" ",
-" MARCHENKO3D - Iterative Green's function and focusing functions retrieval in 3D",
-" ",
-" marchenko3D file_tinv= file_shot= [optional parameters]",
-" ",
-" Required parameters: ",
-" ",
-" file_tinv= ............... direct arrival from focal point: G_d",
-" file_ray= ................ direct arrival from raytimes",
-" 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 ",
-" file_amp= ................ amplitudes for the raytime estimation",
-" 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",
-" HOMOGENEOUS GREEN'S FUNCTION RETRIEVAL OPTIONS ",
-" file_inp= ................ Input source function for the retrieval",
-" scheme=0 ................. Scheme for homogeneous Green's function retrieval",
-" .......................... Scheme for homogeneous Green's function retrieval",
-" kxwfilt=0 ................ Apply a dip filter before integration",
-" 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_imag= ............... output file with image ",
-" file_homg= ............... output file with homogeneous Green's function ",
-" file_ampscl= ............. output file with estimated amplitudes ",
-" 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)",
-" author : Joeri Brackenhoff : 2019 (j.a.brackenhoff@tudelft.nl)",
-" ",
-NULL};
-/**************** end self doc ***********************************/
-
-int main (int argc, char **argv)
-{
- FILE *fp_out, *fp_f1plus, *fp_f1min, *fp_imag;
- FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin, *fp_amp;
- long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
- long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
- long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
- long reci, countmin, mode, n2out, n3out, verbose, ntfft;
- long iter, niter, tracf, *muteW, ampest;
- long hw, smooth, above, shift, *ixpos, npos, ix, nzim, nxim, nyim;
- long 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, *tmpdata;
- float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim, dyim, dxim;
- complex *Refl, *Fop;
- char *file_tinv, *file_shot, *file_green, *file_iter, *file_imag, *file_ampscl;
- char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
- char *file_ray, *file_amp, *file_wav;
- segy *hdrs_out, *hdrs_Nfoc;
-
- 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_ray", &file_ray)) file_ray = NULL;
- if (!getparstring("file_amp", &file_amp)) file_amp = 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 (!getparstring("file_wav", &file_wav)) file_wav = NULL;
- if (!getparstring("file_imag", &file_imag)) file_imag = NULL;
- if (!getparstring("file_ampscl", &file_ampscl)) file_ampscl = NULL;
- if (!getparlong("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 (!getparlong("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 (!getparlong("tap", &tap)) tap = 0;
- if (!getparlong("ntap", &ntap)) ntap = 0;
- if (!getparlong("pad", &pad)) pad = 0;
- if (!getparlong("ampest", &est)) ampest = 0;
-
- if(!getparlong("niter", &niter)) niter = 10;
- if(!getparlong("hw", &hw)) hw = 15;
- if(!getparlong("smooth", &smooth)) smooth = 5;
- if(!getparlong("above", &above)) above = 0;
- if(!getparlong("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 */
- if (file_ray!=NULL) {
- ret = getFileInfo3D(file_ray, &n2, &n1, &n3, &ngath, &d2, &d1, &d3, &f2, &f1, &f3, &scl, &ntraces);
- Nfoc = ngath;
- nxs = n2;
- nys = n3;
- nts = n1;
- dxs = d2*1e3;
- dys = d3;
- fxsb = f2;
- fysb = f3;
- }
- else {
- 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 = loptncr(MAX(nt+pad, nts+pad));
- nfreq = ntfft/2+1;
- nw_low = (long)MIN((fmin*ntfft*dt), nfreq-1);
- nw_low = MAX(nw_low, 1);
- nw_high = MIN((long)(fmax*ntfft*dt), nfreq-1);
- nw = nw_high - nw_low + 1;
- scl = 1.0/((float)ntfft);
- if (!getparlong("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 = (long *)calloc(Nfoc*nys*nxs,sizeof(long));
- 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 = (long *)calloc(Nfoc,sizeof(long)); // number of traces per focal point
- ixpos = (long *)calloc(nys*nxs,sizeof(long)); // 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 = (long *)calloc(nshots,sizeof(long)); // number of traces per shot
-
- if (reci!=0) {
- reci_xsrc = (long *)malloc((nxs*nxs*nys*nys)*sizeof(long));
- reci_xrcv = (long *)malloc((nxs*nxs*nys*nys)*sizeof(long));
- isxcount = (long *)calloc(nxs*nys,sizeof(long));
- 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 */
- if (file_ray!=NULL) {
- makeWindow3D(file_ray, file_amp, file_wav, dt, xrcvsyn, yrcvsyn, xsyn, ysyn, zsyn,
- xnxsyn, Nfoc, nxs, nys, ntfft, muteW, G_d, verbose);
- }
- else {
- 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;
-
- /*Determine the shape of the focal positions*/
- nzim = unique_elements(zsyn,Nfoc);
- if (nzim>1) dzim = zsyn[1]-zsyn[0];
- else dzim = 1.0;
- nyim = unique_elements(ysyn,Nfoc);
- nxim = unique_elements(xsyn,Nfoc);
-
-
- /* 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=%li", 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 = (fxse - fxsb)/(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 = (fyse - fysb)/(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=%li", 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 (%li) and focusing operator (%li) are not equal",nt, nts);
- if (verbose) {
- vmess("Number of focusing operators = %li", Nfoc);
- vmess("Number of receivers in focusop = x:%li y:%li total:%li", nxs, nys, nxs*nys);
- vmess("number of shots = %li", nshots);
- vmess("number of receiver/shot = x:%li y:%li total:%li", 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:%li y:%li", dxi, dyi);
- vmess("number of time samples (nt,nts) = %li (%li,%li)", 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:%li y:%li total:%li", n2out, n3out, n2out*n3out);
- vmess("number of output samples = %li", 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=%li nyshot=%li nshots=%li npos=%li", 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 %li: Ni at iteration %li 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) */
- applyMute3D(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 %li finished ***", iter);
-
- } /* end of iterations */
-
- free(Ni);
- if (ampest < 1) 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];
- }
- }
- }
- applyMute3D(green, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
-
- /* compute upgoing Green's function G^+,- */
- if (file_gmin != NULL || file_imag!= 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 */
- applyMute3D(Gmin, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
- } /* end if Gmin */
-
- /* compute downgoing Green's function G^+,+ */
- if (file_gplus != NULL || ampest > 0) {
- 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];
- }
- }
- }
- /* Apply mute with window for Gplus */
- applyMute3D(Gplus, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
- } /* end if Gplus */
-
- /* Estimate the amplitude of the Marchenko Redatuming */
- if (ampest>0) {
- if (verbose>0) vmess("Estimating amplitude scaling");
-
- // Allocate memory and copy data
- ampscl = (float *)calloc(Nfoc,sizeof(float));
- Gd = (float *)calloc(Nfoc*nxs*nys*ntfft,sizeof(float));
- memcpy(Gd,Gplus,sizeof(float)*Nfoc*nxs*nys*ntfft);
- applyMute3D(Gd, muteW, smooth, 2, Nfoc, nxs*nys, nts, ixpos, npos, shift);
-
- // Determine amplitude and apply scaling
- AmpEst3D(G_d,Gd,ampscl,Nfoc,nxs,nys,ntfft,ixpos,npos,file_wav,dxs,dys,dt);
- for (l=0; l<Nfoc; l++) {
- for (j=0; j<nxs*nys*nts; j++) {
- green[l*nxs*nts+j] *= ampscl[l];
- if (file_gplus != NULL) Gplus[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_gmin != NULL || file_imag != NULL) Gmin[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f2 != NULL) f2p[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_pmin != NULL) pmin[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f1plus != NULL || file_imag != NULL) f1plus[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f1min != NULL) f1min[l*nxs*nys*nts+j] *= ampscl[l];
- }
- if (verbose>1) vmess("Amplitude of focal position %li is equal to %.3e",l,ampscl[l]);
- }
-
- if (file_ampscl!=NULL) { //Write the estimation of the amplitude to file
- hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
- for (l=0; l<nyim; l++){
- for (j=0; j<nxim; j++){
- hdrs_Nfoc[l*nxim+j].ns = nzim;
- hdrs_Nfoc[l*nxim+j].fldr = 1;
- hdrs_Nfoc[l*nxim+j].tracl = 1;
- hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
- hdrs_Nfoc[l*nxim+j].trid = 1;
- hdrs_Nfoc[l*nxim+j].scalco = -1000;
- hdrs_Nfoc[l*nxim+j].scalel = -1000;
- hdrs_Nfoc[l*nxim+j].sx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].gx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].gy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
- hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
- hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
- hdrs_Nfoc[l*nxim+j].d2 = xsyn[1]-xsyn[0];
- hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
- hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
- hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
- }
- }
- // Write the data
- fp_amp = fopen(file_ampscl, "w+");
- if (fp_amp==NULL) verr("error on creating output file %s", file_ampscl);
- ret = writeData3D(fp_amp, (float *)&scl[0], hdrs_Nfoc, nzim, nxim*nyim);
- if (ret < 0 ) verr("error on writing output file.");
- fclose(fp_amp);
- free(hdrs_Nfoc);
- free(ampscl);
- }
- free(Gd);
- if (file_gplus == NULL) free(Gplus);
- }
-
- /* Apply imaging*/
- if (file_imag!=NULL) {
- // Determine Image
- Image = (float *)calloc(Nfoc,sizeof(float));
- imaging3D(Image, Gmin, f1plus, nxs, nys, ntfft, dxs, dys, dt, Nfoc, verbose);
- if (file_gmin==NULL) free(Gmin);
- // Set headers
- hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
- for (l=0; l<nyim; l++){
- for (j=0; j<nxim; j++){
- hdrs_Nfoc[l*nxim+j].ns = nzim;
- hdrs_Nfoc[l*nxim+j].fldr = 1;
- hdrs_Nfoc[l*nxim+j].tracl = 1;
- hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
- hdrs_Nfoc[l*nxim+j].trid = 1;
- hdrs_Nfoc[l*nxim+j].scalco = -1000;
- hdrs_Nfoc[l*nxim+j].scalel = -1000;
- hdrs_Nfoc[l*nxim+j].sx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].gx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].gy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
- hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
- hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
- hdrs_Nfoc[l*nxim+j].d2 = xsyn[1]-xsyn[0];
- hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
- hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
- hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
- }
- }
- // Write out image
- fp_imag = fopen(file_imag, "w+");
- if (fp_imag==NULL) verr("error on creating output file %s", file_imag);
- ret = writeData3D(fp_imag, (float *)&Image[0], hdrs_Nfoc, nzim, nxim*nyim);
- if (ret < 0 ) verr("error on writing output file.");
- fclose(fp_imag);
- free(hdrs_Nfoc);
- free(Image);
- }
-
- 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 = writeData3D(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 = writeData3D(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 = writeData3D(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 = writeData3D(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 = writeData3D(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 = writeData3D(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 = writeData3D(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 %li 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);
-}
-
-long unique_elements(float *arr, long len)
-{
- if (len <= 0) return 0;
- long unique = 1;
- long 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;
-}
diff --git a/marchenko3D/synthesis3Dotavia.c b/marchenko3D/synthesis3Dotavia.c
deleted file mode 100644
index 4e4812ee82b9c69330dfabe25fb37f9339f592b0..0000000000000000000000000000000000000000
--- a/marchenko3D/synthesis3Dotavia.c
+++ /dev/null
@@ -1,319 +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>
-
-//External functions
-int omp_get_max_threads(void);
-int omp_get_num_threads(void);
-void omp_set_num_threads(int num_threads);
-
-//Kernels
-void setup_fops();
-
-
-#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 */
- //TODO: create a FFT kernel
- 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 (i = 0; i < inx; i++) {
- for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- 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;
-}
-
-void setup_fops(complex *Fop, float *Top, int nxys, int Nfoc, int nw, int npos, int ntfft, int *ixpos, int *first, float dxs, float dys, int nshots, int nxy, int nw_low, int *ircv, float *yrcv, float *xrcv, float fxsb, float fysb){
- int ix, idxs, idys, iloop, iw, k, i, l;
- complex *ctrace;
-
- ctrace = (complex *)calloc(ntfft,sizeof(complex));
-
- iloop = (*first ? npos : nxys)
-
- memset(&Fop[Nfoc*nxys*nw].r, 0, nxys*nw*2*sizeof(float));
-
- /* transform muted Ni (Top) to frequency domain, input for next iteration */
- //TODO: create a FFT kernel
- for (i = 0; i < iloop; i++) {
- /* set Fop to zero, so new operator can be defined within ixpos points */
- for (l = 0; l < Nfoc; l++) {
- rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
- ix = (*first ? i : 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;
- }
- }
- }
- if (*first) {
- idxs = 1.0/dxs;
- idys = 1.0/dys;
- ircv = (int *)malloc(nshots*nxy*sizeof(int));
- for (i = 0; i < nxy; i++) {
- for (k=0; k<nshots; k++) {
- ircv[k*nxy+i] = NINT((yrcv[k*nxy+i]-fysb)*idys)*nx+NINT((xrcv[k*nxy+i]-fxsb)*idxs);
- }
- }
- *first = 0;
- }
-
- free(ctrace);
-}
\ No newline at end of file
diff --git a/utils/HomG.c b/utils/HomG.c
index 570bf7ff84a6bfca26244d72a77376285242ca93..736b4f453b36cf35b121609a6c2a0b96166220dc 100755
--- a/utils/HomG.c
+++ b/utils/HomG.c
@@ -13,7 +13,7 @@
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -21,22 +21,21 @@ typedef struct _complexStruct { /* complex number */
} 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 disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+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);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
-int topdet(float *data, int nt);
-void conjugate(float *data, int nsam, int nrec, float dt);
-
-void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout);
-void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout);
-void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift);
-void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt, int shift);
-void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax);
-void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt);
-void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt);
-void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout);
+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 conjugate(float *data, long nsam, long nrec, float dt);
+
+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 corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift);
+void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt);
+void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt);
+void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, float *datout);
char *sdoc[] = {
" ",
@@ -71,20 +70,23 @@ int main (int argc, char **argv)
FILE *fp_in, *fp_shot, *fp_out;
char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
float *indata, *Ghom, *shotdata, *shotdata_jkz, rho, fmin, fmax;
- float dt, dx, t0, x0, xmin, xmax1, sclsxgx, f1, f2, dxrcv, dzrcv;
+ float dt, dy, dx, t0, y0, x0, xmin, xmax1, sclsxgx, f1, f2, f3, dxrcv, dyrcv, dzrcv;
float *conv, *conv2, *tmp1, *tmp2, cp, shift;
- int nshots, nt, nx, ntraces, ret, ix, it, is, ir, ig, file_det, nxs, nzs, verbose;
- int pos1, npos, zmax, inx, numb, dnumb, count, scheme, ntmax, ntshift, shift_num;
+ long nshots, nt, ny, nx, ntraces, ret, ix, iy, it, is, ir, ig, file_det, nxs, nys, nzs, verbose;
+ long pos1, npos, zmax, inx, numb, dnumb, count, scheme, ntmax, ntshift, shift_num;
segy *hdr_in, *hdr_out, *hdr_shot;
initargs(argc, argv);
requestdoc(1);
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
if (!getparstring("fin", &fin)) fin = NULL;
if (!getparstring("fshot", &fshot)) fshot = NULL;
if (!getparstring("fout", &fout)) fout = "out.su";
- if (!getparint("zmax", &zmax)) zmax = 0;
- if (!getparint("inx", &inx)) inx = 0;
+ if (!getparlong("zmax", &zmax)) zmax = 0;
+ if (!getparlong("inx", &inx)) inx = 0;
if (!getparfloat("zrcv", &f1)) f1 = 0;
if (!getparfloat("xrcv", &f2)) f2 = 0;
if (!getparfloat("dzrcv", &dzrcv)) dzrcv = -1;
@@ -94,35 +96,36 @@ int main (int argc, char **argv)
if (!getparfloat("fmin", &fmin)) fmin=0.0;
if (!getparfloat("fmax", &fmax)) fmax=100.0;
if (!getparfloat("shift", &shift)) shift=0.0;
- if (!getparint("numb", &numb)) numb=0;
- if (!getparint("dnumb", &dnumb)) dnumb=1;
- if (!getparint("scheme", &scheme)) scheme = 0;
- if (!getparint("ntmax", &ntmax)) ntmax = 0;
- if (!getparint("verbose", &verbose)) verbose = 0;
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=1;
+ if (!getparlong("scheme", &scheme)) scheme = 0;
+ if (!getparlong("ntmax", &ntmax)) ntmax = 0;
+ if (!getparlong("verbose", &verbose)) verbose = 0;
if (fin == NULL) verr("Incorrect f2 input");
if (fshot == NULL) verr("Incorrect Green input");
+ /*----------------------------------------------------------------------------*
+ * Split the filename so the number can be changed
+ *----------------------------------------------------------------------------*/
if (dnumb == 0) dnumb = 1;
-
- sprintf(fins,"z%d",numb);
-
+ sprintf(fins,"z%li",numb);
fp_in = fopen(fin, "r");
if (fp_in == NULL) {
verr("error on opening basefile=%s", fin);
}
fclose(fp_in);
-
ptr = strstr(fin,fins);
pos1 = ptr - fin + 1;
-
sprintf(fbegin,"%*.*s", pos1-1, pos1-1, fin);
sprintf(fend,"%s", fin+pos1+1);
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files to be read
+ *----------------------------------------------------------------------------*/
file_det = 1;
nzs=0;
-
while (file_det) {
- sprintf(fins,"z%d",nzs*dnumb+numb);
+ sprintf(fins,"z%li",nzs*dnumb+numb);
sprintf(fin,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin, "r");
if (fp_in == NULL) {
@@ -135,7 +138,7 @@ int main (int argc, char **argv)
break;
}
else {
- vmess("%d files detected",nzs);
+ vmess("%li files detected",nzs);
file_det = 0;
break;
}
@@ -155,10 +158,10 @@ int main (int argc, char **argv)
count=0;
npos = nxs*nzs;
- if (verbose) vmess("nxs: %d, nzs: %d",nxs,nzs);
+ if (verbose) vmess("nxs: %li, nzs: %li",nxs,nzs);
nshots = 0;
- getFileInfo(fshot, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax1, &sclsxgx, &ntraces);
+ getFileInfo3D(fshot, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &sclsxgx, &ntraces);
if (dxrcv < 0) dxrcv=dx;
if (dzrcv < 0) dzrcv=dx;
@@ -171,9 +174,9 @@ int main (int argc, char **argv)
if (fp_shot == NULL) {
verr("Could not open file");
}
- vmess("nt: %d nx: %d nshots: %d",nt,nx,nshots);
+ vmess("nt: %li nx: %li nshots: %li",nt,nx,nshots);
fclose(fp_shot);
- readSnapData(fshot, &shotdata[0], &hdr_shot[0], nshots, nx, nt, 0, nx, 0, nt);
+ readSnapData3D(fshot, &shotdata[0], &hdr_shot[0], nshots, nx, ny, nt, 0, nx, 0, ny, 0, nt);
hdr_out = (segy *)calloc(nxs,sizeof(segy));
@@ -215,14 +218,14 @@ int main (int argc, char **argv)
}
#pragma omp for
for (ir = 0; ir < nzs; ir++) {
- sprintf(fins,"z%d",ir*dnumb+numb);
+ sprintf(fins,"z%li",ir*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin2, "r");
if (fp_in == NULL) {
verr("Danger Will Robinson");
}
fclose(fp_in);
- readSnapData(fin2, &indata[0], &hdr_in[0], nxs, nx, nt, 0, nx, 0, nt);
+ readSnapData3D(fin2, &indata[0], &hdr_in[0], nxs, nx, ny, nt, 0, nx, 0, ny, 0, nt);
for (is=0;is<nxs;is++) {
if (scheme==0) { //Marchenko representation
depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
@@ -269,7 +272,7 @@ int main (int argc, char **argv)
for (ig=0; ig<nshots; ig++) {
convol(&shotdata[ig*nx*nt], &indata[is*nx*nt], conv, nx, nt, dt, -2);
timeDiff(conv, nt, nx, dt, fmin, fmax, -2);
- shift_num = ig*((int)(shift/dt));
+ shift_num = ig*((long)(shift/dt));
for (ix = 0; ix < nx; ix++) {
for (it = nt/2+1; it < nt; it++) {
conv[ix*nt+it] = 0.0;
@@ -293,7 +296,7 @@ int main (int argc, char **argv)
}
count+=1;
- if (verbose) vmess("Creating Homogeneous Green's function at depth %d from %d depths",count,nzs);
+ if (verbose) vmess("Creating Homogeneous Green's function at depth %li from %li depths",count,nzs);
}
free(conv); free(indata); free(hdr_in); free(conv2);
if (scheme==2) {
@@ -302,14 +305,14 @@ int main (int argc, char **argv)
}
free(shotdata);
- if (verbose) vmess("nxs: %d nxz: %d f1: %.7f",nxs,nzs,f1);
+ if (verbose) vmess("nxs: %li nxz: %li f1: %.7f",nxs,nzs,f1);
ntshift=0;
if (ntmax > 0) {
if (ntmax < nt) {
ntshift = (nt-ntmax)/2;
- if (verbose) vmess("Time shifted %d samples",ntshift);
+ if (verbose) vmess("Time shifted %li samples",ntshift);
nt=ntmax;
}
else {
@@ -340,7 +343,7 @@ int main (int argc, char **argv)
hdr_out[ix].gx = (int)roundf(f2 + (ix*hdr_out[ix].d2)*1000.0);
hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
}
- ret = writeData(fp_out, &Ghom[(ir+ntshift)*nxs*nzs], hdr_out, nzs, nxs);
+ ret = writeData3D(fp_out, &Ghom[(ir+ntshift)*nxs*nzs], hdr_out, nzs, nxs);
if (ret < 0 ) verr("error on writing output file.");
}
@@ -348,9 +351,9 @@ int main (int argc, char **argv)
return 0;
}
-void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift)
+void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift)
{
- int i, j, n, optn, nfreq, sign;
+ 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;
@@ -436,9 +439,9 @@ void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt
return;
}
-void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout)
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout)
{
- int it,ix;
+ long it,ix;
for (ix=0;ix<nrec;ix++) {
for (it=0;it<nsam;it++)
datout[ix*nsamout+it]=data[ix*nsam+it];
@@ -447,18 +450,18 @@ void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout)
}
}
-void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout)
+void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout)
{
- int it,ix;
+ long it,ix;
for (ix = 0; ix < nrec; ix++) {
for (it = 0 ; it < nsamout ; it++)
datout[ix*nsamout+it] = scl*data[ix*nsam+it];
}
}
-void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt, int shift)
+void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift)
{
- int i, j, n, optn, nfreq, sign;
+ 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;
@@ -537,9 +540,9 @@ void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt,
return;
}
-void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt)
+void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt)
{
- int optn, iom, iomin, iomax, nfreq, ix, sign;
+ long optn, iom, iomin, iomax, nfreq, ix, sign;
float omin, omax, deltom, om, df, *rdata, scl;
complex *cdata, *cdatascl;
@@ -563,9 +566,9 @@ void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax,
deltom = 2.*PI*df;
omin = 2.*PI*fmin;
omax = 2.*PI*fmax;
- iomin = (int)MIN((omin/deltom), (nfreq));
+ iomin = (long)MIN((omin/deltom), (nfreq));
iomin = MAX(iomin, 1);
- iomax = MIN((int)(omax/deltom), (nfreq));
+ iomax = MIN((long)(omax/deltom), (nfreq));
cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
if (cdatascl == NULL) verr("memory allocation error for cdatascl");
@@ -615,9 +618,9 @@ void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax,
return;
}
-void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt)
+void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt)
{
- int optn, iom, iomin, iomax, nfreq, ix, ikx, nkx, ikxmax;
+ long optn, iom, iomin, iomax, nfreq, ix, ikx, nkx, ikxmax;
float omin, omax, deltom, df, dkx, *rdata, kx, scl;
float kx2, kz2, kp2, kp;
complex *cdata, *cdatascl, kz, kzinv;
@@ -643,9 +646,9 @@ void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin,
omin = 2.*PI*fmin;
omax = 2.*PI*fmax;
- iomin = (int)MIN((omin/deltom), nfreq);
+ iomin = (long)MIN((omin/deltom), nfreq);
iomin = MAX(iomin, 0);
- iomax = MIN((int)(omax/deltom), nfreq);
+ iomax = MIN((long)(omax/deltom), nfreq);
cdatascl = (complex *)malloc(nfreq*nkx*sizeof(complex));
if (cdatascl == NULL) verr("memory allocation error for cdatascl");
@@ -667,7 +670,7 @@ void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin,
kp = (iom*deltom)/c;
kp2 = kp*kp;
- ikxmax = MIN((int)(kp/dkx), nkx/2);
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
for (ikx = 0; ikx < ikxmax; ikx++) {
kx = ikx*dkx;
@@ -698,7 +701,7 @@ void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin,
for (iom = iomin ; iom < iomax ; iom++) {
kp = iom*deltom/c;
kp2 = kp*kp;
- ikxmax = MIN((int)(kp/dkx), nkx/2);
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
for (ikx = 0; ikx < ikxmax; ikx++) {
kx = ikx*dkx;
kx2 = kx*kx;
@@ -737,9 +740,9 @@ void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin,
return;
}
-void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout)
+void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, float *datout)
{
- int it,ix;
+ long it,ix;
for (ix=0;ix<nrec;ix++) {
for (it=0;it<nsam;it++)
datout[ix*nsam+it]=data[ix*nsam+it];
@@ -751,9 +754,9 @@ void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float
datout[ix*nsam+it]=0.0;
}
}
-void conjugate(float *data, int nsam, int nrec, float dt)
+void conjugate(float *data, long nsam, long nrec, float dt)
{
- int optn, nfreq, j, ix, it, sign, ntdiff;
+ long optn, nfreq, j, ix, it, sign, ntdiff;
float *rdata, scl;
complex *cdata;
@@ -793,66 +796,4 @@ void conjugate(float *data, int nsam, int nrec, float dt)
free(rdata);
return;
-}
-
-void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax)
-{
- int 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 = (int)MIN((omin/deltom), (nfreq));
- iomax = MIN((int)(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*shift;
- cdatascl[ix*nfreq+iom].r = cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom);
- cdatascl[ix*nfreq+iom].i = cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom);
- }
- }
- 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;
-}
+}
\ No newline at end of file
diff --git a/utils/Makefile b/utils/Makefile
index d828166c229a26feba0ef68700865ead27cc9806..1efdfb8e3121657a60e30b6eefd592ec843b3ea6 100644
--- a/utils/Makefile
+++ b/utils/Makefile
@@ -6,7 +6,7 @@ include ../Make_include
#OPTC += -openmp
#OPTC += -g -O0
-ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d MuteSnap combine combine_induced reshape_su HomG Snap2Shot
+ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap combine combine_induced reshape_su HomG snap2shot
SRCM = \
makemod.c \
@@ -133,56 +133,56 @@ SRCT = ftr1d.c \
docpkge.c \
getpars.c
-SRCMS = MuteSnap.c \
- getFileInfo.c \
- writeData.c \
+SRCMS = mutesnap.c \
+ getFileInfo3D.c \
+ writeData3D.c \
verbosepkg.c \
getpars.c \
wallclock_time.c \
atopkge.c \
docpkge.c \
- readSnapData.c
+ readSnapData3D.c
SRCCO = combine.c \
- getFileInfo.c \
- writeData.c \
+ getFileInfo3D.c \
+ writeData3D.c \
wallclock_time.c \
getpars.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
- readSnapData.c
+ readSnapData3D.c
SRCCI = combine_induced.c \
- getFileInfo.c \
- writeData.c \
+ getFileInfo3D.c \
+ writeData3D.c \
wallclock_time.c \
getpars.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
- readSnapData.c
+ readSnapData3D.c
SRCRS = reshape_su.c \
- getFileInfo.c \
- writeData.c \
+ getFileInfo3D.c \
+ writeData3D.c \
getpars.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
- readSnapData.c
+ readSnapData3D.c
SRCHG = HomG.c \
- getFileInfo.c \
+ getFileInfo3D.c \
readData.c \
- writeData.c \
+ writeData3D.c \
verbosepkg.c \
getpars.c \
wallclock_time.c \
atopkge.c \
docpkge.c \
- readSnapData.c
+ readSnapData3D.c
SRCSS = snap2shot.c \
getFileInfo3D.c \
@@ -252,8 +252,8 @@ ftr1d: $(OBJT)
OBJMS = $(SRCMS:%.c=%.o)
-MuteSnap: $(OBJMS)
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o MuteSnap $(OBJMS) $(LIBS)
+mutesnap: $(OBJMS)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o mutesnap $(OBJMS) $(LIBS)
OBJCO = $(SRCCO:%.c=%.o)
@@ -277,10 +277,10 @@ HomG: $(OBJHG)
OBJSS = $(SRCSS:%.c=%.o)
-Snap2Shot: $(OBJSS)
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o Snap2Shot $(OBJSS) $(LIBS)
+snap2shot: $(OBJSS)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o snap2shot $(OBJSS) $(LIBS)
-install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d MuteSnap combine combine_induced reshape_su HomG Snap2Shot
+install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap combine combine_induced reshape_su HomG snap2shot
cp makemod $B
cp makewave $B
cp extendModel $B
@@ -292,15 +292,15 @@ install: makemod makewave extendModel fconv correigen green green3D basop syn2d
cp syn2d $B
cp mat2su $B
cp ftr1d $B
- cp MuteSnap $B
+ cp mutesnap $B
cp combine $B
cp combine_induced $B
cp reshape_su $B
cp HomG $B
- cp Snap2Shot $B
+ cp snap2shot $B
clean:
- rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) MuteSnap $(OBJMS) combine $(OBJCO) reshape_su $(OBJRS) combine_induced $(OBJCI) HomG $(OBJHG) Snap2Shot $(OBJSS)
+ rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) mutesnap $(OBJMS) combine $(OBJCO) reshape_su $(OBJRS) combine_induced $(OBJCI) HomG $(OBJHG) snap2shot $(OBJSS)
realclean: clean
- rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d $B/MuteSnap $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/Snap2Shot
+ rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d $B/mutesnap $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/snap2shot
diff --git a/utils/combine.c b/utils/combine.c
index a93a8814fcd21d41cee3106d861efeb018ef3959..45fb07d4d48698c796390148ed3d064c5da6c6b8 100755
--- a/utils/combine.c
+++ b/utils/combine.c
@@ -13,7 +13,7 @@
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -21,11 +21,12 @@ typedef struct _complexStruct { /* complex number */
} 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 disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+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);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+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);
char *sdoc[] = {
" ",
@@ -48,151 +49,164 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_in, *fp_out;
- char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
- float *indata, *outdata, *rtrace, fz, fx;
- float dt, dx, t0, x0, xmin, xmax, sclsxgx, dt2, dx2, t02, x02, xmin2, xmax2, sclsxgx2, dxrcv, dzrcv;
- int nshots, nt, nx, ntraces, nshots2, nt2, nx2, ntraces2, ix, it, is, iz, pos, ifile, file_det, nxs, nzs;
- int xcount, numb, dnumb, ret, nzmax, transpose, verbose;
- segy *hdr_in, *hdr_out;
+ FILE *fp_in, *fp_out;
+ char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
+ float *indata, *outdata, dt;
+ float dz, dy, dx, z0, y0, x0, scl;
+ long nt, nz, ny, nx, ntr, ix, iy, it, is, iz, pos, file_det, nzs;
+ long numb, dnumb, ret, nzmax, transpose, verbose, nxyz, sx, sy, sz;
+ segy *hdr_in, *hdr_out;
initargs(argc, argv);
requestdoc(1);
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
if (!getparstring("file_in", &fin)) fin = NULL;
if (!getparstring("file_out", &fout)) fout = "out.su";
- if (!getparint("numb", &numb)) numb=0;
- if (!getparint("dnumb", &dnumb)) dnumb=0;
- if (!getparint("nzmax", &nzmax)) nzmax=0;
- if (!getparint("verbose", &verbose)) verbose=0;
- if (!getparint("transpose", &transpose)) transpose=0;
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=0;
+ if (!getparlong("nzmax", &nzmax)) nzmax=0;
+ if (!getparlong("verbose", &verbose)) verbose=0;
+ if (!getparlong("transpose", &transpose)) transpose=0;
if (fin == NULL) verr("Incorrect downgoing 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,"%d",numb);
-
+ sprintf(numb1,"%li",numb);
ptr = strstr(fin,numb1);
pos = ptr - fin + 1;
-
sprintf(fbegin,"%*.*s", pos-1, pos-1, fin);
sprintf(fend,"%s", fin+pos);
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files that are present
+ *----------------------------------------------------------------------------*/
file_det = 1;
nzs=0;
-
- while (file_det) {
- sprintf(fins,"%d",nzs*dnumb+numb);
+ while (file_det) { // Check for a file with the filename
+ sprintf(fins,"%li",nzs*dnumb+numb);
sprintf(fin,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin, "r");
- if (fp_in == NULL) {
- if (nzs == 0) {
+ if (fp_in == NULL) { // If the filename does not exist
+ if (nzs == 0) { // The filename is wrong to begin with
verr("error on opening basefile=%s", fin);
}
- else if (nzs == 1) {
+ else if (nzs == 1) { // There is only a single file
vmess("1 file detected");
}
- else {
- vmess("%d files detected",nzs);
+ else { // Stop after the final file has been detected
+ vmess("%li files detected",nzs);
file_det = 0;
break;
}
}
fclose(fp_in);
nzs++;
- if (nzmax!=0 && nzs == nzmax) {
- vmess("%d files detected",nzs);
+ if (nzmax!=0 && nzs == nzmax) { // Stop if the amount of files exceed the indicated maximum
+ vmess("%li files detected",nzs);
file_det = 0;
break;
}
}
- sprintf(fins,"%d",numb);
+ /*----------------------------------------------------------------------------*
+ * Read in the first two files and determine the header values
+ * of the output
+ *----------------------------------------------------------------------------*/
+ sprintf(fins,"%li",numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
- nshots = 0;
- getFileInfo(fin2, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
-
- sprintf(fins,"%d",numb+dnumb);
- sprintf(fin2,"%s%s%s",fbegin,fins,fend);
- nshots = 0;
- getFileInfo(fin2, &nt2, &nx2, &nshots2, &dt2, &dx2, &t02, &x02, &xmin2, &xmax2, &sclsxgx2, &ntraces2);
-
- dxrcv=dx*1000;
- dzrcv=t02-t0;
-
- if (nshots==0) nshots=1;
- nxs = ntraces;
+ nt = 1;
+ getFileInfo3D(fin2, &nz, &nx, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
-
+ nxyz = nx*ny*nzs*nz;
- // ngath zijn het aantal schoten
- hdr_out = (segy *)calloc(nxs*nt,sizeof(segy));
- outdata = (float *)calloc(nxs*nzs*nt,sizeof(float));
- hdr_in = (segy *)calloc(nxs*nt,sizeof(segy));
- indata = (float *)calloc(nxs*nt,sizeof(float));
-
- readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
- nshots = hdr_in[nxs-1].fldr;
- if (transpose==0) {
- nxs = hdr_in[nxs-1].tracf;
- }
- else {
- nxs = hdr_in[nxs-1].ns;
+ if (verbose) {
+ vmess("number of time samples: %li", nt);
+ vmess("Number of virtual receivers: %li, x: %li, y: %li, z: %li",nxyz,nx,ny,nzs*nz);
+ vmess("Starting distance for x: %.3f, y: %.3f, z: %.3f",x0,y0,z0);
+ vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
}
- for (iz = 0; iz < nzs; iz++) {
- if (verbose) vmess("Depth:%d out of %d",iz+1,nzs);
- sprintf(fins,"%d",iz*dnumb+numb);
+ /*----------------------------------------------------------------------------*
+ * Read in a single file to determine if the header values match
+ * and allocate the data
+ *----------------------------------------------------------------------------*/
+ hdr_out = (segy *)calloc(nx*ny,sizeof(segy));
+ outdata = (float *)calloc(nxyz*nt,sizeof(float));
+ hdr_in = (segy *)calloc(nx*ny*nt,sizeof(segy));
+ indata = (float *)calloc(nx*ny*nz*nt,sizeof(float));
+
+ /*----------------------------------------------------------------------------*
+ * Combine the separate files
+ *----------------------------------------------------------------------------*/
+ for (is = 0; is < nzs; is++) {
+ if (verbose) vmess("Combining file %li out of %li",iz+1,nzs);
+ sprintf(fins,"%li",iz*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin2, "r");
if (fp_in == NULL) {
verr("Error opening file");
}
fclose(fp_in);
- if (transpose==0) {
- readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
- }
- else {
- readSnapData(fin2, &indata[0], &hdr_in[0], nshots, 1, nxs, 0, 1, 0, nxs);
- }
- if (iz==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
- if (iz==1) dzrcv=hdr_in[0].f1-fz;
- for (ix = 0; ix < nxs; ix++) {
- for (is = 0; is < nshots; is++) {
- outdata[is*nxs*nzs+ix*nzs+iz] = indata[is*nxs+ix];
+ readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
+ for (it = 0; it < nt; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*ny*nx*nz*nzs+iy*nx*nz*nzs+ix*nz*nzs+iz+is] = indata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
+ }
+ }
}
}
}
free(indata);
-
+ sx = hdr_in[0].sx;
+ sy = hdr_in[0].sy;
+ sz = hdr_in[0].sdepth;
+ dt = hdr_in[0].dt;
+ free(hdr_in);
+
+ /*----------------------------------------------------------------------------*
+ * Write out the data
+ *----------------------------------------------------------------------------*/
fp_out = fopen(fout, "w+");
- for (is = 0; is < nshots; is++) {
- for (ix = 0; ix < nxs; ix++) {
- hdr_out[ix].fldr = is+1;
- hdr_out[ix].tracl = is*nxs+ix+1;
- hdr_out[ix].tracf = ix+1;
- hdr_out[ix].scalco = -1000;
- hdr_out[ix].scalel = -1000;
- hdr_out[ix].sdepth = hdr_in[0].sdepth;
- hdr_out[ix].trid = 1;
- hdr_out[ix].ns = nzs;
- hdr_out[ix].trwf = nxs;
- hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
- hdr_out[ix].f1 = fz;
- hdr_out[ix].f2 = fx;
- hdr_out[ix].dt = dt*(1E6);
- hdr_out[ix].d1 = dzrcv;
- hdr_out[ix].d2 = dxrcv;
- hdr_out[ix].sx = (int)roundf(fx + (ix*hdr_out[ix].d2));
- hdr_out[ix].gx = (int)roundf(fx + (ix*hdr_out[ix].d2));
- hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ for (it = 0; it < nt; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ hdr_out[iy*nx+ix].fldr = it+1;
+ hdr_out[iy*nx+ix].tracl = it*ny*nx+iy*nx+ix+1;
+ hdr_out[iy*nx+ix].tracf = iy*nx+ix+1;
+ hdr_out[iy*nx+ix].scalco = -1000;
+ hdr_out[iy*nx+ix].scalel = -1000;
+ hdr_out[iy*nx+ix].sdepth = sz;
+ hdr_out[iy*nx+ix].selev = -sz;
+ hdr_out[iy*nx+ix].trid = 2;
+ hdr_out[iy*nx+ix].ns = nz*nzs;
+ hdr_out[iy*nx+ix].trwf = nx*ny;
+ hdr_out[iy*nx+ix].ntr = hdr_out[iy*nx+ix].fldr*hdr_out[iy*nx+ix].trwf;
+ hdr_out[iy*nx+ix].f1 = z0;
+ hdr_out[iy*nx+ix].f2 = x0;
+ hdr_out[iy*nx+ix].dt = dt;
+ hdr_out[iy*nx+ix].d1 = dz;
+ hdr_out[iy*nx+ix].d2 = dx;
+ hdr_out[iy*nx+ix].sx = sx;
+ hdr_out[iy*nx+ix].gx = (int)roundf(x0 + (ix*dx))*1000;
+ hdr_out[iy*nx+ix].sy = sy;
+ hdr_out[iy*nx+ix].gy = (int)roundf(y0 + (iy*dy))*1000;
+ hdr_out[iy*nx+ix].offset = (hdr_out[iy*nx+ix].gx - hdr_out[iy*nx+ix].sx)/1000.0;
+ }
}
- ret = writeData(fp_out, &outdata[is*nxs*nzs], hdr_out, nzs, nxs);
+ ret = writeData3D(fp_out, &outdata[it*nxyz], hdr_out, nz*nzs, nx*ny);
if (ret < 0 ) verr("error on writing output file.");
}
-
fclose(fp_out);
+ free(outdata); free(hdr_out);
+ vmess("Wrote data");
return 0;
-}
-
+}
\ No newline at end of file
diff --git a/utils/combine_induced.c b/utils/combine_induced.c
index 0b7601e48a0c994bdaad6108e90dca1f6ba71492..4b766d3f2c64b2e732e9268443d0e3354633d45b 100755
--- a/utils/combine_induced.c
+++ b/utils/combine_induced.c
@@ -13,7 +13,7 @@
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -21,11 +21,12 @@ typedef struct _complexStruct { /* complex number */
} 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 disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+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);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+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);
char *sdoc[] = {
" ",
@@ -48,89 +49,103 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_in, *fp_out;
- char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
- float *indata, *outdata, *rtrace, fz, fx, shift, dtshift, dt_time;
- float dt, dx, t0, x0, xmin, xmax, sclsxgx, dt2, dx2, t02, x02, xmin2, xmax2, sclsxgx2, dxrcv, dzrcv;
- int nshots, nt, nx, ntraces, nshots2, nt2, nx2, ntraces2, ix, it, is, iz, pos, ifile, file_det, nxs, nzs;
- int xcount, numb, dnumb, ret, nzmax, verbose, nshot_out, ishift, nshift;
- segy *hdr_in, *hdr_bin, *hdr_out;
+ FILE *fp_in, *fp_out;
+ char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
+ float *indata, *outdata, fz, fy, fx, shift, dtshift, dt_time;
+ float dt, dy, dx, t0, x0, y0, sclsxgx, dt2, dy2, dx2, t02, x02, y02, sclsxgx2, dxrcv, dyrcv, dzrcv;
+ long nshots, nt, ny, nx, ntraces, nshots2, nt2, ny2, nx2, ntraces2, ix, iy, it, is, iz, pos, file_det, nxs, nys, nzs;
+ long numb, dnumb, ret, nzmax, verbose, nshot_out, ishift, nshift;
+ segy *hdr_in, *hdr_bin, *hdr_out;
initargs(argc, argv);
requestdoc(1);
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
if (!getparstring("file_in", &fin)) fin = NULL;
if (!getparstring("file_out", &fout)) fout = "out.su";
- if (!getparint("numb", &numb)) numb=0;
- if (!getparint("dnumb", &dnumb)) dnumb=0;
- if (!getparint("nzmax", &nzmax)) nzmax=0;
- if (!getparint("verbose", &verbose)) verbose=0;
- if (!getparint("nshift", &nshift)) nshift=0;
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=0;
+ if (!getparlong("nzmax", &nzmax)) nzmax=0;
+ if (!getparlong("verbose", &verbose)) verbose=0;
+ if (!getparlong("nshift", &nshift)) nshift=0;
if (!getparfloat("shift", &shift)) shift=0.0;
if (!getparfloat("dtshift", &dtshift)) dtshift=0.0;
if (!getparfloat("dt_time", &dt_time)) dt_time=0.004;
if (fin == NULL) verr("Incorrect downgoing 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,"%d",numb);
-
+ sprintf(numb1,"%li",numb);
ptr = strstr(fin,numb1);
pos = ptr - fin + 1;
-
sprintf(fbegin,"%*.*s", pos-1, pos-1, fin);
sprintf(fend,"%s", fin+pos);
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files that are present
+ *----------------------------------------------------------------------------*/
file_det = 1;
nzs=0;
-
- while (file_det) {
- sprintf(fins,"%d",nzs*dnumb+numb);
+ while (file_det) { // Check for a file with the filename
+ sprintf(fins,"%li",nzs*dnumb+numb);
sprintf(fin,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin, "r");
- if (fp_in == NULL) {
- if (nzs == 0) {
+ if (fp_in == NULL) { // If the filename does not exist
+ if (nzs == 0) { // The filename is wrong to begin with
verr("error on opening basefile=%s", fin);
}
- else if (nzs == 1) {
+ else if (nzs == 1) { // There is only a single file
vmess("1 file detected");
}
- else {
- vmess("%d files detected",nzs);
+ else { // Stop after the final file has been detected
+ vmess("%li files detected",nzs);
file_det = 0;
break;
}
}
fclose(fp_in);
nzs++;
- if (nzmax!=0 && nzs == nzmax) {
- vmess("%d files detected",nzs);
+ if (nzmax!=0 && nzs == nzmax) { // Stop if the amount of files exceed the indicated maximum
+ vmess("%li files detected",nzs);
file_det = 0;
break;
}
}
- sprintf(fins,"%d",numb);
+ /*----------------------------------------------------------------------------*
+ * Read in the first two files and determine the header values
+ * of the output
+ *----------------------------------------------------------------------------*/
+ sprintf(fins,"%li",numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
nshots = 0;
- getFileInfo(fin2, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+ getFileInfo3D(fin2, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &sclsxgx, &ntraces);
- sprintf(fins,"%d",numb+dnumb);
+ sprintf(fins,"%li",numb+dnumb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
nshots = 0;
- getFileInfo(fin2, &nt2, &nx2, &nshots2, &dt2, &dx2, &t02, &x02, &xmin2, &xmax2, &sclsxgx2, &ntraces2);
+ getFileInfo3D(fin2, &nt2, &nx2, &ny2, &nshots2, &dt2, &dx2, &dy2, &t02, &x02, &y02, &sclsxgx2, &ntraces2);
dxrcv=dx*1000;
+ dyrcv=dy*1000;
dzrcv=t02-t0;
if (nshots==0) nshots=1;
nxs = ntraces;
- // ngath zijn het aantal schoten
+ /*----------------------------------------------------------------------------*
+ * Read in a single file to determine if the header values match
+ * and allocate the data
+ *----------------------------------------------------------------------------*/
hdr_bin = (segy *)calloc(nxs,sizeof(segy));
indata = (float *)calloc(nxs*nt,sizeof(float));
- readSnapData(fin2, &indata[0], &hdr_bin[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ readSnapData3D(fin2, &indata[0], &hdr_bin[0], nshots, nxs, ny, nt, 0, nxs, 0, ny, 0, nt);
nshots = hdr_bin[nxs-1].fldr;
nxs = hdr_bin[nxs-1].tracf;
@@ -140,30 +155,36 @@ int main (int argc, char **argv)
hdr_out = (segy *)calloc(nshot_out*nxs,sizeof(segy));
outdata = (float *)calloc(nshot_out*nxs*nt,sizeof(float));
- vmess("nshot_out:%d nxs=%d nt:%d",nshot_out,nxs,nt);
+ /*----------------------------------------------------------------------------*
+ * Write out the file info in case of verbose
+ *----------------------------------------------------------------------------*/
+ if (verbose) vmess("Number of virtual receivers: %li, nx=%li, ny=%li, nz=%li",nx*ny*nt,nx,ny,nt);
+ /*----------------------------------------------------------------------------*
+ * Parallel loop for reading in and combining the various shots
+ *----------------------------------------------------------------------------*/
#pragma omp parallel default(shared) \
- private(indata,iz,hdr_in,fins,fin2,fp_in,is,ix,it,ishift)
+ private(indata,iz,hdr_in,fins,fin2,fp_in,is,ix,iy,it,ishift)
{
indata = (float *)calloc(ntraces*nt,sizeof(float));
hdr_in = (segy *)calloc(ntraces,sizeof(segy));
#pragma omp for
for (iz = 0; iz < nzs; iz++) {
- if (verbose) vmess("Depth:%d out of %d",iz+1,nzs);
- sprintf(fins,"%d",iz*dnumb+numb);
+ if (verbose) vmess("Depth:%li out of %li",iz+1,nzs);
+ sprintf(fins,"%li",iz*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin2, "r");
if (fp_in == NULL) {
verr("Error opening file");
}
fclose(fp_in);
- readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ readSnapData3D(fin2, &indata[0], &hdr_in[0], nshots, nxs, ny, nt, 0, nxs, 0, ny, 0, nt);
if (iz==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
if (iz==1) dzrcv=hdr_in[0].f1-fz;
- //ishift = (int)((shift+(dtshift*((float)iz)))/dt_time);
+
ishift = nshift*iz;
- if (verbose) vmess("Shifting %d timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)iz)));
+ if (verbose) vmess("Shifting %li timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)iz)));
for (is = ishift; is < nshot_out; is++) {
for (ix = 0; ix < nxs; ix++) {
for (it = 0; it < nt; it++) {
@@ -175,6 +196,9 @@ int main (int argc, char **argv)
free(indata);free(hdr_in);
}
+ /*----------------------------------------------------------------------------*
+ * Write out the data to file
+ *----------------------------------------------------------------------------*/
fp_out = fopen(fout, "w+");
for (is = 0; is < nshot_out; is++) {
@@ -198,7 +222,7 @@ int main (int argc, char **argv)
hdr_out[ix].gx = (int)roundf(fx + (ix*hdr_out[ix].d2));
hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
}
- ret = writeData(fp_out, &outdata[is*nxs*nt], hdr_out, nt, nxs);
+ ret = writeData3D(fp_out, &outdata[is*nxs*nt], hdr_out, nt, nxs);
if (ret < 0 ) verr("error on writing output file.");
}
diff --git a/utils/mutesnap.c b/utils/mutesnap.c
new file mode 100644
index 0000000000000000000000000000000000000000..8676a680e401b102cd9f37e47d574f711c6e3137
--- /dev/null
+++ b/utils/mutesnap.c
@@ -0,0 +1,254 @@
+#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);
+long farrdet(float *array, long nt, float tol);
+long topdet(float *array, long nt);
+
+char *sdoc[] = {
+" ",
+" mutesnap - mute a file of snapshots ",
+" ",
+" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" fhom= .................... File containing the snapshot data that will be muted",
+" fsnap= ................... File containing the snapshot data that will determine the mute window",
+" ",
+" Optional parameters: ",
+" ",
+" fout=out.su .............. Filename of the output",
+" shift=5 .................. Shift from the maximum",
+" smooth=5 ................. Length of smoothing taper",
+" mode=0 ................... Determine first arrival by maximum (mode=0), first event above tol (mode=1) or by raytime (mode=2)",
+" tol=1 .................... Tolerance for the determination of first arrival if mode=1",
+" fray ..................... File containing the raytimes of the first arrivals",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_snap, *fp_hom, *fp_out;
+ char *fhom, *fsnap, *fout, *fray;
+ float *homdata, *snapdata, *outdata, *rtrace, *costaper, scl, tol, *timeval, dt;
+ float dxs, dys, dzs, scls, fzs, fxs, fys;
+ float dxh, dyh, dzh, sclh, fzh, fxh, fyh;
+ float dxrcv, dyrcv, dzrcv, dxpos, offset;
+ long nts, nxs, nys, nzs, ntrs, nth, nxh, nyh, nzh, ntrh;
+ long nxyz, nxy, ret, ix, iy, iz, it, ht, indrcv, shift, rmt, mode, smooth, verbose;
+ segy *hdr_hom, *hdr_snap, *hdrs_mute;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
+ if (!getparstring("fhom", &fhom)) fhom = NULL;
+ if (!getparstring("fsnap", &fsnap)) fsnap = NULL;
+ if (!getparstring("fout", &fout)) fout = "out.su";
+ if (!getparstring("fray", &fray)) fray = NULL;
+ if (!getparlong("shift", &shift)) shift = 5;
+ if (!getparlong("smooth", &smooth)) smooth = 5;
+ if (!getparlong("mode", &mode)) mode = 0;
+ if (!getparlong("verbose", &verbose)) verbose = 1;
+ if (!getparfloat("tol", &tol)) tol = 5;
+ if (fhom == NULL) verr("Incorrect G_hom input");
+ if (mode != 2) {
+ if (fsnap == NULL) verr("Incorrect snapshot input");
+ }
+ if (mode == 2) {
+ if (fray == NULL) verr("No filename for raytimes given");
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Configure the taper for smoothing
+ *----------------------------------------------------------------------------*/
+ if (smooth) {
+ costaper = (float *)malloc(smooth*sizeof(float));
+ scl = M_PI/((float)smooth);
+ for (it=0; it<smooth; it++) {
+ costaper[it] = 0.5*(1.0+cos((it+1)*scl));
+ }
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Determine the parameters of the files and determine if they match
+ *----------------------------------------------------------------------------*/
+ getFileInfo3D(fhom, &nzh, &nxh, &nyh, &nth, &dzh, &dxh, &dyh, &fzh, &fxh, &fyh, &sclh, &ntrh);
+
+ if (mode != 2) {
+ getFileInfo3D(fsnap, &nzs, &nxs, &nys, &nts, &dzs, &dxs, &dys, &fzs, &fxs, &fys, &scls, &ntrs);
+ if (nzs != nzh) verr("Unequal number of samples in z-direction, nzsnap=%li nzhom=%li",nzs,nzh);
+ if (nxs != nxh) verr("Unequal number of samples in x-direction, nxsnap=%li nxhom=%li",nxs,nxh);
+ if (nys != nyh) verr("Unequal number of samples in y-direction, nysnap=%li nyhom=%li",nys,nyh);
+ if (NINT(dzs*1000.0) != NINT(dzh*1000.0)) verr("Sampling distance unequal in z-direction, dzsnap=%f dzhom=%f",dzs,dzh);
+ if (NINT(dxs*1000.0) != NINT(dxh*1000.0)) verr("Sampling distance unequal in x-direction, dxsnap=%f dxhom=%f",dxs,dxh);
+ if (NINT(dys*1000.0) != NINT(dyh*1000.0)) verr("Sampling distance unequal in y-direction, dysnap=%f dyhom=%f",dys,dyh);
+ }
+
+ nxyz = nxh*nyh*nzh;
+ nxy = nxh*nyh;
+
+ if (verbose) {
+ vmess("Number of virtual receivers: %li, nz: %li, nx: %li, ny: %li",nxyz,nzh,nxh,nyh);
+ vmess("Sampling distance in direction of z: %.3f, x: %.3f, y: %.3f",dzh,dxh,dyh);
+ vmess("Number of time samples: %li",nth);
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Allocate data and read in raytime if required
+ *----------------------------------------------------------------------------*/
+ if (mode != 2) {
+ snapdata = (float *)malloc(nxyz*nth*sizeof(float));
+ hdr_snap = (segy *)calloc(nxy*nth,sizeof(segy));
+ }
+ homdata = (float *)malloc(nxyz*nth*sizeof(float));
+ hdr_hom = (segy *)calloc(nxy*nth,sizeof(segy));
+ ht = (long)ceil(nth/2);
+ rtrace = (float *)malloc(nth*sizeof(float));
+
+ if (mode != 2) {
+ readSnapData3D(fsnap, snapdata, hdr_snap, nts, nxs, nys, nzs, 0, nxs, 0, nys, 0, nzs);
+ if (verbose>1) vmess("Read file for mute determination");
+ }
+ readSnapData3D(fhom, homdata, hdr_hom, nth, nxh, nyh, nzh, 0, nxh, 0, nyh, 0, nzh);
+ if (verbose>1) vmess("Read file to be muted");
+
+ if (mode == 0) {
+ vmess("First arrival determined through maximum");
+ }
+ else if (mode == 1) {
+ vmess("First arrival determined through tolerance (=%.4f)",tol);
+ }
+ else if (mode == 2) {
+ vmess("First arrival determined through raytimes");
+ fp_snap = fopen(fray,"r");
+ if (fp_snap == NULL) {
+ verr("Could not open file");
+ }
+ fclose(fp_snap);
+ hdrs_mute = (segy *) calloc(nzh*nyh,sizeof(segy));
+ timeval = (float *)calloc(nxyz,sizeof(float));
+ readSnapData3D(fray, timeval, hdrs_mute, nzh, 1, nyh, nxh, 0, 1, 0, nyh, 0, nxh);
+ dt = hdr_hom[0].dt/1E6;
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Apply the muting to the data
+ *----------------------------------------------------------------------------*/
+ for (iy = 0; iy < nyh; iy++) {
+ for (iz = 0; iz < nzh; iz++) {
+ for (ix = 0; ix < nxh; ix++) {
+ if (mode != 2) {
+ for (it = 0; it < nth; it++) {
+ rtrace[it] = snapdata[it*nxyz+iy*nxh*nzh+ix*nzh+iz];
+ }
+ }
+ if (mode == 0) {
+ indrcv = ht - topdet(rtrace,nth);
+ }
+ else if (mode == 1) {
+ indrcv = ht - farrdet(rtrace,nth,tol);
+ }
+ else if (mode == 2) {
+ indrcv = (long)roundf(timeval[iz*nxh*nyh+iy*nxh+ix]/dt);
+ }
+ rmt = MAX(MIN(nth-indrcv,indrcv)-shift-smooth,0);
+ for (it = ht-rmt+1; it < ht+1; it++) {
+ if (it-(ht-rmt+1) < smooth) {
+ homdata[it*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
+ homdata[(nth-it)*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
+ }
+ else{
+ homdata[it*nxh*nzh+ix*nzh+iz] = 0.0;
+ homdata[(nth-it)*nxh*nzh+ix*nzh+iz] = 0.0;
+ }
+ }
+ }
+ if (verbose) vmess("Muting Homogeneous Green's function at depth %li from %li depths",iz+1,nzh);
+ }
+ }
+ free(rtrace);
+ if (smooth) free(costaper);
+ if (mode == 2) {
+ free(hdrs_mute); free(timeval);
+ }
+ if (mode != 2) {
+ free(snapdata);
+ free(hdr_snap);
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Write out the muted data
+ *----------------------------------------------------------------------------*/
+ fp_out = fopen(fout, "w+");
+
+ for (it = 0; it < nth; it++) {
+ ret = writeData3D(fp_out, &homdata[it*nxyz], &hdr_hom[it*nxy], nzh, nxy);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ free(homdata); free(hdr_hom);
+
+ fclose(fp_out);
+ vmess("Wrote Data");
+ return 0;
+}
+
+long topdet(float *array, long nt)
+{
+ float max;
+ long ht, it, ind;
+/* Find the maximum value in the array */
+ ht = nt/2;
+ ind = 0;
+ max = fabs(array[0]);
+ for (it = 1; it < ht; it++) {
+ if (fabs(array[it]) > max) {
+ ind = it;
+ max = fabs(array[it]);
+ }
+ }
+ return ind;
+}
+
+long farrdet(float *array, long nt, float tol)
+{
+ long ht, it, ind;
+/* Find the first value in the array above the tolerance value */
+ ht = nt/2;
+ ind = 0;
+ for (it = 0; it < ht; it++) {
+ if (fabs(array[it]) > tol) {
+ ind = it;
+ break;
+ }
+ }
+ return ind;
+}
\ No newline at end of file
diff --git a/utils/reshape_su.c b/utils/reshape_su.c
index 9c24040b2df3c4a0848a6205c1c783988669b7eb..6649029c3cdc2bbb8df1c2a726825bf6ddda8bab 100755
--- a/utils/reshape_su.c
+++ b/utils/reshape_su.c
@@ -13,7 +13,7 @@
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -21,14 +21,16 @@ typedef struct _complexStruct { /* complex number */
} 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 disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+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);
char *sdoc[] = {
" ",
-" reshape_su - interchange the 1st and 3rd dimension for SU file",
+" reshape_su - interchange the 1st and 4th dimension for SU file",
" ",
" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
" : Jan Thorbecke : (janth@xs4all.nl)",
@@ -44,80 +46,126 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_in, *fp_out;
- char *fin, *fout;
- float *indata, *outdata;
- float dt, dx, t0, x0, xmin, xmax, sclsxgx;
- int nshots, nt, nx, ntraces, ix, it, is, ir, ret, verbose;
- segy *hdr_in, *hdr_out, hdr;
+ FILE *fp_in, *fp_out;
+ char *fin, *fout;
+ float *indata, *outdata;
+ float d1, d2, d3, d4, f1, f2, f3, f4, scl;
+ long n1, n2, n3, n4, ntr, i1, i2, i3, i4, ret, verbose;
+ segy *hdr_in, *hdr_out, hdr;
initargs(argc, argv);
requestdoc(1);
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
if (!getparstring("file_in", &fin)) fin = NULL;
if (!getparstring("file_out", &fout)) fout = "out.su";
- if (!getparint("verbose", &verbose)) verbose = 0;
+ if (!getparlong("verbose", &verbose)) verbose = 1;
if (fin == NULL) verr("No input file specified");
- nshots = 0;
- getFileInfo(fin, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+ /*----------------------------------------------------------------------------*
+ * Determine the parameters of the files
+ *----------------------------------------------------------------------------*/
+ n4=1;
+ getFileInfo3D(fin, &n1, &n2, &n3, &n4, &d1, &d2, &d3, &f1, &f2, &f3, &scl, &ntr);
fp_in = fopen( fin, "r" );
ret = fread( &hdr, 1, TRCBYTES, fp_in );
assert(ret == TRCBYTES);
fclose(fp_in);
- if (nt==0) nt=hdr.ns;
- if (nx==0) nx=hdr.trwf;nshots=ntraces/nx;
- if (nshots==0) nshots=1;
-
- vmess("nx:%d nt:%d nshots:%d ntraces:%d",nx,nt,nshots,ntraces);
-
- // ngath zijn het aantal schoten
- hdr_out = (segy *)calloc(nx,sizeof(segy));
- outdata = (float *)calloc(nshots*nx*nt,sizeof(float));
- hdr_in = (segy *)calloc(nshots*nx,sizeof(segy));
- indata = (float *)calloc(nshots*nx*nt,sizeof(float));
-
- readSnapData(fin, &indata[0], &hdr_in[0], nshots, nx, nt, 0, nx, 0, nt);
-
- for (ir = 0; ir < nshots; ir++) {
- for (is = 0; is < nx; is++) {
- for (it = 0; it < nt; it++) {
- outdata[it*nx*nshots+is*nshots+ir] = indata[ir*nx*nt+is*nt+it];
+ if (hdr.trid==2) {
+ d4 = ((float)hdr.dt)/1E6;
+ f4 = -((float)(n4/2-1))*d4;
+ }
+ else if (hdr.trid==3) {
+ d4 = hdr.ungpow;
+ f4 = hdr.unscale;
+ f1 += ((float)(n1/2-1))*d1;
+ }
+
+ if (verbose) {
+ vmess("******************************* INPUT *******************************");
+ vmess("Number of samples in dimension 1:%li, 2:%li, 3:%li, 4:%li",n1,n2,n3,n4);
+ vmess("Sampling distance in dimension 1:%.3f, 2:%.3f, 3:%.3f, 4:%.3f",d1,d2,d3,d4);
+ vmess("Starting point in dimension 1:%.3f, 2:%.3f, 3:%.3f, 4:%.3f",f1,f2,f3,f4);
+ vmess("******************************* OUTPUT ******************************");
+ vmess("Number of samples in dimension 1:%li, 2:%li, 3:%li, 4:%li",n4,n2,n3,n1);
+ vmess("Sampling distance in dimension 1:%.3f, 2:%.3f, 3:%.3f, 4:%.3f",d4,d2,d3,d1);
+ vmess("Starting point in dimension 1:%.3f, 2:%.3f, 3:%.3f, 4:%.3f",f4,f2,f3,f1);
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Allocate the data
+ *----------------------------------------------------------------------------*/
+ hdr_out = (segy *)calloc(n2*n3,sizeof(segy));
+ outdata = (float *)calloc(n1*n2*n3*n4,sizeof(float));
+ hdr_in = (segy *)calloc(n2*n3*n4,sizeof(segy));
+ indata = (float *)calloc(n1*n2*n3*n4,sizeof(float));
+
+ /*----------------------------------------------------------------------------*
+ * Read in the data
+ *----------------------------------------------------------------------------*/
+ readSnapData3D(fin, indata, hdr_in, n4, n2, n3, n1, 0, n2, 0, n3, 0, n1);
+ if (verbose) vmess("Read data");
+
+ /*----------------------------------------------------------------------------*
+ * Reshape the data
+ *----------------------------------------------------------------------------*/
+ for (i4 = 0; i4 < n4; i4++) {
+ for (i3 = 0; i3 < n3; i3++) {
+ for (i2 = 0; i2 < n2; i2++) {
+ for (i1 = 0; i1 < n1; i1++) {
+ outdata[i1*n4*n2*n3+i3*n4*n2+i2*n4+i4] = indata[i4*n1*n2*n3+i3*n1*n2+i2*n1+i1];
+ }
}
}
- if (verbose) vmess("Reshaping shot %d out of %d shots",ir+1,nshots);
+ if (verbose>1) vmess("Reshaping dimension 4 number %li out of %li",i4+1,n4);
}
free(indata);
+ /*----------------------------------------------------------------------------*
+ * Write out the reshaped data
+ *----------------------------------------------------------------------------*/
fp_out = fopen(fout, "w+");
-
- for (is = 0; is < nt; is++) {
- for (ix = 0; ix < nx; ix++) {
- hdr_out[ix].fldr = is+1;
- hdr_out[ix].tracl = is*nx+ix+1;
- hdr_out[ix].tracf = ix+1;
- hdr_out[ix].scalco = -1000;
- hdr_out[ix].scalel = -1000;
- hdr_out[ix].sdepth = hdr_in[0].sdepth;
- hdr_out[ix].trid = 1;
- hdr_out[ix].ns = nshots;
- hdr_out[ix].trwf = nx;
- hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
- hdr_out[ix].f1 = -((float)(hdr_in[0].dt/1E6))*(nshots/2);
- hdr_out[ix].f2 = hdr_in[0].f2;
- hdr_out[ix].dt = hdr_in[0].dt;
- hdr_out[ix].d1 = ((float)hdr_in[0].dt);
- hdr_out[ix].d2 = (hdr_in[0].d2);
- hdr_out[ix].sx = (int)roundf(hdr_out[ix].f2 + (ix*hdr_out[ix].d2));
- hdr_out[ix].sx = hdr_in[ix].sx;
- hdr_out[ix].gx = (int)roundf(hdr_out[ix].f2 + (ix*hdr_out[ix].d2));
- hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ for (i1 = 0; i1 < n1; i1++) {
+ for (i3 = 0; i3 < n3; i3++) {
+ for (i2 = 0; i2 < n2; i2++) {
+ hdr_out[i3*n2+i2].fldr = i1+1;
+ hdr_out[i3*n2+i2].tracl = i1*n3*n2+i3*n2+i2+1;
+ hdr_out[i3*n2+i2].tracf = i3*n2+i2+1;
+ hdr_out[i3*n2+i2].scalco = -1000;
+ hdr_out[i3*n2+i2].scalel = -1000;
+ hdr_out[i3*n2+i2].sdepth = hdr_in[0].sdepth;
+ hdr_out[i3*n2+i2].ns = n4;
+ hdr_out[i3*n2+i2].trwf = n2*n3;
+ hdr_out[i3*n2+i2].ntr = hdr_out[i3*n2+i2].fldr*hdr_out[i3*n2+i2].trwf;
+ hdr_out[i3*n2+i2].f1 = f4;
+ hdr_out[i3*n2+i2].f2 = f2;
+ hdr_out[i3*n2+i2].dt = hdr_in[0].dt;
+ hdr_out[i3*n2+i2].d1 = d4;
+ hdr_out[i3*n2+i2].d2 = d2;
+ hdr_out[i3*n2+i2].sx = hdr_in[i3*n2+i2].sx;
+ hdr_out[i3*n2+i2].gx = hdr_in[i3*n2+i2].gx;
+ hdr_out[i3*n2+i2].sy = hdr_in[i3*n2+i2].sy;
+ hdr_out[i3*n2+i2].gy = hdr_in[i3*n2+i2].gy;
+ hdr_out[i3*n2+i2].offset = hdr_in[i3*n2+i2].offset;
+ if (hdr.trid==2) {
+ hdr_out[i3*n2+i2].ungpow = d1;
+ hdr_out[i3*n2+i2].unscale = f1;
+ hdr_out[i3*n2+i2].trid = 3;
+ }
+ else {
+ hdr_out[i3*n2+i2].trid = 2;
+ }
+ }
}
- ret = writeData(fp_out, &outdata[is*nx*nshots], hdr_out, nshots, nx);
+ ret = writeData3D(fp_out, &outdata[i1*n2*n3*n4], hdr_out, n4, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
}
+ free(outdata);free(hdr_in);free(hdr_out);
+ if (verbose) vmess("Wrote data");
fclose(fp_out);
return 0;
diff --git a/utils/snap2shot.c b/utils/snap2shot.c
index 227a5762a58e89f7aab7e721618f76348e670c88..66a53ef4e16acbd1e4b4672ea9caa1378d5279ee 100644
--- a/utils/snap2shot.c
+++ b/utils/snap2shot.c
@@ -21,8 +21,8 @@ typedef struct _complexStruct { /* complex number */
} 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);
+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);
@@ -30,7 +30,8 @@ long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
void name_ext(char *filename, char *extension);
-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);
+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);
char *sdoc[] = {
" ",
@@ -41,17 +42,17 @@ char *sdoc[] = {
" ",
" Required parameters: ",
"",
-" file_snap= ............... File containing the snapshot data that will be muted",
-" file_rcv= ................ File containing the snapshot data that will determine the mute window",
+" file_snap= ............... First filename that needs to be reshaped to receiver data",
" ",
" Optional parameters: ",
" ",
-" fout=out.su .............. Filename of the output",
-" shift=5 .................. Shift from the maximum",
-" smooth=5 ................. Length of smoothing taper",
-" mode=0 ................... Determine first arrival by maximum (mode=0), first event above tol (mode=1) or by raytime (mode=2)",
-" tol=1 .................... Tolerance for the determination of first arrival if mode=1",
-" fray ..................... File containing the raytimes of the first arrivals",
+" file_rcv= ................ File containing the snapshot data that will determine the mute window",
+" numb=0 ................... Starting number in the first file",
+" dnumb=1 .................. Increment per file name",
+" nxmax=0 .................. Maximum number of files that are to be reshaped",
+" numb=0 ................... Starting number in the first file",
+" nzstart=0 ................ Starting depth number to be written out",
+" nzend=end ................ Final depth number to be written out",
NULL};
int main (int argc, char **argv)
@@ -69,29 +70,34 @@ int main (int argc, char **argv)
initargs(argc, argv);
requestdoc(1);
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
if (!getparstring("file_rcv", &file_rcv)) file_rcv = "rcv.su";
if (!getparstring("file_snap", &file_snap)) file_snap = NULL;
if (file_snap == NULL) verr("No input data for the snap file given");
if (!getparlong("numb", &numb)) numb=0;
- if (!getparlong("dnumb", &dnumb)) dnumb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=1;
if (!getparlong("nxmax", &nxmax)) nxmax=0;
if (!getparlong("nzstart", &nzstart)) nzstart=0;
if (!getparlong("nzend", &nzend)) nzend=0;
+ /*----------------------------------------------------------------------------*
+ * Split the filename so the number can be changed
+ *----------------------------------------------------------------------------*/
if (dnumb < 1) dnumb = 1;
-
sprintf(numb1,"%li",numb);
-
ptr = strstr(file_snap,numb1);
pos = ptr - file_snap + 1;
-
pf = pos-1;
sprintf(fbegin,"%*.*s", pf, pf, file_snap);
sprintf(fend,"%s", file_snap+pos);
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files to be read
+ *----------------------------------------------------------------------------*/
file_det = 1;
nxr=0;
-
while (file_det) {
sprintf(fins,"%li",nxr*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
@@ -120,6 +126,9 @@ int main (int argc, char **argv)
}
}
+ /*----------------------------------------------------------------------------*
+ * Get the info from the files
+ *----------------------------------------------------------------------------*/
sprintf(fins,"%li",numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
@@ -131,12 +140,18 @@ int main (int argc, char **argv)
dnz = nzend-nzstart;
+ /*----------------------------------------------------------------------------*
+ * Allocate the data
+ *----------------------------------------------------------------------------*/
rcvdata = (float *)malloc(nxr*dnz*nxs*nys*nts*sizeof(float));
snapdata = (float *)malloc(dnz*nxs*nys*nts*sizeof(float));
hdr_snap = (segy *)calloc(nxs*nys*nts,sizeof(segy));
sx = (int *)malloc(nxr*sizeof(int));
sy = (int *)malloc(nxr*sizeof(int));
+ /*----------------------------------------------------------------------------*
+ * Reshape the data
+ *----------------------------------------------------------------------------*/
for (ixr=0; ixr<nxr; ixr++) {
vmess("Reshaping %li out of %li files",ixr+1,nxr);
sprintf(fins,"%li",ixr*dnumb+numb);
@@ -155,11 +170,12 @@ int main (int argc, char **argv)
}
}
}
-
free(snapdata);
+ /*----------------------------------------------------------------------------*
+ * Write out the data to new files
+ *----------------------------------------------------------------------------*/
hdr_rcv = (segy *)calloc(nxs*nys*nxr,sizeof(segy));
-
for (iz=nzstart; iz<nzend; iz++) {
vmess("Writing depth %li out of %li",iz-nzstart+1,dnz);
strcpy(file_tmp, file_rcv);
@@ -193,7 +209,6 @@ int main (int argc, char **argv)
}
}
}
-
// Write out homogeneous Green's function
ret = writeData3D(fp_rcv, (float *)&rcvdata[(iz-nzstart)*nys*nxs*nts*nxr], hdr_rcv, nts, nxr*nxs*nys);
if (ret < 0 ) verr("error on writing output file.");