home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 1 / FFMCD01.bin / bbs / libdisks / d700t799 / disk709.lha / Planets / planet_pos.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-08-05  |  32.6 KB  |  769 lines
  1. /*******************************************************************************
  2. ** Compute the position of the major planets and return the epli value.       **
  3. *******************************************************************************/
  4.  
  5. #include <stdio.h>
  6. #include <math.h>
  7.  
  8. #define RAD 0.01745329251994329651
  9.  
  10. /*******************************************************************************
  11. ** Crude magnitudes of planets (x100) for output filtering by brightness.     **
  12. ***************************************************************************** */
  13.  
  14. #define MAGSOL -9.9
  15. #define MAGMER -1.5
  16. #define MAGVEN -4.0
  17. #define MAGMAR -1.0
  18. #define MAGJUP -2.0
  19. #define MAGSAT  1.0
  20. #define MAGURA  5.9
  21. #define MAGNEP  8.0
  22.  
  23. double planet_pos(jd,T0)
  24. double jd,T0;
  25.    {
  26.    double M0,M1,M2,M4,M5,M6,M7,M8,D,G,H,M,N,P,Q,S,V,W,RA,DEC,ECC;
  27.    double sinQ,cosQ,sin2Q,cos2Q,sin3Q,cos3Q,sin4Q,cos4Q;
  28.    double sinze,cosze,sin2ze,cos2ze,sin3ze,cos3ze,sin4ze,cos4ze;
  29.    double sinps,cosps,sin2ps,cos2ps,sin3ps,cos3ps;
  30.    double sinH,cosH,sin2H,cos2H,sinG,cosG;
  31.    double sinV,cosV,sin2V,cos2V,sinS,cosS,sin2S,cos2S;
  32.    double w1,w2,a,b,e,i,l,ll,r,u,nu,nu2,psi,eta,th,ze,epli,thapp,omeg;
  33.    double l1pert,w1pert;
  34.    double esun,Cen,Stheta,Sr;
  35.  
  36.    double kepler(),truean(),longi(),lati(),poly(),range();
  37.    void   helio_trans(),speak();
  38.  
  39. /*******************************************************************************
  40. ** Calculate orbital elements for Mercury.                                    **
  41. *******************************************************************************/
  42.    l = range(poly(178.179078,149474.07078,0.0003011,0.,T0));
  43.  
  44.    a = 0.3870986;
  45.  
  46.    e = poly(0.20561421,0.00002046,-0.000000030,0.,T0);
  47.    i = poly(7.00288100,0.00186080,-0.000018300,0.,T0);
  48.  
  49.    w1 = range(poly(28.753753,0.3702806,0.0001208,0.,T0));
  50.    w2 = range(poly(47.145944,1.1852083,0.0001739,0.,T0));
  51.  
  52.    M1 = range(poly(102.27938,149472.51529,0.000007,0.,T0));
  53.  
  54. /*******************************************************************************
  55. ** Solving Kepler find the eccentric anomaly ECC.                             **
  56. *******************************************************************************/
  57.    ECC = kepler(e,M1);
  58.    nu  = truean(e,ECC);
  59.    r   = a * (1. - e * cos(ECC*RAD));
  60.    u   = l + nu - M1 - w2;
  61.    ll  = longi(w2,i,u);
  62.    b   = lati(u,i);
  63.  
  64. /*******************************************************************************
  65. ** Now make corrections due to perturbations.                                 **
  66. *******************************************************************************/
  67.    M2 = range(poly(212.60322,58517.80387, 0.001286,0.,T0));
  68.    M4 = range(poly(319.51913,19139.85475, 0.000181,0.,T0));
  69.    M5 = range(poly(225.32833, 3034.69202,-0.000722,0.,T0));
  70.    M6 = range(poly(175.46622, 1221.55147,-0.000502,0.,T0));
  71.  
  72.    ll += 0.00204 * cos((-2. * M1 + 5. * M2 +  12.220) * RAD)
  73.        + 0.00103 * cos((     -M1 + 2. * M2 - 160.692) * RAD)
  74.        + 0.00091 * cos((     -M1 + 2. * M5 -  37.003) * RAD)
  75.        + 0.00078 * cos((-3. * M1 + 5. * M2 +  10.137) * RAD);
  76.  
  77.    r += 0.000007525 * cos((     -M1 + 2. * M5 +  53.013) * RAD)
  78.       + 0.000006802 * cos((-3. * M1 + 5. * M2 - 259.918) * RAD)
  79.       + 0.000005457 * cos((-2. * M1 + 2. * M2 -  71.188) * RAD)
  80.       + 0.000003569 * cos((     -M1 + 5. * M2 -  77.750) * RAD);
  81.  
  82. /*******************************************************************************
  83. ** Now find the Longi and radius vector of the sun.                           **
  84. *******************************************************************************/
  85.    M = range(poly(358.47583,35999.0497500,-0.0001500,-0.0000033,T0));
  86.  
  87.    esun = poly(0.01675104,-0.0000418,-0.000000126,0.,T0);
  88.  
  89.    Cen = poly(1.919460,-0.004789,-0.000014,0.,T0) * sin(     M * RAD) +
  90.          poly(0.020094,-0.000100, 0.      ,0.,T0) * sin(2. * M * RAD) +
  91.          poly(0.000293, 0.      , 0.      ,0.,T0) * sin(3. * M * RAD);
  92.  
  93.    Stheta = range(Cen + range(poly(279.69668,36000.76892,0.0003025,0.,T0)));
  94.  
  95.    Sr = 1.0000002 * (1. - esun * esun) / (1. + esun * cos((M + Cen) * RAD));
  96.  
  97.    omeg  = 259.18 - 1934.142 * T0;
  98.    thapp = Stheta - 0.00569 - 0.00479 * sin(omeg * RAD);
  99.    epli  = poly(23.45229400,-0.013012500,-0.000001640,0.000000503,T0);
  100.  
  101.    N   = cos(epli * RAD) * sin(thapp * RAD);
  102.    D   = cos(thapp * RAD);
  103.    RA  = atan2(N,D) / RAD;
  104.    DEC = asin(sin(epli * RAD) * sin(thapp * RAD)) / RAD;
  105.  
  106. #ifdef EUTSCH
  107.    speak(RA,DEC,Sr,MAGSOL,"PS","Sonne");
  108. #else
  109.    speak(RA,DEC,Sr,MAGSOL,"PS","Sol");
  110. #endif
  111.  
  112. /*******************************************************************************
  113. ** Tansformation of coordinates on Mercury and output.                        **
  114. *******************************************************************************/
  115. #ifdef EUTSCH
  116.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGMER,"PM","Merkur");
  117. #else
  118.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGMER,"PM","Mercury");
  119. #endif
  120.  
  121. /*******************************************************************************
  122. ** Now start on Venus.                                                        **
  123. *******************************************************************************/
  124.    l = range(poly(342.767053,58519.21191,0.0003097,0.,T0));
  125.  
  126.    a = 0.7233316;
  127.  
  128.    e = poly(0.00682069,-0.00004774, 0.000000091,0.,T0);
  129.    i = poly(3.39363100, 0.00100580,-0.000001000,0.,T0);
  130.  
  131.    w1 = range(poly(54.384186,0.5081861,-0.0013864,0.,T0));
  132.    w2 = range(poly(75.779647,0.8998500, 0.0004100,0.,T0));
  133.  
  134. /*******************************************************************************
  135. ** Venus has a long period pert that needs to be added before Kelper.         **
  136. *******************************************************************************/
  137.    M0  = M2 + 0.00077 * sin((237.24 + 150.27 * T0) * RAD);
  138.  
  139.    l  += M0 - M2;
  140.    ECC = kepler(e,M0);
  141.    nu  = truean(e,ECC);
  142.    r   = a * (1. - e * cos(ECC * RAD));
  143.    u   = l + nu - M0 - w2;
  144.    ll  = longi(w2,i,u);
  145.    b   = lati(u,i);
  146.  
  147. /*******************************************************************************
  148. ** Now Venus perturbations.                                                   **
  149. *******************************************************************************/
  150.    ll += 0.00313 * cos((2. * M  - 2. * M2 - 148.225) * RAD)
  151.        + 0.00198 * cos((3. * M  - 3. * M2 +   2.565) * RAD)
  152.        + 0.00136 * cos((     M  -      M2 - 119.107) * RAD)
  153.        + 0.00096 * cos((3. * M  - 2. * M2 - 135.912) * RAD)
  154.        + 0.00082 * cos((     M5 -      M2 - 208.087) * RAD);
  155.  
  156.    r += 0.000022501 * cos((2. * M  - 2. * M2 -  58.208) * RAD)
  157.       + 0.000019045 * cos((3. * M  - 3. * M2 +  92.577) * RAD)
  158.       + 0.000006887 * cos((     M5 -      M2 - 118.090) * RAD)
  159.       + 0.000005172 * cos((     M  -      M2 -  29.110) * RAD)
  160.       + 0.000003620 * cos((5. * M  - 4. * M2 - 104.208) * RAD)
  161.       + 0.000003283 * cos((4. * M  - 4. * M2 +  63.513) * RAD)
  162.       + 0.000003074 * cos((2. * M5 - 2. * M2 -  55.167) * RAD);
  163.  
  164. /*******************************************************************************
  165. ** Tansformation of coordinates on Venus and output.                          **
  166. *******************************************************************************/
  167.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGVEN,"PV","Venus");
  168.  
  169. /*******************************************************************************
  170. ** Now start the planet Mars.                                                 **
  171. *******************************************************************************/
  172.    l = range(poly(293.737334,19141.69551,0.0003107,0.,T0));
  173.  
  174.    a = 1.5236883;
  175.  
  176.    e = poly(0.09331290, 0.000092064,-0.000000077,0.,T0);
  177.    i = poly(1.85033300,-0.000675000, 0.000012600,0.,T0);
  178.  
  179.    w1 = range(poly(285.431761,1.0697667, 0.0001313, 0.00000414,T0));
  180.    w2 = range(poly( 48.786442,0.7709917,-0.0000014,-0.00000533,T0));
  181.  
  182. /*******************************************************************************
  183. ** Mars has a long period perturbation.                                       **
  184. *******************************************************************************/
  185.    M0  = M4 + -0.01133 * sin((3. * M5 - 8. * M4 + 4. * M) * RAD)
  186.               -0.00933 * cos((3. * M5 - 8. * M4 + 4. * M) * RAD);
  187.  
  188.    l  += M0 - M4;
  189.    ECC = kepler(e,M0);
  190.    nu  = truean(e,ECC);
  191.    r   = a * (1. - e * cos(ECC * RAD));
  192.    u   = l + nu - M0 - w2;
  193.    ll  = longi(w2,i,u);
  194.    b   = lati(u,i);
  195.  
  196. /*******************************************************************************
  197. ** Now Mars Perturbations.                                                    **
  198. *******************************************************************************/
  199.    ll += 0.00705 * cos((     M5 -      M4 -  48.958) * RAD)
  200.        + 0.00607 * cos((2. * M5 -      M4 - 188.350) * RAD)
  201.        + 0.00445 * cos((2. * M5 - 2. * M4 - 191.897) * RAD)
  202.        + 0.00388 * cos((     M  - 2. * M4 +  20.495) * RAD)
  203.        + 0.00238 * cos((     M  -      M4 +  35.097) * RAD)
  204.        + 0.00204 * cos((2. * M  - 3. * M4 + 158.638) * RAD)
  205.        + 0.00177 * cos((    -M2 + 3. * M4 -  57.602) * RAD)
  206.        + 0.00136 * cos((2. * M  - 4. * M4 + 154.093) * RAD)
  207.        + 0.00104 * cos((     M5           +  17.618) * RAD);
  208.  
  209.    r += 0.000053227 * cos((     M5 -      M4 +  41.1306) * RAD)
  210.       + 0.000050989 * cos((2. * M5 - 2. * M4 - 101.9847) * RAD)
  211.       + 0.000038278 * cos((2. * M5 -      M4 -  98.3292) * RAD)
  212.       + 0.000015996 * cos((     M  -      M4 -  55.5550) * RAD)
  213.       + 0.000014764 * cos((2. * M  - 3. * M4 +  68.6220) * RAD)
  214.       + 0.000008966 * cos((     M5 - 2. * M4 +  43.6150) * RAD)
  215.       + 0.000007914 * cos((3. * M5 - 2. * M4 - 139.7370) * RAD)
  216.       + 0.000007004 * cos((2. * M5 - 3. * M4 - 102.8880) * RAD)
  217.       + 0.000006620 * cos((     M  - 2. * M4 + 113.2020) * RAD)
  218.       + 0.000004930 * cos((3. * M5 - 3. * M4 -  76.2430) * RAD)
  219.       + 0.000004693 * cos((3. * M  - 5. * M4 + 190.6030) * RAD)
  220.       + 0.000004571 * cos((2. * M  - 4. * M4 + 244.7020) * RAD)
  221.       + 0.000004409 * cos((3. * M5 -      M4 - 115.8280) * RAD);
  222.  
  223. /*******************************************************************************
  224. ** Tansformation of coordinates on Mars and output.                           **
  225. *******************************************************************************/
  226.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGMAR,"Pm","Mars");
  227.  
  228. /*******************************************************************************
  229. ** Now start Jupiter.                                                         **
  230. *******************************************************************************/
  231.    l = range(poly(238.049257,3036.301986,0.0003347,-0.00000165,T0));
  232.  
  233.    a = 5.202561;
  234.  
  235.    e = poly(0.04833475, 0.000164180,-0.0000004676,-0.0000000017,T0);
  236.    i = poly(1.30873600,-0.005696100, 0.0000039000, 0.0000000000,T0);
  237.  
  238.    w1 = range(poly(273.277558,0.5994317,0.00070405, 0.00000508,T0));
  239.    w2 = range(poly( 99.443414,1.0105300,0.00035222,-0.00000851,T0));
  240.  
  241. /*******************************************************************************
  242. ** Now start perturbation calclations.                                        **
  243. *******************************************************************************/
  244.    nu2 = T0 / 5. + 0.1;
  245.  
  246.    P = 237.47555 + 3034.9061 * T0;
  247.    Q = 265.91650 + 1222.1139 * T0;
  248.    S = 243.51721 +  428.4677 * T0;
  249.  
  250.    V = range(5. * Q - 2. * P) * RAD;
  251.    W = range(2. * P - 6. * Q + 3. * S) * RAD;
  252.  
  253.    ze  = range(Q - P) * RAD;
  254.    psi = range(S - Q) * RAD;
  255.  
  256.    P = range(P) * RAD;
  257.    Q = range(Q) * RAD;
  258.    S = range(S) * RAD;
  259.  
  260. /*******************************************************************************
  261. ** An extreme acceleration can be achieved by using the addtion theorems.     **
  262. *******************************************************************************/
  263.    sinQ = sin(Q);
  264.    cosQ = cos(Q);
  265.  
  266.    sin2Q = 2. * sinQ * cosQ;
  267.    cos2Q = cosQ*cosQ - sinQ*sinQ;
  268.  
  269.    sin3Q = (  3. - 4. * sinQ * sinQ) * sinQ;
  270.    cos3Q = (- 3. + 4. * cosQ * cosQ) * cosQ;
  271.  
  272.    sin4Q = (8. * cosQ * cosQ - 4.) * cosQ * sinQ;
  273.    cos4Q =  8. * cosQ * cosQ * ( cosQ * cosQ - 1. ) + 1.;
  274.  
  275. /*******************************************************************************/
  276.    sinze = sin(ze);
  277.    cosze = cos(ze);
  278.  
  279.    sin2ze = 2. * sinze * cosze;
  280.    cos2ze = cosze*cosze - sinze*sinze;
  281.  
  282.    sin3ze = (  3. - 4. * sinze * sinze) * sinze;
  283.    cos3ze = (- 3. + 4. * cosze * cosze) * cosze;
  284.  
  285.    sin4ze = (8. * cosze * cosze - 4.) * cosze * sinze;
  286.    cos4ze =  8. * cosze * cosze * (cosze * cosze - 1. ) + 1.;
  287.  
  288. /*******************************************************************************/
  289.    sinV = sin(V);
  290.    cosV = cos(V);
  291.  
  292.    sin2V = 2. * sinV * cosV;
  293.    cos2V = cosV * cosV - sinV * sinV;
  294.  
  295. /*******************************************************************************/
  296.    sinS  = sin(S);           cosS  = cos(S);
  297.    sin2S = 2. * sinS * cosS; cos2S = cosS * cosS - sinS * sinS;
  298.  
  299. /*******************************************************************************/
  300.    sinps = sin(psi);
  301.    cosps = cos(psi);
  302.  
  303.    sin2ps = 2. * sinps * cosps;
  304.    cos2ps = cosps*cosps - sinps*sinps;
  305.  
  306.    sin3ps = (  3. - 4. * sinps * sinps) * sinps;
  307.    cos3ps = (- 3. + 4. * cosps * cosps) * cosps;
  308.  
  309. /*******************************************************************************
  310. ** Calculating the perturbation terms.                                        **
  311. *******************************************************************************/
  312.    l1pert = poly(0.331364,-0.010281,-0.004692,0.,nu2) * sinV
  313.           + poly(0.003228,-0.064436, 0.002075,0.,nu2) * cosV
  314.           - poly(0.003083, 0.000275,-0.000489,0.,nu2) * sin2V
  315.           +  0.002472                                 * sin(W)
  316.           +  0.013619                   * sinze
  317.           +  0.018472                   * sin2ze
  318.           +  0.006717                   * sin3ze
  319.           +  0.002775                   * sin4ze
  320.           + (0.007275 - 0.001253 * nu2) * sinze  * sinQ
  321.           +  0.006417                   * sin2ze * sinQ
  322.           +  0.002439                   * sin3ze * sinQ
  323.           - (0.033839 + 0.001125 * nu2) * cosze  * sinQ
  324.           -  0.003767                   * cos2ze * sinQ
  325.           - (0.035681 + 0.001208 * nu2) * sinze  * cosQ
  326.           -  0.004261                   * sin2ze * cosQ
  327.           +  0.002178                            * cosQ
  328.           - (0.006333 - 0.001161 * nu2) * cosze  * cosQ
  329.           -  0.006675                   * cos2ze * cosQ
  330.           -  0.002664                   * cos3ze * cosQ
  331.           -  0.002572                   * sinze  * sin2Q
  332.           -  0.003567                   * sin2ze * sin2Q
  333.           +  0.002094                   * cosze  * cos2Q
  334.           +  0.003342                   * cos2ze * cos2Q;
  335.  
  336.    w1pert = (0.007192 -  0.003147 * nu2)             * sinV
  337.           - poly(0.020428,0.000675,-0.000197,0.,nu2) * cosV
  338.           + (0.007269 +  0.000672 * nu2) * sinze  * sinQ
  339.           -  0.004344                             * sinQ
  340.           +  0.034036                    * cosze  * sinQ
  341.           +  0.005614                    * cos2ze * sinQ
  342.           +  0.002964                    * cos3ze * sinQ
  343.           +  0.037761                    * sinze  * cosQ
  344.           +  0.006158                    * sin2ze * cosQ
  345.           -  0.006603                    * cosze  * cosQ
  346.           -  0.005356                    * sinze  * sin2Q
  347.           +  0.002722                    * sin2ze * sin2Q
  348.           +  0.004483                    * cosze  * sin2Q
  349.           -  0.002642                    * cos2ze * sin2Q
  350.           +  0.004403                    * sinze  * cos2Q
  351.           -  0.002536                    * sin2ze * cos2Q
  352.           +  0.005547                    * cosze  * cos2Q
  353.           -  0.002689                    * cos2ze * cos2Q;
  354.  
  355.    l  += l1pert;
  356.    w1 += w1pert;
  357.  
  358.    M0  = M5 + l1pert - (w1pert / e);
  359.  
  360.    e += poly(0.0003606, 0.0000130,-0.0000043,0.,nu2) * sinV
  361.       + poly(0.0001289,-0.0000580,0.        ,0.,nu2) * cosV
  362.       -  0.0006764                    * sinze  * sinQ
  363.       -  0.0001110                    * sin2ze * sinQ
  364.       -  0.0000224                    * sin3ze * sinQ
  365.       -  0.0000204                             * sinQ
  366.       + (0.0001284 + 0.0000116 * nu2) * cosze  * sinQ
  367.       +  0.0000188                    * cos2ze * sinQ
  368.       + (0.0001460 + 0.0000130 * nu2) * sinze  * cosQ
  369.       +  0.0000224                    * sin2ze * cosQ
  370.       -  0.0000817                             * cosQ
  371.       +  0.0006074                    * cosze  * cosQ
  372.       +  0.0000992                    * cos2ze * cosQ
  373.       +  0.0000508                    * cos3ze * cosQ
  374.       +  0.0000230                    * cos4ze * cosQ
  375.       +  0.0000108                    * cos(5. * ze) * cosQ
  376.       - (0.0000956 + 0.0000073 * nu2) * sinze  * sin2Q
  377.       +  0.0000448                    * sin2ze * sin2Q
  378.       +  0.0000137                    * sin3ze * sin2Q
  379.       - (0.0000997 - 0.0000108 * nu2) * cosze  * sin2Q
  380.       +  0.0000480                    * cos2ze * sin2Q
  381.       +  0.0000148                    * cos3ze * sin2Q
  382.       - (0.0000956 - 0.0000099 * nu2) * sinze  * cos2Q
  383.       +  0.0000490                    * sin2ze * cos2Q
  384.       +  0.0000158                    * sin3ze * cos2Q
  385.       +  0.0000179                             * cos2Q
  386.       + (0.0001024 + 0.0000075 * nu2) * cosze  * cos2Q
  387.       -  0.0000437                    * cos2ze * cos2Q
  388.       -  0.0000132                    * cos3ze * cos2Q;
  389.  
  390.    a += -0.000263          * cosV
  391.       +  0.000205 * cosze
  392.       +  0.000693 * cos2ze
  393.       +  0.000312 * cos3ze
  394.       +  0.000147 * cos4ze
  395.       +  0.000299 * sinze  * sinQ
  396.       +  0.000181 * cos2ze * sinQ
  397.       +  0.000204 * sin2ze * cosQ
  398.       +  0.000111 * sin3ze * cosQ
  399.       -  0.000337 * cosze  * cosQ
  400.       -  0.000111 * cos2ze * cosQ;
  401.  
  402.    ECC = kepler(e,M0);
  403.    nu  = truean(e,ECC);
  404.    r   = a * (1. - e * cos(ECC * RAD));
  405.    u   = l + nu - M0 - w2;
  406.    ll  = longi(w2,i,u);
  407.    b   = lati(u,i);
  408.  
  409. /*******************************************************************************
  410. ** Tansformation of coordinates on Jupiter and output.                        **
  411. *******************************************************************************/
  412.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGJUP,"PJ","Jupiter");
  413.  
  414. /*******************************************************************************
  415. ** Now start Saturn.                                                          **
  416. *******************************************************************************/
  417.    l = range(poly(266.564377,1223.509884,0.0003245,-0.0000058,T0));
  418.  
  419.    a = 9.554747;
  420.  
  421.    e = poly(0.05589232,-0.00034550,-0.000000728,0.00000000074,T0);
  422.    i = poly(2.49251900,-0.00391890,-0.000015490,0.00000004000,T0);
  423.  
  424.    w1 = range(poly(338.307800,1.0852207, 0.00097854, 0.00000992,T0));
  425.    w2 = range(poly(112.790414,0.8731951,-0.00015218,-0.00000531,T0));
  426.  
  427. /*******************************************************************************
  428. ** Now Saturn's perturbations.                                                **
  429. *******************************************************************************/
  430.    l1pert = poly(-0.814181,0.018150, 0.016714,0.,nu2) * sinV
  431.           + poly(-0.010497,0.160906,-0.004100,0.,nu2) * cosV
  432.           +   0.007581                                * sin2V
  433.           -   0.007986                                * sin(W)
  434.           -   0.148811                   * sinze
  435.           -   0.040786                   * sin2ze
  436.           -   0.015208                   * sin3ze
  437.           -   0.006339                   * sin4ze
  438.           -   0.006244                            * sinQ
  439.           + ( 0.008931 + 0.002728 * nu2) * sinze  * sinQ
  440.           -   0.016500                   * sin2ze * sinQ
  441.           -   0.005775                   * sin3ze * sinQ
  442.           + ( 0.081344 + 0.003206 * nu2) * cosze  * sinQ
  443.           +   0.015019                   * cos2ze * sinQ
  444.           + ( 0.085581 + 0.002494 * nu2) * sinze  * cosQ
  445.           + ( 0.025328 - 0.003117 * nu2) * cosze  * cosQ
  446.           +   0.014394                   * cos2ze * cosQ
  447.           +   0.006319                   * cos3ze * cosQ
  448.           +   0.006369                   * sinze  * sin2Q
  449.           +   0.009156                   * sin2ze * sin2Q
  450.           +   0.007525                   * sin3ps * sin2Q
  451.           -   0.005236                   * cosze  * cos2Q
  452.           -   0.007736                   * cos2ze * cos2Q
  453.           -   0.007528                   * cos3ps * cos2Q;
  454.  
  455.    w1pert = poly(0.077108, 0.007186,-0.001533,0.,nu2) * sinV
  456.           + poly(0.045803,-0.014766,-0.000536,0.,nu2) * cosV
  457.           -   0.007075                   * sinze
  458.           -   0.075825                   * sinze  * sinQ
  459.           -   0.024839                   * sin2ze * sinQ
  460.           -   0.008631                   * sin3ze * sinQ
  461.           -   0.072586                            * cosQ
  462.           -   0.150383                   * cosze  * cosQ
  463.           +   0.026897                   * cos2ze * cosQ
  464.           +   0.010053                   * cos3ze * cosQ
  465.           - ( 0.013597 + 0.001719 * nu2) * sinze  * sin2Q
  466.           + (-0.007742 + 0.001517 * nu2) * cosze  * sin2Q
  467.           + ( 0.013586 - 0.001375 * nu2) * cos2ze * sin2Q
  468.           + (-0.013667 + 0.001239 * nu2) * sinze  * cos2Q
  469.           +   0.011981                   * sin2ze * cos2Q
  470.           + ( 0.014861 + 0.001136 * nu2) * cosze  * cos2Q
  471.           - ( 0.013064 + 0.001628 * nu2) * cos2ze * cos2Q;
  472.  
  473.    l  += l1pert;
  474.    w1 += w1pert;
  475.  
  476.    M0  = M6 + l1pert - (w1pert / e);
  477.  
  478.    e += poly(-0.0007927,0.0002548, 0.0000091,0.,nu2) * sinV
  479.       + poly( 0.0013381,0.0001226,-0.0000253,0.,nu2) * cosV
  480.       + ( 0.0000248 -  0.0000121 * nu2)              * sin2V
  481.       - ( 0.0000305 +  0.0000091 * nu2)              * cos2V
  482.       +   0.0000412                    * sin2ze
  483.       +   0.0012415                             * sinQ
  484.       + ( 0.0000390 - 0.0000617 * nu2) * sinze  * sinQ
  485.       + ( 0.0000165 - 0.0000204 * nu2) * sin2ze * sinQ
  486.       +   0.0026599                    * cosze  * sinQ
  487.       -   0.0004687                    * cos2ze * sinQ
  488.       -   0.0001870                    * cos3ze * sinQ
  489.       -   0.0000821                    * cos4ze * sinQ
  490.       -   0.0000377                    * cos(5. * ze)  * sinQ
  491.       +   0.0000497                    * cos2ps * sinQ
  492.       + ( 0.0000163 - 0.0000611 * nu2)          * cosQ
  493.       -   0.0012696                    * sinze  * cosQ
  494.       -   0.0004200                    * sin2ze * cosQ
  495.       -   0.0001503                    * sin3ze * cosQ
  496.       -   0.0000619                    * sin4ze * cosQ
  497.       -   0.0000268                    * sin(5. * ze)  * cosQ
  498.       - ( 0.0000282 + 0.0001306 * nu2) * cosze  * cosQ
  499.       + (-0.0000086 + 0.0000230 * nu2) * cos2ze * cosQ
  500.       +   0.0000461                    * sin2ps * cosQ
  501.       -   0.0000350                             * sin2Q
  502.       + ( 0.0002211 - 0.0000286 * nu2) * sinze  * sin2Q
  503.       -   0.0002208                    * sin2ze * sin2Q
  504.       -   0.0000568                    * sin3ze * sin2Q
  505.       -   0.0000346                    * sin4ze * sin2Q
  506.       - ( 0.0002780 + 0.0000222 * nu2) * cosze  * sin2Q
  507.       + ( 0.0002022 + 0.0000263 * nu2) * cos2ze * sin2Q
  508.       +   0.0000248                    * cos3ze * sin2Q
  509.       +   0.0000242                    * sin3ps * sin2Q
  510.       +   0.0000467                    * cos3ps * sin2Q
  511.       -   0.0000490                             * cos2Q
  512.       - ( 0.0002842 + 0.0000279 * nu2) * sinze  * cos2Q
  513.       + ( 0.0000128 + 0.0000226 * nu2) * sin2ze * cos2Q
  514.       +   0.0000224                    * sin3ze * cos2Q
  515.       + (-0.0001594 + 0.0000282 * nu2) * cosze  * cos2Q
  516.       + ( 0.0002162 - 0.0000207 * nu2) * cos2ze * cos2Q
  517.       +   0.0000561                    * cos3ze * cos2Q
  518.       +   0.0000343                    * cos4ze * cos2Q
  519.       +   0.0000469                    * sin3ps * cos2Q
  520.       -   0.0000242                    * cos3ps * cos2Q
  521.       -   0.0000205                    * sinze  * sin3Q
  522.       +   0.0000262                    * sin3ze * sin3Q
  523.       +   0.0000208                    * cosze  * cos3Q
  524.       -   0.0000271                    * cos3ze * cos3Q
  525.       -   0.0000382                    * cos3ze * sin4Q
  526.       -   0.0000376                    * sin3ze * cos4Q;
  527.  
  528.    a +=  0.000572                            * sinV
  529.       -  0.001590                   * sin2ze * cosQ
  530.       +  0.002933                            * cosV
  531.       -  0.000647                   * sin3ze * cosQ
  532.       +  0.033629                   * cosze
  533.       -  0.000344                   * sin4ze * cosQ
  534.       -  0.003081                   * cos2ze
  535.       +  0.002885                   * cosze  * cosQ
  536.       -  0.001423                   * cos3ze
  537.       + (0.002172 + 0.000102 * nu2) * cos2ze * cosQ
  538.       -  0.000671                   * cos4ze
  539.       +  0.000296                   * cos3ze * cosQ
  540.       -  0.000320                   * cos(5. * ze)
  541.       -  0.000267                   * sin2ze * sin2Q
  542.       +  0.001098                            * sinQ
  543.       -  0.000778                   * cosze  * sin2Q
  544.       -  0.002812                   * sinze  * sinQ
  545.       +  0.000495                   * cos2ze * sin2Q
  546.       +  0.000688                   * sin2ze * sinQ
  547.       +  0.000250                   * cos3ze * sin2Q
  548.       -  0.000393                   * sin3ze * sinQ
  549.       -  0.000228                   * sin4ze * sinQ
  550.       +  0.002138                   * cosze  * sinQ
  551.       -  0.000999                   * cos2ze * sinQ
  552.       -  0.000642                   * cos3ze * sinQ
  553.       -  0.000325                   * cos4ze * sinQ
  554.       -  0.000890                            * cosQ
  555.       +  0.002206                   * sinze  * cosQ
  556.       -  0.000856                   * sinze  * cos2Q
  557.       +  0.000441                   * sin2ze * cos2Q
  558.       +  0.000296                   * cos2ze * cos2Q
  559.       +  0.000211                   * cos3ze * cos2Q
  560.       -  0.000427                   * sinze  * sin3Q
  561.       +  0.000398                   * sin3ze * sin3Q
  562.       +  0.000344                   * cosze  * cos3Q
  563.       -  0.000427                   * cos3ze * cos3Q;
  564.  
  565.    ECC = kepler(e,M0);
  566.    nu  = truean(e,ECC);
  567.    r   = a * (1. - e * cos(ECC * RAD));
  568.    u   = l + nu - M0 - w2;
  569.    ll  = longi(w2,i,u);
  570.  
  571.    b = lati(u,i)
  572.      + 0.000747 * cosze  * sinQ
  573.      + 0.001069 * cosze  * cosQ
  574.      + 0.002108 * sin2ze * sin2Q
  575.      + 0.001261 * cos2ze * sin2Q
  576.      + 0.001236 * sin2ze * cos2Q
  577.      - 0.002075 * cos2ze * cos2Q;
  578.  
  579. /*******************************************************************************
  580. ** Tansformation of coordinates on Saturn and output.                         **
  581. *******************************************************************************/
  582.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGSAT,"Ps","Saturn");
  583.  
  584. /*******************************************************************************
  585. ** Now Start on Uranus.                                                       **
  586. *******************************************************************************/
  587.    l = range(poly(244.197470,429.863546,0.0003160,-0.00000060,T0));
  588.  
  589.    a = 19.21814;
  590.  
  591.    e = poly(0.0463444,-0.00002658,0.000000077,0.,T0);
  592.    i = poly(0.7724640, 0.00062530,0.000039500,0.,T0);
  593.  
  594.    w1 = range(poly(98.071581,0.9857650,-0.0010745,-0.00000061,T0));
  595.    w2 = range(poly(73.477111,0.4986678, 0.0013117, 0.00000000,T0));
  596.  
  597.    M7 = range(poly(72.64878,428.37911,0.000079,0.,T0));
  598.  
  599. /*******************************************************************************
  600. ** Now perturbation corrections for Uranus.                                   **
  601. *******************************************************************************/
  602.    G = (83.76922 + 218.4901 * T0) * RAD;
  603.  
  604.    H   = range((2. * G - S) / RAD) * RAD;
  605.    ze  = range((     S - P) / RAD) * RAD;
  606.    eta = range((     S - Q) / RAD) * RAD;
  607.    th  = range((     G - S) / RAD) * RAD;
  608.  
  609.    G = range(G / RAD) * RAD;
  610.  
  611. /*******************************************************************************
  612. ** Addition theorems.                                                         **
  613. *******************************************************************************/
  614.    sinze = sin(ze); cosze = cos(ze);
  615.  
  616.    sinH  = sin(H);           cosH  = cos(H);
  617.    sin2H = 2. * sinH * cosH; cos2H = cosH * cosH - sinH * sinH;
  618.  
  619.    sinG = sin(G); cosG = cos(G);
  620.  
  621. /******************************************************************************/
  622.  
  623.    l1pert = (0.864319 - 0.001583 * nu2) * sinH
  624.           + (0.082222 - 0.006833 * nu2) * cosH
  625.           +  0.036017                   * sin2H
  626.           -  0.003019                   * cos2H
  627.           +  0.008122                   * sin(W);
  628.  
  629.    w1pert = 0.120303                    * sinH
  630.           + (0.019472 - 0.000947 * nu2) * cosH
  631.           + 0.006197                    * sin2H;
  632.  
  633.    l  += l1pert;
  634.    w1 += w1pert;
  635.  
  636.    M0  = M7 + l1pert - (w1pert / e);
  637.  
  638.    e += (-0.0003349 + 0.0000163 * nu2) * sinH
  639.       +   0.0020981                    * cosH
  640.       +   0.0001311                    * cosH;
  641.  
  642.    a -= 0.003825 * cosH;
  643.  
  644.    ECC = kepler(e,M0);
  645.    nu  = truean(e,ECC);
  646.    r   = a * (1. - e * cos(ECC * RAD));
  647.    u   = l + nu - M0 - w2;
  648.    ll  = longi(w2,i,u);
  649.    b   = lati(u,i);
  650.  
  651.   ll += (0.010122 -  0.000988 * nu2)              * sin(     eta +      S)
  652.      + poly(-0.038581, 0.002031,-0.001910,0.,nu2) * cos(     eta +      S)
  653.      + poly( 0.034964,-0.001038, 0.000868,0.,nu2) * cos(     eta + 2. * S)
  654.      +   0.005594                                 * sin(3. * th  +      S)
  655.      -   0.014808                                 * sinze
  656.      -   0.005794                                 * sin(     eta)
  657.      +   0.002347                                 * cos(     eta)
  658.      +   0.009872                                 * sin(     th)
  659.      +   0.008803                                 * sin(2. * th)
  660.      -   0.004308                                 * sin(3. * th);
  661.  
  662.    b += (0.000458 *  sin(     eta)
  663.      -   0.000642 *  cos(     eta)
  664.      -   0.000517 *  cos(4. * th)) * sinS
  665.      -  (0.000347 *  sin(     eta)
  666.      +   0.000853 *  cos(     eta)
  667.      +   0.000517 *  sin(4. * eta)) * cosS
  668.      +   0.000403 * (cos(2. * th)   * sin2S
  669.                   +  sin(2. * th)   * cos2S);
  670.  
  671.    r += -0.025948
  672.      +   0.004985 * cosze
  673.      -   0.001230                 * cosS
  674.      +   0.003354 * cos(     eta)
  675.      +  (0.005795                 * cosS
  676.      -   0.001165                 * sinS
  677.      +   0.001388 * sin(     eta) * cos2S)
  678.      +  (0.001351                 * cosS
  679.      +   0.005702                 * sinS
  680.      +   0.001388 * cos(     eta) * sin2S)
  681.      +   0.000904 * cos(2. * th)
  682.      +   0.000894 * (cos(    th) - cos(3. * th));
  683.  
  684. /*******************************************************************************
  685. ** Tansformation of coordinates on Uranus and output.                         **
  686. *******************************************************************************/
  687.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGURA,"PU","Uranus");
  688.  
  689. /*******************************************************************************
  690. ** Now start Neptune.                                                         **
  691. *******************************************************************************/
  692.    l = range(poly(84.457994,219.885914,0.0003205,-0.00000060,T0));
  693.  
  694.    a = 30.10957;
  695.  
  696.    e = poly(0.00899704, 0.000006330,-0.000000002,0.,T0);
  697.    i = poly(1.77924200,-0.009543600,-0.000009100,0.,T0);
  698.  
  699.    w1 = range(poly(276.045975,0.3256394,0.00014095, 0.000004113,T0));
  700.    w2 = range(poly(130.681389,1.0989350,0.00024987,-0.000004718,T0));
  701.  
  702.    M8 = poly(37.73063,218.46134,-0.000070,0.,T0);
  703.  
  704. /*******************************************************************************
  705. ** Now perturbation corrections for neptune.                                  **
  706. *******************************************************************************/
  707.    ze  = range((G - P) / RAD) * RAD;
  708.    eta = range((G - Q) / RAD) * RAD;
  709.    th  = range((G - S) / RAD) * RAD;
  710.  
  711.    sinze = sin(ze); cosze = cos(ze);
  712.  
  713.    l1pert = (-0.589833 + 0.001089 * nu2) * sinH
  714.           + (-0.056094 + 0.004658 * nu2) * cosH
  715.           -   0.024286                   * sin2H;
  716.  
  717.    w1pert = 0.024039 * sinH
  718.           - 0.025303 * cosH
  719.           + 0.006206 * sin2H
  720.           - 0.005992 * cos2H;
  721.  
  722.    l  += l1pert;
  723.    w1 += w1pert;
  724.  
  725.    M0  = M8 + l1pert - (w1pert / e);
  726.  
  727.    e += 0.0004389 * sinH
  728.       + 0.0004262 * cosH
  729.       + 0.0001129 * sin2H
  730.       + 0.0001089 * cos2H;
  731.  
  732.    a += -0.000817 * sinH
  733.      +   0.008189 * cosH
  734.      +   0.000781 * cos2H;
  735.  
  736.    ECC = kepler(e,M0);
  737.    nu  = truean(e,ECC);
  738.    r   = a * (1. - e * cos(ECC * RAD));
  739.    u   = l + nu - M0 - w2;
  740.    ll  = longi(w2,i,u);
  741.    b   = lati(u,i);
  742.  
  743.    ll += -0.009556 * sinze
  744.       -   0.005178 * sin(     eta)
  745.       +   0.002572 * sin(2. * th)
  746.       -   0.002972 * cos(2. * th) * sinG
  747.       -   0.002833 * sin(2. * th) * cosG;
  748.  
  749.    b += 0.000336 * cos(2. * th) * sinG
  750.      +  0.000364 * sin(2. * th) * cosG;
  751.  
  752.    r += -0.040596
  753.      +   0.004992 * cosze
  754.      +   0.002744 * cos(     eta)
  755.      +   0.002044 * cos(     th)
  756.      +   0.001051 * cos(2. * th);
  757.  
  758. /*******************************************************************************
  759. ** Tansformation of coordinates on Neptune and output.                        **
  760. *******************************************************************************/
  761. #ifdef EUTSCH
  762.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGNEP,"PN","Neptun");
  763. #else
  764.    helio_trans(r,b,ll,Stheta,Sr,epli,MAGNEP,"PN","Neptune");
  765. #endif
  766.  
  767.    return(epli);
  768.    }
  769.