RE: Sunrise/Sunset

David Wright <[email protected]>
Newsgroups gmane.comp.lang.delphi.programming
Message-ID <[email protected]_la_us.celsinc.com> 
My head hurts from just glancing at this unit....

-----Original Message-----
From: Jeff Steinkamp - N7YG [mailto:[email protected]]
Sent: Wednesday, July 16, 2003 1:30 AM
To: Delphi Yahoo Group
Subject: [Delphi] 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	

ADVERTISEMENT
 <http://rd.yahoo.com/M=251812.3170658.4537139.1261774/D=egroupweb/S=1705115362:HM/A=1652964/R=0/SIG=11t2ts2ch/*http://www.netflix.com/Default?mqso=60178276&partid=3170658> click here	
  <http://us.adserver.yahoo.com/l?M=251812.3170658.4537139.1261774/D=egroupmail/S=:HM/A=1652964/rand=367080494> 	

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

Your use of Yahoo! Groups is subject to the Yahoo! Terms of Service <http://docs.yahoo.com/info/terms/> . 




[Non-text portions of this message have been removed]


------------------------ Yahoo! Groups Sponsor ---------------------~-->
Free shipping on all inkjet cartridge & refill kit orders to US & Canada. Low prices up to 80% off. We have your brand: HP, Epson, Lexmark & more.
http://www.c1tracking.com/l.asp?cid=5510
http://us.click.yahoo.com/GHXcIA/n.WGAA/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/
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.