diff --git a/marchenko/fmute.c b/marchenko/fmute.c
index 4d63faf38fa8bb4cb489ee554bca07e52ce9a439..f5c7668a16dc382b586c2981427888c6fa256f66 100644
--- a/marchenko/fmute.c
+++ b/marchenko/fmute.c
@@ -228,6 +228,8 @@ int main (int argc, char **argv)
/* alternative find maximum at source position */
dxrcv = (hdrs_in1[nx1-1].gx - hdrs_in1[0].gx)*sclsxgx/(float)(nx1-1);
imax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
+ /* make sure that the position fits into the receiver array */
+ imax = MIN(MAX(0,imax),nx1-1);
tmax=0.0;
jmax = 0;
xmax=0.0;
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 491a1c0b10c7b333768fb81720759f0a0fa71b13..6a8acbd37b6d7fbef7ae55138ff5e673656ca4eb 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -87,7 +87,7 @@ OBJJ3 = $(SRCJ3:%.c=%.o)
fmute3D: $(OBJJ3)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute3D $(OBJJ3) $(LIBS)
-install: fmute marchenko test fmute3D
+install: fmute marchenko marchenko3D fmute3D
cp fmute $B
cp marchenko $B
cp marchenko3D $B
diff --git a/marchenko3D/demo/oneD/marchenko.scr b/marchenko3D/demo/oneD/marchenko.scr
index 6fd6cacab5f12eda0d2ecafb39c5ac9d78313f46..8538681e22d6a73dbaef67b3c09efeebfbce2916 100755
--- a/marchenko3D/demo/oneD/marchenko.scr
+++ b/marchenko3D/demo/oneD/marchenko.scr
@@ -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
#apply the Marchenko algorithm
-~/OpenSource/marchenko3D/test file_shot=shotsy.su file_tinv=p0y.su nshots=901 verbose=10 \
- tap=0 niter=8 hw=8 shift=12 smooth=3 \
- file_green=pgreen2.su file_gplus=Gplus02.su file_gmin=Gmin02.su \
- file_f1plus=f1plus02.su file_f1min=f1min02.su file_f2=f22.su
+marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
+ tap=0 niter=8 hw=8 shift=12 smooth=3 scale=4 \
+ file_green=pgreen3.su file_gplus=Gplus03.su file_gmin=Gmin03.su \
+ file_f1plus=f1plus03.su file_f1min=f1min03.su file_f2=f23.su
exit
diff --git a/marchenko3D/demo/oneD/p5all.scr b/marchenko3D/demo/oneD/p5all.scr
index 333be5510ec6a203c098595abfdabe5cdba2466b..c749523043cbeb7f11ea3a1a86f72b397b805271 100755
--- a/marchenko3D/demo/oneD/p5all.scr
+++ b/marchenko3D/demo/oneD/p5all.scr
@@ -4,9 +4,9 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
# 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
-nshots=901
+nshots=451
echo $1
@@ -16,16 +16,15 @@ while (( ishot < nshots ))
do
(( xsrc = -2250000 + ${ishot}*${dxshot} ))
- (( tr1 = 901 - ${ishot} ))
- (( tr2 = ${tr1} + 900 ))
+ (( tr1 = $nshots - ${ishot} ))
+ (( tr2 = ${tr1} + $nshots - 1 ))
echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
(( 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 \
- 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
done
-
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 7323516501465457ae78f8bea999ba224fdd5c92..983aa15579a63adbc16529cf527ce878c1068e57 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -24,7 +24,7 @@ void omp_set_num_threads(int num_threads);
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#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 */
} 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);
-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);
-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);
-int unique_elements(float *arr, int len);
+long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, long *xnx, complex *cdata,
+ long nw, long nw_low, long nshots, long nx, long ny, long ntfft, long mode, float scale, long verbose);
+long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
+ long *xnx, long Nfoc, long nx, long ny, long ntfft, long mode, long *maxval, float *tinv, long hw, long verbose);
+// 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 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);
-int readData(FILE *fp, float *data, segy *hdrs, int n1);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-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 getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
+ float *sclsxgxsygy, long *nxm);
+long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
+long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
+long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
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 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,
-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,
-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,
-int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
+void synthesisPositions3D(long nx, long ny, long nxs, long nys, long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc,
+ long *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys, long nshots, long nxsrc, long nysrc,
+ long *ixpos, long *npos, long reci, long verbose);
+void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, long nx, long ny, long nt, long nxs, long nys, long nts, float dt,
+ 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 **********************/
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: ",
" ",
@@ -102,7 +108,8 @@ char *sdoc[] = {
" 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};
/**************** end self doc ***********************************/
@@ -111,13 +118,13 @@ int main (int argc, char **argv)
{
FILE *fp_out, *fp_f1plus, *fp_f1min;
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;
- int size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
- int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
- int reci, countmin, mode, n2out, n3out, verbose, ntfft;
- int iter, niter, tracf, *muteW;
- int hw, smooth, above, shift, *ixpos, npos, ix;
- int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
+ long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
+ long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
+ long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
+ long reci, countmin, mode, n2out, n3out, verbose, ntfft;
+ long iter, niter, tracf, *muteW;
+ long hw, smooth, above, shift, *ixpos, npos, ix;
+ long nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
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)
if (!getparstring("file_f2", &file_f2)) file_f2 = NULL;
if (!getparstring("file_pmin", &file_pmin)) file_pmin = 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)
verr("file_tinv and file_shot cannot be both input pipe");
if (!getparstring("file_green", &file_green)) {
@@ -155,20 +162,20 @@ int main (int argc, char **argv)
}
if (!getparfloat("fmin", &fmin)) fmin = 0.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("tsq", &tsq)) tsq = 0.0;
if (!getparfloat("Q", &Q)) Q = 0.0;
if (!getparfloat("f0", &f0)) f0 = 0.0;
- if (!getparint("tap", &tap)) tap = 0;
- if (!getparint("ntap", &ntap)) ntap = 0;
- if (!getparint("pad", &pad)) pad = 0;
+ if (!getparlong("tap", &tap)) tap = 0;
+ if (!getparlong("ntap", &ntap)) ntap = 0;
+ if (!getparlong("pad", &pad)) pad = 0;
- if(!getparint("niter", &niter)) niter = 10;
- if(!getparint("hw", &hw)) hw = 15;
- if(!getparint("smooth", &smooth)) smooth = 5;
- if(!getparint("above", &above)) above = 0;
- if(!getparint("shift", &shift)) shift=12;
+ if(!getparlong("niter", &niter)) niter = 10;
+ if(!getparlong("hw", &hw)) hw = 15;
+ if(!getparlong("smooth", &smooth)) smooth = 5;
+ if(!getparlong("above", &above)) above = 0;
+ if(!getparlong("shift", &shift)) shift=12;
if (reci && ntap) vwarn("tapering influences the reciprocal result");
@@ -194,12 +201,12 @@ int main (int argc, char **argv)
ntfft = optncr(MAX(nt+pad, nts+pad));
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_high = MIN((int)(fmax*ntfft*dt), nfreq-1);
+ nw_high = MIN((long)(fmax*ntfft*dt), nfreq-1);
nw = nw_high - nw_low + 1;
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 ================*/
@@ -212,7 +219,7 @@ int main (int argc, char **argv)
iRN = (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));
- muteW = (int *)calloc(Nfoc*nys*nxs,sizeof(int));
+ muteW = (long *)calloc(Nfoc*nys*nxs,sizeof(long));
trace = (float *)malloc(ntfft*sizeof(float));
tapersy = (float *)malloc(nxs*sizeof(float));
xrcvsyn = (float *)calloc(Nfoc*nys*nxs,sizeof(float)); // x-rcv postions of focal points
@@ -220,8 +227,8 @@ int main (int argc, char **argv)
xsyn = (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
- xnxsyn = (int *)calloc(Nfoc,sizeof(int)); // 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)
+ xnxsyn = (long *)calloc(Nfoc,sizeof(long)); // number of traces per focal point
+ ixpos = (long *)calloc(nys*nxs,sizeof(long)); // x-position of source of shot in G_d domain (nxs*nys with dxs, dys)
Refl = (complex *)malloc(nw*ny*nx*nshots*sizeof(complex));
tapersh = (float *)malloc(nx*sizeof(float));
@@ -230,12 +237,12 @@ int main (int argc, char **argv)
xsrc = (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
- 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) {
- reci_xsrc = (int *)malloc((nxs*nxs*nys*nys)*sizeof(int));
- reci_xrcv = (int *)malloc((nxs*nxs*nys*nys)*sizeof(int));
- isxcount = (int *)calloc(nxs*nys,sizeof(int));
+ reci_xsrc = (long *)malloc((nxs*nxs*nys*nys)*sizeof(long));
+ reci_xrcv = (long *)malloc((nxs*nxs*nys*nys)*sizeof(long));
+ isxcount = (long *)calloc(nxs*nys,sizeof(long));
ixmask = (float *)calloc(nxs*nys,sizeof(float));
}
@@ -518,9 +525,9 @@ int main (int argc, char **argv)
t3 = wallclock_time();
tsyn += t3 - t2;
- if (file_iter != NULL) {
- writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs*nys, d2, f2, n2out*n3out, Nfoc, xsyn, zsyn, ixpos, npos, iter);
- }
+ // if (file_iter != NULL) {
+ // 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) */
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
for (l = 0; l < Nfoc; l++) {
@@ -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^+,- */
if (file_gmin != NULL) {
@@ -629,7 +637,7 @@ int main (int argc, char **argv)
}
}
/* 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 */
/* compute downgoing Green's function G^+,+ */
@@ -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 */
t2 = wallclock_time();
diff --git a/marchenko3D/readTinvData3D.c b/marchenko3D/readTinvData3D.c
index c6a9e45811c6d24b25dc65ffca7d257444a46e06..845d08076656e6987299f3fdbbbcf5d0aa55ec1f 100644
--- a/marchenko3D/readTinvData3D.c
+++ b/marchenko3D/readTinvData3D.c
@@ -122,25 +122,28 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
xnx[isyn]=itrace;
/* alternative find maximum at source position */
- dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
- dyrcv = (gy1 - gy0)*scl/(float)(ny1-1);
- //imax = NINT(((sx_shot-gx0)*scl)/dxrcv);
+ if (nx1>1) dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
+ else dxrcv = (gx1 - gx0)*scl/(float)(1);
+ if (dxrcv==0.0) dxrcv=1.0;
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);
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;
}
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;
}
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;
}
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;
}
tmax=0.0;
diff --git a/marchenko3D/synthesis3D.c b/marchenko3D/synthesis3D.c
index 8ab12a94806fe3663dc8d4f27d6f3f388d71afd2..19d4c80b2ade811d3afb0b139a35d68cb34aafc5 100644
--- a/marchenko3D/synthesis3D.c
+++ b/marchenko3D/synthesis3D.c
@@ -136,12 +136,21 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
int nfreq, size, inx;
float scl;
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;
int npe;
static int first=1, *ircv;
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;
nxys = nxs*nys;
@@ -149,7 +158,7 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
nfreq = ntfft/2+1;
/* scale factor 1/N for backward FFT,
* 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);
#ifdef _OPENMP
@@ -212,7 +221,7 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
/* Loop over total number of shots */
if (reci == 0 || reci == 1) {
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);
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
/* compute integral over receiver positions */
/* multiply R with Fop and sum over nx */
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];
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;
diff --git a/raytime/Makefile b/raytime/Makefile
index 192a17315ed0dc71f494646357315a4ffcc2dc82..b3aa6f39b0b3a26f00e0c892e352eebbbe223768 100644
--- a/raytime/Makefile
+++ b/raytime/Makefile
@@ -47,16 +47,16 @@ SRCC = $(PRG).c \
OBJC = $(SRCC:%.c=%.o)
$(PRG): $(OBJC) raytime.h
- $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o Raytime $(OBJC) $(LIBS)
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o raytime $(OBJC) $(LIBS)
install: raytime
- cp Raytime $B
+ cp raytime $B
clean:
- rm -f core $(OBJC) $(OBJM) Raytime
+ rm -f core $(OBJC) $(OBJM) raytime
realclean:
- rm -f core $(OBJC) $(OBJM) $(PRG) $B/Raytime
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/raytime
print: Makefile $(SRC)
diff --git a/raytime3d/Makefile b/raytime3d/Makefile
index 1bed95dcdef7cd3fb664ca9017a3bad0b38788b3..be0bb6e043b09ee22c89254f9c54afceb990ad92 100644
--- a/raytime3d/Makefile
+++ b/raytime3d/Makefile
@@ -27,15 +27,15 @@ PRG = raytime3d
SRCC = $(PRG).c \
vidale3d.c \
src3d.c \
- getParameters.c \
+ getParameters3d.c \
getWaveletInfo.c \
writeSrcRecPos.c \
- readModel.c \
+ readModel3d.c \
getWaveletHeaders.c \
verbosepkg.c \
- getModelInfo.c \
+ getModelInfo3d.c \
wallclock_time.c \
- recvPar.c \
+ recvPar3d.c \
writesufile.c \
name_ext.c \
atopkge.c \
@@ -45,17 +45,17 @@ SRCC = $(PRG).c \
OBJC = $(SRCC:%.c=%.o)
-$(PRG): $(OBJC) raytime.h
- $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o Raytime $(OBJC) $(LIBS)
+$(PRG): $(OBJC) raytime3d.h
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o raytime3d $(OBJC) $(LIBS)
-install: raytime
- cp Raytime $B
+install: raytime3d
+ cp raytime3d $B
clean:
- rm -f core $(OBJC) $(OBJM) Raytime
+ rm -f core $(OBJC) $(OBJM) raytime3d
realclean:
- rm -f core $(OBJC) $(OBJM) $(PRG) $B/Raytime
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/raytime3d
print: Makefile $(SRC)
diff --git a/raytime3d/getModelInfo.c b/raytime3d/getModelInfo.c
deleted file mode 100644
index 378a1b50ebac46e5b7b8a8bef4b5365ac15bef9d..0000000000000000000000000000000000000000
--- a/raytime3d/getModelInfo.c
+++ /dev/null
@@ -1,109 +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 "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))
-
-/**
-* reads gridded model file to compute minimum and maximum values and sampling intervals
-*
-* AUTHOR:
-* Jan Thorbecke (janth@xs4all.nl)
-* The Netherlands
-**/
-
-int getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose)
-{
- FILE *fp;
- size_t nread, trace_sz;
- off_t bytes;
- int ret, i, one_shot, ntraces;
- 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 = hdr.d1;
- *d2 = hdr.d2;
- *f1 = hdr.f1;
- *f2 = hdr.f2;
-
- if ( NINT(100.0*((*d1)/(*d2)))!=100 ) {
- verr("dx and dz are different in the model !");
- }
- if ( NINT(1000.0*(*d1))==0 ) {
- if(!getparfloat("dx",d1)) {
- verr("dx is equal to zero use parameter dx= to set value");
- }
- *d2 = *d1;
- }
- trace_sz = sizeof(float)*(*n1)+TRCBYTES;
- ntraces = (int) (bytes/trace_sz);
- *n2 = ntraces;
-
- /* check to find out min and max values gather */
-
- one_shot = 1;
- trace = (float *)malloc(trace_sz);
- fseeko( fp, TRCBYTES, SEEK_SET );
- nread = fread( trace, sizeof(float), hdr.ns, fp );
- assert (nread == hdr.ns);
- fseeko( fp, TRCBYTES, SEEK_SET );
-
- if (hdr.trid == TRID_DEPTH) *axis = 1; /* samples are z-axis */
- else *axis = 0; /* sample direction respresents the x-axis */
-
- i=0; cmin=trace[0];
- while ( ( (cmin==0.0) && zeroch) && (i<hdr.ns) ) cmin=trace[i++];
-
- *max = cmin;
- *min = cmin;
- /* keep on reading traces until there are no more traces (nread==0) */
- while (one_shot) {
- nread = fread( trace, sizeof(float), hdr.ns, fp );
- assert (nread == hdr.ns);
- for (i=0;i<(*n1);i++) {
- *max = MAX(trace[i],*max);
- cmin = MIN(trace[i],*min);
- if (zeroch) {
- if (cmin!=0.0) *min = MIN(*min, cmin);
- }
- else {
- *min = cmin;
- }
- }
- nread = fread( &hdr, 1, TRCBYTES, fp );
- if (nread==0) break;
- }
- fclose(fp);
- free(trace);
-
- if (verbose>2) {
- vmess("For file %s", file_name);
- vmess("nz=%d nx=%d", *n1, *n2);
- vmess("dz=%f dx=%f", *d1, *d2);
- vmess("min=%f max=%f", *min, *max);
- vmess("zstart=%f xstart=%f", *f1, *f2);
- if (*axis) vmess("sample represent z-axis\n");
- else vmess("sample represent x-axis\n");
- }
- return 0;
-}
-
diff --git a/raytime3d/getParameters.c b/raytime3d/getParameters.c
deleted file mode 100644
index 89d1d34585a977f9ba6fd58eff936bf4396e4ea2..0000000000000000000000000000000000000000
--- a/raytime3d/getParameters.c
+++ /dev/null
@@ -1,305 +0,0 @@
-#include<stdlib.h>
-#include<stdio.h>
-#include<math.h>
-#include<assert.h>
-#include"par.h"
-#include"raytime.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))
-
-/**
-*
-* The routine getParameters reads in all parameters to set up a FD modeling.
-* Model and source parameters are used to calculate stability and dispersion relations
-* Source and receiver positions are calculated and checked if they fit into the model.
-*
-* AUTHOR:
-* Jan Thorbecke (janth@xs4all.nl)
-* The Netherlands
-**/
-
-int getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose);
-
-int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx, int nz);
-
-int getParameters(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose)
-{
- int nx, nz, nsrc, ix, axis, is0;
- int idzshot, idxshot;
- int src_ix0, src_iz0, src_ix1, src_iz1;
- float cp_min, cp_max;
- float sub_x0,sub_z0;
- float srcendx, srcendz, dx, dz;
- float xsrc, zsrc, dxshot, dzshot;
- float dxrcv,dzrcv,dxspread,dzspread;
- float xmax, zmax;
- float xsrc1, xsrc2, zsrc1, zsrc2;
- float *xsrca, *zsrca;
- float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
- size_t nsamp;
- int nxsrc, nzsrc;
- int is;
- char *src_positions;
-
- if (!getparint("verbose",&verbose)) verbose=0;
-
- if (!getparstring("file_cp",&mod->file_cp)) {
- verr("parameter file_cp required!");
- }
- if (!getparstring("file_rcv",&rec->file_rcv)) rec->file_rcv="recv.su";
- if (!getparint("src_at_rcv",&src->src_at_rcv)) src->src_at_rcv=1;
-
- /* read model parameters, which are used to set up source and receivers and check stability */
-
- getModelInfo(mod->file_cp, &nz, &nx, &dz, &dx, &sub_z0, &sub_x0, &cp_min, &cp_max, &axis, 1, verbose);
- mod->cp_max = cp_max;
- mod->cp_min = cp_min;
- mod->dz = dz;
- mod->dx = dx;
- mod->nz = nz;
- mod->nx = nx;
-
- /* origin of model in real (non-grid) coordinates */
- mod->x0 = sub_x0;
- mod->z0 = sub_z0;
- xmax = sub_x0+(nx-1)*dx;
- zmax = sub_z0+(nz-1)*dz;
-
- if (verbose) {
- vmess("*******************************************");
- vmess("*************** model info ****************");
- vmess("*******************************************");
- vmess("nz = %8d nx = %8d", nz, nx);
- vmess("dz = %8.4f dx = %8.4f", dz, dx);
- vmess("zmin = %8.4f zmax = %8.4f", sub_z0, zmax);
- vmess("xmin = %8.4f xmax = %8.4f", sub_x0, xmax);
- vmess("min(cp) = %9.3f max(cp) = %9.3f", cp_min, cp_max);
- }
-
- /* define the number and type of shots to model */
- /* each shot can have multiple sources arranged in different ways */
-
- if (!getparfloat("xsrc",&xsrc)) xsrc=sub_x0+((nx-1)*dx)/2.0;
- if (!getparfloat("zsrc",&zsrc)) zsrc=sub_z0;
- if (!getparint("nxshot",&shot->nx)) shot->nx=1;
- if (!getparint("nzshot",&shot->nz)) shot->nz=1;
- if (!getparfloat("dxshot",&dxshot)) dxshot=dx;
- if (!getparfloat("dzshot",&dzshot)) dzshot=dz;
-
- shot->n = (shot->nx)*(shot->nz);
-
- if (shot->nx>1) {
- idxshot=MAX(0,NINT(dxshot/dx));
- }
- else {
- idxshot=0.0;
- }
- if (shot->nz>1) {
- idzshot=MAX(0,NINT(dzshot/dz));
- }
- else {
- idzshot=0.0;
- }
-
- /* calculate the shot positions */
-
- src_ix0=MAX(0,NINT((xsrc-sub_x0)/dx));
- src_ix0=MIN(src_ix0,nx);
- src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
- src_iz0=MIN(src_iz0,nz);
- srcendx=(shot->nx-1)*dxshot+xsrc;
- srcendz=(shot->nz-1)*dzshot+zsrc;
- src_ix1=MAX(0,NINT((srcendx-sub_x0)/dx));
- src_ix1=MIN(src_ix1,nx);
- src_iz1=MAX(0,NINT((srcendz-sub_z0)/dz));
- src_iz1=MIN(src_iz1,nz);
-
- shot->x = (int *)calloc(shot->nx,sizeof(int));
- shot->z = (int *)calloc(shot->nz,sizeof(int));
- for (is=0; is<shot->nx; is++) {
- shot->x[is] = src_ix0+is*idxshot;
- if (shot->x[is] > nx-1) shot->nx = is-1;
- }
- for (is=0; is<shot->nz; is++) {
- shot->z[is] = src_iz0+is*idzshot;
- if (shot->z[is] > nz-1) shot->nz = is-1;
- }
-
- /* check if source array is defined */
-
- nxsrc = countparval("xsrca");
- nzsrc = countparval("zsrca");
- if (nxsrc != nzsrc) {
- verr("Number of sources in array xsrca (%d), zsrca(%d) are not equal",nxsrc, nzsrc);
- }
-
- /* check if sources on a circle are defined */
-
- if (getparfloat("rsrc", &rsrc)) {
- if (!getparfloat("dphisrc",&dphisrc)) dphisrc=2.0;
- if (!getparfloat("oxsrc",&oxsrc)) oxsrc=0.0;
- if (!getparfloat("ozsrc",&ozsrc)) ozsrc=0.0;
- ncsrc = NINT(360.0/dphisrc);
- src->n = nsrc;
-
- src->x = (int *)malloc(ncsrc*sizeof(int));
- src->z = (int *)malloc(ncsrc*sizeof(int));
-
- for (ix=0; ix<ncsrc; ix++) {
- src->x[ix] = NINT((oxsrc-sub_x0+rsrc*cos(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dx);
- src->z[ix] = NINT((ozsrc-sub_z0+rsrc*sin(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dz);
- if (verbose>4) fprintf(stderr,"Source on Circle: xsrc[%d]=%d zsrc=%d\n", ix, src->x[ix], src->z[ix]);
- }
-
- }
-
- /* TO DO propagate src_positions parameter and structure through code */
-
- if (!getparstring("src_positions",&src_positions)) src_positions="single";
- src->random=0;
- src->plane=0;
- src->array=0;
- src->single=0;
- if (strstr(src_positions, "single")) src->single=1;
- else if (strstr(src_positions, "array")) src->array=1;
- else if (strstr(src_positions, "random")) src->random=1;
- else if (strstr(src_positions, "plane")) src->plane=1;
- else src->single=1;
-
- /* to maintain functionality of older parameters usage */
- if (!getparint("src_random",&src->random)) src->random=0;
- if (!getparint("plane_wave",&src->plane)) src->plane=0;
-
- if (src->random) {
- src->plane=0;
- src->array=0;
- src->single=0;
- }
- if (src->plane) {
- src->random=0;
- src->array=0;
- src->single=0;
- }
-
-
- /* number of sources per shot modeling */
-
- if (!getparint("src_window",&src->window)) src->window=0;
- if (!getparint("distribution",&src->distribution)) src->distribution=0;
- if (!getparfloat("amplitude", &src->amplitude)) src->amplitude=0.0;
- if (src->random && nxsrc==0) {
- if (!getparint("nsrc",&nsrc)) nsrc=1;
- if (!getparfloat("xsrc1", &xsrc1)) xsrc1=sub_x0;
- if (!getparfloat("xsrc2", &xsrc2)) xsrc2=xmax;
- if (!getparfloat("zsrc1", &zsrc1)) zsrc1=sub_z0;
- if (!getparfloat("zsrc2", &zsrc2)) zsrc2=zmax;
- dxshot = xsrc2-xsrc1;
- dzshot = zsrc2-zsrc1;
- src->x = (int *)malloc(nsrc*sizeof(int));
- src->z = (int *)malloc(nsrc*sizeof(int));
- nsamp = 0;
-
- }
- else if (nxsrc != 0) {
- /* source array is defined */
- nsrc=nxsrc;
- src->x = (int *)malloc(nsrc*sizeof(int));
- src->z = (int *)malloc(nsrc*sizeof(int));
- xsrca = (float *)malloc(nsrc*sizeof(float));
- zsrca = (float *)malloc(nsrc*sizeof(float));
- getparfloat("xsrca", xsrca);
- getparfloat("zsrca", zsrca);
- for (is=0; is<nsrc; is++) {
- src->x[is] = NINT((xsrca[is]-sub_x0)/dx);
- src->z[is] = NINT((zsrca[is]-sub_z0)/dz);
- if (verbose>3) fprintf(stderr,"Source Array: xsrc[%d]=%f zsrc=%f\n", is, xsrca[is], zsrca[is]);
- }
- src->random = 1;
- free(xsrca);
- free(zsrca);
- }
- else {
- if (src->plane) { if (!getparint("nsrc",&nsrc)) nsrc=1;}
- else nsrc=1;
-
- if (nsrc > nx) {
- vwarn("Number of sources used in plane wave is larger than ");
- vwarn("number of gridpoints in X. Plane wave will be clipped to the edges of the model");
- nsrc = mod->nx;
- }
-
- /* for a source defined on mutliple gridpoint calculate p delay factor */
-
- src->x = (int *)malloc(nsrc*sizeof(int));
- src->z = (int *)malloc(nsrc*sizeof(int));
- is0 = -1*floor((nsrc-1)/2);
- for (is=0; is<nsrc; is++) {
- src->x[is] = is0 + is;
- src->z[is] = 0;
- }
-
- }
-
- src->n=nsrc;
-
- if (verbose) {
- if (src->n>1) {
- vmess("*******************************************");
- vmess("*********** source array info *************");
- vmess("*******************************************");
- vmess("Areal source array is defined with %d sources.",nsrc);
- vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
- }
- if (src->random) {
- vmess("Sources are placed at random locations in domain: ");
- vmess(" x[%.2f : %.2f] z[%.2f : %.2f] ", xsrc1, xsrc2, zsrc1, zsrc2);
- }
- }
-
- /* define receivers */
-
- if (!getparint("sinkdepth",&rec->sinkdepth)) rec->sinkdepth=0;
- if (!getparint("sinkdepth_src",&src->sinkdepth)) src->sinkdepth=0;
- if (!getparint("sinkvel",&rec->sinkvel)) rec->sinkvel=0;
- if (!getparint("max_nrec",&rec->max_nrec)) rec->max_nrec=15000;
- if (!getparfloat("dxspread",&dxspread)) dxspread=0;
- if (!getparfloat("dzspread",&dzspread)) dzspread=0;
- if (!getparint("nt",&rec->nt)) rec->nt=1024;
-
- /* calculates the receiver coordinates */
-
- recvPar(rec, sub_x0, sub_z0, dx, dz, nx, nz);
-
- if (verbose) {
- if (rec->n) {
- dxrcv = rec->xr[MIN(1,rec->n-1)]-rec->xr[0];
- dzrcv = rec->zr[MIN(1,rec->n-1)]-rec->zr[0];
- vmess("*******************************************");
- vmess("************* receiver info ***************");
- vmess("*******************************************");
- vmess("ntrcv = %d nrcv = %d ", rec->nt, rec->n);
- vmess("dzrcv = %f dxrcv = %f ", dzrcv, dxrcv);
- vmess("Receiver array at coordinates: ");
- vmess("zmin = %f zmax = %f ", rec->zr[0]+sub_z0, rec->zr[rec->n-1]+sub_z0);
- vmess("xmin = %f xmax = %f ", rec->xr[0]+sub_x0, rec->xr[rec->n-1]+sub_x0);
- vmess("which are gridpoints: ");
- vmess("izmin = %d izmax = %d ", rec->z[0], rec->z[rec->n-1]);
- vmess("ixmin = %d ixmax = %d ", rec->x[0], rec->x[rec->n-1]);
- fprintf(stderr,"\n");
- }
- else {
- vmess("*************** no receivers **************");
- }
- }
-
- /* Ray tracing parameters */
- if (!getparint("smoothwindow",&ray->smoothwindow)) ray->smoothwindow=0;
- if (!getparint("useT2",&ray->useT2)) ray->useT2=0;
- if (!getparint("geomspread",&ray->geomspread)) ray->geomspread=1;
- if (!getparint("nraystep",&ray->nray)) ray->nray=5;
-
- return 0;
-}
-
diff --git a/raytime3d/raytime.h b/raytime3d/raytime.h
deleted file mode 100644
index 882c78b9a2eb0b945737dac10ae5c6f07aa55123..0000000000000000000000000000000000000000
--- a/raytime3d/raytime.h
+++ /dev/null
@@ -1,135 +0,0 @@
-#include<stdlib.h>
-#include<stdio.h>
-#include<limits.h>
-#include<float.h>
-#include<math.h>
-
-#ifndef COMPLEX
-typedef struct _complexStruct { /* complex number */
- float r,i;
-} complex;
-typedef struct _dcomplexStruct { /* complex number */
- double r,i;
-} dcomplex;
-#endif/* complex */
-
-
-typedef struct _icoord { /* 3D coordinate integer */
- int z;
- int x;
- int y;
-} icoord;
-
-typedef struct _fcoord { /* 3D coordinate float */
- float z;
- float x;
- float y;
-} fcoord;
-
-struct s_ecount {
- int corner,corner_min,side;
-};
-
-typedef struct _receiverPar { /* Receiver Parameters */
- char *file_rcv;
- int n;
- int nt;
- int max_nrec;
- int *z;
- int *x;
- float *zr;
- float *xr;
- int scale;
- int sinkdepth;
- int sinkvel;
- float cp;
- float rho;
-} recPar;
-
-typedef struct _modelPar { /* Model Parameters */
- int sh;
- char *file_cp;
- float dz;
- float dx;
- float dt;
- float z0;
- float x0;
- /* medium max/min values */
- float cp_min;
- float cp_max;
- int nz;
- int nx;
-} modPar;
-
-typedef struct _waveletPar { /* Wavelet Parameters */
- char *file_src; /* general source */
- int nsrcf;
- int nt;
- int ns;
- int nx;
- float dt;
- float ds;
- float fmax;
- int random;
- int seed;
- int nst;
- size_t *nsamp;
-} wavPar;
-
-typedef struct _sourcePar { /* Source Array Parameters */
- int n;
- int type;
- int orient;
- int *z;
- int *x;
- int single;
- int plane;
- int circle;
- int array;
- int random;
- int multiwav;
- float angle;
- float velo;
- float amplitude;
- int distribution;
- int window;
- int injectionrate;
- int sinkdepth;
- int src_at_rcv; /* Indicates that wavefield should be injected at receivers */
-} srcPar;
-
-typedef struct _shotPar { /* Shot Parameters */
- int n;
- int nx;
- int nz;
- int *z;
- int *x;
-} shotPar;
-
-typedef struct _raypar { /* ray-tracing parameters */
- int smoothwindow;
- int useT2;
- int geomspread;
- int nray;
-} rayPar;
-
-#ifndef TRUE
-# define TRUE 1
-#endif
-
-#ifndef FALSE
-# define FALSE 0
-#endif
-
-#define equal(x,y) !strcmp(x,y)
-#define min2(a,b) (((a) < (b)) ? (a) : (b))
-#define max2(a,b) (((a) > (b)) ? (a) : (b))
-
-#define Infinity FLT_MAX
-
-#if __STDC_VERSION__ >= 199901L
- /* "restrict" is a keyword */
-#else
-#define restrict
-#endif
-
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
index 9b0c88d6e59fe63c13275cc139f7352bff266501..73e57f41822da6bec2322f85b0b943e25afad73f 100644
--- a/raytime3d/raytime3d.c
+++ b/raytime3d/raytime3d.c
@@ -1,15 +1,56 @@
-#include <DELPHI_IOc.h>
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
+#include <genfft.h>
+#include"par.h"
+#include"raytime3d.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))
+
+double wallclock_time(void);
+
+void name_ext(char *filename, char *extension);
+
+void threadAffinity(void);
+
+
+int getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose);
+
+int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord r, rayPar ray, fcoord *T, float *Jr);
+
+void applyMovingAverageFilter(float *slowness, icoord size, int window, int dim, float *averageModel);
+
+int readModel3d(char *file_name, float *slowness, int n1, int n2, int n3, int nz, int nx, int ny, float h, int verbose);
+
+int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int reverse, int verbose);
+
+int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
+
+void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int xs, int ys, int zs, int NCUBE);
+
+void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube);
+
/*********************** self documentation **********************/
char *sdoc[] = {
" ",
" RAYTIME3D - modeling of one-way traveltime for operators in 3D media",
" ",
-" raytime3d file_vel= xsrc1= zsrc1= ysrc1= [optional parameters]",
+" raytime3d file_cp= xsrc1= zsrc1= ysrc1= [optional parameters]",
" ",
" Required parameters:",
" ",
-" file_vel= ................ gridded velocity file ",
+" file_cp= ................ gridded velocity file ",
+" file_src= ......... file with source signature",
+" file_rcv=recv.su .. base name for receiver files",
+" file_rcvtime= ..... receiver file in x-t",
+" h= ................ read from model file: if d1==0 then h= can be used to set it",
+" nt=1024 ........... number of time-samples in file_rcvtime",
" xsrc1= ................... x-position of the source (m)",
" ysrc1= ................... y-position of the source (m)",
" zsrc1= ................... z-position of the source (m)",
@@ -18,16 +59,17 @@ char *sdoc[] = {
" ",
" INPUT AND OUTPUT ",
" key=gy ................... input data sorting key",
-" ny=1 ..................... if 2D file number of y traces (2D model)",
-" nxmax=512 ................ maximum number of traces in input files",
-" ntmax=1024 ............... maximum number of samples/trace in input files",
+" nx=1 ..................... if 1D file number of points in x ",
+" ny=1 ..................... if 2D file number of points in y ",
" file_out= ................ output file with traveltime cube",
" file_amp= ................ output file with approximate amplitudes",
-" RAY TRACING ",
-" dT=0 ..................... put traces on one-way time grid with step dT",
-" Tmin=first shot........... minimum time of one-way time grid",
-" Tmax=last shot ........... maximum time of one-way time grid",
-" hom=1 .................... 1: draw straight rays in homogeneous layers",
+" ",
+//" RAY TRACING PARAMETERS:",
+//" dT=0 ..................... put traces on one-way time grid with step dT",
+//" Tmin=first shot........... minimum time of one-way time grid",
+//" Tmax=last shot ........... maximum time of one-way time grid",
+//" hom=1 .................... 1: draw straight rays in homogeneous layers",
+//" ",
" SOURCE POSITIONS ",
" xsrc2=xsrc1 .............. x-position of last source",
" dxsrc=0 .................. step in source x-direction",
@@ -36,11 +78,11 @@ char *sdoc[] = {
" zsrc2=zsrc1 .............. z-position of last source",
" dzsrc=0 .................. step in source z-direction",
" RECEIVER POSITIONS ",
-" xrcv=0,(nx-1)*dx ......... x-position's of receivers (array)",
-" yrcv=0,(ny-1)*dy ......... y-position's of receivers (array)",
+" xrcv=-(nx-1/2)*h,(nx-1/2)*h .. x-position's of receivers (array)",
+" yrcv=-(ny-1)/2*h,(ny-1/2)*h .. y-position's of receivers (array)",
" zrcv=0,0 ................. z-position's of receivers (array)",
-" dxrcv=dx ................. step in receiver x-direction",
-" dyrcv=dy ................. step in receiver y-direction",
+" dxrcv=h ................. step in receiver x-direction",
+" dyrcv=h ................. step in receiver y-direction",
" dzrcv=0 .................. step in receiver z-direction",
" dxspr=0 .................. step of receiver spread in x-direction",
" dyspr=0 .................. step of receiver spread in y-direction",
@@ -65,11 +107,13 @@ NULL};
#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
#define s0(x,y,z) slow0[nxy*(z) + nx*(y) + (x)]
-void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int xs, int ys, int zs, int NCUBE);
-void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube);
-
void main(int argc, char *argv[])
{
+ modPar mod;
+ recPar rec;
+ srcPar src;
+ shotPar shot;
+ rayPar ray;
int
nx, /* x-dimension of mesh (LEFT-TO-RIGHT) */
ny, /* y-dimension of mesh (FRONT-TO-BACK) */
@@ -80,13 +124,11 @@ void main(int argc, char *argv[])
h, /* spatial mesh interval (units consistant with vel) */
*slow0, *time0;
-/* to read the delphi velocity file */
- int32 type, dom1, dom2;
- int error, n1, n2, ret, size, nkeys, verbose;
- int ntmax, nxmax;
- float d1, d2, f1, f2, *tmpdata, c, scl, ox, oz, oy;
- char *file_vel, *file_out, *key;
- segyhdr *hdrs;
+/* to read the velocity file */
+ int error, n1, n2, n3, ret, size, nkeys, verbose;
+ float d1, d2, d3, f1, f2, f3, *tmpdata, c, scl, ox, oz, oy;
+ char *file_cp, *file_out;
+ segy *hdrs;
/*---------------------------------------------------------------------------*
* Read input parameters and query for any that are needed but missing.
@@ -97,122 +139,57 @@ void main(int argc, char *argv[])
if (!getparint("verbose",&verbose)) verbose=0;
if (verbose) {
- samess("Hole, J.A., and B.C. Zelt, 1995. \"3-D finite-difference");
- samess("reflection traveltimes\". Geophys. J. Int., 121, 427-434");
- }
- if(!getparstring("file_out",&file_out)) saerr("file_out not given");
- if(!getparstring("file_vel", &file_vel)) {
- sawarn("file_vel not defined, assuming homogeneous model");
- if(!getparfloat("c",&c)) saerr("c not defined");
- if(!getparint("nx",&nx)) saerr("nx not defined");
- if(!getparint("ny",&ny)) saerr("ny not defined");
- if(!getparint("nz",&nz)) saerr("nz not defined");
- if(!getparfloat("h",&h)) saerr("h not defined");
+ vmess("Hole, J.A., and B.C. Zelt, 1995. \"3-D finite-difference");
+ vmess("reflection traveltimes\". Geophys. J. Int., 121, 427-434");
}
- if(!getparint("ny",&ny)) saerr("ny not defined");
+ if(!getparstring("file_out",&file_out)) verr("file_out not given");
+
+ getParameters3d(&mod, &rec, &src, &shot, &ray, verbose);
+
/*---------------------------------------------------------------------------*
* Open velocity file
*---------------------------------------------------------------------------*/
- if (file_vel != NULL) {
- error = open_file(file_vel, GUESS_TYPE, DELPHI_READ);
- if (error < 0 ) saerr("error in opening file %s", file_vel);
- error = get_dims(file_vel, &n1, &n2, &type);
- if (ret >= 0) {
- if (!getparint("ntmax", &ntmax)) ntmax = n1;
- if (!getparint("nxmax", &nxmax)) nxmax = n2;
- if (verbose>=2 && (ntmax!=n1 || nxmax!=n2))
- samess("dimensions overruled: %d x %d",ntmax,nxmax);
+ if (file_cp != NULL) {
+
+ if (n2==1) { /* 1D model */
+ if(!getparint("nx",&nx)) verr("for 1D medium nx not defined");
+ if(!getparint("ny",&nx)) verr("for 1D medium ny not defined");
+ nz = n1;
+ oz = f1; ox = ((nx-1)/2)*d1; oy = ((ny-1)/2)*d1;
+ dz = d1; dx = d1; dy = d1;
}
- else {
- if (!getparint("ntmax", &ntmax)) ntmax = 1024;
- if (!getparint("nxmax", &nxmax)) nxmax = 512;
- if (verbose>=2) samess("dimensions used: %d x %d",ntmax,nxmax);
+ else if (n3==1) { /* 2D model */
+ if(!getparint("ny",&nx)) verr("for 2D medium ny not defined");
+ nz = n1; nx = n2;
+ oz = f1; ox = f2; oy = ((ny-1)/2)*d1;
+ dz = d1; dx = d1; dy = d1;
}
- size = ntmax * nxmax;
- tmpdata = alloc1float(size);
- hdrs = (segyhdr *) malloc(nxmax*sizeof(segyhdr));
-
- if (!getparstring("key", &key)) {
- ret = get_sort(file_vel);
- if (ret < 0) key = "gy";
- else key = getkey(ret);
+ else { /* Full 3D model */
+ nz = n1; nx = n2; nz = n3;
+ oz = f1; ox = f2; oy = f3;
+ dz = d1; dx = d1; dy = d1;
}
- if (verbose) samess("input sorting key is %s",key);
- set_sukey(key);
-
- ret = read_data(file_vel,tmpdata,size,&n1,&n2,&f1,&f2,&d1,&d2,
- &type,hdrs);
- if (ret < 0) saerr("error in reading data from file %s", file_vel);
- if (hdrs[0].scalco < 0) scl = 1.0/fabs(hdrs[0].scalco);
- else if (hdrs[0].scalco == 0) scl = 1.0;
- else scl = hdrs[0].scalco;
- get_axis(&dom1, &dom2);
- if (d1 != d2)
- saerr("d1 != d2; this is not allowed in the calculation");
+
h = d1;
- if (dom1 == SA_AXIS_Z) {
- nx = n2; nz = n1;
- ox = hdrs[0].gx*scl; oy = hdrs[0].gy*scl; oz = f1;
- }
- else {
- nx = n1; nz = n2;
- ox = f1; oy = hdrs[0].gy*scl; oz = f1;
- }
+ slow0 = (float *)malloc(nz*nx*ny*sizeof(float));
+ if (slow0 == NULL) verr("Out of memory for slow0 array!");
- slow0 = (float *)malloc(ny*n1*n2*sizeof(float));
- if (slow0 == NULL) saerr("Out of memory for slow0 array!");
- nxy = nx * ny;
- if (verbose) samess("h = %.2f nx = %d nz = %d ny = %d", h, nx, nz, ny);
-
- yy = 0;
- while (ret >= 0) {
- if (verbose==2) disp_info(file_vel,n1,n2,f1,f2,d1,d2,type);
-
- if (dom1 == SA_AXIS_Z) {
- if (n2 != nx || n1 != nz) saerr("dimensions changed");
- for (i = 0; i < n2; i++) {
- for (j = 0; j < n1; j++)
- slow0[j*nxy+yy*nx+i] = h/tmpdata[i*n1+j];
- }
- }
- else {
- if (n1 != nx || n2 != nz) saerr("dimensions changed");
- for (i = 0; i < n2; i++) {
- for (j = 0; j < n1; j++)
- slow0[i*nxy+yy*nx+j] = h/tmpdata[i*n1+j];
- }
- }
+ readModel3d(file_cp, slow0, n1, n2, n3, nz, nx, ny, h, verbose);
- yy += 1;
- ret = read_data(file_vel, tmpdata, size, &n1, &n2, &f1, &f2,
- &d1, &d2, &type, hdrs);
- }
- ret = close_file(file_vel);
- if (ret < 0) sawarn("err %d on closing input file",ret);
- free1float(tmpdata);
- if (yy == 1) {
- if(!getparint("ny",&ny)) samess("2D model defined");
- else {
- slow0 = (float *)realloc(slow0, ny*nx*nz*sizeof(float));
- if (slow0 == NULL) saerr("Out of memory for slow0 array!");
-
- samess("3D model defined from 2D model");
- for (zz = 0; zz < nz; zz++) {
- for (yy = 1; yy < ny; yy++) {
- for (xx = 0; xx < nx; xx++)
- slow0[zz*nxy+yy*nx+xx] = slow0[zz*nxy+xx];
- }
- }
- }
- }
+ if (verbose) vmess("h = %.2f nx = %d nz = %d ny = %d", h, nx, nz, ny);
}
else {
nxy = nx * ny;
+ if(!getparint("nx",&nx)) verr("for homogenoeus medium nx not defined");
+ if(!getparint("ny",&nx)) verr("for homogenoeus medium ny not defined");
+ if(!getparint("nz",&nx)) verr("for homogenoeus medium nz not defined");
+ oz = 0; ox = ((nx-1)/2)*d1; oy = ((ny-1)/2)*d1;
+
slow0 = (float *)malloc(nx*nz*ny*sizeof(float));
- if (slow0 == NULL) saerr("Out of memory for slow0 array!");
+ if (slow0 == NULL) verr("Out of memory for slow0 array!");
scl = h/c;
ox = 0; oy = 0; oz = 0;
for (zz = 0; zz < nz; zz++) {
@@ -228,7 +205,7 @@ void main(int argc, char *argv[])
/* ALLOCATE MAIN GRID FOR TIMES */
time0 = (float *) malloc(sizeof(float)*nxyz);
- if(time0 == NULL) saerr("error in allocation of array time0");
+ if(time0 == NULL) verr("error in allocation of array time0");
/*---------------------------------------------------------------------------*
* Input the source locations.
@@ -237,7 +214,7 @@ void main(int argc, char *argv[])
*---------------------------------------------------------------------------*/
src3d(time0, slow0, nz, nx, ny, h, ox, oy, oz, &xs, &ys, &zs, &cube);
- if (verbose) samess("source positions xs = %d ys = %d zs = %d", xs,ys,zs);
+ if (verbose) vmess("source positions xs = %d ys = %d zs = %d", xs,ys,zs);
/* for (zz = 0; zz < nz; zz++) {
for (yy = 0; yy < ny; yy++) {
@@ -255,7 +232,6 @@ void main(int argc, char *argv[])
* Check and set parameters
*---------------------------------------------------------------------------*/
-
/*---------------------------------------------------------------------------*
* Compute traveltimes.
*---------------------------------------------------------------------------*/
@@ -274,17 +250,17 @@ void main(int argc, char *argv[])
}
}
*/
- ret = open_file(file_out, GUESS_TYPE, DELPHI_CREATE);
- if (ret < 0 ) saerr("error in creating output file %s", file_out);
+// ret = open_file(file_out, GUESS_TYPE, DELPHI_CREATE);
+// if (ret < 0 ) verr("error in creating output file %s", file_out);
- hdrs = (segyhdr *) malloc(ny*sizeof(segyhdr));
- tmpdata = alloc1float(nxy);
+ hdrs = (segy *) malloc(ny*sizeof(segy));
+ tmpdata = (float *)malloc(nxy*sizeof(float));
f1 = ox;
f2 = oy;
d1 = h;
d2 = h;
- gen_hdrs(hdrs,nx,ny,f1,f2,d1,d2,TRID_ZX);
+// gen_hdrs(hdrs,nx,ny,f1,f2,d1,d2,TRID_ZX);
for (i = 0; i < ny; i++) {
hdrs[i].scalco = -1000;
hdrs[i].scalel = -1000;
@@ -300,18 +276,17 @@ void main(int argc, char *argv[])
}
}
- if (ret < 0 ) sawarn("error on writing keys.");
- ret = set_axis(dom1, dom2);
- if (ret < 0 ) saerr("error on writing axis.");
+/*
ret = write_data(file_out,tmpdata,nx,ny,f1,f2,d1,d2,type,hdrs);
- if (ret < 0 ) saerr("error on writing output file.");
+ if (ret < 0 ) verr("error on writing output file.");
ret = close_file(file_out);
- if (ret < 0) saerr("err %d on closing output file",ret);
+ if (ret < 0) verr("err %d on closing output file",ret);
+*/
free(time0);
free(slow0);
free(hdrs);
- free1float(tmpdata);
+ free(tmpdata);
exit(0);
diff --git a/raytime3d/readModel.c b/raytime3d/readModel.c
deleted file mode 100644
index 27f1da715a50dbc36eac4827efe68434288ee9e7..0000000000000000000000000000000000000000
--- a/raytime3d/readModel.c
+++ /dev/null
@@ -1,80 +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 "segy.h"
-#include "par.h"
-#include "raytime.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))
-
-/**
-* Reads gridded model files and compute from them medium parameters used in the FD kernels.
-* The files read in contain the P (and S) wave velocity and density.
-* The medium parameters calculated are lambda, mu, lambda+2mu, and 1/ro.
-*
-* AUTHOR:
-* Jan Thorbecke (janth@xs4all.nl)
-* The Netherlands
-**/
-
-
-int readModel(modPar mod, float *velocity, float *slowness, int nw)
-{
- FILE *fpcp;
- size_t nread;
- int i, tracesToDo, j;
- int nz, nx;
- segy hdr;
-
-
- /* grid size and start positions for the components */
- nz = mod.nz;
- nx = mod.nx;
-
-/* open files and read first header */
-
- fpcp = fopen( mod.file_cp, "r" );
- assert( fpcp != NULL);
- nread = fread(&hdr, 1, TRCBYTES, fpcp);
- assert(nread == TRCBYTES);
-
-/* read all traces */
-
- tracesToDo = mod.nx;
- i = 0;
- while (tracesToDo) {
- nread = fread(&velocity[i*nz], sizeof(float), hdr.ns, fpcp);
- assert (nread == hdr.ns);
- for (j=0;j<nz;j++) {
- if (velocity[i*nz+j]!=0.0) {
- slowness[(i+nw)*nz+j+nw] = 1.0/velocity[i*nz+j];
- }
- }
- for (j=0;j<nw;j++) slowness[(i+nw)*nz+j] = slowness[(i+nw)*nz+nw];
- for (j=nz+nw;j<nz+2*nw;j++) slowness[(i+nw)*nz+j] = slowness[(i+nw)*nz+nz+nw-1];
-
- nread = fread(&hdr, 1, TRCBYTES, fpcp);
- if (nread==0) break;
- i++;
- }
- fclose(fpcp);
-
- for (i=0;i<nw;i++) {
- for (j=0;j<nz+2*nw;j++) {
- slowness[(i)*nz+j] = slowness[(nw)*nz+j];
- slowness[(nx+nw+i)*nz+j] = slowness[(nx+nw-1)*nz+j];
- }
- }
-
- return 0;
-}
-
-
diff --git a/raytime3d/recvPar.c b/raytime3d/recvPar.c
deleted file mode 100644
index 4ad9eae9be2f3c7bcbeb1ee061dc97d5a4ff9d77..0000000000000000000000000000000000000000
--- a/raytime3d/recvPar.c
+++ /dev/null
@@ -1,519 +0,0 @@
-#include <stdio.h>
-#include <assert.h>
-#include <math.h>
-
-#include "raytime.h"
-#include "par.h"
-
-#define MIN(x,y) ((x) < (y) ? (x) : (y))
-#define MAX(x,y) ((x) > (y) ? (x) : (y))
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-
-/**
-* Calculates the receiver positions based on the input parameters
-*
-* AUTHOR:
-* Jan Thorbecke (janth@xs4all.nl)
-*
-* Ammendments:
-* Max Holicki changing the allocation receiver array (2-2016)
-* The Netherlands
-**/
-
-
-void name_ext(char *filename, char *extension);
-
-int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx, int nz)
-{
- float *xrcv1, *xrcv2, *zrcv1, *zrcv2;
- int i, ix, ir, verbose;
- float dxrcv, dzrcv, *dxr, *dzr;
- float rrcv, dphi, oxrcv, ozrcv, arcv;
- double circ, h, a, b, e, s, xr, zr, dr, srun, phase;
- float xrange, zrange, sub_x1, sub_z1;
- int Nx1, Nx2, Nz1, Nz2, Ndx, Ndz, iarray, nrec, nh;
- int nxrcv, nzrcv, ncrcv, nrcv, ntrcv, *nlrcv;
- float *xrcva, *zrcva;
- char* rcv_txt;
- FILE *fp;
-
- if (!getparint("verbose", &verbose)) verbose = 0;
-
- /* Calculate Model Dimensions */
- sub_x1=sub_x0+(nx-1)*dx;
- sub_z1=sub_z0+(nz-1)*dz;
-
-/* Compute how many receivers are defined and then allocate the receiver arrays */
-
- /* Receivers on a Circle */
- if (getparfloat("rrcv",&rrcv)) {
- if (!getparfloat("dphi",&dphi)) dphi=2.0;
- ncrcv=NINT(360.0/dphi);
- if (verbose) vmess("Total number of receivers on a circle: %d",ncrcv);
- }
- else {
- ncrcv=0;
- }
-
- /* Receivers from a File */
- ntrcv=0;
- if (!getparstring("rcv_txt",&rcv_txt)) rcv_txt=NULL;
- if (rcv_txt!=NULL) {
- /* Open text file */
- fp=fopen(rcv_txt,"r");
- assert(fp!=NULL);
- /* Get number of lines */
- while (!feof(fp)) if (fgetc(fp)=='\n') ntrcv++;
- fseek(fp,-1,SEEK_CUR);
- if (fgetc(fp)!='\n') ntrcv++; /* Checks if last line terminated by /n */
- if (verbose) vmess("Number of receivers in rcv_txt file: %d",ntrcv);
- rewind(fp);
- }
-
- /* Receiver Array */
- nxrcv=countparval("xrcva");
- nzrcv=countparval("zrcva");
- if (nxrcv!=nzrcv) verr("Number of receivers in array xrcva (%d), zrcva(%d) are not equal",nxrcv,nzrcv);
- if (verbose&&nxrcv) vmess("Total number of array receivers: %d",nxrcv);
-
- /* Linear Receiver Arrays */
- Nx1 = countparval("xrcv1");
- Nx2 = countparval("xrcv2");
- Nz1 = countparval("zrcv1");
- Nz2 = countparval("zrcv2");
- if (Nx1!=Nx2) verr("Number of receivers starting points in 'xrcv1' (%d) and number of endpoint in 'xrcv2' (%d) are not equal",Nx1,Nx2);
- if (Nz1!=Nz2) verr("Number of receivers starting points in 'zrcv1' (%d) and number of endpoint in 'zrcv2' (%d) are not equal",Nz1,Nz2);
- if (Nx1!=Nz2) verr("Number of receivers starting points in 'xrcv1' (%d) and number of endpoint in 'zrcv2' (%d) are not equal",Nx1,Nz2);
-
- rec->max_nrec=ncrcv+ntrcv+nxrcv;
-
- /* no receivers are defined use default linear array of receivers on top of model */
- if (!rec->max_nrec && Nx1==0) Nx1=1; // Default is to use top of model to record data
-
- if (Nx1) {
- /* Allocate Start & End Points of Linear Arrays */
- xrcv1=(float *)malloc(Nx1*sizeof(float));
- xrcv2=(float *)malloc(Nx1*sizeof(float));
- zrcv1=(float *)malloc(Nx1*sizeof(float));
- zrcv2=(float *)malloc(Nx1*sizeof(float));
- if (!getparfloat("xrcv1",xrcv1)) xrcv1[0]=sub_x0;
- if (!getparfloat("xrcv2",xrcv2)) xrcv2[0]=sub_x1;
- if (!getparfloat("zrcv1",zrcv1)) zrcv1[0]=sub_z0;
- if (!getparfloat("zrcv2",zrcv2)) zrcv2[0]=zrcv1[0];
-
- /* check if receiver arrays fit into model */
- for (iarray=0; iarray<Nx1; iarray++) {
- xrcv1[iarray] = MAX(sub_x0, xrcv1[iarray]);
- xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
- xrcv2[iarray] = MAX(sub_x0, xrcv2[iarray]);
- xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
-
- zrcv1[iarray] = MAX(sub_z0, zrcv1[iarray]);
- zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
- zrcv2[iarray] = MAX(sub_z0, zrcv2[iarray]);
- zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
- }
-
- /* Crop to Fit Model */
-/* Max's addtion still have to check if it has the same fucntionality */
- for (iarray=0;iarray<Nx1;iarray++) {
- if (xrcv1[iarray]<sub_x0) {
- if (xrcv2[iarray]<sub_x0) {
- verr("Linear array %d outside model bounds",iarray);
- }
- else {
- vwarn("Cropping element %d of 'xrcv1' (%f) to model bounds (%f)",iarray,xrcv1[iarray],sub_x0);
- xrcv1[iarray]=sub_x0;
- }
- }
- else if (xrcv1[iarray] > sub_x1) {
- verr("Linear array %d outside model bounds",iarray);
- }
- if ( (xrcv2[iarray] < xrcv1[iarray]) ) {
- verr("Ill defined linear array %d, 'xrcv1' (%f) greater than 'xrcv2' (%f)",iarray,xrcv1[iarray],xrcv2[iarray]);
- }
- else if (xrcv2[iarray]>sub_x1) {
- vwarn("Cropping element %d of 'xrcv2' (%f) to model bounds (%f)",iarray,xrcv2[iarray],sub_x1);
- xrcv2[iarray]=sub_x1;
- }
-
- if (zrcv1[iarray] < sub_z0) {
- if (zrcv2[iarray] < sub_z0) {
- verr("Linear array %d outside model bounds",iarray);
- }
- else {
- vwarn("Cropping element %d of 'zrcv1' (%f) to model bounds (%f)",iarray,zrcv1[iarray],sub_z0);
- zrcv1[iarray]=sub_z0;
- }
- }
- else if (zrcv1[iarray] > sub_z1) {
- verr("Linear array %d outside model bounds",iarray);
- }
- if ( (zrcv2[iarray] < zrcv1[iarray]) ) {
- verr("Ill defined linear array %d, 'zrcv1' (%f) greater than 'zrcv2' (%f)",iarray,zrcv1[iarray],zrcv2[iarray]);
- }
- else if (zrcv2[iarray]>sub_z1) {
- vwarn("Cropping element %d of 'xrcv2' (%f) to model bounds (%f)",iarray,zrcv2[iarray],sub_z1);
- zrcv2[iarray]=sub_z1;
- }
- }
-
- /* Get Sampling Rates */
- Ndx = countparval("dxrcv");
- Ndz = countparval("dzrcv");
-
- dxr = (float *)malloc(Nx1*sizeof(float));
- dzr = (float *)malloc(Nx1*sizeof(float));
- if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
- if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
- if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
- for (i=1; i<Nx1; i++) {
- dxr[i] = dxr[0];
- dzr[i] = dzr[0];
- }
- Ndx=1;
- Ndz=1;
- }
- else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
- for (i=1; i<Nx1; i++) {
- dxr[i] = dxr[0];
- dzr[i] = dzr[0];
- }
- Ndz=1;
- Ndx=1;
- }
- else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
- if (Ndx!=Ndz) {
- verr("Number of 'dxrcv' (%d) is not equal to number of 'dzrcv' (%d) or 1",Ndx,Ndz);
- }
- if (Ndx!=Nx1 && Ndx!=1) {
- verr("Number of 'dxrcv' (%d) is not equal to number of starting points in 'xrcv1' (%d) or 1",Ndx,Nx1);
- }
- }
-
- /* check consistency of receiver steps */
- for (iarray=0; iarray<Ndx; iarray++) {
- if (dxr[iarray]<0) {
- dxr[i]=dx;
- vwarn("'dxrcv' element %d (%f) is less than zero, changing it to %f'",iarray,dxr[iarray],dx);
- }
- }
- for (iarray=0;iarray<Ndz;iarray++) {
- if (dzr[iarray]<0) {
- dzr[iarray]=dz;
- vwarn("'dzrcv' element %d (%f) is less than zero, changing it to %f'",iarray,dzr[iarray],dz);
- }
- }
- for (iarray=0;iarray<Ndx;iarray++){
- if (dxr[iarray]==0 && dzr[iarray]==0) {
- xrcv2[iarray]=xrcv1[iarray];
- dxr[iarray]=1.;
- vwarn("'dxrcv' element %d & 'dzrcv' element 1 are both 0.",iarray+1);
- vmess("Placing 1 receiver at (%d,%d)",xrcv1[iarray],zrcv1[iarray]);
- }
- }
- for (iarray=0;iarray<Ndx;iarray++){
- if (xrcv1[iarray]==xrcv2[iarray] && dxr[iarray]!=0) {
- dxr[iarray]=0.;
- vwarn("Linear array %d: 'xrcv1'='xrcv2' and 'dxrcv' is not 0. Setting 'dxrcv'=0",iarray+1);
- }
- }
- for (iarray=0;iarray<Ndx;iarray++){
- if (zrcv1[iarray]==zrcv2[iarray] && dzr[iarray]!=0.){
- dzr[iarray]=0.;
- vwarn("Linear array %d: 'zrcv1'='zrcv2' and 'dzrcv' is not 0. Setting 'dzrcv'=0",iarray+1);
- }
- }
-
- /* Calculate Number of Receivers */
- nrcv = 0;
- nlrcv=(int *)malloc(Nx1*sizeof(int));
- for (iarray=0; iarray<Nx1; iarray++) {
- xrange = (xrcv2[iarray]-xrcv1[iarray]);
- zrange = (zrcv2[iarray]-zrcv1[iarray]);
- if (dxr[iarray] != 0.0) {
- nlrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
- }
- else {
- if (dzr[iarray] == 0) {
- verr("For receiver array %d: receiver distance dzrcv is not given", iarray);
- }
- nlrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
- }
- nrcv+=nlrcv[iarray];
- }
-
- /* Calculate Number of Receivers */
-/*
- nlrcv=(int *)malloc(Nx1*sizeof(int));
- if (!isnan(*xrcv1)) *nlrcv=MIN(NINT((*xrcv2-*xrcv1)/(*dxr)),NINT((*zrcv2-*zrcv1)/(*dzr)))+1;
- else *nlrcv=0;
- nrcv=*nlrcv;
- if (verbose>4 && nlrcv[iarray]!=0) vmess("Linear receiver array 1 has final bounds: (X: %f -> %f,Z: %f ->
-%f)",xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*(*dxr),zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*(*dzr));
- if (Ndx>1) {
- for (iarray=1;iarray<Nx1;iarray++) {
- if (!isnan(xrcv1[iarray])) {
- nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/dxr[iarray]),NINT((zrcv2[iarray]-zrcv1[iarray])/dzr[iarray]))+1;
- }
- else {
- nlrcv[iarray]=0;
- }
- nrcv+=nlrcv[iarray];
- if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %d has final bounds: (X: %f -> %f,Z: %f ->
-%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*dxr[iarray],zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*dzr[iarray]);
- }
- }
- else {
- for (iarray=1;iarray<Nx1;iarray++) {
- if (!isnan(xrcv1[iarray])) nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/(*dxr)),NINT((zrcv2[iarray]-zrcv1[iarray])/(*dzr)))+1;
- else nlrcv[iarray]=0;
- nrcv+=nlrcv[iarray];
- if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %d has final bounds: (X: %f -> %f,Z: %f ->
-%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]**dxr,zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]**dzr);
- }
- }
-*/
- if (verbose) vmess("Total number of linear array receivers: %d",nrcv);
- if (!nrcv) {
- free(xrcv1);
- free(xrcv2);
- free(zrcv1);
- free(zrcv2);
- free(dxr);
- free(dzr);
- free(nlrcv);
- }
- rec->max_nrec+=nrcv;
- }
- else {
- nrcv=0;
- }
-
-/* allocate the receiver arrays */
-
- /* Total Number of Receivers */
- if (verbose) vmess("Total number of receivers: %d",rec->max_nrec);
-
- /* Allocate Arrays */
- rec->x = (int *)calloc(rec->max_nrec,sizeof(int));
- rec->z = (int *)calloc(rec->max_nrec,sizeof(int));
- rec->xr = (float *)calloc(rec->max_nrec,sizeof(float));
- rec->zr = (float *)calloc(rec->max_nrec,sizeof(float));
-
-/* read in the receiver postions */
-
- nrec=0;
- /* Receivers on a Circle */
- if (ncrcv) {
- if (!getparfloat("oxrcv",&oxrcv)) oxrcv=0.0;
- if (!getparfloat("ozrcv",&ozrcv)) ozrcv=0.0;
- if (!getparfloat("arcv",&arcv)) {
- arcv=rrcv;
- for (ix=0; ix<ncrcv; ix++) {
- rec->xr[ix] = oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI));
- rec->zr[ix] = ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI));
- rec->x[ix] = NINT(rec->xr[ix]/dx);
- rec->z[ix] = NINT(rec->zr[ix]/dz);
- //rec->x[ix] = NINT((oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI)))/dx);
- //rec->z[ix] = NINT((ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI)))/dz);
- if (verbose>4) fprintf(stderr,"Receiver Circle: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
- }
- }
- else { /* an ellipse */
- /* simple numerical solution to find equidistant points on an ellipse */
- nh = (ncrcv)*1000; /* should be fine enough for most configurations */
- h = 2.0*M_PI/nh;
- a = MAX(rrcv, arcv);
- b = MIN(rrcv, arcv);
- e = sqrt(a*a-b*b)/a;
- //fprintf(stderr,"a=%f b=%f e=%f\n", a, b, e);
- circ = 0.0;
- for (ir=0; ir<nh; ir++) {
- s = sin(ir*h);
- circ += sqrt(1.0-e*e*s*s);
- }
- circ = a*h*circ;
- //fprintf(stderr,"circ = %f circle=%f\n", circ, 2.0*M_PI*rrcv);
- /* define distance between receivers on ellipse */
- dr = circ/ncrcv;
- ix = 0;
- srun = 0.0;
- if (arcv >= rrcv) phase=0.0;
- else phase=0.5*M_PI;
- for (ir=0; ir<nh; ir++) {
- s = sin(ir*h);
- srun += sqrt(1.0-e*e*s*s);
- if (a*h*srun >= ix*dr ) {
- xr = rrcv*cos(ir*h+phase);
- zr = arcv*sin(ir*h+phase);
- rec->xr[ix] = oxrcv-sub_x0+xr;
- rec->zr[ix] = ozrcv-sub_z0+zr;
- rec->x[ix] = NINT(rec->xr[ix]/dx);
- rec->z[ix] = NINT(rec->zr[ix]/dz);
- if (verbose>4) fprintf(stderr,"Receiver Ellipse: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
- ix++;
- }
- if (ix == ncrcv) break;
- }
- }
-
- /* check if receivers fit into the model otherwise clip to edges */
- for (ix=0; ix<ncrcv; ix++) {
- rec->x[ix] = MIN(nx-1, MAX(rec->x[ix], 0));
- rec->z[ix] = MIN(nz-1, MAX(rec->z[ix], 0));
- }
- nrec += ncrcv;
- }
-
- /* Receiver Text File */
-
- if (ntrcv) {
- /* Allocate arrays */
- xrcva = (float *)malloc(nrcv*sizeof(float));
- zrcva = (float *)malloc(nrcv*sizeof(float));
- /* Read in receiver coordinates */
- for (i=0;i<nrcv;i++) {
- if (fscanf(fp,"%e %e\n",&xrcva[i],&zrcva[i])!=2) vmess("Receiver Text File: Can not parse coordinates on line %d.",i);
- }
- /* Close file */
- fclose(fp);
- /* Process coordinates */
- for (ix=0; ix<nrcv; ix++) {
- rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
- rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
- rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
- rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
- if (verbose>4) vmess("Receiver Text Array: xrcv[%d]=%f zrcv=%f", ix, rec->xr[nrec+ix]+sub_x0, rec->zr[nrec+ix]+sub_z0);
- }
- free(xrcva);
- free(zrcva);
- nrec += ntrcv;
- }
-
- /* Receiver Array */
- if (nxrcv != 0) {
- /* receiver array is defined */
- xrcva = (float *)malloc(nxrcv*sizeof(float));
- zrcva = (float *)malloc(nxrcv*sizeof(float));
- getparfloat("xrcva", xrcva);
- getparfloat("zrcva", zrcva);
- for (ix=0; ix<nxrcv; ix++) {
- rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
- rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
- rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
- rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
- if (verbose>4) fprintf(stderr,"Receiver Array: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[nrec+ix]+sub_x0, rec->zr[nrec+ix]+sub_z0);
- }
- free(xrcva);
- free(zrcva);
- nrec += nxrcv;
- }
-
- /* Linear Receiver Arrays */
- if (nrcv!=0) {
- xrcv1 = (float *)malloc(Nx1*sizeof(float));
- xrcv2 = (float *)malloc(Nx1*sizeof(float));
- zrcv1 = (float *)malloc(Nx1*sizeof(float));
- zrcv2 = (float *)malloc(Nx1*sizeof(float));
-
- if(!getparfloat("xrcv1", xrcv1)) xrcv1[0]=sub_x0;
- if(!getparfloat("xrcv2", xrcv2)) xrcv2[0]=(nx-1)*dx+sub_x0;
- if(!getparfloat("zrcv1", zrcv1)) zrcv1[0]=sub_z0;
- if(!getparfloat("zrcv2", zrcv2)) zrcv2[0]=zrcv1[0];
-
- Ndx = countparval("dxrcv");
- Ndz = countparval("dzrcv");
-
- dxr = (float *)malloc(Nx1*sizeof(float));
- dzr = (float *)malloc(Nx1*sizeof(float));
- if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
- if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
- if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
- for (i=1; i<Nx1; i++) {
- dxr[i] = dxr[0];
- dzr[i] = dzr[0];
- }
- Ndx=1;
- }
- else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
- for (i=1; i<Nx1; i++) {
- dxr[i] = dxr[0];
- dzr[i] = dzr[0];
- }
- Ndz=1;
- }
- else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
- if (Ndx>1) assert(Ndx==Nx1);
- if (Ndz>1) assert(Ndz==Nx1);
- }
-
-/*
- if ( (Ndx!=0) && (Ndz!=0) ) {
- vwarn("Both dzrcv and dxrcv are set: dxrcv value is used");
- Ndz=0;
- for (i=0; i<Nx1; i++) dzr[i] = 0.0;
- }
-*/
- /* check if receiver arrays fit into model */
- for (iarray=0; iarray<Nx1; iarray++) {
- xrcv1[iarray] = MAX(sub_x0, xrcv1[iarray]);
- xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
- xrcv2[iarray] = MAX(sub_x0, xrcv2[iarray]);
- xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
-
- zrcv1[iarray] = MAX(sub_z0, zrcv1[iarray]);
- zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
- zrcv2[iarray] = MAX(sub_z0, zrcv2[iarray]);
- zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
- }
-
- /* calculate receiver array and store into rec->x,z */
-
- for (iarray=0; iarray<Nx1; iarray++) {
- xrange = (xrcv2[iarray]-xrcv1[iarray]);
- zrange = (zrcv2[iarray]-zrcv1[iarray]);
- if (dxr[iarray] != 0.0) {
- nrcv = nlrcv[iarray];
- dxrcv=dxr[iarray];
- dzrcv = zrange/(nrcv-1);
- if (dzrcv != dzr[iarray]) {
- vwarn("For receiver array %d: calculated dzrcv=%f given=%f", iarray, dzrcv, dzr[iarray]);
- vwarn("The calculated receiver distance %f is used", dzrcv);
- }
- }
- else {
- if (dzr[iarray] == 0) {
- verr("For receiver array %d: receiver distance dzrcv is not given", iarray);
- }
- nrcv = nlrcv[iarray];
- dxrcv = xrange/(nrcv-1);
- dzrcv = dzr[iarray];
- if (dxrcv != dxr[iarray]) {
- vwarn("For receiver array %d: calculated dxrcv=%f given=%f", iarray, dxrcv, dxr[iarray]);
- vwarn("The calculated receiver distance %f is used", dxrcv);
- }
- }
-
- // calculate coordinates
- for (ir=0; ir<nrcv; ir++) {
- rec->xr[nrec]=xrcv1[iarray]-sub_x0+ir*dxrcv;
- rec->zr[nrec]=zrcv1[iarray]-sub_z0+ir*dzrcv;
-
- rec->x[nrec]=NINT((rec->xr[nrec])/dx);
- rec->z[nrec]=NINT((rec->zr[nrec])/dz);
- nrec++;
- }
- }
- free(xrcv1);
- free(xrcv2);
- free(zrcv1);
- free(zrcv2);
- free(dxr);
- free(dzr);
- free(nlrcv);
- }
-
- rec->n=rec->max_nrec;
- return 0;
-}
diff --git a/raytime3d/src3d.c b/raytime3d/src3d.c
index 1dadfa5746b00a8bb5504f16de26188c06da5c97..545c1974d14ab60f4bd28b259e68d2d4b4ee0e21 100644
--- a/raytime3d/src3d.c
+++ b/raytime3d/src3d.c
@@ -1,5 +1,8 @@
-#include <cwp.h>
-#include <su.h>
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
#include <fcntl.h>
#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
@@ -7,10 +10,10 @@
#define SQR(x) ((x) * (x))
#define DIST(x,y,z,x1,y1,z1) sqrt(SQR(x-(x1))+SQR(y-(y1)) + SQR(z-(z1)))
-/* definitions from saerrpkge.c */
-extern void saerr(char *fmt, ...);
-extern void sawarn(char *fmt, ...);
-extern void samess(char *fmt, ...);
+/* definitions from verbose.c */
+extern void verr(char *fmt, ...);
+extern void vwarn(char *fmt, ...);
+extern void vmess(char *fmt, ...);
void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube)
{
@@ -48,9 +51,9 @@ void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox
if(!getparint("srctype",&srctype)) srctype=1;
if(srctype==1) {
- if(!getparfloat("xsrc1",&xsrc1)) saerr("xsrc1 not given");
- if(!getparfloat("ysrc1",&ysrc1)) saerr("ysrc1 not given");
- if(!getparfloat("zsrc1",&zsrc1)) saerr("zsrc1 not given");
+ if(!getparfloat("xsrc1",&xsrc1)) verr("xsrc1 not given");
+ if(!getparfloat("ysrc1",&ysrc1)) verr("ysrc1 not given");
+ if(!getparfloat("zsrc1",&zsrc1)) verr("zsrc1 not given");
fxs = (xsrc1-ox)/h;
fys = (ysrc1-oy)/h;
fzs = (zsrc1-oz)/h;
@@ -58,30 +61,30 @@ void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox
ys = (int)(fys + 0.5);
zs = (int)(fzs + 0.5);
if(xs<2 || ys<2 || zs<2 || xs>nx-3 || ys>ny-3 || zs>nz-3){
- sawarn("Source near an edge, beware of traveltime errors");
- sawarn("for raypaths that travel parallel to edge ");
- sawarn("while wavefronts are strongly curved, (JV, 8/17/88)\n");
+ vwarn("Source near an edge, beware of traveltime errors");
+ vwarn("for raypaths that travel parallel to edge ");
+ vwarn("while wavefronts are strongly curved, (JV, 8/17/88)\n");
}
*pxs = xs; *pys = ys, *pzs = zs, *cube = NCUBE;
}
else if (srctype==2) {
- if (!getparint("srcwall",&srcwall)) saerr("srcwall not given");
- if (!getparstring("wallfile",&wallfile)) saerr("wallfile not given");
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
+ if (!getparstring("wallfile",&wallfile)) verr("wallfile not given");
if((wfint=open(wallfile,O_RDONLY,0664))<=1) {
fprintf(stderr,"cannot open %s\n",wallfile);
exit(-1);
}
}
else if (srctype==3) {
- if (!getparint("srcwall",&srcwall)) saerr("srcwall not given");
- if (!getparstring("oldtfile",&oldtfile)) saerr("oldtfile not given");
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
+ if (!getparstring("oldtfile",&oldtfile)) verr("oldtfile not given");
if((ofint=open(oldtfile,O_RDONLY,0664))<=1) {
fprintf(stderr,"cannot open %s\n",oldtfile);
exit(-1);
}
}
else {
- saerr("ERROR: incorrect value of srctype");
+ verr("ERROR: incorrect value of srctype");
}
nxy = nx * ny;
@@ -127,8 +130,8 @@ void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox
/* can't handle asin(1.) ! */
if (fabs(rz1-rzc)>=rho) rho=1.0000001*fabs(rz1-rzc);
theta2 = asin((rz1-rzc)/rho);
- if (rxyc<0) theta1=PI-theta1;
- if (rxyc<rxy1) theta2=PI-theta2;
+ if (rxyc<0) theta1=M_PI-theta1;
+ if (rxyc<rxy1) theta2=M_PI-theta2;
t0(xx,yy,zz) = log(tan(theta2/2)/tan(theta1/2)) / dv;
}
}
diff --git a/raytime3d/vidale3d.c b/raytime3d/vidale3d.c
index 06b4b388274def5752a95ab1f0bbb1b881b47c11..ff8a74c1c4f75f4179ce168c514b6ee1949ab837 100644
--- a/raytime3d/vidale3d.c
+++ b/raytime3d/vidale3d.c
@@ -1,5 +1,9 @@
-#include <cwp.h>
-#include <su.h>
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
+#include"par.h"
#define SQR2 1.414213562
#define SQR3 1.732050808
@@ -7,10 +11,10 @@
#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
#define s0(x,y,z) slow0[nxy*(z) + nx*(y) + (x)]
-/* definitions from saerrpkge.c */
-extern void saerr(char *fmt, ...);
-extern void sawarn(char *fmt, ...);
-extern void samess(char *fmt, ...);
+/* definitions from verbose.c */
+extern void verr(char *fmt, ...);
+extern void vwarn(char *fmt, ...);
+extern void vmess(char *fmt, ...);
struct sorted
{ float time; int i1, i2;};
@@ -82,7 +86,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
/* SET MAXIMUM RADIUS TO COMPUTE */
if (maxoff > 0.) {
maxrad = maxoff/h + 1;
- sawarn("WARNING: Computing only to max radius = %d",maxrad);
+ vwarn("WARNING: Computing only to max radius = %d",maxrad);
}
else maxrad = 99999999;
@@ -106,7 +110,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
}
wall = (float *) malloc(4*nwall);
if(sort == NULL || wall == NULL)
- saerr("error in allocation of arrays sort and wall");
+ verr("error in allocation of arrays sort and wall");
if(!getparint("srctype",&srctype)) srctype=1;
if(srctype==1) {
@@ -126,7 +130,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
else{ z2 = nz; dz2 = 0;}
}
else {
- if (!getparint("srcwall",&srcwall)) saerr("srcwall not given");
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
/* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE */
radius = 1;
if (srcwall == 1) x2=1;
@@ -1245,7 +1249,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
/* UPDATE RADIUS */
radius++;
- if(radius%10 == 0 && verbose) samess("Completed radius = %d",radius);
+ if(radius%10 == 0 && verbose) vmess("Completed radius = %d",radius);
if(radius == maxrad) rad0 = 0;
} /* END BIG LOOP */
@@ -1259,42 +1263,42 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
else {
head=0;
if (headw[1]>0) {
- if(verbose) samess("Head waves found on left: %d",headw[1]);
+ if(verbose) vmess("Head waves found on left: %d",headw[1]);
if (headw[1]>head) {
head = headw[1];
srcwall = 1;
}
}
if (headw[2]>0) {
- if(verbose) samess("Head waves found on right: %d",headw[2]);
+ if(verbose) vmess("Head waves found on right: %d",headw[2]);
if (headw[2]>head) {
head = headw[2];
srcwall = 2;
}
}
if (headw[3]>0) {
- if(verbose) samess("Head waves found on front: %d",headw[3]);
+ if(verbose) vmess("Head waves found on front: %d",headw[3]);
if (headw[3]>head) {
head = headw[3];
srcwall = 3;
}
}
if (headw[4]>0) {
- if(verbose) samess("Head waves found on back: %d",headw[4]);
+ if(verbose) vmess("Head waves found on back: %d",headw[4]);
if (headw[4]>head) {
head = headw[4];
srcwall = 4;
}
}
if (headw[5]>0) {
- if(verbose) samess("Head waves found on top: %d",headw[5]);
+ if(verbose) vmess("Head waves found on top: %d",headw[5]);
if (headw[5]>head) {
head = headw[5];
srcwall = 5;
}
}
if (headw[6]>0) {
- if(verbose) samess("Head waves found on bottom: %d",headw[6]);
+ if(verbose) vmess("Head waves found on bottom: %d",headw[6]);
if (headw[6]>head) {
head = headw[6];
srcwall = 6;
@@ -1302,32 +1306,32 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
}
if (headpref>0 && headw[headpref]>0) {
if(verbose)
- samess("Preference to restart on wall opposite source");
+ vmess("Preference to restart on wall opposite source");
srcwall = headpref;
}
/* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE */
dx1=1; dx2=1; dy1=1; dy2=1; dz1=1; dz2=1; rad0=1;
radius = 1;
if (srcwall == 1) { x2=1;
- samess("RESTART at left side of model"); }
+ vmess("RESTART at left side of model"); }
else { x2=nx; dx2=0; }
if (srcwall == 2) { x1=nx-2;
- samess("RESTART at right side of model"); }
+ vmess("RESTART at right side of model"); }
else { x1= -1; dx1=0; }
if (srcwall == 3) { y2=1;
- samess("RESTART at front side of model"); }
+ vmess("RESTART at front side of model"); }
else { y2=ny; dy2=0; }
if (srcwall == 4) { y1=ny-2;
- samess("RESTART at back side of model"); }
+ vmess("RESTART at back side of model"); }
else { y1= -1; dy1=0; }
if (srcwall == 5) { z2=1;
- samess("RESTART at top side of model"); }
+ vmess("RESTART at top side of model"); }
else { z2=nz; dz2=0; }
if (srcwall == 6) { z1=nz-2;
- samess("RESTART at bottom side of model"); }
+ vmess("RESTART at bottom side of model"); }
else { z1= -1; dz1=0; }
if (reverse == 0)
- sawarn("RESTART CANCELLED by choice of input parameter `reverse`");
+ vwarn("RESTART CANCELLED by choice of input parameter `reverse`");
}
reverse--;
diff --git a/raytime3d/writeSrcRecPos.c b/raytime3d/writeSrcRecPos.c
index faf297811e469314c8de0265f78b87543b27af9b..d9d03da16dae329b7628e3420224bda8d9a6cce4 100644
--- a/raytime3d/writeSrcRecPos.c
+++ b/raytime3d/writeSrcRecPos.c
@@ -3,7 +3,7 @@
#include<math.h>
#include<assert.h>
#include"par.h"
-#include"raytime.h"
+#include"raytime3d.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
diff --git a/raytime3d/writesufile.c b/raytime3d/writesufile.c
index 6eac57d300aaa653e621b99fe2731c9ed3ddb60a..7ebcbbd198858193b67a2b9c2fbb5670feb8b5d8 100644
--- a/raytime3d/writesufile.c
+++ b/raytime3d/writesufile.c
@@ -3,7 +3,7 @@
#include <assert.h>
#include <string.h>
#include "par.h"
-#include "raytime.h"
+#include "raytime3d.h"
#include "SUsegy.h"
#include "segy.h"
diff --git a/utils/Makefile b/utils/Makefile
index 3b59b4f71f844acc56c9b9851deba24861d82e5b..9d47eceee936c6e5dffd52824bbd667ed5a4adea 100644
--- a/utils/Makefile
+++ b/utils/Makefile
@@ -6,7 +6,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC += -openmp
#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 = \
makemod.c \
@@ -84,6 +84,17 @@ SRCG = green.c \
docpkge.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 \
getFileInfo.c \
kxwfilter.c \
@@ -152,6 +163,11 @@ OBJG = $(SRCG:%.c=%.o)
green: $(OBJG)
$(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)
basop: $(OBJB)
@@ -172,23 +188,24 @@ OBJT = $(SRCT:%.c=%.o)
ftr1d: $(OBJT)
$(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 makewave $B
cp extendModel $B
cp fconv $B
cp correigen $B
cp green $B
+ cp green3D $B
cp basop $B
cp syn2d $B
cp mat2su $B
cp ftr1d $B
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
- 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
diff --git a/utils/getrecpos3D.c b/utils/getrecpos3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..85e9623fbc57864eb7be454d66520faf5e88c28b
--- /dev/null
+++ b/utils/getrecpos3D.c
@@ -0,0 +1,135 @@
+#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;
+}
diff --git a/utils/green3D b/utils/green3D
new file mode 100755
index 0000000000000000000000000000000000000000..cc36382df2db967794d33a3b4afe6e4bb4e56e07
Binary files /dev/null and b/utils/green3D differ
diff --git a/utils/green3D.c b/utils/green3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..69daa8ea3552b7083eb3871d87ac1ee87694fe0a
--- /dev/null
+++ b/utils/green3D.c
@@ -0,0 +1,761 @@
+#include <genfft.h>
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+
+#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 COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+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);
+
+void xwgreen3D(float *data, int nt, int nx, int ny, float dt, float fmin, float fmax, float *xi, float xsrc,
+ float dx, float *yi, float ysrc, float dy, float *zi, float zsrc, float c, float cs, float rho,
+ float *wavelet, float dipx, float maxdip, int far, int p_vz, int dip, int verbose);
+
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" green - calculation of 2D Greens function in homogenoeus medium based one exact expressions",
+" ",
+" green c= zsrc1= [optional parameters]",
+" ",
+" Required parameters:",
+" ",
+" c= ....................... P-wave velocity",
+" cs=0.7*c ................. S-wave velocity",
+" zsrc1= ................... depth of source",
+" ",
+" Optional parameters:",
+" ",
+" file_out= ................ output file (default SU-pipe)",
+" RECEIVER POSITIONS ",
+" xrcv=-1500,1500 .......... x-position's of receivers (array)",
+" yrcv=-1500,1500 .......... y-position's of receivers (array)",
+" zrcv=0,0 ................. z-position's of receivers (array)",
+" dxrcv=15 ................. step in receiver x-direction",
+" dyrcv=15 ................. step in receiver y-direction",
+" var=0 .................... variance for irregular sampling (dxrcv +- var)",
+" seed=0 ................... seed for random generator",
+" lint=1 ................... linear interpolate between the rcv points",
+" rrcv= .................... radius for receivers on a circle ",
+" oxrcv=0.0 ................ x-center position of circle",
+" oyrcv=0.0 ................ y-center position of circle",
+" ozrcv=0.0 ................ z-center position of circle",
+" dphi=2 ................... angle between receivers on circle ",
+" SOURCE POSITIONS ",
+" xsrc1=0.0 ................ x-position of first source",
+" xsrc2=xsrc1 .............. x-position of last source",
+" dxsrc=0.0 ................ step in source x-direction",
+" ysrc1=0.0 ................ y-position of first source",
+" ysrc2=ysrc1 .............. y-position of last source",
+" dysrc=0.0 ................ step in source y-direction",
+" zsrc2=zsrc1 .............. depth position (z) of last source",
+" dzsrc=0.0 ................ step in source z-direction",
+" SAMPLING AND SOURCE DEFINITION ",
+" file_src=spike ........... source wavelet (overrules dt)",
+" nt=256 ................... number of samples",
+" dt=0.004 ................. stepsize in time-direction ",
+" fmin=0 ................... minimum frequency",
+" fmax=70 .................. maximum frequency",
+" dipx=0 ................... local dip of the dipole in x-direction",
+" dipy=0 ................... local dip of the dipole in y-direction",
+" dip=1 .................... 1; dipole 0; monopole source",
+" rho=1000 ................. density",
+" FIELD DEFINITION ",
+" far=0 .................... farfield approximation 0=off)",
+" p_vz=0 .................. P or Vz field (0 = P field, 1 = Vz field)",
+" Fz=0 .................... Force source in z with Vz receivers",
+" Fx=0 .................... Force source in x with Vz receivers",
+" maxdip=90 ................ maximum angle (degrees) to be computed ",
+" sum=0 .................... sum all sources",
+" verbose=0 ................ silent option; >0 display info",
+"",
+" The P or Vz field of a dipole source at depth z below the receivers",
+" in a homogeneous 2-D medium is calculated.",
+" ",
+" author : Jan Thorbecke : 23-03-1995 (janth@xs4all.nl)",
+" : revision 2010",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main(int argc, char **argv)
+{
+ FILE *fp_in, *fp_out;
+ int n1, n2, n3, i, j, l, nrx, nry, nrz, dip;
+ int far, p_vz, nt, nx, ny, Nsx, Nsy, is, isy, sum, lint, verbose;
+ int size, ntraces, ngath, Fz, Fx;
+ float scl, xmin, xmax, ymin, ymax;
+ float dx, dy, dt, d1, d2, d3, fmin, fmax, f1, f2, f3, c, cs, rho;
+ float *data, *wavelet, *tmpdata, dipx, dipy, xsrc1, xsrc2, ysrc1, ysrc2;
+ float *xrcv, *yrcv, *zrcv, *xi, *yi, *zi, x0, y0, maxdip;
+ float rrcv, dphi, oxrcv, ozrcv;
+ float zsrc1, zsrc2, dxsrc, dysrc, dzsrc, xsrc, ysrc, zsrc, dxrcv, dyrcv;
+ char *file_src, *file_out;
+ size_t nwrite;
+ segy *hdrs;
+
+/* ========================= Reading parameters ====================== */
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if(!getparint("verbose", &verbose)) verbose = 0;
+ if(!getparstring("file_out", &file_out)){
+ if (verbose) vwarn("parameter file_out not found, assume pipe");
+ file_out = NULL;
+ }
+ if(!getparstring("file_src", &file_src)) file_src = NULL;
+ if(!getparfloat("c", &c)) verr("velocity must be specified.");
+ if(!getparfloat("cs", &cs)) cs=0.7*c;
+ if(!getparfloat("zsrc1", &zsrc1)) verr("zsrc1(depth) must be specified.");
+ if(!getparint("lint", &lint)) lint=1;
+ if(!getparfloat("maxdip", &maxdip)) maxdip=90.0;
+
+ nrx = countparval("xrcv");
+ nry = countparval("yrcv");
+ // nrz = countparval("zrcv");
+ nrz = 0;
+ if(!getparfloat("dxrcv",&dxrcv)) dxrcv = 15;
+ if(!getparfloat("dyrcv",&dyrcv)) dyrcv = 15;
+
+ if (nrx != 0 && nry != 0 && nrz == 0) {
+ if (nrx != 2) verr("xrcv should have only two values");
+ if (nry != 2) verr("yrcv should have only two values");
+ xrcv = (float *)malloc(nrx*sizeof(float));
+ yrcv = (float *)malloc(nry*sizeof(float));
+ getparfloat("xrcv",xrcv);
+ getparfloat("yrcv",yrcv);
+ nx = NINT((xrcv[1] - xrcv[0])/dxrcv) + 1;
+ ny = NINT((yrcv[1] - yrcv[0])/dyrcv) + 1;
+ xi = (float *)malloc(nx*ny*sizeof(float));
+ yi = (float *)malloc(nx*ny*sizeof(float));
+ zi = (float *)malloc(nx*ny*sizeof(float));
+ x0 = xrcv[0];
+ y0 = yrcv[0];
+ for (i = 0; i < ny; i++) {
+ for (j = 0; j < nx; j++) {
+ xi[i*nx+j] = x0 + j*dxrcv;
+ yi[i*nx+j] = y0 + i*dyrcv;
+ zi[i*nx+j] = 0;
+ }
+ }
+ }
+ else if (nrx == 0 && nry == 0 && nrz == 0) {
+ nx = NINT((3000)/dxrcv) + 1;
+ ny = NINT((3000)/dyrcv) + 1;
+ xi = (float *)malloc(nx*ny*sizeof(float));
+ yi = (float *)malloc(nx*ny*sizeof(float));
+ zi = (float *)malloc(nx*ny*sizeof(float));
+ x0 = -1500;
+ y0 = -1500;
+ for (i = 0; i < ny; i++) {
+ for (j = 0; j < nx; j++) {
+ xi[i*nx+j] = x0 + j*dxrcv;
+ yi[i*nx+j] = y0 + i*dyrcv;
+ zi[i*nx+j] = 0;
+ }
+ }
+ }
+ else verr("Number of xrcv and yrcv values are not equal");
+
+ if (verbose) vmess("number of receivers nx = %d, ny = %d total = %d", nx, ny, nx*ny);
+ if (verbose == 13) {
+ for (j = 0; j < ny; j++) {
+ for (i = 0; i < nx; i++) {
+ vmess("xi = %d yi = %d x = %f y=%f z = %f", i, j, xi[j*nx+i], yi[j*nx+i], zi[j*nx+i]);
+ }
+ }
+ }
+
+ if(!getparfloat("xsrc1", &xsrc1)) xsrc1=0;
+ if(!getparfloat("xsrc2", &xsrc2)) xsrc2=xsrc1;
+ if(!getparfloat("dxsrc", &dxsrc)) dxsrc=0.0;
+ if(!getparfloat("ysrc1", &ysrc1)) ysrc1=0;
+ if(!getparfloat("ysrc2", &ysrc2)) ysrc2=ysrc1;
+ if(!getparfloat("dysrc", &dysrc)) dysrc=0.0;
+ if(!getparfloat("zsrc2", &zsrc2)) zsrc2=zsrc1;
+ if(!getparfloat("dzsrc", &dzsrc)) dzsrc=0;
+ if(!getparint("nt", &nt)) nt = 256;
+ if(!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if(!getparfloat("fmax", &fmax)) fmax = 70.0;
+ if(!getparfloat("dipx", &dipx)) dipx = 0.0;
+ if(!getparfloat("dipy", &dipy)) dipy = 0.0;
+ if(!getparfloat("rho", &rho)) rho = 1000.0;
+ if(!getparint("far", &far)) far = 0;
+ if(!getparint("p_vz", &p_vz)) p_vz = 0;
+ if(!getparint("Fz", &Fz)) Fz = 0;
+ if(!getparint("Fx", &Fx)) Fx = 0;
+ if(!getparint("dip", &dip)) dip = 1;
+ if(!getparint("sum", &sum)) sum = 0;
+ if(Fz) p_vz=2;
+ if(Fx) p_vz=3;
+
+/* ========================= Opening wavelet file ====================== */
+
+ if (file_src == NULL){
+ if(!getparfloat("dt", &dt)) dt = 0.004;
+ wavelet = (float *)calloc(nt,sizeof(float));
+ wavelet[0] = 1.0;
+ }
+ else {
+ if (verbose) vmess("Reading wavelet from file %s.", file_src);
+ ngath = 1;
+ getFileInfo(file_src, &n1, &n2, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &scl, &ntraces);
+
+ fp_in = fopen(file_src, "r");
+ if (fp_in == NULL) verr("error on opening input file_src=%s", file_src);
+
+ tmpdata = (float *)calloc(n1*n2,sizeof(float));
+ hdrs = (segy *) calloc(n2,sizeof(segy));
+
+ n2 = readData(fp_in, tmpdata, hdrs, n1);
+ fclose(fp_in);
+ if (verbose) {
+ disp_fileinfo(file_src, n1, n2, f1, f2, d1, d2, hdrs);
+ }
+ dt = d1;
+ wavelet = (float *)calloc(nt,sizeof(float));
+
+ if (n1 <= nt) {
+ for (i = 0; i < n1; i++) wavelet[i] = tmpdata[i];
+ for (i = n1; i < nt; i++) wavelet[i] = 0.0;
+ }
+ else {
+ vwarn("file_src has more samples than output");
+ for (i = 0; i < nt; i++) wavelet[i] = tmpdata[i];
+ }
+ if( tmpdata ) free(tmpdata);
+ if( hdrs ) free( (void *) hdrs);
+ }
+
+/* ============ INITIALIZE AND CHECK PARAMETERS =============== */
+
+ if (xsrc2==xsrc1) Nsx = 1;
+ else Nsx = NINT((xsrc2 - xsrc1)/dxsrc) + 1;
+ if (ysrc2==ysrc1) Nsy = 1;
+ else Nsy = NINT((ysrc2 - ysrc1)/dysrc) + 1;
+
+ if (verbose) vmess("Number of shot records to generate x = %d y = %d", Nsx, Nsy);
+ if (Nsx > 1 && Nsy > 1) {
+ dxsrc = (xsrc2-xsrc1)/(Nsx-1);
+ dysrc = (ysrc2-ysrc1)/(Nsy-1);
+ dzsrc = (zsrc2-zsrc1)/(Nsx-1);
+ if (verbose) {
+ vmess("dxsrc = %f", dxsrc);
+ vmess("dysrc = %f", dysrc);
+ vmess("dzsrc = %f", dzsrc);
+ }
+ }
+
+ size = nt * nx *ny;
+ dx = dxrcv;
+ dy = dyrcv;
+ tmpdata = (float *)calloc(size,sizeof(float));
+ data = (float *)calloc(size,sizeof(float));
+ hdrs = (segy *) calloc(nx*ny,sizeof(segy));
+ for (i = 0; i < ny; i++) {
+ for(j = 0; j < nx; j++) {
+ hdrs[i*nx+j].f1= 0.0;
+ hdrs[i*nx+j].f2= x0;
+ hdrs[i*nx+j].d1= dt;
+ hdrs[i*nx+j].d2= dx;
+ hdrs[i*nx+j].ns= nt;
+ hdrs[i*nx+j].dt= (int)1000000*dt;
+ hdrs[i*nx+j].trwf= nx*ny;
+ hdrs[i*nx+j].tracl= i*nx+j+1;
+ hdrs[i*nx+j].tracf= i*nx+j+1;
+ hdrs[i*nx+j].gx = (x0 + j*dx)*1000;
+ hdrs[i*nx+j].gy = (y0 + i*dy)*1000;
+ hdrs[i*nx+j].scalco = -1000;
+ hdrs[i*nx+j].trid = TREAL;
+ }
+ }
+ if (file_out==NULL) fp_out=stdout;
+ else fp_out = fopen(file_out,"w");
+ if (fp_out == NULL) verr("error in creating output file");
+
+ for (isy = 0; isy < Nsy; isy++) {
+ for (is = 0; is < Nsx; is++) {
+ xsrc = xsrc1 + is*dxsrc;
+ ysrc = ysrc1 + isy*dysrc;
+ zsrc = zsrc1 + is*dzsrc;
+ if (verbose) vmess("xsrc = %f ysrc=%f zsrc = %f", xsrc, ysrc, zsrc);
+
+ xwgreen3D(data,nt,nx,ny,dt,fmin,fmax,xi,xsrc,dx,yi,ysrc,dy,zi,zsrc,c,cs,rho,wavelet,
+ dipx, maxdip, far, p_vz, dip, verbose);
+
+ for (l = 0; l < ny; l++) {
+ for (i = 0; i < nx; i++) {
+ for (j = 0; j < nt; j++) tmpdata[l*nx*nt+i*nt+j] = data[l*nx*nt+i*nt+j];
+ hdrs[l*nx+i].sx = NINT(xsrc*1000);
+ hdrs[l*nx+i].sy = NINT(ysrc*1000);
+ hdrs[l*nx+i].scalco = -1000;
+ hdrs[l*nx+i].offset = xi[l*nx+i]-xsrc;
+ hdrs[l*nx+i].gx = NINT(xi[l*nx+i]*1000);
+ hdrs[l*nx+i].gy = NINT(yi[l*nx+i]*1000);
+ hdrs[l*nx+i].fldr = isy*Nsx+is+1;
+ hdrs[l*nx+i].trwf = nx*ny;
+ nwrite = fwrite( &hdrs[l*nx+i], 1, TRCBYTES, fp_out);
+ assert(nwrite == TRCBYTES);
+ nwrite = fwrite( &tmpdata[l*nx*nt+i*nt], sizeof(float), nt, fp_out);
+ assert(nwrite == nt);
+ }
+ }
+ }
+ }
+
+ if( xi ) free(xi);
+ if( yi ) free(yi);
+ if( zi ) free(zi);
+ if( wavelet ) free( wavelet );
+
+ fclose(fp_out);
+
+ if( data ) free(data);
+ if( tmpdata ) free(tmpdata);
+ if( hdrs ) free( hdrs);
+
+ exit ( 0 );
+}
+
+/***************************************************************************
+*
+* Calculation of pulse response in homogeneous medium
+*
+*
+***************************************************************************/
+
+void xwgreen3D(float *data, int nt, int nx, int ny, float dt, float fmin, float fmax, float *xi, float xsrc, float dx, float *yi, float ysrc, float dy, float *zi, float zsrc, float c, float cs, float rho, float *wavelet, float dipx, float maxdip, int far, int p_vz, int dip, int verbose)
+{
+ int iomin, iomax, iom, ix, iy, nfreq, i, sign, optn;
+ float df, deltom, om, k, r, x, y, invr, phi, phi2, cosphi;
+ float *rwave, *rdata, cos2, scl, z, kp, ks, sclr;
+ complex *cwave, *cdata, tmp, tmp2, sum;
+ complex H02p, H12p, H02s, H12s, Gp, Gs;
+
+ optn = optncr(nt);
+ nfreq = 1+(optn/2);
+ df = 1.0/(dt*optn);
+ deltom = 2.*M_PI*df;
+ iomin = (int)MIN((fmin*dt*optn), (nfreq-1));
+ iomin = MAX(iomin, 1);
+ iomax = MIN((int)(fmax*dt*optn), (nfreq-1));
+
+ rdata = (float *)calloc(optn*nx*ny,sizeof(float));
+ cdata = (complex *)calloc(nfreq*nx*ny,sizeof(complex));
+ rwave = (float *)calloc(optn,sizeof(float));
+ cwave = (complex *)calloc(nfreq,sizeof(complex));
+
+ for (i = 0; i < nt; i++) rwave[i] = wavelet[i]*dt;
+ for (i = nt; i < optn; i++) rwave[i] = 0.0;
+
+ sign = -1;
+ rc1fft(rwave, cwave, optn, sign);
+
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+
+ if (p_vz == 0) {
+ if (far == 0 && dip == 1) {
+ if (verbose) vmess("near and far P field of dipole");
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ phi2 = SGN(x)*phi - (dipx*M_PI/180.0);
+ cosphi = 0.25*cos(phi2)*rho;
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ k = om/c;
+ tmp.r = -k*cosphi*y1(k*r);
+ tmp.i = -k*cosphi*j1(k*r);
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+ }
+ else if (far == 1 && dip == 1){
+ if (verbose) vmess("far P field of dipole");
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ phi2 = SGN(x)*phi - (dipx*M_PI/180.0);
+ cosphi = 0.5*cos(phi2)*rho/sqrt(r);
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ k = om/c;
+ tmp.r = sqrt(k/(2.0*M_PI))*cosphi*cos(k*r-M_PI/4.0);
+ tmp.i = -sqrt(k/(2.0*M_PI))*cosphi*sin(k*r-M_PI/4.0);
+
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+ }
+ else if (far == 0 && dip == 0){
+ if (verbose) vmess("near and far P field of monopole");
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ scl = 0.25*rho;
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ k = om/c;
+ tmp.r = -scl*y0(k*r);
+ tmp.i = -scl*j0(k*r);
+
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+ }
+ else if (far == 1 && dip == 0){
+ if (verbose) vmess("far P field of monopole");
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI*0.5;
+ scl = 0.5*rho/sqrt(r);
+ if (fabs(phi) <= M_PI*(maxdip/180.0)) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ k = om/c;
+ tmp.r = -sqrt(1.0/(2.0*M_PI*k))*scl*sin(k*r-M_PI/4.0);
+ tmp.i = -sqrt(1.0/(2.0*M_PI*k))*scl*cos(k*r-M_PI/4.0);
+
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+ }
+ }
+ else if (p_vz == 1) {
+ if (dip == 1) {
+ if (far == 0) {
+ if (verbose) vmess("near and far Vz field of dipole");
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+ invr = -0.25/(c);
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ phi2 = SGN(x)*phi - (dipx*M_PI/180.0);
+ cosphi = cos(phi2);
+ cos2 = cosphi*cosphi;
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ k = om/c;
+ tmp.r = k*cos2*invr*j0(k*r);
+ tmp.i = -k*cos2*invr*y0(k*r);
+ tmp2.r = k*(1-2*cos2)*invr*j1(k*r)/(k*r);
+ tmp2.i = -k*(1-2*cos2)*invr*y1(k*r)/(k*r);
+ sum.r = tmp.r + tmp2.r;
+ sum.i = tmp.i + tmp2.i;
+
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = sum.r*cwave[iom].r -
+ sum.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = sum.r*cwave[iom].i +
+ sum.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+ }
+ else {
+ if (verbose) vmess("far Vz field of dipole");
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+ invr = -0.5/(c*sqrt(r));
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ phi2 = SGN(x)*phi - (dipx*M_PI/180.0);
+ cosphi = cos(phi2);
+ cos2 = cosphi*cosphi;
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ k = om/c;
+ tmp.r = sqrt(k/(2.0*M_PI))*invr*cos2*cos(k*r-M_PI/4.0);
+ tmp.i = -sqrt(k/(2.0*M_PI))*invr*cos2*sin(k*r-M_PI/4.0);
+
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+ }
+ }
+ else {
+ if (verbose) vmess("near and far Vz field of monopole");
+
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ phi2 = SGN(x)*phi - (dipx*M_PI/180.0);
+ cosphi = cos(phi2);
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ k = om/c;
+ tmp.i = -cosphi*y1(k*r)/(4.0*c);
+ tmp.r = cosphi*j1(k*r)/(4.0*c);
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+ }
+ }
+ else if (p_vz == 2) { /* Fz source with Vz receivers Fz=1 == p_vz=2 */
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ phi2 = SGN(x)*phi - (dipx*M_PI/180.0);
+ cosphi = cos(phi2);
+ sclr = (z*z-x*x-y*y)/(r);
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ kp = om/c;
+ ks = om/cs;
+ H02p.r = j0(kp*r);
+ H02p.i = -y0(kp*r);
+ H12p.r = j1(kp*r);
+ H12p.i = -y1(kp*r);
+
+ H02s.r = j0(ks*r);
+ H02s.i = -y0(ks*r);
+ H12s.r = j1(ks*r);
+ H12s.i = -y1(ks*r);
+
+ Gp.r = kp/(4*om*rho*r*r)*(-z*z*kp*H02p.r + sclr*H12p.r);
+ Gp.i = kp/(4*om*rho*r*r)*(-z*z*kp*H02p.i + sclr*H12p.i);
+
+ Gs.r = ks/(4*om*rho*r*r)*(-z*z*ks*H02s.r + sclr*H12s.r);
+ Gs.i = ks/(4*om*rho*r*r)*(-z*z*ks*H02s.i + sclr*H12s.i);
+
+ tmp.i = (-1.0/om)*((om/(4*rho*cs*cs))*(H02s.r) - Gp.r + Gs.r);
+ tmp.r = ( 1.0/om)*((om/(4*rho*cs*cs))*(H02s.i) - Gp.i + Gs.i);
+
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+
+ }
+ else if (p_vz == 3) { /* Fx source with Vz receivers Fx=1 == p_vz=3 */
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ x = xi[iy*nx+ix] - xsrc;
+ y = yi[iy*nx+ix] - ysrc;
+ z = fabs(zi[iy*nx+ix] - zsrc);
+ r = sqrt(x*x + y*y + z*z);
+
+ if (r != 0) phi = acos(z/r);
+ else phi = M_PI/2;
+ phi2 = SGN(x)*phi - (dipx*M_PI/180.0);
+ cosphi = cos(phi2);
+ scl = (z*x*y)/(4.0*r*r*rho);
+ if (fabs(phi) < maxdip*M_PI/180.0) {
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ kp = om/c;
+ ks = om/cs;
+ H02p.r = kp*kp*j0(kp*r);
+ H02p.i = -kp*kp*y0(kp*r);
+ H12p.r = 2.0*kp*j1(kp*r)/r;
+ H12p.i = -2.0*kp*y1(kp*r)/r;
+
+ H02s.r = ks*ks*j0(ks*r);
+ H02s.i = -ks*ks*y0(ks*r);
+ H12s.r = 2.0*ks*j1(ks*r)/r;
+ H12s.i = -2.0*ks*y1(ks*r)/r;
+
+ tmp.i = (scl/(om*om))*((H02p.r-H12p.r) - (H02s.r-H12s.r));
+ tmp.r = -(scl/(om*om))*((H02p.i-H12p.i) - (H02s.i-H12s.i));
+
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = tmp.r*cwave[iom].r -
+ tmp.i*cwave[iom].i;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = tmp.r*cwave[iom].i +
+ tmp.i*cwave[iom].r;
+ }
+ }
+ else {
+ for (iom = iomin; iom <= iomax; iom++) {
+ cdata[iy*nx*nfreq+ix*nfreq+iom].r = 0.0;
+ cdata[iy*nx*nfreq+ix*nfreq+iom].i = 0.0;
+ }
+ }
+ }
+ }
+
+ }
+
+
+ scl = df;
+ sign = 1;
+ crmfft(&cdata[0], &rdata[0], optn, nx*ny, nfreq, optn, sign);
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (i = 0; i < nt; i++) {
+ data[iy*nx*nt+ix*nt+i] = scl*rdata[iy*nx*optn+ix*optn+i];
+ }
+ }
+ }
+
+ free(cdata);
+ free(cwave);
+ free(rdata);
+ free(rwave);
+
+ return;
+}