function degreedays(x0,h)
	x=x0; maxtime=30; numsteps=100000;
	for n=1:numsteps
    	if x<1
        	x=x+h*0.004*(28+5*cos(2*pi*(n-1)*h)-15);
    	else
        	break
    	end
	end
	
	x
	days=n*h
end

degreedays(0, 0.01)
degreedays(0, 0.001)
degreedays(-12, 0.002)
degreedays(-25, 0.005)