Skip to content
Snippets Groups Projects
Commit 2ecb1b4a authored by joeri.brackenhoff's avatar joeri.brackenhoff
Browse files

all

parent f528715e
No related branches found
No related tags found
No related merge requests found
...@@ -87,7 +87,7 @@ OBJJ3 = $(SRCJ3:%.c=%.o) ...@@ -87,7 +87,7 @@ OBJJ3 = $(SRCJ3:%.c=%.o)
fmute3D: $(OBJJ3) fmute3D: $(OBJJ3)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute3D $(OBJJ3) $(LIBS) $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute3D $(OBJJ3) $(LIBS)
install: fmute marchenko test fmute3D install: fmute marchenko marchenko3D fmute3D
cp fmute $B cp fmute $B
cp marchenko $B cp marchenko $B
cp marchenko3D $B cp marchenko3D $B
......
...@@ -7,10 +7,10 @@ export OMP_NUM_THREADS=1 ...@@ -7,10 +7,10 @@ export OMP_NUM_THREADS=1
fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=8 fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=8
#apply the Marchenko algorithm #apply the Marchenko algorithm
~/OpenSource/marchenko3D/test file_shot=shotsy.su file_tinv=p0y.su nshots=901 verbose=10 \ marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
tap=0 niter=8 hw=8 shift=12 smooth=3 \ tap=0 niter=8 hw=8 shift=12 smooth=3 scale=4 \
file_green=pgreen2.su file_gplus=Gplus02.su file_gmin=Gmin02.su \ file_green=pgreen3.su file_gplus=Gplus03.su file_gmin=Gmin03.su \
file_f1plus=f1plus02.su file_f1min=f1min02.su file_f2=f22.su file_f1plus=f1plus03.su file_f1min=f1min03.su file_f2=f23.su
exit exit
......
...@@ -4,9 +4,9 @@ export PATH=$HOME/src/OpenSource/bin:$PATH: ...@@ -4,9 +4,9 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
# Generate the full R matrix for a fixed spread geometry. # Generate the full R matrix for a fixed spread geometry.
dxshot=5000 # with scalco factor of 1000 dxshot=10000 # with scalco factor of 1000
ishot=0 ishot=0
nshots=901 nshots=451
echo $1 echo $1
...@@ -16,16 +16,15 @@ while (( ishot < nshots )) ...@@ -16,16 +16,15 @@ while (( ishot < nshots ))
do do
(( xsrc = -2250000 + ${ishot}*${dxshot} )) (( xsrc = -2250000 + ${ishot}*${dxshot} ))
(( tr1 = 901 - ${ishot} )) (( tr1 = $nshots - ${ishot} ))
(( tr2 = ${tr1} + 900 )) (( tr2 = ${tr1} + $nshots - 1 ))
echo xsrc=$xsrc tr1=$tr1 tr2=$tr2 echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
(( ishot = $ishot + 1)) (( ishot = $ishot + 1))
suwind < shot5_rp.su key=tracl min=$tr1 max=$tr2 | \ suwind < shot_kxky.su key=tracl min=$tr1 max=$tr2 | \
sushw key=sx,gx,fldr,trwf \ sushw key=sx,gx,fldr,trwf \
a=$xsrc,-2250000,$ishot,901 b=0,5000,0,0 j=0,901,0,0 | \ a=$xsrc,-2250000,$ishot,$nshots b=0,$dxshot,0,0 j=0,$nshots,0,0 | \
suchw key1=offset key2=gx key3=sx c=-1 d=1000 >> shotsdx5_rp.su suchw key1=offset key2=gx key3=sx c=-1 d=1000 >> shotsdx5_rp.su
done done
...@@ -24,7 +24,7 @@ void omp_set_num_threads(int num_threads); ...@@ -24,7 +24,7 @@ void omp_set_num_threads(int num_threads);
#ifndef MIN #ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y)) #define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif #endif
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5)) #define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
...@@ -34,36 +34,42 @@ typedef struct _complexStruct { /* complex number */ ...@@ -34,36 +34,42 @@ typedef struct _complexStruct { /* complex number */
} complex; } complex;
#endif/* complex */ #endif/* complex */
int readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots, int nx, int ny, int ntfft, int mode, float scale, int verbose); long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, long *xnx, complex *cdata,
int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, int *xnx, int Nfoc, int nx, int ny, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose); long nw, long nw_low, long nshots, long nx, long ny, long ntfft, long mode, float scale, long verbose);
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
float *zsyn, int *ixpos, int npos, int iter); long *xnx, long Nfoc, long nx, long ny, long ntfft, long mode, long *maxval, float *tinv, long hw, long verbose);
int unique_elements(float *arr, int len); // int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn,
// float *zsyn, int *ixpos, int npos, int iter);
long unique_elements(float *arr, long len);
void name_ext(char *filename, char *extension); void name_ext(char *filename, char *extension);
void applyMute(float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npos, int shift); void applyMute3D(float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nt, long *xrcvsyn, long npos, long shift);
int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, float *sclsxgxsygy, int *nxm); long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
int readData(FILE *fp, float *data, segy *hdrs, int n1); float *sclsxgxsygy, long *nxm);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2); long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
int disp_fileinfo3D(char *file, int n1, int n2, int n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs); long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
double wallclock_time(void); double wallclock_time(void);
void synthesisPositions3D(int nx, int ny, int nxs, int nys, int Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, int *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys, int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose); void synthesisPositions3D(long nx, long ny, long nxs, long nys, long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc,
void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int ny, int nt, int nxs, int nys, int nts, float dt, float *xsyn, float *ysyn, long *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys, long nshots, long nxsrc, long nysrc,
int Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, int *xnx, float fxse, float fxsb, float fyse, float fysb, float dxs, float dys, float dxsrc, long *ixpos, long *npos, long reci, long verbose);
float dysrc, float dx, float dy, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int nxsrc, int nysrc, void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, long nx, long ny, long nt, long nxs, long nys, long nts, float dt,
int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose); float *xsyn, float *ysyn, long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, long *xnx,
float fxse, float fxsb, float fyse, float fysb, float dxs, float dys, float dxsrc, float dysrc,
float dx, float dy, long ntfft, long nw, long nw_low, long nw_high, long mode, long reci, long nshots, long nxsrc, long nysrc,
long *ixpos, long npos, double *tfft, long *isxcount, long *reci_xsrc, long *reci_xrcv, float *ixmask, long verbose);
int linearsearch(int *array, size_t N, int value); long linearsearch(long *array, size_t N, long value);
/*********************** self documentation **********************/ /*********************** self documentation **********************/
char *sdoc[] = { char *sdoc[] = {
" ", " ",
" MARCHENKO - Iterative Green's function and focusing functions retrieval", " MARCHENKO3D - Iterative Green's function and focusing functions retrieval in 3D",
" ", " ",
" marchenko file_tinv= file_shot= [optional parameters]", " marchenko3D file_tinv= file_shot= [optional parameters]",
" ", " ",
" Required parameters: ", " Required parameters: ",
" ", " ",
...@@ -102,7 +108,8 @@ char *sdoc[] = { ...@@ -102,7 +108,8 @@ char *sdoc[] = {
" verbose=0 ................ silent option; >0 displays info", " verbose=0 ................ silent option; >0 displays info",
" ", " ",
" ", " ",
" author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)", " author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)",
" author : Joeri Brackenhoff : 2019 (j.a.brackenhoff@tudelft.nl)",
" ", " ",
NULL}; NULL};
/**************** end self doc ***********************************/ /**************** end self doc ***********************************/
...@@ -111,13 +118,13 @@ int main (int argc, char **argv) ...@@ -111,13 +118,13 @@ int main (int argc, char **argv)
{ {
FILE *fp_out, *fp_f1plus, *fp_f1min; FILE *fp_out, *fp_f1plus, *fp_f1min;
FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin; FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin;
int i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath; long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
int size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad; long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn; long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, n2out, n3out, verbose, ntfft; long reci, countmin, mode, n2out, n3out, verbose, ntfft;
int iter, niter, tracf, *muteW; long iter, niter, tracf, *muteW;
int hw, smooth, above, shift, *ixpos, npos, ix; long hw, smooth, above, shift, *ixpos, npos, ix;
int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv; long nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn; float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0; double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
float d1, d2, d3, f1, f2, f3, fxsb, fxse, fysb, fyse, ft, fx, fy, *xsyn, *ysyn, dxsrc, dysrc; float d1, d2, d3, f1, f2, f3, fxsb, fxse, fysb, fyse, ft, fx, fy, *xsyn, *ysyn, dxsrc, dysrc;
...@@ -146,7 +153,7 @@ int main (int argc, char **argv) ...@@ -146,7 +153,7 @@ int main (int argc, char **argv)
if (!getparstring("file_f2", &file_f2)) file_f2 = NULL; if (!getparstring("file_f2", &file_f2)) file_f2 = NULL;
if (!getparstring("file_pmin", &file_pmin)) file_pmin = NULL; if (!getparstring("file_pmin", &file_pmin)) file_pmin = NULL;
if (!getparstring("file_iter", &file_iter)) file_iter = NULL; if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
if (!getparint("verbose", &verbose)) verbose = 0; if (!getparlong("verbose", &verbose)) verbose = 0;
if (file_tinv == NULL && file_shot == NULL) if (file_tinv == NULL && file_shot == NULL)
verr("file_tinv and file_shot cannot be both input pipe"); verr("file_tinv and file_shot cannot be both input pipe");
if (!getparstring("file_green", &file_green)) { if (!getparstring("file_green", &file_green)) {
...@@ -155,20 +162,20 @@ int main (int argc, char **argv) ...@@ -155,20 +162,20 @@ int main (int argc, char **argv)
} }
if (!getparfloat("fmin", &fmin)) fmin = 0.0; if (!getparfloat("fmin", &fmin)) fmin = 0.0;
if (!getparfloat("fmax", &fmax)) fmax = 70.0; if (!getparfloat("fmax", &fmax)) fmax = 70.0;
if (!getparint("reci", &reci)) reci = 0; if (!getparlong("reci", &reci)) reci = 0;
if (!getparfloat("scale", &scale)) scale = 2.0; if (!getparfloat("scale", &scale)) scale = 2.0;
if (!getparfloat("tsq", &tsq)) tsq = 0.0; if (!getparfloat("tsq", &tsq)) tsq = 0.0;
if (!getparfloat("Q", &Q)) Q = 0.0; if (!getparfloat("Q", &Q)) Q = 0.0;
if (!getparfloat("f0", &f0)) f0 = 0.0; if (!getparfloat("f0", &f0)) f0 = 0.0;
if (!getparint("tap", &tap)) tap = 0; if (!getparlong("tap", &tap)) tap = 0;
if (!getparint("ntap", &ntap)) ntap = 0; if (!getparlong("ntap", &ntap)) ntap = 0;
if (!getparint("pad", &pad)) pad = 0; if (!getparlong("pad", &pad)) pad = 0;
if(!getparint("niter", &niter)) niter = 10; if(!getparlong("niter", &niter)) niter = 10;
if(!getparint("hw", &hw)) hw = 15; if(!getparlong("hw", &hw)) hw = 15;
if(!getparint("smooth", &smooth)) smooth = 5; if(!getparlong("smooth", &smooth)) smooth = 5;
if(!getparint("above", &above)) above = 0; if(!getparlong("above", &above)) above = 0;
if(!getparint("shift", &shift)) shift=12; if(!getparlong("shift", &shift)) shift=12;
if (reci && ntap) vwarn("tapering influences the reciprocal result"); if (reci && ntap) vwarn("tapering influences the reciprocal result");
...@@ -194,12 +201,12 @@ int main (int argc, char **argv) ...@@ -194,12 +201,12 @@ int main (int argc, char **argv)
ntfft = optncr(MAX(nt+pad, nts+pad)); ntfft = optncr(MAX(nt+pad, nts+pad));
nfreq = ntfft/2+1; nfreq = ntfft/2+1;
nw_low = (int)MIN((fmin*ntfft*dt), nfreq-1); nw_low = (long)MIN((fmin*ntfft*dt), nfreq-1);
nw_low = MAX(nw_low, 1); nw_low = MAX(nw_low, 1);
nw_high = MIN((int)(fmax*ntfft*dt), nfreq-1); nw_high = MIN((long)(fmax*ntfft*dt), nfreq-1);
nw = nw_high - nw_low + 1; nw = nw_high - nw_low + 1;
scl = 1.0/((float)ntfft); scl = 1.0/((float)ntfft);
if (!getparint("countmin", &countmin)) countmin = 0.3*nx*ny; if (!getparlong("countmin", &countmin)) countmin = 0.3*nx*ny;
/*================ Allocating all data arrays ================*/ /*================ Allocating all data arrays ================*/
...@@ -212,7 +219,7 @@ int main (int argc, char **argv) ...@@ -212,7 +219,7 @@ int main (int argc, char **argv)
iRN = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float)); iRN = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
Ni = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float)); Ni = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
G_d = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float)); G_d = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
muteW = (int *)calloc(Nfoc*nys*nxs,sizeof(int)); muteW = (long *)calloc(Nfoc*nys*nxs,sizeof(long));
trace = (float *)malloc(ntfft*sizeof(float)); trace = (float *)malloc(ntfft*sizeof(float));
tapersy = (float *)malloc(nxs*sizeof(float)); tapersy = (float *)malloc(nxs*sizeof(float));
xrcvsyn = (float *)calloc(Nfoc*nys*nxs,sizeof(float)); // x-rcv postions of focal points xrcvsyn = (float *)calloc(Nfoc*nys*nxs,sizeof(float)); // x-rcv postions of focal points
...@@ -220,8 +227,8 @@ int main (int argc, char **argv) ...@@ -220,8 +227,8 @@ int main (int argc, char **argv)
xsyn = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points xsyn = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points
ysyn = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points ysyn = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points
zsyn = (float *)malloc(Nfoc*sizeof(float)); // z-src position of focal points zsyn = (float *)malloc(Nfoc*sizeof(float)); // z-src position of focal points
xnxsyn = (int *)calloc(Nfoc,sizeof(int)); // number of traces per focal point xnxsyn = (long *)calloc(Nfoc,sizeof(long)); // number of traces per focal point
ixpos = (int *)calloc(nys*nxs,sizeof(int)); // x-position of source of shot in G_d domain (nxs*nys with dxs, dys) ixpos = (long *)calloc(nys*nxs,sizeof(long)); // x-position of source of shot in G_d domain (nxs*nys with dxs, dys)
Refl = (complex *)malloc(nw*ny*nx*nshots*sizeof(complex)); Refl = (complex *)malloc(nw*ny*nx*nshots*sizeof(complex));
tapersh = (float *)malloc(nx*sizeof(float)); tapersh = (float *)malloc(nx*sizeof(float));
...@@ -230,12 +237,12 @@ int main (int argc, char **argv) ...@@ -230,12 +237,12 @@ int main (int argc, char **argv)
xsrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots xsrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
ysrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots ysrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
zsrc = (float *)calloc(nshots,sizeof(float)); // z-src position of shots zsrc = (float *)calloc(nshots,sizeof(float)); // z-src position of shots
xnx = (int *)calloc(nshots,sizeof(int)); // number of traces per shot xnx = (long *)calloc(nshots,sizeof(long)); // number of traces per shot
if (reci!=0) { if (reci!=0) {
reci_xsrc = (int *)malloc((nxs*nxs*nys*nys)*sizeof(int)); reci_xsrc = (long *)malloc((nxs*nxs*nys*nys)*sizeof(long));
reci_xrcv = (int *)malloc((nxs*nxs*nys*nys)*sizeof(int)); reci_xrcv = (long *)malloc((nxs*nxs*nys*nys)*sizeof(long));
isxcount = (int *)calloc(nxs*nys,sizeof(int)); isxcount = (long *)calloc(nxs*nys,sizeof(long));
ixmask = (float *)calloc(nxs*nys,sizeof(float)); ixmask = (float *)calloc(nxs*nys,sizeof(float));
} }
...@@ -518,9 +525,9 @@ int main (int argc, char **argv) ...@@ -518,9 +525,9 @@ int main (int argc, char **argv)
t3 = wallclock_time(); t3 = wallclock_time();
tsyn += t3 - t2; tsyn += t3 - t2;
if (file_iter != NULL) { // if (file_iter != NULL) {
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs*nys, d2, f2, n2out*n3out, Nfoc, xsyn, zsyn, ixpos, npos, iter); // writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs*nys, d2, f2, n2out*n3out, Nfoc, xsyn, zsyn, ixpos, npos, iter);
} // }
/* N_k(x,t) = -N_(k-1)(x,-t) */ /* N_k(x,t) = -N_(k-1)(x,-t) */
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */ /* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
for (l = 0; l < Nfoc; l++) { for (l = 0; l < Nfoc; l++) {
...@@ -604,6 +611,7 @@ int main (int argc, char **argv) ...@@ -604,6 +611,7 @@ int main (int argc, char **argv)
} }
} }
} }
applyMute(green, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
/* compute upgoing Green's function G^+,- */ /* compute upgoing Green's function G^+,- */
if (file_gmin != NULL) { if (file_gmin != NULL) {
...@@ -629,7 +637,7 @@ int main (int argc, char **argv) ...@@ -629,7 +637,7 @@ int main (int argc, char **argv)
} }
} }
/* Apply mute with window for Gmin */ /* Apply mute with window for Gmin */
applyMute(Gmin, muteW, smooth, 1, Nfoc, nxs*nys, nts, ixpos, npos, shift); applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
} /* end if Gmin */ } /* end if Gmin */
/* compute downgoing Green's function G^+,+ */ /* compute downgoing Green's function G^+,+ */
...@@ -655,6 +663,8 @@ int main (int argc, char **argv) ...@@ -655,6 +663,8 @@ int main (int argc, char **argv)
} }
} }
} }
/* Apply mute with window for Gplus */
applyMute(Gplus, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
} /* end if Gplus */ } /* end if Gplus */
t2 = wallclock_time(); t2 = wallclock_time();
......
...@@ -122,25 +122,28 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float ...@@ -122,25 +122,28 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
xnx[isyn]=itrace; xnx[isyn]=itrace;
/* alternative find maximum at source position */ /* alternative find maximum at source position */
dxrcv = (gx1 - gx0)*scl/(float)(nx1-1); if (nx1>1) dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
dyrcv = (gy1 - gy0)*scl/(float)(ny1-1); else dxrcv = (gx1 - gx0)*scl/(float)(1);
//imax = NINT(((sx_shot-gx0)*scl)/dxrcv); if (dxrcv==0.0) dxrcv=1.0;
ixmax = NINT(((sx_shot-gx0)*scl)/dxrcv); ixmax = NINT(((sx_shot-gx0)*scl)/dxrcv);
if (ny1>1) dyrcv = (gy1 - gy0)*scl/(float)(ny1-1);
else dyrcv = (gy1 - gy0)*scl/(float)(1);
if (dyrcv==0.0) dyrcv=1.0;
iymax = NINT(((sy_shot-gy0)*scl)/dyrcv); iymax = NINT(((sy_shot-gy0)*scl)/dyrcv);
if (iymax > ny1-1) { if (iymax > ny1-1) {
vmess("source of y is past array, snapping to nearest y"); vmess("source of y (%d) is past array, snapping to nearest y (%d)",iymax,ny1-1);
iymax = ny1-1; iymax = ny1-1;
} }
if (iymax < 0) { if (iymax < 0) {
vmess("source of y is before array, snapping to nearest y"); vmess("source of y (%d) is before array, snapping to nearest y (%d)",iymax,0);
iymax = 0; iymax = 0;
} }
if (ixmax > nx1-1) { if (ixmax > nx1-1) {
vmess("source of x is past array, snapping to nearest x"); vmess("source of x (%d) is past array, snapping to nearest x (%d)",ixmax,nx1-1);
ixmax = nx1-1; ixmax = nx1-1;
} }
if (ixmax < 0) { if (ixmax < 0) {
vmess("source of x is before array, snapping to nearest x"); vmess("source of x (%d) is before array, snapping to nearest x (%d)",ixmax,nx1-1);
ixmax = 0; ixmax = 0;
} }
tmax=0.0; tmax=0.0;
......
...@@ -136,12 +136,21 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr ...@@ -136,12 +136,21 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
int nfreq, size, inx; int nfreq, size, inx;
float scl; float scl;
int i, j, l, m, iw, ix, k, isrc, il, ik, nxy, nxys; int i, j, l, m, iw, ix, k, isrc, il, ik, nxy, nxys;
float *rtrace, idxs, idys; float *rtrace, idxs, idys, fxb, fyb, fxe, fye;
complex *sum, *ctrace; complex *sum, *ctrace;
int npe; int npe;
static int first=1, *ircv; static int first=1, *ircv;
static double t0, t1, t; static double t0, t1, t;
if (fxsb < 0) fxb = 1.001*fxsb;
else fxb = 0.999*fxsb;
if (fysb < 0) fyb = 1.001*fysb;
else fyb = 0.999*fysb;
if (fxse > 0) fxe = 1.001*fxse;
else fxe = 0.999*fxse;
if (fyse > 0) fye = 1.001*fyse;
else fye = 0.999*fyse;
nxy = nx*ny; nxy = nx*ny;
nxys = nxs*nys; nxys = nxs*nys;
...@@ -149,7 +158,7 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr ...@@ -149,7 +158,7 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
nfreq = ntfft/2+1; nfreq = ntfft/2+1;
/* scale factor 1/N for backward FFT, /* scale factor 1/N for backward FFT,
* scale dt for correlation/convolution along time, * scale dt for correlation/convolution along time,
* scale dx (or dxsrc) for integration over receiver (or shot) coordinates */ * scale dx*dy (or dxsrc*dysrc) for integration over receiver (or shot) coordinates */
scl = 1.0*dt/((float)ntfft); scl = 1.0*dt/((float)ntfft);
#ifdef _OPENMP #ifdef _OPENMP
...@@ -212,7 +221,7 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr ...@@ -212,7 +221,7 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
/* Loop over total number of shots */ /* Loop over total number of shots */
if (reci == 0 || reci == 1) { if (reci == 0 || reci == 1) {
for (k=0; k<nshots; k++) { for (k=0; k<nshots; k++) {
if ((xsrc[k] < 0.999*fxsb) || (xsrc[k] > 1.001*fxse) || (ysrc[k] < 0.999*fysb) || (ysrc[k] > 1.001*fyse)) continue; if ((xsrc[k] < fxb) || (xsrc[k] > fxe) || (ysrc[k] < fyb) || (ysrc[k] > fye)) continue;
isrc = NINT((ysrc[k] - fysb)/dys)*nxs+NINT((xsrc[k] - fxsb)/dxs); isrc = NINT((ysrc[k] - fysb)/dys)*nxs+NINT((xsrc[k] - fxsb)/dxs);
inx = xnx[k]; /* number of traces per shot */ inx = xnx[k]; /* number of traces per shot */
...@@ -234,8 +243,8 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr ...@@ -234,8 +243,8 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
/* compute integral over receiver positions */ /* compute integral over receiver positions */
/* multiply R with Fop and sum over nx */ /* multiply R with Fop and sum over nx */
memset(&sum[0].r,0,nfreq*2*sizeof(float)); memset(&sum[0].r,0,nfreq*2*sizeof(float));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) { for (i = 0; i < inx; i++) {
for (i = 0; i < inx; i++) { for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
ix = ircv[k*nxy+i]; ix = ircv[k*nxy+i];
sum[j].r += Refl[k*nw*nxy+m*nxy+i].r*Fop[l*nw*nxys+m*nxys+ix].r - sum[j].r += Refl[k*nw*nxy+m*nxy+i].r*Fop[l*nw*nxys+m*nxys+ix].r -
Refl[k*nw*nxy+m*nxy+i].i*Fop[l*nw*nxys+m*nxys+ix].i; Refl[k*nw*nxy+m*nxy+i].i*Fop[l*nw*nxys+m*nxys+ix].i;
......
...@@ -6,7 +6,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM) ...@@ -6,7 +6,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC += -openmp #OPTC += -openmp
#OPTC += -g -O0 #OPTC += -g -O0
ALL: makemod makewave extendModel fconv correigen green basop syn2d mat2su ftr1d ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d
SRCM = \ SRCM = \
makemod.c \ makemod.c \
...@@ -84,6 +84,17 @@ SRCG = green.c \ ...@@ -84,6 +84,17 @@ SRCG = green.c \
docpkge.c \ docpkge.c \
getpars.c getpars.c
SRCG3 = green3D.c \
getFileInfo.c \
getrecpos3D.c \
readData.c \
writeData.c \
wallclock_time.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
getpars.c
SRCB = basop.c \ SRCB = basop.c \
getFileInfo.c \ getFileInfo.c \
kxwfilter.c \ kxwfilter.c \
...@@ -152,6 +163,11 @@ OBJG = $(SRCG:%.c=%.o) ...@@ -152,6 +163,11 @@ OBJG = $(SRCG:%.c=%.o)
green: $(OBJG) green: $(OBJG)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o green $(OBJG) $(LIBS) $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o green $(OBJG) $(LIBS)
OBJG3 = $(SRCG3:%.c=%.o)
green3D: $(OBJG3)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o green3D $(OBJG3) $(LIBS)
OBJB = $(SRCB:%.c=%.o) OBJB = $(SRCB:%.c=%.o)
basop: $(OBJB) basop: $(OBJB)
...@@ -172,23 +188,24 @@ OBJT = $(SRCT:%.c=%.o) ...@@ -172,23 +188,24 @@ OBJT = $(SRCT:%.c=%.o)
ftr1d: $(OBJT) ftr1d: $(OBJT)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o ftr1d $(OBJT) $(LIBS) $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o ftr1d $(OBJT) $(LIBS)
install: makemod makewave extendModel fconv correigen green basop syn2d mat2su ftr1d install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d
cp makemod $B cp makemod $B
cp makewave $B cp makewave $B
cp extendModel $B cp extendModel $B
cp fconv $B cp fconv $B
cp correigen $B cp correigen $B
cp green $B cp green $B
cp green3D $B
cp basop $B cp basop $B
cp syn2d $B cp syn2d $B
cp mat2su $B cp mat2su $B
cp ftr1d $B cp ftr1d $B
clean: clean:
rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT)
realclean: clean realclean: clean
rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/basop $B/syn2d $B/mat2su $B/ftr1d rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d
......
#include "par.h"
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
/**
* read receiver positions used in green
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define SGN(x) ((x) < 0 ? -1.0 : 1.0)
#ifndef ABS
#define ABS(x) ((x) < 0 ? -(x) : (x))
#endif
void getrecpos3D(float *xi, float *yi, float *zi, int nx, int ny, float *xrcv, float *yrcv, float *zrcv, int verbose)
{
int nrx, nry, i, j, l, ndeltx, ndelty, np, lint, seed;
long idum;
float xprev, yprev, zprev, deltx, delty, deltz, dxrcv, dyrcv, dzrcv, var, irr, maxirr;
float rrcv, dphi, oxrcv, oyrcv, ozrcv;
nrx = countparval("xrcv");
nry = countparval("yrcv");
if(!getparfloat("dxrcv",&dxrcv)) dxrcv = 15;
if(!getparfloat("dyrcv",&dyrcv)) dyrcv = 15;
if(!getparfloat("var", &var)) var=0;
if(!getparint("lint", &lint)) lint=1;
if(!getparint("seed", &seed)) seed=0;
/* check if receiver positions on a circle are defined */
if (getparfloat("rrcv", &rrcv)) {
if (!getparfloat("dphi",&dphi)) dphi=2.0;
if (!getparfloat("oxrcv",&oxrcv)) oxrcv=0.0;
if (!getparfloat("oyrcv",&oyrcv)) oyrcv=0.0;
if (!getparfloat("ozrcv",&ozrcv)) ozrcv=0.0;
np = 0;
for (i=0; i<ny; i++) {
for (j=0; j<ny; j++) {
xi[np] = oxrcv+rrcv*cos(((i*dphi)/360.0)*(2.0*M_PI));
yi[np] = oyrcv+rrcv*cos(((i*dphi)/360.0)*(2.0*M_PI));
zi[np++] = ozrcv+rrcv*sin(((i*dphi)/360.0)*(2.0*M_PI));
if (verbose>4) fprintf(stderr,"Receiver Circle: xrcv[%d]=%f yrcv=%f zrcv=%f\n", i, xi[i], yi[i], zi[i]);
}
}
return;
}
if (var <= 0) {
if (lint == 1) {
xprev = xrcv[0];
yprev = yrcv[0];
zprev = zrcv[0];
np = 0;
for (i = 1; i < nry; i++) {
for (l = 1; l < nrx; l++) {
deltx = xrcv[i] - xprev;
delty = yrcv[i] - yprev;
deltz = zrcv[i] - zprev;
ndeltx = NINT(ABS(deltx/dxrcv));
ndelty = NINT(ABS(delty/dyrcv));
dzrcv = deltz/ndeltx;
for (j = 0; j < ndeltx; j++) {
zi[np] = zprev + j*dzrcv;
yi[np] = yprev + i*dyrcv;
xi[np++] = xprev + j*dxrcv;
}
xprev = xrcv[i*nx+l];
yprev = yrcv[i*nx+l];
zprev = zrcv[i*nx+l];
}
xi[i*nx+nx-1] = xrcv[nrx-1];
yi[i*nx+nx-1] = yrcv[nrx-1];
zi[i*nx+nx-1] = zrcv[nrx-1];
}
}
else {
for (i = 0; i < nry; i++) {
for (l = 0; l < nrx; l++) {
xi[i*nx+l] = xrcv[l];
yi[i*nx+l] = yrcv[i];
zi[i*nx+l] = zrcv[l];
}
}
}
}
else {
xprev = xrcv[0];
yprev = yrcv[0];
zprev = zrcv[0];
np = 0;
maxirr = 0;
idum = (long) seed;
srand48(idum);
for (i = 1; i < nrx; i++) {
deltx = xrcv[i] - xprev;
deltz = zrcv[i] - zprev;
ndeltx = NINT(ABS(deltx/dxrcv));
dzrcv = deltz/ndeltx;
for (j = 0; j < ndeltx; j++) {
irr = var*((float)drand48());
if (fabs(irr) > maxirr) maxirr = fabs(irr);
zi[np] = zprev + j*dzrcv;
xi[np++] = xprev + j*dxrcv + irr;
if (verbose==13)vmess("xrcv %d = %f (%f)",np-1,xi[np-1], irr);
}
xprev = xrcv[i];
zprev = zrcv[i];
}
irr = var*((float)drand48());
if (fabs(irr) > maxirr) maxirr = fabs(irr);
xi[nx-1] = xrcv[nrx-1] + irr;
zi[nx-1] = zrcv[nrx-1];
if (verbose) vmess("maximum error in receiver position %f", maxirr);
if (verbose==13) vmess("xrcv %d = %f (%f)", nx-1, xi[nx-1], irr);
}
if (verbose) vmess("getrecpos number of receivers = %d", np+1);
return;
}
File added
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment