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