diff --git a/fdelmodc3D/.vscode/settings.json b/fdelmodc3D/.vscode/settings.json
new file mode 100644
index 0000000000000000000000000000000000000000..871aedf628d81ee20a28f494175a0b13b0d82556
--- /dev/null
+++ b/fdelmodc3D/.vscode/settings.json
@@ -0,0 +1,5 @@
+{
+ "files.associations": {
+ "*.tcc": "c"
+ }
+}
\ No newline at end of file
diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c
index 8ba30a7c0869efe1884b033f0acbbd5ecbcc1454..1bb3ed3fc890c62e4548ad1867f72114d52647a4 100644
--- a/fdelmodc3D/fdelmodc3D.c
+++ b/fdelmodc3D/fdelmodc3D.c
@@ -639,29 +639,33 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
/* Note that time step it=0 (t=0 for t**-fields t=-1/2 dt for v*-field) is not recorded */
isam = (it-rec.delay-itwritten)/rec.skipdt+1;
/* store time at receiver positions */
- getRecTimes(mod, rec, bnd, it, isam, vx, vz, tzz, txx, txz,
- l2m, rox, roz,
- rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
+ getRecTimes3D(mod, rec, bnd, it, isam, vx, vy, vz,
+ tzz, tyy, txx, txz, txy, tyz,
+ l2m, rox, roy, roz,
+ rec_vx, rec_vy, rec_vz,
+ rec_txx, rec_tyy, rec_tzz, rec_txz, rec_txy, rec_tyz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
/* at the end of modeling a shot, write receiver array to output file(s) */
if (writeToFile && (it+rec.skipdt <= it1-1) ) {
fileno = ( ((it-rec.delay)/rec.skipdt)+1)/rec.nt;
- writeRec(rec, mod, bnd, wav, ixsrc, izsrc, isam+1, ishot, fileno,
- rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
- rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
+ writeRec3D(rec, mod, bnd, wav, ixsrc, iysrc, izsrc, isam+1, ishot, fileno,
+ rec_vx, rec_vy, rec_vz, rec_txx, rec_tyy, rec_tzz, rec_txz, rec_tyz, rec_txy,
+ rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
}
}
/* write snapshots to output file(s) */
if (sna.nsnap) {
- writeSnapTimes(mod, sna, bnd, wav, ixsrc, izsrc, it, vx, vz, tzz, txx, txz, verbose);
+ writeSnapTimes3D(mod, sna, bnd, wav, ixsrc, iysrc, izsrc, it,
+ vx, vy, vz, tzz, tyy, txx, txz, tyz, txy, verbose);
+
}
/* calculate beams */
if(sna.beam) {
- getBeamTimes(mod, sna, vx, vz, tzz, txx, txz,
- beam_vx, beam_vz, beam_txx, beam_tzz, beam_txz,
+ getBeamTimes3D(mod, sna, vx, vy, vz, tzz, tyy, txx, txz, tyz, txy,
+ beam_vx, beam_vy, beam_vz, beam_txx, beam_tyy, beam_tzz, beam_txz, beam_tyz, beam_txy,
beam_p, beam_pp, beam_ss, verbose);
}
}
@@ -686,27 +690,29 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
if (fileno) fileno++;
if (rec.scale==1) { /* scale receiver with distance src-rcv */
- float xsrc, zsrc, Rrec, rdx, rdz;
+ float xsrc, ysrc, zsrc, Rrec, rdx, rdy, rdz;
long irec;
xsrc=mod.x0+mod.dx*ixsrc;
+ ysrc=mod.y0+mod.dy*iysrc;
zsrc=mod.z0+mod.dz*izsrc;
for (irec=0; irec<rec.n; irec++) {
rdx=mod.x0+rec.xr[irec]-xsrc;
+ rdy=mod.y0+rec.yr[irec]-ysrc;
rdz=mod.z0+rec.zr[irec]-zsrc;
- Rrec = sqrt(rdx*rdx+rdz*rdz);
- fprintf(stderr,"Rec %li is scaled with distance %f R=%.2f,%.2f S=%.2f,%.2f\n", irec, Rrec,rdx,rdz,xsrc,zsrc);
+ Rrec = sqrt(rdx*rdx+rdy*rdy+rdz*rdz);
+ fprintf(stderr,"Rec %li is scaled with distance %f R=%.2f,%.2f,%.2f S=%.2f,%.2f,%.2f\n", irec, Rrec,rdx,rdy,rdz,xsrc,ysrc,zsrc);
for (it=0; it<rec.nt; it++) {
rec_p[irec*rec.nt+it] *= sqrt(Rrec);
}
}
}
- writeRec(rec, mod, bnd, wav, ixsrc, izsrc, isam+1, ishot, fileno,
- rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
- rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
+ writeRec3D(rec, mod, bnd, wav, ixsrc, iysrc, izsrc, isam+1, ishot, fileno,
+ rec_vx, rec_vy, rec_vz, rec_txx, rec_tyy, rec_tzz, rec_txz, rec_tyz, rec_txy,
+ rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
- writeBeams(mod, sna, ixsrc, izsrc, ishot, fileno,
- beam_vx, beam_vz, beam_txx, beam_tzz, beam_txz,
- beam_p, beam_pp, beam_ss, verbose);
+ writeBeams3D(mod, sna, ixsrc, iysrc, izsrc, ishot, fileno,
+ beam_vx, beam_vy, beam_vz, beam_txx, beam_tyy, beam_tzz, beam_txz, beam_tyz, beam_txy,
+ beam_p, beam_pp, beam_ss, verbose);
} /* end of loop over number of shots */
@@ -721,20 +727,26 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
initargs(argc,argv); /* this will free the arg arrays declared */
free(rox);
+ free(roy);
free(roz);
free(l2m);
free(src_nwav[0]);
free(src_nwav);
free(vx);
+ free(vy);
free(vz);
free(tzz);
freeStoreSourceOnSurface();
if (rec.type.vz) free(rec_vz);
+ if (rec.type.vy) free(rec_vy);
if (rec.type.vx) free(rec_vx);
if (rec.type.p) free(rec_p);
if (rec.type.txx) free(rec_txx);
+ if (rec.type.tyy) free(rec_tyy);
if (rec.type.tzz) free(rec_tzz);
if (rec.type.txz) free(rec_txz);
+ if (rec.type.tyz) free(rec_tyz);
+ if (rec.type.txy) free(rec_txy);
if (rec.type.pp) free(rec_pp);
if (rec.type.ss) free(rec_ss);
if (rec.type.ud) {
@@ -743,11 +755,15 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
}
if(sna.beam) {
if (sna.type.vz) free(beam_vz);
+ if (sna.type.vy) free(beam_vy);
if (sna.type.vx) free(beam_vx);
if (sna.type.p) free(beam_p);
if (sna.type.txx) free(beam_txx);
+ if (sna.type.tyy) free(beam_tyy);
if (sna.type.tzz) free(beam_tzz);
if (sna.type.txz) free(beam_txz);
+ if (sna.type.tyz) free(beam_tyz);
+ if (sna.type.txy) free(beam_txy);
if (sna.type.pp) free(beam_pp);
if (sna.type.ss) free(beam_ss);
}
@@ -761,7 +777,10 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
free(lam);
free(mul);
free(txz);
+ free(tyz);
+ free(txy);
free(txx);
+ free(tyy);
}
if (mod.ischeme==4) {
free(tss);
diff --git a/fdelmodc3D/getBeamTimes3D.c b/fdelmodc3D/getBeamTimes3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..f23ef7b7a4e4f163cdef3dc4b8d367b06e910743
--- /dev/null
+++ b/fdelmodc3D/getBeamTimes3D.c
@@ -0,0 +1,246 @@
+#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 "segy.h"
+#include "fdelmodc3D.h"
+
+/**
+* getBeamTimes: stores energy fields (beams) in arrays at certain time steps
+* writeBeams: writes the stored fields 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);
+void name_ext(char *filename, char *extension);
+void vmess(char *fmt, ...);
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+long getBeamTimes3D(modPar mod, snaPar sna, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *tyz, float *txy,
+ float *beam_vx, float *beam_vy, float *beam_vz, float *beam_txx, float *beam_tyy, float *beam_tzz, float *beam_txz, float *beam_tyz, float *beam_txy,
+ float *beam_p, float *beam_pp, float *beam_ss, long verbose)
+{
+ long n1, n2, ibndx, ibndy, ibndz, ixs, iys, izs, ize, i, j, l;
+ long ix, iy, iz, ix2, iy2, iz2;
+ float sdx, s, p;
+
+ ibndx = mod.ioPx;
+ ibndy = mod.ioPy;
+ ibndz = mod.ioPz;
+ n1 = mod.naz;
+ n2 = mod.nax;
+ sdx = 1.0/mod.dx;
+ izs = sna.z1+ibndx;
+ ize = sna.z2+ibndz;
+
+ for (iys=sna.y1, l=0; iys<=sna.y2; iys+=sna.skipdy, l++) {
+ iy = iys+ibndy;
+ iy2 = iy+1;
+ for (ixs=sna.x1, i=0; ixs<=sna.x2; ixs+=sna.skipdx, i++) {
+ ix = ixs+ibndx;
+ ix2 = ix+1;
+
+ if (sna.type.vx) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_vx[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(vx[iy*n1*n2+ix2*n1+iz]*vx[iy*n1*n2+ix2*n1+iz]);
+ }
+ }
+ if (sna.type.vy) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_vy[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(vy[iy2*n1*n2+ix*n1+iz]*vx[iy2*n1*n2+ix*n1+iz]);
+ }
+ }
+ if (sna.type.vz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_vz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(vz[iy*n1*n2+ix*n1+iz+1]*vz[iy*n1*n2+ix*n1+iz+1]);
+ }
+ }
+ if (sna.type.p) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_p[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tzz[iy*n1*n2+ix*n1+iz]*tzz[iy*n1*n2+ix*n1+iz]);
+ }
+ }
+ if (sna.type.tzz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_tzz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tzz[iy*n1*n2+ix*n1+iz]*tzz[iy*n1*n2+ix*n1+iz]);
+ }
+ }
+ if (sna.type.tyy) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_tyy[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tyy[iy*n1*n2+ix*n1+iz]*tyy[iy*n1*n2+ix*n1+iz]);
+ }
+ }
+ if (sna.type.txx) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_txx[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(txx[iy*n1*n2+ix*n1+iz]*txx[iy*n1*n2+ix*n1+iz]);
+ }
+ }
+ if (sna.type.txz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_txz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(txz[iy*n1*n2+ix2*n1+iz+1]*txz[iy*n1*n2+ix2*n1+iz+1]);
+ }
+ }
+ if (sna.type.tyz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_tyz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tyz[iy2*n1*n2+ix*n1+iz+1]*tyz[iy2*n1*n2+ix*n1+iz+1]);
+ }
+ }
+ if (sna.type.txz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ beam_txy[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(txy[iy2*n1*n2+ix2*n1+iz]*txy[iy2*n1*n2+ix2*n1+iz]);
+ }
+ }
+ /* calculate divergence of velocity field */
+ if (sna.type.pp) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ iz2 = iz+1;
+ p = sdx*((vx[iy*n1*n2+ix2*n1+iz]-vx[iy*n1*n2+ix*n1+iz])+
+ (vy[iy2*n1*n2+ix*n1+iz]-vy[iy*n1*n2+ix*n1+iz])+
+ (vz[iy*n1*n2+ix*n1+iz2]-vz[iy*n1*n2+ix*n1+iz]));
+ beam_pp[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(p*p);
+ }
+ }
+ /* calculate rotation of velocity field */
+ if (sna.type.ss) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ iz2 = iz+1;
+ s = sdx*((vx[iy2*n1*n2+ix2*n1+iz2]-vx[iy*n1*n2+ix2*n1+iz])-
+ (vy[iy2*n1*n2+ix2*n1+iz2]-vy[iy2*n1*n2+ix*n1+iz])-
+ (vz[iy2*n1*n2+ix2*n1+iz2]-vz[iy*n1*n2+ix*n1+iz2]));
+ beam_ss[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(s*s);
+ }
+ }
+ }
+ }
+ return 0;
+}
+
+
+long writeBeams3D(modPar mod, snaPar sna, long ixsrc, long iysrc, long izsrc, long ishot, long fileno,
+ float *beam_vx, float *beam_vy, float *beam_vz, float *beam_txx, float *beam_tyy, float *beam_tzz, float *beam_txz, float *beam_tyz, float *beam_txy,
+ float *beam_p, float *beam_pp, float *beam_ss, long verbose)
+{
+ FILE *fpvx, *fpvy, *fpvz, *fptxx, *fptyy, *fptzz, *fptxz, *fptyz, *fptxy, *fpp, *fppp, *fpss;
+ long append;
+ long ix, iy;
+ char number[16], filename[1024];
+ segy hdr;
+
+ if (sna.beam==0) return 0;
+ /* all beam snapshots are written to the same output file(s) */
+ if (ishot) append=1;
+ else append=0;
+
+ strcpy(filename, sna.file_beam);
+ if (fileno) {
+ sprintf(number,"_%03d",fileno);
+ name_ext(filename, number);
+ }
+ if (verbose>2) vmess("Writing beam data to file %s", filename);
+
+
+ if (sna.type.vx) fpvx = fileOpen(filename, "_bvx", (int)append);
+ if (sna.type.vy) fpvy = fileOpen(filename, "_bvy", (int)append);
+ if (sna.type.vz) fpvz = fileOpen(filename, "_bvz", (int)append);
+ if (sna.type.p) fpp = fileOpen(filename, "_bp", (int)append);
+ if (sna.type.txx) fptxx = fileOpen(filename, "_btxx", (int)append);
+ if (sna.type.tyy) fptyy = fileOpen(filename, "_btyy", (int)append);
+ if (sna.type.tzz) fptzz = fileOpen(filename, "_btzz", (int)append);
+ if (sna.type.txz) fptxz = fileOpen(filename, "_btxz", (int)append);
+ if (sna.type.tyz) fptyz = fileOpen(filename, "_btyz", (int)append);
+ if (sna.type.txy) fptxy = fileOpen(filename, "_btxy", (int)append);
+ if (sna.type.pp) fppp = fileOpen(filename, "_bpp", (int)append);
+ if (sna.type.ss) fpss = fileOpen(filename, "_bss", (int)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.sy = 1000*(mod.y0+iysrc*mod.dy);
+ hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
+ hdr.fldr = ishot+1;
+ hdr.trid = 1;
+ hdr.ns = sna.nz;
+ hdr.trwf = sna.nx*sna.ny;
+ hdr.ntr = sna.nx*sna.ny;
+ 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;
+
+ for (iy=0; iy<sna.ny; iy++) {
+ for (ix=0; ix<sna.nx; ix++) {
+ hdr.tracf = iy*sna.nx+ix+1;
+ hdr.tracl = iy*sna.nx+ix+1;
+ hdr.gx = 1000*(mod.x0+(sna.x1+ix)*mod.dx);
+ hdr.gy = 1000*(mod.y0+(sna.y1+ix)*mod.dy);
+
+ if (sna.type.vx) {
+ traceWrite( &hdr, &beam_vx[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpvx) ;
+ }
+ if (sna.type.vy) {
+ traceWrite( &hdr, &beam_vy[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpvy) ;
+ }
+ if (sna.type.vz) {
+ traceWrite( &hdr, &beam_vz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpvz) ;
+ }
+ if (sna.type.p) {
+ traceWrite( &hdr, &beam_p[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpp) ;
+ }
+ if (sna.type.tzz) {
+ traceWrite( &hdr, &beam_tzz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptzz) ;
+ }
+ if (sna.type.tyy) {
+ traceWrite( &hdr, &beam_tyy[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptyy) ;
+ }
+ if (sna.type.txx) {
+ traceWrite( &hdr, &beam_txx[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptxx) ;
+ }
+ if (sna.type.txz) {
+ traceWrite( &hdr, &beam_txz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptxz) ;
+ }
+ if (sna.type.txy) {
+ traceWrite( &hdr, &beam_txy[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptxy) ;
+ }
+ if (sna.type.tyz) {
+ traceWrite( &hdr, &beam_tyz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptyz) ;
+ }
+ if (sna.type.pp) {
+ traceWrite( &hdr, &beam_pp[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fppp) ;
+ }
+ if (sna.type.ss) {
+ traceWrite( &hdr, &beam_ss[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpss) ;
+ }
+ }
+ }
+
+ if (sna.type.vx) fclose(fpvx);
+ if (sna.type.vy) fclose(fpvy);
+ if (sna.type.vz) fclose(fpvz);
+ if (sna.type.p) fclose(fpp);
+ if (sna.type.txx) fclose(fptxx);
+ if (sna.type.tyy) fclose(fptyy);
+ if (sna.type.tzz) fclose(fptzz);
+ if (sna.type.txz) fclose(fptxz);
+ if (sna.type.tyz) fclose(fptyz);
+ if (sna.type.txy) fclose(fptxy);
+ if (sna.type.pp) fclose(fppp);
+ if (sna.type.ss) fclose(fpss);
+
+ return 0;
+}
\ No newline at end of file
diff --git a/fdelmodc3D/getRecTimes3D.c b/fdelmodc3D/getRecTimes3D.c
index ea2c276fc9c12e8d0f89bbf36f49dc306ec9cb13..33f5578732303b8951f70b416b89e703fdd1ccbc 100644
--- a/fdelmodc3D/getRecTimes3D.c
+++ b/fdelmodc3D/getRecTimes3D.c
@@ -17,10 +17,10 @@
* The Netherlands
**/
-int getRecTimes3D(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *l2m, float *rox, float *roy, float *roz, float *rec_vx, float *rec_vy, float *rec_vz, float *rec_txx, float *rec_tyy, float *rec_tzz, float *rec_txz, float *rec_txy, float *rec_tyz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose)
+long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *l2m, float *rox, float *roy, float *roz, float *rec_vx, float *rec_vy, float *rec_vz, float *rec_txx, float *rec_tyy, float *rec_tzz, float *rec_txz, float *rec_txy, float *rec_tyz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, long verbose)
{
- int n1, n2, ibndx, ibndy, ibndz;
- int irec, ix, iy, iz, ix2, iy2, iz2, ix1, iy1, iz1;
+ long n1, n2, ibndx, ibndy, ibndz;
+ long irec, ix, iy, iz, ix2, iy2, iz2, ix1, iy1, iz1;
float dvx, dvy, dvz, rdz, rdy, rdx;
float C000, C100, C010, C001, C110, C101, C011, C111;
float *vz_t, c1, c2, lroz, field;
@@ -63,9 +63,9 @@ int getRecTimes3D(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float
/* interpolation to precise (not necessary on a grid point) position */
if ( rec.int_p==3 ) {
- iz = (int)floorf(rec.zr[irec]/mod.dz)+ibndz;
- iy = (int)floorf(rec.yr[irec]/mod.dy)+ibndy;
- ix = (int)floorf(rec.xr[irec]/mod.dx)+ibndx;
+ iz = (long)floorf(rec.zr[irec]/mod.dz)+ibndz;
+ iy = (long)floorf(rec.yr[irec]/mod.dy)+ibndy;
+ ix = (long)floorf(rec.xr[irec]/mod.dx)+ibndx;
rdz = (rec.zr[irec] - (iz-ibndz)*mod.dz)/mod.dz;
rdy = (rec.yr[irec] - (iy-ibndy)*mod.dy)/mod.dy;
rdx = (rec.xr[irec] - (ix-ibndx)*mod.dx)/mod.dx;
@@ -98,198 +98,384 @@ int getRecTimes3D(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float
C010 = tzz[iy*n1*n2+ix*n1+iz+1];
C001 = tzz[(iy+1)*n1*n2+ix*n1+iz];
C110 = tzz[iy*n1*n2+(ix+1)*n1+iz+1];
- C101 = tzz[(iy+1)*n1*n2+ix*n1+iz+1];
- C011 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C101 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = tzz[(iy+1)*n1*n2+ix*n1+iz+1];
C111 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
- rec_p[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdy)*(1.0-rdz) +
- C100*rdx*(1.0-rdy)*(1.0-rdz) +
- C010*(1.0-rdx)*(1.0-rdy)*rdz +
- C11*rdx*rdz;
+ rec_p[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
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;
+ C000 = txx[iy*n1*n2+ix*n1+iz];
+ C100 = txx[iy*n1*n2+(ix+1)*n1+iz];
+ C010 = txx[iy*n1*n2+ix*n1+iz+1];
+ C001 = txx[(iy+1)*n1*n2+ix*n1+iz];
+ C110 = txx[iy*n1*n2+(ix+1)*n1+iz+1];
+ C101 = txx[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = txx[(iy+1)*n1*n2+ix*n1+iz+1];
+ C111 = txx[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_txx[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
+ }
+ if (rec.type.tyy) {
+ C000 = tyy[iy*n1*n2+ix*n1+iz];
+ C100 = tyy[iy*n1*n2+(ix+1)*n1+iz];
+ C010 = tyy[iy*n1*n2+ix*n1+iz+1];
+ C001 = tyy[(iy+1)*n1*n2+ix*n1+iz];
+ C110 = tyy[iy*n1*n2+(ix+1)*n1+iz+1];
+ C101 = tyy[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = tyy[(iy+1)*n1*n2+ix*n1+iz+1];
+ C111 = tyy[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_tyy[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
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;
+ C000 = tzz[iy*n1*n2+ix*n1+iz];
+ C100 = tzz[iy*n1*n2+(ix+1)*n1+iz];
+ C010 = tzz[iy*n1*n2+ix*n1+iz+1];
+ C001 = tzz[(iy+1)*n1*n2+ix*n1+iz];
+ C110 = tzz[iy*n1*n2+(ix+1)*n1+iz+1];
+ C101 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = tzz[(iy+1)*n1*n2+ix*n1+iz+1];
+ C111 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_tzz[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
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;
+ C000 = txz[iy*n1*n2+ix*n1+iz];
+ C100 = txz[iy*n1*n2+(ix+1)*n1+iz];
+ C010 = txz[iy*n1*n2+ix*n1+iz+1];
+ C001 = txz[(iy+1)*n1*n2+ix*n1+iz];
+ C110 = txz[iy*n1*n2+(ix+1)*n1+iz+1];
+ C101 = txz[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = txz[(iy+1)*n1*n2+ix*n1+iz+1];
+ C111 = txz[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_txz[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
+ }
+ if (rec.type.txy) {
+ C000 = txy[iy*n1*n2+ix*n1+iz];
+ C100 = txy[iy*n1*n2+(ix+1)*n1+iz];
+ C010 = txy[iy*n1*n2+ix*n1+iz+1];
+ C001 = txy[(iy+1)*n1*n2+ix*n1+iz];
+ C110 = txy[iy*n1*n2+(ix+1)*n1+iz+1];
+ C101 = txy[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = txy[(iy+1)*n1*n2+ix*n1+iz+1];
+ C111 = txy[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_txy[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
+ }
+ if (rec.type.tyz) {
+ C000 = tyz[iy*n1*n2+ix*n1+iz];
+ C100 = tyz[iy*n1*n2+(ix+1)*n1+iz];
+ C010 = tyz[iy*n1*n2+ix*n1+iz+1];
+ C001 = tyz[(iy+1)*n1*n2+ix*n1+iz];
+ C110 = tyz[iy*n1*n2+(ix+1)*n1+iz+1];
+ C101 = tyz[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = tyz[(iy+1)*n1*n2+ix*n1+iz+1];
+ C111 = tyz[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_tyz[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
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;
+ C000 = (vx[iy*n1*n2+ix2*n1+iz]-vx[iy*n1*n2+ix*n1+iz] +
+ vy[iy2*n1*n2+ix*n1+iz]-vy[iy*n1*n2+ix*n1+iz] +
+ vz[iy*n1*n2+ix*n1+iz2]-vz[iy*n1*n2+ix*n1+iz])/mod.dx;
+ C100 = (vx[iy*n1*n2+(ix2+1)*n1+iz]-vx[iy*n1*n2+(ix+1)*n1+iz] +
+ vy[iy2*n1*n2+(ix+1)*n1+iz]-vy[iy*n1*n2+(ix+1)*n1+iz] +
+ vz[iy*n1*n2+(ix+1)*n1+iz2]-vz[iy*n1*n2+(ix+1)*n1+iz])/mod.dx;
+ C010 = (vx[iy*n1*n2+ix2*n1+iz+1]-vx[iy*n1*n2+ix*n1+iz+1] +
+ vy[iy2*n1*n2+ix*n1+iz+1]-vy[iy*n1*n2+ix*n1+iz+1] +
+ vz[iy*n1*n2+ix*n1+iz2+1]-vz[iy*n1*n2+ix*n1+iz+1])/mod.dx;
+ C001 = (vx[(iy+1)*n1*n2+ix2*n1+iz]-vx[(iy+1)*n1*n2+ix*n1+iz] +
+ vy[(iy2+1)*n1*n2+ix*n1+iz]-vy[(iy+1)*n1*n2+ix*n1+iz] +
+ vz[(iy+1)*n1*n2+ix*n1+iz2]-vz[(iy+1)*n1*n2+ix*n1+iz])/mod.dx;
+ C110 = (vx[iy*n1*n2+(ix2+1)*n1+iz+1]-vx[iy*n1*n2+(ix+1)*n1+iz+1] +
+ vy[iy2*n1*n2+(ix+1)*n1+iz+1]-vy[iy*n1*n2+(ix+1)*n1+iz+1] +
+ vz[iy*n1*n2+(ix+1)*n1+iz2+1]-vz[iy*n1*n2+(ix+1)*n1+iz+1])/mod.dx;
+ C101 = (vx[(iy+1)*n1*n2+(ix2+1)*n1+iz]-vx[(iy+1)*n1*n2+(ix+1)*n1+iz] +
+ vy[(iy2+1)*n1*n2+(ix+1)*n1+iz]-vy[(iy+1)*n1*n2+(ix+1)*n1+iz] +
+ vz[(iy+1)*n1*n2+(ix+1)*n1+iz2]-vz[(iy+1)*n1*n2+(ix+1)*n1+iz])/mod.dx;
+ C011 = (vx[(iy+1)*n1*n2+ix2*n1+iz+1]-vx[(iy+1)*n1*n2+ix*n1+iz+1] +
+ vy[(iy2+1)*n1*n2+ix*n1+iz+1]-vy[(iy+1)*n1*n2+ix*n1+iz+1] +
+ vz[(iy+1)*n1*n2+ix*n1+iz2+1]-vz[(iy+1)*n1*n2+ix*n1+iz+1])/mod.dx;
+ C111 = (vx[(iy+1)*n1*n2+(ix2+1)*n1+iz+1]-vx[(iy+1)*n1*n2+(ix+1)*n1+iz+1] +
+ vy[(iy2+1)*n1*n2+(ix+1)*n1+iz+1]-vy[(iy+1)*n1*n2+(ix+1)*n1+iz+1] +
+ vz[(iy+1)*n1*n2+(ix+1)*n1+iz2+1]-vz[(iy+1)*n1*n2+(ix+1)*n1+iz+1])/mod.dx;
+ rec_pp[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
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;
+ C000 = (vx[iy2*n1*n2+ix2*n1+iz2]-vx[iy*n1*n2+ix2*n1+iz] -
+ (vy[iy2*n1*n2+ix2*n1+iz2]-vy[iy2*n1*n2+ix*n1+iz]) -
+ (vz[iy2*n1*n2+ix2*n1+iz2]-vz[iy*n1*n2+ix*n1+iz2]))/mod.dx;
+ C100 = (vx[iy2*n1*n2+(ix2+1)*n1+iz2]-vx[iy*n1*n2+(ix2+1)*n1+iz] -
+ (vy[iy2*n1*n2+(ix2+1)*n1+iz2]-vy[iy2*n1*n2+(ix+1)*n1+iz]) -
+ (vz[iy2*n1*n2+(ix2+1)*n1+iz2]-vz[iy*n1*n2+(ix+1)*n1+iz2]))/mod.dx;
+ C010 = (vx[iy2*n1*n2+ix2*n1+iz2+1]-vx[iy*n1*n2+ix2*n1+iz+1] -
+ (vy[iy2*n1*n2+ix2*n1+iz2+1]-vy[iy2*n1*n2+ix*n1+iz+1]) -
+ (vz[iy2*n1*n2+ix2*n1+iz2+1]-vz[iy*n1*n2+ix*n1+iz2+1]))/mod.dx;
+ C001 = (vx[(iy2+1)*n1*n2+ix2*n1+iz2]-vx[(iy+1)*n1*n2+ix2*n1+iz] -
+ (vy[(iy2+1)*n1*n2+ix2*n1+iz2]-vy[(iy2+1)*n1*n2+ix*n1+iz]) -
+ (vz[(iy2+1)*n1*n2+ix2*n1+iz2]-vz[(iy+1)*n1*n2+ix*n1+iz2]))/mod.dx;
+ C110 = (vx[iy2*n1*n2+(ix2+1)*n1+iz2+1]-vx[iy*n1*n2+(ix2+1)*n1+iz+1] -
+ (vy[iy2*n1*n2+(ix2+1)*n1+iz2+1]-vy[iy2*n1*n2+(ix+1)*n1+iz+1]) -
+ (vz[iy2*n1*n2+(ix2+1)*n1+iz2+1]-vz[iy*n1*n2+(ix+1)*n1+iz2+1]))/mod.dx;
+ C101 = (vx[(iy2+1)*n1*n2+(ix2+1)*n1+iz2]-vx[(iy+1)*n1*n2+(ix2+1)*n1+iz] -
+ (vy[(iy2+1)*n1*n2+(ix2+1)*n1+iz2]-vy[(iy2+1)*n1*n2+(ix+1)*n1+iz]) -
+ (vz[(iy2+1)*n1*n2+(ix2+1)*n1+iz2]-vz[(iy+1)*n1*n2+(ix+1)*n1+iz2]))/mod.dx;
+ C011 = (vx[(iy2+1)*n1*n2+ix2*n1+iz2+1]-vx[(iy+1)*n1*n2+ix2*n1+iz+1] -
+ (vy[(iy2+1)*n1*n2+ix2*n1+iz2+1]-vy[(iy2+1)*n1*n2+ix*n1+iz+1]) -
+ (vz[(iy2+1)*n1*n2+ix2*n1+iz2+1]-vz[(iy+1)*n1*n2+ix*n1+iz2+1]))/mod.dx;
+ C111 = (vx[(iy2+1)*n1*n2+(ix2+1)*n1+iz2+1]-vx[(iy+1)*n1*n2+(ix2+1)*n1+iz+1] -
+ (vy[(iy2+1)*n1*n2+(ix2+1)*n1+iz2+1]-vy[(iy2+1)*n1*n2+(ix+1)*n1+iz+1]) -
+ (vz[(iy2+1)*n1*n2+(ix2+1)*n1+iz2+1]-vz[(iy+1)*n1*n2+(ix+1)*n1+iz2+1]))/mod.dx;
+ rec_ss[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
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;
+ C000 = vz[iy*n1*n2+ix*n1+iz2];
+ C100 = vz[iy*n1*n2+(ix+1)*n1+iz2];
+ C010 = vz[iy*n1*n2+ix*n1+iz2+1];
+ C001 = vz[(iy+1)*n1*n2+ix*n1+iz2];
+ C110 = vz[iy*n1*n2+(ix+1)*n1+iz2+1];
+ C101 = vz[(iy+1)*n1*n2+(ix+1)*n1+iz2];
+ C011 = vz[(iy+1)*n1*n2+ix*n1+iz2+1];
+ C111 = vz[(iy+1)*n1*n2+(ix+1)*n1+iz2+1];
+ rec_vz[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
+ }
+ if (rec.type.vy) {
+ C000 = vy[iy2*n1*n2+ix*n1+iz];
+ C100 = vy[iy2*n1*n2+(ix+1)*n1+iz];
+ C010 = vy[iy2*n1*n2+ix*n1+iz+1];
+ C001 = vy[(iy2+1)*n1*n2+ix*n1+iz];
+ C110 = vy[iy2*n1*n2+(ix+1)*n1+iz+1];
+ C101 = vy[(iy2+1)*n1*n2+(ix+1)*n1+iz];
+ C011 = vy[(iy2+1)*n1*n2+ix*n1+iz+1];
+ C111 = vy[(iy2+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_vy[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
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;
+ C000 = vy[iy*n1*n2+ix2*n1+iz];
+ C100 = vy[iy*n1*n2+(ix2+1)*n1+iz];
+ C010 = vy[iy*n1*n2+ix2*n1+iz+1];
+ C001 = vy[(iy+1)*n1*n2+ix2*n1+iz];
+ C110 = vy[iy*n1*n2+(ix2+1)*n1+iz+1];
+ C101 = vy[(iy+1)*n1*n2+(ix2+1)*n1+iz];
+ C011 = vy[(iy+1)*n1*n2+ix2*n1+iz+1];
+ C111 = vy[(iy+1)*n1*n2+(ix2+1)*n1+iz+1];
+ rec_vx[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdz)*(1.0-rdy) +
+ C100*rdx*(1.0-rdz)*(1.0-rdy) +
+ C010*(1.0-rdx)*rdz*(1.0-rdy) +
+ C001*(1.0-rdx)*(1.0-rdz)*rdy +
+ C110*rdx*rdz*(1.0-rdy) +
+ C101*rdx*(1.0-rdz)*rdy +
+ C011*(1.0-rdx)*rdz*rdy +
+ C111*rdx*rdz*rdy;
}
}
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);
+ vmess("Receiver %li read at gridpoint ix=%li iy=%li iz=%li",irec, ix, iy, 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);
+ dvx = c1*(vx[iy*n1*n2+(ix+1)*n1+iz] - vx[iy*n1*n2+ix*n1+iz]) +
+ c2*(vx[iy*n1*n2+(ix+2)*n1+iz] - vx[iy*n1*n2+(ix-1)*n1+iz]);
+ dvy = c1*(vy[(iy+1)*n1*n2+ix*n1+iz] - vy[iy*n1*n2+ix*n1+iz]) +
+ c2*(vy[(iy+2)*n1*n2+ix*n1+iz] - vy[(iy-1)*n1*n2+ix*n1+iz]);
+ dvz = c1*(vz[iy*n1*n2+ix*n1+iz+1] - vz[iy*n1*n2+ix*n1+iz]) +
+ c2*(vz[iy*n1*n2+ix*n1+iz+2] - vz[iy*n1*n2+ix*n1+iz-1]);
+ field = tzz[iy*n1*n2+ix*n1+iz] + (1.0/2.0)*l2m[iy*n1*n2+ix*n1+iz]*(dvx+dvy+dvz);
+ dvx = c1*(vx[iy*n1*n2+(ix+1)*n1+iz1] - vx[iy*n1*n2+ix*n1+iz1]) +
+ c2*(vx[iy*n1*n2+(ix+2)*n1+iz1] - vx[iy*n1*n2+(ix-1)*n1+iz1]);
+ dvy = c1*(vy[(iy+1)*n1*n2+ix*n1+iz1] - vy[iy*n1*n2+ix*n1+iz1]) +
+ c2*(vy[(iy+2)*n1*n2+ix*n1+iz1] - vy[(iy-1)*n1*n2+ix*n1+iz1]);
+ dvz = c1*(vz[iy*n1*n2+ix*n1+iz1+1] - vz[iy*n1*n2+ix*n1+iz1]) +
+ c2*(vz[iy*n1*n2+ix*n1+iz1+2] - vz[iy*n1*n2+ix*n1+iz1-1]);
+ field += tzz[iy*n1*n2+ix*n1+iz1] + (1.0/2.0)*l2m[iy*n1*n2+ix*n1+iz1]*(dvx+dvy+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]);
+ rec_p[irec*rec.nt+isam] = 0.5*(tzz[iy*n1*n2+ix*n1+iz1]+tzz[iy*n1*n2+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);
+ dvx = c1*(vx[iy*n1*n2+(ix+1)*n1+iz] - vx[iy*n1*n2+ix*n1+iz]) +
+ c2*(vx[iy*n1*n2+(ix+2)*n1+iz] - vx[iy*n1*n2+(ix-1)*n1+iz]);
+ dvy = c1*(vy[(iy+1)*n1*n2+ix*n1+iz] - vy[iy*n1*n2+ix*n1+iz]) +
+ c2*(vy[(iy+2)*n1*n2+ix*n1+iz] - vy[(iy-1)*n1*n2+ix*n1+iz]);
+ dvz = c1*(vz[iy*n1*n2+ix*n1+iz+1] - vz[iy*n1*n2+ix*n1+iz]) +
+ c2*(vz[iy*n1*n2+ix*n1+iz+2] - vz[iy*n1*n2+ix*n1+iz-1]);
+ field = tzz[iy*n1*n2+ix*n1+iz] + (1.0/2.0)*l2m[iy*n1*n2+ix*n1+iz]*(dvx+dvy+dvz);
+ dvx = c1*(vx[iy*n1*n2+(ix1+1)*n1+iz] - vx[iy*n1*n2+ix1*n1+iz]) +
+ c2*(vx[iy*n1*n2+(ix1+2)*n1+iz] - vx[iy*n1*n2+(ix1-1)*n1+iz]);
+ dvy = c1*(vy[(iy+1)*n1*n2+ix1*n1+iz] - vy[iy*n1*n2+ix1*n1+iz]) +
+ c2*(vy[(iy+2)*n1*n2+ix1*n1+iz] - vy[(iy-1)*n1*n2+ix1*n1+iz]);
+ dvz = c1*(vz[iy*n1*n2+ix1*n1+iz+1] - vz[iy*n1*n2+ix1*n1+iz]) +
+ c2*(vz[iy*n1*n2+ix1*n1+iz+2] - vz[iy*n1*n2+ix1*n1+iz-1]);
+ field += tzz[iy*n1*n2+ix1*n1+iz] + (1.0/2.0)*l2m[iy*n1*n2+ix1*n1+iz]*(dvx+dvy+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]);
+ rec_p[irec*rec.nt+isam] = 0.5*(tzz[iy*n1*n2+ix1*n1+iz]+tzz[iy*n1*n2+ix*n1+iz]);
}
}
else {
- rec_p[irec*rec.nt+isam] = tzz[ix*n1+iz];
+ rec_p[irec*rec.nt+isam] = tzz[iy*n1*n2+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.txx) rec_txx[irec*rec.nt+isam] = txx[iy*n1*n2+ix*n1+iz];
+ if (rec.type.tzz) rec_tzz[irec*rec.nt+isam] = tzz[iy*n1*n2+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]);
+ rec_txz[irec*rec.nt+isam] = 0.125*(
+ txz[iy*n1*n2+ix*n1+iz] + txz[iy*n1*n2+ix*n1+iz2]+
+ txz[iy2*n1*n2+ix*n1+iz] + txz[iy2*n1*n2+ix*n1+iz2]+
+ txz[iy*n1*n2+ix2*n1+iz] + txz[iy*n1*n2+ix2*n1+iz2]+
+ txz[iy2*n1*n2+ix2*n1+iz] + txz[iy2*n1*n2+ix2*n1+iz2]);
}
else {
- rec_txz[irec*rec.nt+isam] = txz[ix2*n1+iz2];
+ rec_txz[irec*rec.nt+isam] = txz[iy2*n1*n2+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;
+ rec_pp[irec*rec.nt+isam] = (vx[iy*n1*n2+ix2*n1+iz]-vx[iy*n1*n2+ix*n1+iz] +
+ vy[iy2*n1*n2+ix*n1+iz]-vy[iy*n1*n2+ix*n1+iz] +
+ vz[iy*n1*n2+ix*n1+iz2]-vz[iy*n1*n2+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;
+ rec_ss[irec*rec.nt+isam] = (vx[iy2*n1*n2+ix2*n1+iz2]-vx[iy*n1*n2+ix2*n1+iz] -
+ (vy[iy2*n1*n2+ix2*n1+iz2]-vy[iy2*n1*n2+ix*n1+iz]) -
+ (vz[iy2*n1*n2+ix2*n1+iz2]-vz[iy*n1*n2+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]);
+ rec_vz[irec*rec.nt+isam] = 0.125*(
+ vz[iy*n1*n2+ix*n1+iz2]+vz[iy*n1*n2+ix1*n1+iz2]+
+ vz[iy*n1*n2+ix*n1+iz] +vz[iy*n1*n2+ix1*n1+iz]+
+ vz[iy1*n1*n2+ix*n1+iz2]+vz[iy1*n1*n2+ix1*n1+iz2]+
+ vz[iy1*n1*n2+ix*n1+iz] +vz[iy1*n1*n2+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]));
+ field = vz[iy*n1*n2+ix*n1+iz] - 0.5*roz[iy*n1*n2+ix*n1+iz]*(
+ c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) +
+ c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+ field += vz[iy*n1*n2+ix*n1+iz2] - 0.5*roz[iy*n1*n2+ix*n1+iz2]*(
+ c1*(tzz[iy*n1*n2+ix*n1+iz2] - tzz[iy*n1*n2+ix*n1+iz2-1]) +
+ c2*(tzz[iy*n1*n2+ix*n1+iz2+1] - tzz[iy*n1*n2+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]);
+ rec_vz[irec*rec.nt+isam] = 0.5*(vz[iy*n1*n2+ix*n1+iz2]+vz[iy*n1*n2+ix*n1+iz]);
}
}
else {
- rec_vz[irec*rec.nt+isam] = vz[ix*n1+iz2];
+ rec_vz[irec*rec.nt+isam] = vz[iy*n1*n2+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]);
+ //fprintf(stderr,"isam=%li vz[%li]=%e vz[%li]=%e vz[%li]=%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]);
+ rec_vx[irec*rec.nt+isam] = 0.125*(
+ vx[iy*n1*n2+ix2*n1+iz]+vx[iy*n1*n2+ix2*n1+iz1]+
+ vx[iy*n1*n2+ix*n1+iz]+vx[iy*n1*n2+ix*n1+iz1]+
+ vx[iy2*n1*n2+ix2*n1+iz]+vx[iy2*n1*n2+ix2*n1+iz1]+
+ vx[iy2*n1*n2+ix*n1+iz]+vx[iy2*n1*n2+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[iy*n1*n2+ix*n1+iz] - 0.5*rox[iy*n1*n2+ix*n1+iz]*(
+ c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) +
+ c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(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]));
+ c1*(tzz[iy*n1*n2+ix2*n1+iz] - tzz[iy*n1*n2+(ix2-1)*n1+iz]) +
+ c2*(tzz[iy*n1*n2+(ix2+1)*n1+iz] - tzz[iy*n1*n2+(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]);
+ rec_vx[irec*rec.nt+isam] = 0.5*(vx[iy*n1*n2+ix2*n1+iz]+vx[iy*n1*n2+ix*n1+iz]);
}
}
else {
- rec_vx[irec*rec.nt+isam] = vx[ix2*n1+iz];
+ rec_vx[irec*rec.nt+isam] = vx[iy*n1*n2+ix2*n1+iz];
}
}
}
@@ -300,22 +486,26 @@ int getRecTimes3D(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float
if (rec.type.ud) {
iz = rec.z[0]+ibndz;
iz2 = iz+1;
- vz_t = (float *)calloc(2*mod.nax,sizeof(float));
+ vz_t = (float *)calloc(2*mod.nax*mod.nay,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];
+ for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
+ for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
+ vz_t[iy*mod.nax+ix] = vz[iy*n1*n2+ix*n1+iz] - lroz*(
+ c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) +
+ c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+ vz_t[iy*mod.nax+mod.nay*mod.nax+ix] = vz[iy*n1*n2+ix*n1+iz2] - lroz*(
+ c1*(tzz[iy*n1*n2+ix*n1+iz2] - tzz[iy*n1*n2+ix*n1+iz2-1]) +
+ c2*(tzz[iy*n1*n2+ix*n1+iz2+1] - tzz[iy*n1*n2+ix*n1+iz2-2]));
+ }
+ }
+ for (iy=0; iy<mod.nay; iy++) {
+ 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[iy*mod.nax*rec.nt+ix*rec.nt+isam] = 0.25*(vz[iy*n1*n2+ix*n1+iz2]+vz[iy*n1*n2+ix*n1+iz]+vz_t[iy*mod.nax+mod.nay*mod.nax+ix]+vz_t[iy*mod.nax+ix]);
+ rec_udp[iy*mod.nax*rec.nt+ix*rec.nt+isam] = tzz[iy*n1*n2+ix*n1+iz];
+ }
}
free(vz_t);
}
diff --git a/fdelmodc3D/writeRec3D.c b/fdelmodc3D/writeRec3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..4dd9cb4ff5debb11049f283506aac02d57946ea4
--- /dev/null
+++ b/fdelmodc3D/writeRec3D.c
@@ -0,0 +1,255 @@
+#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 "segy.h"
+#include "fdelmodc3D.h"
+#include <genfft.h>
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+typedef struct _dcomplexStruct { /* complex number */
+ double r,i;
+} dcomplex;
+#endif/* complex */
+
+/**
+* Writes the receiver array(s) 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) ;
+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 vznorm, int verbose);
+
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+long writeRec3D(recPar rec, modPar mod, bndPar bnd, wavPar wav, long ixsrc, long iysrc, long izsrc, long nsam, long ishot, long fileno,
+ float *rec_vx, float *rec_vy, float *rec_vz, float *rec_txx, float *rec_tyy, float *rec_tzz, float *rec_txz, float *rec_tyz, float *rec_txy,
+ float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, long verbose)
+{
+ FILE *fpvx, *fpvy, *fpvz, *fptxx, *fptyy, *fptzz, *fptxz, *fptyz, *fptxy, *fpp, *fppp, *fpss, *fpup, *fpdown;
+ float *rec_up, *rec_down, *trace, *rec_vze, *rec_pe;
+ float dx, dy, dt, cp, rho, fmin, fmax;
+ complex *crec_vz, *crec_p, *crec_up, *crec_dw;
+ long irec, ntfft, nfreq, nkx, nky, xorig, yorig, ix, iy, iz, it, ibndx, ibndy;
+ long append, vznorm, sx, sy;
+ double ddt;
+ char number[16], filename[1024];
+ segy hdr;
+
+ if (!rec.n) return 0;
+ if (ishot) append=1;
+ else append=0;
+
+ /* if the total number of samples exceeds rec_ntsam then a new (numbered) file is opened */
+ /* fileno has a non-zero value (from fdelmodc.c) if the number of samples exceeds rec_ntsam. */
+ strcpy(filename, rec.file_rcv);
+ if (fileno) {
+ sprintf(number,"_%03d",fileno);
+ name_ext(filename, number);
+ }
+#ifdef MPI
+ sx = (long)mod.x0+ixsrc*mod.dx;
+ sprintf(number,"_%06d",sx);
+ name_ext(filename, number);
+#endif
+
+ if (verbose>2) vmess("Writing receiver data to file %s", filename);
+ if (nsam != rec.nt && verbose) vmess("Number of samples written to last file = %li",nsam);
+
+ memset(&hdr,0,TRCBYTES);
+ ddt = (double)mod.dt;/* to avoid rounding in 32 bit precision */
+ dt = (float)ddt*rec.skipdt;
+ dx = (rec.x[1]-rec.x[0])*mod.dx;
+ dy = (rec.y[1]-rec.y[0])*mod.dy;
+ hdr.dt = (unsigned short)lround((((double)1.0e6*ddt*rec.skipdt)));
+ hdr.scalco = -1000;
+ hdr.scalel = -1000;
+ hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
+ hdr.sy = 1000*(mod.y0+ixsrc*mod.dy);
+ hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
+ hdr.selev = (int)(-1000.0*(mod.z0+izsrc*mod.dz));
+ hdr.fldr = ishot+1;
+ hdr.trid = 1;
+ hdr.ns = nsam;
+ hdr.trwf = rec.n;
+ hdr.ntr = rec.n;
+ if (mod.grid_dir) { /* reverse time modeling */
+ hdr.f1 = (-mod.nt+1)*mod.dt;
+ }
+ else {
+ hdr.f1 = 0.0;
+ }
+ hdr.d1 = mod.dt*rec.skipdt;
+ hdr.d2 = (rec.x[1]-rec.x[0])*mod.dx;
+ hdr.f2 = mod.x0+rec.x[0]*mod.dx;
+
+ if (rec.type.vx) fpvx = fileOpen(filename, "_rvx", (int)append);
+ if (rec.type.vy) fpvy = fileOpen(filename, "_rvy", (int)append);
+ if (rec.type.vz) fpvz = fileOpen(filename, "_rvz", (int)append);
+ if (rec.type.p) fpp = fileOpen(filename, "_rp", (int)append);
+ if (rec.type.txx) fptxx = fileOpen(filename, "_rtxx", (int)append);
+ if (rec.type.tyy) fptyy = fileOpen(filename, "_rtyy", (int)append);
+ if (rec.type.tzz) fptzz = fileOpen(filename, "_rtzz", (int)append);
+ if (rec.type.txz) fptxz = fileOpen(filename, "_rtxz", (int)append);
+ if (rec.type.tyz) fptyz = fileOpen(filename, "_rtyz", (int)append);
+ if (rec.type.txy) fptxy = fileOpen(filename, "_rtxy", (int)append);
+ if (rec.type.pp) fppp = fileOpen(filename, "_rpp", (int)append);
+ if (rec.type.ss) fpss = fileOpen(filename, "_rss", (int)append);
+
+ /* decomposed wavefield */
+ if (rec.type.ud && (mod.ischeme==1 || mod.ischeme==2) ) {
+ fpup = fileOpen(filename, "_ru", (int)append);
+ fpdown = fileOpen(filename, "_rd", (int)append);
+ ntfft = optncr(nsam);
+ nfreq = ntfft/2+1;
+ fmin = 0.0;
+ fmax = wav.fmax;
+ nkx = optncc(2*mod.nax);
+ nky = optncc(2*mod.nay);
+ ibndx = mod.ioPx;
+ ibndy = mod.ioPy;
+ if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+ if (bnd.fro==4 || bnd.fro==2) ibndy += 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*nky,sizeof(float));
+ rec_down= (float *)calloc(ntfft*nkx*nky,sizeof(float));
+ crec_vz = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+ crec_p = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+ crec_up = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+ crec_dw = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+
+ rec_vze = rec_up;
+ rec_pe = rec_down;
+ /* copy input data into extended arrays with padded zeroes */
+ for (iy=0; iy<mod.nay; iy++) {
+ for (ix=0; ix<mod.nax; ix++) {
+ memcpy(&rec_vze[iy*mod.nax*ntfft+ix*ntfft],&rec_udvz[iy*mod.nax*rec.nt+ix*rec.nt],nsam*sizeof(float));
+ memcpy(&rec_pe[iy*mod.nax*ntfft+ix*ntfft], &rec_udp[iy*mod.nax*rec.nt+ix*rec.nt], nsam*sizeof(float));
+ }
+ }
+
+ /* transform from t-x to kx-w */
+ xorig = ixsrc+ibndx;
+ xt2wkx(rec_vze, crec_vz, ntfft, nkx, ntfft, nkx, xorig);
+ xt2wkx(rec_pe, crec_p, ntfft, nkx, ntfft, nkx, xorig);
+
+ /* apply decomposition operators */
+ kxwdecomp(crec_p, crec_vz, crec_up, crec_dw,
+ 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);
+ wkx2xt(crec_dw, rec_down, ntfft, nkx, nkx, ntfft, xorig);
+
+ /* reduce array to rec.nt samples rec.n traces */
+ for (irec=0; irec<rec.n; irec++) {
+ ix = rec.x[irec]+ibndx;
+ iy = rec.y[irec]+ibndy;
+ for (it=0; it<rec.nt; it++) {
+ rec_up[irec*rec.nt+it] = rec_up[iy*nkx*ntfft+ix*ntfft+it];
+ rec_down[irec*rec.nt+it] = rec_down[iy*nkx*ntfft+ix*ntfft+it];
+ }
+ }
+ free(crec_vz);
+ free(crec_p);
+ free(crec_up);
+ free(crec_dw);
+ }
+ if (rec.type.ud && (mod.ischeme==3 || mod.ischeme==4) ) {
+ }
+
+ for (irec=0; irec<rec.n; irec++) {
+ hdr.tracf = irec+1;
+ hdr.tracl = ishot*rec.n+irec+1;
+ hdr.gx = 1000*(mod.x0+rec.x[irec]*mod.dx);
+ hdr.gy = 1000*(mod.y0+rec.y[irec]*mod.dy);
+ hdr.offset = (rec.x[irec]-ixsrc)*mod.dx;
+ hdr.gelev = (int)(-1000*(mod.z0+rec.z[irec]*mod.dz));
+
+ if (rec.type.vx) {
+ traceWrite( &hdr, &rec_vx[irec*rec.nt], (int)nsam, fpvx) ;
+ }
+ if (rec.type.vy) {
+ traceWrite( &hdr, &rec_vy[irec*rec.nt], (int)nsam, fpvy) ;
+ }
+ if (rec.type.vz) {
+ traceWrite( &hdr, &rec_vz[irec*rec.nt], (int)nsam, fpvz) ;
+ }
+ if (rec.type.p) {
+ traceWrite( &hdr, &rec_p[irec*rec.nt], (int)nsam, fpp) ;
+ }
+ if (rec.type.txx) {
+ traceWrite( &hdr, &rec_txx[irec*rec.nt], (int)nsam, fptxx) ;
+ }
+ if (rec.type.tyy) {
+ traceWrite( &hdr, &rec_tyy[irec*rec.nt], (int)nsam, fptyy) ;
+ }
+ if (rec.type.tzz) {
+ traceWrite( &hdr, &rec_tzz[irec*rec.nt], (int)nsam, fptzz) ;
+ }
+ if (rec.type.txz) {
+ traceWrite( &hdr, &rec_txz[irec*rec.nt], (int)nsam, fptxz) ;
+ }
+ if (rec.type.txy) {
+ traceWrite( &hdr, &rec_txy[irec*rec.nt], (int)nsam, fptxy) ;
+ }
+ if (rec.type.tyz) {
+ traceWrite( &hdr, &rec_tyz[irec*rec.nt], (int)nsam, fptyz) ;
+ }
+ if (rec.type.pp) {
+ traceWrite( &hdr, &rec_pp[irec*rec.nt], (int)nsam, fppp) ;
+ }
+ if (rec.type.ss) {
+ traceWrite( &hdr, &rec_ss[irec*rec.nt], (int)nsam, fpss) ;
+ }
+ if (rec.type.ud && mod.ischeme==1) {
+ traceWrite( &hdr, &rec_up[irec*rec.nt], (int)nsam, fpup) ;
+ traceWrite( &hdr, &rec_down[irec*rec.nt], (int)nsam, fpdown) ;
+ }
+ }
+
+ if (rec.type.vx) fclose(fpvx);
+ if (rec.type.vy) fclose(fpvy);
+ if (rec.type.vz) fclose(fpvz);
+ if (rec.type.p) fclose(fpp);
+ if (rec.type.txx) fclose(fptxx);
+ if (rec.type.tyy) fclose(fptyy);
+ if (rec.type.tzz) fclose(fptzz);
+ if (rec.type.txz) fclose(fptxz);
+ if (rec.type.txy) fclose(fptxy);
+ if (rec.type.tyz) fclose(fptyz);
+ if (rec.type.pp) fclose(fppp);
+ if (rec.type.ss) fclose(fpss);
+ if (rec.type.ud) {
+ fclose(fpup);
+ fclose(fpdown);
+ free(rec_up);
+ free(rec_down);
+ }
+
+ return 0;
+}
+
diff --git a/fdelmodc3D/writeSnapTimes3D.c b/fdelmodc3D/writeSnapTimes3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..56c81f5bec79781651efa9a0aaad1bf1e93f86f2
--- /dev/null
+++ b/fdelmodc3D/writeSnapTimes3D.c
@@ -0,0 +1,261 @@
+#define _FILE_OFFSET_BITS 64
+#define _LARGEFILE_SOURCE
+#define _LARGEFILE64_SOURCE
+
+#define ISODD(n) ((n) & 01)
+
+#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 "fdelmodc3D.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) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+long writeSnapTimes3D(modPar mod, snaPar sna, bndPar bnd, wavPar wav, long ixsrc, long iysrc, long izsrc, long itime, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *tyz, float *txy, long verbose)
+{
+ FILE *fpvx, *fpvy, *fpvz, *fptxx, *fptyy, *fptzz, *fptxz, *fptyz, *fptxy, *fpp, *fppp, *fpss;
+ long append, isnap;
+ static long first=1;
+ long n1, n2, ibndx, ibndy, ibndz, ixs, iys, izs, ize, i, j, l;
+ long ix, iy, iz, ix2, iy2;
+ float *snap, sdx, stime;
+ segy hdr;
+
+ if (sna.nsnap==0) return 0;
+
+ ibndx = mod.ioXx;
+ ibndy = mod.ioXy;
+ ibndz = mod.ioXz;
+ n1 = mod.naz;
+ n2 = mod.nax;
+ sdx = 1.0/mod.dx;
+
+ if (sna.withbnd) {
+ sna.nz=mod.naz;
+ sna.z1=0;
+ sna.z2=mod.naz-1;
+ sna.skipdz=1;
+
+ sna.ny=mod.nax;
+ sna.y1=0;
+ sna.y2=mod.nay-1;
+ sna.skipdy=1;
+
+ sna.nx=mod.nax;
+ sna.x1=0;
+ sna.x2=mod.nax-1;
+ sna.skipdx=1;
+ }
+
+ /* 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 (mod.grid_dir) stime = (-wav.nt+1+itime+1)*mod.dt; /* reverse time modeling */
+ else stime = itime*mod.dt;
+ if (verbose) vmess("Writing snapshot(%li) at time=%.4f", isnap+1, stime);
+
+ if (first) {
+ append=0;
+ first=0;
+ }
+ else {
+ append=1;
+ }
+
+ if (sna.type.vx) fpvx = fileOpen(sna.file_snap, "_svx", (int)append);
+ if (sna.type.vy) fpvy = fileOpen(sna.file_snap, "_svy", (int)append);
+ if (sna.type.vz) fpvz = fileOpen(sna.file_snap, "_svz", (int)append);
+ if (sna.type.p) fpp = fileOpen(sna.file_snap, "_sp", (int)append);
+ if (sna.type.txx) fptxx = fileOpen(sna.file_snap, "_stxx", (int)append);
+ if (sna.type.tyy) fptyy = fileOpen(sna.file_snap, "_styy", (int)append);
+ if (sna.type.tzz) fptzz = fileOpen(sna.file_snap, "_stzz", (int)append);
+ if (sna.type.txz) fptxz = fileOpen(sna.file_snap, "_stxz", (int)append);
+ if (sna.type.tyz) fptyz = fileOpen(sna.file_snap, "_styz", (int)append);
+ if (sna.type.txy) fptxy = fileOpen(sna.file_snap, "_stxy", (int)append);
+ if (sna.type.pp) fppp = fileOpen(sna.file_snap, "_spp", (int)append);
+ if (sna.type.ss) fpss = fileOpen(sna.file_snap, "_sss", (int)append);
+
+ memset(&hdr,0,TRCBYTES);
+ hdr.dt = 1000000*(sna.skipdt*mod.dt);
+ hdr.ungpow = (sna.delay*mod.dt);
+ hdr.scalco = -1000;
+ hdr.scalel = -1000;
+ hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
+ hdr.sy = 1000*(mod.y0+iysrc*mod.dy);
+ hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
+ hdr.fldr = isnap+1;
+ hdr.trid = 1;
+ hdr.ns = sna.nz;
+ hdr.trwf = sna.nx*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;
+ if (sna.withbnd) {
+ if ( !ISODD(bnd.top)) hdr.f1 = mod.z0 - bnd.ntap*mod.dz;
+ if ( !ISODD(bnd.lef)) hdr.f2 = mod.x0 - bnd.ntap*mod.dx;
+ //if ( !ISODD(bnd.rig)) ;
+ //if ( !ISODD(bnd.bot)) store=1;
+ }
+
+/***********************************************************************
+* 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 (iys=sna.y1, l=0; iys<=sna.x2; iys+=sna.skipdx, l++) {
+ for (ixs=sna.x1, i=0; ixs<=sna.x2; ixs+=sna.skipdx, i++) {
+ hdr.tracf = l*sna.nx+i+1;
+ hdr.tracl = isnap*sna.nx*sna.ny+l*sna.nx+i+1;
+ hdr.gx = 1000*(mod.x0+ixs*mod.dx);
+ hdr.gy = 1000*(mod.y0+ixs*mod.dy);
+ ix = ixs+ibndx;
+ ix2 = ix+1;
+ iy = iys+ibndy;
+ iy2 = iy+1;
+
+ izs = sna.z1+ibndz;
+ ize = sna.z2+ibndz;
+
+ if (sna.withbnd) {
+ izs = 0;
+ ize = sna.z2;
+ ix = ixs;
+ ix2 = ix;
+ iy = iys;
+ iy2 = iy;
+ if (sna.type.vz || sna.type.txz || sna.type.tyz) izs = -1;
+ if ( !ISODD(bnd.lef)) hdr.gx = 1000*(mod.x0 - bnd.ntap*mod.dx);
+ if ( !ISODD(bnd.fro)) hdr.gy = 1000*(mod.y0 - bnd.ntap*mod.dy);
+ }
+
+ if (sna.type.vx) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = vx[iy*n1*n2+ix2*n1+iz];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fpvx);
+ }
+ if (sna.type.vy) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = vy[iy2*n1*n2+ix*n1+iz];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fpvy);
+ }
+ if (sna.type.vz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = vz[iy*n1*n2+ix*n1+iz+1];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fpvz);
+ }
+ if (sna.type.p) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = tzz[iy*n1*n2+ix*n1+iz];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fpp);
+ }
+ if (sna.type.tzz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = tzz[iy*n1*n2+ix*n1+iz];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fptzz);
+ }
+ if (sna.type.tyy) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = tyy[iy*n1*n2+ix*n1+iz];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fptyy);
+ }
+ if (sna.type.txx) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = txx[iy*n1*n2+ix*n1+iz];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fptxx);
+ }
+ if (sna.type.txz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = txz[iy*n1*n2+ix2*n1+iz+1];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fptxz);
+ }
+ if (sna.type.txy) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = txy[iy2*n1*n2+ix2*n1+iz];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fptxy);
+ }
+ if (sna.type.tyz) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = tyz[iy2*n1*n2+ix*n1+iz+1];
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fptyz);
+ }
+ /* calculate divergence of velocity field */
+ if (sna.type.pp) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = sdx*((vx[iy*n1*n2+(ix+1)*n1+iz]-vx[iy*n1*n2+ix*n1+iz])+
+ (vy[(iy+1)*n1*n2+ix*n1+iz]-vy[iy*n1*n2+ix*n1+iz])+
+ (vz[iy*n1*n2+ix*n1+iz+1]-vz[iy*n1*n2+ix*n1+iz]));
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fppp);
+ }
+ /* calculate rotation of velocity field */
+ if (sna.type.ss) {
+ for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
+ snap[j] = sdx*((vx[iy*n1*n2+ix*n1+iz]-vx[(iy-1)*n1*n2+ix*n1+iz-1])-
+ (vy[iy*n1*n2+ix*n1+iz]-vy[iy*n1*n2+(ix-1)*n1+iz-1])-
+ (vz[iy*n1*n2+ix*n1+iz]-vz[(iy-1)*n1*n2+(ix-1)*n1+iz]));
+ }
+ traceWrite(&hdr, snap, (int)sna.nz, fpss);
+ }
+
+ }
+ }
+
+ if (sna.type.vx) fclose(fpvx);
+ if (sna.type.vy) fclose(fpvy);
+ if (sna.type.vz) fclose(fpvz);
+ if (sna.type.p) fclose(fpp);
+ if (sna.type.txx) fclose(fptxx);
+ if (sna.type.tyy) fclose(fptyy);
+ if (sna.type.tzz) fclose(fptzz);
+ if (sna.type.txz) fclose(fptxz);
+ if (sna.type.tyz) fclose(fptyz);
+ if (sna.type.txy) fclose(fptxy);
+ if (sna.type.pp) fclose(fppp);
+ if (sna.type.ss) fclose(fpss);
+
+ free(snap);
+ }
+
+ return 0;
+}
+