Wednesday, August 25, 2010

Sun Position in C#

I recently needed to calculate the Sun's position in the form of celestial coordinates, in this case Altitude and Azimuth angles. While there are a ton of examples explaining celestial coordinates and giving parts of the calculations, as well as a number of online applications for generating these coordinates, there are very few resources that are understandable to non-astronomers.

Finally after some wrestling I managed to create a function that produces the same output as these two sites:
http://pvcdrom.pveducation.org/SUNLIGHT/SUNCALC.HTM
http://www.usno.navy.mil/USNO/astronomical-applications/data-services/alt-az-world
(In the former, remember to subtract 1 from the hour field if you're calculating the angles for a date during the summer.. this is to account for daylight saving time. I don't remember if that is needed in the latter.

The calculations in my code are closely based on this site. I've tried to give as descriptive variable names as possible, but using the site for reference is still recommended if you want to understand the calculations. Currently the function simply prints the results, so you have to change that part to fit your needs. Note that the hour angle, right ascension and sun declination are also calculated in the function, if you are looking for any of those.

public static class SunPosition
{
    private const double Deg2Rad = Math.PI / 180.0;
    private const double Rad2Deg = 180.0 / Math.PI;

    /*!
     * \brief Calculates the sun light.
     *
     * CalcSunPosition calculates the suns "position" based on a
     * given date and time in local time, latitude and longitude
     * expressed in decimal degrees. It is based on the method
     * found here:
     * http://www.astro.uio.no/~bgranslo/aares/calculate.html
     * The calculation is only satisfiably correct for dates in
     * the range March 1 1900 to February 28 2100.
     * \param dateTime Time and date in local time.
     * \param latitude Latitude expressed in decimal degrees.
     * \param longitude Longitude expressed in decimal degrees.
     */
    public static void CalculateSunPosition(
        DateTime dateTime, double latitude, double longitude)
    {
        // Convert to UTC
        dateTime = dateTime.ToUniversalTime();

        // Number of days from J2000.0.
        double julianDate = 367 * dateTime.Year -
            (int)((7.0 / 4.0) * (dateTime.Year + 
            (int)((dateTime.Month + 9.0) / 12.0))) +
            (int)((275.0 * dateTime.Month) / 9.0) +
            dateTime.Day - 730531.5;
        
        double julianCenturies = julianDate / 36525.0;

        // Sidereal Time
        double siderealTimeHours = 6.6974 + 2400.0513 * julianCenturies;
        
        double siderealTimeUT = siderealTimeHours +
            (366.2422 / 365.2422) * (double)dateTime.TimeOfDay.TotalHours;
        
        double siderealTime = siderealTimeUT * 15 + longitude;

        // Refine to number of days (fractional) to specific time.
        julianDate += (double)dateTime.TimeOfDay.TotalHours / 24.0;
        julianCenturies = julianDate / 36525.0;

        // Solar Coordinates
        double meanLongitude = CorrectAngle(Deg2Rad *
            (280.466 + 36000.77 * julianCenturies));
        
        double meanAnomaly = CorrectAngle(Deg2Rad *
            (357.529 + 35999.05 * julianCenturies));
        
        double equationOfCenter = Deg2Rad * ((1.915 - 0.005 * julianCenturies) * 
            Math.Sin(meanAnomaly) + 0.02 * Math.Sin(2 * meanAnomaly));
        
        double elipticalLongitude =
            CorrectAngle(meanLongitude + equationOfCenter);
        
        double obliquity = (23.439 - 0.013 * julianCenturies) * Deg2Rad;

        // Right Ascension
        double rightAscension = Math.Atan2(
            Math.Cos(obliquity) * Math.Sin(elipticalLongitude),
            Math.Cos(elipticalLongitude));
        
        double declination = Math.Asin(
            Math.Sin(rightAscension) * Math.Sin(obliquity));

        // Horizontal Coordinates
        double hourAngle = CorrectAngle(siderealTime * Deg2Rad) - rightAscension;
        
        if (hourAngle > Math.PI)
        {
            hourAngle -= 2 * Math.PI;
        }

        double altitude = Math.Asin(Math.Sin(latitude * Deg2Rad) *
            Math.Sin(declination) + Math.Cos(latitude * Deg2Rad) *
            Math.Cos(declination) * Math.Cos(hourAngle));

        // Nominator and denominator for calculating Azimuth
        // angle. Needed to test which quadrant the angle is in.
        double aziNom = -Math.Sin(hourAngle);
        double aziDenom =
            Math.Tan(declination) * Math.Cos(latitude * Deg2Rad) -
            Math.Sin(latitude * Deg2Rad) * Math.Cos(hourAngle);
        
        double azimuth = Math.Atan(aziNom / aziDenom);
        
        if (aziDenom < 0) // In 2nd or 3rd quadrant
        {
            azimuth += Math.PI;
        }
        else if (aziNom < 0) // In 4th quadrant
        {
            azimuth += 2 * Math.PI;
        }

        // Altitude
        Console.WriteLine("Altitude: " + altitude * Rad2Deg);

        // Azimut
        Console.WriteLine("Azimuth: " + azimuth * Rad2Deg);
    }

    /*!
    * \brief Corrects an angle.
    *
    * \param angleInRadians An angle expressed in radians.
    * \return An angle in the range 0 to 2*PI.
    */
    private static double CorrectAngle(double angleInRadians)
    {
        if (angleInRadians < 0)
        {
            return 2 * Math.PI - (Math.Abs(angleInRadians) % (2 * Math.PI));
        }
        else if (angleInRadians > 2 * Math.PI)
        {
            return angleInRadians % (2 * Math.PI);
        }
        else
        {
            return angleInRadians;
        }
    }
}

5 comments:

  1. Thanks a lot for this information.
    I took only a quick look so far... perhaps the use of the % modulus operator in CorrectAngle introduces some errors - see http://msdn.microsoft.com/en-us/library/0w4e0fzs.aspx - "Avoid modulus operator with types float and double."

    ReplyDelete
  2. How can we convert these azimuth and altitude to cartesian (x,y,z) coordinates

    ReplyDelete
  3. This comment has been removed by the author.

    ReplyDelete
  4. Under which license do you publish this code? Do you allow usage of this function in closed source software?

    ReplyDelete