#include "SolarCalculation.h"
#include "../Inputs/Inputs.h"

//Note: Based on code used to develop CoolRiver, these functions should be using the date of simulation, yet use Start and End dates
SolarCalculation::SolarCalculation(Inputs* input) :input(input)
{
	int Date_YYYYMMDD;
	//startYear (YYYYMMDD) is set equal to the HydroPlusConfig.xml input StartDate_YYYYMMDD, which is an integer
	Date_YYYYMMDD = input->SimulationNumericalParams["StartDate_YYYYMMDD"];

	//DateParts_Gregorian[4] integer array initialized to hold parts of date for ExtractGregorianDateParts function
	int DateParts_Gregorian[4];
	//ExtractGregorianDateParts function extracts Date_YYYYMMDD year, month, day, hour parts into DateParts_Gregorian vector
	Inputs::ExtractGregorianDateParts(Date_YYYYMMDD, DateParts_Gregorian);
	//startYear is year extracted from DateParts_Gregorian[0]
	startYear = DateParts_Gregorian[0];
	//startMonth is year extracted from DateParts_Gregorian[1]
	startMonth = DateParts_Gregorian[1];
	//startDay is year extracted from DateParts_Gregorian[2]
	startDay = DateParts_Gregorian[2];

	//Initialize variables, defined in SolarCalculation.h and used below
	latDeg = 0;
	latMin = 0;
	lngDeg = 0;
	lngMin = 0;
	Day_JulianDay = 0;

	//stdMeridian defined from input HydroPlusConfig.xml CR_standardMeridian
	//Note: Refactor to replace CR_standardMeridian with GMT offset, and then compute CR_standardMeridian
	stdMeridian = input->CoolRiverVariables["CR_standardMeridian"];
	//SolarConst function called
	SolarConst();
	
	//Date_JulianDay from Date_YYYYMMDD_to_JulianDay function sending Date_YYYYMMDD as YYYYMMDD in first 8 digits
	Day_JulianDay = Inputs::Date_YYYYMMDD_to_JulianDay(Date_YYYYMMDD);

	//SolarDay function called
	SolarDay();
}

//Note: Routine from DHSVM 3.1.2 CalcSolar.c by Wigmosta, function SolarConst()
//Note: Purpose is obtain solar constants and conversion from degree to radians
void SolarCalculation::SolarConst()
{
	latDeg = input->CoolRiverVariables["CR_tarLat"];

	lngDeg = input->CoolRiverVariables["CR_tarLong"];

	Latitude_rad = (latDeg)*M_PI / 180.;
	Longitude = (lngDeg)*M_PI / 180.;
}


//Note: Routine from DHSVM 3.1.2 CalcSolar.c by Wigmosta, function SolarDay()
//Note: Purpose is calculate daily solar values
void SolarCalculation::SolarDay()
{
	/* coefficient for equation of time */
	float B;
	/* adjustment for equation of time (min) */
	float EqnOfTime;
	/* adjustment for longitude (min) */
	float LongitudeAdjust;
	/* cosine of the half-day length */
	float CosineHalfDayLength;
	float HalfDayLength;

	//LongitudeAdjust = (MINPDEG * DEGPRAD) * (stdMeridian - Longitude)
	//Note: Standard Time: Add 4 min per degree away from StandardMeridian (4 min/degree * 180 degree/pi radian) 
	LongitudeAdjust = (MINPDEG * DEGPRAD) * (stdMeridian - Longitude);

	/* equation of time */
	B = (2.0 * M_PI * (Day_JulianDay - 81)) / 364.;
	EqnOfTime = 9.87 * sin(2 * B) - 7.53 * cos(B) - 1.5 * sin(B);

	/* adjustment factor to convert local time to solar time
	solar time = local time + TimeAdjustment  */
	/* for example from GMT to the west coast of the us (PST) */
	/* is a -8 hour shift, i.e. PST = GMT-8 */
	TimeAdjustment = -(LongitudeAdjust + EqnOfTime) / MINPHOUR;

	/* solar Declination  */
	Declination = 0.4098 * sin(2 * M_PI * (284 + Day_JulianDay) / DAYPYEAR);

	/* half-day length  */
	CosineHalfDayLength = -tan(Latitude_rad) * tan(Declination);
	if (CosineHalfDayLength >= 1.0)
		HalfDayLength = M_PI;
	else
		HalfDayLength = acos(CosineHalfDayLength);

	/* convert HalfDayLength from radians to Hours
	1 radian = (180 deg / M_PI) * (1 hr / 15 degrees rotation) */
	HalfDayLength = HalfDayLength / RADPHOUR;

	/* solar time of Sunrise and Sunset  */
	Sunrise = NoonHour - HalfDayLength;
	Sunset = NoonHour + HalfDayLength;

	/* calculate the sun-earth distance */
	SunEarthDist = 1.0 + 0.033 * cos(RADPDEG * (360. * Day_JulianDay / 365));
}

//Note: Routine from DHSVM 3.1.2 Calendar.c by Nijssen, function IsLeapYear()
//Note: Purpose is generic functions to manipulate times and dates
char SolarCalculation::IsLeapYear(int Year)
{
	Year = startYear;
	if ((Year % 4 == 0 && Year % 100 != 0) || Year % 400 == 0)
		return true;
	return false;
}

//Note: Routine header from DHSVM 3.1.2 CalcSolar.c by Wigmosta, function SolarHour()
/*****************************************************************************
Function name: SolarHour()

Purpose: This subroutine calculates position of the sun as a function of
the time of day, the length of time the sun is above the horizon,
and the maximum radiation.

Required:
float Latitude_rad		- site laditude (rad)
float LocalHour		- local time (hr)
float Dt			- length of current timestep (hr)
float NoonHour		- true solar noon (hr)
float Declination		- solar Declination (rad)
float Sunrise		- time of Sunrise (hr)
float Sunset		- time of Sunset (hr)
float TimeAdjustment	- required adjustment to convert local time
to solar time (hr)
float SunEarthDist          - distance from Sun to Earth

Returns: void

Modifies:
float SineSolarAltitude - sine of sun's SolarAltitude
int*Daylight	     - false: measured solar radiation and the sun is
below the horizon.
true: sun is above the horizon
float SolarTimeStep     - fraction of the timestep the sun is above the
horizon
float SunMax            - calculated solar radiation at the top of the
atmosphere (W/m^2)

Comments     : EXECUTE AT START OF EACH TIMESTEP
*****************************************************************************/
//Note: Routine from DHSVM 3.1.2 CalcSolar.c by Wigmosta, function SolarHour()
//Note: Purpose is calculation of position of sun during day
void SolarCalculation::SolarHour(Inputs* input, int timeStep)
{

	//float SolarZenith;		/* sun zenith angle (rads) */
	float StartHour = 0; /* current Hour in solar time (hr) */
	float EndHour = 0;   /* mid-point of current solar Hour (hr) */
	float Hour;			 /* angle of current hour from solar noon
						(rads) */

						/* NOTE THAT HERE Dt IS IN HOURS NOT IN SECONDS */

						/* all calculations based on hour, the solar corrected local time */
	float Dt = 1; // The time step considered 1 hour

	float LocalHour = stof(input->SimulationTime_HMS[timeStep]);
	Hour = LocalHour + TimeAdjustment;

	//Hour = LocalHour + TimeAdjustment - stdMeridian;            //*******************************just test********

	if (Hour < 0) {
		Hour += 24;
	}
	if (Hour > 24) {
		Hour -= 24;
	}

	Daylight = false;
	if ((Hour > Sunrise) && ((Hour - Dt) < Sunset)) {
		Daylight = true;
	}
	else {
		SolarAltitude = 0;
		SolarAzimuth = 0;
	}
	if (Daylight == true) {
		/* sun is above the horizon */

		if (Dt > 0.0) {
			/* compute average solar SolarAltitude over the timestep */
			StartHour = (Hour - Dt > Sunrise) ? Hour - Dt : Sunrise;
			EndHour = (Hour < Sunset) ? Hour : Sunset;

			/*  convert to radians  */
			StartHour = RADPHOUR * (StartHour - NoonHour);
			EndHour = RADPHOUR * (EndHour - NoonHour);
			SolarTimeStep = EndHour - StartHour;

			/*  determine the average geometry of the sun angle  */
			SineSolarAltitude = sin(Latitude_rad) * sin(Declination) + (cos(Latitude_rad) * cos(Declination) * (sin(EndHour) - sin(StartHour)) / SolarTimeStep);
		}
		else {
			Hour = RADPHOUR * (Hour - NoonHour);
			SineSolarAltitude = sin(Latitude_rad) * sin(Declination) + (cos(Latitude_rad) * cos(Declination) * cos(Hour));
		}

		SolarAltitude = asin(SineSolarAltitude);

		SolarZenith = M_PI / 2 - SolarAltitude;

		SolarAzimuth = ((sin(Latitude_rad) * (SineSolarAltitude)-sin(Declination)) / (cos(Latitude_rad) * sin(SolarZenith)));

		if (SolarAzimuth > 1.) {
			SolarAzimuth = 1.;
		}

		SolarAzimuth = acos(-(SolarAzimuth));

		if (Dt > 0.0) {
			if (fabs(EndHour) > fabs(StartHour)) {
				SolarAzimuth = 2 * M_PI - (SolarAzimuth);
			}
		}
		else {
			if (Hour > 0) {
				SolarAzimuth = 2 * M_PI - (SolarAzimuth);
			}
		}
		/*     printf("Solar Azimuth = %g\n", *SolarAzimuth*180./M_PI); */

		SunMax = SOLARCON * SunEarthDist * SineSolarAltitude;
	}
	//Note: End of DHSVM 3.1.2 CalcSolar.c by Wigmosta, function SolarHour()
	//**************************************************************************************

	//Note: Routine from HFlux (Lautz & Glose) hflux_shortwave_refl function and HeatSource (Boyd & Kasper) Sub SubSolarVariables() function
	//Note: Purpose is calculation of position of sun during day
	//Note: Reza used simplified methodology to calculate the dst; need to check the HFlux code's hflux_shortwave_refl routine.

	int Date_YYYYMMDD = input->SimulationDate_YYYYMMDD[timeStep];
	int Date_Year = 0;
	int Date_Month = 0;
	int Date_Day = 0;
	int DateParts_Gregorian[4];

	//ExtractGregorianDateParts function extracts Date_YYYYMMDD year, month, and day date parts from Gregorian date, format YYYYMMDDHH
	Inputs::ExtractGregorianDateParts(Date_YYYYMMDD, DateParts_Gregorian);
	//Date_Year is year extracted from DateParts_Gregorian[0]
	Date_Year = DateParts_Gregorian[0];
	//Date_Month is year extracted from DateParts_Gregorian[1]
	Date_Month = DateParts_Gregorian[1];
	//Date_Day is year extracted from DateParts_Gregorian[2]
	Date_Day = DateParts_Gregorian[2];


	double time_fract = LocalHour * (1.0 / 24.0);

	//t_gmt = time_fract + (-stdMeridian / 24) had been t_gmt = dst + (-stdMeridian / 24), but dst equaled time_fract for Standard Time
	//Note: HFlux had Daylight Saving correction flag and adjustment, which was removed
	//Note: HydroPlus does not use Daylight Saving time in calculations, given Weather.csv is Standard Time
	double t_gmt = time_fract + (-stdMeridian / 24);
	double A1 = round(startYear / 100.0);
	double B1 = 2 - A1 + round(A1 / 4.0);
	double t_jd = round(365.25 * (Date_Year + 4716)) + round(30.6001 * (Date_Month + 1)) + Date_Day + B1 - 1524.5;
	double t_jdc = (t_jd + t_gmt - 2451545) / 36525;

	// Solar position relative to Earth
	double S = 21.448 - t_jdc * (46.815 + t_jdc * (0.00059 - (t_jdc * 0.001813)));
	double O_ob_mean = 23 + ((26 + (S / 60)) / 60);
	double O_ob = O_ob_mean + 0.00256 * cos(125.04 - 1934.136 * t_jdc);
	double e_c = 0.016708634 - t_jdc * (0.000042037 + 0.0000001267 * t_jdc);
	double o_as_mean = 357.52911 + t_jdc * (35999.05029 - 0.0001537 * t_jdc);
	double tempVariable = floor((280.46646 + t_jdc * (36000.76983 + 0.0003032 * t_jdc)) / 360.0);
	double o_ls_mean = (280.46646 + t_jdc * (36000.76983 + 0.0003032 * t_jdc)) - (tempVariable * 360.0);

	double a = (o_as_mean * M_PI / 180);
	double b = sin(a);
	double c = sin(b * 2);
	double d = sin(c * 3);
	double O_cs = b * (1.914602 - t_jdc * (0.004817 + 0.000014 * t_jdc)) + c * (0.019993 - 0.000101 * t_jdc) + d * 0.000289;
	double O_ls = o_ls_mean + O_cs;
	double O_al = O_ls - 0.00569 - (0.00478 * sin((125.04 - 1934.136 * t_jdc) * M_PI / 180));
	double solar_dec = asin(sin(O_ob * M_PI / 180) * sin(O_al * M_PI / 180)) * 180 / M_PI;

	double A = pow((tan(O_ob * M_PI / 360)), 2.0);
	double B = sin(2 * o_ls_mean * M_PI / 180);
	double C = sin(o_as_mean * M_PI / 180);
	double D = cos(2 * o_ls_mean * M_PI / 180);
	double E = sin(4 * o_ls_mean * M_PI / 180);
	double F = sin(2 * o_as_mean * M_PI / 180);

	double e_t = 4 * (A * B - 2 * e_c * C + 4 * e_c * A * C * D - 0.5 * pow(A, 2.0) * E - (4.0 / 3.0) * pow(e_c, 2.0) * F) * (180 / M_PI);
	double h = LocalHour;
	double roundLong = -15 * round(lngDeg / 15);
	double t_s = (h * 60) + 4 * (roundLong + lngDeg) + e_t;
	double o_ha = (t_s / 4 - 180) * M_PI / 180;
	double m = sin(Latitude_rad) * sin(solar_dec * (M_PI / 180)) + (cos(Latitude_rad) * cos(solar_dec * (M_PI / 180)) * cos(o_ha));
	SolarZenithHFlux = acos(m);

	SolAzimuth.push_back(SolarAzimuth);
	SolAltitude.push_back(SolarAltitude);
	SolZenith.push_back(SolarZenithHFlux);

}
