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

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

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

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

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

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

#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