diff --git a/raytime/raytime.c b/raytime/raytime.c
index 414673e43627ba094f6f70bd2b34f64c08907ec1..4f05a55368aa4523e1c19f0cee555184b41c2e18 100644
--- a/raytime/raytime.c
+++ b/raytime/raytime.c
@@ -334,7 +334,7 @@ private (coordgx,irec,Time,Jr)
 						//fprintf(stderr,"ir=%d ix=%d iz=%d velocity=%f\n", ir, ix, iz, velocity[ix*mod.nz+iz]);
 						cp_average += velocity[ix*mod.nz+iz];
 					}
-					cp_average = cp_average/((float)nr);
+					cp_average = cp_average/((float)(nr-1));
             		ampl[ipos] = sqrt(time[ipos]*cp_average);
             		if (verbose>4) vmess("FD: shot=%f,%f receiver at %f(%d),%f(%d) T=%e V=%f Ampl=%f",coordsx.x, coordsx.z, coordgx.x, rec.x[irec], coordgx.z, rec.z[irec], time[ipos], cp_average, ampl[ipos]); 
         		}
diff --git a/raytime3d/Makefile b/raytime3d/Makefile
index 4aff110d5a69b98799cbb3274f7798c464f1f23e..1125552ed09f052253f4cc22fc24ab3cc28949a7 100644
--- a/raytime3d/Makefile
+++ b/raytime3d/Makefile
@@ -26,7 +26,7 @@ PRG = raytime3d
 
 SRCC	= $(PRG).c \
 		vidale3d.c \
-		src3d.c \
+		src3D.c \
 		getParameters3d.c  \
 		getWaveletInfo.c  \
 		writeSrcRecPos.c  \
diff --git a/raytime3d/getParameters3d.c b/raytime3d/getParameters3d.c
index af5b7609591faa6773d52d094ad37e89d5a748f4..13fc3737bbdf5393724e44e31f759dccc41fc6ad 100644
--- a/raytime3d/getParameters3d.c
+++ b/raytime3d/getParameters3d.c
@@ -26,7 +26,7 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0, float dx,
 
 long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, long verbose)
 {
-	long 	nx, nz, ny, nsrc, ix, iy, axis, is0;
+	long 	nx, nz, ny, nsrc, ix, iy, iz, axis, is0;
 	long 	n1, n2, n3; 
 	long	idzshot, idxshot, idyshot;
 	long	src_ix0, src_iz0, src_iy0, src_ix1, src_iz1, src_iy1;
@@ -137,20 +137,28 @@ long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPa
 	src_iz1=MAX(0,NINT((srcendz-sub_z0)/dz));
 	src_iz1=MIN(src_iz1,nz);
 
-	shot->y = (long *)calloc(shot->ny,sizeof(long));
-	shot->x = (long *)calloc(shot->nx,sizeof(long));
-	shot->z = (long *)calloc(shot->nz,sizeof(long));
-	for (is=0; is<shot->ny; is++) {
-		shot->y[is] = src_iy0+is*idyshot;
-		if (shot->y[is] > ny-1) shot->ny = is-1;
-	}
-	for (is=0; is<shot->nx; is++) {
-		shot->x[is] = src_ix0+is*idxshot;
-		if (shot->x[is] > nx-1) shot->nx = is-1;
-	}
-	for (is=0; is<shot->nz; is++) {
-        shot->z[is] = src_iz0+is*idzshot;
-        if (shot->z[is] > nz-1) shot->nz = is-1;
+	shot->y		= (long *)calloc(shot->nx*shot->ny*shot->nz,sizeof(long));
+	shot->x		= (long *)calloc(shot->nx*shot->ny*shot->nz,sizeof(long));
+	shot->z		= (long *)calloc(shot->nx*shot->ny*shot->nz,sizeof(long));
+	shot->ys	= (float *)calloc(shot->nx*shot->ny*shot->nz,sizeof(float));
+	shot->xs	= (float *)calloc(shot->nx*shot->ny*shot->nz,sizeof(float));
+	shot->zs	= (float *)calloc(shot->nx*shot->ny*shot->nz,sizeof(float));
+
+	for (iy=0; iy<shot->ny; iy++) {
+		for (ix=0; ix<shot->nx; ix++) {
+			for (iz=0; iz<shot->nz; iz++) {
+				shot->x[iz*shot->ny*shot->nx+iy*shot->nx+ix]	= src_ix0+ix*idxshot+1;
+				shot->y[iz*shot->ny*shot->nx+iy*shot->nx+ix]	= src_iy0+iy*idyshot+1;
+				shot->z[iz*shot->ny*shot->nx+iy*shot->nx+ix]	= src_iz0+iz*idzshot+1;
+				shot->xs[iz*shot->ny*shot->nx+iy*shot->nx+ix]	= xsrc+ix*dxshot;
+				shot->ys[iz*shot->ny*shot->nx+iy*shot->nx+ix]	= ysrc+iy*dyshot;
+				shot->zs[iz*shot->ny*shot->nx+iy*shot->nx+ix]	= zsrc+iz*dzshot;
+
+				if (shot->x[iz*shot->ny*shot->nx+iy*shot->nx+ix] > nx-1) shot->nx = ix-1;
+				if (shot->y[iz*shot->ny*shot->nx+iy*shot->nx+ix] > ny-1) shot->ny = iy-1;
+				if (shot->z[iz*shot->ny*shot->nx+iy*shot->nx+ix] > nz-1) shot->nz = iz-1;
+			}
+		}
     }
 
 	/* check if source array is defined */
@@ -323,6 +331,7 @@ long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPa
 			vmess("************* receiver info ***************");
 			vmess("*******************************************");
 			vmess("ntrcv   = %li nrcv    = %li ", rec->nt, rec->n);
+			vmess("nxrcv   = %li nyrcv   = %li ", rec->nx, rec->ny);
 			vmess("dzrcv   = %f dxrcv   = %f dyrcv   = %f ", dzrcv, dxrcv, dyrcv);
 			vmess("Receiver array at coordinates: ");
 			vmess("zmin    = %f zmax    = %f ", rec->zr[0]+sub_z0, rec->zr[rec->n-1]+sub_z0);
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
index fc2d685b100d66169c126b6d48855dd29f95df5a..7d5e65a0ead722ade4f4634ffe5bdfe991006eca 100644
--- a/raytime3d/raytime3d.c
+++ b/raytime3d/raytime3d.c
@@ -33,7 +33,7 @@ long writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
 
 void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, long xs, long ys, long zs, long NCUBE);
 
-void src3d(float *time0, float *slow0, long nz, long nx, long ny, float h, float ox, float oy, float oz, long *pxs, long *pys, long *pzs, long *cube);
+void src3D(float *time0, float *slow0, long nz, long nx, long ny, float h, float ox, float oy, float oz, long xs, long ys, long zs, long *cube);
 
 
 /*********************** self documentation **********************/
@@ -119,16 +119,19 @@ void main(int argc, char *argv[])
 		ny,			/* y-dimension of mesh (FRONT-TO-BACK) */
 		nz,			/* z-dimension of mesh  (TOP-TO-BOTTOM) */
 		nxy, nxyz, xs, ys, zs, cube,
-		xx, yy, zz,	i, j, k;
+		nx2, ny2, nz2, nxy2, nxyz2,
+		nrx, nry, nrz, nr, ix, iy, iz,
+		xx, yy, zz,	i, j, k, is, writer;
 	float
 		h,		/* spatial mesh interval (units consistant with vel) */
-		*slow0, *time0;
+		cp_average,
+		*slow0, *time0, *time1, *ampl;
 
 /* to read the velocity file */
 	long    error, n1, n2, n3, ret, size, nkeys, verbose, nwrite;
 	float	d1, d2, d3, f1, f2, f3, *tmpdata, c, scl, ox, oz, oy;
-	char	*file_cp, *file_out;
-	FILE	*fpt;
+	char	*file_cp, *file_out, file_time[100], file_amp[100];
+	FILE	*fpt, *fpa;
 	segy	*hdrs;
 
 /*---------------------------------------------------------------------------*
@@ -182,7 +185,7 @@ void main(int argc, char *argv[])
 		}
 
 		h = mod.dx;
-		slow0 = (float *)malloc(nz*nx*ny*sizeof(float));
+		slow0 = (float *)malloc((nz+2)*(nx+2)*(ny+2)*sizeof(float));
 		if (slow0 == NULL) verr("Out of memory for slow0 array!");
 
 		readModel3d(mod.file_cp, slow0, nz, nx, ny, h, verbose);
@@ -211,141 +214,131 @@ void main(int argc, char *argv[])
 		}
 	}
 
-	nxy = nx * ny;
-	nxyz = nx * ny * nz;
+	nxy = (nx+2) * (ny+2);
+	nxyz = (nx+2) * (ny+2) * (nz+2);
 
-	/* ALLOCATE MAIN GRID FOR TIMES */
-	time0 = (float *) malloc(sizeof(float)*nxyz);
-	if(time0 == NULL) verr("error in allocation of array time0");
-
-/*---------------------------------------------------------------------------*
- *  Input the source locations.
- *          and
- *  Initialize the traveltime array.  Place t=0 at source position.
- *---------------------------------------------------------------------------*/
-
-	src3d(time0, slow0, nz, nx, ny, h, ox, oy, oz, &xs, &ys, &zs, &cube);
-	if (verbose) vmess("source positions xs = %li ys = %li zs = %li", xs,ys,zs);
+	strcpy(file_time, file_out);
+    name_ext(file_time, "_time");
+	fpt = fopen(file_time, "w");
+    assert(fpt != NULL);
 
-/*	for (zz = 0; zz < nz; zz++) {
-		for (yy = 0; yy < ny; yy++) {
-			for (xx = 0; xx < nx; xx++) 
-				if (time0[zz*nxy+yy*nx+xx] != 1e10) fprintf(stderr,"slow[%li,%li,%li] = %f\n", xx,yy,zz, time0[zz*nxy+yy*nx+xx]);
-		}
+	if (ray.geomspread==1) {
+		strcpy(file_amp, file_out);
+		name_ext(file_amp, "_amp");
+		fpa = fopen(file_amp, "w");
+		assert(fpa != NULL);
 	}
-*/
-
-/*---------------------------------------------------------------------------*
- *  Read in receiver positions
- *---------------------------------------------------------------------------*/
 
-/*---------------------------------------------------------------------------*
- *  Check and set parameters
- *---------------------------------------------------------------------------*/
-
-/*---------------------------------------------------------------------------*
- *  Compute traveltimes.
- *---------------------------------------------------------------------------*/
-
-	vidale3d(slow0, time0, nz, nx, ny, h, xs, ys, zs, cube);
-
-/*---------------------------------------------------------------------------*
- *  Write output
- *---------------------------------------------------------------------------*/
+	if (verbose>2 && ray.geomspread==1) vmess("Computing geometrical spreading factor");
+
+	writer = 0;
+
+#pragma omp parallel for schedule(static,1) default(shared) \
+private (is,time0,ampl,nrx,nry,nrz,nr,cp_average,i,j,k,ix,iy,iz,hdrs,tmpdata,nwrite) 
+	for (is = 0; is < shot.n; is++) {
+
+		/* ALLOCATE MAIN GRID FOR TIMES */
+		time0 = (float *) malloc(sizeof(float)*nxyz);
+		if(time0 == NULL) verr("error in allocation of array time0");
+
+		/*---------------------------------------------------------------------------*
+		*  Input the source locations.
+		*          and
+		*  Initialize the traveltime array.  Place t=0 at source position.
+		*---------------------------------------------------------------------------*/
+
+		src3D(time0, slow0, nz+2, nx+2, ny+2, h, shot.xs[is]-ox+d2, shot.ys[is]-oy+d3, shot.zs[is]-oz+d1, shot.x[is], shot.y[is], shot.z[is], &cube);
+
+		/*---------------------------------------------------------------------------*
+		*  Compute traveltimes.
+		*---------------------------------------------------------------------------*/
+
+		vidale3d(slow0, time0, nz+2, nx+2, ny+2, h, shot.x[is], shot.y[is], shot.z[is], cube);
+
+		/*---------------------------------------------------------------------------*
+		*  Compute geometrical spreading.
+		*---------------------------------------------------------------------------*/
+
+		if (ray.geomspread==1) {
+			ampl = (float *) malloc(sizeof(float)*rec.n);
+			for (i = 0; i < rec.n; i++) {
+				/* compute average velocity between source and receiver */
+				nrx = (rec.x[i]-(shot.x[is]-1));
+				nry = (rec.y[i]-(shot.y[is]-1));
+				nrz = (rec.z[i]-(shot.z[is]-1));
+				nr = abs(nrx) + abs(nry) + abs(nrz);
+				cp_average = 0.0;
+				for (j=0; j<nr; j++) {
+					ix = shot.x[is] + floor((j*nrx)/nr);
+					iy = shot.y[is] + floor((j*nry)/nr);
+					iz = shot.z[is] + floor((j*nrz)/nr);
+					cp_average += 1.0/slow0[iz*nxy+iy*(nx+2)+ix];
+				}
+				cp_average = cp_average/((float)(nr-1));
+				ampl[i] = (time0[iz*nxy+iy*(nx+2)+ix]*cp_average);
+			}
+		}
 
+		/*---------------------------------------------------------------------------*
+		*  Write output
+		*---------------------------------------------------------------------------*/
 
-	for (zz = 0; zz < nz; zz++) {
-		for (yy = 0; yy < ny; yy++) {
-			for (xx = 0; xx < nx; xx++) {
-				if (time0[zz*nxy+yy*nx+xx] == 1e10) time0[zz*nxy+yy*nx+xx]=0.0;
-			}
+		while (writer < is) {
+			#pragma omp flush(writer)
 		}
-	}
 
-//	ret = open_file(file_out, GUESS_TYPE, DELPHI_CREATE);
-//	if (ret < 0 ) verr("error in creating output file %s", file_out);
-	fpt = fopen(file_out, "w");
-    assert(fpt != NULL);
+		if (verbose>2) vmess("Writing src %li of %li sources",is+1,shot.n);
+		if (verbose) vmess("xsrc[%li]=%f ysrc[%li]=%f zsrc[%li]=%f",shot.x[is],shot.xs[is],shot.y[is],shot.ys[is],shot.z[is],shot.zs[is]);
 
-	hdrs = (segy *) malloc(1*sizeof(segy));
-	tmpdata = (float *)malloc(nz*sizeof(float));
-	// f1   = ox;
-	// f2   = oy;
-	// d1   = h;
-	// d2   = h;
-
-//	gen_hdrs(hdrs,nx,ny,f1,f2,d1,d2,TRID_ZX);
-	// for (i = 0; i < ny; i++) {
-	// 	hdrs[i].tracl	= i+1;
-	// 	hdrs[i].tracf	= i+1;
-	// 	hdrs[i].scalco	= -1000;
-	// 	hdrs[i].scalel	= -1000;
-	// 	hdrs[i].sx		= (long)(ox+xs*h)*1000;
-	// 	hdrs[i].sy		= (long)(oy+ys*h)*1000;
-	// 	hdrs[i].gy		= (long)(oy+i*d2)*1000;
-	// 	hdrs[i].sdepth	= (long)(oz+zs*h)*1000;
-	// 	hdrs[i].selev	= (long)(oz+zs*h)*-1000;
-	// 	hdrs[i].ns 		= nx;
-	// 	hdrs[i].ntr		= ny;
-	// 	hdrs[i].trwf	= ny;
-	// 	hdrs[i].f1		= mod.x0;
-	// 	hdrs[i].f2		= mod.y0;
-	// 	hdrs[i].dt 		= mod.dx*1e6;
-	// 	hdrs[i].d1 		= mod.dx;
-	// 	hdrs[i].d2 		= mod.dy;
-	// 	hdrs[i].fldr	= 1;
-	// 	hdrs[i].trwf	= ny;
-	// 	for (j = 0; j < nx; j++) {
-	// 		tmpdata[i*nx+j] = time0[i*nx+j];
-	// 	}
-	// 	nwrite = fwrite( &hdrs[i], 1, TRCBYTES, fpt);
-	// 	assert(nwrite == TRCBYTES);
-    //     nwrite = fwrite( &tmpdata[i*nx], sizeof(float), nx, fpt );
-	// 	assert(nwrite == nx);
-	// }
-	for (i = 0; i < ny; i++) {
-		for (j = 0; j < nx; j++) {
-			hdrs[0].tracl	= i*nx+j+1;
-			hdrs[0].tracf	= i*nx+j+1;
+		hdrs = (segy *) calloc(1,sizeof(segy));
+		tmpdata = (float *)malloc((rec.nx)*sizeof(float));
+		
+		for (i = 0; i < rec.ny; i++) {
+			hdrs[0].fldr	= is+1;
+			hdrs[0].tracl	= i+1;
+			hdrs[0].tracf	= i+1;
 			hdrs[0].scalco	= -1000;
 			hdrs[0].scalel	= -1000;
-			hdrs[0].sx		= (long)(f2+j*d2)*1000;
-			hdrs[0].gx		= (long)(f2+j*d2)*1000;
-			hdrs[0].sy		= (long)(f3+i*d3)*1000;
-			hdrs[0].gy		= (long)(f3+i*d3)*1000;
-			hdrs[0].ns 		= nz;
-			hdrs[0].ntr		= ny*nx;
-			hdrs[0].trwf	= nx;
-			hdrs[0].f1		= f1;
-			hdrs[0].f2		= f2;
-			hdrs[0].dt 		= (int)(d1*1e6);
-			hdrs[0].d1 		= d1;
-			hdrs[0].d2 		= d2;
-			hdrs[0].fldr	= 1;
-
-			for (k = 0; k < nz; k++) {
-				tmpdata[k] = time0[k*nxy+i*nx+j];
+			hdrs[0].sx		= (long)(f2+(shot.x[is]-1)*d2)*1000;
+			hdrs[0].sy		= (long)(f3+(shot.y[is]-1)*d3)*1000;
+			hdrs[0].sdepth	= (long)(f1+(shot.z[is]-1)*d1)*1000;
+			hdrs[0].selev	= -(long)(f1+(shot.z[is]-1)*d1)*1000;
+			hdrs[0].gy		= mod.y0+rec.yr[i*rec.nx];
+			hdrs[0].ns 		= rec.nx;
+			hdrs[0].ntr		= rec.ny*shot.n;
+			hdrs[0].trwf	= rec.ny*shot.n;
+			hdrs[0].f1		= mod.x0+rec.xr[0];
+			hdrs[0].f2		= mod.y0+rec.yr[0];
+			hdrs[0].dt 		= (long)(d2*1e6);
+			hdrs[0].d1 		= rec.xr[1]-rec.xr[0];
+			hdrs[0].d2 		= rec.yr[rec.nx]-rec.yr[0];
+
+			for (k = 0; k < rec.nx; k++) {
+				tmpdata[k] = time0[(rec.z[i*rec.nx+k]+1)*nxy+(rec.y[i*rec.nx+k]+1)*(nx+2)+rec.x[i*rec.nx+k]+1];
 			}
 
 			nwrite = fwrite( &hdrs[0], 1, TRCBYTES, fpt);
 			assert(nwrite == TRCBYTES);
-			nwrite = fwrite( tmpdata, sizeof(float), nz, fpt );
-			assert(nwrite == nz);
+			nwrite = fwrite( tmpdata, sizeof(float), rec.nx, fpt );
+			assert(nwrite == rec.nx);
+
+			if (ray.geomspread==1) {
+				nwrite = fwrite( &hdrs[0], 1, TRCBYTES, fpa);
+				assert(nwrite == TRCBYTES);
+				nwrite = fwrite( &ampl[i*rec.nx], sizeof(float), rec.nx, fpa );
+				assert(nwrite == rec.nx);
+			}
 		}
+		writer++;
+		free(time0);
+		free(hdrs);
+		free(tmpdata);
+		if (ray.geomspread==1) free(ampl);
 	}
-	fclose(fpt);
-
-/*
-	ret = write_data(file_out,tmpdata,nx,ny,f1,f2,d1,d2,type,hdrs);
-	if (ret < 0 ) verr("error on writing output file.");
-	ret = close_file(file_out);
-	if (ret < 0) verr("err %li on closing output file",ret);
-*/
 
-	free(time0);
+	if (ray.geomspread==1) fclose(fpa);
+	fclose(fpt);
 	free(slow0);
-	free(hdrs);
-	free(tmpdata);
 
 	exit(0);
 
diff --git a/raytime3d/raytime3d.h b/raytime3d/raytime3d.h
index 39355d4f07d7f5c3ed2baa5734a53893793a5cdc..dbf479ac40df596daebffefb34759c8922b19538 100644
--- a/raytime3d/raytime3d.h
+++ b/raytime3d/raytime3d.h
@@ -115,6 +115,9 @@ typedef struct _shotPar { /* Shot Parameters */
 	long *z;
 	long *x;
 	long *y;
+	float *zs;
+	float *xs;
+	float *ys;
 } shotPar;
 
 typedef struct _raypar { /* ray-tracing parameters */
diff --git a/raytime3d/readModel3d.c b/raytime3d/readModel3d.c
index 7615609574090f903be61113cfd079a671a67f7d..c3e325dd59ef15002140d80329a38d3668d1fad5 100644
--- a/raytime3d/readModel3d.c
+++ b/raytime3d/readModel3d.c
@@ -34,7 +34,7 @@ long readModel3d(char *file_name, float *slowness, long nz, long nx, long ny, fl
 	float *tmp;
     segy hdr;
 
-	nxy = nx * ny;
+	nxy = (nx+2) * (ny+2);
 	tmp = (float *)malloc(nz*sizeof(float));
 
 /* open files and read first header */
@@ -43,7 +43,6 @@ long readModel3d(char *file_name, float *slowness, long nz, long nx, long ny, fl
    	assert( fpcp != NULL);
    	nread = fread(&hdr, 1, TRCBYTES, fpcp);
    	assert(nread == TRCBYTES);
-	vmess("nz=%li hdr.ns=%li",nz,hdr.ns);
 	assert(hdr.ns == nz);
 
 	if (nx==1) { /* 1D model */
@@ -51,7 +50,7 @@ long readModel3d(char *file_name, float *slowness, long nz, long nx, long ny, fl
 		for (j = 0; j < nz; j++) {
 			for (k = 0; k < ny; k++) {
 				for (i = 0; i < nx; i++) {
-					slowness[j*nxy+k*nx+i] = h/tmp[j];
+					slowness[(j+1)*nxy+(k+1)*(nx+2)+(i+1)] = h/tmp[j];
 				}
 			}
 		}
@@ -61,7 +60,7 @@ long readModel3d(char *file_name, float *slowness, long nz, long nx, long ny, fl
        		nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
 			for (j = 0; j < nz; j++) {
 				for (k = 0; k < ny; k++) {
-					slowness[j*nxy+k*nx+i] = h/tmp[j];
+					slowness[(j+1)*nxy+(k+1)*(nx+2)+(i+1)] = h/tmp[j];
 				}
 			}
        		nread = fread(&hdr, 1, TRCBYTES, fpcp);
@@ -73,13 +72,123 @@ long readModel3d(char *file_name, float *slowness, long nz, long nx, long ny, fl
 			for (i = 0; i < nx; i++) {
        			nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
 				for (j = 0; j < nz; j++) {
-					slowness[j*nxy+k*nx+i] = h/tmp[j];
+					slowness[(j+1)*nxy+(k+1)*(nx+2)+(i+1)] = h/tmp[j];
 				}
        			nread = fread(&hdr, 1, TRCBYTES, fpcp);
 			}
 		}
 	}
 
+	/* Fill in the sides */
+	k=0;
+	for (j = 0; j < nz; j++) {
+		for (i = 0; i < nx; i++) {
+			slowness[(j+1)*nxy+k*(nx+2)+(i+1)] = slowness[(j+1)*nxy+(k+1)*(nx+2)+(i+1)];
+		}
+	}
+	k=ny+1;
+	for (j = 0; j < nz; j++) {
+		for (i = 0; i < nx; i++) {
+			slowness[(j+1)*nxy+k*(nx+2)+(i+1)] = slowness[(j+1)*nxy+(k-1)*(nx+2)+(i+1)];
+		}
+	}
+
+	j=0;
+	for (k = 0; k < ny; k++) {
+		for (i = 0; i < nx; i++) {
+			slowness[j*nxy+(k+1)*(nx+2)+(i+1)] = slowness[(j+1)*nxy+(k+1)*(nx+2)+(i+1)];
+		}
+	}
+	j=nz+1;
+	for (k = 0; k < ny; k++) {
+		for (i = 0; i < nx; i++) {
+			slowness[j*nxy+(k+1)*(nx+2)+(i+1)] = slowness[(j-1)*nxy+(k+1)*(nx+2)+(i+1)];
+		}
+	}
+
+	i=0;
+	for (k = 0; k < ny; k++) {
+		for (j = 0; j < nz; j++) {
+			slowness[(j+1)*nxy+(k+1)*(nx+2)+i] = slowness[(j+1)*nxy+(k+1)*(nx+2)+(i+1)];
+		}
+	}
+	i=nx+1;
+	for (k = 0; k < ny; k++) {
+		for (j = 0; j < nz; j++) {
+			slowness[(j+1)*nxy+(k+1)*(nx+2)+i] = slowness[(j+1)*nxy+(k+1)*(nx+2)+(i-1)];
+		}
+	}
+
+	/* Fill in corners of sides */
+	k=0;	j=0;
+	for (i = 0; i < nx; i++) {
+		slowness[j*nxy+k*(nx+2)+(i+1)] = (1.0/2.0)*(slowness[(j+1)*nxy+k*(nx+2)+(i+1)]+slowness[j*nxy+(k+1)*(nx+2)+(i+1)]);
+	}
+	k=ny+1;	j=0;
+	for (i = 0; i < nx; i++) {
+		slowness[j*nxy+k*(nx+2)+(i+1)] = (1.0/2.0)*(slowness[(j+1)*nxy+k*(nx+2)+(i+1)]+slowness[j*nxy+(k-1)*(nx+2)+(i+1)]);
+	}
+	k=0;	j=nz+1;
+	for (i = 0; i < nx; i++) {
+		slowness[j*nxy+k*(nx+2)+(i+1)] = (1.0/2.0)*(slowness[(j-1)*nxy+k*(nx+2)+(i+1)]+slowness[j*nxy+(k+1)*(nx+2)+(i+1)]);
+	}
+	k=ny+1;	j=nz+1;
+	for (i = 0; i < nx; i++) {
+		slowness[j*nxy+k*(nx+2)+(i+1)] = (1.0/2.0)*(slowness[(j-1)*nxy+k*(nx+2)+(i+1)]+slowness[j*nxy+(k-1)*(nx+2)+(i+1)]);
+	}
+
+	k=0;	i=0;
+	for (j = 0; j < nz; j++) {
+		slowness[(j+1)*nxy+k*(nx+2)+i] = (1.0/2.0)*(slowness[(j+1)*nxy+k*(nx+2)+(i+1)]+slowness[(j+1)*nxy+(k+1)*(nx+2)+i]);
+	}
+	k=ny+1;	i=0;
+	for (j = 0; j < nz; j++) {
+		slowness[(j+1)*nxy+k*(nx+2)+i] = (1.0/2.0)*(slowness[(j+1)*nxy+k*(nx+2)+(i+1)]+slowness[(j+1)*nxy+(k-1)*(nx+2)+i]);
+	}
+	k=0;	i=nx+1;
+	for (j = 0; j < nz; j++) {
+		slowness[(j+1)*nxy+k*(nx+2)+i] = (1.0/2.0)*(slowness[(j+1)*nxy+k*(nx+2)+(i-1)]+slowness[(j+1)*nxy+(k+1)*(nx+2)+i]);
+	}
+	k=ny+1;	i=nx+1;
+	for (j = 0; j < nz; j++) {
+		slowness[(j+1)*nxy+k*(nx+2)+i] = (1.0/2.0)*(slowness[(j+1)*nxy+k*(nx+2)+(i-1)]+slowness[(j+1)*nxy+(k-1)*(nx+2)+i]);
+	}
+
+	j=0;	i=0;
+	for (k = 0; k < ny; k++) {
+		slowness[j*nxy+(k+1)*(nx+2)+i] = (1.0/2.0)*(slowness[j*nxy+(k+1)*(nx+2)+(i+1)]+slowness[(j+1)*nxy+(k+1)*(nx+2)+i]);
+	}
+	j=nz+1;	i=0;
+	for (k = 0; k < ny; k++) {
+		slowness[j*nxy+(k+1)*(nx+2)+i] = (1.0/2.0)*(slowness[j*nxy+(k+1)*(nx+2)+(i+1)]+slowness[(j-1)*nxy+(k+1)*(nx+2)+i]);
+	}
+	j=0;	i=nx+1;
+	for (k = 0; k < ny; k++) {
+		slowness[j*nxy+(k+1)*(nx+2)+i] = (1.0/2.0)*(slowness[j*nxy+(k+1)*(nx+2)+(i-1)]+slowness[(j+1)*nxy+(k+1)*(nx+2)+i]);
+	}
+	j=nz+1;	i=nx+1;
+	for (k = 0; k < ny; k++) {
+		slowness[j*nxy+(k+1)*(nx+2)+i] = (1.0/2.0)*(slowness[j*nxy+(k+1)*(nx+2)+(i-1)]+slowness[(j-1)*nxy+(k+1)*(nx+2)+i]);
+	}
+
+	/* Fill in corners of cubes */
+	i=0;	j=0;	k=0;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j+1)*nxy+(k+1)*(nx+2)+i]+slowness[j*nxy+(k+1)*(nx+2)+(i+1)]+slowness[(j+1)*nxy+k*(nx+2)+(i+1)]);
+	i=nx+1;	j=0;	k=0;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j+1)*nxy+(k+1)*(nx+2)+i]+slowness[j*nxy+(k+1)*(nx+2)+(i-1)]+slowness[(j+1)*nxy+k*(nx+2)+(i-1)]);
+	i=0;	j=nz+1;	k=0;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j-1)*nxy+(k+1)*(nx+2)+i]+slowness[j*nxy+(k+1)*(nx+2)+(i+1)]+slowness[(j-1)*nxy+k*(nx+2)+(i+1)]);
+	i=0;	j=0;	k=ny+1;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j+1)*nxy+(k-1)*(nx+2)+i]+slowness[j*nxy+(k-1)*(nx+2)+(i+1)]+slowness[(j+1)*nxy+k*(nx+2)+(i+1)]);
+	i=nx+1;	j=nz+1;	k=0;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j-1)*nxy+(k+1)*(nx+2)+i]+slowness[j*nxy+(k+1)*(nx+2)+(i-1)]+slowness[(j-1)*nxy+k*(nx+2)+(i-1)]);
+	i=nx+1;	j=0;	k=ny+1;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j+1)*nxy+(k-1)*(nx+2)+i]+slowness[j*nxy+(k-1)*(nx+2)+(i-1)]+slowness[(j+1)*nxy+k*(nx+2)+(i-1)]);
+	i=0;	j=nz+1;	k=ny+1;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j-1)*nxy+(k-1)*(nx+2)+i]+slowness[j*nxy+(k-1)*(nx+2)+(i+1)]+slowness[(j-1)*nxy+k*(nx+2)+(i+1)]);
+	i=nx+1;	j=nz+1;	k=ny+1;
+	slowness[j*nxy+k*(nx+2)+i] = (1.0/3.0)*(slowness[(j-1)*nxy+(k-1)*(nx+2)+i]+slowness[j*nxy+(k-1)*(nx+2)+(i-1)]+slowness[(j-1)*nxy+k*(nx+2)+(i-1)]);
+
    	fclose(fpcp);
 
     return 0;
diff --git a/raytime3d/readModel3d_backup.c b/raytime3d/readModel3d_backup.c
deleted file mode 100644
index c180042373fb7c4af9aeeb9d5d632fbe3b2df9b8..0000000000000000000000000000000000000000
--- a/raytime3d/readModel3d_backup.c
+++ /dev/null
@@ -1,87 +0,0 @@
-#define _FILE_OFFSET_BITS 64
-#define _LARGEFILE_SOURCE
-#define _LARGEFILE64_SOURCE
-
-#include <assert.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <errno.h>
-#include <math.h>
-#include "segy.h"
-#include "par.h"
-#include "raytime3d.h"
-
-#define     MAX(x,y) ((x) > (y) ? (x) : (y))
-#define     MIN(x,y) ((x) < (y) ? (x) : (y))
-#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
-
-/**
-*  Reads gridded model files and compute from them medium parameters used in the FD kernels.
-*  The files read in contain the P (and S) wave velocity and density.
-*  The medium parameters calculated are lambda, mu, lambda+2mu, and 1/ro.
-*
-*   AUTHOR:
-*           Jan Thorbecke (janth@xs4all.nl)
-*           The Netherlands 
-**/
-
-long readModel3d(char *file_name, float *slowness, long nz, long nx, long ny, float h, long verbose)
-{
-    FILE    *fpcp;
-    size_t  nread;
-    long i, j, k, tracesToDo;
-	long nxy, nxz;
-	float *tmp;
-    segy hdr;
-
-	nxy = nx * ny;
-	tmp = (float *)malloc(nz*sizeof(float));
-
-/* open files and read first header */
-
-   	fpcp = fopen( file_name, "r" );
-   	assert( fpcp != NULL);
-   	nread = fread(&hdr, 1, TRCBYTES, fpcp);
-   	assert(nread == TRCBYTES);
-	assert(hdr.ns == nz);
-
-	if (nx==1) { /* 1D model */
-       	nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
-		for (j = 0; j < nz; j++) {
-			for (k = 0; k < ny; k++) {
-				for (i = 0; i < nx; i++) {
-					slowness[j*nxy+k*nx+i] = h/tmp[j];
-				}
-			}
-		}
-	}
-	else if (ny==1) { /* 2D model */
-		for (i = 0; i < nx; i++) {
-       		nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
-			for (j = 0; j < nz; j++) {
-				for (k = 0; k < ny; k++) {
-					slowness[j*nxy+k*nx+i] = h/tmp[j];
-				}
-			}
-       		nread = fread(&hdr, 1, TRCBYTES, fpcp);
-		}
-	}
-	else { /* Full 3D model */
-		/* read all traces */
-		for (k = 0; k < ny; k++) {
-			for (i = 0; i < nx; i++) {
-       			nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
-				for (j = 0; j < nz; j++) {
-					slowness[j*nxy+k*nx+i] = h/tmp[j];
-				}
-       			nread = fread(&hdr, 1, TRCBYTES, fpcp);
-			}
-		}
-	}
-
-   	fclose(fpcp);
-
-    return 0;
-}
-
-
diff --git a/raytime3d/recvPar3d.c b/raytime3d/recvPar3d.c
deleted file mode 100644
index 5776d7b49d4f9dcf948c91856551ce480afb1e4b..0000000000000000000000000000000000000000
--- a/raytime3d/recvPar3d.c
+++ /dev/null
@@ -1,617 +0,0 @@
-#include <stdio.h>
-#include <assert.h>
-#include <math.h>
-
-#include "raytime3d.h"
-#include "par.h"
-
-#define MIN(x,y) ((x) < (y) ? (x) : (y))
-#define MAX(x,y) ((x) > (y) ? (x) : (y))
-#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
-
-/**
-*  Calculates the receiver positions based on the input parameters
-*
-*   AUTHOR:
-*           Jan Thorbecke (janth@xs4all.nl)
-*
-*   Ammendments:
-*           Max Holicki changing the allocation receiver array (2-2016)
-*           The Netherlands 
-**/
-
-
-void name_ext(char *filename, char *extension);
-
-long recvPar3d(recPar *rec, float sub_x0, float sub_z0, float sub_y0, 
-float dx, float dz, float dy, long nx, long nz, long ny)
-{
-	float *yrcv1, *yrcv2, *xrcv1, *xrcv2, *zrcv1, *zrcv2;
-	long   i, ix, iy, ir, verbose;
-	float dxrcv, dzrcv, dyrcv, *dxr, *dzr, *dyr;
-	float rrcv, dphi, oxrcv, ozrcv, oyrcv, arcv;
-	double circ, h, a, b, e, s, xr, zr, yr, dr, srun, phase;
-	float xrange, zrange, yrange, sub_x1, sub_z1, sub_y1;
-	long Nx1, Nx2, Nz1, Nz2, Ny1, Ny2, Ndx, Ndz, Ndy, iarray, nrec, nh;
-	long nxrcv, nzrcv, nyrcv, ncrcv, nrcv, ntrcv, *nlrcv;
-	float *xrcva, *zrcva, *yrcva;
-	char* rcv_txt;
-	FILE *fp;
-
-	if (!getparlong("verbose", &verbose)) verbose = 0;
-
-    /* Calculate Model Dimensions */
-    sub_x1=sub_x0+(nx-1)*dx;
-    sub_z1=sub_z0+(nz-1)*dz;
-    sub_y1=sub_y0+(ny-1)*dy;
-
-/* Compute how many receivers are defined and then allocate the receiver arrays */
-
-    /* Receivers on a Circle */
-    if (getparfloat("rrcv",&rrcv)) {
-        if (!getparfloat("dphi",&dphi)) dphi=2.0;
-        ncrcv=NINT(360.0/dphi);
-        if (verbose) vmess("Total number of receivers on a circle: %li",ncrcv);
-    } 
-	else {
-		ncrcv=0;
-	}
-
-    /* Receivers from a File */
-    ntrcv=0;
-    if (!getparstring("rcv_txt",&rcv_txt)) rcv_txt=NULL;
-    if (rcv_txt!=NULL) {
-        /* Open text file */
-        fp=fopen(rcv_txt,"r");
-        assert(fp!=NULL);
-        /* Get number of lines */
-        while (!feof(fp)) if (fgetc(fp)=='\n') ntrcv++;
-        fseek(fp,-1,SEEK_CUR);
-        if (fgetc(fp)!='\n') ntrcv++; /* Checks if last line terminated by /n */
-        if (verbose) vmess("Number of receivers in rcv_txt file: %li",ntrcv);
-        rewind(fp);
-    }
-
-    /* Receiver Array */
-    nyrcv=countparval("yrcva");
-    nxrcv=countparval("xrcva");
-    nzrcv=countparval("zrcva");
-    if (nxrcv!=nzrcv) verr("Number of receivers in array xrcva (%li), zrcva(%li) are not equal",nxrcv,nzrcv);
-    if (nxrcv!=nyrcv) verr("Number of receivers in array xrcva (%li), yrcva(%li) are not equal",nxrcv,nyrcv);
-    if (verbose&&nxrcv) vmess("Total number of array receivers: %li",nxrcv);
-
-    /* Linear Receiver Arrays */
-	Ny1 = countparval("yrcv1");
-	Ny2 = countparval("yrcv2");
-	Nx1 = countparval("xrcv1");
-	Nx2 = countparval("xrcv2");
-	Nz1 = countparval("zrcv1");
-	Nz2 = countparval("zrcv2");
-    if (Ny1!=Ny2) verr("Number of receivers starting points in 'yrcv1' (%li) and number of endpoint in 'yrcv2' (%li) are not equal",Ny1,Ny2);
-    if (Nx1!=Nx2) verr("Number of receivers starting points in 'xrcv1' (%li) and number of endpoint in 'xrcv2' (%li) are not equal",Nx1,Nx2);
-    if (Nz1!=Nz2) verr("Number of receivers starting points in 'zrcv1' (%li) and number of endpoint in 'zrcv2' (%li) are not equal",Nz1,Nz2);
-    if (Nx1!=Nz2) verr("Number of receivers starting points in 'xrcv1' (%li) and number of endpoint in 'zrcv2' (%li) are not equal",Nx1,Nz2);
-    if (Nx1!=Ny2) verr("Number of receivers starting points in 'xrcv1' (%li) and number of endpoint in 'yrcv2' (%li) are not equal",Nx1,Ny2);
-
-    rec->max_nrec=ncrcv+ntrcv+nxrcv;
-
-	/* no receivers are defined use default linear array of receivers on top of model */
-    if (!rec->max_nrec && Nx1==0) Nx1=1; // Default is to use top of model to record data
-
-    if (Nx1) {
-        /* Allocate Start & End Points of Linear Arrays */
-        yrcv1=(float *)malloc(Nx1*sizeof(float));
-        yrcv2=(float *)malloc(Nx1*sizeof(float));
-        xrcv1=(float *)malloc(Nx1*sizeof(float));
-        xrcv2=(float *)malloc(Nx1*sizeof(float));
-        zrcv1=(float *)malloc(Nx1*sizeof(float));
-        zrcv2=(float *)malloc(Nx1*sizeof(float));
-        if (!getparfloat("yrcv1",yrcv1)) yrcv1[0]=sub_y0;
-        if (!getparfloat("yrcv2",yrcv2)) yrcv2[0]=sub_y1;
-        if (!getparfloat("xrcv1",xrcv1)) xrcv1[0]=sub_x0;
-        if (!getparfloat("xrcv2",xrcv2)) xrcv2[0]=sub_x1;
-        if (!getparfloat("zrcv1",zrcv1)) zrcv1[0]=sub_z0;
-        if (!getparfloat("zrcv2",zrcv2)) zrcv2[0]=zrcv1[0];
-
-		/* check if receiver arrays fit into model */
-		for (iarray=0; iarray<Nx1; iarray++) {
-			yrcv1[iarray] = MAX(sub_y0,      yrcv1[iarray]);
-			yrcv1[iarray] = MIN(sub_y0+ny*dy,yrcv1[iarray]);
-			yrcv2[iarray] = MAX(sub_y0,      yrcv2[iarray]);
-			yrcv2[iarray] = MIN(sub_y0+ny*dy,yrcv2[iarray]);
-			
-			xrcv1[iarray] = MAX(sub_x0,      xrcv1[iarray]);
-			xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
-			xrcv2[iarray] = MAX(sub_x0,      xrcv2[iarray]);
-			xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
-			
-			zrcv1[iarray] = MAX(sub_z0,      zrcv1[iarray]);
-			zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
-			zrcv2[iarray] = MAX(sub_z0,      zrcv2[iarray]);
-			zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
-		}
-
-        /* Crop to Fit Model */
-/* Max's addtion still have to check if it has the same fucntionality */
-        for (iarray=0;iarray<Nx1;iarray++) {
-            if (yrcv1[iarray]<sub_y0) {
-                if (yrcv2[iarray]<sub_y0) {
-                    verr("Linear array %li outside model bounds",iarray);
-                }
-				else {
-                    vwarn("Cropping element %li of 'yrcv1' (%f) to model bounds (%f)",iarray,yrcv1[iarray],sub_y0);
-                    yrcv1[iarray]=sub_y0;
-                }
-            } 
-			else if (yrcv1[iarray] > sub_y1) {
-                verr("Linear array %li outside model bounds",iarray);
-            }
-            if ( (yrcv2[iarray] < yrcv1[iarray]) ) {
-                verr("Ill defined linear array %li, 'yrcv1' (%f) greater than 'yrcv2' (%f)",iarray,yrcv1[iarray],yrcv2[iarray]);
-            }
-			else if (yrcv2[iarray]>sub_y1) {
-                vwarn("Cropping element %li of 'yrcv2' (%f) to model bounds (%f)",iarray,yrcv2[iarray],sub_y1);
-                yrcv2[iarray]=sub_y1;
-            }
-
-            if (xrcv1[iarray]<sub_x0) {
-                if (xrcv2[iarray]<sub_x0) {
-                    verr("Linear array %li outside model bounds",iarray);
-                }
-				else {
-                    vwarn("Cropping element %li of 'xrcv1' (%f) to model bounds (%f)",iarray,xrcv1[iarray],sub_x0);
-                    xrcv1[iarray]=sub_x0;
-                }
-            } 
-			else if (xrcv1[iarray] > sub_x1) {
-                verr("Linear array %li outside model bounds",iarray);
-            }
-            if ( (xrcv2[iarray] < xrcv1[iarray]) ) {
-                verr("Ill defined linear array %li, 'xrcv1' (%f) greater than 'xrcv2' (%f)",iarray,xrcv1[iarray],xrcv2[iarray]);
-            }
-			else if (xrcv2[iarray]>sub_x1) {
-                vwarn("Cropping element %li of 'xrcv2' (%f) to model bounds (%f)",iarray,xrcv2[iarray],sub_x1);
-                xrcv2[iarray]=sub_x1;
-            }
-
-            if (zrcv1[iarray] < sub_z0) {
-                if (zrcv2[iarray] < sub_z0) {
-                    verr("Linear array %li outside model bounds",iarray);
-                }
-				else {
-               		vwarn("Cropping element %li of 'zrcv1' (%f) to model bounds (%f)",iarray,zrcv1[iarray],sub_z0);
-                	zrcv1[iarray]=sub_z0;
-                }
-            }
-			else if (zrcv1[iarray] > sub_z1) {
-                verr("Linear array %li outside model bounds",iarray);
-            }
-            if ( (zrcv2[iarray] < zrcv1[iarray]) ) {
-                verr("Ill defined linear array %li, 'zrcv1' (%f) greater than 'zrcv2' (%f)",iarray,zrcv1[iarray],zrcv2[iarray]);
-            }
-			else if (zrcv2[iarray]>sub_z1) {
-                vwarn("Cropping element %li of 'xrcv2' (%f) to model bounds (%f)",iarray,zrcv2[iarray],sub_z1);
-                zrcv2[iarray]=sub_z1;
-            }
-        }
-
-        /* Get Sampling Rates */
-		Ndy = countparval("dyrcv");
-		Ndx = countparval("dxrcv");
-		Ndz = countparval("dzrcv");
-
-		dyr = (float *)malloc(Nx1*sizeof(float));
-		dxr = (float *)malloc(Nx1*sizeof(float));
-		dzr = (float *)malloc(Nx1*sizeof(float));
-		if(!getparfloat("dyrcv", dyr)) dyr[0]=dy;
-		if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
-		if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
-/* TODO logic extended to 3D */
-		if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
-			for (i=1; i<Nx1; i++) {
-				dxr[i] = dxr[0];
-				dzr[i] = dzr[0];
-			}
-			Ndx=1;
-			Ndz=1;
-		}
-		else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
-			for (i=1; i<Nx1; i++) {
-				dxr[i] = dxr[0];
-				dzr[i] = dzr[0];
-			}
-			Ndz=1;
-			Ndx=1;
-		}
-/* end TODO */
-		else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
-			if (Ndx!=Ndz) {
-				verr("Number of 'dxrcv' (%li) is not equal to number of 'dzrcv' (%li) or 1",Ndx,Ndz);
-			}
-			if (Ndx!=Nx1 && Ndx!=1) {
-				verr("Number of 'dxrcv' (%li) is not equal to number of starting points in 'xrcv1' (%li) or 1",Ndx,Nx1);
-			}
-			if (Ndy!=Ndz) {
-				verr("Number of 'dyrcv' (%li) is not equal to number of 'dzrcv' (%li) or 1",Ndy,Ndz);
-			}
-		}
-
-		/* check consistency of receiver steps */
-        for (iarray=0; iarray<Ndy; iarray++) {
-            if (dyr[iarray]<0) {
-				dyr[i]=dy;
-				vwarn("'dyrcv' element %li (%f) is less than zero, changing it to %f'",iarray,dyr[iarray],dy);
-			}
-        }
-        for (iarray=0; iarray<Ndx; iarray++) {
-            if (dxr[iarray]<0) {
-				dxr[i]=dx;
-				vwarn("'dxrcv' element %li (%f) is less than zero, changing it to %f'",iarray,dxr[iarray],dx);
-			}
-        }
-        for (iarray=0;iarray<Ndz;iarray++) {
-            if (dzr[iarray]<0) {
-				dzr[iarray]=dz;
-				vwarn("'dzrcv' element %li (%f) is less than zero, changing it to %f'",iarray,dzr[iarray],dz);
-			}
-        }
-        for (iarray=0;iarray<Ndx;iarray++){
-            if (dxr[iarray]==0 && dzr[iarray]==0) {
-                xrcv2[iarray]=xrcv1[iarray];
-				dxr[iarray]=1.;
-                vwarn("'dxrcv' element %li & 'dzrcv' element 1 are both 0.",iarray+1);
-                vmess("Placing 1 receiver at (%li,%li)",xrcv1[iarray],zrcv1[iarray]);
-            }
-        }
-        for (iarray=0;iarray<Ndx;iarray++){
-            if (xrcv1[iarray]==xrcv2[iarray] && dxr[iarray]!=0) {
-                dxr[iarray]=0.;
-                vwarn("Linear array %li: 'xrcv1'='xrcv2' and 'dxrcv' is not 0. Setting 'dxrcv'=0",iarray+1);
-            }
-        }
-        for (iarray=0;iarray<Ndx;iarray++){
-            if (zrcv1[iarray]==zrcv2[iarray] && dzr[iarray]!=0.){
-                dzr[iarray]=0.;
-                vwarn("Linear array %li: 'zrcv1'='zrcv2' and 'dzrcv' is not 0. Setting 'dzrcv'=0",iarray+1);
-            }
-        }
-
-        /* Calculate Number of Receivers */
-		nrcv = 0;
-        nlrcv=(long *)malloc(Nx1*sizeof(long));
-		for (iarray=0; iarray<Nx1; iarray++) {
-			yrange = (yrcv2[iarray]-yrcv1[iarray]); 
-			xrange = (xrcv2[iarray]-xrcv1[iarray]); 
-			zrange = (zrcv2[iarray]-zrcv1[iarray]); 
-			if (dyr[iarray] != 0.0) {
-				nlrcv[iarray] = NINT(fabs(yrange/dyr[iarray]))+1;
-			}
-			if (dxr[iarray] != 0.0) {
-				nlrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
-			}
-			else {
-				if (dzr[iarray] == 0) {
-					verr("For receiver array %li: receiver distance dzrcv is not given", iarray);
-				}
-				nlrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
-			}
-            nrcv+=nlrcv[iarray];
-		}
-
-        /* Calculate Number of Receivers */
-/*
-        nlrcv=(long *)malloc(Nx1*sizeof(long));
-        if (!isnan(*xrcv1)) *nlrcv=MIN(NINT((*xrcv2-*xrcv1)/(*dxr)),NINT((*zrcv2-*zrcv1)/(*dzr)))+1;
-        else *nlrcv=0;
-        nrcv=*nlrcv;
-        if (verbose>4 && nlrcv[iarray]!=0) vmess("Linear receiver array 1 has final bounds: (X: %f -> %f,Z: %f ->
-%f)",xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*(*dxr),zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*(*dzr));
-        if (Ndx>1) {
-            for (iarray=1;iarray<Nx1;iarray++) {
-                if (!isnan(xrcv1[iarray])) {
-					nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/dxr[iarray]),NINT((zrcv2[iarray]-zrcv1[iarray])/dzr[iarray]))+1;
-				}
-                else {
-					nlrcv[iarray]=0;
-				}
-                nrcv+=nlrcv[iarray];
-                if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %li has final bounds: (X: %f -> %f,Z: %f ->
-%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*dxr[iarray],zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*dzr[iarray]);
-            }
-        }
-		 else {
-            for (iarray=1;iarray<Nx1;iarray++) {
-                if (!isnan(xrcv1[iarray])) nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/(*dxr)),NINT((zrcv2[iarray]-zrcv1[iarray])/(*dzr)))+1;
-                else nlrcv[iarray]=0;
-                nrcv+=nlrcv[iarray];
-                if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %li has final bounds: (X: %f -> %f,Z: %f ->
-%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]**dxr,zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]**dzr);
-            }
-        }
-*/
-        if (verbose) vmess("Total number of linear array receivers: %li",nrcv);
-        if (!nrcv) {
-            free(yrcv1);
-            free(yrcv2);
-            free(xrcv1);
-            free(xrcv2);
-            free(zrcv1);
-            free(zrcv2);
-            free(dyr);
-            free(dxr);
-            free(dzr);
-            free(nlrcv);
-        }
-        rec->max_nrec+=nrcv;
-    } 
-	else {
-		nrcv=0;
-	}
-
-/* allocate the receiver arrays */
-
-    /* Total Number of Receivers */
-    if (verbose) vmess("Total number of receivers: %li",rec->max_nrec);
-
-    /* Allocate Arrays */
-    rec->y  = (long *)calloc(rec->max_nrec,sizeof(long));
-    rec->x  = (long *)calloc(rec->max_nrec,sizeof(long));
-    rec->z  = (long *)calloc(rec->max_nrec,sizeof(long));
-    rec->yr = (float *)calloc(rec->max_nrec,sizeof(float));
-    rec->xr = (float *)calloc(rec->max_nrec,sizeof(float));
-    rec->zr = (float *)calloc(rec->max_nrec,sizeof(float));
-
-/* read in the receiver postions */
-
-	nrec=0;
-    /* Receivers on a Circle */
-    if (ncrcv) {
-		if (!getparfloat("oyrcv",&oxrcv)) oyrcv=0.0;
-		if (!getparfloat("oxrcv",&oxrcv)) oxrcv=0.0;
-		if (!getparfloat("ozrcv",&ozrcv)) ozrcv=0.0;
-		if (!getparfloat("arcv",&arcv)) {
-			arcv=rrcv; 
-			for (ix=0; ix<ncrcv; ix++) {
-				rec->yr[ix] = oyrcv-sub_y0+rrcv*cos(((iy*dphi)/360.0)*(2.0*M_PI));
-				rec->xr[ix] = oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI));
-				rec->zr[ix] = ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI));
-				rec->y[ix] = NINT(rec->yr[ix]/dy);
-				rec->x[ix] = NINT(rec->xr[ix]/dx);
-				rec->z[ix] = NINT(rec->zr[ix]/dz);
-				//rec->x[ix] = NINT((oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI)))/dx);
-				//rec->z[ix] = NINT((ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI)))/dz);
-				if (verbose>4) fprintf(stderr,"Receiver Circle: yrcv[%li]=%f xrcv=%f zrcv=%f\n", ix, rec->yr[ix]+sub_y0, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
-			}
-		}
-		else { /* an ellipse */
-/*TODO fix for 3D */
-			/* simple numerical solution to find equidistant points on an ellipse */
-			nh  = (ncrcv)*1000; /* should be fine enough for most configurations */
-			h = 2.0*M_PI/nh;
-			a = MAX(rrcv, arcv);
-			b = MIN(rrcv, arcv);
-			e = sqrt(a*a-b*b)/a;
-			//fprintf(stderr,"a=%f b=%f e=%f\n", a, b, e);
-			circ = 0.0;
-			for (ir=0; ir<nh; ir++) {
-				s = sin(ir*h);
-				circ += sqrt(1.0-e*e*s*s);
-			}
-			circ = a*h*circ;
-			//fprintf(stderr,"circ = %f circle=%f\n", circ, 2.0*M_PI*rrcv);
-			/* define distance between receivers on ellipse */
-			dr = circ/ncrcv;
-			ix = 0;
-			srun = 0.0;
-			if (arcv >= rrcv) phase=0.0;
-			else phase=0.5*M_PI;
-			for (ir=0; ir<nh; ir++) {
-				s = sin(ir*h);
-				srun += sqrt(1.0-e*e*s*s);
-				if (a*h*srun >= ix*dr ) {
-					yr = rrcv*cos(ir*h+phase);
-					xr = rrcv*cos(ir*h+phase);
-					zr = arcv*sin(ir*h+phase);
-					rec->yr[ix] = oyrcv-sub_y0+xr;
-					rec->xr[ix] = oxrcv-sub_x0+xr;
-					rec->zr[ix] = ozrcv-sub_z0+zr;
-					rec->y[ix] = NINT(rec->yr[ix]/dy);
-					rec->x[ix] = NINT(rec->xr[ix]/dx);
-					rec->z[ix] = NINT(rec->zr[ix]/dz);
-					if (verbose>4) fprintf(stderr,"Receiver Ellipse: xrcv[%li]=%f zrcv=%f\n", ix, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
-					ix++;
-				}
-				if (ix == ncrcv) break;
-			}
-		}
-
-		/* check if receivers fit into the model otherwise clip to edges */
-		for (ix=0; ix<ncrcv; ix++) {
-			rec->y[ix] = MIN(ny-1, MAX(rec->y[ix], 0));
-			rec->x[ix] = MIN(nx-1, MAX(rec->x[ix], 0));
-			rec->z[ix] = MIN(nz-1, MAX(rec->z[ix], 0));
-		}
-		nrec += ncrcv;
-	}
-
-    /* Receiver Text File */
-
-    if (ntrcv) {
-		/* Allocate arrays */
-		yrcva = (float *)malloc(nrcv*sizeof(float));
-		xrcva = (float *)malloc(nrcv*sizeof(float));
-		zrcva = (float *)malloc(nrcv*sizeof(float));
-		/* Read in receiver coordinates */
-		for (i=0;i<nrcv;i++) {
-			if (fscanf(fp,"%e %e %e\n",&yrcva[i],&xrcva[i],&zrcva[i])!=3) vmess("Receiver Text File: Can not parse coordinates on line %li.",i);
-		}
-		/* Close file */
-		fclose(fp);
-		/* Process coordinates */
-		for (ix=0; ix<nrcv; ix++) {
-			rec->yr[nrec+ix] = yrcva[ix]-sub_y0;
-			rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
-			rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
-			rec->y[nrec+ix] = NINT((yrcva[ix]-sub_y0)/dy);
-			rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
-			rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
-			if (verbose>4) fprintf(stderr,"Receiver Text Array: yrcv[%li]=%f xrcv=%f zrcv=%f\n", ix, rec->yr[ix]+sub_y0, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
-		}
-		free(yrcva);
-		free(xrcva);
-		free(zrcva);
-		nrec += ntrcv;
-	}
-
-    /* Receiver Array */
-	if (nxrcv != 0) {
-		/* receiver array is defined */
-		yrcva = (float *)malloc(nxrcv*sizeof(float));
-		xrcva = (float *)malloc(nxrcv*sizeof(float));
-		zrcva = (float *)malloc(nxrcv*sizeof(float));
-		getparfloat("yrcva", yrcva);
-		getparfloat("xrcva", xrcva);
-		getparfloat("zrcva", zrcva);
-		for (ix=0; ix<nxrcv; ix++) {
-			rec->yr[nrec+ix] = yrcva[ix]-sub_y0;
-			rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
-			rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
-			rec->y[nrec+ix] = NINT((yrcva[ix]-sub_y0)/dy);
-			rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
-			rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
-			if (verbose>4) fprintf(stderr,"Receiver Array: yrcv[%li]=%f xrcv=%f zrcv=%f\n", ix, rec->yr[ix]+sub_y0, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
-		}
-		free(yrcva);
-		free(xrcva);
-		free(zrcva);
-		nrec += nxrcv;
-	}
-
-    /* Linear Receiver Arrays */
-    if (nrcv!=0) {
-		yrcv1 = (float *)malloc(Nx1*sizeof(float));
-		yrcv2 = (float *)malloc(Nx1*sizeof(float));
-		xrcv1 = (float *)malloc(Nx1*sizeof(float));
-		xrcv2 = (float *)malloc(Nx1*sizeof(float));
-		zrcv1 = (float *)malloc(Nx1*sizeof(float));
-		zrcv2 = (float *)malloc(Nx1*sizeof(float));
-		
-		if(!getparfloat("yrcv1", yrcv1)) yrcv1[0]=sub_y0;
-		if(!getparfloat("yrcv2", yrcv2)) yrcv2[0]=(ny-1)*dy+sub_y0;
-		if(!getparfloat("xrcv1", xrcv1)) xrcv1[0]=sub_x0;
-		if(!getparfloat("xrcv2", xrcv2)) xrcv2[0]=(nx-1)*dx+sub_x0;
-		if(!getparfloat("zrcv1", zrcv1)) zrcv1[0]=sub_z0;
-		if(!getparfloat("zrcv2", zrcv2)) zrcv2[0]=zrcv1[0];		
-		
-		Ndy = countparval("dyrcv");
-		Ndx = countparval("dxrcv");
-		Ndz = countparval("dzrcv");
-
-		dyr = (float *)malloc(Nx1*sizeof(float));
-		dxr = (float *)malloc(Nx1*sizeof(float));
-		dzr = (float *)malloc(Nx1*sizeof(float));
-		if(!getparfloat("dyrcv", dyr)) dyr[0]=dy;
-		if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
-		if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
-		if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
-			for (i=1; i<Nx1; i++) {
-				dyr[i] = dyr[0];
-				dxr[i] = dxr[0];
-				dzr[i] = dzr[0];
-			}
-			Ndx=1;
-		}
-		else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
-			for (i=1; i<Nx1; i++) {
-				dyr[i] = dyr[0];
-				dxr[i] = dxr[0];
-				dzr[i] = dzr[0];
-			}
-			Ndz=1;
-		}
-		else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
-			if (Ndx>1) assert(Ndx==Nx1);
-			if (Ndz>1) assert(Ndz==Nx1);
-			if (Ndy>1) assert(Ndy==Nx1);
-		}
-		
-/*
-		if ( (Ndx!=0) && (Ndz!=0) ) {
-			vwarn("Both dzrcv and dxrcv are set: dxrcv value is used");
-			Ndz=0;
-			for (i=0; i<Nx1; i++) dzr[i] = 0.0;
-		}
-*/
-		/* check if receiver arrays fit into model */
-		for (iarray=0; iarray<Nx1; iarray++) {
-			yrcv1[iarray] = MAX(sub_y0,      yrcv1[iarray]);
-			yrcv1[iarray] = MIN(sub_y0+ny*dy,yrcv1[iarray]);
-			yrcv2[iarray] = MAX(sub_y0,      yrcv2[iarray]);
-			yrcv2[iarray] = MIN(sub_y0+ny*dy,yrcv2[iarray]);
-
-			xrcv1[iarray] = MAX(sub_x0,      xrcv1[iarray]);
-			xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
-			xrcv2[iarray] = MAX(sub_x0,      xrcv2[iarray]);
-			xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
-			
-			zrcv1[iarray] = MAX(sub_z0,      zrcv1[iarray]);
-			zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
-			zrcv2[iarray] = MAX(sub_z0,      zrcv2[iarray]);
-			zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
-		}
-
-		/* calculate receiver array and store into rec->x,z */
-
-		for (iarray=0; iarray<Nx1; iarray++) {
-			yrange = (yrcv2[iarray]-yrcv1[iarray]); 
-			xrange = (xrcv2[iarray]-xrcv1[iarray]); 
-			zrange = (zrcv2[iarray]-zrcv1[iarray]); 
-			if (dxr[iarray] != 0.0) {
-				nrcv = nlrcv[iarray];
-				dxrcv=dxr[iarray];
-				dzrcv = zrange/(nrcv-1);
-				if (dzrcv != dzr[iarray]) {
-					vwarn("For receiver array %li: calculated dzrcv=%f given=%f", iarray, dzrcv, dzr[iarray]);
-					vwarn("The calculated receiver distance %f is used", dzrcv);
-				}
-			}
-			else {
-				if (dzr[iarray] == 0) {
-					verr("For receiver array %li: receiver distance dzrcv is not given", iarray);
-				}
-				nrcv = nlrcv[iarray];
-				dxrcv = xrange/(nrcv-1);
-				dzrcv = dzr[iarray];
-				if (dxrcv != dxr[iarray]) {
-					vwarn("For receiver array %li: calculated dxrcv=%f given=%f", iarray, dxrcv, dxr[iarray]);
-					vwarn("The calculated receiver distance %f is used", dxrcv);
-				}
-			}
-
-			// calculate coordinates
-			for (ir=0; ir<nrcv; ir++) {
-				rec->yr[nrec]=yrcv1[iarray]-sub_y0+ir*dyrcv;
-				rec->xr[nrec]=xrcv1[iarray]-sub_x0+ir*dxrcv;
-				rec->zr[nrec]=zrcv1[iarray]-sub_z0+ir*dzrcv;
-
-				rec->y[nrec]=NINT((rec->yr[nrec])/dy);
-				rec->x[nrec]=NINT((rec->xr[nrec])/dx);
-				rec->z[nrec]=NINT((rec->zr[nrec])/dz);
-				nrec++;
-			}
-		}
-		free(yrcv1);
-		free(yrcv2);
-		free(xrcv1);
-		free(xrcv2);
-		free(zrcv1);
-		free(zrcv2);
-		free(dyr);
-		free(dxr);
-		free(dzr);
-        free(nlrcv);
-	}
-
-    rec->n=rec->max_nrec;
-	return 0;
-}
diff --git a/raytime3d/src3d.c b/raytime3d/src3D.c
similarity index 72%
rename from raytime3d/src3d.c
rename to raytime3d/src3D.c
index 80a2d2d80c98e1c75c75dccef816b56911c3fc09..369cd993233172b3a0adb14fc2bd525394eef535 100644
--- a/raytime3d/src3d.c
+++ b/raytime3d/src3D.c
@@ -17,7 +17,7 @@ extern void verr(char *fmt, ...);
 extern void vwarn(char *fmt, ...);
 extern void vmess(char *fmt, ...);
 
-void src3d(float *time0, float *slow0, long nz, long nx, long ny, float h, float ox, float oy, float oz, long *pxs, long *pys, long *pzs, long *cube)
+void src3D(float *time0, float *slow0, long nz, long nx, long ny, float h, float fxs, float fys, float fzs, long xs, long ys, long zs, long *cube)
 {
 	long
 		srctype=1,	/* if 1, source is a point;
@@ -26,18 +26,18 @@ void src3d(float *time0, float *slow0, long nz, long nx, long ny, float h, float
 		srcwall,	/* if 1, source on x=0 wall, if 2, on x=nx-1 wall
 						if 3, source on y=0 wall, if 4, on y=ny-1 wall
 						if 5, source on z=0 wall, if 6, on z=nz-1 wall */
-		xs,			/* shot x position (in grid points) */
-		ys,			/* shot y position */
-		zs,			/* shot depth */
+		xs1,			/* shot x position (in grid points) */
+		ys1,			/* shot y position */
+		zs1,			/* shot depth */
 		xx, yy, zz,	/* Used to loop around xs, ys, zs coordinates	*/
 		ii, i, j, k, 
 		wfint, ofint,
 		nxy, nyz, nxz, nxyz, nwall,
 		NCUBE=2;
 	float
-		fxs,	/* shot position in X (in real units)*/
-		fys,	/* shot position in Y (in real units)*/
-		fzs,	/* shot position in Z (in real units)*/
+		fxs1,	/* shot position in X (in real units)*/
+		fys1,	/* shot position in Y (in real units)*/
+		fzs1,	/* shot position in Z (in real units)*/
 		*wall,
 		/* maximum offset (real units) to compute */
 		/* used in linear velocity gradient cube source */
@@ -50,51 +50,15 @@ void src3d(float *time0, float *slow0, long nz, long nx, long ny, float h, float
 
 
 	if(!getparlong("NCUBE",&NCUBE)) NCUBE=2;
+	*cube = NCUBE;
 
 	if(!getparlong("srctype",&srctype)) srctype=1;
-	if(srctype==1) {
-		if(!getparfloat("xsrc1",&xsrc1)) verr("xsrc1 not given");
-		if(!getparfloat("ysrc1",&ysrc1)) verr("ysrc1 not given");
-		if(!getparfloat("zsrc1",&zsrc1)) verr("zsrc1 not given");
-		fxs = (xsrc1-ox)/h;
-		fys = (ysrc1-oy)/h;
-		fzs = (zsrc1-oz)/h;
-		xs = (long)(fxs + 0.5);
-		ys = (long)(fys + 0.5);
-		zs = (long)(fzs + 0.5);
-		if(xs<2 || ys<2 || zs<2 || xs>nx-3 || ys>ny-3 || zs>nz-3){
-			vwarn("Source near an edge, beware of traveltime errors");
-			vwarn("for raypaths that travel parallel to edge ");
-			vwarn("while wavefronts are strongly curved, (JV, 8/17/88)\n");
-		}
-		*pxs = xs; *pys = ys, *pzs = zs, *cube = NCUBE;
-	}
-	else if (srctype==2)  {
-		if (!getparlong("srcwall",&srcwall)) verr("srcwall not given");
-		if (!getparstring("wallfile",&wallfile)) verr("wallfile not given");
-		if((wfint=open(wallfile,O_RDONLY,0664))<=1) {
-			fprintf(stderr,"cannot open %s\n",wallfile);
-			exit(-1);
-		}
-	}
-	else if (srctype==3)  {
-		if (!getparlong("srcwall",&srcwall)) verr("srcwall not given");
-		if (!getparstring("oldtfile",&oldtfile)) verr("oldtfile not given");
-		if((ofint=open(oldtfile,O_RDONLY,0664))<=1) {
-			fprintf(stderr,"cannot open %s\n",oldtfile);
-			exit(-1);
-		}
-	}
-	else  {
-		verr("ERROR: incorrect value of srctype");
-	}
 
 	nxy = nx * ny;
 	nyz = ny * nz;
 	nxz = nx * nz;
 	nxyz = nx * ny * nz;
 
-
 	/* SET TIMES TO DUMMY VALUE */
 	for(i=0;i<nxyz;i++) time0[i] = 1.0e10;
 
diff --git a/raytime3d/vidale3d.c b/raytime3d/vidale3d.c
index 3faaefa7779b72bded75283b63d984f1afc01715..d7e560baa6810a5c645388facb1b4a2b0724c5b5 100644
--- a/raytime3d/vidale3d.c
+++ b/raytime3d/vidale3d.c
@@ -348,7 +348,7 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
 				if (fhead>headtest)  headw[5]++;
 			}
 		}
-		if(z1 == 0) dz1 = 0;
+		if(z1 == 1) dz1 = 0;
 		z1--;
 	}
       }
@@ -527,7 +527,7 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
 				if (fhead>headtest)  headw[6]++;
 			}
 		}
-		if(z2 == nz-1) dz2 = 0;
+		if(z2 == nz-2) dz2 = 0;
 		z2++;
 	}
       }
@@ -706,7 +706,7 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
 				if (fhead>headtest)  headw[3]++;
 			}
 		}
-		if(y1 == 0) dy1 = 0;
+		if(y1 == 1) dy1 = 0;
 		y1--;
 	}
       }
@@ -885,7 +885,7 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
 				if (fhead>headtest)  headw[4]++;
 			}
 		}
-		if(y2 == ny-1) dy2 = 0;
+		if(y2 == ny-2) dy2 = 0;
 		y2++;
 	}
       }
@@ -1064,7 +1064,7 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
 				if (fhead>headtest)  headw[1]++;
 			}
 		}
-		if(x1 == 0) dx1 = 0;
+		if(x1 == 1) dx1 = 0;
 		x1--;
 	}
       }
@@ -1243,14 +1243,14 @@ void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, lo
 				if (fhead>headtest)  headw[2]++;
 			}
 		}
-		if(x2 == nx-1) dx2 = 0;
+		if(x2 == nx-2) dx2 = 0;
 		x2++;
 	}
       }
 
 		/* UPDATE RADIUS */
 		radius++;
-		if(radius%10 == 0 && verbose) vmess("Completed radius = %li",radius);
+		if(radius%10 == 0 && verbose>5) vmess("Completed radius = %li",radius);
         if(radius == maxrad) rad0 = 0;
 
 	}	/* END BIG LOOP */
diff --git a/raytime3d/vidale3d_backup.c b/raytime3d/vidale3d_backup.c
deleted file mode 100644
index 89fb16e1813054960f2f113ffdbde24d74d7c316..0000000000000000000000000000000000000000
--- a/raytime3d/vidale3d_backup.c
+++ /dev/null
@@ -1,1453 +0,0 @@
-#include<stdlib.h>
-#include<stdio.h>
-#include<math.h>
-#include<assert.h>
-#include<string.h>
-#include"par.h"
-
-#define SQR2 1.414213562
-#define SQR3 1.732050808
-#define SQR6 2.449489743
-#define t0(x,y,z)   time0[nxz*(y) + nz*(x) + (z)]
-#define s0(x,y,z)   slow0[nxz*(y) + nz*(x) + (z)]
-
-/* definitions from verbose.c */
-extern void verr(char *fmt, ...);
-extern void vwarn(char *fmt, ...);
-extern void vmess(char *fmt, ...);
-
-struct sorted
-	{ float time; long i1, i2;};
-
-// int compar(struct sorted *a,struct sorted *b);
-int compar(const void * a, const void * b);
-
-float fdhne(float t1,float t2,float t3,float t4,float t5,float ss0,float s1,float s2,float s3);
-float fdh3d(float t1,float t2,float t3,float t4,float t5,float t6,float t7,float ss0,float s1,float s2,float s3,float s4,float s5,float s6,float s7);
-float fdh2d(float t1,float t2,float t3,float ss0,float s1,float s2,float s3);
-float fdhnf(float t1,float t2,float t3,float t4,float t5,float ss0,float s1);
-
-void vidale3d(float *slow0, float *time0, long nz, long nx, long ny, float h, long xs, long ys, long zs, long NCUBE)
-{
-	long
-		srctype=1,	/* if 1, source is a point;
-						2, source is on the walls of the data volume;
-						3, source on wall, time field known; */
-		srcwall,	/* if 1, source on x=0 wall, if 2, on x=nx-1 wall
-						if 3, source on y=0 wall, if 4, on y=ny-1 wall
-						if 5, source on z=0 wall, if 6, on z=nz-1 wall */
-		iplus=1,	/* rate of expansion of "cube" in the */
-		iminus=1,	/*    plus/minus x/y/z direction */
-		jplus=1,
-		jminus=1,
-		kplus=1,
-		kminus=1,
-		igrow,		/* counter for "cube" growth */
-		X1, X2, lasti, index, ii, i, j, k, radius, 
-		nxy, nyz, nxz, nxyz, nwall,
-		/* counters for the position of the sides of current cube */
-		x1, x2, y1, y2, z1, z2,
-		/* flags set to 1 until a side has been reached */
-		dx1=1, dx2=1, dy1=1, dy2=1, dz1=1, dz2=1, rad0=1,
-		maxrad,		/* maximum radius to compute */
-		reverse=1,	/* will automatically do up to this number of
-						reverse propagation steps to fix waves that travel 
-						back into expanding cube */
-		headpref=6,	/* if headpref starts > 0, will determine 
-						model wall closest to source and will prefer to start
-						reverse calculations on opposite wall */
-		/* counters for detecting head waves on sides of current cube */
-		head,headw[7], verbose;
-	float
-		*wall,
-		guess, try,
-		/* maximum offset (real units) to compute */
-		maxoff = -1.,
-		/* used to detect head waves:  if headwave operator decreases 
-		   the previously-computed traveltime by at least 
-		   headtest*<~time_across_cell> then the headwave counter is 
-		   triggered */
-		fhead,headtest=1.e-3;
-
-	/* ARRAY TO ORDER SIDE FOR SOLUTION IN THE RIGHT ORDER */
-	struct sorted *sort;
-
-	if(!getparlong("verbose",&verbose)) verbose=0;
-	if(!getparfloat("maxoff",&maxoff)) maxoff = -1.;
-	if(!getparlong("iminus",&iminus)) iminus=1;
-	if(!getparlong("iplus",&iplus)) iplus=1;
-	if(!getparlong("jminus",&jminus)) jminus=1;
-	if(!getparlong("jplus",&jplus)) jplus=1;
-	if(!getparlong("kminus",&kminus)) kminus=1;
-	if(!getparlong("kplus",&kplus)) kplus=1;
-	if(!getparlong("reverse",&reverse)) reverse=0;
-	if(!getparlong("headpref",&headpref)) headpref=6;
-	if(!getparlong("NCUBE",&NCUBE)) NCUBE=2;
-
-	/* SET MAXIMUM RADIUS TO COMPUTE */
-	if (maxoff > 0.) {
-		maxrad = maxoff/h + 1;
-		vwarn("WARNING: Computing only to max radius = %li",maxrad);
-	}
-	else maxrad = 99999999;
-
-	nxy = nx * ny;
-	nyz = ny * nz;
-	nxz = nx * nz;
-	nxyz = nx * ny * nz;
-
-	/* MAKE ARRAY SORT LARGE ENOUGH FOR ANY SIDE */
-	if(nx <= ny && nx <= nz)  {
-		sort = (struct sorted *) malloc(sizeof(struct sorted)*ny*nz);
-		nwall = nyz;
-	}
-	else if(ny <= nx && ny <= nz)  {
-		sort = (struct sorted *) malloc(sizeof(struct sorted)*nx*nz);
-		nwall = nxz;
-	}
-	else  {
-		sort = (struct sorted *) malloc(sizeof(struct sorted)*nx*ny);
-		nwall = nxy;
-	}
-	wall = (float *) malloc(4*nwall);
-	if(sort == NULL || wall == NULL) 
-		verr("error in allocation of arrays sort and wall");
-
-	if(!getparlong("srctype",&srctype)) srctype=1;
-	if(srctype==1) {
-		/* SETS LOCATION OF THE SIDES OF THE CUBE	*/
-		radius = NCUBE;
-		if(xs > NCUBE) x1 = xs - (NCUBE + 1);
-		else{ x1 = -1; dx1 = 0;}
-		if(xs < nx-(NCUBE + 1)) x2 = xs + (NCUBE + 1);
-		else{ x2 = nx; dx2 = 0;}
-		if(ys > NCUBE) y1 = ys - (NCUBE + 1);
-		else{ y1 = -1; dy1 = 0;}
-		if(ys < ny-(NCUBE + 1)) y2 = ys + (NCUBE + 1);
-		else{ y2 = ny; dy2 = 0;}
-		if(zs > NCUBE) z1 = zs - (NCUBE + 1);
-		else{ z1 = -1; dz1 = 0;}
-		if(zs < nz-(NCUBE + 1)) z2 = zs + (NCUBE + 1);
-		else{ z2 = nz; dz2 = 0;}
-	}
-	else {
-		if (!getparlong("srcwall",&srcwall)) verr("srcwall not given");
-		/* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE  */
-		radius = 1;
-		if (srcwall == 1)	x2=1;
-		else	{  x2=nx;	dx2=0;  }
-		if (srcwall == 2)	x1=nx-2;
-		else	{  x1= -1;	dx1=0;  }
-		if (srcwall == 3)	y2=1;
-		else	{  y2=ny;	dy2=0;  }
-		if (srcwall == 4)	y1=ny-2;
-		else	{  y1= -1;	dy1=0;  }
-		if (srcwall == 5)	z2=1;
-		else	{  z2=nz;	dz2=0;  }
-		if (srcwall == 6)	z1=nz-2;
-		else	{  z1= -1;	dz1=0;  }
-	}
-
-
-	if (headpref>0) {	/* HOLE - PREFERRED REVERSE DIRECTION */
-		head = nx*ny*nz;
-		if (nx>5 && x2<=head)   {headpref=2;  head=x2;}
-		if (nx>5 && (nx-1-x1)<=head)   {headpref=1;  head=nx-1-x1;}
-		if (ny>5 && y2<=head)   {headpref=4;  head=y2;}
-		if (ny>5 && (ny-1-y1)<=head)   {headpref=3;  head=ny-1-y1;}
-		if (nz>5 && z2<=head)   {headpref=6;  head=z2;}
-		if (nz>5 && (nz-1-z1)<=head)   {headpref=5;  head=nz-1-z1;}
-	}
-
-	/* BIGGER LOOP - HOLE - ALLOWS AUTOMATIC REVERSE PROPAGATION IF 
-		HEAD WAVES ARE ENCOUNTERED ON FACES OF EXPANDING CUBE, 
-		ALLOWING WAVES TO TRAVEL BACK INTO THE CUBE */
-
-	while ( reverse > -1 )  {
-
-		headw[1]=0; headw[2]=0; headw[3]=0; headw[4]=0;
-		headw[5]=0; headw[6]=0;
-
-	/* BIG LOOP */
-	while(rad0 && (dx1 || dx2 || dy1 || dy2 || dz1 || dz2))  {
-
-		/* CALCULATE ON PRIMARY (time0) GRID */
-
-		/* TOP SIDE */
-      for (igrow=1;igrow<=kminus;igrow++) {  
-	if(dz1){
-		ii = 0;
-		for(j=y1+1; j<=y2-1; j++){
-			for(i=x1+1; i<=x2-1; i++){
-				sort[ii].time = t0(i,j,z1+1);
-				sort[ii].i1 = i;
-				sort[ii].i2 = j;
-				ii++;
-			}
-		}
-		qsort((char *)sort,ii,sizeof(struct sorted),compar);
-		for(i=0;i<ii;i++){
-			X1 = sort[i].i1;
-			X2 = sort[i].i2;
-			index = z1*nxy + X2*nx + X1;
-			lasti = (z1+1)*nxy + X2*nx + X1;
-			fhead = 0.;
-			guess = time0[index];
-                        if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
-			   && time0[index+nx] < 1.e9 && X2<ny-1  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,X2,z1+1),
-				      t0(X1+1,X2,z1+1),t0(X1+1,X2+1,z1+1),t0(X1,X2+1,z1+1),
-				      t0(X1+1,X2,z1  ),t0(X1+1,X2+1,z1  ),t0(X1,X2+1,z1  ),
-				      s0(X1,X2,z1), s0(X1,X2,z1+1),
-				      s0(X1+1,X2,z1+1),s0(X1+1,X2+1,z1+1),s0(X1,X2+1,z1+1),
-				      s0(X1+1,X2,z1  ),s0(X1+1,X2+1,z1  ),s0(X1,X2+1,z1  ));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
-			   && time0[index+nx] < 1.e9 && X2<ny-1  && X1>0 ) {
-			  try = fdh3d(              t0(X1,X2,z1+1),
-				      t0(X1-1,X2,z1+1),t0(X1-1,X2+1,z1+1),t0(X1,X2+1,z1+1),
-				      t0(X1-1,X2,z1  ),t0(X1-1,X2+1,z1  ),t0(X1,X2+1,z1  ),
-				      s0(X1,X2,z1), s0(X1,X2,z1+1),
-				      s0(X1-1,X2,z1+1),s0(X1-1,X2+1,z1+1),s0(X1,X2+1,z1+1),
-				      s0(X1-1,X2,z1  ),s0(X1-1,X2+1,z1  ),s0(X1,X2+1,z1  ));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
-			   && time0[index-nx] < 1.e9 && X2>0  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,X2,z1+1),
-				      t0(X1+1,X2,z1+1),t0(X1+1,X2-1,z1+1),t0(X1,X2-1,z1+1),
-				      t0(X1+1,X2,z1  ),t0(X1+1,X2-1,z1  ),t0(X1,X2-1,z1  ),
-				      s0(X1,X2,z1), s0(X1,X2,z1+1),
-				      s0(X1+1,X2,z1+1),s0(X1+1,X2-1,z1+1),s0(X1,X2-1,z1+1),
-				      s0(X1+1,X2,z1  ),s0(X1+1,X2-1,z1  ),s0(X1,X2-1,z1  ));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
-			   && time0[index-nx] < 1.e9 && X2>0  && X1>0 ) {
-			  try = fdh3d(              t0(X1,X2,z1+1),
-				      t0(X1-1,X2,z1+1),t0(X1-1,X2-1,z1+1),t0(X1,X2-1,z1+1),
-				      t0(X1-1,X2,z1  ),t0(X1-1,X2-1,z1  ),t0(X1,X2-1,z1  ),
-				      s0(X1,X2,z1), s0(X1,X2,z1+1),
-				      s0(X1-1,X2,z1+1),s0(X1-1,X2-1,z1+1),s0(X1,X2-1,z1+1),
-				      s0(X1-1,X2,z1  ),s0(X1-1,X2-1,z1  ),s0(X1,X2-1,z1  ));
-			  if (try<guess) guess = try;
-			}
-			if(guess > 1.0e9){ 
-			  if(time0[index+1] < 1.e9 && X1<nx-1 && X2>y1+1 && X2<y2-1 )  {
-			      try = fdhne(t0(X1,X2,z1+1),t0(X1+1,X2,z1+1),t0(X1+1,X2,z1),
-					  t0(X1+1,X2-1,z1+1),t0(X1+1,X2+1,z1+1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1+1,X2,z1+1),s0(X1+1,X2,z1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 && X2>y1+1 && X2<y2-1 )  {
-			      try = fdhne(t0(X1,X2,z1+1),t0(X1-1,X2,z1+1),t0(X1-1,X2,z1),
-					  t0(X1-1,X2-1,z1+1),t0(X1-1,X2+1,z1+1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1-1,X2,z1+1),s0(X1-1,X2,z1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nx] < 1.e9 && X2<ny-1 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,X2,z1+1),t0(X1,X2+1,z1+1),t0(X1,X2+1,z1),
-					  t0(X1-1,X2+1,z1+1),t0(X1+1,X2+1,z1+1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1,X2+1,z1+1),s0(X1,X2+1,z1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,X2,z1+1),t0(X1,X2-1,z1+1),t0(X1,X2-1,z1),
-					  t0(X1-1,X2-1,z1+1),t0(X1+1,X2-1,z1+1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1,X2-1,z1+1),s0(X1,X2-1,z1) );
-			    if (try<guess)  guess = try;
-			  }
-		        } 
-			  if(time0[index+1] < 1.e9 && X1<nx-1 )  {
-			    try = fdh2d(t0(X1,X2,z1+1),t0(X1+1,X2,z1+1),t0(X1+1,X2,z1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1+1,X2,z1+1),s0(X1+1,X2,z1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 )  {
-			    try = fdh2d(t0(X1,X2,z1+1),t0(X1-1,X2,z1+1),t0(X1-1,X2,z1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1-1,X2,z1+1),s0(X1-1,X2,z1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nx] < 1.e9 && X2<ny-1 )  {
-			    try = fdh2d(t0(X1,X2,z1+1),t0(X1,X2+1,z1+1),t0(X1,X2+1,z1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1,X2+1,z1+1),s0(X1,X2+1,z1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X2>0 )  {
-			    try = fdh2d(t0(X1,X2,z1+1),t0(X1,X2-1,z1+1),t0(X1,X2-1,z1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1),s0(X1,X2-1,z1+1),s0(X1,X2-1,z1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
-			     && time0[index+nx] < 1.e9 && X2<ny-1  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,X2,z1),t0(X1+1,X2+1,z1),t0(X1,X2+1,z1),
-					s0(X1,X2,z1),
-					s0(X1+1,X2,z1),s0(X1+1,X2+1,z1),s0(X1,X2+1,z1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
-			     && time0[index-nx] < 1.e9 && X2>0  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,X2,z1),t0(X1+1,X2-1,z1),t0(X1,X2-1,z1),
-					s0(X1,X2,z1),
-					s0(X1+1,X2,z1),s0(X1+1,X2-1,z1),s0(X1,X2-1,z1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
-			     && time0[index+nx] < 1.e9 && X2<ny-1  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,X2,z1),t0(X1-1,X2+1,z1),t0(X1,X2+1,z1),
-					s0(X1,X2,z1),
-					s0(X1-1,X2,z1),s0(X1-1,X2+1,z1),s0(X1,X2+1,z1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
-			     && time0[index-nx] < 1.e9 && X2>0  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,X2,z1),t0(X1-1,X2-1,z1),t0(X1,X2-1,z1),
-					s0(X1,X2,z1),
-					s0(X1-1,X2,z1),s0(X1-1,X2-1,z1),s0(X1,X2-1,z1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if(guess > 1.0e9){ 
-			  if ( X1>x1+1 && X1<x2-1 && X2>y1+1 && X2<y2-1 ) {
-			    try = fdhnf(t0(X1,X2,z1+1),
-					  t0(X1+1,X2,z1+1),t0(X1,X2+1,z1+1),
-					  t0(X1-1,X2,z1+1),t0(X1,X2-1,z1+1),
-					  s0(X1,X2,z1),
-					  s0(X1,X2,z1+1) );
-			    if (try<guess)  guess = try;
-			  }
-			} 
-                          try = t0(X1,X2,z1+1) + .5*(s0(X1,X2,z1)+s0(X1,X2,z1+1));
-			  if (try<guess)  guess = try;
-                          if ( time0[index+1]<1.e9 && X1<nx-1 )  {
-			    try = t0(X1+1,X2,z1) + .5*(s0(X1,X2,z1)+s0(X1+1,X2,z1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-1]<1.e9 && X1>0 )  {
-			    try = t0(X1-1,X2,z1) + .5*(s0(X1,X2,z1)+s0(X1-1,X2,z1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index+nx]<1.e9 && X2<ny-1 )  {
-			    try = t0(X1,X2+1,z1) + .5*(s0(X1,X2,z1)+s0(X1,X2+1,z1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nx]<1.e9 && X2>0 )  {
-			    try = t0(X1,X2-1,z1) + .5*(s0(X1,X2,z1)+s0(X1,X2-1,z1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if (guess<time0[index])  {
-				time0[index] = guess;
-				if (fhead>headtest)  headw[5]++;
-			}
-		}
-		if(z1 == 0) dz1 = 0;
-		z1--;
-	}
-      }
-		/* BOTTOM SIDE */
-      for (igrow=1;igrow<=kplus;igrow++) {  
-	if(dz2){
-		ii = 0;
-		for(j=y1+1; j<=y2-1; j++){
-			for(i=x1+1; i<=x2-1; i++){
-				sort[ii].time = t0(i,j,z2-1);
-				sort[ii].i1 = i;
-				sort[ii].i2 = j;
-				ii++;
-			}
-		}
-		qsort((char *)sort,ii,sizeof(struct sorted),compar);
-		for(i=0;i<ii;i++){
-			X1 = sort[i].i1;
-			X2 = sort[i].i2;
-			index = z2*nxy + X2*nx + X1;
-			lasti = (z2-1)*nxy + X2*nx + X1;
-			fhead = 0.;
-			guess = time0[index];
-                        if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
-			   && time0[index+nx] < 1.e9 && X2<ny-1  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,X2,z2-1),
-				      t0(X1+1,X2,z2-1),t0(X1+1,X2+1,z2-1),t0(X1,X2+1,z2-1),
-				      t0(X1+1,X2,z2  ),t0(X1+1,X2+1,z2  ),t0(X1,X2+1,z2  ),
-				      s0(X1,X2,z2), s0(X1,X2,z2-1),
-				      s0(X1+1,X2,z2-1),s0(X1+1,X2+1,z2-1),s0(X1,X2+1,z2-1),
-				      s0(X1+1,X2,z2  ),s0(X1+1,X2+1,z2  ),s0(X1,X2+1,z2  ));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
-			   && time0[index+nx] < 1.e9 && X2<ny-1  && X1>0 ) {
-			  try = fdh3d(              t0(X1,X2,z2-1),
-				      t0(X1-1,X2,z2-1),t0(X1-1,X2+1,z2-1),t0(X1,X2+1,z2-1),
-				      t0(X1-1,X2,z2  ),t0(X1-1,X2+1,z2  ),t0(X1,X2+1,z2  ),
-				      s0(X1,X2,z2), s0(X1,X2,z2-1),
-				      s0(X1-1,X2,z2-1),s0(X1-1,X2+1,z2-1),s0(X1,X2+1,z2-1),
-				      s0(X1-1,X2,z2  ),s0(X1-1,X2+1,z2  ),s0(X1,X2+1,z2  ));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
-			   && time0[index-nx] < 1.e9 && X2>0  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,X2,z2-1),
-				      t0(X1+1,X2,z2-1),t0(X1+1,X2-1,z2-1),t0(X1,X2-1,z2-1),
-				      t0(X1+1,X2,z2  ),t0(X1+1,X2-1,z2  ),t0(X1,X2-1,z2  ),
-				      s0(X1,X2,z2), s0(X1,X2,z2-1),
-				      s0(X1+1,X2,z2-1),s0(X1+1,X2-1,z2-1),s0(X1,X2-1,z2-1),
-				      s0(X1+1,X2,z2  ),s0(X1+1,X2-1,z2  ),s0(X1,X2-1,z2  ));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
-			   && time0[index-nx] < 1.e9 && X2>0  && X1>0 ) {
-			  try = fdh3d(              t0(X1,X2,z2-1),
-				      t0(X1-1,X2,z2-1),t0(X1-1,X2-1,z2-1),t0(X1,X2-1,z2-1),
-				      t0(X1-1,X2,z2  ),t0(X1-1,X2-1,z2  ),t0(X1,X2-1,z2  ),
-				      s0(X1,X2,z2), s0(X1,X2,z2-1),
-				      s0(X1-1,X2,z2-1),s0(X1-1,X2-1,z2-1),s0(X1,X2-1,z2-1),
-				      s0(X1-1,X2,z2  ),s0(X1-1,X2-1,z2  ),s0(X1,X2-1,z2  ));
-			  if (try<guess) guess = try;
-			}
-                        if(guess > 1.0e9){ 
-			  if(time0[index+1] < 1.e9 && X1<nx-1 && X2>y1+1 && X2<y2-1 )  {
-			      try = fdhne(t0(X1,X2,z2-1),t0(X1+1,X2,z2-1),t0(X1+1,X2,z2),
-					  t0(X1+1,X2-1,z2-1),t0(X1+1,X2+1,z2-1),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1+1,X2,z2-1),s0(X1+1,X2,z2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 && X2>y1+1 && X2<y2-1 )  {
-			      try = fdhne(t0(X1,X2,z2-1),t0(X1-1,X2,z2-1),t0(X1-1,X2,z2),
-					  t0(X1-1,X2-1,z2-1),t0(X1-1,X2+1,z2-1),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1-1,X2,z2-1),s0(X1-1,X2,z2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nx] < 1.e9 && X2<ny-1 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,X2,z2-1),t0(X1,X2+1,z2-1),t0(X1,X2+1,z2),
-					  t0(X1-1,X2+1,z2-1),t0(X1+1,X2+1,z2-1),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1,X2+1,z2-1),s0(X1,X2+1,z2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,X2,z2-1),t0(X1,X2-1,z2-1),t0(X1,X2-1,z2),
-					  t0(X1-1,X2-1,z2-1),t0(X1+1,X2-1,z2-1),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1,X2-1,z2-1),s0(X1,X2-1,z2) );
-			    if (try<guess)  guess = try;
-			  }
-		        }
-			  if(time0[index+1] < 1.e9 && X1<nx-1 )  {
-			    try = fdh2d(t0(X1,X2,z2-1),t0(X1+1,X2,z2-1),t0(X1+1,X2,z2),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1+1,X2,z2-1),s0(X1+1,X2,z2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 )  {
-			    try = fdh2d(t0(X1,X2,z2-1),t0(X1-1,X2,z2-1),t0(X1-1,X2,z2),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1-1,X2,z2-1),s0(X1-1,X2,z2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nx] < 1.e9 && X2<ny-1 )  {
-			    try = fdh2d(t0(X1,X2,z2-1),t0(X1,X2+1,z2-1),t0(X1,X2+1,z2),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1,X2+1,z2-1),s0(X1,X2+1,z2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X2>0 )  {
-			    try = fdh2d(t0(X1,X2,z2-1),t0(X1,X2-1,z2-1),t0(X1,X2-1,z2),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1),s0(X1,X2-1,z2-1),s0(X1,X2-1,z2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index+nx+1] < 1.e9
-			     && time0[index+nx] < 1.e9 && X2<ny-1  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,X2,z2),t0(X1+1,X2+1,z2),t0(X1,X2+1,z2),
-					s0(X1,X2,z2),
-					s0(X1+1,X2,z2),s0(X1+1,X2+1,z2),s0(X1,X2+1,z2) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index-nx+1] < 1.e9
-			     && time0[index-nx] < 1.e9 && X2>0  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,X2,z2),t0(X1+1,X2-1,z2),t0(X1,X2-1,z2),
-					s0(X1,X2,z2),
-					s0(X1+1,X2,z2),s0(X1+1,X2-1,z2),s0(X1,X2-1,z2) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index+nx-1] < 1.e9
-			     && time0[index+nx] < 1.e9 && X2<ny-1  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,X2,z2),t0(X1-1,X2+1,z2),t0(X1,X2+1,z2),
-					s0(X1,X2,z2),
-					s0(X1-1,X2,z2),s0(X1-1,X2+1,z2),s0(X1,X2+1,z2) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index-nx-1] < 1.e9
-			     && time0[index-nx] < 1.e9 && X2>0  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,X2,z2),t0(X1-1,X2-1,z2),t0(X1,X2-1,z2),
-					s0(X1,X2,z2),
-					s0(X1-1,X2,z2),s0(X1-1,X2-1,z2),s0(X1,X2-1,z2) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if(guess > 1.0e9){ 
-			  if ( X1>x1+1 && X1<x2-1 && X2>y1+1 && X2<y2-1 ) {
-			    try = fdhnf(t0(X1,X2,z2-1),
-					  t0(X1+1,X2,z2-1),t0(X1,X2+1,z2-1),
-					  t0(X1-1,X2,z2-1),t0(X1,X2-1,z2-1),
-					  s0(X1,X2,z2),
-					  s0(X1,X2,z2-1) );
-			    if (try<guess)  guess = try;
-			  }
-			} 
-			  try = t0(X1,X2,z2-1) + .5*(s0(X1,X2,z2)+s0(X1,X2,z2-1));
-			  if (try<guess)  guess = try;
-                          if ( time0[index+1]<1.e9 && X1<nx-1 )  {
-			    try = t0(X1+1,X2,z2) + .5*(s0(X1,X2,z2)+s0(X1+1,X2,z2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-1]<1.e9 && X1>0 )  {
-			    try = t0(X1-1,X2,z2) + .5*(s0(X1,X2,z2)+s0(X1-1,X2,z2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index+nx]<1.e9 && X2<ny-1 )  {
-			    try = t0(X1,X2+1,z2) + .5*(s0(X1,X2,z2)+s0(X1,X2+1,z2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nx]<1.e9 && X2>0 )  {
-			    try = t0(X1,X2-1,z2) + .5*(s0(X1,X2,z2)+s0(X1,X2-1,z2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if (guess<time0[index]) {
-				time0[index] = guess;
-				if (fhead>headtest)  headw[6]++;
-			}
-		}
-		if(z2 == nz-1) dz2 = 0;
-		z2++;
-	}
-      }
-		/* FRONT SIDE */
-      for (igrow=1;igrow<=jminus;igrow++) {  
-	if(dy1){
-		ii = 0;
-		for(k=z1+1; k<=z2-1; k++){
-			for(i=x1+1; i<=x2-1; i++){
-				sort[ii].time = t0(i,y1+1,k);
-				sort[ii].i1 = i;
-				sort[ii].i2 = k;
-				ii++;
-			}
-		}
-		qsort((char *)sort,ii,sizeof(struct sorted),compar);
-		for(i=0;i<ii;i++){
-			X1 = sort[i].i1;
-			X2 = sort[i].i2;
-			index = X2*nxy + y1*nx + X1;
-			lasti = X2*nxy + (y1+1)*nx + X1;
-			fhead = 0.;
-			guess = time0[index];
-			if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,y1+1,X2),
-				      t0(X1+1,y1+1,X2),t0(X1+1,y1+1,X2+1),t0(X1,y1+1,X2+1),
-				      t0(X1+1,y1  ,X2),t0(X1+1,y1  ,X2+1),t0(X1,y1  ,X2+1),
-				      s0(X1,y1,X2), s0(X1,y1+1,X2),
-				      s0(X1+1,y1+1,X2),s0(X1+1,y1+1,X2+1),s0(X1,y1+1,X2+1),
-				      s0(X1+1,y1  ,X2),s0(X1+1,y1  ,X2+1),s0(X1,y1  ,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			  try = fdh3d(              t0(X1,y1+1,X2),
-				      t0(X1-1,y1+1,X2),t0(X1-1,y1+1,X2+1),t0(X1,y1+1,X2+1),
-				      t0(X1-1,y1  ,X2),t0(X1-1,y1  ,X2+1),t0(X1,y1  ,X2+1),
-				      s0(X1,y1,X2), s0(X1,y1+1,X2),
-				      s0(X1-1,y1+1,X2),s0(X1-1,y1+1,X2+1),s0(X1,y1+1,X2+1),
-				      s0(X1-1,y1  ,X2),s0(X1-1,y1  ,X2+1),s0(X1,y1  ,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,y1+1,X2),
-				      t0(X1+1,y1+1,X2),t0(X1+1,y1+1,X2-1),t0(X1,y1+1,X2-1),
-				      t0(X1+1,y1  ,X2),t0(X1+1,y1  ,X2-1),t0(X1,y1  ,X2-1),
-				      s0(X1,y1,X2), s0(X1,y1+1,X2),
-				      s0(X1+1,y1+1,X2),s0(X1+1,y1+1,X2-1),s0(X1,y1+1,X2-1),
-				      s0(X1+1,y1  ,X2),s0(X1+1,y1  ,X2-1),s0(X1,y1  ,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			  try = fdh3d(              t0(X1,y1+1,X2),
-				      t0(X1-1,y1+1,X2),t0(X1-1,y1+1,X2-1),t0(X1,y1+1,X2-1),
-				      t0(X1-1,y1  ,X2),t0(X1-1,y1  ,X2-1),t0(X1,y1  ,X2-1),
-				      s0(X1,y1,X2), s0(X1,y1+1,X2),
-				      s0(X1-1,y1+1,X2),s0(X1-1,y1+1,X2-1),s0(X1,y1+1,X2-1),
-				      s0(X1-1,y1  ,X2),s0(X1-1,y1  ,X2-1),s0(X1,y1  ,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(guess > 1.0e9){ 
-			  if(time0[index+1] < 1.e9 && X1<nx-1 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(X1,y1+1,X2),t0(X1+1,y1+1,X2),t0(X1+1,y1,X2),
-					  t0(X1+1,y1+1,X2-1),t0(X1+1,y1+1,X2+1),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1+1,y1+1,X2),s0(X1+1,y1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(X1,y1+1,X2),t0(X1-1,y1+1,X2),t0(X1-1,y1,X2),
-					  t0(X1-1,y1+1,X2-1),t0(X1-1,y1+1,X2+1),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1-1,y1+1,X2),s0(X1-1,y1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,y1+1,X2),t0(X1,y1+1,X2+1),t0(X1,y1,X2+1),
-					  t0(X1-1,y1+1,X2+1),t0(X1+1,y1+1,X2+1),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1,y1+1,X2+1),s0(X1,y1,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,y1+1,X2),t0(X1,y1+1,X2-1),t0(X1,y1,X2-1),
-					  t0(X1-1,y1+1,X2-1),t0(X1+1,y1+1,X2-1),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1,y1+1,X2-1),s0(X1,y1,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-		        } 
-			  if(time0[index+1] < 1.e9 && X1<nx-1 )  {
-			    try = fdh2d(t0(X1,y1+1,X2),t0(X1+1,y1+1,X2),t0(X1+1,y1,X2),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1+1,y1+1,X2),s0(X1+1,y1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 )  {
-			    try = fdh2d(t0(X1,y1+1,X2),t0(X1-1,y1+1,X2),t0(X1-1,y1,X2),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1-1,y1+1,X2),s0(X1-1,y1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 )  {
-			    try = fdh2d(t0(X1,y1+1,X2),t0(X1,y1+1,X2+1),t0(X1,y1,X2+1),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1,y1+1,X2+1),s0(X1,y1,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 )  {
-			    try = fdh2d(t0(X1,y1+1,X2),t0(X1,y1+1,X2-1),t0(X1,y1,X2-1),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2),s0(X1,y1+1,X2-1),s0(X1,y1,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,y1,X2),t0(X1+1,y1,X2+1),t0(X1,y1,X2+1),
-					s0(X1,y1,X2),
-					s0(X1+1,y1,X2),s0(X1+1,y1,X2+1),s0(X1,y1,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,y1,X2),t0(X1+1,y1,X2-1),t0(X1,y1,X2-1),
-					s0(X1,y1,X2),
-					s0(X1+1,y1,X2),s0(X1+1,y1,X2-1),s0(X1,y1,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,y1,X2),t0(X1-1,y1,X2+1),t0(X1,y1,X2+1),
-					s0(X1,y1,X2),
-					s0(X1-1,y1,X2),s0(X1-1,y1,X2+1),s0(X1,y1,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,y1,X2),t0(X1-1,y1,X2-1),t0(X1,y1,X2-1),
-					s0(X1,y1,X2),
-					s0(X1-1,y1,X2),s0(X1-1,y1,X2-1),s0(X1,y1,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if(guess > 1.0e9){ 
-			  if ( X1>x1+1 && X1<x2-1 && X2>z1+1 && X2<z2-1 ) {
-			    try = fdhnf(t0(X1,y1+1,X2),
-					  t0(X1+1,y1+1,X2),t0(X1,y1+1,X2+1),
-					  t0(X1-1,y1+1,X2),t0(X1,y1+1,X2-1),
-					  s0(X1,y1,X2),
-					  s0(X1,y1+1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			} 
-			  try = t0(X1,y1+1,X2) + .5*(s0(X1,y1,X2)+s0(X1,y1+1,X2));
-			  if (try<guess)  guess = try;
-                          if ( time0[index+1]<1.e9 && X1<nx-1 )  {
-			    try = t0(X1+1,y1,X2) + .5*(s0(X1,y1,X2)+s0(X1+1,y1,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-1]<1.e9 && X1>0 )  {
-			    try = t0(X1-1,y1,X2) + .5*(s0(X1,y1,X2)+s0(X1-1,y1,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index+nxy]<1.e9 && X2<nz-1 )  {
-			    try = t0(X1,y1,X2+1) + .5*(s0(X1,y1,X2)+s0(X1,y1,X2+1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nxy]<1.e9 && X2>0 )  {
-			    try = t0(X1,y1,X2-1) + .5*(s0(X1,y1,X2)+s0(X1,y1,X2-1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if (guess<time0[index]) {
-				time0[index] = guess;
-				if (fhead>headtest)  headw[3]++;
-			}
-		}
-		if(y1 == 0) dy1 = 0;
-		y1--;
-	}
-      }
-		/* BACK SIDE */
-      for (igrow=1;igrow<=jplus;igrow++) {  
-	if(dy2){
-		ii = 0;
-		for(k=z1+1; k<=z2-1; k++){
-			for(i=x1+1; i<=x2-1; i++){
-				sort[ii].time = t0(i,y2-1,k);
-				sort[ii].i1 = i;
-				sort[ii].i2 = k;
-				ii++;
-			}
-		}
-		qsort((char *)sort,ii,sizeof(struct sorted),compar);
-		for(i=0;i<ii;i++){
-			X1 = sort[i].i1;
-			X2 = sort[i].i2;
-			index = X2*nxy + y2*nx + X1;
-			lasti = X2*nxy + (y2-1)*nx + X1;
-			fhead = 0.;
-			guess = time0[index];
-			if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,y2-1,X2),
-				      t0(X1+1,y2-1,X2),t0(X1+1,y2-1,X2+1),t0(X1,y2-1,X2+1),
-				      t0(X1+1,y2  ,X2),t0(X1+1,y2  ,X2+1),t0(X1,y2  ,X2+1),
-				      s0(X1,y2,X2), s0(X1,y2-1,X2),
-				      s0(X1+1,y2-1,X2),s0(X1+1,y2-1,X2+1),s0(X1,y2-1,X2+1),
-				      s0(X1+1,y2  ,X2),s0(X1+1,y2  ,X2+1),s0(X1,y2  ,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			  try = fdh3d(              t0(X1,y2-1,X2),
-				      t0(X1-1,y2-1,X2),t0(X1-1,y2-1,X2+1),t0(X1,y2-1,X2+1),
-				      t0(X1-1,y2  ,X2),t0(X1-1,y2  ,X2+1),t0(X1,y2  ,X2+1),
-				      s0(X1,y2,X2), s0(X1,y2-1,X2),
-				      s0(X1-1,y2-1,X2),s0(X1-1,y2-1,X2+1),s0(X1,y2-1,X2+1),
-				      s0(X1-1,y2  ,X2),s0(X1-1,y2  ,X2+1),s0(X1,y2  ,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1<nx-1 ) {
-			  try = fdh3d(              t0(X1,y2-1,X2),
-				      t0(X1+1,y2-1,X2),t0(X1+1,y2-1,X2-1),t0(X1,y2-1,X2-1),
-				      t0(X1+1,y2  ,X2),t0(X1+1,y2  ,X2-1),t0(X1,y2  ,X2-1),
-				      s0(X1,y2,X2), s0(X1,y2-1,X2),
-				      s0(X1+1,y2-1,X2),s0(X1+1,y2-1,X2-1),s0(X1,y2-1,X2-1),
-				      s0(X1+1,y2  ,X2),s0(X1+1,y2  ,X2-1),s0(X1,y2  ,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			  try = fdh3d(              t0(X1,y2-1,X2),
-				      t0(X1-1,y2-1,X2),t0(X1-1,y2-1,X2-1),t0(X1,y2-1,X2-1),
-				      t0(X1-1,y2  ,X2),t0(X1-1,y2  ,X2-1),t0(X1,y2  ,X2-1),
-				      s0(X1,y2,X2), s0(X1,y2-1,X2),
-				      s0(X1-1,y2-1,X2),s0(X1-1,y2-1,X2-1),s0(X1,y2-1,X2-1),
-				      s0(X1-1,y2  ,X2),s0(X1-1,y2  ,X2-1),s0(X1,y2  ,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(guess > 1.0e9){ 
-			  if(time0[index+1] < 1.e9 && X1<nx-1 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(X1,y2-1,X2),t0(X1+1,y2-1,X2),t0(X1+1,y2,X2),
-					  t0(X1+1,y2-1,X2-1),t0(X1+1,y2-1,X2+1),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1+1,y2-1,X2),s0(X1+1,y2,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(X1,y2-1,X2),t0(X1-1,y2-1,X2),t0(X1-1,y2,X2),
-					  t0(X1-1,y2-1,X2-1),t0(X1-1,y2-1,X2+1),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1-1,y2-1,X2),s0(X1-1,y2,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,y2-1,X2),t0(X1,y2-1,X2+1),t0(X1,y2,X2+1),
-					  t0(X1-1,y2-1,X2+1),t0(X1+1,y2-1,X2+1),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1,y2-1,X2+1),s0(X1,y2,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 && X1>x1+1 && X1<x2-1 )  {
-			      try = fdhne(t0(X1,y2-1,X2),t0(X1,y2-1,X2-1),t0(X1,y2,X2-1),
-					  t0(X1-1,y2-1,X2-1),t0(X1+1,y2-1,X2-1),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1,y2-1,X2-1),s0(X1,y2,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-		        } 
-			  if(time0[index+1] < 1.e9 && X1<nx-1 )  {
-			    try = fdh2d(t0(X1,y2-1,X2),t0(X1+1,y2-1,X2),t0(X1+1,y2,X2),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1+1,y2-1,X2),s0(X1+1,y2,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-1] < 1.e9 && X1>0 )  {
-			    try = fdh2d(t0(X1,y2-1,X2),t0(X1-1,y2-1,X2),t0(X1-1,y2,X2),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1-1,y2-1,X2),s0(X1-1,y2,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 )  {
-			    try = fdh2d(t0(X1,y2-1,X2),t0(X1,y2-1,X2+1),t0(X1,y2,X2+1),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1,y2-1,X2+1),s0(X1,y2,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 )  {
-			    try = fdh2d(t0(X1,y2-1,X2),t0(X1,y2-1,X2-1),t0(X1,y2,X2-1),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2),s0(X1,y2-1,X2-1),s0(X1,y2,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index+nxy+1] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,y2,X2),t0(X1+1,y2,X2+1),t0(X1,y2,X2+1),
-					s0(X1,y2,X2),
-					s0(X1+1,y2,X2),s0(X1+1,y2,X2+1),s0(X1,y2,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index+1] < 1.e9 && time0[index-nxy+1] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1<nx-1 ) {
-			    try = fdh2d(t0(X1+1,y2,X2),t0(X1+1,y2,X2-1),t0(X1,y2,X2-1),
-					s0(X1,y2,X2),
-					s0(X1+1,y2,X2),s0(X1+1,y2,X2-1),s0(X1,y2,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index+nxy-1] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,y2,X2),t0(X1-1,y2,X2+1),t0(X1,y2,X2+1),
-					s0(X1,y2,X2),
-					s0(X1-1,y2,X2),s0(X1-1,y2,X2+1),s0(X1,y2,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-1] < 1.e9 && time0[index-nxy-1] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			    try = fdh2d(t0(X1-1,y2,X2),t0(X1-1,y2,X2-1),t0(X1,y2,X2-1),
-					s0(X1,y2,X2),
-					s0(X1-1,y2,X2),s0(X1-1,y2,X2-1),s0(X1,y2,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if(guess > 1.0e9){ 
-			  if ( X1>x1+1 && X1<x2-1 && X2>z1+1 && X2<z2-1 ) {
-			    try = fdhnf(t0(X1,y2-1,X2),
-					  t0(X1+1,y2-1,X2),t0(X1,y2-1,X2+1),
-					  t0(X1-1,y2-1,X2),t0(X1,y2-1,X2-1),
-					  s0(X1,y2,X2),
-					  s0(X1,y2-1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			} 
-			  try = t0(X1,y2-1,X2) + .5*(s0(X1,y2,X2)+s0(X1,y2-1,X2));
-			  if (try<guess)  guess = try;
-                          if ( time0[index+1]<1.e9 && X1<nx-1 )  {
-			    try = t0(X1+1,y2,X2) + .5*(s0(X1,y2,X2)+s0(X1+1,y2,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-1]<1.e9 && X1>0 )  {
-			    try = t0(X1-1,y2,X2) + .5*(s0(X1,y2,X2)+s0(X1-1,y2,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index+nxy]<1.e9 && X2<nz-1 )  {
-			    try = t0(X1,y2,X2+1) + .5*(s0(X1,y2,X2)+s0(X1,y2,X2+1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nxy]<1.e9 && X2>0 )  {
-			    try = t0(X1,y2,X2-1) + .5*(s0(X1,y2,X2)+s0(X1,y2,X2-1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if (guess<time0[index]) {
-				time0[index] = guess;
-				if (fhead>headtest)  headw[4]++;
-			}
-		}
-		if(y2 == ny-1) dy2 = 0;
-		y2++;
-	}
-      }
-		/* LEFT SIDE */
-      for (igrow=1;igrow<=iminus;igrow++) {  
-	if(dx1){
-		ii = 0;
-		for(k=z1+1; k<=z2-1; k++){
-			for(j=y1+1; j<=y2-1; j++){
-				sort[ii].time = t0(x1+1,j,k);
-				sort[ii].i1 = j;
-				sort[ii].i2 = k;
-				ii++;
-			}
-		}
-		qsort((char *)sort,ii,sizeof(struct sorted),compar);
-		for(i=0;i<ii;i++){
-			X1 = sort[i].i1;
-			X2 = sort[i].i2;
-			index = X2*nxy + X1*nx + x1;
-			lasti = X2*nxy + X1*nx + (x1+1);
-			fhead = 0.;
-			guess = time0[index];
-			if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<ny-1 ) {
-			  try = fdh3d(              t0(x1+1,X1,X2),
-				      t0(x1+1,X1+1,X2),t0(x1+1,X1+1,X2+1),t0(x1+1,X1,X2+1),
-				      t0(x1  ,X1+1,X2),t0(x1  ,X1+1,X2+1),t0(x1  ,X1,X2+1),
-				      s0(x1,X1,X2), s0(x1+1,X1,X2),
-				      s0(x1+1,X1+1,X2),s0(x1+1,X1+1,X2+1),s0(x1+1,X1,X2+1),
-				      s0(x1  ,X1+1,X2),s0(x1  ,X1+1,X2+1),s0(x1  ,X1,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			  try = fdh3d(              t0(x1+1,X1,X2),
-				      t0(x1+1,X1-1,X2),t0(x1+1,X1-1,X2+1),t0(x1+1,X1,X2+1),
-				      t0(x1  ,X1-1,X2),t0(x1  ,X1-1,X2+1),t0(x1  ,X1,X2+1),
-				      s0(x1,X1,X2), s0(x1+1,X1,X2),
-				      s0(x1+1,X1-1,X2),s0(x1+1,X1-1,X2+1),s0(x1+1,X1,X2+1),
-				      s0(x1  ,X1-1,X2),s0(x1  ,X1-1,X2+1),s0(x1  ,X1,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1<ny-1 ) {
-			  try = fdh3d(              t0(x1+1,X1,X2),
-				      t0(x1+1,X1+1,X2),t0(x1+1,X1+1,X2-1),t0(x1+1,X1,X2-1),
-				      t0(x1  ,X1+1,X2),t0(x1  ,X1+1,X2-1),t0(x1  ,X1,X2-1),
-				      s0(x1,X1,X2), s0(x1+1,X1,X2),
-				      s0(x1+1,X1+1,X2),s0(x1+1,X1+1,X2-1),s0(x1+1,X1,X2-1),
-				      s0(x1  ,X1+1,X2),s0(x1  ,X1+1,X2-1),s0(x1  ,X1,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			  try = fdh3d(              t0(x1+1,X1,X2),
-				      t0(x1+1,X1-1,X2),t0(x1+1,X1-1,X2-1),t0(x1+1,X1,X2-1),
-				      t0(x1  ,X1-1,X2),t0(x1  ,X1-1,X2-1),t0(x1  ,X1,X2-1),
-				      s0(x1,X1,X2), s0(x1+1,X1,X2),
-				      s0(x1+1,X1-1,X2),s0(x1+1,X1-1,X2-1),s0(x1+1,X1,X2-1),
-				      s0(x1  ,X1-1,X2),s0(x1  ,X1-1,X2-1),s0(x1  ,X1,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(guess > 1.0e9){ 
-			  if(time0[index+nx] < 1.e9 && X1<ny-1 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1+1,X2),t0(x1,X1+1,X2),
-					  t0(x1+1,X1+1,X2-1),t0(x1+1,X1+1,X2+1),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1+1,X2),s0(x1,X1+1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1-1,X2),t0(x1,X1-1,X2),
-					  t0(x1+1,X1-1,X2-1),t0(x1+1,X1-1,X2+1),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1-1,X2),s0(x1,X1-1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>y1+1 && X1<y2-1 )  {
-			      try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1,X2+1),t0(x1,X1,X2+1),
-					  t0(x1+1,X1-1,X2+1),t0(x1+1,X1+1,X2+1),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1,X2+1),s0(x1,X1,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 && X1>y1+1 && X1<y2-1 )  {
-			      try = fdhne(t0(x1+1,X1,X2),t0(x1+1,X1,X2-1),t0(x1,X1,X2-1),
-					  t0(x1+1,X1-1,X2-1),t0(x1+1,X1+1,X2-1),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1,X2-1),s0(x1,X1,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-		        } 
-			  if(time0[index+nx] < 1.e9 && X1<ny-1 )  {
-			    try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1+1,X2),t0(x1,X1+1,X2),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1+1,X2),s0(x1,X1+1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X1>0 )  {
-			    try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1-1,X2),t0(x1,X1-1,X2),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1-1,X2),s0(x1,X1-1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 )  {
-			    try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1,X2+1),t0(x1,X1,X2+1),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1,X2+1),s0(x1,X1,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 )  {
-			    try = fdh2d(t0(x1+1,X1,X2),t0(x1+1,X1,X2-1),t0(x1,X1,X2-1),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2),s0(x1+1,X1,X2-1),s0(x1,X1,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<ny-1 ) {
-			    try = fdh2d(t0(x1,X1+1,X2),t0(x1,X1+1,X2+1),t0(x1,X1,X2+1),
-					s0(x1,X1,X2),
-					s0(x1,X1+1,X2),s0(x1,X1+1,X2+1),s0(x1,X1,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1<ny-1 ) {
-			    try = fdh2d(t0(x1,X1+1,X2),t0(x1,X1+1,X2-1),t0(x1,X1,X2-1),
-					s0(x1,X1,X2),
-					s0(x1,X1+1,X2),s0(x1,X1+1,X2-1),s0(x1,X1,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			    try = fdh2d(t0(x1,X1-1,X2),t0(x1,X1-1,X2+1),t0(x1,X1,X2+1),
-					s0(x1,X1,X2),
-					s0(x1,X1-1,X2),s0(x1,X1-1,X2+1),s0(x1,X1,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			    try = fdh2d(t0(x1,X1-1,X2),t0(x1,X1-1,X2-1),t0(x1,X1,X2-1),
-					s0(x1,X1,X2),
-					s0(x1,X1-1,X2),s0(x1,X1-1,X2-1),s0(x1,X1,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if(guess > 1.0e9){ 
-			  if ( X1>y1+1 && X1<y2-1 && X2>z1+1 && X2<z2-1 ) {
-			    try = fdhnf(t0(x1+1,X1,X2),
-					  t0(x1+1,X1+1,X2),t0(x1+1,X1,X2+1),
-					  t0(x1+1,X1-1,X2),t0(x1+1,X1,X2-1),
-					  s0(x1,X1,X2),
-					  s0(x1+1,X1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			} 
-			  try = t0(x1+1,X1,X2) + .5*(s0(x1,X1,X2)+s0(x1+1,X1,X2));
-			  if (try<guess)  guess = try;
-                          if ( time0[index+nx]<1.e9 && X1<ny-1 )  {
-			    try = t0(x1,X1+1,X2) + .5*(s0(x1,X1,X2)+s0(x1,X1+1,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nx]<1.e9 && X1>0 )  {
-			    try = t0(x1,X1-1,X2) + .5*(s0(x1,X1,X2)+s0(x1,X1-1,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index+nxy]<1.e9 && X2<nz-1 )  {
-			    try = t0(x1,X1,X2+1) + .5*(s0(x1,X1,X2)+s0(x1,X1,X2+1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nxy]<1.e9 && X2>0 )  {
-			    try = t0(x1,X1,X2-1) + .5*(s0(x1,X1,X2)+s0(x1,X1,X2-1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if (guess<time0[index]) {
-				time0[index] = guess;
-				if (fhead>headtest)  headw[1]++;
-			}
-		}
-		if(x1 == 0) dx1 = 0;
-		x1--;
-	}
-      }
-		/* RIGHT SIDE */
-      for (igrow=1;igrow<=iplus;igrow++) {  
-	if(dx2){
-		ii = 0;
-		for(k=z1+1; k<=z2-1; k++){
-			for(j=y1+1; j<=y2-1; j++){
-				sort[ii].time = t0(x2-1,j,k);
-				sort[ii].i1 = j;
-				sort[ii].i2 = k;
-				ii++;
-			}
-		}
-		qsort((char *)sort,ii,sizeof(struct sorted),compar);
-		for(i=0;i<ii;i++){
-			X1 = sort[i].i1;
-			X2 = sort[i].i2;
-			index = X2*nxy + X1*nx + x2;
-			lasti = X2*nxy + X1*nx + (x2-1);
-			fhead = 0.;
-			guess = time0[index];
-			if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<ny-1 ) {
-			  try = fdh3d(              t0(x2-1,X1,X2),
-				      t0(x2-1,X1+1,X2),t0(x2-1,X1+1,X2+1),t0(x2-1,X1,X2+1),
-				      t0(x2  ,X1+1,X2),t0(x2  ,X1+1,X2+1),t0(x2  ,X1,X2+1),
-				      s0(x2,X1,X2), s0(x2-1,X1,X2),
-				      s0(x2-1,X1+1,X2),s0(x2-1,X1+1,X2+1),s0(x2-1,X1,X2+1),
-				      s0(x2  ,X1+1,X2),s0(x2  ,X1+1,X2+1),s0(x2  ,X1,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
-			   && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			  try = fdh3d(              t0(x2-1,X1,X2),
-				      t0(x2-1,X1-1,X2),t0(x2-1,X1-1,X2+1),t0(x2-1,X1,X2+1),
-				      t0(x2  ,X1-1,X2),t0(x2  ,X1-1,X2+1),t0(x2  ,X1,X2+1),
-				      s0(x2,X1,X2), s0(x2-1,X1,X2),
-				      s0(x2-1,X1-1,X2),s0(x2-1,X1-1,X2+1),s0(x2-1,X1,X2+1),
-				      s0(x2  ,X1-1,X2),s0(x2  ,X1-1,X2+1),s0(x2  ,X1,X2+1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1<ny-1 ) {
-			  try = fdh3d(              t0(x2-1,X1,X2),
-				      t0(x2-1,X1+1,X2),t0(x2-1,X1+1,X2-1),t0(x2-1,X1,X2-1),
-				      t0(x2  ,X1+1,X2),t0(x2  ,X1+1,X2-1),t0(x2  ,X1,X2-1),
-				      s0(x2,X1,X2), s0(x2-1,X1,X2),
-				      s0(x2-1,X1+1,X2),s0(x2-1,X1+1,X2-1),s0(x2-1,X1,X2-1),
-				      s0(x2  ,X1+1,X2),s0(x2  ,X1+1,X2-1),s0(x2  ,X1,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
-			   && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			  try = fdh3d(              t0(x2-1,X1,X2),
-				      t0(x2-1,X1-1,X2),t0(x2-1,X1-1,X2-1),t0(x2-1,X1,X2-1),
-				      t0(x2  ,X1-1,X2),t0(x2  ,X1-1,X2-1),t0(x2  ,X1,X2-1),
-				      s0(x2,X1,X2), s0(x2-1,X1,X2),
-				      s0(x2-1,X1-1,X2),s0(x2-1,X1-1,X2-1),s0(x2-1,X1,X2-1),
-				      s0(x2  ,X1-1,X2),s0(x2  ,X1-1,X2-1),s0(x2  ,X1,X2-1));
-			  if (try<guess) guess = try;
-			}
-			if(guess > 1.0e9){ 
-			  if(time0[index+nx] < 1.e9 && X1<ny-1 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1+1,X2),t0(x2,X1+1,X2),
-					  t0(x2-1,X1+1,X2-1),t0(x2-1,X1+1,X2+1),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1+1,X2),s0(x2,X1+1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X1>0 && X2>z1+1 && X2<z2-1 )  {
-			      try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1-1,X2),t0(x2,X1-1,X2),
-					  t0(x2-1,X1-1,X2-1),t0(x2-1,X1-1,X2+1),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1-1,X2),s0(x2,X1-1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 && X1>y1+1 && X1<y2-1 )  {
-			      try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1,X2+1),t0(x2,X1,X2+1),
-					  t0(x2-1,X1-1,X2+1),t0(x2-1,X1+1,X2+1),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1,X2+1),s0(x2,X1,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 && X1>y1+1 && X1<y2-1 )  {
-			      try = fdhne(t0(x2-1,X1,X2),t0(x2-1,X1,X2-1),t0(x2,X1,X2-1),
-					  t0(x2-1,X1-1,X2-1),t0(x2-1,X1+1,X2-1),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1,X2-1),s0(x2,X1,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-		        } 
-			  if(time0[index+nx] < 1.e9 && X1<ny-1 )  {
-			    try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1+1,X2),t0(x2,X1+1,X2),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1+1,X2),s0(x2,X1+1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nx] < 1.e9 && X1>0 )  {
-			    try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1-1,X2),t0(x2,X1-1,X2),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1-1,X2),s0(x2,X1-1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nxy] < 1.e9 && X2<nz-1 )  {
-			    try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1,X2+1),t0(x2,X1,X2+1),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1,X2+1),s0(x2,X1,X2+1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index-nxy] < 1.e9 && X2>0 )  {
-			    try = fdh2d(t0(x2-1,X1,X2),t0(x2-1,X1,X2-1),t0(x2,X1,X2-1),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2),s0(x2-1,X1,X2-1),s0(x2,X1,X2-1) );
-			    if (try<guess)  guess = try;
-			  }
-			  if(time0[index+nx] < 1.e9 && time0[index+nxy+nx] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1<ny-1 ) {
-			    try = fdh2d(t0(x2,X1+1,X2),t0(x2,X1+1,X2+1),t0(x2,X1,X2+1),
-					s0(x2,X1,X2),
-					s0(x2,X1+1,X2),s0(x2,X1+1,X2+1),s0(x2,X1,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index+nx] < 1.e9 && time0[index-nxy+nx] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1<ny-1 ) {
-			    try = fdh2d(t0(x2,X1+1,X2),t0(x2,X1+1,X2-1),t0(x2,X1,X2-1),
-					s0(x2,X1,X2),
-					s0(x2,X1+1,X2),s0(x2,X1+1,X2-1),s0(x2,X1,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-nx] < 1.e9 && time0[index+nxy-nx] < 1.e9
-			     && time0[index+nxy] < 1.e9 && X2<nz-1  && X1>0 ) {
-			    try = fdh2d(t0(x2,X1-1,X2),t0(x2,X1-1,X2+1),t0(x2,X1,X2+1),
-					s0(x2,X1,X2),
-					s0(x2,X1-1,X2),s0(x2,X1-1,X2+1),s0(x2,X1,X2+1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if(time0[index-nx] < 1.e9 && time0[index-nxy-nx] < 1.e9
-			     && time0[index-nxy] < 1.e9 && X2>0  && X1>0 ) {
-			    try = fdh2d(t0(x2,X1-1,X2),t0(x2,X1-1,X2-1),t0(x2,X1,X2-1),
-					s0(x2,X1,X2),
-					s0(x2,X1-1,X2),s0(x2,X1-1,X2-1),s0(x2,X1,X2-1) );
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if(guess > 1.0e9){ 
-			  if ( X1>y1+1 && X1<y2-1 && X2>z1+1 && X2<z2-1 ) {
-			    try = fdhnf(t0(x2-1,X1,X2),
-					  t0(x2-1,X1+1,X2),t0(x2-1,X1,X2+1),
-					  t0(x2-1,X1-1,X2),t0(x2-1,X1,X2-1),
-					  s0(x2,X1,X2),
-					  s0(x2-1,X1,X2) );
-			    if (try<guess)  guess = try;
-			  }
-			} 
-			  try = t0(x2-1,X1,X2) + .5*(s0(x2,X1,X2)+s0(x2-1,X1,X2));
-			  if (try<guess)  guess = try;
-                          if ( time0[index+nx]<1.e9 && X1<ny-1 )  {
-			    try = t0(x2,X1+1,X2) + .5*(s0(x2,X1,X2)+s0(x2,X1+1,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nx]<1.e9 && X1>0 )  {
-			    try = t0(x2,X1-1,X2) + .5*(s0(x2,X1,X2)+s0(x2,X1-1,X2));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index+nxy]<1.e9 && X2<nz-1 )  {
-			    try = t0(x2,X1,X2+1) + .5*(s0(x2,X1,X2)+s0(x2,X1,X2+1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			  if ( time0[index-nxy]<1.e9 && X2>0 )  {
-			    try = t0(x2,X1,X2-1) + .5*(s0(x2,X1,X2)+s0(x2,X1,X2-1));
-			    if (try<guess)  {fhead=(guess-try)/slow0[index]; guess=try;}
-			  }
-			if (guess<time0[index]) {
-				time0[index] = guess;
-				if (fhead>headtest)  headw[2]++;
-			}
-		}
-		if(x2 == nx-1) dx2 = 0;
-		x2++;
-	}
-      }
-
-		/* UPDATE RADIUS */
-		radius++;
-		if(radius%10 == 0 && verbose) vmess("Completed radius = %li",radius);
-        if(radius == maxrad) rad0 = 0;
-
-	}	/* END BIG LOOP */
-
-
-	/* TEST IF REVERSE PROPAGATION IS NEEDED */
-
-	if (headw[1]==0 && headw[2]==0 && headw[3]==0 && headw[4]==0 
-		     && headw[5]==0 && headw[6]==0)
-		reverse=0;
-	else {
-		head=0;
-		if (headw[1]>0) {
-			if(verbose) vmess("Head waves found on left: %li",headw[1]);
-			if (headw[1]>head)  {
-				head = headw[1];
-				srcwall = 1;
-			}
-		}
-		if (headw[2]>0) {
-			if(verbose) vmess("Head waves found on right: %li",headw[2]);
-			if (headw[2]>head)  {
-				head = headw[2];
-				srcwall = 2;
-			}
-		}
-		if (headw[3]>0) {
-			if(verbose) vmess("Head waves found on front: %li",headw[3]);
-			if (headw[3]>head)  {
-				head = headw[3];
-				srcwall = 3;
-			}
-		}
-		if (headw[4]>0) {
-			if(verbose) vmess("Head waves found on back: %li",headw[4]);
-			if (headw[4]>head)  {
-				head = headw[4];
-				srcwall = 4;
-			}
-		}
-		if (headw[5]>0) {
-			if(verbose) vmess("Head waves found on top: %li",headw[5]);
-			if (headw[5]>head)  {
-				head = headw[5];
-				srcwall = 5;
-			}
-		}
-		if (headw[6]>0) {
-			if(verbose) vmess("Head waves found on bottom: %li",headw[6]);
-			if (headw[6]>head)  {
-				head = headw[6];
-				srcwall = 6;
-			}
-		}
-		if (headpref>0 && headw[headpref]>0) {
-			if(verbose) 
-				vmess("Preference to restart on wall opposite source");
-			srcwall = headpref;
-		}
-		/* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE */
-		dx1=1; dx2=1; dy1=1; dy2=1; dz1=1; dz2=1; rad0=1;
-		radius = 1;
-		if (srcwall == 1)	{  x2=1;
-			vmess("RESTART at left side of model");  }
-		else	{  x2=nx;	dx2=0;  }
-		if (srcwall == 2)	{ x1=nx-2;
-			vmess("RESTART at right side of model");  }
-		else	{  x1= -1;	dx1=0;  }
-		if (srcwall == 3)	{ y2=1;
-			vmess("RESTART at front side of model");  }
-		else	{  y2=ny;	dy2=0;  }
-		if (srcwall == 4)	{ y1=ny-2;
-			vmess("RESTART at back side of model");  }
-		else	{  y1= -1;	dy1=0;  }
-		if (srcwall == 5)	{ z2=1;
-			vmess("RESTART at top side of model");  }
-		else	{  z2=nz;	dz2=0;  }
-		if (srcwall == 6)	{ z1=nz-2;
-			vmess("RESTART at bottom side of model");  }
-		else	{  z1= -1;	dz1=0;  }
-		if (reverse == 0)  
-			vwarn("RESTART CANCELLED by choice of input parameter `reverse`");
-	}
-	reverse--;
-
-	}	/* END BIGGER LOOP - HOLE */
-
-	free(sort);
-	free(wall);
-}
-
-// int compar(struct sorted *a,struct sorted *b)
-// {
-// 	if(a->time > b->time) return(1);
-// 	if(b->time > a->time) return(-1);
-// 	else return(0);
-// }
-int compar(const void * a, const void * b)
-{
-	struct sorted *A = (struct sorted *)a;
-	struct sorted *B = (struct sorted *)b;
-	if(A->time > B->time) return(1);
-	if(B->time > A->time) return(-1);
-	else return(0);
-}
-
-
-/* 3D TRANSMISSION STENCIL
-   STENCIL FROM VIDALE; CONDITIONS AND OTHER OPTIONS FROM HOLE
-   JAH 11/91 */
-float fdh3d(t1,t2,t3,t4,t5,t6,t7,ss0,s1,s2,s3,s4,s5,s6,s7)
-     float  t1,t2,t3,t4,t5,t6,t7,ss0,s1,s2,s3,s4,s5,s6,s7;
-     /* ss0 at newpoint; s1,t1 adjacent on oldface;
-	s2,t2 and s4,t4 on oldface adjacent to s1;
-	s3,t3 on oldface diametrically opposite newpoint;
-	s5,t5 on newface adjacent to newpoint AND to s2;
-	s6,t6 on newface diagonal to newpoint (adjacent to s3);
-	s7,t7 on newface adjacent to newpoint AND to s4
-	*/
-{
-  float x,slo;
-  double sqrt();
-  slo = .125*(ss0+s1+s2+s3+s4+s5+s6+s7);
-  x = 6.*slo*slo - (t4-t2)*(t4-t2) - (t2-t6)*(t2-t6) - (t6-t4)*(t6-t4)
-                 - (t7-t5)*(t7-t5) - (t5-t1)*(t5-t1) - (t1-t7)*(t1-t7);
-  if (x>=0.)  {
-    x = t3 + sqrt(.5*x);
-    if ( (x<t1) || (x<t2) || (x<t4) || (x<t5) || (x<t6) || (x<t7) )  
-      x = 1.e11;   /* ACAUSAL; ABORT */
-  }
-  else  x = 1.e11;   /* SQRT IMAGINARY; ABORT */
-  return(x);
-}
-
-/* 3D STENCIL FOR NEW EDGE
-   STENCIL FROM VIDALE; CONDITIONS AND OTHER OPTIONS FROM HOLE
-   JAH 11/91 */
-float fdhne(t1,t2,t3,t4,t5,ss0,s1,s2,s3)
-     float  t1,t2,t3,t4,t5,ss0,s1,s2,s3;
-     /* ss0 at newpoint; s1,t1 adjacent on oldface;
-	s2,t2 diagonal on oldface; s3,t3 adjacent on newface;
-	t4,t5 beside t2 on old face opposite each other */
-{
-  float x,slo;
-  double sqrt();
-  slo = .25*(ss0+s1+s2+s3);
-  x = 2.*slo*slo - (t3-t1)*(t3-t1) - .5*(t5-t4)*(t5-t4);
-  if (x>=0.)  {
-    x = t2 + sqrt(x);
-    if ( (x<t1) || (x<t3) || (x<t4) || (x<t5) )     /* ACAUSAL; ABORT */
-      x = 1.e11;
-  }
-  else  x = 1.e11;   /* SQRT IMAGINARY; ABORT */
-  return(x);
-}
-
-/* 2D TRANSMISSION STENCIL (FOR HEAD WAVES ON FACES OF GRID CELLS)
-   STENCIL FROM VIDALE (1988 2D PAPER); CONDITIONS AND OTHER OPTIONS FROM HOLE
-   JAH 11/91 */
-float fdh2d(t1,t2,t3,ss0,s1,s2,s3)
-     float  t1,t2,t3,ss0,s1,s2,s3;
-     /* ss0 at newpoint; s1,t1 & s3,t3 adjacent; s2,t2 diagonal
-      */
-{
-  float x,slo;
-  double sqrt();
-  slo = .25*(ss0+s1+s2+s3);
-  x = 2.*slo*slo - (t3-t1)*(t3-t1);
-  if (x>=0.)  {
-    x = t2 + sqrt(x);
-    if ( (x<t1) || (x<t3) )  x = 1.e11;   /* ACAUSAL; ABORT */
-  }
-  else  x = 1.e11;   /* SQRT IMAGINARY; ABORT */
-  return(x);
-}
-
-/* 3D STENCIL FOR NEW FACE
-   STENCIL FROM VIDALE; CONDITIONS AND OTHER OPTIONS FROM HOLE
-   JAH 11/91 */
-float fdhnf(t1,t2,t3,t4,t5,ss0,s1)
-     float  t1,t2,t3,t4,t5,ss0,s1;
-     /* ss0 at newpoint; s1,t1 adjacent on old face;
-	t2,t4 beside t1 on old face and opposite each other;
-	t3,t5 beside t1 on old face and opposite each other
-	*/
-{
-  float x,slo;
-  double sqrt();
-  slo = .5*(ss0+s1);
-  x = slo*slo - .25*( (t4-t2)*(t4-t2) + (t5-t3)*(t5-t3) );
-  if (x>=0.)  {
-    x = t1 + sqrt(x);
-    if ( (x<t2) || (x<t3) || (x<t4) || (x<t5) )     /* ACAUSAL; ABORT */
-      x = 1.e11;
-  }
-  else  x = 1.e11;   /* SQRT IMAGINARY; ABORT */
-  return(x);
-}
-
-