开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > OpenGL > 查看源码
customrotation.cpp
资源名称：celestia-1.6.0.tar.gz　[点击查看]
上传用户：center1979
上传日期：2022-07-26
资源大小：50633k
文件大小：33k
源码类别：

OpenGL

开发平台：

Visual C++

customrotation.cpp：源码内容
  1. // customrotation.cpp
  2. //
  3. // Custom rotation models for Solar System bodies.
  4. //
  5. // Copyright (C) 2008, the Celestia Development Team
  6. // Initial version by Chris Laurel, claurel@gmail.com
  7. //
  8. // This program is free software; you can redistribute it and/or
  9. // modify it under the terms of the GNU General Public License
  10. // as published by the Free Software Foundation; either version 2
  11. // of the License, or (at your option) any later version.
  13. #include <map>
  14. #include <string>
  15. #include <celengine/customrotation.h>
  16. #include <celengine/rotation.h>
  17. #include <celengine/astro.h>
  18. #include <celengine/precession.h>
  20. using namespace std;
  23. static map<string, RotationModel*> CustomRotationModels;
  24. static bool CustomRotationModelsInitialized = false;
  27. // Clamp secular terms in IAU rotation models to this number of centuries
  28. // from J2000. Extrapolating much further can lead to ridiculous results,
  29. // such as planets 'tipping over' Periodic terms are not clamped; their
  30. // validity over long time ranges is questionable, but extrapolating them
  31. // doesn't produce obviously absurd results.
  32. static const double IAU_SECULAR_TERM_VALID_CENTURIES = 50.0;
  34. // The P03 long period precession theory for Earth is valid for a one
  35. // million year time span centered on J2000. For dates outside far outside
  36. // that range, the polynomial terms produce absurd results.
  37. static const double P03LP_VALID_CENTURIES = 5000.0;
  39. /*! Base class for IAU rotation models. All IAU rotation models are in the
  40.  *  J2000.0 Earth equatorial frame.
  41.  */
  42. class IAURotationModel : public CachingRotationModel
  43. {
  44. public:
  45.     IAURotationModel(double _period) :
  46.         period(_period), flipped(false)
  47.     {
  48.     }
  49.     
  50.     virtual ~IAURotationModel() {}
  51.     
  52.     bool isPeriodic() const { return true; }
  53.     double getPeriod() const { return period; }
  54.     
  55.     virtual Quatd computeSpin(double t) const
  56.     {
  57.         // Time argument of IAU rotation models is actually day since J2000.0 TT, but
  58.         // Celestia uses TDB. The difference should be so minute as to be irrelevant.
  59.         t = t - astro::J2000;
  60.         if (flipped)
  61.             return Quatd::yrotation( degToRad(180.0 + meridian(t)));
  62.         else
  63.             return Quatd::yrotation(-degToRad(180.0 + meridian(t)));
  64.     }
  65.     
  66.     virtual Quatd computeEquatorOrientation(double t) const
  67.     {
  68.         double poleRA = 0.0;
  69.         double poleDec = 0.0;
  70.         
  71.         t = t - astro::J2000;
  72.         pole(t, poleRA, poleDec);
  73.         double node = poleRA + 90.0;
  74.         double inclination = 90.0 - poleDec;
  75.         
  76.         if (flipped)
  77.             return Quatd::xrotation(PI) * Quatd::xrotation(degToRad(-inclination)) * Quatd::yrotation(degToRad(-node));
  78.         else
  79.             return Quatd::xrotation(degToRad(-inclination)) * Quatd::yrotation(degToRad(-node));
  80.     }
  81.     
  82.     // Return the RA and declination (in degrees) of the rotation axis
  83.     virtual void pole(double t, double& ra, double& dec) const = 0;
  84.     
  85.     virtual double meridian(double t) const = 0;
  87. protected:
  88.     static void clamp_centuries(double& T)
  89.     {
  90.         if (T < -IAU_SECULAR_TERM_VALID_CENTURIES)
  91.             T = -IAU_SECULAR_TERM_VALID_CENTURIES;
  92.         else if (T > IAU_SECULAR_TERM_VALID_CENTURIES)
  93.             T = IAU_SECULAR_TERM_VALID_CENTURIES;
  94.     }
  96.     void setFlipped(bool _flipped)
  97.     {
  98.         flipped = _flipped;
  99.     }
  100.     
  101. private:
  102.     double period;
  103.     bool flipped;
  104. };
  108. /******* Earth rotation model *******/
  110. class EarthRotationModel : public CachingRotationModel
  111. {
  112. public:
  113.     EarthRotationModel()
  114.     {
  115.     }
  116.     
  117.     ~EarthRotationModel()
  118.     {
  119.     }
  120.     
  121.     Quatd computeSpin(double tjd) const
  122.     {
  123.         // TODO: Use a more accurate model for sidereal time
  124.         double t = tjd - astro::J2000;
  125.         double theta = 2 * PI * (t * 24.0 / 23.9344694 - 259.853 / 360.0);
  127.         return Quatd::yrotation(-theta);
  128.     }
  129.     
  130.     Quatd computeEquatorOrientation(double tjd) const
  131.     {
  132.         double T = (tjd - astro::J2000) / 36525.0;
  133.         
  134.         // Clamp T to the valid time range of the precession theory.
  135.         if (T < -P03LP_VALID_CENTURIES)
  136.             T = -P03LP_VALID_CENTURIES;
  137.         else if (T > P03LP_VALID_CENTURIES)
  138.             T = P03LP_VALID_CENTURIES;
  139.         
  140.         astro::PrecessionAngles prec = astro::PrecObliquity_P03LP(T);
  141.         astro::EclipticPole pole = astro::EclipticPrecession_P03LP(T);
  142.         
  143.         double obliquity = degToRad(prec.epsA / 3600);
  144.         double precession = degToRad(prec.pA / 3600);
  145.         
  146.         // Calculate the angles pi and Pi from the ecliptic pole coordinates
  147.         // P and Q:
  148.         //   P = sin(pi)*sin(Pi)
  149.         //   Q = sin(pi)*cos(Pi)
  150.         double P = pole.PA * 2.0 * PI / 1296000;
  151.         double Q = pole.QA * 2.0 * PI / 1296000;
  152.         double piA = asin(sqrt(P * P + Q * Q));
  153.         double PiA = atan2(P, Q);
  155.         // Calculate the rotation from the J2000 ecliptic to the ecliptic
  156.         // of date.
  157.         Quatd RPi = Quatd::zrotation(PiA);
  158.         Quatd rpi = Quatd::xrotation(piA);
  159.         Quatd eclRotation = ~RPi * rpi * RPi;
  161.         Quatd q = Quatd::xrotation(obliquity) * Quatd::zrotation(-precession) * ~eclRotation;
  162.         
  163.         // convert to Celestia's coordinate system
  164.         return Quatd::xrotation(PI / 2.0) * q * Quatd::xrotation(-PI / 2.0);
  165.     }
  166.     
  167.     double getPeriod() const
  168.     {
  169.         return 23.9344694 / 24.0;
  170.     }
  172.     bool isPeriodic() const
  173.     {
  174.         return true;
  175.     }
  176. };
  180. /******* IAU rotation models for the planets *******/
  183. /*! IAUPrecessingRotationModel is a rotation model with uniform rotation about
  184.  *  a pole that precesses linearly in RA and declination.
  185.  */
  186. class IAUPrecessingRotationModel : public IAURotationModel
  187. {
  188. public:
  189.     
  190.     /*! rotationRate is in degrees per Julian day
  191.      *  pole precession is in degrees per Julian century
  192.      */
  193.     IAUPrecessingRotationModel(double _poleRA,
  194.                                double _poleRARate,
  195.                                double _poleDec,
  196.                                double _poleDecRate,
  197.                                double _meridianAtEpoch,
  198.                                double _rotationRate) :
  199.     IAURotationModel(std::abs(360.0 / _rotationRate)),
  200.         poleRA(_poleRA),
  201.         poleRARate(_poleRARate),
  202.         poleDec(_poleDec),
  203.         poleDecRate(_poleDecRate),
  204.         meridianAtEpoch(_meridianAtEpoch),
  205.         rotationRate(_rotationRate)
  206.     {
  207.         if (rotationRate < 0.0)
  208.             setFlipped(true);
  209.     }
  210.     
  211.     void pole(double d, double& ra, double &dec) const
  212.     {
  213.         double T = d / 36525.0;
  214.         clamp_centuries(T);
  215.         ra = poleRA + poleRARate * T;
  216.         dec = poleDec + poleDecRate * T;
  217.     }
  218.     
  219.     double meridian(double d) const
  220.     {
  221.         return meridianAtEpoch + rotationRate * d;
  222.     }
  223.     
  224. private:
  225.     double poleRA;
  226.     double poleRARate;
  227.     double poleDec;
  228.     double poleDecRate;
  229.     double meridianAtEpoch;
  230.     double rotationRate;
  231. };
  234. class IAUNeptuneRotationModel : public IAURotationModel
  235. {
  236. public:
  237.     IAUNeptuneRotationModel() : IAURotationModel(360.0 / 536.3128492) {}
  238.     
  239.     void pole(double d, double& ra, double &dec) const
  240.     {
  241.         double T = d / 36525.0;
  242.         double N = degToRad(357.85 + 52.316 * T);
  243.         ra = 299.36 + 0.70 * sin(N);
  244.         dec = 43.46 - 0.51 * cos(N);
  245.     }
  246.     
  247.     double meridian(double d) const
  248.     {
  249.         double T = d / 36525.0;
  250.         double N = degToRad(357.85 + 52.316 * T);
  251.         return 253.18 + 536.3128492 * d - 0.48 * sin(N);
  252.     }    
  253. };
  256. /*! IAU rotation model for the Moon.
  257.  *  From the IAU/IAG Working Group on Cartographic Coordinates and Rotational Elements:
  258.  *  http://astrogeology.usgs.gov/Projects/WGCCRE/constants/iau2000_table2.html
  259.  */
  260. class IAULunarRotationModel : public IAURotationModel
  261. {
  262. public:
  263.     IAULunarRotationModel() : IAURotationModel(360.0 / 13.17635815) {}
  264.     
  265.     void calcArgs(double d, double E[14]) const
  266.     {
  267.         E[1]  = degToRad(125.045 -  0.0529921 * d);
  268.         E[2]  = degToRad(250.089 -  0.1059842 * d);
  269.         E[3]  = degToRad(260.008 + 13.012009 * d);
  270.         E[4]  = degToRad(176.625 + 13.3407154 * d);
  271.         E[5]  = degToRad(357.529 +  0.9856993 * d);
  272.         E[6]  = degToRad(311.589 + 26.4057084 * d);
  273.         E[7]  = degToRad(134.963 + 13.0649930 * d);
  274.         E[8]  = degToRad(276.617 +  0.3287146 * d);
  275.         E[9]  = degToRad( 34.226 +  1.7484877 * d);
  276.         E[10] = degToRad( 15.134 -  0.1589763 * d);
  277.         E[11] = degToRad(119.743 +  0.0036096 * d);
  278.         E[12] = degToRad(239.961 +  0.1643573 * d);
  279.         E[13] = degToRad( 25.053 + 12.9590088 * d);
  280.     }
  281.     
  282.     void pole(double d, double& ra, double& dec) const
  283.     {
  284.         double T = d / 36525.0;
  285.         clamp_centuries(T);
  287.         double E[14];
  288.         calcArgs(d, E);
  289.         
  290.         ra = 269.9949
  291.             + 0.0013*T
  292.             - 3.8787 * sin(E[1])
  293.             - 0.1204 * sin(E[2])
  294.             + 0.0700 * sin(E[3])
  295.             - 0.0172 * sin(E[4])
  296.             + 0.0072 * sin(E[6])
  297.             - 0.0052 * sin(E[10])
  298.             + 0.0043 * sin(E[13]);
  299.             
  300.         dec = 66.5392
  301.             + 0.0130 * T
  302.             + 1.5419 * cos(E[1])
  303.             + 0.0239 * cos(E[2])
  304.             - 0.0278 * cos(E[3])
  305.             + 0.0068 * cos(E[4])
  306.             - 0.0029 * cos(E[6])
  307.             + 0.0009 * cos(E[7])
  308.             + 0.0008 * cos(E[10])
  309.             - 0.0009 * cos(E[13]);
  310.             
  311.             
  312.     }
  313.     
  314.     double meridian(double d) const
  315.     {
  316.         double E[14];
  317.         calcArgs(d, E);
  319.         // d^2 represents slowing of lunar rotation as the Moon recedes
  320.         // from the Earth. This may need to be clamped at some very large
  321.         // time range (1 Gy?)
  323.         return (38.3213
  324.                 + 13.17635815 * d
  325.                 - 1.4e-12 * d * d
  326.                 + 3.5610 * sin(E[1])
  327.                 + 0.1208 * sin(E[2])
  328.                 - 0.0642 * sin(E[3])
  329.                 + 0.0158 * sin(E[4])
  330.                 + 0.0252 * sin(E[5])
  331.                 - 0.0066 * sin(E[6])
  332.                 - 0.0047 * sin(E[7])
  333.                 - 0.0046 * sin(E[8])
  334.                 + 0.0028 * sin(E[9])
  335.                 + 0.0052 * sin(E[10])
  336.                 + 0.0040 * sin(E[11])
  337.                 + 0.0019 * sin(E[12])
  338.                 - 0.0044 * sin(E[13]));
  339.     }
  340. };
  343. // Rotations of Martian, Jovian, and Uranian satellites from IAU/IAG Working group 
  344. // on Cartographic Coordinates and Rotational Elements (Corrected for known errata
  345. // as of 17 Feb 2006)
  346. // See: http://astrogeology.usgs.gov/Projects/WGCCRE/constants/iau2000_table2.html
  348. class IAUPhobosRotationModel : public IAURotationModel
  349. {
  350. public:
  351.     IAUPhobosRotationModel() : IAURotationModel(360.0 / 1128.8445850) {}
  352.     
  353.     void pole(double t, double& ra, double& dec) const
  354.     {
  355.         double T = t / 36525.0;
  356.         double M1 = degToRad(169.51 - 0.04357640 * t);
  357.         clamp_centuries(T);
  358.         ra  = 317.68 - 0.108 * T + 1.79 * sin(M1);
  359.         dec =  52.90 - 0.061 * T - 1.08 * cos(M1);
  360.     }
  361.     
  362.     double meridian(double t) const
  363.     {
  364.         // Note: negative coefficient of T^2 term for meridian angle indicates faster
  365.         // rotation as Phobos's orbit evolves inward toward Mars
  366.         double T = t / 36525.0;
  367.         double M1 = degToRad(169.51 - 0.04357640 * t);
  368.         double M2 = degToRad(192.93 + 1128.4096700 * t + 8.864 * T * T);
  369.         return 35.06 + 1128.8445850 * t + 8.864 * T * T - 1.42 * sin(M1) - 0.78 * sin(M2);
  370.     }
  371. };
  374. class IAUDeimosRotationModel : public IAURotationModel
  375. {
  376. public:
  377.     IAUDeimosRotationModel() : IAURotationModel(360.0 / 285.1618970) {}
  378.     
  379.     void pole(double t, double& ra, double& dec) const
  380.     {
  381.         double T = t / 36525.0;
  382.         double M3 = degToRad(53.47 - 0.0181510 * t);
  383.         clamp_centuries(T);
  384.         ra  = 316.65 - 0.108 * T + 2.98 * sin(M3);
  385.         dec =  53.52 - 0.061 * T - 1.78 * cos(M3);
  386.     }
  387.     
  388.     double meridian(double t) const
  389.     {
  390.         // Note: positive coefficient of T^2 term for meridian angle indicates slowing
  391.         // rotation as Deimos's orbit evolves outward from Mars
  392.         double T = t / 36525.0;
  393.         double M3 = degToRad(53.47 - 0.0181510 * t);
  394.         return 79.41 + 285.1618970 * t + 0.520 * T * T - 2.58 * sin(M3) + 0.19 * cos(M3);
  395.     }
  396. };
  399. /****** Satellites of Jupiter ******/
  400.         
  401. class IAUAmaltheaRotationModel : public IAURotationModel
  402. {
  403. public:
  404.     IAUAmaltheaRotationModel() : IAURotationModel(360.0 / 722.6314560) {}
  405.     
  406.     void pole(double t, double& ra, double& dec) const
  407.     {
  408.         double T = t / 36525.0;
  409.         double J1 = degToRad(73.32 + 91472.9 * T);
  410.         clamp_centuries(T);
  411.         ra = 268.05 - 0.009 * T - 0.84 * sin(J1) + 0.01 * sin(2.0 * J1);
  412.         dec = 64.49 + 0.003 * T - 0.36 * cos(J1);
  413.     }
  414.     
  415.     double meridian(double t) const
  416.     {
  417.         double T = t / 36525.0;
  418.         double J1 = degToRad(73.32 + 91472.9 * T);
  419.         return 231.67 + 722.6314560 * t + 0.76 * sin(J1) - 0.01 * sin(2.0 * J1);
  420.     }
  421. };
  424. class IAUThebeRotationModel : public IAURotationModel
  425. {
  426. public:
  427.     IAUThebeRotationModel() : IAURotationModel(360.0 / 533.7004100) {}
  428.     
  429.     void pole(double t, double& ra, double& dec) const
  430.     {
  431.         double T = t / 36525.0;
  432.         double J2 = degToRad(24.62 + 45137.2 * T);
  433.         clamp_centuries(T);
  434.         ra = 268.05 - 0.009 * T - 2.11 * sin(J2) + 0.04 * sin(2.0 * J2);
  435.         dec = 64.49 + 0.003 * T - 0.91 * cos(J2) + 0.01 * cos(2.0 * J2);        
  436.     }
  437.     
  438.     double meridian(double t) const
  439.     {
  440.         double T = t / 36525.0;
  441.         double J2 = degToRad(24.62 + 45137.2 * T);
  442.         return 8.56 + 533.7004100 * t + 1.91 * sin(J2) - 0.04 * sin(2.0 * J2);
  443.     }    
  444. };
  447. class IAUIoRotationModel : public IAURotationModel
  448. {
  449. public:
  450.     IAUIoRotationModel() : IAURotationModel(360.0 / 203.4889538) {}
  451.     
  452.     void pole(double t, double& ra, double& dec) const
  453.     {
  454.         double T = t / 36525.0;
  455.         double J3 = degToRad(283.90 + 4850.7 * T);
  456.         double J4 = degToRad(355.80 + 1191.3 * T);
  457.         clamp_centuries(T);
  458.         ra = 268.05 - 0.009 * T + 0.094 * sin(J3) + 0.024 * sin(J4);
  459.         dec = 64.49 + 0.003 * T + 0.040 * cos(J3) + 0.011 * cos(J4);        
  460.     }
  461.     
  462.     double meridian(double t) const
  463.     {
  464.         double T = t / 36525.0;
  465.         double J3 = degToRad(283.90 + 4850.7 * T);
  466.         double J4 = degToRad(355.80 + 1191.3 * T);
  467.         return 200.39 + 203.4889538 * t - 0.085 * sin(J3) - 0.022 * sin(J4);
  468.     }    
  469. };
  472. class IAUEuropaRotationModel : public IAURotationModel
  473. {
  474. public:
  475.     IAUEuropaRotationModel() : IAURotationModel(360.0 / 101.3747235) {}
  476.     
  477.     void pole(double t, double& ra, double& dec) const
  478.     {
  479.         double T = t / 36525.0;
  480.         double J4 = degToRad(355.80 + 1191.3 * T);
  481.         double J5 = degToRad(119.90 + 262.1 * T);
  482.         double J6 = degToRad(229.80 + 64.3 * T);
  483.         double J7 = degToRad(352.35 + 2382.6 * T);
  484.         clamp_centuries(T);
  485.         ra = 268.05 - 0.009 * T + 1.086 * sin(J4) + 0.060 * sin(J5) + 0.015 * sin(J6) + 0.009 * sin(J7);
  486.         dec = 64.49 + 0.003 * T + 0.486 * cos(J4) + 0.026 * cos(J5) + 0.007 * cos(J6) + 0.002 * cos(J7);        
  487.     }
  488.     
  489.     double meridian(double t) const
  490.     {
  491.         double T = t / 36525.0;
  492.         double J4 = degToRad(355.80 + 1191.3 * T);
  493.         double J5 = degToRad(119.90 + 262.1 * T);
  494.         double J6 = degToRad(229.80 + 64.3 * T);
  495.         double J7 = degToRad(352.35 + 2382.6 * T);
  496.         return 36.022 + 101.3747235 * t - 0.980 * sin(J4) - 0.054 * sin(J5) - 0.014 * sin(J6) - 0.008 * sin(J7);
  497.     }    
  498. };
  501. class IAUGanymedeRotationModel : public IAURotationModel
  502. {
  503. public:
  504.     IAUGanymedeRotationModel() : IAURotationModel(360.0 / 50.3176081) {}
  505.     
  506.     void pole(double t, double& ra, double& dec) const
  507.     {
  508.         double T = t / 36525.0;
  509.         double J4 = degToRad(355.80 + 1191.3 * T);
  510.         double J5 = degToRad(119.90 + 262.1 * T);
  511.         double J6 = degToRad(229.80 + 64.3 * T);
  512.         clamp_centuries(T);
  513.         ra = 268.05 - 0.009 * T - 0.037 * sin(J4) + 0.431 * sin(J5) + 0.091 * sin(J6);
  514.         dec = 64.49 + 0.003 * T - 0.016 * cos(J4) + 0.186 * cos(J5) + 0.039 * cos(J6);        
  515.     }
  516.     
  517.     double meridian(double t) const
  518.     {
  519.         double T = t / 36525.0;
  520.         double J4 = degToRad(355.80 + 1191.3 * T);
  521.         double J5 = degToRad(119.90 + 262.1 * T);
  522.         double J6 = degToRad(229.80 + 64.3 * T);
  523.         return 44.064 + 50.3176081 * t + 0.033 * sin(J4) - 0.389 * sin(J5) - 0.082 * sin(J6);
  524.     }    
  525. };
  528. class IAUCallistoRotationModel : public IAURotationModel
  529. {
  530. public:
  531.     IAUCallistoRotationModel() : IAURotationModel(360.0 / 21.5710715) {}
  532.     
  533.     void pole(double t, double& ra, double& dec) const
  534.     {
  535.         double T = t / 36525.0;
  536.         double J5 = degToRad(119.90 + 262.1 * T);
  537.         double J6 = degToRad(229.80 + 64.3 * T);
  538.         double J8 = degToRad(113.35 + 6070.0 * T);
  539.         clamp_centuries(T);
  540.         ra = 268.05 - 0.009 * T - 0.068 * sin(J5) + 0.590 * sin(J6) + 0.010 * sin(J8);
  541.         dec = 64.49 + 0.003 * T - 0.029 * cos(J5) + 0.254 * cos(J6) - 0.004 * cos(J8);        
  542.     }
  543.     
  544.     double meridian(double t) const
  545.     {
  546.         double T = t / 36525.0;
  547.         double J5 = degToRad(119.90 + 262.1 * T);
  548.         double J6 = degToRad(229.80 + 64.3 * T);
  549.         double J8 = degToRad(113.35 + 6070.0 * T);
  550.         return 259.51 + 21.5710715 * t + 0.061 * sin(J5) - 0.533 * sin(J6) - 0.009 * sin(J8);
  551.     }    
  552. };
  555. /*
  556. S1 = 353.32 + 75706.7 * T
  557. S2 = 28.72 + 75706.7 * T
  558. S3 = 177.40 - 36505.5 * T
  559. S4 = 300.00 - 7225.9 * T
  560. S5 = 53.59 - 8968.6 * T
  561. S6 = 143.38 - 10553.5 * T
  562. S7 = 345.20 - 1016.3 * T
  563. S8 = 29.80 - 52.1 * T
  564. S9 = 316.45 + 506.2 * T
  565. */
  567. // Rotations of Saturnian satellites from Seidelmann, _Explanatory Supplement to the
  568. // Astronomical Almanac_ (1992).
  570. class IAUMimasRotationModel : public IAURotationModel
  571. {
  572. public:
  573.     IAUMimasRotationModel() : IAURotationModel(360.0 / 381.9945550) {}
  574.     
  575.     void pole(double t, double& ra, double& dec) const
  576.     {
  577.         double T = t / 36525.0;
  578.         double S3 = degToRad(177.40 - 36505.5 * T);
  579.         clamp_centuries(T);
  580.         ra = 40.66 - 0.036 * T + 13.56 * sin(S3);
  581.         dec = 83.52 - 0.004 * T - 1.53 * cos(S3);
  582.     }
  583.     
  584.     double meridian(double t) const
  585.     {
  586.         double T = t / 36525.0;
  587.         double S3 = degToRad(177.40 - 36505.5 * T);
  588.         double S9 = degToRad(316.45 + 506.2 * T);
  589.         return 337.46 + 381.9945550 * t - 13.48 * sin(S3) - 44.85 * sin(S9);        
  590.     }
  591. };
  594. class IAUEnceladusRotationModel : public IAURotationModel
  595. {
  596. public:
  597.     IAUEnceladusRotationModel() : IAURotationModel(360.0 / 262.7318996) {}
  598.     
  599.     void pole(double t, double& ra, double& dec) const
  600.     {
  601.         double T = t / 36525.0;
  602.         clamp_centuries(T);
  603.         ra = 40.66 - 0.036 * T;
  604.         dec = 83.52 - 0.004 * T;
  605.     }
  606.     
  607.     double meridian(double t) const
  608.     {
  609.         return 2.82 + 262.7318996 * t;
  610.     }
  611. };
  614. class IAUTethysRotationModel : public IAURotationModel
  615. {
  616. public:
  617.     IAUTethysRotationModel() : IAURotationModel(360.0 / 190.6979085) {}
  618.     
  619.     void pole(double t, double& ra, double& dec) const
  620.     {
  621.         double T = t / 36525.0;
  622.         double S4 = degToRad(300.00 - 7225.9 * T);
  623.         clamp_centuries(T);
  624.         ra = 40.66 - 0.036 * T - 9.66 * sin(S4);
  625.         dec = 83.52 - 0.004 * T - 1.09 * cos(S4);
  626.     }
  627.     
  628.     double meridian(double t) const
  629.     {
  630.         double T = t / 36525.0;
  631.         double S4 = degToRad(300.00 - 7225.9 * T);
  632.         double S9 = degToRad(316.45 + 506.2 * T);
  633.         return 10.45 + 190.6979085 * t - 9.60 * sin(S4) + 2.23 * sin(S9);        
  634.     }
  635. };
  638. class IAUTelestoRotationModel : public IAURotationModel
  639. {
  640. public:
  641.     IAUTelestoRotationModel() : IAURotationModel(360.0 / 190.6979330) {}
  642.     
  643.     void pole(double t, double& ra, double& dec) const
  644.     {
  645.         double T = t / 36525.0;
  646.         clamp_centuries(T);
  647.         ra = 50.50 - 0.036 * T;
  648.         dec = 84.06 - 0.004 * T;
  649.     }
  650.     
  651.     double meridian(double t) const
  652.     {
  653.         return 56.88 + 190.6979330 * t;
  654.     }
  655. };
  658. class IAUCalypsoRotationModel : public IAURotationModel
  659. {
  660. public:
  661.     IAUCalypsoRotationModel() : IAURotationModel(360.0 / 190.6742373) {}
  662.     
  663.     void pole(double t, double& ra, double& dec) const
  664.     {
  665.         double T = t / 36525.0;
  666.         double S5 = degToRad(53.59 - 8968.6 * T);
  667.         clamp_centuries(T);
  668.         ra = 40.58 - 0.036 * T - 13.943 * sin(S5) - 1.686 * sin(2.0 * S5);
  669.         dec = 83.43 - 0.004 * T - 1.572 * cos(S5) + 0.095 * cos(2.0 * S5);
  670.     }
  671.     
  672.     double meridian(double t) const
  673.     {
  674.         double T = t / 36525.0;
  675.         double S5 = degToRad(53.59 - 8968.6 * T);
  676.         return 149.36 + 190.6742373 * t - 13.849 * sin(S5) + 1.685 * sin(2.0 * S5);
  677.     }
  678. };
  681. class IAUDioneRotationModel : public IAURotationModel
  682. {
  683. public:
  684.     IAUDioneRotationModel() : IAURotationModel(360.0 / 131.5349316) {}
  685.     
  686.     void pole(double t, double& ra, double& dec) const
  687.     {
  688.         double T = t / 36525.0;
  689.         clamp_centuries(T);
  690.         ra = 40.66 - 0.036 * T;
  691.         dec = 83.52 - 0.004 * T;
  692.     }
  693.     
  694.     double meridian(double t) const
  695.     {
  696.         return 357.00 + 131.5349316 * t;
  697.     }
  698. };
  701. class IAUHeleneRotationModel : public IAURotationModel
  702. {
  703. public:
  704.     IAUHeleneRotationModel() : IAURotationModel(360.0 / 131.6174056) {}
  705.     
  706.     void pole(double t, double& ra, double& dec) const
  707.     {
  708.         double T = t / 36525.0;
  709.         double S6 = degToRad(143.38 - 10553.5 * T);
  710.         clamp_centuries(T);
  711.         ra = 40.58 - 0.036 * T + 1.662 * sin(S6) + 0.024 * sin(2.0 * S6);
  712.         dec = 83.52 - 0.004 * T - 0.187 * cos(S6) + 0.095 * cos(2.0 * S6);
  713.     }
  714.     
  715.     double meridian(double t) const
  716.     {
  717.         double T = t / 36525.0;
  718.         double S6 = degToRad(143.38 - 10553.5 * T);
  719.         return 245.39 + 131.6174056 * t - 1.651 * sin(S6) + 0.024 * sin(2.0 * S6);
  720.     }
  721. };
  724. class IAURheaRotationModel : public IAURotationModel
  725. {
  726. public:
  727.     IAURheaRotationModel() : IAURotationModel(360.0 / 79.6900478) {}
  728.     
  729.     void pole(double t, double& ra, double& dec) const
  730.     {
  731.         double T = t / 36525.0;
  732.         double S7 = degToRad(345.20 - 1016.3 * T);
  733.         clamp_centuries(T);
  734.         ra = 40.38 - 0.036 * T + 3.10 * sin(S7);
  735.         dec = 83.55 - 0.004 * T - 0.35 * cos(S7);
  736.     }
  737.     
  738.     double meridian(double t) const
  739.     {
  740.         double T = t / 36525.0;
  741.         double S7 = degToRad(345.20 - 1016.3 * T);
  742.         return 235.16 + 79.6900478 * t - 1.651 - 3.08 * sin(S7);
  743.     }
  744. };
  747. class IAUTitanRotationModel : public IAURotationModel
  748. {
  749. public:
  750.     IAUTitanRotationModel() : IAURotationModel(360.0 / 22.5769768) {}
  751.     
  752.     void pole(double t, double& ra, double& dec) const
  753.     {
  754.         double T = t / 36525.0;
  755.         double S8 = degToRad(29.80 - 52.1 * T);
  756.         clamp_centuries(T);
  757.         ra = 36.41 - 0.036 * T + 2.66 * sin(S8);
  758.         dec = 83.94 - 0.004 * T - 0.30 * cos(S8);
  759.     }
  760.     
  761.     double meridian(double t) const
  762.     {
  763.         double T = t / 36525.0;
  764.         double S8 = degToRad(29.80 - 52.1 * T);
  765.         return 189.64 + 22.5769768 * t - 2.64 * sin(S8);
  766.     }
  767. };
  770. class IAUIapetusRotationModel : public IAURotationModel
  771. {
  772. public:
  773.     IAUIapetusRotationModel() : IAURotationModel(360.0 / 4.5379572) {}
  774.     
  775.     void pole(double t, double& ra, double& dec) const
  776.     {
  777.         double T = t / 36525.0;
  778.         clamp_centuries(T);
  779.         ra = 318.16 - 3.949 * T;
  780.         dec = 75.03 - 1.142 * T;
  781.     }
  782.     
  783.     double meridian(double t) const
  784.     {
  785.         return 350.20 + 4.5379572 * t;
  786.     }
  787. };
  790. class IAUPhoebeRotationModel : public IAURotationModel
  791. {
  792. public:
  793.     IAUPhoebeRotationModel() : IAURotationModel(360.0 / 22.5769768) {}
  794.     
  795.     void pole(double t, double& ra, double& dec) const
  796.     {
  797.         double T = t / 36525.0;
  798.         clamp_centuries(T);
  799.         ra = 355.16;
  800.         dec = 68.70 - 1.143 * T;
  801.     }
  802.     
  803.     double meridian(double t) const
  804.     {
  805.         return 304.70 + 930.8338720 * t;
  806.     }
  807. };
  810. class IAUMirandaRotationModel : public IAURotationModel
  811. {
  812. public:
  813.     IAUMirandaRotationModel() : IAURotationModel(360.0 / 254.6906892) { setFlipped(true); }
  814.     
  815.     void pole(double t, double& ra, double& dec) const
  816.     {
  817.         double T = t / 36525.0;
  818.         double U11 = degToRad(102.23 - 2024.22 * T);
  819.         ra =  257.43 + 4.41 * sin(U11) - 0.04 * sin(2.0 * U11);
  820.         dec = -15.08 + 4.25 * cos(U11) - 0.02 * cos(2.0 * U11);
  821.     }
  822.     
  823.     double meridian(double t) const
  824.     {
  825.         double T = t / 36525.0;
  826.         double U11 = degToRad(102.23 - 2024.22 * T);
  827.         double U12 = degToRad(316.41 + 2863.96 * T);
  828.         return 30.70 - 254.6906892 * t
  829.                - 1.27 * sin(U12) + 0.15 * sin(2.0 * U12)
  830.                + 1.15 * sin(U11) - 0.09 * sin(2.0 * U11);
  831.     }
  832. };
  835. class IAUArielRotationModel : public IAURotationModel
  836. {
  837. public:
  838.     IAUArielRotationModel() : IAURotationModel(360.0 / 142.8356681) { setFlipped(true); }
  839.     
  840.     void pole(double t, double& ra, double& dec) const
  841.     {
  842.         double T = t / 36525.0;
  843.         double U13 = degToRad(304.01 - 51.94 * T);
  844.         ra =  257.43 + 0.29 * sin(U13);
  845.         dec = -15.10 + 0.28 * cos(U13);
  846.     }
  847.     
  848.     double meridian(double t) const
  849.     {
  850.         double T = t / 36525.0;
  851.         double U12 = degToRad(316.41 + 2863.96 * T);
  852.         double U13 = degToRad(304.01 - 51.94 * T);
  853.         return 156.22 - 142.8356681 * t + 0.05 * sin(U12) + 0.08 * sin(U13);
  854.     }
  855. };
  858. class IAUUmbrielRotationModel : public IAURotationModel
  859. {
  860. public:
  861.     IAUUmbrielRotationModel() : IAURotationModel(360.0 / 86.8688923) { setFlipped(true); }
  862.     
  863.     void pole(double t, double& ra, double& dec) const
  864.     {
  865.         double T = t / 36525.0;
  866.         double U14 = degToRad(308.71 - 93.17 * T);
  867.         ra =  257.43 + 0.21 * sin(U14);
  868.         dec = -15.10 + 0.20 * cos(U14);
  869.     }
  870.     
  871.     double meridian(double t) const
  872.     {
  873.         double T = t / 36525.0;
  874.         double U12 = degToRad(316.41 + 2863.96 * T);
  875.         double U14 = degToRad(308.71 - 93.17 * T);
  876.         return 108.05 - 86.8688923 * t - 0.09 * sin(U12) + 0.06 * sin(U14);
  877.     }
  878. };
  881. class IAUTitaniaRotationModel : public IAURotationModel
  882. {
  883. public:
  884.     IAUTitaniaRotationModel() : IAURotationModel(360.0 / 41.351431) { setFlipped(true); }
  885.     
  886.     void pole(double t, double& ra, double& dec) const
  887.     {
  888.         double T = t / 36525.0;
  889.         double U15 = degToRad(340.82 - 75.32 * T);
  890.         ra =  257.43 + 0.29 * sin(U15);
  891.         dec = -15.10 + 0.28 * cos(U15);
  892.     }
  893.     
  894.     double meridian(double t) const
  895.     {
  896.         double T = t / 36525.0;
  897.         double U15 = degToRad(340.82 - 75.32 * T);
  898.         return 77.74 - 41.351431 * t + 0.08 * sin(U15);
  899.     }
  900. };
  903. class IAUOberonRotationModel : public IAURotationModel
  904. {
  905. public:
  906.     IAUOberonRotationModel() : IAURotationModel(360.0 / 26.7394932) { setFlipped(true); }
  907.     
  908.     void pole(double t, double& ra, double& dec) const
  909.     {
  910.         double T = t / 36525.0;
  911.         double U16 = degToRad(259.14 - 504.81 * T);
  912.         ra =  257.43 + 0.16 * sin(U16);
  913.         dec = -15.10 + 0.16 * cos(U16);
  914.     }
  915.     
  916.     double meridian(double t) const
  917.     {
  918.         double T = t / 36525.0;
  919.         double U16 = degToRad(259.14 - 504.81 * T);
  920.         return 6.77 - 26.7394932 * t + 0.04 * sin(U16);
  921.     }
  922. };
  926. RotationModel*
  927. GetCustomRotationModel(const std::string& name)
  928. {
  929.     // Initialize the custom rotation model table.
  930.     if (!CustomRotationModelsInitialized)
  931.     {
  932.         CustomRotationModelsInitialized = true;
  933.         
  934.         CustomRotationModels["earth-p03lp"] = new EarthRotationModel();
  935.         
  936.         // IAU rotation elements for the planets
  937.         CustomRotationModels["iau-mercury"] = new IAUPrecessingRotationModel(281.01, -0.033,
  938.                                                                              61.45, -0.005,
  939.                                                                              329.548, 6.1385025);
  940.         CustomRotationModels["iau-venus"]   = new IAUPrecessingRotationModel(272.76, 0.0,
  941.                                                                              67.16, 0.0,
  942.                                                                              160.20, -1.4813688);
  943.         CustomRotationModels["iau-earth"]   = new IAUPrecessingRotationModel(0.0, -0.641,
  944.                                                                              90.0, -0.557,
  945.                                                                              190.147, 360.9856235);
  946.         CustomRotationModels["iau-mars"]    = new IAUPrecessingRotationModel(317.68143, -0.1061,
  947.                                                                              52.88650, -0.0609,
  948.                                                                              176.630, 350.89198226);
  949.         CustomRotationModels["iau-jupiter"] = new IAUPrecessingRotationModel(268.05, -0.009,
  950.                                                                              64.49, -0.003,
  951.                                                                              284.95, 870.5366420);
  952.         CustomRotationModels["iau-saturn"]  = new IAUPrecessingRotationModel(40.589, -0.036,
  953.                                                                              83.537, -0.004,
  954.                                                                              38.90, 810.7939024);
  955.         CustomRotationModels["iau-uranus"]  = new IAUPrecessingRotationModel(257.311, 0.0,
  956.                                                                              -15.175, 0.0,
  957.                                                                              203.81, -501.1600928);
  958.         CustomRotationModels["iau-neptune"] = new IAUNeptuneRotationModel();
  959.         CustomRotationModels["iau-pluto"]   = new IAUPrecessingRotationModel(313.02, 0.0,
  960.                                                                              9.09, 0.0,
  961.                                                                              236.77, -56.3623195);
  962.        
  963.         // IAU elements for satellite of Earth
  964.         CustomRotationModels["iau-moon"] = new IAULunarRotationModel();
  965.         
  966.         // IAU elements for satellites of Mars
  967.         CustomRotationModels["iau-phobos"] = new IAUPhobosRotationModel();
  968.         CustomRotationModels["iau-deimos"] = new IAUDeimosRotationModel();
  969.         
  970.         // IAU elements for satellites of Jupiter
  971.         CustomRotationModels["iau-metis"]    = new IAUPrecessingRotationModel(268.05, -0.009,
  972.                                                                               64.49, 0.003,
  973.                                                                               346.09, 1221.2547301);
  974.         CustomRotationModels["iau-adrastea"] = new IAUPrecessingRotationModel(268.05, -0.009,
  975.                                                                               64.49, 0.003,
  976.                                                                               33.29, 1206.9986602);
  977.         CustomRotationModels["iau-amalthea"] = new IAUAmaltheaRotationModel();
  978.         CustomRotationModels["iau-thebe"]    = new IAUThebeRotationModel();
  979.         CustomRotationModels["iau-io"]       = new IAUIoRotationModel();
  980.         CustomRotationModels["iau-europa"]   = new IAUEuropaRotationModel();
  981.         CustomRotationModels["iau-ganymede"] = new IAUGanymedeRotationModel();
  982.         CustomRotationModels["iau-callisto"] = new IAUCallistoRotationModel();
  983.         
  984.         // IAU elements for satellites of Saturn
  985.         CustomRotationModels["iau-pan"]        = new IAUPrecessingRotationModel(40.6, -0.036,
  986.                                                                                 83.5, -0.004,
  987.                                                                                 48.8, 626.0440000);
  988.         CustomRotationModels["iau-atlas"]      = new IAUPrecessingRotationModel(40.6, -0.036,
  989.                                                                                 83.5, -0.004,
  990.                                                                                 137.88, 598.3060000);
  991.         CustomRotationModels["iau-prometheus"] = new IAUPrecessingRotationModel(40.6, -0.036,
  992.                                                                                 83.5, -0.004,
  993.                                                                                 296.14, 587.289000);
  994.         CustomRotationModels["iau-pandora"]    = new IAUPrecessingRotationModel(40.6, -0.036,
  995.                                                                                 83.5, -0.004,
  996.                                                                                 162.92, 572.7891000);
  997.         CustomRotationModels["iau-mimas"]     = new IAUMimasRotationModel();
  998.         CustomRotationModels["iau-enceladus"] = new IAUEnceladusRotationModel();
  999.         CustomRotationModels["iau-tethys"]    = new IAUTethysRotationModel();
  1000.         CustomRotationModels["iau-telesto"]   = new IAUTelestoRotationModel();
  1001.         CustomRotationModels["iau-calypso"]   = new IAUCalypsoRotationModel();
  1002.         CustomRotationModels["iau-dione"]     = new IAUDioneRotationModel();
  1003.         CustomRotationModels["iau-helene"]    = new IAUHeleneRotationModel();
  1004.         CustomRotationModels["iau-rhea"]      = new IAURheaRotationModel();
  1005.         CustomRotationModels["iau-titan"]     = new IAUTitanRotationModel();
  1006.         CustomRotationModels["iau-iapetus"]   = new IAUIapetusRotationModel();
  1007.         CustomRotationModels["iau-phoebe"]    = new IAUPhoebeRotationModel();
  1008.         
  1009.         CustomRotationModels["iau-miranda"]   = new IAUMirandaRotationModel();
  1010.         CustomRotationModels["iau-ariel"]     = new IAUArielRotationModel();
  1011.         CustomRotationModels["iau-umbriel"]   = new IAUUmbrielRotationModel();
  1012.         CustomRotationModels["iau-titania"]   = new IAUTitaniaRotationModel();
  1013.         CustomRotationModels["iau-oberon"]    = new IAUOberonRotationModel();
  1014.     }
  1016.     if (CustomRotationModels.count(name) > 0)
  1017.         return CustomRotationModels[name];
  1018.     else
  1019.         return NULL;
  1020. }