//  -*- Mode: C++;  -*-
//  File: grid.cc
//  Author: Dino Bellugi (dino@geomorph.berkeley.edu)
//  Copyright  Dino Bellugi, BlueG SoftWear, U.C. Berkeley, 1999 (C)
// *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// * FUNCTION: implementation of all grid related functions
// *
// * RELATED PACKAGES: all
// *
// * HISTORY:
// * Created:  Mon Nov  6 1995 (dino, based on rob reiss' version)
// * Modified: Sat Jun 26 1999 (dino)
// *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

// the following include file (notice.h) contains important copyright 
// information, as well as other legal statements. it must be present 
// in all files and in any distribution. removal or modification is
// strictly forbidden and violates u.s. and international copyright law.
#include "notice.h"
//
#include <cstdlib>
#include <fstream>
#include <ctime>
#include <cmath>
#include "macros.h"
#include "grid.h"
#include "display.h"
#include <iostream>
using namespace std;

#ifdef COUNT_SIDES
extern double lcount, rcount, ucount, dcount;
#endif

Grid::Grid(int x, int y, double xo, double yo, double spac, 
	double nodat, int flds) 
{
	// Creation of Grid
	xdim = x;
  	ydim = y;
  	xorg = xo;
  	yorg = yo;
  	xcur = ycur = xprev1 = yprev1 = xprev2 = yprev2 = 0;
  	spacing = xsize = ysize = spac;
  	nodata = nodat;
	fields = flds;

	for (int k = 0; k < fields; ++k)
	  	arrValues[k] = new double[x * y];

	timeseed();
}

Grid::Grid(Grid &z) 
{
	int i, j, k;

  	// Creation of Grid
  	xdim = z.xdim;
  	ydim = z.ydim;
  	xorg = z.xorg;
  	yorg = z.yorg;
  	xcur = ycur = xprev1 = yprev1 = xprev2 = yprev2 = 0;
  	spacing = z.spacing;
  	xsize = z.xsize;
  	ysize = z.ysize;
  	nodata = z.nodata;
	fields = z.fields;

 	for (k = 0; k < fields; ++k)
 		arrValues[k] = new double[z.xdim * z.ydim];

	for (k = 0; k < fields; ++k)
	{
	  	for (j = 0; j < ydim; j++) 
	  	{
	    	for(i = 0; i < xdim; i++) 
			{
	      		at(i, j, k) = z.at(i, j, k);
	    	}
	  	}
	}
	timeseed();
}

Grid::~Grid() 
{
	for (int k = 0; k < fields; ++k)
  		delete[] arrValues[k];
}

int Grid::addField(int srcfield)	// adds a field at the end of arrValues
{
	if (fields == MAXFIELDS ||
		(arrValues[fields] = new double[xdim * ydim]) == (double *)NULL)
	{
		cerr << "Not enough memory for " << fields + 1 << " fields ...\n";
		return(0);
	}
	++fields;
	if (srcfield >= 0)
		setValue(fields - 1, srcfield);
	return(1);
}

int Grid::delField()	// deletes last field from arrValues
{
	if (fields < 1 || arrValues[fields - 1] == (double *)NULL)
		return(0);

  	delete[] arrValues[--fields];
	return(1);
}

int Grid::setSize(int rows, int cols, int init, int flds)
{
	int k;

	for (int k = 0; k < fields; ++k)
  		delete[] arrValues[k];

	xdim = cols;
	ydim = rows;
	fields = flds;

 	for (k = 0; k < fields; ++k)
 		if ((arrValues[k] = new double[xdim * ydim]) == (double *)NULL)
			return(0);

	if (init)
	{
  		xorg = yorg = xcur = ycur = 
		xprev1 = yprev1 = xprev2 = yprev2 = 0;
		spacing = xsize = ysize = DEF_SIZE;
	}
	return(1);
}

int Grid::setSize(Grid &z)
{
	int k;

	for (int k = 0; k < fields; ++k)
  		delete[] arrValues[k];

	xdim = z.xdim;
	ydim = z.ydim;
  	xorg = z.xorg;
  	yorg = z.yorg;
  	xcur = ycur = xprev1 = yprev1 = xprev2 = yprev2 = 0;
  	spacing = z.spacing;
  	xsize = z.xsize;
  	ysize = z.ysize;
	nodata = z.nodata;
	fields = z.fields;

 	for (k = 0; k < fields; ++k)
 		if ((arrValues[k] = new double[xdim * ydim]) == (double *)NULL)
			return(0);

	return(1);
}


void Grid::setNodata(double value)	{nodata = value;}

void Grid::chgNodata(double value, int set, int field)	
{
	chgValue(nodata, value, field);
	if (set)
		setNodata(value);
}

void Grid::setValue(double value, int field)
{
        int i, j;
 
        for (j = 0; j < ydim; j++)
                for(i = 0; i < xdim; i++)
                        at(i, j, field) = value;
}
 
void Grid::setValue(int field1, int field2)
{
	int i, j;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++) 
			at(i, j, field1) = at(i, j, field2);
}

void Grid::setValue(int field1, int field2, double value)
{
	int i, j;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++)
			if (isValue(i, j, field2))
				at(i, j, field1) = at(i, j, field2) + value;
			else
				at(i, j, field1) = nodata;
}

void Grid::setValue(int field1, double value, int field2)
{
	int i, j;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++) 
			if (isValue(i, j, field2))
				at(i, j, field1) = value;
			else
				at(i, j, field1) = nodata;
}

long Grid::countValue(double value, int field)
{
	int i, j;
	long cnt = 0L;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++) 
			if (at(i, j, field) == value)
				++cnt;
	return(cnt);
}

long Grid::countValue(int field)
{
	int i, j;
	long cnt = 0L;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++) 
			if (isValue(i, j, field))
				++cnt;
	return(cnt);
}

void Grid::chgValue(double old_value, double new_value, int field)
{
	int i, j;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++) 
			if (at(i, j, field) == old_value)
				at(i, j, field) = new_value;
}

void Grid::chgValue(int nvalues, double *old_value, double new_value, int field)
{
	int i, j, n;

	for (j = 0; j < ydim; ++j) 
	{
		for(i = 0; i < xdim; ++i) 
		{
			for (n = 0; n < nvalues; ++n)
			{
				if (at(i, j, field) == old_value[n])
				{
					at(i, j, field) = new_value;
					break;
				}
			}
		}
	}
}

void Grid::addValue(int field1, int field2)
{
	int i, j;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++)
			if (isValue(i, j, field1) && isValue(i, j, field2))
				at(i, j, field1) += at(i, j, field2);
}

void Grid::subValue(int field1, int field2)
{
	int i, j;

	for (j = 0; j < ydim; j++) 
		for(i = 0; i < xdim; i++) 
			if (isValue(i, j, field1) && isValue(i, j, field2))
				at(i, j, field1) -= at(i, j, field2);
}

void Grid::mkBinary(double threshold, int reverse, int infield, int outfield)
{
        int i, j;
 
        for (j = 0; j < ydim; j++)
	{
                for(i = 0; i < xdim; i++)
		{
			if (isValue(i, j, infield))
			{
                        	if (at(i, j, infield) > threshold)
					at(i, j, outfield) = (reverse ? 0 : 1);
				else
					at(i, j, outfield) = (reverse ? 1 : 0);
			}
			else
				at(i, j, outfield) = at(i, j, infield);
		}
	}
}

double Grid::getQoverT(double theta, double aoverb, double ps, double phi)
{
	if (theta == nodata || aoverb == nodata || ps <= 0.0 || phi <= 0.0)
		return(nodata);

	double ps_pw = ps / 1000.0, pw_ps = 1000.0 / ps, q_t;

	phi = deg2rad(phi);
 
	if (tan(theta) >= tan(phi))
	{
		q_t = DEF_STEEP;
	}
	else if (tan(theta) < ((1.0 - pw_ps) * tan(phi)))
	{
		q_t = DEF_FLAT;
	}
	else
	{
		q_t = log10((sin(theta) / aoverb) * ps_pw *
			(1.0 - (tan(theta) / tan(phi))));
	}

	return(q_t);
}

double Grid::getQoverT(int i, int j, int theta_field, int aoverb_field,
	double ps, double phi)
{
	double theta, aoverb;

        if ((theta = at(i, j, theta_field)) == nodata ||
		(aoverb = at(i, j, aoverb_field)) == nodata ||
		ps <= 0.0 || phi <= 0.0)
                return(nodata);
 
        double ps_pw = ps / 1000.0, pw_ps = 1000.0 / ps, q_t;
 
        phi = deg2rad(phi);
 
        if (tan(theta) >= tan(phi))
        {
                q_t = DEF_STEEP;
        }
        else if (tan(theta) < ((1.0 - pw_ps) * tan(phi)))
        {
                q_t = DEF_FLAT;
        }
        else
        {
                q_t = log10((sin(theta) / aoverb) * ps_pw *
                        (1.0 - (tan(theta) / tan(phi))));
        }
 
        return(q_t);
}
 
double Grid::getQoverT(int i, int j, int theta_field, double aoverb,
        double ps, double phi)
{
        double theta;
 
        if ((theta = at(i, j, theta_field)) == nodata ||
                aoverb == nodata || ps <= 0.0 || phi <= 0.0)
                return(nodata);
 
        double ps_pw = ps / 1000.0, pw_ps = 1000.0 / ps, q_t;
 
        phi = deg2rad(phi);
 
        if (tan(theta) >= tan(phi))
        {
                q_t = DEF_STEEP;
        }
        else if (tan(theta) < ((1.0 - pw_ps) * tan(phi)))
        {
                q_t = DEF_FLAT;
        }
        else
        {
                q_t = log10((sin(theta) / aoverb) * ps_pw *
                        (1.0 - (tan(theta) / tan(phi))));
        }
 
        return(q_t);
}
 
double Grid::getQoverT(int i, int j, double theta, int aoverb_field,
        double ps, double phi)
{
        double aoverb;
 
        if ((aoverb = at(i, j, aoverb_field)) == nodata ||
                theta == nodata || ps <= 0.0 || phi <= 0.0)
                return(nodata);
 
        double ps_pw = ps / 1000.0, pw_ps = 1000.0 / ps, q_t;
 
        phi = deg2rad(phi);
 
        if (tan(theta) >= tan(phi))
        {
                q_t = DEF_STEEP;
        }
        else if (tan(theta) < ((1.0 - pw_ps) * tan(phi)))
        {
                q_t = DEF_FLAT;
        }
        else
        {
                q_t = log10((sin(theta) / aoverb) * ps_pw *
                        (1.0 - (tan(theta) / tan(phi))));
        }
 
        return(q_t);
}
 
double Grid::getQoverT2(int i, int j, int theta_field, int aoverb_field,
	double ps, double phi) // returns q/T (not log) with no checks
{
	double theta, aoverb;

        if ((theta = at(i, j, theta_field)) == nodata ||
		(aoverb = at(i, j, aoverb_field)) == nodata ||
		ps <= 0.0 || phi <= 0.0)
                return(nodata);
 
        double ps_pw = ps / 1000.0, pw_ps = 1000.0 / ps;
 
        phi = deg2rad(phi);
 
        return((sin(theta) / aoverb) * ps_pw * (1.0 - (tan(theta) / tan(phi))));
 }
 
double Grid::get3dFS(int i, int j, int slope_field, int aspect_field, int soil_field, 
	int cb_field, int cl_field, double root_depth, int sat_field, double ps, double phi)
{
	double theta, cl, cb, z, a, b, c, d, z1, z2, z3, z4, z12, z22, z32, z42, h, ka, kp, ko, 
		frc, fs, basal, llateral, rlateral, upslope, dnslope, shalstab, aspect;
	const double g = 9.81, pw = 1000.0;

	// check data
        if ((theta = at(i, j, slope_field)) == nodata ||
		(cl = at(i, j, cl_field)) == nodata ||
		(cb = at(i, j, cb_field)) == nodata ||
		(h = at(i, j, sat_field)) == nodata ||
		(z = at(i, j, soil_field)) == nodata ||
		(aspect = at(i, j, aspect_field)) == nodata ||
		(a = at(i, j + 1, soil_field)) == nodata ||
		(b = at(i + 1, j, soil_field)) == nodata ||
		(c = at(i, j - 1, soil_field)) == nodata ||
		(d = at(i - 1, j, soil_field)) == nodata)
                return(nodata);
	else if (z <= 0)
		return(DEF_FLAT);

	// z1 = upslope, z3 = downslope, z2 = left side (facing downslope), z4 = right side
	if ((aspect > 315 && aspect <= 360) || (aspect >= 0 && aspect <= 45))
	{
		z1 = c; z2 = d; z3 = a; z4 = b;
	}
	else if (aspect > 45 && aspect <= 135)
	{
		z1 = d; z2 = a; z3 = b; z4 = c;
	}
	else if (aspect > 135 && aspect <= 235)
	{
		z1 = a; z2 = b; z3 = c; z4 = d;
	}
	else if (aspect > 235 && aspect <= 315)
	{
		z1 = b; z2 = c; z3 = d; z4 = a;
	}
	else
	{
		cerr << "invalid aspect value: " << aspect << "\n";
		return(nodata);
	}	

	// check ranges
	if (rad2deg(theta) <= 0 || phi >= 90)
		return(DEF_FLAT);
	if (rad2deg(theta) >= 90 || phi <= 0)
		return(DEF_STEEP);
 
	// convert angles to radians
        phi = deg2rad(phi);

	// earth pressure coefficients
	ka = sqr(tan(deg2rad(45) - (phi / 2.0)));
	kp = sqr(tan(deg2rad(45) + (phi / 2.0)));
	ko = 1 - sin(phi);

	// limit soil columns to rooting depth
	z12 = min(z1, root_depth);
	z22 = min(z2, root_depth);
	z32 = min(z3, root_depth);
	z42 = min(z4, root_depth);

#ifdef NEW_3DSTAB
	basal = (z <= root_depth ? (cb * spacing * spacing * sec(theta)) : 0);
	
	if (phi > theta)
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) + 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);
	else
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) - 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);

	llateral = (cl * spacing * z22) + 
		(0.5 * g * ko * (ps * z2 * z2 - pw * h * h) * spacing * cos(phi) * tan(phi));

	rlateral = (cl * spacing * z42) + 
		(0.5 * g * ko * (ps * z4 * z4 - pw * h * h) * spacing * cos(phi) * tan(phi));

	upslope = cl * spacing * z12; 

	dnslope = (cl * spacing * z32) + 
		(0.5 * kp * g * (ps * z3 * z3 - pw * h * h) * spacing * cos(phi - theta));

	frc = (ps * g * z * sin(theta) * spacing * spacing) +
		(0.5 * ka * g * (ps * z * z - pw * h * h) * spacing * cos(phi - theta));
#else	
	basal = (z <= root_depth ? (cb * spacing * spacing / cos(theta)) : 0);
	// shalstab = ((ps * g * z * spacing * spacing * cos(theta) * cos(theta) / cos(theta)) - 
	//	(pw * g * h * spacing * spacing * cos(theta) * cos(theta) / cos(theta))) * tan(phi);
	shalstab = ((ps * g * z * spacing * spacing * cos(theta)) - 
		(pw * g * h * spacing * spacing * cos(theta))) * tan(phi);
	// llateral = (cl * spacing * z22 * cos(theta) / cos(theta)) + 
		// (0.5 * g * z2 * cos(theta) * z2 * cos(theta) * spacing / cos(theta) * 
		// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	llateral = (cl * spacing * z22) + 
		(0.5 * g * z2 * cos(theta) * z2 * spacing  * 
		((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	// rlateral = (cl * spacing * z42 * cos(theta) / cos(theta)) + 
		// (0.5 * g * z4 * cos(theta) * z4 * cos(theta) * spacing / cos(theta) * 
		// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	rlateral = (cl * spacing * z42) + 
		(0.5 * g * z4 * cos(theta) * z4 * spacing  * 
		((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	upslope = (cl * spacing * z12 * cos(theta)) + 
		(0.5 * g * z1 * cos(theta) * z1 * cos(theta) * spacing * 
		((ps - pw) * ka + pw) * 
		cos(phi) * tan(phi));
	dnslope = (cl * spacing * z32 * cos(theta)) + 
		(0.5 * g * z3 * cos(theta) * z3 * cos(theta) * spacing * 
		((ps - pw) * kp + pw) * 
		cos(phi) * tan(phi));

	// frc = ps * g * z * sin(theta) * cos(theta) * spacing * spacing / cos(theta);
	frc = ps * g * z * sin(theta) * spacing * spacing;
#endif

	fs = (basal + shalstab + llateral + rlateral + upslope + dnslope) / frc;

#ifdef PRINT_FS
	cerr << "\nka\tkp\tko\tbasal\tshalstab\tllateral\trlateral\tupslope\tdnslope\tsum\tforce\tfs\n" <<
		ka << "\t" << kp << "\t" << ko << "\t" << basal << "\t" << shalstab << "\t" << llateral <<
		"\t" << rlateral << "\t" << upslope << "\t" << dnslope << "\t" << 
		basal + shalstab + llateral + rlateral + upslope + dnslope << "\t" << 
		frc << "\t" << fs << "\n";
#endif
 
        return(fs);
}

double Grid::get3dFS(int i, int j, int slope_field, int aspect_field, int soil_field, 
	int cb_field, int cl_field, double root_depth, double sat_ratio, double ps, double phi)
{
	double theta, cl, cb, z, a, b, c, d, z1, z2, z3, z4, z12, z22, z32, z42, h, ka, kp, ko,
		frc, fs, basal, llateral, rlateral, upslope, dnslope, shalstab, aspect;
	const double g = 9.81, pw = 1000.0;

	// check data
        if ((theta = at(i, j, slope_field)) == nodata ||
		(cl = at(i, j, cl_field)) == nodata ||
		(cb = at(i, j, cb_field)) == nodata ||
		(z = at(i, j, soil_field)) == nodata ||
		(h = z * sat_ratio) < 0 ||
		(aspect = at(i, j, aspect_field)) == nodata ||
		(a = at(i, j + 1, soil_field)) == nodata ||
		(b = at(i + 1, j, soil_field)) == nodata ||
		(c = at(i, j - 1, soil_field)) == nodata ||
		(d = at(i - 1, j, soil_field)) == nodata)
                return(nodata);
	else if (z <= 0)
		return(DEF_FLAT);

	// z1 = upslope, z3 = downslope, z2 = left side (facing downslope), z4 = right side
	if ((aspect > 315 && aspect <= 360) || (aspect >= 0 && aspect <= 45))
	{
		z1 = c; z2 = d; z3 = a; z4 = b;
	}
	else if (aspect > 45 && aspect <= 135)
	{
		z1 = d; z2 = a; z3 = b; z4 = c;
	}
	else if (aspect > 135 && aspect <= 235)
	{
		z1 = a; z2 = b; z3 = c; z4 = d;
	}
	else if (aspect > 235 && aspect <= 315)
	{
		z1 = b; z2 = c; z3 = d; z4 = a;
	}
	else
	{
		cerr << "invalid aspect value: " << aspect << "\n";
		return(nodata);
	}	

	// check ranges
	if (rad2deg(theta) <= 0 || phi >= 90)
		return(DEF_FLAT);
	if (rad2deg(theta) >= 90 || phi <= 0)
		return(DEF_STEEP);
 
	// convert angles to radians
        phi = deg2rad(phi);

	// earth pressure coefficients
	ka = sqr(tan(deg2rad(45) - (phi / 2.0)));
	kp = sqr(tan(deg2rad(45) + (phi / 2.0)));
	ko = 1 - sin(phi);

	// limit neighboring soil columns to rooting depth
	z12 = min(z1, root_depth);
	z22 = min(z2, root_depth);
	z32 = min(z3, root_depth);
	z42 = min(z4, root_depth);

#ifdef NEW_3DSTAB
	basal = (z <= root_depth ? (cb * spacing * spacing * sec(theta)) : 0);
	
	if (phi > theta)
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) + 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);
	else
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) - 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);

	llateral = (cl * spacing * z22) + 
		(0.5 * g * ko * (ps * z2 * z2 - pw * h * h) * spacing * cos(phi) * tan(phi));

	rlateral = (cl * spacing * z42) + 
		(0.5 * g * ko * (ps * z4 * z4 - pw * h * h) * spacing * cos(phi) * tan(phi));

	upslope = cl * spacing * z12; 

	dnslope = (cl * spacing * z32) + 
		(0.5 * kp * g * (ps * z3 * z3 - pw * h * h) * spacing * cos(phi - theta));

	frc = (ps * g * z * sin(theta) * spacing * spacing) +
		(0.5 * ka * g * (ps * z * z - pw * h * h) * spacing * cos(phi - theta));
#else	
	basal = (z <= root_depth ? (cb * spacing * spacing / cos(theta)) : 0);
	// shalstab = ((ps * g * z * spacing * spacing * cos(theta) * cos(theta) / cos(theta)) - 
	//	(pw * g * h * spacing * spacing * cos(theta) * cos(theta) / cos(theta))) * tan(phi);
	shalstab = ((ps * g * z * spacing * spacing * cos(theta)) - 
		(pw * g * h * spacing * spacing * cos(theta))) * tan(phi);
	// llateral = (cl * spacing * z22 * cos(theta) / cos(theta)) + 
		// (0.5 * g * z2 * cos(theta) * z2 * cos(theta) * spacing / cos(theta) * 
		// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	llateral = (cl * spacing * z22) + 
		(0.5 * g * z2 * cos(theta) * z2 * spacing  * 
		((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	// rlateral = (cl * spacing * z42 * cos(theta) / cos(theta)) + 
		// (0.5 * g * z4 * cos(theta) * z4 * cos(theta) * spacing / cos(theta) * 
		// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	rlateral = (cl * spacing * z42) + 
		(0.5 * g * z4 * cos(theta) * z4 * spacing  * 
		((ps - pw) * ko + pw) * cos(phi) * tan(phi));
	upslope = (cl * spacing * z12 * cos(theta)) + 
		(0.5 * g * z1 * cos(theta) * z1 * cos(theta) * spacing * 
		((ps - pw) * ka + pw) * 
		cos(phi) * tan(phi));
	dnslope = (cl * spacing * z32 * cos(theta)) + 
		(0.5 * g * z3 * cos(theta) * z3 * cos(theta) * spacing * 
		((ps - pw) * kp + pw) * 
		cos(phi) * tan(phi));

	// frc = ps * g * z * sin(theta) * cos(theta) * spacing * spacing / cos(theta);
	frc = ps * g * z * sin(theta) * spacing * spacing;
#endif

	fs = (basal + shalstab + llateral + rlateral + upslope + dnslope) / frc;
 
	//cerr << "basal = " << basal << ", shalstab = " << shalstab << "\n";
	//cerr << "llateral = " << llateral << ", rlateral = " << rlateral << "\n";
	//cerr << "upslope = " << upslope << ", dnslope = " << dnslope << "\n";
 	//cerr << "force = " << frc << ", FS = " << fs << "\n"; 

#ifdef PRINT_FS
	cerr << "\nka\tkp\tko\tbasal\tshalstab\tllateral\trlateral\tupslope\tdnslope\tsum\tforce\tfs\n" <<
		ka << "\t" << kp << "\t" << ko << "\t" << basal << "\t" << shalstab << "\t" << llateral <<
		"\t" << rlateral << "\t" << upslope << "\t" << dnslope << "\t" << 
		basal + shalstab + llateral + rlateral + upslope + dnslope << "\t" << 
		frc << "\t" << fs << "\n";
#endif

        return(fs);
}

// version that computes FS based on the definition of a rigid block in the control grid
double Grid::get3dFS(int i, int j, int slope_field, int aspect_field, int soil_field, 
	int cb_field, int cl_field, int ctrl_field, double root_depth, int sat_field, 
	double ps, double phi, double *resis, double *force)
{
	double theta, cl, cb, z, h, a1, b1, c1, d1, a2, b2, c2, d2,  
		z1, z2, z3, z4, z12, z22, z32, z42, n1, n2, n3, n4,  
		ka, kp, ko, fs, frc, res = 0, basal, shalstab, aspect, control, 
		llateral = 0, rlateral = 0, upslope = 0, dnslope = 0;
	const double g = 9.81, pw = 1000.0;

	// check data
        if ((control = at(i, j, ctrl_field)) == nodata ||
		(theta = at(i, j, slope_field)) == nodata ||
		(cl = at(i, j, cl_field)) == nodata ||
		(cb = at(i, j, cb_field)) == nodata ||
		(h = at(i, j, sat_field)) == nodata ||
		(z = at(i, j, soil_field)) == nodata ||
		(aspect = at(i, j, aspect_field)) == nodata ||
		(a1 = at(i, j + 1, soil_field)) == nodata ||
		(b1 = at(i + 1, j, soil_field)) == nodata ||
		(c1 = at(i, j - 1, soil_field)) == nodata ||
		(d1 = at(i - 1, j, soil_field)) == nodata)
                return(nodata);
	else if (z <= 0)
		return(DEF_FLAT);

	a2 = at(i, j + 1, ctrl_field);
	b2 = at(i + 1, j, ctrl_field);
	c2 = at(i, j - 1, ctrl_field);
	d2 = at(i - 1, j, ctrl_field);

	// z1 = upslope, z3 = downslope, z2 = left side (facing downslope), z4 = right side
	if ((aspect > 315 && aspect <= 360) || (aspect >= 0 && aspect <= 45))
	{
		z1 = c1; z2 = d1; z3 = a1; z4 = b1;
		n1 = c2; n2 = d2; n3 = a2; n4 = b2;
	}
	else if (aspect > 45 && aspect <= 135)
	{
		z1 = d1; z2 = a1; z3 = b1; z4 = c1;
		n1 = d2; n2 = a2; n3 = b2; n4 = c2;
	}
	else if (aspect > 135 && aspect <= 235)
	{
		z1 = a1; z2 = b1; z3 = c1; z4 = d1;
		n1 = a2; n2 = b2; n3 = c2; n4 = d2;
	}
	else if (aspect > 235 && aspect <= 315)
	{
		z1 = b1; z2 = c1; z3 = d1; z4 = a1;
		n1 = b2; n2 = c2; n3 = d2; n4 = a2;
	}
	else
	{
		cerr << "invalid aspect value: " << aspect << "\n";
		return(nodata);
	}	

	// check ranges
	if (rad2deg(theta) <= 0 || phi >= 90)
		return(DEF_FLAT);
	if (rad2deg(theta) >= 90 || phi <= 0)
		return(DEF_STEEP);
 
	// convert angles to radians
        phi = deg2rad(phi);

	// earth pressure coefficients
	ka = sqr(tan(deg2rad(45) - (phi / 2.0)));
	kp = sqr(tan(deg2rad(45) + (phi / 2.0)));
	ko = 1 - sin(phi);

	// limit soil columns to overlap and rooting depth
	// was z1 = min(z1, root_depth), changed to take the minimum overlapping height
	// then was z1 = ((v1 + v) / 2) - (((v1 - z1) + (v - z)) / 2);
	z1 = avg(z, z1);
	z2 = avg(z, z2);
	z3 = avg(z, z3);
	z4 = avg(z, z4);
	z12 = min(z1, root_depth);
	z22 = min(z2, root_depth);
	z32 = min(z3, root_depth);
	z42 = min(z4, root_depth);

	// everywhere inside the shape
#ifdef NEW_3DSTAB
	basal = (z <= root_depth ? (cb * spacing * spacing * sec(theta)) : 0);
	
	if (phi > theta)
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) + 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);
	else
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) - 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);
#else
	basal = (z <= root_depth ? (cb * spacing * spacing / cos(theta)) : 0);
	// shalstab = ((ps * g * z * spacing * spacing * cos(theta) * cos(theta) / cos(theta)) - 
	//	(pw * g * h * spacing * spacing * cos(theta) * cos(theta) / cos(theta))) * tan(phi);
	shalstab = ((ps * g * z * spacing * spacing * cos(theta)) - 
		(pw * g * h * spacing * spacing * cos(theta))) * tan(phi);
#endif

	res += basal;
	res += shalstab;

	// perimeter of the shape
	if (control != 2)
	{

		if (n2 == nodata)
		{
#ifdef NEW_3DSTAB
			llateral = (cl * spacing * z22) + 
				(0.5 * g * ko * (ps * z2 * z2 - pw * h * h) * spacing * cos(phi) * tan(phi));
#ifdef PRINT_FS
			cerr << "\nllateral = (cl * spacing * z22) + (0.5 * g * ko * (ps * z2 * z2 - pw * h * h) * spacing * cos(phi) * tan(phi));\n";
			cerr << llateral << " = (" << cl << " * " << spacing << " * " << z22 << ") + (0.5 * " << g << " * " << ko << " * (" << ps << " * " << z2 << " * " << z2 << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi) << " * " << tan(phi) << ");\n";
#endif
				
#else
			// llateral = (cl * spacing * z22 * cos(theta) / cos(theta)) + 
			// (0.5 * g * z2 * cos(theta) * z2 * cos(theta) * spacing / cos(theta) * 
			// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
			llateral = (cl * spacing * z22) + 
			(0.5 * g * z2 * cos(theta) * z2 * spacing  * 
			((ps - pw) * ko + pw) * cos(phi) * tan(phi));
#endif
			res += llateral;
#ifdef COUNT_SIDES
			lcount += spacing;
#endif
		}
		if (n4 == nodata)
		{
#ifdef NEW_3DSTAB
			rlateral = (cl * spacing * z42) + 
				(0.5 * g * ko * (ps * z4 * z4 - pw * h * h) * spacing * cos(phi) * tan(phi));
#ifdef PRINT_FS
			cerr << "\nrlateral = (cl * spacing * z42) + (0.5 * g * ko * (ps * z4 * z4 - pw * h * h) * spacing * cos(phi) * tan(phi));\n";
			cerr << rlateral << " = (" << cl << " * " << spacing << " * " << z42 << ") + (0.5 * " << g << " * " << ko << " * (" << ps << " * " << z4 << " * " << z4 << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi) << " * " << tan(phi) << ");\n";
#endif
#else
			// rlateral = (cl * spacing * z42 * cos(theta) / cos(theta)) + 
			// (0.5 * g * z4 * cos(theta) * z4 * cos(theta) * spacing / cos(theta) * 
			// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
			rlateral = (cl * spacing * z42) + 
			(0.5 * g * z4 * cos(theta) * z4 * spacing  * 
			((ps - pw) * ko + pw) * cos(phi) * tan(phi));
#endif
			res += rlateral;
#ifdef COUNT_SIDES
			rcount += spacing;
#endif
		}
		if (n1 == nodata)
		{
#ifdef NEW_3DSTAB
#ifdef SAME_FORCE
			upslope = (cl * spacing * z12) + 
				(0.5 * g * ko * (ps * z1 * z1 - pw * h * h) * spacing * cos(phi) * tan(phi));
#else
			upslope = cl * spacing * z12; 
#endif
#ifdef PRINT_FS
			cerr << "\nupslope = cl * spacing * z12;\n";
			cerr << upslope << " = " << cl << " * " << spacing << " * " << z12 << ";\n";
#endif
#else
			upslope = (cl * spacing * z12 * cos(theta)) + 
			(0.5 * g * z1 * cos(theta) * z1 * cos(theta) * spacing * 
			((ps - pw) * ka + pw) * 
			cos(phi) * tan(phi));
#endif
			res += upslope;
#ifdef COUNT_SIDES
			ucount += spacing;
#endif
		}
		if (n3 == nodata)
		{
#ifdef NEW_3DSTAB
#ifdef SAME_FORCE
			dnslope = (cl * spacing * z32) + 
				(0.5 * g * ko * (ps * z3 * z3 - pw * h * h) * spacing * cos(phi) * tan(phi));
#else
			dnslope = (cl * spacing * z32) + 
				(0.5 * kp * g * (ps * z3 * z3 - pw * h * h) * spacing * cos(phi - theta));
#endif
#ifdef PRINT_FS
			cerr << "\ndnslope = (cl * spacing * z32) + (0.5 * kp * g * (ps * z3 * z3 - pw * h * h) * spacing * cos(phi - theta));\n";
			cerr << dnslope << " = (" << cl << " * " << spacing << " * " << z32 << ") + (0.5 * " << kp << " * " << g << " * (" << ps << " * " << z3 << " * " << z3 << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi - theta) << ");\n";
#endif
#else
			dnslope = (cl * spacing * z32 * cos(theta)) + 
			(0.5 * g * z3 * cos(theta) * z3 * cos(theta) * spacing * 
			((ps - pw) * kp + pw) * 
			cos(phi) * tan(phi));
#endif
			res += dnslope;
#ifdef COUNT_SIDES
			dcount += spacing;
#endif
		}
	}

#ifdef NEW_3DSTAB
	frc = (ps * g * z * sin(theta) * spacing * spacing) +
		(0.5 * ka * g * (ps * z * z - pw * h * h) * spacing * cos(phi - theta));
#ifdef PRINT_FS
	cerr << "\nfrc = (ps * g * z * sin(theta) * spacing * spacing) + (0.5 * ka * g * (ps * z * z - pw * h * h) * spacing * cos(phi - theta));\n";
	cerr << frc << " = (" << ps << " * " << g << " * " << z << " * " << sin(theta) << " * " << spacing << " * " << spacing << ") + (0.5 * " << ka << " * " << g << " * (" << ps << " * " << z << " * " << z << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi - theta) << ");\n";
#endif
#else
	// frc = ps * g * z * sin(theta) * cos(theta) * spacing * spacing / cos(theta);
	frc = ps * g * z * sin(theta) * spacing * spacing;
#endif
	fs = (basal + shalstab + llateral + rlateral + upslope + dnslope) / frc;

	// cerr << "FS = " << fs << "\n"; 

	*resis += res;
	*force += frc;

#ifdef PRINT_FS
	cerr << "\nka\tkp\tko\tbasal\tshalstab\tllateral\trlateral\tupslope\tdnslope\tsum\tforce\tfs\n" <<
		ka << "\t" << kp << "\t" << ko << "\t" << basal << "\t" << shalstab << "\t" << llateral <<
		"\t" << rlateral << "\t" << upslope << "\t" << dnslope << "\t" << 
		basal + shalstab + llateral + rlateral + upslope + dnslope << "\t" << 
		frc << "\t" << fs << "\n";
#endif
 
	return(fs);

}

double Grid::get3dFS(int i, int j, int slope_field, int aspect_field, int soil_field, 
	int cb_field, int cl_field, int ctrl_field, double root_depth, double sat_ratio, 
	double ps, double phi, double *resis, double *force)
{
	double theta, cl, cb, z, h, a1, b1, c1, d1, a2, b2, c2, d2,  
		z1, z2, z3, z4, z12, z22, z32, z42, n1, n2, n3, n4,  
		ka, kp, ko, fs, frc, res = 0, basal, shalstab, aspect, control, 
		llateral = 0, rlateral = 0, upslope = 0, dnslope = 0;
	const double g = 9.81, pw = 1000.0;

	// check data
        if ((control = at(i, j, ctrl_field)) == nodata ||
		(theta = at(i, j, slope_field)) == nodata ||
		(cl = at(i, j, cl_field)) == nodata ||
		(cb = at(i, j, cb_field)) == nodata ||
		(z = at(i, j, soil_field)) == nodata ||
		(h = z * sat_ratio) < 0 ||
		(aspect = at(i, j, aspect_field)) == nodata ||
		(a1 = at(i, j + 1, soil_field)) == nodata ||
		(b1 = at(i + 1, j, soil_field)) == nodata ||
		(c1 = at(i, j - 1, soil_field)) == nodata ||
		(d1 = at(i - 1, j, soil_field)) == nodata)
                return(nodata);
	else if (z <= 0)
		return(DEF_FLAT);

	a2 = at(i, j + 1, ctrl_field);
	b2 = at(i + 1, j, ctrl_field);
	c2 = at(i, j - 1, ctrl_field);
	d2 = at(i - 1, j, ctrl_field);

	// z1 = upslope, z3 = downslope, z2 = left side (facing downslope), z4 = right side
	if ((aspect > 315 && aspect <= 360) || (aspect >= 0 && aspect <= 45))
	{
		z1 = c1; z2 = d1; z3 = a1; z4 = b1;
		n1 = c2; n2 = d2; n3 = a2; n4 = b2;
	}
	else if (aspect > 45 && aspect <= 135)
	{
		z1 = d1; z2 = a1; z3 = b1; z4 = c1;
		n1 = d2; n2 = a2; n3 = b2; n4 = c2;
	}
	else if (aspect > 135 && aspect <= 235)
	{
		z1 = a1; z2 = b1; z3 = c1; z4 = d1;
		n1 = a2; n2 = b2; n3 = c2; n4 = d2;
	}
	else if (aspect > 235 && aspect <= 315)
	{
		z1 = b1; z2 = c1; z3 = d1; z4 = a1;
		n1 = b2; n2 = c2; n3 = d2; n4 = a2;
	}
	else
	{
		cerr << "invalid aspect value: " << aspect << "\n";
		return(nodata);
	}	

	// check ranges
	if (rad2deg(theta) <= 0 || phi >= 90)
		return(DEF_FLAT);
	if (rad2deg(theta) >= 90 || phi <= 0)
		return(DEF_STEEP);
 
	// convert angles to radians
        phi = deg2rad(phi);

	// earth pressure coefficients
	ka = sqr(tan(deg2rad(45) - (phi / 2.0)));
	kp = sqr(tan(deg2rad(45) + (phi / 2.0)));
	ko = 1 - sin(phi);

	// limit soil columns to overlap and rooting depth
	// was z1 = min(z1, root_depth), changed to take the minimum overlapping height
	// then was z1 = ((v1 + v) / 2) - (((v1 - z1) + (v - z)) / 2);
	z1 = avg(z, z1);
	z2 = avg(z, z2);
	z3 = avg(z, z3);
	z4 = avg(z, z4);
	z12 = min(z1, root_depth);
	z22 = min(z2, root_depth);
	z32 = min(z3, root_depth);
	z42 = min(z4, root_depth);

	// cerr << "aspect = " << aspect << ", z1 = " << z1 << ", z2 = " << z2 << ", z3 = " << z3 << ", z4 = " << z4 << "\n";

	// everywhere inside the shape
#ifdef NEW_3DSTAB
	basal = (z <= root_depth ? (cb * spacing * spacing * sec(theta)) : 0);
	
	if (phi > theta)
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) + 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);
	else
		shalstab = (((ps * z - pw * h) * g * spacing * spacing) - 
		(0.5 * g * (ka -kp) * (ps * z * z - pw * h * h) * spacing * sin(phi - theta))) * tan(phi);
#else
	basal = (z <= root_depth ? (cb * spacing * spacing / cos(theta)) : 0);
	// shalstab = ((ps * g * z * spacing * spacing * cos(theta) * cos(theta) / cos(theta)) - 
	//	(pw * g * h * spacing * spacing * cos(theta) * cos(theta) / cos(theta))) * tan(phi);
	shalstab = ((ps * g * z * spacing * spacing * cos(theta)) - 
		(pw * g * h * spacing * spacing * cos(theta))) * tan(phi);
#endif
	res += basal;
	res += shalstab;

	// perimeter of the shape
	if (control != 2)
	{

		if (n2 == nodata)
		{
#ifdef NEW_3DSTAB
			llateral = (cl * spacing * z22) + 
				(0.5 * g * ko * (ps * z2 * z2 - pw * h * h) * spacing * cos(phi) * tan(phi));
#ifdef PRINT_FS
			cerr << "\nllateral = (cl * spacing * z22) + (0.5 * g * ko * (ps * z2 * z2 - pw * h * h) * spacing * cos(phi) * tan(phi));\n";
			cerr << llateral << " = (" << cl << " * " << spacing << " * " << z22 << ") + (0.5 * " << g << " * " << ko << " * (" << ps << " * " << z2 << " * " << z2 << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi) << " * " << tan(phi) << ");\n";
#endif
#else
			// llateral = (cl * spacing * z22 * cos(theta) / cos(theta)) + 
			// (0.5 * g * z2 * cos(theta) * z2 * cos(theta) * spacing / cos(theta) * 
			// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
			llateral = (cl * spacing * z22) + 
			(0.5 * g * z2 * cos(theta) * z2 * spacing  * 
			((ps - pw) * ko + pw) * cos(phi) * tan(phi));
#endif
			res += llateral;
		}
		if (n4 == nodata)
		{
#ifdef NEW_3DSTAB
			rlateral = (cl * spacing * z42) + 
				(0.5 * g * ko * (ps * z4 * z4 - pw * h * h) * spacing * cos(phi) * tan(phi));
#ifdef PRINT_FS
			cerr << "\nrlateral = (cl * spacing * z42) + (0.5 * g * ko * (ps * z4 * z4 - pw * h * h) * spacing * cos(phi) * tan(phi));\n";
			cerr << rlateral << " = (" << cl << " * " << spacing << " * " << z42 << ") + (0.5 * " << g << " * " << ko << " * (" << ps << " * " << z4 << " * " << z4 << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi) << " * " << tan(phi) << ");\n";
#endif
#else
			// rlateral = (cl * spacing * z42 * cos(theta) / cos(theta)) + 
			// (0.5 * g * z4 * cos(theta) * z4 * cos(theta) * spacing / cos(theta) * 
			// ((ps - pw) * ko + pw) * cos(phi) * tan(phi));
			rlateral = (cl * spacing * z42) + 
			(0.5 * g * z4 * cos(theta) * z4 * spacing  * 
			((ps - pw) * ko + pw) * cos(phi) * tan(phi));
#endif
			res += rlateral;
		}
		if (n1 == nodata)
		{
#ifdef NEW_3DSTAB
			upslope = cl * spacing * z12; 
#ifdef PRINT_FS
			cerr << "\nupslope = cl * spacing * z12;\n";
			cerr << upslope << " = " << cl << " * " << spacing << " * " << z12 << ";\n";
#endif
#else
			upslope = (cl * spacing * z12 * cos(theta)) + 
			(0.5 * g * z1 * cos(theta) * z1 * cos(theta) * spacing * 
			((ps - pw) * ka + pw) * 
			cos(phi) * tan(phi));
#endif
			res += upslope;
		}
		if (n3 == nodata)
		{
#ifdef NEW_3DSTAB
			dnslope = (cl * spacing * z32) + 
				(0.5 * kp * g * (ps * z3 * z3 - pw * h * h) * spacing * cos(phi - theta));
#ifdef PRINT_FS
			cerr << "\ndnslope = (cl * spacing * z32) + (0.5 * kp * g * (ps * z3 * z3 - pw * h * h) * spacing * cos(phi - theta));\n";
			cerr << dnslope << " = (" << cl << " * " << spacing << " * " << z32 << ") + (0.5 * " << kp << " * " << g << " * (" << ps << " * " << z3 << " * " << z3 << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi - theta) << ");\n";
#endif
#else
			dnslope = (cl * spacing * z32 * cos(theta)) + 
			(0.5 * g * z3 * cos(theta) * z3 * cos(theta) * spacing * 
			((ps - pw) * kp + pw) * 
			cos(phi) * tan(phi));
#endif
			res += dnslope;
		}
	}

#ifdef NEW_3DSTAB
	frc = (ps * g * z * sin(theta) * spacing * spacing) +
		(0.5 * ka * g * (ps * z * z - pw * h * h) * spacing * cos(phi - theta));
#ifdef PRINT_FS
	cerr << "\nfrc = (ps * g * z * sin(theta) * spacing * spacing) + (0.5 * ka * g * (ps * z * z - pw * h * h) * spacing * cos(phi - theta));\n";
	cerr << frc << " = (" << ps << " * " << g << " * " << z << " * " << sin(theta) << " * " << spacing << " * " << spacing << ") + (0.5 * " << ka << " * " << g << " * (" << ps << " * " << z << " * " << z << " - " << pw << " * " << h << " * " << h << ") * " << spacing << " * " << cos(phi - theta) << ");\n";
#endif
#else
	// frc = ps * g * z * sin(theta) * cos(theta) * spacing * spacing / cos(theta);
	frc = ps * g * z * sin(theta) * spacing * spacing;
#endif
	fs = (basal + shalstab + llateral + rlateral + upslope + dnslope) / frc;

	// cerr << "basal = " << basal << ", shalstab = " << shalstab << "\n";
	// cerr << "llateral = " << llateral << ", rlateral = " << rlateral << "\n";
	// cerr << "upslope = " << upslope << ", dnslope = " << dnslope << "\n";
 	// cerr << "force = " << frc << ", FS = " << fs << "\n"; 

	*resis += res;
	*force += frc;

#ifdef PRINT_FS
	cerr << "\nka\tkp\tko\tbasal\tshalstab\tllateral\trlateral\tupslope\tdnslope\tsum\tforce\tfs\n" <<
		ka << "\t" << kp << "\t" << ko << "\t" << basal << "\t" << shalstab << "\t" << llateral <<
		"\t" << rlateral << "\t" << upslope << "\t" << dnslope << "\t" << 
		basal + shalstab + llateral + rlateral + upslope + dnslope << "\t" << 
		frc << "\t" << fs << "\n";
#endif
 
	return(fs);

}

double Grid::getSlope(int i, int j, int field, int mode) 
// Our Method
{
	double dz_dx, dz_dy, slope, a, b, c, d, e, f, g, h, l;

 	if ((e = at(i,     j, field)) == nodata || 
		(d = at(i-1,   j, field)) == nodata || 
		(f = at(i+1,   j, field)) == nodata || 
		(h = at(i,   j-1, field)) == nodata || 
		(b = at(i,   j+1, field)) == nodata || 
		(c = at(i+1, j+1, field)) == nodata || 
		(g = at(i-1, j-1, field)) == nodata || 
		(a = at(i-1, j+1, field)) == nodata || 
		(l = at(i+1, j-1, field)) == nodata) 
	{
		return(nodata);
	}

	slope = (sqrt(sqr((f - d) / (2.0 * spacing)) +
		sqr((b - h) / (2.0 * spacing))) +
		sqrt(sqr((c - g) / (SQRT2_x2 * spacing)) +
		sqr((a - l) / (SQRT2_x2 * spacing)))) / 2.0;

	//NOTE: TAN_DEG makes no sense and the radians are used only in the first case
	switch (mode) 
	{
		case ATAN_RAD:	// theta, radians
			return(atan(slope));
		case TAN_RAD:	// tan(theta), radians
			return(slope);
		case SIN_RAD:	// sin(theta), radians
			return(sin(atan(slope)));
		case COS_RAD:	// cos(theta), radians
			return(cos(atan(slope)));
                case ATAN_DEG:  // theta, degrees
                        return(rad2deg(atan(slope)));
                case TAN_DEG:   // tan(theta), degrees
                        return(rad2deg(slope));
		default:	// theta, radians
			return(atan(slope));
	}
}

double Grid::getSlope2(int i, int j, int field, int mode)
// Arc/Info Method (Horn)
{
        double dz_dx, dz_dy, slope, a, b, c, d, e, f, g, h, l;
 
        if ((e = at(i,     j, field)) == nodata ||
                (d = at(i-1,   j, field)) == nodata ||
                (f = at(i+1,   j, field)) == nodata ||
                (h = at(i,   j-1, field)) == nodata ||
                (b = at(i,   j+1, field)) == nodata ||
                (c = at(i+1, j+1, field)) == nodata ||
                (g = at(i-1, j-1, field)) == nodata ||
                (a = at(i-1, j+1, field)) == nodata ||
                (l = at(i+1, j-1, field)) == nodata)
        {
                return(nodata);
        }
 
        dz_dx = ((a + 2 * d + g) - (c + 2 * f + l)) / (8 * xsize);
        dz_dy = ((a + 2 * b + c) - (g + 2 * h + l)) / (8 * ysize);
        slope = sqrt(sqr(dz_dx) + sqr(dz_dy));
 
        switch (mode)
        {
                case ATAN_RAD:  // theta, radians
                        return(atan(slope));
                case TAN_RAD:   // tan(theta), radians
                        return(slope);
                case SIN_RAD:   // sin(theta), radians
                        return(sin(atan(slope)));
                case COS_RAD:   // cos(theta), radians
                        return(cos(atan(slope)));
                case ATAN_DEG:  // theta, degrees
                        return(rad2deg(atan(slope)));
                case TAN_DEG:   // tan(theta), degrees
                        return(rad2deg(slope));
		default:	// theta, radians
			return(atan(slope));
        }
}
 
double Grid::getSlope3(int i, int j, int field, int mode) 
// Downhill slope
{
	double a, b, c, d, e, f, g, h, l, val, dsize, slope = 0.0;
	int ndir = 0;

 	if ((e = at(i,     j, field)) == nodata || 
		(d = at(i-1,   j, field)) == nodata || 
		(f = at(i+1,   j, field)) == nodata || 
		(h = at(i,   j-1, field)) == nodata || 
		(b = at(i,   j+1, field)) == nodata || 
		(c = at(i+1, j+1, field)) == nodata || 
		(g = at(i-1, j-1, field)) == nodata || 
		(a = at(i-1, j+1, field)) == nodata || 
		(l = at(i+1, j-1, field)) == nodata) 
	{
		return(nodata);
	}

	dsize = sqrt(sqr(xsize) + sqr(ysize));
	
	if ((val = e - d) > 0.0)
	{
		slope += val / xsize;
		ndir++;
	}

	if ((val = e - f) > 0.0)
	{
		slope += val / xsize;
		ndir++;
	}

	if ((val = e - b) > 0.0)
	{
		slope += val / ysize;
		ndir++;
	}

	if ((val = e - h) > 0.0)
	{
		slope += val / ysize;
		ndir++;
	}

	if ((val = e - a) > 0.0)
	{
		slope += val / dsize;
		ndir++;
	}

	if ((val = e - c) > 0.0)
	{
		slope += val / dsize;
		ndir++;
	}

	if ((val = e - g) > 0.0)
	{
		slope += val / dsize;
		ndir++;
	}

	if ((val = e - l) > 0.0)
	{
		slope += val / dsize;
		ndir++;
	}
	
	slope /= ndir;

	switch (mode)
	{
                case ATAN_RAD:  // theta, radians
                        return(atan(slope));
                case TAN_RAD:   // tan(theta), radians
                        return(slope);
                case SIN_RAD:   // sin(theta), radians
                        return(sin(atan(slope)));
                case COS_RAD:   // cos(theta), radians
                        return(cos(atan(slope)));
                case ATAN_DEG:  // theta, degrees
                        return(rad2deg(atan(slope)));
                case TAN_DEG:   // tan(theta), degrees
                        return(rad2deg(slope));
		default:	// theta, radians
			return(atan(slope));
	}
}

double Grid::getSlope4(int i, int j, int field, int mode) 
// positive average (1D use)
{
	double a, b, c, d, e, f, g, h, l, val, dsize, slope = 0.0;
	int ndir = 0;

 	if ((e = at(i, j, field)) == nodata) 
		return(nodata);

	dsize = sqrt(sqr(xsize) + sqr(ysize));
	
	if ((d = at(i-1, j, field)) != nodata)
	{
		val = e - d;
		slope += abs(val) / xsize;
		ndir++;
	}

	if ((f = at(i+1, j, field)) != nodata)
	{
		val = e - f;
		slope += abs(val) / xsize;
		ndir++;
	}

	if ((b = at(i, j+1, field)) != nodata)
	{
		val = e - b;
		slope += abs(val) / ysize;
		ndir++;
	}

	if ((h = at(i, j-1, field)) != nodata)
	{
		val = e - h;
		slope += abs(val) / ysize;
		ndir++;
	}

	if ((a = at(i-1, j+1, field)) != nodata)
	{
		val = e - a;
		slope += abs(val) / dsize;
		ndir++;
	}

	if ((c = at(i+1, j+1, field)) != nodata)
	{
		val = e - c;
		slope += abs(val) / dsize;
		ndir++;
	}

	if ((g = at(i-1, j-1, field)) != nodata)
	{
		val = e - g;
		slope += abs(val) / dsize;
		ndir++;
	}

	if ((l = at(i+1, j-1, field)) != nodata)
	{
		val = e - l;
		slope += abs(val) / dsize;
		ndir++;
	}

	if (!ndir)
		return(nodata);
	
	slope /= ndir;

	switch (mode)
	{
                case ATAN_RAD:  // theta, radians
                        return(atan(slope));
                case TAN_RAD:   // tan(theta), radians
                        return(slope);
                case SIN_RAD:   // sin(theta), radians
                        return(sin(atan(slope)));
                case COS_RAD:   // cos(theta), radians
                        return(cos(atan(slope)));
                case ATAN_DEG:  // theta, degrees
                        return(rad2deg(atan(slope)));
                case TAN_DEG:   // tan(theta), degrees
                        return(rad2deg(slope));
		default:	// theta, radians
			return(atan(slope));
	}
}

double Grid::getSlope5(int i, int j, int field, int mode) 
// Our Method modified for nodata in the 3x3 neighborhood (use center cell instead)
{
	double dz_dx, dz_dy, slope, a, b, c, d, e, f, g, h, l;

 	if ((e = at(i,     j, field)) == nodata)
		return(nodata);
	if ((a = at(i-1, j+1, field)) == nodata)
		a = e;
	if ((b = at(i,   j+1, field)) == nodata)
		b = e;
	if ((c = at(i+1, j+1, field)) == nodata)
		c = e;
	if ((d = at(i-1,   j, field)) == nodata)
		d = e;
	if ((f = at(i+1,   j, field)) == nodata)
		f = e;
	if ((g = at(i-1, j-1, field)) == nodata)
		g = e;
	if ((h = at(i,   j-1, field)) == nodata)
		h = e;
	if ((l = at(i+1, j-1, field)) == nodata)
		l = e;

	slope = (sqrt(sqr((f - d) / (2.0 * spacing)) +
		sqr((b - h) / (2.0 * spacing))) +
		sqrt(sqr((c - g) / (SQRT2_x2 * spacing)) +
		sqr((a - l) / (SQRT2_x2 * spacing)))) / 2.0;

	switch (mode) //NOTE: TAN_DEG makes no sense
	{
		case ATAN_RAD:	// theta, radians
			return(atan(slope));
		case TAN_RAD:	// tan(theta), radians
			return(slope);
		case SIN_RAD:	// sin(theta), radians
			return(sin(atan(slope)));
		case COS_RAD:	// cos(theta), radians
			return(cos(atan(slope)));
                case ATAN_DEG:  // theta, degrees
                        return(rad2deg(atan(slope)));
                case TAN_DEG:   // tan(theta), degrees
                        return(rad2deg(slope));
		default:	// theta, radians
			return(atan(slope));
	}
}

double Grid::getAspect(int i, int j, int field) 
// from Zebenbergen and Thorne
{
	double aspect, a, b, c, d, e, f, g, h, l, tmp;

 	if ((e = at(i,     j, field)) == nodata || 
		(d = at(i-1,   j, field)) == nodata || 
		(f = at(i+1,   j, field)) == nodata || 
		(h = at(i,   j-1, field)) == nodata || 
		(b = at(i,   j+1, field)) == nodata)
	{
		return(nodata);
	}

	if (d == f && h == b)
		aspect = nodata;
	else if (d == f)
		aspect = (b < h ? 180 : 0);
	else if (b == h)
		aspect = (f < d ? 90 : 270);
	else
	{
		tmp = rad2deg(atan2((b - h) / (2.0 * spacing), (d - f) / (2.0 * spacing)));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}

	return(aspect);
}

Direction Grid::getAspect2(int i, int j, int field, int mode) 
// from Zebenbergen and Thorne, returns direction
{
	double aspect, a, b, c, d, e, f, g, h, l, tmp;

 	if ((e = at(i,     j, field)) == nodata || 
		(d = at(i-1,   j, field)) == nodata || 
		(f = at(i+1,   j, field)) == nodata || 
		(h = at(i,   j-1, field)) == nodata || 
		(b = at(i,   j+1, field)) == nodata)
	{
		return(NODIR);
	}

	// compute aspect
	if (d == f && h == b)
		return(NODIR);
	else if (d == f)
		aspect = (b < h ? 180 : 0);
	else if (b == h)
		aspect = (f < d ? 90 : 270);
	else
	{
		tmp = rad2deg(atan2((b - h) / (2.0 * spacing), (d - f) / (2.0 * spacing)));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}

	return(mode ? deg2dir(aspect, i, j) : deg2dir(aspect));
}

double Grid::getAspect3(int i, int j, int field) 
// from Zebenbergen and Thorne, modified for nodata in the 3x3 neighborhood (use center cell instead)
{
	double aspect, a, b, c, d, e, f, g, h, l, tmp;

 	if ((e = at(i,     j, field)) == nodata)
		return(nodata);
	if ((b = at(i,   j+1, field)) == nodata)
		b = e;
	if ((d = at(i-1,   j, field)) == nodata)
		d = e;
	if ((f = at(i+1,   j, field)) == nodata)
		f = e;
	if ((h = at(i,   j-1, field)) == nodata)
		h = e;

	if (d == f && h == b)
		aspect = nodata;
	else if (d == f)
		aspect = (b < h ? 180 : 0);
	else if (b == h)
		aspect = (f < d ? 90 : 270);
	else
	{
		tmp = rad2deg(atan2((b - h) / (2.0 * spacing), (d - f) / (2.0 * spacing)));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}
	return(aspect);
}

Direction Grid::getAspect4(int i, int j, int field, int mode) 
// from Zebenbergen and Thorne, returns direction, modified for nodata in the 3x3 neighborhood (use center cell instead)
{
	double aspect, a, b, c, d, e, f, g, h, l, tmp;

 	if ((e = at(i,     j, field)) == nodata)
		return(NODIR);
	if ((b = at(i,   j+1, field)) == nodata)
		b = e;
	if ((d = at(i-1,   j, field)) == nodata)
		d = e;
	if ((f = at(i+1,   j, field)) == nodata)
		f = e;
	if ((h = at(i,   j-1, field)) == nodata)
		h = e;

	// compute aspect
	if (d == f && h == b)
		return(NODIR);
	else if (d == f)
		aspect = (b < h ? 180 : 0);
	else if (b == h)
		aspect = (f < d ? 90 : 270);
	else
	{
		tmp = rad2deg(atan2((b - h) / (2.0 * spacing), (d - f) / (2.0 * spacing)));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}

	return(mode ? deg2dir(aspect, i, j) : deg2dir(aspect));
}

double Grid::getAspect5(int i, int j, int field)
// Arc/Info Method (Horn) (UNTESTED)
{
        double dz_dx, dz_dy, aspect, a, b, c, d, e, f, g, h, l, tmp;
 
        if ((e = at(i,     j, field)) == nodata ||
                (d = at(i-1,   j, field)) == nodata ||
                (f = at(i+1,   j, field)) == nodata ||
                (h = at(i,   j-1, field)) == nodata ||
                (b = at(i,   j+1, field)) == nodata ||
                (c = at(i+1, j+1, field)) == nodata ||
                (g = at(i-1, j-1, field)) == nodata ||
                (a = at(i-1, j+1, field)) == nodata ||
                (l = at(i+1, j-1, field)) == nodata)
        {
                return(nodata);
        }
 
        dz_dx = ((c + 2 * f + l) - (a + 2 * d + g)) / (8 * xsize);
        dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * ysize);

	// compute aspect
	if (dz_dx == 0 && dz_dy == 0)
		aspect = nodata;
	else if (dz_dx == 0)
		aspect = (dz_dy > 0 ? 0 : 180);
	else if (dz_dy == 0)
		aspect = (dz_dx > 0 ? 270 : 90);
	else
	{
		tmp = rad2deg(atan2(dz_dy, dz_dx));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}
	
	return(aspect);

}

Direction Grid::getAspect6(int i, int j, int field, int mode)
// Arc/Info Method (Horn), returns direction (UNTESTED)
{
        double dz_dx, dz_dy, aspect, a, b, c, d, e, f, g, h, l, tmp;
 
        if ((e = at(i,     j, field)) == nodata ||
                (d = at(i-1,   j, field)) == nodata ||
                (f = at(i+1,   j, field)) == nodata ||
                (h = at(i,   j-1, field)) == nodata ||
                (b = at(i,   j+1, field)) == nodata ||
                (c = at(i+1, j+1, field)) == nodata ||
                (g = at(i-1, j-1, field)) == nodata ||
                (a = at(i-1, j+1, field)) == nodata ||
                (l = at(i+1, j-1, field)) == nodata)
        {
                return(NODIR);
        }
 
        dz_dx = ((c + 2 * f + l) - (a + 2 * d + g)) / (8 * xsize);
        dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * ysize);

	// compute aspect
	if (dz_dx == 0 && dz_dy == 0)
		return(NODIR);
	else if (dz_dx == 0)
		aspect = (dz_dy > 0 ? 0 : 180);
	else if (dz_dy == 0)
		aspect = (dz_dx > 0 ? 270 : 90);
	else
	{
		tmp = rad2deg(atan2(dz_dy, dz_dx));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}

	return(mode ? deg2dir(aspect, i, j) : deg2dir(aspect));
}

double Grid::getAspect7(int i, int j, int field)
// Arc/Info Method (Horn), modified for nodata cases (UNTESTED)
{
        double dz_dx, dz_dy, aspect, a, b, c, d, e, f, g, h, l, tmp;
 
 	if ((e = at(i,     j, field)) == nodata)
		return(nodata);
        if ((d = at(i-1,   j, field)) == nodata)
		d = e;
	if ((f = at(i+1,   j, field)) == nodata)
		f = e;
	if ((h = at(i,   j-1, field)) == nodata)
		h = e;
	if ((b = at(i,   j+1, field)) == nodata)
		b = e;
	if ((c = at(i+1, j+1, field)) == nodata)
		c = e;
	if ((g = at(i-1, j-1, field)) == nodata)
		g = e;
	if ((a = at(i-1, j+1, field)) == nodata)
		a = e;
	if ((l = at(i+1, j-1, field)) == nodata)
		l = e;
 
        dz_dx = ((c + 2 * f + l) - (a + 2 * d + g)) / (8 * xsize);
        dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * ysize);

	// compute aspect
	if (dz_dx == 0 && dz_dy == 0)
		aspect = nodata;
	else if (dz_dx == 0)
		aspect = (dz_dy > 0 ? 0 : 180);
	else if (dz_dy == 0)
		aspect = (dz_dx > 0 ? 270 : 90);
	else
	{
		tmp = rad2deg(atan2(dz_dy, dz_dx));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}

	return(aspect);
}

Direction Grid::getAspect8(int i, int j, int field, int mode)
// Arc/Info Method (Horn), returns direction, modified for nodata cases (UNTESTED)
{
        double dz_dx, dz_dy, aspect, a, b, c, d, e, f, g, h, l, tmp;
 
 	if ((e = at(i,     j, field)) == nodata)
		return(NODIR);
        if ((d = at(i-1,   j, field)) == nodata)
		d = e;
	if ((f = at(i+1,   j, field)) == nodata)
		f = e;
	if ((h = at(i,   j-1, field)) == nodata)
		h = e;
	if ((b = at(i,   j+1, field)) == nodata)
		b = e;
	if ((c = at(i+1, j+1, field)) == nodata)
		c = e;
	if ((g = at(i-1, j-1, field)) == nodata)
		g = e;
	if ((a = at(i-1, j+1, field)) == nodata)
		a = e;
	if ((l = at(i+1, j-1, field)) == nodata)
		l = e;
 
        dz_dx = ((c + 2 * f + l) - (a + 2 * d + g)) / (8 * xsize);
        dz_dy = ((g + 2 * h + i) - (a + 2 * b + c)) / (8 * ysize);

	// compute aspect
	if (dz_dx == 0 && dz_dy == 0)
		return(NODIR);
	else if (dz_dx == 0)
		aspect = (dz_dy > 0 ? 0 : 180);
	else if (dz_dy == 0)
		aspect = (dz_dx > 0 ? 270 : 90);
	else
	{
		tmp = rad2deg(atan2(dz_dy, dz_dx));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}

	return(mode ? deg2dir(aspect, i, j) : deg2dir(aspect));
}

Direction Grid::getAspect9(int i, int j, int field, int mode) 
// from Zebenbergen and Thorne, returns direction, modified for nodata in the 3x3 neighborhood 
// like getAspect4 but instead of center cell uses the opposite cell (specific for tunnel of love)
{
	double aspect, a, b, c, d, e, f, g, h, l, tmp;

 	if ((e = at(i,     j, field)) == nodata)
	{
cerr << "getAspect: center cell is nodata (" << i << ", " << j << ") ...\n";
		return(NODIR);
	}

	b = at(i,   j+1, field);
	d = at(i-1,   j, field);
	f = at(i+1,   j, field);
	h = at(i,   j-1, field);

	if (b == nodata)
		b = h + RS_INCREMENT;
	if (d == nodata)
		d = f + RS_INCREMENT;
	if (f == nodata)
		f = d + RS_INCREMENT;
	if (h == nodata)
		h = b + RS_INCREMENT;

	// compute aspect
	if (d == f && h == b)
	{
// return random direction
		Direction d = NODIR;
		int n = RS_MAXPASSES;
		do
			d = (Direction)getRandInt(7);
		while (!isValue(i, j, d) && --n);
		return(d);
	}
	else if (d == f)
		aspect = (b < h ? 180 : 0);
	else if (b == h)
		aspect = (f < d ? 90 : 270);
	else
	{
		tmp = rad2deg(atan2((b - h) / (2.0 * spacing), (d - f) / (2.0 * spacing)));
		if (tmp > 90)
			aspect = 360 - tmp + 90;
		else
			aspect = 90 - tmp;		
	}

	return(mode ? deg2dir(aspect, i, j) : deg2dir(aspect));
}

double Grid::getCurvature(int i, int j, int field) 
{
	double a, b, c, d, e, f, g, h, l;

 	if ((e = at(i,     j, field)) == nodata || 
		(d = at(i-1,   j, field)) == nodata || 
		(f = at(i+1,   j, field)) == nodata || 
		(h = at(i,   j-1, field)) == nodata || 
		(b = at(i,   j+1, field)) == nodata || 
		(c = at(i+1, j+1, field)) == nodata || 
		(g = at(i-1, j-1, field)) == nodata || 
		(a = at(i-1, j+1, field)) == nodata || 
		(l = at(i+1, j-1, field)) == nodata) 
	{
		return(nodata);
	}

	return((2 * (f + d + b + h) + (c + a + g + l) - 12 * e) / 
		(4 * xsize * ysize));
}


double Grid::getMaxDHSlope(int i, int j, int field, Direction *dir, double *elev)
{
	int k;
	Direction d = NODIR;
	float slope = -1000000, val1, val2, elv, slp;

	if (isValue(i, j, field))
	{
		val1 = at(i, j, field);

		for (k = 0; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val2 = at(i, j, (Direction)k, field);
				if (val2 > val1)
					continue;

				slp = (val1 - val2) / 
					(isDiag((Direction)k) ? spacing * SQRT2 : spacing);
				if (slp > slope)
				{
					slope = slp;
					d = (Direction)k;
					elv = val2;
				}
			}
		}
	}


	if (dir != (Direction *)NULL)
		*dir = d;
	if (elev != (double *)NULL)
		*elev = elv;
	return(slope);  // negative means no downhill
}



double Grid::getMaxUHSlope(int i, int j, int field, Direction *dir, double *elev)
{
	int k;
	Direction d = NODIR;
	float slope = -1000000, val1, val2, elv, slp;

	if (isValue(i, j, field))
	{
		val1 = at(i, j, field);


		for (k = 0; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val2 = at(i, j, (Direction)k, field);
				if (val1 > val2)
					continue;
				slp = (val2 - val1) / 
					(isDiag((Direction)k) ? spacing * SQRT2 : spacing);
				if (slp > slope)
				{
					slope = slp;
					d = (Direction)k;
					elv = val2;
				}
			}
		}
	}


	if (dir != (Direction *)NULL)
		*dir = d;
	if (elev != (double *)NULL)
		*elev = elv;
	return(slope);  // negative means no uphill
}

Direction Grid::getMaxDHDir(int i, int j, int field)
{
	int k;
	Direction d = NODIR;
	double val1, val2;

	if (isValue(i, j, field))
	{
		val1 = at(i, j, field);

		for (k = 0; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val2 = at(i, j, (Direction)k, field);
				if (val2 >= val1)
					continue;
				d = (Direction)k;
				val1 = val2;
			}
		}
	}


	return(d);  // negative means no downhill
}

Direction Grid::getMaxDHDir2(int i, int j, int field)
// picks the lower elevation of the neighbors even if not downhill from center
{
	int k;
	Direction d = NODIR;
	double val1, val2;

	if (isValue(i, j, field))
	{
		// pick a starting neighboring elevation for comparison
		for (k = 0; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val1 = at(i, j, (Direction)k, field);
				d = (Direction)k;
				break;
			}
		}

		// any valid neighboring elevation?
		if (d == NODIR)
		{
			cerr << "getMaxDHDir2: all surrounding neighbors are NODATA\n";
			return(d);
		}

		for (k = 0, d = NODIR; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val2 = at(i, j, (Direction)k, field);
				if (val2 > val1)
					continue;
				d = (Direction)k;
				val1 = val2;
			}
		}
	}


	return(d);  // negative means no downhill
}

Direction Grid::getMaxUHDir(int i, int j, int field)
{
	int k;
	Direction d = NODIR;
	double val1, val2;

	if (isValue(i, j, field))
	{
		val1 = at(i, j, field);

		for (k = 0; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val2 = at(i, j, (Direction)k, field);
				if (val2 <= val1)
					continue;
				d = (Direction)k;
				val1 = val2;
			}
		}
	}


	return(d);  // negative means no uphill
}

Direction Grid::getMaxUHDir2(int i, int j, int field)
// picks the highest elevation of the neighbors even if not uphill from center
{
	int k;
	Direction d = NODIR;
	double val1, val2;

	if (isValue(i, j, field))
	{
		// pick a starting neighboring elevation for comparison
		for (k = 0; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val1 = at(i, j, (Direction)k, field);
				d = (Direction)k;
				break;
			}
		}

		// any valid neighboring elevation?
		if (d == NODIR)
		{
			cerr << "getMaxUHDir2: all surrounding neighbors are NODATA\n";
			return(d);
		}

		for (k = 0, d = NODIR; k < 8; ++k)
		{
			if (isValue(i, j, (Direction)k, field))
			{
				val2 = at(i, j, (Direction)k, field);
				if (val2 < val1)
					continue;
				d = (Direction)k;
				val1 = val2;
			}
		}
	}


	return(d);  // negative means no uphill
}


double Grid::getMaxRel(int i, int j, int efield, int vrfield, int rfield, int dfield)
// expects sink removed elevations, control grid with channels set to true
// and ridges set to false and the rest nodata, writes relief grid and optional distance
{
	Direction dir;
	double elv, slp, rel;
	int i0, j0;

	// nodata case
	if (!isValue(i, j, efield))
	{
		at(i, j, rfield) = nodata;
		if (dfield >= 0)
			at(i, j, dfield) = nodata;
		return(nodata);
	}

	// already computed
	if (isValue(i, j, rfield))
	{
		return(at(i, j, rfield));
	}

	// ridge case
	if (isValue(i, j, vrfield) && at(i, j, vrfield) == FALSE)
	{
		at(i, j, rfield) = 0;
		if (dfield >= 0)
			at(i, j, dfield) = 0;
		return(0);
	}

	// find steepset ascent
	slp = getMaxUHSlope(i, j, efield, &dir, &elv);

	// at the top
	if (slp < 0)
	{
		at(i, j, rfield) = 0;
		if (dfield >= 0)
			at(i, j, dfield) = 0;
		return(0);
	}

	// start recursion
	i0 = get_i(i, dir);
	j0 = get_j(j, dir);
	rel = getMaxRel(i0, j0, efield, vrfield, rfield, dfield);

	// add this cell
	rel += elv - at(i, j, efield);
	at(i, j, rfield) = rel;

	// keep track of distance (optional)
	if (dfield >= 0)
		at(i, j, dfield) = at(i0, j0, dfield) + 
			isDiag(dir) ? spacing * SQRT2 : spacing;
	return(rel);
}

Direction Grid::getFlowDir(int i, int j, int field)
{
	return(getMaxDHDir(i, j, field));
}

double Grid::chDist2Out(int i, int j, int flow_field, int dist_field)
{
	Direction fdir = (Direction)(at(i, j, flow_field)); 
	double dist = at(i, j, dist_field);

	// check: shouldn't it be 0 to 7??
	if (fdir < (Direction)1 || fdir > (Direction)8)	// can't follow
	{
		at(i, j, dist_field) = 0;
		return(0);
	}

	if (dist == nodata)		// follow (not yet computed)
		dist = chDist2Out(get_i(i, fdir), get_j(j, fdir), 
			flow_field, dist_field);
	dist += (isDiag(fdir) ? spacing * SQRT2 : spacing);
	at(i, j, dist_field) = dist;
	return(dist);
}


//
// Input/Output
//

void Grid::coutGrid(int field, int header, double zmin, double zmax) 
{
	cout.precision(15);
	switch (header)
	{
		case ARCINFO:
			cout << "ncols "	<< xdim		<< "\n";
			cout << "nrows "	<< ydim		<< "\n";
			cout << "xllcorner "	<< xorg		<< "\n";
			cout << "yllcorner "	<< yorg		<< "\n";
			cout << "cellsize "	<< spacing	<< "\n";
			cout << "NODATA_value " << nodata	<< "\n";
			break;

		case SURFER:
			if (zmin == 0.0 && zmax == 0.0)
			{
				zmin = findMin(field);
				zmax = findMax(field);
#if (GR_MESSAGES >= 1)
				cerr << "min: " << zmin << ", max: " << 
					zmax << "\n";
#endif
			}
			cout << "DSAA\n";
			cout << xdim << SEPARATOR << ydim				 << "\n";
			cout << xorg << SEPARATOR << xorg + xdim * xsize - 1 << "\n";
			cout << yorg << SEPARATOR << yorg + ydim * ysize - 1 << "\n";
			cout << zmin << SEPARATOR << zmax				 << "\n";
			break;

		default:
			break;
	}
	cout.precision(8); 

	for(int j = 0; j < ydim; ++j)  
	{
		for(int i = 0; i < xdim; ++i) 
		{
			cout << at(i, j, field) << "\t";
		}
		cout << "\n";	
	}
}

void Grid::cerrGrid(int field, int header, double zmin, double zmax)
{
        cerr.precision(15);
        switch (header)
        {
                case ARCINFO:
                        cerr << "ncols "        << xdim         << "\n";
                        cerr << "nrows "        << ydim         << "\n";
                        cerr << "xllcorner "    << xorg         << "\n";
                        cerr << "yllcorner "    << yorg         << "\n";
                        cerr << "cellsize "     << spacing      << "\n";
                        cerr << "NODATA_value " << nodata       << "\n";
                        break;
 
                case SURFER:
                        if (zmin == 0.0 && zmax == 0.0)
                        {
                                zmin = findMin(field);
                                zmax = findMax(field);
                        }
                        cerr << "DSAA\n";
                        cerr << xdim << SEPARATOR << ydim << "\n";
                        cerr << xorg << SEPARATOR << xorg + xdim * xsize - 1 << "\n";
                        cerr << yorg << SEPARATOR << yorg + ydim * ysize - 1 << "\n";
                        cerr << zmin << SEPARATOR << zmax << "\n";
                        break;
 
                default:
                        break;
        }
        cerr.precision(8);
 
        for(int j = 0; j < ydim; ++j)
        {
                for(int i = 0; i < xdim; ++i)
                {
                        cerr << at(i, j, field) << "\t";
                }
                cerr << "\n";
        }
}
 
void Grid::coutPoints(int field) 
{
	int i, j, n;

	cout.precision(15);

	// >= 0: write a single field
	if (field >= 0)
	{
		for(j = 0; j < ydim; ++j)  
		{
			for(i = 0; i < xdim; ++i) 
			{
				cout << xorg + i * xsize << SEPARATOR <<
					yorg + j * ysize << SEPARATOR <<
					at(i, j, field) << "\n";
			}
		}
	}
	// -1: write all fields
	else
	{
		for(j = 0; j < ydim; ++j)  
		{
			for(i = 0; i < xdim; ++i) 
			{
				cout << xorg + i * xsize << SEPARATOR <<
					yorg + j * ysize << SEPARATOR;
				for(n = 0; n < fields; ++n) 
					cout << at(i, j, n) << SEPARATOR;
				cout << "\n";
			}
		}
	}
	cout.precision(8); 
}

int Grid::writeGrid(char* filename, int field, int header, 
	double zmin, double zmax) 
{
	//  Saves a Grid to the Unix File "filename".
	ofstream outfile(filename);
	if (!outfile) 
	{
#if (GR_MESSAGES >= 1)
		cerr << "writeGrid: Can't open " << filename << " ...\n";
		return(0);
#endif
	}

	outfile.precision(15);
	switch (header)
	{
		case ARCINFO:
			outfile << "ncols "		<< xdim		<< "\n";
			outfile << "nrows "		<< ydim		<< "\n";
			outfile << "xllcorner "		<< xorg		<< "\n";
			outfile << "yllcorner "		<< yorg		<< "\n";
			outfile << "cellsize "		<< spacing	<< "\n";
			outfile << "NODATA_value "	<< nodata	<< "\n";
			break;

		case SURFER:
			if (zmin == 0.0 && zmax == 0.0)
			{
				zmin = findMin(field);
				zmax = findMax(field);
#if (GR_MESSAGES >= 1)
				cerr << "min: " << zmin << ", max: " << 
					zmax << "\n";
#endif
			}
			outfile << "DSAA\n";
			outfile << xdim << SEPARATOR << ydim				<< "\n";
			outfile << xorg << SEPARATOR << xorg + xdim * xsize - 1 << "\n";
			outfile << yorg << SEPARATOR << yorg + ydim * ysize - 1 << "\n";
			outfile << zmin << SEPARATOR << zmax				<< "\n";
			break;

		default:
			break;
	}
	outfile.precision(8);  

	for(int j = 0; j < ydim; j++)  
	{
		for(int i = 0; i < xdim; i++) 
		{
			outfile << at(i, j, field) << "\t";
		}
		outfile << "\n";
	}
	outfile.close();
	return(1);
}
 
int Grid::readGrid(char* filename, int field, int header) 
{
	//  Reads in a Grid from the Unix file "filename"
	double tmp, txm, tym;

	ifstream infile(filename);
	if (!infile) 
	{
#if (GR_MESSAGES >= 1)
		cerr << "readGrid: Can't open " << filename << " ...\n";
#endif
		return(0);

	}

	char tempstr[80];

	infile.precision(20);
	infile.width(20);

	if (!field)
	{
		switch (header)
		{
			case ARCINFO:
				infile >> tempstr >> xdim;
				infile >> tempstr >> ydim;
				infile >> tempstr >> xorg;
				infile >> tempstr >> yorg;
				infile >> tempstr >> spacing;
				infile >> tempstr >> nodata;
				xsize = ysize = spacing;
				break;

			case SURFER:
				infile >> tempstr;
				infile >> xdim >> ydim;
				infile >> xorg >> txm;
				infile >> yorg >> tym;
				infile >> tmp  >> tmp;
				xsize = ((txm - xorg + 1) / xdim);
				ysize = ((tym - yorg + 1) / ydim);
				spacing = (xsize + ysize) / 2.0;
				nodata = SURF_VAL;
				break;

			default:
				return(0);
		}

		// Initialize all the data structures to the appropriate size
		if (!setSize(ydim,xdim))
			return(0); 
	}
	else
	{
		int txdim, tydim;
		double txorg, tyorg, tspacing, tnodata, txsize, tysize;

		switch (header)
		{
			case ARCINFO:
				infile >> tempstr >> txdim;
				infile >> tempstr >> tydim;
				infile >> tempstr >> txorg;
				infile >> tempstr >> tyorg;
				infile >> tempstr >> tspacing;
				infile >> tempstr >> tnodata;
				txsize = tysize = tspacing;
				break;

			case SURFER:
				infile >> tempstr;
				infile >> txdim >> tydim;
				infile >> txorg >> txm;
				infile >> tyorg >> tym;
				infile >> tmp   >> tmp;
				txsize = ((txm - txorg + 1) / txdim);
				tysize = ((tym - tyorg + 1) / tydim);
				tspacing = (txsize + tysize) / 2.0;
				tnodata = nodata;
				break;

			default:
				return(0);
		}
		if (txdim != xdim || tydim != ydim ||
			txorg != xorg || tyorg != yorg ||
			txsize != xsize || tysize != ysize ||
			tspacing != spacing || tnodata != nodata)
		{
#if (GR_MESSAGES >= 1)
			cerr << "readGrid: grid " << filename << 
				" has a different header ...\n";
#endif
			return(0);
		}
	}

	// Read In The Data
	for(int j = 0; j < ydim; j++) 
	{
		for(int i = 0; i < xdim; i++) 
		{
			infile >> at(i, j, field);
		}
	}
	infile.close();
	if (fields <= 0)
		fields = 1;
	return(1);
}

//
// Initialization
//

int Grid::initGrid(int rows, int cols, int init, double csize, 
	double value, int field) 
{
	if (setSize(rows, cols, init))
	{
		setValue(value, field);
		spacing = xsize = ysize = csize;
		timeseed();
		return(1);
	}
	return(0);
}

int Grid::initGrid(Grid &z, double value, int field) 
{
	if (setSize(z))
	{
		setValue(value, field);
		timeseed();
		return(1);
	}
	return(0);
}

int Grid::initGrid(Grid &z) 
{
	int i, j, k;
	if (setSize(z))
	{
		for (k = 0; k < fields; ++k)
	  		for (i = 0; i < xdim; ++i) 
	    			for (j = 0; j < ydim; ++j)
	      				at(i,j,k) = z.at(i,j,k);
		timeseed();
		return(1);
	}
	return(0);
}

//
// Access
//

double& Grid::at(int i, int j, int field) 
{
	if (i >= xdim || i < 0 || j >= ydim || j < 0 || field < 0 || field >= fields)
		return(nodata);
	return(arrValues[field][i + (j * xdim)]);
}

double& Grid::at(int field) 
{
	return(arrValues[field][xcur + (ycur * xdim)]);
}

double& Grid::at(Direction dir, int field) 
{
	return(at(xcur, ycur, dir, field));
}

double& Grid::at(int i, int j, Direction dir, int field) 
{
	switch (dir) 
  	{
  		case R:
    			return(at(i+1, j, field));
    			break;
  		case L:
    			return(at(i-1, j, field));
    			break;
  		case U: 
    			return(at(i, j+1, field));
    			break;
  		case D:
    			return(at(i, j-1, field));
    			break;
  		case UR: 
    			return(at(i+1, j+1, field));
    			break;
  		case DL:
    			return(at(i-1, j-1, field));
    			break;
  		case DR: 
    			return(at(i+1, j-1, field));
    			break;
  		case UL:
    			return(at(i-1, j+1, field));
    			break;
		default:
#if (GR_MESSAGES >= 1)
			cerr << "at: invalid direction ...\n";
#endif
			return(at(i, j, field));
			break;
  	}
}

double& Grid::get(int i, int j, Direction dir, int field) 	// For compatibility
{
	return(at(i, j, dir, field));
}

double& Grid::goTo(int i, int j, int field) 
{
	setCurPos(i, j); 
	return(at(field));
}

double& Grid::goTo(Direction dir, int field) 
{
	return(goTo(xcur, ycur, dir, field));
}

double& Grid::goTo(int i, int j, Direction dir, int field) 
{
	switch (dir) 
  	{
  		case R:
    		return(goTo(i+1, j));
		case L:
			return(goTo(i-1, j));
		case U: 
			return(goTo(i, j+1));
		case D:
    		return(goTo(i, j-1));
		case UR: 
    		return(goTo(i+1, j+1));
		case DL:
			return(goTo(i-1, j-1));
		case DR: 
			return(goTo(i+1, j-1));
		case UL:
			return(goTo(i-1, j+1));
		default:
#if (GR_MESSAGES >= 1)
			cerr << "goTo: invalid direction ...\n";
			return(nodata);
			break;
#endif
	}
}

double& Grid::goToPrev(int field) 
{
	// ATTENTION: do we need a setcurpos since we have prev2?
	setCurPos(xcur, ycur);
	return(goTo(xprev1, yprev1, field));
}

int Grid::get_i(int i, Direction direc) 
{
  	switch (direc) 
  	{
  		case R:
    		return(i+1);
    		break;
  		case L:
    		return(i-1);
    		break;
  		case U: 
    		return(i);
    		break;
  		case D:
    		return(i);
    		break;
  		case UR: 
    		return(i+1);
    		break;
  		case DL:
    		return(i-1);
    		break;
  		case DR: 
    		return(i+1);
    		break;
  		case UL:
    		return(i-1);
    		break;
		default:
#if (GR_MESSAGES >= 1)
			cerr << "get_i: invalid direction ...\n";
			return(ERROR);
			break;
#endif
  	}
}

int Grid::get_i(Direction direc) 
{
  	switch (direc) 
  	{
  		case R:
    		return(xcur+1);
    		break;
  		case L:
    		return(xcur-1);
    		break;
  		case U: 
    		return(xcur);
    		break;
  		case D:
    		return(xcur);
    		break;
  		case UR: 
    		return(xcur+1);
    		break;
  		case DL:
    		return(xcur-1);
    		break;
  		case DR: 
    		return(xcur+1);
    		break;
  		case UL:
    		return(xcur-1);
    		break;
#if (GR_MESSAGES >= 1)
				cerr << "get_i: invalid direction ...\n";
				return(ERROR);
				break;
#endif
  	}
}

int Grid::get_j(int j, Direction direc) 
{
  	switch (direc) 
  	{
  		case R:
    		return(j);
    		break;
  		case L:
    		return(j);
    		break;
  		case U: 
    		return(j+1);
    		break;
  		case D:
    		return(j-1);
    		break;
  		case UR: 
    		return(j+1);
    		break;
  		case DL:
    		return(j-1);
    		break;
  		case DR: 
    		return(j-1);
    		break;
  		case UL:
    		return(j+1);
   			break;
#if (GR_MESSAGES >= 1)
			cerr << "get_j: invalid direction ...\n";
			return(ERROR);
			break;
#endif
  	}
}

int Grid::get_j(Direction direc) 
{
  	switch (direc) 
  	{
  		case R:
    		return(ycur);
    		break;
  		case L:
    		return(ycur);
    		break;
  		case U: 
    		return(ycur+1);
    		break;
  		case D:
    		return(ycur-1);
    		break;
  		case UR: 
    		return(ycur+1);
    		break;
  		case DL:
    		return(ycur-1);
    		break;
  		case DR: 
    		return(ycur-1);
    		break;
  		case UL:
    		return(ycur+1);
    		break;
#if (GR_MESSAGES >= 1)
			cerr << "get_j: invalid direction ...\n";
			return(ERROR);
			break;
#endif
  	}
}

int Grid::saveCurPos(POS *pos) 
{
	if (pos == (POS *)NULL)
	{
#if (GR_MESSAGES >= 1)
		cerr << "storeCurPos: position not stored!\n";
#endif
		return(0);	// must deal with error messages
	}
	pos->xcur = xcur; 
	pos->ycur = ycur;
	pos->xprev1 = xprev1; 
	pos->yprev1 = yprev1; 
	pos->xprev2 = xprev2;
	pos->yprev2 = yprev2;
	return(1);
}

int Grid::restoreCurPos(POS *pos) 
{
	if (pos == (POS *)NULL)
	{
#if (GR_MESSAGES >= 1)
		cerr << "restoreCurPos: position not restored!\n";
#endif
		return(0);	// must deal with error messages
	}
	xcur = pos->xcur; 
	ycur = pos->ycur;
	xprev1 = pos->xprev1; 
	yprev1 = pos->yprev1; 
	xprev2 = pos->xprev2;
	yprev2 = pos->yprev2;
	return(1);
}

int Grid::setCurPos(int i, int j) 
{
	if (i < 0 || i >= xdim || j < 0 || j >= ydim)
		return(0);	// must deal with error messages
	xprev2 = xprev1;
	yprev2 = yprev1;
	xprev1 = xcur; 
	yprev1 = ycur; 
	xcur = i; 
	ycur = j;
	return(1);
}

int Grid::setCurPos(Grid *grid) 
{
	if (grid->xcur < 0 || grid->xcur >= xdim || 
		grid->ycur < 0 || grid->ycur >= ydim)
		return(0);	// must deal with error messages
	xprev2 = xprev1;
	yprev2 = yprev1;
	xprev1 = xcur; 
	yprev1 = ycur; 
	xcur = grid->xcur; 
	ycur = grid->ycur;
	return(1);
}

double Grid::getLocMean(int i, int j, int field)
// returns mean val in 8 directions (not including center!)
{
        double val = 0.0;
        int h, n = 0;

        for (h = 0; h < 8; ++h)
        {
                if (isValue(i, j, (Direction)h, field))
		{
                        val += at(i, j, (Direction)h, field);
                        ++n;
		}
        }
        return(n ? val / n : nodata);
}

double Grid::getLocMin(int i, int j, int field)	// returns min val in 8 directions
{
	double val, val0 = at(i, j, field);
	int h;

	for (h = 0; h < 8; ++h)
	{
		if (isValue(i, j, (Direction)h, field) &&
			(val = at(i, j, (Direction)h, field)) < val0)
			val0 = val;
	}
	return(val);
}

double Grid::getLocMin(int field)
{
	return(getLocMin(xcur, ycur, field));
}

double Grid::getLocMax(int i, int j, int field) // returns max val in 8 directions
{
        double val, val0 = at(i, j, field);
        int h;
 
        for (h = 0; h < 8; ++h)
        {
                if (isValue(i, j, (Direction)h, field) &&
                        (val = at(i, j, (Direction)h, field)) > val0)
                        val0 = val;
        }
        return(val);
}

double Grid::getLocMax(int field)
{
	return(getLocMax(xcur, ycur, field));
}

double Grid::getWinRel(int i, int j, int dist, int field) // returns max relief
{
        double min, max, val; 
        int a, b, c, d, n, m;

	min = max = at(i, j, field);
	a = max(0, j - dist);
	b = max(0, i - dist);
	c = min(j + dist, ydim);
	d = min(i + dist, xdim);
 
        for (n = a; n <= c; ++n)
	{
        	for (m = b; m <= d; ++m)
        	{
                	if (isValue(m, n, field))
			{
				val = at(m, n, field);
                        	if (val > max)
					max = val;
				else if (val < min)
					min = val;
			}
		}
        }
        return(max - min);
}

double Grid::getWinRel(int i, int j, int dist, int field, 
	WIN *prev, WIN *cur) // returns max relief (faster)
{
        double val; 
        int xmin, xmax, ymin, ymax, minmax, m, n;

	cur->min = cur->max = at(i, j, field);
	cur->xmin = max(0, i - dist);
	cur->ymin = max(0, j - dist);
	cur->xmax = min(i + dist, xdim);
	cur->ymax = min(j + dist, ydim);

	// check overlap
	if (((cur->xmax <= prev->xmax && cur->xmax >= prev->xmin) ||
		(cur->xmin <= prev->xmax && cur->xmin >= prev->xmin)) &&
		((cur->ymax <= prev->ymax && cur->ymax >= prev->ymin) ||
		(cur->ymin <= prev->ymax && cur->ymin >= prev->ymin)))
	{
		xmin = max(cur->xmin, prev->xmin);
		ymin = max(cur->ymin, prev->ymin);
		xmax = min(cur->xmax, prev->xmax);
		ymax = min(cur->ymax, prev->ymax);

		// check max and min positions
		if (prev->mini >= xmin && 
			prev->mini <= xmax &&
			prev->minj >= ymin &&
			prev->minj <= ymin &&
			prev->maxi >= xmin && 
			prev->maxi <= xmax &&
			prev->maxj >= ymin &&
			prev->maxj <= ymin)
		{
			cur->max = prev->max;
			cur->min = prev->min;
			minmax = TRUE;
		}
		else
			minmax = FALSE;
	}
	else
		minmax = FALSE;

	// special case
	if (minmax)
	{
		// case 1
		for (n = cur->ymin; n <= ymin; ++n)
		{
	   	     	for (m = cur->xmin; m <= xmin; ++m)
	    	    	{
   	        	     	if (isValue(m, n, field))
				{
					val = at(m, n, field);
	                        	if (val > cur->max)
					{
						cur->max = val;
						cur->maxi = m;
						cur->maxj = n;
					}
					else if (val < cur->min)
					{
						cur->min = val;
						cur->mini = m;
						cur->minj = n;
					}
				}
			}
		}

		// case 2
		for (n = ymax; n <= cur->ymax; ++n)
		{
	   	     	for (m = cur->xmin; m <= xmin; ++m)
	    	    	{
   	        	     	if (isValue(m, n, field))
				{
					val = at(m, n, field);
	                        	if (val > cur->max)
					{
						cur->max = val;
						cur->maxi = m;
						cur->maxj = n;
					}
					else if (val < cur->min)
					{
						cur->min = val;
						cur->mini = m;
						cur->minj = n;
					}
				}
			}
	        }

		// case 3
		for (n = cur->ymin; n <= ymin; ++n)
		{
	   	     	for (m = xmax; m <= cur->xmin; ++m)
	    	    	{
   	        	     	if (isValue(m, n, field))
				{
					val = at(m, n, field);
	                        	if (val > cur->max)
					{
						cur->max = val;
						cur->maxi = m;
						cur->maxj = n;
					}
					else if (val < cur->min)
					{
						cur->min = val;
						cur->mini = m;
						cur->minj = n;
					}
				}
			}
		}
		// case 4
		for (n = ymax; n <= cur->ymax; ++n)
		{
	   	     	for (m = xmax; m <= cur->xmax; ++m)
	    	    	{
   	        	     	if (isValue(m, n, field))
				{
					val = at(m, n, field);
	                        	if (val > cur->max)
					{
						cur->max = val;
						cur->maxi = m;
						cur->maxj = n;
					}
					else if (val < cur->min)
					{
						cur->min = val;
						cur->mini = m;
						cur->minj = n;
					}
				}
			}
	        }
	}

	// general case
	else
	{
	        for (n = cur->ymin; n <= cur->ymax; ++n)
		{
	        	for (m = cur->xmin; m <= cur->xmax; ++m)
	        	{
	                	if (isValue(m, n, field))
				{
					val = at(m, n, field);
	                        	if (val > cur->max)
					{
						cur->max = val;
						cur->maxi = m;
						cur->maxj = n;
					}
					else if (val < cur->min)
					{
						cur->min = val;
						cur->mini = m;
						cur->minj = n;
					}
				}
			}
	        }
	}

        return(cur->max - cur->min);
}


double Grid::findMax(int field, int xinit, int yinit, int backw)
{
	int i, j;
	double value = nodata, maxval = -10.0e+100;

	if (!setCurPos(xinit, yinit))
		return(nodata);
	for (i = xcur; backw ? i >= 0 : i < xdim; backw ? --i : ++i)
	{
		for (j = ycur; backw ? j >= 0 : j < ydim; backw ? --j : ++j)
		{
			if ((value = at(i, j, field)) != nodata && 
				value > maxval)
			{
				maxval = value;
				setCurPos(i, j);
			}
		}
	}
	return(maxval);
}

/*
double Grid::findMax(int field, int *x, int *y, int xinit, int yinit, int backw)
{
	int i, j;
	double value = nodata, maxval = -10.0e+100;

	if (!setCurPos(xinit, yinit))
		return(nodata);
	for (i = xcur; backw ? i >= 0 : i < xdim; backw ? --i : ++i)
	{
		for (j = ycur; backw ? j >= 0 : j < ydim; backw ? --j : ++j)
		{
			if ((value = at(i, j, field)) != nodata && 
				value > maxval)
			{
				maxval = value;
				setCurPos(i, j);
				if (x != (int *)NULL && y != (int *)NULL)
				{
					*x = i; 
					*y = j;
				}
			}
		}
	}
	return(maxval);
}
*/

double Grid::findMax(int field, int *x, int *y)
{
	int i, j;
	double value = nodata, maxval = -10.0e+100;

	for (i = 0; i < xdim; ++i)
	{
		for (j = 0; j < ydim; ++j)
		{
			if ((value = at(i, j, field)) != nodata && 
				value > maxval)
			{
				maxval = value;
				if (x != (int *)NULL && y != (int *)NULL)
				{
					*x = i; 
					*y = j;
				}
			}
		}
	}
	return(maxval);
}

double Grid::findMin(int field, int xinit, int yinit, int backw)
{
	int i, j;
	double value = nodata, minval = 10.0e+100;

	if (!setCurPos(xinit, yinit))
		return(nodata);
	for (i = xcur; backw ? i >= 0 : i < xdim; backw ? --i : ++i)
	{
		for (j = ycur; backw ? j >= 0 : j < ydim; backw ? --j : ++j)
		{
			if ((value = at(i, j, field)) != nodata && 
				value < minval)
			{
				minval = value;
				setCurPos(i, j);
			}
		}
	}
	
	return(minval);
}

/*
double Grid::findMin(int field, int *x, int *y, int xinit, int yinit, int backw)
{
	int i, j;
	double value = nodata, minval = 10.0e+100;

	if (!setCurPos(xinit, yinit))
		return(nodata);
	for (i = xcur; backw ? i >= 0 : i < xdim; backw ? --i : ++i)
	{
		for (j = ycur; backw ? j >= 0 : j < ydim; backw ? --j : ++j)
		{
			if ((value = at(i, j, field)) != nodata && 
				value < minval)
			{
				minval = value;
				setCurPos(i, j);
				if (x != (int *)NULL && y != (int *)NULL)
				{
					*x = i; 
					*y = j;
				}
			}
		}
	}
	
	return(minval);
}
*/

double Grid::findMin(int field, int *x, int *y)
{
	int i, j;
	double value = nodata, minval = 10.0e+100;

	for (i = 0; i < xdim; ++i)
	{
		for (j = 0; j < ydim; ++j)
		{
			if ((value = at(i, j, field)) != nodata && 
				value < minval)
			{
				minval = value;
				if (x != (int *)NULL && y != (int *)NULL)
				{
					*x = i; 
					*y = j;
				}
			}
		}
	}
	
	return(minval);
}


int Grid::findCell(double value, int backw, int xinit, int yinit, int field)
{
	int i, j;

	if (xinit < 0 || xinit >= xdim || yinit < 0 || yinit >= ydim ||
		field < 0 || field >= fields)
		return(0);	// must deal with error messages

	for (i = xinit; backw ? i >= 0 : i < xdim; backw ? --i : ++i)
	{
		for (j = yinit; backw ? j >= 0 : j < ydim; backw ? --j : ++j)
		{
			if (at(i, j, field) == value)
			{
				return(setCurPos(i, j));
			}
		}
	}
	
	return(0);
}

int Grid::findNext(double value, int backw, int field)
{
	return(findCell(value, backw, xcur, ycur, field));
}

int Grid::isDiag(Direction dir, int field)
{
	return((dir == U || dir == D || dir == R || dir == L) ? 
		FALSE : TRUE);
}

int Grid::isValue(int field)
{
	return(at(field) == nodata ? 0 : 1);
}

int Grid::isValue(Direction dir, int field)
{
	return(at(dir, field) == nodata ? 0 : 1);
}

int Grid::isValue(int i, int j, int field)
{
	return(at(i, j, field) == nodata ? 0 : 1);
}

int Grid::isValue(int i, int j, Direction dir, int field)
{
	return(at(i, j, dir, field) == nodata ? 0 : 1);
}

//
// Operators
//

Grid& Grid::operator = (Grid& v) 
{
	int i, j, k;

  	xdim = v.xdim;
  	ydim = v.ydim;
  	xorg = v.xorg;
  	yorg = v.yorg;
  	spacing = v.spacing;
  	xsize = v.xsize;
  	ysize = v.ysize;
  	nodata = v.nodata;
	setSize(v);
	for (k = 0; k < fields; ++k)
  		for(i = 0; i < xdim; i++) 
    			for(j = 0; j < ydim; j++) 
      				at(i,j,k) = v.at(i,j,k);
  	return(*this);
}

Grid& Grid::operator = (double num) 
{
  	for(int i = 0; i < xdim; i++) 
    		for(int j = 0; j < ydim; j++) 
      			at(i,j) = num;
	return(*this);
}

Grid& Grid::operator -= (double num) 
{
  	for(int i = 0; i < xdim; i++) 
    		for(int j = 0; j < ydim; j++) 
      			at(i,j) = at(i,j) - num;
  	return(*this);
}

Grid& Grid::operator += (double num) 
{
  	for(int i = 0; i < xdim; i++) 
    		for(int j = 0; j < ydim; j++) 
      			at(i,j) = at(i,j) + num;
  	return(*this);
}

Grid& Grid::operator *= (double num) 
{
  	for(int i = 0; i < xdim; i++) 
    		for(int j = 0; j < ydim; j++) 
      			at(i,j) = at(i,j) * num;
  	return(*this);
}

Grid& Grid::operator /= (double num) 
{
  	for(int i = 0; i < xdim; i++) 
    		for(int j = 0; j < ydim; j++) 
      			at(i,j) = at(i,j) / num;
  	return(*this);
}

Grid& Grid::operator - (Grid& v) 
{
  	if(xdim != v.xdim || ydim != v.ydim) 
  	{
    		cerr << "Error- Attempt to perform subtraction on grids of unequal dimension.\n";
    		exit(1);
  	}
  	for (int k = 0; k < fields; ++k) 
	{
  		for (int i = 0; i < xdim; i++) 
  		{
    		for (int j = 0; j < ydim; j++) 
    		{
      			if (at(i,j) == nodata || v.at(i,j) == v.nodata) 
      			{
					at(i,j) = nodata;
      			}
      			else 
      			{
					at(i,j) = at(i,j) - v.at(i,j);
      			}
    		}
  		}
	}
  	return(*this);
}

Grid& Grid::operator + (Grid& v) 
{
  	if(xdim != v.xdim || ydim != v.ydim) 
  	{
    		cerr << "Error- Attempt to perform subtraction on grids of unequal dimension.\n";
    		exit(1);
  	}
	for (int k = 0; k < fields; ++k)
	{
  		for(int i = 0; i < xdim; i++) 
  		{
    		for(int j = 0; j < ydim; j++) 
    		{
      			if(at(i,j,k) == nodata || v.at(i,j,k) == v.nodata) 
      			{
					at(i,j,k) = nodata;
      			}
      			else 
      			{
					at(i,j,k) = at(i,j,k) + v.at(i,j,k);
      			}
    		}
  		}
	}
  	return(*this);
}

Grid& Grid::setBorders(double val, int field) 
{
  	for(int i = 0; i < xdim; i++) 
  	{
    	at(i, 0, field) = val;
    	at(i, ydim - 1, field) = val;
  	}
  	for(int j = 0; j < ydim; j++) 
  	{
    	at(0, j, field) = val;
    	at(xdim - 1, j, field) = val;
  	}
  	return(*this);
}

// Random number generation
#define RANDOM_IA		16807
#define RANDOM_IM		2147483647
#define RANDOM_AM		(1.0 / RANDOM_IM)
#define RANDOM_IQ		127773
#define RANDOM_IR		2836
#define RANDOM_RMASK	123459876

double Grid::dRandino()
{
	long l;
	double num;

	seedino ^= RANDOM_RMASK;
	l = seedino / RANDOM_IQ;
	seedino = RANDOM_IA * (seedino - l * RANDOM_IQ) - RANDOM_IR * l;
	if (seedino < 0L)
		seedino += RANDOM_IM;
	num = RANDOM_AM * seedino;
	seedino ^= RANDOM_RMASK;
	return(num);
}

int  Grid::getRandInt(int n)   
{
	return(int((dRandino()  * (double)n) + 0.5));
}

long Grid::getRandLong(long n) 
{
	return(long((dRandino() * (double)n) + 0.5));
}

int  Grid::doRndGrd(int hits, int range, int unique, int field)
{
	int n, rx, ry, rn;

	if (!hits || (unique && (hits > xdim * ydim)))
		hits = xdim * ydim;

	for (n = 0; n < hits; ++n)
	{
		// get a random number between 0 and xdim - 1
		rx = getRandInt(xdim - 1);
		// get a random number between 0 and ydim - 1
		ry = getRandInt(ydim - 1);
		// already assigned?
		if (unique)
		{
			if (isValue(rx, ry, field))
			{
				--n;
				continue;
			}
		}
		// get a random number between 1 and range
		rn = getRandInt(range - 1) + 1;
		// assign value
		at(rx, ry, field) = rn;
	}
	return(1);
}

int Grid::doRndGrd(int hits, int range, int unique, int infield, int outfield)
{
	int n, rx, ry, rn;

	for (n = 0; n < hits; ++n)
	{
		// get a random number between 0 and xdim - 1
		rx = getRandInt(xdim - 1);
		// get a random number between 0 and ydim - 1
		ry = getRandInt(ydim - 1);
		// already assigned?
		if (!(isValue(rx, ry, infield)) || 
			(unique && isValue(rx, ry, outfield)))
		{
			--n;
			continue;
		}
		// get a random number between 1 and range
		rn = getRandInt(range - 1) + 1;
		// assign value
		at(rx, ry, outfield) = rn;
	}
	return(1);
}

// Sink removal functions

void Grid::removeSinks(double inc, int maxpasses, int neighbors,
	int field, int disp)
{
	int cnt = 0, snk = 0;

	if (disp >= 1)
		cerr << "Removing sinks (max. passes: " <<
			maxpasses << ") ...\n";

	while (maxpasses--)
	{
		if (disp >= 1)
			cerr << "Removing sinks, pass number: " <<
				++cnt << "\n";

		snk = remSinksPass(inc, neighbors, field);
#if (RS_CHECKSINKS == 1)
 		if (!snk)
			break;
#else 
#if (RS_CHECKSINKS == 2)
		if (disp >= 1)
			cerr << "Checking for sinks ...\n";
		if (checkForSinks(field, FALSE))
			break;
#endif
#endif
	}
	if (disp >= 1)
		cerr << "Done removing sinks ...\n";
}

int Grid::remSinksPass(double addinc, int neighbors, int field, int disp) 
{
	// This is a simple Sink removal procedure.
	// Remove sinks by adding elevation to all sinks
	// untill there are no sinks left in the grid.
  	int i, j, stillSinks = 0, perc_comp = 0, last_comp = 0, 
		cnt = 0, tot = xdim * ydim;

	if (disp >= 1)
		cerr << "Completed with: " << perc_comp << " %\r";

    	for (j = 0; j < ydim; ++j) 
	{ 
      		for (i = 0; i < xdim; ++i) 
		{
			if (disp >= 1)
			{
				// Let the user know how far we are
				perc_comp = int(((double)(++cnt) /
					(double)tot) * 100.0);
				if (last_comp != perc_comp)
				{
					last_comp = perc_comp;
					cerr << "Completed with: " <<
						perc_comp << " %\r";
				}
			}
			if (isSink(i, j, field)) 
			{
	  			++stillSinks;
	 			at(i, j, field) = sinkAddValue(i, j, field) + addinc;
				if (neighbors)
	  				checkSinkNeighbors(i, j, addinc, field);
			}
      		}
    	}
	if (disp >= 1)
    		cerr << "\nNumber of sinks: " << stillSinks << "\n";

	return(stillSinks);
}

void Grid::checkSinkNeighbors(int i, int j, double addinc, int field) 
{
	if (isSink(i - 1, j, field)) 
	{
		at(i - 1, j, field) = sinkAddValue(i - 1, j, field) + addinc;
		checkSinkNeighbors(i - 1, j, addinc, field);
	}
	if (isSink(i, j - 1, field)) 
	{
		at(i, j - 1, field) = sinkAddValue(i, j - 1, field) + addinc;
		checkSinkNeighbors(i, j - 1, addinc, field);
	}
	if (isSink(i - 1, j - 1, field)) 
	{
		at(i - 1, j - 1, field) = sinkAddValue(i - 1, j - 1, field) + addinc;
		checkSinkNeighbors(i - 1, j - 1, addinc, field);
	}
	if (isSink(i + 1, j + 1, field)) 
	{
		at(i + 1, j + 1, field) = sinkAddValue(i + 1, j + 1, field) + addinc;
		checkSinkNeighbors(i + 1, j + 1, addinc, field);
	}
	if (isSink(i + 1, j, field)) 
	{
		at(i + 1, j, field) = sinkAddValue(i + 1, j, field) + addinc;
		checkSinkNeighbors(i + 1, j, addinc, field);
	}
	if (isSink(i, j + 1, field)) 
	{
		at(i, j + 1, field) = sinkAddValue(i, j + 1, field) + addinc;
		checkSinkNeighbors(i, j + 1,addinc, field);
	}
	if (isSink(i + 1, j - 1, field)) 
	{
		at(i + 1, j - 1, field) = sinkAddValue(i + 1, j - 1, field) + addinc;
		checkSinkNeighbors(i + 1, j - 1, addinc, field);
	}
	if (isSink(i - 1, j + 1, field)) 
	{
		at(i - 1, j + 1, field) = sinkAddValue(i - 1, j + 1, field) + addinc;
		checkSinkNeighbors(i - 1, j + 1, addinc, field);
	}
}


double Grid::sinkAddValue(int i, int j, int field) 
{
	double minneigh = get(i, j, L, field);
	for (int k = 0; k < 8; ++k) 
		if (get(i, j, (Direction)k, field) < minneigh) 
      			minneigh = get(i, j, (Direction)k, field);
  	return(minneigh);
}

int Grid::isSink(int i, int j, int field)
{
        if (at(i,     j, field) == nodata ||
                at(i-1,   j, field) == nodata ||
                at(i+1,   j, field) == nodata ||
                at(i,   j-1, field) == nodata ||
                at(i,   j+1, field) == nodata ||
                at(i+1, j+1, field)     == nodata ||
                at(i-1, j-1, field) == nodata ||
                at(i-1, j+1, field) == nodata ||
                at(i+1, j-1, field) == nodata)
        {
                return(0);
        }
 
        if (at(i, j, field) <= at(i-1,   j, field) &&
                at(i, j, field) <= at(i+1,   j, field) &&
                at(i, j, field) <= at(i,   j-1, field) &&
                at(i, j, field) <= at(i,   j+1, field) &&
                at(i, j, field) <= at(i+1, j+1, field) &&
                at(i, j, field) <= at(i-1, j+1, field) &&
                at(i, j, field) <= at(i+1, j-1, field) &&
                at(i, j, field) <= at(i-1, j-1, field))
        {
                return(1);
        }
        else
        {
                return(0);
        }
}


int Grid::checkForSinks(int field, int disp)
{
        int totsinks = 0, perc_comp = 0, last_comp = 0,
                cnt = 0, tot = xdim * ydim;
 
        for (int j = 0; j < ydim; j++)
        {
                for(int i = 0; i < xdim; i++)
                {
			if (disp >= 2)
			{
                        	// Let the user know how far we are
                        	perc_comp = int(((double)(++cnt) /
					(double)tot) * 100.0);

                        	if (last_comp != perc_comp)
                        	{
                               		last_comp = perc_comp;
                               	 	cerr << "Completed with: "
						<< perc_comp << " %\r";
                        	}
			}

                        if (isSink(i, j, field))
                        {
                                if (!disp)
                                        return(0);
                                ++totsinks;
                                if (disp == 3)
                                        cerr << "Sink found in elevation grid at (" << i
                                                << ", " << j << ") ...\n";
                        }
                }
        }
 
	if (disp >= 2)
		cerr << "\n";

        if (disp && totsinks)
        {
                cerr << "checkForSinks: " << totsinks <<
                        " sinks found in elevation grid!\n";
                return(0);
        }
 
        return(1);
}

// Sink removal functions for channels (modified to consider nodata as a "wall")

void Grid::removeSinks2(double inc, int maxpasses, int neighbors,
	int field, int disp)
{
	int cnt = 0, snk = 0;

	if (disp >= 1)
		cerr << "Removing sinks (max. passes: " <<
			maxpasses << ") ...\n";

	while (maxpasses--)
	{
		if (disp >= 1)
			cerr << "Removing sinks, pass number: " <<
				++cnt << "\n";

		snk = remSinksPass2(inc, neighbors, field);
#if (RS_CHECKSINKS == 1)
 		if (!snk)
			break;
#else 
#if (RS_CHECKSINKS == 2)
		if (disp >= 1)
			cerr << "Checking for sinks ...\n";
		if (checkForSinks2(field, FALSE))
			break;
#endif
#endif
	}
	if (disp >= 1)
		cerr << "Done removing sinks ...\n";
}

int Grid::remSinksPass2(double addinc, int neighbors, int field, int disp) 
{
	// This is a simple Sink removal procedure.
	// Remove sinks by adding elevation to all sinks
	// untill there are no sinks left in the grid.
  	int i, j, stillSinks = 0, perc_comp = 0, last_comp = 0, 
		cnt = 0, tot = xdim * ydim;

	if (disp >= 1)
		cerr << "Completed with: " << perc_comp << " %\r";

    	for (j = 0; j < ydim; ++j) 
	{ 
      		for (i = 0; i < xdim; ++i) 
		{
			if (disp >= 1)
			{
				// Let the user know how far we are
				perc_comp = int(((double)(++cnt) /
					(double)tot) * 100.0);
				if (last_comp != perc_comp)
				{
					last_comp = perc_comp;
					cerr << "Completed with: " <<
						perc_comp << " %\r";
				}
			}
			if (isSink2(i, j, field)) 
			{
	  			++stillSinks;
	 			at(i, j, field) = sinkAddValue(i, j, field) + addinc;
				if (neighbors)
	  				checkSinkNeighbors2(i, j, addinc, field);
			}
      		}
    	}
	if (disp >= 1)
    		cerr << "\nNumber of sinks: " << stillSinks << "\n";

	return(stillSinks);
}

void Grid::checkSinkNeighbors2(int i, int j, double addinc, int field) 
{
	if (isSink2(i - 1, j, field)) 
	{
		at(i - 1, j, field) = sinkAddValue(i - 1, j, field) + addinc;
		checkSinkNeighbors2(i - 1, j, addinc, field);
	}
	if (isSink2(i, j - 1, field)) 
	{
		at(i, j - 1, field) = sinkAddValue(i, j - 1, field) + addinc;
		checkSinkNeighbors2(i, j - 1, addinc, field);
	}
	if (isSink2(i - 1, j - 1, field)) 
	{
		at(i - 1, j - 1, field) = sinkAddValue(i - 1, j - 1, field) + addinc;
		checkSinkNeighbors2(i - 1, j - 1, addinc, field);
	}
	if (isSink2(i + 1, j + 1, field)) 
	{
		at(i + 1, j + 1, field) = sinkAddValue(i + 1, j + 1, field) + addinc;
		checkSinkNeighbors2(i + 1, j + 1, addinc, field);
	}
	if (isSink2(i + 1, j, field)) 
	{
		at(i + 1, j, field) = sinkAddValue(i + 1, j, field) + addinc;
		checkSinkNeighbors2(i + 1, j, addinc, field);
	}
	if (isSink2(i, j + 1, field)) 
	{
		at(i, j + 1, field) = sinkAddValue(i, j + 1, field) + addinc;
		checkSinkNeighbors2(i, j + 1,addinc, field);
	}
	if (isSink2(i + 1, j - 1, field)) 
	{
		at(i + 1, j - 1, field) = sinkAddValue(i + 1, j - 1, field) + addinc;
		checkSinkNeighbors2(i + 1, j - 1, addinc, field);
	}
	if (isSink2(i - 1, j + 1, field)) 
	{
		at(i - 1, j + 1, field) = sinkAddValue(i - 1, j + 1, field) + addinc;
		checkSinkNeighbors2(i - 1, j + 1, addinc, field);
	}
}

int Grid::checkForSinks2(int field, int disp)
{
        int totsinks = 0, perc_comp = 0, last_comp = 0,
                cnt = 0, tot = xdim * ydim;
 
        for (int j = 0; j < ydim; j++)
        {
                for(int i = 0; i < xdim; i++)
                {
			if (disp >= 2)
			{
                        	// Let the user know how far we are
                        	perc_comp = int(((double)(++cnt) /
					(double)tot) * 100.0);

                        	if (last_comp != perc_comp)
                        	{
                               		last_comp = perc_comp;
                               	 	cerr << "Completed with: "
						<< perc_comp << " %\r";
                        	}
			}

                        if (isSink2(i, j, field))
                        {
                                if (!disp)
                                        return(0);
                                ++totsinks;
                                if (disp == 3)
                                        cerr << "Sink found in elevation grid at (" << i
                                                << ", " << j << ") ...\n";
                        }
                }
        }
 
	if (disp >= 2)
		cerr << "\n";

        if (disp && totsinks)
        {
                cerr << "checkForSinks2: " << totsinks <<
                        " sinks found in elevation grid!\n";
                return(0);
        }
 
        return(1);
}

int Grid::isSink2(int i, int j, int field) // modified for edges
{
	int k, vals = 0;

	if (!isValue(i, j, field))
		return(0);

	for (k = 0; k < 8; ++k)
	{
		if (isValue(i, j, (Direction)k, field))
		{
			if (at(i, j, (Direction)k, field) < at(i, j, field))
				return(0);
			++vals;
		}
	}
	return(vals > 1 ? 1 : 0); // try to exclude outlet
}

// Other miscellaneous functions

// converts from degrees to direction
Direction Grid::deg2dir(double deg)
{
	// convert to direction
	if (deg >= 337.5 || deg < 22.5)
		return(U);
	else if (deg >= 22.5 && deg < 67.5)
		return(UR);
	else if (deg >= 67.5 && deg < 112.5)
		return(R);
	else if (deg >= 112.5 && deg < 157.5)
		return(DR);
	else if (deg >= 157.5 && deg < 202.5)
		return(D);
	else if (deg >= 202.5 && deg < 247.5)
		return(DL);
	else if (deg >= 247.5 && deg < 292.5)
		return(L);
	else if (deg >= 292.5 && deg < 337.5)
		return(UL);
	return(NODIR);
}

// converts from degrees to direction, returns the next best direction when nodata
Direction Grid::deg2dir(double deg, int i, int j)
{
	// convert to direction
	if (deg >= 337.5 || deg < 22.5)
	{
		if (isValue(i, j, U))
			return(U);
		else if (deg >= 337.5)
		{
			if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, D))
				return(D);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, D))
				return(D);
			else
				return(NODIR);
		}
	}
	else if (deg >= 22.5 && deg < 67.5)
	{
		if (isValue(i, j, UR))
			return(UR);
		else if (deg >= 45)
		{
			if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, DL))
				return(DL);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, DL))
				return(DL);
			else
				return(NODIR);
		}
	}
	else if (deg >= 67.5 && deg < 112.5)
	{
		if (isValue(i, j, R))
			return(R);
		else if (deg >= 90)
		{
			if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, L))
				return(L);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, L))
				return(L);
			else
				return(NODIR);
		}
	}
	else if (deg >= 112.5 && deg < 157.5)
	{
		if (isValue(i, j, DR))
			return(DR);
		else if (deg >= 135)
		{
			if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, UL))
				return(UL);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, UL))
				return(UL);
			else
				return(NODIR);
		}
	}
	else if (deg >= 157.5 && deg < 202.5)
	{
		if (isValue(i, j, D))
			return(D);
		else if (deg >= 180)
		{
			if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, U))
				return(U);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, U))
				return(U);
			else
				return(NODIR);
		}
	}
	else if (deg >= 202.5 && deg < 247.5)
	{
		if (isValue(i, j, DL))
			return(DL);
		else if (deg >= 225)
		{
			if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, UR))
				return(UR);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, UR))
				return(UR);
			else
				return(NODIR);
		}
	}
	else if (deg >= 247.5 && deg < 292.5)
	{
		if (isValue(i, j, L))
			return(L);
		else if (deg >= 270)
		{
			if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, R))
				return(R);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, UL))
				return(UL);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, DR))
				return(DR);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, R))
				return(R);
			else
				return(NODIR);
		}
	}
	else if (deg >= 292.5 && deg < 337.5)
	{
		if (isValue(i, j, UL))
			return(UL);
		else if (deg >= 315)
		{
			if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, DR))
				return(DR);
			else
				return(NODIR);
		}
		else
		{
			if (isValue(i, j, L))
				return(L);
			else if (isValue(i, j, U))
				return(U);
			else if (isValue(i, j, DL))
				return(DL);
			else if (isValue(i, j, UR))
				return(UR);
			else if (isValue(i, j, D))
				return(D);
			else if (isValue(i, j, R))
				return(R);
			else if (isValue(i, j, DR))
				return(DR);
			else
				return(NODIR);
		}
	}
	return(NODIR);
}

// converts from direction to degrees
double Grid::dir2deg(Direction dir)
{
	switch (dir)
	{
		case U:
			return(0);
		case UR:
			return(45);
		case R:
			return(90);
		case DR:
			return(135);
		case D:
			return(180);
		case DL:
			return(225);
		case L:
			return(270);
		case UL:
			return(315);
		default:
			return(DEF_VAL);
	}
}