load'local_defs.ijs'
NB. Precess coordinates from one epoch to another.
NB. Usage: (d1 m1 y1; d2 m2 y2) precess h m s; d ' " --> result    
NB. where d1 m1 y1 is day,month(integer),year of first epoch, etc. 
NB. coordinates to be precessed are RA= h m s & Dec = d ' "
NB. Example: precess RA = 3h 45m 0s, Dec = 30deg 40min 0sec
NB. from 0 Jan 1950 to 4 Jul 2004. Enter:
NB.     (0 1 1950;4 7 2004) precess 3 45 0; 30 40 0 
NB. ---> result: 3 48 23.7964 | 30 50 0.272944
precess=: 4 : 0 
r=. (htr >{.y),(dtr >{:y)
v=.((cos {:r)*(cos {.r),(sin {.r)),sin {:r
'eF eD'=. (mjd >0{x),(mjd>1{x)
NB. Base epoch is J2000.0 = 2451545.0 <-> mjd 1.5 1 2000 = 51544.5
epoch=. 36525%~(eF-51544.5),(eD-eF)    NB. T & t in Julian centuries
x=. {. v =. (Pr_Mat epoch) mm v 
a=. rth (o.(x<0))+atan 1{v%x
NB. As a 2X3 matrix:  a=. 2 3$a, Asin {:v
a; Asin {:v
)

NB. Precession Matrix (see Seidelmann pp. 103-104)
Pr_Mat=: 3 : 0  NB. Th = theta, Z = z_A, Z0 = Zeta_A 
'T t'=. y       NB. T = Fixed - 2000, t = Date - Fixed
Z0=. 0,(2306.2181 1.39656 _0.000139 p. T),(0.30188+_0.000344*T),0.017998
cZ0=. cos Z0 =. dtr 0,0,(Z0 p. t)
Z=. 0,(2306.2181 1.39656 _0.000139 p. T),(1.09468+0.000066*T),0.018203
cZ=. cos Z =. dtr 0,0,(Z p. t)
Th=. 0,(2004.3109 _0.85330 _0.000217 p. T),(_0.42665+_0.000217*T),_0.041833
cT=. cos Th =. dtr 0,0,(Th p. t)
'sZ0 sZ sT'=. sin Z0, Z, Th
f=.(cZ0*cZ*cT)-sZ0*sZ
f=.f,(-(cZ0*sZ)+sZ0*cZ*cT),(-cZ*sT)
f=.f,(sZ0*cZ)+cZ0*sZ*cT
f=.f,(cZ0*cZ)-sZ0*sZ*cT
3 3$f,(-sZ*sT),(cZ0*sT),(-sZ0*sT),cT
)

NB. Modified Julian date: Usage: mjd d m yr ---> julian date
mjd=: 3 : 0  
m1=. m+12*i=. 3<~ m=. 1{y
y1=. (-i)+(yr<0)+ yr=. {:y
i=. 0<: #. |. *y -15 10 1582
b=. i*(2+(-a)+<. 4%~a=. <.y1%100)
c=. (*c)* <. |c=.(365.25*y1)-0.75*(y1<0)
_679006+b+c+({.y)+<.30.6*m1+1
NB. *modified* JD = JD - 2400000.5 (Green p 251, who dropped a 0)
)