home *** CD-ROM | disk | FTP | other *** search
/ Gambler 19 / GAMBLERCD19.BIN / UTILS / RIVA128 / Uzytki / ssystem-1.2 / SOURCE / positions.c < prev    next >
C/C++ Source or Header  |  1998-05-02  |  8KB  |  252 lines
  1. /*
  2.  
  3.     Calculate planets positions in a given date.
  4.     
  5.     The main algorithm uses elements dated Julian 2450680.5 (8/20/1997) so 
  6.     it's quite accurate for dates near this.
  7.     
  8.     To be improved in the future ...
  9.     
  10.     Galilean and Saturn satellites computed using code extracted with
  11.     permission from XEphem, (c) 1997 Elwood Charles Downey.
  12.     
  13.     Keith Burnett's planet positioning algorithm used with persmission.
  14.     
  15. */
  16. #include "planets.h"
  17.  
  18. static void GalileanSat(void),SaturnSat(void),NeptuneSat(void),PlutoSat(void);
  19. static double SolveKepler(double, double);
  20.  
  21. /* Iterative solution of Kepler equation */
  22. static double SolveKepler(double M, double ec)
  23. {
  24.  static double E,d,delta;
  25.  
  26.  E=M;
  27.  d=E-ec*sin(E)-M;
  28.  while (fabs(d)>1.0E-10) {
  29.     delta=d/(1.0-ec*cos(E));
  30.     E=E-delta;
  31.     d=E-ec*sin(E)-M;
  32.  }
  33.  return E;
  34. }
  35.  
  36.  
  37. /* This function is nearly a clone of the jupiter_data function in xephem 3.0 
  38.    uses elements dated 12/31/1899
  39. */
  40. static void GalileanSat()
  41. {
  42. #define dsin(x) sin(DEG2RAD(x))
  43. #define dcos(x) cos(DEG2RAD(x))
  44.     static double A, B, Del, J, K, M, N, R, V;
  45.     static double cor_u1, cor_u2, cor_u3, cor_u4;
  46.     static double tmp, GG, H, psi, r, r1, r2, r3, r4;
  47.     static double u1, u2, u3, u4;
  48.     static double lam, Ds;
  49.     static double z1, z2, z3,  z4;
  50.     static double De, dsinDe,d;
  51.  
  52.     d=days+35295.0;
  53.     V = 134.63 + 0.00111587 * d;
  54.  
  55.     M = (358.47583 + 0.98560003*d);
  56.     N = (225.32833 + 0.0830853*d) + 0.33 * dsin (V);
  57.  
  58.     J = 221.647 + 0.9025179*d - 0.33 * dsin(V);;
  59.  
  60.     A = 1.916*dsin(M) + 0.02*dsin(2*M);
  61.     B = 5.552*dsin(N) + 0.167*dsin(2*N);
  62.     K = (J+A-B);
  63.     R = 1.00014 - 0.01672 * dcos(M) - 0.00014 * dcos(2*M);
  64.     r = 5.20867 - 0.25192 * dcos(N) - 0.00610 * dcos(2*N);
  65.     Del = sqrt (R*R + r*r - 2*R*r*dcos(K));
  66.     psi = RAD2DEG (asin(R/Del*dsin(K)));
  67.  
  68.     tmp = psi - B;
  69.  
  70.     u1 = 84.5506 + 203.4058630 * d + tmp;
  71.     u2 = 41.5015 + 101.2916323 * d + tmp;
  72.     u3 = 109.9770 + 50.2345169 * d + tmp;
  73.     u4 = 176.3586 + 21.4879802 * d + tmp;
  74.  
  75.     GG = 187.3 + 50.310674 * d;
  76.     H = 311.1 + 21.569229 * d;
  77.       
  78.     cor_u1 =  0.472 * dsin (2*(u1-u2));
  79.     cor_u2 =  1.073 * dsin (2*(u2-u3));
  80.     cor_u3 =  0.174 * dsin (GG);
  81.     cor_u4 =  0.845 * dsin (H);
  82.       
  83.     r1 = RADIUSSCALE(planets[JUPITER].Radius)*5.9061 - 0.0244 * dcos (2*(u1-u2));
  84.     r2 = RADIUSSCALE(planets[JUPITER].Radius)*9.3972 - 0.0889 * dcos (2*(u2-u3));
  85.     r3 = RADIUSSCALE(planets[JUPITER].Radius)*14.9894 - 0.0227 * dcos (GG);
  86.     r4 = RADIUSSCALE(planets[JUPITER].Radius)*26.3649 - 0.1944 * dcos (H);
  87.  
  88.     planets[IO].posx = r1 * dsin (u1+cor_u1);
  89.     planets[EUROPA].posx = r2 * dsin (u2+cor_u2);
  90.     planets[GANYMEDE].posx = r3 * dsin (u3+cor_u3);
  91.     planets[CALLISTO].posx = r4 * dsin (u4+cor_u4);
  92.  
  93.     lam = 238.05 + 0.083091*d + 0.33*dsin(V) + B;
  94.     Ds = 3.07*dsin(lam + 44.5);
  95.     De = Ds - 2.15*dsin(psi)*dcos(lam+24.)
  96.         - 1.31*(r-Del)/Del*dsin(lam-99.4);
  97.     dsinDe = dsin(De);
  98.  
  99.     z1 = r1 * dcos(u1+cor_u1);
  100.     z2 = r2 * dcos(u2+cor_u2);
  101.     z3 = r3 * dcos(u3+cor_u3);
  102.     z4 = r4 * dcos(u4+cor_u4);
  103.  
  104.     planets[IO].posy = z1*dsinDe;
  105.     planets[EUROPA].posy = z2*dsinDe;
  106.     planets[GANYMEDE].posy = z3*dsinDe;
  107.     planets[CALLISTO].posy = z4*dsinDe;
  108.  
  109.     planets[IO].posz = z1;
  110.     planets[EUROPA].posz = z2;
  111.     planets[GANYMEDE].posz = z3;
  112.     planets[CALLISTO].posz = z4;
  113.     
  114. #undef dsin
  115. #undef dcos
  116. }
  117.  
  118.  
  119. /* Another xephem based algorithm
  120.    Elements dated Jan 0.0 1980 = December 31,1979 0:0:0 UT
  121. */
  122. static void SaturnSat()
  123. {
  124.      static double SMA[11], U[11], U0[11], PD[11];
  125.      static double NN;
  126.      static int i;
  127.  
  128.      /*  Semimajor Axis of the Satellites' Orbit in planet radius */
  129.      SMA[1] = 3.08; SMA[2] = 3.95; SMA[3] = 4.89; SMA[4] = 6.26;
  130.      SMA[5] = 8.74; SMA[6] = 20.27; SMA[7] = 24.58; SMA[8] = 59.09;
  131.      /*  Eccentricity of Satellites' Orbit [Program uses 0] */
  132.      /*  Synodic Orbital Period of Moons in Days */
  133.      PD[1] = .9425049;        /* Mimas */
  134.      PD[2] = 1.3703731;        /* Enceladus */
  135.      PD[3] = 1.8880926;        /* Tethys */
  136.      PD[4] = 2.7375218;        /* Dione */
  137.      PD[5] = 4.5191631;        /* Rhea */
  138.      PD[6] = 15.9669028;    /* Titan */
  139.      PD[7] = 21.3174647;    /* Hyperion */
  140.      PD[8] = 79.9190206;    /* Iapetus */
  141.     
  142.      U0[1] = DEG2RAD(18.2919);
  143.      U0[2] = DEG2RAD(174.2135);
  144.      U0[3] = DEG2RAD(172.8546);
  145.      U0[4] = DEG2RAD(76.8438);
  146.      U0[5] = DEG2RAD(37.2555);
  147.      U0[6] = DEG2RAD(57.7005);
  148.      U0[7] = DEG2RAD(266.6977);
  149.      U0[8] = DEG2RAD(195.3513);
  150.  
  151.      NN=days+7305.0;
  152.  
  153.      for (i=TETHYS;i<=TITAN;i++) {
  154.         U[i-TETHYS+1] = U0[i-TETHYS+1] + (NN * 2 * PI / PD[i-TETHYS+1]);
  155.         U[i-TETHYS+1] = ((U[i-TETHYS+1] / (2 * PI)) - (int)(U[i-TETHYS+1] / (2 * PI))) * 2 * PI;
  156.          planets[i].posx = SMA[i-TETHYS+1] * sin(U[i-TETHYS+1]) * RADIUSSCALE(planets[SATURN].Radius);
  157.          planets[i].posz = SMA[i-TETHYS+1] * cos(U[i-TETHYS+1]) * RADIUSSCALE(planets[SATURN].Radius);
  158.          planets[i].posy = -SMA[i-TETHYS+1] * cos(U[i-TETHYS+1]) * sin(DEG2RAD(1.0)) * RADIUSSCALE(planets[SATURN].Radius);
  159.      }
  160. }
  161.  
  162.  
  163. static void NeptuneSat()
  164. {
  165.      static double SMA[11], U[11], U0[11], PD[11];
  166.      static double NN;
  167.  
  168.      /*  Semimajor Axis of the Satellites' Orbit in planet radius */
  169.      SMA[1] = 14.33;
  170.      /*  Eccentricity of Satellites' Orbit [Program uses 0] */
  171.      /*  Synodic Orbital Period of Moons in Days */
  172.      PD[1] = 5.876854;        /* Triton */
  173.     
  174.      U0[1] = 0.0; /* ? */
  175.  
  176.      NN=days;
  177.  
  178.      U[1] = U0[1] + (NN * 2 * PI / PD[1]);
  179.      U[1] = ((U[1] / (2 * PI)) - (int)(U[1] / (2 * PI))) * 2 * PI;
  180.      planets[TRITON].posx = SMA[1] * sin(U[1]) * RADIUSSCALE(planets[NEPTUNE].Radius);
  181.      planets[TRITON].posz = SMA[1] * cos(U[1]) * RADIUSSCALE(planets[NEPTUNE].Radius);
  182.      planets[TRITON].posy = -SMA[1] * cos(U[1]) * sin(DEG2RAD(1.0)) * RADIUSSCALE(planets[NEPTUNE].Radius);
  183. }
  184.  
  185.  
  186. static void PlutoSat()
  187. {
  188.      static double SMA[11], U[11], U0[11], PD[11];
  189.      static double NN;
  190.  
  191.      /*  Semimajor Axis of the Satellites' Orbit in planet radius */
  192.      SMA[1] = 17.0;
  193.      /*  Eccentricity of Satellites' Orbit [Program uses 0] */
  194.      /*  Synodic Orbital Period of Moons in Days */
  195.      PD[1] = 6.38725;        /* Charon */
  196.     
  197.      U0[1] = 0.0; /* ? */
  198.  
  199.      NN=days;
  200.  
  201.      U[1] = U0[1] + (NN * 2 * PI / PD[1]);
  202.      U[1] = ((U[1] / (2 * PI)) - (int)(U[1] / (2 * PI))) * 2 * PI;
  203.      planets[CHARON].posx = SMA[1] * RADIUSSCALE(planets[PLUTO].Radius) 
  204.                      * sin(U[1]);
  205.      planets[CHARON].posz = SMA[1] * RADIUSSCALE(planets[PLUTO].Radius)
  206.                      * cos(U[1]);
  207.      planets[CHARON].posy = -SMA[1] * RADIUSSCALE(planets[PLUTO].Radius) 
  208.                      * cos(U[1]) * sin(DEG2RAD(1.0));
  209. }
  210.  
  211.  
  212. /* Based on Keith Burnett's QBASIC code found here:
  213.     http://www.xylem.demon.co.uk/kepler/
  214. */
  215. void UpdatePositions(void)
  216. {
  217.  static int j;
  218.  static double e,M,E,r,v,o,p,i,x,y,z;
  219.  
  220.  for (j=1;j<=MOON;j++) {
  221.      planets[j].DeltaRotation+=planets[j].Rotation*timefactor;
  222.      e=planets[j].Eccentricity;
  223.      M=planets[j].DailyMotion*days+planets[j].MeanLongitude-planets[j].Perihelion;
  224.      E=SolveKepler(M,e);
  225.      r=planets[j].MeanDistance*(1.0-e*cos(E));
  226.      v=2.0*atan(sqrt((1+e)/(1-e))*tan(E/2.0));
  227.      o=planets[j].AscendingNode;
  228.      p=planets[j].Perihelion;
  229.      i=planets[j].Inclination;
  230.      planets[j].posx=r*(cos(o)*cos(v+p-o)-sin(o)*sin(v+p-o)*cos(i));
  231.      planets[j].posz=-r*(sin(o)*cos(v+p-o)+cos(o)*sin(v+p-o)*cos(i));
  232.      planets[j].posy=r*(sin(v+p-o)*sin(i));
  233.  }
  234.  GalileanSat();
  235.  SaturnSat();
  236.  NeptuneSat();
  237.  PlutoSat();
  238.  for (j=MOON;j<=CHARON;j++) {
  239.     x=planets[j].posx;
  240.     y=planets[j].posy;
  241.     z=planets[j].posz;
  242.      o=cos(DEG2RAD(planets[planets[j].Sat].Degrees)+planets[j].Inclination);
  243.      p=sin(DEG2RAD(planets[planets[j].Sat].Degrees)+planets[j].Inclination);
  244.      planets[j].posx = planets[planets[j].Sat].posx+x;
  245.          planets[j].posy = planets[planets[j].Sat].posy+y*o+z*p;
  246.          planets[j].posz = planets[planets[j].Sat].posz-y*p+z*o;
  247.  }
  248.  i=atan2(planets[EARTH].posz,planets[EARTH].posx)*180.0/PI;
  249.  planets[EARTH].DeltaRotation=(days-floor(days))*planets[EARTH].Rotation-i+90.0;
  250. }
  251.  
  252.