S2_TD5_python_notebook_fr_student.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "F1WI4b8ShUv1"
   },
   "source": [
    "# <span style=\"color:blue\"> <center>Student / TD5 : 4TPU279U $-$ Bachelor 1st year $-$ spring 2023</center></span>\n",
    "# <center>Introduction to python programming</center>\n",
    "# <hr style=\"border:1px solid black\"><center>  Integration of movement equations </center><hr style=\"border:1px solid black\">\n",
    "</br>\n",
    "\n",
    "<div style=\"text-align: right\"> Credits: R. Boisgard, L. Truflandier, Philippe Paillou, Julien Burgin, Sara Zein, Leo Delmarre, Simon Villain-Guillot </div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following topics will be cover in this notebook:\n",
    "- Equations de la dynamique pour un mouvement plan\n",
    "- Champs gravitationnel et mouvement circulaire uniforme \n",
    "- Champs gravitationnel et mouvement elliptique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\">  Equations de la dynamique pour un mouvement plan <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'objectif ici est de  construire les trajectoires de points matériels soumis à diverses forces (ressort, poids, frottement, gravité) de façon numérique, c'est-à-dire sans résoudre analytiquement des équations différentielles de la dynamique :\n",
    "\n",
    "$$\\begin{align}\n",
    "\\vec f &=m \\vec a =m \\frac {d \\vec v}{dt}\\quad\\quad\\textcolor{green}{(8)}\\\\\n",
    "\\vec v &= \\frac {d \\vec r}{dt}\n",
    "\\end{align}$$\n",
    "\n",
    "avec $\\vec{v}$ le vecteur vitesse de coordonnées $(v_x,v_y,v_z)$ et $\\vec{r}$ le vecteur position de coordonnées $(x,y,z)$. A partir de la méthode d'Euler, cf. ***équation*** $\\textcolor{green}{(6)}$, le calcul de $\\vec{v}$ et $\\vec{r}$ peut s'obtenir par récurrence suivant :\n",
    "\n",
    "$$\\begin{align}\n",
    "\\vec{r}(t_{i+1}) &= \\vec{r}(t_i) + h\\ \\vec{v}(t_i)\\\\\n",
    "\\vec{v}(t_{i+1}) &= \\vec{v}(t_i) + \\frac{h}{m}\\vec{f}(t_i)\n",
    "\\end{align}$$\n",
    "\n",
    "avec $h=t_{i+1}-t_{i}$ le pas d'intégration des équations du mouvement. Dans le cas de mouvements  plans, \n",
    "on a :\n",
    "\n",
    "\n",
    "$$\\begin{align}\n",
    "  \\vec{r}= & \\begin{pmatrix}\n",
    "           x \\\\\n",
    "           y \\\\\n",
    "           0 \\\\\n",
    "         \\end{pmatrix},\\quad\n",
    "  \\vec{v}=  \\begin{pmatrix}\n",
    "           v_x \\\\\n",
    "           v_y \\\\\n",
    "           0 \\\\\n",
    "         \\end{pmatrix},\\quad\n",
    "  \\vec{f}=  \\begin{pmatrix}\n",
    "           f_x \\\\\n",
    "           f_y \\\\\n",
    "           0 \\\\\n",
    "         \\end{pmatrix}\n",
    "  \\end{align}$$ \n",
    "\n",
    "A partir de d'intégration numérique des équations du mouvement traité dans le TD 4, plus particulièrement \n",
    "de l'algorithme de Hooke-Newton, on obtient les 2 relations de récurrence suivantes pour la position : \n",
    "\n",
    "$$\\begin{align}\n",
    "x(t_{i+1})   &= x(t_i) + h\\;v_x(t_i)\\\\ \n",
    "y(t_{i+1})   &= y(t_i) + h\\;v_y(t_i)\\\\ \n",
    "\\end{align}$$\n",
    "\n",
    "et 2 relations pour la vitesse :\n",
    "\n",
    "$$\\begin{align}\n",
    "v_x(t_{i+1}) &= v_x(t_i)+\\frac{h}m f_x\\left(x(t_{i+1}),y(t_{i+1})\\right) \\\\ \n",
    "v_y(t_{i+1}) &= v_y(t_i)+\\frac{h}m f_y\\left(x(t_{i+1}),y(t_{i+1})\\right) \\\\ \n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\">  Champs gravitationnel et mouvement circulaire uniforme <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> On souhaite étudier le mouvement d'un objet de masse $m$ en orbite autour du soleil et soumis à la force gravitationnelle de Newton :\n",
    "$$\\vec f=-G\\frac{Mm}{r^3}\\vec r\\quad\\quad\\textcolor{green}{(1)}$$\n",
    "avec $G$ la constante gravitationnelle (m$^3$kg$^{-1}$s$^{-2}$) et $M$ la masse du soleil.\n",
    "On prendra comme unité de temps une année $T$ et comme unité de distance l'unité astronomique (ua), avec 1 ua = distance terre - soleil (m). Dans ce sytème d'unité :\n",
    "$$\\vec f=-\\frac{km}{r^3}\\vec r  \\quad\\text{avec}\\quad k=\\frac{GM}{\\text{ua}^3}T^2\\simeq (2\\pi)^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 5.1 : constante de Gauss <hr style=\"border:1px solid grey\">\n",
    "\n",
    "Démontrer numériquement que $\\sqrt{k}\\simeq 2\\pi$. On utilisera les constantes physiques `gravitational_constant` (ou `G`) et `astronomical_unit` (ou `au`) de la librairie `scipy.constants` avec la masse du soleil $M=1.9884\\times 10^{30}$ kg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.constants import au, G, pi\n",
    "from scipy.constants import Julian_year as T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Pour avoir une trajectoire circulaire, on prendra des conditions initiales telles que :\n",
    "$$Ec=\\frac{|Ep|}2=\\frac 12 \\frac{km}{r}=\\frac 12 \\frac{(2\\pi)^2m}r$$\n",
    "avec, pour $m=1$ kg : $r_0=1$ ua et $v_0=2\\pi$ ua/année. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 5.2 : orbite circulaire <hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Utiliser l'algorithme de Hook-Newton pour simuler la trajectoire de l'objet autour du soleil.\n",
    "- On discrétisera $t$ sur une année avec un pas $h=0.01$, ie. $t\\in[t_\\text{min},t_\\text{max}]$ avec\n",
    "$t_\\text{min} = 0$ et $t_\\text{max} = 1$\n",
    "- Les composantes de la positions $(x_0,y_0)=(1,0)$\n",
    "- Les composantes de la vitesse initiale seront $(v^0_x,v^0_y)=(0,2\\pi)$\n",
    "- Les composantes de la force $(f_x(t),f_y(t))$ sont données par\n",
    "$$f_x\\left(x,y\\right)=-\\frac{km}{r^3}x\\quad\\mathrm{et}\\quad f_y\\left(x,y\\right)=-\\frac{km}{r^3}y\\quad\\mathrm{avec}\\quad r=\\sqrt{\\left(x^2+y^2\\right)}$$\n",
    "On écrira 2 fonctions `fx`et `fy` à cet effet. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import arange, pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import pi, arange, zeros, sqrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Représenter graphiquement la variation de $v_x$ et $v_y$ en fonction du temps et comparer aux solutions anaytiques $v_x(t) = -r\\omega_0\\sin(\\omega_0 t)$ et $v_x(t) = +r\\omega_0\\cos(\\omega_0 t)$ avec $w = v_0/r_0 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "from numpy import cos, sin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Représenter graphiquement la trajectoire $(x(t),y(t))$ et les conditions initiales $(x_0(t),y_0(t))$. Ajouter un cercle de rayon r = 1 ainsi qu'un marker décrivant la position du soleil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "from numpy import linspace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Représenter graphiquement $x$, $y$ et $r=\\sqrt{x^2+y^2}$ en fonction de $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Représenter graphiquement $v_x$, $v_y$ et $v=\\sqrt{v_x^2+v_y^2}$ en fonction de $t$. Calculer la valeur moyenne temporelle de $v$ et son écart-type. Que ce passera t-il si on réduit le pas de temps $h$ et que l'on recommence la simulation ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import mean, std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. Représenter graphiquement l'énergie potentielle $E_p$, l'énergie cinétique $E_c$, et l'énergie totale $E=E_p+E_c$ en fonction de $t$, avec : $$E_p = -\\frac{km}{r}\\quad\\text{et}\\quad E_c =\\frac{1}{2}mv^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. Calculer les moyennes temporelles $\\langle E_p\\rangle$, $\\langle E_c\\rangle$ et montrer que $\\langle E_p\\rangle/\\langle E_c\\rangle=-2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <hr style=\"border:1px solid black\">  Champs gravitationnel et mouvement elliptique <hr style=\"border:1px solid black\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avec un repère centré sur le foyer, l'équation en coordonnée polaire d'une l'ellipse est donnée par :\n",
    "$$r(\\theta) = \\frac{(a-c)(a+c)}{a-c\\cos(\\theta)}\\quad\\textcolor{green}{(2)}$$\n",
    "avec \n",
    "- $a$ le demi-grand axe correspondant au diamètre de longueur maximale de l'ellipse\n",
    "- $c$ la distance entre le foyer et le centre de l'ellipse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Calcul du rayon $r$ de l'ellipse en fonction de $\\theta$ pour une période de révolution complète, ie. $\\theta=[0,2\\pi]$ avec $a=3$ et $c=2$ ua."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linspace, pi, cos\n",
    "\n",
    "a = 3.0\n",
    "c = 1.0\n",
    "theta = linspace(0,2*pi,256)\n",
    "r = (a - c)*(a + c)/(a - c*cos(theta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Représentation graphique d'une ellipse en coordonnées polaires"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "\n",
    "# Représentation de l'orbite\n",
    "plt.polar(theta,r)\n",
    "\n",
    "# Représentation de la distance a\n",
    "plt.polar([0,0],[c,c+a],marker='',color='black')\n",
    "plt.text(-0.0,c+a/2,'$a$',fontsize=12)\n",
    "\n",
    "# Représentation de la distance c\n",
    "plt.polar([0,0],[0,c],marker='',color='black')\n",
    "plt.text(-0.0,c/2,'$c$',fontsize=12)\n",
    "\n",
    "# Localisation du foyer de l'ellipse\n",
    "plt.polar(0,0,marker='o', color='orange',label='foyer', linestyle='')\n",
    "\n",
    "# Localisation du centre de l'ellipse\n",
    "plt.polar(0,c,marker='o', color='green', label='centre',linestyle='')\n",
    "\n",
    "# Localisation de l'apogée\n",
    "plt.polar(0,c+a ,marker='x',color='red', label='apogée', linestyle='')\n",
    "\n",
    "# Localisation du périgée\n",
    "plt.polar(pi,a-c,marker='x',color='blue',label='périgée',linestyle='')\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarques :\n",
    "- Le ***périgée*** correspond au point de l'orbite le plus proche du foyer. La distance périgée-foyer sera notée $d_{\\text{min}}$.\n",
    "- L'***apogée*** correspond au point de l'orbite le plus éloigné du foyer. La distance apogée-foyer sera noté $d_{\\text{max}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 5.3 : première loi de Kepler <hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reprendre l'***Exercise 5.2*** avec les conditions suivantes :\n",
    "- On discrétisera $t$ sur 4 années avec un pas $h=0.001$\n",
    "- Les composantes de la positions $(x_0,y_0)=(1,0)$\n",
    "- Les composantes de la vitesse initiale seront selectionnées de telle sorte que \n",
    "$$\\frac{|Ep|}2<Ec<|Ep|$$ \n",
    "soit $r_0=1$ et $2\\pi <v_0 <\\sqrt 8\\pi$. Par exemple : $(v^0_x,v^0_y)=(0,2.5\\pi)$.\n",
    "\n",
    "Comparer les résultats avec ceux obtenu pour l'orbite circulaire. Et vérifier visuellement la première loi de Kepler : *Les planètes tournent autour du Soleil en suivant des orbites en forme d'ellipse dont le Soleil occupe un des foyers.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 5.4 : troisième loi de Kepler  <hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Ecrire un algorithme permettant de calculer la distance maximum $d_{\\text{max}}$ entre le foyer (= le soleil) et l'ellipse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Ecrire un algorithme permettant de calculer la distance minimum $d_{\\text{min}}$ entre le foyer et l'ellipse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. En utilisant $d_{\\text{min}}$ et $d_{\\text{min}}$ calculer la longueur du demi-grand axe $a$ de l'ellipse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. A partir de La troisième loi de Kepler qui relit la période $P$ de l'orbite au demi-grand axe par\n",
    "$a^3/P^2=k$, calculer $P$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. En utilisant la valeur de la période $P$ pour discrétiser le temps, c'est-à-dire, $t\\in[0,P]$, tout en conservant les conditions initiales précédente reprendre l'***Exercise 5.3***. Pourquoi le théorème du viriel est de nouveau vérifié ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"border:1px solid grey\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> La deuxième loi de Kepler stipule que le rayon vecteur planète-Soleil balaie des aires égales en des temps égaux. Cette loi peut sécrire de la manière suivante :\n",
    "$$A(\\Delta t)=\\frac{||\\vec{L}||}{2m}\\Delta t=\\frac{1}{2}\\int^{t_1}_{t_0}r^2(t)\\dot\\theta(t)dt\\quad\\textcolor{green}{(3)}$$\n",
    "avec $\\vec{L}$ le moment cinétique orbital de l'objet de masse $m$ et $\\Delta t = t_1 -t_0$.\n",
    "\n",
    ">Pour vérifier cet énoncé nous devons calculer la vitesse angulaire $\\dot\\theta(t)={d\\theta}/{dt}$ à chaque instant $t$, puis réaliser l'intégration de l'***équation*** $\\textcolor{green}{(3)}$ numériquement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Exercise 5.6 : deuxième loi de Kepler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Répresenter sur un graphique l'ellipse numérique ainsi que les 4 vecteurs $(r_0,r_1,r_2,r_3)$ associée aux temps $t_0=0$, $t_1=P/16$, $t_2=P/2$ et $t_3=7P/16$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Entre l'apogée et le périgée, calculer $\\theta(t)$ à partir des coordonnées $x$ et $y$ et représenter les résultats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. En utilisant l'approximation des différences finies *centrées* calculer $\\dot\\theta(t)$. Les valeurs de $\\dot\\theta(t_0)$ et $\\dot\\theta(t_\\text{max})$ seront calculées à l'aide des différences finies *progressives*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Intégrer $r^2d\\dot\\theta$ entre $t_0=0$ et $t_1=P/16$. Comparer à l'intégrale obtenue pour $t_2=P/2$ et $t_3=7P/16$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Diaporama",
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}