unit uMoon;

 {$i ah_def.inc }

{ Copyright 1997-2001 Andreas H鰎stemeier            Version 2.0 2001-07-07 }
{ this component is public domain - please check the file moon.hlp for       }
{ more detailed info on usage and distributing                               }

{ Algorithms taken from the book "Astronomical Algorithms" by Jean Meeus     }

(*$b-*)   { I may make use of the shortcut boolean eval }

(*@/// interface *)
interface


(*@/// uses *)
uses
  uAhMath,
  uVsop,
  sysutils;
(*@\\\*)

type
  TMoonPhase=(Newmoon,WaxingCrescrent,FirstQuarter,WaxingGibbous,
              Fullmoon,WaningGibbous,LastQuarter,WaningCrescent);
  TSeason=(Winter,Spring,Summer,Autumn);
  TEclipse=(none, partial, noncentral, circular, circulartotal, total, halfshadow);
  E_NoRiseSet=class(Exception);
  E_OutOfAlgorithmRange=class(Exception);
  TSolarTerm=(st_z2,st_j3,st_z3,st_j4,st_z4,st_j5,st_z5,st_j6,st_z6,
              st_j7,st_z7,st_j8,st_z8,st_j9,st_z9,st_j10,st_z10,
              st_j11,st_z11,st_j12,st_z12,st_j1,st_z1,st_j2);

  //屬性 '鼠','牛','虎','兔','龍','蛇','馬','羊','猴','雞','狗','豬'
  TChineseZodiac=(ch_rat,ch_ox,ch_tiger,ch_rabbit,ch_dragon,ch_snake,
                      ch_horse,ch_goat,ch_monkey,ch_chicken,ch_dog,ch_pig);
  //'甲','乙','丙','丁','戊','己','庚','辛','寅','癸'
  TChineseStem=(ch_jia,ch_yi,ch_bing,ch_ding,ch_wu,ch_ji,
                    ch_geng,ch_xin,ch_ren,ch_gui);
  (*@/// TChineseCycle= record *)
  TChineseCycle=record
    zodiac: TChineseZodiac;
    stem:   TChineseStem;
    end;
  (*@\\\*)
  (*@/// TChineseDate = record *)
  TChineseDate = record
    cycle: integer;
    year: integer;
    epoch_years: integer;
    month: integer;
    leap: boolean;
    leapyear: boolean;
    day: integer;
    yearcycle: TChineseCycle;
    daycycle: TChineseCycle;
    monthcycle: TChineseCycle;
    end;
  (*@\\\*)

const
  (* Date of calendar reformation - start of gregorian calendar *)
  CHINESEZODIACEX:array[0..11] of string = ('鼠','牛','虎','兔','龍','蛇','馬',
        '羊','猴','雞','狗','豬');

  CHINESEZODIAC:array[0..11] of string = ('子','丑','寅','卯','辰','巳','午',
        '未','申','酉','戌','亥');

  CHINESEDAYNAME:array[0..29] of string = ('初一','初二','初三','初四',
        '初五','初六','初七','初八','初九','初十','十一','十二','十三',
        '十四','十五','十六','十七','十八','十九','二十','廿一','廿二',
        '廿三','廿四','廿五','廿六','廿七','廿八','廿九','卅');
        
  CHINESEMONTHNAME:array[0..11] of string = ('正月','二月','三月','四月',
        '五月','六月','七月','八月','九月','十月','十一月','十二月');

  CHINESESTEM:array[0..9] of string = ('甲','乙','丙','丁','戊','己','庚','辛','寅','癸');

  CHINESEWEEKDAY:array[1..7] of string = ('一','二','三','四','五','六','日');
  CHINESELONGWEEKDAY:array[1..7] of string = ('星期一','二','三','四','五','六','日');
  calendar_change_standard: extended = 2299160.5;
  calendar_change_russia: extended = 2421638.5;
  calendar_change_england: extended = 2361221.5;
  calendar_change_sweden: extended = 2361389.5;
  (*@/// Jewish_Month_Name:array[1..13] of string *)
  Jewish_Month_Name:array[1..13] of string = (
    'Nisan',
    'Iyar',
    'Sivan',
    'Tammuz',
    'Av',
    'Elul',
    'Tishri',
    'Heshvan',
    'Kislev',
    'Tevet',
    'Shevat',
    'Adar',
    'Adar 2'
    );
  (*@\\\*)


{ Calendar algorithms }
function julian_date(date:TDateTime):extended;
function delphi_date(juldat:extended):TDateTime;
function EasterDate(year:integer):TDateTime;
function EasterDateJulian(year:integer):TDateTime;
function PesachDate(year:integer):TDateTime;
procedure DecodeDateJewish(date: TDateTime; var year,month,day: word);
function EncodeDateJewish(year,month,day: word):TDateTime;
function WeekNumber(date:TDateTime):integer;

{ Convert date to julian date and back }
function Calc_Julian_date_julian(year,month,day:word):extended;
function Calc_Julian_date_gregorian(year,month,day:word):extended;
function Calc_Julian_date_switch(year,month,day:word; switch_date:extended):extended;
function Calc_Julian_date(year,month,day:word):extended;
procedure Calc_Calendar_date_julian(juldat:extended; var year,month,day:word);
procedure Calc_Calendar_date_gregorian(juldat:extended; var year,month,day:word);
procedure Calc_Calendar_date_switch(juldat:extended; var year,month,day:word; switch_date:extended);
procedure Calc_Calendar_date(juldat:extended; var year,month,day:word);

{ corrected TDateTime functions }
function isleapyearcorrect(year:word):boolean;
function EncodedateCorrect(year,month,day: word):TDateTime;
procedure DecodedateCorrect(date:TDateTime; var year,month,day: word);
procedure DecodetimeCorrect(date:TDateTime; var hour,min,sec,msec: word);
function FalsifyTdateTime(date:TDateTime):TdateTime;

{ Sun and Moon }
function sun_distance(date:TDateTime): extended;
function moon_distance(date:TDateTime): extended;
function age_of_moon(date:TDateTime): extended;

function last_phase(date:TDateTime; phase:TMoonPhase):TDateTime;
function next_phase(date:TDateTime; phase:TMoonPhase):TDateTime;
function nearest_phase(date: TDateTime):TMoonPhase;
function next_blue_moon(date: TDateTime):TDateTime;
function is_blue_moon(lunation: integer):boolean;

function moon_phase_angle(date: TDateTime):extended;
function current_phase(date:TDateTime):extended;
function lunation(date:TDateTime):integer;

function sun_diameter(date:TDateTime):extended;
function moon_diameter(date:TDateTime):extended;

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 Morning_Twilight_Civil(date:TDateTime; latitude, longitude:extended):TDateTime;
function Evening_Twilight_Civil(date:TDateTime; latitude, longitude:extended):TDateTime;
function Morning_Twilight_Nautical(date:TDateTime; latitude, longitude:extended):TDateTime;
function Evening_Twilight_Nautical(date:TDateTime; latitude, longitude:extended):TDateTime;
function Morning_Twilight_Astronomical(date:TDateTime; latitude, longitude:extended):TDateTime;
function Evening_Twilight_Astronomical(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 nextperigee(date:TDateTime):TDateTime;
function nextapogee(date:TDateTime):TDateTime;
function nextperihel(date:TDateTime):TDateTime;
function nextaphel(date:TDateTime):TDateTime;

function NextEclipse(var date:TDateTime; sun:boolean):TEclipse;

procedure Moon_Position_Horizontal(date:TdateTime; longitude,latitude: extended; var elevation,azimuth: extended);
procedure Sun_Position_Horizontal(date:TdateTime; longitude,latitude: extended; var elevation,azimuth: extended);

{ Further useful functions }
function star_time(date:TDateTime):extended;
function StartSeason(year: integer; season:TSeason):TDateTime;
function CalcSolarTerm(year: integer; term: TSolarTerm):TDateTime;

{ Chinese calendar }
function ChineseNewYear(year: integer): TDateTime;
function ChineseDate(date: TdateTime): TChineseDate;
function EncodeDateChinese(date:TChineseDate):TDateTime;
(*@\\\0000000301*)
(*@/// implementation *)
implementation

(*$undef low_accuracy *)

const
  AU=149597869;             (* astronomical unit in km *)
  mean_lunation=29.530589;  (* Mean length of a month *)
  tropic_year=365.242190;   (* Tropic year length *)
  earth_radius=6378.15;     (* Radius of the earth *)

(*$ifdef delphi_ge_3 *)
var
(*$else *)
const
(*$endif *)
  (* Shortcuts to avoid calling Encodedate too often *)
  datetime_2000_01_01: extended = 0;
  datetime_1999_01_01: extended = 0;
  datetime_chinese_epoch: extended = 0;
  datetime_first_lunation: extended = 0;
  julian_offset: extended = 0;
  (* How broken is the TDateTime? *)
  negative_dates_broken: boolean = false;
  calendar_reform_supported: boolean = true;
  julian_calendar_before_1582: boolean = true;

const
  beijing_longitude = -(116+25/60);

type
  (*@/// t_coord = record *)
  t_coord = record
    longitude, latitude, radius: extended; (* lambda, beta, R *)
    rektaszension, declination: extended;  (* alpha, delta *)
    parallax: extended;
    elevation, azimuth: extended;          (* h, A *)
    end;
  (*@\\\*)
  T_RiseSet=(_rise,_set,_transit,_rise_civil,_rise_nautical,_rise_astro,_set_civil,_set_nautical,_set_astro);
  TJewishYearStyle=(ys_common_deficient,ys_common_regular,ys_common_complete,
                    ys_leap_deficient,ys_leap_regular,ys_leap_complete);
const
  (*@/// Jewish_Month_length:array[1..13,TJewishYearStyle] of shortint *)
  Jewish_Month_length:array[1..13,TJewishYearStyle] of word = (
   ( 30,30,30,30,30,30),
   ( 29,29,29,29,29,29),
   ( 30,30,30,30,30,30),
   ( 29,29,29,29,29,29),
   ( 30,30,30,30,30,30),
   ( 29,29,29,29,29,29),
   ( 30,30,30,30,30,30),
   ( 29,29,30,29,29,30),
   ( 29,30,30,29,30,30),
   ( 29,29,29,29,29,29),
   ( 30,30,30,30,30,30),
   ( 29,29,29,30,30,30),
   (  0, 0, 0,29,29,29)
   );
  (*@\\\*)
  (*@/// Jewish_Month_Name_short:array[1..13] of string *)
  Jewish_Month_Name_short:array[1..13] of string = (
    'Nis',
    'Iya',
    'Siv',
    'Tam',
    'Av' ,
    'Elu',
    'Tis',
    'Hes',
    'Kis',
    'Tev',
    'She',
    'Ada',
    'Ad2'
    );
  (*@\\\*)
  Jewish_year_length:array[TJewishYearStyle] of integer = (353,354,355,383,384,385);

{ Julian date }
(*@/// function julian_date(date:TDateTime):extended; *)
function julian_date(date:TDateTime):extended;
begin
  julian_date:=julian_offset+date
  end;
(*@\\\*)
(*@/// function delphi_date(juldat:extended):TDateTime; *)
function delphi_date(juldat:extended):TDateTime;
begin
  delphi_date:=juldat-julian_offset;
  end;
(*@\\\*)
(*@/// function isleapyearcorrect(year:word):boolean; *)
function isleapyearcorrect(year:word):boolean;
begin
  if year<=1582 then
    result:=((year mod 4)=0)
  else
    result:=(((year mod 4)=0) and ((year mod 100)<>0)) or
             ((year mod 400)=0);
  end;
(*@\\\*)
(*@/// function Calc_Julian_date_julian(year,month,day:word):extended; *)
function Calc_Julian_date_julian(year,month,day:word):extended;
begin
  if (year<1) or (year>9999) then
    raise EConvertError.Create('Invalid year');
  if month<3 then begin
    month:=month+12;
    year:=year-1;
    end;
  case month of
    3,5,7,8,10,12,13: if (day<1) or (day>31) then EConvertError.Create('Invalid day');
    4,6,9,11:         if (day<1) or (day>30) then EConvertError.Create('Invalid day');
    14: case day of
          1..28: ;
          29: if (year+1) mod 4<>0 then EConvertError.Create('Invalid day');
          else EConvertError.Create('Invalid day');
          end;
     else raise EConvertError.Create('Invalid month');
     end;
  result:=trunc(365.25*(year+4716))+trunc(30.6001*(month+1))+day-1524.5;
  end;
(*@\\\*)
(*@/// function Calc_Julian_date_gregorian(year,month,day:word):extended; *)
function Calc_Julian_date_gregorian(year,month,day:word):extended;
var
  a,b: longint;
begin
  if (year<1) or (year>9999) then
    raise EConvertError.Create('Invalid year');
  if month<3 then begin
    month:=month+12;
    year:=year-1;
    end;
  a:=year div 100;
  case month of
    3,5,7,8,10,12,13: if (day<1) or (day>31) then EConvertError.Create('Invalid day');
    4,6,9,11:         if (day<1) or (day>30) then EConvertError.Create('Invalid day');
    14: case day of
          1..28: ;
          29: if (((year mod 4)<>0) or ((year mod 100)=0)) and
                 ((year mod 400)<>0) then EConvertError.Create('Invalid day');
          else EConvertError.Create('Invalid day');
          end;
     else raise EConvertError.Create('Invalid month');
     end;
  b:=2-a+(a div 4);
  result:=trunc(365.25*(year+4716))+trunc(30.6001*(month+1))+day+b-1524.5;
  end;
(*@\\\*)
(*@/// function Calc_Julian_date_switch(year,month,day:word; switch_date:extended):extended; *)
function Calc_Julian_date_switch(year,month,day:word; switch_date:extended):extended;
begin
  result:=Calc_Julian_date_julian(year,month,day);
  if result>=switch_date then begin
    result:=Calc_Julian_date_gregorian(year,month,day);
    if result<switch_date then
      raise EConvertError.Create('Date invalid due to calendar change');
    end;
  end;
(*@\\\*)
(*@/// function Calc_Julian_date(year,month,day:word):extended; *)
function Calc_Julian_date(year,month,day:word):extended;
begin
  result:=Calc_Julian_date_switch(year,month,day,calendar_change_standard);
  end;
(*@\\\*)
(*@/// procedure Calc_Calendar_date_julian(juldat:extended; var year,month,day:word); *)
procedure Calc_Calendar_date_julian(juldat:extended; var year,month,day:word);
var
  z,a,b,c,d,e: longint;
begin
  if juldat<0 then
    raise EConvertError.Create('Negative julian dates not supported');
  juldat:=juldat+0.5;
  z:=trunc(juldat);
  a:=z;
  b:=a+1524;
  c:=trunc((b-122.1)/365.25);
  d:=trunc(365.25*c);
  e:=trunc((b-d)/30.6001);
  day:=b-d-trunc(30.6001*e);
  year:=c-4716;
  if e<14 then
    month:=e-1
  else begin
    month:=e-13;
    year:=year+1;
    end;
  end;
(*@\\\*)
(*@/// procedure Calc_Calendar_date_gregorian(juldat:extended; var year,month,day:word); *)
procedure Calc_Calendar_date_gregorian(juldat:extended; var year,month,day:word);
var
  alpha,z,a,b,c,d,e: longint;
begin
  if juldat<0 then
    raise EConvertError.Create('Negative julian dates not supported');
  juldat:=juldat+0.5;
  z:=trunc(juldat);
  alpha:=trunc((z-1867216.25)/36524.25);
  a:=z+1+alpha-trunc(alpha/4);
  b:=a+1524;
  c:=trunc((b-122.1)/365.25);
  d:=trunc(365.25*c);
  e:=trunc((b-d)/30.6001);
  day:=b-d-trunc(30.6001*e);
  year:=c-4716;
  if e<14 then
    month:=e-1
  else begin
    month:=e-13;
    year:=year+1;
    end;
  end;
(*@\\\*)
(*@/// procedure Calc_Calendar_date_switch(juldat:extended; var year,month,day:word; switch_date:extended); *)
procedure Calc_Calendar_date_switch(juldat:extended; var year,month,day:word; switch_date:extended);
begin
  if juldat<0 then
    raise EConvertError.Create('Negative julian dates not supported');
  if juldat<switch_date then
    Calc_Calendar_date_julian(juldat,year,month,day)
  else
    Calc_Calendar_date_gregorian(juldat,year,month,day);
  end;
(*@\\\*)
(*@/// procedure Calc_Calendar_date(juldat:extended; var year,month,day:word); *)
procedure Calc_Calendar_date(juldat:extended; var year,month,day:word);
begin
  Calc_Calendar_date_switch(juldat,year,month,day,calendar_change_standard);
  end;
(*@\\\*)

{ TDateTime correction }
(*@/// procedure check_TDatetime; *)
(* Check how many bugs the TDateTime has compare to julian date *)
procedure check_TDatetime;
var
  h,m,s,ms: word;
  d1,d2: TDateTime;
begin
  (*$ifndef delphi_1 *)        { Delphi 1 did not allow negative values }
  decodetime(-1.9,h,m,s,ms);
  negative_dates_broken:=h>12;
  (*$endif delphi_1 *)
  d1:=EncodeDate(1582,10,15);
  d2:=EncodeDate(1582,10,4);
  calendar_reform_supported:=((d1-d2)=1);
  d1:=EncodeDate(1500,3,1);
  d2:=EncodeDate(1500,2,28);