#include <sys/stat.h>					 
#include "TerrainProcessor.h"
#include "Inputs.h"

#include "../Buffer/BufferSpatialCalc.h"

//all the static variables need to be redefined in cpp files
vector<vector<double>> TerrainProcessor::FlowDirection_Map, TerrainProcessor::FlowAccumulation_Map, TerrainProcessor::NeighborAnalysis_Map, TerrainProcessor::InFlowDirection_Map;
int TerrainProcessor::nRows, TerrainProcessor::nCols, TerrainProcessor::NATB;
double TerrainProcessor::Length_Pixel_Side_m, TerrainProcessor::Area_Catchment_1_frac, TerrainProcessor::NODATA_code;
vector<vector<double>> TerrainProcessor::elevation_m_2D, TerrainProcessor::elevation_working_m_2D, TerrainProcessor::elevation_depression_filled_m_2D, TerrainProcessor::elevation_flats_filled_m_2D;
int TerrainProcessor::currentDepth, TerrainProcessor::maxDepth;
vector<vector<double>> TerrainProcessor::Slope_Ground_Horn_rad_Map, TerrainProcessor:: Slope_Ground_rise_run_mpm_Map;// slope of each organizer folder
vector<double> TerrainProcessor::slopeorganizer_riserun1D;
vector<double> TerrainProcessor::DX_vec;
vector<double> TerrainProcessor::AreaFlowingToPervious_frac_vec;
vector<double> TerrainProcessor::elevation_m_1D, TerrainProcessor::elevation_m_negative_1D;
vector<double> TerrainProcessor::flowAccum1D;
vector<vector<pair<int, int>>>  TerrainProcessor::ContributingAreaPixels_FlowIntoAllPixels_vec, TerrainProcessor::DispersalAreaPixels_FlowFromAllPixels_vec;
vector<Cell> TerrainProcessor::undefined_flow_cells;
vector<pair<int, int>>  TerrainProcessor::possible_outlets;
vector< vector <int>> TerrainProcessor::flats;
bool TerrainProcessor::saveContributingDispersalArea = false;
//Computes topographic index using exponential decay or power decay approaches to TopographicIndex_Log_to_Power ratio of flow accumulation area per contour width to slope
//TopographicIndex_Log_to_Power function converts topographic index (TI) values in map and histogram from exponential decay ln formulation to power decay formulation
//Note: TI transformed by taking exponential of ln(a/tanB), raised to power (1/n), units change from ln([L]) to [L]^(1/n)

//Computes flow accumulation area using D8 (O'Callaghan and Mark, 1984), MFD (Quinn et al., 1991), or D-Infinity (Tarboton, 1997) flow direction approaches
//Horn, B. K. P. (1981). Hill shading and the reflectance map. IEEE Proceedings, 69, 14-47. North = 0, Clockwise to 360 degrees
//O'Callaghan, J. F., & Mark, D. M. (1984). The Extraction of Drainage Networks from Digital Elevation Data. Computer Vision, Graphics, and Image Processing, 28, 323-344. 
//Quinn, P., Beven, K., Chevallier, P., & Planchon, O. (1991). The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling using Digital Terrain Models. Hydrological Processes, 5, 59-79. 
//Tarboton, D. G. (1997). A New Method for the Determination of Flow Directions and Upslope Areas in Grid Digital Elevation Models. Water Resources Research, 33, 309-319. 

//Note: HydroPlusConfig.xml Flag_TI_PerviousOnly=1 excludes impervious area in statistical or spatial mode using TopUrban theory of Valeo and Moin (2000)
//Note: Statistical and Spatial Hydro use different approaches to update FlowAccumulation_Map with pervious area
//Note: Another option is to use pervious area and impervious area draining to pervious, e.g., (1-[DCIA*impervious_fraction])
//Valeo, C., and Moin, S. M. A. (2000). Variable source area modelling in urbanizing watersheds. Journal of Hydrology, 228(1-2), 68-81. 

//Topographic index algorithm taken from TOPMODEL code developed by Keith Beven and used by USGS D.M. Wolock https://pubs.er.usgs.gov/publication/wri934124
//Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124).
//Note: The algorithm requires the topographic index (TI) bins (ia) are sorted from highest (wettest) to lowest TI values
//Note: The algorithm requires when ia = 0, the TI bin TI_Value[0] should be the upper limit and its TI_Area_frac[0] = 0
//Note: The algorithm presumes for each ia, the area of catchment between discrete TI values TI_Value[ia] and TI_Value[ia+1] equals TI_Area_frac[ia]
//Note: The algorithm presumes for each ia, the area of catchment bounding TI_Value[ia] equals (TI_Area_frac[ia] + TI_Area_frac[ia+1])/2 
//Note: The algorithm in C++ has TI values for ia = 0 to ia = [Size_TI_Value_Vector - 1], yet requires use of TI_Area_frac[Size_TI_Value_Vector] = 0 in the above fractional area equation

//TerrainProcessor constructor for Algorithm_SoilKsatDecay = PowerDecay or ExponentialDecay with Inputs pointer
TerrainProcessor::TerrainProcessor(Inputs* input)
{
	//Directory_Input_CLArg is input->SimulationStringParams["InputDirectory"] with ending backslash
	Directory_Input_CLArg = input->SimulationStringParams["InputDirectory"];

	//ImperviousArea_total_frac (frac) is Area_Impervious_frac, defined differently for SpatialTemperatureHydro and StatisticalHydro
	double ImperviousArea_total_frac = input->SimulationLocationParams["Area_Impervious_frac"];
	//stat_PC_frac (fraction) is pervious cover fraction in StatisticalHydro, defined as 1 - ImperviousArea_total_frac
	stat_PC_frac = 1 - ImperviousArea_total_frac;

	//initialize NAC as number of topographic index bins TopographicIndexBins
	NAC = int(input->SimulationNumericalParams["TopographicIndexBins"]);
	
	//TerrainProcessor::init() function in TerrainProcessor.h initializes the following private data members ...
	//... nRows = 0; nCols = 0;  NATB = 0; Length_Pixel_Side_m = 0.0; Area_Catchment_1_frac = 1.0; totalFImpArea = 0.0;
	TerrainProcessor::init();
}

//TopographicIndex_Histogram_Read function reads Topographic Index (TI) as histogram with bins TI values with fractional areas summing to 1
//Note: If z_TI_PowerDecay.csv or z_TI_ExponentialDecay.csv do not exist, or Flag_Recompute_TopographicIndex true, then TI is computed
//Note: TopographicIndex_Histogram_Read for statistically distributed applications, i.e., Hydro, with histogram based TI values
//Note: TopographicIndex_Histogram_Read not called by the spatially distributed applications using the TI map
//Note: Algorithm taken from TOPMODEL code of by Keith Beven and used by USGS D.M. Wolock https://pubs.er.usgs.gov/publication/wri934124
//Note: The algorithm requires the topographic index (TI) bins (ia) are sorted from highest (wettest) to lowest TI values
//Note: The algorithm requires when ia = 0, the TI bin TI_Value[0] should be the upper limit and its TI_Area_frac[0] = 0
//Note: The algorithm for each ia has the area of catchment between discrete TI values TI_Value[ia] and TI_Value[ia+1] equal TI_Area_frac[ia]
//Note: The algorithm for each ia has the area of catchment bounding TI_Value[ia] equal (TI_Area_frac[ia] + TI_Area_frac[ia+1])/2 
//Reference: Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124).
void TerrainProcessor::TopographicIndex_Histogram_Read(Inputs* input)
{
	//file_TI_Histogram defined 
	string file_TI_Histogram;
	//If (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "ExponentialDecay") then file to read is z_TI_ExponentialDecay.csv
	if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "ExponentialDecay") {
		//TI histogram is Exponential Decay, which has much smaller values than Power Decay
		file_TI_Histogram = "z_TI_ExponentialDecay.csv";
	}
	//Else If (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") then file to read is z_TI_PowerDecay.csv
	else if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") {
		//TI histogram is Power Decay, which has much larger values than Exponential Decay
		file_TI_Histogram = "z_TI_PowerDecay.csv";
	}

	//If (!fileExists(Directory_Input_CLArg + file_TI_Histogram) || input->SimulationNumericalParams["Flag_Recompute_TopographicIndex"] == 1)
	//Note: File does not exist or Flag_Recompute_TopographicIndex is true then create file
	//Note: Not re-computing TI saves time, but re-computing allows changes to HydroPlusConfig.xml parameters
	if (!fileExists(Directory_Input_CLArg + file_TI_Histogram) || input->SimulationNumericalParams["Flag_Recompute_TopographicIndex"] == 1) {
		//TerrainProcessor::TopographicIndex_Map_Create(input) called to create TI map
		TerrainProcessor::TopographicIndex_Map_Create(input);
		//TerrainProcessor::TopographicIndex_Map_Create(input) called to create TI histogram
		TerrainProcessor::TopographicIndex_Histogram_Create(input);
	}
	//ifstream readIndex(Directory_Input_CLArg + file_TI_Histogram); open readIndex for processing
	ifstream readIndex(Directory_Input_CLArg + file_TI_Histogram);
	//If (!readIndex) then alert user. 
	if (!readIndex) {
		cout << "Warning: Unable to read this topographic index file within input directory, " << file_TI_Histogram << "." << endl;
	}

	//TerrainProcessor::clearVectorsRelatedToTopographicIndex() to clear vector of values
	TerrainProcessor::clearVectorsRelatedToTopographicIndex();
	//TI_Area_temp and TI_Value_temp assigned arbitrary value of -1111.0
	double TI_Area_temp = -1111.0;
	double TI_Value_temp = -1111.0;
	string TI_file_line;

	//getline reads first TI_file_line of header from readIndex file
	getline(readIndex, TI_file_line);

	//>> reads next TI_file_line from readIndex file, as variable indexSize, and subFarea
	getline(readIndex, TI_file_line);
	//stringstream headerLine(TI_file_line); create headerline from TI_file_line
	stringstream headerLine(TI_file_line);
	//char comma_symbol created to store comma symbol
	char comma_symbol;
	//headerLine >> indexSize >> subFarea; send headerline to indexSize and subFarea
	headerLine >> indexSize >> subFarea;

	//While (getline(readIndex, TI_file_line)); read TI_Area_temp and TI_Value_temp values TI_file_line-by-TI_file_line from CSV
	while (getline(readIndex, TI_file_line)) {
		//stringstream dataLine(TI_file_line)
		stringstream dataLine(TI_file_line);
		//dataLine >> TI_Area_temp >> comma_symbol >> TI_Value_temp; read in values
		dataLine >> TI_Area_temp >> comma_symbol >> TI_Value_temp;
		//TI_Value.push_back(TI_Value_temp); appends value of TI 
		TI_Value.push_back(TI_Value_temp);
		///TI_Area_frac.push_back(TI_Area_temp); appends value of fractional area
		TI_Area_frac.push_back(TI_Area_temp);
	}
	//TI_Area_frac vector extended to hold a 0 value at TI_Value=NAC, after reading last item in histogram NAC-1 given C++ starts vector at 0
	//Note: This final TI_Area_frac value is used in algorithms computing TI_binArea_frac (ACF in Topmodel)
	TI_Area_frac.push_back(0);
	//readIndex.close(); close file
	readIndex.close();

	TI_Value_Average = 0.0;
	double totalArea = 0.0;
	//For (size_t i = 0; i < TI_Value.size(); i++); loop through TI values
	for (size_t i = 0; i < TI_Value.size(); i++) {
		//totalArea += TI_Area_frac[i]; accumulate totalArea
		totalArea += TI_Area_frac[i];
	}

	//Error handling: Guard against zero or NaN totalArea
	if (totalArea <= 0.0 || isnan(totalArea)) {
		cerr << "ERROR: TopographicIndex_Histogram_Read: totalArea <= 0 or NaN." << endl;
		cerr << "  Check TI histogram file for TI_Area entries." << endl;
		cerr << "  TI_Value.size()      = " << TI_Value.size() << endl;
		cerr << "  TI_Area_frac.size()  = " << TI_Area_frac.size() << endl;
		abort();
	}

	//TI_Area_frac[0] = TI_Area_frac[0] / totalArea for first bin
	TI_Area_frac[0] = TI_Area_frac[0] / totalArea;
	//For (size_t i = 1; i < TI_Value.size(); i++); loop through TI values
	for (size_t i = 1; i < TI_Value.size(); i++) {
		//TI_Area_frac[i] = TI_Area_frac[i] / totalArea; compute fraction of area
		//Note: algorithm works for values originally stored as fractions (totalArea = 1) or as area units
		TI_Area_frac[i] = TI_Area_frac[i] / totalArea;
		//TI_Value_Average is the final sum of the (average of neighboring discrete TI values defining bin edges), weighted by the (fractional area within each bin)
		//TI_Value_Average stored as input->TI_Value_Average
		TI_Value_Average = TI_Value_Average + TI_Area_frac[i] * (TI_Value[i] + TI_Value[i - 1]) / 2;
		//Original TOPMODEL Eq TL = TL + AC(J) * (ST(J) + ST(J - 1)) / 2; find algorithm by searching Entire Solution for "CALCULATE AREAL INTEGRAL OF LN(A / TANB)"
	}
}

//TopographicIndex_Map_Read function reads Topographic Index (TI) as map with each cell containing a TI value
//Note: If z_TIorganizer.asc or z_FDorganizer.asc do not exist, or Flag_Recompute_TopographicIndex true, then TI is computed
//Note: TopographicIndex_Map_Read for spatially distributed applications, i.e., CoolAir, with map based TI values
//Note: TopographicIndex_Map_Read not called by the statistically distributed applications using the TI histogram
void TerrainProcessor::TopographicIndex_Map_Read(Inputs* input)
{
	//subFarea is set to 1 for each folder
	subFarea = 1.0;
	//TI_Value_Average for map will be computed in distAveIndex function using map topographic index values
	TI_Value_Average = 0.0;
	indexSize = 0;
	indexAreaSize = 0;

	double temp;
	string space;
	string line, temp2;

	//If z_TIorganizer.asc and z_FDorganizer.asc do not exist or Flag_Recompute_TopographicIndex is true, then compute TI 
	//Note: Not re-computing TI saves time, but re-computing allows changes to HydroPlusConfig.xml parameters
	if (!fileExists(Directory_Input_CLArg + "z_TIorganizer.asc") && !fileExists(Directory_Input_CLArg + "z_FDorganizer.asc") ||
		input->SimulationNumericalParams["Flag_Recompute_TopographicIndex"] == 1) {

		//TopographicIndex_Map_Create function computes topographic index map from elevation data
		TerrainProcessor::TopographicIndex_Map_Create(input);
	}
	//Else If z_TIorganizer.asc and z_FDorganizer.asc exist or Flag_Recompute_TopographicIndex is false, then read DEM
	else {
		//DEM_Read_and_Prepare function called 
		DEM_Read_and_Prepare(input);
	}

	//getElevation1D called with elevation_m_1D
	TerrainProcessor::getElevation1D(elevation_m_1D);
	//getflowAccum1D called with flowAccum1D
	TerrainProcessor::getflowAccum1D(flowAccum1D);

	//ifstream function readTIorganizer directory path and file name to read file, which contains topographic index values for each map pixel
	ifstream readTIorganizer(Directory_Input_CLArg + "z_TIorganizer.asc");

	//TerrainProcessor::clearVectorsRelatedToTopographicIndex() function clears vectors before they are used in the index calculations
	TerrainProcessor::clearVectorsRelatedToTopographicIndex();
	//if readTIorganizer function returns a False value, meaning it did not work, then
	if (!readTIorganizer) {
		cout << "Warning: z_TIorganizer.asc could not be opened." << endl;
	}

	//getline readTIorganizer is reading in the 6 header lines of data from map file
	getline(readTIorganizer, line);
	istringstream(line) >> space >> nCols;
	getline(readTIorganizer, line, '\n');
	istringstream(line) >> space >> nRows;
	getline(readTIorganizer, line, '\n');
	istringstream(line) >> space >> xllcorner_m;
	getline(readTIorganizer, line, '\n');
	istringstream(line) >> space >> yllcorner_m;
	getline(readTIorganizer, line, '\n');
	istringstream(line) >> space >> cellsize_m;
	getline(readTIorganizer, line, '\n');
	istringstream(line) >> space >> NODATA_code;

	//while loop to readTIorganizer all pixel data values, assigned to temp variable
	while (readTIorganizer >> temp) {
		//if pixel value is real data then
		if (temp != NODATA_code) {
			//TI_Value vector has temp value added
			TI_Value.push_back(temp);
		}
		//if pixel is NODATA_code then
		else {
			//TI_Value vector has NODATA_code value added
			TI_Value.push_back(NODATA_code);
		}
	}
	//TI_Area_frac vector stores value of 1 for each pixel; only 1 pixel of area to be used by all spatially distributed map folders
	TI_Area_frac.push_back(1.0);

	//close the readTIorganizer file
	readTIorganizer.close();

	//ifstream function readFDorganizer directory path and file name to read file 
	ifstream readFDorganizer(Directory_Input_CLArg + "z_FDorganizer.asc");

	//if readFDorganizer function returns a False value, meaning it did not work, then
	if (!readFDorganizer) {
		cout << "Warning: z_FDorganizer.asc could not be opened." << endl;
	}

	//getline readFDorganizer is reading in the 6 header lines of data from map file
	getline(readFDorganizer, line);
	istringstream(line) >> space >> nCols;
	getline(readFDorganizer, line, '\n');
	istringstream(line) >> space >> nRows;
	getline(readFDorganizer, line, '\n');
	istringstream(line) >> space >> xllcorner_m;
	getline(readFDorganizer, line, '\n');
	istringstream(line) >> space >> yllcorner_m;
	getline(readFDorganizer, line, '\n');
	istringstream(line) >> space >> cellsize_m;
	getline(readFDorganizer, line, '\n');
	istringstream(line) >> space >> NODATA_code;

	//while loop to readFDorganizer all pixel data values, assigned to temp variable
	while (readFDorganizer >> temp) {
		//if pixel value is real data then
		if (temp != NODATA_code) {
			//FDorganizer vector has temp value added
			FDorganizer_1D.push_back(temp);
		}
		//if pixel is NODATA_code then
		else {
			//FDorganizer vector has NODATA_code value added
			FDorganizer_1D.push_back(NODATA_code);
		}
	}

	//close readFDorganizer file
	readFDorganizer.close();

	//if TI_Value.size is less than input->nCols * input->nRows, z_TIorganizer.asc is different than HydroPlusConfig.xml rows and cols
	if (TI_Value.size() < static_cast<size_t>(Inputs::nCols * Inputs::nRows)) {
		//Warn user, likely end simulation
		cout << "Warning: TI calculations encountered an error." << endl;
	}
	//initialize variables before summing
	int cnt = 0;
	TI_Value_Average = 0;

	//loop through all pixels, i.e., folders of map
	for (size_t i = 0; i < TI_Value.size(); i++) {
		//if pixel is real data, not NODATA_code
		if (TI_Value[i] != NODATA_code) {
			//TI_Value_Average is sum of all pixel topographic index value
			TI_Value_Average = TI_Value_Average + TI_Value[i];
			//cnt counter is advanced + 1
			cnt++;
		}
	}
	//if all pixels were NODATA_code and cnt = 0
	if (cnt == 0) {
		cout << "Warning: Error calculating distributed average index. \n All map pixels are NODATA." << endl;
	}
	//if some pixels were real data and cnt > 0
	else {
		//TI_Value_Average is quotient of summed TI_Value values and summed number of pixels
		//TI_Value_Average stored as input->TI_Value_Average
		TI_Value_Average = TI_Value_Average / cnt;
	}
}

void TerrainProcessor::TopographicIndex_Map_Create(Inputs* input)
{
	//Resizing the vectors to contain NAC number of topographic index bins
	//TI_Values_In_HistogramBin is counter of occurrences of TI values in each topographic index bin
	TI_Values_In_HistogramBin.resize(NAC);
	//TI_Area_frac is area of catchment between two topographic index (TI) values
	TI_Area_frac.resize(NAC);
	//TI_Value is TI, defined as Ln(a/tanB) or (a/tanB)^n 
	TI_Value.resize(NAC);

	//DEM_Read_and_Prepare function called 
	DEM_Read_and_Prepare(input);

	cout << "Computing the slopes across the domain." << endl;
	//DEM_to_Slope_Calc function calculates DEM slope and aspect w/ 3rd order difference algorithm published by Horn (1981)
	//Note: Slope not used for flow direction algorithms (D8, MFD, DInf) which compute their own slopes, e.g., calcFlowDirection_D8
	DEM_to_Slope_Calc();

	cout << "Computing the flow directions across the domain." << endl;
	//DEM_to_FlowDirection_Manager computes DEM flow direction used for flow accumulation, needed in TI calculation
	DEM_to_FlowDirection_Manager();

	cout << "Computing the specific contributing area across the domain." << endl;
	//DEM_to_FlowAccumulation_Manager computes DEM flow accumulation needed in TI calculation
	DEM_to_FlowAccumulation_Manager();

	cout << "Computing the topographic index values across the domain." << endl;
	//TopographicIndex_Map_Compute function called to compute topographic index using different flow routing options
	TopographicIndex_Map_Compute();

	//If Model_Selection is StatisticalHydro or SpatialBufferwGI then enter
	//Note: Consider refactor to allow for TopographicIndex_Map_to_Histogram with Model_Selection = SpatialTemperatureHydro
	if (input->SimulationStringParams["Model_Selection"] == "StatisticalHydro" || input->SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {
		//TopographicIndex_Map_to_Histogram function converts TI map to histogram (log space), with bins holding magnitude of each TI value
		TopographicIndex_Map_to_Histogram();
	}
	//If Algorithm_SoilKsatDecay is PowerDecay then enter
	if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") {
		//TopographicIndex_Log_to_Power function converts topographic index (TI) values from exponential decay to power decay formulation
		//Note: TI transformed by taking exponential of ln(a/tanB), raised to power (1/n), units change from ln([L]) to [L]^(1/n)
		TopographicIndex_Log_to_Power();
	}

	//TopographicCharacteristics_Vector_to_Map function creates contributing area and slope maps to compute Topographic Index
	TopographicCharacteristics_Vector_to_Map();
}

//TopographicIndex_Histogram_Create function writes topographic index (TI) histogram, either as Power Decay or Exponential Decay
void TerrainProcessor::TopographicIndex_Histogram_Create(Inputs* input)
{
	//determine file_TI_Histogram of TI histogram output file, initialize string file_TI_Histogram 
	string file_TI_Histogram;

	//If (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "ExponentialDecay") then write z_TI_ExponentialDecay.csv
	if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "ExponentialDecay") {
		file_TI_Histogram = "z_TI_ExponentialDecay.csv";
	}
	//Else If (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") then write z_TI_PowerDecay.csv
	else if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") {
		file_TI_Histogram = "z_TI_PowerDecay.csv";
	}

	//outfile has Directory_Input_CLArg path and file file_TI_Histogram
	ofstream outfile(Directory_Input_CLArg + file_TI_Histogram);

	//if outfile.good returns a False value, meaning it did not work, then
	if (!outfile.good()) {
		cout << "Warning: Could not write topographic index (TI) histogram file " << file_TI_Histogram << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: Without the TI histogram file, the model cannot create lateral transfers of water." << endl;
		cout << "Correction: Check your DEM data file, dem.asc, perhaps using a Test Case dem.asc to troubleshoot for any DEM errors." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//outfile << "Area_in_TI_Value_Bin_frac,TI_Value_Histogram_Bin" << endl; write header rows
	outfile << "Area_in_TI_Value_Bin_frac,TI_Value_Histogram_Bin" << endl;
	//outfile << NAC << "," << Area_Catchment_1_frac << endl; write NAC which is count of TI bins
	outfile << NAC << "," << Area_Catchment_1_frac << endl;

	//for (int TI_HistogramBin_ID = 0; TI_HistogramBin_ID < NAC; TI_HistogramBin_ID++) 
	//Note: TI bins are sorted from largest TI value to smallest
	for (int TI_HistogramBin_ID = 0; TI_HistogramBin_ID < NAC; TI_HistogramBin_ID++) {
		//outfile << TI_Area_frac[TI_HistogramBin_ID] << "," << TI_Value[TI_HistogramBin_ID] << endl
		outfile << TI_Area_frac[TI_HistogramBin_ID] << "," << TI_Value[TI_HistogramBin_ID] << endl;
	}
	//outfile.close() to close file
	outfile.close();
}

//getAveSoilMoistureDeficit function called for exponential or power decay to initializes folder average soil moisture deficit (AveSMD)
//Note: Theory of updating AveSMD (m) with drainage from vadose zone and runoff from saturated zone based on Wolock (1993).
//Note: Soil volumetric water has 3 zones ranging from wilting point, 0wp to field capacity, 0fc, to saturation, 0sat
//Note: 1) evapotranspiration dynamically contains 0wp to 0sat; 0sat-0fc (gravitational water) releases to lower vadose or saturated zone
//Note: 2) vadose zone dynamically contains 0sat to 0fc; 0sat-0fc drains to saturated zone; it is ground surface to water table
//Note: 3) saturated zone generates subsurface runoff; it can rise into vadose and evapotranspiration zones to fill gravitational voids
//Note: The evapotranspiration and vadose zones are conceptually situated alongside each other, each filling a different volumetric component
//Note: The evapotranspiration zone extends to the limit of evapotranspiration, only processing capillary held water
//Note: The vadose zone is represented by average soil moisture deficit, AveSMD, StorageDeficit_VadoseZone_m, or StorageDeficit_VadoseZone_TI_m
//Note: AveSMD and StorageDeficit_VadoseZone_m are catchment averages, while StorageDeficit_VadoseZone_TI_m varies with topographic index
//Note: For ts=0, AveSMD (m) defined in TerrainProcessor::getAveSoilMoistureDeficit using Q0_mph, T0_m2ph, mM, nN with Exponential or Power Decay
//Note: For ts>0, StorageDeficit_VadoseZone_m & StorageDeficit_VadoseZone_TI_m defined in Flux_LateralGroundWater::calculate as above
//Note: StorageDeficit_FieldCapacity_SoilEvapZone_Max_m (m) = Evapotranspiration_Depth_m * (Soil_FieldCapacity_m3pm3 - Soil_WiltingPoint_m3pm3)
//Note: StorageDeficit_VadoseZone_m (m) = (Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3) * Depth_SoilColumn_Ave_m
//Reference: Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 
double TerrainProcessor::getAveSoilMoistureDeficit(Inputs* input)
{
	//AveSMD (m) is computed below as average soil moisture deficit
	//Note: AveSMD (m) = (Soil_SaturationPoint_m3pm3 - Soil_FieldCapacity_m3pm3) * Depth_SoilColumn_Ave_m
	double AveSMD = 0.0;
	//Q0_mpdt (m/h) is initial subsurface discharge from catchment, used to derive initial AveSMD, from HydroPlusConfig.xml Discharge_Subsurface_Initial_mph
	//Note: BuildDataOrganizer::Check_HydroPlusConfig_Inputs checkes that Discharge_Subsurface_Initial_mph and Q0_mph > 0
	double Q0_mph = input->InputXml["Discharge_Subsurface_Initial_mph"];
	//T0_m2ph (m2/h) is soil transmissivity at saturation or AveSMD = 0 (soil fully saturated), from HydroPlusCofig.xml VadoseZone_Transmissivity_Max_m2ph 
	//Note: BuildDataOrganizer::Check_HydroPlusConfig_Inputs checkes that VadoseZone_Transmissivity_Max_m2ph and T0_m2ph > 0
	double T0_m2ph = input->InputXml["VadoseZone_Transmissivity_Max_m2ph"];
	//mM is Parameter_m_KsatExpDecay, scaling parameter for decay in saturated hydraulic conductivity
	double mM = input->InputXml["Parameter_m_KsatExpDecay"];
	//nN is Parameter_n_KsatPowerDecay, power exponent for decay in saturated hydraulic conductivity
	double nN = input->InputXml["Parameter_n_KsatPowerDecay"];

	//Note: Exponential decay uses mM (m) to describe decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Note: Theory of Exponential function decay explained by Wolock (1993).
	//Note: Beven, Lamb, et al. (1995) give a physical interpretation of the decay parameter m (m) controling effective depth, z (m) of catchment soil profile ...
	//Note: ... This it does interactively in Eq T = T0_m2pdt * exp (-f*z), where f = (Theta_sat - Theta_fc)/m. 
	//Note: ... See Figure 18.1 in Beven, Lamb et al. (1995) for values of m ranging from 0.02 to 0.05 creating effective depths z ranging from 1 to 2 m.
	//Note: Williamson et al. (2009) has Eqs 6, 7, and 8 to derive parameter m (m) from SSURGO data across a thickness where hydraulic conductivity Ksat > 1 um/sec
	//Note: ... where Eq 8 m = (Theta_sat - Theta_fc)/f, Eq 7 f = ln(ConMult)/SoilThickness, and Eq 6 ConMult = Ksat_surface/Ksat_bottom
	//Note: ... The parameter SoilThickness is the distance between the surface and bottom soil layers where Ksat > 1 um/sec

	//Note: Power decay uses mM (m) and nN to describe decay the upper limit for the soil moisture deficit used in Topmodel theory
	//Note: Theory of power function decay explained between Eq 2 and 4 of Wang et al. (2005) and after Eq 5 of Wang et al. (2006)
	//Note: Power function decay theory requires m <= S (AveSMD) based on constraint governing decay with depth of transmissivity, Tz = T0_m2pdt*(1-S/m)^n  
	//Note: nN is coefficient for power function decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
	//Note: Wang et al (2006) after Eq 5 T = T0_m2pdt(1-S/m)^2 explains in power decay the m parameter is upper limit for soil moisture deficit S (m) or AveSMD (m) 
	//Note: ... The Eq 5 a parabolic form of decay in transmissivity, and is generalized in Eq 7, T = T0_m2pdt(1-S/m)^n, Fig 1 n power decay value ranges 0.25 and 16  
	//Note: ... Fig 1 shows how T(z) becomes a smaller fraction of T0_m2pdt (i.e., normalized transmissivity or T(z)/T0_m2pdt) as the storage deficit ratio, S/m, increases and as power function decay rate decrease (N.B. typo in manuscript says increases)

	//References: 
	//Beven, K., Lamb, R., Quinn, P., Romanowics, R., & Freer, J. (1995). TOPMODEL. In V. P. Singh (Ed.), Computer models of watershed hydrology (pp. 627-688). Colorado: Water Resources Publications.
	//Hornberger, G. M., Beven, K. J., Cosby, B. J., & Sappington, D. E. (1985). Shenandoah Watershed Study: Calibration of a Topography-Based, Variable Contributing Area Hydrological Model to a Small Forested Catchment. Water Resources Research, 21(12), 1841-1850.
	//Wang, J., Hassett, J.M., & Endreny, T.A. (2005). An Object Oriented Approach to the Description& Simulation of Watershed Scale Hydrologic Processes.Computersand Geosciences, 31(4), 425 - 435.
	//Wang, J., Endreny, T. A., & Hassett, J. M. (2006). Power function decay of hydraulic conductivity for a TOPMODEL-based infiltration routine Hydrological Processes, 20(18), 3825-3834. 
	//Williamson, T. N., Odom, K. R., Newson, J. K., Downs, A. C., Nelson, H. L., Cinotto, P. J., & Ayers, M. A. (2009). The Water Availability Tool for Environmental Resources (WATER): A Water-budget Modeling Approach for Managing Water-supply Resources in Kentucky–Phase I: Data Processing, Model Development, and Application to Non-karst Areas: US Geological Survey Scientific Investigations Report 2009-5248.
	//Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL (USGS Water-Resources Investigation Report 93-4124). 

	//If (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") thene Power function decay computation of AveSMD 
	//Note: Soil moisture deficits are average for catchment and local for TI bin, computed for both approaches as:
	//Note: PowerDecay (timeStep=0) AveSMD = mM * (1.0 - SZQ) and SZQ = pow(Q0_mph / T0_m2ph, 1 / nN)*(avgIndex);
	//Note: PowerDecay StorageDeficit_VadoseZone_TI_m = AveSMD + mM * pow(Q0_mph/T0_m2ph, 1/nN)*(TI_Value_Average - TI_Value[ia]);
	//Note: Calculations of AveSMD are insensitive to the model time step used if Q0 and T0 are in the same time unit, typically per hour
	if (input->SimulationStringParams["Algorithm_SoilKsatDecay"] == "PowerDecay") {

		//SZQ (m) initialilzed to 0, maximum potential subsurface flow rate, when groundwater table is at the ground surface
		//Note: SZQ (m) is maximum potential subsurface flow rate for watershed pervious area, when groundwater table is at the ground surface
		double SZQ = 0.0;

		//SZQ = pow(Q0_mph / T0_m2ph, 1 / nN) * (input->TI_Value_Average)
		//Note: Eq for SZQ is RHS of Eq for S_bar, AveSMD, given in text below Eq 9 of Wang et. al. (2006)
		//Note: From Wang et al. (2006) S_bar = m[1-(R/T0_m2ph)^(1/n)*Lambda], where Lambda = TI_Value_Average, R = Q0_mph below; and SZQ = (R/T0_m2ph)^(1/n)*Lambda
		//Note: From Wang et al. (2006) S_bar = m[1-(R/T0_m2ph)^(1/n)*Lambda], where Lambda = TI_Value_Average, R = Q0_mph below; and SZQ = (R/T0_m2ph)^(1/n)*Lambda
		//Note: Q0_mph or T0_m2ph must be in same time units and be non-zero values in HydroPlusConfig.xml
		//Note: BuildDataOrganizer::Check_HydroPlusConfig_Inputs checkes that Q0_mph and T0_m2ph > 0; 
		//Note: nN is coefficient for power function decay of saturated hydraulic conductivity with depth, due to decreasing macropores with depth
		SZQ = pow(Q0_mph / T0_m2ph, 1 / nN) * (input->TI_Value_Average);

		//If "Type" != RainBarrel or RoofDisconnect and Q0_mph or T0_m2ph <= 0 in HydroPlusConfig.xml, then model will abort in BuildDataOrganizer to avoid average SMD going to infinity
		//if SZQ greater than 1 set to unity, which is maximum for PowerDecay functional approach
		if (SZQ > 1.0) { SZQ = 1.0; }
		//if SZQ less than 0 set to 0, which is minimum for PowerDecay functional approach
		if (SZQ < 0.0) { SZQ = 0.0; }

		//Compute initial average soil moisture deficit, AveSMD (m) for watershed
		//Eq for AveSMD is below Eq 9 of Wang, J., Endreny, T. A., & Hassett, J. M. (2006). Power function decay of hydraulic conductivity for a TOPMODEL-based infiltration routine Hydrological Processes, 20(18), 3825-3834. 
		//Note: From Wang et al. (2006) S_bar = m[1-(R/T0_m2ph)^(1/n)*Lambda], where SZQ = (R/T0_m2ph)^(1/n)*Lambda, Lambda = TI_Value_Average, R = Q0_mph above
		//Eq for AveSMD is below Eq 5 of Wang, J., Endreny, T. A., & Hassett, J. M. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
		//Note: From Wang et al. (2005) S_bar = m[1-(R/T0_m2ph)^(1/n)*Lambda], where SZQ = (R/T0_m2ph)^(1/n)*Lambda, Lambda = TI_Value_Average, R = Q0_mph above
		AveSMD = mM * (1.0 - SZQ);
	}

	//Else if input->SimulationStringParams["Algorithm_SoilKsatDecay"] != "PowerDecay" then Exponential decay computation of AveSMD
	//Note: Soil moisture deficits are average for catchment and local for TI bin, computed for both approaches as:
	//Note: ExponentialDecay (timeStep=0) AveSMD = -mM * log(Q0_mph / SZQ) and SZQ = T0_m2ph * exp(-avgIndex); or 
	//Note: ExponentialDecay StorageDeficit_VadoseZone_TI_m = AveSMD + mM * (TI_Value_Average - TI_Value[ia]) or 
	else {
		//SZQ (m) initialilzed to 0, maximum potential subsurface flow rate, when groundwater table is at the ground surface
		//Note: SZQ (m) is maximum potential subsurface flow rate for watershed pervious area, when groundwater table is at the ground surface
		double SZQ = 0.0;

		//SZQ = T0_m2ph * exp(-TI_Value_Average); Units of T0_m2ph will balance with Q0_mph when computing AveSMD
		//Note: Eq 41 from Wolock (1993) gives T0_m2ph = product of z (m) and K (m/dt), then adjusted upwards by 2 orders of magnitude
		//Note: Eq for SZQ similar to Eq 31b for Q0_mph in Famiglietti, J. S., & Wood, E. F. (1994). Multiscale Modeling of Spatially Variable Water and Energy Balance Processes. Water Resources Research, 30(11), 3061-3078. 
		//Note: Eq similar to RHS of Eq for S_bar, AveSMD, given in text below Eq 3 of Wang, J., Endreny, T. A., & Hassett, J. M. (2005). Flexible Modeling Package for Topographically Based Watershed Hydrology. Journal of Hydrology, 314(1-4), 78-91. 
		//Note: From Wang et al. (2005) S_bar = -m*ln(R/T0_m2ph) -m*Lambda, where Lambda = TI_Value_Average, R = Q0_mph below 
		//Note: From Wang et al. (2005) given AveSMD Eq below on RHS takes log of SZQ, SZQ becomes T0_m2ph * exp(-Lambda)
		//Note: Eq 32 from Wolock (1993) gives SZQ (m) with per dt implicit, represented as SZQ or q_submax (mm/day), and uses T0_m2ph (m^2/dt) was represented as TRANS or T0_m2ph (mm^2/day), and avgIndex was represented as TOPMEAN or lambda ln(meters)^2
		//Note: Unit balance: exp(-avgIndex) has units of 1/m, given Topographic Index = ln (a / tanB) and a is units of m. 
		SZQ = T0_m2ph * exp(-input->TI_Value_Average);
		//TOPMOEL Code:  SZQ = EXP(T0DT - TL), where TL is average topographic index, and T0DT is natural log of transmissivity, ln(T0_m2ph)

		//If (SZQ > 0 and Q0_mph > 0) then enter ...
		//Note: SZQ <= 0 occurs when T0_m2ph was not defined in HydroPlusConfig.xml with element VadoseZone_Transmissivity_Max_m2ph
		//Note: Q0_mph <= 0 occurs when Q0_mph was not defined in HydroPlusConfig.xml with element Discharge_Subsurface_Initial_mph
		//Note: TI computed and exp(-input->TI_Value_Average) can never yield a negative number.
		if (SZQ > 0 and Q0_mph > 0) {
			//AveSMD (m) is initial average soil moisture deficit for watershed pervious area
			//Note: Eq 4 in Hornberger et al. (1985) Qb = SZQ * EXP(-S/SZM)
			//Note; Eq given as SBAR=-SZM*ALOG(Q0_mph/SZQ) in TMOD9502.FOR TOPMODEL DEMONSTRATION PROGRAM VERSION 95.02, By Keith Beven or ...
			//Note: Eq given as misc.S_mean[0] = - params.m * log(misc.qs0 / misc.qss) in parameter.init.c ...
			//Note: ... https://cran.r-project.org/web/packages/topmodel/index.html
			//Note: Eq similar to AveSMD below Eq 3 in Wang et al. (2005) S_bar = -m*ln(R/T0_m2ph) -m*Lambda, = -m[log(Q0_mph/T0_m2ph) / exp(-Lambda)] where Lambda = TI_Value_Average, R = Q0_mph below
			//Note: Similar to Beven and Kirby (1979), Eq 3, Q0_mph = SZQ exp (AveSMD/mM) or Qb = Q0 exp (S/m), where Q0 is flow when S = 0 ...
			//Note: ... where S = storage deficit and Q0 = SZQ or Q at saturated S
			AveSMD = -mM * log(Q0_mph / SZQ);
		}
		
		//Else If SZQ or Q0_mph <= 0, then allowed in BuildDataOrganizer::Check_HydroPlusConfig_Inputs for certain Model_Selection and Type
		//Note: Model_Selection settings with Hydro using Type GI of RainBarrel or RoofDisconnect permits AveSMD = 0
		//Note: Model_Selection SpatialBuffer permits AveSMD = 0
		else {
			//AveSMD (m) set to 0
			AveSMD = 0;
			//If Model_Selection is SpatialTemperatureModel or StatisticalHydro then alert user
			if (input->SimulationStringParams["Model_Selection"] == "SpatialTemperatureModel" || input->SimulationStringParams["Model_Selection"] == "StatisticalHydro") {
				cout << "Note: Average soil moisture deficit initialized to 0 for model type " << input->SimulationStringParams["Model_Selection"] << "." << endl;
				cout << "Note: To initialize to another value provide in the HydroPlusConfig.xmls the following elements:" << endl;
				cout << "Discharge_Subsurface_Initial_mph, VadoseZone_Transmissivity_Max_m2ph, ..." << endl;
				cout << "... Parameter_m_KsatExpDecay, Parameter_n_KsatPowerDecay." << endl;
			}
		}
	}
	//AveSMD (m) returned as double to call within Inputs::ProcessTopographyData function
	return AveSMD;
}

//DEM_Read_Header function check header
//Note: Header information is needed prior to reading other map inputs, but the DEM data cannot be processed until other inputs are read
void TerrainProcessor::DEM_Read_Header(Inputs* input)
{
	//file name is dem.asc, format recognized by QGIS and ArcGIS as input type for given header
	string location = Directory_Input_CLArg + "dem.asc";
	bool datExists = fileExists(location);
	//MAX is capacity of buffer vector, which contains characters along a row of DEM file
	const int MAX = 80;
	//initialize buffer vector, which contains characters along a row of DEM file, separated by commas
	char buffer[MAX];

	//Note: Consider refactor to use one set of fileExists functionality, such as infile.good from Inputs
	//if file does not exist
	if (!fileExists(location)) {
		cout << "Warning: The input file dem.asc or dem.asc is missing." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: Without DEM data the model cannot simulate the water balance." << endl;
		cout << "Correction: Provide a DEM file in ArcGIS ASCII form with 6 header rows, named dem.asc or dem.asc." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}
	//ifstream function reads infile from location
	ifstream infile(location);

	//Note: Consider refactor to use one set of fileExists functionality, such as infile.good from Inputs
	//Tokenizer function, may check every line of data is followed by a carriage return without spaces
	if (!infile.good()) {
		cout << "Warning: The input file dem.asc or dem.asc is missing." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: Without DEM data the model cannot simulate the water balance." << endl;
		cout << "Correction: Provide a DEM file in ArcGIS ASCII form with 6 header rows, named dem.asc or dem.asc." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//Console message with program status
	cout << "Checking map header information ..." << endl;

	//buffer vector takes on values of infile, up to first 80 columns set by MAX, separated by commas
	infile >> buffer;
	//Ascii_Header_str used to hold map header data read from infile
	Ascii_Header_str = buffer;

	//If (Inputs::Convert_String_to_Lower_Case(Ascii_Header_str) == "ncols"), then read all infile header values
	if (Inputs::Convert_String_to_Lower_Case(Ascii_Header_str) == "ncols") {
		//infile defines nRows, columns, lower left xllcorner_m, yllcorner_m, cell size, nodata value
		infile >> Inputs::nCols >> buffer >> Inputs::nRows >> buffer >> Inputs::xllcorner_m >> buffer >> Inputs::yllcorner_m >> buffer >> Inputs::Length_Pixel_Side_m >> buffer >> Inputs::NODATA_code;
	}
	//Else alert user that map is not properly formatted
	else {
		cout << "Warning: Input map dem.asc does not have nCols as 1st item in header." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model is designed for ESRI formatted ASCII maps." << endl;
		cout << "Correction: Provide maps with ESRI formatted ASCII headers." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}
	//close DEM file
	infile.close();
}

//DEM_Read_and_Prepare function reads in the DEMDATA
void TerrainProcessor::DEM_Read_and_Prepare(Inputs* input)
{
	//tempElv temporarily stores the read in elevation data (m)
	double tempElv;  
	//initialize size
	int size = 0;
	//MAX is capacity of buffer vector, which contains characters along a row of DEM file
	const int MAX = 80;
	//initialize buffer vector, which contains characters along a row of DEM file, separated by commas
	char buffer[MAX];
	//elevation_m_2D.clear and Elevation_DEM_m.clear to initialize vector
	elevation_m_2D.clear();
	input->Elevation_DEM_m.clear();

	//file name is dem.asc, format recognized by QGIS and ArcGIS as input type for given header
	string location = Directory_Input_CLArg + "dem.asc";
	bool datExists = fileExists(location);

	//Note: Consider refactor to use one set of fileExists functionality, such as infile.good from Inputs
	//if file does not exist
	if (!fileExists(location))	{
		cout << "Warning: The input file dem.asc or dem.asc is missing." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: Without DEM data the model cannot simulate the water balance." << endl;
		cout << "Correction: Provide a DEM file in ArcGIS ASCII form with 6 header rows, named dem.asc or dem.asc." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}
	//ifstream function reads infile from location
	ifstream infile(location);    

	//Note: Consider refactor to use one set of fileExists functionality, such as infile.good from Inputs
	//Tokenizer function, may check every line of data is followed by a carriage return without spaces
	if (!infile.good())	{
		cout << "Warning: The input file dem.asc or dem.asc is missing." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: Without DEM data the model cannot simulate the water balance." << endl;
		cout << "Correction: Provide a DEM file in ArcGIS ASCII form with 6 header rows, named dem.asc or dem.asc." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//Console message with program status
	cout << "Processing terrain elevation data ..." << endl;

	//buffer vector takes on values of infile, up to first 80 columns set by MAX, separated by commas
	infile >> buffer;
	//Ascii_Header_str used to hold map header data read from infile
	Ascii_Header_str = buffer;

	//If (Inputs::Convert_String_to_Lower_Case(Ascii_Header_str) == "ncols"), then read all infile header values
	if (Inputs::Convert_String_to_Lower_Case(Ascii_Header_str) == "ncols") {
		//infile defines nRows, columns, lower left xllcorner_m, yllcorner_m, cell size, nodata value
		infile >> nCols >> buffer >> nRows >> buffer >> xllcorner_m >> buffer >> yllcorner_m >> buffer >> Length_Pixel_Side_m >> buffer >> NODATA_code;

		//Define for Inputs::nCols, Inputs::nRows, Inputs::Length_Pixel_Side_m, Inputs::NODATA_code, Inputs::xllcorner_m, Inputs::yllcorner_m
		Inputs::nCols = nCols;
		Inputs::nRows = nRows;
		Inputs::xllcorner_m = xllcorner_m;
		Inputs::yllcorner_m = yllcorner_m;
		Inputs::Length_Pixel_Side_m = Length_Pixel_Side_m;
		Inputs::NODATA_code = NODATA_code;
	}
	//Else alert user that map is not properly formatted
	else {
		cout << "Warning: Input map dem.asc does not have nCols as 1st item in header." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: The HydroPlus model is designed for ESRI formatted ASCII maps." << endl;
		cout << "Correction: Provide maps with ESRI formatted ASCII headers." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//initialize rows in files that will hold data extracted from DEM file
	elevation_m_2D.reserve(nRows);
	Slope_Ground_Horn_rad_Map.reserve(nRows);
	Slope_Ground_rise_run_mpm_Map.reserve(nRows);
	Aspect_Ground_N_0_rad_Map.reserve(nRows);
	FlowDirection_Map.reserve(nRows);
	FlowAccumulation_Map.reserve(nRows);
	NeighborAnalysis_Map.reserve(nRows);
	InFlowDirection_Map.reserve(nRows);
	//TopographicIndex_Map is log form, not power form of TI
	TopographicIndex_Map.reserve(nRows);
	input->Elevation_DEM_m.reserve(nRows * nCols);

	//initialize cols in files that will hold data extracted from DEM file
	for (int row = 0; row < nRows; row++)	{
		//if (Inputs::SimulationNumericalParams["Flag_DEM_RemovePitsFlats"] == 1), only allocate memory space 
		if (Inputs::SimulationNumericalParams["Flag_DEM_RemovePitsFlats"] == 1) {
			//original elevation that read from input
			elevation_working_m_2D.push_back(*new vector<double>(nCols));
			//elevation after filling depression
			elevation_depression_filled_m_2D.push_back(*new vector<double>(nCols));
			//elevation after filling flats
			elevation_flats_filled_m_2D.push_back(*new vector<double>(nCols));
		}
		
		//elevation[row][col] is the final elevation_m_2D map we will use in the following calculations
		elevation_m_2D.push_back(*new vector<double>(nCols));

		Slope_Ground_Horn_rad_Map.push_back(*new vector<double>(nCols));
		Slope_Ground_rise_run_mpm_Map.push_back(*new vector<double>(nCols));
		Aspect_Ground_N_0_rad_Map.push_back(*new vector<double>(nCols));
		FlowDirection_Map.push_back(*new vector<double>(nCols));
		FlowAccumulation_Map.push_back(*new vector<double>(nCols));
		NeighborAnalysis_Map.push_back(*new vector<double>(nCols));
		InFlowDirection_Map.push_back(*new vector<double>(nCols));
		TopographicIndex_Map.push_back(*new vector<double>(nCols));
	}

	double temp = 0.0;
	
	//loop through rows in DEM file along row counter
	for (int row = 0; row < nRows; row++) {
		//loop through cols in DEM file along col counter
		for (int col = 0; col < nCols; col++) {
			//Tokenizer function, may test if every line of data is followed by a carriage return without spaces
			if (!infile.good()) {
				infile.clear();
				//break from reading line
				break;
			}

			//read elevation data from DEM file
			infile >> tempElv;

			if (Inputs::SimulationNumericalParams["Flag_DEM_RemovePitsFlats"] == 1) {
				//initialize elevation_working_m_2D map as original dem
				elevation_working_m_2D[row][col] = tempElv;
				//initialize elevation_depression_filled_m_2D map as original dem
				elevation_depression_filled_m_2D[row][col] = tempElv;
				//initialize elevation_depression_filled_m_2D map as original dem
				elevation_flats_filled_m_2D[row][col] = tempElv;
			}
			//elevation_m_2D (m) 2D vector takes tempElv elevation as original dem; this will be modified by any pit filling
			elevation_m_2D[row][col] = tempElv;
			//Elevation_DEM_m (m) 1D vector takes tempElv elevation as original dem; this will not be modified by any pit filling
			input->Elevation_DEM_m.push_back(tempElv);

			//Initializing new organizers using NODATA_code, and real values later assigned
			Slope_Ground_Horn_rad_Map[row][col] = NODATA_code;
			Slope_Ground_rise_run_mpm_Map[row][col] = NODATA_code;
			Aspect_Ground_N_0_rad_Map[row][col] = NODATA_code;
			FlowDirection_Map[row][col] = NODATA_code;
			NeighborAnalysis_Map[row][col] = NODATA_code;
			InFlowDirection_Map[row][col] = NODATA_code;
			TopographicIndex_Map[row][col] = NODATA_code;
			//FlowAccumulation_Map not set all cells to NODATA_code, and previously set to 1 if dem.asc has data
			//the function flowAccumorganizer_initialize() will be put into FlowAccum() to initialize before accumulating
		}
	}
	//close DEM file
	infile.close();

	// ----------------------------------------------------------------------------------
	// PIT AND FLAT FILLING SECTION
	// ----------------------------------------------------------------------------------
	//This block handles hydrologic conditioning of the DEM by removing pits and resolving flats.
	// It ensures that water flow is continuous and not trapped due to depressions or flat surfaces. It is only triggered when Flag_DEM_RemovePitsFlats = 1 in the config file.

	//Set initial flags and parameters for pit and flat filling
	bool flag_fillDepressions = false;
	bool flag_fillFlats = false;
	int counter = 1;

	//Elevation increment used to force slope in flat areas (in meters)
	//fill_elevation_flats_increment_m (m) defined as 0.001 m 
	//Note: fill_elevation_flats_increment_m (m) should be small (<0.01m) to make the slope enforcement efficient; if > 0.01m then it creates new pits
	double fill_elevation_flats_increment_m = 0.001;
	//maximum allowed iterations to fill depressions and resolve flats set to 5000
	int iteration_max = 5000;

	//turn the flags on if the SimulationStringParams["Flag_DEM_RemovePitsFlats"] == "1"
	if (Inputs::SimulationNumericalParams["Flag_DEM_RemovePitsFlats"] == 1) {
		//identify initial pit cells in DEM
		Find_pit(elevation_working_m_2D);
		//Command Prompt notification
		cout << "Terrain elevation data has " << undefined_flow_cells.size() << " pit cells. Flag_DEM_RemovePitsFlats=1 triggers pit search." << endl;
		
		//Only enable filling process if pits are actually found
		if (!undefined_flow_cells.empty()) {
			flag_fillDepressions = true;
			flag_fillFlats = true;
		}
	}

	//If pit filling is enabled, begin iterative filling loop
	if (flag_fillDepressions) {
		//While loop if undefined_flow_cells is not empty, and pits or flats remain 
		while (!undefined_flow_cells.empty()) {
			//track number of pits before this iteration, pits_count_prior_ts defined as undefined_flow_cells.size() at start of loop
			int pits_count_prior_ts = undefined_flow_cells.size();
			
			// Copy the current elevation state into the temporary working maps
			//Note elevation_depression_filled_m_2D and elevation_flats_filled_m_2D are local temp variable
			//Note elevation_working_m_2D is the latest working copy of the filled elevation, is actively updated during iterations, local temp variable
			for (int row = 0; row < nRows; row++) {
				for (int col = 0; col < nCols; col++) {
					elevation_depression_filled_m_2D[row][col] = elevation_working_m_2D[row][col];
					elevation_flats_filled_m_2D[row][col] = elevation_working_m_2D[row][col];
				}
			}

			//If flag_fillFlats is false then fill depressions but not flats; this is never accessed unless a developer changes default settings above
			//Note: flag_fillFlats is automatically set to true by HydroPlusConfig.xml element Flag_DEM_RemovePitsFlats set to 1
			if (!flag_fillFlats) {
				//elevation_working_m_2D = elevation_depression_filled_m_2D;
				FillDepressions(elevation_working_m_2D);
				//For loop with row through nRows 
				for (int row = 0; row < nRows; row++) {
					//For loop with col through nCols
					for (int col = 0; col < nCols; col++) {
						//Update working elevation map after pit filling
						elevation_working_m_2D[row][col] = elevation_depression_filled_m_2D[row][col];
					}
				}

				//Recompute pits after updating elevations
				Find_pit(elevation_working_m_2D);
			}
			//Else If flag_fillFlats is true then fill depressions and enforce slope in flat areas
			else if (flag_fillFlats) {
				//FillDepressions function called sending elevation_working_m_2D (m) and creating elevation_flats_filled_m_2D (m);
				FillDepressions(elevation_working_m_2D);
				//FixFlats function sends elevation_working_m_2D (m) data and fill_elevation_flats_increment_m (m) value to incrementally update flats
				FixFlats(elevation_working_m_2D, fill_elevation_flats_increment_m);
				//For loop with row through nRows 
				for (int row = 0; row < nRows; row++) {
					//For loop with col through nCols
					for (int col = 0; col < nCols; col++) {
						//Update working elevation map after pit filling
						elevation_working_m_2D[row][col] = elevation_flats_filled_m_2D[row][col];
					}
				}

				//Recompute pits after filling
				Find_pit(elevation_working_m_2D);
			}

			//If counter modulus 100 is equal to zero, then write update
			if (counter % 1000 == 0) {
				//Command Prompt notification of status
				cout << "Adjusting terrain elevation on iteration " << counter << " to remove pits. " << undefined_flow_cells.size() << " pits remaining. " << endl;
			}

			//If counter > iteration_max or undefined_flow_cells.empty() then explain pits were either not filled and iteration_max was exceeded
			if (undefined_flow_cells.empty() || counter > iteration_max) {
				//If counter > iteration_max
				if (counter > iteration_max) {
						cout << "Notice: Rerun the simulation replacing dem.asc with the dem_filled.asc file saved in the output directory." << endl;
						cout << "Explanation: This simulation did not fill all dem.asc elevation pits within the max " << iteration_max << " iterations." << endl;
				}
				break;
			}
			counter++;
		}
		//Command Prompt Notification
		cout << "Flow routing uses the pit-filled elevation data, and writes these data to dem_filled.asc in output folder." << endl;
		cout << "Energy balances and adiabatic temperature adjustment use the original elevation data." << endl;
	}
	
	//loop through rows in DEM file along row counter
	//Note: After filling depressions, save the updated DEM to elevation_m_1D and elevation_m_negative_1D
	for (int row = 0; row < nRows; row++) {
		//loop through cols in DEM file along col counter
		for (int col = 0; col < nCols; col++) {
			//if fill pits or flats
			if (flag_fillDepressions || flag_fillFlats) {
				// Copy the final working DEM (after pit/flat filling) to the main elevation grid
				elevation_m_2D[row][col] = elevation_working_m_2D[row][col];
			}
			
			//old logic  selected different maps depending on pit/flat flags
			/*if (flag_fillDepressions && !flag_fillFlats) {
				elevation_m_2D[row][col] = elevation_depression_filled_m_2D[row][col];
			}
			//if filling depression and filling flats, elevation_m_2D = elevation_flats_filled_m_2D
			else if (flag_fillDepressions && flag_fillFlats) {
				elevation_m_2D[row][col] = elevation_flats_filled_m_2D[row][col];
			}*/

			//elevation_m_1D (m) 1D vector takes values of elevation_m_2D (m) 2D vector
			elevation_m_1D.push_back(elevation_m_2D[row][col]);
			//elevation_m_negative_1D (m) 1D vector takes negative values of elevation_m_2D (m) 2D vector; inverting elevation
			elevation_m_negative_1D.push_back(-elevation_m_2D[row][col]);
		}
	}

	//only write out the filled dem if the Flag_DEM_RemovePitsFlats == 1 and filled depressions or flats
	if (Inputs::SimulationNumericalParams["Flag_DEM_RemovePitsFlats"] == 1 && (flag_fillDepressions == true || flag_fillFlats == true)) {
		//calling the output writer in Buffer to write out the filled dem
		BufferSpatialCalc::BufferModel_IntermediateOutputs_writer("dem_filled");
	}
}

//DEM_to_Slope_Calc function calculates DEM slope and aspect w/ 3rd order difference algorithm published by Horn (1981)
//Note: Nearest points weighted more than diagonal neighbors, center point has no influence; also known as the Sobel operator (Richards, 1986)
//Note: Slope used computing topographic index and for cosine incidence angle with solar radiation
//Note: Slope not used for flow direction algorithms (D8, MFD, DInf) which compute their own slopes, e.g., calcFlowDirection_D8
void TerrainProcessor::DEM_to_Slope_Calc()
{
	//DEM folder layout: 	|6|7|0|
	//						|5|8|1|
	//						|4|3|2|
	double z;
	double organizer[8];
	double fx, fy;
	double Slope_Ground_Horn_temp_rad;
	double progress = 0;

	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int Dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int Dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	//loop through all map rows 
	for (int row = 0; row < nRows; row++) {
		//loop through all map columns
		for (int col = 0; col < nCols; col++) {
			//z is elevation_m_2D value for pixel
			z = getElevation(row, col);

			//if z is real data
			if (z != NODATA_code) {
				//loop through all 8 neighbors
				for (int k = 0; k < 8; k++) {

					//checkfolder function confirms pixel is inside the map
					if (checkfolder(row + Dy[k], col + Dx[k])) {
						//organizer updated for neighbor
						organizer[k] = getElevation(row + Dy[k], col + Dx[k]);
					}
					else {
						organizer[k] = NODATA_code;
					}
					//if organizer is NODATA_code
					if (organizer[k] == NODATA_code) {
						//organizer equal to elevation_m_2D of center pixel
						organizer[k] = z;
					}
				}

				//Calculating the folder's slope based on Horn (1981), with clockwise rotation from NorthEast pixel = organizer[0]
				//New form of Horn Eqn (1981) for fx and fy with sorted addition and subtraction terms
				fx = (organizer[2] + 2 * organizer[1] + organizer[0] - organizer[4] - 2 * organizer[5] - organizer[6]) / (8 * Length_Pixel_Side_m);
				fy = (organizer[6] + 2 * organizer[7] + organizer[0] - organizer[4] - 2 * organizer[3] - organizer[2]) / (8 * Length_Pixel_Side_m);
				//Old form of Horn Eqn (1981) for fx and fy with intermixed addition and subtraction terms
				//fx = (organizer[2] - organizer[4] + 2 * (organizer[1] - organizer[5]) + organizer[0] - organizer[6]) / (8 * Length_Pixel_Side_m);
				//fy = (organizer[6] - organizer[4] + 2 * (organizer[7] - organizer[3]) + organizer[0] - organizer[2]) / (8 * Length_Pixel_Side_m);

				//Slope_Ground_Horn_temp_rad (rad) is slope based on Horn (1981), taking atan of (rise over run) in degrees
				Slope_Ground_Horn_temp_rad = atan(sqrt(fx * fx + fy * fy));

				//fx forced to minimum value when 0 for computation of aspect
				if (fx == 0) {
					fx = 0.0000001;
				}
				//fx forced to minimum value when 0
				if (fy == 0) {
					fy = 0.0000001;
				}
				//Aspect_Ground_N_0_rad_Map is aspect based on Horn (1981) and WhiteBox, w/ PI ~ 3.1415926 and PI/2 ~ 1.5707963 to properly rotate aspect
				//Note: Horn (1981) convention is N=0, 360 clockwise rotation
				Aspect_Ground_N_0_rad_Map[row][col] = M_PI - atan(fy / fx) + M_PI / 2 * fx / abs(fx);

				double Slope_minimum_deg = 0.001;
				//Slope_Ground_Horn_temp_rad converted to deg and checked against minimum value to limit range of topographic index
				if ((Slope_Ground_Horn_temp_rad * Ratio_Degree_to_Radian) <= Slope_minimum_deg) {
					Slope_Ground_Horn_rad_Map[row][col] = Slope_minimum_deg * Ratio_Radian_to_Degree;
				}
				else {
					Slope_Ground_Horn_rad_Map[row][col] = Slope_Ground_Horn_temp_rad;
				}
			}
			else {
				Slope_Ground_Horn_rad_Map[row][col] = NODATA_code;
			}
		}
	}
	/*
	//For loop it from Slope_Ground_Horn_rad_Map.begin to .end 
	for (auto it = Slope_Ground_Horn_rad_Map.begin(); it != Slope_Ground_Horn_rad_Map.end(); it++) {
		for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
			input->SlopeGround_rad.push_back(*it2);
			input->AspectGround_N_0_rad.push_back(*it2);
		}
	}
	*/
}

//DEM_to_FlowDirection_Manager computes DEM flow direction used for flow accumulation, needed in TI calculation
void TerrainProcessor::DEM_to_FlowDirection_Manager()
{
	//calcFlowDirection_DInifity computes the direction DEM cells will flow, using D-Infinity algorithm (Tarboton, 1997)
	//Note: Slope is explicitly computed within the flow direction algorithm to find steepest path facet
	if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "DInfinity") {
		calcFlowDirection_DInf();
	}
	//calcFlowDirection_D8 computes the direction DEM cells will flow, using D8 algorithm (O'Callaghan and Mark, 1984)
	//Note: Slope is explicitly computed within the flow direction algorithm to find steepest single path
	else if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "D8") {
		calcFlowDirection_D8();
	}
	//calcFlowDirection_MFD computes the direction DEM cells will flow, using multiple flow direction algorithm (Quinn et al. 1991)
	//Note: Slope is implicitly computed within the flow direction algorithm to find all downslope paths
	else if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "MFD") {
		calcFlowDirection_MFD();
	}
}

//DEM_to_FlowAccumulation_Manager computes DEM flow accumulation needed in TI calculation
void TerrainProcessor::DEM_to_FlowAccumulation_Manager()
{
	//DX_vec (m) vector used for flow accumulation when it is contributing area per contour length
	DX_vec.resize(nRows * nCols);

	//For loop of i through nRows * nCols, starting at 0
	for (int i = 0; i < nRows * nCols; i++) {
		//If elevation_m_1D is not NODATA_code then enter
		if (elevation_m_1D[i] != NODATA_code) {
			//DX_vec vector stores Length_Pixel_Side_m variable, size of pixel length (m)
			DX_vec[i] = Length_Pixel_Side_m;
		}
		//Else If NODATA_code then 
		else {
			DX_vec[i] = NODATA_code;
		}
	}

	if (Inputs::SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {
		AreaFlowingToPervious_frac_vec.resize(nRows * nCols);
		for (int i = 0; i < nRows * nCols; i++) {
			if (elevation_m_1D[i] != NODATA_code) {
				AreaFlowingToPervious_frac_vec[i] = Inputs::PerCoverdata[i] + (1 - Inputs::InputXml["DirectlyConnectedImperviousArea_frac"]) * Inputs::ImperviousCover_frac[i];
			}
			else {
				AreaFlowingToPervious_frac_vec[i] = NODATA_code;
			}
		}
	}

	//Initialize accumulation organizer to accumulate pervious/impervous area
	//If Model_Selection = SpatialTemperatureHydro then can work with spatial maps of in_weight_raster_frac_vec 
	if (Inputs::SimulationStringParams["Model_Selection"] == "SpatialTemperatureHydro") {
		//If Flag_TI_PerviousOnly=1, then use only pervious area in FlowAccumulation_Map when computing topographic index in spatial mode
		//Note: HydroPlusConfig.xml Flag_TI_PerviousOnly=1 excludes impervious area in statistical or spatial mode using TopUrban theory of Valeo and Moin (2000)
		if (Inputs::SimulationNumericalParams["Flag_TI_PerviousOnly"] == 1) {
			//FlowAccum function (which contains many function calls) is called with two vectors, in_weight_raster_frac_vec and in_weight_raster_vec
			//Note: in_weight_raster_frac_vec set to PerCoverdata, a vector of pervious area (fraction) in each cell, constructed in Inputs.cpp
			//Note: in_weight_raster_vec set to DX_vec, a vector of area per contour length (m) in each cell, constructed above
			FlowAccum(Inputs::PerCoverdata, DX_vec);
		}
		else {
			//FlowAccum function (which contains many function calls) is called with one vector, in_weight_raster_vec
			//Note: in_weight_raster_vec set to DX_vec, a vector of total area (m2) in each cell, constructed above
			FlowAccum(DX_vec);
		}
	}
	//Else if SimulationStringParams != SpatialTemperatureHydro then only 
	else {
		//FlowAccum function (which contains many function calls) is called with one vector, in_weight_raster_vec
		//Note: in_weight_raster_vec set to DX_vec, a vector of total area (m2) in each cell, constructed above
		FlowAccum(DX_vec);
	}

	flowAccum1D.clear();
	//For loop row from FlowAccumulation_Map.begin to FlowAccumulation_Map.end is outer container of 2D vector
	for (auto row = FlowAccumulation_Map.begin(); row != FlowAccumulation_Map.end(); row++) {
		//For loop col from *row.begin to *row.end is inner container of 2D vector
		for (auto col = (*row).begin(); col != (*row).end(); col++) {
			//flowAccum1D is 1D vector initialized to size *col
			flowAccum1D.push_back(*col);
		}
	}
}

//FlowAccum function called for spatial HydroPlus models, with 2 inputs: in_weight_raster_frac_vec and in_weight_raster_vec. 
//Note: Calls 4 functions, 1st initialilize, 2nd pre-accumulate weights, 3rd calculate accumulation, 4th post accumulate weights
void TerrainProcessor::FlowAccum(vector<double>in_weight_raster_frac_vec, vector<double>in_weight_raster_vec) {
	//flowAccumorganizer_initialize initializes FlowAccumOrganizer vector to NODATA_code or 1 
	flowAccumorganizer_initialize();
	//preAccumWeighting function takes product of FlowAccumOrganizer and in_weight_raster_frac_vec; where the weight is impervious or pervious cover fraction, or river cells as 0 or 1
	preAccumWeighting(in_weight_raster_frac_vec);
	//start accumulating
	//calcFlowAccum_Dinf computes the contributing area to DEM cell will flow, using D-Infinity algorithm (Tarboton, 1997)
	if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "DInfinity")
		calcFlowAccum_Dinf();
	//calcFlowAccum_D8 computes the contributing area to DEM cell will flow, using D8 algorithm (O'Callaghan and Mark, 1984)
	else if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "D8")
		calcFlowAccum_D8();
	//calcSCA_DMFD computes the contributing area to DEM cell will flow, using multiple flow direction algorithm (Quinn et al., 1991)
	else if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "MFD")
		calcFlowAccum_MFD();
	
	//Accumulate for Area per contour (Length_Pixel_Side_m) or Area (Length_Pixel_Side_m*Length_Pixel_Side_m) or rise or waters0
	postAccumWeighting(in_weight_raster_vec);
}


//FlowAccum function called for spatial HydroPlus models, with 1 input: in_weight_raster_vec. 
//Note: Calls 4 functions, 1st initialilize, 2nd calculate accumulation, 3rd post accumulate weights
void TerrainProcessor::FlowAccum(vector<double> in_weight_raster_vec) {
	//flowAccumorganizer_initialize initializes FlowAccumOrganizer vector to NODATA_code or 1 
	flowAccumorganizer_initialize();
	//start accumulating
	//calcFlowAccum_Dinf computes the contributing area to DEM cell will flow, using D-Infinity algorithm (Tarboton, 1997)
	if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "DInfinity")
		calcFlowAccum_Dinf();
	//calcFlowAccum_D8 computes the contributing area to DEM cell will flow, using D8 algorithm (O'Callaghan and Mark, 1984)
	else if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "D8") 
		calcFlowAccum_D8();
	//calcSCA_DMFD computes the contributing area to DEM cell will flow, using multiple flow direction algorithm (Quinn et al., 1991)
	else if (Inputs::SimulationStringParams["Algorithm_FlowDirection"] == "MFD")
		calcFlowAccum_MFD();
	//postAccumWeighting function takes accumulates for each pixel the in_weight_raster_vec, which may be area per contour length (Length_Pixel_Side_m) or total area (Length_Pixel_Side_m*Length_Pixel_Side_m), or elevation rise or water ID of 0
	postAccumWeighting(in_weight_raster_vec);
}

//flowAccumorganizer_initialize function initializes FlowAccumulation_Map to NODATA_code or 1
void TerrainProcessor::flowAccumorganizer_initialize() {
	//loop through rows in DEM file along row counter
	for (int row = 0; row < nRows; row++) {
		//loop through cols in DEM file along col counter
		for (int col = 0; col < nCols; col++) {
			//Testing if DEM is actually NODATA_code
			if (elevation_m_2D[row][col] == NODATA_code) {
				//FlowAccumulation_Map is set to NODATA_code
				FlowAccumulation_Map[row][col] = NODATA_code;
			}
			else {
				//FlowAccumulation_Map is initialized to values of 1, unweighted flow accumulation
				FlowAccumulation_Map[row][col] = 1;
			}
		}
	}
}

// Function to set flags for saving contributing/dispersal areas based on conditions
void TerrainProcessor::setContributingDispersalFlags() {
	// Check BufferSpatialWeightingStringParams for "PrintRowCol"
	for (const auto& i : Inputs::BufferSpatialWeightingStringParams) {
		if (i.first.find("PrintRowCol") != string::npos) {
			int writeout_row = atoi(i.second.substr(0, i.second.find(",")).c_str()) - 1;
			int writeout_col = atoi(i.second.substr(i.second.find(",") + 1, i.second.size() - i.second.find(",")).c_str()) - 1;
			if (writeout_row >= 0 && writeout_col >= 0) {
				saveContributingDispersalArea = true;
			}
			break;
		}
	}

	// Check if SpatialBufferGI is a condition to save contributing/dispersal areas
	if (Inputs::SimulationStringParams["Model_Selection"] == "SpatialBufferwGI") {
		saveContributingDispersalArea = true;
	}
	
	if (saveContributingDispersalArea){
		// The safety check is based on the actual number of valid DEM cells(ignores NODATA).
		// Count valid cells once here (O(n))
		int valid_cells = 0;
		for (int i = 0; i < nRows; ++i) {
			for (int j = 0; j < nCols; ++j) {
				if (elevation_m_2D[i][j] != NODATA_code) ++valid_cells;
			}
		}
		// Use a cell-count threshold roughly equivalent to 500x500
		const int max_contrib_pixels = 250000; // 500 * 500

		if (valid_cells > max_contrib_pixels) {
			cout << "Warning: Valid DEM cells (" << valid_cells << ") exceed threshold (" << max_contrib_pixels << "). This may slow the run or crash on constrained systems." << endl;
			cout << "Explanation: Contributing/dispersal pixel write-out skipped to avoid heavy I/O/memory. " << endl;
			//cout << "            The HydroPlusConfig.xml BufferSpatialWeightingStringParams <PrintRowCol#> were ignored." << endl;
			cout << "Suggestion: Simulate a smaller area, or raise this threshold if resources allow. " << endl;
			// When exceeded, we log a warning and turn saveContributingDispersalArea off to prevent potential crashes from heavy memory.
			//saveContributingDispersalArea = false;
		}
	}	
}

//calcFlowDirection_DInf function calculates upslope contributing areas using the D-Infinity algorithm (Tarboton, 1997) 
//Reference: Tarboton, D. G. (1997). A New Method for the Determination of Flow Directions and Upslope Areas in Grid Digital Elevation Models. Water Resources Research, 33, 309-319. 
void TerrainProcessor::calcFlowDirection_DInf()
{
	//Length_Pixel_Diagonal_m is distance between diagonal folder centers
	double Length_Pixel_Diagonal_m = Length_Pixel_Side_m * sqrt(2.0);
	//Elevation_Pixel_m is elevation
	double Elevation_Pixel_m;
	//Slope_Pixel_Max_mpm is maximum slope
	double Slope_Pixel_Max_mpm;
	//flowDir is flow direction
	double flowDir = 0;
	//constant arctan of 1
	double atanof1 = atan(1.0);
	//constants
	double Facet_DInfinity_Multiplier_af;
	double Facet_DInfinity_Constant_ac;
	double Elevation_PixelNeighbor_1_m, ArcTangent_Slope_1_to_Slope_2, Slope_RiseOverRun_1_mpm, Slope_RiseOverRun_2_mpm, Slope_RiseOverRun_mpm, Elevation_PixelNeighbor_2_m;
	int current_row1, current_col1;
	int current_row2, current_col2;

	//DEM facet pixel coordinates for 8 neighbors, as defined in Table 1 of Tarboton (1997) describing D-Infinity
	int acVals[8] = { 0, 1, 1, 2, 2, 3, 3, 4 };
	int afVals[8] = { 1, -1, 1, -1, 1, -1, 1, -1 };

	//DEM 1 pixel coordinates for 8 neighbors, as defined in Table 1 of Tarboton (1997) describing D-Infinity
	int e1Col[8] = { 1, 0, 0, -1, -1, 0, 0, 1 };
	int e1Row[8] = { 0, -1, -1, 0, 0, 1, 1, 0 };

	//DEM 2 pixel coordinates for 8 neighbors, as defined in Table 1 of Tarboton (1997) describing D-Infinity
	int e2Col[8] = { 1, 1, -1, -1, -1, -1, 1, 1 };
	int e2Row[8] = { -1, -1, -1, -1, 1, 1, 1, 1 };

	//loop through all map rows 
	for (int row = 0; row < nRows; row++) {
		//loop through all map columns
		for (int col = 0; col < nCols; col++) {
			//Elevation_Pixel_m is elevation from DEM map 
			Elevation_Pixel_m = elevation_m_2D[row][col];
			//if Elevation_Pixel_m is real data
			if (Elevation_Pixel_m != NODATA_code) {
				//initialize variables for algorithm
				Slope_Pixel_Max_mpm = -9999999999;
				flowDir = 360; //Li: Initialize to -1 instead of 360?
				//loop through all 8 neighbors
				for (int k = 0; k < 8; k++) {
					//Facet_DInfinity_Constant_ac and Facet_DInfinity_Multiplier_af are updated for neighbor
					Facet_DInfinity_Constant_ac = acVals[k];
					Facet_DInfinity_Multiplier_af = afVals[k];
					current_row1 = row + e1Row[k];
					current_col1 = col + e1Col[k];
					current_row2 = row + e2Row[k];
					current_col2 = col + e2Col[k];

					//checkfolder function confirms pixel is inside the map
					Elevation_PixelNeighbor_1_m = checkfolder(current_row1, current_col1) ? elevation_m_2D[current_row1][current_col1] : NODATA_code;
					Elevation_PixelNeighbor_2_m = checkfolder(current_row2, current_col2) ? elevation_m_2D[current_row2][current_col2] : NODATA_code;

					//if Elevation_PixelNeighbor_1_m and Elevation_PixelNeighbor_2_m are real data
					if (Elevation_PixelNeighbor_1_m != NODATA_code && Elevation_PixelNeighbor_2_m != NODATA_code) {

						//if Elevation_Pixel_m is higher than Elevation_PixelNeighbor_1_m and Elevation_PixelNeighbor_2_m
						//Note: calcFlowDirection_DInf requires Epsilon_Tolerance_1E_negative15 for peformance agreement between Linux and Windows builds
						if (Inputs::isGreaterThan(Elevation_Pixel_m, Elevation_PixelNeighbor_1_m,Epsilon_Tolerance_1E_negative15) && Inputs::isGreaterThan(Elevation_Pixel_m,Elevation_PixelNeighbor_2_m, Epsilon_Tolerance_1E_negative15))	{
							//Equations from Tarboton (1997)
							Slope_RiseOverRun_1_mpm = (Elevation_Pixel_m - Elevation_PixelNeighbor_1_m) / Length_Pixel_Side_m;
							if (Slope_RiseOverRun_1_mpm == 0) {
								Slope_RiseOverRun_1_mpm = 0.00001;
							}
							Slope_RiseOverRun_2_mpm = (Elevation_PixelNeighbor_1_m - Elevation_PixelNeighbor_2_m) / Length_Pixel_Side_m;
							ArcTangent_Slope_1_to_Slope_2 = atan(Slope_RiseOverRun_2_mpm / Slope_RiseOverRun_1_mpm);
							Slope_RiseOverRun_mpm = sqrt(Slope_RiseOverRun_1_mpm * Slope_RiseOverRun_1_mpm + Slope_RiseOverRun_2_mpm * Slope_RiseOverRun_2_mpm);
							
							if (Slope_RiseOverRun_1_mpm < 0 && Slope_RiseOverRun_2_mpm < 0) {
								Slope_RiseOverRun_mpm = -1 * Slope_RiseOverRun_mpm;
							}
							if (Slope_RiseOverRun_1_mpm < 0 && Slope_RiseOverRun_2_mpm == 0) {
								Slope_RiseOverRun_mpm = -1 * Slope_RiseOverRun_mpm;
							}
							if (Slope_RiseOverRun_1_mpm == 0 && Slope_RiseOverRun_2_mpm < 0) {
								Slope_RiseOverRun_mpm = -1 * Slope_RiseOverRun_mpm;
							}
							if (Slope_RiseOverRun_1_mpm == 0.00001 && Slope_RiseOverRun_2_mpm < 0) {
								Slope_RiseOverRun_mpm = -1 * Slope_RiseOverRun_mpm;
							}
							if (ArcTangent_Slope_1_to_Slope_2 < 0 || ArcTangent_Slope_1_to_Slope_2 > atanof1)
							{
								if (ArcTangent_Slope_1_to_Slope_2 < 0) {
									ArcTangent_Slope_1_to_Slope_2 = 0;
									Slope_RiseOverRun_mpm = Slope_RiseOverRun_1_mpm;
								}
								else {
									ArcTangent_Slope_1_to_Slope_2 = atanof1;
									Slope_RiseOverRun_mpm = (Elevation_Pixel_m - Elevation_PixelNeighbor_2_m) / Length_Pixel_Diagonal_m;
								}
							}
							if (Inputs::isGreaterThanOrAlmostEqual(Slope_RiseOverRun_mpm, Slope_Pixel_Max_mpm, Epsilon_Tolerance_1E_negative15) && Slope_RiseOverRun_mpm != 0.00001) {
								Slope_Pixel_Max_mpm = Slope_RiseOverRun_mpm;
								flowDir = Facet_DInfinity_Multiplier_af * ArcTangent_Slope_1_to_Slope_2 + Facet_DInfinity_Constant_ac * (M_PI / 2);
							}
						}
						else if (Inputs::isGreaterThan(Elevation_Pixel_m, Elevation_PixelNeighbor_1_m, Epsilon_Tolerance_1E_negative15) || Inputs::isGreaterThan(Elevation_Pixel_m,Elevation_PixelNeighbor_2_m,Epsilon_Tolerance_1E_negative15)) {
							if (Elevation_Pixel_m > Elevation_PixelNeighbor_1_m) {
								ArcTangent_Slope_1_to_Slope_2 = 0;
								Slope_RiseOverRun_mpm = (Elevation_Pixel_m - Elevation_PixelNeighbor_1_m) / Length_Pixel_Side_m;
							}
							else {
								ArcTangent_Slope_1_to_Slope_2 = atanof1;
								Slope_RiseOverRun_mpm = (Elevation_Pixel_m - Elevation_PixelNeighbor_2_m) / Length_Pixel_Diagonal_m;
							}
							if (Inputs::isGreaterThanOrAlmostEqual(Slope_RiseOverRun_mpm, Slope_Pixel_Max_mpm, Epsilon_Tolerance_1E_negative15) && Slope_RiseOverRun_mpm != 0.00001)
							{
								Slope_Pixel_Max_mpm = Slope_RiseOverRun_mpm;
								flowDir = Facet_DInfinity_Multiplier_af * ArcTangent_Slope_1_to_Slope_2 + Facet_DInfinity_Constant_ac * (M_PI / 2);
							}
						}
					}
				}
				//saving the slope to Slope_Ground_rise_run_mpm_Map
				//Dinf assign flow from each pixel to one of its eight neighbors, in the direction with steepest downward slope
				Slope_Ground_rise_run_mpm_Map[row][col] = Slope_Pixel_Max_mpm;

				//if Slope_Pixel_Max_mpm <= 0, then there is no flow to a neighboring cell and no flow direction 
				if (Slope_Pixel_Max_mpm <= 0) {
					//FlowDirection_Map set to -1 to indicate no flow direction
					//Note: No flow direction will require inference to find drainage network
					FlowDirection_Map[row][col] = -1;
				}
				else {
					//re-orient the flow direction angles to match North = 0/360 degrees
					flowDir = floor((flowDir * (180 / M_PI)) * 10) / 10;
					flowDir = 360 - flowDir + 90;
					if (flowDir > 360) {
						flowDir = flowDir - 360;
					}
					FlowDirection_Map[row][col] = flowDir;
				}
			}
			else {
				FlowDirection_Map[row][col] = NODATA_code;
				Slope_Ground_rise_run_mpm_Map[row][col] = NODATA_code;
			}
		}
	}
	//saving the 2D slope organizer to a 1D vector
	for (auto it = Slope_Ground_rise_run_mpm_Map.begin(); it != Slope_Ground_rise_run_mpm_Map.end(); it++) {
		for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
			slopeorganizer_riserun1D.push_back(*it2);
		}
	}
}


//calcFlowDirection_D8 function calculates upslope contributing areas using the D8 (O'Callaghan and Mark, 1984) algorithm from Whitebox GAT
//O'Callaghan, J. F., & Mark, D. M. (1984). The Extraction of Drainage Networks from Digital Elevation Data. Computer Vision, Graphics, and Image Processing, 28, 323-344. 
//Note: HydroPlus WhiteBox flowDir is 1 for NE corner, 2 for E, 4 for SE, 8 for S, 16 for SW, 32 for W, 64 for NW, 128 for N
//Note: ArcGIS ESRI flowDir is 128 for NE corner, 1 for E, 2 for SE, 4 for S, 8 for SW, 16 for W, 32 for NW, 64 for N
void TerrainProcessor::calcFlowDirection_D8()
{
	//e0 and e1 are elevation values
	double e0, e1;
	double slope;
	double Slope_Pixel_Max_mpm;
	//initialize flowDir flow direction
	double flowDir = 0;
	int current_row, current_col;
	double flowDir_ESRI;

	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	//Note: ESRI ArcGIS flow direction vectors arranged clockwise: x=1y=0 in E (right middle), x=1,y=-1 in NE (top right)
	//int dx[8] = { 1, 1, 0, -1, -1, -1, 0, 1 };
	//int dy[8] = { 0, 1, 1, 1, 0, -1, -1, -1 };

	//Length_Pixel_Diagonal_m is distance between diagonal folder centers
	double Length_Pixel_Diagonal_m = Length_Pixel_Side_m * sqrt(2.0);
	//dist vector contains distance to neighbors, along diagonal and cardinal directions
	double dist[8] = { Length_Pixel_Diagonal_m, Length_Pixel_Side_m, Length_Pixel_Diagonal_m, Length_Pixel_Side_m, Length_Pixel_Diagonal_m, Length_Pixel_Side_m, Length_Pixel_Diagonal_m, Length_Pixel_Side_m };

	//loop through all map rows 
	for (int row = 0; row < nRows; row++) {
		//loop through all map columns
		for (int col = 0; col < nCols; col++) {
			//e0 is elevation_m_2D of pixel
			e0 = elevation_m_2D[row][col];

			//if e0 is real data
			if (e0 != NODATA_code) {
				//initialize variables for algorithm
				//Slope_Pixel_Max_mpm initialized to negative value; if no larger slope found when looking at neighbors, no flow direction is the result
				Slope_Pixel_Max_mpm = -9999999999;
				flowDir = -9999999999;
				//loop through all 8 neighbors
				for (int k = 0; k < 8; k++) {
					//current_row and current_col are updated for neighbor
					current_row = row + dy[k];
					current_col = col + dx[k];

					//checkfolder function confirms pixel is inside the map
					if (checkfolder(current_row, current_col)) {
						//e1 updated for neighbor
						e1 = elevation_m_2D[current_row][current_col];
					}
					else {
						e1 = NODATA_code;
					}
					//if e1 is real data then compute slope
					if (e1 != NODATA_code) {
						//slope is simply rise over run
						slope = (e0 - e1) / dist[k];
						//if slope > Slope_Pixel_Max_mpm identify new Slope_Pixel_Max_mpm
						if (Inputs::isGreaterThan(slope, Slope_Pixel_Max_mpm, Epsilon_Tolerance_1E_negative15)) {
							Slope_Pixel_Max_mpm = slope;
							//Note: HydroPlus flowDir is 1 for NE corner, 2 for E, 4 for SE, 8 for S, 16 for SW, 32 for W, 64 for NW, 128 for N
							//Note: ArcGIS ESRI flowDir is 128 for NE corner, 1 for E, 2 for SE, 4 for S, 8 for SW, 16 for W, 32 for NW, 64 for N
							flowDir = pow(2, k);
							//saving the slope to Slope_Ground_rise_run_mpm_Map
							Slope_Ground_rise_run_mpm_Map[row][col] = slope;
							//Note: To compute ESRI flow direction, this code was used, but accuracy is not known
							/*
							flowDir_ESRI = pow(2, k-1);
							if (k == 0) {
								flowDir_ESRI = 128;
							}
							*/
						}
					}
				}
				//if Slope_Pixel_Max_mpm <= 0, then there is no flow to a neighboring cell and no flow direction 
				if (Slope_Pixel_Max_mpm <= 0) {
					//FlowDirection_Map set to -1 to indicate no flow direction
					//Note: No flow direction will require inference to find drainage network
					FlowDirection_Map[row][col] = -1;
				}
				else {
					FlowDirection_Map[row][col] = flowDir;
				}
			}
			else {
				FlowDirection_Map[row][col] = NODATA_code;
				Slope_Ground_rise_run_mpm_Map[row][col] = NODATA_code;
			}
		}
	}

	//saving the 2D slope organizer to a 1D vector
	for (auto it = Slope_Ground_rise_run_mpm_Map.begin(); it != Slope_Ground_rise_run_mpm_Map.end(); it++) {
		for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
			slopeorganizer_riserun1D.push_back(*it2);
		}
	}
}

//calcFlowDirection_MFD function calculates upslope contributing areas using the multiple flow direction algorithm (Quinn et al., 1991) from Whitebox GAT
//Quinn, P., Beven, K., Chevallier, P., & Planchon, O. (1991). The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling using Digital Terrain Models. Hydrological Processes, 5, 59-79. 
//Note: HydroPlus flowDir is 1 for NE corner, 2 for E, 4 for SE, 8 for S, 16 for SW, 32 for W, 64 for NW, 128 for N
//Note: ArcGIS ESRI flowDir is 128 for NE corner, 1 for E, 2 for SE, 4 for S, 8 for SW, 16 for W, 32 for NW, 64 for N
void TerrainProcessor::calcFlowDirection_MFD()
{
	//e0 and e1 are elevation values
	double e0, e1;
	//initialize flowDir flow direction
	double flowDir = 0;
	int current_row, current_col;
	double slope, total_slope;
	int slope_downslope_count;
	double slope_avg;
	//Note: Whitebox (HydroPlus) flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	//Note: ESRI ArcGIS flow direction vectors arranged clockwise:, x=1y=0 in E (right middle), x=1,y=-1 in NE (top right)
	//int dx[8] = { 1, 1, 0, -1, -1, -1, 0, 1 };
	//int dy[8] = { 0, 1, 1, 1, 0, -1, -1, -1 };

	//Length_Pixel_Diagonal_m is distance between diagonal folder centers
	double Length_Pixel_Diagonal_m = Length_Pixel_Side_m * sqrt(2.0); 
	//dist vector contains distance to neighbors, along diagonal and cardinal directions
	double dist[8] = { Length_Pixel_Diagonal_m, Length_Pixel_Side_m, Length_Pixel_Diagonal_m, Length_Pixel_Side_m, Length_Pixel_Diagonal_m, Length_Pixel_Side_m, Length_Pixel_Diagonal_m, Length_Pixel_Side_m };

	for (int row = 0; row < nRows; row++) {
		for (int col = 0; col < nCols; col++) {

			//e0 is pixel elevation_m_2D
			e0 = elevation_m_2D[row][col];
			//if e0 is real data
			if (e0 != NODATA_code) {
				//Note: No flow direction will require inference to find drainage network
				flowDir = 0;
				//loop through all 8 neighbors
				for (int k = 0; k < 8; k++)	{
					//current_row is neighbor row
					current_row = row + dy[k];
					//current_col is neighbor column
					current_col = col + dx[k];

					//checkfolder function confirms pixel is inside the map
					if (checkfolder(row + dy[k], col + dx[k]))
						//e1 updated for neighbor
						e1 = elevation_m_2D[current_row][current_col];
					else
						e1 = NODATA_code;

					// Beginning of the slope calculation loop
					//Note: HydroPlus flowDir is 1 for NE corner, 2 for E, 4 for SE, 8 for S, 16 for SW, 32 for W, 64 for NW, 128 for N
					//Note: ArcGIS ESRI flowDir is 128 for NE corner, 1 for E, 2 for SE, 4 for S, 8 for SW, 16 for W, 32 for NW, 64 for N
					if (e1 != NODATA_code && Inputs::isLessThan(e1, e0 ,Epsilon_Tolerance_1E_negative15))
					{
						//Using the base-2 naming 
						flowDir += pow(2, k);
						//slope is simply rise over run
						slope = (e0 - e1) / dist[k];
						//the count of downslope +1 
						slope_downslope_count++;
						//total_slope for average
						total_slope += slope;
					}
				}
				//saving flowDir to FlowDirection_Map
				FlowDirection_Map[row][col] = flowDir;

				if (slope_downslope_count > 0) {
					// slope_avg is the averge slope of all the downslpre cells
					slope_avg = total_slope / slope_downslope_count;
					//saving the slope to Slope_Ground_rise_run_mpm_Map
					Slope_Ground_rise_run_mpm_Map[row][col] = slope_avg;
				}
				else {
					Slope_Ground_rise_run_mpm_Map[row][col] = 0;
				}
				//reset total_slope and slope_downslope_count to 0
				total_slope = 0;
				slope_downslope_count = 0;
			}
			else {
				FlowDirection_Map[row][col] = NODATA_code;
				Slope_Ground_rise_run_mpm_Map[row][col] = NODATA_code;
			}
		}
	}
	//saving the 2D slope organizer to a 1D vector
	for (auto it = Slope_Ground_rise_run_mpm_Map.begin(); it != Slope_Ground_rise_run_mpm_Map.end(); it++) {
		for (auto it2 = (*it).begin(); it2 != (*it).end(); it2++) {
			slopeorganizer_riserun1D.push_back(*it2);
		}
	}
}


//preAccumWeighting function for spatial HydroPlus models defines FlowAccumulation_Map = in_weight_raster_frac_vec
//Note: weight_vector is impervious or pervious cover fraction, or river cells as 0 or 1
void TerrainProcessor::preAccumWeighting(vector<double> weight_vector) {
	//loop through rows in DEM file along row counter
	for (int row = 0; row < nRows; row++) {
		//loop through cols in DEM file along col counter
		for (int col = 0; col < nCols; col++) {
			//MapPixel_ID is map vector location, i.e., DataDrawer_ID, computed as a function of col_max and pixel row, column ID
			//Note: Conversion Eq: MapPixel_ID = row * col_max + col, with MapPixel_ID starting at 0, row & col starting at 1
			//Note: Conversion Eq: row = [int(MapPixel_ID/col_max)+1]; col = [mod(MapPixel_ID,col_max)+1]
			int col_max = nCols;
			//MapPixel_ID is defined as row * col_max + col, a function that converts row col array index to position in vector
			int MapPixel_ID = row * col_max + col;
			//FlowAccumulation_Map was set to 1 for all data other than NODATA_code
			if (FlowAccumulation_Map[row][col] != NODATA_code && weight_vector[MapPixel_ID] != NODATA_code) {
				//FlowAccumulation_Map[row][col] used to initialize to 1 per valid cell. Since in this case we apply weight_vector which is typically fraction of pervious area,
				//FlowAccumulation_Map is re-initialized as weight_vector
				//Note: Calculation peformed before FlowAccumulation_Map accumulates all drainage area pixels
				FlowAccumulation_Map[row][col] = weight_vector[MapPixel_ID];
			}
		}
	}
}

//postAccumWeighting function for spatial HydroPlus models defines FlowAccumulation_Map = in_weight_raster_frac_vec
//Note: weight_vector is typically cell length or area, as area per contour (m2/m = m) or area (m2)
void TerrainProcessor::postAccumWeighting(vector<double> weight_vector) {
	//loop through rows in DEM file along row counter
	for (int row = 0; row < nRows; row++) {
		//loop through cols in DEM file along col counter
		for (int col = 0; col < nCols; col++) {
			//accumulate Area perContour
			int col_max = nCols;
			//MapPixel_ID is defined as row * col_max + col, a function that converts row col array index to position in vector
			int MapPixel_ID = row * col_max + col;
			//FlowAccumulation_Map was set to 1 for all data other than NODATA_code
			if (FlowAccumulation_Map[row][col] != NODATA_code && weight_vector[MapPixel_ID] != NODATA_code) {
				//FlowAccumulation_Map is updated as product of FlowAccumulation_Map and weight_vector, which is typically Length_Pixel_Side_m, the area per contour length
				//Note: Calculation is peformed after FlowAccumulation_Map accumulated all drainage area pixels
				FlowAccumulation_Map[row][col] = FlowAccumulation_Map[row][col] * weight_vector[MapPixel_ID];
			}
		}
	}
}

//calcFlowAccum_Dinf function calculates specific contributing area (SCA) based on the D-Infinity (Tarboton, 1997)
//Reference: Tarboton, D. G. (1997). A New Method for the Determination of Flow Directions and Upslope Areas in Grid Digital Elevation Models. Water Resources Research, 33, 309-319. 
void TerrainProcessor::calcFlowAccum_Dinf()
{
	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dCol[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dRow[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	double organizerRes = 1;
	double startFD[8] = { 180, 225, 270, 315, 0, 45, 90, 135 };
	double endFD[8] = { 270, 315, 360, 45, 90, 135, 180, 225 };

	double flowDir, flowAccum;
	int row, col, newRow, newCol;
	int counter_int;
	int loopNum = 0;

	//maxDepth is maximum number of pixels in flow path length; 1000 used in Whitebox GAT (John Linsdsay)
	maxDepth = 10000;

	//somethingDone is boolean flag used to enter and exit do-while loop tracking flow path 
	bool somethingDone;

	//For loop row 0 to nRows
	for (row = 0; row < nRows; row++) {
		//For loop col 0 to nCols
		for (col = 0; col < nCols; col++) {
			//flowDir is FlowDirection_Map map value
			flowDir = FlowDirection_Map[row][col];
			//If flowDir is not NODATA_code
			if (flowDir != NODATA_code) {
				//counter_int set to 0
				counter_int = 0;
				//For loop k through all 8 cell neighbors
				for (int k = 0; k < 8; k++)	{
					//newCol is neighbor col
					newCol = col + dCol[k];
					//newRow is neighbor row
					newRow = row + dRow[k];

					//If checkfolder function returns true then neighbor pixel is inside the map
					if (checkfolder(newRow, newCol)) {
						//flowDir updated for neighbor value
						flowDir = FlowDirection_Map[newRow][newCol];
					}
					//Else continue to next neighbor
					else {
						continue;
					}
					//If flowDir is within allowable range of 0 to 360 deg
					if (flowDir >= 0 && flowDir <= 360) {
						//If neighbor is not directly south, which has endFD < startFD
						if (k != 3) {
							//If flowDir > startFD AND < endFD then increase counter
							if (flowDir > startFD[k] && flowDir < endFD[k]) {
								//increment flow direction counter
								counter_int++;
							}
						}
						//Else If neighbor is directly south, which has endFD > startFD, divided across north
						else {
							//If flowDir > startFD OR < endFD then increase counter
							if (flowDir > startFD[k] || flowDir < endFD[k]) {
								//increment flow direction counter 
								counter_int++;
							}
						}
					}
				}
				//assign flow direction counter to maps
				NeighborAnalysis_Map[row][col] = counter_int;
				InFlowDirection_Map[row][col] = counter_int;
			}
			//Else If flowDir is NODATA_code
			else {
				FlowAccumulation_Map[row][col] = NODATA_code;
			}
		}
	}

	//do loop initiates when somethingDone = false, and continues if somethingDone becomes true 
	//Note: to accumulate flow from upslope cells draining to pixel 
	do {
		loopNum++;
		//somethingDone flag to False
		somethingDone = false;
		//For loop row from 0 to nRows 
		for (row = 0; row < nRows; row++) {
			//For loop col from 0 to nCols 
			for (col = 0; col < nCols; col++) {
				//If getElevation is not NODATA_code 
				if (getElevation(row, col) != NODATA_code) {
					//if flow has not been computed for pixel
					if (NeighborAnalysis_Map[row][col] == 0) {   
						//initialize currentDepth to 0 to accumulate pixels
						currentDepth = 0;
						//somethingDone set to True
						somethingDone = true;
						//DInfAccum function called to find flow accumulation area
						FlowAccumulation_DInfinity(row, col);
					}
				}
				//Else If getElevation is NODATA_code 
				else {
					//continue to next pixel
					continue;
				}
			}
		}
	} while (somethingDone);
	// End of Dinf Accum main function
}


//calcFlowAccum_D8 function calculates specific contributing area (SCA) based on the D8 (O'Callaghan and Mark, 1984) 
//Reference: O'Callaghan, J. F., & Mark, D. M. (1984). The Extraction of Drainage Networks from Digital Elevation Data. Computer Vision, Graphics, and Image Processing, 28, 323-344. 
void TerrainProcessor::calcFlowAccum_D8()
{	
	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	//inflowDirVals inflow direction vector starts Northeast (NE=16), rotating clockwise to E=32, SE=64, S=128, SW=1, W=2, NW=4, N=8
	//Note: HydroPlus or Whitebox uses this inflowDirVals structure
	double inflowDirVals[8] = { 16, 32, 64, 128, 1, 2, 4, 8 };

	//Note: HydroPlus or WhiteBox flowDir is 1 for NE corner, 2 for E, 4 for SE, 8 for S, 16 for SW, 32 for W, 64 for NW, 128 for N
	//Whitebox FDR		Whitebox -FDR	Pairs of x,y 
	//64	128	1		4	8	16		-1,-1	0,-1	1,-1
	//32		2		2		32		-1,0			1,0
	//16	8	4		1	128	64		-1,1	0,1		1,1

	//If saveContributingDispersalArea is true initialize the ContributingAreaPixels_FlowIntoAllPixels_vec vector
	if (saveContributingDispersalArea) {
		//ContributingAreaPixels_FlowIntoAllPixels_vec clear and resize to nRows * nCols
		ContributingAreaPixels_FlowIntoAllPixels_vec.clear();
		ContributingAreaPixels_FlowIntoAllPixels_vec.resize(nRows * nCols);
	}
	
	//For loop row 0 to nRows
	for (int row = 0; row < nRows; ++row) {
		//For loop col 0 to nCols
		for (int col = 0; col < nCols; ++col) {
			//If elevation_m_2D is not NODATA_code 
			if (elevation_m_2D[row][col] != NODATA_code) {
				//If saveContributingDispersalArea is true add the pixel itself to ContributingAreaPixels_FlowIntoAllPixels_vec
				if (saveContributingDispersalArea) {
					ContributingAreaPixels_FlowIntoAllPixels_vec[(row * nCols) + col].push_back(make_pair(row, col));
				}
			}
		}
	}

	//Warning: ESRI ArcGIS method not used in HydroPlus, remove or leave commented out
	//Note: ESRI ArcGIS flow direction vectors arranged clockwise: x=1y=0 in E (right middle), x=1,y=-1 in NE (top right)
	//int dx[8] = { 1, 1, 0, -1, -1, -1, 0, 1 };
	//int dy[8] = { 0, 1, 1, 1, 0, -1, -1, -1 };

	//ESRI ArcGIS Pairs of x, y
	//0, -1		0, -1	1, 0
	//-1, -1			1, 1
	//-1, 0		-1, 1	0, 1

	//inflowDirVals inflow direction vector starts Northeast (NE=32), rotating clockwise to E=64, SE=128, S=1, SW=2, W=4, NW=8, N=16
	//Note: ESRI ArcGIS uses this inflowDirVals structure
	//double inflowDirVals[8] = {32, 64, 128, 1, 2, 4, 8, 16 };

	//Note: ArcGIS ESRI flowDir is 128 for NE corner, 1 for E, 2 for SE, 4 for S, 8 for SW, 16 for W, 32 for NW, 64 for N
	//ESRI FDR			ESRI -FDR		Pairs of x,y 
	//32	64	128		2	4	8		-1,-1	0,-1	1,-1
	//16		1		1		16		-1,0			1,0
	//8		4	2		128	64	32		-1,1	0,1		1,1

	double flowDir, flowAccum;
	int row, col, newRow, newCol, temporganizerVal, FlowDirection_ID;
	int z = 0;
	int loopNum = 0;

	//somethingDone is boolean flag used to enter and exit do-while loop tracking flow path 
	bool somethingDone;

	//For loop row from 0 to nRows
	for (int row = 0; row < nRows; row++)	{
		//For loop col from 0 to nCols
		for (int col = 0; col < nCols; col++)	{
			//flowDir is FlowDirection_Map map value
			flowDir = FlowDirection_Map[row][col];
			//If flowDir is not NODATA_code
			if (flowDir != NODATA_code) {
				//z set to 0
				z = 0;
				//For loop k through all 8 cell neighbors
				for (int k = 0; k < 8; k++) {
					//newCol is neighbor col
					newCol = col + dx[k];
					//newRow is neighbor row
					newRow = row + dy[k];

					//If checkfolder function returns true then neighbor pixel is inside the map
					if (checkfolder(newRow, newCol)) {
						//flowDir updated for neighbor value
						flowDir = FlowDirection_Map[newRow][newCol];
						//if flowDir with neighbor value equals inflow direction, then the neighbor flow into this cell
						if (flowDir == inflowDirVals[k]) {
							//z is advanced, counting how many neighboring cells flow into cell
							z++;
						}
					}
				}
				//The counts are stored in the NeighborAnalysis_Map and InFlowDirection_Map maps.
				NeighborAnalysis_Map[row][col] = z;
				InFlowDirection_Map[row][col] = z;
			}
			//If flowDir is NODATA_code
			else {
				NeighborAnalysis_Map[row][col] = NODATA_code;
				InFlowDirection_Map[row][col] = NODATA_code;
				FlowAccumulation_Map[row][col] = NODATA_code;
			}
		}
	}

	//For loop row from 0 to nRows 
	for (row = 0; row < nRows; row++) {
		//For loop col from 0 to nCols 
		for (col = 0; col < nCols; col++) {
			//temporganizerVal is int of NeighborAnalysis_Map[row][col]
			temporganizerVal = int(NeighborAnalysis_Map[row][col]);
			//If the cell has a flow direction counter of 0, indicating that flow accumulation needs to be calculated start from here
			//Note: For cells with temporganizerVal > 0, they have already contributed flow to downstream cells.
			if (temporganizerVal == 0) {
				//NeighborAnalysis_Map for rwo and col set to -1
				NeighborAnalysis_Map[row][col] = -1.0;
				//somethingDone initialized to false, which does not prevent do loop from running once
				somethingDone = false;
				//newCol and newRow initialized to cell, before searching neighbors
				newCol = col;
				newRow = row;
				
				//do loop initiates when somethingDone = false, and continues if somethingDone becomes true 
				//Note: to accumulate flow from upslope cells draining to pixel 
				do {
					//The flow accumulation value and flow direction at the current cell are retrieved.
					flowAccum = FlowAccumulation_Map[newRow][newCol];
					flowDir = FlowDirection_Map[newRow][newCol];

					//ContributingAreaPixels_temp initalized 
					vector<pair<int, int>> ContributingAreaPixels_temp;
					//If saveContributingDispersalArea is true
					if (saveContributingDispersalArea) {
						//ContributingAreaPixels_temp inserted with ContributingAreaPixels_FlowIntoAllPixels_vec.begin to .end 
						//Note: contributing cells from the current cell are brought into the downstream cell
						ContributingAreaPixels_temp.insert(
							ContributingAreaPixels_temp.end(),
							ContributingAreaPixels_FlowIntoAllPixels_vec[(newRow * nCols) + newCol].begin(),
							ContributingAreaPixels_FlowIntoAllPixels_vec[(newRow * nCols) + newCol].end());
					}
					//If flowDir >= 0 then use flow direction to identify the neighboring cell that receives flow from the current cell.
					if (flowDir >= 0) {
						//FlowDirection_ID is int of log2 of flowDir, which is value of 0 to 8, for accessing dx and dy vector
						FlowDirection_ID = int(log2(flowDir));
						//prevCol and prevRow saved prior to going to neighbor
						int prevCol = newCol;
						int prevRow = newRow;
						//newCol and newRow take on dx and dy vector neighbor values
						newCol += dx[FlowDirection_ID];
						newRow += dy[FlowDirection_ID];

						//FlowAccumulation_Map map updated with inflow of flowAccum from neighboring cell
						FlowAccumulation_Map[newRow][newCol] = FlowAccumulation_Map[newRow][newCol] + flowAccum;
										
						//If saveContributingDispersalArea is true
						if (saveContributingDispersalArea) {
							//ContributingAreaPixels_temp inserted with ContributingAreaPixels_FlowIntoAllPixels_vec.begin to .end 
							//Note: contributing cells from the current cell are brought into the downstream cell
							ContributingAreaPixels_FlowIntoAllPixels_vec[(newRow * nCols) + newCol].insert(
								ContributingAreaPixels_FlowIntoAllPixels_vec[(newRow * nCols) + newCol].end(),
								ContributingAreaPixels_temp.begin(),
								ContributingAreaPixels_temp.end());
						}

						//NeighborAnalysis_Map for newRow and newCol neighbor decreased by 1 
						NeighborAnalysis_Map[newRow][newCol] = NeighborAnalysis_Map[newRow][newCol] - 1; 
						//If NeighborAnalysis_Map for neighbor equals 0 then it has no more contributing neighbors, it is marked as -1, 
						//Note: it means this downstream cell has received all possible inflow from its upstream neighbors.
						//Note: this cell is now ready to contribute its accumulated flow to its own downstream neighbor.
						if (NeighborAnalysis_Map[newRow][newCol] == 0) {
							NeighborAnalysis_Map[newRow][newCol] = -1;
							//set somethingDone to true to continue the loop to process this new cell
							//Note: this will become an "upstream" cell for further downstream cells.
							somethingDone = true;
						}
						else {
							//Otherwise, set somethingDone to false to exit the loop.
							somethingDone = false;
						}
					}
					else {
						somethingDone = false;
					}
				} while (somethingDone);
			}
		}
	}
	//If saveContributingDispersalArea is true then 
	if (saveContributingDispersalArea) {
		//trackDispersalAreas_D8 called to track dispersal areas
		trackDispersalAreas_D8();
	}
}

void TerrainProcessor::trackDispersalAreas_D8() {
	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	DispersalAreaPixels_FlowFromAllPixels_vec.clear();
	DispersalAreaPixels_FlowFromAllPixels_vec.resize(nRows * nCols);

	// Loop through all cells to update their dispersal areas
	for (int row = 0; row < nRows; row++) {
		for (int col = 0; col < nCols; col++) {
			int index = (row * nCols) + col;
			// For each contributing cell of the current cell, add the current cell to its dispersal area
			for (const auto& contributingCell : ContributingAreaPixels_FlowIntoAllPixels_vec[index]) {
				int contributingIndex = (contributingCell.first * nCols) + contributingCell.second;
				// Add the current cell to the dispersal area of the contributing cell
				DispersalAreaPixels_FlowFromAllPixels_vec[contributingIndex].push_back(make_pair(row, col));
				// Also add all current cell's dispersal areas to the contributing cell's dispersal area
				/*DispersalAreaPixels_FlowFromAllPixels_vec[contributingIndex].insert(
					DispersalAreaPixels_FlowFromAllPixels_vec[contributingIndex].end(),
					DispersalAreaPixels_FlowFromAllPixels_vec[index].begin(),
					DispersalAreaPixels_FlowFromAllPixels_vec[index].end()
				);*/
			}
		}
	}
}

//calcFlowAccum_MFD function calculates specific contributing area (SCA) based on the multiple flow direction (Quinn et al., 1991)
//Reference: Quinn, P., Beven, K., Chevallier, P., & Planchon, O. (1991). The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling using Digital Terrain Models. Hydrological Processes, 5, 59-79. 
void TerrainProcessor::calcFlowAccum_MFD()
{
	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	//inflowDirVals inflow direction vector starts Northeast (NE=16), rotating clockwise to E=32, SE=64, S=128, SW=1, W=2, NW=4, N=8
	//Note: HydroPlus or Whitebox uses this inflowDirVals structure
	double inflowDirVals[8] = { 16, 32, 64, 128, 1, 2, 4, 8 };

	//Note: HydroPlus or WhiteBox flowDir is 1 for NE corner, 2 for E, 4 for SE, 8 for S, 16 for SW, 32 for W, 64 for NW, 128 for N
	//Whitebox FDR		Whitebox -FDR	Pairs of x,y 
	//64	128	1		4	8	16		-1,-1	0,-1	1,-1
	//32		2		2		32		-1,0			1,0
	//16	8	4		1	128	64		-1,1	0,1		1,1

	//Warning: ESRI ArcGIS method not used in HydroPlus, remove or leave commented out
	//Note: ESRI ArcGIS flow direction vectors are arranged with clockwise indexing
	//int dx[8] = { 1, 1, 0, -1, -1, -1, 0, 1 };
	//int dy[8] = { 0, 1, 1, 1, 0, -1, -1, -1 };

	//ESRI ArcGIS Pairs of x, y
	//0, -1		0, -1	1, 0
	//-1, -1			1, 1
	//-1, 0		-1, 1	0, 1

	//inflowDirVals inflow direction vector starts Northeast (NE=32), rotating clockwise to E=64, SE=128, S=1, SW=2, W=4, NW=8, N=16
	//Note: ESRI ArcGIS uses this inflowDirVals structure
	//double inflowDirVals[8] = {32, 64, 128, 1, 2, 4, 8, 16 };

	//Note: ArcGIS ESRI flowDir is 128 for NE corner, 1 for E, 2 for SE, 4 for S, 8 for SW, 16 for W, 32 for NW, 64 for N
	//ESRI FDR			ESRI -FDR		Pairs of x,y 
	//32	64	128		2	4	8		-1,-1	0,-1	1,-1
	//16		1		1		16		-1,0			1,0
	//8		4	2		128	64	32		-1,1	0,1		1,1

	double e0, e1;
	int row, col, newRow, newCol;
	int numNeighbors = 0;
	int loopNum = 0;

	//maxDepth is maximum number of pixels in flow path length; 1000 used in Whitebox GAT (John Linsdsay)
	maxDepth = 10000;

	//somethingDone is boolean flag used to enter and exit do-while loop tracking flow path 
	bool somethingDone;

	//loop through all map rows 
	for (int row = 0; row < nRows; row++) {
		//loop through all map columns
		for (int col = 0; col < nCols; col++) {
			//e0 is DEM elevation_m_2D
			e0 = elevation_m_2D[row][col];
			//if e0 is real data
			if (e0 != NODATA_code) {
				//initialize numNeighbors counter to 0
				numNeighbors = 0;

				//loop through all 8 neighbors
				for (int k = 0; k < 8; k++)	{
					//newCol is neighbor col
					newCol = col + dx[k];
					//newRow is neighbor row
					newRow = row + dy[k];
					//checkfolder function confirms pixel is inside the map
					if (checkfolder(newRow, newCol)) {
						//e1 is elevation_m_2D of neighbor pixel
						e1 = elevation_m_2D[newRow][newCol];
						//if e1 neighbor is lower than e0 center pixel, and is real data
						if (e1 != NODATA_code && (e1 > e0)) {
							//increment numNeighbors counter
							numNeighbors++;
						}
					}
					else {
						//continue to next neighbor
						continue;
					}
				}
				//assign numNeighbor counter of downslope neighbors
				NeighborAnalysis_Map[row][col] = numNeighbors;
				InFlowDirection_Map[row][col] = numNeighbors;

			}
			else {
				NeighborAnalysis_Map[row][col] = NODATA_code;
				InFlowDirection_Map[row][col] = NODATA_code;
			}
		}
	}
	
	//Do loop while somethingDone = True to accumulate flow from upslope cells draining to pixel 
	do {
		//increment loopNum counter
		loopNum++;
		//reset somethingDone = False to exit loop if needed
		somethingDone = false;
		//loop through all map rows 
		for (row = 0; row < nRows; row++) {
			//loop through all map columns
			for (col = 0; col < nCols; col++) {
				//if getElevation is real data
				if (getElevation(row, col) != NODATA_code) {
					//if NeighborAnalysis_Map = 0, then pixel has not been checked
					if (NeighborAnalysis_Map[row][col] == 0) {   
						//initialize currentDepth to 0 to accumulate pixels
						currentDepth = 0;
						//call FlowAccumulation_MultipleFlowDirection function to find flow accumulation area
						FlowAccumulation_MultipleFlowDirection(row, col);
						//somethingDone set to True
						somethingDone = true;
					}
				}
				else {
					//continue to next pixel
					continue;
				}
			}
		}
	} while (somethingDone);
}

//FlowAccumulation_DInfinity function computes flow accumulation values for map with D-Infinity algorithm (Tarboton, 1997), taken from Whitebox GAT
//Reference: Tarboton, D. G. (1997). A New Method for the Determination of Flow Directions and Upslope Areas in Grid Digital Elevation Models. Water Resources Research, 33, 309-319. 
void TerrainProcessor::FlowAccumulation_DInfinity(int row, int col)
{
	double flowAccumVal = FlowAccumulation_Map[row][col];
	//flowDir_deg (deg) is FlowDirection_Map[row][col]; taken from 2D vector created by DInfinity
	double flowDir_deg = FlowDirection_Map[row][col];
	double tempAccum = 0;
	double proportion1 = 0;
	double proportion2 = 0;
	int a1 = 0;
	int b1 = 0;
	int a2 = 0;
	int b2 = 0;

	//currentDepth is incremented by 1 as flow path is evaluated
	currentDepth++;
	//if currentDepth exceeds maxDepth, then recursion exceeds limit and error likely
	if (currentDepth > maxDepth) {
		//break from routine
		return;
	}

	//NeighborAnalysis_Map set to -1 to indicate pixel flow path was evaluated; used to avoid revisiting same pixel
	NeighborAnalysis_Map[row][col] = -1; 

	//If flowDir_deg (deg) is greater than 0, search within the 360 deg array
	//Note: search for flow direction facets that receive flow, and compute proportion of flow received
	if (flowDir_deg >= 0) {

		if (flowDir_deg >= 0 && flowDir_deg < 45) {
			//proportion1 is fraction of flow from 45 degree facet
			proportion1 = (45 - flowDir_deg) / 45;
			//a1 is neighbor column
			a1 = col;
			//b1 is neighbor row
			b1 = row - 1;
			//proportion2 is fraction of flow from 45 degree facet
			proportion2 = flowDir_deg / 45;
			a2 = col + 1;
			b2 = row - 1;
		}
		else if (flowDir_deg >= 45 && flowDir_deg < 90) {
			proportion1 = (90 - flowDir_deg) / 45;
			a1 = col + 1;
			b1 = row - 1;
			proportion2 = (flowDir_deg - 45) / 45;
			a2 = col + 1;
			b2 = row;
		}
		else if (flowDir_deg >= 90 && flowDir_deg < 135) {
			proportion1 = (135 - flowDir_deg) / 45;
			a1 = col + 1;
			b1 = row;
			proportion2 = (flowDir_deg - 90) / 45;
			a2 = col + 1;
			b2 = row + 1;
		}
		else if (flowDir_deg >= 135 && flowDir_deg < 180) {
			proportion1 = (180 - flowDir_deg) / 45;
			a1 = col + 1;
			b1 = row + 1;
			proportion2 = (flowDir_deg - 135) / 45;
			a2 = col;
			b2 = row + 1;
		}
		else if (flowDir_deg >= 180 && flowDir_deg < 225) {
			proportion1 = (225 - flowDir_deg) / 45;
			a1 = col;
			b1 = row + 1;
			proportion2 = (flowDir_deg - 180) / 45;
			a2 = col - 1;
			b2 = row + 1;
		}
		else if (flowDir_deg >= 225 && flowDir_deg < 270) {
			proportion1 = (270 - flowDir_deg) / 45;
			a1 = col - 1;
			b1 = row + 1;
			proportion2 = (flowDir_deg - 225) / 45;
			a2 = col - 1;
			b2 = row;
		}
		else if (flowDir_deg >= 270 && flowDir_deg < 315) {
			proportion1 = (315 - flowDir_deg) / 45;
			a1 = col - 1;
			b1 = row;
			proportion2 = (flowDir_deg - 270) / 45;
			a2 = col - 1;
			b2 = row - 1;
		}
		else if (flowDir_deg >= 315 && flowDir_deg <= 360) {
			proportion1 = (360 - flowDir_deg) / 45;
			a1 = col - 1;
			b1 = row - 1;
			proportion2 = (flowDir_deg - 315) / 45;
			a2 = col;
			b2 = row - 1;
		}

		//checkfolder function confirms pixel is inside the map
		if (checkfolder(b1, a1)) {
			//if a proportion of flow is assigned to downslope cell along facet of flow path
			if (proportion1 > 0 && FlowAccumulation_Map[b1][a1] != NODATA_code) {
				tempAccum = FlowAccumulation_Map[b1][a1];

				//FlowAccumulation_Map map updated with inflow of flowAccum from neighboring cell
				FlowAccumulation_Map[b1][a1] = (tempAccum + (flowAccumVal * proportion1));

				//NeighborAnalysis_Map reduced by 1; used as indicator pixel flow path was searched
				NeighborAnalysis_Map[b1][a1]--;

				//recursive call to FlowAccumulation_DInfinity if NeighborAnalysis_Map is 0
				if (NeighborAnalysis_Map[b1][a1] == 0) {
					FlowAccumulation_DInfinity(b1, a1);
				}
			}
		}
		//checkfolder function confirms pixel is inside the map
		if (checkfolder(b2, a2)) {
			//if a proportion of flow is assigned to downslope cell along facet of flow path
			if (proportion2 > 0 && FlowAccumulation_Map[b2][a2] != NODATA_code) {
				tempAccum = FlowAccumulation_Map[b2][a2];

				//FlowAccumulation_Map map updated with inflow of flowAccum from neighboring cell
				FlowAccumulation_Map[b2][a2] = (tempAccum + (flowAccumVal * proportion2));

				//NeighborAnalysis_Map reduced by 1; used as indicator pixel flow path was searched
				NeighborAnalysis_Map[b2][a2]--;
				//NeighborAnalysis_Map[b2][a2] = (NeighborAnalysis_Map[b2][a2] - 1);

				//recursive call to FlowAccumulation_DInfinity if NeighborAnalysis_Map is 0
				if (NeighborAnalysis_Map[b2][a2] == 0) {
					FlowAccumulation_DInfinity(b2, a2);
				}
			}
		}
	}
	//currentDepth is reduced by 1
	currentDepth--;
}

//FlowAccumulation_MultipleFlowDirection function computes flow accumulation values with multiple flow direction algorithm (Quinn et al, 1991)
//Reference: Quinn, P., Beven, K., Chevallier, P., & Planchon, O. (1991). The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling using Digital Terrain Models. Hydrological Processes, 5, 59-79. 
//Reference: Quinn, P. F., Beven, K. J., & Lamb, R. (1995). The ln(a/Tan-Beta) Index - How To Calculate It and How To Use It Within The TOPMODEL Framework. Hydrological Processes, 9(2), 161-182. doi:10.1002/hyp.3360090204
void TerrainProcessor::FlowAccumulation_MultipleFlowDirection(int row, int col)
{

	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	//inflowDirVals inflow direction vector starts Northeast (NE=16), rotating clockwise to E=32, SE=64, S=128, SW=1, W=2, NW=4, N=8
	//Note: HydroPlus or Whitebox uses this inflowDirVals structure
	double inflowDirVals[8] = { 16, 32, 64, 128, 1, 2, 4, 8 };

	//Note: HydroPlus or WhiteBox flowDir is 1 for NE corner, 2 for E, 4 for SE, 8 for S, 16 for SW, 32 for W, 64 for NW, 128 for N
	//Whitebox FDR		Whitebox -FDR	Pairs of x,y 
	//64	128	1		4	8	16		-1,-1	0,-1	1,-1
	//32		2		2		32		-1,0			1,0
	//16	8	4		1	128	64		-1,1	0,1		1,1

	//Warning: ESRI ArcGIS method not used in HydroPlus, remove or leave commented out
	//Note: ESRI ArcGIS flow direction vectors arranged clockwise: x=1y=0 in E (right middle), x=1,y=-1 in NE (top right)
	//int dx[8] = { 1, 1, 0, -1, -1, -1, 0, 1 };
	//int dy[8] = { 0, 1, 1, 1, 0, -1, -1, -1 };

	//ESRI ArcGIS Pairs of x, y
	//0, -1		0, -1	1, 0
	//-1, -1			1, 1
	//-1, 0		-1, 1	0, 1

	//inflowDirVals inflow direction vector starts Northeast (NE=32), rotating clockwise to E=64, SE=128, S=1, SW=2, W=4, NW=8, N=16
	//Note: ESRI ArcGIS uses this inflowDirVals structure
	//double inflowDirVals[8] = {32, 64, 128, 1, 2, 4, 8, 16 };

	//Note: ArcGIS ESRI flowDir is 128 for NE corner, 1 for E, 2 for SE, 4 for S, 8 for SW, 16 for W, 32 for NW, 64 for N
	//ESRI FDR			ESRI -FDR		Pairs of x,y 
	//32	64	128		2	4	8		-1,-1	0,-1	1,-1
	//16		1		1		16		-1,0			1,0
	//8		4	2		128	64	32		-1,1	0,1		1,1

	double flowAccumVal = FlowAccumulation_Map[row][col];
	double flowDir = FlowDirection_Map[row][col];
	double tempAccum = 0;
	double temporganizerVal = 0;
	double proportion = 0;
	double totalRelief = 0;
	double power = 1.0;
	int newRow, newCol = 0;
	double e0 = 0;
	double e1 = 0;

	//currentDepth is incremented by 1 as flow path is evaluated
	currentDepth++;
	//if currentDepth exceeds maxDepth, then recursion exceeds limit and DEM may have error 
	if (currentDepth > maxDepth) {
		return;
	}

	//NeighborAnalysis_Map set to -1 to indicate pixel flow path was evaluated; used to avoid revisiting same pixel
	NeighborAnalysis_Map[row][col] = -1;

	e0 = elevation_m_2D[row][col];

	//search for flow direction facets that receive flow, and compute proportion of flow received
	//loop through all 8 neighbors
	for (int k = 0; k < 8; k++)	{
		//newRow is neighbor row
		newRow = row + dy[k];
		//newCol is neighbor column
		newCol = col + dx[k];

		//checkfolder function confirms pixel is inside the map
		if (checkfolder(newRow, newCol)) {
			//e1 is neighbor elevation_m_2D
			e1 = elevation_m_2D[newRow][newCol];
			//if e1 neighbor is lower than e0 center pixel, and is real data
			if ((e0 > e1) && (e1 != NODATA_code)) {
				//accumulate totalRelief
				totalRelief += pow((e0 - e1), power);
			}
		}
		else {
			//continue to next neighbor
			continue;
		}
	}
	//loop through all 8 neighbors
	for (int k = 0; k < 8; k++) {
		//newRow is neighbor row
		newRow = row + dy[k];
		//newCol is neighbor column
		newCol = col + dx[k];
		double new_FA;

		//checkfolder function confirms pixel is inside the map
		if (checkfolder(newRow, newCol)) {

			//e1 is neighbor elevation_m_2D
			e1 = elevation_m_2D[newRow][newCol];
			//if e1 neighbor is lower than e0 center pixel, and is real data
			if ((e0 > e1) && (e1 != NODATA_code)) {
				//Check algorithm with Quinn et al. (1991, 1995)
				//Quinn, P., Beven, K., Chevallier, P., & Planchon, O. (1991). The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling using Digital Terrain Models. Hydrological Processes, 5, 59-79. 
				//Quinn, P. F., Beven, K. J., & Lamb, R. (1995). The ln(a/Tan-Beta) Index - How To Calculate It and How To Use It Within The TOPMODEL Framework. Hydrological Processes, 9(2), 161-182. doi:10.1002/hyp.3360090204
				proportion = pow((e0 - e1), power) / totalRelief;
				tempAccum = FlowAccumulation_Map[newRow][newCol];

				//new_FA is fraction of cell flowing into neighbor as determined by product of flowAccumVal and proportion
				new_FA = (tempAccum + (flowAccumVal * proportion));

				//FlowAccumulation_Map is increased by proportion of accumulated flow 
				FlowAccumulation_Map[newRow][newCol] = new_FA;

				temporganizerVal = NeighborAnalysis_Map[newRow][newCol];
				//NeighborAnalysis_Map reduced by 1; used as indicator pixel flow path was searched
				NeighborAnalysis_Map[newRow][newCol]--;

				//recursive call to FlowAccumulation_MultipleFlowDirection if NeighborAnalysis_Map is 0
				if (NeighborAnalysis_Map[newRow][newCol] == 0)	{
					FlowAccumulation_MultipleFlowDirection(newRow, newCol);
				}
			}
		}
		else {
			//continue to next neighbor
			continue;
		}
	}
	//currentDepth is reduced by 1
	currentDepth--;
}

//TopographicIndex_Map_Compute function for exponential and power decay, but TopographicIndex_Log_to_Power() function modifies power decay values
//Note: Computes topographic index for each map pixel, using RHS terms computed in TopographicIndex_Map_Compute function
void TerrainProcessor::TopographicIndex_Map_Compute()
{
	//loop through all map rows 
	for (int row = 0; row < nRows; row++) {
		//loop through all map columns
		for (int col = 0; col < nCols; col++) {

			//If FlowAccumulation_Map is not NODATA_code and FlowAccumulation_Map is not 0 then enter functions to divide by FlowAccumulation_Map
			if (FlowAccumulation_Map[row][col] != NODATA_code && FlowAccumulation_Map[row][col] != 0) {
				//If stat_PC_frac < 0 then set to 0.01 to prevent mathematical error of log(0)
				if (stat_PC_frac < 0) {
					stat_PC_frac = 0.01;
				}
				//If Flag_TI_PerviousOnly=1, then use only pervious area in FlowAccumulation_Map when computing topographic index in statistical mode
				//Note: HydroPlusConfig.xml Flag_TI_PerviousOnly=1 excludes impervious area in statistical or spatial mode using TopUrban theory of Valeo and Moin (2000)
				if (Inputs::SimulationNumericalParams["Flag_TI_PerviousOnly"] == 1) {
					//TopographicIndex_Map is natural log of quotient of contributing area and tangent of slope; C++ log is base e
					TopographicIndex_Map[row][col] = log((FlowAccumulation_Map[row][col] * stat_PC_frac) / (tan(Slope_Ground_Horn_rad_Map[row][col])));
				}
				//Else if Flag_TI_Pervious_Accumulated not equal to one then accumulate impervious and pervious area
				else {
					//TopographicIndex_Map is natural log of quotient of contributing area and tangent of slope; C++ log is base e
					TopographicIndex_Map[row][col] = log((FlowAccumulation_Map[row][col]) / (tan(Slope_Ground_Horn_rad_Map[row][col])));
				}
			}
			//Else If FlowAccumulation_Map equals 0 then TopographicIndex_Map set to 0
			else if(FlowAccumulation_Map[row][col] == 0)  {
				TopographicIndex_Map[row][col] = 0;
			}
			//Else If FlowAccumulation_Map is NODATA_code then TopographicIndex_Map set to NODATA_code
			else {
				TopographicIndex_Map[row][col] = NODATA_code;
			}
		}
	}
}

//Topographic index algorithm taken from TOPMODEL code developed by Keith Beven and used by USGS D.M. Wolock https://pubs.er.usgs.gov/publication/wri934124
//Reference: Wolock, D. M. (1993). Simulating the variable-source-area concept of watershed hydrology with TOPMODEL, USGS WRI Report 93-4124).
//Note: The algorithm requires the topographic index (TI) bins (ia) are sorted from highest (wettest) to lowest TI values
//Note: The algorithm requires when ia = 0, the TI bin TI_Value[0] should be the upper limit and its TI_Area_frac[0] = 0
//Note: The algorithm presumes for each ia, the area of catchment between discrete TI values TI_Value[ia] and TI_Value[ia+1] equals TI_Area_frac[ia]
//Note: The algorithm presumes for each ia, the area of catchment bounding TI_Value[ia] equals (TI_Area_frac[ia] + TI_Area_frac[ia+1])/2 

//TopographicIndex_Map_to_Histogram function will generate histogram of topographic index
void TerrainProcessor::TopographicIndex_Map_to_Histogram()
{
	//initialize TI_Value as bin number
	int TI_Value_int = 0;
	//initialize TI_Value_max as topographic wetness TI_Value maximum, where TWI is another name for TI
	double TI_Value_max = 0.0;
	//initialize TI_Value_min as topographic wetness TI_Value minimum, where TWI is another name for TI
	double TI_Value_min = 10000;
	double TI_Values_In_Map = 0.0;

	cout << "Generating the Topographic Index histogram with topographic bins = " << NAC << endl;

	//TI_Values_In_HistogramBin initialized, NAC is number of TI bins
	TI_Values_In_HistogramBin.assign(NAC, 0.0);
	//TI_Area_frac initialized, NAC is number of TI bins
	TI_Area_frac.assign(NAC, 0.0);
	//TI_Value initialized, NAC is number of TI bins
	TI_Value.assign(NAC, 0.0);

	//TI_validCount initialized ---- Find min / max and count valid pixels ----
	int TI_validCount = 0;
	//For row to nRows, size of TopographicIndex_Map array
	for (int row = 0; row < nRows; row++) {
		//For col to nCols, size of TopographicIndex_Map array
		for (int col = 0; col < nCols; col++) {
			//If TopographicIndex_Map[row][col] != NODATA_code then enter
			if (TopographicIndex_Map[row][col] != NODATA_code) {
				//TI_Value_row_col = TopographicIndex_Map[row][col], assign map value to double
				double TI_Value_row_col = TopographicIndex_Map[row][col];
				//TI_validCount advanced 1
				TI_validCount++;
				//Find Max value: If (TI_Value_row_col > TI_Value_max) then TI_Value_max = TI_Value_row_col
				if (TI_Value_row_col > TI_Value_max) {
					TI_Value_max = TI_Value_row_col;
				}
				//Find Min value: If (TI_Value_row_col < TI_Value_min) then TI_Value_min = TI_Value_row_col
				if (TI_Value_row_col < TI_Value_min) {
					TI_Value_min = TI_Value_row_col;
				}
			}
		}
	}
	//Error handling if TI_validCount == 0
	if (TI_validCount == 0) {
		cerr << "ERROR: No valid TI values found in TopographicIndex_Map.\n";
		abort();
	}

	//Special Case
	//If fabs(TI_Value_max - TI_Value_min) < 1e-9 then uniform TI map, all TI values the same
	if (fabs(TI_Value_max - TI_Value_min) < 1e-9) {
		//Alert user, in case this was not intended
		cout << "Note: TI map has uniform value " << TI_Value_max
			<< ". Using equal-area TI histogram.\n";

		//TI_Value for bin 0 is reserved for TI_Value_max with zero area
		TI_Value[0] = TI_Value_max;
		TI_Area_frac[0] = 0.0;

		//Distribute area equally among remaining bins
		if (NAC > 1) {
			double TI_Area_i_frac = 1.0 / static_cast<double>(NAC - 1);
			for (int k = 1; k < NAC; ++k) {
				//All TI_Value are equal; TI_Value[k] = TI_Value_max
				TI_Value[k] = TI_Value_max;  
				//All TI_Area_frac are equal; TI_Area_frac[k] = TI_Area_i_frac
				TI_Area_frac[k] = TI_Area_i_frac;
			}
		}

		//Average index is just that uniform value; TI_Value_Average = TI_Value_max
		TI_Value_Average = TI_Value_max;
		return;    // skip the usual incremental-bin logic
	}

	//Normal Case; different values across TI map and histogram
	//TI_HistogramBin_Increment is incremental distance between TWI max and min, across the number of bins
	//TI_HistogramBin_Increment will create breakpoints in TI values within histogram
	int TI_HistogramBins_MinusOne = NAC - 1;
	double TI_HistogramBin_Increment = (TI_Value_max - TI_Value_min) / (TI_HistogramBins_MinusOne);

	//create TI histogram
	//loop through rows of topographic index map, along row
	for (int row = 0; row < nRows; row++) {
		//loop through cols of topographic index map, along col
		for (int col = 0; col < nCols; col++) {
			//if real topographic index map values exist, and not NODATA_code
			if (TopographicIndex_Map[row][col] != NODATA_code) {
				//TI_Value is integer of scaled TI values used to put TI values into the NAC topographic bins, ranging from 1 to NAC
				//Note: Ex w/ NAC=30 for Exponential Decay: TWIOrganizer[row][col] = 15, TI_Value_max = 30, TI_Value_min = 3, TI_HistogramBin_Increment = 0.931, then TI_Value = int(13.89) = 14
				//Note: Ex w/ NAC=30 for Power Decay: TWIOrganizer[row][col] = 1500, TI_Value_max = 3000, TI_Value_min = 3, TI_HistogramBin_Increment = 103.345, then TI_Value = int(15.48) = 15
				//Note: Any TopographicIndex_Map value equal to TI_Value_min becomes an TI_Value value of 1
				//Note: Any TopographicIndex_Map value equal to TI_Value_max becomes an TI_Value value of TI_Value_max
				//Note: Puts entire map of TI values into discrete bins
				TI_Value_int = static_cast<int>(((TopographicIndex_Map[row][col] - TI_Value_min) / TI_HistogramBin_Increment) + 1);

				//TI_Value should never be larger than NAC-1 (equal to NAC) which is reserved for TI_Value_max
				if (TI_Value_int > NAC - 1) {
					TI_Value_int = NAC - 1;
				}

				//TI_Values_In_HistogramBin[TI_Value_int] counts number of TopographicIndex_Map values going into each bin
				TI_Values_In_HistogramBin[TI_Value_int] += 1;
				//TI_Values_In_Map counts all TopographicIndex_Map values, used to scale TI_Values_In_HistogramBin[TI_Value] into fractions of TI_Values_In_Map
				TI_Values_In_Map += 1;
			}
		}
	}

	//TI_Value[0] takes the maximum TWI value, which can be used later in TI algorithms 
	TI_Value[0] = TI_Value_max;

	//loop through all TI bins in histogram
	for (int TI_HistogramBin_ID = 1; TI_HistogramBin_ID < NAC; TI_HistogramBin_ID++) {
		int TI_HistogramBins_MinusOne = NAC - TI_HistogramBin_ID;
		//TI_Area_frac[TI_HistogramBin_ID] is area (fraction) of catchment for each bin, TI_Values_In_HistogramBin[] is counter of occurrence in map of TI values in each bin
		//Note: N.B. TI_Area_frac[0] = 0, where 0 is bin with TI_Value_max value
		TI_Area_frac[TI_HistogramBin_ID] = TI_Values_In_HistogramBin[TI_HistogramBins_MinusOne] / TI_Values_In_Map;
		//TI_Value[TI_HistogramBin_ID] is TI value for each bin, constructed with TI_HistogramBin_Increment based on incremental distance between TI_Value_max and TI_Value_min, 
		TI_Value[TI_HistogramBin_ID] = TI_Value_max - (TI_HistogramBin_ID)*TI_HistogramBin_Increment;
	}
}

//TopographicIndex_Log_to_Power function converts topographic index (TI) values in map and histogram from exponential decay ln formulation to power decay formulation
//Note: TI transformed by taking exponential of ln(a/tanB), raised to power (1/n), units change from ln([L]) to [L]^(1/n)
void TerrainProcessor::TopographicIndex_Log_to_Power()
{
	double nN = input->InputXml["Parameter_n_KsatPowerDecay"];
	//loop through all rows of map
	for (int row = 0; row < nRows; row++) {
		//loop through all cols of map
		for (int col = 0; col < nCols; col++) {
			//for pixels with NODATA_code, keep NODATA_code
			if (TopographicIndex_Map[row][col] == NODATA_code) {
				TopographicIndex_Map[row][col] = NODATA_code;
			}
			//for pixels with data, TopographicIndex_Log_to_Power 
			else {
				//TopographicIndex_Map is TI map value transformed value by taking exponential of ln(a/tanB), raised to power (1/n)
				TopographicIndex_Map[row][col] = pow(exp(TopographicIndex_Map[row][col]), 1 / nN);
			}
			//remove zero values
			if (TopographicIndex_Map[row][col] < 0) {
				TopographicIndex_Map[row][col] = 0.01;
			}

		}
	}

	//loop through all TI bins in histogram
	for (int TI_HistogramBin_ID = 0; TI_HistogramBin_ID < NAC; TI_HistogramBin_ID++) {
		//TI_Value is TI histogram value transformed value by taking exponential of ln(a/tanB), raised to power (1/n)
		TI_Value[TI_HistogramBin_ID] = pow(exp(TI_Value[TI_HistogramBin_ID]), 1 / nN);
	}
}

//TopographicCharacteristics_Vector_to_Map function converts DEM byproducts (TI, slope, aspect, FD, FA) from vectors to maps
//Note: Maps called z_TIorganizer.asc, z_SlopeGround_rad.asc, z_AspectGround_N_0_rad.asc, z_FDorganizer.asc, z_FlowAccum.asc
void TerrainProcessor::TopographicCharacteristics_Vector_to_Map()
{
	cout << "Writing or confirming files with processed topographic data." << endl;

	// Outputting the z_TIorganizer.asc
	ofstream outfile_TopographicIndex(Directory_Input_CLArg + "z_TIorganizer.asc");
	if (!outfile_TopographicIndex.good()) {
		cout << "Warning: Could not create " << Directory_Input_CLArg << "z_TIorganizer.asc" << endl;
	}

	outfile_TopographicIndex << "ncols" << setw(20) << Inputs::nCols << endl;
	outfile_TopographicIndex << "nrows" << setw(20) << Inputs::nRows << endl;
	outfile_TopographicIndex << "xllcorner" << setw(16) << Inputs::xllcorner_m << endl;
	outfile_TopographicIndex << "yllcorner" << setw(16) << Inputs::yllcorner_m << endl;
	outfile_TopographicIndex << "cellsize" << setw(17) << Inputs::Length_Pixel_Side_m << endl;
	outfile_TopographicIndex << "NODATA_value" << setw(13) << Inputs::NODATA_code << endl;

	//TI_sum initialized as zero
	float TI_sum = 0;
	//rem initialized as modulus of nCols divided by 8
	int rem = nCols % 8;
	//for row = 0 to row < nRows, incrementing row
	for (int row = 0; row < nRows; row++) {
		//for col = 0 to col < nCols, incrementing col
		for (int col = 0; col < nCols; col++) {
			//If elevation_m_2D is not NODATA_code
			if (elevation_m_2D[row][col] != NODATA_code) {
				//outfile_TopographicIndex has written TopographicIndex_Map value
				outfile_TopographicIndex << " " << setw(10) << TopographicIndex_Map[row][col];
				//TI_sum is incremented by value of TopographicIndex_Map[row][col]
				TI_sum = TI_sum + TopographicIndex_Map[row][col];
			}
			//Else outfile_TopographicIndex has written NODATA_code
			else {
				outfile_TopographicIndex << " " << setw(10) << NODATA_code;
			}
		}
		outfile_TopographicIndex << endl;
	}
	// End of outputting the z_TIorganizer.asc

	//If TI_sum <= 0 then no TI values
	if (TI_sum <= 0) {
		cout << "Warning: All Topographic Index Values are zero in the z_TIorganizer.asc file." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: When the z_TIorganizer.asc has all zero values, there is no subsurface hydrological funcationality." << endl;
		cout << "Correction: Check your DEM.asc file or modify HydroPlusConfig.xml parameters such as Flag_TI_PerviousOnly = 0." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	ofstream outfile_SlopeGround_rad(Directory_Input_CLArg + "z_SlopeGround_rad.asc");
	if (!outfile_SlopeGround_rad.good()) {
		cout << "Warning: Could not create " << Directory_Input_CLArg << "z_SlopeGround_rad.asc" << endl;
	}
	outfile_SlopeGround_rad << left << setw(14) << "ncols" << setw(15) << Inputs::nCols << endl;
	outfile_SlopeGround_rad << left << setw(14) << "nrows" << setw(15) << Inputs::nRows << endl;
	outfile_SlopeGround_rad << left << setw(14) << "xllcorner" << setw(15) << setprecision(15) << Inputs::xllcorner_m << endl;
	outfile_SlopeGround_rad << left << setw(14) << "yllcorner" << setw(15) << setprecision(15) << Inputs::yllcorner_m << endl;
	outfile_SlopeGround_rad << left << setw(14) << "cellsize" << setw(15) << setprecision(15) << Inputs::Length_Pixel_Side_m << endl;
	outfile_SlopeGround_rad << left << setw(14) << "NODATA_value" << setw(15) << Inputs::NODATA_code << endl;

	for (int row = 0; row < nRows; row++) {
		for (int col = 0; col < nCols; col++) {
			if (elevation_m_2D[row][col] != NODATA_code) {
				//Slope_Ground_Horn_rad_Map takes Slope_Ground_Horn_rad_Map value
				outfile_SlopeGround_rad << Slope_Ground_Horn_rad_Map[row][col] << " ";
			}
			else {
				outfile_SlopeGround_rad << NODATA_code << " ";
			}
		}
		outfile_SlopeGround_rad << endl;
	}
	outfile_SlopeGround_rad.close();
	// End of outputting the slope organizer

	// Outputting the aspect organizer 
	ofstream outfile_AspectGround_N_0_rad(Directory_Input_CLArg + "z_AspectGround_N_0_rad.asc");
	if (!outfile_AspectGround_N_0_rad.good())	{
		cout << "Warning: Could not create " << Directory_Input_CLArg << "z_AspectGround_N_0_rad.asc" << endl;
	}
	outfile_AspectGround_N_0_rad << left << setw(14) << "ncols" << setw(15) << Inputs::nCols << endl;
	outfile_AspectGround_N_0_rad << left << setw(14) << "nrows" << setw(15) << Inputs::nRows << endl;
	outfile_AspectGround_N_0_rad << left << setw(14) << "xllcorner" << setw(15) << setprecision(15) << Inputs::xllcorner_m << endl;
	outfile_AspectGround_N_0_rad << left << setw(14) << "yllcorner" << setw(15) << setprecision(15) << Inputs::yllcorner_m << endl;
	outfile_AspectGround_N_0_rad << left << setw(14) << "cellsize" << setw(15) << setprecision(15) << Inputs::Length_Pixel_Side_m << endl;
	outfile_AspectGround_N_0_rad << left << setw(14) << "NODATA_value" << setw(15) << Inputs::NODATA_code << endl;

	for (int row = 0; row < nRows; row++) {
		for (int col = 0; col < nCols; col++) {
			if (elevation_m_2D[row][col] != NODATA_code) {
				outfile_AspectGround_N_0_rad << Aspect_Ground_N_0_rad_Map[row][col] << " ";
			}
			else {
				outfile_AspectGround_N_0_rad << NODATA_code << " ";
			}
		}
		outfile_AspectGround_N_0_rad << endl;
	}
	outfile_AspectGround_N_0_rad.close();
	// End of outputting the aspect organizer

	// outputting the flowdirection organizer	
	ofstream outfile_FlowDirection(Directory_Input_CLArg + "z_FDorganizer.asc");
	if (!outfile_FlowDirection.good()) {
		cout << "Warning: Could not create " << Directory_Input_CLArg << "z_FDorganizer.asc" << endl;
	}
	outfile_FlowDirection << left << setw(14) << "ncols" << setw(15) << Inputs::nCols << endl;
	outfile_FlowDirection << left << setw(14) << "nrows" << setw(15) << Inputs::nRows << endl;
	outfile_FlowDirection << left << setw(14) << "xllcorner" << setw(15) << setprecision(15) << Inputs::xllcorner_m << endl;
	outfile_FlowDirection << left << setw(14) << "yllcorner" << setw(15) << setprecision(15) << Inputs::yllcorner_m << endl;
	outfile_FlowDirection << left << setw(14) << "cellsize" << setw(15) << setprecision(15) << Inputs::Length_Pixel_Side_m << endl;
	outfile_FlowDirection << left << setw(14) << "NODATA_value" << setw(15) << Inputs::NODATA_code << endl;

	for (int row = 0; row < nRows; row++) {
		for (int col = 0; col < nCols; col++) {
			if (elevation_m_2D[row][col] != NODATA_code) {
				outfile_FlowDirection << FlowDirection_Map[row][col] << " ";
			}
			else {
				outfile_FlowDirection << NODATA_code << " ";
			}
		}
		outfile_FlowDirection << endl;
	}
	outfile_FlowDirection.close();
	// End of outputting the z_FDorganizer.asc

	// Outputting the flow accumulation organizer	
	ofstream outfile_FlowAccumulation(Directory_Input_CLArg + "z_FlowAccum.asc");
	if (!outfile_FlowAccumulation.good()) {
		cout << "Warning: Could not create " << Directory_Input_CLArg << "z_FlowAccum.asc" << endl;
	}
	outfile_FlowAccumulation << left << setw(14) << "ncols" << setw(15) << Inputs::nCols << endl;
	outfile_FlowAccumulation << left << setw(14) << "nrows" << setw(15) << Inputs::nRows << endl;
	outfile_FlowAccumulation << left << setw(14) << "xllcorner" << setw(15) << setprecision(15) << Inputs::xllcorner_m << endl;
	outfile_FlowAccumulation << left << setw(14) << "yllcorner" << setw(15) << setprecision(15) << Inputs::yllcorner_m << endl;
	outfile_FlowAccumulation << left << setw(14) << "cellsize" << setw(15) << setprecision(15) << Inputs::Length_Pixel_Side_m << endl;
	outfile_FlowAccumulation << left << setw(14) << "NODATA_value" << setw(15) << Inputs::NODATA_code << endl;

	for (int row = 0; row < nRows; row++) {
		for (int col = 0; col < nCols; col++) {
			if (elevation_m_2D[row][col] != NODATA_code) {
				outfile_FlowAccumulation << FlowAccumulation_Map[row][col] << " ";
			}
			else {
				outfile_FlowAccumulation << NODATA_code << " ";
			}
		}
		outfile_FlowAccumulation << endl;
	}
	outfile_FlowAccumulation.close();
	// End of outputting the flow accumulation organizer
}

//To find the maximum value in a 2D vector while excluding specific values like -9999 (nodata), -nan, inf, and -inf
double TerrainProcessor::max_2d_element(vector<vector<double>> vec2D) {
	double max = DBL_MIN;

	for (const auto& row : vec2D) {
		for (const double& value : row) {
			if (value != NODATA_code && !isnan(value) && isfinite(value)) {
				if (value > max) {
					max = value;
				}
			}
		}
	}
	/*for (auto it = vec2D.begin(); it != vec2D.end(); it++) {
		double max_temp = *max_element((*it).begin(), (*it).end());
		if (max_temp != NODATA_code && !isnan(max_temp) && isfinite(max_temp)) {
			if (max_temp > max) {
				max = max_temp;
			}	
		}
	}*/
	return max;
}
//To find the minimum value in a 2D vector while excluding specific values like -9999 (nodata), -nan, inf, and -inf
double TerrainProcessor::min_2d_element(vector<vector<double>> vec2D) {
	double min = DBL_MAX;

	for (const auto& row : vec2D) {
		for (const double& value : row) {
			if (value != NODATA_code && !isnan(value) && isfinite(value)) {
				if (value < min) {
					min = value;
				}
			}
		}
	}
	/*for (auto it = vec2D.begin(); it != vec2D.end(); it++) {
		double min_temp = *min_element((*it).begin(), (*it).end());
		if (min_temp != NODATA_code && !isnan(min_temp) && isfinite(min_temp)) {
			if (min_temp < min) {
				min = min_temp;
			}
		}
	}*/
	return min;
}

template <class Q>
void clearQueue(Q& q) {
	q = Q();
}

void TerrainProcessor::Find_pit(vector<vector<double>> elevation_map) {

	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	// Find pit cells. This step is parallelized.
	undefined_flow_cells.clear();
	for (int row = 1; row < nRows - 1; row++) {
		for (int col = 1; col < nCols - 1; col++) {
			double z = elevation_map[row][col];
			if (z != NODATA_code) {
				bool flag_findpit = true;

				for (int n = 0; n < 8; n++) {
					double zn = elevation_map[row + dy[n]][col + dx[n]];
					if (zn < z || zn == NODATA_code) {
						// It has a lower neighbour or is an edge cell.
						flag_findpit = false;
						break;
					}
				}
				if (flag_findpit) {
					// it's a cell with undefined flow
					undefined_flow_cells.push_back({ row, col, z });
				}
			}
		}
	}
}


void TerrainProcessor::FillDepressions(vector<vector<double>> elevation_working_m_2D) {

	//Note: HydroPlus or Whitebox flow direction vectors arranged clockwise: x=1,y=-1 in NE (top right), x=0,y=-1 in N (top middle)
	int dx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 };
	int dy[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//HydroPlus or Whitebox Pairs of x, y
	//- 1, -1	0, -1	1, -1
	//- 1, 0			1, 0
	//- 1, 1	0, 1	1, 1

	//Each pitncell is visited from highest to lowest
	sort(undefined_flow_cells.begin(), undefined_flow_cells.end(),
		[](const auto& a, const auto& b) {
			return a.z > b.z;
		});
	
	//Optional maximum depression depth to fill.
	double max_depth = max_2d_element(elevation_working_m_2D) - min_2d_element(elevation_working_m_2D);
	//cout << " max_depth is " << max_depth << " = " << max_2d_element(elevation_working_m_2D) << " - " << min_2d_element(elevation_working_m_2D) << endl;
	// Now we need to perform an in-place depression filling
	vector< vector <int>> visited(nRows, vector<int>(nCols, 0));
	deque<pair<int, int>> queue;
	flats.clear();
	flats.resize(nRows, vector<int>(nCols, 0));

	//cout << "Filling depressions..." << endl;
	possible_outlets.clear();
	//go through the undefined_flow_cells list
	while (!undefined_flow_cells.empty()) {
		//take the value of this cell
		auto cell = undefined_flow_cells.back();
		//pop out this cell from the vector
		undefined_flow_cells.pop_back();

		int row = cell.row;
		int col = cell.col;
		//cout << "Filling cell row: " << row  + 1 << ", col: " << col + 1<< ", flat[][] is " << flats[row][col] << endl;
		//if it's not a flat region
		if (flats[row][col] != 1) {
			//cout << "row: " << row + 1 << ", col : " << col + 1 << " is not a flat region" << endl;
			// First there is a priority region-growing operation to find the outlets.
			double z_pit = elevation_depression_filled_m_2D[row][col];
			//minheap is a structure to contain a group of numbers, row, column, z
			priority_queue<GridCell> minheap;
			//push the row, col and elevation to minheap
			minheap.push({ row, col, z_pit });
			//cout << "minheap pushing back " << row << ", " << col << endl;
			//mark this cell as visited
			visited[row][col] = 1;
			//initilize the outlet_found as false and outlet_z as a infinity number
			bool outlet_found = false;
			double outlet_z = numeric_limits<double>::infinity();

			if (!queue.empty()) {
				queue.clear();
			}
			// if this minheap is not empty 
			while (!minheap.empty()) {
				GridCell cell2 = minheap.top();
				minheap.pop();
				double z = cell2.priority;

				if (outlet_found && z > outlet_z) {
					break;
				}

				if (z - z_pit > max_depth) {
					// No outlet could be found that was low enough.
					//cout << "z - z_pit is "<< z - z_pit <<  ",> max_depth. No outlet could be found that was low enough." << endl;
					break;
				}
				//if haven't found the outlet
				if (!outlet_found) {
					//cout << "Outlet haven't found, loop through the neighbour cells of " << cell2.row << ", " << cell2.column << endl;
					//loop through its neighbour cells
					for (int n = 0; n < 8; n++) {
						int cn = cell2.column + dx[n];
						int rn = cell2.row + dy[n];
						//cout << "Verifying cn: " << cn << "; rn: " << rn << endl;
						
						//if this neighbour cell haven't been visited
						if (checkfolder(rn, cn)) {
							if (visited[rn][cn] == 0) {
								double zn = elevation_depression_filled_m_2D[rn][cn];
							
								if (!outlet_found) {
									//we haven't found outlet
									//cout << "Outlet not found " << endl;
									if (zn >= z && zn != NODATA_code) {
										minheap.push({ rn, cn, zn });
										//cout << "minheap pushing back " << rn << ", " << cn << endl;
										visited[rn][cn] = 1;
									}
									else if (zn != NODATA_code) {
										// zn < z
										// 'cell' has a lower neighbour that hasn't already passed through minheap.
										// Therefore, 'cell' is a pour point cell.
										outlet_found = true;
										outlet_z = z;
										queue.push_back({ cell2.row, cell2.column });
										possible_outlets.push_back({ cell2.row, cell2.column });
									}
								}
								else if (zn == outlet_z) {
									// We've found the outlet but are still looking for additional depression cells.
									//cout << "2367 Outlet found but are still looking for additional depression cells" << endl;
									minheap.push({ rn, cn, zn });
									// cout << "minheap pushing back " << rn << ", " << cn << endl;
									visited[rn][cn] = 1;
								}
							}
						}	
					}
				}
				else {
					// We've found the outlet but are still looking for additional depression cells and potential outlets.
					//cout << "2378 Outlet found but are still looking for additional depression cells" << endl;
					if (z == outlet_z) {
						bool flag = false;

						for (int n = 0; n < 8; n++) {
							int cn = cell2.column + dx[n];
							int rn = cell2.row + dy[n];

							if (checkfolder(rn, cn)) {
								if (visited[rn][cn] == 0) {
									double zn = elevation_depression_filled_m_2D[rn][cn];

									if (zn < z) {
										flag = true;
									}
									else if (zn == outlet_z) {
										minheap.push({ rn, cn, zn });
										visited[rn][cn] = 1;
									}
								}
							}	
						}

						if (flag) {
							// it's an outlet
							queue.push_back({ cell2.row, cell2.column });
							possible_outlets.push_back({ cell2.row, cell2.column });
						}
						else {
							visited[cell2.row][cell2.column] = 1;
						}
					}
				}
			}

			if (outlet_found) {
				//cout << "Outlet found" << endl;
				// Now that we have the outlets, raise the interior of the depression.
				// Start from the outlets.
				priority_queue<GridCell> minheap;
				while (!queue.empty()) {
					auto cell2 = queue.front();
					queue.pop_front();
					for (int n = 0; n < 8; n++) {
						int rn = cell2.first + dy[n];
						int cn = cell2.second + dx[n];

						if (checkfolder(rn, cn)) {
							if (visited[rn][cn] == 1) {
								visited[rn][cn] = 0;
								queue.push_back({ rn, cn });
								double z = elevation_depression_filled_m_2D[rn][cn];
								//cout << "z is " << z << ", outlet_z is " << outlet_z << endl;
								if (z < outlet_z) {
									elevation_depression_filled_m_2D[rn][cn] = outlet_z;
									//initialize elevation_flats_filled_m_2D as elevation_depression_filled_m_2D
									elevation_flats_filled_m_2D[rn][cn] = outlet_z;
									flats[rn][cn] = 1;
									//cout << "Pit cell row: " << rn + 1 << ", col: " << cn + 1 << ", increase elevation from " << z << " to: " << outlet_z << endl;
									//cout << "[4,8] is: " << elevation_flats_filled_m_2D[3][7] << endl;
									//cout << "Flat change to 1, row: " << rn + 1 << ", col: " << cn + 1 << endl;
								}
								else if (z == outlet_z) {
									flats[rn][cn] = 1;
									//cout << "Flat found, row: " << rn + 1 << ", col: " << cn + 1 << endl;
									//cout << "Flat change to 1, row: " << rn + 1 << ", col: " << cn + 1 << endl;
								}
							}
						}	
					}
				}
			}
			else {
				queue.push_back({ row, col }); // start at the pit cell and clean up visited

				while (!queue.empty()) {
					auto cell2 = queue.front();
					queue.pop_front();

					for (int n = 0; n < 8; n++) {
						int rn = cell2.first + dy[n];
						int cn = cell2.second + dx[n];

						if (checkfolder(rn, cn)) {
							if (visited[rn][cn] == 1) {
								visited[rn][cn] = 0;
								queue.push_back({ rn, cn });
							}
						}			
					}
				}
			}
		}
	}
}
	
void TerrainProcessor::FixFlats(vector<vector<double>> elevation_working_m_2D, double small_num) {

	//cout << "Fixing flow on flats..." << endl;

	vector<int> dx = { 1, 1, 1, 0, -1, -1, -1, 0 };
	vector<int> dy = { -1, 0, 1, 1, 1, 0, -1, -1 };

	//cout << "[4,8] is: " << elevation_flats_filled_m_2D[3][7] << endl;
	//redefine minheap, since priority_queue doesn't have clear() in c++
	priority_queue<GridCell> minheap;
	// Fix the flats
	//some of the potential outlets have lower cells
	while (!possible_outlets.empty()) {
		auto cell = possible_outlets.back();
		possible_outlets.pop_back();
		double z = elevation_flats_filled_m_2D[cell.first][cell.second];
		bool flag = false;

		for (int n = 0; n < 8; n++) {
			int rn = cell.first + dy[n];
			int cn = cell.second + dx[n];

			if (checkfolder(rn, cn)) {
				double zn = elevation_flats_filled_m_2D[rn][cn];

				if (zn < z && zn != NODATA_code) {
					flag = true;
					break;
				}
			}
		}

		if (flag) {
			//it's confirmed as an outlet
			minheap.push({ cell.first, cell.second, z });
		}
	}

	int num_outlets = minheap.size();
	vector<GridCell> outlets;

	while (!minheap.empty()) {
		auto cell = minheap.top();
		minheap.pop();

		if (flats[cell.row][cell.column] != 3) {
			double z = elevation_flats_filled_m_2D[cell.row][cell.column];
			flats[cell.row][cell.column] = 3;

			if (!outlets.empty()) {
				outlets.clear();
			}

			outlets.push_back(cell);
			//are there any other outlet cells at the same elevation (likely for the same feathre)
			bool flag = true;
			while (flag) {
				if (!minheap.empty()) {
					auto cell2 = minheap.top();
					if (cell2.priority == z) {
						flats[cell2.row][cell2.column] = 3;
						outlets.push_back(cell2);
						minheap.pop();
					}
					else {
						flag = false;
					}
				}
				else {
					flag = false;
				}
			}

			priority_queue<GridCell2> minheap2;
			
			for (const auto& cell2 : outlets) {
				double z = elevation_flats_filled_m_2D[cell2.row][cell2.column];
				for (int n = 0; n < 8; n++) {
					int rn = cell2.row + dy[n];
					int cn = cell2.column + dx[n];

					if (checkfolder(rn, cn)) {
						if (flats[rn][cn] != 3) {
							double zn = elevation_flats_filled_m_2D[rn][cn];

							if (zn == z && zn != NODATA_code) {
								minheap2.push({ rn, cn, z, elevation_working_m_2D[rn][cn] });
								elevation_flats_filled_m_2D[rn][cn] = z + small_num;
								//cout << "Flat found, row: " << rn + 1 << ", col: " << cn + 1 << ", increase elevation to: " << z + small_num << endl;
								flats[rn][cn] = 3;
							}
						}
					}	
				}
			}

			while (!minheap2.empty()) {
				auto cell2 = minheap2.top();
				minheap2.pop();
				double z = elevation_flats_filled_m_2D[cell2.row][cell2.column];

				for (int n = 0; n < 8; n++) {
					int rn = cell2.row + dy[n];
					int cn = cell2.column + dx[n];

					if (checkfolder(rn, cn)) {
						if (flats[rn][cn] != 3) {
							double zn = elevation_flats_filled_m_2D[rn][cn];

							if (zn < z + small_num && zn >= cell2.z && fabs(zn - elevation_depression_filled_m_2D[rn][cn]) > 1e-6) {
								minheap2.push({ rn, cn, cell2.z, elevation_working_m_2D[rn][cn] });
								elevation_flats_filled_m_2D[rn][cn] = z + small_num;
								//cout << "Flat found, row: " << rn + 1 << ", col: " << cn + 1 << ", increase elevation to: " << z + small_num << endl;
								flats[rn][cn] = 3;
							}
						}
					}
				}
			}
		}
	}
}
//Note: Consider refactor to simply convert Slope_Ground_Horn_rad_Map and Aspect_Ground_N_0_rad_Map to 1D vectors
//Ground_Slope_and_Aspect_File_Read function reads prepared topographic data
void TerrainProcessor::Ground_Slope_and_Aspect_File_Read()
{
	double temp;
	string space, line;

	//ifstream function readSlopeFile directory path and file name to read file 
	ifstream readSlopeFile(Directory_Input_CLArg + "z_SlopeGround_rad.asc");

	//if readSlopeFile function returns a False value, meaning it did not work, then
	if (!readSlopeFile) {
		cout << "Ground_Slope_and_Aspect_File_Read: z_SlopeGround_rad.asc could not be opened. Delete z_TIorganizer.asc and z_FDorganizer.asc and try again." << endl;
	}

	//getline readFDorganizer is reading in the 6 header lines of data from map file
	int xllcorner_m, yllcorner_m;
	getline(readSlopeFile, line);
	istringstream(line) >> space >> nCols;
	getline(readSlopeFile, line, '\n');
	istringstream(line) >> space >> nRows;
	getline(readSlopeFile, line, '\n');
	istringstream(line) >> space >> xllcorner_m;
	getline(readSlopeFile, line, '\n');
	istringstream(line) >> space >> yllcorner_m;
	getline(readSlopeFile, line, '\n');
	istringstream(line) >> space >> cellsize_m;
	getline(readSlopeFile, line, '\n');
	istringstream(line) >> space >> NODATA_code;

	//while loop to readSlopeFile all pixel data values, assigned to temp variable
	while (readSlopeFile >> temp) {
		//if pixel value is real data then
		if (temp != NODATA_code) {
			//slopeGround_1D_rad vector has temp value added
			slopeGround_1D_rad.push_back(temp);
		}
		//if pixel is NODATA then
		else {
			//slopeGround_1D_rad vector has NODATA_code value added
			slopeGround_1D_rad.push_back(NODATA_code);
		}
	}
	//close readSlopeFile file
	readSlopeFile.close();

	//ifstream function readAspectFile directory path and file name to read file 
	ifstream readAspectFile(Directory_Input_CLArg + "z_AspectGround_N_0_rad.asc");

	//if readAspectFile function returns a False value, meaning it did not work, then
	if (!readAspectFile) {
		cout << "Ground_Slope_and_Aspect_File_Read: z_AspectGround_N_0_rad.asc could not be opened. Delete z_TIorganizer.asc and z_FDorganizer.asc and try again." << endl;
	}

	//getline readFDorganizer is reading in the 6 header lines of data from map file
	getline(readAspectFile, line);
	istringstream(line) >> space >> nCols;
	getline(readAspectFile, line, '\n');
	istringstream(line) >> space >> nRows;
	getline(readAspectFile, line, '\n');
	istringstream(line) >> space >> xllcorner_m;
	getline(readAspectFile, line, '\n');
	istringstream(line) >> space >> yllcorner_m;
	getline(readAspectFile, line, '\n');
	istringstream(line) >> space >> cellsize_m;
	getline(readAspectFile, line, '\n');
	istringstream(line) >> space >> NODATA_code;

	//while loop to readAspectFile all pixel data values, assigned to temp variable
	while (readAspectFile >> temp) {
		//if pixel value is real data then
		if (temp != NODATA_code) {
			//aspectGround_1D_N_0_rad vector has temp value added
			aspectGround_1D_N_0_rad.push_back(temp);
		}
		//if pixel is NODATA_code then
		else {
			//aspectGround_1D_N_0_rad vector has NODATA_code value added
			aspectGround_1D_N_0_rad.push_back(NODATA_code);
		}
	}
	//close readAspectFile file
	readAspectFile.close();
}

//elevation function gets value of map at specific rows and column for current folder
double TerrainProcessor::getElevation(int row, int col)
{
	return elevation_m_2D[row][col];
}

//TerrainProcessor::clear() function clears all the vectors before they are used in the index calculations
void TerrainProcessor::clearVectorsRelatedToTopographicIndex()
{
	//clear variable vectors
	ATB.clear();
	Area.clear();
	TI_Values_In_HistogramBin.clear();
	TI_Area_frac.clear();
	TI_Value.clear();
}

//checkfolder function tests if row or col are beyond range of 0 to nRows and nCols
bool TerrainProcessor::checkfolder(int row, int col)
{
	if (row < 0 || row >(nRows - 1) || col < 0 || col >(nCols - 1)) {
		return false;
	}
	else {
		return true;
	}
}

//Note: Consider rrefactor to replace fileExists with EnsureFileExists or other singular function
//fileExists function tests if filename exists
bool TerrainProcessor::fileExists(const string& filename)
{
	struct stat buf;
	if (stat(filename.c_str(), &buf) != -1)
	{
		return true;
	}
	return false;
}
