sys/src/cmd/astro/helio.c

   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 }