stastro.pas

条码控件： 一维条码控件 二维条码控件 PDF417Barcode MaxiCodeBarcode

end;

procedure PhasePrim(LD : TStDate; var PhaseArray : TStPhaseArray);
  {-primitive phase calculation}
var
  I,
  D,
  M,
  Y      : Integer;
  K, TD,
  LYear  : Double;

begin
  StDateToDMY(LD, D, M, Y);
  LYear := Y - 0.05;
  K := (LYear - 2000) * 12.3685;
  K := Int(K);
  TD := K / 12.3685 + 2000;
  if TD > Y then
    K := K-1;

{compute phases for each lunar cycle throughout the year}
  for I := 0 to 13 do begin
    GetPhases(K, PhaseArray[I]);
    K := K + 1;
  end;
end;

function GenSearchPhase(SD : TStDate; PV : Byte) : TStLunarRecord;
{searches for the specified phase in the given month/year expressed by SD}
var
  I,
  C,
  LD,
  LM,
  LY,
  TD,
  TM,
  TY    : Integer;
  ADate : TStDate;
  JD    : Double;
  PhaseArray : TStPhaseArray;
begin
  C := 0;
  FillChar(Result, SizeOf(Result), $FF);

  StDateToDMY(SD, LD, LM, LY);
  PhasePrim(SD, PhaseArray);
  for I := LM-1 to LM+1 do begin
    if (PV = 0) then
      JD := PhaseArray[I].NMDate
    else if (PV = 1) then
      JD := PhaseArray[I].FQDate
    else if (PV = 2) then
      JD := PhaseArray[I].FMDate
    else
      JD := PhaseArray[I].LQDate;
    ADate := AstJulianDateToStDate(JD, True);

    StDateToDMY(ADate, TD, TM, TY);
    if TM < LM then
      continue
    else if TM = LM then begin
      Result.T[C].D := ADate;
      Result.T[C].T := Trunc((Frac(JD) + 0.5) * 86400);
      if Result.T[C].T >= SecondsInDay then
        Dec(Result.T[C].T, SecondsInDay);
      Inc(C);
    end;
  end;
end;

function FirstQuarter(D : TStDate) : TStLunarRecord;
  {-compute date/time of FirstQuarter(s)}
begin
  Result := GenSearchPhase(D, 1);
end;

function FullMoon(D : TStDate) : TStLunarRecord;
  {-compute the date/time of FullMoon(s)}
begin
  Result := GenSearchPhase(D, 2);
end;

function LastQuarter(D: TStDate) : TStLunarRecord;
  {-compute the date/time of LastQuarter(s)}
begin
  Result := GenSearchPhase(D, 3);
end;

function NewMoon(D : TStDate) : TStLunarRecord;
  {-compute the date/time of NewMoon(s)}
begin
  Result := GenSearchPhase(D, 0);
end;

function NextPrevPhase(D : TStDate; Ph : Byte;
                       FindNext : Boolean) : TStDateTimeRec;
var
  LD,
  LM,
  LY  : Integer;
  K,
  JD,
  TJD : Double;
  PR  : TStPhaseRecord;
  OK  : Boolean;
begin
   if (D < MinDate) or (D > MaxDate) then begin
     Result.D := BadDate;
     Result.T := BadTime;
     Exit;
   end;

   StDateToDMY(D, LD, LM, LY);
   K := ((LY + LM/12 + LD/365.25) - 2000) * 12.3685 - 0.5;
   if FindNext then
     K := Round(K)-1
   else
     K := Round(K)-2;

   OK := False;
   TJD := AstJulianDate(D);
   repeat
     GetPhases(K, PR);

     if (Ph = 0) then
       JD := PR.NMDate
     else if (Ph = 1) then
       JD := PR.FQDate
     else if (Ph = 2) then
       JD := PR.FMDate
     else
       JD := PR.LQDate;

     if (FindNext) then begin
       if (JD > TJD) then
         OK := True
       else
         K := K + 1;
     end else begin
       if (JD < TJD) then
         OK := True
       else
         K := K - 1;
     end;
   until OK;

   Result.D := AstJulianDateToStDate(JD, True);
   if (Result.D <> BadDate) then begin
     Result.T := Trunc((Frac(JD) + 0.5) * 86400);
     if Result.T >= SecondsInDay then
       Dec(Result.T, SecondsInDay);
   end else
     Result.T := BadTime;
end;

function NextFirstQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest FirstQuarter}
begin
   Result := NextPrevPhase(D, 1, True);
end;

function NextFullMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest FullMoon}
begin
  Result := NextPrevPhase(D, 2, True);
end;

function NextLastQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest LastQuarter}
begin
  Result := NextPrevPhase(D, 3, True);
end;

function NextNewMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest NewMoon}
begin
  Result := NextPrevPhase(D, 0, True);
end;

function PrevFirstQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest FirstQuarter}
begin
  Result := NextPrevPhase(D, 1, False);
end;

function PrevFullMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest FullMoon}
begin
  Result := NextPrevPhase(D, 2, False);
end;

function PrevLastQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest LastQuarter}
begin
  Result := NextPrevPhase(D, 3, False);
end;

function PrevNewMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest NewMoon}
begin
  Result := NextPrevPhase(D, 0, False);
end;

{Calendar procedures/functions}

function SolEqPrim(Y : Integer; K : Byte) : TStDateTimeRec;
{primitive routine for finding equinoxes and solstices}
var
  JD, TD, LY,
  JCent, MA,
  SLong       : Double;

begin
  JD := 0;
  Result.D := BadDate;
  Result.T := BadTime;

{the following algorithm is valid only in the range of [1000..3000 AD]}
{but is limited to returning dates in [MinYear..MaxYear]}
  if (Y < MinYear) or (Y > MaxYear) then
    Exit;

{compute approximate date/time for specified event}
  LY := (Y - 2000) / 1000;
  case K of
    0 : JD := 2451623.80984
            + 365242.37404 * LY
            + 0.05169 * sqr(LY)
            - 0.00411 * LY * sqr(LY)
            - 0.00057 * sqr(sqr(LY));

    1 : JD := 2451716.56767
            + 365241.62603 * LY
            + 0.00325 * sqr(LY)
            + 0.00888 * LY * sqr(LY)
            - 0.00030 * sqr(sqr(LY));

    2 : JD := 2451810.21715
            + 365242.01767 * LY
            - 0.11575 * sqr(LY)
            + 0.00337 * sqr(LY) * LY
            + 0.00078 * sqr(sqr(LY));

    3 : JD := 2451900.05952
            + 365242.74049 * LY
            - 0.06223 * sqr(LY)
            - 0.00823 * LY * sqr(LY)
            + 0.00032 * sqr(sqr(LY));
  end;

{refine date/time by computing corrections due to solar longitude,}
{nutation and abberation. Iterate using the corrected time until}
{correction is less than one minute}
  repeat
    Result.D := AstJulianDateToStDate(JD, True);
    Result.T := Trunc((Frac(JD) + 0.5) * 86400);
    if Result.T >= SecondsInDay then
      Dec(Result.T, SecondsInDay);
    JCent := (JD - 2451545.0)/36525.0;

{approximate solar longitude - no FK5 correction}
    SLong := 280.46645
           + 36000.76983 * JCent
           + 0.0003032 * sqr(JCent);
    SLong := Frac((SLong)/360.0) * 360.0;
    if SLong < 0 then
      SLong := SLong + 360;

{Equation of the center correction}
    MA := 357.52910
        + 35999.05030 * JCent;
    MA := MA/radcor;
    SLong := SLong
           + (1.914600 - 0.004817 * JCent - 0.000014 * sqr(JCent)) * sin(MA)
           + (0.019993 - 0.000101 * JCent) * sin(2*MA);

{approximate nutation}
    TD := 125.04452
        - 1934.136261 * JCent
        + 0.0020708 * sqr(JCent);
    TD := TD/radcor;
    TD := (-17.20 * sin(TD) - 1.32 * sin(2*SLong/radcor))/3600;

{approximate abberation - solar distance is assumed to be 1 A.U.}
    SLong := SLong - (20.4989/3600) + TD;

{correction to compute Julian Date for event}
    TD := 58 * sin((K*90 - SLong)/radcor);
    if abs(TD) >= 0.0005 then
      JD := JD + TD;
  until abs(TD) < 0.0005;
end;

function Solstice(Y, Epoch : Integer; Summer : Boolean) : TStDateTimeRec;
  {-compute the date/time of the summer or winter solstice}
  {if Summer = True,  compute astronomical summer solstice (summer in N. Hem.)}
  {          = False, compute astronomical winter solstice (winter in N. Hem.)}
begin
  Y := CheckYear(Y, Epoch);
  if Summer then
    Result := SolEqPrim(Y, 1)
  else
    Result := SolEqPrim(Y, 3);
end;

function Equinox(Y, Epoch : Integer; Vernal : Boolean) : TStDateTimeRec;
  {-compute the date/time of the vernal/autumnal equinox}
  {if Vernal = True,  compute astronomical vernal equinox (spring in N. Hem.)}
  {          = False, compute astronomical autumnal equinox (fall in N. Hem.)}
begin
  Y := CheckYear(Y, Epoch);
  if Vernal then
    Result := SolEqPrim(Y, 0)
  else
    Result := SolEqPrim(Y, 2);
end;

function Easter(Y : Integer; Epoch : Integer) : TStDate;
  {-compute the date of Easter}
var
  A, B,
  C, D,
  E, F,
  G, H,
  I, K,
  L, M,
  N, P : LongInt;
begin
  Y := CheckYear(Y, Epoch);

  if (Y < MinYear) or (Y > MaxYear) then begin
    Result := BadDate;
    Exit;
  end;

  A := Y mod 19;
  B := Y div 100;
  C := Y mod 100;
  D := B div 4;
  E := B mod 4;
  F := (B+8) div 25;
  G := (B-F+1) div 3;
  H := (19*A + B - D - G + 15) mod 30;
  I := C div 4;
  K := C mod 4;
  L := (32 + 2*E + 2*I - H - K) mod 7;
  M := (A + 11*H + 22*L) div 451;
  N := (H + L - 7*M + 114) div 31;
  P := (H + L - 7*M + 114) mod 31 + 1;

  Result := DMYToStDate(P, N, Y, Epoch);
end;

{Conversion routines}
function HoursMin(RA : Double) : ShortString;
  {-convert RA to formatted hh:mm:ss string}
var
  HR, MR  : Double;
  HS, MS  : string[12];

begin
  if abs(RA) >= 360.0 then
    RA := Frac(RA / 360.0) * 360;
  if RA < 0 then
    RA := RA + 360.0;

  HR := Int(RA / 15.0);
  MR  := Frac(RA / 15.0) * 60;

  Str(MR:5:2, MS);
  if MS = '60.00' then begin
    MS := '00.00';
    HR := HR + 1;
    if HR = 24 then
      HS := '0'
    else
      Str(HR:2:0, HS);
  end else begin
    if MS[1] = ' ' then
      MS[1] := '0';
    Str(Hr:2:0, HS);
  end;
  Result := HS + ' ' + MS;
end;

function DegsMin(DC : Double) : ShortString;
  {-convert Declination to formatted +/-dd:mm:ss string}
var
  DR, MR  : Double;
  DS, MS  : string[12];
begin
  if abs(DC) > 90 then
    DC := Frac(DC / 90.0) * 90.0;

  DR := Int(DC);
  MR := Frac(abs(DC)) * 60;

  Str(MR:4:1, MS);
  if MS = '60.0' then begin
    MS := '00.0';
    if DC >= 0 then
      DR := DR + 1
    else
      DR := DR - 1;
  end;

  if abs(DC) < 10 then begin
    Str(DR:2:0, DS);
    DS := TrimLeadS(DS);
    if DC < 0 then begin
      if DC > -1 then
        DS := '- 0'
      else
        DS := '- ' + DS[2];
    end else
      DS := '+ ' + DS;
  end else begin
    Str(DR:3:0, DS);
    DS := TrimLeadS(DS);
    if DC < 0 then begin
      Delete(DS,1,1);
      DS := '-' + DS;
    end else
      DS := '+' + DS;
  end;
  if MS[1] = ' ' then
    MS[1] := '0';
  Result := DS + ' ' + MS;
end;

function DateTimeToAJD(D : TDateTime) : Double;
begin
  Result := D + AJDOffset;
end;

function AJDToDateTime(D : Double) : TDateTime;
begin
  Result := D - AJDOffset;
end;


initialization
  AJDOffSet := AstJulianDatePrim(1600, 1, 1, 0) - EncodeDate(1600, 1, 1);
end.