L'objectif de cette séance est de:
Nous allons utiliser, dans la bibliothèque standard, notamment:
datetime pour la gestion intuitive du temps;Parmi les bibliothèques non officielles mais au statut standard:
numpy pour la gestion des tableaux de données;matplotlib pour l'affichage;requests pour la récupération de données sur Internet;Nouveauté pour cette séance, le module skyfield spécialisé dans les traitements appliqués à l'astronomie.
In [ ]:
from skyfield.api import load, utc
ts = load.timescale()
# chargement des éphémérides
planets = load('de421.bsp')
In [ ]:
earth = planets['earth']
sun = planets['sun']
moon = planets['moon']
Des méthodes peuvent être enchaînées sur les corps célestes présents dans planets.
In [ ]:
# Position de la Terre au 1er janvier 2017
earth.at(ts.utc(2017, 1, 1))
La méthode utc permet d'entrer des données temporelles manuellement, mais aussi d'utiliser le module datetime de Python.
Dans ce cas, il faut spécifier le fuseau horaire!
In [ ]:
import datetime
now = datetime.datetime.now(utc)
now
In [ ]:
# Position relative du soleil par rapport à la Terre au 1er janvier 2017, 12h10
earth.at(ts.utc(now)).observe(sun)
Pour manipuler des objets spatiaux, il convient d'être familier avec différents systèmes de coordonnées. Nous allons manipuler trois systèmes de coordonnées:
(source Wikipedia)
Le système de coordonnées équatoriales est un système de coordonnées célestes dont les valeurs sont indépendantes de la position de l'observateur. Ce système utilise comme plan de référence la projection, sur la sphère céleste, de l'équateur de la Terre. Cette projection est l'équateur céleste, qui divise le ciel en deux hémisphères, chacun ayant comme axe de référence la projection d'un pôle terrestre, perpendiculaire à l'équateur céleste. À partir de ces divisions, le système permet d'établir deux coordonnées angulaires : l'ascension droite et la déclinaison.
In [ ]:
earth.at(ts.utc(now)).observe(sun).radec()
In [ ]:
In [ ]:
# %load solutions/declinaison_nulle.py
In [ ]:
In [ ]:
# %load solutions/tropiques.py
(source Wikipedia)
Le système de coordonnées horizontales, également appelé système local ou système de coordonnées alt-azimutales, est un système de coordonnées célestes utilisé en astronomie par un observateur au sol. Le système sépare le ciel en deux hémisphères : l'un situé au-dessus de l'observateur et l'autre situé au-dessous, caché par le sol. Le grand cercle séparant les deux hémisphères situe le plan horizontal, à partir duquel est établi une altitude et un azimut, qui constituent les deux principales coordonnées de ce système.
Pour l'heure d'été, nous pourrons effectuer les mêmes calculs le 1er juillet.
In [ ]:
from collections import namedtuple
from skyfield.api import Topos
import pytz
city = namedtuple('city', ['coords', 'winter', 'summer'])
fmt = '%H:%M:%S %Z %z'
msg = "Midi à {} ({}): azimuth de {:0.02f} deg"
timezones = {
# ORY
'Europe/Paris': city(coords=(48.725278, 2.359444), winter=1, summer=2),
# SXF
'Europe/Berlin': city(coords=(52.380001, 13.52258), winter=1, summer=2),
# LHR
'Europe/London': city(coords=(51.4775, -0.461389), winter=0, summer=1),
# HND
'Asia/Tokyo': city(coords=(35.552258, 139.779694), winter=9, summer=9),
# EZE
'America/Buenos_Aires': city(coords=(-34.822222, -58.535833), winter=-3, summer=-3)
}
print("Heures d'hiver:")
for tz, city in timezones.items():
noon = datetime.datetime(2017, 1, 1, 12 - city.winter, tzinfo=utc)
_, az, _ = (earth + Topos(*city.coords)).at(ts.utc(noon)).observe(sun).apparent().altaz()
print(msg.format(tz, noon.astimezone(pytz.timezone(tz)).strftime(fmt), az.degrees))
print()
print("Heures d'été:")
for tz, city in timezones.items():
if 'Europe' not in tz: continue
noon = datetime.datetime(2017, 7, 1, 12 - city.summer, tzinfo=utc)
_, az, _ = (earth + Topos(*city.coords)).at(ts.utc(noon)).observe(sun).apparent().altaz()
print(msg.format(tz, noon.astimezone(pytz.timezone(tz)).strftime(fmt), az.degrees))
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [ ]:
today = datetime.date.today()
# Combien y a-t-il de jours ce mois-ci ?
import calendar
_, total_days = calendar.monthrange(today.year, today.month)
# La Lune observée à Paris
paris = earth + Topos(*timezones['Europe/Paris'].coords)
stack=[]
for day in range(1, total_days+1):
dt = ts.utc(datetime.date(today.year, today.month, day))
_, moon_az, _ = paris.at(dt).observe(moon).apparent().altaz()
_, sun_az, _ = paris.at(dt).observe(sun).apparent().altaz()
stack.append(moon_az.radians - sun_az.radians)
# Revenir entre -π et π
stack = np.angle(np.exp(1j*np.array(stack)))
# Détecter le premier passage de -π à π
wh = np.where(np.abs(stack[1:]-stack[:-1]) > np.pi)[0][0]
fig = plt.figure(figsize=(10, 7))
ax = fig.gca()
# Un trait vertical par jour
for i in range(total_days):
ax.plot([i,i], [-np.pi, np.pi], color='#eeeeee')
# Un trait horizontal par phase principale
pi = np.zeros(stack.size)
phase = ['Pleine lune', 'Premier quartier', 'Nouvelle lune',
'Dernier quartier', 'Pleine lune']
for i in range(5):
plt.plot((i-2)*np.pi/2 + pi, '--', color="#aaaaaa")
plt.annotate(phase[i], (5, (i-2)*np.pi/2 + .1), )
# L'angle d'éclairage sur la Lune, vue de la Terre
plt.plot(list(range(wh+1)), stack[:wh+1], color="#f13a31")
plt.plot(list(range(wh+1, total_days)), stack[wh+1:], color="#f13a31")
# Les axes
ax.set_xticks(list(range(total_days)))
ax.set_xticklabels(list(range(1, total_days+1)))
ax.xaxis.set_ticks_position('bottom')
ax.set_yticks([(i-2)*np.pi/2 for i in range(5)])
ax.set_yticklabels(["- π", "- π/2", "0", "π/2", "π"])
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('axes', -0.02))
# Le titre
month_name = [None, "janvier", "février", "mars",
"avril", "mai", "juin","juillet",
"août", "septembre", "octobre",
"novembre", "décembre"]
ax.set_title("Les phases de la Lune en {} {}".format(month_name[today.month], today.year))
ax.set_frame_on(False)
Le système global GCRS (geocentric coordinates referential system) spécifie la position en coordonnées cartésiennes $(x, y, z)$ d'objets proches de la Terre (comme les satellites). Les coordonnées $x$ et $y$ sont projetées sur le plan équatorial, avec $x$ qui point vers le point vernal. La coordonnées $z$ pointe elle vers le pôle Nord.
Il est relativement facile de calculer une latitude, une longitude et une altitude à partir de coordonnées $(x, y, z)$ mais il est important de connaître une heure de référence pour se repérer par rapport à une coordonnée longitudinale terrestre. N'oublions pas que le système GCRS étant indépendant de la rotation de la Terre sur elle-même.
On manipule alors le temps sidéral qui est à un instant et en un lieu donné l'angle horaire du point vernal.
Pour nos calculs, nous allons utiliser l'interface de la bibliothèque geodesy que nous avons utilisée lors de la séance précédente, ainsi que skyfield pour le calcul du temps sidéral. Chaque objet de type Time a un attribut gmst pour Greenwich Mean Sidereal Time, soit l'angle horaire au méridien de Greenwich par rapport au point vernal.
In [ ]:
import geodesy.wgs84 as geo
from skyfield.api import utc
now = datetime.datetime.now(utc)
t = ts.utc(now)
sun_pos = earth.at(t).observe(planets['sun'])
print(sun_pos.position) # by default in astronomic units
print(sun_pos.position.m)
lat, lon, _ = geo.cartesian_to_geodesic(*sun_pos.position.m)
print("lat: {:.4f} deg lon: {:.4f} deg".format(lat, (lon - 15 * t.gmst))) # 15 * 24 h = 360 deg
In [ ]:
In [ ]:
# %load solutions/between_m180_180.py
Le NORAD et la NASA mesurent et calculent la position de satellites, stations spatiales et débris orbitaux. Ils publient leurs observations sous la forme de Two-Line Elements (TLE), qui peuvent être représentés sur le modèle de la ligne suivante (pour l'exemple de la station spatiale internationale):
ISS (ZARYA)
1 25544U 98067A 14273.50403866 .00012237 00000-0 21631-3 0 1790
2 25544 51.6467 297.5710 0002045 126.1182 27.2142 15.50748592907666Le format est décrit sur la page Wikipedia relative aux TLE. Le détail est intéressant mais n'aura pas grande importance pour nous: nous allons utiliser une bibliothèque qui lit ce format pour propager la trajectoire des objets et donner leur position à un moment donné.
Le minimum à savoir est que:
ISS pour International Space Station);1, ici 14273.50403866 donne le $t_0$ pour la propagation. Les deux premiers chiffres (14) font référence à l'année (2014), la suite est le jour de l'année (273.50403866), soit le 273e jour de l'année (30 septembre).Nous allons utiliser un code de propagation de trajectoires qui utilise un modèle numérique nommé SGP4. Celui-ci est précis pendant quelques heures à partir du $t_0$; il est alors important de mettre régulièrement à jour les données depuis cette page.
In [ ]:
# for the example
from skyfield.api import EarthSatellite
iss_text = """
ISS (ZARYA)
1 25544U 98067A 14273.50403866 .00012237 00000-0 21631-3 0 1790
2 25544 51.6467 297.5710 0002045 126.1182 27.2142 15.50748592907666
"""
line1, line2 = iss_text.splitlines()[-2:]
iss = EarthSatellite(line1, line2)
iss
In [ ]:
# compute present position (makes no sense after so many years...)
iss.at(ts.utc(now)).position.m
In [ ]:
In [ ]:
# %load solutions/requests.py
In [ ]:
In [ ]:
# %load solutions/iss_track.py
Le terminateur est la ligne fictive et mobile qui sépare la face éclairée et la face illuminée de la Terre. On trouve sur la page Wikipedia la formule qui donne ses coordonnées ainsi qu'une vidéo qui montre la traversée du terminateur par la station spatiale internationale. Cette formule prend en compte notamment la déclinaison du soleil à l'instant donné (équivalente à l'inclinaison de l'axe de rotation de la planète par rapport à l'écliptique).
La station spatiale internationale peut être visible sur Terre la nuit quand ses panneaux solaires reflètent les rayons solaires vers la Terre. Pour essayer de modéliser ce phénomène, nous allons nous pencher sur une série de définitions du crépuscule.
Le terminateur trace la ligne des points sur Terre qui voient le soleil sur l'horizon (élévation nulle), mais:
Nous considérerons alors que la station spatiale internationale est visible de notre position actuelle si :
Bonus: Ajouter les limites de nuits civile, aéronautique et astronomique.
Ce bonus est intéressant mais n'empêche pas d'aller plus loin. On pourra comparer le résultat avec cette carte.
In [ ]:
In [ ]:
# %load solutions/terminator.py
In [ ]:
In [ ]:
# %load solutions/compute_next_pass.py
In [ ]:
In [ ]:
# %load solutions/map_next_pass.py
La photo suivante a été prise dans la nuit du 13 au 14 janvier 2017 depuis le site l'ESA Madrid.
L'anecdote est racontée ici. On y apprend que le passage de l'ISS devant la Lune a duré 0,56 secondes !
Données :
pour la taille de la Lune dans le ciel, on considère souvent qu'elle occupe un demi-degré angulaire ;
les coordonnées GPS de l'ESA Madrid :
(40.4438, -3.9529)
la dernière TLE publiée le 13 janvier 2017 :
1 25544U 98067A 17013.66453757 .00002774 00000-0 49270-4 0 9991
2 25544 51.6436 88.6266 0007381 79.9762 16.7314 15.54061850 37728
In [ ]:
# Date en question
now = datetime.datetime(2017, 1, 14, tzinfo=utc)
# Coordonnées GPS depuis leur site
esac_madrid = earth + Topos(40.4438, -3.9529)
# Archives TLE
iss = EarthSatellite(
"1 25544U 98067A 17013.66453757 .00002774 00000-0 49270-4 0 9991",
"2 25544 51.6436 88.6266 0007381 79.9762 16.7314 15.54061850 37728"
)
visible = passing_over(
now, esac_madrid,
lambda iss_alt, sun_alt: iss_alt.degrees > 10,
horizon=datetime.timedelta(days=1),
timestep=datetime.timedelta(minutes=1)
)
In [ ]:
# Compute a new track (and keep additional parameters)
MoonPoint = namedtuple(
"MoonPoint", ["iss_alt", "iss_az", "moon_alt", "moon_az", "localtime"]
)
def moon_track(start, position):
track = []
moon = planets["moon"]
for k in range(1200): # 10 minutes à 0.5 secondes
t = ts.utc(start + k * datetime.timedelta(seconds=.5))
iss_alt, iss_az, _ = position.at(t).observe(earth + iss).apparent().altaz()
moon_alt, moon_az, _ = position.at(t).observe(moon).apparent().altaz()
if iss_alt.degrees > 10:
point = MoonPoint(
iss_alt,
iss_az,
moon_alt,
moon_az,
t.astimezone(pytz.timezone("Europe/Madrid")).strftime("%H:%M:%S"),
)
track.append(point)
return track
# Compute the track of the pass over
track = moon_track(visible[0][0] - datetime.timedelta(minutes=3), esac_madrid)
small_track = track[368:372]
fig = plt.figure()
ax = fig.gca()
plt.plot(
[t.iss_az.degrees for t in small_track],
[t.iss_alt.degrees for t in small_track],
"-o",
color="#aaaaaa",
)
for t in small_track:
c = plt.Circle(
(t.moon_az.degrees, t.moon_alt.degrees),
radius=0.25, # pertinence du 0.25 ?
facecolor="#d4cf6a",
edgecolor="#7d7813",
)
ax.add_patch(c)
ax.annotate(t.localtime, (t.iss_az.degrees, t.iss_alt.degrees + 0.1))
ax.axis("scaled")
ax.set_xlim((154, 157.5))
ax.set_ylim((61, 63))
ax.set_frame_on(False)
Le calcul de l'élévation du soleil est donné par la formule suivante :
$$h_{sun} = \arcsin\left(\sin{\varphi} \cdot \sin{\delta_{sun}} + \cos{\varphi}\cdot\cos{\delta_{sun}}\cdot\cos\left(\lambda + \tau\right)\right)$$avec:
Pour le terminateur simple, avec $h_{sun} = 0$:
\begin{eqnarray} \sin{\varphi} \cdot \sin{\delta_{sun}} &=& -\cos{\varphi} \cdot \cos{\delta_{sun}} \cdot \cos\left(\lambda + \tau\right)\\ \tan{\varphi} &=& - \frac{\cos\left(\lambda + \tau\right)}{\tan{\delta_{sun}}} \end{eqnarray}Sinon, on peut réécrire l'équation
$$\sin{h_{sun}} = a \sin{\varphi} + b \cos{\varphi}$$avec $a = \sin{\delta_{sun}}$ et $b = \cos{\delta_{sun}}\cos\left(\lambda + \tau\right)$.
On résout cette équation en posant: $$y = \arctan\frac{b}{a} = \arctan\frac{\cos\left(\lambda + \tau\right)}{\tan\delta_{sun}}$$
Alors:
\begin{eqnarray} \sin{h_{sun}} &=& a\sin{\varphi}+b\cos{\varphi}\\ &=& \sqrt{a^2+b^2}\left(\frac{a}{\sqrt{a^2+b^2}}\sin{\varphi}+\frac{b}{\sqrt{a^2+b^2}}\cos{\varphi}\right)\\ &=& \sqrt{a^2+b^2}\Big(\sin{\varphi}\cos{y}+\cos{\varphi}\sin{y}\Big)\\ &=& \sqrt{a^2+b^2}\sin(\varphi + y) \end{eqnarray}$$ \varphi = - \arctan\left(\dfrac{\cos\left(\lambda + \tau\right)}{\tan{\delta_{sun}}}\right) + \arcsin\left(\frac{\sin{h_{sun}}}{\sqrt{\sin^2 {\delta_{sun}} + \cos^2{\delta_{sun}}\cos^2 \left(\lambda + \tau\right)}}\right)$$
In [ ]: