#include "HeatMetrics_Calc.h"

//HeatMetrics_Calc::HeatMetrics_Calc Constructor for use here with initialization list or pointer
//Note: Initialization list : or colon after function name initializes SolarCalculator and resistanceVar objects
HeatMetrics_Calc::HeatMetrics_Calc(Inputs* input)
{
	//Cp_Ratio_DryAir_to_WaterVapor =  SpecificHeat_DryAir_J_p_kg_K * Molecular_Weight_DryAir_kg_p_kmol / Molecular_Weight_WaterVapor_kg_p_kmol 
	//Note: Ratio relates the thermal properties of air to water vapor
	Cp_Ratio_DryAir_to_WaterVapor = SpecificHeat_DryAir_J_p_kg_K * Molecular_Weight_DryAir_kg_p_kmol / Molecular_Weight_WaterVapor_kg_p_kmol;
	//Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K (m2/s2K) is R_Universal_Gas_Constant_J_p_kmol_K / Molecular_Weight_DryAir_kg_p_kmol; 
	//Note: Units  m2/(s2 K) = J/(kg K) = [J/(kmol K)] / [kg/kmol]
	Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K = R_Universal_Gas_Constant_J_p_kmol_K / Molecular_Weight_DryAir_kg_p_kmol;
	//Prandtl_Number_wind is modified Prandtl number relating momentum diffusivity or viscosity to thermal diffusivity
	Prandtl_Number_wind = SpecificHeat_DryAir_J_p_kg_K / (SpecificHeat_DryAir_J_p_kg_K + 1.25 * Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K); 

	Wick_Emissivity_frac = 0.95;
	Wick_Albedo_frac = 0.4;
	Wick_Diameter_m = 0.007;
	Wick_Length_m = 0.0254;
	Globe_Emissivity_frac = 0.95;
	Globe_Albedo_frac = 0.05;
	Globe_Diameter_m = 0.0508;
	SurfaceGround_Emissivity_frac = 0.999;
	SurfaceGround_Albedo_frac = 0.45;
}

//HeatIndex_K_Calc function will manage calculation of the heat index temperature and return value
double HeatMetrics_Calc::HeatIndex_K_Calc(Inputs* input, double Tair_K, double Tdew_K) {
	//Note: Heat Index from Rothfusz NWS (1990) algorithm at: https://www.wpc.ncep.noaa.gov/html/heatindex_equation.shtml and Steadman (1979)
	//Note: Heat Index algorithm below does not necessarily handle extreme temperature and relative humidity
	//References: 
	//Rothfusz NWS Technical Attachment SR 90-23 (1990) and adapted to fit Steadman (1979) results for all temperature and humidity levels
	//Steadman, R. G. (1979). "The Assessment of Sultriness. Part I: A Temperature-Humidity Index Based on Human Physiology and Clothing Science." Journal of Applied Meteorology 18(7): 861-873.

	//Tair_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tair_C = Tair_K - 273.15;
	//Tdew_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tdew_C = Tdew_K - 273.15;
	//Tair_F (F) for NOAA NWS HeatIndex Formula is converted from Tair_C using standard equation, Tair_C * 9/5 + 32
	double Tair_F = Tair_C * 9.0 / 5.0 + 32;
	//VaporPressure_Saturated_kPa (kPa) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
	double VaporPressure_Saturated_kPa = 0.6108 * exp((17.27 * Tair_C) / (237.3 + Tair_C));
	//VaporPressure_Actual_kPa (kPa) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
	double VaporPressure_Actual_kPa = 0.6108 * exp((17.27 * Tdew_C) / (237.3 + Tdew_C));
	//RelativeHumidity_percent (%) is ratio of VaporPressure_Actual_kPa / VaporPressure_Saturated_kPa , Eq 5.6 from Stull (2000)
	double RelativeHumidity_percent = VaporPressure_Actual_kPa / VaporPressure_Saturated_kPa * 100.0;
	//HeatIndex_NWS_F (F) and HeatIndex_NWS_K (K) defined as empty
	double HeatIndex_NWS_F = 0;
	double HeatIndex_NWS_K = 0;

	//HeatIndex_NWS_F (F) is 1st defined using Eq 2 of Rothfusz NWS (1990) algorithm
	//Note: HeatIndex_NWS_F = 0.5 * (Tair_F + 61.0 + ((Tair_F - 68.0) * 1.2) + (RelativeHumidity_percent * 0.094))
	HeatIndex_NWS_F = 0.5 * (Tair_F + 61.0 + ((Tair_F - 68.0) * 1.2) + (RelativeHumidity_percent * 0.094));
	//HeatIndex_NWS_F (F) is 2nd defined using text after Eq 2 of Rothfusz NWS (1990) algorithm, averaging HeatIndex_NWS_F and Tair_F
	//Note: HeatIndex_NWS_F = 0.5 * (HeatIndex_NWS_F + Tair_F)
	HeatIndex_NWS_F = 0.5 * (HeatIndex_NWS_F + Tair_F);

	//If HeatIndex_NWS_F >= 80.0 F then utilize Rothfusz (1990) regression, but consider Else If terms from Rothfusz NWS (1990) algorithm
	if (HeatIndex_NWS_F >= 80.0) {
		//HeatIndex_NWS_F (F) is 3rd defined using Eq 1 of the Rothfusz NWS (1990) algorithm
		HeatIndex_NWS_F = -42.379 + (2.04901523 * Tair_F) + (10.14333127 * RelativeHumidity_percent) - (0.22475541 * Tair_F * RelativeHumidity_percent) -
			(0.00683783 * Tair_F * Tair_F) - (0.05481717 * RelativeHumidity_percent * RelativeHumidity_percent) + (0.00122874 * Tair_F * Tair_F * RelativeHumidity_percent) + (0.00085282 * Tair_F * RelativeHumidity_percent * RelativeHumidity_percent) - (0.00000199 * Tair_F * Tair_F * RelativeHumidity_percent * RelativeHumidity_percent);
	}
	//Else If RelativeHumidity_percent < 13% and 80 F < Tair_F < 112 F then utilize adjustment from Rothfusz NWS (1990) algorithm
	else if (RelativeHumidity_percent < 13.0 && Tair_F > 80.0 && Tair_F < 112.0) {
		//HeatIndex_Adjust_F (F) is 1st adjustment below Eq 1 of Rothfusz NWS (1990) algorithm
		double HeatIndex_Adjust_F = ((13.0 - RelativeHumidity_percent) / 4.0) * pow(((17.0 - abs(Tair_F - 95.0)) / 17.0), 0.5);
		//HeatIndex_NWS_F (F) is decreased by HeatIndex_Adjust_F (F)
		HeatIndex_NWS_F = HeatIndex_NWS_F - HeatIndex_Adjust_F;
	}
	//Else If RelativeHumidity_percent > 85% and 80 F < Tair_F < 87 F then utilize adjustment from Rothfusz NWS (1990) algorithm
	else if (RelativeHumidity_percent > 85.0 && Tair_F > 80.0 && Tair_F < 87.0) {
		//HeatIndex_Adjust_F (F) is 2nd adjustment below Eq 1 of Rothfusz NWS (1990) algorithm
		double HeatIndex_Adjust_F = ((RelativeHumidity_percent - 85.0) / 10.0) * ((87.0 - Tair_F) / 5.0);
		//HeatIndex_NWS_F (F) is increased by HeatIndex_Adjust_F (F)
		HeatIndex_NWS_F = HeatIndex_NWS_F + HeatIndex_Adjust_F;
	}
	//Else HeatIndex_NWS_F < 80 F then use Tair_K
	else {
		//HeatIndex_NWS_F (F) defined as Tair_F (F)
		HeatIndex_NWS_F = Tair_F;
	}

	//HeatIndex_NWS_K (K) is derived from HeatIndex_NWS_F using conversion to degrees C, then to degrees K
	HeatIndex_NWS_K = ((HeatIndex_NWS_F - 32) * 5.0 / 9.0) + 273.15;

	//return HeatIndex_NWS_K
	return HeatIndex_NWS_K;
}

//Humidex_K_Calc_Calc function will manage calculation of the humidex temperature and return value
double HeatMetrics_Calc::Humidex_K_Calc(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID) {
	//Note: Humidex from Diaconescu et al. (2023) Eq 1 and 2, where Eq 2 provided vapor pressure, e, and was replaced with Eq 4.2.2 of Shuttleworth (1993)
	//Note: Humidex values were within 2% of each other when computed using e from Eq 2 of Diaconescu et al. (2023) or Eq 4.2.2 of Shuttleworth (1993)
	//References: 
	//Diaconescu, E., H. Sankare, K. Chow, T. Q. Murdock and A. J. Cannon (2023). "A short note on the use of daily climate data to calculate Humidex heat-stress indices." International Journal of Climatology 43(2): 837-849.

	//Tair_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tair_C = beC->by_key(MapPixel_ID, DataFolder_ID, "Tair_K") - 273.15;
	//Tdew_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tdew_C = beC->by_key(MapPixel_ID, DataFolder_ID, "Tdew_K") - 273.15;
	//VaporPressure_Actual_kPa (kPa) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
	double VaporPressure_Actual_kPa = 0.6108 * exp((17.27 * Tdew_C) / (237.3 + Tdew_C));
	
	//VaporPressure_Actual_hPa (hPa) is VaporPressure_Actual_kPa * 10
	double VaporPressure_Actual_hPa = VaporPressure_Actual_kPa * 10;
	//Humidex_C (C) is Tair_C + (5.0/9.0) * (VaporPressure_Actual_hPa - 10); 
	//Note: 5/9 is sometimes given as 0.555, and is scaling VaporPressure_Actual_hPa
	//Note: From Diaconescu et al. (2023) Eq 1 and 2, where Eq 2 provided vapor pressure, e, and was replaced with Eq 4.2.2 of Shuttleworth (1993)
	double Humidex_C = Tair_C + (5.0 / 9.0) * (VaporPressure_Actual_hPa - 10);
	//Humidex_K (K) is Humidex_C + 273.15
	double Humidex_K = Humidex_C + 273.15;

	//return Humidex_K
	return Humidex_K;
}

//WetBulbGlobeTemperature_C_Calc function will manage calculation of the wet bulb globe temperature (WBGT) and return value
//Note: WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Note: Radiation_Shortwave_Direct_W_p_m2 specific to MapPixel_ID when Flag_LatitudeLongitude_Computed is true
//References:
//EPA (U.S.Environmental Protection Agency). (2000). Meteorological data processing. In Meteorological Monitoring Guidance for Regulatory Modeling Applications (EPA-454 / R-99-005). Washington, D.C.
//Liljegren, J. C., R. A. Carhart, P. Lawday, S. Tschopp, and R. Sharp. (2008). Modeling the Wet Bulb Globe Temperature Using Standard Meteorological Measurements. The Journal of Occupational and Environmental Hygiene, vol. 5:10, pp. 645-655. DOI: 10.1080/15459620802310770 
//Murray, F.W. 1967. On the computation of saturation vapour pressure. J. Applied Meteorology 6: 203-204. https://journals.ametsoc.org/view/journals/apme/6/1/1520-0450_1967_006_0203_otcosv_2_0_co_2.xml
//Shuttleworth, J. W. (1993). Evaporation. In D. R. Maidment (Ed.), Handbook of Hydrology (pp. 4.1-4.5.3). New York: McGraw-Hill.
void HeatMetrics_Calc::WetBulbGlobeTemperature_C_Calc(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID, int timeStep, double Tair_K, double Tdew_K, double Radiation_Shortwave_Direct_W_p_m2, double Radiation_Shortwave_Diffuse_W_p_m2, double Cos_ZenithAngle_Solar, double& Twetbulbglobe_C, double& Twetbulbnatural_C, double& Twetbulbpsychrometric_C, double& UTCI_C) {

	//YYYYMMDD is input->SimulationDate_YYYYMMDD[timeStep]
	int Date_YYYYMMDD = input->SimulationDate_YYYYMMDD[timeStep];
	//HH_MM_SS is input->SimulationTime_HMS[timeStep]
	string HH_MM_SS = input->SimulationTime_HMS[timeStep];

	//Year_YYYY extracted from YYYYMMDD, first 4 spaces
	Year_YYYY = stoi(to_string(Date_YYYYMMDD).substr(0, 4));
	//Month_0_Julian_or_1_to_12_MM extracted from YYYYMMDD, next 2 spaces
	Month_MM = stoi(to_string(Date_YYYYMMDD).substr(4, 2));
	//Day_of_Month_or_Julian_if_Month_0_DD extracted from YYYYMMDD, next 2 spaces
	Day_DD = stoi(to_string(Date_YYYYMMDD).substr(6, 2));
	//Hour_HH extracted from HH_MM_SS; values must be 1 to 24
	Hour_HH = stoi(HH_MM_SS.substr(0, 2)) + 1;
	//Minute_Mn extracted from HH_MM_SS
	Minute_Mn = stoi(HH_MM_SS.substr(3, 2));
	//Tair_C = Tair_K - 273.15
	Tair_C = Tair_K - 273.15;
	//Tdew_C = Tdew_K - 273.15
	Tdew_C = Tdew_K - 273.15;

	//Time_averaging_MeteorologicalInputs_minutes (min) is set to thirty (30) min by Liljegren et al. (2008) code
	Time_averaging_MeteorologicalInputs_minutes = 30;
	//AtmosphericPressure_hPa = input->AtmPres_kPa[timeStep] * 10.0, convert kPa to hPa
	AtmosphericPressure_hPa = input->AtmPres_kPa[timeStep] * 10.0;

	//VaporPressure_Saturated_kPa (kPa) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
	double VaporPressure_Saturated_kPa = 0.6108 * exp((17.27 * Tair_C) / (237.3 + Tair_C));
	//VaporPressure_Actual_kPa (kPa) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
	double VaporPressure_Actual_kPa = 0.6108 * exp((17.27 * Tdew_C) / (237.3 + Tdew_C));
	//RelativeHumidity_percent computed as VaporPressure_Actual_kPa / VaporPressure_Saturated_kPa * 100
	RelativeHumidity_percent = VaporPressure_Actual_kPa / VaporPressure_Saturated_kPa * 100;

	//WindSpeed_mps
	WindSpeed_mps = input->WindSpd_mps[timeStep];
	//Height_WindSpeed_Actual_m (m) is SimulationNumericalParams["Height_Sensor_Wind_m"], reference height for wind speed measurements
	Height_WindSpeed_Actual_m = input->SimulationNumericalParams["Height_Sensor_Wind_m"];
	//Delta_Height_Tair_m (m) is set to zero (0) by Liljegren et al. (2008) code
	Delta_Height_Tair_m = 0;
	//Flag_urban is set to one (1) by Liljegren et al. (2008) code
	Flag_urban = 1;
	//Flag_urban is initialized to one (1) by Liljegren et al. (2008) code
	Flag_daytime = 1;

	//Ratio_LAI_BAI_Tree_Actual_to_Max_frac (frac) is  ratio of LAI_BAI_Tree_m2_p_m2[timeStep] to LAI_BAI_map["LAI_BAI_Tree_max_m2pm2"]
	double Ratio_LAI_BAI_Tree_Actual_to_Max_frac = input->LAI_BAI_Tree_m2_p_m2[timeStep] / input->LAI_BAI_map["LAI_BAI_Tree_max_m2pm2"];

	//Radiation_Shortwave_Direct_W_p_m2 (W/m2) is Radiation_Shortwave_Direct_W_p_m2 * (1 - (TreeCover_frac * Ratio_LAI_BAI_Tree_Actual_to_Max_frac ))
	//Note: Radiation_Shortwave_Direct_W_p_m2 is reduced based on tree cover shade, with a maximum when Ratio_LAI_BAI_Tree_Actual_to_Max_frac = 1
	//Note: Consider Refactor to create shade based on NLCD Class 21-24 HeatFluxCall SkyViewFactor_Ground_Open_frac or ...
	//Note: ... Area_Roof_Urban21_frac, or ratio Height_Avg_Urban22_m to Width_Road_Urban22_m; but Radiation_Shortwave_Direct_W_p_m2 may include this influence
	Radiation_Shortwave_Direct_W_p_m2 = Radiation_Shortwave_Direct_W_p_m2 * (1 - (beC->by_key(MapPixel_ID, DataFolder_ID, "TreeCover_frac") * Ratio_LAI_BAI_Tree_Actual_to_Max_frac));
	//Radiation_Shortwave_Wpm2 is Radiation_Shortwave_Direct_W_p_m2 + Radiation_Shortwave_Diffuse_W_p_m2, combining direct and diffuse
	Radiation_Shortwave_Wpm2 = Radiation_Shortwave_Direct_W_p_m2 + Radiation_Shortwave_Diffuse_W_p_m2;

	//If Radiation_Shortwave_Wpm2 > 0 then 
	if (Radiation_Shortwave_Wpm2 > 0) {
		//Radiation_Shortwave_direct_frac is Radiation_Shortwave_Direct_W_p_m2 / Radiation_Shortwave_Wpm2
		//Note: Liljegren et al. (2008) code called calc_solar_parameters and solarposition functions to obtain radiation parameters
		//Note : Liljegren et al. (2008) explain Eq 13 is empirically derived from 1 year of data from Oklahoma, shown in Figure 
		Radiation_Shortwave_direct_frac = Radiation_Shortwave_Direct_W_p_m2 / Radiation_Shortwave_Wpm2;
	}
	//Else Radiation_Shortwave_direct_frac is 0
	else {
		Radiation_Shortwave_direct_frac = 0;
	}

	//If Height_WindSpeed_Actual_m from HydroPlusConfig.xml not equal to Height_WindSpeed_Reference_m for Twbg calculations, then enter
	//Note: Height_WindSpeed_Actual_m = 10 m for most NOAA weather stations, while Height_WindSpeed_Reference_m = 2 m, so lower wind speed
	if (Height_WindSpeed_Actual_m != Height_WindSpeed_Reference_m) {
		//Flag_daytime set to True if Cos_ZenithAngle_Solar > 0, or False otherwise
		Flag_daytime = (Cos_ZenithAngle_Solar > 0) ? 1 : 0;
		//StabilityClass_Atmosphere_int is determined from function StabilityClass_Estimator
		StabilityClass_Atmosphere_int = HeatMetrics_Calc::StabilityClass_Estimator(Flag_daytime, WindSpeed_mps, Radiation_Shortwave_Wpm2, Delta_Height_Tair_m);
		//WindSpeed_estimated_mps (m/s) computed in WindSpeed_Estimate_mps_Calc to adjust from measured height to Twbg height
		//Note: Height_WindSpeed_Actual_m = 10 m for most NOAA weather stations, while Height_WindSpeed_Reference_m = 2 m, so lower wind speed
		WindSpeed_estimated_mps = HeatMetrics_Calc::WindSpeed_Estimate_mps_Calc(WindSpeed_mps, Height_WindSpeed_Actual_m, Height_WindSpeed_Reference_m, StabilityClass_Atmosphere_int, Flag_urban);
		//WindSpeed_mps = WindSpeed_estimated_mps, updated based on stability
		WindSpeed_mps = WindSpeed_estimated_mps;
	}

	//WindSpeed_mps limited to minimum wind speed of 0.13 m/s used by Liljegren et al. (2008) code
	//WindSpeed_min_mps = 0.13
	double WindSpeed_min_mps = 0.13;
	//WindSpeed_mps = max(WindSpeed_mps, WindSpeed_min_mps)
	WindSpeed_mps = max(WindSpeed_mps, WindSpeed_min_mps);

	//RelativeHumidity_frac(frac) converted from RelativeHumidity_percent * 0.01
	double RelativeHumidity_frac = RelativeHumidity_percent * 0.01;

	//Calculate globe and natural wet bulb temperatures
	//Tglobe_C computed in Tglobe_C_Calc
	Tglobe_C = HeatMetrics_Calc::Tglobe_C_Calc(input, Tair_K, RelativeHumidity_frac, AtmosphericPressure_hPa, WindSpeed_mps, Radiation_Shortwave_Wpm2, Radiation_Shortwave_direct_frac, Cos_ZenithAngle_Solar, MapPixel_ID, timeStep);

	//Tglobe_K = Tglobe_C + 273.15; convert from C to K
	double Tglobe_K = Tglobe_C + 273.15;
	//Tglobe_K (K) is scaled to 2 digits precision
	//Note: RHS function divides by 100 after taking integer of the term that takes product with 100 to round up
	beC->by_key(MapPixel_ID, DataFolder_ID, "Tglobe_K") = Inputs::roundToDecimalPlaces(Tglobe_K, 2);

	//Twetbulbnatural_C computed in Twetbulb_C_Calc with last term sent = 1
	Twetbulbnatural_C = HeatMetrics_Calc::Twetbulb_C_Calc(Tair_K, RelativeHumidity_frac, AtmosphericPressure_hPa, WindSpeed_mps, Radiation_Shortwave_Wpm2, Radiation_Shortwave_direct_frac, Cos_ZenithAngle_Solar, 1);

	//Twetbulbnatural_C = min(Twetbulbnatural_C, Tair_C); Twetbulbnatural_C should never exceed Tair_C
	Twetbulbnatural_C = min(Twetbulbnatural_C, Tair_C);

	//Twetbulbpsychrometric_C computed in Twetbulb_C_Calc with last term sent = 0
	Twetbulbpsychrometric_C = HeatMetrics_Calc::Twetbulb_C_Calc(Tair_K, RelativeHumidity_frac, AtmosphericPressure_hPa, WindSpeed_mps, Radiation_Shortwave_Wpm2, Radiation_Shortwave_direct_frac, Cos_ZenithAngle_Solar, 0);

	//Twetbulbpsychrometric_C = min(Twetbulbpsychrometric_C , Tair_C); Twetbulbpsychrometric_C should never exceed Tair_C
	Twetbulbpsychrometric_C = min(Twetbulbpsychrometric_C, Tair_C);

	//Twetbulbglobe_C = 0.7 * Twetbulbnatural_C + 0.2 * Tglobe_C + 0.1 * Tair_C
	Twetbulbglobe_C = 0.7 * Twetbulbnatural_C + 0.2 * Tglobe_C + 0.1 * Tair_C;

	//UTCI_C returned from call to UniversalThermalClimateIndex_C_Calc after Tglobe_C is determined, which is needed for thermal mean radiant temperature 
	//Note: UniversalThermalClimateIndex_K or UTCI is additional heat metric; efficiently called from within WetBulbGlobeTemperature_C_Calc
	UTCI_C = UniversalThermalClimateIndex_C_Calc(input, MapPixel_ID, timeStep, Tair_C, WindSpeed_mps, VaporPressure_Actual_kPa, Tglobe_C);
}

//Tglobe_C_Calc function
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
double HeatMetrics_Calc::Tglobe_C_Calc(Inputs* input, double Tair_K, double RelativeHumidity_frac, double AtmosphericPressure_hPa, double WindSpeed_mps, double Radiation_Shortwave_Wpm2, double Radiation_Shortwave_direct_frac, double Cos_ZenithAngle_Solar, int MapPixel_ID, int timeStep) {

	//Tsfc_K initialized as Tair_K
	double Tsfc_K = Tair_K;
	//Tglobe_prev_K initialized as Tair_K
	double Tglobe_prev_K = Tair_K;
	//Flag_converged = false
	bool Flag_converged = false;
	int Iteration_Count = 0;

	//while (!Flag_converged && Iteration_Count < Iterations_Maximum, iterative loop
	while (!Flag_converged && Iteration_Count < Iterations_Maximum) {

		//Iteration_Count is advanced
		Iteration_Count++;
		//Tref_K = 0.5 * (Tglobe_prev_K + Tair_K)
		double Tref_K = 0.5 * (Tglobe_prev_K + Tair_K);

		//Coefficient_ConvectiveHeatTransfer computed from function HeatTransfer_Air_Sphere_Wpm2K_Calc from Eq 10 and following paragraph in Liljegren et al. (2008)
		double Coefficient_ConvectiveHeatTransfer = HeatMetrics_Calc::HeatTransfer_Air_Sphere_Wpm2K_Calc(Globe_Diameter_m, Tref_K, AtmosphericPressure_hPa, WindSpeed_mps);

		//Tglobe_new_K from Eq 17 in Liljegren et al. (2008), rearranged to solve for Tglobe_new_K by taking power 0.25 of RHS
		//Note: Eq 17 uses Radiation_Shortwave_direct_frac with Globe_Albedo_frac
		Tglobe_new_K = pow((0.5 * (HeatMetrics_Calc::Emissivity_Atmosphere_Calc(Tair_K, RelativeHumidity_frac) * pow(Tair_K, 4) + SurfaceGround_Emissivity_frac * pow(Tsfc_K, 4)) - Coefficient_ConvectiveHeatTransfer / (Sigma_Stefan_Boltzmann * Globe_Emissivity_frac) * (Tglobe_prev_K - Tair_K) + Radiation_Shortwave_Wpm2 / (2 * Sigma_Stefan_Boltzmann * Globe_Emissivity_frac) * (1 - Globe_Albedo_frac) * (1 + (1 / (2 * Cos_ZenithAngle_Solar) - 1) * Radiation_Shortwave_direct_frac + SurfaceGround_Albedo_frac)), 0.25);

		//If abs(Tglobe_new_K - Tglobe_prev_K) < Tolerance_Convergence then converged
		if (abs(Tglobe_new_K - Tglobe_prev_K) < Tolerance_Convergence) {
			Flag_converged = true;
		}

		//Tglobe_prev_K = 0.9 * Tglobe_prev_K + 0.1 * Tglobe_new_K, based on relaxation weights (0.9 and 0.1)
		//Note: Liljegren et al. (2008) describe this as an iterative, relaxation solution
		Tglobe_prev_K = 0.9 * Tglobe_prev_K + 0.1 * Tglobe_new_K;
	}
	//Tglobe_new_C (C) is Tglobe_new_K - 273.15
	double Tglobe_new_C = Tglobe_new_K - 273.15;

	//If Flag_converged = false then alert user
	if (Flag_converged == false) {
		cout << "Warning: The i-Tree Cool Air model cannot converge on a globe temperature" << endl;
		cout << " ... which is currently " << Tglobe_new_C << " C after " << Iterations_Maximum << " iterations." << endl;
		cout << "Aborting: This warning triggers the HydroPlus simulation to abort." << endl;
		cout << "Explanation: This condition suggests errors with input maps, parameters, calculated resistances, etc.." << endl;
		cout << "Explanation: This occurred at timestep " << timeStep << " and date hour " << input->SimulationDate_GDH[timeStep] << endl;
		cout << "Explanation: This occurred at MapPixel_ID " << MapPixel_ID << endl;
		cout << "Explanation: LandCover_Class " << input->LandCover_NLCD_Class[MapPixel_ID] << endl;
		cout << "Explanation: TreeCover_frac " << input->TreeCover_frac[MapPixel_ID] << endl;
		cout << "Explanation: ImperviousCover_frac " << input->ImperviousCover_frac[MapPixel_ID] << endl;
		cout << "Explanation: DEM_m " << input->Elevation_DEM_m[MapPixel_ID] << endl;
		cout << "Explanation: Tair_K " << Tair_K << endl;
		cout << "Explanation: RelativeHumidity_frac " << RelativeHumidity_frac << endl;
		cout << "Explanation: AtmosphericPressure_hPa " << AtmosphericPressure_hPa << endl;
		cout << "Explanation: WindSpeed_mps " << WindSpeed_mps << endl;
		cout << "Explanation: Radiation_Shortwave_Wpm2 " << Radiation_Shortwave_Wpm2 << endl;
		cout << "Explanation: Radiation_Shortwave_direct_frac " << Radiation_Shortwave_direct_frac << endl;
		cout << "Explanation: Cos_ZenithAngle_Solar " << Cos_ZenithAngle_Solar << endl;
		cout << "Correction: Improve inputs or modify code to allow for such large Tair variation." << endl;
		//Call abort function, which ends the HydroPlus.exe simulation
		abort();
	}

	//Return Tglobe_new_C
	return Tglobe_new_C;
}

//Twetbulb_C_Calc function
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
double HeatMetrics_Calc::Twetbulb_C_Calc(double Tair_K, double RelativeHumidity_frac, double AtmosphericPressure_hPa, double WindSpeed_mps,
	double Radiation_Shortwave_Wpm2, double Radiation_Shortwave_direct_frac, double Cos_ZenithAngle_Solar, int rad_flag) {

	//Tsfc_K initialized as Tair_K
	double Tsfc_K = Tair_K;
	//ZenithAngle_rad = acos(Cos_ZenithAngle_Solar)
	double ZenithAngle_rad = acos(Cos_ZenithAngle_Solar);
	//VaporPres_actual_mb = RelativeHumidity_frac * HeatMetrics_Calc::VaporPressure_haPa_Calc(Tair_K, 0)
	double VaporPres_actual_mb = RelativeHumidity_frac * HeatMetrics_Calc::VaporPressure_haPa_Calc(Tair_K, 0);

	//Tdew_K = Tdew_C + 273.15, where Tdew_C was defined as beC->by_key(MapPixel_ID, DataFolder_ID, "Tdew_K") in HeatMetrics_Calc::WetBulbGlobeTemperature_C_Calc
	//Note: Liljegren et al. (2008) computed Tdew_K with HeatMetrics_Calc::Dewpoint_Temperature_K_Calc(VaporPres_actual_mb, 0)
	double Tdew_K = Tdew_C + 273.15; 

	//Twetbulb_prev_K initialized as Tdew_K, but will be calculated given these are not equivalent terms
	double Twetbulb_prev_K = Tdew_K;
	//Flag_converged = false
	bool Flag_converged = false;
	int Iteration_Count = 0;

	//while (!Flag_converged && Iteration_Count < Iterations_Maximum), iterative loop
	while (!Flag_converged && Iteration_Count < Iterations_Maximum) {

		//Iteration_Count increased
		Iteration_Count++;
		//Tref_K = 0.5 * (Twetbulb_prev_K + Tair_K)
		double Tref_K = 0.5 * (Twetbulb_prev_K + Tair_K);

		//DeltaFnet_HeatFluxToWick_K_p_m2 from Eq 12 in Liljegren et al. (2008)
		double DeltaFnet_HeatFluxToWick_K_p_m2 = (Sigma_Stefan_Boltzmann * Wick_Emissivity_frac *
			(0.5 * (HeatMetrics_Calc::Emissivity_Atmosphere_Calc(Tair_K, RelativeHumidity_frac) * pow(Tair_K, 4) +
				SurfaceGround_Emissivity_frac * pow(Tsfc_K, 4)) -
				pow(Twetbulb_prev_K, 4)) +
			(1 - Wick_Albedo_frac) * Radiation_Shortwave_Wpm2 *
			((1 - Radiation_Shortwave_direct_frac) *
				(1 + 0.25 * Wick_Diameter_m / Wick_Length_m) +
				Radiation_Shortwave_direct_frac *
				((tan(ZenithAngle_rad) / M_PI) +
					0.25 * Wick_Diameter_m / Wick_Length_m) +
				SurfaceGround_Albedo_frac));

		//VaporPressure_wick_hPa (hPa or mb) from HeatMetrics_Calc::VaporPressure_haPa_Calc
		double VaporPressure_wick_hPa = HeatMetrics_Calc::VaporPressure_haPa_Calc(Twetbulb_prev_K, 0);
		//Air_Density_kg_p_m3 = AtmosphericPressure_hPa * 100 / (287.05 * Tref_K) using Ideal gas law
		double Air_Density_kg_p_m3 = AtmosphericPressure_hPa * 100 / (287.05 * Tref_K); 
		//Schmidt_Number_mass from Viscosity_Air_Pa_s_Calc function / (Air_Density_kg_p_m3 * Diffusivity_WaterVapor_m2ps_Calc function)
		double Schmidt_Number_mass = HeatMetrics_Calc::Viscosity_Air_Pa_s_Calc(Tref_K) / (Air_Density_kg_p_m3 * HeatMetrics_Calc::Diffusivity_WaterVapor_m2ps_Calc(Tref_K, AtmosphericPressure_hPa));
		//Coefficient_ConvectiveHeatTransfer from HeatMetrics_Calc::HeatTransfer_Air_Cylinder_Wpm2K_Calc function
		double Coefficient_ConvectiveHeatTransfer = HeatMetrics_Calc::HeatTransfer_Air_Cylinder_Wpm2K_Calc(Wick_Diameter_m, Wick_Length_m, Tref_K, AtmosphericPressure_hPa, WindSpeed_mps);

		//Twetbulb_new_K from Eq 9 in Liljegren et al. (2008)
		Twetbulb_new_K = Tair_K -
			HeatMetrics_Calc::LatentHeat_Vaporization_JpkgK_Calc(Tref_K) / Cp_Ratio_DryAir_to_WaterVapor *
			(VaporPressure_wick_hPa - VaporPres_actual_mb) /
			(AtmosphericPressure_hPa - VaporPressure_wick_hPa) *
			pow(Prandtl_Number_wind / Schmidt_Number_mass, 0.56) +
			DeltaFnet_HeatFluxToWick_K_p_m2 / Coefficient_ConvectiveHeatTransfer * rad_flag;

		//if (abs(Twetbulb_new_K - Twetbulb_prev_K) < Tolerance_Convergence then converged
		if (abs(Twetbulb_new_K - Twetbulb_prev_K) < Tolerance_Convergence) {
			Flag_converged = true;
		}

		//Twetbulb_prev_K = 0.9 * Twetbulb_prev_K + 0.1 * Twetbulb_new_K, updating with relaxation weights
		Twetbulb_prev_K = 0.9 * Twetbulb_prev_K + 0.1 * Twetbulb_new_K;
	}

	//Twetbulb_new_C = Twetbulb_new_K - 273.15
	Twetbulb_new_C = Twetbulb_new_K - 273.15;

	//If Flag_converged = false then alert user
	if (Flag_converged = false) {
		cout << "Warning: Convergence not achieved for calculation of Twetbulb_new_K" << endl;
	}

	//Return Twetbulb_new_C
	return Twetbulb_new_C;
}

//StabilityClass_Estimator to estimate the stability class based on Flag_daytime, wind WindSpeed_mps, Radiation_Shortwave_Wpm2 irradiance, and temperature difference
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Note: Reference: EPA (2000) section 6.2.5
int HeatMetrics_Calc::StabilityClass_Estimator(bool Flag_daytime, double WindSpeed_mps, double Radiation_Shortwave_Wpm2, double Delta_Height_Tair_m) {

	//Pasquill_Gifford_StabilityClasses_array is stability class lookup table
	//Note: The Pasquill - Gifford stability classes range from Value A (1) to F (6) or G (7), based on solar radiation, cloud cover, and wind speed
	//Note: Numbers, not letters, represent atmospheric stability classes in this application
	//Note: 1 indicating very unstable conditions and 6 or 7 representing very stable conditions
	//Note: Numbers taken from array based on identifying column j and row i, using 0-based indexing
	//Note: Column j = 4 (5th from left) and j = 7 (last) are not used in calculations
	const vector<vector<int>> Pasquill_Gifford_StabilityClasses_array = {
		{1, 1, 2, 4, 0, 5, 6, 0},
		{1, 2, 3, 4, 0, 5, 6, 0},
		{2, 2, 3, 4, 0, 4, 4, 0},
		{3, 3, 4, 4, 0, 0, 0, 0},
		{3, 4, 4, 4, 0, 0, 0, 0},
		{0, 0, 0, 0, 0, 0, 0, 0}
	};

	//Initialize j column and i row to 0
	int j = 0, i = 0;

	//Determine the column index j based on Flag_daytime conditions and Radiation_Shortwave_Wpm2 irradiance
	//If Flag_daytime then use Radiation_Shortwave_Wpm2 to determine 4 possible j column indices
	if (Flag_daytime) {
		if (Radiation_Shortwave_Wpm2 >= 925.0) {
			//j is column index
			j = 0;
		}
		else if (Radiation_Shortwave_Wpm2 >= 675.0) {
			j = 1;
		}
		else if (Radiation_Shortwave_Wpm2 >= 175.0) {
			j = 2;
		}
		else {
			j = 3;
		}
	}
	//Else not Flag_daytime then use Delta_Height_Tair_m to determine 2 other j column indices
	else {
		j = (Delta_Height_Tair_m >= 0.0) ? 6 : 5;
	}

	// Determine the row index `i` based on wind speed
	//If Flag_daytime then use WindSpeed_mps to determine 5 i row indices
	if (Flag_daytime) {
		if (WindSpeed_mps >= 6.0) {
			//i is row index 
			i = 4;
		}
		else if (WindSpeed_mps >= 5.0) {
			i = 3;
		}
		else if (WindSpeed_mps >= 3.0) {
			i = 2;
		}
		else if (WindSpeed_mps >= 2.0) {
			i = 1;
		}
		else {
			i = 0;
		}
	}
	//Else not Flag_daytime then use WindSpeed_mps to determine 3 i row indices
	else {
		if (WindSpeed_mps >= 2.5) {
			i = 2;
		}
		else if (WindSpeed_mps >= 2.0) {
			i = 1;
		}
		else {
			i = 0;
		}
	}

	//Return the stability class from the lookup table Pasquill_Gifford_StabilityClasses_array
	return Pasquill_Gifford_StabilityClasses_array[i][j];
}

//WindSpeed_Estimate_mps_Calc to calculate the wind WindSpeed_mps at a height for Twbg calculations using stability class and Flag_urban/rural factors.
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Note: Reference: EPA-454/5-99-005, 2000, section 6.2.5
//Reference: U.S. Environmental Protection Agency (EPA): Meteorological data processing. In Meteorological Monitoring Guidance for Regulatory Modeling Applications (EPA-454/R-99-005). Washington, D.C.: EPA, 2000
double HeatMetrics_Calc::WindSpeed_Estimate_mps_Calc(double WindSpeed_mps, double Height_WindSpeed_Actual_m, double Height_WindSpeed_Target_m, int StabilityClass_Atmosphere_int, bool Flag_urban) {

	//Define the Exponent_WindSpeed_AdjustForHeight_Eqn arrays for Flag_urban and rural areas
	const vector<double> urban_exp = { 0.15, 0.15, 0.20, 0.25, 0.30, 0.30 };
	const vector<double> rural_exp = { 0.07, 0.07, 0.10, 0.15, 0.35, 0.55 };

	//Choose the appropriate Exponent_WindSpeed_AdjustForHeight_Eqn based on Flag_urban or rural setting
	double Exponent_WindSpeed_AdjustForHeight_Eqn = Flag_urban
		? urban_exp[StabilityClass_Atmosphere_int - 1]
		: rural_exp[StabilityClass_Atmosphere_int - 1];

	//Estimate the wind WindSpeed_mps adjusted for the reference height
	//WindSpeed_estimated_mps (m/s) computed in WindSpeed_Estimate_mps_Calc to adjust from measured height to Twbg height
	//Note: Height_WindSpeed_Actual_m = 10 m for most NOAA weather stations, while Height_WindSpeed_Reference_m = 2 m, so lower wind speed
	double WindSpeed_estimated_mps = WindSpeed_mps * pow(Height_WindSpeed_Reference_m / Height_WindSpeed_Actual_m, Exponent_WindSpeed_AdjustForHeight_Eqn);

	//Ensure WindSpeed_mps does not fall below the minimum WindSpeed_mps
	WindSpeed_estimated_mps = max(WindSpeed_estimated_mps, WindSpeed_Minimum_mps);

	return WindSpeed_estimated_mps;
}

//HeatTransfer_Air_Sphere_Wpm2K_Calc to calculate the convective heat transfer coefficient for flow around a sphere
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
double HeatMetrics_Calc::HeatTransfer_Air_Sphere_Wpm2K_Calc(double Sphere_Diameter_m, double Tair_K, double AtmosphericPressure_hPa, double WindSpeed_mps) {

	//Air_Density_kg_p_m3 = AtmosphericPressure_hPa * 100 / (Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K * Tair_K), using Ideal gas law
	//Note: *100 converts hPa to Pa, ensures correct units
	double Air_Density_kg_p_m3 = AtmosphericPressure_hPa * 100 / (Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K * Tair_K);

	//Calculate Reynolds number, ensuring speed is above the minimum using Viscosity_Air_Pa_s_Calc function
	double Reynolds_Number_wind = MAX(WindSpeed_mps, WindSpeed_Minimum_mps) * Air_Density_kg_p_m3 * Sphere_Diameter_m / HeatMetrics_Calc::Viscosity_Air_Pa_s_Calc(Tair_K);

	//Nusselt_Number_wind = 2.0 + 0.6 * sqrt(Reynolds_Number_wind) * pow(Prandtl_Number_wind, 0.3333) as Eq 16 from Liljegren et al (2008)
	double Nusselt_Number_wind = 2.0 + 0.6 * sqrt(Reynolds_Number_wind) * pow(Prandtl_Number_wind, 0.3333);

	//HeatTransfer_Convective_Wpm2K W/(m2·K) = Nusselt_Number_wind * HeatMetrics_Calc::ThermalConductivity_DryAir_WpmK_Calc(Tair_K) / Sphere_Diameter_m
	double HeatTransfer_Wpm2K = Nusselt_Number_wind * HeatMetrics_Calc::ThermalConductivity_DryAir_WpmK_Calc(Tair_K) / Sphere_Diameter_m;

	//return HeatTransfer_Convective_Wpm2K 
	return HeatTransfer_Wpm2K;
}

//HeatTransfer_Air_Cylinder_Wpm2K_Calc: convective heat transfer coefficient calculation
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Reference: Bird, Stewart, and Lightfoot (BSL), page 409.
//Reference: Bird, R.B., W.E. Stewart, and E.N. Lightfoot: Transport Phenomena. New York: John Wiley & Sons, 1960. pp. 649–652.
//Coefficient_ConvectiveHeatTransfer from Eq 10 and following paragraph in Liljegren et al. (2008)
double HeatMetrics_Calc::HeatTransfer_Air_Cylinder_Wpm2K_Calc(double Diameter_cylinder_m, double Length_cylinder_m, double Tair_K, double AtmosphericPressure_hPa, double WindSpeed_mps) {

	//Air_Density_kg_p_m3 = AtmosphericPressure_hPa * 100 / (Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K * Tair_K)
	double Air_Density_kg_p_m3 = AtmosphericPressure_hPa * 100 / (Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K * Tair_K);

	//Reynolds_Number_wind = MAX(WindSpeed_mps, WindSpeed_Minimum_mps) * Air_Density_kg_p_m3 * Diameter_cylinder_m / HeatMetrics_Calc::Viscosity_Air_Pa_s_Calc(Tair_K)
	double Reynolds_Number_wind = MAX(WindSpeed_mps, WindSpeed_Minimum_mps) * Air_Density_kg_p_m3 * Diameter_cylinder_m / HeatMetrics_Calc::Viscosity_Air_Pa_s_Calc(Tair_K);

	//Nusselt number estimator from Eq 10 and following paragraph in Liljegren et al. (2008)
	//Note: Reynolds_Number_wind^0.6 from Nusselt_Number_wind * Reynolds_Number_wind^(-1) * Prandtl_Number_wind^(a-1) = b * Reynolds_Number_wind^(-c), ..
	//Note: ... where a = 0.56, b = 0.281, c = 0.4, hence combine Reynolds_Number_wind terms on RHS, = Reynolds_Number_wind^(1-c) = Reynolds_Number_wind^0.6
	double Nusselt_Number_wind = 0.281 * pow(Reynolds_Number_wind, 0.6) * pow(Prandtl_Number_wind, 0.44);

	//HeatTransfer_Convective_Wpm2K W/(m2·K) = Nusselt_Number_wind * HeatMetrics_Calc::ThermalConductivity_DryAir_WpmK_Calc(Tair_K) / Sphere_Diameter_m
	double HeatTransfer_Convective_Wpm2K = Nusselt_Number_wind * HeatMetrics_Calc::ThermalConductivity_DryAir_WpmK_Calc(Tair_K) / Diameter_cylinder_m;

	//return HeatTransfer_Convective_Wpm2K
	return HeatTransfer_Convective_Wpm2K;
}

//Emissivity_Atmosphere_Calc to calculate the atmospheric emissivity
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Reference: Oke (2nd edition), page 373.
//Reference: Oke, T.R.: Boundary Layer Climates, 2nd Ed. New York: Methuen & Co., 1978. pp. 373.
double HeatMetrics_Calc::Emissivity_Atmosphere_Calc(double Tair_K, double RelativeHumidity_frac) {

	//VaporPres_hPa = RelativeHumidity_frac * HeatMetrics_Calc::VaporPressure_haPa_Calc(Tair_K, 0)
	double VaporPres_hPa = RelativeHumidity_frac * HeatMetrics_Calc::VaporPressure_haPa_Calc(Tair_K, 0);
	//Emissivity_Atmosphere_frac = 0.575 * pow(VaporPres_hPa, 0.143) using the formula from Oke (1978) page 373
	double Emissivity_Atmosphere_frac = 0.575 * pow(VaporPres_hPa, 0.143);

	//return Emissivity_Atmosphere_frac
	return Emissivity_Atmosphere_frac;
}

//VaporPressure_haPa_Calc to calculate saturation vapor pressure over water as liquid or ice
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Note: Difference of VaporPressure_mb values for Tair_C < 0 when using Flag_WaterPhase_is_Ice = 0 or 1 is ~ 10% at -5 C
//Note: Lilijegren et al. (2008) code used Buck's (1981) approximation (eqn 3) of Wexler's (1976) formulae ...
//Note: Lilijegren et al. (2008) code had several commented out sections to adjust for atmospheric pressure
//Note: HydroPlus removed the uncommented Lilijegren et al. (2008) code adjustment to correct for pressure given it had nominal effect
//Note: HydroPlus replaced the Lilijegren et al. (2008) vapor pressure Eqs with Eq 4.2.2 in Shuttleworth (1993) for code consistency
//Reference:
//Buck, A.L.: New equations for computing vapor pressure and enhancement factor. J. Appl. Meteorol. 20:1527–1532 (1981).
double HeatMetrics_Calc::VaporPressure_haPa_Calc(double Tair_K, int Flag_WaterPhase_is_Ice) {
	
	double VaporPres_saturated_haPa;
	//Tair_C = Tair_K - 273.15
	double Tair_C = Tair_K - 273.15;

	//If Flag_WaterPhase_is_Ice = 0 then saturation vapor pressure over water 
	if (Flag_WaterPhase_is_Ice == 0) {
		//VaporPres_saturated_haPa (mb) as Tetens' formula in Eq 4.2.2 Shuttleworth (1993), e_kPa=0.6108*exp[(17.27*T_C)/(237.3+T_C)]
		//Note: Coefficient 6.108 mb is replacement of 0.6108 kPa used in Shuttleworth (1993) 
		VaporPres_saturated_haPa = 6.108 * exp((17.27 * Tair_C) / (237.3 + Tair_C));
	}
	//Else If Flag_WaterPhase_is_Ice = 1 then saturation vapor pressure over ice
	else {
		//VaporPressure_mb (mb) as Tetens' over ice from Murray (1967) in Wikipedia, https://en.wikipedia.org/wiki/Tetens_equation
		VaporPres_saturated_haPa = 6.108 * exp((21.875 * Tair_C) / (265.5 + Tair_C));
	}

	//return VaporPres_saturated_haPa, note mb = haPa
	return VaporPres_saturated_haPa;
}

//Dewpoint_Temperature_K_Calc to calculate the dew point temperature over water (Flag_WaterPhase_is_Ice=0) or ice (Flag_WaterPhase_is_Ice=1).
double HeatMetrics_Calc::Dewpoint_Temperature_K_Calc(double VaporPres_mb, int Flag_WaterPhase_is_Ice) {
	double Tdew_C;
	//If Flag_WaterPhase_is_Ice is 0 then dewpoint temperature over water
	if (Flag_WaterPhase_is_Ice == 0) {
		//T_C = 237.3 * ln(e_mb / 6.108) / (17.27 - ln(e_mb / 6.108)) derived from Tetens' formula in Eq 4.2.2 Shuttleworth (1993) ...
		//Note: Step 1: e_mb = 6.108 * exp[(17.27 * T_C) / (237.3 + T_C)], for vapor pressure in mb, not kPa as in Shuttleworth (1993) ...
		//Note: Step 2: ln(e_mb / 6.108) = (17.27 * T_C) / (237.3 + T_C)
		//Note: Step 3: (237.3 + T_C) * ln(e_mb / 6.108) = (17.27 * T_C)
		//Note: Step 4: 237.3 * ln(e_mb / 6.108) + T_C * ln(e_mb / 6.108) = (17.27 * T_C)
		//Note: Step 5: 237.3 * ln(e_mb / 6.108) = (17.27 * T_C) - T_C * ln(e_mb / 6.108)
		//Note: Step 6: 237.3 * ln(e_mb / 6.108) = T_C * (17.27 - ln(e_mb / 6.108))
		//Note: Step 7: T_C = 237.3 * ln(e_mb / 6.108) / (17.27 - ln(e_mb / 6.108))
		Tdew_C = 237.3 * log(VaporPres_mb / 6.108) / (17.27 - log(VaporPres_mb / 6.108));
	}
	//Else dewpoint temperature over ice
	else {
		//T_C = (265.5 * ln(e_mb / 6.108)) / (21.875 - ln(e_mb / 6.108)) derived from Tetens' formula over ice ...
		//Note: Step 1: e_mb = 6.108 * exp((21.875 * T_C) / (265.5 + T_C)) from Murray (1967) in Wikipedia, https://en.wikipedia.org/wiki/Tetens_equation
		//Note: Step 2 to 7 follow the steps for over water
		Tdew_C = 265.5 * log(VaporPres_mb / 6.108) / (21.875 - log(VaporPres_mb / 6.108));
	}

	//Tdew_K = Tdew_C + 273.15
	double Tdew_K = Tdew_C + 273.15;
	//return Tdew_K
	return Tdew_K;
}

//ThermalConductivity_DryAir_WpmK_Calc to calculate the thermal conductivity of dry air, W/(m K)
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Note: Bird, Stewart, and Lightfoot (BSL), page 257.
//Reference: 
//Bird, R.B., W.E. Stewart, and E.N. Lightfoot: Transport Phenomena. New York: John Wiley & Sons, 1960. pp. 649–652.
double HeatMetrics_Calc::ThermalConductivity_DryAir_WpmK_Calc(double Tair_K) {

	//ThermalConductivity_DryAir_W_p_m_K (W/(m·K)) is (SpecificHeat_DryAir_J_p_kg_K + 1.25 * Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K) * HeatMetrics_Calc::Viscosity_Air_Pa_s_Calc(Tair_K)
	double ThermalConductivity_DryAir_W_p_m_K = (SpecificHeat_DryAir_J_p_kg_K + 1.25 * Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K) * HeatMetrics_Calc::Viscosity_Air_Pa_s_Calc(Tair_K);
	
	//return ThermalConductivity_DryAir_W_p_m_K
	return ThermalConductivity_DryAir_W_p_m_K;
}

//SpecificHeat_HumidAir_JpkgK_Calc to calculate the specific heat of air of any humidity, J/(kg K)
double HeatMetrics_Calc::SpecificHeat_HumidAir_JpkgK_Calc(double Tair_K, double Tdew_K, double AtmPres_kpa) {

	//Tair_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tair_C = Tair_K - 273.15;
	//Tdew_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tdew_C = Tdew_K - 273.15;

	//SpecificHeat_DryAir_Tair_J_p_kg_K (J/kg K) is 0.0002 * pow(Tair_C, 2) + 0.033 * Tair_C + 1003.8, capturing variation with Tair, based on:
	//Note: Data points to derive Eq from Bohren and Albrecht, Smithsonian Table 256 within Stull (2017) Figure 3I.1 Emperical Estimates of Cpd
	//Note: Data Vectors: Tair_C: -100,-30,0,20,40,60,100; Cpd_J_p_kg_K: 1002.5,1002.99,1003.8,1004.54,1005.44,1006.5,1009.1
	//Note: See Info Specific Heat Cp for Air. Emperical Estimates of Cpd, Stull (2017), Practical Meteorology, Chpt. Thermodynamics (p56)
	double SpecificHeat_DryAir_Tair_J_p_kg_K = 0.0002 * pow(Tair_C, 2) + 0.033 * Tair_C + 1003.8;

	int Flag_WaterPhase_is_Ice = 0;
	//VaporPres_actual_mb = HeatMetrics_Calc::VaporPressure_haPa_Calc(Tdew_K, Flag_WaterPhase_is_Ice), sending Tdew_K and 0 as flag indicating liquid
	double VaporPres_actual_mb = HeatMetrics_Calc::VaporPressure_haPa_Calc(Tdew_K, Flag_WaterPhase_is_Ice);
	//VaporPres_actual_kPa = VaporPres_actual_mb * 0.1, converting from mb or hPa to kPa
	double VaporPres_actual_kPa = VaporPres_actual_mb * 0.1; 
	//MixingRatio_g_p_g is (Ratio_MolecularWeightVapor_to_MolecularWeightDryAir * VaporPres_actual_kPa) / (AtmPres_kpa - VaporPres_actual_kPa)
	//Note: From Eq 4.4 of Stull (2017)
	double MixingRatio_g_p_g = (Ratio_MolecularWeightVapor_to_MolecularWeightDryAir * VaporPres_actual_kPa) / (AtmPres_kpa - VaporPres_actual_kPa);

	//SpecificHeat_HumidAir_J_p_kg_K (W/(m·K)) is SpecificHeat_DryAir_Tair_J_p_kg_K * (1 + 1.84 * MixingRatio_g_p_g)
	//Note: From Eq 3.3 of Stull (2017)
	double SpecificHeat_HumidAir_J_p_kg_K = SpecificHeat_DryAir_Tair_J_p_kg_K * (1 + 1.84 * MixingRatio_g_p_g);

	//return SpecificHeat_HumidAir_J_p_kg_K
	return SpecificHeat_HumidAir_J_p_kg_K;
}

//LatentHeat_Vaporization_JpkgK_Calc to calculate the heat of evaporation, J/(kg K), for temperature in the range 283-313 K.
//Note: Van Wylen and Sonntag (1986), Table A.1.1
//Reference: 
//Van Wylen and Sonntag. (1986) Fundamentals of Classical Thermodynamics, 3rd ed., John Wiley & Sons, New York.
double HeatMetrics_Calc::LatentHeat_Vaporization_JpkgK_Calc(double Tair_K) {
	//LatentHeat_Evaporation_J_p_kg_K (J/kg/K) is (313.15 - Tair_K) / 30.0 * -71100.0 + 2.4073E6
	double LatentHeat_Evaporation_J_p_kg_K = (313.15 - Tair_K) / 30.0 * -71100.0 + 2.4073E6;
	return LatentHeat_Evaporation_J_p_kg_K;
}

//Viscosity_Air_Pa_s_Calc to calculate viscosity
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Reference: Bird, Stewart, and Lightfoot (BSL), page 23
//Reference: Bird, R.B., W.E. Stewart, and E.N. Lightfoot: Transport Phenomena. New York: John Wiley & Sons, 1960. pp. 649–652.
double HeatMetrics_Calc::Viscosity_Air_Pa_s_Calc(double Tair_K) {

	//Viscosity_Air_Pa_s is 2.6693E-6 * sqrt(Molecular_Weight_DryAir_kg_p_kmol * Tair_K) / (3.617 * 3.617 * (((Tair_K / 97.0 - 2.9) / 0.4) * -0.034 + 1.048))
	double Viscosity_Air_Pa_s = 2.6693E-6 * sqrt(Molecular_Weight_DryAir_kg_p_kmol * Tair_K) / (3.617 * 3.617 * (((Tair_K / 97.0 - 2.9) / 0.4) * -0.034 + 1.048));
	//return Viscosity_Air_Pa_s
	return Viscosity_Air_Pa_s;
}

//Diffusivity_WaterVapor_m2ps_Calc to calculate the diffusivity of water vapor in air, m²/s
//Note: Part of WetBulbGlobeTemperature algorithm modified from Liljegren et al. (2008) code: https://github.com/mdljts/wbgt/blob/master/src/wbgt.c
//Note: Bird, Stewart, and Lightfoot (BSL), page 505.
//Reference: 
//Bird, R.B., W.E. Stewart, and E.N. Lightfoot: Transport Phenomena. New York: John Wiley & Sons, 1960. pp. 649–652.
double HeatMetrics_Calc::Diffusivity_WaterVapor_m2ps_Calc(double Tair_K, double AtmosphericPressure_hPa) {
	
	//Critical properties of air and water
	const double Pcrit_air = 36.4, Pcrit_h2o = 218.0; // Critical pressures (bar)
	const double Tcrit_air = 132.0, Tcrit_h2o = 647.3; // Critical temperatures (K)
	const double Coeff_a = 3.640E-4, Power_b = 2.334; // Coefficients for Diffusivity_WaterVapor_m2ps_Calc calculation

	//Intermediate calculations
	double Pcrit_air_water_1p3 = pow(Pcrit_air * Pcrit_h2o, 1.0 / 3.0);
	double Tcrit_air_water_5p12 = pow(Tcrit_air * Tcrit_h2o, 5.0 / 12.0);
	double Tcrit_air_water_1p2 = sqrt(Tcrit_air * Tcrit_h2o);
	double MolecularWeightMix_air_water_1p2 = sqrt(1.0 / Molecular_Weight_DryAir_kg_p_kmol + 1.0 / Molecular_Weight_WaterVapor_kg_p_kmol);
	double Ratio_AtmPressure_Actual_to_SeaLevelConstant = AtmosphericPressure_hPa / 1013.25;

	//Diffusivity_WaterVapor_m2_p_s (m²/s) is (Coeff_a * pow(Tair_K / Tcrit_air_water_1p2, Power_b) * Pcrit_air_water_1p3 * Tcrit_air_water_5p12 * MolecularWeightMix_air_water_1p2 / Ratio_AtmPressure_Actual_to_SeaLevelConstant * 1E-4)
	double Diffusivity_WaterVapor_m2_p_s = (Coeff_a * pow(Tair_K / Tcrit_air_water_1p2, Power_b) * Pcrit_air_water_1p3 * Tcrit_air_water_5p12 * MolecularWeightMix_air_water_1p2 / Ratio_AtmPressure_Actual_to_SeaLevelConstant * 1E-4);

	//return Diffusivity_WaterVapor_m2_p_s
	return Diffusivity_WaterVapor_m2_p_s;
}

//Density_Water_kgpm3_Calc to calculate the density of liquid sater, kg/m3
//Note: From Equation in Figure 11.1.1 McCutcheon (1993)
//Reference: 
//McCutcheon, S.C., Martin, J. and T.O. Barnwell Jr. (1993). Water Quality. In D. R. Maidment (Ed.), Handbook of Hydrology (pp. 11.1-11.73). New York: McGraw-Hill.
double HeatMetrics_Calc::Density_Water_kgpm3_Calc(double Tair_K) {

	//Temperature_Air_C = Tair_K - 273.15
	double Temperature_Air_C = Tair_K - 273.15;

	//Density_Water_kgpm3 is density of water (kg/m^3) from Equation in Figure 11.1.1 McCutcheon (1993)
	double Density_Water_kgpm3 = 1000 * (1 - ((Temperature_Air_C + 288.9414) / (508929.2 * (Temperature_Air_C + 68.12963))) *
		pow((Temperature_Air_C - 3.9863), 2));

	//return Density_Water_kgpm3
	return Density_Water_kgpm3;
}

//UniversalThermalClimateIndex_C_Calc function will manage calculation of the UTCI temperature and return value
//Note: Consider clarification if function needs wind speed at 10 m or height of Tglobe_C and Tmrt, which is 2 m height
double HeatMetrics_Calc::UniversalThermalClimateIndex_C_Calc(Inputs* input, int MapPixel_ID, int timeStep, double Tair_C, double WindSpeed_mps, double VaporPressure_Actual_kPa, double Tglobe_C) {
	//Note: UTCI from https://www.utci.org/utci_calc.php and COST Action 730 UTCI UTCI_a002.f90 code, written by Peter Broede et al. (2012)
	//Note: COST Action 730 UTCI calculation uses 6th order polynomial with meteorological inputs including Thermal mean radiant temperature (Tmrt)
	//Note: Thermal mean radiant temperature computed with Eq 2 forced convection of Li et al. (2022), which is the same Eq used in degree F at ...
	//Note: ... https://www.engineersedge.com/calculators/mean_radiant_temperature_15726.htm 
	//Note: UTCI sensitivity discussed at: https://compute.procedural.build/docs/comfort/utci_sensitivity, which references Broede et al. (2012)
	//References: 
	//Bröde, P., Fiala, D., Błażejczyk, K., Holmér, I., Jendritzky, G., Kampmann, B., et al. (2012). Deriving the operational procedure for the Universal Thermal Climate Index (UTCI). International Journal of Biometeorology, 56(3), 481-494. https://doi.org/10.1007/s00484-011-0454-1
	//Li, Z., Feng, X., & Fang, Z. (2022). A modified method to measure outdoor mean radiant temperature: Comparison between two-hemisphere method and six-direction integral method. Building and Environment, 221, 109292. https://www.sciencedirect.com/science/article/pii/S0360132322005273

	//ThermalMeanRadiantTemp_C is pow(pow(Tglobe_C + 273.15, 4) + (1.1E8 * pow(WindSpeed_mps, 0.6)) / (Globe_Emissivity_frac * pow(Globe_Diameter_m, 0.4)) * (Tglobe_C - Tair_C), 0.25) - 273.15
	//Note: Based on Eq 2 of Li et al. (2012); equivalent to https://www.engineersedge.com/calculators/mean_radiant_temperature_15726.htm 
	ThermalMeanRadiantTemp_C = pow(pow(Tglobe_C + 273.15, 4) + (1.1E8 * pow(WindSpeed_mps, 0.6)) / (Globe_Emissivity_frac * pow(Globe_Diameter_m, 0.4)) * (Tglobe_C - Tair_C), 0.25) - 273.15;

	//Delta_Tmrt_Tair_C = ThermalMeanRadiantTemp_C - Tair_C; term within the UTCI calculation
	double Delta_Tmrt_Tair_C = ThermalMeanRadiantTemp_C - Tair_C;

	//UniversalThermalClimateIndex_C is a 6th order polynomial from COST Action 730 UTCI UTCI_a002.f90 code, written by Peter Broede et al. (2012)
	//Note: The equation has 210 coefficients (!), and was copied directly the source code, tested against https://www.utci.org/utci_calc.php
	UniversalThermalClimateIndex_C = Tair_C +
		(6.07562052E-01) + 
		(-2.27712343E-2) * Tair_C + 
		(8.06470249E-04) * Tair_C * Tair_C + 
		(-1.54271372E-04) * Tair_C * Tair_C * Tair_C + 
		(-3.24651735E-06) * Tair_C * Tair_C * Tair_C * Tair_C + 
		(7.32602852E-08) * Tair_C * Tair_C * Tair_C * Tair_C * Tair_C + 
		(1.35959073E-09) * Tair_C * Tair_C * Tair_C * Tair_C * Tair_C * Tair_C + 
		(-2.25836520E+00) * WindSpeed_mps + 
		(8.80326035E-02) * Tair_C * WindSpeed_mps + 
		(2.16844454E-03) * Tair_C * Tair_C * WindSpeed_mps + 
		(-1.53347087E-05) * Tair_C * Tair_C * Tair_C * WindSpeed_mps + 
		(-5.72983704E-07) * Tair_C * Tair_C * Tair_C * Tair_C * WindSpeed_mps + 
		(-2.55090145E-09) * Tair_C * Tair_C * Tair_C * Tair_C * Tair_C * WindSpeed_mps + 
		(-7.51269505E-01) * WindSpeed_mps * WindSpeed_mps + 
		(-4.08350271E-03) * Tair_C * WindSpeed_mps * WindSpeed_mps + 
		(-5.21670675E-05) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps + 
		(1.94544667E-06) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps + 
		(1.14099531E-08) * Tair_C * Tair_C * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps + 
		(1.58137256E-01) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(-6.57263143E-05) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(2.22697524E-07) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(-4.16117031E-08) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(-1.27762753E-02) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(9.66891875E-06) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(2.52785852E-09) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(4.56306672E-04) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(-1.74202546E-07) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(-5.91491269E-06) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps + 
		(3.98374029E-01) * Delta_Tmrt_Tair_C + 
		(1.83945314E-04) * Tair_C * Delta_Tmrt_Tair_C + 
		(-1.73754510E-04) * Tair_C * Tair_C * Delta_Tmrt_Tair_C + 
		(-7.60781159E-07) * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C + 
		(3.77830287E-08) * Tair_C * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C + 
		(5.43079673E-10) * Tair_C * Tair_C * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C + 
		(-2.00518269E-02) * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(8.92859837E-04) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(3.45433048E-06) * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(-3.77925774E-07) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(-1.69699377E-09) * Tair_C * Tair_C * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(1.69992415E-04) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(-4.99204314E-05) * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(2.47417178E-07) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(1.07596466E-08) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(8.49242932E-05) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(1.35191328E-06) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(-6.21531254E-09) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(-4.99410301E-06) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(-1.89489258E-08) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(8.15300114E-08) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C + 
		(7.55043090E-04) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-5.65095215E-05) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-4.52166564E-07) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(2.46688878E-08) * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(2.42674348E-10) * Tair_C * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(1.54547250E-04) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(5.24110970E-06) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-8.75874982E-08) * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-1.50743064E-09) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-1.56236307E-05) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-1.33895614E-07) * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(2.49709824E-09) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(6.51711721E-07) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(1.94960053E-09) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-1.00361113E-08) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-1.21206673E-05) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-2.18203660E-07) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(7.51269482E-09) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(9.79063848E-11) * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(1.25006734E-06) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-1.81584736E-09) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-3.52197671E-10) * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-3.36514630E-08) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(1.35908359E-10) * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(4.17032620E-10) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-1.30369025E-09) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(4.13908461E-10) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(9.22652254E-12) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-5.08220384E-09) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-2.24730961E-11) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(1.17139133E-10) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(6.62154879E-10) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(4.03863260E-13) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(1.95087203E-12) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(-4.73602469E-12) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C + 
		(5.12733497E+00) * VaporPressure_Actual_kPa + 
		(-3.12788561E-01) * Tair_C * VaporPressure_Actual_kPa + 
		(-1.96701861E-02) * Tair_C * Tair_C * VaporPressure_Actual_kPa + 
		(9.99690870E-04) * Tair_C * Tair_C * Tair_C * VaporPressure_Actual_kPa + 
		(9.51738512E-06) * Tair_C * Tair_C * Tair_C * Tair_C * VaporPressure_Actual_kPa + 
		(-4.66426341E-07) * Tair_C * Tair_C * Tair_C * Tair_C * Tair_C * VaporPressure_Actual_kPa + 
		(5.48050612E-01) * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-3.30552823E-03) * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-1.64119440E-03) * Tair_C * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-5.16670694E-06) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(9.52692432E-07) * Tair_C * Tair_C * Tair_C * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-4.29223622E-02) * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(5.00845667E-03) * Tair_C * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(1.00601257E-06) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-1.81748644E-06) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-1.25813502E-03) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-1.79330391E-04) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(2.34994441E-06) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(1.29735808E-04) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(1.29064870E-06) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-2.28558686E-06) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa + 
		(-3.69476348E-02) * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(1.62325322E-03) * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-3.14279680E-05) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(2.59835559E-06) * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-4.77136523E-08)* Tair_C* Tair_C* Tair_C* Tair_C* Delta_Tmrt_Tair_C* VaporPressure_Actual_kPa + 
		(8.64203390E-03) * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-6.87405181E-04) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-9.13863872E-06) * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(5.15916806E-07) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-3.59217476E-05) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(3.28696511E-05) * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-7.10542454E-07) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-1.24382300E-05) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-7.38584400E-09) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(2.20609296E-07) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-7.32469180E-04) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-1.87381964E-05) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(4.80925239E-06) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-8.75492040E-08) * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(2.77862930E-05) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-5.06004592E-06) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(1.14325367E-07) * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(2.53016723E-06) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-1.72857035E-08) * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-3.95079398E-08) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-3.59413173E-07) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(7.04388046E-07) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-1.89309167E-08) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-4.79768731E-07) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(7.96079978E-09) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(1.62897058E-09) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(3.94367674E-08) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-1.18566247E-09) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(3.34678041E-10) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-1.15606447E-10) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa + 
		(-2.80626406E+00) * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(5.48712484E-01) * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-3.99428410E-03) * Tair_C * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-9.54009191E-04) * Tair_C * Tair_C * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.93090978E-05) * Tair_C * Tair_C * Tair_C * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-3.08806365E-01) * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.16952364E-02) * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(4.95271903E-04) * Tair_C * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-1.90710882E-05) * Tair_C * Tair_C * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.10787756E-03) * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-6.98445738E-04) * Tair_C * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.30109073E-05) * Tair_C * Tair_C * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(4.17856590E-04) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-1.27043871E-05) * Tair_C * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-3.04620472E-06) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(5.14507424E-02) * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-4.32510997E-03) * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(8.99281156E-05) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-7.14663943E-07) * Tair_C * Tair_C * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-2.66016305E-04) * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.63789586E-04) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-7.01199003E-06) * Tair_C * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-1.06823306E-04) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(3.61341136E-06) * Tair_C * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.29748967E-07) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(3.04788893E-04) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-6.42070836E-05) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.16257971E-06) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(7.68023384E-06) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-5.47446896E-07) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-3.59937910E-08) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-4.36497725E-06) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.68737969E-07) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.67489271E-08) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(3.23926897E-09) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-3.53874123E-02) * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-2.21201190E-01) * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.55126038E-02) * Tair_C * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-2.63917279E-04) * Tair_C * Tair_C * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(4.53433455E-02) * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-4.32943862E-03) * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.45389826E-04) * Tair_C * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.17508610E-04) * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-6.66724702E-05) * Tair_C * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(3.33217140E-05) * WindSpeed_mps * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-2.26921615E-03) * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(3.80261982E-04) * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-5.45314314E-09) * Tair_C * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-7.96355448E-04) * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.53458034E-05) * Tair_C * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-6.31223658E-06) * WindSpeed_mps * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(3.02122035E-04) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-4.77403547E-06) * Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.73825715E-06) * WindSpeed_mps * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-4.09087898E-07) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(6.14155345E-01) * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-6.16755931E-02) * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.33374846E-03) * Tair_C * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(3.55375387E-03) * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-5.13027851E-04) * Tair_C * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.02449757E-04) * WindSpeed_mps * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-1.48526421E-03) * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-4.11469183E-05) * Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-6.80434415E-06) * WindSpeed_mps * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-9.77675906E-06) * Delta_Tmrt_Tair_C * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(8.82773108E-02) * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(-3.01859306E-03) * Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.04452989E-03) * WindSpeed_mps * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(2.47090539E-04) * Delta_Tmrt_Tair_C * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa + 
		(1.48348065E-03) * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa * VaporPressure_Actual_kPa;
	//End of polynomial
	
	//return UniversalThermalClimateIndex_C
	return UniversalThermalClimateIndex_C;

	/*
	//If conditional to print output for verification 
	//Note: UniversalThermalClimateIndex_C=53.911 when Tiar_C=36.68, Delta_Tmrt_Tair_C=65.7684, VaporPressure_Actual_hPa=22.4025 WindSpeed_mps=1.75573
	if (timeStep == 12 && MapPixel_ID == 0) {
		cout << "HeatMetrics: UniversalThermalClimateIndex_C " << UniversalThermalClimateIndex_C << endl;
		cout << "HeatMetrics: Tair_C " << Tair_C << " Delta_Tmrt_Tair_C " << Delta_Tmrt_Tair_C << endl;
		cout << "HeatMetrics: VaporPressure_Actual_hPa " << VaporPressure_Actual_kPa * 10 << " WindSpeed_mps " << WindSpeed_mps << endl;
		cout << "HeatMetrics: ThermalMeanRadiantTemp_C " << ThermalMeanRadiantTemp_C << " Tglobe_C " << Tglobe_C << endl;
	}
	*/
}

//WindChill_K_Calc function computes wind chill temperature using algorithm and equations from Eq 3.64a and 3.64b of Stull (2017)
//Note: Alternative The NOAA NWS Wind Chill Temperature Index updated in 2001, https://www.weather.gov/safety/cold-wind-chill-chart
//Note: Formula for wind chill is on bottom of chart, and windspeed is mph or miles per hour ...
//Note: WindChill_F = 35.74 + 0.6215(Tair_K) - [35.75 * (Wspd_miles_p_hour ^ 0.16)] + [(0.4275 * Tair_K * (Wspd_miles_p_hour ^ 0.16)]
void HeatMetrics_Calc::WindChill_K_Calc(Inputs* input, int timeStep, double Tair_K, double Tdew_K, double& WindChill_K) {
	//Tair_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tair_C = Tair_K - 273.15;
	//Tdew_C (C) is converted from Tair_K (K) to degrees C by subtracting 273.15
	double Tdew_C = Tdew_K - 273.15;
	//Windspeed_Reference_kmph = 4.8 km/h from Stull (2017), above Eq 3.64a in section 3.7 on Apparent Temperature Indices
	double Windspeed_Reference_kmph = 4.8;
	//Windspeed_Reference_mps = Windspeed_Reference_kmph * Ratio_m_to_km * Ratio_hr_to_s
	double Windspeed_Reference_mps = Windspeed_Reference_kmph * Ratio_m_to_km * Ratio_hr_to_s;

	double WindSpeed_10m_mps = input->WindSpd_mps[timeStep]; 
	//Height_WindSpeed_10_m defined at 10 m according to algorithm and equations from Eq 3.64a and 3.64b of Stull (2017)
	double Height_WindSpeed_10_m = 10;
	//If Height_WindSpeed_Actual_m from HydroPlusConfig.xml not equal to Height_WindSpeed_10_m for Twbg calculations, then enter
	//Note: Height_WindSpeed_Actual_m = 10 m for most NOAA weather stations
	if (Height_WindSpeed_Actual_m != Height_WindSpeed_10_m) {
		//Flag_daytime set to True if Cos_ZenithAngle_Solar > 0, or False otherwise
		Flag_daytime = (Cos_ZenithAngle_Solar > 0) ? 1 : 0;
		//StabilityClass_Atmosphere_int is determined from function StabilityClass_Estimator
		StabilityClass_Atmosphere_int = HeatMetrics_Calc::StabilityClass_Estimator(Flag_daytime, WindSpeed_10m_mps, Radiation_Shortwave_Wpm2, Delta_Height_Tair_m);
		//WindSpeed_wc_mps (m/s) computed in WindSpeed_Estimate_mps_Calc to adjust from measured height to Twbg height
		WindSpeed_10m_mps = HeatMetrics_Calc::WindSpeed_Estimate_mps_Calc(WindSpeed_10m_mps, Height_WindSpeed_Actual_m, Height_WindSpeed_10_m, StabilityClass_Atmosphere_int, Flag_urban);
	}


	//Coefficients from Eq 3.64 of Stull (2017) 
	double Coeff_a = 0.62;
	double Coeff_b = 0.51;
	double Tair_Coeff_1_C = 13.1;
	double Tair_Coeff_2_C = 14.6;
	//Twind_chill_C initiatialized to Tair_C 
	double WindChill_C = Tair_C;

	//If WindSpeed_10m_mps > Windspeed_Reference_mps then
	//Note: Algorithm and equations from Eq 3.64a and 3.64b of Stull (2017)
	if (WindSpeed_10m_mps > Windspeed_Reference_mps) {
		//WindChill_C = (Coeff_a * Tair_C + Tair_Coeff_1_C) + (Coeff_b * Tair_C - Tair_Coeff_2_C) * pow((WindSpeed_10m_mps / Windspeed_Reference_mps), 0.16)
		WindChill_C = (Coeff_a * Tair_C + Tair_Coeff_1_C) + (Coeff_b * Tair_C - Tair_Coeff_2_C) * pow((WindSpeed_10m_mps / Windspeed_Reference_mps), 0.16);
	}
	//Else WindSpeed_10m_mps <= Windspeed_Reference_mps 
	else {
		//WindChill_C = Tair_C
		WindChill_C = Tair_C;
	}

	//If (WindChill_C > Tair_C) then set to Tair_C
	if (WindChill_C > Tair_C) {
		//WindChill_C = Tair_C
		WindChill_C = Tair_C;
	}

	//WindChill_K = WindChill_C + 273.15; convert from C to K
	WindChill_K = WindChill_C + 273.15;
}