from pylab import *
from numpy import *
def sunsetangle(latitude,J):
    ds=0.4093*sin(2*pi/365*J-1.405)      # to be replaced by class
    sha=arccos(-tan(latitude)*tan(ds))
    return sha
latitude=48
J=arange(1,365,1)
plot(J,sunsetangle(latitude,J))
ytext = ylabel('sunset hour angle')
xtext = xlabel('Julian day')
show()