uran.c 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154
  1. #include "astro.h"
  2. static double elem[] =
  3. {
  4. 36525.0, // [0] eday of epoc
  5. 19.19126393, // [1] semi major axis (au)
  6. 0.04716771, // [2] eccentricity
  7. 0.76986, // [3] inclination (deg)
  8. 74.22988, // [4] longitude of the ascending node (deg)
  9. 170.96424, // [5] longitude of perihelion (deg)
  10. 313.23218, // [6] mean longitude (deg)
  11. 0.00152025, // [1+6] (au/julian century)
  12. -0.00019150, // [2+6] (e/julian century)
  13. -2.09, // [3+6] (arcsec/julian century)
  14. -1681.40, // [4+6] (arcsec/julian century)
  15. 1312.56, // [5+6] (arcsec/julian century)
  16. 1542547.79, // [6+6] (arcsec/julian century)
  17. };
  18. void
  19. uran(void)
  20. {
  21. double pturbl, pturbb, pturbr;
  22. double lograd;
  23. double dele, enom, vnom, nd, sl;
  24. double capj, capn, eye, comg, omg;
  25. double sb, su, cu, u, b, up;
  26. double sd, ca, sa;
  27. double cy;
  28. cy = (eday - elem[0]) / 36525.; // per julian century
  29. mrad = elem[1] + elem[1+6]*cy;
  30. ecc = elem[2] + elem[2+6]*cy;
  31. cy = cy / 3600; // arcsec/deg per julian century
  32. incl = elem[3] + elem[3+6]*cy;
  33. node = elem[4] + elem[4+6]*cy;
  34. argp = elem[5] + elem[5+6]*cy;
  35. anom = elem[6] + elem[6+6]*cy - argp;
  36. motion = elem[6+6] / 36525. / 3600;
  37. incl *= radian;
  38. node *= radian;
  39. argp *= radian;
  40. anom = fmod(anom,360.)*radian;
  41. enom = anom + ecc*sin(anom);
  42. do {
  43. dele = (anom - enom + ecc * sin(enom)) /
  44. (1. - ecc*cos(enom));
  45. enom += dele;
  46. } while(fabs(dele) > converge);
  47. vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
  48. cos(enom/2.));
  49. rad = mrad*(1. - ecc*cos(enom));
  50. lambda = vnom + argp;
  51. pturbl = 0.;
  52. lambda += pturbl*radsec;
  53. pturbb = 0.;
  54. pturbr = 0.;
  55. /*
  56. * reduce to the ecliptic
  57. */
  58. nd = lambda - node;
  59. lambda = node + atan2(sin(nd)*cos(incl),cos(nd));
  60. sl = sin(incl)*sin(nd) + pturbb*radsec;
  61. beta = atan2(sl, pyth(sl));
  62. lograd = pturbr*2.30258509;
  63. rad *= 1. + lograd;
  64. lambda -= 1185.*radsec;
  65. beta -= 51.*radsec;
  66. motion *= radian*mrad*mrad/(rad*rad);
  67. semi = 83.33;
  68. /*
  69. * here begins the computation of magnitude
  70. * first find the geocentric equatorial coordinates of Saturn
  71. */
  72. sd = rad*(cos(beta)*sin(lambda)*sin(obliq) +
  73. sin(beta)*cos(obliq));
  74. sa = rad*(cos(beta)*sin(lambda)*cos(obliq) -
  75. sin(beta)*sin(obliq));
  76. ca = rad*cos(beta)*cos(lambda);
  77. sd += zms;
  78. sa += yms;
  79. ca += xms;
  80. alpha = atan2(sa,ca);
  81. delta = atan2(sd,sqrt(sa*sa+ca*ca));
  82. /*
  83. * here are the necessary elements of Saturn's rings
  84. * cf. Exp. Supp. p. 363ff.
  85. */
  86. capj = 6.9056 - 0.4322*capt;
  87. capn = 126.3615 + 3.9894*capt + 0.2403*capt2;
  88. eye = 28.0743 - 0.0128*capt;
  89. comg = 168.1179 + 1.3936*capt;
  90. omg = 42.9236 - 2.7390*capt - 0.2344*capt2;
  91. capj *= radian;
  92. capn *= radian;
  93. eye *= radian;
  94. comg *= radian;
  95. omg *= radian;
  96. /*
  97. * now find saturnicentric ring-plane coords of the earth
  98. */
  99. sb = sin(capj)*cos(delta)*sin(alpha-capn) -
  100. cos(capj)*sin(delta);
  101. su = cos(capj)*cos(delta)*sin(alpha-capn) +
  102. sin(capj)*sin(delta);
  103. cu = cos(delta)*cos(alpha-capn);
  104. u = atan2(su,cu);
  105. b = atan2(sb,sqrt(su*su+cu*cu));
  106. /*
  107. * and then the saturnicentric ring-plane coords of the sun
  108. */
  109. su = sin(eye)*sin(beta) +
  110. cos(eye)*cos(beta)*sin(lambda-comg);
  111. cu = cos(beta)*cos(lambda-comg);
  112. up = atan2(su,cu);
  113. /*
  114. * at last, the magnitude
  115. */
  116. sb = sin(b);
  117. mag = -8.68 +2.52*fabs(up+omg-u)-
  118. 2.60*fabs(sb) + 1.25*(sb*sb);
  119. helio();
  120. geo();
  121. }