//  -*- Mode: C++;  -*-
//  File: grid.h
//  Author: Dino Bellugi (dino@geomorph.berkeley.edu)
//  Copyright  Dino Bellugi, BlueG SoftWear, U.C. Berkeley, 1999 (C)
// *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// * FUNCTION: all the grid 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"
//
#define NEW_3DSTAB	// activates new model
//#define PRINT_FS	// prints values at each cell
#define COUNT_SIDES	// counts perimeter
#define SAME_FORCE	// lateral everywhere for testing


#ifndef	DEF_GRID
#define	DEF_GRID
#define CKCELLS		100
#define DEF_VAL		-9999.0
#define SURF_VAL	1.70141e+38
#define DEF_SIZE	1.0
#define DEF_STEEP	-10.0
#define DEF_FLAT	10.0
#define DEF_PSPW	1.7
#define DEF_PHI		45.0
#define TRUE 		1
#define FALSE 		0
#define GR_MESSAGES 	1

#include "gridedges.h"

#define DEFFIELDS	1	// default number of fields stored in grid
#define MAXFIELDS	20	// max number of fields stored in grid
#define DEF_FIELD	0
#define ARCINFO		0
#define SURFER		1
#define SEPARATOR	" "
#define ATAN_RAD	0
#define TAN_RAD		1
#define SIN_RAD		2
#define COS_RAD		3
#define ATAN_DEG	4
#define TAN_DEG		5
#define SLP_MODES	6

#ifndef DEF_REMOVESINKS
#define DEF_REMOVESINKS
#define RS_NEIGHBORS    1       // 0 = don't check, 1 = check neighbors
#define RS_MAXPASSES    100     // maximum number of passes
#define RS_INCREMENT    0.1     // soil increment
#define RS_MESSAGES     2       // message level
#define RS_CHECKSINKS   1       // 0 = don't check, 1 = test, 2 = check
#endif

#ifndef DEF_QOVERT
#define DEF_QOVERT
#define QT_PS		1600
#define QT_PHI		45
#endif

// Macros
#define srandino(l)	seedino = l
#define timeseed()	{time_t tt; time(&tt); srandino(tt);}
//#define getRandInt(n)	 int((dRandino() * (double)n) + 0.5)
//#define getRandLong(n) long((dRandino() * (double)n) + 0.5)

typedef struct position
{
	int xcur, ycur, xprev1, yprev1, xprev2, yprev2;
} POS;

typedef struct window
{
	int xmin, ymin, xmax, ymax, mini, minj, maxi, maxj;
	double min, max;
} WIN;

#define copywin(d,s)	{(d)->xmin=(s)->xmin;(d)->xmax=(s)->xmax; \
			(d)->ymin=(s)->ymin;(d)->ymax=(s)->ymax; \
			(d)->mini=(s)->mini;(d)->minj=(s)->minj; \
			(d)->maxi=(s)->maxi;(d)->maxj=(s)->maxj; \
			(d)->max=(s)->max;(d)->min=(s)->min;} \

class Grid 
{
private:
	int setSize(int i, int j, int init = FALSE, int flds = DEFFIELDS);
	int setSize(Grid &z);
	double *arrValues[MAXFIELDS];
	int xcur, ycur, xprev1, yprev1, xprev2, yprev2;
	long seedino;		// Random number seed
public:
	double xorg, yorg;
	double spacing, xsize, ysize, nodata;
	int xdim, ydim, fields;

	// Creation
	Grid(int x = 0, int y = 0, double xorg = 0.0, double yorg = 0.0, 
		double csize = DEF_SIZE, double nodat = DEF_VAL, 
		int flds = DEFFIELDS);
  	Grid(Grid &z);
  	~Grid();

	// Initialization
  	int initGrid(int nrows, int ncols, int init, double csize = DEF_SIZE, 
		double value = DEF_VAL, int field = 0);
  	int initGrid(Grid &z, double value, int field = 0);
  	int initGrid(Grid &z);
	void setNodata(double value);
	void chgNodata(double value, int set = 0, int field = 0);
	void chgValue(double old_value, double new_value, int field = 0);
	void chgValue(int nvalues, double *old_value, double new_value, int field = 0);
  	void setValue(double value = DEF_VAL, int field = 0);
  	void setValue(int field1, int field2);
  	void setValue(int field1, int field2, double value);
  	void setValue(int field1, double value = DEF_VAL, int field2 = 0);
  	void addValue(int field1, int field2);
  	void subValue(int field1, int field2);
  	int addField(int copy = -1);
  	int delField();
	long countValue(double value, int field = 0);
	long countValue(int field = 0);
	void mkBinary(double threshold, int reverse = 0, int infield = 0, int outfield = 0);

	// Access
  	double& at(int field = 0);
  	double& at(Direction dir, int field = 0);
  	double& at(int i, int j, int field = 0);
  	double& at(int i, int j, Direction dir, int field = 0);
  	double& get(int i, int j, Direction dir, int field = 0);
  	double& goTo(int i, int j, int field = 0);
  	double& goTo(Direction dir, int field = 0);
  	double& goTo(int i, int j, Direction dir, int field = 0);
  	double& goToPrev(int field = 0);
  	double getLocMin(int field = 0);
  	double getLocMin(int i, int j, int field = 0);
  	double getLocMean(int i, int j, int field = 0);
  	double getLocMax(int field = 0);
  	double getLocMax(int i, int j, int field = 0);
  	double getWinRel(int i, int j, int dist, int field = 0);
  	double getWinRel(int i, int j, int dist, int field,
		WIN *prev, WIN *cur);
  	double findMax(int field = 0, int xinit = 0, int yinit = 0,
		int backw = FALSE);
  	//double findMax(int field, int *x, int *y, 
	//	int xinit = 0, int yinit = 0, int backw = FALSE);
  	double findMax(int field, int *x, int *y); 
  	double findMin(int field = 0, int xinit = 0, int yinit = 0,
		int backw = FALSE);
  	//double findMin(int field, int *x, int *y, 
	//	int xinit = 0, int yinit = 0, int backw = FALSE);
  	double findMin(int field, int *x, int *y);
  	int saveCurPos(POS *pos);
  	int restoreCurPos(POS *pos);
  	int setCurPos(int i, int j);
  	int setCurPos(Grid *grid);
  	int findCell(double value, int backw = FALSE, int xinit = 0, 
		int yinit = 0, int field = 0);
  	int findNext(double value, int backw = FALSE, int field = 0);
 	int isValue(int field = 0);
  	int isDiag(Direction dir, int field = 0);
  	int isValue(Direction dir, int field = 0);
  	int isValue(int i, int j, int field = 0);
  	int isValue(int i, int j, Direction dir, int field = 0);
  	int nextValue(int field = 0);
	double getArrVal(int index, int field) {return(arrValues[field][index]);}

	// q/T, slope, curvature, distance
	double getQoverT(double theta, double aoverb, double ps = DEF_PSPW,
		double phi = DEF_PHI);
	double getQoverT(int i, int j, int theta_field, int aoverb_field,
        	double ps = DEF_PSPW, double phi = DEF_PHI);
	double getQoverT(int i, int j, int theta_field, double aoverb,
        	double ps = DEF_PSPW, double phi = DEF_PHI);
	double getQoverT(int i, int j, double theta, int aoverb_field,
        	double ps = DEF_PSPW, double phi = DEF_PHI);
	double getQoverT2(int i, int j, int theta_field, int aoverb_field,
        	double ps = DEF_PSPW, double phi = DEF_PHI);
	double 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 = DEF_PSPW, double phi = DEF_PHI);
 	double 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 = DEF_PSPW, double phi = DEF_PHI);
	double 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 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 getSlope(int i, int j, int field = 0, int mode = 0);
  	double getSlope2(int i, int j, int field = 0, int mode = 0);
  	double getSlope3(int i, int j, int field = 0, int mode = 0);
  	double getSlope4(int i, int j, int field = 0, int mode = 0);
  	double getSlope5(int i, int j, int field = 0, int mode = 0);
  	double getAspect(int i, int j, int field = 0);
  	Direction getAspect2(int i, int j, int field = 0, int mode = 0);
  	double getAspect3(int i, int j, int field = 0);
  	Direction getAspect4(int i, int j, int field = 0, int mode = 0);
  	double getAspect5(int i, int j, int field = 0);
  	Direction getAspect6(int i, int j, int field = 0, int mode = 0);
  	double getAspect7(int i, int j, int field = 0);
  	Direction getAspect8(int i, int j, int field = 0, int mode = 0);
  	Direction getAspect9(int i, int j, int field = 0, int mode = 0);
  	double getCurvature(int i, int j, int field = 0);
	double getMaxDHSlope(int i, int j, int field = 0, 
		Direction *dir = (Direction *)NULL, 
		double *elev = (double *)NULL);
	double getMaxUHSlope(int i, int j, int field = 0, 
		Direction *dir = (Direction *)NULL, 
		double *elev = (double *)NULL);
	Direction getMaxDHDir(int i, int j, int field = 0);
	Direction getMaxUHDir(int i, int j, int field = 0);
	Direction getMaxDHDir2(int i, int j, int field = 0);
	Direction getMaxUHDir2(int i, int j, int field = 0);
	double getMaxRel(int i, int j, int efield, int vrfield, 
		int rfield, int dfield = -1);
	Direction getFlowDir(int i, int j, int field = 0);
	double chDist2Out(int i, int j, int flow_field, int dist_field);

	// I.O.
  	int writeGrid(char *name, int field = 0, int header = ARCINFO,
		double zmin = 0.0, double zmax = 0.0);
  	void coutGrid(int field = 0, int header = ARCINFO,
		double zmin = 0.0, double zmax = 0.0);
  	void cerrGrid(int field = 0, int header = ARCINFO,
		double zmin = 0.0, double zmax = 0.0);
  	void coutPoints(int field = 0);
  	int readGrid(char *name, int field = 0, int header = ARCINFO);

  	int get_i(int i, Direction k);
  	int get_j(int j, Direction k);
  	int get_i(Direction k);
  	int get_j(Direction k);

	// Sink removal
	void removeSinks(double addinc, int maxpasses = RS_MAXPASSES,
		int neighbors = RS_NEIGHBORS, int field = 0,
		int disp = RS_MESSAGES);
  	int isSink(int i, int j, int field = 0); 
  	int checkForSinks(int field = 0, int disp = RS_MESSAGES);
	int remSinksPass(double addinc, int neighbors = RS_NEIGHBORS,
		int field = 0, int disp = RS_MESSAGES);
	void checkSinkNeighbors(int i, int j, double addinc, int field = 0);
	double sinkAddValue(int i, int j, int field = 0);

	// Sink removal for channels
	void removeSinks2(double addinc, int maxpasses = RS_MAXPASSES,
		int neighbors = RS_NEIGHBORS, int field = 0,
		int disp = RS_MESSAGES);
  	int isSink2(int i, int j, int field = 0); 
  	int checkForSinks2(int field = 0, int disp = RS_MESSAGES);
	int remSinksPass2(double addinc, int neighbors = RS_NEIGHBORS,
		int field = 0, int disp = RS_MESSAGES);
	void checkSinkNeighbors2(int i, int j, double addinc, int field = 0);

	// Operators
  	Grid& operator = (Grid& a);
  	Grid& operator = (double num);
  	Grid& operator -= (double num);
  	Grid& operator += (double num);
  	Grid& operator *= (double num);
  	Grid& operator /= (double num);
  	Grid& operator - (Grid& a);
  	Grid& operator + (Grid& a);

	// Random number generator
	double dRandino(void);
	int getRandInt(int);
	long getRandLong(long);
	int doRndGrd(int hits = 0, int range = 100, int unique = 0, 
		int field = 0);
	int doRndGrd(int hits, int range, int unique, 
		int infield, int outfield);

	// Other
  	Grid& setBorders(double val, int field = 0);
	Direction deg2dir(double deg);
	Direction deg2dir(double deg, int i, int j);
	double dir2deg(Direction dir);
};

#endif