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 intscan_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 intscan_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 intscan_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]