/* 2D ELASTIC FINITE DIFFERENCE MODELING PROGRAM */
/* AUTHORS: Robert W. Clayton and John E. Vidale */ 
/*                    and Art Frankel            */ 
/* revised and tested by Richard Stead, 4-18-86  */
/*                                    , 10-22-86 */
/*                                    , 4-29-87  */
/*          all at:  Caltech, Pasadena CA  91125 */
/* modfified for restart capability by Doug Dreger 2-24-92*/

#include   <stdio.h>
#include   <sys/file.h>
#ifdef USAGE
/* compile with -DUSAGE if system is bsd4.2 */
#include <sys/time.h>
#include <sys/resource.h>
#endif
int REPORT = 50;                  /* usage reporting interval (timesteps) */

#define    HETCHK       1.0e-6    /* min. diff. for material comparison */
int NTSIZE = 1000;                /* starting number of table entries */

/* table entries for material parameters, each grid point */
#define    ta(x)      na[pt[x]+0]
#define    tb(x)      na[pt[x]+1]
#define    tc(x)      na[pt[x]+2]
#define    td(x)      na[pt[x]+3]
#define    te(x)      na[pt[x]+4]

/* table entries for fourth order homogeneous */
#define    tcc        na[pt[0]+3]
#define    tb16       na[pt[0]+4]
#define    tbb        na[pt[0]+5]
#define    ta16       na[pt[0]+6]
#define    taa        na[pt[0]+7]
#define    tcd        na[pt[0]+8]
#define    tcd10      na[pt[0]+9]

/* table entries for fourth order heterogeneous */
#define    tlu2m2     na[pt[0]+3]
#define    tlu2m1     na[pt[0]+4]
#define    tlu2       na[pt[0]+5]
#define    tlu2p1     na[pt[0]+6]
#define    tlu2p2     na[pt[0]+7]
#define    tlu1m2     na[pt[0]+8]
#define    tlu1m1     na[pt[0]+9]
#define    tlu1       na[pt[0]+10]
#define    tlu1p1     na[pt[0]+11]
#define    tlu1p2     na[pt[0]+12]
#define    tl00m2     na[pt[0]+13]
#define    tl00m1     na[pt[0]+14]
#define    tl00       na[pt[0]+15]
#define    tl00p1     na[pt[0]+16]
#define    tl00p2     na[pt[0]+17]
#define    tld1m2     na[pt[0]+18]
#define    tld1m1     na[pt[0]+19]
#define    tld1       na[pt[0]+20]
#define    tld1p1     na[pt[0]+21]
#define    tld1p2     na[pt[0]+22]
#define    tld2m2     na[pt[0]+23]
#define    tld2m1     na[pt[0]+24]
#define    tld2       na[pt[0]+25]
#define    tld2p1     na[pt[0]+26]
#define    tld2p2     na[pt[0]+27]
#define    tmu2m2     na[pt[0]+28]
#define    tmu2m1     na[pt[0]+29]
#define    tmu2       na[pt[0]+30]
#define    tmu2p1     na[pt[0]+31]
#define    tmu2p2     na[pt[0]+32]
#define    tmu1m2     na[pt[0]+33]
#define    tmu1m1     na[pt[0]+34]
#define    tmu1       na[pt[0]+35]
#define    tmu1p1     na[pt[0]+36]
#define    tmu1p2     na[pt[0]+37]
#define    tm00m2     na[pt[0]+38]
#define    tm00m1     na[pt[0]+39]
#define    tm00       na[pt[0]+40]
#define    tm00p1     na[pt[0]+41]
#define    tm00p2     na[pt[0]+42]
#define    tmd1m2     na[pt[0]+43]
#define    tmd1m1     na[pt[0]+44]
#define    tmd1       na[pt[0]+45]
#define    tmd1p1     na[pt[0]+46]
#define    tmd1p2     na[pt[0]+47]
#define    tmd2m2     na[pt[0]+48]
#define    tmd2m1     na[pt[0]+49]
#define    tmd2       na[pt[0]+50]
#define    tmd2p1     na[pt[0]+51]
#define    tmd2p2     na[pt[0]+52]


int   homono   =0;   /* number of homogeneous grid points */
int   ns       =1;   /* number of shots */
int   ds       =1;   /* mesh spacing between shots */
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   idt      =1;   /* record at every idt time steps */
int   s0;            /* initial shot position */
int   zs       =1;   /* shot depth */
int   zg       =1;   /* receiver depth */
int   style    =0;   /* 1= compress. explo., 0 and 2=double couple */
                     /*             3= shear explosion */
int   ww      =10;   /* width of space-specified source */
int   tslice   =0;   /* number of timesteps between timeslice output */
int   ostyle   =0;   /* ostyle=1 for divergence and curl */
                     /* ostyle=0 for u (x-disp) and w (z-disp) */
int   freesurf =0;   /* 1=freesurf, 0=absorb at top */
int   rigid    =0;   /* 1 makes all boundaries rigid */
int   ioff0    =0;   /* distance from shot to first receiver */
int   idoff    =1;   /* mesh increment between receivers */
int   noff     =0;   /* number of receivers */
int   nang     =0;   /* # of angles for rad output */
float dang;          /* angular distance bewteen angles for rad output */
float radius;        /* radius to angular output */
float ang0;          /* starting angle for angular output */
float grid_off=0.;   /* value to offset grid by in zap_fill */
int   pressure =0;   /* pressure instead of u & w disp. is output */
int   kdx      =1;   /* kirchoff output at every kdx grid points */
int   kz;            /* kirchoff depth */

/* following for analytical source */
int   rdsource =1;   /* 0 for quick, 1 for source from ntde code */
float pvel;          /* alpha for quick init plane wave */
float svel;          /* beta for quick init plane wave */
int   ntend;         /* how many time points source stays on */
int   source   =9;   /* inside diameter of source box */
                     /* note that arrays are hard wired to this value */
		     /* If a larger value of source is desired then */
		     /* redimension or dynamically allocate */

/* filenames */
char name[40];
int  tsfile=0;    /* filename on if set to 1 */
int  seisfile=1;  
int  kirfile=0;  

int   it;            /* timestep counter */

struct disp
{
   float u;
   float w;
};
#define WAVE struct disp      /*variable type for wave field */
#define NW   2                /*number of floats in WAVE */
struct disp *p1, *p2;

/* following for free surface */
float *fsg;
float fsu[6], fsw[6];

/* following for table */
/* float ntable [53*NTSIZE]; */
float *ntable;
float *pa, *pb, *pc, *buf;
int   *ptable;

/* following for analytical source */
char source1[40];
int source1in;
float ring1[160],ring2[176],ring3[192],ring4[208];
struct disp sinner1[160],sinner2[176],outer1[192],outer2[208];
struct disp p3outer1[192],p3outer2[208],p3inner[900];

/* following for record ouput */
struct disp *comhold;
float *comheld, *comkhold;

/* following for restart */
int restart=0, save=0, getabc(), table_size;
void save_run(), get_savefile(), chkpar(), error();
int tsout=0, seisout=0, kirout=0; /*Moved from local in main to global*/

main(ac,av)
int ac; 
char **av;
{
   struct disp *pp1, *pp2, *temp;
   int i, is, s, len, knx;
#ifdef USAGE
   struct rusage use1, use2;
   float useout;
#else
   int time0, time1, time2;
#endif

#ifndef USAGE
   time(&time0);
#endif

/* Restart Block - Dreger */
fprintf(stderr,"Welcome to %s ...\n",av[0]);
   setpar(ac,av);
    
      getpar("restart",  "d", &restart);
        if (restart == 1) {           /* restarting saved run */
	get_savefile(&s,&is,&it,&len);
	fprintf(stderr,"s=%d is=%d it=%d len=%d\n",s,is,it,len);
	pp1 = p1;
	pp2 = p2;
#ifdef USAGE
	getrusage(RUSAGE_SELF,&use1);
	fprintf(stderr,
	"status given at initialization and every %d timesteps\n",REPORT);
	fprintf(stderr,"\nAll times in seconds\n");
	fprintf(stderr,
	"Timestep   System    User   Cum. Sys.  Cum. User  Faults  Swaps\n");          
	fprintf(stderr,
	"  %5d                      %5d.     %5d.     %5d  %5d\n",1,
	use1.ru_stime.tv_sec,use1.ru_utime.tv_sec,use1.ru_majflt,use1.ru_nswap);
#else
	time(&time1);
	fprintf(stderr,"initialization took %d seconds\n",time1-time0);
	time2=time1;
#endif
	endpar();
	goto start_point;
   }


/* general parameters */
 /*  setpar(ac,av); */
   mstpar("name","s",name);
   getpar("save","d",&save);
   getpar("REPORT", "d", &REPORT);
   getpar("NTSIZE", "d", &NTSIZE);
   mstpar("nx",     "d", &nx);
   mstpar("nz",     "d", &nz);
   mstpar("nt",     "d", &nt);
   getpar("ns",     "d", &ns);
   getpar("ds",     "d", &ds);
   getpar("s0",     "d", &s0);
   getpar("zs",     "d", &zs);
   zs +=1;
   getpar("zg",     "d", &zg);
   getpar("style",  "d", &style);
   getpar("ww",     "d", &ww);
   mstpar("h",      "f", &h);
   mstpar("dt",     "f", &dt);
   getpar("idt",    "d", &idt);
   getpar("tslice", "d", &tslice);
   getpar("ostyle", "d", &ostyle);
   getpar("freesurf", "d", &freesurf);
   getpar("pressure", "d", &pressure);
   getpar("grid_off", "f", &grid_off);
   getpar("rigid",  "d", &rigid);

   /* for seismograms */
   getpar("ioff0",  "d", &ioff0);
   getpar("idoff",  "d", &idoff);
   getpar("noff",   "d", &noff);
   getpar("nang",   "d", &nang);
   getpar("dang",   "f", &dang);
   getpar("radius", "f", &radius);
   getpar("ang0",   "f", &ang0);

   /* for kirchoff */
   getpar("kdx",    "d", &kdx);
   getpar("kz",     "d", &kz);

   /* following for analytic source */
   getpar("source", "d", &source);
   getpar("pvel",   "f", &pvel);
   getpar("svel",   "f", &svel);
   if(source>10){
      fprintf(stderr,"source too large: redimension arrays");
      exit();
   }
   getpar("ntend",  "d", &ntend);
   getpar("rdsource", "d", &rdsource);
   if(ntend==0) ntend=nt;
   if(getpar("source1","s",source1)) {
      if((source1in=open(source1,O_RDONLY,0664)) <= 1)
         fprintf(stderr,"cannot open %s\n",source1);
      fprintf(stderr,"File %s assigned descriptor %d\n",source1,source1in);
   }

   /* create output files */
   getpar("seisfile","d",&seisfile);
   getpar("kirfile","d",&kirfile);
   getpar("tsfile","d",&tsfile);
   if(seisfile)
      seisout= chkcreat(name,"seis");
   if(kirfile)
      kirout= chkcreat(name,"kir");
   if(tsfile)
      tsout= chkcreat(name,"ts");

   /* tabulate material properties and stencil coefficients */
   len= nx*(nz+1)*sizeof(float);
   pa=     (float *) malloc(len);
   pb=     (float *) malloc(len);
   pc=     (float *) malloc(len);
   ptable= (int   *) malloc(nx*(nz+1)*sizeof(int));
   buf=    (float *) malloc(8*nx);
   fsg=    (float *) malloc(4*nx);
   if(pa == NULL || pb == NULL || pc == NULL || ptable == NULL)
   { 
      fprintf(stderr,"cannot allocate a, b, c memory\n"); 
      exit(-1); 
   }
   table_size=getabc(pa,pb,pc,nx,nz,ptable,&ntable);
   fprintf(stderr,"table_size=%d\n",table_size);
   free(pa);
   free(pb);
   free(pc);

   /* allocate output arrays */
   if (seisout) {
	comhold= (struct disp *) malloc(noff*8);
   	comheld= (float *) malloc(noff*4);
   }
   if (kirout) {
   	knx = nx/kdx;
   	comkhold= (float *) malloc(4*knx);
   }

   /* allocate displacement arrays */
   p1= (struct disp *) malloc(2*len);
   p2= (struct disp *) malloc(2*len);
   if(p1==NULL||p2==NULL)
   { 
      fprintf(stderr,"cannot allocate p1, p2 memory\n"); 
      exit(-1); 
   }


   /* tabulate free surface parameters */
   fsfill(fsg,fsu,fsw,ptable,ntable,nx);

   endpar();

   /* enter source(s) */
   for(is=0; is < ns; is++)
   {
      s= s0 + is*ds;
      pp1= p1;
      pp2= p2;
      zap_fill(p1,2*nx*(nz+1),&grid_off);
      zap_fill(p2,2*nx*(nz+1),&grid_off);
      if(rdsource)
      {
            dread(source1in,ring1,16*source);
            dread(source1in,ring2,16*(source+1));
            dread(source1in,ring3,16*(source+2));
            dread(source1in,ring4,16*(source+3));
            mkring2(p2,source,s0,zs,nx,sinner1);
            addeqsource(p2,source,s0,zs,nx,ring1);
            mkring2(p2,source+1,s0,zs,nx,sinner2);
            addeqsource(p2,source+1,s0,zs,nx,ring2);

         for(i=0; i<8*(source+2); i++) 
         {
            outer1[i].u =0.;
            outer1[i].w =0.;
         }
         for(i=0; i<8*(source+3); i++) 
         {
            outer2[i].u =0.;
            outer2[i].w =0.;
         }
      }


      /* generate source */
      if(rdsource==0){
         quickinit(pp1,pp2,nx,s0,zs,style,ww,pvel,svel);
      }

      /* output any points necessary from source */
      if(seisout) {
         if(noff)record(pp1,nx,zg,s,seisout);
         if(nang)radrec(pp1,nx,zs,s,seisout);
      }
      if(seisout && idt==1){
         if(noff)record(pp2,nx,zg,s,seisout);
         if(nang)radrec(pp1,nx,zs,s,seisout);
      }
      if(kirout) kirrecord(pp1,kirout);
      if(kirout && idt==1) kirrecord(pp2,kirout);
      if(tslice)    slice(pp1,nx,nz,tsout,buf);
      if(tslice==1) slice(pp2,nx,nz,tsout,buf);

      /* report progress */

#ifdef USAGE
      /* Only for bsd4.2 systems: */
      getrusage(RUSAGE_SELF,&use1);
      fprintf(stderr,"status given at initialization and every %d timesteps\n",REPORT);
      fprintf(stderr,"\nAll times in seconds\n");
      fprintf(stderr,"Timestep   System    User   Cum. Sys.  Cum. User  Faults  Swaps\n");
      fprintf(stderr,"  %5d                      %5d.     %5d.     %5d  %5d\n",1,use1.ru_stime.tv_sec,use1.ru_utime.tv_sec,use1.ru_majflt,use1.ru_nswap);
#else
      /* Otherwise: */
      time(&time1);
      fprintf(stderr,"initialization took %d seconds\n",time1-time0);
      time2=time1;
#endif

      /* begin time iteration */
      for(it=2; it< nt; it++)
      {

         /* report progress */

#ifdef USAGE
         /* Only for bsd4.2 systems: */
         if (it%(REPORT)==(REPORT-1)) getrusage(RUSAGE_SELF,&use2);
         if (it%(REPORT)==0) {
            getrusage(RUSAGE_SELF,&use1);
            fprintf(stderr,"  %5d   ",it);
 	    useout = use1.ru_stime.tv_sec-use2.ru_stime.tv_sec+(use1.ru_stime.tv_usec-use2.ru_stime.tv_usec)/1000000.;
            fprintf(stderr,"%7.3f  ",useout);
 	    useout = use1.ru_utime.tv_sec-use2.ru_utime.tv_sec+(use1.ru_utime.tv_usec-use2.ru_utime.tv_usec)/1000000.;
            fprintf(stderr,"%7.3f   ",useout);
            fprintf(stderr,"%5d.     %5d.     %5d  %5d\n",use1.ru_stime.tv_sec,use1.ru_utime.tv_sec,use1.ru_majflt,use1.ru_nswap);
         }
#else
	 /* Otherwise: */
         time1=time2;
         time(&time2);
         if(it%(REPORT)==0) fprintf(stderr,"time step %d took %d seconds\n",it,time2-time1);
#endif


	 /* Before call of tstep, pp2 is the previous timestep, pp1 is
	    the timestep before pp2.  pp1 will be overwritten in tstep. */
         tstep(pp1,pp2,ptable,ntable,nx,nz);

	 /* Output required values of current timestep (pp1) */
         if(seisout && it%idt == 0){
            if(noff)record(pp1,nx,zg,s,seisout);
            if(nang)radrec(pp1,nx,zs,s,seisout);
         }
         if(kirout && it%idt == 0) kirrecord(pp1,kirout);
         if(tslice && it%tslice == 0) slice(pp1,nx,nz,tsout,buf);

	 /* Output runstate if save is requested */
		  if (save && it%REPORT==0) save_run(is,it,len,pp1,pp2);
start_point:

	 /* make pp2 current timestep, saving previous in pp1 */
         temp= pp1; 
         pp1= pp2; 
         pp2= temp;
      }
      if ((int)lseek(source1in,(long)0,0) < 0) {
	       fprintf(stderr,"failure to seek in source file - exit\n");
			exit(-1);
      }
   }
}

tstep(p1,p2,ptable,ntable,nx,nz)

/*
	   Subroutine TSTEP is responsible for generating the next
	timestep from the previous two timesteps, p1 and p2,   for
	the grid and material parameters represented by    ptable,
	ntable, nx and nz.  It takes care of entering and   saving
	the source, applying boundary conditions, and calling  the
	finite difference subroutines.
*/

int *ptable;
float *ntable;
int nx ,nz;
struct disp *p1, *p2;
{
   int ix, iz;
   register struct disp *pp1, *pp2;
   register int *pptable;
   int offset, cnt;

   if(it<ntend && rdsource) savit(p1,source,s0,zs,nx,p3inner);
   /* use fourth order for most grid points */
   for(iz=2; iz < nz-1; iz++)
   {
      offset=2+iz*nx;
      pp1    = p1     + offset;
      pp2    = p2     + offset;
      pptable= ptable + offset;
      inner4(pp1,pp2,pptable,ntable,nx-4,nx);
   }
   if(rigid==0 || rigid==1){
      /* apply second order code to edges */

      offset = nx + 1 ;
      pp1     = p1     + offset;
      pp2     = p2     + offset;
      pptable = ptable + offset;
      inner(pp1,pp2,pptable,ntable,nx-2);

      offset = nx*(nz-1)+1;
      pp1     = p1     + offset;
      pp2     = p2     + offset;
      pptable = ptable + offset;
      inner(pp1,pp2,pptable,ntable,nx-2);
      for ( ix = 2 ; ix < nz-1 ; ix++ )
      {
         offset = nx*ix+1;
         pp1     = p1     + offset;
         pp2     = p2     + offset;
         pptable = ptable + offset;
         inner(pp1,pp2,pptable,ntable,1);
         offset = nx*(ix+1)-2;
         pp1     = p1     + offset;
         pp2     = p2     + offset;
         pptable = ptable + offset;
         inner(pp1,pp2,pptable,ntable,1);
      }
   }

   /*  for analytical source */
   if(it<ntend && rdsource)
   {
      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(p1,source,s0,zs,nx,p3inner);
      offset= s0-source-1+(zs-source-1)*nx;
      pp1= p1+offset;
      pp2= p2+offset;
      pptable= ptable+offset;
      cnt= 3+ 2*source;
      for (iz=zs-source-1; iz<= zs+source+1; iz++)
      {
         inner4(pp1,pp2,pptable,ntable,cnt,nx);
         pp1 +=nx; 
         pp2 +=nx; 
         pptable +=nx;
      }
      if(rdsource==1){
         dread(source1in,ring1,16*source);
         dread(source1in,ring2,16*(source+1));
         dread(source1in,ring3,16*(source+2));
         dread(source1in,ring4,16*(source+3));
         mkring2(p1,source,s0,zs,nx,sinner1);
         mkring2(p1,source+1,s0,zs,nx,sinner2);
         mkeqring(p1,source+2,s0,zs,nx,ring3,outer1);
         mkeqring(p1,source+3,s0,zs,nx,ring4,outer2);
         addeqsource(p1,source,s0,zs,nx,ring1);
         addeqsource(p1,source+1,s0,zs,nx,ring2);
      }
      bkring(p2,source+2,s0,zs,nx,p3outer1);
      bkring(p2,source+3,s0,zs,nx,p3outer2);
   }

   /* apply boundary conditions here */
   if(rigid==0){
      if(freesurf==0) absorb(p1,p2,ptable,ntable,0,1,nx,nx);
      absorb(p1,p2,ptable,ntable, nx*nz+1, 1,-nx,  nx-2);
      absorb(p1,p2,ptable,ntable, nx     ,nx,  1,  nz);
      absorb(p1,p2,ptable,ntable, 2*nx-1 ,nx, -1,  nz);
   }
   if(freesurf!=0) frob(p1,p2,nx,fsg,fsu,fsw);
}

inner(p1,p2,pt,na,cnt)

/*
	   Subroutine INNER applies a second order stencil to the
	edges of the finite difference grid, accounting for   the
	points between the interior (fourth order) and the  boun-
	dary.
*/

struct disp *p1, *p2;
int *pt;
float *na;
int cnt;
{
   float a, b, cc, cd;
   float piece0, piece1, piece2, piece3;
   double fabs();

   p1--; 
   p2--;
   pt--;
   while(cnt--)
   {
      p1++;
      p2++;
      pt++;

      /* homogeneous case */
      if(tc(0)<0){
         a = -ta(0)/tc(0);
         b = -(tb(0)+2.*ta(0))/tc(0);
         cd= (b-a)/4.;
         cc= 2.0*(1.0-(a+b));
         p1->u=  cc*p2[0].u - p1[0].u
             +b*(p2[1].u +p2[-1].u)
            +a*(p2[nx].u+p2[-nx].u)
               +cd*(p2[nx+1].w+p2[-nx-1].w-p2[nx-1].w-p2[-nx+1].w);

         p1->w=  cc*p2[0].w - p1[0].w
             +a*(p2[1].w +p2[-1].w)
            +b*(p2[nx].w+p2[-nx].w)
               +cd*(p2[nx+1].u +p2[-nx-1].u-p2[nx-1].u-p2[-nx+1].u);
      }
      else{

      /* inhomogeneous case */
         piece0=(tb( 1)+tb(0))*(p2[ 1].u-p2[0].u)/2.+
             tb( 1)*(p2[ nx+1].w-p2[1-nx].w)/4.+
             (ta( 1)+ta(0))*(p2[ 1].u-p2[0].u);

         piece1=(tb(-1)+tb(0))*(p2[-1].u-p2[0].u)/2.+
             tb(-1)*(p2[-nx-1].w-p2[nx-1].w)/4.+
             (ta(-1)+ta(0))*(p2[-1].u-p2[0].u);

         piece2=(ta( nx)+ta(0))*(p2[ nx].u-p2[0].u)/2.+
	     ta( nx)*(p2[ nx+1].w-p2[nx-1].w)/4.;

         piece3=(ta(-nx)+ta(0))*(p2[-nx].u-p2[0].u)/2.+
	     ta(-nx)*(p2[-nx-1].w-p2[1-nx].w)/4.;

         p1->u=(piece0+piece2+piece1+piece3)/
            tc(0)+2.*p2[0].u-p1[0].u;

         piece0=(tb( nx)+tb(0))*(p2[ nx].w-p2[0].w)/2.+
             tb( nx)*(p2[ nx+1].u-p2[nx-1].u)/4.+
             (ta( nx)+ta(0))*(p2[ nx].w-p2[0].w);

         piece1=(tb(-nx)+tb(0))*(p2[-nx].w-p2[0].w)/2.+
             tb(-nx)*(p2[-nx-1].u-p2[1-nx].u)/4.+
             (ta(-nx)+ta(0))*(p2[-nx].w-p2[0].w);

         piece2=(ta( 1)+ta(0))*(p2[ 1].w-p2[0].w)/2.+
	     ta( 1)*(p2[ nx+1].u-p2[1-nx].u)/4.;

         piece3=(ta(-1)+ta(0))*(p2[-1].w-p2[0].w)/2.+
	     ta(-1)*(p2[-nx-1].u-p2[nx-1].u)/4.;

         p1->w=(piece0+piece2+piece1+piece3)/
            tc(0)+2.*p2[0].w-p1[0].w;
      }
   }
}

slice(p1,nx,nz,tsout,buf)

/*
	   Subroutine SLICE writes a timeslice to tsout every
	tslice timesteps.  The options are:  ostyle=0,  ouput
	is displacements; ostyle=1, output is divergence  and
	curl.
*/

register float *buf;
register struct disp *p1;
register int nx, nz, tsout;
{
   int ix,iz;
   if(ostyle==0)write(tsout,&p1[nx].u,nx*nz*8);

   /* following for div and curl */
   if(ostyle==1){
      for(ix=0; ix<2*nx; ix++) buf[ix]= 0.0;
      write(tsout,buf,8*nx);
      for(iz=1; iz<nz-1; iz++)
      {
         for(ix=1; ix<nx-1; ix++)
         {

	    /* divergence */
            buf[2*ix] =(p1[1+ix+nx*iz].u-p1[-1+ix+nx*iz].u
                +p1[nx+ix+nx*iz].w-p1[-nx+ix+nx*iz].w);

	    /* curl */
            buf[2*ix+1]= -(p1[nx+ix+nx*iz].u-p1[-nx+ix+nx*iz].u
                -p1[1+ix+nx*iz].w+p1[-1+ix+nx*iz].w);
         }
         write(tsout,buf,8*nx);
      }
      for(ix=0; ix<2*nx; ix++) buf[ix]= 0.0;
      write(tsout,buf,8*nx);
   }
}

chkcreat(name1,type)

/*
	   Function CHKCREAT attempts to create filename, returning
	the resulting file descriptor if successful, or printing an
	error and aborting the program if not.
*/

char *name1, *type;
{
   char filename[100];
   int f;
   sprintf(filename,"%s.%s",name1,type);
   if( (f=open(filename,O_WRONLY | O_CREAT | O_TRUNC,0664)) < 0)
   {
      fprintf(stderr,"cannot creat %s\n",filename);
      exit(-1);
   }
   fprintf(stderr,"File %s assigned file descriptor %d\n",filename,f);
   return(f);
}

dread(file,buf,n)

/*
	   Subroutine DREAD reads n float values from file file
	into memory starting at pointer buf, printing an  error
	and aborting if unsuccessful.
*/

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

zap_fill(v,cnt,value)

/*
	   Subroutine ZAP_FILL fills cnt spaces in memory with *value,
	starting at pointer v.
*/

register int *v;
register int cnt;
int *value;
{
   register int x;

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

zap(v,cnt)

/*
	   Subroutine ZAP zeroes cnt spaces in memory of type float,
	starting at pointer v.
*/

register float *v;
register int cnt;
{
   for(;cnt--;) *v++ = 0.0;
}

quickinit(p0,p1,nx,s,zs,style,ww,pvel,svel)
struct disp *p0, *p1;
float pvel, svel;
int nx, s, zs, style,ww;
{
   float alpha = 1./25.;
   float x, zr, x1, x2, amp, ampmax, theta;
   float zrsq, ampzr, ampz, ampx, i;
   int ix, iz;
   double atan2(),sin(), cos(), exp(), fabs();

   ampmax= 0.0;
   if(style==5||style==6){
      for(ix=0;ix<50;ix++){
         i = ix - 25;
         x = i * exp( -i*i *alpha);
         if(style==5){
            i = ix - 25 - pvel*dt/h;
            x1 = i * exp( -i*i*alpha);
            fprintf(stderr,"style 5 x= %g x1= %g \n",x,x1);
            for(iz=0;iz<nz+1;iz++){
               p0[iz*nx+ix].u=x;
               p1[iz*nx+ix].u=x1;
            }
         }
         if(style==6){
            i = ix - 25 - svel*dt/h;
            x1 = i * exp( -i*i*alpha);
            fprintf(stderr," style 6 x= %g x1= %g \n",x,x1);
            for(iz=0;iz<nz+1;iz++){
               p0[iz*nx+ix].w=x;
               p1[iz*nx+ix].w=x1;
            }
         }
      }
      return;
   }
   for( iz= zs-ww; iz <= zs +ww; iz++)
   {
      zr= (float) (iz-zs);
      zrsq= zr*zr;
      ampzr= exp(-alpha*zrsq);
      ampz= ampzr;

      for( ix= s-ww; ix <= s + ww; ix++)
      {
         x= (float) (s-ix);
         x2= x*x;
         ampx= exp( -alpha*x2 );
         amp= ampz*ampx;
         if(fabs(amp) > ampmax) ampmax= fabs(amp);
         theta=atan2((double)(iz-zs),(double)(ix-s));
         if(style==0){
            p0[ix + iz*nx].u= p1[ix + iz*nx].u= -amp*sin(theta);
            p0[ix + iz*nx].w= p1[ix + iz*nx].w= amp*cos(theta);
         }
         if(style==1){
            p0[ix + iz*nx].u= p1[ix + iz*nx].u= amp*cos(theta);
            p0[ix + iz*nx].w= p1[ix + iz*nx].w= -amp*sin(theta);
         }
         if(style==2){
            p0[ix + iz*nx].u= p1[ix + iz*nx].u= amp*cos(theta);
            p0[ix + iz*nx].w= p1[ix + iz*nx].w= amp*sin(theta);
         }
         if(style==3){
            p0[ix + iz*nx].u= p1[ix + iz*nx].u= -amp*sin(theta);
            p0[ix + iz*nx].w= p1[ix + iz*nx].w= -amp*cos(theta);
         }
         p0[s+zs*nx].u=0.;
         p0[s+zs*nx].w=0.;
         p1[s+zs*nx].u=0.;
         p1[s+zs*nx].w=0.;
      }
   }
   for( iz= zs-ww; iz <= zs +ww; iz++)
   {
      if(iz < 0 ) continue;
      for( ix= s-ww; ix <= s +ww; ix++)
      {
         p0[ix + iz*nx].u /= ampmax;
         p1[ix + iz*nx].u /= ampmax;
         p0[ix + iz*nx].w /= ampmax;
         p1[ix + iz*nx].w /= ampmax;
      }
   }

   /* diagnostic loop, style = 10 */
   if(style==10)for( iz= zs-ww; iz <= zs +ww; iz++)
   {
      for( ix= s-ww; ix <= s +ww; ix++)
      {
         p0[ix + iz*nx].u = 0.;
         p1[ix + iz*nx].u = 0.;
         p0[ix + iz*nx].w = 0.;
         p1[ix + iz*nx].w = 0.;
      }
      p1[s+zs*nx].u=1.0;
      p1[s+zs*nx].w=1.0;
   }
}

#define D    2   /* density */
#define SV   4   /* s-wave velocity */
#define PV   8   /* p-wave velocity */
#define NOTFOUND    0
#define NUMBER      1
#define MATRIX      2
int getabc(a,b,c,nx,nz,ptable,ntable)

/*
	   Subroutine GETABC evalutes the media parameters at
	each grid point and tabulates them, along with    the
	necessary combinations for the finite      difference
	stencils and boundary conditions (except for the free
	surface.)  Comparisons eliminate redundant table  en-
	tries to save space.
*/

float *a, *b, *c;
int nx,nz;
int *ptable;
float **ntable;
{
   float *na;
   int key, cnt, side;
   float dval, svval, pvval, tdtest, tetest;
   int   dkey, svkey, pvkey;
   int   df,   svf,   pvf;
   int len, i, j, k, l, test, y[2], iz, ix, nx2, j2;
   int pt[2];
   int syze;
   float *pa, *pb, *pc;
   float vsmax, vpmax, x, fac, facu, facw;
   float aa, bb, lu2m2, lu2m1, lu2, lu2p1, lu2p2;
   float lu1m2, lu1m1, lu1, lu1p1, lu1p2;
   float l00m2, l00m1, l00, l00p1, l00p2;
   float ld1m2, ld1m1, ld1, ld1p1, ld1p2;
   float ld2m2, ld2m1, ld2, ld2p1, ld2p2;
   float mu2m2, mu2m1, mu2, mu2p1, mu2p2;
   float mu1m2, mu1m1, mu1, mu1p1, mu1p2;
   float m00m2, m00m1, m00, m00p1, m00p2;
   float md1m2, md1m1, md1, md1p1, md1p2;
   float md2m2, md2m1, md2, md2p1, md2p2;
   double sqrt(), fabs();

   /* see what input medium parameters are specified
      currently, 3 and only 3 media parameters   are
      available and must be specified.               */
   dkey= svkey= pvkey= key=0;
   if((svkey=getparm("svelocity", &svf,&svval)))  key |= SV;
   if((pvkey=getparm("pvelocity", &pvf,&pvval)))  key |= PV;
   if((dkey=getparm("density",  &df,&dval)))  key |= D;

   /* construct a=beta**2, b=alpha**2, c=density */

   len= nx*nz;
   nx2 = 2 * nx;
   switch(key)
   {
   case 0: 
   case D: 
   case SV: 
   case PV:
      fprintf(stderr,"must specify three medium parameter files\n");
      exit(-1);
      break;
   case SV|PV: 
   case SV|D: 
   case PV|D: 
      fprintf(stderr,"must specify three medium parameter files\n");
      exit(-1);
      break;
   case PV|SV|D:/* medium specified by 2 velocities and density */
      if(svkey==MATRIX) dread(svf,&a[nx],nx*nz);
      if(pvkey==MATRIX) dread(pvf,&b[nx],nx*nz);
      if(dkey==MATRIX) dread(df,&c[nx],nx*nz);
      if(svkey==MATRIX) close(svf);
      if(pvkey==MATRIX) close(pvf);
      if(dkey==MATRIX) close(df);
      for(pa=a+nx,pb=b+nx,pc=c+nx,cnt=len; cnt--; pa++, pb++, pc++)
      {
         if(svkey==MATRIX) svval= (*pa);
         if(pvkey==MATRIX) pvval= (*pb);
         if(dkey==MATRIX) dval= *pc;
         *pa= svval * svval;
         *pb= pvval * pvval;
         *pc= dval;
      }
      break;
   }
   for(i=0;i<nx;i++){
      a[i]=a[i+nx];
      b[i]=b[i+nx];
      c[i]=c[i+nx];
   }

   /* check for stability  */
   vsmax=0.;
   vpmax=0.;
   for(i=0;i<nx*(nz+1);i++){
      if(sqrt(b[i])>vpmax)vpmax=sqrt(b[i]);
      if(sqrt(a[i])>vsmax)vsmax=sqrt(a[i]);
   }
   fac=sqrt(4./3.)*dt*sqrt(vpmax*vpmax+vsmax*vsmax)/h;
   if(fac>=1.0){
      fprintf(stderr,"fac= %4.2f< instability !\n",fac);
      exit();
   }
   else{
      fprintf(stderr,"fac= %4.2f, stable\n",fac);
   }

   /* convert to b = lambda, a = mu, and multiply by grid factor */
   fac=dt*dt/(h*h);
   for(i=0;i<nx*(nz+1);i++){
      b[i]=(b[i]-2*a[i])*c[i]*fac;
      a[i]*=c[i]*fac;
   }

   /* check for homogeneous regions */
   for(iz=2; iz <  (nz-1);    iz++)
      for(ix=2; ix < (nx-2); ix++)
      {
         len= iz*nx + ix;
         pb= b + len;
         x=pb[0];
         if(fabs((x-pb[-1]  )/x) > HETCHK) continue;
         if(fabs((x-pb[1]   )/x) > HETCHK) continue;
         if(fabs((x-pb[-nx] )/x) > HETCHK) continue;
         if(fabs((x-pb[nx]  )/x) > HETCHK) continue;
         if(fabs((x-pb[-nx2])/x) > HETCHK) continue;
         if(fabs((x-pb[-2]  )/x) > HETCHK) continue;
         if(fabs((x-pb[2]   )/x) > HETCHK) continue;
         if(fabs((x-pb[nx2] )/x) > HETCHK) continue;

	 /* at this point, homogeneous in mu */
         pa= a + len;
         x=pa[0];
         if(fabs(x-pa[-1]  ) > HETCHK*x) continue;
         if(fabs(x-pa[1]   ) > HETCHK*x) continue;
         if(fabs(x-pa[-nx] ) > HETCHK*x) continue;
         if(fabs(x-pa[nx]  ) > HETCHK*x) continue;
         if(fabs(x-pa[-nx2]) > HETCHK*x) continue;
         if(fabs(x-pa[-2]  ) > HETCHK*x) continue;
         if(fabs(x-pa[2]   ) > HETCHK*x) continue;
         if(fabs(x-pa[nx2] ) > HETCHK*x) continue;

	 /* at this point, homogeneous in lambda, also */
         homono+=1;
         c[len]*= (-1);        /* negative density is the homogeneous flag */
      }

   /* check edges for second order homogeneous case */
   for( ix = 1 ; ix < nx-1 ; ix++)
   {
      y[0] = nx+ix;              /* top edge */
      y[1] = nx*(nz-1)+ix;       /* bottom edge */
      for(iz=0;iz<2;iz++)
      {
         if(c[y[iz]]<0)continue;
         x=b[y[iz]];
         if(fabs((x-b[y[iz] -1])/x)>HETCHK)continue;
         if(fabs((x-b[y[iz] +1])/x)>HETCHK)continue;
         if(fabs((x-b[y[iz]-nx])/x)>HETCHK)continue;
         if(fabs((x-b[y[iz]+nx])/x)>HETCHK)continue;

	 /* at this point, homogeneous in mu */
         x=a[y[iz]];
         if(fabs(x-a[y[iz] -1])>HETCHK*x)continue;
         if(fabs(x-a[y[iz] +1])>HETCHK*x)continue;
         if(fabs(x-a[y[iz]-nx])>HETCHK*x)continue;
         if(fabs(x-a[y[iz]+nx])>HETCHK*x)continue;

	 /* at this point, homogeneous in lambda, also */
         homono+=1;
         c[y[iz]]*= (-1);
      }
   }
   for(iz=2;iz<nz-1;iz++)
   {
      y[0] = nx*iz+1;            /* left edge */
      y[1] = nx*(iz+1)-2;        /* right edge */
      for(ix=0;ix<2;ix++)
      {
         if(c[y[ix]]<0)continue;
         x=b[y[ix]];
         if(fabs((x-b[y[ix] -1])/x)>HETCHK)continue;
         if(fabs((x-b[y[ix] +1])/x)>HETCHK)continue;
         if(fabs((x-b[y[ix]-nx])/x)>HETCHK)continue;
         if(fabs((x-b[y[ix]+nx])/x)>HETCHK)continue;
         x=a[y[ix]];
         if(fabs(x-a[y[ix] -1])>HETCHK*x)continue;
         if(fabs(x-a[y[ix] +1])>HETCHK*x)continue;
         if(fabs(x-a[y[ix]-nx])>HETCHK*x)continue;
         if(fabs(x-a[y[ix]+nx])>HETCHK*x)continue;
         homono+=1;
         c[y[ix]]*= (-1);
      }
   }
   fprintf(stderr,"homo= %d, hetero= %d \n",homono, nx*nz-homono);

   /* encode media parameters in table */
   fprintf(stderr,"encode media parameters\n");
   *ntable = (float *) malloc(53*NTSIZE*sizeof(float));
   na = *ntable;
   fprintf(stderr,"starting ntable size: %d entries. (%d bytes, %.2f Mbytes)\n"
	  ,NTSIZE,NTSIZE*53*sizeof(float),NTSIZE*53*sizeof(float)/1000000.);
   fflush(stderr);

   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;
         pc = c + i;
         test=0;
	 if (pc[0] < 0) {
	    aa = -pa[0] / (pc[0] * 12.);
	    bb = -(pb[0] + 2. * pa[0]) / (pc[0] * 12.);
	 }
	 else {
	    lu2m2 =  (pb[-2] + pa[-nx2]) / (pc[0] * 144.);
	    lu2m1 = -(pb[-1] + pa[-nx2]) / (pc[0] * 18.);
	    lu2   = -(pa[ 0] + pa[-nx2]) / (pc[0] * 24.);
	    lu2p1 =  (pb[ 1] + pa[-nx2]) / (pc[0] * 18.);
	    lu2p2 = -(pb[ 2] + pa[-nx2]) / (pc[0] * 144.);
	    lu1m2 = -(pb[-2] + pa[-nx ]) / (pc[0] * 18.);
	    lu1m1 =  (pb[-1] + pa[-nx ]) / (pc[0] * 2.25);
	    lu1   =  (pa[ 0] + pa[-nx ]) / (pc[0] * 1.5);
	    lu1p1 = -(pb[ 1] + pa[-nx ]) / (pc[0] * 2.25);
	    lu1p2 =  (pb[ 2] + pa[-nx ]) / (pc[0] * 18.);
	    l00m2 = -(pb[0] + 2*pa[0] + pb[-2] + 2*pa[-2]) / (pc[0] * 24.);
	    l00m1 =  (pb[0] + 2*pa[0] + pb[-1] + 2*pa[-1]) / (pc[0] * 1.5);
	/*  l00   = see below                                     */
	    l00p1 =  (pb[0] + 2*pa[0] + pb[ 1] + 2*pa[ 1]) / (pc[0] * 1.5);
	    l00p2 = -(pb[0] + 2*pa[0] + pb[ 2] + 2*pa[ 2]) / (pc[0] * 24.);
	    ld1m2 =  (pb[-2] + pa[ nx ]) / (pc[0] * 18.);
	    ld1m1 = -(pb[-1] + pa[ nx ]) / (pc[0] * 2.25);
	    ld1   =  (pa[ 0] + pa[ nx ]) / (pc[0] * 1.5);
	    ld1p1 =  (pb[ 1] + pa[ nx ]) / (pc[0] * 2.25);
	    ld1p2 = -(pb[ 2] + pa[ nx ]) / (pc[0] * 18.);
	    ld2m2 = -(pb[-2] + pa[ nx2]) / (pc[0] * 144.);
	    ld2m1 =  (pb[-1] + pa[ nx2]) / (pc[0] * 18.);
	    ld2   = -(pa[ 0] + pa[ nx2]) / (pc[0] * 24.);
	    ld2p1 = -(pb[ 1] + pa[ nx2]) / (pc[0] * 18.);
	    ld2p2 =  (pb[ 2] + pa[ nx2]) / (pc[0] * 144.);
	    l00   = -(lu2 + lu1 + l00m2 + l00m1 + l00p1 + l00p2 + ld1 + ld2);
	    mu2m2 =  (pa[-2] + pb[-nx2]) / (pc[0] * 144.);
	    mu2m1 = -(pa[-1] + pb[-nx2]) / (pc[0] * 18.);
	    mu2   = -(pb[0] + 2*pa[0] + pb[-nx2] + 2*pa[-nx2]) / (pc[0] * 24.);
	    mu2p1 =  (pa[ 1] + pb[-nx2]) / (pc[0] * 18.);
	    mu2p2 = -(pa[ 2] + pb[-nx2]) / (pc[0] * 144.);
	    mu1m2 = -(pa[-2] + pb[-nx ]) / (pc[0] * 18.);
	    mu1m1 =  (pa[-1] + pb[-nx ]) / (pc[0] * 2.25);
	    mu1   =  (pb[0] + 2*pa[0] + pb[-nx ] + 2*pa[-nx ]) / (pc[0] * 1.5);
	    mu1p1 = -(pa[ 1] + pb[-nx ]) / (pc[0] * 2.25);
	    mu1p2 =  (pa[ 2] + pb[-nx ]) / (pc[0] * 18.);
	    m00m2 = -(pa[ 0] + pa[-2]) / (pc[0] * 24.);
	    m00m1 =  (pa[ 0] + pa[-1]) / (pc[0] * 1.5);
	/*  m00   = see below                                     */
	    m00p1 =  (pa[ 0] + pa[ 1]) / (pc[0] * 1.5);
	    m00p2 = -(pa[ 0] + pa[ 2]) / (pc[0] * 24.);
	    md1m2 =  (pa[-2] + pb[ nx ]) / (pc[0] * 18.);
	    md1m1 = -(pa[-1] + pb[ nx ]) / (pc[0] * 2.25);
	    md1   =  (pb[0] + 2*pa[0] + pb[ nx ] + 2*pa[ nx ]) / (pc[0] * 1.5);
	    md1p1 =  (pa[ 1] + pb[ nx ]) / (pc[0] * 2.25);
	    md1p2 = -(pa[ 2] + pb[ nx ]) / (pc[0] * 18.);
	    md2m2 = -(pa[-2] + pb[ nx2]) / (pc[0] * 144.);
	    md2m1 =  (pa[-1] + pb[ nx2]) / (pc[0] * 18.);
	    md2   = -(pb[0] + 2*pa[0] + pb[ nx2] + 2*pa[ nx2]) / (pc[0] * 24.);
	    md2p1 = -(pa[ 1] + pb[ nx2]) / (pc[0] * 18.);
	    md2p2 =  (pa[ 2] + pb[ nx2]) / (pc[0] * 144.);
	    m00   = -(mu2 + mu1 + m00m2 + m00m1 + m00p1 + m00p2 + md1 + md2);
	 }
	 if (pc[0] < 0) {
            for(j=len-53;j>=0 && test==0;j-=53){
               pt[0]=j;
               if(pa[0]==ta(0) && pb[0]==tb(0) && pc[0]==tc(0)) 
               {
                  test=1;
                  ptable[i]=pt[0];
               }
	    }
         }
	 else {

	    /* first check the four nearest neighbors already tabulated */
	    for(j2=0;j2<4 && test==0;j2++) {

	       if (j2==0 && k==2) continue;
	       if (j2==1 && (k==2 || l==2)) continue;
	       if (j2==2 && l==2) continue;
	       if (j2==3 && (k==nx-2 || l==2)) continue;

	       if (j2==0) pt[0] = ptable[i-1];
	       else       pt[0] = ptable[i-nx+j2-2];
	       if (l00==tl00 && m00==tm00) {
	          if (lu2m2==tlu2m2 && ld2p2==tld2p2 && 
		      mu2m2==tmu2m2 && md2p2==tmd2p2) {
		     if (lu2m1==tlu2m1 && lu2==tlu2 && lu1m2==tlu1m2
			 && lu1m1==tlu1m1 && lu1==tlu1 && l00m2==tl00m2
		         && l00m1==tl00m1 && l00p1==tl00p1 && l00p2==tl00p2) {
		     if (ld1==tld1 && ld1p1==tld1p1 && ld1p2==tld1p2
			 && ld2==tld2 && ld2p1==tld2p1 && mu2m1==tmu2m1
		         && mu2==tmu2 && mu1m2==tmu1m2 && mu1m1==tmu1m1
			 && mu1==tmu1) {
		     if (m00m2==tm00m2 && m00m1==tm00m1 && m00p1==tm00p1
			 && m00p2==tm00p2 && md1==tmd1 && md1p1==tmd1p1
			 && md1p2==tmd1p2 && md2==tmd2 && md2p1==tmd2p1) {
                        test=1;
                        ptable[i]=pt[0];
                     }
                     }
                     }
		  }
	       }
	    }
	    if (test==0) {

	       /* check all remaining table entries */
               for(j=len-106;j>=0 && test==0;j-=53){
                  pt[0]=j;
	          if (l00==tl00 && m00==tm00) {
	             if (lu2m2==tlu2m2 && ld2p2==tld2p2 && 
		         mu2m2==tmu2m2 && md2p2==tmd2p2) {
		        if (lu2m1==tlu2m1 && lu2==tlu2 && lu1m2==tlu1m2
			    && lu1m1==tlu1m1 && lu1==tlu1 && l00m2==tl00m2
		            && l00m1==tl00m1 && l00p1==tl00p1 && l00p2==tl00p2) {
		        if (ld1==tld1 && ld1p1==tld1p1 && ld1p2==tld1p2
			    && ld2==tld2 && ld2p1==tld2p1 && mu2m1==tmu2m1
		            && mu2==tmu2 && mu1m2==tmu1m2 && mu1m1==tmu1m1
			    && mu1==tmu1) {
		        if (m00m2==tm00m2 && m00m1==tm00m1 && m00p1==tm00p1
			    && m00p2==tm00p2 && md1==tmd1 && md1p1==tmd1p1
			    && md1p2==tmd1p2 && md2==tmd2 && md2p1==tmd2p1) {
                           test=1;
                           ptable[i]=pt[0];
                        }
                        }
                        }
		     }
	          }
	       }
	    }
	 }
         if(test==0){                   /* homogeneous */
            pt[0]=len;
            if(pc[0]<0){
               ta(0)=pa[0];
               tb(0)=pb[0];
               tc(0)=pc[0];
               taa= aa;
               tbb= bb;
               tcd=(bb-aa)/2.;
               tcc=2.-30.*(aa+bb);
               ta16=16.*aa;
               tb16=16.*bb;
               tcd10=10.*tcd;
            }
            else{                       /* heterogeneous */
	       ta(0)  = pa[0];
	       tb(0)  = pb[0];
	       tc(0)  = pc[0];
	       tlu2m2 = lu2m2;
	       tlu2m1 = lu2m1;
	       tlu2   = lu2  ;
	       tlu2p1 = lu2p1;
	       tlu2p2 = lu2p2;
	       tlu1m2 = lu1m2;
	       tlu1m1 = lu1m1;
	       tlu1   = lu1  ;
	       tlu1p1 = lu1p1;
	       tlu1p2 = lu1p2;
	       tl00m2 = l00m2;
	       tl00m1 = l00m1;
	       tl00   = l00  ;
	       tl00p1 = l00p1;
	       tl00p2 = l00p2;
	       tld1m2 = ld1m2;
	       tld1m1 = ld1m1;
	       tld1   = ld1  ;
	       tld1p1 = ld1p1;
	       tld1p2 = ld1p2;
	       tld2m2 = ld2m2;
	       tld2m1 = ld2m1;
	       tld2   = ld2  ;
	       tld2p1 = ld2p1;
	       tld2p2 = ld2p2;
	       tmu2m2 = mu2m2;
	       tmu2m1 = mu2m1;
	       tmu2m1 = mu2m1;
	       tmu2   = mu2  ;
	       tmu2p1 = mu2p1;
	       tmu2p2 = mu2p2;
	       tmu1m2 = mu1m2;
	       tmu1m1 = mu1m1;
	       tmu1   = mu1  ;
	       tmu1p1 = mu1p1;
	       tmu1p2 = mu1p2;
	       tm00m2 = m00m2;
	       tm00m1 = m00m1;
	       tm00   = m00  ;
	       tm00p1 = m00p1;
	       tm00p2 = m00p2;
	       tmd1m2 = md1m2;
	       tmd1m1 = md1m1;
	       tmd1   = md1  ;
	       tmd1p1 = md1p1;
	       tmd1p2 = md1p2;
	       tmd2m2 = md2m2;
	       tmd2m1 = md2m1;
	       tmd2   = md2  ;
	       tmd2p1 = md2p1;
	       tmd2p2 = md2p2;
            }
            ptable[i]=pt[0];
            len+=53;
            if(len>=53*NTSIZE){
	       NTSIZE *= 2;
	       syze = NTSIZE * 53 * sizeof(float);
	       fprintf(stderr,
	       "redimensioning ntable to %d entries. (%d bytes, %.2f Mbytes)\n"
	       ,NTSIZE,syze,syze/1000000.);
	       *ntable = (float *) realloc((char *)na,syze);
	       if (*ntable==NULL) {
		  fprintf(stderr,"Memory re-allocation for ntable fails! (try again)\n");
		  NTSIZE *= 0.55;
	          syze = NTSIZE * 53 * sizeof(float);
	          fprintf(stderr,
	          "redimensioning ntable to %d entries. (%d bytes, %.2f Mbytes)\n"
	          ,NTSIZE,syze,syze/1000000.);
	          *ntable = (float *) realloc((char *)na,syze);
	          if (*ntable==NULL) {
		     fprintf(stderr,"Memory re-allocation for ntable fails! (quit)\n");
		     exit(-3);
		  }
	       }
	       if (*ntable != na) copy(na,*ntable,len);
	       na = *ntable;
	       fflush(stderr);
            }
         }
      }

   /* tabulate second order portion */
   fprintf(stderr,"do second order portion\n");
   for(k=0;k<2*(nx+nz)-10;k++)
   {
      if(k<nx-3) i = k + nx + 1;
      else if( k<nx+nz-5) i = nx*(k-nx+5) - 2;
      else if( k<nx2+nz-8) i= k+nx*(nz-2)-nz+7;
      else i = nx*(k-nx2-nz+10) + 1;
      pa = a + i;
      pb = b + i;
      pc = c + i;
      test=0;
      for(j=len-53;j>=0;j-=53){
         pt[0]=j;
         if(pa[0]==ta(0)&&pb[0]==tb(0)&&pc[0]==tc(0)){ 
            test=1;
            ptable[i]=pt[0];
         }
      }
      if(test==0){
         pt[0]=len;
         ta(0)=pa[0];
         tb(0)=pb[0];
         tc(0)=pc[0];
         ptable[i]=pt[0];
         len+=53;
         if(len>=53*NTSIZE){
	    NTSIZE += 2*(nx+nz)-10+nx2+2*nz-2;
	    syze = NTSIZE * 53 * sizeof(float);
	    fprintf(stderr,
	    "redimensioning ntable to %d entries. (%d bytes, %.2f Mbytes)\n"
	    ,NTSIZE,syze,syze/1000000.);
	    *ntable = (float *) realloc((char *)na,syze);
	    if (*ntable==NULL) {
	       fprintf(stderr,"Memory re-allocation for ntable fails!\n");
	       exit(-3);
	    }
	    if (*ntable != na) copy(na,*ntable,len);
	    na = *ntable;
	    fflush(stderr);
         }
      }
   }

   /* tabulate edges and                                      */
   /* set the edges to (1-v*dt/h)/(1+v*dt/h) for the absorbing*/
   /* boundary conditions, a[] for u comp, b[] for w comp     */
   fprintf(stderr,"do edges\n");
   for(k=0;k<nx2+2*nz-2;k++)
   {
      if(k==0) {
         i = k;
         side=1;
      }
      else if(k<nx-1) {
         i = k;
         side=0;
      }
      else if(k<nx+nz){ 
         i = nx*(k-nx+2)-1;
         side=1;
      }
      else if(k<nx2+nz-2){ 
         i = k-nx-nz+nx*nz+1;
         side=0;
      }
      else { 
         i = (k-nx2-nz+3)*nx;
         side=1;
      }
      pa = a + i;
      pb = b + i;
      pc = c + i;
      test=0;
      if(side==0){
         facu=sqrt(a[i]/c[i]);
         facw=sqrt((b[i]+2.*a[i])/c[i]);
      }
      else{
         facu=sqrt((b[i]+2.*a[i])/c[i]);
         facw=sqrt(a[i]/c[i]);
      }
      tdtest= (1.0-facu)/(1.0+facu);
      tetest= (1.0-facw)/(1.0+facw);
      for(j=len-53;j>=0;j-=53){
         pt[0]=j;
         if(pa[0]==ta(0)&&pb[0]==tb(0)&&pc[0]==tc(0)&&tdtest==td(0)&&tetest==te(0)) 
         {
            test=1;
            ptable[i]=pt[0];
         }
      }
      if(test==0){
         pt[0]=len;
         ta(0)=pa[0];
         tb(0)=pb[0];
         tc(0)=pc[0];
         td(0)=tdtest;
         te(0)=tetest;
         ptable[i]=pt[0];
         len+=53;
         if(len>53*NTSIZE){
	    NTSIZE += nx2+2*nz-2;
	    syze = NTSIZE * 53 * sizeof(float);
	    fprintf(stderr,
	    "redimensioning ntable to %d entries. (%d bytes, %.2f Mbytes)\n"
	    ,NTSIZE,syze,syze/1000000.);
	    *ntable = (float *) realloc((char *)na,syze);
	    if (*ntable==NULL) {
	       fprintf(stderr,"Memory re-allocation for ntable fails!\n");
	       exit(-3);
	    }
	    if (*ntable != na) copy(na,*ntable,len);
	    na = *ntable;
	    fflush(stderr);
         }
      }
   }
   fprintf(stderr,"done with table formation, %d types of media! \n",len/53);
   NTSIZE = len/53 + 1;
   syze = NTSIZE * 53 * sizeof(float);
   fprintf(stderr,
   "final ntable size: %d entries. (%d bytes, %.2f Mbytes)\n"
   ,NTSIZE,syze,syze/1000000.);
   *ntable = (float *) realloc((char *)na,syze);
   if (*ntable != na) copy(na,*ntable,len);
   fflush(stderr);
   return(syze);
}

copy(a,b,n)
register float *a, *b;
register int n;
{
   while (n--) *b++ = *a++;
}

getparm(name2,file,val)

/* getparm checks the value assigned to the parameter 'name'
identifying it as a number or string, and opening a file if 
the value is a string */

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

   *file= chkopen(buf,name2,0);
   return(MATRIX);
}

chkopen(filename,filetype,n)

/* attempts to open the file 'filename' for reading */
int n;
char *filename, *filetype;
{
   int f;
   char name3[100];
  if(n==0)
   {
   if( (f=open(filename,O_RDONLY,0664)) <= 1)
   {
      fprintf(stderr,"cannot open %s\n",filename);
      exit(-1);
   }
   fprintf(stderr,"File %s assigned file descriptor %d\n",filename,f);
   }
  if(n==1)
   {
   sprintf(name3,"%s.%s",filename,filetype);
   if( (f=open(name3,O_RDONLY,0664)) <= 1)
   {
      fprintf(stderr,"cannot open %s\n",name3);
      exit(-1);
   }
   fprintf(stderr,"File %s assigned file descriptor %d\n",name3,f);
   }
   return(f);
}

isnumber(str)

/* isnumber checks 'str', returning 1 if 'str' is a floating-point number */
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);
}

record(p,nx,zg,s,seisout)
struct disp *p;
int nx, zg, s, seisout;
{
   register int i, i0, i1;
   register float *held;
   register struct disp *pp, *hold;
   if(pressure!=1){
      /*
      hold= (struct disp *) malloc(noff*8);
      */
      hold = comhold;
      if(zg==0 || zg==1){
         pp= p + ioff0 + s;
         for(i=0; i< noff; i++) {
            i0 = i*idoff;
            i1 = i0 + nx;
            hold[i].u= 0.5*(pp[i0].u + pp[i1].u);
            hold[i].w= 0.5*(pp[i0].w + pp[i1].w);
         }
      }
      else{
         pp= p+ (zg+1)*nx + ioff0 + s;
         for(i=0; i< noff; i++) {
            hold[i].u= pp[i*idoff].u;
            hold[i].w= pp[i*idoff].w;
         }
      }
      write(seisout,hold,8*noff);
      /*
      free(hold);
      */
   }
   else{
      /*
      held= (float *) malloc(noff*4);
      */
      held = comheld;
      pp= p+ (zg+1)*nx + ioff0 + s;
      for(i=0; i< noff; i++)
         held[i] = pp[i*idoff+1].u  - pp[i*idoff-1].u 
             + pp[i*idoff+nx].w - pp[i*idoff-nx].w;
      write(seisout,held,4*noff);
      /*
      free(held);
      */
   }
}
kirrecord(p,kirout)
struct disp *p;
int kirout;
{
	register int i;
	register struct disp *pp;
	register float *hold;
	float h2, h3;
	int knx;

	knx = nx/kdx;
	h2 = h*h;
	h3 = h2*h;
	/*
	hold= (float *) malloc(4*knx);
	*/
	hold = comkhold;
	*hold = 0.0;
	hold[knx-1] = 0.0;
	pp= p+ (kz+1)*nx;
	for(i=1; i< (knx-1); i++)		/* (div div u)x */
	{
		hold[i] = ((pp[1+nx+i*kdx].w-pp[nx-1+i*kdx].w
			    +pp[-1-nx+i*kdx].w-pp[1-nx+i*kdx].w)/4.
			    +pp[1+i*kdx].u-2.*pp[i*kdx].u+pp[-1+i*kdx].u)/h2;
	}
	write(kirout,hold,4*knx);
	for(i=1; i< (knx-1); i++)		/* z-derivative */
	{
		hold[i] = (pp[1+i*kdx+nx].u+pp[-1+i*kdx+nx].u
			    -2.*pp[nx+i*kdx].u+2.*pp[-nx+i*kdx].u
		            -pp[1+i*kdx-nx].u-pp[-1+i*kdx-nx].u
			    +pp[1+i*kdx+nx].w-pp[-1+i*kdx+nx].w
			    -2.*pp[1+i*kdx].w+2.*pp[-1+i*kdx].w
		            +pp[1+i*kdx-nx].w-pp[-1+i*kdx-nx].w)/(2*h3);
	}
	write(kirout,hold,4*knx);
	for(i=1; i< (knx-1); i++)		/* (div div u)z */
	{
		hold[i] = ((pp[1+nx+i*kdx].u-pp[nx-1+i*kdx].u
			    +pp[-1-nx+i*kdx].u-pp[1-nx+i*kdx].u)/4.
			    +pp[nx+i*kdx].w-2.*pp[i*kdx].w+pp[-nx+i*kdx].w)/h2;
	}
	write(kirout,hold,4*knx);
	for(i=1; i< (knx-1); i++)		/* z-derivative */
	{
		hold[i] = (pp[1+i*kdx+nx].u-pp[-1+i*kdx+nx].u
			    -2.*pp[1+i*kdx].u+2.*pp[-1+i*kdx].u
		            +pp[1+i*kdx-nx].u-pp[-1+i*kdx-nx].u
			    +pp[i*kdx+2*nx].w-pp[i*kdx-2*nx].w
			    -2.*pp[nx+i*kdx].w+2.*pp[-nx+i*kdx].w)/(2*h3);
	}
	write(kirout,hold,4*knx);
	for(i=1; i< (knx-1); i++)		/* (curl curl u)x */
	{
		hold[i] = -((pp[1+nx+i*kdx].w-pp[nx-1+i*kdx].w
			    +pp[-1-nx+i*kdx].w-pp[1-nx+i*kdx].w)/4.
			    -pp[nx+i*kdx].u+2.*pp[i*kdx].u-pp[-nx+i*kdx].u)/h2;
	}
	write(kirout,hold,4*knx);
	for(i=1; i< (knx-1); i++)		/* z-derivative */
	{
		hold[i] = -(pp[1+i*kdx+nx].w-pp[-1+i*kdx+nx].w
			    -2.*pp[1+i*kdx].w+2.*pp[-1+i*kdx].w
		            +pp[1+i*kdx-nx].w-pp[-1+i*kdx-nx].w
			    -pp[i*kdx+2*nx].u+pp[i*kdx-2*nx].u
			    +2.*pp[nx+i*kdx].u-2.*pp[-nx+i*kdx].u)/(2*h3);
	}
	write(kirout,hold,4*knx);
	for(i=1; i< (knx-1); i++)		/* (curl curl u)z */
	{
		hold[i] = -((pp[1+nx+i*kdx].u-pp[nx-1+i*kdx].u
			    +pp[-1-nx+i*kdx].u-pp[1-nx+i*kdx].u)/4.
			    -pp[1+i*kdx].w+2.*pp[i*kdx].w-pp[-1+i*kdx].w)/h2;
	}
	write(kirout,hold,4*knx);
	for(i=1; i< (knx-1); i++)		/* z-derivative */
	{
		hold[i] = -(-pp[1+i*kdx+nx].w-pp[-1+i*kdx+nx].w
			    +2.*pp[nx+i*kdx].w-2.*pp[-nx+i*kdx].w
		            +pp[1+i*kdx-nx].w+pp[-1+i*kdx-nx].w
			    +pp[1+i*kdx+nx].u-pp[-1+i*kdx+nx].u
			    -2.*pp[1+i*kdx].u+2.*pp[-1+i*kdx].u
		            +pp[1+i*kdx-nx].u-pp[-1+i*kdx-nx].u)/(2*h3);
	}
	write(kirout,hold,4*knx);
}
absorb(p1,p2,pt,na,offset,tang,norm,cnt)
register struct disp *p1, *p2;
register float *na;
register int *pt;
register int tang, norm, cnt;
int offset;
{
   p1 += offset;
   p2 += offset;
   pt += offset;
   while(cnt--)
   {
      p1[0].u= p2[norm].u -td(0)*(p1[norm].u -p2[0].u);
      p1[0].w= p2[norm].w -te(0)*(p1[norm].w -p2[0].w);
      p1 += tang;
      p2 += tang;
      pt += tang;
   }
}
radrec(p,nx,zs,s0,seisout)
struct disp *p;
int nx, zs, s0, seisout;
{
   float *out;
   int iang, ix, iz, i;
   float s, c, ang, fac;
   double sin(), cos();
   fac= 3.14159/180.0;
   out= (float *)malloc(8*nang);

   for(iang=0; iang<nang; iang++)
   {
      ang= ang0 + iang*dang;
      s= sin(fac*ang);
      c= cos(fac*ang);
      ix= (int)(s0 + c*radius + 0.5);
      iz= (int)(zs + s*radius + 0.5);
      if(ix<1 || ix >=nx-1 || iz<1 || iz >=nz)
      {
         out[iang]= 0.0;
         out[nang+iang]= 0.0;
         continue;
      }
      i= iz*nx + ix;
      if( ostyle == 1){
         out[iang] = p[1+i].u - p[-1+i].u + p[nx+i].w - p[-nx+i].w;
         out[nang+iang] = -(p[nx+i].u - p[-nx+i].u - p[1+i].w + p[-1+i].w);
      }
      else{
	 out[iang] = p[i].u;
	 out[iang + nang] = p[i].w;
      }
   }
   write(seisout,out,8*nang);
   free(out);
}

/*Restart Subs */

int op[100];

void save_run(is,it,len,p1,p2)
 int is, it, len;
 WAVE *p1, *p2;
 {
 static int not_called=1, fdr=0;

 if(not_called)
   {
   char filename[128];

   not_called=0;
   op[0]=homono;   /* number of homogeneous grid points */
   op[1]=ns;       /* number of shots */
   op[2]=ds;       /* mesh spacing between shots */
   op[3]=nx;       /* x-dimension of mesh */
   op[4]=nz;       /* z-dimension of mesh */
   op[5]=nt;       /* number of time steps */
   op[6]= *((int *) &h);           /* spatial mesh interval */
   op[7]= *((int *) &dt);          /* time digitization interval */
   op[8]=idt;      /* record at every idt time steps */
   op[9]=s0;       /* initial shot position */
   op[10]=zs;       /* shot depth */
   op[11]=zg;       /* receiver depth */
   op[12]=style;    /* 1= compress. explo., 0 and 2=double couple */
                   /*             3= shear explosion */
   op[13]=ww;       /* width of space-specified source */
   op[14]=tslice;   /* number of timesteps between timeslice output */
   op[15]=ostyle;   /* ostyle=1 for divergence and curl */
                   /* ostyle=0 for u (x-disp) and w (z-disp) */
   op[16]=freesurf; /* 1=freesurf, 0=absorb at top */
   op[17]=rigid;    /* 1 makes all boundaries rigid */
   op[18]=ioff0;    /* distance from shot to first receiver */
   op[19]=idoff;    /* mesh increment between receivers */
   op[20]=noff;     /* number of receivers */
   op[21]=nang;     /* # of angles for rad output */
   op[22]=dang;     /* angular distance bewteen angles for rad output */
   op[23]=radius;   /* radius to angular output */
   op[24]=ang0;     /* starting angle for angular output */
   op[25]=grid_off; /* value to offset grid by in zap_fill */
   op[26]=pressure; /* pressure instead of u & w disp. is output */
   op[27]=kdx;      /* kirchoff output at every kdx grid points */
   op[28]=kz;       /* kirchoff depth */

/* following for analytical source */
   op[29]=rdsource; /* 0 for quick, 1 for mksource, 2 for source */
		   /*	generated by ntde code  */
   op[30]= *((int *) &pvel);   /* alpha for quick init plane wave */
   op[31]= *((int *) &svel);   /* beta for quick init plane wave */
   op[32]=ntend;    /* how many time points source stays on */
   op[33]=source;   /* inside diameter of source box */

/* filenames */
   op[34]=seisfile;
   op[35]=tsfile;
   op[36]=kirfile;
   op[37]=REPORT;

   sprintf(filename,"%s.%s",name,"save");

   if((fdr=open(filename,O_RDWR | O_CREAT | O_TRUNC,0644)) < 2)
     fprintf(stderr,"Cannot creat %s...\n running without save file\n",
	     filename);

   fprintf(stderr,"File %s assigned file descriptor %d\n",filename,fdr);

   } /*end not_called test*/

   if(fdr < 2) return;
   if((int)lseek(fdr,(long)0,0)<0) return;
   write(fdr,p1,NW*len);
   write(fdr,p2,NW*len);
   write(fdr,op,400);
   if(rdsource)
     {
     write(fdr,ring1,   4*NW*8*(source  ));
     write(fdr,ring2,   4*NW*8*(source+1));
     write(fdr,ring3,   4*NW*8*(source+2));
     write(fdr,ring4,   4*NW*8*(source+3));
     write(fdr,sinner1, 4*NW*8*(source  ));
     write(fdr,sinner2, 4*NW*8*(source+1));
     write(fdr,outer1,  4*NW*8*(source+2));
     write(fdr,outer2,  4*NW*8*(source+3));
     write(fdr,p3outer1,4*NW*8*(source+2));
     write(fdr,p3outer2,4*NW*8*(source+3));
     write(fdr,p3inner, 4*NW*(2*source+3)*(2*source+3));
     }
   write(fdr,&is,4);
   write(fdr,&it,4);
   write(fdr,&table_size,4);
   write(fdr,ptable,len);
   write(fdr,ntable,table_size);

   }

   /*get save file */

void get_savefile(s2,is2,it2,len2)
 int *s2, *is2, *it2, *len2;
 {
 int fdr, len, offset, is, s, it, knx;
 mstpar("name","s",name);
 if((fdr=chkopen(name,"save",1)) == 0) exit(-1);
 mstpar("nx","d",&nx);
 mstpar("nz","d",&nz);
 len=nx*(nz+1)*sizeof(float);

 /*allocate displacement arrays*/
 p1=(WAVE *) malloc(NW*len);
 p2=(WAVE *) malloc(NW*len);
 if(p1 == NULL || p2== NULL)
   {
   fprintf(stderr,"Cannot allocate p1, p2 memory %d bytes\n",2*NW*len);
   exit(-1);
   }

 if(read(fdr,p1,NW*len) != NW*len) error("timeslice 1");
 if(read(fdr,p2,NW*len) != NW*len) error("timeslice 2");
 if(read(fdr,op,400) != 400)       error("old parameters");

 getpar("zs","d",&zs);

chkpar(0,"homono","d",&homono);
chkpar(1,"ns","d",&ns);
chkpar(2,"ds","d",&ds);
chkpar(3,"nx","d",&nx);
chkpar(4,"nz","d",&nz);
chkpar(5,"nt","d",&nt);
chkpar(6,"h","f",&h);
chkpar(7,"dt","f",&dt);
chkpar(8,"idt","d",&idt);
chkpar(9,"s0","d",&s0);
chkpar(10,"zs","d",&zs);
chkpar(11,"zg","d",&zg);
chkpar(12,"style","d",&style);
chkpar(13,"ww","d",&ww);
chkpar(14,"tslice","d",&tslice);
chkpar(15,"ostyle","d",&ostyle);
chkpar(16,"freesurf","d",&freesurf);
chkpar(17,"rigid","d",&rigid);
chkpar(18,"ioff0","d",&ioff0);
chkpar(19,"idoff","d",&idoff);
chkpar(20,"noff","d",&noff);
chkpar(21,"nang","d",&nang);
chkpar(22,"dang","f",&dang);
chkpar(23,"radius","f",&radius);
chkpar(24,"ang0","f",&ang0);
chkpar(25,"grid_off","f",&grid_off);
chkpar(26,"pressure","d",&pressure);
chkpar(27,"kdx","d",&kdx);
chkpar(28,"kz","d",&kz);
chkpar(29,"rdsource","d",&rdsource);
chkpar(30,"pvel","f",&pvel);
chkpar(31,"svel","f",&svel);
chkpar(32,"ntend","d",&ntend);
chkpar(33,"source","d",&source);
chkpar(34,"seisfile","d",&seisfile);
chkpar(35,"tsfile","d",&tsfile);
chkpar(36,"kirfile","d",&kirfile);
chkpar(37,"REPORT","d",&REPORT);
/*endpar();*/
zs++;

if(op[29])  /* rdsource = 1*/
  {
  /*
  ring1=(float *)    malloc(4*NW*8*(source  ));
  ring2=(float *)    malloc(4*NW*8*(source+1));
  ring3=(float *)    malloc(4*NW*8*(source+2));
  ring4=(float *)    malloc(4*NW*8*(source+3));
  sinner1  =(WAVE *) malloc(4*NW*8*(source  ));
  sinner2  =(WAVE *) malloc(4*NW*8*(source+1));
  outer1   =(WAVE *) malloc(4*NW*8*(source+2));
  outer2   =(WAVE *) malloc(4*NW*8*(source+3));
  p3outer1 =(WAVE *) malloc(4*NW*8*(source+2));
  p3outer2 =(WAVE *) malloc(4*NW*8*(source+3));
  p3inner  =(WAVE *) malloc(4*NW*(2*source+3)*(2*source+3));
  */

  /* Get source state */
  if(read(fdr,ring1,4*NW*8*(source  ))
     != 4*NW*8*(source  )) error("ring1");
  if(read(fdr,ring2,4*NW*8*(source+1))
     != 4*NW*8*(source+1)) error("ring2");
  if(read(fdr,ring3,4*NW*8*(source+2))
     != 4*NW*8*(source+2)) error("ring3");
  if(read(fdr,ring4,4*NW*8*(source+3))
     != 4*NW*8*(source+3)) error("ring4");
  if(read(fdr,sinner1,4*NW*8*(source  ))
     != 4*NW*8*(source  )) error("sinner1");
  if(read(fdr,sinner2,4*NW*8*(source+1))
     != 4*NW*8*(source+1)) error("sinner2");
  if(read(fdr,outer1,4*NW*8*(source+2))
     != 4*NW*8*(source+2)) error("outer1");
  if(read(fdr,outer2,4*NW*8*(source+3))
     != 4*NW*8*(source+3)) error("outer2");
  if(read(fdr,p3outer1,4*NW*8*(source+2))
     != 4*NW*8*(source+2)) error("p3outer1");
  if(read(fdr,p3outer2,4*NW*8*(source+3))
     != 4*NW*8*(source+3)) error("p3outer2");
  if(read(fdr,p3inner,4*NW*(2*source+3)*(2*source+3))
     != 4*NW*(2*source+3)*(2*source+3)) error("p3inner");
  }

  if(ntend <= 0) ntend=nt;

  if(read(fdr,&is,4) != 4) error("source number");
  if(read(fdr,&it,4) != 4) error("timestep");
  fprintf(stderr,"it=%d\n",it);
  getpar("it","d",&it);/* Option to change it variable */
  endpar();
  fprintf(stderr,"it=%d\n",it);
  if(read(fdr,&table_size,4) != 4) error("table_size");
  fprintf(stderr,"table_size=%d\n",table_size);

  /* get table */
  ptable= (int *) malloc(nx*(nz+1)*sizeof(int));
  ntable= (float *) malloc(table_size);
  if(ptable == NULL || ntable == NULL)
   fprintf(stderr,"cannot allocate memory %d bytes\n",
	  nx*(nz+1)*sizeof(int)+table_size);

  if(read(fdr,ptable,nx*(nz+1)*sizeof(int))
     != nx*(nz+1)*sizeof(int)) error("table pointers");

  if(read(fdr,ntable,table_size)
     != table_size) error("table");

  /* reopen output files */
  if(seisfile)
   {
   if(noff)
   offset=NW*noff*4*(it/idt +1);
   if(nang)
   offset=NW*noff*4*(it/idt +1);
   seisout=chkreopen(name,"seis",offset);
   }
  if(kirfile)
   {
   offset=NW*(nx/kdx)*4*(it/idt +1);
   kirout=chkreopen(name,"kir",offset);
   }
  if(tsfile)
   {
   offset=(it/tslice +1)*(len-nx*sizeof(float))*NW;
   tsout=chkreopen(name,"ts",offset);
   }

  /* allocate output arrays */
  if (seisout)
    {
    comhold= (WAVE *) malloc(noff*8);
    comheld= (float *) malloc(noff*4);
    }
  if (kirout)
    {
    knx = nx/kdx;
    comkhold= (float *) malloc(4*knx);
    }


  /*tabulate free surface parameters*/
  if((fsg=(float *) malloc(sizeof(float)*nx))==NULL)
    {
    fprintf(stderr,"cannot allocate %d bytes for free surface\n",4*nx);
    exit(-1);
    }
  fsfill(fsg,fsu,fsw,ptable,ntable,nx);

  s=s0 + is*ds;
  if(rdsource==0)
    {
    quickinit(p1,p2,nx,s0,zs,style,ww,pvel,svel);
    }
    else
    {
    source1in=chkopen(name,"src",1);
    if(source1in ==0) exit(-1);
    offset=4*it*8*NW*(4*source+6);
    if((int)lseek(source1in,(long)0,2) < offset)
      {
      fprintf(stderr,"Existing file less than %d bytes - can't proceed\n",              offset);
      exit(-1);
      }
    if((int)lseek(source1in,(long)offset,0) != offset)
      {
      fprintf(stderr,"failure to seek in source file on open - exit\n");
      exit(-1);
      }
    }

    *s2=s;
    *is2=is;
    *it2=it;
    *len2=len;
    close(fdr);
}


void
chkpar(n,name,type,pt)
int n, *pt;
char *name, *type;
{
   getpar(name,type,pt);
   if (op[n] != *pt) {
      fprintf(stderr,"warning: parameter %s has changed since run was saved\n",
              name);
      fprintf(stderr," - new value will be used if needed\n");
   }
}

void
error(varname)
char *varname;
{
   fprintf(stderr,"read error in get_savefile:\n");
   fprintf(stderr,"     fatal read error:  reading %s\nquit\n",varname);
}

int
chkreopen(name3,type,offset)

/*
           Function CHKREOPEN attempts to reopen filename, returning
        the resulting file descriptor if successful, or printing  an
        error and aborting the program if not.  It requires the file
        to exist and to contain >= offset bytes.
*/

char *name3, *type;
int offset;
{
   int f;
   char filename[128];

      sprintf(filename,"%s.%s",name,type);
   if ((f = open(filename,O_RDWR,0644)) < 0) {
      fprintf(stderr,"cannot reopen %s\n",filename);
      fprintf(stderr,"trying to create...\n");
      return(chkcreat(name3,type));
   }
   fprintf(stderr,"File %s reopened on file descriptor %d\n",filename,f);
   if ((int)lseek(f,(long)0,2) < offset) {
      fprintf(stderr,"Existing file less than %d bytes long - cannot proceed\n"
              ,offset);
      exit(-1);
   }
   if ((int)lseek(f,(long)offset,0) != offset) {
      fprintf(stderr,"failure to seek in file %s on reopen - exit\n",filename);
      exit(-1);
   }
   return(f);
}