diff --git a/index.html b/index.html
new file mode 100644
index 0000000000000000000000000000000000000000..7f1b66e8c5edeea520527b9380e5659ad3fa8231
--- /dev/null
+++ b/index.html
@@ -0,0 +1,18 @@
+
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
+<HTML>
+<HEAD>
+
+<!-- Global site tag (gtag.js) - Google Analytics -->
+<script async src="https://www.googletagmanager.com/gtag/js?id=UA-133845553-1"></script>
+<script>
+ window.dataLayer = window.dataLayer || [];
+ function gtag(){dataLayer.push(arguments);}
+ gtag('js', new Date());
+
+ gtag('config', 'UA-133845553-1');
+</script>
+
+</HEAD>
+<BODY BGCOLOR="#ffffff">
+<P>
diff --git a/marchenko/fmute.c b/marchenko/fmute.c
index 4d63faf38fa8bb4cb489ee554bca07e52ce9a439..f5c7668a16dc382b586c2981427888c6fa256f66 100644
--- a/marchenko/fmute.c
+++ b/marchenko/fmute.c
@@ -228,6 +228,8 @@ int main (int argc, char **argv)
/* alternative find maximum at source position */
dxrcv = (hdrs_in1[nx1-1].gx - hdrs_in1[0].gx)*sclsxgx/(float)(nx1-1);
imax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
+ /* make sure that the position fits into the receiver array */
+ imax = MIN(MAX(0,imax),nx1-1);
tmax=0.0;
jmax = 0;
xmax=0.0;
diff --git a/raytime/Makefile b/raytime/Makefile
index 192a17315ed0dc71f494646357315a4ffcc2dc82..b3aa6f39b0b3a26f00e0c892e352eebbbe223768 100644
--- a/raytime/Makefile
+++ b/raytime/Makefile
@@ -47,16 +47,16 @@ SRCC = $(PRG).c \
OBJC = $(SRCC:%.c=%.o)
$(PRG): $(OBJC) raytime.h
- $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o Raytime $(OBJC) $(LIBS)
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o raytime $(OBJC) $(LIBS)
install: raytime
- cp Raytime $B
+ cp raytime $B
clean:
- rm -f core $(OBJC) $(OBJM) Raytime
+ rm -f core $(OBJC) $(OBJM) raytime
realclean:
- rm -f core $(OBJC) $(OBJM) $(PRG) $B/Raytime
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/raytime
print: Makefile $(SRC)
diff --git a/raytime/raytime.c.veryold b/raytime/raytime.c.veryold
deleted file mode 100644
index e5e01b50b514111a8d45a27deb3421131244f40c..0000000000000000000000000000000000000000
--- a/raytime/raytime.c.veryold
+++ /dev/null
@@ -1,704 +0,0 @@
-#include <DELPHI_IOc.h>
-#include <errno.h>
-#include <memory.h>
-#include <time.h>
-#include <raytime.h>
-
-/* Plane-wave modeling externals */
-void plane_wave(float *tfinal,float *slowness, char *candidate, struct i_xyz *ndim, float scale, int order);
-
-/* Vidale modeling externals */
-extern void vidale(float *ttime, float *slow, struct i_xyz *nm, struct i_xyz *isrc, struct f_xyz *scale, int order);
-
-/* Graph_Theory modeling externals */
-extern void dijkstra(float *tfinal, float *slowness, char *candidate, int *raypath, XYZIndex *ndim, int order, int big, struct s_stencil *template);
-
-extern struct s_stencil *make_stencil(int order, XYZIndex *nm, XYZPosition *scale);
-
-extern void getrecpos(int *xi, int *zi, int *nrec, int nx, float h, float ox, float oz, int verbose);
-
-extern float setzsrc(int nb, int *boundary, float **inter, float *slow, int ni, float zsrc1, float dzsrc, float dz, float oz, int nx, int nz, float xsrc, float dx, float ox, int id, int verbose);
-
-void rm_head(float *slow, struct i_xyz *ndim, struct i_xyz *isrc, int mzrcv, struct f_xyz *scale, float **inter, int ni, int *nzm);
-
-extern void draw_rays(char *file_eps, int *raypath, struct i_xyz *ircv, int nx, int nz, float *slow, float dx, float dz, int nrec, float **inter, int ni, int hom, int verbose);
-
-extern void opint(float **data, int nrec, int Ns, int ix, float **dataT, float Tmin, float Tmax, float dT, int nT);
-
-/*********************** self documentation **********************/
-char *sdoc[] = {
-" ",
-" raytime - modeling of one-way traveltime for CFP operators",
-" ",
-" raytime file_vel= xsrc1= zsrc1= [optional parameters]",
-" ",
-" Required parameters:",
-" ",
-" file_vel= ................ gridded velocity file (DELPHI format)",
-" file_svel= ............... gridded velocity file (DELPHI format)",
-" xsrc1= ................... x-position of the source (m)",
-" zsrc1= ................... z-position of the source (m)",
-" ",
-" Optional parameters:",
-" ",
-" INPUT AND OUTPUT ",
-" file_out= ................ output file with traveltimes",
-" file_amp= ................ output file with approximate amplitudes",
-" file_int= ................ input file describing the interfaces (makemod)",
-" file_ray= ................ postscript file with rays (only method=graph)",
-" onegath=0 ................ 1; writes operators in one gather",
-" RAY TRACING ",
-" method=fd ................ calculation method (fd, plane or graph) ",
-" sbox=1 ................... radius of inner straight ray (fd method)",
-" order= ................... accuracy plane(=2)[0-2] and graph(=8)[1-10]",
-" radius=0 ................. radius in plane method",
-" dT=0 ..................... put traces on one-way time grid with step dT",
-" Tmin=0 ................... minimum time of one-way time grid (0 not used)",
-" Tmax=Tmin ................ maximum time of one-way time grid",
-" hom=1 .................... 1: draw straight rays in homogeneous layers",
-" SOURCE POSITIONS ",
-" xsrc2=xsrc1 .............. x-position of last source",
-" dxsrc=0 .................. step in source x-direction",
-" zsrc2=zsrc1 .............. z-position of last source",
-" dzsrc=0 .................. step in source z-direction",
-" boundary=0 ............... boundary to place the sources (overrules zsrc)",
-" add=0 .................... 1: adds all defined sources",
-" head=0 ................... 1: calculate head waves near source position",
-" RECEIVER POSITIONS ",
-" xrcv=0,(nx-1)*dx ......... x-position's of receivers (array)",
-" zrcv=0,0 ................. z-position's of receivers (array)",
-" Rboundary=0 .............. boundary to place the receivers (overrules zrcv)",
-" dxrcv=dx ................. step in receiver x-direction",
-" dzrcv=0 .................. step in receiver z-direction",
-" dxspr=0 .................. step of receiver spread in x-direction",
-" lint=1 ................... linear interpolate between the rcv points",
-" verbose=0 ................ verbose option",
-" ",
-" raytime calculates the first arrival time at the defined receiver array ",
-" for the defined shots at different depth and lateral positions.",
-" Note that one output gather consists of shots which are defined at the",
-" same lateral position but with different depth positions.",
-" Every new lateral position (with dxsrc) gives a new output gather.",
-" The parameter dT defines the one-way time between two shot records. The",
-" defined shot records (at depth) are used to get the interpolated values ",
-" at time steps of dT. ",
-" ",
-" PROGRAM TO CALCULATE TRAVEL TIMES IN 2D MEDIA ",
-" AUTHOR: Joseph R. Matarese ",
-" Copyright (c) 1993: Joseph R. Matarese and ",
-" Massachusetts Institute of Technology ",
-" ",
-" Translated to DELPHI environment: Jan Thorbecke 16-02-1996",
-" ",
-NULL};
-/**************** end self doc ***********************************/
-
-int main(int argc, char *argv[])
-{
- float32 *slowness, *slowness_S;
- uint8 *candidate;
- float32 *ttime, *ttime_p, slow_src;
- struct i_xyz *ndim, *isrc, *ircv;
- struct f_xyz scale;
- int id, nd, iz, ix, is, ir, ie, i;
- int iz_min, iz_max, ix_min, ix_max;
- int node_src, idz, idx, idp, sign;
- int *raypath;
- struct s_stencil *stencil;
-
- intn seqnr[MAX_KEYS];
- int32 type, dom1, dom2;
- int error, n1, n2, ret, size, verbose, nkeys, a;
- int k, nx, nz, sbox, Ns, nrec, ni, add, hom, ib, Nb;
- int *xi, *zi, j, ispr, ik, nzm, head, nT, mzrcv;
- int Nd, nb, *boundary, onegath, order, radius, Rboundary;
- float xsrc1, xsrc2, dxsrc, zsrc1=0, zsrc2, dzsrc, sx, sz, sl;
- float d1, d2, f1, f2, *tmpdata, dx, dz, dxspr, **data;
- float xsrc, *zsrc, drcv, dxrcv, dzrcv, t0, t1, t2, x, z, r, signz;
- float ox, oz, **inter, dT, Tmin, Tmax, **dataT, *trueslow;
- char *file_vel, *file_out, *file_int, *file_amp, *file_svel;
- char *keys[MAX_KEYS], *method, *file_ray, file_base[256], *pf;
- segyhdr *hdrs, *hdrsT;
-
- t0 = cputime();
- initargs(argc, argv);
- requestdoc(1);
-
-/*---------------------------------------------------------------------------*
- * Read input parameters and query for any that are needed but missing.
- *---------------------------------------------------------------------------*/
-
- if(!getparint("verbose", &verbose)) verbose = 0;
- if(!getparstring("file_vel", &file_vel)) saerr("file_vel not defined");
- if(!getparstring("file_svel", &file_svel)) file_svel=NULL;
- if(!getparstring("file_out", &file_out)) file_out=NULL;
- if(!getparstring("file_int", &file_int)) file_int=NULL;
- if(!getparstring("file_ray", &file_ray)) file_ray=NULL;
- if(!getparstring("file_amp", &file_amp)) file_amp=NULL;
- if(!getparstring("method", &method)) method="fd";
- if(!getparfloat("xsrc1", &xsrc1)) saerr("xsrc1 not defined");
- if(!getparfloat("xsrc2", &xsrc2)) xsrc2=xsrc1;
- if(!getparfloat("dxsrc", &dxsrc)) dxsrc=0;
- if(!getparfloat("Tmin", &Tmin)) Tmin=0;
- if(!getparint("Rboundary", &Rboundary)) Rboundary=0;
- if (Rboundary) {
- if(file_int == NULL) saerr("file_int must be specified for Rboundary");
- }
- nb = countparval("boundary");
- if(nb == 0 && Tmin == 0) {
- if(!getparfloat("zsrc1", &zsrc1))
- saerr("zsrc1 and boundary not defined, one must be defined");
- }
- else if (Tmin == 0) {
- if(file_int == NULL) saerr("file_int must be specified for boundary");
- boundary = alloc1int(nb);
- getparint("boundary", boundary);
- if (verbose) samess("source definition on boundary");
- }
- if(!getparfloat("zsrc2", &zsrc2)) zsrc2=zsrc1;
- if(!getparfloat("dzsrc", &dzsrc)) dzsrc=0;
- if(!getparint("head", &head)) head = 0;
- if(!getparint("sbox", &sbox)) sbox = 1;
- if(!getparint("onegath", &onegath)) onegath = 0;
- if(!getparint("add", &add)) add = 0;
- if(!getparint("hom", &hom)) hom = 1;
-
- if(equal(method,"fd")) {
- if (verbose)
- samess("finite_difference (Vidale, 1988, BSSA V. 78 #6, p. 2062)");
- }
- else if(equal(method,"graph")) {
- if (verbose)
- samess("graph_theory (Moser, 1991, Geophysics V. 56 #1, p. 59)");
- }
- else if(equal(method,"plane")) {
- if (verbose)
- samess("plane_wave (Matarese, 1993, Ph.D. Thesis, MIT)");
- }
- else {
- samess("unknown method: %s",method);
- samess("Possible choices include:\n\n");
- samess("\tfd based on Vidale's finite difference method\n");
- samess("\t (Vidale, 1988, BSSA V. 78 #6, p. 2062)\n");
- samess("\tgraph based on Moser's graph theoretical method\n");
- samess("\t (Moser, 1991, Geophysics V. 56 #1, p. 59)\n");
- samess("\tplane based on Matarese's plane wave extrapolation\n");
- samess("\t (Matarese, 1993, Ph.D. Thesis, MIT)\n");
- return(0);
- }
-
- if (file_ray != NULL) {
- pf = strrchr(file_ray, '.');
- *pf = '\0';
- if(!equal(method,"graph")) {
- sawarn("If file_ray is defined then method=graph");
- sawarn("So, method is set to graph");
- method = "graph";
- }
- }
-
- if(equal(method,"plane")) {
- if(!getparint("order", &order)) order = 2;
- if(order < 0 || order > 2) {
- sawarn("order must be within [0-2]");
- sawarn("order set to 2");
- order = 2;
- }
- if(!getparint("radius", &radius)) radius = 0;
- }
- else if(equal(method,"graph")) {
- if(!getparint("order", &order)) order = 8;
- if(order < 1) saerr("order must be within [1-10]");
- if(order > 10) sawarn("Attempting an order > 10. Good luck!");
- }
-
- if(equal(method,"plane") && add) {
- if (radius) {
- sawarn("Plane wave method with non-zero radius requires one source.");
- samess("Continuing with radius = 0.");
- }
- radius = 0;
- }
-
- if(add && equal(method,"fd"))
- saerr("Finite difference methods don't support extended source.");
-
-/*---------------------------------------------------------------------------*
- * Input the slowness grid.
- * It's gotta be 2-D and the variable is named "slowness".
- *---------------------------------------------------------------------------*
- * Open velocity file
- *---------------------------------------------------------------------------*/
-
- error = open_file(file_vel, GUESS_TYPE, DELPHI_READ);
- if (error < 0 ) saerr("error in opening file %s", file_vel);
- error = get_dims(file_vel, &n1, &n2, &type);
- size = n1 * n2;
- tmpdata = alloc1float(size);
- hdrs = (segyhdr *) malloc(n2*sizeof(segyhdr));
- read_data(file_vel,tmpdata,size,&n1,&n2,&f1,&f2,&d1,&d2,&type,hdrs);
- get_axis(&dom1, &dom2);
- if (verbose) disp_info(file_vel,n1,n2,f1,f2,d1,d2,type);
- ret = close_file(file_vel);
- if (ret < 0) sawarn("err %d on closing input file",ret);
- slowness = alloc1float(n1*n2);
-
- if (dom2 == SA_AXIS_X) {
- nx = n2; nz = n1;
- dx = d2; dz = d1;
- ox = f2; oz = f1;
- /* look at the coordinates of gx is the f2-axis is not defined */
- if (dx < 1e-9) {
- sawarn("f2 and d2 axis not defined, use the gx values");
- if (hdrs[0].scalco < 0) sl = 1.0/fabs(hdrs[0].scalco);
- else if (hdrs[0].scalco == 0) sl = 1.0;
- else sl = hdrs[0].scalco;
- ox = hdrs[0].gx*sl;
- dx = (hdrs[1].gx-hdrs[0].gx)*sl;
- }
-
- if (verbose) samess("Input model is transposed");
- for(ix=0; ix<nx; ix++) {
- for(iz=0; iz<nz; iz++) slowness[iz*nx+ix] = 1.0/tmpdata[ix*nz+iz];
- }
- }
- else {
- nx = n1; nz = n2;
- dx = d1; dz = d2;
- ox = f1; oz = f2;
- for(iz=0; iz<nz; iz++) {
- for(ix=0; ix<nx; ix++) slowness[iz*nx+ix] = 1.0/tmpdata[iz*nx+ix];
- }
- }
- free1float(tmpdata);
- free(hdrs);
-
-/*---------------------------------------------------------------------------*
- * Open S-wave velocity file
- *---------------------------------------------------------------------------*/
- if (file_svel!= NULL) {
- error = open_file(file_svel, GUESS_TYPE, DELPHI_READ);
- if (error < 0 ) saerr("error in opening file %s", file_svel);
- error = get_dims(file_svel, &n1, &n2, &type);
- size = n1 * n2;
- tmpdata = alloc1float(size);
- hdrs = (segyhdr *) malloc(n2*sizeof(segyhdr));
- read_data(file_svel,tmpdata,size,&n1,&n2,&f1,&f2,&d1,&d2,&type,hdrs);
- get_axis(&dom1, &dom2);
- if (verbose) disp_info(file_svel,n1,n2,f1,f2,d1,d2,type);
- ret = close_file(file_svel);
- if (ret < 0) sawarn("err %d on closing input file",ret);
- slowness_S = alloc1float(n1*n2);
-
- if (dom2 == SA_AXIS_X) {
- if (n2 != nx) saerr("nx of file %s (%d) != nx of file %s (%d)", file_vel, nx, file_svel, n2);
- if (n1 != nz) saerr("nz of file %s (%d) != nz of file %s (%d)", file_vel, nz, file_svel, n1);
-
- if (verbose) samess("Input S-model is transposed");
- for(ix=0; ix<nx; ix++) {
- for(iz=0; iz<nz; iz++) slowness_S[iz*nx+ix] = 1.0/tmpdata[ix*nz+iz];
- }
- }
- else {
- if (n1 != nx) saerr("nx of file %s (%d) != nx of file %s (%d)", file_vel, nx, file_svel, n1);
- if (n2 != nz) saerr("nz of file %s (%d) != nz of file %s (%d)", file_vel, nz, file_svel, n2);
- for(iz=0; iz<nz; iz++) {
- for(ix=0; ix<nx; ix++) slowness_S[iz*nx+ix] = 1.0/tmpdata[iz*nx+ix];
- }
- }
- free1float(tmpdata);
- free(hdrs);
- }
- else {
- slowness_S = alloc1float(nx*nz);
- for(iz=0; iz<nz; iz++) {
- for(ix=0; ix<nx; ix++)
- slowness_S[iz*nx+ix] = slowness[iz*nx+ix];
- }
- }
-
- if (NINT(dx*1000) != NINT(dz*1000)) saerr("dx must be equal to dz");
-
-/*---------------------------------------------------------------------------*
- * Open interface file (if available)
- *---------------------------------------------------------------------------*/
-
- if (file_int != NULL) {
- error = open_file(file_int, GUESS_TYPE, DELPHI_READ);
- if (error < 0 ) saerr("error in opening file %s", file_int);
- error = get_dims(file_int, &n1, &n2, &type);
- size = n1 * n2;
- tmpdata = alloc1float(size);
- hdrs = (segyhdr *) malloc(n2*sizeof(segyhdr));
- read_data(file_int,tmpdata,size,&n1,&n2,&f1,&f2,&d1,&d2,&type,hdrs);
- if (verbose>=2) disp_info(file_int,n1,n2,f1,f2,d1,d2,type);
- ret = close_file(file_int);
- free(hdrs);
- if (ret < 0) sawarn("err %d on closing input file",ret);
- ni = n2;
- if (n1 != nx) saerr("n1 != nx; wrong interface file");
-
- inter = alloc2float(nx, ni);
- for(i=0; i<ni; i++) {
- for(j=0; j<nx; j++) inter[i][j] = tmpdata[i*nx+j];
- }
- free1float(tmpdata);
- }
- else ni = 0;
-
-/*================ Read in receiver positions ================*/
-
- zi = alloc1int(nx+nz);
- xi = alloc1int(nx+nz);
- if(Rboundary<=0) {
- getrecpos(xi, zi, &nrec, nx, dz, ox, oz, verbose);
- }
- else {
- if (verbose) samess("Placing receivers on boundary %d.",Rboundary);
- if (verbose>=3) samess("receiver positions are:");
- if(!getparfloat("dxrcv",&dxrcv)) dxrcv = dx;
- nrec = NINT((nx-1)*dx/dxrcv)+1;
- for (ir = 0; ir < nrec; ir++) {
- xi[ir] = NINT(ir*dxrcv/dx);
- zi[ir] = NINT(inter[Rboundary-1][xi[ir]]/dz);
- if (verbose>=3) fprintf(stderr,"x=%f z=%f\n",(ox+xi[ir]*dx),(oz+zi[ir]*dz));
- }
- }
- if(!getparfloat("dxspr",&dxspr)) dxspr= 0;
- if(verbose) samess("nrec = %d", nrec);
-
-/* ============ Check and set parameters =============== */
-
- ispr = NINT(dxspr/dx);
- if (NINT(ispr*dx) != NINT(dxspr))
- saerr("dxspr is not a multiple of dx; this is not allowed");
-
- mzrcv = 0;
- for (ir = 0; ir < nrec; ir++) mzrcv = MAX(zi[ir], mzrcv);
- if (mzrcv > (nz-1)) saerr("deepest receiver outside model");
-
- if (nb) {dzsrc = 0.0; zsrc1 = inter[boundary[0]-1][0]; Nd = nb;}
- else if (dzsrc == 0) Nd = 1;
- else if (dzsrc != 0) Nd = NINT((zsrc2 - zsrc1)/dzsrc) + 1;
- if (dxsrc == 0) Ns = 1;
- else if (dxsrc != 0) Ns = NINT((xsrc2 - xsrc1)/dxsrc) + 1;
-
- if ((zsrc1+(Nd-1)*dzsrc-oz) > nz*dz) {
- sawarn("Deepest source outside model; last shot(s) not calculated");
- Nd -= 1;
- while( (zsrc1+(Nd-1)*dzsrc-oz ) > nz*dz) Nd--;
- }
- if (xi[nrec-1]*dx + (Ns-1)*dxspr > nx*dx)
- saerr("Moving spread moves outside model");
-
- if(!getparfloat("Tmin", &Tmin)) Tmin=0;
- if(!getparfloat("Tmax", &Tmax)) Tmax=Tmin;
- if(!getparfloat("dT", &dT)) dT=0;
- if (NINT(1000*dT) != 0) Nd = NINT((Tmax - Tmin)/dT) + 1;
- else if (Tmin) Nd = 1;
-
- if (verbose) {
- samess("Number of shot positions to generate = %d", Ns);
- samess("For every shot postion %d depth positions", Nd);
- samess("orig of model (x, z) = %.2f, %.2f", ox, oz);
- }
-
- ndim = (struct i_xyz *)jm_alloc(1,sizeof(struct i_xyz),1);
- ndim->z = nz; ndim->y = 1; ndim->x = nx;
- scale.x = dx; scale.y = 0.; scale.z = dz;
- nd = nz*nx;
-
-/*---------------------------------------------------------------------------*
- * If not finite difference method, allocate traveltime mask array(candidate).
- *---------------------------------------------------------------------------
- * Allocate the traveltime asrray.
- * If graph method, allocate the raypath and stencil arrays.
- *---------------------------------------------------------------------------*/
-
- candidate = (uint8 *)NULL;
- ttime = (float32 *)jm_alloc(nd, sizeof(float32), 0);
-
- if(equal(method,"graph")) {
- candidate = (uint8 *)jm_alloc(nd, sizeof(uint8), 0);
- raypath = (int *)jm_alloc(nd, sizeof(int), 0);
- stencil = make_stencil(order, ndim, &scale);
- }
-
-/* ============ Initializations =============== */
-
-/*---------------------------------------------------------------------------*
- * Input the source locations.
- * and
- * Initialize the traveltime array. Place t=0 @ source position.
- *---------------------------------------------------------------------------*/
-
- if (add) {
- data = alloc2float(nrec, 1);
- isrc = (struct i_xyz *)jm_alloc(Ns*Nd+1,sizeof(struct i_xyz),1);
- zsrc = alloc1float(Nd);
- for(id=0, ttime_p=ttime; id<nd; id++, ttime_p++)
- *ttime_p = Infinity;
-
- ie = 0;
- for (is = 0; is < Ns; is++) {
- xsrc = xsrc1 + is*dxsrc - ox;
- for (id = 0; id < Nd; id++) {
- zsrc[id] = setzsrc(nb,boundary,inter,slowness_S,ni,zsrc1,dzsrc,
- dz,oz,nx,nz,xsrc,dx,ox,id,verbose);
-
- isrc[ie].x = NINT(xsrc/dx);
- isrc[ie].y = 0;
- isrc[ie].z = NINT(zsrc[id]/dz);
- node_src = isrc[ie].z*nx + isrc[ie].x;
- sx = isrc[ie].x*dx-xsrc;
- sz = isrc[ie].z*dz-zsrc[id];
- sign = -1;
- if (sz < 0) sign = 1;
- ttime[node_src] = sign*sqrt(sx*sx+sz*sz)*slowness[node_src];
- if((isrc[ie].x > nx-1) || (isrc[ie].x < 0) ||
- (isrc[ie].z > nz-1) || (isrc[ie].z < 0))
- { saerr("source %d out of bounds ix=%d iz=%d", ie, isrc[ie].x, isrc[ie].z); }
- ie++;
- }
- }
- Ns = 1;
- Nd = 1;
- }
- else {
- data = alloc2float(nrec, Nd);
- isrc = (struct i_xyz *)jm_alloc(2,sizeof(struct i_xyz),1);
- zsrc = alloc1float(Nd);
- }
-
-/* ============ Initializations (2) =============== */
-
- if(!getparfloat("dxrcv",&dxrcv)) dxrcv = dx;
- if(!getparfloat("dzrcv",&dzrcv)) dzrcv = 0;
- drcv = sqrt(dxrcv*dxrcv+dzrcv*dzrcv);
- keys[0] = (char *) malloc(MAX_KEY_LENGTH);
- nkeys = 1;
- keys[0] = SA_OPER;
- seqnr[0] = 1;
- type = SA_TYPE_REAL;
- dom1 = SA_AXIS_X;
- if (dT > 0) dom2 = SA_AXIS_TIME;
- else dom2 = SA_AXIS_Z;
- if (file_ray != NULL) {
- ircv = (struct i_xyz *)jm_alloc(nrec+1,sizeof(struct i_xyz),1);
- }
- if (head == 0) {
- trueslow = alloc1float(nd);
- for (k = 0; k < nd; k++) trueslow[k] = slowness[k];
- }
-
- if (file_amp != NULL) {
- ret = open_file(file_amp, GUESS_TYPE, DELPHI_CREATE);
- if (ret < 0 ) saerr("error in creating output file %s", file_amp);
- }
- ret = open_file(file_out, GUESS_TYPE, DELPHI_CREATE);
- if (ret < 0 ) saerr("error in creating output file %s", file_out);
-
-/*---------------------------------------------------------------------------*
- * Compute traveltimes and (if applicable) raypaths.
- *---------------------------------------------------------------------------*/
-
- for (is = 0; is < Ns; is++) {
- xsrc = xsrc1 + is*dxsrc - ox;
- if (verbose) samess("**** gather %d ****", is+1);
-
- for (id = 0, ib=0; id < Nd; id++) {
- if (nb) {
- if (inter[boundary[id]-1][NINT(xsrc/dx)] == 0) continue;
- }
- zsrc[ib] = setzsrc(nb,boundary,inter,slowness_S,ni,zsrc1,dzsrc,dz,
- oz,nx,nz,xsrc,dx,ox,id,verbose);
- if (verbose) samess("xsrc = %f zsrc = %f", xsrc+ox, zsrc[ib]+oz);
-
- if (!add) {
- for(i=0, ttime_p=ttime; i<nd; i++, ttime_p++)
- *ttime_p = Infinity;
-
- isrc[0].x = NINT(xsrc/dx);
- isrc[0].y = 0;
- isrc[0].z = NINT(zsrc[ib]/dz);
- if((isrc[0].x > nx-1) || (isrc[0].x < 0) ||
- (isrc[0].z > nz-1) || (isrc[0].z < 0)) {
- { saerr("source %d out of bounds ix=%d iz=%d", is, isrc[is].x, isrc[is].z); }
- }
- node_src = isrc[0].z*nx + isrc[0].x;
- sx = isrc[0].x*dx-xsrc;
- sz = isrc[0].z*dz-zsrc[ib];
- sign = -1;
- if (sz < 0) sign = 1;
- ttime[node_src] = sign*sqrt(sx*sx+sz*sz)*slowness[node_src];
-
- if(equal(method,"plane")) {
- iz_min = max2(0,isrc[0].z-radius);
- iz_max = min2(ndim->z-1,isrc[0].z+radius);
- ix_min = max2(0,isrc[0].x-radius);
- ix_max = min2(ndim->x-1,isrc[0].x+radius);
-
- node_src = isrc[0].z*ndim->x + isrc[0].x;
- slow_src = slowness[node_src];
- for(iz=iz_min;iz<=iz_max;iz++) {
- idz = iz - isrc[0].z;
- for(ix=ix_min;ix<=ix_max;ix++) {
- idx = ix - isrc[0].x;
- idp = iz*ndim->x + ix;
- ttime[idp] = 0.5 * scale.x * sqrt(1.*idx*idx +
- idz*idz) * (slow_src + slowness[idp]);
- }
- }
- }
-
- /*=== avoid calculation of head waves from below zsrc ===*/
- if (head == 0) {
- ndim->z = nz;
- for (k = 0; k < nd; k++) slowness[k] = trueslow[k];
- rm_head(slowness,ndim,isrc,mzrcv,&scale,inter,ni,&nzm);
- ndim->z = MAX(nzm,2);
- }
-
- }
-
- t1 = cputime();
-
- if(equal(method,"plane")) {
- plane_wave(ttime,slowness,(char *)candidate,ndim,
- scale.x,order);
- }
- else if(equal(method,"fd")) {
- vidale(ttime,slowness,ndim,isrc,&scale,sbox);
- }
- else if(equal(method,"graph")) {
- dijkstra(ttime,slowness,(char *)candidate,raypath,ndim,
- order,0,stencil);
- }
-
- t2 = cputime();
- if (verbose>=3)
- samess("CPU-time computing traveltimes = %.2f s",t2-t1);
-
- for (ir = 0; ir < nrec; ir++) {
- ik = xi[ir] + is*ispr;
- data[ib][ir] = ttime[zi[ir]*nx+ik];
- }
-
- if(equal(method,"graph") && (file_ray != NULL)) {
- sprintf(file_base, "%s_s%02dd%02d.eps", file_ray, is+1, id+1);
-
- for(ir=0;ir<nrec;ir++) {
- ircv[ir].x = xi[ir];
- ircv[ir].y = 0;
- ircv[ir].z = zi[ir];
- }
- draw_rays(file_base,raypath,ircv,nx,nz,slowness,dx,dz,nrec,
- inter,ni,hom,verbose);
- if(verbose>=2)
- samess("finished plotting raypaths in postscript file %s", file_base);
- }
- ib++;
- }
- Nb = ib;
-
-/* ================ write to output file ================*/
-
- hdrs = (segyhdr *) malloc(Nd*sizeof(segyhdr));
- f1 = xi[0]*dx + ox;
- f2 = zsrc1;
- d2 = dxsrc;
- if (nb) f2 = 1.0;
- if (onegath) seqnr[0] = 1;
- else seqnr[0] = is+1;
-
- gen_hdrs(hdrs,nrec,Nd,f1,f2,drcv,dzsrc,TRID_ZX);
- for (i = 0; i < Nb; i++) {
- hdrs[i].scalco = -1000;
- hdrs[i].scalel = -1000;
- hdrs[i].offset = NINT(xi[0]*dx + is*ispr*dx - xsrc);
- hdrs[i].sx = NINT((xsrc+ox)*1000.0);
- hdrs[i].sdepth = NINT((zsrc[i]+oz)*1000.0);
- if (onegath) {
- hdrs[i].fldr = 1;
- hdrs[i].trwf = Ns*Nd;
- }
- else {
- hdrs[i].fldr = is+1;
- hdrs[i].trwf = Nd;
- }
- }
-
- ret = set_keys(keys, seqnr, nkeys);
- if (ret < 0 ) sawarn("error on writing keys.");
- ret = set_axis(dom1, dom2);
- if (ret < 0 ) saerr("error on writing axis.");
-
- if (verbose>1) disp_info(file_out,nrec,Nb,f1,f2,drcv,dzsrc,type);
- ret = write_data(file_out,*data,nrec,Nb,f1,f2,
- drcv,dzsrc,type,hdrs);
- if (ret < 0 ) saerr("error on writing output file.");
-
- free(hdrs);
-
- if (file_amp != NULL) {
- hdrs = (segyhdr *) malloc(Nd*sizeof(segyhdr));
- f1 = xi[0]*dx + ox;
- f2 = zsrc1;
- d2 = dxsrc;
- if (nb) f2 = 1.0;
- if (onegath) seqnr[0] = 1;
- else seqnr[0] = is+1;
-
- gen_hdrs(hdrs,nrec,Nd,f1,f2,drcv,dzsrc,TRID_ZX);
- for (i = 0; i < Nb; i++) {
- hdrs[i].scalco = -1000;
- hdrs[i].scalel = -1000;
- hdrs[i].offset = NINT(xi[0]*dx + is*ispr*dx - xsrc);
- hdrs[i].sx = NINT((xsrc+ox)*1000.0);
- hdrs[i].sdepth = NINT((zsrc[i]+oz)*1000.0);
- if (onegath) {
- hdrs[i].fldr = 1;
- hdrs[i].trwf = Ns*Nd;
- }
- else {
- hdrs[i].fldr = is+1;
- hdrs[i].trwf = Nd;
- }
- for (ir = 0; ir < nrec; ir++) {
- x = xsrc - (xi[ir] + is*ispr)*dx;
- z = zsrc[i] - (zi[ir] + is*ispr)*dzrcv;
- r = sqrt(x*x+z*z); /* cos(phi) = z/r */
- if (r != 0) data[i][ir] = fabs(z)/(r*sqrt(r));
- else data[i][ir] = 1.0;
- }
-
- }
-
- ret = set_keys(keys, seqnr, nkeys);
- if (ret < 0 ) sawarn("error on writing keys.");
- ret = set_axis(dom1, dom2);
- if (ret < 0 ) saerr("error on writing axis.");
- ret = write_data(file_amp,*data,nrec,Nb,f1,f2,
- drcv,dzsrc,type,hdrs);
- if (ret < 0 ) saerr("error on writing output file.");
-
- free(hdrs);
- }
- }
-
-/*---------------------------------------------------------------------------*
- * Output traveltime array to a file.
- *---------------------------------------------------------------------------*/
-
- ret = close_file(file_out);
- if (ret < 0) saerr("err %d on closing output file",ret);
- if (file_amp != NULL) {
- ret = close_file(file_amp);
- if (ret < 0) saerr("err %d on closing output file",ret);
- }
-
- t1 = cputime();
- if (verbose) samess("Total CPU-time = %f",t1-t0);
-
- return(0);
-}
diff --git a/raytime3d/Makefile b/raytime3d/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..be0bb6e043b09ee22c89254f9c54afceb990ad92
--- /dev/null
+++ b/raytime3d/Makefile
@@ -0,0 +1,72 @@
+# Makefile
+
+include ../Make_include
+
+########################################################################
+# define general include and system library
+ALLINC = -I.
+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 \
+# Grid2Time1.c \
+
+all: raytime3d
+
+PRG = raytime3d
+
+SRCC = $(PRG).c \
+ vidale3d.c \
+ src3d.c \
+ getParameters3d.c \
+ getWaveletInfo.c \
+ writeSrcRecPos.c \
+ readModel3d.c \
+ getWaveletHeaders.c \
+ verbosepkg.c \
+ getModelInfo3d.c \
+ wallclock_time.c \
+ recvPar3d.c \
+ writesufile.c \
+ name_ext.c \
+ atopkge.c \
+ docpkge.c \
+ threadAffinity.c \
+ getpars.c
+
+OBJC = $(SRCC:%.c=%.o)
+
+$(PRG): $(OBJC) raytime3d.h
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o raytime3d $(OBJC) $(LIBS)
+
+install: raytime3d
+ cp raytime3d $B
+
+clean:
+ rm -f core $(OBJC) $(OBJM) raytime3d
+
+realclean:
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/raytime3d
+
+
+print: Makefile $(SRC)
+ $(PRINT) $?
+ @touch print
+
+count:
+ @wc $(SRC)
+
+tar:
+ @tar cf $(PRG).tar Makefile $(SRC) && compress $(PRG).tar
+
+
+
diff --git a/raytime3d/SUsegy.h b/raytime3d/SUsegy.h
new file mode 100644
index 0000000000000000000000000000000000000000..a9133b999320911b29505909c2e4cd5f33b83dc5
--- /dev/null
+++ b/raytime3d/SUsegy.h
@@ -0,0 +1,391 @@
+/* This file is property of the Colorado School of Mines.
+
+ Copyright © 2007, Colorado School of Mines,
+ All rights reserved.
+
+
+ Redistribution and use in source and binary forms, with or
+ without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the following
+ disclaimer in the documentation and/or other materials provided
+ with the distribution.
+ * Neither the name of the Colorado School of Mines nor the names of
+ its contributors may be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ Warranty Disclaimer:
+ THIS SOFTWARE IS PROVIDED BY THE COLORADO SCHOOL OF MINES AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ COLORADO SCHOOL OF MINES OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+ STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+ IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+
+
+ Export Restriction Disclaimer:
+ We believe that CWP/SU: Seismic Un*x is a low technology product that does
+ not appear on the Department of Commerce CCL list of restricted exports.
+ Accordingly, we believe that our product meets the qualifications of
+ an ECCN (export control classification number) of EAR99 and we believe
+ it fits the qualifications of NRR (no restrictions required), and
+ is thus not subject to export restrictions of any variety.
+
+ Approved Reference Format:
+ In publications, please refer to SU as per the following example:
+ Cohen, J. K. and Stockwell, Jr. J. W., (200_), CWP/SU: Seismic Un*x
+ Release No. __: an open source software package for seismic
+ research and processing,
+ Center for Wave Phenomena, Colorado School of Mines.
+
+ Articles about SU in peer-reviewed journals:
+ Saeki, T., (1999), A guide to Seismic Un*x (SU)(2)---examples of data processing (part 1), data input and preparation of headers, Butsuri-Tansa (Geophysical Exploration), vol. 52, no. 5, 465-477.
+ Stockwell, Jr. J. W. (1999), The CWP/SU: Seismic Un*x Package, Computers and Geosciences, May 1999.
+ Stockwell, Jr. J. W. (1997), Free Software in Education: A case study of CWP/SU: Seismic Un*x, The Leading Edge, July 1997.
+ Templeton, M. E., Gough, C.A., (1998), Web Seismic Un*x: Making seismic reflection processing more accessible, Computers and Geosciences.
+
+ Acknowledgements:
+ SU stands for CWP/SU:Seismic Un*x, a processing line developed at Colorado
+ School of Mines, partially based on Stanford Exploration Project (SEP)
+ software.
+ */
+
+/* segy.h - include file for SEGY traces
+ *
+ * declarations for:
+ * typedef struct {} segy - the trace identification header
+ * typedef struct {} bhed - binary header
+ *
+ * Note:
+ * If header words are added, run the makefile in this directory
+ * to recreate hdr.h.
+ *
+ * Reference:
+ * K. M. Barry, D. A. Cavers and C. W. Kneale, "Special Report:
+ * Recommended Standards for Digital Tape Formats",
+ * Geophysics, vol. 40, no. 2 (April 1975), P. 344-352.
+ *
+ * $Author: john $
+ * $Source: /usr/local/cwp/src/su/include/RCS/segy.h,v $
+ * $Revision: 1.23 $ ; $Date: 1998/03/26 23:48:18 $
+ */
+
+/* #define SU_NFLTS 800000 Arbitrary limit on data array size */
+
+/**
+* This segyhdr has been redefined and uses an integer (32 bit) for number of samples (ns)
+* Jan Thorbecke
+**/
+
+
+/* TYPEDEFS */
+typedef struct { /* segy - trace identification header */
+
+ int tracl; /* trace sequence number within line */
+
+ int tracr; /* trace sequence number within reel */
+
+ int fldr; /* field record number */
+
+ int tracf; /* trace number within field record */
+
+ int ep; /* energy source point number */
+
+ int cdp; /* CDP ensemble number */
+
+ int cdpt; /* trace number within CDP ensemble */
+
+ short trid; /* trace identification code:
+ 1 = seismic data
+ 2 = dead
+ 3 = dummy
+ 4 = time break
+ 5 = uphole
+ 6 = sweep
+ 7 = timing
+ 8 = water break
+ 9---, N = optional use (N = 32,767)
+
+ Following are CWP id flags:
+
+ 9 = autocorrelation
+
+ 10 = Fourier transformed - no packing
+ xr[0],xi[0], ..., xr[N-1],xi[N-1]
+
+ 11 = Fourier transformed - unpacked Nyquist
+ xr[0],xi[0],...,xr[N/2],xi[N/2]
+
+ 12 = Fourier transformed - packed Nyquist
+ even N:
+ xr[0],xr[N/2],xr[1],xi[1], ...,
+ xr[N/2 -1],xi[N/2 -1]
+ (note the exceptional second entry)
+ odd N:
+ xr[0],xr[(N-1)/2],xr[1],xi[1], ...,
+ xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2]
+ (note the exceptional second & last entries)
+
+ 13 = Complex signal in the time domain
+ xr[0],xi[0], ..., xr[N-1],xi[N-1]
+
+ 14 = Fourier transformed - amplitude/phase
+ a[0],p[0], ..., a[N-1],p[N-1]
+
+ 15 = Complex time signal - amplitude/phase
+ a[0],p[0], ..., a[N-1],p[N-1]
+
+ 16 = Real part of complex trace from 0 to Nyquist
+
+ 17 = Imag part of complex trace from 0 to Nyquist
+
+ 18 = Amplitude of complex trace from 0 to Nyquist
+
+ 19 = Phase of complex trace from 0 to Nyquist
+
+ 21 = Wavenumber time domain (k-t)
+
+ 22 = Wavenumber frequency (k-omega)
+
+ 23 = Envelope of the complex time trace
+
+ 24 = Phase of the complex time trace
+
+ 25 = Frequency of the complex time trace
+
+ 30 = Depth-Range (z-x) traces
+
+ 43 = Seismic Data, Vertical Component
+
+ 44 = Seismic Data, Horizontal Component 1
+
+ 45 = Seismic Data, Horizontal Component 2
+
+ 46 = Seismic Data, Radial Component
+
+ 47 = Seismic Data, Transverse Component
+
+ 101 = Seismic data packed to bytes (by supack1)
+
+ 102 = Seismic data packed to 2 bytes (by supack2)
+ */
+
+ short nvs; /* number of vertically summed traces (see vscode
+ in bhed structure) */
+
+ short nhs; /* number of horizontally summed traces (see vscode
+ in bhed structure) */
+
+ short duse; /* data use:
+ 1 = production
+ 2 = test */
+
+ int offset; /* distance from source point to receiver
+ group (negative if opposite to direction
+ in which the line was shot) */
+
+ int gelev; /* receiver group elevation from sea level
+ (above sea level is positive) */
+
+ int selev; /* source elevation from sea level
+ (above sea level is positive) */
+
+ int sdepth; /* source depth (positive) */
+
+ int gdel; /* datum elevation at receiver group */
+
+ int sdel; /* datum elevation at source */
+
+ int swdep; /* water depth at source */
+
+ int gwdep; /* water depth at receiver group */
+
+ short scalel; /* scale factor for previous 7 entries
+ with value plus or minus 10 to the
+ power 0, 1, 2, 3, or 4 (if positive,
+ multiply, if negative divide) */
+
+ short scalco; /* scale factor for next 4 entries
+ with value plus or minus 10 to the
+ power 0, 1, 2, 3, or 4 (if positive,
+ multiply, if negative divide) */
+
+ int sx; /* X source coordinate */
+
+ int sy; /* Y source coordinate */
+
+ int gx; /* X group coordinate */
+
+ int gy; /* Y group coordinate */
+
+ short counit; /* coordinate units code:
+ for previous four entries
+ 1 = length (meters or feet)
+ 2 = seconds of arc (in this case, the
+ X values are longitude and the Y values
+ are latitude, a positive value designates
+ the number of seconds east of Greenwich
+ or north of the equator */
+
+ short wevel; /* weathering velocity */
+
+ short swevel; /* subweathering velocity */
+
+ short sut; /* uphole time at source */
+
+ short gut; /* uphole time at receiver group */
+
+ short sstat; /* source static correction */
+
+ short gstat; /* group static correction */
+
+ short tstat; /* total static applied */
+
+ short laga; /* lag time A, time in ms between end of 240-
+ byte trace identification header and time
+ break, positive if time break occurs after
+ end of header, time break is defined as
+ the initiation pulse which maybe recorded
+ on an auxiliary trace or as otherwise
+ specified by the recording system */
+
+ short lagb; /* lag time B, time in ms between the time break
+ and the initiation time of the energy source,
+ may be positive or negative */
+
+ short delrt; /* delay recording time, time in ms between
+ initiation time of energy source and time
+ when recording of data samples begins
+ (for deep water work if recording does not
+ start at zero time) */
+
+ short muts; /* mute time--start */
+
+ short mute; /* mute time--end */
+
+ short igc; /* instrument gain constant */
+
+ int ns; /* number of samples in this trace */
+
+ unsigned short dt; /* sample interval; in micro-seconds */
+
+ short gain; /* gain type of field instruments code:
+ 1 = fixed
+ 2 = binary
+ 3 = floating point
+ 4 ---- N = optional use */
+
+ short igi; /* instrument early or initial gain */
+
+ short corr; /* correlated:
+ 1 = no
+ 2 = yes */
+
+ short sfs; /* sweep frequency at start */
+
+ short sfe; /* sweep frequency at end */
+
+ short slen; /* sweep length in ms */
+
+ short styp; /* sweep type code:
+ 1 = linear
+ 2 = cos-squared
+ 3 = other */
+
+ short stas; /* sweep trace length at start in ms */
+
+ short stae; /* sweep trace length at end in ms */
+
+ short tatyp; /* taper type: 1=linear, 2=cos^2, 3=other */
+
+ short afilf; /* alias filter frequency if used */
+
+ short afils; /* alias filter slope */
+
+ short nofilf; /* notch filter frequency if used */
+
+ short nofils; /* notch filter slope */
+
+ short lcf; /* low cut frequency if used */
+
+ short hcf; /* high cut frequncy if used */
+
+ short lcs; /* low cut slope */
+
+ short hcs; /* high cut slope */
+
+ short year; /* year data recorded */
+
+ short day; /* day of year */
+
+ short hour; /* hour of day (24 hour clock) */
+
+ short minute; /* minute of hour */
+
+ short sec; /* second of minute */
+
+ short timbas; /* time basis code:
+ 1 = local
+ 2 = GMT
+ 3 = other */
+
+ short trwf; /* trace weighting factor, defined as 1/2^N
+ volts for the least sigificant bit */
+
+ short grnors; /* geophone group number of roll switch
+ position one */
+
+ short grnofr; /* geophone group number of trace one within
+ original field record */
+
+ short grnlof; /* geophone group number of last trace within
+ original field record */
+
+ short gaps; /* gap size (total number of groups dropped) */
+
+ short otrav; /* overtravel taper code:
+ 1 = down (or behind)
+ 2 = up (or ahead) */
+
+ /* local assignments */
+
+ short mark; /* mark selected traces */
+
+ float d1; /* sample spacing for non-seismic data */
+
+ float f1; /* first sample location for non-seismic data */
+
+ float d2; /* sample spacing between traces */
+
+ float f2; /* first trace location */
+
+ float ungpow; /* negative of power used for dynamic
+ range compression */
+
+ float unscale; /* reciprocal of scaling factor to normalize
+ range */
+
+ int ntr; /* number of traces */
+
+/* short shortpad; alignment padding */
+
+ short unass[14]; /* unassigned--NOTE: last entry causes
+ a break in the word alignment, if we REALLY
+ want to maintain 240 bytes, the following
+ entry should be an odd number of short/UINT2
+ OR do the insertion above the "mark" keyword
+ entry */
+
+} SUsegy;
+
+
diff --git a/raytime3d/atopkge.c b/raytime3d/atopkge.c
new file mode 100644
index 0000000000000000000000000000000000000000..ef0b21803a52fb2ee021b6d73938d44bcefb0f92
--- /dev/null
+++ b/raytime3d/atopkge.c
@@ -0,0 +1,444 @@
+/*
+
+ This file is property of the Colorado School of Mines.
+
+ Copyright (C) 2007, Colorado School of Mines,
+ All rights reserved.
+
+
+ Redistribution and use in source and binary forms, with or
+ without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the following
+ disclaimer in the documentation and/or other materials provided
+ with the distribution.
+ * Neither the name of the Colorado School of Mines nor the names of
+ its contributors may be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ Warranty Disclaimer:
+ THIS SOFTWARE IS PROVIDED BY THE COLORADO SCHOOL OF MINES AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ COLORADO SCHOOL OF MINES OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+ STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+ IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+
+
+ Export Restriction Disclaimer:
+ We believe that CWP/SU: Seismic Un*x is a low technology product that does
+ not appear on the Department of Commerce CCL list of restricted exports.
+ Accordingly, we believe that our product meets the qualifications of
+ an ECCN (export control classification number) of EAR99 and we believe
+ it fits the qualifications of NRR (no restrictions required), and
+ is thus not subject to export restrictions of any variety.
+
+ Approved Reference Format:
+ In publications, please refer to SU as per the following example:
+ Cohen, J. K. and Stockwell, Jr. J. W., (200_), CWP/SU: Seismic Un*x
+ Release No. __: an open source software package for seismic
+ research and processing,
+ Center for Wave Phenomena, Colorado School of Mines.
+
+ Articles about SU in peer-reviewed journals:
+ Saeki, T., (1999), A guide to Seismic Un*x (SU)(2)---examples of data processing (part 1), data input and preparation of headers, Butsuri-Tansa (Geophysical Exploration), vol. 52, no. 5, 465-477.
+ Stockwell, Jr. J. W. (1999), The CWP/SU: Seismic Un*x Package, Computers and Geosciences, May 1999.
+ Stockwell, Jr. J. W. (1997), Free Software in Education: A case study of CWP/SU: Seismic Un*x, The Leading Edge, July 1997.
+ Templeton, M. E., Gough, C.A., (1998), Web Seismic Un*x: Making seismic reflection processing more accessible, Computers and Geosciences.
+
+ Acknowledgements:
+ SU stands for CWP/SU:Seismic Un*x, a processing line developed at Colorado
+ School of Mines, partially based on Stanford Exploration Project (SEP)
+ software.
+ */
+
+/*********************** self documentation **********************/
+/***************************************************************************
+ATOPKGE - convert ascii to arithmetic and with error checking
+
+
+eatoh ascii to short
+eatou ascii to unsigned short
+eatoi ascii to int
+eatop ascii to unsigned
+eatol ascii to long
+eatov ascii to unsigned long
+eatof ascii to float
+eatod ascii to double
+
+****************************************************************************
+Function Prototypes:
+short eatoh(char *s);
+unsigned short eatou(char *s);
+int eatoi(char *s);
+unsigned int eatop(char *s);
+long eatol(char *s);
+unsigned long eatov(char *s);
+float eatof(char *s);
+double eatod(char *s);
+
+****************************************************************************
+Input:
+s string
+
+Returned: type indicated
+
+****************************************************************************
+Notes:
+Each of these routines acts like atoi, but has error checking:
+
+This is a major revision of the tedious code required before
+vendors implemented the ANSI C strtol, strtoul and strtod.
+
+In addition to the size checks for each integer type, a
+specific test on errno is required. For example, INT_MAX
+may (and probably does) equal LONG_MAX. In this case,
+if fed a number exceeding INT_MAX (and LONG_MAX), strtol
+will make a quiet return with the wrong answer and it is up
+to the user to check if errno == ERANGE.
+
+Size limits are machine dependent and are read from the
+ANSI C include files limits.h and float.h.
+
+Bug Report: With NeXT c and Gnucc, when x > DBL_MAX (x <-DBL_MAX),
+the return value from strtod was +Infinity (-Infinity), not HUGE_VAL
+and more important, errno was not set to ERANGE. To cope with this,
+I put explicit size checks in eatod (which would not be needed if
+errno were set as it should be in ANSI C. jkc 01/29/94
+
+On IBM RS6000, the return value from strtod was +-Inf on
+overflow, but errno was set correctly.
+
+****************************************************************************
+References:
+For old code:
+Plum: Reliable Data Structures in C, p. 2-17.
+Kernighan and Ritchie: The C Programming Language, p. 58.
+
+CWP: Jack K. Cohen, Brian Sumner
+
+For new code:
+ANSI C routines with a little help from Jack
+
+****************************************************************************
+Author: Jack Cohen, Center for Wave Phenomena, 1994.
+***************************************************************************/
+/**************** end self doc ********************************/
+
+#include "par.h"
+#include <float.h>
+#include <limits.h>
+#include <stdarg.h>
+#include <errno.h>
+
+/* eatoh - convert string s to short integer {SHRT_MIN:SHRT_MAX} */
+short eatoh(char *s)
+{
+ long n = strtol(s, NULL, 10);
+
+ if ( (n > SHRT_MAX) || (n < SHRT_MIN) || (errno == ERANGE) )
+ err("%s: eatoh: overflow", __FILE__);
+
+ return (short) n;
+}
+
+
+/* eatou - convert string s to unsigned short integer {0:USHRT_MAX} */
+unsigned short eatou(char *s)
+{
+ unsigned long n = strtoul(s, NULL, 10);
+
+ if ( (n > USHRT_MAX) || (errno == ERANGE) )
+ err("%s: eatou: overflow", __FILE__);
+
+ return (unsigned short) n;
+}
+
+
+/* eatoi - convert string s to integer {INT_MIN:INT_MAX} */
+int eatoi(char *s)
+{
+ long n = strtol(s, NULL, 10);
+
+ if ( (n > INT_MAX) || (n < INT_MIN) || (errno == ERANGE) )
+ err("%s: eatoi: overflow", __FILE__);
+
+ return (int) n;
+}
+
+
+/* eatop - convert string s to unsigned integer {0:UINT_MAX} */
+unsigned int eatop(char *s)
+{
+ unsigned long n = strtoul(s, NULL, 10);
+
+ if ( (n > UINT_MAX) || (errno == ERANGE) )
+ err("%s: eatop: overflow", __FILE__);
+
+ return (unsigned int) n;
+}
+
+
+/* eatol - convert string s to long integer {LONG_MIN:LONG_MAX} */
+long eatol(char *s)
+{
+ long n = strtol(s, NULL, 10);
+
+ if (errno == ERANGE)
+ err("%s: eatol: overflow", __FILE__);
+
+ return n;
+}
+
+
+/* eatov - convert string s to unsigned long {0:ULONG_MAX} */
+unsigned long eatov(char *s)
+{
+ unsigned long n = strtoul(s, NULL, 10);
+
+ if (errno == ERANGE)
+ err("%s: eatov: overflow", __FILE__);
+
+ return n;
+}
+
+
+/* eatof - convert string s to float {-FLT_MAX:FLT_MAX} */
+float eatof(char *s)
+{
+ float x = strtod(s, NULL);
+
+ if ( (x > FLT_MAX) || (x < -FLT_MAX) || (errno == ERANGE) )
+ err("%s: eatof: overflow", __FILE__);
+
+ return (float) x;
+}
+
+
+/* eatod - convert string s to double {-DBL_MAX:DBL_MAX} */
+double eatod(char *s)
+{
+ double x = strtod(s, NULL);
+
+ /* errno == ERANGE suffices if compiler sets errno on overflow */
+ if ( (errno == ERANGE) || (x > DBL_MAX) || (x < -DBL_MAX) )
+ err("%s: eatod: overflow", __FILE__);
+
+ return x;
+}
+
+
+/**************************************************************************
+ERRPKGE - routines for reporting errors
+
+err print warning on application program error and die
+warn print warning on application program error
+syserr print warning on application program error using errno and die
+
+***************************************************************************
+Function Prototypes:
+void err (char *fmt, ...);
+void warn (char *fmt, ...);
+void syserr (char *fmt, ...);
+
+***************************************************************************
+Return: void
+
+***************************************************************************
+Notes:
+fmt a printf format string ("\n" not needed)
+... the variables referenced in the format string
+
+Examples:
+ err("Cannot divide %f by %f", x, y);
+ warn("fmax = %f exceeds half nyquist= %f", fmax, 0.25/dt);
+
+ if (NULL == (fp = fopen(xargv[1], "r")))
+ err("can't open %s", xargv[1]);
+ ...
+ if (-1 == close(fd))
+ err("close failed");
+
+***************************************************************************
+References:
+Kernighan and Pike, "The UNIX Programming Environment", page 207.
+Also Rochkind, "Advanced UNIX Programming", page 13.
+
+***************************************************************************
+Authors:SEP: Jeff Thorson, Stew Levin CWP: Shuki Ronen, Jack Cohen
+**************************************************************************/
+
+
+void err(char *fmt, ...)
+{
+ va_list args;
+
+
+ if (EOF == fflush(stdout)) {
+ fprintf(stderr, "\nerr: fflush failed on stdout");
+ }
+ fprintf(stderr, "\n%s: ", xargv[0]);
+ va_start(args,fmt);
+ vfprintf(stderr, fmt, args);
+ va_end(args);
+ fprintf(stderr, "\n");
+ exit(EXIT_FAILURE);
+}
+
+
+void warn(char *fmt, ...)
+{
+ va_list args;
+
+ if (EOF == fflush(stdout)) {
+ fprintf(stderr, "\nwarn: fflush failed on stdout");
+ }
+ fprintf(stderr, "\n%s: ", xargv[0]);
+ va_start(args,fmt);
+ vfprintf(stderr, fmt, args);
+ va_end(args);
+ fprintf(stderr, "\n");
+ return;
+}
+
+
+void syserr(char *fmt, ...)
+{
+ va_list args;
+
+ if (EOF == fflush(stdout)) {
+ fprintf(stderr, "\nsyserr: fflush failed on stdout");
+ }
+ fprintf(stderr, "\n%s: ", xargv[0]);
+ va_start(args,fmt);
+ vfprintf(stderr, fmt, args);
+ va_end(args);
+ fprintf(stderr, " (%s)\n", strerror(errno));
+ exit(EXIT_FAILURE);
+}
+
+#ifdef TEST
+main(int argc, char **argv)
+{
+ char s[BUFSIZ];
+ short nh;
+ unsigned short nu;
+ int ni;
+ unsigned int np;
+ long nl;
+ unsigned long nv;
+
+ initargs(argc, argv);
+
+
+ /* Test code for eatoh */
+ if (SHRT_MAX == LONG_MAX) {
+ warn("Warning: eatoh not used on this machine.\n");
+ } else {
+ warn("\n");
+ }
+ strcpy(s, "0");
+ nh = eatoh(s);
+ warn("eatoh(%s) = %hd\n", s, nh);
+
+ strcpy(s, "32767");
+ nh = eatoh(s);
+ warn("eatoh(%s) = %hd\n", s, nh);
+
+ strcpy(s, "-32768");
+ nh = eatoh(s);
+ warn("eatoh(%s) = %hd\n", s, nh);
+
+
+ /* Test code for eatou */
+ if (USHRT_MAX == ULONG_MAX) {
+ warn("Warning: eatou not used on this machine.\n");
+ } else {
+ warn("\n");
+ }
+ strcpy(s, "0");
+ nu = eatou(s);
+ warn("eatou(%s) = %hu\n", s, nu);
+
+ strcpy(s, "65535");
+ nu = eatou(s);
+ warn("eatou(%s) = %hu\n", s, nu);
+
+
+ /* Test code for eatoi */
+ if (INT_MAX == LONG_MAX) {
+ warn("Warning: eatoi not used on this machine.\n");
+ } else {
+ warn("\n");
+ }
+ strcpy(s, "0");
+ ni = eatoi(s);
+ warn("eatoi(%s) = %d\n", s, ni);
+
+ strcpy(s, "2147483647");
+ ni = eatoi(s);
+ warn("eatoi(%s) = %d\n", s, ni);
+
+
+ strcpy(s, "-2147483648");
+ ni = eatoi(s);
+ warn("eatoi(%s) = %d\n", s, ni);
+
+
+ /* Test code for eatop */
+ if (INT_MAX == LONG_MAX) {
+ warn("Warning: eatop not used on this machine.\n");
+ } else {
+ warn("\n");
+ }
+ strcpy(s, "0");
+ np = eatop(s);
+ warn("eatop(%s) = %lu\n", s, np);
+
+ strcpy(s, "4294967295");
+ np = eatop(s);
+ warn("eatop(%s) = %lu\n", s, np);
+
+
+ /* Test code for eatol */
+ warn("\n");
+ strcpy(s, "0");
+ nl = eatol(s);
+ warn("eatol(%s) = %ld\n", s, nl);
+
+ strcpy(s, "2147483647");
+ nl = eatol(s);
+ warn("eatol(%s) = %ld\n", s, nl);
+
+ strcpy(s, "-2147483648");
+ nl = eatol(s);
+ warn("eatol(%s) = %ld\n", s, nl);
+
+
+ /* Test code for eatov */
+ strcpy(s, "0");
+ nv = eatov(s);
+ warn("eatov(%s) = %lu\n", s, nv);
+
+ strcpy(s, "4294967295");
+ nv = eatov(s);
+ warn("eatov(%s) = %lu\n", s, nv);
+
+ warn("Now we feed in 4294967296, expecting fatal error exit\n");
+ strcpy(s, "4294967296");
+ nv = eatov(s);
+ warn("eatov(%s) = %lu\n", s, nv);
+
+ return EXIT_SUCCESS;
+}
+#endif
diff --git a/raytime3d/docpkge.c b/raytime3d/docpkge.c
new file mode 100644
index 0000000000000000000000000000000000000000..a74b4c331cbb3c882fd311ef5dacbda776c14e04
--- /dev/null
+++ b/raytime3d/docpkge.c
@@ -0,0 +1,188 @@
+/*
+ This file is property of the Colorado School of Mines.
+
+ Copyright (C) 2007, Colorado School of Mines,
+ All rights reserved.
+
+
+ Redistribution and use in source and binary forms, with or
+ without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the following
+ disclaimer in the documentation and/or other materials provided
+ with the distribution.
+ * Neither the name of the Colorado School of Mines nor the names of
+ its contributors may be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ Warranty Disclaimer:
+ THIS SOFTWARE IS PROVIDED BY THE COLORADO SCHOOL OF MINES AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ COLORADO SCHOOL OF MINES OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+ STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+ IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+
+
+ Export Restriction Disclaimer:
+ We believe that CWP/SU: Seismic Un*x is a low technology product that does
+ not appear on the Department of Commerce CCL list of restricted exports.
+ Accordingly, we believe that our product meets the qualifications of
+ an ECCN (export control classification number) of EAR99 and we believe
+ it fits the qualifications of NRR (no restrictions required), and
+ is thus not subject to export restrictions of any variety.
+
+ Approved Reference Format:
+ In publications, please refer to SU as per the following example:
+ Cohen, J. K. and Stockwell, Jr. J. W., (200_), CWP/SU: Seismic Un*x
+ Release No. __: an open source software package for seismic
+ research and processing,
+ Center for Wave Phenomena, Colorado School of Mines.
+
+ Articles about SU in peer-reviewed journals:
+ Saeki, T., (1999), A guide to Seismic Un*x (SU)(2)---examples of data processing (part 1), data input and preparation of headers, Butsuri-Tansa (Geophysical Exploration), vol. 52, no. 5, 465-477.
+ Stockwell, Jr. J. W. (1999), The CWP/SU: Seismic Un*x Package, Computers and Geosciences, May 1999.
+ Stockwell, Jr. J. W. (1997), Free Software in Education: A case study of CWP/SU: Seismic Un*x, The Leading Edge, July 1997.
+ Templeton, M. E., Gough, C.A., (1998), Web Seismic Un*x: Making seismic reflection processing more accessible, Computers and Geosciences.
+
+ Acknowledgements:
+ SU stands for CWP/SU:Seismic Un*x, a processing line developed at Colorado
+ School of Mines, partially based on Stanford Exploration Project (SEP)
+ software.
+ */
+
+/*********************** self documentation **********************/
+/***************************************************************************
+DOCPKGE - Function to implement the CWP self-documentation facility
+
+requestdoc give selfdoc on user request (i.e. when name of main is typed)
+pagedoc print self documentation string
+
+****************************************************************************
+Function Prototypes:
+void requestdoc(flag);
+void pagedoc();
+
+****************************************************************************
+requestoc:
+Input:
+flag integer specifying i.o. cases
+
+pagedoc():
+Returns: the self-documentation, an array of strings
+
+****************************************************************************
+Notes:
+requestdoc:
+In the usual case, stdin is used to pass in data. However,
+some programs (eg. synthetic data generators) don't use stdin
+to pass in data and some programs require two or more arguments
+besides the command itself (eg. sudiff) and don't use stdin.
+In this last case, we give selfdoc whenever too few arguments
+are given, since these usages violate the usual SU syntax.
+In all cases, selfdoc can be requested by giving only the
+program name.
+
+The flag argument distinguishes these cases:
+ flag = 0; fully defaulted, no stdin
+ flag = 1; usual case
+ flag = n > 1; no stdin and n extra args required
+
+pagedoc:
+Intended to be called by requesdoc(), but conceivably could be
+used directly as in:
+ if (xargc != 3) selfdoc();
+
+Based on earlier versions by:
+SEP: Einar Kjartansson, Stew Levin CWP: Jack Cohen, Shuki Ronen
+HRC: Lyle
+
+****************************************************************************
+Author: Jack K. Cohen, Center for Wave Phenomena
+****************************************************************************/
+/**************** end self doc ********************************/
+
+#include "par.h"
+
+#ifndef EXIT_FAILURE
+#define EXIT_FAILURE (1)
+#endif
+#ifndef EXIT_SUCCESS
+#define EXIT_SUCCESS (0)
+#endif
+
+
+/* definitions of global variables */
+int xargc; char **xargv;
+
+
+void requestdoc(int flag)
+/***************************************************************************
+print selfdocumentation as directed by the user-specified flag
+****************************************************************************
+Notes:
+In the usual case, stdin is used to pass in data. However,
+some programs (eg. synthetic data generators) don't use stdin
+to pass in data and some programs require two or more arguments
+besides the command itself (eg. sudiff) and don't use stdin.
+In this last case, we give selfdoc whenever too few arguments
+are given, since these usages violate the usual SU syntax.
+In all cases, selfdoc can be requested by giving only the
+program name.
+
+The flag argument distinguishes these cases:
+ flag = 0; fully defaulted, no stdin
+ flag = 1; usual case
+ flag = n > 1; no stdin and n extra args required
+
+pagedoc:
+Intended to be called by pagedoc(), but conceivably could be
+used directly as in:
+ if (xargc != 3) selfdoc();
+
+****************************************************************************
+Authors: Jack Cohen, Center for Wave Phenomena, 1993, based on on earlier
+versions by:
+SEP: Einar Kjartansson, Stew Levin CWP: Jack Cohen, Shuki Ronen
+HRC: Lyle
+****************************************************************************/
+{
+ switch(flag) {
+ case 1:
+ if (xargc == 1 && isatty(STDIN)) pagedoc();
+ break;
+ case 0:
+ if (xargc == 1 && isatty(STDIN) && isatty(STDOUT)) pagedoc();
+ break;
+ default:
+ if (xargc <= flag) pagedoc();
+ break;
+ }
+ return;
+}
+
+
+void pagedoc(void)
+{
+ extern char *sdoc[];
+ char **p = sdoc;
+ FILE *fp;
+
+ fflush(stdout);
+ fp = popen("more -22 1>&2", "w");
+ while(*p) (void)fprintf(fp, "%s\n", *p++);
+ pclose(fp);
+
+ exit(EXIT_FAILURE);
+}
+/*----------------------End of Package--------------------------------*/
diff --git a/raytime3d/getWaveletHeaders.c b/raytime3d/getWaveletHeaders.c
new file mode 100644
index 0000000000000000000000000000000000000000..5bff37528015722251741fcdb434db218e06ed90
--- /dev/null
+++ b/raytime3d/getWaveletHeaders.c
@@ -0,0 +1,52 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <math.h>
+#include "segy.h"
+
+/**
+* reads file which contain the source wavelets and reads receiver positions
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int getWaveletHeaders(char *file_src, int n1, int n2, float *gx, float *sx, float *gelev, float *selev, int verbose)
+{
+ FILE *fp;
+ size_t nread;
+ int ix;
+ size_t trace_sz;
+ off_t offset;
+ float scl, scll;
+ segy hdr;
+
+ if (file_src == NULL) return 0; /* Input pipe can not be handled */
+ else fp = fopen( file_src, "r" );
+ assert( fp != NULL);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+ if (hdr.scalel < 0) scll = 1.0/fabs(hdr.scalel);
+ else if (hdr.scalel == 0) scll = 1.0;
+ else scll = hdr.scalel;
+ trace_sz = (size_t)sizeof(float)*(n1)+TRCBYTES;
+
+ for (ix=0; ix<n2; ix++) {
+ offset = ix*trace_sz;
+ fseeko( fp, offset, SEEK_SET );
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ gx[ix] = hdr.gx*scl;
+ sx[ix] = hdr.sx*scl;
+ gelev[ix] = -1.0*hdr.gelev*scll;
+ selev[ix] = -1.0*hdr.selev*scll;
+ }
+ fclose(fp);
+ return 0;
+}
+
diff --git a/raytime3d/getWaveletInfo.c b/raytime3d/getWaveletInfo.c
new file mode 100644
index 0000000000000000000000000000000000000000..2f3734aae6c38e54653fab909ec5e936a157d8ce
--- /dev/null
+++ b/raytime3d/getWaveletInfo.c
@@ -0,0 +1,138 @@
+#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 "segy.h"
+
+/**
+* reads file which contain the source wavelets and computes sampling interval
+* and tries to estimate the maximum frequency content.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+typedef struct _dcomplexStruct { /* complex number */
+ double r,i;
+} dcomplex;
+#endif/* complex */
+
+int optncr(int n);
+void rc1fft(float *rdata, complex *cdata, int n, int sign);
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *fmax, int *nxm, int verbose)
+{
+ FILE *fp;
+ size_t nread, trace_sz;
+ off_t bytes;
+ int ret, one_shot, ntraces;
+ int optn, nfreq, i, iwmax;
+ float *trace;
+ float ampl, amplmax, tampl, tamplmax;
+ complex *ctrace;
+ segy hdr;
+
+ if (file_src == NULL) return 0; /* Input pipe can not be handled */
+ else fp = fopen( file_src, "r" );
+ assert( fp != NULL);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ ret = fseeko( fp, 0, SEEK_END );
+ if (ret<0) perror("fseeko");
+ bytes = ftello( fp );
+
+ *n1 = hdr.ns;
+ if (hdr.trid == 1 || hdr.dt != 0) {
+ *d1 = ((float) hdr.dt)*1.e-6;
+ *f1 = ((float) hdr.delrt)/1000.;
+ if (*d1 == 0.0) *d1 = hdr.d1;
+ }
+ else {
+ *d1 = hdr.d1;
+ *f1 = hdr.f1;
+ }
+ *f2 = hdr.f2;
+
+ trace_sz = (size_t)(sizeof(float)*(*n1)+TRCBYTES);
+ ntraces = (int) (bytes/trace_sz);
+ *n2 = ntraces;
+
+ /* check to find out number of traces in shot gather */
+
+ optn = optncr(*n1);
+ nfreq = optn/2 + 1;
+ ctrace = (complex *)malloc(nfreq*sizeof(complex));
+ one_shot = 1;
+ trace = (float *)malloc(optn*sizeof(float));
+ fseeko( fp, TRCBYTES, SEEK_SET );
+
+ while (one_shot) {
+ memset(trace,0,optn*sizeof(float));
+ nread = fread( trace, sizeof(float), *n1, fp );
+ assert (nread == *n1);
+ tamplmax = 0.0;
+ for (i=0;i<(*n1);i++) {
+ tampl = fabsf(trace[i]);
+ if (tampl > tamplmax) tamplmax = tampl;
+ }
+ if (trace[0]*1e-3 > tamplmax) {
+ fprintf(stderr,"WARNING: file_src has a large amplitude %f at t=0\n", trace[0]);
+ fprintf(stderr,"This will introduce high frequencies and can cause dispersion.\n");
+ }
+
+ /* estimate maximum frequency assuming amplitude spectrum is smooth */
+ rc1fft(trace,ctrace,optn,1);
+
+ /* find maximum amplitude */
+ amplmax = 0.0;
+ iwmax = 0;
+ for (i=0;i<nfreq;i++) {
+ ampl = sqrt(ctrace[i].r*ctrace[i].r+ctrace[i].i*ctrace[i].i);
+ if (ampl > amplmax) {
+ amplmax = ampl;
+ iwmax = i;
+ }
+ }
+ /* from the maximum amplitude position look for the largest frequency
+ * which has an amplitude 400 times weaker than the maximum amplitude */
+ for (i=iwmax;i<nfreq;i++) {
+ ampl = sqrt(ctrace[i].r*ctrace[i].r+ctrace[i].i*ctrace[i].i);
+ if (400*ampl < amplmax) {
+ *fmax = (i-1)*(1.0/(optn*(*d1)));
+ break;
+ }
+ }
+
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread==0) break;
+ }
+ *nxm = (int)ntraces;
+
+ if (verbose>2) {
+ vmess("For file %s", file_src);
+ vmess("nt=%d nx=%d", *n1, *n2);
+ vmess("dt=%f dx=%f", *d1, *d2);
+ vmess("fmax=%f", *fmax);
+ vmess("tstart=%f", *f1);
+ }
+
+ fclose(fp);
+ free(trace);
+ free(ctrace);
+
+ return 0;
+}
diff --git a/raytime3d/getpars.c b/raytime3d/getpars.c
new file mode 100644
index 0000000000000000000000000000000000000000..5099c5801bef214253daf07667c1b3c55b1008b1
--- /dev/null
+++ b/raytime3d/getpars.c
@@ -0,0 +1,732 @@
+/* This file is property of the Colorado School of Mines.
+
+ Copyright (C) 2007, Colorado School of Mines,
+ All rights reserved.
+
+
+ Redistribution and use in source and binary forms, with or
+ without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the following
+ disclaimer in the documentation and/or other materials provided
+ with the distribution.
+ * Neither the name of the Colorado School of Mines nor the names of
+ its contributors may be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ Warranty Disclaimer:
+ THIS SOFTWARE IS PROVIDED BY THE COLORADO SCHOOL OF MINES AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ COLORADO SCHOOL OF MINES OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+ STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+ IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+
+
+ Export Restriction Disclaimer:
+ We believe that CWP/SU: Seismic Un*x is a low technology product that does
+ not appear on the Department of Commerce CCL list of restricted exports.
+ Accordingly, we believe that our product meets the qualifications of
+ an ECCN (export control classification number) of EAR99 and we believe
+ it fits the qualifications of NRR (no restrictions required), and
+ is thus not subject to export restrictions of any variety.
+
+ Approved Reference Format:
+ In publications, please refer to SU as per the following example:
+ Cohen, J. K. and Stockwell, Jr. J. W., (200_), CWP/SU: Seismic Un*x
+ Release No. __: an open source software package for seismic
+ research and processing,
+ Center for Wave Phenomena, Colorado School of Mines.
+
+ Articles about SU in peer-reviewed journals:
+ Saeki, T., (1999), A guide to Seismic Un*x (SU)(2)---examples of data processing (part 1), data input and preparation of headers, Butsuri-Tansa (Geophysical Exploration), vol. 52, no. 5, 465-477.
+ Stockwell, Jr. J. W. (1999), The CWP/SU: Seismic Un*x Package, Computers and Geosciences, May 1999.
+ Stockwell, Jr. J. W. (1997), Free Software in Education: A case study of CWP/SU: Seismic Un*x, The Leading Edge, July 1997.
+ Templeton, M. E., Gough, C.A., (1998), Web Seismic Un*x: Making seismic reflection processing more accessible, Computers and Geosciences.
+
+ Acknowledgements:
+ SU stands for CWP/SU:Seismic Un*x, a processing line developed at Colorado
+ School of Mines, partially based on Stanford Exploration Project (SEP)
+ software.
+ */
+
+/*********************** self documentation **********************/
+/*****************************************************************************
+GETPARS - Functions to GET PARameterS from the command line. Numeric
+ parameters may be single values or arrays of int, uint,
+ short, ushort, long, ulong, float, or double. Single character
+ strings (type string or char *) may also be gotten.
+ Arrays of strings, delimited by, but not containing
+ commas are permitted.
+
+The functions are:
+
+initargs Makes command line args available to subroutines (re-entrant).
+ Every par program starts with this call!
+getparint get integers
+getparuint get unsigned integers
+getparshort get short integers
+getparushort get unsigned short integers
+getparlong get long integers
+getparulong get unsigned long integers
+getparfloat get float
+getpardouble get double
+getparstring get a single string
+getparstringarray get string array (fields delimited by commas)
+getpar get parameter by type
+getnparint get n'th occurrence of integer
+getnparuint get n'th occurrence of unsigned int
+getnparshort get n'th occurrence of short integer
+getnparushort get n'th occurrence of unsigned short int
+getnparlong get n'th occurrence of long integer
+getnparulong get n'th occurrence of unsigned long int
+getnparfloat get n'th occurrence of float integer
+getnpardouble get n'th occurrence of double integer
+getnparstring get n'th occurrence of string integer
+getnparstringarray get n'th occurrence of string integer array
+getnpar get n'th occurrence by type
+countparname return the number of times a parameter names is used
+countparval return the number of values in the last occurrence
+ of a parameter
+countnparval return the number of values in the n'th occurrence
+ of a parameter
+getPar Promax compatible version of getpar
+
+******************************************************************************
+Function Prototypes:
+void initargs (int argc, char **argv);
+int getparint (char *name, int *p);
+int getparuint (char *name, unsigned int *p);
+int getparshort (char *name, short *p);
+int getparushort (char *name, unsigned short *p);
+int getparlong (char *name, long *p);
+int getparulong (char *name, unsigned long *p);
+int getparfloat (char *name, float *p);
+int getpardouble (char *name, double *p);
+int getparstring (char *name, char **p);
+int getparstringarray (char *name, char **p);
+int getnparint (int n, char *name, int *p);
+int getnparuint (int n, char *name, unsigned int *p);
+int getnparshort (int n, char *name, short *p);
+int getnparushort (int n, char *name, unsigned short *p);
+int getnparlong (int n, char *name, long *p);
+int getnparulong (int n, char *name, unsigned long *p);
+int getnparfloat (int n, char *name, float *p);
+int getnpardouble (int n, char *name, double *p);
+int getnparstring (int n, char *name, char **p);
+int getnparstringarray (int n, char *name, char **p);
+int getnpar (int n, char *name, char *type, void *ptr);
+int countparname (char *name);
+int countparval (char *name);
+int countnparval (int n, char *name);
+void getPar(char *name, char *type, void *ptr);
+
+******************************************************************************
+Notes:
+Here are some usage examples:
+
+ ... if integer n not specified, then default to zero.
+ if (!getparint("n", &n)) n = 0;
+
+ ... if array of floats vx is specified, then
+ if (nx=countparval("vx")) {
+ ... allocate space for array
+ vx = (float *)malloc(nx*sizeof(float));
+ ... and get the floats
+ getparfloat("vx",vx);
+ }
+
+The command line for the above examples might look like:
+ progname n=35 vx=3.21,4,9.5
+ Every par program starts with this call!
+
+More examples are provided in the DTEST code at the end of this file.
+
+The functions: eatoh, eatou, eatol, eatov, eatoi, eatop used
+below are versions of atoi that check for overflow. The source
+file for these functions is atopkge.c.
+
+******************************************************************************
+Authors:
+Rob Clayton & Jon Claerbout, Stanford University, 1979-1985
+Shuki Ronen & Jack Cohen, Colorado School of Mines, 1985-1990
+Dave Hale, Colorado School of Mines, 05/29/90
+Credit to John E. Anderson for re-entrant initargs 03/03/94
+*****************************************************************************/
+/**************** end self doc ********************************/
+
+#include "par.h"
+
+#ifndef TRUE
+#define TRUE (1)
+#endif
+#ifndef FALSE
+#define FALSE (0)
+#endif
+
+/* parameter table */
+typedef struct {
+ char *name; /* external name of parameter */
+ char *asciival; /* ascii value of parameter */
+} pointer_table;
+
+/* global variables declared and used internally */
+static pointer_table *argtbl; /* parameter table */
+static int nargs; /* number of args that parse */
+static int tabled = FALSE; /* true when parameters tabled */
+static int targc; /* total number of args */
+static char **targv; /* pointer to arg strings */
+static char *argstr; /* storage for command line */
+
+/* functions declared and used internally */
+static int getparindex (int n, char *name);
+static void getparinit(void);
+static void tabulate (int argc, char **argv);
+static char *getpfname (void);
+static int white2null (char *str, int len);
+static int ccount (char c, char *s);
+static void strchop(char *s, char *t);
+
+/* make command line args available to subroutines -- re-entrant version */
+void initargs(int argc, char **argv)
+{
+ xargc = argc; xargv = argv;
+ if(tabled==TRUE){
+ free(argstr);
+ free(targv);
+ free(argtbl);
+ }
+ tabled = FALSE;
+ return;
+}
+
+/* functions to get values for the last occurrence of a parameter name */
+int getparint (char *name, int *ptr)
+{
+ return getnpar(0,name,"i",ptr);
+}
+int getparuint (char *name, unsigned int *ptr)
+{
+ return getnpar(0,name,"p",ptr);
+}
+int getparshort (char *name, short *ptr)
+{
+ return getnpar(0,name,"h",ptr);
+}
+int getparushort (char *name, unsigned short *ptr)
+{
+ return getnpar(0,name,"u",ptr);
+}
+int getparlong (char *name, long *ptr)
+{
+ return getnpar(0,name,"l",ptr);
+}
+int getparulong (char *name, unsigned long *ptr)
+{
+ return getnpar(0,name,"v",ptr);
+}
+int getparfloat (char *name, float *ptr)
+{
+ return getnpar(0,name,"f",ptr);
+}
+int getpardouble (char *name, double *ptr)
+{
+ return getnpar(0,name,"d",ptr);
+}
+int getparstring (char *name, char **ptr)
+{
+ return getnpar(0,name,"s",ptr);
+}
+int getparstringarray (char *name, char **ptr)
+{
+ return getnpar(0,name,"a",ptr);
+}
+int getpar (char *name, char *type, void *ptr)
+{
+ return getnpar(0,name,type,ptr);
+}
+
+/* functions to get values for the n'th occurrence of a parameter name */
+int getnparint (int n, char *name, int *ptr)
+{
+ return getnpar(n,name,"i",ptr);
+}
+int getnparuint (int n, char *name, unsigned int *ptr)
+{
+ return getnpar(n,name,"p",ptr);
+}
+int getnparshort (int n, char *name, short *ptr)
+{
+ return getnpar(n,name,"h",ptr);
+}
+int getnparushort (int n, char *name, unsigned short *ptr)
+{
+ return getnpar(n,name,"u",ptr);
+}
+int getnparlong (int n, char *name, long *ptr)
+{
+ return getnpar(n,name,"l",ptr);
+}
+int getnparulong (int n, char *name, unsigned long *ptr)
+{
+ return getnpar(n,name,"v",ptr);
+}
+int getnparfloat (int n, char *name, float *ptr)
+{
+ return getnpar(n,name,"f",ptr);
+}
+int getnpardouble (int n, char *name, double *ptr)
+{
+ return getnpar(n,name,"d",ptr);
+}
+int getnparstring (int n, char *name, char **ptr)
+{
+ return getnpar(n,name,"s",ptr);
+}
+int getnparstringarray (int n, char *name, char **ptr)
+{
+ return getnpar(n,name,"a",ptr);
+}
+int getnpar (int n, char *name, char *type, void *ptr)
+{
+ int i; /* index of name in symbol table */
+ int nval; /* number of parameter values found */
+ char *aval; /* ascii field of symbol */
+
+ if (xargc == 1) return 0;
+ if (!tabled) getparinit();/* Tabulate command line and parfile */
+ i = getparindex(n,name);/* Get parameter index */
+ if (i < 0) return 0; /* Not there */
+
+ /*
+ * handle string type as a special case, since a string
+ * may contain commas.
+ */
+ if (type[0]=='s') {
+ *((char**)ptr) = argtbl[i].asciival;
+ return 1;
+ }
+
+ /* convert vector of ascii values to numeric values */
+ for (nval=0,aval=argtbl[i].asciival; *aval; nval++) {
+ switch (type[0]) {
+ case 'i':
+ *(int*)ptr = eatoi(aval);
+ ptr = (int*)ptr+1;
+ break;
+ case 'p':
+ *(unsigned int*)ptr = eatop(aval);
+ ptr = (unsigned int*)ptr+1;
+ break;
+ case 'h':
+ *(short*)ptr = eatoh(aval);
+ ptr = (short*)ptr+1;
+ break;
+ case 'u':
+ *(unsigned short*)ptr = eatou(aval);
+ ptr = (unsigned short*)ptr+1;
+ break;
+ case 'l':
+ *(long*)ptr = eatol(aval);
+ ptr = (long*)ptr+1;
+ break;
+ case 'v':
+ *(unsigned long*)ptr = eatov(aval);
+ ptr = (unsigned long*)ptr+1;
+ break;
+ case 'f':
+ *(float*)ptr = eatof(aval);
+ ptr = (float*)ptr+1;
+ break;
+ case 'd':
+ *(double*)ptr = eatod(aval);
+ ptr = (double*)ptr+1;
+ break;
+ case 'a':
+ { char *tmpstr="";
+ tmpstr = (char *)calloc(strlen(aval)+1,1);
+
+ strchop(aval,tmpstr);
+ *(char**)ptr = tmpstr;
+ ptr=(char **)ptr + 1;
+ }
+ break;
+ default:
+ err("%s: invalid parameter type = %s",
+ __FILE__,type);
+ }
+ while (*aval++ != ',') {
+ if (!*aval) break;
+ }
+ }
+ return nval;
+}
+/* Promax compatible version of getnpar */
+void getPar(char *name, char *type, void *ptr)
+{
+ (void) getnpar(0,name,type,ptr);
+ return;
+}
+
+/* return number of occurrences of parameter name */
+int countparname (char *name)
+{
+ int i,nname;
+
+ if (xargc == 1) return 0;
+ if (!tabled) getparinit();
+ for (i=0,nname=0; i<nargs; ++i)
+ if (!strcmp(name,argtbl[i].name)) ++nname;
+ return nname;
+}
+
+/* return number of values in n'th occurrence of parameter name */
+int countnparval (int n, char *name)
+{
+ int i;
+
+ if (xargc == 1) return 0;
+ if (!tabled) getparinit();
+ i = getparindex(n,name);
+ if (i>=0)
+ return ccount(',',argtbl[i].asciival) + 1;
+ else
+ return 0;
+}
+
+/* return number of values in last occurrence of parameter name */
+int countparval (char *name)
+{
+ return countnparval(0,name);
+}
+
+
+
+/*
+ * Return the index of the n'th occurrence of a parameter name,
+ * except if n==0, return the index of the last occurrence.
+ * Return -1 if the specified occurrence does not exist.
+ */
+static int getparindex (int n, char *name)
+{
+ int i;
+ if (n==0) {
+ for (i=nargs-1; i>=0; --i)
+ if (!strcmp(name,argtbl[i].name)) break;
+ return i;
+ } else {
+ for (i=0; i<nargs; ++i)
+ if (!strcmp(name,argtbl[i].name))
+ if (--n==0) break;
+ if (i<nargs)
+ return i;
+ else
+ return -1;
+ }
+}
+
+/* Initialize getpar */
+static void getparinit (void)
+{
+ static char *pfname; /* name of parameter file */
+ FILE *pffd=NULL; /* file id of parameter file */
+ int pflen; /* length of parameter file in bytes */
+ static int pfargc; /* arg count from parameter file */
+ int parfile; /* parfile existence flag */
+ int argstrlen;
+ char *pargstr; /* storage for parameter file args */
+ int nread; /* bytes fread */
+ int i, j; /* counters */
+
+
+ tabled = TRUE; /* remember table is built */
+
+ /* Check if xargc was initiated */
+ if(!xargc)
+ err("%s: xargc=%d -- not initiated in main", __FILE__, xargc);
+
+ /* Space needed for command lines */
+ for (i = 1, argstrlen = 0; i < xargc; i++) {
+ argstrlen += strlen(xargv[i]) + 1;
+ }
+
+ /* Get parfile name if there is one */
+ /* parfile = (pfname = getpfname()) ? TRUE : FALSE; */
+ if ((pfname = getpfname())) {
+ parfile = TRUE;
+ } else {
+ parfile = FALSE;
+ }
+
+ if (parfile) {
+ pffd = fopen(pfname, "r");
+
+ /* Get the length */
+ fseek(pffd, 0, SEEK_END);
+ pflen = ftell(pffd);
+ rewind(pffd);
+ argstrlen += pflen;
+ } else {
+ pflen = 0;
+ }
+
+ /* Allocate space for command line and parameter file
+ plus nulls at the ends to help with parsing. */
+ /* argstr = (char *) calloc((size_t) (1+argstrlen+1), 1); */
+ /*argstr = (char *) ealloc1(1+argstrlen+1, 1);*/
+ argstr = (char *) calloc((size_t) (1+argstrlen+1), 1);
+
+ if (parfile) {
+ /* Read the parfile */
+ nread = fread(argstr + 1, 1, pflen, pffd);
+ if (nread != pflen) {
+ err("%s: fread only %d bytes out of %d from %s",
+ __FILE__, nread, pflen, pfname);
+ }
+ fclose(pffd);
+
+ /* Zap whites in parfile to help in parsing */
+ pfargc = white2null(argstr, pflen);
+
+ } else {
+ pfargc = 0;
+ }
+
+ /* Total arg count */
+ targc = pfargc + xargc - 1;
+
+ /* Allocate space for total arg pointers */
+ targv = (char **) calloc(targc, sizeof(char*));
+
+ if (parfile) {
+ /* Parse the parfile. Skip over multiple NULLs */
+ for (j = 1, i = 0; j < pflen; j++) {
+ if (argstr[j] && !argstr[j-1]) {
+ targv[i++] = argstr + j;
+ }
+ }
+ } else {
+ i = 0;
+ }
+
+ /* Copy command line arguments */
+ for (j = 1, pargstr = argstr + pflen + 2; j < xargc; j++) {
+ strcpy(pargstr,xargv[j]);
+ targv[i++] = pargstr;
+ pargstr += strlen(xargv[j]) + 1;
+ }
+
+ /* Allocate space for the pointer table */
+ argtbl = (pointer_table*) calloc(targc, sizeof(pointer_table));
+
+ /* Tabulate targv */
+ tabulate(targc, targv);
+
+ return;
+}
+
+#define PFNAME "par="
+/* Get name of parameter file */
+static char *getpfname (void)
+{
+ int i;
+ int pfnamelen;
+
+ pfnamelen = strlen(PFNAME);
+ for (i = xargc-1 ; i > 0 ; i--) {
+ if(!strncmp(PFNAME, xargv[i], pfnamelen)
+ && strlen(xargv[i]) != pfnamelen) {
+ return xargv[i] + pfnamelen;
+ }
+ }
+ return NULL;
+}
+
+#define iswhite(c) ((c) == ' ' || (c) == '\t' || (c) == '\n')
+
+/*
+ * Replace the whites by (possibly multiple) nulls. If we see a non-white
+ * and the previous char is a null, this signals the start of a string
+ * and we bump the count. This routine returns a count of the strings.
+ */
+static int white2null (char *str, int len)
+{
+ int i;
+ int count;
+ int inquote = FALSE;
+
+ str[0] = '\0'; /* This line added by Dave Hale, 1/30/96. */
+ for (i = 1, count = 0; i < len; i++) {
+ if (str[i]=='"') inquote=(inquote==TRUE)?FALSE:TRUE;
+ if (!inquote) {
+ if (iswhite(str[i])) { /* Is this a new word ? */
+ str[i] = '\0';
+ } else if (!str[i-1]) { /* multiple whites */
+ count++;
+ }
+ }
+ }
+ for (i = 1, inquote=FALSE; i < len; i++) {
+ if (str[i]=='"') inquote=(inquote==TRUE)?FALSE:TRUE;
+ if (inquote) {
+ if (str[i+1]!='"') {
+ str[i] = str[i+1];
+ } else {
+ str[i] = '\0';
+ str[i+1] = '\0';
+ inquote = FALSE;
+ }
+ }
+ }
+ str[len] = '\0';
+ return count;
+}
+
+/* Install symbol table */
+static void tabulate (int argc, char **argv)
+{
+ int i;
+ char *eqptr;
+
+ for (i = 0, nargs = 0 ; i < argc; i++) {
+ eqptr = (char *)strchr(argv[i], '=');
+ if (eqptr) {
+ argtbl[nargs].name = argv[i];
+ argtbl[nargs].asciival = eqptr + 1;
+ *eqptr = (char)0;
+
+ /* Debugging dump */
+/* fprintf(stderr, */
+/* "argtbl[%d]: name=%s asciival=%s\n", */
+/* nargs,argtbl[nargs].name,argtbl[nargs].asciival); */
+
+ nargs++;
+ }
+ }
+ return;
+}
+
+/* Count characters in a string */
+static int ccount (char c, char *s)
+{
+ int i, count;
+ for (i = 0, count = 0; s[i] != 0; i++)
+ if(s[i] == c) count++;
+ return count;
+}
+
+static void strchop(char *s, char *t)
+/***********************************************************************
+strchop - chop off the tail end of a string "s" after a "," returning
+ the front part of "s" as "t".
+************************************************************************
+Notes:
+Based on strcpy in Kernighan and Ritchie's C [ANSI C] book, p. 106.
+************************************************************************
+Author: CWP: John Stockwell and Jack K. Cohen, July 1995
+***********************************************************************/
+{
+
+ while ( (*s != ',') && (*s != '\0') ) {
+ *t++ = *s++;
+ }
+ *t='\0';
+}
+
+
+#ifdef TEST
+#define N 100
+main(int argc, char **argv)
+{
+ char *s;
+ short h, vh[N];
+ unsigned short u, vu[N];
+ long l, vl[N];
+ unsigned long v, vv[N];
+ int i, vi[N], ipar, npar, nval;
+ unsigned int p, vp[N];
+ float f, vf[N];
+ double d, vd[N];
+
+ initargs(argc, argv);
+
+ /* int parameters */
+ npar = countparname("i");
+ printf("\nnumber of i pars = %d\n",npar);
+ for (ipar=1; ipar<=npar; ++ipar) {
+ getnparint(ipar,"i",&i);
+ printf("occurrence %d of i=%d\n",ipar,i);
+ }
+ if (getparint("i", &i))
+ printf("last occurrence of i=%d\n",i);
+ npar = countparname("vi");
+ printf("number of vi pars = %d\n",npar);
+ for (ipar=1; ipar<=npar; ++ipar) {
+ nval = countnparval(ipar,"vi");
+ printf("occurrence %d has %d values\n",ipar,nval);
+ nval = getnparint(ipar,"vi",vi);
+ printf("vi=");
+ for (i=0; i<nval; i++)
+ printf("%d%c",vi[i],i==nval-1?'\n':',');
+ }
+ if (npar>0) {
+ nval = countparval("vi");
+ printf("last occurrence has %d values\n",nval);
+ getparint("vi",vi);
+ printf("vi=");
+ for (i=0; i<nval; i++)
+ printf("%d%c",vi[i],i==nval-1?'\n':',');
+ }
+
+ /* float parameters */
+ npar = countparname("f");
+ printf("\nnumber of f pars = %d\n",npar);
+ for (ipar=1; ipar<=npar; ++ipar) {
+ getnparfloat(ipar,"f",&f);
+ printf("occurrence %d of f=%g\n",ipar,f);
+ }
+ if (getparfloat("f", &f))
+ printf("last occurrence of f=%g\n",f);
+ npar = countparname("vf");
+ printf("number of vf pars = %d\n",npar);
+ for (ipar=1; ipar<=npar; ++ipar) {
+ nval = countnparval(ipar,"vf");
+ printf("occurrence %d has %d values\n",ipar,nval);
+ nval = getnparfloat(ipar,"vf",vf);
+ printf("vf=");
+ for (i=0; i<nval; i++)
+ printf("%g%c",vf[i],i==nval-1?'\n':',');
+ }
+ if (npar>0) {
+ nval = countparval("vf");
+ printf("last occurrence has %d values\n",nval);
+ getparfloat("vf",vf);
+ printf("vf=");
+ for (i=0; i<nval; i++)
+ printf("%g%c",vf[i],i==nval-1?'\n':',');
+ }
+
+ /* string parameters */
+ npar = countparname("s");
+ printf("\nnumber of s pars = %d\n",npar);
+ for (ipar=1; ipar<=npar; ++ipar) {
+ getnparstring(ipar,"s",&s);
+ printf("occurrence %d of s=%s\n",ipar,s);
+ }
+ if (getparstring("s", &s))
+ printf("last occurrence of s=%s\n",s);
+
+ return EXIT_SUCCESS;
+}
+#endif
+
diff --git a/raytime3d/name_ext.c b/raytime3d/name_ext.c
new file mode 100644
index 0000000000000000000000000000000000000000..8fa1e09d254153d4d783ab040cebfada4d82d3b7
--- /dev/null
+++ b/raytime3d/name_ext.c
@@ -0,0 +1,44 @@
+#include<stdlib.h>
+#include<string.h>
+#include<stdio.h>
+
+/**
+* inserts a character string after the filename, before the extension
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+
+void name_ext(char *filename, char *extension)
+{
+ char ext[100];
+
+ if (strstr(filename, ".su") != NULL) {
+ sprintf(ext,"%s.su", extension);
+ strcpy(strstr(filename, ".su"), ext);
+ }
+ else if (strstr(filename, ".segy") != NULL) {
+ sprintf(ext,"%s.segy", extension);
+ strcpy(strstr(filename, ".segy"), ext);
+ }
+ else if (strstr(filename, ".mat") != NULL) {
+ sprintf(ext,"%s.mat", extension);
+ strcpy(strstr(filename, ".mat"), ext);
+ }
+ else if (strstr(filename, ".hdf") != NULL) {
+ sprintf(ext,"%s.hdf", extension);
+ strcpy(strstr(filename, ".hdf"), ext);
+ }
+ else if (strrchr(filename, '.') != NULL) {
+ sprintf(ext,"%s.su", extension);
+ strcpy(strrchr(filename, '.'), ext);
+ }
+ else {
+ sprintf(ext,"%s.su", extension);
+ strcat(filename, ext);
+ }
+
+ return;
+}
diff --git a/raytime3d/par.h b/raytime3d/par.h
new file mode 100644
index 0000000000000000000000000000000000000000..fce76ed344382c6d5719737e5395fd8fb3ad0a5b
--- /dev/null
+++ b/raytime3d/par.h
@@ -0,0 +1,217 @@
+/* This file is property of the Colorado School of Mines.
+
+ Copyright © 2007, Colorado School of Mines,
+ All rights reserved.
+
+
+ Redistribution and use in source and binary forms, with or
+ without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ * Redistributions in binary form must reproduce the above
+ copyright notice, this list of conditions and the following
+ disclaimer in the documentation and/or other materials provided
+ with the distribution.
+ * Neither the name of the Colorado School of Mines nor the names of
+ its contributors may be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ Warranty Disclaimer:
+ THIS SOFTWARE IS PROVIDED BY THE COLORADO SCHOOL OF MINES AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ COLORADO SCHOOL OF MINES OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+ STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
+ IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+
+
+ Export Restriction Disclaimer:
+ We believe that CWP/SU: Seismic Un*x is a low technology product that does
+ not appear on the Department of Commerce CCL list of restricted exports.
+ Accordingly, we believe that our product meets the qualifications of
+ an ECCN (export control classification number) of EAR99 and we believe
+ it fits the qualifications of NRR (no restrictions required), and
+ is thus not subject to export restrictions of any variety.
+
+ Approved Reference Format:
+ In publications, please refer to SU as per the following example:
+ Cohen, J. K. and Stockwell, Jr. J. W., (200_), CWP/SU: Seismic Un*x
+ Release No. __: an open source software package for seismic
+ research and processing,
+ Center for Wave Phenomena, Colorado School of Mines.
+
+ Articles about SU in peer-reviewed journals:
+ Saeki, T., (1999), A guide to Seismic Un*x (SU)(2)---examples of data processing (part 1), data input and preparation of headers, Butsuri-Tansa (Geophysical Exploration), vol. 52, no. 5, 465-477.
+ Stockwell, Jr. J. W. (1999), The CWP/SU: Seismic Un*x Package, Computers and Geosciences, May 1999.
+ Stockwell, Jr. J. W. (1997), Free Software in Education: A case study of CWP/SU: Seismic Un*x, The Leading Edge, July 1997.
+ Templeton, M. E., Gough, C.A., (1998), Web Seismic Un*x: Making seismic reflection processing more accessible, Computers and Geosciences.
+
+ Acknowledgements:
+ SU stands for CWP/SU:Seismic Un*x, a processing line developed at Colorado
+ School of Mines, partially based on Stanford Exploration Project (SEP)
+ software.
+ . */
+/* par.h - include file for getpar, selfdoc, and error handling functions */
+
+#ifndef PAR_H
+#define PAR_H
+void verr(char *fmt, ...);
+void vmess(char *fmt, ...);
+void vwarn(char *fmt, ...);
+void vsyserr(char *fmt, ...);
+
+/* TYPEDEFS */
+typedef union { /* storage for arbitrary type */
+ char s[8];
+ short h;
+ unsigned short u;
+ long l;
+ unsigned long v;
+ int i;
+ unsigned int p;
+ float f;
+ double d;
+ unsigned int U:16;
+ unsigned int P:32;
+} Value;
+
+/* INCLUDES */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <fcntl.h> /* non-ANSI */
+#include <unistd.h> /* non-ANSI */
+#include <sys/types.h> /* non-ANSI */
+#include<string.h>
+
+
+/* GLOBAL DECLARATIONS */
+extern int xargc; extern char **xargv;
+
+
+/* TYPEDEFS */
+typedef char *cwp_String;
+
+typedef enum {BADFILETYPE = -1,
+ TTY, DISK, DIRECTORY, TAPE, PIPE, FIFO, SOCKET, SYMLINK} FileType;
+
+/* DEFINES */
+
+/* getpar macros */
+#define MUSTGETPARINT(x,y) if(!getparint(x,y)) err("must specify %s=",x)
+#define MUSTGETPARFLOAT(x,y) if(!getparfloat(x,y)) err("must specify %s=",x)
+#define MUSTGETPARSTRING(x,y) if(!getparstring(x,y)) err("must specify %s=",x)
+
+#define STDIN (0)
+#define STDOUT (1)
+#define STDERR (2)
+
+/* FUNCTION PROTOTYPES */
+
+#ifdef __cplusplus /* if C++, specify external C linkage */
+extern "C" {
+#endif
+
+/* getpar parameter parsing */
+void initargs (int argc, char **argv);
+int getparint (char *name, int *p);
+int getparuint (char *name, unsigned int *p);
+int getparshort (char *name, short *p);
+int getparushort (char *name, unsigned short *p);
+int getparlong (char *name, long *p);
+int getparulong (char *name, unsigned long *p);
+int getparfloat (char *name, float *p);
+int getpardouble (char *name, double *p);
+int getparstring (char *name, char **p);
+int getparstringarray (char *name, char **p);
+int getnparint (int n, char *name, int *p);
+int getnparuint (int n, char *name, unsigned int *p);
+int getnparshort (int n, char *name, short *p);
+int getnparushort (int n, char *name, unsigned short *p);
+int getnparlong (int n, char *name, long *p);
+int getnparulong (int n, char *name, unsigned long *p);
+int getnparfloat (int n, char *name, float *p);
+int getnpardouble (int n, char *name, double *p);
+int getnparstring (int n, char *name, char **p);
+int getnparstringarray (int n, char *name, char **p);
+int getnpar (int n, char *name, char *type, void *ptr);
+int countparname (char *name);
+int countparval (char *name);
+int countnparval (int n, char *name);
+
+/* For ProMAX */
+void getPar(char *name, char *type, void *ptr);
+
+/* errors and warnings */
+void err (char *fmt, ...);
+void syserr (char *fmt, ...);
+void warn (char *fmt, ...);
+
+/* self documentation */
+void pagedoc (void);
+void requestdoc (int i);
+
+/* system calls with error trapping */
+int ecreat(char *path, int perms);
+int efork(void);
+int eopen(char *path, int flags, int perms);
+int eclose(int fd);
+int eunlink(char *path);
+long elseek(int fd, long offset, int origin);
+int epipe(int fd[2]);
+int eread(int fd, char *buf, int nbytes);
+int ewrite(int fd, char *buf, int nbytes);
+
+/* system subroutine calls with error trapping */
+FILE *efopen(const char *file, const char *mode);
+FILE *efreopen(const char *file, const char *mode, FILE *stream1);
+FILE *efdopen(int fd, const char *mode);
+FILE *epopen(char *command, char *type);
+int efclose(FILE *stream);
+int epclose(FILE *stream);
+int efflush(FILE *stream);
+int eremove(const char *file);
+int erename(const char *oldfile, const char* newfile);
+int efseek(FILE *stream, long offset, int origin);
+void erewind(FILE *stream);
+long eftell(FILE *stream);
+FILE *etmpfile(void);
+char *etmpnam(char *namebuffer);
+void *emalloc(size_t size);
+void *erealloc(void *memptr, size_t size);
+void *ecalloc(size_t count, size_t size);
+size_t efread(void *bufptr, size_t size, size_t count, FILE *stream);
+size_t efwrite(void *bufptr, size_t size, size_t count, FILE *stream);
+
+#ifndef SUN_A
+int efgetpos(FILE *stream, fpos_t *position);
+int efsetpos(FILE *stream, const fpos_t *position);
+#endif
+
+/* string to numeric conversion with error checking */
+short eatoh(char *s);
+unsigned short eatou(char *s);
+int eatoi(char *s);
+unsigned int eatop(char *s);
+long eatol(char *s);
+unsigned long eatov(char *s);
+float eatof(char *s);
+double eatod(char *s);
+
+/* file type checking */
+FileType filestat(int fd);
+char *printstat(int fd);
+
+#ifdef __cplusplus /* if C++ (external C linkage is being specified) */
+}
+#endif
+
+#endif /* PAR_H */
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..73e57f41822da6bec2322f85b0b943e25afad73f
--- /dev/null
+++ b/raytime3d/raytime3d.c
@@ -0,0 +1,293 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
+#include <genfft.h>
+#include"par.h"
+#include"raytime3d.h"
+#include "segy.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+double wallclock_time(void);
+
+void name_ext(char *filename, char *extension);
+
+void threadAffinity(void);
+
+
+int getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose);
+
+int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord r, rayPar ray, fcoord *T, float *Jr);
+
+void applyMovingAverageFilter(float *slowness, icoord size, int window, int dim, float *averageModel);
+
+int readModel3d(char *file_name, float *slowness, int n1, int n2, int n3, int nz, int nx, int ny, float h, int verbose);
+
+int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int reverse, int verbose);
+
+int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
+
+void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int xs, int ys, int zs, int NCUBE);
+
+void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube);
+
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" RAYTIME3D - modeling of one-way traveltime for operators in 3D media",
+" ",
+" raytime3d file_cp= xsrc1= zsrc1= ysrc1= [optional parameters]",
+" ",
+" Required parameters:",
+" ",
+" file_cp= ................ gridded velocity file ",
+" file_src= ......... file with source signature",
+" file_rcv=recv.su .. base name for receiver files",
+" file_rcvtime= ..... receiver file in x-t",
+" h= ................ read from model file: if d1==0 then h= can be used to set it",
+" nt=1024 ........... number of time-samples in file_rcvtime",
+" xsrc1= ................... x-position of the source (m)",
+" ysrc1= ................... y-position of the source (m)",
+" zsrc1= ................... z-position of the source (m)",
+" ",
+" Optional parameters:",
+" ",
+" INPUT AND OUTPUT ",
+" key=gy ................... input data sorting key",
+" nx=1 ..................... if 1D file number of points in x ",
+" ny=1 ..................... if 2D file number of points in y ",
+" file_out= ................ output file with traveltime cube",
+" file_amp= ................ output file with approximate amplitudes",
+" ",
+//" RAY TRACING PARAMETERS:",
+//" dT=0 ..................... put traces on one-way time grid with step dT",
+//" Tmin=first shot........... minimum time of one-way time grid",
+//" Tmax=last shot ........... maximum time of one-way time grid",
+//" hom=1 .................... 1: draw straight rays in homogeneous layers",
+//" ",
+" SOURCE POSITIONS ",
+" xsrc2=xsrc1 .............. x-position of last source",
+" dxsrc=0 .................. step in source x-direction",
+" ysrc2=ysrc1 .............. y-position of last source",
+" dysrc=0 .................. step in source y-direction",
+" zsrc2=zsrc1 .............. z-position of last source",
+" dzsrc=0 .................. step in source z-direction",
+" RECEIVER POSITIONS ",
+" xrcv=-(nx-1/2)*h,(nx-1/2)*h .. x-position's of receivers (array)",
+" yrcv=-(ny-1)/2*h,(ny-1/2)*h .. y-position's of receivers (array)",
+" zrcv=0,0 ................. z-position's of receivers (array)",
+" dxrcv=h ................. step in receiver x-direction",
+" dyrcv=h ................. step in receiver y-direction",
+" dzrcv=0 .................. step in receiver z-direction",
+" dxspr=0 .................. step of receiver spread in x-direction",
+" dyspr=0 .................. step of receiver spread in y-direction",
+" lint=1 ................... linear interpolate between the rcv points",
+" verbose=0 ................ verbose option",
+" ",
+" raytime3d calculates the first arrival time at the defined receiver array ",
+" for the defined shots at different depth and lateral positions.",
+" Every new lateral position (with dxsrc) gives a new output gather.",
+" ",
+" PROGRAM TO CALCULATE TRAVEL TIMES IN 3D MEDIA ",
+" AUTHORS: John E. Vidale(1986-1989), J. Hole(1990-1995) ",
+" ",
+" Translated to DELPHI environment: Jan Thorbecke 17-04-1996",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+#define SQR2 1.414213562
+#define SQR3 1.732050808
+#define SQR6 2.449489743
+#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
+#define s0(x,y,z) slow0[nxy*(z) + nx*(y) + (x)]
+
+void main(int argc, char *argv[])
+{
+ modPar mod;
+ recPar rec;
+ srcPar src;
+ shotPar shot;
+ rayPar ray;
+ int
+ nx, /* x-dimension of mesh (LEFT-TO-RIGHT) */
+ ny, /* y-dimension of mesh (FRONT-TO-BACK) */
+ nz, /* z-dimension of mesh (TOP-TO-BOTTOM) */
+ nxy, nxyz, xs, ys, zs, cube,
+ xx, yy, zz, i, j;
+ float
+ h, /* spatial mesh interval (units consistant with vel) */
+ *slow0, *time0;
+
+/* to read the velocity file */
+ int error, n1, n2, n3, ret, size, nkeys, verbose;
+ float d1, d2, d3, f1, f2, f3, *tmpdata, c, scl, ox, oz, oy;
+ char *file_cp, *file_out;
+ segy *hdrs;
+
+/*---------------------------------------------------------------------------*
+ * Read input parameters and query for any that are needed but missing.
+ *---------------------------------------------------------------------------*/
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparint("verbose",&verbose)) verbose=0;
+ if (verbose) {
+ vmess("Hole, J.A., and B.C. Zelt, 1995. \"3-D finite-difference");
+ vmess("reflection traveltimes\". Geophys. J. Int., 121, 427-434");
+ }
+ if(!getparstring("file_out",&file_out)) verr("file_out not given");
+
+ getParameters3d(&mod, &rec, &src, &shot, &ray, verbose);
+
+
+/*---------------------------------------------------------------------------*
+ * Open velocity file
+ *---------------------------------------------------------------------------*/
+
+ if (file_cp != NULL) {
+
+ if (n2==1) { /* 1D model */
+ if(!getparint("nx",&nx)) verr("for 1D medium nx not defined");
+ if(!getparint("ny",&nx)) verr("for 1D medium ny not defined");
+ nz = n1;
+ oz = f1; ox = ((nx-1)/2)*d1; oy = ((ny-1)/2)*d1;
+ dz = d1; dx = d1; dy = d1;
+ }
+ else if (n3==1) { /* 2D model */
+ if(!getparint("ny",&nx)) verr("for 2D medium ny not defined");
+ nz = n1; nx = n2;
+ oz = f1; ox = f2; oy = ((ny-1)/2)*d1;
+ dz = d1; dx = d1; dy = d1;
+ }
+ else { /* Full 3D model */
+ nz = n1; nx = n2; nz = n3;
+ oz = f1; ox = f2; oy = f3;
+ dz = d1; dx = d1; dy = d1;
+ }
+
+ h = d1;
+ slow0 = (float *)malloc(nz*nx*ny*sizeof(float));
+ if (slow0 == NULL) verr("Out of memory for slow0 array!");
+
+ readModel3d(file_cp, slow0, n1, n2, n3, nz, nx, ny, h, verbose);
+
+ if (verbose) vmess("h = %.2f nx = %d nz = %d ny = %d", h, nx, nz, ny);
+
+ }
+ else {
+ nxy = nx * ny;
+ if(!getparint("nx",&nx)) verr("for homogenoeus medium nx not defined");
+ if(!getparint("ny",&nx)) verr("for homogenoeus medium ny not defined");
+ if(!getparint("nz",&nx)) verr("for homogenoeus medium nz not defined");
+ oz = 0; ox = ((nx-1)/2)*d1; oy = ((ny-1)/2)*d1;
+
+ slow0 = (float *)malloc(nx*nz*ny*sizeof(float));
+ if (slow0 == NULL) verr("Out of memory for slow0 array!");
+ scl = h/c;
+ ox = 0; oy = 0; oz = 0;
+ for (zz = 0; zz < nz; zz++) {
+ for (yy = 0; yy < ny; yy++) {
+ for (xx = 0; xx < nx; xx++)
+ slow0[zz*nxy+yy*nx+xx] = scl;
+ }
+ }
+ }
+
+ nxy = nx * ny;
+ nxyz = nx * ny * nz;
+
+ /* ALLOCATE MAIN GRID FOR TIMES */
+ time0 = (float *) malloc(sizeof(float)*nxyz);
+ if(time0 == NULL) verr("error in allocation of array time0");
+
+/*---------------------------------------------------------------------------*
+ * Input the source locations.
+ * and
+ * Initialize the traveltime array. Place t=0 at source position.
+ *---------------------------------------------------------------------------*/
+
+ src3d(time0, slow0, nz, nx, ny, h, ox, oy, oz, &xs, &ys, &zs, &cube);
+ if (verbose) vmess("source positions xs = %d ys = %d zs = %d", xs,ys,zs);
+
+/* for (zz = 0; zz < nz; zz++) {
+ for (yy = 0; yy < ny; yy++) {
+ for (xx = 0; xx < nx; xx++)
+ if (time0[zz*nxy+yy*nx+xx] != 1e10) fprintf(stderr,"slow[%d,%d,%d] = %f\n", xx,yy,zz, time0[zz*nxy+yy*nx+xx]);
+ }
+ }
+*/
+
+/*---------------------------------------------------------------------------*
+ * Read in receiver positions
+ *---------------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------------------*
+ * Check and set parameters
+ *---------------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------------------*
+ * Compute traveltimes.
+ *---------------------------------------------------------------------------*/
+
+ vidale3d(slow0, time0, nz, nx, ny, h, xs, ys, zs, cube);
+
+/*---------------------------------------------------------------------------*
+ * Write output
+ *---------------------------------------------------------------------------*/
+
+/*
+ for (zz = 0; zz < nz; zz++) {
+ for (yy = 0; yy < ny; yy++) {
+ for (xx = 0; xx < nx; xx++)
+ if (time0[zz*nxy+yy*nx+xx] != 1e10) fprintf(stderr,"slow[%d,%d,%d] = %f\n", xx,yy,zz, time0[zz*nxy+yy*nx+xx]);
+ }
+ }
+*/
+// ret = open_file(file_out, GUESS_TYPE, DELPHI_CREATE);
+// if (ret < 0 ) verr("error in creating output file %s", file_out);
+
+ hdrs = (segy *) malloc(ny*sizeof(segy));
+ tmpdata = (float *)malloc(nxy*sizeof(float));
+ f1 = ox;
+ f2 = oy;
+ d1 = h;
+ d2 = h;
+
+// gen_hdrs(hdrs,nx,ny,f1,f2,d1,d2,TRID_ZX);
+ for (i = 0; i < ny; i++) {
+ hdrs[i].scalco = -1000;
+ hdrs[i].scalel = -1000;
+/* hdrs[i].offset = xi[0]*dx + is*ispr*dx - xsrc;*/
+ hdrs[i].sx = (int)(ox+xs*h)*1000;
+ hdrs[i].sy = (int)(oy+ys*h)*1000;
+ hdrs[i].gy = (int)(oy+i*d2)*1000;
+ hdrs[i].sdepth = (int)(oz+zs*h)*1000;
+ hdrs[i].fldr = 1;
+ hdrs[i].trwf = ny;
+ for (j = 0; j < nx; j++) {
+ tmpdata[i*nx+j] = time0[i*nx+j];
+ }
+ }
+
+/*
+ ret = write_data(file_out,tmpdata,nx,ny,f1,f2,d1,d2,type,hdrs);
+ if (ret < 0 ) verr("error on writing output file.");
+ ret = close_file(file_out);
+ if (ret < 0) verr("err %d on closing output file",ret);
+*/
+
+ free(time0);
+ free(slow0);
+ free(hdrs);
+ free(tmpdata);
+
+ exit(0);
+
+}
diff --git a/raytime3d/segy.h b/raytime3d/segy.h
new file mode 100644
index 0000000000000000000000000000000000000000..d0a0d769d1548115b04396076b4a15d5be0ee687
--- /dev/null
+++ b/raytime3d/segy.h
@@ -0,0 +1,849 @@
+/* Copyright (c) Colorado School of Mines, 2011.*/
+/* All rights reserved. */
+
+/* segy.h - include file for SEGY traces
+ *
+ * declarations for:
+ * typedef struct {} segy - the trace identification header
+ * typedef struct {} bhed - binary header
+ *
+ * Note:
+ * If header words are added, run the makefile in this directory
+ * to recreate hdr.h.
+ *
+ * Reference:
+ * K. M. Barry, D. A. Cavers and C. W. Kneale, "Special Report:
+ * Recommended Standards for Digital Tape Formats",
+ * Geophysics, vol. 40, no. 2 (April 1975), P. 344-352.
+ *
+ * $Author: john $
+ * $Source: /usr/local/cwp/src/su/include/RCS/segy.h,v $
+ * $Revision: 1.33 $ ; $Date: 2011/11/11 23:56:14 $
+ */
+
+#include <limits.h>
+#include "par.h"
+
+#ifndef SEGY_H
+#define SEGY_H
+#define TRCBYTES 240
+
+#define SU_NFLTS 32767 /* Arbitrary limit on data array size */
+
+
+/* TYPEDEFS */
+typedef struct { /* segy - trace identification header */
+
+ int tracl; /* Trace sequence number within line
+ --numbers continue to increase if the
+ same line continues across multiple
+ SEG Y files.
+ byte# 1-4
+ */
+
+ int tracr; /* Trace sequence number within SEG Y file
+ ---each file starts with trace sequence
+ one
+ byte# 5-8
+ */
+
+ int fldr; /* Original field record number
+ byte# 9-12
+ */
+
+ int tracf; /* Trace number within original field record
+ byte# 13-16
+ */
+
+ int ep; /* energy source point number
+ ---Used when more than one record occurs
+ at the same effective surface location.
+ byte# 17-20
+ */
+
+ int cdp; /* Ensemble number (i.e. CDP, CMP, CRP,...)
+ byte# 21-24
+ */
+
+ int cdpt; /* trace number within the ensemble
+ ---each ensemble starts with trace number one.
+ byte# 25-28
+ */
+
+ short trid; /* trace identification code:
+ -1 = Other
+ 0 = Unknown
+ 1 = Seismic data
+ 2 = Dead
+ 3 = Dummy
+ 4 = Time break
+ 5 = Uphole
+ 6 = Sweep
+ 7 = Timing
+ 8 = Water break
+ 9 = Near-field gun signature
+ 10 = Far-field gun signature
+ 11 = Seismic pressure sensor
+ 12 = Multicomponent seismic sensor
+ - Vertical component
+ 13 = Multicomponent seismic sensor
+ - Cross-line component
+ 14 = Multicomponent seismic sensor
+ - in-line component
+ 15 = Rotated multicomponent seismic sensor
+ - Vertical component
+ 16 = Rotated multicomponent seismic sensor
+ - Transverse component
+ 17 = Rotated multicomponent seismic sensor
+ - Radial component
+ 18 = Vibrator reaction mass
+ 19 = Vibrator baseplate
+ 20 = Vibrator estimated ground force
+ 21 = Vibrator reference
+ 22 = Time-velocity pairs
+ 23 ... N = optional use
+ (maximum N = 32,767)
+
+ Following are CWP id flags:
+
+ 109 = autocorrelation
+ 110 = Fourier transformed - no packing
+ xr[0],xi[0], ..., xr[N-1],xi[N-1]
+ 111 = Fourier transformed - unpacked Nyquist
+ xr[0],xi[0],...,xr[N/2],xi[N/2]
+ 112 = Fourier transformed - packed Nyquist
+ even N:
+ xr[0],xr[N/2],xr[1],xi[1], ...,
+ xr[N/2 -1],xi[N/2 -1]
+ (note the exceptional second entry)
+ odd N:
+ xr[0],xr[(N-1)/2],xr[1],xi[1], ...,
+ xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2]
+ (note the exceptional second & last entries)
+ 113 = Complex signal in the time domain
+ xr[0],xi[0], ..., xr[N-1],xi[N-1]
+ 114 = Fourier transformed - amplitude/phase
+ a[0],p[0], ..., a[N-1],p[N-1]
+ 115 = Complex time signal - amplitude/phase
+ a[0],p[0], ..., a[N-1],p[N-1]
+ 116 = Real part of complex trace from 0 to Nyquist
+ 117 = Imag part of complex trace from 0 to Nyquist
+ 118 = Amplitude of complex trace from 0 to Nyquist
+ 119 = Phase of complex trace from 0 to Nyquist
+ 121 = Wavenumber time domain (k-t)
+ 122 = Wavenumber frequency (k-omega)
+ 123 = Envelope of the complex time trace
+ 124 = Phase of the complex time trace
+ 125 = Frequency of the complex time trace
+ 130 = Depth-Range (z-x) traces
+ 201 = Seismic data packed to bytes (by supack1)
+ 202 = Seismic data packed to 2 bytes (by supack2)
+ byte# 29-30
+ */
+
+ short nvs; /* Number of vertically summed traces yielding
+ this trace. (1 is one trace,
+ 2 is two summed traces, etc.)
+ byte# 31-32
+ */
+
+ short nhs; /* Number of horizontally summed traces yielding
+ this trace. (1 is one trace
+ 2 is two summed traces, etc.)
+ byte# 33-34
+ */
+
+ short duse; /* Data use:
+ 1 = Production
+ 2 = Test
+ byte# 35-36
+ */
+
+ int offset; /* Distance from the center of the source point
+ to the center of the receiver group
+ (negative if opposite to direction in which
+ the line was shot).
+ byte# 37-40
+ */
+
+ int gelev; /* Receiver group elevation from sea level
+ (all elevations above the Vertical datum are
+ positive and below are negative).
+ byte# 41-44
+ */
+
+ int selev; /* Surface elevation at source.
+ byte# 45-48
+ */
+
+ int sdepth; /* Source depth below surface (a positive number).
+ byte# 49-52
+ */
+
+ int gdel; /* Datum elevation at receiver group.
+ byte# 53-56
+ */
+
+ int sdel; /* Datum elevation at source.
+ byte# 57-60
+ */
+
+ int swdep; /* Water depth at source.
+ byte# 61-64
+ */
+
+ int gwdep; /* Water depth at receiver group.
+ byte# 65-68
+ */
+
+ short scalel; /* Scalar to be applied to the previous 7 entries
+ to give the real value.
+ Scalar = 1, +10, +100, +1000, +10000.
+ If positive, scalar is used as a multiplier,
+ if negative, scalar is used as a divisor.
+ byte# 69-70
+ */
+
+ short scalco; /* Scalar to be applied to the next 4 entries
+ to give the real value.
+ Scalar = 1, +10, +100, +1000, +10000.
+ If positive, scalar is used as a multiplier,
+ if negative, scalar is used as a divisor.
+ byte# 71-72
+ */
+
+ int sx; /* Source coordinate - X
+ byte# 73-76
+ */
+
+ int sy; /* Source coordinate - Y
+ byte# 77-80
+ */
+
+ int gx; /* Group coordinate - X
+ byte# 81-84
+ */
+
+ int gy; /* Group coordinate - Y
+ byte# 85-88
+ */
+
+ short counit; /* Coordinate units: (for previous 4 entries and
+ for the 7 entries before scalel)
+ 1 = Length (meters or feet)
+ 2 = Seconds of arc
+ 3 = Decimal degrees
+ 4 = Degrees, minutes, seconds (DMS)
+
+ In case 2, the X values are longitude and
+ the Y values are latitude, a positive value designates
+ the number of seconds east of Greenwich
+ or north of the equator
+
+ In case 4, to encode +-DDDMMSS
+ counit = +-DDD*10^4 + MM*10^2 + SS,
+ with scalco = 1. To encode +-DDDMMSS.ss
+ counit = +-DDD*10^6 + MM*10^4 + SS*10^2
+ with scalco = -100.
+ byte# 89-90
+ */
+
+ short wevel; /* Weathering velocity.
+ byte# 91-92
+ */
+
+ short swevel; /* Subweathering velocity.
+ byte# 93-94
+ */
+
+ short sut; /* Uphole time at source in milliseconds.
+ byte# 95-96
+ */
+
+ short gut; /* Uphole time at receiver group in milliseconds.
+ byte# 97-98
+ */
+
+ short sstat; /* Source static correction in milliseconds.
+ byte# 99-100
+ */
+
+ short gstat; /* Group static correction in milliseconds.
+ byte# 101-102
+ */
+
+ short tstat; /* Total static applied in milliseconds.
+ (Zero if no static has been applied.)
+ byte# 103-104
+ */
+
+ short laga; /* Lag time A, time in ms between end of 240-
+ byte trace identification header and time
+ break, positive if time break occurs after
+ end of header, time break is defined as
+ the initiation pulse which maybe recorded
+ on an auxiliary trace or as otherwise
+ specified by the recording system
+ byte# 105-106
+ */
+
+ short lagb; /* lag time B, time in ms between the time break
+ and the initiation time of the energy source,
+ may be positive or negative
+ byte# 107-108
+ */
+
+ short delrt; /* delay recording time, time in ms between
+ initiation time of energy source and time
+ when recording of data samples begins
+ (for deep water work if recording does not
+ start at zero time)
+ byte# 109-110
+ */
+
+ short muts; /* mute time--start
+ byte# 111-112
+ */
+
+ short mute; /* mute time--end
+ byte# 113-114
+ */
+
+ unsigned short ns; /* number of samples in this trace
+ byte# 115-116
+ */
+
+ unsigned short dt; /* sample interval; in micro-seconds
+ byte# 117-118
+ */
+
+ short gain; /* gain type of field instruments code:
+ 1 = fixed
+ 2 = binary
+ 3 = floating point
+ 4 ---- N = optional use
+ byte# 119-120
+ */
+
+ short igc; /* instrument gain constant
+ byte# 121-122
+ */
+
+ short igi; /* instrument early or initial gain
+ byte# 123-124
+ */
+
+ short corr; /* correlated:
+ 1 = no
+ 2 = yes
+ byte# 125-126
+ */
+
+ short sfs; /* sweep frequency at start
+ byte# 127-128
+ */
+
+ short sfe; /* sweep frequency at end
+ byte# 129-130
+ */
+
+ short slen; /* sweep length in ms
+ byte# 131-132
+ */
+
+ short styp; /* sweep type code:
+ 1 = linear
+ 2 = cos-squared
+ 3 = other
+ byte# 133-134
+ */
+
+ short stas; /* sweep trace length at start in ms
+ byte# 135-136
+ */
+
+ short stae; /* sweep trace length at end in ms
+ byte# 137-138
+ */
+
+ short tatyp; /* taper type: 1=linear, 2=cos^2, 3=other
+ byte# 139-140
+ */
+
+ short afilf; /* alias filter frequency if used
+ byte# 141-142
+ */
+
+ short afils; /* alias filter slope
+ byte# 143-144
+ */
+
+ short nofilf; /* notch filter frequency if used
+ byte# 145-146
+ */
+
+ short nofils; /* notch filter slope
+ byte# 147-148
+ */
+
+ short lcf; /* low cut frequency if used
+ byte# 149-150
+ */
+
+ short hcf; /* high cut frequncy if used
+ byte# 151-152
+ */
+
+ short lcs; /* low cut slope
+ byte# 153-154
+ */
+
+ short hcs; /* high cut slope
+ byte# 155-156
+ */
+
+ short year; /* year data recorded
+ byte# 157-158
+ */
+
+ short day; /* day of year
+ byte# 159-160
+ */
+
+ short hour; /* hour of day (24 hour clock)
+ byte# 161-162
+ */
+
+ short minute; /* minute of hour
+ byte# 163-164
+ */
+
+ short sec; /* second of minute
+ byte# 165-166
+ */
+
+ short timbas; /* time basis code:
+ 1 = local
+ 2 = GMT
+ 3 = other
+ byte# 167-168
+ */
+
+ short trwf; /* trace weighting factor, defined as 1/2^N
+ volts for the least sigificant bit
+ byte# 169-170
+ */
+
+ short grnors; /* geophone group number of roll switch
+ position one
+ byte# 171-172
+ */
+
+ short grnofr; /* geophone group number of trace one within
+ original field record
+ byte# 173-174
+ */
+
+ short grnlof; /* geophone group number of last trace within
+ original field record
+ byte# 175-176
+ */
+
+ short gaps; /* gap size (total number of groups dropped)
+ byte# 177-178
+ */
+
+ short otrav; /* overtravel taper code:
+ 1 = down (or behind)
+ 2 = up (or ahead)
+ byte# 179-180
+ */
+
+#ifdef SLTSU_SEGY_H /* begin Unocal SU segy.h differences */
+
+
+ /* cwp local assignments */
+ float d1; /* sample spacing for non-seismic data
+ byte# 181-184
+ */
+
+ float f1; /* first sample location for non-seismic data
+ byte# 185-188
+ */
+
+ float d2; /* sample spacing between traces
+ byte# 189-192
+ */
+
+ float f2; /* first trace location
+ byte# 193-196
+ */
+
+ float ungpow; /* negative of power used for dynamic
+ range compression
+ byte# 197-200
+ */
+
+ float unscale; /* reciprocal of scaling factor to normalize
+ range
+ byte# 201-204
+ */
+
+ short mark; /* mark selected traces
+ byte# 205-206
+ */
+
+ /* SLTSU local assignments */
+ short mutb; /* mute time at bottom (start time)
+ bottom mute ends at last sample
+ byte# 207-208
+ */
+ float dz; /* depth sampling interval in (m or ft)
+ if =0.0, input are time samples
+ byte# 209-212
+ */
+
+ float fz; /* depth of first sample in (m or ft)
+ byte# 213-116
+ */
+
+ short n2; /* number of traces per cdp or per shot
+ byte# 217-218
+ */
+
+ short shortpad; /* alignment padding
+ byte# 219-220
+ */
+
+ int ntr; /* number of traces
+ byte# 221-224
+ */
+
+ /* SLTSU local assignments end */
+
+ short unass[8]; /* unassigned
+ byte# 225-240
+ */
+
+#else
+
+ /* cwp local assignments */
+ float d1; /* sample spacing for non-seismic data
+ byte# 181-184
+ */
+
+ float f1; /* first sample location for non-seismic data
+ byte# 185-188
+ */
+
+ float d2; /* sample spacing between traces
+ byte# 189-192
+ */
+
+ float f2; /* first trace location
+ byte# 193-196
+ */
+
+ float ungpow; /* negative of power used for dynamic
+ range compression
+ byte# 197-200
+ */
+
+ float unscale; /* reciprocal of scaling factor to normalize
+ range
+ byte# 201-204
+ */
+
+ int ntr; /* number of traces
+ byte# 205-208
+ */
+
+ short mark; /* mark selected traces
+ byte# 209-210
+ */
+
+ short shortpad; /* alignment padding
+ byte# 211-212
+ */
+
+
+ short unass[14]; /* unassigned--NOTE: last entry causes
+ a break in the word alignment, if we REALLY
+ want to maintain 240 bytes, the following
+ entry should be an odd number of short/UINT2
+ OR do the insertion above the "mark" keyword
+ entry
+ byte# 213-240
+ */
+#endif
+
+} segy;
+
+
+typedef struct { /* bhed - binary header */
+
+ int jobid; /* job identification number */
+
+ int lino; /* line number (only one line per reel) */
+
+ int reno; /* reel number */
+
+ short ntrpr; /* number of data traces per record */
+
+ short nart; /* number of auxiliary traces per record */
+
+ unsigned short hdt; /* sample interval in micro secs for this reel */
+
+ unsigned short dto; /* same for original field recording */
+
+ unsigned short hns; /* number of samples per trace for this reel */
+
+ unsigned short nso; /* same for original field recording */
+
+ short format; /* data sample format code:
+ 1 = floating point, 4 byte (32 bits)
+ 2 = fixed point, 4 byte (32 bits)
+ 3 = fixed point, 2 byte (16 bits)
+ 4 = fixed point w/gain code, 4 byte (32 bits)
+ 5 = IEEE floating point, 4 byte (32 bits)
+ 8 = two's complement integer, 1 byte (8 bits)
+ */
+
+ short fold; /* CDP fold expected per CDP ensemble */
+
+ short tsort; /* trace sorting code:
+ 1 = as recorded (no sorting)
+ 2 = CDP ensemble
+ 3 = single fold continuous profile
+ 4 = horizontally stacked */
+
+ short vscode; /* vertical sum code:
+ 1 = no sum
+ 2 = two sum ...
+ N = N sum (N = 32,767) */
+
+ short hsfs; /* sweep frequency at start */
+
+ short hsfe; /* sweep frequency at end */
+
+ short hslen; /* sweep length (ms) */
+
+ short hstyp; /* sweep type code:
+ 1 = linear
+ 2 = parabolic
+ 3 = exponential
+ 4 = other */
+
+ short schn; /* trace number of sweep channel */
+
+ short hstas; /* sweep trace taper length at start if
+ tapered (the taper starts at zero time
+ and is effective for this length) */
+
+ short hstae; /* sweep trace taper length at end (the ending
+ taper starts at sweep length minus the taper
+ length at end) */
+
+ short htatyp; /* sweep trace taper type code:
+ 1 = linear
+ 2 = cos-squared
+ 3 = other */
+
+ short hcorr; /* correlated data traces code:
+ 1 = no
+ 2 = yes */
+
+ short bgrcv; /* binary gain recovered code:
+ 1 = yes
+ 2 = no */
+
+ short rcvm; /* amplitude recovery method code:
+ 1 = none
+ 2 = spherical divergence
+ 3 = AGC
+ 4 = other */
+
+ short mfeet; /* measurement system code:
+ 1 = meters
+ 2 = feet */
+
+ short polyt; /* impulse signal polarity code:
+ 1 = increase in pressure or upward
+ geophone case movement gives
+ negative number on tape
+ 2 = increase in pressure or upward
+ geophone case movement gives
+ positive number on tape */
+
+ short vpol; /* vibratory polarity code:
+ code seismic signal lags pilot by
+ 1 337.5 to 22.5 degrees
+ 2 22.5 to 67.5 degrees
+ 3 67.5 to 112.5 degrees
+ 4 112.5 to 157.5 degrees
+ 5 157.5 to 202.5 degrees
+ 6 202.5 to 247.5 degrees
+ 7 247.5 to 292.5 degrees
+ 8 293.5 to 337.5 degrees */
+
+ short hunass[170]; /* unassigned */
+
+} bhed;
+
+/* DEFINES */
+#define gettr(x) fgettr(stdin, (x))
+#define vgettr(x) fvgettr(stdin, (x))
+#define puttr(x) fputtr(stdout, (x))
+#define vputtr(x) fvputtr(stdout, (x))
+#define gettra(x, y) fgettra(stdin, (x), (y))
+
+
+/* TOTHER represents "other" */
+#define TOTHER -1
+/* TUNK represents time traces of an unknown type */
+#define TUNK 0
+/* TREAL represents real time traces */
+#define TREAL 1
+/* TDEAD represents dead time traces */
+#define TDEAD 2
+/* TDUMMY represents dummy time traces */
+#define TDUMMY 3
+/* TBREAK represents time break traces */
+#define TBREAK 4
+/* UPHOLE represents uphole traces */
+#define UPHOLE 5
+/* SWEEP represents sweep traces */
+#define SWEEP 6
+/* TIMING represents timing traces */
+#define TIMING 7
+/* WBREAK represents timing traces */
+#define WBREAK 8
+/* NFGUNSIG represents near field gun signature */
+#define NFGUNSIG 9
+/* FFGUNSIG represents far field gun signature */
+#define FFGUNSIG 10
+/* SPSENSOR represents seismic pressure sensor */
+#define SPSENSOR 11
+/* TVERT represents multicomponent seismic sensor
+ - vertical component */
+#define TVERT 12
+/* TXLIN represents multicomponent seismic sensor
+ - cross-line component */
+#define TXLIN 13
+/* TINLIN represents multicomponent seismic sensor
+ - in-line component */
+#define TINLIN 14
+/* ROTVERT represents rotated multicomponent seismic sensor
+ - vertical component */
+#define ROTVERT 15
+/* TTRANS represents rotated multicomponent seismic sensor
+ - transverse component */
+#define TTRANS 16
+/* TRADIAL represents rotated multicomponent seismic sensor
+ - radial component */
+#define TRADIAL 17
+/* VRMASS represents vibrator reaction mass */
+#define VRMASS 18
+/* VBASS represents vibrator baseplate */
+#define VBASS 19
+/* VEGF represents vibrator estimated ground force */
+#define VEGF 20
+/* VREF represents vibrator reference */
+#define VREF 21
+
+/*** CWP trid assignments ***/
+/* ACOR represents autocorrelation */
+#define ACOR 109
+/* FCMPLX represents fourier transformed - no packing
+ xr[0],xi[0], ..., xr[N-1],xi[N-1] */
+#define FCMPLX 110
+/* FUNPACKNYQ represents fourier transformed - unpacked Nyquist
+ xr[0],xi[0],...,xr[N/2],xi[N/2] */
+#define FUNPACKNYQ 111
+/* FTPACK represents fourier transformed - packed Nyquist
+ even N: xr[0],xr[N/2],xr[1],xi[1], ...,
+ xr[N/2 -1],xi[N/2 -1]
+ (note the exceptional second entry)
+ odd N:
+ xr[0],xr[(N-1)/2],xr[1],xi[1], ...,
+ xr[(N-1)/2 -1],xi[(N-1)/2 -1],xi[(N-1)/2]
+ (note the exceptional second & last entries)
+*/
+#define FTPACK 112
+/* TCMPLX represents complex time traces */
+#define TCMPLX 113
+/* FAMPH represents freq domain data in amplitude/phase form */
+#define FAMPH 114
+/* TAMPH represents time domain data in amplitude/phase form */
+#define TAMPH 115
+/* REALPART represents the real part of a trace to Nyquist */
+#define REALPART 116
+/* IMAGPART represents the real part of a trace to Nyquist */
+#define IMAGPART 117
+/* AMPLITUDE represents the amplitude of a trace to Nyquist */
+#define AMPLITUDE 118
+/* PHASE represents the phase of a trace to Nyquist */
+#define PHASE 119
+/* KT represents wavenumber-time domain data */
+#define KT 121
+/* KOMEGA represents wavenumber-frequency domain data */
+#define KOMEGA 122
+/* ENVELOPE represents the envelope of the complex time trace */
+#define ENVELOPE 123
+/* INSTPHASE represents the phase of the complex time trace */
+#define INSTPHASE 124
+/* INSTFREQ represents the frequency of the complex time trace */
+#define INSTFREQ 125
+/* DEPTH represents traces in depth-range (z-x) */
+#define TRID_DEPTH 130
+/* 3C data... v,h1,h2=(11,12,13)+32 so a bitmask will convert */
+/* between conventions */
+/* CHARPACK represents byte packed seismic data from supack1 */
+#define CHARPACK 201
+/* SHORTPACK represents 2 byte packed seismic data from supack2 */
+#define SHORTPACK 202
+
+
+#define ISSEISMIC(id) (( (id)==TUNK || (id)==TREAL || (id)==TDEAD || (id)==TDUMMY || (id)==TBREAK || (id)==UPHOLE || (id)==SWEEP || (id)==TIMING || (id)==WBREAK || (id)==NFGUNSIG || (id)==FFGUNSIG || (id)==SPSENSOR || (id)==TVERT || (id)==TXLIN || (id)==TINLIN || (id)==ROTVERT || (id)==TTRANS || (id)==TRADIAL || (id)==ACOR ) ? cwp_true : cwp_false )
+
+/* FUNCTION PROTOTYPES */
+#ifdef __cplusplus /* if C++, specify external linkage to C functions */
+extern "C" {
+#endif
+
+/* get trace and put trace */
+int fgettr(FILE *fp, segy *tp);
+int fvgettr(FILE *fp, segy *tp);
+void fputtr(FILE *fp, segy *tp);
+void fvputtr(FILE *fp, segy *tp);
+int fgettra(FILE *fp, segy *tp, int itr);
+
+/* get gather and put gather */
+segy **fget_gather(FILE *fp, cwp_String *key,cwp_String *type,Value *n_val,
+ int *nt,int *ntr, float *dt,int *first);
+segy **get_gather(cwp_String *key, cwp_String *type, Value *n_val,
+ int *nt, int *ntr, float *dt, int *first);
+segy **fput_gather(FILE *fp, segy **rec,int *nt, int *ntr);
+segy **put_gather(segy **rec,int *nt, int *ntr);
+
+/* hdrpkge */
+void gethval(const segy *tp, int index, Value *valp);
+void puthval(segy *tp, int index, Value *valp);
+void getbhval(const bhed *bhp, int index, Value *valp);
+void putbhval(bhed *bhp, int index, Value *valp);
+void gethdval(const segy *tp, char *key, Value *valp);
+void puthdval(segy *tp, char *key, Value *valp);
+char *hdtype(const char *key);
+char *getkey(const int index);
+int getindex(const char *key);
+void swaphval(segy *tp, int index);
+void swapbhval(bhed *bhp, int index);
+void printheader(const segy *tp);
+
+void tabplot(segy *tp, int itmin, int itmax);
+
+#ifdef __cplusplus /* if C++, end external linkage specification */
+}
+#endif
+
+#endif
diff --git a/raytime3d/src3d.c b/raytime3d/src3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..545c1974d14ab60f4bd28b259e68d2d4b4ee0e21
--- /dev/null
+++ b/raytime3d/src3d.c
@@ -0,0 +1,211 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
+#include <fcntl.h>
+
+#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
+#define s0(x,y,z) slow0[nxy*(z) + nx*(y) + (x)]
+#define SQR(x) ((x) * (x))
+#define DIST(x,y,z,x1,y1,z1) sqrt(SQR(x-(x1))+SQR(y-(y1)) + SQR(z-(z1)))
+
+/* definitions from verbose.c */
+extern void verr(char *fmt, ...);
+extern void vwarn(char *fmt, ...);
+extern void vmess(char *fmt, ...);
+
+void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube)
+{
+ int
+ srctype=1, /* if 1, source is a point;
+ 2, source is on the walls of the data volume;
+ 3, source on wall, time field known; */
+ srcwall, /* if 1, source on x=0 wall, if 2, on x=nx-1 wall
+ if 3, source on y=0 wall, if 4, on y=ny-1 wall
+ if 5, source on z=0 wall, if 6, on z=nz-1 wall */
+ xs, /* shot x position (in grid points) */
+ ys, /* shot y position */
+ zs, /* shot depth */
+ xx, yy, zz, /* Used to loop around xs, ys, zs coordinates */
+ ii, i, j, k,
+ wfint, ofint,
+ nxy, nyz, nxz, nxyz, nwall,
+ NCUBE=2;
+ float
+ fxs, /* shot position in X (in real units)*/
+ fys, /* shot position in Y (in real units)*/
+ fzs, /* shot position in Z (in real units)*/
+ *wall,
+ /* maximum offset (real units) to compute */
+ /* used in linear velocity gradient cube source */
+ rx, ry, rz, dvz, dv, v0,
+ rzc, rxyc, rz1, rxy1, rho, theta1, theta2,
+ xsrc1, ysrc1, zsrc1;
+ char
+ *oldtfile, /* file through which old travel times are input */
+ *wallfile; /* file containing input wall values of traveltimes */
+
+
+ if(!getparint("NCUBE",&NCUBE)) NCUBE=2;
+
+ if(!getparint("srctype",&srctype)) srctype=1;
+ if(srctype==1) {
+ if(!getparfloat("xsrc1",&xsrc1)) verr("xsrc1 not given");
+ if(!getparfloat("ysrc1",&ysrc1)) verr("ysrc1 not given");
+ if(!getparfloat("zsrc1",&zsrc1)) verr("zsrc1 not given");
+ fxs = (xsrc1-ox)/h;
+ fys = (ysrc1-oy)/h;
+ fzs = (zsrc1-oz)/h;
+ xs = (int)(fxs + 0.5);
+ ys = (int)(fys + 0.5);
+ zs = (int)(fzs + 0.5);
+ if(xs<2 || ys<2 || zs<2 || xs>nx-3 || ys>ny-3 || zs>nz-3){
+ vwarn("Source near an edge, beware of traveltime errors");
+ vwarn("for raypaths that travel parallel to edge ");
+ vwarn("while wavefronts are strongly curved, (JV, 8/17/88)\n");
+ }
+ *pxs = xs; *pys = ys, *pzs = zs, *cube = NCUBE;
+ }
+ else if (srctype==2) {
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
+ if (!getparstring("wallfile",&wallfile)) verr("wallfile not given");
+ if((wfint=open(wallfile,O_RDONLY,0664))<=1) {
+ fprintf(stderr,"cannot open %s\n",wallfile);
+ exit(-1);
+ }
+ }
+ else if (srctype==3) {
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
+ if (!getparstring("oldtfile",&oldtfile)) verr("oldtfile not given");
+ if((ofint=open(oldtfile,O_RDONLY,0664))<=1) {
+ fprintf(stderr,"cannot open %s\n",oldtfile);
+ exit(-1);
+ }
+ }
+ else {
+ verr("ERROR: incorrect value of srctype");
+ }
+
+ nxy = nx * ny;
+ nyz = ny * nz;
+ nxz = nx * nz;
+ nxyz = nx * ny * nz;
+
+
+ /* SET TIMES TO DUMMY VALUE */
+ for(i=0;i<nxyz;i++) time0[i] = 1.0e10;
+
+ if (srctype == 1) { /* VIDALE'S POINT SOURCE */
+ /* FILL IN CUBE AROUND SOURCE POINT */
+ /* HOLE'S NEW LINEAR VELOCITY GRADIENT CUBE (APRIL 1991)*/
+ v0 = h/s0(xs,ys,zs);
+ for (xx = xs-NCUBE; xx <= xs+NCUBE; xx++) {
+ if (xx < 0 || xx >= nx) continue;
+ for (yy = ys-NCUBE; yy <= ys+NCUBE; yy++) {
+ if (yy < 0 || yy >= ny) continue;
+ for (zz = zs-NCUBE; zz <= zs+NCUBE; zz++) {
+ if (zz < 0 || zz >= nz) continue;
+ if (zz == zs)
+ dvz = 1/s0(xx,yy,zz+1)-1/s0(xs,ys,zs);
+ else
+ dvz = (1/s0(xx,yy,zz)-1/s0(xs,ys,zs))/(zz-zs);
+ dv = fabs(dvz);
+ if (dv == 0.) {
+ t0(xx,yy,zz) = s0(xs,ys,zs)*DIST(fxs,fys,fzs,xx,yy,zz);
+ continue;
+ }
+ rzc = -v0/dv;
+ rx = h*(xx - fxs);
+ ry = h*(yy - fys);
+ rz = h*(zz - fzs);
+ rz1 = rz*dvz/dv;
+ rxy1 = sqrt(rx*rx+ry*ry+rz*rz-rz1*rz1);
+ if (rxy1<=h/1.e6)
+ t0(xx,yy,zz) = fabs(log((v0+dv*rz1)/v0)/dv);
+ else {
+ rxyc = (rz1*rz1+rxy1*rxy1-2*rz1*rzc)/(2*rxy1);
+ rho = sqrt(rzc*rzc+rxyc*rxyc);
+ theta1 = asin(-rzc/rho);
+ /* can't handle asin(1.) ! */
+ if (fabs(rz1-rzc)>=rho) rho=1.0000001*fabs(rz1-rzc);
+ theta2 = asin((rz1-rzc)/rho);
+ if (rxyc<0) theta1=M_PI-theta1;
+ if (rxyc<rxy1) theta2=M_PI-theta2;
+ t0(xx,yy,zz) = log(tan(theta2/2)/tan(theta1/2)) / dv;
+ }
+ }
+ }
+ }
+ }
+ else if (srctype == 2) { /* HOLE'S EXTERNAL SOURCE */
+
+ /* FILL IN WALLS' TIMES FROM EXTERNAL DATAFILE */
+ read (wfint,wall,4*nwall); /* READ X=0 WALL */
+ if (wall[0]>-1.e-20) {
+ ii = 0;
+ for (k=0; k<nz; k++) {
+ for (j=0; j<ny; j++) {
+ t0(0,j,k) = wall[ii];
+ ii++;
+ }
+ }
+ }
+ read (wfint,wall,4*nwall); /* READ X=NX-1 WALL */
+ if (wall[0]>-1.e-20) {
+ ii = 0;
+ for (k=0; k<nz; k++) {
+ for (j=0; j<ny; j++) {
+ t0(nx-1,j,k) = wall[ii];
+ ii++;
+ }
+ }
+ }
+ read (wfint,wall,4*nwall); /* READ Y=0 WALL */
+ if (wall[0]>-1.e-20) {
+ ii = 0;
+ for (k=0; k<nz; k++) {
+ for (i=0; i<nx; i++) {
+ t0(i,0,k) = wall[ii];
+ ii++;
+ }
+ }
+ }
+ read (wfint,wall,4*nwall); /* READ Y=NY-1 WALL */
+ if (wall[0]>-1.e-20) {
+ ii = 0;
+ for (k=0; k<nz; k++) {
+ for (i=0; i<nx; i++) {
+ t0(i,ny-1,k) = wall[ii];
+ ii++;
+ }
+ }
+ }
+ read (wfint,wall,4*nwall); /* READ Z=0 WALL */
+ if (wall[0]>-1.e-20) {
+ ii = 0;
+ for (j=0; j<ny; j++) {
+ for (i=0; i<nx; i++) {
+ t0(i,j,0) = wall[ii];
+ ii++;
+ }
+ }
+ }
+ read (wfint,wall,4*nwall); /* READ Z=NZ-1 WALL */
+ if (wall[0]>-1.e-20) {
+ ii = 0;
+ for (j=0; j<ny; j++) {
+ for (i=0; i<nx; i++) {
+ t0(i,j,nz-1) = wall[ii];
+ ii++;
+ }
+ }
+ }
+ }
+ else if (srctype == 3) { /* HOLE'S REDO OLD TIMES */
+ /* READ IN OLD TIME FILE */
+ if (srctype == 3) read(ofint,time0,nxyz*4);
+ }
+
+ return;
+}
diff --git a/raytime3d/threadAffinity.c b/raytime3d/threadAffinity.c
new file mode 100644
index 0000000000000000000000000000000000000000..49ca7e9d45bb953c9e63601c217d2deef69afcd1
--- /dev/null
+++ b/raytime3d/threadAffinity.c
@@ -0,0 +1,109 @@
+#define _GNU_SOURCE
+
+#include <stdio.h>
+#include <unistd.h>
+#include <string.h>
+#ifdef __USE_GNU
+#include <omp.h>
+#include <sched.h>
+#else /* for OSX */
+#include <sched.h>
+#include <sys/types.h>
+#include <sys/sysctl.h>
+
+#define CPU_SETSIZE 1024
+#define SYSCTL_CORE_COUNT "machdep.cpu.core_count"
+void vmess(char *fmt, ...);
+
+typedef struct cpu_set {
+ uint32_t count;
+} cpu_set_t;
+
+static inline void
+CPU_ZERO(cpu_set_t *cs) { cs->count = 0; }
+
+static inline void
+CPU_SET(int num, cpu_set_t *cs) { cs->count |= (1 << num); }
+
+static inline int
+CPU_ISSET(int num, cpu_set_t *cs) { return (cs->count & (1 << num)); }
+
+int sched_getaffinity(pid_t pid, size_t cpu_size, cpu_set_t *cpu_set)
+{
+ int32_t core_count = 0;
+ size_t len = sizeof(core_count);
+ int ret = sysctlbyname(SYSCTL_CORE_COUNT, &core_count, &len, 0, 0);
+ if (ret) {
+ printf("error while get core count %d\n", ret);
+ return -1;
+ }
+ cpu_set->count = 0;
+ for (int i = 0; i < core_count; i++) {
+ cpu_set->count |= (1 << i);
+ }
+
+ return 0;
+}
+#endif
+
+/* Borrowed from util-linux-2.13-pre7/schedutils/taskset.c */
+
+static char *cpuset_to_cstr(cpu_set_t *mask, char *str)
+{
+ char *ptr = str;
+ int i, j, entry_made = 0;
+ for (i = 0; i < CPU_SETSIZE; i++) {
+ if (CPU_ISSET(i, mask)) {
+ int run = 0;
+ entry_made = 1;
+ for (j = i + 1; j < CPU_SETSIZE; j++) {
+ if (CPU_ISSET(j, mask)) run++;
+ else break;
+ }
+ if (!run)
+ sprintf(ptr, "%d,", i);
+ else if (run == 1) {
+ sprintf(ptr, "%d,%d,", i, i + 1);
+ i++;
+ } else {
+ sprintf(ptr, "%d-%d,", i, i + run);
+ i += run;
+ }
+ while (*ptr != 0) ptr++;
+ }
+ }
+ ptr -= entry_made;
+ *ptr = 0;
+ return(str);
+}
+
+void threadAffinity(void)
+{
+ int thread;
+ cpu_set_t coremask;
+ char clbuf[7 * CPU_SETSIZE], hnbuf[64];
+ char prefix[200];
+
+ memset(clbuf, 0, sizeof(clbuf));
+ memset(hnbuf, 0, sizeof(hnbuf));
+ (void)gethostname(hnbuf, sizeof(hnbuf));
+
+ strcpy(prefix,"Hello world from");
+
+// #pragma omp parallel private(thread, coremask, clbuf)
+/* for use inside parallel region */
+ #pragma omp critical
+ {
+#ifdef __USE_GNU
+ thread = omp_get_thread_num();
+#else
+ thread = 1;
+#endif
+ (void)sched_getaffinity(0, sizeof(coremask), &coremask);
+ cpuset_to_cstr(&coremask, clbuf);
+ vmess("%s thread %d, on %s. (core affinity = %s)", prefix, thread, hnbuf, clbuf);
+
+ }
+ return;
+}
+
diff --git a/raytime3d/verbosepkg.c b/raytime3d/verbosepkg.c
new file mode 100644
index 0000000000000000000000000000000000000000..483e5f92bd5e1c1a495c66b7a63b9e8113943897
--- /dev/null
+++ b/raytime3d/verbosepkg.c
@@ -0,0 +1,77 @@
+#include <stdio.h>
+#include <stdarg.h>
+#include "par.h"
+#include <string.h>
+#ifdef _CRAYMPP
+#include <intrinsics.h>
+#endif
+
+/**
+* functions to print out verbose, error and warning messages to stderr.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+
+void verr(char *fmt, ...)
+{
+ va_list args;
+
+ if (EOF == fflush(stderr)) {
+ fprintf(stderr, "\nverr: fflush failed on stderr");
+ }
+ fprintf(stderr, " Error in %s: ", xargv[0]);
+#ifdef _CRAYMPP
+ fprintf(stderr, "PE %d: ", _my_pe());
+#elif defined(SGI)
+ fprintf(stderr, "PE %d: ", mp_my_threadnum());
+#endif
+ va_start(args,fmt);
+ vfprintf(stderr, fmt, args);
+ va_end(args);
+ fprintf(stderr, "\n");
+
+ exit(EXIT_FAILURE);
+}
+
+void vwarn(char *fmt, ...)
+{
+ va_list args;
+
+ if (EOF == fflush(stderr)) {
+ fprintf(stderr, "\nvwarn: fflush failed on stderr");
+ }
+ fprintf(stderr, " Warning in %s: ", xargv[0]);
+#ifdef _CRAYMPP
+ fprintf(stderr, "PE %d: ", _my_pe());
+#elif defined(SGI)
+ fprintf(stderr, "PE %d: ", mp_my_threadnum());
+#endif
+ va_start(args,fmt);
+ vfprintf(stderr, fmt, args);
+ va_end(args);
+ fprintf(stderr, "\n");
+ return;
+}
+
+void vmess(char *fmt, ...)
+{
+ va_list args;
+
+ if (EOF == fflush(stderr)) {
+ fprintf(stderr, "\nvmess: fflush failed on stderr");
+ }
+ fprintf(stderr, " %s: ", xargv[0]);
+#ifdef _CRAYMPP
+ fprintf(stderr, "PE %d: ", _my_pe());
+#elif defined(SGI)
+ fprintf(stderr, "PE %d: ", mp_my_threadnum());
+#endif
+ va_start(args,fmt);
+ vfprintf(stderr, fmt, args);
+ va_end(args);
+ fprintf(stderr, "\n");
+ return;
+}
diff --git a/raytime3d/vidale3d.c b/raytime3d/vidale3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..ff8a74c1c4f75f4179ce168c514b6ee1949ab837
--- /dev/null
+++ b/raytime3d/vidale3d.c
@@ -0,0 +1,1444 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
+#include"par.h"
+
+#define SQR2 1.414213562
+#define SQR3 1.732050808
+#define SQR6 2.449489743
+#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
+#define s0(x,y,z) slow0[nxy*(z) + nx*(y) + (x)]
+
+/* definitions from verbose.c */
+extern void verr(char *fmt, ...);
+extern void vwarn(char *fmt, ...);
+extern void vmess(char *fmt, ...);
+
+struct sorted
+ { float time; int i1, i2;};
+
+int compar(struct sorted *a,struct sorted *b);
+
+float fdhne(float t1,float t2,float t3,float t4,float t5,float ss0,float s1,float s2,float s3);
+float fdh3d(float t1,float t2,float t3,float t4,float t5,float t6,float t7,float ss0,float s1,float s2,float s3,float s4,float s5,float s6,float s7);
+float fdh2d(float t1,float t2,float t3,float ss0,float s1,float s2,float s3);
+float fdhnf(float t1,float t2,float t3,float t4,float t5,float ss0,float s1);
+
+void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int xs, int ys, int zs, int NCUBE)
+{
+ int
+ srctype=1, /* if 1, source is a point;
+ 2, source is on the walls of the data volume;
+ 3, source on wall, time field known; */
+ srcwall, /* if 1, source on x=0 wall, if 2, on x=nx-1 wall
+ if 3, source on y=0 wall, if 4, on y=ny-1 wall
+ if 5, source on z=0 wall, if 6, on z=nz-1 wall */
+ iplus=1, /* rate of expansion of "cube" in the */
+ iminus=1, /* plus/minus x/y/z direction */
+ jplus=1,
+ jminus=1,
+ kplus=1,
+ kminus=1,
+ igrow, /* counter for "cube" growth */
+ X1, X2, lasti, index, ii, i, j, k, radius,
+ nxy, nyz, nxz, nxyz, nwall,
+ /* counters for the position of the sides of current cube */
+ x1, x2, y1, y2, z1, z2,
+ /* flags set to 1 until a side has been reached */
+ dx1=1, dx2=1, dy1=1, dy2=1, dz1=1, dz2=1, rad0=1,
+ maxrad, /* maximum radius to compute */
+ reverse=1, /* will automatically do up to this number of
+ reverse propagation steps to fix waves that travel
+ back into expanding cube */
+ headpref=6, /* if headpref starts > 0, will determine
+ model wall closest to source and will prefer to start
+ reverse calculations on opposite wall */
+ /* counters for detecting head waves on sides of current cube */
+ head,headw[7], verbose;
+ float
+ *wall,
+ guess, try,
+ /* maximum offset (real units) to compute */
+ maxoff = -1.,
+ /* used to detect head waves: if headwave operator decreases
+ the previously-computed traveltime by at least
+ headtest*<~time_across_cell> then the headwave counter is
+ triggered */
+ fhead,headtest=1.e-3;
+
+ /* ARRAY TO ORDER SIDE FOR SOLUTION IN THE RIGHT ORDER */
+ struct sorted *sort;
+
+ if (!getparint("verbose",&verbose)) verbose=0;
+ if(!getparfloat("maxoff",&maxoff)) maxoff = -1.;
+ if(!getparint("iminus",&iminus)) iminus=1;
+ if(!getparint("iplus",&iplus)) iplus=1;
+ if(!getparint("jminus",&jminus)) jminus=1;
+ if(!getparint("jplus",&jplus)) jplus=1;
+ if(!getparint("kminus",&kminus)) kminus=1;
+ if(!getparint("kplus",&kplus)) kplus=1;
+ if(!getparint("reverse",&reverse)) reverse=0;
+ if(!getparint("headpref",&headpref)) headpref=6;
+ if(!getparint("NCUBE",&NCUBE)) NCUBE=2;
+
+ /* SET MAXIMUM RADIUS TO COMPUTE */
+ if (maxoff > 0.) {
+ maxrad = maxoff/h + 1;
+ vwarn("WARNING: Computing only to max radius = %d",maxrad);
+ }
+ else maxrad = 99999999;
+
+ nxy = nx * ny;
+ nyz = ny * nz;
+ nxz = nx * nz;
+ nxyz = nx * ny * nz;
+
+ /* MAKE ARRAY SORT LARGE ENOUGH FOR ANY SIDE */
+ if(nx <= ny && nx <= nz) {
+ sort = (struct sorted *) malloc(sizeof(struct sorted)*ny*nz);
+ nwall = nyz;
+ }
+ else if(ny <= nx && ny <= nz) {
+ sort = (struct sorted *) malloc(sizeof(struct sorted)*nx*nz);
+ nwall = nxz;
+ }
+ else {
+ sort = (struct sorted *) malloc(sizeof(struct sorted)*nx*ny);
+ nwall = nxy;
+ }
+ wall = (float *) malloc(4*nwall);
+ if(sort == NULL || wall == NULL)
+ verr("error in allocation of arrays sort and wall");
+
+ if(!getparint("srctype",&srctype)) srctype=1;
+ if(srctype==1) {
+ /* SETS LOCATION OF THE SIDES OF THE CUBE */
+ radius = NCUBE;
+ if(xs > NCUBE) x1 = xs - (NCUBE + 1);
+ else{ x1 = -1; dx1 = 0;}
+ if(xs < nx-(NCUBE + 1)) x2 = xs + (NCUBE + 1);
+ else{ x2 = nx; dx2 = 0;}
+ if(ys > NCUBE) y1 = ys - (NCUBE + 1);
+ else{ y1 = -1; dy1 = 0;}
+ if(ys < ny-(NCUBE + 1)) y2 = ys + (NCUBE + 1);
+ else{ y2 = ny; dy2 = 0;}
+ if(zs > NCUBE) z1 = zs - (NCUBE + 1);
+ else{ z1 = -1; dz1 = 0;}
+ if(zs < nz-(NCUBE + 1)) z2 = zs + (NCUBE + 1);
+ else{ z2 = nz; dz2 = 0;}
+ }
+ else {
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
+ /* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE */
+ radius = 1;
+ if (srcwall == 1) x2=1;
+ else { x2=nx; dx2=0; }
+ if (srcwall == 2) x1=nx-2;
+ else { x1= -1; dx1=0; }
+ if (srcwall == 3) y2=1;
+ else { y2=ny; dy2=0; }
+ if (srcwall == 4) y1=ny-2;
+ else { y1= -1; dy1=0; }
+ if (srcwall == 5) z2=1;
+ else { z2=nz; dz2=0; }
+ if (srcwall == 6) z1=nz-2;
+ else { z1= -1; dz1=0; }
+ }
+
+
+ if (headpref>0) { /* HOLE - PREFERRED REVERSE DIRECTION */
+ head = nx*ny*nz;
+ if (nx>5 && x2<=head) {headpref=2; head=x2;}
+ if (nx>5 && (nx-1-x1)<=head) {headpref=1; head=nx-1-x1;}
+ if (ny>5 && y2<=head) {headpref=4; head=y2;}
+ if (ny>5 && (ny-1-y1)<=head) {headpref=3; head=ny-1-y1;}
+ if (nz>5 && z2<=head) {headpref=6; head=z2;}
+ if (nz>5 && (nz-1-z1)<=head) {headpref=5; head=nz-1-z1;}
+ }
+
+ /* BIGGER LOOP - HOLE - ALLOWS AUTOMATIC REVERSE PROPAGATION IF
+ HEAD WAVES ARE ENCOUNTERED ON FACES OF EXPANDING CUBE,
+ ALLOWING WAVES TO TRAVEL BACK INTO THE CUBE */
+
+ while ( reverse > -1 ) {
+
+ headw[1]=0; headw[2]=0; headw[3]=0; headw[4]=0;
+ headw[5]=0; headw[6]=0;
+
+ /* BIG LOOP */
+ while(rad0 && (dx1 || dx2 || dy1 || dy2 || dz1 || dz2)) {
+
+ /* CALCULATE ON PRIMARY (time0) GRID */
+
+ /* TOP SIDE */
+ for (igrow=1;igrow<=kminus;igrow++) {
+ if(dz1){
+ ii = 0;
+ for(j=y1+1; j<=y2-1; j++){
+ for(i=x1+1; i<=x2-1; i++){
+ sort[ii].time = t0(i,j,z1+1);
+ sort[ii].i1 = i;
+ sort[ii].i2 = j;
+ ii++;
+ }
+ }
+ qsort((char *)sort,ii,sizeof(struct sorted),compar);
+ for(i=0;i<ii;i++){
+ X1 = sort[i].i1;
+ X2 = sort[i].i2;
+ index = z1*nxy + X2*nx + X1;
+ lasti = (z1+1)*nxy + X2*nx + X1;
+ fhead = 0.;
+ guess = time0[index];
+ if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,X2,z1+1),
+ t0(X1+1,X2,z1+1),t0(X1+1,X2+1,z1+1),t0(X1,X2+1,z1+1),
+ t0(X1+1,X2,z1 ),t0(X1+1,X2+1,z1 ),t0(X1,X2+1,z1 ),
+ s0(X1,X2,z1), s0(X1,X2,z1+1),
+ s0(X1+1,X2,z1+1),s0(X1+1,X2+1,z1+1),s0(X1,X2+1,z1+1),
+ s0(X1+1,X2,z1 ),s0(X1+1,X2+1,z1 ),s0(X1,X2+1,z1 ));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1>0 ) {
+ try = fdh3d( t0(X1,X2,z1+1),
+ t0(X1-1,X2,z1+1),t0(X1-1,X2+1,z1+1),t0(X1,X2+1,z1+1),
+ t0(X1-1,X2,z1 ),t0(X1-1,X2+1,z1 ),t0(X1,X2+1,z1 ),
+ s0(X1,X2,z1), s0(X1,X2,z1+1),
+ s0(X1-1,X2,z1+1),s0(X1-1,X2+1,z1+1),s0(X1,X2+1,z1+1),
+ s0(X1-1,X2,z1 ),s0(X1-1,X2+1,z1 ),s0(X1,X2+1,z1 ));
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,X2,z1+1),
+ t0(X1+1,X2,z1+1),t0(X1+1,X2-1,z1+1),t0(X1,X2-1,z1+1),
+ t0(X1+1,X2,z1 ),t0(X1+1,X2-1,z1 ),t0(X1,X2-1,z1 ),
+ s0(X1,X2,z1), s0(X1,X2,z1+1),
+ s0(X1+1,X2,z1+1),s0(X1+1,X2-1,z1+1),s0(X1,X2-1,z1+1),
+ s0(X1+1,X2,z1 ),s0(X1+1,X2-1,z1 ),s0(X1,X2-1,z1 ));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh3d( t0(X1,X2,z1+1),
+ t0(X1-1,X2,z1+1),t0(X1-1,X2-1,z1+1),t0(X1,X2-1,z1+1),
+ t0(X1-1,X2,z1 ),t0(X1-1,X2-1,z1 ),t0(X1,X2-1,z1 ),
+ s0(X1,X2,z1), s0(X1,X2,z1+1),
+ s0(X1-1,X2,z1+1),s0(X1-1,X2-1,z1+1),s0(X1,X2-1,z1+1),
+ s0(X1-1,X2,z1 ),s0(X1-1,X2-1,z1 ),s0(X1,X2-1,z1 ));
+ if (try<guess) guess = try;
+ }
+ if(guess > 1.0e9){
+ if(time0[index+1] < 1.e9 && X1<nx-1 && X2>y1+1 && X2<y2-1 ) {
+ try = fdhne(t0(X1,X2,z1+1),t0(X1+1,X2,z1+1),t0(X1+1,X2,z1),
+ t0(X1+1,X2-1,z1+1),t0(X1+1,X2+1,z1+1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1+1,X2,z1+1),s0(X1+1,X2,z1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 && X2>y1+1 && X2<y2-1 ) {
+ try = fdhne(t0(X1,X2,z1+1),t0(X1-1,X2,z1+1),t0(X1-1,X2,z1),
+ t0(X1-1,X2-1,z1+1),t0(X1-1,X2+1,z1+1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1-1,X2,z1+1),s0(X1-1,X2,z1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && X2<ny-1 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,X2,z1+1),t0(X1,X2+1,z1+1),t0(X1,X2+1,z1),
+ t0(X1-1,X2+1,z1+1),t0(X1+1,X2+1,z1+1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1,X2+1,z1+1),s0(X1,X2+1,z1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,X2,z1+1),t0(X1,X2-1,z1+1),t0(X1,X2-1,z1),
+ t0(X1-1,X2-1,z1+1),t0(X1+1,X2-1,z1+1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1,X2-1,z1+1),s0(X1,X2-1,z1) );
+ if (try<guess) guess = try;
+ }
+ }
+ if(time0[index+1] < 1.e9 && X1<nx-1 ) {
+ try = fdh2d(t0(X1,X2,z1+1),t0(X1+1,X2,z1+1),t0(X1+1,X2,z1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1+1,X2,z1+1),s0(X1+1,X2,z1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 ) {
+ try = fdh2d(t0(X1,X2,z1+1),t0(X1-1,X2,z1+1),t0(X1-1,X2,z1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1-1,X2,z1+1),s0(X1-1,X2,z1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && X2<ny-1 ) {
+ try = fdh2d(t0(X1,X2,z1+1),t0(X1,X2+1,z1+1),t0(X1,X2+1,z1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1,X2+1,z1+1),s0(X1,X2+1,z1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X2>0 ) {
+ try = fdh2d(t0(X1,X2,z1+1),t0(X1,X2-1,z1+1),t0(X1,X2-1,z1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1),s0(X1,X2-1,z1+1),s0(X1,X2-1,z1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,X2,z1),t0(X1+1,X2+1,z1),t0(X1,X2+1,z1),
+ s0(X1,X2,z1),
+ s0(X1+1,X2,z1),s0(X1+1,X2+1,z1),s0(X1,X2+1,z1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,X2,z1),t0(X1+1,X2-1,z1),t0(X1,X2-1,z1),
+ s0(X1,X2,z1),
+ s0(X1+1,X2,z1),s0(X1+1,X2-1,z1),s0(X1,X2-1,z1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1>0 ) {
+ try = fdh2d(t0(X1-1,X2,z1),t0(X1-1,X2+1,z1),t0(X1,X2+1,z1),
+ s0(X1,X2,z1),
+ s0(X1-1,X2,z1),s0(X1-1,X2+1,z1),s0(X1,X2+1,z1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh2d(t0(X1-1,X2,z1),t0(X1-1,X2-1,z1),t0(X1,X2-1,z1),
+ s0(X1,X2,z1),
+ s0(X1-1,X2,z1),s0(X1-1,X2-1,z1),s0(X1,X2-1,z1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(guess > 1.0e9){
+ if ( X1>x1+1 && X1<x2-1 && X2>y1+1 && X2<y2-1 ) {
+ try = fdhnf(t0(X1,X2,z1+1),
+ t0(X1+1,X2,z1+1),t0(X1,X2+1,z1+1),
+ t0(X1-1,X2,z1+1),t0(X1,X2-1,z1+1),
+ s0(X1,X2,z1),
+ s0(X1,X2,z1+1) );
+ if (try<guess) guess = try;
+ }
+ }
+ try = t0(X1,X2,z1+1) + .5*(s0(X1,X2,z1)+s0(X1,X2,z1+1));
+ if (try<guess) guess = try;
+ if ( time0[index+1]<1.e9 && X1<nx-1 ) {
+ try = t0(X1+1,X2,z1) + .5*(s0(X1,X2,z1)+s0(X1+1,X2,z1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-1]<1.e9 && X1>0 ) {
+ try = t0(X1-1,X2,z1) + .5*(s0(X1,X2,z1)+s0(X1-1,X2,z1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index+nx]<1.e9 && X2<ny-1 ) {
+ try = t0(X1,X2+1,z1) + .5*(s0(X1,X2,z1)+s0(X1,X2+1,z1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nx]<1.e9 && X2>0 ) {
+ try = t0(X1,X2-1,z1) + .5*(s0(X1,X2,z1)+s0(X1,X2-1,z1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if (guess<time0[index]) {
+ time0[index] = guess;
+ if (fhead>headtest) headw[5]++;
+ }
+ }
+ if(z1 == 0) dz1 = 0;
+ z1--;
+ }
+ }
+ /* BOTTOM SIDE */
+ for (igrow=1;igrow<=kplus;igrow++) {
+ if(dz2){
+ ii = 0;
+ for(j=y1+1; j<=y2-1; j++){
+ for(i=x1+1; i<=x2-1; i++){
+ sort[ii].time = t0(i,j,z2-1);
+ sort[ii].i1 = i;
+ sort[ii].i2 = j;
+ ii++;
+ }
+ }
+ qsort((char *)sort,ii,sizeof(struct sorted),compar);
+ for(i=0;i<ii;i++){
+ X1 = sort[i].i1;
+ X2 = sort[i].i2;
+ index = z2*nxy + X2*nx + X1;
+ lasti = (z2-1)*nxy + X2*nx + X1;
+ fhead = 0.;
+ guess = time0[index];
+ if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,X2,z2-1),
+ t0(X1+1,X2,z2-1),t0(X1+1,X2+1,z2-1),t0(X1,X2+1,z2-1),
+ t0(X1+1,X2,z2 ),t0(X1+1,X2+1,z2 ),t0(X1,X2+1,z2 ),
+ s0(X1,X2,z2), s0(X1,X2,z2-1),
+ s0(X1+1,X2,z2-1),s0(X1+1,X2+1,z2-1),s0(X1,X2+1,z2-1),
+ s0(X1+1,X2,z2 ),s0(X1+1,X2+1,z2 ),s0(X1,X2+1,z2 ));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1>0 ) {
+ try = fdh3d( t0(X1,X2,z2-1),
+ t0(X1-1,X2,z2-1),t0(X1-1,X2+1,z2-1),t0(X1,X2+1,z2-1),
+ t0(X1-1,X2,z2 ),t0(X1-1,X2+1,z2 ),t0(X1,X2+1,z2 ),
+ s0(X1,X2,z2), s0(X1,X2,z2-1),
+ s0(X1-1,X2,z2-1),s0(X1-1,X2+1,z2-1),s0(X1,X2+1,z2-1),
+ s0(X1-1,X2,z2 ),s0(X1-1,X2+1,z2 ),s0(X1,X2+1,z2 ));
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,X2,z2-1),
+ t0(X1+1,X2,z2-1),t0(X1+1,X2-1,z2-1),t0(X1,X2-1,z2-1),
+ t0(X1+1,X2,z2 ),t0(X1+1,X2-1,z2 ),t0(X1,X2-1,z2 ),
+ s0(X1,X2,z2), s0(X1,X2,z2-1),
+ s0(X1+1,X2,z2-1),s0(X1+1,X2-1,z2-1),s0(X1,X2-1,z2-1),
+ s0(X1+1,X2,z2 ),s0(X1+1,X2-1,z2 ),s0(X1,X2-1,z2 ));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh3d( t0(X1,X2,z2-1),
+ t0(X1-1,X2,z2-1),t0(X1-1,X2-1,z2-1),t0(X1,X2-1,z2-1),
+ t0(X1-1,X2,z2 ),t0(X1-1,X2-1,z2 ),t0(X1,X2-1,z2 ),
+ s0(X1,X2,z2), s0(X1,X2,z2-1),
+ s0(X1-1,X2,z2-1),s0(X1-1,X2-1,z2-1),s0(X1,X2-1,z2-1),
+ s0(X1-1,X2,z2 ),s0(X1-1,X2-1,z2 ),s0(X1,X2-1,z2 ));
+ if (try<guess) guess = try;
+ }
+ if(guess > 1.0e9){
+ if(time0[index+1] < 1.e9 && X1<nx-1 && X2>y1+1 && X2<y2-1 ) {
+ try = fdhne(t0(X1,X2,z2-1),t0(X1+1,X2,z2-1),t0(X1+1,X2,z2),
+ t0(X1+1,X2-1,z2-1),t0(X1+1,X2+1,z2-1),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1+1,X2,z2-1),s0(X1+1,X2,z2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 && X2>y1+1 && X2<y2-1 ) {
+ try = fdhne(t0(X1,X2,z2-1),t0(X1-1,X2,z2-1),t0(X1-1,X2,z2),
+ t0(X1-1,X2-1,z2-1),t0(X1-1,X2+1,z2-1),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1-1,X2,z2-1),s0(X1-1,X2,z2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && X2<ny-1 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,X2,z2-1),t0(X1,X2+1,z2-1),t0(X1,X2+1,z2),
+ t0(X1-1,X2+1,z2-1),t0(X1+1,X2+1,z2-1),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1,X2+1,z2-1),s0(X1,X2+1,z2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,X2,z2-1),t0(X1,X2-1,z2-1),t0(X1,X2-1,z2),
+ t0(X1-1,X2-1,z2-1),t0(X1+1,X2-1,z2-1),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1,X2-1,z2-1),s0(X1,X2-1,z2) );
+ if (try<guess) guess = try;
+ }
+ }
+ if(time0[index+1] < 1.e9 && X1<nx-1 ) {
+ try = fdh2d(t0(X1,X2,z2-1),t0(X1+1,X2,z2-1),t0(X1+1,X2,z2),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1+1,X2,z2-1),s0(X1+1,X2,z2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 ) {
+ try = fdh2d(t0(X1,X2,z2-1),t0(X1-1,X2,z2-1),t0(X1-1,X2,z2),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1-1,X2,z2-1),s0(X1-1,X2,z2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && X2<ny-1 ) {
+ try = fdh2d(t0(X1,X2,z2-1),t0(X1,X2+1,z2-1),t0(X1,X2+1,z2),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1,X2+1,z2-1),s0(X1,X2+1,z2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X2>0 ) {
+ try = fdh2d(t0(X1,X2,z2-1),t0(X1,X2-1,z2-1),t0(X1,X2-1,z2),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1),s0(X1,X2-1,z2-1),s0(X1,X2-1,z2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,X2,z2),t0(X1+1,X2+1,z2),t0(X1,X2+1,z2),
+ s0(X1,X2,z2),
+ s0(X1+1,X2,z2),s0(X1+1,X2+1,z2),s0(X1,X2+1,z2) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,X2,z2),t0(X1+1,X2-1,z2),t0(X1,X2-1,z2),
+ s0(X1,X2,z2),
+ s0(X1+1,X2,z2),s0(X1+1,X2-1,z2),s0(X1,X2-1,z2) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
+ && time0[index+nx] < 1.e9 && X2<ny-1 && X1>0 ) {
+ try = fdh2d(t0(X1-1,X2,z2),t0(X1-1,X2+1,z2),t0(X1,X2+1,z2),
+ s0(X1,X2,z2),
+ s0(X1-1,X2,z2),s0(X1-1,X2+1,z2),s0(X1,X2+1,z2) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
+ && time0[index-nx] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh2d(t0(X1-1,X2,z2),t0(X1-1,X2-1,z2),t0(X1,X2-1,z2),
+ s0(X1,X2,z2),
+ s0(X1-1,X2,z2),s0(X1-1,X2-1,z2),s0(X1,X2-1,z2) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(guess > 1.0e9){
+ if ( X1>x1+1 && X1<x2-1 && X2>y1+1 && X2<y2-1 ) {
+ try = fdhnf(t0(X1,X2,z2-1),
+ t0(X1+1,X2,z2-1),t0(X1,X2+1,z2-1),
+ t0(X1-1,X2,z2-1),t0(X1,X2-1,z2-1),
+ s0(X1,X2,z2),
+ s0(X1,X2,z2-1) );
+ if (try<guess) guess = try;
+ }
+ }
+ try = t0(X1,X2,z2-1) + .5*(s0(X1,X2,z2)+s0(X1,X2,z2-1));
+ if (try<guess) guess = try;
+ if ( time0[index+1]<1.e9 && X1<nx-1 ) {
+ try = t0(X1+1,X2,z2) + .5*(s0(X1,X2,z2)+s0(X1+1,X2,z2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-1]<1.e9 && X1>0 ) {
+ try = t0(X1-1,X2,z2) + .5*(s0(X1,X2,z2)+s0(X1-1,X2,z2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index+nx]<1.e9 && X2<ny-1 ) {
+ try = t0(X1,X2+1,z2) + .5*(s0(X1,X2,z2)+s0(X1,X2+1,z2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nx]<1.e9 && X2>0 ) {
+ try = t0(X1,X2-1,z2) + .5*(s0(X1,X2,z2)+s0(X1,X2-1,z2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if (guess<time0[index]) {
+ time0[index] = guess;
+ if (fhead>headtest) headw[6]++;
+ }
+ }
+ if(z2 == nz-1) dz2 = 0;
+ z2++;
+ }
+ }
+ /* FRONT SIDE */
+ for (igrow=1;igrow<=jminus;igrow++) {
+ if(dy1){
+ ii = 0;
+ for(k=z1+1; k<=z2-1; k++){
+ for(i=x1+1; i<=x2-1; i++){
+ sort[ii].time = t0(i,y1+1,k);
+ sort[ii].i1 = i;
+ sort[ii].i2 = k;
+ ii++;
+ }
+ }
+ qsort((char *)sort,ii,sizeof(struct sorted),compar);
+ for(i=0;i<ii;i++){
+ X1 = sort[i].i1;
+ X2 = sort[i].i2;
+ index = X2*nxy + y1*nx + X1;
+ lasti = X2*nxy + (y1+1)*nx + X1;
+ fhead = 0.;
+ guess = time0[index];
+ if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,y1+1,X2),
+ t0(X1+1,y1+1,X2),t0(X1+1,y1+1,X2+1),t0(X1,y1+1,X2+1),
+ t0(X1+1,y1 ,X2),t0(X1+1,y1 ,X2+1),t0(X1,y1 ,X2+1),
+ s0(X1,y1,X2), s0(X1,y1+1,X2),
+ s0(X1+1,y1+1,X2),s0(X1+1,y1+1,X2+1),s0(X1,y1+1,X2+1),
+ s0(X1+1,y1 ,X2),s0(X1+1,y1 ,X2+1),s0(X1,y1 ,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh3d( t0(X1,y1+1,X2),
+ t0(X1-1,y1+1,X2),t0(X1-1,y1+1,X2+1),t0(X1,y1+1,X2+1),
+ t0(X1-1,y1 ,X2),t0(X1-1,y1 ,X2+1),t0(X1,y1 ,X2+1),
+ s0(X1,y1,X2), s0(X1,y1+1,X2),
+ s0(X1-1,y1+1,X2),s0(X1-1,y1+1,X2+1),s0(X1,y1+1,X2+1),
+ s0(X1-1,y1 ,X2),s0(X1-1,y1 ,X2+1),s0(X1,y1 ,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,y1+1,X2),
+ t0(X1+1,y1+1,X2),t0(X1+1,y1+1,X2-1),t0(X1,y1+1,X2-1),
+ t0(X1+1,y1 ,X2),t0(X1+1,y1 ,X2-1),t0(X1,y1 ,X2-1),
+ s0(X1,y1,X2), s0(X1,y1+1,X2),
+ s0(X1+1,y1+1,X2),s0(X1+1,y1+1,X2-1),s0(X1,y1+1,X2-1),
+ s0(X1+1,y1 ,X2),s0(X1+1,y1 ,X2-1),s0(X1,y1 ,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh3d( t0(X1,y1+1,X2),
+ t0(X1-1,y1+1,X2),t0(X1-1,y1+1,X2-1),t0(X1,y1+1,X2-1),
+ t0(X1-1,y1 ,X2),t0(X1-1,y1 ,X2-1),t0(X1,y1 ,X2-1),
+ s0(X1,y1,X2), s0(X1,y1+1,X2),
+ s0(X1-1,y1+1,X2),s0(X1-1,y1+1,X2-1),s0(X1,y1+1,X2-1),
+ s0(X1-1,y1 ,X2),s0(X1-1,y1 ,X2-1),s0(X1,y1 ,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(guess > 1.0e9){
+ if(time0[index+1] < 1.e9 && X1<nx-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(X1,y1+1,X2),t0(X1+1,y1+1,X2),t0(X1+1,y1,X2),
+ t0(X1+1,y1+1,X2-1),t0(X1+1,y1+1,X2+1),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1+1,y1+1,X2),s0(X1+1,y1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(X1,y1+1,X2),t0(X1-1,y1+1,X2),t0(X1-1,y1,X2),
+ t0(X1-1,y1+1,X2-1),t0(X1-1,y1+1,X2+1),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1-1,y1+1,X2),s0(X1-1,y1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,y1+1,X2),t0(X1,y1+1,X2+1),t0(X1,y1,X2+1),
+ t0(X1-1,y1+1,X2+1),t0(X1+1,y1+1,X2+1),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1,y1+1,X2+1),s0(X1,y1,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,y1+1,X2),t0(X1,y1+1,X2-1),t0(X1,y1,X2-1),
+ t0(X1-1,y1+1,X2-1),t0(X1+1,y1+1,X2-1),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1,y1+1,X2-1),s0(X1,y1,X2-1) );
+ if (try<guess) guess = try;
+ }
+ }
+ if(time0[index+1] < 1.e9 && X1<nx-1 ) {
+ try = fdh2d(t0(X1,y1+1,X2),t0(X1+1,y1+1,X2),t0(X1+1,y1,X2),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1+1,y1+1,X2),s0(X1+1,y1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 ) {
+ try = fdh2d(t0(X1,y1+1,X2),t0(X1-1,y1+1,X2),t0(X1-1,y1,X2),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1-1,y1+1,X2),s0(X1-1,y1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 ) {
+ try = fdh2d(t0(X1,y1+1,X2),t0(X1,y1+1,X2+1),t0(X1,y1,X2+1),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1,y1+1,X2+1),s0(X1,y1,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 ) {
+ try = fdh2d(t0(X1,y1+1,X2),t0(X1,y1+1,X2-1),t0(X1,y1,X2-1),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2),s0(X1,y1+1,X2-1),s0(X1,y1,X2-1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,y1,X2),t0(X1+1,y1,X2+1),t0(X1,y1,X2+1),
+ s0(X1,y1,X2),
+ s0(X1+1,y1,X2),s0(X1+1,y1,X2+1),s0(X1,y1,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,y1,X2),t0(X1+1,y1,X2-1),t0(X1,y1,X2-1),
+ s0(X1,y1,X2),
+ s0(X1+1,y1,X2),s0(X1+1,y1,X2-1),s0(X1,y1,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh2d(t0(X1-1,y1,X2),t0(X1-1,y1,X2+1),t0(X1,y1,X2+1),
+ s0(X1,y1,X2),
+ s0(X1-1,y1,X2),s0(X1-1,y1,X2+1),s0(X1,y1,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh2d(t0(X1-1,y1,X2),t0(X1-1,y1,X2-1),t0(X1,y1,X2-1),
+ s0(X1,y1,X2),
+ s0(X1-1,y1,X2),s0(X1-1,y1,X2-1),s0(X1,y1,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(guess > 1.0e9){
+ if ( X1>x1+1 && X1<x2-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhnf(t0(X1,y1+1,X2),
+ t0(X1+1,y1+1,X2),t0(X1,y1+1,X2+1),
+ t0(X1-1,y1+1,X2),t0(X1,y1+1,X2-1),
+ s0(X1,y1,X2),
+ s0(X1,y1+1,X2) );
+ if (try<guess) guess = try;
+ }
+ }
+ try = t0(X1,y1+1,X2) + .5*(s0(X1,y1,X2)+s0(X1,y1+1,X2));
+ if (try<guess) guess = try;
+ if ( time0[index+1]<1.e9 && X1<nx-1 ) {
+ try = t0(X1+1,y1,X2) + .5*(s0(X1,y1,X2)+s0(X1+1,y1,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-1]<1.e9 && X1>0 ) {
+ try = t0(X1-1,y1,X2) + .5*(s0(X1,y1,X2)+s0(X1-1,y1,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index+nxy]<1.e9 && X2<nz-1 ) {
+ try = t0(X1,y1,X2+1) + .5*(s0(X1,y1,X2)+s0(X1,y1,X2+1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nxy]<1.e9 && X2>0 ) {
+ try = t0(X1,y1,X2-1) + .5*(s0(X1,y1,X2)+s0(X1,y1,X2-1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if (guess<time0[index]) {
+ time0[index] = guess;
+ if (fhead>headtest) headw[3]++;
+ }
+ }
+ if(y1 == 0) dy1 = 0;
+ y1--;
+ }
+ }
+ /* BACK SIDE */
+ for (igrow=1;igrow<=jplus;igrow++) {
+ if(dy2){
+ ii = 0;
+ for(k=z1+1; k<=z2-1; k++){
+ for(i=x1+1; i<=x2-1; i++){
+ sort[ii].time = t0(i,y2-1,k);
+ sort[ii].i1 = i;
+ sort[ii].i2 = k;
+ ii++;
+ }
+ }
+ qsort((char *)sort,ii,sizeof(struct sorted),compar);
+ for(i=0;i<ii;i++){
+ X1 = sort[i].i1;
+ X2 = sort[i].i2;
+ index = X2*nxy + y2*nx + X1;
+ lasti = X2*nxy + (y2-1)*nx + X1;
+ fhead = 0.;
+ guess = time0[index];
+ if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,y2-1,X2),
+ t0(X1+1,y2-1,X2),t0(X1+1,y2-1,X2+1),t0(X1,y2-1,X2+1),
+ t0(X1+1,y2 ,X2),t0(X1+1,y2 ,X2+1),t0(X1,y2 ,X2+1),
+ s0(X1,y2,X2), s0(X1,y2-1,X2),
+ s0(X1+1,y2-1,X2),s0(X1+1,y2-1,X2+1),s0(X1,y2-1,X2+1),
+ s0(X1+1,y2 ,X2),s0(X1+1,y2 ,X2+1),s0(X1,y2 ,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh3d( t0(X1,y2-1,X2),
+ t0(X1-1,y2-1,X2),t0(X1-1,y2-1,X2+1),t0(X1,y2-1,X2+1),
+ t0(X1-1,y2 ,X2),t0(X1-1,y2 ,X2+1),t0(X1,y2 ,X2+1),
+ s0(X1,y2,X2), s0(X1,y2-1,X2),
+ s0(X1-1,y2-1,X2),s0(X1-1,y2-1,X2+1),s0(X1,y2-1,X2+1),
+ s0(X1-1,y2 ,X2),s0(X1-1,y2 ,X2+1),s0(X1,y2 ,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh3d( t0(X1,y2-1,X2),
+ t0(X1+1,y2-1,X2),t0(X1+1,y2-1,X2-1),t0(X1,y2-1,X2-1),
+ t0(X1+1,y2 ,X2),t0(X1+1,y2 ,X2-1),t0(X1,y2 ,X2-1),
+ s0(X1,y2,X2), s0(X1,y2-1,X2),
+ s0(X1+1,y2-1,X2),s0(X1+1,y2-1,X2-1),s0(X1,y2-1,X2-1),
+ s0(X1+1,y2 ,X2),s0(X1+1,y2 ,X2-1),s0(X1,y2 ,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh3d( t0(X1,y2-1,X2),
+ t0(X1-1,y2-1,X2),t0(X1-1,y2-1,X2-1),t0(X1,y2-1,X2-1),
+ t0(X1-1,y2 ,X2),t0(X1-1,y2 ,X2-1),t0(X1,y2 ,X2-1),
+ s0(X1,y2,X2), s0(X1,y2-1,X2),
+ s0(X1-1,y2-1,X2),s0(X1-1,y2-1,X2-1),s0(X1,y2-1,X2-1),
+ s0(X1-1,y2 ,X2),s0(X1-1,y2 ,X2-1),s0(X1,y2 ,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(guess > 1.0e9){
+ if(time0[index+1] < 1.e9 && X1<nx-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(X1,y2-1,X2),t0(X1+1,y2-1,X2),t0(X1+1,y2,X2),
+ t0(X1+1,y2-1,X2-1),t0(X1+1,y2-1,X2+1),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1+1,y2-1,X2),s0(X1+1,y2,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(X1,y2-1,X2),t0(X1-1,y2-1,X2),t0(X1-1,y2,X2),
+ t0(X1-1,y2-1,X2-1),t0(X1-1,y2-1,X2+1),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1-1,y2-1,X2),s0(X1-1,y2,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,y2-1,X2),t0(X1,y2-1,X2+1),t0(X1,y2,X2+1),
+ t0(X1-1,y2-1,X2+1),t0(X1+1,y2-1,X2+1),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1,y2-1,X2+1),s0(X1,y2,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 ) {
+ try = fdhne(t0(X1,y2-1,X2),t0(X1,y2-1,X2-1),t0(X1,y2,X2-1),
+ t0(X1-1,y2-1,X2-1),t0(X1+1,y2-1,X2-1),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1,y2-1,X2-1),s0(X1,y2,X2-1) );
+ if (try<guess) guess = try;
+ }
+ }
+ if(time0[index+1] < 1.e9 && X1<nx-1 ) {
+ try = fdh2d(t0(X1,y2-1,X2),t0(X1+1,y2-1,X2),t0(X1+1,y2,X2),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1+1,y2-1,X2),s0(X1+1,y2,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-1] < 1.e9 && X1>0 ) {
+ try = fdh2d(t0(X1,y2-1,X2),t0(X1-1,y2-1,X2),t0(X1-1,y2,X2),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1-1,y2-1,X2),s0(X1-1,y2,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 ) {
+ try = fdh2d(t0(X1,y2-1,X2),t0(X1,y2-1,X2+1),t0(X1,y2,X2+1),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1,y2-1,X2+1),s0(X1,y2,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 ) {
+ try = fdh2d(t0(X1,y2-1,X2),t0(X1,y2-1,X2-1),t0(X1,y2,X2-1),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2),s0(X1,y2-1,X2-1),s0(X1,y2,X2-1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,y2,X2),t0(X1+1,y2,X2+1),t0(X1,y2,X2+1),
+ s0(X1,y2,X2),
+ s0(X1+1,y2,X2),s0(X1+1,y2,X2+1),s0(X1,y2,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<nx-1 ) {
+ try = fdh2d(t0(X1+1,y2,X2),t0(X1+1,y2,X2-1),t0(X1,y2,X2-1),
+ s0(X1,y2,X2),
+ s0(X1+1,y2,X2),s0(X1+1,y2,X2-1),s0(X1,y2,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh2d(t0(X1-1,y2,X2),t0(X1-1,y2,X2+1),t0(X1,y2,X2+1),
+ s0(X1,y2,X2),
+ s0(X1-1,y2,X2),s0(X1-1,y2,X2+1),s0(X1,y2,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh2d(t0(X1-1,y2,X2),t0(X1-1,y2,X2-1),t0(X1,y2,X2-1),
+ s0(X1,y2,X2),
+ s0(X1-1,y2,X2),s0(X1-1,y2,X2-1),s0(X1,y2,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(guess > 1.0e9){
+ if ( X1>x1+1 && X1<x2-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhnf(t0(X1,y2-1,X2),
+ t0(X1+1,y2-1,X2),t0(X1,y2-1,X2+1),
+ t0(X1-1,y2-1,X2),t0(X1,y2-1,X2-1),
+ s0(X1,y2,X2),
+ s0(X1,y2-1,X2) );
+ if (try<guess) guess = try;
+ }
+ }
+ try = t0(X1,y2-1,X2) + .5*(s0(X1,y2,X2)+s0(X1,y2-1,X2));
+ if (try<guess) guess = try;
+ if ( time0[index+1]<1.e9 && X1<nx-1 ) {
+ try = t0(X1+1,y2,X2) + .5*(s0(X1,y2,X2)+s0(X1+1,y2,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-1]<1.e9 && X1>0 ) {
+ try = t0(X1-1,y2,X2) + .5*(s0(X1,y2,X2)+s0(X1-1,y2,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index+nxy]<1.e9 && X2<nz-1 ) {
+ try = t0(X1,y2,X2+1) + .5*(s0(X1,y2,X2)+s0(X1,y2,X2+1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nxy]<1.e9 && X2>0 ) {
+ try = t0(X1,y2,X2-1) + .5*(s0(X1,y2,X2)+s0(X1,y2,X2-1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if (guess<time0[index]) {
+ time0[index] = guess;
+ if (fhead>headtest) headw[4]++;
+ }
+ }
+ if(y2 == ny-1) dy2 = 0;
+ y2++;
+ }
+ }
+ /* LEFT SIDE */
+ for (igrow=1;igrow<=iminus;igrow++) {
+ if(dx1){
+ ii = 0;
+ for(k=z1+1; k<=z2-1; k++){
+ for(j=y1+1; j<=y2-1; j++){
+ sort[ii].time = t0(x1+1,j,k);
+ sort[ii].i1 = j;
+ sort[ii].i2 = k;
+ ii++;
+ }
+ }
+ qsort((char *)sort,ii,sizeof(struct sorted),compar);
+ for(i=0;i<ii;i++){
+ X1 = sort[i].i1;
+ X2 = sort[i].i2;
+ index = X2*nxy + X1*nx + x1;
+ lasti = X2*nxy + X1*nx + (x1+1);
+ fhead = 0.;
+ guess = time0[index];
+ if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<ny-1 ) {
+ try = fdh3d( t0(x1+1,X1,X2),
+ t0(x1+1,X1+1,X2),t0(x1+1,X1+1,X2+1),t0(x1+1,X1,X2+1),
+ t0(x1 ,X1+1,X2),t0(x1 ,X1+1,X2+1),t0(x1 ,X1,X2+1),
+ s0(x1,X1,X2), s0(x1+1,X1,X2),
+ s0(x1+1,X1+1,X2),s0(x1+1,X1+1,X2+1),s0(x1+1,X1,X2+1),
+ s0(x1 ,X1+1,X2),s0(x1 ,X1+1,X2+1),s0(x1 ,X1,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh3d( t0(x1+1,X1,X2),
+ t0(x1+1,X1-1,X2),t0(x1+1,X1-1,X2+1),t0(x1+1,X1,X2+1),
+ t0(x1 ,X1-1,X2),t0(x1 ,X1-1,X2+1),t0(x1 ,X1,X2+1),
+ s0(x1,X1,X2), s0(x1+1,X1,X2),
+ s0(x1+1,X1-1,X2),s0(x1+1,X1-1,X2+1),s0(x1+1,X1,X2+1),
+ s0(x1 ,X1-1,X2),s0(x1 ,X1-1,X2+1),s0(x1 ,X1,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<ny-1 ) {
+ try = fdh3d( t0(x1+1,X1,X2),
+ t0(x1+1,X1+1,X2),t0(x1+1,X1+1,X2-1),t0(x1+1,X1,X2-1),
+ t0(x1 ,X1+1,X2),t0(x1 ,X1+1,X2-1),t0(x1 ,X1,X2-1),
+ s0(x1,X1,X2), s0(x1+1,X1,X2),
+ s0(x1+1,X1+1,X2),s0(x1+1,X1+1,X2-1),s0(x1+1,X1,X2-1),
+ s0(x1 ,X1+1,X2),s0(x1 ,X1+1,X2-1),s0(x1 ,X1,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh3d( t0(x1+1,X1,X2),
+ t0(x1+1,X1-1,X2),t0(x1+1,X1-1,X2-1),t0(x1+1,X1,X2-1),
+ t0(x1 ,X1-1,X2),t0(x1 ,X1-1,X2-1),t0(x1 ,X1,X2-1),
+ s0(x1,X1,X2), s0(x1+1,X1,X2),
+ s0(x1+1,X1-1,X2),s0(x1+1,X1-1,X2-1),s0(x1+1,X1,X2-1),
+ s0(x1 ,X1-1,X2),s0(x1 ,X1-1,X2-1),s0(x1 ,X1,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(guess > 1.0e9){
+ if(time0[index+nx] < 1.e9 && X1<ny-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1+1,X2),t0(x1,X1+1,X2),
+ t0(x1+1,X1+1,X2-1),t0(x1+1,X1+1,X2+1),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1+1,X2),s0(x1,X1+1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1-1,X2),t0(x1,X1-1,X2),
+ t0(x1+1,X1-1,X2-1),t0(x1+1,X1-1,X2+1),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1-1,X2),s0(x1,X1-1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>y1+1 && X1<y2-1 ) {
+ try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1,X2+1),t0(x1,X1,X2+1),
+ t0(x1+1,X1-1,X2+1),t0(x1+1,X1+1,X2+1),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1,X2+1),s0(x1,X1,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 && X1>y1+1 && X1<y2-1 ) {
+ try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1,X2-1),t0(x1,X1,X2-1),
+ t0(x1+1,X1-1,X2-1),t0(x1+1,X1+1,X2-1),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1,X2-1),s0(x1,X1,X2-1) );
+ if (try<guess) guess = try;
+ }
+ }
+ if(time0[index+nx] < 1.e9 && X1<ny-1 ) {
+ try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1+1,X2),t0(x1,X1+1,X2),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1+1,X2),s0(x1,X1+1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X1>0 ) {
+ try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1-1,X2),t0(x1,X1-1,X2),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1-1,X2),s0(x1,X1-1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 ) {
+ try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1,X2+1),t0(x1,X1,X2+1),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1,X2+1),s0(x1,X1,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 ) {
+ try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1,X2-1),t0(x1,X1,X2-1),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2),s0(x1+1,X1,X2-1),s0(x1,X1,X2-1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<ny-1 ) {
+ try = fdh2d(t0(x1,X1+1,X2),t0(x1,X1+1,X2+1),t0(x1,X1,X2+1),
+ s0(x1,X1,X2),
+ s0(x1,X1+1,X2),s0(x1,X1+1,X2+1),s0(x1,X1,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<ny-1 ) {
+ try = fdh2d(t0(x1,X1+1,X2),t0(x1,X1+1,X2-1),t0(x1,X1,X2-1),
+ s0(x1,X1,X2),
+ s0(x1,X1+1,X2),s0(x1,X1+1,X2-1),s0(x1,X1,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh2d(t0(x1,X1-1,X2),t0(x1,X1-1,X2+1),t0(x1,X1,X2+1),
+ s0(x1,X1,X2),
+ s0(x1,X1-1,X2),s0(x1,X1-1,X2+1),s0(x1,X1,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh2d(t0(x1,X1-1,X2),t0(x1,X1-1,X2-1),t0(x1,X1,X2-1),
+ s0(x1,X1,X2),
+ s0(x1,X1-1,X2),s0(x1,X1-1,X2-1),s0(x1,X1,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(guess > 1.0e9){
+ if ( X1>y1+1 && X1<y2-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhnf(t0(x1+1,X1,X2),
+ t0(x1+1,X1+1,X2),t0(x1+1,X1,X2+1),
+ t0(x1+1,X1-1,X2),t0(x1+1,X1,X2-1),
+ s0(x1,X1,X2),
+ s0(x1+1,X1,X2) );
+ if (try<guess) guess = try;
+ }
+ }
+ try = t0(x1+1,X1,X2) + .5*(s0(x1,X1,X2)+s0(x1+1,X1,X2));
+ if (try<guess) guess = try;
+ if ( time0[index+nx]<1.e9 && X1<ny-1 ) {
+ try = t0(x1,X1+1,X2) + .5*(s0(x1,X1,X2)+s0(x1,X1+1,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nx]<1.e9 && X1>0 ) {
+ try = t0(x1,X1-1,X2) + .5*(s0(x1,X1,X2)+s0(x1,X1-1,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index+nxy]<1.e9 && X2<nz-1 ) {
+ try = t0(x1,X1,X2+1) + .5*(s0(x1,X1,X2)+s0(x1,X1,X2+1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nxy]<1.e9 && X2>0 ) {
+ try = t0(x1,X1,X2-1) + .5*(s0(x1,X1,X2)+s0(x1,X1,X2-1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if (guess<time0[index]) {
+ time0[index] = guess;
+ if (fhead>headtest) headw[1]++;
+ }
+ }
+ if(x1 == 0) dx1 = 0;
+ x1--;
+ }
+ }
+ /* RIGHT SIDE */
+ for (igrow=1;igrow<=iplus;igrow++) {
+ if(dx2){
+ ii = 0;
+ for(k=z1+1; k<=z2-1; k++){
+ for(j=y1+1; j<=y2-1; j++){
+ sort[ii].time = t0(x2-1,j,k);
+ sort[ii].i1 = j;
+ sort[ii].i2 = k;
+ ii++;
+ }
+ }
+ qsort((char *)sort,ii,sizeof(struct sorted),compar);
+ for(i=0;i<ii;i++){
+ X1 = sort[i].i1;
+ X2 = sort[i].i2;
+ index = X2*nxy + X1*nx + x2;
+ lasti = X2*nxy + X1*nx + (x2-1);
+ fhead = 0.;
+ guess = time0[index];
+ if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<ny-1 ) {
+ try = fdh3d( t0(x2-1,X1,X2),
+ t0(x2-1,X1+1,X2),t0(x2-1,X1+1,X2+1),t0(x2-1,X1,X2+1),
+ t0(x2 ,X1+1,X2),t0(x2 ,X1+1,X2+1),t0(x2 ,X1,X2+1),
+ s0(x2,X1,X2), s0(x2-1,X1,X2),
+ s0(x2-1,X1+1,X2),s0(x2-1,X1+1,X2+1),s0(x2-1,X1,X2+1),
+ s0(x2 ,X1+1,X2),s0(x2 ,X1+1,X2+1),s0(x2 ,X1,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh3d( t0(x2-1,X1,X2),
+ t0(x2-1,X1-1,X2),t0(x2-1,X1-1,X2+1),t0(x2-1,X1,X2+1),
+ t0(x2 ,X1-1,X2),t0(x2 ,X1-1,X2+1),t0(x2 ,X1,X2+1),
+ s0(x2,X1,X2), s0(x2-1,X1,X2),
+ s0(x2-1,X1-1,X2),s0(x2-1,X1-1,X2+1),s0(x2-1,X1,X2+1),
+ s0(x2 ,X1-1,X2),s0(x2 ,X1-1,X2+1),s0(x2 ,X1,X2+1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<ny-1 ) {
+ try = fdh3d( t0(x2-1,X1,X2),
+ t0(x2-1,X1+1,X2),t0(x2-1,X1+1,X2-1),t0(x2-1,X1,X2-1),
+ t0(x2 ,X1+1,X2),t0(x2 ,X1+1,X2-1),t0(x2 ,X1,X2-1),
+ s0(x2,X1,X2), s0(x2-1,X1,X2),
+ s0(x2-1,X1+1,X2),s0(x2-1,X1+1,X2-1),s0(x2-1,X1,X2-1),
+ s0(x2 ,X1+1,X2),s0(x2 ,X1+1,X2-1),s0(x2 ,X1,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh3d( t0(x2-1,X1,X2),
+ t0(x2-1,X1-1,X2),t0(x2-1,X1-1,X2-1),t0(x2-1,X1,X2-1),
+ t0(x2 ,X1-1,X2),t0(x2 ,X1-1,X2-1),t0(x2 ,X1,X2-1),
+ s0(x2,X1,X2), s0(x2-1,X1,X2),
+ s0(x2-1,X1-1,X2),s0(x2-1,X1-1,X2-1),s0(x2-1,X1,X2-1),
+ s0(x2 ,X1-1,X2),s0(x2 ,X1-1,X2-1),s0(x2 ,X1,X2-1));
+ if (try<guess) guess = try;
+ }
+ if(guess > 1.0e9){
+ if(time0[index+nx] < 1.e9 && X1<ny-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1+1,X2),t0(x2,X1+1,X2),
+ t0(x2-1,X1+1,X2-1),t0(x2-1,X1+1,X2+1),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1+1,X2),s0(x2,X1+1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1-1,X2),t0(x2,X1-1,X2),
+ t0(x2-1,X1-1,X2-1),t0(x2-1,X1-1,X2+1),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1-1,X2),s0(x2,X1-1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>y1+1 && X1<y2-1 ) {
+ try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1,X2+1),t0(x2,X1,X2+1),
+ t0(x2-1,X1-1,X2+1),t0(x2-1,X1+1,X2+1),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1,X2+1),s0(x2,X1,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 && X1>y1+1 && X1<y2-1 ) {
+ try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1,X2-1),t0(x2,X1,X2-1),
+ t0(x2-1,X1-1,X2-1),t0(x2-1,X1+1,X2-1),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1,X2-1),s0(x2,X1,X2-1) );
+ if (try<guess) guess = try;
+ }
+ }
+ if(time0[index+nx] < 1.e9 && X1<ny-1 ) {
+ try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1+1,X2),t0(x2,X1+1,X2),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1+1,X2),s0(x2,X1+1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nx] < 1.e9 && X1>0 ) {
+ try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1-1,X2),t0(x2,X1-1,X2),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1-1,X2),s0(x2,X1-1,X2) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nxy] < 1.e9 && X2<nz-1 ) {
+ try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1,X2+1),t0(x2,X1,X2+1),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1,X2+1),s0(x2,X1,X2+1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index-nxy] < 1.e9 && X2>0 ) {
+ try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1,X2-1),t0(x2,X1,X2-1),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2),s0(x2-1,X1,X2-1),s0(x2,X1,X2-1) );
+ if (try<guess) guess = try;
+ }
+ if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1<ny-1 ) {
+ try = fdh2d(t0(x2,X1+1,X2),t0(x2,X1+1,X2+1),t0(x2,X1,X2+1),
+ s0(x2,X1,X2),
+ s0(x2,X1+1,X2),s0(x2,X1+1,X2+1),s0(x2,X1,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1<ny-1 ) {
+ try = fdh2d(t0(x2,X1+1,X2),t0(x2,X1+1,X2-1),t0(x2,X1,X2-1),
+ s0(x2,X1,X2),
+ s0(x2,X1+1,X2),s0(x2,X1+1,X2-1),s0(x2,X1,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
+ && time0[index+nxy] < 1.e9 && X2<nz-1 && X1>0 ) {
+ try = fdh2d(t0(x2,X1-1,X2),t0(x2,X1-1,X2+1),t0(x2,X1,X2+1),
+ s0(x2,X1,X2),
+ s0(x2,X1-1,X2),s0(x2,X1-1,X2+1),s0(x2,X1,X2+1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
+ && time0[index-nxy] < 1.e9 && X2>0 && X1>0 ) {
+ try = fdh2d(t0(x2,X1-1,X2),t0(x2,X1-1,X2-1),t0(x2,X1,X2-1),
+ s0(x2,X1,X2),
+ s0(x2,X1-1,X2),s0(x2,X1-1,X2-1),s0(x2,X1,X2-1) );
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if(guess > 1.0e9){
+ if ( X1>y1+1 && X1<y2-1 && X2>z1+1 && X2<z2-1 ) {
+ try = fdhnf(t0(x2-1,X1,X2),
+ t0(x2-1,X1+1,X2),t0(x2-1,X1,X2+1),
+ t0(x2-1,X1-1,X2),t0(x2-1,X1,X2-1),
+ s0(x2,X1,X2),
+ s0(x2-1,X1,X2) );
+ if (try<guess) guess = try;
+ }
+ }
+ try = t0(x2-1,X1,X2) + .5*(s0(x2,X1,X2)+s0(x2-1,X1,X2));
+ if (try<guess) guess = try;
+ if ( time0[index+nx]<1.e9 && X1<ny-1 ) {
+ try = t0(x2,X1+1,X2) + .5*(s0(x2,X1,X2)+s0(x2,X1+1,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nx]<1.e9 && X1>0 ) {
+ try = t0(x2,X1-1,X2) + .5*(s0(x2,X1,X2)+s0(x2,X1-1,X2));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index+nxy]<1.e9 && X2<nz-1 ) {
+ try = t0(x2,X1,X2+1) + .5*(s0(x2,X1,X2)+s0(x2,X1,X2+1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if ( time0[index-nxy]<1.e9 && X2>0 ) {
+ try = t0(x2,X1,X2-1) + .5*(s0(x2,X1,X2)+s0(x2,X1,X2-1));
+ if (try<guess) {fhead=(guess-try)/slow0[index]; guess=try;}
+ }
+ if (guess<time0[index]) {
+ time0[index] = guess;
+ if (fhead>headtest) headw[2]++;
+ }
+ }
+ if(x2 == nx-1) dx2 = 0;
+ x2++;
+ }
+ }
+
+ /* UPDATE RADIUS */
+ radius++;
+ if(radius%10 == 0 && verbose) vmess("Completed radius = %d",radius);
+ if(radius == maxrad) rad0 = 0;
+
+ } /* END BIG LOOP */
+
+
+ /* TEST IF REVERSE PROPAGATION IS NEEDED */
+
+ if (headw[1]==0 && headw[2]==0 && headw[3]==0 && headw[4]==0
+ && headw[5]==0 && headw[6]==0)
+ reverse=0;
+ else {
+ head=0;
+ if (headw[1]>0) {
+ if(verbose) vmess("Head waves found on left: %d",headw[1]);
+ if (headw[1]>head) {
+ head = headw[1];
+ srcwall = 1;
+ }
+ }
+ if (headw[2]>0) {
+ if(verbose) vmess("Head waves found on right: %d",headw[2]);
+ if (headw[2]>head) {
+ head = headw[2];
+ srcwall = 2;
+ }
+ }
+ if (headw[3]>0) {
+ if(verbose) vmess("Head waves found on front: %d",headw[3]);
+ if (headw[3]>head) {
+ head = headw[3];
+ srcwall = 3;
+ }
+ }
+ if (headw[4]>0) {
+ if(verbose) vmess("Head waves found on back: %d",headw[4]);
+ if (headw[4]>head) {
+ head = headw[4];
+ srcwall = 4;
+ }
+ }
+ if (headw[5]>0) {
+ if(verbose) vmess("Head waves found on top: %d",headw[5]);
+ if (headw[5]>head) {
+ head = headw[5];
+ srcwall = 5;
+ }
+ }
+ if (headw[6]>0) {
+ if(verbose) vmess("Head waves found on bottom: %d",headw[6]);
+ if (headw[6]>head) {
+ head = headw[6];
+ srcwall = 6;
+ }
+ }
+ if (headpref>0 && headw[headpref]>0) {
+ if(verbose)
+ vmess("Preference to restart on wall opposite source");
+ srcwall = headpref;
+ }
+ /* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE */
+ dx1=1; dx2=1; dy1=1; dy2=1; dz1=1; dz2=1; rad0=1;
+ radius = 1;
+ if (srcwall == 1) { x2=1;
+ vmess("RESTART at left side of model"); }
+ else { x2=nx; dx2=0; }
+ if (srcwall == 2) { x1=nx-2;
+ vmess("RESTART at right side of model"); }
+ else { x1= -1; dx1=0; }
+ if (srcwall == 3) { y2=1;
+ vmess("RESTART at front side of model"); }
+ else { y2=ny; dy2=0; }
+ if (srcwall == 4) { y1=ny-2;
+ vmess("RESTART at back side of model"); }
+ else { y1= -1; dy1=0; }
+ if (srcwall == 5) { z2=1;
+ vmess("RESTART at top side of model"); }
+ else { z2=nz; dz2=0; }
+ if (srcwall == 6) { z1=nz-2;
+ vmess("RESTART at bottom side of model"); }
+ else { z1= -1; dz1=0; }
+ if (reverse == 0)
+ vwarn("RESTART CANCELLED by choice of input parameter `reverse`");
+ }
+ reverse--;
+
+ } /* END BIGGER LOOP - HOLE */
+
+ free(sort);
+ free(wall);
+}
+
+int compar(struct sorted *a,struct sorted *b)
+{
+ if(a->time > b->time) return(1);
+ if(b->time > a->time) return(-1);
+ else return(0);
+}
+
+
+/* 3D TRANSMISSION STENCIL
+ STENCIL FROM VIDALE; CONDITIONS AND OTHER OPTIONS FROM HOLE
+ JAH 11/91 */
+float fdh3d(t1,t2,t3,t4,t5,t6,t7,ss0,s1,s2,s3,s4,s5,s6,s7)
+ float t1,t2,t3,t4,t5,t6,t7,ss0,s1,s2,s3,s4,s5,s6,s7;
+ /* ss0 at newpoint; s1,t1 adjacent on oldface;
+ s2,t2 and s4,t4 on oldface adjacent to s1;
+ s3,t3 on oldface diametrically opposite newpoint;
+ s5,t5 on newface adjacent to newpoint AND to s2;
+ s6,t6 on newface diagonal to newpoint (adjacent to s3);
+ s7,t7 on newface adjacent to newpoint AND to s4
+ */
+{
+ float x,slo;
+ double sqrt();
+ slo = .125*(ss0+s1+s2+s3+s4+s5+s6+s7);
+ x = 6.*slo*slo - (t4-t2)*(t4-t2) - (t2-t6)*(t2-t6) - (t6-t4)*(t6-t4)
+ - (t7-t5)*(t7-t5) - (t5-t1)*(t5-t1) - (t1-t7)*(t1-t7);
+ if (x>=0.) {
+ x = t3 + sqrt(.5*x);
+ if ( (x<t1) || (x<t2) || (x<t4) || (x<t5) || (x<t6) || (x<t7) )
+ x = 1.e11; /* ACAUSAL; ABORT */
+ }
+ else x = 1.e11; /* SQRT IMAGINARY; ABORT */
+ return(x);
+}
+
+/* 3D STENCIL FOR NEW EDGE
+ STENCIL FROM VIDALE; CONDITIONS AND OTHER OPTIONS FROM HOLE
+ JAH 11/91 */
+float fdhne(t1,t2,t3,t4,t5,ss0,s1,s2,s3)
+ float t1,t2,t3,t4,t5,ss0,s1,s2,s3;
+ /* ss0 at newpoint; s1,t1 adjacent on oldface;
+ s2,t2 diagonal on oldface; s3,t3 adjacent on newface;
+ t4,t5 beside t2 on old face opposite each other */
+{
+ float x,slo;
+ double sqrt();
+ slo = .25*(ss0+s1+s2+s3);
+ x = 2.*slo*slo - (t3-t1)*(t3-t1) - .5*(t5-t4)*(t5-t4);
+ if (x>=0.) {
+ x = t2 + sqrt(x);
+ if ( (x<t1) || (x<t3) || (x<t4) || (x<t5) ) /* ACAUSAL; ABORT */
+ x = 1.e11;
+ }
+ else x = 1.e11; /* SQRT IMAGINARY; ABORT */
+ return(x);
+}
+
+/* 2D TRANSMISSION STENCIL (FOR HEAD WAVES ON FACES OF GRID CELLS)
+ STENCIL FROM VIDALE (1988 2D PAPER); CONDITIONS AND OTHER OPTIONS FROM HOLE
+ JAH 11/91 */
+float fdh2d(t1,t2,t3,ss0,s1,s2,s3)
+ float t1,t2,t3,ss0,s1,s2,s3;
+ /* ss0 at newpoint; s1,t1 & s3,t3 adjacent; s2,t2 diagonal
+ */
+{
+ float x,slo;
+ double sqrt();
+ slo = .25*(ss0+s1+s2+s3);
+ x = 2.*slo*slo - (t3-t1)*(t3-t1);
+ if (x>=0.) {
+ x = t2 + sqrt(x);
+ if ( (x<t1) || (x<t3) ) x = 1.e11; /* ACAUSAL; ABORT */
+ }
+ else x = 1.e11; /* SQRT IMAGINARY; ABORT */
+ return(x);
+}
+
+/* 3D STENCIL FOR NEW FACE
+ STENCIL FROM VIDALE; CONDITIONS AND OTHER OPTIONS FROM HOLE
+ JAH 11/91 */
+float fdhnf(t1,t2,t3,t4,t5,ss0,s1)
+ float t1,t2,t3,t4,t5,ss0,s1;
+ /* ss0 at newpoint; s1,t1 adjacent on old face;
+ t2,t4 beside t1 on old face and opposite each other;
+ t3,t5 beside t1 on old face and opposite each other
+ */
+{
+ float x,slo;
+ double sqrt();
+ slo = .5*(ss0+s1);
+ x = slo*slo - .25*( (t4-t2)*(t4-t2) + (t5-t3)*(t5-t3) );
+ if (x>=0.) {
+ x = t1 + sqrt(x);
+ if ( (x<t2) || (x<t3) || (x<t4) || (x<t5) ) /* ACAUSAL; ABORT */
+ x = 1.e11;
+ }
+ else x = 1.e11; /* SQRT IMAGINARY; ABORT */
+ return(x);
+}
+
+
diff --git a/raytime3d/wallclock_time.c b/raytime3d/wallclock_time.c
new file mode 100644
index 0000000000000000000000000000000000000000..1e75530ccee3215724badefdd6144c2c59246dbc
--- /dev/null
+++ b/raytime3d/wallclock_time.c
@@ -0,0 +1,33 @@
+#include <time.h>
+#include <sys/time.h>
+#include <stdio.h>
+
+/**
+* function used to calculate wallclock times
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+double wallclock_time(void)
+{
+ struct timeval s_val;
+ static struct timeval b_val;
+ double time;
+ static int base=0;
+
+ gettimeofday(&s_val,0);
+
+ if (!base) {
+ b_val = s_val;
+ base = 1;
+ return 0.0;
+ }
+
+ time = (double)(s_val.tv_sec-b_val.tv_sec) +
+ (double)(1e-6*((double)s_val.tv_usec-(double)b_val.tv_usec));
+
+ return (double)time;
+}
+
diff --git a/raytime3d/writeSrcRecPos.c b/raytime3d/writeSrcRecPos.c
new file mode 100644
index 0000000000000000000000000000000000000000..d9d03da16dae329b7628e3420224bda8d9a6cce4
--- /dev/null
+++ b/raytime3d/writeSrcRecPos.c
@@ -0,0 +1,136 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"par.h"
+#include"raytime3d.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* 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
+**/
+
+int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2, float d1, float d2);
+
+int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot)
+{
+ FILE *fp;
+ float *dum, sub_x0, sub_z0, dx, dz;
+ int is, nx, nz, is0, ish, ix, iz, ndot, idx, idz;
+ char tmpname[1024];
+
+ ndot = 2;
+ nx = mod->nx;
+ nz = mod->nz;
+ dx = mod->dx;
+ dz = mod->dz;
+ sub_x0 = mod->x0;
+ sub_z0 = mod->z0;
+
+// ibndx = mod.ioPx;
+// ibndz = mod.ioPz;
+// if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+// if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+
+ /* write velocity field with positions of the sources */
+ dum = (float *)calloc(nx*nz, sizeof(float));
+ vmess("Positions: shot=%d src=%d rec=%d", shot->n, src->n, rec->n);
+ /* source positions for random shots */
+ if (src->random) {
+ sprintf(tmpname,"SrcPositions%d.txt",src->n);
+ fp = fopen(tmpname, "w+");
+ for (is=0; is<src->n; is++) {
+ for (idx=0; idx<=ndot; idx++) {
+ for (idz=0; idz<=ndot; idz++) {
+ dum[(MAX(0,src->x[is]-idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+ dum[(MAX(0,src->x[is]-idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+ dum[(MIN(nx-1,src->x[is]+idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+ dum[(MIN(nx-1,src->x[is]+idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+ }
+ }
+ fprintf(fp, "%f %f\n", src->z[is]*dz+sub_z0, 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%d.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;
+ iz = shot->z[ish] + 1;
+ dum[ix*nz+iz] = 1.0;
+ dum[(MAX(0,ix-1))*nz+iz] = 1.0;
+ dum[(MIN(nx-1,ix+1))*nz+iz] = 1.0;
+ dum[ix*nz+MAX(0,iz-1)] = 1.0;
+ dum[ix*nz+MIN(nz-1,iz+1)] = 1.0;
+ fprintf(fp, "(%f, %f)\n", ix*dx+sub_x0, iz*dz+sub_z0);
+ }
+ }
+ fclose(fp);
+ }
+ else if (src->multiwav) {
+ /* source positions for single shot sources with multiple wavelets */
+ sprintf(tmpname,"SrcPositions%d.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];
+ iz = src->z[is];
+ dum[ix*nz+iz] = 1.0;
+ dum[(MAX(0,ix-1))*nz+iz] = 1.0;
+ dum[(MIN(nx-1,ix+1))*nz+iz] = 1.0;
+ dum[ix*nz+MAX(0,iz-1)] = 1.0;
+ dum[ix*nz+MIN(nz-1,iz+1)] = 1.0;
+ fprintf(fp, "(%f, %f)\n", ix*dx+sub_x0, iz*dz+sub_z0);
+ }
+ }
+ fclose(fp);
+ }
+ else {
+ sprintf(tmpname,"SrcPositions%d.txt",shot->n);
+ fp = fopen(tmpname, "w+");
+ for (is=0; is<shot->n; is++) {
+ for (idx=0; idx<=ndot; idx++) {
+ for (idz=0; idz<=ndot; idz++) {
+ dum[(MAX(0,shot->x[is]-idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+ dum[(MAX(0,shot->x[is]-idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+ dum[(MIN(nx-1,shot->x[is]+idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+ dum[(MIN(nx-1,shot->x[is]+idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+ }
+ }
+ fprintf(fp, "%f %f\n", shot->z[is]*dz+sub_z0, shot->x[is]*dx+sub_x0);
+ }
+ fclose(fp);
+ }
+
+ /* receiver positions */
+ sprintf(tmpname,"RcvPositions%d.txt",rec->n);
+ fp = fopen(tmpname, "w+");
+ for (is=0; is<rec->n; is++) {
+ dum[rec->x[is]*nz+rec->z[is]] = -1.0;
+ dum[(MAX(0,rec->x[is]-1))*nz+rec->z[is]] = -1.0;
+ dum[(MIN(nx-1,rec->x[is]+1))*nz+rec->z[is]] = -1.0;
+ dum[rec->x[is]*nz+MAX(0,rec->z[is]-1)] = -1.0;
+ dum[rec->x[is]*nz+MIN(nz-1,rec->z[is]+1)] = -1.0;
+
+// vmess("receiver position %d at grid[ix=%d, iz=%d] = (x=%f z=%f)", ir, ix+ioPx, rec.z[ir]+ioPz, rec.xr[ir]+mod.x0, rec.zr[ir]+mod.z0);
+ fprintf(fp, "(%f, %f)\n", rec->x[is]*dx+sub_x0, rec->z[is]*dz+sub_z0);
+ }
+ fclose(fp);
+ writesufile("SrcRecPositions.su", dum, nz, nx, sub_z0, sub_x0, dz, dx);
+ free(dum);
+
+ return 0;
+}
diff --git a/raytime3d/writesufile.c b/raytime3d/writesufile.c
new file mode 100644
index 0000000000000000000000000000000000000000..7ebcbbd198858193b67a2b9c2fbb5670feb8b5d8
--- /dev/null
+++ b/raytime3d/writesufile.c
@@ -0,0 +1,169 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <string.h>
+#include "par.h"
+#include "raytime3d.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) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2, float d1, float d2)
+{
+ FILE *file_out;
+ size_t nwrite, itrace;
+ int ns;
+ segy *hdr;
+// char *ptr;
+
+/* Read in parameters */
+
+// ptr = strstr(filename, " ");
+// *ptr = '\0';
+
+
+ if (n1 > USHRT_MAX) {
+ vwarn("Output file %s: number of samples is truncated from %d 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
+*
+**/
+
+int writesufilesrcnwav(char *filename, float **src_nwav, wavPar wav, int n1, int n2, float f1, float f2, float d1, float d2)
+{
+ FILE *file_out;
+ size_t nwrite, itrace;
+ float *trace;
+ int ns;
+ segy *hdr;
+// char *ptr;
+
+/* Read in parameters */
+
+// ptr = strstr(filename, " ");
+// *ptr = '\0';
+
+ if (n1 > USHRT_MAX) {
+ vwarn("Output file %s: number of samples is truncated from %d 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.
+*
+**/
+
+int writeSUfile(char *filename, float *data, int n1, int 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;
+}
+
diff --git a/utils/makemod.c b/utils/makemod.c
index d6bfb9782ed83dfbe77ab2b04740e5c539907e85..574c1355e646a5e06b536f5e07ef3174914ca88b 100644
--- a/utils/makemod.c
+++ b/utils/makemod.c
@@ -646,6 +646,7 @@ int main(int argc, char **argv)
float Wsmooth[5][5], C, iC, xx, xz, *dataS, smooth;
int ixi, izi, nxout, nzout;
+ C=0;
sigma = -1.0*log(sigma)/(dx*(powf(0.25*2.0,2.0)));
for(ixi = -2; ixi < 3; ixi++) {
for(izi = -2; izi < 3; izi++) {