time_3d.c

射线追踪程序

/* copy args (with a few preliminary tests) to internal variables */
    hs_buf=HS;
    t_buf=T;
    nx=NX;
    ny=NY;
    nz=NZ;
    fxs=XS;
    fys=YS;
    fzs=ZS;
    hs_eps_init=HS_EPS_INIT;
    if(hs_eps_init<0.0 || hs_eps_init>1.0) {
        error(ERR_HS_EPS);
        return ERR_HS_EPS;
    }
    if(MSG<0){
        no_init=1;
        MSG= -MSG;
    }/* this trick inhibits search for homogeneous region around the source */
    messages=MSG;

/* compute */
    if((signal=pre_init())==NO_ERROR){
        signal=propagate_point(init_point());
        free_ptrs(nx);
    }
    if(init_stage==0 || signal!=NO_ERROR) error(signal);
    return signal;
}

/*---------------------------- Time_3d_(): FORTRAN INTERFACE ---------------*/
 
int time_3d_(float *HS, float *T, int *NZ, int *NY, int *NX,
             float *ZS, float *YS, float *XS, float *HS_EPS_INIT, int *MSG)
 
/* All FORTRAN arguments are pointers. Dimensions X and Z are */
/* swapped in order to mimic FORTRAN mapping conventions.     */
 
{
    int     signal;
 
#ifdef DEBUG_ARGS
    fprintf(stderr,"Entering function time_3d using Fortran-style interface\n");
    fprintf(stderr,"Dimensions nz=%d ny=%d nx=%d\n",*NZ,*NY,*NX);
    fprintf(stderr,"(arrays stored in memory as Fortran-array(nz,ny,nx))\n");
    fprintf(stderr,"Source location zs=%g ys=%g xs=%g\n",*ZS,*YS,*XS);
    fprintf(stderr,"Tolerance used (HS_EPS_INIT) %g\n",*HS_EPS_INIT);
    fprintf(stderr,"Verbosity (MSG) %d (%s)\n",*MSG,(*MSG)?"verbose":"silent");
    fflush(stderr);
#endif

    hs_buf=HS;
    t_buf=T;
    nx= *NX;
    ny= *NY;
    nz= *NZ;   
    fxs= *XS;
    fys= *YS;
    fzs= *ZS;
    hs_eps_init= *HS_EPS_INIT;
    if(hs_eps_init<0.0 || hs_eps_init>1.0) {
        error(ERR_HS_EPS);
        return ERR_HS_EPS;
    }
    if(*MSG<0){
        no_init=1;
        messages= -(*MSG);
    }
    else messages= *MSG;
 
    if((signal=pre_init())==NO_ERROR){
        signal=propagate_point(init_point());
        free_ptrs(nx);
    }
    if(init_stage==0 || signal!=NO_ERROR) error(signal);
    return signal;
}

/*------------------------------------------------Pre_init()----------------*/

static int pre_init(void)

{
    int
        x,y,z,
        np,nt,
        n0,n1;
    float
        *pf;

    nmesh_x=nx-1;
    nmesh_y=ny-1;
    nmesh_z=nz-1;
    np=ny*nz;
    nt=nx*np;
    n1=max(nx,ny);
    n0=max(nx,nz);
    if(n1==n0) n0=max(ny,nz);
    n1*=n0;

/* allocate pointers */
    if(!(hs=(float ***)malloc((unsigned)nx*sizeof(float **))))
        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;
	timeshift=0.0;
    }

/* 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){
            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(void)

{
    int
        signal=NO_ERROR,
        x,y,z,
        test,
        t_X0,t_X1,t_Y0,t_Y1,t_Z0,t_Z1;
    float
        min_t, max_t,
        hs0=0.0,           /* initialization required by gcc -Wall, unused */
        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=max_t=t[0][0][0];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(t[x][y][z]>max_t)
	                max_t=t[x][y][z];
	        }
        if(ISINF(min_t)) return ERR_MULTINF;
        if(min_t==max_t) return ERR_MULTNONE;
	if(init_stage==0){
	  if(min_t<0.0){
	    timeshift=min_t;
	    for(x=0;x<nx;x++)
	      for(y=0;y<ny;y++)
	        for(z=0;z<nz;z++)
	          if(!(ISINF(t[x][y][z]))) t[x][y][z]-=timeshift;
	  }/* shift all times to non-negative values */
	  else timeshift=0.0;
	}/* (done because fuzzy tests require all times to be non-negative) */
        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--;