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.app