sys/src/cmd/astro/star.c

   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 }