/****************************************************************************
 *
 * MODULE:       r.simwe.sediment: main program for sediment transport
 *               simulation (SIMWE)
 *
 * AUTHOR(S):    L. Mitas,  H. Mitasova, J. Hofierka
 * PURPOSE:      Sediment transport simulation (SIMWE)
 *
 * COPYRIGHT:    (C) 2002 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *
 *****************************************************************************/

/*-
 * r.simwe.sediment: main program for hydrologic and sediment transport
 * simulation (SIMWE)
 *
 * Original program (2002) and various modifications:
 * Lubos Mitas, Helena Mitasova
 *
 * GRASS5.0 version of the program:
 * J. Hofierka
 *
 * Copyright (C) 2002 L. Mitas,  H. Mitasova, J. Hofierka
 *
 *This program is free software; you can redistribute it and/or
 *modify it under the terms of the GNU General Public License
 *as published by the Free Software Foundation; either version 2
 *of the License, or (at your option) any later version.
 *
 *This program is distributed in the hope that it will be useful,
 *but WITHOUT ANY WARRANTY; without even the implied warranty of
 *MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *GNU General Public License for more details.
 *
 *You should have received a copy of the GNU General Public License
 *along with this program; if not, write to the Free Software
 *Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA
 *
 *
 * Notes on modifications:
 * v. 1.0 June 2002
 *
 */

#define NWALK	"2000000"
#define DIFFC	"0.8"
#define NITER   "1200"
#define ITEROUT "300"
#define DENSITY "200"

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

#include <grass/config.h>
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/linkm.h>
#include <grass/bitmap.h>
#include <grass/site.h>
#include <grass/glocale.h>

#define MAIN
#include <grass/waterglobs.h>

char fncdsm[32];
char filnam[10];

/*struct BM *bitmask;*/
/*struct Cell_head cellhd;*/
struct GModule *module;
struct Map_info Map;

char msg[1024];

int main ( int argc, char *argv[]) 
{
  int i,ii,j,l;
  int ret_val;
  double x_orig, y_orig;
  static int rand1 = 12345;
  static int rand2 = 67891;

  G_gisinit (argv[0]);

  module = G_define_module();
  module->keywords = _("raster, sediment flow, erosion, deposition");
  module->description =        
	_("Sediment transport and erosion/deposition simulation "
	  "using path sampling method (SIMWE)");
                  
  if (G_get_set_window (&cellhd) == -1)
    exit (EXIT_FAILURE);

  conv = G_database_units_to_meters_factor();

  mixx = cellhd.west * conv;
  maxx = cellhd.east * conv;
  miyy = cellhd.south * conv;
  mayy = cellhd.north * conv;

  stepx = cellhd.ew_res * conv;
  stepy = cellhd.ns_res * conv;
/*  step = amin1(stepx,stepy);*/
  step = (stepx + stepy) / 2.;
  mx = cellhd.cols;
  my = cellhd.rows;
  x_orig = cellhd.west * conv;
  y_orig = cellhd.south * conv;/* do we need this? */
    xmin = 0.;
    ymin = 0.;
    xp0 = xmin + stepx / 2.;
    yp0 = ymin + stepy / 2.;
    xmax = xmin + stepx * (float) mx;
    ymax = ymin + stepy * (float) my;
  G_message(_("x,y %f, %f"), xmax, ymax);
    /*	
    bxmi=2093113. * conv;
    bymi=731331. * conv;
    bxma=2093461. * conv;
    byma=731529. * conv;
    bresx=2. * conv;
    bresy=2. * conv;
    maxwab=100000;

    hhc = hhmax = 0.;

    mx2o= (int)((bxma-bxmi)/bresx);
    my2o= (int)((byma-bymi)/bresy);
    */
/* relative small box coordinates: leave 1 grid layer for overlap */
    /*
    bxmi = bxmi - mixx + stepx;
    bymi = bymi - miyy + stepy;
    bxma = bxma - mixx - stepx;
    byma = byma - miyy - stepy;
    mx2 = mx2o - 2*((int) (stepx / bresx));
    my2 = my2o - 2*((int) (stepy / bresy)); 
  G_message(_("bresx,bresy %f, %f"), bresx, bresy);
   */

  parm.elevin = G_define_option();
  parm.elevin->key = "elevin";
  parm.elevin->type = TYPE_STRING;
  parm.elevin->required = YES;
  parm.elevin->gisprompt = "old,cell,raster";
  parm.elevin->description = _("Name of the elevation raster map [m]");
  parm.elevin->guisection  = _("Input_options");

  parm.wdepth = G_define_option();
  parm.wdepth->key = "wdepth";
  parm.wdepth->type = TYPE_STRING;
  parm.wdepth->required = YES;
  parm.wdepth->gisprompt = "old,cell,raster";
  parm.wdepth->description = _("Name of the water depth raster map [m]");
  parm.wdepth->guisection  = _("Input_options");

  parm.dxin = G_define_option();
  parm.dxin->key = "dxin";
  parm.dxin->type = TYPE_STRING;
  parm.dxin->required = YES;
  parm.dxin->gisprompt = "old,cell,raster";
  parm.dxin->description = _("Name of the x-derivatives raster map");
  parm.dxin->guisection  = _("Input_options");

  parm.dyin = G_define_option();
  parm.dyin->key = "dyin";
  parm.dyin->type = TYPE_STRING;
  parm.dyin->required = YES;
  parm.dyin->gisprompt = "old,cell,raster";
  parm.dyin->description = _("Name of the y-derivatives raster map");
  parm.dyin->guisection  = _("Input_options");

  parm.detin = G_define_option();
  parm.detin->key = "detin";
  parm.detin->type = TYPE_STRING;
  parm.detin->required = YES;
  parm.detin->gisprompt = "old,cell,raster";
  parm.detin->description =
	_("Name of the detachment capacity coefficient raster map [s/m]");
  parm.detin->guisection  = _("Input_options");

  parm.tranin = G_define_option();
  parm.tranin->key = "tranin";
  parm.tranin->type = TYPE_STRING;
  parm.tranin->required = YES;
  parm.tranin->gisprompt = "old,cell,raster";
  parm.tranin->description =
	_("Name of the transport capacity coefficient raster map [s]");
  parm.tranin->guisection  = _("Input_options");

  parm.tauin = G_define_option();
  parm.tauin->key = "tauin";
  parm.tauin->type = TYPE_STRING;
  parm.tauin->required = YES;
  parm.tauin->gisprompt = "old,cell,raster";
  parm.tauin->description = _("Name of the critical shear stress raster map [Pa]");
  parm.tauin->guisection  = _("Input_options");

  parm.manin = G_define_option();
  parm.manin->key = "manin";
  parm.manin->type = TYPE_STRING;
  parm.manin->required = YES;
  parm.manin->gisprompt = "old,cell,raster";
  parm.manin->description = _("Name of the Mannings n raster map");
  parm.manin->guisection  = _("Input_options");

/* needs to be updated to GRASS 6 vector format !! */
  parm.sfile = G_define_option ();
  parm.sfile->key = "vector";
  parm.sfile->type = TYPE_STRING;
  parm.sfile->required = NO;
  parm.sfile->gisprompt = "old,vector,vector";
  parm.sfile->description = _("Name of the sampling locations vector points map");
  parm.sfile->guisection  = _("Input_options");

  parm.tc = G_define_option();
  parm.tc->key = "tc";
  parm.tc->type = TYPE_STRING;
  parm.tc->required = NO;
  parm.tc->gisprompt = "new,cell,raster";
  parm.tc->description = _("Output transport capacity raster map [kg/ms]");
  parm.tc->guisection  = _("Output_options");

  parm.et = G_define_option();
  parm.et->key = "et";
  parm.et->type = TYPE_STRING;
  parm.et->required = NO;
  parm.et->gisprompt = "new,cell,raster";
  parm.et->description =
	_("Output transp.limited erosion-deposition raster map [kg/m2s]");
  parm.et->guisection  = _("Output_options");

  parm.conc = G_define_option();
  parm.conc->key = "conc";
  parm.conc->type = TYPE_STRING;
  parm.conc->required = NO;
  parm.conc->gisprompt = "new,cell,raster";
  parm.conc->description = _("Output sediment concentration raster map [particle/m3]");
  parm.conc->guisection  = _("Output_options");

  parm.flux = G_define_option();
  parm.flux->key = "flux";
  parm.flux->type = TYPE_STRING;
  parm.flux->required = NO;
  parm.flux->gisprompt = "new,cell,raster";
  parm.flux->description = _("Output sediment flux raster map [kg/ms]");
  parm.flux->guisection  = _("Output_options");

  parm.erdep = G_define_option();
  parm.erdep->key = "erdep";
  parm.erdep->type = TYPE_STRING;
  parm.erdep->required = NO;
  parm.erdep->gisprompt = "new,cell,raster";
  parm.erdep->description = _("Output erosion-deposition raster map [kg/m2s]");
  parm.erdep->guisection  = _("Output_options");

  parm.nwalk = G_define_option();
  parm.nwalk->key = "nwalk";
  parm.nwalk->type = TYPE_INTEGER;
  parm.nwalk->answer = NWALK;
  parm.nwalk->required = NO;
  parm.nwalk->description = _("Number of walkers");
  parm.nwalk->guisection  = _("Parameters");

  parm.niter = G_define_option();
  parm.niter->key = "niter";
  parm.niter->type = TYPE_INTEGER;
  parm.niter->answer = NITER;
  parm.niter->required = NO;
  parm.niter->description = _("Number of time iterations (minutes)");
  parm.niter->guisection  = _("Parameters");

  parm.outiter = G_define_option();
  parm.outiter->key = "outiter";
  parm.outiter->type = TYPE_INTEGER;
  parm.outiter->answer = ITEROUT;
  parm.outiter->required = NO;
  parm.outiter->description = _("Time step for saving output maps (minutes)");
  parm.outiter->guisection  = _("Parameters");

  parm.density = G_define_option();
  parm.density->key = "density";
  parm.density->type = TYPE_INTEGER;
  parm.density->answer = DENSITY;
  parm.density->required = NO;
  parm.density->description = _("Density of output walkers");
  parm.density->guisection  = _("Parameters");

  parm.diffc = G_define_option();
  parm.diffc->key = "diffc";
  parm.diffc->type = TYPE_DOUBLE;
  parm.diffc->answer = DIFFC;
  parm.diffc->required = NO;
  parm.diffc->description = _("Water diffusion constant");
  parm.diffc->guisection  = _("Parameters");

  /*
  flag.mscale = G_define_flag ();
  flag.mscale->key = 'm';
  flag.mscale->description = _("Multiscale simulation");
  */
  flag.tserie = G_define_flag ();
  flag.tserie->key = 't';
  flag.tserie->description = _("Time-series (dynamic) output");

  if (G_parser (argc, argv))
    exit (EXIT_FAILURE);

/*  mscale=flag.mscale->answer ? "m" : NULL; */
/* water has this ts=flag.tserie->answer; */
  tserie=flag.tserie->answer ? "t" : NULL;

  elevin = parm.elevin->answer;
  wdepth = parm.wdepth->answer;
  dxin = parm.dxin->answer;
  dyin = parm.dyin->answer;
/*  maskmap = parm.maskmap->answer;*/
  detin = parm.detin->answer;
  tranin = parm.tranin->answer;
  tauin = parm.tauin->answer;
  manin = parm.manin->answer;
  tc = parm.tc->answer;
  et = parm.et->answer;
  conc = parm.conc->answer;
  flux = parm.flux->answer;
  erdep = parm.erdep->answer;
  sfile = parm.sfile->answer;

  sscanf(parm.nwalk->answer, "%d", &maxwa);
  sscanf(parm.niter->answer, "%d", &timesec);
  sscanf(parm.outiter->answer, "%d", &iterout);
  sscanf(parm.density->answer, "%d", &ldemo);
  sscanf(parm.diffc->answer, "%lf", &frac);

  timesec = timesec * 60.0;
  iterout = iterout * 60.0;
  G_message(_("Converted time %f, %f"), timesec, timesec);

    rwalk = (double) maxwa;
  G_message(_("rwalk %f"), rwalk);

  if (conv != 1.0) 
    G_message(_("Using metric conversion factor %f, step=%f"), conv, step);

  /*
   * G_set_embedded_null_value_mode(1);
   */

  if ((tc == NULL) && (et == NULL) && (conc == NULL) && (flux == NULL) && (erdep == NULL))
    G_warning(_("You are not outputing any raster or site files"));

  ret_val = input_data();
  if (ret_val != 1)
    G_fatal_error(_("Input failed"));

/* mandatory for si,sigma */

   si = (double **)G_malloc (sizeof(double *)*(my));
       for(l=0;l<my;l++)
          {
            si[l]   = (double*)G_malloc (sizeof(double)*(mx));
           }
       for (j = 0; j < my; j++)
           {
            for (i = 0; i < mx; i++)
              si[j][i] = 0.;
             }

   sigma = (double **)G_malloc (sizeof(double *)*(my));
       for(l=0;l<my;l++)
          {
            sigma[l]   = (double*)G_malloc (sizeof(double)*(mx));
           }
       for (j = 0; j < my; j++)
           {
            for (i = 0; i < mx; i++)
              sigma[j][i] = 0.;
             }

/* memory allocation for output grids */

        dif = (float **)G_malloc (sizeof(float *)*(my));
           for(l=0;l<my;l++)
              {
                dif[l]   = (float*)G_malloc (sizeof(float)*(mx));
              }
           for (j = 0; j < my; j++)
              {
                for (i = 0; i < mx; i++)
                   dif[j][i] = 0.;
               }

   if (erdep != NULL || et != NULL)
      {
        er = (float **)G_malloc (sizeof(float *)*(my));
           for(l=0;l<my;l++)
              {
                er[l]   = (float*)G_malloc (sizeof(float)*(mx));
              }
           for (j = 0; j < my; j++)
              {
                for (i = 0; i < mx; i++)
                   er[j][i] = 0.;
               }
       }

/*  if (maskmap != NULL)
    bitmask = BM_create (cols, rows);

  IL_create_bitmask (&params, bitmask);
*/

  seeds(rand1,rand2);

  grad_check();

  if (et != NULL)
	erod(si);
/* treba dat output pre topoerdep */
/*  fprintf (stdout, "\n Percent complete: ");*/

  main_loop();

  if (tserie == NULL) {
  ii=output_data (0,1.);
  if (ii != 1)
    G_fatal_error (_("Cannot write raster maps"));
  }

  if (fdwalkers != NULL)
    fclose (fdwalkers);

  if (sfile != NULL)
                G_sites_close(fw);

  exit (EXIT_SUCCESS);
}