diff --git a/fdelmodc/.DS_Store b/fdelmodc/.DS_Store
index c0cb59e77110d2339cdcace7ece85b15b8ff66a3..f32ee3de07f0d73c742fdadceb960e8afa6cc66c 100644
Binary files a/fdelmodc/.DS_Store and b/fdelmodc/.DS_Store differ
diff --git a/fdelmodc/demo/demoStaircase.scr b/fdelmodc/demo/demoStaircase.scr
new file mode 100755
index 0000000000000000000000000000000000000000..4fd328a3cc41f996de5639e46763e2894db0341b
--- /dev/null
+++ b/fdelmodc/demo/demoStaircase.scr
@@ -0,0 +1,85 @@
+#!/bin/bash
+
+#export PATH=$HOME/bin:$HOME/src/OpenSource/utils:$PATH:
+
+cp=2000
+rho=2500
+dx=2.5
+dt=0.0005
+
+#dx=5
+#dt=0.0005
+
+makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
+ orig=-2500,0 file_base=tilt.su \
+ intt=def x=-2500,0,2500 z=2000,1200,400 cp=2500 ro=5500
+
+makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
+ orig=-2500,0 file_base=hom.su
+
+makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
+
+export OMP_NUM_THREADS=2
+
+zsrc=1100
+zsrc=0
+
+../fdelmodc \
+ file_cp=tilt_cp.su ischeme=1 iorder=4 \
+ file_den=tilt_ro.su \
+ file_src=wave.su \
+ file_rcv=shot_fd.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=0 \
+ rec_type_vz=1 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.01 \
+ dxrcv=10.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=$zsrc \
+ file_snap=snapF_$zsrc \
+ tsnap1=0.1 tsnap2=2.0 dtsnap=0.05 dxsnap=5 dzsnap=5 \
+ ntaper=101 \
+ snapwithbnd=1 \
+ left=2 right=2 top=2 bottom=2
+
+#suxmovie < snapF_${zsrc}_svz.su loop=1 clip=1e-13
+sugain scale=-1 < snapF_0_svz.su | sustrip > snapF_0_dx${dx}_svz.bin
+
+exit;
+../fdelmodc \
+ file_cp=hom_cp.su ischeme=1 iorder=4 \
+ file_den=hom_ro.su \
+ file_src=wave.su \
+ file_rcv=shoth_fd.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=0 \
+ rec_type_vz=1 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.01 \
+ dxrcv=10.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=$zsrc \
+ file_snap=snapH_$zsrc \
+ tsnap1=0.1 tsnap2=2.0 dtsnap=0.05 dxsnap=5 dzsnap=5 \
+ ntaper=101 \
+ snapwithbnd=1 \
+ left=2 right=2 top=2 bottom=2
+
+sudiff snapF_${zsrc}_svz.su snapH_${zsrc}_svz.su > snap_svz.su
+
+sugain scale=-1 < snap_svz.su | sustrip > snap_svz.bin
+
+suxmovie < snap_svz.su loop=1 clip=1e-13
diff --git a/fdelmodc/demo/fdelmodc_taper.scr b/fdelmodc/demo/fdelmodc_taper.scr
index e6dfeeaf79cebdd597bca0e25e5fcbe9222299ed..cdb8aed52b6996561755ac5a171097d9becb81d4 100755
--- a/fdelmodc/demo/fdelmodc_taper.scr
+++ b/fdelmodc/demo/fdelmodc_taper.scr
@@ -11,6 +11,7 @@
dt=0.0002
dx=1
cp=1500
+clip=1
makewave file_out=wavelet.su dt=$dt nt=1024 fp=85 shift=1 w=g2 verbose=1
@@ -22,7 +23,36 @@ export filecp=model_cp.su
export filecs=model_cs.su
export filero=model_ro.su
-for ntaper in 0 200;
+ntaper=0
+../fdelmodc \
+ file_cp=$filecp file_cs=$filecs file_den=$filero \
+ ischeme=1 \
+ file_src=wavelet.su verbose=4 \
+ file_rcv=rec.su \
+ file_snap=snap_nodisp.su \
+ xrcv1=0 xrcv2=1000 dxrcv=5 \
+ zrcv1=300 zrcv2=300 \
+ rec_type_vx=1 rec_type_pp=1 rec_type_ss=1 rec_int_vx=1 \
+ dtrcv=0.004 \
+ xsrc=500 zsrc=500 nshot=1 \
+ src_type=1 tmod=1.0 \
+ ntaper=$ntaper \
+ left=3 right=3 bottom=3 top=3 \
+ tsnap1=0.2 tsnap2=1.0 dtsnap=0.1 \
+ sna_type_ss=1 sna_type_pp=1 fmax=25
+
+suwind key=fldr min=3 max=3 < snap_nodisp_sp.su | \
+ supsimage \
+ wbox=4 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
+ d2=$dx f2=0 clip=$clip \
+ label1="depth [m]" label2="lateral position [m]" > snap_tap${ntaper}.eps
+
+supsimage < rec_rp.su \
+ wbox=3 hbox=4 titlesize=-1 labelsize=10 clip=$clip verbose=1 \
+ label1="time [s]" label2="lateral position [m]" > rec_tap${ntaper}_rp.eps
+
+
+for ntaper in 50 100;
do
../fdelmodc \
file_cp=$filecp file_cs=$filecs file_den=$filero \
@@ -44,14 +74,46 @@ do
suwind key=fldr min=3 max=3 < snap_nodisp_sp.su | \
supsimage \
wbox=4 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
- d2=$dx f2=0 clip=0.003 \
+ d2=$dx f2=0 clip=$clip \
label1="depth [m]" label2="lateral position [m]" > snap_tap${ntaper}.eps
supsimage < rec_rp.su \
- wbox=3 hbox=4 titlesize=-1 labelsize=10 clip=0.003 verbose=1 \
+ wbox=3 hbox=4 titlesize=-1 labelsize=10 clip=$clip verbose=1 \
label1="time [s]" label2="lateral position [m]" > rec_tap${ntaper}_rp.eps
done
+
+for npml in 5 10 20;
+do
+../fdelmodc \
+ file_cp=$filecp file_cs=$filecs file_den=$filero \
+ ischeme=1 \
+ file_src=wavelet.su verbose=4 \
+ file_rcv=rec.su \
+ file_snap=snap_nodisp.su \
+ xrcv1=0 xrcv2=1000 dxrcv=5 \
+ zrcv1=300 zrcv2=300 \
+ rec_type_vx=1 rec_type_pp=1 rec_type_ss=1 rec_int_vx=1 \
+ dtrcv=0.004 \
+ xsrc=500 zsrc=500 nshot=1 \
+ src_type=1 tmod=1.0 \
+ npml=$npml \
+ left=2 right=2 bottom=2 top=2 \
+ tsnap1=0.2 tsnap2=1.0 dtsnap=0.1 \
+ sna_type_ss=1 sna_type_pp=1 fmax=25
+
+suwind key=fldr min=3 max=3 < snap_nodisp_sp.su | \
+ supsimage \
+ wbox=4 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
+ d2=$dx f2=0 clip=$clip \
+ label1="depth [m]" label2="lateral position [m]" > snap_pml${npml}.eps
+
+supsimage < rec_rp.su \
+ wbox=3 hbox=4 titlesize=-1 labelsize=10 clip=$clip verbose=1 \
+ label1="time [s]" label2="lateral position [m]" > rec_pml${npml}_rp.eps
+
+done
+
exit;
diff --git a/fdemmodc/em4.c b/fdemmodc/em4.c
index 760399c1ea12f990e7959a404c40c8065efafaf5..fb6e6475bc17aa5fb373ff6fef485ed666113015 100644
--- a/fdemmodc/em4.c
+++ b/fdemmodc/em4.c
@@ -5,16 +5,15 @@
#include<assert.h>
#include"fdelmodc.h"
-int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float *vx, float *vz, float *tzz,
-float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, int verbose);
+int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, int verbose);
int storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
int reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
-int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int verbose);
+int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
-int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int verbose);
+int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *hz, float *hx, float *Ey, float *eprs, float *ksigma, float *mu, int verbose)
{
@@ -93,6 +92,7 @@ int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, in
c1*(Ey[ix*n1+iz] - Ey[(ix-1)*n1+iz]) +
c2*(Ey[(ix+1)*n1+iz] - Ey[(ix-2)*n1+iz]));
//if (hz[ix*n1+iz] > 0.1*FLT_MAX) fprintf(stderr,"%d: hz[%d %d] = %e\n", itime, ix, iz, hz[ix*n1+iz]);
+// if (hz[ix*n1+iz] != 0.0) fprintf(stderr,"%d: hz[%d %d] = %e\n", itime, ix, iz, hz[ix*n1+iz]);
}
}
@@ -105,16 +105,19 @@ int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, in
c1*(Ey[ix*n1+iz] - Ey[ix*n1+iz-1]) +
c2*(Ey[ix*n1+iz+1] - Ey[ix*n1+iz-2]));
//if (hx[ix*n1+iz] > 0.1*FLT_MAX) fprintf(stderr,"%d: hx[%d %d] = %e\n", itime, ix, iz, hx[ix*n1+iz]);
+// if (fabs(hx[ix*n1+iz]) != 0.0) fprintf(stderr,"%d: hx[%d %d] = %e\n", itime, ix, iz, hx[ix*n1+iz]);
}
}
/* Add force source */
if (src.type > 5) {
applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, hz, hx, Ey, NULL, NULL, mu, mu, eprs, src_nwav, verbose);
+// fprintf(stderr,"%d: hx[%d %d] = %e hz=%e \n", itime, ixsrc, izsrc, hx[ixsrc*n1+izsrc], hz[ixsrc*n1+izsrc]);
+//fprintf(stderr,"src.n = %d\n", src.n);
}
/* boundary condition clears velocities on boundaries */
- boundariesP(mod, bnd, hz, hx, Ey, NULL, NULL, mu, mu, eprs, NULL, NULL, verbose);
+ boundariesP(mod, bnd, hz, hx, Ey, NULL, NULL, mu, mu, eprs, NULL, NULL, itime, verbose);
/* calculate p/tzz for all grid points except on the virtual boundary */
for (ix=mod.ioPx; ix<mod.iePx; ix++) {
@@ -137,6 +140,7 @@ int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, in
for (iz=mod.ioPz; iz<mod.iePz; iz++) {
Ey[ix*n1+iz] -= eprs[ix*n1+iz]*(dzhx[iz]+dxhz[iz]) + ksigma[ix*n1+iz]*Ey[ix*n1+iz];
//if (Ey[ix*n1+iz] > 0.1*FLT_MAX) fprintf(stderr,"%d: Ey[%d %d] = %e\n", itime, ix, iz, Ey[ix*n1+iz]);
+// if (Ey[ix*n1+iz] != 0.0) fprintf(stderr,"%d: Ey[%d %d] = %e\n", itime, ix, iz, Ey[ix*n1+iz]);
}
}
diff --git a/fdemmodc/fdemmodc.c b/fdemmodc/fdemmodc.c
index fbd2ffb3f9a29ec3e8e275112161a6bbb525e9da..85d876c2f8dbd1aa164f4df4186b3c465dbdcff3 100644
--- a/fdemmodc/fdemmodc.c
+++ b/fdemmodc/fdemmodc.c
@@ -16,7 +16,7 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
int readEmModel(modPar mod, bndPar bnd, float *eprs, float *ksigma, float *mu);
-int defineSource(wavPar wav, srcPar src, float **src_nwav, int reverse, 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);
@@ -273,7 +273,7 @@ int main(int argc, char **argv)
}
}
- defineSource(wav, src, src_nwav, mod.grid_dir, verbose);
+ defineSource(wav, src, mod, src_nwav, mod.grid_dir, verbose);
/* allocate arrays for wavefield and receiver arrays */
diff --git a/fdemmodc/getEmParameters.c b/fdemmodc/getEmParameters.c
index 6e49d7eae143b782305953654b149198d8c3394b..95bb85f65b0034a34e7733e66568a95244c097c7 100644
--- a/fdemmodc/getEmParameters.c
+++ b/fdemmodc/getEmParameters.c
@@ -20,6 +20,8 @@
* The Netherlands
**/
+int optncr(int n);
+
float gaussGen();
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);
@@ -36,7 +38,7 @@ int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2,
int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *src, shotPar *shot, bndPar *bnd, int verbose)
{
int isnapmax1, isnapmax2, isnapmax, sna_nrsna;
- int n1, n2, nx, nz, nsrc, ix, axis, ioPz, is0;
+ int n1, n2, nx, nz, nsrc, ix, axis, ioPz, is0, optn;
int idzshot, idxshot;
int src_ix0, src_iz0, src_ix1, src_iz1;
int disable_check;
@@ -49,7 +51,7 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
float tsnap1, tsnap2, dtsnap, dxsnap, dzsnap, dtrcv;
float xsnap1, xsnap2, zsnap1, zsnap2, xmax, zmax;
float xsrc1, xsrc2, zsrc1, zsrc2, tsrc1, tsrc2, tlength, tactive;
- float src_angle, src_velo, p, grad2rad, rdelay;
+ float src_angle, src_velo, p, grad2rad, rdelay, scaledt;
float *xsrca, *zsrca, rrcv;
float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
size_t nsamp;
@@ -102,13 +104,28 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
/* define wavelet(s), modeling time and wavelet maximum frequency */
if (wav->file_src!=NULL) {
- getWaveletInfo(wav->file_src, &wav->nt, &wav->nx, &wav->dt, &d2, &f1, &f2, &fmax, &ntraces, verbose);
- if (wav->dt <= 0.0) {
+ getWaveletInfo(wav->file_src, &wav->ns, &wav->nx, &wav->ds, &d2, &f1, &f2, &fmax, &ntraces, verbose);
+ if (wav->ds <= 0.0) {
vwarn("dt in wavelet (file_src) equal to 0.0 or negative.");
vwarn("Use parameter dt= to overule dt from file_src.");
}
+ wav->nt = wav->ns;
+ wav->dt = wav->ds;
if(!getparfloat("tmod",&mod->tmod)) mod->tmod = (wav->nt-1)*wav->dt;
if(!getparfloat("dt",&mod->dt)) mod->dt=wav->dt;
+ if (NINT(wav->ds*1000000) != NINT(mod->dt*1000000)) {
+ if (wav->dt > mod->dt) {
+ scaledt = wav->dt/mod->dt;
+ scaledt = floorf(wav->dt/mod->dt);
+ optn = optncr(wav->ns);
+ wav->nt = floorf(scaledt*optn);
+ vmess("file_src dt-scalefactor=%f : wav.dt=%e ==interpolated==> mod.dt=%e", scaledt, wav->dt, mod->dt);
+ wav->dt = mod->dt;
+ }
+ else {
+ wav->dt = mod->dt; /* in case if wav.dt is smaller than 1e-7 and can not be read by SU-getpar */
+ }
+ }
if(!getparfloat("fmax",&wav->fmax)) wav->fmax=fmax;
}
else {
@@ -714,7 +731,7 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
vmess("*******************************************");
vmess("************* wavelet info ****************");
vmess("*******************************************");
- vmess("wav_nt = %6d wav_nx = %d", wav->nt, wav->nx);
+ vmess("wav_nt = %6d wav_nx = %d", wav->ns, wav->nx);
vmess("src_type = %6d src_orient = %d", src->type, src->orient);
vmess("fmax = %10.3e", fmax);
fprintf(stderr," %s: Source type : ",xargv[0]);
diff --git a/raytime/Makefile b/raytime/Makefile
index 11948fc94c7a5e51bcb170adc4bca2f9d2dbe801..bc429087cf0e991413e7a471300071241132b0e5 100644
--- a/raytime/Makefile
+++ b/raytime/Makefile
@@ -26,6 +26,7 @@ PRG = raytime
SRCC = $(PRG).c \
JespersRayTracer.c \
vidale.c \
+ Grid2Time1.c \
getParameters.c \
getWaveletInfo.c \
writeSrcRecPos.c \
diff --git a/utils/Makefile b/utils/Makefile
index 8f6b1a5670dd97e30248a9c8780a61751851fd71..3b59b4f71f844acc56c9b9851deba24861d82e5b 100644
--- a/utils/Makefile
+++ b/utils/Makefile
@@ -86,6 +86,7 @@ SRCG = green.c \
SRCB = basop.c \
getFileInfo.c \
+ kxwfilter.c \
readData.c \
writeData.c \
wallclock_time.c \
diff --git a/utils/basop.c b/utils/basop.c
index 98a10f82e53cf667e53e0982ac7d76be2d23426d..4456acbaaf94921101fd404e8eef7307690a6ee5 100644
--- a/utils/basop.c
+++ b/utils/basop.c
@@ -39,6 +39,7 @@ void timeRotate(float *data, int nsam, int nrec, int nrot);
void rma(float *data, int nsam, int nrec);
void rmat(float *data, int nsam, int nrec);
void decompAcoustic(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, float rho, int opt);
+void kxwfilter(float *data, int nt, int nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc);
int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
@@ -70,6 +71,7 @@ char *sdoc[] = {
" shift=0 .................. time shift in seconds",
" rot=0 .................... phase rotation in parts of PI",
" c=1500 ................... velocity of the medium (for kz and ghost)",
+" alpha=65 ................. angle for kfilt option: k > (w/c)*cos(alpha) is filtered out",
" rho=1000 ................. density in kg/m^3 of the medium (for decomposition operator)",
" dz=9 ..................... depth of receivers for deghosting (m)",
" eps=0.01 ................. stabilization for deghosting",
@@ -105,6 +107,7 @@ char *sdoc[] = {
" - 23 - som = multiply with sqrt(w) in frequency domain",
" - 24 - deca1 = acoustic decompostion multiply with sqrt(2 kz/(w rho)) in kx-w domain",
" - 25 - deca2 = acoustic decompostion multiply with sqrt((2 w rho)/(kz)) in kx-w domain",
+" - 26 - kfilt = dipfilter in the kx-w domain given alpha and cp",
" ",
" author : Jan Thorbecke : 12-12-1994 (janth@xs4all.nl)",
" Alexander Koek (E.A.Koek@CTG.TuDelft.NL)",
@@ -124,7 +127,7 @@ int main(int argc, char **argv)
float dt, dx, dy, c, rho, d1, d2, f1, f2;
float scl, xmin, xmax, trot;
double t0, t1, t2;
- float fmin, fmax, *data, shift, rot, dz, eps, *tmpdata;
+ float fmin, fmax, *data, shift, rot, alpha, perc, dz, eps, *tmpdata;
char *file_in, *file_out;
char choice[10], *choicepar;
segy *hdrs;
@@ -148,6 +151,7 @@ int main(int argc, char **argv)
if(!getparfloat("fmin", &fmin)) fmin = 0.0;
if(!getparfloat("shift", &shift)) shift = 0.0;
if(!getparfloat("c", &c)) c = 1500.0;
+ if(!getparfloat("alpha", &alpha)) alpha = 65.0;
if(!getparfloat("rho", &rho)) rho = 1000.0;
if(!getparfloat("dz", &dz)) dz = 9;
if(!getparfloat("eps", &eps)) eps = 0.01;
@@ -379,7 +383,11 @@ int main(int argc, char **argv)
opt = 2;
decompAcoustic(data, nsam, nrec, dt, dx, fmin, fmax, c, rho, opt);
}
-
+ else if (!strcmp(choice, "kfilt") || !strcmp(choice, "26")) {
+ if (verbose) vmess("kx-w filtering for angles > alpha");
+ perc=0.15;
+ kxwfilter(data, nsam, nrec, dt, dx, fmin, fmax, alpha, c, perc);
+ }
t2 = wallclock_time();
diff --git a/utils/demo/green.scr b/utils/demo/green.scr
index 36d4c6d2d98599858f26c4efda6c1c1e986f6397..d573a519748a84172d608704f3ef5633a8c57f6d 100755
--- a/utils/demo/green.scr
+++ b/utils/demo/green.scr
@@ -1,9 +1,9 @@
#!/bin/bash
-makewave nt=512 w=g2 file_out=wave.su w=g2 dt=0.004
+makewave nt=512 file_out=wave.su w=g2 dt=0.004
green c=2000 nt=512 dt=0.004 zsrc1=500 xsrc1=0 file_src=wave.su file_out=green.su
-fconv file_in1=green.su file_in2=wave.su mode=dec eps=0.0 reps=0.01 file_out=dec.su mode=dec verbose=1
+../fconv file_in2=green.su file_in1=wave.su eps=0.01 reps=0.01 file_out=dec.su mode=dec verbose=2
diff --git a/utils/fconv.c b/utils/fconv.c
index ed409fbbb60ea0db6a7497da7a724493575d1d49..6afc67774ad68e9bf7c6a7d9fee23be0df667765 100644
--- a/utils/fconv.c
+++ b/utils/fconv.c
@@ -647,9 +647,11 @@ void deconv(float *data1, float *data2, float *decon, int nrec, int nsam,
qr = (float *) &cdec[0].r;
p1i = p1r + 1;
p2i = p2r + 1;
+ qi = qr + 1;
leps = reps*maxden+eps;
// fprintf(stderr,"eps=%e reps=%e max=%e => leps=%e\n", eps, reps, maxden, leps);
for (j = 0; j < n; j++) {
+
if (fabs(*p2r)>=fabs(*p2i)) {
*qr = (*p2r**p1r+*p2i**p1i)/(den[j]+leps);
*qi = (*p2r**p1i-*p2i**p1r)/(den[j]+leps);
diff --git a/utils/kxwfilter.c b/utils/kxwfilter.c
new file mode 100644
index 0000000000000000000000000000000000000000..30396ec31f8ca29348ac7d42f0af9d553a8939d3
--- /dev/null
+++ b/utils/kxwfilter.c
@@ -0,0 +1,137 @@
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <genfft.h>
+
+#define ISODD(n) ((n) & 01)
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+void kxwfilt(complex *data, float k, float dx, int nkx, float a1, float perc);
+
+void kxwfilter(float *data, int nt, int nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc)
+{
+ int ntfft, nfreq, nkx, ix, it, iomin, iomax, xorig, iom, ikx;
+ float df, dkx, deltom, omin, omax, kp, om;
+ float *pdata, *filter;
+ complex *cdata;
+
+ ntfft = optncr(nt);
+ nfreq = ntfft/2+1;
+ nkx = optncc(2*nx);
+
+ df = 1.0/((float)ntfft*dt);
+ dkx = 2.0*M_PI/(nkx*dx);
+ deltom = 2.*M_PI*df;
+ omin = 2.*M_PI*fmin;
+ omax = 2.*M_PI*fmax;
+
+ iomin = (int)MIN((omin/deltom), (nfreq-1));
+ iomin = MAX(iomin, 1);
+ iomax = MIN((int)(omax/deltom), (nfreq-1));
+
+ pdata = (float *)calloc(ntfft*nkx,sizeof(float));
+ cdata = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ filter = (float *)malloc(nkx*sizeof(float));
+
+ /* copy input data into extended arrays with padded zeroes */
+ for (ix=0; ix<nx; ix++) {
+ memcpy(&pdata[ix*ntfft],&data[ix*nt],nt*sizeof(float));
+ }
+
+ /* transform from t-x to kx-w */
+ xorig = nkx/2;
+ xt2wkx(pdata, cdata, ntfft, nkx, ntfft, nkx, xorig);
+
+
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+ kp = om/cp;
+
+ kxwfilt(&cdata[iom*nkx], kp, dx, nkx, angle, perc);
+
+ for (ikx = 0; ikx < nkx; ikx++) {
+ cdata[ikx].r *= filter[ikx];
+ cdata[ikx].i *= filter[ikx];
+ }
+ }
+
+ /* transform back to t-x */
+ wkx2xt(cdata, pdata, ntfft, nkx, nkx, ntfft, xorig);
+
+ /* reduce array to nt samples nx traces */
+ for (ix=0; ix<nx; ix++) {
+ memcpy(&data[ix*nt],&pdata[ix*ntfft],nt*sizeof(float));
+ }
+
+ free(pdata);
+ free(cdata);
+
+ return;
+}
+
+
+void kxwfilt(complex *data, float k, float dx, int nkx, float a1, float perc)
+{
+ int ikx, ik1, ik2, ntap;
+ float kxnyq, dkx, kxfmax, kfilt;
+ float kpos, band, li, *filter;
+
+ kxnyq = M_PI/dx;
+ dkx = 2.0*M_PI/(nkx*dx);
+ if (a1 > 90.0) kpos = kxnyq;
+ else kpos = k*sin(M_PI*a1/180.0);
+
+ filter = (float *)malloc(nkx*sizeof(float));
+
+ band = fabs(2*kpos);
+ ntap = (int)fabs((int)(perc*band/dkx));
+ kfilt = fabs(dkx*ntap);
+
+ if (perc > 0) {
+ if (kpos+kfilt < kxnyq) {
+ kxfmax = kpos+kfilt;
+ }
+ else {
+ kxfmax = kxnyq;
+ ntap = (int)(0.15*nkx/2);
+ }
+ }
+ else {
+ kxfmax = MIN(kpos, kxnyq);
+ }
+
+ ik1 = (int)(kxfmax/dkx);
+ ik2 = ik1 - ntap;
+
+ if (perc < -0.5 || perc > 1.0) {
+ if (kpos > 0.85*kxnyq) {
+ kpos = 0.85*kxnyq;
+ }
+ ik1 = nkx/2-1;
+ ik2 = (int)(kpos/dkx);
+ }
+
+ li = 1.0/(ik1-ik2);
+ for (ikx = 0; ikx < ik2; ikx++)
+ filter[ikx] = 1.0;
+ for (ikx = ik2; ikx < ik1; ikx++)
+ filter[ikx] = 0.5*(cos(M_PI*(ikx-ik2)*li)+1);
+ for (ikx = ik1; ikx <= nkx/2; ikx++)
+ filter[ikx] = 0.0;
+
+ for (ikx = 0; ikx <= (nkx/2); ikx++) {
+ data[ikx].r *= filter[ikx];
+ data[ikx].i *= filter[ikx];
+ }
+ for (ikx = (nkx/2+1); ikx < nkx; ikx++) {
+ data[ikx].r *= filter[nkx-ikx];
+ data[ikx].i *= filter[nkx-ikx];
+ }
+
+ free(filter);
+ return;
+}
+
+