6 double xm, ym, zm, dxm, dym, dzm;
7 double xx, yx, zx, yy, zy, zz, tau;
8 double capt0, capt1, capt12, capt13, sl, sb, cl;
11 * remove E-terms of aberration
12 * except when finding catalog mean places
15 alpha += (.341/(3600.*15.))*sin((alpha+11.26)*15.*radian)
17 delta += (.341/3600.)*cos((alpha+11.26)*15.*radian)
18 *sin(delta*radian) - (.029/3600.)*cos(delta*radian);
21 * correct for proper motion
24 tau = (eday - epoch)/365.24220;
25 alpha += tau*da/3600.;
26 delta += tau*dd/3600.;
31 * convert to rectangular coordinates merely for convenience
34 xm = cos(delta)*cos(alpha);
35 ym = cos(delta)*sin(alpha);
39 * convert mean places at epoch of startable to current
40 * epoch (i.e. compute relevant precession)
43 capt0 = (epoch - 18262.427)/36524.220e0;
44 capt1 = (eday - epoch)/36524.220;
46 capt13 = capt12*capt1;
48 xx = - (.00029696+26.e-8*capt0)*capt12
50 yx = -(.02234941+1355.e-8*capt0)*capt1
51 - 676.e-8*capt12 + 221.e-8*capt13;
52 zx = -(.00971690-414.e-8*capt0)*capt1
53 + 207.e-8*capt12 + 96.e-8*capt13;
54 yy = - (.00024975+30.e-8*capt0)*capt12
56 zy = -(.00010858+2.e-8*capt0)*capt12;
57 zz = - (.00004721-4.e-8*capt0)*capt12;
59 dxm = xx*xm + yx*ym + zx*zm;
60 dym = - yx*xm + yy*ym + zy*zm;
61 dzm = - zx*xm + zy*ym + zz*zm;
68 * convert to mean ecliptic system of date
71 alpha = atan2(ym, xm);
72 delta = atan2(zm, sqrt(xm*xm+ym*ym));
73 cl = cos(delta)*cos(alpha);
74 sl = cos(delta)*sin(alpha)*cos(obliq) + sin(delta)*sin(obliq);
75 sb = -cos(delta)*sin(alpha)*sin(obliq) + sin(delta)*cos(obliq);
76 lambda = atan2(sl, cl);
77 beta = atan2(sb, sqrt(cl*cl+sl*sl));