time_3d.c

射线追踪程序

/* if relevant, retrieve INFINITY-masked hs values at model boundaries */
    if(init_stage==0 && hs_keep){
        pf=hs_keep;
        for(x=0;x<nx;x++){
            for(y=0;y<ny;y++) hs[x][y][nmesh_z]= *pf++;
            for(z=0;z<nmesh_z;z++) hs[x][nmesh_y][z]= *pf++;
        }
        for(y=0;y<nmesh_y;y++)
            for(z=0;z<nmesh_z;z++) hs[nmesh_x][y][z]= *pf++;
        free((char *)hs_keep);
    }

/* if relevant, undo global time-shift before returning */
    if(init_stage==0 && timeshift<0.0){
        for(x=0;x<nx;x++)
            for(y=0;y<ny;y++)
                for(z=0;z<nz;z++) t[x][y][z]+=timeshift;
    }

/* free pointers */
    for(x=0;x<max_x;x++){
        free((char *)hs[x]);
        free((char *)t[x]);
    }
    free((char *)hs);
    free((char *)t);
    free((char *)longflags);

}
/****end mail1/3****/
/*--------------------LOCAL 3-D STENCILS (FUNCTIONS AND MACROS)---------------*/
/* Counts of stencils refer to the local inwards isotropic formulation of the */
/* local finite difference computation function (170 stencils for a given     */
/* point, taken as the center of a 2*2*2 cube, with 26 timed neighbours).     */
/* See Podvin and Lecomte, 1991.                                              */
/* In this implementation, this function is never fully computed, because we  */
/* sequentially select only a limited number of relevant directions of propa- */
/* gation, according to the recent history of the signal.                     */
/* As a consequence, only 28 stencils are generally tested at each point.     */
/* The corresponding count is indicated between brackets.                     */
/*----------------------------------------------------------------------------*/
/* Code was restructured in order to minimize the number of actually computed */
/* sqrt()s. Some redundancy was introduced in tests in order to cope properly */
/* with problems arising from numerical errors (e.g. sqrt(x*x) may be lower   */
/* than x, a fact leading to infinite recursions in some awkward cases).      */
/*----------------------------------------------------------------------------*/

/*----------------------------------------- exact_delay() ------------------- */

static float exact_delay(float vx, float vy, float vz, int xm, int ym, int zm)

{
    float estimate;

    if(xm<0 || xm>=nmesh_x || ym<0 || ym>=nmesh_y || zm<0 || zm>=nmesh_z)
        return INFINITY;
    estimate=(vx*vx+vy*vy+vz*vz)*hs[xm][ym][zm]*hs[xm][ym][zm];
    return sqrt(estimate);
}

/*----------------------------------------- 1-D transmission : 6 stencils [3] */
/*------------------------------------- (Direct arrival from first neighbour) */

static int t_1d(int x, int y, int z,
                float t0, float hs0, float hs1, float hs2, float hs3)

{
    float estimate,dt;
    estimate=t0+min4(hs0,hs1,hs2,hs3);
    /* 12 July 2003: FUZZIFIED COMPARISON */
    dt=t[x][y][z]-estimate;
    if(dt>EPS_FUZZY*t[x][y][z]){
        t[x][y][z]=estimate;
        return 1;
    }
    return 0;
}

/*----------------------------------------- 2-D diffraction : 12 stencils [3] */
/*------------------------------------ (Direct arrival from second neighbour) */

static int diff_2d(int x, int y, int z, float t0, float hs0, float hs1)

{
    float estimate,dt;
    estimate=t0+M_SQRT2*min(hs0,hs1);
    /* 12 July 2003: FUZZIFIED COMPARISON */
    dt=t[x][y][z]-estimate;
    if(dt>EPS_FUZZY*t[x][y][z]){
        t[x][y][z]=estimate;
        return 1;
    }
    return 0;
}

/*------------------------------------ 3-D point diffraction : 8 stencils [1] */
/*------------------------------------- (Direct arrival from third neighbour) */

static int point_diff(int x, int y, int z, float t0, float hs0)

{
    float estimate;
    estimate=t0+hs0*M_SQRT3;
    if(estimate<t[x][y][z]){
        t[x][y][z]=estimate;
        return 1;
    }
    return 0;
}

/*---------------------------------------- 2-D transmission : 24 stencils [6] */
/*----------------------------------------- (Arrival from coplanar mesh edge) */

static int t_2d(int x, int y, int z, float t0, float t1, float hs0, float hs1)

{   float estimate,dt,hsm,test2,u2;
    dt=t1-t0;
    test2=t[x][y][z]-t1;
    if(dt<0.0 || test2<0.0) return 0;
    test2*=test2;
    hsm=min(hs0,hs1);
    u2=hsm*hsm-dt*dt;
    if(dt<=hsm/M_SQRT2 && u2<=test2){
        estimate=t1+sqrt(u2);
        /* 12 July 2003: FUZZIFIED COMPARISON */
        dt=t[x][y][z]-estimate;
        if(dt>EPS_FUZZY*t[x][y][z]){
            t[x][y][z]=estimate;
            return 1;
        }
    }
    return 0;
}

/*------------------------------------ 3-D edge diffraction : 24 stencils [3] */
/*------------------------------------- (Arrival from non-coplanar mesh edge) */

static int edge_diff(int x, int y, int z, float t0, float t1, float hs0)

{   float estimate,u2,test2,dt;
    dt=t1-t0;
    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/M_SQRT3 && 2.0*u2<=test2){
        estimate=t1+M_SQRT2*sqrt(u2);
        /* 12 July 2003: FUZZIFIED COMPARISON (not required...) */
        dt=t[x][y][z]-estimate;
        if(dt>EPS_FUZZY*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_(int x, int y, int z, float t0, float tl, float tr, float td,
                 float hs0, int 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(int x, int y, int z,
                      float t0, float tl, float tr, float 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(int y_begin, int y_end, int z_begin, int z_end,
                  int x, int 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;