/*----------------------------------------------------------------------------*//*  TIME_3D: FINITE DIFFERENCE COMPUTATION OF FIRST ARRIVAL TIMES IN 3D.      *//*----------------------------------------------------------------------------*//*  P.PODVIN, CGGM/Geophysique, Ecole des Mines de Paris, Fontainebleau.      *//*  E-MAIL: podvin@cggm.ensmp.fr            Tel: 33-(1) 64 69 49 25.          *//*  April 1991, last revision: 7 May 1993                                     *//*                                                                            *//*  This is a plain K&R C (not ANSI-C) code.                                  *//*                                                                            *//*  USAGE: (from a program written in C or in FORTRAN)                        *//*       (int)i_status=time_3d(HS,T,NX,NY,NZ,XS,YS,ZS,HS_EPS_INIT,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.            *//*                          e.g. 0.0<=XS<=(float)(NX-1)                       *//*                          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.        *//*      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) *//*                                                                            *//*  BUGS :                                                                    *//*       All known bugs are ascribed to bad argument list  or 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 !)         *//*                                                                            *//*  REFERENCE: Podvin & Lecomte, Geophys.J.Intern. 105, 271-284, 1991.        *//*----------------------------------------------------------------------------*/#include <stdio.h>#include <math.h>#define MM_SQRT2     1.41421356237309504880#define MM_SQRT3     1.732050807568877076#define INFINITY    1.0e+19#define ISINF(x)    ((x)>1.0e+18)#define NINT(x)     (int)floor((x)+0.5)#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))))/*-------------------------------------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(/* int,int,int,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 *//* 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    2#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=2,                       /* 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_MULT         (-1)#define ERR_INF          (-2)#define ERR_RECURS       (-3)#define ERR_MALLOC       (-4)#define ERR_HS_EPS       (-5)#define ERR_NONPHYSICAL  (-6)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"            };/*-------------------------------------------------Error()------------------*/static voiderror(flag)int flag;{    fflush(stdout);    if(messages || flag) fprintf(stderr,"%s",err_msg[-flag]);}/*-------------------------------------------------Time_3d()----------------*/inttime_3d(HS,T,NX,NY,NZ,XS,YS,ZS,HS_EPS_INIT,MSG)float *HS,*T,HS_EPS_INIT,XS,YS,ZS;int   NX,NY,NZ,MSG;{    int    signal;    hs_buf=HS;    t_buf=T;    nx=NX;    ny=NY;    nz=NZ;    fxs=XS;    fys=YS;    fzs=ZS;    hs_eps_init=HS_EPS_INIT;    if(hs_eps_init<0.0 || hs_eps_init>1.0) {        error(ERR_HS_EPS);        return(ERR_HS_EPS);    }    if(MSG<0){        no_init=1;        MSG= -MSG;    }    messages=MSG;    if((signal=pre_init())==NO_ERROR){        signal=propagate_point(init_point());        free_ptrs(nx);    }    if(init_stage==0 || signal!=NO_ERROR) error(signal);    return(signal);} /*---------------------------- Time_3d_(): FORTRAN INTERFACE ---------------*/ inttime_3d_(HS,T,NZ,NY,NX,ZS,YS,XS,HS_EPS_INIT,MSG) /* All FORTRAN arguments are pointers. Dimensions X and Z are *//* swapped in order to fit FORTRAN mapping conventions.       */float *HS,*T,*HS_EPS_INIT,*XS,*YS,*ZS;int   *NX,*NY,*NZ,*MSG; {    int     signal;     hs_buf=HS;    t_buf=T;    nx= *NX;    ny= *NY;    nz= *NZ;       fxs= *XS;    fys= *YS;    fzs= *ZS;    hs_eps_init= *HS_EPS_INIT;    if(hs_eps_init<0.0 || hs_eps_init>1.0) {        error(ERR_HS_EPS);        return(ERR_HS_EPS);    }    if(*MSG<0){        no_init=1;        messages= -(*MSG);    }    else messages= *MSG;     if((signal=pre_init())==NO_ERROR){        signal=propagate_point(init_point());        free_ptrs(nx);    }    if(init_stage==0 || signal!=NO_ERROR) error(signal);    return(signal);}/*------------------------------------------------Pre_init()----------------*/static intpre_init(){    int        x,y,z,        np,nt,        n0,n1;    float        *pf;    nmesh_x=nx-1;    nmesh_y=ny-1;    nmesh_z=nz-1;    np=ny*nz;    nt=nx*np;    n1=max(nx,ny);    n0=max(nx,nz);    if(n1==n0) n0=max(ny,nz);    n1*=n0;/* allocate pointers */    if(!(hs=(float ***)malloc((unsigned)nx*sizeof(float **))))