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.

1. 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.

2. 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----

-----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.

Answer 1

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.