diff --git a/.gitignore b/.gitignore
index 72d36c2b470b91c15c8c9c8819a4ad75540a02e3..64ba6df2394e1590a987ec1ecae0ab82345fc1b3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,3 +23,7 @@ utils/syn2d
fdelmodc/fdelmodc
include/genfft.h
Make_include
+marchenko_applications/fmute
+marchenko_applications/marchenko
+marchenko_full/fmute
+marchenko_full/marchenko
diff --git a/Makefile b/Makefile
index 62c2fe254f7cebadc153b1a4ddf23d4fd9b5bbe4..cda1ad6cbee3775abbeeec2980fa403a9a7f2bae 100644
--- a/Makefile
+++ b/Makefile
@@ -5,6 +5,7 @@ all: mkdirs
cd fdelmodc ; $(MAKE) install
cd utils ; $(MAKE) install
cd marchenko ; $(MAKE) install
+ cd corrvir ; $(MAKE) install
mkdirs:
-mkdir -p lib
@@ -16,12 +17,14 @@ clean:
cd fdelmodc ; $(MAKE) $@
cd utils ; $(MAKE) $@
cd marchenko ; $(MAKE) $@
+ cd corrvir ; $(MAKE) $@
realclean:
cd FFTlib ; $(MAKE) $@
cd fdelmodc ; $(MAKE) $@
cd utils ; $(MAKE) $@
cd marchenko ; $(MAKE) $@
+ cd corrvir ; $(MAKE) $@
rm -f lib/*
rm -f include/*
rm -f bin/*
diff --git a/corrvir/Makefile b/corrvir/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..d39cd50d15a15954f5ddd50ff5d69ecbad4eecd3
--- /dev/null
+++ b/corrvir/Makefile
@@ -0,0 +1,47 @@
+# Makefile for corrvir
+
+include ../Make_include
+
+ALLINC = -I.
+LIBS += -L$L -lgenfft -lm
+#OPTC += -g
+
+all: corrvir
+
+PRG = corrvir
+
+SRCC = $(PRG).c \
+ getFileInfo.c \
+ correlate.c \
+ atopkge.c \
+ docpkge.c \
+ wallclock_time.c \
+ getpars.c
+
+OBJC = $(SRCC:%.c=%.o)
+
+$(PRG): $(OBJC) corr.h
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o $(PRG) $(OBJC) $(LIBS)
+
+install: $(PRG)
+ cp $(PRG) $B
+
+clean:
+ rm -f core $(OBJC) $(OBJM) $(PRG)
+
+realclean:
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/$(PRG)
+
+
+print: Makefile $(SRC)
+ $(PRINT) $?
+ @touch print
+
+count:
+ @wc $(SRC)
+
+tar:
+ @tar cf $(PRG).tar Makefile $(SRC) && compress $(PRG).tar
+
+
+
diff --git a/corrvir/atopkge.c b/corrvir/atopkge.c
new file mode 120000
index 0000000000000000000000000000000000000000..5107e2b2ccd382ede29d397838d8fad88126a516
--- /dev/null
+++ b/corrvir/atopkge.c
@@ -0,0 +1 @@
+../utils/atopkge.c
\ No newline at end of file
diff --git a/corrvir/corr.h b/corrvir/corr.h
new file mode 100644
index 0000000000000000000000000000000000000000..9030ee50d60e9ab5d3be94393172aee3cbee9a15
--- /dev/null
+++ b/corrvir/corr.h
@@ -0,0 +1,26 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+
+typedef struct _tracePos { /* Type */
+ int gx;
+ int gy;
+ int gelev;
+} tracePos;
+
+typedef struct _traceCoord { /* Type */
+ int x;
+ int y;
+ int peg;
+ int fpos;
+} traceCoord;
+
+typedef struct _crgPos { /* Type */
+ int gx;
+ int gy;
+ int rpeg;
+ int gelev;
+ int nsrc;
+ traceCoord *src;
+} crgPos;
+
diff --git a/corrvir/correlate.c b/corrvir/correlate.c
new file mode 100644
index 0000000000000000000000000000000000000000..1a4ae8de7879ad6520fb56e116fb23346f14e2c7
--- /dev/null
+++ b/corrvir/correlate.c
@@ -0,0 +1,97 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+
+double wallclock_time(void);
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+int correlate(complex *cmaster, complex *cslaves, int nfreq, int ncor, double *tcorr, int verbose)
+{
+ int j, istation, icc;
+ double t0, t1;
+ complex cC;
+
+ t0 = wallclock_time();
+ for (istation=0; istation<ncor; istation++) {
+ icc = istation*nfreq;
+
+#pragma ivdep
+ for (j=0; j<nfreq; j++) {
+ /* A*B */
+ cC.r = cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i;
+ cC.i = cmaster[j+icc].r*cslaves[icc+j].i-cmaster[j+icc].i*cslaves[icc+j].r;
+ /* AB* */
+// cC.r = cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i;
+// cC.i = -cmaster[j+icc].r*cslaves[icc+j].i+cmaster[j+icc].i*cslaves[icc+j].r;
+ /* A*B + AB* */
+// cC.r = cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i;
+// cC.i = 0.0;
+
+ cslaves[icc+j] = cC;
+ }
+
+ }
+ t1 = wallclock_time();
+ *tcorr += t1-t0;
+
+ return 0;
+}
+
+int coherence(complex *cmaster, complex *cslaves, int nfreq, int ncor, float reps, float epsmax, double *tcorr, int verbose)
+{
+ int j, istation, icc;
+ float maxden, *den, leps, am1, am2, scl;
+ double t0, t1;
+ complex cC;
+
+ t0 = wallclock_time();
+ den = (float *)malloc(nfreq*ncor*sizeof(float));
+ assert(den != NULL);
+
+ for (istation=0; istation<ncor; istation++) {
+ icc = istation*nfreq;
+ for (j=0; j<nfreq; j++) {
+ am1 = sqrt(cmaster[j+icc].r*cmaster[j+icc].r+cmaster[j+icc].i*cmaster[j+icc].i);
+ am2 = sqrt(cslaves[j+icc].r*cslaves[j+icc].r+cslaves[j+icc].i*cslaves[j+icc].i);
+ den[j+icc] = am1*am2;
+ }
+ }
+ leps = reps*maxden;
+
+ for (istation=0; istation<ncor; istation++) {
+ icc = istation*nfreq;
+
+ maxden=0.0;
+ for (j=0; j<nfreq; j++) {
+ maxden = MAX(maxden, den[j+icc]);
+ }
+ leps = reps*maxden;
+#pragma ivdep
+ for (j=0; j<nfreq; j++) {
+ /* A*B */
+ if (den[j+icc]>epsmax*maxden) scl = 1.0/(den[j+icc]);
+ else if (den[j+icc]<epsmax*maxden && den[j+icc]!=0) scl = 1.0/(den[j+icc]+leps);
+ else if (den[j+icc]==0) scl = 1.0;
+
+ cC.r = (cmaster[j+icc].r*cslaves[icc+j].r+cmaster[j+icc].i*cslaves[icc+j].i)*scl;
+ cC.i = (cmaster[j+icc].r*cslaves[icc+j].i-cmaster[j+icc].i*cslaves[icc+j].r)*scl;
+
+ cslaves[icc+j] = cC;
+ }
+
+ }
+ free(den);
+
+ t1 = wallclock_time();
+ *tcorr += t1-t0;
+
+ return 0;
+}
diff --git a/corrvir/corrvir.c b/corrvir/corrvir.c
new file mode 100644
index 0000000000000000000000000000000000000000..44a0025ce907aeec456c9d1794c58e6fa8689cc4
--- /dev/null
+++ b/corrvir/corrvir.c
@@ -0,0 +1,704 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+#include "par.h"
+#include "corr.h"
+#include "segy.h"
+
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+double wallclock_time(void);
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+void read_sutrace_at_position(FILE *fp, int itrace, complex *tracedata, complex *tracebuffer, int *trace_in_buffer, size_t *nbuffer, int ntfft, int nt, int nfreq, size_t nbufmax, size_t trace_sz, double *tread, double *tfft);
+
+int optncr(int n);
+void cr1fft(complex *cdata, float *data, int n, int sign);
+void rc1fft(float *rdata, complex *cdata, int n, int sign);
+
+int getFileInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, int verbose);
+
+int computeFFT(float *datain, int ld1, int nStart, int ntfft, int nfreq, int nstation, complex *cA, complex *cB, double *tfft, int verbose);
+
+int correlate(complex *cmaster, complex *cslaves, int nfreq, int ncor, double *tcorr, int verbose);
+int coherence(complex *cmaster, complex *cslaves, int nfreq, int ncor, float reps, float epsmax, double *tcorr, int verbose);
+
+/**************
+* ntc output samples of correlation result
+* note that nt (the number of samples read by the IO routine)
+* should be 2*ntc and a number efficient for FFT's
+*/
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" corrvir - correlation to compute virtual shot records",
+" ",
+" corrvir file_shots= file_out= [optional parameters]",
+" ",
+" Required parameters: ",
+" ",
+" file_shots= ....... filename of shots used to calculate virtual shots",
+" file_out= ......... filename of the output file with the virtual shot records",
+" nsources= ......... maximum number of actual shots in a common-receiver gather",
+" nreceivers= ....... maximum number of receiver positions",
+" xsrcv= ............ defines virtual-source x-positions",
+" ysrcv= ............ defines virtual-source y-positions",
+" ",
+" Optional parameters: ",
+" ",
+" derr = 0 ............ maximum tolerated error distance (m) between defined virtual-source position and receiver positions",
+" fmax = 70 ........... maximum frequency to use in correlation",
+" nbm = 4 ............. amount of memory (in GB) to buffer traces for reusage",
+" nsrc_ctr = 1 ...... . minumum number of sources to contribute to Fresnel stack",
+" normsrc = 0 ....... . normalize each correlated trace (before source summation)",
+" normalize = 1 ....... normalize the virtual trace with the number of contributing sources",
+" cohr = 0 ............ use cross-coherence: {f1(w)*f2(w)/(|f1(w)|*|f2(w)|+eps)} ",
+" reps = 0.0 .......... relative stabilization factor for cohr; eps=reps*|f1(w)|*|f2(w)|",
+" epsmax = 0.1 ........ cut off range above which spectrum is flattened",
+" src_sel = 0 ......... use all sources that are contributing ",
+" = 1 ......... use only sources that are outside square area enclosing master and slave positions",
+" = 2 ......... use only sources that are in donut",
+" bmin = 1 ............ define inner circle of donut",
+" bmax = 4 ............ define outer circle of donut",
+" ",
+" OUTPUT DEFINITION ",
+" causal=1 .......... output causal(1), non-causal(2), both(3), or summed(4)",
+" verbose=0 ......... silent option; >0 displays info",
+" ",
+" Notes: ",
+" ntc is the number of output samples of correlation result",
+" nt is the number of samples per trace read by the IO routine, should be 2*ntc and a number efficient for FFT's",
+" For causal=3: t=0 is sample nt/2 and trace goes from [-nt/2*dt - 0 - nt/2*dt]",
+" ",
+" Authors : Jan Thorbecke : 2013 (j.w.thorbecke@tudelft.nl)",
+" Boris Boullenger : 2015 (b.boullenger@tudelft.nl)",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+ int i, j, k, ret, nshots, maxnfiles, nfiles;
+ int size, n1, n2, ntfft, nf, ifreq, nfreq;
+ int verbose, causal;
+ int nt, nstation, ntc, nb, istation, jstation;
+ int pgsz, istep, nsrc_ctr, nsrc, nsrcMaster, nsrcSlave;
+ size_t outputSize, nwrite, cdatainSize, datainSize, cdataoutSize;
+ size_t jstep, lsz, offset, trace_sz, nread;
+ int ierr, itrace_out, done, cohr, normalize, normsrc;
+ int storedrcv, storedsrc, src_sel, read_trace;
+ int nsources, nreceivers, jsource, ivirtual, ivrcv, ivsrc, isrcs, isrcm, ivs;
+ int sxv, syv, sx, sy, svelev, sxm, sym, sxs, sys;
+ int cx3, cy3, gxv, gyv, gvelev, gx, gy, gelev;
+ double cx1, cy1, cx2, cy2;
+ int distmc2, distsc2, distsc2min, distsc2max;
+ int bmin, bmax;
+ int rpeg, speg, vrpeg, vspeg;
+ int xmin, xmax, ymin, ymax;
+ int *trace_in_buffer;
+ size_t nbuffer, nbufmax, nbufmemory;
+ size_t ntrace;
+ int nbm;
+ crgPos *rcvSrc;
+ int ircv, *isrc, *nsrc_crg, *vsrc;
+ float dx, fmax, df, d1, d2, scl, fscalco, sclfft;
+ float dt, reps, epsmax;
+ float derr;
+ float *r, *vtrace;
+ float diffx, diffy, dd, ddmin, ddtol;
+ float *xsrcv, *ysrcv;
+ float rms2, sclrms;
+ complex tmpcc;
+ int nxsrcv, nysrcv, nvsrc;
+ int scalco;
+ complex *c, *cmaster, *cslaves, *tracebuffer;
+ double t0, t1, t2, t3, tinit, tread, tfft, tcorr, twrite, tlogic, trdloc, tftloc, tcorrloc;
+ double tread_ivs, tfft_ivs, tcorr_ivs, twrite_ivs, tlogic_ivs;
+ char *file_shots, *file_out, filename[1024], basename[1200], base[1200], *pbase;
+ int pe=0, root_pe=0, npes=1, ipe, size_s;
+ FILE *fpin, *fpout;
+ segy hdr;
+ segy *segy_hdrs;
+ segy outputhdr;
+
+ t0 = wallclock_time();
+ tinit = twrite = tread = tcorr = tfft = tlogic = 0.0;
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparint("verbose", &verbose)) verbose = 0;
+ if (!getparstring("file_shots", &file_shots)) file_shots=NULL;
+ assert(file_shots != NULL);
+ if (!getparstring("file_out", &file_out)) file_out=NULL;
+ if (!getparint("nsources", &nsources)) nsources = 0;
+ assert(nsources!=0);
+ if (!getparint("nreceivers", &nreceivers)) nreceivers = 0;
+ assert(nreceivers!=0);
+ if (!getparint("nbm", &nbm)) nbm = 4;
+ nbufmemory=nbm;
+ if (!getparint("nsrc_ctr", &nsrc_ctr)) nsrc_ctr = 1;
+ if (!getparint("normsrc", &normsrc)) normsrc = 0;
+ if (!getparint("causal", &causal)) causal = 1;
+ if (!getparint("normalize", &normalize)) normalize = 1;
+ if (!getparint("cohr", &cohr)) cohr = 0;
+ if (!getparint("src_sel", &src_sel)) src_sel = 0;
+ if (!getparint("bmin", &bmin)) bmin = 1;
+ if (!getparint("bmax", &bmax)) bmax = 4;
+ if (!getparfloat("reps", &reps)) reps = 0.0;
+ if (!getparfloat("epsmax", &epsmax)) epsmax = 0.1;
+ if (!getparfloat("fmax", &fmax)) fmax = 70;
+ if (!getparfloat("derr", &derr)) derr = 0;
+ nxsrcv = countparval("xsrcv");
+ nysrcv = countparval("ysrcv");
+ if (nxsrcv != nysrcv) {
+ verr("Number of sources in array xsrcv (%d), ysrcv(%d) are not equal",nxsrcv, nysrcv);
+ }
+ xsrcv = (float *)malloc(nxsrcv*sizeof(float));
+ ysrcv = (float *)malloc(nxsrcv*sizeof(float));
+
+ getparfloat("xsrcv", xsrcv);
+ getparfloat("ysrcv", ysrcv);
+
+ getFileInfo(file_shots, &n1, &n2, &d1, &d2, verbose);
+ ntrace = n2;
+ nt = n1;
+ dt = d1;
+ fpin = fopen(file_shots, "r");
+ assert( fpin != NULL);
+ nread = fread( &hdr, 1, TRCBYTES, fpin );
+ assert(nread == TRCBYTES);
+ fpout = fopen(file_out, "w+");
+
+ fprintf(stderr,"nt=%d dt=%f ntrace=%ld\n",nt, dt, (long) ntrace);
+ t1 = wallclock_time();
+ tinit += t1-t0;
+
+ /*================ Read geometry of all traces in file_shots ================*/
+
+ trace_sz = sizeof(float)*(n1)+TRCBYTES;
+ segy_hdrs = (segy *)calloc(ntrace,sizeof(segy));
+ assert(segy_hdrs != NULL);
+ for (i=0; i<ntrace; i++) {
+ offset = i*trace_sz;
+ ret = fseeko(fpin , offset, SEEK_SET);
+ if (ret<0) perror("fseeko");
+ nread = fread( &segy_hdrs[i], 1, TRCBYTES, fpin );
+ assert(nread == TRCBYTES);
+ }
+ t2 = wallclock_time();
+ tread += t2-t1;
+
+
+ /*================ Sort traces in common receiver gathers ================*/
+
+ /* nreceivers is total number of receiver positions which is (much) smaller than the number of traces in the file */
+ /* Every receiver position has a number of sources: common receiver gather */
+
+ rcvSrc = (crgPos *)calloc(nreceivers,sizeof(crgPos));
+ for (i=0; i<nreceivers; i++) {
+ rcvSrc[i].src = (traceCoord *)calloc(nsources,sizeof(traceCoord));
+ }
+
+ rcvSrc[0].gx = segy_hdrs[0].gx;
+ rcvSrc[0].gy = segy_hdrs[0].gy;
+ rcvSrc[0].rpeg = segy_hdrs[0].gdel;
+ rcvSrc[0].gelev = segy_hdrs[0].gelev;
+ rcvSrc[0].src[0].x = segy_hdrs[0].sx;
+ rcvSrc[0].src[0].y = segy_hdrs[0].sy;
+ rcvSrc[0].src[0].peg = segy_hdrs[0].sdel;
+ scalco = segy_hdrs[0].scalco;
+ rcvSrc[0].src[0].fpos = 0;
+ rcvSrc[0].nsrc = 1;
+ ircv=1;
+ if (scalco < 0) fscalco = 1.0/(-1.0*scalco);
+ else fscalco = (float)scalco;
+
+ for (jstation=0; jstation<ntrace; jstation++) {
+ gx = segy_hdrs[jstation].gx;
+ gy = segy_hdrs[jstation].gy;
+ rpeg = segy_hdrs[jstation].gdel;
+ gelev = segy_hdrs[jstation].gelev;
+ sx = segy_hdrs[jstation].sx;
+ sy = segy_hdrs[jstation].sy;
+ speg = segy_hdrs[jstation].sdel;
+
+ /* bookkeeping: find out how many sources are contributing to each receiver position */
+ storedrcv = 0;
+ for (i=0; i<ircv; i++) {
+ if ( (rcvSrc[i].gx == gx) && (rcvSrc[i].gy == gy) ) {
+ storedrcv = 1;
+ storedsrc = 0;
+ for (j=0; j<rcvSrc[i].nsrc; j++) {
+ if ( (rcvSrc[i].src[j].x == sx) && (rcvSrc[i].src[j].y == sy) ) {
+ storedsrc=1;
+ break;
+ }
+ }
+ if ( !storedsrc ) {
+ rcvSrc[i].src[rcvSrc[i].nsrc].x = sx;
+ rcvSrc[i].src[rcvSrc[i].nsrc].y = sy;
+ rcvSrc[i].src[rcvSrc[i].nsrc].peg = speg;
+ rcvSrc[i].src[rcvSrc[i].nsrc].fpos = jstation;
+ rcvSrc[i].nsrc += 1;
+ }
+ break;
+ }
+ }
+ if ( !storedrcv ) {
+ if (verbose>=2) fprintf(stderr,"for jstation = %d number of receiver positions %d\n",jstation, ircv+1);
+ rcvSrc[ircv].gx = gx;
+ rcvSrc[ircv].gy = gy;
+ rcvSrc[ircv].rpeg = rpeg;
+ rcvSrc[ircv].gelev = gelev;
+ rcvSrc[ircv].src[0].x = sx;
+ rcvSrc[ircv].src[0].y = sy;
+ rcvSrc[ircv].src[0].peg = speg;
+ rcvSrc[ircv].src[0].fpos = jstation;
+ rcvSrc[ircv].nsrc = 1;
+ ircv++;
+ }
+ }
+ fprintf(stderr,"number of receiver positions in file = %d \n",ircv);
+ nreceivers=ircv;
+ nsources = 0;
+ nsrc_crg = (int *)calloc(nreceivers,sizeof(int));
+ for (i=0; i<nreceivers; i++) {
+ nsrc_crg[i] = rcvSrc[i].nsrc;
+ nsources = MAX(nsources,rcvSrc[i].nsrc);
+ if(verbose>=2) fprintf(stderr,"number of sources in receiver position %d = %d \n",i+1, rcvSrc[i].nsrc);
+ }
+
+
+
+ /*================ Find receiver positions corresponding to requested virtual-sources ================*/
+ if (nxsrcv) {
+ nvsrc = nxsrcv;
+ vsrc = (int *)malloc(nvsrc*sizeof(int));
+ ddtol= derr*derr;
+ for (j=0; j<nxsrcv; j++) {
+ diffx = xsrcv[j]- rcvSrc[0].gx*fscalco;
+ diffy = ysrcv[j]- rcvSrc[0].gy*fscalco;
+ ddmin = diffx*diffx+diffy*diffy;
+ vsrc[j] = 0;
+ for (i=1; i<nreceivers; i++) {
+ diffx = xsrcv[j]- rcvSrc[i].gx*fscalco;
+ diffy = ysrcv[j]- rcvSrc[i].gy*fscalco;
+ dd = diffx*diffx+diffy*diffy;
+ if (dd < ddmin ) {
+ ddmin = dd;
+ vsrc[j] = i;
+ }
+ }
+ fprintf(stderr,"ddmin for requested vsrc position (%.2f,%.2f) = %.2f\n",xsrcv[j],ysrcv[j],ddmin);
+ if (ddmin <= ddtol) {
+ xsrcv[j]=rcvSrc[vsrc[j]].gx*fscalco;
+ ysrcv[j]=rcvSrc[vsrc[j]].gy*fscalco;
+ fprintf(stderr,"Found vsrc position: (%.2f,%.2f)\n",xsrcv[j],ysrcv[j]);
+ }
+ else {
+ fprintf(stderr,"Error: Requested vsrc position (%.2f,%.2f) do not coincide with a rcv position\n",xsrcv[j],ysrcv[j]);
+ }
+ }
+ }
+ else { /* all receivers are made virtual sources */
+ nvsrc = nreceivers;
+ vsrc = (int *)malloc(nvsrc*sizeof(int));
+ for (i=0; i<nreceivers; i++) {
+ vsrc[i] = i;
+ }
+ }
+
+
+ /*================ initializations ================*/
+
+ ntfft = optncr(2*nt);
+ nf = ntfft/2+1;
+ df = 1.0/(ntfft*dt);
+ nfreq = MIN(nf,(int)(fmax/df)+1);
+ sclfft= 1.0/((float)ntfft);
+
+// for (j=0; j<MIN(rcvSrc[0].nsrc,rcvSrc[1].nsrc); j++) fprintf(stderr,"before allocs sx=%d sy=%d +1 sx=%d sx=%d\n",rcvSrc[0].src[j].x, rcvSrc[0].src[j].y, rcvSrc[1].src[j].x, rcvSrc[1].src[j].y);
+
+ nbuffer = 0;
+ nbufmax = (nbufmemory*1024*1024*1024/(nfreq*sizeof(complex))); /* number of traces in buffer */
+ tracebuffer = (complex *)malloc(nbufmax*nfreq*sizeof(complex));
+ trace_in_buffer = (int *)malloc(nbufmax*sizeof(int));
+ for (i=0; i<nbufmax; i++) trace_in_buffer[i] = -1;
+
+
+ /*================ for nvirtual shots find suitable receivers ================*/
+
+ cmaster = (complex *)malloc(nsources*nfreq*sizeof(complex));
+ cslaves = (complex *)malloc(nsources*nfreq*sizeof(complex));
+
+ vtrace = (float *)malloc(nt*sizeof(float));
+ r = (float *)malloc(ntfft*sizeof(float));
+ c = (complex *)malloc(nf*sizeof(complex));
+
+ t1 = wallclock_time();
+ tinit += t1-t2;
+
+ itrace_out=0;
+ for (ivs=0; ivs<nvsrc; ivs++) {/* loop over the number of virtual source positions to be created */
+
+ t1 = wallclock_time();
+ tread_ivs = tfft_ivs = tcorr_ivs = twrite_ivs = tlogic_ivs = 0.0;
+
+ ivsrc=vsrc[ivs];
+ sxv = rcvSrc[ivsrc].gx;
+ syv = rcvSrc[ivsrc].gy;
+ vspeg = rcvSrc[ivsrc].rpeg;
+ svelev = rcvSrc[ivsrc].gelev;
+ nsrcMaster = nsrc_crg[ivsrc];
+
+ for (ivrcv=0; ivrcv<nreceivers; ivrcv++) {
+ gxv = rcvSrc[ivrcv].gx;
+ gyv = rcvSrc[ivrcv].gy;
+ vrpeg = rcvSrc[ivrcv].rpeg;
+ gvelev = rcvSrc[ivrcv].gelev;
+ nsrcSlave = nsrc_crg[ivrcv];
+
+ if (src_sel == 1) {
+ xmin = MIN(gxv,sxv);
+ ymin = MIN(gyv,syv);
+ xmax = MAX(gxv,sxv);
+ ymax = MAX(gyv,syv);
+ }
+
+ if (src_sel == 2) {
+ cx1 = 0.5*(sxv+gxv);
+ cy1 = 0.5*(syv+gyv);
+ cx2 = ceil(cx1);
+ cy2 = ceil(cy1);
+ cx3 = cx2;
+ cy3 = cy2;
+ distmc2 = (sxv-cx3)*(sxv-cx3)+(syv-cy3)*(syv-cy3);
+ distsc2min = distmc2+bmin*bmin*distmc2;
+ distsc2max = distmc2+bmax*bmax*distmc2;
+ distsc2 = (sxv-cx3)*(sxv-cx3)+(syv-cy3)*(syv-cy3);
+ }
+
+ //fprintf(stderr,"ivsrc=%d sxv=%d nsrcMaster=%d ivrcv=%d gxv=%d nsrcSlave=%d\n", ivsrc, sxv, nsrcMaster, ivrcv, gxv, nsrcSlave);
+
+ nsrc = 0;
+ memset(cslaves,0,nfreq*nsources*sizeof(complex));
+ trdloc = tftloc = tcorrloc = 0.0;
+ for (isrcm=0; isrcm<nsrcMaster; isrcm++) { /* for all traces find which traces are recorded by both virtual source and receiver position */
+ sxm = rcvSrc[ivsrc].src[isrcm].x;
+ sym = rcvSrc[ivsrc].src[isrcm].y;
+ for (isrcs=0; isrcs<nsrcSlave; isrcs++) {
+ sxs = rcvSrc[ivrcv].src[isrcs].x;
+ sys = rcvSrc[ivrcv].src[isrcs].y;
+ if ( (sxm == sxs) && (sym == sys) ) { /* master and slave have same contributing source */
+
+ /* if the source is outside the inclosing area of the virtual-receiver and virtual-source coordinates
+ then the source is accepted for contribution to the final summation */
+ if (src_sel == 0) {
+ read_trace=1;
+ }
+ else if (src_sel == 1) {
+ read_trace = 0;
+ if ( !( (sxs>xmin) && (sxs<xmax) ) ) {
+ if ( !( (sys>ymin) && (sys<ymax) ) ) {
+ read_trace=1;
+ }
+ }
+ }
+ else if (src_sel == 2) {
+ read_trace = 0;
+ if ( !( (distsc2 < distsc2min) && (distsc2 > distsc2max) ) ) {
+ read_trace = 1;
+ }
+ }
+ else if (src_sel == 3) {
+ read_trace = 0;
+ if ( (sxv < gxv) && (syv < gyv) ) {
+ if ( (sxm < sxv) && (sym < syv) ) {
+ read_trace = 1;
+ }
+ }
+ else if ( (sxv < gxv) && (syv > gyv) ) {
+ if ( (sxm < sxv) && (sym > syv) ) {
+ read_trace = 1;
+ }
+ }
+ else if ( (sxv > gxv) && (syv < gyv) ) {
+ if ( (sxm > sxv) && (sym < syv) ) {
+ read_trace = 1;
+ }
+ }
+ else if ( (sxv > gxv) && (syv > gyv) ) {
+ if ( (sxm > sxv) && (sym > syv) ) {
+ read_trace = 1;
+ }
+ }
+
+ }
+ else if (src_sel == 4) {
+ read_trace = 0;
+ if ( (sxv < gxv) && (syv < gyv) ) {
+ if ( (sxm < gxv) && (sym < gyv) ) {
+ read_trace = 1;
+ }
+ }
+ else if ( (sxv < gxv) && (syv > gyv) ) {
+ if ( (sxm < gxv) && (sym > gyv) ) {
+ read_trace = 1;
+ }
+ }
+ else if ( (sxv > gxv) && (syv < gyv) ) {
+ if ( (sxm > gxv) && (sym < gyv) ) {
+ read_trace = 1;
+ }
+ }
+ else if ( (sxv > gxv) && (syv > gyv) ) {
+ if ( (sxm > gxv) && (sym > gyv) ) {
+ read_trace = 1;
+ }
+ }
+
+ }
+ /* read data from file for sources that passed the tests above */
+ if (read_trace) {
+ /* read master trace */
+ read_sutrace_at_position(fpin, rcvSrc[ivsrc].src[isrcm].fpos, &cmaster[nsrc*nfreq], tracebuffer, trace_in_buffer, &nbuffer, ntfft, nt, nfreq, nbufmax, trace_sz, &trdloc, &tftloc);
+
+ /* read slave trace */
+ read_sutrace_at_position(fpin, rcvSrc[ivrcv].src[isrcs].fpos, &cslaves[nsrc*nfreq], tracebuffer, trace_in_buffer, &nbuffer, ntfft, nt, nfreq, nbufmax, trace_sz, &trdloc, &tftloc);
+ nsrc++;
+ }
+
+ } /* end of if test that slave and master have same contributing source */
+ } /* end of nsrcSlave loop */
+ } /* end of nsrcMaster loop */
+ tread_ivs += trdloc;
+ tread += trdloc;
+ tfft_ivs += tftloc;
+ tfft += tftloc;
+
+ if (nsrc < nsrc_ctr) nsrc = -1; /* only compute virtual receiver when there are sufficient active sources contributing */
+
+ /* correlation of virtual source ivsrc with virtual receiver ivcrv */
+ if (cohr) {
+ coherence(cmaster, cslaves, nfreq, nsrc, reps, epsmax, &tcorrloc, verbose);
+ }
+ else {
+ correlate(cmaster, cslaves, nfreq, nsrc, &tcorrloc, verbose);
+ }
+ tcorr_ivs += tcorrloc;
+ tcorr += tcorrloc;
+
+ /* summation over the contributing sources */
+ memset(c,0,nf*sizeof(complex));
+ for (jsource=0; jsource<nsrc; jsource++) {
+ if (normsrc == 1) {
+ rms2 = 0.0;
+ for (ifreq=0; ifreq<nfreq; ifreq++) {
+ tmpcc = cslaves[jsource*nfreq+ifreq];
+ rms2 += tmpcc.r*tmpcc.r+tmpcc.i*tmpcc.i;
+ }
+ if (rms2 > 0) sclrms = 1.0/sqrt(rms2);
+ else sclrms = 1.0;
+ }
+ else (sclrms = 1.0);
+ for (ifreq=0; ifreq<nfreq; ifreq++) {
+ c[ifreq].r += cslaves[jsource*nfreq+ifreq].r*sclrms;
+ c[ifreq].i += cslaves[jsource*nfreq+ifreq].i*sclrms;
+ }
+ }
+
+ t2 = wallclock_time();
+
+ cr1fft(&c[0],r,ntfft,1); /* transform virtual trace back to time */
+
+ t3 = wallclock_time();
+ tfft_ivs += t3-t2;
+ tfft += t3-t2;
+
+ if (normalize && nsrc>0) scl=sclfft/nsrc; /* normalize by the number of contributing sources */
+ else scl = sclfft;
+
+ if (causal==1) {
+ for (j=0;j<nt; j++) {
+ vtrace[j] = r[j]*scl;
+ }
+ }
+ else if (causal==2) {
+ vtrace[0] = r[0]*scl;
+ for (j=1;j<nt; j++) {
+ vtrace[j] = r[ntfft-j]*scl;
+ }
+ }
+ else if (causal==3) {
+ for (j=1;j<=(nt/2); j++) {
+ vtrace[nt/2-j] = r[ntfft-j]*scl;
+ }
+ for (j=nt/2;j<nt; j++) {
+ vtrace[j] = r[j-nt/2]*scl;
+ }
+ }
+ else if (causal==4) {
+ vtrace[0] = r[0]*scl;
+ for (j=1;j<nt; j++) {
+ vtrace[j] = 0.5*(r[ntfft-j] + r[j])*scl;
+ }
+ } /* one virtual source and virtual receiver are now calculated; write the common virtual shot gather to disk */
+
+ t2 = wallclock_time();
+
+ memcpy(&outputhdr,&segy_hdrs[(rcvSrc[ivrcv].src[0].fpos)],TRCBYTES);
+ outputhdr.fldr = ivsrc+1;
+ outputhdr.tracl = itrace_out+1;
+ outputhdr.tracf = ivrcv+1;
+ outputhdr.sx = sxv;
+ outputhdr.sy = syv;
+ outputhdr.selev = svelev;
+ outputhdr.gx = gxv;
+ outputhdr.gy = gyv;
+ outputhdr.sdel = vspeg;
+ outputhdr.gdel = vrpeg;
+ outputhdr.gelev = gvelev;
+ outputhdr.ns = nt;
+ // outputhdr.offset = (gxv - sxv)/1000;
+ nwrite = fwrite( &outputhdr, 1, TRCBYTES, fpout );
+ assert (nwrite == TRCBYTES);
+ nwrite = fwrite( &vtrace[0], sizeof(float), nt, fpout );
+ assert (nwrite == nt);
+ fflush(fpout);
+
+ t3 = wallclock_time();
+ twrite_ivs += t3-t2;
+ twrite += t3-t2;
+
+ itrace_out++;
+ } /* end of virtual receiver loop */
+
+ if (verbose>=3) {
+ t3 = wallclock_time();
+ tlogic_ivs = ((t3-t1)-tread_ivs-tfft_ivs-tcorr_ivs-twrite_ivs);
+ fprintf(stderr,"************* Timings ************* vshot= %.d \n", vspeg);
+ fprintf(stderr,"CPU-time read data = %.3f\n", tread_ivs);
+ fprintf(stderr,"CPU-time FFT's = %.3f\n", tfft_ivs);
+ fprintf(stderr,"CPU-time correlation = %.3f\n", tcorr_ivs);
+ fprintf(stderr,"CPU-time write data = %.3f\n", twrite_ivs);
+ fprintf(stderr,"CPU-time logic = %.3f\n", tlogic_ivs);
+ fprintf(stderr,"Total CPU-time this shot = %.3f\n", t3-t1);
+ }
+ } /* end of virtual source loop */
+
+ fclose(fpin);
+ fclose(fpout);
+
+ if (verbose) {
+ t3 = wallclock_time();
+ tlogic = ((t3-t0)-tread-tfft-tcorr-twrite);
+ fprintf(stderr,"************* Total Timings ************* \n");
+ fprintf(stderr,"CPU-time initilization = %.3f\n", tinit);
+ fprintf(stderr,"CPU-time read data = %.3f\n", tread);
+ fprintf(stderr,"CPU-time FFT's = %.3f\n", tfft);
+ fprintf(stderr,"CPU-time correlation = %.3f\n", tcorr);
+ fprintf(stderr,"CPU-time write data = %.3f\n", twrite);
+ fprintf(stderr,"CPU-time logic = %.3f\n", tlogic);
+ fprintf(stderr,"Total CPU-time = %.3f\n", t3-t0);
+ }
+
+ free(cmaster);
+ free(cslaves);
+ free(vtrace);
+ free(c);
+ free(r);
+
+/*================ end ================*/
+
+ return 0;
+}
+
+
+void read_sutrace_at_position(FILE *fp, int itrace, complex *tracedata, complex *tracebuffer, int *trace_in_buffer, size_t *nbuffer,
+int ntfft, int nt, int nfreq, size_t nbufmax, size_t trace_sz, double *tread, double *tfft)
+{
+ static size_t freebuffer=0;
+ size_t i, inbuffer, offset, nread, nwrite;
+ int ret;
+ float *trace;
+ float *r;
+ double t0, t1, t2, t3;
+ complex *c;
+ FILE *fpt;
+ char filename[1024];
+
+ inbuffer = -1;
+// fprintf(stderr,"reading trace %d\n", itrace);
+
+ for (i=0; i<nbufmax; i++) {
+ if (itrace == trace_in_buffer[i]) {
+ inbuffer = i;
+ break;
+ }
+ }
+
+// fprintf(stderr,"inbuffer = %d bufmax=%d\n", inbuffer,nbufmax);
+ if (inbuffer == -1) {
+ t0 = wallclock_time();
+ trace = (float *)calloc(nt, sizeof(float));
+ r = (float *)malloc(ntfft*sizeof(float));
+ c = (complex *)malloc((ntfft/2+1)*sizeof(complex));
+
+ /* read trace from file */
+ offset = itrace*trace_sz+TRCBYTES;
+ ret = fseeko(fp , offset, SEEK_SET);
+ if (ret<0) perror("fseeko");
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert(nread == nt);
+
+ /* transform to frequency domain */
+ memset(r,0,ntfft*sizeof(float));
+ memcpy(r,&trace[0],nt*sizeof(float)); /* might not be needed if segy_read_traces_c is nice */
+ t1 = wallclock_time();
+ *tread += t1-t0;
+
+ rc1fft(r,c,ntfft,-1);
+
+ t2 = wallclock_time();
+ *tfft += t2-t1;
+
+ /* reset buffer counter when buffer is full */
+ if (freebuffer<nbufmax-1) {
+ freebuffer++;
+ }
+ else freebuffer=0;
+
+// fprintf(stderr,"freebuffer = %d\n", freebuffer);
+// fprintf(stderr,"nbuffer = %d\n", *nbuffer);
+ memcpy(&tracebuffer[freebuffer*nfreq].r,&c[0].r,nfreq*sizeof(complex));
+ memcpy(&tracedata[0].r,&tracebuffer[freebuffer*nfreq].r,nfreq*sizeof(complex));
+ trace_in_buffer[freebuffer]=itrace;
+
+ free(trace);
+ free(r);
+ free(c);
+ t1 = wallclock_time();
+ *tread += t1-t2;
+
+ }
+ else {
+ /* read trace from buffer */
+ t0 = wallclock_time();
+// tracedata = (complex *)&tracebuffer[inbuffer*nfreq]; /* this is better, but need some more changes TODO */
+
+// fprintf(stderr,"read from inbuffer = %d bufmax=%d\n", inbuffer,nbufmax);
+
+ memcpy(&tracedata[0].r,&tracebuffer[inbuffer*nfreq].r,nfreq*sizeof(complex));
+ t1 = wallclock_time();
+ *tread += t1-t0;
+ }
+
+ return;
+}
+
diff --git a/corrvir/demo/compvir_all.sh b/corrvir/demo/compvir_all.sh
new file mode 100755
index 0000000000000000000000000000000000000000..08c66a3b08588739aa0f6aee86dd2a132e0604f0
--- /dev/null
+++ b/corrvir/demo/compvir_all.sh
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+filename_in=/vardim/home/boulleng/Canada_Data/Lalor_shots_clean_4000ms_agc_sortgdelsdel_sdel136out.su
+filename_out=virshots_fmax160_FROM_all.su
+
+
+ulimit -c unlimited
+../corrvir \
+ file_shots=${filename_in} \
+ file_out=${filename_out} \
+ nsources=908 \
+ src_sel=0 \
+ fmax=160 \
+ nbm=1 \
+ nreceivers=2685 \
+ normsrc=0 \
+ normalize=0 \
+ cohr=0 \
+ causal=1 \
+ verbose=4
diff --git a/corrvir/docpkge.c b/corrvir/docpkge.c
new file mode 120000
index 0000000000000000000000000000000000000000..5384bb3801703c3f0db8fcc032235ca6130fa08b
--- /dev/null
+++ b/corrvir/docpkge.c
@@ -0,0 +1 @@
+../utils/docpkge.c
\ No newline at end of file
diff --git a/corrvir/getFileInfo.c b/corrvir/getFileInfo.c
new file mode 100644
index 0000000000000000000000000000000000000000..152c5e290c4632d81b86f2d9a543ff4efddb3612
--- /dev/null
+++ b/corrvir/getFileInfo.c
@@ -0,0 +1,43 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include "par.h"
+#include "segy.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+int getFileInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, int verbose)
+{
+ FILE *fp;
+ size_t nread;
+ off_t bytes, ret, trace_sz, ntraces;
+ int i, itrace, one_shot, igath;
+ float *trace, cmin;
+ segy hdr;
+
+ fp = fopen( file_name, "r" );
+ assert( fp != NULL);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ ret = fseeko( fp, 0, SEEK_END );
+ if (ret<0) perror("fseeko");
+ bytes = ftello( fp );
+
+ *n1 = hdr.ns;
+ *d1 =1e-6*hdr.dt;
+ *d2 = hdr.d2;
+
+ trace_sz = sizeof(float)*(*n1)+TRCBYTES;
+ ntraces = (int) (bytes/trace_sz);
+ *n2 = ntraces;
+
+ if (verbose) {
+ fprintf(stderr,"file_base=%s has nt=%d dt=%f and %d traces\n", file_name, *n1, *d1, *n2);
+ }
+
+ return 0;
+}
+
diff --git a/corrvir/getpars.c b/corrvir/getpars.c
new file mode 120000
index 0000000000000000000000000000000000000000..fa7dc3355428e8ea9013fafad6e319dde3a48ebb
--- /dev/null
+++ b/corrvir/getpars.c
@@ -0,0 +1 @@
+../utils/getpars.c
\ No newline at end of file
diff --git a/corrvir/par.h b/corrvir/par.h
new file mode 120000
index 0000000000000000000000000000000000000000..0fa273cea748f9ead16e0e231201941174a3dd46
--- /dev/null
+++ b/corrvir/par.h
@@ -0,0 +1 @@
+../utils/par.h
\ No newline at end of file
diff --git a/corrvir/segy.h b/corrvir/segy.h
new file mode 100644
index 0000000000000000000000000000000000000000..326f285caf8bc14f17331d168904553858cab030
--- /dev/null
+++ b/corrvir/segy.h
@@ -0,0 +1,326 @@
+/* Copyright (c) Colorado School of Mines, 2001.*/
+/* All rights reserved. */
+
+/* segy.h - include file for SEGY traces
+ *
+ * declarations for:
+ * typedef struct {} segy - the trace identification header
+ * typedef struct {} bhed - binary header
+ *
+ * Note:
+ * If header words are added, run the makefile in this directory
+ * to recreate hdr.h.
+ *
+ * Reference:
+ * K. M. Barry, D. A. Cavers and C. W. Kneale, "Special Report:
+ * Recommended Standards for Digital Tape Formats",
+ * Geophysics, vol. 40, no. 2 (April 1975), P. 344-352.
+ *
+ * $Author: john $
+ * $Source: /usr/local/cwp/src/su/include/RCS/segy.h,v $
+ * $Revision: 1.23 $ ; $Date: 1998/03/26 23:48:18 $
+ */
+
+#define TRCBYTES 240
+
+
+/* TYPEDEFS */
+typedef struct { /* segy - trace identification header */
+
+ int tracl; /* trace sequence number within line */
+
+ int tracr; /* trace sequence number within reel */
+
+ int fldr; /* field record number */
+
+ int tracf; /* trace number within field record */
+
+ int ep; /* energy source point number */
+
+ int cdp; /* CDP ensemble number */
+
+ int cdpt; /* trace number within CDP ensemble */
+
+ short trid; /* trace identification code:
+ 1 = seismic data
+ 2 = dead
+ 3 = dummy
+ 4 = time break
+ 5 = uphole
+ 6 = sweep
+ 7 = timing
+ 8 = water break
+ 9---, N = optional use (N = 32,767)
+
+ Following are CWP id flags:
+
+ 9 = autocorrelation
+
+ 10 = Fourier transformed - no packing
+ xr[0],xi[0], ..., xr[N-1],xi[N-1]
+
+ 11 = Fourier transformed - unpacked Nyquist
+ xr[0],xi[0],...,xr[N/2],xi[N/2]
+
+ 12 = Fourier transformed - packed Nyquist
+ even N:
+ xr[0],xr[N/2],xr[1],xi[1], ...,
+ xr[N/2 -1],xi[N/2 -1]
+ (note the exceptional second entry)
+ odd N:
+ xr[0],xr[(N-1)/2],xr[1],xi[1], ...,
+ xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2]
+ (note the exceptional second & last entries)
+
+ 13 = Complex signal in the time domain
+ xr[0],xi[0], ..., xr[N-1],xi[N-1]
+
+ 14 = Fourier transformed - amplitude/phase
+ a[0],p[0], ..., a[N-1],p[N-1]
+
+ 15 = Complex time signal - amplitude/phase
+ a[0],p[0], ..., a[N-1],p[N-1]
+
+ 16 = Real part of complex trace from 0 to Nyquist
+
+ 17 = Imag part of complex trace from 0 to Nyquist
+
+ 18 = Amplitude of complex trace from 0 to Nyquist
+
+ 19 = Phase of complex trace from 0 to Nyquist
+
+ 21 = Wavenumber time domain (k-t)
+
+ 22 = Wavenumber frequency (k-omega)
+
+ 23 = Envelope of the complex time trace
+
+ 24 = Phase of the complex time trace
+
+ 25 = Frequency of the complex time trace
+
+ 30 = Depth-Range (z-x) traces
+
+ 43 = Seismic Data, Vertical Component
+
+ 44 = Seismic Data, Horizontal Component 1
+
+ 45 = Seismic Data, Horizontal Component 2
+
+ 46 = Seismic Data, Radial Component
+
+ 47 = Seismic Data, Transverse Component
+
+ 101 = Seismic data packed to bytes (by supack1)
+
+ 102 = Seismic data packed to 2 bytes (by supack2)
+ */
+
+ short nvs; /* number of vertically summed traces (see vscode
+ in bhed structure) */
+
+ short nhs; /* number of horizontally summed traces (see vscode
+ in bhed structure) */
+
+ short duse; /* data use:
+ 1 = production
+ 2 = test */
+
+ int offset; /* distance from source point to receiver
+ group (negative if opposite to direction
+ in which the line was shot) */
+
+ int gelev; /* receiver group elevation from sea level
+ (above sea level is positive) */
+
+ int selev; /* source elevation from sea level
+ (above sea level is positive) */
+
+ int sdepth; /* source depth (positive) */
+
+ int gdel; /* datum elevation at receiver group */
+
+ int sdel; /* datum elevation at source */
+
+ int swdep; /* water depth at source */
+
+ int gwdep; /* water depth at receiver group */
+
+ short scalel; /* scale factor for previous 7 entries
+ with value plus or minus 10 to the
+ power 0, 1, 2, 3, or 4 (if positive,
+ multiply, if negative divide) */
+
+ short scalco; /* scale factor for next 4 entries
+ with value plus or minus 10 to the
+ power 0, 1, 2, 3, or 4 (if positive,
+ multiply, if negative divide) */
+
+ int sx; /* X source coordinate */
+
+ int sy; /* Y source coordinate */
+
+ int gx; /* X group coordinate */
+
+ int gy; /* Y group coordinate */
+
+ short counit; /* coordinate units code:
+ for previous four entries
+ 1 = length (meters or feet)
+ 2 = seconds of arc (in this case, the
+ X values are longitude and the Y values
+ are latitude, a positive value designates
+ the number of seconds east of Greenwich
+ or north of the equator */
+
+ short wevel; /* weathering velocity */
+
+ short swevel; /* subweathering velocity */
+
+ short sut; /* uphole time at source */
+
+ short gut; /* uphole time at receiver group */
+
+ short sstat; /* source static correction */
+
+ short gstat; /* group static correction */
+
+ short tstat; /* total static applied */
+
+ short laga; /* lag time A, time in ms between end of 240-
+ byte trace identification header and time
+ break, positive if time break occurs after
+ end of header, time break is defined as
+ the initiation pulse which maybe recorded
+ on an auxiliary trace or as otherwise
+ specified by the recording system */
+
+ short lagb; /* lag time B, time in ms between the time break
+ and the initiation time of the energy source,
+ may be positive or negative */
+
+ short delrt; /* delay recording time, time in ms between
+ initiation time of energy source and time
+ when recording of data samples begins
+ (for deep water work if recording does not
+ start at zero time) */
+
+ short muts; /* mute time--start */
+
+ short mute; /* mute time--end */
+
+ unsigned short ns; /* number of samples in this trace */
+
+ unsigned short dt; /* sample interval; in micro-seconds */
+
+ short gain; /* gain type of field instruments code:
+ 1 = fixed
+ 2 = binary
+ 3 = floating point
+ 4 ---- N = optional use */
+
+ short igc; /* instrument gain constant */
+
+ short igi; /* instrument early or initial gain */
+
+ short corr; /* correlated:
+ 1 = no
+ 2 = yes */
+
+ short sfs; /* sweep frequency at start */
+
+ short sfe; /* sweep frequency at end */
+
+ short slen; /* sweep length in ms */
+
+ short styp; /* sweep type code:
+ 1 = linear
+ 2 = cos-squared
+ 3 = other */
+
+ short stas; /* sweep trace length at start in ms */
+
+ short stae; /* sweep trace length at end in ms */
+
+ short tatyp; /* taper type: 1=linear, 2=cos^2, 3=other */
+
+ short afilf; /* alias filter frequency if used */
+
+ short afils; /* alias filter slope */
+
+ short nofilf; /* notch filter frequency if used */
+
+ short nofils; /* notch filter slope */
+
+ short lcf; /* low cut frequency if used */
+
+ short hcf; /* high cut frequncy if used */
+
+ short lcs; /* low cut slope */
+
+ short hcs; /* high cut slope */
+
+ short year; /* year data recorded */
+
+ short day; /* day of year */
+
+ short hour; /* hour of day (24 hour clock) */
+
+ short minute; /* minute of hour */
+
+ short sec; /* second of minute */
+
+ short timbas; /* time basis code:
+ 1 = local
+ 2 = GMT
+ 3 = other */
+
+ short trwf; /* trace weighting factor, defined as 1/2^N
+ volts for the least sigificant bit */
+
+ short grnors; /* geophone group number of roll switch
+ position one */
+
+ short grnofr; /* geophone group number of trace one within
+ original field record */
+
+ short grnlof; /* geophone group number of last trace within
+ original field record */
+
+ short gaps; /* gap size (total number of groups dropped) */
+
+ short otrav; /* overtravel taper code:
+ 1 = down (or behind)
+ 2 = up (or ahead) */
+
+ /* local assignments */
+ float d1; /* sample spacing for non-seismic data */
+
+ float f1; /* first sample location for non-seismic data */
+
+ float d2; /* sample spacing between traces */
+
+ float f2; /* first trace location */
+
+ float ungpow; /* negative of power used for dynamic
+ range compression */
+
+ float unscale; /* reciprocal of scaling factor to normalize
+ range */
+
+ int ntr; /* number of traces */
+
+ short mark; /* mark selected traces */
+
+ short shortpad; /* alignment padding */
+
+
+ short unass[14]; /* unassigned--NOTE: last entry causes
+ a break in the word alignment, if we REALLY
+ want to maintain 240 bytes, the following
+ entry should be an odd number of short/UINT2
+ OR do the insertion above the "mark" keyword
+ entry */
+
+} segy;
+
diff --git a/corrvir/wallclock_time.c b/corrvir/wallclock_time.c
new file mode 120000
index 0000000000000000000000000000000000000000..0bd00b4c2878f007a8dc398f0af7c7cb44f50717
--- /dev/null
+++ b/corrvir/wallclock_time.c
@@ -0,0 +1 @@
+../utils/wallclock_time.c
\ No newline at end of file
diff --git a/doc/fdelmodcManual.pdf b/doc/fdelmodcManual.pdf
index a6ec5a39a708ac0c4636ba578292c176ec484a3e..3cebcbe08381ab0575edb73d8192bfa034eea3e5 100644
Binary files a/doc/fdelmodcManual.pdf and b/doc/fdelmodcManual.pdf differ
diff --git a/fdelmodc/decomposition.c b/fdelmodc/decomposition.c
index 6558de07b5a75d82e89429a7c5853703993a96dd..3b487787b5323c633c5dfcb15dfbcc83f888f53d 100644
--- a/fdelmodc/decomposition.c
+++ b/fdelmodc/decomposition.c
@@ -5,7 +5,6 @@
* GEOPHYSICS, VOL. 63, NO. 4 (JULY-AUGUST 1998); P. 1795–1798
*
* Created by Jan Thorbecke on 17/03/2014.
- * Copyright 2014 TU Delft. All rights reserved.
*
*/
diff --git a/fdelmodc/decomposition.c.new b/fdelmodc/decomposition.c.new
deleted file mode 100644
index 3b03ba5a6c5956578ee28c89977bdde47c9d073f..0000000000000000000000000000000000000000
--- a/fdelmodc/decomposition.c.new
+++ /dev/null
@@ -1,247 +0,0 @@
-/*
- * decomposition.c
- *
- * Kees Wapenaar "Reciprocity properties of one-way propagators"
- * GEOPHYSICS, VOL. 63, NO. 4 (JULY-AUGUST 1998); P. 1795–1798
- *
- * Created by Jan Thorbecke on 17/03/2014.
- * Copyright 2014 TU Delft. All rights reserved.
- *
- */
-
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <errno.h>
-#include <math.h>
-#include <string.h>
-
-#define MAX(x,y) ((x) > (y) ? (x) : (y))
-#define MIN(x,y) ((x) < (y) ? (x) : (y))
-
-#ifndef COMPLEX
-typedef struct _complexStruct { /* complex number */
- float r,i;
-} complex;
-typedef struct _dcomplexStruct { /* complex number */
- double r,i;
-} dcomplex;
-#endif/* complex */
-
-complex firoot(float x, float stab)
-complex ciroot(complex x, float stab);
-complex cwp_csqrt(complex z);
-
-void decud(float om, float rho, float cp, float dx, int nkx, complex *pu);
-
-void kxwfilter(complex *data, float k, float dx, int nkx,
- float alfa1, float alfa2, float perc);
-
-void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
- int nkx, float dx, int nt, float dt, float fmin, float fmax,
- float cp, float rho)
-{
- int iom, iomin, iomax, ikx, nfreq;
- float omin, omax, deltom, om, kp, df;
- float alpha, eps;
- complex *pu, w;
- complex ax, az;
-
- df = 1.0/((float)nt*dt);
- deltom = 2.*M_PI*df;
- omin = 2.*M_PI*fmin;
- omax = 2.*M_PI*fmax;
- nfreq = nt/2+1;
- eps = 0.01;
- alpha = 0.1;
-
- iomin = (int)MIN((omin/deltom), (nfreq-1));
- iomin = MAX(iomin, 1);
- iomax = MIN((int)(omax/deltom), (nfreq-1));
-
- pu = (complex *)malloc(nkx*sizeof(complex));
-
- for (iom = iomin; iom <= iomax; iom++) {
- om = iom*deltom;
-
- decud(om, rho, cp, dx, nkx, alpha, eps, pu);
-
-/*
- kxwfilter(dpux, kp, dx, nkx, alfa1, alfa2, perc);
- kxwfilter(dpuz, kp, dx, nkx, alfa1, alfa2, perc);
-*/
-
- for (ikx = 0; ikx < nkx; ikx++) {
- ax.r = 0.5*rp[iom*nkx+ikx].r;
- ax.i = 0.5*rp[iom*nkx+ikx].i;
- az.r = 0.5*(rvz[iom*nkx+ikx].r*pu[ikx].r - rvz[iom*nkx+ikx].i*pu[ikx].i);
- az.i = 0.5*(rvz[iom*nkx+ikx].i*pu[ikx].r + rvz[iom*nkx+ikx].r*pu[ikx].i);
-
- down[iom*nkx+ikx].r = ax.r + az.r;
- down[iom*nkx+ikx].i = ax.i + az.i;
- up[iom*nkx+ikx].r = ax.r - az.r;
- up[iom*nkx+ikx].i = ax.i - az.i;
- }
-
- }
-
- free(pu);
-
- return;
-}
-
-void decud(float om, float rho, float cp, float dx, int nkx, float alpha, float eps, complex *pu);
-{
- int ikx, ikxmax1, ikxmax2, filterpoints, filterppos;
- float mu, kp, kp2, ks, ks2, ksk;
- float kx, kx2, kzp2, kzs2, dkx;
- float kxfmax, kxnyq, kpos, kneg, alfa, kfilt, perc, band, *filter;
- complex kzp, kzs, cste, ckp, ckp2, ckzp2;
-
-/* with complex frequency
- wom.r=om;
- wom.i=alpha;
-
- ckp.r = wom.r/cp;
- ckp.i = wom.i/cp;
- ckp2.r = ckp.r*ckp.r-ckp.i*ckp.i;
- ckp2.i = 2.0*ckp.r*ckp.i;
- stab = eps*eps*(ckp.r*ckp.r+ckp.i*ckp.i);
-*/
-
-
- kp = om/cp;
- kp2 = kp*kp;
- dkx = 2.0*M_PI/(nkx*dx);
- stab = eps*eps*kp*kp;
-
- /* make kw filter at maximum angle alfa */
- alfa = 90.0;
- perc = 0.10; /* percentage of band to use for smooth filter */
- filter = (float *)malloc(nkx*sizeof(float));
- kpos = kp*sin(M_PI*alfa/180.0);
- kneg = -kpos;
- kxnyq = M_PI/dx;
- if (kpos > kxnyq) kpos = kxnyq;
- band = kpos;
- filterpoints = (int)fabs((int)(perc*band/dkx));
- kfilt = fabs(dkx*filterpoints);
- if (kpos+kfilt < kxnyq) {
- kxfmax = kpos+kfilt;
- filterppos = filterpoints;
- }
- else {
- kxfmax = kxnyq;
- filterppos = (int)(0.15*nkx/2);
- }
- ikxmax1 = (int) (kxfmax/dkx);
- ikxmax2 = ikxmax1 - filterppos;
- // fprintf(stderr,"ikxmax1=%d ikxmax2=%d nkp=%d nkx=%d\n", ikxmax1, ikxmax2, (int)(kp/dkx), nkx);
-
- for (ikx = 0; ikx < ikxmax2; ikx++)
- filter[ikx]=1.0;
- for (ikx = ikxmax2; ikx < ikxmax1; ikx++)
- filter[ikx] =(cos(M_PI*(ikx-ikxmax2)/(ikxmax1-ikxmax2))+1)/2.0;
- for (ikx = ikxmax1; ikx <= nkx/2; ikx++)
- filter[ikx] = 0.0;
- /* end of kxfilter */
-
- for (ikx = 0; ikx <= (nkx/2); ikx++) {
- kx = ikx*dkx;
- kx2 = kx*kx;
- kzp2 = kp2 - kx2;
- kzp = firoot(kzp2);
-
-/* with complex frequency
- kzp2.r = kp2.r - kx2;
- kzp2.i = kp2.i;
- kzp = ciroot(kzp2, stab);
-*/
- if (kzp2 != 0) {
- pu[ikx].r = filter[ikx]*om*rho*kzp.r;
- pu[ikx].i = filter[ikx]*om*rho*kzp.i;
-// pu[ikx].r = om*rho*kzp.r;
-// pu[ikx].i = om*rho*kzp.i;
- }
- else {
- pu[ikx].r = 0.0;
- pu[ikx].i = 0.0;
- }
- }
-
-/* operators are symmetric in kx-w domain */
- for (ikx = (nkx/2+1); ikx < nkx; ikx++) {
- pu[ikx] = pu[nkx-ikx];
- }
- free(filter);
-
- return;
-}
-
-/* compute 1/x */
-complex firoot(float x, float stab)
-{
- complex z;
- if (x > 0.0) {
- z.r = 1.0/sqrt(x+stab);
- z.i = 0.0;
- }
- else if (x == 0.0) {
- z.r = 0.0;
- z.i = 0.0;
- }
- else {
- z.r = 0.0;
- z.i = 1.0/sqrt(-x+stab);
- }
- return z;
-}
-
-complex ciroot(complex x, float stab)
-{
- complex z, kz, kzz;
- float kd;
-
- if (x.r == 0.0) {
- z.r = 0.0;
- z.i = 0.0;
- }
- else {
- kzz = cwp_csqrt(x);
- kz.r = kzz.r;
- kz.i = -abs(kzz.i);
- kd = kz.r*kz.r+kz.i*kz.i+stab;
- z.r = kz.r/kd;
- z.i = -kz.i/kd;
- }
- return z;
-}
-
-complex cwp_csqrt(complex z)
-{
- complex c;
- float x,y,w,r;
- if (z.r==0.0 && z.i==0.0) {
- c.r = c.i = 0.0;
- return c;
- } else {
- x = fabs(z.r);
- y = fabs(z.i);
- if (x>=y) {
- r = y/x;
- w = sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
- } else {
- r = x/y;
- w = sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
- }
- if (z.r>=0.0) {
- c.r = w;
- c.i = z.i/(2.0*w);
- } else {
- c.i = (z.i>=0.0) ? w : -w;
- c.r = z.i/(2.0*c.i);
- }
- return c;
- }
-}
-
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 899f78c67894af737539e63a513f1061351f4fee..66b19a5a96540be86a800298eddf1e1586e3bfc1 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -138,6 +138,7 @@ char *sdoc[] = {
" dzshot=0 .......... if nshot > 1: z-shift in shot locations",
" xsrca= ............ defines source array x-positions",
" zsrca= ............ defines source array z-positions",
+" src_txt=........... text file with source coordinates. Col 1: x, Col. 2: z",
" wav_random=1 ...... 1 generates (band limited by fmax) noise signatures ",
" fmax=from_src ..... maximum frequency in wavelet",
" src_multiwav=0 .... use traces in file_src as areal source",
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index bef3b7952088e86f67b4f549ba02e8b08f941b16..7c6527079281c77d0e0110c7a500bf48f0de3844 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -39,7 +39,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
{
int isnapmax1, isnapmax2, isnapmax, sna_nrsna;
int n1, n2, nx, nz, nsrc, ix, axis, ioPz, is0, optn;
- int idzshot, idxshot;
+ int idzshot, idxshot, nsrctext;
int src_ix0, src_iz0, src_ix1, src_iz1;
int disable_check;
float cp_min, cp_max, cs_min, cs_max, ro_min, ro_max;
@@ -63,6 +63,8 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
int is,ntraces,length_random;
float rand;
char *src_positions, tmpname[1024];
+ char* src_txt;
+ FILE *fp;
if (!getparint("verbose",&verbose)) verbose=0;
if (!getparint("disable_check",&disable_check)) disable_check=0;
@@ -584,6 +586,9 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
verr("Number of sources in array xsrca (%d), zsrca(%d) are not equal",nxsrc, nzsrc);
}
+ /* source positions defined through txt file */
+ if (!getparstring("src_txt",&src_txt)) src_txt=NULL;
+
/* check if sources on a circle are defined */
if (getparfloat("rsrc", &rsrc)) {
@@ -729,17 +734,45 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
}
}
- else if (nxsrc != 0) {
+ else if ( (nxsrc != 0) || (src_txt != NULL) ) {
/* source array is defined */
- nsrc=nxsrc;
+ if (src_txt!=NULL) {
+ /* Sources from a Text File */
+ /* Open text file */
+ nsrctext=0;
+ fp=fopen(src_txt,"r");
+ assert(fp!=NULL);
+ /* Get number of lines */
+ while (!feof(fp)) if (fgetc(fp)=='\n') nsrctext++;
+ fseek(fp,-1,SEEK_CUR);
+ if (fgetc(fp)!='\n') nsrctext++; /* Checks if last line terminated by /n */
+ if (verbose) vmess("Number of sources in src_txt file: %d",nsrctext);
+ rewind(fp);
+ nsrc=nsrctext;
+ }
+ else {
+ nsrc=nxsrc;
+ }
+ /* Allocate arrays */
src->x = (int *)malloc(nsrc*sizeof(int));
src->z = (int *)malloc(nsrc*sizeof(int));
src->tbeg = (float *)malloc(nsrc*sizeof(float));
src->tend = (float *)malloc(nsrc*sizeof(float));
xsrca = (float *)malloc(nsrc*sizeof(float));
zsrca = (float *)malloc(nsrc*sizeof(float));
- getparfloat("xsrca", xsrca);
- getparfloat("zsrca", zsrca);
+ if (src_txt!=NULL) {
+ /* Read in source coordinates */
+ for (i=0;i<nsrc;i++) {
+ if (fscanf(fp,"%e %e\n",&xsrca[i],&zsrca[i])!=2) vmess("Source Text File: Can not parse coordinates on line %d.",i);
+ }
+ /* Close file */
+ fclose(fp);
+ }
+ else {
+ getparfloat("xsrca", xsrca);
+ getparfloat("zsrca", zsrca);
+ }
+ /* Process coordinates */
for (is=0; is<nsrc; is++) {
src->x[is] = NINT((xsrca[is]-sub_x0)/dx);
src->z[is] = NINT((zsrca[is]-sub_z0)/dz);
@@ -747,6 +780,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
src->tend[is] = (wav->nt-1)*wav->dt;
if (verbose>3) fprintf(stderr,"Source Array: xsrc[%d]=%f zsrc=%f\n", is, xsrca[is], zsrca[is]);
}
+
src->random = 1;
wav->nsamp = (size_t *)malloc((nsrc+1)*sizeof(size_t));
if (wav->random) {
@@ -1064,7 +1098,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
rec->skipdt=NINT(dtrcv/dt);
dtrcv = mod->dt*rec->skipdt;
if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
- if (!getparint("rec_ntsam",&rec->nt)) rec->nt=NINT((mod->tmod)/dtrcv)+1;
+ if (!getparint("rec_ntsam",&rec->nt)) rec->nt=NINT((mod->tmod-rdelay)/dtrcv)+1;
if (!getparint("rec_int_p",&rec->int_p)) rec->int_p=0;
if (!getparint("rec_int_vx",&rec->int_vx)) rec->int_vx=0;
if (!getparint("rec_int_vz",&rec->int_vz)) rec->int_vz=0;
@@ -1072,7 +1106,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparint("scale",&rec->scale)) rec->scale=0;
if (!getparfloat("dxspread",&dxspread)) dxspread=0;
if (!getparfloat("dzspread",&dzspread)) dzspread=0;
- rec->nt=MIN(rec->nt, NINT((mod->tmod)/dtrcv)+1);
+ rec->nt=MIN(rec->nt, NINT((mod->tmod-rdelay)/dtrcv)+1);
/* allocation of receiver arrays is done in recvPar */
/*
diff --git a/fdelmodc/writeSnapTimes.c.org b/fdelmodc/writeSnapTimes.c.org
deleted file mode 100644
index 69120ea518a619ede9580aec09e0c92f8f42a26c..0000000000000000000000000000000000000000
--- a/fdelmodc/writeSnapTimes.c.org
+++ /dev/null
@@ -1,172 +0,0 @@
-#define _FILE_OFFSET_BITS 64
-#define _LARGEFILE_SOURCE
-#define _LARGEFILE64_SOURCE
-
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <errno.h>
-#include <math.h>
-#include <string.h>
-#include "par.h"
-#include "segy.h"
-#include "fdelmodc.h"
-
-/**
-* Writes gridded wavefield(s) at a desired time to output file(s)
-*
-* AUTHOR:
-* Jan Thorbecke (janth@xs4all.nl)
-* The Netherlands
-**/
-
-
-FILE *fileOpen(char *file, char *ext, int append);
-int traceWrite(segy *hdr, float *data, int n, FILE *fp);
-
-#define MAX(x,y) ((x) > (y) ? (x) : (y))
-#define MIN(x,y) ((x) < (y) ? (x) : (y))
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-
-int writeSnapTimes(modPar mod, snaPar sna, int ixsrc, int izsrc, int itime, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose)
-{
- FILE *fpvx, *fpvz, *fptxx, *fptzz, *fptxz, *fpp, *fppp, *fpss;
- int append, isnap;
- int n1, ibnd, ixs, izs, ize, i, j;
- int ix, iz, ix2, iz2;
- float *snap, sdx;
- segy hdr;
-
- if (sna.nsnap==0) return 0;
-
- ibnd = mod.iorder/2-1;
- n1 = mod.naz;
- sdx = 1.0/mod.dx;
-
- /* check if this itime is a desired snapshot time */
- if ( (((itime-sna.delay) % sna.skipdt)==0) &&
- (itime >= sna.delay) &&
- (itime <= sna.delay+(sna.nsnap-1)*sna.skipdt) ) {
-
- isnap = NINT((itime-sna.delay)/sna.skipdt);
- if (verbose) vmess("Writing snapshot(%d) at time=%.3f", isnap+1, itime*mod.dt);
-
- if (isnap) append=1;
- else append=0;
-
- if (sna.type.vx) fpvx = fileOpen(sna.file_snap, "_svx", append);
- if (sna.type.vz) fpvz = fileOpen(sna.file_snap, "_svz", append);
- if (sna.type.p) fpp = fileOpen(sna.file_snap, "_sp", append);
- if (sna.type.txx) fptxx = fileOpen(sna.file_snap, "_stxx", append);
- if (sna.type.tzz) fptzz = fileOpen(sna.file_snap, "_stzz", append);
- if (sna.type.txz) fptxz = fileOpen(sna.file_snap, "_stxz", append);
- if (sna.type.pp) fppp = fileOpen(sna.file_snap, "_spp", append);
- if (sna.type.ss) fpss = fileOpen(sna.file_snap, "_sss", append);
-
- memset(&hdr,0,TRCBYTES);
- hdr.dt = 1000000*(mod.dt);
- hdr.scalco = -1000;
- hdr.scalel = -1000;
- hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
- hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
- hdr.fldr = isnap+1;
- hdr.trid = 1;
- hdr.ns = sna.nz;
- hdr.trwf = sna.nx;
- hdr.ntr = (isnap+1)*sna.nx;
- hdr.f1 = sna.z1*mod.dz+mod.z0;
- hdr.f2 = sna.x1*mod.dx+mod.x0;
- hdr.d1 = mod.dz*sna.skipdz;
- hdr.d2 = mod.dx*sna.skipdx;
-
-/***********************************************************************
-* vx velocities have one sample less in x-direction
-* vz velocities have one sample less in z-direction
-* txz stresses have one sample less in z-direction and x-direction
-***********************************************************************/
-
- snap = (float *)malloc(sna.nz*sizeof(float));
-
- /* Decimate, with skipdx and skipdz, the number of gridpoints written to file
- and write to file. */
- for (ixs=sna.x1, i=0; ixs<=sna.x2; ixs+=sna.skipdx, i++) {
- hdr.tracf = i+1;
- hdr.tracl = isnap*sna.nx+i+1;
- hdr.gx = 1000*(mod.x0+ixs*mod.dx);
- ix = ixs+ibnd;
- ix2 = ix+1;
-
- izs = sna.z1+ibnd;
- ize = sna.z2+ibnd;
- if (sna.type.vx) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = vx[ix2*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fpvx);
- }
- if (sna.type.vz) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = vz[ix*n1+iz+1];
- }
- traceWrite(&hdr, snap, sna.nz, fpvz);
- }
- if (sna.type.p) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = tzz[ix*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fpp);
- }
- if (sna.type.tzz) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = tzz[ix*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fptzz);
- }
- if (sna.type.txx) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = txx[ix*n1+iz];
- }
- traceWrite(&hdr, snap, sna.nz, fptxx);
- }
- if (sna.type.txz) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- snap[j] = txz[ix2*n1+iz+1];
- }
- traceWrite(&hdr, snap, sna.nz, fptxz);
- }
- /* calculate divergence of velocity field */
- if (sna.type.pp) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- iz2 = iz+1;
- snap[j] = sdx*((vx[ix2*n1+iz]-vx[ix*n1+iz])+
- (vz[ix*n1+iz2]-vz[ix*n1+iz]));
- }
- traceWrite(&hdr, snap, sna.nz, fppp);
- }
- /* calculate rotation of velocity field */
- if (sna.type.ss) {
- for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
- iz2 = iz+1;
- snap[j] = sdx*((vx[ix2*n1+iz2]-vx[ix2*n1+iz])-
- (vz[ix2*n1+iz2]-vz[ix*n1+iz2]));
- }
- traceWrite(&hdr, snap, sna.nz, fpss);
- }
-
- }
-
- if (sna.type.vx) fclose(fpvx);
- if (sna.type.vz) fclose(fpvz);
- if (sna.type.p) fclose(fpp);
- if (sna.type.txx) fclose(fptxx);
- if (sna.type.tzz) fclose(fptzz);
- if (sna.type.txz) fclose(fptxz);
- if (sna.type.pp) fclose(fppp);
- if (sna.type.ss) fclose(fpss);
-
- free(snap);
- }
-
- return 0;
-}
-
diff --git a/marchenko_applications/fmute b/marchenko_applications/fmute
deleted file mode 100755
index d035dc3d1e8d3256fb2f6b0a7ffa3429f5681361..0000000000000000000000000000000000000000
Binary files a/marchenko_applications/fmute and /dev/null differ
diff --git a/marchenko_applications/marchenko b/marchenko_applications/marchenko
deleted file mode 100755
index 33fa897e36ac1c3d0d02bba1e71bba3e3a3c7f6a..0000000000000000000000000000000000000000
Binary files a/marchenko_applications/marchenko and /dev/null differ
diff --git a/marchenko_full/fmute b/marchenko_full/fmute
deleted file mode 100755
index c87afaa36d11069199cf8405451d18a19a98b6fd..0000000000000000000000000000000000000000
Binary files a/marchenko_full/fmute and /dev/null differ
diff --git a/utils/freqwave.c b/utils/freqwave.c
index 156299e18ff222943be57597663f57e906da725a..d5c894c7481eeee668e0959ebe7bdfb0f971ebd2 100644
--- a/utils/freqwave.c
+++ b/utils/freqwave.c
@@ -33,470 +33,468 @@ void hilbertTrans(float *data, int nsam);
void freqwave(float *wave, int nt, float dt, float fp, float fmin, float flef, float frig, float fmax, float t0, float db, int shift, int cm, int cn, char *w, float scale, int scfft, int inverse, float eps, int verbose)
{
- int it, iof, nfreq, nf, i, j, sign, optn, stored;
- int ifmin1, ifmin2, ifmax1, ifmax2;
- float df, fact, alfa, f, max, freq, att, ampl, phase;
- float tt, dum;
- float *rwave, *amplitude;
- complex *cwave, tmp, *mpwave;
-
- optn = optncr(nt);
- nfreq = 1+(optn/2);
- df = 1.0/(dt*optn);
- iof = MAX(NINT(fmax/df), NINT(fp/df));
- att = pow(10.0, db/20.0);
-
- if (iof > nfreq) verr("characterizing frequency aliased");
-
- cwave = (complex *)malloc(nfreq*sizeof(complex));
- rwave = (float *)malloc((optn+2)*sizeof(float));
-
- stored = 0;
-
- if (strstr(w, "g0") != NULL) {
- i = NINT(fmax/df);
- for (iof = i; iof > 0; iof--) {
- f = iof*df;
- if((gauss0freq(fmax, f) < att)&&(stored != 1)) {
- freq = f;
- stored = 1;
- }
- }
- if (stored == 0) verr("No valid wavelet found.");
- stored = 0;
- if (shift == 1) {
- for (i = 0; i < optn; i++) {
- if ((fabs(gauss0time((float)i*dt,freq))<1e-3)&&(stored != 1)) {
- t0 = (float)i*dt;
- stored = 1;
- }
- }
- }
- for (iof = 0; iof < nfreq; iof++) {
- f = iof*df;
- fact = f*f/(freq*freq);
- fact = exp(-fact);
- cwave[iof].r = fact*cos(2.0*M_PI*f*t0);
- cwave[iof].i = -fact*sin(2.0*M_PI*f*t0);
- }
- if (verbose >= 1) {
- vmess("Gaussian wavelet");
- vmess("----------------");
- vmess("Number of time samples .. = %d", nt);
- vmess("time step ............... = %f (s)", dt);
- vmess("maximum frequency at ... = %f (Hz)",fmax);
- vmess("with attenutation ....... = %f", att);
- vmess("time shift .............. = %f (s)", t0);
- }
- }
- else if (strstr(w, "g1") != NULL) {
- if (fp < 0.0) {
- i = NINT(fmax/df);
- for (iof = i; iof > 0; iof--) {
- f = iof*df;
- if((gauss1freq(fmax, f) < att)&&(stored != 1)) {
- freq = f;
- stored = 1;
- }
- }
- if (stored == 0) verr("No valid wavelet found.");
- }
- else freq = fp;
- alfa = sqrt(2.0)*freq;
- stored = 0;
-
- if (shift == 1) {
- for (i = 1; i < optn; i++) {
- tt=(float)i*dt;
- dum = fabs(gauss1time(tt,freq))+fabs(gauss1time((tt+dt),freq));
- if ((dum<1e-4)&&(stored != 1)) {
- t0 = (float)i*dt;
- stored = 1;
- }
- }
- }
-
- for (iof = 0; iof < nfreq; iof++) {
- f = iof*df;
- fact = f*f/(alfa*alfa);
- fact = f*exp(-fact)/alfa;
- cwave[iof].r = fact*sin(2.0*M_PI*f*t0);
- cwave[iof].i = fact*cos(2.0*M_PI*f*t0);
- }
- if (verbose >= 1) {
- vmess("Derivative of Gaussian wavelet");
- vmess("------------------------------");
- vmess("Number of time samples .. = %d",nt);
- vmess("time step ............... = %f (s)",dt);
- if (fp < 0) {
- vmess("maximum frequency at ... = %f (Hz)",fmax);
- vmess("with attenutation ....... = %f", att);
- }
- vmess("frequency peak at ....... = %f Hz", freq);
- vmess("time shift .............. = %f (s)", t0);
- }
- }
- else if (strstr(w, "cs") != NULL) {
- freq = acos((float)(cm-cn)/(float)(cm+cn))/(2.0*M_PI*dt);
- fact = 1.0/(cos(freq*2.0*M_PI*dt)+sin(freq*2.0*M_PI*dt));
- if (shift == 1) t0 = (cn+cm)*dt;
-
- for (iof = 0; iof < nfreq; iof++) {
- f = 2.0*M_PI*iof/(float)nt;
- ampl = pow((1.0-cos(f)), cn/2.0)*pow((1.0+cos(f)), cm/2.0);
- phase = atan((fact*sin(f))/(1.0+fact*cos(f)));
- cwave[iof].r = ampl*cos(phase-f*nt*t0);
- cwave[iof].i = ampl*sin(phase-f*nt*t0);
- }
- if (verbose >= 1) {
- vmess("Neidell Type of wavelet");
- vmess("-----------------------");
- vmess("Number of time samples .. = %d", nt);
- vmess("time step ............... = %f (s)", dt);
- vmess("frequency peak at ....... = %f Hz", freq);
- vmess("time shift .............. = %f (s)", t0);
- }
- }
- else if (strstr(w, "fw") != NULL) {
- ifmin1 = (int) (fmin/df);
- ifmin2 = (int) (flef/df);
- ifmax2 = (int) (frig/df);
- ifmax1 = (int) (fmax/df);
- for (j = 0; j < ifmin1; j++) {
- cwave[j].r = 0.0;
- cwave[j].i = 0.0;
- }
- for (j = ifmin1; j < ifmin2; j++) {
- cwave[j].r = (cos(M_PI*(j-ifmin2)/(ifmin1-ifmin2))+1.0)/2.0;
- cwave[j].i = 0.0;
- }
- for (j = ifmin2; j < ifmax2; j++) {
- cwave[j].r = 1.0;
- cwave[j].i = 0.0;
- }
- for (j = ifmax2; j < ifmax1; j++) {
- cwave[j].r =(cos(M_PI*(j-ifmax2)/(ifmax1-ifmax2))+1.0)/2.0;
- cwave[j].i = 0.0;
- }
- for (j = ifmax1; j < nfreq; j++) {
- cwave[j].r = 0.0;
- cwave[j].i = 0.0;
- }
- if (shift == 1) {
- sign = 1;
- cr1fft(cwave, rwave, optn, sign);
+ int it, iof, nfreq, nf, i, j, sign, optn, stored;
+ int ifmin1, ifmin2, ifmax1, ifmax2;
+ float df, fact, alfa, f, max, freq, att, ampl, phase;
+ float tt, dum;
+ float *rwave, *amplitude;
+ complex *cwave, tmp, *mpwave;
+
+ optn = optncr(nt);
+ nfreq = 1+(optn/2);
+ df = 1.0/(dt*optn);
+ iof = MAX(NINT(fmax/df), NINT(fp/df));
+ att = pow(10.0, db/20.0);
+
+ if (iof > nfreq) verr("characterizing frequency aliased");
+
+ cwave = (complex *)malloc(nfreq*sizeof(complex));
+ rwave = (float *)malloc((optn+2)*sizeof(float));
+
+ stored = 0;
+
+ if (strstr(w, "g0") != NULL) {
+ i = NINT(fmax/df);
+ for (iof = i; iof > 0; iof--) {
+ f = iof*df;
+ if((gauss0freq(fmax, f) < att)&&(stored != 1)) {
+ freq = f;
+ stored = 1;
+ }
+ }
+ if (stored == 0) verr("No valid wavelet found.");
+ stored = 0;
+ if (shift == 1) {
+ for (i = 0; i < optn; i++) {
+ if ((fabs(gauss0time((float)i*dt,freq))<1e-3)&&(stored != 1)) {
+ t0 = (float)i*dt;
+ stored = 1;
+ }
+ }
+ }
+ for (iof = 0; iof < nfreq; iof++) {
+ f = iof*df;
+ fact = f*f/(freq*freq);
+ fact = exp(-fact);
+ cwave[iof].r = fact*cos(2.0*M_PI*f*t0);
+ cwave[iof].i = -fact*sin(2.0*M_PI*f*t0);
+ }
+ if (verbose >= 1) {
+ vmess("Gaussian wavelet");
+ vmess("----------------");
+ vmess("Number of time samples .. = %d", nt);
+ vmess("time step ............... = %f (s)", dt);
+ vmess("maximum frequency at ... = %f (Hz)",fmax);
+ vmess("with attenutation ....... = %f", att);
+ vmess("time shift .............. = %f (s)", t0);
+ }
+ }
+ else if (strstr(w, "g1") != NULL) {
+ if (fp < 0.0) {
+ i = NINT(fmax/df);
+ for (iof = i; iof > 0; iof--) {
+ f = iof*df;
+ if((gauss1freq(fmax, f) < att)&&(stored != 1)) {
+ freq = f;
+ stored = 1;
+ }
+ }
+ if (stored == 0) verr("No valid wavelet found.");
+ }
+ else freq = fp;
+ alfa = sqrt(2.0)*freq;
+ stored = 0;
+
+ if (shift == 1) {
+ for (i = 1; i < optn; i++) {
+ tt=(float)i*dt;
+ dum = fabs(gauss1time(tt,freq))+fabs(gauss1time((tt+dt),freq));
+ if ((dum<1e-4)&&(stored != 1)) {
+ t0 = (float)i*dt;
+ stored = 1;
+ }
+ }
+ }
+
+ for (iof = 0; iof < nfreq; iof++) {
+ f = iof*df;
+ fact = f*f/(alfa*alfa);
+ fact = f*exp(-fact)/alfa;
+ cwave[iof].r = fact*sin(2.0*M_PI*f*t0);
+ cwave[iof].i = fact*cos(2.0*M_PI*f*t0);
+ }
+ if (verbose >= 1) {
+ vmess("Derivative of Gaussian wavelet");
+ vmess("------------------------------");
+ vmess("Number of time samples .. = %d",nt);
+ vmess("time step ............... = %f (s)",dt);
+ if (fp < 0) {
+ vmess("maximum frequency at ... = %f (Hz)",fmax);
+ vmess("with attenutation ....... = %f", att);
+ }
+ vmess("frequency peak at ....... = %f Hz", freq);
+ vmess("time shift .............. = %f (s)", t0);
+ }
+ }
+ else if (strstr(w, "cs") != NULL) {
+ freq = acos((float)(cm-cn)/(float)(cm+cn))/(2.0*M_PI*dt);
+ fact = 1.0/(cos(freq*2.0*M_PI*dt)+sin(freq*2.0*M_PI*dt));
+ if (shift == 1) t0 = (cn+cm)*dt;
+
+ for (iof = 0; iof < nfreq; iof++) {
+ f = 2.0*M_PI*iof/(float)nt;
+ ampl = pow((1.0-cos(f)), cn/2.0)*pow((1.0+cos(f)), cm/2.0);
+ phase = atan((fact*sin(f))/(1.0+fact*cos(f)));
+ cwave[iof].r = ampl*cos(phase-f*nt*t0);
+ cwave[iof].i = ampl*sin(phase-f*nt*t0);
+ }
+ if (verbose >= 1) {
+ vmess("Neidell Type of wavelet");
+ vmess("-----------------------");
+ vmess("Number of time samples .. = %d", nt);
+ vmess("time step ............... = %f (s)", dt);
+ vmess("frequency peak at ....... = %f Hz", freq);
+ vmess("time shift .............. = %f (s)", t0);
+ }
+ }
+ else if (strstr(w, "fw") != NULL) {
+ ifmin1 = (int) (fmin/df);
+ ifmin2 = (int) (flef/df);
+ ifmax2 = (int) (frig/df);
+ ifmax1 = (int) (fmax/df);
+ for (j = 0; j < ifmin1; j++) {
+ cwave[j].r = 0.0;
+ cwave[j].i = 0.0;
+ }
+ for (j = ifmin1; j < ifmin2; j++) {
+ cwave[j].r = (cos(M_PI*(j-ifmin2)/(ifmin1-ifmin2))+1.0)/2.0;
+ cwave[j].i = 0.0;
+ }
+ for (j = ifmin2; j < ifmax2; j++) {
+ cwave[j].r = 1.0;
+ cwave[j].i = 0.0;
+ }
+ for (j = ifmax2; j < ifmax1; j++) {
+ cwave[j].r =(cos(M_PI*(j-ifmax2)/(ifmax1-ifmax2))+1.0)/2.0;
+ cwave[j].i = 0.0;
+ }
+ for (j = ifmax1; j < nfreq; j++) {
+ cwave[j].r = 0.0;
+ cwave[j].i = 0.0;
+ }
+ if (shift == 1) {
+ sign = 1;
+ cr1fft(cwave, rwave, optn, sign);
max = rwave[0]*rwave[0];
- for (i = 1; i < optn; i++) {
+ for (i = 1; i < optn; i++) {
if (rwave[i] < 0.0) { it=i; break;}
}
it = it;
dum = 0.0;
- for (j = 0; j < it; j++)
+ for (j = 0; j < it; j++)
dum += fabs(rwave[j]);
max = dum;
- for (i = 1; i < optn/2; i++) {
+ for (i = 1; i < optn/2; i++) {
dum = dum - fabs(rwave[i-1]) + fabs(rwave[i+it-1]);
- //fprintf(stderr,"dum=%f i=%d\n", dum, i);
- if ((dum)>(0.03*max*att/it)) {
- t0 = (float)i*dt;
- }
- }
- }
- for (iof = 0; iof < nfreq; iof++) {
- f = iof*df;
- tmp.r = cwave[iof].r*cos(2.0*M_PI*f*t0);
- tmp.i = -cwave[iof].r*sin(2.0*M_PI*f*t0);
- cwave[iof].r = tmp.r;
- cwave[iof].i = tmp.i;
+ if ((dum)>(0.03*max*att/it)) {
+ t0 = (float)i*dt;
+ }
+ }
+ }
+ for (iof = 0; iof < nfreq; iof++) {
+ f = iof*df;
+ tmp.r = cwave[iof].r*cos(2.0*M_PI*f*t0);
+ tmp.i = -cwave[iof].r*sin(2.0*M_PI*f*t0);
+ cwave[iof].r = tmp.r;
+ cwave[iof].i = tmp.i;
/* older version has multiplication with dt changed in april 2014
- cwave[iof].r = dt*tmp.r;
- cwave[iof].i = dt*tmp.i;
-*/
- }
-
- if (verbose >= 1) {
- vmess("Flat spectrum wavelet");
- vmess("---------------------");
- vmess("Number of time samples .. = %d", nt);
- vmess("time step ............... = %f (s)", dt);
- vmess("maximum frequency ....... = %f Hz", fmax);
- vmess("left cut-off frequency .. = %f Hz", flef);
- vmess("right cut-off frequency . = %f Hz", frig);
- vmess("minimum frequency ....... = %f Hz", fmin);
- vmess("time shift .............. = %f (s)", t0);
- }
-
- }
- else if (strstr(w, "mon") != NULL) {
- for (j = 0; j < nfreq; j++) {
- cwave[j].r = 0.0;
- cwave[j].i = 0.0;
- }
- i = NINT(fp/df);
- cwave[i].r = 0.5*cos(2.0*M_PI*i*df*t0);
- cwave[i].i = -0.5*sin(2.0*M_PI*i*df*t0);
-
- if (verbose >= 1) {
- vmess("Monochromatic wavelet");
- vmess("---------------------");
- vmess("Number of time samples .. = %d", nt);
- vmess("time step ............... = %e (s)", dt);
- vmess("frequency ............... = %f Hz", i*df);
- vmess("time shift .............. = %e (s)", t0);
- }
- }
- else if (strstr(w, "sqrtg2") != NULL) {
- if (fp < 0.0) {
- i = NINT(fmax/(2.0*df));
- for (iof = i; iof > 0; iof--) {
- f = iof*df;
- if((gauss2freq(fmax, f) < att)&&(stored != 1)) {
- freq = f;
- stored = 1;
- }
- }
- if (stored == 0) verr("No valid wavelet found.");
- }
- else freq = fp;
- stored = 0;
-
- if (shift == 1) {
- for (i = 0; i < optn; i++) {
- tt=(float)i*dt;
- dum = fabs(gauss2time(tt,freq))+fabs(gauss2time((tt+dt),freq));
- if ((dum<1e-3)&&(stored != 1)) {
- t0 = (float)i*dt;
- stored = 1;
- }
- }
- }
-
- for (iof = 0; iof < nfreq; iof++) {
- f = iof*df;
- fact = f/(freq);
- fact *= exp(-0.5*fact*fact);
- cwave[iof].r = fact*cos(2.0*M_PI*f*t0);
- cwave[iof].i = -fact*sin(2.0*M_PI*f*t0);
- }
- if (verbose >= 1) {
- vmess("Sqrt of Second derivative of Gaussian wavelet");
- vmess("-------------------------------------");
- vmess("Number of time samples .. = %d", nt);
- vmess("time step ............... = %f (s)", dt);
- if (fp < 0) {
- vmess("maximum frequency at ... = %f (Hz)",fmax);
- vmess("with attenutation ....... = %f", att);
- }
- vmess("frequency peak at ....... = %f Hz", freq);
- vmess("time shift .............. = %f (s)", t0);
- }
- }
- else {
- if (fp < 0.0) {
- i = NINT(fmax/(2.0*df));
- for (iof = i; iof > 0; iof--) {
- f = iof*df;
- if((gauss2freq(fmax, f) < att)&&(stored != 1)) {
- freq = f;
- stored = 1;
- }
- }
- if (stored == 0) verr("No valid wavelet found.");
- }
- else freq = fp;
- stored = 0;
-
- if (shift == 1) {
- for (i = 0; i < optn; i++) {
- tt=(float)i*dt;
- dum = fabs(gauss2time(tt,freq))+fabs(gauss2time((tt+dt),freq));
- if ((dum<1e-3)&&(stored != 1)) {
- t0 = (float)i*dt;
- stored = 1;
- }
- }
- }
+ cwave[iof].r = dt*tmp.r;
+ cwave[iof].i = dt*tmp.i; */
+ }
+
+ if (verbose >= 1) {
+ vmess("Flat spectrum wavelet");
+ vmess("---------------------");
+ vmess("Number of time samples .. = %d", nt);
+ vmess("time step ............... = %f (s)", dt);
+ vmess("maximum frequency ....... = %f Hz", fmax);
+ vmess("left cut-off frequency .. = %f Hz", flef);
+ vmess("right cut-off frequency . = %f Hz", frig);
+ vmess("minimum frequency ....... = %f Hz", fmin);
+ vmess("time shift .............. = %f (s)", t0);
+ }
+
+ }
+ else if (strstr(w, "mon") != NULL) {
+ for (j = 0; j < nfreq; j++) {
+ cwave[j].r = 0.0;
+ cwave[j].i = 0.0;
+ }
+ i = NINT(fp/df);
+ cwave[i].r = 0.5*cos(2.0*M_PI*i*df*t0);
+ cwave[i].i = -0.5*sin(2.0*M_PI*i*df*t0);
+
+ if (verbose >= 1) {
+ vmess("Monochromatic wavelet");
+ vmess("---------------------");
+ vmess("Number of time samples .. = %d", nt);
+ vmess("time step ............... = %e (s)", dt);
+ vmess("frequency ............... = %f Hz", i*df);
+ vmess("time shift .............. = %e (s)", t0);
+ }
+ }
+ else if (strstr(w, "sqrtg2") != NULL) {
+ if (fp < 0.0) {
+ i = NINT(fmax/(2.0*df));
+ for (iof = i; iof > 0; iof--) {
+ f = iof*df;
+ if((gauss2freq(fmax, f) < att)&&(stored != 1)) {
+ freq = f;
+ stored = 1;
+ }
+ }
+ if (stored == 0) verr("No valid wavelet found.");
+ }
+ else freq = fp;
+ stored = 0;
+
+ if (shift == 1) {
+ for (i = 0; i < optn; i++) {
+ tt=(float)i*dt;
+ dum = fabs(gauss2time(tt,freq))+fabs(gauss2time((tt+dt),freq));
+ if ((dum<1e-3)&&(stored != 1)) {
+ t0 = (float)i*dt;
+ stored = 1;
+ }
+ }
+ }
+
for (iof = 0; iof < nfreq; iof++) {
- f = iof*df;
- fact = f*f/(freq*freq);
- fact *= exp(-fact);
- cwave[iof].r = fact*cos(2.0*M_PI*f*t0);
- cwave[iof].i = -fact*sin(2.0*M_PI*f*t0);
- }
- if (verbose >= 1) {
- vmess("Second derivative of Gaussian wavelet");
- vmess("-------------------------------------");
- vmess("Number of time samples .. = %d", nt);
- vmess("time step ............... = %f (s)", dt);
- if (fp < 0) {
- vmess("maximum frequency at ... = %f (Hz)",fmax);
- vmess("with attenutation ....... = %f", att);
- }
- vmess("frequency peak at ....... = %f Hz", freq);
- vmess("time shift .............. = %f (s)", t0);
- }
- }
+ f = iof*df;
+ fact = f/(freq);
+ fact *= exp(-0.5*fact*fact);
+ cwave[iof].r = fact*cos(2.0*M_PI*f*t0);
+ cwave[iof].i = -fact*sin(2.0*M_PI*f*t0);
+ }
+ if (verbose >= 1) {
+ vmess("Sqrt of Second derivative of Gaussian wavelet");
+ vmess("-------------------------------------");
+ vmess("Number of time samples .. = %d", nt);
+ vmess("time step ............... = %f (s)", dt);
+ if (fp < 0) {
+ vmess("maximum frequency at ... = %f (Hz)",fmax);
+ vmess("with attenutation ....... = %f", att);
+ }
+ vmess("frequency peak at ....... = %f Hz", freq);
+ vmess("time shift .............. = %f (s)", t0);
+ }
+ }
+ else {
+ if (fp < 0.0) {
+ i = NINT(fmax/(2.0*df));
+ for (iof = i; iof > 0; iof--) {
+ f = iof*df;
+ if((gauss2freq(fmax, f) < att)&&(stored != 1)) {
+ freq = f;
+ stored = 1;
+ }
+ }
+ if (stored == 0) verr("No valid wavelet found.");
+ }
+ else freq = fp;
+ stored = 0;
+
+ if (shift == 1) {
+ for (i = 0; i < optn; i++) {
+ tt=(float)i*dt;
+ dum = fabs(gauss2time(tt,freq))+fabs(gauss2time((tt+dt),freq));
+ if ((dum<1e-3)&&(stored != 1)) {
+ t0 = (float)i*dt;
+ stored = 1;
+ }
+ }
+ }
+ for (iof = 0; iof < nfreq; iof++) {
+ f = iof*df;
+ fact = f*f/(freq*freq);
+ fact *= exp(-fact);
+ cwave[iof].r = fact*cos(2.0*M_PI*f*t0);
+ cwave[iof].i = -fact*sin(2.0*M_PI*f*t0);
+ }
+ if (verbose >= 1) {
+ vmess("Second derivative of Gaussian wavelet");
+ vmess("-------------------------------------");
+ vmess("Number of time samples .. = %d", nt);
+ vmess("time step ............... = %f (s)", dt);
+ if (fp < 0) {
+ vmess("maximum frequency at ... = %f (Hz)",fmax);
+ vmess("with attenutation ....... = %f", att);
+ }
+ vmess("frequency peak at ....... = %f Hz", freq);
+ vmess("time shift .............. = %f (s)", t0);
+ }
+ }
if (inverse==1) {
vmess("inverse with eps ....... = %f (s)", eps);
for (iof = 1; iof < nfreq; iof++) {
fact = cwave[iof].r*cwave[iof].r + cwave[iof].i*cwave[iof].i;
- cwave[iof].r = cwave[iof].r/(fact+eps);
- cwave[iof].i = -cwave[iof].i/(fact+eps);
+ cwave[iof].r = cwave[iof].r/(fact+eps);
+ cwave[iof].i = -cwave[iof].i/(fact+eps);
}
cwave[0].r = 0.0;
cwave[0].i = 0.0;
}
- /* minimum phase calculation */
+ /* minimum phase calculation */
if (inverse==2) {
vmess("minimum phase calculation ");
- nf = (2*(nfreq-1));
- mpwave = (complex *)calloc(nf,sizeof(complex));
-
- fprintf(stderr,"nf=%d\n", nf);
- amplitude = (float *)calloc(2*nf,sizeof(float));
+ nf = (2*(nfreq-1));
+ mpwave = (complex *)calloc(nf,sizeof(complex));
+
+ fprintf(stderr,"nf=%d\n", nf);
+ amplitude = (float *)calloc(2*nf,sizeof(float));
for (iof = 0; iof < nfreq; iof++) {
- fact = sqrt(cwave[iof].r*cwave[iof].r + cwave[iof].i*cwave[iof].i);
- if (fact > 0.0) amplitude[iof] = log(fact);
- else amplitude[iof] = 0.0;
+ fact = sqrt(cwave[iof].r*cwave[iof].r + cwave[iof].i*cwave[iof].i);
+ if (fact > 0.0) amplitude[iof] = log(fact);
+ else amplitude[iof] = 0.0;
amplitude[nf+iof] = fact;
- }
- hilbertTrans(amplitude, nf);
+ }
+ hilbertTrans(amplitude, nf);
for (iof = 0; iof < nfreq; iof++) {
- fact = amplitude[nf+iof];
- fprintf(stderr,"amplitude[%d] = %f phase = %f\n", iof, fact, amplitude[iof]);
- if (fact != 0.0) {
- mpwave[iof].r = (float) fact*cos(amplitude[iof]);
- mpwave[iof].i = (float) -fact*sin(amplitude[iof]);
- }
- else {
- mpwave[iof].r=0.0;
- mpwave[iof].i=0.0;
- }
- }
+ fact = amplitude[nf+iof];
+ fprintf(stderr,"amplitude[%d] = %f phase = %f\n", iof, fact, amplitude[iof]);
+ if (fact != 0.0) {
+ mpwave[iof].r = (float) fact*cos(amplitude[iof]);
+ mpwave[iof].i = (float) -fact*sin(amplitude[iof]);
+ }
+ else {
+ mpwave[iof].r=0.0;
+ mpwave[iof].i=0.0;
+ }
+ }
for (iof = nf-1; iof > nfreq; iof--) {
- mpwave[iof].r=mpwave[nf-iof].r;
- mpwave[iof].i=-1.0*mpwave[nf-iof].i;
- }
- cc1fft(mpwave, nf, 1);
- for (i = 0; i < nt; i++) rwave[i] = mpwave[i].r;
-
- free(amplitude);
- free(mpwave);
+ mpwave[iof].r=mpwave[nf-iof].r;
+ mpwave[iof].i=-1.0*mpwave[nf-iof].i;
+ }
+ cc1fft(mpwave, nf, 1);
+ for (i = 0; i < nt; i++) rwave[i] = mpwave[i].r;
+
+ free(amplitude);
+ free(mpwave);
}
- else {
- sign = 1;
- cr1fft(cwave, rwave, optn, sign);
- }
-
- max = rwave[0];
- for (i = 0; i < nt; i++) if (rwave[i] > max) max = rwave[i];
- max = scale/max;
- if (scale == 0) {
- if (scfft == 0) max = 1.0/(float)nt;
- else max = df;
- }
- //fprintf(stderr,"scaling factor back FFT=%e\n", max);
- for (i = 0; i < nt; i++) wave[i]= rwave[i]*max;
-
- free(cwave);
- free(rwave);
-
- return;
+ else {
+ sign = 1;
+ cr1fft(cwave, rwave, optn, sign);
+ }
+
+ max = rwave[0];
+ for (i = 0; i < nt; i++) if (rwave[i] > max) max = rwave[i];
+ max = scale/max;
+ if (scale == 0) {
+ if (scfft == 0) max = 1.0/(float)nt;
+ else max = df;
+ }
+ //fprintf(stderr,"scaling factor back FFT=%e\n", max);
+ for (i = 0; i < nt; i++) wave[i]= rwave[i]*max;
+
+ free(cwave);
+ free(rwave);
+
+ return;
}
float gauss2time(float t, float f)
{
- float value;
+ float value;
- value = ((1.0-2.0*M_PI*M_PI*f*f*t*t)*exp(-M_PI*M_PI*f*f*t*t));
- return value;
+ value = ((1.0-2.0*M_PI*M_PI*f*f*t*t)*exp(-M_PI*M_PI*f*f*t*t));
+ return value;
}
float gauss1time(float t, float f)
{
- float value;
+ float value;
- value = (-t*exp(-M_PI*M_PI*f*f*t*t))*(sqrt(2.0)*f*M_PI*exp(0.5));
- return value;
+ value = (-t*exp(-M_PI*M_PI*f*f*t*t))*(sqrt(2.0)*f*M_PI*exp(0.5));
+ return value;
}
float gauss0time(float t, float f)
{
- float value;
+ float value;
- value = exp(-M_PI*M_PI*f*f*t*t);
- return value;
+ value = exp(-M_PI*M_PI*f*f*t*t);
+ return value;
}
float gauss2freq(float f, float freq)
{
- float value;
+ float value;
- value = f*f/(freq*freq);
- value *= exp(1.0)*exp(-value);
+ value = f*f/(freq*freq);
+ value *= exp(1.0)*exp(-value);
- return value;
+ return value;
}
float gauss1freq(float f, float freq)
{
- float value;
+ float value;
- value = f*f/(2.0*freq*freq);
- value = sqrt(2.0*exp(1))*f*exp(-value)/(sqrt(2.0)*freq);
+ value = f*f/(2.0*freq*freq);
+ value = sqrt(2.0*exp(1))*f*exp(-value)/(sqrt(2.0)*freq);
- return value;
+ return value;
}
float gauss0freq(float f, float freq)
{
- float value;
+ float value;
- value = f*f/(freq*freq);
- value = exp(-value);
+ value = f*f/(freq*freq);
+ value = exp(-value);
- return value;
+ return value;
}
void hilbertTrans(float *data, int nsam)
{
- int optn, j, sign, nfreq;
- float scale;
- complex *cdata;
-
- optn = optncr(nsam);
- nfreq = optn/2+1;
- fprintf(stderr,"Hilbert optn=%d nsam=%d nfreq=%d\n", optn, nsam, nfreq);
- cdata = (complex *)malloc(optn*sizeof(complex));
- if (cdata == NULL) verr("memory allocation error for cdata");
-
- for(j = 0; j < nsam; j++){
- cdata[j].r = data[j];
- cdata[j].i = 0.0;
- }
- for(j = nsam; j < optn; j++){
- cdata[j].r = 0.0;
- cdata[j].i = 0.0;
- }
- sign = -1;
- cc1fft(&cdata[0], optn, sign);
-
- for(j = nfreq; j < optn; j++){
- cdata[j].r = 0.0;
- cdata[j].i = 0.0;
- }
-
- sign = 1;
- cc1fft(&cdata[0], optn, sign);
-
- scale= 1.0/(float)optn;
- for (j = 0 ; j < nsam ; j++) data[j] = cdata[j].i*scale;
-
- free(cdata);
-
- return;
+ int optn, j, sign, nfreq;
+ float scale;
+ complex *cdata;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ fprintf(stderr,"Hilbert optn=%d nsam=%d nfreq=%d\n", optn, nsam, nfreq);
+ cdata = (complex *)malloc(optn*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ for(j = 0; j < nsam; j++){
+ cdata[j].r = data[j];
+ cdata[j].i = 0.0;
+ }
+ for(j = nsam; j < optn; j++){
+ cdata[j].r = 0.0;
+ cdata[j].i = 0.0;
+ }
+ sign = -1;
+ cc1fft(&cdata[0], optn, sign);
+
+ for(j = nfreq; j < optn; j++){
+ cdata[j].r = 0.0;
+ cdata[j].i = 0.0;
+ }
+
+ sign = 1;
+ cc1fft(&cdata[0], optn, sign);
+
+ scale= 1.0/(float)optn;
+ for (j = 0 ; j < nsam ; j++) data[j] = cdata[j].i*scale;
+
+ free(cdata);
+
+ return;
}