thelio.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       thelio.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 }