#ifndef TerrainProcessor_H
#define TerrainProcessor_H
/*
* HydroPlus is source code and models developed and managed by Theodore Endreny and his students at SUNY ESF, te@esf.edu
* The US Forest Service and Davey Tree have provided strategic software development support through funding and personnel
* Attribution for use of software and code is requested. Please see Executive.cpp for documentation.
*/
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <sstream>
#include <utility>
#include <iomanip>
#include <map>
#include <limits>
#include <algorithm>   
#include <queue>
//_USE_MATH_DEFINES gives access to Pi as M_PI, and must precede #include <cmath> or <math.h>
#define _USE_MATH_DEFINES
#include <cmath>
#include <cfloat>

//using namespace std allows for use of standard namespace without the preface 
using namespace std;

//Forward declaration of class Inputs
class Inputs;

//ratio of degrees to radians
#define Ratio_Degree_to_Radian 180 / M_PI
//ratio of radians to degrees
#define Ratio_Radian_to_Degree M_PI / 180



//using namespace std;                    

struct Cell {
	int row;
	int col;
	double z;
};

struct GridCell {
	int row;
	int column;
	double priority;

	bool operator<(const GridCell& other) const {
		return priority > other.priority; // Use '>' for min-heap behavior
	}
};

struct GridCell2 {
	int row;
	int column;
	double z;
	double priority;

	bool operator<(const GridCell2& other) const {
		return priority < other.priority;
	}
};



class TerrainProcessor
{
public:

	TerrainProcessor(Inputs* input);
	void DEM_Read_Header(Inputs* input);
	void TopographicIndex_Histogram_Read(Inputs* input);
	void TopographicIndex_Map_Read(Inputs* input);
	double getAveSoilMoistureDeficit(Inputs* input);
	void TopographicIndex_Map_Create(Inputs* input);
	void TopographicIndex_Histogram_Create(Inputs* input);
	//TopographicCharacteristics_Vector_to_Map function converts DEM byproducts (TI, slope, aspect, FD, FA) from vectors to maps
	void TopographicCharacteristics_Vector_to_Map();
	void Ground_Slope_and_Aspect_File_Read();

	bool fileExists(const string& filename);
	        
	void getElevation(vector<vector<double>> &elevation_m_2D){elevation_m_2D = elevation_m_2D;}
	void getElevation1D(vector<double> &elevation_m_2D){elevation_m_2D = elevation_m_1D;}
	void getflowAccum1D(vector<double>& flowAccum) { flowAccum = flowAccum1D; }

	int getcols() { return nCols; }
	int getrows() { return nRows; }
	double getNODATA() { return NODATA_code; }
	double getxllcorner() { return xllcorner_m; }
	double getyllcorner() { return yllcorner_m; }
	void setcols(int x) { nCols = x; }
	void setrows(int x) { nRows = x; }
	void setxllcorner(double x) { xllcorner_m = x; }
	void setyllcorner(double x) { yllcorner_m = x; }

	static void calcFlowAccum_Dinf();
	static void calcFlowAccum_D8();
	static void calcFlowAccum_MFD();
	static void FlowAccum(vector<double>in_weight_raster_frac_vec, vector<double> in_weight_raster_vec);
	static void FlowAccum(vector<double> in_weight_raster_vec);
	static void preAccumWeighting(vector<double> weight_vector);
	static void postAccumWeighting(vector<double> weight_vector);
	static void flowAccumorganizer_initialize();
	
	static vector<vector<double>> FlowAccumulation_Map;// Dinf flow accumulation organizer
	static vector<double> flowAccum1D;
	static vector<double> slopeorganizer_riserun1D;
	//vector<double> slopeorganizer_Horn1D;
	static vector<vector<double>> FlowDirection_Map;	// Dinf flow direction organizer
	static vector<double> DX_vec;
	static vector<double> AreaFlowingToPervious_frac_vec;
	//elevation of each folder, elevation_m_2D(i,j)
	static vector<vector<double>> elevation_m_2D;
	static vector<double> elevation_m_1D, elevation_m_negative_1D;
	static void calcFlowDirection_DInf();
	static void calcFlowDirection_MFD();
	static void calcFlowDirection_D8();
	//static vector<double> flowDir1D, neg_flowDir1D;
	static vector<vector<pair<int, int>>> ContributingAreaPixels_FlowIntoAllPixels_vec; // 2D vector to store the contributing path [row, col] pairs for each cell
	static vector<vector<pair<int, int>>> DispersalAreaPixels_FlowFromAllPixels_vec;// 2D vector to store the dispersal path [row, col] pairs for each cell
	static double max_2d_element(vector<vector<double>> vec2D);
	static double min_2d_element(vector<vector<double>> vec2D);

	//checkfolder function checks if a folder exists, returns 1 if True, NODATA_code values if False
	static bool checkfolder(int row, int col);
	static void trackDispersalAreas_D8();
	static bool saveContributingDispersalArea;
	static void setContributingDispersalFlags();


	static vector<vector<double>> Slope_Ground_rise_run_mpm_Map, Slope_Ground_Horn_rad_Map;

	double TI_Value_Average = 0;
	double indexSize = 0;
	int indexAreaSize = 0;
	vector<double> slopeGround_1D_rad;
	vector<double> aspectGround_1D_N_0_rad;
	double stat_PC_frac;
	vector<double> TI_Area_frac;
	vector<double> TI_Value;

private:
	//Inputs pointer stored
	Inputs* input;

	static vector<vector<double>> NeighborAnalysis_Map;	// temporary organizer for the flow accumulation function
	static vector<vector<double>> InFlowDirection_Map;	// permanent organizer for the checking the flow accumulation function
	vector<vector<double>> TopographicIndex_Map;// final topographic wetness index  
	vector<vector<double>> Aspect_Ground_N_0_rad_Map;

	vector<double> FDorganizer_1D;
	//ATB array (i,j), as acronym for "a tangent Beta", the topographic index (exponential = ln(a/tanB) or power = (a/tanB)^n) 
	vector<vector<double>> ATB;
	//Area array (i,j) is upslope area draining per contour (pixel) width
	vector<vector<double>> Area;
	//TI_Values_In_HistogramBin is counter of index bins
	vector<double> TI_Values_In_HistogramBin;

	//nRows, nCols are size of map input; NATB is number of TI_Value bins and not used
	static int nRows, nCols, NATB;
	//xll and yll are lower left map coordinates in xllcorner_m (easting) and yllcorner_m (northing) directions
	double xllcorner_m, yllcorner_m;
	double cellsize_m = 0;

	string simulationMode;
	string routingAlgorithm;
	string Flag_TI_Pervious_Accumulated;
	string SimulationStringParams;

	// simulation method (0;1): 1: need to load impervious coverage data
	int simMethod;                      
	// Number of ln(a/tanB) or (a/tanB)^n histogram ordinates
	string Ascii_Header_str;                    
	// Length_Pixel_Side_m, organizer spacing in elevation_m_2D array, Area_Catchment_1_frac: fractional area of the subcatchment
	static double Length_Pixel_Side_m;
	static double Area_Catchment_1_frac;
	static double NODATA_code;
	// total fractional impervious area of the simulated subcatchment
	double totalFImpArea;
    string Directory_Input_CLArg;
	//NAC is number of topographic index bins
	int NAC;
	//currentDepth is parameter for flow accumulation routing, used in interation through DEM cells
	static int currentDepth;
	//maxDepth is parameter threshold for flow accumulation routing
	static int maxDepth;
	double subFarea = 0;
	static vector<vector<double>> elevation_working_m_2D, elevation_depression_filled_m_2D, elevation_flats_filled_m_2D;
	//vector <double>ImperviousCover_frac;
	//Functions with TerrainProcessor

	//init() function will set initial values for the following variables used in the topographic index
	void init() { nRows = 0; nCols = 0;  NATB = 0; Length_Pixel_Side_m = 0.0; Area_Catchment_1_frac = 1.0; totalFImpArea = 0.0; }
	void clearVectorsRelatedToTopographicIndex();
	//DEM_Read_and_Prepare function checks if DEM was read, returns 1 if True, 0 if False
	void DEM_Read_and_Prepare(Inputs* input);
	static double getElevation(int row, int col);
	void DEM_to_Slope_Calc();
	void FillDepressions(vector<vector<double>> elevation_working_m_2D);
	void FixFlats(vector<vector<double>> elevation_working_m_2D, double small_num);
	void Find_pit(vector<vector<double>> elevation_map);
	//DEM_to_FlowDirection_Manager computes DEM flow direction used for flow accumulation, needed in TI calculation
	void DEM_to_FlowDirection_Manager();
	//DEM_to_FlowAccumulation_Manager computes DEM flow accumulation needed in TI calculation
	void DEM_to_FlowAccumulation_Manager();
	//TopographicIndex_Map_Compute computes topographic index (TI) 
	void TopographicIndex_Map_Compute();
	//TopographicIndex_Map_to_Histogram function converts TI map to histogram (log space), with bins holding magnitude of each TI value
	void TopographicIndex_Map_to_Histogram();
	//TopographicIndex_Log_to_Power function converts topographic index values from log to power magnitude
	void TopographicIndex_Log_to_Power();

	//FlowAccumulation_DInfinity function is recursive Dinfinity flow accumulation algorithm 
	static void FlowAccumulation_DInfinity(int row, int col); 
	//FlowAccumulation_MultipleFlowDirection function is recursive MFD (multiple flow direction) flow accumulation algorithm
	static void FlowAccumulation_MultipleFlowDirection(int row, int col);  

	static vector<Cell> undefined_flow_cells;
	static vector<pair<int, int>> possible_outlets;
	static vector< vector <int>> flats; 
};

#endif

