3D

fdelmodc3D/acoustic4_3D.c

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+long applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, long verbose);
+
+long storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, long verbose);
+
+long reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, long verbose);
+
+long boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
+
+long boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
+
+long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, float *p, float *rox, float *roy, float *roz, float *l2m, long verbose)
+{
+/*********************************************************************
+       COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID: 
+
+  The captial symbols T (=P) Txz,Vx,Vz represent the actual grid
+  The indices ix,iz are related to the T grid, so the capital T 
+  symbols represent the actual modelled grid.
+
+  one cel (iz,ix)
+       |
+       V                              extra column of vx,txz
+                                                      |
+    -------                                           V
+   | txz vz| txz vz  txz vz  txz vz  txz vz  txz vz txz
+   |       |      
+   | vx  t | vx  t   vx  t   vx  t   vx  t   vx  t  vx
+    -------
+     txz vz  txz vz  txz vz  txz vz  txz vz  txz vz  txz
+
+     vx  t   vx  T---Vx--T---Vx--T---Vx--T   vx  t   vx
+                 |   |   |   |   |   |   | 
+     txz vz  txz Vz--Txz-Vz--Txz-Vz  Txz-Vz  txz vz  txz
+                 |   |   |   |   |   |   |
+     vx  t   vx  T---Vx--T---Vx--T---Vx--T   vx  t   vx
+                 |   |   |   |   |   |   |
+     txz vz  txz Vz  Txz-Vz  Txz-Vz  Txz-Vz  txz vz  txz
+                 |   |   |   |   |   |   |
+     vx  t   vx  T---Vx--T---Vx--T---Vx--T   vx  t   vx
+                 |   |   |   |   |   |   |
+     txz vz  txz Vz  Txz-Vz  Txz-Vz  Txz-Vz  txz vz  txz
+                 |   |   |   |   |   |   |
+     vx  t   vx  T---Vx--T---Vx--T---Vx--T   vx  t   vx
+
+     txz vz  txz vz  txz vz  txz vz  txz vz  txz vz  txz
+
+     vx  t   vx  t   vx  t   vx  t   vx  t   vx  t  vx
+
+     txz vz  txz vz  txz vz  txz vz  txz vz  txz vz  txz  <--| 
+                                                             |
+                                         extra row of txz/vz |
+
+   AUTHOR:
+           Jan Thorbecke (janth@xs4all.nl)
+           The Netherlands 
+
+***********************************************************************/
+
+	float   c1, c2;
+	long    ix, iy, iz;
+	long    n1, n2;
+	long    ioXx, ioXy, ioXz, ioYx, ioYy, ioYz, ioZx, ioZy, ioZz, ioPx, ioPy, ioPz;
+
+	c1 = 9.0/8.0; 
+	c2 = -1.0/24.0;
+	n1  = mod.naz;
+    n2  = mod.nay;
+
+	/* calculate vx for all grid points except on the virtual boundary*/
+#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
+	for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
+#pragma ivdep
+        for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
+            for (iz=mod.ioXz; iz<mod.ieXz; iz++) {
+                vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*(
+                    c1*(p[iy*n2*n1+ix*n1+iz]     - p[iy*n2*n1+(ix-1)*n1+iz]) +
+                    c2*(p[iy*n2*n1+(ix+1)*n1+iz] - p[iy*n2*n1+(ix-2)*n1+iz]));
+            }
+        }
+    }
+
+	/* calculate vy for all grid points except on the virtual boundary*/
+#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
+	for (iy=mod.ioYy; iy<mod.ieYy; iy++) {
+#pragma ivdep
+        for (ix=mod.ioYx; ix<mod.ieYx; ix++) {
+            for (iz=mod.ioYz; iz<mod.ieYz; iz++) {
+                vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*(
+                    c1*(p[iy*n2*n1+ix*n1+iz]     - p[(iy-1)*n2*n1+ix*n1+iz]) +
+                    c2*(p[(iy+1)*n2*n1+ix*n1+iz] - p[(iy-2)*n2*n1+ix*n1+iz]));
+            }
+        }
+    }
+
+	/* calculate vz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1) 
+	for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
+#pragma ivdep
+        for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
+            for (iz=mod.ioZz; iz<mod.ieZz; iz++) {
+                vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*(
+                    c1*(p[iy*n2*n1+ix*n1+iz]   - p[iy*n2*n1+ix*n1+iz-1]) +
+                    c2*(p[iy*n2*n1+ix*n1+iz+1] - p[iy*n2*n1+ix*n1+iz-2]));
+            }
+        }
+    }
+        
+	/* boundary condition clears velocities on boundaries */
+    boundariesP3D(mod, bnd, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, NULL, NULL, itime, verbose);
+
+	/* Add force source */
+	if (src.type > 5) {
+         applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, src_nwav, verbose);
+	}
+
+	/* this is needed because the P fields are not using tapered boundaries (bnd....=4) */
+    if (bnd.top==2) mod.ioPz += bnd.npml;
+    if (bnd.bot==2) mod.iePz -= bnd.npml;
+    if (bnd.fro==2) mod.ioPy += bnd.npml;
+    if (bnd.bac==2) mod.iePy -= bnd.npml;
+    if (bnd.lef==2) mod.ioPx += bnd.npml;
+    if (bnd.rig==2) mod.iePx -= bnd.npml;
+
+	/* calculate p/tzz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1) 
+//#pragma omp for private (ix, iz) schedule(dynamic) 
+#pragma ivdep
+	for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+#pragma ivdep
+		for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+            for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+                p[iy*n1*n2+ix*n1+iz] -= l2m[iy*n1*n2+ix*n1+iz]*(
+                            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]) +
+                            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]) +
+                            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]));
+            }
+        }
+	}
+    if (bnd.top==2) mod.ioPz -= bnd.npml;
+    if (bnd.bot==2) mod.iePz += bnd.npml;
+    if (bnd.fro==2) mod.ioPy -= bnd.npml;
+    if (bnd.bac==2) mod.iePy += bnd.npml;
+    if (bnd.lef==2) mod.ioPx -= bnd.npml;
+    if (bnd.rig==2) mod.iePx += bnd.npml;
+
+	/* Add stress source */
+	if (src.type < 6) {
+        applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, src_nwav, verbose);
+	}
+    
+/* Free surface: calculate free surface conditions for stresses */
+
+	/* check if there are sources placed on the free surface */
+    storeSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, verbose);
+
+	/* Free surface: calculate free surface conditions for stresses */
+    boundariesV3D(mod, bnd, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, NULL, NULL, itime, verbose);
+
+	/* restore source positions on the edge */
+    reStoreSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, verbose);
+
+	return 0;
+}

fdelmodc3D/applySource3D.c

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+
+void vmess(char *fmt, ...);
+
+#define c1 (9.0/8.0)
+#define c2 (-1.0/24.0)
+
+/*********************************************************************
+ * 
+ * Add's the source amplitude(s) to the grid.
+ * 
+ * For the acoustic schemes, the source-type must not be txx tzz or txz.
+ *
+ *   AUTHOR:
+ *           Jan Thorbecke (janth@xs4all.nl)
+ *           The Netherlands 
+ *
+ **********************************************************************/
+
+long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose)
+{
+	long is0, ibndz, ibndy, ibndx;
+	long isrc, ix, iy, iz, n1, n2;
+	long id1, id2, id3;
+	float src_ampl, time, scl, dt, sdx;
+	float Mxx, Myy, Mzz, Mxz, Myz, Mxy, rake;
+	static long first=1;
+
+	if (src.type==6) {
+    	ibndz = mod.ioXz;
+    	ibndy = mod.ioXy;
+    	ibndx = mod.ioXx;
+	}
+	else if (src.type==7) {
+    	ibndz = mod.ioZz;
+    	ibndy = mod.ioZy;
+    	ibndx = mod.ioZx;
+	}
+	else if (src.type==2) {
+    	ibndz = mod.ioTz;
+    	ibndy = mod.ioTy;
+    	ibndx = mod.ioTx;
+    	if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+    	if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
+    	if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+	}
+	else {	
+    	ibndz = mod.ioPz;
+    	ibndy = mod.ioPy;
+    	ibndx = mod.ioPx;
+    	if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+    	if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
+    	if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+	}
+
+	n1   = mod.naz;
+    n2   = mod.nax;
+	dt   = mod.dt;
+	sdx  = 1.0/mod.dx;
+
+	/* special txz source activated? */
+
+	if ((bnd.top==1) && (src.type==2)) {
+		iz = izsrc + ibndz;
+		if (iz==ibndz) {
+            if (src.orient != 1) {
+				if (first) {
+					vmess("Only monopole Txz source allowed at surface. Reset to monopole");
+					first = 0;
+				}
+				src.orient=1;
+			}
+		}
+	}
+             
+/*
+* for plane wave sources the sources are placed 
+* around the central shot position 
+* the first source position has an offset in x of is0
+*
+* itime = 0 corresponds with time=0
+* itime = 1 corresponds with time=dt
+* src[0] (the first sample) corresponds with time = 0
+*/
+
+	is0 = -1*floor((src.n-1)/2);
+#pragma omp	for private (isrc, src_ampl, ix, iy, iz, time, id1, id2, id3, scl) 
+	for (isrc=0; isrc<src.n; isrc++) {
+		src_ampl=0.0;
+		/* calculate the source position */
+		if (src.random || src.multiwav) {
+			ix = src.x[isrc] + ibndx;
+			iy = src.y[isrc] + ibndy;
+			iz = src.z[isrc] + ibndz;
+		}
+		else { /* plane wave and point sources */
+            ix = ixsrc + ibndx + is0 + isrc;
+            iy = iysrc + ibndy + is0 + isrc;
+            iz = izsrc + ibndz;
+		}
+		time = itime*dt - src.tbeg[isrc];
+		id1 = floor(time/dt);
+		id2 = id1+1;
+        
+		/* delay not reached or no samples left in source wavelet? */
+		if ( (time < 0.0) || ( (itime*dt) >= src.tend[isrc]) ) continue;
+
+//		fprintf(stderr,"isrc=%li ix=%li iz=%li src.x=%li src.z=%li\n", isrc, ix, iz, src.x[isrc], src.z[isrc]);
+
+		if (!src.multiwav) { /* only one wavelet for all sources */
+			src_ampl = src_nwav[0][id1]*(id2-time/dt) + src_nwav[0][id2]*(time/dt-id1);
+		}
+		else { /* multi-wavelet sources */
+			src_ampl = src_nwav[isrc][id1]*(id2-time/dt) + src_nwav[isrc][id2]*(time/dt-id1);
+		}
+
+		if (src_ampl==0.0) continue;
+		if ( ((ix-ibndx)<0) || ((ix-ibndx)>mod.nx) ) continue; /* source outside grid */
+        if ( ((iy-ibndy)<0) || ((iy-ibndy)>mod.ny) ) continue; /* source outside grid */
+
+		if (verbose>=4 && itime==0) {
+			vmess("Source %li positioned at grid ix=%li iy=%li iz=%li",isrc, ix, iy, iz);
+		}
+
+		/* cosine squared windowing to reduce edge effects on shot arrays */
+		if ( (src.n>1) && src.window) {
+            scl = 1.0;
+			if (isrc < src.window) {
+				scl = cos(0.5*M_PI*(src.window - isrc)/src.window);
+			}
+			else if (isrc > src.n-src.window+1) {
+				scl = cos(0.5*M_PI*(src.window - (src.n-isrc+1))/src.window);
+			}
+			src_ampl *= scl*scl;
+		}
+
+		/* source scaling factor to compensate for discretisation */
+
+		/* old amplitude setting does not obey reciprocity */
+		// src_ampl *= rox[ix*n1+iz]*l2m[ix*n1+iz]/(dt);
+
+/* in older version added factor 2.0 to be compliant with defined Green's functions in Marchenko algorithm */
+/* this is now set to 1.0 */
+		src_ampl *= (1.0/mod.dx)*l2m[iy*n1*n2+ix*n1+iz];
+
+		if (verbose>5) {
+			vmess("Source %li at grid [ix=%li,iy=%li,iz=%li] at itime %li has value %e",isrc, ix, iy, iz, itime, src_ampl);
+		}
+
+		/* Force source */
+
+		if (src.type == 6) {
+			vx[iy*n1*n2+ix*n1+iz] += src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+			/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
+			//vx[ix*n1+iz] = 0.5*(vx[(ix+1)*n1+iz]+vx[(ix-1)*n1+iz])+src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
+		}
+		else if (src.type == 7) {
+			vz[iy*n1*n2+ix*n1+iz] += src_ampl*roz[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+			/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
+			/* stable implementation changes amplitude and more work is needed */
+			//vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
+			//vz[ix*n1+iz] = 0.25*(vz[ix*n1+iz-2]+vz[ix*n1+iz-1]+vz[ix*n1+iz]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
+        } /* src.type */
+
+        
+		/* Stress source */
+
+		if (mod.ischeme <= 2) { /* Acoustic scheme */
+			/* Compressional source */
+			if (src.type == 1) {
+				if (src.orient != 1) src_ampl=src_ampl/mod.dx;
+
+				if (src.orient==1) { /* monopole */
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+				}
+				else if (src.orient==2) { /* dipole +/- */
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+					tzz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+				}
+				else if (src.orient==3) { /* dipole - + */
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+					tzz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+				}
+				else if (src.orient==4) { /* dipole +/0/- */
+					if (iz > ibndz) 
+						tzz[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
+					if (iz < mod.nz+ibndz-1) 
+						tzz[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
+				}
+				else if (src.orient==5) { /* dipole + - */
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+					tzz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+				}
+			}
+		}
+		else { /* Elastic scheme */
+			/* Compressional source */
+			if (src.type == 1) {
+				if (src.orient==1) { /* monopole */
+					txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+				}
+				else if (src.orient==2) { /* dipole +/- */
+					txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+					txx[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+					tzz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+				}
+				else if (src.orient==3) { /* dipole - + */
+					txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+					txx[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+					tzz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+				}
+				else if (src.orient==4) { /* dipole +/0/- */
+					if (iz > ibndz) {
+						txx[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
+						tzz[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
+					}
+					if (iz < mod.nz+ibndz-1) {
+						txx[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
+						tzz[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
+					}
+				}
+				else if (src.orient==5) { /* dipole + - */
+					txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+					tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+					txx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+					tzz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+				}
+			}
+			else if (src.type == 2) {
+				/* Txz source */
+				if ((iz == ibndz) && bnd.top==1) {
+					txz[iy*n1*n2+(ix-1)*n1+iz-1] += src_ampl;
+					txz[iy*n1*n2+ix*n1+iz-1] += src_ampl;
+				}
+				else {
+					txz[iy*n1*n2+ix*n1+iz] += src_ampl;
+				}
+				/* possible dipole orientations for a txz source */
+				if (src.orient == 2) { /* dipole +/- */
+					txz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+				}
+				else if (src.orient == 3) { /* dipole - + */
+					txz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+				}
+				else if (src.orient == 4) { /*  dipole +/O/- */
+					/* correction: subtrace previous value to prevent z-1 values. */
+					txz[iy*n1*n2+ix*n1+iz] -= 2.0*src_ampl;
+					txz[iy*n1*n2+ix*n1+iz+1] += src_ampl;
+				}
+				else if (src.orient == 5) { /* dipole + - */
+					txz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+				}
+			}
+			/* Tzz source */
+			else if(src.type == 3) {
+				tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+			} 
+			/* Txx source */
+			else if(src.type == 4) {
+				txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+			} 
+
+/***********************************************************************
+* pure potential shear S source (experimental)
+* Curl S-pot = CURL(F) = dF_x/dz - dF_z/dx
+***********************************************************************/
+			else if(src.type == 5) {
+				src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+				if (src.orient == 3) src_ampl = -src_ampl;
+                /* first order derivatives */
+				vx[iy*n1*n2+ix*n1+iz]         += src_ampl*sdx;
+				vx[(iy-1)*n1*n2+ix*n1+iz-1]   -= src_ampl*sdx;
+				vy[iy*n1*n2+ix*n1+iz]         += src_ampl*sdx;
+				vy[iy*n1*n2+(ix-1)*n1+iz-1]   -= src_ampl*sdx;
+				vz[iy*n1*n2+ix*n1+iz]         -= src_ampl*sdx;
+				vz[(iy-1)*n1*n2+(ix-1)*n1+iz] += src_ampl*sdx;
+                
+                /* second order derivatives */
+                /*
+				vx[ix*n1+iz]     += c1*src_ampl*sdx;
+                vx[ix*n1+iz-1]   -= c1*src_ampl*sdx;
+				vx[ix*n1+iz+1]   += c2*src_ampl*sdx;
+                vx[ix*n1+iz-2]   -= c2*src_ampl*sdx;
+
+                vz[ix*n1+iz]     -= c1*src_ampl*sdx;
+				vz[(ix-1)*n1+iz] += c1*src_ampl*sdx;
+				vz[(ix+1)*n1+iz] -= c2*src_ampl*sdx;
+				vz[(ix-2)*n1+iz] += c2*src_ampl*sdx;
+                 */
+
+				/* determine second position of dipole */
+				if (src.orient == 2) { /* dipole +/- vertical */
+					iz += 1;
+                    vx[iy*n1*n2+ix*n1+iz]         -= src_ampl*sdx;
+                    vx[(iy-1)*n1*n2+ix*n1+iz-1]   += src_ampl*sdx;
+                    vy[iy*n1*n2+ix*n1+iz]         -= src_ampl*sdx;
+                    vy[iy*n1*n2+(ix-1)*n1+iz-1]   += src_ampl*sdx;
+                    vz[iy*n1*n2+ix*n1+iz]         += src_ampl*sdx;
+                    vz[(iy-1)*n1*n2+(ix-1)*n1+iz] -= src_ampl*sdx;
+				}
+				else if (src.orient == 3) { /* dipole - + horizontal */
+					ix += 1;
+                    vx[iy*n1*n2+ix*n1+iz]         -= src_ampl*sdx;
+                    vx[(iy-1)*n1*n2+ix*n1+iz-1]   += src_ampl*sdx;
+                    vy[iy*n1*n2+ix*n1+iz]         -= src_ampl*sdx;
+                    vy[iy*n1*n2+(ix-1)*n1+iz-1]   += src_ampl*sdx;
+                    vz[iy*n1*n2+ix*n1+iz]         += src_ampl*sdx;
+                    vz[(iy-1)*n1*n2+(ix-1)*n1+iz] -= src_ampl*sdx;
+				}
+            }
+/***********************************************************************
+* pure potential pressure P source (experimental)
+* Divergence P-pot = DIV(F) = dF_x/dx + dF_z/dz
+***********************************************************************/
+            else if(src.type == 8) {
+			    src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+                if (src.orient == 3) src_ampl = -src_ampl;
+                vx[iy*n1*n2+(ix+1)*n1+iz] += src_ampl*sdx;
+                vx[iy*n1*n2+ix*n1+iz]     -= src_ampl*sdx;
+                vy[(iy+1)*n1*n2+ix*n1+iz] += src_ampl*sdx;
+                vy[iy*n1*n2+ix*n1+iz]     -= src_ampl*sdx;
+                vz[iy*n1*n2+ix*n1+iz+1]   += src_ampl*sdx;
+                vz[iy*n1*n2+ix*n1+iz]     -= src_ampl*sdx;
+                /* determine second position of dipole */
+                if (src.orient == 2) { /* dipole +/- */
+                    iz += 1;
+                    vx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl*sdx;
+                    vx[iy*n1*n2+ix*n1+iz]     += src_ampl*sdx;
+                    vy[(iy+1)*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+                    vy[iy*n1*n2+ix*n1+iz]     += src_ampl*sdx;
+                    vz[iy*n1*n2+ix*n1+iz+1]   -= src_ampl*sdx;
+                    vz[iy*n1*n2+ix*n1+iz]     += src_ampl*sdx;
+                }
+                else if (src.orient == 3) { /* dipole - + */
+                    ix += 1;
+                    vx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl*sdx;
+                    vx[iy*n1*n2+ix*n1+iz]     += src_ampl*sdx;
+                    vy[(iy+1)*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+                    vy[iy*n1*n2+ix*n1+iz]     += src_ampl*sdx;
+                    vz[iy*n1*n2+ix*n1+iz+1]   -= src_ampl*sdx;
+                    vz[iy*n1*n2+ix*n1+iz]     += src_ampl*sdx;
+                }
+			}
+            else if(src.type == 9) {
+				rake = 0.5*M_PI;
+				Mxx = -1.0*(sin(src.dip)*cos(rake)*sin(2.0*src.strike)+sin(src.dip*2.0)*sin(rake)*sin(src.strike)*sin(src.strike));
+				Mxz = -1.0*(cos(src.dip)*cos(rake)*cos(src.strike)+cos(src.dip*2.0)*sin(rake)*sin(src.strike));
+				Mzz = sin(src.dip*2.0)*sin(rake);
+
+				txx[iy*n1*n2+ix*n1+iz] -= Mxx*src_ampl;
+				tzz[iy*n1*n2+ix*n1+iz] -= Mzz*src_ampl;
+				txz[iy*n1*n2+ix*n1+iz] -= Mxz*src_ampl;
+			} /* src.type */
+		} /* ischeme */
+	} /* loop over isrc */
+
+	return 0;
+}

fdelmodc3D/fdelmodc3D.c

@@ -589,51 +589,41 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
 //						vx, vz, tzz, rox, roz, l2m, verbose);
 //					break;
 				case -1 : /* Acoustic dissipative media FD kernel */
-					acoustic4_qr(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-						vx, vz, tzz, rox, roz, l2m, verbose);
+					vmess("Acoustic dissipative not yet available");
 					break;
 				case 1 : /* Acoustic FD kernel */
 					if (mod.iorder==2) {
-						acoustic2(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-							vx, vz, tzz, rox, roz, l2m, verbose);
+						vmess("Acoustic order 2 not yet available");
 					}
 					else if (mod.iorder==4) {
                         if (mod.sh) {
-                            acousticSH4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-                                  vx, vz, tzz, rox, roz, l2m, verbose);
+                            vmess("SH order 4 not yet available");
                         }
                         else {
-                            acoustic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-                                      vx, vz, tzz, rox, roz, l2m, verbose);
+							acoustic4_3D(mod, src, wav, bnd, it, ixsrc, iysrc, izsrc, src_nwav,
+									vx, vy, vz, tzz, rox, roy, roz, l2m, verbose);
                         }
 					}
 					else if (mod.iorder==6) {
-						acoustic6(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-							vx, vz, tzz, rox, roz, l2m, verbose);
+						vmess("Acoustic order 6 not yet available");
 					}
 					break;
 				case 2 : /* Visco-Acoustic FD kernel */
-					viscoacoustic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-							vx, vz, tzz, rox, roz, l2m, tss, tep, q, verbose);
+						vmess("Visco-Acoustic order 4 not yet available");
 					break;
 				case 3 : /* Elastic FD kernel */
                     if (mod.iorder==4) {
-                        elastic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-                            vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
+						vmess("Elastic order 4 not yet available");
 					}
 					else if (mod.iorder==6) {
-                        elastic6(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-							vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
+						vmess("Elastic order 6 not yet available");
                     }
 					break;
 				case 4 : /* Visco-Elastic FD kernel */
-					viscoelastic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-						vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, 
-						tss, tep, tes, r, q, p, verbose);
+						vmess("Visco-Elastic order 4 not yet available");
 					break;
 				case 5 : /* Elastic FD kernel with S-velocity set to zero*/
-                     elastic4dc(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav, 
-                            vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
+						vmess("DC-Elastic order 4 not yet available");
 					break;
 			}
 

fdelmodc3D/getRecTimes3D.c

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.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 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)
+{
+	int n1, n2, ibndx, ibndy, ibndz;
+	int 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;
+
+    ibndx = mod.ioPx;
+    ibndy = mod.ioPy;
+    ibndz = mod.ioPz;
+    if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+    if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
+    if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+	n1    = mod.naz;
+    n2    = mod.nax;
+	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;
+		iy = rec.y[irec]+ibndy;
+		ix = rec.x[irec]+ibndx;
+		iz1 = iz-1;
+		iy1 = iy-1;
+		ix1 = ix-1;
+		iz2 = iz+1;
+		iy2 = iy+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;
+			iy = (int)floorf(rec.yr[irec]/mod.dy)+ibndy;
+			ix = (int)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;
+			iz1 = iz-1;
+			iy1 = iy-1;
+			ix1 = ix-1;
+			iz2 = iz+1;
+			iy2 = iy+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 */
+				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*n1+iz+1];
+				C011 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz];
+				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;
+			}
+			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;
+}

fdelmodc3D/writeSrcRecPos3D.c

+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"par.h"
+#include"fdelmodc3D.h"
+
+#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))
+
+/**
+* Writes the source and receiver positions into a gridded file,
+* which has the same size as the input gridded model files. 
+* Source positions have a value +1 and receivers -1.
+*
+*   AUTHOR:
+*           Jan Thorbecke (janth@xs4all.nl)
+*           The Netherlands 
+**/
+
+long writesufile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2);
+
+long writeSrcRecPos3D(modPar *mod, recPar *rec, srcPar *src, shotPar *shot)
+{
+	FILE *fp;
+	float *dum, sub_x0, sub_y0, sub_z0, dx, dy, dz;
+	long is, nx, ny, nz, is0, ish, ix, iy, iz, ndot, idx, idy, idz;
+	char tmpname[1024];
+
+ 	ndot = 2;
+	nx = mod->nx;
+    ny = mod->ny;
+	nz = mod->nz;
+	dx = mod->dx;
+    dy = mod->dy;
+	dz = mod->dz;
+	sub_x0 = mod->x0;
+	sub_y0 = mod->y0;
+	sub_z0 = mod->z0;
+
+	/* write velocity field with positions of the sources */
+	dum = (float *)calloc(nx*ny*nz, sizeof(float));
+	vmess("Positions: shot=%li src=%li rec=%li", shot->n, src->n, rec->n);
+	/* source positions for random shots */
+	if (src->random) {
+		sprintf(tmpname,"SrcPositions%li.txt",src->n);
+		fp = fopen(tmpname, "w+");
+		for (is=0; is<src->n; is++) {
+            for (idy=0; idy<=ndot; idy++) {
+                for (idx=0; idx<=ndot; idx++) {
+                    for (idz=0; idz<=ndot; idz++) {
+                        dum[(MAX(0,src->y[is]-idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+                        dum[(MAX(0,src->y[is]-idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+                        dum[(MAX(0,src->y[is]-idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+                        dum[(MAX(0,src->y[is]-idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+                        dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+                        dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+                        dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+                        dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+                    }
+                }
+            }
+			fprintf(fp, "%f %f %f\n", src->z[is]*dz+sub_z0, src->y[is]*dy+sub_y0, src->x[is]*dx+sub_x0);
+		}
+		fclose(fp);
+	}
+	/* source positions for single shot sources with plane waves */
+	else if (src->plane) {
+    	is0 = -1*floor((src->n-1)/2);
+		sprintf(tmpname,"SrcPositions%li.txt",shot->n);
+		fp = fopen(tmpname, "w+");
+		for (ish=0; ish<shot->n; ish++) {
+			for (is=0; is<src->n; is++) {
+				ix = shot->x[ish] + 1 + is0 + is;
+				iy = shot->y[ish] + 1 + is0 + is;
+				iz = shot->z[ish] + 1;
+				dum[iy*nx*nz+ix*nz+iz] = 1.0;
+
+                dum[(MAX(0,iy-1))*nx*nz+ix*nz+iz] = 1.0;
+				dum[(MIN(ny-1,iy+1))*nx*nz+ix*nz+iz] = 1.0;
+				dum[iy*nx*nz+(MAX(0,ix-1))*nz+iz] = 1.0;
+				dum[iy*nx*nz+(MIN(nx-1,ix+1))*nz+iz] = 1.0;
+				dum[iy*nx*nz+ix*nz+MAX(0,iz-1)] = 1.0;
+				dum[iy*nx*nz+ix*nz+MIN(nz-1,iz+1)] = 1.0;
+				fprintf(fp, "(%f, %f, %f)\n", ix*dx+sub_x0, iy*dy+sub_y0, iz*dz+sub_z0);
+			}
+		}
+		fclose(fp);
+	}
+	else if (src->multiwav) {
+	/* source positions for single shot sources with multiple wavelets */
+		sprintf(tmpname,"SrcPositions%li.txt",shot->n);
+		fp = fopen(tmpname, "w+");
+		for (ish=0; ish<shot->n; ish++) {
+			for (is=0; is<src->n; is++) {
+				ix = src->x[is];
+                iy = src->y[is];
+				iz = src->z[is];
+				dum[iy*nx*nz+ix*nz+iz] = 1.0;
+
+                dum[(MAX(0,iy-1))*nx*nz+ix*nz+iz] = 1.0;
+				dum[(MIN(ny-1,iy+1))*nx*nz+ix*nz+iz] = 1.0;
+				dum[iy*nx*nz+(MAX(0,ix-1))*nz+iz] = 1.0;
+				dum[iy*nx*nz+(MIN(nx-1,ix+1))*nz+iz] = 1.0;
+				dum[iy*nx*nz+ix*nz+MAX(0,iz-1)] = 1.0;
+				dum[iy*nx*nz+ix*nz+MIN(nz-1,iz+1)] = 1.0;
+				fprintf(fp, "(%f, %f, %f)\n", ix*dx+sub_x0, iy*dy+sub_y0, iz*dz+sub_z0);
+			}
+		}
+		fclose(fp);
+	}
+	else {
+		sprintf(tmpname,"SrcPositions%li.txt",shot->n);
+		fp = fopen(tmpname, "w+");
+		for (is=0; is<shot->n; is++) {
+            for (idy=0; idy<=ndot; idy++) {
+                for (idx=0; idx<=ndot; idx++) {
+                    for (idz=0; idz<=ndot; idz++) {
+                        dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+                        dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+                        dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+                        dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+                        dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+                        dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+                        dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+                        dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+                    }
+                }
+            }
+			fprintf(fp, "%f %f %f\n", shot->z[is]*dz+sub_z0, shot->y[is]*dy+sub_y0, shot->x[is]*dx+sub_x0);
+		}
+		fclose(fp);
+	}
+
+	/* receiver positions */
+	sprintf(tmpname,"RcvPositions%li.txt",rec->n);
+	fp = fopen(tmpname, "w+");
+	for (is=0; is<rec->n; is++) {
+		dum[rec->y[is]*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
+		dum[(MAX(0,rec->y[is]-1))*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
+		dum[(MIN(ny-1,rec->y[is]+1))*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
+		dum[rec->y[is]*nx*nz+(MAX(0,rec->x[is]-1))*nz+rec->z[is]] = -1.0;
+		dum[rec->y[is]*nx*nz+(MIN(nx-1,rec->x[is]+1))*nz+rec->z[is]] = -1.0;
+		dum[rec->y[is]*nx*nz+rec->x[is]*nz+MAX(0,rec->z[is]-1)] = -1.0;
+		dum[rec->y[is]*nx*nz+rec->x[is]*nz+MIN(nz-1,rec->z[is]+1)] = -1.0;
+
+		if (rec->int_vx==3) {
+			fprintf(fp, "(%f, %f, %f)\n", rec->xr[is]*dx+sub_x0, rec->yr[is]*dy+sub_y0, rec->zr[is]*dz+sub_z0);
+		}
+		else {
+			fprintf(fp, "(%f, %f, %f)\n", rec->x[is]*dx+sub_x0, rec->y[is]*dy+sub_y0, rec->z[is]*dz+sub_z0);
+		}
+	}
+	fclose(fp);
+	writesufile3D("SrcRecPositions.su", dum, nz, nx*ny, sub_z0, sub_x0, dz, dx);
+	free(dum);
+
+	return 0;
+}

fdelmodc3D/writesufile3D.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <string.h>
+#include "par.h"
+#include "fdelmodc3D.h"
+#include "SUsegy.h"
+#include "segy.h"
+
+/**
+*  Writes an 2D array to a SU file
+*
+*   AUTHOR:
+*           Jan Thorbecke (janth@xs4all.nl)
+*           The Netherlands 
+**/
+
+#define TRCBYTES 240
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define ISODD(n) ((n) & 01)
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+long writesufile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2)
+{
+	FILE    *file_out;
+	size_t  nwrite, itrace;
+	long     ns;
+	segy    *hdr;
+
+/* Read in parameters */
+
+
+	
+	if (n1 > USHRT_MAX) {
+		vwarn("Output file %s: number of samples is truncated from %li to USHRT_MAX.", filename, n1);
+	}
+	ns = MIN(n1,USHRT_MAX);
+
+	file_out = fopen( filename, "w+" );
+	assert( file_out );
+
+	hdr = (segy *)calloc(1,TRCBYTES);
+	hdr->ns = ns;
+	hdr->dt = NINT(1000000*(d1));
+	hdr->d1 = d1;
+	hdr->d2 = d2;
+	hdr->f1 = f1;
+	hdr->f2 = f2;
+	hdr->fldr = 1;
+	hdr->trwf = n2;
+
+	for (itrace=0; itrace<n2; itrace++) {
+		hdr->tracl = itrace+1;
+		nwrite = fwrite( hdr, 1, TRCBYTES, file_out );
+		assert (nwrite == TRCBYTES);
+		nwrite = fwrite( &data[itrace*n1], sizeof(float), ns, file_out );
+		assert (nwrite == ns);
+	} 
+	fclose(file_out);
+	free(hdr);
+
+	return 0;
+}
+
+/**
+*  Writes an 2D array to a SU file
+*  special routine for src_nwav array which has a different number of samples for each shot 
+*
+**/
+
+long writesufilesrcnwav3D(char *filename, float **src_nwav, wavPar wav, long n1, long n2, float f1, float f2, float d1, float d2)
+{
+	FILE    *file_out;
+	size_t  nwrite, itrace;
+	float   *trace;
+	long     ns;
+	segy    *hdr;
+
+/* Read in parameters */
+
+	if (n1 > USHRT_MAX) {
+		vwarn("Output file %s: number of samples is truncated from %li to USHRT_MAX.", filename, n1);
+	}
+	ns = MIN(n1,USHRT_MAX);
+
+	file_out = fopen( filename, "w+" );
+	assert( file_out );
+
+	trace = (float *)malloc(n1*sizeof(float));
+	hdr = (segy *)calloc(1,TRCBYTES);
+	hdr->ns = ns;
+	hdr->dt = NINT(1000000*(d1));
+	hdr->d1 = d1;
+	hdr->d2 = d2;
+	hdr->f1 = f1;
+	hdr->f2 = f2;
+	hdr->fldr = 1;
+	hdr->trwf = n2;
+
+	for (itrace=0; itrace<n2; itrace++) {
+		hdr->tracl = itrace+1;
+		nwrite = fwrite( hdr, 1, TRCBYTES, file_out );
+		assert (nwrite == TRCBYTES);
+		memset(trace, 0, n1*sizeof(float));
+		memcpy(trace, &src_nwav[itrace][0], wav.nsamp[itrace]*sizeof(float));
+		nwrite = fwrite( &trace[0], sizeof(float), ns, file_out );
+		assert (nwrite == ns);
+	} 
+	fclose(file_out);
+	free(hdr);
+	free(trace);
+
+	return 0;
+}
+
+/**
+*  Writes an 2D array to a SU file
+*  special routine which used segyhdrs which have ns defined as integer (32 bit)
+*  to handle more than 2^16 samples per trace.
+*
+**/
+
+long writeSUfile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2)
+{
+	FILE    *file_out;
+	size_t  nwrite, itrace;
+	SUsegy  *SUhdr;
+	char    *ptr;
+
+/* Read in parameters */
+
+    ptr = strstr(filename, " ");
+    *ptr = '\0';
+
+	file_out = fopen( filename, "w+" );
+	assert( file_out );
+
+	SUhdr = (SUsegy *)calloc(1,TRCBYTES);
+	SUhdr->ns = n1;
+	SUhdr->dt = NINT(1000000*(d1));
+	SUhdr->d1 = d1;
+	SUhdr->d2 = d2;
+	SUhdr->f1 = f1;
+	SUhdr->f2 = f2;
+	SUhdr->fldr = 1;
+	SUhdr->trwf = n2;
+
+	for (itrace=0; itrace<n2; itrace++) {
+		SUhdr->tracl = itrace+1;
+		nwrite = fwrite( SUhdr, 1, TRCBYTES, file_out );
+		assert (nwrite == TRCBYTES);
+		nwrite = fwrite( &data[itrace*n1], sizeof(float), n1, file_out );
+		assert (nwrite == n1);
+	} 
+	fclose(file_out);
+	free(SUhdr);
+
+	return 0;
+}
+