/* 2D ACOUSTIC AND SH FINITE DIFFERENCE MODEL PROGRAM	*/
/* AUTHORS: John E. Vidale, Robert W. Clayton,		*/
/*  Richard Stead and Art Frankel, 			*/
/* Last modified: 3-17-87 to add gliding window option	*/
/*  		Caltech, Pasadena CA  91125  		*/

#include	<stdio.h>
#include	<sys/file.h>
#include	<math.h>

#define HETCHK 1.0e-10

/* DEFINITIONS FOR ELEMENTS OF TABLE
   (STENCILS ARE TABLE DRIVEN.)      */
#define ta0		na[pt[0]]
#define tac		na[pt[0]+1]
#define tai		na[pt[0]+2]
#define tao		na[pt[0]+3]
#define tb0		na[pt[0]+1]
#define tbp1		na[pt[0]+2]
#define tbm1		na[pt[0]+3]
#define tbu1		na[pt[0]+4]
#define tbd1		na[pt[0]+5]
#define tbp2		na[pt[0]+6]
#define tbm2		na[pt[0]+7]
#define tbu2		na[pt[0]+8]
#define tbd2		na[pt[0]+9]

/* DEFINITIONS FOR GLIDING ROUTINES */
#define bw1		ptmgm[0]
#define cnt1		ptmgm[1]
#define bw2		ptmgm[2]
#define cnt2		ptmgm[3]

/* ALLOW FOR LARGE NEGATIVE EXPONENTS TO BE 0.0 */
#define badexp(x)  (x< -50? 0.0 : exp((double)x))

/* GRID SPECIFICATIONS (SPACE, TIME AND BOUNDARY CONDITIONS) */
int	nx;		/* x-dimension of mesh */
int	nz;		/* z-dimension of mesh */
int	nt;		/* number of time steps */
float	h;		/* spatial mesh interval */
float	dt;		/* time digitization interval */
int	SHdisp	=1;     /* SHdisp = 1 for SH, = 0 for acoustic */
int	freesurf=1;	/* = 1 for freesurf, = 0 for absorbing top */
float	glide	=0;	/* compute only glide seconds after first arrival */
float	grease	=0;	/* compute grease seconds ahead of first arrival */

/* SOURCE SPECIFICATIONS */
int 	ns	=1;	/* number of shots */
int 	ds	=1;	/* mesh spacing between shots */
int 	s0;		/* initial shot position */
int	zs	=1;	/* shot depth */
int     rdsource=1;	/* rdsource=1, read source from input file */
int     style	=0;	/* 0=explosion,1=sin(theta),2=cos(theta) */
int	source	=7;	/* inner radius of source */
float   alpha=1./12.5;	/* width of gaussian*/
int	plane	=0;	/* plane=0, no plane wave.
			   plane=1, horizontal plane wave.
			   plane=2, vertical plane wave.
			   plane=3, plane wave at an angle.  */
float	plane_ang=0.0;	/* angle for plane=3 */
float	srcvel=0.0;	/* velocity to offset plane=3 */

/* POSITIONS AND OUTPUT SPECIFICATIONS OF	 */
/* TIME SLICES, RECEIVERS AND KIRCHHOFF ELEMENTS */
int	tslice	=0;	/* timestep interval between timeslices */
int	idt	=1;	/* record at every idt time steps */
int	zg	=1;	/* receiver depth */
int	ioff0	=0;	/* distance from shot to first receiver */
int	idoff	=1;	/* mesh increment between receivers */
int	noff	=0;	/* number of receivers */
int	zeropad	=1;	/* =1, fill in zeroes for receivers off the mesh */
int	allrecord=0;	/* allrecord=1 is equivalent to ioff0= -s0,
			   idoff=1, noff=nx (an entire horizontal line) */
int	nrcv    =0;	/* num. receivers for pointrecord */
int	*pts_x, *pts_z, *grid_pts;   /* positions of rcvs. for pointrecord */
int 	xvsp	=0;	/* x-position of a vertical line of recievers
			   xvsp<0, recievers positioned at source */
int	kz;		/* kirchhoff surface depth */
int	kdx	=1;	/* kir. out at every kdx grid points */


/* INPUT AND OUTPUT FILES */
char	source1[40];	/* source file */
char	tsfile[40];	/* timeslices */
char	tfile[40];	/* timing of first arrival (read in) */
char	seisfile[40];	/* seismograms (along the horizontal) */
char	zfile[40];	/* seismograms (along the vertical) */
char	pointfile[40];	/* seismograms (at arbitrary points) */
char	kirfile[40];	/* records for kirchhoff */
char	initfile[40];	/* input file:  first two timeslices */
char	finalfile[40];	/* output file containing final two timeslices
			   for restart purposes */

/* FLAGS FOR I/O REQUIREMENTS */
int	final	=0;	/* require output of final two timeslices */
int	back	=0;	/* initfile timeslices to be input in reverse
			   for reverse propagation */

/* FILE DESCRIPTORS */
int	zout	=0;	/* zfile file descriptor */
int	tsout	=0;	/* timeslice output file descriptor */
int	source1in;	/* source file descriptor */

/* FOR SOURCES */
float ring1[30],ring2[30],ring3[30],ring4[30];
float sinner1[300],sinner2[300],outer1[300],outer2[300];
float p3outer1[96],p3outer2[104];
float p3inner[450];

/* FOR TABLE */
unsigned short *ptable;
float ntable[400000];
float *a, *b;

/* OTHER INFO */
int 	it;		/* current timestep */
float	ct;		/* current time (it*dt) */
float	*p2, *p3;	/* current and previous timeslices */
float	*tmg;		/* stores first arrival time for grid */
int	*tmgm;		/* stores windows for a particular row */
float	*buf;		/* buffer for subroutine 'record' */
float	*o;		/* buffer for subroutine 'kirrecord' */
float   *pbuf;		/* buffer for subroutine 'pointrecord' */
int	restart	=0;	/* restart=1 for continuation of previous run */
			/* requires initfile from last run's finalfile */
int	homono	=0;	/* number of homogeneous grid points */
float   grid_off=0;	/* grid offset, !=0 to avoid underflow */

main(ac,av)
int ac; 
char **av;
{
	int is, s;				/* for sources */
	int in, kirout, seisout, pointout;	/* file descriptors */
	int isstart;				/* for restart */
	int i, len;
	float *pp2, *pp3, *temp;
	char pararg[16];

	setpar(ac,av);
	getpar("restart","d", &restart);
	getpar("nx",     "d", &nx);
	getpar("nz",     "d", &nz);
	getpar("nt",     "d", &nt);
	getpar("ns",     "d", &ns);
	getpar("ds",     "d", &ds);
	getpar("s0",     "d", &s0);
	getpar("zs",     "d", &zs);
	getpar("zg",     "d", &zg);
	getpar("h",      "f", &h);
	getpar("dt",     "f", &dt);
	getpar("kz",     "d", &kz);
	getpar("kdx",    "d", &kdx);
	getpar("alpha",  "f", &alpha);
	getpar("style",  "d", &style);
	getpar("idt",    "d", &idt);
	getpar("tslice", "d", &tslice);
	getpar("SHdisp", "d", &SHdisp);
	if(SHdisp)
		fprintf(stdout,"SH RUN!!\n");
	else
		fprintf(stdout,"ACOUSTIC RUN!!\n");
	getpar("plane",  "d", &plane);
	getpar("plane_ang","f",&plane_ang);
	getpar("srcvel","f",&srcvel);
	getpar("grid_off","f",&grid_off);
	getpar("allrecord", "d", &allrecord);
	getpar("freesurf", "d", &freesurf);
	getpar("ioff0",  "d", &ioff0);
	getpar("idoff",  "d", &idoff);
	getpar("noff",   "d", &noff);
	getpar("zeropad","d", &zeropad);
	getpar("source","d",&source);
	getpar("rdsource","d",&rdsource);

	/* BUFFERS FOR SOME I/O ROUTINES */
	buf = (float *) malloc(4*nx);
	o = (float *) malloc(8*nx);

	if(getpar("source1","s",source1))
	{
		if((source1in=open(source1,O_RDONLY,0644))<=1)
		{
			fprintf(stderr,"cannot open %s\n",source1);
		}
	}
	if(getpar("seisfile","s",seisfile))
	{
		if(restart)
		{
			int offset, ntrace, nshot;
			if((seisout= open(seisfile,O_RDONLY,0644)) <= 1)
			{
				fprintf(stderr,"cannot creat %s\n",seisfile);
				exit(-1);
			}
			offset= lseek(seisout,0L,2)/4;
			ntrace= offset/(nt/idt);
			nshot= ntrace/ noff;
			fprintf(stdout,
			"offset=%d ntrace=%d nshot=%d\n",
			offset,ntrace,nshot);
			isstart=nshot;
			lseek(seisout,nshot*noff*(nt/idt)*4,0);
		}
		else
		{
			if((seisout= open(seisfile,O_WRONLY|O_TRUNC|O_CREAT,0644)) <= 1)
			{
				fprintf(stderr,"cannot creat %s\n",seisfile);
				exit(-1);
			}
			isstart=0;
		}
	}
	if(getpar("kirfile","s",kirfile))
	{
		if(restart)
		{
			int offset, ntrace, nshot;
			if((kirout= open(kirfile,1)) == -1)
			{
				fprintf(stderr,"cannot creat %s\n",kirfile);
				exit(-1);
			}
			offset= lseek(kirout,0L,2)/8;
			ntrace= offset/(nt/idt);
			nshot= ntrace/ nx;
			fprintf(stderr,
			"offset=%d ntrace=%d nshot=%d\n",
			offset,ntrace,nshot);
			lseek(kirout,nshot*nx*(nt/idt)*8,0);
		}
		else
		{
			if((kirout= creat(kirfile,0644)) == -1)
			{
				fprintf(stderr,"cannot creat %s\n",kirfile);
				exit(-1);
			}
		}
	}
	if(getpar("pointfile","s",pointfile))
	{
		if((pointout= open(pointfile,O_WRONLY|O_TRUNC|O_CREAT,0644)) <= 1)
		{
			fprintf(stderr,"cannot creat %s\n",pointfile);
			exit(-1);
		}
		mstpar("nrcv","d",&nrcv);
		pts_x = grid_pts = (int *) malloc(4*nrcv);
		pts_z =            (int *) malloc(4*nrcv);
		pbuf  =          (float *) malloc(4*nrcv);
		sprintf(pararg,"vd(%d)",nrcv);
		mstpar("points_x",pararg,pts_x);
		mstpar("points_z",pararg,pts_z);
		for (i=0; i<nrcv; i++)
			grid_pts[i] += nx * (pts_z[i] + 1);
	}
	if(getpar("tsfile","s",tsfile))
	{
		if((tsout= open(tsfile,O_WRONLY|O_TRUNC|O_CREAT,0644)) <= 1)
		{
			fprintf(stderr,"cannot creat %s\n",tsfile);
			exit(-1);
		}
	}
	if(getpar("zfile","s",zfile))
	{
		xvsp= -999;
		getpar("xvsp","d",&xvsp);
		if((zout= open(zfile,O_WRONLY|O_TRUNC|O_CREAT,0644)) <= 1)
		{
			fprintf(stderr,"cannot creat %s\n",zfile);
			exit(-1);
		}
	}
	if(getpar("final","s",finalfile))
	{
		if((final= open(finalfile,O_WRONLY|O_TRUNC|O_CREAT,0644)) <= 1)
		{
			fprintf(stderr,"cannot creat %s\n",finalfile);
			exit(-1);
		}
	}
	getpar("back","d",&back);
        if(getpar("glide","f",&glide)) {
        	getpar("grease","f",&grease);
		tmg = (float *) malloc(4*nx*nz);
		tmgm = (int *) malloc(4*4*nz);
		for(i=0;i<nz;i++){
			tmgm[4*i] = -3;
			tmgm[4*i+1] = -1;
			tmgm[4*i+2] = -1;
			tmgm[4*i+3] = -1;
		}
		tmgm[4*zs-4] = tmgm[4*zs+4] = tmgm[4*zs] = s0 - 2;
		tmgm[4*zs-3] = tmgm[4*zs+5] = tmgm[4*zs+1] = 1;
		tmgm[4*zs-8] = -2;
		tmgm[4*zs+8] = -2;
		mstpar("tfile","s",tfile);
		if((in= open(tfile,O_RDONLY,0644)) <= 1) {
			fprintf(stderr,"cannot open tfile %s\n",tfile);
			exit(-1);
		}
		dread(in,tmg,nx*nz);
	}

	/* CONSTRUCT TABLE */
	len= nx*(nz+1)*4;
	a=  (float *) malloc(len);
	b=  (float *) malloc(len);
	ptable= (unsigned short *) malloc(len/2);
	if( a == NULL || b == NULL || ptable == NULL )
	{
		fprintf(stderr,"cannot allocate memory\n");
		exit(-1);
	}
	getab(a,b,nx,nz,ptable,ntable);
	free(a); 
	free(b);

	/* ALLOCATE SPACE FOR TIMESLICES */
	p2= (float *) malloc(len);
	p3= (float *) malloc(len);

	/* LOOP OVER SOURCES */
	for(is=isstart; is < ns; is++)
	{
		s= s0 + is*ds;
		pp2= p2; 
		pp3= p3;
		zap_fill(p2,nx*(nz+1),&grid_off);
		zap_fill(p3,nx*(nz+1),&grid_off);
		if (rdsource)
		{
			if(source>25){
				fprintf(stderr,"REDIMENSION SOURCE ARRAYS!");
				exit(-1);
			}
			dread(source1in,ring1,source+1);
			dread(source1in,ring2,source+2);
			dread(source1in,ring3,source+3);
			dread(source1in,ring4,source+4);
			mkring2(p2,source,s0,zs,nx,sinner1);
			addsource(p2,source,s0,zs,nx,ring1);
			mkring2(p2,source+1,s0,zs,nx,sinner2);
			addsource(p2,source+1,s0,zs,nx,ring2);
			for(i=0; i<8*(source+2); i++) outer1[i]=0.;
			for(i=0; i<8*(source+3); i++) outer2[i]=0.;
		}
		if (rdsource==0)
			init(p3,p2,nx,s,zs,SHdisp,alpha);
		if(seisout) record(p3,nx,zg,s,seisout);
		if(kirout) kirrecord(p3,o,kz,kirout);
		if(seisout && idt==1) record(p2,nx,zg,s,seisout);
		if(kirout && idt==1) kirrecord(p2,o,kz,kirout);
		if(pointout) pointrecord(p3,nrcv,grid_pts,pointout);
		if(pointout && idt==1)pointrecord(p2,nrcv,grid_pts,pointout);
		if(zout) zrecord(p3,nx,nz,s,zout);
		if(zout && idt==1) zrecord(p2,nx,nz,s,zout);
		if(tslice) slice(p3,nx,nz,tsout);

		/* LOOP OVER TIME */
		for(it=2; it< nt; it++)
		{
			ct = it * dt + grease;
			tstep(pp2,pp3,nx,nz,ptable,ntable);
			if(seisout && it%idt == 0) record(pp3,nx,zg,s,seisout);
			if(kirout && it%idt == 0) kirrecord(pp3,o,kz,kirout);
			if(pointout && it%idt == 0) pointrecord(pp3,nrcv,grid_pts,pointout);
			if(zout && it%idt == 0) zrecord(pp3,nx,nz,s,zout);
			if(tslice && it%tslice == 0) slice(pp3,nx,nz,tsout);
			if(final && it== nt-1)
			{
				write(final,&pp3[nx],nx*nz*4);
				write(final,&pp2[nx],nx*nz*4);
			}

			temp= pp3; 
			pp3= pp2; 
			pp2= temp;
		}
		close(source1in);
		if(getpar("source1","s",source1))
		{
			if((source1in=open(source1,O_RDONLY,0644))<=1)
			{
				fprintf(stderr,"cannot open %s\n",source1);
			}
		}
	}
	endpar();
}

tstep(p2,p3,nx,nz,ptable,ntable)
int nx ,nz;
float *p2, *p3, *ntable;
unsigned short *ptable;
{
	int iz;
	float srcoff= -1;
	register float *pp2, *pp3, *ptmg;
	register int *ptmgm;
	register unsigned short *pptable;
	int offset;
	int cnt;

	offset= 2*nx + 2;
	cnt= nx - 4;
	if(glide>0.0) srcoff = ct - 2.*glide;
	pp2 = p2 +offset;
	pp3 = p3 +offset;
	ptmg = tmg + offset;
	ptmgm = tmgm + 8;
	pptable = ptable + offset;

	if(rdsource && srcoff<0) savit(p3,source,s0,zs,nx,p3inner);

	/* EXECUTE FOURTH-ORDER EQUATION */
	for(iz=2; iz < nz-1; iz++)
	{
		if(glide!=0){
			if( ptmgm[0] >= -2){
				re_tmg(ptmgm,ptmg);
				if( ptmgm[0] >= 0 )
					inner4(pp2+bw1,pp3+bw1,ntable,nx,cnt1,pptable+bw1);
				if( ptmgm[2] >= 0 )
					inner4(pp2+bw2,pp3+bw2,ntable,nx,cnt2,pptable+bw2);
			}
		}
		else {inner4(pp2,pp3,ntable,nx,cnt,pptable);}

		pp2 += nx; 
		pp3 += nx;
		pptable += nx;
		ptmg += nx;
		ptmgm += 4;
	}

	/* DO SECOND-ORDER EDGE CALCULATION */

	ord2(p3,p2,ptable,ntable,nx+2,1,nx-4);
	ord2(p3,p2,ptable,ntable,nx*(nz-1)+2,1,nx-4);
	ord2(p3,p2,ptable,ntable,nx+1,nx,nz-1);
	ord2(p3,p2,ptable,ntable,2*nx-2,nx,nz-1);

	if(rdsource && srcoff<0)
	{
		mkring2(p2,source+2,s0,zs,nx,p3outer1);
		mkring2(p2,source+3,s0,zs,nx,p3outer2);
		bkring(p2,source,s0,zs,nx,sinner1);
		bkring(p2,source+1,s0,zs,nx,sinner2);
		bkring(p2,source+2,s0,zs,nx,outer1);
		bkring(p2,source+3,s0,zs,nx,outer2);
		recall(p3,source,s0,zs,nx,p3inner);
		offset= s0-source-1+(zs-source-1)*nx;
		pp2= p2+offset;
		pp3= p3+offset;
		pptable = ptable + offset;
		cnt= 3+ 2*source;
		for (iz=zs-source-1; iz<= zs+source+1; iz++)
		{
			inner4(pp2,pp3,ntable,nx,cnt,pptable);
			pp2 +=nx; 
			pp3 +=nx;
			pptable += nx;
		}
		dread(source1in,ring1,source+1);
		dread(source1in,ring2,source+2);
		dread(source1in,ring3,source+3);
		dread(source1in,ring4,source+4);
		mkring2(p3,source,s0,zs,nx,sinner1);
		mkring2(p3,source+1,s0,zs,nx,sinner2);
		mkring(p3,source+2,s0,zs,nx,ring3,outer1);
		mkring(p3,source+3,s0,zs,nx,ring4,outer2);
		addsource(p3,source,s0,zs,nx,ring1);
		addsource(p3,source+1,s0,zs,nx,ring2);
		bkring(p2,source+2,s0,zs,nx,p3outer1);
		bkring(p2,source+3,s0,zs,nx,p3outer2);
	}

	switch (plane) {
	case 0:			/* no plane-wave */
		/* FREE SURFACE */
		if(freesurf) {
			if(SHdisp==0) zap(p3,nx); 
			if(SHdisp==1) zip(p3,nx); 
		}
		else
			absorb(p3,p2,ptable,ntable,1,1,nx,nx-2);
		absorb(p3,p2,ptable,ntable,nx*nz+1,1,-nx,nx-2);
		absorb(p3,p2,ptable,ntable,0,nx,1,nz+1);
		absorb(p3,p2,ptable,ntable,nx-1,nx,-1,nz+1);
		break;
	case 1:			/* horz. plane-wave */
		/* FREE SURFACE */
		if(freesurf) {
			if(SHdisp==0) zap(p3,nx); 
			if(SHdisp==1) zip(p3,nx); 
		}
		else
			absorb(p3,p2,ptable,ntable,1,1,nx,nx-2);
		absorb(p3,p2,ptable,ntable,nx*nz+1,1,-nx,nx-2);
		reflect(p3,0,nx,1,nz+1);
		reflect(p3,nx-1,nx,-1,nz+1);
		break;
	case 2:			/* vert. plane-wave */
		reflect(p3,1,1,nx,nx-2);
		reflect(p3,nx*nz+1,1,-nx,nx-2);
		absorb(p3,p2,ptable,ntable,0,nx,1,nz+1);
		absorb(p3,p2,ptable,ntable,nx-1,nx,-1,nz+1);
		break;
	case 3:			/* angle plane-wave */
		absorb(p3,p2,ptable,ntable,1,1,nx,nx-2);
		absorb(p3,p2,ptable,ntable,nx*nz+1,1,-nx,nx-2);
		absorb(p3,p2,ptable,ntable,0,nx,1,nz+1);
		absorb(p3,p2,ptable,ntable,nx-1,nx,-1,nz+1);
		break;
	}
}

reflect(p3,offset,tang,norm,cnt)
register float *p3;
register int tang, norm, cnt;
int offset;
{
	p3 += offset;
	while(cnt--)
	{
		*p3= p3[norm];
		p3 += tang;
	}
}
ord2(p3,p2,pt,na,offset,tang,cnt)
register float *p3, *p2, *na;
unsigned short *pt;
register int tang, cnt;
int offset;
{
	p3 += offset;
	p2 += offset;
	pt += offset;
	while(cnt--)
	{
		inner(p2,p3,na,nx,1,pt);
		p3 += tang;
		p2 += tang;
		pt += tang;
	}
}
absorb(p3,p2,pt,na,offset,tang,norm,cnt)
register float *p3, *p2, *na;
unsigned short *pt;
register int tang, norm, cnt;
int offset;
{
	p3 += offset;
	p2 += offset;
	pt += offset;
	while(cnt--)
	{
		*p3= p2[norm] -ta0*(p3[norm] -p2[0]);
		p3 += tang;
		p2 += tang;
		pt += tang;
	}
}

slice(p1,nx,nz,tsout)
register float *p1;
register int nx, nz, tsout;
{
	p1 += nx ;
	write(tsout,p1,nx*nz*4);
}

record(p,nx,zg,s,seisout)
float *p;
int nx, zg, s, seisout;
{
	register int n, i, index;
	register float *pp;
	if(allrecord)
	{
		write(seisout,&p[zg*nx],nx*4);
		return;
	}
	n=0;
	pp= p+ zg*nx;
	for(i=0; i< noff; i++)
	{
		index= (ioff0+s) + i*idoff;
		if(index < 0 || index >= nx)
		{
			if(zeropad) buf[n++]= 0.0;
		}
		else
		{
			buf[n++]= pp[index];
		}
	}
	write(seisout,buf,4*n);
}

kirrecord(p,o,kz,kirout)
float *p;
float *o;
int kz, kirout;
{
	register int i;
	register float *pp;
	float *out;
	int knx;

	pp= p+ kz*nx;
	out = o;
	knx = nx / kdx;
	for(i=0; i< knx; i++)
	{
		*out++ = *pp;
		pp += kdx;
	}
	write(kirout,o,4*knx);
	pp= p+ kz*nx;
	out = o;
	for(i=0; i< knx; i++)
	{
		*out++ = (pp[nx] - pp[-nx])/(2*h);
		pp += kdx;
	}
	write(kirout,o,4*knx);
}

pointrecord(p,nrcv,grid_pts,pointout)
float *p;
int nrcv, pointout, *grid_pts;
{
	register int i;

	for(i=0; i< nrcv; i++)
		pbuf[i]= p[grid_pts[i]];
	write(pointout,pbuf,4*nrcv);
}

zrecord(p,nx,nz,s,zout)
float *p;
int nx, nz, s, zout;
{
	register int i;
	register float *pp, *junk;
	junk = (float *) malloc(4*nz);
	if(xvsp < 0)	pp= p+ nx+ s;
	else		pp= p+ nx+ xvsp;
	for(i=0; i< nz; i++)
	{
		junk[i]= *pp;
		pp += nx;
	}
	write(zout,junk,4*nz);
	free(junk);
}

init(p0,p1,nx,s,zs,SHdisp,alpha)
float *p0, *p1, alpha;
int nx, s, zs, SHdisp;
{
	float x,zr, zi,x2,amp;
	float zrsq, zisq, ampzr, ampzi, ampz, ampx;
	float check, c_pa, s_pa, rp, rpsq, taper[4], fac;
	float s0b, zsb, *pp;
	int ix0, index;
	int in, ix, iz, j, need_taper;
  
        if(getpar("init","s",initfile))
	   {
		if((in= open(initfile,O_RDONLY,0644)) <= 1)
		   {
			fprintf(stderr,"cannot open %s\n",initfile);
			exit(-1);
                   }
                 if(back)
		    {
			dread(in,&p0[nx],nx*nz);
			dread(in,&p1[nx],nx*nz);
                    }
                 else
		    {
			dread(in,&p1[nx],nx*nz);
			dread(in,&p0[nx],nx*nz);
                    }
                  return;
            }
	if(plane) goto planewave;
	for( iz= 0; iz <= nz; iz++)
	{
		zr= iz - zs;
		zi= iz + zs;
		zrsq= zr*zr;
		zisq= zi*zi;
		ampzr= badexp(-alpha*zrsq);
		ampzi= badexp(-alpha*zisq);
		if(freesurf==1) ampz= (SHdisp? ampzr+ampzi : ampzr-ampzi );
		else ampz= ampzr ;

		for( ix= 0; ix < nx; ix++)
		{
			x= s - ix;
			x2= x*x;
			ampx= badexp( -alpha*x2 );
			amp= ampz*ampx;
			p0[ix + iz*nx]= p1[ix + iz*nx]= amp;
		}
	}
	return;
planewave:
	if(plane==1){
	for( iz= 0; iz <= nz; iz++)
	{
		zr= (float) (iz-zs);
		zrsq= zr*zr;
		amp= badexp(-alpha*zrsq);

		for( ix= 0; ix < nx; ix++)
			p0[ix + iz*nx]= p1[ix + iz*nx]= amp;
	}
	}
	if(plane==2){
	for( ix= 0; ix < nx; ix++)
	{
		zi= (float) (s0-ix);
		zisq= zi*zi;
		amp= badexp(-alpha*zisq);

		for( iz= 0; iz <= nz; iz++)
			p0[ix + iz*nx]= p1[ix + iz*nx]= amp;
	}
	}
	if (plane==3){
	   need_taper = 1;
	   check = (float)(((int)(plane_ang / 180.0 + 100.5)) - 100);
	   plane_ang -= check * 180.0;
	   plane_ang *= 3.1415926 / 180.0;
	   c_pa = cos(plane_ang);
	   s_pa = sin(plane_ang);
	   pp = p0;
           for( iz= 0; iz <= nz; iz++) {
              for( ix = 0; ix < nx; ix++) {
                 rp= fabs(s_pa * (s0 - ix) + c_pa * (zs - iz));
                 rpsq= rp*rp;
                 amp= badexp(-alpha*rpsq);
                 *pp++ += amp;
              }
	   }
	   s0b = s0 - s_pa * (srcvel * dt);
	   zsb = zs - c_pa * (srcvel * dt);
	   pp = p1;
           for( iz= 0; iz <= nz; iz++) {
              for( ix = 0; ix < nx; ix++) {
                 rp= fabs(s_pa * (s0b - ix) + c_pa * (zsb - iz));
                 rpsq= rp*rp;
                 amp= badexp(-alpha*rpsq);
                 *pp++ += amp;
              }
	   }
	}
	if (need_taper) {
	   taper[0] = 0.1; taper[1] = 0.4; taper[2] = 0.6; taper[3] = 0.9;
	   for (j=0;j<4;j++) {
	      fac = taper[j];
	      /* Top */
	      ix0 = j*nx;
	      index = (j+1)*nx;
	      for(ix=ix0;ix<index;ix++) { p0[ix] *= fac; p1[ix] *= fac; }
	      /* Bottom */
	      ix0 = (nz-j)*nx;
	      index = (nz+1-j)*nx;
	      for(ix=ix0;ix<index;ix++) { p0[ix] *= fac; p1[ix] *= fac; }
	      /* Left */
	      ix0 = j;
	      index = j + nz*nx;
	      for(ix=ix0;ix<index;ix += nx) { p0[ix] *= fac; p1[ix] *= fac; }
	      /* Right */
	      ix0 = nx-1-j;
	      index = nx-1-j + nz*nx;
	      for(ix=ix0;ix<index;ix += nx) { p0[ix] *= fac; p1[ix] *= fac; }
	   }
	}
	return;
}

#define M	1	/* modulus */
#define D	2	/* density */
#define V	4	/* velocity */
#define R	8	/* impedence */
#define NOTFOUND	0
#define NUMBER		1
#define MATRIX		2

getab(a,b,nx,nz,ptable,na)
unsigned short *ptable;
float *a, *b, *na;
int nx, nz;
{
	int pt[2], key, len, cnt, iz, ix, i, j, k, l, test;
	float a0,ac, ai, ao, b0, bp1, bm1, bp2, bm2, bu1, bd1, bu2, bd2;
	float *pa, *pb;
	float vmax, vsq, vsqmax, fac;
	float mval, dval, vval, rval;
	int check[2];
	int   mkey, dkey, vkey, rkey;
	int   mf,   df,   vf,   rf;
	double sqrt();

	/* SEE WHAT INPUT MEDIUM PARAMETERS ARE SPECIFIED */
	dkey= vkey= rkey= mkey= key=0;
	if((mkey=getparm("modulus",  &mf,&mval)))  key |= M;
	if((vkey=getparm("velocity", &vf,&vval)))  key |= V;
	if((dkey=getparm("density",  &df,&dval)))  key |= D;
	if((rkey=getparm("impedance",&rf,&rval)))  key |= R;

	/* IF SHdisp, CONSTRUCT THE 'b' (=MODULUS) AND 'a' (=1.0/DENSITY) */
	/* OTHERWISE, CONSTRUCT THE 'a' (=MODULUS) AND 'b' (=1.0/DENSITY) */
	len= nx*nz;
	switch(key)
	{
	case 0: 
	case M: 
	case D: 
	case V: 
	case R:
	case M|D|V: 
	case M|D|R: 
	case D|V|R: 
	case M|D|V|R:
		fprintf(stderr,"TWO MEDIUM PARAMETER FILES, PLEASE!\n");
		exit(-1);
		break;

	case M|D:	/* medium specified by modulus and density */
		if(mkey==MATRIX) dread(mf,&a[nx],nx*nz);
		if(dkey==MATRIX) dread(df,&b[nx],nx*nz);
		if(mkey==MATRIX) close(mf);
		if(dkey==MATRIX) close(df);
		vsqmax= 0.0;
		for(pa=a+nx,pb=b+nx,cnt=len; cnt--; pa++, pb++)
		{
			if(mkey==MATRIX) mval= *pa;
			if(dkey==MATRIX) dval= *pb;
			vsq= mval/dval;
			if(vsq > vsqmax) vsqmax= vsq;
			if(SHdisp==1){
				*pa= 1.0/dval;
				*pb= mval;
			}
			if(SHdisp==0){
				*pb=1.0/dval;
				*pa=mval;
			}
		}
		vmax= sqrt(vsqmax);
		break;
	case M|V:	/* medium specified by modulus and velocity */
		if(mkey==MATRIX) dread(mf,&a[nx],nx*nz);
		if(vkey==MATRIX) dread(vf,&b[nx],nx*nz);
		if(mkey==MATRIX) close(mf);
		if(vkey==MATRIX) close(vf);
		vmax= 0.0;
		for(pa=a+nx,pb=b+nx,cnt=len; cnt--; pa++, pb++)
		{
			if(mkey==MATRIX) mval= *pa;
			if(vkey==MATRIX) vval= *pb;
			dval= mval/(vval*vval);
			if(vval > vmax) vmax= vval;
			if(SHdisp==1){
				*pa= 1.0/dval;
				*pb= mval;
			}
			if(SHdisp==0){
				*pb=1.0/dval;
				*pa=mval;
			}
		}
		break;
	case M|R:	/* medium specified by modulus and impedence */
		if(mkey==MATRIX) dread(mf,&a[nx],nx*nz);
		if(rkey==MATRIX) dread(rf,&b[nx],nx*nz);
		if(mkey==MATRIX) close(mf);
		if(rkey==MATRIX) close(rf);
		vsqmax= 0.0;
		for(pa=a+nx,pb=b+nx,cnt=len; cnt--; pa++, pb++)
		{
			if(mkey==MATRIX) mval= *pa;
			if(rkey==MATRIX) rval= *pb;
			dval= rval*rval/mval;
			vsq= mval/dval;
			if(vsq > vsqmax) vsqmax= vsq;
			if(SHdisp==1){
				*pa= 1.0/dval;
				*pb= mval;
			}
			if(SHdisp==0){
				*pb=1.0/dval;
				*pa=mval;
			}
		}
		vmax= sqrt(vsqmax);
		break;
	case D|V:	/* medium specified by density and velocity */
		if(vkey==MATRIX) dread(vf,&a[nx],nx*nz);
		if(dkey==MATRIX) dread(df,&b[nx],nx*nz);
		if(vkey==MATRIX) close(vf);
		if(dkey==MATRIX) close(df);
		vmax= 0.0;
		for(pa=a+nx,pb=b+nx,cnt=len; cnt--; pa++, pb++)
		{
			if(vkey==MATRIX) vval= *pa;
			if(dkey==MATRIX) dval= *pb;
			mval= dval*vval*vval;
			if(vval > vmax) vmax= vval;
			if(SHdisp==1){
				*pa= 1.0/dval;
				*pb= mval;
			}
			if(SHdisp==0){
				*pb=1.0/dval;
				*pa=mval;
			}
		}
		break;
	case D|R:	/* medium specified by density and impedence */
		if(rkey==MATRIX) dread(rf,&a[nx],nx*nz);
		if(dkey==MATRIX) dread(df,&b[nx],nx*nz);
		if(rkey==MATRIX) close(rf);
		if(dkey==MATRIX) close(df);
		vsqmax= 0.0;
		for(pa=a+nx,pb=b+nx,cnt=len; cnt--; pa++, pb++)
		{
			if(rkey==MATRIX) rval= *pa;
			if(dkey==MATRIX) dval= *pb;
			mval= rval*rval/dval;
			vsq= mval/dval;
			if(vsq > vsqmax) vsqmax= vsq;
			if(SHdisp==1){
				*pa= 1.0/dval;
				*pb= mval;
			}
			if(SHdisp==0){
				*pb=1.0/dval;
				*pa=mval;
			}
		}
		vmax= sqrt(vsqmax);
		break;
	case V|R:	/* medium specified by velocity and impedence */
		if(vkey==MATRIX) dread(vf,&a[nx],nx*nz);
		if(rkey==MATRIX) dread(rf,&b[nx],nx*nz);
		if(vkey==MATRIX) close(vf);
		if(rkey==MATRIX) close(rf);
		vmax= 0.0;
		for(pa=a+nx,pb=b+nx,cnt=len; cnt--; pa++, pb++)
		{
			if(vkey==MATRIX) vval= *pa;
			if(rkey==MATRIX) rval= *pb;
			dval= rval/vval;
			mval= rval*vval;
			if(vval > vmax) vmax= vval;
			if(SHdisp==1){
				*pa= 1.0/dval;
				*pb= mval;
			}
			if(SHdisp==0){
				*pb=1.0/dval;
				*pa=mval;
			}
		}
		break;
	}

	fac= sqrt(8./3.)*dt*vmax/h;
	if(fac >= 1.0)
	{
		fprintf(stderr,"**instability** fac=%8.5f  (>1.0) --stop\n",fac);
		exit();
	}
	fprintf(stdout,"**stability** fac=%8.5f  (<1.0) --go\n",fac);
	/* MULTIPLY THE 'a' FILE BY THE FINITE DIFFERENCE CONSTANTS */
	fac= 0.5*dt*dt/(h*h);
	for(pa=a+nx,cnt=len; cnt--; pa++) *pa= fac* *pa;

	/* INSTALL THE TOP ROW */
	for(pa=a,cnt=nx; cnt--; pa++) *pa= pa[nx];
	for(pb=b,cnt=nx; cnt--; pb++) *pb= pb[nx];

	/* SCAN THE 'b' FILE FOR REGIONS OF CONSTANT VALUE (OVER THE FINITE */
	/* DIFFERENCE STENCIL. ENCODE THESE REGIONS BY SETTING THE 'a' FILE */
	/* VALUE NEGATIVE (AND MULTIPLYING BY 2*a). */
	for(iz=2; iz <  (nz-1);    iz++)
		for(ix=2; ix < (nx-2); ix++)
		{
			len= iz*nx + ix;
			pb= b + len;
			if( fabs((pb[0] - pb[-1] )  /pb[0])>HETCHK) continue;
			if( fabs((pb[0] - pb[1]  )  /pb[0])>HETCHK) continue;
			if( fabs((pb[0] - pb[-nx])  /pb[0])>HETCHK) continue;
			if( fabs((pb[0] - pb[nx] )  /pb[0])>HETCHK) continue;
			if( fabs((pb[0] - pb[-2*nx])/pb[0])>HETCHK) continue;
			if( fabs((pb[0] - pb[-2])   /pb[0])>HETCHK) continue;
			if( fabs((pb[0] - pb[2])    /pb[0])>HETCHK) continue;
			if( fabs((pb[0] - pb[2*nx]) /pb[0])>HETCHK) continue;
			pa= a + len;
			pa[0]= -2.0*pa[0]*pb[0];
			homono+=1;
		}

	/* RECHECK EDGES FOR SECOND ORDER HOMOGENEOUS CASE */
	for( ix = 1 ; ix < nx-1 ; ix++)
	{
		check[0] = nx+ix;
		check[1] = nx*(nz-1)+ix;
		for(iz=0;iz<2;iz++)
		{
			if(a[check[iz]]<0)continue;
			if(fabs((b[check[iz]]-b[check[iz] -1])/b[check[iz]])>HETCHK)continue;
			if(fabs((b[check[iz]]-b[check[iz] +1])/b[check[iz]])>HETCHK)continue;
			if(fabs((b[check[iz]]-b[check[iz]-nx])/b[check[iz]])>HETCHK)continue;
			if(fabs((b[check[iz]]-b[check[iz]+nx])/b[check[iz]])>HETCHK)continue;
			a[check[iz]]*= (-2.0)*b[check[iz]];
			homono+=1;
		}
	}
	for(iz=2;iz<nz-1;iz++)
	{
		check[0] = nx*iz+1;
		check[1] = nx*(iz+1)-2;
		for(ix=0;ix<2;ix++)
		{
			if(a[check[ix]]<0)continue;
			if(fabs((b[check[ix]]-b[check[ix] -1])/b[check[ix]])>HETCHK)continue;
			if(fabs((b[check[ix]]-b[check[ix] +1])/b[check[ix]])>HETCHK)continue;
			if(fabs((b[check[ix]]-b[check[ix]-nx])/b[check[ix]])>HETCHK)continue;
			if(fabs((b[check[ix]]-b[check[ix]+nx])/b[check[ix]])>HETCHK)continue;
			a[check[ix]]*= (-2.0)*b[check[ix]];
			homono+=1;
		}
	}
	/* SET THE SIDE EDGES TO (1-v*dt/h)/(1+v*dt/h) FOR THE ABSORBING */
	/* BOUNDARY CONDITIONS. */
	for(ix=0; ix<nx; ix++)
	{
		len= ix;
		fac= sqrt(2.0*a[len]*b[len]);
		a[len]= (1.0-fac)/(1.0+fac);
		len= nx*nz +ix;
		fac= sqrt(2.0*a[len]*b[len]);
		a[len]= (1.0-fac)/(1.0+fac);
	}
	for(iz=0; iz<=nz; iz++)
	{
		len= iz*nx;
		fac= sqrt(2.0*a[len]*b[len]);
		a[len]= (1.0-fac)/(1.0+fac);
		len= nx*iz + (nx-1);
		fac= sqrt(2.0*a[len]*b[len]);
		a[len]= (1.0-fac)/(1.0+fac);
	}
	fprintf(stdout,"homo= %d, hetero= %d \n",homono, nx*nz-homono);
	/* ENCODE MEDIA PARAMETERS IN TABLE */
	fprintf(stdout,"encode media parameters\n");
	len=0;
	/* CHECK FOURTH ORDER PORTION */
	for(l=2;l<nz-1;l++)
		for(k=2;k<nx-2;k++){
			i = nx*l + k;
			pa = a + i;
			pb = b + i;
			test=0;
			a0= pa[0];
			if(a0<0){
				ac=5.*a0+2.;
				ai=(-a0*4./3.);
				ao=a0*1./12.;
			}
			else{
				bp1=  (pb[    1]+pb[0])*a0*4./3.;
				bm1=  (pb[   -1]+pb[0])*a0*4./3.;
				bd1=  (pb[   nx]+pb[0])*a0*4./3.;
				bu1=  (pb[  -nx]+pb[0])*a0*4./3.;
				bp2= -(pb[    2]+pb[0])*a0/12.;
				bm2= -(pb[   -2]+pb[0])*a0/12.;
				bd2= -(pb[ 2*nx]+pb[0])*a0/12.;
				bu2= -(pb[-2*nx]+pb[0])*a0/12.;
				b0 = 2.-(bp1+bm1+bd1+bu1+bp2+bm2+bd2+bu2);
			}
			for(j=len-1;j>=0;j--){
				pt[0]=j*10;
				if(a0<0){
					if(ac==tac&&ai==tai&&ao==tao
					    &&a0==ta0) 
					{
						test=1;
						ptable[i]=pt[0];
					}
				}
				else{
					if(bp1==tbp1 && b0==tb0 &&
					    bm1==tbm1 && bm1==tbm1 &&
					    bp2==tbp2 && bm2==tbm2 &&
					    bd1==tbd1 && bu1==tbu1 &&
					    bd2==tbd2 && bu2==tbu2 &&
					    a0 ==ta0)
					{
						test=1;
						ptable[i]=pt[0];
					}
				}
			}
			if(test==0){
				pt[0]=len*10;
				ta0=a0;
				if(a0<0){
					tac=ac; 
					tai=ai; 
					tao=ao;
				}
				else{
					tb0=b0; 
					tbp1=bp1; 
					tbm1=bm1;
					tbp2=bp2; 
					tbm2=bm2; 
					tbd1=bd1;
					tbu1=bu1; 
					tbd2=bd2; 
					tbu2=bu2;
				}
				ptable[i]=pt[0];
				len+=1;
				if(len>6000){
					fprintf(stderr,"redimension ntable or debug program!\n");
					fprintf(stderr,"fourth order z=%d x=%d\n",l,k);
					exit();
				}
			}
		}
	/* TABULATE SECOND ORDER PORTION */
	fprintf(stdout,"tabulate second order portion\n");
	for(k=0;k<2*nx+2*nz-10;k++)
	{
		if(k<nx-3) i = k + nx + 1;
		if(k>=nx-3 && k<nx+nz-4) i = nx*((k-(nx-3))+2) - 2;
		if(k>=nx+nz-4 && k<2*nx+nz-7) i=(k-(nx+nz-4))+ nx*(nz-1)+1;
		if(k>= 2*nx+nz-7 ) i = nx*((k-(2*nx+nz-7))+2) + 1;
		pa = a + i;
		pb = b + i;
		test=0;
		a0= pa[0];
		ac=2.+4.*a0;
		bp1= (pb[  1] + pb[0])*a0;
		bm1= (pb[ -1] + pb[0])*a0;
		bd1= (pb[ nx] + pb[0])*a0;
		bu1= (pb[-nx] + pb[0])*a0;
		b0= 2. - (bp1 + bm1 + bd1 + bu1);
		for(j=len-1;j>=0;j--){
			pt[0]=j*10;
			if(a0<0){	
				if(a0==ta0 && ac==tac)
				{
					test=1;
					ptable[i]=pt[0];
				}
			}
			else{
				if(b0==tb0 && bp1==tbp1  && bm1==tbm1 &&
				    bd1==tbd1  && bu1==tbu1 && a0==ta0)
				{
					test=1;
					ptable[i]=pt[0];
				}
			}
		}
		if(test==0){
			pt[0]=len*10;
			ta0=a0;
			if(a0<0){
				tac=ac;
			}
			else{
				tb0=b0;
				tbp1=bp1;
				tbm1=bm1;
				tbd1=bd1;
				tbu1=bu1;
			}
			ptable[i]=pt[0];
			len+=1;
			if(len>4000){
				fprintf(stderr,"redimension ntable or debug program!\n");
					fprintf(stderr,"second order k=%d\n",k);
				exit();
			}
		}
	}
	/* TABULATE EDGES */
	fprintf(stdout,"tabulate edges\n");
	for(k=0;k<2*nx+2*nz-2;k++)
	{
		if(k<nx-1) i = k;
		if(k>=nx-1&&k<nx+nz-1) i = nx*((k-(nx-1))+1) - 1;
		if(k>=nx+nz-1&&k<2*nx+nz-2) i = (k-(nx+nz-1))+nx*nz + 1;
		if(k>=2*nx+nz-2) i = nx*((k-(2*nx+nz-2))+1);
		pa = a + i;
		pb = b + i;
		test=0;
		a0= pa[0];
		b0= pb[0];
		for(j=len-1;j>=0;j--){
			pt[0]=j*10;
			if(a0==ta0 && b0==tb0) 
			{
				test=1;
				ptable[i]=pt[0];
			}
		}
		if(test==0){
			pt[0]=len*10;
			ta0=a0; 
			tb0=b0;
			ptable[i]=pt[0];
			len+=1;
			if(len>4000){
				fprintf(stderr,"redimension ntable or debug program!\n");
					fprintf(stderr,"edge\n");
				exit();
			}
		}
	}
	fprintf(stdout,"done with table formation, %d types of media! \n",len);
}

getparm(name,file,val)
char *name;
float *val;
int *file;
{
	char junk[128];
	double atof();
	if(getpar(name,"s",junk) == 0) return(NOTFOUND);
	if(isnumber(junk))
	{
		*val= atof(junk);
		return(NUMBER);
	}

	*file= chkopen(junk,name);
	return(MATRIX);
}

chkopen(filename,filetype)
char *filename, *filetype;
{
	int f;
	if( (f=open(filename,O_RDONLY,0644)) <= 1)
	{
		fprintf(stderr,"cannot open %s=%s\n",filetype,filename);
		exit(-1);
	}
	return(f);
}

dread(file,junk,n)
register int file, n;
float *junk;
{
	if(read(file,junk,4*n)==4*n) return;
	fprintf(stderr,"dread: bad read\n");
	exit(-1);
}

zip(v,cnt)
register float *v;
register int cnt;
{
	int nx;
	nx=cnt;
	while(cnt--) {
		*v = v[2*nx];
		v++;
	}
}

zap(v,cnt)
register float *v;
register int cnt;
{
	while(cnt--) *v++ = 0.0;
}

zap_fill(v,cnt,value)
register int *v;
register int cnt;
int *value;
{
   register int x;

   x = *value;
   for(;cnt--;) *v++ = x;
}

isnumber(str)
char *str;
{
	if(*str == '\0') return(0);
	do
	    {
		if(*str=='+' || *str=='-' || *str=='.' || *str=='e') continue;
		if(*str <= '9' && *str >= '0') continue;
		return(0);
	} 
	while(*(++str));
	return(1);
}

inner(p2,p3,na,nx,cnt,pt)
unsigned short *pt;
register float *p2, *p3, *na;
register int nx, cnt;
{

	while(cnt--)
	{
		if(ta0 < 0.0)
		{
			*p3 = tac*p2[0] -p3[0]
			    - ta0*(p2[1] + p2[-1] + p2[nx] + p2[-nx]);
		}
		else
		{
			*p3 = tbp1*p2[1]  + tbm1*p2[-1]
			    + tbd1*p2[nx] + tbu1*p2[-nx]
			    + tb0*p2[0]   - p3[0];
		}

		p2++; 
		p3++; 
		pt++;
	}
}

inner4(p2,p3,na,nx,cnt,pt)
register float *p2, *p3, *na;
register unsigned short *pt;
register int nx, cnt;
{

	while(cnt--)
	{
		if(ta0 < 0.0)
		{
			*p3 = tac*p2[0] -p3[0]
			    + tai*( p2[1] + p2[-1] + p2[nx] + p2[-nx] )
				+ tao*( p2[2] + p2[-2] + p2[2*nx] + p2[-2*nx] );
		}
		else
		{
			*p3 = tbp1*p2[1   ]  + tbm1*p2[-1   ]
			    + tbp2*p2[2   ]  + tbm2*p2[-2   ]
			    + tbd1*p2[nx  ]  + tbu1*p2[-nx  ]
			    + tbd2*p2[2*nx]  + tbu2*p2[-2*nx]
			    + tb0 * p2[0] - p3[0];

		}

		p2++; 
		p3++; 
		pt++;
	}
}
re_tmg(ptmgm,ptmg)
float *ptmg;
int *ptmgm;
{
	int workL, workR;
	/* TWO HALVES ABOVE AND BELOW */
	if(ptmgm[-2] >= 0 && ptmgm[6] >= 0 ){
		workR = fillinL(&ptmgm[2],ptmg);
		if(workR){
			fillinR(&ptmgm[2],ptmg);
			workL = fillinL(&ptmgm[0],ptmg);
			if(workL){
				fillinR(&ptmgm[0],ptmg);
				return;
			}
			/* LEFT HALF MOVED OFF GRID */
			else{
				ptmgm[0] = ptmgm[2];
				ptmgm[1] = ptmgm[3];
				ptmgm[2] = -1;
				ptmgm[3] = -1;
			}
		}
		/* RIGHT HALF MOVED OFF GRID */
		else{
			ptmgm[2] = -1;
			ptmgm[3] = -1;
			workL = fillinL(&ptmgm[0],ptmg);
			if(workL){
				fillinR(&ptmgm[0],ptmg);
				return;
			}
			/* LEFT AND RIGHT HALVES MOVED OFF GRID */
			else{
				ptmgm[0] = -4;
				ptmgm[1] = -1;
			}
		}
	}
	/* SOMETHING ABOVE AND BELOW */
	else if(ptmgm[-4] >= 0 && ptmgm[4] >= 0 ){
		/* NOT BRACKETED BY TWO HALVES */
		workL = fillinL(&ptmgm[0],ptmg);
		if(workL){
			fillinx(&ptmgm[0],ptmg);
			/* SEE WHETHER IT SPLITS */
			split(ptmgm,ptmg);
		}
		else{
			ptmgm[0] = -4;
			ptmgm[1] = ptmgm[2] = ptmgm[3] = -1;
		}
		return;
	}
	/* ROW ABOVE OR BELOW NOT HIT */
	else {
		start(ptmgm,ptmg);
		return;
	}
}
fillinL(ptmgm,ptmg)
float *ptmg;
int *ptmgm;
{
	int i, x0, x1, junk;
	x0 = ptmgm[-4];
	x1 = ptmgm[4];
	if( x0 > x1 ) { junk = x1; x1 = x0; x0 = junk;}
	if( x0 != 0) x0 -= 1;
	if( x1 != nx - 5) x1 += 1;
	for( i = x0 ; i <= x1 ; i++ )
		if(ct > ptmg[i] && ct < ptmg[i] + glide){
			ptmgm[0] = i;
			return(1);
		}
	return(0);
}
fillinR(ptmgm,ptmg)
float *ptmg;
int *ptmgm;
{
	int i, x0, x1, junk;
	x0 = ptmgm[-4] + ptmgm[-3] - 1;
	x1 = ptmgm[4] + ptmgm[5] - 1;
	if( x0 > x1 ) { junk = x1; x1 = x0; x0 = junk;}
	if( x0 != 0) x0 -= 1;
	if( x1 != nx - 5) x1 += 1;
	for( i = x1 ; i >= x0 ; i-- )
		if(ct > ptmg[i] && ct < ptmg[i] + glide){
			ptmgm[1] = i - ptmgm[0] + 1;
			return;
		}
	fprintf(stderr,"FILLINR FAILED\n CALL THE REPAIRMAN: J. VIDALE\n");
}
fillinx(ptmgm,ptmg)
float *ptmg;
int *ptmgm;
{
	int i, x0= -1, x1= -1, junk;
	if(ptmgm[-2] >= 0) x0 = ptmgm[-2] + ptmgm[-1] -1;
	else { x0 = ptmgm[-4] + ptmgm[-3] - 1;}
	if(ptmgm[6] >= 0) x1 = ptmgm[6] + ptmgm[7] -1;
	else { x1 = ptmgm[4] + ptmgm[5] - 1;}
	if( x0 > x1 ){junk = x0; x0 = x1; x1 = junk;}
	if( x0 != 0) x0 -= 1;
	if( x1 != nx - 5) x1 += 1;
	for( i = x1 ; i >= x0 ; i-- )
		if(ct > ptmg[i] && ct < ptmg[i] + glide){
			ptmgm[1] = i - ptmgm[0] + 1;
			ptmgm[2] = -1;
			ptmgm[3] = -1;
			return;
		}
	fprintf(stderr,"FILLINX FAILED\n CALL THE REPAIRMAN: J. VIDALE\n");
}
start(ptmgm,ptmg)
float *ptmg;
int *ptmgm;
{
	int i, beg, end, wbeg, wend;
	/* DOWNWARD WAVE */
	if(ptmgm[-4] >= 0){
		beg = ptmgm[-4];
		if( beg != 0) beg -= 1;
		end = ptmgm[-4] + ptmgm[-3] - 1;
		if(ptmgm[-2] >= 0) end = ptmgm[-2] + ptmgm[-1] - 1;
		if( end != nx - 5) end += 1;
	}
	/* UPWARD WAVE */
	else{
		beg = ptmgm[4];
		if( beg != 0) beg -= 1;
		end = ptmgm[4] + ptmgm[5] - 1;
		if(ptmgm[6] >= 0) end = ptmgm[6] + ptmgm[7] - 1;
		if( end != nx - 5) end += 1;
	}
	wend = -1;
	wbeg = nx - 4;
	for(i=beg;i<=end;i++)
		if(ct > ptmg[i] && ct < ptmg[i] + glide){
			if(wend<i) wend = i;
			if(wbeg>i) wbeg = i;
		}
	/* ROW IS ACTIVE */
	if(wend >= 0){
		ptmgm[0] = wbeg;
		ptmgm[1] = wend - wbeg + 1;
		if(ptmgm[-4] == -3) ptmgm[-4] = -2;
		if(ptmgm[ 4] == -3) ptmgm[ 4] = -2;
	}
	/* ROW IS NOT ACTIVE */
	else{
		ptmgm[0] = -2;
		ptmgm[1] = -1;
		ptmgm[2] = -1;
		ptmgm[3] = -1;
		if(ptmgm[-4] == -4) ptmgm[ 0] = -4;
		if(ptmgm[ 4] == -4) ptmgm[ 0] = -4;
	}
	return(1);
}
split(ptmgm,ptmg)
float *ptmg;
int *ptmgm;
{
	int i, nend, nbeg, oend, obeg;
	obeg = ptmgm[0];
	oend = ptmgm[0] + ptmgm[1] - 1;
	i = obeg;
	while(ct > ptmg[i] && ct < ptmg[i] + glide && i <= oend) i++;
	if( i == oend + 1 ) return(1);
	nend = i - 1;
	while((ct <= ptmg[i] || ct >= ptmg[i] + glide) && i <= oend) i++;
	if( i == oend + 1 ) fprintf(stderr,"SCREW-UP IN SPLIT: DEBUG!!\n");
	nbeg = i;
	ptmgm[1] = nend - obeg + 1; 
	ptmgm[2] = nbeg;
	ptmgm[3] = oend - nbeg + 1; 
	return(1);
}