time_3d.c

3D ray trace of elstic wave in homogenious medium

        hs_bf,hs_bb,hs_fb,
        hs_ube,hs_ubb,hs_ueb;

    x_sf=x_s+x-x0;

/* We first propagate headwaves along the two borders of the current zone */
/* These headwaves are usually relevant only when a local minimum valley  */
/* is present along these borders on the preceding x_side (x=x0).         */
/* This is analogous with timing local minima in the 2-D implementation.  */

/* interface waves along y_side: 1 1D transmission and 1 2D transmission */
    hs_bb=hs_ubb=hs_ube=INFINITY;
    for(y=y_start,z=z_start;y<y_end;y++){
        hs_bf=hs[x_s][y][z];
        if(z) hs_bb=hs[x_s][y][z-1];
        if(x_sf>=0 && x_sf<nmesh_x){
            hs_ube=hs[x_sf][y][z];
            if(z) hs_ubb=hs[x_sf][y][z-1];
        }
        alert1=t_1d(x,y+1,z,t[x][y][z],hs_bb,hs_bf,hs_ubb,hs_ube);
        alert0=t_2d(x,y+1,z,t[x0][y][z],t[x][y][z],hs_bb,hs_bf);
        updated+=alert0+alert1;
        if(alert1) longflags[y*nz+nz+z]=1;
        if(alert0) longflags[y*nz+nz+z]=0;
    }

/* interface waves along z_side: 1 1D transmission and 1 2D transmission */
    hs_bb=hs_ubb=hs_ueb=INFINITY;
    for(y=y_start,z=z_start;z<z_end;z++){
        hs_fb=hs[x_s][y][z];
        if(y) hs_bb=hs[x_s][y-1][z];
        if(x_sf>=0 && x_sf<nmesh_x){
            hs_ueb=hs[x_sf][y][z];
            if(y) hs_ubb=hs[x_sf][y-1][z];
        }
        alert1=t_1d(x,y,z+1,t[x][y][z],hs_bb,hs_fb,hs_ubb,hs_ueb);
        alert0=t_2d(x,y,z+1,t[x0][y][z],t[x][y][z],hs_bb,hs_fb);
        updated+=alert0+alert1;
        if(alert1) longflags[y*nz+z+1]=1;
        if(alert0) longflags[y*nz+z+1]=0;
    }

/* We now propagate bulk and head waves into the region of interest. */

    for(y=y_start;y<y_end;y++){
        for(z=z_start;z<z_end;z++){

            hs_bb=hs[x_s][y][z];
            hs_bf=hs[x_s][y][z+1];
            hs_fb=hs[x_s][y+1][z];
            if(x_sf>=0 && x_sf<nmesh_x){
                hs_ubb=hs[x_sf][y][z];
                hs_ube=hs[x_sf][y][z+1];
                hs_ueb=hs[x_sf][y+1][z];
            }
            else hs_ubb=hs_ube=hs_ueb=INFINITY;

/* bulk waves: 1 3D edge diffraction and 2 (*4) 3D transmission */
            alert0=edge_diff(x,y+1,z+1,t[x0][y][z],t[x][y][z],hs_bb)
                  +t_3d(x,y+1,z+1,t[x0][y][z],
                    t[x0][y+1][z],t[x][y][z],t[x][y+1][z],hs_bb)
                  +t_3d(x,y+1,z+1,t[x0][y][z],
                    t[x0][y][z+1],t[x][y][z],t[x][y][z+1],hs_bb);
            if(alert0){
                updated+=alert0;
                longflags[y*nz+nz+z+1]=0;
            }

/* interface waves along y_side: 1 1D transmission and 1 2D transmission */
            alert1=t_1d(x,y+1,z+1,t[x][y][z+1],hs_bb,hs_bf,hs_ubb,hs_ube);
            alert0=t_2d(x,y+1,z+1,t[x0][y][z+1],t[x][y][z+1],hs_bb,hs_bf);
            if(alert0+alert1){
                updated+=alert0+alert1;
                flag_fb++; /* this scan must be (re-)examined */
                if(y_start_fb>y) y_start_fb=y;
                if(z_start_fb<z+1) z_start_fb=z+1;
                if(alert1) longflags[y*nz+nz+z+1]=1;
                else       longflags[y*nz+nz+z+1]=0;
            }

/* interface waves along z_side: 1 1D transmission and 1 2D transmission */
            alert1=t_1d(x,y+1,z+1,t[x][y+1][z],hs_bb,hs_fb,hs_ubb,hs_ueb);
            alert0=t_2d(x,y+1,z+1,t[x0][y+1][z],t[x][y+1][z],hs_bb,hs_fb);
            if(alert0+alert1){
                updated+=alert0+alert1;
                flag_bf++; /* this scan must be (re-)examined */
                if(y_start_bf<y+1) y_start_bf=y+1;
                if(z_start_bf>z) z_start_bf=z;
                if(alert1) longflags[y*nz+nz+z+1]=1;
                else       longflags[y*nz+nz+z+1]=0;
            }

/* interface waves along x_side : 2 2D transmission and 1 2D diffraction */
            alert1=diff_2d(x,y+1,z+1,t[x][y][z],hs_bb,hs_ubb)
                  +t_2d(x,y+1,z+1,t[x][y][z],t[x][y+1][z],hs_bb,hs_ubb)
                  +t_2d(x,y+1,z+1,t[x][y][z],t[x][y][z+1],hs_bb,hs_ubb);
            if(alert1){
                updated+=alert1;
                longflags[y*nz+nz+z+1]=1;
            }
        }
    }

    flag_ff=0;
    y_start_ff=y_end;
    z_start_ff=z_end;
    /* this direction has been examined: unset corresponding flag */

    return(updated);
}

/*--------------------------------------X_SIDE() : SCAN_X_BE()--------------*/

static int
scan_x_bf(y_begin,y_start,z_start,z_end,x0,x,x_s)
int    y_start,y_begin,z_start,z_end,x0,x,x_s;

{    int
        updated=0,
        alert0,alert1,
        y,z,
        x_sf;
    float
        hs_ff,hs_bb,hs_fb,
        hs_uee,hs_ubb,hs_ueb;

    x_sf=x_s+x-x0;

    hs_fb=hs_uee=hs_ueb=INFINITY;
    for(y=y_start,z=z_start;y>y_begin;y--){
        hs_ff=hs[x_s][y-1][z];
        if(z) hs_fb=hs[x_s][y-1][z-1];
        if(x_sf>=0 && x_sf<nmesh_x){
            hs_uee=hs[x_sf][y-1][z];
            if(z) hs_ueb=hs[x_sf][y-1][z-1];
        }
        alert1=t_1d(x,y-1,z,t[x][y][z],hs_fb,hs_ff,hs_ueb,hs_uee);
        alert0=t_2d(x,y-1,z,t[x0][y][z],t[x][y][z],hs_fb,hs_ff);
        updated+=alert0+alert1;
        if(alert1) longflags[y*nz-nz+z]=1;
        if(alert0) longflags[y*nz-nz+z]=0;
    }

    hs_bb=hs_ubb=hs_ueb=INFINITY;
    for(y=y_start,z=z_start;z<z_end;z++){
        if(y) hs_bb=hs[x_s][y-1][z];
        hs_fb=hs[x_s][y][z];
        if(x_sf>=0 && x_sf<nmesh_x){
            if(y) hs_ubb=hs[x_sf][y-1][z];
            hs_ueb=hs[x_sf][y][z];
        }
        alert1=t_1d(x,y,z+1,t[x][y][z],hs_bb,hs_fb,hs_ubb,hs_ueb);
        alert0=t_2d(x,y,z+1,t[x0][y][z],t[x][y][z],hs_bb,hs_fb);
        updated+=alert0+alert1;
        if(alert1) longflags[y*nz+z+1]=1;
        if(alert0) longflags[y*nz+z+1]=0;
    }

    for(y=y_start;y>y_begin;y--){
        for(z=z_start;z<z_end;z++){

            hs_ff=hs[x_s][y-1][z+1];
            if(y>1) hs_bb=hs[x_s][y-2][z];
            else hs_bb=INFINITY;
            hs_fb=hs[x_s][y-1][z];
            if(x_sf>=0 && x_sf<nmesh_x){
                hs_uee=hs[x_sf][y-1][z+1];
                if(y>1) hs_ubb=hs[x_sf][y-2][z];
                else hs_ubb=INFINITY;
                hs_ueb=hs[x_sf][y-1][z];
            }
            else hs_ubb=hs_uee=hs_ueb=INFINITY;

/* bulk waves: 1 3D edge diffraction and 2 (*4) 3D transmission */
            alert0=edge_diff(x,y-1,z+1,t[x0][y][z],t[x][y][z],hs_fb)
                  +t_3d(x,y-1,z+1,t[x0][y][z],
                    t[x0][y-1][z],t[x][y][z],t[x][y-1][z],hs_fb)
                  +t_3d(x,y-1,z+1,t[x0][y][z],
                    t[x0][y][z+1],t[x][y][z],t[x][y][z+1],hs_fb);
            if(alert0){
                updated+=alert0;
                longflags[y*nz-nz+z+1]=0;
            }

/* interface waves along y_side: 1 1D transmission and 1 2D transmission */
            alert1=t_1d(x,y-1,z+1,t[x][y][z+1],hs_ff,hs_fb,hs_uee,hs_ueb);
            alert0=t_2d(x,y-1,z+1,t[x0][y][z+1],t[x][y][z+1],hs_ff,hs_fb);
            if(alert0+alert1){
                updated+alert0+alert1;
                flag_bb++; /* this scan must be (re-)examined */
                if(y_start_bb<y) y_start_bb=y;
                if(z_start_bb<z+1) z_start_bb=z+1;
                if(alert1) longflags[y*nz-nz+z+1]=1;
                else       longflags[y*nz-nz+z+1]=0;
            }

/* interface waves along z_side: 1 1D transmission and 1 2D transmission */
            alert1=t_1d(x,y-1,z+1,t[x][y-1][z],hs_bb,hs_fb,hs_ubb,hs_ueb);
            alert0=t_2d(x,y-1,z+1,t[x0][y-1][z],t[x][y-1][z],hs_bb,hs_fb);
            if(alert0+alert1){
                updated+alert0+alert1;
                flag_ff++; /* this scan must be (re-)examined */
                if(y_start_ff>y-1) y_start_ff=y-1;
                if(z_start_ff>z) z_start_ff=z;
                if(alert1) longflags[y*nz-nz+z+1]=1;
                else       longflags[y*nz-nz+z+1]=0;
            }

/* interface waves along x_side : 2 2D transmission and 1 2D diffraction */
            alert1=diff_2d(x,y-1,z+1,t[x][y][z],hs_fb,hs_ueb)
                  +t_2d(x,y-1,z+1,t[x][y][z],t[x][y-1][z],hs_fb,hs_ueb)
                  +t_2d(x,y-1,z+1,t[x][y][z],t[x][y][z+1],hs_fb,hs_ueb);
            if(alert1){
                updated+=alert1;
                longflags[y*nz-nz+z+1]=1;
            }
        }
    }

    flag_bf=0;
    y_start_bf=y_begin;
    z_start_bf=z_end;
    /* this direction has been examined: unset corresponding flag */

    return(updated);
}

/*--------------------------------------X_SIDE() : SCAN_X_BB()--------------*/

static int
scan_x_bb(y_begin,y_start,z_begin,z_start,x0,x,x_s)
int    y_start,y_begin,z_begin,z_start,x0,x,x_s;

{    int
        updated=0,
        alert0,alert1,
        y,z,
        x_sf;
    float
        hs_ff,hs_bf,hs_fb,
        hs_uee,hs_ube,hs_ueb;

    x_sf=x_s+x-x0;

    hs_ff=hs_uee=hs_ueb=INFINITY;
    for(y=y_start,z=z_start;y>y_begin;y--){
        if(z) hs_ff=hs[x_s][y-1][z-1];
        hs_fb=hs[x_s][y-1][z];
        if(x_sf>=0 && x_sf<nmesh_x){
            if(z) hs_uee=hs[x_sf][y-1][z-1];
            hs_ueb=hs[x_sf][y-1][z];
        }
        alert1=t_1d(x,y-1,z,t[x][y][z],hs_fb,hs_ff,hs_ueb,hs_uee);
        alert0=t_2d(x,y-1,z,t[x0][y][z],t[x][y][z],hs_fb,hs_ff);
        updated+=alert0+alert1;
        if(alert1) longflags[y*nz-nz+z]=1;
        if(alert0) longflags[y*nz-nz+z]=0;
    }

    hs_ff=hs_uee=hs_ube=INFINITY;
    for(y=y_start,z=z_start;z>z_begin;z--){
        if(y) hs_ff=hs[x_s][y-1][z-1];
        hs_bf=hs[x_s][y][z-1];
        if(x_sf>=0 && x_sf<nmesh_x){
            if(y) hs_uee=hs[x_sf][y-1][z-1];
            hs_ube=hs[x_sf][y][z-1];
        }
        alert1=t_1d(x,y,z-1,t[x][y][z],hs_bf,hs_ff,hs_ube,hs_uee);
        alert0=t_2d(x,y,z-1,t[x0][y][z],t[x][y][z],hs_bf,hs_ff);
        updated+=alert0+alert1;
        if(alert1) longflags[y*nz+z-1]=1;
        if(alert0) longflags[y*nz+z-1]=0;
    }

    for(y=y_start;y>y_begin;y--){
        for(z=z_start;z>z_begin;z--){

            hs_ff=hs[x_s][y-1][z-1];
            if(y>1) hs_bf=hs[x_s][y-2][z-1];
            else hs_bf=INFINITY;
            if(z>1) hs_fb=hs[x_s][y-1][z-2];
            else hs_fb=INFINITY;
            if(x_sf>=0 && x_sf<nmesh_x){
                hs_uee=hs[x_sf][y-1][z-1];
                if(y>1) hs_ube=hs[x_sf][y-2][z-1];
                else hs_ube=INFINITY;
                if(z>1) hs_ueb=hs[x_sf][y-1][z-2];
                else hs_ueb=INFINITY;
            }
            else hs_uee=hs_ube=hs_ueb=INFINITY;

/* bulk waves: 1 3D edge diffraction and 2 (*4) 3D transmission */
            alert0=edge_diff(x,y-1,z-1,t[x0][y][z],t[x][y][z],hs_ff)
                  +t_3d(x,y-1,z-1,t[x0][y][z],
                        t[x0][y-1][z],t[x][y][z],t[x][y-1][z],hs_ff)
                  +t_3d(x,y-1,z-1,t[x0][y][z],
                        t[x0][y][z-1],t[x][y][z],t[x][y][z-1],hs_ff);
            if(alert0){
                updated+=alert0;
                longflags[y*nz-nz+z-1]=0;
            }

/* interface waves along y_side: 1 1D transmission and 1 2D transmission */
            alert1=t_1d(x,y-1,z-1,t[x][y][z-1],hs_ff,hs_fb,hs_uee,hs_ueb);
            alert0=t_2d(x,y-1,z-1,t[x0][y][z-1],t[x][y][z-1],hs_ff,hs_fb);
            if(alert0+alert1){
                updated+=alert0+alert1;
                flag_bf++; /* this scan must be (re-)examined */
                if(y_start_bf<y) y_start_bf=y;
                if(z_start_bf>z-1) z_start_bf=z-1;
                if(alert1) longflags[y*nz-nz+z-1]=1;
                else       longflags[y*nz-nz+z-1]=0;
            }

/* interface waves along z_side: 1 1D transmission and 1 2D transmission */
            alert1=t_1d(x,y-1,z-1,t[x][y-1][z],hs_ff,hs_bf,hs_uee,hs_ube);
            alert0=t_2d(x,y-1,z-1,t[x0][y-1][z],t[x][y-1][z],hs_ff,hs_bf);
            if(alert0+alert1){
                updated+=alert0+alert1;
                flag_fb++; /* this scan must be (re-)examined */
                if(y_start_fb>y-1) y_start_fb=y-1;
                if(z_start_fb<z) z_start_fb=z;
                if(alert1) longflags[y*nz-nz+z-1]=1;
                else       longflags[y*nz-nz+z-1]=0;
            }

/* interface waves along x_side : 2 2D transmission and 1 2D diffraction */
            alert1=diff_2d(x,y-1,z-1,t[x][y][z],hs_ff,hs_uee)
                  +t_2d(x,y-1,z-1,t[x][y][z],t[x][y-1][z],hs_ff,hs_uee)
                  +t_2d(x,y-1,z-1,t[x][y][z],t[x][y][z-1],hs_ff,hs_uee);
            if(alert1){
                updated+=alert1;
                longflags[y*nz-nz+z-1]=1;
            }
        }
    }

    flag_bb=0;
    y_start_bb=y_begin;
    z_start_bb=z_begin;
    /* this direction has been examined: unset corresponding flag */

    return(updated);
}

/*--------------------------------------X_SIDE() : SCAN_X_EB()--------------*/

static int
scan_x_fb(y_start,y_end,z_begin,z_start,x0,x,x_s)
int    y_start,y_end,z_begin,z_start,x0,x,x_s;

{    int
        updated=0,
        alert0,alert1,
        y,z,
        x_sf;
    float
        hs_ff,hs_bb,hs_bf,
        hs_uee,hs_ubb,hs_ube;

    x_sf=x_s+x-x0;

    hs_bf=hs_ubb=hs_ube=INFINITY;
    for(y=y_start,z=z_start;y<y_end;y++){
        if(z) hs_bf=hs[x_s][y][z-1];
        hs_bb=hs[x_s][y][z];
        if(x_sf>=0 && x_sf<nmesh_x){
            if(z) hs_ube=hs[x_sf]