xref: /plan9/sys/src/cmd/astro/star.c (revision 59cc4ca53493a3c6d2349fe2b7f7c40f7dce7294)
13e12c5d1SDavid du Colombier #include "astro.h"
23e12c5d1SDavid du Colombier 
33e12c5d1SDavid du Colombier void
star(void)43e12c5d1SDavid du Colombier star(void)
53e12c5d1SDavid du Colombier {
63e12c5d1SDavid du Colombier 	double xm, ym, zm, dxm, dym, dzm;
73e12c5d1SDavid du Colombier 	double xx, yx, zx, yy, zy, zz, tau;
83e12c5d1SDavid du Colombier 	double capt0, capt1, capt12, capt13, sl, sb, cl;
93e12c5d1SDavid du Colombier 
103e12c5d1SDavid du Colombier /*
113e12c5d1SDavid du Colombier  *	remove E-terms of aberration
123e12c5d1SDavid du Colombier  *	except when finding catalog mean places
133e12c5d1SDavid du Colombier  */
143e12c5d1SDavid du Colombier 
153e12c5d1SDavid du Colombier 	alpha += (.341/(3600.*15.))*sin((alpha+11.26)*15.*radian)
163e12c5d1SDavid du Colombier 		  /cos(delta*radian);
173e12c5d1SDavid du Colombier 	delta += (.341/3600.)*cos((alpha+11.26)*15.*radian)
183e12c5d1SDavid du Colombier 		  *sin(delta*radian) - (.029/3600.)*cos(delta*radian);
193e12c5d1SDavid du Colombier 
203e12c5d1SDavid du Colombier /*
213e12c5d1SDavid du Colombier  *	correct for proper motion
223e12c5d1SDavid du Colombier  */
233e12c5d1SDavid du Colombier 
24*59cc4ca5SDavid du Colombier 	tau = (eday - epoch)/365.24220;
253e12c5d1SDavid du Colombier 	alpha += tau*da/3600.;
263e12c5d1SDavid du Colombier 	delta += tau*dd/3600.;
273e12c5d1SDavid du Colombier 	alpha *= 15.*radian;
283e12c5d1SDavid du Colombier 	delta *= radian;
293e12c5d1SDavid du Colombier 
303e12c5d1SDavid du Colombier /*
313e12c5d1SDavid du Colombier  *	convert to rectangular coordinates merely for convenience
323e12c5d1SDavid du Colombier  */
333e12c5d1SDavid du Colombier 
343e12c5d1SDavid du Colombier 	xm = cos(delta)*cos(alpha);
353e12c5d1SDavid du Colombier 	ym = cos(delta)*sin(alpha);
363e12c5d1SDavid du Colombier 	zm = sin(delta);
373e12c5d1SDavid du Colombier 
383e12c5d1SDavid du Colombier /*
393e12c5d1SDavid du Colombier  *	convert mean places at epoch of startable to current
403e12c5d1SDavid du Colombier  *	epoch (i.e. compute relevant precession)
413e12c5d1SDavid du Colombier  */
423e12c5d1SDavid du Colombier 
433e12c5d1SDavid du Colombier 	capt0 = (epoch - 18262.427)/36524.220e0;
443e12c5d1SDavid du Colombier 	capt1 = (eday - epoch)/36524.220;
453e12c5d1SDavid du Colombier 	capt12 = capt1*capt1;
463e12c5d1SDavid du Colombier 	capt13 = capt12*capt1;
473e12c5d1SDavid du Colombier 
483e12c5d1SDavid du Colombier 	xx = - (.00029696+26.e-8*capt0)*capt12
493e12c5d1SDavid du Colombier 		  - 13.e-8*capt13;
503e12c5d1SDavid du Colombier 	yx =  -(.02234941+1355.e-8*capt0)*capt1
513e12c5d1SDavid du Colombier 		  - 676.e-8*capt12 + 221.e-8*capt13;
523e12c5d1SDavid du Colombier 	zx = -(.00971690-414.e-8*capt0)*capt1
533e12c5d1SDavid du Colombier 		  + 207.e-8*capt12 + 96.e-8*capt13;
543e12c5d1SDavid du Colombier 	yy = - (.00024975+30.e-8*capt0)*capt12
553e12c5d1SDavid du Colombier 		  - 15.e-8*capt13;
563e12c5d1SDavid du Colombier 	zy = -(.00010858+2.e-8*capt0)*capt12;
573e12c5d1SDavid du Colombier 	zz = - (.00004721-4.e-8*capt0)*capt12;
583e12c5d1SDavid du Colombier 
593e12c5d1SDavid du Colombier 	dxm =  xx*xm + yx*ym + zx*zm;
603e12c5d1SDavid du Colombier 	dym = - yx*xm + yy*ym + zy*zm;
613e12c5d1SDavid du Colombier 	dzm = - zx*xm + zy*ym + zz*zm;
623e12c5d1SDavid du Colombier 
633e12c5d1SDavid du Colombier 	xm = xm + dxm;
643e12c5d1SDavid du Colombier 	ym = ym + dym;
653e12c5d1SDavid du Colombier 	zm = zm + dzm;
663e12c5d1SDavid du Colombier 
673e12c5d1SDavid du Colombier /*
683e12c5d1SDavid du Colombier  *	convert to mean ecliptic system of date
693e12c5d1SDavid du Colombier  */
703e12c5d1SDavid du Colombier 
713e12c5d1SDavid du Colombier 	alpha = atan2(ym, xm);
723e12c5d1SDavid du Colombier 	delta = atan2(zm, sqrt(xm*xm+ym*ym));
733e12c5d1SDavid du Colombier 	cl = cos(delta)*cos(alpha);
743e12c5d1SDavid du Colombier 	sl = cos(delta)*sin(alpha)*cos(obliq) + sin(delta)*sin(obliq);
753e12c5d1SDavid du Colombier 	sb = -cos(delta)*sin(alpha)*sin(obliq) + sin(delta)*cos(obliq);
763e12c5d1SDavid du Colombier 	lambda = atan2(sl, cl);
773e12c5d1SDavid du Colombier 	beta = atan2(sb, sqrt(cl*cl+sl*sl));
783e12c5d1SDavid du Colombier 	rad = 1.e9;
793e12c5d1SDavid du Colombier 	if(px != 0)
803e12c5d1SDavid du Colombier 		rad = 20600/px;
813e12c5d1SDavid du Colombier 	motion = 0;
823e12c5d1SDavid du Colombier 	semi = 0;
833e12c5d1SDavid du Colombier 
843e12c5d1SDavid du Colombier 	helio();
853e12c5d1SDavid du Colombier 	geo();
863e12c5d1SDavid du Colombier }
87