/***********************************************************
 *
 *  pgmglacier.c - a heightfield glaciation tool
 *
 *  Mark J. Stock, September 24, 1998, version 0.2
 *
 *  0.0  1998-09-03  Initial program
 *  0.1  1998-09-03  Improved ice flow algorithm
 *  0.2  1998-09-24  Added default writing to stdout
 *
 ***********************************************************/


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "const.h"

/* The heightfield array, values between 0k and 65k */
unsigned short int hf[MAX_SIZE][MAX_SIZE];

/* The ice/glacier thickness */
unsigned short int gt[MAX_SIZE][MAX_SIZE];

int inheight, inwidth;		/* size of the input file */

void main(int argc,char **argv) {

   int i,j,step,temp_i;
   int interval = 50;
   int upgt,uphf,dngt,dnhf,ltgt,lthf,rtgt,rthf;
   int num_steps = 200;
   int use_alt = 0;
   int zero_alt = 500;
   double temp_d;
   double rate = 0.2;		/* rate of deposition/sublimation, 1.0 for 250x250 file */
   char infile[80];
   char tempfile[80];
   char outfile[80] = "-";	/* unless reset, write to stdout */
   char progname[80];

   inheight = 256;		/* default values */
   inwidth = 256;

   /*
    * Parse command-line args
    */
   (void) strcpy(progname,argv[0]);
   if (argc < 2) (void) Usage(progname,0);
   if (strncmp(argv[1], "-help", 2) == 0)
      (void) Usage(progname,0);
   (void) strcpy(infile,argv[1]);
   for (i=2; i<argc; i++) {
      if (strncmp(argv[i], "-help", 2) == 0)
         (void) Usage(progname,0);
      else if (strncmp(argv[i], "-n", 2) == 0)
         num_steps = atoi(argv[++i]); 
      else if (strncmp(argv[i], "-a", 2) == 0)
         use_alt = 1;
      else if (strncmp(argv[i], "-outfile", 2) == 0)
         (void) strcpy(outfile,argv[++i]); 
      else
         (void) Usage(progname,0);
   }
   fprintf(stderr,"Writing file to %s\n",outfile);


   /* read in the input heightfield file */
   (void)read_pgm_ushort(infile,hf);

   /* Set initial ice thickness by altitude or evenly */
   if (!use_alt) {
      for (i=0; i<inheight; i++) {
         for (j=0;j<inwidth;j++) {
            gt[i][j] = 20;
         }
      }
   } else {
      for (i=0; i<inheight; i++) {
         for (j=0;j<inwidth;j++) {
            if (hf[i][j] > 7000)
               gt[i][j] = 100;
            else if (hf[i][j] > 3000)
               gt[i][j] = (hf[i][j]-3000)/40;
            else
               gt[i][j] = 0;
         }
      }
   }

   /* then, iterate through steps */
   fprintf(stderr,"Stepping through the simulation ");
   for (step=0; step<num_steps; step++) {

      fprintf(stderr,".");
      fflush(stderr);

      /* add/remove snow based on altitude */
      for (i=0;i<inheight;i++) {
         for (j=0;j<inwidth;j++) {
            /* temp_i = gt[i][j] + hf[i][j]/1000 - 3; */
            /* if (temp_i > 0) gt[i][j] = temp_i; */
            temp_d = rate*(hf[i][j]-zero_alt)/1000;
            if (temp_d < 0.0) {
               if (((int)(temp_d) - temp_d) > rand()/(RAND_MAX+1.0)) temp_i = (int)temp_d - 1;
               else temp_i = (int)temp_d;
            } else {
               if ((temp_d - (int)temp_d) > rand()/(RAND_MAX+1.0)) temp_i = (int)temp_d + 1;
               else temp_i = (int)temp_d;
            }
            temp_i += gt[i][j];
            if (temp_i < 0) gt[i][j] = 0;
            else gt[i][j] = temp_i;
         }
      }

      /* calculate new values from old values */
      for (i=0;i<inheight;i++) {
         for (j=0;j<inwidth;j++) {
            if (i==0) {
               upgt = gt[0][j];
               uphf = hf[0][j];
               dngt = gt[1][j];
               dnhf = hf[1][j];
            } else if (i==inheight-1) {
               upgt = gt[inheight-2][j];
               uphf = hf[inheight-2][j];
               dngt = gt[inheight-1][j];
               dnhf = hf[inheight-1][j];
            } else {
               upgt = gt[i-1][j];
               uphf = hf[i-1][j];
               dngt = gt[i+1][j];
               dnhf = hf[i+1][j];
            }
            if (j==0) {
               ltgt = gt[i][0];
               lthf = hf[i][0];
               rtgt = gt[i][1];
               rthf = hf[i][1];
            } else if (j==inwidth-1) {
               ltgt = gt[i][inwidth-2];
               lthf = hf[i][inwidth-2];
               rtgt = gt[i][inwidth-1];
               rthf = hf[i][inwidth-1];
            } else {
               ltgt = gt[i][j-1];
               lthf = hf[i][j-1];
               rtgt = gt[i][j+1];
               rthf = hf[i][j+1];
            }
            temp_i = (upgt+uphf+dngt+dnhf+ltgt+lthf+rtgt+rthf)/4;
            if (temp_i > 65535) {
               fprintf(stderr,"Value too high for short int, exiting.\n");
               fprintf(stderr,"i= %d, j= %d, temp_i= %d\n\n",i,j,temp_i);
               exit(0);
            } else if (temp_i > hf[i][j]) gt[i][j] = temp_i-hf[i][j];
            else gt[i][j] = 0;
         }
      }

      /* periodically, write a temporary gt file */
      if (step/interval == (step+interval-1)/interval && step != 0) {
         fprintf(stderr,"gt%03d.pgm",step);
         fflush(stderr);
         sprintf(tempfile,"/tmp/gt%03d.pgm",step);
         (void) write_pgm_ushort(tempfile,gt,inheight,inwidth);
      }

   }
   fprintf(stderr,"\nDone.\n");

   /* write out a frame at the end */

   /* Add the glacier thickness to the elevation file */
   for (i=0; i<inheight; i++) {
      for (j=0;j<inwidth;j++) {
         /* hf[i][j] += gt[i][j]; */
      }
   }

   /* write the total elevation (land and ice) */
   /* (void) write_pgm_ushort(stdout,hf,inheight,inwidth); */

   /* write the glacier thickness file */
   (void) write_pgm_ushort(outfile,gt,inheight,inwidth);

}


/*
 * This function writes basic usage information to stderr,
 * and then quits. Too bad.
 */
int Usage(char progname[80],int status) {

   /* Usage for pgmglacier */
   static char **cpp, *help_message[] =
   {
       "where options include:",
       "  -help          writes this help information",
       "  -outfile name  glacier thickness output file (ASCII PGM)",
       "  -n val         number of cycles, def= 50",
       "  -a             vary the initial ice level by elevation",
       " ",
       "Options may be abbreviated to an unambiguous length.",
       NULL
   };

   fprintf(stderr, "usage:\n  %s infile [-options ...]\n\n", progname);
   for (cpp = help_message; *cpp; cpp++)
      fprintf(stderr, "%s\n", *cpp);
   exit(status);
   return(0);
}

