from pylab import *
def declination(J):
    ds=0.4093*sin(2*pi/365*J-1.405)
    return ds
J = arange(1,365,1)
plot(J,declination(J))
ytext = ylabel('solar declination in radians')
xtext = xlabel('Julian day')
show()