plan9port

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

helio.c (1629B)


      1 #include "astro.h"
      2 
      3 void
      4 helio(void)
      5 {
      6 /*
      7  *	uses lambda, beta, rad, motion
      8  *	sets alpha, delta, rp
      9  */
     10 
     11 /*
     12  *	helio converts from ecliptic heliocentric coordinates
     13  *	referred to the mean equinox of date
     14  *	to equatorial geocentric coordinates referred to
     15  *	the true equator and equinox
     16  */
     17 
     18 	double xmp, ymp, zmp;
     19 	double beta2;
     20 
     21 /*
     22  *	compute geocentric distance of object and
     23  *	compute light-time correction (i.i. planetary aberration)
     24  */
     25 
     26 	xmp = rad*cos(beta)*cos(lambda);
     27 	ymp = rad*cos(beta)*sin(lambda);
     28 	zmp = rad*sin(beta);
     29 	rp = sqrt((xmp+xms)*(xmp+xms) +
     30 		(ymp+yms)*(ymp+yms) +
     31 		(zmp+zms)*(zmp+zms));
     32 	lmb2 = lambda - .0057756e0*rp*motion;
     33 
     34 	xmp = rad*cos(beta)*cos(lmb2);
     35 	ymp = rad*cos(beta)*sin(lmb2);
     36 	zmp = rad*sin(beta);
     37 
     38 /*
     39  *	compute annual parallax from the position of the sun
     40  */
     41 
     42 	xmp += xms;
     43 	ymp += yms;
     44 	zmp += zms;
     45 	rp = sqrt(xmp*xmp + ymp*ymp + zmp*zmp);
     46 
     47 /*
     48  *	compute annual (i.e. stellar) aberration
     49  *	from the orbital velocity of the earth
     50  *	(by an incorrect method)
     51  */
     52 
     53 	xmp -= xdot*rp;
     54 	ymp -= ydot*rp;
     55 	zmp -= zdot*rp;
     56 
     57 /*
     58  *	perform the nutation and so convert from the mean
     59  *	equator and equinox to the true
     60  */
     61 
     62 	lmb2 = atan2(ymp, xmp);
     63 	beta2 = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
     64 	lmb2 += phi;
     65 
     66 /*
     67  *	change to equatorial coordinates
     68  */
     69 
     70 	xmp = rp*cos(lmb2)*cos(beta2);
     71 	ymp = rp*(sin(lmb2)*cos(beta2)*cos(tobliq) - sin(tobliq)*sin(beta2));
     72 	zmp = rp*(sin(lmb2)*cos(beta2)*sin(tobliq) + cos(tobliq)*sin(beta2));
     73 
     74 	alpha = atan2(ymp, xmp);
     75 	delta = atan2(zmp, sqrt(xmp*xmp+ymp*ymp));
     76 
     77 	hp = 8.794e0*radsec/rp;
     78 	semi /= rp;
     79 	if(rad > 0 && rad < 2.e5)
     80 		mag += 2.17*log(rad*rp);
     81 }