/**********************************************************************/
/***                            Cathy.C                             ***/
/***                      Modified Voter Model                      ***/
/***								    ***/
/***	Extensively modified original program array algorithms	    ***/
/***	to be 1 to 2 orders of magnitude faster for larger array    ***/
/***	sizes.  Also modified output graphics to allow for much	    ***/
/***	higher speeds.						    ***/
/***					Garr Updegraff, 1/21/99	    ***/
/***					School of Biol, UC Irvine   ***/
/***					garru@uci.edu		    ***/
/***								    ***/
/***	6/1/99: Added option for reflective boundaries, instead of  ***/
/***		wrap around boundaries.				    ***/
/***								    ***/
/**********************************************************************/

#define YES 1
#define NO  0

#define PC
#define GRAPHICS

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <math.h>

#ifdef PC
    #include <alloc.h>
    #include <conio.h>
    #include <dos.h>
#endif

#define WRTDETAIL
#define USER_INPUT

#define MAX_GRIDS 250	/* Desire max grid size of 250x250 */


/***  global variables  ***/
FILE	*inptr1, *outptr1;

unsigned char cm[2] = {0, 9},	/* color palatte */
	      bAlphaPos;	/* 1 (true) if alpha is pos, 0 (false) if alpha is neg */

float	fAlpha,		/* a part of dispersal equation */
	propA;		/* defines initial proportion of A */

long	lDispRegion,	/* actual dispersal region: sqr(2*dispersal-1) */
	lInteractRegion;/* actual interaction region: sqr(2*interaction-1) */

short	dispersal,	/* defined integer for dispersal region, the n in (2n-1)^2 */
	interaction;	/* defined integer for interaction region, the n in (2n-1)^2 */

unsigned char bDisplayGraph,	/* switch: 1 (true) if graphics should be displayed */
	      bReflectiveBoundary, /* switch: 0 for wrap-around, 1 for reflective */
	      bPaintedCenter,	/* switch: 1 (true) for painted center */
	      bBlockFormStart,  /* switch for starting in block formation */
	      iPaintBits = 4;   /* number of bits wide for each plot dot */

long	lStopTime,      /* stopping time */
	lTime;		/* time step */

int	iReplicateSeed,	/* random seed, incr's with each Replicate */
	iGph_Interval,   /* timestep interval for graphing */
	nFinalReplicates, /* final number of main program loop cycles */
	nReplicates,	/* current number of main program loop cycles */
	nGridX = 100,	/* grid length in x dimension */
	nGridY = 100;	/* grid length in y dimension */

char	sTemp[119];	/* just a temporary string */
	sOutFile[149];	/* name of output data file */

unsigned char bFullProbab;  /* 1 (true) if afProbab[] contains full probability, */
			    /* 0 (false) if it just contains habitat fraction partial calc's. */
#ifdef MF
    unsigned char *a[MAX_GRIDS],  /*grid, initial and grid to graph from*/
		  *b[MAX_GRIDS];  /*temp holding grid for dispersal sums during voter procedure*/
		  *c[MAX_GRIDS];  /*temp holding grid for interaction sums during " */
    float *afProbab;	/* If user specifies dispersal=interaction region size, */
			/*	then this contains pre-calculated probabilities */
			/*	for any dispersal=interaction regions sums. */
			/* Otherwise, it contains habitat fractions based on */
			/*	interaction region sums as partial pre-calculations */
#endif

#ifdef PC
    unsigned char far *a[MAX_GRIDS], far *b[MAX_GRIDS], far *c[MAX_GRIDS];
    float far *afProbab; /* See comments in MF section above. */

    static unsigned char far *egaptr =	(unsigned char far *)0xa0000000L;
    static unsigned long far *INT43 = (unsigned long far *)0x00000010c;
#endif

/***  Procedures  ***/
void   cluster(int *);	/* calculate cluster and % occupied */
void   initial(int nReplicates);/* initialize the grid */
void   input(void);		/* ask user for input */
void   run_init(void);		/* allocate memory for program */
void   r250_init(int);		/* fire up random procedure dr250() */
void   output(int *, int nReplicates);	/* write to output file */
void   voter( void );		/* do voter calculations */
double dr250(void);		/* random number procedure */

#ifdef PC
    void  graphit(int *);	/* graph the grid to screen*/
    void  PaintRow( short iRow, unsigned int wColor );
    void  setvmode(char mode);
    void  wrch(char b, char c, int x, int y);
    void  wrstr(char phrase[], char c, int  x, int y);
#endif


/**********************************************************************/
/***                          Main Program                          ***/
/**********************************************************************/

main()
{
  long	/* dr,	/*dispersal region loop */
	ir,	/*interaction region */
	lGph_Time;	/* time at which graphing occurs to screen */

  int	aiStats[2];	/* aiStats[0]=cluster; aiStats[1]=occupied */

  char   bKb;		/* becomes true when keyboard char is hit */

  double da,		/* proportion of A in dispersal region */
	 db,		/* proportion of B in dispersal region */
	 fa,		/* frequency of A in interaction region */
	 dHabitatQ,	/* Habitat Quality fraction, based on Interaction region sum */
	 dDenom;	/* denominator for fraction */

  input();             /* input run parameters */

  bFullProbab = ( dispersal==interaction );  /* 1 (true) if they match, 0 otherwise */
  if ( nGridY <= 110 ) iPaintBits = 4;
  else if ( nGridY <= 220 ) iPaintBits = 2;
  else iPaintBits = 1;

  run_init();     /***  set up run environment using input()'s specs. ***/

  /* Precalculate probabilities, or at least habitat quality fractions
  /*	for interaction region sums. */
  for (ir=0; ir<=lInteractRegion; ++ir)  {  /*cycles through interaction region */
    printf (".");
    fa = 0.5 - ((double)ir) / lInteractRegion;
    /* calculate habitat quality (between 0 and 1) for all possible interaction sizes. */
    if ( bAlphaPos )  dHabitatQ = 0.5 - fAlpha*fa;
    else dHabitatQ = 0.5 + fAlpha*fa;
    afProbab[(int)ir] = dHabitatQ;  /* store habitat calc, in case not bFullProbab */

    /* if interact region = dispersal region, then pre-calc full probabilities */
    if ( bFullProbab ) {  /* interact region=dispersal region */
      da = ((double)ir) / lDispRegion;	/* fraction present in Disp region */
      db = 1.0 - da;		/* fraction not present in Disp region (species 2) */
      dDenom = da * dHabitatQ + db * (1.00-dHabitatQ);
      if (dDenom == 0) afProbab[(int)ir] = 0.0;
      else afProbab[(int)ir] = dHabitatQ*da / dDenom;
      }
    }    /* end for dr */

  #ifdef GRAPHICS
    if ( bDisplayGraph )  setvmode(0x12);  /* activate 16-color 640x480 VGA graphics mode */
  #endif

  /* Perform main program loops -- default is 100 replicates. */
  /* ... stop if any keyboard activity */
  for ( nReplicates=1;  !kbhit() && (nReplicates <= nFinalReplicates);  nReplicates++, iReplicateSeed++ )  {
    r250_init(iReplicateSeed);

    lGph_Time = 1;     /* graph start time       */
    initial(nReplicates); /***  initialize grid, reset lTime  ***/
    cluster(aiStats);    /*** calculate statistic for clustering  ***/

    if ( bDisplayGraph )  {
      #ifdef GRAPHICS
      graphit(aiStats);	/*** put initial grid on screen  ***/
      #endif
      }
    else printf( "Rep %d ", nReplicates );  /* no graphics, just display Replicate */

    /* printf("lTime=%d, lStopTime=%d, aiStats[1]=%d, Hit any key...\n", lTime, lStopTime, aiStats[1]);
    /* scanf("%20s", sTemp);  /* debug */

    /***  loop thru generations for this replicate  ***/
    #ifdef PC
    while (lTime<lStopTime && aiStats[1]!=0 && aiStats[1]!=100)  {
    #endif
    #ifdef MF
    while (lTime<lStopTime && aiStats[1]!=0 && aiStats[1]!=100)  {
    #endif
      lTime++;
      voter();	/***  calculate values in next time step  ***/
      cluster(aiStats);

      #ifdef GRAPHICS
	if (lTime == lGph_Time)  {  /* update graphics display */
	  lGph_Time += iGph_Interval;
	  if ( bDisplayGraph )  graphit(aiStats);
	  else printf (".");  /* no graphics */
	  /*	else printf( "g%d, ", nReplicates );  /* debug */
	  }
      #endif
      #ifdef PC
	bKb = kbhit();
	if ( bKb ) break;
      #endif
      }     /* end "while" generation loop  */
      output(aiStats, nReplicates);  /* write results for this replicate to output file */

      /* printf( "lTime=%ld, lStopTime=%ld, aiStats[1]=%d, bKb=%d\n", lTime, lStopTime, aiStats[1], (int)bKb );
      /* printf( "Hit 1 and Enter for next loop..." );
      /* scanf( "%f", &fAlpha );  /* */

    }   /* end for nReplicates loop */


  #ifdef GRAPHICS
    if (bDisplayGraph )
      setvmode(0x03);	/* this resets screen to text mode */
  #endif

  #ifdef WRTDETAIL
    fclose (outptr1);
  #endif

  return( 0 );
  } /* end for main program*/


/**********************************************************************/
/***                      Input Run Parameters                      ***/
/**********************************************************************/

void input(void)
{
  int	L;
  short bErr;	/* 1 (true) on input errors */

  iGph_Interval = 1;

  #ifdef USER_INPUT
    printf("*** Voter Model ver 1.00 ***\n\n");

    printf( "Name of output data file: " );
    scanf( "%148s", sOutFile );
    if ( strlen( (char *)sOutFile ) <= 0 )  {
      strcpy( (char *)sOutFile, "Coexist.dat" );
      printf( "   Using Default Output Name '%s'\n", sOutFile );
      }

    do {
      printf ("Display graphics while running (1=yes; 0=no)? ");
      scanf ( "%d", &L );	bDisplayGraph = (unsigned char)L;
      bErr = (bDisplayGraph & 1) != bDisplayGraph;  /* error if not 0 or 1 */
      if ( bErr )  /*** force user to enter 0 or 1 for graphics ***/
	printf ("*** Please enter 0 or 1.\n");
      } while ( bErr );

    do {
      printf( "Grid size (50 to 250) = " );
      scanf( "%d", &nGridX );
      bErr = ( (nGridX < 50) || (nGridX > 250) );
      if ( bErr )  /*** force user to enter proper range ***/
	printf ("*** Grid size not within bounds.\n");
      } while ( bErr );
    nGridY = nGridX;	/* use same bounds for Y axis */

    printf("Number of replicates (main program loops) = ");
    scanf( "%d", &nFinalReplicates );

    printf("Seed for first replicate (normally 101) = ");
    scanf( "%d", &iReplicateSeed );

    printf("Stopping time = ");
    scanf("%ld", &lStopTime);

    printf("Graph interval = ");
    scanf("%d", &iGph_Interval);

    printf("Use the following scheme for dispersal and interaction regions \n");
    printf("   0: 5 cell star, for this case, Disperal=Interact Region \n");
    printf("   1: 1 cell region\n");
    printf("   2: 9 cell region\n");
    printf("   3: 25 cell region\n");
    printf("    . \n");
    printf("    . \n");
    printf("    . \n");
    printf("   %d:", ((100 + 1)/2));
    printf(" entire grid area (max allowable n) \n");

    do {
      printf("Enter dispersal region: ");
      scanf("%hd", &dispersal);
      bErr = ( (dispersal < 0) || (dispersal > 50) );
      if ( bErr )  /*** force user to enter valid dispersal value ***/
	printf( "*** dispersal is not within correct limits.  Try again.\n" );
      } while ( bErr );
    if (dispersal <= 0)  lDispRegion = 5;
    else {
      L = 2 * dispersal - 1;
      lDispRegion = L * L;
      }

    do {
      printf("Enter interaction region: ");
      scanf("%hd", &interaction);
      bErr = ((interaction < 0) || (interaction > 50));
      if ( bErr )  /*** force user to enter value between 0 and 50 ***/
	printf( "*** Interaction is not within correct limits.  Try again.\n" );
      } while ( bErr );
    if ( interaction <= 0 )  lInteractRegion = 5;
    else {
      L = 2 * interaction - 1;
      lInteractRegion = L * L;
      }

    do {
      printf("\nEnter alpha between -1 and 1: " );
      scanf("%f", &fAlpha);
      bAlphaPos = ( fAlpha >= 0 ); /* 0 for neg fAlpha, 1 for pos */
      fAlpha = fabs(fAlpha);	/* store absolute value */
      bErr = fAlpha > 1.0;
      if ( bErr )  /*** force user to enter valid Alpha value ***/
	printf ("*** alpha is not within correct limits; try again.");
      } while ( bErr );

    printf ("\nDefine the simulation boundary behavior:");
    do {
      printf("\n""Wrap-around"" or ""Reflective"" boundaries (0 or 1)? ");
      scanf ( "%d", &L );	bReflectiveBoundary = (unsigned char)L;
      bErr = (bReflectiveBoundary & 1) != bReflectiveBoundary;  /* error if not 0 or 1 */
      if ( bErr )  /*** force user to enter 0 or 1 for boundary conditions ***/
	printf ("*** Please enter 0 or 1.\n");
      } while ( bErr );

    printf ("\nDefine the starting conditions:\n");
    do {
      printf ("Would you like a 10x10 painted center (0=no 1=yes)? ");
      scanf ( "%d", &L );	bPaintedCenter = (unsigned char)L;
      bErr = (bPaintedCenter & 1) != bPaintedCenter;    /* error if not 0 or 1 */
      if ( bErr )   /*** force user to enter 0 or 1 ***/
	printf ( "*** Please enter 0 or 1.\n" );
      } while ( bErr );

    if ( !bPaintedCenter )  {  /* no painted center */
      do  {
	printf( "Enter the proportion (0.0 to 1.0) of species A\n");
	printf( "    in the initial grid: " );
	scanf ( "%f", &propA );
	bErr = (propA <  0.0) || (propA > 1.0);
	if ( bErr )	/*** for user to enter value between 0 and 1 ***/
	  printf( "*** Specified initial proportion of A is out of bounds.  Try again.\n" );
	} while ( bErr );

      if (propA != 0.50 )  bBlockFormStart=1;
      else do {
	printf( "Would you like a random placement or a block design? \n");
	printf( "    Enter 0 for block, or 1 for random placement: ");
	scanf ( "%d", &L );	bBlockFormStart = (unsigned char)L;
	bErr = (bBlockFormStart & 1) != bBlockFormStart;  /* error if not 0 or 1 */
	if ( bErr )  /*** force user to enter 0 or 1 ***/
	  printf ("*** Invalid response.  Try again.\n" );
	} while ( bErr );
      }   /*end for bPaintedCenter*/
  #endif

  }


/**********************************************************************/
/***                     Setup Run Environment                      ***/
/**********************************************************************/

void run_init(void)
{
  short i;
  unsigned int iSizeY;

  /***  open output file  ***/
  #ifdef WRTDETAIL
    /* outptr1=fopen ("100.dat", "wb");  /* old */
    #ifdef MF
      outptr1 = fopen( (char *)sOutFile, "wb" );  /* binary output */
    #endif
    #ifdef PC
      outptr1 = fopen( (char *)sOutFile, "wt" );  /* text output */
    #endif

    if (outptr1 == (FILE *) NULL)  {
      printf ("*** Cannot open %s\n", sOutFile );
      exit(2);
      }
  #endif

  /***  calculate storage requirements for the arrays  ***/
  /* If dispersal=interaction, then work can pre-calculcate full */
  /*	probabilities for any possible region sums in probab[]. */
  /* Otherwise, must do partial pre-calculation where probab[] */
  /*	contains habitat fractions. */

  #ifdef MF
    for (i = 0; i < nGridX; i++)  {
      a[i] = (unsigned char *) malloc( nGridY * sizeof(a[1]) );
      b[i] = (unsigned char *) malloc( nGridY * sizeof(b[1]) );
      }
    afProbab = (float *) malloc( (lInteractRegion+1) * sizeof(afProbab[1]) );
  #endif

  #ifdef PC
    /* Actual array size must be divisible by 8 for graphing routines; */
    /*     add 7 and mask off low bits to get larger array (if necessary) */
    /*     divisible by 8. */
    iSizeY = 7;		/* set 3 low bits */
    iSizeY = ~iSizeY & (nGridY + 7); /* reverse bits for AND mask */

    for (i = 0; i < nGridX; i++)  {
      a[i] = (unsigned char far *) farmalloc( iSizeY * sizeof(a[0][0]));
      memset( a[i], 0, iSizeY * sizeof(a[0][0]));
      b[i] = (unsigned char far *) farmalloc( iSizeY * sizeof(b[0][0]));
      memset( b[i], 0, iSizeY * sizeof(b[0][0]));
      c[i] = (unsigned char far *) farmalloc( iSizeY * sizeof(c[0][0]));
      memset( c[i], 0, iSizeY * sizeof(c[0][0]));
      }
    afProbab = (float far *) farmalloc( (lInteractRegion+1) * sizeof(afProbab[1]) );
  #endif
  }


/**********************************************************************/
/***                        Initialize Grid                         ***/
/**********************************************************************/

void initial(int nReplicates)
{
  static char sGrid[] =  {"Grid Size   =         "};
  static char sDispl[] = {"Dispersion  =         "};
  static char sInter[] = {"Interaction =         "};
  static char sTime[]  = {"Time        =         "};
  static char sStat1[] = {"% Clustered =         "};
  static char sStat2[] = {"% Occupied  =         "};
  static char sSeed[]  = {"Random Seed =         "};
  static char sReplicate[] = {"Replicate   =         "};

  /* char cdispl[5], cinter[5]; */
  float prob;
  short i, j;
  char  S[31];

  lTime = 0;       /*** initialize time  ***/

  if (bBlockFormStart != 0)  {  /*** initialize working grid at random  ***/
    for (i = 0; i < nGridX; i++)
      for (j = 0; j < nGridY; j++)  {
	prob = dr250();
	a[i][j] = (prob < propA);  /* 1 if true, 0 otherwise */
	}
    }

  if (bBlockFormStart != 1)  {
    /*** initialize working grid as a 50/50 block  ***/
    for (i = 0; i < nGridX; i++)
      for (j = 0; j < nGridY; j++)
	a[i][j] = (j > nGridY/2);  /* 1 if true, 0 otherwise */
    }

  if ( bPaintedCenter )  {
    for (i=0; i<nGridX; ++i)
      for (j=0; j<nGridY; ++j)
	a[i][j] = (i>43 && j>43 && i<55 && j<55);  /* 1 if true, 0 otherwise */
    }

  /***  open screen for graphics  ***/
  #ifdef GRAPHICS
  if ( bDisplayGraph )  {
    /* setvmode(0x12); */

    wrstr(sGrid,  12, 125, 57);
    wrstr(sInter, 12, 150, 57);
    wrstr(sDispl, 12, 175, 57);
    wrstr(sTime,  12, 200, 57);
    wrstr(sStat1, 12, 250, 57);
    wrstr(sStat2, 12, 275, 57);
    wrstr(sReplicate, 12, 300, 57);
    wrstr(sSeed, 12, 325, 57);

    /* write values to screen that won't change for this replicate */

    /* write grid size to the screen */
    itoa(nGridY, S, 10);	  wrstr(S, 12, 125, 71);
    /* write interaction region size to the screen */
    itoa(interaction, S, 10);	  wrstr(S, 12, 150, 71);
    /* write dispersal region size to the screen */
    itoa(dispersal, S, 10);	  wrstr(S, 12, 175, 71);
    /* write Replicate number to the screen  ***/
    itoa(nReplicates, S, 10);	  wrstr(S, 12, 300, 71); /* */
    /* write Replicate Seed to the screen  ***/
    itoa(iReplicateSeed, S, 10);  wrstr(S, 12, 325, 71); /* */

    }
  #endif
  }


/**********************************************************************/
/***             Calculate Presence Based On Randomness             ***/
/***             Dispersal and Interaction Region Sums              ***/
/**********************************************************************/
char cPresence( long lDisperSum, long lInterSum )
{
  double da, db, dDenom, dHabitatQ, dProb;  /* variables used in calc'ing probability */
  float  compare;

  compare = dr250();
  if ( bFullProbab )  /* Dispersal region same as Interaction region */
    return (compare <= afProbab[(int)lDisperSum]);   /* 1 if true, 0 otherwise */

  else {	/* Dispersal region not same as interaction region. */
    /* Need to finish probability calculations. */
    dHabitatQ = afProbab[(int)lInterSum];  /* this much already calculated */
    da = ((double)lDisperSum) / lDispRegion;  /* fraction present in Disp region */
    db = 1.0 - da;		/* fraction not present in Disp region (species 2) */
    dDenom = da * dHabitatQ + db * (1.00-dHabitatQ);
    if (dDenom == 0) dProb = 0.0;
    else dProb = dHabitatQ*da / dDenom;
    return (compare <= dProb);   /* 1 if true, 0 otherwise */
    }
  }


/**********************************************************************/
/***                          Voter Step                            ***/
/**********************************************************************/

void voter()
{
  int i, j, k, l, k2, l2;
  long	lDisperSum, lInterSum;  /* Respective sums for Dispersal and Interaction regions */

  #ifdef MF
    unsigned char *pTemp[1];	  /* temporary pointer for swapping arrays */
  #endif
  #ifdef PC
    unsigned char far *pTemp[1];
  #endif

  /* Sum surrounding region for each cell */
  if (dispersal <= 0)  {  /*** voter step for dispersal=0 (star) ***/
    /* Wrap indexes if less or greater than array bounds. */
    /* ...Star pattern involves one cell above, below, left, and right of center. */
    for ( i=0;  i<nGridX;  i++ )  {
      k = i - 1;	/* set left star index k */
      if ( k < 0 )	/* beyond left boundary */
	if ( bReflectiveBoundary ) k = 1;	/* reflect */
	else k += nGridX;			/* wrap-around */
      k2 = i + 1;	   /* set right star index k2 */
      if ( k2 >= nGridX )  /* beyond right boundary */
	if ( bReflectiveBoundary ) k2 -= 2;	/* reflect */
	else k2 = 0;				/* wrap-around */
      for ( j=0;  j<nGridY;  j++ )  {
	lDisperSum=0;
	l = j - 1;	/* set upper star index L */
	if ( l < 0 )	/* beyond upper boundary */
	  if ( bReflectiveBoundary ) l = 1;	/* reflect */
	  else l += nGridY;			/* wrap-around */
	l2 = j + 1;		/* set bottom star index L2 */
	if ( l2 >= nGridY )	/* beyond bottom boundary */
	  if ( bReflectiveBoundary ) l2 -= 2;	/* reflect */
	  else l2 = 0;				/* wrap-around */
	/* sum the five cells within "star" pattern */
	lDisperSum = a[k][j] + a[i][l] + a[i][j] + a[k2][j] + a[i][l2];
	lInterSum=lDisperSum;
	b[i][j] = cPresence( lDisperSum, lInterSum );
	}
      }
    /* fast copy b[] array values to a[] by swapping pointers */
    for (i = 0; i < nGridX; i++)  {
      pTemp[0] = a[i];	a[i] = b[i];	b[i] = pTemp[0];
      }
    }		/* end star case   */

  else {  /* Start summing regions for non-star cases (dispersal > 0) */
    for ( j=0;  j<nGridY;  j++ )  {
      /* sum displacement region along each column and store in b[][] */
      lDisperSum = 0;		lInterSum = 0;
      /* do column's top-most dispersal region first */
      for ( i=-dispersal+1;  i<dispersal; i++ )  {
	k = i;
	if ( k < 0 )  /* beyond left boundary */
	  if ( bReflectiveBoundary ) k = -k;	/* reflect */
	  else k += nGridX;			/* wrap-around */
	lDisperSum += a[k][j];  /* add cell, if present */
	}
      if ( !bFullProbab )	/* need to sum different sized interaction region */
	for ( i=-interaction+1;  i<interaction; i++ )  {
	  k = i;
	  if ( k < 0 )  /* beyond left boundary */
	    if ( bReflectiveBoundary ) k = -k;	/* reflect */
	    else k += nGridX;			/* wrap-around */
	  lInterSum += a[k][j];  /* add cell, if present */
	  }

      b[0][j] = (int)lDisperSum;   /* store value for left most row region */
      if ( !bFullProbab )  /* need to sum different sized interaction region */
	c[0][j] = (int)lInterSum;  /* store value for left most row region */

      /* Work across row, subtracting left-most cell from old Sum */
      /*	and adding right-most cell in new region */
      for ( i=1;  i<nGridX;  i++ )  {
	k = i - dispersal;
	if ( k < 0 )	/* beyond left boundary */
	  if ( bReflectiveBoundary )  k = -k;	/* reflect */
	  else k += nGridX;			/* wrap-around */
	lDisperSum -= a[k][j];
	k = i + dispersal - 1;
	if ( k >= nGridX )  /* beyond right boundary */
	  if ( bReflectiveBoundary )  k = 2*nGridX - k - 2;  /* reflect */
	  else k -= nGridX;			/* wrap-around */
	lDisperSum += a[k][j];
	b[i][j] = (int)lDisperSum;  /* store value for row region */
	if ( !bFullProbab )  {  /* need to sum different sized interaction regions */
	  k = i - interaction;
	  if ( k < 0 )	/* beyond left boundary */
	    if ( bReflectiveBoundary )  k = -k;	/* reflect */
	    else k += nGridX; 			/* wrap-around */
	  lInterSum -= a[k][j];

	  k = i + interaction - 1;
	  if ( k >= nGridX )  /* beyond right boundary */
	    if ( bReflectiveBoundary )  k = 2*nGridX - k - 2;  /* reflect */
	    else k -= nGridX;			/* wrap-around */
	  lInterSum += a[k][j];
	  c[i][j] = (int)lInterSum;  /* store value for row region */
	  }
	}
      }

    /* now sum horizontal regions from b[] vertical sums */
    for ( i=0;  i<nGridX;  i++ )  {
      lDisperSum = 0;		lInterSum = 0;
      /* do column's upper-most region first */
      for ( j=-dispersal+1;  j<dispersal;  j++ )  {
	k = j;
	if ( k < 0 )	/* beyond top boundary */
	  if ( bReflectiveBoundary ) k = -k;	/* reflect */
	  else k += nGridY;			/* wrap-around */
	lDisperSum += b[i][k];  /* add cell, if present */
	}
      if ( bFullProbab )  lInterSum = lDisperSum;
      else      /* need to sum different sized interaction region */
	for ( j=-interaction+1;  j<interaction;  j++ )  {
	  k = j;
	  if ( k < 0 )  /* beyond top boundary */
	    if ( bReflectiveBoundary ) k = -k;	/* reflect */
	    else k += nGridY;			/* wrap-around */
	  lInterSum += c[i][k];  /* add cell, if present */
	  }
      a[i][0] = cPresence( lDisperSum, lInterSum );

      /* Work down column, subtracting upper-most cell from old Sum */
      /*	and adding bottom-most cell in new region */
      for ( j=1;  j<nGridY;  j++ )  {
	k = j - dispersal;
	if ( k < 0 )		/* beyond top boundary */
	  if ( bReflectiveBoundary )  k = -k;	/* reflect */
	  else k += nGridY;			/* wrap-around */
	lDisperSum -= b[i][k];

	k = j + dispersal - 1;
	if ( k >= nGridY )	/* beyond bottom boundary */
	  if ( bReflectiveBoundary )  k = 2*nGridY - k - 2;  /* reflect */
	  else k -= nGridY;			/* wrap-around */
	lDisperSum += b[i][k];

	if ( bFullProbab)  lInterSum = lDisperSum;
	else  {  /* need to sum different sized interaction region */
	  k = j - interaction;
	  if ( k < 0 )		/* beyond top boundary */
	    if ( bReflectiveBoundary )  k = -k;	/* reflect */
	    else k += nGridY;			/* wrap-around */
	  lInterSum -= c[i][k];

	  k = j + interaction - 1;
	  if ( k >= nGridY )	/* beyond bottom boundary */
	    if ( bReflectiveBoundary ) k = 2*nGridY - k - 2;  /* reflect */
	    else k -= nGridY;			/* wrap-around */
	  lInterSum += c[i][k];
	  }
	a[i][j] = cPresence( lDisperSum, lInterSum );
	}
      }
    }  /* end non-star cases */
  }

/**********************************************************************/
/***                    Random Number Generator                     ***/
/**********************************************************************/
static unsigned long r250_buffer[250];
static int r250_index;

void r250_init(int seed)
{	/* initialize random number generator */
  int j,k;
  unsigned long mask;
  unsigned long msb;
  unsigned long token;
  srand(seed);
  r250_index=0;

  for (j=0;  j<250;  j++)  {
    r250_buffer[j] = rand();
    token = rand();
    token <<= 16;
    r250_buffer[j] ^= token;
    }

  for (j=0;  j<250;  j++)  {
    if (rand() > 16348)  r250_buffer[j]|=0x00008000;
    }

  msb=0x40000000L;
  mask=0x7fffffffL;
  for (j=0;j<31;j++)  {
    k = 7 * j + 3;
    r250_buffer[k] &= mask;
    r250_buffer[k] |= msb;
    mask >>= 1;
    msb >>= 1;
    }
  }


double dr250(void)
{	/* get number from random number generator */
  register int j;
  register unsigned long new_rand;

  if (r250_index >= 147)  j=r250_index - 147;
  else  j=r250_index + 103;

  new_rand = r250_buffer[r250_index] ^ r250_buffer[j];  /* bitwise XOR */
  r250_buffer[r250_index] = new_rand;
  if (r250_index>=249)  r250_index=0;
  else r250_index++;

  return new_rand/(double)0x7fffffffL;
  }



/**********************************************************************/
/***                   Measure of Clustering                        ***/
/**********************************************************************/

void cluster(int *aiStats)
{
  int i, j;
  long  cluster = 0,
	occ = 0,
	lGridXY;

  lGridXY = nGridX * (long)nGridY;

  /***  calculate clustering statistics    ***/
  /* count number of cells equal to 1 */
  for (i = 0;  i<nGridX;  i++)
    for (j = 0;  j<nGridY;  j++)
      occ += a[i][j];	/* add value of 1 or 0 */
  aiStats[1] = (int)((occ * 100) / lGridXY);

  /* count number of horizontally adjacent cells that are the same */
  for (i = 0;  i<nGridX;  i++)
    for (j = 1;  j<nGridY;  j++)
      cluster += (a[i][j] == a[i][j-1]);  /* adds 1 (true) if they're equal */
  aiStats[0] = (int)((cluster * 100) / lGridXY);
  }


/**********************************************************************/
/***                      Output Binary File                        ***/
/**********************************************************************/

void output(int *aiStats, int nReplicates)
{
  /* save state of system at each time step   */
  #ifdef PC
    if ( nReplicates <= 1 )	/* print column headings for output file */
      fprintf( outptr1,
	    "grid inter displ  init blk alpha  seed  time  clust  occ\n" );

    fprintf(outptr1,
	  "%4d %5ld %5ld %5.2f %3d %5.2f %5d %5ld %5d %5d\n",
	  nGridY, lInteractRegion,lDispRegion,propA,bBlockFormStart,fAlpha,
	  iReplicateSeed,lTime,aiStats[0],aiStats[1]);
    /* "ir= %5.2f dr= %3.2f initial= %3.2f bBlockFormStart= %d alpha= %3.2f seed= %d time= %d cluster= %d occ=%d\n", */
  #endif

  #ifdef MF
    char *pstate, *state_ptr;
    short i, j;
    long grid;

    /***  prepare data for output in binary format  ***/
    grid = nGridX * nGridY;
    pstate = (char *) malloc(grid);
    state_ptr = pstate;
    for (i=0;  i<nGridX;  i++)
      for (j=0;  j<nGridY;  j++)  {
	*state_ptr = a[i][j];
	state_ptr++;
	}

    fwrite(&lTime, sizeof(short), 1, outptr1);
    fwrite(pstate, sizeof(char), grid, outptr1);
  #endif
  }

/**********************************************************************/
/***                        Setvmode for PC                         ***/
/**********************************************************************/

void setvmode(char mode)
{
  union REGS regs;
  regs.h.ah = 0;
  regs.h.al = mode;
  int86(0x10, &regs, &regs);
  }


/**********************************************************************/
/***                        graph the grid                          ***/
/**********************************************************************/

void graphit(int *aiStats)
{
  char  S[19];
  int i; /* */

  /*** write time step to the screen  ***/
  sprintf( S, "%ld", lTime );	    wrstr(S, 12, 200, 71);

  /*** write % clustered to the screen  ***/
  sprintf( S, "%2d", aiStats[0] );  wrstr(S, 12, 250, 71);

  /*** write % occupied to the screen  ***/
  sprintf( S, "%2d", aiStats[1] );  wrstr(S, 12, 275, 71);

  for (i=0;  i<nGridX;  i++)
    PaintRow( i, 9 ); /* */
  }


/**********************************************************************/
/***                      PaintRow for PC                           ***/
/**********************************************************************/

void PaintRow( short iRow, unsigned int wColor )
{	/* Paint row from array a[row][0] */
  static unsigned char
	acBuf[80],	/* stores copy of compressed screen bytes */
	*pacBuf;	/* ptr to above */
  int j,k;
  unsigned char C;
  unsigned char far *pcArow,	/* points within row of A[iRow][] */
		    *pcEga;	/* points to row of 640x480 EGA screen */

  /* Set Bit Mask to 0xFF */
  outport( 0x03ce, 0xFF08);     /* enable all mask bits in high byte, "8" is bit mask index */
  outport( 0x03ce, 0x0003 );	/* select EGA write mode 0, "3" is index for Rotate/Mode reg */
  outport( 0x03ce, 0x0001 );	/* CPU byte maps to pixels, "1" is index for Enable Set/Reset reg */
  outport( 0x03C4, (wColor << 8) + 2 );	/* enable "blue" planes (9), "2" is index for Color Map Mask reg */

  pcEga = egaptr + 80 * iRow * iPaintBits;     /* calc screen byte offset */
  pcArow = a[iRow];	/* set ptr to row a[iRow][0] */
  pacBuf = acBuf;	/* set ptr to acBuf */

  switch ( iPaintBits )  {
    case 1: /* Each byte of 0 or 1 maps to one screen pixel. */
	/* Compress 8 bytes to 8 bits and write directly screen. */
	for ( j=0;  j<nGridY;  j+=8 )  {
	   C = *(pcArow++);	/* assign a[] row value, incr row ptr */
	   for ( k=1;  k<8;  k++ )
	     C = C*2 + *(pcArow++);	/* add row value, advance row ptr */
	   *(pcEga++) = C;	/* store value to screen, advance ptr */
	   }
	break;		/* finished writing to screen */

    case 2: /* Each byte of 0/1 maps to two screen pixels. */
	/* Compress 4 bytes to 8 bits and store in acScreenRow[] for double write. */
	/* Sub 1 from each source byte to make 0xFF or 0, then mask with 3 to get inversed bits. */
	for ( j=0;  j<nGridY;  j+=4 )  {
	   C = ( *(pcArow++) - 1 );  /* sub 1 for inverse, advance ptr */
	   for ( k=1;  k<4;  k++ )
	     C = C*4 + ((*(pcArow++) - 1) & 3);	/* sub 1 for inverse, mask, add value, advance row ptr */
	   *(pacBuf++) = ~C;  /* complement to restore value, store in buffer array, advance ptr */
	   }
	for ( j=0;  j<2;  j++ )  {	/* copy buffer to screen */
	    memmove( pcEga, acBuf, (nGridY+3) / 4);
	    pcEga += 80;  /* advance ptr to next row */
	    }
	break;
    case 4: /* Each byte of 0/1 maps to four screen pixels */
	/* Compress 2 bytes to 8 bits and store in acScreenRow[] for quad write. */
	/* Sub 1 from each source byte to make 0xFF or 0, then mask with 0Fh to get inversed bits. */
	for ( j=0;  j<nGridY;  j+=2 )  {
	   C = ( *(pcArow++) - 1 );  /* sub 1 for inverse, advance ptr */
	   C = C*16 + ((*(pcArow++) - 1) & 0x0F); /* sub 1 for inverse, mask, add value, advance row ptr */
	   *(pacBuf++) = ~C;  /* complement to restore value, store in buffer array, advance ptr */
	   }
	for ( j=0;  j<4;  j++ )  {	/* copy buffer to screen */
	    memmove( pcEga, acBuf, (nGridY+1) / 2);
	    pcEga += 80;  /* advance ptr to next row */
	    }
	break;
    }

  outport( 0x03C4, 0xFF02 );	/* enable all color planes (FF), "2" is index for Color Map Mask reg */
  }


/**********************************************************************/
/***                    Write a Character                           ***/
/**********************************************************************/

void wrch(char b, char c, int x, int y)
{
  int i;
  unsigned char far *ptr;
  unsigned char far *chptr;
  unsigned char ch_size;

  ch_size = *( (char far *)(0x00400085L) ); /* current font size */

  /***  first erase what is there now  ***/

  /* reset data rotate register so new data replaces old */
  outportb(0x03ce, 0x03);
  outportb(0x03cf, 0x00);

  /* Put color black in set/reset register */
  outportb(0x03ce, 0x00);
  outportb(0x03cf, 0x00);

  /* enable set rest register */
  outportb(0x03ce, 0x01);
  outportb(0x03cf, 0x0f);

  /* reset bit mask register */
  outportb(0x03ce, 0x08);
  outportb(0x03cf, 0xff);

  ptr = egaptr + 80*x + y;
  for (i = 0; i < ch_size; i++)  {
    *ptr &= 0xff;
    ptr += 80;
    }

  /* Put color c in set/reset register */
  outportb(0x03ce, 0);
  outportb(0x03cf, c);

  /* enable set rest register */
  outportb(0x03ce, 0x01);
  outportb(0x03cf, 0x0f);

  chptr = (unsigned char far *) (*(INT43) + (b*ch_size));
  ptr = egaptr + 80*x + y;
  for (i = 0; i < ch_size; i++)  {
    /* reset bit mask register with next line of character */
    outportb(0x03ce, 8);
    outportb(0x03cf, *(chptr++));
    *ptr &= 0xff;
    ptr  += 80;
    }
  }


/**********************************************************************/
/***                    Write a Phrase                              ***/
/**********************************************************************/

void wrstr(char phrase[], char c, int  x, int y)
{
  int  i;

  for (i=0; phrase[i] != '\0'; i++)
    wrch(phrase[i], c, x, y+i);
  }