From 70c95cc4659bf9cb44123d44742f79314ee4872f Mon Sep 17 00:00:00 2001
From: JBrackenhoff <J.A.Brackenhoff@tudelft.nl>
Date: Mon, 13 Jan 2020 10:21:35 +0100
Subject: [PATCH] gmin
---
marchenko/applyMute.c | 20 +
marchenko/marchenko.c | 2 +
marchenko3D/Makefile | 2 +
marchenko3D/TWtransform.c | 184 +++++++--
marchenko3D/applyMute3D.c | 23 +-
marchenko3D/getFileInfo3Dzfp.c | 81 ++++
marchenko3D/homogeneousg3D.c | 189 ++++++++-
marchenko3D/homogeneousg3D_backup.c | 577 ++++++++++++++++++++++++++++
marchenko3D/makeWindow3D.c | 6 +-
marchenko3D/marchenko3D.c | 109 +++---
marchenko3D/readShotData3D.c | 2 +
marchenko3D/readShotData3Dzfp.c | 167 ++++++++
utils/Makefile | 23 +-
utils/combine.c | 130 +++++--
14 files changed, 1385 insertions(+), 130 deletions(-)
create mode 100644 marchenko3D/getFileInfo3Dzfp.c
create mode 100644 marchenko3D/homogeneousg3D_backup.c
create mode 100644 marchenko3D/readShotData3Dzfp.c
diff --git a/marchenko/applyMute.c b/marchenko/applyMute.c
index acbbeea..64f8068 100644
--- a/marchenko/applyMute.c
+++ b/marchenko/applyMute.c
@@ -60,6 +60,26 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
}
}
}
+ else if (above==-2){
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ ts = tsynW[isyn*nxs+imute];
+ if (tmute >= nt/2) {
+ memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,-2*ts+tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ }
+ }
else if (above==-1){
for (i = 0; i < npos; i++) {
imute = ixpos[i];
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 637ac8f..99514fd 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -508,6 +508,7 @@ sqrt(energyNi/energyN0));
}
}
}
+ if (above==-2) applyMute(f1min, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
}
else { /* plane wave scheme */
for (l = 0; l < Nfoc; l++) {
@@ -521,6 +522,7 @@ sqrt(energyNi/energyN0));
}
}
}
+ if (above==-2) applyMute(f1min, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
}
}
else {/* odd iterations update: => f_1^+(t) */
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 800660b..3afe49a 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -12,8 +12,10 @@ ALL: marchenko3D fmute3D TWtransform
SRCT = marchenko3D.c \
getFileInfo3D.c \
getFileInfo3DW.c \
+ getFileInfo3Dzfp.c \
readData3D.c \
readShotData3D.c \
+ readShotData3Dzfp.c \
readTinvData3D.c \
synthesis3D.c \
applyMute3D.c \
diff --git a/marchenko3D/TWtransform.c b/marchenko3D/TWtransform.c
index 4c8dba0..8895271 100644
--- a/marchenko3D/TWtransform.c
+++ b/marchenko3D/TWtransform.c
@@ -3,7 +3,9 @@
#include <string.h>
#include <math.h>
#include "segy.h"
+#include "zfpmar.h"
#include <assert.h>
+#include <zfp.h>
typedef struct { /* complex number */
float r,i;
@@ -14,6 +16,7 @@ typedef struct { /* complex number */
#define MIN(x,y) ((x) < (y) ? (x) : (y))
int optncr(int n);
+long zfpcompress(float* data, long nx, long ny, long nz, double tolerance, zfpmar zfpm, FILE *file);
void cc1fft(complex *data, int n, int sign);
void rc1fft(float *rdata, complex *cdata, int n, int sign);
long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
@@ -38,6 +41,9 @@ char *sdoc[] = {
" fmin=0 ................... minimum frequency in the output",
" fmax=70 .................. maximum frequency in the output",
" mode=1 ................... sign of the frequency transform",
+" zfp=0 .................... compress the transformed data using zfp",
+" tolerance=1e-3 ........... accuracy of the zfp compression",
+" weight=2.0 ............... scaling of the reflection data",
" ",
" ",
" author : Joeri Brackenhoff : (j.a.brackenhoff@tudelft.nl)",
@@ -49,14 +55,18 @@ NULL};
int main (int argc, char **argv)
{
FILE *fp, *fp_out;
- segy hdr;
+ segy hdr, *hdr_out;
size_t nread;
long fldr_shot, sx_shot, sy_shot, itrace, one_shot, igath, iw;
long end_of_file, nt, ntfft, nxy, nx, ny, nshots, ret, ntraces;
- long k, l, m, j, nfreq, nw_low, nw_high, nw, mode, verbose;
- float scl, scel, *trace, dt, dx, dy, ft, fx, fy, fmin, fmax, *cdata;
+ long nfreq, nw_low, nw_high, nw, mode, verbose, zfp;
+ long inx, iny, gy;
+ float scl, *trace, dt, dx, dy, ft, fx, fy, fmin, fmax, *cdata, scale;
+ double tolerance;
char *file_T, *file_W;
complex *ctrace;
+ zfptop zfpt;
+ zfpmar zfpm;
initargs(argc, argv);
requestdoc(1);
@@ -67,9 +77,11 @@ int main (int argc, char **argv)
if (file_W==NULL) verr("No file_W is given");
if (!getparfloat("fmin", &fmin)) fmin = 0;
if (!getparfloat("fmax", &fmax)) fmax = 70;
- if (!getparfloat("weight", &scl)) scl = 2.0;
+ if (!getpardouble("tolerance", &tolerance)) tolerance = 1e-3;
+ if (!getparfloat("weight", &scale)) scale = 2.0;
if (!getparlong("mode", &mode)) mode = 1;
if (!getparlong("verbose", &verbose)) verbose = 1;
+ if (!getparlong("zfp", &zfp)) zfp = 0;
nshots = 0;
ret = getFileInfo3D(file_T, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
@@ -112,30 +124,69 @@ int main (int argc, char **argv)
nt = hdr.ns;
dt = hdr.dt/(1E6);
- trace = (float *)calloc(ntfft,sizeof(float));
- ctrace = (complex *)malloc(ntfft*sizeof(complex));
- cdata = (float *)calloc(nw*2,sizeof(float));
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+ cdata = (float *)calloc(nw*2*nxy,sizeof(float));
+ hdr_out = (segy *) calloc(nxy,sizeof(segy));
end_of_file = 0;
one_shot = 1;
igath = 0;
+ if (zfp) {
+ vmess("zfp compression applied");
+ zfpt.dx = dx;
+ zfpt.dy = dy;
+ zfpt.dz = dt;
+ zfpt.ndim = 3;
+ zfpt.ns = nshots;
+ zfpt.scale = hdr.scalco;
+ zfpt.nt = ntfft;
+ zfpt.fmin = fmin;
+ zfpt.fmax = fmax;
+ zfpt.nz = 2*nw;
+ zfpt.tolerance = tolerance;
+ zfpt.fz = ft;
+ zfpt.fx = fx;
+ zfpt.fy = fy;
+
+ nread = fwrite(&zfpt, 1, TOPBYTES, fp_out);
+ assert(nread == TOPBYTES);
+ }
+ else {
+ vmess("no zfp compression applied");
+ }
+
/* Read shots in file */
while (!end_of_file) {
/* start reading data (shot records) */
- itrace = 0;
+ itrace = 0;
+ iny = 1;
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread != TRCBYTES) { /* no more data in file */
break;
}
- sx_shot = hdr.sx;
- sy_shot = hdr.sy;
- fldr_shot = hdr.fldr;
+ sx_shot = hdr.sx;
+ sy_shot = hdr.sy;
+ gy = hdr.gy;
+ fldr_shot = hdr.fldr;
+ if (zfp) {
+ zfpm.gx = hdr.gx;
+ zfpm.gy = hdr.gy;
+ zfpm.sx = hdr.sx;
+ zfpm.sy = hdr.sy;
+ zfpm.sz = hdr.selev;
+ }
while (one_shot) {
+ if (hdr.gy != gy) {
+ iny++;
+ gy = hdr.gy;
+ }
+
nread = fread( trace, sizeof(float), nt, fp );
assert (nread == hdr.ns);
@@ -145,18 +196,16 @@ int main (int argc, char **argv)
rc1fft(trace,ctrace,(int)ntfft,-1);
for (iw=0; iw<nw; iw++) {
- cdata[iw*2] = scl*ctrace[nw_low+iw].r;
- cdata[(iw*2)+1] = scl*mode*ctrace[nw_low+iw].i;
+ cdata[itrace*nw*2+(iw*2)] = scale*ctrace[nw_low+iw].r;
+ cdata[itrace*nw*2+(iw*2)+1] = scale*mode*ctrace[nw_low+iw].i;
}
itrace++;
- hdr.ep = ntfft;
- hdr.ns = 2*nw;
- hdr.unscale = fmin;
- hdr.ungpow = fmax;
-
- ret = writeData3D(fp_out, (float *)&cdata[0], &hdr, 2*nw, 1);
- if (ret < 0 ) verr("error on writing output file.");
+ hdr.ep = ntfft;
+ hdr.ns = 2*nw;
+ hdr.unscale = fmin;
+ hdr.ungpow = fmax;
+ hdr_out[itrace-1] = hdr;
/* read next hdr of next trace */
nread = fread( &hdr, 1, TRCBYTES, fp );
@@ -167,8 +216,17 @@ int main (int argc, char **argv)
}
if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
}
+ inx = itrace/iny;
if (verbose) {
- vmess("finished reading shot x=%li y=%li (%li) with %li traces",sx_shot,sy_shot,igath,itrace);
+ vmess("finished reading shot x=%li y=%li (%li) with %li traces (nx=%li ny=%li) and weight=%.3f",sx_shot,sy_shot,igath,itrace,inx,iny,scale);
+ }
+
+ if (zfp) {
+ zfpcompress(cdata,nx,ny,2*nw,tolerance,zfpm,fp_out);
+ }
+ else {
+ ret = writeData3D(fp_out, (float *)&cdata[0], hdr_out, 2*nw, itrace);
+ if (ret < 0 ) verr("error on writing output file.");
}
if (itrace != 0) { /* end of shot record */
@@ -179,9 +237,91 @@ int main (int argc, char **argv)
end_of_file = 1;
}
}
-
+ fclose(fp_out);
free(ctrace);
free(trace);
exit(0);
+}
+
+long zfpcompress(float* data, long nx, long ny, long nz, double tolerance, zfpmar zfpm, FILE *file)
+{
+
+ zfp_field* field = NULL;
+ zfp_stream* zfp = NULL;
+ bitstream* stream = NULL;
+ void* fi = NULL;
+ void* fo = NULL;
+ void* buffer = NULL;
+ size_t rawsize = 0;
+ size_t zfpsize = 0;
+ size_t bufsize = 0;
+ size_t nwrite;
+ zfp_exec_policy exec = zfp_exec_serial;
+
+ zfp = zfp_stream_open(NULL);
+ field = zfp_field_alloc();
+
+ zfp_field_set_pointer(field, (void *)data);
+
+ zfp_field_set_type(field, zfp_type_float);
+ zfp_field_set_size_3d(field, (uint)nz, (uint)nx, (uint)ny);
+
+ zfp_stream_set_accuracy(zfp, tolerance);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ bufsize = zfp_stream_maximum_size(zfp, field);
+ if (!bufsize) {
+ fprintf(stderr, "invalid compression parameters\n");
+ return EXIT_FAILURE;
+ }
+
+ buffer = malloc(bufsize);
+ if (!buffer) {
+ fprintf(stderr, "cannot allocate memory\n");
+ return EXIT_FAILURE;
+ }
+
+ stream = stream_open(buffer, bufsize);
+ if (!stream) {
+ fprintf(stderr, "cannot open compressed stream\n");
+ return EXIT_FAILURE;
+ }
+ zfp_stream_set_bit_stream(zfp, stream);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ zfpsize = zfp_compress(zfp, field);
+ if (zfpsize == 0) {
+ fprintf(stderr, "compression failed\n");
+ return EXIT_FAILURE;
+ }
+
+ zfpm.nx = nx;
+ zfpm.ny = ny;
+ zfpm.compsize = zfpsize;
+
+ // file = fopen(zfppath, "wb");
+ // if (file==NULL) {
+ // fprintf(stderr,"input file %s has an error\n", zfppath);
+ // perror("error in opening file: ");
+ // fflush(stderr);
+ // return -1;
+ // }
+ nwrite = fwrite(&zfpm, 1, MARBYTES, file);
+ assert(nwrite == MARBYTES);
+ if (fwrite(buffer, 1, zfpsize, file) != zfpsize) {
+ fprintf(stderr, "cannot write compressed file\n");
+ return EXIT_FAILURE;
+ }
+ // fclose(file);
+
+ return 1;
}
\ No newline at end of file
diff --git a/marchenko3D/applyMute3D.c b/marchenko3D/applyMute3D.c
index b50fd44..03097a2 100644
--- a/marchenko3D/applyMute3D.c
+++ b/marchenko3D/applyMute3D.c
@@ -41,10 +41,10 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
}
}
}
- else if (above==0){
+ else if (above==0){ //Classic Theta window removes Gd
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
- tmute = mute[isyn*nxs+imute];
+ tmute = mute[isyn*nxys+imute];
if (tmute >= nt/2) {
memset(&data[isyn*nxys*nt+i*nt],0, sizeof(float)*nt);
continue;
@@ -60,6 +60,25 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
}
}
}
+ else if (above==-2){ //New Theta window keeps Gd
+ for (i = 0; i < npos; i++) {
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ if (tmute >= nt/2) {
+ memset(&data[isyn*nxys*nt+i*nt],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,tmute+shift),l=0; j < MAX(0,tmute+shift+smooth); j++,l++) {
+ data[isyn*nxys*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
+ data[isyn*nxys*nt+i*nt+j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ data[isyn*nxys*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ }
+ }
else if (above==-1){
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
diff --git a/marchenko3D/getFileInfo3Dzfp.c b/marchenko3D/getFileInfo3Dzfp.c
new file mode 100644
index 0000000..36854f5
--- /dev/null
+++ b/marchenko3D/getFileInfo3Dzfp.c
@@ -0,0 +1,81 @@
+#define _FILE_OFFSET_BITS 64
+#define _LARGEFILE_SOURCE
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <math.h>
+#include "segy.h"
+#include "zfpmar.h"
+#include <zfp.h>
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* gets sizes, sampling and min/max values of a SU file
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+void vmess(char *fmt, ...);
+void verr(char *fmt, ...);
+int optncr(int n);
+
+long getFileInfo3Dzfp(char *filename, long *n1, long *n2, long *n3, long *ngath,
+ float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
+ float *fmin, float *fmax, float *scl, long *nxm)
+{
+ FILE *fp_in;
+ size_t nread;
+ zfpmar zfpm;
+ zfptop zfpt;
+ long ishot, compsize;
+
+ fp_in = fopen(filename, "r");
+ if (fp_in==NULL) {
+ fprintf(stderr,"input file %s has an error\n", fp_in);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ nread = fread(&zfpt, 1, TOPBYTES, fp_in);
+ assert(nread == TOPBYTES);
+ *ngath = zfpt.ns;
+ *n1 = zfpt.nt;
+ *n2 = 1;
+ *n3 = 1;
+ *fmin = zfpt.fmin;
+ *fmax = zfpt.fmax;
+ *f1 = zfpt.fz;
+ *f2 = zfpt.fx;
+ *f3 = zfpt.fy;
+ *d1 = zfpt.dz;
+ *d2 = zfpt.dx;
+ *d3 = zfpt.dy;
+ *nxm = 1;
+
+ if (zfpt.scale < 0.0) *scl = 1.0/fabs((float)zfpt.scale);
+ else if (zfpt.scale == 0.0) *scl = 1.0;
+ else *scl = zfpt.scale;
+
+ compsize = 0;
+
+ for (ishot=0; ishot<zfpt.ns; ishot++) {
+ fseek(fp_in, compsize, SEEK_CUR);
+ nread = fread(&zfpm, 1, MARBYTES, fp_in);
+ assert(nread == MARBYTES);
+ *n2 = MAX(*n2,zfpm.nx);
+ *n3 = MAX(*n3,zfpm.ny);
+ compsize = zfpm.compsize;
+ }
+
+ *nxm = (*n2)*(*n3);
+
+ return 0;
+}
\ No newline at end of file
diff --git a/marchenko3D/homogeneousg3D.c b/marchenko3D/homogeneousg3D.c
index e5b2cfd..5da29e0 100644
--- a/marchenko3D/homogeneousg3D.c
+++ b/marchenko3D/homogeneousg3D.c
@@ -24,6 +24,7 @@ void conjugate(float *data, long nsam, long nrec, float dt);
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 convol2(float *data1, float *data2, float *con, long nrec, long nsam, float dt, float fmin, float fmax, long opt);
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 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);
@@ -32,7 +33,7 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
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 homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx, long ny, long nt, float dx, float dy,
+void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1m, float *zsyn, long nx, long ny, long nt, float dx, float dy,
float dt, long Nfoc, long *sx, long *sy, long *sz, long verbose)
{
char *file_inp;
@@ -127,6 +128,36 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx,
else if (scheme==5) {
if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with f2 source");
}
+ else if (scheme==8) { // 0=f1p 1=f1m
+ if (verbose) vmess("f1+ redatuming");
+ if (n_source<2) verr("Not enough input for the homogeneous Green's function");
+ for (k = 0; k < ny; k++) {
+ depthDiff(&input[0*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ conjugate(&input[0*ny*nx*nt+k*nx*nt], nt, nx, dt);
+ depthDiff(&input[1*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ conjugate(&input[1*ny*nx*nt+k*nx*nt], nt, nx, dt);
+ }
+ }
+ else if (scheme==9) { // 0=f1p 1=f1m
+ if (verbose) vmess("f1- redatuming");
+ if (n_source<2) verr("Not enough input for the homogeneous Green's function");
+ for (k = 0; k < ny; k++) {
+ depthDiff(&input[0*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ depthDiff(&input[1*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ }
+ }
+ else if (scheme==10) {
+ if (verbose) vmess("2i IM(f1) redatuming");
+ inputjkz = (float *)calloc(n_source*nx*ny*nt,sizeof(float));
+ for (k = 0; k < ny; k++) {
+ for (l = 0; l < nx*nt; l++) {
+ inputjkz[k*nx*nt+l] = input[k*nx*nt+l];
+ }
+ conjugate(&inputjkz[k*nx*nt], nt, nx, dt);
+ depthDiff(&inputjkz[k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ depthDiff(&input[k*nx*nt] , nt, nx, dt, dx, fmin, fmax, cp, 1);
+ }
+ }
else {
if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with G source");
}
@@ -139,7 +170,7 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx,
tmp1 = (float *)calloc(nx*nt,sizeof(float));
tmp2 = (float *)calloc(nx*nt,sizeof(float));
}
- if (scheme==3) tmp1 = (float *)calloc(nx*nt,sizeof(float));
+ if (scheme==3 || scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nx*nt,sizeof(float));
#pragma omp for schedule(guided,1)
for (l = 0; l < Nfoc; l++) {
@@ -181,6 +212,42 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx,
}
}
}
+ else if (scheme==8) { // f1+ redatuming 0=f1p 1=f1m
+ for (k = 0; k < ny; k++) {
+ convol2(&input[0*ny*nx*nt+k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
+ convol2(&input[1*ny*nx*nt+k*nx*nt], &f1m[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] -= 2.0*(conv[i*nt+j] + tmp1[i*nt+j])/rho;
+ HomG[j*Nfoc+l] -= 2.0*(conv[i*nt+(j+nt/2)] + tmp1[i*nt+(j+nt/2)])/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==9) { // f1- redatuming 0=f1p 1=f1m
+ for (k = 0; k < ny; k++) {
+ convol2(&input[0*ny*nx*nt+k*nx*nt], &f1m[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
+ convol2(&input[1*ny*nx*nt+k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*(conv[i*nt+j] + tmp1[i*nt+j])/rho;
+ HomG[j*Nfoc+l] += 2.0*(conv[i*nt+(j+nt/2)] + tmp1[i*nt+(j+nt/2)])/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==10) { // 2i IM(f1) redatuming
+ for (k = 0; k < ny; k++) {
+ convol2(&input[k*nx*nt] , &f1m[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 2);
+ convol2(&inputjkz[k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 2);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 4.0*(conv[i*nt+j] - tmp1[i*nt+j])/rho;
+ HomG[j*Nfoc+l] += 4.0*(conv[i*nt+(j+nt/2)] - tmp1[i*nt+(j+nt/2)])/rho;
+ }
+ }
+ }
+ }
else if (scheme==1) { //classical representation
for (k = 0; k < ny; k++) {
convol(&greenjkz[l*ny*nx*nt+k*nx*nt], &input[k*nx*nt], tmp1, nx, nt, dt, 0);
@@ -202,12 +269,11 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx,
else if (scheme==2) { //Marchenko representation without time-reversal G source
for (k = 0; k < ny; k++) {
depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- convol(&input[k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
- timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ convol2(&input[k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 2*conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += 2*conv[i*nt+(j+nt/2)]/rho;
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2.0*conv[i*nt+(j+nt/2)]/rho;
}
}
}
@@ -219,8 +285,8 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx,
timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 2*conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += 2*conv[i*nt+(j+nt/2)]/rho;
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2.0*conv[i*nt+(j+nt/2)]/rho;
}
}
}
@@ -242,8 +308,8 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx,
timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 2*conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += 2*conv[i*nt+(j+nt/2)]/rho;
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2.0*conv[i*nt+(j+nt/2)]/rho;
}
}
}
@@ -575,3 +641,106 @@ void conjugate(float *data, long nsam, long nrec, float dt)
return;
}
+
+void convol2(float *data1, float *data2, float *con, long nrec, long nsam, float dt, float fmin, float fmax, long opt)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, sign, i, j, n;
+ float omin, omax, deltom, om, df, dw, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccon, tmp, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ 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);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (long)MIN((omin/deltom), (nfreq));
+ iomin = MAX(iomin, 1);
+ iomax = MIN((long)(omax/deltom), (nfreq));
+
+ cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ if (opt==1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*ccon[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = -om*ccon[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt==2) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = -om*ccon[ix*nfreq+iom].r;
+ }
+ }
+ }
+ free(ccon);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,con,nsam);
+
+ free(cdatascl);
+ free(rdata1);
+ free(rdata2);
+
+ return;
+ return;
+}
\ No newline at end of file
diff --git a/marchenko3D/homogeneousg3D_backup.c b/marchenko3D/homogeneousg3D_backup.c
new file mode 100644
index 0000000..e5b2cfd
--- /dev/null
+++ b/marchenko3D/homogeneousg3D_backup.c
@@ -0,0 +1,577 @@
+#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))
+
+double wallclock_time(void);
+
+void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
+void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, float *datout);
+void conjugate(float *data, long nsam, long nrec, float dt);
+
+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);
+
+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 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);
+
+//void kxwfilter(float *data, long nt, long nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc);
+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 homogeneousg3D(float *HomG, float *green, float *f2p, float *zsyn, long nx, long ny, long nt, float dx, float dy,
+ float dt, long Nfoc, long *sx, long *sy, long *sz, long verbose)
+{
+ char *file_inp;
+ long i, j, k, l, ret, scheme, count=0, n_source;
+ long is, kxwfilt, n1, n2, n3, ngath, ntraces, zsrc;
+ float d1, d2, d3, f1, f2, f3, scl;
+ float *conv, fmin, fmax, alpha, cp, rho, perc;
+ float *greenjkz, *input, *inputjkz, *tmp1, *tmp2;
+ double t0, t2, tfft;
+ segy *hdr_inp;
+
+ if (!getparstring("file_inp", &file_inp)) file_inp = NULL;
+ if (!getparlong("scheme", &scheme)) scheme = 0;
+ if (!getparlong("kxwfilt", &kxwfilt)) kxwfilt = 0;
+ if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (!getparfloat("fmax", &fmax)) fmax = 100.0;
+ if (!getparfloat("alpha", &alpha)) alpha = 65.0;
+ if (!getparfloat("cp", &cp)) cp = 1000.0;
+ if (!getparfloat("rho", &rho)) rho = 1000.0;
+ if (!getparfloat("perc", &perc)) perc = 0.15;
+
+ tfft = 0.0;
+ ret = 0;
+ t0 = wallclock_time();
+
+ /* Read in the source function input and check the size */
+ n_source=0;
+ ret = getFileInfo3D(file_inp, &n1, &n2, &n3, &n_source, &d1, &d2, &d3, &f1, &f2, &f3, &scl, &ntraces);
+
+ if (n1!=nt) verr("number of t-samples between input (%li) and retrieved (%li) is not equal",n1,nt);
+ if (n2!=nx) verr("number of x-samples between input (%li) and retrieved (%li) is not equal",n2,nx);
+ if (n3!=ny) verr("number of y-samples between input (%li) and retrieved (%li) is not equal",n3,ny);
+ if (NINT(d1*1000.0)!=NINT(dt*1000.0)) verr("dt sampling between input (%.3e) and retrieved (%.3e) is not equal",d1,dt);
+ if (NINT(d2*1000.0)!=NINT(dx*1000.0)) verr("dx sampling between input (%.3e) and retrieved (%.3e) is not equal",d2,dx);
+ if (NINT(d3*1000.0)!=NINT(dy*1000.0)) verr("dy sampling between input (%.3e) and retrieved (%.3e) is not equal",d3,dy);
+
+ input = (float *)calloc(n_source*nx*ny*nt,sizeof(float));
+ hdr_inp = (segy *)calloc(n_source*nx*ny,sizeof(segy));
+ readSnapData3D(file_inp, input, hdr_inp, n_source, nx, ny, nt, 0, nx, 0, ny, 0, nt);
+
+ zsrc = labs(hdr_inp[0].selev);
+ for (l = 0; l < nx*ny; l++) {
+ sx[l] = hdr_inp[0].sx;
+ sy[l] = hdr_inp[0].sy;
+ sz[l] = hdr_inp[0].sdepth;
+ }
+ for (l = 0; l < n_source; l++) {
+ sx[l] = hdr_inp[l*nx*ny].sx;
+ sy[l] = hdr_inp[l*nx*ny].sy;
+ sz[l] = hdr_inp[l*nx*ny].sdepth;
+ }
+
+ if (scheme==1) { // Scale the Green's functions if the classical scheme is used
+ if (verbose) vmess("Classical Homogeneous Green's function retrieval");
+ greenjkz = (float *)calloc(Nfoc*nx*ny*nt,sizeof(float));
+ inputjkz = (float *)calloc(n_source*nx*ny*nt,sizeof(float));
+
+ for (l = 0; l < Nfoc; l++) {
+ for (k = 0; k < ny; k++) {
+ depthDiff(&greenjkz[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ }
+ }
+
+ for (l = 0; l < n_source; l++) {
+ for (i = 0; i < ny*nx*nt; i++) {
+ inputjkz[l*ny*nx*nt+i] = input[l*ny*nx*nt+i];
+ }
+ conjugate(&inputjkz[l*ny*nx*nt], nt, nx*ny, dt);
+ conjugate(&input[l*ny*nx*nt], nt, nx*ny, dt);
+ for (k = 0; k < ny; k++) {
+ depthDiff(&inputjkz[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ }
+ }
+ }
+ else if (scheme==2) {
+ if (verbose) vmess("Marchenko Green's function retrieval with G source");
+ }
+ else if (scheme==6) {
+ if (verbose) vmess("Marchenko Green's function retrieval with f2 source");
+ }
+ else if (scheme==7) {
+ if (verbose) vmess("Marchenko Green's function retrieval with source depending on position");
+ }
+ else if (scheme==3) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with multiple sources");
+ if (verbose) vmess("Looping over %li source positions",n_source);
+ }
+ else if (scheme==4) {
+ if (verbose) vmess("Back propagation with multiple sources");
+ if (verbose) vmess("Looping over %li source positions",n_source);
+ }
+ else if (scheme==5) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with f2 source");
+ }
+ else {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with G source");
+ }
+
+#pragma omp parallel default(shared) \
+ private(i,j,k,is,conv,tmp1,tmp2,zsrc)
+{
+ conv = (float *)calloc(nx*nt,sizeof(float));
+ if (scheme==1) {
+ tmp1 = (float *)calloc(nx*nt,sizeof(float));
+ tmp2 = (float *)calloc(nx*nt,sizeof(float));
+ }
+ if (scheme==3) tmp1 = (float *)calloc(nx*nt,sizeof(float));
+
+#pragma omp for schedule(guided,1)
+ for (l = 0; l < Nfoc; l++) {
+
+ count+=1;
+ zsrc = hdr_inp[0].selev;
+
+ if (verbose>2) vmess("Creating Homogeneous G at location %li out of %li",l,Nfoc);
+
+ if (scheme==0) { //Marchenko representation with G source
+ for (k = 0; k < ny; k++) {
+ depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(&input[k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -3);
+ if (kxwfilt) {
+ //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
+ }
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==5) { //Marchenko representation with f2 source
+ for (k = 0; k < ny; k++) {
+ depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(&input[k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -3);
+ if (kxwfilt) {
+ //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
+ }
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==1) { //classical representation
+ for (k = 0; k < ny; k++) {
+ convol(&greenjkz[l*ny*nx*nt+k*nx*nt], &input[k*nx*nt], tmp1, nx, nt, dt, 0);
+ convol(&green[l*ny*nx*nt+k*nx*nt], &inputjkz[k*nx*nt], tmp2, nx, nt, dt, 0);
+ for (i = 0; i < nx; i++) {
+ for (j = 0; j < nt; j++) {
+ conv[i*nt+j] = tmp1[i*nt+j]+tmp2[i*nt+j];
+ }
+ }
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==2) { //Marchenko representation without time-reversal G source
+ for (k = 0; k < ny; k++) {
+ depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(&input[k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 2*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2*conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==6) { //Marchenko representation without time-reversal f2 source
+ for (k = 0; k < ny; k++) {
+ depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(&input[k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 2*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2*conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==7) { //Marchenko representation without time-reversal using varying sources
+ for (k = 0; k < ny; k++) {
+ depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<(n_source/2); is+=2) {
+ zsrc = labs(hdr_inp[is].selev);
+ if (zsrc > NINT(1000.0*zsyn[l])) {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses G source (zsrc=%li)",NINT(1000.0*zsyn[l]));
+ convol(&input[is*ny*nx*nt+k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ }
+ else {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses f2 source (zsrc=%li)",NINT(1000.0*zsyn[l]));
+ convol(&input[(is+1)*ny*nx*nt+k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ }
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 2*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2*conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ }
+ else if (scheme==3) { //Marchenko representation with multiple shot gathers
+ for (k = 0; k < ny; k++) {
+ depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<n_source; is++) {
+ convol(&input[is*ny*nx*nt+k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -3);
+ if (kxwfilt) {
+ //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
+ }
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[is*nt*Nfoc+(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
+ HomG[is*nt*Nfoc+j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ }
+ else if (scheme==4) { //Marchenko representation with multiple shot gathers without time-reversal
+ for (k = 0; k < ny; k++) {
+ depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<n_source; is++) {
+ convol(&input[is*ny*nx*nt+k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ if (kxwfilt) {
+ //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
+ }
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[is*nt*Nfoc+(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
+ HomG[is*nt*Nfoc+j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ }
+
+ }
+ free(conv);
+ if (scheme==1) {
+ free(tmp1);
+ free(tmp2);
+ }
+ if (scheme==3) free(tmp1);
+ if (scheme==4) free(tmp1);
+}
+ if (scheme==1) {
+ free(input);
+ free(inputjkz);
+ free(greenjkz);
+ }
+ free(hdr_inp);
+
+ t2 = wallclock_time();
+ if (verbose) {
+ vmess("Total Homogeneous G time = %.3f", t2-t0);
+ }
+
+ return;
+}
+
+void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (long)MIN((omin/deltom), (nfreq));
+ iomin = MAX(iomin, 1);
+ iomax = MIN((long)(omax/deltom), (nfreq));
+
+ cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ if (opt == 1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ cdatascl[ix*nfreq+iom].r = -om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt == -1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = -om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt == -2) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 4.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].r;
+ cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].i;
+ }
+ }
+ else if (opt == -3) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = 2*om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ free(cdata);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt)
+{
+ 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;
+
+ optn = optncr(nsam);
+ nfreq = optncr(nsam)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nrec);
+ dkx = 2.0*PI/(nkx*dx);
+ cdata = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad2d_data(data,nsam,nrec,optn,nkx,rdata);
+
+ /* double forward FFT */
+ xt2wkx(&rdata[0], &cdata[0], optn, nkx, optn, nkx, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (long)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((long)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (iom = 0; iom < iomin; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
+
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+ }
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkx2xt(&cdatascl[0], &rdata[0], optn, nkx, nkx, optn, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, float *datout)
+{
+ long it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsam+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+ for (ix=nrec;ix<nrecout;ix++) {
+ for (it=0;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+}
+
+void conjugate(float *data, long nsam, long nrec, float dt)
+{
+ long optn, nfreq, j, ix, it, sign, ntdiff;
+ float *rdata, scl;
+ complex *cdata;
+
+ optn = optncr(nsam);
+ ntdiff = optn-nsam;
+ nfreq = optn/2+1;
+
+ 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);
+
+ /* take complex conjugate */
+ for(ix = 0; ix < nrec; ix++) {
+ for(j = 0; j < nfreq; j++) cdata[ix*nfreq+j].i = -cdata[ix*nfreq+j].i;
+ }
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdata[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsam ; it++)
+ data[ix*nsam+it] = scl*rdata[ix*optn+it+ntdiff];
+ }
+ //scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdata);
+ free(rdata);
+
+ return;
+}
diff --git a/marchenko3D/makeWindow3D.c b/marchenko3D/makeWindow3D.c
index 2809e63..c2f6f64 100644
--- a/marchenko3D/makeWindow3D.c
+++ b/marchenko3D/makeWindow3D.c
@@ -64,6 +64,7 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
rc1fft(wavelet,cwav,ntfft,-1);
free(wavtmp);
free(wavelet);
+ vmess("Read the wavelet");
}
@@ -150,7 +151,10 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
}
}
free(timeval); free(hdrs_mute);
- if (file_amp!=NULL) free(amp); free(hdrs_amp);
+ if (file_amp!=NULL) {
+ free(amp);
+ free(hdrs_amp);
+ }
return;
}
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 02047cb..fdfb590 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -36,6 +36,8 @@ typedef struct _complexStruct { /* complex number */
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 readShotData3Dzfp(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
+ long *xnx, complex *cdata, long nw, long nshots, long nx, long ny, 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);
@@ -52,6 +54,9 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath,
long getFileInfo3DW(char *filename, long *n1, long *n2, long *n3, long *ngath,
float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, float *fmin, float *fmax,
float *sclsxgxsygy, long *nxm);
+long getFileInfo3Dzfp(char *filename, long *n1, long *n2, long *n3, long *ngath,
+ float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
+ float *fmin, float *fmax, float *scl, 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);
@@ -76,7 +81,7 @@ long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2,
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 *f2p, float *zsyn, long nx, long ny, long nt, float dx, float dy,
+void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1m, float *zsyn, long nx, long ny, long nt, float dx, float dy,
float dt, long Nfoc, long *sx, long *sy, long *sz, long verbose);
long linearsearch(long *array, size_t N, long value);
@@ -169,17 +174,17 @@ int main (int argc, char **argv)
long hw, smooth, above, shift, *ixpos, *iypos, npos, ix, iy, nzim, nxim, nyim, compact;
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;
+ 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, *Nig, *trace, *Gmin, *Gplus, *HomG;
float xmin, xmax, ymin, ymax, scale, tsq, Q, f0, *tmpdata;
float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim;
- float grad2rad, p, src_angle, src_velo;
+ float grad2rad, p, src_angle, src_velo, *mutetest;
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter, *file_imag, *file_homg, *file_ampscl;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin, *file_inp;
- char *file_ray, *file_amp, *file_wav, *file_shotw;
+ char *file_ray, *file_amp, *file_wav, *file_shotw, *file_shotzfp;
segy *hdrs_out, *hdrs_Nfoc, *hdrs_iter;
initargs(argc, argv);
@@ -190,7 +195,8 @@ int main (int argc, char **argv)
if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
if (!getparstring("file_shotw", &file_shotw)) file_shotw = NULL;
- if (file_shot==NULL && file_shotw==NULL) verr("No input for the shot data given");
+ if (!getparstring("file_shotzfp", &file_shotzfp)) file_shotzfp = NULL;
+ if (file_shot==NULL && file_shotw==NULL && file_shotzfp==NULL) verr("No input for the shot data given");
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;
@@ -209,8 +215,8 @@ int main (int argc, char **argv)
if (file_homg!=NULL && file_inp==NULL) verr("Cannot create HomG if no file_inp is given");
if (!getparstring("file_ampscl", &file_ampscl)) file_ampscl = NULL;
if (!getparlong("verbose", &verbose)) verbose = 0;
- if (file_tinv == NULL && file_shot == NULL && file_shotw == NULL)
- verr("file_tinv, file_shotw and file_shot cannot be both input pipe");
+ if (file_tinv == NULL && file_shot == NULL && file_shotw == NULL && file_shotzfp==NULL)
+ verr("file_tinv, file_shotw, file_shotzfp and file_shot cannot all be input pipe");
if (!getparstring("file_green", &file_green)) {
if (verbose) vwarn("parameter file_green not found, assume pipe");
file_green = NULL;
@@ -274,12 +280,18 @@ int main (int argc, char **argv)
if (verbose) vmess("Retrieved file info of the first arrivals");
ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
- if (file_shot!=NULL) {
- ret = getFileInfo3D(file_shot, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+ if (file_shotzfp!=NULL) {
+ vmess("zfp shot data");
+ ret = getFileInfo3Dzfp(file_shotzfp, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &fmin, &fmax, &scl, &ntraces);
}
else if (file_shotw!=NULL) {
+ vmess("Frequency shot data");
ret = getFileInfo3DW(file_shotw, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &fmin, &fmax, &scl, &ntraces);
}
+ else if (file_shot!=NULL) {
+ vmess("Time shot data");
+ ret = getFileInfo3D(file_shot, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+ }
if (verbose) vmess("Retrieved file info of the shot data");
nshots = ngath;
assert (nxs*nys >= nshots);
@@ -319,6 +331,7 @@ int main (int argc, char **argv)
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)
iypos = (long *)calloc(nys*nxs,sizeof(long)); // y-position of source of shot in G_d domain (nxs*nys with dxs, dys)
+ energyN0= (double *)calloc(Nfoc,sizeof(double)); // minimum energy for each focal position
Refl = (complex *)malloc(nw*ny*nx*nshots*sizeof(complex));
tapersh = (float *)malloc(nx*sizeof(float));
@@ -380,6 +393,18 @@ int main (int argc, char **argv)
}
}
+ // mutetest = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
+ // for (i=0; i<Nfoc*nys*nxs*ntfft; i++) {
+ // mutetest[i] = 1.0;
+ // }
+ // applyMute3D(mutetest, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
+
+ // fp_out = fopen( "mute.bin", "w+" );
+ // for (i=0; i<nx*ny; i++) {
+ // fprintf(fp_out,"%.7f\n",mutetest[i]);
+ // }
+ // fclose(fp_out);
+
/* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
@@ -444,16 +469,20 @@ int main (int argc, char **argv)
/*================ Reading shot records ================*/
- if (file_shot!=NULL) {
- mode=1;
- readShotData3D(file_shot, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw,
- nw_low, nshots, nx, ny, ntfft, mode, scale, verbose);
+ if (file_shotzfp!=NULL){
+ readShotData3Dzfp(file_shotzfp, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl,
+ nw, nshots, nx, ny, scl, verbose);
}
- else {
+ else if (file_shotw!=NULL) {
mode=0;
readShotData3D(file_shotw, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw,
nw_low, nshots, nx, ny, ntfft, mode, scale, verbose);
}
+ else if (file_shot!=NULL) {
+ mode=1;
+ readShotData3D(file_shot, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw,
+ nw_low, nshots, nx, ny, ntfft, mode, scale, verbose);
+ }
mode=1;
tapersh = (float *)malloc(nx*sizeof(float));
@@ -488,7 +517,7 @@ int main (int argc, char **argv)
nyshot = nshots/nxshot;
fxf = xsrc[0];
- if (nx > 1) dxf = (xrcv[ny*nx-1] - xrcv[0])/(float)(nx-1);
+ if (nx > 1) dxf = xrcv[1] - xrcv[0];
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);
@@ -497,7 +526,7 @@ int main (int argc, char **argv)
vmess("dx used => %f", dx);
}
fyf = ysrc[0];
- if (ny > 1) dyf = (yrcv[ny*nx-1] - yrcv[0])/(float)(ny-1);
+ if (ny > 1) dyf = yrcv[nx] - yrcv[0];
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);
@@ -536,7 +565,7 @@ int main (int argc, char **argv)
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, x:%li, y%li, z:%li", Nfoc, nxim, nyim, nzim);
+ vmess("Number of focusing operators = %li, x:%li, y:%li, z:%li", Nfoc, nxim, nyim, nzim);
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);
@@ -548,6 +577,7 @@ int main (int argc, char **argv)
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("frequency cutoffs = min:%.3f max:%.3f",fmin,fmax);
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);
@@ -717,9 +747,9 @@ int main (int argc, char **argv)
}
}
}
- if (iter==0) energyN0 = energyNi;
+ if (iter==0) energyN0[l] = 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));
+ l, iter, sqrt(energyNi), sqrt(energyNi/energyN0[l]));
}
/* apply mute window based on times of direct arrival (in muteW) */
@@ -1013,7 +1043,7 @@ int main (int argc, char **argv)
// Determine Homogeneous Green's function
HomG = (float *)calloc(Nfoc*ntfft,sizeof(float));
- homogeneousg3D(HomG, green, f2p, zsyn, nxs, nys, ntfft, dxs, dys, dt, Nfoc, sx, sy, sz, verbose);
+ homogeneousg3D(HomG, green, f2p, f1plus, f1min, zsyn, nxs, nys, ntfft, dxs, dys, dt, Nfoc, sx, sy, sz, verbose);
// Set headers and write out the data
fp_homg = fopen(file_homg, "w+");
@@ -1095,9 +1125,10 @@ int main (int argc, char **argv)
/*================ write output files ================*/
-
- fp_out = fopen(file_green, "w+");
- if (fp_out==NULL) verr("error on creating output file %s", file_green);
+ if (file_green != NULL) {
+ 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);
@@ -1148,8 +1179,10 @@ int main (int argc, char **argv)
}
}
- ret = writeData3D(fp_out, (float *)&green[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
+ if (file_green != NULL) {
+ 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);
@@ -1168,37 +1201,15 @@ int main (int argc, char **argv)
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_green != NULL) {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);}
@@ -1231,7 +1242,7 @@ long unique_elements(float *arr, long len)
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 ((arr[inner] >= arr[outer]-0.1) && (arr[inner] <= arr[outer]+0.1)) is_unique = 0;
}
if (is_unique) ++unique;
}
diff --git a/marchenko3D/readShotData3D.c b/marchenko3D/readShotData3D.c
index 3c477da..6bedbe8 100644
--- a/marchenko3D/readShotData3D.c
+++ b/marchenko3D/readShotData3D.c
@@ -59,10 +59,12 @@ long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
dt = hdr.dt/(1E6);
if (mode==0){
+ if (verbose) vmess("Reading in frequency traces");
trace = (float *)calloc(2*nw,sizeof(float));
ctrace = (complex *)malloc(nw*sizeof(complex));
}
else {
+ if (verbose) vmess("Reading in time traces");
trace = (float *)calloc(ntfft,sizeof(float));
ctrace = (complex *)malloc(ntfft*sizeof(complex));
}
diff --git a/marchenko3D/readShotData3Dzfp.c b/marchenko3D/readShotData3Dzfp.c
new file mode 100644
index 0000000..8de346f
--- /dev/null
+++ b/marchenko3D/readShotData3Dzfp.c
@@ -0,0 +1,167 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+#include "zfpmar.h"
+#include <zfp.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+int optncr(int n);
+void cc1fft(complex *data, int n, int sign);
+void rc1fft(float *rdata, complex *cdata, int n, int sign);
+long zfpdecompress(float* data, long nx, long ny, long nz, long comp, double tolerance, FILE *file);
+
+long readShotData3Dzfp(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
+ long *xnx, complex *cdata, long nw, long nshots, long nx, long ny, float scale, long verbose)
+{
+ FILE *fp;
+ size_t nread;
+ long sx_shot, sy_shot, iw, nxy, ishot, inx, iny, ix, iy;
+ float *trace, dx, dy, scl;
+ double tolerance;
+ zfpmar zfpm;
+ zfptop zfpt;
+
+ nxy = nx*ny;
+
+ /* Reading first header */
+
+ if (filename == NULL) fp = stdin;
+ else fp = fopen( filename, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", filename);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &zfpt, 1, TOPBYTES, fp );
+ assert(nread == TOPBYTES);
+
+ if (verbose) vmess("Reading in zfp frequency traces");
+
+ tolerance = zfpt.tolerance;
+ dx = zfpt.dx;
+ dy = zfpt.dy;
+
+ if (zfpt.scale < 0.0) scl = 1.0/fabs((float)zfpt.scale);
+ else if (zfpt.scale == 0.0) scl = 1.0;
+ else scl = zfpt.scale;
+
+ trace = (float *)calloc(2*nw,sizeof(float));
+
+ /* Read shots in file */
+
+ for (ishot=0; ishot<nshots; ishot++) {
+
+ /* start reading data (shot records) */
+ nread = fread( &zfpm, 1, MARBYTES, fp );
+ if (nread != MARBYTES) { /* no more data in file */
+ break;
+ }
+
+ inx = zfpm.nx;
+ iny = zfpm.ny;
+ sx_shot = zfpm.sx;
+ sy_shot = zfpm.sy;
+ xsrc[ishot] = sx_shot*scl;
+ ysrc[ishot] = sy_shot*scl;
+ zsrc[ishot] = zfpm.sz*scl;
+ xnx[ishot] = inx*iny;
+
+ trace = (float *)realloc(trace, inx*iny*2*nw*sizeof(float));
+ zfpdecompress(trace, inx, iny, 2*nw, zfpm.compsize , tolerance, fp);
+
+ for (iy = 0; iy < iny; iy++) {
+ for (ix = 0; ix < inx; ix++) {
+ for (iw = 0; iw < nw; iw++) {
+ cdata[ishot*nxy*nw+iw*nxy+iy*inx+ix].r = trace[iy*inx*nw*2+ix*nw*2+(iw*2)];
+ cdata[ishot*nxy*nw+iw*nxy+iy*inx+ix].i = trace[iy*inx*nw*2+ix*nw*2+(iw*2)+1];
+ }
+
+ xrcv[ishot*nxy+iy*inx+ix] = (zfpm.gx*scl)+(ix*dx);
+ yrcv[ishot*nxy+iy*inx+ix] = (zfpm.gy*scl)+(iy*dy);
+ }
+ }
+ if (verbose>2) {
+ vmess("finished reading shot x=%li y=%li (%li) with %li traces (nx=%li ny=%li)",sx_shot,sy_shot,ishot,inx*iny,inx,iny);
+ }
+ }
+
+ free(trace);
+
+ return 0;
+}
+
+long zfpdecompress(float* data, long nx, long ny, long nz, long comp, double tolerance, FILE *file)
+{
+ zfp_field* field = NULL;
+ zfp_stream* zfp = NULL;
+ bitstream* stream = NULL;
+ zfp_exec_policy exec = zfp_exec_serial;
+ size_t nread, compsize;
+ void *buffer;
+
+ zfp = zfp_stream_open(NULL);
+ field = zfp_field_alloc();
+ compsize = comp;
+
+ buffer = malloc(compsize);
+ if (!buffer) {
+ fprintf(stderr, "cannot allocate memory\n");
+ return EXIT_FAILURE;
+ }
+ nread = fread((uchar*)buffer, 1, compsize, file);
+ assert(nread==compsize);
+
+ stream = stream_open(buffer, compsize);
+ if (!stream) {
+ fprintf(stderr, "cannot open compressed stream\n");
+ return EXIT_FAILURE;
+ }
+ zfp_stream_set_bit_stream(zfp, stream);
+
+ zfp_field_set_type(field, zfp_type_float);
+ zfp_field_set_size_3d(field, (uint)nz, (uint)nx, (uint)ny);
+
+ zfp_stream_set_accuracy(zfp, tolerance);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ zfp_stream_rewind(zfp);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ zfp_field_set_pointer(field, (void *)data);
+
+ while (!zfp_decompress(zfp, field)) {
+ /* fall back on serial decompression if execution policy not supported */
+ if (zfp_stream_execution(zfp) != zfp_exec_serial) {
+ if (!zfp_stream_set_execution(zfp, zfp_exec_serial)) {
+ fprintf(stderr, "cannot change execution policy\n");
+ return EXIT_FAILURE;
+ }
+ }
+ else {
+ fprintf(stderr, "decompression failed\n");
+ return EXIT_FAILURE;
+ }
+ }
+
+ return 1;
+}
\ No newline at end of file
diff --git a/utils/Makefile b/utils/Makefile
index 1efdfb8..4caf5ee 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 makeR1D
SRCM = \
makemod.c \
@@ -153,6 +153,15 @@ SRCCO = combine.c \
docpkge.c \
readSnapData3D.c
+SRCMR = makeR1D.c \
+ getFileInfo3D.c \
+ writeData3D.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData3D.c
SRCCI = combine_induced.c \
getFileInfo3D.c \
@@ -260,6 +269,11 @@ OBJCO = $(SRCCO:%.c=%.o)
combine: $(OBJCO)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine $(OBJCO) $(LIBS)
+OBJMR = $(SRCMR:%.c=%.o)
+
+makeR1D: $(OBJMR)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o makeR1D $(OBJMR) $(LIBS)
+
OBJCI = $(SRCCI:%.c=%.o)
combine_induced: $(OBJCI)
@@ -280,7 +294,7 @@ OBJSS = $(SRCSS:%.c=%.o)
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 makeR1D
cp makemod $B
cp makewave $B
cp extendModel $B
@@ -298,9 +312,10 @@ install: makemod makewave extendModel fconv correigen green green3D basop syn2d
cp reshape_su $B
cp HomG $B
cp snap2shot $B
+ cp makeR1D $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) makeR1D $(OBJMR) 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 $B/makeR1D
diff --git a/utils/combine.c b/utils/combine.c
index 4c9da58..b2bf1e4 100755
--- a/utils/combine.c
+++ b/utils/combine.c
@@ -48,17 +48,19 @@ char *sdoc[] = {
" numb=0 ................... integer number of first file",
" dnumb=1 .................. integer number of increment in files",
" nzmax=0 .................. Maximum number of files read",
+" direction=z .............. The direction over which the data is stacked",
" verbose=1 ................ Give detailed information of process",
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];
+ char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100], *direction;
float *indata, *outdata;
float dz, dy, dx, z0, y0, x0, scl;
long nt, nz, ny, nx, ntr, ix, iy, it, is, iz, pos, file_det, nzs, dt;
long numb, dnumb, ret, nzmax, compact, verbose, nxyz, sx, sy, sz;
+ long nxout, nyout, nzout, ixout, iyout, izout;
segy *hdr_in, *hdr_out;
initargs(argc, argv);
@@ -69,6 +71,7 @@ int main (int argc, char **argv)
*----------------------------------------------------------------------------*/
if (!getparstring("file_in", &fin)) fin = NULL;
if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparstring("direction", &direction)) direction = "z";
if (!getparlong("numb", &numb)) numb=0;
if (!getparlong("dnumb", &dnumb)) dnumb=1;
if (!getparlong("nzmax", &nzmax)) nzmax=0;
@@ -76,16 +79,22 @@ int main (int argc, char **argv)
if (!getparlong("compact", &compact)) compact=0;
if (fin == NULL) verr("Incorrect downgoing input");
+
+ vmess("Direction given is %s",direction);
+ if (strcmp(direction,"x") != 0 && strcmp(direction,"y") != 0 && strcmp(direction,"z") != 0) {
+ verr("Direction needs to be either x, y or z");
+ }
+
/*----------------------------------------------------------------------------*
* 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,"%li",numb);
+ sprintf(numb1,"%s%li",direction,numb);
ptr = strstr(fin,numb1);
pos = ptr - fin + 1;
sprintf(fbegin,"%*.*s", pos-1, pos-1, fin);
- sprintf(fend,"%s", fin+pos);
+ sprintf(fend,"%s", fin+pos+1);
/*----------------------------------------------------------------------------*
* Determine the amount of files that are present
@@ -93,7 +102,7 @@ int main (int argc, char **argv)
file_det = 1;
nzs=0;
while (file_det) { // Check for a file with the filename
- sprintf(fins,"%li",nzs*dnumb+numb);
+ sprintf(fins,"%s%li",direction,nzs*dnumb+numb);
sprintf(fin,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin, "r");
if (fp_in == NULL) { // If the filename does not exist
@@ -124,7 +133,7 @@ int main (int argc, char **argv)
* Read in the first two files and determine the header values
* of the output
*----------------------------------------------------------------------------*/
- sprintf(fins,"%li",numb);
+ sprintf(fins,"%s%li",direction,numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
nt = 1;
if (compact==0) {
@@ -140,11 +149,27 @@ int main (int argc, char **argv)
readCompact(fin2, &nx, &ny, &nz, &nt, &dx, &dy, &dz, &dt, &x0, &y0, &z0, &sx, &sy, &sz, &scl);
}
- nxyz = nx*ny*nzs*nz;
+ if (strcmp(direction,"z") == 0) {
+ nxout = nx;
+ nyout = ny;
+ nzout = nz*nzs;
+ }
+ if (strcmp(direction,"y") == 0) {
+ nxout = nx;
+ nyout = ny*nzs;
+ nzout = nz;
+ }
+ if (strcmp(direction,"x") == 0) {
+ nxout = nx*nzs;
+ nyout = ny;
+ nzout = nz;
+ }
+
+ nxyz = nxout*nyout*nzout;
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("Number of virtual receivers: %li, x: %li, y: %li, z: %li",nxyz,nxout,nyout,nzout);
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);
}
@@ -153,7 +178,7 @@ int main (int argc, char **argv)
* Read in a single file to determine if the header values match
* and allocate the data
*----------------------------------------------------------------------------*/
- hdr_out = (segy *)calloc(nx*ny,sizeof(segy));
+ hdr_out = (segy *)calloc(nxout*nyout,sizeof(segy));
outdata = (float *)calloc(nxyz*nt,sizeof(float));
indata = (float *)calloc(nx*ny*nz*nt,sizeof(float));
if (compact != 2) {
@@ -168,7 +193,7 @@ int main (int argc, char **argv)
*----------------------------------------------------------------------------*/
for (is = 0; is < nzs; is++) {
if (verbose) vmess("Combining file %li out of %li",is+1,nzs);
- sprintf(fins,"%li",is*dnumb+numb);
+ sprintf(fins,"%s%li",direction,is*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin2, "r");
if (fp_in == NULL) {
@@ -177,14 +202,21 @@ int main (int argc, char **argv)
fclose(fp_in);
if (compact==0) {
readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
- if (dz == 1.0) {
- if (is==1) dz = hdr_in[0].f1-z0;
+ if (is==1) {
+ if (dz == 1.0) dz = hdr_in[0].f1-z0;
+ if (dx == 1.0) dx = hdr_in[0].f2-x0;
}
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
+ if (strcmp(direction,"y") == 0) iyout = iy+is;
+ else iyout = iy;
for (ix = 0; ix < nx; ix++) {
+ if (strcmp(direction,"x") == 0) ixout = ix+is;
+ else ixout = 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];
+ if (strcmp(direction,"z") == 0) izout = iz+is;
+ else izout = iz;
+ outdata[it*nyout*nxout*nzout+iyout*nxout*nzout+ixout*nzout+izout] = indata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
}
}
}
@@ -192,14 +224,21 @@ int main (int argc, char **argv)
}
else if (compact==1) {
readSnapData3D(fin2, indata, hdr_in, nt, nz, ny, nx, 0, nz, 0, ny, 0, nx);
- if (dz == 1.0) {
- if (is==1) dz = hdr_in[0].f1-z0;
+ if (is==1) {
+ if (dz == 1.0) dz = hdr_in[0].f1-z0;
+ if (dx == 1.0) dx = hdr_in[0].f2-x0;
}
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
+ if (strcmp(direction,"y") == 0) iyout = iy+is;
+ else iyout = iy;
for (ix = 0; ix < nx; ix++) {
+ if (strcmp(direction,"x") == 0) ixout = ix+is;
+ else ixout = 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*nz*nx+iy*nz*nx+iz*nx+ix];
+ if (strcmp(direction,"z") == 0) izout = iz+is;
+ else izout = iz;
+ outdata[it*nyout*nxout*nzout+iyout*nxout*nzout+ixout*nzout+izout] = indata[it*ny*nz*nx+iy*nz*nx+iz*nx+ix];
}
}
}
@@ -207,14 +246,21 @@ int main (int argc, char **argv)
}
else if (compact==2) {
readSnapData3D(fin2, indata, hdr_in, 1, 1, 1, nx*ny*nz*nt, 0, 1, 0, 1, 0, nx*ny*nz*nt);
- if (dz == 1.0) {
- if (is==1) dz = hdr_in[0].f1-z0;
+ if (is==1) {
+ if (dz == 1.0) dz = hdr_in[0].f1-z0;
+ if (dx == 1.0) dx = hdr_in[0].f2-x0;
}
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
+ if (strcmp(direction,"y") == 0) iyout = iy+is;
+ else iyout = iy;
for (ix = 0; ix < nx; ix++) {
+ if (strcmp(direction,"x") == 0) ixout = ix+is;
+ else ixout = 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];
+ if (strcmp(direction,"z") == 0) izout = iz+is;
+ else izout = iz;
+ outdata[it*nyout*nxout*nzout+iyout*nxout*nzout+ixout*nzout+izout] = indata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
}
}
}
@@ -236,32 +282,32 @@ int main (int argc, char **argv)
fp_out = fopen(fout, "w+");
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 = roundf(z0*1000.0)/1000.0;
- hdr_out[iy*nx+ix].f2 = roundf(x0*1000.0)/1000.0;
- hdr_out[iy*nx+ix].dt = dt;
- hdr_out[iy*nx+ix].d1 = roundf(dz*1000.0)/1000.0;
- hdr_out[iy*nx+ix].d2 = roundf(dx*1000.0)/1000.0;
- 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;
+ for (iy = 0; iy < nyout; iy++) {
+ for (ix = 0; ix < nxout; ix++) {
+ hdr_out[iy*nxout+ix].fldr = it+1;
+ hdr_out[iy*nxout+ix].tracl = it*nyout*nxout+iy*nxout+ix+1;
+ hdr_out[iy*nxout+ix].tracf = iy*nxout+ix+1;
+ hdr_out[iy*nxout+ix].scalco = -1000;
+ hdr_out[iy*nxout+ix].scalel = -1000;
+ hdr_out[iy*nxout+ix].sdepth = sz;
+ hdr_out[iy*nxout+ix].selev = -sz;
+ hdr_out[iy*nxout+ix].trid = 2;
+ hdr_out[iy*nxout+ix].ns = nzout;
+ hdr_out[iy*nxout+ix].trwf = nxout*nyout;
+ hdr_out[iy*nxout+ix].ntr = hdr_out[iy*nxout+ix].fldr*hdr_out[iy*nxout+ix].trwf;
+ hdr_out[iy*nxout+ix].f1 = roundf(z0*1000.0)/1000.0;
+ hdr_out[iy*nxout+ix].f2 = roundf(x0*1000.0)/1000.0;
+ hdr_out[iy*nxout+ix].dt = dt;
+ hdr_out[iy*nxout+ix].d1 = roundf(dz*1000.0)/1000.0;
+ hdr_out[iy*nxout+ix].d2 = roundf(dx*1000.0)/1000.0;
+ hdr_out[iy*nxout+ix].sx = sx;
+ hdr_out[iy*nxout+ix].gx = (int)roundf(x0 + (ix*dx))*1000;
+ hdr_out[iy*nxout+ix].sy = sy;
+ hdr_out[iy*nxout+ix].gy = (int)roundf(y0 + (iy*dy))*1000;
+ hdr_out[iy*nxout+ix].offset = (hdr_out[iy*nxout+ix].gx - hdr_out[iy*nxout+ix].sx)/1000.0;
}
}
- ret = writeData3D(fp_out, &outdata[it*nxyz], hdr_out, nz*nzs, nx*ny);
+ ret = writeData3D(fp_out, &outdata[it*nxyz], hdr_out, nzout, nxout*nyout);
if (ret < 0 ) verr("error on writing output file.");
}
fclose(fp_out);
--
GitLab