-
Jan Thorbecke authored
removing unused variables in all code, bug fixes in marchenko to handle different acquisition geometries, re-organizing calls within marchenko algorithm in preperation for software release
Jan Thorbecke authoredremoving unused variables in all code, bug fixes in marchenko to handle different acquisition geometries, re-organizing calls within marchenko algorithm in preperation for software release
getRecTimes.c 11.77 KiB
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<assert.h>
#include"fdelmodc.h"
#include"par.h"
/**
* Stores the wavefield at the receiver positions.
*
* On a staggered grid the fields are all on different positions,
* to compensate for that the rec.int_vx and rec.int_vz options
* can be set.
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *vx, float *vz, float *tzz, float *txx, float *txz, float *l2m, float *rox, float *roz, float *rec_vx, float *rec_vz, float *rec_txx, float *rec_tzz, float *rec_txz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose)
{
int n1, ibndx, ibndz;
int irec, ix, iz, ix2, iz2, ix1, iz1;
float dvx, dvz, rdz, rdx, C00, C10, C01, C11;
float *vz_t, c1, c2, lroz, field;
ibndx = mod.ioPx;
ibndz = mod.ioPz;
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
n1 = mod.naz;
c1 = 9.0/8.0;
c2 = -1.0/24.0;
if (!rec.n) return 0;
/***********************************************************************
* velocity or txz or potential registrations issues:
* rec_x and rec_z are related to actual txx/tzz/p positions.
* offsets from virtual boundaries must be taken into account.
*
* 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
*
* Note, in the acoustic scheme P is stored in the Tzz array.
***********************************************************************/
for (irec=0; irec<rec.n; irec++) {
iz = rec.z[irec]+ibndz;
ix = rec.x[irec]+ibndx;
iz1 = iz-1;
ix1 = ix-1;
iz2 = iz+1;
ix2 = ix+1;
/* interpolation to precise (not necessary on a grid point) position */
if ( rec.int_p==3 ) {
iz = (int)floorf(rec.zr[irec]/mod.dz)+ibndz;
ix = (int)floorf(rec.xr[irec]/mod.dx)+ibndx;
rdz = (rec.zr[irec] - (iz-ibndz)*mod.dz)/mod.dz;
rdx = (rec.xr[irec] - (ix-ibndx)*mod.dx)/mod.dx;
iz1 = iz-1;
ix1 = ix-1;
iz2 = iz+1;
ix2 = ix+1;
/*
// Interpolate according to Dirk Kraaijpool's scheme
// Reference: "Seismic ray fields and ray field maps : theory and algorithms" , // PhD thesis Utrecht University,Faculty of Geosciences, 2003)
C00 = tzz[ix*n1+iz] + 0.5*((tzz[(ix+1)*n1+iz] +tzz[(ix-1)*n1+iz]+
tzz[(ix )*n1+iz+1] +tzz[(ix )*n1+iz-1])/(2.0*mod.dx));
C10 = tzz[(ix+1)*n1+iz] + 0.5*((tzz[(ix+2)*n1+iz] +tzz[(ix )*n1+iz]+
tzz[(ix+1)*n1+iz+1] +tzz[(ix+1)*n1+iz-1])/(2.0*mod.dz));
C01 = tzz[ix*n1+iz+1] + 0.5*((tzz[(ix+1)*n1+iz+1] +tzz[(ix-1)*n1+iz+1]+
tzz[(ix)*n1+iz+2] +tzz[(ix )*n1+iz])/(2.0*mod.dx));
C11 = tzz[(ix+1)*n1+iz+1]+ 0.5*((tzz[(ix+2)*n1+iz+1] +tzz[(ix )*n1+iz+1]+
tzz[(ix+1)*n1+iz+2] +tzz[(ix+1)*n1+iz])/(2.0*mod.dz));
*/
if (rec.type.p){
/* bi-linear interpolation */
C00 = tzz[ix*n1+iz];
C10 = tzz[(ix+1)*n1+iz];
C01 = tzz[ix*n1+iz+1];
C11 = tzz[(ix+1)*n1+iz+1];
rec_p[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.txx) {
C00 = txx[ix*n1+iz];
C10 = txx[(ix+1)*n1+iz];
C01 = txx[ix*n1+iz+1];
C11 = txx[(ix+1)*n1+iz+1];
rec_txx[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.tzz) {
C00 = tzz[ix*n1+iz];
C10 = tzz[(ix+1)*n1+iz];
C01 = tzz[ix*n1+iz+1];
C11 = tzz[(ix+1)*n1+iz+1];
rec_tzz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.txz) {
C00 = txz[ix2*n1+iz2];
C10 = txz[(ix2+1)*n1+iz2];
C01 = txz[ix2*n1+iz2+1];
C11 = txz[(ix2+1)*n1+iz2+1];
rec_txz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.pp) {
C00 = (vx[ix2*n1+iz]-vx[ix*n1+iz] +
vz[ix*n1+iz2]-vz[ix*n1+iz])/mod.dx;
C10 = (vx[(ix2+1)*n1+iz]-vx[(ix+1)*n1+iz] +
vz[(ix+1)*n1+iz2]-vz[(ix+1)*n1+iz])/mod.dx;
C01 = (vx[ix2*n1+iz+1]-vx[ix*n1+iz+1] +
vz[ix*n1+iz2+1]-vz[ix*n1+iz+1])/mod.dx;
C11 = (vx[(ix2+1)*n1+iz+1]-vx[(ix+1)*n1+iz+1] +
vz[(ix+1)*n1+iz2+1]-vz[(ix+1)*n1+iz+1])/mod.dx;
rec_pp[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.ss) {
C00 = (vx[ix2*n1+iz2]-vx[ix2*n1+iz] -
(vz[ix2*n1+iz2]-vz[ix*n1+iz2]))/mod.dx;
C10 = (vx[(ix2+1)*n1+iz2]-vx[(ix2+1)*n1+iz] -
(vz[(ix2+1)*n1+iz2]-vz[(ix+1)*n1+iz2]))/mod.dx;
C01 = (vx[ix2*n1+iz2+1]-vx[ix2*n1+iz+1] -
(vz[ix2*n1+iz2+1]-vz[ix*n1+iz2+1]))/mod.dx;;
C11 = (vx[(ix2+1)*n1+iz2+1]-vx[(ix2+1)*n1+iz+1] -
(vz[(ix2+1)*n1+iz2+1]-vz[(ix+1)*n1+iz2+1]))/mod.dx;
rec_ss[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.vz) { C00 = vz[ix*n1+iz2];
C10 = vz[(ix+1)*n1+iz2];
C01 = vz[ix*n1+iz2+1];
C11 = vz[(ix+1)*n1+iz2+1];
rec_vz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
if (rec.type.vx) {
C00 = vx[ix2*n1+iz];
C10 = vx[(ix2+1)*n1+iz];
C01 = vx[ix2*n1+iz+1];
C11 = vx[(ix2+1)*n1+iz+1];
rec_vx[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
C01*(1.0-rdx)*rdz + C11*rdx*rdz;
}
}
else { /* read values directly from the grid points */
if (verbose>=4 && isam==0) {
vmess("Receiver %d read at gridpoint ix=%d iz=%d",irec, ix, iz);
}
/* interpolation of receivers to same time step is only done for acoustic scheme */
if (rec.type.p) {
if (rec.int_p == 1) {
if (mod.ischeme == 1) { /* interpolate Tzz times -1/2 Dt backward to Vz times */
dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
field = tzz[ix*n1+iz] + 0.5*l2m[ix*n1+iz]*(dvx+dvz);
dvx = c1*(vx[(ix+1)*n1+iz1] - vx[ix*n1+iz1]) +
c2*(vx[(ix+2)*n1+iz1] - vx[(ix-1)*n1+iz1]);
dvz = c1*(vz[ix*n1+iz1+1] - vz[ix*n1+iz1]) +
c2*(vz[ix*n1+iz1+2] - vz[ix*n1+iz1-1]);
field += tzz[ix*n1+iz1] + 0.5*l2m[ix*n1+iz1]*(dvx+dvz);
rec_p[irec*rec.nt+isam] = 0.5*field;
}
else {
rec_p[irec*rec.nt+isam] = 0.5*(tzz[ix*n1+iz1]+tzz[ix*n1+iz]);
}
}
else if (rec.int_p == 2) {
if (mod.ischeme == 1) { /* interpolate Tzz times -1/2 Dt backward to Vx times */
dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
field = tzz[ix*n1+iz] + 0.5*l2m[ix*n1+iz]*(dvx+dvz);
dvx = c1*(vx[(ix1+1)*n1+iz] - vx[ix1*n1+iz]) +
c2*(vx[(ix+2)*n1+iz] - vx[(ix1-1)*n1+iz]);
dvz = c1*(vz[ix1*n1+iz+1] - vz[ix1*n1+iz]) +
c2*(vz[ix1*n1+iz+2] - vz[ix1*n1+iz-1]);
field += tzz[ix1*n1+iz] + 0.5*l2m[ix1*n1+iz]*(dvx+dvz);
rec_p[irec*rec.nt+isam] = 0.5*field;
}
else {
rec_p[irec*rec.nt+isam] = 0.5*(tzz[ix1*n1+iz]+tzz[ix*n1+iz]);
}
}
else {
rec_p[irec*rec.nt+isam] = tzz[ix*n1+iz];
}
}
if (rec.type.txx) rec_txx[irec*rec.nt+isam] = txx[ix*n1+iz];
if (rec.type.tzz) rec_tzz[irec*rec.nt+isam] = tzz[ix*n1+iz];
if (rec.type.txz) { /* time interpolation to be done */
if (rec.int_vz == 2 || rec.int_vx == 2) {
rec_txz[irec*rec.nt+isam] = 0.25*(
txz[ix*n1+iz2]+txz[ix2*n1+iz2]+
txz[ix*n1+iz]+txz[ix2*n1+iz]); }
else {
rec_txz[irec*rec.nt+isam] = txz[ix2*n1+iz2];
}
}
if (rec.type.pp) {
rec_pp[irec*rec.nt+isam] = (vx[ix2*n1+iz]-vx[ix*n1+iz] +
vz[ix*n1+iz2]-vz[ix*n1+iz])/mod.dx;
}
if (rec.type.ss) {
rec_ss[irec*rec.nt+isam] = (vx[ix2*n1+iz2]-vx[ix2*n1+iz] -
(vz[ix2*n1+iz2]-vz[ix*n1+iz2]))/mod.dx;
}
if (rec.type.vz) {
/* interpolate vz to vx position to the right and above of vz */
if (rec.int_vz == 1) {
rec_vz[irec*rec.nt+isam] = 0.25*(
vz[ix*n1+iz2]+vz[ix1*n1+iz2]+
vz[ix*n1+iz] +vz[ix1*n1+iz]);
}
/* interpolate vz to Txx/Tzz position by taking the mean of 2 values */
else if (rec.int_vz == 2) {
if (mod.ischeme == 1) { /* interpolate Vz times +1/2 Dt forward to P times */
field = vz[ix*n1+iz] - 0.5*roz[ix*n1+iz]*(
c1*(tzz[ix*n1+iz] - tzz[ix*n1+iz-1]) +
c2*(tzz[ix*n1+iz+1] - tzz[ix*n1+iz-2]));
field += vz[ix*n1+iz2] - 0.5*roz[ix*n1+iz2]*(
c1*(tzz[ix*n1+iz2] - tzz[ix*n1+iz2-1]) +
c2*(tzz[ix*n1+iz2+1] - tzz[ix*n1+iz2-2]));
rec_vz[irec*rec.nt+isam] = 0.5*field;
}
else {
rec_vz[irec*rec.nt+isam] = 0.5*(vz[ix*n1+iz2]+vz[ix*n1+iz]);
}
}
else {
rec_vz[irec*rec.nt+isam] = vz[ix*n1+iz2];
//rec_vz[irec*rec.nt+isam] = vz[ix*n1+iz];
//fprintf(stderr,"isam=%d vz[%d]=%e vz[%d]=%e vz[%d]=%e \n",isam, iz-1,vz[ix*n1+iz-1],iz,vz[ix*n1+iz], iz+1, vz[ix*n1+iz+1]);
}
}
if (rec.type.vx) {
/* interpolate vx to vz position to the left and below of vx */
if (rec.int_vx == 1) {
rec_vx[irec*rec.nt+isam] = 0.25*(
vx[ix2*n1+iz]+vx[ix2*n1+iz1]+
vx[ix*n1+iz]+vx[ix*n1+iz1]);
}
/* interpolate vx to Txx/Tzz position by taking the mean of 2 values */
else if (rec.int_vx == 2) {
if (mod.ischeme == 1) { /* interpolate Vx times +1/2 Dt forward to P times */
field = vx[ix*n1+iz] - 0.5*rox[ix*n1+iz]*(
c1*(tzz[ix*n1+iz] - tzz[(ix-1)*n1+iz]) +
c2*(tzz[(ix+1)*n1+iz] - tzz[(ix-2)*n1+iz]));
field += vx[ix2*n1+iz] - 0.5*rox[ix2*n1+iz]*(
c1*(tzz[ix2*n1+iz] - tzz[(ix2-1)*n1+iz]) +
c2*(tzz[(ix2+1)*n1+iz] - tzz[(ix2-2)*n1+iz]));
rec_vx[irec*rec.nt+isam] = 0.5*field;
}
else {
rec_vx[irec*rec.nt+isam] = 0.5*(vx[ix2*n1+iz]+vx[ix*n1+iz]);
}
}
else {
rec_vx[irec*rec.nt+isam] = vx[ix2*n1+iz];
}
}
}
} /* end of irec loop */
/* store all x-values on z-level for P Vz for up-down decomposition */
if (rec.type.ud) {
iz = rec.z[0]+ibndz;
iz2 = iz+1;
vz_t = (float *)calloc(2*mod.nax,sizeof(float));
/* P and Vz are staggered in time and need to correct for this */
/* -1- compute Vz at next time step and average with current time step */
lroz = mod.dt/(mod.dx*rec.rho);
for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
vz_t[ix] = vz[ix*n1+iz] - lroz*(
c1*(tzz[ix*n1+iz] - tzz[ix*n1+iz-1]) +
c2*(tzz[ix*n1+iz+1] - tzz[ix*n1+iz-2]));
vz_t[mod.nax+ix] = vz[ix*n1+iz2] - lroz*(
c1*(tzz[ix*n1+iz2] - tzz[ix*n1+iz2-1]) +
c2*(tzz[ix*n1+iz2+1] - tzz[ix*n1+iz2-2]));
}
for (ix=0; ix<mod.nax; ix++) {
/* -2- compute average in time and depth to get Vz at same depth and time as P */
rec_udvz[ix*rec.nt+isam] = 0.25*(vz[ix*n1+iz2]+vz[ix*n1+iz]+vz_t[mod.nax+ix]+vz_t[ix]);
rec_udp[ix*rec.nt+isam] = tzz[ix*n1+iz];
}
free(vz_t);
}
return 0;
}