#ifndef HeatMetrics_Calc_H
#define HeatMetrics_Calc_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 <vector>
#include <string>
#include <iostream>
//_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 "../Inputs/Inputs.h"
#include "../DataFolder.h"
#include "../Inputs/SolarCalculator.h"

//Ratio_Radian_to_Degree is 0.017453292519943295
#define Ratio_Radian_to_Degree (M_PI / 180.0) 
//Radiation_Shortwave_SolarConstant_Wpm2 is 1367 W/m2
#define Radiation_Shortwave_SolarConstant_Wpm2 1367.0
//Stefan-Boltzmann constant is 5.670374419E-8 W/(m^2*K^4)
#define Sigma_Stefan_Boltzmann 5.670374419E-8
//Ratio_MolecularWeightVapor_to_MolecularWeightDryAir is ratio of molecular weight for water vapor to dry air (g/g)
#define Ratio_MolecularWeightVapor_to_MolecularWeightDryAir 0.622
//Specific heat capacity of air at constant pressure (J/kg/K) 1006.0 at ~290 to 315 K, 
//Note: 1003.5 used in Liljegren et al. (2008), but this does not seem to make physical sense
#define SpecificHeat_DryAir_J_p_kg_K 1006.0
//Molecular_Weight_DryAir_kg_p_kmol is 28.97 kg/kilomole or g/mole of dry air
#define Molecular_Weight_DryAir_kg_p_kmol 28.97 
//Molecular_Weight_WaterVapor_kg_p_kmol is 18.015 kg/kilomole or g/mole of water vapor
#define Molecular_Weight_WaterVapor_kg_p_kmol 18.015 
//R_Universal_Gas_Constant_J_p_kmol_K is universal gas constant 8314.34 J/(kilomole K) 
#define R_Universal_Gas_Constant_J_p_kmol_K 8314.34 
//Cos_ZenithAngle_Minimum minimum cos zenith angle, cos(89.5 deg) = cos(1.562 rad) Eq 13 and Eq 14 Liljegren et al. (2008)
#define Cos_ZenithAngle_Minimum 0.00873 
//Ratio_RadiationShortwaveActual_to_RadiationShortwaveMax_Minimum_frac is 0.85 from Liljegren et al. (2008) github code
#define Ratio_RadiationShortwaveActual_to_RadiationShortwaveMax_Minimum_frac 0.85
//Height_WindSpeed_Reference_m (m) set to 2.0 from Liljegren et al. (2008) github code; height experienced for Twbg
#define Height_WindSpeed_Reference_m 2.0
//WindSpeed_Minimum_mps (m/s) is 0.13 set by Liljegren et al. (2008) github code
#define WindSpeed_Minimum_mps 0.13
//Ratio_m_to_km is 1000.
#define Ratio_m_to_km 1000.0
//Ratio_hr_to_s is 1/3600.
#define Ratio_hr_to_s 1/3600.0
//Tolerance_Convergence is 0.02 set by Liljegren et al. (2008) github code
#define Tolerance_Convergence 0.02
//Iterations_Maximum is 500 set by Liljegren et al. (2008) github code
#define Iterations_Maximum 500

//MAX function for finding maximum 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 MAX(x_var,y_var) (((x_var)>=(y_var)) ? (x_var) : (y_var))        
#define MIN(x_var,y_var) (((x_var)<=(y_var)) ? (x_var) : (y_var))        

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

//class HeatMetrics_Calc
class HeatMetrics_Calc
{
private:

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

	//vector<vector<pair<double, double>>> Latitude_Longitude_vectorPair_dd vector created to store latitude and longitude as decimal degree
	vector<vector<pair<double, double>>> Latitude_Longitude_vectorPair_dd;
	double Cp_Ratio_DryAir_to_WaterVapor = 0;
	double Ratio_GasConstant_to_MolecularWeightAir_m2_p_s2_K = 0;
	double Prandtl_Number_wind = 0;
	double Wick_Emissivity_frac = 0;
	double Wick_Albedo_frac = 0;
	double Wick_Diameter_m = 0;
	double Wick_Length_m = 0;
	double Globe_Emissivity_frac = 0;
	double Globe_Albedo_frac = 0;
	double Globe_Diameter_m = 0;
	double SurfaceGround_Emissivity_frac = 0;
	double SurfaceGround_Albedo_frac = 0;

	int Year_YYYY = 0;
	int Month_MM = 0;
	int Day_DD = 0;
	int Hour_HH = 0;
	int Minute_Mn = 0;
	int GMT_Offset_int = 0;
	int Time_averaging_MeteorologicalInputs_minutes = 0;
	double Latitude_NH_pos_dd = 0;
	double Longitude_WH_pos_dd = 0;
	double AtmosphericPressure_hPa = 0;
	double Tair_C = 0;
	double Tdew_C = 0;
	double RelativeHumidity_percent = 0;
	double WindSpeed_mps = 0;
	double Height_WindSpeed_Actual_m = 0;
	double Delta_Height_Tair_m = 0;
	int Flag_urban = 0;
	int Flag_daytime = 0;

	double Radiation_Shortwave_Wpm2 = 0;
	double Radiation_Shortwave_TopOfAtm_Wpm2 = 0;
	double Radiation_Shortwave_adjusted_Wpm2 = 0;
	double Ratio_RadiationShortwaveActual_to_RadiationShortwaveMax_frac = 0;
	double Radiation_Shortwave_direct_frac = 0;

	int Hour_GreenwichMeanTime_int = 0;
	double Day_and_Hour_decimal = 0;
	double Cos_ZenithAngle_Solar = 0;
	double days_1900 = 0;
	double Solar_RightAscension_Apparent_hr = 0;
	double Solar_Declination_Apparent_rad = 0; 
	double Solar_AltitudeAngle_rad = 0; 
	double Solar_AltitudeAngle_deg = 0;
	double Solar_AzimuthAngle_rad = 0;
	double Distance_Earth_to_Sun_astrologicalunits_frac = 0;
	double Solar_AltitudeAngle_RefractionCorrection_rad = 0;
	double Solar_AltitudeAngle_RefractionCorrection_deg = 0;
	int StabilityClass_Atmosphere_int = 0;
	double WindSpeed_estimated_mps = 0;

	double Solar_HourAngle_Local_rad = 0;
	double Tangent_Solar_AltitudeAngle = 0;

	double Cos_Solar_AltitudeAngle = 0;
	double Cos_Solar_AzimuthAngle = 0;
	double Sine_Solar_AzimuthAngle = 0;

	double Tglobe_new_K = 0;
	double Twetbulb_new_K = 0;
	double Twetbulb_new_C = 0;

	double Twetbulbglobe_C = 0;
	double Tglobe_C = 0;
	double Twetbulbnatural_C = 0; 
	double Twetbulbpsychrometric_C = 0;
	
	double ThermalMeanRadiantTemp_C = 0;
	double UniversalThermalClimateIndex_C = 0;

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

	//WetBulbGlobeTemperature_C_Calc function will manage computation of the wet bulb globe temperature and its return
	void 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);
	double 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, int rad_flag);
	double 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, int MapPixel_ID, int timeStep);
	double HeatTransfer_Air_Sphere_Wpm2K_Calc(double Sphere_Diameter_m, double Tair_K, double AtmosphericPressure_hPa, double WindSpeed_mps);
	double HeatTransfer_Air_Cylinder_Wpm2K_Calc(double Diameter_cylinder_m, double Length_cylinder_m, double Tair_K, double AtmosphericPressure_hPa, double WindSpeed_mps);
	int StabilityClass_Estimator(bool Flag_daytime, double WindSpeed_mps, double Radiation_Shortwave_Wpm2, double Delta_Height_Tair_m);
	double WindSpeed_Estimate_mps_Calc(double WindSpeed_mps, double Height_WindSpeed_Actual_m, double Height_WindSpeed_Target_m, int StabilityClass_Atmosphere_int, bool Flag_urban);
	double Diffusivity_WaterVapor_m2ps_Calc(double Tair_K, double AtmosphericPressure_hPa);
	double Dewpoint_Temperature_K_Calc(double VaporPres_hPa, int phase);
	double ThermalConductivity_DryAir_WpmK_Calc(double Tair_K);
	double SpecificHeat_HumidAir_JpkgK_Calc(double Tair_K, double Tdew_K, double AtmPres_kpa);
	double VaporPressure_haPa_Calc(double Tair_K, int Flag_WaterPhase_is_Ice);
	double Viscosity_Air_Pa_s_Calc(double Tair_K);
	double Emissivity_Atmosphere_Calc(double Tair_K, double RelativeHumidity_frac);
	double LatentHeat_Vaporization_JpkgK_Calc(double Tair_K);
	double Density_Water_kgpm3_Calc(double Tair_K);

	//HeatIndex_K_Calc function will manage computation of the heat index temperature and its return
	double HeatIndex_K_Calc(Inputs* input, double Tair_K, double Tdew_K);
	//Humidex_K_Calc function will manage computation of the humidex temperature and its return
	double Humidex_K_Calc(Inputs* input, CompactRagged* beC, int MapPixel_ID, int DataFolder_ID);
	double UniversalThermalClimateIndex_C_Calc(Inputs* input, int MapPixel_ID, int timeStep, double Tair_C, double WindSpeed_mps, double VaporPressure_Actual_kPa, double Tglobe_C);
	void WindChill_K_Calc(Inputs* input, int timeStep, double Tair_K, double Tdew_K, double& WindChill_K);
};

#endif