analemma.m position of sun in sky

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Octave (Matlab) script to calculate position of sun in the sky %
% Peter Stallinga, February 2012                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
hold off;
clc; home;

% Part I: Inspired by http://www.stargazing.net/kepler/sun.html
%1. Find the days before J2000.0 (d) from the table

function [azi,alti] = azialti(y, m, d, h, mins, north, east)
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 = fmod(280.461 + 0.9856474 * d, 360.0);
L = 280.461 + 0.9856474 * d;
while (L<0)
L = L+360.0;
endwhile

%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);

%   (note that the sin(g) and sin(2*g) terms constitute an
%    approximation to the 'equation of centre' for the orbit
%    of the Sun)

% beta = 0;

%   (by definition as the Sun's orbit defines the
%             ecliptic plane. This results in a simplification
%             of the formulas below)

%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

Y = cos(epsilon*pi/180.0) * sin(lambda*pi/180.0);
X = cos(lambda*pi/180.0);

a = (180.0/pi)*atan(Y/X);

if (X < 0)
alpha = a + 180;
else
if (Y < 0) && (X > 0)
alpha = a + 360;
else
alpha = a;
endif
endif
alpha;

delta = (180.0/pi)*asin(sin(epsilon*pi/180.0)*sin(lambda*pi/180.0));

% Part II:
From http://www.geoastro.de/elevaz/basics/index.htm
%compute Sidereal time at Greenwich:
T = d/ 36525.0;
theta0 = 280.46061837 + 360.98564736629*d + 0.000387933*T*T - T*T*T/38710000.0;
while (theta0<0)
theta0 = theta0+360.0;
endwhile
theta0;
theta = theta0 + east;
tau = theta - alpha;

beta = north;
beta = pi*beta/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;
endif
if azi>180.0
azi=azi-360.0;
endif
endfunction

maanden = [31,28,31,30,31,30,31,31,30,31,30,31]; % months

y = 2011.0;
% Position of: Faro Vale da Amoreira
north = 37.028176;
east = -7.931132;

% part 1: analemmas:

for h=1:23
i = 0;
for n=1:365
d = n;
m = 1;
while (d>maanden(m))
d = d - maanden(m);
m = m+1;
endwhile

[az, alt] = azialti(y, m, d, h, 0, north, east);
%if ((alt>0) && (az<-30.0))
if (alt>0)
i=i+1;
az0h(h, i)=az;
alt0h(h, i)=alt;
endif
endfor
aserieslen(h)=i;
endfor

% part 2. Trajectory on day 21 of each month:

for m=1:12
n=0;
d = 21;
h=23.00;
mins=0;
alt = 100;
while (h>0)
[az,alt] = azialti(y, m, d, h, mins, north, east);
%if ((alt>0) && (az<-30.0))
if (alt>0)
n=n+1;
az21m(m,n)=az;
alt21m(m,n)=alt;
endif
mins=mins-1;
if (mins<0)
mins=mins+60;
h = h-1;
endif
endwhile
bserieslen(m)=n;
endfor

hold off;
for h=1:23
if aserieslen(h)>0
plot(az0h(h,1:aserieslen(h)), alt0h(h,1:aserieslen(h)),'bo');
hold on;
endif
endfor
for m=1:12
if bserieslen(m)>0
plot(az21m(m,1:bserieslen(m)), alt21m(m,1:bserieslen(m)),'k-');
endif
endfor
xlabel('Azimuth (degrees)');
ylabel('Altitude (degrees)');
axis('tight');

title('Position of Sun 2011. 0^o is South, -90^o is East, +90^o is West');