solar_system.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['jh', 'hj l;']\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "00001    3.34  0.12 K172G 266.80164   72.90774   80.30986   10.59240  0.0756828  0.21403411   2.7678809  0 MPO389571  6651 113 1801-2016 0.60 M-v 30h MPCLINUX   0000      (1) Ceres              20161108\n",
      "\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "Bogus: 00001    3.34  0.12 K172G 266.80164   72.90774   80.30986   10.59240  0.0756828  0.21403411   2.7678809  0 MPO389571  6651 113 1801-2016 0.60 M-v 30h MPCLINUX   0000      (1) Ceres              20161108\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-94188d3cd923>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     11\u001b[0m         \u001b[0mline\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdb\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreadline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m         \u001b[0;32mprint\u001b[0m \u001b[0mline\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m         \u001b[0mephem\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreaddb\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: Bogus: 00001    3.34  0.12 K172G 266.80164   72.90774   80.30986   10.59240  0.0756828  0.21403411   2.7678809  0 MPO389571  6651 113 1801-2016 0.60 M-v 30h MPCLINUX   0000      (1) Ceres              20161108\n"
     ]
    }
   ],
   "source": [
    "import ephem\n",
    "\n",
    "fnm = 'side/MPCORB.DAT'\n",
    "with open(fnm, 'r+') as db:\n",
    "    line = db.readline()\n",
    "    while len(line) > 0:\n",
    "        if line.startswith('---'):\n",
    "            break\n",
    "        line = db.readline()\n",
    "    for i in xrange(1):\n",
    "        line = db.readline()\n",
    "        print line\n",
    "        ephem.readdb(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ -2.27891559e+08]\n",
      " [ -7.42380374e+07]\n",
      " [ -2.78883236e+07]\n",
      " [  7.63362697e+05]\n",
      " [ -1.62758469e+06]\n",
      " [ -7.67220749e+05]]\n"
     ]
    }
   ],
   "source": [
    "import de423\n",
    "from jplephem import Ephemeris\n",
    "#https://pypi.python.org/pypi/jplephem\n",
    "\n",
    "e = Ephemeris(de423)\n",
    "jed = 2444391.5  # 1980.06.01\n",
    "print e.compute('mars', jed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[16, 11, 2015]\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'sv' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-2-4a03d7451c8e>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     96\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mg\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     97\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 98\u001b[1;33m \u001b[0msun_a\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msun_b\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlinpar\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m12\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m13\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     99\u001b[0m \u001b[0mf_sun\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0msun_a\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0msun_b\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    100\u001b[0m \u001b[0mt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m24.\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'sv' is not defined"
     ]
    }
   ],
   "source": [
    "import ephem\n",
    "import numpy as np\n",
    "from utils import parse_coord\n",
    "\n",
    "def dpm(year):\n",
    "    ds = np.array([30]*12)\n",
    "    ds[ :7:2] += 1\n",
    "    ds[7: :2] += 1\n",
    "    if year % 4 == 0:\n",
    "        ds[1] = 29\n",
    "    else:\n",
    "        ds[1] = 28 \n",
    "    return ds\n",
    "\n",
    "\n",
    "def prev_date(date):\n",
    "    ds = dpm(date[2])\n",
    "    \n",
    "    if date[0] == 1:\n",
    "        if date[1] == 1:\n",
    "            nd = [31, 12, date[2] - 1]\n",
    "        else:\n",
    "            nd = [ds[date[1]-2], date[1]-1, date[2]]\n",
    "    else:\n",
    "        nd = [date[0]-1, date[1], date[2]]\n",
    "    return nd\n",
    "\n",
    "def next_date(date):\n",
    "    ds = dpm(date[2])\n",
    "    \n",
    "    if date[0] == ds[date[1]-1]:\n",
    "        if date[1] == 12:\n",
    "            nd = [1, 1, date[2] + 1]\n",
    "        else:\n",
    "            nd = [1, date[1]+1, date[2]]\n",
    "    else:\n",
    "        nd = [date[0]+1, date[1], date[2]]\n",
    "    return nd\n",
    "\n",
    "def ephsun(date):\n",
    "    # Moscow coordinates\n",
    "    (lat, lon) = (55.7522200, 37.6155600) \n",
    "\n",
    "    date2 = prev_date(date)\n",
    "    \n",
    "    time = '21:00' #\n",
    "    sdate1 = str(date2[2]) + '/' + str(date2[1]) + '/' + str(date2[0]) + ' '\n",
    "    sdate2 = str(date[2]) + '/' + str(date[1]) + '/' + str(date[0]) + ' '\n",
    "    sun = ephem.Sun()\n",
    "    #gatech = ephem.Observer()\n",
    "    #gatech.lon, gatech.lat = str(lon), str(lat)\n",
    "    \n",
    "    gatech = sdate1 + time\n",
    "    #gatech.date = sdate1 + time\n",
    "\n",
    "    sun.compute(gatech)\n",
    "    sun1 = [parse_coord(str(sun.a_ra), 'ra'), parse_coord(str(sun.a_dec), 'deg')]\n",
    "    \n",
    "    gatech = sdate2 + time #.date\n",
    "    #gatech.date = sdate2 + time\n",
    "\n",
    "    sun.compute(gatech)\n",
    "    sun2 = [parse_coord(str(sun.a_ra), 'ra'), parse_coord(str(sun.a_dec), 'deg')]\n",
    "    return sun1, sun2\n",
    "\n",
    "def ephsun_iv(date1, date2):\n",
    "    # Moscow coordinates\n",
    "    (lat, lon) = (55.7522200, 37.6155600) \n",
    "\n",
    "    resra = []\n",
    "    resdec = []\n",
    "    curdate = prev_date(date1)\n",
    "    \n",
    "    time = '21:00' #\n",
    "    sun = ephem.Sun()\n",
    "    while curdate != date2:\n",
    "        sdate = str(curdate[2]) + '/' + str(curdate[1]) + '/' + str(curdate[0]) + ' '\n",
    "        gatech = sdate + time\n",
    "        sun.compute(gatech)\n",
    "        sun1 = [parse_coord(str(sun.a_ra), 'ra'), parse_coord(str(sun.a_dec), 'deg')]\n",
    "        resra.append(sun1[0].f2hd())\n",
    "        resdec.append(sun1[1].f2hd())\n",
    "        curdate = next_date(curdate)\n",
    "    return resra, resdec\n",
    "\n",
    "\n",
    "date1 = [15, 11, 2015]\n",
    "date2 = [1, 1, 2016]\n",
    "print next_date(date1)\n",
    "svra, svdec = ephsun_iv(date1, date2)    \n",
    "\n",
    "def linpar(ra1, ra2):\n",
    "    g = np.array([[ 1., 0.],\n",
    "                  [-1., 1.]])\n",
    "    t = np.array([ra1, ra2]).reshape(-1,1)\n",
    "    return np.dot(g, t)\n",
    "\n",
    "sun_a, sun_b = linpar(sv[12],sv[13])\n",
    "f_sun = lambda x: sun_a + sun_b*x\n",
    "t = np.linspace(0,1,24.)\n",
    "print svdec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Month ready: 11 649\n",
      "Month ready: 12 649\n",
      "File type DAF/SPK and format LTL-IEEE with 18 segments:\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Janus (610)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Epimetheus (611)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Helene (612)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Telesto (613)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Calypso (614)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Atlas (615)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Prometheus (616)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Pandora (617)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Methone (632)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Pallene (633)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Polydeuces (634)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Anthe (649)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Aegaeon (653)\n",
      "2433282.50..2469807.50  Saturn Barycenter (6) -> Saturn (699)\n",
      "2433282.50..2469807.50  Solar System Barycenter (0) -> Earth Barycenter (3)\n",
      "2433282.50..2469807.50  Solar System Barycenter (0) -> Saturn Barycenter (6)\n",
      "2433282.50..2469807.50  Solar System Barycenter (0) -> Sun (10)\n",
      "2433282.50..2469807.50  Earth Barycenter (3) -> Earth (399)\n"
     ]
    }
   ],
   "source": [
    "from jplephem.spk import SPK\n",
    "from coord import coord\n",
    "import numpy as np\n",
    "import ephem\n",
    "\n",
    "def bary2eq(pos):\n",
    "    [a, b, c] = pos\n",
    "    Dec = 180./np.pi*np.arctan(c / np.sqrt( a*a + b*b ))\n",
    "    RA  = 180. + 180./np.pi*np.arctan(b / a)\n",
    "    return RA / 15., Dec\n",
    "\n",
    "def date2jed(date):\n",
    "    st = ''\n",
    "    for i in xrange(len(date)):\n",
    "        x = date[i]\n",
    "        st += str(x)\n",
    "        if i < 2:\n",
    "            st += '/'\n",
    "        if i == 2:\n",
    "            st += ' '\n",
    "        if (i > 2) and (i < 5):\n",
    "            if i != len(date) - 1:\n",
    "                st += ':'\n",
    "    djd = float(ephem.Date(st)) # DJD, dublin jd\n",
    "    return st, djd + 2415020\n",
    "\n",
    "def dmy2(date):\n",
    "    return [date[2], date[1], date[0]]\n",
    "\n",
    "\n",
    "def planetpos2(a, b, jed):\n",
    "    # Sattelite-Saturn\n",
    "    position1, _  = kernel[6, 699].compute_and_differentiate(jed) \n",
    "    # Satellite-Saturn + Saturn-Solar -> Satellite-Sollar\n",
    "    position2, _  = kernel[0,   6].compute_and_differentiate(jed) \n",
    "    # Satellite-Sollar - Earth-Solar -> Satellite-Earth\n",
    "    position3, _  = kernel[0,   3].compute_and_differentiate(jed) \n",
    "    # Satellite-Earth - Earth-Earth -> Fin!\n",
    "    position4, _  = kernel[3, 399].compute_and_differentiate(jed) \n",
    "\n",
    "    pos = position1[:3] + position2 - position3 - position4 #\n",
    "    return pos\n",
    "\n",
    "def planetpos(barc, obj, jed):\n",
    "    # Sattelite-Saturn\n",
    "    position1, _  = kernel[barc, obj].compute_and_differentiate(jed) \n",
    "    # Satellite-Saturn + Saturn-Solar -> Satellite-Sollar\n",
    "    position2, _  = kernel[0,   barc].compute_and_differentiate(jed) \n",
    "    # Satellite-Sollar - Earth-Solar -> Satellite-Earth\n",
    "    position3, _  = kernel[0,   3].compute_and_differentiate(jed) \n",
    "    # Satellite-Earth - Earth-Earth -> Fin!\n",
    "    position4, _  = kernel[3, 399].compute_and_differentiate(jed) \n",
    "\n",
    "    pos = position1[:3] + position2 - position3 - position4 #\n",
    "    return pos\n",
    "\n",
    "def ephplan(barc, obj, kernel, date):\n",
    "    date1 = dmy2(date) + [21, 0, 0]\n",
    "    st, jed = date2jed(date1)#2457061.5 # 8.02.2015 00:00\n",
    "    pos = planetpos(barc, obj, jed)\n",
    "    eqpos = bary2eq(pos)\n",
    "    return eqpos\n",
    "\n",
    "def ephplan_iv(barc, obj, kernel, date1, date2):\n",
    "    # Moscow coordinates\n",
    "    (lat, lon) = (55.7522200, 37.6155600) \n",
    "\n",
    "    resdec = []\n",
    "    resra = []\n",
    "    curdate = prev_date(date1)\n",
    "    while curdate != date2:\n",
    "        rd = ephplan(barc, obj, kernel, curdate)\n",
    "        resra.append(rd[0])\n",
    "        resdec.append(rd[1])\n",
    "        curdate = next_date(curdate)\n",
    "    return resra, resdec\n",
    "\n",
    "def linpar(ra1, ra2):\n",
    "    g = np.array([[ 1., 0.],\n",
    "                  [-1., 1.]])\n",
    "    t = np.array([ra1, ra2]).reshape(-1,1)\n",
    "    par = np.dot(g, t)\n",
    "    f = lambda x: par[0] + par[1]*x\n",
    "    return f\n",
    "\n",
    "dirnm = './saturn/'\n",
    "fnm = 'sat382.bsp'\n",
    "\n",
    "barc = 6\n",
    "obj = 649\n",
    "\n",
    "#date1 = [15, 11, 2015]\n",
    "#date2 = [1, 1, 2016]\n",
    "dm = dpm(year)\n",
    "kernel = SPK.open(dirnm+fnm)\n",
    "t = np.linspace(0,1,24.)\n",
    "\n",
    "year = 2015\n",
    "\n",
    "for month in range(11,13):\n",
    "    fnm = './csv/' + str(month)+str(year)+'.csv'\n",
    "    f = open(fnm, 'w+')\n",
    "    f.write('Date,Time,Dist(RA),Dist(DEC),RA,DEC\\n')\n",
    "\n",
    "    date1 = [1, month, year]\n",
    "    date2 = [dm[month-1], month, year]\n",
    "    date2 = next_date(date2)\n",
    "    svra, svdec = ephsun_iv(date1, date2)    \n",
    "    pvra, pvdec = ephplan_iv(barc, obj, kernel, date1, date2)\n",
    "    \n",
    "    for i in xrange(len(pvra)-1):\n",
    "        f_pl = linpar(pvra[i], pvra[i+1])\n",
    "        f_pl_dec = linpar(pvdec[i], pvdec[i+1])\n",
    "        f_sun = linpar(svra[i],svra[i+1])\n",
    "        f_sun_dec = linpar(svdec[i],svdec[i+1])\n",
    "        str_date = str(i+1) + '/' + str(month) + '/' + str(year)\n",
    "        for k in xrange(t.size):\n",
    "            tk = t[k]\n",
    "            dra = f_sun(tk) - f_pl(tk)\n",
    "            ddec = f_sun_dec(tk) - f_pl_dec(tk)\n",
    "            #sdra = coord(0, 0, dra[0]*60*60, np.sign(dra), 'ra')\n",
    "            #sddec = coord(0, 0, ddec[0]*60*60, np.sign(ddec), 'deg')\n",
    "            sra_pl = coord(0, 0, f_pl(tk)[0]*60*60, np.sign(f_pl(tk)[0]), 'ra')\n",
    "            sdec_pl = coord(0, 0, f_pl_dec(tk)[0]*60*60, np.sign(f_pl_dec(tk)), 'deg')                \n",
    "            sra_sun = coord(0, 0, f_sun(tk)[0]*60*60, np.sign(f_sun(tk)[0]), 'ra')\n",
    "            sdec_sun = coord(0, 0, f_sun_dec(tk)[0]*60*60, np.sign(f_sun_dec(tk)), 'deg')   \n",
    "            #dra = sra_sun - sra_pl\n",
    "            #ddec = sdec_sun - sdec_pl\n",
    "            #dra = dra.f2hd()\n",
    "            #dra = dra[0] % 12\n",
    "            #ddec = ddec[0] % 180\n",
    "            dra = dra[0]\n",
    "            ddec = ddec[0]\n",
    "            strtime = str(k) + ':' + '00'\n",
    "            f.write(str_date + ',' + strtime + ',' + str(dra) + ',' + str(ddec) + ',' +  str(sra) + ',' +str(sdec)+'\\n')\n",
    "        f.write('\\n')\n",
    "    f.close()\n",
    "    print \"Month ready:\", month, obj\n",
    "\n",
    "print kernel\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File type DAF/SPK and format LTL-IEEE with 300 segments:\n",
      "2378495.00..2524575.00  Sun (10) -> Ceres (2000001)\n",
      "2378495.00..2524575.00  Sun (10) -> Pallas (2000002)\n",
      "2378495.00..2524575.00  Sun (10) -> Vesta (2000004)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000029)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000016)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000015)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000010)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000007)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000006)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000003)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000345)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000012)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000442)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000027)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000287)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000043)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000084)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000030)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000536)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000163)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000105)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000554)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000313)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000115)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000230)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000337)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000080)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000009)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000063)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000025)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000192)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000304)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000020)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000566)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000654)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000135)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000618)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000635)\n",
      "2378495.00..2524575.00  Sun (10) -> Lutetia (2000021)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000042)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000019)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000011)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000168)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000017)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2001467)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000329)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000046)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000751)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000089)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000449)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000334)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000134)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000535)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000626)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000005)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000712)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000107)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000409)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000362)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000405)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000032)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000014)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000091)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000111)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2001015)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000404)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000344)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000849)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000389)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000347)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000070)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000194)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000053)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000078)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000407)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000454)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000023)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000124)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000164)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000037)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000072)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000224)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000369)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000050)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000476)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000713)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000085)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000026)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000455)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000144)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000488)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000466)\n",
      "2378495.00..2524575.00  Sun (10) -> Davida (2000511)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000240)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2001036)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000141)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000097)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000326)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000077)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000694)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000145)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000075)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000225)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000201)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000324)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000034)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000098)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000031)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000175)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000109)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000058)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000103)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000059)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000054)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000146)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000045)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000209)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000410)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000160)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000140)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000156)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000110)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000187)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000018)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000387)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000200)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000185)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000057)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000247)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000372)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000980)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000481)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000521)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000206)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000038)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000173)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000036)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000171)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000490)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000381)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000128)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000093)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000071)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000356)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000127)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000041)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000143)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000094)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000039)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000088)\n",
      "2378495.00..2524575.00  Sun (10) -> Eros (2000433)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000008)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000532)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000092)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000275)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000065)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000444)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000393)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000147)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000205)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000028)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000139)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000068)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000322)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000416)\n",
      "2378495.00..2524575.00  Sun (10) -> Kleopatra (2000216)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000354)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000346)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000236)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000365)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000266)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000804)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000385)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000081)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000909)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000129)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000508)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000195)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000047)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000104)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000426)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000471)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000431)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000386)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000238)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000338)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000022)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000024)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000674)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000702)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000040)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000705)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000776)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000375)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000069)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000150)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000117)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000035)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000747)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000772)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000360)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000388)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000250)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000176)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000162)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000350)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000506)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000211)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000514)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000241)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000740)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000096)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000663)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000451)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000895)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000423)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000095)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000049)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000602)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000052)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000159)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000357)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000545)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000469)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000373)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000048)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000137)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000086)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000121)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000196)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000062)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000090)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000791)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000120)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000683)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000212)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000283)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000691)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000596)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000709)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000191)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000074)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000203)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000377)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000112)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000165)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000363)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000276)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000739)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2001093)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000268)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2001021)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000210)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000505)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000099)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000690)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000366)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000814)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000788)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000303)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000780)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000259)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000489)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000568)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000358)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000424)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000412)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000762)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000308)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000102)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000056)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000154)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000083)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000328)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000769)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000181)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000221)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000491)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000349)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000148)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000704)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000595)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000213)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000419)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000335)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000114)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000233)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000498)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000914)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000044)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000790)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000076)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000420)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000106)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000130)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000013)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000087)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000051)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000773)\n",
      "2378495.00..2524575.00  Sun (10) -> Unknown Target (2000336)\n"
     ]
    }
   ],
   "source": [
    "from jplephem.spk import SPK\n",
    "from coord import coord\n",
    "import numpy as np\n",
    "import ephem\n",
    "\n",
    "dirnm = './'\n",
    "fnm = 'codes_300ast_20100725.bsp'\n",
    "fnm_planet = 'de432s.bsp'\n",
    "year = 2015\n",
    "\n",
    "dm = dpm(year)\n",
    "kernel = SPK.open(dirnm+fnm)\n",
    "kernel.planet = SPK\n",
    "t = np.linspace(0,1,24.)\n",
    "\n",
    "print kernel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.728155339806\n",
      "0.970873786408\n",
      "1.21359223301\n",
      "1.69902912621\n",
      "2.42718446602\n",
      "2.66990291262\n",
      "2.91262135922\n",
      "3.39805825243\n",
      "3.64077669903\n",
      "4.12621359223\n",
      "4.36893203883\n",
      "4.61165048544\n",
      "4.85436893204\n",
      "5.09708737864\n",
      "5.33980582524\n",
      "6.06796116505\n",
      "6.31067961165\n",
      "6.79611650485\n",
      "7.28155339806\n",
      "7.52427184466\n",
      "8.00970873786\n",
      "8.25242718447\n",
      "8.49514563107\n",
      "8.73786407767\n",
      "8.98058252427\n",
      "9.22330097087\n",
      "9.46601941748\n",
      "9.70873786408\n",
      "9.95145631068\n",
      "10.1941747573\n",
      "10.4368932039\n",
      "10.6796116505\n",
      "10.9223300971\n",
      "11.4077669903\n",
      "11.6504854369\n",
      "11.8932038835\n",
      "12.1359223301\n",
      "12.6213592233\n",
      "12.8640776699\n",
      "13.1067961165\n",
      "13.3495145631\n",
      "13.5922330097\n",
      "14.0776699029\n",
      "14.3203883495\n",
      "14.5631067961\n",
      "15.0485436893\n",
      "15.5339805825\n",
      "15.7766990291\n",
      "16.0194174757\n",
      "16.2621359223\n",
      "16.5048543689\n",
      "16.7475728155\n",
      "16.9902912621\n",
      "17.2330097087\n",
      "17.4757281553\n",
      "17.7184466019\n",
      "17.9611650485\n",
      "18.4466019417\n",
      "18.6893203883\n",
      "18.932038835\n",
      "19.1747572816\n",
      "19.4174757282\n",
      "19.6601941748\n",
      "19.9029126214\n",
      "20.145631068\n",
      "20.3883495146\n",
      "20.6310679612\n",
      "21.1165048544\n",
      "21.359223301\n",
      "21.6019417476\n",
      "21.8446601942\n",
      "22.3300970874\n",
      "22.572815534\n",
      "22.8155339806\n",
      "23.0582524272\n",
      "23.3009708738\n",
      "23.5436893204\n",
      "23.786407767\n",
      "24.0291262136\n",
      "24.2718446602\n",
      "24.7572815534\n",
      "25.0\n",
      "25.2427184466\n",
      "25.4854368932\n",
      "25.7281553398\n",
      "25.9708737864\n",
      "26.213592233\n",
      "26.4563106796\n",
      "27.1844660194\n",
      "27.427184466\n",
      "27.6699029126\n",
      "27.9126213592\n",
      "28.1553398058\n",
      "28.3980582524\n",
      "28.640776699\n",
      "29.1262135922\n",
      "29.3689320388\n",
      "29.854368932\n",
      "30.0970873786\n",
      "30.3398058252\n",
      "30.5825242718\n",
      "31.3106796117\n",
      "31.5533980583\n",
      "32.0388349515\n",
      "32.2815533981\n",
      "32.5242718447\n",
      "32.7669902913\n",
      "33.0097087379\n",
      "33.2524271845\n",
      "34.2233009709\n",
      "34.4660194175\n",
      "34.7087378641\n",
      "34.9514563107\n",
      "35.1941747573\n",
      "35.4368932039\n",
      "35.9223300971\n",
      "36.1650485437\n",
      "36.4077669903\n",
      "36.6504854369\n",
      "36.8932038835\n",
      "37.1359223301\n",
      "37.3786407767\n",
      "37.6213592233\n",
      "37.8640776699\n",
      "38.1067961165\n",
      "38.3495145631\n",
      "38.5922330097\n",
      "38.8349514563\n",
      "39.0776699029\n",
      "39.5631067961\n",
      "39.8058252427\n",
      "40.0485436893\n",
      "40.7766990291\n",
      "41.0194174757\n",
      "41.2621359223\n",
      "41.7475728155\n",
      "42.2330097087\n",
      "42.7184466019\n",
      "42.9611650485\n",
      "43.4466019417\n",
      "43.6893203883\n",
      "43.932038835\n",
      "44.4174757282\n",
      "44.6601941748\n",
      "45.145631068\n",
      "45.6310679612\n",
      "45.8737864078\n",
      "46.1165048544\n",
      "46.359223301\n",
      "46.8446601942\n",
      "47.0873786408\n",
      "47.3300970874\n",
      "47.572815534\n",
      "47.8155339806\n",
      "48.5436893204\n",
      "48.786407767\n",
      "49.0291262136\n",
      "49.2718446602\n",
      "49.5145631068\n",
      "49.7572815534\n",
      "50.0\n",
      "50.2427184466\n",
      "50.4854368932\n",
      "50.7281553398\n",
      "50.9708737864\n",
      "51.213592233\n",
      "51.4563106796\n",
      "51.6990291262\n",
      "52.1844660194\n",
      "52.427184466\n",
      "52.6699029126\n",
      "52.9126213592\n",
      "53.1553398058\n",
      "53.640776699\n",
      "53.8834951456\n",
      "54.1262135922\n",
      "54.3689320388\n",
      "54.6116504854\n",
      "55.3398058252\n",
      "55.8252427184\n",
      "56.067961165\n",
      "56.3106796117\n",
      "56.5533980583\n",
      "56.7961165049\n",
      "57.2815533981\n",
      "57.5242718447\n",
      "57.7669902913\n",
      "58.0097087379\n",
      "58.2524271845\n",
      "58.4951456311\n",
      "58.7378640777\n",
      "59.2233009709\n",
      "59.7087378641\n",
      "59.9514563107\n",
      "60.4368932039\n",
      "60.6796116505\n",
      "60.9223300971\n",
      "61.1650485437\n",
      "61.6504854369\n",
      "61.8932038835\n",
      "62.1359223301\n",
      "62.3786407767\n",
      "62.6213592233\n",
      "62.8640776699\n",
      "63.1067961165\n",
      "63.3495145631\n",
      "63.5922330097\n",
      "64.0776699029\n",
      "64.3203883495\n",
      "64.5631067961\n",
      "65.0485436893\n",
      "65.2912621359\n",
      "65.5339805825\n",
      "65.7766990291\n",
      "66.0194174757\n",
      "66.2621359223\n",
      "66.9902912621\n",
      "67.2330097087\n",
      "67.4757281553\n",
      "67.7184466019\n",
      "67.9611650485\n",
      "68.2038834951\n",
      "68.4466019417\n",
      "68.6893203883\n",
      "68.932038835\n",
      "69.1747572816\n",
      "69.6601941748\n",
      "69.9029126214\n",
      "70.3883495146\n",
      "70.6310679612\n",
      "70.8737864078\n",
      "71.1165048544\n",
      "71.6019417476\n",
      "71.8446601942\n",
      "72.0873786408\n",
      "72.3300970874\n",
      "72.8155339806\n",
      "73.3009708738\n",
      "73.786407767\n",
      "74.2718446602\n",
      "74.5145631068\n",
      "75.2427184466\n",
      "75.4854368932\n",
      "75.7281553398\n",
      "75.9708737864\n",
      "76.213592233\n",
      "76.4563106796\n",
      "76.9417475728\n",
      "77.1844660194\n",
      "77.6699029126\n",
      "77.9126213592\n",
      "78.1553398058\n",
      "78.640776699\n",
      "79.1262135922\n",
      "79.3689320388\n",
      "79.854368932\n",
      "80.0970873786\n",
      "80.3398058252\n",
      "80.8252427184\n",
      "81.067961165\n",
      "81.3106796117\n",
      "81.5533980583\n",
      "81.7961165049\n",
      "82.0388349515\n",
      "82.2815533981\n",
      "82.5242718447\n",
      "82.7669902913\n",
      "83.0097087379\n",
      "83.2524271845\n",
      "83.4951456311\n",
      "83.7378640777\n",
      "83.9805825243\n",
      "84.4660194175\n",
      "84.7087378641\n",
      "85.1941747573\n",
      "85.4368932039\n",
      "85.6796116505\n",
      "85.9223300971\n",
      "86.1650485437\n",
      "86.4077669903\n",
      "86.6504854369\n",
      "86.8932038835\n",
      "87.1359223301\n",
      "87.3786407767\n",
      "87.6213592233\n",
      "87.8640776699\n",
      "88.5922330097\n",
      "88.8349514563\n",
      "89.0776699029\n",
      "89.3203883495\n",
      "89.5631067961\n",
      "89.8058252427\n",
      "90.0485436893\n",
      "90.5339805825\n",
      "91.2621359223\n",
      "91.5048543689\n",
      "91.7475728155\n",
      "91.9902912621\n",
      "92.2330097087\n",
      "92.7184466019\n",
      "92.9611650485\n",
      "93.2038834951\n",
      "93.4466019417\n",
      "94.1747572816\n",
      "94.4174757282\n",
      "94.6601941748\n",
      "95.145631068\n",
      "95.3883495146\n",
      "95.8737864078\n",
      "96.1165048544\n",
      "96.359223301\n",
      "96.6019417476\n",
      "96.8446601942\n",
      "97.0873786408\n",
      "97.3300970874\n",
      "97.572815534\n",
      "97.8155339806\n",
      "98.0582524272\n",
      "98.3009708738\n",
      "98.5436893204\n",
      "98.786407767\n",
      "99.0291262136\n",
      "99.2718446602\n",
      "Month  1 ready\n",
      "Month  2 ready\n",
      "Month  3 ready\n",
      "Month  4 ready\n",
      "Month  5 ready\n",
      "Month  6 ready\n",
      "Month  7 ready\n",
      "Month  8 ready\n",
      "Month  9 ready\n",
      "Month  10 ready\n",
      "Month  11 ready\n",
      "Month  12 ready\n",
      "Check it!\n"
     ]
    }
   ],
   "source": [
    "import ephem\n",
    "from coord import coord\n",
    "\n",
    "import numpy as np\n",
    "import sys, requests\n",
    "import math\n",
    "from scipy.optimize import brentq\n",
    "\n",
    "def dpm(year):\n",
    "    ds = np.array([30]*12)\n",
    "    ds[ :7:2] += 1\n",
    "    ds[7: :2] += 1\n",
    "    if year % 4 == 0:\n",
    "        ds[1] = 29\n",
    "    else:\n",
    "        ds[1] = 28 \n",
    "    return ds\n",
    "\n",
    "def next_date(date):\n",
    "    ds = dpm(date[2])\n",
    "    \n",
    "    if date[0] == ds[date[1]-1]:\n",
    "        if date[1] == 12:\n",
    "            nd = [1, 1, date[2] + 1]\n",
    "        else:\n",
    "            nd = [1, date[1]+1, date[2]]\n",
    "    else:\n",
    "        nd = [date[0]+1, date[1], date[2]]\n",
    "    return nd\n",
    "\n",
    "\n",
    "def ephsun_year(year):\n",
    "    # Moscow coordinates\n",
    "    (lat, lon) = (55.7522200, 37.6155600) \n",
    "\n",
    "    time = '21:00' #\n",
    "    days_per_month = dpm(year)\n",
    "    gatech = str(year-1) + '/12/31 ' + time\n",
    "    sun = ephem.Sun()\n",
    "    sun.compute(gatech)\n",
    "    rv = [[gatech, sun.a_ra, sun.a_dec]]\n",
    "    \n",
    "    for i, x in enumerate(days_per_month):\n",
    "        for j in xrange(x):\n",
    "            gatech = (str(year) + '/' + str(i+1) + '/' + str(j+1) + ' ' + time)\n",
    "            sun.compute(gatech)\n",
    "            rv.append([gatech, sun.a_ra, sun.a_dec])\n",
    "    return rv\n",
    "\n",
    "def mpc_data(r):\n",
    "    r = r.split(',')\n",
    "    params = {}\n",
    "    #if num >= uplim:\n",
    "    #    d = r[0].replace(' ', '+')\n",
    "    #    params['designation'] = d\n",
    "    #else:\n",
    "    d = r[0].split(' ')\n",
    "    params['number'] = d[0]\n",
    "    params['json'] = 1\n",
    "    \n",
    "    req = requests.post(url, params, auth = ('mpc_ws', 'mpc!!ws'))\n",
    "    y = req.text\n",
    "    if y == '[]':\n",
    "        return ''\n",
    "    \n",
    "    y = y.replace('{', '')\n",
    "    y = y.replace('}', '')\n",
    "    y = y.replace('[', '')\n",
    "    y = y.replace(']', '')\n",
    "    y = y.replace('\"','')\n",
    "    y = y.split(',')\n",
    "\n",
    "    res = {}\n",
    "    for x in y:\n",
    "        z = x.split(':', 1)\n",
    "        res[z[0]] = z[1]\n",
    "\n",
    "    diam = res['diameter_neowise']\n",
    "    if diam == 'null':\n",
    "        diam = ''\n",
    "    #name = res['designation']\n",
    "    #while name.startswith('0'):\n",
    "    #    name = name[1:]\n",
    "    #if num < uplim:\n",
    "    #name = '(' + name + ') '\n",
    "    #tmp = res['name']\n",
    "    #if tmp == 'null':\n",
    "    #    tmp = res['']\n",
    "    #name = name + res['name']\n",
    "    return diam\n",
    "\n",
    "year = 2016\n",
    "DECLIM = 16.\n",
    "\n",
    "days_per_month = dpm(year)\n",
    "result = [[]]*len(days_per_month)\n",
    "\n",
    "sun = ephem.Sun()\n",
    "sun_eph = ephsun_year(year)\n",
    "#print sun_eph\n",
    "\n",
    "dbfnm = 'large.edb'#'MPCORB.DAT.edb'\n",
    "url = 'http://mpc.cfa.harvard.edu/ws/search'\n",
    "\n",
    "max_entry = 412#455144\n",
    "save_int = 100\n",
    "offset = 0\n",
    "\n",
    "with open(dbfnm, 'r+') as f:\n",
    "    r = '#'\n",
    "    while r.startswith('#'):\n",
    "        r = f.readline()\n",
    "    for idx in xrange(max_entry):\n",
    "        if idx > 0:\n",
    "            r = f.readline()\n",
    "        else:\n",
    "            for p in xrange(offset):\n",
    "                r = f.readline()\n",
    "        yh = ephem.readdb(r)\n",
    "        for i in xrange(len(sun_eph) - 1):\n",
    "            if idx % save_int == 0:\n",
    "                np.savez('asteroids', result)\n",
    "            [gatech1, sun1_ra, sun1_dec] = sun_eph[i]\n",
    "            yh.compute(gatech1)\n",
    "            dif_ra1 = yh.a_ra - sun1_ra\n",
    "\n",
    "            [gatech2, sun2_ra, sun2_dec] = sun_eph[i+1]\n",
    "            yh.compute(gatech2)\n",
    "            dif_ra2 = yh.a_ra - sun2_ra\n",
    "\n",
    "            if (dif_ra1 * dif_ra2) < 0:\n",
    "                a = ephem.date(gatech1)\n",
    "                b = ephem.date(gatech2)\n",
    "                def func(x):\n",
    "                    sun.compute(x)\n",
    "                    yh.compute(x)\n",
    "                    dist = yh.a_ra - sun.a_ra\n",
    "                    dist *= 180./np.pi\n",
    "                    return dist\n",
    "                t0 = brentq(func, a, b)\n",
    "                sun.compute(t0)\n",
    "                yh.compute(t0)\n",
    "                \n",
    "                if abs(sun.a_ra - yh.a_ra)* 180./np.pi > 1e-3:\n",
    "                    continue\n",
    "                \n",
    "                date = gatech1.split(' ')[0]\n",
    "                date = date.split('/')\n",
    "                \n",
    "                date = [int(date[2]), int(date[1]), int(date[0])]\n",
    "                date = next_date(date)\n",
    "                date = [str(x) for x in date]\n",
    "                \n",
    "                ind_month = int(date[1]) - 1\n",
    "                if len(date[1]) == 1:\n",
    "                    date[1] = ' ' + date[1]\n",
    "                if len(date[0]) == 1:\n",
    "                    date[0] = ' ' + date[0]\n",
    "                date = date[0] + '/' + date[1] + '/' + date[2]\n",
    "                time = coord(0, 0, (t0 - a)*24.*60.*60., 1, 'ra')\n",
    "                diam = mpc_data(r)\n",
    "                q = r.split(',')\n",
    "                name = q[0]\n",
    "                \n",
    "                yh_ra = str(yh.a_ra)\n",
    "                yh_dec = str(yh.a_dec)\n",
    "                dist = ephem.degrees(yh.a_dec - sun.a_dec) \n",
    "                if abs(dist) * 180. / np.pi > DECLIM:\n",
    "                    continue\n",
    "                print (idx / float(max_entry))*100\n",
    "\n",
    "                #print date, sun.a_ra, yh.a_ra, sun.a_ra - yh.a_ra\n",
    "\n",
    "                dist = str(dist)\n",
    "                \n",
    "                output = date + ',' + str(time) + ',' + name + ',' + diam + ',' + dist + ',' + yh_ra + ',' + yh_dec\n",
    "                result[ind_month] = result[ind_month] + [[t0, output]]\n",
    "np.savez('asteroids', result)\n",
    "\n",
    "dirnm = './csv/'\n",
    "\n",
    "head = 'Date,Time,Name,Diameter,Distance,Ra,Dec\\n'\n",
    "for k in xrange(len(result)):\n",
    "    f = open(dirnm + str(k+1)+'.csv', 'w+')\n",
    "    res = sorted(result[k], key = lambda x : x[0])\n",
    "    f.writelines(head)\n",
    "    for p, g in res:\n",
    "        f.writelines(g + '\\n')\n",
    "    f.close()\n",
    "    print 'Month ', k+1, 'ready'\n",
    "print 'Check it!'\n",
    "            \n",
    "\n",
    "\n",
    "#gatech = ephem.Observer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.70546926818 269:36:12.7 269.603529695\n",
      "11730\n",
      "Month  1 ready\n",
      "Month  2 ready\n",
      "Month  3 ready\n",
      "Month  4 ready\n",
      "Month  5 ready\n",
      "Month  6 ready\n",
      "Month  7 ready\n",
      "Month  8 ready\n",
      "Month  9 ready\n",
      "Month  10 ready\n",
      "Month  11 ready\n",
      "Month  12 ready\n",
      "Check it!\n"
     ]
    }
   ],
   "source": [
    "q  = sun.a_ra + 0.2\n",
    "print q, ephem.degrees(q), q*180./np.pi\n",
    "print idx\n",
    "\n",
    "dirnm = './csv/'\n",
    "\n",
    "head = 'Date,Time,Name,Diameter,Distance,Ra,Dec\\n'\n",
    "for k in xrange(len(result)):\n",
    "    f = open(dirnm + str(k+1)+'.csv', 'w+')\n",
    "    res = sorted(result[k], key = lambda x : x[0])\n",
    "    f.writelines(head)\n",
    "    for p, g in res:\n",
    "        f.writelines(g + '\\n')\n",
    "    f.close()\n",
    "    print 'Month ', k+1, 'ready'\n",
    "print 'Check it!'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = []\n",
    "y == []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2008 YA,e,31.20091,119.32641,201.3196,2.651856,0,0.8316239,237.1706,1/13/2016,2000.0,18,0.15\n",
      "kkk\n",
      "[u'']\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-151-994d4426bad8>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     63\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdiam\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     64\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 65\u001b[1;33m \u001b[0mrv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmpc_data\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mI\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     66\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mrv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     67\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-151-994d4426bad8>\u001b[0m in \u001b[0;36mmpc_data\u001b[1;34m(r, num, uplim)\u001b[0m\n\u001b[0;32m     50\u001b[0m     \u001b[1;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[1;32min\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     51\u001b[0m         \u001b[0mz\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msplit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m':'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 52\u001b[1;33m         \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mz\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mz\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     53\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     54\u001b[0m     \u001b[0mdiam\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'diameter_neowise'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "import ephem\n",
    "import numpy as np\n",
    "import sys, requests\n",
    "import math\n",
    "\n",
    "def sun_eph():\n",
    "    pass\n",
    "\n",
    "\n",
    "dbfnm = 'MPCORB.DAT.edb'\n",
    "url = 'http://mpc.cfa.harvard.edu/ws/search'\n",
    "\n",
    "f = open(dbfnm, 'r+')\n",
    "\n",
    "r = '#'\n",
    "while r.startswith('#'):\n",
    "    r = f.readline()\n",
    "    #print \n",
    "    \n",
    "for i in xrange(I):\n",
    "    r = f.readline()\n",
    "print r, 'kkk'\n",
    "#while '  ' in r:\n",
    "#    r = r.replace('  ', ' ')\n",
    "#r = r.replace(' ', ',')\n",
    "    \n",
    "\n",
    "\n",
    "rv = mpc_data(r, I)\n",
    "print rv\n",
    "\n",
    "yh = ephem.readdb(r)\n",
    "dt = '2015/09/28 22:10'\n",
    "yh.compute(dt)\n",
    "print(\"%s %s\" % (yh.a_ra, yh.a_dec))\n",
    "print(\"%s %s\" % (ephem.constellation(yh), yh.mag))\n",
    "\n",
    "s = ephem.Sun()\n",
    "s.compute(dt)\n",
    "print(\"%s %s\" % (s.a_ra, s.a_dec))\n",
    "print(\"%s %s\" % (ephem.constellation(s), s.mag))\n",
    "\n",
    "def estimate_diameter(obj):\n",
    "    angle = obj.size / 3600. * np.pi/180.\n",
    "    distance = obj.earth_distance * ephem.meters_per_au / 1000.\n",
    "    return 2 * distance * np.sin(angle / 2.)\n",
    "\n",
    "print estimate_diameter(yh)\n",
    "print ephem.degrees(yh.a_dec - s.a_dec) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = 'fff'\n",
    "a.split?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "''.split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 []\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'o' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-25154905f569>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     50\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtext\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtext\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     51\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 52\u001b[1;33m \u001b[0mtu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mo\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     53\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mobj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mr\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     54\u001b[0m     \u001b[0mN\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mobj\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'ascending_node'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'o' is not defined"
     ]
    }
   ],
   "source": [
    "import ephem\n",
    "\n",
    "\n",
    "fnm = 'MPCORB.DAT'\n",
    "\n",
    "if False:\n",
    "    f = open(fnm, 'r+')\n",
    "\n",
    "    r = ''\n",
    "    while not r.startswith('-'):\n",
    "        r = f.readline()\n",
    "        #print r\n",
    "    for i in xrange(180000):\n",
    "        r = f.readline()\n",
    "    #while '  ' in r:\n",
    "    #    r = r.replace('  ', ' ')\n",
    "    #r = r.replace(' ', ',')\n",
    "    r = r.split()\n",
    "    for rr in r:\n",
    "        print rr\n",
    "\n",
    "    #nb = ephem.readdb('SRPTNGPALO Object,f|M|sp,'\n",
    "    #                  +str(ephem.hours(ephem.degrees(str(coord.RA))))+\n",
    "    #                  ','\n",
    "    #                  +str(ephem.degrees(str(coord.DEC)))\n",
    "    #                  +',0,\n",
    "    #                  '+'2000')\n",
    "    yh = ephem.readdb('C/2002 Y1 (Juels-Holvorcem),e,,,,,,,,,,,')\n",
    "    yh.compute('2003/4/11')     \n",
    "    print(\"%s %s\" % (yh.ra, yh.dec))\n",
    "    print(\"%s %s\" % (ephem.constellation(yh), yh.mag))\n",
    "\n",
    "\n",
    "import sys, requests\n",
    "import math\n",
    "\n",
    "url = 'http://mpc.cfa.harvard.edu/ws/search'\n",
    "\n",
    "params = {}\n",
    "#params['diameter_neowise_min'] = 20.\n",
    "params['json'] = 1\n",
    "#params['inclination_max']= 1.0\n",
    "\n",
    "params['name'] = '1995%20UL40'\n",
    "#params['number'] = 9160\n",
    "\n",
    "\n",
    "r = requests.post(url, params, auth = ('mpc_ws', 'mpc!!ws'))\n",
    "\n",
    "print len(r.text), r.text\n",
    "\n",
    "tu = o\n",
    "for obj in r:\n",
    "    N = obj['ascending_node']\n",
    "    i = obg['inclination']\n",
    "    w = obj['argument_of_perihelion']\n",
    "    a = obj['semimajor_axis']\n",
    "    e = obj['eccentricity']\n",
    "    M = obj['mean_anomaly']\n",
    "    \n",
    "    w1 = N + w # longitude of perihelion\n",
    "    L = M + w1 # mean longitude\n",
    "    q = a*(1-e) # perihelion distance\n",
    "    Q =a*(1+e) # aphelion distance\n",
    "    P = a**1.5 # orbitalperiod\n",
    "    T = (M / 360.) / P # time of perihelion\n",
    "    #v = \n",
    "    #e\n",
    "    \n",
    "    # eccentric anomaly\n",
    "    E = M + e*(180/np.pi) * np.sin(M) * ( 1.0 + e * np.cos(M) )\n",
    "    if e > 0.05:\n",
    "        E0 = E\n",
    "        E1 = E0*10.\n",
    "        while abs(E0 - E1)>1e-5:\n",
    "            E0 = E1\n",
    "            E1 = E0 - ( E0 - e*(180/np.pi) * np.sin(E0) - M ) / ( 1 - e * np.cos(E0) )\n",
    "        E = E1    \n",
    "    xv = a * ( np.cos(E) - e )\n",
    "    yv = a * ( np.sqrt(1.0 - e*e) * np.sin(E) )\n",
    "\n",
    "    v = math.atan2( yv, xv )\n",
    "    r = np.sqrt( xv*xv + yv*yv )\n",
    "    \n",
    "    # heliocentric\n",
    "    xh = r * ( np.cos(N) * np.cos(v+w) - np.sin(N) * np.sin(v+w) * np.cos(i) )\n",
    "    yh = r * ( np.sin(N) * np.cos(v+w) + np.cos(N) * np.sin(v+w) * np.cos(i) )\n",
    "    zh = r * ( np.sin(v+w) * np.sin(i) )\n",
    "        \n",
    "    \n",
    "########### Primary orbital elements\n",
    "    #N = longitude of the ascending node\n",
    "    #i = inclination to the ecliptic (plane of the Earth's orbit)\n",
    "    #w = argument of perihelion\n",
    "    #a = semi-major axis, or mean distance from Sun\n",
    "    #e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)\n",
    "    #M = mean anomaly (0 at perihelion; increases uniformly with time)\n",
    "    \n",
    "    \n",
    "######### Related orbital elements\n",
    "    #w1 = N + w   = longitude of perihelion\n",
    "    #L  = M + w1  = mean longitude\n",
    "    #q  = a*(1-e) = perihelion distance\n",
    "    #Q  = a*(1+e) = aphelion distance\n",
    "    #P  = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)\n",
    "    #T  = Epoch_of_M - (M(deg)/360_deg) / P  = time of perihelion\n",
    "    #v  = true anomaly (angle between position and perihelion)\n",
    "    #E  = eccentric anomaly\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "{\"absolute_magnitude\":\"4.13\",\n",
    " \"albedo_neowise\":0.142,\n",
    " \"albedo_neowise_unc\":0.046,\n",
    " \"albedo_updated_at\":\"2014-07-23T19:43:54Z\",\n",
    " \"aphelion_distance\":\"3.413\",\n",
    " \"apparition_end_date\":null,\n",
    " \"apparition_end_magn\":null,\n",
    " \"apparition_start_date\":null,\n",
    " \"apparition_start_magn\":null,\n",
    " \"arc_length\":70495,\n",
    " \"argument_of_perihelion\":\"309.92586\",\n",
    " \"ascending_node\":\"173.0968983\",\n",
    " \"b_minus_v\":0.635,\n",
    " \"b_minus_v_source\":\"M\",\n",
    " \"binary_object\":null,\n",
    " \"color_ignore\":\"???\",\n",
    " \"color_updated_at\":\"2014-07-23T19:44:17Z\",\n",
    " \"critical_list_numbered_object\":false,\n",
    " \"delta_v\":14.2,\"designation\":\"0000002\",\n",
    " \"diameter_neowise\":544.0,\n",
    " \"diameter_neowise_unc\":60.71,\n",
    " \"earth_moid\":null,\n",
    " \"eccentricity\":\"0.2313681\",\n",
    " \"epoch\":\"2014-05-23.0\",\n",
    " \"epoch_jd\":\"2456800.5\",\n",
    " \"first_observation_date_used\":\"1821-05-16.0\",\n",
    " \"first_opposition_used\":\"1821\",\n",
    " \"g_neowise\":0.11,\n",
    " \"greatest_elong\":null,\n",
    " \"greatest_elong_date\":null,\n",
    " \"greatest_elong_decl\":null,\n",
    " \"greatest_elong_magn\":null,\n",
    " \"h_neowise\":4.06,\n",
    " \"inclination\":\"34.84103\",\n",
    " \"jupiter_moid\":null,\n",
    " \"km_neo\":false,\n",
    " \"last_observation_date_used\":\"2014-05-20.0\",\n",
    " \"last_opposition_used\":\"2014\",\n",
    " \"lightcurve_notes\":\".\",\n",
    " \"lightcurve_quality\":\"3\",\n",
    " \"mars_moid\":null,\n",
    " \"mean_anomaly\":\"35.50312\",\n",
    " \"mean_daily_motion\":\"0.2135891\",\n",
    " \"mercury_moid\":null,\n",
    " \"name\":\"Pallas\",\n",
    " \"neo\":false,\n",
    " \"number\":2,\n",
    " \"observations\":7714,\n",
    " \"oppositions\":103,\n",
    " \"orbit_type\":0,\n",
    " \"orbit_uncertainty\":\"0\",\n",
    " \"orbit_updated_at\":\"2014-07-10T16:08:38Z\",\n",
    " \"p_vector_x\":\"-0.56149472\",\n",
    " \"p_vector_y\":\"0.81832779\",\n",
    " \"p_vector_z\":\"-0.1227327\",\n",
    " \"panstarrs_v_minus_gprime\":-0.35,\n",
    " \"panstarrs_v_minus_gprime_source\":\"P\",\n",
    " \"panstarrs_v_minus_iprime\":0.32,\n",
    " \"panstarrs_v_minus_iprime_source\":\"P\",\n",
    " \"panstarrs_v_minus_rprime\":0.14,\n",
    " \"panstarrs_v_minus_rprime_source\":\"P\",\n",
    " \"panstarrs_v_minus_uprime\":0.0,\n",
    " \"panstarrs_v_minus_wprime\":-0.13,\n",
    " \"panstarrs_v_minus_wprime_source\":\"P\",\n",
    " \"panstarrs_v_minus_yprime\":0.32,\n",
    " \"panstarrs_v_minus_yprime_source\":\"P\",\n",
    " \"panstarrs_v_minus_zprime\":0.26,\n",
    " \"panstarrs_v_minus_zprime_source\":\"P\",\n",
    " \"perihelion_date\":\"2013-12-07.77843\",\n",
    " \"perihelion_date_jd\":\"2456634.27843\",\n",
    " \"perihelion_distance\":\"2.1304346\",\n",
    " \"period\":\"4.61\",\n",
    " \"pha\":false,\n",
    " \"phase_slope\":\"0.11\",\n",
    " \"q_vector_x\":\"-0.82462646\",\n",
    " \"q_vector_y\":\"-0.54106127\",\n",
    " \"q_vector_z\":\"0.16505729\",\n",
    " \"rc_minus_ic\":0.4,\n",
    " \"rc_minus_ic_source\":\"D\",\n",
    " \"residual_rms\":\"0.58\",\n",
    " \"saturn_moid\":null,\n",
    " \"semimajor_axis\":\"2.771723\",\n",
    " \"spin_amplitude_flag\":\"-\",\n",
    " \"spin_max_amplitude\":0.16,\n",
    " \"spin_min_amplitude\":0.03,\n",
    " \"spin_period\":7.81,\n",
    " \"spin_period_description\":null,\n",
    " \"spin_period_flag\":null,\n",
    " \"spin_updated_at\":\"2014-10-27T13:51:07Z\",\n",
    " \"taxonomy_class\":\"B\",\n",
    " \"taxonomy_class_source\":\"B\",\n",
    " \"tisserand_jupiter\":3.0,\n",
    " \"u_minus_b\":0.284,\n",
    " \"u_minus_b_source\":\"M\",\n",
    " \"updated_at\":\"2014-10-27T13:51:07Z\",\n",
    " \"uranus_moid\":null,\n",
    " \"v_minus_gprime\":-0.261,\n",
    " \"v_minus_gprime_source\":\"m\",\n",
    " \"v_minus_iprime\":0.321,\n",
    " \"v_minus_iprime_source\":\"d\",\n",
    " \"v_minus_rc\":0.361,\n",
    " \"v_minus_rc_source\":\"T\",\n",
    " \"v_minus_rprime\":0.157,\n",
    " \"v_minus_rprime_source\":\"m\",\n",
    " \"v_minus_uprime\":0.0,\n",
    " \"v_minus_wprime\":0.0,\n",
    " \"v_minus_yprime\":0.0,\n",
    " \"v_minus_zprime\":0.435,\n",
    " \"v_minus_zprime_source\":\"d\",\n",
    " \"venus_moid\":null}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(pv, label = 'Saturn')\n",
    "plt.plot(sv, label = 'Sun')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.xlim(14.5, 16.5)\n",
    "plt.show()\n",
    "plt.plot(pv, label = 'Saturn')\n",
    "plt.plot(sv, label = 'Sun')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.xlim(12, 20)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[31 28 31 30 31 30 31 31 30 31 30 31]\n"
     ]
    }
   ],
   "source": [
    "print dm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16:16:46.49\n"
     ]
    }
   ],
   "source": [
    "from compute import *\n",
    "from utils import parse_date, parse_coord\n",
    "import time\n",
    "\n",
    "#from vdb import process_vdb, handle_date\n",
    "\n",
    "#day = int(time.strftime(\"%d\"))\n",
    "month = 1#int(time.strftime(\"%m\"))\n",
    "year = int('20' + time.strftime('%y'))\n",
    "dm = dpm(year)\n",
    "vphi = 20\n",
    "phi= coord(vphi,0,0,1,'deg')\n",
    "\n",
    "for month in range(1,13):\n",
    "    days = range(1, dm[month-1]+1)\n",
    "    fnm = './csv/' + str(month)+str(year)+'.csv'\n",
    "    f = open(fnm, 'w+')\n",
    "    f.write('Date,Time,Dist,RA,DEC\\n')\n",
    "\n",
    "\n",
    "    for i in xrange(len(days)):\n",
    "        curdate = [days[i], month, year]\n",
    "        [[ra1, dec1], [ra2, dec2]] = ephsun(curdate)\n",
    "\n",
    "        if dec1.f2s() < dec2.f2s():\n",
    "            c = dec1\n",
    "            d = dec2\n",
    "        else:\n",
    "            c = dec2\n",
    "            d = dec1\n",
    "\n",
    "        ra_interest = [ra1.f2hd(), ra2.f2hd()]\n",
    "        dec_interest = [(c - phi).f2hd(), (d + phi).f2hd()]\n",
    "\n",
    "        #fields = ['RA', 'DEC', 'flux(V)']\n",
    "        customSimbad = Simbad()#columns = fields)\n",
    "\n",
    "        query = \"cat in ('\" + cat + \"')\" + ' & '\n",
    "        query += 'rah > '+ str(ra_interest[0])  + ' & ' + 'rah < ' + str(ra_interest[1])+ ' & '\n",
    "        query += 'dec > '+ str(dec_interest[0]) + ' & ' + 'dec < ' + str(dec_interest[1])\n",
    "\n",
    "        if cat == 'hip':\n",
    "            customSimbad.add_votable_fields('flux(V)', 'id(hip)')\n",
    "        elif cat == 'ngc':\n",
    "            customSimbad.add_votable_fields('id(m)')\n",
    "        customSimbad.ROW_LIMIT = -1\n",
    "\n",
    "\n",
    "        result = customSimbad.query_criteria(query)\n",
    "\n",
    "\n",
    "        if result is not None:\n",
    "            ll = len(result)\n",
    "            result.sort('RA')\n",
    "        else:\n",
    "            continue\n",
    "\n",
    "        for i in xrange(ll):\n",
    "            row = result[i]\n",
    "            [ra, dec] = [row['RA'], row['DEC']]\n",
    "            if cat == 'hip':\n",
    "                vmag = row['FLUX_V']\n",
    "            ra = parse_coord( ra, 'ra', ' ')\n",
    "            dec = parse_coord( dec, 'deg', ' ')\n",
    "            try:\n",
    "                t = comp_time(ra1, ra2, ra)\n",
    "            except:\n",
    "                print ra1, ra2, ra, row['MAIN_ID']\n",
    "                continue\n",
    "                \n",
    "            tmp = str(curdate[0])\n",
    "            if len(tmp) <  2:\n",
    "                tmp = ' ' + tmp\n",
    "            date = '' + tmp + '/'\n",
    "            tmp = str(curdate[1])\n",
    "            if len(tmp) < 2:\n",
    "                tmp = ' ' + tmp\n",
    "            date = date + tmp + '/' + str(curdate[2])\n",
    "            dist = comp_dist(curdate[0], curdate[1], t, dec)\n",
    "\n",
    "            hdist = abs(dist.f2s() / 3600.)\n",
    "            if (cat == 'hip')and False:\n",
    "                if (str(vmag)!='--'):\n",
    "                    if (vmag > 3):\n",
    "                        if (vmag <= 4):\n",
    "                            if hdist > 15:\n",
    "                                continue\n",
    "                        elif (vmag <= 5):\n",
    "                            if hdist > 14:\n",
    "                                continue\n",
    "                        elif (vmag <= 6):\n",
    "                            if hdist > 13:\n",
    "                                continue\n",
    "                        elif (vmag <= 7):\n",
    "                            if hdist > 12:\n",
    "                                continue\n",
    "                        elif (vmag <= 8):\n",
    "                            if hdist > 8:\n",
    "                                continue\n",
    "                        elif (vmag <= 9):\n",
    "                            if hdist > 4:\n",
    "                                continue\n",
    "                        else: #(vmag <= 10):\n",
    "                            if hdist > 2:\n",
    "                                continue\n",
    "                        #else:\n",
    "                            #print \"Strange! \" + row['MAIN_ID'], vmag, hdist\n",
    "                            #continue\n",
    "                            \n",
    "            printer = '' + date + ',' + t.__str__(rpar = 0) \n",
    "\n",
    "            stri = [row['MAIN_ID'], str(dist)]\n",
    "            if cat == 'hip':\n",
    "                alt_names = Simbad.query_objectids(row['MAIN_ID'])\n",
    "                stnms = ''\n",
    "                dl = ''\n",
    "                hd = ''\n",
    "                for an in alt_names:\n",
    "                    if an['ID'].startswith('HD'):\n",
    "                        if len(hd) > 0:\n",
    "                            hd = hd + ' / '\n",
    "                        hd = hd + an['ID']\n",
    "                    if an['ID'].startswith('* '):\n",
    "                        if len(dl) > 0:\n",
    "                            dl = dl + ' / '\n",
    "                        dl = dl + an['ID'].replace(\"*\", \"\")\n",
    "                    if an['ID'].startswith('NAME '):\n",
    "                        if len(stnms) > 0:\n",
    "                            stnms = stnms + ' / '\n",
    "                        tmp = an['ID'].replace(\"NAME\", '')\n",
    "                        stnms = stnms + tmp\n",
    "                stri = [row['ID_hip'], hd, dl, stnms, str(dist), str(vmag)]\n",
    "            elif cat == 'ngc':\n",
    "                mes = row['ID_m']\n",
    "                if (len(mes) == 0) and (hdist > 15):\n",
    "                    continue\n",
    "                alt_names = Simbad.query_objectids(row['MAIN_ID'])\n",
    "                stnms = ''\n",
    "                ngc = ''\n",
    "                for an in alt_names:\n",
    "                    if an['ID'].startswith('NGC '):\n",
    "                        if len(ngc) > 0:\n",
    "                            ngc = ngc + ' / '\n",
    "                        ngc = ngc + an['ID']\n",
    "                    if an['ID'].startswith('NAME '):\n",
    "                        if len(stnms) > 0:\n",
    "                            stnms = stnms + ' / '\n",
    "                        tmp = an['ID'].replace(\"NAME\", '')\n",
    "                        stnms = stnms + tmp    \n",
    "                stri = [ngc, mes, stnms, str(dist)]\n",
    "            stri += [str(ra), str(dec)]\n",
    "\n",
    "            for x in stri:\n",
    "                printer = printer + ',' + x\n",
    "            f.writelines(printer+'\\n')\n",
    "            #print printer\n",
    "        f.write('\\n')\n",
    "    f.close()\n",
    "    print \"Month\", month, \"is ready!\"\n",
    "print \"Check it!\""
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}