{
 "metadata": {
  "name": "Stellar-Color-Index"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import stats\n",
      "from scipy.optimize import curve_fit\n",
      "\n",
      "tau = 2.0 * math.pi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with open('data/hygfull.csv') as f:\n",
      "    print f.read(250) + '...'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import csv\n",
      "with open('data/hygfull.csv') as f:\n",
      "    stars = list(csv.DictReader(f))\n",
      "\n",
      "print 'The database holds', len(stars), 'stars'\n",
      "print 'An example star:', stars[123]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from spectral_classification import build_color_chart\n",
      "color_chart = build_color_chart('starcolors.txt')\n",
      "\n",
      "print 'The RGB color for class G2(V) is:\\n',\n",
      "print color_chart['G2(V)']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "usefuls = [\n",
      "    star for star in stars\n",
      "    if star['Spectrum'] in color_chart and\n",
      "       star['ColorIndex'] != ''\n",
      "    ]\n",
      "\n",
      "print 'There are', len(usefuls), '/', len(stars), 'useful stars'\n",
      "print 'Which is', 100.0 * float(len(usefuls)) / len(stars), '%'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "color_index = np.array(\n",
      "    [ float(star['ColorIndex']) for star in usefuls ])\n",
      "spectral_code = np.array(\n",
      "    [ star['Spectrum'] for star in usefuls ])\n",
      "\n",
      "print 'First few spectral codes:', spectral_code[:5]\n",
      "print 'First few color indexes:', color_index[:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rgb_color = np.array(\n",
      "    [ color_chart[code] for code in spectral_code ])\n",
      "\n",
      "print 'First few RGB colors:\\n', rgb_color[:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "red = rgb_color[:,0]\n",
      "green = rgb_color[:,1]\n",
      "blue = rgb_color[:,2]\n",
      "\n",
      "big = (8.0, 16.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(color_index, red, s=1.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(color_index, green, s=1.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(color_index, blue, s=1.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "result = stats.linregress(color_index, red)\n",
      "slope, intercept, r_value, p_value, std_err = result\n",
      "\n",
      "print slope, intercept, r_value, p_value, std_err"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(color_index, red, s=1.0)\n",
      "x0, x1 = -1, 1.5\n",
      "plot([x0, x1], [intercept + x0 * slope,\n",
      "                intercept + x1 * slope])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def fancy_regression(color):\n",
      "    result = stats.linregress(color_index[color < 1.0],\n",
      "                              color[color < 1.0])\n",
      "    slope, intercept, r_value, p_value, std_err = result\n",
      "\n",
      "    print slope, intercept, r_value, p_value, std_err\n",
      "\n",
      "    scatter(color_index, color, s=1.0)\n",
      "    x0, x1 = axis()[:2]\n",
      "    plot([x0, x1], [intercept + x0 * slope,\n",
      "                    intercept + x1 * slope])\n",
      "    axis([None, None, 0.3, 1.1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fancy_regression(red)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fancy_regression(blue)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def func(x, a, b, c):\n",
      "    return a * exp(-b*x) + c\n",
      "\n",
      "x = np.linspace(0,4,50)\n",
      "y = func(x, 2.5, 1.3, 0.5)\n",
      "yn = y + 0.2*np.random.normal(size=len(x))\n",
      "popt, pcov = curve_fit(func, x, yn)\n",
      "print popt\n",
      "print pcov"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def parabola(x, a, b, c):\n",
      "    return a*x*x + b*x + c\n",
      "\n",
      "popt, pcov = curve_fit(parabola, color_index, green)\n",
      "print popt\n",
      "print pcov"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "scatter(color_index, green, s=1.0)\n",
      "x0, x1, y0, y1 = axis()\n",
      "xspace = np.linspace(x0, x1, 100)\n",
      "yvalue = parabola(xspace, popt[0], popt[1], popt[2])\n",
      "plot(xspace, yvalue)\n",
      "axis([None, None, 0.6, 1.1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def indexed_color(color_index):\n",
      "    x = color_index\n",
      "    r =  0.448834253080 * x + 0.697688812105\n",
      "    g = -0.22906088 * x * x + 0.38220694 * x + 0.77957636\n",
      "    b = -0.349557269948 * x + 1.170004977010\n",
      "    return np.array((r, g, b)).clip(0.0, 1.0)\n",
      "\n",
      "print indexed_color(0.0)\n",
      "print indexed_color(1.0)\n",
      "print indexed_color(2.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xspace = np.linspace(-1.0, 3.0, 21)[2:-2]\n",
      "colors = indexed_color(xspace).transpose()\n",
      "scatter(xspace, xspace * 0.0 + 0.5,\n",
      "        s=100.0, c=colors)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}