From 7b61a45cbc85271bc7d24ac253346382cf4fe9e7 Mon Sep 17 00:00:00 2001
From: JBrackenhoff <J.A.Brackenhoff@tudelft.nl>
Date: Thu, 21 Dec 2017 13:00:27 +0100
Subject: [PATCH] fast
---
marchenko_applications/AmpEst.c | 2 +-
marchenko_applications/Cost.c | 2 +-
marchenko_applications/Makefile | 26 ++++++++++++++++++++-----
marchenko_applications/marchenko.c | 31 ++++++++++++++++--------------
marchenko_full/AmpEst.c | 2 +-
marchenko_full/Cost.c | 2 +-
marchenko_full/JespersRayTracer.c | 4 ++--
marchenko_full/Makefile | 2 +-
marchenko_full/marchenko.c | 4 ++--
raytime/JespersRayTracer.c | 4 ++--
10 files changed, 49 insertions(+), 30 deletions(-)
diff --git a/marchenko_applications/AmpEst.c b/marchenko_applications/AmpEst.c
index d8dbe82..7b04e27 100755
--- a/marchenko_applications/AmpEst.c
+++ b/marchenko_applications/AmpEst.c
@@ -19,7 +19,7 @@ int optncr(int n);
int maxest(float *data, int nt);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
-int AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav)
+void AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav)
{
int l, i, ix, iw, nfreq;
diff --git a/marchenko_applications/Cost.c b/marchenko_applications/Cost.c
index aa15b3a..1c09013 100755
--- a/marchenko_applications/Cost.c
+++ b/marchenko_applications/Cost.c
@@ -15,7 +15,7 @@ typedef struct _complexStruct { /* complex number */
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-int Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn)
+void Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn)
{
int l, i, ix, iw, nfreq;
diff --git a/marchenko_applications/Makefile b/marchenko_applications/Makefile
index 7fbfaca..06b2acd 100644
--- a/marchenko_applications/Makefile
+++ b/marchenko_applications/Makefile
@@ -3,11 +3,11 @@
include ../Make_include
LIBS += -L$L -lgenfft -lm $(LIBSM)
-#OPTC += -g -O0 -Wall
+#OPTC += -g -Wall
#ALL: fmute marchenko marchenko2
-ALL: fmute marchenko_app combine
+ALL: fmute marchenko_app combine gmshift
SRCJ = fmute.c \
getFileInfo.c \
@@ -65,6 +65,16 @@ SRCC = combine.c \
docpkge.c \
readSnapData.c
+SRCG = gmshift.c \
+ getFileInfo.c \
+ writeData.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
+
OBJJ = $(SRCJ:%.c=%.o)
fmute: $(OBJJ)
@@ -80,18 +90,24 @@ OBJC = $(SRCC:%.c=%.o)
combine: $(OBJC)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine $(OBJC) $(LIBS)
-install: fmute marchenko_app combine
+OBJG = $(SRCG:%.c=%.o)
+
+gmshift: $(OBJG)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o gmshift $(OBJG) $(LIBS)
+
+install: fmute marchenko_app combine gmshift
cp fmute $B
cp marchenko_app $B
cp combine $B
+ cp gmshift $B
# cp marchenko2 $B
clean:
- rm -f core fmute $(OBJJ) marchenko_app $(OBJH) combine $(OBJC)
+ rm -f core fmute $(OBJJ) marchenko_app $(OBJH) combine $(OBJC) gmshift $(OBJG)
realclean: clean
- rm -f $B/fmute $B/marchenko_app $B/combine
+ rm -f $B/fmute $B/marchenko_app $B/combine $B/gmshift
diff --git a/marchenko_applications/marchenko.c b/marchenko_applications/marchenko.c
index 7a975bf..b8f0bea 100644
--- a/marchenko_applications/marchenko.c
+++ b/marchenko_applications/marchenko.c
@@ -41,17 +41,17 @@ int readSnapData(char *filename, float *data, segy *hdrs, int nsnaps, int nx, in
int readTinvData(char *filename, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nsyn, float *xsyn, float *zsyn, int iter);
void name_ext(char *filename, char *extension);
-int Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn);
+void Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn);
void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift, int pad, int nt0);
-int AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav);
+void AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav);
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 readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
-int makeWindow(WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
+void makeWindow(WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
void iterations (complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *first, int verbose);
-void imaging (float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int *first, int verbose);
+void imaging (float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose);
void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose);
void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, double *tfft, int *first, int verbose);
@@ -659,6 +659,7 @@ int main (int argc, char **argv)
f2p[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
f1plus[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
green[i*nts+j] = hG_d[ix*nts+nts-j];
+ //vmess("%.8f",f1plus[l*nxs*nts+i*nts+j]);
}
}
}
@@ -675,7 +676,7 @@ int main (int argc, char **argv)
first=0;
imaging(Image,WPs,Refl,nx,nt,nxs,nts,dt,xsyn,Nsyn,xrcv,xsrc,fxs2,fxs,dxs,dxsrc,dx,ixa,ixb,
ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
- f2p,G_d,muteW,smooth,shift,above,pad,nt0,synpos,&first,verbose);
+ f2p,G_d,muteW,smooth,shift,above,pad,nt0,synpos,verbose);
/*============= write output files ================*/
@@ -797,8 +798,17 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
/* reset output data to zero */
memset(&iRN[0], 0, Nsyn*nxs*nts*sizeof(float));
+ idxs = 1.0/dxs;
+ if (ixrcv == NULL) {
+ ixrcv = (int *)malloc(nshots*nx*sizeof(int));
+ }
+ for (k=0; k<nshots; k++) {
+ for (i = 0; i < nx; i++) {
+ ixrcv[k*nx+i] = NINT((xrcv[k*nx+i]-fxs)*idxs);
+ }
+ }
ctrace = (complex *)calloc(ntfft,sizeof(complex));
- if (!first) {
+ if (!*first) {
/* transform muted Ni (Top) to frequency domain, input for next iteration */
for (l = 0; l < Nsyn; l++) {
/* set Fop to zero, so new operator can be defined within ixpossyn points */
@@ -815,7 +825,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
}
else { /* only for first call to synthesis */
/* transform G_d to frequency domain, over all nxs traces */
- first=0;
+ *first=0;
for (l = 0; l < Nsyn; l++) {
/* set Fop to zero, so new operator can be defined within all ix points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
@@ -827,13 +837,6 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
}
}
}
- idxs = 1.0/dxs;
- ixrcv = (int *)malloc(nshots*nx*sizeof(int));
- for (k=0; k<nshots; k++) {
- for (i = 0; i < nx; i++) {
- ixrcv[k*nx+i] = NINT((xrcv[k*nx+i]-fxs)*idxs);
- }
- }
}
free(ctrace);
t1 = wallclock_time();
diff --git a/marchenko_full/AmpEst.c b/marchenko_full/AmpEst.c
index d8dbe82..7b04e27 100755
--- a/marchenko_full/AmpEst.c
+++ b/marchenko_full/AmpEst.c
@@ -19,7 +19,7 @@ int optncr(int n);
int maxest(float *data, int nt);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
-int AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav)
+void AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav)
{
int l, i, ix, iw, nfreq;
diff --git a/marchenko_full/Cost.c b/marchenko_full/Cost.c
index aa15b3a..1c09013 100755
--- a/marchenko_full/Cost.c
+++ b/marchenko_full/Cost.c
@@ -15,7 +15,7 @@ typedef struct _complexStruct { /* complex number */
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-int Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn)
+void Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn)
{
int l, i, ix, iw, nfreq;
diff --git a/marchenko_full/JespersRayTracer.c b/marchenko_full/JespersRayTracer.c
index bd4f9da..ba6ae35 100644
--- a/marchenko_full/JespersRayTracer.c
+++ b/marchenko_full/JespersRayTracer.c
@@ -123,7 +123,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
if (ray.useT2 != 0)
T2 += getdT2(x, z, so, angle, ds, nRayTmp, s, r, rayReference3D, slowness, size);
- if (ray.geomspread != 0) {
+ /*if (ray.geomspread != 0) {
if (so <= 0) {
dQdPhi = 0;
}
@@ -131,7 +131,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
greentmp = greenIntP(lengthRefRay, so, lengthRefRay, slowness, size, nRayTmp, r, s);
dQdPhi += greentmp*secondDerivativeU1(slowness, size, x, z, angle, r, s)*ds/so;
}
- }
+ }*/
}
if (ray.useT2)
diff --git a/marchenko_full/Makefile b/marchenko_full/Makefile
index 931555d..5ffeddc 100644
--- a/marchenko_full/Makefile
+++ b/marchenko_full/Makefile
@@ -3,7 +3,7 @@
include ../Make_include
LIBS += -L$L -lgenfft -lm $(LIBSM)
-OPTC += -g -O0 -Wall
+#OPTC += -g -O0 -Wall
#ALL: fmute marchenko marchenko2
diff --git a/marchenko_full/marchenko.c b/marchenko_full/marchenko.c
index 26bf5dd..82878a1 100644
--- a/marchenko_full/marchenko.c
+++ b/marchenko_full/marchenko.c
@@ -41,9 +41,9 @@ int readSnapData(char *filename, float *data, segy *hdrs, int nsnaps, int nx, in
int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nsyn, float *xsyn, float *zsyn, int iter);
void name_ext(char *filename, char *extension);
-int Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn);
+void Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn);
void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift, int pad, int nt0);
-int AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav);
+void AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav);
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 readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
diff --git a/raytime/JespersRayTracer.c b/raytime/JespersRayTracer.c
index bd4f9da..ba6ae35 100644
--- a/raytime/JespersRayTracer.c
+++ b/raytime/JespersRayTracer.c
@@ -123,7 +123,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
if (ray.useT2 != 0)
T2 += getdT2(x, z, so, angle, ds, nRayTmp, s, r, rayReference3D, slowness, size);
- if (ray.geomspread != 0) {
+ /*if (ray.geomspread != 0) {
if (so <= 0) {
dQdPhi = 0;
}
@@ -131,7 +131,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
greentmp = greenIntP(lengthRefRay, so, lengthRefRay, slowness, size, nRayTmp, r, s);
dQdPhi += greentmp*secondDerivativeU1(slowness, size, x, z, angle, r, s)*ds/so;
}
- }
+ }*/
}
if (ray.useT2)
--
GitLab