time_3d.c

3D ray trace of elstic wave in homogenious medium

    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()
{
    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];
/*    FILE*   *tmpfile;   */

/* 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 */
/*sh    tmpfile = fopen("vel-src.tmp", "w");  sh*/
    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;
/*sh    fprintf(tmpfile,"\nsrc pos. xs_: %d  ys_: %d  zs_:  %d",xs_,ys_,zs_); */
    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;
    }
/*sh    fprintf(tmpfile,"\nnew dim. nx: %d ny: %d  nz: %d",nx,ny,nz);
    fprintf(tmpfile,"\n ihs: %d  jhs: %d  khs: %d",ihs,jhs,khs);  */
    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];
/*		fprintf(tmpfile,"\n %d %d %d %f",i,j,k,0.5/HS[n]); */
                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(start)

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(max_x)

int max_x;

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

/* 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);
    }

/* 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);

}

/*--------------------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(vx,vy,vz,xm,ym,zm)

float vx,vy,vz;
int xm,ym,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(x,y,z,t0,hs0,hs1,hs2,hs3)
int x,y,z;
float t0,hs0,hs1,hs2,hs3;
{
    float estimate;
    estimate=t0+min4(hs0,hs1,hs2,hs3);
    if(estimate<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(x,y,z,t0,hs0,hs1)
int x,y,z;
float t0,hs0,hs1;
{
    float estimate;
    estimate=t0+MM_SQRT2*min(hs0,hs1);
    if(estimate<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(x,y,z,t0,hs0)
int x,y,z;
float t0,hs0;
{
    float estimate;
    estimate=t0+hs0*MM_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(x,y,z,t0,t1,hs0,hs1)
int x,y,z;
float t0,t1,hs0,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/MM_SQRT2 && u2<=test2){
        estimate=t1+sqrt(u2);
        if(estimate<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(x,y,z,t0,t1,hs0)
int x,y,z;
float t0,t1,hs0;

{   float estimate,u2,test2,dt;
    dt=t1-t0;