r.b. goldstein  r.w. king  april 1978
routine to calculate the zenith angle and/or zenith
angle rate of a source
the routine needs to be called separately for each
site.

izctl:
1 bit=1: calculate za
2 bit=1: calculate zar
4 bit=1: calculate sitnrm (otherwise it is assumed to exis

note: the zenith angle is calculated using the site
normal, whereas the rate is calculated using xsite/r.
these two are different in that the site normal is
calculated on the basis of an elipsoidal shape.

x(6):  coordinates of the observed source w.r.t. the observing
u(3):  unit vector of x
r   :  magnitude of x
s(3):  planetocentric site coordinates
sr  :  magnitude of s
su(3):  coordinates of site normal