PRO lp_azel,ha,dec,az,el
;computes the azel of the sun in radians given ha in hours and dec in degrees
;  latitude of obs.
	la=28.758*!dtor
dr=dec*!dtor		;radian form of declination
hr=ha*(360.D0/24.D0)*!dtor	;and of hour angle
;use Ken's formula from gdr.ana
s1=sin(la)*sin(dr)
c1=cos(la)*cos(dr)*cos(hr)
xq=s1+c1
xq=xq<1.0   &	xq=xq>(-1.0)
el=asin(xq)
s1=sin(dr)-sin(la)*xq
c1=cos(la)*sqrt(1.-xq*xq)
xq=s1/c1
xq=xq<1.0   &	xq=xq>(-1.0)
az=acos(xq)
;Ken's formula loses the az sign, we can restore it using the ha sign, at
;least for La Palma
az=az+(ha gt 0)*(2.*!pi-2.*az)
az=az-!pi
end