{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fonctions d'une variable complexe (IV)\n",
    "\n",
    "# La sphère de Riemann\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "10 avril 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import cmath\n",
    "import math\n",
    "import colorsys\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Une quasi-bijection entre la sphère et le plan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 La sphère de Riemann $\\mathcal S^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $R>0$ un réel fixé. Soit $\\mathcal S^2$ l'ensemble des points $(X,Y,Z)\\in\\mathbb R^3$ tels que $X^2+Y^2+Z^2=R^2$. L'ensemble $\\mathcal S^2$ est la sphère de centre $O$ et de rayon R, nous l'appellerons la __sphère de Riemann__.\n",
    "\n",
    "__Remarque__ : Dans tous les livres parlant de fonctions de la variable complexe, le réel $R$ vaut 1. Pour des raisons qui apparaîtront plus tard, nous généralisons un peu ...\n",
    "\n",
    "Notons $N=(0,0,R)$ le __pôle nord__ de $\\mathcal S^2$. Soi $M=(X,Y,Z)\\in\\mathcal S^2$ un point différent de $N$. La droite $(NM)$ coupe le plan $xOy$, d'équation $z=0$, en un unique point. Quel est ce point ? Notons le $P=(x,y,0)$. Il existe un réel $\\lambda$ tel que $\\overrightarrow{NP}=\\lambda\\overrightarrow{NM}$ ce qui donne, en termes de coordonnées,\n",
    "\n",
    "$$x-0=\\lambda (X-0), y-0=\\lambda (Y-0), 0-R=\\lambda(Z-R)$$\n",
    "\n",
    "De la dernière égalité, on déduit \n",
    "\n",
    "$$\\lambda = \\frac{R}{R-Z}$$\n",
    "\n",
    "Puis, en reportant :\n",
    "\n",
    "$$x=\\frac{RX}{R-Z}, y=\\frac{RY}{R-Z}$$\n",
    "\n",
    "En identifiant le plan $xOy$ au plan complexe $\\mathbb C$, nous venons donc de définir une fonction $\\varphi:\\mathcal S^2\\setminus\\{N\\}\\to\\mathbb C$, définie par\n",
    "\n",
    "$$\\varphi(X,Y,Z)=\\frac{R}{R-Z}(X+iY)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def projection_C(X, Y, Z, R):\n",
    "    x = R * X / (R - Z)\n",
    "    y = R * Y / (R - Z)\n",
    "    return x + 1j * y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projection_C(1, 0, 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 La réciproque de $\\varphi$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IL est facile de voir que la fonction $\\varphi$ est bijective. Pour cela, prenons $z=x+iy\\in\\mathbb C$. Déterminons tous les points $M=(X,Y,Z)\\in\\mathcal S^2$, $M\\ne N$, tels que $\\varphi(M)=z$. Un tel point convient si et seulement si $(R-Z)x = RX$ et $(R-Z)y = RY$. Élevons au carré, additionnons, et rappelons-nous que $X^2+Y^2+Z^2=R^2$. il vient\n",
    "\n",
    "$$(R-Z)^2(x^2+y^2)=R^2(R^2-Z^2)$$\n",
    "\n",
    "ou encore, puisque $Z\\ne R$,\n",
    "\n",
    "$$(R-Z)(x^2+y^2)=R^2(R+Z)$$\n",
    "\n",
    "On en tire facilement \n",
    "\n",
    "$$Z=R\\frac{x^2+y^2-R^2}{x^2+y^2+R^2}$$\n",
    "\n",
    "Calculer $X$ et $Y$ est maintenant une formalité. Remarquons que $R-Z=\\frac {2R^3} {x^2+y^2+R^2}$. On en déduit\n",
    "\n",
    "$$X=\\frac{2R^2x}{x^2+y^2+R^2}\\quad{\\rm et}\\quad Y=\\frac{2R^2y}{x^2+y^2+R^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def anti_projection(z, R):\n",
    "    x, y = z.real, z.imag\n",
    "    r = x ** 2 + y ** 2 + R ** 2\n",
    "    return (2 * R ** 2 * x / r, 2 * R ** 2 * y / r, R - 2 * R ** 3 / r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si on compose `anti_projection` et `projection` on doit trouver l'identité."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 7\n",
    "X, Y, Z = anti_projection(1.234 - 5.678 * 1j, R)\n",
    "projection_C(X, Y, Z, R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enfin, disons l'identité à $\\varepsilon$ près :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Le point à l'infini du plan complexe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le pauvre pôle nord de la sphère n'a pas d'image par $\\varphi$. Soit $M=(X,Y,Z)\\in\\mathcal S^2$ un point différent de $N$. Que fait $\\varphi(M)$ lorsque $M$ s'approche indéfiniment de $N$, c'est à dire lorsque $Z\\to R$ ? On a\n",
    "\n",
    "$$|\\varphi(M)|^2=\\frac {R^2}{(R-Z)^2}(X^2+Y^2)=\\frac {R^2(R^2-Z^2)}{(R-Z)^2}=\\frac {R^2(R+Z)}{R-Z}$$\n",
    "\n",
    "Le numérateur de cette fraction tend vers $2R^3$ et son dénominateur tend vers 0. Ainsi, $|\\varphi(M)|$ tend vers $+\\infty$ lorsque $M$ tend vers $N$ : Le nombre complexe $\\varphi(M)$ \"part à l'infini\".\n",
    "\n",
    "Avez-vous noté les guillemets ? Nous allons les supprimer :-). Rajoutons au plan complexe un point que nous noterons $\\infty$. Notons $\\hat{\\mathbb C}=\\mathbb C\\cup\\{\\infty\\}$ notre \"nouveau\" plan : nous l'appellerons le __plan complexe étendu__. Maintenant, il est facile de prolonger $\\varphi$ en une bijection de $\\mathcal S^2$ sur $\\hat{\\mathbb C}$ en posant $\\varphi(N)=\\infty$. Maintenant, lorsque $M$ tend vers $N$, nous avons $\\varphi(M)$ qui tend vers l'infini. Qui tend ? Il faudrait pour cela pouvoir mesurer la distance entre $\\varphi(M)$ et $\\infty$, non ? Eh bien nous allons pour cela définir une distance adéquate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Une distance bornée sur le plan complexe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour $z,z'\\in\\hat{\\mathbb C}$, posons $\\delta(z,z')=d(M,M')$ où $M$ et $M'$ sont les antécédents de $z$ et $z'$ par $\\varphi$ et $d$ est la distance usuelle dans l'espace. Si vous préférez :\n",
    "\n",
    "$$\\delta(z, z')=d(\\varphi^{-1}(z),\\varphi^{-1}(z'))$$\n",
    "\n",
    "\n",
    "\n",
    "Comme $d$ est une distance et $\\varphi$ est bijective, la fonction $\\delta$ est une distance dans le plan complexe étendu. Et comme $M$ et $M'$ sont sont la sphère $\\mathcal S^2$, leur distance est inférieure ou égale au diamètre de la sphère, à savoir $2R$.\n",
    "\n",
    "Si vous le désirez, nous pouvons calculer explicitement $\\delta(z, z')$. Il suffit de reprendre les expressions de $M$ et $M'$ déjà trouvées. Je vous épargne les calculs. On obtient les résultats suivants :\n",
    "\n",
    "- Si $z,z'\\in\\mathbb C$, alors\n",
    "\n",
    "$$\\delta(z,z')=\\frac{2R^2|z'-z|}{\\sqrt{|z|^2+R^2}\\sqrt{|z'|^2+R^2}}$$\n",
    "\n",
    "- Si $z\\in\\mathbb C$, alors\n",
    "\n",
    "$$\\delta(z,\\infty)=\\frac{2R^2}{\\sqrt{|z|^2+R^2}}$$\n",
    "\n",
    "Remarquons avec la dernière égalité que si $|z|$ tend vers $+\\infty$, alors $\\delta(z,\\infty)$ tend vers 0. Nous pourrons donc enlever les modules et tout simplement dire que $z\\to\\infty$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_sphere(z1, z2, R):\n",
    "    r1 = math.sqrt(abs(z1) ** 2 + R ** 2)\n",
    "    r2 = math.sqrt(abs(z2) ** 2 + R ** 2)\n",
    "    return 2 * R ** 2 * abs(z2 - z1) / (r1 * r2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_sphere(1, 1j, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notez l'exemple ci-dessous :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_sphere(-1000, 1000, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous constatons que $\\delta(-1000,1000)$ est petit. Pourquoi ? Chacun de ces deux nombres est \"grand\", donc proche de $\\infty$. Par l'inégalité triangulaire appliquée à $\\delta$, ils sont donc proches l'un de l'autre."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Remarque__ : Pour les amateurs de topologie, nous disposons maintenant de __deux__ distances dans le plan complexe : la distance usuelle $d(z,z')=|z'-z|$ et la distance $\\delta$. On peut montrer que ces deux distances sont __équivalentes__, elles définissent la même topologie sur $\\mathbb C$.\n",
    "\n",
    "Par ailleurs, $\\hat{\\mathbb C}$, muni de la distance $\\delta$, devient un __espace métrique__. La sphère $\\mathcal S^2$, elle-aussi, est un espace métrique. On peut démontrer que la bijection $\\varphi:\\mathcal S^2\\to \\hat{\\mathbb C}$ est continue, et que sa réciproque est aussi continue. C'est ce que l'on appelle un __homéomorphisme__ entre $\\mathcal S^2$ et $\\hat{\\mathbb C}$. Nous pourrions aussi bien dire maintenant que $\\hat{\\mathbb C}$ __EST__ la sphère de Riemann.\n",
    "\n",
    "Je n'entrerai pas, excepté dans une ou deux occasions, dans les détails topologiques dans ce qui va suivre."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Dessiner la sphère"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Lancer de rayons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment dessiner la sphère de Riemann ? Imaginez un appareil photo placé sur l'axe $OX$, disons au point $\\Omega=(\\omega, 0, 0)$ où $\\omega>0$ est un réel fixé. L'appareil photo pointe vers la sphère de Riemann. À l'intérieur de l'appareil se trouve un ensemble de capteurs que nous représenterons par un plan vertical, parallèle à $yOz$, formé de pixels. Les rayons lumineux issus de la sphère frappent les pixels qui enregistrent la couleur de ces rayons. Nous avons obtenu une image.\n",
    "\n",
    "En fait, nous allons faire exactement le contraire : c'est la technique du __lancer de rayons__. Nous allons lancer des rayons (c'est à dire des demi-droites) partant du point $\\Omega$ vers tous les pixels de l'appareil. Si un tel  rayon frappe la sphère, on colorie le pixel avec la couleur du point frappé sur la sphère.\n",
    "\n",
    "Facile à dire, facile à faire ? Au travail. Mais commençons par un peu de géométrie."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Intersection d'un rayon et de la sphère"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour simplifier un peu nos calculs, nous supposons que le plan des pixels est le plan d'équation $x=0$. Nous appellerons ce plan le __plan de projection__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rappelons que $\\Omega=(\\omega, 0, 0)$. Soit $M=(0,y,z)$ un point dans le plan de projection. Soit $\\mathcal D$ la droite $(\\Omega M)$. Un point $P$ est sur $\\mathcal D$ si et seulement si il existe un réel $t$ tel que $\\overrightarrow{\\Omega P}=t\\overrightarrow{\\Omega M}$, ou encore $P=\\Omega + t\\overrightarrow{\\Omega M}$.\n",
    "\n",
    "\n",
    "En termes de coordonnées, cela donne $P=(\\omega-t\\omega, ty, tz)$ où $t\\in\\mathbb R$. Ce point est sur la sphère $\\mathcal S^2$ si et seulement si\n",
    "\n",
    "$$\\omega^2(1-t)^2+ t^2y^2+t^2z^2=R^2$$\n",
    "\n",
    "ou encore\n",
    "\n",
    "$$(\\omega^2 + y^2+z^2)t^2-2\\omega^2 t +\\omega^2-R^2=0$$\n",
    "\n",
    "Nous avons là une équation du second degré en $t$. Son discriminant est\n",
    "\n",
    "$$D=4(\\omega^4+(\\omega^2 + y^2+z^2)(R^2-\\omega^2))=4((y^2+z^2)(R^2-\\omega^2)+R^2\\omega^2)$$\n",
    "\n",
    "Évidemment, l'écran d'un appareil photo n'est pas un plan infini avec une densité infinie de pixels. Dans la fonction ci-dessous, nous supposons que l'écran est le carré $[-\\frac 6 5 R, \\frac 6 5 R[^2$, et qu'il contient $n$ pixels par ligne et par colonne. Le pixel ligne $i$, colonne $j$ est d'abord converti en un couple $(y, z)$ de réels, qui sont des coordonnées dans le plan d'équation $x=0$. Puis on résout l'équation ci-dessus. \n",
    "\n",
    "- Si le discriminant est strictement négatif, le rayon ne frappe pas la sphère et on renvoie la liste vide.\n",
    "- Sinon, le rayon frappe la sphère en un ou deux points. On renvoie une liste contenant un unique triplet de réels qui est le point frappé par le rayon : il s'agit du point qui est le plus proche de l'observateur, c'est à dire celui pour lequel $t$ est le plus petit.\n",
    "\n",
    "Comme j'imagine que vous savez résoudre une équation du second degré, je ne détaille pas plus avant la fonction `intersect_sphere`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def intersect_sphere(i, j, n, R, omega):\n",
    "    y = -1.2 * R + 2.4 * R * j / n\n",
    "    z = 1.2 * R - 2.4 * R * i / n\n",
    "    D = (y ** 2 + z ** 2) * (R ** 2 - omega ** 2) + R ** 2 * omega ** 2\n",
    "    if D < 0: return []\n",
    "    else:\n",
    "        d = math.sqrt(D)\n",
    "        t = (omega ** 2 - d) / (y ** 2 + z ** 2 + omega ** 2)\n",
    "        return [((1 - t) * omega, t * y, t * z)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Dessiner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction de dessin est maintenant immédiate. Elle prend un paramètre `methode` qui est une fonction qui renvoie une couleur (une couleur doit être un triplet $(r,g,b)$ de réels entre 0 et 1). Plus les paramètres nécessaires $n$, $R$ et `omega`. Si vous voulez éviter les surprises, prenz toujours `omega`$>R$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_riemann0(methode, n=500, R=5, omega=20):\n",
    "    ws = [[(0, 0, 0) for i in range(n)] for j in range(n)]\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            s = intersect_sphere(i, j, n, R, omega)\n",
    "            if s != []:\n",
    "                X, Y, Z = s[0]\n",
    "                ws[i][j] = methode(X, Y, Z, R)\n",
    "    plt.imshow(ws, interpolation='bicubic', extent=[-1.2, 1.2, -1.2, 1.2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tentons avec une méthode bête de coloration. On renvoie une couleur dépendant de $X$ et $Y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def methode_naive(X, Y, Z, R):\n",
    "    c1 = abs(math.sin(10 * (X + math.sin(5 * Y)) / R))\n",
    "    c2 = 1.\n",
    "    c3 = 1.\n",
    "    return (c1, 0., 0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_riemann0(methode_naive)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Rotations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voilà une bonne chose de faite, mais ce serait bien si nous pouvions faire tourner la sphère. Dans l'espace, on tourne autour d'une droite, d'un certain angle.\n",
    "\n",
    "Voici tout d'abord une fonction `rotation_Oz`. Elle prend en paramètres trois réels $x,y,z$ et un angle $t$. Il renvoie l'image du point $(x, y, z)$ par la rotation d'axe $Oz$ et d'angle $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotation_Oz(x, y, z, t):\n",
    "    c, s = math.cos(t), math.sin(t)\n",
    "    x1 = c * x - s * y\n",
    "    y1 = s * x + c * y\n",
    "    return (x1, y1, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation_Oz(1, 1, 1, math.pi / 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et voici de même une fonction qui fait tourner autour de l'axe $Oy$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotation_Oy(x, y, z, t):\n",
    "    c, s = math.cos(t), math.sin(t)\n",
    "    z1 = c * z - s * x\n",
    "    x1 = s * z + c * x\n",
    "    return (x1, y, z1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation_Oy(1, 1, 1, math.pi / 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Dessiner (bis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Réécrivons notre fonction de tracé en lui ajoutant deux paramètres `roty` et `rotz`. Ces deux paramètres contiennent les angles dont on veut faire tourner la sphère, __exprimés en degrés__. Soyons plus précis.\n",
    "\n",
    "1. On dessine la sphère en la coloriant avec la fonction `methode`.\n",
    "2. On la fait tourner autour de $Oz$ d'un angle `rotz`.\n",
    "3. On la fait tourner autour de $Oy$ d'un angle `roty`.\n",
    "\n",
    "Bon, ça c'est ce qu'on veut, mais on se rend vite compte que notre algorithme de lancer de rayons ne va pas du tout apprécier. Alors nous allons procéder ... à l'envers !\n",
    "\n",
    "\n",
    "1. On fait tourner la sphère autour de $Oy$ d'un angle $-$`roty`.\n",
    "\n",
    "2. On fait tourner la sphère autour de $Oz$ d'un angle $-$`rotz`.\n",
    "3. On dessine la sphère en la coloriant avec la fonction `methode`.\n",
    "\n",
    "Voici notre nouvelle fonction de tracé."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_riemann1(methode, roty=0, rotz=0, n=500, R=5, omega=20):\n",
    "    ws = [[(0, 0, 0) for i in range(n)] for j in range(n)]\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            s = intersect_sphere(i, j, n, R, omega)\n",
    "            if s != []:\n",
    "                X0, Y0, Z0 = s[0]\n",
    "                X, Y, Z = rotation_Oy(X0, Y0, Z0, -roty * math.pi / 180)\n",
    "                X, Y, Z = rotation_Oz(X, Y, Z, -rotz * math.pi / 180)\n",
    "                ws[i][j] = methode(X, Y, Z, R)\n",
    "    plt.imshow(ws, interpolation='bicubic', extent=[-1.2, 1.2, -1.2, 1.2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_riemann1(methode_naive, roty=45, rotz=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maintenant que nous savons dessiner la sphère $\\mathcal S^2$ nous allons utiliser notre tout nouvel outil pour dessiner le graphe de fonctions de $\\mathbb C$ vers $\\mathbb C$.\n",
    "\n",
    "À partir de maintenant, je vais supposer que vous avez lu au moins le premier notebook sur les fonctions holomorphes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6 La grille de coordonnées"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ce serait bien de pouvoir tracer sur la sphère une grille de coordonnées, un peu comme le fait `pyplot` avec la fonction `pyplot.grid`. Mais c'est quoi une grille de coordonnées sur la sphère de Riemann ? Voici l'idée : nous allons prendre une grille de coordonnées (droites horizontales et verticales) dans le plan complexe $\\mathbb C$. Puis nous allons dé-projeter cette grille sur la sphère $\\mathcal S^2$. Que devient une droite de $\\mathbb C$ lorsqu'on la dé-projette ? Eh bien une droite passe par $\\infty$. On devrait donc obtenir une courbe qui passe par le pôle nord de la sphère.\n",
    "\n",
    "On peut montrer, et je ne le ferai pas ici, que ces courbes sont des cercles passant par le pôle nord $N$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tout d'abord, comment tracer une \"droite\" ? La fonction `tracer_droite` ci-dessous fait le travail et trace la \"droite\" passant par $z\\in\\mathbb C$ et dirigée par $u\\in\\mathbb C^*$. Elle prend en paramètres\n",
    "\n",
    "- Une matrice `ws` représentant les pixels de l'écran.\n",
    "- Deux nombres complexes $z$ et $u\\ne 0$ : $z$ est un point de la droite et $u$ est un vecteur directeur de celle-ci.\n",
    "- Deux paramètres `roty` et `rotz` de rotation.\n",
    "- $R$, le rayon de la sphère, et `omega`, l'abscisse de l'observateur.\n",
    "\n",
    "La fonction déprojette un certain nombre de points de la droite puis projette les dé-projetés sur \"l'écran\" `ws` en mettant le pixel correspondant à la couleur blanche.\n",
    "\n",
    "__Exercice__ : Soit $\\Omega=(\\omega,0,0)$. Soit $M=(X,Y,Z)$. Que vaut le produit scalaire $<\\overrightarrow{\\Omega O},\\overrightarrow{OM}>$ ? Expliquez la raison d'être du test `if X * (X - omega) + Y ** 2 + Z ** 2 < 0:`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `scale` est une bijection de $]-1,1[$ sur $\\mathbb R$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scale(t):\n",
    "    return math.tan(math.pi * t / 2) ** 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tracer_droite(ws, z, u, roty, rotz, R, omega, n):\n",
    "    for t in range(-999, 1000):\n",
    "        w = z + u * scale(t / 1000)\n",
    "        X, Y, Z = anti_projection(w, R)\n",
    "        X, Y, Z = rotation_Oz(X, Y, Z, rotz * math.pi / 180)\n",
    "        X, Y, Z = rotation_Oy(X, Y, Z, roty * math.pi / 180)\n",
    "        if X * (X - omega) + Y ** 2 + Z ** 2 <= 0:\n",
    "            i, j = pixel(X, Y, Z, R, omega, n)\n",
    "            if i >= 0 and i < n and j >= 0 and j < n:\n",
    "                color = 1.\n",
    "                ws[i][j] = (color, color, color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pixel(X, Y, Z, R, omega, n):\n",
    "    y = omega * Y / (omega - X)\n",
    "    z = omega * Z / (omega - X)\n",
    "    j = int(0.5 + n * (y + 1.2 * R) / (2.4 * R))\n",
    "    i = int(0.5 + n * (1.2 * R - z) / (2.4 * R))\n",
    "    return (i, j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tracer la grille c'est simplement tracer quelques droites verticales et quelques droites horizontales."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tracer_grille(ws, roty, rotz, R, omega, n):\n",
    "    for x in range(8):\n",
    "        tracer_droite(ws, 2 ** x, 1j, roty, rotz, R, omega, n)\n",
    "        tracer_droite(ws, -2 ** x, 1j, roty, rotz, R, omega, n)\n",
    "    for y in range(8):\n",
    "        tracer_droite(ws, -1j * 2 ** y, 1, roty, rotz, R, omega, n)\n",
    "        tracer_droite(ws, 1j * 2 ** y, 1, roty, rotz, R, omega, n)\n",
    "    tracer_droite(ws, 0, 1j, roty, rotz, R, omega, n)\n",
    "    tracer_droite(ws, 0, 1, roty, rotz, R, omega, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et voici notre fonction ultime de tracé de $\\mathcal S^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_riemann(methode, roty=0, rotz=0, n=500, R=5, omega=20, axes=True):\n",
    "    ws = [[(0, 0, 0) for i in range(n)] for j in range(n)]\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            s = intersect_sphere(i, j, n, R, omega)\n",
    "            if s != []:\n",
    "                X0, Y0, Z0 = s[0]\n",
    "                X, Y, Z = rotation_Oy(X0, Y0, Z0, -roty * math.pi / 180)\n",
    "                X, Y, Z = rotation_Oz(X, Y, Z, -rotz * math.pi / 180)\n",
    "                ws[i][j] = methode(X, Y, Z, R)\n",
    "    if axes: tracer_grille(ws, roty, rotz, R, omega, n)\n",
    "    plt.imshow(ws, interpolation='bicubic', extent=[-1.2, 1.2, -1.2, 1.2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allons-y ! Voici la sphère vue du dessus. La grille de coordonnées est composée de deux familles de cercles passant par le pôle nord. Dans chacune des familles les cercles sont tangents entre-eux au pôle nord. Remarquons également (sans preuve ici) que les cercle d'une famille coupent les cercles de l'autre famlille à angle droit en __DEUX__ points : le pôle nord, et un autre point. Si vous pensez que ces cercles sont des droites dans le plan complexe, vous obtenez le résultat suivant : \n",
    "\n",
    "- Deux droites non parallèles se coupent en __DEUX__ points, l'un des deux étant $\\infty$.\n",
    "- Deux droites parallèles se coupent en __UN__ point, $\\infty$, et sont tangentes en ce point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_riemann(methode_naive, roty=90, rotz=0, axes=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notre méthode de tracé de la grille n'est pas idéale. Par exemple, certaines parties des \"droites\" sont en pointillés. Il faudrait beaucoup plus de travail pour obtenir une représentation plus esthétique. Nous nous contenterons de ce que nous avons fait :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Graphe d'une fonction de variable complexe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Rappels du premier notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction ci-dessous renvoie $\\frac 1 {2\\pi}\\arg z$ où l'argument de $z$ est dans $[0,2\\pi[$. Elle renvoie donc un réel entre 0 et 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def arg01(z):\n",
    "    theta = cmath.phase(z) / (2 * math.pi)\n",
    "    if theta < 0: return theta + 1\n",
    "    else: return theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J'ai parlé de la fonction `psi` dans le premier notebook sur les fonctions holomorphes. Elle prend un réel $x>0$ et deux réels $m$ et $d$, et renvoie un nombre entre 0 et 1 particulièrement bien adapté pour colorier le graphe d'une fonction complexe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def psi(x, m, d):\n",
    "    if x == 0: return 0\n",
    "    else:\n",
    "        return abs(math.sin(math.pi * math.log(x, 2) / m)) ** d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 La méthode de tracé"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les sphères que nous avons tracées pour l'instant étaient toutes grises. L'idée, maintenant, est la suivante. Étant donnée une fonction $f:\\mathbb C\\to\\mathbb C$, chaque point $M$ de la sphère de Riemann est colorié en fonction de la valeur de $\\varphi(M)$. \n",
    "\n",
    "- On projette le point $M$ sur le plan complexe étendu. Soit $z$ son projeté.\n",
    "- On calcule $w=f(z)$.\n",
    "- À partir de $w$ on calcule un triplet représentant une couleur dans le modèle HSV (voir le premier notebook).\n",
    "- On renvoie ce triplet, ce sera la couleur choisie pour $M$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def methode_hsv(X, Y, Z, R, f, m, d):\n",
    "    try:\n",
    "        z = projection_C(X, Y, Z, R)\n",
    "        w = f(z)\n",
    "        r = psi(abs(w), m, d)\n",
    "        ph = arg01(w)\n",
    "        return colorsys.hsv_to_rgb(ph, 1, r)\n",
    "    except: return (1., 1., 1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "À quoi sert le `try ... except` ? Les raisons pour que les calculs ci-dessous échouent sont nombreuses. L'une d'entre-elles, et non des moindres, est que nous calculons des quantités du genre $f(z)$ pour $z$ voisin de l'infini. Essayez donc d'évaluer `math.exp(1000)`. Et pourtant, 1000 n'est pas un bien grand nombre."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 La fonction de tracé"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction de tracé d'une fonction $f:\\mathbb C\\to \\mathbb C$ est maintenant immédiate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_complex(f, roty=0, rotz=0, m=1, d=0.2, n=500, R=5, omega=20, axes=True):\n",
    "    plot_riemann(lambda X, Y, Z, R:methode_hsv(X, Y, Z, R, f, m, d), roty, rotz, n, R, omega, axes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Un premier test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons la fonction $f:z\\mapsto z^2-1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:z ** 2 - 1, roty=0, rotz=0, R=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La sphère de Riemann est coloriée par des bandes horizontales sur chacune desquelles ont lieu des cycle de couleurs. Peut-on voir les deux racines de $f$ en même temps ? Oui, il suffit de prendre $R$ un peu plus grand et effectuer une rotation par rapport à $Oy$ de $-\\frac \\pi 2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:z ** 2 - 1, roty=-90, rotz=0, R=5, omega=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ah mais oui, c'est ce que nous avions vu dans le premier notebook !\n",
    "\n",
    "Ah mais non, les couleurs tournent dans le mauvais sens ???\n",
    "\n",
    "__Réfléchissez un peu : la projection de la sphère de Riemann sur le plan complexe renverse les sens de rotation__\n",
    "\n",
    "__Conséquence__ :\n",
    "\n",
    "- Racine : les couleurs tournent dans le sens trigonométrique inverse.\n",
    "- Pôle : les couleurs tournent dans le sens trigonométrique direct.\n",
    "\n",
    "Que se passe -t-il en $\\infty$ ? Tournons de $\\frac \\pi 2$ autour de $Oy$, l'infini sera alors juste devant nos yeux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:z ** 2 - 1, roty=90, rotz=0, R=5, omega=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Deux__ cycles de couleurs, parcourus dans le sens trigonométrique ? Oserons-nous ?\n",
    "\n",
    "__Propriété__ : $\\infty$ est un pôle double de $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 L'infini comme pôle ou racine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Définition__ : On dit que $\\infty$ est un pôle de $f$ lorsque 0 est un pôle de $g:z\\mapsto f(\\frac 1 z)$. L'ordre de $\\infty$ comme pôle de $f$ est alors l'ordre de 0 en tant que pôle de $g$. On définit de même la notion de racine en $\\infty$ et l'ordre de multiplicité de la racine.\n",
    "\n",
    "Prenons par exemple $f(z)=z^2-1$, prolongée à $\\hat{\\mathbb C}$ par $f(\\infty)=\\infty$. On a \n",
    "\n",
    "$$g(z)=\\frac 1 {z^2} - 1$$\n",
    "\n",
    "et 0 est un pôle d'ordre 2 de $g$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.6 Un exemple où l'infini est une racine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons $f(z)=\\frac 1 {z^3-1}$.\n",
    "\n",
    "$1$, $j$ et $j^2$ sont maintenant trois pôles de $f$ d'ordre 1 : notez le cycle de couleurs dans le sens trigonométrique autour de ces deux pôles. Et $-i$ et $i$ sont deux racines d'ordre 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:(z ** 2 + 1) / (z ** 3 - 1), roty=-90, rotz=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et l'infini ? Eh bien c'est une racine d'ordre 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:(z ** 2 + 1) / (z ** 3 - 1), roty=90, rotz=0, axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une preuve ? On a \n",
    "\n",
    "$$f(\\frac 1 z)=\\frac{\\frac 1 {z^2}+1}{\\frac 1 {z^3}-1}=\\frac{z(1+z^2)}{1-z^3}$$\n",
    "\n",
    "et 0 est bien racine simple de $z\\mapsto f(\\frac 1 z)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Et maintenant ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Je ne vais pas reprendre tous les exemples de fonctions que nous avons tracées dans les trois premiers notebooks : \n",
    "\n",
    "- polynômes\n",
    "- fractions rationnelles\n",
    "- exponentielle\n",
    "- logarithmes\n",
    "- puissances\n",
    "- fonctions trigonométriques\n",
    "\n",
    "__Faites-le__. Et regardez l'infini :-).\n",
    "\n",
    "Histoire de vous donner envie :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Logarithmes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici une vue \"standard\" du logarithme, avec sa racine simple 1, sa singularité non isolée $0$ et une partie de la coupure $\\mathbb R_-$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:cmath.log(z), roty=-90, rotz=0, m=0.25, axes=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici une vue centrée sur (presque tous) les réels négatifs. Eh oui, malgré toute notre bonne volonté il nous est très difficile de voir en même temps le pôle sud et le pôle nord de $\\mathcal S^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:cmath.log(z), roty=180, rotz=0, m=0.25, axes=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et enfin, une vue centrée sur $\\infty$. Mettez `axes` à `True` si vous désirez l'affichage de la grille."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:cmath.log(z), roty=90, rotz=0, m=0.25, axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le logarithme est holomorphe sur $\\mathcal S^2$ privée d'un demi-cercle, allant du pôle sud (0) au pôle nord ($\\infty$). Ceci explique les remarques que j'avais faites en leur temps sur l'existence de déterminations du logarithme. C'est à ce moment que j'avais parlé de logarithmes fous, définis sur $\\mathbb C$ privé d'une spirale logarithmique. Eh bien cette spirale joint 0 à $\\infty$. Dessinons-le, ce logarithme bizarre."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logfou(z, a):\n",
    "    r = abs(z)\n",
    "    theta = cmath.phase(z)\n",
    "    k = math.floor((math.log(r, a)- theta) / (2 * math.pi))\n",
    "    return math.log(r) + 1j * (theta + 2 * k * math.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici le logarithme fou autour de $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:logfou(z, 1.1), roty=-90, rotz=0, m=0.25, axes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et le voici autour de $\\infty$ !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z:logfou(z, 1.1), roty=90, rotz=0, m=0.2, axes=False, d=0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Sinus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Terminons en beauté par un zoom à l'infini de la fonction sinus : $\\infty$ est une singularité essentielle de $\\sin$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_complex(lambda z: cmath.sin(z), roty=90, rotz=0, m=5, R=2, d=0.2, axes=False, n=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 À vous de jouer ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
