Solar Position

How to locate the sun

Referring to the Wikipedia articles Solar Zenith Angle and Solar Azimuth Angle, we have two equations. 'A' comes before 'Z', so lets attack azimuth first. Put into code, we have...


Azimuth = atan ( sin(HourAngle) / (cos(HourAngle) * sin(Latitude) - tan(Declination) * cos(Latitude)))

We need 3 inputs to this equation. The first input we need is HourAngle, and... what the heck is that? HourAngle can be defined as the current angle of rotation of the Earth, with respect to a solar day. It ranges from -180 to +180 degrees, with 0 degrees being solar noon. Luckily, AVR-Libc provides the solar_noon() function, from which we can derive the HourAngle for a point in time...
HourAngle

	#include <math.h>
	#include <time.h>
	time_t * point_in_time;
	time_t noon = solar_noon(point_in_time) % ONE_DAY; // time of noon mod one day
	time_t tday = *point_in_time % ONE_DAY;            // current time mod one day
	long r = tday - noon;                              // offset from noon, -43200 to + 43200 seconds
	double HourAngle = r / 86400.0;                    // hour angle as percentage of full day
	HourAngle = HourAngle * 2.0 * M_PI;                // hour angle in radians

The second thing we need is Declination, which is the current angle of the Sun above Earths equatorial plane. Again, AVR-Libc provides a function for that...
Declination


   double Declination = solar_declination(point_in_time);

Finally, we need the Latitude. AVR-Libc has a global long, __latitude, which we need to convert from seconds of arc into radians...
Latitude


  extern long __latitude;

  #define RAD_TO_DEGREE 57.29577951308233

  double Latitude =  __latitude / ( ONE_DEGREE * RAD_TO_DEGREE); // latitude in radians

Lets write a function to compute the azimuth...

Solar Azimuth

#include <math.h>
#include <time.h>
#define RAD_TO_DEGREE 57.29577951308233
double SolarAzimuth(time_t * point_in_time){

  time_t noon = solar_noon(point_in_time) % ONE_DAY;      // time of noon mod one day
  time_t tday = *point_in_time % ONE_DAY;                 // current time mod one day
  long r = tday - noon;                                   // offset from noon, -43200 to + 43200 seconds
  double HourAngle = r / 86400.0;                         // hour angle as percentage of full day
  HourAngle = HourAngle * 2.0 * M_PI;                     // hour angle in radians

  double Declination = solar_declination(point_in_time);  // declination in radians

  extern long  __latitude;
  double Latitude =  __latitude / ( ONE_DEGREE * RAD_TO_DEGREE);          // latitude in radians

  double Azimuth = atan ( sin(HourAngle) / (cos(HourAngle) * sin(Latitude) - tan(Declination) * cos(Latitude)));

  return Azimuth; // return azimuth in radians
}

For altitude, the equation is...


Altitude = asin( cos(HourAngle) * cos(Declination) * cos(Latitude) + sin(Declination) * sin(Latitude));

We can write the function...

Solar Altitude

double SolarAltitude(time_t * point_in_time){

  time_t noon = solar_noon(point_in_time) % ONE_DAY;      // time of noon mod one day
  time_t tday = *point_in_time % ONE_DAY;                 // current time mod one day
  long r = tday - noon;                                   // offset from noon, -43200 to + 43200 seconds
  double HourAngle = r / 86400.0;                         // hour angle as percentage of full day
  HourAngle = HourAngle * 2.0 * M_PI;                     // hour angle in radians

  double Declination = solar_declination(point_in_time);  // declination in radians

  extern long  __latitude;
  Latitude =  __latitude / ( ONE_DEGREE * RAD_TO_DEGREE);          // latitude in radians

  Altitude = asin( cos(HourAngle) * cos(Declination) * cos(Latitude) + sin(Declination) * sin(Latitude));

  return Altitude; // altitude in radians

}

Now, generally whenever we need the altitude, we also need the azimuth, & vice versa. There is a lot of common code between the two functions, so we could merge them...

Solar Position

#include <math.h>
#include <time.h>
#define RAD_TO_DEGREE 57.29577951308233
void SolarPosition(time_t * point_in_time, double * Azimuth, double * Altitude){


  time_t noon = solar_noon(point_in_time) % ONE_DAY;      // time of noon mod one day
  time_t tday = *point_in_time % ONE_DAY;                 // current time mod one day
  long r = tday - noon;                                   // offset from noon, -43200 to + 43200 seconds
  double HourAngle = r / 86400.0;                         // hour angle as percentage of full day
  HourAngle = HourAngle * 2.0 * M_PI;                     // hour angle in radians

  double Declination = solar_declination(point_in_time);  // declination in radians

  extern long  __latitude;
  double Latitude =  __latitude / ( ONE_DEGREE * RAD_TO_DEGREE);          // latitude in radians

  *Azimuth = atan ( sin(HourAngle) / (cos(HourAngle) * sin(Latitude) - tan(Declination) * cos(Latitude)));
  *Altitude = asin( cos(HourAngle) * cos(Declination) * cos(Latitude) + sin(Declination) * sin(Latitude));
}

Now we only have to compute HourAngle, Declination, and Latitude once. Next, lets look at all those trig functions. We compute many of these twice. Lets re-arrange to compute everything just once...

Solar Position (optimized)

#include <math.h>
#include <time.h>
#define RAD_TO_DEGREE 57.29577951308233
void SolarPosition(time_t * point_in_time, double * Azimuth, double * Altitude){


  time_t noon = solar_noon(point_in_time) % ONE_DAY;      // time of noon mod one day
  time_t tday = *point_in_time % ONE_DAY;                 // current time mod one day
  long r = tday - noon;                                   // offset from noon, -43200 to + 43200 seconds
  double HourAngle = r / 86400.0;                         // hour angle as percentage of full day
  HourAngle = HourAngle * 2.0 * M_PI;                     // hour angle in radians

  double Declination = solar_declination(point_in_time);  // declination in radians

  extern long  __latitude;
  double Latitude =  __latitude / ( ONE_DEGREE * RAD_TO_DEGREE);          // latitude in radians

  double CosHourAngle = cos(HourAngle);
  double SinLatitude = sin(Latitude);
  double CosLatitude = cos(Latitude);

  *Azimuth = atan ( sin(HourAngle) / ( CosHourAngle * SinLatitude - tan(Declination ) * CosLatitude ) );
  *Altitude = asin( CosHourAngle * cos(Declination) * CosLatitude + sin(Declination) * SinLatitude );
}

As noted in the Wiki articles, these algorithms are approximations. However, they do come quite close, the error being less than 1 degree over most of the Earth, at most times.

Error generally increases with latitude. In particular, where Latitude == 90 degrees - Declination, an attempted division by zero will produce wildly inaccurate results. We could test for this condition and correct for it, but we leave that as an exercise for the reader.

This algorithm also does not take into account atmospheric refraction. Causing light to bend around the horizon, the sun can appear above the horizion, when geometrically it is actually below.

The azimuth computed by this function uses the 'traditional' definition, where the azimuth is negative in the morning after sunrise, zero at solar noon, and positive in the afternoon... in other words the azimuth approximates the HourAngle.

Questions and comments can be sent to