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 intinit_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 voidinit_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 voidinit_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;