#ifndef Inputs_H
#define Inputs_H
/*
* HydroPlus is source code and models developed and managed by Theodore Endreny and his students at SUNY ESF, te@esf.edu
* The US Forest Service and Davey Tree have provided strategic software development support through funding and personnel
* Attribution for use of software and code is requested. Please see Executive.cpp for documentation.
*/
#include <string>
#include <map>
//_USE_MATH_DEFINES gives access to Pi as M_PI, and must precede #include <cmath> or <math.h>
#define _USE_MATH_DEFINES
#include <cmath>
#include <math.h>
#include <fstream>
#include <sstream>
#include <unordered_set>
#include <numeric>
#include <iomanip>
#include <unordered_map>
#include <functional>
#include <stdexcept>
#include <utility>
#include <algorithm>
#include "../XMLReader/XML_Document.h"
#include "WeatherProcessor.h"
#include <cstring>
#include <span>   // C++20

//Forward declaration of class TerrainProcessor
class TerrainProcessor;

//Earth_Radius_m 6378137.0 m
#define Earth_Radius_m 6378137.0
//Earth_Eccentricity is 0.081819191; related to Flattening_WGS84_Ellipsoid = 1 / 298.257223563
#define Earth_Eccentricity 0.081819191

//ratio of degrees to radians
#define Ratio_Degree_to_Radian 180 / M_PI
//ratio of radians to degrees
#define Ratio_Radian_to_Degree M_PI / 180
//ratio of hour to seconds
#define Ratio_Hour_to_Second 1.0/3600.0
//ratio of seconds to hour
#define Ratio_Second_to_Hour 3600
//ratio of hour to day
#define Ratio_Hour_to_Day 24
//ratio of metric tonne to kg 
#define Ratio_Tonne_to_kg 1000
//Conversion between fraction and percent
#define Ratio_Decimal_to_Percent 0.01
//Define constant effective zero value for sums of water
#define Constant_1E_negative10 1.E-10     
//Define constant effective zero value for fractions of area
#define Constant_1E_negative5 1.E-5     
//ratio of kg to mg
#define Ratio_kg_to_mg 1000000
//ratio of mm*ha to L
#define Ratio_mmHa_to_L 10000
//Epsilon_Tolerance_1E_negative15 a tolerance-based threshold
//Note: calcFlowDirection_DInf requires Epsilon_Tolerance_1E_negative15 for peformance agreement between Linux and Windows builds
#define Epsilon_Tolerance_1E_negative15 1e-15

//MapMetadata Struct initialized for Albers or Lamberts or UTM
struct MapMetadata {
	double lowerLeftX_m;
	double lowerLeftY_m;
	double cellSize_m;
	int nCols;
	int nRows;
};

// ProjectionParams_Albers struct contains WKID 102008 North America Albers Equal Area Conic projection parameters
// Note: Source: Snyder, J.P., (1987). Map projections: A working manual. USGS Professional Paper 1395, https://pubs.usgs.gov/publication/pp1395
// Note: WKID 102008 is North America Albers Equal Area Conic projection, with predefined parameters:
// Note: ... Latitude of Origin (rad): 23∘ N = 23×(π/180) = 0.4014 radians
// Note: ... Central Meridian (rad): −96∘ W = −96×(π/180) = −1.6755 radians
// Note: ... Standard Parallel 1 (rad​): 29.5∘ N = 29.5×(π/180) = 0.5149 radians
// Note: ... Standard Parallel 2 (rad​): 45.5∘ N = 45.5×(π/180) = 0.7941 radians
// 
//Li- this is the parameters for NAD 1983 Contiguous USA Albers WKID 102965/5070
struct ProjectionParams_Albers {
	double phi0;	// phi0 = Latitude of origin (radians)
	double lambda0; // lambda0 = Central meridian (radians) 
	double phi1;	// phi1 = Standard parallel 1 (radians)
	double phi2;	// phi2 = Standard parallel 2 (radians)
	double falseEasting;
	double falseNorthing;
};

struct ProjectionParams_Lambert {
	double phi0;      // Latitude of Origin (in radians)
	double lambda0;   // Central Meridian (in radians)
	double phi1;      // Standard Parallel 1 (in radians)
	double phi2;      // Standard Parallel 2 (in radians)
	double falseEasting;
	double falseNorthing;
};

// ProjectionParams_UTM struct contains WKID 32600 and 32700 UTM northern and southern hemisphere projection parameters
// Note: Source: Snyder, J.P., (1987). Map projections: A working manual. USGS Professional Paper 1395, https://pubs.usgs.gov/publication/pp1395
// Note: WKID 102008 is North America Albers Equal Area Conic projection, with predefined parameters:
struct ProjectionParams_UTM {
	int zone;// UTM zone number (1–60)
	double centralMeridian;// Central meridian of the UTM zone (in radians)
	double falseEasting;// False easting (usually 500,000 meters)
	double falseNorthing;// False northing (0 meters for north, 10,000,000 meters for south
	double scaleFactor;// Scale factor at the central meridian (0.9996)
};

//EnsureTrailingBackslash function to ensure it has no parenthesis and ends with backslash
string EnsureTrailingBackslash(const string& directoryPath);
bool EnsureFileExists(const string& filePath);

// CompactRagged
//   A single flat (contiguous) array for all values in the model.
//	 Storage is organized as:
//   slot = (drawer, folder) pair
//   index = slot * Nvars + variable
//   This way, all values are stored in one contiguous array but
//   access is still O(1) by drawer/folder/var.
//==================================================================
// A helper to map variable names (strings) to integer indices
struct NameTable {
	vector<string> names;                 // ordered list of variable names
	unordered_map<string, int> to_index;   // lookup table: name -> column index

	// Add a new variable if missing, return its index
	int add(const string& key) {
		auto it = to_index.find(key);
		if (it != to_index.end()) return it->second; // already exists
		int id = (int)names.size();                 // new index = current size
		names.push_back(key);                       // record name in list
		to_index.emplace(names.back(), id);         // record mapping in hash table
		return id;
	}
	int size() const { return (int)names.size(); }
};

struct CompactRagged {
	NameTable varNames;     // registry for variable names
	// ----------------- Drawer/Folder dimensions -----------------
	int Ndrawers = 0;                       // number of drawers (pixels or statistical drawers)
	vector<int> folders_per_drawer;    // how many folders belong to each drawer
	vector<size_t> drawer_offset;      // prefix sum array: offset of first slot for drawer d
	// Example:
	//   folders_per_drawer = [2,3,1]: Drawer 0 has 2 folders, Drawer 1 has 3 folders, Drawer 2 has 1 folder
	//   drawer_offset      = [0,2,5,6]
	//   Drawer 0 starts at index 0, Drawer 1 starts after Drawer 0’s 2 folders → 0 + 2 = 2, Drawer 2 starts after Drawer 1’s 3 folders → 2 + 3 = 5, The final value (6) is a sentinel that equals the total number of folders (2 + 3 + 1).
	size_t Nslots = 0;                      // total number of (drawer,folder) slots = sum(folders_per_drawer)

	// ----------------- Variables and data -----------------
	int Nvars = 0;                          // number of variable columns
	vector<double> data;               // contiguous buffer, size = Nslots * Nvars
	// Layout: slot-major, then variable
	//   data[ slot * Nvars + v ]  => value of variable v for (drawer,folder)

	// ======================================================
	// Initialize with per-drawer folder counts
	//   per_drawer_folders[d] = number of folders in drawer d
	//   nvars = number of variables to pre-allocate
	//   init_val = fill value for all entries
	// ======================================================
	void init_ragged(const vector<int>& per_drawer_folders, int nvars, double init_val)
	{
		folders_per_drawer = per_drawer_folders;
		Ndrawers = (int)folders_per_drawer.size();
		// Build prefix sum array so we can map (drawer,folder) -> slot
		drawer_offset.resize(Ndrawers + 1);
		drawer_offset[0] = 0;
		for (int d = 0; d < Ndrawers; ++d) {
			drawer_offset[d + 1] = drawer_offset[d] + (size_t)folders_per_drawer[d];
		}
		Nslots = drawer_offset[Ndrawers];   // total slots = last prefix value

		// Initialize variable name table
		varNames.names.clear(); varNames.to_index.clear();
		varNames.names.reserve(nvars);
		varNames.to_index.reserve((size_t)nvars * 2);
		// Pre-populate generic variable names (var_0001, var_0002, ...)
		for (int v = 0; v < nvars; ++v) {
			ostringstream k; k << "var_" << setw(4)
				<< setfill('0') << (v + 1);
			varNames.add(k.str());
		}
		Nvars = nvars;

		// Allocate data buffer: Nslots rows × Nvars columns
		data.assign(Nslots * (size_t)Nvars, init_val);
	}

	// Uniform init: every drawer has exactly `folders` folders
	void init_uniform(int n_drawers, int folders, int nvars, double init_val) {
		vector<int> per_drawer(n_drawers, folders);
		init_ragged(per_drawer, nvars, init_val);   // calls your current ragged init
	}
	// ======================================================
	// Indexing helpers
	// ======================================================

	// Map (drawer,folder) -> flat slot index
	inline size_t slot(int drawer, int folder) const {
		// drawer_offset[d] tells you the starting slot for drawer d
		// Add folder index to jump inside this drawer’s range
		return drawer_offset[drawer] + (size_t)folder;
	}

	// Map (drawer,folder,var) -> final index in data[]
	inline size_t idx(int drawer, int folder, int v) const {
		return slot(drawer, folder) * (size_t)Nvars + (size_t)v;
	}

	// Direct reference by integer variable index
	inline double& by_index(int drawer, int folder, int v) {
		return data[idx(drawer, folder, v)];
	}

	// Direct reference by variable name
	inline double& by_key(int drawer, int folder, const string& name) {
		auto it = varNames.to_index.find(name);
		int v = (it == varNames.to_index.end())
			? add_var(name, 0)   // reallocates once, returns new index
			: it->second;
		return by_index(drawer, folder, v);
	}

	// ======================================================
	// Add new variables (columns)
	// ======================================================

	// Add one new variable column; return its index.
	int add_var(const string& name, double init_val = 0.0) {
		auto it = varNames.to_index.find(name);
		if (it != varNames.to_index.end()) return it->second; // already exists

		const int oldN = Nvars;
		const int newIdx = varNames.add(name);
		const int newN = oldN + 1;

		// Create a new buffer with one extra column per slot
		vector<double> newData;
		newData.resize(Nslots * (size_t)newN);

		if (oldN > 0) {
			// For each slot, copy old columns and set new column = init_val
			for (size_t s = 0; s < Nslots; ++s) {
				const size_t baseOld = s * (size_t)oldN;
				const size_t baseNew = s * (size_t)newN;
				memcpy(&newData[baseNew], &data[baseOld], oldN * sizeof(double));
				newData[baseNew + (size_t)newIdx] = init_val;
			}
		}
		else {
			// First variable: just fill the new column with init_val
			for (size_t s = 0; s < Nslots; ++s) {
				newData[s] = init_val;
			}
		}

		// Swap buffers and update count
		data.swap(newData);
		Nvars = newN;
		return newIdx;
	}

	// Does beC contain a column with this name?
	inline bool has_var(const string& varName) const {
		return varNames.to_index.find(varName) != varNames.to_index.end();
	}

	// --- Optional dimension helpers (read-only) ---
	// for debugging, usage: beC -> print_var_dims(beC, "Storage_RainSnow_TreeCanopy_m");
	inline int n_drawers(const CompactRagged* beC) {
		return static_cast<int>(beC->drawer_offset.size()) - 1;  // assuming sentinel at end
	}
	inline int folders_in_drawer(const CompactRagged* beC, int d) {
		return static_cast<int>(beC->drawer_offset[d + 1] - beC->drawer_offset[d]);
	}
	inline size_t total_slots(const CompactRagged* beC) { return beC->Nslots; }

	inline bool in_range_drawer(const CompactRagged* beC, int d) {
		return d >= 0 && d + 1 < static_cast<int>(beC->drawer_offset.size());
	}
	inline bool in_range_folder(const CompactRagged* beC, int d, int f) {
		return in_range_drawer(beC, d) && f >= 0 && f < folders_in_drawer(beC, d);
	}

	// Print beC shape + a variable’s presence/index
	inline void print_var_dims(const CompactRagged* beC, const string& var) {
		cout << "beC dims: drawers=" << n_drawers(beC)
			<< ", slots=" << total_slots(beC)
			<< ", Nvars=" << beC->Nvars << "\n";

		/*cout << "folders per drawer: ";
		for (int d = 0; d < n_drawers(beC); ++d) {
			if (d) cout << ", ";
			cout << folders_in_drawer(beC, d);
		}
		cout << "\n";*/

		if (beC->has_var(var)) {
			int v = beC->varNames.to_index.at(var);
			cout << "var \"" << var << "\" exists, index=" << v << "\n";
		}
		else {
			cout << "var \"" << var << "\" NOT registered\n";
		}
	}

	inline void print_vec_dims(const CompactRagged* beC, const string& vec_name)
	{
		auto it = beC->vec_name_to_idx.find(vec_name);
		if (it == beC->vec_name_to_idx.end()) {
			cout << "vec \"" << vec_name << "\" NOT found.\n";
			return;
		}

		const auto& col = beC->vec_cols[it->second];
		cout << "vec \"" << vec_name << "\" exists.\n";
		cout << "  total slots: " << beC->Nslots
			<< ", total values: " << col.values.size() << "\n";

		// Print the first few slot lengths
		for (size_t s = 0; s < min<size_t>(beC->Nslots, 5); ++s) {
			uint32_t a = col.offsets[s];
			uint32_t b = col.offsets[s + 1];
			uint32_t len = b - a;

			cout << "  slot " << s << " length = " << len << " values = [";

			if (len == 0) {
				cout << "]\n";
				continue;
			}

			// Print a few values (cap to avoid spamming)
			uint32_t printCount = min<uint32_t>(len, 10); // print up to 10 values per slot
			for (uint32_t i = 0; i < printCount; ++i) {
				cout << col.values[a + i];
				if (i + 1 < printCount) cout << ", ";
			}
			if (len > printCount) cout << ", ...";
			cout << "]\n";
		}
	}
	// ----------------------------------------------------
	//  Ragged (vector-per-slot) support
	// ----------------------------------------------------
	struct RaggedColumn {
		vector<double>   values;    // concatenated data for all slots
		vector<uint32_t> offsets;   // length = Nslots + 1; slot s => [offsets[s], offsets[s+1])
	};

	unordered_map<string, uint32_t> vec_name_to_idx; // variable name → column index
	vector<RaggedColumn>                 vec_cols;         // actual column storage

	inline bool has_vec(const string& name) const {
		return vec_name_to_idx.find(name) != vec_name_to_idx.end();
	}

	// Ensure a vector column exists; initially all slot vectors are empty
	uint32_t ensure_vec(const string& name) {
		auto it = vec_name_to_idx.find(name);
		if (it != vec_name_to_idx.end()) return it->second;

		uint32_t idx = (uint32_t)vec_cols.size();
		vec_name_to_idx[name] = idx;
		RaggedColumn rc;
		rc.offsets.assign(Nslots + 1, 0);
		vec_cols.push_back(move(rc));
		return idx;
	}

	// Mutable span (content only, no length change)
	//modify existing content only. Get a reference (span) to an already-existing vector slice for a given (drawer, folder, var).
	inline span<double> vec_by_key(int drawer, int folder, const string& name) {
		auto it = vec_name_to_idx.find(name);
		if (it == vec_name_to_idx.end()) it = vec_name_to_idx.emplace(name, ensure_vec(name)).first;
		RaggedColumn& col = vec_cols[it->second];
		const size_t s = slot(drawer, folder);
		uint32_t a = col.offsets[s], b = col.offsets[s + 1];
		return span<double>(col.values.data() + a, b - a);
	}

	
	// Resize the vector for a specific (drawer, folder, variable) slot. This changes only ONE slot's vector length, preserving existing values.
	// ... Expanding adds new elements initialized to fill_val. Shrinking removes tail elements.
	inline void vec_resize_slot(int drawer, int folder, const string& name, uint32_t new_length, double fill_val = 0.0)
	{	
		// Ensure the vector-column exists; get its column index
		uint32_t col_idx = ensure_vec(name);
		RaggedColumn& col = vec_cols[col_idx];
		// Convert (drawer, folder) → slot index
		const size_t slot_index = slot(drawer, folder);
		// Current slice range for this slot in the flattened `values` array
		uint32_t offset_begin = col.offsets[slot_index], offset_end = col.offsets[slot_index + 1];
		uint32_t current_length = offset_end - offset_begin;
		// No change needed
		if (current_length == new_length) return;
		// EXPAND slot vector
		if (current_length < new_length) {
			// Expand: insert extra values
			uint32_t add_count = new_length - current_length;
			col.values.insert(col.values.begin() + offset_end, add_count, fill_val);
			// All following offsets must be increased
			for (size_t i = slot_index + 1; i < col.offsets.size(); ++i)
				col.offsets[i] += add_count;
		}
		//SHRINK slot vector
		else {
			// Shrink: erase tail
			uint32_t remove_count = current_length - new_length;
			col.values.erase(col.values.begin() + (offset_begin + new_length), col.values.begin() + offset_end);

			for (size_t i = slot_index + 1; i < col.offsets.size(); ++i)
				col.offsets[i] -= remove_count;
		}
	}

	// Ensure a vector for this slot has *at least* want_length elements. If expansion is needed, new values use fill_val. Returns a fresh span to the updated slot.
	inline span<double> vec_ensure_length(int drawer, int folder, const string& name, uint32_t want_length, double fill_val = 0.0)
	{
		// Get current vector slice (creates empty one if missing)
		auto slice = vec_by_key(drawer, folder, name);
		if (slice.size() < want_length) {
			vec_resize_slot(drawer, folder, name, want_length, fill_val);
			slice = vec_by_key(drawer, folder, name);
		}
		return slice;
	}

	// ----------------------------------------------------
	// String variable storage (parallel to numeric)
	// ----------------------------------------------------
	int Nvars_str = 0;
	unordered_map<string, vector<string>> data_str;  // name → string values

	// Add or access string variable
	inline string& by_key_str(int drawer, int folder, const string& name) {
		auto it = data_str.find(name);
		size_t s = slot(drawer, folder);

		// If variable not exist, allocate it for all slots with empty string
		if (it == data_str.end()) {
			data_str[name] = vector<string>(Nslots, "");
			++Nvars_str;
		}

		return data_str[name][s];
	}

	// Ensure a string column exists; fill all slots with init
	inline void ensure_var_str(const string& name, const string& init = "") {
		if (!has_var_str(name)) {
			data_str[name] = vector<string>(Nslots, init);
			++Nvars_str;
		}
	}

	// Check if string variable exists
	inline bool has_var_str(const string& varName) const {
		return data_str.find(varName) != data_str.end();
	}
	// ============================================================================
	// NOTE ON ACCESS VIOLATION BUG (beC::by_key reallocation issue)
	// ----------------------------------------------------------------------------
	// When a variable name does not exist in beC, calling by_key() auto-triggers add_var(), which reallocates the entire internal data buffer. 
	// If multiple by_key() calls appear in the same expression (e.g. a + b), and one of them causes reallocation, the others may access invalidated references or memory, leading to an access violation.
	//
	// Example (original problematic line):
	//     beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_RainSnow_TreeCanopy_m") = beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") + beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Snow_TreeCanopy_m");
	//
	// Fix / Prevention:
	// 1. Always ensure all variables used by by_key() are registered before use. For safety, if a variable might not exist, assign a default value (e.g. 0.0).
	//		e.g. if (!beC->has_var("Storage_Rain_TreeCanopy_m")) beC->by_key(MapPixel_ID, DataFolder_ID, "Storage_Rain_TreeCanopy_m") = 0.0;
	//
	// 2. Calling by_key() in separate statements or switching evaluation order to avoid simultaneous reallocation.
	//
	// (Documented by Li Zhang, Nov 2025)
	// ============================================================================

};
//==================================================================

class Inputs
{
public:

	Inputs(string &Directory_Input_CLArg);
	CompactRagged beC;   // main (grid or statistical drawers)
	CompactRagged beBG;  // 1 slot (folder 0) per unique Block Group
	vector<int> Unique_BlockGroup_IDs; // BG IDs in drawer order
	unordered_map<int, int> bg_to_drawer; // BG ID -> drawer index [0..Unique_BlockGroup_IDs.size()-1]
	// Build once after blockgroup.asc has been read and nCols/nRows are known
	void buildBlockGroupIndex(double NODATA_code);

	// Convenience accessors
	inline double& by_BlockGroup(int blockGroupID, const string& varName) {
		int d = bg_to_drawer.at(blockGroupID);    // throws out_of_range if unknown BG
		return beBG.by_key(/*drawer=*/d, /*folder=*/0, varName);
	}
	inline bool has_BlockGroup(int blockGroupID) const {
		return bg_to_drawer.find(blockGroupID) != bg_to_drawer.end();
	}
	

	void ProcessTopographyData();
	static void InitializeLAIVariables();
	static void CalcTreeLAI(int JulianDay);
	static void CalcShortVegLAI(int JulianDay);
	//Date_YYYYMMDD_to_JulianDay function calculates the Julian day based on inputs of year, month, and days
	static int Date_YYYYMMDD_to_JulianDay(int Date_YYYYMMDD);
	//Month_Year_to_Days_in_Month function calculates the days in a month given inputs of year and month
	static int Month_Year_to_Days_in_Month(int Year_YYYY, int Month_MM);
	//ExtractGregorianDateParts function extracts YYYY MM DD and if available HH into an array
	static void ExtractGregorianDateParts(int Date_YYYYMMDDHH, int DateParts_Gregorian[]);
	//GregorianToJulianDate function converts Gregorian date YYYYMMDD into an integer number based on Julian Date System
	static int GregorianToJulianDate(int DateParts_Gregorian[]);
	//JulianDateToWeekday function converts Julian Date System integer into weekday integer, ranging from 0=Sunday to 6=Saturday
	static int JulianDateToWeekday(int Date_JulianDate_base_4713BCE_int);
	//get reversed flow direction 
	static double D8_flowdir_reverse(double flowdir);
	//get manning n and Hydraulic radius
	static double getManningn(double nlcd_categories);
	static double getHydraulicradius(double nlcd_categories);
	static double get_P_release_mean(double nlcd_categories);
	//convert the code of fdr from esri to whitebox which is using in HydroPlus
	static double flowdircode_esri_to_whitebox(double flowdir);
	static double flowdircode_whitebox_to_esri(double flowdir);
	//calculates the mean value of a vector while excluding a specified NODATA value
	static double calculateMean(const vector<double>& data);
	//isReferenceStationPixel determines if MapPixel_ID is reference station
	bool isReferenceStationPixel(int MapPixel_ID) const;
	int getReferenceStationIndex(int MapPixel_ID);

	//Inputs::ProcessTemperatureExecutionParams function reads HydroPlusConfig.xml inputs for Spatial Buffer with GI
	void ProcessSpatialBufferGIExecutionParams();
	//update bulk area for each GI device
	double BufferGIContributingAreaCalc(int GI_Index);
	double BufferGIAreaCalc(int GI_Index);
	//Check_Map_Size_and_NoData overloaded for map vectors that are double vs double
	void Check_Map_Size_and_NoData(const vector<double>& map1, const vector<double>& map2, const string& map1Name, const string& map2Name);
	//Check_Map_Size_and_NoData overloaded for map vectors that are double vs int
	void Check_Map_Size_and_NoData(const vector<double>& map1, const vector<int>& map2, const string& map1Name, const string& map2Name);
	void ProcessTemperatureExecutionParams();
	//Process_Scalar_TemperatureExecutionParams to read HydroPlusConfig.xml scalars of map values from TemperatureExecutionParams
	void Process_Scalar_TemperatureExecutionParams(const string& TemperatureExecutionParams_str);
	//Process_Vector_TemperatureExecutionParams to read HydroPlusConfig.xml vectors of anthropogenic heat and degree hour from TemperatureExecutionParams
	void Process_Vector_TemperatureExecutionParams(const string& input, vector<double>& targetVector);
	//Convert_String_to_Lower_Case uses pointer to string that is modified
	static string Convert_String_to_Lower_Case(const string& string_Unknown_Case);
	//validateRatioVector to validate length of vectors for CoolAir Anthropogenic Heat Flux modifications
	void validateRatioVector(vector<double>& RatioVector, int expected_size, const string& ratio_name, const string& default_csv);
	//IsNearUrban will search the neighbors around a pixel to determine if they are an urban class
	bool IsNearUrban(int mapPixel_ID, int searchRadius = 1); 
	//IsAppropriateLandCoverForScenario will return true if MapPixel_ID is appropriate for scenario, e.g., irrigation, no transpiration
	bool IsAppropriateLandCoverForScenario(int MapPixel_ID, int Scenario_CoolAir_LandCover_Class);
	//IsAppropriateDay_or_NightForScenario will return true if ElevationAngle_Solar_rad is appropriate for scenario, e.g., day or night
	bool IsAppropriateDay_or_NightForScenario(int timeStep, double ElevationAngle_Solar_rad, int Scenario_CoolAir_Irrigate_Day_or_Night);
	//IrrigationCountByDay map contains the MapPixel_ID, CurrentDate, LandCover_Type and then the count
	map<int, map<int, map<string, int>>> IrrigationCountByDay;
	//IsAppropriateTimesPerDay will return true if IrrigationCountByDay <= Scenario_CoolAir_Irrigate_MaxPerDay, or Scenario_CoolAir_Irrigate_MaxPerDay is 0
	bool IsAppropriateTimesPerDay(int MapPixel_ID, int currentDate, const string& landCoverType, int Scenario_CoolAir_Irrigate_MaxPerDay);
	//IncreaseTimesPerDay will increase the count once an action such as irrigation took place
	void IncreaseTimesPerDay(int MapPixel_ID, int currentDate, const string& landCoverType);
	//CleanTimesPerDayMap function will clear map if current date is not the map date for maps such as IrrigationCountByDay
	void CleanTimesPerDayMap(map<int, map<int, map<string, int>>>& IrrigationCountByDay, int currentDate);
	//Check_Map_Header called to compare with dem.asc header
	static void Check_Map_Header(string map_name, int nCols, int nRows, double NODATA_code, double Length_Pixel_Side_m, double xllcorner, double yllcorner);
	//SafeDivide function avoids division by 0 and returns 0 rather than nan; relies on .h version containing Eq_defaultValue = 0.0
	static double SafeDivide(double Eq_numerator, double Eq_denominator, double Eq_defaultValue = 0.0);
	//roundToDecimalPlaces rounds input to 3 decimal places
	static double roundToDecimalPlaces(double value, int decimal_places);

	int get_Q_count() { return (totalDischargeCount - 1); }
	int getcols() { return nDemDataCols; }
	int getrows() { return nDemDataRows; }
	double getNODATA() { return NODATA_code; }
	double getxllcorner() { return nDemDataxllcorner; }
	double getyllcorner() { return nDemDatayllcorner; }

	void setcols(int x) { nDemDataCols = x; }
	void setrows(int x) { nDemDataRows = x; }
	void setxllcorner(double x) { nDemDataxllcorner = x; }
	void setyllcorner(double x) { nDemDatayllcorner = x; }

	int getCount() { return (totalCount); }

	//vector<vector<pair<double, double>>> Latitude_Longitude_vectorPair_dd vector created to store latitude and longitude as decimal degree
	vector<vector<pair<double, double>>> Latitude_Longitude_vectorPair_dd;
	//Flag_LatitudeLongitude_Computed flag initialized to false, to track if Latitude_Longitude_vectorPair_dd computed
	bool Flag_LatitudeLongitude_Computed = false;

	string sanitizeLineEndings(const string& variable_string);

	static bool isGreaterThanOrAlmostEqual(double a, double b, double eps = Epsilon_Tolerance_1E_negative15);
	static bool isLessThanOrAlmostEqual(double a, double b, double eps = Epsilon_Tolerance_1E_negative15);
	static bool isGreaterThan(double Variable_a, double Variable_b, double eps = Epsilon_Tolerance_1E_negative15);
	static bool isLessThan(double Variable_a, double Variable_b, double eps = Epsilon_Tolerance_1E_negative15);
	static double forceZeroIfLessThanOrAlmostEqualZero(double value, double eps = Epsilon_Tolerance_1E_negative15);
	static size_t findIndex_Meteorological_Min_or_Max(const vector<double>& Meteo_values, string Max_or_Min);
	static size_t findIndex_Min_or_Max_in_Window(const vector<double>& Meteo_values, size_t Index_Vector, string Max_or_Min, int window_hours, double timestep_hours = 1.0);

	//HydroPlusConfig.xml parameters
	static unordered_map<string, string> SimulationStringParams;
	static unordered_map<string, double> SimulationNumericalParams;
	static unordered_map<string, double> SimulationScenarios;
	static unordered_map<string, double> VariableInputMaps;
	unordered_map<string, double> LandCoverParams;
	//map<string, double> VegetationParams;
	//map<string, double> SoilParams;
	static unordered_map<string, double> SimulationLocationParams;

	static unordered_map<string, double> InputXml;
	unordered_map<string, string> InputXml_string;
	//InputXml_StandardFolder and InputXmlString_StandardFolder are default read from HydroPlusConfig.xml for all folders
	static unordered_map<string, double> InputXml_StandardFolder;
	unordered_map<string, string> InputXmlString_StandardFolder;
	//InputXml_ReferenceStationFolder and InputXmlString_ReferenceStationFolder are optionally read from HydroPlusConfig.xml for reference folder
	static unordered_map<string, double> InputXml_ReferenceStationFolder;
	unordered_map<string, string> InputXmlString_ReferenceStationFolder;
	//Flag_InputXml_ReferenceStationFolder is default false, and changed to true if InputXml_ReferenceStationFolder exists
	bool Flag_InputXml_ReferenceStationFolder = false;

	vector<int> FolderType_List;
	vector<vector<unordered_map<string, double>>> DataDrawers;
	vector<vector<unordered_map<string, string>>> StringDataDrawers;

	unordered_map<string, double> TemperatureCalculationParams;
	unordered_map<string, double> TemperatureCalculationParams_ReferenceStation;
	bool Flag_TemperatureCalculationParams_ReferenceStation_Exists = false;
	unordered_map<string, string> TemperatureExecutionParams;

	unordered_map<string, double> CoolRiverVariables;
	unordered_map<string, string> CoolBuildingStringParams;
	unordered_map<string, double> CoolBuildingNumericalParams;
	unordered_map<string, double> BufferTemporalLoadingNumericalParams;
	unordered_map<string, double> BufferSpatialWeightingNumericalParams;
	static unordered_map<string, string> BufferSpatialWeightingStringParams;
	unordered_map<string, string>BufferNutrientsStringParams;

	//Time parameters for internal weather and radiation related calculations, based on standard time
	vector<string> SimulationTime_HMS;
	vector<int> SimulationDate_YYYYMMDD;
	vector<int> SimulationDate_GDH;
	//Time parameters for writing outputs to Daylight Saving or standard time
	vector<string> SimulationTime_Output_HMS;
	vector<int> SimulationDate_Output_GD;
	vector<int> SimulationDate_Output_GDH;
	
	//Vector of time or duration for each timeStep or time step
	vector<int> SimulationTimeStep_Duration_sec;
	int SimulationTimePeriod_seconds; 
	int SimulationTimePeriod_timeSteps;
	int SimulationDay_Duration_sec{};
	int SimulationDateStart_YYYYMMDDHH{};
	int SimulationDateStop_YYYYMMDDHH{};

	//vector<int> UserOutput_SimulationStartHour_HH, UserOutput_SimulationEndHour_HH;
	//AtmPres_kPa (kPa)
	vector<double> AtmPres_kPa;
	vector<double> WindSpd_mps;
	//WindDir_deg (degrees,N=0=360); used in EnergyLoad model
	vector<double> WindDir_deg;
	//Precip_mpdt (meters, per time step, dt; typically m/h)
	vector<double> Precip_mpdt;
	vector<double> RadiationNet_Wpm2;
	vector<double> Tair_C;
	vector<double> Tdew_C;
	vector<double> Rain_m;
	vector<double> Snow_m;
	double TotalRain_m;
	double TotalPrecip_m;
	double TotalSnow_m;
	//Defining vectors for meteorological inputs
	double Precipitation_Sum_m, Rain_Sum_m, Snow_Sum_m;
	vector<double> TemperatureAir_TS_F, TemperatureDewPoint_TS_F, NetRadiation_TS_Wpm2, WindSpeed_TS_mps, WindDirection_TS_deg, PressureAtm_TS_kPa, 
		Precipitation_TS_mpdt, Rain_TS_mpdt, Snow_TS_mpdt, SnowDepth_TS_m;
	vector<double> Runoff_Observed_m;
	//MonthlyAvgTair_C computed in CreateFile_EPW.cpp and used in CreateFile_IDF 
	vector<float> MonthlyAvgTair_C;

	//Flag_PrintBlockGroupDailyAndHourly_str to control use of Block Group data
	string Flag_PrintBlockGroupDailyAndHourly_str = "0";

	//Evaporation.csv parameters
	bool Flag_ComputePET;
	string evapXMLVal;
	vector<double> PotentialEvaporation_WaterOnGround_m;
	vector<double> PotentialEvaporation_WaterOnTree_m;
	vector<double> PotentialEvapotranspiration_TreeCover_m;
	vector<double> PotentialSublimation_SnowOnGround_m;
	vector<double> PotentialSublimation_SnowOnTree_m;
	vector<double> PotentialEvaporation_Ground_m;
	vector<double> PotentialEvaporation_Tree_m;
	vector<double> PotentialEvapotranspiration_m;
	vector<double> PotentialSublimation_Ground_m;
	vector<double> PotentialSublimation_Tree_m;

	//Radiation.csv parameters
	vector<double> Radiation_Shortwave_Direct_Wpm2;  
	vector<double> Radiation_Shortwave_Diffuse_Wpm2; 
	vector<double> Radiation_Longwave_Downwelling_Wpm2;
	vector<double> Radiation_Longwave_Upwelling_Wpm2;
	double Latitude_rad;
	vector<double> HourAngle_Solar_rad;
	vector<double> DeclinationAngle_Solar_rad;

	//Distributed cover maps: landcover.asc, imperviouscover.asc, treecover.asc, dem.asc, blockgroup.asc, AH_flux_Qtot_avg_Wpm2.asc
	//Note: File extension .asc is recognized by ArcGIS and QGIS as a standard input type when using the 6 header rows
	static vector<double> LandCover_NLCD_Class;
	static vector<double> ImperviousCover_frac, PerCoverdata;
	vector<double> TreeCover_frac;
	vector<double> TreeCover_base_frac;
	vector<double> TotalMapCover_frac;
	vector<double> SoilCoverNoTC_frac;
	vector<double> WaterCoverNoTC_frac;
	vector<double> ImperviousCoverNoTreeCover_frac;
	vector<double> ShortVegCoverNoTC_frac;
	vector<double> ConnectedImpervious;
	vector<double> TreeEvergreen_frac;
	vector<double> ShortVegEvergreen_frac;
	vector<double> TreeCoverOnPerviousCover_frac;
	vector<double> TreeCoverOnImperviousCover_frac;
	vector<double> PermeablePavementCover_frac;

	static vector<double> Waters01, Waters0;
	static vector<double> NHDfdr, NHDfdr_reverse, DEMfdr, negDEMfdr;
	static vector<double> s_trans;
	static vector<double>  s_Ksat;
	static double s_trans_avg, s_depth_avg;
	static vector<double> s_depth;
	static vector<double> s_wtable;
	static vector<double> GI_Area_m2, GI_Contributing_Area_m2;

	//AnthropogenicHeat_Flux_Qtot_Avg_Wpm2 is anthropogenic heat as total (Qtot) average energy (W/m2) from AH4GUC
	vector<double> AnthropogenicHeat_Flux_Qtot_Avg_Wpm2;
	string AnthropogenicHeat_Flux_Qtot_Avg_Wpm2_str;
	//AnthropogenicHeat_Flux_Qcr_Avg_Wpm2 is anthropogenic heat as commercial residential (Qcr) average energy (W/m2) from AH4GUC
	vector<double> AnthropogenicHeat_Flux_Qcr_Avg_Wpm2;
	string AnthropogenicHeat_Flux_Qcr_Avg_Wpm2_str;
	//AnthropogenicHeat_Flux_Qncr_Avg_Wpm2 is anthropogenic heat as non-commercial-residential (Qncr) average energy (W/m2) from AH4GUC
	vector<double> AnthropogenicHeat_Flux_Qncr_Avg_Wpm2;
	string AnthropogenicHeat_Flux_Qncr_Avg_Wpm2_str;
	//AnthropogenicHeat_Flux_Qtot_Map_Sum_Wpm2, Qcr (commercial residential) and Qncr stores the sum values 
	double AnthropogenicHeat_Flux_Qtot_Map_Sum_Wpm2;
	double AnthropogenicHeat_Flux_Qcr_Map_Sum_Wpm2;
	double AnthropogenicHeat_Flux_Qncr_Map_Sum_Wpm2;
	//ImperviousCover_Sum_frac stores the sum value to divide into the AnthropogenicHeat_Flux_Qtot_Map_Sum_Wpm2 etc. for AH_Flux/IC_frac
	double ImperviousCover_Sum_frac;
	//AH_Flux_Hour_Ratio defined as vector of AHE average hour value to average annual value
	vector<double> AH_Flux_Hour_Ratio;
	//AH_Flux_Month_Ratio defined as vector of AHE average month value to average annual value
	vector<double> AH_Flux_Month_Ratio;

	//map_Irrigation_Area_frac is fraction of pixel that received irrigation within Flag_Scenario_CoolAir_Irrigate
	vector<double> map_Irrigation_Area_frac;
	//map_Permeable_Area_of_IC_frac is fraction of imperviouscover.asc (IC) that becomes PermeablePavementCover_frac
	vector<double> map_Permeable_Area_of_IC_frac;


	unordered_set<int> LocationIdsToPrintTempOutput;
	vector<int> StartDatesForRowColTempOutput;
	vector<int> EndDatesForRowColTempOutput;
	vector<int> StartDatesForTimeBasedTempOutput;
	vector<int> EndDatesForTimeBasedTempOutput;
	vector<int> StartingBlockIdsToPrintTempOutput;
	vector<int> EndingBlockIdsToPrintTempOutput;
	vector<int> BlockGroup_ID;

	string demXMLVal;
	double AveSMD;
	double Size_TI_Value_Vector;
	double Size_TI_Area_Vector;
	vector<double> TI_Area_frac;
	double TI_Value_Average;

	//Define TI_Value vector double
	vector<double> TI_Value;
	//Define DEM vector double
	vector<double> Elevation_DEM_m;
	vector<double> FlowDir;
	//Define SlopeGround_rad vector double
	vector<double> SlopeGround_rad;
	//Define AspectGround_N_0_rad vector double
	vector<double> AspectGround_N_0_rad;

	//Define RepoDict map double
	unordered_map<string, double> RepoDict;
	//Define RepoVecDict map vector
	unordered_map<string, vector<double>> RepoVecDict;

	//LAI calculations
	static unordered_map<string, double> LAI_BAI_map;
	static vector<double> LAI_BAI_Tree_m2_p_m2;
	static vector<double> LAI_BAI_SVeg_m2_p_m2;
	static vector<double> LAI_Tree_m2_p_m2;
	static vector<double> LAI_SVeg_m2_p_m2;

	double MinAllowableWindSpeed;
	//Define non-routed runoff (m) vectors, pervious, directly connected impervious, sub-surface, and total
	//PAQi_m (m) is pervious area runoff
	vector<double> PAQi_m;
	//DCIAQi_m (m) is directly connected impervious runoff, the impervious runoff not sent to pervious area
	vector<double> DCIAQi_m;
	//WAQi_m (m) is water area runoff
	vector<double> WAQi_m;
	//SurQi_m (m) is suurface area runoff
	vector<double> SurfQi_m;
	//SubQi_m (m) is subsurface pervious area runoff
	vector<double> SubQi_m;
	//Qi_m (m) is total runoff
	vector<double> Qi_m;

	//Define routed runoff (m) vectors, directly connected impervious, pervious, water, sub-surface, and total
	vector<double> WAQi_Routed_m;
	vector<double> PAQi_Routed_m;
	vector<double> DCIAQi_Routed_m;
	vector<double> SurfQi_Routed_m;
	vector<double> SubQi_Routed_m;
	vector<double> Qi_Routed_m;

	// For CoolRiver RiverTimeSeries parameters
	vector<int> DatesRiverTime;
	vector<string> RiverTimeSeries;

	vector<double> Sediment_T_C;
	vector<double> Cloudiness_F;
	vector<double> Humidity_Per;

	vector<vector<double>> SW_Q_cms;
	vector<vector<double>> SW_T_C;

	map<int, vector<double>> SW_Q_cms1;
	map<int, vector<double>> SW_T_C1;

	map<int, vector<double>> GW_Q_cms1;
	map<int, vector<double>> GW_T_C1;

	vector<vector<double>> GW_Q_cms;
	vector<vector<double>> GW_T_C;
	vector<int> number;
	vector<double> dischargeDist;
	vector<double> dischargeValue1;

	// For CoolRiver RiverStationSeries parameters
	vector<int> StationNumber;
	vector<int> StationDistanceBedData;
	vector<double> DepthMesurments;

	vector<string> TypeOfRiver;
	vector<double> HorizontalBedConnectivity;

	vector<double> BedParticleSize;
	vector<double> Embeddedness;

	vector<int> StationDistanceMorphology;
	vector<double> AreaR_m2;
	vector<double> WidthR_m;

	vector<double> DepthR_m;
	vector<double> SlopeR_f;

	vector<int> RowRiver;
	vector<int> ColRiver;
	vector<double> LongitudeR_frac;
	vector<double> LatitudeR_frac;
	vector<double> ZRiver_m;

	vector<int> StationDistanceShading;
	vector<double> EastBankHR_m;
	vector<double> EastTreeHR_m;
	vector<double> EastBuildingHR_m;
	vector<double> EastBankDistR_m;
	vector<double> EastCanDistR_m;
	vector<double> EastBuildingDistR_m;
	vector<double> EastBufferWR_m;

	vector<double> WestBBankHR_m;
	vector<double> WestBankHR_m;
	vector<double> WestTreeHR_m;
	vector<double> WestBuildingHR_m;
	vector<double> WestBankDistR_m;
	vector<double> WestCanDistR_m;
	vector<double> WestBuildingDistR_m;
	vector<double> WestBufferWR_m;

	vector<int> StationDistanceShadingPercent;
	vector<double> ElevationR_m;
	vector<double> StreamAzimuthR_deg;
	vector<double> ShadeFactorR_f;
	vector<double> ViewToSkyR_f;
	//HecRas.dat file
	vector<double> downStreamDist;
	vector<double> totalQ;
	vector<double> chElev;
	vector<double> waterSurElev;
	vector<double> velChnl;
	vector<double> xArea;
	vector<double> topWidth;
	vector<double> waterSurSlope;
	vector<double> wettedPerimeter;

	//unordered_map<string, int> MultipleStationObjectID_to_AttributeVector_Index_map;
	unordered_map<string, int> MultipleStationObjectID_to_OrganizerDataDrawersIndex_map;
	unordered_map<string, int> MultipleStationObjectID_to_MapPixelID_map;

	//CoolRiver variable
	int totalCount;
	//ECdynamic cell number
	double cell_count = 0.0, cell_count_Nowater = 0.0;
	vector<double> Cell_number;
	map <string, vector<double>> EC_vec_map;
	map <string, vector<double>> EMC_vec_map;
	map <string, vector<double>> water_volume_vec_map;
	vector<double> DischargeRecord;
	vector<int> DischargeMonth, DischargeDay, DischargeYear;
	double foldersize;
	//GI Water Quality
	//Pollutants.csv pollutant types TSS, BOD, .., Zn are read into vector Pollutant_Name[], used to ID HydroPlusConfig.xml removal efficiencies
	//Note: The HydroPlusConfig.xml elements TSS, BOD, .., Zn are removal efficiencies, units percent
	vector<string> Pollutant_Name;
	vector<double> Pollutant_ConcentrationMean_mg_p_L;
	vector<double> Pollutant_ConcentrationMedian_mg_p_L;

	//Initialize NLCD land cover types to zero for summing in Inputs::readNLCDLandCoverFile function
	int LC21 = 0, LC22 = 0, LC23 = 0, LC24 = 0, LC41 = 0, LC42 = 0, LC43 = 0, LC71 = 0, LC81 = 0, LC82 = 0;
	int LC31 = 0, LC52 = 0;
	//water area
	int LC11 = 0, LC90 = 0, LC95 = 0;

	int run_ECDynamicModel = 1;
	//Flag_Scenario_CoolAir_LandCover initialized to 0 to result in no change to land cover
	int Flag_Scenario_CoolAir_LandCover = 0;
	//Scenario_CoolAir_Base_TCI_to_TC_frac initialized to -1 to result in no change to land cover
	double Scenario_CoolAir_Base_TCI_to_TC_frac = -1;
	//Scenario_CoolAir_Alternative_TC_Min_frac initialized to 0.0 to result in no change to land cover
	double Scenario_CoolAir_Alternative_TC_Min_frac = 0.0;
	//Scenario_CoolAir_Alternative_TC_Max_frac initialized to 1.0 to result in no change to land cover
	double Scenario_CoolAir_Alternative_TC_Max_frac = 1.0;
	//Scenario_CoolAir_Alternative_IC_replacing_TC_frac initialized to 0.0 to result in no change to land cover
	double Scenario_CoolAir_Alternative_IC_replacing_TC_frac = 0.0;
	//Scenario_CoolAir_Alternative_IC_Max_frac initialized to 1.0 to result in no change to land cover
	double Scenario_CoolAir_Alternative_IC_Max_frac = 1.0;
	//Scenario_CoolAir_Alternative_IC_Min_frac initialized to 0.0 to result in no change to land cover
	double Scenario_CoolAir_Alternative_IC_Min_frac = 0.0;
	//Scenario_CoolAir_Alternative_TC_Change_frac initialized to 0.0 to result in no change to land cover
	double Scenario_CoolAir_Alternative_TC_Change_frac = 0.0;
	//Scenario_CoolAir_Alternative_IC_Change_frac initialized to 0.0 to result in no change to land cover
	double Scenario_CoolAir_Alternative_IC_Change_frac = 0.0;
	//Scenario_CoolAir_Alternative_Change_type initialized to 0 for Absolute, 1 for Relative
	//Note: Absolute adds a fixed fraction to the existing value, Relative adds a fraction of the existing value to the existing value
	int Scenario_CoolAir_Alternative_Change_type = 0;

	int row_ReferenceStation;
	int col_ReferenceStation;
	//int formulation used when reference station was always a single pixel inside map area
	int MapPixel_ReferenceStation_ID;
	vector<int> MapPixel_ReferenceStation_IDs;
	static int nCols;
	static int nRows;
	static double Length_Pixel_Side_m;
	static double cellsize_m;
	static double xllcorner_m;
	static double yllcorner_m;
	static double NODATA_code;

	static int nCols_temp;
	static int nRows_temp;
	static double cellsize_m_temp;
	static double xllcorner_m_temp;
	static double yllcorner_m_temp;
	static double NODATA_temp;

	vector<string> Station_OBJECTID_vec, Station_USAF_vec, Station_WBAN_vec, Station_Path_folder_vec, Station_TI_option_vec;
	vector<long double> Station_NLCD_Class_vec, Station_Elev_m_vec, Station_Slope_deg_vec, Station_Aspect_deg_vec, Station_TC_per_vec, Station_IC_per_vec, Station_SVeg_in_Gap_frac_vec, Station_Soil_per_vec, Station_Water_per_vec, Station_AH_Flux_Wpm2_vec;
	vector<long double>  Station_Easting_m_vec, Station_Northing_m_vec;
	vector<int> Station_BlockGroup_ID_vec;

	double conductivity_min, soildepth_min;
	vector<string> GI_list;
	vector<double> GI_ContributingArea_list;
	vector<int> GI_ulcorner_rows, GI_ulcorner_cols, GI_lrcorner_rows, GI_lrcorner_cols;
	static pair<double, double> Projection_LatitudeLongitude_to_Albers(double Latitude_MapPixel_dd, double Longitude_MapPixel_dd, const ProjectionParams_Albers& ProjectionParams_Albers);
	static pair<double, double> Projection_LatitudeLongitude_to_Lambert(double Latitude_MapPixel_dd, double Longitude_MapPixel_dd, const ProjectionParams_Lambert& ProjectionParams_Lambert);
	static pair<double, double>  Projection_LatitudeLongitude_to_UTM(double latitude_deg, double longitude_deg, const ProjectionParams_UTM& utmParams);
	// Declaration only (no initialization here non-const static members cannot be initialized inside the class definition)
	// Actual definition and initialization in cpp file
	static unordered_map<int, ProjectionParams_Albers> AlbersParamMap;
	static unordered_map<int, ProjectionParams_Lambert> LambertParamMap;
	static ProjectionParams_UTM initializeUTMParams(int MapProjection_WKID);

private:
	string Directory_Input_CLArg;
	int nDemDataCols;
	int nDemDataRows;
	double nDemDataxllcorner;
	double nDemDatayllcorner;
	int totalDischargeCount;

	//readRunoff_ObservedFile function reads Runoff_Observed_for_Calibration.csv file
	void readRunoff_ObservedFile();

	//RataDieNumber function converts Gregorian date YYYYMMDD into an integer number based on Rata Die System 
	int RataDieNumber(int y, int m, int d);

	//HydroPlusConfigFile.xml
	void readHydroPlusConfigFile();

	//Distributed cover maps: landcover.asc, imperviouscover.asc, treecover.asc, dem.asc, blockgroup.asc, AH_flux_Qtot_avg_Wpm2.asc
	//Note: File extension .asc is recognized by ArcGIS and QGIS as a standard input type when using the 6 header rows
	bool readNLCDLandCoverFile(string &Directory_Input_CLArg, double NODATA_code);
	void readTreeCoverFile(string &Directory_Input_CLArg, double NODATA_code);
	//bool readImperviousFile(string &Directory_Input_CLArg, double NODATA_code);
	void readImperviousFile(string& Directory_Input_CLArg, double NODATA_code);
	void readBlockGroupOutputFileMap(string& Directory_Input_CLArg, double NODATA_code);
	//readAnthropogenicHeatMap bool function to read anthropogenic heat maps from AH4GUC
	void readAnthropogenicHeatMap(string& Directory_Input_CLArg, double NODATA_code);
	void readHydroPlusConfigVariablesAsMap(string& Directory_Input_CLArg, double NODATA_code, string Variable_name, vector<double> &Variable_map);
	void LandCoverAssignment_by_NLCD_Class(int NumberOfLocations, double NODATA_code);

	void readRiverTimeSeriesFile(double startDate, double endDate);
	void readRiverGWTimeSeriesFile(double startDate, double endDate);
	void readRiverSWTimeSeriesFile(double startDate, double endDate);
	void readRiverStationSeriesFile(string &Directory_Input_CLArg, double NODATA_code);
	void readHecRasFile();

	//ECDynamic
	bool readDischargeData(string& Directory_Input_CLArg);
	void readPollutionData(string& Directory_Input_CLArg);
	void readExportCoefficientTable();
	void readSoilDepthMap_cm();
	void readHydraulicConductivityMap_umps();
	void calculateSoilTransmissivityChangeUnits();
	bool readWatertableMap_cm();
	//readWaterFeaturesMap function for SpatialBuffer inputs
	void readWaterFeaturesMap();
	void readNHDFlowDirectionMap();

	//writeAsciiMap takes a variableName and MapData_vec to create ArcGIS ready cover maps
	void writeAsciiMap(const string& variableName, const vector<double>& MapData_vec);
	//generate warning of printing out incomplete EC data
	void printIncompleteECDataWarning(const string& landCoverType, const string& landCoverRange, const string& huc);
	//LatitudeLongitude_Projection_Manager will create latitude longitude pair by calling Projection_Albers_to_LatitudeLongitude function
	void LatitudeLongitude_Projection_Manager(vector<vector<pair<double, double>>>& Latitude_Longitude_vectorPair_dd);
	void ValidateAndExtractDrawerFolderID(const string& label, const string& param, const vector<vector<unordered_map<string, double>>>& DataDrawers, const vector<vector<unordered_map<string, string>>>& StringDataDrawers, int& outDrawerID, int& outFolderID);
	//LoadScenarioParametersFromSimulationScenarios function checks HydroPlusConfig.xml for scenario data
	void LoadScenarioParametersFromSimulationScenarios();
};

#endif