#ifndef HeatFluxCal_H
#define HeatFluxCal_H
/*
* HydroPlus is source code and models developed and managed by Theodore Endreny and his students at SUNY ESF, te@esf.edu
* The US Forest Service and Davey Tree have provided strategic software development support through funding and personnel
* Attribution for use of software and code is requested. Please see Executive.cpp for documentation.
*/
#include <iostream>
#include <fstream>
#include <vector>
//_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 <stdio.h>
#include <algorithm> 
#include "../Inputs/Inputs.h"
#include "../DataFolder.h"
#include "ResistanceCal.h"
#include "../Inputs/WeatherProcessor.h"
#include "../Inputs/SolarCalculator.h"

//Note: Consider refactor HydroPlusConfig.xml to use a shared header file with constants, to ensure they are same across the .sln
//Note: Use these lines: #ifndef SHARED_CONSTANTS_H, #define SHARED_CONSTANTS_H, #define Variable_01 10, #endif

//Latent heat of vaporization for water (J/kg)
#define LatentHeatVaporization_Water_J_p_kg 2257000.0
//Stefan-Boltzmann constant (W/m^2/K^4) as 5.670374419E-8 W/(m^2*K^4)
#define Sigma_Stefan_Boltzmann 5.670374419E-8
//Specific heat capacity of air at constant pressure (J/kg/K) 1006.0 at ~290 to 315 K, 
#define SpecificHeat_DryAir_J_p_kg_K 1006.0
//Constant SpecificHeat_MoistAir_J_p_kg_p_C is specific heat of moist air at constant pressure, defined by Chin (2021) w Eq 13.37
#define SpecificHeat_MoistAir_J_p_kg_p_C 1013
//Ratio_MolecularWeightVapor_to_MolecularWeightDryAir is ratio of molecular weight for water vapor to dry air (g/g)
#define Ratio_MolecularWeightVapor_to_MolecularWeightDryAir 0.622
//Constant von Karman number; Range of values reported in research, recently tending towward 0.4
#define vonKarman_Constant 0.41
//Density of air (kg/m3)
#define Density_Air_kg_p_m3 1.205
//Conversion m to mm
#define Ratio_m_to_mm 1.0/1000.0
//Conversion km to m
#define Ratio_km_to_m 1.0/1000.0
//Ratio of Pascall to kiloPascall
#define Ratio_Pa_to_kPa 1000.0
//Ratio of Joules to megaJoules
#define Ratio_J_to_MJ 1000000.0
//Ratio of megaJoules to Joules
#define Ratio_MJ_to_J 1.0 / 1000000.0
//Conversion hour to day
#define Ratio_Hour_to_Day 24
//Constant environmental lapse rate dry (ELR) air (from Stull, 2000, estimated as Gamma_elr = 5.5 C/km), ...
//Note: ... there is no ELR equation but is often between 5 and 6 C/km, hence the estimate of 5.5 C.km
#define Gamma_EnvironmentalLapseRate_C_p_km 5.5
//Gas constant for pure water vapor is 0.461 or 4.61E-4 * 1000 to convert from kPa/(K*m3*g) to kPa/(K*m3*kg) 
#define Gas_constant_water_vapor_kPa_p_K_m3_kg 0.461
//Constant adiabatic lapse rate dry air (from Stull, 2000, defined after Eq 5.13, Gamma_d = 9.8 C/km)
//Note: Adiabatic lapse rate dry air is not used for estimating air 2 m above ground when summiting a mountain 
#define Gamma_dry_C_p_km 9.8
//Constant adiabatic lapse rate dew point air (from Stull, 2000, derived from Eq 5.8; 1/0.125 = 1/Z_lcl * (T - Td), if Z_lcl=1 km, T=9.8 then Td=1.8)
#define Gamma_dew_C_p_km 1.8
//Define constant threshold for fractional area
#define Constant_1E_negative4 1.E-4  
//Define constant threshold for temperature accuracy 0.001 or 0.0005
//Note: 0.0005 works better to ensure agreement between base and alternative case simulations when 2 digit precision
//Note: Consider refactor to include in HydroPlusConfig.xml in place of element IterationCount 
#define Temperature_ErrorThreshold_K 0.0005 
//Define coefficient c=0.08 Eq 28 of Yang et al (2013)
#define Coefficient_c_Eq28_Yang_etal_2013 0.08
//Define constant effective zero value
#define Constant_1E_negative6 1.E-6
//ZeroPlaneDisplacementHeight_frac (fraction) from Eq 13.3 Chin (2021) adjusts cover height to get zero plane displacement height
#define ZeroPlaneDisplacementHeight_frac 0.67
//Zom_RoughnessLengthMomentumTransferCoefficient roughness length coefficient controlling momentum transfer defined Eq 13.3 Chin (2021)
#define Zom_RoughnessLengthMomentumTransferCoefficient 0.123
//Zom_RoughnessLengthHeatVaporTransferCoefficient roughness length coefficient controlling heat and vapor transfer defined Eq 13.4 and following text Chin (2021)
#define Zom_RoughnessLengthHeatVaporTransferCoefficient 0.0123
//roughness length (m) for classification open, e.g. airport, from Davenport-Wieringa classification, from Table 4.1 in Stull (2000)
#define RoughnessLength_Airport_m 0.03
//RoughnessLength_Water_m = 0.00137 m, based on roughness length (m) for water set to z0 = 0.00137 m Chin (2021) and Thom and Oliver (1977)
//Note: Chin, D. A. (2021). Water Resources Engineering, Fourth Edition. Hoboken, NJ: Pearson Education.
//Note: Thom, A. S., & Oliver, H. R. (1977). On Penman's equation for estimating regional evaporation. Quarterly Journal of the Royal Meteorological Society, 103(436), 345-357. doi:https://doi.org/10.1002/qj.49710343610
#define RoughnessLength_Water_m 0.00137
//Roughness Lengths from Davenport-Wieringa Zo (m) to Classification and Landscape Table 4.1 in Stull (2000)
//Note: Classification smooth = Snow; Classification Open = Airport; Note UCAR Community Land Model uses 0.024 
#define RoughnessLength_Snow_m 0.005

//SGN function for finding sign of var using conditional ternary operator (?:)
#define SGN(var) ((var < 0) ? -1.0 : (var > 0) ? 1.0 : 0.0);
//MIN function for finding minimum of x_var and y_var using conditional ternary operator (?:)
//Note: Ternary operator if x_var>=y_var is true, then x_var is taken, otherwise y_var taken.
#define MIN(x_var,y_var) (((x_var)<=(y_var)) ? (x_var) : (y_var))        
//MAX function for finding maximum of x_var and y_var using conditional ternary operator (?:)
#define MAX(x_var,y_var) (((x_var)>=(y_var)) ? (x_var) : (y_var))        

//class HeatFluxCal
class HeatFluxCal
{
private:

	//SolarCalculator class pointer solarCalculator accesses function
	//Note: solarCalculator is owned by SimulationCoordinator so is not deleted with Destructor
	SolarCalculator* solarCalculator;
	Inputs* input;

	//ResistanceCal class pointer resistanceVar accesses values using resistanceVar->Resistance_Aerodynamic_Tree_s_p_m
	//Note: Resistances are: Resistance_Aerodynamic_Tree_s_p_m, Resistance_Aerodynamic_SVeg_s_p_m, Resistance_Aerodynamic_Impervious_s_p_m, Resistance_Aerodynamic_Soil_s_p_m, Resistance_Aerodynamic_CanopyToMesoLayer_s_p_m, Resistance_Surface_Soil_s_p_m, Resistance_Surface_Tree_s_p_m, Resistance_Surface_SVeg_s_p_m
	ResistanceCal* resistanceVar;

	int Diurnal_Hour_HH{};
	int Annual_Month_MM{};
	int LandCover_NLCD_Class{};
	double TreeCover_frac{};
	double ImperviousCover_frac{};
	double SoilCover_noTreeCanopy_frac{};
	double PermeablePavementCover_noTreeCanopy_frac{};
	double ShortVegCover_noTreeCanopy_frac{};
	double ImperviousCover_noTreeCanopy_frac{};
	double WaterCover_noTreeCanopy_frac{};
	double TreeCanopyCover_overImpervious_frac{};

	double Ratio_Road_to_Ground{};
	double Ratio_Garden_to_Ground{};
	double Direction_Road_rad{};
	int Symmetry_Wall_int{};
	double Glazing_Window{};

	//Elevation_lcl_km (km) is elevation of lifting condensation level
	double Elevation_lcl_km{};
	//Elevation_delta_km (km) is elevation difference, Elevation_pixel - Elevation_referenceStation
	double Elevation_delta_km{};

	//Flag_UrbanCanyonResistance = 1 to simulate urban canyon, = 0 to ignore urban canyon 
	int Flag_UrbanCanyonResistance{};
	double Count_CanyonWalls{};


	//Tair_mesoScale_Prior_K (K) mesoscale air temperature from previous iteration in numerical method search 
	double Tair_mesoScale_Prior_K{};
	//AbsHumidity_mesoScale_Prior_kg_p_m3 (K) mesoscale absolute humidity from previous iteration in numerical method search 
	double AbsHumidity_mesoScale_Prior_kg_p_m3{};

	//Tair_localCell_K (K) canopy air temperature
	double Tair_localCell_K{};
	//Tair_referenceCell_K (K) used in error handling
	double Tair_referenceCell_K{};
	//Tdew_localCell_K (K) canopy dew point temperature
	double Tdew_localCell_K{};
	//WindSpeed_Reference_m_p_s (m/s) wind speed
	double WindSpeed_Reference_m_p_s{};

	// saturated vapor pressure (absolute humidity Kg/m^3) at Ta
	double AbsHumiditySat_localCell_kg_p_m3{};
	// real vapor pressure (absolute humidity Kg/m^3 at Tdew_localCell_K)
	double AbsHumidityAct_localCell_kg_p_m3{};
	double AbsHumidityGradient{};
	double Radiation_Longwave_Down_W_p_m2{};
	double Radiation_Longwave_Up_W_p_m2{};
	double Cos_IncidenceAngle_Solar{};
	double Radiation_Shortwave_W_p_m2{};
	double Radiation_Shortwave_Direct_W_p_m2{};
	double Radiation_Shortwave_Diffuse_W_p_m2{};

	double Radiation_Shortwave_Direct_Diffuse_TotalCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_Pavement_UrbanCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_Garden_UrbanCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_Tree_UrbanCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_SVeg_UrbanCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_Soil_UrbanCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_Water_UrbanCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_Wall_UrbanCanyon_Absorbed_W_p_m2{};
	double Radiation_Shortwave_Direct_Diffuse_Roof_UrbanCanyon_Absorbed_W_p_m2{};

	double Radiation_Longwave_Down_to_Wall_W_p_m2{};
	double Radiation_Longwave_Down_to_Ground_W_p_m2{};
	double Radiation_Longwave_Net_Pavement_W_p_m2{};
	double Radiation_Longwave_Net_Tree_W_p_m2{};
	double Radiation_Longwave_Net_SVeg_W_p_m2{};
	double Radiation_Longwave_Net_Soil_W_p_m2{};
	double Radiation_Longwave_Net_Water_W_p_m2{};
	double Radiation_Longwave_Net_Wall_W_p_m2{};
	double Temperature_UrbanCanyonRoad_K{};
	double Temperature_UrbanCanyonTree_K{};
	double Temperature_UrbanCanyonSVeg_K{};
	double Temperature_UrbanCanyonSoil_K{};
	double Temperature_UrbanCanyonWater_K{};
	double Temperature_UrbanCanyonWall_K{};

	// available water on Impervious for evaporate
	double Saturation_ImperviousStorage_frac{};
	double Saturation_PerviousStorage_frac{};
	double Saturation_TreeStorage_frac{};
	double Saturation_SVegStorage_frac{};
	//Saturation_SoilStorage_frac (frac) ranges from 0=dry, 1=saturated, available water in the soil to evaporate
	double Saturation_SoilStorage_frac{};
	//Saturation_WaterStorage_frac (frac) is always presumed 1=saturated, term used to maintain equation form
	double Saturation_WaterStorage_frac{};

	double LAI_BAI_Tree_m2pm2{};
	double LAI_BAI_SVeg_m2pm2{};
	double LAI_Tree_m2pm2{};
	double LAI_SVeg_m2pm2{};

	//Height_LandCover_Average_m is average or blended height of all cover types
	double Height_LandCover_Average_m{};
	//GroundHeatFlux_NightTimeLag_hr default value is 1
	int GroundHeatFlux_NightTimeLag_hr = 1;

	double Albedo_Pavement_frac{};
	double Emissivity_Pavement_frac{};
	double Albedo_Roof_frac{};
	double Emissivity_Roof_frac{};
	double Albedo_Road_frac{};
	double Emissivity_Road_frac{};
	double Albedo_Garden_frac{};
	double Emissivity_Garden_frac{};

	double Albedo_Wall_frac{};
	double Emissivity_Wall_frac{};
	double Albedo_Window_frac{};
	double Emissivity_Window_frac{};

	double VaporPressure_Actual_kPa{};
	double Emissivity_Sky_frac{};

	// impervious intermediate heat flux values
	double ImpNR_W_p_m2{};
	double ImpDeltaQ_W_p_m2{};
	double ImpC{};
	double ImpH_W_p_m2{};
	double ImpLE_W_p_m2{};

	double RoofNR_W_p_m2{};
	double RoofDeltaQ_W_p_m2{};
	double RoofC{};
	double RoofH_W_p_m2{};
	double RoofLE_W_p_m2{};

	double SoilNR_W_p_m2{};
	double SoilDeltaQ_W_p_m2{};
	double SoilC{};
	double SoilH_W_p_m2{};
	double SoilLE_W_p_m2{};

	double TreeNR_W_p_m2{};
	double TreeDeltaQ_W_p_m2{};
	double TreeCE{};
	double TreeCT{};
	double TreeC{};
	double TreeH_W_p_m2{};
	double TreeLEE_W_p_m2{};
	double TreeLET_W_p_m2{};
	double TreeLE_W_p_m2{};

	double SVegNR_W_p_m2{};
	double SVegDeltaQ_W_p_m2{};
	double SVegCE{};
	double SVegCT{};
	double SVegC{};
	double SVegH_W_p_m2{};
	double SVegLEE_W_p_m2{};
	double SVegLET_W_p_m2{};
	double SVegLE_W_p_m2{};

	double WaterNR_W_p_m2{};
	double WaterDeltaQ_W_p_m2{};
	double WaterC{};
	double WaterH_W_p_m2{};
	double WaterLE_W_p_m2{};
	double LE_total_W_p_m2{};
	double NR_total_W_p_m2{};
	double DeltaQ_W_p_m2{};
	double AH_total_W_p_m2{};

	double AnthropogenicHeat_Flux_Total_Wpm2{};
	double Imp_AH_W_p_m2{};
	double Soil_AH_W_p_m2{};
	double Tree_AH_W_p_m2{};
	double SVeg_AH_W_p_m2{};
	double Water_AH_W_p_m2{};

	//NR_priorTS_W_p_m2_Map map used with StationInfoMap for DeltaQ calculations across Multiple Weather Stations 
	map<string, double> NR_priorTS_W_p_m2_Map;
	double ImpNR_priorTS_W_p_m2{};
	double TreeNR_priorTS_W_p_m2{};
	double SVegNR_priorTS_W_p_m2{};
	double SoilNR_priorTS_W_p_m2{};
	double WaterNR_priorTS_W_p_m2{};
	double RoofNR_priorTS_W_p_m2{};

	// soil parameters
	// average soil height
	double Height_Soil_m{};
	// soil effective albedo
	double Albedo_Soil_frac{};
	// soil emissivity
	double Emissivity_Soil_frac{};

	double a1_OHM_Soil_frac{};
	double a2_OHM_Soil_hr{};
	double a3_OHM_Soil_W_p_m2{};

	// urban parameters
	double Albedo_Urban_frac{};
	double Emissivity_Urban_frac{};
	double Height_UrbanBuilding_m{};
	double Width_UrbanCanyon_m{};
	double a1_OHM_Pavement_frac{};
	double a2_OHM_Pavement_hr{};
	double a3_OHM_Pavement_W_p_m2{};
	double a1_OHM_Building_frac{};
	double a2_OHM_Building_hr{};
	double a3_OHM_Building_W_p_m2{};

	double a1_OHM_Urban_frac{};
	double a2_OHM_Urban_hr{};
	double a3_OHM_Urban_W_p_m2{};
	double Building_Area_frac{};
	double Storage_Depression_Max_Impervious_m{};
	double Storage_Depression_Max_Pervious_m{};

	// parameters calculated based on the parameters above
	double AspectRatio_Canyon_frac{};
	double SkyViewFactor_Road_frac{};
	double SkyViewFactor_Wall_frac{};
	double SkyViewFactor_Garden_frac{};
	double SkyViewFactor_Ground_Open_frac{};
	double SkyViewFactor_Ground_Wall_frac{};

	//tree parameters

	double a1_OHM_Tree_frac{};
	double a2_OHM_Tree_hr{};
	double a3_OHM_Tree_W_p_m2{};

	double Height_Tree_m{};
	double Albedo_Tree_frac{};
	double Emissivity_Tree_frac{};
	double Storage_per_LAI_BAI_Tree_m{};

	// short vegetation parameters
	double a1_OHM_SVeg_frac{};
	double a2_OHM_SVeg_hr{};
	double a3_OHM_SVeg_W_p_m2{};

	double Height_SVeg_m{};
	double Albedo_SVeg_frac{};
	double Emissivity_SVeg_frac{};
	double Storage_per_LAI_BAI_SVeg_m{};

	// water parameters

	double Height_Water_m{};
	double Albedo_Water_frac{};
	double Emissivity_Water_frac{};
	double a1_OHM_Water_frac{};
	double a2_OHM_Water_hr{};
	double a3_OHM_Water_W_p_m2{};

	double ZenithAngle_Solar_rad{};
	double Cos_ZenithAngle_Solar{};
	double AzimuthAngle_Solar_N_0_rad{};
	//Ratio_AlbedoDirCZA_to_AlbedoDirCZA60  from Eq 15 of Yang et al. (2008)
	double Ratio_AlbedoDirCZA_to_AlbedoDirCZA60{};
	double AtmosphericBoundaryLayer_Thickness_km{};
	double ABL_adjust_frac{};
	double Resistance_Aerodynamic_CanopyToMesoLayer_max_s_p_m{};
	double Tair_mesoScale_final_pointer_K{};
	double H_total_pointer_W_p_m2{};

public:
	
	//HeatFluxCal(Inputs* input, SolarCalculator* solarCalc) constructor must match .h class name
	HeatFluxCal(Inputs* input, SolarCalculator* solarCalc);
	//Note: Variables defined in the constructor are shared across all void functions, otherwise limited to void function
	//HeatFluxCal(Inputs* input, SolarCalculator* solarCalc) : solarCalculator(solarCalc) {}
	//~HeatFluxCal() destructor (without body {}) that ensures resources used by the object can be released 
	~HeatFluxCal();

	//Load_TemperatureCalculationParams function initializes HydroPlusConfig TemperatureCalculationParams for folders with and without reference station 
	void Load_TemperatureCalculationParams(Inputs* input, int MapPixel_ID, int DataFolder_ID, int timeStep);

	//EnergyFlux_WaterFlux_DefineTerms function calculates energy balance once for each folder at each timestep, before heading into the water balance calculations
	void EnergyFlux_WaterFlux_DefineTerms(Inputs* input, CompactRagged* beC, bool Flag_simulateReferenceStation, int timeStep, WeatherProcessor* WeatherPro, string StationID_string, int MapPixel_ID, int DataFolder_ID);

	//UrbanCanyon_RadiationLongwave function calculates urban canyon longwave energy balance once for each folder at each timestep
	void UrbanCanyon_RadiationLongwave(Inputs* input, CompactRagged* beC, bool Flag_UrbanCanyonResistance, int MapPixel_ID, int DataFolder_ID, int timeStep);
	//UrbanCanyon_RadiationShortwave function calculates urban canyon direct and diffuse shortwave energy balance once for each folder at each timestep
	void UrbanCanyon_RadiationShortwave(Inputs* input, bool Flag_UrbanCanyonResistance, int MapPixel_ID, int DataFolder_ID, int timeStep);

	//Tair_AbsHumidity_LayerBalance function balances air temperature and absolute humidity between the urban canopy and mesoscale layers
	void Tair_AbsHumidity_LayerBalance(Inputs* input, CompactRagged* beC, DataOrganizer* organizer, WeatherProcessor* WeatherPro, int MapPixel_ID, int DataFolder_ID, int timeStep);
	//Tair_AbsHumidity_LayerBalance_NumericalMethod function only called for non-reference cells 
	void Tair_AbsHumidity_LayerBalance_NumericalMethod(Inputs* input, CompactRagged* beC, DataOrganizer* organizer, WeatherProcessor* WeatherPro,
		int MapPixel_ID, int DataFolder_ID, int timeStep, int iterationCount);

	//AdiabaticLapseRates function will adjust air and dew point temperatures based on elevation, modifying assumption of homogeneous mesoclimate
	void AdiabaticLapseRates(Inputs* input, int MapPixel_ID, int DataFolder_ID, int timestep, WeatherProcessor* WeatherPro);
	//RescaleVariables_CoolAir function will rescale latent and sensible energy fluxes, compute corresponding mesoScale values for Tair and AbsoluteHumidity
	void RescaleVariables_CoolAir(Inputs* input, CompactRagged* beC, DataOrganizer* organizer, int MapPixel_ID, int DataFolder_ID, int timestep);
	//Rescale_LE_H_CoolAir
	void Rescale_LE_H_CoolAir(Inputs* input, int MapPixel_ID, int DataFolder_ID, int timeStep);
	//CollectVariables_CoolAir function will store all values to memory
	void CollectVariables_CoolAir(Inputs* input, CompactRagged* beC, SolarCalculator* solarCalc, int MapPixel_ID, int DataFolder_ID, int timestep);
	//used to initialize the AtmosphericBoundaryLayer parameters
	void Initialize_AtmosphericBoundaryLayer_param();
	//Tair_mesoScale_K (K) mesoscale air temperature at current calculation
	double Tair_mesoScale_K{};
	double Tair_mesoScale_final_K{};
	double Tdew_mesoScale_K{};
	double H_total_W_p_m2{};
	//AbsHumidity_mesoScale_kg_p_m3 (kg/m3) mescoscale air absolute humidity 
	double AbsHumidity_mesoScale_kg_p_m3{};
	double AbsHumidity_mesoScale_final_kg_p_m3{};
	
	//HeatFlux_InitializeVectors function to initialize vectors
	static void HeatFlux_InitializeVectors(Inputs* input); 
	//NetRadiation_Prior_Calc function to update prior values of net radation for OHM; NetRadiation_W_p_m2_vec stores pair of hours-ago and NR
	double NetRadiation_Prior_Calc(Inputs* input, int timeStep, double NR_vec_hrs, vector<pair<double, double>>& NetRadiation_W_p_m2_vec, double NetRadiation_W_p_m2);
	//to make sure all vectors only been initialized once while having multiple runs of CollectVariables_CoolAir
	static bool vectors_initialized;
	//Note: static vectors need to also be defined at the start of the .cpp file scoped to HeatFluxCal::
	static vector<double> AH_total_W_p_m2_vec;
	static vector<double> Imp_AH_W_p_m2_vec;
	static vector<double> Tree_AH_W_p_m2_vec;
	static vector<double> SVeg_AH_W_p_m2_vec;
	static vector<double> Soil_AH_W_p_m2_vec;
	static vector<double> Water_AH_W_p_m2_vec;
	static vector<double> NR_total_W_p_m2_vec;
	static vector<double> DeltaQ_W_p_m2_vec;
	static vector<double> WaterDeltaQ_W_p_m2_vec;
	static vector<double> SoilDeltaQ_W_p_m2_vec;
	static vector<double> SVegDeltaQ_W_p_m2_vec;
	static vector<double> TreeDeltaQ_W_p_m2_vec;
	static vector<double> ImpDeltaQ_W_p_m2_vec;
	//Note: static vectors need to also be defined at the start of the .cpp file scoped to HeatFluxCal::
	static vector<pair<double, double>> ImpNR_priorTS_W_p_m2_vec;
	static vector<pair<double, double>> TreeNR_priorTS_W_p_m2_vec;
	static vector<pair<double, double>> SVegNR_priorTS_W_p_m2_vec;
	static vector<pair<double, double>> SoilNR_priorTS_W_p_m2_vec;
	static vector<pair<double, double>> WaterNR_priorTS_W_p_m2_vec;
	static vector<pair<double, double>> RoofNR_priorTS_W_p_m2_vec;
};

#endif