Trying to calculate sunrise...ain't getting the right answer

Question

This is my current code:

    public class Sunpos {
   final private double Pi = Math.PI;
   final private double eul = 2.71828182845904523552 ;
   final private double sonauf = 90;
   final private double RAD = 0.017453292519943295769236907684886;

public double sunrisefinal (double Breitengrad, double Laengengrad, int tagzahl, int sommerzeit, int nacht) {

    double lngHour = Laengengrad/15;
    double t = tagzahl + ((6 - lngHour)/24);

    // double ab = tagzahl + ((18 - lngHour)/24);
    double M = (0.9856 * t) - 3.289;
    double L = M + (1.916 * Math.sin(M)) + (0.020 * Math.sin(2 * M)) + 282.634;

        if (L >= 359) { L -= 360; }
        else if (L < 0) { L += 360; }  

    double RA = (Math.atan(0.91764 * Math.tan(Pi/180)*L));

        if (RA >= 359) { RA -= 360; }
        else if (RA < 0) { RA += 360; }

    double Lquadrant = (Math.floor(L/90)*90);
    double RAquadrant = (Math.floor(RA/90))*90;
    RA = RA + (Lquadrant - RAquadrant);
    RA = RA/15;

    double sinDec = 0.39782 * Math.sin((Pi/180)*L);
    double cosDec = (180/Pi)*(Math.cos(Math.asin(sinDec)));
    double cosH = (Math.cos((Pi/180)*sonauf)-(sinDec*Math.sin((Pi/180)*Breitengrad)))/(cosDec * Math.cos((Pi/180)*Breitengrad));
    double H = 360 - Math.acos(cosH);
    H /= 15;

    double T = H + RA -(0.06571 * t) - 6.622;
    double UTC = T - lngHour;

        if (UTC >= 23) { UTC -= 24; }
        else if (UTC < 0) { UTC += 24; }

    double locTime = UTC; // Fuer die schweiz!
    System.out.println(locTime);
    return(0);


}

The inputs are the following: ( 50, 10, 294, 1, 0). The last 2 can be ignored.

Now I am basing this on the following page: http://williams.best.vwh.net/sunrise_sunset_algorithm.htm

The code should be complete according to the site, but I don't get anywhere near the supposed results. I should get around 7.5 for today but I'm getting a 9.358.

Now, that might be because something with radiants/degrees? I can't quite get my Mind into that, as I've been trying to insert those converters (Pi/180) into the code, without any usable result.

Can anyone tell me where to put them or point me in the right direction? I've spent waaaay too much time on this already, and now I'm so close.

Answer 1

Everone seems to link to this http://williams.best.vwh.net/sunrise_sunset_algorithm.htm which doesn't exist anymore. Why not try something that gets updated once in a while like https://en.wikipedia.org/wiki/Sunrise_equation Then if you like you could help edit it to make it better.

Answer 2

I'll just post my implementation here in case people need it (ported from the same source as yours)

https://gist.github.com/zhong-j-yu/2232343b14a5b5ef5b9d

    public class SunRiseSetAlgo
    {
        static double calcSunrise(int dayOfYear, double localOffset, double latitude, double longitude)
        {
            return calc(dayOfYear, localOffset, latitude, longitude, true);
        }
        static double calcSunset(int dayOfYear, double localOffset, double latitude, double longitude)
        {
            return calc(dayOfYear, localOffset, latitude, longitude, false);
        }

        // http://williams.best.vwh.net/sunrise_sunset_algorithm.htm
        static double calc(int dayOfYear, double localOffset, double latitude, double longitude, boolean rise)
        {
            //1. first calculate the day of the year

    //        int N1 = floor(275 * month / 9.0);
    //        int N2 = floor((month + 9) / 12.0);
    //        int N3 = (1 + floor((year - 4 * floor(year / 4.0) + 2) / 3.0));
    //        int N = N1 - (N2 * N3) + day - 30;
            int N = dayOfYear;

            //2. convert the longitude to hour value and calculate an approximate time

            double lngHour = longitude / 15;
            double t = rise?
                N + (( 6 - lngHour) / 24) :
                N + ((18 - lngHour) / 24);

            //3. calculate the Sun's mean anomaly

            double M = (0.9856 * t) - 3.289;

            //4. calculate the Sun's true longitude

            double L = M + (1.916 * sin(M)) + (0.020 * sin(2 * M)) + 282.634;
            L = mod(L, 360);

            //5a. calculate the Sun's right ascension
            double RA = atan(0.91764 * tan(L));
            RA = mod(RA, 360);

            //5b. right ascension value needs to be in the same quadrant as L
            double Lquadrant  = (floor( L/90)) * 90;
            double 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
            double sinDec = 0.39782 * sin(L);
            double cosDec = cos(asin(sinDec));

            //7a. calculate the Sun's local hour angle
            double zenith = 90 + 50.0/60;
            double cosH = (cos(zenith) - (sinDec * sin(latitude))) / (cosDec * cos(latitude));

            if (cosH >  1)
              throw new Error("the sun never rises on this location (on the specified date");
            if (cosH < -1)
              throw new Error("the sun never sets on this location (on the specified date");

            //7b. finish calculating H and convert into hours
            double H = rise?
                360 - acos(cosH) :
                acos(cosH);
            H = H / 15;

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

            //9. adjust back to UTC
            double UT = T - lngHour;

            //10. convert UT value to local time zone of latitude/longitude
            double localT = UT + localOffset;
            localT = mod(localT, 24);
            return localT;
        }

        static int floor(double d){ return (int)Math.floor(d); }

        static double sin(double degree)
        {
            return Math.sin(degree*Math.PI/180);
        }
        static double cos(double degree)
        {
            return Math.cos(degree*Math.PI/180);
        }
        static double tan(double degree)
        {
            return Math.tan(degree*Math.PI/180);
        }
        static double atan(double x)
        {
            return Math.atan(x) *180/Math.PI;
        }
        static double asin(double x)
        {
            return Math.asin(x) *180/Math.PI;
        }
        static double acos(double x)
        {
            return Math.acos(x) *180/Math.PI;
        }

        static double mod(double x, double lim)
        {
            return x - lim * floor(x/lim);
        }

    }