time_3d.c

射线追踪程序

/*----------------------------------------------------------------------------*/
/*  TIME_3D: FINITE DIFFERENCE COMPUTATION OF FIRST ARRIVAL TIMES IN 3D.      */
/*----------------------------------------------------------------------------*/
/*  P.PODVIN, Geophysique, Ecole des Mines de Paris, Fontainebleau.           */
/*  E-MAIL: [e-mail]          Tel: [phone].          */
/*  April 1991, previous revision: 7 May 1993                                 */
/*  This "robustified and ANSIfied" version is dated 12 July 2003             */
/*                                                                            */
/*  PROTOTYPE :                                                               */
/*    int time_3d(float *HS, float *T, int NX, int NY, int NZ,                */
/*                float XS, float YS, float ZS, float HS_EPS_INIT, int MSG);  */
/*                                                                            */
/*  ARGUMENTS: (C description. All FORTRAN arguments are pointers)            */
/*    (int)      NX,NY,NZ : Dimensions of the timefield. Cells are cubic.     */
/*                          Real model dimensions are NX-1,NY-1,NZ-1.         */
/*    (float *)      HS,T : 1-D arrays. Length: NX*NY*NZ                      */
/*                          HS[] is interpreted as slowness*mesh parameter,   */
/*                          so that associated physical dimension is Time.    */
/*                          Both arrays must be arranged as a succession of   */
/*                          planes x (resp. z)=Ct, each plane consisting of   */
/*                          a suite of y=Ct vectors when this program is      */
/*                          called from a program written in C (resp.FORTRAN).*/
/*                          hs and t are 3D pointers allocated in order       */
/*                          to address points and cells straightforwardly.    */
/*                          Implicitly staggered grids are used: t[][][]      */
/*                          represent arrival times at grid-points, whereas   */
/*                          hs[][][] describe (constant) slownesses in cells. */
/*                          Cells such that x=NX-1 or y=NY-1 or z=NZ-1 are    */
/*                          thus dummy cells located OUT of the model.        */
/*                          Their associated slowness is treated as INFINITY  */
/*                          (absorbing boundary conditions).                  */
/*                          Negative hs values are unphysical and forbidden.  */
/*    (float)    XS,YS,ZS : real point source coordinates, i.e. distances     */
/*                          along axes between source point and corner 0,0,0  */
/*                          of the timefield, measured in h units.            */
/*                          i.e. 0.0<=XS<=(float)NX                           */
/*                          If source coordinates are illicit (out of model   */
/*                          bounds), the program uses T as a pre-initialized  */
/*                          timefield (multiple source, e.g. exploding        */
/*                          reflector). At least one element of T must then   */
/*                          be finite.                                        */
/*    (float) HS_EPS_INIT : fraction (typically 1.0E-3) defining the toler-   */
/*                          ated model inhomogeneity for exact initialization.*/
/*                          A tolerance larger than 0.01 will potentially     */
/*                          create errors larger than those involved by the   */
/*                          F.D. scheme without any exact initialization.     */
/*                          (this arg. is used only in the point source case) */
/*    (int)           MSG : Message flag (0:silent,1:few messages,2:verbose)  */
/*                          A negative value inhibits "clever" initialization.*/
/*                                                                            */
/*  VALUE: time_3d() returns a nonzero integer if an error occurred.          */
/*      In such case, an explicit message is printed on 'stderr'.             */
/*                                                                            */
/*  FORTRAN INTERFACE: see function time_3d_(). Standard FORTRAN conventions  */
/*      are adopted: indexes start at 1 (but this is irrelevant here!),       */
/*      and arrays are transposed to fit FORTRAN-style memory mapping.        */
/*      Please note that time_3d is an INTEGER FUNCTION which name DOES       */
/*      NOT BEGIN WITH [I,J,K,L,M,N] ! You must declare it explicitly !       */
/*      e.g.: MAIN program                                                    */
/*            parameter (MAXNODES=...)                                        */
/*            real*4 hs(MAXNODES),t(MAXNODES)                                 */
/*            integer*4 time_3d                                               */
/*            nx=... ny=... nz=...                        ! #grid-points      */
/*            xs=...(0<=xs<=nx-1) ys=... zs=...           ! source point      */
/*            eps_hs=0.001                                ! tolerance         */
/*            call build_model(hs,nx,ny,nz)                                   */
/*            itime_3d_status=time_3d(hs,t,nx,ny,nz,xs,ys,zs,eps_hs,msg)      */
/*            if(itime_3d_status.ne.0) ----> an error occurred !              */
/*            ...                                                             */
/*            stop                                                            */
/*            end                                                             */
/*                                                                            */
/*            SUBROUTINE build_model(hs,nx,ny,nz)                             */
/*            real*4 hs(nx,ny,nz)                                             */
/*            nested loops: k=...(1<=k<=nz-1) j=... i=... ! cell coordinates  */
/*                  hs(i,j,k)=...                                             */
/*            end loops                                                       */
/*            return                                                          */
/*            end                                                             */
/*                                                                            */
/*      With SunOS, the FORTRAN compiler appends the final '_' to FORTRAN     */
/*      function name 'time_3d'. Full name 'time_3d_' may be necessary on     */
/*      other systems. Remember that NO standard C/Fortran interface exists ! */
/*                                                                            */
/*  PREPROCESSOR OPTIONS:       (/bin/cpp options applicable at compile time) */
/*      -DNO_IEEE_PROTOCOL : for systems where IEEE standards are not defined */
/*      -DINIT_MIN=value   :         these two options control initialization */
/*      -DINIT_RECURS_LIMIT=value     (see comments after #ifndef directives) */
/*      -DDEBUG_ARGS       (to check what function you called with what args) */
/*                                                                            */
/*  COMMON ERRORS:                                                            */
/*       Incomprehensible behaviours in general result from C/Fortran         */
/*       interface problems . Most frequent sources of confusing errors :     */
/*       -incorrect argument types (NX MUST be an integer, XS a real number)  */
/*       -incorrect initialization of array T when using a multiple source    */
/*       (e.g. T being initialized with zeroes,if source coordinates are out  */
/*        of model boundaries, then array T will be left unchanged ! This is  */
/*        now checked as an error condition)                                  */
/*                                                                            */
/*  LAST CHANGES ("robustification") :                                        */
/*    This routine has be used intensively by many people in the past ten     */
/*    years, but recently, I was informed of incomprehensible problems        */
/*    encountered pseudo-randomly (and very rarely) by some users...          */
/*    Here is what I found out and the solution I designed :                  */
/*    Everyone should know (;=)) that no program should critically depend     */
/*    on "exact" comparisons between floating point numbers, because these    */
/*    numbers are NOT represented exactly on a computer. I only recently      */
/*    became aware of the fact that this was a serious (and vicious) flaw     */
/*    in this source code, that only reveals on some platforms (Intel-based   */
/*    PCs in particular !) Technically, recursive calls in this program are   */
/*    triggered by tests of the form (a<b) and, with no precaution, such      */
/*    tests may be evaluated as true when a and b are evaluated with exactly  */
/*    the same set of instructions yet in two different parts of the source   */
/*    code. Clearly, a and b should then be exactly equal, but... they are    */
/*    not and this results into infinite recursions crashing the program !    */
/*    The problem is illustrated and explained as follows in gcc.info :       */
/*        "On 68000 and x86 systems, for instance, you can get paradoxical    */
/*        results if you test the precise values of floating point numbers.   */
/*        For example, you can find that a floating point value which is not  */
/*        a NaN is not equal to itself.  This results from the fact that the  */
/*        floating point registers hold a few more bits of precision than     */
/*        fit in a `double' in memory.  Compiled code moves values between    */
/*        memory and floating point registers at its convenience, and moving  */
/*        them into memory truncates them."                                   */
/*    In this "robustified" version, all critical tests (the ones triggering  */
/*    recursive calls) have been "fuzzified" (i.e., a<b became b-a>eps) to    */
/*    avoid this pathological behaviour. Because this depends on floating     */
/*    point encoding, "eps" also is. It is in fact evaluated in accordance    */
/*    with the imprecision of standard encoding schemes such as IEEE754.      */
/*    Fuzzification is controlled by the macro parameter EPS_FUZZY.           */
/*    If you get into trouble (i.e., encounter cases where time_3d crashes    */
/*    after having entered an infinite loop of recursive calls), you might    */
/*    consider changing this macro to a higher value, but PLEASE, be kind     */
/*    enough to also inform me by e-mail, thanks!                             */
/*                                                                            */
/*  REFERENCE: Podvin & Lecomte, Geophys.J.Intern. 105, 271-284, 1991.        */
/*----------------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#ifndef M_SQRT2
#define M_SQRT2     1.41421356237309504880
#endif
#define M_SQRT3     1.732050807568877076

#ifdef NO_IEEE_PROTOCOL
#define INFINITY    1.0e+19
    /* INFINITY should be < sqrt(MAX_FLOAT_VALUE), machine-dependent */
#define ISINF(x)    ((x)>1.0e+18)
#define NINT(x)     (int)floor((x)+0.5)
    /* NINT = "Nearest INTeger", used here only with positive numbers */
#else
#ifdef sun
#ifdef BSD
/* Note : do NOT use option -Xs with Sun cc as this option UNDEFINES "sun" ! */
#include <sunmath.h>
#endif
#endif
#define INFINITY    (float)infinity()
#define ISINF(x)    isinf(x)
#define NINT(x)     nint(x)
#endif /* NO_IEEE_PROTOCOL */

#define min(x,y) (((x)<(y)) ? (x):(y))
#define max(x,y) (((x)>(y)) ? (x):(y))
#define min3(x,y,z) (min(x,min(y,z)))
#define min4(x,y,z,t) (min(x,min(y,min(z,t))))

#include <malloc.h>

#define EPS_FUZZY	1.2e-07
/* slightly more than 1/2^23 (4-byte-float mantissas are encoded on 23 bits) */

/*-------------------------------------Static functions-----------------------*/

static int
    pre_init(void),
    init_point(void),
    recursive_init(void),
    propagate_point(int),
    x_side(int,int,int,int,int,int),
    y_side(int,int,int,int,int,int),
    z_side(int,int,int,int,int,int),
    scan_x_ff(int,int,int,int,int,int,int),
    scan_x_fb(int,int,int,int,int,int,int),
    scan_x_bf(int,int,int,int,int,int,int),
    scan_x_bb(int,int,int,int,int,int,int),
    scan_y_ff(int,int,int,int,int,int,int),
    scan_y_fb(int,int,int,int,int,int,int),
    scan_y_bf(int,int,int,int,int,int,int),
    scan_y_bb(int,int,int,int,int,int,int),
    scan_z_ff(int,int,int,int,int,int,int),
    scan_z_fb(int,int,int,int,int,int,int),
    scan_z_bf(int,int,int,int,int,int,int),
    scan_z_bb(int,int,int,int,int,int,int);
    /* the only fully commented "side" functions are x_side() ans scan_x_ff() */

static void
    error(int),
    init_nearest(void),
    init_cell(float,float,float,int,int,int),
    free_ptrs(int);

static float
    exact_delay(float,float,float,int,int,int);

static int
    t_1d(int,int,int,float,float,float,float,float),
    t_2d(int,int,int,float,float,float,float),
    diff_2d(int,int,int,float,float,float),
    t_3d_(int,int,int,float,float,float,float,float,int),
    t_3d_part1(int,int,int,float,float,float,float),
    point_diff(int,int,int,float,float),
    edge_diff(int,int,int,float,float,float);

/*-------------------------------------Static variables-----------------------*/

/* MODEL */

static int
    nmesh_x,nmesh_y,nmesh_z;          /* Model dimensions (cells) */
static float
    ***hs,*hs_buf,                    /* 1D and 3D arrays */
    *hs_keep=(float *)NULL;           /* to save boundary values */

/* TIMEFIELD */

static int
    nx,ny,nz;                         /* Timefield dimensions (nodes) */
static float
    ***t,*t_buf;                      /* 1D and 3D arrays */
static float
    timeshift;                        /* used to render all times >=0.0 */

/* SOURCE */

static float
    fxs,fys,fzs;                      /* Point source coordinates */
static int
    xs,ys,zs;                         /* Nearest node */
static int
    mult=0;                           /* Flag used for multiple source */

/* PARAMETERS */

#ifndef INIT_MIN
#define INIT_MIN             7
#endif /* INIT_MIN */
#define N_INIT_X    (4*INIT_MIN+3)
#define N_INIT      N_INIT_X*N_INIT_X*N_INIT_X
                                      /* This conventional value defines  */
                                      /* the maximum size of the box that */
                                      /* will be initialized recursively. */
                                      /* (ADJUSTABLE at compile time).    */
                                      /* Default value is reasonable.     */
                                      /* Cost is already high: 3584 to    */
                                      /* 26416 more points according to   */
                                      /* source position.                 */

#ifndef INIT_RECURS_LIMIT
#define INIT_RECURS_LIMIT    1
#endif /* INIT_RECURS_LIMIT */
                                      /* This parameter defines the maximal */
                                      /* level of recursivity during init.  */
                                      /* (ADJUSTABLE at compile time).      */
                                      /* Value zero would practically rest- */
                                      /* rain initialization to the source  */
                                      /* point in heterogeneous models.     */
                                      /* Value 2 is advisable only when     */
                                      /* VERY severe heterogeneities are    */
                                      /* located close to the source point. */

static int
    messages,                         /* message flag (0:silent)              */
    source_at_node=0,                 /* are source coordinate int's ? (0/1)  */
    no_init=0,                        /* 1: inhibition of "clever" init.      */
    init_stage=0,                     /* level of recursivity during init.    */
    current_side_limit,               /* actual boundary of computations      */
    X0,X1,Y0,Y1,Z0,Z1,                /* inclusive boundaries of timed region */
    sum_updated,                      /* total count of adopted FD stencils   */
    reverse_order,                    /* level of recursivity in FD scheme    */
    *longflags,                       /* local headwave flags.                */
    flag_fb,x_start_fb,y_start_fb,z_start_fb,
    flag_bf,x_start_bf,y_start_bf,z_start_bf,
    flag_ff,x_start_ff,y_start_ff,z_start_ff,
    flag_bb,x_start_bb,y_start_bb,z_start_bb;
                                      /* control current side scanning.       */

static float
    hs_eps_init;                      /* tolerance on homogeneity
                                       (fraction of slowness at source point) */

#define SMALLTALK        messages
#define VERBOSE          messages==2
/*------------------------------------------------Error flags---------------*/

#define NO_ERROR         0
#define ERR_MULTINF      (-1)
#define ERR_INF          (-2)
#define ERR_RECURS       (-3)
#define ERR_MALLOC       (-4)
#define ERR_HS_EPS       (-5)
#define ERR_NONPHYSICAL  (-6)
#define ERR_MULTNONE     (-7)

static char *err_msg[]=
            {
                "\ntime_3d: Computations terminated normally.\n",
                "\ntime_3d: Multiple source but no source at finite time.\n",
                "\ntime_3d: Source point is in a zero velocity zone.\n",
                "\ntime_3d: Fatal error during recursive init.\n",
                "\ntime_3d: Memory Allocation failed.\n",
                "\ntime_3d: Init: Illegal tolerance on inhomogeneity.\n",
                "\ntime_3d: Illegal negative slowness value.\n",
                "\ntime_3d: Multiple source but input time map is uniform.\n"
            };

/*-------------------------------------------------Error()------------------*/

static void error(int flag)

{
    fflush(stdout);
    if(messages || flag) fprintf(stderr,"%s",err_msg[-flag]);
}

/*-------------------------------------------------Time_3d()----------------*/

int time_3d(float *HS, float *T, int NX, int NY, int NZ,
            float XS, float YS, float ZS, float HS_EPS_INIT, int MSG)

/* This function merely does nothing else than copying its arguments  */
/* to internal (static) variables. This allows you to alter the user  */
/* interface (e.g., pass on vector N[3] instead of NX,NY,NZ, or use   */
/* members of a struct or whatever else fits your needs) very easily. */
/* Arrays are passed as 1D vectors in order to make life easier with  */
/* Fortran calling programs...                                        */

{
    int    signal;

#ifdef DEBUG_ARGS
    fprintf(stderr,"Entering function time_3d using C-style interface\n");
    fprintf(stderr,"Dimensions nx=%d ny=%d nz=%d\n",NX,NY,NZ);
    fprintf(stderr,"(arrays stored in memory as C-array[nx][ny][nz])\n");
    fprintf(stderr,"Source location xs=%g ys=%g zs=%g\n",XS,YS,ZS);
    fprintf(stderr,"Tolerance used (HS_EPS_INIT) %g\n",HS_EPS_INIT);
    fprintf(stderr,"Verbosity (MSG) %d (%s)\n",MSG,(MSG)?"verbose":"silent");
    fflush(stderr);
#endif