## Sunrise/Sunset

A place to discuss the implementation and style of computer programs.

Moderators: phlip, Moderators General, Prelates

Paragon99
Posts: 16
Joined: Thu Sep 20, 2012 4:10 am UTC

### Sunrise/Sunset

I found an algorithm on line for calculating the Sunrise/Sunset... but I can't seem to get it implemented correctly.
The original algorithm is found
[url]williams.best.vwh.net/sunrise_sunset_algorithm.htm[/url]
It keeps telling me that sunrise is after noon....

Code: Select all

`    public static class DMath    {        //Extension method... Ie  Double d=5; d.ToRadian()==DMath.ToRadian(d);        public static double ToRadian(this double D)        {            return D * (double)Math.PI / 180;        }        public static double ToDegrees(this double R)        {            return R * (180 / (double)Math.PI);        }        public static double Sin(double D)        {            return Math.Sin(D.ToRadian()).ToDegrees();        }        public static double Cos(double D)        {            return Math.Cos(D.ToRadian()).ToDegrees();        }        public static double Tan(double D)        {            return Math.Tan(D.ToRadian()).ToDegrees();        }        public static double Atan(Double D)        {            return Math.Atan(D.ToRadian()).ToDegrees();        }        public static double Asin(Double D)        {            return Math.Asin(D.ToRadian()).ToDegrees();        }        public static double Acos(Double D)        {            return Math.Cos(D.ToRadian()).ToDegrees();        }    }`

Algorith

Code: Select all

`       public static DateTime Sunrise(Point LatLong)        {            double lat = LatLong.X;            double lng = LatLong.Y;            double localOffSet = -6;            double lngHour = lng / 15;                double trise = DateTime.Now.DayOfYear + ((6-lngHour)/24);                double M = (0.9856 * trise) - 3.289;                double L=M + (1.916* DMath.Sin(M)) + (.020*DMath.Sin(2*M)) + 282.634;                if (L>360)L-=360;                if (L<0) L+=360;                double RA = DMath.Atan(0.91764 * DMath.Tan(L));                if (RA>360) RA-=360;                if (RA<0) RA+=360;                double LQ = (Math.Floor(L / 90)) * 90;                double RAQ = (Math.Floor(RA / 90)) * 90;                RA=RA + (LQ-RAQ);                RA=RA/15;               double sinDec = 0.39782 * DMath.Sin(L);                double cosDec = DMath.Cos(DMath.Asin(sinDec));                double zenith=96;                double   cosH = (DMath.Cos(zenith) - (sinDec * DMath.Sin(lat))) / (cosDec * DMath.Cos(lat));               //Currently ignoring these cases               //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)                double H = 360 - DMath.Acos(cosH);//cosH.DegCos();               H = H / 15;                double T = H + RA - (0.06571 * trise) - 6.622;                 double UT = T - lngHour;                 if (UT > 24) UT -= 24;                 if (UT < 0) UT += 24;                 double localT = UT + localOffSet;                return DateTime.Now;        }`

PM 2Ring
Posts: 3700
Joined: Mon Jan 26, 2009 3:19 pm UTC
Location: Sydney, Australia

### Re: Sunrise/Sunset

Sorry, I don't know what language this is, but I think I see a few errors.
Paragon99 wrote:

Code: Select all

`    public static class DMath    {        //Extension method... Ie  Double d=5; d.ToRadian()==DMath.ToRadian(d);        public static double ToRadian(this double D)        {            return D * (double)Math.PI / 180;        }        public static double ToDegrees(this double R)        {            return R * (180 / (double)Math.PI);        }        public static double Sin(double D)        {            return Math.Sin(D.ToRadian()).ToDegrees();        }        public static double Atan(Double D)        {            return Math.Atan(D.ToRadian()).ToDegrees();        }`

ToDegrees() and ToRadian() are fine, but your trig functions look wrong. Eg sin takes an angle but returns a number, so

Code: Select all

`        public static double Sin(double D)        {            return Math.Sin(D.ToRadian());        }`

Similarly, atan takes a number and returns an angle, so

Code: Select all

`        public static double Atan(Double D)        {            return Math.Atan().ToDegrees();        }`

Also, your Sunrise() function looks like it returns the current time rather than the calculated time of sunrise.

Xenomortis
Not actually a special flower.
Posts: 1445
Joined: Thu Oct 11, 2012 8:47 am UTC

### Re: Sunrise/Sunset

It looks like C#, judging from the class names.

The first obvious bug is that you return DateTime.Now regardless of your input, rather than localT.

Secondly, to avoid confusion, I would immediately convert your longitude and latitude into radians and not bother with wrappers around the Math functions. You have errors in all of those (already pointed out).
Glancing through it, nothing appears to depend on using degrees over radians, other than the line:

Code: Select all

`double H = 360 - DMath.Acos(cosH);//cosH.DegCos();`
fixed by changing the 360 to 2*Math.Pi

There's a Math.Atan2 method which adjusts for the correct quadrant for you.

And naming a variable cosH is asking for trouble.

Edit:
Actually, a number of things do depend on it being degrees. Although I'm pretty sure it'll work out if you fix the quadrant calculation (90 -> 0.5 * Math.PI).
Edit 2:
L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
I don't know how this line is supposed to work right now, but that last constant looks suspiciously degree dependent. (It's supposed to give an angle by adding some (non-angle numbers) to an angle?)

PM 2Ring
Posts: 3700
Joined: Mon Jan 26, 2009 3:19 pm UTC
Location: Sydney, Australia

### Re: Sunrise/Sunset

Xenomortis wrote:L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634
I don't know how this line is supposed to work right now, but that last constant looks suspiciously degree dependent. (It's supposed to give an angle by adding some (non-angle numbers) to an angle?)

That line looks like it calculates the Sun's true ecliptic longitude from its mean anomaly. And yes, it is degree-dependent.

FWIW, I wrote a C program in 2000 to calculate various solar data including the Equation of Time, RA & declination, altitude & azimuth, rising & setting times, and shadow lengths, but the algorithm in my code is slightly different to that of the OP. I got my algorithm & data (for the 0/1/1983 epoch) from an old Astronomical Almanac of the USNO.

Paragon99
Posts: 16
Joined: Thu Sep 20, 2012 4:10 am UTC

### Re: Sunrise/Sunset

Thanks guys, the major issue was the incorrect Trig functions.

DeGuerre
Posts: 51
Joined: Mon Feb 04, 2008 6:41 am UTC

### Re: Sunrise/Sunset

Yeah, I wrote some code to do this a couple of years ago too, just to see how it's done.

I translated it into Java last night. Here it is in all of its glory. Note that my knowledge of java's Date and Calendar stuff dates from pre-1.1 days, so I may not have used this idiomatically. Note that this method is broken out into pieces that are easily googlable if you don't know what words like "nutation" mean.

Incidentally, you've made one seemingly minor (but actually very important) mistake: you've assigned latitude to be the X coordinate and longitude to be the Y coordinate.

Yes, I'm aware that latitude is traditionally given first. That's because a) it was easier to measure, and b) it was unambiguous since the coordinates are suffixed by N/S and E/W.

Yes, I'm aware that Google Maps does it that way. Google Maps also decrees that the Earth is spherical, presumably because they're Google, so they can always fix the Earth later. Google Maps was not designed by geographers, and got just about everything wrong that it's possible to get wrong.

There is a right way. ISO/IEC 13249 (which was written by geographers) does it the right way. Every GIS system does it the right way. That's the way you should do it too.

End rant. Here's the code:

Code: Select all

`package com.futplex.astronomy;import java.util.Date;import java.util.Calendar;import java.util.GregorianCalendar;import java.util.SimpleTimeZone;import java.util.TimeZone;public class SunriseSunset {   LonLat pos;   GregorianCalendar date;   final static double d2r(double degrees) {      return degrees * (Math.PI / 180.0);   }   final static double r2d(double radians) {      return radians * (180.0 / Math.PI);   }         /**    * Convert to local time.    *     * @param the day    * @param time the time in minutes since 0 UTC    * @return the local time    */   Date localTime(double t) {      long ms = Math.round(t * 60000);      ms += date.getTimeZone().getRawOffset();      int year = date.get(Calendar.YEAR);      int month = date.get(Calendar.MONTH);      int dom = date.get(Calendar.DAY_OF_MONTH);      int hour = (int)(ms / (1000 * 60 * 60));      ms %= 1000 * 60 * 60;      int minute = (int)(ms / (1000 * 60));      ms %= 1000 * 60;      int second = (int)(ms / 1000);      ms %= 1000;      Calendar time = new GregorianCalendar();      time.setTimeZone(date.getTimeZone());      time.set(year, month, dom, hour, minute, second);      time.set(Calendar.MILLISECOND, 0); // Not accurate to ms      return time.getTime();   }      /**    * Calculate the Julian Day Number from a Date.    */   public static double julianDayNumber(GregorianCalendar day) {      long year = day.get(Calendar.YEAR);      long month = day.get(Calendar.MONTH);      long dom = day.get(Calendar.DAY_OF_MONTH);      long gregorianCorrection = 2 - (year / 100) + (year / 400);            return (1461 * (year + 4716) / 4            + Math.floor(30.6001 * (month + 1))            + dom            + gregorianCorrection) - 1524.5;   }   /**    * Calculate the Julian century relative to 1 Jan 2000 UTC.    */   public static double julianCentury(double julianDay) {      final double j2000 = 2451545;       return (julianDay - j2000)/36525.0;   }      /**    * Constructor.    */    public SunriseSunset(GregorianCalendar date, LonLat pos) {       this.date = date;       this.pos = pos;    }    /**     * Mean obliquity of ecliptic.     *      * This equation is from:     *   J. Laskar, "New Formulas for the Precession, Valid Over 10000 years",     *   Astronomy and Astrophysics, Vol 157, p68 (1986).     */    static double meanObliquityOfEcliptic(double jc) {       double t = jc / 100.0;       if (t < -1.0 || t > 1.0) {          throw new RuntimeException("Need up-to-date astronomical data");       }       return (84381.448 + t *             (-4680.93 + t *             (-1.55 + t *             (1999.25 + t *             (-51.38 + t *             (-249.67 + t *             (-39.05 + t *              (-7.12 + t *              (28.87 + t *              (5.79 + t * 2.45)))))))))) / 3600.0;    }        /**     * True obliquity of the ecliptic.     */    static double trueObliquityOfEcliptic(double jc) {       double nutation = 0.00256 * Math.cos(d2r(125.04 - 1934.136 * jc));       return meanObliquityOfEcliptic(jc) + nutation;    }        static double adjustAngle(double angle) {       while (angle < 0) angle += 360;       while (angle >= 360) angle -= 360;       return angle;    }    static double solarGeometricMeanLongitude(double jc) {       return adjustAngle(280.46646 + jc * (36000.76983 + jc * 0.0003032));    }    static double solarGeometricMeanAnomaly(double jc) {       return 357.52911 + jc * (35999.05029 + jc * -0.0001537);    }        static double earthOrbitalEccentricity(double jc) {       return 0.016708634 + jc * (-0.000042037 + jc * -0.0000001267);           }        static double solarEquationOfCenter(double jc) {       double m = d2r(solarGeometricMeanAnomaly(jc));       return Math.sin(m) * (1.914602 + jc*(-0.004817 + jc * -0.000014))             + Math.sin(2*m) * (0.019993 + jc * -0.000101)             + Math.sin(3*m) * 0.00028;    }    static double solarTrueLongitude(double jc) {       return solarGeometricMeanLongitude(jc) + solarEquationOfCenter(jc);    }        static double solarApparentLongitude(double jc) {       double omega = 125.04 - 1934.136 * jc;       return solarTrueLongitude(jc) - 0.00569 - 0.00478 * Math.sin(d2r(omega));    }        /**     * Calculates the discrepancy between apparent solar time and mean solar time.     *      * (May be more familiar as the East-West component of the analemma.)     */    public static double equationOfTime(double jc) {       double l0r = d2r(solarGeometricMeanLongitude(jc));       double e = earthOrbitalEccentricity(jc);       double mr = d2r(solarGeometricMeanAnomaly(jc)), smr = Math.sin(mr);       double y = Math.tan(d2r(trueObliquityOfEcliptic(jc)) * 0.5);       double y2 = y*y;       return 4.0 * r2d(y2 * Math.sin(2.0 * l0r)             - 2.0 * e * smr             + 4.0 * e * y2 * smr * Math.cos(2.0 * l0r)             - 0.5 * y2 * y2 * Math.sin(4.0 * l0r)             - 1.25 * e * e * Math.sin(2.0 * mr)       );    }        static double solarDeclination(double jc) {       double apparentLong = d2r(solarApparentLongitude(jc));       double e = d2r(trueObliquityOfEcliptic(jc));       return r2d(Math.asin(Math.sin(e) * Math.sin(apparentLong)));    }        static double hourAngleSunrise(double lat, double solarDec) {       double latr = d2r(lat);       double sdr = d2r(solarDec);       return r2d(Math.acos(Math.cos(d2r(90.833)) / (Math.cos(latr) * Math.cos(sdr))             - Math.tan(latr) * Math.tan(sdr)));    }    /**     * Calculate sunrise or sunset in minutes since 0 UTC.     *      * @param jd The Julian day     * @param factor 1 for sunrise, -1 for sunset.     * @return the time     */    double calcSunriseSunsetUTC(double jd, double factor) {       final int ITERATIONS = 2;       double sunriseSunsetTime = 0;       for (int i = 0; i < ITERATIONS; ++i) {          double jc = julianCentury(jd + sunriseSunsetTime / 1440.0);          double hourAngle = hourAngleSunrise(pos.lat, solarDeclination(jc));          double et = equationOfTime(jc);          sunriseSunsetTime = 720.0 - 4 * (pos.lon + factor * hourAngle) - et;       }       return sunriseSunsetTime;    }    Date sunrise() {       return localTime(calcSunriseSunsetUTC(julianDayNumber(date), 1));    }    Date sunset() {       return localTime(calcSunriseSunsetUTC(julianDayNumber(date), -1));    }         public static void main(String argv[]) {       TimeZone tz = new SimpleTimeZone(10 * 1000 * 3600, "GMT+10");       GregorianCalendar cal = new GregorianCalendar();       cal.setTimeZone(tz);       cal.set(2013, 5, 21);       LonLat pos = new LonLat(145.2, -37.48);       SunriseSunset sun = new SunriseSunset(cal, pos);       System.out.println("Sunrise: " + sun.sunrise());       System.out.println("Sunset: " + sun.sunset());    }}`

thoughtfully
Posts: 2253
Joined: Thu Nov 01, 2007 12:25 am UTC
Location: Minneapolis, MN
Contact:

### Re: Sunrise/Sunset

If anyone is looking for a ready-made library, I've used NOVAS in the past. It's got Fortran, C, and now Python versions. I had some fun playing with SWIG to make Python bindings when I was using it. It's good to see it's being actively maintained.

Source code is available. It's the US gov't, so it's probably under pretty liberal terms.

Perfection is achieved, not when there is nothing more to add, but when there is nothing left to take away.
-- Antoine de Saint-Exupery

Xanthir
My HERO!!!
Posts: 5400
Joined: Tue Feb 20, 2007 12:49 am UTC
Contact:

### Re: Sunrise/Sunset

Everything produced by the US govt is, in fact, in the public domain. It's part of copyright law.
(defun fibs (n &optional (a 1) (b 1)) (take n (unfold '+ a b)))

thoughtfully
Posts: 2253
Joined: Thu Nov 01, 2007 12:25 am UTC
Location: Minneapolis, MN
Contact:

### Re: Sunrise/Sunset

I try to avoid using the term. It doesn't mean what you probably think it means. But in the case of works created by the gov't, it probably does.

Perfection is achieved, not when there is nothing more to add, but when there is nothing left to take away.
-- Antoine de Saint-Exupery