From 67f11bde7f5ce3d9597a21794d34d253d196d2f3 Mon Sep 17 00:00:00 2001
From: Jan Thorbecke <janth@xs4all.nl>
Date: Thu, 12 Apr 2018 14:45:33 +0200
Subject: [PATCH] added Velocity normalised decomposition
---
fdelmodc/decomposition.c | 149 +++++++++++++++++++++++++++-----
fdelmodc/demo/decomposition.scr | 31 ++++++-
fdelmodc/fdelmodc.c | 3 +-
fdelmodc/writeRec.c | 8 +-
4 files changed, 161 insertions(+), 30 deletions(-)
diff --git a/fdelmodc/decomposition.c b/fdelmodc/decomposition.c
index 3b48778..130d52f 100644
--- a/fdelmodc/decomposition.c
+++ b/fdelmodc/decomposition.c
@@ -31,10 +31,13 @@ typedef struct _dcomplexStruct { /* complex number */
void vmess(char *fmt, ...);
complex firoot(float x, float stab);
+complex froot(float x);
complex ciroot(complex x, float stab);
complex cwp_csqrt(complex z);
-void decud(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu);
+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);
+
int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2, float d1, float d2);
void kxwfilter(complex *data, float k, float dx, int nkx,
@@ -42,7 +45,7 @@ void kxwfilter(complex *data, float k, float dx, int nkx,
void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
int nkx, float dx, int nt, float dt, float fmin, float fmax,
- float cp, float rho, int verbose)
+ float cp, float rho, int vznorm, int verbose)
{
int iom, iomin, iomax, ikx, nfreq, a, av;
float omin, omax, deltom, om, df, dkx;
@@ -126,27 +129,48 @@ 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);
- for (iom = iomin; iom <= iomax; iom++) {
- om = iom*deltom;
+ if (vznorm) { /* Vz normalised decompostion */
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
+
+ decudVz(om, rho, cp, dx, nkx, kangle, alpha, eps, pu);
+ /*
+ kxwfilter(dpux, kp, dx, nkx, alfa1, alfa2, perc);
+ kxwfilter(dpuz, kp, dx, nkx, alfa1, alfa2, perc);
+ */
+ for (ikx = 0; ikx < nkx; ikx++) {
+ ax.r = 0.5*rvz[iom*nkx+ikx].r;
+ ax.i = 0.5*rvz[iom*nkx+ikx].i;
+ az.r = 0.5*(rp[iom*nkx+ikx].r*pu[ikx].r - rp[iom*nkx+ikx].i*pu[ikx].i);
+ az.i = 0.5*(rp[iom*nkx+ikx].i*pu[ikx].r + rp[iom*nkx+ikx].r*pu[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;
+ }
+ }
+
+ }
+ else { /* P normalised decompostion */
+ for (iom = iomin; iom <= iomax; iom++) {
+ om = iom*deltom;
- decud(om, rho, cp, dx, nkx, kangle, alpha, eps, pu);
-/*
- kxwfilter(dpux, kp, dx, nkx, alfa1, alfa2, perc);
- kxwfilter(dpuz, kp, dx, nkx, alfa1, alfa2, perc);
-*/
- for (ikx = 0; ikx < nkx; ikx++) {
- ax.r = 0.5*rp[iom*nkx+ikx].r;
- ax.i = 0.5*rp[iom*nkx+ikx].i;
- az.r = 0.5*(rvz[iom*nkx+ikx].r*pu[ikx].r - rvz[iom*nkx+ikx].i*pu[ikx].i);
- az.i = 0.5*(rvz[iom*nkx+ikx].i*pu[ikx].r + rvz[iom*nkx+ikx].r*pu[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;
- }
+ decudP(om, rho, cp, dx, nkx, kangle, alpha, eps, pu);
- }
+ for (ikx = 0; ikx < nkx; ikx++) {
+ ax.r = 0.5*rp[iom*nkx+ikx].r;
+ ax.i = 0.5*rp[iom*nkx+ikx].i;
+ az.r = 0.5*(rvz[iom*nkx+ikx].r*pu[ikx].r - rvz[iom*nkx+ikx].i*pu[ikx].i);
+ az.i = 0.5*(rvz[iom*nkx+ikx].i*pu[ikx].r + rvz[iom*nkx+ikx].r*pu[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;
+ }
+ }
+ }
free(pu);
free(angle);
@@ -154,7 +178,9 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
return;
}
-void decud(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu)
+/* Pressure normalised decompostion */
+
+void decudP(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu)
{
int ikx, ikxmax1, ikxmax2, filterpoints, filterppos;
float kp, kp2;
@@ -240,6 +266,85 @@ void decud(float om, float rho, float cp, float dx, int nkx, float kangle, float
return;
}
+/* Particle Velocity normalised decompostion */
+
+void decudVz(float om, float rho, float cp, float dx, int nkx, float kangle, float alpha, float eps, complex *pu)
+{
+ int ikx, ikxmax1, ikxmax2, filterpoints, filterppos;
+ float kp, kp2;
+ float kx, kx2, kzp2, dkx, stab;
+ float kxfmax, kxnyq, kpos, kfilt, perc, band, *filter;
+ complex kzp, ckp, ckp2;
+
+ kp = om/cp;
+ kp2 = kp*kp;
+ dkx = 2.0*M_PI/(nkx*dx);
+ stab = eps*eps*kp*kp;
+
+ /* 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;
+
+ 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;
+ kzp = froot(kzp2);
+
+ pu[ikx].r = filter[ikx]*kzp.r/(om*rho);
+ pu[ikx].i = filter[ikx]*kzp.i/(om*rho);
+ }
+
+ /* operators are symmetric in kx-w domain */
+ for (ikx = (nkx/2+1); ikx < nkx; ikx++) {
+ pu[ikx] = pu[nkx-ikx];
+ }
+ free(filter);
+
+ return;
+}
+
+
+complex froot(float x)
+{
+ complex z;
+ if (x >= 0.0) {
+ z.r = sqrt(x);
+ z.i = 0.0;
+ return z;
+ }
+ else {
+ z.r = 0.0;
+ z.i = -sqrt(-x);
+ return z;
+ }
+}
+
+
/* compute 1/x */
complex firoot(float x, float stab)
{
diff --git a/fdelmodc/demo/decomposition.scr b/fdelmodc/demo/decomposition.scr
index 5077fb0..384422d 100755
--- a/fdelmodc/demo/decomposition.scr
+++ b/fdelmodc/demo/decomposition.scr
@@ -27,13 +27,13 @@ makemod sizex=6000 sizez=4000 dx=$dx dz=$dx cp0=$cp cs0=$cs ro0=$rho \
makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-
+# pressure Normalised decompostion
../fdelmodc \
file_cp=syncl_cp.su ischeme=1 \
file_den=syncl_ro.su \
file_src=wave.su \
file_rcv=shot_fd_Fz_zsrc1150.su \
- nshot=2 dxshot=500 \
+ nshot=1 dxshot=500 \
src_type=7 \
src_orient=1 \
src_injectionrate=1 \
@@ -46,8 +46,8 @@ makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
xrcv1=-2250 xrcv2=2250 \
zrcv1=1150 zrcv2=1150 \
xsrc=0 zsrc=0 \
- ntaper=200 \
- left=4 right=4 top=4 bottom=4
+ ntaper=50 \
+ left=2 right=2 top=2 bottom=2
supsgraph < anglerp.su hbox=2 wbox=4 style=normal \
@@ -70,3 +70,26 @@ supsimage < shot_fd_Fz_zsrc1150_rd.su \
wbox=2 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
label1="time [s]" label2="lateral position [m]" > shot_fd_Fz_zsrc1150_rd.eps
+# Particle Velocity Normalised decompostion
+../fdelmodc \
+ file_cp=syncl_cp.su ischeme=1 \
+ file_den=syncl_ro.su \
+ file_src=wave.su \
+ file_rcv=shot_fd_Fz_zsrc1150_vznorm.su \
+ nshot=1 dxshot=500 \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_ud=2 \
+ dtrcv=0.0040 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=3.10 \
+ dxrcv=10 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=1150 zrcv2=1150 \
+ xsrc=0 zsrc=0 \
+ ntaper=50 \
+ left=2 right=2 top=2 bottom=2
+
+
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 66b19a5..ebbf8d7 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -223,7 +223,8 @@ char *sdoc[] = {
" rec_type_txz=0 .... Txz registration _rtxz",
" rec_type_pp=0 ..... P (divergence) registration _rP",
" rec_type_ss=0 ..... S (curl) registration _rS",
-" rec_type_ud=0 ..... decomposition in up and downgoing waves _ru, _rd",
+" 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",
" kangle= ........... maximum wavenumber angle for decomposition",
" rec_int_vx=0 ..... interpolation of Vx receivers",
" - 0=Vx->Vx (no interpolation)",
diff --git a/fdelmodc/writeRec.c b/fdelmodc/writeRec.c
index 167e5f7..6c2ce08 100644
--- a/fdelmodc/writeRec.c
+++ b/fdelmodc/writeRec.c
@@ -34,7 +34,7 @@ int traceWrite(segy *hdr, float *data, int n, FILE *fp) ;
void name_ext(char *filename, char *extension);
void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
int nkx, float dx, int nt, float dt, float fmin, float fmax,
- float cp, float rho, int verbose);
+ float cp, float rho, int vznorm, int verbose);
#define MAX(x,y) ((x) > (y) ? (x) : (y))
@@ -50,7 +50,7 @@ int writeRec(recPar rec, modPar mod, bndPar bnd, wavPar wav, int ixsrc, int izsr
float dx, dt, cp, rho, fmin, fmax;
complex *crec_vz, *crec_p, *crec_up, *crec_dw;
int irec, ntfft, nfreq, nkx, xorig, ix, iz, it, ibndx;
- int append;
+ int append, vznorm;
double ddt;
char number[16], filename[1024];
segy hdr;
@@ -117,6 +117,8 @@ 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 (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));
@@ -140,7 +142,7 @@ int writeRec(recPar rec, modPar mod, bndPar bnd, wavPar wav, int ixsrc, int izsr
/* apply decomposition operators */
kxwdecomp(crec_p, crec_vz, crec_up, crec_dw,
- nkx, mod.dx, nsam, dt, fmin, fmax, cp, rho, verbose);
+ nkx, mod.dx, nsam, dt, fmin, fmax, cp, rho, vznorm, verbose);
/* transform back to t-x */
wkx2xt(crec_up, rec_up, ntfft, nkx, nkx, ntfft, xorig);
--
GitLab