function GetSunTime( const Dt: TDate; const Latitude, Longitude, Zenith: Extended;
  const LocalOffset: Byte; const SunTime: TSunTime ): TTime;
var
  N: Word;
  LngHour, M, L, RA, Lquadrant, RAquadrant,
  sinDec, cosDec, HCos, H, LocalT, UT: Double;
  t: TTime;
begin
   // 1. first calculate the day of the year
   N := Trunc( GetDays( Dt ) );

// ----------

   // 2. convert the longitude to hour value and calculate an approximate time
   LngHour := Longitude / 15;
   // if rising time is desired:
	 // t := N + ( ( 6 - LngHour ) / 24 );
	 // if setting time is desired:
	 // t := N + ( ( 18 - LngHour ) / 24 );
   case SunTime of
      Sunrise: t := N + ( ( 6 - LngHour ) / 24 );
      Sunset: t := N + ( ( 18 - LngHour ) / 24 );
   end;

// ----------

   // 3. calculate the Sun's mean anomaly
   M := ( 0.9856 * t ) - 3.289;

// ----------

   // 4. calculate the Sun's true longitude
   //    [Note throughout the arguments of the trig functions
   //    (sin, tan) are in degrees. It will likely be necessary to
   //    convert to radians. eg sin(170.626 deg) =sin(170.626*pi/180
   //    radians)=0.16287]
   L := M + ( 1.916 * Sin( DegToRad( M ) ) ) +
        ( 0.020 * Sin( DegToRad( 2 * M ) ) ) + 282.634;
   // NOTE: L potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
   Adjust( L, 360 );

// ----------

   // 5a. calculate the Sun's right ascension
   RA := RadToDeg( ArcTan( 0.91764 * Tan( DegToRad( L ) ) ) );
	 // NOTE: RA potentially needs to be adjusted into the range [0,360) by adding/subtracting 360
   Adjust( RA, 360 );

// ----------

   // 5b. right ascension value needs to be in the same quadrant as L
   Lquadrant := Floor( L / 90 ) * 90;
   RAquadrant := Floor( RA / 90 ) * 90;
   RA := RA + ( Lquadrant - RAquadrant );

// ----------

   // 5c. right ascension value needs to be converted into hours
   RA := RA / 15;

// ----------

   // 6. calculate the Sun's declination
   sinDec := 0.39782 * Sin( DegToRad( L ) );
   cosDec := Cos( ArcSin( sinDec ) );

// ----------

   // 7a. calculate the Sun's local hour angle
   HCos := ( Cos( DegToRad( Zenith ) ) - ( sinDec * Sin( DegToRad( Latitude ) ) ) ) /
           ( cosDec * Cos( DegToRad( Latitude ) ) );
   // if (cosH >  1)
	 //    the sun never rises on this location (on the specified date)
	 // if (cosH < -1)
	 //    the sun never sets on this location (on the specified date)
   if ( HCos > 1 ) or ( HCos < -1 ) then
      Exit;

// ----------

   // 7b. finish calculating H and convert into hours
   // if  rising time is desired:
   // H := 360 - RadToDeg( ArcCos( HCos ) );
   // if setting time is desired:
   // H := RadToDeg( ArcCos( HCos ) );
   case SunTime of
      Sunrise: H := 360 - RadToDeg( ArcCos( HCos ) );
      Sunset: H := RadToDeg( ArcCos( HCos ) );
   end;

   H := H / 15;

// ----------

   // 8. calculate local mean time of rising/setting
   LocalT := H + RA - ( 0.06571 * t ) - 6.622;

// ----------

   // 9. adjust back to UTC

   UT := LocalT - LngHour;
   // NOTE: UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24
   Adjust( UT, 24 );

// ----------

   // 10. convert UT value to local time zone of latitude/longitude
   Result := UT + LocalOffset;
   Adjust( Double( Result ), 24 );
end;

ST = UTC + долгота в часовых единицах + EOT (значение уравнения времени)