[gelöst] Kennt jemand Mond/Sonnen-Komponenten?

[catweasel]

Moin

Ich bin auf der Suche nach einer (freien) Komponente um in meinem Kalenderprogramm noch Sonnen auf-/ Untergangszeiten sowie Mondphasen eintragen zu können.
Kennt jemand eine entsprechende download Seite?

Gruß
Michael

[six1]

Hi,
habe "irgendwo" die AstroUnit heruntergeladen.
Dazu gab es keine Copyright Hinweise.
EDIT: kam von hier: https://www.delphipraxis.net/145690-astro-daten.html

Leider ist die nur für Windows, ich versuche gerade die für Linux fit zu machen.

Code: Alles auswählen

unit AstroUnit;

interface

uses
  {Windows, }SysUtils, Classes, Controls, StdCtrls, DateUtils
  {$ifdef Unix}
//   , LibC
   , BaseUnix
   , unixutil
  {$else}
   , Windows
  {$endif}
  ;


type
  E_NoRiseSet=class(Exception);
  E_OutOfAlgorithRange=class(Exception);
  TGeoAngle = Extended;

type
  TEclipse=(ecNone, ecPartial, ecNoncentral, ecCircular, ecCirculartotal, ecTotal, ecHalfshadow);
  TdeState = (de_Baden_Wuerttemberg, de_Bayern, de_Berlin, de_Brandenburg,
              de_Bremen, de_Hamburg, de_Hessen, de_Mecklenburg_Vorpommern,
              de_Niedersachsen, de_Nordrhein_Westfalen, de_Rheinland_Pfalz,
              de_Saarland, de_Sachsen, de_Sachsen_Anhalt, de_Schleswig_Holstein,
              de_Thueringen, de_None);

  TMoonPhase = (mpNewmoon, mpFirstQuarter, mpFullmoon, mpLastQuarter);
  TSeason = (seWinter, seSpring, seSummer, seAutumn, seNone);
  TMoonPhase8 = (mp_NewMoon, mp_WaxingMoonFirstEighth, mp_WaxingMoonFirstQuarter, mp_WaxingMoonThirdEighth,
                               mp_FullMoon, mp_WaningMoonFifthEighth, mp_WaningMoonLastQuarter, mp_WaningMoonSevensEighth);

const
  Days : array[1..7] of string =
   ('Sonntag', 'Montag', 'Dienstag', 'Mittwoch', 'Donnerstag', 'Freitag', 'Sonnabend');
  Bundesland : array [TdeState] of string =
   ('Baden-Württemberg', 'Bayern', 'Berlin', 'Brandenburg', 'Bremen', 'Hamburg',
    'Hessen', 'Mecklenburg-Vorpommern', 'Niedersachsen', 'Nordrhein-Westfalen',
    'Rheinland-Pfalz', 'Saarland', 'Sachsen', 'Sachsen-Anhalt', 'Schleswig-Holstein',
    'Thüringen', 'Benutzerdefiniert');
  EclipseName : array [TEclipse] of string =
   ('', 'Teilfinsternis', 'Teilfinsternis', 'Randfinsternis',
    'Teilfinsternis mit Rand', 'Totale Finsternis', 'Halbschatten');
  GermanStateLong : array [0..15] of extended =
   (-9, -11.5, -13.4, -13.4, -8.8, -10, -8.7, -12.2, -8.8, -7.5, -7.3, -7, -14, -11.7, -10.2, -11);
  GermanStateLat : array [0..15] of extended =
   (48.6, 48.8, 52.5, 52.5, 53.1, 53.5, 50.5, 53.7, 53.1, 51.6, 50.2, 49.2, 51, 52, 54.3, 51);
  GeogrZentrum : array [TdeState] of string =
   ('Tübingen', 'Ingolstadt', 'Berlin', 'Berlin', 'Bremen', 'Hamburg', 'Gießen',
    'Güstrow', 'Bremen', 'Dortmund', 'Cochem', 'Saarbrücken', 'Dresden', 'Magdeburg',
    'Kiel', 'Erfurt', '');
  Jahreszeit : array [TSeason] of string =
    ('Winter','Frühling', 'Sommer', 'Herbst', '');
  Sternzeichen : array [1..12] of string =
   ('Wassermann', 'Fische', 'Widder', 'Stier', 'Zwilling', 'Krebs', 'Löwe', 'Jungfrau',
    'Waage', 'Skorpion', 'Schütze', 'Steinbock');
  Feiertage : array [1..17] of string =
   ('Neujahr', 'Maifeiertag', 'Tag der deutschen Einheit', 'Allerheiligen',
    '1. Weihnachtstag', '2. Weihnachtstag', 'Karfreitag', 'Ostersonntag',
    'Ostermontag', 'Christi Himmelfahrt', 'Pfingstsonntag', 'Pfingstmontag',
    'Fronleichnam', 'Heilige 3 Könige', 'Mariä Himmelfahrt', 'Reformationstag',
    'Buß- und Bettag');
  Sondertage : array [1..26] of string =
   ('Mariä Lichtmeß', 'Valentinstag', 'Weiberfastnacht', 'Rosenmontag', 'Fastnacht',
    'Aschermittwoch', 'Mariä Verkündigung', 'Palmsonntag', 'Gründonnerstag', 'Muttertag',
    'Peter und Paul', 'Mariä Geburt', 'Erntedankfest', 'Mariä Empfängnis', 'Silvester',
    '1. Advent', '2. Advent', '3. Advent', '4. Advent', 'Heiligabend', 'Frühlingsanfang',
    'Sommmeranfang', 'Herbstanfang', 'Winteranfang', 'Totensonntag', 'Volkstrauertag');

 var
    FAutumnDate          : longint;
    FdeState             : TdeState;
    FHoliday             : string;
    FHolidayNr           : integer;
    FMoonDistance        : extended;
    FMoonPhase           : TMoonPhase;
    FSeason              : TSeason;
    FSpringDate          : longint;
    FSummerDate          : longint;
    FWinterDate          : longint;
    eX                   : array[0..2000] of Extended; // rel. Position in X
    eY                   : array[0..2000] of extended; // rel. Position in Y
                                                       // Bezug Mittelpunkt im
                                                       // Kreis mit R=1
    EUp, EDown, ECurrent : Integer;                    // Konst. f. Bahn


function GetDayOfYear(CDate:TDateTime):Integer;
function GetDaysPerYear(CDate:TDateTime):Integer;
function IsLeapYear(CDate:TDateTime):boolean;
function GetWeekOfYear(CDate:TDateTime):Integer;
function GetWeeksPerYear(CDate:TDateTime):Integer;
function GetDaysPerMonth(CDate:TDateTime):Integer;
function GetFirstDayOfWeek(CDate:TDateTime; DayIndex:Integer):TDateTime;
function IsSummertime(CDate:TDateTime):boolean;
function GetHolidayIndex(CDate:TDateTime; Land:TdeState):integer;
function last_phase(date:TDateTime; phase:TMoonPhase):TDateTime;
function next_phase(date:TDateTime; phase:TMoonPhase):TDateTime;
function age_of_moon(date: TDateTime):extended;
function current_phase(date:TDateTime):extended;
function moon_distance(date: TDateTime): extended;
function lunation(date:TDateTime):integer;
function sun_distance(date: TDateTime): extended;
function nextperigee(date:TDateTime):TDateTime;
function nextapogee(date:TDateTime):TDateTime;
function StartSeason(year: integer; season:TSeason):TDateTime;
function Eclipse(var date:TDateTime; sun:boolean):TEclipse;
function NextEclipse(var date:TDateTime; sun:boolean):TEclipse;
function Sun_Rise(date:TDateTime; latitude, longitude:extended):TDateTime;
function Sun_Set(date:TDateTime; latitude, longitude:extended):TDateTime;
function Sun_Transit(date:TDateTime; latitude, longitude:extended):TDateTime;
function Moon_Rise(date:TDateTime; latitude, longitude:extended):TDateTime;
function Moon_Set(date:TDateTime; latitude, longitude:extended):TDateTime;
function Moon_Transit(date:TDateTime; latitude, longitude:extended):TDateTime;
function GetSternzeichen(const CDat:TDateTime):Integer;
function GetSeason(CDat:TDate):String;
function StartDateSeason(CDat:TDate; Season:TSeason): TDate;
function GetMoonPhase(const Dat:TDateTime): String;
procedure SetRelatPosition;  
procedure EDrive(DtTm : TDateTime; Lat, Long : TGeoAngle; Sun : boolean);
function GetMoonPhaseNum(Age : extended):Integer;
function GetMoonPhaseConst(Age : extended):TMoonPhase8;
function DateTimeReal(DtTm : TDateTime; Long : extended):TDateTime;
function DateTimeTZ(DtTm : TDateTime; Hour : ShortInt):TDateTime;


implementation

uses Math;

const
  AU=149597869;
  mean_lunation=29.530589;
  earth_radius=6378.15;

type
t_coord = record
  longitude, latitude, radius: extended;
  rektaszension, declination: extended;
  parallax: extended;
  end;
  T_RiseSet=(_rise,_set,_transit);


var
  FirstWeekDay      : Integer = 2;  // 2 : montag (nach DIN 1355)
  FirstWeekDate     : Integer = 4;  // 1 : beginnt am ersten januar
                                    // 4 : erste-4 tage-woche (nach DIN 1355)
                                    // 7 : erste volle woche
  LocaleIDate       : Integer;
  LocaleILDate      : Integer;
  CurrentYear2Digit : Integer;
  CurrentCentury    : Integer;

function TimeZoneBias:longint;
var
  tz_info: TTimeZoneInformation;
begin
  case GetTimeZoneInformation(tz_info) of
    1: result:=tz_info.StandardBias+tz_info.Bias;
    2: result:=tz_info.DaylightBias+tz_info.Bias;
    else result:=0;
  end;
end;

function TransformDate(CDatum:TDatetime): TDateTime;
var
  bias, jetzt,
  start_time,
  first_now   : TDateTime;
begin
  first_now:=now;

  bias:=TimeZoneBias/(60*24);
  start_time:=CDatum-(now-first_now)-bias;

  jetzt:=(now-first_now)+start_time;
  result:=jetzt;
end;

function put_in_360(x:extended):extended;
begin
  result:=x-round(x/360)*360;
  while result<0 do result:=result+360;
  end;

function sin_d(x:extended):extended;
begin
  sin_d:=sin(put_in_360(x)*pi/180);
end;

function cos_d(x:extended):extended;
begin
  cos_d:=cos(put_in_360(x)*pi/180);
end;

function tan_d(x:extended):extended;
begin
  tan_d:=tan(put_in_360(x)*pi/180);
end;

function arctan2_d(a,b:extended):extended;
begin
  result:=arctan2(a,b)*180/pi;
end;

function arcsin_d(x:extended):extended;
begin
  result:=arcsin(x)*180/pi;
end;

function arccos_d(x:extended):extended;
begin
  result:=arccos(x)*180/pi;
end;

function arctan_d(x:extended):extended;
begin
  result:=arctan(x)*180/pi
end;

function Julian_Date(CDate:TDateTime):extended;
var FirstJulian,
    First2k : TDateTime;
begin
  FirstJulian := EncodeDate(1582, 10, 15);
  if CDate >= FirstJulian then begin
    First2k := EncodeDate(2000, 1, 1);
    Result := 2451544.5 - First2k + CDate;
  end
  else
    Result:=0;
end;

function Delphi_Date(JulDat:extended):TDateTime;
var FirstJulian,
    First2k : TDateTime;
begin
  FirstJulian := EncodeDate(1582, 10, 15);
  if JulDat >= Julian_Date(FirstJulian) then
   begin
    First2k := EncodeDate(2000, 1, 1);
    Result := JulDat - 2451544.5 + First2k;
   end
    else
     Result:=0;
end;

function Star_Time(CDate:TDateTime):extended;
var
  jd, t: extended;
begin
  jd := Julian_Date(CDate);
  t := (jd - 2451545.0) / 36525;
  Result := Put_in_360(280.46061837 + 360.98564736629 * (jd - 2451545.0)
                       + t * t * (0.000387933 - t / 38710000));
end;

function GetDayOfYear(CDate:TDateTime):Integer;
var T,M,J     : word;
    CFirstDay : TDateTime;
begin
  try
    DecodeDate(CDate,J,M,T);
    CFirstDay :=EncodeDate(J,1,1);
    Result :=Trunc(CDate-CFirstDay+1);
  except
    Result:=0;
  end;
end;

function GetDaysPerYear(CDate:TDateTime):Integer;
var nYear, nMonth, nDays : word;
begin
  try
    DecodeDate(CDate, nYear, nMonth, nDays);
    CDate := EncodeDate(nYear, 12, 31);
    Result :=GetDayOfYear(CDate);
  except
    Result:=0;
  end;
end;

function IsLeapYear(CDate:TDateTime):boolean;
var nYear, nMonth, nDay : word;
begin
  DecodeDate(CDate, nYear, nMonth, nDay);
  Result:=(nYear mod 4 = 0) and ((nYear mod 100 <> 0) or (nYear mod 400 = 0));
end;

function GetWeekOfYear(CDate:TDateTime):Integer;
var nYear, nMonth, nDay : word;
begin
  CDate:=CDate-((DayOfWeek(CDate)-FirstWeekDay+7) mod 7)+ 7-FirstWeekDate;
  DecodeDate(CDate, nYear, nMonth, nDay);
  Result:=(Trunc(CDate-EncodeDate(nYear, 1, 1)) div 7)+1;
end;

function GetWeeksPerYear(CDate:TDateTime):Integer;
var cWeek : byte;
    nYear, nMonth, nDay : word;
begin
  DecodeDate(CDate, nYear, nMonth, nDay);
  cWeek:=GetWeekOfYear(EncodeDate(nYear,12,31));
  if cWeek=1 then
    Result:=52
  else
    Result:=cWeek;
end;

function GetDaysPerMonth(CDate:TDateTime):Integer;
const
  DaysInMonth: array [1..12] of Integer =
   (31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
var nYear, nMonth, nDay : word;
begin
  DecodeDate(CDate, nYear, nMonth, nDay);
  Result:=DaysInMonth[nMonth];
  if (nMonth=2) and IsLeapYear(nYear) then
    Inc(Result);
end;

function GetFirstDayOfWeek(CDate:TDateTime; DayIndex:Integer):TDateTime;
begin
  if DayIndex<1 then
   DayIndex := FirstWeekDay;
  while DayOfWeek(CDate) <> DayIndex do
   CDate := CDate-1;
    Result := CDate;
end;

function GetDay(CDate:TDate):Integer;
begin
 Result := DayOfWeek(CDate);
end;

function IsSummertime(CDate:TDateTime):boolean;
var nYear, nMonth, nDay   : word;
    nBeginn, nEnde   : TDateTime;
begin
  try
    CDate := Trunc(CDate);
    DecodeDate(CDate, nYear, nMonth, nDay);
    if nYear < 1980 then // keine SZ vor 1980
      Result:=false
    else
     begin // Beginn SZ
      nBeginn := EncodeDate(nYear, 3, 31);
      while DayOfWeek(nBeginn) <> 1 do
        nBeginn := nBeginn-1;
      // Ende SZ
      if nYear <= 1995
       then // ... 1995: letzter So in 09.
        nEnde := EncodeDate(nYear, 9, 30)
       else // 1996 ...: letzter So in 10
        nEnde:=EncodeDate(nYear, 10, 31);
      while DayOfWeek(nEnde) <> 1 do
        nEnde := nEnde-1;
      Result:= (CDate >= nBeginn) and (CDate < nEnde);
    end;
  except
    Result:=false;
  end;
end;

procedure InitLocale;
var SystemTime: TSystemTime;

  function GetLocaleInt(CInfo:LCTYPE):Integer;
  var
    Buffer: array[0..1] of Char;
  begin
    if GetLocaleInfo(GetThreadLocale, CInfo, Buffer, 2) > 0
     then  Result:=Ord(Buffer[0])-Ord('0')
      else Result := -1;
  end;

begin
  LocaleIDate  := GetLocaleInt(LOCALE_IDATE);
  LocaleILDate := GetLocaleInt(LOCALE_ILDATE);
  GetLocalTime(SystemTime);
  CurrentYear2Digit := SystemTime.wYear mod 100;
  CurrentCentury    := SystemTime.wYear - CurrentYear2Digit;
end;

function GetHolidayIndex(CDate:TDateTime; Land:TdeState):integer;
var CurrY,Y,M,
    D,dw,OM,aw : word;
    Dat,
    Ostern,
    Weihnacht  : TDateTime;

  function OsterSonntag(Y:word):TDateTime;
  var oa,ob,oc,od,oe,
      Tag,Monat : integer;
  begin
    oa := y mod 19 ;
    ob := y mod 4;
    oc := y mod 7;
    od := (19 * oa + 24) mod 30;
    oe := (2 * ob + 4 * oc + 6 * od + 5) mod 7;
    Tag := 22 + od + oe;
    Monat := 3;
    if Tag > 31 then
     begin
      Tag := od + oe -9;
      Monat := 4;
     end;
    if (Tag = 26) and (Monat = 4) then Tag := 19;
    if (Tag = 25) and (Monat = 4) and (od = 28) and (oe = 6) and (oa > 10) then Tag := 18;
    try
      Result:= EncodeDate(y, Monat, Tag);
    except
      Result:=0;
    end;
  end; // ende ostersonntag

begin
  Result:=0;
  try
    DecodeDate(CDate, Y, M, D);
  except
    Exit;
  end;
  CurrY:=GetDayOfYear(CDate);
  if CurrY=GetDayOfYear(FSpringDate) then Result := -21;
  if CurrY=GetDayOfYear(FSummerDate) then Result := -22;
  if CurrY=GetDayOfYear(FAutumnDate) then Result := -23;
  if CurrY=GetDayOfYear(FWinterDate) then Result := -24;
  if (D >= 1) and (M >= 1) and (M <= 12) and (Y >= 1900) then
   begin
    Ostern:=OsterSonntag(Y);
    try
      DecodeDate(Ostern, Y, OM, D);
    except
      OM := 4;
    end;
    try
      Weihnacht:=EncodeDate(Y, 12, 25);
      if (DayOfWeek(Weihnacht) -1) =0
       then dw:=7
        else dw:=DayOfWeek(Weihnacht)-1;
    except
      Weihnacht:=-1;
      dw:=0;
    end;

    Dat := EncodeDate(Y, 2, 2);
    if CurrY=GetDayOfYear(Dat) then Result := -1;
    Dat := Encodedate(Y, 2, 14);
    if CurrY=GetDayOfYear(Dat) then Result := -2;
    Dat := Ostern -45;
    while DayOfWeek(Dat) <> 2 do Dat := Dat -1;
    if CurrY = GetDayOfYear(Dat-4) then Result := -3;
    if CurrY=GetDayOfYear(Dat) then Result := -4;
    if CurrY=GetDayOfYear(Dat+1) then Result := -5;
      if CurrY=GetDayOfYear(Dat+2) then Result := -6;
    Dat:=Encodedate(Y, 3, 25);
    if CurrY=GetDayOfYear(Dat) then Result := -7;
    if CurrY=GetDayOfYear(Ostern-7) then Result := -8;
    if CurrY=GetDayOfYear(Ostern-3) then Result := -9;
    Dat:=EncodeDate(y, 4, 30);
    aw:=DayOfWeek(Dat)-1;
    Dat:=Dat - aw + 14;
    if Dat=(Ostern + 49) then Dat := Dat-7;
    if CurrY=GetDayOfYear(Dat) then Result := -10;
    Dat:=Encodedate(Y, 6, 29);
    if CurrY=GetDayOfYear(Dat) then Result := -11;
    Dat:=Encodedate(Y, 9, 8);
    if CurrY=GetDayOfYear(Dat) then Result := -12;
    Dat:=Encodedate(Y, 9, 29);
    while DayOfWeek(Dat)<>1 do Dat := Dat + 1;
    if CurrY=GetDayOfYear(Dat) then Result := -13;
    Dat:=Encodedate(Y, 12, 8);
    if CurrY=GetDayOfYear(Dat) then Result := -14;
    Dat:=Encodedate(Y, 12, 31);
    if CurrY=GetDayOfYear(Dat) then Result := -15;
    Dat:=Weihnacht-1;
    while DayOfWeek(Dat)<>1 do Dat := Dat - 1;
    if CurrY=GetDayOfYear(Dat-21) then Result := -16;
    if CurrY=GetDayOfYear(Dat-14) then Result := -17;
    if CurrY=GetDayOfYear(Dat-7) then Result := -18;
    if CurrY=GetDayOfYear(Dat) then Result := -19;
    if CurrY=GetDayOfYear(Weihnacht-1) then Result := -20;
    if (Weihnacht >= 0) and (CurrY=GetDayOfYear(Weihnacht - dw - 28)) then Result := -25;
    if (Weihnacht >= 0) and (CurrY=GetDayOfYear(Weihnacht - dw - 35)) then Result := -26;

    if CurrY=1 then Result := 1;
    Dat:=EncodeDate(Y, 5, 1);
    if CurrY=GetDayOfYear(Dat) then Result := 2;
    Dat:=EncodeDate(Y, 10, 3);
    if CurrY=GetDayOfYear(Dat) then Result := 3;
    if (Land=de_Baden_Wuerttemberg) or (Land=de_Bayern) or (Land=de_Nordrhein_Westfalen)
     or (Land=de_Rheinland_Pfalz) or (Land=de_Saarland) then
      begin
       Dat:=EncodeDate(Y, 11, 1);
       if CurrY=GetDayOfYear(Dat) then Result := 4;
     end;
    if (Weihnacht >= 0) and (CurrY=GetDayOfYear(Weihnacht)) then Result := 5;
    if (Weihnacht >= 0) and (CurrY=GetDayOfYear(Weihnacht+1)) then Result := 6;
    if CurrY=GetDayOfYear(Ostern-2)  then Result := 7;
    if CurrY=GetDayOfYear(Ostern)    then Result := 8;
    if CurrY=GetDayOfYear(Ostern+1)  then Result := 9;
    if CurrY=GetDayOfYear(Ostern+39) then Result := 10;
    if CurrY=GetDayOfYear(Ostern+49) then Result := 11;
    if CurrY=GetDayOfYear(Ostern+50) then Result := 12;
    if (Land<>de_None) and ((Land<de_Berlin) or (Land=de_Hessen) or
     ((Land>=de_Nordrhein_Westfalen) and (Land<de_Sachsen)) or
     (Land=de_Thueringen)) then
      if CurrY=GetDayOfYear(Ostern+60) then Result := 13;
    if (Land=de_Baden_Wuerttemberg) or (Land=de_Bayern) or (Land=de_Sachsen_Anhalt) then
      if CurrY=6 then Result := 14;
    if (Land=de_Bayern) or (Land=de_Saarland) then
     begin
      Dat:=EncodeDate(Y, 8, 15);
      if CurrY=GetDayOfYear(Dat) then Result := 15;
     end;
    if (Land=de_Brandenburg) or (Land=de_Mecklenburg_Vorpommern) or
     (Land=de_Sachsen) or (Land=de_Sachsen_Anhalt) or (Land=de_Thueringen) then
      begin
       Dat:=Encodedate(Y, 10, 31);
       if CurrY=GetDayOfYear(Dat) then Result := 16;
     end;
    if (Weihnacht >= 0) and (Land=de_Sachsen)
     and (CurrY=GetDayOfYear(Weihnacht - dw - 32)) then Result := 17;
  end;
end;

procedure calc_geocentric(var coord:t_coord; date:TDateTime);
const
arg_mul:array[0..30, 0..4] of shortint =
  (( 0, 0, 0, 0, 1),
   (-2, 0, 0, 2, 2),
   ( 0, 0, 0, 2, 2),
   ( 0, 0, 0, 0, 2),
   ( 0, 1, 0, 0, 0),
   ( 0, 0, 1, 0, 0),
   (-2, 1, 0, 2, 2),
   ( 0, 0, 0, 2, 1),
   ( 0, 0, 1, 2, 2),
   (-2,-1, 0, 2, 2),
   (-2, 0, 1, 0, 0),
   (-2, 0, 0, 2, 1),
   ( 0, 0,-1, 2, 2),
   ( 2, 0, 0, 0, 0),
   ( 0, 0, 1, 0, 1),
   ( 2, 0,-1, 2, 2),
   ( 0, 0,-1, 0, 1),
   ( 0, 0, 1, 2, 1),
   (-2, 0, 2, 0, 0),
   ( 0, 0,-2, 2, 1),
   ( 2, 0, 0, 2, 2),
   ( 0, 0, 2, 2, 2),
   ( 0, 0, 2, 0, 0),
   (-2, 0, 1, 2, 2),
   ( 0, 0, 0, 2, 0),
   (-2, 0, 0, 2, 0),
   ( 0, 0,-1, 2, 1),
   ( 0, 2, 0, 0, 0),
   ( 2, 0,-1, 0, 1),
   (-2, 2, 0, 2, 2),
   ( 0, 1, 0, 0, 1));

arg_phi:array[0..30, 0..1] of longint =
  ((-171996,-1742),
   ( -13187,  -16),
   (  -2274,   -2),
   (   2062,    2),
   (   1426,  -34),
   (    712,    1),
   (   -517,   12),
   (   -386,   -4),
   (   -301,    0),
   (    217,   -5),
   (   -158,    0),
   (    129,    1),
   (    123,    0),
   (     63,    0),
   (     63,    1),
   (    -59,    0),
   (    -58,   -1),
   (    -51,    0),
   (     48,    0),
   (     46,    0),
   (    -38,    0),
   (    -31,    0),
   (     29,    0),
   (     29,    0),
   (     26,    0),
   (    -22,    0),
   (     21,    0),
   (     17,   -1),
   (     16,    0),
   (    -16,    1),
   (    -15,    0));

arg_eps:array[0..30, 0..1] of longint =
  (( 92025,   89),
   (  5736,  -31),
   (   977,   -5),
   (  -895,    5),
   (    54,   -1),
   (    -7,    0),
   (   224,   -6),
   (   200,    0),
   (   129,   -1),
   (   -95,    3),
   (     0,    0),
   (   -70,    0),
   (   -53,    0),
   (     0,    0),
   (   -33,    0),
   (    26,    0),
   (    32,    0),
   (    27,    0),
   (     0,    0),
   (   -24,    0),
   (    16,    0),
   (    13,    0),
   (     0,    0),
   (   -12,    0),
   (     0,    0),
   (     0,    0),
   (   -10,    0),
   (     0,    0),
   (    -8,    0),
   (     7,    0),
   (     9,    0));

var
  t,omega: extended;
  d,m,ms,f,s: extended;
  i: integer;
  epsilon,epsilon_0,delta_epsilon: extended;
  delta_phi: extended;
  alpha,delta: extended;
begin
  t:=(julian_date(date) - 2451545.0) / 36525;
  omega := put_in_360(125.04452 + (-1934.136261 + (0.0020708 + 1 / 450000 * t) * t) * t);
  d := put_in_360(297.85036 + (445267.111480+ (-0.0019142 + t / 189474) * t ) * t);
  m := put_in_360(357.52772 + (35999.050340 + (-0.0001603 - t / 300000) * t ) * t);
  ms := put_in_360(134.96298+ (477198.867398+ (0.0086972 + t / 56250) * t) * t);
  f := put_in_360(93.27191+ (483202.017538+ (-0.0036825 + t / 327270) * t) * t);

  delta_phi:=0;
  delta_epsilon:=0;

  for i := 0 to 30 do begin
    s:= arg_mul[i, 0] * d
       +arg_mul[i, 1] * m
       +arg_mul[i, 2] * ms
       +arg_mul[i, 3] * f
       +arg_mul[i, 4] * omega;
    delta_phi := delta_phi + (arg_phi[i, 0] + arg_phi[i, 1] * t * 0.1) * sin_d(s);
    delta_epsilon := delta_epsilon+(arg_eps[i, 0] + arg_eps[i, 1] * t * 0.1) * cos_d(s);
    end;

  delta_phi := delta_phi * 0.0001 / 3600;
  delta_epsilon := delta_epsilon * 0.0001 / 3600;

  epsilon_0 := 84381.448+ (-46.8150 + (-0.00059 + 0.001813 * t) * t) * t;
  epsilon := (epsilon_0 + delta_epsilon) / 3600;
  coord.longitude := put_in_360(coord.longitude + delta_phi);

  alpha := arctan2_d( sin_d(coord.longitude) * cos_d(epsilon)
                   - tan_d(coord.latitude) * sin_d(epsilon)
                  , cos_d(coord.longitude));
  delta := arcsin_d( sin_d(coord.latitude) * cos_d(epsilon)
                   + cos_d(coord.latitude) * sin_d(epsilon) * sin_d(coord.longitude));

  coord.rektaszension := alpha;
  coord.declination := delta;
end;

function sun_coordinate(date:TDateTime):t_coord;
var
  t,e,m,c,nu: extended;
  l0,o,omega,lambda: extended;
begin
  t := (julian_date(date) - 2451545.0) / 36525;
  l0 := 280.46645 + (36000.76983 + 0.0003032 * t) * t;
  e := 0.016708617 + (-0.000042037 - 0.0000001236 * t) * t;
  m := 357.52910 + (35999.05030 - (0.0001559 + 0.00000048 * t) * t) *t;
  c := (1.914600 + (-0.004817 - 0.000014 * t) * t) * sin_d(m)
       + (0.019993 - 0.000101 * t) * sin_d(2*m)
       + 0.000290 * sin_d(3*m);
  o := put_in_360(l0 + c);
  nu := m + c;
  result.radius := (1.000001018 * (1 - e * e)) / (1 + e * cos_d(nu)) *AU;
  omega := 125.04452 + (-1934.136261 + (0.0020708 + 1 / 450000 * t) * t) * t;
  lambda := put_in_360(o - 0.00569 - 0.00478 * sin_d(omega)
                      - 20.4898 / 3600 / (result.radius / AU));
  result.longitude := lambda;
  result.latitude := 0;
  calc_geocentric(result, date);
end;

function moon_coordinate(date:TDateTime):t_coord;
const
arg_lr:array[0..59, 0..3] of shortint =
  (( 0, 0, 1, 0),
   ( 2, 0,-1, 0),
   ( 2, 0, 0, 0),
   ( 0, 0, 2, 0),
   ( 0, 1, 0, 0),
   ( 0, 0, 0, 2),
   ( 2, 0,-2, 0),
   ( 2,-1,-1, 0),
   ( 2, 0, 1, 0),
   ( 2,-1, 0, 0),
   ( 0, 1,-1, 0),
   ( 1, 0, 0, 0),
   ( 0, 1, 1, 0),
   ( 2, 0, 0,-2),
   ( 0, 0, 1, 2),
   ( 0, 0, 1,-2),
   ( 4, 0,-1, 0),
   ( 0, 0, 3, 0),
   ( 4, 0,-2, 0),
   ( 2, 1,-1, 0),
   ( 2, 1, 0, 0),
   ( 1, 0,-1, 0),
   ( 1, 1, 0, 0),
   ( 2,-1, 1, 0),
   ( 2, 0, 2, 0),
   ( 4, 0, 0, 0),
   ( 2, 0,-3, 0),
   ( 0, 1,-2, 0),
   ( 2, 0,-1, 2),
   ( 2,-1,-2, 0),
   ( 1, 0, 1, 0),
   ( 2,-2, 0, 0),
   ( 0, 1, 2, 0),
   ( 0, 2, 0, 0),
   ( 2,-2,-1, 0),
   ( 2, 0, 1,-2),
   ( 2, 0, 0, 2),
   ( 4,-1,-1, 0),
   ( 0, 0, 2, 2),
   ( 3, 0,-1, 0),
   ( 2, 1, 1, 0),
   ( 4,-1,-2, 0),
   ( 0, 2,-1, 0),
   ( 2, 2,-1, 0),
   ( 2, 1,-2, 0),
   ( 2,-1, 0,-2),
   ( 4, 0, 1, 0),
   ( 0, 0, 4, 0),
   ( 4,-1, 0, 0),
   ( 1, 0,-2, 0),
   ( 2, 1, 0,-2),
   ( 0, 0, 2,-2),
   ( 1, 1, 1, 0),
   ( 3, 0,-2, 0),
   ( 4, 0,-3, 0),
   ( 2,-1, 2, 0),
   ( 0, 2, 1, 0),
   ( 1, 1,-1, 0),
   ( 2, 0, 3, 0),
   ( 2, 0,-1,-2));

arg_b:array[0..59, 0..3] of shortint =
  (( 0, 0, 0, 1),
   ( 0, 0, 1, 1),
   ( 0, 0, 1,-1),
   ( 2, 0, 0,-1),
   ( 2, 0,-1, 1),
   ( 2, 0,-1,-1),
   ( 2, 0, 0, 1),
   ( 0, 0, 2, 1),
   ( 2, 0, 1,-1),
   ( 0, 0, 2,-1),
   ( 2,-1, 0,-1),
   ( 2, 0,-2,-1),
   ( 2, 0, 1, 1),
   ( 2, 1, 0,-1),
   ( 2,-1,-1, 1),
   ( 2,-1, 0, 1),
   ( 2,-1,-1,-1),
   ( 0, 1,-1,-1),
   ( 4, 0,-1,-1),
   ( 0, 1, 0, 1),
   ( 0, 0, 0, 3),
   ( 0, 1,-1, 1),
   ( 1, 0, 0, 1),
   ( 0, 1, 1, 1),
   ( 0, 1, 1,-1),
   ( 0, 1, 0,-1),
   ( 1, 0, 0,-1),
   ( 0, 0, 3, 1),
   ( 4, 0, 0,-1),
   ( 4, 0,-1, 1),
   ( 0, 0, 1,-3),
   ( 4, 0,-2, 1),
   ( 2, 0, 0,-3),
   ( 2, 0, 2,-1),
   ( 2,-1, 1,-1),
   ( 2, 0,-2, 1),
   ( 0, 0, 3,-1),
   ( 2, 0, 2, 1),
   ( 2, 0,-3,-1),
   ( 2, 1,-1, 1),
   ( 2, 1, 0, 1),
   ( 4, 0, 0, 1),
   ( 2,-1, 1, 1),
   ( 2,-2, 0,-1),
   ( 0, 0, 1, 3),
   ( 2, 1, 1,-1),
   ( 1, 1, 0,-1),
   ( 1, 1, 0, 1),
   ( 0, 1,-2,-1),
   ( 2, 1,-1,-1),
   ( 1, 0, 1, 1),
   ( 2,-1,-2,-1),
   ( 0, 1, 2, 1),
   ( 4, 0,-2,-1),
   ( 4,-1,-1,-1),
   ( 1, 0, 1,-1),
   ( 4, 0, 1,-1),
   ( 1, 0,-1,-1),
   ( 4,-1, 0,-1),
   ( 2,-2, 0, 1));

sigma_r: array[0..59] of longint =
 (-20905355,
   -3699111,
   -2955968,
    -569925,
      48888,
      -3149,
     246158,
    -152138,
    -170733,
    -204586,
    -129620,
     108743,
     104755,
      10321,
          0,
      79661,
     -34782,
     -23210,
     -21636,
      24208,
      30824,
      -8379,
     -16675,
     -12831,
     -10445,
     -11650,
      14403,
      -7003,
          0,
      10056,
       6322,
      -9884,
       5751,
          0,
      -4950,
       4130,
          0,
      -3958,
          0,
       3258,
       2616,
      -1897,
      -2117,
       2354,
          0,
          0,
      -1423,
      -1117,
      -1571,
      -1739,
          0,
      -4421,
          0,
          0,
          0,
          0,
       1165,
          0,
          0,
       8752);

sigma_l: array[0..59] of longint =
 (6288774,
  1274027,
   658314,
   213618,
  -185116,
  -114332,
    58793,
    57066,
    53322,
    45758,
   -40923,
   -34720,
   -30383,
    15327,
   -12528,
    10980,
    10675,
    10034,
     8548,
    -7888,
    -6766,
    -5163,
     4987,
     4036,
     3994,
     3861,
     3665,
    -2689,
    -2602,
     2390,
    -2348,
     2236,
    -2120,
    -2069,
     2048,
    -1773,
    -1595,
     1215,
    -1110,
     -892,
     -810,
      759,
     -713,
     -700,
      691,
      596,
      549,
      537,
      520,
     -487,
     -399,
     -381,
      351,
     -340,
      330,
      327,
     -323,
      299,
      294,
        0);

sigma_b: array[0..59] of longint =
 (5128122,
   280602,
   277693,
   173237,
    55413,
    46271,
    32573,
    17198,
     9266,
     8822,
     8216,
     4324,
     4200,
    -3359,
     2463,
     2211,
     2065,
    -1870,
     1828,
    -1794,
    -1749,
    -1565,
    -1491,
    -1475,
    -1410,
    -1344,
    -1335,
     1107,
     1021,
      833,
      777,
      671,
      607,
      596,
      491,
     -451,
      439,
      422,
      421,
     -366,
     -351,
      331,
      315,
      302,
     -283,
     -229,
      223,
      223,
     -220,
     -220,
     -185,
      181,
     -177,
      176,
      166,
     -164,
      132,
     -119,
      115,
      107);

var
  t,d,m,ms,f,e,ls : extended;
  sr,sl,sb,temp: extended;
  a1,a2,a3: extended;
  lambda,beta,delta: extended;
  i: integer;
begin
  t := (julian_date(date) - 2451545) / 36525;
  d := 297.8502042 + (445267.1115168 + (-0.0016300 + (1 / 545868 - 1 / 113065000 * t) * t) * t) * t;
  m := 357.5291092 + (35999.0502909 + (-0.0001536 + 1 / 24490000 * t) * t) * t;
  ms := 134.9634114 + (477198.8676313 + (0.0089970 + (1 / 69699 - 1 / 1471200 * t) * t) * t) * t;
  f := 93.2720993 + (483202.0175273 + (-0.0034029 + (-1 / 3526000 + 1 / 863310000 * t) * t) * t) * t;
  e := 1.0 + (-0.002516 - 0.0000074 * t) * t;
  ls := 218.3164591 + (481267.88134236 + (-0.0013268 + (1 / 538841 - 1 / 65194000 * t) * t) * t) * t;
  a1 := 119.75 + 131.849 * t;
  a2 := 53.09 + 479264.290 * t;
  a3 := 313.45 + 481266.484 * t;

  sr := 0;
  for i := 0 to 59 do begin
    temp := sigma_r[i] * cos_d( arg_lr[i, 0] *d
                              + arg_lr[i, 1] * m
                              + arg_lr[i, 2] * ms
                              + arg_lr[i, 3] * f);
    if abs(arg_lr[i, 1]) = 1 then temp := temp * e;
    if abs(arg_lr[i, 1]) = 2 then temp := temp * e;
    sr := sr + temp;
    end;

  sl := 0;
  for i := 0 to 59 do begin
    temp := sigma_l[i] * sin_d( arg_lr[i, 0] * d
                              + arg_lr[i, 1] * m
                              + arg_lr[i, 2] * ms
                              + arg_lr[i, 3] * f);
    if abs(arg_lr[i,1]) = 1 then temp := temp * e;
    if abs(arg_lr[i,1]) = 2 then temp := temp * e;
    sl := sl + temp;
  end;

  sl := sl + 3958 * sin_d(a1)
           + 1962 * sin_d(ls-f)
           + 318 * sin_d(a2);

  sb := 0;
  for i := 0 to 59 do begin
    temp := sigma_b[i] * sin_d( arg_b[i, 0] * d
                              + arg_b[i, 1] * m
                              + arg_b[i, 2] * ms
                              + arg_b[i, 3] * f);
    if abs(arg_b[i,1]) = 1 then temp := temp * e;
    if abs(arg_b[i,1]) = 2 then temp := temp * e;
    sb := sb + temp;
  end;

  sb := sb -2235 * sin_d(ls)
           + 382 * sin_d(a3)
           + 175 * sin_d(a1-f)
           + 175 * sin_d(a1+f)
           + 127 * sin_d(ls-ms)
           - 115 * sin_d(ls+ms);

  lambda := ls + sl / 1000000;
  beta := sb / 1000000;
  delta:=385000.56+sr/1000;

  result.radius := delta;
  result.longitude := lambda;
  result.latitude := beta;

  calc_geocentric(result,date);
end;

procedure calc_phase_data(date:TDateTime; phase:TMoonPhase; var jde,kk,m,ms,f,o,e: extended);
var
  t: extended;
  k: longint;
  ts: extended;
begin
  k   := round((date - encodedate(2000, 1, 1)) / 36525.0 * 1236.85);
  ts  := (date - encodedate(2000, 1, 1)) / 36525.0;
  kk  := int(k) + ord(phase) / 4.0;
  t   := kk / 1236.85;
  jde := 2451550.09765 + 29.530588853 * kk
         + t * t * (0.0001337 - t * (0.000000150 -0.00000000073 * t));
  m   := 2.5534 + 29.10535669 * kk - t * t * (0.0000218 + 0.00000011 * t);
  ms  := 201.5643 + 385.81693528 * kk + t * t * (0.1017438 + t * (0.00001239 - t * 0.000000058));
  f   := 160.7108 + 390.67050274 * kk - t * t * (0.0016341 + t * (0.00000227 - t * 0.000000011));
  o   := 124.7746 - 1.56375580 * kk + t * t * (0.0020691 + t * 0.00000215);
  e   := 1 - ts * (0.002516 + ts * 0.0000074);
end;

function nextphase(date:TDateTime; phase:TMoonPhase):TDateTime;
var
  t: extended;
  kk: extended;
  jde: extended;
  m,ms,f,o,e: extended;
  korr,w,akorr: extended;
  a:array[1..14] of extended;
begin
  calc_phase_data(date, phase, jde, kk, m, ms, f, o, e);
  t := kk / 1236.85;
  case phase of
    mpNewmoon: begin
     korr := - 0.40720 * sin_d(ms)
             + 0.17241 * e * sin_d(m)
             + 0.01608 * sin_d(2 * ms)
             + 0.01039 * sin_d(2 * f)
             + 0.00739 * e * sin_d(ms-m)
             - 0.00514 * e * sin_d(ms+m)
             + 0.00208 * e * e * sin_d(2 * m)
             - 0.00111 * sin_d(ms-2 * f)
             - 0.00057 * sin_d(ms + 2 * f)
             + 0.00056 * e * sin_d(2 * ms + m)
             - 0.00042 * sin_d(3 * ms)
             + 0.00042 * e * sin_d(m + 2 * f)
             + 0.00038 * e * sin_d(m - 2 * f)
             - 0.00024 * e * sin_d(2 * ms - m)
             - 0.00017 * sin_d(o)
             - 0.00007 * sin_d(ms + 2 * m)
             + 0.00004 * sin_d(2 * ms - 2 * f)
             + 0.00004 * sin_d(3 * m)
             + 0.00003 * sin_d(ms + m - 2 * f)
             + 0.00003 * sin_d(2 * ms + 2 * f)
             - 0.00003 * sin_d(ms + m + 2 * f)
             + 0.00003 * sin_d(ms - m + 2 * f)
             - 0.00002 * sin_d(ms - m - 2 * f)
             - 0.00002 * sin_d(3 * ms + m)
             + 0.00002 * sin_d(4 * ms);
    end;
    mpFirstQuarter, mpLastQuarter: begin
     korr := - 0.62801 * sin_d(ms)
             + 0.17172 * e * sin_d(m)
             - 0.01183 * e * sin_d(ms+m)
             + 0.00862 * sin_d(2 * ms)
             + 0.00804 * sin_d(2 * f)
             + 0.00454 * e * sin_d(ms - m)
             + 0.00204 * e * e * sin_d(2 * m)
             - 0.00180 * sin_d(ms - 2 * f)
             - 0.00070 * sin_d(ms + 2 * f)
             - 0.00040 * sin_d(3 * ms)
             - 0.00034 * e * sin_d(2 * ms - m)
             + 0.00032 * e * sin_d(m + 2 * f)
             + 0.00032 * e * sin_d(m - 2 * f)
             - 0.00028 * e * e * sin_d(ms + 2 * m)
             + 0.00027 * e * sin_d(2 * ms + m)
             - 0.00017 * sin_d(o)
             - 0.00005 * sin_d(ms - m - 2 * f)
             + 0.00004 * sin_d(2 * ms + 2 * f)
             - 0.00004 * sin_d(ms + m + 2 * f)
             + 0.00004 * sin_d(ms - 2 * m)
             + 0.00003 * sin_d(ms + m - 2 * f)
             + 0.00003 * sin_d(3 * m)
             + 0.00002 * sin_d(2 * ms - 2 * f)
             + 0.00002 * sin_d(ms - m + 2 * f)
             - 0.00002 * sin_d(3 * ms + m);
      w := 0.00306 - 0.00038 * e * cos_d(m)
                   + 0.00026 * cos_d(ms)
                   - 0.00002 * cos_d(ms - m)
                   + 0.00002 * cos_d(ms + m)
                   + 0.00002 * cos_d(2 * f);
      if phase = mpFirstQuarter then korr := korr + w else korr := korr - w;
    end;
    mpFullmoon: begin
     korr := - 0.40614 * sin_d(ms)
             + 0.17302 * e * sin_d(m)
             + 0.01614 * sin_d(2 * ms)
             + 0.01043 * sin_d(2 * f)
             + 0.00734 * e * sin_d(ms-m)
             - 0.00515 * e * sin_d(ms+m)
             + 0.00209 * e * e * sin_d(2 * m)
             - 0.00111 * sin_d(ms - 2 * f)
             - 0.00057 * sin_d(ms + 2 * f)
             + 0.00056 * e * sin_d(2 * ms + m)
             - 0.00042 * sin_d(3 * ms)
             + 0.00042 * e * sin_d(m + 2 * f)
             + 0.00038 * e * sin_d(m - 2 * f)
             - 0.00024 * e * sin_d(2 * ms - m)
             - 0.00017 * sin_d(o)
             - 0.00007 * sin_d(ms + 2 * m)
             + 0.00004 * sin_d(2 * ms - 2 * f)
             + 0.00004 * sin_d(3 * m)
             + 0.00003 * sin_d(ms + m - 2 * f)
             + 0.00003 * sin_d(2 * ms + 2 * f)
             - 0.00003 * sin_d(ms + m + 2 * f)
             + 0.00003 * sin_d(ms - m + 2 * f)
             - 0.00002 * sin_d(ms - m - 2 * f)
             - 0.00002 * sin_d(3 * ms + m)
             + 0.00002 * sin_d(4 * ms);
    end;
    else
      korr := 0;
  end;

  a[1]  := 299.77 + 0.107408 * kk - 0.009173 * t * t;
  a[2]  := 251.88 + 0.016321 * kk;
  a[3]  := 251.83 + 26.651886 * kk;
  a[4]  := 349.42 + 36.412478 * kk;
  a[5]  := 84.66 + 18.206239 * kk;
  a[6]  := 141.74 + 53.303771 * kk;
  a[7]  := 207.14 + 2.453732 * kk;
  a[8]  := 154.84 + 7.306860 * kk;
  a[9]  := 34.52+27.261239 * kk;
  a[10] := 207.19 + 0.121824 * kk;
  a[11] := 291.34 + 1.844379 * kk;
  a[12] := 161.72 + 24.198154 * kk;
  a[13] := 239.56 + 25.513099 * kk;
  a[14] := 331.55 + 3.592518 * kk;
  akorr :=   + 0.000325 * sin_d(a[1])
             + 0.000165 * sin_d(a[2])
             + 0.000164 * sin_d(a[3])
             + 0.000126 * sin_d(a[4])
             + 0.000110 * sin_d(a[5])
             + 0.000062 * sin_d(a[6])
             + 0.000060 * sin_d(a[7])
             + 0.000056 * sin_d(a[8])
             + 0.000047 * sin_d(a[9])
             + 0.000042 * sin_d(a[10])
             + 0.000040 * sin_d(a[11])
             + 0.000037 * sin_d(a[12])
             + 0.000035 * sin_d(a[13])
             + 0.000023 * sin_d(a[14]);
  korr := korr+ akorr;
  nextphase := delphi_date(jde + korr);
end;

function last_phase(date:TDateTime; phase:TMoonPhase):TDateTime;
var
  temp_date: TDateTime;
begin
  temp_date := date + 28;
  result := temp_date;
  while result > date do
   begin
    result := nextphase(temp_date, phase);
    temp_date := temp_date - 28;
   end;
end;

function next_phase(date:TDateTime; phase:TMoonPhase):TDateTime;
var
  temp_date: TDateTime;
begin
  temp_date := date - 28;
  result := temp_date;
  while result < date do
   begin
    result := nextphase(temp_date, phase);
    temp_date := temp_date + 28;
   end;
end;

function moon_phase_angle(date: TDateTime):extended;
var
  sun_coord,moon_coord: t_coord;
  phi,i: extended;
begin
  sun_coord  := sun_coordinate(date);
  moon_coord := moon_coordinate(date);
  phi := arccos(cos_d(moon_coord.latitude)
               * cos_d(moon_coord.longitude - sun_coord.longitude));
  i := arctan(sun_coord.radius*sin(phi) /
             (moon_coord.radius - sun_coord.radius*cos(phi)));
  if i < 0 then  result := i / pi * 180 + 180
            else  result := i / pi * 180;

  if put_in_360(moon_coord.longitude-sun_coord.longitude) > 180
   then result := -result;

end;

function age_of_moon(date: TDateTime):extended;
var
  sun_coord,moon_coord: t_coord;
begin
  sun_coord  := sun_coordinate(date);
  moon_coord := moon_coordinate(date);
  result := put_in_360(moon_coord.longitude - sun_coord.longitude) / 360 * mean_lunation;
end;

function current_phase(date:TDateTime):extended;
begin
  result := (1 + cos_d(moon_phase_angle(date))) / 2;
end;

function moon_distance(date: TDateTime): extended;
begin
  result := moon_coordinate(date).radius;
end;

function lunation(date:TDateTime):integer;
begin
  result := round((last_phase(date, mpNewMoon) - delphi_date(2423436)) / mean_lunation) + 1;
end;

function sun_distance(date: TDateTime): extended;
begin
  result := sun_coordinate(date).radius / AU;
end;

function sun_diameter(date:TDateTime):extended;
begin
  result := 959.63 / (sun_coordinate(date).radius / AU) * 2;
end;

function moon_diameter(date:TDateTime):extended;
begin
  result := 358473400 / moon_coordinate(date).radius * 2;
end;

function nextXXXgee(date:TDateTime; apo: boolean):TDateTime;
const
arg_apo:array[0..31, 0..2] of shortint =
   { D  F  M }
  (( 2, 0, 0),
   ( 4, 0, 0),
   ( 0, 0, 1),
   ( 2, 0,-1),
   ( 0, 2, 0),
   ( 1, 0, 0),
   ( 6, 0, 0),
   ( 4, 0,-1),
   ( 2, 2, 0),
   ( 1, 0, 1),
   ( 8, 0, 0),
   ( 6, 0,-1),
   ( 2,-2, 0),
   ( 2, 0,-2),
   ( 3, 0, 0),
   ( 4, 2, 0),
   ( 8, 0,-1),
   ( 4, 0,-2),
   (10, 0, 0),
   ( 3, 0, 1),
   ( 0, 0, 2),
   ( 2, 0, 1),
   ( 2, 0, 2),
   ( 6, 2, 0),
   ( 6, 0,-2),
   (10, 0,-1),
   ( 5, 0, 0),
   ( 4,-2, 0),
   ( 0, 2, 1),
   (12, 0, 0),
   ( 2, 2,-1),
   ( 1, 0,-1));

arg_per:array[0..59, 0..2] of shortint =
   { D  F  M }
  (( 2, 0, 0),
   ( 4, 0, 0),
   ( 6, 0, 0),
   ( 8, 0, 0),
   ( 2, 0,-1),
   ( 0, 0, 1),
   (10, 0, 0),
   ( 4, 0,-1),
   ( 6, 0,-1),
   (12, 0, 0),
   ( 1, 0, 0),
   ( 8, 0,-1),
   (14, 0, 0),
   ( 0, 2, 0),
   ( 3, 0, 0),
   (10, 0,-1),
   (16, 0, 0),
   (12, 0,-1),
   ( 5, 0, 0),
   ( 2, 2, 0),
   (18, 0, 0),
   (14, 0,-1),
   ( 7, 0, 0),
   ( 2, 1, 0),
   (20, 0, 0),
   ( 1, 0, 1),
   (16, 0,-1),
   ( 4, 0, 1),
   ( 2, 0,-2),
   ( 4, 0,-2),
   ( 6, 0,-2),
   (22, 0, 0),
   (18, 0,-1),
   ( 6, 0, 1),
   (11, 0, 0),
   ( 8, 0, 1),
   ( 4,-2, 0),
   ( 6, 2, 0),
   ( 3, 0, 1),
   ( 5, 0, 1),
   (13, 0, 0),
   (20, 0,-1),
   ( 3, 0, 2),
   ( 4, 2,-2),
   ( 1, 0, 2),
   (22, 0,-1),
   ( 0, 4, 0),
   ( 6,-2, 0),
   ( 2,-2, 1),
   ( 0, 0, 2),
   ( 0, 2,-1),
   ( 2, 4, 0),
   ( 0, 2,-2),
   ( 2,-2, 2),
   (24, 0, 0),
   ( 4,-4, 0),
   ( 9, 0, 0),
   ( 4, 2, 0),
   ( 2, 0, 2),
   ( 1, 0,-1));

koe_apo:array[0..31, 0..1] of longint =
   {    1   T }
  (( 4392,  0),
   (  684,  0),
   (  456,-11),
   (  426,-11),
   (  212,  0),
   ( -189,  0),
   (  144,  0),
   (  113,  0),
   (   47,  0),
   (   36,  0),
   (   35,  0),
   (   34,  0),
   (  -34,  0),
   (   22,  0),
   (  -17,  0),
   (   13,  0),
   (   11,  0),
   (   10,  0),
   (    9,  0),
   (    7,  0),
   (    6,  0),
   (    5,  0),
   (    5,  0),
   (    4,  0),
   (    4,  0),
   (    4,  0),
   (   -4,  0),
   (   -4,  0),
   (    3,  0),
   (    3,  0),
   (    3,  0),
   (   -3,  0));

koe_per:array[0..59, 0..1] of longint =
   {     1   T }
  ((-16769,  0),
   (  4589,  0),
   ( -1856,  0),
   (   883,  0),
   (  -773, 19),
   (   502,-13),
   (  -460,  0),
   (   422,-11),
   (  -256,  0),
   (   253,  0),
   (   237,  0),
   (   162,  0),
   (  -145,  0),
   (   129,  0),
   (  -112,  0),
   (  -104,  0),
   (    86,  0),
   (    69,  0),
   (    66,  0),
   (   -53,  0),
   (   -52,  0),
   (   -46,  0),
   (   -41,  0),
   (    40,  0),
   (    32,  0),
   (   -32,  0),
   (    31,  0),
   (   -29,  0),
   (   -27,  0),
   (    24,  0),
   (   -21,  0),
   (   -21,  0),
   (   -21,  0),
   (    19,  0),
   (   -18,  0),
   (   -14,  0),
   (   -14,  0),
   (   -14,  0),
   (    14,  0),
   (   -14,  0),
   (    13,  0),
   (    13,  0),
   (    11,  0),
   (   -11,  0),
   (   -10,  0),
   (    -9,  0),
   (    -8,  0),
   (     8,  0),
   (     8,  0),
   (     7,  0),
   (     7,  0),
   (     7,  0),
   (    -6,  0),
   (    -6,  0),
   (     6,  0),
   (     5,  0),
   (    27,  0),
   (    27,  0),
   (     5,  0),
   (    -4,  0));

var
  k, jde, t: extended;
  d,m,f,v: extended;
  i: integer;
begin
  k := round(((date - encodedate(1999, 1, 1)) / 365.25-0.97) * 13.2555);
  if apo then k := k + 0.5;
  t := k / 1325.55;
  jde := 2451534.6698 + 27.55454988 * k + (-0.0006886 +
        (-0.000001098 + 0.0000000052 * t) * t) * t * t;
  d := 171.9179 + 335.9106046 * k + (-0.0100250 + (-0.00001156 + 0.000000055 * t) * t) * t * t;
  m := 347.3477 + 27.1577721 * k + (-0.0008323 - 0.0000010 * t) * t * t;
  f := 316.6109 + 364.5287911 * k + (-0.0125131 - 0.0000148 * t) * t * t;
  v := 0;
  if apo then
    for i := 0 to 31 do
      v := v + sin_d(arg_apo[i, 0] * d + arg_apo[i, 1] * f + arg_apo[i, 2] * m) *
          (koe_apo[i, 0] * 0.0001 + koe_apo[i, 1] * 0.00001 * t)
  else
    for i := 0 to 59 do
      v := v + sin_d(arg_per[i, 0] * d + arg_per[i, 1] * f + arg_per[i, 2] * m) *
          (koe_per[i, 0] * 0.0001 + koe_per[i, 1] * 0.00001 * t);
  result := delphi_date(jde + v);
end;

function nextperigee(date:TDateTime):TDateTime;
var
  temp_date: TDateTime;
begin
  temp_date := date - 28;
  result := temp_date;
  while result < date do
   begin
    result := nextXXXgee(temp_date, false);
    temp_date := temp_date + 28;
   end;
end;

function nextapogee(date:TDateTime):TDateTime;
var
  temp_date: TDateTime;
begin
  temp_date := date - 28;
  result := temp_date;
  while result < date do
   begin
    result := nextXXXgee(temp_date, true);
    temp_date := temp_date + 28;
   end;
end;

function StartSeason(year: integer; season:TSeason):TDateTime;
var
  y: extended;
  jde0: extended;
  t, w, dl, s: extended;
  i: integer;
const
a: array[0..23] of integer =
  (485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, 50,
   45, 44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8 );

bc : array[0..23, 1..2] of extended =
  (( 324.96,   1934.136 ),
   ( 337.23,  32964.467 ),
   ( 342.08,     20.186 ),
   (  27.85, 445267.112 ),
   (  73.14,  45036.886 ),
   ( 171.52,  22518.443 ),
   ( 222.54,  65928.934 ),
   ( 296.72,   3034.906 ),
   ( 243.58,   9037.513 ),
   ( 119.81,  33718.147 ),
   ( 297.17,    150.678 ),
   (  21.02,   2281.226 ),
   ( 247.54,  29929.562 ),
   ( 325.15,  31555.956 ),
   (  60.93,   4443.417 ),
   ( 155.12,  67555.328 ),
   ( 288.79,   4562.452 ),
   ( 198.04,  62894.029 ),
   ( 199.76,  31436.921 ),
   (  95.39,  14577.848 ),
   ( 287.11,  31931.756 ),
   ( 320.81,  34777.259 ),
   ( 227.73,   1222.114 ),
   (  15.45,  16859.074 ));

begin
  case year of
    -1000..+999: begin
      y := year / 1000;
      case season of
        seSpring: jde0 := 1721139.29189 + (365242.13740 + ( 0.06134 + ( 0.00111-0.00071 * y) * y) * y) * y;
        seSummer: jde0 := 1721223.25401 + (365241.72562 + (-0.05323 + ( 0.00907+0.00025 * y) * y) * y) * y;
        seAutumn: jde0 := 1721325.70455 + (365242.49558 + (-0.11677 + (-0.00297+0.00074 * y) * y) * y) * y;
        seWinter: jde0 := 1721414.39987 + (365242.88257 + (-0.00769 + (-0.00933-0.00006 * y) * y) * y) * y;
        else    jde0 := 0;
      end;
    end;
    +1000..+3000: begin
      y := (year - 2000) / 1000;
      case season of
        seSpring: jde0 := 2451623.80984 + (365242.37404 + ( 0.05169 + (-0.00411 - 0.00057 * y) * y) * y) * y;
        seSummer: jde0 := 2451716.56767 + (365241.62603 + ( 0.00325 + ( 0.00888 - 0.00030 * y) * y) * y) * y;
        seAutumn: jde0 := 2451810.21715 + (365242.01767 + (-0.11575 + ( 0.00337 + 0.00078 * y) * y) * y) * y;
        seWinter: jde0 := 2451900.05952 + (365242.74049 + (-0.06223 + (-0.00823 + 0.00032 * y) * y) * y) * y;
        else    jde0 := 0;
      end;
    end;
    else raise E_OutOfAlgorithRange.Create('algorithmischer fehler');
  end;
  t := (jde0 - 2451545.0) / 36525;
  w := 35999.373 * t - 2.47;
  dl := 1 + 0.0334 * cos_d(w) + 0.0007 * cos_d(2 * w);
  s := 0;
  for i := 0 to 23 do
    s := s + a[i] * cos_d(bc[i, 1] + bc[i, 2] * t);

  Result := TransformDate(Delphi_Date(jde0 + (0.00001 * s) / dl));
end;

function Eclipse(var date:TDateTime; sun:boolean):TEclipse;
var
  jde,kk,m,ms,f,o,e: extended;
  t,f1,a1: extended;
  p,q,w,gamma,u: extended;
begin
  if sun then
    calc_phase_data(date, mpNewMoon, jde, kk, m, ms, f, o, e)
  else
    calc_phase_data(date, mpFullMoon, jde, kk, m, ms, f, o, e);
  t := kk / 1236.85;
  if abs(sin_d(f)) > 0.36 then
    result := ecNone
  else
   begin
    f1 := f - 0.02665 * sin_d(o);
    a1 := 299.77 + 0.107408 * kk - 0.009173 * t * t;
    if sun then
      jde := jde - 0.4075     * sin_d(ms)
                 + 0.1721 * e * sin_d(m)
    else
      jde := jde - 0.4065     * sin_d(ms)
                 + 0.1727 * e * sin_d(m);
    jde := jde   + 0.0161     * sin_d(2*ms)
                - 0.0097     * sin_d(2*f1)
                + 0.0073 * e * sin_d(ms-m)
                - 0.0050 * e * sin_d(ms+m)
                - 0.0023     * sin_d(ms-2*f1)
                + 0.0021 * e * sin_d(2*m)
                + 0.0012     * sin_d(ms+2*f1)
                + 0.0006 * e * sin_d(2*ms+m)
                - 0.0004     * sin_d(3*ms)
                - 0.0003 * e * sin_d(m+2*f1)
                + 0.0003     * sin_d(a1)
                - 0.0002 * e * sin_d(m-2*f1)
                - 0.0002 * e * sin_d(2*ms-m)
                - 0.0002     * sin_d(o);
    p :=        + 0.2070 * e * sin_d(m)
                + 0.0024 * e * sin_d(2*m)
                - 0.0392     * sin_d(ms)
                + 0.0116     * sin_d(2*ms)
                - 0.0073 * e * sin_d(ms+m)
                + 0.0067 * e * sin_d(ms-m)
                + 0.0118     * sin_d(2*f1);
    q :=        + 5.2207
                - 0.0048 * e * cos_d(m)
                + 0.0020 * e * cos_d(2*m)
                - 0.3299     * cos_d(ms)
                - 0.0060 * e * cos_d(ms+m)
                + 0.0041 * e * cos_d(ms-m);
    w := abs(cos_d(f1));
    gamma:=(p * cos_d(f1) + q * sin_d(f1)) * (1 - 0.0048 * w);
    u := + 0.0059
         + 0.0046 * e * cos_d(m)
         - 0.0182     * cos_d(ms)
         + 0.0004     * cos_d(2*ms)
         - 0.0005     * cos_d(m+ms);

    if sun then
     begin
      if abs(gamma) < 0.9972 then
       begin
        if u < 0 then
          result := ecTotal
        else if u > 0.0047 then
          result := ecCircular
        else if u < 0.00464 * sqrt(1 - gamma * gamma) then
          result := ecCirculartotal
        else
          result := ecCircular;
      end
      else if abs(gamma) > 1.5433 + u then
        result:= ecNone
      else if abs(gamma) < 0.9972 + abs(u) then
        result := ecNoncentral
      else
        result := ecPartial;
    end
    else
     begin
      if (1.0128 - u - abs(gamma)) / 0.5450 > 0 then
        result := ecTotal
      else if (1.5573 + u - abs(gamma)) / 0.5450 > 0 then
        result := ecHalfshadow
      else
        result := ecNone;
    end;
  end;
  date := delphi_date(jde);
end;

function NextEclipse(var date:TDateTime; sun:boolean):TEclipse;
var
  temp_date: TDateTime;
begin
  result := ecNone;
  temp_date := date - 28 * 2;
  while temp_date < date do
   begin
    temp_date := temp_date + 28;
    result := Eclipse(temp_date,sun);
   end;
  date := temp_date;
end;

procedure correct_position(var position:t_coord; date:TDateTime;
                           latitude,longitude,height:extended);
var
  u,h,delta_alpha: extended;
  rho_sin, rho_cos: extended;
const
  b_a = 0.99664719;
begin
  u := arctan_d(b_a * b_a * tan_d(longitude));
  rho_sin := b_a * sin_d(u) + height / 6378140 * sin_d(longitude);
  rho_cos := cos_d(u) + height / 6378140 * cos_d(longitude);

  position.parallax := arcsin_d(sin_d(8.794 / 3600)/(moon_distance(date) / AU));
  h := star_time(date) - longitude - position.rektaszension;
  delta_alpha := arctan_d(
                        (-rho_cos * sin_d(position.parallax) * sin_d(h)) /
                        (cos_d(position.declination) -
                         rho_cos * sin_d(position.parallax) * cos_d(h)));
  position.rektaszension := position.rektaszension + delta_alpha;
  position.declination := arctan_d(
                                  (( sin_d(position.declination)
                                   - rho_sin * sin_d(position.parallax))
                                  * cos_d(delta_alpha)) /
                                  ( cos_d(position.declination)
                                  - rho_cos * sin_d(position.parallax) * cos_d(h)));
end;

function Calc_Set_Rise(date:TDateTime; latitude, longitude:extended;
                       sun: boolean; kind: T_RiseSet):TDateTime;
var
  h: Extended;
  pos1, pos2, pos3: t_coord;
  h0, theta0, cos_h0, cap_h0: extended;
  m0,m1,m2: extended;

  function interpolation(y1, y2, y3, n: extended):extended;
  var
    a, b, c: extended;
  begin
    a := y2 - y1;
    b := y3 - y2;
    if a > 100 then  a := a - 360;
    if a < -100 then  a := a + 360;
    if b > 100 then  b := b - 360;
    if b < -100 then  b := b + 360;
    c := b - a;
    result := y2 + 0.5 * n * (a + b + n * c);
  end; // end interpolation

  function correction(m:extended; kind:integer):extended;
  var
    alpha,delta,h, height: extended;
  begin
    alpha := interpolation(pos1.rektaszension,
                           pos2.rektaszension,
                           pos3.rektaszension,
                           m);
    delta := interpolation(pos1.declination,
                           pos2.declination,
                           pos3.declination,
                           m);
    h := put_in_360((theta0 + 360.985647 * m) - longitude - alpha);
    if h > 180 then h := h - 360;

    height := arcsin_d(sin_d(latitude) * sin_d(delta)
                      + cos_d(latitude) * cos_d(delta) * cos_d(h));

    case kind of
      0    : result := -h / 360;
      1, 2 : result := (height - h0) / (360 * cos_d(delta) * cos_d(latitude) * sin_d(h));
      else result := 0;
    end;
  end; // end correction

begin
  if sun then
    h0 := -0.8333
  else
   begin
    pos1 := moon_coordinate(date);
    correct_position(pos1, date, latitude, longitude, 0);
    h0 := 0.7275 * pos1.parallax - 34 / 60;
   end;

  h := int(date);
  theta0 := star_time(h);
  if sun then
   begin
    pos1:=sun_coordinate(h-1);
    pos2:=sun_coordinate(h);
    pos3:=sun_coordinate(h+1);
   end
    else
     begin
      pos1 := moon_coordinate(h - 1);
      correct_position(pos1, h-1, latitude, longitude, 0);
      pos2 := moon_coordinate(h);
      correct_position(pos2, h, latitude, longitude, 0);
      pos3 := moon_coordinate(h + 1);
      correct_position(pos3, h + 1, latitude, longitude, 0);
     end;

  cos_h0 := (sin_d(h0) - sin_d(latitude) * sin_d(pos2.declination)) /
            (cos_d(latitude) * cos_d(pos2.declination));
  if (cos_h0 < -1) or (cos_h0 > 1) then
    raise E_NoRiseSet.Create('Werte nicht errechenbar');
  cap_h0 := arccos_d(cos_h0);

  m0 := (pos2.rektaszension + longitude - theta0) / 360;
  m1 := m0 - cap_h0 / 360;
  m2 := m0 + cap_h0 / 360;

  m0 := frac(m0);
  if m0 < 0 then m0 := m0 + 1;
  m1 := frac(m1);
  if m1 < 0 then m1 := m1 + 1;
  m2 := frac(m2);
  if m2 < 0 then m2 := m2+1;

  m0 := m0 + correction(m0, 0);
  m1 := m1 + correction(m1, 1);
  m2 := m2 + correction(m2, 2);

  case kind of
    _rise:    result := h + m1;
    _set:     result := h + m2;
    _transit: result := h + m0;
    else      result := 0;
  end;
end;

function Sun_Rise(date:TDateTime; latitude, longitude:extended):TDateTime;
begin
  result := Calc_Set_Rise(date, latitude, longitude, true, _rise);
end;

function Sun_Set(date:TDateTime; latitude, longitude:extended):TDateTime;
begin
  result := Calc_Set_Rise(date, latitude, longitude, true, _set);
end;

function Sun_Transit(date:TDateTime; latitude, longitude:extended):TDateTime;
begin
  result := Calc_Set_Rise(date, latitude, longitude, true, _transit);
end;

function Moon_Rise(date:TDateTime; latitude, longitude:extended):TDateTime;
begin
  result := Calc_Set_Rise(date, latitude, longitude, false, _rise);
end;

function Moon_Set(date:TDateTime; latitude, longitude:extended):TDateTime;
begin
  result := Calc_Set_Rise(date, latitude, longitude, false, _set);
end;

function Moon_Transit(date:TDateTime; latitude, longitude:extended):TDateTime;
begin
  result:=Calc_Set_Rise(date, latitude, longitude, false, _transit);
end;

procedure EDrive(DtTm : TDateTime; Lat, Long : TGeoAngle; Sun : boolean);
var
 dtJr, dtJt, dtJs, dtJc : extended;
 Y, M, D, th, tm, ts : word;
 ss : word;
 s : String;
 eDate : TDate;
 eTime : TTime;
begin
 eDate := DtTm;
 eTime := DtTm;
 dtJc := Julian_Date(DtTm);
 if Sun then dtJr := Julian_Date(Sun_Rise(DtTm, Lat, Long))
         else dtJr := Julian_Date(Moon_Rise(DtTm, Lat, Long));
 if Sun then dtJt := Julian_Date(Sun_Transit(DtTm, Lat, Long))
         else dtJt := Julian_Date(Moon_Transit(DtTm, Lat, Long));
 if Sun then dtJs := Julian_Date(Sun_Set(DtTm, Lat, Long))
         else dtJs := Julian_Date(Moon_Set(DtTm, Lat, Long));
   if (dtJs < dtJr) and not Sun then
    begin
     DecodeDate(eDate + 1, Y, M, D);
     DecodeTime(eTime, th, tm, ts, ss);
     dtJs := Julian_Date(Moon_Set((EncodeDateTime(Y, M, D, th, tm, ts, ss)), Lat, Long));
    end;
   if (dtJt < dtJr) and not Sun then
    begin
     DecodeDate(eDate + 1, Y, M, D);
     DecodeTime(eTime, th, tm, ts, ss);
     dtJt := Julian_Date(Moon_Transit((EncodeDateTime(Y, M, D, th, tm, ts, ss)), Lat, Long));
    end;

 s := Copy(FloatToStr(dtJr), 7, Length(FloatToStr(dtJr)) -6);
 dtJr := Trunc(StrToFloat(s) * 10000);
 s := Copy(FloatToStr(dtJt), 7, Length(FloatToStr(dtJt)) -6);
 dtJt := Trunc(StrToFloat(s) * 10000);
 s := Copy(FloatToStr(dtJs), 7, Length(FloatToStr(dtJs)) -6);
 dtJs := Trunc(StrToFloat(s) * 10000);
 s := Copy(FloatToStr(dtJc), 7, Length(FloatToStr(dtJc)) -6);
 dtJc := Trunc(StrToFloat(s) * 10000);

 EUp   := Round(dtJt - dtJr);
 EDown := Round(dtJs - dtJt);

end;


function GetSternzeichen(const CDat:TDateTime):Integer;
var Y12 : word;
begin
  Result := 0;
  Y12:=GetDayOfYear(CDat);
  case Y12 of
   1..20      : Result := 12;
   21 .. 49   : Result := 1;
   50 .. 79   : Result := 2;
   80 .. 111  : Result := 3;
   112 .. 141 : Result := 4;
   142 .. 172 : Result := 5;
   173 .. 203 : Result := 6;
   204 .. 235 : Result := 7;
   236 .. 266 : Result := 8;
   267 .. 296 : Result := 9;
   297 .. 326 : Result := 10;
   327 .. 355 : Result := 11;
   356 .. 365 : Result := 12;
  end;
end;

function GetSeason(CDat:TDate):String;
var Y, M, D : word;
begin
  DecodeDate(CDat, Y, M, D);
  FWinterDate := Trunc(StartSeason(Y, seWinter));
  FSpringDate := Trunc(StartSeason(Y, seSpring));
  FSummerDate := Trunc(StartSeason(Y, seSummer));
  FAutumnDate := Trunc(StartSeason(Y, seAutumn));
  CDat := Trunc(CDat);
  Result := '';
  if (CDat >= FWinterDate) or (CDat < FSpringDate) then
    Result := 'Winter';
  if (CDat >= FSpringDate) and (CDat < FSummerDate) then
    Result := 'Frühling';
  if (CDat >= FSummerDate) and (CDat < FAutumnDate) then
    Result := 'Sommer';
  if (CDat >= FAutumnDate) and (CDat < FWinterDate) then
    Result := 'Herbst';
end;

function StartDateSeason(CDat:TDate; Season:TSeason): TDate;
var Y, M, D : word;
begin
 DecodeDate(CDat, Y, M, D);
 case Season of
  seWinter : Result := Trunc(StartSeason(Y, seWinter));
  seSpring : Result := Trunc(StartSeason(Y, seSpring));
  seSummer : Result := Trunc(StartSeason(Y, seSummer));
  seAutumn : Result := Trunc(StartSeason(Y, seAutumn));
 end;
end;

function GetMoonPhase(const Dat:TDateTime): String;
var TimeDiff         : extended;

  function LowestPhase(Dat:TDateTime):extended;
  var Phase   : extended;
      Std     : byte;
  begin
    Result := Current_Phase(trunc(Dat));
    for Std := 1 to 23 do begin
      Phase := Current_Phase(trunc(Dat) + Std / 24);
      if Phase < Result then Result := Phase;
    end;
  end; { LowestPhase }

begin
  FMoonDistance:=Moon_Distance(Dat);
  if LowestPhase(Dat-1)>LowestPhase(Dat) then
   begin
    if LowestPhase(Dat+1)>LowestPhase(Dat)
     then Result := 'Neumond'
      else Result := 'abnehmender Mond';
   end
    else
     begin
      if LowestPhase(Dat+1)<LowestPhase(Dat)
       then Result := 'Vollmond'
        else Result := 'zunehmender Mond';
     end;
end;

procedure SetRelatPosition;
var
 i : Integer;
 x, y : extended;
 alpha : extended;
 radius : extended;
begin
 alpha := 0;
 radius := 1;
 for i := 0 to 1999 do
  begin
   eY[i] := sin(alpha*Pi/180)*radius;
   eX[i] := sin((alpha-90)*Pi/180)*radius;
   alpha := alpha + 0.09;
  end;
 eY[2000] := 0;
 eX[2000] := 1;
end;

function GetMoonPhaseNum(Age : extended):Integer;
begin
 Result := 0;
 case Trunc(Age * 10000) of
  0 .. 5000        : Result := 1;
  5001  .. 50833   : Result := 2;
  50834 .. 96666   : Result := 3;
  96667 .. 142499  : Result := 4;
  142500 .. 152500 : Result := 5;
  152501 .. 198333 : Result := 6;
  198334 .. 244166 : Result := 7;
  244167 .. 289999 : Result := 8;
  290000 .. 295000 : Result := 1;
 end; // end case of
end;

function GetMoonPhaseConst(Age : extended):TMoonPhase8;
begin
 case Trunc(Age * 10000) of
  0 .. 5000        : Result := mp_NewMoon;
  5001  .. 50833   : Result := mp_WaxingMoonFirstEighth;
  50834 .. 96666   : Result := mp_WaxingMoonFirstQuarter;
  96667 .. 142499  : Result := mp_WaxingMoonThirdEighth;
  142500 .. 152500 : Result := mp_FullMoon;
  152501 .. 198333 : Result := mp_WaningMoonFifthEighth;
  198334 .. 244166 : Result := mp_WaningMoonLastQuarter;
  244167 .. 289999 : Result := mp_WaningMoonSevensEighth;
  290000 .. 295000 : Result := mp_NewMoon;
 end; // end case of
end;

function DateTimeReal(DtTm : TDateTime; Long : extended):TDateTime;
var
 h, m, s, ss : word;
 Tm : TTime;
 nSec, nMin, nHour : Integer;
begin
 Result := DtTm;
 nSec := Trunc(240 * Long);
 if nSec > 60 then
  begin
   nMin := Trunc(nSec / 60);
   nSec := (nSec - (nMin * 60));
  end;
 if nMin > 60 then
  begin
   nHour := Trunc(nMin / 60);
   nMin := (nMin - (nHour * 60));
  end
   else
    nHour := 0;
  h  := nHour;
  m  := nMin;
  s  := nSec;
  ss := 0;
  Tm := EncodeTime(h, m, s, ss);
  Result := DtTm + Tm;
end;

function DateTimeTZ(DtTm : TDateTime; Hour : ShortInt):TDateTime;
var
 h, m, s, ss : word;
 Tm : TTime;
begin
 if Hour >= 0 then h  := Hour else h := Hour * -1;
 m  := 0;
 s  := 0;
 ss := 0;
 Tm := EncodeTime(h, m, s, ss);
 if Hour > 0
  then Result := DtTm + Tm
   else Result := DtTm - Tm;
end;

end.                                 

Es gibt über OnlinePackageManager auch die Komponente Mondphasen.

[wp_xyz]

Im Online-Package-Manager findest du DelphiMoon von Andreas Hörstemeier, das auch unter Lazarus funktioniert und auch Sonne und Planeten behandelt. Die Mondphasen können durch ein entsprechendes Mond-Bild dargestellt werden.

In meinem kleinen "Sunrise" Projekt (https://github.com/wp-xyz/sunrise) kommt die Unit solar von Jeff Steinkamp zur Anwendung. Kann allerdings nur Sonne.

[catweasel]

Danke an alle für die Tips.
Das schaue ich mir mal an.

Schöne Woche
Michael

[six1]

@wp_xyz,
ich habe die solar.pas mal probiert.

Code: Alles auswählen

procedure TForm1.SetData;
var
  sun: TSun;
begin
    sun := SolarStuff(trunc(now()), Latitude.Value, Longitude.value);
end;

sun.srise und sun.sset liegen einiges neben dran...
Mach ich da was falsch?

https://www.laengengrad-breitengrad.de/
lat/lon Frankfurt/M
(50.109185, 8.680712)

https://www.sunrise-and-sunset.com/de/s ... rt-am-main
Sonnenaufgang heute 07:52
Sonnenuntergang heute 16:43

sun.srise = 08:18:01
sun.sset = 17:16:17


Dresden (51.052874, 13.738285)
https://www.sunrise-and-sunset.com/de/s ... nd/Dresden
Sonnenaufgang heute 08:08
Sonnenuntergang heute 17:06

sun.srise = 08:41:36
sun.sset = 17:33:11

[sstvmaster]

Bei mir verrechnet sich Solar.pas um jeweils ca. 50 minuten, zuviel.

Soll:
Aufgang: 07:52
Untergang: 16:43

Errechnet:
Aufgang: 08:41
Untergang: 17:33

[wp_xyz]

Meine Ergebnisse:

Frankfurt (50.110922, 8.682127): solar sunrise 8:08:37, sunset 17:06:44;
https://www.sunrise-and-sunset.com/de/s ... rt-am-main: 8:08 bzw 17:06

Dresden (51.050409, 13.737262): solar sunrise 7:51:47, sunset 16:43:09
https://www.sunrise-and-sunset.com/de/s ... nd/Dresden: 7:52, 16:43

Habt ihr die Bogen-Minuten richtig eingegeben? Das war eine bescheuerte idee Grad und Minuten zu trennen... Hab'nun die Minuten rausgenommen und die Grad zu einer Floating-Point-Zahl gemacht. Außerdem muss die Erdhälfte north/south west/east angeklickt sein (ich weiß, das ist umständlich...)

Das Programm hatte ich seinerzeit geschrieben, um mich davon zu überzeugen, dass am 21. Dezember zwar der kürzeste Tag ist, der Tag des spätesten Sonnenaufgangs aber erst Anfang Januar und der Tag des frühesten Sonnenuntergangs schon früh im Dezember ist.

[EDIT]
Ich sehe gerade bei den Koordinaten der eingebauten Orte, dass das Vorzeichen der Längengrade offenbar anders herum definiert wird: positiv nach Westen, negativ nach Osten. Ich meine mich zu erinnern, dass ich über diese (für mich unübliche) Konvention schon einmal gestolpert bin.

[sstvmaster]

@wp_xyz

setze mal bitte das Kampatilitäts flag für ältere Versionen, Danke.

[sstvmaster]

Wo kommmen dann bei dem Testprogramm die 10 Minuten unterschied her? Julian Date?

Code: Alles auswählen

program project1;

uses
  solar, sysutils, dateutils;

var
  sun: TSun;

begin
  sun := SolarStuff(trunc(now()), 51.116665, 13.750000);

  WriteLn('Latitude: 51.116665');
  WriteLn('Logitude: 13.750000');

  WriteLn('Sunrise : ' + TimeToStr( LocalTimeToUniversal(sun.srise) ));
  WriteLn('Noon    : ' + TimeToStr( LocalTimeToUniversal(sun.snoon) ));
  WriteLn('Sunset  : ' + TimeToStr( LocalTimeToUniversal(sun.sset) ));

  ReadLn;
end.  

[wp_xyz]

Mit was vergleichst du?

Trage ich die geographischen Koordinaten in https://keisan.casio.com/exec/system/1224686065 ein, wobei zu berücksichtigen ist, dass die solar-Unit die Längengrade mit dem umgekehrten Vorzeichen versieht (wie oben bemerkt), komme ich genau auf dieselben Werte (in deinem Programm stimmt übrigens die Zeitzonenumrechnung nicht. solar berechnet UTC-Zeiten und daher musst du UniverstalTimeToLocal aufrufen, nicht LocalToUniversalTime):

Code: Alles auswählen

program project1;

uses
  solar, sysutils, dateutils;

var
  sun: TSun;

begin
  sun := SolarStuff(trunc(now()), 51.116665, -13.750000);

  WriteLn('Latitude: 51.116665');
  WriteLn('Logitude: 13.750000');

  WriteLn('Sunrise : UTC ' + TimeToStr(sun.srise) + ', local ' + TimeToStr( UniversalTimeToLocal(sun.srise) ));
  WriteLn('Noon    : UTC ' + TimetoStr(sun.snoon) + ', local ' + TimeToStr( UniversalTimeToLocal(sun.snoon) ));
  WriteLn('Sunset  : UTC ' + TimeToStr(sun.sset) + ', local ' + TimeToStr( UniversalTimeToLocal(sun.sset) ));

  ReadLn;
end. 

[sstvmaster]

Mit was vergleichst du?

Mit deinem Programm?

[sstvmaster]

Hmm, komisch, jetzt habe ich das noch mal so gemacht. Jetzt passt das mit der Zeit. Programm ist das gleiche geblieben.
Sehr seltsam. Danke trotzdem.

[wp_xyz]

Programm ist das gleiche geblieben.

Naja, das Programm, das ich als Antwort auf deines gepostet habe, hat ein Minuszeichen vor dem Längengrad. Damit erhalte ich dasselbe Ergebnis wie auf der Webseite, und daraus schließe ich, dass solar die Längengrade entgegen der üblichen Konvention nach Westen positiv rechnet.

Zwischen dem direkten Aufruf von SolarStuff bei dir und meinem Programm ist der Unterschied, dass ich die üblichen Längengrade gleich negiere, um sie in der von solar benötigten Konvention zu übergeben.

[sstvmaster]

ABER -13° ist wo ganz anders das wäre ja West und nicht Ost.
Dann sollte das Vorzeichen in der Funktion geändert werden.

Siehe hier: https://www.google.de/maps/place/51°06' ... 169,3.25z/

Ost ist immer positiv vom Nullmeridian und positiv vom Äquator noch Norden.

[wp_xyz]

Ost ist immer positiv vom Nullmeridian

Nein, offenbar nicht für alle Leute, z.B. hier: https://www.esrl.noaa.gov/gmd/grad/solcalc/azel.html ("Long: East=- West=+"). Ich könnte mir durchaus vorstellen, dass man es im Land jenseits des Atlantik nicht gerne sieht wenn ihnen "negative" Koordinaten ("pfui") zugeordnet werden und man das Vorzeichen einfach umgedreht hat.