/*
 * SUMMARY:      SlopeAspect.c - Calculate slope and aspect of each pixel
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       William A Perkins
 * ORG:          Battelle Memorial Institute Pacific Northwest Laboratory
 * E-MAIL:       perk@clio.muse.pnl.gov
 * ORIG-DATE:    21-May-96
 * Last Change: Wed Feb  3 23:11:31 1999 by Bart Nijssen <nijssen@u.washington.edu>
 * DESCRIPTION:  This module contains two routines to compute "slope" and
 *               "aspect"  direction of slope): one which uses only terrain
 *               elevations and another which uses water table elevations.
 * DESCRIP-END.
 * FUNCTIONS:    valid_cell()
 *               slope_aspect()
 *               flow_fractions()
 *               ElevationSlopeAspect()
 *               HeadSlopeAspect()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: SlopeAspect.c,v 1.2 2001/03/31 00:10:36 waichler Exp $";
#endif /* lint */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include "constants.h"
#include "settings.h"
#include "data.h"
#include "slopeaspect.h"
#include "fileio.h"

/* These indices are so neighbors can be looked up quickly */
/*
int xneighbor[NDIRS] = {
#if NDIRS == 4
  0, 1, 0, -1
#elif NDIRS == 8
  -1, 0, 1, 1, 1, 0, -1, -1
#endif
};
int yneighbor[NDIRS] = {
#if NDIRS == 4
  -1, 0, 1, 0
#elif NDIRS == 8
  1, 1, 1, 0, -1, -1, -1, 0
#endif
};
*/

  

/* -------------------------------------------------------------
   valid_cell
   Checks to see if grid indices, x and y, are within the grid 
   defined by the specified Map
   ------------------------------------------------------------- */
int 
valid_cell(MAPSIZE *Map, int x, int y) 
{
  return (x >= 0 && y >= 0 && x < Map->NX && y < Map->NY);
}



/* -------------------------------------------------------------
   slope_aspect
   Calculation of slope and aspect given elevations of cell and neighbors
   ------------------------------------------------------------- */
static void slope_aspect(float dx, float dy, float celev, float
                         nelev[NDIRS], float* slope, float* aspect)
{
  int n;
  float dzdx = 0, dzdy = 0;

  *slope = 0.0;
  *aspect = 0.0;
/* test  for (n = 0; n < NDIRS; n++) nelev[n] = 100.0; */
  
  switch (NDIRS) {
  case 8:
    /* for eight neighbors, this is
       exactly the same algorithm that
       Arc/Info uses */
    
    for (n = 0; n < NDIRS; n++) {
      if (nelev[n] == OUTSIDEBASIN) {
        nelev[n] = celev;
      }
    }
    
    dzdx = ((nelev[7] + 2*nelev[6] + nelev[5]) -
            (nelev[1] + 2*nelev[2] + nelev[3]))/(8*dx);
    dzdy = ((nelev[7] + 2*nelev[0] + nelev[1]) -
            (nelev[5] + 2*nelev[4] + nelev[3]))/(8*dy);
    break;
    
    
  case 4:
    if (nelev[1] == (float)OUTSIDEBASIN && nelev[3] == (float)OUTSIDEBASIN) {
      dzdx = 0.0;
    } else if (nelev[1] == (float)OUTSIDEBASIN) {
      dzdx = (nelev[3] - celev)/dx;
    } else if (nelev[3] == (float)OUTSIDEBASIN) {
      dzdx = (celev - nelev[1])/dx;
    } else {
      dzdx = (nelev[3] - nelev[1])/(2*dx);
    }
    
    if (nelev[0] == (float)OUTSIDEBASIN && nelev[2] == (float)OUTSIDEBASIN) {
      dzdy = 0.0;
    } else if (nelev[2] == (float)OUTSIDEBASIN) {
      dzdy = (celev - nelev[0])/dy;
    } else if (nelev[0] == (float)OUTSIDEBASIN) {
      dzdy = (nelev[2] - celev)/dy;
    } else {
      dzdy = (nelev[2] - nelev[0])/(2*dy);
    }
    break;
  default:
    assert(0);                  /* nothing else works */ 
    break;
  }
  
  *slope = sqrt(dzdx*dzdx + dzdy*dzdy); if (dzdx == 0.0 && dzdy == 0.0) {
    *aspect = 0.0;
  } else {
    *aspect = atan2(dzdx, dzdy);
  }
  
  return;
}

/* -------------------------------------------------------------
   flow_fractions
   Computes subsurface flow fractions given the slope and aspect 
------------------------------------------------------------- */
static void flow_fractions(float dx, float dy, float slope, float aspect,
                           float nelev[NDIRS], float *grad, 
                           unsigned char dir[NDIRS], unsigned int *total_dir)
{
  float cosine = 0.0;
  float sine = 0.0;
  float total_width = 0.0, effective_width = 0.0;
  int n;

  cosine = cos(aspect);
  sine = sin(aspect);

  switch (NDIRS) {
  case 4:
    
    /* fudge any cells which flow outside
       the basin by just pointing the
       aspect in the opposite direction */
    
    if ( (cosine > 0 && nelev[0] == (float)OUTSIDEBASIN) ||
         (cosine < 0 && nelev[2] == (float)OUTSIDEBASIN) )
      cosine = -cosine;
    if ( (sine > 0 && nelev[1] == (float)OUTSIDEBASIN) ||
         (sine < 0 && nelev[3] == (float)OUTSIDEBASIN) )
      sine = -sine;
    
    /* compute flow widths */
    total_width = fabs(cosine)*dx + fabs(sine)*dy; 
    *grad = slope*total_width;
    *total_dir = 0;
    for (n = 0; n < NDIRS; n++) {
      switch (n) {
      case 0:
        effective_width = (cosine > 0 ? cosine*dx : 0.0); break;
      case 2:
        effective_width = (cosine < 0 ? -cosine*dx : 0.0); break;
      case 1:
        effective_width = (sine > 0 ? sine*dx : 0.0); break;
      case 3:
        effective_width = (sine < 0 ? -sine*dx : 0.0); break;
      default:
        assert(0);              /* How can this happen? */
      }
      dir[n] = (int)((effective_width/total_width)*255.0 + 0.5);
      *total_dir += dir[n];
    }
    break;
  case 8:
    assert(0);                  /* can't handle this */ break;
  default:
    assert(0);                  /* other cases don't work either */
  }
  return;
}

/* -------------------------------------------------------------
   ElevationSlopeAspect
   ------------------------------------------------------------- */
void ElevationSlopeAspect(MAPSIZE *Map, TOPOPIX **TopoMap, int ReadFlowDir,
                          char FlowDirFileName[BUFSIZE+1])
{
  int x = 0, xn = 0;
  int y = 0, yn = 0;
  int n, NumLoops = 0;
  float neighbor_elev[NDIRS];
  
  for (n = 0; n < NDIRS; n++) neighbor_elev[n] = 0.0;

  /* fill neighbor array */
  for (x = 0; x < Map->NX; x++) {
    for (y = 0; y < Map->NY; y++) {
      if (TopoMap[y][x].Mask) {
        
        for (n = 0; n < NDIRS; n++) {
          xn = x + Map->xneighbor[n];
          yn = y + Map->yneighbor[n];
          if (valid_cell(Map, xn, yn)) {
            neighbor_elev[n] =
              ((TopoMap[yn][xn].Mask) ? TopoMap[yn][xn].Dem : (float)OUTSIDEBASIN);
          } else {
            neighbor_elev[n] = OUTSIDEBASIN;
          }
        }
        slope_aspect(Map->DX, Map->DY, TopoMap[y][x].Dem, neighbor_elev,
                     &(TopoMap[y][x].Slope), &(TopoMap[y][x].Aspect));
        
        /* fill Dirs in TopoMap too */
        if( ReadFlowDir == FALSE ) {
          flow_fractions(Map->DX, Map->DY, TopoMap[y][x].Slope,
                         TopoMap[y][x].Aspect,
                         neighbor_elev, &(TopoMap[y][x].FlowGrad),
                         TopoMap[y][x].Dir, &(TopoMap[y][x].TotalDir));
        }
      }
    }
  }

  if( ReadFlowDir == TRUE )
  {
    printf("Reading Flow Fractions from File\n");
    ReadFlowDirs( *Map, TopoMap, FlowDirFileName);
  }

  /* check for flow loops and flow out of the basin */
  NumLoops = CheckForLoops( *Map, TopoMap);

  return;
}

/* -------------------------------------------------------------
   ReadFlowDirs
   Reads in flow directions flow user specified file. 
   ------------------------------------------------------------- */
void ReadFlowDirs( MAPSIZE Map, TOPOPIX **TopoMap, char FlowDirFileName[BUFSIZE+1])
{

  FILE *DirFile;
	unsigned char *Dir;
	int n, x, y, NumRead, InBasin = 0, row = 0, TotalDir = 0;
	
	Dir = Build1DUChar(NDIRS * Map.NY * Map.NX);	/* Allocate temp memory  */

  /* open flow direction file */ 
  OpenFile(&DirFile, FlowDirFileName, "r", FALSE);

	NumRead = fread( Dir, sizeof(unsigned char), NDIRS * Map.NX * Map.NY, DirFile);
	if( NumRead != NDIRS * Map.NY * Map.NX ){
		fprintf(stderr, "Error: %d Flow Dir Values Read\n", NumRead);
	}

	fclose( DirFile );

	/* Assign values */
    row = 0;
	for( y = 0; y < Map.NY; y++ ){
		for( x = 0; x < Map.NX; x++ ){

      TotalDir = 0;
      for( n = 0; n < NDIRS; n++ ){
          TopoMap[y][x].Dir[n] = Dir[row];
          TotalDir += Dir[row];
          row++;
      }
      TopoMap[y][x].TotalDir = TotalDir;

      if( TopoMap[y][x].Slope > 0 )
        TopoMap[y][x].FlowGrad = TopoMap[y][x].Slope * Map.DX;
      else
        TopoMap[y][x].FlowGrad = MIN_SLOPE * Map.DX;

      if( TopoMap[y][x].Mask && TotalDir == 0 )
      {
        printf("Error: cell %d %d has no flow directions\n", y, x);
        exit(1);
      }
      if( TotalDir > 255 )
      {
        printf("Error: TotalDir for cell %d %d > 255\n", y, x);
        exit(1);
      }
    }
  }

  if( row != NDIRS * Map.NY * Map.NX )
  {
    printf("Error reading flow directions\n");
    exit(1);
  }

  Free1DUChar( Dir );

}

/* -------------------------------------------------------------
   CheckForLoops
   Check for loops in flow directions. 
   ------------------------------------------------------------- */
int CheckForLoops( MAPSIZE Map, TOPOPIX **TopoMap)
{
  int x, y, n, nn, NumLoops = 0, NumFlowOutOfBasin = 0;
    
    /* fill neighbor array */
    for(x = 0; x < Map.NX; x++)
    {
      for(y = 0; y < Map.NY; y++)
      {
        if(TopoMap[y][x].Mask)
        {
           for(n = 0; n < NDIRS; n++)
           {
              int     xn = x + Map.xneighbor[n];
              int     yn = y + Map.yneighbor[n];
                    
              if(xn >= 0 && xn < Map.NX && yn >= 0 && yn < Map.NY && TopoMap[yn][xn].Mask)
              {
                if( n < (NDIRS / 2) )
                  nn = n + (NDIRS / 2);
                else
                  nn = n - (NDIRS / 2);

                if( TopoMap[y][x].Dir[n] > 0 && TopoMap[yn][xn].Dir[nn] > 0 )
                {
                  NumLoops++;
                    printf("Loop between cell %d %d and %d %d in the %d direction\n",
                      y, x, yn, xn, n);
                }

              }
              else
              {
                if( TopoMap[y][x].Dir[n] > 0 )
                {
                  printf("Cell %d %d has flow outside the basin in the %d direction\n",
                  y, x, n);
                  NumFlowOutOfBasin++;
                }
                  
              }
           }
              
        }
      }
    }

  printf("Number of flow loops: %d\n", NumLoops);
  printf("Cells with flow out of the basin: %d\n", NumFlowOutOfBasin);
    /* exit(1); */

  return(NumLoops);
}

/* -------------------------------------------------------------
   HeadSlopeAspect
   This computes slope and aspect using the water table elevation. 
   ------------------------------------------------------------- */
void HeadSlopeAspect(MAPSIZE *Map, TOPOPIX **TopoMap, SOILPIX **SoilMap)
{

  int x;
  int y;
  int n;
  float neighbor_elev[NDIRS];

  /* let's assume for now that WaterLevel is the SOILPIX map is
     computed elsewhere */
  for (x = 0; x < Map->NX; x++) {
    for (y = 0; y < Map->NY; y++) {
      if (TopoMap[y][x].Mask) {
        
        float slope, aspect;
        
        for (n = 0; n < NDIRS; n++) {
          int xn = x + Map->xneighbor[n];
          int yn = y + Map->yneighbor[n];
          
          if (valid_cell(Map, xn, yn)) {
            neighbor_elev[n] =
              ((TopoMap[yn][xn].Mask) ? SoilMap[yn][xn].WaterLevel :
               OUTSIDEBASIN);
          } else {
            neighbor_elev[n] = OUTSIDEBASIN;
          }
        }
        
        slope_aspect(Map->DX, Map->DY, SoilMap[y][x].WaterLevel, neighbor_elev,
                     &slope, &aspect);
        flow_fractions(Map->DX, Map->DY, slope, aspect, neighbor_elev,
                       &(TopoMap[y][x].FlowGrad), TopoMap[y][x].Dir, 
                       &(TopoMap[y][x].TotalDir));
      }
    }
  }
  return;
}


