time_3d.c

射线追踪程序

                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(int y_start, int y_end, int z_start, int z_end,
          int x0, int x, int 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
        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(int y_begin, int y_start, int z_start, int z_end,
          int x0, int x, int 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){