{
 "metadata": {
  "name": "ex15a.ipynb"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Introduction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook runs through example 15a in [1]. I used this to debug the sun.py code and step through the errors in the tests. \n",
      "\n",
      "[1] Meeus, Jean. 1998. Astronomical Algorithms. 2nd edition. Richmond, Va: Willmann-Bell."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import astropy.units as u \n",
      "import astropy.coordinates as c\n",
      "import astropy.time as t\n",
      "import numpy as np\n",
      "import math"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "boston = c.EarthLocation.from_geodetic(-71.0833*u.deg,42.3333*u.deg)\n",
      "midnight = t.Time( '1988-03-20', scale='utc')\n",
      "h0 = -0.5667 *u.deg\n",
      "venus = c.SkyCoord(ra=41.73129*u.deg, dec=18.44092*u.deg)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The first problem is that the apparent sidereal times at greenwich for midnight of the day in question do not match. They match to within one arcsecond, but the precision reported in the book is 1/100th of an arcsecond."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sidereal_time = c.Angle( (11,50,58.10), unit=u.hour)\n",
      "print sidereal_time - midnight.sidereal_time('apparent','greenwich')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-0h00m00.2357s\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using the provided values, equation 15.2 computes the time of the transit. This agrees with the example to 0.4s."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m0 = (venus.ra - boston.longitude - sidereal_time)/(360*u.deg/u.sday)\n",
      "m0+=1*u.sday\n",
      "print m0 \n",
      "print (m0 - 0.81965*u.sday).to(u.s)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.819645851852 sday\n",
        "-0.357421412573 s\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculating the approximate rise and set times, we find that H0 and cos_H0 agree to within the precision reported by the Astronomical Algorithms book. m1 and m2 disagree by 0.8e-5 sidereal days, or 0.73s.\n",
      "\n",
      "To summarize, all of the m variables (time of the event expressed in fractional days) agree with the example to within 0.75s (8.5e-6 sidereal days)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def fix_day(day) : \n",
      "    if day > 1 * u.sday : \n",
      "        day -= 1 *u.sday\n",
      "    elif day < 0 * u.sday : \n",
      "        day += 1 * u.sday\n",
      "    return day\n",
      "\n",
      "\n",
      "cos_H0 = (np.sin(h0) - np.sin(boston.latitude)*np.sin(venus.dec))/(np.cos(boston.latitude)*np.cos(venus.dec))\n",
      "print cos_H0 # -0.3178735\n",
      "H0 = np.arccos(cos_H0)\n",
      "m1 = m0 - H0/(360*u.deg/u.sday)\n",
      "m2 = m0 + H0/(360*u.deg/u.sday)\n",
      "m = u.Quantity([fix_day(m0),fix_day(m1),fix_day(m2)])\n",
      "ref_m = [0.81965, 0.51817, 0.12113] * u.sday\n",
      "print H0.to(u.deg) # 108.5344\n",
      "print m\n",
      "print (m - ref_m)\n",
      "print (m - ref_m).to(u.s)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-0.317873477695\n",
        "108.534370473 deg\n",
        "[ 0.81964585  0.51816149  0.12113021] sday\n",
        "[ -4.14814815e-06  -8.51057216e-06   2.14275860e-07] sday\n",
        "[-0.35742141 -0.73330571  0.01846288] s\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now to calculate the corrections to the event times...Note that right off the bat, the sidereal times for the events (at least the transit), are greater than 360 degrees. To get the answer in the example, one needs to wrap the angle into 0-360 degrees.\n",
      "\n",
      "Calculated values agree to about the 10-15 arcsecond level. Note that from here on out, an array represents the three events. In index order, the events are: transit, rise, set."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#m_solar =  m * 1.0027379093604878 * u.day/u.sday\n",
      "m_solar =  m * 360.985647 * u.deg/u.sday \n",
      "#m_solar = m_solar * u.day/u.sday\n",
      "theta0 = sidereal_time + m_solar\n",
      "theta0.wrap_at(360*u.deg, inplace=True)\n",
      "print theta0.deg\n",
      "ref_theta0 = [113.62397, 4.79401, 221.46827] * u.deg\n",
      "print (theta0 - ref_theta0).deg\n",
      "print (theta0 - ref_theta0).to(u.arcsec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 113.62247147    4.79094384  221.4683521 ]\n",
        "[ -1.49852506e-03  -3.06615507e-03   8.21049532e-05]\n",
        "[u'-5.39469arcsec' u'-11.0382arcsec' u'0.295578arcsec']\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we repeat the above calculation using the given values for \"m\", the error is reduced substantially. Note that this error is accumulated over only one multiply and one addition.\n",
      "\n",
      "We establish the convention of prefixing variable names with \"cr\\_\" when that variable is calculated using the reference values taken from the example."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cr_m_solar = ref_m * 360.985647 * u.deg/u.sday\n",
      "cr_theta0 = sidereal_time + cr_m_solar\n",
      "cr_theta0.wrap_at(360*u.deg, inplace=True)\n",
      "print (cr_theta0-ref_theta0).deg\n",
      "print (cr_theta0-ref_theta0).to(u.arcsec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ -1.10311673e-06   6.03932333e-06   4.75444332e-06]\n",
        "[u'-0.00397122arcsec' u'0.0217416arcsec' u'0.017116arcsec']\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "test_theta0 = sidereal_time + 360.985647 * 0.12113 * u.deg\n",
      "print (test_theta0 - 221.46827 * u.deg).to(u.arcsec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.017116arcsec\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we calculate the local hour angle of the body at the times of the events. We use the given values for the interpolated sky coordinates. The error is exactly the same as for the error in our calculated theta0."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "venus_interp = c.SkyCoord(ra=[42.59324,42.27648,41.85927], dec=[0,18.64229,18.48835], unit=u.deg)\n",
      "print venus_interp\n",
      "H = theta0 + boston.longitude - venus_interp.ra\n",
      "print H.deg\n",
      "ref_H = [-0.05257, -108.56577, 108.52570] * u.deg\n",
      "print (H - ref_H).deg\n",
      "print (H - ref_H).to(u.arcsec)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<SkyCoord (ICRS): (ra, dec) in deg\n",
        "    [(42.59324, 0.0), (42.27648, 18.64229), (41.85927, 18.48835)]>\n",
        "[ -5.40685251e-02  -1.08568836e+02   1.08525782e+02]\n",
        "[ -1.49852506e-03  -3.06615507e-03   8.21049532e-05]\n",
        "[u'-5.39469arcsec' u'-11.0382arcsec' u'0.295578arcsec']\n"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Again, we repeat using the reference values, to see how much error is attributable to propagating errors from theta0. The error is reduced in all cases to essentially 0. Therefore, this calculation of hour angle does not introduce significant error."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cr_H = c.Angle(ref_theta0 + boston.longitude - venus_interp.ra)\n",
      "print cr_H.deg\n",
      "print (cr_H - ref_H).deg\n",
      "print (cr_H - ref_H).to(u.arcsec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ -5.25700000e-02  -1.08565770e+02   1.08525700e+02]\n",
        "[  4.21190860e-15   1.42108547e-14   0.00000000e+00]\n",
        "[u'1.51629e-11arcsec' u'5.11591e-11arcsec' u'0arcsec']\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we calculate venus's altitude above the horizon using eqn 13.6. This is not a relevant parameter for a transit, so ignore the first element in the array. Errors here seem to be influenced by the errors from theta0 as well. So we repeat, using the reference values. Again, this reduces the error to virtually zero."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sin_h = np.sin(boston.latitude)*np.sin(venus_interp.dec) + (np.cos(boston.latitude)*np.cos(venus_interp.dec)*np.cos(H))\n",
      "h = np.arcsin(sin_h)\n",
      "print h.to(u.deg)\n",
      "ref_h = [90, -0.44393, -0.52711] * u.deg\n",
      "print (h - ref_h).to(u.deg)\n",
      "print (h - ref_h).to(u.arcsec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 47.666672    -0.44596351  -0.52716223] deg\n",
        "[ -4.23333280e+01  -2.03351418e-03  -5.22255658e-05] deg\n",
        "[ -1.52399981e+05  -7.32065104e+00  -1.88012037e-01] arcsec\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cr_sin_h = np.sin(boston.latitude)*np.sin(venus_interp.dec) + (np.cos(boston.latitude)*np.cos(venus_interp.dec)*np.cos(ref_H))\n",
      "cr_h = np.arcsin(cr_sin_h)\n",
      "print cr_h.to(u.deg)\n",
      "print (cr_h - ref_h).to(u.deg)\n",
      "print (cr_h - ref_h).to(u.arcsec)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 47.66667353  -0.44392754  -0.52710764] deg\n",
        "[ -4.23333265e+01   2.46034892e-06   2.35656200e-06] deg\n",
        "[ -1.52399975e+05   8.85725610e-03   8.48362319e-03] arcsec\n"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we calculate first order corrections to the event times. Same thing: first with the values we've been calculating ourselves, next with the given values. In this case, however, the error in transit time correction is one order of magnitude higher when using the provided values. Rise and set time errors are approximately the same order of magnitude."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dm0 = - (H[0] / (360*(u.deg/u.sday)))\n",
      "dm1 = (h[1]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[1].dec) * np.cos(boston.latitude) * np.sin(H[1]))\n",
      "dm2 = (h[2]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[2].dec) * np.cos(boston.latitude) * np.sin(H[2]))\n",
      "delta_m = u.Quantity([dm0, dm1, dm2])\n",
      "ref_delta_m = [0.00015, -0.00051, 0.00017] *u.sday\n",
      "print delta_m.to(u.sday)\n",
      "print (delta_m - ref_delta_m).to(u.min)\n",
      "print (delta_m - ref_delta_m).to(u.s)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.00015019 -0.0005051   0.00016521] sday\n",
        "[ 0.00027335  0.00704102 -0.00687288] min\n",
        "[ 0.01640111  0.42246093 -0.41237303] s\n"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dm0 = - (ref_H[0] / (360*(u.deg/u.sday)))\n",
      "dm1 = (ref_h[1]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[1].dec) * np.cos(boston.latitude) * np.sin(ref_H[1]))\n",
      "dm2 = (ref_h[2]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[2].dec) * np.cos(boston.latitude) * np.sin(ref_H[2]))\n",
      "cr_delta_m = u.Quantity([dm0, dm1, dm2])\n",
      "print cr_delta_m.to(u.min)\n",
      "print (cr_delta_m - ref_delta_m).to(u.min)\n",
      "print (cr_delta_m - ref_delta_m).to(u.s)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.20970584 -0.73755733  0.23757199] min\n",
        "[-0.00570438 -0.00516256 -0.0065596 ] min\n",
        "[-0.34226292 -0.30975339 -0.39357613] s\n"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    }
   ],
   "metadata": {}
  }
 ]
}