time_3d.c

3D ray trace of elstic wave in homogenious medium

/*----------------------------------------------------------------------------*/
/*  TIME_3D: FINITE DIFFERENCE COMPUTATION OF FIRST ARRIVAL TIMES IN 3D.      */
/*----------------------------------------------------------------------------*/
/*  P.PODVIN, CGGM/Geophysique, Ecole des Mines de Paris, Fontainebleau.      */
/*  E-MAIL: [e-mail]            Tel: [phone].          */
/*  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 void
error(flag)

int flag;

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

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

int
time_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 ---------------*/
 
int
time_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 int
pre_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 **))))