time_3d.c

3D ray trace of elstic wave in homogenious medium

        return(ERR_MALLOC);
    if(!(t =(float ***)malloc((unsigned)nx*sizeof(float **)))){
        free((char *)hs);
        return(ERR_MALLOC);
    }
    if(!(longflags =(int *)malloc((unsigned)n1*sizeof(int)))){
        free((char *)t);
        free((char *)hs);
        return(ERR_MALLOC);
    }/* size of the largest side of the model */
    for(x=0;x<nx;x++)
        if(   !(hs[x]=(float **)malloc((unsigned)ny*sizeof(float *)))
           || !(t[x] =(float **)malloc((unsigned)ny*sizeof(float *)))
            ){
                free_ptrs(x);
                return(ERR_MALLOC);
            }
    for(x=0;x<nx;x++)
        for(y=0;y<ny;y++){
            hs[x][y]=hs_buf+x*np+y*nz;
            t[x][y]=t_buf+x*np+y*nz;
        }

/* stop here if recursive call */
    if(init_stage) return(NO_ERROR);

/* initialize all times as INFINITY if licit point source */
    if(fxs>=0.0 && fxs<=nx-1 && fys>=0 && fys<=ny-1 && fzs>=0 && fzs<=nz-1)
        for(x=0,pf=t_buf;x<nt;x++) *pf++=INFINITY;

/* assign INFINITY to hs in dummy meshes (x=nmesh_x|y=nmesh_y|z=nmesh_z) */
/* and keep masked values in hs_keep[].                                  */
    x=((nx+1)*(ny+1)+(nx+1)*nz+nz*ny)*sizeof(float);
    if(!(hs_keep=(float *)malloc((unsigned)x))) {
        free_ptrs(nx);
        return(ERR_MALLOC);
    }
    pf=hs_keep;
    for(x=0;x<nx;x++){
        for(y=0;y<ny;y++) {
            *pf++=hs[x][y][nmesh_z];
            hs[x][y][nmesh_z]=INFINITY;
        }
        for(z=0;z<nmesh_z;z++) {
            *pf++=hs[x][nmesh_y][z];
            hs[x][nmesh_y][z]=INFINITY;
        }
    }
    for(y=0;y<nmesh_y;y++)
        for(z=0;z<nmesh_z;z++) {
            *pf++=hs[nmesh_x][y][z];
            hs[nmesh_x][y][z]=INFINITY;
        }

/* test for negative slowness value */
    for(x=0,pf=hs_buf;x<nx*ny*nz;x++,pf++)
        if(*pf<=0.0){
	    printf("\n%d   %f",x,*pf);
            free_ptrs(nx);
            return(ERR_NONPHYSICAL);
        }/* a negative value would provoke an infinitely recursive call */
         /* and act as a "black hole" driving all times to -INFINITY !! */

    return(NO_ERROR);
}

/*------------------------------------------------Init_point()--------------*/

static int
init_point()

{
    int
        signal=NO_ERROR,
        x,y,z,
        test,
        t_X0,t_X1,t_Y0,t_Y1,t_Z0,t_Z1;
    float
        min_t,
        hs0,
        allowed_delta_hs,
        dist;

/* test relevance of source position or locate minimum time source point */
    if(fxs<0.0 || fxs>nx-1 || fys<0 || fys>ny-1 || fzs<0 || fzs>nz-1){
        for(x=0,min_t=INFINITY;x<nx;x++)
            for(y=0;y<ny;y++)
                for(z=0;z<nz;z++)
                    if(t[x][y][z]<min_t){
                        min_t=t[x][y][z];
                        xs=x;
                        ys=y;
                        zs=z;
                    }
        if(ISINF(min_t)) return(ERR_MULT);
        source_at_node=1;
        mult=1;
        if(SMALLTALK)
            printf("\nMultiple source starting at node [%d,%d,%d] at time %g.",
                xs,ys,zs,min_t);
    }

    else {
/* locate node closest to source */
        xs=NINT(fxs);
        ys=NINT(fys);
        zs=NINT(fzs);
        if(xs==fxs && ys==fys && zs==fzs){
            source_at_node=1;
            if(SMALLTALK) printf("\nSource located exactly at node [%d,%d,%d].",
                    xs,ys,zs);
        }
        mult=0;
    }

/* test relevance of slowness at the vicinity of the source */
/* (not tested in the case of multiple source)              */
    if(source_at_node){
        hs0=hs[xs][ys][zs];
        if(ISINF(hs0) && zs) hs0=hs[xs][ys][zs-1];
        if(ISINF(hs0) && ys) hs0=hs[xs][ys-1][zs];
        if(ISINF(hs0) && xs) hs0=hs[xs-1][ys][zs];
        if(ISINF(hs0) && zs && ys) hs0=hs[xs][ys-1][zs-1];
        if(ISINF(hs0) && zs && xs) hs0=hs[xs-1][ys][zs-1];
        if(ISINF(hs0) && ys && xs) hs0=hs[xs-1][ys-1][zs];
        if(ISINF(hs0) && zs && ys && xs) hs0=hs[xs-1][ys-1][zs-1];
    }
    else if(!mult){
        x=(fxs<xs) ? xs-1:xs;
        y=(fys<ys) ? ys-1:ys;
        z=(fzs<zs) ? zs-1:zs;
        hs0=hs[x][y][z];
        if(ISINF(hs0) && fxs==xs && xs){
            hs0=hs[x-1][y][z];
            if(ISINF(hs0) && fys==ys && ys) hs0=hs[x-1][y-1][z];
            if(ISINF(hs0) && fzs==zs && zs) hs0=hs[x-1][y][z-1];
        }
        if(ISINF(hs0) && fys==ys && ys){
            hs0=hs[x][y-1][z];
            if(ISINF(hs0) && fzs==zs && zs) hs0=hs[x][y-1][z-1];
        }
        if(ISINF(hs0) && fzs==zs && zs) hs0=hs[x][y][z-1];
    }
    if(ISINF(hs0)) return(ERR_INF);

/* if source is multiple, do not initialize at all the timefield */
    if(mult){
        X0=X1=xs;
        Y0=Y1=ys;
        Z0=Z1=zs;
        return(NO_ERROR);
    }/* this case has higher priority than the no_init directive */

/* if asked for, use minimal initialization : t=0.0 at the source point */
    if(no_init){
        X0=X1=xs;
        Y0=Y1=ys;
        Z0=Z1=zs;
        init_nearest();
        return(NO_ERROR);
    }

/* initialize inclusive boundaries of explored region */
    X0=max(xs-1,0);
    Y0=max(ys-1,0);
    Z0=max(zs-1,0);
    X1=min(xs+1,nmesh_x-1);
    Y1=min(ys+1,nmesh_y-1);
    Z1=min(zs+1,nmesh_z-1);

/* search largest parallelepipedic homogeneous box centered on the source */
    t_X0=t_X1=t_Y0=t_Y1=t_Z0=t_Z1=0;
    /* these flags will signal that a heterogeneity has been reached */
    allowed_delta_hs=hs0*hs_eps_init;
    /* defines tolerated inhomogeneity for exact initialization */
    do{
        test=0;
        if(X0 && !t_X0){
            test++;
            x=X0;
            for(y=Y0;y<=Y1 && y<nmesh_y && !t_X0;y++)
                for(z=Z0;z<=Z1 && z<nmesh_z && !t_X0;z++)
                    if(fabs(hs[x][y][z]-hs0)>allowed_delta_hs) t_X0=1;
            if(!t_X0) X0--;
        }
        if(Y0 && !t_Y0){
            test++;
            y=Y0;
            for(x=X0;x<=X1 && x<nmesh_x && !t_Y0;x++)
                for(z=Z0;z<=Z1 && z<nmesh_z && !t_Y0;z++)
                    if(fabs(hs[x][y][z]-hs0)>allowed_delta_hs) t_Y0=1;
            if(!t_Y0) Y0--;
        }
        if(Z0 && !t_Z0){
            test++;
            z=Z0;
            for(x=X0;x<=X1 && x<nmesh_x && !t_Z0;x++)
                for(y=Y0;y<=Y1 && y<nmesh_y && !t_Z0;y++)
                    if(fabs(hs[x][y][z]-hs0)>allowed_delta_hs) t_Z0=1;
            if(!t_Z0) Z0--;
        }
        if(X1<nmesh_x && !t_X1){
            test++;
            X1++;
            x=X1;
            for(y=Y0;y<=Y1 && y<nmesh_y && !t_X1;y++)
                for(z=Z0;z<=Z1 && z<nmesh_z && !t_X1;z++)
                    if(fabs(hs[x][y][z]-hs0)>allowed_delta_hs) t_X1=1;
        }
        if(Y1<nmesh_y && !t_Y1){
            test++;
            Y1++;
            y=Y1;
            for(x=X0;x<=X1 && x<nmesh_x && !t_Y1;x++)
                for(z=Z0;z<=Z1 && z<nmesh_z && !t_Y1;z++)
                    if(fabs(hs[x][y][z]-hs0)>allowed_delta_hs) t_Y1=1;
        }
        if(Z1<nmesh_z && !t_Z1){
            test++;
            Z1++;
            z=Z1;
            for(x=X0;x<=X1 && x<nmesh_x && !t_Z1;x++)
                for(y=Y0;y<=Y1 && y<nmesh_y && !t_Z1;y++)
                    if(fabs(hs[x][y][z]-hs0)>allowed_delta_hs) t_Z1=1;
        }
    } while(test);

    if(X0) X0++;
    if(Y0) Y0++;
    if(Z0) Z0++;
    if(X1<nmesh_x) X1--;
    if(Y1<nmesh_y) Y1--;
    if(Z1<nmesh_z) Z1--;
    /* 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()

/* 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(vx,vy,vz,xl,yl,zl)

float vx,vy,vz;
int xl,yl,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;