{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Interpolation of the Trajectory of the Asteroid Pallas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Motivated by the discovery of the asteroids Ceres (1801) and Pallas (1802), Carl Friedrich Gauss\n", "studied the computation of planet trajectories in the beginning of the 19th century. There, he\n", "was faced with the following problem of trigonometric interpolation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpolation of the Asteroid's Trajectory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following data of the trajectory have been available to Gauss:\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "
 Ascension $\\theta$ (in degrees) Declination $X$ (in minutes) 0 30 60 90 120 150 180 210 240 270 300 330 408 89 -66 10 338 807 1238 1511 1583 1462 1183 804
\n", "\n", "Since the declination $X$ is periodic with regard to $\\theta$, the given\n", "trajectory data should be interpolated by the following trigonometric function:\n", "\n", "\\begin{equation}\n", " X(\\theta) = a_0\n", " + \\sum\\limits_{k=1}^5 \\left(\n", " a_k \\cos\\left( \\frac{2\\pi k \\theta}{360}\\right)\n", " +b_k \\sin\\left( \\frac{2\\pi k \\theta}{360}\\right)\n", " \\right)\n", " + a_6 \\cos \\left( \\frac{2\\pi \\cdot 6 \\theta}{360}\\right) \n", "\\end{equation}\n", "\n", "The data $X_l$ and $\\theta_l = 30l$ have to satisfy $X(\\theta_l) = X_l$\n", "for all $l = 0,\\ldots,11$. Thus,\n", "\n", "\\begin{equation} \\label{eq:interpol}\n", " X_l = a_0\n", " + \\sum\\limits_{k=1}^5 \\left(\n", " a_k \\cos\\left( \\frac{\\pi k l}{6}\\right)\n", " +b_k \\sin\\left( \\frac{\\pi k l}{6}\\right)\n", " \\right)\n", " + a_6 \\cos \\left( \\pi l \\right). \n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python Demo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a Python script (using IPython Notebook) to compute the\n", "coefficients $a_k$ and $b_k$.\n", "Plot the graph of the interpolated trajectory." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n", "import math\n", "import matplotlib.pyplot as plt\n", "\n", "X = [408, 89, -66, 10, 338, 807, 1238, 1511, 1583, 1462, 1183, 804]\n", "#X = [0, 89, -66, 10, 338, 807, 0, -807, -338, -10, 66, -89]\n", "#X = [50, 89, -66, 10, 338, 807, 22, 807, 338, 10, -66, 89]\n", "\n", "# find coefficients for interpolation function\n", "A = np.array([[math.cos(math.pi*k*l/6) for k in range(0,7)] for l in range(0,12)])\n", "B = np.array([[math.sin(math.pi*k*l/6) for k in range(1,6)] for l in range(0,12)])\n", "\n", "M = np.concatenate([A,B], axis=1);\n", "\n", "#print \"A = \" + str(A)\n", "print (\"A = \")\n", "for Ai in A:\n", " print ('%s' % ''.join(['%7.3f' % Aij for Aij in Ai]))\n", " \n", "#print \"B = \" + str(B)\n", "print (\"B = \")\n", "for Bi in B:\n", " print ('%s' % ''.join(['%7.3f' % Bij for Bij in Bi]))\n", "\n", "#print \"M = \" + str(M)\n", "print (\"M = \")\n", "for Mi in M:\n", " print ('%s' % ''.join(['%7.3f' % Mij for Mij in Mi]))\n", "\n", "#solve the system to get ab\n", "ab = np.linalg.solve(M,X)\n", "#separate a and b\n", "a = ab[0:7];\n", "b = ab[7:12];\n", "\n", "print (\"a:\", a)\n", "print (\"b:\", b)\n", "\n", "#generate interpolation point\n", "def Xint_gen(a,b,t):\n", " return np.sum(np.array([a[k]*math.cos(math.pi*k*t/180) for k in range(0,7)])) + \\\n", " np.sum(np.array([b[k-1]*math.sin(math.pi*k*t/180) for k in range(1,6)]))\n", " \n", "\n", "# plotting interpolation results and initial data\n", "resolution = 2\n", "plt.plot(range(0,720,resolution), [Xint_gen(a,b,t) for t in range(0,720,resolution)])\n", "plt.plot(range(0,360,30), X, 'o')\n", "plt.xlabel('ascension')\n", "plt.ylabel('declination')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Show that the interpolation problem in the previously used interpolation equation is equivalent to\n", "\n", "\\begin{equation} \\label{eq:complinter}\n", " X_l = \\sum\\limits_{k=-5}^6 c_k e^{i2\\pi kl/12},\n", "\\end{equation}\n", "\n", "if for $k=1,\\ldots,5$ $a_k$ and $b_k$ are chosen as $a_k = 2\\mathrm{Re}(\\{c_k\\})$\n", "and $b_k = -2\\mathrm{Im}(\\{c_k\\})$, while $c_0 = a_0$ and $c_6 = a_6$.\n", "\n", "Use the fact that all $X_l \\in \\mathbb{R}$ and, thus, that $c_{-k} = c_k^*$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python Demo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The equation from Excercise 1 also results from an interpolation problem\n", "with the complex interpolation function\n", "\\begin{equation}\n", " C(x) = \\sum\\limits_{k=-5}^6 c_k e^{ikx}\n", "\\end{equation}\n", "and the supporting points $x_n = 2\\pi n/N$.\n", "\n", "Use Python to compute and plot the interpolation function $C(x)$. Use the $a_k$ and $b_k$\n", "from Excercise 1 and construct the $C_k$ for all $k=-\\frac{N}{2}+1,\\ldots, \\frac{N}{2}$.\n", "\n", "Can $C(x)$ be used to describe the asteroid's trajectory?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#Note: we are in the same python environment as the previous python demo so X is still defined\n", "import cmath\n", "\n", "#compute c_k using the IDFT formula, only the c_k with k>=0 are required\n", "c = np.array([np.sum([X[l]*cmath.exp(-2*math.pi*k*l*1J/12)/12 for l in range(0,12)]) for k in range(0,7)])\n", "\n", "#compute a and b coefficients using the formula from exercise 1\n", "aFromC = c.real\n", "aFromC[1:6] = np.product([,aFromC[1:6]])\n", "bFromC = np.product([[-2],c[1:6].imag])\n", "\n", "# plotting interpolation results and initial data (reuse previously defined X_int_gen)\n", "plt.plot(range(0,720,resolution), [Xint_gen(aFromC,bFromC,t) for t in range(0,720,resolution)])\n", "plt.plot(range(0,360,30), X, 'o');\n", "plt.xlabel('ascension')\n", "plt.ylabel('declination')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Small Dictionary of Astronomy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Declination \n", "\n", " The angle between the celestial object and the celestial equator (projection\n", " of the earth equator on the celestial sphere).\n", "\n", " Ascension \n", "\n", " The angle between the First Point of Aries (The point where the ecliptic intersects\n", " the celestial equator) and the intersection point of the meridian of a celestial object\n", " and the celestial equator. It is equvalent to the geographical longitude but is\n", " measured to the east on the celestial equator. The units are usually hours, minutes\n", " and seconds, where 24 hours are equal to 360$^\\circ$.\n" ] } ], "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.0" } }, "nbformat": 4, "nbformat_minor": 0 }