unit Calculate;

interface

  uses Globals;

  type

    SunTime = record Dawn,Dusk:Float; end;

    Daytime = record
      Astronomical,
      Nautical,
      Civil,
      Actual:SunTime;
      end;

  procedure Phenomina(lat,long,zone:float; y:longint;m,d:integer;
                      var Day:DayTime);


implementation

  const
        Sunup    = pi/2;
        Sundown  = pi*3/2;


  function Normalize(z:float):float;
    begin
      z := 2*pi * Frac(z/(2*pi));
      if z<0 then
        z := z + 2*pi;
      Normalize := z;
    end;

  procedure CalcSolar(SunAt,lat,long:float; y:longint; m,d:integer;
                  var ApproxTime,
                      SolarRightAscension, SolarDeclination:float);
    var
      SolarMeanAnomaly,
      SolarTrueLong
        :Float;
      n :integer;
    begin
      ApproxTime := DayOfYear(y,m,d) + (SunAt+long)/(2*pi);
      SolarMeanAnomaly := ApproxTime * 0.017202 - 0.0574039;
      SolarTrueLong := Normalize(
                        SolarMeanAnomaly
                        + 0.0334405 * sin(SolarMeanAnomaly)
                        + 3.49066E-4 * sin(2*SolarMeanAnomaly)
                        + 4.93289
                       );
      {Quadrant Determination}
        if Frac(SolarTrueLong*2/pi) = 0.0 then
          SolarTrueLong := SolarTrueLong + 4.84814e-6;
        n := 2;
        if SolarTrueLong <= pi/2 then
          n := 0
        else if SolarTrueLong <= pi*3/2 then
          n := 1;
      SolarRightAscension := ArcTan(0.91746 * Tan(SolarTrueLong));
      {Quadrant Adjustment}
        if n = 1 then
          SolarRightAscension := SolarRightAscension + pi
        else if n = 2 then
          SolarRightAscension := SolarRightAscension + 2*pi;
      {Solar Declination}
        SolarDeclination := 0.39782 * Sin(SolarTrueLong);
        SolarDeclination := ArcTan(SolarDeclination
                                  / Sqrt(1-Sqr(SolarDeclination))
                                  );
    end;

  procedure ConvertCoords(r,q,lat:float; var s:float);
    begin
      s := (r-(sin(q)*sin(lat))) / (cos(q)*cos(lat));
    end;

  function SettingAdjustment(var SolarDeclination:float):float;
    begin
      SettingAdjustment := pi/2 - ArcTan( SolarDeclination
                                        / Sqrt(1-sqr(SolarDeclination))
                                        );
    end;

  function RisingAdjustment(var SolarDeclination:float):float;
    begin
      RisingAdjustment := 2*pi - SettingAdjustment(SolarDeclination);
    end;

  procedure GetTime(long,zone:float; y,m,d: integer;
                    ApproxTime,SolarRightAscension, SolarDeclination:float;
                var t:float);
    begin
      t := {Loval Apparent Time}
             SolarDeclination + SolarRightAscension
                 - ApproxTime * 2*pi/365.2425 - 1.73364
           {Universal Time}
               + long
           {Wall Clock Time}
               - zone;
      t := Normalize(t) * 24{hours}/(2*pi{radians});
    end;


  procedure Phenomina(lat,long,zone:float; y:longint;m,d:integer;
                      var Day:DayTime);
    var time,ascension,declination:float;
    procedure RisingTime(r:float; var t:SunTime);
      var s:float;
          i:integer;
      begin;
        ConvertCoords(r,declination,lat,s);
        if abs(s)>1.0 then
          t.Dawn := 0
        else
          GetTime(long,zone,y,m,d,time,ascension,RisingAdjustment(s),t.Dawn);
      end;
    procedure SettingTime(r:float; var t:SunTime);
      var s:float;
          i:integer;
      begin;
        ConvertCoords(r,declination,lat,s);
        if abs(s)>1.0 then
          t.Dusk := 0
        else
          GetTime(long,zone,y,m,d,time,ascension,SettingAdjustment(s),t.Dusk);
      end;
    const
       Astro   = -0.309017;
       Nautic  = -0.207912;
       Civil   = -0.104528;
       Actual  = -0.0145439;
    begin
      lat  := Radians(lat);
      long := Radians(long);
      zone := Radians(zone)*15.0; {15 degrees/Time Zone}

      CalcSolar(Sunup,lat,long, y,m,d, time,ascension,declination);
      RisingTime(Astro  ,Day.Astronomical);
      RisingTime(Nautic ,Day.Nautical);
      RisingTime(Civil  ,Day.Civil);
      RisingTime(Actual ,Day.Actual);

      CalcSolar(Sundown,lat,long, y,m,d, time,ascension,declination);
      SettingTime(Actual ,Day.Actual);
      SettingTime(Civil  ,Day.Civil);
      SettingTime(Nautic ,Day.Nautical);
      SettingTime(Astro  ,Day.Astronomical);
    end;

end.