1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
| #!/usr/env python3
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
##############################################
# Constantes physiques
##############################################
g = 9.81 # m/s2
lat = 48*2*np.pi/360 # rad, latitude (à Paris : 48°)
Tjour = 3600*24 # s, durée d'un jour
#Tjour = 360 #s, durée de dix minutes, qui rendent l'observation rapide
Omega = 2*np.pi/Tjour
##############################################
# Paramètres du pendule
##############################################
t0 = 0 # s, temps pour lancer le scipt en retard
l = 67 # m, longueur du pendule
x0 = 15 # m, position initiale historique
omega = np.sqrt(g/l)
T = 2*np.pi/omega
print(r'PĂ©riode du pendule: {:.2f}'.format(T))
##############################################
# Calculs
##############################################
def xy(t):
x0cosomegat = x0 * np.cos(omega*t)
x = x0cosomegat * np.sin(Omega*t)
y = - x0cosomegat * np.cos(Omega*t)
return x, y
##############################################
# Animation
##############################################
dt = 50 # ms, delta entre les temps d'affichage
fig = plt.figure()
plt.title(r'Simulation du pendule du pantheon Ă Paris')
ttxt = plt.text(0.05, 0.95, r't = 0', transform=plt.gca().transAxes, va='top')
# Limites
max = x0 * 1.1
plt.xlim(-max, max)
plt.ylim(-max, max)
plt.axes().set_aspect('equal')
plt.xticks([])
plt.yticks([])
# Axe de référence et point initial
plt.axvline(x=0, color='k')
plt.plot(*xy(0), linestyle='none', marker='o', color='k')
# Point
point, = plt.plot([], [], linestyle='none', marker='o', color='b')
# Variables pour la détermination de la période lors de l'animation
# Le pendule Ă parcouru une demi-oscillation si
# dr = r(t) - r(t-dt)
# augmente alors qu'il diminuait juste avant
rprev = 0
drprev = 0
dr_up = False
def animate(i):
""" Fonction pour l'animation """
# Calcul
ti = dt*i*1e-3 + t0
x, y = xy(ti)
# Mise Ă jour du plot
point.set_data(x, y)
ttxt.set_text(r'$t = {:.2f}$ s'.format(ti))
# Détermination de la période du pendule
r = np.sqrt(np.power(x, 2) + np.power(y, 2))
global rprev, drprev, dr_up
dr = r - rprev
if dr_up and dr < drprev:
print(r'Demi-période à {:.2f} +/- {:.2f}'.format(ti, dt*1e-3))
print(r' x = {:.4f}'.format(x))
dr_up = False
if dr_up == False and dr > drprev:
dr_up = True
drprev = dr
anim = FuncAnimation(fig, animate, interval=dt)
# Affichage
plt.show()
|