{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Root finding methods\n",
    "\n",
    "#### scipy.optimize documentation at https://docs.scipy.org/doc/scipy/reference/optimize.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1) Relaxation method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "def relax(func,x0,tol):\n",
    "    \n",
    "    # Finds x such that func(x) = x, \n",
    "    # via relaxation method.\n",
    "    # start is the starting point of the recursion\n",
    "    # tol is the tolerance\n",
    "    \n",
    "    ynew = func(x0)\n",
    "    yold = start\n",
    "    eps = abs(ynew - yold)\n",
    "        \n",
    "    while (eps > tol):\n",
    "        y = func(ynew)\n",
    "        eps = abs(y - ynew)\n",
    "        ynew, yold = y, ynew\n",
    "        \n",
    "        #print(ynew)\n",
    "        \n",
    "    return ynew"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise: solve  $2 - e^{-x} = x$\n",
    "\n",
    "1) Make a plot to find the total number of roots of the equation\n",
    "\n",
    "2) Find all roots using the relaxation method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x): return 2.0 - np.exp(-x)\n",
    "\n",
    "start = 4.0\n",
    "tol = 1e-6\n",
    "\n",
    "# Plotting\n",
    "f = [ ]\n",
    "g = [ ]\n",
    "x = np.arange(-3.0,4.0,0.1)\n",
    "for t in x:\n",
    "    f.append(func(t))\n",
    "    g.append(t)\n",
    "plt.xlim(-2,4)\n",
    "plt.ylim(-6,5)\n",
    "plt.plot(x,f,label = 'y=f(x)')\n",
    "plt.plot(x,g,label = 'y=x')\n",
    "plt.xlabel('x', size=16)\n",
    "\n",
    "# Solving\n",
    "root = relax(func,start,tol)\n",
    "# Plotting solution\n",
    "plt.axvline(x=root,ls='dashed',color='r',lw=1,label='roots')\n",
    "\n",
    "print('Root 1 = ',  \"{:.6f}\".format(root))\n",
    "\n",
    "# Since |df/dx| > 1 in the negative solution, the relaxation algorithm does \n",
    "# not converge there. Need to find a suitable transformation and write the\n",
    "# equation in the form g(x) = x, where |dg/dx| < 1.\n",
    "# This can be done by rearranging the equation as  x = - log(2-x)\n",
    "\n",
    "def func(x): return - np.log(2-x)\n",
    "\n",
    "start = -1.0\n",
    "tol = 1e-6\n",
    "\n",
    "# Solving\n",
    "root = relax(func,start,tol)\n",
    "# Plotting solution\n",
    "plt.axvline(x=root,ls='dashed',color='r',lw=1)\n",
    "plt.legend()\n",
    "\n",
    "print('Root 2 = ', \"{:.6f}\".format(root))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 6.10. Solve $x = 1 - e^{-cx}$ for c from 0 to 3 in steps of 0.01 and plot x(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using relaxation method\n",
    "\n",
    "# Input function\n",
    "def func(x,c):\n",
    "    f = 1 - np.exp(-c*x)\n",
    "    return f\n",
    "\n",
    "# Problem parameters\n",
    "start = 10.0\n",
    "tol = 1e-6\n",
    "cmin = 0.\n",
    "cmax = 3.\n",
    "step = 0.01\n",
    "\n",
    "roots = [ ]\n",
    "\n",
    "for c in np.arange(cmin,cmax,step):\n",
    "\n",
    "    # Calling function func for fixed parameter c\n",
    "    # To be passed to function relax\n",
    "    f = lambda x: func(x,c) \n",
    "    # Finding root \n",
    "    root = relax(f,start,tol)\n",
    "    # Updating list of roots for varying c\n",
    "    roots.append(root)\n",
    "\n",
    "plt.plot(np.arange(cmin,cmax,step),roots)\n",
    "plt.xlabel('c', size = 16)\n",
    "plt.ylabel('x(c)', size = 16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2) Bisection mehod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "def bisection(func,a,b,tol):\n",
    "    \n",
    "    # Finds a solution of func(x) = 0, via\n",
    "    # bisection method\n",
    "    # a and b are the boundaries of the search interval [a,b]\n",
    "    # tol is the tolerance\n",
    "    \n",
    "    if (a > b):\n",
    "        a,b = b,a\n",
    "  \n",
    "    x1 = a\n",
    "    x2 = b\n",
    "    f1 = func(x1)\n",
    "    f2 = func(x2)\n",
    "    \n",
    "    if (f1/f2 > 0.):\n",
    "        sys.exit('f(x) does not have opposite signs at the boundaries')\n",
    "    else:  \n",
    "        while (abs(x1-x2) > tol):\n",
    "            x3 = (x1 + x2)/2.0\n",
    "            f3 = func(x3)\n",
    "            if (f3/f1 < 0.):\n",
    "                x2 = x3\n",
    "            else:\n",
    "                x1 = x3\n",
    "                \n",
    "    return ((x1+x2)/2.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solve again $2 - e^{-x} = x$, but now use the bisection method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Testing\n",
    "\n",
    "def func(x): return 2.0 - np.exp(-x) - x\n",
    "\n",
    "tol = 1e-6\n",
    "\n",
    "# Solving for root 1\n",
    "a = 0.1\n",
    "b = 10\n",
    "root1 = bisection(func,a,b,tol)\n",
    "\n",
    "# Solving for root 2\n",
    "a = -20.\n",
    "b = 1.5\n",
    "root2 = bisection(func,a,b,tol)\n",
    "\n",
    "print('Root 1 = ', root1)\n",
    "print('Root 2 = ', root2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting\n",
    "f = [ ]\n",
    "x = np.arange(-4.0,10.0,0.1)\n",
    "for t in x:\n",
    "    f.append(func(t))\n",
    "plt.plot(x,f)\n",
    "plt.xlim(-3,6)\n",
    "plt.ylim(-20,5)\n",
    "plt.axhline(linestyle='dashed',color = 'red', lw = 1)\n",
    "plt.axvline(x=root1,ls='dashed',color='r',lw=1)\n",
    "plt.axvline(x=root2,ls='dashed',color='r',lw=1)\n",
    "plt.xlabel('x',size=16)\n",
    "plt.ylabel('f(x)',size=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 6.13: Wien's displacement constant\n",
    "\n",
    "Blackbody law: $I(\\lambda) =\\frac{2 \\pi h c^2 \\lambda^{-5}}{e^{{hc /\\lambda k_B T} } - 1}$.\n",
    "\n",
    "1) Differentiate and show that the maximum of emitted radiation is the solution of the equation $5 e^{-x} + x - 5 = 0$, where \n",
    "$x = hc/\\lambda k_B T$. Therefore the peak wavelenght is $\\lambda = b/T$, with $b = hc/k_B x$.\n",
    "\n",
    "2) Solve the equation above for x with binary search (bisection), with accuracy 1e-6.\n",
    "\n",
    "3) Estimate the temperature of the sun, knowing that the peak wavelength of emitted radiation is $\\lambda = 502$ nm.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x): return 5.*np.exp(-x) + x - 5.\n",
    "\n",
    "a = 2.0\n",
    "b = 10.\n",
    "tol = 1e-6\n",
    "\n",
    "# Solving\n",
    "root = bisection(func,a,b,tol)\n",
    "\n",
    "print('Root =', root, '   f(root)=', func(root))\n",
    "\n",
    "# Estimating sun's temperature\n",
    "# Parameters (SI units)\n",
    "h = 6.6261e-34\n",
    "c = 2.998e8\n",
    "kB = 1.381e-23\n",
    "wavelength = 5.02e-7\n",
    "b = (h*c)/(kB*root)\n",
    "T = b/wavelength\n",
    "\n",
    "print('Estimated temperature of the sun: T = ',\"{:.0f}\".format(T), 'K')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting\n",
    "xlist = [ ]\n",
    "flist = [ ]\n",
    "for x in np.arange(2.,10.,0.1):\n",
    "    xlist.append(x)\n",
    "    flist.append(func(x))\n",
    "    \n",
    "plt.plot(xlist,flist)\n",
    "plt.xlabel('x',size=16)\n",
    "plt.ylabel('$f(x) = 5e^{-x} + x -5$',size=14)\n",
    "plt.axhline(linestyle='dashed',color = 'red',lw=1)\n",
    "plt.axvline(x=root,linestyle='dashed',color='red',lw=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below I am solving the same exercise, using functions from scipy.optimize, instad of my own"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.optimize as opt\n",
    "\n",
    "def func(x): return 5.*np.exp(-x) + x - 5.\n",
    "\n",
    "a = 2.0\n",
    "b = 10.\n",
    "root = opt.root_scalar(func, method='bisect', bracket=[a,b])\n",
    "\n",
    "print(root)\n",
    "\n",
    "# If you want specific items from the output object root,\n",
    "# do as follows (some examples; you can uncomment to see the result)\n",
    "\n",
    "#print(root.root)\n",
    "#x = root.iterations\n",
    "#print(x)\n",
    "#print(root.flag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another syntax\n",
    "root = opt.bisect(func,a,b)\n",
    "print('Root =', root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3) Newton's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton(func,deriv,x0,tol=1e-4):\n",
    "    \n",
    "    # Finds a solution of func(x)=0 using Newton's method\n",
    "    # deriv is an input function which computes df/dx\n",
    "    # x0 is the starting point of the search\n",
    "    # tol is the tolerance, set by default at 1e-4\n",
    "    \n",
    "    demomode = False\n",
    "    \n",
    "    x = x0\n",
    "    acc = tol + 1.\n",
    "    \n",
    "    while (acc > tol):\n",
    "        delta = func(x)/deriv(x)\n",
    "        x1 = x - delta\n",
    "        if demomode:\n",
    "            print('root = ', x1, 'acc = ', abs(delta) )\n",
    "        acc = abs(x1 - x)\n",
    "        x = x1\n",
    "    \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise: assuming circular orbits and that the Earth is much more massive than the moon, find the L1 point of a satellite orbiting between the two\n",
    "\n",
    "1) Show that, in L1: $\\frac{GM}{r^2} - \\frac{Gm}{(R-r)^2} = \\omega^2 r$, where $R$ is the moon-earth distance, M and m are the respective masses, $G$ is Newton's constant and $\\omega$ is the angular velocity of both the moon and the satellite \n",
    "\n",
    "2) Write a program that solves the equation above with at least four significant figures, using Newton's or secant's method. Parameters (SI units):\n",
    "   $G = 6.674 \\times 10^{-11}$, $M= 5.974 \\times 10^{24}$, $m = 7.348 \\times 10^{22}$, $R = 3.844 \\times 10^8$, \n",
    "   $\\omega =    2.662 \\times 10^{-6}$.\n",
    "\n",
    "3) Make a plot to verify that your code computed the correct solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Parameters\n",
    "G = 6.674e-11\n",
    "M = 5.974e+24\n",
    "m = 7.348e+22\n",
    "R = 3.844e+8\n",
    "omega = 2.662e-6\n",
    "\n",
    "# Function\n",
    "f = lambda r: (G*M)/(r**2.) - (G*m)/(R-r)**2. - omega**2*r \n",
    "# Derivative\n",
    "deriv = lambda r: -(2.*G*M)/(r**3.) - (2.*G*m)/(R-r)**3. - omega**2.\n",
    "\n",
    "# Solving for position of L1, using Newton's method\n",
    "start = 2.*R/3.\n",
    "L1 = newton(f,deriv,start,tol=1e-6)\n",
    "\n",
    "print(' ')\n",
    "print('Position of L1 in units of R = ', L1/R)\n",
    "print(' ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below I am solving the same exercise, now using functions from scipy optimize instead of my own"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize as opt\n",
    "\n",
    "# Parameters\n",
    "G = 6.674e-11\n",
    "M = 5.974e+24\n",
    "m = 7.348e+22\n",
    "R = 3.844e+8\n",
    "omega = 2.662e-6\n",
    "\n",
    "# Function\n",
    "f = lambda r: (G*M)/(r**2.) - (G*m)/(R-r)**2. - omega**2*r \n",
    "# Derivative\n",
    "deriv = lambda r: -(2.*G*M)/(r**3.) - (2.*G*m)/(R-r)**3. - omega**2.\n",
    "\n",
    "start = 2.*R/3.\n",
    "\n",
    "# Newton's method\n",
    "root = opt.root_scalar(f, method='newton', fprime=deriv, x0 = start)\n",
    "\n",
    "print(root)\n",
    "print(' ')\n",
    "\n",
    "# Extracting items from object root, example:\n",
    "print('The solution in units of R is: x/R = ', root.root/R )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another syntax\n",
    "root = opt.newton(f,start,deriv)\n",
    "\n",
    "print('The solution in units of R is: x/R = ', root/R )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting\n",
    "fig,ax = plt.subplots()\n",
    "pos = np.arange(0.5*R,0.95*R,0.01*R)\n",
    "x = np.arange(0.5,0.95,0.01)\n",
    "fx = f(pos)\n",
    "ax.axhline(linestyle='dashed',color = 'red', linewidth=1)\n",
    "ax.axvline(x=L1/R, linestyle='dashed',color='red',linewidth=1)\n",
    "ax.set_xlabel('x/R',size=16)\n",
    "ax.set_ylabel('f(x)', size=16)\n",
    "ax.annotate('$L_1$',xy=(2,0.),xycoords='axes points',xytext=(270,-15),color='red',size=14)\n",
    "                                                    \n",
    "ax.plot(x,fx)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4) Secant method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton_secant(func,x0,deriv=None,x1=None,tol=1e-4):\n",
    "    \n",
    "    # Finds a solution of func(x)=0 using Newton's method\n",
    "    # deriv is an input function which computes df/dx\n",
    "    # x0 is the starting point of the search\n",
    "    # tol is the tolerance, set by default at 1e-4\n",
    "    \n",
    "    demomode = True\n",
    "    \n",
    "    acc = tol + 1.\n",
    "    \n",
    "    # If deriv is not passed, use secant\n",
    "    if (deriv == None):\n",
    "        \n",
    "        if demomode:\n",
    "            print('Using secant method')\n",
    "        \n",
    "        # If x1 is not passed for secant method,\n",
    "        # choose x1 = x0 + 1 by default\n",
    "        if (x1 == None):\n",
    "            x1 = x0 + 1.\n",
    "            \n",
    "        while (acc > tol):\n",
    "            fprime = (func(x1) - func(x0))/(x1 - x0)\n",
    "            delta = func(x1)/fprime\n",
    "            x2 = x1 - delta\n",
    "            if demomode:\n",
    "                print('root = ', x2, 'acc = ', abs(delta))\n",
    "            acc = abs(x2 - x1)\n",
    "            x0,x1 = x1,x2\n",
    "       \n",
    "        x = x2\n",
    "    \n",
    "    # If deriv is passed, use Newton\n",
    "    else:\n",
    "        \n",
    "        if demomode:\n",
    "            print('Using Newton method')\n",
    "        \n",
    "        x = x0\n",
    "    \n",
    "        while (acc > tol):\n",
    "            delta = func(x)/deriv(x)\n",
    "            x1 = x - delta\n",
    "            if demomode:\n",
    "                print('root = ', x1, 'acc = ', abs(delta))\n",
    "            acc = abs(x1 - x)\n",
    "            x = x1\n",
    "            \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Finding moon-earth L1 point, using secant method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "G = 6.674e-11\n",
    "M = 5.974e+24\n",
    "m = 7.348e+22\n",
    "R = 3.844e+8\n",
    "omega = 2.662e-6\n",
    "\n",
    "# Function\n",
    "f = lambda r: (G*M)/(r**2.) - (G*m)/(R-r)**2. - omega**2*r \n",
    "# Derivative\n",
    "deriv = lambda r: -(2.*G*M)/(r**3.) - (2.*G*m)/(R-r)**3. - omega**2.\n",
    "\n",
    "# Solving for position of L1, using Newton's method\n",
    "start = R/100.\n",
    "\n",
    "# Uses secant method with default x1\n",
    "#L1 = newton_secant(f,start,tol=1e-6)\n",
    "# Uses secant method with x1 passed by user\n",
    "L1 = newton_secant(f,start,x1=start+R/10.,tol=1e-6)\n",
    "# Uses Newton's method\n",
    "#L1 = newton_secant(f,start,deriv,tol=1e-6)\n",
    "\n",
    "print(' ')\n",
    "print('The root in units of R is: x/R = ', L1/R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize as opt\n",
    "\n",
    "# Parameters\n",
    "G = 6.674e-11\n",
    "M = 5.974e+24\n",
    "m = 7.348e+22\n",
    "R = 3.844e+8\n",
    "omega = 2.662e-6\n",
    "\n",
    "# Function\n",
    "f = lambda r: (G*M)/(r**2.) - (G*m)/(R-r)**2. - omega**2*r \n",
    "# Derivative\n",
    "deriv = lambda r: -(2.*G*M)/(r**3.) - (2.*G*m)/(R-r)**3. - omega**2.\n",
    "\n",
    "start = 2.*R/3.\n",
    "\n",
    "# Secant method, using scipy\n",
    "root = opt.root_scalar(f, method='secant', x0 = start, x1 = start + 1.)\n",
    "\n",
    "print(root)\n",
    "print(' ')\n",
    "print('The solution in units of R is: x/R = ', root.root/R )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another syntax.\n",
    "# If you do not pass fprime to opt.newton, it automatically uses the secant method\n",
    "root = opt.newton(f,start)\n",
    "\n",
    "print('The solution in units of R is: x/R = ', root/R )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary exercise from slides: eccentric anomaly\n",
    "\n",
    "Solve the exercise using relaxation, bisection and Newton's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Relaxation, bisection and Newton's root finding functions, as coded above \n",
    "###########################################################################\n",
    "def relax(func,x0,tol):\n",
    "    \n",
    "    demomode = False\n",
    "        \n",
    "    # Finds x such that func(x) = x, \n",
    "    # via relaxation method.\n",
    "    \n",
    "    ynew = func(x0)\n",
    "    yold = start\n",
    "    eps = abs(ynew - yold)\n",
    "        \n",
    "    while (eps > tol):\n",
    "        y = func(ynew)\n",
    "        eps = abs(y - ynew)\n",
    "        ynew, yold = y, ynew\n",
    "        \n",
    "        if demomode:\n",
    "            print(ynew, eps)\n",
    "            \n",
    "    return ynew\n",
    "\n",
    "def bisection(func,a,b,tol):\n",
    "    \n",
    "    # Finds solution of f(x) = 0, via\n",
    "    # bisection method\n",
    "    \n",
    "    if (a > b):\n",
    "        a,b = b,a\n",
    "  \n",
    "    x1 = a\n",
    "    x2 = b\n",
    "    f1 = func(x1)\n",
    "    f2 = func(x2)\n",
    "    \n",
    "    if (f1/f2 > 0.):\n",
    "        sys.exit('f(x) does not have opposite signs at the boundaries')\n",
    "    else:  \n",
    "        while (abs(x1-x2) > tol):\n",
    "            x3 = (x1 + x2)/2.0\n",
    "            f3 = func(x3)\n",
    "            if (f3/f1 < 0.):\n",
    "                x2 = x3\n",
    "            else:\n",
    "                x1 = x3\n",
    "                \n",
    "    return ((x1+x2)/2.0)\n",
    "\n",
    "\n",
    "def newton(func,deriv,x0,tol=1e-4):\n",
    "    \n",
    "    demomode = False\n",
    "    \n",
    "    x = x0\n",
    "    acc = tol + 1.\n",
    "    \n",
    "    while (acc > tol):\n",
    "        delta = func(x)/deriv(x)\n",
    "        x1 = x - delta\n",
    "        if demomode:\n",
    "            print('root = ', x1, 'acc = ', delta)\n",
    "        acc = abs(x1 - x)\n",
    "        x = x1\n",
    "    \n",
    "    return x\n",
    "#######################################################################\n",
    "\n",
    "\n",
    "\n",
    "# Solving the exercise using all three methods above\n",
    "#######################################################\n",
    "\n",
    "mean_anomaly = [np.pi,np.pi/3.]\n",
    "eccentricity = [0.1,0.7,0.9]\n",
    "\n",
    "# Solving: Newton\n",
    "\n",
    "def func(E,ma,ecc): return E - ecc*np.sin(E) - ma\n",
    "def deriv(E,ecc): return 1 - ecc*np.cos(E)    \n",
    "    \n",
    "newt = [ ]\n",
    "start = np.pi/2\n",
    "\n",
    "print(' ')\n",
    "print('NEWTONS METHOD, SOLUTIONS:')\n",
    "print('==========================')\n",
    "print('ecc = eccentricty')\n",
    "print('ma = mean anomaly') \n",
    "print('E = eccentric anomaly')\n",
    "print('--------------------------------------------')\n",
    "\n",
    "for ma in mean_anomaly:\n",
    "    \n",
    "    der = lambda x: deriv(x,ecc)\n",
    "    \n",
    "    for ecc in eccentricity:\n",
    "            f = lambda x: func(x,ma,ecc)\n",
    "            newt.append(newton(f,der,start,tol=1e-6))\n",
    "            \n",
    "            print('ecc =', (\"{:.6f}\".format(ecc)), ' ma =', (\"{:.6f}\".format(ma)), ' E =', (\"{:.6f}\".format(newt[-1]) ))\n",
    "\n",
    "                        \n",
    "# Solving: bisection\n",
    "\n",
    "def func(E,ma,ecc): return E - ecc*np.sin(E) - ma\n",
    "\n",
    "a = -100.\n",
    "b = 100.\n",
    "tol = 1e-6\n",
    "bisec = [ ]\n",
    "\n",
    "print(' ')\n",
    "print('BISECTION METHOD, SOLUTIONS:')\n",
    "print('==========================')\n",
    "print('ecc = eccentricty')\n",
    "print('ma = mean anomaly') \n",
    "print('E = eccentric anomaly')\n",
    "print('--------------------------------------------')\n",
    "\n",
    "for ma in mean_anomaly:\n",
    "    for ecc in eccentricity:\n",
    "        \n",
    "            f = lambda x: func(x,ma,ecc)\n",
    "            bisec.append(bisection(f,a,b,tol))\n",
    "            \n",
    "            print('ecc =', (\"{:.6f}\".format(ecc)), ' ma =', (\"{:.6f}\".format(ma)), ' E =', (\"{:.6f}\".format(bisec[-1]) ))\n",
    "\n",
    "            \n",
    "# Solving: relaxation\n",
    "\n",
    "start = np.pi/2.\n",
    "tol = 1e-6\n",
    "rel = [ ]\n",
    "\n",
    "def func(E,ma,ecc): return ma + ecc*np.sin(E)\n",
    "\n",
    "print(' ')\n",
    "print('RELAXATION METHOD, SOLUTIONS:')\n",
    "print('==========================')\n",
    "print('ecc = eccentricty')\n",
    "print('ma = mean anomaly') \n",
    "print('E = eccentric anomaly')\n",
    "print('--------------------------------------------')\n",
    "\n",
    "for ma in mean_anomaly:\n",
    "    for ecc in eccentricity:\n",
    "        \n",
    "            f = lambda x: func(x,ma,ecc)\n",
    "            rel.append(relax(f,start,tol))\n",
    "            \n",
    "            print('ecc =', (\"{:.6f}\".format(ecc)), ' ma =', (\"{:.6f}\".format(ma)), ' E =', (\"{:.6f}\".format(rel[-1]) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}