« Python pour le calcul scientifique/Éléments de programmation » : différence entre les versions
Contenu supprimé Contenu ajouté
Ligne 696 :
{{clear}}
{{Boîte déroulante/début |titre=Analyse d'optique géométrique}}
[[Fichier:Lentille hemispherique analyse geometrique.svg|vignette|Analyse géométrique du problème.]]
Ligne 706 :
: ''h'' = ''d'' ⋅ tan θ.
L'angle d'incidence vaut θ. D'après la loi de Snell-Descartes, l'angle de réfraction θ<sub>2</sub> vaut :
: θ<sub>2</sub> = arcsin(
Le rayon réfracté passe par le points de coordonnées (0, ''h''). L'équation de la droite est donc :
: ''y'' = a ⋅ ''x'' + ''h''
Ligne 720 :
ce qui nous permet de calculer cet angle :
: <math>\theta_\mathrm{i} = \operatorname{arcos} \left ( \frac{x_\mathrm{M} + a \cdot y_\mathrm{M}}{\mathrm{R} \cdot \sqrt{1^2 + a^2} } \right )</math>
Comme nous passons vers un milieu d'indice plus faible, il y a un risque de réflexion totale. L'angle limite est :
: θ<sub>max</sub> = arcsin(1/''n'').
Si l'on a θ<sub>i</sub> > θ<sub>max</sub>, le rayon repart vers l'intérieur. Nous ne traçons pas le rayon car cela nous emmènerait trop loin dans l'analyse. En revanche, si θ<sub>i</sub> ≤ θ<sub>max</sub>, alors nous pouvons appliquer la loi de Snell-Descartes pour avoir l'angle de réfraction θ<sub>e</sub> :
: θ<sub>e</sub> = arcsin(''n'' ⋅ sin θ<sub>i</sub>).
Pour tracer le rayon sortant, il nous faut l'angle θ<sub>3</sub> par rapport à l'horizontale. L'angle du rayon [OM] par rapport à l'horizontal vaut arctan(''y''<sub>M</sub> / ''x''<sub>M</sub>), nous avons donc
: θ<sub>3</sub> = arctan(''y''<sub>M</sub> / ''x''<sub>M</sub>) + θ<sub>e</sub>.
{{Boîte déroulante/fin}}
{{Boîte déroulante/début |titre=Solution}}
Nous créons une fonction <code>refrac()</code> qui permet de calculer l'angle réfracté à partir de l'angle d'incidence et des indices de réfraction. S'il y a rélexion totale, alors nous générons une erreur.
La fonction <code>lanceRayon()</code> calcule les différents points de passage du rayon. Elle appelle pour cela la fonction <code>refrac()</code>. Si un appel de la commande <code>refrac()</code> génère une erreur, alors nous générons également une erreur.
Nous déterminons l'angle d'émision du rayon <code>thetaLimite</code> qui provoque une réflecxion totale. Pour cela, nous créons une fonction <code>rechercheLimite()</code> qui cherche par dichotomie :
* nous partons de l'angle maximum possible, lorsque le rayon frappe le sommet de la lentille, et nous appelons <code>lanceRayon()</code> ; si cela ne génère pas d'erreur, alors nous pouvons aller jusqu'à cette valeur, la recherche est terminée ; si cela génère une erreur, alors nous divisons la valeur par deux ;
* à une étape de la recherche donnée, si <code>lanceRayon()</code> ne génère pas d'erreur avec l'angle testé, alors nous savons que l'angle limite est supérieur à cette valeur ; cette valeur est donc la valeur basse recherchée ; si au contraire <code>lanceRayon()</code> génère une erreur, alors c'est que l'angle est trop important, nous nous rapprochons donc de la valeur basse ;
* nous nous arrêtons lorsque les valeurs haute et basse sont suffisamment proche.
Nous traçons un rayon tous les 5° jusqu'à la valeur limite.
<syntaxhighlight lang="python">
#!/usr/bin/env python3
# coding: utf-8
# ******************************************************
# ******************************************************
# ** Lancer de rayons pour une lentille hémisphérique **
# ******************************************************
# ******************************************************
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as nppol
# **************
# * Constantes *
# **************
# Pour la conversion degrés ↔ radians
radversdeg = 180/np.pi
degversrad = 1/radversdeg
epsilon = np.finfo(float).eps # précision maximale
# *************
# * Fonctions *
# *************
def boucleEntreeNombre(messageSaisie, valeurDefaut):
"""Permet de s’assurer que l’utilisateur·rice a bien entré un nombre.
Entrée :
— message à afficher (chaîne de caractères) ;
— valeur par défaut (réel à virgule flottante).
Sortie : nombre (réel à virgule flottante)."""
messageErreur = "Veuillez entrer une valeur numérique (ou vide pour accepter la valeur par défaut).\n"
execute = True
while execute:
strNombre = input(messageSaisie+f" (valeur par défaut {valeurDefaut}) : ")
if strNombre == "":
nombre = valeurDefaut
execute = False
else:
try:
nombre = float(strNombre)
except:
print(messageErreur)
else:
execute = False
return nombre
def initialisation():
"""L’utilisateur·rice entre les variables du problème.
Entrées : aucune.
Sorties :
— R1 (mm) : rayon de la lentille ;
— d (mm) : distance de la source au dioptre plan ;
— n (sans dimension) : indice de réfraction du verre."""
R1 = boucleEntreeNombre("Rayon de la lentille en mm", 20.0)
d = boucleEntreeNombre("Distance de la source au dioptre plan en mm", 20.0)
n = boucleEntreeNombre("Indice de réfraction (sans dimension)", 1.5)
return (R1, d, n)
def refrac(n1, n2, theta1):
"""Calcule l’angle de réfraction theta2 (radians) en fonction
— de l’angle d’incidence theta1 (radians);
— de l’indice de réfraction n1 du premier milieu ;
— de l’indice de réfraction n2 du second milieu."""
reflexionTotale=False
rapport=n2/n1
rapportinv=np.reciprocal(rapport)
if n1 > n2:
thetal = np.arcsin(rapport) # angle limite pour la réflexion totale
if theta1 >= thetal:
reflexionTotale=True
if reflexionTotale:
print("Réflexion totale")
raise ValueError
else:
return np.arcsin(rapportinv*np.sin(theta1))
def lanceRayon(n1, n2, d, R, theta1):
"""Détermine le rayon issu de la source
située à une distance d (mm) du bareau
et avec une élévation de theta (radians),
en fonction des indices de réfraction n1 et n2.
Les éléments retournés sont :
— la hauteur h (mm) à laquelle le rayon frappe le barreau ;
— l’angle de réfraction theta2 (radians)) dans le barreau ;
— l’angle de réfraction theta3 (radians) à la sortie du barreau
— le point M(x, y) (mm) auquel le rayon sort du barreau."""
h = d*np.tan(theta1)
if h >= R:
print("Le rayon est au-dessus du barreau")
raise ValueError
else:
theta2 = refrac(n1, n2, theta1)
a = np.tan(theta2)
x = max(nppol.polyroots([h*h - R*R, 2*a*h, 1+a*a])) # recherche de l'intersection du rayon avec le cercle
y = a*x + h
M = np.array([x, y])
thetaint = np.arccos((x + a*y)/(R*np.sqrt(1 + a*a)))
theta3 = np.arctan(y/x) - refrac(n2, n1, thetaint)
return (h, theta2, theta3, M)
def rechercheLimite(n1, n2, d, R):
"""Recherche l’angle limite pour la réflexion totale.
Entrée :
— indice de réfraction des milieux 1 et 2, n1 et n2 ;
— distance au barreau, d(mm).
Sortie : angle limite theta (radians)"""
angleHaut = np.arctan(R/d)
angleBas = 0
angleTest = angleHaut
try:
lanceRayon(n1, n2, d, angleTest, R)
except:
condition = True # il y a réflexion total en haut de la lentille
else:
condition = False # il n’y a jamais réflexion totale dans la lentille
while condition: #dichotomie
angleTest = np.mean([angleHaut, angleBas]) # on ajuste la valeur de test
try:
lanceRayon(n1, n2, d, R, angleTest)
except:
angleHaut = angleTest # réflexion totale : on abaisse la valeur maximale
else:
angleBas = angleTest # pas de réflexion totale : on monte la valeur minimale
condition = ((angleHaut - angleBas) >= 0.001) # on a cerné la limite à 0,001 rad près
if not condition:
angleTest = angleBas
return angleTest
# ***********************
# * Programme principal *
# ***********************
(R1, d, n) = initialisation()
xmax = round(R1 + d)
thetaLimite = rechercheLimite(1, n, d, R1)
thetaLimiteDeg = thetaLimite*radversdeg
print(f"Angle limite pour la réflexion totale : {thetaLimiteDeg:.1f}°.\n")
anglesDeg = np.arange(0, thetaLimiteDeg, 5)[1:] # trace un rayon tous les 5°
anglesRad = anglesDeg*degversrad
nb = len(anglesDeg)
h = np.zeros(nb) # initialisation des vecteurs de valeurs
theta2 = np.zeros(nb)
theta3 = np.zeros(nb)
M = np.zeros((nb, 2))
for i in range(nb):
(h[i], theta2[i], theta3[i], M[i, :]) = lanceRayon(1, n, d, R1, anglesRad[i])
(h_lim, theta2_lim, theta3_lim, M_lim) = lanceRayon(1, n, d, R1, thetaLimite)
# tracé
anglesCercle = 0.5*np.pi*(np.linspace(1, 0, 20))
x_cercle = R1*np.cos(anglesCercle) # coordonnées des pints du cercle
y_cercle = R1*np.sin(anglesCercle)
fig = plt.plot([-d,xmax], [0, 0], "k-.", linewidth="0.5") # tracé de l'axe optique
for i in range(nb):
plt.plot([-d, 0, M[i, 0], xmax], [0, h[i], M[i, 1], M[i, 1] + (xmax - M[i, 0])*np.tan(theta3[i])],
label=f"{anglesDeg[i]:.0f}°")
plt.plot([-d, 0, M_lim[0], xmax], [0, h_lim, M_lim[1], M_lim[1] + (xmax - M_lim[0])*np.tan(theta3_lim)],
label=f"{0.1*int(np.trunc(10*thetaLimite*radversdeg)):.1f}°")
plt.plot(x_cercle, y_cercle, "k", linewidth="0.5") # tracé du cercle
plt.plot([0,0], [0, R1], "k", linewidth="0.5") # tracé du premier dioptre
#plt.axis("square")
plt.gca().set_aspect("equal", adjustable="box")
plt.xlabel("x (mm)")
plt.ylabel("y (mm)")
plt.title("Lentille hémisphérique, lancer de rayons")
plt.legend()
plt.savefig("lentille_hemispherique_lancer_rayon.svg", format="svg")
plt.show()
</syntaxhighlight>
{{Boîte déroulante/fin}}
| |||