ash/forni  june 1968   subroutine occult
revisions:
slade 1971,  ng 1976,  chandler 1977
type of observation is specified by ncodf and further by ib
arguments

ncodet= ncodf - 8
ncodet=   -1   occultation
ncodet=    0   transit
ncodet.ge. 1   photographs as black dot crosses sun

ib       event                    logical center
ib=1 moon occults star                    observer
ib=2 planet occults star                  observer
ib=3 moon occults planet                  observer
ib=4 planet transits sun                  sun
ib=5 (not used)
ib=6 planet occults planet (mutual event) distant body
ib=7 (reserved for black dot photos)
ib=8 planet occults planet (mid-time)     observer
ib=9 planet eclipses planet (mid-time)    sun

xo   = vector pointing from earth to observed body
rr1  = xo(1)**2+xo(2)**2+xo(3)**2
r    = SQRT(rr1)
xz = vector pointing to second (distant) body
rz = SQRT(xz(1)**2+xz(2)**2+xz(3)**2 )

set up star direction unit vector
N.B.  This assumes star positions are given for mean epoch 1950,
which is not necessarily true.  More code will be needed if
lunar occultations are seriously used.

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*  start=3000
ncodet=    0   transit
vector   r     sub e =  xe(1-3)
vector   r dot sub e =  xe(4-6)
vector   r     sub m = xps(1-3)
vector   r dot sub m = xps(4-6)
transit of one planet across another (mutual occultation)
substitute coordinates wrt occulted body

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*  start=4000
mid-time of mutual occultation or eclipse
set up arrays for closest approach computation
mid-time occultation
phase correction to position of occulted body

rotate rhomn to selenographic frame and read watt's limb corect
note: some special treatment needed here for planet case
call prodct(rot,rhomn,rhomn0,3,3,1)
compute correction to limb
mrad=watlim(rhomn0)*mrad0
compute correction to occultation time