Calculations of the path of the Sun

Question

I'm writing several methods necessary to calculate the path of the Sun across a specific point. I have written the code using two different sources for my calculations and neither is producing the desired result. The sources are: http://www.pveducation.org/pvcdrom/properties-of-sunlight/suns-position and http://www.esrl.noaa.gov/gmd/grad/solcalc/solareqns.PDF

Note: Degrees to arcminutes is Deg * 60 min.

localSolartime: I have converted the longitude to 'minutes', the local standard time meridian(lstm) derived from the localStandardTimeMeridian method returns a value that is in 'minutes', and the equationOfTime which is also returned in 'minutes'. Using the equation from pveducation, I've calculated the time correction which accounts for the small time variations within a given time zone. When I apply this result and the localTime, each in minutes, to the local solar time (lst) equation, the result is 676.515 (at this moment), which does not make any sense to me. The local solar time, as I understand it, represents the time with respect to the Sun and when it is at its highest point in the sky, locally, is considered solar noon. 676.515 does not make sense. Does anybody understand what might be causing this.
HourAngle: I'm hoping that once I fix the localSolarTime method, this will not need to be corrected.

I've chosen Washington DC for the latitude and longitude. Both the Zenith and Azimuth readings should be positive values, and for my region at this moment, are 66 and 201 respectively.

public class PathOfSun {
    static LocalTime localTime = LocalTime.now();
    static double dcLat = 38.83;
    static double dcLong =  -77.02;
    static DecimalFormat df = new DecimalFormat("#.0");

    public static void main(String [] args) {
        int day = dayOfYear();
        double equationOfTime = equationOfTime(day);
        double lstm = localTimeMeridian();
        double lst = localSolarTime(equationOfTime, dcLong, lstm);
        double declination = declination(day);
        double hourAngle = hourAngle(lst);

        double zenith = zenith(dcLat, declination, hourAngle);
        double azimuth = azimuth(dcLong, declination, zenith, hourAngle); 

    }

    //Longitude of timezone meridian
    public static double localTimeMeridian() {
        TimeZone gmt = TimeZone.getTimeZone("GMT");
        TimeZone est = TimeZone.getTimeZone("EST");
        int td = gmt.getRawOffset() - est.getRawOffset();
        double localStandardTimeMeridian = 15 * (td/(1000*60*60)); //convert td to hours
        //System.out.println("Local Time Meridian: " + localStandardTimeMeridian);
        return localStandardTimeMeridian;
    }

    //Get the number of days since Jan. 1
    public static int dayOfYear() {
        Calendar localCalendar = Calendar.getInstance(TimeZone.getDefault());
        int dayOfYear = localCalendar.get(Calendar.DAY_OF_YEAR); 
        //System.out.println("Day: " + dayOfYear);
        return dayOfYear;
    }

    //Emperical equation to correct the eccentricity of Earth's orbit and axial tilt
    public static double equationOfTime (double day) {
        double d =(360.0/365.0)*(day - 81);
        d = Math.toRadians(d);
        double equationTime = 9.87*sin(2*d)-7.53*cos(d)-1.54*sin(d); 
        //System.out.println("Equation Of Time: " + equationTime);
        return equationTime;
    }
    //The angle between the equator and a line drawn from the center of the Sun(degrees)
    public static double declination(int dayOfYear) {
        double declination = 23.5*sin((Math.toRadians(360.0/365.0))*(dayOfYear - 81));
        //System.out.println("Declination: " + df.format(declination));
        return declination;
    }

    //Add the number of minutes past midnight localtime//
    public static double hourAngle(double localSolarTime) {
        double hourAngle = 15 * (localSolarTime - 13); 
        System.out.println("Hour Angle: " + df.format(hourAngle)); //(degrees)
        return hourAngle;
    }

    //Account for the variation within timezone - increases accuracy
    public static double localSolarTime(double equationOfTime, double longitude, double lstm) { 
        //LocalSolarTime = 4min * (longitude + localStandardTimeMeridian) + equationOfTime
        //Time Correction is time variation within given time zone (minutes)
        //longitude = longitude/60; //convert degrees to arcminutes
        double localStandardTimeMeridian = lstm;
        double timeCorrection = (4 * (longitude + localStandardTimeMeridian) + equationOfTime);
        System.out.println("Time Correction: " + timeCorrection); //(in minutes)
        //localSolarTime represents solar time where noon represents sun's is highest position 
        // in sky and the hour angle is 0 -- hour angle is negative in morning, and positive after solar noon.
        double localSolarTime = (localTime.toSecondOfDay() + (timeCorrection*60)); //(seconds)
        localSolarTime = localSolarTime/(60*60);  //convert from seconds to hours
        //Convert double to Time (HH:mm:ss) for console output
        int hours = (int) Math.floor(localSolarTime);
        int minutes = (int) ((localSolarTime - hours) * 60);
        //-1 for the daylight savings
        Time solarTime = new Time((hours-1), minutes, 0);
        System.out.println("Local Solar Time: " + solarTime); //hours

        return localSolarTime;
    }

    public static double azimuth(double lat, double declination, double zenith, double hourAngle) {
        double azimuthDegree = 0;
        double elevation = 90 - zenith;
        elevation = Math.toRadians(elevation);
        zenith = Math.toRadians(zenith);
        lat = Math.toRadians(lat);
        declination = Math.toRadians(declination);
        hourAngle = Math.round(hourAngle);
        hourAngle = Math.toRadians(hourAngle);

        //double azimuthRadian = -sin(hourAngle)*cos(declination) / cos(elevation);
        double azimuthRadian = ((sin(declination)*cos(lat)) - (cos(hourAngle)*cos(declination)*
                sin(lat)))/cos(elevation);

        //Account for time quadrants
        Calendar cal = Calendar.getInstance();
        int hour = cal.get(Calendar.HOUR_OF_DAY);
        if(hour > 0 && hour < 6) {
        azimuthDegree =  Math.toDegrees(acos(azimuthRadian));
        }
        else if(hour >= 6 && hour < 12) {
            azimuthDegree = Math.toDegrees(acos(azimuthRadian));
            azimuthDegree = 180 - azimuthDegree;
        } else if (hour >= 12 && hour < 18) {
            azimuthDegree = Math.toDegrees(acos(azimuthRadian));
            azimuthDegree = azimuthDegree - 180;
        } else if (hour >= 18 && hour < 24) {
            azimuthDegree = Math.toDegrees(acos(azimuthRadian));
            azimuthDegree = 360 - azimuthDegree;
        }

        System.out.println("Azimuth: " + df.format(azimuthDegree));
        return azimuthDegree;
    }

    public static double zenith(double lat, double declination, double hourAngle) {
        lat = Math.toRadians(lat);
        declination = Math.toRadians(declination);
        hourAngle = Math.round(hourAngle);
        hourAngle = Math.toRadians(hourAngle);
        //Solar Zenith Angle 
        double zenith = Math.toDegrees(acos(sin(lat)*sin(declination) + (cos(lat)*cos(declination)*cos(hourAngle))));
        //Solar Elevation Angle
        double elevation = Math.toDegrees(asin(sin(lat)*sin(declination) + (cos(lat)*cos(declination)*cos(hourAngle))));
        System.out.println("Elevation: " + df.format(elevation));
        System.out.println("Zenith: " + df.format(zenith));
        return zenith;
    }
}

Just to reiterate, the day, local time meridian are exactly correct, and the equation of time and declination are accurate but not exact. ----UPDATE OUTPUT---- new output

sensor program

-----UPDATE----- Used the scatterchart to display the sun's elevation/azimuth throughout day. I am still having trouble figuring out the azimuth output. It is correct for long time, but then it will change from increasing and start to decrease (~270-->0). I will be sure to update the code once I finally get the output right.

Tough to understand what the exact problem is due to the wording. You say it's supposed to produce minutes and then give the results in minutes; why?? You say the two should be positive and then give positive results; why?? Have you debugged this to relate the inputs and outputs? — ChiefTwoPencils, May 04 '15 at 19:42
Your zenith() and azimuth() functions do not reflect the formulas given in the citations you provide. Did you try to simplify them? — BadZen, May 04 '15 at 19:48
Also, please note that you will have to mod the results [0,2*Pi] to get a "normalized" (non-negative) angle. — BadZen, May 04 '15 at 19:50
@BadZen, sorry about that, I think I must have rewritten the equations wrong. I will update the post to clarify. — wellington, May 04 '15 at 20:01
There are also another approaches to this task [calculate the time when the sun is X degrees below/above the Horizon](http://stackoverflow.com/a/25722336/2521214) and also this might help a bit [Calculate whether it is close to dawn/dusk based on sunrise/sunset?](http://stackoverflow.com/a/27038880/2521214) — Spektre, May 05 '15 at 07:57

erickson · Accepted Answer · 2015-05-07T15:38:35.313

3

You pass the longitude to localSolarTime() as degrees, and then you divide that by 60, with a comment claiming this is in order to convert to minutes of arc. This is wrong; your later calculations require degrees, and even if you needed minutes of arc, you'd multiply by 60, not divide.

This mistaken division results in a longitude of -1.3°, and when you find the angle between your local time meridian and your position, you get a large angle (about 75°). It should be a small angle, generally ±7.5°. The large angle results in a large time correction, and throws everything off.

Update: In the updated version of the azimuth() method, the quadrant selection should be based on the hour angle of the sun, or, equivalently, on local solar time, rather than standard wall clock time. And, the hour angle used in all calculations should not be rounded. Rather than testing four different quadrants, the method could look like this:

public static double azimuth(double lat, double declination, double zenith, double hourAngle)
{
  double elevation = Math.toRadians(90 - zenith);
  lat = Math.toRadians(lat);
  declination = Math.toRadians(declination);
  hourAngle = Math.toRadians(hourAngle);
  double azimuthRadian = acos(((sin(declination) * cos(lat)) - (cos(hourAngle) * cos(declination) * sin(lat))) / cos(elevation));
  double azimuthDegree = Math.toDegrees(azimuthRadian);
  if (hourAngle > 0)
    azimuthDegree = 360 - azimuthDegree;
  System.out.println("Azimuth: " + df.format(azimuthDegree));
  return azimuthDegree;
}

Finally, you are passing dcLong in as the lat parameter of the azimuth() method; this should be dcLat.

I'd recommend using radians internally throughout, and only converting from and to degrees on input and output. This will help prevent mistakes, and cut down on rounding errors and unnecessary clutter.

edited May 07 '15 at 15:38

answered May 05 '15 at 18:40

erickson

249,448
50
371
469

Thanks erickson for pointing that out. I will adjust the code and hopefully it will produce something a little more pleasant! – wellington May 05 '15 at 18:52
I've worked through the problem and the methods work correctly now, except the final calculation - azimuth angle is incorrect. In the morning, the range of the azimuth angle is 0-180 (from 0-12) and 180-360 in the afternoon(12-24hr). Given that, I need to adjust the value within the acos( ) to consider the change in sign for each quadrant of the unit circle. I've tried using two different equations, both from http://en.wikipedia.org/wiki/Solar_azimuth_angle, but output of the second equations produces a '?'. The first equation is fairly close but is about 10 degrees different than actual. – wellington May 06 '15 at 15:59
@wellington If you need additional help, add an update to your post with your current azimuth method and sample inputs and expected results. – erickson May 06 '15 at 17:17
Do you have any recommendations on how I could reduce the amount of code? I feel like a caveman and that there isn't a lot of finesse in my coding. :) I appreciate your help. – wellington May 06 '15 at 17:30
1

Your last problem is that you are passing longitude to your azimuth method instead of the expected latitude. As for reducing the amount of code, you are doing work you don't need to. For example, the four if-else cases in azimuth really should be `return (hour > 0) ? 360 - azimuth : azimuth;` And you should handle everything as radians internally, converting only on input and output to degrees if needed. Use `java.time` APIs consistently: `OffsetDateTime`, `LocalTime`, and `Duration`. Write some unit tests. Write your tests to use a hard-coded time rather than `now()`. – erickson May 06 '15 at 23:53
Yeah, I found the mistake of using the wrong time value (i.e. hourAngle) for the azimuth degree angle. That def helped to produce a better value. I would like to say that, atm, I use 180 not 360, in order to get the azimuth value that matches the actual. This might be bc of the direction(?): http://www.pveducation.org/pvcdrom/properties-of-sunlight/azimuth-angle So far this is working, but we will see at the critical time (when hourAngle < 0). Just doing the degree conversions for troubleshooting. Really helpful info erickson, appreciate all of your help! – wellington May 07 '15 at 20:15

Calculations of the path of the Sun

1 Answers1

Linked