time_3d.c

射线追踪程序

    /* limits are decremented so that interfaces where heterogeneities     */
    /* were detected are dealt with by finite differences (cf. headwaves). */
    /* (but this is not necessary when located at model boundaries !)      */

    if( init_stage>=INIT_RECURS_LIMIT ||
        (   (X0==0 || (xs-X0)>=INIT_MIN) &&
            (Y0==0 || (ys-Y0)>=INIT_MIN) &&
            (Z0==0 || (zs-Z0)>=INIT_MIN) &&
            (X1==nmesh_x || (X1-xs)>=INIT_MIN) &&
            (Y1==nmesh_y || (Y1-ys)>=INIT_MIN) &&
            (Z1==nmesh_z || (Z1-zs)>=INIT_MIN)     )   ) {
        if((X1-X0+1)*(Y1-Y0+1)*(Z1-Z0+1)==1)
            init_nearest();
        else for(x=X0;x<=X1;x++)
            for(y=Y0;y<=Y1;y++)
                for(z=Z0;z<=Z1;z++){
                    dist=(x-fxs)*(x-fxs)+(y-fys)*(y-fys)+(z-fzs)*(z-fzs);
                    t[x][y][z]=hs0*sqrt(dist);
                }
        if(SMALLTALK)
            printf("\nHomogeneous region: x[%d->%d] y[%d->%d] z[%d->%d]\n",
                    X0,X1,Y0,Y1,Z0,Z1);
    } /* if smallest distance from source to boundaries of the homogeneous */
      /* box exceeds conventional limit INIT_MIN, OR if no further recursi-*/
      /* vity is allowed, then exact arrivals are computed in this region. */

    else {
        if((signal=recursive_init())!=NO_ERROR) return signal;
        X0=max(xs-INIT_MIN,0);
        Y0=max(ys-INIT_MIN,0);
        Z0=max(zs-INIT_MIN,0);
        X1=min(xs+INIT_MIN,nmesh_x);
        Y1=min(ys+INIT_MIN,nmesh_y);
        Z1=min(zs+INIT_MIN,nmesh_z);
    } /* otherwise, time_3d() is used recursively   */
      /* on a re-discretized (2*INIT_MIN+1)^3 cube. */
 
    return signal;

}

/*------------------------------------------------Init_nearest()------------*/

static void init_nearest(void)

/* initialize the 8|12|18 nearest neighbour nodes of the source    */
/* according to source position (inside a mesh or at a boundary).  */
/* WARNING: errors are maximal when the source is located close to */
/* a grid-point. Best configurations are close to the centre of a  */
/* mesh face, or of a mesh. Errors increase (anisotropically) when */
/* the source gets close to a grid-point. Better use the grid-     */
/* point itself as the source in such case...                      */
{
    int   x,y,z;
    float distx,disty,distz;

    if(source_at_node){
        t[xs][ys][zs]=0.0;
        return;
    }
    x=(fxs<xs) ? xs-1:xs;
    y=(fys<ys) ? ys-1:ys;
    z=(fzs<zs) ? zs-1:zs;
    /* x,y,z : coordinates of current cell */
    distx=fabs(fxs-x);
    disty=fabs(fys-y);
    distz=fabs(fzs-z);
    /* dist* : distances from source to node minx,miny,minz of current cell */

    init_cell(distx,disty,distz,x,y,z);
    /* this is enough if the source is strictly located */
    /* within the current cell (init: 8 neighbours).    */

    if(fxs==xs){
        if(fys==ys){
            if(x) init_cell(1.,0.,distz,x-1,y,z);
            if(y) init_cell(0.,1.,distz,x,y-1,z);
            if(x && y) init_cell(1.,1.,distz,x-1,y-1,z);
        }/* source located on cell edge parallel to z (18 neighbours) */
        else
        if(fzs==zs){
            if(x) init_cell(1.,disty,0.,x-1,y,z);
            if(z) init_cell(0.,disty,1.,x,y,z-1);
            if(z && x) init_cell(1.,disty,1.,x-1,y,z-1);
        }/* source located on cell edge parallel to y (18 neighbours) */
        else{
            if(x) init_cell(1.,disty,distz,x-1,y,z);
        }/* source located on cell face perpendicular to x (12 neighbours) */
    }
    else
    if(fys==ys){
        if(fzs==zs){
            if(y) init_cell(distx,1.,0.,x,y-1,z);
            if(z) init_cell(distz,0.,1.,x,y,z-1);
            if(y && z) init_cell(distx,1.,1.,x,y-1,z-1);
        }/* source located on cell edge parallel to x (18 neighbours) */
        else {
            if(y) init_cell(distx,1.,distz,x,y-1,z);
        }/* source located on cell face perpendicular to y (12 neighbours) */
    }
    else
    if(fzs==zs){
        if(y) init_cell(distx,disty,1.,x,y,z-1);
    }/* source located on cell face perpendicular to z (12 neighbours) */

}
/*------------------------------------------------Init_cell()---------------*/

static void init_cell(float vx, float vy, float vz, int xl, int yl, int zl)

/* compute delays between floating source and nodes of current cell     */
/* xl,yl,zl are current cell coordinates,                               */
/* vx,vy,vz are distances from source to node xl,yl,zl (0<=vx<=1.0,...) */

{
    float est;
    est=exact_delay(vx,vy,vz,xl,yl,zl);
    if(est<t[xl][yl][zl]) t[xl][yl][zl]=est;
    est=exact_delay(1.0-vx,vy,vz,xl,yl,zl);
    if(est<t[xl+1][yl][zl]) t[xl+1][yl][zl]=est;
    est=exact_delay(vx,1.0-vy,vz,xl,yl,zl);
    if(est<t[xl][yl+1][zl]) t[xl][yl+1][zl]=est;
    est=exact_delay(vx,vy,1.0-vz,xl,yl,zl);
    if(est<t[xl][yl][zl+1]) t[xl][yl][zl+1]=est;
    est=exact_delay(1.0-vx,1.0-vy,vz,xl,yl,zl);
    if(est<t[xl+1][yl+1][zl]) t[xl+1][yl+1][zl]=est;
    est=exact_delay(1.0-vx,vy,1.0-vz,xl,yl,zl);
    if(est<t[xl+1][yl][zl+1]) t[xl+1][yl][zl+1]=est;
    est=exact_delay(vx,1.0-vy,1.0-vz,xl,yl,zl);
    if(est<t[xl][yl+1][zl+1]) t[xl][yl+1][zl+1]=est;
    est=exact_delay(1.0-vx,1.0-vy,1.0-vz,xl,yl,zl);
    if(est<t[xl+1][yl+1][zl+1]) t[xl+1][yl+1][zl+1]=est;
}

/*------------------------------------------------Recursive_init()----------*/

static int recursive_init(void)

{
    int
            signal,
            nx_,ny_,nz_,
            xs_,ys_,zs_,
            X0_,X1_,Y0_,Y1_,Z0_,Z1_,
            n,d,
            i,ii,ihs,i0,
            j,jj,jhs,j0,
            k,kk,khs,k0;
    float
            *hs_buf_,*t_buf_,
            fxs_,fys_,fzs_,
            HS[N_INIT],T[N_INIT];

/* increment count of recursivity level */
    init_stage++;
    if(SMALLTALK)
        printf("\nRecursive initialization: level %d",init_stage);

/* free locally allocated pointers (float ***) */
    free_ptrs(nx);

/* save static parameters at this stage */
    nx_=nx;
    ny_=ny;
    nz_=nz;
    hs_buf_=hs_buf;
    t_buf_=t_buf;
    xs_=xs;
    ys_=ys;
    zs_=zs;
    fxs_=fxs;
    fys_=fys;
    fzs_=fzs;
    X0_=X0;
    X1_=X1;
    Y0_=Y0;
    Y1_=Y1;
    Z0_=Z0;
    Z1_=Z1;

/* build the re-discretized local model and the associated source position */
    for(i=0;i<N_INIT;i++) HS[i]=T[i]=INFINITY;
    nx=ny=nz=N_INIT_X;
    xs=ys=zs=2*INIT_MIN+1;
    i0=j0=k0=1;
    ihs=xs_-INIT_MIN-1;
    if((d=INIT_MIN-xs_)>=0){
        ihs+=d+1;
        d=1+2*d;
        nx-=d;
        xs-=d;
        i0=0;
    }
    if((d=xs_+INIT_MIN-nx_+1)>=0) nx-=1+2*d;
    jhs=ys_-INIT_MIN-1;
    if((d=INIT_MIN-ys_)>=0){
        jhs+=d+1;
        d=1+2*d;
        ny-=d;
        ys-=d;
        j0=0;
    }
    if((d=ys_+INIT_MIN-ny_+1)>=0) ny-=1+2*d;
    khs=zs_-INIT_MIN-1;
    if((d=INIT_MIN-zs_)>=0){
        khs+=d+1;
        d=1+2*d;
        nz-=d;
        zs-=d;
        k0=0;
    }
    if((d=zs_+INIT_MIN-nz_+1)>=0) nz-=1+2*d;
    for(i=ihs,n=ii=0;ii<nx;ii++){
        for(j=jhs,jj=0;jj<ny;jj++){
            for(k=khs,kk=0;kk<nz;kk++,n++){
                HS[n]=0.5*hs_buf_[i*ny_*nz_+j*nz_+k];
                if(kk%2!=k0) k++;
            }
            if(jj%2!=j0) j++;
        }
        if(ii%2!=i0) i++;
    }/* No smoothing is associated with this re-discretization */
    fxs=xs+2.0*(fxs_-xs_);
    fys=ys+2.0*(fys_-ys_);
    fzs=zs+2.0*(fzs_-zs_);

    if(VERBOSE)
        printf("\nRediscretized timefield dimensions: %d %d %d",nx,ny,nz);

/* recursively compute times on this rediscretized model */
    signal=time_3d(HS,T,nx,ny,nz,fxs,fys,fzs,hs_eps_init,messages);

/* assign relevant times to parent timefield */
    if(signal==NO_ERROR){
        for(i=ihs+i0,ii=i0;ii<nx;ii+=2,i++)
            for(j=jhs+j0,jj=j0;jj<ny;jj+=2,j++)
                for(k=khs+k0,kk=k0;kk<nz;kk+=2,k++)
                    t_buf_[i*ny_*nz_+j*nz_+k]=T[ii*ny*nz+jj*nz+kk];
    }
 
/* retrieve initial static parameters */
    nx=nx_;
    ny=ny_;
    nz=nz_;
    hs_buf=hs_buf_;
    t_buf=t_buf_;
    xs=xs_;
    ys=ys_;
    zs=zs_;
    fxs=fxs_;
    fys=fys_;
    fzs=fzs_;
    X0=X0_;
    X1=X1_;
    Y0=Y0_;
    Y1=Y1_;
    Z0=Z0_;
    Z1=Z1_;
 
/* reallocate pointers (but do not re-initialize!) */
    signal=pre_init();

/* decrement count of recursivity level */
    init_stage--;
 
    return signal;
 
}

/*------------------------------------------------Propagate_point()---------*/

static int propagate_point(int start)

{
    int
        msg,test;

    if(start!=NO_ERROR) return start; /* Initialization failed */

    sum_updated=0;

/* Make recursive_init silent */
    if(SMALLTALK) printf("\nStarting F.D. computation...");
    msg=messages;
    if(init_stage) messages=0;

/* Increment boundaries of timed zone as long as necessary... */
/* (Outwards propagation is adopted as an initial guess).     */
    do {
        test=0;

        if(X0>0){
            X0--;
            if(VERBOSE) printf("\nx_side %d->%d: ",X0+1,X0);
            x_side(Y0,Y1,Z0,Z1,X0,-1);
            test++;
        }

        if(Y0>0){
            Y0--;
            if(VERBOSE) printf("\ny_side %d->%d: ",Y0+1,Y0);
            y_side(X0,X1,Z0,Z1,Y0,-1);
            test++;
        }

        if(Z0>0){
            Z0--;
            if(VERBOSE) printf("\nz_side %d->%d: ",Z0+1,Z0);
            z_side(X0,X1,Y0,Y1,Z0,-1);
            test++;
        }

        if(X1<nmesh_x){
            X1++;
            if(VERBOSE) printf("\nx_side %d->%d: ",X1-1,X1);
            x_side(Y0,Y1,Z0,Z1,X1,1);
            test++;
        }

        if(Y1<nmesh_y){
            Y1++;
            if(VERBOSE) printf("\ny_side %d->%d: ",Y1-1,Y1);
            y_side(X0,X1,Z0,Z1,Y1,1);
            test++;
        }

        if(Z1<nmesh_z){
            Z1++;
            if(VERBOSE) printf("\nz_side %d->%d: ",Z1-1,Z1);
            z_side(X0,X1,Y0,Y1,Z1,1);
            test++;
        }

    } while(test);

    messages=msg;

    return NO_ERROR;

}

/*---------------------------------------------- Free_ptrs()------------------*/

static void free_ptrs(int max_x)

{
    int x,y,z;
    float *pf;