/* #pragma warning( disable : 4244 ) */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <malloc.h>
#include "settings.h"
#include "data.h"
#include "DHSVMerror.h"
#include "functions.h"
#include "fileio.h"
#include "getinit.h"
#include "elevationband.h"

enum KEY {input_file, output_file};

/****  MAIN FUNCTIONS   **************/

/*****************************************************************************
  InitElevBand()
*****************************************************************************/
void
InitElevBand (LISTPTR Input, MAPSIZE Map, BANDSTRUCT ** ElevBand, BANDINFO * ElevBandInfo)
{
  const char *Routine = "InitElevBand()";
  FILE *ElevBandFile;
  char data_line[120];

  int i, x, y, Huc;
  STRINIENTRY StrEnv[] = {
    {"ELEVBAND", "ELEVBAND INPUT FILE"    , ""  , NULL},
    {"ELEVBAND", "ELEVBAND OUTPUT FILE", ""  , NULL},
    {NULL     , NULL                  , ""  , NULL}
  };

  printf("Initializing elevation bands\n");
  ElevBandInfo->ReadBinFiles = FALSE;
  
  /* Read the key-entry pairs from the [ELEVBAND] section of the input file */
  for (i = 0; StrEnv[i].SectionName; i++) {
    GetInitString(StrEnv[i].SectionName, StrEnv[i].KeyName, StrEnv[i].Default,
                  StrEnv[i].VarStr, (unsigned long) BUFSIZE, Input);
    if (!StrEnv[i].VarStr)
      ReportError(StrEnv[i].KeyName, 51);
  }


  /* open elevation band file and read number of bands */
  OpenFile(&ElevBandFile, StrEnv[input_file].VarStr, "r", FALSE);

  if( fgets( data_line, 120, ElevBandFile ) == NULL ){
      printf("Error reading ElevBand File: exit program\n");
      exit(1);
  }
  else
  {
    if( sscanf(data_line, "%d %d",  &(ElevBandInfo->NumBands)  
          ,&(ElevBandInfo->PrintOutflows) ) != 2 ) {
       printf("Incorrect number of parameters on first line of ElevBand File\n");
       printf("Expecting 2 parametrs: Num of Bands and PrintOutflows (1 or 0)\n");
       printf("Exiting program\n");
       exit(1);
     }
  }
  

  if (!(*ElevBand = (BANDSTRUCT *) calloc (ElevBandInfo->NumBands, sizeof (BANDSTRUCT))))
    fprintf (stderr, "Could not allocate memory for ElevBand\n");

  ElevBandInfo->TotalNumCells = 0;

  /* read band data */
  for (i = 0; i < ElevBandInfo->NumBands; i++) {
    if( fgets( data_line, 120, ElevBandFile ) == NULL ){
        printf("Error reading ElevBand File: exit program\n");
        exit(1);
    }
    else {
      if( sscanf(data_line, "%d %g %g %g %lg %g %g %d %d %f\n", &Huc, &(*ElevBand)[i].MeanElev,
	      &(*ElevBand)[i].MeanSlope, &(*ElevBand)[i].MeanAspect,
	      &(*ElevBand)[i].Area, &(*ElevBand)[i].MeanLatitude,
	      &(*ElevBand)[i].MeanLongitude,  &(*ElevBand)[i].VegType,
        &(*ElevBand)[i].SoilType, &(*ElevBand)[i].SoilDepth) != 10 ) {
          ElevBandInfo->ReadBinFiles = TRUE;
      }

      (*ElevBand)[i].NumCells = (int)( ( (*ElevBand)[i].Area / (Map.DX * Map.DY) ) + 0.5 );
      
      ElevBandInfo->TotalNumCells += (*ElevBand)[i].NumCells;
    }
  }

  fclose (ElevBandFile);

  fprintf (stderr, "Total Number of Cells for Elevation Bands: %d\n"
           ,ElevBandInfo->TotalNumCells);


  /* if flows from each band need to be printed, set pointer to output file */
  if( ElevBandInfo->PrintOutflows == 1)
    OpenFile(&(ElevBandInfo->ElevBandOut), StrEnv[output_file].VarStr, "w", TRUE);

}

/*****************************************************************************
  ElevBandUpdate()
*****************************************************************************/
void
ElevBandUpdate(MAPSIZE Map, TOPOPIX *** TopoMap, VEGPIX ***VegMap,
              BANDSTRUCT ** ElevBand, BANDINFO * ElevBandInfo, SOILPIX ***SoilMap)
{
  /* assign values to TopoMap and VegMap*/
  int x, y, n, i = 0;
  for (x = 0; x < Map.NX; x++) {
    for (y = 0; y < Map.NY; y++) {
      if (i > ElevBandInfo->NumBands) {
	      printf ("ERROR IN ELEVATION BANDS: NumBands %d Map.NY %d\n", ElevBandInfo->NumBands, Map.NY);
        exit(1);
      }
      else {
        (*TopoMap)[y][x].Dem = (*ElevBand)[i].MeanElev;
        (*TopoMap)[y][x].Slope = (*ElevBand)[i].MeanSlope;
        (*TopoMap)[y][x].Aspect = (*ElevBand)[i].MeanAspect;
        (*TopoMap)[y][x].Area =  (*ElevBand)[i].Area;

        if( ElevBandInfo->ReadBinFiles == FALSE ) {
          (*TopoMap)[y][x].FlowGrad = (*TopoMap)[y][x].Slope 
                                      * sqrt((*TopoMap)[y][x].Area);
          (*TopoMap)[y][x].TotalDir = 255;
          for (n = 0; n < NDIRS; n++)(*TopoMap)[y][x].Dir[n] = 0;
          (*TopoMap)[y][x].Dir[0] = 255;
          (*VegMap)[y][x].Veg = (*ElevBand)[i].VegType;
          (*SoilMap)[y][x].Soil = (*ElevBand)[i].SoilType;
          (*SoilMap)[y][x].Depth = (*ElevBand)[i].SoilDepth;
        }
      i++;
     }
    }
  }

}

/*****************************************************************************
  ElevBandBuildGrids()
*****************************************************************************/
void
ElevBandBuildGrids(MAPSIZE Map, TOPOPIX *** TopoMap, VEGPIX ***VegMap,
              BANDSTRUCT ** ElevBand, BANDINFO * ElevBandInfo, SOILPIX ***SoilMap, 
              LAYER *Soil)
{
  const char *Routine = "ElevBandBuildGrids";
  int x, y, n, Type, i = 0;

  /* allocate memory for TopoMap and set initial values  -------------------------*/
  if (!(*TopoMap = (TOPOPIX **) calloc(Map.NY, sizeof(TOPOPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map.NY; y++) {
    if (!((*TopoMap)[y] = (TOPOPIX *) calloc(Map.NX, sizeof(TOPOPIX))))
      ReportError((char *) Routine, 1);
  }

  for (y = 0, i = 0; y < Map.NY; y++) {
    for (x = 0; x < Map.NX; x++, i++) {
      (*TopoMap)[y][x].Dem = (*ElevBand)[i].MeanElev;
      (*TopoMap)[y][x].Mask = 1;
      (*TopoMap)[y][x].Slope = (*ElevBand)[i].MeanSlope;
      (*TopoMap)[y][x].Aspect = (*ElevBand)[i].MeanAspect;
      (*TopoMap)[y][x].Area =  (*ElevBand)[i].Area;
      (*TopoMap)[y][x].Travel = 0;
      (*TopoMap)[y][x].FlowGrad = (*TopoMap)[y][x].Slope * sqrt((*TopoMap)[y][x].Area);
      (*TopoMap)[y][x].TotalDir = 255;
      for (n = 0; n < NDIRS; n++)(*TopoMap)[y][x].Dir[n] = 0;
      (*TopoMap)[y][x].Dir[0] = 255;
    }
  }

  /* allocate memory for SoilMap and set initial values ------------------------*/
  if (!(*SoilMap = (SOILPIX **) calloc(Map.NY, sizeof(SOILPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map.NY; y++) {
    if (!((*SoilMap)[y] = (SOILPIX *) calloc(Map.NX, sizeof(SOILPIX))))
      ReportError((char *) Routine, 1);
  }
  
  for (y = 0, i = 0; y < Map.NY; y++) {
    for (x = 0; x < Map.NX; x++, i++) {

      Type = (*ElevBand)[i].SoilType;

      if (Type > Soil->NTypes) 
	        ReportError((char *) Routine, 32);

      (*SoilMap)[y][x].Soil = Type;
      (*SoilMap)[y][x].Depth = (*ElevBand)[i].SoilDepth;

      /* allocate memory for the number of root layers, plus an additional 
	       layer below the deepest root layer */
      if (!((*SoilMap)[y][x].Moist = 
              (float *) calloc((Soil->NLayers[Type-1]+1), sizeof(float))))
          ReportError((char *) Routine, 1);
      if (!((*SoilMap)[y][x].Perc = 
              (float *) calloc(Soil->NLayers[Type-1], sizeof(float))))
          ReportError((char *) Routine, 1);
      if (!((*SoilMap)[y][x].Temp = 
              (float *) calloc(Soil->NLayers[Type-1], sizeof(float))))
         ReportError((char *) Routine, 1);
    }
  }

  /* allocate memory for VegMap and set initial values ------------------------*/
  if (!(*VegMap = (VEGPIX **) calloc(Map.NY, sizeof(VEGPIX *))))
    ReportError((char *) Routine, 1);
  for (y = 0; y < Map.NY; y++) {
    if (!((*VegMap)[y] = (VEGPIX *) calloc(Map.NX, sizeof(VEGPIX))))
      ReportError((char *) Routine, 1);
  }
  for (y = 0, i = 0; y < Map.NY; y++) {
    for (x = 0; x < Map.NX; x++, i++) {
      (*VegMap)[y][x].Veg = (*ElevBand)[i].VegType;
      (*VegMap)[y][x].Tcanopy = 0.0;
    }
  }

}

/*****************************************************************************
  InitElevBandDump()
*****************************************************************************/
void
InitElevBandDump (BANDINFO * ElevBandInfo, char *DumpPath)
{

  char buffer[NAMESIZE];

  /* Open elevation band file for future use */
  /*MakeFileName (DumpPath, "ElevBand.Out", "", buffer);*/
  OpenFile (&(ElevBandInfo->ElevBandOut), buffer, "w", TRUE);

}



/*****************************************************************************
  PrintElevBand()
*****************************************************************************/
void
PrintElevBand (DATE *Current, BANDINFO ElevBandInfo, Channel *net, BANDSTRUCT *ElevBand)
{
  int BandCount = 0;
  int order_count;
  static int first_time = 1;
  unsigned int order;
  Channel *current;

  if( first_time ){
    first_time = 0;
    fprintf (ElevBandInfo.ElevBandOut, "Order of Segments: ");
    for (order = 1; ; order += 1) {
      order_count = 0;
      current = net;
      while (current != NULL) {
        if (current->order == order) {
          fprintf(ElevBandInfo.ElevBandOut, " %d ", current->id );
          order_count += 1;
        }
        current = current->next;
      }
      if (order_count == 0) break;
    }
    fprintf (ElevBandInfo.ElevBandOut, "\n");
  }
  
  
  /* quick_flow */
  PrintDate(Current, ElevBandInfo.ElevBandOut);
  for (order = 1; ; order += 1) {
    order_count = 0;
    current = net;
    while (current != NULL) {
      if (current->order == order) {
        fprintf(ElevBandInfo.ElevBandOut, " %g ", current->quick_flow );
        order_count += 1;
      }
      current = current->next;
    }
    if (order_count == 0) break;
  }
  fprintf (ElevBandInfo.ElevBandOut, "\n");

  /* matrix_flow */
  PrintDate(Current, ElevBandInfo.ElevBandOut);
  for (order = 1; ; order += 1) {
    order_count = 0;
    current = net;
    while (current != NULL) {
      if (current->order == order) {
        fprintf(ElevBandInfo.ElevBandOut, " %g ", 
               (current->lateral_inflow - current->quick_flow) );
        order_count += 1;
      }
      current = current->next;
    }
    if (order_count == 0) break;
  }
  fprintf (ElevBandInfo.ElevBandOut, "\n");
}


