stastro.pas

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

(* ***** BEGIN LICENSE BLOCK *****
 * Version: MPL 1.1
 *
 * The contents of this file are subject to the Mozilla Public License Version
 * 1.1 (the "License"); you may not use this file except in compliance with
 * the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * Software distributed under the License is distributed on an "AS IS" basis,
 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
 * for the specific language governing rights and limitations under the
 * License.
 *
 * The Original Code is TurboPower SysTools
 *
 * The Initial Developer of the Original Code is
 * TurboPower Software
 *
 * Portions created by the Initial Developer are Copyright (C) 1996-2002
 * the Initial Developer. All Rights Reserved.
 *
 * Contributor(s):
 *
 * ***** END LICENSE BLOCK ***** *)

{*********************************************************}
{* SysTools: StAstro.pas 4.03                            *}
{*********************************************************}
{* SysTools: Astronomical Routines                       *}
{*********************************************************}

{$I StDefine.inc}

{ ************************************************************** }
{ Sources:                                                       }
{   1. Astronomical Algorithms, Jean Meeus, Willmann-Bell, 1991. }
{                                                                }
{   2. Planetary and Lunar Coordinates (1984-2000), U.S. Govt,   }
{         1983.                                                  }
{                                                                }
{   3. Supplement to the American Ephemeris and Nautical Almanac,}
{        U.S. Govt, 1964.                                        }
{                                                                }
{   4. MPO96 source files, Brian D. Warner, 1995.                }
{                                                                }
{ ************************************************************** }

unit StAstro;

interface

uses
  Windows, SysUtils,
  StConst, StBase, StDate, StStrS, StDateSt, StMath;

type
  TStTwilight = (ttCivil, ttNautical, ttAstronomical);

  TStRiseSetRec = record
    ORise : TStTime;
    OSet  : TStTime;
  end;

  TStPosRec = record
    RA,
    DC : Double;
  end;
  TStPosRecArray = array[1..3] of TStPosRec;

  TStSunXYZRec = record
    SunX,
    SunY,
    SunZ,
    RV,
    SLong,
    SLat   : Double;
  end;

  TStLunarRecord = record
    T : array[0..1] of TStDateTimeRec;
  end;

  TStPhaseRecord = packed record
    NMDate,
    FQDate,
    FMDate,
    LQDate  : Double;
  end;
  TStPhaseArray = array[0..13] of TStPhaseRecord;

  TStDLSDateRec = record
    Starts : TStDate;
    Ends   : TStDate;
  end;

  TStMoonPosRec = record
    RA,
    DC,
    Phase,
    Dia,
    Plx,
    Elong   : Double;
  end;


const
  radcor  = 57.29577951308232;    {number of degrees in a radian}
  StdDate = 2451545.0;            {Ast. Julian Date for J2000 Epoch}
  OB2000  = 0.409092804;          {J2000 obliquity of the ecliptic (radians)}


{Sun procedures/functions}
function AmountOfSunlight(LD : TStDate; Longitude, Latitude : Double) : TStTime;
  {-compute the hours, min, sec of sunlight on a given date}
function FixedRiseSet(LD : TStDate; RA, DC, Longitude, Latitude : Double) : TStRiseSetRec;
  {-compute the rise and set time for a fixed object, e.g., a star}
function SunPos(UT : TStDateTimeRec) : TStPosRec;
  {-compute the J2000 RA/Declination of the Sun}
function SunPosPrim(UT : TStDateTimeRec) : TStSunXYZRec;
  {-compute the J2000 geocentric rectangular & ecliptic coordinates of the sun}
function SunRiseSet(LD : TStDate; Longitude, Latitude : Double) : TStRiseSetRec;
  {-compute the Sun rise or set time}
function Twilight(LD : TStDate; Longitude, Latitude : Double; TwiType : TStTwilight) : TStRiseSetRec;
  {-compute the beginning and end of twilight (civil, nautical, or astron.)}

{Lunar procedures/functions}
function LunarPhase(UT : TStDateTimeRec) : Double;
  {-compute the phase of the moon}
function MoonPos(UT : TStDateTimeRec) : TStMoonPosRec;
  {-compute the J2000 RA/Declination of the moon}
function MoonRiseSet(LD : TStDate; Longitude, Latitude : Double) : TStRiseSetRec;
  {-compute the Moon rise and set time}
function FirstQuarter(D : TStDate) : TStLunarRecord;
  {-compute date/time of FirstQuarter(s)}
function FullMoon(D : TStDate) : TStLunarRecord;
  {-compute the date/time of FullMoon(s)}
function LastQuarter(D : TStDate) : TStLunarRecord;
  {-compute the date/time of LastQuarter(s)}
function NewMoon(D : TStDate) : TStLunarRecord;
  {-compute the date/time of NewMoon(s)}
function NextFirstQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest FirstQuarter}
function NextFullMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest FullMoon}
function NextLastQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest LastQuarter}
function NextNewMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the next closest NewMoon}
function PrevFirstQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest FirstQuarter}
function PrevFullMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest FullMoon}
function PrevLastQuarter(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest LastQuarter}
function PrevNewMoon(D : TStDate) : TStDateTimeRec;
  {-compute the date/time of the prev closest NewMoon}

{Calendar procedures/functions}
function SiderealTime(UT : TStDateTimeRec) : Double;
  {-compute Sidereal Time at Greenwich}
function Solstice(Y, Epoch : Integer; Summer : Boolean) : TStDateTimeRec;
  {-compute the date/time of the summer or winter solstice}
function Equinox(Y, Epoch : Integer; Vernal : Boolean) : TStDateTimeRec;
  {-compute the date/time of the vernal/autumnal equinox}
function Easter(Y, Epoch : Integer) : TStDate;
  {-compute the date of Easter (astronmical calendar)}

{Astronomical Julian Date Conversions}
function DateTimeToAJD(D : TDateTime) : Double;

{Conversion routines}
function HoursMin(RA : Double) : ShortString;
  {-convert RA to hh:mm:ss string}
function DegsMin(DC : Double) : ShortString;
  {-convert Declination to +/-dd:mm:ss string}
function AJDToDateTime(D : Double) : TDateTime;


implementation

var
  AJDOffset : Double;

function CheckDate(UT : TStDateTimeRec) : Boolean;
begin
  with UT do begin
    if (D < MinDate) or (D > MaxDate) or
       (T < 0) or (T > MaxTime) then
      Result := False
    else
      Result := True;
  end;
end;

function CheckYear(Y, Epoch : Integer) : Integer;
begin
  if Y < 100 then begin
    if Y >= (Epoch mod 100) then
      Result := ((Epoch div 100) * 100) + Y
    else
      Result := ((Epoch div 100) * 100) + 100 + Y;
  end else
    Result := Y;
end;

function SiderealTime(UT : TStDateTimeRec) : Double;
  {-compute Sidereal Time at Greenwich in degrees}
var
  T,
  JD : Double;
begin
  if not CheckDate(UT) then begin
    Result := -1;
    Exit;
  end;

  JD := AstJulianDate(UT.D) + UT.T/86400;

  T := (JD - 2451545.0) / 36525.0;

  Result := 280.46061837
          + 360.98564736629 * (JD - 2451545.0)
          + 0.000387933 * sqr(T)
          - (sqr(T) * T / 38710000);
  Result := Frac(Result/360.0) * 360.0;
  if Result < 0 then
    Result := 360 + Result;
end;

function SunPosPrim(UT : TStDateTimeRec) : TStSunXYZRec;
  {-compute J2000 XYZ coordinates of the Sun}
var
  JD,
  T0,
  A,
  L,
  B,
  X,Y,Z  : Double;

begin
  if not CheckDate(UT) then begin
    with Result do begin
      SunX := -99;
      SunY := -99;
      SunZ := -99;
      RV := -99;
      SLong := -99;
      SLat := -99;
    end;
    Exit;
  end;

  JD := AstJulianDate(UT.D) + UT.T/86400;
  T0 := (JD - StdDate) / 365250;

{solar longitude}
  L := 175347046
     +   3341656 * cos(4.6692568 +  6283.07585*T0)
     +     34894 * cos(4.6261000 + 12566.1517*T0)
     +      3497 * cos(2.7441000 +  5753.3849*T0)
     +      3418 * cos(2.8289000 +     3.5231*T0)
     +      3136 * cos(3.6277000 + 77713.7715*T0)
     +      2676 * cos(4.4181000 +  7860.4194*T0)
     +      2343 * cos(6.1352000 +  3930.2097*T0)
     +      1324 * cos(0.7425000 + 11506.7698*T0)
     +      1273 * cos(2.0371000 +   529.6910*T0)
     +      1199 * cos(1.1096000 +  1577.3435*T0)
     +       990 * cos(5.2330000 +  5884.9270*T0)
     +       902 * cos(2.0450000 +    26.1490*T0)
     +       857 * cos(3.5080000 +   398.149*T0)
     +       780 * cos(1.1790000 +  5223.694*T0)
     +       753 * cos(2.5330000 +  5507.553*T0)
     +       505 * cos(4.5830000 + 18849.228*T0)
     +       492 * cos(4.2050000 +   775.523*T0)
     +       357 * cos(2.9200000 +     0.067*T0)
     +       317 * cos(5.8490000 + 11790.626*T0)
     +       284 * cos(1.8990000 +   796.298*T0)
     +       271 * cos(0.3150000 + 10977.079*T0)
     +       243 * cos(0.3450000 +  5486.778*T0)
     +       206 * cos(4.8060000 +  2544.314*T0)
     +       205 * cos(1.8690000 +  5573.143*T0)
     +       202 * cos(2.4580000 +  6069.777*T0)
     +       156 * cos(0.8330000 +   213.299*T0)
     +       132 * cos(3.4110000 +  2942.463*T0)
     +       126 * cos(1.0830000 +    20.775*T0)
     +       115 * cos(0.6450000 +     0.980*T0)
     +       103 * cos(0.6360000 +  4694.003*T0)
     +       102 * cos(0.9760000 + 15720.839*T0)
     +       102 * cos(4.2670000 +     7.114*T0)
     +        99 * cos(6.2100000 +  2146.170*T0)
     +        98 * cos(0.6800000 +   155.420*T0)
     +        86 * cos(5.9800000 +161000.690*T0)
     +        85 * cos(1.3000000 +  6275.960*T0)
     +        85 * cos(3.6700000 + 71430.700*T0)
     +        80 * cos(1.8100000 + 17260.150*T0);

  A := 628307584999.0
     + 206059 * cos(2.678235 +  6283.07585*T0)
     +   4303 * cos(2.635100 + 12566.1517*T0)
     +    425 * cos(1.590000 +     3.523*T0)
     +    119 * cos(5.796000 +    26.298*T0)
     +    109 * cos(2.966000 +  1577.344*T0)
     +     93 * cos(2.590000 + 18849.23*T0)
     +     72 * cos(1.140000 +   529.69*T0)
     +     68 * cos(1.870000 +   398.15*T0)
     +     67 * cos(4.410000 +  5507.55*T0)
     +     59 * cos(2.890000 +  5223.69*T0)
     +     56 * cos(2.170000 +   155.42*T0)
     +     45 * cos(0.400000 +   796.30*T0)
     +     36 * cos(0.470000 +   775.52*T0)
     +     29 * cos(2.650000 +     7.11*T0)
     +     21 * cos(5.340000 +     0.98*T0)
     +     19 * cos(1.850000 +  5486.78*T0)
     +     19 * cos(4.970000 +   213.30*T0)
     +     17 * cos(2.990000 +  6275.96*T0)
     +     16 * cos(0.030000 +  2544.31*T0);
  L := L + (A * T0);

  A := 8722 * cos(1.0725 +  6283.0758*T0)
     +  991 * cos(3.1416)
     +  295 * cos(0.437  + 12566.1520*T0)
     +   27 * cos(0.050  +     3.52*T0)
     +   16 * cos(5.190  +    26.30*T0)
     +   16 * cos(3.69   +   155.42*T0)
     +    9 * cos(0.30   + 18849.23*T0)
     +    9 * cos(2.06   + 77713.77*T0);
  L := L + (A * sqr(T0));

  A := 289 * cos(5.842 +  6283.076*T0)
     +  21 * cos(6.05  + 12566.15*T0)
     +   3 * cos(5.20  +   155.42*T0)
     +   3 * cos(3.14);
  L := L + (A * sqr(T0) * T0);
  L := L / 1.0E+8;


{solar latitude}
  B := 280 * cos(3.199 + 84334.662*T0)
     + 102 * cos(5.422 +  5507.553*T0)
     +  80 * cos(3.88  +  5223.69*T0)
     +  44 * cos(3.70  +  2352.87*T0)
     +  32 * cos(4.00  +  1577.34*T0);
  B := B / 1.0E+8;

  A := 227778 * cos(3.413766 + 6283.07585*T0)
     +   3806 * cos(3.3706 + 12566.1517*T0)
     +   3620
     +     72 * cos(3.33 + 18849.23*T0)
     +      8 * cos(3.89 +  5507.55*T0)
     +      8 * cos(1.79 +  5223.69*T0)
     +      6 * cos(5.20 +  2352.87*T0);
  B := B + (A * T0 / 1.0E+8);

  A := 9721 * cos(5.1519 + 6283.07585*T0)
     +  233 * cos(3.1416)
     +  134 * cos(0.644  + 12566.152*T0)
     +    7 * cos(1.07   + 18849.23*T0);
  B := B + (A * sqr(T0) / 1.0E+8);

  A := 276 * cos(0.595 + 6283.076*T0)
     +  17 * cos(3.14)
     +   4 * cos(0.12  + 12566.15*T0);
  B := B + (A * sqr(T0) * T0 / 1.0E+8);


{solar radius vector (astronomical units)}
  Result.RV := 100013989
             +   1670700 * cos(3.0984635 +  6283.07585*T0)
             +     13956 * cos(3.05525   + 12566.15170*T0)
             +      3084 * cos(5.1985    + 77713.7715*T0)
             +      1628 * cos(1.1739    +  5753.3849*T0)
             +      1576 * cos(2.8649    +  7860.4194*T0)
             +       925 * cos(5.453     + 11506.770*T0)
             +       542 * cos(4.564     +  3930.210*T0)
             +       472 * cos(3.661     +  5884.927*T0)
             +       346 * cos(0.964     +  5507.553*T0)
             +       329 * cos(5.900     +  5223.694*T0)
             +       307 * cos(0.299     +  5573.143*T0)
             +       243 * cos(4.273     + 11790.629*T0)
             +       212 * cos(5.847     +  1577.344*T0)
             +       186 * cos(5.022     + 10977.079*T0)
             +       175 * cos(3.012     + 18849.228*T0)
             +       110 * cos(5.055     +  5486.778*T0)
             +        98 * cos(0.89      +  6069.78*T0)
             +        86 * cos(5.69      + 15720.84*T0)
             +        86 * cos(1.27      +161000.69*T0)
             +        65 * cos(0.27      + 17260.15*T0)
             +        63 * cos(0.92      +   529.69*T0)
             +        57 * cos(2.01      + 83996.85*T0)
             +        56 * cos(5.24      + 71430.70*T0)
             +        49 * cos(3.25      +  2544.31*T0)
             +        47 * cos(2.58      +   775.52*T0)
             +        45 * cos(5.54      +  9437.76*T0)
             +        43 * cos(6.01      +  6275.96*T0)
             +        39 * cos(5.36      +  4694.00*T0)
             +        38 * cos(2.39      +  8827.39*T0)
             +        37 * cos(0.83      + 19651.05*T0)
             +        37 * cos(4.90      + 12139.55*T0)
             +        36 * cos(1.67      + 12036.46*T0)
             +        35 * cos(1.84      +  2942.46*T0)
             +        33 * cos(0.24      +  7084.90*T0)
             +        32 * cos(0.18      +  5088.63*T0)
             +        32 * cos(1.78      +   398.15*T0)
             +        28 * cos(1.21      +  6286.60*T0)
             +        28 * cos(1.90      +  6279.55*T0)
             +        26 * cos(4.59      + 10447.39*T0);
  Result.RV := Result.RV / 1.0E+8;

  A := 103019 * cos(1.107490 +  6283.075850*T0)
     +   1721 * cos(1.0644   + 12566.1517*T0)
     +    702 * cos(3.142)
     +     32 * cos(1.02     + 18849.23*T0)
     +     31 * cos(2.84     +  5507.55*T0)
     +     25 * cos(1.32     +  5223.69*T0)
     +     18 * cos(1.42     +  1577.34*T0)
     +     10 * cos(5.91     + 10977.08*T0)
     +      9 * cos(1.42     +  6275.96*T0)
     +      9 * cos(0.27     +  5486.78*T0);
  Result.RV := Result.RV + (A * T0 / 1.0E+8);

  A := 4359 * cos(5.7846 +  6283.0758*T0)
     +  124 * cos(5.579  + 12566.152*T0)
     +   12 * cos(3.14)
     +    9 * cos(3.63   + 77713.77*T0)
     +    6 * cos(1.87   +  5573.14*T0)
     +    3 * cos(5.47   + 18849.23*T0);
  Result.RV := Result.RV + (A * sqr(T0) / 1.0E+8);

  L := (L + PI);
  L := Frac(L / 2.0 / PI) * 2.0 * Pi;
  if L < 0 then
    L := L + (2.0*PI);
  B := -B;

  Result.SLong := L * radcor;
  Result.SLat  := B * radcor;

  X := Result.RV * cos(B) * cos(L);
  Y := Result.RV * cos(B) * sin(L);
  Z := Result.RV * sin(B);

  Result.SunX :=               X +     4.40360E-7 * Y - 1.90919E-7 * Z;
  Result.SunY := -4.79966E-7 * X + 0.917482137087 * Y - 0.397776982902 * Z;
  Result.SunZ :=                   0.397776982902 * Y + 0.917482137087 * Z;
end;

function MoonPosPrim(UT : TStDateTimeRec) : TStMoonPosRec;
  {-computes J2000 coordinates of the moon}