solar.luau

--!native
--!strict

-- Source of calculations: https://gml.noaa.gov/grad/solcalc/calcdetails.html

export type SolarData = {
	JulianDay: number,
	JulianCentury: number,
	GeometricMeanLongitudeSunDeg: number,
	GeometricMeanAnomSunDeg: number,
	EccentricEarthOrbit: number,
	SunEqOfCenter: number,
	SunTrueLongitudeDeg: number,
	SunTrueAnomDeg: number,
	SunRadVectorAus: number,
	SunAppLongDeg: number,
	MeanObliqueEclipticDeg: number,
	ObliqueCorrectionDeg: number,
	SunRtAscentDeg: number,
	SunDeclineDeg: number,
	EqOfTimeMinutes: number,
	HASunriseDeg: number,
	SolarNoon: DateTime,
	SunriseTime: DateTime,
	SunsetTime: DateTime,
	SunlightDuration: number,
	TrueSolarTime: number,
	HourAngleDeg: number,
	SolarZenithAngleDeg: number,
	SolarElevationAngleDeg: number,
	ApproxAtmosphericRefractionDeg: number,
	SolarElevationCorrectedForAtmosphericRefractionDeg: number,
	SolarAzimuthAngleDegreeClockwiseFromNorth: number,
}

export type SunriseSunsetData = {
	SunriseTime: DateTime,
	SunsetTime: DateTime,
}

type DateTimeInfo = {
	Year: number,
	Month: number,
	Day: number,
	UTCOffset: number, -- UTC east offset in seconds (e.g. EST-5 would be -7200)
}

local function calculateSolarData(date: DateTimeInfo, lat: number, lon: number): SolarData
	-- Calculation source: https://gml.noaa.gov/grad/solcalc/calcdetails.html

	local timezoneOffset = date.UTCOffset / 3600
	local t = DateTime.fromLocal(date.Year, date.Month, date.Day)

	-- Letters below indicate the spreadsheet column in the calculation source above

	-- E
	local timePastLocalMidnight = 0.5

	-- F
	local julianDay = t.Timestamp / 86400 + 2440587.5 + 0.5

	-- G
	local julianCentury = (julianDay - 2451545) / 36525

	-- I
	local geomMeanLongSunDeg = (280.46646 + julianCentury * (36000.76983 + julianCentury * 0.0003032)) % 360

	-- J
	local geomMeanAnomSunDeg = 357.52911 + julianCentury * (35999.05029 - 0.0001537 * julianCentury)

	-- K
	local eccentEarthOrbit = 0.016708634 - julianCentury * (0.000042037 + 0.0000001267 * julianCentury)

	-- L
	local sunEqOfCtr = math.sin(math.rad(geomMeanAnomSunDeg))
			* (1.914602 - julianCentury * (0.004817 + 0.000014 * julianCentury))
		+ math.sin(math.rad(2.0 * geomMeanAnomSunDeg)) * (0.019993 - 0.000101 * julianCentury)
		+ math.sin(math.rad(3.0 * geomMeanAnomSunDeg)) * 0.000289

	-- M
	local sunTrueLongDeg = geomMeanLongSunDeg + sunEqOfCtr

	-- N
	local sunTrueAnomDeg = geomMeanAnomSunDeg + sunEqOfCtr

	-- O
	local sunRadVectorAus = (1.000001018 * (1 - eccentEarthOrbit * eccentEarthOrbit))
		/ (1 + eccentEarthOrbit * math.cos(math.rad(sunTrueAnomDeg)))

	-- P
	local sunAppLongDeg = sunTrueLongDeg - 0.00569 - 0.00478 * math.sin(math.rad(125.04 - 1934.136 * julianCentury))

	-- Q
	local meanObliqEclipticDeg = 23
		+ (26 + (21.448 - julianCentury * (46.815 + julianCentury * (0.00059 - julianCentury * 0.001813))) / 60)
			/ 60

	-- R
	local obliqCorrDeg = meanObliqEclipticDeg + 0.00256 * math.cos(math.rad(125.04 - 1934.136 * julianCentury))

	-- S
	local sunRtAscentDeg = math.deg(
		math.atan2(
			math.cos(math.rad(obliqCorrDeg)) * math.sin(math.rad(sunAppLongDeg)),
			math.cos(math.rad(sunAppLongDeg))
		)
	)

	-- T
	local sunDeclinDeg = math.deg(math.asin(math.sin(math.rad(obliqCorrDeg)) * math.sin(math.rad(sunAppLongDeg))))

	-- U
	local varY = math.tan(math.rad(obliqCorrDeg / 2)) * math.tan(math.rad(obliqCorrDeg / 2))

	-- V
	local eqOfTimeMinutes = 4
		* math.deg(
			varY * math.sin(2 * math.rad(geomMeanLongSunDeg))
				- 2 * eccentEarthOrbit * math.sin(math.rad(geomMeanAnomSunDeg))
				+ 4 * eccentEarthOrbit * varY * math.sin(math.rad(geomMeanAnomSunDeg)) * math.cos(
					2 * math.rad(geomMeanLongSunDeg)
				)
				- 0.5 * varY * varY * math.sin(4 * math.rad(geomMeanLongSunDeg))
				- 1.25 * eccentEarthOrbit * eccentEarthOrbit * math.sin(2 * math.rad(geomMeanAnomSunDeg))
		)

	-- W
	local haSunriseDeg = math.deg(
		math.acos(
			math.cos(math.rad(90.833)) / (math.cos(math.rad(lat)) * math.cos(math.rad(sunDeclinDeg)))
				- math.tan(math.rad(lat)) * math.tan(math.rad(sunDeclinDeg))
		)
	)

	-- X
	local solarNoon = (720 - 4 * lon - eqOfTimeMinutes + timezoneOffset * 60) / 1440

	-- Y
	local sunriseTime = (solarNoon * 1440 - haSunriseDeg * 4) / 1440

	-- Z
	local sunsetTime = (solarNoon * 1440 + haSunriseDeg * 4) / 1440

	-- AA
	local sunlightDurationMinutes = 8 * haSunriseDeg

	-- AB
	local trueSolarTimeMinutes = (timePastLocalMidnight * 1440 + eqOfTimeMinutes + 4 * lon - 60 * timezoneOffset) % 1440

	-- AC
	local hourAngleDeg: number
	if trueSolarTimeMinutes / 4 < 0 then
		hourAngleDeg = trueSolarTimeMinutes / 4 + 180
	else
		hourAngleDeg = trueSolarTimeMinutes / 4 - 180
	end

	-- AD
	local solarZenithAngleDeg = math.deg(
		math.acos(
			math.sin(math.rad(lat)) * math.sin(math.rad(sunDeclinDeg))
				+ math.cos(math.rad(lat)) * math.cos(math.rad(sunDeclinDeg)) * math.cos(math.rad(hourAngleDeg))
		)
	)

	-- AE
	local solarElevationAngleDeg = 90 - solarZenithAngleDeg

	-- AF
	local approxAtmosphericRefractionDeg: number
	if solarElevationAngleDeg > 85 then
		approxAtmosphericRefractionDeg = 0
	elseif solarElevationAngleDeg > 5 then
		approxAtmosphericRefractionDeg = 58.1 / math.tan(math.rad(solarElevationAngleDeg))
			- 0.07 / math.pow(math.tan(math.rad(solarElevationAngleDeg)), 3)
			+ 0.000086 / math.pow(math.tan(math.rad(solarElevationAngleDeg)), 5)
	elseif solarElevationAngleDeg > -0.575 then
		approxAtmosphericRefractionDeg = 1735
			+ solarElevationAngleDeg
				* (-518.2 + solarElevationAngleDeg * (103.4 + solarElevationAngleDeg * (-12.79 + solarElevationAngleDeg * 0.711)))
	else
		approxAtmosphericRefractionDeg = -20.772 / math.tan(math.rad(solarElevationAngleDeg))
	end
	approxAtmosphericRefractionDeg /= 3600

	-- AG
	local solarElevationCorrectedForAtmRefractionDeg = solarElevationAngleDeg + approxAtmosphericRefractionDeg

	-- AH
	local solarAzimuthAngleDegCwFromNorth: number
	if hourAngleDeg > 0 then
		solarAzimuthAngleDegCwFromNorth = math.deg(
			math.acos(
				((math.sin(math.rad(lat)) * math.cos(math.rad(solarZenithAngleDeg))) - math.sin(math.rad(sunDeclinDeg)))
					/ (math.cos(math.rad(lat)) * math.sin(math.rad(solarZenithAngleDeg)))
			)
		) + 180
	else
		solarAzimuthAngleDegCwFromNorth = 540
			- math.deg(
				math.acos(
					(
						(math.sin(math.rad(lat)) * math.cos(math.rad(solarZenithAngleDeg)))
						- math.sin(math.rad(sunDeclinDeg))
					) / (math.cos(math.rad(lat)) * math.sin(math.rad(solarZenithAngleDeg)))
				)
			)
	end
	solarAzimuthAngleDegCwFromNorth %= 360

	return table.freeze({
		JulianDay = julianDay,
		JulianCentury = julianCentury,
		GeometricMeanLongitudeSunDeg = geomMeanLongSunDeg,
		GeometricMeanAnomSunDeg = geomMeanAnomSunDeg,
		EccentricEarthOrbit = eccentEarthOrbit,
		SunEqOfCenter = sunEqOfCtr,
		SunTrueLongitudeDeg = sunTrueLongDeg,
		SunTrueAnomDeg = sunTrueAnomDeg,
		SunRadVectorAus = sunRadVectorAus,
		SunAppLongDeg = sunAppLongDeg,
		MeanObliqueEclipticDeg = meanObliqEclipticDeg,
		ObliqueCorrectionDeg = obliqCorrDeg,
		SunRtAscentDeg = sunRtAscentDeg,
		SunDeclineDeg = sunDeclinDeg,
		EqOfTimeMinutes = eqOfTimeMinutes,
		HASunriseDeg = haSunriseDeg,
		SolarNoon = DateTime.fromTimestamp(t.Timestamp + solarNoon * 86400),
		SunriseTime = DateTime.fromTimestamp(t.Timestamp + sunriseTime * 86400),
		SunsetTime = DateTime.fromTimestamp(t.Timestamp + sunsetTime * 86400),
		SunlightDuration = sunlightDurationMinutes,
		TrueSolarTime = trueSolarTimeMinutes,
		HourAngleDeg = hourAngleDeg,
		SolarZenithAngleDeg = solarZenithAngleDeg,
		SolarElevationAngleDeg = solarElevationAngleDeg,
		ApproxAtmosphericRefractionDeg = approxAtmosphericRefractionDeg,
		SolarElevationCorrectedForAtmosphericRefractionDeg = solarElevationCorrectedForAtmRefractionDeg,
		SolarAzimuthAngleDegreeClockwiseFromNorth = solarAzimuthAngleDegCwFromNorth,
	})
end

return {
	calculateSolarData = calculateSolarData,
}