From df87bff00484c7b31707e8ccefad6a00628dc08e Mon Sep 17 00:00:00 2001
From: Jan Thorbecke <jan@cray.com>
Date: Tue, 24 Jan 2017 14:17:29 +0100
Subject: [PATCH] cosmetic changes and marchenko code cleaning
---
FFTlib/Makefile | 2 +-
fdelmodc/defineSource.c | 442 +++++++++++------------
fdelmodc/getModelInfo.c | 96 ++---
fdelmodc/getWaveletInfo.c | 22 +-
fdelmodc/replacetab.scr | 3 +
fdelmodc/writeRec.c | 314 ++++++++--------
marchenko/applyMute.c | 2 +-
marchenko/demo/oneD/README | 7 +
marchenko/demo/oneD/epsMarchenkoIter.scr | 8 +
marchenko/demo/oneD/initialFocus.scr | 2 +-
marchenko/fmute.c | 19 +-
marchenko/marchenko.c | 163 +++------
utils/ftr1d.c | 2 +-
utils/syn2d.c | 13 +
14 files changed, 542 insertions(+), 553 deletions(-)
create mode 100755 fdelmodc/replacetab.scr
diff --git a/FFTlib/Makefile b/FFTlib/Makefile
index f239468..7585a21 100644
--- a/FFTlib/Makefile
+++ b/FFTlib/Makefile
@@ -43,7 +43,7 @@ $(LIB) : $(ARCH)
mkdirs:
-mkdir -p $(ROOT)/lib
-mkdir -p $(ROOT)/include
- -ln -sf $(ROOT)/FFTlib/genfft.h $I/genfft.h
+ -ln -sf genfft.h ../include/genfft.h
bins:
cd test ; $(MAKE)
diff --git a/fdelmodc/defineSource.c b/fdelmodc/defineSource.c
index e120b05..27d0f43 100644
--- a/fdelmodc/defineSource.c
+++ b/fdelmodc/defineSource.c
@@ -57,71 +57,71 @@ int defineSource(wavPar wav, srcPar src, float **src_nwav, int reverse, int verb
FILE *fp;
size_t nread;
int optn, nfreq, i, j, k, iwmax, tracesToDo;
- int iw, n1, namp;
+ int iw, n1, namp;
float scl, d1, df, deltom, om, tshift;
- float amp1, amp2, amp3;
- float *trace, maxampl;
+ float amp1, amp2, amp3;
+ float *trace, maxampl;
complex *ctrace, tmp;
segy hdr;
- n1 = wav.nt;
- if (wav.random) { /* initialize random sequence */
- srand48(wav.seed+1);
- seedCMWC4096();
- for (i=0; i<8192; i++) {
- amp1 = dcmwc4096();
- }
+ n1 = wav.nt;
+ if (wav.random) { /* initialize random sequence */
+ srand48(wav.seed+1);
+ seedCMWC4096();
+ for (i=0; i<8192; i++) {
+ amp1 = dcmwc4096();
+ }
- }
- else {
+ }
+ else {
/* read first header and last byte to get file size */
- fp = fopen( wav.file_src, "r" );
- assert( fp != NULL);
- nread = fread( &hdr, 1, TRCBYTES, fp );
- assert(nread == TRCBYTES);
-
+ fp = fopen( wav.file_src, "r" );
+ assert( fp != NULL);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+
/* read all traces */
- tracesToDo = wav.nx;
- i = 0;
- while (tracesToDo) {
- memset(&src_nwav[i][0],0,n1*sizeof(float));
- nread = fread(&src_nwav[i][0], sizeof(float), hdr.ns, fp);
- assert (nread == hdr.ns);
-
- nread = fread( &hdr, 1, TRCBYTES, fp );
- if (nread==0) break;
- tracesToDo--;
- i++;
- }
- fclose(fp);
- }
-
-
- optn = optncr(2*n1);
- nfreq = optn/2 + 1;
- ctrace = (complex *)calloc(nfreq,sizeof(complex));
- trace = (float *)calloc(optn,sizeof(float));
-
- df = 1.0/(optn*wav.dt);
+ tracesToDo = wav.nx;
+ i = 0;
+ while (tracesToDo) {
+ memset(&src_nwav[i][0],0,n1*sizeof(float));
+ nread = fread(&src_nwav[i][0], sizeof(float), hdr.ns, fp);
+ assert (nread == hdr.ns);
+
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread==0) break;
+ tracesToDo--;
+ i++;
+ }
+ fclose(fp);
+ }
+
+
+ optn = optncr(2*n1);
+ nfreq = optn/2 + 1;
+ ctrace = (complex *)calloc(nfreq,sizeof(complex));
+ trace = (float *)calloc(optn,sizeof(float));
+
+ df = 1.0/(optn*wav.dt);
deltom = 2.*M_PI*df;
- scl = 1.0/optn;
- maxampl=0.0;
- iwmax = nfreq;
-
- for (i=0; i<wav.nx; i++) {
- if (wav.random) {
- randomWavelet(wav, src, &src_nwav[i][0], src.tbeg[i], src.tend[i], verbose);
- }
- else {
- memset(&ctrace[0].r,0,nfreq*sizeof(complex));
- memset(&trace[0],0,optn*sizeof(float));
- memcpy(&trace[0],&src_nwav[i][0],n1*sizeof(float));
- rc1fft(trace,ctrace,optn,-1);
- /* Scale source from file with -j/w (=1/(jw)) for volume source injections
- no scaling is applied for volume source injection rates */
+ scl = 1.0/optn;
+ maxampl=0.0;
+ iwmax = nfreq;
+
+ for (i=0; i<wav.nx; i++) {
+ if (wav.random) {
+ randomWavelet(wav, src, &src_nwav[i][0], src.tbeg[i], src.tend[i], verbose);
+ }
+ else {
+ memset(&ctrace[0].r,0,nfreq*sizeof(complex));
+ memset(&trace[0],0,optn*sizeof(float));
+ memcpy(&trace[0],&src_nwav[i][0],n1*sizeof(float));
+ rc1fft(trace,ctrace,optn,-1);
+ /* Scale source from file with -j/w (=1/(jw)) for volume source injections
+ no scaling is applied for volume source injection rates */
if (src.injectionrate==0) {
for (iw=1;iw<iwmax;iw++) {
om = 1.0/(deltom*iw);
@@ -134,99 +134,99 @@ int defineSource(wavPar wav, srcPar src, float **src_nwav, int reverse, int verb
/*
*/
if (src.type < 6) { // shift wavelet with +1/2 DeltaT due to staggered in time
- tshift=0.5*wav.dt;
+ tshift=0.5*wav.dt;
for (iw=1;iw<iwmax;iw++) {
- om = deltom*iw*tshift;
- tmp.r = ctrace[iw].r*cos(-om) - ctrace[iw].i*sin(-om);
- tmp.i = ctrace[iw].i*cos(-om) + ctrace[iw].r*sin(-om);
+ om = deltom*iw*tshift;
+ tmp.r = ctrace[iw].r*cos(-om) - ctrace[iw].i*sin(-om);
+ tmp.i = ctrace[iw].i*cos(-om) + ctrace[iw].r*sin(-om);
ctrace[iw].r = tmp.r;
ctrace[iw].i = tmp.i;
}
- }
-
- /* zero frequency iw=0 set to 0 if the next sample has amplitude==0*/
- amp1 = sqrt(ctrace[1].r*ctrace[1].r+ctrace[1].i*ctrace[1].i);
- if (amp1 == 0.0) {
- ctrace[0].r = ctrace[0].i = 0.0;
- }
- else { /* stabilization for w=0: extrapolate amplitudes to 0 */
- amp2 = sqrt(ctrace[2].r*ctrace[2].r+ctrace[2].i*ctrace[2].i);
- amp3 = sqrt(ctrace[3].r*ctrace[3].r+ctrace[3].i*ctrace[3].i);
- ctrace[0].r = amp1+(2.0*(amp1-amp2)-(amp2-amp3));
- ctrace[0].i = 0.0;
- if (ctrace[1].r < 0.0) {
- ctrace[0].r *= -1.0;
- }
- }
- for (iw=iwmax;iw<nfreq;iw++) {
- ctrace[iw].r = 0.0;
- ctrace[iw].i = 0.0;
- }
- cr1fft(ctrace,trace,optn,1);
- /* avoid a (small) spike in the last sample
- this is done to avoid diffraction from last wavelet sample
- which will act as a pulse */
- if (reverse) {
- for (j=0; j<n1; j++) src_nwav[i][j] = scl*(trace[n1-j-1]-trace[0]);
-// for (j=0; j<n1; j++) src_nwav[i][j] = scl*(trace[j]-trace[optn-1]);
- }
- else {
- for (j=0; j<n1; j++) src_nwav[i][j] = scl*(trace[j]-trace[optn-1]);
- }
- }
+ }
+
+ /* zero frequency iw=0 set to 0 if the next sample has amplitude==0*/
+ amp1 = sqrt(ctrace[1].r*ctrace[1].r+ctrace[1].i*ctrace[1].i);
+ if (amp1 == 0.0) {
+ ctrace[0].r = ctrace[0].i = 0.0;
+ }
+ else { /* stabilization for w=0: extrapolate amplitudes to 0 */
+ amp2 = sqrt(ctrace[2].r*ctrace[2].r+ctrace[2].i*ctrace[2].i);
+ amp3 = sqrt(ctrace[3].r*ctrace[3].r+ctrace[3].i*ctrace[3].i);
+ ctrace[0].r = amp1+(2.0*(amp1-amp2)-(amp2-amp3));
+ ctrace[0].i = 0.0;
+ if (ctrace[1].r < 0.0) {
+ ctrace[0].r *= -1.0;
+ }
+ }
+ for (iw=iwmax;iw<nfreq;iw++) {
+ ctrace[iw].r = 0.0;
+ ctrace[iw].i = 0.0;
+ }
+ cr1fft(ctrace,trace,optn,1);
+ /* avoid a (small) spike in the last sample
+ this is done to avoid diffraction from last wavelet sample
+ which will act as a pulse */
+ if (reverse) {
+ for (j=0; j<n1; j++) src_nwav[i][j] = scl*(trace[n1-j-1]-trace[0]);
+// for (j=0; j<n1; j++) src_nwav[i][j] = scl*(trace[j]-trace[optn-1]);
+ }
+ else {
+ for (j=0; j<n1; j++) src_nwav[i][j] = scl*(trace[j]-trace[optn-1]);
+ }
+ }
}
free(ctrace);
free(trace);
/* use random amplitude gain factor for each source */
- if (src.amplitude > 0.0) {
- namp=wav.nx*10;
- trace = (float *)calloc(2*namp,sizeof(float));
- for (i=0; i<wav.nx; i++) {
- if (src.distribution) {
- scl = gaussGen()*src.amplitude;
- k = (int)MAX(MIN(namp*(scl+5*src.amplitude)/(10*src.amplitude),namp-1),0);
- d1 = 10.0*src.amplitude/(namp-1);
- }
- else {
- scl = (float)(drand48()-0.5)*src.amplitude;
- k = (int)MAX(MIN(namp*(scl+1*src.amplitude)/(2*src.amplitude),namp-1),0);
- d1 = 2.0*src.amplitude/(namp-1);
- }
-
- trace[k] += 1.0;
-/* trace[i] = scl; */
- if (wav.random) n1 = wav.nsamp[i];
- else n1 = wav.nt;
- for (j=0; j<n1; j++) {
- src_nwav[i][j] *= scl;
- }
- }
- if (verbose>2) writesufile("src_ampl.su", trace, namp, 1, -5*src.amplitude, 0.0, d1, 1);
+ if (src.amplitude > 0.0) {
+ namp=wav.nx*10;
+ trace = (float *)calloc(2*namp,sizeof(float));
+ for (i=0; i<wav.nx; i++) {
+ if (src.distribution) {
+ scl = gaussGen()*src.amplitude;
+ k = (int)MAX(MIN(namp*(scl+5*src.amplitude)/(10*src.amplitude),namp-1),0);
+ d1 = 10.0*src.amplitude/(namp-1);
+ }
+ else {
+ scl = (float)(drand48()-0.5)*src.amplitude;
+ k = (int)MAX(MIN(namp*(scl+1*src.amplitude)/(2*src.amplitude),namp-1),0);
+ d1 = 2.0*src.amplitude/(namp-1);
+ }
+
+ trace[k] += 1.0;
+/* trace[i] = scl; */
+ if (wav.random) n1 = wav.nsamp[i];
+ else n1 = wav.nt;
+ for (j=0; j<n1; j++) {
+ src_nwav[i][j] *= scl;
+ }
+ }
+ if (verbose>2) writesufile("src_ampl.su", trace, namp, 1, -5*src.amplitude, 0.0, d1, 1);
/*
- qsort(trace,wav.nx,sizeof(float), comp);
- for (i=0; i<wav.nx; i++) {
- scl = trace[i];
- trace[i] = normal(scl, 0.0, src.amplitude);
- }
- if (verbose>2) writesufile("src_ampl.su", trace, wav.nx, 1, -5*src.amplitude, 0.0, d1, 1);
+ qsort(trace,wav.nx,sizeof(float), comp);
+ for (i=0; i<wav.nx; i++) {
+ scl = trace[i];
+ trace[i] = normal(scl, 0.0, src.amplitude);
+ }
+ if (verbose>2) writesufile("src_ampl.su", trace, wav.nx, 1, -5*src.amplitude, 0.0, d1, 1);
*/
- free(trace);
- }
+ free(trace);
+ }
- if (verbose>3) writesufilesrcnwav("src_nwav.su", src_nwav, wav, wav.nt, wav.nx, 0.0, 0.0, wav.dt, 1);
+ if (verbose>3) writesufilesrcnwav("src_nwav.su", src_nwav, wav, wav.nt, wav.nx, 0.0, 0.0, wav.dt, 1);
/* set maximum amplitude in source file to 1.0 */
/*
- assert(maxampl != 0.0);
- scl = wav.dt/(maxampl);
- scl = 1.0/(maxampl);
- for (i=0; i<wav.nx; i++) {
- for (j=0; j<n1; j++) {
- src_nwav[i*n1+j] *= scl;
- }
+ assert(maxampl != 0.0);
+ scl = wav.dt/(maxampl);
+ scl = 1.0/(maxampl);
+ for (i=0; i<wav.nx; i++) {
+ for (j=0; j<n1; j++) {
+ src_nwav[i*n1+j] *= scl;
+ }
}
*/
@@ -237,93 +237,93 @@ int defineSource(wavPar wav, srcPar src, float **src_nwav, int reverse, int verb
int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend, int verbose)
{
int optn, nfreq, j, iwmax;
- int iw, n1, itbeg, itmax, nsmth;
+ int iw, n1, itbeg, itmax, nsmth;
float df, amp1;
- float *rtrace;
- float x1, x2, z1, z2, dzdx1, dzdx2, a, b, c, d, t;
+ float *rtrace;
+ float x1, x2, z1, z2, dzdx1, dzdx2, a, b, c, d, t;
complex *ctrace;
- n1 = wav.nt; /* this is set to the maximum length (tlength/dt) */
-
- optn = optncr(2*n1);
- nfreq = optn/2 + 1;
- ctrace = (complex *)calloc(nfreq,sizeof(complex));
- rtrace = (float *)calloc(optn,sizeof(float));
-
- df = 1.0/(optn*wav.dt);
+ n1 = wav.nt; /* this is set to the maximum length (tlength/dt) */
+
+ optn = optncr(2*n1);
+ nfreq = optn/2 + 1;
+ ctrace = (complex *)calloc(nfreq,sizeof(complex));
+ rtrace = (float *)calloc(optn,sizeof(float));
+
+ df = 1.0/(optn*wav.dt);
+
+ iwmax = MIN(NINT(wav.fmax/df),nfreq);
+
+ for (iw=1;iw<iwmax;iw++) {
+ ctrace[iw].r = (float)(drand48()-0.5);
+ ctrace[iw].i = (float)(drand48()-0.5);
+ }
+ for (iw=iwmax;iw<nfreq;iw++) {
+ ctrace[iw].r = 0.0;
+ ctrace[iw].i = 0.0;
+ }
+ cr1fft(ctrace,rtrace,optn,1);
+
+ /* find first zero crossing in wavelet */
+ amp1 = rtrace[0];
+ j = 1;
+ if (amp1 < 0.0) {
+ while (rtrace[j] < 0.0) j++;
+ }
+ else {
+ while (rtrace[j] > 0.0) j++;
+ }
+ itbeg = j;
+
+ /* find last zero crossing in wavelet */
+// itmax = itbeg+MIN(NINT((tend-tbeg)/wav.dt),n1);
+ itmax = MIN(NINT(itbeg+(tend-tbeg)/wav.dt),n1);
+
+ amp1 = rtrace[itmax-1];
+ j = itmax;
+ if (amp1 < 0.0) {
+ while (rtrace[j] < 0.0 && j>itbeg) j--;
+ }
+ else {
+ while (rtrace[j] > 0.0 && j>itbeg) j--;
+ }
+ itmax = j;
+
+ /* make smooth transitions to zero aamplitude */
+ nsmth=MIN(10,itmax);
+ x1 = 0.0;
+ z1 = 0.0;
+ dzdx1 = 0.0;
+ x2 = nsmth;
+ z2 = rtrace[itbeg+nsmth];
+// dzdx2 = (rtrace[itbeg+(nsmth+1)]-rtrace[itbeg+(nsmth-1)])/(2.0);
+ dzdx2 = (rtrace[itbeg+nsmth-2]-8.0*rtrace[itbeg+nsmth-1]+
+ 8.0*rtrace[itbeg+nsmth+1]-rtrace[itbeg+nsmth+2])/(12.0);
+ spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
+ for (j=0; j<nsmth; j++) {
+ t = j;
+ rtrace[itbeg+j] = a*t*t*t+b*t*t+c*t+d;
+ }
+
+ x1 = 0.0;
+ z1 = rtrace[itmax-nsmth];
+// dzdx1 = (rtrace[itmax-(nsmth-1)]-rtrace[itmax-(nsmth+1)])/(2.0);
+ dzdx1 = (rtrace[itmax-nsmth-2]-8.0*rtrace[itmax-nsmth-1]+
+ 8.0*rtrace[itmax-nsmth+1]-rtrace[itmax-nsmth+2])/(12.0);
+ x2 = nsmth;
+ z2 = 0.0;
+ dzdx2 = 0.0;
+
+// fprintf(stderr,"x1=%f z1=%f d=%f x2=%f, z2=%f d=%f\n",x1,z1,dzdx1,x2,z2,dzdx2);
+ spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
+ for (j=0; j<nsmth; j++) {
+ t = j;
+ rtrace[itmax-nsmth+j] = a*t*t*t+b*t*t+c*t+d;
+// fprintf(stderr,"a=%f b=%f c=%f d=%f rtrace%d=%f \n",a,b,c,d,j,rtrace[itmax-nsmth+j]);
+ }
+
+ for (j=itbeg; j<itmax; j++) trace[j-itbeg] = rtrace[j];
- iwmax = MIN(NINT(wav.fmax/df),nfreq);
-
- for (iw=1;iw<iwmax;iw++) {
- ctrace[iw].r = (float)(drand48()-0.5);
- ctrace[iw].i = (float)(drand48()-0.5);
- }
- for (iw=iwmax;iw<nfreq;iw++) {
- ctrace[iw].r = 0.0;
- ctrace[iw].i = 0.0;
- }
- cr1fft(ctrace,rtrace,optn,1);
-
- /* find first zero crossing in wavelet */
- amp1 = rtrace[0];
- j = 1;
- if (amp1 < 0.0) {
- while (rtrace[j] < 0.0) j++;
- }
- else {
- while (rtrace[j] > 0.0) j++;
- }
- itbeg = j;
-
- /* find last zero crossing in wavelet */
-// itmax = itbeg+MIN(NINT((tend-tbeg)/wav.dt),n1);
- itmax = MIN(NINT(itbeg+(tend-tbeg)/wav.dt),n1);
-
- amp1 = rtrace[itmax-1];
- j = itmax;
- if (amp1 < 0.0) {
- while (rtrace[j] < 0.0 && j>itbeg) j--;
- }
- else {
- while (rtrace[j] > 0.0 && j>itbeg) j--;
- }
- itmax = j;
-
- /* make smooth transitions to zero aamplitude */
- nsmth=MIN(10,itmax);
- x1 = 0.0;
- z1 = 0.0;
- dzdx1 = 0.0;
- x2 = nsmth;
- z2 = rtrace[itbeg+nsmth];
-// dzdx2 = (rtrace[itbeg+(nsmth+1)]-rtrace[itbeg+(nsmth-1)])/(2.0);
- dzdx2 = (rtrace[itbeg+nsmth-2]-8.0*rtrace[itbeg+nsmth-1]+
- 8.0*rtrace[itbeg+nsmth+1]-rtrace[itbeg+nsmth+2])/(12.0);
- spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
- for (j=0; j<nsmth; j++) {
- t = j;
- rtrace[itbeg+j] = a*t*t*t+b*t*t+c*t+d;
- }
-
- x1 = 0.0;
- z1 = rtrace[itmax-nsmth];
-// dzdx1 = (rtrace[itmax-(nsmth-1)]-rtrace[itmax-(nsmth+1)])/(2.0);
- dzdx1 = (rtrace[itmax-nsmth-2]-8.0*rtrace[itmax-nsmth-1]+
- 8.0*rtrace[itmax-nsmth+1]-rtrace[itmax-nsmth+2])/(12.0);
- x2 = nsmth;
- z2 = 0.0;
- dzdx2 = 0.0;
-
-// fprintf(stderr,"x1=%f z1=%f d=%f x2=%f, z2=%f d=%f\n",x1,z1,dzdx1,x2,z2,dzdx2);
- spline3(x1, x2, z1, z2, dzdx1, dzdx2, &a, &b, &c, &d);
- for (j=0; j<nsmth; j++) {
- t = j;
- rtrace[itmax-nsmth+j] = a*t*t*t+b*t*t+c*t+d;
-// fprintf(stderr,"a=%f b=%f c=%f d=%f rtrace%d=%f \n",a,b,c,d,j,rtrace[itmax-nsmth+j]);
- }
-
- for (j=itbeg; j<itmax; j++) trace[j-itbeg] = rtrace[j];
-
free(ctrace);
free(rtrace);
@@ -332,16 +332,16 @@ int randomWavelet(wavPar wav, srcPar src, float *trace, float tbeg, float tend,
float normal(double x,double mu,double sigma)
{
- return (float)(1.0/(2.0*M_PI*sigma*sigma))*exp(-1.0*(((x-mu)*(x-mu))/(2.0*sigma*sigma)) );
+ return (float)(1.0/(2.0*M_PI*sigma*sigma))*exp(-1.0*(((x-mu)*(x-mu))/(2.0*sigma*sigma)) );
}
int comp (const float *a, const float *b)
{
- if (*a==*b)
- return 0;
- else
- if (*a < *b)
- return -1;
- else
- return 1;
+ if (*a==*b)
+ return 0;
+ else
+ if (*a < *b)
+ return -1;
+ else
+ return 1;
}
diff --git a/fdelmodc/getModelInfo.c b/fdelmodc/getModelInfo.c
index 20ed279..378a1b5 100644
--- a/fdelmodc/getModelInfo.c
+++ b/fdelmodc/getModelInfo.c
@@ -25,85 +25,85 @@
int getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose)
{
FILE *fp;
- size_t nread;
- off_t bytes, ret, trace_sz, ntraces;
- int i, one_shot;
- float *trace, cmin;
- segy hdr;
+ size_t nread, trace_sz;
+ off_t bytes;
+ int ret, i, one_shot, ntraces;
+ float *trace, cmin;
+ segy hdr;
fp = fopen( file_name, "r" );
assert( fp != NULL);
nread = fread( &hdr, 1, TRCBYTES, fp );
assert(nread == TRCBYTES);
ret = fseeko( fp, 0, SEEK_END );
- if (ret<0) perror("fseeko");
+ if (ret<0) perror("fseeko");
bytes = ftello( fp );
- *n1 = hdr.ns;
- *d1 = hdr.d1;
- *d2 = hdr.d2;
- *f1 = hdr.f1;
- *f2 = hdr.f2;
+ *n1 = hdr.ns;
+ *d1 = hdr.d1;
+ *d2 = hdr.d2;
+ *f1 = hdr.f1;
+ *f2 = hdr.f2;
- if ( NINT(100.0*((*d1)/(*d2)))!=100 ) {
- verr("dx and dz are different in the model !");
- }
- if ( NINT(1000.0*(*d1))==0 ) {
- if(!getparfloat("dx",d1)) {
- verr("dx is equal to zero use parameter dx= to set value");
- }
- *d2 = *d1;
- }
+ if ( NINT(100.0*((*d1)/(*d2)))!=100 ) {
+ verr("dx and dz are different in the model !");
+ }
+ if ( NINT(1000.0*(*d1))==0 ) {
+ if(!getparfloat("dx",d1)) {
+ verr("dx is equal to zero use parameter dx= to set value");
+ }
+ *d2 = *d1;
+ }
trace_sz = sizeof(float)*(*n1)+TRCBYTES;
ntraces = (int) (bytes/trace_sz);
- *n2 = ntraces;
+ *n2 = ntraces;
/* check to find out min and max values gather */
one_shot = 1;
trace = (float *)malloc(trace_sz);
fseeko( fp, TRCBYTES, SEEK_SET );
- nread = fread( trace, sizeof(float), hdr.ns, fp );
- assert (nread == hdr.ns);
+ nread = fread( trace, sizeof(float), hdr.ns, fp );
+ assert (nread == hdr.ns);
fseeko( fp, TRCBYTES, SEEK_SET );
- if (hdr.trid == TRID_DEPTH) *axis = 1; /* samples are z-axis */
- else *axis = 0; /* sample direction respresents the x-axis */
+ if (hdr.trid == TRID_DEPTH) *axis = 1; /* samples are z-axis */
+ else *axis = 0; /* sample direction respresents the x-axis */
- i=0; cmin=trace[0];
- while ( ( (cmin==0.0) && zeroch) && (i<hdr.ns) ) cmin=trace[i++];
+ i=0; cmin=trace[0];
+ while ( ( (cmin==0.0) && zeroch) && (i<hdr.ns) ) cmin=trace[i++];
- *max = cmin;
- *min = cmin;
- /* keep on reading traces until there are no more traces (nread==0) */
+ *max = cmin;
+ *min = cmin;
+ /* keep on reading traces until there are no more traces (nread==0) */
while (one_shot) {
nread = fread( trace, sizeof(float), hdr.ns, fp );
assert (nread == hdr.ns);
- for (i=0;i<hdr.ns;i++) {
- *max = MAX(trace[i],*max);
- cmin = MIN(trace[i],*min);
+ for (i=0;i<(*n1);i++) {
+ *max = MAX(trace[i],*max);
+ cmin = MIN(trace[i],*min);
if (zeroch) {
- if (cmin!=0.0) *min = MIN(*min, cmin);
- }
- else {
- *min = cmin;
- }
- }
+ if (cmin!=0.0) *min = MIN(*min, cmin);
+ }
+ else {
+ *min = cmin;
+ }
+ }
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread==0) break;
}
fclose(fp);
free(trace);
- if (verbose>2) {
- vmess("For file %s", file_name);
- vmess("nz=%d nx=%d", *n1, *n2);
- vmess("dz=%f dx=%f", *d1, *d2);
- vmess("min=%f max=%f", *min, *max);
- vmess("zstart=%f xstart=%f", *f1, *f2);
- if (*axis) vmess("sample represent z-axis\n");
- else vmess("sample represent x-axis\n");
- }
+ if (verbose>2) {
+ vmess("For file %s", file_name);
+ vmess("nz=%d nx=%d", *n1, *n2);
+ vmess("dz=%f dx=%f", *d1, *d2);
+ vmess("min=%f max=%f", *min, *max);
+ vmess("zstart=%f xstart=%f", *f1, *f2);
+ if (*axis) vmess("sample represent z-axis\n");
+ else vmess("sample represent x-axis\n");
+ }
return 0;
}
diff --git a/fdelmodc/getWaveletInfo.c b/fdelmodc/getWaveletInfo.c
index bcf2a77..4fcfbc0 100644
--- a/fdelmodc/getWaveletInfo.c
+++ b/fdelmodc/getWaveletInfo.c
@@ -37,12 +37,12 @@ void rc1fft(float *rdata, complex *cdata, int n, int sign);
int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *fmax, int *nxm, int verbose)
{
FILE *fp;
- size_t nread;
- off_t bytes, ret, trace_sz, ntraces;
- int one_shot;
- int optn, nfreq, i, iwmax;
- float *trace;
- float ampl, amplmax, tampl, tamplmax;
+ size_t nread, trace_sz;
+ off_t bytes;
+ int ret, one_shot, ntraces;
+ int optn, nfreq, i, iwmax;
+ float *trace;
+ float ampl, amplmax, tampl, tamplmax;
complex *ctrace;
segy hdr;
@@ -67,13 +67,13 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
}
*f2 = hdr.f2;
- trace_sz = sizeof(float)*(*n1)+TRCBYTES;
+ trace_sz = (size_t)(sizeof(float)*(*n1)+TRCBYTES);
ntraces = (int) (bytes/trace_sz);
*n2 = ntraces;
/* check to find out number of traces in shot gather */
- optn = optncr(hdr.ns);
+ optn = optncr(*n1);
nfreq = optn/2 + 1;
ctrace = (complex *)malloc(nfreq*sizeof(complex));
one_shot = 1;
@@ -82,10 +82,10 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
while (one_shot) {
memset(trace,0,optn*sizeof(float));
- nread = fread( trace, sizeof(float), hdr.ns, fp );
- assert (nread == hdr.ns);
+ nread = fread( trace, sizeof(float), *n1, fp );
+ assert (nread == *n1);
tamplmax = 0.0;
- for (i=0;i<hdr.ns;i++) {
+ for (i=0;i<(*n1);i++) {
tampl = fabsf(trace[i]);
if (tampl > tamplmax) tamplmax = tampl;
}
diff --git a/fdelmodc/replacetab.scr b/fdelmodc/replacetab.scr
new file mode 100755
index 0000000..48fc0b5
--- /dev/null
+++ b/fdelmodc/replacetab.scr
@@ -0,0 +1,3 @@
+#!/bin/bash
+sed 's/ / /g' $1 > nep
+mv nep $1
diff --git a/fdelmodc/writeRec.c b/fdelmodc/writeRec.c
index e15ccaf..877bbe9 100644
--- a/fdelmodc/writeRec.c
+++ b/fdelmodc/writeRec.c
@@ -14,10 +14,10 @@
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
- float r,i;
+ float r,i;
} complex;
typedef struct _dcomplexStruct { /* complex number */
- double r,i;
+ double r,i;
} dcomplex;
#endif/* complex */
@@ -42,172 +42,172 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
int writeRec(recPar rec, modPar mod, bndPar bnd, wavPar wav, int ixsrc, int izsrc, int nsam, int ishot, int fileno,
- float *rec_vx, float *rec_vz, float *rec_txx, float *rec_tzz, float *rec_txz,
- float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose)
+ float *rec_vx, float *rec_vz, float *rec_txx, float *rec_tzz, float *rec_txz,
+ float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose)
{
FILE *fpvx, *fpvz, *fptxx, *fptzz, *fptxz, *fpp, *fppp, *fpss, *fpup, *fpdown;
- float *rec_up, *rec_down, *trace, *rec_vze, *rec_pe;
- float dx, dt, cp, rho, fmin, fmax;
- complex *crec_vz, *crec_p, *crec_up, *crec_dw;
+ float *rec_up, *rec_down, *trace, *rec_vze, *rec_pe;
+ float dx, dt, cp, rho, fmin, fmax;
+ complex *crec_vz, *crec_p, *crec_up, *crec_dw;
int irec, ntfft, nfreq, nkx, xorig, ix, iz, it, ibndx;
- int append;
- double ddt;
- char number[16], filename[1024];
+ int append;
+ double ddt;
+ char number[16], filename[1024];
segy hdr;
- if (!rec.n) return 0;
- if (ishot) append=1;
- else append=0;
-
- /* if the total number of samples exceeds rec_ntsam then a new (numbered) file is opened */
- /* fileno has a non-zero value (from fdelmodc.c) if the number of samples exceeds rec_ntsam. */
- strcpy(filename, rec.file_rcv);
- if (fileno) {
- sprintf(number,"_%03d",fileno);
- name_ext(filename, number);
- }
-
- if (verbose>2) vmess("Writing receiver data to file %s", filename);
- if (nsam != rec.nt && verbose) vmess("Number of samples written to last file = %d",nsam);
-
- memset(&hdr,0,TRCBYTES);
- ddt = (double)mod.dt;/* to avoid rounding in 32 bit precision */
+ if (!rec.n) return 0;
+ if (ishot) append=1;
+ else append=0;
+
+ /* if the total number of samples exceeds rec_ntsam then a new (numbered) file is opened */
+ /* fileno has a non-zero value (from fdelmodc.c) if the number of samples exceeds rec_ntsam. */
+ strcpy(filename, rec.file_rcv);
+ if (fileno) {
+ sprintf(number,"_%03d",fileno);
+ name_ext(filename, number);
+ }
+
+ if (verbose>2) vmess("Writing receiver data to file %s", filename);
+ if (nsam != rec.nt && verbose) vmess("Number of samples written to last file = %d",nsam);
+
+ memset(&hdr,0,TRCBYTES);
+ ddt = (double)mod.dt;/* to avoid rounding in 32 bit precision */
dt = (float)ddt*rec.skipdt;
- dx = (rec.x[1]-rec.x[0])*mod.dx;
- hdr.dt = (unsigned short)lround((((double)1.0e6*ddt*rec.skipdt)));
- hdr.scalco = -1000;
- hdr.scalel = -1000;
- hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
- hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
+ dx = (rec.x[1]-rec.x[0])*mod.dx;
+ hdr.dt = (unsigned short)lround((((double)1.0e6*ddt*rec.skipdt)));
+ hdr.scalco = -1000;
+ hdr.scalel = -1000;
+ hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
+ hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
hdr.selev = (int)(-1000.0*(mod.z0+izsrc*mod.dz));
- hdr.fldr = ishot+1;
- hdr.trid = 1;
- hdr.ns = nsam;
- hdr.trwf = rec.n;
- hdr.ntr = (ishot+1)*rec.n;
- hdr.f1 = 0.0;
- hdr.d1 = mod.dt*rec.skipdt;
- hdr.d2 = (rec.x[1]-rec.x[0])*mod.dx;
- hdr.f2 = mod.x0+rec.x[0]*mod.dx;
-
- if (rec.type.vx) fpvx = fileOpen(filename, "_rvx", append);
- if (rec.type.vz) fpvz = fileOpen(filename, "_rvz", append);
- if (rec.type.p) fpp = fileOpen(filename, "_rp", append);
- if (rec.type.txx) fptxx = fileOpen(filename, "_rtxx", append);
- if (rec.type.tzz) fptzz = fileOpen(filename, "_rtzz", append);
- if (rec.type.txz) fptxz = fileOpen(filename, "_rtxz", append);
- if (rec.type.pp) fppp = fileOpen(filename, "_rpp", append);
- if (rec.type.ss) fpss = fileOpen(filename, "_rss", append);
-
- /* decomposed wavefield */
- if (rec.type.ud && (mod.ischeme==1 || mod.ischeme==2) ) {
- fpup = fileOpen(filename, "_ru", append);
- fpdown = fileOpen(filename, "_rd", append);
- ntfft = optncr(nsam);
- nfreq = ntfft/2+1;
- fmin = 0.0;
- fmax = wav.fmax;
- nkx = optncc(2*mod.nax);
- ibndx = mod.ioPx;
- if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
- cp = rec.cp;
- rho = rec.rho;
- if (verbose) vmess("Decomposition array at z=%.2f with cp=%.2f rho=%.2f", rec.zr[0]+mod.z0, cp, rho);
- rec_up = (float *)calloc(ntfft*nkx,sizeof(float));
- rec_down= (float *)calloc(ntfft*nkx,sizeof(float));
- crec_vz = (complex *)malloc(nfreq*nkx*sizeof(complex));
- crec_p = (complex *)malloc(nfreq*nkx*sizeof(complex));
- crec_up = (complex *)malloc(nfreq*nkx*sizeof(complex));
- crec_dw = (complex *)malloc(nfreq*nkx*sizeof(complex));
-
- rec_vze = rec_up;
- rec_pe = rec_down;
- /* copy input data into extended arrays with padded zeroes */
- for (ix=0; ix<mod.nax; ix++) {
- memcpy(&rec_vze[ix*ntfft],&rec_udvz[ix*rec.nt],nsam*sizeof(float));
- memcpy(&rec_pe[ix*ntfft], &rec_udp[ix*rec.nt], nsam*sizeof(float));
- }
-
- /* transform from t-x to kx-w */
- xorig = ixsrc+ibndx;
- xt2wkx(rec_vze, crec_vz, ntfft, nkx, ntfft, nkx, xorig);
- xt2wkx(rec_pe, crec_p, ntfft, nkx, ntfft, nkx, xorig);
-
- /* apply decomposition operators */
- kxwdecomp(crec_p, crec_vz, crec_up, crec_dw,
+ hdr.fldr = ishot+1;
+ hdr.trid = 1;
+ hdr.ns = nsam;
+ hdr.trwf = rec.n;
+ hdr.ntr = (ishot+1)*rec.n;
+ hdr.f1 = 0.0;
+ hdr.d1 = mod.dt*rec.skipdt;
+ hdr.d2 = (rec.x[1]-rec.x[0])*mod.dx;
+ hdr.f2 = mod.x0+rec.x[0]*mod.dx;
+
+ if (rec.type.vx) fpvx = fileOpen(filename, "_rvx", append);
+ if (rec.type.vz) fpvz = fileOpen(filename, "_rvz", append);
+ if (rec.type.p) fpp = fileOpen(filename, "_rp", append);
+ if (rec.type.txx) fptxx = fileOpen(filename, "_rtxx", append);
+ if (rec.type.tzz) fptzz = fileOpen(filename, "_rtzz", append);
+ if (rec.type.txz) fptxz = fileOpen(filename, "_rtxz", append);
+ if (rec.type.pp) fppp = fileOpen(filename, "_rpp", append);
+ if (rec.type.ss) fpss = fileOpen(filename, "_rss", append);
+
+ /* decomposed wavefield */
+ if (rec.type.ud && (mod.ischeme==1 || mod.ischeme==2) ) {
+ fpup = fileOpen(filename, "_ru", append);
+ fpdown = fileOpen(filename, "_rd", append);
+ ntfft = optncr(nsam);
+ nfreq = ntfft/2+1;
+ fmin = 0.0;
+ fmax = wav.fmax;
+ nkx = optncc(2*mod.nax);
+ ibndx = mod.ioPx;
+ if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+ cp = rec.cp;
+ rho = rec.rho;
+ if (verbose) vmess("Decomposition array at z=%.2f with cp=%.2f rho=%.2f", rec.zr[0]+mod.z0, cp, rho);
+ rec_up = (float *)calloc(ntfft*nkx,sizeof(float));
+ rec_down= (float *)calloc(ntfft*nkx,sizeof(float));
+ crec_vz = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ crec_p = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ crec_up = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ crec_dw = (complex *)malloc(nfreq*nkx*sizeof(complex));
+
+ rec_vze = rec_up;
+ rec_pe = rec_down;
+ /* copy input data into extended arrays with padded zeroes */
+ for (ix=0; ix<mod.nax; ix++) {
+ memcpy(&rec_vze[ix*ntfft],&rec_udvz[ix*rec.nt],nsam*sizeof(float));
+ memcpy(&rec_pe[ix*ntfft], &rec_udp[ix*rec.nt], nsam*sizeof(float));
+ }
+
+ /* transform from t-x to kx-w */
+ xorig = ixsrc+ibndx;
+ xt2wkx(rec_vze, crec_vz, ntfft, nkx, ntfft, nkx, xorig);
+ xt2wkx(rec_pe, crec_p, ntfft, nkx, ntfft, nkx, xorig);
+
+ /* apply decomposition operators */
+ kxwdecomp(crec_p, crec_vz, crec_up, crec_dw,
nkx, mod.dx, nsam, dt, fmin, fmax, cp, rho, verbose);
- /* transform back to t-x */
- wkx2xt(crec_up, rec_up, ntfft, nkx, nkx, ntfft, xorig);
- wkx2xt(crec_dw, rec_down, ntfft, nkx, nkx, ntfft, xorig);
-
- /* reduce array to rec.nt samples rec.n traces */
- for (irec=0; irec<rec.n; irec++) {
- ix = rec.x[irec]+ibndx;
- for (it=0; it<rec.nt; it++) {
- rec_up[irec*rec.nt+it] = rec_up[ix*ntfft+it];
- rec_down[irec*rec.nt+it] = rec_down[ix*ntfft+it];
- }
- }
- free(crec_vz);
- free(crec_p);
- free(crec_up);
- free(crec_dw);
- }
- if (rec.type.ud && (mod.ischeme==3 || mod.ischeme==4) ) {
- }
-
- for (irec=0; irec<rec.n; irec++) {
- hdr.tracf = irec+1;
- hdr.tracl = ishot*rec.n+irec+1;
- hdr.gx = 1000*(mod.x0+rec.x[irec]*mod.dx);
- hdr.offset = (rec.x[irec]-ixsrc)*mod.dx;
+ /* transform back to t-x */
+ wkx2xt(crec_up, rec_up, ntfft, nkx, nkx, ntfft, xorig);
+ wkx2xt(crec_dw, rec_down, ntfft, nkx, nkx, ntfft, xorig);
+
+ /* reduce array to rec.nt samples rec.n traces */
+ for (irec=0; irec<rec.n; irec++) {
+ ix = rec.x[irec]+ibndx;
+ for (it=0; it<rec.nt; it++) {
+ rec_up[irec*rec.nt+it] = rec_up[ix*ntfft+it];
+ rec_down[irec*rec.nt+it] = rec_down[ix*ntfft+it];
+ }
+ }
+ free(crec_vz);
+ free(crec_p);
+ free(crec_up);
+ free(crec_dw);
+ }
+ if (rec.type.ud && (mod.ischeme==3 || mod.ischeme==4) ) {
+ }
+
+ for (irec=0; irec<rec.n; irec++) {
+ hdr.tracf = irec+1;
+ hdr.tracl = ishot*rec.n+irec+1;
+ hdr.gx = 1000*(mod.x0+rec.x[irec]*mod.dx);
+ hdr.offset = (rec.x[irec]-ixsrc)*mod.dx;
hdr.gelev = (int)(-1000*(mod.z0+rec.z[irec]*mod.dz));
- if (rec.type.vx) {
- traceWrite( &hdr, &rec_vx[irec*rec.nt], nsam, fpvx) ;
- }
- if (rec.type.vz) {
- traceWrite( &hdr, &rec_vz[irec*rec.nt], nsam, fpvz) ;
- }
- if (rec.type.p) {
- traceWrite( &hdr, &rec_p[irec*rec.nt], nsam, fpp) ;
- }
- if (rec.type.txx) {
- traceWrite( &hdr, &rec_txx[irec*rec.nt], nsam, fptxx) ;
- }
- if (rec.type.tzz) {
- traceWrite( &hdr, &rec_tzz[irec*rec.nt], nsam, fptzz) ;
- }
- if (rec.type.txz) {
- traceWrite( &hdr, &rec_txz[irec*rec.nt], nsam, fptxz) ;
- }
- if (rec.type.pp) {
- traceWrite( &hdr, &rec_pp[irec*rec.nt], nsam, fppp) ;
- }
- if (rec.type.ss) {
- traceWrite( &hdr, &rec_ss[irec*rec.nt], nsam, fpss) ;
- }
- if (rec.type.ud && mod.ischeme==1) {
- traceWrite( &hdr, &rec_up[irec*rec.nt], nsam, fpup) ;
- traceWrite( &hdr, &rec_down[irec*rec.nt], nsam, fpdown) ;
- }
- }
-
- if (rec.type.vx) fclose(fpvx);
- if (rec.type.vz) fclose(fpvz);
- if (rec.type.p) fclose(fpp);
- if (rec.type.txx) fclose(fptxx);
- if (rec.type.tzz) fclose(fptzz);
- if (rec.type.txz) fclose(fptxz);
- if (rec.type.pp) fclose(fppp);
- if (rec.type.ss) fclose(fpss);
- if (rec.type.ud) {
- fclose(fpup);
- fclose(fpdown);
- free(rec_up);
- free(rec_down);
- }
+ if (rec.type.vx) {
+ traceWrite( &hdr, &rec_vx[irec*rec.nt], nsam, fpvx) ;
+ }
+ if (rec.type.vz) {
+ traceWrite( &hdr, &rec_vz[irec*rec.nt], nsam, fpvz) ;
+ }
+ if (rec.type.p) {
+ traceWrite( &hdr, &rec_p[irec*rec.nt], nsam, fpp) ;
+ }
+ if (rec.type.txx) {
+ traceWrite( &hdr, &rec_txx[irec*rec.nt], nsam, fptxx) ;
+ }
+ if (rec.type.tzz) {
+ traceWrite( &hdr, &rec_tzz[irec*rec.nt], nsam, fptzz) ;
+ }
+ if (rec.type.txz) {
+ traceWrite( &hdr, &rec_txz[irec*rec.nt], nsam, fptxz) ;
+ }
+ if (rec.type.pp) {
+ traceWrite( &hdr, &rec_pp[irec*rec.nt], nsam, fppp) ;
+ }
+ if (rec.type.ss) {
+ traceWrite( &hdr, &rec_ss[irec*rec.nt], nsam, fpss) ;
+ }
+ if (rec.type.ud && mod.ischeme==1) {
+ traceWrite( &hdr, &rec_up[irec*rec.nt], nsam, fpup) ;
+ traceWrite( &hdr, &rec_down[irec*rec.nt], nsam, fpdown) ;
+ }
+ }
+
+ if (rec.type.vx) fclose(fpvx);
+ if (rec.type.vz) fclose(fpvz);
+ if (rec.type.p) fclose(fpp);
+ if (rec.type.txx) fclose(fptxx);
+ if (rec.type.tzz) fclose(fptzz);
+ if (rec.type.txz) fclose(fptxz);
+ if (rec.type.pp) fclose(fppp);
+ if (rec.type.ss) fclose(fpss);
+ if (rec.type.ud) {
+ fclose(fpup);
+ fclose(fpdown);
+ free(rec_up);
+ free(rec_down);
+ }
return 0;
}
diff --git a/marchenko/applyMute.c b/marchenko/applyMute.c
index f15ad90..4ba1798 100644
--- a/marchenko/applyMute.c
+++ b/marchenko/applyMute.c
@@ -66,7 +66,7 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
}
}
}
- else if (above==4) { //Psi gate which is the inverse of the Theta gate
+ else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
for (i = 0; i < npossyn; i++) {
imute = xrcvsyn[i];
tmute = mute[isyn*nxs+imute];
diff --git a/marchenko/demo/oneD/README b/marchenko/demo/oneD/README
index 31e3323..22d2b75 100644
--- a/marchenko/demo/oneD/README
+++ b/marchenko/demo/oneD/README
@@ -178,3 +178,10 @@ backpropf2sum_0.15.eps
backpropf2sum_0.30.eps
+The figures in the appendix, to explain the different options in the programs, are reproduced by
+
+==> run figAppendi.scr
+
+-- Figure A-1
+
+-- Figure A-2
diff --git a/marchenko/demo/oneD/epsMarchenkoIter.scr b/marchenko/demo/oneD/epsMarchenkoIter.scr
index 35197df..0c795ed 100755
--- a/marchenko/demo/oneD/epsMarchenkoIter.scr
+++ b/marchenko/demo/oneD/epsMarchenkoIter.scr
@@ -98,3 +98,11 @@ scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
done
+
+#special treatment of f1+ zero-iteration: which is zero, to make a nice gray plot (and not black)
+file=f1plus_001.su
+file_base=${file%.su}
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 bclip=1 wclip=-1> $file_base.eps
+
diff --git a/marchenko/demo/oneD/initialFocus.scr b/marchenko/demo/oneD/initialFocus.scr
index 1bbacbf..4d4fd68 100755
--- a/marchenko/demo/oneD/initialFocus.scr
+++ b/marchenko/demo/oneD/initialFocus.scr
@@ -11,7 +11,7 @@ makemod sizex=10000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=1000 \
intt=def x=-5000,5000 z=400,400 poly=0 cp=2300 ro=3000 \
intt=def x=-5000,5000 z=700,700 poly=0 cp=2000 ro=1100
-#makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
export OMP_NUM_THREADS=1
diff --git a/marchenko/fmute.c b/marchenko/fmute.c
index 9a01568..5abbae6 100644
--- a/marchenko/fmute.c
+++ b/marchenko/fmute.c
@@ -25,7 +25,7 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
-void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int nx, int shift);
+void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift);
double wallclock_time(void);
/*********************** self documentation **********************/
@@ -33,12 +33,12 @@ char *sdoc[] = {
" ",
" fmute - mute in time domain file_shot along curve of maximum amplitude in file_mute ",
" ",
-" fmute file_mute= {file_shot=} [optional parameters]",
+" fmute file_shot= {file_mute=} [optional parameters]",
" ",
" Required parameters: ",
" ",
-" file_mute= ................ input file 1",
-" file_shot= ................ input file 2",
+" file_mute= ................ input file with event that defines the mute line",
+" file_shot= ................ input data thats is muted",
" ",
" Optional parameters: ",
" ",
@@ -47,10 +47,10 @@ char *sdoc[] = {
" shift=0 .................. number of points above(positive) / below(negative) maximum time for mute",
" check=0 .................. plots muting window on top of file_mute: output file check.su",
" scale=0 .................. scale data by dividing through maximum",
-" hw=15 .................... window in time samples to look for maximum in next trace",
+" hw=15 .................... window in time samples to look for maximum in next trace of file_mute",
" smooth=0 ................. number of points to smooth mute with cosine window",
-" nxmax=512 ................ maximum number of traces in input file",
-" ntmax=1024 ............... maximum number of samples/trace in input file",
+//" nxmax=512 ................ maximum number of traces in input file",
+//" ntmax=1024 ............... maximum number of samples/trace in input file",
" verbose=0 ................ silent option; >0 display info",
" ",
" author : Jan Thorbecke : 2012 (janth@xs4all.nl)",
@@ -160,6 +160,7 @@ int main (int argc, char **argv)
f2=f2b;
tmpdata = tmpdata2;
hdrs_in1 = hdrs_in2;
+ sclsxgx = sclshot;
}
if (verbose) vmess("sampling file_mute=%d, file_shot=%d", nt1, nt2);
@@ -227,6 +228,7 @@ int main (int argc, char **argv)
imax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
tmax=0.0;
jmax = 0;
+ xmax=0.0;
for (j = 0; j < nt1; j++) {
lmax = fabs(tmpdata[imax*nt1+j]);
if (lmax > tmax) {
@@ -276,13 +278,14 @@ int main (int argc, char **argv)
if (scale==1) {
for (i = 0; i < nx2; i++) {
lmax = fabs(tmpdata2[i*nt2+maxval[i]]);
- xrcv[i] = i;
for (j = 0; j < nt2; j++) {
tmpdata2[i*nt2+j] = tmpdata2[i*nt2+j]/lmax;
}
}
}
+ for (i = 0; i < nx2; i++) xrcv[i] = i;
+
/*================ apply mute window ================*/
applyMute(tmpdata2, maxval, smooth, above, 1, nx2, nt2, xrcv, nx2, shift);
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 589161f..d0af00d 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -45,14 +45,14 @@ void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn
/*********************** self documentation **********************/
char *sdoc[] = {
" ",
-" MARCHENKO - Iterative Green's functions retrieval in frequency domain",
+" MARCHENKO - Iterative Green's function and focusing functions retrieval",
" ",
-" marchenko file_tinv= file_shot= nshots= [optional parameters]",
+" marchenko file_tinv= file_shot= [optional parameters]",
" ",
" Required parameters: ",
" ",
-" file_tinv= ............... focusing operator(s)",
-" file_shot= ............... shot records with Reflection data",
+" file_tinv= ............... direct arrival from focal point: G_d",
+" file_shot= ............... Reflection response: R",
" ",
" Optional parameters: ",
" ",
@@ -68,17 +68,17 @@ char *sdoc[] = {
" 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",
-" weight=1 ................. weight factor for summation of muted field with Tinv",
+" weight=1 ................. weight factor of R for summation of Ni with G_d",
" 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_pplus= .............. output file with p+ ",
" file_f2= ................. output file with f2 (=p+) ",
" file_pmin= ............... output file with p- ",
-" file_iter= ............... output file with N for each iteration",
+" file_pplus= .............. output file with p+ ",
+" file_iter= ............... output file with -Ni(-t) for each iteration",
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
@@ -135,7 +135,8 @@ int main (int argc, char **argv)
if (!getparfloat("fmax", &fmax)) fmax = 70.0;
if (!getparint("ixa", &ixa)) ixa = 0;
if (!getparint("ixb", &ixb)) ixb = ixa;
- if (!getparint("reci", &reci)) reci = 0;
+// if (!getparint("reci", &reci)) reci = 0;
+ reci=0; // source-receiver reciprocity is not yet fully build into the code
if (!getparfloat("weight", &weight)) weight = 1.0;
if (!getparint("tap", &tap)) tap = 0;
if (!getparint("ntap", &ntap)) ntap = 0;
@@ -382,10 +383,7 @@ int main (int argc, char **argv)
t2 = wallclock_time();
-/*================ construction of Ni(-t) = - \int R(x,t) Fop(t) ================*/
-
- /* set Fop to zero, so new operator can be defined within ixpossyn points */
- //memset(&Fop[0].r, 0, Nsyn*nxs*nw*2*sizeof(float));
+/*================ construction of Ni(-t) = - \int R(x,t) Ni(t) ================*/
synthesis(Refl, Fop, Ni, iRN, nx, nt, nxs, nts, dt, xsyn, Nsyn,
xrcv, xsrc, fxs2, fxs, dxs, dxsrc, dx, ixa, ixb, ntfft, nw, nw_low, nw_high, mode,
@@ -398,34 +396,27 @@ int main (int argc, char **argv)
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nsyn, xsyn, zsyn, iter);
}
/* N_k(x,t) = -N_(k-1)(x,-t) */
+ /* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
+ pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
+ pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+i*nts+j];
}
}
}
+ /* apply mute window based on times of direct arrival (in muteW) */
+ applyMute(Ni, muteW, smooth, above, Nsyn, nxs, nts, ixpossyn, npossyn, shift);
+
/* initialization */
if (iter==0) {
/* N_0(t) = M_0(t) = -p0^-(x,-t) = -(R * T_d^inv)(-t) */
- /* p0^-(x,t) = iRN = (R * T_d^inv)(t) */
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j = 0;
- pmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- pmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+i*nts+j];
- }
- }
- }
-
- applyMute(Ni, muteW, smooth, above, Nsyn, nxs, nts, ixpossyn, npossyn, shift);
-
- /* even iterations: => - f_1^-(-t) = windowed(iRN) */
+ /* zero iteration: => f_1^-(t) = windowed(iRN = -(Ni(-t)) */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -436,6 +427,7 @@ int main (int argc, char **argv)
}
}
+ /* Initialize f2 */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -446,42 +438,11 @@ int main (int argc, char **argv)
}
}
}
- /* Pressure based scheme */
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j=0;
- ix = ixpossyn[i];
- green[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j] + pmin[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- green[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts-j]+ pmin[l*nxs*nts+i*nts+j];
- }
- }
- }
}
else if (iter==1) {
/* Ni(x,t) = -\int R(x,t) M_0(x,-t) dxdt*/
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j = 0;
- pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
- }
- }
- }
- applyMute(Ni, muteW, smooth, above, Nsyn, nxs, nts, ixpossyn, npossyn, shift);
-
- /* Pressure based scheme */
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j=0;
- green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
- }
- }
- }
+ /* Update f2 */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -491,7 +452,8 @@ int main (int argc, char **argv)
}
}
}
- /* odd iterations: M_m^+ */
+
+ /* first iteration: => f_1^+(t) = G_d + windowed(iRN) */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -507,31 +469,7 @@ int main (int argc, char **argv)
/* next iterations */
/* N_k(x,t) = -N_(k-1)(x,-t) */
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j = 0;
- pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- pmin[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
- }
- }
- }
- applyMute(Ni, muteW, smooth, above, Nsyn, nxs, nts, ixpossyn, npossyn, shift);
-
- /* compute full Green's function G = p^+(-t) + p^-(t) */
- if (iter == niter-1) {
- /* Pressure based scheme */
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j=0;
- green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
- }
- }
- }
- } /* end if for last iteration */
-
+ /* update f2 */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -542,8 +480,7 @@ int main (int argc, char **argv)
}
}
-
- if (iter % 2 == 0) { /* even iterations: => - f_1^- (-t) = pmin(t) */
+ if (iter % 2 == 0) { /* even iterations: => f_1^-(t) */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -554,7 +491,7 @@ int main (int argc, char **argv)
}
}
}
- else {/* odd iterations: M_m^+ */
+ else {/* odd iterations: => f_1^+(t) */
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
@@ -568,6 +505,7 @@ int main (int argc, char **argv)
} /* end else (iter!=0) branch */
+
t2 = wallclock_time();
tcopy += t2 - t3;
@@ -577,6 +515,17 @@ int main (int argc, char **argv)
free(Ni);
free(G_d);
+ /* compute full Green's function G = int R * f2(t) + f2(-t) */
+ for (l = 0; l < Nsyn; l++) {
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+
/* compute upgoing Green's function G^+,- */
if (file_gmin != NULL) {
Gmin = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
@@ -623,7 +572,6 @@ int main (int argc, char **argv)
}
} /* end if Gplus */
-
t2 = wallclock_time();
if (verbose) {
vmess("Total CPU-time marchenko = %.3f", t2-t0);
@@ -773,11 +721,11 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
{
int nfreq, size, iox, inx;
float scl;
- int i, j, l, m, iw, ix, ixrcv, k;
- float *rtrace;
- complex *sum, tmp, *ctrace;
+ int i, j, l, m, iw, ix, k;
+ float *rtrace, idxs;
+ complex *sum, *ctrace;
int npe;
- static int first=1;
+ static int first=1, *ixrcv;
static double t0, t1, t;
size = nxs*nts;
@@ -805,6 +753,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
if (!first) {
/* transform muted Ni (Top) to frequency domain, input for next iteration */
for (l = 0; l < Nsyn; l++) {
+ /* set Fop to zero, so new operator can be defined within ixpossyn points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < npossyn; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
@@ -820,6 +769,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
/* transform G_d to frequency domain, over all nxs traces */
first=0;
for (l = 0; l < Nsyn; l++) {
+ /* set Fop to zero, so new operator can be defined within all ix points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < nxs; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
@@ -829,6 +779,13 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
}
}
}
+ idxs = 1.0/dxs;
+ ixrcv = (int *)malloc(nshots*nx*sizeof(int));
+ for (k=0; k<nshots; k++) {
+ for (i = 0; i < nx; i++) {
+ ixrcv[k*nx+i] = NINT((xrcv[k*nx+i]-fxs)*idxs);
+ }
+ }
}
free(ctrace);
t1 = wallclock_time();
@@ -860,8 +817,8 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
shared(iRN, dx, npe, nw, verbose) \
shared(Refl, Nsyn, reci, xrcv, xsrc, xsyn, fxs, nxs, dxs) \
shared(nx, ixa, ixb, dxsrc, iox, inx, k, nfreq, nw_low, nw_high) \
- shared(Fop, size, nts, ntfft, scl, stderr) \
- private(l, ix, j, m, i, ixrcv, sum, rtrace, tmp)
+ shared(Fop, size, nts, ntfft, scl, ixrcv, stderr) \
+ private(l, ix, j, m, i, sum, rtrace)
{ /* start of parallel region */
sum = (complex *)malloc(nfreq*sizeof(complex));
rtrace = (float *)calloc(ntfft,sizeof(float));
@@ -874,14 +831,12 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
/* multiply R with Fop and sum over nx */
memset(&sum[0].r,0,nfreq*2*sizeof(float));
//for (j = 0; j < nfreq; j++) sum[j].r = sum[j].i = 0.0;
- for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- for (i = iox; i < inx; i++) {
- ixrcv = NINT((xrcv[k*nx+i]-fxs)/dxs);
- tmp = Fop[l*nw*nxs+m*nxs+ixrcv];
- sum[j].r += Refl[k*nw*nx+m*nx+i].r*tmp.r -
- Refl[k*nw*nx+m*nx+i].i*tmp.i;
- sum[j].i += Refl[k*nw*nx+m*nx+i].i*tmp.r +
- Refl[k*nw*nx+m*nx+i].r*tmp.i;
+ for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
+ for (i = iox; i < inx; i++) {
+ sum[j].r += Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ixrcv[k*nx+i]].r -
+ Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ixrcv[k*nx+i]].i;
+ sum[j].i += Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ixrcv[k*nx+i]].r +
+ Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ixrcv[k*nx+i]].i;
}
}
diff --git a/utils/ftr1d.c b/utils/ftr1d.c
index 4406c88..05889ca 100644
--- a/utils/ftr1d.c
+++ b/utils/ftr1d.c
@@ -58,7 +58,7 @@ char *sdoc[] = {
NULL};
/**************** end self doc ***********************************/
-void main(int argc, char *argv[])
+int main(int argc, char *argv[])
{
FILE *fp_in, *fp_out;
short trid, trid_out;
diff --git a/utils/syn2d.c b/utils/syn2d.c
index 02ea06b..81aed4a 100644
--- a/utils/syn2d.c
+++ b/utils/syn2d.c
@@ -7,6 +7,10 @@
#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
@@ -543,6 +547,15 @@ void synthesis(float *shotdata, float *syndata, int nx, int nt, int nxs, int nts
complex *sum, *cdata, *shotcdata, tmp, ts, to;
int npe;
+#ifdef _OPENMP
+ npe = omp_get_max_threads();
+ /* parallelisation is over number of virtual source positions (Nsyn) */
+ if (npe > Nsyn) {
+ vmess("Number of OpenMP threads set to %d (was %d)", Nsyn, npe);
+ omp_set_num_threads(Nsyn);
+ }
+#endif
+
/*============= Transform synthesis operator(s), only one time =============*/
--
GitLab