Sunrise/Sunset

I want to thank all those you provided information on my request for
sunrise and sunset calculations.  While the component and packages
looked quite nice, I felt it was a bit overboard.  I did find some Java
code at the USGS site that I was able to convert.  This is pretty
accurate to within about 30 seconds for latitude below 70 degrees and 10
minutes for latitudes above 70;

Here is the unit for those that would like to add this to their arsenal.

Jeff Steinkamp


----------------------------------astor unit--------------------

unit astro;


interface
uses math,sysutils,dateutils;

type tSun = record
     srise : Tdatetime;
     sset : Tdatetime;
     snoon : tdatetime;

end;

Function calcJD(ayear,amonth,aday : integer) : extended;
Function SolarStuff(adate : TdateTime; alat,alon : extended) : Tsun;
Function calcJDfromJulianCent(t : extended): Extended;
Function calcTimeJulianCent(jd : extended) : extended;
Function calcSunSetUTC(jd,latitude,longitude : extended) :extended;
function calcSolarNoonUTC(t,longitude : Extended) : Extended;
Function calcSunRiseUTC(jd,latitude,longitude : extended) :extended;
Function calcHourAngleSunrise(lat,solarDec : extended) : extended;
function calcHourAngleSunset(lat,solarDec : extended) : extended;
Function calcEquationOfTime(t:extended) : Extended;
Function calcObliquityCorrection(t:extended) : extended;
Function calcMeanObliquityOfEcliptic(t: Extended) : extended;
Function calcGeoMeanLongSun(t:extended) : extended;
Function calcEccentricityEarthOrbit(t : extended) : extended;
Function calcJulianDaytoCent(jd : integer) : extended;
Function CalcSunDeclination(t:extended): Extended;
function calcSunAapparentLong(t:extended) : extended;
function calcSunEqOfCenter(t:extended) : extended;
function calcGeoMeanAnomalySun(t : extended) : extended;
function calcSunTrueLong(t:extended) : extended;

implementation

Function calcEquationOfTime(t:extended) : Extended;
//***************************************************************
// Purpose:  calculate the difference between the ture solar time
//           and mean solar time
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     Equation of time in minutes of time
//****************************************************************
var epsilon,l0,e,m,y,Etime: extended;
    Sin2l0,sinm,cos2l0,sin4l0,sin2m : extended;
Begin
    epsilon := calcObliquityCorrection(t);
    l0 := calcGeoMeanLongSun(t);
    e := calcEccentricityEarthOrbit(t);
    m := calcGeoMeanAnomalySun(t);
    y := tan(degtorad(epsilon)/2.0);
    y := y*y;
    sin2l0 := sin(2.0 *degtorad(l0));
    sinm := sin(degtorad(m));
    cos2l0 := cos(2 * degtorad(l0));
    sin4l0 := sin(4 * degtorad(l0));
    sin2m := sin(2 * degtorad(m));
    Etime := y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * -0.5 * y
*y * sin4l0 - 1.25 * e * e * sin2m;
    result := radtodeg(Etime) *4.0;  // in minutes of time

end;

Function calcObliquityCorrection(t:extended) : extended;
//***************************************************************
// Purpose:  calculate the correcrted obliquity of the ecliptic
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     corrected obliquity in degrees
//****************************************************************
var e0,omega : extended;
begin
    e0 := calcMeanObliquityOfEcliptic(t);
    omega := 125.04 - 1934.136 * t;
    result := e0 + 0.00256 * cos(degtorad(omega));  //in degres
end;

Function  calcMeanObliquityOfEcliptic(t: Extended) : extended;
//***************************************************************
// Purpose:  calculate the mean obliquity of the ecliptic
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     mean obliquity in degrees
//****************************************************************
var sec : extended;
begin
    sec := 21.448 - t*(46.82150 + t*(0.00059-t*(0.001813)));
    result := 23.0 + (26.0 - sec/60.0)/60.0;  //in degrees
end;

Function calcGeoMeanLongSun(t:extended) : extended;
//***************************************************************
// Purpose:  calculate the Geometric Mean Longitude of the Sun
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     Geometric Mean Longitude of the sun in degrees
//****************************************************************
var l : extended;
begin
    l := 280.46646 + t *(36000.76983 + 0.0003032 *t);
    if l > 360.0 then l := l - 360;
    if l < 0 then l := l + 360;
    result := l; //in degrees;
end;

Function calcEccentricityEarthOrbit(t : extended) : extended;
//***************************************************************
// Purpose:  calculate the eccentricity of earth's orbig
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     unitless eccentricity
//****************************************************************
var e : extended;
begin
    e := 0.016708634 - t * (0.000042037 + 0.0000002367 * t);
    result := e; //unitless
end;

function calcGeoMeanAnomalySun(t : extended) : extended;
//***************************************************************
// Purpose:  calculate the Geometric Mean Anomoly of the Sun
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     Geometric Mean Anomoly of the sun in degrees
//****************************************************************
var M : extended;
begin
    M := 357.52911 + t * (35999.05029 - 0.0001537 * t);
    result := M; //in degrees
end;

Function calcJulianDaytoCent(jd : integer) : extended;
//***************************************************************
// Purpose:  convert the Julian dat to centuries since J2000.0
//
//Arguments:  jd: Julian day of the year;
//
//Return:     number of Julian centuries since J2000.0
//****************************************************************
var t: extended;
begin
    t := (jd-2451545.0/36525.0);
    result := t;
end;



Function CalcSunDeclination(t:extended): Extended;
//***************************************************************
// Purpose:  calculate the declination of the sun
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     suns's declination in degrees
//****************************************************************
var e,lambda,sint,theta : extended;
begin
    e := calcObliquityCorrection(t);
    lambda := calcSunAapparentLong(t);
    sint := sin(degtorad(e)) * sin(degtorad(lambda));
    theta := radtodeg(arcsin(sint));
    result := theta; //in degrees
end;

function calcSunAapparentLong(t:extended) : extended;
//***************************************************************
// Purpose:  calculate the apparent longitude of the sun
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     suns's apparent longitude in degrees
//****************************************************************
var alpha,omega,lambda : extended;
begin
    alpha := calcSunTrueLong(t);
    omega := 125.05 - 1934.136 * t;
    lambda := alpha - 0.00569 - 0.00478 * sin(degtorad(omega));
    result := lambda; // in degrees
end;

function calcSunTrueLong(t:extended) : extended;
//***************************************************************
// Purpose:  calculate the true longitude of the sun
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     suns's true longitude in degrees
//****************************************************************
var l0,c0 : extended;
begin
    l0 := calcGeoMeanLongSun(t);
    c0 := calcSunEqOfCenter(t);
    result := l0 + c0; //in degrees
end;

function calcSunEqOfCenter(t:extended) : extended;
//***************************************************************
// Purpose:  calculate the equation of the center of the sun
//
//Arguments:  t: number of Julian centuries since J2000.0
//
//Return:     suns's equation of the center in degrees
//****************************************************************
var m,mrad,sinm,sin2m,sin3m,c : extended;
begin
    m := calcGeoMeanAnomalySun(t);
    mrad := degtorad(m);
    sinm := sin(mrad);
    sin2m := sin(2*mrad);
    sin3m := sin(3*mrad);
    c := sinm *(1.914602 - t *(0.004817 + 0.000014 *t)) +
         sin2m * (0.019993 - 0.000101 *t) + sin3m * 0.000289;
         result := c; //in degrees
end;

Function calcHourAngleSunrise(lat,solarDec : extended) : extended;
//***************************************************************
// Purpose:  calculate the hour angle of the sun at sunrise for the
latitude
//
//Arguments:  lat : latitude in degrees
//            solarDec : declination angle of sun in degreees
//
//Return:     hour angle of sunrise in radians
//****************************************************************
var latrad,sdrad,haarg,ha : extended;
begin
    latrad := degtorad(lat);
    sdrad := degtorad(solardec);
    haarg :=
(cos(degtorad(90.833))/cos(latrad)*cos(sdrad)-tan(latrad)*tan(sdrad));
    ha :=
(arccos(cos(degtorad(90.833))/cos(latrad)*cos(sdrad)-tan(latrad)*tan(sdr
ad)));
    result := ha;
end;

function calcHourAngleSunset(lat,solarDec : extended) : extended;
//***************************************************************
// Purpose:  calculate the hour angle of the sun at sunset for the
latitude
//
//Arguments:  lat : latitude in degrees
//            solarDec : declination angle of sun in degreees
//
//Return:     hour angle of sunset in radians
//****************************************************************
var latrad,sdrad,haarg,ha : extended;
begin
    latrad := degtorad(lat);
    sdrad := degtorad(solardec);
    haarg :=
(cos(degtorad(90.833))/cos(latrad)*cos(sdrad)-tan(latrad)*tan(sdrad));
    ha :=
(arccos(cos(degtorad(90.833))/cos(latrad)*cos(sdrad)-tan(latrad)*tan(sdr
ad)));
    result := -ha;
end;

Function calcSunRiseUTC(jd,latitude,longitude : extended) :extended;
//***************************************************************
// Purpose:  calculate the UTC time of the sun at sunrise for the
latitude
//
//Arguments:  lat : latitude in degrees
//            solarDec : declination angle of sun in degreees
//
//Return:     time in minutes from zero
//****************************************************************
var
t,noonmin,tnoon,eqtime,solardec,hourangle,delta,timediff,timeutc,newt :
extended;
begin
    t := calcTimeJulianCent(JD);
    //find the time of solar noon at the location, and use
    //that declination.  This is better than start of the julian day
    noonmin := calcSolarNoonUTC(t,Longitude);
    tnoon := calcTimeJulianCent(jd+noonmin/1440);

    //** First pass to approximate sunrise using solar noon

    eqtime := calcEquationofTime(tnoon);
    solardec := calcSunDeclination(tnoon);
    hourAngle := calcHourAngleSunrise(latitude,solarDec);

    delta := longitude - radtodeg(hourangle);
    timediff := 4 * delta; //in minutes of time
    timeUTC := 720 + timeDiff - eqtime; // in minutes

    //** Second pass including the fractional jday in gamma calc

    newt := calctimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440);
    eqtime := calcEquationOfTime(newt);
    solarDec := calcSunDeclination(newt);
    hourangle := calcHourAngleSunrise(latitude,solarDec);
    delta := longitude - radtodeg(hourangle);
    timediff := 4 * delta;
    timeUTC := 720 + timeDiff - eqtime; //in minutes
    result := timeUTC;
end;

function calcSolarNoonUTC(t,longitude:extended) : Extended;
//***************************************************************
// Purpose:  calculate the UTC time of solar noon for a given
//           day and at the givin location on earth
//
//Arguments:  t : number of Julian centuries since J2000.0
//            longitude : longitude of boserver in degrees
//
//Return:     time in minutes from zero
//****************************************************************
var newt,eqtime,noonUTC : extended;
begin
    newt := calcTimeJulianCent(calcJDFromJulianCent(t) + 0.5 +
longitude/360);
    eqtime := calcEquationOfTime(newt);
    NoonUTC := 720 + (longitude *4) - eqtime; //minutes
    result := NoonUTC;
end;

Function calcSunSetUTC(jd,latitude,longitude : extended) :extended;
//***************************************************************
// Purpose:  calculate the UTC time of the sun at sunset for the
latitude
//
//Arguments:  lat : latitude in degrees
//            solarDec : declination angle of sun in degreees
//
//Return:     time in minutes from zero
//****************************************************************
var
t,noonmin,tnoon,eqtime,solardec,hourangle,delta,timediff,timeutc,newt :
extended;
begin
    t := calcTimeJulianCent(JD);
    //find the time of solar noon at the location, and use
    //that declination.  This is better than start of the julian day
    noonmin := calcSolarNoonUTC(t,Longitude);
    tnoon := calcTimeJulianCent(jd+noonmin/1440);

    //** First pass to approximate sunrise using solar noon

    eqtime := calcEquationofTime(tnoon);
    solardec := calcSunDeclination(tnoon);
    hourAngle := calcHourAngleSunset(latitude,solarDec);

    delta := longitude - radtodeg(hourangle);
    timediff := 4 * delta; //in minutes of time
    timeUTC := 720 + timeDiff - eqtime; // in minutes

    //** Second pass including the fractional jday in gamma calc

    newt := calctimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440);
    eqtime := calcEquationOfTime(newt);
    solarDec := calcSunDeclination(newt);
    hourangle := calcHourAngleSunset(latitude,solarDec);
    delta := longitude - radtodeg(hourangle);
    timediff := 4 * delta;
    timeUTC := 720 + timeDiff - eqtime; //in minutes
    result := timeUTC;
end;

Function calcTimeJulianCent(jd : extended) : extended;
//***************************************************************
// Purpose:  Convert JulianDay to Centruies since J2000.0
//
//Arguments:  jd : the julian day to convert
//
//
//Return:     the T value corresponding to the julian day
//****************************************************************
var t: extended;
begin
    t := (jd - 2451545.0)/36525.0;
    result := t;
end;

Function calcJDfromJulianCent(t : extended): Extended;
//***************************************************************
// Purpose:  Convert Centruies since J2000.0 to Julian Day
//
//Arguments:  t : number of Julian centruies since J2000;
//
//
//Return:     the Julian Day corresponding to the T value
//****************************************************************
var jd : extended;
begin
     jd := t * 36525.0 + 2451545.0;
     result := jd;
end;
Function SolarStuff(adate : TDateTime; alat,alon : extended) : Tsun;
//***************************************************************
// Purpose:  calculate the sunrise, sunset and solar noon for a given
//           date and location on the face of the earth
//
//Arguments:  adate : the date in TDatetime format
//            alat : latitude of location in degrees
//            alon : longitude of location in degrees
//
//Return:     Solar Sturcture containg sunrise, sunset and
//            high noon in TdateTime format;
//****************************************************************
var sunrise,sunset,noon,jd : extended;
    sun : Tsun;
begin

    jd := calcJD(yearof(adate),monthof(adate),dayof(adate));
    noon := calcSolarNoonUTC(calcTimeJulianCent(JD),alon);
    sunrise := calcSunRiseUTC(jd,alat,alon);
    sunset := calcSunSetUTC(jd,alat,alon);
    sun.srise := adate + (sunrise/1440);
    sun.sset := adate + (sunset/1440);
    sun.snoon := adate + (noon/1440);
    result := sun;
end;

Function calcJD(ayear,amonth,aday : integer) : extended;
//***************************************************************
// Purpose:  Convert calandar day to Julian
//
//Arguments:  ayear : 4 digit year
//            amonth : 2 digit month
//            aday : 2 digit dat
//
//Return:     the Julian Day
//****************************************************************
var a,b : integer;
    jd : extended;
begin
    if amonth <= 2 then
    begin
       ayear := ayear -1;
       amonth := amonth + 12;
    end;
    a := trunc(ayear/100);
    b := 2-a+trunc (a/4);
    jd := trunc(365.25 * (ayear + 4716)) + trunc(30.6001 *(amonth +1)) +
aday + b -1524.5;
    result := jd;
end;
end.



------------------------ Yahoo! Groups Sponsor ---------------------~-->
Buy Naturally Painless & Spray Away Backaches & Joint Pain. $19.97
http://www.challengerone.com/t/l.asp?cid=2867&lp=m331.html
http://us.click.yahoo.com/tJIe0D/79VGAA/ySSFAA/i7folB/TM
---------------------------------------------------------------------~->

CodeCoffer - new Delphi code protection Tool!
http://www.delphicollection.com/public/CodeCoffer.htm
---------------------------------------------------------------
Unsubscribe:[email protected]
List owner:[email protected]
--------------------------------------------------------------- 

Your use of Yahoo! Groups is subject to http://docs.yahoo.com/info/terms/