-module(sonnenstand).

-moduledoc """
Solar position calculations.

Given a date/time and a location (latitude/longitude in degrees), this module
computes the position of the sun as **elevation** (height above the horizon)
and **azimuth** (compass direction).

The formulas follow the algorithm described at
<http://www.geoastro.de/SME/tk/index.htm>. They are approximations intended
for everyday use (shading, solar-panel orientation, photography), not for
high-precision astronomy.

## Conventions

- Angles are in **degrees**.
- Latitude is positive **north**, longitude positive **east**.
- Times are interpreted as **UTC**.
- Azimuth is measured **clockwise from north**: `0.0` = north, `90.0` = east,
  `180.0` = south, `270.0` = west.

## Example

```erlang
1> DateTime = {{2024, 6, 21}, {12, 0, 0}}.
2> sonnenstand:elevation(DateTime, 52.5, 13.4).
60.91...
3> sonnenstand:azimuth(DateTime, 52.5, 13.4).
176.13...
```
""".

-export([day_of_year/1,
         declination/1,
         equation_of_time/1,
         hour_angle/2,
         elevation/3,
         azimuth/3]).

-type latitude() :: number().
%% Geographic latitude in degrees, positive north (range -90.0..90.0).

-type longitude() :: number().
%% Geographic longitude in degrees, positive east (range -180.0..180.0).

-type degrees() :: float().
%% An angle expressed in degrees.

-export_type([latitude/0, longitude/0, degrees/0]).

%%====================================================================
%% Public API
%%====================================================================

-doc """
Returns the ordinal day of the year for `Date` (1 for 1 January).

```erlang
1> sonnenstand:day_of_year({2024, 1, 1}).
1
2> sonnenstand:day_of_year({2024, 12, 31}).
366
```
""".
-spec day_of_year(calendar:date()) -> pos_integer().
day_of_year({Year, _Month, _Day} = Date) ->
    FirstOfYear = {Year, 1, 1},
    Midnight = {0, 0, 0},
    {Days, _} = calendar:time_difference({FirstOfYear, Midnight},
                                         {Date, Midnight}),
    Days + 1.

-doc """
Returns the sun's declination in degrees for `Date`.

The declination is the angle between the sun's rays and the equatorial plane,
ranging from about `-23.45` at the December solstice to `+23.45` at the June
solstice.
""".
-spec declination(calendar:date()) -> degrees().
declination(Date) ->
    -23.45 * math:cos(deg2rad(360 * (day_of_year(Date) + 10) / 365)).

-doc """
Returns the equation of time in **minutes** for `Date`.

This is the difference between apparent solar time and mean solar time, caused
by the eccentricity of Earth's orbit and the tilt of its axis.
""".
-spec equation_of_time(calendar:date()) -> float().
equation_of_time(Date) ->
    N = day_of_year(Date),
    60 * (-0.171 * math:sin(0.0337 * N + 0.465)
          - 0.1299 * math:sin(0.01787 * N - 0.168)).

-doc """
Returns the sun's hour angle in degrees for `DateTime` (UTC) at `Longitude`.

The hour angle is negative before solar noon (morning), zero at solar noon, and
positive afterwards (afternoon).
""".
-spec hour_angle(calendar:datetime(), longitude()) -> degrees().
hour_angle({Date, {Hour, Minute, Second}}, Longitude) ->
    Clock = Hour + Minute / 60 + Second / 3600,
    15 * (Clock + Longitude / 15 - 12 + equation_of_time(Date) / 60).

-doc """
Returns the sun's elevation in degrees for `DateTime` at the given location.

The elevation (or altitude) is the angle of the sun above the horizon:
positive when the sun is up, negative when it is below the horizon.
""".
-spec elevation(calendar:datetime(), latitude(), longitude()) -> degrees().
elevation(DateTime, Latitude, Longitude) ->
    rad2deg(math:asin(sin_elevation(DateTime, Latitude, Longitude))).

-doc """
Returns the sun's azimuth in degrees for `DateTime` at the given location.

The azimuth is measured clockwise from north (`0.0` = north, `90.0` = east,
`180.0` = south, `270.0` = west). Values before solar noon fall in the eastern
half (`0.0`..`180.0`); values afterwards fall in the western half
(`180.0`..`360.0`).
""".
-spec azimuth(calendar:datetime(), latitude(), longitude()) -> degrees().
azimuth(DateTime, Latitude, Longitude) ->
    Base = rad2deg(math:acos(cos_azimuth(DateTime, Latitude, Longitude))),
    %% acos only yields 0..180 degrees, which covers the eastern (morning) sky.
    %% In the afternoon (positive hour angle) the sun is in the west, so the
    %% azimuth is mirrored across the north-south meridian.
    case hour_angle(DateTime, Longitude) > 0 of
        true -> 360.0 - Base;
        false -> Base
    end.

%%====================================================================
%% Internal functions
%%====================================================================

%% sin(elevation): the third coordinate of the sun's unit vector.
-spec sin_elevation(calendar:datetime(), latitude(), longitude()) -> float().
sin_elevation({Date, _Time} = DateTime, Latitude, Longitude) ->
    Decl = declination(Date),
    math:sin(deg2rad(Latitude)) * math:sin(deg2rad(Decl))
        + math:cos(deg2rad(Latitude)) * math:cos(deg2rad(Decl))
              * math:cos(deg2rad(hour_angle(DateTime, Longitude))).

%% cos(azimuth), clamped to [-1.0, 1.0] to keep acos well-defined against
%% floating-point rounding near the horizon.
-spec cos_azimuth(calendar:datetime(), latitude(), longitude()) -> float().
cos_azimuth({Date, _Time} = DateTime, Latitude, Longitude) ->
    Decl = declination(Date),
    Elev = elevation(DateTime, Latitude, Longitude),
    Y = -(math:sin(deg2rad(Latitude)) * math:sin(deg2rad(Elev))
          - math:sin(deg2rad(Decl)))
        / (math:cos(deg2rad(Latitude))
           * math:sin(math:acos(math:sin(deg2rad(Elev))))),
    max(-1.0, min(1.0, Y)).

-spec deg2rad(number()) -> float().
deg2rad(Degrees) ->
    Degrees * math:pi() / 180.

-spec rad2deg(number()) -> float().
rad2deg(Radians) ->
    Radians * 180 / math:pi().