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; } } }

Thanks a lot for this information.

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

Thank you!

ReplyDeleteHow can we convert these azimuth and altitude to cartesian (x,y,z) coordinates

ReplyDeleteThis comment has been removed by the author.

ReplyDeleteUnder which license do you publish this code? Do you allow usage of this function in closed source software?

ReplyDelete