From e624e8d2abf3984a2ebdca06d731b1faff911e55 Mon Sep 17 00:00:00 2001
From: JBrackenhoff <J.A.Brackenhoff@tudelft.nl>
Date: Mon, 27 May 2019 16:51:48 +0200
Subject: [PATCH] Joeri update
---
marchenko3D/Makefile | 1 +
marchenko3D/fmute3D.c | 4 +-
marchenko3D/homogeneousg3D.c | 565 +++++++++++++++++++++++
marchenko3D/marchenko3D.c | 124 +++--
raytime/Makefile | 2 +-
utils/HomG.c | 858 +++++++++++++++++++++++++++++++++++
utils/Makefile | 106 ++++-
utils/MuteSnap.c | 227 +++++++++
utils/combine.c | 198 ++++++++
utils/combine_induced.c | 208 +++++++++
utils/getFileInfo3D.c | 244 ++++++++++
utils/readSnapData.c | 58 +++
utils/readSnapData3D.c | 56 +++
utils/reshape_su.c | 125 +++++
utils/snap2shot.c | 205 +++++++++
utils/writeData3D.c | 28 ++
16 files changed, 2974 insertions(+), 35 deletions(-)
create mode 100644 marchenko3D/homogeneousg3D.c
create mode 100755 utils/HomG.c
create mode 100644 utils/MuteSnap.c
create mode 100755 utils/combine.c
create mode 100755 utils/combine_induced.c
create mode 100644 utils/getFileInfo3D.c
create mode 100755 utils/readSnapData.c
create mode 100644 utils/readSnapData3D.c
create mode 100755 utils/reshape_su.c
create mode 100644 utils/snap2shot.c
create mode 100644 utils/writeData3D.c
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 691bfc1..a2aad66 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -20,6 +20,7 @@ SRCT = marchenko3D.c \
makeWindow3D.c \
ampest3D.c \
imaging3D.c \
+ homogeneousg3D.c \
readSnapData3D.c \
writeDataIter.c \
wallclock_time.c \
diff --git a/marchenko3D/fmute3D.c b/marchenko3D/fmute3D.c
index 9c52073..627c745 100644
--- a/marchenko3D/fmute3D.c
+++ b/marchenko3D/fmute3D.c
@@ -442,9 +442,9 @@ long main (int argc, char **argv)
for (l=0; l<ny1; l++) {
for (i=0; i<nx1; i++) {
jmax = maxval[l*nx1+i]-shift;
- ret = fprintf(fp_psline1, "%.5f %.5f \n",jmax*dt,hdrs_in1[l*nx1+i].gx*sclshot,hdrs_in1[l*nx1+i].gy*sclshot);
+ ret = fprintf(fp_psline1, "%.5f %.5f %.5f \n",jmax*dt,hdrs_in1[l*nx1+i].gx*sclshot,hdrs_in1[l*nx1+i].gy*sclshot);
jmax =-maxval[l*nx1+i]+shift;
- ret = fprintf(fp_psline2, "%.5f %.5f \n",jmax*dt,hdrs_in1[l*nx1+i].gx*sclshot,hdrs_in1[l*nx1+i].gy*sclshot);
+ ret = fprintf(fp_psline2, "%.5f %.5f %.5f \n",jmax*dt,hdrs_in1[l*nx1+i].gx*sclshot,hdrs_in1[l*nx1+i].gy*sclshot);
}
}
}
diff --git a/marchenko3D/homogeneousg3D.c b/marchenko3D/homogeneousg3D.c
new file mode 100644
index 0000000..22c887e
--- /dev/null
+++ b/marchenko3D/homogeneousg3D.c
@@ -0,0 +1,565 @@
+#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);
+
+//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 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);
+
+ 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;
+}
\ No newline at end of file
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 21e71c1..17bd022 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -70,7 +70,7 @@ void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, long nx, l
void imaging3D(float *Image, float *Gmin, float *f1plus, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
-void homogeneousg3D(float *HomG, float *green, float *f2, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
+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 verbose);
long linearsearch(long *array, size_t N, long value);
@@ -90,6 +90,7 @@ char *sdoc[] = {
" Optional parameters: ",
" ",
" INTEGRATION ",
+" ampest=0 ................. Estimate a scalar amplitude correction with depth (=1)",
" tap=0 .................... lateral taper focusing(1), shot(2) or both(3)",
" ntap=0 ................... number of taper points at boundaries",
" fmin=0 ................... minimum frequency in the Fourier transform",
@@ -98,6 +99,7 @@ char *sdoc[] = {
" niter=10 ................. number of iterations",
" MUTE-WINDOW ",
" file_amp= ................ amplitudes for the raytime estimation",
+" file_wav= ................ Wavelet applied to the raytime data",
" above=0 .................. mute above(1), around(0) or below(-1) the first travel times of file_tinv",
" shift=12 ................. number of points above(positive) / below(negative) travel time for mute",
" hw=8 ..................... window in time samples to look for maximum in next trace",
@@ -110,8 +112,16 @@ char *sdoc[] = {
" HOMOGENEOUS GREEN'S FUNCTION RETRIEVAL OPTIONS ",
" file_inp= ................ Input source function for the retrieval",
" scheme=0 ................. Scheme for homogeneous Green's function retrieval",
-" .......................... Scheme for homogeneous Green's function retrieval",
+" .......................... scheme=0 Marchenko HomG retrieval with time-reversal",
+" .......................... scheme=1 Classical retrieval scheme",
+" .......................... scheme=2 Marchenko HomG retrieval without time-reversal",
+" .......................... scheme=3 Back propagation with multiple sources",
+" .......................... scheme=4 Marchenko HomG retrieval with multiple sources",
" kxwfilt=0 ................ Apply a dip filter before integration",
+" alpha=65.0 ............... Alpha filter for the kxwfilter",
+" perc=0.15 ................ Percentage for the kxwfilter",
+" cp=1000.0 ................ Velocity of upper layer for certain operations",
+" rho=1000.0 ............... Density of upper layer for certain operations",
" OUTPUT DEFINITION ",
" file_green= .............. output file with full Green function(s)",
" file_gplus= .............. output file with G+ ",
@@ -136,7 +146,7 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_out, *fp_f1plus, *fp_f1min, *fp_imag;
+ FILE *fp_out, *fp_f1plus, *fp_f1min, *fp_imag, *fp_homg;
FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin, *fp_amp;
long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
@@ -149,12 +159,12 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
float d1, d2, d3, f1, f2, f3, fxsb, fxse, fysb, fyse, ft, fx, fy, *xsyn, *ysyn, dxsrc, dysrc;
float *green, *f2p, *pmin, *G_d, dt, dx, dy, dxs, dys, scl, mem;
- float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
+ float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus, *HomG;
float xmin, xmax, ymin, ymax, scale, tsq, Q, f0, *tmpdata;
float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim, dyim, dxim;
complex *Refl, *Fop;
- char *file_tinv, *file_shot, *file_green, *file_iter, *file_imag, *file_ampscl;
- char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
+ char *file_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;
segy *hdrs_out, *hdrs_Nfoc;
@@ -178,6 +188,9 @@ int main (int argc, char **argv)
if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
if (!getparstring("file_wav", &file_wav)) file_wav = NULL;
if (!getparstring("file_imag", &file_imag)) file_imag = NULL;
+ if (!getparstring("file_homg", &file_homg)) file_homg = NULL;
+ if (!getparstring("file_inp", &file_inp)) file_inp = NULL;
+ 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)
@@ -206,6 +219,12 @@ int main (int argc, char **argv)
if (reci && ntap) vwarn("tapering influences the reciprocal result");
+ if (file_inp!=NULL) {
+ fp_out = fopen( file_inp, "r" );
+ if (fp_out == NULL) verr("File %s does not exist or cannot be opened", file_inp);
+ fclose(fp_out);
+ }
+
/*================ Reading info about shot and initial operator sizes ================*/
ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
@@ -303,11 +322,10 @@ int main (int argc, char **argv)
/*Determine the shape of the focal positions*/
nzim = unique_elements(zsyn,Nfoc);
- if (nzim>1) dzim = zsyn[1]-zsyn[0];
- else dzim = 1.0;
nyim = unique_elements(ysyn,Nfoc);
nxim = unique_elements(xsyn,Nfoc);
-
+ if (nzim>1) dzim = zsyn[nyim*nxim]-zsyn[0];
+ else dzim = 1.0;
/* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
@@ -456,7 +474,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", Nfoc);
+ 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);
@@ -749,12 +767,12 @@ int main (int argc, char **argv)
for (l=0; l<Nfoc; l++) {
for (j=0; j<nxs*nys*nts; j++) {
green[l*nxs*nts+j] *= ampscl[l];
+ f2p[l*nxs*nys*nts+j] *= ampscl[l];
+ pmin[l*nxs*nys*nts+j] *= ampscl[l];
+ f1plus[l*nxs*nys*nts+j] *= ampscl[l];
+ f1min[l*nxs*nys*nts+j] *= ampscl[l];
if (file_gplus != NULL) Gplus[l*nxs*nys*nts+j] *= ampscl[l];
if (file_gmin != NULL || file_imag != NULL) Gmin[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f2 != NULL) f2p[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_pmin != NULL) pmin[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f1plus != NULL || file_imag != NULL) f1plus[l*nxs*nys*nts+j] *= ampscl[l];
- if (file_f1min != NULL) f1min[l*nxs*nys*nts+j] *= ampscl[l];
}
if (verbose>1) vmess("Amplitude of focal position %li is equal to %.3e",l,ampscl[l]);
}
@@ -770,15 +788,15 @@ int main (int argc, char **argv)
hdrs_Nfoc[l*nxim+j].trid = 1;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
- hdrs_Nfoc[l*nxim+j].sx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].gx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].gy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l*nxim+j]*(1e3);
hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
- hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
- hdrs_Nfoc[l*nxim+j].d2 = xsyn[1]-xsyn[0];
+ hdrs_Nfoc[l*nxim+j].d1 = dzim;
+ hdrs_Nfoc[l*nxim+j].d2 = dxs;
hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
@@ -799,10 +817,12 @@ int main (int argc, char **argv)
/* Apply imaging*/
if (file_imag!=NULL) {
+
// Determine Image
Image = (float *)calloc(Nfoc,sizeof(float));
imaging3D(Image, Gmin, f1plus, nxs, nys, ntfft, dxs, dys, dt, Nfoc, verbose);
if (file_gmin==NULL) free(Gmin);
+
// Set headers
hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
for (l=0; l<nyim; l++){
@@ -814,20 +834,21 @@ int main (int argc, char **argv)
hdrs_Nfoc[l*nxim+j].trid = 1;
hdrs_Nfoc[l*nxim+j].scalco = -1000;
hdrs_Nfoc[l*nxim+j].scalel = -1000;
- hdrs_Nfoc[l*nxim+j].sx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].gx = xsyn[j]*(1e3);
- hdrs_Nfoc[l*nxim+j].gy = ysyn[l]*(1e3);
- hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l*nxim+j]*(1e3);
hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
- hdrs_Nfoc[l*nxim+j].d1 = zsyn[1]-zsyn[0];
- hdrs_Nfoc[l*nxim+j].d2 = xsyn[1]-xsyn[0];
+ hdrs_Nfoc[l*nxim+j].d1 = dzim;
+ hdrs_Nfoc[l*nxim+j].d2 = dxs;
hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
}
}
+
// Write out image
fp_imag = fopen(file_imag, "w+");
if (fp_imag==NULL) verr("error on creating output file %s", file_imag);
@@ -838,6 +859,53 @@ int main (int argc, char **argv)
free(Image);
}
+ /* Determine homogeneous Green's function*/
+ if (file_homg!=NULL) {
+
+ // Determine Image
+ HomG = (float *)calloc(Nfoc*ntfft,sizeof(float));
+ homogeneousg3D(HomG, green, f2p, zsyn, nxs, nys, ntfft, dxs, dys, dt, Nfoc, verbose);
+
+ fp_homg = fopen(file_homg, "w+");
+ if (fp_homg==NULL) verr("error on creating output file %s", file_homg);
+ hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
+
+ for (i=0; i<ntfft; i++) {
+ // Set headers
+ for (l=0; l<nyim; l++){
+ for (j=0; j<nxim; j++){
+ hdrs_Nfoc[l*nxim+j].ns = nzim;
+ hdrs_Nfoc[l*nxim+j].fldr = i+1;
+ hdrs_Nfoc[l*nxim+j].tracl = 1;
+ hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
+ hdrs_Nfoc[l*nxim+j].trid = 1;
+ hdrs_Nfoc[l*nxim+j].scalco = -1000;
+ hdrs_Nfoc[l*nxim+j].scalel = -1000;
+ hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
+ hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
+ hdrs_Nfoc[l*nxim+j].d1 = dzim;
+ hdrs_Nfoc[l*nxim+j].d2 = dxs;
+ hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
+ hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
+ hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
+ }
+ }
+
+ // Write out homogeneous Green's function
+ ret = writeData3D(fp_homg, (float *)&HomG[i*Nfoc], hdrs_Nfoc, nzim, nxim*nyim);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_homg);
+ free(hdrs_Nfoc);
+ free(HomG);
+ }
+
t2 = wallclock_time();
if (verbose) {
vmess("Total CPU-time marchenko = %.3f", t2-t0);
diff --git a/raytime/Makefile b/raytime/Makefile
index 398bd1f..bde8b02 100644
--- a/raytime/Makefile
+++ b/raytime/Makefile
@@ -49,7 +49,7 @@ OBJC = $(SRCC:%.c=%.o)
$(PRG): $(OBJC) raytime.h
$(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o raytime $(OBJC) $(LIBS)
-install: raytime
+install: Raytime
cp raytime $B
clean:
diff --git a/utils/HomG.c b/utils/HomG.c
new file mode 100755
index 0000000..570bf7f
--- /dev/null
+++ b/utils/HomG.c
@@ -0,0 +1,858 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+double wallclock_time(void);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+int topdet(float *data, int nt);
+void conjugate(float *data, int nsam, int nrec, float dt);
+
+void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout);
+void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout);
+void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift);
+void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt, int shift);
+void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax);
+void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt);
+void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt);
+void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout);
+
+char *sdoc[] = {
+" ",
+" HomG - Calculate a Homogeneous Green's function ",
+" ",
+" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. First file of the array of receivers",
+" file_shot= ............... File containing the shot location",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first snapshot file",
+" dnumb= ................... integer number of increment in snapshot files",
+" zmax= .................... Integer number of maximum depth level",
+" inx= ..................... Number of sources per depth level",
+" zrcv= .................... z-coordinate of first receiver location",
+" xrcv= .................... x-coordinate of first receiver location",
+" dzrcv= ................... z-spacing of receivers",
+" dxrcv= ................... x-spacing of receivers",
+" shift=0.0 ................ shift per shot",
+" scheme=0 ................. Scheme used for retrieval. 0=Marchenko,",
+" 1=Marchenko with multiple sources, 2=classical",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_shot, *fp_out;
+ char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
+ float *indata, *Ghom, *shotdata, *shotdata_jkz, rho, fmin, fmax;
+ float dt, dx, t0, x0, xmin, xmax1, sclsxgx, f1, f2, dxrcv, dzrcv;
+ float *conv, *conv2, *tmp1, *tmp2, cp, shift;
+ int nshots, nt, nx, ntraces, ret, ix, it, is, ir, ig, file_det, nxs, nzs, verbose;
+ int pos1, npos, zmax, inx, numb, dnumb, count, scheme, ntmax, ntshift, shift_num;
+ segy *hdr_in, *hdr_out, *hdr_shot;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("fin", &fin)) fin = NULL;
+ if (!getparstring("fshot", &fshot)) fshot = NULL;
+ if (!getparstring("fout", &fout)) fout = "out.su";
+ if (!getparint("zmax", &zmax)) zmax = 0;
+ if (!getparint("inx", &inx)) inx = 0;
+ if (!getparfloat("zrcv", &f1)) f1 = 0;
+ if (!getparfloat("xrcv", &f2)) f2 = 0;
+ if (!getparfloat("dzrcv", &dzrcv)) dzrcv = -1;
+ if (!getparfloat("dxrcv", &dxrcv)) dxrcv = -1;
+ if (!getparfloat("rho", &rho)) rho=1000.0;
+ if (!getparfloat("cp", &cp)) cp = 1500.0;
+ if (!getparfloat("fmin", &fmin)) fmin=0.0;
+ if (!getparfloat("fmax", &fmax)) fmax=100.0;
+ if (!getparfloat("shift", &shift)) shift=0.0;
+ if (!getparint("numb", &numb)) numb=0;
+ if (!getparint("dnumb", &dnumb)) dnumb=1;
+ if (!getparint("scheme", &scheme)) scheme = 0;
+ if (!getparint("ntmax", &ntmax)) ntmax = 0;
+ if (!getparint("verbose", &verbose)) verbose = 0;
+ if (fin == NULL) verr("Incorrect f2 input");
+ if (fshot == NULL) verr("Incorrect Green input");
+
+ if (dnumb == 0) dnumb = 1;
+
+ sprintf(fins,"z%d",numb);
+
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ verr("error on opening basefile=%s", fin);
+ }
+ fclose(fp_in);
+
+ ptr = strstr(fin,fins);
+ pos1 = ptr - fin + 1;
+
+ sprintf(fbegin,"%*.*s", pos1-1, pos1-1, fin);
+ sprintf(fend,"%s", fin+pos1+1);
+
+ file_det = 1;
+ nzs=0;
+
+ while (file_det) {
+ sprintf(fins,"z%d",nzs*dnumb+numb);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ if (nzs == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzs == 1) {
+ vmess("1 file detected");
+ file_det = 0;
+ break;
+ }
+ else {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzs++;
+ }
+
+ if (inx < 1) {
+ inx = 1;
+ }
+
+ if (zmax < 1) zmax=1;
+ if (zmax < nzs) nzs=zmax;
+
+ nxs = inx;
+ count=0;
+ npos = nxs*nzs;
+
+ if (verbose) vmess("nxs: %d, nzs: %d",nxs,nzs);
+
+ nshots = 0;
+ getFileInfo(fshot, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax1, &sclsxgx, &ntraces);
+
+ if (dxrcv < 0) dxrcv=dx;
+ if (dzrcv < 0) dzrcv=dx;
+
+ // ngath zijn het aantal schoten
+ shotdata = (float *)malloc(nt*nx*nshots*sizeof(float));
+ hdr_shot = (segy *)calloc(nx*nshots,sizeof(segy));
+
+ fp_shot = fopen(fshot,"r");
+ if (fp_shot == NULL) {
+ verr("Could not open file");
+ }
+ vmess("nt: %d nx: %d nshots: %d",nt,nx,nshots);
+ fclose(fp_shot);
+ readSnapData(fshot, &shotdata[0], &hdr_shot[0], nshots, nx, nt, 0, nx, 0, nt);
+
+
+ hdr_out = (segy *)calloc(nxs,sizeof(segy));
+ Ghom = (float *)malloc(nt*npos*sizeof(float));
+
+ if (scheme==2) {
+ vmess("Classical representation");
+ shotdata_jkz = (float *)malloc(nt*nx*nshots*sizeof(float));
+ for (ix = 0; ix < nx; ix++) {
+ for (it = 0; it < nt; it++) {
+ shotdata_jkz[ix*nt+it] = shotdata[ix*nt+it];
+ }
+ }
+ conjugate(shotdata_jkz, nt, nx, dt);
+ conjugate(shotdata, nt, nx, dt);
+ depthDiff(shotdata_jkz, nt, nx, dt, dx, fmin, fmax, cp, 1);
+ if (verbose) vmess("Applied jkz to source data");
+ }
+ else if (scheme==0) {
+ vmess("Marchenko representation");
+ }
+ else if (scheme==1) {
+ vmess("Marchenko representation with multiple sources");
+ }
+ else if (scheme==3) {
+ vmess("Marchenko representation with multiple shot gathers");
+ }
+
+#pragma omp parallel default(shared) \
+ private(ix,it,is,indata, hdr_in,fins,fin2,fp_in,conv,ig,conv2,tmp1,tmp2)
+{
+ indata = (float *)malloc(nt*nx*nxs*sizeof(float));
+ hdr_in = (segy *)calloc(nx*nxs,sizeof(segy));
+ conv = (float *)calloc(nx*nt,sizeof(float));
+ conv2 = (float *)calloc(nx*nt,sizeof(float));
+ if (scheme==2) {
+ tmp1 = (float *)calloc(nx*nt,sizeof(float));
+ tmp2 = (float *)calloc(nx*nt,sizeof(float));
+ }
+#pragma omp for
+ for (ir = 0; ir < nzs; ir++) {
+ sprintf(fins,"z%d",ir*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin2, "r");
+ if (fp_in == NULL) {
+ verr("Danger Will Robinson");
+ }
+ fclose(fp_in);
+ readSnapData(fin2, &indata[0], &hdr_in[0], nxs, nx, nt, 0, nx, 0, nt);
+ for (is=0;is<nxs;is++) {
+ if (scheme==0) { //Marchenko representation
+ depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(shotdata, &indata[is*nx*nt], conv, nx, nt, dt, -2);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -2);
+ for (ix=0; ix<nx; ix++) {
+ for (it=0; it<nt/2; it++) {
+ Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv[ix*nt+it]/rho;
+ Ghom[it*nxs*nzs+is*nzs+ir] += conv[ix*nt+(it+nt/2)]/rho;
+ }
+ }
+ }
+ else if (scheme==1) { //Marchenko representation with multiple sources
+ depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(shotdata, &indata[is*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ for (ix=0; ix<nx; ix++) {
+ for (it=0; it<nt/2; it++) {
+ Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += 2*conv[ix*nt+it]/rho;
+ Ghom[it*nxs*nzs+is*nzs+ir] += 2*conv[ix*nt+(it+nt/2)]/rho;
+ }
+ }
+ }
+ else if (scheme==2) { //classical representation
+ convol(&indata[is*nx*nt], shotdata_jkz, tmp1, nx, nt, dt, 0);
+ depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(&indata[is*nx*nt], shotdata, tmp2, nx, nt, dt, 0);
+ //corr(&indata[is*nx*nt], shotdata, tmp2, nx, nt, dt, 0);
+ for (ix = 0; ix < nx; ix++) {
+ for (it = 0; it < nt; it++) {
+ conv[ix*nt+it] = tmp2[ix*nt+it]+tmp1[ix*nt+it];
+ }
+ }
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ for (ix=0; ix<nx; ix++) {
+ for (it=0; it<nt/2; it++) {
+ Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv[ix*nt+it]/rho;
+ Ghom[it*nxs*nzs+is*nzs+ir] += conv[ix*nt+(it+nt/2)]/rho;
+ }
+ }
+ }
+ if (scheme==3) { //Marchenko representation with multiple shot gathers
+ depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ for (ig=0; ig<nshots; ig++) {
+ convol(&shotdata[ig*nx*nt], &indata[is*nx*nt], conv, nx, nt, dt, -2);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -2);
+ shift_num = ig*((int)(shift/dt));
+ for (ix = 0; ix < nx; ix++) {
+ for (it = nt/2+1; it < nt; it++) {
+ conv[ix*nt+it] = 0.0;
+ }
+ for (it = shift_num; it < nt; it++) {
+ conv2[ix*nt+it] = conv[ix*nt+it-shift_num];
+ }
+ for (it = 0; it < shift_num; it++) {
+ conv2[ix*nt+it] = conv[ix*nt+nt-shift_num+it];
+ }
+ }
+ for (ix=0; ix<nx; ix++) {
+ Ghom[(-1+nt/2)*nxs*nzs+is*nzs+ir] += conv2[ix*nt+nt-1]/rho;
+ for (it=0; it<nt/2; it++) {
+ Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv2[ix*nt+it]/rho;
+ //Ghom[it*nxs*nzs+is*nzs+ir] += conv2[ix*nt+(it+nt/2)]/rho;
+ }
+ }
+ }
+ }
+ }
+
+ count+=1;
+ if (verbose) vmess("Creating Homogeneous Green's function at depth %d from %d depths",count,nzs);
+ }
+ free(conv); free(indata); free(hdr_in); free(conv2);
+ if (scheme==2) {
+ free(tmp1);free(tmp2);
+ }
+}
+ free(shotdata);
+
+ if (verbose) vmess("nxs: %d nxz: %d f1: %.7f",nxs,nzs,f1);
+
+ ntshift=0;
+
+ if (ntmax > 0) {
+ if (ntmax < nt) {
+ ntshift = (nt-ntmax)/2;
+ if (verbose) vmess("Time shifted %d samples",ntshift);
+ nt=ntmax;
+ }
+ else {
+ if (verbose) vmess("Max time samples larger than original samples");
+ }
+ }
+
+ fp_out = fopen(fout, "w+");
+
+ for (ir = 0; ir < nt; ir++) {
+ for (ix = 0; ix < nxs; ix++) {
+ hdr_out[ix].fldr = ir+1;
+ hdr_out[ix].tracl = ir*nxs+ix+1;
+ hdr_out[ix].tracf = ix+1;
+ hdr_out[ix].scalco = hdr_shot[0].scalco;
+ hdr_out[ix].scalel = hdr_shot[0].scalel;
+ hdr_out[ix].sdepth = hdr_shot[0].sdepth;
+ hdr_out[ix].trid = 1;
+ hdr_out[ix].ns = nzs;
+ hdr_out[ix].trwf = nxs;
+ hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
+ hdr_out[ix].f1 = f1;
+ hdr_out[ix].f2 = f2/1000;
+ hdr_out[ix].dt = dt*(1E6);
+ hdr_out[ix].d1 = dzrcv;
+ hdr_out[ix].d2 = dxrcv;
+ hdr_out[ix].sx = hdr_shot[0].sx;
+ hdr_out[ix].gx = (int)roundf(f2 + (ix*hdr_out[ix].d2)*1000.0);
+ hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &Ghom[(ir+ntshift)*nxs*nzs], hdr_out, nzs, nxs);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ return 0;
+}
+
+void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt, int shift)
+{
+ int i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccon, tmp;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccon = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (ccon == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
+ rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
+
+ /* apply convolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccon[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ n = nrec*nfreq;
+ for (j = 0; j < n; j++) {
+ *qr = (*p2r**p1r-*p2i**p1i);
+ *qi = (*p2r**p1i+*p2i**p1r);
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ if (shift==1) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+ tau = dt*(nsam/2);
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = ccon[j*nfreq+i].r*cos(om*tau) + ccon[j*nfreq+i].i*sin(om*tau);
+ tmp.i = ccon[j*nfreq+i].i*cos(om*tau) - ccon[j*nfreq+i].r*sin(om*tau);
+ ccon[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+ if (shift==-2) {
+ for (j = 0; j < nrec; j++) {
+ for (i = 0; i < nfreq; i++) {
+ ccon[j*nfreq+i].r = ccon[j*nfreq+i].i;
+ ccon[j*nfreq+i].i = 0.0;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/((float)(optn));
+ crmfft(&ccon[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,con,nsam);
+
+ free(ccon);
+ free(rdata1);
+ free(rdata2);
+ return;
+}
+
+void pad_data(float *data, int nsam, int nrec, int nsamout, float *datout)
+{
+ int it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsamout+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsamout+it]=0.0;
+ }
+}
+
+void scl_data(float *data, int nsam, int nrec, float scl, float *datout, int nsamout)
+{
+ int it,ix;
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsamout ; it++)
+ datout[ix*nsamout+it] = scl*data[ix*nsam+it];
+ }
+}
+
+void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt, int shift)
+{
+ int i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccov, tmp;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccov = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (ccov == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
+ rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
+
+ /* apply correlation */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccov[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ n = nrec*nfreq;
+ for (j = 0; j < n; j++) {
+ *qr = (*p1r * *p2r + *p1i * *p2i);
+ *qi = (*p1i * *p2r - *p1r * *p2i);
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ /* shift t=0 to middle of time window (nsam/2)*/
+ if (shift) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+ tau = dt*(nsam/2);
+
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = ccov[j*nfreq+i].r*cos(om*tau) + ccov[j*nfreq+i].i*sin(om*tau);
+ tmp.i = ccov[j*nfreq+i].i*cos(om*tau) - ccov[j*nfreq+i].r*sin(om*tau);
+ ccov[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&ccov[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,cov,nsam);
+
+ free(ccov);
+ free(rdata1);
+ free(rdata2);
+ return;
+}
+
+void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt)
+{
+ int 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 = (int)MIN((omin/deltom), (nfreq));
+ iomin = MAX(iomin, 1);
+ iomax = MIN((int)(omax/deltom), (nfreq));
+
+ cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ 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;
+ }
+ }
+ }
+ 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, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt)
+{
+ int 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 = (int)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((int)(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((int)(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((int)(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, int nsam, int nrec, int nsamout, int nrecout, float *datout)
+{
+ int 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, int nsam, int nrec, float dt)
+{
+ int 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;
+}
+
+void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax)
+{
+ int optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, tom, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (int)MIN((omin/deltom), (nfreq));
+ iomax = MIN((int)(omax/deltom), (nfreq));
+
+ cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ tom = om*shift;
+ cdatascl[ix*nfreq+iom].r = cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom);
+ cdatascl[ix*nfreq+iom].i = cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom);
+ }
+ }
+ free(cdata);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
diff --git a/utils/Makefile b/utils/Makefile
index 7d84e9b..d828166 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
+ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d MuteSnap combine combine_induced reshape_su HomG Snap2Shot
SRCM = \
makemod.c \
@@ -133,6 +133,68 @@ SRCT = ftr1d.c \
docpkge.c \
getpars.c
+SRCMS = MuteSnap.c \
+ getFileInfo.c \
+ writeData.c \
+ verbosepkg.c \
+ getpars.c \
+ wallclock_time.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
+
+SRCCO = combine.c \
+ getFileInfo.c \
+ writeData.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
+
+
+SRCCI = combine_induced.c \
+ getFileInfo.c \
+ writeData.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
+
+SRCRS = reshape_su.c \
+ getFileInfo.c \
+ writeData.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
+
+SRCHG = HomG.c \
+ getFileInfo.c \
+ readData.c \
+ writeData.c \
+ verbosepkg.c \
+ getpars.c \
+ wallclock_time.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
+
+SRCSS = snap2shot.c \
+ getFileInfo3D.c \
+ writeData3D.c \
+ verbosepkg.c \
+ getpars.c \
+ wallclock_time.c \
+ atopkge.c \
+ docpkge.c \
+ name_ext.c \
+ readSnapData3D.c
+
OBJM = $(SRCM:%.c=%.o)
makemod: $(OBJM)
@@ -188,7 +250,37 @@ OBJT = $(SRCT:%.c=%.o)
ftr1d: $(OBJT)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o ftr1d $(OBJT) $(LIBS)
-install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d
+OBJMS = $(SRCMS:%.c=%.o)
+
+MuteSnap: $(OBJMS)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o MuteSnap $(OBJMS) $(LIBS)
+
+OBJCO = $(SRCCO:%.c=%.o)
+
+combine: $(OBJCO)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine $(OBJCO) $(LIBS)
+
+OBJCI = $(SRCCI:%.c=%.o)
+
+combine_induced: $(OBJCI)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine_induced $(OBJCI) $(LIBS)
+
+OBJRS = $(SRCRS:%.c=%.o)
+
+reshape_su: $(OBJRS)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o reshape_su $(OBJRS) $(LIBS)
+
+OBJHG = $(SRCHG:%.c=%.o)
+
+HomG: $(OBJHG)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o HomG $(OBJHG) $(LIBS)
+
+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
cp makemod $B
cp makewave $B
cp extendModel $B
@@ -200,9 +292,15 @@ install: makemod makewave extendModel fconv correigen green green3D basop syn2d
cp syn2d $B
cp mat2su $B
cp ftr1d $B
+ cp MuteSnap $B
+ cp combine $B
+ cp combine_induced $B
+ cp reshape_su $B
+ cp HomG $B
+ cp Snap2Shot $B
clean:
- rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT)
+ rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) MuteSnap $(OBJMS) combine $(OBJCO) reshape_su $(OBJRS) combine_induced $(OBJCI) HomG $(OBJHG) Snap2Shot $(OBJSS)
realclean: clean
- rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d
+ rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d $B/MuteSnap $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/Snap2Shot
diff --git a/utils/MuteSnap.c b/utils/MuteSnap.c
new file mode 100644
index 0000000..b90681e
--- /dev/null
+++ b/utils/MuteSnap.c
@@ -0,0 +1,227 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+double wallclock_time(void);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+int topdet(float *data, int nt);
+int farrdet(float *data, int nt, float tol);
+
+char *sdoc[] = {
+" ",
+" HomG - Calculate a Homogeneous Green's function ",
+" ",
+" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" fhom= .................... File containing the snapshot data that will be muted",
+" fsnap= ................... File containing the snapshot data that will determine the mute window",
+" ",
+" Optional parameters: ",
+" ",
+" fout=out.su .............. Filename of the output",
+" shift=5 .................. Shift from the maximum",
+" smooth=5 ................. Length of smoothing taper",
+" mode=0 ................... Determine first arrival by maximum (mode=0), first event above tol (mode=1) or by raytime (mode=2)",
+" tol=1 .................... Tolerance for the determination of first arrival if mode=1",
+" fray ..................... File containing the raytimes of the first arrivals",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_snap, *fp_hom, *fp_out;
+ char *fhom, *fsnap, *fout, *fray;
+ float *homdata, *snapdata, *outdata, *rtrace, *costaper, scl, tol, *timeval, dt;
+ float dx, dz, z0, x0, xmin, xmax, sclsxgx, f1, f2, dxrcv, dzrcv, dxpos, offset;
+ int nt, nts, nx, nxs, nxh, ntraces, ret, ix, it, is, ir, nzs, nzh, nz, ht, indrcv, shift;
+ int rmt, smooth, mode, nzh1, nzs1, nxh1, nxs1, nts1, nt1;
+ segy *hdr_hom, *hdr_snap, *hdrs_mute;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("fhom", &fhom)) fhom = NULL;
+ if (!getparstring("fsnap", &fsnap)) fsnap = NULL;
+ if (!getparstring("fout", &fout)) fout = "out.su";
+ if (!getparstring("fray", &fray)) fray = NULL;
+ if (!getparint("shift", &shift)) shift = 5;
+ if (!getparint("smooth", &smooth)) smooth = 5;
+ if (!getparint("mode", &mode)) mode = 0;
+ if (!getparfloat("tol", &tol)) tol = 5;
+ if (fhom == NULL) verr("Incorrect G_hom input");
+ if (fsnap == NULL) verr("Incorrect snapshot input");
+ if (mode == 2) {
+ if (fray == NULL) verr("No filename for raytimes given");
+ }
+ if (!getparint("nxs1", &nxs1)) nxs1 = 0;
+ if (!getparint("nxh1", &nxh1)) nxh1 = 0;
+ if (!getparint("nzs1", &nzs1)) nzs1 = 0;
+ if (!getparint("nzh1", &nzh1)) nzh1 = 0;
+ if (!getparint("nts1", &nts1)) nts1 = 0;
+ if (!getparint("nt1", &nt1)) nt1 = 0;
+
+ if (smooth) {
+ costaper = (float *)malloc(smooth*sizeof(float));
+ scl = M_PI/((float)smooth);
+ for (is=0; is<smooth; is++) {
+ costaper[is] = 0.5*(1.0+cos((is+1)*scl));
+ }
+ }
+
+ getFileInfo(fsnap, &nzs, &nxs, &nts, &dz, &dx, &z0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+ getFileInfo(fhom, &nzh, &nxh, &nt, &dz, &dx, &z0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+
+ if (nxh1 != 0) nxh = nxh1;
+ if (nxs1 != 0) nxs = nxs1;
+ if (nzh1 != 0) nzh = nzh1;
+ if (nzs1 != 0) nzs = nzs1;
+ if (nt1 != 0) nt = nt1;
+ if (nts1 != 0) nts = nts1;
+
+ if (nzs != nzh || nxs != nxh) {
+ verr("sampling in x or z direction is incorrect, nzs=%d nzh=%d, nxs=%d nxh=%d",nzs,nzh,nxs,nxh);
+ }
+
+ vmess("nzs=%d nzh=%d, nxs=%d nxh=%d, nts=%d nt=%d",nzs,nzh,nxs,nxh,nts,nt);
+
+ nz = nzh;
+ nx = nxh;
+
+ snapdata = (float *)malloc(nz*nx*nts*sizeof(float));
+ hdr_snap = (segy *)calloc(nx*nts,sizeof(segy));
+ homdata = (float *)malloc(nz*nx*nt*sizeof(float));
+ hdr_hom = (segy *)calloc(nx*nt,sizeof(segy));
+ ht = (int)ceil(nt/2);
+ rtrace = (float *)malloc(nts*sizeof(float));
+
+ if (mode != 2) {
+ readSnapData(fsnap, &snapdata[0], &hdr_snap[0], nts, nx, nz, 0, nx, 0, nz);
+ vmess("Read Snapshot data");
+ }
+ readSnapData(fhom, &homdata[0], &hdr_hom[0], nt, nx, nz, 0, nx, 0, nz);
+ vmess("Read G_hom");
+
+ if (mode == 0) {
+ vmess("First arrival determined through maximum");
+ }
+ else if (mode == 1) {
+ vmess("First arrival determined through tolerance (=%.4f)",tol);
+ }
+ else if (mode == 2) {
+ vmess("First arrival determined through raytimes");
+ fp_snap = fopen(fray,"r");
+ if (fp_snap == NULL) {
+ verr("Could not open file");
+ }
+ fclose(fp_snap);
+ hdrs_mute = (segy *) calloc(nz,sizeof(segy));
+ timeval = (float *)calloc(nz*nx,sizeof(float));
+ readSnapData(fray, timeval, hdrs_mute, nz, 1, nx, 0, 1, 0, nx);
+ dt = hdr_hom[0].dt/1E6;
+ }
+
+ for (ir = 0; ir < nz; ir++) {
+ for (is = 0; is < nx; is++) {
+ for (it = 0; it < nts; it++) {
+ rtrace[it] = snapdata[it*nxs*nzs+is*nzs+ir];
+ }
+ if (mode == 0) {
+ indrcv = topdet(&rtrace[0],nts);
+ }
+ else if (mode == 1) {
+ indrcv = farrdet(&rtrace[0],nts,tol);
+ }
+ else if (mode == 2) {
+ indrcv = (int)roundf(timeval[ir*nx+is]/dt);
+ }
+ rmt = MIN(nt-indrcv,indrcv)-shift;
+ for (it = ht-rmt+1; it < ht; it++) {
+ if (it-(ht-rmt) < smooth) {
+ homdata[it*nxs*nzs+is*nzs+ir] *= costaper[it-(ht-rmt)];
+ }
+ else {
+ homdata[it*nxs*nzs+is*nzs+ir] = 0.0;
+ }
+ }
+ for (it = ht; it < ht+rmt; it++) {
+ if (it-(ht+rmt)+smooth > 0) {
+ homdata[it*nxs*nzs+is*nzs+ir] *= costaper[smooth-(it-(ht+rmt)+smooth)];
+ }
+ else {
+ homdata[it*nxs*nzs+is*nzs+ir] = 0.0;
+ }
+ }
+ }
+ vmess("Muting Homogeneous Green's function at depth %d from %d depths",ir+1,nzs);
+ }
+ free(snapdata);free(hdr_snap);
+
+ fp_out = fopen(fout, "w+");
+
+ for (ir = 0; ir < nt; ir++) {
+ ret = writeData(fp_out, &homdata[ir*nxs*nzs], &hdr_hom[ir*nx], nz, nx);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ vmess("Wrote Data");
+ return 0;
+}
+
+int topdet(float *data, int nt)
+{
+ int it,ind;
+ float maxval;
+
+ maxval = data[0];
+ ind = 0;
+
+ for (it = 1; it < nt; it++) {
+ if (fabs(data[it]) > maxval) {
+ maxval = data[it];
+ ind = it;
+ }
+ }
+
+ return ind;
+}
+
+int farrdet(float *data, int nt, float tol)
+{
+ int it,ind;
+
+ ind = 0;
+
+ for (it = 0; it < nt; it++) {
+ if (fabs(data[it]) > tol) {
+ ind = it;
+ break;
+ }
+ }
+
+ return ind;
+}
diff --git a/utils/combine.c b/utils/combine.c
new file mode 100755
index 0000000..a93a881
--- /dev/null
+++ b/utils/combine.c
@@ -0,0 +1,198 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+double wallclock_time(void);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+
+char *sdoc[] = {
+" ",
+" combine - Combine results into a single result ",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. File containing the first data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first file",
+" dnumb= ................... integer number of increment in files",
+" nzmax= ................... Maximum number of files read",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_out;
+ char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
+ float *indata, *outdata, *rtrace, fz, fx;
+ float dt, dx, t0, x0, xmin, xmax, sclsxgx, dt2, dx2, t02, x02, xmin2, xmax2, sclsxgx2, dxrcv, dzrcv;
+ int nshots, nt, nx, ntraces, nshots2, nt2, nx2, ntraces2, ix, it, is, iz, pos, ifile, file_det, nxs, nzs;
+ int xcount, numb, dnumb, ret, nzmax, transpose, verbose;
+ segy *hdr_in, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparint("numb", &numb)) numb=0;
+ if (!getparint("dnumb", &dnumb)) dnumb=0;
+ if (!getparint("nzmax", &nzmax)) nzmax=0;
+ if (!getparint("verbose", &verbose)) verbose=0;
+ if (!getparint("transpose", &transpose)) transpose=0;
+ if (fin == NULL) verr("Incorrect downgoing input");
+
+ if (dnumb < 1) dnumb = 1;
+
+ sprintf(numb1,"%d",numb);
+
+ ptr = strstr(fin,numb1);
+ pos = ptr - fin + 1;
+
+ sprintf(fbegin,"%*.*s", pos-1, pos-1, fin);
+ sprintf(fend,"%s", fin+pos);
+
+ file_det = 1;
+ nzs=0;
+
+ while (file_det) {
+ sprintf(fins,"%d",nzs*dnumb+numb);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ if (nzs == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzs == 1) {
+ vmess("1 file detected");
+ }
+ else {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzs++;
+ if (nzmax!=0 && nzs == nzmax) {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+
+ sprintf(fins,"%d",numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+
+ sprintf(fins,"%d",numb+dnumb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt2, &nx2, &nshots2, &dt2, &dx2, &t02, &x02, &xmin2, &xmax2, &sclsxgx2, &ntraces2);
+
+ dxrcv=dx*1000;
+ dzrcv=t02-t0;
+
+ if (nshots==0) nshots=1;
+ nxs = ntraces;
+
+
+
+ // ngath zijn het aantal schoten
+ hdr_out = (segy *)calloc(nxs*nt,sizeof(segy));
+ outdata = (float *)calloc(nxs*nzs*nt,sizeof(float));
+ hdr_in = (segy *)calloc(nxs*nt,sizeof(segy));
+ indata = (float *)calloc(nxs*nt,sizeof(float));
+
+ readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ nshots = hdr_in[nxs-1].fldr;
+ if (transpose==0) {
+ nxs = hdr_in[nxs-1].tracf;
+ }
+ else {
+ nxs = hdr_in[nxs-1].ns;
+ }
+
+ for (iz = 0; iz < nzs; iz++) {
+ if (verbose) vmess("Depth:%d out of %d",iz+1,nzs);
+ sprintf(fins,"%d",iz*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin2, "r");
+ if (fp_in == NULL) {
+ verr("Error opening file");
+ }
+ fclose(fp_in);
+ if (transpose==0) {
+ readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ }
+ else {
+ readSnapData(fin2, &indata[0], &hdr_in[0], nshots, 1, nxs, 0, 1, 0, nxs);
+ }
+ if (iz==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
+ if (iz==1) dzrcv=hdr_in[0].f1-fz;
+ for (ix = 0; ix < nxs; ix++) {
+ for (is = 0; is < nshots; is++) {
+ outdata[is*nxs*nzs+ix*nzs+iz] = indata[is*nxs+ix];
+ }
+ }
+ }
+ free(indata);
+
+ fp_out = fopen(fout, "w+");
+
+ for (is = 0; is < nshots; is++) {
+ for (ix = 0; ix < nxs; ix++) {
+ hdr_out[ix].fldr = is+1;
+ hdr_out[ix].tracl = is*nxs+ix+1;
+ hdr_out[ix].tracf = ix+1;
+ hdr_out[ix].scalco = -1000;
+ hdr_out[ix].scalel = -1000;
+ hdr_out[ix].sdepth = hdr_in[0].sdepth;
+ hdr_out[ix].trid = 1;
+ hdr_out[ix].ns = nzs;
+ hdr_out[ix].trwf = nxs;
+ hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
+ hdr_out[ix].f1 = fz;
+ hdr_out[ix].f2 = fx;
+ hdr_out[ix].dt = dt*(1E6);
+ hdr_out[ix].d1 = dzrcv;
+ hdr_out[ix].d2 = dxrcv;
+ hdr_out[ix].sx = (int)roundf(fx + (ix*hdr_out[ix].d2));
+ hdr_out[ix].gx = (int)roundf(fx + (ix*hdr_out[ix].d2));
+ hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &outdata[is*nxs*nzs], hdr_out, nzs, nxs);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ return 0;
+}
+
diff --git a/utils/combine_induced.c b/utils/combine_induced.c
new file mode 100755
index 0000000..0b7601e
--- /dev/null
+++ b/utils/combine_induced.c
@@ -0,0 +1,208 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+double wallclock_time(void);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+
+char *sdoc[] = {
+" ",
+" combine_induced - Combine induced seismicity results together ",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. File containing the first data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first file",
+" dnumb= ................... integer number of increment in files",
+" nzmax= ................... Maximum number of files read",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_out;
+ char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100];
+ float *indata, *outdata, *rtrace, fz, fx, shift, dtshift, dt_time;
+ float dt, dx, t0, x0, xmin, xmax, sclsxgx, dt2, dx2, t02, x02, xmin2, xmax2, sclsxgx2, dxrcv, dzrcv;
+ int nshots, nt, nx, ntraces, nshots2, nt2, nx2, ntraces2, ix, it, is, iz, pos, ifile, file_det, nxs, nzs;
+ int xcount, numb, dnumb, ret, nzmax, verbose, nshot_out, ishift, nshift;
+ segy *hdr_in, *hdr_bin, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparint("numb", &numb)) numb=0;
+ if (!getparint("dnumb", &dnumb)) dnumb=0;
+ if (!getparint("nzmax", &nzmax)) nzmax=0;
+ if (!getparint("verbose", &verbose)) verbose=0;
+ if (!getparint("nshift", &nshift)) nshift=0;
+ if (!getparfloat("shift", &shift)) shift=0.0;
+ if (!getparfloat("dtshift", &dtshift)) dtshift=0.0;
+ if (!getparfloat("dt_time", &dt_time)) dt_time=0.004;
+ if (fin == NULL) verr("Incorrect downgoing input");
+
+ if (dnumb < 1) dnumb = 1;
+
+ sprintf(numb1,"%d",numb);
+
+ ptr = strstr(fin,numb1);
+ pos = ptr - fin + 1;
+
+ sprintf(fbegin,"%*.*s", pos-1, pos-1, fin);
+ sprintf(fend,"%s", fin+pos);
+
+ file_det = 1;
+ nzs=0;
+
+ while (file_det) {
+ sprintf(fins,"%d",nzs*dnumb+numb);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ if (nzs == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzs == 1) {
+ vmess("1 file detected");
+ }
+ else {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzs++;
+ if (nzmax!=0 && nzs == nzmax) {
+ vmess("%d files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+
+ sprintf(fins,"%d",numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+
+ sprintf(fins,"%d",numb+dnumb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ nshots = 0;
+ getFileInfo(fin2, &nt2, &nx2, &nshots2, &dt2, &dx2, &t02, &x02, &xmin2, &xmax2, &sclsxgx2, &ntraces2);
+
+ dxrcv=dx*1000;
+ dzrcv=t02-t0;
+
+ if (nshots==0) nshots=1;
+ nxs = ntraces;
+
+ // ngath zijn het aantal schoten
+ hdr_bin = (segy *)calloc(nxs,sizeof(segy));
+ indata = (float *)calloc(nxs*nt,sizeof(float));
+
+ readSnapData(fin2, &indata[0], &hdr_bin[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ nshots = hdr_bin[nxs-1].fldr;
+ nxs = hdr_bin[nxs-1].tracf;
+
+ nshot_out = nshots/2;
+ free(indata);
+
+ hdr_out = (segy *)calloc(nshot_out*nxs,sizeof(segy));
+ outdata = (float *)calloc(nshot_out*nxs*nt,sizeof(float));
+
+ vmess("nshot_out:%d nxs=%d nt:%d",nshot_out,nxs,nt);
+
+#pragma omp parallel default(shared) \
+ private(indata,iz,hdr_in,fins,fin2,fp_in,is,ix,it,ishift)
+{
+ indata = (float *)calloc(ntraces*nt,sizeof(float));
+ hdr_in = (segy *)calloc(ntraces,sizeof(segy));
+
+#pragma omp for
+ for (iz = 0; iz < nzs; iz++) {
+ if (verbose) vmess("Depth:%d out of %d",iz+1,nzs);
+ sprintf(fins,"%d",iz*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin2, "r");
+ if (fp_in == NULL) {
+ verr("Error opening file");
+ }
+ fclose(fp_in);
+ readSnapData(fin2, &indata[0], &hdr_in[0], nshots, nxs, nt, 0, nxs, 0, nt);
+ if (iz==0) fz=hdr_in[0].f1; fx=hdr_in[0].f2;
+ if (iz==1) dzrcv=hdr_in[0].f1-fz;
+ //ishift = (int)((shift+(dtshift*((float)iz)))/dt_time);
+ ishift = nshift*iz;
+ if (verbose) vmess("Shifting %d timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)iz)));
+ for (is = ishift; is < nshot_out; is++) {
+ for (ix = 0; ix < nxs; ix++) {
+ for (it = 0; it < nt; it++) {
+ outdata[is*nxs*nt+ix*nt+it] += indata[(is-ishift+(nshots/2))*nxs*nt+ix*nt+it];
+ }
+ }
+ }
+ }
+ free(indata);free(hdr_in);
+}
+
+ fp_out = fopen(fout, "w+");
+
+ for (is = 0; is < nshot_out; is++) {
+ for (ix = 0; ix < nxs; ix++) {
+ hdr_out[ix].fldr = is+1;
+ hdr_out[ix].tracl = is*nxs+ix+1;
+ hdr_out[ix].tracf = ix+1;
+ hdr_out[ix].scalco = -1000;
+ hdr_out[ix].scalel = -1000;
+ hdr_out[ix].sdepth = hdr_bin[0].sdepth;
+ hdr_out[ix].trid = 1;
+ hdr_out[ix].ns = nt;
+ hdr_out[ix].trwf = nxs;
+ hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
+ hdr_out[ix].f1 = fz;
+ hdr_out[ix].f2 = fx;
+ hdr_out[ix].dt = dt_time*(1E6);
+ hdr_out[ix].d1 = dzrcv;
+ hdr_out[ix].d2 = dxrcv;
+ hdr_out[ix].sx = (int)roundf(fx + (ix*hdr_out[ix].d2));
+ hdr_out[ix].gx = (int)roundf(fx + (ix*hdr_out[ix].d2));
+ hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &outdata[is*nxs*nt], hdr_out, nt, nxs);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ return 0;
+}
+
diff --git a/utils/getFileInfo3D.c b/utils/getFileInfo3D.c
new file mode 100644
index 0000000..fd00cee
--- /dev/null
+++ b/utils/getFileInfo3D.c
@@ -0,0 +1,244 @@
+#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"
+
+#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 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)
+{
+ FILE *fp;
+ size_t nread, data_sz;
+ off_t bytes, ret, trace_sz, ntraces;
+ long sx_shot, sy_shot, gx_start, gx_end, gy_start, gy_end, itrace, one_shot, igath, end_of_file, fldr_shot;
+ long verbose=1, igy, nsx, nsy;
+ float scl, *trace, dxsrc, dxrcv, dysrc, dyrcv;
+ segy hdr;
+
+ if (filename == NULL) { /* read from input pipe */
+ *n1=0;
+ *n2=0;
+ *n3=0;
+ return -1; /* Input pipe */
+ }
+ else fp = fopen( filename, "r" );
+ if (fp == NULL) verr("File %s does not exist or cannot be opened", filename);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ ret = fseeko( fp, 0, SEEK_END );
+ if (ret<0) perror("fseeko");
+ bytes = ftello( fp );
+
+ *n1 = hdr.ns;
+ if ( (hdr.trid == 1) && (hdr.dt != 0) ) {
+ *d1 = ((float) hdr.dt)*1.e-6;
+ *f1 = ((float) hdr.delrt)/1000.;
+ }
+ else {
+ *d1 = hdr.d1;
+ *f1 = hdr.f1;
+ }
+ *f2 = hdr.f2;
+ *f3 = hdr.gy;
+
+ data_sz = sizeof(float)*(*n1);
+ trace_sz = sizeof(float)*(*n1)+TRCBYTES;
+ ntraces = (long) (bytes/trace_sz);
+
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+
+ *sclsxgxsygy = scl;
+ /* check to find out number of traces in shot gather */
+
+ one_shot = 1;
+ itrace = 1;
+ igy = 1;
+ fldr_shot = hdr.fldr;
+ sx_shot = hdr.sx;
+ sy_shot = hdr.sy;
+ gx_start = hdr.gx;
+ gy_start = hdr.gy;
+ gy_end = gy_start;
+ trace = (float *)malloc(hdr.ns*sizeof(float));
+ fseeko( fp, TRCBYTES, SEEK_SET );
+
+ while (one_shot) {
+ nread = fread( trace, sizeof(float), hdr.ns, fp );
+ assert (nread == hdr.ns);
+ if (hdr.gy != gy_end) {
+ gy_end = hdr.gy;
+ igy++;
+ }
+ gx_end = hdr.gx;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread==0) break;
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr) ) break;
+ itrace++;
+ }
+
+ if (itrace>1) {
+ *n2 = itrace/igy;
+ *n3 = igy;
+ if (*n2>1) {
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ }
+ else {
+ dxrcv = 1.0/scl;
+ }
+ if (*n3>1) {
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ }
+ else {
+ dyrcv = 1.0/scl;
+ }
+ *d2 = fabs(dxrcv)*scl;
+ *d3 = fabs(dyrcv)*scl;
+ if (NINT(dxrcv*1e3) != NINT(fabs(hdr.d2)*1e3)) {
+ if (dxrcv != 0) *d2 = fabs(dxrcv)*scl;
+ else *d2 = hdr.d2;
+ }
+ }
+ else {
+ *n2 = MAX(hdr.trwf, 1);
+ *n3 = 1;
+ *d2 = hdr.d2;
+ *d3 = 1.0;
+ dxrcv = hdr.d2;
+ dyrcv = 0.0;
+ }
+
+/* check if the total number of traces (ntraces) is correct */
+
+/* expensive way to find out how many gathers there are */
+
+// fprintf(stderr, "ngath = %li dxrcv=%f d2=%f scl=%f \n", *ngath, dxrcv, *d2, scl);
+ if (*ngath == 0) {
+ *n2 = 0;
+ *n3 = 0;
+
+ end_of_file = 0;
+ one_shot = 1;
+ igath = 0;
+ fseeko( fp, 0, SEEK_SET );
+ dxrcv = *d2;
+ dyrcv = *d3;
+
+ while (!end_of_file) {
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { break; }
+ fldr_shot = hdr.fldr;
+ sx_shot = hdr.sx;
+ gx_start = hdr.gx;
+ gx_end = hdr.gx;
+ sy_shot = hdr.sy;
+ gy_start = hdr.gy;
+ gy_end = hdr.gy;
+
+ itrace = 1;
+ igy = 1;
+ while (one_shot) {
+ fseeko( fp, data_sz, SEEK_CUR );
+ if (hdr.gx != gx_end) dxrcv = MIN(dxrcv,abs(hdr.gx-gx_end));
+ if (hdr.gy != gy_end) {
+ igy++;
+ gy_end = hdr.gy;
+ dyrcv = MIN(dyrcv,abs(hdr.gy-gy_end));
+ }
+ gx_end = hdr.gx;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
+ itrace++;
+ }
+ if (itrace>1) {
+ *n2 = MAX(itrace/igy,*n2);
+ *n3 = igy;
+ if (*n2>1) {
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ }
+ else {
+ dxrcv = 1.0/scl;
+ }
+ if (*n3>1) {
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ }
+ else {
+ dyrcv = 1.0/scl;
+ }
+ dxsrc = (float)(hdr.sx - sx_shot)*scl;
+ dysrc = (float)(hdr.sy - sy_shot)*scl;
+ }
+ else {
+ *n2 = MAX(MAX(hdr.trwf, 1),*n2);
+ *n3 = 1;
+ *d2 = hdr.d2;
+ *d3 = 1.0;
+ dxrcv = hdr.d2/scl;
+ dyrcv = 1.0/scl;
+ }
+ if (verbose>1) {
+ fprintf(stderr," . Scanning shot %li (%li) with %li traces dxrcv=%.2f dxsrc=%.2f %li %li dyrcv=%.2f dysrc=%.2f %li %li\n",sx_shot,igath,itrace,dxrcv*scl,dxsrc,gx_end,gx_start,dyrcv*scl,dysrc,gy_end,gy_start);
+ }
+ if (itrace != 0) { /* end of shot record */
+ fseeko( fp, -TRCBYTES, SEEK_CUR );
+ igath++;
+ }
+ else {
+ end_of_file = 1;
+ }
+ }
+ *ngath = igath;
+ *d2 = dxrcv*scl;
+ *d3 = dyrcv*scl;
+ }
+ else {
+ /* read last trace header */
+
+ fseeko( fp, -trace_sz, SEEK_END );
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ *ngath = ntraces/((*n2)*(*n3));
+ }
+// *nxm = NINT((*xmax-*xmin)/dxrcv)+1;
+ *nxm = (long)ntraces;
+
+ fclose( fp );
+ free(trace);
+
+ return 0;
+}
+
+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)
+{
+ vmess("file %s contains", file);
+ vmess("*** n1 = %li n2 = %li n3 = %li ntftt=%li", n1, n2, n3, (long)optncr((int)n1));
+ vmess("*** d1 = %.5f d2 = %.5f d3 = %.5f", d1, d2, d3);
+ vmess("*** f1 = %.5f f2 = %.5f f3 = %.5f", f1, f2, f3);
+ vmess("*** fldr = %li sx = %li sy = %li", hdrs[0].fldr, hdrs[0].sx, hdrs[0].sy);
+
+ return 0;
+}
diff --git a/utils/readSnapData.c b/utils/readSnapData.c
new file mode 100755
index 0000000..f7b953f
--- /dev/null
+++ b/utils/readSnapData.c
@@ -0,0 +1,58 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+int optncr(int n);
+
+
+int readSnapData(char *filename, float *data, segy *hdrs, int nsnaps, int nx, int nz, int sx, int ex, int sz, int ez)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ int nt, it, ix, iz, dx, dz;
+ float *tmpdata;
+
+ tmpdata = (float *)malloc(nsnaps*nx*nz*sizeof(float));
+ /* 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;
+ }
+ //nread = fread(&hdr, 1, TRCBYTES, fp);
+ for (it = 0; it < nsnaps*nx; it++) {
+ nread = fread(&hdr, 1, TRCBYTES, fp);
+ if (nread != TRCBYTES) {
+ break;
+ }
+ assert(nread == TRCBYTES);
+ nread = fread(&tmpdata[it*nz], sizeof(float), nz, fp);
+ assert (nread == nz);
+ memcpy(&hdrs[it], &hdr, TRCBYTES);
+ }
+ dx = ex-sx;
+ dz = ez-sz;
+ for (iz = sz; iz < ez; iz++) {
+ for (ix = sx; ix < ex; ix++) {
+ for (it = 0; it < nsnaps; it++) {
+ data[it*dx*dz+(ix-sx)*dz+iz-sz]=tmpdata[it*nx*nz+ix*nz+iz];
+ }
+ }
+ }
+ fclose(fp);
+ free(tmpdata);
+ return 0;
+}
diff --git a/utils/readSnapData3D.c b/utils/readSnapData3D.c
new file mode 100644
index 0000000..503d490
--- /dev/null
+++ b/utils/readSnapData3D.c
@@ -0,0 +1,56 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+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)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ long nt, it, ix, iy, iz, dx, dy, dz;
+ float *tmpdata;
+
+ tmpdata = (float *)malloc(nsnaps*nx*ny*nz*sizeof(float));
+ /* 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;
+ }
+ //nread = fread(&hdr, 1, TRCBYTES, fp);
+ for (it = 0; it < nsnaps*nx*ny; it++) {
+ nread = fread(&hdr, 1, TRCBYTES, fp);
+ if (nread != TRCBYTES) {
+ break;
+ }
+ assert(nread == TRCBYTES);
+ nread = fread(&tmpdata[it*nz], sizeof(float), nz, fp);
+ assert (nread == nz);
+ memcpy(&hdrs[it], &hdr, TRCBYTES);
+ }
+ dx = ex-sx;
+ dy = ey-sy;
+ dz = ez-sz;
+ for (iz = sz; iz < ez; iz++) {
+ for (iy = sy; iy < ey; iy++) {
+ for (ix = sx; ix < ex; ix++) {
+ for (it = 0; it < nsnaps; it++) {
+ data[it*dy*dx*dz+(iy-sy)*dx*dz+(ix-sx)*dz+iz-sz]=tmpdata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
+ }
+ }
+ }
+ }
+ fclose(fp);
+ free(tmpdata);
+ return 0;
+}
diff --git a/utils/reshape_su.c b/utils/reshape_su.c
new file mode 100755
index 0000000..9c24040
--- /dev/null
+++ b/utils/reshape_su.c
@@ -0,0 +1,125 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
+
+char *sdoc[] = {
+" ",
+" reshape_su - interchange the 1st and 3rd dimension for SU file",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. File containing the first data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_out;
+ char *fin, *fout;
+ float *indata, *outdata;
+ float dt, dx, t0, x0, xmin, xmax, sclsxgx;
+ int nshots, nt, nx, ntraces, ix, it, is, ir, ret, verbose;
+ segy *hdr_in, *hdr_out, hdr;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparint("verbose", &verbose)) verbose = 0;
+ if (fin == NULL) verr("No input file specified");
+
+ nshots = 0;
+ getFileInfo(fin, &nt, &nx, &nshots, &dt, &dx, &t0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
+
+ fp_in = fopen( fin, "r" );
+ ret = fread( &hdr, 1, TRCBYTES, fp_in );
+ assert(ret == TRCBYTES);
+ fclose(fp_in);
+
+ if (nt==0) nt=hdr.ns;
+ if (nx==0) nx=hdr.trwf;nshots=ntraces/nx;
+ if (nshots==0) nshots=1;
+
+ vmess("nx:%d nt:%d nshots:%d ntraces:%d",nx,nt,nshots,ntraces);
+
+ // ngath zijn het aantal schoten
+ hdr_out = (segy *)calloc(nx,sizeof(segy));
+ outdata = (float *)calloc(nshots*nx*nt,sizeof(float));
+ hdr_in = (segy *)calloc(nshots*nx,sizeof(segy));
+ indata = (float *)calloc(nshots*nx*nt,sizeof(float));
+
+ readSnapData(fin, &indata[0], &hdr_in[0], nshots, nx, nt, 0, nx, 0, nt);
+
+ for (ir = 0; ir < nshots; ir++) {
+ for (is = 0; is < nx; is++) {
+ for (it = 0; it < nt; it++) {
+ outdata[it*nx*nshots+is*nshots+ir] = indata[ir*nx*nt+is*nt+it];
+ }
+ }
+ if (verbose) vmess("Reshaping shot %d out of %d shots",ir+1,nshots);
+ }
+ free(indata);
+
+ fp_out = fopen(fout, "w+");
+
+ for (is = 0; is < nt; is++) {
+ for (ix = 0; ix < nx; ix++) {
+ hdr_out[ix].fldr = is+1;
+ hdr_out[ix].tracl = is*nx+ix+1;
+ hdr_out[ix].tracf = ix+1;
+ hdr_out[ix].scalco = -1000;
+ hdr_out[ix].scalel = -1000;
+ hdr_out[ix].sdepth = hdr_in[0].sdepth;
+ hdr_out[ix].trid = 1;
+ hdr_out[ix].ns = nshots;
+ hdr_out[ix].trwf = nx;
+ hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
+ hdr_out[ix].f1 = -((float)(hdr_in[0].dt/1E6))*(nshots/2);
+ hdr_out[ix].f2 = hdr_in[0].f2;
+ hdr_out[ix].dt = hdr_in[0].dt;
+ hdr_out[ix].d1 = ((float)hdr_in[0].dt);
+ hdr_out[ix].d2 = (hdr_in[0].d2);
+ hdr_out[ix].sx = (int)roundf(hdr_out[ix].f2 + (ix*hdr_out[ix].d2));
+ hdr_out[ix].sx = hdr_in[ix].sx;
+ hdr_out[ix].gx = (int)roundf(hdr_out[ix].f2 + (ix*hdr_out[ix].d2));
+ hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &outdata[is*nx*nshots], hdr_out, nshots, nx);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ return 0;
+}
+
diff --git a/utils/snap2shot.c b/utils/snap2shot.c
new file mode 100644
index 0000000..4218a2b
--- /dev/null
+++ b/utils/snap2shot.c
@@ -0,0 +1,205 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
+ float *sclsxgxsygy, long *nxm);
+
+double wallclock_time(void);
+
+long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
+
+void name_ext(char *filename, char *extension);
+
+long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps, long nx, long ny, long nz, long sx, long ex, long sy, long ey, long sz, long ez);
+
+char *sdoc[] = {
+" ",
+" snap2shot - Reshape snapshot data to receiver data",
+" ",
+" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_snap= ............... File containing the snapshot data that will be muted",
+" file_rcv= ................ File containing the snapshot data that will determine the mute window",
+" ",
+" Optional parameters: ",
+" ",
+" fout=out.su .............. Filename of the output",
+" shift=5 .................. Shift from the maximum",
+" smooth=5 ................. Length of smoothing taper",
+" mode=0 ................... Determine first arrival by maximum (mode=0), first event above tol (mode=1) or by raytime (mode=2)",
+" tol=1 .................... Tolerance for the determination of first arrival if mode=1",
+" fray ..................... File containing the raytimes of the first arrivals",
+NULL};
+
+void main (int argc, char **argv)
+{
+ FILE *fp_snap, *fp_rcv;
+ char *file_snap, *file_rcv, file_tmp[150], ins[100], fbegin[100], fend[100], fins[100], fin2[100], numb1[100], *ptr;
+ float *rcvdata, *snapdata;
+ float dxs, dys, dzs, fxs, fys, fzs, sclr;
+ long nxs, nys, nzs, nts, ntrs, ret, file_det;
+ long it, ix, iy, iz, ixr, nxr, dnumb, numb, pos, nxmax;
+ long nzstart, nzend, dnz;
+ int *sx, *sy;
+ segy *hdr_snap, *hdr_rcv;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_rcv", &file_rcv)) file_rcv = "rcv.su";
+ if (!getparstring("file_snap", &file_snap)) file_snap = NULL;
+ if (file_snap == NULL) verr("No input data for the snap file given");
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=0;
+ if (!getparlong("nxmax", &nxmax)) nxmax=0;
+ if (!getparlong("nzstart", &nzstart)) nzstart=0;
+ if (!getparlong("nzend", &nzend)) nzend=0;
+
+ if (dnumb < 1) dnumb = 1;
+
+ sprintf(numb1,"%li",numb);
+
+ ptr = strstr(file_snap,numb1);
+ pos = ptr - file_snap + 1;
+
+ sprintf(fbegin,"%*.*s", pos-1, pos-1, file_snap);
+ sprintf(fend,"%s", file_snap+pos);
+
+ file_det = 1;
+ nxr=0;
+
+ while (file_det) {
+ sprintf(fins,"%li",nxr*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ fp_snap = fopen(fin2, "r");
+ if (fp_snap == NULL) {
+ if (nxr == 0) {
+ verr("error on opening basefile=%s", fin2);
+ }
+ else if (nxr == 1) {
+ vmess("1 file detected");
+ file_det = 0;
+ break;
+ }
+ else {
+ vmess("%d files detected",nxr);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_snap);
+ nxr++;
+ if (nxr!=0 && nxr == nxmax) {
+ vmess("%li files detected",nxr);
+ file_det = 0;
+ break;
+ }
+ }
+
+ sprintf(fins,"%li",numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+
+ getFileInfo3D(fin2, &nzs, &nxs, &nys, &nts, &dzs, &dxs, &dys, &fzs, &fxs, &fys, &sclr, &ntrs);
+
+ if (nzend > nzs || nzend==0) nzend=nzs;
+ if (nzstart < 0) nzstart=0;
+ if (nzstart > nzend) verr("The value of nzstart (%li) is greater than nzend (%li)",nzstart,nzend);
+
+ dnz = nzend-nzstart;
+
+ rcvdata = (float *)malloc(nxr*dnz*nxs*nys*nts*sizeof(float));
+ snapdata = (float *)malloc(dnz*nxs*nys*nts*sizeof(float));
+ hdr_snap = (segy *)calloc(nxs*nys*nts,sizeof(segy));
+ sx = (int *)malloc(nxr*sizeof(int));
+ sy = (int *)malloc(nxr*sizeof(int));
+
+ for (ixr=0; ixr<nxr; ixr++) {
+ vmess("Reshaping %li out of %li files",ixr+1,nxr);
+ sprintf(fins,"%li",ixr*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ readSnapData3D(fin2, snapdata, hdr_snap, nts, nxs, nys, nzs, 0, nxs, 0, nys, nzstart, nzend);
+ sx[ixr] = hdr_snap[0].sx;
+ sy[ixr] = hdr_snap[0].sy;
+ if (ixr==0) dzs = hdr_snap[0].d1;
+ for (iz=0; iz<dnz; iz++) {
+ for (iy=0; iy<nys; iy++) {
+ for (ix=0; ix<nxs; ix++) {
+ for (it=0; it<nts; it++) {
+ rcvdata[iz*nys*nxs*nxr*nts+iy*nxs*nxr*nts+ix*nxr*nts+ixr*nts+it] = snapdata[it*nys*nxs*dnz+iy*nxs*dnz+ix*dnz+iz];
+ }
+ }
+ }
+ }
+ }
+
+ free(snapdata);
+
+ hdr_rcv = (segy *)calloc(nxs*nys*nxr,sizeof(segy));
+
+ for (iz=nzstart; iz<nzend; iz++) {
+ vmess("Writing depth %li out of %li",iz-nzstart+1,dnz);
+ strcpy(file_tmp, file_rcv);
+ sprintf(ins,"_z%li", iz);
+ name_ext(file_tmp, ins);
+ fp_rcv = fopen(file_tmp, "w+");
+ if (fp_rcv==NULL) verr("error on creating output file %s", file_rcv);
+ // Set headers
+ for (iy=0; iy<nys; iy++){
+ for (ix=0; ix<nxs; ix++){
+ for (ixr=0; ixr<nxr; ixr++) {
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].ns = nts;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].fldr = iy*nxs+ix+1;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].tracf = ixr+1;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].tracl = iy*nxs*nxr+ix*nxr+ixr+1;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].trid = 1;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].scalco = -1000;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].scalel = -1000;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].sx = (int)((fxs+dxs*ix)*(1e3));
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].sy = (int)((fys+dys*iy)*(1e3));
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].gx = sx[ixr];
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].gy = sy[ixr];
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].sdepth = (int)((fzs+dzs*iz)*(1e3));
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].f1 = 0.0;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].f2 = sx[0]/(1e3);
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].d1 = (float)(hdr_snap[0].dt/1e6);
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].d2 = dxs;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].dt = hdr_snap[0].dt;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].trwf = nxs*nys*nxr;
+ hdr_rcv[iy*nxs*nxr+ix*nxr+ixr].ntr = nxs*nys*nxr;
+ }
+ }
+ }
+
+ // Write out homogeneous Green's function
+ ret = writeData3D(fp_rcv, (float *)&rcvdata[(iz-nzstart)*nys*nxs*nts*nxr], hdr_rcv, nts, nxr*nxs*nys);
+ if (ret < 0 ) verr("error on writing output file.");
+ fclose(fp_rcv);
+ }
+
+ free(rcvdata); free(hdr_rcv); free(hdr_snap);
+
+ return;
+}
\ No newline at end of file
diff --git a/utils/writeData3D.c b/utils/writeData3D.c
new file mode 100644
index 0000000..d8c3e47
--- /dev/null
+++ b/utils/writeData3D.c
@@ -0,0 +1,28 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "segy.h"
+
+/**
+* writes an 2D array to a SU file
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2)
+{
+ size_t nwrite;
+ long i;
+
+ for (i=0; i<n2; i++) {
+ nwrite = fwrite(&hdrs[i], 1, TRCBYTES, fp);
+ assert(nwrite == TRCBYTES);
+ nwrite = fwrite(&data[i*n1], sizeof(float), n1, fp);
+ assert (nwrite == n1);
+ }
+
+ return 0;
+}
+
--
GitLab