import numpy as np
import scipy as sc
import matplotlib.pyplot as plt
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)
import datetime
import pandas as pd
import time
#Option d'affichage et taille de police sur les figures:
fs=20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
#plt.rc('text', usetex=True)
plt.rcParams["font.family"] = "Droid Sans Mono"
%matplotlib inline
On se propose dans ce TN, d'utiliser les connaissances acquises lors des 2 premiers TN pour construire les principaux éléments d'un astrolabe. Cet instrument ancien datant de l'antiquité, a été utilisé jusqu'au 17ème siècle pour les calculs d'éphéméride astronomique et également pour connaitre l'heure locale en mesurant la hauteur du soleil dans le ciel ou celle d'une étoile visible la nuit.
Je vous conseille très fortement de visiter le musée des Arts et Métiers ou de nombreux astrolabes d'origine diverses sont exposés (Ci-contre l'Astrolabe d'Arsenius de 1569).
Ce TN doit vous permettre de mettre en oeuvre les différents éléments abordé dans le TN1 et le TN2. L'idée étant de construire succéssivement l'équation du temps rapportée à l'année civile, la projection stéréographique gravée sur les Tympans pour une latitude donnée ainsi que l'Araignée représentant les étoiles les plus brillante de l'hémisphère Nord.
En reprenant les fonction du TN2 représenter dans un diagramme polaire l'équation du temps en minute en fonction des mois de l'année civile.
#Lattitude du point d'observation
Lat=48.8
# Longitude pour la correction horaire
long=2.47
# Calendrier
N=np.arange(1,366)
N1=np.cumsum([1,31,28,31,30,31,30,31,31,30,31,30]);
#Longitude Ecliptique
long_eq=2*np.pi*(N-80)/365
long_eq[long_eq<0]=2*np.pi+long_eq[long_eq<0] # !! Il faut fournir un vecteur de donnée positives pour polar()
long_eq1=2*np.pi*(N1-80)/365
long_eq1[long_eq1<0]=2*np.pi+long_eq1[long_eq1<0]
def Hsol(N,L):
long_ec=2*np.pi*(N+10)/365
delta=23.44*np.cos(long_ec)
return 90-L-delta
def Real_time_eq(L,N):
rad=np.pi/180
# Ecart au 1er Janvier de l'époque J2000:
d0=(datetime.date(2022,1,1)-datetime.date(2000,1,1)).days
N=N+d0
M=357.5291+0.98560028*N
#Excentricité de l'orbite Terrestre
e=0.01671
#Obliquité
eps=23.43929
C=(2*e-0.25*e**3)*np.sin(M*rad)+1.25*e**2*np.sin(2*M*rad)+(13/12)*e**3*np.sin(3*M*rad)
l=280.4665+0.98564736*N+C
y=np.tan(eps*rad/2)
R=-y**2*np.sin(2*l*rad)+0.5*y**4*np.sin(4*l*rad)-(y**6/3)*np.sin(6*l*rad)
T=60*(C+R)
#Hauteur solaire
P=Hsol(N,L)
return 4*T,P
month=['Jan','Fev','Mar','Avr','Mai','Juin','Juil','Aout','Sep','Oct','Nov','Dec']
Lfigure=30#cm
fig=plt.figure(figsize=(Lfigure/2.54,Lfigure/2.54), tight_layout=True, linewidth=2)
ax=fig.add_subplot(111,projection='polar')
ax.set_theta_zero_location('E')#Angle 0 en haut
plt.xticks(long_eq1, month, fontsize=fs)
plt.yticks(np.arange(-18,16,5), fontsize=fs)
plt.plot(...,...)
ax.xaxis.set_minor_locator(AutoMinorLocator())
ax.tick_params(which='both', width=2)
ax.tick_params(which='major', length=7)
ax.tick_params(which='minor', length=4, color='r')
ax.xaxis.grid(True, which='major')
#plt.ylim(-30,30)
La seconde partie est le coeur de l'astrolabe. C'est le Tympan, il indique les graduations des déclinaisons en fonction de l'azimute et donc de l'heure de la journée. Pour représenter ces graduations on utilise la projection stéréographique qui permet de représenter une sphère sur un plan.
Un point dans le ciel observé depuis une latitude $\phi$ ayant des coordonées locales $(hauteur, azimute)=(h,\Theta)$ se projette sur la cercle de rayon $R$ et de centre $(0,y_c)$:
$$ (R,y_c)=\left(\dfrac{\cos(h)}{\sin(\phi)+\sin(h)},\dfrac{\cos(\phi)}{\sin(\phi)+\sin(h)}\right) $$En utilisant les liens fourni, retrouver la formule précédente pour la projection stéréographique
Représenter les cercles de projection de hauteur et d'azimute tous les 5 degrès
# Vecteur des hauteur 0...90
high=np.arange(...)*np.pi/180
# Vecteur des azimutes 5..360
azimut=np.arange(...)*np.pi/180
# Vecteur radial pour la représentation
theta=np.linspace(0,2*np.pi,150)
# Latitude en radians
phi=Lat*np.pi/180
# Fonction tracant le cerle de centre *center* de rayon *radius*
# angle_pos correspond à la description radiale
def circle(ax,center, radius,angle_pos,color='k',lw=1.):
x = radius * np.cos( angle_pos )
y = radius * np.sin( angle_pos )
ax.plot( x-center[0], center[1]-y,color,linewidth=lw )
Lfigure=30#cm
fig=plt.figure(figsize=(Lfigure/2.54,Lfigure/2.54), tight_layout=True, linewidth=2)
ax=fig.add_subplot(111)
for h in high:
yc=np.cos(phi)/(np.sin(phi)+np.sin(h))
r=np.cos(h)/(np.sin(phi)+np.sin(h))
circle(ax,[0,yc], r,theta)
for a in azimut:
xc=-np.cos(a)/(np.sin(a)*np.cos(phi))
yc=-np.tan(phi)
r=1/(np.sin(a)*np.cos(phi))
circle(ax,[xc,yc], r,theta,'gray',lw=0.2)
plt.xlim(-1,1)
plt.ylim(-1,2)
La dernière étape consiste à représenter l'araignée, c'est à dire la carte du ciel et l'écliptique. Pour cela on pourra reprendre le fichier Simbad_starsV6.txt du TN1:
VStar=pd.read_csv('Simbad_starsV6.txt',delimiter='|')
Tracer les 100 étoiles les plus brillantes du ciel à l'aide de leurs coordonnées équatoriales et représenter l'écliptique en rouge
selec=VStar.sort_values(...)[...]
Ecliptique=...
Lfigure=30#cm
fig=plt.figure(figsize=(Lfigure/2.54,Lfigure/2.54), tight_layout=True, linewidth=2)
ax=fig.add_subplot(111,projection='polar')
plt.plot(selec['RA'],selec['dec'],'*')
plt.plot(long_eq,Ecliptique,'r')
#plt.grid()
Assembler les 3 modules créer précedemment pour permettre une lecture de l'heure en fonction de la hauteur du soleil.
def astrolabe(h,Lat,long,date):
...
...
...
return heure
Vous pouvez trouver un simulateur d'Astrolabe ICI, à la latitude de Madrid.
Brigitte ALIX, réalise des Astrolabes modernes (c'est à dire avec la version actuelle du zodiaque et le point vernal dans les Poisson) [Presentation], son site commercialise de nombreux objets.
from IPython.core.display import HTML
style=open('notebooks.css', "r").read()
HTML(style)