%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);
}
For more information, contact me at The University of The
Algarve,
Prof. Peter Stallinga
http://www.stallinga.org