analemma.c position of sun in sky

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Arduino (C++) script to calculate position of sun in the sky   %
% Peter Stallinga, February 2012                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

void FloatToString(float f, char c[], int numdec){
int i, j;
unsigned char d;
float rond;

if (f<0){
c[0]='-';
c++;  // now isn't this an amazing coincidence?! C plus plus
f = abs(f);
}
rond = 0.4999999999;
for (j=0; j<numdec; j++) rond /= 10.0;
f += rond;
itoa((int) f, c, 10);
f -= (int) f;
while (*c!=0) c++;
if (numdec>0){
*c = '.';
c++;
for (j=0; j<numdec; j++){
f = 10.0*f;
d = (unsigned char) f;
*c = d + '0';
c++;
f = f-d;
}
*c = 0; // terminate string
}
}

void azialti(double y, double m, double d, double h, double mins,
double north, double east, double *azi, double *alti){
/*****************************************************************
*   Finds position of the sun in the sky: azimuth and altitude  *
*   returned in variables azi and alti                          *
*    azi: 0=South, -90=East, 90=West                            *
*    alti: 0=horizon, 90=Zenith                                 *
*   Input: y=year, m=month, d=day, h=hour, mins=minutes         *
*          north=latitude, east=longitude (east=positive)       *
*             unit: degrees:                                    *
*****************************************************************/
double L, xx, yy, a, g, t, epsilon, alpha, beta, lambda, theta0,
delta, theta, tau, alt, az;
double pi = 3.14159265358979323846;

// PART I: Calculate Sun position in starmap
// (based on http://www.stargazing.net/kepler/sun.html)
h = h + mins/60.0;
//1. Find the day:
d =  367 * y - floor(7 * (y +floor( (m + 9)/12))/ 4) + floor(275 * m / 9) + d - 730531.5 + h / 24;
//2. Find the Mean Longitude (L) of the Sun
L = 280.461 + 0.9856474 * d;
while (L<0)
L = L+360.0;
//3. Find the Mean anomaly (g) of the Sun
g = 357.528 + 0.9856003 * d;
//4. Find the ecliptic longitude (lambda) of the sun
lambda = L + 1.915 * sin(g*pi/180.0) + 0.020 * sin(2*g*pi/180.0);
//5. Find the obliquity of the ecliptic plane (epsilon)
epsilon = 23.439 - 0.0000004 * d;
//6. Find the Right Ascension (alpha) and Declination (delta) of  the Sun
yy = cos(epsilon*pi/180.0) * sin(lambda*pi/180.0);
xx = cos(lambda*pi/180.0);
a = (180.0/pi)*atan(yy/xx);
if (xx < 0)
alpha = a + 180.0;
else {
if ((yy < 0.0) && (xx > 0.0))
alpha = a + 360;
else
alpha = a;
}
delta = (180.0/pi)*asin(sin(epsilon*pi/180.0)*sin(lambda*pi/180.0));
// PART II: Calculate of position in sky
//compute Sidereal time at Greenwich:
//(based on http://www.geoastro.de/elevaz/basics/)
t = d/36525.0;
theta0 = 280.46061837 + 0.98564736629*d + 360.0*(d-floor(d)) + 0.000387933*t*t - t*t*t/38710000.0;
while (theta0<0.0)
theta0 = theta0+360.0;
while (theta0>360.0)
theta0 = theta0-360.0;
theta = theta0 + east;
tau = theta - alpha;
beta = pi*north/180.0;
delta = pi*delta/180.0;
tau = pi*tau/180.0;
alt = asin(sin(beta)*sin(delta) + cos(beta)*cos(delta)*cos(tau));
az = atan2(- sin(tau) , (cos(beta)*tan(delta) - sin(beta)*cos(tau)));
*alti = 180.0*alt/pi;
*azi = 180.0*az/pi-180.0;
if (*azi<-180.0)
*azi = *azi+360.0;
if (*azi>180.0)
*azi = *azi-360.0;
}

void setup()   {
Serial.begin(9600);
}

void loop()
{
double y = 2011.0;
double m = 12.0;
double d = 30.0;
double h = 9.0;
double mins = 0.0;
double north = 37.028176;
double east = -7.931132;

double azi, alti;
char s[80];

azialti(y, m, d, h, mins, north, east, &azi, &alti);
Serial.print("Azimuth: ");
FloatToString(azi, s, 5);
Serial.println(s);
Serial.print("Altitude: ");
FloatToString(alti, s, 5);
Serial.println(s);
delay(10000);
}