procedure Sun_RiseSet2 (const DT: TDateTime; const Position: TESBPosition; const TimeOfs: Extended;
    out SunRise, SunSet: TDateTime);
var
    Gamma, Declination, EqTime, HourAngle: Extended;
begin
    Gamma := TwoPi * (DayOfYear (DT) - 1.5);
    EqTime := 229.18 * (0.000075 + 0.001868 * Cos (Gamma) - 0.032077 * Sin (Gamma)
         - 0.014615 * cos (2 * Gamma) - 0.040849 * sin (2 * Gamma)); // In Minutes

    Declination := 0.006918 - 0.399912 * cos (Gamma) + 0.070257 * sin (Gamma)
         - 0.006758 * cos (2 * Gamma) + 0.000907 * sin (2 * Gamma)
         - 0.002697 * cos (3 * Gamma) + 0.00148 * sin (3 * Gamma); // In Radians

    HourAngle := ESBArcCosDeg (ESBCosDeg (90.833) / (EsbCosDeg (Position.Latitude)
         * Cos (Declination)) - ESBTanDeg (Position.Latitude) * ESBTan (Declination)); // In Degrees

    SunRise := (720 + 4 * (-1 * Position.Longitude - HourAngle) - EqTime) / 60 + TimeOfs;
    SunSet := (720 + 4 * (-1 * Position.Longitude + HourAngle) - EqTime) / 60 + TimeOfs;

    if SunRise > 24 then
         SunRise := SunRise - 24;
    if SunRise < 0 then
         SunRise := SunRise + 24;
    if SunSet > 24 then
         SunSet := SunSet - 24;
    if SunSet < 0 then
         SunSet := SunSet + 24;

    SunRise := Int (DT) + SunRise / 24;
    SunSet := Int (DT) + SunSet / 24;
end;