/************************************************************************
 *This program is a modification of model8.c to include the subroutine  *                                                                      
 *rrandom.c.  Random.c generates random integers from 0 to the argument *
 *Thus, the call to random in modiy8.c has been changed to              *
 *rrandom(idim-1).  See line 367.                                       *
 *The random number generator has been changed to ran1.c                *
 *Changed input to read x and y only.  Other parameters #defined        *
 *Changed output to append summary only to argv[2]                      *
 *                                                                      *
 *      FILE            MODEL84.C                                       *
 *      PROGRAM         CIRCLE GROWTH MODEL FOR PLANTS                  *
 *      DATE            August 1996                                     *
 *      AUTHORS         BILL PRICE, ROBERT DUBEAU                       *
 *      LOCATION        UNIVERSITY OF IDAHO AGRICULTURE                 *
 *      HARDWARE        IBM PC COMPATIBLE                               *
 *      SOFTWARE        WATCOM C COMPILER.                              *
 *      MODIFICATIONS   H. J. PRICE                                     *                                                                      *
 *      DESCRIPTION:    Reads plant data from input file and            *
 *        simulates growth and competition among plants according       *
 *        to a circle growth model.  Plant resource usage is            *
 *        stored in a dynamic array, and after completion               *
 *        of the simulation, final plant data are written to output     *
 *        file.  Input and output file pathnames must be specified in   *
 *        the command line.                                             *
 *                                                                      *
 *      USAGE:          model8 infile_pathname outfile_pathname         *
 *                                                                      *
 *      CONSTANTS:                                                      *
 *                                                                      *
 *  ESC                 escape key code                                 *
 *  b0                  constant in the growth equation                 *
 *  b1                  constant in the growth equation                 *
 *  MIN_NEEDED          percent of requested resources necessary to     *
 *                        maintain a plant's health                     *
 *  RESOURCE_MIN        percent of resources NOT available to plants.   *
 *                        can be set higher than 0 to simulate resource *
 *                        scarcity                                      *
 *  FIELD_X             field size x                                    *
 *  FIELD_Y             field size y                                    *
 *                                                                      *
 *  enum {              indices to array of plant characteristics       *
 *        X               x coordinate of plant                         *
 *        Y               y coordinate of plant                         *
 *        MAX_R           maximum attainable radius of plant            *
 *        ALPHA           intensity index for resource uptake           *
 *        CUR_R           current plant radius                          *
 *        HEALTH          current plant health (decremented if percent  *
 *                          resource uptake is below MIN_NEEDED)        *
 *        PCT             percent of requested resources received       *
 *        CUM_PCT         cumulative percent of requested resources     *
 *                        received                                      *
 *        NUM_CHARS       number of different plant characteristics     *
 *                          for each plant                              *
 *       }                                                              *
 *                                                                      *
 ************************************************************************/

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

#define ESC 27

#define b0 1.955244
#define b1 0.103715

#define MIN_NEEDED 33
#define RESOURCE_MIN 0
#define FIELD_X 525
#define FIELD_Y 525

#define BIG_R 30
#define SMALL_r 2
#define ALPH 0.06
#define HLTH 51

enum {X, Y, MAX_R, ALPHA, CUR_R, HEALTH, PCT, CUM_PCT, NUM_CHARS};

typedef unsigned char byte;     /* for small int data storage */

byte **dist;      /* array of distance values from center of a plant */
byte **plt_vals;  /* array of plant uptake values indexed by alpha (0..10) */
byte **field;      /* array of plant locations */

byte **create_byte_array(size_t idim, size_t jdim);     /* array of zeros */
int **create_int_array(int idim, int jdim);     /* array of zeros */
int **copy_int_ptrs(int **array, int idim);     /* for after shuffle */

void shuffle(int **array, int idim);    /* shuffles array pointers */
void pause(void); /* check for user interrupt */
void iterate_plant(int **plant_data, byte **field, int index);
        /* biggest time-user. does the pixel calculations for each plant */

int file_size(char *filename);
void read_file(char *filename, int **matrix, int num, int *max_plantsize);
void write_file(char *filename, int **matrix, int num);

/************************************************************************
 *                                                                      *
 *      main()                                                          *
 *                                                                      *
 *      DESCRIPTION:    Calls data input from file; initializes byte    *
 *        tables for iterate_plant; runs simulation; calls file output  *
 *                                                                      *
 ************************************************************************/
void main(int argc, char* argv[])
{
  int **plant_data;             /* array of plant characteristics */
  int **plant_data_save;        /* saves original order ofplant_data */
/*  byte **field;   */              /* this is where the plants live */
  int plant_num = 0;            /* number of plants found in infile */
  int cycle;                    /* main loop index */
  int plant;                    /* current plant in cycle */
  int i, j;                     /* index counters */
  int max_plantsize;            /* largest plant radius found by read_file */
  int r;                        /* radius beyond max_plantsize for dist */
  int rrandom();                 /* random returns a random integer */

        /***** exit if I/O not specified in command line *****/
/*  printf("start of program\n");*/
  if (argc != 3) {
    printf("Usage: %s infile outfile\n", argv[0]);
    exit(1);
  }

        /***** get input data *****/
  plant_num = file_size(argv[1]);  /* get number of plants in infile */
  printf("\nThe file size is: %d\n", plant_num);

  plant_data = create_int_array(plant_num, NUM_CHARS); /* for read_file */
  plant_data_save = copy_int_ptrs(plant_data, plant_num);

        /* retrieves order after shuffling */
  read_file(argv[1], plant_data, plant_num, &max_plantsize);
  /* print max_plantsize */
  printf("max_plantsize is %5d\n",max_plantsize);
  printf("plant_num is %5d\n",plant_num);

        /* infile data goes into plant_data, max_plantsize is returned */
        /* ***** if data out of bounds, read_file will exit program ****/

        /***** create array of distances for iterate_plant() *****/

  r = max_plantsize + 2;  /* must leave room beyond plant boundary */
  dist = create_byte_array(r, r);
  for (i=0; i<r; i++)
    for (j=0; j<r; j++)
      dist[i][j] = (byte) (sqrt((i*i) + (j*j))+.5);  /* round decimals */

        /***** create array of exp(-alpha d) for iterate_plant *****/

  plt_vals = create_byte_array(10, max_plantsize+2);  /* 10 alpha values */

  for (i=0; i<10; i++)  /* i is 100 alpha (0.00 to 0.09) */
    for (j=0; j<max_plantsize; j++){  /* j is just d */
      plt_vals[i][j] = (byte) (exp(-i/100.0*j)*9 + 1);
    }

        /***** create array of field elements for plants to live in *****/

  field = create_byte_array((int) FIELD_X,(int) FIELD_Y);

        /***** loop through 50 cycles times plant_num plants *****/

  for (cycle = 1; cycle <= 50; ++cycle){ 
        
    /* reset the field */

    for(i=0;i<FIELD_X;i++){
       for(j=0;j<FIELD_Y;j++){
         field[i][j]=(byte)(10);
         }
    }

    shuffle(plant_data, plant_num);  /* reorder plants for each iteration */
    for (plant=0; plant<plant_num; plant++) /* step through all plants */
      if (plant_data[plant][HEALTH] > 0) /* if plant is still alive */
        iterate_plant(plant_data, field, plant); /* feed and grow the plant */
  } /* for cycle */
  printf("\n");

        /***** write final plant data to output file *****/

  write_file(argv[2], plant_data_save, plant_num);
} /* main() */

/************************************************************************
 *                                                                      *
 *      file_size()                                                     *
 *                                                                      *
 *      DESCRIPTION:    opens infile; counts the number of data lines   *
 *        by counting carriage returns in the file; closes infile.      *
 *                                                                      *
 ************************************************************************/
int file_size(char *filename)
{
  int size=0;
  char c;
  FILE *infile;

  infile = fopen(filename, "r");
  while((c=getc(infile)) != EOF) {
    if (c == '\n'){
      ++size;
    }
  }
  fclose(infile);
  return size;
}

/************************************************************************
 *                                                                      *
 *      read_file()                                                     *
 *                                                                      *
 *      DESCRIPTION:    opens file; reads data and stores in array;     *
 *        closes file; returns largest plant size in *max_plantsize;    *
 *        exits program if plant data will not fit in field array       *
 *                                                                      *
 ************************************************************************/
void read_file(char *filename, int **m, int num, int *max_plantsize)
{
  int plant;
  int x, y, R = BIG_R;
  int min_x = 1000, min_y = 1000, max_x = 0, max_y = 0;
  float alpha;
  FILE *infile;

  infile = fopen(filename, "r");

  *max_plantsize = 0;

  for (plant=0; plant<num; ++plant) {
    fscanf(infile, "%d %d",
                &x, &y);
    m[plant][X] = x - 225;
    m[plant][Y] = y - 75;
    m[plant][MAX_R] = BIG_R;
    m[plant][ALPHA] = 100*ALPH+.5;
    m[plant][CUR_R] = SMALL_r;
    m[plant][HEALTH] = HLTH;

    if (R > *max_plantsize) *max_plantsize = R;
    if (x - 225 < min_x) min_x = x - 225;
    if (y - 75 < min_y) min_y = y - 75;
    if (x - 225 > max_x) max_x = x - 225;
    if (y - 75 > max_y) max_y = y - 75;
  }
  fclose(infile);
  if ((min_x < 25) || (min_y < 25) || (max_x > 500) || (max_y > 500)) {
    printf(" min_x: %d\n min_y: %d\n max_x: %d\n max_y: %d\n",
              min_x, min_y, max_x, max_y);
    printf(" max_plantsize: %d\n", *max_plantsize);
    printf("Plant data for this file is out of bounds for the\n");
    printf("field size of this program.\n");
    exit(1);
  }
}

/************************************************************************
 *                                                                      *
 *      write_file()                                                    *
 *                                                                      *
 *      DESCRIPTION:    opens output file; writes array to file;        *
 *        closes file                                                   *
 *                                                                      *
 ************************************************************************/
void write_file(char *filename, int **m, int num)
{
  int i;
  int one_m_n=0;
  float cumm=0.0;
  float avg_bio;

  FILE *outfile;

  outfile = fopen(filename, "a");

  for (i=0; i<num; i++){
    if ((m[i][0]>=90) && (m[i][1]>=85) && 
       (m[i][0]<=190) && (m[i][1]<=185)){
        cumm = cumm + (float)(m[i][4]*m[i][4]);
        one_m_n = one_m_n + 1;
    }
  }
  avg_bio = (cumm/(float)one_m_n);

    fprintf(outfile, "%3d  %6.2f  %6.2f\n", one_m_n, cumm, avg_bio);

  fclose(outfile);
}

/************************************************************************
 *                                                                      *
 *      create_byte_array()                                             *
 *                                                                      *
 *      DESCRIPTION:    Uses calloc() to allocate storage for idim by   *
 *        jdim array of byte and a vector of pointers to the array;     *
 *        initializes pointers to first elements of rows in the array;  *
 *        returns array pointer                                         *
 *                                                                      *
 ************************************************************************/
byte **create_byte_array(size_t idim, size_t jdim)
{
  byte **array;
  size_t i;
  static int j;

  array = (byte **) calloc(idim, sizeof(byte *));
  /* printf("Size of array is %u\n",_msize(array)); */
  *array = (byte *) calloc(idim * jdim,sizeof(byte));
  /* printf("Size of *array is %u\n",_msize(*array)); */

  for (i=0; i<idim; i++)
    *(array+i) = *array + jdim*i;

  return array;
}

/************************************************************************
 *                                                                      *
 *      create_int_array()                                              *
 *                                                                      *
 *      DESCRIPTION:    Uses calloc() to allocate storage for idim by   *
 *        jdim array of int and a vector of pointers to the array;      *
 *        initializes pointers to first elements of rows in the array;  *
 *        returns array pointer                                         *
 *                                                                      *
 ************************************************************************/
int **create_int_array(int idim, int jdim)
{
  int **array;
  int i;

  array = (int **) calloc(idim, sizeof(int *));
  *array = (int *) calloc(idim * jdim, sizeof(int));

  for (i=0; i<idim; i++)
    *(array+i) = *array + jdim*i;

  return array;
}

/************************************************************************
 *                                                                      *
 *      copy_int_pointers()                                             *
 *                                                                      *
 *      DESCRIPTION:    allocates space for and copies only the vector  *
 *        of pointers to int array; returns pointer to the copy         *
 *                                                                      *
 ************************************************************************/
int **copy_int_ptrs(int **array, int idim)
{
  int **copy;
  int i;

  copy = (int **) calloc(idim, sizeof(int *));

  for (i=0; i<idim; i++)
    *(copy+i) = *(array+i);

  return copy;
}

/************************************************************************
 *                                                                      *
 *      shuffle()                                                       *
 *                                                                      *
 *      DESCRIPTION:    reorder the array randomly by shuffling the     *
 *        pointers in *array.  permanently reorders.  copy the pointers *
 *        before using this                                             *
 *                                                                      *
 ************************************************************************/
void shuffle(int **array, int idim)
{
  int rnum, *temp;
  int num = idim;
  /*
  randomize();
  */

  while(--num >= 0) {
    rnum = rrandom(idim-1);
    temp = *(array+num);
    *(array+num) = *(array+rnum);
    *(array+rnum) = temp;
  }
}


/************************************************************************
 *                                                                      *
 *      iterate_plant()                                                 *
 *                                                                      *
 *      DESCRIPTION:    feeds and grows a plant.  somewhat repetitive   *
 *        code is designed for speed, not conciseness.  calculations    *
 *        on field elements are done by integer table lookup.           *
 *                                                                      *
 ************************************************************************/
void iterate_plant(int **plant_data, byte **field, int index)
{
  extern byte **dist;  /* array of calculated distances within a plant */
  extern byte **plt_vals;  /* array of calculated resource uptake */
  int cur_r;  /* from **plant_data */
  int distance, plt_value;  /* looked up in **dist and **plt_vals */
  long actual_tot=0;  /* actual total food intake (note long int) */
  long potent_tot=0;  /* potential total food intake (note long int) */
  int cell_value;  /* amount of resource left in each cell (pixel) */
  int pct_got;  /* percent of needed resources actually gotten */
  int i, j;  /* indices */
  int x, y;  /* from **plant_data */
  byte *pvals;  /* pointer to intake values for given alpha */

  pvals = plt_vals[plant_data[index][ALPHA]];
        /* get pointer to plt_vals for this plant's alpha */

  x = plant_data[index][X];
  y = plant_data[index][Y];
  cur_r = plant_data[index][CUR_R];

/***********  pixel calculations  **************/

        /***** do center pixel *****/

  plt_value = pvals[0];  /* resources needed at zero distance */
  actual_tot += plt_value;  /* pretend they're all here */
  potent_tot += plt_value;  /* store needed amount */
  cell_value = field[x][y] - plt_value; /* allocate resource to plant */
  if (cell_value < RESOURCE_MIN) { /* ran out of resource */
    actual_tot -= (RESOURCE_MIN - cell_value); /* take some back */
    cell_value = RESOURCE_MIN; /* and reset cell to minimum */
  }
  field[x][y] = cell_value; /* put remaining resource back in field */

        /***** do horizontal midline *****/

  for (i=1; i <= cur_r; i++) { /* repeat for left and right of center */
    plt_value = pvals[i];
    actual_tot += 2*plt_value;
    potent_tot += 2*plt_value;

    cell_value = field[x+i][y] - plt_value;
    if (cell_value < RESOURCE_MIN) {
      actual_tot -= (RESOURCE_MIN - cell_value);
      cell_value = RESOURCE_MIN;
    }
    field[x+i][y] = cell_value;

    cell_value = field[x-i][y] - plt_value;
    if (cell_value < RESOURCE_MIN) {
      actual_tot -= (RESOURCE_MIN - cell_value);
      cell_value = RESOURCE_MIN;
    }
    field[x-i][y] = cell_value;
  }

        /***** do vertical midline *****/

  for (j=1; j <= cur_r; j++) { /* repeat for above and below center */
    plt_value = pvals[j];
    actual_tot += 2*plt_value;
    potent_tot += 2*plt_value;

    cell_value = field[x][y+j] - plt_value;
    if (cell_value < RESOURCE_MIN) {
      actual_tot -= (RESOURCE_MIN - cell_value);
      cell_value = RESOURCE_MIN;
    }
    field[x][y+j] = cell_value;

    cell_value = field[x][y-j] - plt_value;
    if (cell_value < RESOURCE_MIN) {
      actual_tot -= (RESOURCE_MIN - cell_value);
      cell_value = RESOURCE_MIN;
    }
    field[x][y-j] = cell_value;
  }

        /***** do non-central pixels *****/

  for (i=1; dist[i][0] <= cur_r; i++) { /* repeat now for four quadrants */
    for (j=1; (distance = dist[i][j]) <= cur_r; j++) {

      plt_value = pvals[distance];
      actual_tot += 4*plt_value;
      potent_tot += 4*plt_value;

      cell_value = field[x+i][y+j] - plt_value;
      if (cell_value < RESOURCE_MIN) {
        actual_tot -= (RESOURCE_MIN - cell_value);
        cell_value = RESOURCE_MIN;
      }
      field[x+i][y+j] = cell_value;

      cell_value = field[x-i][y+j] - plt_value;
      if (cell_value < RESOURCE_MIN) {
        actual_tot -= (RESOURCE_MIN - cell_value);
        cell_value = RESOURCE_MIN;
      }
      field[x-i][y+j] = cell_value;

      cell_value = field[x+i][y-j] - plt_value;
      if (cell_value < RESOURCE_MIN) {
        actual_tot -= (RESOURCE_MIN - cell_value);
        cell_value = RESOURCE_MIN;
      }
      field[x+i][y-j] = cell_value;

      cell_value = field[x-i][y-j] - plt_value;
      if (cell_value < RESOURCE_MIN) {
        actual_tot -= (RESOURCE_MIN - cell_value);
        cell_value = RESOURCE_MIN;
      }
      field[x-i][y-j] = cell_value;
      }
    }

        /***** grow the plant *****/

    pct_got = ((1000*actual_tot)/potent_tot + 5)/10;
        /* gives rounded integer percent */
    plant_data[index][PCT] = pct_got; /* for output file */
    plant_data[index][CUM_PCT] += pct_got; /* for output */
    plant_data[index][CUR_R] = /* the growth function */
                1.0/(1.0+exp(b0-b1*plant_data[index][CUM_PCT]/100.0))
                *(double)plant_data[index][MAX_R] + 0.5;

    if (plant_data[index][CUR_R] < 1) /* maintenance stuff */
      plant_data[index][CUR_R] = 1;
    if (pct_got <= MIN_NEEDED)
      plant_data[index][HEALTH] -= 1;
    if (plant_data[index][HEALTH] < 0)
      plant_data[index][HEALTH] = 0;
}

/***************************************************************** 
This routine returns a random integer between 0 and imax.  
******************************************************************/

#include <stdlib.h>
#include <math.h>

int rrandom(int imax)
{
     static long idum = (-1);
     float ran1(),r;
     int imax_r;

     r=ran1(&idum);
     imax_r=imax*r;

     return(((imax*r-imax_r)>0.5)?(imax_r+1):imax_r);

}/* end rrandom() */
/************************************************************* 
This is the random number generator ran1() from "Numerical
Recipes in C"
**************************************************************/

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran1(long *idum)
/* Minimum random number generator of Park an Miller with Bays-Durham
shuffle and added safeguards.  Returns uniform random deviate between
0.0 and 1.0 (exclusive of the endpoint values).  Call with idum a
negative integer to initialize; thereafter, do not alter idum between
successive deviates in a sequence.  RNMX should approximate the largest
floating value that is less than 1. */
{
     int j;
     long k;
     static long iy=0;
     static long iv[NTAB];
     float temp;

     if(*idum <= 0 || !iy){
        if(-(*idum) < 1) *idum=1;
        else *idum = -(*idum);
        for(j = NTAB+7;j >= 0;j--){
           k = (*idum)/IQ;
           *idum = IA*(*idum - k*IQ) - IR*k;
           if(*idum < 0) *idum += IM;
           if(j < NTAB) iv[j] = *idum;
        }
        iy = iv[0];
     }
     k = (*idum)/IQ; /* Start here when not initializing */
     *idum = IA*(*idum - k*IQ) - IR*k; /* Compute idum = (IA*idum) % IM
                                          without overflows by Schrage's
                                          method. */
     if(*idum < 0) *idum += IM;
     j = iy/NDIV; /* Will be in the range 0..NTAB-1. */
     iy = iv[j]; /* Output previously stored value and refill the
                    shuffle table. */
     iv[j] = *idum;
     if((temp = AM*iy) > RNMX) return RNMX; /* Because users don't expect
                                               endpoint values. */
     else return temp;
} /* end ran1 */