steclpse.pas

来自「条码控件： 一维条码控件 二维条码控件 PDF417Barcode Ma」· PAS 代码 · 共 732 行 · 第 1/2 页

 non-central   Axis of moon's umbral shadow does not touch earth's surface
               Eclipse is usually partial but can be one of possibilities
               for central eclipse if very near one of the earth's poles

 Partial eclipses
 ----------------
 partial       None of the moon's umbral shadow touches the earth's surface
}

    if (abs(Gamma) <= (0.9972 + abs(Mu))) then
      NonPartialSolarEclipse(CentralJD, Mu, Gamma)
    else begin
      if (abs(Gamma) < 1.5433 + Mu) then
        PartialSolarEclipse(CentralJD, Mu, Gamma);
    end;
  end;
end;

procedure TStEclipses.TotalLunarEclipse(CentralJD, MoonAnom, Mu,
                                        PMag, UMag, Gamma : Double);
var
  TLE              : PStEclipseRecord;
  PartialSemiDur,
  TotalSemiDur,
  Dbl1             : Double;
begin
  New(TLE);
  TLE^.Magnitude  := UMag;
  TLE^.Hemisphere := htNone;
  TLE^.EType      := etLunarTotal;
  TLE^.Path       := nil;
  CentralJD := AJDToDateTime(CentralJD);

  PartialSemiDur := 1.0128 - Mu;
  TotalSemiDur   := 0.4678 - Mu;
  Dbl1 := 0.5458 + 0.0400 * cos(MoonAnom);

  PartialSemiDur := 60/Dbl1 * sqrt(sqr(PartialSemiDur) - sqr(Gamma)) / 1440;
  TotalSemiDur   := 60/Dbl1 * sqrt(sqr(TotalSemiDur) - sqr(Gamma)) / 1440;

  TLE^.LContacts.FirstContact  := CentralJD - PartialSemiDur;
  TLE^.LContacts.SecondContact := CentralJD - TotalSemiDur;
  TLE^.LContacts.MidEclipse    := CentralJD;
  TLE^.LContacts.ThirdContact  := CentralJD + TotalSemiDur;
  TLE^.LContacts.FourthContact := CentralJD + PartialSemiDur;

  PartialSemiDur := 60/Dbl1 * sqrt(sqr(1.5573 + Mu) - sqr(Gamma)) / 1440;
  TLE^.LContacts.UT1 := CentralJD - PartialSemiDur;
  TLE^.LContacts.UT2 := CentralJD + PartialSemiDur;

  Self.Append(TLE);
end;

procedure TStEclipses.PartialLunarEclipse(CentralJD, MoonAnom, Mu,
                                          PMag, UMag, Gamma : Double);
var
  TLE              : PStEclipseRecord;
  PartialSemiDur,
  Dbl1             : Double;
begin
  New(TLE);
  TLE^.Magnitude := UMag;
  TLE^.Hemisphere := htNone;
  TLE^.EType      := etLunarPartial;
  TLE^.Path       := nil;
  CentralJD := AJDToDateTime(CentralJD);

  PartialSemiDur := 1.0128 - Mu;
  Dbl1 := 0.5458 + 0.0400 * cos(MoonAnom);

  PartialSemiDur := 60/Dbl1 * sqrt(sqr(PartialSemiDur) - sqr(Gamma)) / 1440;

  TLE^.LContacts.FirstContact  := CentralJD - PartialSemiDur;
  TLE^.LContacts.SecondContact := 0;
  TLE^.LContacts.MidEclipse    := CentralJD;
  TLE^.LContacts.ThirdContact  := 0;
  TLE^.LContacts.FourthContact := CentralJD + PartialSemiDur;

  PartialSemiDur := 60/Dbl1 * sqrt(sqr(1.5573 + Mu) - sqr(Gamma)) / 1440;
  TLE^.LContacts.UT1 := CentralJD - PartialSemiDur;
  TLE^.LContacts.UT2 := CentralJD + PartialSemiDur;

  Self.Append(TLE);
end;

procedure TStEclipses.PenumbralLunarEclipse(CentralJD, MoonAnom, Mu,
                                            PMag, UMag, Gamma : Double);
var
  TLE              : PStEclipseRecord;
  PartialSemiDur,
  Dbl1             : Double;
begin
  New(TLE);
  TLE^.Magnitude := PMag;
  TLE^.Hemisphere := htNone;
  TLE^.EType      := etLunarPenumbral;
  TLE^.Path       := nil;
  CentralJD := AJDToDateTime(CentralJD);

  TLE^.LContacts.FirstContact  := 0;
  TLE^.LContacts.SecondContact := 0;
  TLE^.LContacts.MidEclipse    := CentralJD;
  TLE^.LContacts.ThirdContact  := 0;
  TLE^.LContacts.FourthContact := 0;

  Dbl1 := 0.5458 + 0.0400 * cos(MoonAnom);
  PartialSemiDur := 60/Dbl1 * sqrt(sqr(1.5573 + Mu) - sqr(Gamma)) / 1440;
  TLE^.LContacts.UT1 := CentralJD - PartialSemiDur;
  TLE^.LContacts.UT2 := CentralJD + PartialSemiDur;

  Self.Append(TLE);
end;

procedure TStEclipses.GetBesselianElements(CentralJD : Double);
var
  I,
  Mins       : LongInt;
  CurJD,
  SidTime,
  SunDist,
  MoonDist,
  DistRatio,
  Alpha,
  Theta,
  Sun1,
  Sun2,
  Moon1,
  Moon2,
  Dbl3,
  F1, F2     : Double;
  DTRec      : TStDateTimeRec;
  SunXYZ     : TStSunXYZRec;
  Sun        : TStPosRec;
  Moon       : TStMoonPosRec;
begin
{compute BesselianElements every 10 minutes starting 2 hours prior to CentralJD}
{but forcing positions to be at multiple of ten minutes}
  for I := 1 to 25 do begin
    CurJD   := CentralJD + ((I-13) * (10/1440));
    DTRec.D := AstJulianDateToStDate(CurJD, True);
    if ((Frac(CurJD) + 0.5) >= 1) then
      Mins := Trunc(((Frac(CurJD) + 0.5)-1) * 1440)
    else
      Mins    := Trunc((Frac(CurJD) + 0.5) * 1440);
   {changed because, for example, both 25 and 35 rounded to 30}
    Mins    := ((Mins + 5) div 10) * 10;
    if (Mins = 1440) then
      Mins := 0;
    DTRec.T := Mins * 60;

    SidTime := SiderealTime(DTRec) / radcor;
    SunXYZ  := SunPosPrim(DTRec);
    Sun     := SunPos(DTRec);
    Moon    := MoonPos(DTRec);

    Sun1   := Sun.RA/radcor;
    Sun2   := Sun.DC/radcor;
    Moon1  := Moon.RA/radcor;
    Moon2  := Moon.DC/radcor;

    FBesselianElements[I].JD := StDateToDateTime(DTRec.D)
                              + StTimeToDateTime(DTRec.T);

    SunDist   := 1.0 / sin(8.794/SunXYZ.RV/3600.0/radcor);
    MoonDist  := 1.0 / sin(Moon.Plx/radcor);
    DistRatio := MoonDist / SunDist;

    Theta := DistRatio/cos(Sun2) * cos(Moon2) * (Moon1 - Sun1);
    Theta := Theta/(1.0-DistRatio);
    Alpha := Sun1 - Theta;

    Theta := DistRatio/(1.0 - DistRatio) * (Moon2 - Sun2);

    FBesselianElements[I].Delta := Sun2 - Theta;
    FBesselianElements[I].Angle := SidTime - Alpha;
    FBesselianElements[I].XAxis := MoonDist * cos(Moon2) * (sin(Moon1 - Alpha));

    Dbl3 := FBesselianElements[I].Delta;
    FBesselianElements[I].YAxis := MoonDist * (sin(Moon2) * cos(Dbl3)
                                 - cos(Moon2) * sin(Dbl3) * cos(Moon1 - Alpha));

    Theta := DistRatio * SunXYZ.RV;
    Theta := SunXYZ.RV - Theta;
    F1 := StInvSin(0.004664012/Theta);
    F2 := StInvSin(0.004640787/Theta);

    Theta := MoonDist * (sin(Moon2) * sin(Dbl3) + cos(Moon2)
             * cos(Dbl3) * cos(Moon1 - Alpha));
    FBesselianElements[I].L1 := (Theta + 0.272453/sin(F1)) * StTan(F1);
    FBesselianElements[I].L2 := (Theta - 0.272453/sin(F2)) * StTan(F2);

    if (I = 13) then
      F0 := StTan(F2);

    if (I = 16) then begin
      FUPrime := FBesselianElements[16].Angle - FBesselianElements[10].Angle;
      FDPrime := FBesselianElements[16].Delta - FBesselianElements[10].Delta;
    end;
  end;
end;

procedure TStEclipses.GetShadowPath(I1, I2 : Integer; Path : TStList);
var
  J,
  I3, I4,
  I5, I6         : integer;

  Delta,
  Dbl1,
  Dbl2,
  P1,
  XAxis,
  YAxis,
  Eta,
  R1, R2,
  D1, D2,
  D3, D4,
  V3, V4,
  V5, V6, V7,
  XPrime,
  YPrime      : Double;

  Position    : PStLongLat;
begin
  for J := I1 to I2 do begin
    Eta := 0.00669454;
    Delta := FBesselianElements[J].Delta;
    XAxis := FBesselianElements[J].XAxis;
    YAxis := FBesselianElements[J].YAxis;

    R1 := sqrt(1.0 - Eta * sqr(cos(Delta)));
    R2 := sqrt(1.0 - Eta * sqr(sin(Delta)));

    D1 := sin(Delta) / R1;
    D2 := sqrt(1.0 - Eta) * cos(Delta) / R1;
    D3 := Eta * sin(Delta) * cos(Delta) / R1 / R2;
    D4 := sqrt(1.0 - Eta) / R1 / R2;

    V3 := YAxis / R1;
    V4 := sqrt(1.0 - sqr(XAxis) - sqr(V3));
    V5 := R2 * (V4 * D4 - V3 * D3);
    V6 := FUPrime * (-YAxis * sin(Delta) + V5 * cos(Delta));
    V7 := FUPrime * XAxis * sin(Delta) - FDPrime * V5;

    if ((I2-I1) div 2) >= 4 then begin
      I3 := (I2-I1) div 2;
      I4 := I1 + I3;
      I5 := I4 - 3;
      I6 := I4 + 3;
      XPrime := FBesselianElements[I6].XAxis
              - FBesselianElements[I5].XAxis;
      YPrime := FBesselianElements[I6].YAxis
              - FBesselianElements[I5].YAxis;
    end else begin
      XPrime := (FBesselianElements[J+1].XAxis -
                 FBesselianElements[J-1].XAxis) * 3;
      YPrime := (FBesselianElements[J+1].YAxis -
                 FBesselianElements[J-1].YAxis) * 3;
    end;

    New(Position);
    Dbl1 := sqrt(sqr(XPrime - V6) + sqr(YPrime - V7));
    Position^.JD := FBesselianElements[J].JD;

    Dbl2 := (FBesselianElements[J].L2 - V5 * F0) / Dbl1;
    Dbl2 := abs(Dbl2 * 120);
    Position^.Duration := int(Dbl2) + frac(Dbl2) * 0.6;

    Dbl1 := -V3 * D1 + V4 * D2;
    P1 := StInvTan2(Dbl1, XAxis);

    Dbl2 := (FBesselianElements[J].Angle - P1) * radcor;
    Dbl2 := frac(Dbl2 / 360.0) * 360;
    if (Dbl2 < 0) then
      Dbl2 := Dbl2 + 360.0;
    if (Dbl2 > 180.0) and (Dbl2 < 360.0) then
       Dbl2 := Dbl2 - 360.0;
    Position^.Longitude := Dbl2;

    Dbl1 := StInvSin(V3 * D2 + V4 * D1);
    Dbl2 := ArcTan(1.003364 * StTan(Dbl1)) * radcor;
    Position^.Latitude := Dbl2;

    Path.Append(Position);
  end;
end;

procedure TStEclipses.NonPartialSolarEclipse(CentralJD, Mu, Gamma : Double);
var
  SolEc   : PStEclipseRecord;
  I1, I2  : Integer;
begin
  New(SolEc);
  if (Mu < 0) then
    SolEc^.EType := etSolarTotal
  else if (Mu > 0.0047) then
    SolEc^.EType := etSolarAnnular
  else begin
    if (Mu < (0.00464 * sqrt(1 - sqr(Gamma)))) then
      SolEc^.EType := etSolarAnnularTotal
    else
      SolEc^.EType := etSolarTotal;
  end;

  SolEc^.Magnitude := -1;
  if (Gamma > 0) then
    SolEc^.Hemisphere := htNorthern
  else
    SolEc^.Hemisphere := htSouthern;
  FillChar(SolEc^.LContacts, SizeOf(TStContactTimes), #0);
  SolEc^.LContacts.MidEclipse := AJDtoDateTime(CentralJD);

  GetBesselianElements(CentralJD);

{find limits - then go one step inside at each end}
  I1 := 1;
  while (sqr(FBesselianElements[I1].XAxis) +
         sqr(FBesselianElements[I1].YAxis) >= 1.0) and (I1 < 25) do
    Inc(I1);
  Inc(I1);

  I2 := I1;
  repeat
    if (I2 < 25) then begin
      if (sqr(FBesselianElements[I2+1].XAxis) +
          sqr(FBesselianElements[I2+1].YAxis) < 1) then
        Inc(I2)
      else
        break;
    end;
  until (sqr(FBesselianElements[I2].XAxis) +
        sqr(FBesselianElements[I2].YAxis) >= 1) or (I2 >= 25);
  Dec(I2);

{this test accounts for non-central eclipses, i.e., those that are total}
{and/or annular but the axis of the moon's shadow does not touch the earth}
  if (I1 <> I2) and (I1 < 26) and (I2 < 26) then begin
    SolEc^.Path := TStList.Create(TStListNode);
    SolEc^.Path.DisposeData := DisposePathData;
    GetShadowPath(I1, I2, SolEc^.Path);
  end else
    SolEc^.Path := nil;
  Self.Append(SolEc);
end;

procedure TStEclipses.PartialSolarEclipse(CentralJD, Mu, Gamma : Double);
var
  SolEc   : PStEclipseRecord;
begin
  New(SolEc);
  SolEc^.EType := etSolarPartial;
  SolEc^.Magnitude := (1.5433 + Mu - abs(Gamma)) / (0.5461 + 2*Mu);
  if (Gamma > 0) then
    SolEc^.Hemisphere := htNorthern
  else
    SolEc^.Hemisphere := htSouthern;
  FillChar(SolEc^.LContacts, SizeOf(TStContactTimes), #0);
  SolEc^.LContacts.MidEclipse := AJDToDateTime(CentralJD);
  SolEc^.Path := nil;
  Self.Append(SolEc);
end;


end.