from math import pi, sin, cos, fmod, tan, atan, fabs, atan2, asin, acos, sqrt, log10
import doctest
from celme.dates import Date
from celme.angles import Angle
from celme.coords import Coords

# ========================================================
# ========================================================
# === Mechanics
# ========================================================
# ========================================================

class Mechanics:
    """ Class to compute celestial mechanics
    """
# ========================================================
# === attributs
# ========================================================

    _DR = pi/180
    _PI = pi
    _PISUR2 = pi/2
    
    # analytical fit methods
    _AFM_VFP79 = 0
    _AFM_VSOP87 = 1
    _AFM_GMS86 = 2
    _AFM_ELP2000_82B = 3
    _AFM_MELP98 = 4
    
    # planet names
    _PLN_OTHERPLANET = -2
    _PLN_ALLPLANETS = -1
    _PLN_SOLEIL = 0
    _PLN_MERCURE = 1
    _PLN_VENUS = 2
    _PLN_TERRE = 3
    _PLN_MARS = 4
    _PLN_JUPITER = 5
    _PLN_SATURNE = 6
    _PLN_URANUS = 7
    _PLN_NEPTUNE = 8
    _PLN_PLUTON = 9
    _PLN_LUNE_ELP = 10
    _PLN_LUNE = 11
    _PLN_NB_PLANETES = 12
    
    # astronomical constants
    _CST_UA = 1.49597870691e11 ; # (m)
    _CST_CLIGHT = 2.99792458e8 ; # (m/s)
    _CST_EARTH_SEMI_MAJOR_RADIUS = 6378137 ; # (m) WGS 84
    _CST_EARTH_INVERSE_FLATTENING = 298.257223563 ; # WGS 84
    
# ========================================================
# === internal methods : Generals
# ========================================================

    def _init_mechanics(self):
        """ Object initialization where planet_name is the name of the planet
        
        :param planet_name: A string (cf. help(PLanet))
        :type planet_name: string
        """
        pass
    
# ========================================================
# === internal methods : Planet orbits
# ========================================================
 
    def _mc_name2planetnum(self, planet_name):
        name = str(planet_name[0:3]).upper()
        planet_num = self._PLN_OTHERPLANET
        if (name == "SUN"):
            planet_num = self._PLN_SOLEIL
        elif (name == "MER"):
            planet_num = self._PLN_MERCURE
        elif (name == "VEN"):
            planet_num = self._PLN_VENUS
        elif (name == "EAR"):
            planet_num = self._PLN_TERRE
        elif (name == "MAR"):
            planet_num = self._PLN_MARS
        elif (name == "JUP"):
            planet_num = self._PLN_JUPITER
        elif (name == "SAT"):
            planet_num = self._PLN_SATURNE
        elif (name == "URA"):
            planet_num = self._PLN_URANUS
        elif (name == "NEP"):
            planet_num = self._PLN_NEPTUNE
        elif (name == "PLU"):
            planet_num = self._PLN_PLUTON
        elif (name == "ELP"):
            planet_num = self._PLN_LUNE_ELP
        elif (name == "MOO"):
            planet_num = self._PLN_LUNE
        return planet_num
        
    def _mc_jd2lbr1a(self, jj):
        """
        /***************************************************************************/
        /* Retourne les valeurs des tableaux L, M, U necessaires pour tenir        */
        /* compte des principales perturbations planetaires a entrer dans          */
        /* la fonction mc_jd2lbr1b                                                 */
        /***************************************************************************/
        /* Les tableaux *l, *m et *u doivent etre dimensiones chacun avec          */
        /* 9 elements (0 a 8)                                                      */
        /***************************************************************************/
        """
        #  l, m, u
        t=(jj-2415020.0)/36525;
        l=[0]*self._PLN_NB_PLANETES
        m=[0]*self._PLN_NB_PLANETES
        u=[0]*self._PLN_NB_PLANETES
        l[0]=(279.6964027+36000.7695173*t)*self._DR;
        m[0]=(358.4758635+35999.0494965*t)*self._DR;
        u[0]=(270.435377+481267.880863*t)*self._DR;
        l[1]=(178.178814+149474.071386*t)*self._DR;
        m[1]=(102.279426+149472.515334*t)*self._DR;
        u[1]=(131.032888+149472.885872*t)*self._DR;
        l[2]=(342.766738+58519.212542*t)*self._DR;
        m[2]=(212.601892+58517.806388*t)*self._DR;
        u[2]=(266.987445+58518.311835*t)*self._DR;
        l[3]=(293.747201+19141.699879*t)*self._DR;
        m[3]=(319.529273+19139.858887*t)*self._DR;
        u[3]=(244.960887+19140.928953*t)*self._DR;
        l[4]=(237.352259+3034.906621*t)*self._DR;
        m[4]=(225.444539+3034.906621*t)*self._DR;
        u[4]=(138.419219+3034.906621*t)*self._DR;
        l[5]=(265.869357+1222.116843*t)*self._DR;
        m[5]=(175.758477+1222.116843*t)*self._DR;
        u[5]=(153.521637+1222.116843*t)*self._DR;
        l[6]=(243.362437+429.898403*t)*self._DR;
        m[6]=( 74.313637+429.898403*t)*self._DR;
        u[6]=(169.872293+429.388747*t)*self._DR;
        l[7]=( 85.024943+219.863377*t)*self._DR;
        m[7]=( 41.269103+219.863377*t)*self._DR;
        u[7]=(314.346275+218.761885*t)*self._DR;
        l[8]=( 92.312712+146.674728*t)*self._DR;
        m[8]=(229.488633+145.278567*t)*self._DR;
        u[8]=(343.369233+145.278567*t)*self._DR;
        res = (l, m, u)
        return res
    
    def _mc_jd2lbr1b(self, jj, planete, l, m, u, afm):
        """
        /***************************************************************************/
        /* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr  */
        /* pour l'equinoxe moyen de la date. ()                                    */
        /***************************************************************************/
        /* Les tableaux l, m et u doivent etre dimensiones chacun avec             */
        /* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a     */
        /***************************************************************************/
        """
        coords = Coords((0,0,0))
        if afm == self._AFM_VFP79:
            t=(jj-2415020.0)/36525
            if ((planete==self._PLN_SOLEIL) or (planete==self._PLN_TERRE)):
                l0=(1.91944-.004722*t)*sin(m[0])+.02*sin(2*m[0])-.001944*cos(m[0]-m[4])
                +.001666*sin(u[0]-l[0])+.001388*sin(4*m[0]-8*m[3]+m[4])-.001388*cos(m[0]-m[2])
                l0+=(-.001111*sin(m[0]-m[2])+.001111*cos(4*m[0]-8*m[3]+m[4])+.000833
                *sin(2*m[0]-m[2])-.000833*sin(m[4])-.000833*sin(2*m[0]-m[4]))
                l0=l[0]+l0*self._DR
                b0=0
                e=.01675104-.418e-4*t-.126e-6*t*t;
                r0=1.0000002*(1+.5*e*e-(e-3./8*e*e*e)*cos(m[0])-.5*e*e*cos(2*m[0])-3.*e*e*e/8*cos(3*m[0]))
                if (planete==self._PLN_TERRE):
                    l0+=self._PI
                    b0=-b0
              
            elif (planete==self._PLN_MERCURE):
                l0=8*t*sin(m[1])+84378*sin(m[1])+10733*sin(2*m[1])+1892*sin(3*m[1])
                +381*sin(4*m[1])+83*sin(5*m[1])+19*sin(6*m[1])-646*sin(2*u[1])-306*sin(m[1]-2*u[1])
                -274*sin(m[1]-2*u[1])-92*sin(2*m[1]+2*u[1])-28*sin(3*m[1]+2*u[1])+25*sin(2*m[1]-2*u[1])
                l0+=(-9*sin(4*m[1]+2*u[1])+7*cos(2*m[1]-5*m[2]))
                l0=l[1]+l0/3600*self._DR
                b0=24134*sin(u[1])-10*sin(3*u[1])+5180*sin(m[1]-u[1])+4910*sin(m[1]+u[1])
                +1124*sin(2*m[1]+u[1])+271*sin(3*m[1]+u[1])+132*sin(2*m[1]-u[1])
                +67*sin(4*m[1]+u[1])+18*sin(3*m[1]-u[1])+17*sin(5*m[1]+u[1])-9*sin(m[1]-3*u[1])
                b0=b0/3600*self._DR
                r0=0.39528-.07834*cos(m[1])-.00795*cos(2*m[1])-.00121*cos(3*m[1])
                -.00022*cos(4*m[1])
            
            elif (planete==self._PLN_VENUS):
                l0=20*t*sin(m[2])+2814*sin(m[2])+12*sin(2*m[2])-181*sin(2*u[2])
                -10*cos(2*m[0]-2*m[2])+7*cos(3*m[0]-3*m[2]);
                l0=l[2]+l0/3600*self._DR;
                b0=12215*sin(u[2])+83*sin(m[2]+u[2])+83*sin(m[2]-u[2]);
                b0=b0/3600*self._DR;
                r0=.72335-0.00493*cos(m[2]);
            
            elif (planete==self._PLN_MARS):
                l0=37*t*sin(m[3])+4*t*sin(2*m[3])+38451*sin(m[3])+2238*sin(2*m[3])
                +181*sin(3*m[3])+17*sin(4*m[3])-52*sin(2*u[3])-22*cos(m[3]-2*m[4])-19*sin(m[3]-m[4])
                +17*cos(m[3]-m[4])-16*cos(2*m[3]-2*m[4])+13*cos(m[0]-2*m[3])-10*sin(m[3]-2*u[3]);
                l0+=(-10*sin(m[3]+2*u[3])+7*cos(m[0]-m[3])-7*cos(2*m[0]-3*m[3])-5*sin(m[2]-3*m[3])
                -5*sin(m[0]-m[3])-5*sin(m[0]-2*m[3])-4*cos(2*m[0]-4*m[3])+4*cos(m[4])+3*cos(m[2]-3*m[3])
                +3*sin(2*m[3]-2*m[4]));
                l0=l[3]+l0/3600*self._DR
                b0=6603*sin(u[3])+622*sin(m[3]-u[3])+615*sin(m[3]+u[3])+64*sin(2*m[3]+u[3])
                b0=b0/3600*self._DR
                r0=1.53031-.14170*cos(m[3])-.0066*cos(2*m[3])-.00047*cos(3*m[3])
                
            elif (planete==self._PLN_JUPITER):
                l0,b0,r0 = self._mc_jd2lbr1c(jj,l,m,u)
            
            elif (planete==self._PLN_SATURNE):
                l0,b0,r0 = self._mc_jd2lbr1e(jj,l,m,u)
            
            elif (planete==self._PLN_URANUS):
                l0,b0,r0 = self._mc_jd2lbr1f(jj,l,m,u)
            
            elif (planete==self._PLN_NEPTUNE):
                l0,b0,r0 = self._mc_jd2lbr1g(jj,l,m,u)
            
            elif (planete==self._PLN_PLUTON):
                l0,b0,r0 = self._mc_jd2lbr1h(jj,l,m,u)
                
            elif (planete==self._PLN_LUNE):
                l0,b0,r0 = self._mc_jd2lbr1d(jj)
            
            ll = Angle(l0/self._DR)
            bb = Angle(b0/self._DR)
            coords = Coords((r0, ll, bb))
        return coords

    def _mc_jd2lbr1c(self, jj, l, m, u):
        """
        /***************************************************************************/
        /* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr  */
        /* pour l'equinoxe moyen de la date. (longitudes vraies)                   */
        /* JUPITER                                                                 */
        /***************************************************************************/
        /* Les tableaux l, m et u doivent etre dimensiones chacun avec             */
        /* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a     */
        /***************************************************************************/
        """
        u[0]*=1.;
        t=(jj-2415020.0)/36525;
        l0=2511+5023*t+19934*sin(m[4])+601*sin(2*m[4])+1093*cos(2*m[4]-5*m[5])
        -479*sin(2*m[4]-5*m[5])-185*sin(2*m[4]-2*m[5])+137*sin(3*m[4]-5*m[5])-131*sin(m[4]-2*m[5])
        +79*cos(m[4]-m[5])-76*cos(2*m[4]-2*m[5])-74*t*cos(m[4])+68*t*sin(m[4]);
        l0+=(+66*cos(2*m[4]-3*m[5])+63*cos(3*m[4]-5*m[5])+53*cos(m[4]-5*m[5])+49*sin(2*m[4]-3*m[5])
        -43*t*sin(2*m[4]-5*m[5])-37*cos(m[4])+25*sin(2*l[4])+25*sin(3*m[4])-23*sin(m[4]-5*m[5])
        -19*t*cos(2*m[4]-5*m[5])+17*cos(2*m[4]-4*m[5]));
        l0+=(+17*cos(3*m[4]-3*m[5])-14*sin(m[4]-m[5]));
        l0+=(-13*sin(3*m[4]-4*m[5])-9*cos(2*l[4])+9*cos(m[5])-9*sin(m[5])-9*sin(3*m[4]-2*m[5])
        +9*sin(4*m[4]-5*m[5])+9*sin(2*m[4]-6*m[5]+3*m[6])-8*cos(4*m[4]-10*m[5])+7*cos(3*m[4]-4*m[5])
        -7*cos(m[4]-3*m[5])-7*sin(4*m[4]-10*m[5]));
        l0+=(-7*sin(m[4]-3*m[5])+6*cos(4*m[4]-5*m[5]));
        l0+=(-6*sin(3*m[4]-3*m[5])+5*cos(2*m[5])-4*sin(4*m[4]-4*m[5])-4*cos(3*m[5])+4*cos(2*m[4]-m[5])
        -4*cos(3*m[4]-2*m[5])-4*t*cos(2*m[4])+3*t*sin(2*m[4])+3*cos(5*m[5])+3*cos(5*m[4]-10*m[5])
        +3*sin(2*m[5])-2*sin(2*l[4]-m[4])+2*sin(2*l[4]+m[4]));
        l0+=(-2*t*sin(3*m[4]-5*m[5])-2*t*sin(m[4]-5*m[5]));
        l0=l[4]+l0/3600*self._DR;
        b0=-4692*cos(m[4])+259*sin(m[4])+227-227*cos(2*m[4])+30*t*sin(m[4])+21*t*cos(m[4])
        +16*sin(3*m[4]-5*m[5])-13*sin(m[4]-5*m[5])-12*cos(3*m[4])+12*sin(2*m[4])+7*cos(3*m[4]-5*m[5])
        -5*cos(m[4]-5*m[5]);
        b0=b0/3600*self._DR;
        r0=5.20883-.25122*cos(m[4])-.00604*cos(2*m[4])+.0026*cos(2*m[4]-2*m[5])
        -.00170*cos(3*m[4]-5*m[5])-.0016*sin(2*m[4]-2*m[5])-.00091*t*sin(m[4])
        -.00084*t*cos(m[4])+.00069*sin(2*m[4]-3*m[5])-.00067*sin(m[4]-5*m[5]);
        r0+=(.00066*sin(3*m[4]-5*m[5])+.00063*sin(m[4]-m[5])-.00051*cos(2*m[4]-3*m[5])
        -.00046*sin(m[4])-.00029*cos(m[4]-5*m[5])+.00027*cos(m[4]-2*m[5])
        -.00022*cos(3*m[4])-.00021*sin(2*m[4]-5*m[5]));
        ll0=l0
        bb0=b0
        rr0=r0
        return ll0, bb0, rr0
        
    def _mc_jd2lbr1e(self, jj, l, m, u):
        """
        /***************************************************************************/
        /* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr  */
        /* pour l'equinoxe moyen de la date. (longitudes vraies)                   */
        /* SATURNE                                                                 */
        /***************************************************************************/
        /* Les tableaux l, m et u doivent etre dimensiones chacun avec             */
        /* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a     */
        /***************************************************************************/
        """
        u[0]*=1.;
        t=(jj-2415020.0)/36525;
        l0=2507+5014*t+23043*sin(m[5])-2689*cos(2*m[4]-5*m[5])+1177*sin(2*m[4]-5*m[5])-826*cos(2*m[4]-4*m[5])+802*sin(2*m[5])+425*sin(m[4]-2*m[5])-229*t*cos(m[5])-142*t*sin(m[5])-153*cos(2*m[4]-6*m[5])-114*cos(m[5])+101*t*sin(2*m[4]-5*m[5]);
        l0+=-70*cos(2*l[5])+67*sin(2*l[5])+66*sin(2*m[4]-6*m[5])+60*t*cos(2*m[4]-5*m[5])+41*sin(m[4]-3*m[5])+39*sin(3*m[5])+31*sin(m[4]-m[5])+31*sin(2*m[4]-2*m[5])-29*cos(2*m[4]-3*m[5])-28*sin(2*m[4]-6*m[5]+3*m[6])+28*cos(m[4]-3*m[5]);
        l0+=+22*t*sin(2*m[4]-4*m[5])-22*sin(m[5]-3*m[6])+20*sin(2*m[4]-3*m[5])+20*cos(4*m[4]-10*m[5])+19*cos(2*m[5]-3*m[6])+19*sin(4*m[4]-10*m[5])-17*t*cos(2*m[5])-16*cos(m[5]-3*m[6])-12*sin(2*m[4]-4*m[5])+12*cos(m[4]);
        l0+=-12*sin(2*m[5]-2*m[6])-11*t*sin(2*m[5])-11*cos(2*m[4])-12*sin(2*m[5]-2*m[6])-11*t*sin(2*m[5])-11*cos(2*m[4]-7*m[5])+10*sin(2*m[5]-3*m[6])+10*cos(2*m[4]-2*m[5])+9*sin(4*m[4]-9*m[5])-8*sin(m[5]-2*m[6])-8*cos(2*l[5]+m[5]);
        l0+=+8*cos(2*l[5]-m[5])+8*cos(m[5]-m[6])-8*sin(2*l[5]-m[5])+7*sin(2*l[5]+m[5])-7*cos(m[4]-2*m[5])-7*cos(2*m[5])-6*t*sin(4*m[4]-10*m[5])+6*t*cos(4*m[4]-10*m[5])+6*t*(2*m[4]-6*m[5])-5*sin(3*m[4]-7*m[5])-5*cos(3*m[4]-3*m[5]);
        l0+=-5*cos(2*m[5]-2*m[6])+5*sin(3*m[4]-4*m[5])+5*sin(2*m[4]-7*m[5])+4*sin(3*m[4]-3*m[5])+4*sin(3*m[4]-5*m[5])+4*t*cos(m[4]-3*m[5])+3*t*cos(2*m[4]-4*m[5])+3*cos(2*m[4]-6*m[5]+3*m[6])-3*t*sin(2*l[5]);
        l0+=+3*t*cos(2*m[4]-6*m[5])-3*t*cos(2*l[5])+3*cos(3*m[4]-7*m[5])+3*cos(4*m[4]-9*m[5])+3*sin(3*m[4]-6*m[5])+3*sin(2*m[4]-m[5])+3*sin(m[4]-4*m[5])+2*cos(3*m[5]-3*m[6])+2*t*sin(m[4]-2*m[5])+2*sin(4*m[5])-2*cos(3*m[4]-4*m[5])-2*cos(2*m[4]-m[5]);
        l0+=-2*sin(2*m[4]-7*m[5]+3*m[6])+2*cos(m[4]-4*m[5])+2*cos(4*m[4]-11*m[5])-2*sin(m[5]-m[6]);
        l0=l[5]+l0/3600*self._DR;
        b0=185+8297*sin(m[5])-3346*cos(m[5])+462*sin(2*m[5])-189*cos(2*m[5])+79*t*cos(m[5])-71*cos(2*m[4]-4*m[5])+46*sin(2*m[4]-6*m[5])-45*cos(2*m[4]-6*m[5])+29*sin(3*m[5])-20*cos(2*m[4]-3*m[5])+18*t*sin(m[5]);
        b0+=-14*cos(2*m[4]-5*m[5])-11*cos(3*m[5])-10*t+9*sin(m[4]-3*m[5])+8*sin(m[4]-m[5])-6*sin(2*m[4]-3*m[5])+5*sin(2*m[4]-7*m[5])-5*cos(2*m[4]-7*m[5])+4*sin(2*m[4]-5*m[5])-4*t*sin(2*m[5])-4*cos(m[4]-m[5])+3*cos(m[4]-3*m[5])+3*t*sin(2*m[4]-4*m[5]);
        b0+=+3*sin(m[4]-2*m[5])+2*sin(4*m[5])-2*cos(2*m[4]-2*m[5]);
        b0=b0/3600*self._DR;
        r0=9.55774-.00028*t-.53252*cos(m[5])-.01878*sin(2*m[4]-4*m[5])-.01482*cos(2*m[5])+.00817*sin(m[4]-m[5])-.00539*cos(m[4]-2*m[5])-.00524*t*sin(m[5])+.00349*sin(2*m[4]-5*m[5])+.00347*sin(2*m[4]-6*m[5]);
        r0+=+.00328*t*cos(m[5])-.00225*sin(m[5])+.00149*cos(2*m[4]-6*m[5])-.00126*cos(2*m[4]-2*m[5])+.00104*cos(m[4]-m[5])+.00101*cos(2*m[4]-5*m[5])+.00098*cos(m[4]-3*m[5])-.00073*cos(2*m[4]-3*m[5])-.00062*cos(3*m[5]);
        r0+=+.00043*sin(2*m[5]-3*m[6])+.00041*sin(2*m[4]-2*m[5])-.00040*sin(m[4]-3*m[5])+.0004*cos(2*m[4]-4*m[5])-.00023*sin(m[4])+.0002*sin(2*m[4]-7*m[5]);
        ll0=l0
        bb0=b0
        rr0=r0
        return ll0, bb0, rr0
        
    def _mc_jd2lbr1f(self, jj, l, m, u):
        """
        /***************************************************************************/
        /* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr  */
        /* pour l'equinoxe moyen de la date. (longitudes vraies)                   */
        /* URANUS                                                                  */
        /***************************************************************************/
        /* Les tableaux l, m et u doivent etre dimensiones chacun avec             */
        /* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a     */
        /***************************************************************************/
        """
        u[0]*=1.;
        t=(jj-2415020.0)/36525;
        l0=32*t*t+19397*sin(m[6])+570*sin(2*m[6])-536*t*cos(m[6])+143*sin(m[5]-2*m[6])+110*t*sin(m[6])+102*sin(m[5]-3*m[6])+76*cos(m[5]-3*m[6])-49*sin(m[4]-m[6])-30*t*cos(2*m[6])+29*sin(2*m[4]-6*m[5]+3*m[6])+29*cos(2*m[6]-2*m[7]);
        l0+=-28*cos(m[6]-m[7])+23*sin(3*m[6])-21*cos(m[4]-m[6])+20*sin(m[6]-m[7])+20*cos(m[5]-m[6])-12*t*t*cos(m[6])-12*cos(m[6])+10*sin(2*m[6]-2*m[7])-9*sin(2*u[6])-9*t*t*sin(m[6])+9*cos(2*m[6]-3*m[7])+8*t*cos(m[5]-2*m[6]);
        l0+=+7*t*cos(m[5]-3*m[6])-7*t*sin(m[5]-3*m[6])+7*t*sin(2*m[6])+6*sin(2*m[4]-6*m[5]+2*m[6])+6*cos(2*m[4]-6*m[5]+2*m[6])+5*sin(m[5]-4*m[6])-4*sin(3*m[6]-4*m[7])+4*cos(3*m[6]-3*m[7])-3*cos(m[7])-2*sin(m[7]);
        l0=l[6]+l0/3600*self._DR;
        b0=2775*sin(u[6])+131*sin(m[6]-u[6])+130*sin(m[6]+u[6]);
        b0=b0/3600*self._DR;
        r0=19.21216-.90154*cos(m[6])-.02488*t*sin(m[6])-.02121*cos(2*m[6])-.00585*cos(m[5]-2*m[6])-.00508*t*cos(m[6])-.00451*cos(m[4]-m[6])+.00336*sin(m[5]-m[6])+.00198*sin(m[4]-m[6])+.00118*cos(m[5]-3*m[6])+.00107*sin(m[5]-2*m[6]);
        r0+=-.00103*t*sin(2*m[6])-.00081*cos(3*m[6]-3*m[7]);
        ll0=l0
        bb0=b0
        rr0=r0
        return ll0, bb0, rr0
        
    def _mc_jd2lbr1g(self, jj, l, m, u):
        """
        /***************************************************************************/
        /* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr  */
        /* pour l'equinoxe moyen de la date. (longitudes vraies)                   */
        /* NEPTUNE                                                                 */
        /***************************************************************************/
        /* Les tableaux l, m et u doivent etre dimensiones chacun avec             */
        /* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a     */
        /***************************************************************************/
        """
        u[0]*=1.;
        t=(jj-2415020.0)/36525;
        l0=3523*sin(m[7])-50*sin(2*u[7])-43*t*cos(m[7])+29*sin(m[4]-m[7])+19*sin(2*m[7])-18*cos(m[4]-m[7])+13*cos(m[5]-m[7])+13*sin(m[5]-m[7])-9*sin(2*m[6]-3*m[7])+9*cos(2*m[6]-2*m[7])-5*cos(2*m[6]-3*m[7]);
        l0+=-4*t*sin(m[7])+4*cos(m[6]-2*m[7])+4*t*t*sin(m[7]);
        l0=l[7]+l0/3600*self._DR;
        b0=6404*sin(u[7])+55*sin(m[7]+u[7])+55*sin(m[7]-u[7])-33*t*sin(u[7]);
        b0=b0/3600*self._DR;
        r0=30.07175-.22701*cos(m[7])-.00787*cos(2*l[6]-m[6]-2*l[7])+.00409*cos(m[4]-m[7])-.00314*t*sin(m[7])+.0025*sin(m[4]-m[7])-.00194*sin(m[5]-m[7])+.00185*cos(m[5]-m[7]);
        ll0=l0
        bb0=b0
        rr0=r0
        return ll0, bb0, rr0

    def _mc_jd2lbr1h(self, jj, l, m, u):
        """
        /***************************************************************************/
        /* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr  */
        /* pour l'equinoxe J2000. ()                                               */
        /* PLUTON                                                                  */
        /***************************************************************************/
        /* Les tableaux l, m et u doivent etre dimensiones chacun avec             */
        /* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a     */
        /* mais ne servent a rien car on se sert de l'algo de Meeus                */
        /* Astronomical Algorithms page 247                                        */
        /***************************************************************************/
        """
        meeus=[
        1,0,0,1,-19798886,19848454,-5453098,-14974876,66867334,68955876,
        2,0,0,2,897499,-4955707,3527363,1672673,-11826086,-333765,
        3,0,0,3,610820,1210521,-1050939,327763,1593657,-1439953,
        4,0,0,4,-341639,-189719,178691,-291925,-18948,482443,
        5,0,0,5,129027,-34863,18763,100448,-66634,-85576,
        6,0,0,6,-38215,31061,-30594,-25838,30841,-5765,
        7,0,1,-1,20349,-9886,4965,11263,-6140,22254,
        8,0,1,0,-4045,-4904,310,-132,4434,4443,
        9,0,1,1,-5885,-3238,2036,-947,-1518,641,
        10,0,1,2,-3812,3011,-2,-674,-5,792,
        11,0,1,3,-601,3468,-329,-563,518,518,
        12,0,2,-2,1237,463,-64,39,-13,-221,
        13,0,2,-1,1086,-911,-94,210,837,-494,
        14,0,2,0,595,-1229,-8,-160,-281,616,
        15,1,-1,0,2484,-485,-177,259,260,-395,
        16,1,-1,1,839,-1414,17,234,-191,-396,
        17,1,0,-3,-964,1059,582,-285,-3218,370,
        18,1,0,-2,-2303,-1038,-298,692,8019,-7869,
        19,1,0,-1,7049,747,157,201,105,45637,
        20,1,0,0,1179,-358,304,825,8623,8444,
        21,1,0,1,393,-63,-124,-29,-896,-801,
        22,1,0,2,111,-268,15,8,208,-122,
        23,1,0,3,-52,-154,7,15,-133,65,
        24,1,0,4,-78,-30,2,2,-16,1,
        25,1,1,-3,-34,-26,4,2,-22,7,
        26,1,1,-2,-43,1,3,0,-8,16,
        27,1,1,-1,-15,21,1,-1,2,9,
        28,1,1,0,-1,15,0,-2,12,5,
        29,1,1,1,4,7,1,0,1,-3,
        30,1,1,3,1,5,1,-1,1,0,
        31,2,0,-6,8,3,-2,-3,9,5,
        32,2,0,-5,-3,6,1,2,2,-1,
        33,2,0,-4,6,-13,-8,2,14,10,
        34,2,0,-3,10,22,10,-7,-65,12,
        35,2,0,-2,-57,-32,0,21,126,-233,
        36,2,0,-1,157,-46,8,5,270,1068,
        37,2,0,0,12,-18,13,16,254,155,
        38,2,0,1,-4,8,-2,-3,-26,-2,
        39,2,0,2,-5,0,0,0,7,0,
        40,2,0,3,3,4,0,1,-11,4,
        41,3,0,-2,-1,-1,0,0,4,-14,
        42,3,0,-1,6,-3,0,0,18,35,
        43,3,0,0,-1,-2,0,1,13,3];
        u[0]*=1.;
        t=(jj-2451545.0)/36525;
        j=34.35+3034.9057*t;
        s=50.08+1222.1138*t;
        p=238.96+144.9600*t;
        l0=0.;
        b0=0.;
        r0=0.;
        for k in range(0,43):
            ij=meeus[k*10+1];
            it=meeus[k*10+2];
            ip=meeus[k*10+3];
            a=(self._DR)*(j*ij+s*it+p*ip);
            sina=sin(a);
            cosa=cos(a);
            A=meeus[k*10+4];
            B=meeus[k*10+5];
            l0+=(A*sina+B*cosa);
            A=meeus[k*10+6];
            B=meeus[k*10+7];
            b0+=(A*sina+B*cosa);
            A=meeus[k*10+8];
            B=meeus[k*10+9];
            r0+=(A*sina+B*cosa);
        l0=(238.956785+144.96*t+l0*1e-6)*self._DR;
        b0=(-3.908202+b0*1e-6)*self._DR;
        r0=40.7247248+r0*1e-7;
        ll0=l0
        bb0=b0
        rr0=r0
        return ll0, bb0, rr0

    def _mc_jd2lbr1d(self, jj):
        """
        /***************************************************************************/
        /* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr  */
        /* pour l'equinoxe J2000.0.                                                */
        /* LUNE                                                                    */
        /***************************************************************************/
        /* J. Meeus, "Astronomical Algorithms", chapter 45, 307-317                */
        /***************************************************************************/
        """
        arg_lr=[
        0,0,1,0,
        2,0,-1,0,
        2,0,0,0,
        0,0,2,0,
        0,1,0,0,
        0,0,0,2,
        2,0,-2,0,
        2,-1,-1,0,
        2,0,1,0,
        2,-1,0,0,
        0,1,-1,0,
        1,0,0,0,
        0,1,1,0,
        2,0,0,-2,
        0,0,1,2,
        0,0,1,-2,
        4,0,-1,0,
        0,0,3,0,
        4,0,-2,0,
        2,1,-1,0,
        2,1,0,0,
        1,0,-1,0,
        1,1,0,0,
        2,-1,1,0,
        2,0,2,0,
        4,0,0,0,
        2,0,-3,0,
        0,1,-2,0,
        2,0,-1,2,
        2,-1,-2,0,
        1,0,1,0,
        2,-2,0,0,
        0,1,2,0,
        0,2,0,0,
        2,-2,-1,0,
        2,0,1,-2,
        2,0,0,2,
        4,-1,-1,0,
        0,0,2,2,
        3,0,-1,0,
        2,1,1,0,
        4,-1,-2,0,
        0,2,-1,0,
        2,2,-1,0,
        2,1,-2,0,
        2,-1,0,-2,
        4,0,1,0,
        0,0,4,0,
        4,-1,0,0,
        1,0,-2,0,
        2,1,0,-2,
        0,0,2,-2,
        1,1,1,0,
        3,0,-2,0,
        4,0,-3,0,
        2,-1,2,0,
        0,2,1,0,
        1,1,-1,0,
        2,0,3,0,
        2,0,-1,-2]
        
        sinl=[
        6288774,
        1274027,
        658314,
        213618,
        -185116,
        -114332,
        58793,
        57066,
        53322,
        45758,
        -40923,
        -34720,
        -30383,
        15327,
        -12528,
        10980,
        10675,
        10034,
        8548,
        -7888,
        -6766,
        -5163,
        4987,
        4036,
        3994,
        3861,
        3665,
        -2689,
        -2602,
        2390,
        -2348,
        2236,
        -2120,
        -2069,
        2048,
        -1773,
        -1595,
        1215,
        -1110,
        -892,
        -810,
        759,
        -713,
        -700,
        691,
        596,
        549,
        537,
        520,
        -487,
        -399,
        -381,
        351,
        -340,
        330,
        327,
        -323,
        299,
        294,
        0]
        
        cosr=[
        -20905355,
        -3699111,
        -2955968,
        -569925,
        48888,
        -3149,
        246158,
        -152138,
        -170733,
        -204586,
        -129620,
        108743,
        104755,
        10321,
        0,
        79661,
        -34782,
        -23210,
        -21636,
        24208,
        30824,
        -8379,
        -16675,
        -12831,
        -10445,
        -11650,
        14403,
        -7003,
        0,
        10056,
        6322,
        -9884,
        5751,
        0,
        -4950,
        4130,
        0,
        -3958,
        0,
        3258,
        2616,
        -1897,
        -2117,
        2354,
        0,
        0,
        -1423,
        -1117,
        -1571,
        -1739,
        0,
        -4421,
        0,
        0,
        0,
        0,
        1165,
        0,
        0,
        8752]
        
        arg_b=[
        0,0,0,1,
        0,0,1,1,
        0,0,1,-1,
        2,0,0,-1,
        2,0,-1,1,
        2,0,-1,-1,
        2,0,0,1,
        0,0,2,1,
        2,0,1,	-1,
        0,0,2,	-1,
        2,-1,0,-1,
        2,0,-2,-1,
        2,0,1,1,
        2,1,0,-1,
        2,-1,-1,1,
        2,-1,0,1,
        2,-1,-1,-1,
        0,1,-1,-1,
        4,0,-1,-1,
        0,1,0,1,
        0,0,0,3,
        0,1,-1,1,
        1,0,0,1,
        0,1,1,1,
        0,1,1,-1,
        0,1,0,-1,
        1,0,0,-1,
        0,0,3,1,
        4,0,0,-1,
        4,0,-1,1,
        0,0,1,-3,
        4,0,-2,1,
        2,0,0,-3,
        2,0,2,-1,
        2,-1,1,-1,
        2,0,-2,1,
        0,0,3,-1,
        2,0,2,1,
        2,0,-3,-1,
        2,1,-1,1,
        2,1,0,1,
        4,0,0,1,
        2,-1,1,1,
        2,-2,0,-1,
        0,0,1,3,
        2,1,1,-1,
        1,1,0,-1,
        1,1,0,1,
        0,1,-2,-1,
        2,1,-1,-1,
        1,0,1,1,
        2,-1,-2,-1,
        0,1,2,1,
        4,0,-2,-1,
        4,-1,-1,-1,
        1,0,1,-1,
        4,0,1,-1,
        1,0,-1,-1,
        4,-1,0,-1,
        2,-2,0,1]
        
        sinb=[
        5128122,
        280602,
        277693,
        173237,
        55413,
        46271,
        32573,
        17198,
        9266,
        8822,
        8216,
        4324,
        4200,
        -3359,
        2463,
        2211,
        2065,
        -1870,
        1828,
        -1794,
        -1749,
        -1565,
        -1491,
        -1475,
        -1410,
        -1344,
        -1335,
        1107,
        1021,
        833,
        777,
        671,
        607,
        596,
        491,
        -451,
        439,
        422,
        421,
        -366,
        -351,
        331,
        315,
        302,
        -283,
        -229,
        223,
        223,
        -220,
        -220,
        -185,
        181,
        -177,
        176,
        166,
        -164,
        132,
        -119,
        115,
        107]
        
        T=(jj-2451545.)/36525.;
        lp=(218.3164591+481267.88134236*T-.0013268*T*T+T*T*T/538841.-T*T*T*T/65194000)*(self._DR);
        d=(297.8502042+445267.1115168*T-.00016300*T*T+T*T*T/545868-T*T*T*T/113065000)*(self._DR);
        m=(357.5291092+35999.0502909*T-.0001536*T*T+T*T*T*T/24490000)*(self._DR);
        mp=(134.9634114+477198.8676313*T+.0089970*T*T+T*T*T/69699.-T*T*T*T/14712000)*(self._DR);
        f=(93.2720993+483202.0175273*T-.0034029*T*T+T*T*T/3526000+T*T*T*T/863310000)*(self._DR);
        a1=(119.75+131.849*T)*(self._DR);
        a2=(53.09+479264.290*T)*(self._DR);
        a3=(313.45+481266.484*T)*(self._DR);
        e=1-.002516*T-.0000074*T*T;
        e2=e*e;
        
        # --- longitude & radius ---
        l=0.;
        r=0.;
        for k in range(0,60):
            xe=1.;
            if (fabs(arg_lr[k*4+1])==1):
                xe=e;
            elif (fabs(arg_lr[k*4+1])==2):
                xe=e2
            angle=1.*arg_lr[k*4+0]*d+arg_lr[k*4+1]*m+arg_lr[k*4+2]*mp+arg_lr[k*4+3]*f
            sina=sin(angle)
            cosa=cos(angle)
            l+=sinl[k]*sina*xe
            r+=cosr[k]*cosa*xe
        l+=3958.*sin(a1)+1962.*sin(lp-f)+318.*sin(a2)
        
        # --- latitude ---
        b=0.
        for k in range(0,60):
            xe=1.;
            if (fabs(arg_b[k*4+1])==1):
                xe=e
            elif (fabs(arg_b[k*4+1])==2):
                xe=e2
            angle=1.*arg_b[k*4+0]*d+arg_b[k*4+1]*m+arg_b[k*4+2]*mp+arg_b[k*4+3]*f
            sina=sin(angle)
            b+=sinb[k]*sina*xe
        b+=-2235.*sin(lp)+382.*sin(a3)+175.*sin(a1-f)+175.*sin(a1+f)+127.*sin(lp-mp)-115.*sin(lp+mp)
        
        l=lp+(l*1.0e-6)*(self._DR)
        b=(b*1.0e-6)*(self._DR)
        r=(385000.56e3+r)/(self._CST_UA)
        
        ll0=l;
        bb0=b;
        rr0=r;
        return ll0, bb0, rr0


# ========================================================
# === internal methods : Coordinate transformations
# ========================================================
        
    def _mc_obliqmoy(self, jj):
        """
        /***************************************************************************/
        /* Retourne la valeur de l'obliquite terrestre moyenne pour jj             */
        /***************************************************************************/
        /* formule de Laskar (JM)                                                  */
        /***************************************************************************/
        """
        t=(jj-2451545.0)/36525
        u=t/100
        eps0=u*(-4680.93-u*(1.55+u*(1999.25-u*(51.38-u*(249.67-u*(39.05+u*(7.12+u*(27.87+u*(5.79+u*(2.45))))))))))
        eps0=(23.4392911111+eps0/3600)*self._DR
        return eps0

    def _mc_xyzec2eq(self, xec, yec, zec, eps):
        """
        /***************************************************************************/
        /* Transforme les coord. cart. ecliptiques vers equatoriales               */
        /***************************************************************************/
        /***************************************************************************/
        """
        xeq0=xec
        yeq0=yec*cos(eps)-zec*sin(eps)
        zeq0=yec*sin(eps)+zec*cos(eps)
        return xeq0, yeq0, zeq0

    def _mc_xyzeq2ec(self, xeq, yeq, zeq, eps):
        """
        /***************************************************************************/
        /* Transforme les coord. cart. equatoriales vers ecliptiques               */
        /***************************************************************************/
        /***************************************************************************/
        """
        eps=-eps
        xec0=xeq
        yec0=yeq*cos(eps)-zeq*sin(eps)
        zec0=yeq*sin(eps)+zeq*cos(eps)
        return xec0, yec0, zec0

    def _mc_paraldxyzeq(self, jj, longuai, rhocosphip, rhosinphip):
        """
        /***************************************************************************/
        /* Calcul des corrections cartesiennes equatoriales de la parallaxe        */
        /* Xtopo = Xgeo - *dxeq etc...                                             */
        /***************************************************************************/
        /* ref : Danby J.M.A. "Fundamentals of Celestial Mechanics" 2nd ed. 1992   */
        /*       Willmann Bell                                                     */
        /* Formules (6.17.1) pp 208                                                */
        /***************************************************************************/
        """
        cst=(self._CST_EARTH_SEMI_MAJOR_RADIUS)/(self._CST_UA); # equatorial radius of the Earth in U.A.
        tsl = self._mc_tsl(jj,-longuai)
        dxeq=(cst*rhocosphip*cos(tsl))
        dyeq=(cst*rhocosphip*sin(tsl))
        dzeq=(cst*rhosinphip)
        return dxeq, dyeq, dzeq

    def _mc_latalt2rhophi(self, latitude, altitude):
        """
        /***************************************************************************/
        /* Retourne les valeurs de rhocosphi' et rhosinphi' (en rayons equatorial  */
        /* terrestres) a partir de la latitude et de l'altitude.                   */
        /***************************************************************************/
        /* Latitude en degres decimaux.                                            */
        /* Altitude en metres.                                                     */
        /* Algo : Meeus "Astronomical Algorithms" p78                              */
        /***************************************************************************/
        """
        lat=latitude*(self._DR);
        alt=altitude;
        aa=self._CST_EARTH_SEMI_MAJOR_RADIUS;
        ff=1./self._CST_EARTH_INVERSE_FLATTENING;
        bb=aa*(1-ff);
        u=atan(bb/aa*tan(lat));
        a=bb/aa*sin(u)+alt/aa*sin(lat);
        b=      cos(u)+alt/aa*cos(lat);
        rhocosphip=b;
        rhosinphip=a;
        return rhocosphip, rhosinphip

    def _mc_rhophi2latalt(self, rhosinphip,rhocosphip):
        """
        /***************************************************************************/
        /* Retourne les valeurs de la latitude et de l'altitude a partir de        */
        /* rhocosphi' et rhosinphi' (en rayons equatorial terrestres)              */
        /***************************************************************************/
        /* Latitude en degres decimaux.                                            */
        /* Altitude en metres.                                                     */
        /* Algo : Meeus "Astronomical Algorithms" p78                              */
        /***************************************************************************/
        """
        aa=self._CST_EARTH_SEMI_MAJOR_RADIUS;
        ff=1./self._CST_EARTH_INVERSE_FLATTENING;
        bb=aa*(1-ff);
        rho=sqrt(rhosinphip*rhosinphip+rhocosphip*rhocosphip);
        if (rho==0.):
            latitude=0.;
            altitude=-aa;
            return latitude, altitude
        phip=atan2(rhosinphip,rhocosphip)
        phi0=atan(aa*aa/bb/bb*tan(phip));  # alt=0
        u0=atan(bb/aa*tan(phi0));
        sinu0=sin(u0);
        cosu0=cos(u0);
        sinphi0=sin(phi0);
        cosphi0=cos(phi0);
        alt=-1000
        while (alt<20000.):
            rhosinphip0 = bb/aa*sinu0 + alt/aa*sinphi0 ;
            rhocosphip0 =       cosu0 + alt/aa*cosphi0 ;
            rho0=sqrt(rhosinphip0*rhosinphip0+rhocosphip0*rhocosphip0);
            if ((rho-rho0)<0):
                break;
            alt += 0.1
        lat=phi0;
        alt-=0.1;
        latitude=lat/(self._DR);
        altitude=alt;
        return latitude, altitude

    def _mc_nutation(self, jj, precision):
        """    
        /***************************************************************************/
        /* Retourne la valeur de la nutation pour jj donne                         */
        /***************************************************************************/
        /* precision=0 pour une incertitude de 0"5 sur dpsi 0"1 sur deps           */
        /* precision=1 pour de la tres haute precision (environ 0.01")             */
        /***************************************************************************/
        """
        t=(jj-2451545.0)/36525;
        o=125.04452-1934.136261*t+.0020708*t*t+t*t*t/450000;
        o=fmod(o*self._DR,2*self._PI);
        if (precision==0):
            l=280.4665+36000.7698*t;
            lp=218.3165+481267.8813*t;
            l=fmod(l*self._DR,2*self._PI);
            lp=fmod(lp*self._DR,2*self._PI);
            dpsi0=-17.20*sin(o)-1.32*sin(2*l)-.23*sin(2*lp)+.21*sin(2*o);
            dpsi0=dpsi0/3600*self._DR;
            deps0=9.20*cos(o)+.57*cos(2*l)+.10*cos(2*lp)-.09*cos(2*o);
            deps0=deps0/3600*self._DR;
        else: 
            # if (precision==1)
            d=297.85036+445267.111480*t-.0019142*t*t+t*t*t/189474
            m=357.52772+35999.050340*t-.0001603*t*t-t*t*t/300000
            mp=134.96298+477198.867398*t+.0086972*t*t+t*t*t/56250
            f=93.27191+483202.017538*t-.0036825*t*t+t*t*t/327270
            d=fmod(d*self._DR,2*self._PI)
            m=fmod(m*self._DR,2*self._PI)
            mp=fmod(mp*self._DR,2*self._PI)
            f=fmod(f*self._DR,2*self._PI)
            dpsi0=(-171996-174.2*t)*sin(o)+(-13187-1.6*t)*sin(-2*d+2*f+2*o)
            +(-2274-.02*t)*sin(2*f+2*o)+(2062+.2*t)*sin(2*o)
            +(1426-3.4*t)*sin(m)+(712+.1*t)*sin(mp)
            dpsi0+=((-517+1.2*t)*sin(-2*d+m+2*f+2*o)-(386-.4*t)*sin(2*f+o)
            -301*sin(mp+2*f+2*o)+(217-.5*t)*sin(-2*d-m+2*f+2*o)
            -158*sin(-2*d+mp)+(129+.1*t)*sin(-2*d+2*f+o))
            dpsi0+=(123*sin(-mp+2*f+2*o)+63*sin(2*d)+(63+.1*t)*sin(mp+o)
            -59*sin(2*d-mp+2*f+2*o)+(-58-.1*t)*sin(-mp+o)-51*sin(mp+2*f+o)
            +48*sin(-2*d+2*mp)+46*sin(-2*mp+2*f+o)-38*sin(2*d+2*f+2*o))
            dpsi0=dpsi0*1e-4/3600*self._DR
            deps0=(92025+8.9*t)*cos(o)+(5736-3.1*t)*cos(-2*d+2*f+2*o)
            +(977-.5*t)*cos(2*f+2*o)+(-895+.5*t)*cos(2*o)+(54-.1*t)*cos(m)
            -7*cos(mp)+(224-.6*t)*cos(-2*d+m+2*f+2*o)+200*cos(2*f+o);
            deps0+=((129-.1*t)*cos(mp+2*f+2*o)+(-95+.3*t)*cos(-2*d-m+2*f+2*o)
            -70*cos(-2*d+2*f+o)-53*cos(-mp+2*f+2*o)-33*cos(mp+o)
            +26*cos(2*d-mp+2*f+2*o)+32*cos(-mp+o)+27*cos(mp+2*f+o))
            deps0=deps0*1e-4/3600*self._DR
        dpsi=dpsi0;
        deps=deps0;
        return dpsi, deps          

    def _mc_nutradec(self, jj, asd1, dec1, signe):
        """
        /***************************************************************************/
        /* Corrige asd1,dec1 de la nutation et retourne asd2 et dec2               */
        /***************************************************************************/
        /* Trueblood & Genet : Telescop Control ed. Willmann Bell (1997) p71       */
        /***************************************************************************/
        """
        methode=1;
        
        # --- obliquite moyenne --- */
        eps = self._mc_obliqmoy(jj)

        # --- nutation ---*/
        dpsi, deps = self._mc_nutation(jj,1)

        if (fabs(dec1) >= (self._PI)/2.):
            asd2=asd1
            dec2=dec1
        else:
            if (methode==0):
                tand=tan(dec1);
                dasd=(cos(eps)+sin(eps)*sin(asd1)*tand)*dpsi-cos(asd1)*tand*deps;
                ddec=sin(eps)*cos(asd1)*dpsi+sin(asd1)*deps;
                dasd*=signe;
                ddec*=signe;
                asd1+=dasd;
                dec1+=ddec;
            else:
                # eq->ecl
                l=asd1 ; b=dec1 ; r=1;                
                bb = Angle(b/self._DR)
                ll = Angle(l/self._DR)
                coords = Coords((r,ll,bb))
                xeq,yeq,zeq = coords.cart()
                xec,yec,zec = self._mc_xyzeq2ec(xeq,yeq,zeq,eps)
                coords = Coords((xec,yec,zec))
                r,l,b = coords.sphe("rad","rad")
                l=l+dpsi*signe;
                eps=eps+deps*signe;
                bb = Angle(b/self._DR)
                ll = Angle(l/self._DR)
                coords = Coords((r,ll,bb))
                xeq,yeq,zeq = coords.cart()
                # ecl->eq
                xeq,yeq,zeq = self._mc_xyzec2eq(xec,yec,zec,eps)
                coords = Coords((xeq,yeq,zeq))
                r,asd1,dec1 = coords.sphe("rad","rad")
            asd1=fmod(4*self._PI+asd1,2*self._PI);
            asd2=asd1;
            dec2=dec1;
        return asd2, dec2


    def _mc_he2ge(self, xh,yh,zh,xs,ys,zs):
        """
        /***************************************************************************/
        /* Translation du repere heliocentrique en geocentrique (cartesien)        */
        /***************************************************************************/
        /* ENTREES :                                                               */
        /* xh,yh,zh : coordonnees heliocentriques de l'astre                       */
        /* xs,ys,zs : coordonnees geocentriques du Soleil                          */
        /* SORTIES :                                                               */
        /* *xg,*yg,*zg : coordonnees geocentriques de l'astre                      */
        /***************************************************************************/
        """
        xg=xh+xs
        yg=yh+ys
        zg=zh+zs
        return xg, yg, zg

    def _mc_aberpla(self, jj1, delta):
        """
        /***************************************************************************/
        /* Corrige jj1 de l'aberration de la lumiere et retourne jj2               */
        /***************************************************************************/
        /***************************************************************************/
        """
        jj2=jj1-0.0057755*delta
        return jj2
    
    def _mc_parallaxe_stellaire(self, jj, asd1, dec1, plx_mas):
        """
        /***************************************************************************/
        /* Corrige asd1,dec1 de la parallaxe stellaire et retourne asd2 et dec2    */
        /***************************************************************************/
        /* A. Danjon : Astronomie Generale ed. A. Blanchard (1980) p130            */
        /***************************************************************************/
        """
        afm = self._AFM_VFP79
        # --- obliquite moyenne --- */
        eps = self._mc_obliqmoy(jj)        
        # --- longitude vraie du soleil ---*/
        llp, mmp, uup = self._mc_jd2lbr1a(jj);
        coords = self._mc_jd2lbr1b(jj, self._PLN_SOLEIL, llp, mmp, uup, afm)
        rs,ls,bs = coords.sphe("rad","rad")        
        dpsi, deps = self._mc_nutation(jj,1)
        ls += dpsi
        plxrad=plx_mas*1e-3/3600.*(self._DR);
        secd=cos(dec1);
        if (secd==0):
            asd2=asd1;
            dec2=dec1;
        else:
            secd=1./secd;
            dasd=(cos(eps)*cos(asd1)*sin(ls)-sin(asd1)*cos(ls))*secd;
            ddec=(sin(eps)*cos(dec1)*sin(ls)-sin(dec1)*cos(asd1)*cos(ls)-cos(eps)*sin(dec1)*sin(asd1)*sin(ls));
            asd1+=plxrad*dasd;
            dec1+=plxrad*ddec;
            asd1=fmod(4*self._PI+asd1,2*self._PI);
            asd2=asd1;
            dec2=dec1;
        return asd2, dec2

    def _mc_aberration_annuelle(self, jj, asd1, dec1, signe):
        """
        /***************************************************************************/
        /* Corrige asd1,dec1 de l'aberration annuelle et retourne asd2 et dec2     */
        /***************************************************************************/
        /* Trueblood & Genet : Telescop Control ed. Willmann Bell (1997) p81-82    */
        /* Formule sans les E-terms                                                */
        /***************************************************************************/
        """
        k=20.49552; # constant of annual aberration        
        # --- computation method
        afm = self._AFM_VFP79
        # --- obliquite moyenne ---
        eps = self._mc_obliqmoy(jj)    
        # --- longitude vraie du soleil ---
        llp, mmp, uup = self._mc_jd2lbr1a(jj);
        coords = self._mc_jd2lbr1b(jj, self._PLN_SOLEIL, llp, mmp, uup, afm)
        rs,ls,bs = coords.sphe("rad","rad")
        dpsi, deps = self._mc_nutation(jj,1)
        ls += dpsi;
        # ---
        secd=cos(dec1)
        if (secd==0):
            asd2=asd1
            dec2=dec1
            return asd2, dec2
        secd=1./secd
        if (fabs(secd) < 100):
            c = cos(asd1)*secd;
            d = sin(asd1)*secd;
            cp = tan(eps)*cos(dec1)-sin(asd1)*sin(dec1);
            dp = cos(asd1)*sin(dec1);
            cc = -k*cos(eps)*cos(ls);
            dd = -k*sin(ls);
            dasd = (cc*c+dd*d)/3600.*(self._DR);
            ddec = (cc*cp+dd*dp)/3600.*(self._DR);
            dasd *= float(signe);
            ddec *= float(signe);
            asd1 += dasd;
            dec1 += ddec;
        asd1 = fmod(4*self._PI+asd1, 2*self._PI);
        asd2 = asd1;
        dec2 = dec1;
        return asd2, dec2

    def _mc_aberration_diurne(self, jj, asd1, dec1, longuai, rhocosphip, rhosinphip, signe):
        """
        /***************************************************************************/
        /* Corrige asd1,dec1 de l'aberration diurne et retourne asd2 et dec2       */
        /***************************************************************************/
        /* Trueblood & Genet : Telescop Control ed. Willmann Bell (1997) p83-84    */
        /***************************************************************************/
        """
        a=(self._CST_EARTH_SEMI_MAJOR_RADIUS)*1e-3;
        tsl = self._mc_tsl(jj,-longuai)
        h=tsl-asd1;
        latitude, altitude = self._mc_rhophi2latalt(rhosinphip, rhocosphip)
        if (rhocosphip==0.):
            if (rhosinphip>0):
                phip=(self._PI)/2.
            else:
                phip=-(self._PI)/2.
        else:
            phip=atan2(rhosinphip,rhocosphip)         
        sinphi=sin(latitude)
        r=(21*sinphi*sinphi+a)/a;
        secd=cos(dec1);
        if (secd==0):
            asd2=asd1
            dec2=dec1
        else:
            secd=1./secd;
            if (fabs(secd)<100):
                dasd=(0.320*r*cos(phip)*cos(h)*secd)/3600.*(self._DR);
                ddec=(0.320*r*cos(phip)*sin(h)*sin(dec1))/3600.*(self._DR);
            else:
                dasd=0
                ddec=0
            dasd*=signe;
            ddec*=signe;
            asd1+=dasd;
            dec1+=ddec;
            asd1=fmod(4*self._PI+asd1,2*self._PI);
        asd2=asd1;
        dec2=dec1;
        return asd2, dec2

    def _mc_precad(self, jd1, asd1, dec1, jd2):
        """
        /***************************************************************************/
        /* Passage des coordonnees spheri. equatoriales d'un equinoxe a un autre   */
        /***************************************************************************/
        /***************************************************************************/
        """
        t=(jd2-jd1)/36525.
        tt=(jd1-2451545.0)/36525.
        dz=(2306.2181+1.39656*tt-.000139*tt*tt)*t+(0.30188-.000344*tt)*t*t+.017998*t*t*t
        dz=dz/3600*self._DR
        zz=(2306.2181+1.39656*tt-.000139*tt*tt)*t+(1.09468+.000066*tt)*t*t+.018203*t*t*t
        zz=zz/3600*self._DR
        th=(2004.3109-0.85330*tt-.000217*tt*tt)*t-(0.42665+.000217*tt)*t*t-.041833*t*t*t
        th=th/3600*self._DR
        cosasddz=cos(asd1+dz)
        sinasddz=sin(asd1+dz)
        costh=cos(th)
        sinth=sin(th)
        cosdec=cos(dec1)
        sindec=sin(dec1)
        a=cosdec*sinasddz
        b=costh*cosdec*cosasddz-sinth*sindec
        c=sinth*cosdec*cosasddz+costh*sindec
        a=atan2(a,b)+zz
        asd2=fmod(4*self._PI+a,2*self._PI)
        dec2=asin(c)
        return asd2, dec2    
    
    def _mc_elonphas(self, r, rsol, delta):
        """
        /***************************************************************************/
        /* Calcul des angles d'elongation et de phase.                             */
        /***************************************************************************/
        /***************************************************************************/
        """
        elong=0.
        phase=0.
        if (delta!=0.):
            if (rsol!=0.):
                elong=acos((rsol*rsol+delta*delta-r*r)/(2*rsol*delta))
            if (r!=0.):
                phase=acos((r*r+delta*delta-rsol*rsol)/(2*r*delta))
        return elong, phase    

    def _mc_elonphaslimb(self, asd, dec, asds, decs, r, delta):
        """    
        /***************************************************************************/
        /* Calcul des angles d'elongation de phase et de position du limbe.        */
        /***************************************************************************/
        /* Meeus page 316 (46.2) (46.3) (46.5)                                     */
        /* asd,dec : planete.                                                      */
        /* asds,decs : soleil                                                      */
        /* r : distance Terre-Soleil (corrigee de l'abberation).                   */
        /* delta : distance Terre-Planete (corrigee de l'abberation).              */
        /***************************************************************************/
        """
        elong=0.
        phase=0.
        elong=acos(sin(decs)*sin(dec)+cos(decs)*cos(dec)*cos(asds-asd))
        i=atan2(r*sin(elong),delta-r*cos(elong))
        if (i<0):
            i=i+self._PI
        phase=i;
        i=atan2(cos(decs)*sin(asds-asd),sin(decs)*cos(dec)-cos(decs)*sin(dec)*cos(asds-asd));
        if (i<0):
            i=2*self._PI+i
        posang_brightlimb=fmod(4*self._PI+i,2*self._PI)
        return elong, phase, posang_brightlimb

    def _mc_hd2ah(self, ha, dec, latitude):
        """
        /***************************************************************************/
        /* Transforme les coord. sph. equatoriales vers sph. azinuth hauteur       */
        /***************************************************************************/
        /* ha est l'angle horaire local.                                           */
        /* Tout en radian                                                          */
        /***************************************************************************/
        """
        if (dec>=self._PISUR2):
            aa=(self._PI)
            hh=latitude
        elif (dec<=-self._PISUR2):
            aa=0.
            hh=-latitude
        else:
            aa=atan2(sin(ha),cos(ha)*sin(latitude)-tan(dec)*cos(latitude))
            hh=asin(sin(latitude)*sin(dec)+cos(latitude)*cos(dec)*cos(ha))
        az=fmod(4*self._PI+aa,2*self._PI);
        h=hh;
        return az, h

    def _mc_ad2hd(self, jd, longuai, asd):
        """
        /***************************************************************************/
        /* Transforme l'ascension droite en angle horaire                          */
        /***************************************************************************/
        /***************************************************************************/
        """
        # --- calcul du TSL en radians ---*/
        tsl = self._mc_tsl(jd,-longuai)
        h=tsl-asd;
        ha=fmod(4*self._PI+h,2*self._PI);
        return ha

    def _mc_ah2hd(self, az, h, latitude):
        """
        /***************************************************************************/
        /* Transforme les coord. sph. azinuth hauteur vers sph. equatoriales       */
        /***************************************************************************/
        /* ha est l'angle horaire local.                                           */
        /***************************************************************************/
        """
        if (h>=self._PISUR2):
            ahh=0.
            decc=latitude
        elif (h<=-self._PISUR2):
            ahh=0.
            decc=-latitude
        else:
            ahh=atan2(sin(az),cos(az)*sin(latitude)+tan(h)*cos(latitude))
            decc=asin(sin(latitude)*sin(h)-cos(latitude)*cos(h)*cos(az))
        ha=fmod(4*self._PI+ahh,2*self._PI)
        dec=decc
        return ha, dec

    def _mc_hd2ad(self, jd, longuai, ha):
        """
        /***************************************************************************/
        /* Transforme l'angle horaire en ascension droite                          */
        /***************************************************************************/
        /***************************************************************************/
        """
        #--- calcul du TSL en radians ---*/
        tsl = self._mc_tsl(jd,-longuai)
        ra=tsl-ha
        asd=fmod(4*self._PI+ra,2*self._PI)
        return asd
    
    def _mc_hd2parallactic(self, ha, dec, latitude):
        """
        /***************************************************************************/
        /* Transforme les coord. sph. equatoriales vers angle parallactic          */
        /***************************************************************************/
        /* ha est l'angle horaire local.                                           */
        /***************************************************************************/
        """
        if (fabs(latitude)>=self._PISUR2):
            aa=0
        elif ((ha==0.) and (dec==latitude)):
            aa=0.
        else:
            aa=atan2(sin(ha),cos(dec)*tan(latitude)-sin(dec)*cos(ha))
        parallactic=aa;
        return parallactic

# ========================================================
# === internal methods : Refraction
# ========================================================

# ========================================================
# === internal methods : Time
# ========================================================
        
    def _mc_tsl(self, jj, longitude):
        """
        /***************************************************************************/
        /* Calcul du temps sideral local apparent (en radian)                      */
        /* La longitude est comptee en radian positive vers l'ouest                */
        /* jj = UT1  (on a toujours |UTC-UT1|<1 sec)                               */
        /***************************************************************************/
        /* Formules (11.1) et (11.4) de Meeus                                      */
        /***************************************************************************/
        """
        j=(jj-2451545.0)
        t=j/36525
        theta0=280.460618375+360.98564736629*j+.000387933*t*t-t*t*t/38710000        
        eps = self._mc_obliqmoy(jj)
        dpsi, deps = self._mc_nutation(jj,1)
        theta0+=(dpsi*cos(eps+deps)/(self._DR))
        theta0-=longitude/(self._DR)
        theta0=fmod(theta0,360.)
        theta0=fmod(theta0+720.,360.)*self._DR
        tsl=theta0
        return tsl

    def _mc_td2tu(self, jjtd):
        """
        /***************************************************************************/
        /* Retourne la valeur de jj en TU a partir de jj en temps dynamique        */
        /***************************************************************************/
        /* Algo : Meeus "Astronomical Algorithms" p73 (9.1)                        */
        /***************************************************************************/
        """
        dt=0.
        dt = self._mc_tdminusut(jjtd)
        jjtu=jjtd-dt/86400.
        return jjtu


    def _mc_tu2td(self, jjtu):
        """
        /***************************************************************************/
        /* Retourne la valeur de jj en temps dynamique a partir de jj en TU        */
        /* UTC -> TT                                                               */
        /***************************************************************************/
        /* Algo : Meeus "Astronomical Algorithms" p73 (9.1)                        */
        /***************************************************************************/
        """
        dt=0.;
        dt = self._mc_tdminusut(jjtu)
        jjtd=jjtu+dt/86400.
        return jjtd

    def _mc_tdminusut(self, jj):
        """
        /***************************************************************************/
        /* Retourne la valeur dt(sec)=TT-UT partir de jj en TU                     */
        /***************************************************************************/
        /* Algo : Meeus "Astronomical Algorithms" p72                              */
        /* and ftp://maia.usno.navy.mil/ser7/tai-utc.dat                           */
        /*                                                ET 1960-1983             */
        /*                                                TDT 1984-2000            */
        /* UTC 1972-  GPS 1980-    TAI 1958-               TT 2001-                */
        /*----+---------+-------------+-------------------------+-----             */
        /*    |         |             |                         |                  */
        /*    |<------ TAI-UTC ------>|<-----   TT-TAI    ----->|                  */
        /*    |         |             |      32.184s fixed      |                  */
        /*    |<GPS-UTC>|<- TAI-GPS ->|                         |                  */
        /*    |         |  19s fixed  |                         |                  */
        /*    |                                                 |                  */
        /*    <> delta-UT = UT1-UTC                             |                  */
        /*     | (max 0.9 sec)                                  |                  */
        /*-----+------------------------------------------------+-----             */
        /*     |<-------------- delta-T = TT-UT1 -------------->|                  */
        /*    UT1 (UT)                                       TT/TDT/ET             */
        /*                                                                         */
        /* http://stjarnhimlen.se/comp/time.html                                   */
        /***************************************************************************/
        """
        ds=0.
        t=0.
        indexmax=201
        table = []
        table.append(2312752.5) ; table.append(  +124.) ; # 1620 
        table.append(2313483.5) ; table.append(   +115.)
        table.append(2314213.5) ; table.append(   +106.)
        table.append(2314944.5) ; table.append(   +98.)
        table.append(2315674.5) ; table.append(   +91.)
        table.append(2316405.5) ; table.append(   +85.)
        table.append(2317135.5) ; table.append(   +79.)
        table.append(2317866.5) ; table.append(   +74.)
        table.append(2318596.5) ; table.append(   +70.)
        table.append(2319327.5) ; table.append(   +65.)
        table.append(2320057.5) ; table.append(   +62.)
        table.append(2320788.5) ; table.append(   +58.)
        table.append(2321518.5) ; table.append(   +55.)
        table.append(2322249.5) ; table.append(   +53.)
        table.append(2322979.5) ; table.append(   +50.)
        table.append(2323710.5) ; table.append(   +48.)
        table.append(2324440.5) ; table.append(   +46.)
        table.append(2325171.5) ; table.append(   +44.)
        table.append(2325901.5) ; table.append(   +42.)
        table.append(2326632.5) ; table.append(   +40.)
        table.append(2327362.5) ; table.append(   +37.)
        table.append(2328093.5) ; table.append(   +35.)
        table.append(2328823.5) ; table.append(   +33.)
        table.append(2329554.5) ; table.append(   +31.)
        table.append(2330284.5) ; table.append(   +28.)
        table.append(2331015.5) ; table.append(   +26.)
        table.append(2331745.5) ; table.append(   +24.)
        table.append(2332476.5) ; table.append(   +22.)
        table.append(2333206.5) ; table.append(   +20.)
        table.append(2333937.5) ; table.append(   +18.)
        table.append(2334667.5) ; table.append(   +16.)
        table.append(2335398.5) ; table.append(   +14.)
        table.append(2336128.5) ; table.append(   +13.)
        table.append(2336859.5) ; table.append(   +12.)
        table.append(2337589.5) ; table.append(   +11.)
        table.append(2338320.5) ; table.append(   +10.)
        table.append(2339050.5) ; table.append(    +9.)
        table.append(2339781.5) ; table.append(    +9.)
        table.append(2340511.5) ; table.append(    +9.)
        table.append(2341242.5) ; table.append(    +9.)
        table.append(2341972.5) ; table.append(     9.)   ; # 1700 
        table.append(2342702.5) ; table.append(     9.)   ; # 1702 
        table.append(2343432.5) ; table.append(     9.)   ; # 1704 
        table.append(2344163.5) ; table.append(     9.)   ; # 1706 
        table.append(2344893.5) ; table.append(    10.)   ; # 1708 
        table.append(2345624.5) ; table.append(    10.)   ; # 1710 
        table.append(2346354.5) ; table.append(    10.)   ; # 1712 
        table.append(2347085.5) ; table.append(    10.)   ; # 1714 
        table.append(2347815.5) ; table.append(    10.)   ; # 1716 
        table.append(2348546.5) ; table.append(    11.)   ; # 1718 
        table.append(2349276.5) ; table.append(    11.)   ; # 1720 
        table.append(2350007.5) ; table.append(    11.)   ; # 1722 
        table.append(2350737.5) ; table.append(    11.)   ; # 1724 
        table.append(2351468.5) ; table.append(    11.)   ; # 1726 
        table.append(2352198.5) ; table.append(    11.)   ; # 1728 
        table.append(2352929.5) ; table.append(    11.)   ; # 1730 
        table.append(2353659.5) ; table.append(    11.)   ; # 1732 
        table.append(2354390.5) ; table.append(    12.)   ; # 1734 
        table.append(2355120.5) ; table.append(    12.)   ; # 1736 
        table.append(2355851.5) ; table.append(    12.)   ; # 1738 
        table.append(2356581.5) ; table.append(    12.)   ; # 1740 
        table.append(2357312.5) ; table.append(    12.)   ; # 1742 
        table.append(2358042.5) ; table.append(    13.)   ; # 1744 
        table.append(2358773.5) ; table.append(    13.)   ; # 1746 
        table.append(2359503.5) ; table.append(    13.)   ; # 1748 
        table.append(2360234.5) ; table.append(    13.)   ; # 1750 
        table.append(2360964.5) ; table.append(    14.)   ; # 1752 
        table.append(2361695.5) ; table.append(    14.)   ; # 1754 
        table.append(2362425.5) ; table.append(    14.)   ; # 1756 
        table.append(2363156.5) ; table.append(    15.)   ; # 1758 
        table.append(2363886.5) ; table.append(    15.)   ; # 1760 
        table.append(2364617.5) ; table.append(    15.)   ; # 1762 
        table.append(2365347.5) ; table.append(    15.)   ; # 1764 
        table.append(2366078.5) ; table.append(    16.)   ; # 1766 
        table.append(2366808.5) ; table.append(    16.)   ; # 1768 
        table.append(2367539.5) ; table.append(    16.)   ; # 1770 
        table.append(2368269.5) ; table.append(    16.)   ; # 1772 
        table.append(2369000.5) ; table.append(    16.)   ; # 1774 
        table.append(2369730.5) ; table.append(    17.)   ; # 1776 
        table.append(2370461.5) ; table.append(    17.)   ; # 1778 
        table.append(2371191.5) ; table.append(    17.)   ; # 1780 
        table.append(2371922.5) ; table.append(    17.)   ; # 1782 
        table.append(2372652.5) ; table.append(    17.)   ; # 1784 
        table.append(2373383.5) ; table.append(    17.)   ; # 1786 
        table.append(2374113.5) ; table.append(    17.)   ; # 1788 
        table.append(2374844.5) ; table.append(    17.)   ; # 1790 
        table.append(2375574.5) ; table.append(    16.)   ; # 1792 
        table.append(2376305.5) ; table.append(    16.)   ; # 1794 
        table.append(2377035.5) ; table.append(    15.)   ; # 1796 
        table.append(2377766.5) ; table.append(    14.)   ; # 1798 
        table.append(2378496.5) ; table.append(    13.7)  ; # 1800 
        table.append(2379226.5) ; table.append(    13.1)  ; # 1802 
        table.append(2379956.5) ; table.append(    12.7)  ; # 1804 
        table.append(2380687.5) ; table.append(    12.5)  ; # 1806 
        table.append(2381417.5) ; table.append(    12.5)  ; # 1808 
        table.append(2382148.5) ; table.append(    12.5)  ; # 1810 
        table.append(2382878.5) ; table.append(    12.5)  ; # 1812 
        table.append(2383609.5) ; table.append(    12.5)  ; # 1814 
        table.append(2384339.5) ; table.append(    12.5)  ; # 1816 
        table.append(2385070.5) ; table.append(    12.3)  ; # 1818 
        table.append(2385800.5) ; table.append(    12.0)  ; # 1820 
        table.append(2386531.5) ; table.append(    11.4)  ; # 1822 
        table.append(2387261.5) ; table.append(    10.6)  ; # 1824 
        table.append(2387992.5) ; table.append(     9.6)  ; # 1826 
        table.append(2388722.5) ; table.append(     8.6)  ; # 1828 
        table.append(2389453.5) ; table.append(     7.5)  ; # 1830 
        table.append(2390183.5) ; table.append(     6.6)  ; # 1832 
        table.append(2390914.5) ; table.append(     6.0)  ; # 1834 
        table.append(2391644.5) ; table.append(     5.7)  ; # 1836 
        table.append(2392375.5) ; table.append(     5.6)  ; # 1838 
        table.append(2393105.5) ; table.append(     5.7)  ; # 1840 
        table.append(2393836.5) ; table.append(     5.9)  ; # 1842 
        table.append(2394566.5) ; table.append(     6.2)  ; # 1844 
        table.append(2395297.5) ; table.append(     6.5)  ; # 1846 
        table.append(2396027.5) ; table.append(     6.8)  ; # 1848 
        table.append(2396758.5) ; table.append(     7.1)  ; # 1850 
        table.append(2397488.5) ; table.append(     7.3)  ; # 1852 
        table.append(2398219.5) ; table.append(     7.5)  ; # 1854 
        table.append(2398949.5) ; table.append(     7.7)  ; # 1856 
        table.append(2399680.5) ; table.append(     7.8)  ; # 1858 
        table.append(2400410.5) ; table.append(     7.9)  ; # 1860 
        table.append(2401141.5) ; table.append(     7.5)  ; # 1862 
        table.append(2401871.5) ; table.append(     6.4)  ; # 1864 
        table.append(2402602.5) ; table.append(     5.4)  ; # 1866 
        table.append(2403332.5) ; table.append(     2.9)  ; # 1868 
        table.append(2404063.5) ; table.append(     1.6)  ; # 1870 
        table.append(2404793.5) ; table.append(    -1.0)  ; # 1872 
        table.append(2405524.5) ; table.append(    -2.7)  ; # 1874 
        table.append(2406254.5) ; table.append(    -3.6)  ; # 1876 
        table.append(2406985.5) ; table.append(    -4.7)  ; # 1878 
        table.append(2407715.5) ; table.append(    -5.4)  ; # 1880 
        table.append(2408446.5) ; table.append(    -5.2)  ; # 1882 
        table.append(2409176.5) ; table.append(    -5.5)  ; # 1884 
        table.append(2409907.5) ; table.append(    -5.6)  ; # 1886 
        table.append(2410637.5) ; table.append(    -5.8)  ; # 1888 
        table.append(2411368.5) ; table.append(    -5.9)  ; # 1890 
        table.append(2412098.5) ; table.append(    -6.2)  ; # 1892 
        table.append(2412829.5) ; table.append(    -6.4)  ; # 1894 
        table.append(2413559.5) ; table.append(    -6.1)  ; # 1896 
        table.append(2414290.5) ; table.append(    -4.7)  ; # 1898 
        table.append(2415020.5) ; table.append(    -2.7)  ; # 1900 
        table.append(2415750.5) ; table.append(    -0.0)  ; # 1902 
        table.append(2416480.5) ; table.append(    +2.6)  ; # 1904 
        table.append(2417211.5) ; table.append(     5.4)  ; # 1906 
        table.append(2417941.5) ; table.append(     7.7)  ; # 1908 
        table.append(2418672.5) ; table.append(    10.5)  ; # 1910 
        table.append(2419402.5) ; table.append(    13.4)  ; # 1912 
        table.append(2420133.5) ; table.append(    16.0)  ; # 1914 
        table.append(2420863.5) ; table.append(    18.2)  ; # 1916 
        table.append(2421594.5) ; table.append(    20.2)  ; # 1918 
        table.append(2422324.5) ; table.append(    21.2)  ; # 1920 
        table.append(2423055.5) ; table.append(    22.4)  ; # 1922 
        table.append(2423785.5) ; table.append(    23.5)  ; # 1924 
        table.append(2424516.5) ; table.append(    23.9)  ; # 1926 
        table.append(2425246.5) ; table.append(    24.3)  ; # 1928 
        table.append(2425977.5) ; table.append(    24.0)  ; # 1930 
        table.append(2426707.5) ; table.append(    23.9)  ; # 1932 
        table.append(2427438.5) ; table.append(    23.9)  ; # 1934 
        table.append(2428168.5) ; table.append(    23.7)  ; # 1936 
        table.append(2428899.5) ; table.append(    24.0)  ; # 1938 
        table.append(2429629.5) ; table.append(    24.3)  ; # 1940 
        table.append(2430360.5) ; table.append(    25.3)  ; # 1942 
        table.append(2431090.5) ; table.append(    26.2)  ; # 1944 
        table.append(2431821.5) ; table.append(    27.3)  ; # 1946 
        table.append(2432551.5) ; table.append(    28.2)  ; # 1948 
        table.append(2433282.5) ; table.append(    29.1)  ; # 1950 
        table.append(2434012.5) ; table.append(    30.0)  ; # 1952 
        table.append(2434743.5) ; table.append(    30.7)  ; # 1954 
        table.append(2435473.5) ; table.append(    31.4)  ; # 1956 
        table.append(2436204.5) ; table.append(    32.2)  ; # 1958 
        table.append(2436934.5) ; table.append(    33.1)  ; # 1960 
        table.append(2437665.5) ; table.append(    34.0)  ; # 1962 
        table.append(2438395.5) ; table.append(    35.0)  ; # 1964 
        table.append(2439126.5) ; table.append(    36.5)  ; # 1966 
        table.append(2439856.5) ; table.append(    38.3)  ; # 1968 
        table.append(2440587.5) ; table.append(    40.2)  ; # 1970 
        table.append(2441317.5) ; table.append(  42.184)  ; # 1972 JAN  1 
        table.append(2441499.5) ; table.append(  43.184)  ; # 1972 JUL  1 
        table.append(2441683.5) ; table.append(  44.184)  ; # 1973 JAN  1 
        table.append(2442048.5) ; table.append(  45.184)  ; # 1974 JAN  1 
        table.append(2442413.5) ; table.append(  46.184)  ; # 1975 JAN  1 
        table.append(2442778.5) ; table.append(  47.184)  ; # 1976 JAN  1 
        table.append(2443144.5) ; table.append(  48.184)  ; # 1977 JAN  1 
        table.append(2443509.5) ; table.append(  49.184)  ; # 1978 JAN  1 
        table.append(2443874.5) ; table.append(  50.184)  ; # 1979 JAN  1 
        table.append(2444239.5) ; table.append(  51.184)  ; # 1980 JAN  1 
        table.append(2444786.5) ; table.append(  52.184)  ; # 1981 JUL  1 
        table.append(2445151.5) ; table.append(  53.184)  ; # 1982 JUL  1 
        table.append(2445516.5) ; table.append(  54.184)  ; # 1983 JUL  1 
        table.append(2446247.5) ; table.append(  55.184)  ; # 1985 JUL  1 
        table.append(2447161.5) ; table.append(  56.184)  ; # 1988 JAN  1 
        table.append(2447892.5) ; table.append(  57.184)  ; # 1990 JAN  1 
        table.append(2448257.5) ; table.append(  58.184)  ; # 1991 JAN  1 
        table.append(2448804.5) ; table.append(  59.184)  ; # 1992 JUL  1 
        table.append(2449169.5) ; table.append(  60.184)  ; # 1993 JUL  1 
        table.append(2449534.5) ; table.append(  61.184)  ; # 1994 JUL  1 
        table.append(2450083.5) ; table.append(  62.184)  ; # 1996 JAN  1 
        table.append(2450630.5) ; table.append(  63.184)  ; # 1997 JUL  1 
        table.append(2451179.5) ; table.append(  64.184)  ; # 1999 JAN  1 
        table.append(2453736.5) ; table.append(  65.184)  ; # 2006 JAN  1 
        table.append(2454832.5) ; table.append(  66.184)  ; # 2009 JAN  1 
        table.append(2456109.5) ; table.append(  67.184)  ; # 2012 JUL  1 UTC-TAI = - 35s 
        table.append(2457204.5) ; table.append(  68.184)  ; # 2015 JUL  1 UTC-TAI = - 36s
        table.append(2457754.5) ; table.append(  69.184)  ; # 2017 JAN  1 UTC-TAI = - 37s
        jjmax=table[(indexmax-1)*2]
        if (jj<=2067314.5):
            # --- date <= anne948 ---
            t=(jj- 2451545.0)/36525.
            ds=2715.6+573.36*t+46.5*t*t
            dt=ds;
            return dt
        if (jj<=2312752.5):
            #--- date <= anne1620 ---
            t=(jj- 2451545.0)/36525.
            ds=50.6+67.5*t+22.5*t*t
            dt=ds
            return dt
        if (jj<=jjmax):
            # --- date <= indexmax ---
            k2=indexmax
            for k in range(1,indexmax):
                k2=k
                jj2=table[k2*2]
                if (jj<=jj2):
                    break
            jj2=table[k2*2]
            ds2=table[k2*2+1]
            jj1=table[(k2-1)*2]
            ds1=table[(k2-1)*2+1]
            ds=ds1+(jj-jj1)/(jj2-jj1)*(ds2-ds1)
            dt=ds
            return dt
        # --- extrapolation ---
        ds=table[2*indexmax+1]
        dt=ds
        return dt

# ========================================================
# === internal methods : Physical data for planets
# ========================================================

    def _mc_magplanet(self, r, delta, planete_num, phase, l, b):
        """
        /***************************************************************************/
        /* Calcule la magnitude d'une planete                                      */
        /***************************************************************************/
        /* Meuus page 269 et 302-303 pour Saturne                                  */
        /***************************************************************************/
        """
        i=phase/(self._DR)
        mag=0.
        diamapp=0.
        if (r==0.):
            r=1e-10
        if (delta==0.):
            delta=1e-10
        if (r*delta>0):
            mag=5*log10(r*delta)
            if (planete_num == self._PLN_SOLEIL):
                mag=2.5*log10(delta)
        if (planete_num == self._PLN_MERCURE):
                mag+=(-0.42+.0380*i-.000273*i*i+.000002*i*i*i)
                diamapp=2.*atan(1.6303e-5/delta)
        elif (planete_num == self._PLN_VENUS):
            mag+=(-4.40+.0009*i+.000239*i*i-.00000065*i*i*i)
            diamapp=2.*atan(4.0455e-5/delta)
        elif (planete_num == self._PLN_MARS):
            mag+=(-1.52+.016*i)
            diamapp=2.*atan(2.2694e-5/delta);
        elif (planete_num == self._PLN_JUPITER):
            mag+=(-9.40+.005)
            diamapp=2.*atan(4.7741e-4/delta)
        elif (planete_num == self._PLN_URANUS):
            mag+=(-7.19)
            diamapp=2.*atan(1.6979e-4/delta)
        elif (planete_num == self._PLN_NEPTUNE):
            mag+=(-6.87)
            diamapp=2.*atan(1.6243e-4/delta)
        elif (planete_num == self._PLN_PLUTON):
            mag+=(-1.00)
            diamapp=2.*atan(1.5608e-5/delta)
        elif (planete_num == self._PLN_SOLEIL):
            mag+=(-26.86)
            diamapp=2.*atan(4.6541e-3/delta)
        elif ((planete_num == self._PLN_LUNE) or (planete_num == self._PLN_LUNE_ELP)):
            mag+=(0.38+2.97*(i/100.)-0.78*(i/100.)*(i/100.)+.90*(i/100.)*(i/100.)*(i/100.))
            diamapp=2.*atan(1.1617e-5/delta)
        elif (planete_num == self._PLN_SATURNE):
            i=28.08*(self._DR)
            o=169.51*(self._DR)
            n=113.67*(self._DR)
            bb=asin(sin(i)*cos(b)*sin(l-o)-cos(i)*sin(b))
            lp=l-(.01759/r)*(self._DR)
            bp=b-(.000764*cos(l-n)/r)*(self._DR)
            u1=atan2(sin(i)*sin(bp)+cos(i)*cos(bp)*sin(lp-o),cos(bp)*cos(lp-o))
            u2=atan2(sin(i)*sin(b)+cos(i)*cos(b)*sin(l-o),cos(b)*cos(l-o))
            bb=sin(fabs(bb))
            mag+=(-8.68+.044*fabs(u1-u2)/(self._DR)-2.60*bb+1.25*bb*bb)
            diamapp=2.*atan(4.0395e-4/delta)
        return mag, diamapp

    def _mc_libration(self, jj, planete, longmpc, rhocosphip, rhosinphip, equinoxe, astrometric, asd, dec, delta, r, ls, bs):
        """
        /***************************************************************************/
        /* Calcul de la libration apparentes de la Lune a jj donne.                */
        /***************************************************************************/
        /* D'apres Meeus "Astronomical Algorithms" p341-347                        */
        /* lonc : longitude selenographique du centre topocentrique                */
        /* latc : latitude selenographique du centre topocentrique                 */
        /* p    : position de l'angle de l'axe de rotation (axe des poles)         */
        /* lons : longitude selenographique du centre heliocentrique               */
        /* lats : latitude selenographique du centre heliocentrique                */
        /***************************************************************************/
        inputs asd, dec a l'equinoxe de la date
        """
        
        # --- nutation ---*/
        eps = self._mc_obliqmoy(jj)
        dpsi, deps = self._mc_nutation(jj,1)
        eps += deps

        # --- Moon arguments ---*/
        T=(jj-2451545.)/36525.;
        d=(297.8502042+445267.1115168*T-.00016300*T*T+T*T*T/545868-T*T*T*T/113065000)*(self._DR);
        m=(357.5291092+35999.0502909*T-.0001536*T*T+T*T*T*T/24490000)*(self._DR);
        mp=(134.9634114+477198.8676313*T+.0089970*T*T+T*T*T/69699.-T*T*T*T/14712000)*(self._DR);
        f=(93.2720993+483202.0175273*T-.0034029*T*T+T*T*T/3526000+T*T*T*T/863310000)*(self._DR);
        e=1-.002516*T-.0000074*T*T;
        f=(93.2720993+483202.0175273*T-.0034029*T*T+T*T*T/3526000+T*T*T*T/863310000)*(self._DR);
        o=(125.0445550-1934.1361849*T+0.0020762*T*T+T*T*T/467410.-T*T*T*T/60616000.)*(self._DR);
        i=1.54242*(self._DR);
        f=fmod(4*self._PI+f,2*self._PI);
        # --- parameters of physical libration ---*/
        k1=(119.75+131.849*T)*(self._DR);
        k2=(72.56+20.186*T)*(self._DR);
        rho=-0.02752*cos(mp)-0.02245*sin(f)+0.00684*cos(mp-2*f)-0.00293*cos(2*f)-0.00085*cos(2*f-2*d);
        rho=rho-0.00054*cos(mp-2*d)-0.00020*cos(mp+f)-0.00020*cos(mp+2*f)-0.00020*cos(mp-f);
        rho=rho+0.00014*cos(mp+2*f-2*d);
        sigma=-0.02816*sin(mp)+0.02244*cos(f)-0.00682*sin(mp-2*f)-0.00279*sin(2*f)-0.00083*sin(2*f-2*d);
        sigma=sigma+0.00069*sin(mp-2*d)+0.00040*cos(mp+f)-0.00025*sin(2*mp);
        sigma=sigma-0.00023*sin(mp+2*f)+0.00020*cos(mp-f)+0.00019*sin(mp-f)+0.00013*sin(mp+2*f-2*d);
        sigma=sigma-0.00010*cos(mp-3*f);
        tau=0.02520*e*sin(m)+0.00473*sin(2*mp-2*f)-0.00467*sin(mp)+0.00396*sin(k1)+0.00276*sin(2*mp-2*d);
        tau=tau+0.00196*sin(o)-0.00183*cos(mp-f)+0.00115*sin(mp-2*d)-0.00096*sin(mp-d);
        tau=tau+0.00046*sin(2*f-2*d)-0.00039*sin(mp-f)-0.00032*sin(mp-m-d)+0.00027*sin(2*mp-m-2*d);
        tau=tau+0.00023*sin(k2)-0.00014*sin(2*d)+0.00014*cos(2*mp-2*f)-0.00012*sin(mp-2*f);
        tau=tau-0.00012*sin(2*mp)+0.00011*sin(2*mp-2*m-2*d);
        tau*=(self._DR);
        sigma*=(self._DR);
        rho*=(self._DR);
        # =============== earth lon,lat =====================*/
        # --- coord ecl */
        bb = Angle(dec/self._PI)
        ll = Angle(asd/self._PI)
        coords = Coords((1,ll,bb))
        xeq,yeq,zeq = coords.cart()
        xec,yec,zec = self._mc_xyzeq2ec(xeq,yeq,zeq,eps)
        coords = Coords((xec,yec,zec))
        rr,l,b = coords.sphe("rad","rad")
        
        # --- optical topocentric libration ---*/
        w=l-dpsi-o;
        a=atan2(sin(w)*cos(b)*cos(i)-sin(b)*sin(i),cos(w)*cos(b));
        lp=a-f;
        bp=asin(-sin(w)*cos(b)*sin(i)-sin(b)*cos(i));
        # --- physical topocentric libration ---*/
        lpp=-tau+(rho*cos(a)+sigma*sin(a))*tan(bp);
        bpp=sigma*cos(a)-rho*sin(a);
        # --- topocentric libration  ---*/
        v=fmod(4*self._PI+lp+lpp,2*self._PI);
        if (v>(self._PI)):
            v-=(2*(self._PI))
        lonc=v;
        latc=(bp+bpp);
        # --- pole axis position ---*/
        v=o+deps+sigma/sin(i);
        x=sin(i+rho)*sin(v);
        y=sin(i+rho)*cos(v)*cos(eps)-cos(i+rho)*sin(eps);
        w=atan2(x,y);
        v=asin(sqrt(x*x+y*y)*cos(asd-w)/cos(bp+bpp));
        p=fmod(4*self._PI+v,2*self._PI);
        
        # =============== earth lon,lat =====================*/
        # --- sun ---*/
        l=ls+(self._PI)+delta/r*cos(b)*sin(ls-l);
        b=delta/r*b;
                
        # --- optical topocentric libration ---*/
        w=l-dpsi-o;
        a=atan2(sin(w)*cos(b)*cos(i)-sin(b)*sin(i),cos(w)*cos(b));
        lp=a-f;
        bp=asin(-sin(w)*cos(b)*sin(i)-sin(b)*cos(i));
        # --- physical topocentric libration ---*/
        lpp=-tau+(rho*cos(a)+sigma*sin(a))*tan(bp);
        bpp=sigma*cos(a)-rho*sin(a);
        # --- topocentric libration  ---*/
        v=fmod(4*self._PI+lp+lpp,2*self._PI);
        if (v>(self._PI)):
            v-=(2*(self._PI))
        lons=v;
        lats=(bp+bpp)
        return lonc, latc, p, lons, lats

    def _mc_physephem_moon(self, jj,planete,xg,yg,zg,longmpc,rhocosphip,rhosinphip, equinoxe, astrometric, asd, dec, delta, r, ls, bs):
        """
        /***************************************************************************/
        /* Cacul des parametres d'observation physique des planetes                */
        /***************************************************************************/
        /* D'apres "Astronomy with computer".                                      */
        /* pages 135-154 (chapitre 7)                                              */
        /***************************************************************************/
        """
        if ((planete == self._PLN_LUNE) or (planete == self._PLN_LUNE_ELP)):
            diamapp_equ=0.;
            diamapp_pol=0.;
            long1=0.;
            long2=0.;
            long3=0.;
            lati=0.;
            posangle_sun=0.;
            posangle_north=0.;
            # ---
            long1,lati,posangle_north,lons,lats = self._mc_libration(jj,planete,longmpc,rhocosphip,rhosinphip,equinoxe, astrometric, asd, dec)
            long2=long1;
            long3=long1;
            r=sqrt(xg*xg+yg*yg+zg*zg);
            diamapp_equ=2.*atan(1.1617e-5/r);
            diamapp_pol=diamapp_equ;
            posangle_sun=0.;
            lati_sun=lats;
            long1_sun=lons;
        return diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun

    def _mc_physephem_nomoon(self, jj,planete,xg,yg,zg,x,y,z):
        """
        /***************************************************************************/
        /* Cacul des parametres d'observation physique des planetes                */
        /***************************************************************************/
        /* D'apres "Astronomy with computer".                                      */
        /* pages 135-154 (chapitre 7)                                              */
        /***************************************************************************/
        """
        if ((planete == self._PLN_LUNE) or (planete == self._PLN_LUNE_ELP)):
            return
        a0=0.;d0=0.;w1=0.;w2=0.;w3=0.;f=0.;req=0.;
        sense=1.;
        diamapp_equ=0.;
        diamapp_pol=0.;
        long1=0.;
        long2=0.;
        long3=0.;
        lati=0.;
        posangle_sun=0.;
        posangle_north=0.;
        t=(jj-2451545.)/36525.;
        d=(jj-2451545.);
        if (planete == self._PLN_MERCURE):
            a0=281.01-0.033*t;
            d0=61.45-0.005*t;
            req=2439.;
            f=0.;
            w1=w2=w3=(329.68+6.1385025*d);
            sense=1.;
        elif (planete == self._PLN_VENUS):
            a0=272.76;
            d0=67.16;
            req=6051.;
            f=0.;
            w1=w2=w3=(160.20-1.4813688*d);
            sense=-1.;
        elif (planete == self._PLN_TERRE):
            a0=0.-0.641*t;
            d0=90.0-0.557*t;
            req=6378.14;
            f=0.00335281;
            w1=w2=w3=(190.16+360.9856235*d);
            sense=1.;
        elif (planete == self._PLN_MARS):
            a0=317.681-0.108*t;
            d0=52.886-0.061*t;
            req=3393.4;
            f=0.0051865;
            w1=w2=w3=(176.901+350.8919830*d);
            sense=1.;
        elif (planete == self._PLN_JUPITER):
            a0=268.05-0.009*t;
            d0=64.49+0.003*t;
            req=71398.;
            f=0.0648088;
            w1=(67.1+877.900*d);
            w2=(43.3+870.270*d);
            w3=(284.695+870.536*d);
            sense=1.;
        elif (planete == self._PLN_SATURNE):
            a0=40.589-0.036*t;
            d0=83.537-0.004*t;
            req=60000.;
            f=0.1076209;
            w1=w2=(227.2037+844.300*d);
            w3=(38.90+810.7939024*d);
            sense=1.;
        elif (planete == self._PLN_URANUS):
            a0=257.311;
            d0=-15.175;
            req=25400.;
            f=0.030;
            w1=w2=w3=(203.81-501.1600928*d);
            sense=-1.;
        elif (planete == self._PLN_NEPTUNE):
            n=(357.85+52.316*t)*(self._DR);
            a0=299.36+0.7*sin(n);
            d0=43.46-0.51*cos(n);
            req=24300.;
            f=0.0259;
            w1=w2=w3=(253.18+536.3128492*d-0.48*sin(n));
            sense=1.;
        elif (planete == self._PLN_PLUTON):
            a0=313.02;
            d0=9.09;
            req=1500.;
            f=0.;
            w1=w2=w3=(236.77-56.3623195*d);
            sense=-1.;
        elif (planete == self._PLN_SOLEIL):
            a0=286.13;
            d0=63.87;
            req=696000.;
            f=0.;
            w1=w2=w3=(84.182+14.1844000*d);
            sense=-1.;
        a0*=(self._DR);
        d0*=(self._DR);
        w1*=(self._DR);
        w2*=(self._DR);
        w3*=(self._DR);
        req/=(self._UA*1e-3);
        sina=sin(a0);
        cosa=cos(a0);
        sind=sin(d0);
        cosd=cos(d0);
        dx=cosa*cosd;
        dy=sina*cosd;
        dz=sind;
        rho=sqrt(xg*xg+yg*yg);
        r=sqrt(xg*xg+yg*yg+zg*zg);
        diamapp_equ=2*asin(req/r);
        costh=((-xg*zg)*dx+(-yg*zg)*dy+(xg*xg+yg*yg)*dz)/(r*rho);
        sinth=(-yg*dx+xg*dy)/rho;
        th=atan2(sinth,costh);
        posangle_north=fmod(4*self._PI+th,2*self._PI);
        # /* w1 */
        cosw=cos(w1);
        sinw=sin(w1);
        e1x=-cosw*sina-sinw*sind*cosa;
        e1y= cosw*cosa-sinw*sind*sina;
        e1z= sinw*cosd;
        e2x= sinw*sina-cosw*sind*cosa;
        e2y=-sinw*cosa-cosw*sind*sina;
        e2z= cosw*cosd;
        e3x= cosd*cosa;
        e3y= cosd*sina;
        e3z= sind;
        sx=-(e1x*xg+e1y*yg+e1z*zg);
        sy=-(e2x*xg+e2y*yg+e2z*zg);
        sz=-(e3x*xg+e3y*yg+e3z*zg);
        phip=atan2(sz,sqrt(sx*sx+sy*sy));
        phi=atan2(tan(phip),((1-f)*(1-f)));
        lambdp=atan2(sy,sx);
        lambd=-1.*sense*lambdp;
        lati=phi;
        diamapp_pol=(diamapp_equ)*(1-f*cos(phip)*cos(phip));
        long1=fmod(4*self._PI+lambd,2*self._PI);
        # /* w1 sun */
        sx=-(e1x*x+e1y*y+e1z*z);
        sy=-(e2x*x+e2y*y+e2z*z);
        sz=-(e3x*x+e3y*y+e3z*z);
        phip=atan2(sz,sqrt(sx*sx+sy*sy));
        phi=atan2(tan(phip),((1-f)*(1-f)));
        lambdp=atan2(sy,sx);
        lambd=-1.*sense*lambdp;
        lati_sun=phi;
        long1_sun=fmod(4*self._PI+lambd,2*self._PI);
        # /* w2 */
        cosw=cos(w2);
        sinw=sin(w2);
        e1x=-cosw*sina-sinw*sind*cosa;
        e1y= cosw*cosa-sinw*sind*sina;
        e1z= sinw*cosd;
        e2x= sinw*sina-cosw*sind*cosa;
        e2y=-sinw*cosa-cosw*sind*sina;
        e2z= cosw*cosd;
        e3x= cosd*cosa;
        e3y= cosd*sina;
        e3z= sind;
        sx=-(e1x*xg+e1y*yg+e1z*zg);
        sy=-(e2x*xg+e2y*yg+e2z*zg);
        lambdp=atan2(sy,sx);
        lambd=-1.*sense*lambdp;
        long2=fmod(4*self._PI+lambd,2*self._PI);
        # /* w3 */
        cosw=cos(w3);
        sinw=sin(w3);
        e1x=-cosw*sina-sinw*sind*cosa;
        e1y= cosw*cosa-sinw*sind*sina;
        e1z= sinw*cosd;
        e2x= sinw*sina-cosw*sind*cosa;
        e2y=-sinw*cosa-cosw*sind*sina;
        e2z= cosw*cosd;
        e3x= cosd*cosa;
        e3y= cosd*sina;
        e3z= sind;
        sx=-(e1x*xg+e1y*yg+e1z*zg);
        sy=-(e2x*xg+e2y*yg+e2z*zg);
        sz=-(e3x*xg+e3y*yg+e3z*zg);
        lambdp=atan2(sy,sx);
        lambd=-1.*sense*lambdp;
        long3=fmod(4*self._PI+lambd,2*self._PI);
        # /* sun */
        d=sqrt(x*x+y*y+z*z);
        dx=-x;
        dy=-y;
        dz=-z;
        costh=((-xg*zg)*dx+(-yg*zg)*dy+(xg*xg+yg*yg)*dz)/(r*rho);
        sinth=(-yg*dx+xg*dy)/rho;
        th=atan2(sinth,costh);
        posangle_sun=fmod(4*self._PI+th,2*self._PI);
        if (planete == self._PLN_SOLEIL):
            posangle_sun=0.;
            lati_sun=lati;
            long1_sun=long1;
        return diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun

# ========================================================
# === internal methods : Misc
# ========================================================

    def _mc_altitude2tp(self, altitude_m, pressure_0m_pascal=101325):
        """
        /****************************************************************************/
        /* Copute the ICAO temperature & pressure from a given altitude             */
        /****************************************************************************/
        /****************************************************************************/
        """
        alti=altitude_m
        p0m=pressure_0m_pascal
        tk0m=273.15+15;
        if (alti<11000):
            tk=tk0m-0.0065*alti
            p=p0m*pow(tk/tk0m,5.255)
        elif (alti<15000):
            tk0m=273.15+15
            tk=tk0m-0.0065*11000
            p=p0m*pow((tk0m-0.0065*alti)/tk0m,5.255)
        elif ((alti>=15000) and (alti<20000)):
            h1=15000; p1=p0m*pow((tk0m-0.0065*h1)/tk0m,5.255); t1=tk0m-0.0065*11000
            h2=20000; p2=5500; t2=273.15-46
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=20000) and (alti<30000)):
            h1=20000; p1=5500; t1=273.15-46
            h2=30000; p2=1100; t2=273.15-38
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=30000) and (alti<40000)):
            h1=30000; p1=1100; t1=273.15-38
            h2=40000; p2=300; t2=273.15-5
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=40000) and (alti<50000)):
            h1=40000; p1=300; t1=273.15-5
            h2=50000; p2=90; t2=273.15+1
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=50000) and (alti<60000)):
            h1=50000; p1=90; t1=273.15+1
            h2=60000; p2=25; t2=273.15-20
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=60000) and (alti<100000)):
            h1=60000; p1=25; t1=273.15-20
            h2=100000; p2=0.04; t2=273.15-64
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=100000) and (alti<200000)):
            h1=100000; p1=0.04; t1=273.15-64
            h2=200000; p2=1.3e-4; t2=273.15-82.2
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=200000) and (alti<400000)):
            h1=200000; p1=1.3e-4; t1=273.15-82.2
            h2=400000; p2=4.4e-6; t2=273.15-97.3
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif ((alti>=200000) and (alti<400000)):
            h1=400000; p1=4.4e-6; t1=273.15-97.3
            h2=500000; p2=0; t2=273.15-97.7
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        else:
            p=0
            tk=273.15-97.7
        return tk,p

# ========================================================
# === internal methods : Ephemeris
# ========================================================

    def _mc_radec2app(self, date_utc, longmpc, rhocosphip, rhosinphip, ra, dec, equinox, epoch, mura, mudec, plx, tk,ppa,hump,lnm,outputs):
        # --- date        
        jjutc = date_utc.jd()
        jds = [ jjutc ]

        # --- home
        longmpc *= self._DR
        latitude, altitude = self._mc_rhophi2latalt(rhosinphip, rhocosphip)
        latrad = latitude * self._DR

        # --- outputs (list of keywords, such as RA, DEC, ELONG)
        outputs = [ output.upper() for output in outputs ]
        results = []
        
        #print("-1. ra={} dec={}".format(ra,dec))
        ra*=self._DR ;
        dec*=self._DR;
        #print("0. ra={} dec={}".format(ra/self._DR,dec/self._DR))
        muramas=mura;
        mudecmas=mudec;
        for jd in jds:
            
            #print("1. ra={} dec={}".format(ra/self._DR,dec/self._DR))
            cosdec=cos(dec);
            mura=muramas*1e-3/86400/cosdec;
            mudec=mudecmas*1e-3/86400;
            parallax=plx;
                
            # --- aberration annuelle ---*/
            asd2, dec2 = self._mc_aberration_annuelle(jd,ra,dec,1)
            ra=asd2
            dec=dec2
            #print("2. ra={} dec={}".format(ra/self._DR,dec/self._DR))
            
            #--- calcul de mouvement propre ---*/
            ra+=(jd-epoch)/365.25*mura;
            dec+=(jd-epoch)/365.25*mudec;
            #print("3. ra={} dec={}".format(ra/self._DR,dec/self._DR))
            
				# --- calcul de la precession ---*/
            asd2, dec2 = self._mc_precad(equinox,ra,dec,jd)
            ra=asd2;
            dec=dec2;
            #print("4. ra={} dec={}".format(ra/self._DR,dec/self._DR))
            
            # --- correction de parallaxe stellaire*/
            if (parallax>0):
                asd2, dec2 = self._mc_parallaxe_stellaire(jd,ra,dec,parallax)
                ra=asd2;
                dec=dec2;
            #print("5. ra={} dec={}".format(ra/self._DR,dec/self._DR))

            # --- correction de nutation */
            asd2, dec2 = self._mc_nutradec(jd,ra,dec,1);
            ra=asd2
            dec=dec2
            #print("6. ra={} dec={}".format(ra/self._DR,dec/self._DR))

            # --- aberration de l'aberration diurne*/
            asd2, dec2 = self._mc_aberration_diurne(jd,ra,dec,longmpc,rhocosphip,rhosinphip,1)
            ra=asd2;
            dec=dec2;
            #print("7. ra={} dec={}".format(ra/self._DR,dec/self._DR))

				# --- coordonnees horizontales---*/
            ha = self._mc_ad2hd(jd,longmpc,ra)
            az, h = self._mc_hd2ah(ha,dec,latrad)

				# --- refraction ---*/
            #lnm=590,
            #refraction = self._mc_refraction2(h,1,tk,ppa,lnm,hump,latitude,altitude)
            #h += refraction;     
            #ha, dec = self._mc_ah2hd(az,h,latrad)

            parallactic = self._mc_hd2parallactic(ha,dec,latrad)
            ra = self._mc_hd2ad(jd,longmpc,ha)

        if 'RA_APP' in outputs or 'RA' in outputs:
            results.append(['RA_APP',ra/self._DR])
        if 'DEC_APP' in outputs or 'DEC' in outputs :
            results.append(['DEC_APP',dec/self._DR])
        if 'HA' in outputs:
            results.append(['HA',ha/self._DR])
        if 'AZ' in outputs:
            results.append(['AZ',az/self._DR])
        if 'ELEV' in outputs or 'ALT' in outputs or 'ALTITUDE' in outputs:
            results.append(['ELEV',h/self._DR])
        if 'PARALLACTIC' in outputs:
            results.append(['PARALLACTIC',parallactic/self._DR])
        return dict(results)
    
    def _mc_ephem_planet(self, planet_name, date_utc, longmpc, rhocosphip, rhosinphip, astrometric, jd_equinox, outputs):
        """ Compute a planet ephemeris.
        
        :param planet_name: The name of a planet.
        :type planet_name: string
        :param date_utc: The Date (UTC)
        :type date_utc: Date
        :param coords_home: The Coords of the home place (for parallax)
        :type coords_home: Coords
        :param astrometric: Flag to compute astrometric (=1) or true coordinates (=0)
        :type astrometric: int
        :param jd_equinox: The equinox julian day if astrometric flag = 1
        :type jd_equinox: float
        """
        # --- planet
        planet_num = self._mc_name2planetnum(planet_name)
        
        # --- The case of a false planet
        if (planet_num == self._PLN_OTHERPLANET):
            raise Exception
            return ""
        
        # --- date        
        jjutc = date_utc.jd()
        jj = self._mc_tu2td(jjutc)
        jjd=jj

        # --- home
        #altitude, longmpc , latitude = coords_home.sphe("rad","deg")
        latitude, altitude = self._mc_rhophi2latalt(rhosinphip, rhocosphip)
        longmpc*=self._DR
        
        # --- computation method
        afm = self._AFM_VFP79
        
        # --- outputs (list of keywords, such as RA, DEC, ELONG)
        outputs = [ output.upper() for output in outputs ]
        results = []
        
        if (planet_num == self._PLN_LUNE) or (planet_num == self._PLN_LUNE_ELP):
            # =====================================
            # === The moon only
            # =====================================
            
            ttminusutc=jj-jjutc
            # --- soleil ---
            delta = 1. ; # correction de l'abberation planetaire a priori pour la Lune            
            jjds = self._mc_aberpla(jjd,delta)
            llp, mmp, uup = self._mc_jd2lbr1a(jjds);
            coords = self._mc_jd2lbr1b(jjds, self._PLN_SOLEIL, llp, mmp, uup, afm)
            xs,ys,zs = coords.cart()
            rs,ls,bs = coords.sphe("rad","rad")            
            
            # --- correction de la parallaxe ---
            eps = self._mc_obliqmoy(jjds)
            xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
            dxeq,dyeq,dzeq = self._mc_paraldxyzeq(jjds-ttminusutc,longmpc,rhocosphip,rhosinphip)
            xs-=dxeq
            ys-=dyeq
            zs-=dzeq
            xs,ys,zs = self._mc_xyzeq2ec(xs,ys,zs,eps)

            # --- soleil : coordonnes asd,dec du Soleil ---
            if (astrometric==0):
                coords = Coords((xs,ys,zs))
                rs,ls,bs = coords.sphe("rad","rad")
                dpsi, deps = self._mc_nutation(jjd,1)
                ls += dpsi
                eps += deps
                bb = Angle(bs/self._DR)
                ll = Angle(ls/self._DR)
                coords = Coords((rs,ll,bb))
                xs,ys,zs = coords.cart()
            xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
            coords = Coords((xs,ys,zs))
            delta,asds,decs = coords.sphe("rad","rad")

            # a ce niveau on retient (asds,decs) pour le calcul de la phase */
            # On recommence tout le calcul sans abberation pour la position */
            #dans le repere heliocentrique */

            # --- Terre ---
            llp, mmp, uup = self._mc_jd2lbr1a(jjd);
            coords = self._mc_jd2lbr1b(jjd, self._PLN_SOLEIL, llp, mmp, uup, afm)
            xs,ys,zs = coords.cart()
            r,ls,bs = coords.sphe("rad","rad")            
            
            # --- Terre : correction de la parallaxe ---
            eps = self._mc_obliqmoy(jjd)
            xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
            dxeq,dyeq,dzeq = self._mc_paraldxyzeq(jjd-ttminusutc,longmpc,rhocosphip,rhosinphip)
            xs-=dxeq
            ys-=dyeq
            zs-=dzeq
            xs,ys,zs = self._mc_xyzeq2ec(xs,ys,zs,eps)
            
            # --- Terre : coordonnes asd,dec du Soleil ---
            if (astrometric==0):
                coords = Coords((xs,ys,zs))
                rs,ls,bs = coords.sphe("rad","rad")
                dpsi, deps = self._mc_nutation(jjd,1)
                ls += dpsi
                eps += deps
                bb = Angle(bs/self._DR)
                ll = Angle(ls/self._DR)
                coords = Coords((r,ll,bb))
                xs,ys,zs = coords.cart()
            xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
            coords = Coords((xs,ys,zs))
            delta,asd,dec = coords.sphe("rad","rad")
            
            # ici on ne garde que (xs,ys,zs,delta) pour la position geocentrique

            # --- LUNE : planete (coordonnees ecliptiques astrometriques a la date) ---
            llp, mmp, uup = self._mc_jd2lbr1a(jjd);
            coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
            x,y,z = coords.cart()
            eps = self._mc_obliqmoy(jjd)
            x,y,z = self._mc_xyzec2eq(x,y,z,eps)
            coords = Coords((x,y,z))
            delta,asd,dec = coords.sphe("rad","rad")

            # --- planete corrigee de l'aberration de la lumiere ---
            jjd = self._mc_aberpla(jjd,delta)
            llp, mmp, uup = self._mc_jd2lbr1a(jjd);
            coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
            x,y,z = coords.cart()
            eps = self._mc_obliqmoy(jjd)
            
            # --- correction de la nutation ---
            r,l,b = coords.sphe()
            if (astrometric==0):
                dpsi, deps = self._mc_nutation(jjd,1)
                l += dpsi
                eps += deps
                bb = Angle(b/self._DR)
                ll = Angle(l/self._DR)
                coords = Coords((r,ll,bb))
                x,y,z = coords.cart()
                
            # --- correction de la parallaxe ---
            x,y,z = self._mc_xyzec2eq(x,y,z,eps)
            x-=dxeq
            y-=dyeq
            z-=dzeq
            
            # --- coord. spheriques ---
            coords = Coords((x,y,z))
            delta,asd,dec = coords.sphe("rad","rad")

            #--- parametres physiques --- (not implemented)
            if ('APPDIAMEQU' in outputs) or ('APPDIAMPOL' in outputs) or ('LONGI' in outputs) or ('LONGII' in outputs) or ('LONGIII' in outputs) or ('LATI' in outputs) or ('POSNORTH' in outputs) or ('POSSUN' in outputs) or ('LONGI_SUN' in outputs) or ('LATI_SUN' in outputs):
                diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun = self._mc_physephem_moon(jjd,planet_num,x,y,z,longmpc, rhocosphip, rhosinphip, asd, dec, delta, r, ls, bs)
            
            # --- parametres elong et magnitude ---
            x,y,z = self._mc_he2ge(x,y,z,-xs,-ys,-zs)
            r=sqrt(x*x+y*y+z*z)
            if ('ELONG' in outputs) or ('PHASE' in outputs) or ('POSSUN' in outputs) or ('MAG' in outputs) or ('APPDIAM' in outputs):
                elong, phase, limb = self._mc_elonphaslimb(asd,dec,asds,decs,rs,delta);
                posangle_sun = limb
            if ('MAG' in outputs) or ('APPDIAM' in outputs):
                mag, diamapp = self._mc_magplanet(r, delta, planet_num,phase, l, b)

            # --- correction de l'aberration annuelle ---
            if (astrometric==0):
                asd, dec = self._mc_aberration_annuelle(jjd, asd, dec, 1)
    
            # --- correction de la precession ---
            asd,dec = self._mc_precad(jjd,asd,dec,jd_equinox); # equatorial J2000

            # --- parametres physiques ---*/
            
        else:
            # =====================================
            # === All planets and Sun
            # =====================================
            
            # --- soleil ---
            llp, mmp, uup = self._mc_jd2lbr1a(jjd);
            coords = self._mc_jd2lbr1b(jjd, self._PLN_SOLEIL, llp, mmp, uup, afm)
            xs,ys,zs = coords.cart()
            rs,ls,bs = coords.sphe("rad","rad")
            if (planet_num == self._PLN_SOLEIL):
                # --- special case of Sun
                jjd = self._mc_aberpla(jjd,rs)
                llp, mmp, uup = self._mc_jd2lbr1a(jjd);
                coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
                xs,ys,zs = coords.cart()
                rs,ls,bs = coords.sphe("rad","rad")            
            
            # --- correction de la parallaxe ---
            eps = self._mc_obliqmoy(jjd)
            xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
            dxeq,dyeq,dzeq = self._mc_paraldxyzeq(jjd,longmpc,rhocosphip,rhosinphip)
            xs-=dxeq
            ys-=dyeq
            zs-=dzeq
            xs,ys,zs = self._mc_xyzeq2ec(xs,ys,zs,eps)
    
            if (planet_num == self._PLN_SOLEIL):
                xg=xs; yg=ys; zg=zs
            else:
                # --- planet ---
                coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
                x,y,z = coords.cart()   
                x,y,z = self._mc_he2ge(x,y,z,xs,ys,zs)
                x,y,z = self._mc_xyzec2eq(x,y,z,eps)
                coords = Coords((x,y,z))
                delta,asd,dec = coords.sphe("rad","rad")
        
                # --- planete corrigee de l'aberration de la lumiere ---
                jjd = self._mc_aberpla(jjd,delta)
                llp, mmp, uup = self._mc_jd2lbr1a(jjd)
                coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
                x,y,z = coords.cart()
                rr,l,b = coords.sphe("rad","rad")
                xg,yg,zg = self._mc_he2ge(x,y,z,xs,ys,zs)
                eps = self._mc_obliqmoy(jjd)
    
            # --- correction de la nutation ---
            if (astrometric==0):
                coords = Coords((xg,yg,zg))
                r,l,b = coords.sphe("rad","rad")
                dpsi, deps = self._mc_nutation(jjd,1)
                l += dpsi
                eps += deps
                bb = Angle(b/self._DR)
                ll = Angle(l/self._DR)
                coords = Coords((r,ll,bb))
                xg,yg,zg = coords.cart()
    
            # --- coord. spheriques ---
            xg,yg,zg = self._mc_xyzec2eq(xg,yg,zg,eps)
            coords = Coords((xg,yg,zg))
            delta,asd,dec = coords.sphe("rad","rad")
    
            # --- correction de l'aberration annuelle ---
            if (astrometric==0):
                asd, dec = self._mc_aberration_annuelle(jjd, asd, dec, 1)
    
            # --- correction de la precession ---
            if (planet_num != self._PLN_PLUTON):
                asd,dec = self._mc_precad(jjd,asd,dec,jd_equinox); # equatorial J2000
    
            # --- parametres elong et magnitude ---*/
            if (planet_num != self._PLN_SOLEIL):
                r = rr
                if ('ELONG' in outputs) or ('PHASE' in outputs) or ('MAG' in outputs) or ('DIAMAPP' in outputs):
                    elong, phase = self._mc_elonphas(r, rs, delta)
            else:
                elong = 0
                phase = 0
            if ('MAG' in outputs) or ('APPDIAM' in outputs):
                mag, diamapp = self._mc_magplanet(r, delta, planet_num, phase, l, b)

            # --- parametres physiques ---*/
            if ('APPDIAMEQU' in outputs) or ('APPDIAMPOL' in outputs) or ('LONGI' in outputs) or ('LONGII' in outputs) or ('LONGIII' in outputs) or ('LATI' in outputs) or ('POSNORTH' in outputs) or ('POSSUN' in outputs) or ('LONGI_SUN' in outputs) or ('LATI_SUN' in outputs):
                x,y,z = self._mc_xyzec2eq(x,y,z,eps) ; # coord. helio equatoriales
                diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun = self._mc_physephem_nommon(jjd,planet_num,xg,yg,zg,x,y,z)
        
        # --- outputs 
        if 'NAME' in outputs:
            results.append(['NAME',planet_name])
        if 'RA' in outputs:
            results.append(['RA',asd/self._DR])
        if 'DEC' in outputs:
            results.append(['DEC',dec/self._DR])
        if 'EQUINOX' in outputs:
            results.append(['EQUINOX',(Date(jd_equinox)).equinox()])
        if 'ELONG' in outputs:
            results.append(['ELONG',elong/self._DR])
        if 'PHASE' in outputs:
            results.append(['PHASE',phase/self._DR])
        if 'DELTA' in outputs:
            results.append(['DELTA',delta])
        if 'R' in outputs:
            results.append(['R',rs])
        if 'MAG' in outputs:
            results.append(['MAG',mag])
        if 'APPDIAM' in outputs:
            results.append(['APPDIAM',diamapp/self._DR])
        if 'APPDIAMEQU' in outputs:
            results.append(['APPDIAMEQU',diamapp_equ/self._DR])
        if 'APPDIAMPOL' in outputs:
            results.append(['APPDIAMPOL',diamapp_pol/self._DR])
        if 'LONGI' in outputs:
            results.append(['LONGI',long1/self._DR])
        if 'LONGII' in outputs:
            results.append(['LONGII',long2/self._DR])
        if 'LONGIII' in outputs:
            results.append(['LONGIII',long3/self._DR])
        if 'LATI' in outputs:
            results.append(['LATI',lati/self._DR])
        if 'POSNORTH' in outputs:
            results.append(['POSNORTH',posangle_north/self._DR])
        if 'POSSUN' in outputs:
            results.append(['POSNORTH',posangle_sun/self._DR])
        if 'LONGI_SUN' in outputs:
            results.append(['LONGI_SUN',long1_sun/self._DR])
        if 'LATI_SUN' in outputs:
            results.append(['LATI_SUN',lati_sun/self._DR])
        return dict(results)
        
# ========================================================
# === special methods
# ========================================================
        
    def __init__(self):
        """ Object initialization where planet_name is the name of the planet
        """
        self._init_mechanics()