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
# ========================================================
# ========================================================
[docs]class Mechanics:
""" Class to compute celestial mechanics
"""
# ========================================================
# === attributs
# ========================================================
_DR = pi/180
_PI = pi
_PISUR2 = pi/2
# analytical fit methods
_AFM_VFP79 = 0
_AFM_VSOP87 = 1
_AFM_GMS86 = 2
_AFM_ELP2000_82B = 3
_AFM_MELP98 = 4
# planet names
_PLN_OTHERPLANET = -2
_PLN_ALLPLANETS = -1
_PLN_SOLEIL = 0
_PLN_MERCURE = 1
_PLN_VENUS = 2
_PLN_TERRE = 3
_PLN_MARS = 4
_PLN_JUPITER = 5
_PLN_SATURNE = 6
_PLN_URANUS = 7
_PLN_NEPTUNE = 8
_PLN_PLUTON = 9
_PLN_LUNE_ELP = 10
_PLN_LUNE = 11
_PLN_NB_PLANETES = 12
# astronomical constants
_CST_UA = 1.49597870691e11 ; # (m)
_CST_CLIGHT = 2.99792458e8 ; # (m/s)
_CST_EARTH_SEMI_MAJOR_RADIUS = 6378137 ; # (m) WGS 84
_CST_EARTH_INVERSE_FLATTENING = 298.257223563 ; # WGS 84
# ========================================================
# === internal methods : Generals
# ========================================================
def _init_mechanics(self):
""" Object initialization where planet_name is the name of the planet
:param planet_name: A string (cf. help(PLanet))
:type planet_name: string
"""
pass
# ========================================================
# === internal methods : Planet orbits
# ========================================================
def _mc_name2planetnum(self, planet_name):
name = str(planet_name[0:3]).upper()
planet_num = self._PLN_OTHERPLANET
if (name == "SUN"):
planet_num = self._PLN_SOLEIL
elif (name == "MER"):
planet_num = self._PLN_MERCURE
elif (name == "VEN"):
planet_num = self._PLN_VENUS
elif (name == "EAR"):
planet_num = self._PLN_TERRE
elif (name == "MAR"):
planet_num = self._PLN_MARS
elif (name == "JUP"):
planet_num = self._PLN_JUPITER
elif (name == "SAT"):
planet_num = self._PLN_SATURNE
elif (name == "URA"):
planet_num = self._PLN_URANUS
elif (name == "NEP"):
planet_num = self._PLN_NEPTUNE
elif (name == "PLU"):
planet_num = self._PLN_PLUTON
elif (name == "ELP"):
planet_num = self._PLN_LUNE_ELP
elif (name == "MOO"):
planet_num = self._PLN_LUNE
return planet_num
def _mc_jd2lbr1a(self, jj):
"""
/***************************************************************************/
/* Retourne les valeurs des tableaux L, M, U necessaires pour tenir */
/* compte des principales perturbations planetaires a entrer dans */
/* la fonction mc_jd2lbr1b */
/***************************************************************************/
/* Les tableaux *l, *m et *u doivent etre dimensiones chacun avec */
/* 9 elements (0 a 8) */
/***************************************************************************/
"""
# l, m, u
t=(jj-2415020.0)/36525;
l=[0]*self._PLN_NB_PLANETES
m=[0]*self._PLN_NB_PLANETES
u=[0]*self._PLN_NB_PLANETES
l[0]=(279.6964027+36000.7695173*t)*self._DR;
m[0]=(358.4758635+35999.0494965*t)*self._DR;
u[0]=(270.435377+481267.880863*t)*self._DR;
l[1]=(178.178814+149474.071386*t)*self._DR;
m[1]=(102.279426+149472.515334*t)*self._DR;
u[1]=(131.032888+149472.885872*t)*self._DR;
l[2]=(342.766738+58519.212542*t)*self._DR;
m[2]=(212.601892+58517.806388*t)*self._DR;
u[2]=(266.987445+58518.311835*t)*self._DR;
l[3]=(293.747201+19141.699879*t)*self._DR;
m[3]=(319.529273+19139.858887*t)*self._DR;
u[3]=(244.960887+19140.928953*t)*self._DR;
l[4]=(237.352259+3034.906621*t)*self._DR;
m[4]=(225.444539+3034.906621*t)*self._DR;
u[4]=(138.419219+3034.906621*t)*self._DR;
l[5]=(265.869357+1222.116843*t)*self._DR;
m[5]=(175.758477+1222.116843*t)*self._DR;
u[5]=(153.521637+1222.116843*t)*self._DR;
l[6]=(243.362437+429.898403*t)*self._DR;
m[6]=( 74.313637+429.898403*t)*self._DR;
u[6]=(169.872293+429.388747*t)*self._DR;
l[7]=( 85.024943+219.863377*t)*self._DR;
m[7]=( 41.269103+219.863377*t)*self._DR;
u[7]=(314.346275+218.761885*t)*self._DR;
l[8]=( 92.312712+146.674728*t)*self._DR;
m[8]=(229.488633+145.278567*t)*self._DR;
u[8]=(343.369233+145.278567*t)*self._DR;
res = (l, m, u)
return res
def _mc_jd2lbr1b(self, jj, planete, l, m, u, afm):
"""
/***************************************************************************/
/* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr */
/* pour l'equinoxe moyen de la date. () */
/***************************************************************************/
/* Les tableaux l, m et u doivent etre dimensiones chacun avec */
/* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a */
/***************************************************************************/
"""
coords = Coords((0,0,0))
if afm == self._AFM_VFP79:
t=(jj-2415020.0)/36525
if ((planete==self._PLN_SOLEIL) or (planete==self._PLN_TERRE)):
l0=(1.91944-.004722*t)*sin(m[0])+.02*sin(2*m[0])-.001944*cos(m[0]-m[4])
+.001666*sin(u[0]-l[0])+.001388*sin(4*m[0]-8*m[3]+m[4])-.001388*cos(m[0]-m[2])
l0+=(-.001111*sin(m[0]-m[2])+.001111*cos(4*m[0]-8*m[3]+m[4])+.000833
*sin(2*m[0]-m[2])-.000833*sin(m[4])-.000833*sin(2*m[0]-m[4]))
l0=l[0]+l0*self._DR
b0=0
e=.01675104-.418e-4*t-.126e-6*t*t;
r0=1.0000002*(1+.5*e*e-(e-3./8*e*e*e)*cos(m[0])-.5*e*e*cos(2*m[0])-3.*e*e*e/8*cos(3*m[0]))
if (planete==self._PLN_TERRE):
l0+=self._PI
b0=-b0
elif (planete==self._PLN_MERCURE):
l0=8*t*sin(m[1])+84378*sin(m[1])+10733*sin(2*m[1])+1892*sin(3*m[1])
+381*sin(4*m[1])+83*sin(5*m[1])+19*sin(6*m[1])-646*sin(2*u[1])-306*sin(m[1]-2*u[1])
-274*sin(m[1]-2*u[1])-92*sin(2*m[1]+2*u[1])-28*sin(3*m[1]+2*u[1])+25*sin(2*m[1]-2*u[1])
l0+=(-9*sin(4*m[1]+2*u[1])+7*cos(2*m[1]-5*m[2]))
l0=l[1]+l0/3600*self._DR
b0=24134*sin(u[1])-10*sin(3*u[1])+5180*sin(m[1]-u[1])+4910*sin(m[1]+u[1])
+1124*sin(2*m[1]+u[1])+271*sin(3*m[1]+u[1])+132*sin(2*m[1]-u[1])
+67*sin(4*m[1]+u[1])+18*sin(3*m[1]-u[1])+17*sin(5*m[1]+u[1])-9*sin(m[1]-3*u[1])
b0=b0/3600*self._DR
r0=0.39528-.07834*cos(m[1])-.00795*cos(2*m[1])-.00121*cos(3*m[1])
-.00022*cos(4*m[1])
elif (planete==self._PLN_VENUS):
l0=20*t*sin(m[2])+2814*sin(m[2])+12*sin(2*m[2])-181*sin(2*u[2])
-10*cos(2*m[0]-2*m[2])+7*cos(3*m[0]-3*m[2]);
l0=l[2]+l0/3600*self._DR;
b0=12215*sin(u[2])+83*sin(m[2]+u[2])+83*sin(m[2]-u[2]);
b0=b0/3600*self._DR;
r0=.72335-0.00493*cos(m[2]);
elif (planete==self._PLN_MARS):
l0=37*t*sin(m[3])+4*t*sin(2*m[3])+38451*sin(m[3])+2238*sin(2*m[3])
+181*sin(3*m[3])+17*sin(4*m[3])-52*sin(2*u[3])-22*cos(m[3]-2*m[4])-19*sin(m[3]-m[4])
+17*cos(m[3]-m[4])-16*cos(2*m[3]-2*m[4])+13*cos(m[0]-2*m[3])-10*sin(m[3]-2*u[3]);
l0+=(-10*sin(m[3]+2*u[3])+7*cos(m[0]-m[3])-7*cos(2*m[0]-3*m[3])-5*sin(m[2]-3*m[3])
-5*sin(m[0]-m[3])-5*sin(m[0]-2*m[3])-4*cos(2*m[0]-4*m[3])+4*cos(m[4])+3*cos(m[2]-3*m[3])
+3*sin(2*m[3]-2*m[4]));
l0=l[3]+l0/3600*self._DR
b0=6603*sin(u[3])+622*sin(m[3]-u[3])+615*sin(m[3]+u[3])+64*sin(2*m[3]+u[3])
b0=b0/3600*self._DR
r0=1.53031-.14170*cos(m[3])-.0066*cos(2*m[3])-.00047*cos(3*m[3])
elif (planete==self._PLN_JUPITER):
l0,b0,r0 = self._mc_jd2lbr1c(jj,l,m,u)
elif (planete==self._PLN_SATURNE):
l0,b0,r0 = self._mc_jd2lbr1e(jj,l,m,u)
elif (planete==self._PLN_URANUS):
l0,b0,r0 = self._mc_jd2lbr1f(jj,l,m,u)
elif (planete==self._PLN_NEPTUNE):
l0,b0,r0 = self._mc_jd2lbr1g(jj,l,m,u)
elif (planete==self._PLN_PLUTON):
l0,b0,r0 = self._mc_jd2lbr1h(jj,l,m,u)
elif (planete==self._PLN_LUNE):
l0,b0,r0 = self._mc_jd2lbr1d(jj)
ll = Angle(l0/self._DR)
bb = Angle(b0/self._DR)
coords = Coords((r0, ll, bb))
return coords
def _mc_jd2lbr1c(self, jj, l, m, u):
"""
/***************************************************************************/
/* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr */
/* pour l'equinoxe moyen de la date. (longitudes vraies) */
/* JUPITER */
/***************************************************************************/
/* Les tableaux l, m et u doivent etre dimensiones chacun avec */
/* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a */
/***************************************************************************/
"""
u[0]*=1.;
t=(jj-2415020.0)/36525;
l0=2511+5023*t+19934*sin(m[4])+601*sin(2*m[4])+1093*cos(2*m[4]-5*m[5])
-479*sin(2*m[4]-5*m[5])-185*sin(2*m[4]-2*m[5])+137*sin(3*m[4]-5*m[5])-131*sin(m[4]-2*m[5])
+79*cos(m[4]-m[5])-76*cos(2*m[4]-2*m[5])-74*t*cos(m[4])+68*t*sin(m[4]);
l0+=(+66*cos(2*m[4]-3*m[5])+63*cos(3*m[4]-5*m[5])+53*cos(m[4]-5*m[5])+49*sin(2*m[4]-3*m[5])
-43*t*sin(2*m[4]-5*m[5])-37*cos(m[4])+25*sin(2*l[4])+25*sin(3*m[4])-23*sin(m[4]-5*m[5])
-19*t*cos(2*m[4]-5*m[5])+17*cos(2*m[4]-4*m[5]));
l0+=(+17*cos(3*m[4]-3*m[5])-14*sin(m[4]-m[5]));
l0+=(-13*sin(3*m[4]-4*m[5])-9*cos(2*l[4])+9*cos(m[5])-9*sin(m[5])-9*sin(3*m[4]-2*m[5])
+9*sin(4*m[4]-5*m[5])+9*sin(2*m[4]-6*m[5]+3*m[6])-8*cos(4*m[4]-10*m[5])+7*cos(3*m[4]-4*m[5])
-7*cos(m[4]-3*m[5])-7*sin(4*m[4]-10*m[5]));
l0+=(-7*sin(m[4]-3*m[5])+6*cos(4*m[4]-5*m[5]));
l0+=(-6*sin(3*m[4]-3*m[5])+5*cos(2*m[5])-4*sin(4*m[4]-4*m[5])-4*cos(3*m[5])+4*cos(2*m[4]-m[5])
-4*cos(3*m[4]-2*m[5])-4*t*cos(2*m[4])+3*t*sin(2*m[4])+3*cos(5*m[5])+3*cos(5*m[4]-10*m[5])
+3*sin(2*m[5])-2*sin(2*l[4]-m[4])+2*sin(2*l[4]+m[4]));
l0+=(-2*t*sin(3*m[4]-5*m[5])-2*t*sin(m[4]-5*m[5]));
l0=l[4]+l0/3600*self._DR;
b0=-4692*cos(m[4])+259*sin(m[4])+227-227*cos(2*m[4])+30*t*sin(m[4])+21*t*cos(m[4])
+16*sin(3*m[4]-5*m[5])-13*sin(m[4]-5*m[5])-12*cos(3*m[4])+12*sin(2*m[4])+7*cos(3*m[4]-5*m[5])
-5*cos(m[4]-5*m[5]);
b0=b0/3600*self._DR;
r0=5.20883-.25122*cos(m[4])-.00604*cos(2*m[4])+.0026*cos(2*m[4]-2*m[5])
-.00170*cos(3*m[4]-5*m[5])-.0016*sin(2*m[4]-2*m[5])-.00091*t*sin(m[4])
-.00084*t*cos(m[4])+.00069*sin(2*m[4]-3*m[5])-.00067*sin(m[4]-5*m[5]);
r0+=(.00066*sin(3*m[4]-5*m[5])+.00063*sin(m[4]-m[5])-.00051*cos(2*m[4]-3*m[5])
-.00046*sin(m[4])-.00029*cos(m[4]-5*m[5])+.00027*cos(m[4]-2*m[5])
-.00022*cos(3*m[4])-.00021*sin(2*m[4]-5*m[5]));
ll0=l0
bb0=b0
rr0=r0
return ll0, bb0, rr0
def _mc_jd2lbr1e(self, jj, l, m, u):
"""
/***************************************************************************/
/* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr */
/* pour l'equinoxe moyen de la date. (longitudes vraies) */
/* SATURNE */
/***************************************************************************/
/* Les tableaux l, m et u doivent etre dimensiones chacun avec */
/* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a */
/***************************************************************************/
"""
u[0]*=1.;
t=(jj-2415020.0)/36525;
l0=2507+5014*t+23043*sin(m[5])-2689*cos(2*m[4]-5*m[5])+1177*sin(2*m[4]-5*m[5])-826*cos(2*m[4]-4*m[5])+802*sin(2*m[5])+425*sin(m[4]-2*m[5])-229*t*cos(m[5])-142*t*sin(m[5])-153*cos(2*m[4]-6*m[5])-114*cos(m[5])+101*t*sin(2*m[4]-5*m[5]);
l0+=-70*cos(2*l[5])+67*sin(2*l[5])+66*sin(2*m[4]-6*m[5])+60*t*cos(2*m[4]-5*m[5])+41*sin(m[4]-3*m[5])+39*sin(3*m[5])+31*sin(m[4]-m[5])+31*sin(2*m[4]-2*m[5])-29*cos(2*m[4]-3*m[5])-28*sin(2*m[4]-6*m[5]+3*m[6])+28*cos(m[4]-3*m[5]);
l0+=+22*t*sin(2*m[4]-4*m[5])-22*sin(m[5]-3*m[6])+20*sin(2*m[4]-3*m[5])+20*cos(4*m[4]-10*m[5])+19*cos(2*m[5]-3*m[6])+19*sin(4*m[4]-10*m[5])-17*t*cos(2*m[5])-16*cos(m[5]-3*m[6])-12*sin(2*m[4]-4*m[5])+12*cos(m[4]);
l0+=-12*sin(2*m[5]-2*m[6])-11*t*sin(2*m[5])-11*cos(2*m[4])-12*sin(2*m[5]-2*m[6])-11*t*sin(2*m[5])-11*cos(2*m[4]-7*m[5])+10*sin(2*m[5]-3*m[6])+10*cos(2*m[4]-2*m[5])+9*sin(4*m[4]-9*m[5])-8*sin(m[5]-2*m[6])-8*cos(2*l[5]+m[5]);
l0+=+8*cos(2*l[5]-m[5])+8*cos(m[5]-m[6])-8*sin(2*l[5]-m[5])+7*sin(2*l[5]+m[5])-7*cos(m[4]-2*m[5])-7*cos(2*m[5])-6*t*sin(4*m[4]-10*m[5])+6*t*cos(4*m[4]-10*m[5])+6*t*(2*m[4]-6*m[5])-5*sin(3*m[4]-7*m[5])-5*cos(3*m[4]-3*m[5]);
l0+=-5*cos(2*m[5]-2*m[6])+5*sin(3*m[4]-4*m[5])+5*sin(2*m[4]-7*m[5])+4*sin(3*m[4]-3*m[5])+4*sin(3*m[4]-5*m[5])+4*t*cos(m[4]-3*m[5])+3*t*cos(2*m[4]-4*m[5])+3*cos(2*m[4]-6*m[5]+3*m[6])-3*t*sin(2*l[5]);
l0+=+3*t*cos(2*m[4]-6*m[5])-3*t*cos(2*l[5])+3*cos(3*m[4]-7*m[5])+3*cos(4*m[4]-9*m[5])+3*sin(3*m[4]-6*m[5])+3*sin(2*m[4]-m[5])+3*sin(m[4]-4*m[5])+2*cos(3*m[5]-3*m[6])+2*t*sin(m[4]-2*m[5])+2*sin(4*m[5])-2*cos(3*m[4]-4*m[5])-2*cos(2*m[4]-m[5]);
l0+=-2*sin(2*m[4]-7*m[5]+3*m[6])+2*cos(m[4]-4*m[5])+2*cos(4*m[4]-11*m[5])-2*sin(m[5]-m[6]);
l0=l[5]+l0/3600*self._DR;
b0=185+8297*sin(m[5])-3346*cos(m[5])+462*sin(2*m[5])-189*cos(2*m[5])+79*t*cos(m[5])-71*cos(2*m[4]-4*m[5])+46*sin(2*m[4]-6*m[5])-45*cos(2*m[4]-6*m[5])+29*sin(3*m[5])-20*cos(2*m[4]-3*m[5])+18*t*sin(m[5]);
b0+=-14*cos(2*m[4]-5*m[5])-11*cos(3*m[5])-10*t+9*sin(m[4]-3*m[5])+8*sin(m[4]-m[5])-6*sin(2*m[4]-3*m[5])+5*sin(2*m[4]-7*m[5])-5*cos(2*m[4]-7*m[5])+4*sin(2*m[4]-5*m[5])-4*t*sin(2*m[5])-4*cos(m[4]-m[5])+3*cos(m[4]-3*m[5])+3*t*sin(2*m[4]-4*m[5]);
b0+=+3*sin(m[4]-2*m[5])+2*sin(4*m[5])-2*cos(2*m[4]-2*m[5]);
b0=b0/3600*self._DR;
r0=9.55774-.00028*t-.53252*cos(m[5])-.01878*sin(2*m[4]-4*m[5])-.01482*cos(2*m[5])+.00817*sin(m[4]-m[5])-.00539*cos(m[4]-2*m[5])-.00524*t*sin(m[5])+.00349*sin(2*m[4]-5*m[5])+.00347*sin(2*m[4]-6*m[5]);
r0+=+.00328*t*cos(m[5])-.00225*sin(m[5])+.00149*cos(2*m[4]-6*m[5])-.00126*cos(2*m[4]-2*m[5])+.00104*cos(m[4]-m[5])+.00101*cos(2*m[4]-5*m[5])+.00098*cos(m[4]-3*m[5])-.00073*cos(2*m[4]-3*m[5])-.00062*cos(3*m[5]);
r0+=+.00043*sin(2*m[5]-3*m[6])+.00041*sin(2*m[4]-2*m[5])-.00040*sin(m[4]-3*m[5])+.0004*cos(2*m[4]-4*m[5])-.00023*sin(m[4])+.0002*sin(2*m[4]-7*m[5]);
ll0=l0
bb0=b0
rr0=r0
return ll0, bb0, rr0
def _mc_jd2lbr1f(self, jj, l, m, u):
"""
/***************************************************************************/
/* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr */
/* pour l'equinoxe moyen de la date. (longitudes vraies) */
/* URANUS */
/***************************************************************************/
/* Les tableaux l, m et u doivent etre dimensiones chacun avec */
/* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a */
/***************************************************************************/
"""
u[0]*=1.;
t=(jj-2415020.0)/36525;
l0=32*t*t+19397*sin(m[6])+570*sin(2*m[6])-536*t*cos(m[6])+143*sin(m[5]-2*m[6])+110*t*sin(m[6])+102*sin(m[5]-3*m[6])+76*cos(m[5]-3*m[6])-49*sin(m[4]-m[6])-30*t*cos(2*m[6])+29*sin(2*m[4]-6*m[5]+3*m[6])+29*cos(2*m[6]-2*m[7]);
l0+=-28*cos(m[6]-m[7])+23*sin(3*m[6])-21*cos(m[4]-m[6])+20*sin(m[6]-m[7])+20*cos(m[5]-m[6])-12*t*t*cos(m[6])-12*cos(m[6])+10*sin(2*m[6]-2*m[7])-9*sin(2*u[6])-9*t*t*sin(m[6])+9*cos(2*m[6]-3*m[7])+8*t*cos(m[5]-2*m[6]);
l0+=+7*t*cos(m[5]-3*m[6])-7*t*sin(m[5]-3*m[6])+7*t*sin(2*m[6])+6*sin(2*m[4]-6*m[5]+2*m[6])+6*cos(2*m[4]-6*m[5]+2*m[6])+5*sin(m[5]-4*m[6])-4*sin(3*m[6]-4*m[7])+4*cos(3*m[6]-3*m[7])-3*cos(m[7])-2*sin(m[7]);
l0=l[6]+l0/3600*self._DR;
b0=2775*sin(u[6])+131*sin(m[6]-u[6])+130*sin(m[6]+u[6]);
b0=b0/3600*self._DR;
r0=19.21216-.90154*cos(m[6])-.02488*t*sin(m[6])-.02121*cos(2*m[6])-.00585*cos(m[5]-2*m[6])-.00508*t*cos(m[6])-.00451*cos(m[4]-m[6])+.00336*sin(m[5]-m[6])+.00198*sin(m[4]-m[6])+.00118*cos(m[5]-3*m[6])+.00107*sin(m[5]-2*m[6]);
r0+=-.00103*t*sin(2*m[6])-.00081*cos(3*m[6]-3*m[7]);
ll0=l0
bb0=b0
rr0=r0
return ll0, bb0, rr0
def _mc_jd2lbr1g(self, jj, l, m, u):
"""
/***************************************************************************/
/* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr */
/* pour l'equinoxe moyen de la date. (longitudes vraies) */
/* NEPTUNE */
/***************************************************************************/
/* Les tableaux l, m et u doivent etre dimensiones chacun avec */
/* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a */
/***************************************************************************/
"""
u[0]*=1.;
t=(jj-2415020.0)/36525;
l0=3523*sin(m[7])-50*sin(2*u[7])-43*t*cos(m[7])+29*sin(m[4]-m[7])+19*sin(2*m[7])-18*cos(m[4]-m[7])+13*cos(m[5]-m[7])+13*sin(m[5]-m[7])-9*sin(2*m[6]-3*m[7])+9*cos(2*m[6]-2*m[7])-5*cos(2*m[6]-3*m[7]);
l0+=-4*t*sin(m[7])+4*cos(m[6]-2*m[7])+4*t*t*sin(m[7]);
l0=l[7]+l0/3600*self._DR;
b0=6404*sin(u[7])+55*sin(m[7]+u[7])+55*sin(m[7]-u[7])-33*t*sin(u[7]);
b0=b0/3600*self._DR;
r0=30.07175-.22701*cos(m[7])-.00787*cos(2*l[6]-m[6]-2*l[7])+.00409*cos(m[4]-m[7])-.00314*t*sin(m[7])+.0025*sin(m[4]-m[7])-.00194*sin(m[5]-m[7])+.00185*cos(m[5]-m[7]);
ll0=l0
bb0=b0
rr0=r0
return ll0, bb0, rr0
def _mc_jd2lbr1h(self, jj, l, m, u):
"""
/***************************************************************************/
/* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr */
/* pour l'equinoxe J2000. () */
/* PLUTON */
/***************************************************************************/
/* Les tableaux l, m et u doivent etre dimensiones chacun avec */
/* 9 elements (0 a 8) et sont initialises dans la fonction mc_jd2lbr1a */
/* mais ne servent a rien car on se sert de l'algo de Meeus */
/* Astronomical Algorithms page 247 */
/***************************************************************************/
"""
meeus=[
1,0,0,1,-19798886,19848454,-5453098,-14974876,66867334,68955876,
2,0,0,2,897499,-4955707,3527363,1672673,-11826086,-333765,
3,0,0,3,610820,1210521,-1050939,327763,1593657,-1439953,
4,0,0,4,-341639,-189719,178691,-291925,-18948,482443,
5,0,0,5,129027,-34863,18763,100448,-66634,-85576,
6,0,0,6,-38215,31061,-30594,-25838,30841,-5765,
7,0,1,-1,20349,-9886,4965,11263,-6140,22254,
8,0,1,0,-4045,-4904,310,-132,4434,4443,
9,0,1,1,-5885,-3238,2036,-947,-1518,641,
10,0,1,2,-3812,3011,-2,-674,-5,792,
11,0,1,3,-601,3468,-329,-563,518,518,
12,0,2,-2,1237,463,-64,39,-13,-221,
13,0,2,-1,1086,-911,-94,210,837,-494,
14,0,2,0,595,-1229,-8,-160,-281,616,
15,1,-1,0,2484,-485,-177,259,260,-395,
16,1,-1,1,839,-1414,17,234,-191,-396,
17,1,0,-3,-964,1059,582,-285,-3218,370,
18,1,0,-2,-2303,-1038,-298,692,8019,-7869,
19,1,0,-1,7049,747,157,201,105,45637,
20,1,0,0,1179,-358,304,825,8623,8444,
21,1,0,1,393,-63,-124,-29,-896,-801,
22,1,0,2,111,-268,15,8,208,-122,
23,1,0,3,-52,-154,7,15,-133,65,
24,1,0,4,-78,-30,2,2,-16,1,
25,1,1,-3,-34,-26,4,2,-22,7,
26,1,1,-2,-43,1,3,0,-8,16,
27,1,1,-1,-15,21,1,-1,2,9,
28,1,1,0,-1,15,0,-2,12,5,
29,1,1,1,4,7,1,0,1,-3,
30,1,1,3,1,5,1,-1,1,0,
31,2,0,-6,8,3,-2,-3,9,5,
32,2,0,-5,-3,6,1,2,2,-1,
33,2,0,-4,6,-13,-8,2,14,10,
34,2,0,-3,10,22,10,-7,-65,12,
35,2,0,-2,-57,-32,0,21,126,-233,
36,2,0,-1,157,-46,8,5,270,1068,
37,2,0,0,12,-18,13,16,254,155,
38,2,0,1,-4,8,-2,-3,-26,-2,
39,2,0,2,-5,0,0,0,7,0,
40,2,0,3,3,4,0,1,-11,4,
41,3,0,-2,-1,-1,0,0,4,-14,
42,3,0,-1,6,-3,0,0,18,35,
43,3,0,0,-1,-2,0,1,13,3];
u[0]*=1.;
t=(jj-2451545.0)/36525;
j=34.35+3034.9057*t;
s=50.08+1222.1138*t;
p=238.96+144.9600*t;
l0=0.;
b0=0.;
r0=0.;
for k in range(0,43):
ij=meeus[k*10+1];
it=meeus[k*10+2];
ip=meeus[k*10+3];
a=(self._DR)*(j*ij+s*it+p*ip);
sina=sin(a);
cosa=cos(a);
A=meeus[k*10+4];
B=meeus[k*10+5];
l0+=(A*sina+B*cosa);
A=meeus[k*10+6];
B=meeus[k*10+7];
b0+=(A*sina+B*cosa);
A=meeus[k*10+8];
B=meeus[k*10+9];
r0+=(A*sina+B*cosa);
l0=(238.956785+144.96*t+l0*1e-6)*self._DR;
b0=(-3.908202+b0*1e-6)*self._DR;
r0=40.7247248+r0*1e-7;
ll0=l0
bb0=b0
rr0=r0
return ll0, bb0, rr0
def _mc_jd2lbr1d(self, jj):
"""
/***************************************************************************/
/* Retourne les valeurs de longitude *ll, latitude *bb, rayon vecteur *rr */
/* pour l'equinoxe J2000.0. */
/* LUNE */
/***************************************************************************/
/* J. Meeus, "Astronomical Algorithms", chapter 45, 307-317 */
/***************************************************************************/
"""
arg_lr=[
0,0,1,0,
2,0,-1,0,
2,0,0,0,
0,0,2,0,
0,1,0,0,
0,0,0,2,
2,0,-2,0,
2,-1,-1,0,
2,0,1,0,
2,-1,0,0,
0,1,-1,0,
1,0,0,0,
0,1,1,0,
2,0,0,-2,
0,0,1,2,
0,0,1,-2,
4,0,-1,0,
0,0,3,0,
4,0,-2,0,
2,1,-1,0,
2,1,0,0,
1,0,-1,0,
1,1,0,0,
2,-1,1,0,
2,0,2,0,
4,0,0,0,
2,0,-3,0,
0,1,-2,0,
2,0,-1,2,
2,-1,-2,0,
1,0,1,0,
2,-2,0,0,
0,1,2,0,
0,2,0,0,
2,-2,-1,0,
2,0,1,-2,
2,0,0,2,
4,-1,-1,0,
0,0,2,2,
3,0,-1,0,
2,1,1,0,
4,-1,-2,0,
0,2,-1,0,
2,2,-1,0,
2,1,-2,0,
2,-1,0,-2,
4,0,1,0,
0,0,4,0,
4,-1,0,0,
1,0,-2,0,
2,1,0,-2,
0,0,2,-2,
1,1,1,0,
3,0,-2,0,
4,0,-3,0,
2,-1,2,0,
0,2,1,0,
1,1,-1,0,
2,0,3,0,
2,0,-1,-2]
sinl=[
6288774,
1274027,
658314,
213618,
-185116,
-114332,
58793,
57066,
53322,
45758,
-40923,
-34720,
-30383,
15327,
-12528,
10980,
10675,
10034,
8548,
-7888,
-6766,
-5163,
4987,
4036,
3994,
3861,
3665,
-2689,
-2602,
2390,
-2348,
2236,
-2120,
-2069,
2048,
-1773,
-1595,
1215,
-1110,
-892,
-810,
759,
-713,
-700,
691,
596,
549,
537,
520,
-487,
-399,
-381,
351,
-340,
330,
327,
-323,
299,
294,
0]
cosr=[
-20905355,
-3699111,
-2955968,
-569925,
48888,
-3149,
246158,
-152138,
-170733,
-204586,
-129620,
108743,
104755,
10321,
0,
79661,
-34782,
-23210,
-21636,
24208,
30824,
-8379,
-16675,
-12831,
-10445,
-11650,
14403,
-7003,
0,
10056,
6322,
-9884,
5751,
0,
-4950,
4130,
0,
-3958,
0,
3258,
2616,
-1897,
-2117,
2354,
0,
0,
-1423,
-1117,
-1571,
-1739,
0,
-4421,
0,
0,
0,
0,
1165,
0,
0,
8752]
arg_b=[
0,0,0,1,
0,0,1,1,
0,0,1,-1,
2,0,0,-1,
2,0,-1,1,
2,0,-1,-1,
2,0,0,1,
0,0,2,1,
2,0,1, -1,
0,0,2, -1,
2,-1,0,-1,
2,0,-2,-1,
2,0,1,1,
2,1,0,-1,
2,-1,-1,1,
2,-1,0,1,
2,-1,-1,-1,
0,1,-1,-1,
4,0,-1,-1,
0,1,0,1,
0,0,0,3,
0,1,-1,1,
1,0,0,1,
0,1,1,1,
0,1,1,-1,
0,1,0,-1,
1,0,0,-1,
0,0,3,1,
4,0,0,-1,
4,0,-1,1,
0,0,1,-3,
4,0,-2,1,
2,0,0,-3,
2,0,2,-1,
2,-1,1,-1,
2,0,-2,1,
0,0,3,-1,
2,0,2,1,
2,0,-3,-1,
2,1,-1,1,
2,1,0,1,
4,0,0,1,
2,-1,1,1,
2,-2,0,-1,
0,0,1,3,
2,1,1,-1,
1,1,0,-1,
1,1,0,1,
0,1,-2,-1,
2,1,-1,-1,
1,0,1,1,
2,-1,-2,-1,
0,1,2,1,
4,0,-2,-1,
4,-1,-1,-1,
1,0,1,-1,
4,0,1,-1,
1,0,-1,-1,
4,-1,0,-1,
2,-2,0,1]
sinb=[
5128122,
280602,
277693,
173237,
55413,
46271,
32573,
17198,
9266,
8822,
8216,
4324,
4200,
-3359,
2463,
2211,
2065,
-1870,
1828,
-1794,
-1749,
-1565,
-1491,
-1475,
-1410,
-1344,
-1335,
1107,
1021,
833,
777,
671,
607,
596,
491,
-451,
439,
422,
421,
-366,
-351,
331,
315,
302,
-283,
-229,
223,
223,
-220,
-220,
-185,
181,
-177,
176,
166,
-164,
132,
-119,
115,
107]
T=(jj-2451545.)/36525.;
lp=(218.3164591+481267.88134236*T-.0013268*T*T+T*T*T/538841.-T*T*T*T/65194000)*(self._DR);
d=(297.8502042+445267.1115168*T-.00016300*T*T+T*T*T/545868-T*T*T*T/113065000)*(self._DR);
m=(357.5291092+35999.0502909*T-.0001536*T*T+T*T*T*T/24490000)*(self._DR);
mp=(134.9634114+477198.8676313*T+.0089970*T*T+T*T*T/69699.-T*T*T*T/14712000)*(self._DR);
f=(93.2720993+483202.0175273*T-.0034029*T*T+T*T*T/3526000+T*T*T*T/863310000)*(self._DR);
a1=(119.75+131.849*T)*(self._DR);
a2=(53.09+479264.290*T)*(self._DR);
a3=(313.45+481266.484*T)*(self._DR);
e=1-.002516*T-.0000074*T*T;
e2=e*e;
# --- longitude & radius ---
l=0.;
r=0.;
for k in range(0,60):
xe=1.;
if (fabs(arg_lr[k*4+1])==1):
xe=e;
elif (fabs(arg_lr[k*4+1])==2):
xe=e2
angle=1.*arg_lr[k*4+0]*d+arg_lr[k*4+1]*m+arg_lr[k*4+2]*mp+arg_lr[k*4+3]*f
sina=sin(angle)
cosa=cos(angle)
l+=sinl[k]*sina*xe
r+=cosr[k]*cosa*xe
l+=3958.*sin(a1)+1962.*sin(lp-f)+318.*sin(a2)
# --- latitude ---
b=0.
for k in range(0,60):
xe=1.;
if (fabs(arg_b[k*4+1])==1):
xe=e
elif (fabs(arg_b[k*4+1])==2):
xe=e2
angle=1.*arg_b[k*4+0]*d+arg_b[k*4+1]*m+arg_b[k*4+2]*mp+arg_b[k*4+3]*f
sina=sin(angle)
b+=sinb[k]*sina*xe
b+=-2235.*sin(lp)+382.*sin(a3)+175.*sin(a1-f)+175.*sin(a1+f)+127.*sin(lp-mp)-115.*sin(lp+mp)
l=lp+(l*1.0e-6)*(self._DR)
b=(b*1.0e-6)*(self._DR)
r=(385000.56e3+r)/(self._CST_UA)
ll0=l;
bb0=b;
rr0=r;
return ll0, bb0, rr0
# ========================================================
# === internal methods : Coordinate transformations
# ========================================================
def _mc_obliqmoy(self, jj):
"""
/***************************************************************************/
/* Retourne la valeur de l'obliquite terrestre moyenne pour jj */
/***************************************************************************/
/* formule de Laskar (JM) */
/***************************************************************************/
"""
t=(jj-2451545.0)/36525
u=t/100
eps0=u*(-4680.93-u*(1.55+u*(1999.25-u*(51.38-u*(249.67-u*(39.05+u*(7.12+u*(27.87+u*(5.79+u*(2.45))))))))))
eps0=(23.4392911111+eps0/3600)*self._DR
return eps0
def _mc_xyzec2eq(self, xec, yec, zec, eps):
"""
/***************************************************************************/
/* Transforme les coord. cart. ecliptiques vers equatoriales */
/***************************************************************************/
/***************************************************************************/
"""
xeq0=xec
yeq0=yec*cos(eps)-zec*sin(eps)
zeq0=yec*sin(eps)+zec*cos(eps)
return xeq0, yeq0, zeq0
def _mc_xyzeq2ec(self, xeq, yeq, zeq, eps):
"""
/***************************************************************************/
/* Transforme les coord. cart. equatoriales vers ecliptiques */
/***************************************************************************/
/***************************************************************************/
"""
eps=-eps
xec0=xeq
yec0=yeq*cos(eps)-zeq*sin(eps)
zec0=yeq*sin(eps)+zeq*cos(eps)
return xec0, yec0, zec0
def _mc_paraldxyzeq(self, jj, longuai, rhocosphip, rhosinphip):
"""
/***************************************************************************/
/* Calcul des corrections cartesiennes equatoriales de la parallaxe */
/* Xtopo = Xgeo - *dxeq etc... */
/***************************************************************************/
/* ref : Danby J.M.A. "Fundamentals of Celestial Mechanics" 2nd ed. 1992 */
/* Willmann Bell */
/* Formules (6.17.1) pp 208 */
/***************************************************************************/
"""
cst=(self._CST_EARTH_SEMI_MAJOR_RADIUS)/(self._CST_UA); # equatorial radius of the Earth in U.A.
tsl = self._mc_tsl(jj,-longuai)
dxeq=(cst*rhocosphip*cos(tsl))
dyeq=(cst*rhocosphip*sin(tsl))
dzeq=(cst*rhosinphip)
return dxeq, dyeq, dzeq
def _mc_latalt2rhophi(self, latitude, altitude):
"""
/***************************************************************************/
/* Retourne les valeurs de rhocosphi' et rhosinphi' (en rayons equatorial */
/* terrestres) a partir de la latitude et de l'altitude. */
/***************************************************************************/
/* Latitude en degres decimaux. */
/* Altitude en metres. */
/* Algo : Meeus "Astronomical Algorithms" p78 */
/***************************************************************************/
"""
lat=latitude*(self._DR);
alt=altitude;
aa=self._CST_EARTH_SEMI_MAJOR_RADIUS;
ff=1./self._CST_EARTH_INVERSE_FLATTENING;
bb=aa*(1-ff);
u=atan(bb/aa*tan(lat));
a=bb/aa*sin(u)+alt/aa*sin(lat);
b= cos(u)+alt/aa*cos(lat);
rhocosphip=b;
rhosinphip=a;
return rhocosphip, rhosinphip
def _mc_rhophi2latalt(self, rhosinphip,rhocosphip):
"""
/***************************************************************************/
/* Retourne les valeurs de la latitude et de l'altitude a partir de */
/* rhocosphi' et rhosinphi' (en rayons equatorial terrestres) */
/***************************************************************************/
/* Latitude en degres decimaux. */
/* Altitude en metres. */
/* Algo : Meeus "Astronomical Algorithms" p78 */
/***************************************************************************/
"""
aa=self._CST_EARTH_SEMI_MAJOR_RADIUS;
ff=1./self._CST_EARTH_INVERSE_FLATTENING;
bb=aa*(1-ff);
rho=sqrt(rhosinphip*rhosinphip+rhocosphip*rhocosphip);
if (rho==0.):
latitude=0.;
altitude=-aa;
return latitude, altitude
phip=atan2(rhosinphip,rhocosphip)
phi0=atan(aa*aa/bb/bb*tan(phip)); # alt=0
u0=atan(bb/aa*tan(phi0));
sinu0=sin(u0);
cosu0=cos(u0);
sinphi0=sin(phi0);
cosphi0=cos(phi0);
alt=-1000
while (alt<20000.):
rhosinphip0 = bb/aa*sinu0 + alt/aa*sinphi0 ;
rhocosphip0 = cosu0 + alt/aa*cosphi0 ;
rho0=sqrt(rhosinphip0*rhosinphip0+rhocosphip0*rhocosphip0);
if ((rho-rho0)<0):
break;
alt += 0.1
lat=phi0;
alt-=0.1;
latitude=lat/(self._DR);
altitude=alt;
return latitude, altitude
def _mc_nutation(self, jj, precision):
"""
/***************************************************************************/
/* Retourne la valeur de la nutation pour jj donne */
/***************************************************************************/
/* precision=0 pour une incertitude de 0"5 sur dpsi 0"1 sur deps */
/* precision=1 pour de la tres haute precision (environ 0.01") */
/***************************************************************************/
"""
t=(jj-2451545.0)/36525;
o=125.04452-1934.136261*t+.0020708*t*t+t*t*t/450000;
o=fmod(o*self._DR,2*self._PI);
if (precision==0):
l=280.4665+36000.7698*t;
lp=218.3165+481267.8813*t;
l=fmod(l*self._DR,2*self._PI);
lp=fmod(lp*self._DR,2*self._PI);
dpsi0=-17.20*sin(o)-1.32*sin(2*l)-.23*sin(2*lp)+.21*sin(2*o);
dpsi0=dpsi0/3600*self._DR;
deps0=9.20*cos(o)+.57*cos(2*l)+.10*cos(2*lp)-.09*cos(2*o);
deps0=deps0/3600*self._DR;
else:
# if (precision==1)
d=297.85036+445267.111480*t-.0019142*t*t+t*t*t/189474
m=357.52772+35999.050340*t-.0001603*t*t-t*t*t/300000
mp=134.96298+477198.867398*t+.0086972*t*t+t*t*t/56250
f=93.27191+483202.017538*t-.0036825*t*t+t*t*t/327270
d=fmod(d*self._DR,2*self._PI)
m=fmod(m*self._DR,2*self._PI)
mp=fmod(mp*self._DR,2*self._PI)
f=fmod(f*self._DR,2*self._PI)
dpsi0=(-171996-174.2*t)*sin(o)+(-13187-1.6*t)*sin(-2*d+2*f+2*o)
+(-2274-.02*t)*sin(2*f+2*o)+(2062+.2*t)*sin(2*o)
+(1426-3.4*t)*sin(m)+(712+.1*t)*sin(mp)
dpsi0+=((-517+1.2*t)*sin(-2*d+m+2*f+2*o)-(386-.4*t)*sin(2*f+o)
-301*sin(mp+2*f+2*o)+(217-.5*t)*sin(-2*d-m+2*f+2*o)
-158*sin(-2*d+mp)+(129+.1*t)*sin(-2*d+2*f+o))
dpsi0+=(123*sin(-mp+2*f+2*o)+63*sin(2*d)+(63+.1*t)*sin(mp+o)
-59*sin(2*d-mp+2*f+2*o)+(-58-.1*t)*sin(-mp+o)-51*sin(mp+2*f+o)
+48*sin(-2*d+2*mp)+46*sin(-2*mp+2*f+o)-38*sin(2*d+2*f+2*o))
dpsi0=dpsi0*1e-4/3600*self._DR
deps0=(92025+8.9*t)*cos(o)+(5736-3.1*t)*cos(-2*d+2*f+2*o)
+(977-.5*t)*cos(2*f+2*o)+(-895+.5*t)*cos(2*o)+(54-.1*t)*cos(m)
-7*cos(mp)+(224-.6*t)*cos(-2*d+m+2*f+2*o)+200*cos(2*f+o);
deps0+=((129-.1*t)*cos(mp+2*f+2*o)+(-95+.3*t)*cos(-2*d-m+2*f+2*o)
-70*cos(-2*d+2*f+o)-53*cos(-mp+2*f+2*o)-33*cos(mp+o)
+26*cos(2*d-mp+2*f+2*o)+32*cos(-mp+o)+27*cos(mp+2*f+o))
deps0=deps0*1e-4/3600*self._DR
dpsi=dpsi0;
deps=deps0;
return dpsi, deps
def _mc_nutradec(self, jj, asd1, dec1, signe):
"""
/***************************************************************************/
/* Corrige asd1,dec1 de la nutation et retourne asd2 et dec2 */
/***************************************************************************/
/* Trueblood & Genet : Telescop Control ed. Willmann Bell (1997) p71 */
/***************************************************************************/
"""
methode=1;
# --- obliquite moyenne --- */
eps = self._mc_obliqmoy(jj)
# --- nutation ---*/
dpsi, deps = self._mc_nutation(jj,1)
if (fabs(dec1) >= (self._PI)/2.):
asd2=asd1
dec2=dec1
else:
if (methode==0):
tand=tan(dec1);
dasd=(cos(eps)+sin(eps)*sin(asd1)*tand)*dpsi-cos(asd1)*tand*deps;
ddec=sin(eps)*cos(asd1)*dpsi+sin(asd1)*deps;
dasd*=signe;
ddec*=signe;
asd1+=dasd;
dec1+=ddec;
else:
# eq->ecl
l=asd1 ; b=dec1 ; r=1;
bb = Angle(b/self._DR)
ll = Angle(l/self._DR)
coords = Coords((r,ll,bb))
xeq,yeq,zeq = coords.cart()
xec,yec,zec = self._mc_xyzeq2ec(xeq,yeq,zeq,eps)
coords = Coords((xec,yec,zec))
r,l,b = coords.sphe("rad","rad")
l=l+dpsi*signe;
eps=eps+deps*signe;
bb = Angle(b/self._DR)
ll = Angle(l/self._DR)
coords = Coords((r,ll,bb))
xeq,yeq,zeq = coords.cart()
# ecl->eq
xeq,yeq,zeq = self._mc_xyzec2eq(xec,yec,zec,eps)
coords = Coords((xeq,yeq,zeq))
r,asd1,dec1 = coords.sphe("rad","rad")
asd1=fmod(4*self._PI+asd1,2*self._PI);
asd2=asd1;
dec2=dec1;
return asd2, dec2
def _mc_he2ge(self, xh,yh,zh,xs,ys,zs):
"""
/***************************************************************************/
/* Translation du repere heliocentrique en geocentrique (cartesien) */
/***************************************************************************/
/* ENTREES : */
/* xh,yh,zh : coordonnees heliocentriques de l'astre */
/* xs,ys,zs : coordonnees geocentriques du Soleil */
/* SORTIES : */
/* *xg,*yg,*zg : coordonnees geocentriques de l'astre */
/***************************************************************************/
"""
xg=xh+xs
yg=yh+ys
zg=zh+zs
return xg, yg, zg
def _mc_aberpla(self, jj1, delta):
"""
/***************************************************************************/
/* Corrige jj1 de l'aberration de la lumiere et retourne jj2 */
/***************************************************************************/
/***************************************************************************/
"""
jj2=jj1-0.0057755*delta
return jj2
def _mc_parallaxe_stellaire(self, jj, asd1, dec1, plx_mas):
"""
/***************************************************************************/
/* Corrige asd1,dec1 de la parallaxe stellaire et retourne asd2 et dec2 */
/***************************************************************************/
/* A. Danjon : Astronomie Generale ed. A. Blanchard (1980) p130 */
/***************************************************************************/
"""
afm = self._AFM_VFP79
# --- obliquite moyenne --- */
eps = self._mc_obliqmoy(jj)
# --- longitude vraie du soleil ---*/
llp, mmp, uup = self._mc_jd2lbr1a(jj);
coords = self._mc_jd2lbr1b(jj, self._PLN_SOLEIL, llp, mmp, uup, afm)
rs,ls,bs = coords.sphe("rad","rad")
dpsi, deps = self._mc_nutation(jj,1)
ls += dpsi
plxrad=plx_mas*1e-3/3600.*(self._DR);
secd=cos(dec1);
if (secd==0):
asd2=asd1;
dec2=dec1;
else:
secd=1./secd;
dasd=(cos(eps)*cos(asd1)*sin(ls)-sin(asd1)*cos(ls))*secd;
ddec=(sin(eps)*cos(dec1)*sin(ls)-sin(dec1)*cos(asd1)*cos(ls)-cos(eps)*sin(dec1)*sin(asd1)*sin(ls));
asd1+=plxrad*dasd;
dec1+=plxrad*ddec;
asd1=fmod(4*self._PI+asd1,2*self._PI);
asd2=asd1;
dec2=dec1;
return asd2, dec2
def _mc_aberration_annuelle(self, jj, asd1, dec1, signe):
"""
/***************************************************************************/
/* Corrige asd1,dec1 de l'aberration annuelle et retourne asd2 et dec2 */
/***************************************************************************/
/* Trueblood & Genet : Telescop Control ed. Willmann Bell (1997) p81-82 */
/* Formule sans les E-terms */
/***************************************************************************/
"""
k=20.49552; # constant of annual aberration
# --- computation method
afm = self._AFM_VFP79
# --- obliquite moyenne ---
eps = self._mc_obliqmoy(jj)
# --- longitude vraie du soleil ---
llp, mmp, uup = self._mc_jd2lbr1a(jj);
coords = self._mc_jd2lbr1b(jj, self._PLN_SOLEIL, llp, mmp, uup, afm)
rs,ls,bs = coords.sphe("rad","rad")
dpsi, deps = self._mc_nutation(jj,1)
ls += dpsi;
# ---
secd=cos(dec1)
if (secd==0):
asd2=asd1
dec2=dec1
return asd2, dec2
secd=1./secd
if (fabs(secd) < 100):
c = cos(asd1)*secd;
d = sin(asd1)*secd;
cp = tan(eps)*cos(dec1)-sin(asd1)*sin(dec1);
dp = cos(asd1)*sin(dec1);
cc = -k*cos(eps)*cos(ls);
dd = -k*sin(ls);
dasd = (cc*c+dd*d)/3600.*(self._DR);
ddec = (cc*cp+dd*dp)/3600.*(self._DR);
dasd *= float(signe);
ddec *= float(signe);
asd1 += dasd;
dec1 += ddec;
asd1 = fmod(4*self._PI+asd1, 2*self._PI);
asd2 = asd1;
dec2 = dec1;
return asd2, dec2
def _mc_aberration_diurne(self, jj, asd1, dec1, longuai, rhocosphip, rhosinphip, signe):
"""
/***************************************************************************/
/* Corrige asd1,dec1 de l'aberration diurne et retourne asd2 et dec2 */
/***************************************************************************/
/* Trueblood & Genet : Telescop Control ed. Willmann Bell (1997) p83-84 */
/***************************************************************************/
"""
a=(self._CST_EARTH_SEMI_MAJOR_RADIUS)*1e-3;
tsl = self._mc_tsl(jj,-longuai)
h=tsl-asd1;
latitude, altitude = self._mc_rhophi2latalt(rhosinphip, rhocosphip)
if (rhocosphip==0.):
if (rhosinphip>0):
phip=(self._PI)/2.
else:
phip=-(self._PI)/2.
else:
phip=atan2(rhosinphip,rhocosphip)
sinphi=sin(latitude)
r=(21*sinphi*sinphi+a)/a;
secd=cos(dec1);
if (secd==0):
asd2=asd1
dec2=dec1
else:
secd=1./secd;
if (fabs(secd)<100):
dasd=(0.320*r*cos(phip)*cos(h)*secd)/3600.*(self._DR);
ddec=(0.320*r*cos(phip)*sin(h)*sin(dec1))/3600.*(self._DR);
else:
dasd=0
ddec=0
dasd*=signe;
ddec*=signe;
asd1+=dasd;
dec1+=ddec;
asd1=fmod(4*self._PI+asd1,2*self._PI);
asd2=asd1;
dec2=dec1;
return asd2, dec2
def _mc_precad(self, jd1, asd1, dec1, jd2):
"""
/***************************************************************************/
/* Passage des coordonnees spheri. equatoriales d'un equinoxe a un autre */
/***************************************************************************/
/***************************************************************************/
"""
t=(jd2-jd1)/36525.
tt=(jd1-2451545.0)/36525.
dz=(2306.2181+1.39656*tt-.000139*tt*tt)*t+(0.30188-.000344*tt)*t*t+.017998*t*t*t
dz=dz/3600*self._DR
zz=(2306.2181+1.39656*tt-.000139*tt*tt)*t+(1.09468+.000066*tt)*t*t+.018203*t*t*t
zz=zz/3600*self._DR
th=(2004.3109-0.85330*tt-.000217*tt*tt)*t-(0.42665+.000217*tt)*t*t-.041833*t*t*t
th=th/3600*self._DR
cosasddz=cos(asd1+dz)
sinasddz=sin(asd1+dz)
costh=cos(th)
sinth=sin(th)
cosdec=cos(dec1)
sindec=sin(dec1)
a=cosdec*sinasddz
b=costh*cosdec*cosasddz-sinth*sindec
c=sinth*cosdec*cosasddz+costh*sindec
a=atan2(a,b)+zz
asd2=fmod(4*self._PI+a,2*self._PI)
dec2=asin(c)
return asd2, dec2
def _mc_elonphas(self, r, rsol, delta):
"""
/***************************************************************************/
/* Calcul des angles d'elongation et de phase. */
/***************************************************************************/
/***************************************************************************/
"""
elong=0.
phase=0.
if (delta!=0.):
if (rsol!=0.):
elong=acos((rsol*rsol+delta*delta-r*r)/(2*rsol*delta))
if (r!=0.):
phase=acos((r*r+delta*delta-rsol*rsol)/(2*r*delta))
return elong, phase
def _mc_elonphaslimb(self, asd, dec, asds, decs, r, delta):
"""
/***************************************************************************/
/* Calcul des angles d'elongation de phase et de position du limbe. */
/***************************************************************************/
/* Meeus page 316 (46.2) (46.3) (46.5) */
/* asd,dec : planete. */
/* asds,decs : soleil */
/* r : distance Terre-Soleil (corrigee de l'abberation). */
/* delta : distance Terre-Planete (corrigee de l'abberation). */
/***************************************************************************/
"""
elong=0.
phase=0.
elong=acos(sin(decs)*sin(dec)+cos(decs)*cos(dec)*cos(asds-asd))
i=atan2(r*sin(elong),delta-r*cos(elong))
if (i<0):
i=i+self._PI
phase=i;
i=atan2(cos(decs)*sin(asds-asd),sin(decs)*cos(dec)-cos(decs)*sin(dec)*cos(asds-asd));
if (i<0):
i=2*self._PI+i
posang_brightlimb=fmod(4*self._PI+i,2*self._PI)
return elong, phase, posang_brightlimb
def _mc_hd2ah(self, ha, dec, latitude):
"""
/***************************************************************************/
/* Transforme les coord. sph. equatoriales vers sph. azinuth hauteur */
/***************************************************************************/
/* ha est l'angle horaire local. */
/* Tout en radian */
/***************************************************************************/
"""
if (dec>=self._PISUR2):
aa=(self._PI)
hh=latitude
elif (dec<=-self._PISUR2):
aa=0.
hh=-latitude
else:
aa=atan2(sin(ha),cos(ha)*sin(latitude)-tan(dec)*cos(latitude))
hh=asin(sin(latitude)*sin(dec)+cos(latitude)*cos(dec)*cos(ha))
az=fmod(4*self._PI+aa,2*self._PI);
h=hh;
return az, h
def _mc_ad2hd(self, jd, longuai, asd):
"""
/***************************************************************************/
/* Transforme l'ascension droite en angle horaire */
/***************************************************************************/
/***************************************************************************/
"""
# --- calcul du TSL en radians ---*/
tsl = self._mc_tsl(jd,-longuai)
h=tsl-asd;
ha=fmod(4*self._PI+h,2*self._PI);
return ha
def _mc_ah2hd(self, az, h, latitude):
"""
/***************************************************************************/
/* Transforme les coord. sph. azinuth hauteur vers sph. equatoriales */
/***************************************************************************/
/* ha est l'angle horaire local. */
/***************************************************************************/
"""
if (h>=self._PISUR2):
ahh=0.
decc=latitude
elif (h<=-self._PISUR2):
ahh=0.
decc=-latitude
else:
ahh=atan2(sin(az),cos(az)*sin(latitude)+tan(h)*cos(latitude))
decc=asin(sin(latitude)*sin(h)-cos(latitude)*cos(h)*cos(az))
ha=fmod(4*self._PI+ahh,2*self._PI)
dec=decc
return ha, dec
def _mc_hd2ad(self, jd, longuai, ha):
"""
/***************************************************************************/
/* Transforme l'angle horaire en ascension droite */
/***************************************************************************/
/***************************************************************************/
"""
#--- calcul du TSL en radians ---*/
tsl = self._mc_tsl(jd,-longuai)
ra=tsl-ha
asd=fmod(4*self._PI+ra,2*self._PI)
return asd
def _mc_hd2parallactic(self, ha, dec, latitude):
"""
/***************************************************************************/
/* Transforme les coord. sph. equatoriales vers angle parallactic */
/***************************************************************************/
/* ha est l'angle horaire local. */
/***************************************************************************/
"""
if (fabs(latitude)>=self._PISUR2):
aa=0
elif ((ha==0.) and (dec==latitude)):
aa=0.
else:
aa=atan2(sin(ha),cos(dec)*tan(latitude)-sin(dec)*cos(ha))
parallactic=aa;
return parallactic
# ========================================================
# === internal methods : Refraction
# ========================================================
# ========================================================
# === internal methods : Time
# ========================================================
def _mc_tsl(self, jj, longitude):
"""
/***************************************************************************/
/* Calcul du temps sideral local apparent (en radian) */
/* La longitude est comptee en radian positive vers l'ouest */
/* jj = UT1 (on a toujours |UTC-UT1|<1 sec) */
/***************************************************************************/
/* Formules (11.1) et (11.4) de Meeus */
/***************************************************************************/
"""
j=(jj-2451545.0)
t=j/36525
theta0=280.460618375+360.98564736629*j+.000387933*t*t-t*t*t/38710000
eps = self._mc_obliqmoy(jj)
dpsi, deps = self._mc_nutation(jj,1)
theta0+=(dpsi*cos(eps+deps)/(self._DR))
theta0-=longitude/(self._DR)
theta0=fmod(theta0,360.)
theta0=fmod(theta0+720.,360.)*self._DR
tsl=theta0
return tsl
def _mc_td2tu(self, jjtd):
"""
/***************************************************************************/
/* Retourne la valeur de jj en TU a partir de jj en temps dynamique */
/***************************************************************************/
/* Algo : Meeus "Astronomical Algorithms" p73 (9.1) */
/***************************************************************************/
"""
dt=0.
dt = self._mc_tdminusut(jjtd)
jjtu=jjtd-dt/86400.
return jjtu
def _mc_tu2td(self, jjtu):
"""
/***************************************************************************/
/* Retourne la valeur de jj en temps dynamique a partir de jj en TU */
/* UTC -> TT */
/***************************************************************************/
/* Algo : Meeus "Astronomical Algorithms" p73 (9.1) */
/***************************************************************************/
"""
dt=0.;
dt = self._mc_tdminusut(jjtu)
jjtd=jjtu+dt/86400.
return jjtd
def _mc_tdminusut(self, jj):
"""
/***************************************************************************/
/* Retourne la valeur dt(sec)=TT-UT partir de jj en TU */
/***************************************************************************/
/* Algo : Meeus "Astronomical Algorithms" p72 */
/* and ftp://maia.usno.navy.mil/ser7/tai-utc.dat */
/* ET 1960-1983 */
/* TDT 1984-2000 */
/* UTC 1972- GPS 1980- TAI 1958- TT 2001- */
/*----+---------+-------------+-------------------------+----- */
/* | | | | */
/* |<------ TAI-UTC ------>|<----- TT-TAI ----->| */
/* | | | 32.184s fixed | */
/* |<GPS-UTC>|<- TAI-GPS ->| | */
/* | | 19s fixed | | */
/* | | */
/* <> delta-UT = UT1-UTC | */
/* | (max 0.9 sec) | */
/*-----+------------------------------------------------+----- */
/* |<-------------- delta-T = TT-UT1 -------------->| */
/* UT1 (UT) TT/TDT/ET */
/* */
/* http://stjarnhimlen.se/comp/time.html */
/***************************************************************************/
"""
ds=0.
t=0.
indexmax=201
table = []
table.append(2312752.5) ; table.append( +124.) ; # 1620
table.append(2313483.5) ; table.append( +115.)
table.append(2314213.5) ; table.append( +106.)
table.append(2314944.5) ; table.append( +98.)
table.append(2315674.5) ; table.append( +91.)
table.append(2316405.5) ; table.append( +85.)
table.append(2317135.5) ; table.append( +79.)
table.append(2317866.5) ; table.append( +74.)
table.append(2318596.5) ; table.append( +70.)
table.append(2319327.5) ; table.append( +65.)
table.append(2320057.5) ; table.append( +62.)
table.append(2320788.5) ; table.append( +58.)
table.append(2321518.5) ; table.append( +55.)
table.append(2322249.5) ; table.append( +53.)
table.append(2322979.5) ; table.append( +50.)
table.append(2323710.5) ; table.append( +48.)
table.append(2324440.5) ; table.append( +46.)
table.append(2325171.5) ; table.append( +44.)
table.append(2325901.5) ; table.append( +42.)
table.append(2326632.5) ; table.append( +40.)
table.append(2327362.5) ; table.append( +37.)
table.append(2328093.5) ; table.append( +35.)
table.append(2328823.5) ; table.append( +33.)
table.append(2329554.5) ; table.append( +31.)
table.append(2330284.5) ; table.append( +28.)
table.append(2331015.5) ; table.append( +26.)
table.append(2331745.5) ; table.append( +24.)
table.append(2332476.5) ; table.append( +22.)
table.append(2333206.5) ; table.append( +20.)
table.append(2333937.5) ; table.append( +18.)
table.append(2334667.5) ; table.append( +16.)
table.append(2335398.5) ; table.append( +14.)
table.append(2336128.5) ; table.append( +13.)
table.append(2336859.5) ; table.append( +12.)
table.append(2337589.5) ; table.append( +11.)
table.append(2338320.5) ; table.append( +10.)
table.append(2339050.5) ; table.append( +9.)
table.append(2339781.5) ; table.append( +9.)
table.append(2340511.5) ; table.append( +9.)
table.append(2341242.5) ; table.append( +9.)
table.append(2341972.5) ; table.append( 9.) ; # 1700
table.append(2342702.5) ; table.append( 9.) ; # 1702
table.append(2343432.5) ; table.append( 9.) ; # 1704
table.append(2344163.5) ; table.append( 9.) ; # 1706
table.append(2344893.5) ; table.append( 10.) ; # 1708
table.append(2345624.5) ; table.append( 10.) ; # 1710
table.append(2346354.5) ; table.append( 10.) ; # 1712
table.append(2347085.5) ; table.append( 10.) ; # 1714
table.append(2347815.5) ; table.append( 10.) ; # 1716
table.append(2348546.5) ; table.append( 11.) ; # 1718
table.append(2349276.5) ; table.append( 11.) ; # 1720
table.append(2350007.5) ; table.append( 11.) ; # 1722
table.append(2350737.5) ; table.append( 11.) ; # 1724
table.append(2351468.5) ; table.append( 11.) ; # 1726
table.append(2352198.5) ; table.append( 11.) ; # 1728
table.append(2352929.5) ; table.append( 11.) ; # 1730
table.append(2353659.5) ; table.append( 11.) ; # 1732
table.append(2354390.5) ; table.append( 12.) ; # 1734
table.append(2355120.5) ; table.append( 12.) ; # 1736
table.append(2355851.5) ; table.append( 12.) ; # 1738
table.append(2356581.5) ; table.append( 12.) ; # 1740
table.append(2357312.5) ; table.append( 12.) ; # 1742
table.append(2358042.5) ; table.append( 13.) ; # 1744
table.append(2358773.5) ; table.append( 13.) ; # 1746
table.append(2359503.5) ; table.append( 13.) ; # 1748
table.append(2360234.5) ; table.append( 13.) ; # 1750
table.append(2360964.5) ; table.append( 14.) ; # 1752
table.append(2361695.5) ; table.append( 14.) ; # 1754
table.append(2362425.5) ; table.append( 14.) ; # 1756
table.append(2363156.5) ; table.append( 15.) ; # 1758
table.append(2363886.5) ; table.append( 15.) ; # 1760
table.append(2364617.5) ; table.append( 15.) ; # 1762
table.append(2365347.5) ; table.append( 15.) ; # 1764
table.append(2366078.5) ; table.append( 16.) ; # 1766
table.append(2366808.5) ; table.append( 16.) ; # 1768
table.append(2367539.5) ; table.append( 16.) ; # 1770
table.append(2368269.5) ; table.append( 16.) ; # 1772
table.append(2369000.5) ; table.append( 16.) ; # 1774
table.append(2369730.5) ; table.append( 17.) ; # 1776
table.append(2370461.5) ; table.append( 17.) ; # 1778
table.append(2371191.5) ; table.append( 17.) ; # 1780
table.append(2371922.5) ; table.append( 17.) ; # 1782
table.append(2372652.5) ; table.append( 17.) ; # 1784
table.append(2373383.5) ; table.append( 17.) ; # 1786
table.append(2374113.5) ; table.append( 17.) ; # 1788
table.append(2374844.5) ; table.append( 17.) ; # 1790
table.append(2375574.5) ; table.append( 16.) ; # 1792
table.append(2376305.5) ; table.append( 16.) ; # 1794
table.append(2377035.5) ; table.append( 15.) ; # 1796
table.append(2377766.5) ; table.append( 14.) ; # 1798
table.append(2378496.5) ; table.append( 13.7) ; # 1800
table.append(2379226.5) ; table.append( 13.1) ; # 1802
table.append(2379956.5) ; table.append( 12.7) ; # 1804
table.append(2380687.5) ; table.append( 12.5) ; # 1806
table.append(2381417.5) ; table.append( 12.5) ; # 1808
table.append(2382148.5) ; table.append( 12.5) ; # 1810
table.append(2382878.5) ; table.append( 12.5) ; # 1812
table.append(2383609.5) ; table.append( 12.5) ; # 1814
table.append(2384339.5) ; table.append( 12.5) ; # 1816
table.append(2385070.5) ; table.append( 12.3) ; # 1818
table.append(2385800.5) ; table.append( 12.0) ; # 1820
table.append(2386531.5) ; table.append( 11.4) ; # 1822
table.append(2387261.5) ; table.append( 10.6) ; # 1824
table.append(2387992.5) ; table.append( 9.6) ; # 1826
table.append(2388722.5) ; table.append( 8.6) ; # 1828
table.append(2389453.5) ; table.append( 7.5) ; # 1830
table.append(2390183.5) ; table.append( 6.6) ; # 1832
table.append(2390914.5) ; table.append( 6.0) ; # 1834
table.append(2391644.5) ; table.append( 5.7) ; # 1836
table.append(2392375.5) ; table.append( 5.6) ; # 1838
table.append(2393105.5) ; table.append( 5.7) ; # 1840
table.append(2393836.5) ; table.append( 5.9) ; # 1842
table.append(2394566.5) ; table.append( 6.2) ; # 1844
table.append(2395297.5) ; table.append( 6.5) ; # 1846
table.append(2396027.5) ; table.append( 6.8) ; # 1848
table.append(2396758.5) ; table.append( 7.1) ; # 1850
table.append(2397488.5) ; table.append( 7.3) ; # 1852
table.append(2398219.5) ; table.append( 7.5) ; # 1854
table.append(2398949.5) ; table.append( 7.7) ; # 1856
table.append(2399680.5) ; table.append( 7.8) ; # 1858
table.append(2400410.5) ; table.append( 7.9) ; # 1860
table.append(2401141.5) ; table.append( 7.5) ; # 1862
table.append(2401871.5) ; table.append( 6.4) ; # 1864
table.append(2402602.5) ; table.append( 5.4) ; # 1866
table.append(2403332.5) ; table.append( 2.9) ; # 1868
table.append(2404063.5) ; table.append( 1.6) ; # 1870
table.append(2404793.5) ; table.append( -1.0) ; # 1872
table.append(2405524.5) ; table.append( -2.7) ; # 1874
table.append(2406254.5) ; table.append( -3.6) ; # 1876
table.append(2406985.5) ; table.append( -4.7) ; # 1878
table.append(2407715.5) ; table.append( -5.4) ; # 1880
table.append(2408446.5) ; table.append( -5.2) ; # 1882
table.append(2409176.5) ; table.append( -5.5) ; # 1884
table.append(2409907.5) ; table.append( -5.6) ; # 1886
table.append(2410637.5) ; table.append( -5.8) ; # 1888
table.append(2411368.5) ; table.append( -5.9) ; # 1890
table.append(2412098.5) ; table.append( -6.2) ; # 1892
table.append(2412829.5) ; table.append( -6.4) ; # 1894
table.append(2413559.5) ; table.append( -6.1) ; # 1896
table.append(2414290.5) ; table.append( -4.7) ; # 1898
table.append(2415020.5) ; table.append( -2.7) ; # 1900
table.append(2415750.5) ; table.append( -0.0) ; # 1902
table.append(2416480.5) ; table.append( +2.6) ; # 1904
table.append(2417211.5) ; table.append( 5.4) ; # 1906
table.append(2417941.5) ; table.append( 7.7) ; # 1908
table.append(2418672.5) ; table.append( 10.5) ; # 1910
table.append(2419402.5) ; table.append( 13.4) ; # 1912
table.append(2420133.5) ; table.append( 16.0) ; # 1914
table.append(2420863.5) ; table.append( 18.2) ; # 1916
table.append(2421594.5) ; table.append( 20.2) ; # 1918
table.append(2422324.5) ; table.append( 21.2) ; # 1920
table.append(2423055.5) ; table.append( 22.4) ; # 1922
table.append(2423785.5) ; table.append( 23.5) ; # 1924
table.append(2424516.5) ; table.append( 23.9) ; # 1926
table.append(2425246.5) ; table.append( 24.3) ; # 1928
table.append(2425977.5) ; table.append( 24.0) ; # 1930
table.append(2426707.5) ; table.append( 23.9) ; # 1932
table.append(2427438.5) ; table.append( 23.9) ; # 1934
table.append(2428168.5) ; table.append( 23.7) ; # 1936
table.append(2428899.5) ; table.append( 24.0) ; # 1938
table.append(2429629.5) ; table.append( 24.3) ; # 1940
table.append(2430360.5) ; table.append( 25.3) ; # 1942
table.append(2431090.5) ; table.append( 26.2) ; # 1944
table.append(2431821.5) ; table.append( 27.3) ; # 1946
table.append(2432551.5) ; table.append( 28.2) ; # 1948
table.append(2433282.5) ; table.append( 29.1) ; # 1950
table.append(2434012.5) ; table.append( 30.0) ; # 1952
table.append(2434743.5) ; table.append( 30.7) ; # 1954
table.append(2435473.5) ; table.append( 31.4) ; # 1956
table.append(2436204.5) ; table.append( 32.2) ; # 1958
table.append(2436934.5) ; table.append( 33.1) ; # 1960
table.append(2437665.5) ; table.append( 34.0) ; # 1962
table.append(2438395.5) ; table.append( 35.0) ; # 1964
table.append(2439126.5) ; table.append( 36.5) ; # 1966
table.append(2439856.5) ; table.append( 38.3) ; # 1968
table.append(2440587.5) ; table.append( 40.2) ; # 1970
table.append(2441317.5) ; table.append( 42.184) ; # 1972 JAN 1
table.append(2441499.5) ; table.append( 43.184) ; # 1972 JUL 1
table.append(2441683.5) ; table.append( 44.184) ; # 1973 JAN 1
table.append(2442048.5) ; table.append( 45.184) ; # 1974 JAN 1
table.append(2442413.5) ; table.append( 46.184) ; # 1975 JAN 1
table.append(2442778.5) ; table.append( 47.184) ; # 1976 JAN 1
table.append(2443144.5) ; table.append( 48.184) ; # 1977 JAN 1
table.append(2443509.5) ; table.append( 49.184) ; # 1978 JAN 1
table.append(2443874.5) ; table.append( 50.184) ; # 1979 JAN 1
table.append(2444239.5) ; table.append( 51.184) ; # 1980 JAN 1
table.append(2444786.5) ; table.append( 52.184) ; # 1981 JUL 1
table.append(2445151.5) ; table.append( 53.184) ; # 1982 JUL 1
table.append(2445516.5) ; table.append( 54.184) ; # 1983 JUL 1
table.append(2446247.5) ; table.append( 55.184) ; # 1985 JUL 1
table.append(2447161.5) ; table.append( 56.184) ; # 1988 JAN 1
table.append(2447892.5) ; table.append( 57.184) ; # 1990 JAN 1
table.append(2448257.5) ; table.append( 58.184) ; # 1991 JAN 1
table.append(2448804.5) ; table.append( 59.184) ; # 1992 JUL 1
table.append(2449169.5) ; table.append( 60.184) ; # 1993 JUL 1
table.append(2449534.5) ; table.append( 61.184) ; # 1994 JUL 1
table.append(2450083.5) ; table.append( 62.184) ; # 1996 JAN 1
table.append(2450630.5) ; table.append( 63.184) ; # 1997 JUL 1
table.append(2451179.5) ; table.append( 64.184) ; # 1999 JAN 1
table.append(2453736.5) ; table.append( 65.184) ; # 2006 JAN 1
table.append(2454832.5) ; table.append( 66.184) ; # 2009 JAN 1
table.append(2456109.5) ; table.append( 67.184) ; # 2012 JUL 1 UTC-TAI = - 35s
table.append(2457204.5) ; table.append( 68.184) ; # 2015 JUL 1 UTC-TAI = - 36s
table.append(2457754.5) ; table.append( 69.184) ; # 2017 JAN 1 UTC-TAI = - 37s
jjmax=table[(indexmax-1)*2]
if (jj<=2067314.5):
# --- date <= anne948 ---
t=(jj- 2451545.0)/36525.
ds=2715.6+573.36*t+46.5*t*t
dt=ds;
return dt
if (jj<=2312752.5):
#--- date <= anne1620 ---
t=(jj- 2451545.0)/36525.
ds=50.6+67.5*t+22.5*t*t
dt=ds
return dt
if (jj<=jjmax):
# --- date <= indexmax ---
k2=indexmax
for k in range(1,indexmax):
k2=k
jj2=table[k2*2]
if (jj<=jj2):
break
jj2=table[k2*2]
ds2=table[k2*2+1]
jj1=table[(k2-1)*2]
ds1=table[(k2-1)*2+1]
ds=ds1+(jj-jj1)/(jj2-jj1)*(ds2-ds1)
dt=ds
return dt
# --- extrapolation ---
ds=table[2*indexmax+1]
dt=ds
return dt
# ========================================================
# === internal methods : Physical data for planets
# ========================================================
def _mc_magplanet(self, r, delta, planete_num, phase, l, b):
"""
/***************************************************************************/
/* Calcule la magnitude d'une planete */
/***************************************************************************/
/* Meuus page 269 et 302-303 pour Saturne */
/***************************************************************************/
"""
i=phase/(self._DR)
mag=0.
diamapp=0.
if (r==0.):
r=1e-10
if (delta==0.):
delta=1e-10
if (r*delta>0):
mag=5*log10(r*delta)
if (planete_num == self._PLN_SOLEIL):
mag=2.5*log10(delta)
if (planete_num == self._PLN_MERCURE):
mag+=(-0.42+.0380*i-.000273*i*i+.000002*i*i*i)
diamapp=2.*atan(1.6303e-5/delta)
elif (planete_num == self._PLN_VENUS):
mag+=(-4.40+.0009*i+.000239*i*i-.00000065*i*i*i)
diamapp=2.*atan(4.0455e-5/delta)
elif (planete_num == self._PLN_MARS):
mag+=(-1.52+.016*i)
diamapp=2.*atan(2.2694e-5/delta);
elif (planete_num == self._PLN_JUPITER):
mag+=(-9.40+.005)
diamapp=2.*atan(4.7741e-4/delta)
elif (planete_num == self._PLN_URANUS):
mag+=(-7.19)
diamapp=2.*atan(1.6979e-4/delta)
elif (planete_num == self._PLN_NEPTUNE):
mag+=(-6.87)
diamapp=2.*atan(1.6243e-4/delta)
elif (planete_num == self._PLN_PLUTON):
mag+=(-1.00)
diamapp=2.*atan(1.5608e-5/delta)
elif (planete_num == self._PLN_SOLEIL):
mag+=(-26.86)
diamapp=2.*atan(4.6541e-3/delta)
elif ((planete_num == self._PLN_LUNE) or (planete_num == self._PLN_LUNE_ELP)):
mag+=(0.38+2.97*(i/100.)-0.78*(i/100.)*(i/100.)+.90*(i/100.)*(i/100.)*(i/100.))
diamapp=2.*atan(1.1617e-5/delta)
elif (planete_num == self._PLN_SATURNE):
i=28.08*(self._DR)
o=169.51*(self._DR)
n=113.67*(self._DR)
bb=asin(sin(i)*cos(b)*sin(l-o)-cos(i)*sin(b))
lp=l-(.01759/r)*(self._DR)
bp=b-(.000764*cos(l-n)/r)*(self._DR)
u1=atan2(sin(i)*sin(bp)+cos(i)*cos(bp)*sin(lp-o),cos(bp)*cos(lp-o))
u2=atan2(sin(i)*sin(b)+cos(i)*cos(b)*sin(l-o),cos(b)*cos(l-o))
bb=sin(fabs(bb))
mag+=(-8.68+.044*fabs(u1-u2)/(self._DR)-2.60*bb+1.25*bb*bb)
diamapp=2.*atan(4.0395e-4/delta)
return mag, diamapp
def _mc_libration(self, jj, planete, longmpc, rhocosphip, rhosinphip, equinoxe, astrometric, asd, dec, delta, r, ls, bs):
"""
/***************************************************************************/
/* Calcul de la libration apparentes de la Lune a jj donne. */
/***************************************************************************/
/* D'apres Meeus "Astronomical Algorithms" p341-347 */
/* lonc : longitude selenographique du centre topocentrique */
/* latc : latitude selenographique du centre topocentrique */
/* p : position de l'angle de l'axe de rotation (axe des poles) */
/* lons : longitude selenographique du centre heliocentrique */
/* lats : latitude selenographique du centre heliocentrique */
/***************************************************************************/
inputs asd, dec a l'equinoxe de la date
"""
# --- nutation ---*/
eps = self._mc_obliqmoy(jj)
dpsi, deps = self._mc_nutation(jj,1)
eps += deps
# --- Moon arguments ---*/
T=(jj-2451545.)/36525.;
d=(297.8502042+445267.1115168*T-.00016300*T*T+T*T*T/545868-T*T*T*T/113065000)*(self._DR);
m=(357.5291092+35999.0502909*T-.0001536*T*T+T*T*T*T/24490000)*(self._DR);
mp=(134.9634114+477198.8676313*T+.0089970*T*T+T*T*T/69699.-T*T*T*T/14712000)*(self._DR);
f=(93.2720993+483202.0175273*T-.0034029*T*T+T*T*T/3526000+T*T*T*T/863310000)*(self._DR);
e=1-.002516*T-.0000074*T*T;
f=(93.2720993+483202.0175273*T-.0034029*T*T+T*T*T/3526000+T*T*T*T/863310000)*(self._DR);
o=(125.0445550-1934.1361849*T+0.0020762*T*T+T*T*T/467410.-T*T*T*T/60616000.)*(self._DR);
i=1.54242*(self._DR);
f=fmod(4*self._PI+f,2*self._PI);
# --- parameters of physical libration ---*/
k1=(119.75+131.849*T)*(self._DR);
k2=(72.56+20.186*T)*(self._DR);
rho=-0.02752*cos(mp)-0.02245*sin(f)+0.00684*cos(mp-2*f)-0.00293*cos(2*f)-0.00085*cos(2*f-2*d);
rho=rho-0.00054*cos(mp-2*d)-0.00020*cos(mp+f)-0.00020*cos(mp+2*f)-0.00020*cos(mp-f);
rho=rho+0.00014*cos(mp+2*f-2*d);
sigma=-0.02816*sin(mp)+0.02244*cos(f)-0.00682*sin(mp-2*f)-0.00279*sin(2*f)-0.00083*sin(2*f-2*d);
sigma=sigma+0.00069*sin(mp-2*d)+0.00040*cos(mp+f)-0.00025*sin(2*mp);
sigma=sigma-0.00023*sin(mp+2*f)+0.00020*cos(mp-f)+0.00019*sin(mp-f)+0.00013*sin(mp+2*f-2*d);
sigma=sigma-0.00010*cos(mp-3*f);
tau=0.02520*e*sin(m)+0.00473*sin(2*mp-2*f)-0.00467*sin(mp)+0.00396*sin(k1)+0.00276*sin(2*mp-2*d);
tau=tau+0.00196*sin(o)-0.00183*cos(mp-f)+0.00115*sin(mp-2*d)-0.00096*sin(mp-d);
tau=tau+0.00046*sin(2*f-2*d)-0.00039*sin(mp-f)-0.00032*sin(mp-m-d)+0.00027*sin(2*mp-m-2*d);
tau=tau+0.00023*sin(k2)-0.00014*sin(2*d)+0.00014*cos(2*mp-2*f)-0.00012*sin(mp-2*f);
tau=tau-0.00012*sin(2*mp)+0.00011*sin(2*mp-2*m-2*d);
tau*=(self._DR);
sigma*=(self._DR);
rho*=(self._DR);
# =============== earth lon,lat =====================*/
# --- coord ecl */
bb = Angle(dec/self._PI)
ll = Angle(asd/self._PI)
coords = Coords((1,ll,bb))
xeq,yeq,zeq = coords.cart()
xec,yec,zec = self._mc_xyzeq2ec(xeq,yeq,zeq,eps)
coords = Coords((xec,yec,zec))
rr,l,b = coords.sphe("rad","rad")
# --- optical topocentric libration ---*/
w=l-dpsi-o;
a=atan2(sin(w)*cos(b)*cos(i)-sin(b)*sin(i),cos(w)*cos(b));
lp=a-f;
bp=asin(-sin(w)*cos(b)*sin(i)-sin(b)*cos(i));
# --- physical topocentric libration ---*/
lpp=-tau+(rho*cos(a)+sigma*sin(a))*tan(bp);
bpp=sigma*cos(a)-rho*sin(a);
# --- topocentric libration ---*/
v=fmod(4*self._PI+lp+lpp,2*self._PI);
if (v>(self._PI)):
v-=(2*(self._PI))
lonc=v;
latc=(bp+bpp);
# --- pole axis position ---*/
v=o+deps+sigma/sin(i);
x=sin(i+rho)*sin(v);
y=sin(i+rho)*cos(v)*cos(eps)-cos(i+rho)*sin(eps);
w=atan2(x,y);
v=asin(sqrt(x*x+y*y)*cos(asd-w)/cos(bp+bpp));
p=fmod(4*self._PI+v,2*self._PI);
# =============== earth lon,lat =====================*/
# --- sun ---*/
l=ls+(self._PI)+delta/r*cos(b)*sin(ls-l);
b=delta/r*b;
# --- optical topocentric libration ---*/
w=l-dpsi-o;
a=atan2(sin(w)*cos(b)*cos(i)-sin(b)*sin(i),cos(w)*cos(b));
lp=a-f;
bp=asin(-sin(w)*cos(b)*sin(i)-sin(b)*cos(i));
# --- physical topocentric libration ---*/
lpp=-tau+(rho*cos(a)+sigma*sin(a))*tan(bp);
bpp=sigma*cos(a)-rho*sin(a);
# --- topocentric libration ---*/
v=fmod(4*self._PI+lp+lpp,2*self._PI);
if (v>(self._PI)):
v-=(2*(self._PI))
lons=v;
lats=(bp+bpp)
return lonc, latc, p, lons, lats
def _mc_physephem_moon(self, jj,planete,xg,yg,zg,longmpc,rhocosphip,rhosinphip, equinoxe, astrometric, asd, dec, delta, r, ls, bs):
"""
/***************************************************************************/
/* Cacul des parametres d'observation physique des planetes */
/***************************************************************************/
/* D'apres "Astronomy with computer". */
/* pages 135-154 (chapitre 7) */
/***************************************************************************/
"""
if ((planete == self._PLN_LUNE) or (planete == self._PLN_LUNE_ELP)):
diamapp_equ=0.;
diamapp_pol=0.;
long1=0.;
long2=0.;
long3=0.;
lati=0.;
posangle_sun=0.;
posangle_north=0.;
# ---
long1,lati,posangle_north,lons,lats = self._mc_libration(jj,planete,longmpc,rhocosphip,rhosinphip,equinoxe, astrometric, asd, dec)
long2=long1;
long3=long1;
r=sqrt(xg*xg+yg*yg+zg*zg);
diamapp_equ=2.*atan(1.1617e-5/r);
diamapp_pol=diamapp_equ;
posangle_sun=0.;
lati_sun=lats;
long1_sun=lons;
return diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun
def _mc_physephem_nomoon(self, jj,planete,xg,yg,zg,x,y,z):
"""
/***************************************************************************/
/* Cacul des parametres d'observation physique des planetes */
/***************************************************************************/
/* D'apres "Astronomy with computer". */
/* pages 135-154 (chapitre 7) */
/***************************************************************************/
"""
if ((planete == self._PLN_LUNE) or (planete == self._PLN_LUNE_ELP)):
return
a0=0.;d0=0.;w1=0.;w2=0.;w3=0.;f=0.;req=0.;
sense=1.;
diamapp_equ=0.;
diamapp_pol=0.;
long1=0.;
long2=0.;
long3=0.;
lati=0.;
posangle_sun=0.;
posangle_north=0.;
t=(jj-2451545.)/36525.;
d=(jj-2451545.);
if (planete == self._PLN_MERCURE):
a0=281.01-0.033*t;
d0=61.45-0.005*t;
req=2439.;
f=0.;
w1=w2=w3=(329.68+6.1385025*d);
sense=1.;
elif (planete == self._PLN_VENUS):
a0=272.76;
d0=67.16;
req=6051.;
f=0.;
w1=w2=w3=(160.20-1.4813688*d);
sense=-1.;
elif (planete == self._PLN_TERRE):
a0=0.-0.641*t;
d0=90.0-0.557*t;
req=6378.14;
f=0.00335281;
w1=w2=w3=(190.16+360.9856235*d);
sense=1.;
elif (planete == self._PLN_MARS):
a0=317.681-0.108*t;
d0=52.886-0.061*t;
req=3393.4;
f=0.0051865;
w1=w2=w3=(176.901+350.8919830*d);
sense=1.;
elif (planete == self._PLN_JUPITER):
a0=268.05-0.009*t;
d0=64.49+0.003*t;
req=71398.;
f=0.0648088;
w1=(67.1+877.900*d);
w2=(43.3+870.270*d);
w3=(284.695+870.536*d);
sense=1.;
elif (planete == self._PLN_SATURNE):
a0=40.589-0.036*t;
d0=83.537-0.004*t;
req=60000.;
f=0.1076209;
w1=w2=(227.2037+844.300*d);
w3=(38.90+810.7939024*d);
sense=1.;
elif (planete == self._PLN_URANUS):
a0=257.311;
d0=-15.175;
req=25400.;
f=0.030;
w1=w2=w3=(203.81-501.1600928*d);
sense=-1.;
elif (planete == self._PLN_NEPTUNE):
n=(357.85+52.316*t)*(self._DR);
a0=299.36+0.7*sin(n);
d0=43.46-0.51*cos(n);
req=24300.;
f=0.0259;
w1=w2=w3=(253.18+536.3128492*d-0.48*sin(n));
sense=1.;
elif (planete == self._PLN_PLUTON):
a0=313.02;
d0=9.09;
req=1500.;
f=0.;
w1=w2=w3=(236.77-56.3623195*d);
sense=-1.;
elif (planete == self._PLN_SOLEIL):
a0=286.13;
d0=63.87;
req=696000.;
f=0.;
w1=w2=w3=(84.182+14.1844000*d);
sense=-1.;
a0*=(self._DR);
d0*=(self._DR);
w1*=(self._DR);
w2*=(self._DR);
w3*=(self._DR);
req/=(self._UA*1e-3);
sina=sin(a0);
cosa=cos(a0);
sind=sin(d0);
cosd=cos(d0);
dx=cosa*cosd;
dy=sina*cosd;
dz=sind;
rho=sqrt(xg*xg+yg*yg);
r=sqrt(xg*xg+yg*yg+zg*zg);
diamapp_equ=2*asin(req/r);
costh=((-xg*zg)*dx+(-yg*zg)*dy+(xg*xg+yg*yg)*dz)/(r*rho);
sinth=(-yg*dx+xg*dy)/rho;
th=atan2(sinth,costh);
posangle_north=fmod(4*self._PI+th,2*self._PI);
# /* w1 */
cosw=cos(w1);
sinw=sin(w1);
e1x=-cosw*sina-sinw*sind*cosa;
e1y= cosw*cosa-sinw*sind*sina;
e1z= sinw*cosd;
e2x= sinw*sina-cosw*sind*cosa;
e2y=-sinw*cosa-cosw*sind*sina;
e2z= cosw*cosd;
e3x= cosd*cosa;
e3y= cosd*sina;
e3z= sind;
sx=-(e1x*xg+e1y*yg+e1z*zg);
sy=-(e2x*xg+e2y*yg+e2z*zg);
sz=-(e3x*xg+e3y*yg+e3z*zg);
phip=atan2(sz,sqrt(sx*sx+sy*sy));
phi=atan2(tan(phip),((1-f)*(1-f)));
lambdp=atan2(sy,sx);
lambd=-1.*sense*lambdp;
lati=phi;
diamapp_pol=(diamapp_equ)*(1-f*cos(phip)*cos(phip));
long1=fmod(4*self._PI+lambd,2*self._PI);
# /* w1 sun */
sx=-(e1x*x+e1y*y+e1z*z);
sy=-(e2x*x+e2y*y+e2z*z);
sz=-(e3x*x+e3y*y+e3z*z);
phip=atan2(sz,sqrt(sx*sx+sy*sy));
phi=atan2(tan(phip),((1-f)*(1-f)));
lambdp=atan2(sy,sx);
lambd=-1.*sense*lambdp;
lati_sun=phi;
long1_sun=fmod(4*self._PI+lambd,2*self._PI);
# /* w2 */
cosw=cos(w2);
sinw=sin(w2);
e1x=-cosw*sina-sinw*sind*cosa;
e1y= cosw*cosa-sinw*sind*sina;
e1z= sinw*cosd;
e2x= sinw*sina-cosw*sind*cosa;
e2y=-sinw*cosa-cosw*sind*sina;
e2z= cosw*cosd;
e3x= cosd*cosa;
e3y= cosd*sina;
e3z= sind;
sx=-(e1x*xg+e1y*yg+e1z*zg);
sy=-(e2x*xg+e2y*yg+e2z*zg);
lambdp=atan2(sy,sx);
lambd=-1.*sense*lambdp;
long2=fmod(4*self._PI+lambd,2*self._PI);
# /* w3 */
cosw=cos(w3);
sinw=sin(w3);
e1x=-cosw*sina-sinw*sind*cosa;
e1y= cosw*cosa-sinw*sind*sina;
e1z= sinw*cosd;
e2x= sinw*sina-cosw*sind*cosa;
e2y=-sinw*cosa-cosw*sind*sina;
e2z= cosw*cosd;
e3x= cosd*cosa;
e3y= cosd*sina;
e3z= sind;
sx=-(e1x*xg+e1y*yg+e1z*zg);
sy=-(e2x*xg+e2y*yg+e2z*zg);
sz=-(e3x*xg+e3y*yg+e3z*zg);
lambdp=atan2(sy,sx);
lambd=-1.*sense*lambdp;
long3=fmod(4*self._PI+lambd,2*self._PI);
# /* sun */
d=sqrt(x*x+y*y+z*z);
dx=-x;
dy=-y;
dz=-z;
costh=((-xg*zg)*dx+(-yg*zg)*dy+(xg*xg+yg*yg)*dz)/(r*rho);
sinth=(-yg*dx+xg*dy)/rho;
th=atan2(sinth,costh);
posangle_sun=fmod(4*self._PI+th,2*self._PI);
if (planete == self._PLN_SOLEIL):
posangle_sun=0.;
lati_sun=lati;
long1_sun=long1;
return diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun
# ========================================================
# === internal methods : Misc
# ========================================================
def _mc_altitude2tp(self, altitude_m, pressure_0m_pascal=101325):
"""
/****************************************************************************/
/* Copute the ICAO temperature & pressure from a given altitude */
/****************************************************************************/
/****************************************************************************/
"""
alti=altitude_m
p0m=pressure_0m_pascal
tk0m=273.15+15;
if (alti<11000):
tk=tk0m-0.0065*alti
p=p0m*pow(tk/tk0m,5.255)
elif (alti<15000):
tk0m=273.15+15
tk=tk0m-0.0065*11000
p=p0m*pow((tk0m-0.0065*alti)/tk0m,5.255)
elif ((alti>=15000) and (alti<20000)):
h1=15000; p1=p0m*pow((tk0m-0.0065*h1)/tk0m,5.255); t1=tk0m-0.0065*11000
h2=20000; p2=5500; t2=273.15-46
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=20000) and (alti<30000)):
h1=20000; p1=5500; t1=273.15-46
h2=30000; p2=1100; t2=273.15-38
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=30000) and (alti<40000)):
h1=30000; p1=1100; t1=273.15-38
h2=40000; p2=300; t2=273.15-5
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=40000) and (alti<50000)):
h1=40000; p1=300; t1=273.15-5
h2=50000; p2=90; t2=273.15+1
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=50000) and (alti<60000)):
h1=50000; p1=90; t1=273.15+1
h2=60000; p2=25; t2=273.15-20
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=60000) and (alti<100000)):
h1=60000; p1=25; t1=273.15-20
h2=100000; p2=0.04; t2=273.15-64
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=100000) and (alti<200000)):
h1=100000; p1=0.04; t1=273.15-64
h2=200000; p2=1.3e-4; t2=273.15-82.2
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=200000) and (alti<400000)):
h1=200000; p1=1.3e-4; t1=273.15-82.2
h2=400000; p2=4.4e-6; t2=273.15-97.3
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
elif ((alti>=200000) and (alti<400000)):
h1=400000; p1=4.4e-6; t1=273.15-97.3
h2=500000; p2=0; t2=273.15-97.7
frac=(alti-h1)/(h2-h1)
p=p1+frac*(p2-p1)
tk=t1+frac*(t2-t1)
else:
p=0
tk=273.15-97.7
return tk,p
# ========================================================
# === internal methods : Ephemeris
# ========================================================
def _mc_radec2app(self, date_utc, longmpc, rhocosphip, rhosinphip, ra, dec, equinox, epoch, mura, mudec, plx, tk,ppa,hump,lnm,outputs):
# --- date
jjutc = date_utc.jd()
jds = [ jjutc ]
# --- home
longmpc *= self._DR
latitude, altitude = self._mc_rhophi2latalt(rhosinphip, rhocosphip)
latrad = latitude * self._DR
# --- outputs (list of keywords, such as RA, DEC, ELONG)
outputs = [ output.upper() for output in outputs ]
results = []
#print("-1. ra={} dec={}".format(ra,dec))
ra*=self._DR ;
dec*=self._DR;
#print("0. ra={} dec={}".format(ra/self._DR,dec/self._DR))
muramas=mura;
mudecmas=mudec;
for jd in jds:
#print("1. ra={} dec={}".format(ra/self._DR,dec/self._DR))
cosdec=cos(dec);
mura=muramas*1e-3/86400/cosdec;
mudec=mudecmas*1e-3/86400;
parallax=plx;
# --- aberration annuelle ---*/
asd2, dec2 = self._mc_aberration_annuelle(jd,ra,dec,1)
ra=asd2
dec=dec2
#print("2. ra={} dec={}".format(ra/self._DR,dec/self._DR))
#--- calcul de mouvement propre ---*/
ra+=(jd-epoch)/365.25*mura;
dec+=(jd-epoch)/365.25*mudec;
#print("3. ra={} dec={}".format(ra/self._DR,dec/self._DR))
# --- calcul de la precession ---*/
asd2, dec2 = self._mc_precad(equinox,ra,dec,jd)
ra=asd2;
dec=dec2;
#print("4. ra={} dec={}".format(ra/self._DR,dec/self._DR))
# --- correction de parallaxe stellaire*/
if (parallax>0):
asd2, dec2 = self._mc_parallaxe_stellaire(jd,ra,dec,parallax)
ra=asd2;
dec=dec2;
#print("5. ra={} dec={}".format(ra/self._DR,dec/self._DR))
# --- correction de nutation */
asd2, dec2 = self._mc_nutradec(jd,ra,dec,1);
ra=asd2
dec=dec2
#print("6. ra={} dec={}".format(ra/self._DR,dec/self._DR))
# --- aberration de l'aberration diurne*/
asd2, dec2 = self._mc_aberration_diurne(jd,ra,dec,longmpc,rhocosphip,rhosinphip,1)
ra=asd2;
dec=dec2;
#print("7. ra={} dec={}".format(ra/self._DR,dec/self._DR))
# --- coordonnees horizontales---*/
ha = self._mc_ad2hd(jd,longmpc,ra)
az, h = self._mc_hd2ah(ha,dec,latrad)
# --- refraction ---*/
#lnm=590,
#refraction = self._mc_refraction2(h,1,tk,ppa,lnm,hump,latitude,altitude)
#h += refraction;
#ha, dec = self._mc_ah2hd(az,h,latrad)
parallactic = self._mc_hd2parallactic(ha,dec,latrad)
ra = self._mc_hd2ad(jd,longmpc,ha)
if 'RA_APP' in outputs or 'RA' in outputs:
results.append(['RA_APP',ra/self._DR])
if 'DEC_APP' in outputs or 'DEC' in outputs :
results.append(['DEC_APP',dec/self._DR])
if 'HA' in outputs:
results.append(['HA',ha/self._DR])
if 'AZ' in outputs:
results.append(['AZ',az/self._DR])
if 'ELEV' in outputs or 'ALT' in outputs or 'ALTITUDE' in outputs:
results.append(['ELEV',h/self._DR])
if 'PARALLACTIC' in outputs:
results.append(['PARALLACTIC',parallactic/self._DR])
return dict(results)
def _mc_ephem_planet(self, planet_name, date_utc, longmpc, rhocosphip, rhosinphip, astrometric, jd_equinox, outputs):
""" Compute a planet ephemeris.
:param planet_name: The name of a planet.
:type planet_name: string
:param date_utc: The Date (UTC)
:type date_utc: Date
:param coords_home: The Coords of the home place (for parallax)
:type coords_home: Coords
:param astrometric: Flag to compute astrometric (=1) or true coordinates (=0)
:type astrometric: int
:param jd_equinox: The equinox julian day if astrometric flag = 1
:type jd_equinox: float
"""
# --- planet
planet_num = self._mc_name2planetnum(planet_name)
# --- The case of a false planet
if (planet_num == self._PLN_OTHERPLANET):
raise Exception
return ""
# --- date
jjutc = date_utc.jd()
jj = self._mc_tu2td(jjutc)
jjd=jj
# --- home
#altitude, longmpc , latitude = coords_home.sphe("rad","deg")
latitude, altitude = self._mc_rhophi2latalt(rhosinphip, rhocosphip)
longmpc*=self._DR
# --- computation method
afm = self._AFM_VFP79
# --- outputs (list of keywords, such as RA, DEC, ELONG)
outputs = [ output.upper() for output in outputs ]
results = []
if (planet_num == self._PLN_LUNE) or (planet_num == self._PLN_LUNE_ELP):
# =====================================
# === The moon only
# =====================================
ttminusutc=jj-jjutc
# --- soleil ---
delta = 1. ; # correction de l'abberation planetaire a priori pour la Lune
jjds = self._mc_aberpla(jjd,delta)
llp, mmp, uup = self._mc_jd2lbr1a(jjds);
coords = self._mc_jd2lbr1b(jjds, self._PLN_SOLEIL, llp, mmp, uup, afm)
xs,ys,zs = coords.cart()
rs,ls,bs = coords.sphe("rad","rad")
# --- correction de la parallaxe ---
eps = self._mc_obliqmoy(jjds)
xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
dxeq,dyeq,dzeq = self._mc_paraldxyzeq(jjds-ttminusutc,longmpc,rhocosphip,rhosinphip)
xs-=dxeq
ys-=dyeq
zs-=dzeq
xs,ys,zs = self._mc_xyzeq2ec(xs,ys,zs,eps)
# --- soleil : coordonnes asd,dec du Soleil ---
if (astrometric==0):
coords = Coords((xs,ys,zs))
rs,ls,bs = coords.sphe("rad","rad")
dpsi, deps = self._mc_nutation(jjd,1)
ls += dpsi
eps += deps
bb = Angle(bs/self._DR)
ll = Angle(ls/self._DR)
coords = Coords((rs,ll,bb))
xs,ys,zs = coords.cart()
xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
coords = Coords((xs,ys,zs))
delta,asds,decs = coords.sphe("rad","rad")
# a ce niveau on retient (asds,decs) pour le calcul de la phase */
# On recommence tout le calcul sans abberation pour la position */
#dans le repere heliocentrique */
# --- Terre ---
llp, mmp, uup = self._mc_jd2lbr1a(jjd);
coords = self._mc_jd2lbr1b(jjd, self._PLN_SOLEIL, llp, mmp, uup, afm)
xs,ys,zs = coords.cart()
r,ls,bs = coords.sphe("rad","rad")
# --- Terre : correction de la parallaxe ---
eps = self._mc_obliqmoy(jjd)
xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
dxeq,dyeq,dzeq = self._mc_paraldxyzeq(jjd-ttminusutc,longmpc,rhocosphip,rhosinphip)
xs-=dxeq
ys-=dyeq
zs-=dzeq
xs,ys,zs = self._mc_xyzeq2ec(xs,ys,zs,eps)
# --- Terre : coordonnes asd,dec du Soleil ---
if (astrometric==0):
coords = Coords((xs,ys,zs))
rs,ls,bs = coords.sphe("rad","rad")
dpsi, deps = self._mc_nutation(jjd,1)
ls += dpsi
eps += deps
bb = Angle(bs/self._DR)
ll = Angle(ls/self._DR)
coords = Coords((r,ll,bb))
xs,ys,zs = coords.cart()
xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
coords = Coords((xs,ys,zs))
delta,asd,dec = coords.sphe("rad","rad")
# ici on ne garde que (xs,ys,zs,delta) pour la position geocentrique
# --- LUNE : planete (coordonnees ecliptiques astrometriques a la date) ---
llp, mmp, uup = self._mc_jd2lbr1a(jjd);
coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
x,y,z = coords.cart()
eps = self._mc_obliqmoy(jjd)
x,y,z = self._mc_xyzec2eq(x,y,z,eps)
coords = Coords((x,y,z))
delta,asd,dec = coords.sphe("rad","rad")
# --- planete corrigee de l'aberration de la lumiere ---
jjd = self._mc_aberpla(jjd,delta)
llp, mmp, uup = self._mc_jd2lbr1a(jjd);
coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
x,y,z = coords.cart()
eps = self._mc_obliqmoy(jjd)
# --- correction de la nutation ---
r,l,b = coords.sphe()
if (astrometric==0):
dpsi, deps = self._mc_nutation(jjd,1)
l += dpsi
eps += deps
bb = Angle(b/self._DR)
ll = Angle(l/self._DR)
coords = Coords((r,ll,bb))
x,y,z = coords.cart()
# --- correction de la parallaxe ---
x,y,z = self._mc_xyzec2eq(x,y,z,eps)
x-=dxeq
y-=dyeq
z-=dzeq
# --- coord. spheriques ---
coords = Coords((x,y,z))
delta,asd,dec = coords.sphe("rad","rad")
#--- parametres physiques --- (not implemented)
if ('APPDIAMEQU' in outputs) or ('APPDIAMPOL' in outputs) or ('LONGI' in outputs) or ('LONGII' in outputs) or ('LONGIII' in outputs) or ('LATI' in outputs) or ('POSNORTH' in outputs) or ('POSSUN' in outputs) or ('LONGI_SUN' in outputs) or ('LATI_SUN' in outputs):
diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun = self._mc_physephem_moon(jjd,planet_num,x,y,z,longmpc, rhocosphip, rhosinphip, asd, dec, delta, r, ls, bs)
# --- parametres elong et magnitude ---
x,y,z = self._mc_he2ge(x,y,z,-xs,-ys,-zs)
r=sqrt(x*x+y*y+z*z)
if ('ELONG' in outputs) or ('PHASE' in outputs) or ('POSSUN' in outputs) or ('MAG' in outputs) or ('APPDIAM' in outputs):
elong, phase, limb = self._mc_elonphaslimb(asd,dec,asds,decs,rs,delta);
posangle_sun = limb
if ('MAG' in outputs) or ('APPDIAM' in outputs):
mag, diamapp = self._mc_magplanet(r, delta, planet_num,phase, l, b)
# --- correction de l'aberration annuelle ---
if (astrometric==0):
asd, dec = self._mc_aberration_annuelle(jjd, asd, dec, 1)
# --- correction de la precession ---
asd,dec = self._mc_precad(jjd,asd,dec,jd_equinox); # equatorial J2000
# --- parametres physiques ---*/
else:
# =====================================
# === All planets and Sun
# =====================================
# --- soleil ---
llp, mmp, uup = self._mc_jd2lbr1a(jjd);
coords = self._mc_jd2lbr1b(jjd, self._PLN_SOLEIL, llp, mmp, uup, afm)
xs,ys,zs = coords.cart()
rs,ls,bs = coords.sphe("rad","rad")
if (planet_num == self._PLN_SOLEIL):
# --- special case of Sun
jjd = self._mc_aberpla(jjd,rs)
llp, mmp, uup = self._mc_jd2lbr1a(jjd);
coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
xs,ys,zs = coords.cart()
rs,ls,bs = coords.sphe("rad","rad")
# --- correction de la parallaxe ---
eps = self._mc_obliqmoy(jjd)
xs,ys,zs = self._mc_xyzec2eq(xs,ys,zs,eps)
dxeq,dyeq,dzeq = self._mc_paraldxyzeq(jjd,longmpc,rhocosphip,rhosinphip)
xs-=dxeq
ys-=dyeq
zs-=dzeq
xs,ys,zs = self._mc_xyzeq2ec(xs,ys,zs,eps)
if (planet_num == self._PLN_SOLEIL):
xg=xs; yg=ys; zg=zs
else:
# --- planet ---
coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
x,y,z = coords.cart()
x,y,z = self._mc_he2ge(x,y,z,xs,ys,zs)
x,y,z = self._mc_xyzec2eq(x,y,z,eps)
coords = Coords((x,y,z))
delta,asd,dec = coords.sphe("rad","rad")
# --- planete corrigee de l'aberration de la lumiere ---
jjd = self._mc_aberpla(jjd,delta)
llp, mmp, uup = self._mc_jd2lbr1a(jjd)
coords = self._mc_jd2lbr1b(jjd, planet_num, llp, mmp, uup, afm)
x,y,z = coords.cart()
rr,l,b = coords.sphe("rad","rad")
xg,yg,zg = self._mc_he2ge(x,y,z,xs,ys,zs)
eps = self._mc_obliqmoy(jjd)
# --- correction de la nutation ---
if (astrometric==0):
coords = Coords((xg,yg,zg))
r,l,b = coords.sphe("rad","rad")
dpsi, deps = self._mc_nutation(jjd,1)
l += dpsi
eps += deps
bb = Angle(b/self._DR)
ll = Angle(l/self._DR)
coords = Coords((r,ll,bb))
xg,yg,zg = coords.cart()
# --- coord. spheriques ---
xg,yg,zg = self._mc_xyzec2eq(xg,yg,zg,eps)
coords = Coords((xg,yg,zg))
delta,asd,dec = coords.sphe("rad","rad")
# --- correction de l'aberration annuelle ---
if (astrometric==0):
asd, dec = self._mc_aberration_annuelle(jjd, asd, dec, 1)
# --- correction de la precession ---
if (planet_num != self._PLN_PLUTON):
asd,dec = self._mc_precad(jjd,asd,dec,jd_equinox); # equatorial J2000
# --- parametres elong et magnitude ---*/
if (planet_num != self._PLN_SOLEIL):
r = rr
if ('ELONG' in outputs) or ('PHASE' in outputs) or ('MAG' in outputs) or ('DIAMAPP' in outputs):
elong, phase = self._mc_elonphas(r, rs, delta)
else:
elong = 0
phase = 0
if ('MAG' in outputs) or ('APPDIAM' in outputs):
mag, diamapp = self._mc_magplanet(r, delta, planet_num, phase, l, b)
# --- parametres physiques ---*/
if ('APPDIAMEQU' in outputs) or ('APPDIAMPOL' in outputs) or ('LONGI' in outputs) or ('LONGII' in outputs) or ('LONGIII' in outputs) or ('LATI' in outputs) or ('POSNORTH' in outputs) or ('POSSUN' in outputs) or ('LONGI_SUN' in outputs) or ('LATI_SUN' in outputs):
x,y,z = self._mc_xyzec2eq(x,y,z,eps) ; # coord. helio equatoriales
diamapp_equ,diamapp_pol,long1,long2,long3,lati,posangle_north,posangle_sun,long1_sun,lati_sun = self._mc_physephem_nommon(jjd,planet_num,xg,yg,zg,x,y,z)
# --- outputs
if 'NAME' in outputs:
results.append(['NAME',planet_name])
if 'RA' in outputs:
results.append(['RA',asd/self._DR])
if 'DEC' in outputs:
results.append(['DEC',dec/self._DR])
if 'EQUINOX' in outputs:
results.append(['EQUINOX',(Date(jd_equinox)).equinox()])
if 'ELONG' in outputs:
results.append(['ELONG',elong/self._DR])
if 'PHASE' in outputs:
results.append(['PHASE',phase/self._DR])
if 'DELTA' in outputs:
results.append(['DELTA',delta])
if 'R' in outputs:
results.append(['R',rs])
if 'MAG' in outputs:
results.append(['MAG',mag])
if 'APPDIAM' in outputs:
results.append(['APPDIAM',diamapp/self._DR])
if 'APPDIAMEQU' in outputs:
results.append(['APPDIAMEQU',diamapp_equ/self._DR])
if 'APPDIAMPOL' in outputs:
results.append(['APPDIAMPOL',diamapp_pol/self._DR])
if 'LONGI' in outputs:
results.append(['LONGI',long1/self._DR])
if 'LONGII' in outputs:
results.append(['LONGII',long2/self._DR])
if 'LONGIII' in outputs:
results.append(['LONGIII',long3/self._DR])
if 'LATI' in outputs:
results.append(['LATI',lati/self._DR])
if 'POSNORTH' in outputs:
results.append(['POSNORTH',posangle_north/self._DR])
if 'POSSUN' in outputs:
results.append(['POSNORTH',posangle_sun/self._DR])
if 'LONGI_SUN' in outputs:
results.append(['LONGI_SUN',long1_sun/self._DR])
if 'LATI_SUN' in outputs:
results.append(['LATI_SUN',lati_sun/self._DR])
return dict(results)
# ========================================================
# === special methods
# ========================================================
def __init__(self):
""" Object initialization where planet_name is the name of the planet
"""
self._init_mechanics()