function delta = declination_angle(doy)

% declination angle [Hartmann, Appendix A]
  theta_d = 2*pi*doy/365;
  delta   =   .006918 ...
	    - .399912 * cos(theta_d)   + .070257 * sin(theta_d) ...
	    - .006758 * cos(2*theta_d) + .000907 * sin(2*theta_d) ...
	    - .002697 * cos(3*theta_d) + .001480 * sin(3*theta_d);