plan9port

fork of plan9port with libvec, libstr and libsdb
Log | Files | Refs | README | LICENSE

star.c (1904B)


      1 #include "astro.h"
      2 
      3 void
      4 star(void)
      5 {
      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;
      9 
     10 /*
     11  *	remove E-terms of aberration
     12  *	except when finding catalog mean places
     13  */
     14 
     15 	alpha += (.341/(3600.*15.))*sin((alpha+11.26)*15.*radian)
     16 		  /cos(delta*radian);
     17 	delta += (.341/3600.)*cos((alpha+11.26)*15.*radian)
     18 		  *sin(delta*radian) - (.029/3600.)*cos(delta*radian);
     19 
     20 /*
     21  *	correct for proper motion
     22  */
     23 
     24 	tau = (eday - epoch)/365.24220;
     25 	alpha += tau*da/3600.;
     26 	delta += tau*dd/3600.;
     27 	alpha *= 15.*radian;
     28 	delta *= radian;
     29 
     30 /*
     31  *	convert to rectangular coordinates merely for convenience
     32  */
     33 
     34 	xm = cos(delta)*cos(alpha);
     35 	ym = cos(delta)*sin(alpha);
     36 	zm = sin(delta);
     37 
     38 /*
     39  *	convert mean places at epoch of startable to current
     40  *	epoch (i.e. compute relevant precession)
     41  */
     42 
     43 	capt0 = (epoch - 18262.427)/36524.220e0;
     44 	capt1 = (eday - epoch)/36524.220;
     45 	capt12 = capt1*capt1;
     46 	capt13 = capt12*capt1;
     47 
     48 	xx = - (.00029696+26.e-8*capt0)*capt12
     49 		  - 13.e-8*capt13;
     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
     55 		  - 15.e-8*capt13;
     56 	zy = -(.00010858+2.e-8*capt0)*capt12;
     57 	zz = - (.00004721-4.e-8*capt0)*capt12;
     58 
     59 	dxm =  xx*xm + yx*ym + zx*zm;
     60 	dym = - yx*xm + yy*ym + zy*zm;
     61 	dzm = - zx*xm + zy*ym + zz*zm;
     62 
     63 	xm = xm + dxm;
     64 	ym = ym + dym;
     65 	zm = zm + dzm;
     66 
     67 /*
     68  *	convert to mean ecliptic system of date
     69  */
     70 
     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));
     78 	rad = 1.e9;
     79 	if(px != 0)
     80 		rad = 20600/px;
     81 	motion = 0;
     82 	semi = 0;
     83 
     84 	helio();
     85 	geo();
     86 }