From fa095294a6cc84479c4840c75e7f9b514b8d4179 Mon Sep 17 00:00:00 2001
From: Jan at TU-Delft <J.W.Thorbecke@tudelft.nl>
Date: Thu, 18 Jan 2018 13:36:36 +0100
Subject: [PATCH] adding vidale raytracing method and bug fix in Jesper's
---
raytime/JespersRayTracer.c | 12 -
raytime/Makefile | 6 +
raytime/raytime.c | 101 +++---
raytime/raytime.h | 32 ++
raytime/vidale.c | 607 +++++++++++++++++++++++++++++++++++++
5 files changed, 709 insertions(+), 49 deletions(-)
create mode 100644 raytime/vidale.c
diff --git a/raytime/JespersRayTracer.c b/raytime/JespersRayTracer.c
index 255c4d8..bedb681 100644
--- a/raytime/JespersRayTracer.c
+++ b/raytime/JespersRayTracer.c
@@ -18,18 +18,6 @@
static float H, L, W, iH, iL, iW;
-typedef struct _icoord { /* 3D coordinate integer */
- int z;
- int x;
- int y;
-} icoord;
-
-typedef struct _fcoord { /* 3D coordinate float */
- float z;
- float x;
- float y;
-} fcoord;
-
int getnRay(icoord size, fcoord s, fcoord r, float dx, int nRayStep);
int traceTwoPoint(fcoord s, fcoord r, int nRay, fcoord *rayReference3D);
float takeOffAngle(fcoord s, fcoord r);
diff --git a/raytime/Makefile b/raytime/Makefile
index 6d454e7..2a4d790 100644
--- a/raytime/Makefile
+++ b/raytime/Makefile
@@ -9,17 +9,23 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#LIBS += -L$L -lgenfft -lm -lc
#OPTC = -g -Wall -fsignaling-nans -O0
#OPTC += -fopenmp -Waddress
+OPTC += -g -O0
#OPTC := $(subst -O3 -ffast-math, -O1 -g ,$(OPTC))
#PGI options for compiler feedback
#OPTC += -Mprof=lines
#LDFLAGS += -Mprof=lines
+# side.c \
+# corner.c \
+# near_source.c \
+
all: raytime
PRG = raytime
SRCC = $(PRG).c \
JespersRayTracer.c \
+ vidale.c \
getParameters.c \
getWaveletInfo.c \
writeSrcRecPos.c \
diff --git a/raytime/raytime.c b/raytime/raytime.c
index cf010ac..eea8786 100644
--- a/raytime/raytime.c
+++ b/raytime/raytime.c
@@ -11,18 +11,6 @@
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-typedef struct _icoord { /* 3D coordinate integer */
- int z;
- int x;
- int y;
-} icoord;
-
-typedef struct _fcoord { /* 3D coordinate float */
- float z;
- float x;
- float y;
-} fcoord;
-
double wallclock_time(void);
void name_ext(char *filename, char *extension);
@@ -41,6 +29,8 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
+void vidale(float *ttime, float *slow, icoord *isrc, icoord grid, float dx, int order, int mzrcv);
+
/* Self documentation */
char *sdoc[] = {
" ",
@@ -54,12 +44,13 @@ char *sdoc[] = {
" dz= ............... read from model file: if dz==0 then dz= can be used to set it",
"" ,
" RAY TRACING PARAMETERS:",
+" method=jesper ..... calculation method (jesper, fd) ",
" smoothwindow=0 .... if set lenght of 2/3D smoothing window on slowness",
" useT2=0 ........... 1: compute more accurate T2 pertubation correction",
" geomspread=1 ...... 1: compute Geometrical Spreading Factor",
" nraystep=5 ........ number of points on ray",
+" sbox=1 ............ radius of inner straight ray (fd method)",
" OPTIONAL PARAMETERS:",
-" ischeme=3 ......... 1=acoustic, 2=visco-acoustic 3=elastic, 4=visco-elastic",
" sinkdepth=0 ....... receiver grid points below topography (defined bij cp=0.0)",
" sinkdepth_src=0 ... source grid points below topography (defined bij cp=0.0)",
" sinkvel=0 ......... use velocity of first receiver to sink through to next layer",
@@ -122,16 +113,17 @@ int main(int argc, char **argv)
srcPar src;
shotPar shot;
rayPar ray;
- float *velocity, *slowness, *smooth;
+ float *velocity, *slowness, *smooth, *trueslow, **inter;
double t0, t1, t2, tinit, tray, tio;
size_t size;
int nw, n1, ix, iz, ir, ixshot, izshot;
- int irec;
+ int irec, sbox, ipos;
fcoord coordsx, coordgx, Time;
- icoord grid;
- float Jr, *ampl, *time;
+ icoord grid, isrc;
+ float Jr, *ampl, *time, *ttime, *ttime_p;
+ float dxrcv, dzrcv;
segy hdr;
- char filetime[1024], fileamp[1024];
+ char filetime[1024], fileamp[1024], *method;
size_t nwrite;
int verbose;
FILE *fpt, *fpa;
@@ -141,19 +133,27 @@ int main(int argc, char **argv)
requestdoc(0);
if (!getparint("verbose",&verbose)) verbose=0;
+ if(!getparint("sbox", &sbox)) sbox = 1;
+ if(!getparstring("method", &method)) method="jesper";
getParameters(&mod, &rec, &src, &shot, &ray, verbose);
/* allocate arrays for model parameters: the different schemes use different arrays */
n1 = mod.nz;
- nw = ray.smoothwindow;
+ if(!strcmp(method,"fd")) nw = 0;
+ else nw = ray.smoothwindow;
velocity = (float *)calloc(mod.nx*mod.nz,sizeof(float));
slowness = (float *)calloc((mod.nx+2*nw)*(mod.nz+2*nw),sizeof(float));
-// slowness = (float *)calloc(mod.nx*mod.nz,sizeof(float));
+ trueslow = (float *)calloc(mod.nx*mod.nz,sizeof(float));
- /* read velocity and density files */
+ if(!strcmp(method,"fd")) {
+ ttime = (float *)calloc(mod.nx*mod.nz,sizeof(float));
+ for(ir=0, ttime_p=ttime; ir<mod.nx*mod.nz; ir++, ttime_p++)
+ *ttime_p = Infinity;
+ }
+ /* read velocity and density files */
readModel(mod, velocity, slowness, nw);
/* allocate arrays for wavefield and receiver arrays */
@@ -198,9 +198,9 @@ int main(int argc, char **argv)
grid.x = mod.nx;
grid.z = mod.nz;
grid.y = 1;
- if ( (ray.smoothwindow) != 0 ) { /* smooth slowness */
+ if ( nw != 0 ) { /* smooth slowness */
smooth = (float *)calloc(grid.x*grid.z,sizeof(float));
- applyMovingAverageFilter(slowness, grid, ray.smoothwindow, 2, smooth);
+ applyMovingAverageFilter(slowness, grid, nw, 2, smooth);
memcpy(slowness,smooth,grid.x*grid.z*sizeof(float));
free(smooth);
}
@@ -218,6 +218,7 @@ int main(int argc, char **argv)
assert(fpa != NULL);
}
+ memset(&hdr,0,sizeof(hdr));
hdr.dt = (unsigned short)1;
hdr.scalco = -1000;
hdr.scalel = -1000;
@@ -244,25 +245,50 @@ int main(int argc, char **argv)
vmess(" which are actual positions z=%.2f - %.2f", mod.z0+rec.zr[0], mod.z0+rec.zr[rec.n-1]);
}
- coordsx.x = mod.x0+shot.x[ixshot]*mod.dx;
- coordsx.z = mod.z0+shot.z[izshot]*mod.dz;
+ coordsx.x = shot.x[ixshot]*mod.dx;
+ coordsx.z = shot.z[izshot]*mod.dz;
coordsx.y = 0;
t1=wallclock_time();
tio += t1-t2;
+
+ if (!strcmp(method,"jesper")) {
#pragma omp parallel for default(shared) \
private (coordgx,irec,Time,Jr)
- for (irec=0; irec<rec.n; irec++) {
- coordgx.x=mod.x0+rec.xr[irec];
- coordgx.z=mod.z0+rec.zr[irec];
- coordgx.y = 0;
+ for (irec=0; irec<rec.n; irec++) {
+ coordgx.x=rec.xr[irec];
+ coordgx.z=rec.zr[irec];
+ coordgx.y = 0;
+
+ getWaveParameter(slowness, grid, mod.dx, coordsx, coordgx, ray, &Time, &Jr);
- getWaveParameter(slowness, grid, mod.dx, coordsx, coordgx, ray, &Time, &Jr);
-
- time[((izshot*shot.nx)+ixshot)*rec.n + irec] = Time.x + Time.y + Time.z;
- ampl[((izshot*shot.nx)+ixshot)*rec.n + irec] = Jr;
- if (verbose>4) vmess("shot=%f,%f receiver at %f,%f T0=%f T1=%f T2=%f Jr=%f",coordsx.x, coordsx.z, coordgx.x, coordgx.z, Time.x, Time.y, Time.z, Jr);
- }
+ time[((izshot*shot.nx)+ixshot)*rec.n + irec] = Time.x + Time.y + Time.z;
+ ampl[((izshot*shot.nx)+ixshot)*rec.n + irec] = Jr;
+ if (verbose>4) vmess("JS: shot=%f,%f receiver at %f,%f T0=%f T1=%f T2=%f Jr=%f",coordsx.x, coordsx.z, coordgx.x, coordgx.z, Time.x, Time.y, Time.z, Jr);
+ }
+ }
+ else if(!strcmp(method,"fd")) {
+ int mzrcv;
+
+ isrc.x = shot.x[ixshot];
+ isrc.y = 0;
+ isrc.z = shot.z[izshot];
+
+ mzrcv = 0;
+ for (irec = 0; irec < rec.n; irec++) mzrcv = MAX(rec.z[irec], mzrcv);
+
+ vidale(ttime,slowness,&isrc,grid,mod.dx,sbox, mzrcv);
+ for (irec=0; irec<rec.n; irec++) {
+ coordgx.x=mod.x0+rec.xr[irec];
+ coordgx.z=mod.z0+rec.zr[irec];
+ coordgx.y = 0;
+ ipos = ((izshot*shot.nx)+ixshot)*rec.n + irec;
+
+ time[ipos] = ttime[rec.z[irec]*mod.nx+rec.x[irec]];
+ ampl[ipos] = sqrt(time[ipos]);
+ if (verbose>4) vmess("FD: shot=%f,%f receiver at %f(%d),%f(%d) T=%e Ampl=%f",coordsx.x, coordsx.z, coordgx.x, rec.x[irec], coordgx.z, rec.z[irec], time[ipos], ampl[ipos]);
+ }
+ }
t2=wallclock_time();
tray += t2-t1;
@@ -273,9 +299,10 @@ private (coordgx,irec,Time,Jr)
hdr.tracl = ((izshot*shot.nx)+ixshot)+1;
hdr.tracf = ((izshot*shot.nx)+ixshot)+1;
hdr.ntr = shot.n;
- hdr.d1 = (rec.x[1]-rec.x[0])*mod.dx;
+ //hdr.d1 = (rec.x[1]-rec.x[0])*mod.dx; // discrete
+ hdr.d1 = (rec.xr[1]-rec.xr[0]);
hdr.f1 = mod.x0+rec.x[0]*mod.dx;
- hdr.d2 = (shot.x[1]-shot.x[0])*mod.dx;
+ hdr.d2 = (shot.x[MIN(shot.n-1,1)]-shot.x[0])*mod.dx;
hdr.f2 = mod.x0+shot.x[0]*mod.dx;
nwrite = fwrite( &hdr, 1, TRCBYTES, fpt);
diff --git a/raytime/raytime.h b/raytime/raytime.h
index 045e99e..fc38be6 100644
--- a/raytime/raytime.h
+++ b/raytime/raytime.h
@@ -1,7 +1,25 @@
#include<stdlib.h>
#include<stdio.h>
+#include<limits.h>
+#include<float.h>
#include<math.h>
+typedef struct _icoord { /* 3D coordinate integer */
+ int z;
+ int x;
+ int y;
+} icoord;
+
+typedef struct _fcoord { /* 3D coordinate float */
+ float z;
+ float x;
+ float y;
+} fcoord;
+
+struct s_ecount {
+ int corner,corner_min,side;
+};
+
typedef struct _receiverPar { /* Receiver Parameters */
char *file_rcv;
int n;
@@ -85,6 +103,20 @@ typedef struct _raypar { /* ray-tracing parameters */
int nray;
} rayPar;
+#ifndef TRUE
+# define TRUE 1
+#endif
+
+#ifndef FALSE
+# define FALSE 0
+#endif
+
+#define equal(x,y) !strcmp(x,y)
+#define min2(a,b) (((a) < (b)) ? (a) : (b))
+#define max2(a,b) (((a) > (b)) ? (a) : (b))
+
+#define Infinity FLT_MAX
+
#if __STDC_VERSION__ >= 199901L
/* "restrict" is a keyword */
#else
diff --git a/raytime/vidale.c b/raytime/vidale.c
new file mode 100644
index 0000000..a347169
--- /dev/null
+++ b/raytime/vidale.c
@@ -0,0 +1,607 @@
+/* FILE: vidale.c
+ * AUTHOR: Joseph R. Matarese
+ * DATE: Decenber 18, 1992
+ * Copyright (c) 1992 Joseph R. Matarese and Massachusetts Institute of
+ * Technology
+ *
+ * Permission is granted to copy or modify this code under the condition
+ * that the above copyright notice is maintained. The author and MIT
+ * make ABSOLUTELY NO WARRANTY regarding the fitness of this code for
+ * any purpose.
+ *
+ */
+
+//#include <raytime2.h>
+#include <raytime.h>
+
+extern void near_source(float *ttime, float *slow, icoord *nm, icoord *isrc, fcoord *scale, int order);
+
+extern void corner(float *ttime, float *slow,icoord *nm, icoord *iin, icoord *iout, fcoord *scale, struct s_ecount *ecount);
+
+extern void side(float *ttime, float *slow, icoord *nm, icoord *iin, icoord *iout, fcoord *scale, struct s_ecount *ecount);
+
+void rm_head(float *slow, icoord *ndim, icoord *isrc, int mzrcv, float dx, int *nzm);
+
+void vidale(float *ttime, float *slow1, icoord *isrc, icoord grid, float dx, int order, int mzrcv)
+{
+ int iz_lo, iz_hi, ix_lo, ix_hi, iz, ix, node_src;
+ int nzm;
+ icoord iin, iout, dnm, *nm;
+ fcoord dscale, *scale;
+ short ascending, finished;
+ struct s_ecount ecount;
+ float sx, sz, sign, *trueslow, *slow;
+
+ dscale.x = dx; dscale.y = 0.; dscale.z = dx;
+ scale = &dscale;
+ dnm.x = grid.x; dnm.y = 1; dnm.z = grid.z;
+ nm = &dnm;
+
+/* transpose slowness field for usage in vidale */
+ trueslow = (float *)calloc(grid.x*grid.z,sizeof(float));
+ for (ix=0; ix<grid.x; ix++) {
+ for (iz=0; iz<grid.z; iz++) {
+ trueslow[iz*grid.x + ix] = slow1[ix*grid.z + iz];
+ }
+ }
+ slow = trueslow;
+ rm_head(slow,&grid,isrc,mzrcv,dx,&nzm);
+
+ //for actual position of source not on grid: current no-use
+ node_src = isrc->z*grid.x + isrc->x;
+ sx = isrc->x*dx-isrc->x*dx;
+ sz = isrc->z*dx-isrc->z*dx;
+ if (sz < 0) sign = 1;
+ else sign = -1;
+ ttime[node_src] = sign*sqrt(sx*sx+sz*sz)*slow[node_src];
+
+ if(nm->y != 1)
+ verr("only 2D models implemented");
+
+ ecount.corner = ecount.corner_min = ecount.side = 0;
+
+/*---------------------------------------------------------------------------*
+ * Do near source region.
+ *---------------------------------------------------------------------------*/
+
+ near_source(ttime,slow,nm,isrc,scale,order);
+
+ iz_hi = min2(isrc->z+order,nm->z-1);
+ iz_lo = max2(isrc->z-order,0);
+ ix_hi = min2(isrc->x+order,nm->x-1);
+ ix_lo = max2(isrc->x-order,0);
+
+/*---------------------------------------------------------------------------*
+ * Loop over outer ring - verticals, then horizontals.
+ *---------------------------------------------------------------------------*/
+
+ finished = FALSE;
+ while(!finished) {
+
+ finished = TRUE;
+ if(ix_hi < nm->x-1) {
+ finished = FALSE;
+ iin.x = ix_hi; /* right-hand edge */
+ iout.x = ix_hi+1;
+
+ ascending = FALSE;
+ if(ttime[iz_lo*nm->x+iin.x] <= ttime[(iz_lo+1)*nm->x+iin.x]) {
+ ttime[iz_lo*nm->x+iout.x] = ttime[iz_lo*nm->x+iin.x] + 0.5 *
+ scale->x * (slow[iz_lo*nm->x+iout.x] + slow[iz_lo*nm->x+iin.x]);
+ ecount.corner_min++;
+ }
+ if(ttime[iz_lo*nm->x+iin.x] < ttime[(iz_lo+1)*nm->x+iin.x]) {
+ ascending = TRUE;
+ iin.z = iz_lo;
+ iout.z = iz_lo+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ for(iz=iz_lo+1; iz<=iz_hi-1; iz++) {
+ iin.z = iz;
+ if(ttime[iz*nm->x+iin.x] == ttime[(iz+1)*nm->x+iin.x]) {
+ iout.z = iz;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ if(ttime[iz*nm->x+iin.x] < ttime[(iz+1)*nm->x+iin.x]) {
+ if(!ascending) {
+ iout.z = iz;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ ascending = TRUE;
+ iout.z = iz+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ } else ascending = FALSE;
+ }
+
+ if(ttime[iz_hi*nm->x+iin.x] <= ttime[(iz_hi-1)*nm->x+iin.x]) {
+ ttime[iz_hi*nm->x+iout.x] = ttime[iz_hi*nm->x+iin.x] + 0.5 *
+ scale->x * (slow[iz_hi*nm->x+iout.x] + slow[iz_hi*nm->x+iin.x]);
+ ecount.corner_min++;
+ }
+ for(iz=iz_hi; iz>=iz_lo+1; iz--) {
+ if(ttime[iz*nm->x+iin.x] < ttime[(iz-1)*nm->x+iin.x]) {
+ iin.z = iz;
+ iout.z = iz-1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ }
+
+ }
+
+ if(ix_lo > 0) {
+ finished = FALSE;
+ iin.x = ix_lo; /* left-hand edge */
+ iout.x = ix_lo-1;
+
+ ascending = FALSE;
+ if(ttime[iz_lo*nm->x+iin.x] <= ttime[(iz_lo+1)*nm->x+iin.x]) {
+ ttime[iz_lo*nm->x+iout.x] = ttime[iz_lo*nm->x+iin.x] + 0.5 *
+ scale->x * (slow[iz_lo*nm->x+iout.x] + slow[iz_lo*nm->x+iin.x]);
+ ecount.corner_min++;
+ }
+ if(ttime[iz_lo*nm->x+iin.x] < ttime[(iz_lo+1)*nm->x+iin.x]) {
+ ascending = TRUE;
+ iin.z = iz_lo;
+ iout.z = iz_lo+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ for(iz=iz_lo+1; iz<=iz_hi-1; iz++) {
+ iin.z = iz;
+ if(ttime[iz*nm->x+iin.x] == ttime[(iz+1)*nm->x+iin.x]) {
+ iout.z = iz;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ if(ttime[iz*nm->x+iin.x] < ttime[(iz+1)*nm->x+iin.x]) {
+ if(!ascending) {
+ iout.z = iz;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ ascending = TRUE;
+ iout.z = iz+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ } else ascending = FALSE;
+ }
+
+ if(ttime[iz_hi*nm->x+iin.x] <= ttime[(iz_hi-1)*nm->x+iin.x]) {
+ ttime[iz_hi*nm->x+iout.x] = ttime[iz_hi*nm->x+iin.x] + 0.5 *
+ scale->x * (slow[iz_hi*nm->x+iout.x] + slow[iz_hi*nm->x+iin.x]);
+ ecount.corner_min++;
+ }
+ for(iz=iz_hi; iz>=iz_lo+1; iz--) {
+ if(ttime[iz*nm->x+iin.x] < ttime[(iz-1)*nm->x+iin.x]) {
+ iin.z = iz;
+ iout.z = iz-1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ }
+
+ }
+
+ if(iz_hi < nm->z-1) {
+ finished = FALSE;
+ iin.z = iz_hi; /* bottom edge */
+ iout.z = iz_hi+1;
+
+ ascending = FALSE;
+ if(ttime[iin.z*nm->x+ix_lo] <= ttime[iin.z*nm->x+ix_lo+1]) {
+ ttime[iout.z*nm->x+ix_lo] = ttime[iin.z*nm->x+ix_lo] + 0.5 *
+ scale->x * (slow[iout.z*nm->x+ix_lo] + slow[iin.z*nm->x+ix_lo]);
+ ecount.corner_min++;
+ }
+ if(ttime[iin.z*nm->x+ix_lo] < ttime[iin.z*nm->x+ix_lo+1]) {
+ ascending = TRUE;
+ iin.x = ix_lo;
+ iout.x = ix_lo+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ for(ix=ix_lo+1; ix<=ix_hi-1; ix++) {
+ iin.x = ix;
+ if(ttime[iin.z*nm->x+ix] == ttime[iin.z*nm->x+ix+1]) {
+ iout.x = ix;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ if(ttime[iin.z*nm->x+ix] < ttime[iin.z*nm->x+ix+1]) {
+ if(!ascending) {
+ iout.x = ix;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ ascending = TRUE;
+ iout.x = ix+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ } else ascending = FALSE;
+ }
+
+ if(ttime[iin.z*nm->x+ix_hi] <= ttime[iin.z*nm->x+ix_hi-1]) {
+ ttime[iout.z*nm->x+ix_hi] = ttime[iin.z*nm->x+ix_hi] + 0.5 *
+ scale->x * (slow[iout.z*nm->x+ix_hi] + slow[iin.z*nm->x+ix_hi]);
+ ecount.corner_min++;
+ }
+ for(ix=ix_hi; ix>=ix_lo+1; ix--) {
+ if(ttime[iin.z*nm->x+ix] < ttime[iin.z*nm->x+ix-1]) {
+ iin.x = ix;
+ iout.x = ix-1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ }
+
+ }
+
+ if(iz_lo > 0) {
+ finished = FALSE;
+ iin.z = iz_lo; /* top edge */
+ iout.z = iz_lo-1;
+
+ ascending = FALSE;
+ if(ttime[iin.z*nm->x+ix_lo] <= ttime[iin.z*nm->x+ix_lo+1]) {
+ ttime[iout.z*nm->x+ix_lo] = ttime[iin.z*nm->x+ix_lo] + 0.5 *
+ scale->x * (slow[iout.z*nm->x+ix_lo] + slow[iin.z*nm->x+ix_lo]);
+ ecount.corner_min++;
+ }
+ if(ttime[iin.z*nm->x+ix_lo] < ttime[iin.z*nm->x+ix_lo+1]) {
+ ascending = TRUE;
+ iin.x = ix_lo;
+ iout.x = ix_lo+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ for(ix=ix_lo+1; ix<=ix_hi-1; ix++) {
+ iin.x = ix;
+ if(ttime[iin.z*nm->x+ix] == ttime[iin.z*nm->x+ix+1]) {
+ iout.x = ix;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ if(ttime[iin.z*nm->x+ix] < ttime[iin.z*nm->x+ix+1]) {
+ if(!ascending) {
+ iout.x = ix;
+ side(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ ascending = TRUE;
+ iout.x = ix+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ } else ascending = FALSE;
+ }
+
+ if(ttime[iin.z*nm->x+ix_hi] <= ttime[iin.z*nm->x+ix_hi-1]) {
+ ttime[iout.z*nm->x+ix_hi] = ttime[iin.z*nm->x+ix_hi] + 0.5 *
+ scale->x * (slow[iout.z*nm->x+ix_hi] + slow[iin.z*nm->x+ix_hi]);
+ ecount.corner_min++;
+ }
+ for(ix=ix_hi; ix>=ix_lo+1; ix--) {
+ if(ttime[iin.z*nm->x+ix] < ttime[iin.z*nm->x+ix-1]) {
+ iin.x = ix;
+ iout.x = ix-1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+ }
+
+ }
+
+ if((ix_hi < nm->x-1) || (iz_hi < nm->z-1)) { /* bottom right corner */
+ iin.x = (ix_hi < nm->x-1) ? ix_hi : nm->x-2;
+ iout.x = iin.x+1;
+ iin.z = (iz_hi < nm->z-1) ? iz_hi : nm->z-2;
+ iout.z = iin.z+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+
+ if((ix_lo > 0) || (iz_hi < nm->z-1)) { /* bottom left corner */
+ iin.x = (ix_lo > 0) ? ix_lo : 1;
+ iout.x = iin.x-1;
+ iin.z = (iz_hi < nm->z-1) ? iz_hi : nm->z-2;
+ iout.z = iin.z+1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+
+ if((ix_hi < nm->x-1) || (iz_lo > 0)) { /* top right corner */
+ iin.x = (ix_hi < nm->x-1) ? ix_hi : nm->x-2;
+ iout.x = iin.x+1;
+ iin.z = (iz_lo > 0) ? iz_lo : 1;
+ iout.z = iin.z-1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+
+ if((ix_lo > 0) || (iz_lo > 0)) { /* top left corner */
+ iin.x = (ix_lo > 0) ? ix_lo : 1;
+ iout.x = iin.x-1;
+ iin.z = (iz_lo > 0) ? iz_lo : 1;
+ iout.z = iin.z-1;
+ corner(ttime,slow,nm,&iin,&iout,scale,&ecount);
+ }
+
+ ix_hi = min2(ix_hi+1,nm->x-1);
+ ix_lo = max2(ix_lo-1,0);
+ iz_hi = min2(iz_hi+1,nm->z-1);
+ iz_lo = max2(iz_lo-1,0);
+
+ } /* while(!finished) */
+
+/*
+ jm_message("debug",__FILE__,__LINE__,
+ "errors:\n %s: %d\n %s: %d\n %s: %d",
+ "corner minima",ecount.corner_min,
+ "corner negative sqrts",ecount.corner,
+ "side negative sqrts",ecount.side);
+ */
+ return;
+}
+
+
+
+/* FILE: near_source.c
+ * AUTHOR: Joseph R. Matarese
+ * DATE: Decenber 18, 1992
+ * Copyright (c) 1992 Joseph R. Matarese and Massachusetts Institute of
+ * Technology
+ *
+ * Permission is granted to copy or modify this code under the condition
+ * that the above copyright notice is maintained. The author and MIT
+ * make ABSOLUTELY NO WARRANTY regarding the fitness of this code for
+ * any purpose.
+ *
+ */
+
+/*---------------------------------------------------------------------------*
+ *
+ * Calculate near source traveltimes.
+ *
+ *
+ *
+ *---------------------------------------------------------------------------*/
+
+
+void near_source(float *ttime, float *slow, icoord *nm, icoord *isrc, fcoord *scale, int order)
+{
+ int ix_lo, ix_hi, iz_lo, iz_hi, ix, iz;
+ float slow_0, ttime_0, dist;
+
+ if(nm->y != 1)
+ verr("only 2D models implemented");
+
+/*---------------------------------------------------------------------------*
+ * Boundaries of source region.
+ *---------------------------------------------------------------------------*/
+
+ ix_lo = max2(isrc->x-order,0);
+ ix_hi = min2(isrc->x+order,nm->x-1);
+ iz_lo = max2(isrc->z-order,0);
+ iz_hi = min2(isrc->z+order,nm->z-1);
+
+/*---------------------------------------------------------------------------*
+ * Calculate traveltimes for source region.
+ * We assume this region to nearly homogeneous.
+ *---------------------------------------------------------------------------*/
+
+ slow_0 = slow[isrc->z*nm->x+isrc->x];
+ ttime_0 = ttime[isrc->z*nm->x+isrc->x];
+ for(iz=iz_lo; iz<=iz_hi; iz++) {
+ for(ix=ix_lo; ix<=ix_hi; ix++) {
+ dist = hypot((ix-isrc->x)*scale->x,(iz-isrc->z)*scale->z);
+ ttime[iz*nm->x+ix] = 0.5*dist*(slow[iz*nm->x+ix] + slow_0) + ttime_0;
+ }
+ }
+
+ return;
+}
+
+
+
+
+/* FILE: side.c
+ * AUTHOR: Joseph R. Matarese
+ * DATE: Decenber 18, 1992
+ * Copyright (c) 1992 Joseph R. Matarese and Massachusetts Institute of
+ * Technology
+ *
+ * Permission is granted to copy or modify this code under the condition
+ * that the above copyright notice is maintained. The author and MIT
+ * make ABSOLUTELY NO WARRANTY regarding the fitness of this code for
+ * any purpose.
+ *
+ */
+
+/*---------------------------------------------------------------------------*
+ *
+ * Apply the side stencil.
+ *
+ *---------------------------------------------------------------------------*/
+
+enum e_orientation { vertical, horizontal };
+
+void side(float *ttime, float *slow, icoord *nm, icoord *iin, icoord *iout, fcoord *scale, struct s_ecount *ecount)
+{
+ enum e_orientation orientation;
+ float dt, slow_0, operand, ttnew, tttmp;
+ static double c_sqrt2;
+ static char initialized = (char)0;
+
+ if(!initialized) {
+ if(nm->y != 1)
+ verr("only 2D models implemented");
+ if(scale->x != scale->z)
+ verr("only uniform grid implemented");
+ c_sqrt2 = sqrt(2.);
+ initialized = (char)1;
+ }
+
+ orientation = (iin->x == iout->x) ? horizontal : vertical;
+
+ if(orientation == vertical) {
+ dt = ttime[(iin->z+1)*nm->x+iin->x] - ttime[(iin->z-1)*nm->x+iin->x];
+ slow_0 = 0.25*(slow[(iin->z+1)*nm->x+iin->x] +
+ slow[(iin->z-1)*nm->x+iin->x] +
+ slow[iin->z*nm->x+iin->x] + slow[iin->z*nm->x+iout->x]);
+
+ if((operand = scale->x*scale->x*slow_0*slow_0 - 0.25*dt*dt) < 0.) {
+ slow_0 = 0.5*(slow[(iin->z-1)*nm->x+iin->x]+slow[iin->z*nm->x+iout->x]);
+ ttnew = ttime[(iin->z-1)*nm->x+iin->x] + c_sqrt2*scale->x*slow_0;
+
+ slow_0 = 0.5*(slow[(iin->z+1)*nm->x+iin->x]+slow[iin->z*nm->x+iout->x]);
+ tttmp = ttime[(iin->z+1)*nm->x+iin->x] + c_sqrt2*scale->x*slow_0;
+ if(tttmp < ttnew) ttnew = tttmp;
+
+ slow_0 = 0.5*(slow[iin->z*nm->x+iin->x] + slow[iin->z*nm->x+iout->x]);
+ tttmp = ttime[iin->z*nm->x+iin->x] + scale->x*slow_0;
+ if(tttmp < ttnew) ttnew = tttmp;
+
+ ecount->side++;
+ } else ttnew = ttime[iin->z*nm->x+iin->x] + sqrt(operand);
+ ttime[iin->z*nm->x+iout->x] = ttnew;
+
+ } else {
+ dt = ttime[iin->z*nm->x+iin->x+1] - ttime[iin->z*nm->x+iin->x-1];
+ slow_0 = 0.25*(slow[iin->z*nm->x+iin->x+1] + slow[iin->z*nm->x+iin->x-1] +
+ slow[iin->z*nm->x+iin->x] + slow[iout->z*nm->x+iin->x]);
+
+ if((operand = scale->x*scale->x*slow_0*slow_0 - 0.25*dt*dt) < 0.) {
+ slow_0 = 0.5*(slow[iin->z*nm->x+iin->x-1] + slow[iout->z*nm->x+iin->x]);
+ ttnew = ttime[iin->z*nm->x+iin->x-1] + c_sqrt2*scale->x*slow_0;
+
+ slow_0 = 0.5*(slow[iin->z*nm->x+iin->x+1] + slow[iout->z*nm->x+iin->x]);
+ tttmp = ttime[iin->z*nm->x+iin->x+1] + c_sqrt2*scale->x*slow_0;
+ if(tttmp < ttnew) ttnew = tttmp;
+
+ slow_0 = 0.5*(slow[iin->z*nm->x+iin->x] + slow[iout->z*nm->x+iin->x]);
+ tttmp = ttime[iin->z*nm->x+iin->x] + scale->x*slow_0;
+ if(tttmp < ttnew) ttnew = tttmp;
+
+ ecount->side++;
+ } else ttnew = ttime[iin->z*nm->x+iin->x] + sqrt(operand);
+ ttime[iout->z*nm->x+iin->x] = ttnew;
+
+ }
+
+ return;
+}
+
+
+
+
+/* FILE: corner.c
+ * AUTHOR: Joseph R. Matarese
+ * DATE: Decenber 18, 1992
+ * Copyright (c) 1992 Joseph R. Matarese and Massachusetts Institute of
+ * Technology
+ *
+ * Permission is granted to copy or modify this code under the condition
+ * that the above copyright notice is maintained. The author and MIT
+ * make ABSOLUTELY NO WARRANTY regarding the fitness of this code for
+ * any purpose.
+ *
+ */
+
+/*---------------------------------------------------------------------------*
+ *
+ * Apply the corner stencil.
+ *
+ *---------------------------------------------------------------------------*/
+
+
+void corner(float *ttime, float *slow,icoord *nm, icoord *iin, icoord *iout, fcoord *scale, struct s_ecount *ecount)
+{
+ float dt, slow_0, operand, ttnew, tttmp;
+ static double c_sqrt2;
+ static char initialized = (char)0;
+
+ if(!initialized) {
+ if(nm->y != 1)
+ verr("only 2D models implemented");
+ if(scale->x != scale->z)
+ verr("only uniform grid implemented");
+ c_sqrt2 = sqrt(2.);
+ initialized = (char)1;
+ }
+
+ dt = ttime[iout->z*nm->x+iin->x] - ttime[iin->z*nm->x+iout->x];
+ slow_0 = 0.25*(slow[iout->z*nm->x+iin->x] + slow[iin->z*nm->x+iout->x] +
+ slow[iout->z*nm->x+iout->x] + slow[iin->z*nm->x+iin->x]);
+
+ if((operand = 2*scale->x*scale->x*slow_0*slow_0 - dt*dt) < 0.) {
+ slow_0 = 0.5*(slow[iin->z*nm->x+iout->x] + slow[iout->z*nm->x+iout->x]);
+ ttnew = ttime[iin->z*nm->x+iout->x] + scale->x*slow_0;
+
+ slow_0 = 0.5*(slow[iout->z*nm->x+iin->x] + slow[iout->z*nm->x+iout->x]);
+ tttmp = ttime[iout->z*nm->x+iin->x] + scale->x*slow_0;
+ if(tttmp < ttnew) ttnew = tttmp;
+
+ slow_0 = 0.5*(slow[iout->z*nm->x+iout->x] + slow[iin->z*nm->x+iin->x]);
+ tttmp = ttime[iin->z*nm->x+iin->x] + c_sqrt2*scale->x*slow_0;
+ if(tttmp < ttnew) ttnew = tttmp;
+
+ ecount->corner++;
+ } else ttnew = ttime[iin->z*nm->x+iin->x] + sqrt(operand);
+
+/*---------------------------------------------------------------------------*
+ * Nonuniform grid?
+ *
+ scale2->x = scale->x * scale->x;
+ scale2->z = scale->z * scale->z;
+ scale2_plus = scale2->z + scale2->x;
+ scale2_minus = scale2->z + scale2->x;
+
+ ttnew = ttime[iin->z*nm->x+iin->x] +
+ (scale2_minus * dt + 2 * scale->z * scale->x *
+ sqrt(scale2_plus * slow_0 * slow_0 - dt * dt)) / scale2_plus;
+ *---------------------------------------------------------------------------*/
+
+ if(ttnew < ttime[iout->z*nm->x+iout->x])
+ ttime[iout->z*nm->x+iout->x] = ttnew;
+
+ return;
+}
+
+
+void rm_head(float *slow, icoord *ndim, icoord *isrc, int mzrcv, float dx, int *nzm)
+{
+ int iz, ix, k, l, i, izn, lprev, nz, nx;
+ float brd, val, dz;
+
+ iz = isrc->z;
+ ix = isrc->x;
+ nz = ndim->z;
+ nx = ndim->x;
+ dz = dx;
+
+ if (iz == 0) { ndim->z = 2; return;}
+
+ if (iz >= mzrcv) {
+ brd = slow[iz*nx+ix];
+ val = slow[(iz-1)*nx+ix];
+ while ((brd == val) && iz < nz-1) brd = slow[++iz*nx+ix];
+
+ lprev = iz;
+ for (l = iz; l > 0; l--) {
+ for (k = 0; k < nx; k++) {
+ if (slow[l*nx+k] == brd) {
+ slow[l*nx+k] = val;
+ lprev = l;
+ }
+ }
+ if (lprev-l) l = 0;
+ }
+
+ iz = isrc->z;
+ }
+ else {
+ brd = slow[iz*nx+ix];
+ val = slow[(iz+1)*nx+ix];
+ while (brd == val && iz > 0) brd = slow[--iz*nx+ix];
+
+ if (iz < nz) {
+ lprev = iz;
+ for (l = iz; l < nz; l++) {
+ for (k = 0; k < nx; k++) {
+ if (slow[l*nx+k] == brd) {
+ slow[l*nx+k] = val;
+ lprev = l;
+ }
+ }
+ if (lprev-l) l = nz;
+ }
+ }
+ iz = mzrcv;
+ }
+
+ *nzm = iz+1;
+
+ return;
+}
+
--
GitLab