clear;

%calculate the equation of time for a full year

d2r = pi/180;
r2d = 180/pi;
w = 360/365.24;
d = 0:1:365;%days of the year
a = w * (d + 10);
b = a +(360/pi) * 0.0167 * sin(d2r * w * (d-2));
c = (a - r2d*atan(tan(d2r*b)/cos(23.44*d2r)))/180;
eoT = 720 * (c - round(c));
subplot(2,1,1);plot(d,eoT,'go-');hold on

%declination of the sun for a full year

%dec = r2d*sin(d2/365*2*pi);
dec = 23.45 * sin(d2r * (360 / 365 * (284 + d)));
subplot(2,1,1);plot(d,dec,'ro-');

%analemma = plot of eoT vs declination
subplot(2,1,2);plot(dec,eoT,'ko-');





