time_3d.c

3D ray trace of elstic wave in homogenious medium

    test2=t[x][y][z]-t1;
    if(dt<0.0 || test2<0.0) return(0);
    test2*=test2;
    u2=hs0*hs0-dt*dt;
    if(dt<=hs0/MM_SQRT3 && 2.0*u2<=test2){
        estimate=t1+MM_SQRT2*sqrt(u2);
        if(estimate<t[x][y][z]){
            t[x][y][z]=estimate;
            return(1);
        }
    }
    return(0);
}

/*--------------------------------------- 3-D transmission : 96 stencils [12] */
/*------------------------------------- (Arrival from non-coplanar interface) */
/* 4 stencils per function call or 1+3 using two function calls. */

#define t_3d(x,y,z,a,b,c,d,e)        t_3d_(x,y,z,a,b,c,d,e,0)
#define t_3d_part2(x,y,z,a,b,c,d,e)  t_3d_(x,y,z,a,b,c,d,e,1)

static int
t_3d_(x,y,z,t0,tl,tr,td,hs0,redundant)
float t0,tl,tr,td,hs0;
int x,y,z,redundant;
/* The current point is in diagonal position with respect to t0     */
/* and it is a first neighbour of td. tl,tr are second neighbours.  */
/* One of these estimators is redundant during first step of *_side */
/* functions. See t_3d_part1() which also computes it.              */
/* This function is always called through macros t_3d or t_3d_part2 */
{   float test2,r2,s2,t2,u2,dta,dtb,dta2,dtb2,estimate;
    int action;
    action=0;
    hs0*=hs0;

    dta=tl-t0;
    dtb=tr-t0;
    dta2=dta*dta;
    dtb2=dtb*dtb;
    if(dta>=0.0 && dtb>=0.0 && dta2+dtb2+dta*dtb>=0.5*hs0
        && 2.0*dta2+dtb2<=hs0 && 2.0*dtb2+dta2<=hs0){
        test2=t[x][y][z]-tr-tl+t0;
        if(test2>=0.0){
            test2*=test2;
            r2=hs0-dta2-dtb2;
            if(r2<test2){
                estimate=tr+tl-t0+sqrt(r2);
                if(estimate<t[x][y][z]){
                    t[x][y][z]=estimate;
                    action++;
                }
            }
        }
    }

    test2=t[x][y][z]-td;
    if(test2<0.0) return(action);
    test2*=test2;
    s2=t2=u2=INFINITY;

    dtb=td-tl;
    dtb2=dtb*dtb;
    if(dta>=0.0 && dtb>=dta && 2.0*dtb2+dta2<=hs0)
        s2=hs0-dta2-dtb2;

    dta=td-tr;
    dta2=dta*dta;
    if(!redundant && dta>=0.0 && dtb>=0.0 && dta2+dtb2+dta*dtb<=0.5*hs0)
        t2=hs0-dta2-dtb2;

    dtb=tr-t0;
    dtb2=dtb*dtb;
    if(dtb>=0.0 && dta>=dtb && 2.0*dta2+dtb2<=hs0)
        u2=hs0-dta2-dtb2;

    u2=min3(s2,t2,u2);
    if(u2<test2){
        estimate=td+sqrt(u2);
        if(estimate<t[x][y][z]){
            t[x][y][z]=estimate;
            action++;
        }
    }
    return(action);
}

/* 3-D transmission: partial stencil, introduced because initial scheme */
/* failed to fullfill the exhaustivity condition requested by Fermat's  */
/* principle. (See "a-causal" step in *_side() functions; 18/07/91)     */
static int
t_3d_part1(x,y,z,t0,tl,tr,hs0)
int x,y,z;
float t0,tl,tr,hs0;
/* The current point is a first neighbour of t0; tl,tr are two other */
/* first neighbours of t0. Transmission through 0-l-r is tested.     */
{
    float dtl,dtr,s2,u2,estimate,test2;
    dtl=t0-tl;
    dtr=t0-tr;
    test2=t[x][y][z]-t0;
    if(test2<0.0 || dtl<0.0 || dtr<0.0) return(0);
    test2*=test2;
    hs0*=hs0;
    s2=dtl*dtl+dtr*dtr;
    if(s2+dtl*dtr>0.5*hs0) return(0);
    /* illumination condition */
    u2=hs0-s2;
    if(u2<test2){
        estimate=t0+sqrt(u2);
        if(estimate<t[x][y][z]){
            t[x][y][z]=estimate;
            return(1);
        }
    }
    return(0);
}

/*----------------------------------------------X_SIDE()--------------------*/

static int
x_side(y_begin,y_end,z_begin,z_end,x,future)

int y_begin,y_end,z_begin,z_end,x,future;

/* Propagates computations from side x-future to current side x */
/* between *_begin and *_end coordinates. Returns a nonzero     */
/* integer if something actually happened (a time was lowered). */
/* Extensions _bb, _fb etc... define simple orientation rules:  */
/* _bf means backwards along y axis and forwards along z axis.  */
/* So-called "longitudinal" headwaves refer to first arrivals   */
/* due to a headwave propagating along the current side.        */

{
    int
        updated,                           /* counts adopted FD stencils */
        longhead,                          /* counts "longitudinal" headwaves */
        x0,                                /* past side coordinate */
        x_s,                               /* current meshes coordinate */
        y,z,                               /* current point coordinate */
        sign_ff,sign_bf,sign_bb,sign_fb,   /* sign flags for time differences */
        past,                              /* opposite to future ! */
        test;
    float
        hs_ff,hs_bf,hs_bb,hs_fb;           /* local slownesses */

    if(reverse_order==0)
        current_side_limit=x+future;
    updated=0;
    x0=x-future;
    if(future==1) x_s=x0;
    else          x_s=x;

    flag_fb=flag_bf=flag_ff=flag_bb=0;
    y_start_ff=y_start_fb=y_end;
    y_start_bf=y_start_bb=y_begin;
    z_start_ff=z_start_bf=z_end;
    z_start_fb=z_start_bb=z_begin;

/* First,  Compute all stencils using only nodes of side x0.   */
/* As only times on side x will be changed, these stencils     */
/* are computed initially, once for all, in any order (no      */
/* causality problem involved).                                */
/* During this first pass, future directions of propagation    */
/* are diagnosed, according to the time pattern on side x0.    */
/* Borders of zones to be scanned are computed.                */
/* This part may be seen as the explicit part of the FD scheme.*/

    for(y=y_begin;y<=y_end;y++){
        for(z=z_begin;z<=z_end;z++){

            hs_ff=hs[x_s][y][z];
            if(y>0) hs_bf=hs[x_s][y-1][z];
            else hs_bf=INFINITY;
            if(z>0 && y>0) hs_bb=hs[x_s][y-1][z-1];
            else hs_bb=INFINITY;
            if(z>0) hs_fb=hs[x_s][y][z-1];
            else hs_fb=INFINITY;
            sign_fb=sign_bf=sign_ff=sign_bb=0;

/* illuminate first neighbours */
/* 1 1D transmission and 4 partial 3D transmission */
            updated+=t_1d(x,y,z,t[x0][y][z],hs_ff,hs_bf,hs_bb,hs_fb);
            if(y<y_end && z<z_end)
                updated+=t_3d_part1(x,y,z,
                            t[x0][y][z],t[x0][y+1][z],t[x0][y][z+1],hs_ff);
            if(y>y_begin && z<z_end)
                updated+=t_3d_part1(x,y,z,
                            t[x0][y][z],t[x0][y-1][z],t[x0][y][z+1],hs_bf);
            if(y>y_begin && z>z_begin)
                updated+=t_3d_part1(x,y,z,
                            t[x0][y][z],t[x0][y-1][z],t[x0][y][z-1],hs_bb);
            if(y<y_end && z>z_begin)
                updated+=t_3d_part1(x,y,z,
                            t[x0][y][z],t[x0][y+1][z],t[x0][y][z-1],hs_fb);

/* illuminate second neighbours (if necessary)    */
/* 4 2D diffraction and 4 2D transmission */
            if(y<y_end && t[x0][y][z]<=t[x0][y+1][z]){
                sign_fb++;
                sign_ff++;
                if(y<y_start_ff) y_start_ff=y;
                if(y<y_start_fb) y_start_fb=y;
                updated+=diff_2d(x,y+1,z,t[x0][y][z],hs_ff,hs_fb);
                updated+=t_2d(x,y+1,z,t[x0][y][z],t[x0][y+1][z],hs_ff,hs_fb);
            }
            if(y>y_begin && t[x0][y][z]<=t[x0][y-1][z]){
                sign_bb++; 
                sign_bf++;
                if(y>y_start_bf) y_start_bf=y;
                if(y>y_start_bb) y_start_bb=y;
                updated+=diff_2d(x,y-1,z,t[x0][y][z],hs_bf,hs_bb);
                updated+=t_2d(x,y-1,z,t[x0][y][z],t[x0][y-1][z],hs_bf,hs_bb);
            }
            if(z<z_end && t[x0][y][z]<=t[x0][y][z+1]){
                sign_bf++; 
                sign_ff++;
                if(z<z_start_ff) z_start_ff=z;
                if(z<z_start_bf) z_start_bf=z;
                updated+=diff_2d(x,y,z+1,t[x0][y][z],hs_ff,hs_bf);
                updated+=t_2d(x,y,z+1,t[x0][y][z],t[x0][y][z+1],hs_ff,hs_bf);
            }
            if(z>z_begin && t[x0][y][z]<=t[x0][y][z-1]){
                sign_bb++; 
                sign_fb++;
                if(z>z_start_fb) z_start_fb=z;
                if(z>z_start_bb) z_start_bb=z;
                updated+=diff_2d(x,y,z-1,t[x0][y][z],hs_bb,hs_fb);
                updated+=t_2d(x,y,z-1,t[x0][y][z],t[x0][y][z-1],hs_bb,hs_fb);
            }

/* illuminate third neighbours (if necessary) */
/* 4 3D point diffraction, 8 3D edge diffraction and 12 3D transmission */
            if(sign_ff==2){
                flag_ff=1;
                updated+=point_diff(x,y+1,z+1,t[x0][y][z],hs_ff);
                updated+=edge_diff(x,y+1,z+1,t[x0][y][z],t[x0][y+1][z],hs_ff);
                updated+=edge_diff(x,y+1,z+1,t[x0][y][z],t[x0][y][z+1],hs_ff);
                updated+=t_3d_part2(x,y+1,z+1,t[x0][y][z],
                    t[x0][y+1][z],t[x0][y][z+1],t[x0][y+1][z+1],hs_ff);
            }
            if(sign_bf==2){
                flag_bf=1;
                updated+=point_diff(x,y-1,z+1,t[x0][y][z],hs_bf);
                updated+=edge_diff(x,y-1,z+1,t[x0][y][z],t[x0][y-1][z],hs_bf);
                updated+=edge_diff(x,y-1,z+1,t[x0][y][z],t[x0][y][z+1],hs_bf);
                updated+=t_3d_part2(x,y-1,z+1,t[x0][y][z],
                    t[x0][y-1][z],t[x0][y][z+1],t[x0][y-1][z+1],hs_bf);
            }
            if(sign_bb==2){
                flag_bb=1;
                updated+=point_diff(x,y-1,z-1,t[x0][y][z],hs_bb);
                updated+=edge_diff(x,y-1,z-1,t[x0][y][z],t[x0][y-1][z],hs_bb);
                updated+=edge_diff(x,y-1,z-1,t[x0][y][z],t[x0][y][z-1],hs_bb);
                updated+=t_3d_part2(x,y-1,z-1,t[x0][y][z],
                    t[x0][y-1][z],t[x0][y][z-1],t[x0][y-1][z-1],hs_bb);
            }
            if(sign_fb==2){
                flag_fb=1;
                updated+=point_diff(x,y+1,z-1,t[x0][y][z],hs_fb);
                updated+=edge_diff(x,y+1,z-1,t[x0][y][z],t[x0][y+1][z],hs_fb);
                updated+=edge_diff(x,y+1,z-1,t[x0][y][z],t[x0][y][z-1],hs_fb);
                updated+=t_3d_part2(x,y+1,z-1,t[x0][y][z],
                    t[x0][y+1][z],t[x0][y][z-1],t[x0][y+1][z-1],hs_fb);
            }
        }
    }

/* Now, all remaining stencils depend on nodes located on the current  */
/* side. They must be propagated causally, and occurrences of critical */
/* conditions must be diagnosed, in order to propagate associated head */
/* waves exhaustively. Four independent scanning directions are succes-*/
/* sively explored, and headwave flags are used to generate the supple-*/
/* mentary scans requested by exhaustivity ("Reverse" Propagations).   */
/* flag_* flag  is non-zero while the corresponding direction remains  */
/* to be examined (or reexamined). Its examination may detect critical */
/* conditions which in their turn will make another scan necessary.    */
/* Initialization of this process was achieved during first step.      */
/* This second step may be seen as the implicit part of the FD scheme. */

/* initialize local headwave flags */
    for(y=0;y<ny*nz;y++) longflags[y]=0;

/* enforce all scans if current side has a null surface */
/* (This may only be encountered near the source point) */
    if(y_begin==y_end || z_begin==z_end){
        flag_ff=flag_fb=flag_bf=flag_bb=1;
        y_start_ff=y_start_fb=y_begin;
        y_start_bf=y_start_bb=y_end;
        z_start_ff=z_start_bf=z_begin;
        z_start_fb=z_start_bb=z_end;
    }

/* Reexamine each direction, while necessary */
    do{
        test=0;
        if(flag_ff){
            test++;
            if(VERBOSE) printf("ff ");
            updated+=scan_x_ff(y_start_ff,y_end,z_start_ff,z_end,x0,x,x_s);
        }
        if(flag_fb){
            test++;
            if(VERBOSE) printf("fb ");
            updated+=scan_x_fb(y_start_fb,y_end,z_begin,z_start_fb,x0,x,x_s);
        }
        if(flag_bb){
            test++;
            if(VERBOSE) printf("bb ");
            updated+=scan_x_bb(y_begin,y_start_bb,z_begin,z_start_bb,x0,x,x_s);
        }
        if(flag_bf){
            test++;
            if(VERBOSE) printf("bf ");
            updated+=scan_x_bf(y_begin,y_start_bf,z_start_bf,z_end,x0,x,x_s);
        }
    } while(test);

    sum_updated+=updated;

/* At this stage, all points of the current side have been timed.     */
/* Now, Reverse propagation must be invoked if a headwave propagating */
/* along the current side was generated, because a new branch of the  */
/* timefield (conical wave) propagates towards the already timed box. */

    for(y=longhead=0;y<ny*nz;y++) longhead+=longflags[y];

    if(longhead){

        reverse_order++;

        if(VERBOSE) printf("\nReverse#%d from x_side %d",reverse_order,x);
        past= -future;
        for(x=x0;x!=current_side_limit;x+=past){
            if(x<0 || x>=nx) break;
            if(VERBOSE) printf("\nupdate side x=%d: ",x);
            if(x_side(y_begin,y_end,z_begin,z_end,x,past)==0) break;
            if(VERBOSE) printf("x=%d <R#%d>updated.",x,reverse_order);
        }
        if(VERBOSE) printf("\nEnd Reverse#%d\n",reverse_order);

        reverse_order--;

    }

    return(updated);

}

/*--------------------------------------X_SIDE() : SCAN_X_EE()--------------*/

static int
scan_x_ff(y_start,y_end,z_start,z_end,x0,x,x_s)
int    y_start,y_end,z_start,z_end,x0,x,x_s;

/* scan x_side by increasing y and z ("ff"=forwards, forwards)      */
/* propagating causal stencils with provisional a-priori that some  */
/* significant wavefronts propagate in this direction in the zone   */
/* defined by y_start,y_end,z_start,z_end.                          */
/* Critical conditions on any interface are detected so that other  */
/* relevant directions of propagation due to headwave generation    */
/* may be exhaustively taken into account at a later stage.         */

{    int
        updated=0,
        alert0,alert1,
        y,z,
        x_sf;
    float