function [ delta ] = delta( n )
%DELTA declination angle at day n of the year
%   Uses a third degree Fourier expansion to compute the declination angle
%   of the earth on a given day. January 1 is day 1. The formula is taken 
%   from the book Global physical climatology by Dennis L. Hartmann (p.
%   348).

thetad = (2*pi*(n-1))./365;

% Fourier-coefficients
a0 = 0.006918;
a1 = -0.399912; b1 = 0.070257;
a2 = -0.006758; b2 = 0.000907;
a3 = -0.002697; b3 = 0.001480;

delta = (180/pi)*(a0 + a1*cos(thetad) + b1*sin(thetad) + a2*cos(2*thetad) + b2*sin(2*thetad) + a3*cos(3*thetad) + b3*sin(3*thetad));
