from math import pi, sin, cos, fmod, tan, atan, fabs, atan2, asin, acos, sqrt, log10 import doctest from .dates import Date from .angles import Angle from .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 | */ /* ||<- 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()