Skip to content

Instantly share code, notes, and snippets.

@Sleitnick

Sleitnick/solar.luau

Last active January 23, 2025 04:58
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save Sleitnick/d616319ac9d3b8e6aa1d1a6042c7f1f9 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save Sleitnick/d616319ac9d3b8e6aa1d1a6042c7f1f9 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. Sleitnick revised this gist Jan 23, 2025. 1 changed file with 3 additions and 3 deletions.
    6 changes: 3 additions & 3 deletions solar.luau
    View file
    Original file line number Diff line number Diff line change
    @@ -225,9 +225,9 @@ local function calculateSolarData(date: DateTimeInfo, lat: number, lon: number):
    SunDeclineDeg = sunDeclinDeg,
    EqOfTimeMinutes = eqOfTimeMinutes,
    HASunriseDeg = haSunriseDeg,
    SolarNoon = DateTime.fromTimestamp(t.Timestamp + solarNoon * 86400), --t.Add(time.Duration(solarNoon * 24.0 * float64(time.Hour))),
    SunriseTime = DateTime.fromTimestamp(t.Timestamp + sunriseTime * 86400), --t.Add(time.Duration(sunriseTime * 24.0 * float64(time.Hour))),
    SunsetTime = DateTime.fromTimestamp(t.Timestamp + sunsetTime * 86400), --t.Add(time.Duration(sunsetTime * 24.0 * float64(time.Hour))),
    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,
  2. Sleitnick created this gist Jan 23, 2025.
    244 changes: 244 additions & 0 deletions solar.luau
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,244 @@
    --!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), --t.Add(time.Duration(solarNoon * 24.0 * float64(time.Hour))),
    SunriseTime = DateTime.fromTimestamp(t.Timestamp + sunriseTime * 86400), --t.Add(time.Duration(sunriseTime * 24.0 * float64(time.Hour))),
    SunsetTime = DateTime.fromTimestamp(t.Timestamp + sunsetTime * 86400), --t.Add(time.Duration(sunsetTime * 24.0 * float64(time.Hour))),
    SunlightDuration = sunlightDurationMinutes,
    TrueSolarTime = trueSolarTimeMinutes,
    HourAngleDeg = hourAngleDeg,
    SolarZenithAngleDeg = solarZenithAngleDeg,
    SolarElevationAngleDeg = solarElevationAngleDeg,
    ApproxAtmosphericRefractionDeg = approxAtmosphericRefractionDeg,
    SolarElevationCorrectedForAtmosphericRefractionDeg = solarElevationCorrectedForAtmRefractionDeg,
    SolarAzimuthAngleDegreeClockwiseFromNorth = solarAzimuthAngleDegCwFromNorth,
    })
    end

    return {
    calculateSolarData = calculateSolarData,
    }