diff --git a/README.md b/README.md
index e31ebf72ba99e49263e9ff91331cf8454f532f6d..28e547d72dfdaf18fde68e4d9074d73b5cd58c96 100644
--- a/README.md
+++ b/README.md
@@ -180,6 +180,17 @@ Other make commands which can be useful:
make clean : removes all object files, but leaves libraries and executables
make realclean: removes also object files, libraries and executables.
+COMPILATION PROBLEMS
+--------------------
+If you encounter missing trigometric / mathematical functions during compilation, for example;
+
+defineSource.c:(.text+0x2356): undefined reference to sin getParameters.o: In function getParameters:
+
+add '-lm -lc' around line 109 in Make_include:
+
+LIBS= -L$L -lgenfft $(BLAS) -lm -lc
+
+
UPDATES AND LATEST VERSION
--------------------------
diff --git a/fdelmodc/decomposition.c b/fdelmodc/decomposition.c
index 0b4a3e7cd040eaa7cac60748ecd323a6e438dc1f..a47b41e1cc7292522fbc553d9fdc2a461a3d4260 100644
--- a/fdelmodc/decomposition.c
+++ b/fdelmodc/decomposition.c
@@ -37,6 +37,7 @@ complex cwp_csqrt(complex z);
void decudP(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu);
void decudVz(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu);
+void decudF(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu, complex *ipu);
int writesufile(char *filename, float *data, size_t n1, size_t n2, float f1, float f2, float d1, float d2);
@@ -50,8 +51,8 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
int iom, iomin, iomax, ikx, nfreq, a, av;
float omin, omax, deltom, om, df, dkx;
float alpha, eps, *angle, avrp, avrvz, maxrp, maxrvz;
- float fangle, pangle, vangle, kangle;
- complex *pu;
+ float fangle, pangle, vangle, kangle, scl;
+ complex *pu, *ipu;
complex ax, az;
df = 1.0/((float)nt*dt);
@@ -68,6 +69,7 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
iomax = MIN((int)(omax/deltom), (nfreq-1));
pu = (complex *)malloc(nkx*sizeof(complex));
+ ipu = (complex *)malloc(nkx*sizeof(complex));
angle = (float *)calloc(2*90,sizeof(float));
/* estimate maximum propagation angle in wavefields P and Vz */
@@ -129,7 +131,7 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
if(!getparfloat("kangle",&kangle)) kangle=fangle;
if (verbose>=2) vmess("Up-down going: maximum angle in decomposition= %f", kangle);
- if (vznorm) { /* Vz normalised decompostion */
+ if (vznorm==1) { /* Vz normalised decompostion */
for (iom = iomin; iom <= iomax; iom++) {
om = iom*deltom;
@@ -152,7 +154,7 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
}
}
- else { /* P normalised decompostion */
+ else if (vznorm==0) { /* P normalised decompostion */
for (iom = iomin; iom <= iomax; iom++) {
om = iom*deltom;
@@ -171,8 +173,31 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
}
}
}
+ else if (vznorm==-1) { /* Flux normalised decompostion */
+ scl = 1.0/sqrt(2.0);
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+
+ decudF(om, rho, cp, dx, nkx, kangle, alpha, eps, pu, ipu);
+
+ for (ikx = 0; ikx < nkx; ikx++) {
+ ax.r = scl*(rp[iom*nkx+ikx].r*pu[ikx].r - rp[iom*nkx+ikx].i*pu[ikx].i);
+ ax.i = scl*(rp[iom*nkx+ikx].i*pu[ikx].r + rp[iom*nkx+ikx].r*pu[ikx].i);
+ az.r = scl*(rvz[iom*nkx+ikx].r*ipu[ikx].r - rvz[iom*nkx+ikx].i*ipu[ikx].i);
+ az.i = scl*(rvz[iom*nkx+ikx].i*ipu[ikx].r + rvz[iom*nkx+ikx].r*ipu[ikx].i);
+
+ down[iom*nkx+ikx].r = ax.r + az.r;
+ down[iom*nkx+ikx].i = ax.i + az.i;
+ up[iom*nkx+ikx].r = ax.r - az.r;
+ up[iom*nkx+ikx].i = ax.i - az.i;
+ //down[iom*nkx+ikx] = ax;
+ //up[iom*nkx+ikx] = az;
+ }
+ }
+ }
free(pu);
+ free(ipu);
free(angle);
return;
@@ -328,6 +353,116 @@ void decudVz(float om, float rho, float cp, float dx, int nkx, float kangle, flo
return;
}
+/* Flux normalised decompostion */
+
+void decudF(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu, complex *ipu)
+{
+ int ikx, ikxmax1, ikxmax2, filterpoints, filterppos;
+ float kp, kp2, scl, kzpr;
+ float kx, kx2, kzp2, dkx, stab;
+ float kxfmax, kxnyq, kpos, kfilt, perc, band, *filter;
+ complex kzp, ckp, ckp2, kzpq, kzpi;
+
+ kp = om/cp;
+ kp2 = kp*kp;
+ dkx = 2.0*M_PI/(nkx*dx);
+ stab = eps*eps*kp*kp;
+ scl = (om*rho);
+
+ /* make kw filter at maximum angle alfa */
+ perc = 0.15; /* percentage of band to use for smooth filter */
+ filter = (float *)malloc(nkx*sizeof(float));
+ kpos = kp*sin(M_PI*kangle/180.0);
+ kxnyq = M_PI/dx;
+ if (kpos > kxnyq) kpos = kxnyq;
+ band = kpos;
+ filterpoints = (int)abs((int)(perc*band/dkx));
+ kfilt = fabsf(dkx*filterpoints);
+ if (kpos+kfilt < kxnyq) {
+ kxfmax = kpos+kfilt;
+ filterppos = filterpoints;
+ }
+ else {
+ kxfmax = kxnyq;
+ filterppos = (int)(0.15*nkx/2);
+ }
+ ikxmax1 = (int) (kxfmax/dkx);
+ ikxmax2 = ikxmax1 - filterppos;
+ // fprintf(stderr,"ikxmax1=%d ikxmax2=%d nkp=%d nkx=%d\n", ikxmax1, ikxmax2, (int)(kp/dkx), nkx);
+
+ for (ikx = 0; ikx < ikxmax2; ikx++)
+ filter[ikx]=1.0;
+ for (ikx = ikxmax2; ikx < ikxmax1; ikx++)
+ filter[ikx] =(cos(M_PI*(ikx-ikxmax2)/(ikxmax1-ikxmax2))+1)/2.0;
+ for (ikx = ikxmax1; ikx <= nkx/2; ikx++)
+ filter[ikx] = 0.0;
+ /* end of kxfilter */
+
+ for (ikx = 0; ikx <= (nkx/2); ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kzp2 = (kp2 - kx2);
+ if (kzp2 <= 0.0) continue;
+
+ kzpr = sqrt(kzp2);
+
+ if (kzp2 != 0) {
+ pu[ikx].r = filter[ikx]*sqrt(kzpr/scl);
+ pu[ikx].i = filter[ikx]*sqrt(kzpr/scl);
+ ipu[ikx].r = filter[ikx]*sqrt(scl/kzpr);
+ ipu[ikx].i = filter[ikx]*sqrt(scl/kzpr);
+ }
+ else {
+ ipu[ikx].r = 0.0;
+ ipu[ikx].i = 0.0;
+ }
+ }
+
+/* operators are symmetric in kx-w domain */
+ for (ikx = (nkx/2+1); ikx < nkx; ikx++) {
+ pu[ikx] = pu[nkx-ikx];
+ ipu[ikx] = ipu[nkx-ikx];
+ }
+
+ free(filter);
+
+ return;
+}
+
+void findMatch(complex *A, complex *B, int nx, float om, float trange, float arange, float *topt, float *aopt)
+{
+ int ix, it, ia;
+ float match, tb, ab;
+ float tshift, ampl, corr, *shifted;
+
+
+ shifted = (float *)malloc(nx*sizeof(float));
+ corr = 0.0;
+ for (ix=0; ix<nx; ix++) {
+ corr += A[ix].r*B[ix].r - A[ix].i*B[ix].i;
+ corr += A[ix].i*B[ix].r + A[ix].r*B[ix].i;
+ }
+
+ for (it=0; it<100; it++) {
+ tshift = -0.5*trange+it*(trange/100.0);
+ tom = om*shift;
+ cdatascl[ix*nfreq+iom].r = cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom);
+ cdatascl[ix*nfreq+iom].i = cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom);
+
+ for (ia=0; ia<100; ia++) {
+ ampl = -0.5*arange+it*(arange/100.0);
+ match = 0.0;
+ for (ix=0; ix<nx; ix++) {
+ match += A[ix].r*B[ix].r + A[ix].i*B[ix].i;
+ match += A[ix].r*B[ix].r + A[ix].i*B[ix].i;
+ }
+ if (match > corr) { tb = tshift; ab = ampl; }
+ }
+ }
+
+ free(shifted);
+ return;
+}
complex froot(float x)
{
@@ -345,7 +480,7 @@ complex froot(float x)
}
-/* compute 1/x */
+/* compute 1/sqrt(x) */
complex firoot(float x, float stab)
{
complex z;
diff --git a/fdelmodc/demo/modelOilGas.scr b/fdelmodc/demo/modelOilGas.scr
index d7aa0687d2d92cb86badf6a06d9d1a7f25281ea0..7bf177e779a39c09d01fb3345b029dc0e2c64833 100755
--- a/fdelmodc/demo/modelOilGas.scr
+++ b/fdelmodc/demo/modelOilGas.scr
@@ -29,11 +29,12 @@ which fdelmodc
file_den=syncl_ro.su \
file_src=wave.su \
file_rcv=shot_fd.su \
- src_type=7 \
+ src_type=10 \
src_orient=1 \
src_injectionrate=0 \
- rec_type_vz=1 \
+ rec_type_vz=2 \
rec_type_p=1 \
+ rec_type_ud=3 \
rec_int_vz=2 \
dtrcv=0.004 \
rec_delay=0.1 \
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 1b8f0eadd3e162ec45ec20837d557a93d68da7fb..76229d47e7f019f30f47035cf2ca6ed8e79c0a56 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -234,6 +234,7 @@ char *sdoc[] = {
" rec_type_ss=0 ..... S (curl) registration _rS",
" rec_type_ud=0 ..... 1:pressure normalized decomposition in up and downgoing waves _ru, _rd",
" ................... 2:particle velocity normalized decomposition in up and downgoing waves _ru, _rd",
+" ................... 3:flux normalized decomposition in up and downgoing waves _flup, _flip",
" kangle= ........... maximum wavenumber angle for decomposition",
" rec_int_vx=0 ..... interpolation of Vx receivers",
" - 0=Vx->Vx (no interpolation)",
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index 598267d5a2c6ccb1e90f86b9d2586ec6abfb7f13..61983b0d2b39f2ed4ea83a5006cac6e9bf77d2f4 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -935,6 +935,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
}
}
if (src->type==7) { /* set also src_injectionrate=1 */
+ vwarn("For src_type=7 injectionrate is always set to 1");
src->injectionrate=1;
}
@@ -1233,7 +1234,9 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (rec->type.vz) fprintf(stderr,"Vz ");
if (rec->type.vx) fprintf(stderr,"Vx ");
if (rec->type.p) fprintf(stderr,"p ");
- if (rec->type.ud) fprintf(stderr,"P+ P- ");
+ if (rec->type.ud==1) fprintf(stderr,"P+ P- Pressure normalized");
+ if (rec->type.ud==2) fprintf(stderr,"P+ P- Particle Velocity normalized");
+ if (rec->type.ud==3) fprintf(stderr,"P+ P- Flux normalized");
if (mod->ischeme>2) {
if (rec->type.txx) fprintf(stderr,"Txx ");
if (rec->type.tzz) fprintf(stderr,"Tzz ");
diff --git a/fdelmodc/writeRec.c b/fdelmodc/writeRec.c
index 982cc6bcdcaadd3eb0e515087ddef39f955d58d4..c8ab3b7a2c5fa892b262448b059e4821f5264ca9 100644
--- a/fdelmodc/writeRec.c
+++ b/fdelmodc/writeRec.c
@@ -122,8 +122,9 @@ int writeRec(recPar rec, modPar mod, bndPar bnd, wavPar wav, int ixsrc, int izsr
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
cp = rec.cp;
rho = rec.rho;
- if (rec.type.ud==2) vznorm=1;
- else vznorm=0;
+ if (rec.type.ud==1) vznorm=0;
+ else if (rec.type.ud==2) vznorm=1;
+ else if (rec.type.ud==3) vznorm=-1;
if (verbose) vmess("Decomposition array at z=%.2f with cp=%.2f rho=%.2f", rec.zr[0]+mod.z0, cp, rho);
rec_up = (float *)calloc(ntfft*nkx,sizeof(float));
rec_down= (float *)calloc(ntfft*nkx,sizeof(float));
@@ -167,6 +168,7 @@ int writeRec(recPar rec, modPar mod, bndPar bnd, wavPar wav, int ixsrc, int izsr
free(crec_dw);
}
if (rec.type.ud && (mod.ischeme==3 || mod.ischeme==4) ) {
+// Not yet implemented
}
for (irec=0; irec<rec.n; irec++) {
diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr
index 94b37e1ca0081efa4b5b1581d24b6e444478d768..55ff6255c733218bb040ac5add82e5153c93f7af 100755
--- a/marchenko/demo/oneD/primaries.scr
+++ b/marchenko/demo/oneD/primaries.scr
@@ -19,7 +19,7 @@ select=451
makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
-../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
+../../marchenko_primaries1 file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
niter=31 shift=20 smooth=10 niterec=2 niterskip=50 file_rr=pred_rr.su T=0 file_update=update.su
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 0738c465101ebb72d5582d08989358ae297f80f1..3fd1f9bca86f7c5267bead3e0da3d75f8997b3cf 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -87,6 +87,8 @@ char *sdoc[] = {
" hw=8 ..................... window in time samples to look for maximum in next trace",
" smooth=5 ................. number of points to smooth mute with cosine window",
" plane_wave=0 ............. enable plane-wave illumination function",
+" src_angle=0 .............. angle of plane source array",
+" src_velo=1500 ............ velocity to use in src_angle definition",
" REFLECTION RESPONSE CORRECTION ",
" tsq=0.0 .................. scale factor n for t^n for true amplitude recovery",
" Q=0.0 .......,............ Q correction factor",
@@ -132,6 +134,7 @@ int main (int argc, char **argv)
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
float xmin, xmax, scale, tsq, Q, f0;
float *ixmask;
+ float grad2rad, p, src_angle, src_velo, tshift;
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
@@ -176,6 +179,8 @@ int main (int argc, char **argv)
if(!getparint("shift", &shift)) shift=12;
if (!getparint("plane_wave", &plane_wave)) plane_wave = 0;
+ if (!getparfloat("src_angle",&src_angle)) src_angle=0.;
+ if (!getparfloat("src_velo",&src_velo)) src_velo=1500.;
if (reci && ntap) vwarn("tapering influences the reciprocal result");
@@ -251,9 +256,18 @@ int main (int argc, char **argv)
/* compute time shift for tilted plane waves */
if (plane_wave==1) {
itmin = nt;
+ /* compute time shift for shifted plane waves */
+ grad2rad = 17.453292e-3;
+ p = sin(src_angle*grad2rad)/src_velo;
+ tshift = fabs((nxs-1)*dxs*p);
+
+ /* compute mute window for plane waves */
+ //for (i=0; i<nxs; i++) fprintf(stderr,"i=%d window=%d\n", i, muteW[i]);
for (i=0; i<nxs; i++) itmin = MIN (itmin, muteW[i]);
for (i=0; i<nxs; i++) tsynW[i] = muteW[i]-itmin;
if (Nfoc!=1) verr("For plane-wave focusing only one function can be computed at the same time");
+ //fprintf(stderr,"itmin=%d\n", itmin);
+ //for (i=0; i<nxs; i++) fprintf(stderr,"i=%d window=%f\n", i, tsynW[i]*dt);
}
else { /* just fill with zero's */
itmin=0;
@@ -568,19 +582,20 @@ int main (int argc, char **argv)
/* Apply mute with window for Gmin */
if ( plane_wave==1 ) {
applyMute_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 0, tsynW);
- /* for plane wave with angle shift itmin downward */
+ itmin = 0;
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
memcpy(&trace[0],&Gmin[l*nxs*nts+i*nts],nts*sizeof(float));
for (j = 0; j < itmin; j++) {
Gmin[l*nxs*nts+i*nts+j] = 0.0;
}
- for (j = 0; j < nts-itmin; j++) {
- Gmin[l*nxs*nts+i*nts+j+itmin] = trace[j];
+ for (j = 0; j < nts; j++) {
+ Gmin[l*nxs*nts+i*nts+j] = trace[j];
}
}
}
- //timeShift(Gmin, nts, npos, dt, itmin*dt, fmin, fmax);
+ /* for plane wave with angle shift itmin downward */
+ timeShift(Gmin, nts, npos, dt, tshift, fmin, fmax);
}
else {
applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index 4ddfb85b2bfe27ec59fcc913c759248db7c9aaed..1eddba67121cafede753f57e4f0a3b7670dca991 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -639,12 +639,33 @@ int main (int argc, char **argv)
}
free(SRC);
+/* For testing different window functions */
+
+/*
+ char filename[1024];
+ int A, W;
+ A=shift*3; W=10;
+ sprintf(filename,"windowA%dW%d.txt",A, W);
+ fp_up = fopen(filename, "w+");
+ ii=250;
+ for (i = 0; i < npos; i++) {
+ twplane[i] = dt*A*sin(2*M_PI*i*W/npos);
+ fprintf(fp_up,"time-shift= %d %d %f\n", i, NINT(twplane[i]/dt), twplane[i]);
+ }
+ fclose(fp_up);
+*/
+
/*================ start loop over number of time-samples ================*/
for (ii=istart; ii<iend; ii++) {
/*================ initialization ================*/
+/*
+ for (i = 0; i < npos; i++) {
+ twplane[i] = sqrt(dxs*dxs*(i-npos/2)*(i-npos/2)+ii*ii*dt*dt*2200*2200)/2200-ii*dt;
+ }
+*/
t5 = wallclock_time();
memset(M0, 0, Nfoc*nxs*ntfft*sizeof(float));
memset(v1plus, 0, Nfoc*nxs*ntfft*sizeof(float));
diff --git a/zfp/Makefile b/zfp/Makefile
index 1f73d1047f2d720d0fa03fb8136b22248467a9f4..4b83d4fb44f73d49599a083df0e305fec4fda7aa 100644
--- a/zfp/Makefile
+++ b/zfp/Makefile
@@ -44,7 +44,7 @@ ifneq ($(BUILD_EXAMPLES),0)
@cd examples; $(MAKE) clean all
endif
-install:
+install: all
-cp -rp include/* $I
-cp -rp lib/* $L