Lab 2.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.interpolate as interpol"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f(theta, tau = 365.25636, ecen = 0.01671123):\n",
      "    return tau/(2*pi)*(2*arctan2((1-ecen)**0.5*sin(theta/2),(1+ecen)**0.5*cos(theta/2)))-e*sin(theta)*(1-ecen**2)**0.5/(1+ecen*cos(theta))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = linspace(0,2*pi)\n",
      "y = interpol.InterpolatedUnivariateSpline(x,f(x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xx = linspace(0,2*pi,5000)\n",
      "vv = y(xx,nu = 1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(xx,f(xx))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 61,
       "text": [
        "[<matplotlib.lines.Line2D at 0x10e6a7990>]"
       ]
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def t(psi, w=2*pi/365.25636, ecen = 0.01671123):\n",
      "    return (psi - ecen*sin(psi))/w\n",
      "def theta(psi, ecen= 0.01671123):\n",
      "    thta = arccos((cos(psi)-ecen)/(1-ecen*cos(psi)))\n",
      "#    if 0<psi and psi<=pi:\n",
      " #      return thta\n",
      "  #  if pi < psi and psi <= 2*pi:\n",
      "   #     return 2*pi-thta\n",
      "    #else:\n",
      "     #   print \"Invalid Domain\"\n",
      "      #  return 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 111
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def q(a):\n",
      "    print (a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 119
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q(3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3\n"
       ]
      }
     ],
     "prompt_number": 120
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "zeros_like?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 121
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "clf()\n",
      "Tau = 365.25636\n",
      "e=0.01671123\n",
      "omega=2*pi/Tau\n",
      "from scipy.interpolate import InterpolatedUnivariateSpline\n",
      "def f(theta, tau = 365.25636, ecen = e):\n",
      "    return tau/(2*pi)*((2*arctan2((1-ecen)**0.5*sin(theta/2),(1+ecen)**0.5*cos(theta/2)))-e*sin(theta)*(1-ecen**2)**0.5/(1+ecen*cos(theta)))\n",
      "theta = linspace(0,2*pi,5000)\n",
      "time = linspace(0,Tau,5000)\n",
      "t2 = InterpolatedUnivariateSpline(f(theta),theta)\n",
      "plot(time,t2(2*pi*time/Tau, nu=1)*180/pi,'b--',label = \"Conservation of Momentum\")\n",
      "\n",
      "def ti(psi, w=omega, ecen = e):\n",
      "    return (psi - ecen*sin(psi))/w\n",
      "Psispline = InterpolatedUnivariateSpline(ti(theta),theta,k=3)\n",
      "def thta0(x, ecen = e):\n",
      "    return arccos((cos(x)-ecen)/(1-ecen*cos(x)))\n",
      "def thtapi(x, ecen = e):\n",
      "    return 2*pi - arccos((cos(2*pi - x)-ecen)/(1-ecen*cos(2*pi - x)))\n",
      "Psifill = zeros_like(theta)\n",
      "Psifill[:theta.size/2]=thta0(Psispline(time))[:theta.size/2]\n",
      "Psifill[theta.size/2:]=thtapi(Psispline(time))[theta.size/2:]\n",
      "\n",
      "omegaarray = linspace(omega, omega, 5000)\n",
      "t3 = InterpolatedUnivariateSpline(Psifill,theta)\n",
      "time1 = linspace(0,Tau,4999)\n",
      "plot(time,omegaarray*180/pi,'g:',label = 'Average Angular Velocity')\n",
      "plot(time1, t3(2*pi*time1/Tau,nu=1),'r-',label = \"Kepler's Equation\")\n",
      "xlim(0,366)\n",
      "legend()\n",
      "xlabel('Time (Days)')\n",
      "ylabel('Angular Speed (degrees/day)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 49,
       "text": [
        "<matplotlib.text.Text at 0x110214d50>"
       ]
      }
     ],
     "prompt_number": 49
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "clf()\n",
      "semilogy(time1,abs(t2(2*pi*time1/Tau, nu=1)*180/pi -t3(2*pi*time1/Tau,nu=1)))\n",
      "xlabel('Time (Days)')\n",
      "ylabel('Discrepency in Angular Speed (degrees/day)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 54,
       "text": [
        "<matplotlib.text.Text at 0x10b36c890>"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "semilogy?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}