Skip to content
Snippets Groups Projects
Commit 3f07253d authored by Jan Thorbecke's avatar Jan Thorbecke
Browse files

cleaning

parent 972852ff
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fdelmodc/decomposition.c.new deleted 100644 → 0
+ 0
− 247
View file @ 972852ff
/*
* decomposition.c
*
* Kees Wapenaar "Reciprocity properties of one-way propagators"
* GEOPHYSICS, VOL. 63, NO. 4 (JULY-AUGUST 1998); P. 1795–1798
*
* Created by Jan Thorbecke on 17/03/2014.
* Copyright 2014 TU Delft. All rights reserved.
*
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include <string.h>
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
float r,i;
} complex;
typedef struct _dcomplexStruct { /* complex number */
double r,i;
} dcomplex;
#endif/* complex */
complex firoot(float x, float stab)
complex ciroot(complex x, float stab);
complex cwp_csqrt(complex z);
void decud(float om, float rho, float cp, float dx, int nkx, complex *pu);
void kxwfilter(complex *data, float k, float dx, int nkx,
float alfa1, float alfa2, float perc);
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 iom, iomin, iomax, ikx, nfreq;
float omin, omax, deltom, om, kp, df;
float alpha, eps;
complex *pu, w;
complex ax, az;
df = 1.0/((float)nt*dt);
deltom = 2.*M_PI*df;
omin = 2.*M_PI*fmin;
omax = 2.*M_PI*fmax;
nfreq = nt/2+1;
eps = 0.01;
alpha = 0.1;
iomin = (int)MIN((omin/deltom), (nfreq-1));
iomin = MAX(iomin, 1);
iomax = MIN((int)(omax/deltom), (nfreq-1));
pu = (complex *)malloc(nkx*sizeof(complex));
for (iom = iomin; iom <= iomax; iom++) {
om = iom*deltom;
decud(om, rho, cp, dx, nkx, 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;
}
}
free(pu);
return;
}
void decud(float om, float rho, float cp, float dx, int nkx, float alpha, float eps, complex *pu);
{
int ikx, ikxmax1, ikxmax2, filterpoints, filterppos;
float mu, kp, kp2, ks, ks2, ksk;
float kx, kx2, kzp2, kzs2, dkx;
float kxfmax, kxnyq, kpos, kneg, alfa, kfilt, perc, band, *filter;
complex kzp, kzs, cste, ckp, ckp2, ckzp2;
/* with complex frequency
wom.r=om;
wom.i=alpha;
ckp.r = wom.r/cp;
ckp.i = wom.i/cp;
ckp2.r = ckp.r*ckp.r-ckp.i*ckp.i;
ckp2.i = 2.0*ckp.r*ckp.i;
stab = eps*eps*(ckp.r*ckp.r+ckp.i*ckp.i);
*/
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 */
alfa = 90.0;
perc = 0.10; /* percentage of band to use for smooth filter */
filter = (float *)malloc(nkx*sizeof(float));
kpos = kp*sin(M_PI*alfa/180.0);
kneg = -kpos;
kxnyq = M_PI/dx;
if (kpos > kxnyq) kpos = kxnyq;
band = kpos;
filterpoints = (int)fabs((int)(perc*band/dkx));
kfilt = fabs(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;
kzp = firoot(kzp2);
/* with complex frequency
kzp2.r = kp2.r - kx2;
kzp2.i = kp2.i;
kzp = ciroot(kzp2, stab);
*/
if (kzp2 != 0) {
pu[ikx].r = filter[ikx]*om*rho*kzp.r;
pu[ikx].i = filter[ikx]*om*rho*kzp.i;
// pu[ikx].r = om*rho*kzp.r;
// pu[ikx].i = om*rho*kzp.i;
}
else {
pu[ikx].r = 0.0;
pu[ikx].i = 0.0;
}
}
/* operators are symmetric in kx-w domain */
for (ikx = (nkx/2+1); ikx < nkx; ikx++) {
pu[ikx] = pu[nkx-ikx];
}
free(filter);
return;
}
/* compute 1/x */
complex firoot(float x, float stab)
{
complex z;
if (x > 0.0) {
z.r = 1.0/sqrt(x+stab);
z.i = 0.0;
}
else if (x == 0.0) {
z.r = 0.0;
z.i = 0.0;
}
else {
z.r = 0.0;
z.i = 1.0/sqrt(-x+stab);
}
return z;
}
complex ciroot(complex x, float stab)
{
complex z, kz, kzz;
float kd;
if (x.r == 0.0) {
z.r = 0.0;
z.i = 0.0;
}
else {
kzz = cwp_csqrt(x);
kz.r = kzz.r;
kz.i = -abs(kzz.i);
kd = kz.r*kz.r+kz.i*kz.i+stab;
z.r = kz.r/kd;
z.i = -kz.i/kd;
}
return z;
}
complex cwp_csqrt(complex z)
{
complex c;
float x,y,w,r;
if (z.r==0.0 && z.i==0.0) {
c.r = c.i = 0.0;
return c;
} else {
x = fabs(z.r);
y = fabs(z.i);
if (x>=y) {
r = y/x;
w = sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
} else {
r = x/y;
w = sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
}
if (z.r>=0.0) {
c.r = w;
c.i = z.i/(2.0*w);
} else {
c.i = (z.i>=0.0) ? w : -w;
c.r = z.i/(2.0*c.i);
}
return c;
}
}
fdelmodc/writeSnapTimes.c.org deleted 100644 → 0
+ 0
− 172
View file @ 972852ff
#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 <string.h>
#include "par.h"
#include "segy.h"
#include "fdelmodc.h"
/**
* Writes gridded wavefield(s) at a desired time to output file(s)
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
FILE *fileOpen(char *file, char *ext, int append);
int traceWrite(segy *hdr, float *data, int n, FILE *fp);
#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))
int writeSnapTimes(modPar mod, snaPar sna, int ixsrc, int izsrc, int itime, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose)
{
FILE *fpvx, *fpvz, *fptxx, *fptzz, *fptxz, *fpp, *fppp, *fpss;
int append, isnap;
int n1, ibnd, ixs, izs, ize, i, j;
int ix, iz, ix2, iz2;
float *snap, sdx;
segy hdr;
if (sna.nsnap==0) return 0;
ibnd = mod.iorder/2-1;
n1 = mod.naz;
sdx = 1.0/mod.dx;
/* check if this itime is a desired snapshot time */
if ( (((itime-sna.delay) % sna.skipdt)==0) &&
(itime >= sna.delay) &&
(itime <= sna.delay+(sna.nsnap-1)*sna.skipdt) ) {
isnap = NINT((itime-sna.delay)/sna.skipdt);
if (verbose) vmess("Writing snapshot(%d) at time=%.3f", isnap+1, itime*mod.dt);
if (isnap) append=1;
else append=0;
if (sna.type.vx) fpvx = fileOpen(sna.file_snap, "_svx", append);
if (sna.type.vz) fpvz = fileOpen(sna.file_snap, "_svz", append);
if (sna.type.p) fpp = fileOpen(sna.file_snap, "_sp", append);
if (sna.type.txx) fptxx = fileOpen(sna.file_snap, "_stxx", append);
if (sna.type.tzz) fptzz = fileOpen(sna.file_snap, "_stzz", append);
if (sna.type.txz) fptxz = fileOpen(sna.file_snap, "_stxz", append);
if (sna.type.pp) fppp = fileOpen(sna.file_snap, "_spp", append);
if (sna.type.ss) fpss = fileOpen(sna.file_snap, "_sss", append);
memset(&hdr,0,TRCBYTES);
hdr.dt = 1000000*(mod.dt);
hdr.scalco = -1000;
hdr.scalel = -1000;
hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
hdr.fldr = isnap+1;
hdr.trid = 1;
hdr.ns = sna.nz;
hdr.trwf = sna.nx;
hdr.ntr = (isnap+1)*sna.nx;
hdr.f1 = sna.z1*mod.dz+mod.z0;
hdr.f2 = sna.x1*mod.dx+mod.x0;
hdr.d1 = mod.dz*sna.skipdz;
hdr.d2 = mod.dx*sna.skipdx;
/***********************************************************************
* vx velocities have one sample less in x-direction
* vz velocities have one sample less in z-direction
* txz stresses have one sample less in z-direction and x-direction
***********************************************************************/
snap = (float *)malloc(sna.nz*sizeof(float));
/* Decimate, with skipdx and skipdz, the number of gridpoints written to file
and write to file. */
for (ixs=sna.x1, i=0; ixs<=sna.x2; ixs+=sna.skipdx, i++) {
hdr.tracf = i+1;
hdr.tracl = isnap*sna.nx+i+1;
hdr.gx = 1000*(mod.x0+ixs*mod.dx);
ix = ixs+ibnd;
ix2 = ix+1;
izs = sna.z1+ibnd;
ize = sna.z2+ibnd;
if (sna.type.vx) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = vx[ix2*n1+iz];
}
traceWrite(&hdr, snap, sna.nz, fpvx);
}
if (sna.type.vz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = vz[ix*n1+iz+1];
}
traceWrite(&hdr, snap, sna.nz, fpvz);
}
if (sna.type.p) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = tzz[ix*n1+iz];
}
traceWrite(&hdr, snap, sna.nz, fpp);
}
if (sna.type.tzz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = tzz[ix*n1+iz];
}
traceWrite(&hdr, snap, sna.nz, fptzz);
}
if (sna.type.txx) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = txx[ix*n1+iz];
}
traceWrite(&hdr, snap, sna.nz, fptxx);
}
if (sna.type.txz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = txz[ix2*n1+iz+1];
}
traceWrite(&hdr, snap, sna.nz, fptxz);
}
/* calculate divergence of velocity field */
if (sna.type.pp) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
iz2 = iz+1;
snap[j] = sdx*((vx[ix2*n1+iz]-vx[ix*n1+iz])+
(vz[ix*n1+iz2]-vz[ix*n1+iz]));
}
traceWrite(&hdr, snap, sna.nz, fppp);
}
/* calculate rotation of velocity field */
if (sna.type.ss) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
iz2 = iz+1;
snap[j] = sdx*((vx[ix2*n1+iz2]-vx[ix2*n1+iz])-
(vz[ix2*n1+iz2]-vz[ix*n1+iz2]));
}
traceWrite(&hdr, snap, sna.nz, fpss);
}
}
if (sna.type.vx) fclose(fpvx);
if (sna.type.vz) fclose(fpvz);
if (sna.type.p) fclose(fpp);
if (sna.type.txx) fclose(fptxx);
if (sna.type.tzz) fclose(fptzz);
if (sna.type.txz) fclose(fptxz);
if (sna.type.pp) fclose(fppp);
if (sna.type.ss) fclose(fpss);
free(snap);
}
return 0;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment