diff --git a/raytime3d/Makefile b/raytime3d/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..1bed95dcdef7cd3fb664ca9017a3bad0b38788b3
--- /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 \
+ getParameters.c \
+ getWaveletInfo.c \
+ writeSrcRecPos.c \
+ readModel.c \
+ getWaveletHeaders.c \
+ verbosepkg.c \
+ getModelInfo.c \
+ wallclock_time.c \
+ recvPar.c \
+ writesufile.c \
+ name_ext.c \
+ atopkge.c \
+ docpkge.c \
+ threadAffinity.c \
+ getpars.c
+
+OBJC = $(SRCC:%.c=%.o)
+
+$(PRG): $(OBJC) raytime.h
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o Raytime $(OBJC) $(LIBS)
+
+install: raytime
+ cp Raytime $B
+
+clean:
+ rm -f core $(OBJC) $(OBJM) Raytime
+
+realclean:
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/Raytime
+
+
+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/getModelInfo.c b/raytime3d/getModelInfo.c
new file mode 100644
index 0000000000000000000000000000000000000000..378a1b50ebac46e5b7b8a8bef4b5365ac15bef9d
--- /dev/null
+++ b/raytime3d/getModelInfo.c
@@ -0,0 +1,109 @@
+#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 "par.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))
+
+/**
+* reads gridded model file to compute minimum and maximum values and sampling intervals
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose)
+{
+ FILE *fp;
+ size_t nread, trace_sz;
+ off_t bytes;
+ int ret, i, one_shot, ntraces;
+ float *trace, cmin;
+ segy hdr;
+
+ fp = fopen( file_name, "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;
+ *d1 = hdr.d1;
+ *d2 = hdr.d2;
+ *f1 = hdr.f1;
+ *f2 = hdr.f2;
+
+ if ( NINT(100.0*((*d1)/(*d2)))!=100 ) {
+ verr("dx and dz are different in the model !");
+ }
+ if ( NINT(1000.0*(*d1))==0 ) {
+ if(!getparfloat("dx",d1)) {
+ verr("dx is equal to zero use parameter dx= to set value");
+ }
+ *d2 = *d1;
+ }
+ trace_sz = sizeof(float)*(*n1)+TRCBYTES;
+ ntraces = (int) (bytes/trace_sz);
+ *n2 = ntraces;
+
+ /* check to find out min and max values gather */
+
+ one_shot = 1;
+ trace = (float *)malloc(trace_sz);
+ fseeko( fp, TRCBYTES, SEEK_SET );
+ nread = fread( trace, sizeof(float), hdr.ns, fp );
+ assert (nread == hdr.ns);
+ fseeko( fp, TRCBYTES, SEEK_SET );
+
+ if (hdr.trid == TRID_DEPTH) *axis = 1; /* samples are z-axis */
+ else *axis = 0; /* sample direction respresents the x-axis */
+
+ i=0; cmin=trace[0];
+ while ( ( (cmin==0.0) && zeroch) && (i<hdr.ns) ) cmin=trace[i++];
+
+ *max = cmin;
+ *min = cmin;
+ /* keep on reading traces until there are no more traces (nread==0) */
+ while (one_shot) {
+ nread = fread( trace, sizeof(float), hdr.ns, fp );
+ assert (nread == hdr.ns);
+ for (i=0;i<(*n1);i++) {
+ *max = MAX(trace[i],*max);
+ cmin = MIN(trace[i],*min);
+ if (zeroch) {
+ if (cmin!=0.0) *min = MIN(*min, cmin);
+ }
+ else {
+ *min = cmin;
+ }
+ }
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread==0) break;
+ }
+ fclose(fp);
+ free(trace);
+
+ if (verbose>2) {
+ vmess("For file %s", file_name);
+ vmess("nz=%d nx=%d", *n1, *n2);
+ vmess("dz=%f dx=%f", *d1, *d2);
+ vmess("min=%f max=%f", *min, *max);
+ vmess("zstart=%f xstart=%f", *f1, *f2);
+ if (*axis) vmess("sample represent z-axis\n");
+ else vmess("sample represent x-axis\n");
+ }
+ return 0;
+}
+
diff --git a/raytime3d/getParameters.c b/raytime3d/getParameters.c
new file mode 100644
index 0000000000000000000000000000000000000000..89d1d34585a977f9ba6fd58eff936bf4396e4ea2
--- /dev/null
+++ b/raytime3d/getParameters.c
@@ -0,0 +1,305 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"par.h"
+#include"raytime.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))
+
+/**
+*
+* The routine getParameters reads in all parameters to set up a FD modeling.
+* Model and source parameters are used to calculate stability and dispersion relations
+* Source and receiver positions are calculated and checked if they fit into the model.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose);
+
+int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx, int nz);
+
+int getParameters(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose)
+{
+ int nx, nz, nsrc, ix, axis, is0;
+ int idzshot, idxshot;
+ int src_ix0, src_iz0, src_ix1, src_iz1;
+ float cp_min, cp_max;
+ float sub_x0,sub_z0;
+ float srcendx, srcendz, dx, dz;
+ float xsrc, zsrc, dxshot, dzshot;
+ float dxrcv,dzrcv,dxspread,dzspread;
+ float xmax, zmax;
+ float xsrc1, xsrc2, zsrc1, zsrc2;
+ float *xsrca, *zsrca;
+ float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
+ size_t nsamp;
+ int nxsrc, nzsrc;
+ int is;
+ char *src_positions;
+
+ if (!getparint("verbose",&verbose)) verbose=0;
+
+ if (!getparstring("file_cp",&mod->file_cp)) {
+ verr("parameter file_cp required!");
+ }
+ if (!getparstring("file_rcv",&rec->file_rcv)) rec->file_rcv="recv.su";
+ if (!getparint("src_at_rcv",&src->src_at_rcv)) src->src_at_rcv=1;
+
+ /* read model parameters, which are used to set up source and receivers and check stability */
+
+ getModelInfo(mod->file_cp, &nz, &nx, &dz, &dx, &sub_z0, &sub_x0, &cp_min, &cp_max, &axis, 1, verbose);
+ mod->cp_max = cp_max;
+ mod->cp_min = cp_min;
+ mod->dz = dz;
+ mod->dx = dx;
+ mod->nz = nz;
+ mod->nx = nx;
+
+ /* origin of model in real (non-grid) coordinates */
+ mod->x0 = sub_x0;
+ mod->z0 = sub_z0;
+ xmax = sub_x0+(nx-1)*dx;
+ zmax = sub_z0+(nz-1)*dz;
+
+ if (verbose) {
+ vmess("*******************************************");
+ vmess("*************** model info ****************");
+ vmess("*******************************************");
+ vmess("nz = %8d nx = %8d", nz, nx);
+ vmess("dz = %8.4f dx = %8.4f", dz, dx);
+ vmess("zmin = %8.4f zmax = %8.4f", sub_z0, zmax);
+ vmess("xmin = %8.4f xmax = %8.4f", sub_x0, xmax);
+ vmess("min(cp) = %9.3f max(cp) = %9.3f", cp_min, cp_max);
+ }
+
+ /* define the number and type of shots to model */
+ /* each shot can have multiple sources arranged in different ways */
+
+ if (!getparfloat("xsrc",&xsrc)) xsrc=sub_x0+((nx-1)*dx)/2.0;
+ if (!getparfloat("zsrc",&zsrc)) zsrc=sub_z0;
+ if (!getparint("nxshot",&shot->nx)) shot->nx=1;
+ if (!getparint("nzshot",&shot->nz)) shot->nz=1;
+ if (!getparfloat("dxshot",&dxshot)) dxshot=dx;
+ if (!getparfloat("dzshot",&dzshot)) dzshot=dz;
+
+ shot->n = (shot->nx)*(shot->nz);
+
+ if (shot->nx>1) {
+ idxshot=MAX(0,NINT(dxshot/dx));
+ }
+ else {
+ idxshot=0.0;
+ }
+ if (shot->nz>1) {
+ idzshot=MAX(0,NINT(dzshot/dz));
+ }
+ else {
+ idzshot=0.0;
+ }
+
+ /* calculate the shot positions */
+
+ src_ix0=MAX(0,NINT((xsrc-sub_x0)/dx));
+ src_ix0=MIN(src_ix0,nx);
+ src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
+ src_iz0=MIN(src_iz0,nz);
+ srcendx=(shot->nx-1)*dxshot+xsrc;
+ srcendz=(shot->nz-1)*dzshot+zsrc;
+ src_ix1=MAX(0,NINT((srcendx-sub_x0)/dx));
+ src_ix1=MIN(src_ix1,nx);
+ src_iz1=MAX(0,NINT((srcendz-sub_z0)/dz));
+ src_iz1=MIN(src_iz1,nz);
+
+ shot->x = (int *)calloc(shot->nx,sizeof(int));
+ shot->z = (int *)calloc(shot->nz,sizeof(int));
+ for (is=0; is<shot->nx; is++) {
+ shot->x[is] = src_ix0+is*idxshot;
+ if (shot->x[is] > nx-1) shot->nx = is-1;
+ }
+ for (is=0; is<shot->nz; is++) {
+ shot->z[is] = src_iz0+is*idzshot;
+ if (shot->z[is] > nz-1) shot->nz = is-1;
+ }
+
+ /* check if source array is defined */
+
+ nxsrc = countparval("xsrca");
+ nzsrc = countparval("zsrca");
+ if (nxsrc != nzsrc) {
+ verr("Number of sources in array xsrca (%d), zsrca(%d) are not equal",nxsrc, nzsrc);
+ }
+
+ /* check if sources on a circle are defined */
+
+ if (getparfloat("rsrc", &rsrc)) {
+ if (!getparfloat("dphisrc",&dphisrc)) dphisrc=2.0;
+ if (!getparfloat("oxsrc",&oxsrc)) oxsrc=0.0;
+ if (!getparfloat("ozsrc",&ozsrc)) ozsrc=0.0;
+ ncsrc = NINT(360.0/dphisrc);
+ src->n = nsrc;
+
+ src->x = (int *)malloc(ncsrc*sizeof(int));
+ src->z = (int *)malloc(ncsrc*sizeof(int));
+
+ for (ix=0; ix<ncsrc; ix++) {
+ src->x[ix] = NINT((oxsrc-sub_x0+rsrc*cos(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dx);
+ src->z[ix] = NINT((ozsrc-sub_z0+rsrc*sin(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dz);
+ if (verbose>4) fprintf(stderr,"Source on Circle: xsrc[%d]=%d zsrc=%d\n", ix, src->x[ix], src->z[ix]);
+ }
+
+ }
+
+ /* TO DO propagate src_positions parameter and structure through code */
+
+ if (!getparstring("src_positions",&src_positions)) src_positions="single";
+ src->random=0;
+ src->plane=0;
+ src->array=0;
+ src->single=0;
+ if (strstr(src_positions, "single")) src->single=1;
+ else if (strstr(src_positions, "array")) src->array=1;
+ else if (strstr(src_positions, "random")) src->random=1;
+ else if (strstr(src_positions, "plane")) src->plane=1;
+ else src->single=1;
+
+ /* to maintain functionality of older parameters usage */
+ if (!getparint("src_random",&src->random)) src->random=0;
+ if (!getparint("plane_wave",&src->plane)) src->plane=0;
+
+ if (src->random) {
+ src->plane=0;
+ src->array=0;
+ src->single=0;
+ }
+ if (src->plane) {
+ src->random=0;
+ src->array=0;
+ src->single=0;
+ }
+
+
+ /* number of sources per shot modeling */
+
+ if (!getparint("src_window",&src->window)) src->window=0;
+ if (!getparint("distribution",&src->distribution)) src->distribution=0;
+ if (!getparfloat("amplitude", &src->amplitude)) src->amplitude=0.0;
+ if (src->random && nxsrc==0) {
+ if (!getparint("nsrc",&nsrc)) nsrc=1;
+ if (!getparfloat("xsrc1", &xsrc1)) xsrc1=sub_x0;
+ if (!getparfloat("xsrc2", &xsrc2)) xsrc2=xmax;
+ if (!getparfloat("zsrc1", &zsrc1)) zsrc1=sub_z0;
+ if (!getparfloat("zsrc2", &zsrc2)) zsrc2=zmax;
+ dxshot = xsrc2-xsrc1;
+ dzshot = zsrc2-zsrc1;
+ src->x = (int *)malloc(nsrc*sizeof(int));
+ src->z = (int *)malloc(nsrc*sizeof(int));
+ nsamp = 0;
+
+ }
+ else if (nxsrc != 0) {
+ /* source array is defined */
+ nsrc=nxsrc;
+ src->x = (int *)malloc(nsrc*sizeof(int));
+ src->z = (int *)malloc(nsrc*sizeof(int));
+ xsrca = (float *)malloc(nsrc*sizeof(float));
+ zsrca = (float *)malloc(nsrc*sizeof(float));
+ getparfloat("xsrca", xsrca);
+ getparfloat("zsrca", zsrca);
+ for (is=0; is<nsrc; is++) {
+ src->x[is] = NINT((xsrca[is]-sub_x0)/dx);
+ src->z[is] = NINT((zsrca[is]-sub_z0)/dz);
+ if (verbose>3) fprintf(stderr,"Source Array: xsrc[%d]=%f zsrc=%f\n", is, xsrca[is], zsrca[is]);
+ }
+ src->random = 1;
+ free(xsrca);
+ free(zsrca);
+ }
+ else {
+ if (src->plane) { if (!getparint("nsrc",&nsrc)) nsrc=1;}
+ else nsrc=1;
+
+ if (nsrc > nx) {
+ vwarn("Number of sources used in plane wave is larger than ");
+ vwarn("number of gridpoints in X. Plane wave will be clipped to the edges of the model");
+ nsrc = mod->nx;
+ }
+
+ /* for a source defined on mutliple gridpoint calculate p delay factor */
+
+ src->x = (int *)malloc(nsrc*sizeof(int));
+ src->z = (int *)malloc(nsrc*sizeof(int));
+ is0 = -1*floor((nsrc-1)/2);
+ for (is=0; is<nsrc; is++) {
+ src->x[is] = is0 + is;
+ src->z[is] = 0;
+ }
+
+ }
+
+ src->n=nsrc;
+
+ if (verbose) {
+ if (src->n>1) {
+ vmess("*******************************************");
+ vmess("*********** source array info *************");
+ vmess("*******************************************");
+ vmess("Areal source array is defined with %d sources.",nsrc);
+ vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
+ }
+ if (src->random) {
+ vmess("Sources are placed at random locations in domain: ");
+ vmess(" x[%.2f : %.2f] z[%.2f : %.2f] ", xsrc1, xsrc2, zsrc1, zsrc2);
+ }
+ }
+
+ /* define receivers */
+
+ if (!getparint("sinkdepth",&rec->sinkdepth)) rec->sinkdepth=0;
+ if (!getparint("sinkdepth_src",&src->sinkdepth)) src->sinkdepth=0;
+ if (!getparint("sinkvel",&rec->sinkvel)) rec->sinkvel=0;
+ if (!getparint("max_nrec",&rec->max_nrec)) rec->max_nrec=15000;
+ if (!getparfloat("dxspread",&dxspread)) dxspread=0;
+ if (!getparfloat("dzspread",&dzspread)) dzspread=0;
+ if (!getparint("nt",&rec->nt)) rec->nt=1024;
+
+ /* calculates the receiver coordinates */
+
+ recvPar(rec, sub_x0, sub_z0, dx, dz, nx, nz);
+
+ if (verbose) {
+ if (rec->n) {
+ dxrcv = rec->xr[MIN(1,rec->n-1)]-rec->xr[0];
+ dzrcv = rec->zr[MIN(1,rec->n-1)]-rec->zr[0];
+ vmess("*******************************************");
+ vmess("************* receiver info ***************");
+ vmess("*******************************************");
+ vmess("ntrcv = %d nrcv = %d ", rec->nt, rec->n);
+ vmess("dzrcv = %f dxrcv = %f ", dzrcv, dxrcv);
+ vmess("Receiver array at coordinates: ");
+ vmess("zmin = %f zmax = %f ", rec->zr[0]+sub_z0, rec->zr[rec->n-1]+sub_z0);
+ vmess("xmin = %f xmax = %f ", rec->xr[0]+sub_x0, rec->xr[rec->n-1]+sub_x0);
+ vmess("which are gridpoints: ");
+ vmess("izmin = %d izmax = %d ", rec->z[0], rec->z[rec->n-1]);
+ vmess("ixmin = %d ixmax = %d ", rec->x[0], rec->x[rec->n-1]);
+ fprintf(stderr,"\n");
+ }
+ else {
+ vmess("*************** no receivers **************");
+ }
+ }
+
+ /* Ray tracing parameters */
+ if (!getparint("smoothwindow",&ray->smoothwindow)) ray->smoothwindow=0;
+ if (!getparint("useT2",&ray->useT2)) ray->useT2=0;
+ if (!getparint("geomspread",&ray->geomspread)) ray->geomspread=1;
+ if (!getparint("nraystep",&ray->nray)) ray->nray=5;
+
+ return 0;
+}
+
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/raytime.h b/raytime3d/raytime.h
new file mode 100644
index 0000000000000000000000000000000000000000..882c78b9a2eb0b945737dac10ae5c6f07aa55123
--- /dev/null
+++ b/raytime3d/raytime.h
@@ -0,0 +1,135 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<limits.h>
+#include<float.h>
+#include<math.h>
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+typedef struct _dcomplexStruct { /* complex number */
+ double r,i;
+} dcomplex;
+#endif/* complex */
+
+
+typedef struct _icoord { /* 3D coordinate integer */
+ int z;
+ int x;
+ int y;
+} icoord;
+
+typedef struct _fcoord { /* 3D coordinate float */
+ float z;
+ float x;
+ float y;
+} fcoord;
+
+struct s_ecount {
+ int corner,corner_min,side;
+};
+
+typedef struct _receiverPar { /* Receiver Parameters */
+ char *file_rcv;
+ int n;
+ int nt;
+ int max_nrec;
+ int *z;
+ int *x;
+ float *zr;
+ float *xr;
+ int scale;
+ int sinkdepth;
+ int sinkvel;
+ float cp;
+ float rho;
+} recPar;
+
+typedef struct _modelPar { /* Model Parameters */
+ int sh;
+ char *file_cp;
+ float dz;
+ float dx;
+ float dt;
+ float z0;
+ float x0;
+ /* medium max/min values */
+ float cp_min;
+ float cp_max;
+ int nz;
+ int nx;
+} modPar;
+
+typedef struct _waveletPar { /* Wavelet Parameters */
+ char *file_src; /* general source */
+ int nsrcf;
+ int nt;
+ int ns;
+ int nx;
+ float dt;
+ float ds;
+ float fmax;
+ int random;
+ int seed;
+ int nst;
+ size_t *nsamp;
+} wavPar;
+
+typedef struct _sourcePar { /* Source Array Parameters */
+ int n;
+ int type;
+ int orient;
+ int *z;
+ int *x;
+ int single;
+ int plane;
+ int circle;
+ int array;
+ int random;
+ int multiwav;
+ float angle;
+ float velo;
+ float amplitude;
+ int distribution;
+ int window;
+ int injectionrate;
+ int sinkdepth;
+ int src_at_rcv; /* Indicates that wavefield should be injected at receivers */
+} srcPar;
+
+typedef struct _shotPar { /* Shot Parameters */
+ int n;
+ int nx;
+ int nz;
+ int *z;
+ int *x;
+} shotPar;
+
+typedef struct _raypar { /* ray-tracing parameters */
+ int smoothwindow;
+ int useT2;
+ int geomspread;
+ int nray;
+} rayPar;
+
+#ifndef TRUE
+# define TRUE 1
+#endif
+
+#ifndef FALSE
+# define FALSE 0
+#endif
+
+#define equal(x,y) !strcmp(x,y)
+#define min2(a,b) (((a) < (b)) ? (a) : (b))
+#define max2(a,b) (((a) > (b)) ? (a) : (b))
+
+#define Infinity FLT_MAX
+
+#if __STDC_VERSION__ >= 199901L
+ /* "restrict" is a keyword */
+#else
+#define restrict
+#endif
+
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..9b0c88d6e59fe63c13275cc139f7352bff266501
--- /dev/null
+++ b/raytime3d/raytime3d.c
@@ -0,0 +1,318 @@
+#include <DELPHI_IOc.h>
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" RAYTIME3D - modeling of one-way traveltime for operators in 3D media",
+" ",
+" raytime3d file_vel= xsrc1= zsrc1= ysrc1= [optional parameters]",
+" ",
+" Required parameters:",
+" ",
+" file_vel= ................ gridded velocity file ",
+" 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",
+" ny=1 ..................... if 2D file number of y traces (2D model)",
+" nxmax=512 ................ maximum number of traces in input files",
+" ntmax=1024 ............... maximum number of samples/trace in input files",
+" file_out= ................ output file with traveltime cube",
+" file_amp= ................ output file with approximate amplitudes",
+" RAY TRACING ",
+" 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=0,(nx-1)*dx ......... x-position's of receivers (array)",
+" yrcv=0,(ny-1)*dy ......... y-position's of receivers (array)",
+" zrcv=0,0 ................. z-position's of receivers (array)",
+" dxrcv=dx ................. step in receiver x-direction",
+" dyrcv=dy ................. 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 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);
+
+void main(int argc, char *argv[])
+{
+ 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 delphi velocity file */
+ int32 type, dom1, dom2;
+ int error, n1, n2, ret, size, nkeys, verbose;
+ int ntmax, nxmax;
+ float d1, d2, f1, f2, *tmpdata, c, scl, ox, oz, oy;
+ char *file_vel, *file_out, *key;
+ segyhdr *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) {
+ samess("Hole, J.A., and B.C. Zelt, 1995. \"3-D finite-difference");
+ samess("reflection traveltimes\". Geophys. J. Int., 121, 427-434");
+ }
+ if(!getparstring("file_out",&file_out)) saerr("file_out not given");
+ if(!getparstring("file_vel", &file_vel)) {
+ sawarn("file_vel not defined, assuming homogeneous model");
+ if(!getparfloat("c",&c)) saerr("c not defined");
+ if(!getparint("nx",&nx)) saerr("nx not defined");
+ if(!getparint("ny",&ny)) saerr("ny not defined");
+ if(!getparint("nz",&nz)) saerr("nz not defined");
+ if(!getparfloat("h",&h)) saerr("h not defined");
+ }
+ if(!getparint("ny",&ny)) saerr("ny not defined");
+
+/*---------------------------------------------------------------------------*
+ * Open velocity file
+ *---------------------------------------------------------------------------*/
+
+ if (file_vel != NULL) {
+ 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);
+ if (ret >= 0) {
+ if (!getparint("ntmax", &ntmax)) ntmax = n1;
+ if (!getparint("nxmax", &nxmax)) nxmax = n2;
+ if (verbose>=2 && (ntmax!=n1 || nxmax!=n2))
+ samess("dimensions overruled: %d x %d",ntmax,nxmax);
+ }
+ else {
+ if (!getparint("ntmax", &ntmax)) ntmax = 1024;
+ if (!getparint("nxmax", &nxmax)) nxmax = 512;
+ if (verbose>=2) samess("dimensions used: %d x %d",ntmax,nxmax);
+ }
+ size = ntmax * nxmax;
+ tmpdata = alloc1float(size);
+ hdrs = (segyhdr *) malloc(nxmax*sizeof(segyhdr));
+
+ if (!getparstring("key", &key)) {
+ ret = get_sort(file_vel);
+ if (ret < 0) key = "gy";
+ else key = getkey(ret);
+ }
+ if (verbose) samess("input sorting key is %s",key);
+ set_sukey(key);
+
+ ret = read_data(file_vel,tmpdata,size,&n1,&n2,&f1,&f2,&d1,&d2,
+ &type,hdrs);
+ if (ret < 0) saerr("error in reading data from file %s", file_vel);
+ if (hdrs[0].scalco < 0) scl = 1.0/fabs(hdrs[0].scalco);
+ else if (hdrs[0].scalco == 0) scl = 1.0;
+ else scl = hdrs[0].scalco;
+ get_axis(&dom1, &dom2);
+ if (d1 != d2)
+ saerr("d1 != d2; this is not allowed in the calculation");
+ h = d1;
+ if (dom1 == SA_AXIS_Z) {
+ nx = n2; nz = n1;
+ ox = hdrs[0].gx*scl; oy = hdrs[0].gy*scl; oz = f1;
+ }
+ else {
+ nx = n1; nz = n2;
+ ox = f1; oy = hdrs[0].gy*scl; oz = f1;
+ }
+
+ slow0 = (float *)malloc(ny*n1*n2*sizeof(float));
+ if (slow0 == NULL) saerr("Out of memory for slow0 array!");
+ nxy = nx * ny;
+ if (verbose) samess("h = %.2f nx = %d nz = %d ny = %d", h, nx, nz, ny);
+
+ yy = 0;
+ while (ret >= 0) {
+ if (verbose==2) disp_info(file_vel,n1,n2,f1,f2,d1,d2,type);
+
+ if (dom1 == SA_AXIS_Z) {
+ if (n2 != nx || n1 != nz) saerr("dimensions changed");
+ for (i = 0; i < n2; i++) {
+ for (j = 0; j < n1; j++)
+ slow0[j*nxy+yy*nx+i] = h/tmpdata[i*n1+j];
+ }
+ }
+ else {
+ if (n1 != nx || n2 != nz) saerr("dimensions changed");
+ for (i = 0; i < n2; i++) {
+ for (j = 0; j < n1; j++)
+ slow0[i*nxy+yy*nx+j] = h/tmpdata[i*n1+j];
+ }
+ }
+
+ yy += 1;
+ ret = read_data(file_vel, tmpdata, size, &n1, &n2, &f1, &f2,
+ &d1, &d2, &type, hdrs);
+ }
+ ret = close_file(file_vel);
+ if (ret < 0) sawarn("err %d on closing input file",ret);
+ free1float(tmpdata);
+ if (yy == 1) {
+ if(!getparint("ny",&ny)) samess("2D model defined");
+ else {
+ slow0 = (float *)realloc(slow0, ny*nx*nz*sizeof(float));
+ if (slow0 == NULL) saerr("Out of memory for slow0 array!");
+
+ samess("3D model defined from 2D model");
+ for (zz = 0; zz < nz; zz++) {
+ for (yy = 1; yy < ny; yy++) {
+ for (xx = 0; xx < nx; xx++)
+ slow0[zz*nxy+yy*nx+xx] = slow0[zz*nxy+xx];
+ }
+ }
+ }
+ }
+
+ }
+ else {
+ nxy = nx * ny;
+ slow0 = (float *)malloc(nx*nz*ny*sizeof(float));
+ if (slow0 == NULL) saerr("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) saerr("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) samess("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 ) saerr("error in creating output file %s", file_out);
+
+ hdrs = (segyhdr *) malloc(ny*sizeof(segyhdr));
+ tmpdata = alloc1float(nxy);
+ 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];
+ }
+ }
+
+ 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_out,tmpdata,nx,ny,f1,f2,d1,d2,type,hdrs);
+ if (ret < 0 ) saerr("error on writing output file.");
+ ret = close_file(file_out);
+ if (ret < 0) saerr("err %d on closing output file",ret);
+
+ free(time0);
+ free(slow0);
+ free(hdrs);
+ free1float(tmpdata);
+
+ exit(0);
+
+}
diff --git a/raytime3d/readModel.c b/raytime3d/readModel.c
new file mode 100644
index 0000000000000000000000000000000000000000..27f1da715a50dbc36eac4827efe68434288ee9e7
--- /dev/null
+++ b/raytime3d/readModel.c
@@ -0,0 +1,80 @@
+#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"
+#include "par.h"
+#include "raytime.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))
+
+/**
+* Reads gridded model files and compute from them medium parameters used in the FD kernels.
+* The files read in contain the P (and S) wave velocity and density.
+* The medium parameters calculated are lambda, mu, lambda+2mu, and 1/ro.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+
+int readModel(modPar mod, float *velocity, float *slowness, int nw)
+{
+ FILE *fpcp;
+ size_t nread;
+ int i, tracesToDo, j;
+ int nz, nx;
+ segy hdr;
+
+
+ /* grid size and start positions for the components */
+ nz = mod.nz;
+ nx = mod.nx;
+
+/* open files and read first header */
+
+ fpcp = fopen( mod.file_cp, "r" );
+ assert( fpcp != NULL);
+ nread = fread(&hdr, 1, TRCBYTES, fpcp);
+ assert(nread == TRCBYTES);
+
+/* read all traces */
+
+ tracesToDo = mod.nx;
+ i = 0;
+ while (tracesToDo) {
+ nread = fread(&velocity[i*nz], sizeof(float), hdr.ns, fpcp);
+ assert (nread == hdr.ns);
+ for (j=0;j<nz;j++) {
+ if (velocity[i*nz+j]!=0.0) {
+ slowness[(i+nw)*nz+j+nw] = 1.0/velocity[i*nz+j];
+ }
+ }
+ for (j=0;j<nw;j++) slowness[(i+nw)*nz+j] = slowness[(i+nw)*nz+nw];
+ for (j=nz+nw;j<nz+2*nw;j++) slowness[(i+nw)*nz+j] = slowness[(i+nw)*nz+nz+nw-1];
+
+ nread = fread(&hdr, 1, TRCBYTES, fpcp);
+ if (nread==0) break;
+ i++;
+ }
+ fclose(fpcp);
+
+ for (i=0;i<nw;i++) {
+ for (j=0;j<nz+2*nw;j++) {
+ slowness[(i)*nz+j] = slowness[(nw)*nz+j];
+ slowness[(nx+nw+i)*nz+j] = slowness[(nx+nw-1)*nz+j];
+ }
+ }
+
+ return 0;
+}
+
+
diff --git a/raytime3d/recvPar.c b/raytime3d/recvPar.c
new file mode 100644
index 0000000000000000000000000000000000000000..4ad9eae9be2f3c7bcbeb1ee061dc97d5a4ff9d77
--- /dev/null
+++ b/raytime3d/recvPar.c
@@ -0,0 +1,519 @@
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#include "raytime.h"
+#include "par.h"
+
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* Calculates the receiver positions based on the input parameters
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+*
+* Ammendments:
+* Max Holicki changing the allocation receiver array (2-2016)
+* The Netherlands
+**/
+
+
+void name_ext(char *filename, char *extension);
+
+int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx, int nz)
+{
+ float *xrcv1, *xrcv2, *zrcv1, *zrcv2;
+ int i, ix, ir, verbose;
+ float dxrcv, dzrcv, *dxr, *dzr;
+ float rrcv, dphi, oxrcv, ozrcv, arcv;
+ double circ, h, a, b, e, s, xr, zr, dr, srun, phase;
+ float xrange, zrange, sub_x1, sub_z1;
+ int Nx1, Nx2, Nz1, Nz2, Ndx, Ndz, iarray, nrec, nh;
+ int nxrcv, nzrcv, ncrcv, nrcv, ntrcv, *nlrcv;
+ float *xrcva, *zrcva;
+ char* rcv_txt;
+ FILE *fp;
+
+ if (!getparint("verbose", &verbose)) verbose = 0;
+
+ /* Calculate Model Dimensions */
+ sub_x1=sub_x0+(nx-1)*dx;
+ sub_z1=sub_z0+(nz-1)*dz;
+
+/* Compute how many receivers are defined and then allocate the receiver arrays */
+
+ /* Receivers on a Circle */
+ if (getparfloat("rrcv",&rrcv)) {
+ if (!getparfloat("dphi",&dphi)) dphi=2.0;
+ ncrcv=NINT(360.0/dphi);
+ if (verbose) vmess("Total number of receivers on a circle: %d",ncrcv);
+ }
+ else {
+ ncrcv=0;
+ }
+
+ /* Receivers from a File */
+ ntrcv=0;
+ if (!getparstring("rcv_txt",&rcv_txt)) rcv_txt=NULL;
+ if (rcv_txt!=NULL) {
+ /* Open text file */
+ fp=fopen(rcv_txt,"r");
+ assert(fp!=NULL);
+ /* Get number of lines */
+ while (!feof(fp)) if (fgetc(fp)=='\n') ntrcv++;
+ fseek(fp,-1,SEEK_CUR);
+ if (fgetc(fp)!='\n') ntrcv++; /* Checks if last line terminated by /n */
+ if (verbose) vmess("Number of receivers in rcv_txt file: %d",ntrcv);
+ rewind(fp);
+ }
+
+ /* Receiver Array */
+ nxrcv=countparval("xrcva");
+ nzrcv=countparval("zrcva");
+ if (nxrcv!=nzrcv) verr("Number of receivers in array xrcva (%d), zrcva(%d) are not equal",nxrcv,nzrcv);
+ if (verbose&&nxrcv) vmess("Total number of array receivers: %d",nxrcv);
+
+ /* Linear Receiver Arrays */
+ Nx1 = countparval("xrcv1");
+ Nx2 = countparval("xrcv2");
+ Nz1 = countparval("zrcv1");
+ Nz2 = countparval("zrcv2");
+ if (Nx1!=Nx2) verr("Number of receivers starting points in 'xrcv1' (%d) and number of endpoint in 'xrcv2' (%d) are not equal",Nx1,Nx2);
+ if (Nz1!=Nz2) verr("Number of receivers starting points in 'zrcv1' (%d) and number of endpoint in 'zrcv2' (%d) are not equal",Nz1,Nz2);
+ if (Nx1!=Nz2) verr("Number of receivers starting points in 'xrcv1' (%d) and number of endpoint in 'zrcv2' (%d) are not equal",Nx1,Nz2);
+
+ rec->max_nrec=ncrcv+ntrcv+nxrcv;
+
+ /* no receivers are defined use default linear array of receivers on top of model */
+ if (!rec->max_nrec && Nx1==0) Nx1=1; // Default is to use top of model to record data
+
+ if (Nx1) {
+ /* Allocate Start & End Points of Linear Arrays */
+ xrcv1=(float *)malloc(Nx1*sizeof(float));
+ xrcv2=(float *)malloc(Nx1*sizeof(float));
+ zrcv1=(float *)malloc(Nx1*sizeof(float));
+ zrcv2=(float *)malloc(Nx1*sizeof(float));
+ if (!getparfloat("xrcv1",xrcv1)) xrcv1[0]=sub_x0;
+ if (!getparfloat("xrcv2",xrcv2)) xrcv2[0]=sub_x1;
+ if (!getparfloat("zrcv1",zrcv1)) zrcv1[0]=sub_z0;
+ if (!getparfloat("zrcv2",zrcv2)) zrcv2[0]=zrcv1[0];
+
+ /* check if receiver arrays fit into model */
+ for (iarray=0; iarray<Nx1; iarray++) {
+ xrcv1[iarray] = MAX(sub_x0, xrcv1[iarray]);
+ xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
+ xrcv2[iarray] = MAX(sub_x0, xrcv2[iarray]);
+ xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
+
+ zrcv1[iarray] = MAX(sub_z0, zrcv1[iarray]);
+ zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
+ zrcv2[iarray] = MAX(sub_z0, zrcv2[iarray]);
+ zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
+ }
+
+ /* Crop to Fit Model */
+/* Max's addtion still have to check if it has the same fucntionality */
+ for (iarray=0;iarray<Nx1;iarray++) {
+ if (xrcv1[iarray]<sub_x0) {
+ if (xrcv2[iarray]<sub_x0) {
+ verr("Linear array %d outside model bounds",iarray);
+ }
+ else {
+ vwarn("Cropping element %d of 'xrcv1' (%f) to model bounds (%f)",iarray,xrcv1[iarray],sub_x0);
+ xrcv1[iarray]=sub_x0;
+ }
+ }
+ else if (xrcv1[iarray] > sub_x1) {
+ verr("Linear array %d outside model bounds",iarray);
+ }
+ if ( (xrcv2[iarray] < xrcv1[iarray]) ) {
+ verr("Ill defined linear array %d, 'xrcv1' (%f) greater than 'xrcv2' (%f)",iarray,xrcv1[iarray],xrcv2[iarray]);
+ }
+ else if (xrcv2[iarray]>sub_x1) {
+ vwarn("Cropping element %d of 'xrcv2' (%f) to model bounds (%f)",iarray,xrcv2[iarray],sub_x1);
+ xrcv2[iarray]=sub_x1;
+ }
+
+ if (zrcv1[iarray] < sub_z0) {
+ if (zrcv2[iarray] < sub_z0) {
+ verr("Linear array %d outside model bounds",iarray);
+ }
+ else {
+ vwarn("Cropping element %d of 'zrcv1' (%f) to model bounds (%f)",iarray,zrcv1[iarray],sub_z0);
+ zrcv1[iarray]=sub_z0;
+ }
+ }
+ else if (zrcv1[iarray] > sub_z1) {
+ verr("Linear array %d outside model bounds",iarray);
+ }
+ if ( (zrcv2[iarray] < zrcv1[iarray]) ) {
+ verr("Ill defined linear array %d, 'zrcv1' (%f) greater than 'zrcv2' (%f)",iarray,zrcv1[iarray],zrcv2[iarray]);
+ }
+ else if (zrcv2[iarray]>sub_z1) {
+ vwarn("Cropping element %d of 'xrcv2' (%f) to model bounds (%f)",iarray,zrcv2[iarray],sub_z1);
+ zrcv2[iarray]=sub_z1;
+ }
+ }
+
+ /* Get Sampling Rates */
+ Ndx = countparval("dxrcv");
+ Ndz = countparval("dzrcv");
+
+ dxr = (float *)malloc(Nx1*sizeof(float));
+ dzr = (float *)malloc(Nx1*sizeof(float));
+ if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
+ if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
+ if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
+ for (i=1; i<Nx1; i++) {
+ dxr[i] = dxr[0];
+ dzr[i] = dzr[0];
+ }
+ Ndx=1;
+ Ndz=1;
+ }
+ else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
+ for (i=1; i<Nx1; i++) {
+ dxr[i] = dxr[0];
+ dzr[i] = dzr[0];
+ }
+ Ndz=1;
+ Ndx=1;
+ }
+ else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
+ if (Ndx!=Ndz) {
+ verr("Number of 'dxrcv' (%d) is not equal to number of 'dzrcv' (%d) or 1",Ndx,Ndz);
+ }
+ if (Ndx!=Nx1 && Ndx!=1) {
+ verr("Number of 'dxrcv' (%d) is not equal to number of starting points in 'xrcv1' (%d) or 1",Ndx,Nx1);
+ }
+ }
+
+ /* check consistency of receiver steps */
+ for (iarray=0; iarray<Ndx; iarray++) {
+ if (dxr[iarray]<0) {
+ dxr[i]=dx;
+ vwarn("'dxrcv' element %d (%f) is less than zero, changing it to %f'",iarray,dxr[iarray],dx);
+ }
+ }
+ for (iarray=0;iarray<Ndz;iarray++) {
+ if (dzr[iarray]<0) {
+ dzr[iarray]=dz;
+ vwarn("'dzrcv' element %d (%f) is less than zero, changing it to %f'",iarray,dzr[iarray],dz);
+ }
+ }
+ for (iarray=0;iarray<Ndx;iarray++){
+ if (dxr[iarray]==0 && dzr[iarray]==0) {
+ xrcv2[iarray]=xrcv1[iarray];
+ dxr[iarray]=1.;
+ vwarn("'dxrcv' element %d & 'dzrcv' element 1 are both 0.",iarray+1);
+ vmess("Placing 1 receiver at (%d,%d)",xrcv1[iarray],zrcv1[iarray]);
+ }
+ }
+ for (iarray=0;iarray<Ndx;iarray++){
+ if (xrcv1[iarray]==xrcv2[iarray] && dxr[iarray]!=0) {
+ dxr[iarray]=0.;
+ vwarn("Linear array %d: 'xrcv1'='xrcv2' and 'dxrcv' is not 0. Setting 'dxrcv'=0",iarray+1);
+ }
+ }
+ for (iarray=0;iarray<Ndx;iarray++){
+ if (zrcv1[iarray]==zrcv2[iarray] && dzr[iarray]!=0.){
+ dzr[iarray]=0.;
+ vwarn("Linear array %d: 'zrcv1'='zrcv2' and 'dzrcv' is not 0. Setting 'dzrcv'=0",iarray+1);
+ }
+ }
+
+ /* Calculate Number of Receivers */
+ nrcv = 0;
+ nlrcv=(int *)malloc(Nx1*sizeof(int));
+ for (iarray=0; iarray<Nx1; iarray++) {
+ xrange = (xrcv2[iarray]-xrcv1[iarray]);
+ zrange = (zrcv2[iarray]-zrcv1[iarray]);
+ if (dxr[iarray] != 0.0) {
+ nlrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
+ }
+ else {
+ if (dzr[iarray] == 0) {
+ verr("For receiver array %d: receiver distance dzrcv is not given", iarray);
+ }
+ nlrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
+ }
+ nrcv+=nlrcv[iarray];
+ }
+
+ /* Calculate Number of Receivers */
+/*
+ nlrcv=(int *)malloc(Nx1*sizeof(int));
+ if (!isnan(*xrcv1)) *nlrcv=MIN(NINT((*xrcv2-*xrcv1)/(*dxr)),NINT((*zrcv2-*zrcv1)/(*dzr)))+1;
+ else *nlrcv=0;
+ nrcv=*nlrcv;
+ if (verbose>4 && nlrcv[iarray]!=0) vmess("Linear receiver array 1 has final bounds: (X: %f -> %f,Z: %f ->
+%f)",xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*(*dxr),zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*(*dzr));
+ if (Ndx>1) {
+ for (iarray=1;iarray<Nx1;iarray++) {
+ if (!isnan(xrcv1[iarray])) {
+ nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/dxr[iarray]),NINT((zrcv2[iarray]-zrcv1[iarray])/dzr[iarray]))+1;
+ }
+ else {
+ nlrcv[iarray]=0;
+ }
+ nrcv+=nlrcv[iarray];
+ if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %d has final bounds: (X: %f -> %f,Z: %f ->
+%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*dxr[iarray],zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*dzr[iarray]);
+ }
+ }
+ else {
+ for (iarray=1;iarray<Nx1;iarray++) {
+ if (!isnan(xrcv1[iarray])) nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/(*dxr)),NINT((zrcv2[iarray]-zrcv1[iarray])/(*dzr)))+1;
+ else nlrcv[iarray]=0;
+ nrcv+=nlrcv[iarray];
+ if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %d has final bounds: (X: %f -> %f,Z: %f ->
+%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]**dxr,zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]**dzr);
+ }
+ }
+*/
+ if (verbose) vmess("Total number of linear array receivers: %d",nrcv);
+ if (!nrcv) {
+ free(xrcv1);
+ free(xrcv2);
+ free(zrcv1);
+ free(zrcv2);
+ free(dxr);
+ free(dzr);
+ free(nlrcv);
+ }
+ rec->max_nrec+=nrcv;
+ }
+ else {
+ nrcv=0;
+ }
+
+/* allocate the receiver arrays */
+
+ /* Total Number of Receivers */
+ if (verbose) vmess("Total number of receivers: %d",rec->max_nrec);
+
+ /* Allocate Arrays */
+ rec->x = (int *)calloc(rec->max_nrec,sizeof(int));
+ rec->z = (int *)calloc(rec->max_nrec,sizeof(int));
+ rec->xr = (float *)calloc(rec->max_nrec,sizeof(float));
+ rec->zr = (float *)calloc(rec->max_nrec,sizeof(float));
+
+/* read in the receiver postions */
+
+ nrec=0;
+ /* Receivers on a Circle */
+ if (ncrcv) {
+ if (!getparfloat("oxrcv",&oxrcv)) oxrcv=0.0;
+ if (!getparfloat("ozrcv",&ozrcv)) ozrcv=0.0;
+ if (!getparfloat("arcv",&arcv)) {
+ arcv=rrcv;
+ for (ix=0; ix<ncrcv; ix++) {
+ rec->xr[ix] = oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI));
+ rec->zr[ix] = ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI));
+ rec->x[ix] = NINT(rec->xr[ix]/dx);
+ rec->z[ix] = NINT(rec->zr[ix]/dz);
+ //rec->x[ix] = NINT((oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI)))/dx);
+ //rec->z[ix] = NINT((ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI)))/dz);
+ if (verbose>4) fprintf(stderr,"Receiver Circle: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
+ }
+ }
+ else { /* an ellipse */
+ /* simple numerical solution to find equidistant points on an ellipse */
+ nh = (ncrcv)*1000; /* should be fine enough for most configurations */
+ h = 2.0*M_PI/nh;
+ a = MAX(rrcv, arcv);
+ b = MIN(rrcv, arcv);
+ e = sqrt(a*a-b*b)/a;
+ //fprintf(stderr,"a=%f b=%f e=%f\n", a, b, e);
+ circ = 0.0;
+ for (ir=0; ir<nh; ir++) {
+ s = sin(ir*h);
+ circ += sqrt(1.0-e*e*s*s);
+ }
+ circ = a*h*circ;
+ //fprintf(stderr,"circ = %f circle=%f\n", circ, 2.0*M_PI*rrcv);
+ /* define distance between receivers on ellipse */
+ dr = circ/ncrcv;
+ ix = 0;
+ srun = 0.0;
+ if (arcv >= rrcv) phase=0.0;
+ else phase=0.5*M_PI;
+ for (ir=0; ir<nh; ir++) {
+ s = sin(ir*h);
+ srun += sqrt(1.0-e*e*s*s);
+ if (a*h*srun >= ix*dr ) {
+ xr = rrcv*cos(ir*h+phase);
+ zr = arcv*sin(ir*h+phase);
+ rec->xr[ix] = oxrcv-sub_x0+xr;
+ rec->zr[ix] = ozrcv-sub_z0+zr;
+ rec->x[ix] = NINT(rec->xr[ix]/dx);
+ rec->z[ix] = NINT(rec->zr[ix]/dz);
+ if (verbose>4) fprintf(stderr,"Receiver Ellipse: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
+ ix++;
+ }
+ if (ix == ncrcv) break;
+ }
+ }
+
+ /* check if receivers fit into the model otherwise clip to edges */
+ for (ix=0; ix<ncrcv; ix++) {
+ rec->x[ix] = MIN(nx-1, MAX(rec->x[ix], 0));
+ rec->z[ix] = MIN(nz-1, MAX(rec->z[ix], 0));
+ }
+ nrec += ncrcv;
+ }
+
+ /* Receiver Text File */
+
+ if (ntrcv) {
+ /* Allocate arrays */
+ xrcva = (float *)malloc(nrcv*sizeof(float));
+ zrcva = (float *)malloc(nrcv*sizeof(float));
+ /* Read in receiver coordinates */
+ for (i=0;i<nrcv;i++) {
+ if (fscanf(fp,"%e %e\n",&xrcva[i],&zrcva[i])!=2) vmess("Receiver Text File: Can not parse coordinates on line %d.",i);
+ }
+ /* Close file */
+ fclose(fp);
+ /* Process coordinates */
+ for (ix=0; ix<nrcv; ix++) {
+ rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
+ rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
+ rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
+ rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
+ if (verbose>4) vmess("Receiver Text Array: xrcv[%d]=%f zrcv=%f", ix, rec->xr[nrec+ix]+sub_x0, rec->zr[nrec+ix]+sub_z0);
+ }
+ free(xrcva);
+ free(zrcva);
+ nrec += ntrcv;
+ }
+
+ /* Receiver Array */
+ if (nxrcv != 0) {
+ /* receiver array is defined */
+ xrcva = (float *)malloc(nxrcv*sizeof(float));
+ zrcva = (float *)malloc(nxrcv*sizeof(float));
+ getparfloat("xrcva", xrcva);
+ getparfloat("zrcva", zrcva);
+ for (ix=0; ix<nxrcv; ix++) {
+ rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
+ rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
+ rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
+ rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
+ if (verbose>4) fprintf(stderr,"Receiver Array: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[nrec+ix]+sub_x0, rec->zr[nrec+ix]+sub_z0);
+ }
+ free(xrcva);
+ free(zrcva);
+ nrec += nxrcv;
+ }
+
+ /* Linear Receiver Arrays */
+ if (nrcv!=0) {
+ xrcv1 = (float *)malloc(Nx1*sizeof(float));
+ xrcv2 = (float *)malloc(Nx1*sizeof(float));
+ zrcv1 = (float *)malloc(Nx1*sizeof(float));
+ zrcv2 = (float *)malloc(Nx1*sizeof(float));
+
+ if(!getparfloat("xrcv1", xrcv1)) xrcv1[0]=sub_x0;
+ if(!getparfloat("xrcv2", xrcv2)) xrcv2[0]=(nx-1)*dx+sub_x0;
+ if(!getparfloat("zrcv1", zrcv1)) zrcv1[0]=sub_z0;
+ if(!getparfloat("zrcv2", zrcv2)) zrcv2[0]=zrcv1[0];
+
+ Ndx = countparval("dxrcv");
+ Ndz = countparval("dzrcv");
+
+ dxr = (float *)malloc(Nx1*sizeof(float));
+ dzr = (float *)malloc(Nx1*sizeof(float));
+ if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
+ if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
+ if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
+ for (i=1; i<Nx1; i++) {
+ dxr[i] = dxr[0];
+ dzr[i] = dzr[0];
+ }
+ Ndx=1;
+ }
+ else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
+ for (i=1; i<Nx1; i++) {
+ dxr[i] = dxr[0];
+ dzr[i] = dzr[0];
+ }
+ Ndz=1;
+ }
+ else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
+ if (Ndx>1) assert(Ndx==Nx1);
+ if (Ndz>1) assert(Ndz==Nx1);
+ }
+
+/*
+ if ( (Ndx!=0) && (Ndz!=0) ) {
+ vwarn("Both dzrcv and dxrcv are set: dxrcv value is used");
+ Ndz=0;
+ for (i=0; i<Nx1; i++) dzr[i] = 0.0;
+ }
+*/
+ /* check if receiver arrays fit into model */
+ for (iarray=0; iarray<Nx1; iarray++) {
+ xrcv1[iarray] = MAX(sub_x0, xrcv1[iarray]);
+ xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
+ xrcv2[iarray] = MAX(sub_x0, xrcv2[iarray]);
+ xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
+
+ zrcv1[iarray] = MAX(sub_z0, zrcv1[iarray]);
+ zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
+ zrcv2[iarray] = MAX(sub_z0, zrcv2[iarray]);
+ zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
+ }
+
+ /* calculate receiver array and store into rec->x,z */
+
+ for (iarray=0; iarray<Nx1; iarray++) {
+ xrange = (xrcv2[iarray]-xrcv1[iarray]);
+ zrange = (zrcv2[iarray]-zrcv1[iarray]);
+ if (dxr[iarray] != 0.0) {
+ nrcv = nlrcv[iarray];
+ dxrcv=dxr[iarray];
+ dzrcv = zrange/(nrcv-1);
+ if (dzrcv != dzr[iarray]) {
+ vwarn("For receiver array %d: calculated dzrcv=%f given=%f", iarray, dzrcv, dzr[iarray]);
+ vwarn("The calculated receiver distance %f is used", dzrcv);
+ }
+ }
+ else {
+ if (dzr[iarray] == 0) {
+ verr("For receiver array %d: receiver distance dzrcv is not given", iarray);
+ }
+ nrcv = nlrcv[iarray];
+ dxrcv = xrange/(nrcv-1);
+ dzrcv = dzr[iarray];
+ if (dxrcv != dxr[iarray]) {
+ vwarn("For receiver array %d: calculated dxrcv=%f given=%f", iarray, dxrcv, dxr[iarray]);
+ vwarn("The calculated receiver distance %f is used", dxrcv);
+ }
+ }
+
+ // calculate coordinates
+ for (ir=0; ir<nrcv; ir++) {
+ rec->xr[nrec]=xrcv1[iarray]-sub_x0+ir*dxrcv;
+ rec->zr[nrec]=zrcv1[iarray]-sub_z0+ir*dzrcv;
+
+ rec->x[nrec]=NINT((rec->xr[nrec])/dx);
+ rec->z[nrec]=NINT((rec->zr[nrec])/dz);
+ nrec++;
+ }
+ }
+ free(xrcv1);
+ free(xrcv2);
+ free(zrcv1);
+ free(zrcv2);
+ free(dxr);
+ free(dzr);
+ free(nlrcv);
+ }
+
+ rec->n=rec->max_nrec;
+ return 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..1dadfa5746b00a8bb5504f16de26188c06da5c97
--- /dev/null
+++ b/raytime3d/src3d.c
@@ -0,0 +1,208 @@
+#include <cwp.h>
+#include <su.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 saerrpkge.c */
+extern void saerr(char *fmt, ...);
+extern void sawarn(char *fmt, ...);
+extern void samess(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)) saerr("xsrc1 not given");
+ if(!getparfloat("ysrc1",&ysrc1)) saerr("ysrc1 not given");
+ if(!getparfloat("zsrc1",&zsrc1)) saerr("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){
+ sawarn("Source near an edge, beware of traveltime errors");
+ sawarn("for raypaths that travel parallel to edge ");
+ sawarn("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)) saerr("srcwall not given");
+ if (!getparstring("wallfile",&wallfile)) saerr("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)) saerr("srcwall not given");
+ if (!getparstring("oldtfile",&oldtfile)) saerr("oldtfile not given");
+ if((ofint=open(oldtfile,O_RDONLY,0664))<=1) {
+ fprintf(stderr,"cannot open %s\n",oldtfile);
+ exit(-1);
+ }
+ }
+ else {
+ saerr("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=PI-theta1;
+ if (rxyc<rxy1) theta2=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..06b4b388274def5752a95ab1f0bbb1b881b47c11
--- /dev/null
+++ b/raytime3d/vidale3d.c
@@ -0,0 +1,1440 @@
+#include <cwp.h>
+#include <su.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 saerrpkge.c */
+extern void saerr(char *fmt, ...);
+extern void sawarn(char *fmt, ...);
+extern void samess(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;
+ sawarn("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)
+ saerr("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)) saerr("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) samess("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) samess("Head waves found on left: %d",headw[1]);
+ if (headw[1]>head) {
+ head = headw[1];
+ srcwall = 1;
+ }
+ }
+ if (headw[2]>0) {
+ if(verbose) samess("Head waves found on right: %d",headw[2]);
+ if (headw[2]>head) {
+ head = headw[2];
+ srcwall = 2;
+ }
+ }
+ if (headw[3]>0) {
+ if(verbose) samess("Head waves found on front: %d",headw[3]);
+ if (headw[3]>head) {
+ head = headw[3];
+ srcwall = 3;
+ }
+ }
+ if (headw[4]>0) {
+ if(verbose) samess("Head waves found on back: %d",headw[4]);
+ if (headw[4]>head) {
+ head = headw[4];
+ srcwall = 4;
+ }
+ }
+ if (headw[5]>0) {
+ if(verbose) samess("Head waves found on top: %d",headw[5]);
+ if (headw[5]>head) {
+ head = headw[5];
+ srcwall = 5;
+ }
+ }
+ if (headw[6]>0) {
+ if(verbose) samess("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)
+ samess("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;
+ samess("RESTART at left side of model"); }
+ else { x2=nx; dx2=0; }
+ if (srcwall == 2) { x1=nx-2;
+ samess("RESTART at right side of model"); }
+ else { x1= -1; dx1=0; }
+ if (srcwall == 3) { y2=1;
+ samess("RESTART at front side of model"); }
+ else { y2=ny; dy2=0; }
+ if (srcwall == 4) { y1=ny-2;
+ samess("RESTART at back side of model"); }
+ else { y1= -1; dy1=0; }
+ if (srcwall == 5) { z2=1;
+ samess("RESTART at top side of model"); }
+ else { z2=nz; dz2=0; }
+ if (srcwall == 6) { z1=nz-2;
+ samess("RESTART at bottom side of model"); }
+ else { z1= -1; dz1=0; }
+ if (reverse == 0)
+ sawarn("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..faf297811e469314c8de0265f78b87543b27af9b
--- /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"raytime.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..6eac57d300aaa653e621b99fe2731c9ed3ddb60a
--- /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 "raytime.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;
+}
+