Step 5 - First steps towards fitting models to data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 5 - First steps towards fitting models to data\n",
    "\n",
    "In the first three steps we've built up a framework that allows us to load stellar models, filter curves and combine them at any observed redshift, briefly we'll re-input that code, each of the three functions below are condensed versions of the ground we covered in the first three examples respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "\n",
    "cosmo = FlatLambdaCDM(H0=70., Om0=0.3)\n",
    "\n",
    "\n",
    "def get_model_grid():\n",
    "    \"\"\" Loads up the BPASS grid of stellar models and \n",
    "    resamples it onto a coarser wavelength grid. See Step 1. \"\"\"\n",
    "    \n",
    "    model_path = \"data/spectra-bin-imf135_300.z020.dat\"\n",
    "    raw_wavelengths = np.loadtxt(model_path, usecols=0)\n",
    "    raw_grid = np.loadtxt(model_path)[:,1:]\n",
    "    \n",
    "    grid = np.zeros((wavelengths.shape[0], raw_grid.shape[1]))\n",
    "\n",
    "    for i in range(grid.shape[1]):\n",
    "        grid[:,i] = np.interp(wavelengths, raw_wavelengths, raw_grid[:,i])\n",
    "    \n",
    "    grid *= (3.827*10**33)/(10**6)\n",
    "\n",
    "    return grid\n",
    "\n",
    "    \n",
    "def blueshift_filters(redshift):\n",
    "    \"\"\" A function that resamples filters onto the same wavelength\n",
    "    basis as the model spectrum at the specified redshift. See Step 2. \"\"\"\n",
    "        \n",
    "    resampled_filter_curves = []\n",
    "\n",
    "    for filt in filter_curves:\n",
    "        blueshifted_wavs = filt[:, 0]/(1 + redshift)\n",
    "        \n",
    "        resampled_filt = np.interp(wavelengths,\n",
    "                                   blueshifted_wavs, filt[:, 1],\n",
    "                                   left=0, right=0)\n",
    "        \n",
    "        resampled_filter_curves.append(resampled_filt)\n",
    "        \n",
    "    return resampled_filter_curves\n",
    "\n",
    "\n",
    "def get_model_photometry(redshift, age_index):\n",
    "    \"\"\" For a row in the model grid return model\n",
    "    photometry at the specified redshifts. See Step 3. \"\"\"\n",
    "\n",
    "    ssp_model = np.copy(grid[:, age_index])\n",
    "    \n",
    "    luminosity_distance = (3.086*10**24)*cosmo.luminosity_distance(redshift).value\n",
    "\n",
    "    ssp_model /= 4*np.pi*(1 + redshift)*luminosity_distance**2\n",
    "\n",
    "    filter_curves_z = blueshift_filters(redshift)\n",
    "\n",
    "    photometry = np.zeros(len(filter_curves))\n",
    "\n",
    "    redshifted_wavs = wavelengths*(1 + redshift)\n",
    "    \n",
    "    for i in range(photometry.shape[0]):\n",
    "        flux_contributions = filter_curves_z[i]*ssp_model\n",
    "        photometry[i] = np.trapz(flux_contributions, x=redshifted_wavs)\n",
    "        photometry[i] /= np.trapz(filter_curves_z[i], x=redshifted_wavs)\n",
    "        \n",
    "    return photometry\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also re-input the function we wrote in Step 4 to load up observational data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data(row_no):\n",
    "    \"\"\" Load UltraVISTA photometry from catalogue. See Step 4. \"\"\"\n",
    "\n",
    "    # load up the relevant columns from the catalogue.\n",
    "    catalogue = np.loadtxt(\"data/UltraVISTA_catalogue.cat\",\n",
    "                         usecols=(0,3,4,5,6,7,8,9,10,11,12,13,14,15,\n",
    "                                  16,17,18,19,20,21,22,23,24,25,26))\n",
    "    \n",
    "    # Extract the object we're interested in from the catalogue.\n",
    "    fluxes = catalogue[row_no, 1:13]\n",
    "    fluxerrs = catalogue[row_no, 13:25]\n",
    "\n",
    "    # Convert to microjanskys\n",
    "    fluxes = fluxes*10**29\n",
    "    fluxerrs = fluxerrs*10**29\n",
    "    \n",
    "    # Put an error floor into the data at the 20 sigma level.\n",
    "    # This is normally done to account for systematic uncertainties.\n",
    "    for i in range(fluxes.shape[0]):\n",
    "        if fluxerrs[i] < fluxes[i]/20:\n",
    "            fluxerrs[i] = fluxes[i]/20\n",
    "\n",
    "    return fluxes, fluxerrs\n",
    "\n",
    "fluxes, fluxerrs = load_data(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's also some other basic setting up we should do to define the variables we're going to need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define our basic quantities.\n",
    "wavelengths = np.arange(1000., 60000., 10.)\n",
    "ages = np.arange(2, 53)\n",
    "ages = 10**(6+0.1*(ages-2))\n",
    "\n",
    "grid = get_model_grid()\n",
    "\n",
    "# Load the curves up.\n",
    "filter_names = [\"data/filters/CFHT_u.txt\",\n",
    "                \"data/filters/CFHT_g.txt\",\n",
    "                \"data/filters/CFHT_r.txt\",\n",
    "                \"data/filters/CFHT_i+i2.txt\",\n",
    "                \"data/filters/CFHT_z.txt\",\n",
    "                \"data/filters/subaru_z\",\n",
    "                \"data/filters/VISTA_Y.txt\",\n",
    "                \"data/filters/VISTA_J.txt\",\n",
    "                \"data/filters/VISTA_H.txt\",\n",
    "                \"data/filters/VISTA_Ks.txt\",\n",
    "                \"data/filters/IRAC1\",\n",
    "                \"data/filters/IRAC2\"]\n",
    "\n",
    "filter_curves = []\n",
    "\n",
    "for name in filter_names:\n",
    "    filter_curves.append(np.loadtxt(name))\n",
    "    \n",
    "eff_wavs = np.zeros(len(filter_curves))\n",
    "\n",
    "# Calculate the effective wavelengths of the filter curves\n",
    "for i in range(len(filter_curves)):\n",
    "    filt = filter_curves[i]\n",
    "\n",
    "    wav_weights = filt[:,1]*filt[:,0]\n",
    "    \n",
    "    eff_wavs[i] = np.trapz(wav_weights, x=filt[:,0])\n",
    "\n",
    "    eff_wavs[i] /= np.trapz(filt[:, 1], x=filt[:,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now that that's over with...\n",
    "\n",
    "Having got all that out of the way let's think about how to go about fitting to data. The simplest way to do this is grid-based fitting. Basically we're going to pick a bunch of redshift, age and mass values making up a 3D grid, then calculate chi-squared values at each point and see which one is smallest. Let's start with a single point to develop the idea.\n",
    "\n",
    "We can only use discrete fixed age values corresponding to the ages we have spectra for in the BPASS model grid, so let's just leave things in terms of the column index in the grid instead of actual ages in years for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mass = 10**10  # Solar masses\n",
    "redshift = 1.\n",
    "age_index = 25\n",
    "\n",
    "# Our photometry is in erg/s/cm^2/A per Solar mass of stars\n",
    "model_photometry = get_model_photometry(redshift, age_index)\n",
    "\n",
    "# Convert to erg/s/cm^2/A by multiplying by the total mass we want.\n",
    "model_photometry *= mass\n",
    "\n",
    "# Our observational data is in microJanskys.\n",
    "fluxes, fluxerrs = load_data(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Would you look at that, we have a not-at-all-contrived unit mismatch between our model and our photometric data. Converting between flux per unit Angstrom and flux per unit frequency is the bane of my (and I presume other people in the field's too) existance, I always have to work it out from scratch every time.\n",
    "\n",
    "Let's turn our microJansky observational data into erg/s/cm^2/A and quickly plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert from muJy to erg/s/cm^2/Hz.\n",
    "fluxes = fluxes*10**-29\n",
    "fluxerrs = fluxerrs*10**-29\n",
    "\n",
    "# Convert from erg/s/cm^2/Hz to erg/s/cm^2/A.\n",
    "fluxes = fluxes*2.9979*10**18/eff_wavs**2\n",
    "fluxerrs = fluxerrs*2.9979*10**18/eff_wavs**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAE9CAYAAAB+ykFQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHJJJREFUeJzt3UFoXFeaL/D/EdmIEVidnYTshRpM4qXTT1Kvx8owgqFpUPptHmbaEHtp6ECat0pl1XjAgc5OCoR52fVEDDQDBW251y2JJLvOdALthYmlnVsGDd4EnbdQyVYrkixbUtUt1e8HhVW3rup+pSi69a9z7ndKrTUAAAD03lCvCwAAAGCHgAYAANAQAhoAAEBDCGgAAAANIaABAAA0xEAGtFLKaK9rAAAA2G/gAlop5WqSL/dtu1ZKuVpKudmjsgAAAAYvoNVav0qyuXu/E9gmO9sflFIme1YcAAAw0Bof0EopywdsmyylvN8Z+Xr/JFMWO8FstnOcyVrrg5PUCwAA8Kpe63UBhymlXEsymeTaAQ9/Xmt9q7PfF0k+SfJO5/5B0xQf11qXDjnOZJLfJXmQ5E4p5YtOaAMAAOiqUmvt7gFLuVlrXdy3bTTJtYNCVCml1lrLnvtXk9yptc7u2fa3WuuPXqKGL/cEvPeTLNZaNzt13Ky1/tvLvzIAAICT6cUUxy9KKXd273RC0SdJ7h/z+yez5xqyjsed4PZCuyNznX+TZCnJLzr3f5Jk8dBvBgAAOENdn+JYa/2qlJJOSPtNkjtJ3q217g9dh3n9hMe/n+RHe+4/iFAGAAA0QE+uQeuEtNez0+7+rZcIZ0nyOMn+piAnCm0AAABN0JMujp1pje8k+XWS//uS3/4gBwQyjT0AAIB+16s2+3eS/LrTFOR3e69Je5H9QazThfG4168BAAA0Vi+6ON5J8pu90xo7DT5+sre7Y2fbteyEuX9Lsty5fmzvY18luZpOF8buvQoAAIDT1/WABgAAwMF6NcURAACAfboe0EophuwAAICBdVQm6kmbfSENAADgh3oyxbHW6nYKt3fffbfnNXit/VdzL+vq9rHP+nhn8fz+Rrq9yq2pf2/68TZIP8t+fK1NrLnXNXXz+P14Xq21mefWo7gGrY/9y7/8S69L6Jp+fK1NrbmXdXX72Gd9vKb+N2bw+F08PYP0s+zH19rEmntdUzeP77zaHb1os1+7fUyAflFKeeEnawDA8TXx3NqpqRz0mBE0gAb54IMPel0CAJwr/XZuNYIGAADQRUbQAAAA+oCABgAA0BACGgAAQEMIaAAAAA0hoAEAADSEgAYAANAQAhpAA7VarZRSfnBrtVq9Lg0AOEPWQQNouM5aKb0uAwA4JdZBAzglRrYAgLNkBA3gFXVrZMsIGgCcL0bQAAAA+oCABgAA0BACGgAAQEMIaAAAAA0hoAE0VK01KysrSZKVlRWNQgBgAAhoAA3Ubrdz6dKVzM7eSJLMzv4yly5dSbvd7nFlAMBZEtAAXtJZj2y12+3Mz9/Id9/9Nltbf06SbG19ne+++23m528IaQBwjlkHDeAltNvt3Lr1XjY3S7a2/jsjI29kdDRZWLibubm5Ez9/rTUXL76ZR48+TvJ2Z2tJsvt3814mJm7n4cOvU8qBy6cAAA1nHTSAU9CNka3V1dU8eTKUZPaQPWazuZmsra2d+FgAQPMIaADHUGvNzZu/ytOnn2VnZGv3Q6+S5O08ffpZbt1678TTHTc2NjI0dHnP8+9XMjR0Oevr6yc6DgDQTAIawDF0a2RrbGws29vf5PmUxv1qtre/zfj4+ImOAwA0k4AGcAzdGtmanp7OhQs1yfIheyxndDSZmpo60XEAgGYS0ACOoVsjW6WULC5+lOHh60nu7TleTXIvw8PXs7BwV4MQADinBDSAY+jmyNbc3FyWlj7NxMTtjIxcSZKMjFzJxMTtLC19eirdIgGAZtJmH+CYdrs47jQKmc3OZ1zbSZYzPHz91MNTrTVra2uZmZnJyspKpqamjJwBwDlwVJt9AQ3gJTxfBy3Z2vrLqa+DdpDOH/EzeW4AoPsENIBT1O2RLQENAM4XC1UDnECr1Uop5dltaGgoMzMzSXauTTPtEAA4LUbQAF7C3tGsbo1sGUEDgPPFCBoAAEAfENAAAAAa4oUBrZRytZRys5Tyfinl81LK5BH73iml1FLK30opy0ftC8Dh9l73luTZ161Wq7eFAQBn6shr0Eopo0l+UWtd7Ny/lmSh1vrjQ/a/ubvvEc/pGjSgb/XiGjQA4Hw5yTVok0l+vef+F0kmO8HtRAX5NBjoV0a2AICXtX92zGFe2MWxlDJZa33Q+fpaks9rrT86ZN/3kzxIsplkNslvaq2b+/Yxggb0nVprVldX89Of/jR/+tOftNcHAF7Zibo47oazjltJ3j1i98Va61Kt9X6S3yX540tVCtBA7XY7ly5dyezsjSTJ7Owvc+nSlbTb7R5X1nv714gzoggAJ3PsddBKKTeTPK61Lh37yUupSX60dxTNCBrQT9rtdubnb+Tp08+yMzFgKMl2kuUMD1/P0tKnmZub622RDeGaPAA4nqNG0I4V0DpTG9MZGTtsn6tJPqm1vrVnW91/YAEN6Be11ly8+GYePfo4ydudrSXJ7t+we5mYuJ2HD7823TECGgAc14mmOHaC1+PdcFZKmd/z2N6GIQ+SLOx57FqSY4+2ATTN6upqnjwZys7I2UFms7mZrK2tdbMsAOAce+2oBzvrmP0xyeieT4cf5HnwupNkOTvXnm2WUh50pkImyY9z9PVqAI22sbGRoaHL2Rk1O0jJ0NDlrK+vd7MsAOAcOzKgdRqEHNixsfP4O/vuHzoFEqDfjI2NZXv7m+xMaTwopNVsb3+b8fHxLlcGAJxXL5ziCDCopqenc+FCzc5EgYMsZ3Q0mZqa6mZZAMA5JqABHKKUksXFjzI8fD3JvTxvDlKT3Mvw8PUsLNwd+AYhtdasrKwkSVZWVjQKAYATENAAjjA3N5elpU8zMXE7IyNXkiQjI1cyMXFbi/1YIw4ATtux10E7tQNqsw/0oVpr1tbWMjMzk5WVlUxNTQ38yJk14gDg1Zx4HbRTLkZAA/qWtb52WCMOAF7didZBA4D9rBEHAGdDQAPgpVkjDgDOhoAGwEv7+zXiDmKNOAB4FQIaAC/NGnEAcDYENIAXaLVaKaU8a3ax+3Wr1eptYT1kjTgAOBu6OALwytrtdm7dei+bm8nW1l8yMvJGRkeThYW7WuwDwCG02QfgzFgjDgBejoAGwJmzRhwAHI910AAAAPqAgAYAANAQAhoAAEBDCGgAAAANIaAB8MqsEQcAp0sXRwAAgC7SxREAAKAPCGgAAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAAAA0hoAEAADSEgAYAANAQAhoAAEBDCGgAAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADfHai3YopVxN8pMko0n+V5Jf11ofHLLvZJL5JF8luZpksda6eXrlAgAAnF9HBrRSymiSn9RaFzv3ryVZTvLjQ77l81rrW519v0jySZJ3Tq9cAACA8+tFUxwnk/x6z/0vkkx2gtvf6Yy0Pd693xk5u3YaRQIAAAyCIwNarfWrJLN7Nv0kyeYh0xYnk+zf/rgT3P5OKeXZrdVqvWTJAAAA/aXVaj3LQEcptdZjP2kp5fMkv6u1Lh3w2M0ks7XWd/Zs+2uSdzpBb3dbfZljAgAAnCellNRaD0xqx+7i2AlgB4azjsfZaSSy1+vHfX4AAIBB98Iujsmz5iAPaq33j9jtQQ4IZHtHzwAAADjcC0fQdpt/7IazUsr8nseeNQzZH8Q6LfePCnQAAADsceQ1aJ2Q9WX+furig1rrjzuPf55keU8b/qvZ6dx46DporkEDAAAG2VHXoL1Uk5BTKkZAAwAABtapNAkBAADgbAloAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAAAA0hoAHAAVqtVkopP7i1Wq1elwbAOWYdNKAvtFqtfPjhhz/Y/sEHH3jDzJnrrFfT6zIAOCcsVA2cK94s021+5wA4TRaqBoBXUGvNyspKkmRlZUVIA+DMCWhA3/BmmW5qt9u5dOlKZmdvJElmZ3+ZS5eupN1u97gyAM4zAQ3oC94s003tdjvz8zfy3Xe/zdbWn5MkW1tf57vvfpv5+Rt+7wA4M65BAxpv983y06efJZnNzmdL20mWMzx8PUtLn2Zubq63RXJu1Fpz8eKbefTo4yRvd7aWJLvnrnuZmLidhw+/TikHXj4AAEdyDRrQt2qtuXnzV51w9nZ23iin8+/befr0s9y69Z7pjpya1dXVPHkylJ0PAw4ym83NZG1trZtlATAgBDSg0bxZpts2NjYyNHQ5zz8M2K9kaOhy1tfXu1kWAANCQAMazZtlum1sbCzb29/k+ZTG/Wq2t7/N+Ph4N8sCYEAIaECjebNMt01PT+fChZpk+ZA9ljM6mkxNTXWzLAAGhIAGNJo3y3RbKSWLix9lePh6knt5/uFATXIvw8PXs7BwV4MQAM6EgAY0mjfL9MLc3FyWlj7NxMTtjIxcSZKMjFzJxMRtXUMBOFPa7AN9od1u59at97K5mWxt/SUjI29kdDRZWLjrzTJnptaatbW1zMzMZGVlJVNTUz4MAODEjmqzL6ABfcObZXqlcyLtdRkAnBMCGnCueLNMt/mdA+A0WagaAACgDwhoAHCAVquVUsqzabS7X7dard4WBsC5Zooj0BdarVY+/PDDH2z/4IMPvGEGAPqKa9AAAAAawjVoAAAAfUBAAwAAaAgBDQAAoCEENAAAgIYQ0AAAABpCQAMAAGgIAQ0AAKAhBDQAAICGENAAAAAaQkADAABoCAENAACgIY4V0Eopy8fY504ppZZS/lZKWS6lTJ68PAAAgMHx2lEPllKuJZlMcu0Yz/XXWms5laoAAAAG0JEjaLXW+7XWxW4VAwAAMMhO8xq00VLKfCnlWme64+gpPjcAAMC5d5oBbbHWulRrvZ/kd0n+eNiOpZRnt1ardYolAAAANE+r1XqWgY5Saq0vfLJSSn3Z68tKKTXJj2qtmwc818s8FQAAwLlRSslh+epURtBKKVdLKV/u374/nAEAAHC4Vw5opZTJPdeZPUiysOexa0mWTlgbAADAQHlRm/2r6bTYL6XcSbLcucYsSe4kWc7OtWebpZQHpZSbncd+nOTdM6oZAADgXDrWNWinekDXoAEAAAPszK9BAwAA4OQENAAYIHvbPFv2BqB5THEEgAHVmWLT6zIABo4pjgDAM7XWrKysJElWVlaENIAGEdAAYIC02+1cunQls7M3kiSzs7/MpUtX0m63e1wZAIkpjgAwMNrtdubnb+Tp08+SzGbnc9rtJMsZHr6epaVPMzc319siAQbAUVMcBTQAGAC11ly8+GYePfo4ydudrSXJ7jn5XiYmbufhw69TyoHvGQA4Ja5BA4ABt7q6midPhrIzcnaQ2WxuJmtra90sC4B9BDQAGAAbGxsZGrqcnVGzg5QMDV3O+vp6N8sCYB8BDQAGwNjYWLa3v8nzKY371Wxvf5vx8fFulgXAPgIaAAyA6enpXLhQkywfssdyRkeTqampbpYFwD4CGgAMgFJKFhc/yvDw9ST38nwkrSa5l+Hh61lYuKtBCECP6eIIAAOk3W7n1q33srmZbG39JSMjb2R0NFlYuKvFPkCXaLMPADxTa83a2lpmZmaysrKSqakpI2cAXaTNPgCQJGm1WhkaGsrMzEySZGZmJkNDQ2m1Wr0tDIAkRtAAAAC6yggaAABAHxDQAAAAGkJAAwAAaAgBDQAAoCEENAAAgIYQ0AAAABpCQAMAAGgIAQ0AAKAhBDQAAICGENAAAAAaQkADAABoCAENAACgIQQ0AACAhhDQAAAAGkJAAwAAaAgBDQAAoCEENAAAgIYQ0AAAABpCQAMAAGgIAQ0AAKAhBDQAAICGENAAAAAa4lgBrZSyfIx9Jksp75dSrnX+HT15eQAAAIOj1FoPf7CUa0kmkyzUWsuRT1TKl7XWtzpfjyb5pNb6zgH71aOOCQAAcJ6VUnJYvjoyoO15gnpUQCulXE1yp9Y6u2fb32qtPzrkuY5XOQAAwDlzVEA7rWvQJpNs7tv2uBPcAAAAOIbTCmivv8zOpZRnt1ardUolAAAANFOr1XqWgY5yWgHtcZL9TUEODW211mc3AQ16b+8fDB+gAACcvlar9SwDHeW1UzregxwQyGqtX53S8wNnaGpqKhMTb2Rzs2Rr678zMvJGRkd3tgMA0D2vPILWaas/mvwwiJVSJpPcP2FtQBe02+3Mz9/Id9/9Nltbf06SbG19ne+++23m52+k3W73uEIAgMHxojb7V5NcS3Inyb8lWa613u889nnn/uK+fb9KcjXJYq11f+MQXRyhQWqtuXjxzTx69HGStztbS5Ld/0fvZWLidh4+/PqF86UBADieE7fZP+ViBDRoiJWVlczO3uiMnO3+jdgb0GpGRq7k/v1/z/T0dG+KBAA4Z7rRZh/oQxsbGxkaupzn4Wy/kqGhy1lfX+9mWQAAA0tAgwE2NjaW7e1v8nzEbL+a7e1vMz4+3s2yAAAGloAGA2x6ejoXLtQky4fssaybIwBAFwloMMBKKVlc/CjDw9eT3Mvea8+Sexkevp6FhbsahAAAdImABgNubm4uP//5Pyb5pzz/kzCU5J/y85//Y+bm5npXHADAgNHFEUiy03J/bW0t6+vrGR8fz9TUlJEzAIAzoM0+AABAQ2izDwAA0AcENAAAgIYQ0AAAABpCQAMAAGgIAQ0AAKAhBDQAAICGENAAAAAaQkADAABoCAENAACgIQQ0AACAhhDQAAAAGkJAAwAAaAgBDQAAoCEENAAAgIYQ0AAAABpCQAMAAGgIAQ0AAKAhBDQAAICGENAAABqi1WqllPKDW6vV6nVpQJeUWmt3D1hK7fYxAQD6Ra01q6ur+elPf5o//elPmZ6eTiml12UBp6iUklrrgf9jG0EDAGiIdrudS5euZHb2RpJkdvaXuXTpStrtdo8rA7rFCBoAQAO02+3Mz9/I06efJZnNzufo20mWMzx8PUtLn2Zubq63RQKn4qgRNAENAKDHaq25ePHNPHr0cZK3O1tLkt33TPcyMXE7Dx9+bbojnAOmOAIANNjq6mqePBnKzsjZQWazuZmsra11syygBwQ0AIAe29jYyNDQ5eyMmh2kZGjoctbX17tZFtADAhoAQI+NjY1le/ubPJ/SuF/N9va3GR8f72ZZQA8IaAAAPTY9PZ0LF2qS5UP2WM7oaDI1NdXNsoAeENAAAHqslJLFxY8yPHw9yb08H0mrSe5lePh6FhbuahACA0AXRwCAhmi327l1671sbiZbW3/JyMgbGR1NFhbuarEP54g2+wAAfaLWmrW1tczMzGRlZSVTU1NGzuCc0WYfAKAPtFqtDA0NZWZmJkkyMzOToaGhtFqt3hYGdI0RNDiBWmtWV1ezsbGRsbGxTE9P+5QTAIAjHTWC9toxvnkyyXySr5JcTbJYa908ZN87Sd5PspnkiyS3aq0PXrVwaLLn1wnsrE2zvf2N6wQAADiRF46glVK+rLW+1fl6NMkntdZ3Dtn3Zq118QXPZwSNvh95arfb+dnPfpHvv/+fHzz22mv/kN///j+ENAAADvTK16CVUq4mebx7vzNydu00Ctq9mVM9eNrtdi5dupLZ2Rv513/9f5md/WUuXbqSdrvd69KOpdaamzd/le+//8/stD/e2wq55vvv/zO3br0XH0QAALCr1Wo9y0BHOXIErZQyn+R/7x0xK6X8Nck7tdavDtj//SQPsjPFcTbJb/ZPhzSCNtja7Xbm52/k6dPPsvMrUrITbJYzPHw9S0ufNn7kaWVlJbOzN7K19efs1J88fx1JUjMyciX37/97pqene1MkAACNdZJr0F5/yWM9uz6tlPI4yR+TvPWSz8E5tTvy9PTpbJJ/+sHjT5/+n9y69V4ePvznRk933NjYyNDQ5TwPZ/vtXJO2vr7ezbIAADgHXtRm/3GS0X3bDg1te0fLOiNsVzvXrUFWV1fz5MlQks9y0NTA5LNsbiZra2u9KvFYxsbGsr39TZ7Xv1/N9va3GR8f72ZZAACcAy8KaA9yQCA7ZHrj1VLKlwfse2DHRwbPeRl5mp6ezoULNcn17LyWvdMcS5LrGR1NpqamelUiAAB96siAtj+IdVru3997f88I2YMkC3seu5Zk6fRKpd+dl5GnUkoWFz/K8PBykj8k2c7Oa9pO8ocMDy9nYeFuo6dpAgDQTMdps381O50bf7AOWinl8yTLu631O6FssvOtP44mIexRa83Fi2/m0aOPk7zd2bq3uca9TEzczsOHX/dFuHm+Dlo666B9ax00AABe6KgmIS8MaGdQjIA2wM7b+mG11qytrWV9fT3j4+OZmprqi3AJAEDvCGg0ipEnAAAGmYBG4xh5AgBgUAloAAAADXFUQHtRm30AAAC6READAABoCAENAACgIQQ0AACAhhDQAAAAGkJAAwAAaAgBDQAAoCEENAAAgIYQ0AAAABritV4XAADA+VFrzerqajY2NjI2Npbp6emUUnpdFvQNAQ0AgFPRbrdz69Z72dwsGRq6nO3tbzI6miws3M3c3Fyvy4O+UGqt3T1gKbXbxwQA4Gy12+387Ge/yPff/88PHnvttX/I73//H0IadJRSUms9cGhZQAMA4ERqrbl48c08evRxkrc7W0uS3fd89zIxcTsPH35tuiPk6ICmSQgAACeyurqaJ0+GkswessdsNjeTtbW1bpYFfUlAAwDgRDY2NjI0dDk7o2YH2bkmbX19vZtlQV8S0AAAOJGxsbFsb3+T51Ma96vZ3v424+Pj3SwL+pKABgDAiUxPT+fChZrkenZG0XZH0na/vp7R0WRqaqpXJULfENAAADiRUkoWFz/K8PBykj8k2c7OaNp2kj9keHg5Cwt3NQiBY9DFEQCAU/F8HbR01kH71jpocABt9gEA6Ipaa9bW1rK+vp7x8fFMTU0ZOYN9BDQAAICGsA4aAABAHxDQAAAAGkJAAwAAaAgBDQAAoCEENAAAgIYQ0AAAABpCQANokFar1esSAOBc6bdzq4DWx/7rv/6r1yV0TT++1qbW3Mu6un3ssz7eWTz/hx9+eOrPyfnX1L83/WiQfpb9+FqbWHOva+rm8fvxvJr037lVQOtjvf6D0E39+FqbWrOA1j/PD8fld/H0DNLPsh9faxNr7nVNAtr5U2qt3T1gKd09IAAAQMPUWstB27se0AAAADiYKY4AAAANIaABAAA0hIAG0AdKKaOllIVe1wEA50Ep5ctSynIp5f1e17Lfa70uAIBjuZZkstdFAMA58eta6/1eF3EQI2gAXVZKWT5g22Qp5f1SyrXOv6N7HruWpJEnEQBogpc9tya52tl+p4tlHosujgBd0glak0kW9rfWLaV8WWt9q/P1aJJPaq3vlFImk6TW+qCUslxrne164QDQUK9ybt23z/tJNmuti92q+UWMoAF0Sa31/kEngFLK1SSP9+y3mZ0pjUlyNclkKWU+yeudfQGAvNq5tZQy3wl2ux6nQQQ0gN6bTLK5b9vjUsrVWuvSnjnyowEAjuPQc2s6lw3shrRa61KXazuSJiEAvff6i3bonDwadQIBgAY79NzaGU3b/fCzcdd4G0ED6L3H+eHo2AtDGwBwqL49twpoAL33IAecNGqtX/WgFgA4D/r23CqgAfTY/pNFp3Nj46ZcAEC/6Odzq2vQALqkc2HybgepO0mW9zQAebfT6ver7HRufLc3VQJA/ziP51broAEAADSEKY4AAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAAAA3x/wGYCmOfOjKxTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.errorbar(eff_wavs, fluxes, yerr=fluxerrs, lw=1.0, linestyle=\" \", capsize=3, capthick=1, color=\"black\")\n",
    "plt.scatter(eff_wavs, fluxes, s=75, linewidth=1, facecolor=\"blue\", edgecolor=\"black\")\n",
    "plt.xscale(\"log\")\n",
    "plt.xlim(1000.*(1 + redshift), 60000.*(1 + redshift))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, now we have our units the same let's calculate a chi-squared value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi-squared value: 3505.00879547\n"
     ]
    }
   ],
   "source": [
    "diffs = fluxes - model_photometry\n",
    "\n",
    "chisq = np.sum(diffs**2/fluxerrs**2)\n",
    "\n",
    "print \"Chi-squared value:\", chisq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not great, let's see if we can iterate towards a better solution. First plot this model over the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAE9CAYAAAB+ykFQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3V9opeedJ/jvI1eckaemS0llYVRW1UWl8TqGmYFKLFUgsJAupWlBr8msPA2dYDoGu+7GbBuc2auc7E3jgMPad5IhZE13Q49rcEIvB8ZV3bAQiKTYvlnijkN3XVSXSnuRclRBazWOrWcvzlFZUelflVTnvNL5fOAgnfd96pyfHEXnfM/zPL+31FoDAABA/w31uwAAAAA6BDQAAICGENAAAAAaQkADAABoCAENAACgIQYyoJVSRvpdAwAAwGYDF9BKKeeSvL3p2IVSyrlSyrN9KgsAAGDwAlqt9Z0ky+v3u4HtbPf41VLK2b4VBwAADLTGB7RSyuUtjp0tpbzQnfl6YT9LFrvBbLL7PGdrrVf3Uy8AAMC9OtbvArZTSrmQ5GySC1ucfr3W+sXuuLeSvJrkye79rZYpvl9rvbTN85xN8jdJriZ5sZTyVje0AQAA9FSptfb2CUt5ttY6u+nYSJILW4WoUkqttZYN988lebHWOrnh2K9rrZ+5ixre3hDwXkgyW2td7tbxbK31e3f/kwEAAOxPP5Y4vlVKeXH9TjcUvZrkyh7//dls2EPW9X43uO1qfWau+zVJLiX5T937X0oyu+0/BgAAuI96vsSx1vpOKSXdkPYXSV5M8kytdXPo2s5n9/n8V5J8ZsP9qxHKAACABujLHrRuSPtsOu3uv3gX4SxJ3k+yuSnIvkIbAABAE/Sli2N3WeOTSb6d5H+7y39+NVsEMo09AACAw65fbfZfTPLtblOQv9m4J203m4NYtwvjXvevAQAANFY/uji+mOQvNi5r7Db4+NLG7o7dYxfSCXPfS3K5u39s47l3kpxLtwtj734KAACAg9fzgAYAAMDW+rXEEQAAgE0ENAAAgIboeZv9Uoo1lQAAwECrtZatjvdlBq3W6nYAt2eeeabvNfhZD1/N/ayr1899v5/vfjy+v5Fu93Jr6t+bw3gbpP+Wh/FnbWLN/a6pl89/GF9Xa23ma+tOLHE8xP74j/+43yX0zGH8WZtacz/r6vVz3+/na+r/xgwev4sHZ5D+Wx7Gn7WJNfe7pl4+v9fV3uhHm/3a6+cEOCxKKbt+sgYA7F0TX1u7NTVniSMAW/vOd77T7xIA4Eg5bK+tZtAAAAB6yAwaAADAISCgAQAANISABgAA0BACGgAAQEMIaAAAAA0hoAEAADSEgAbQQK1WK6WUO26tVqvfpQEA95HroAE0XPdaKf0uAwA4IK6DBnBAzGwBAPeTGTSAe9SrmS0zaABwtJhBAwAAOAQENAAAgIYQ0AAAABpCQANoqFpr5ubmkiRzc3P2oQHAABDQABqo3W7nzJnHMjn5dJJkcvJbOXPmsbTb7T5XBgDcTwIawF263zNb7XY709NP5/r1l7Oy8vMkycrKu7l+/eVMTz8tpAHAEabNPsBdaLfbuXjx+Swvl6ys/EOOH380IyPJzMxLmZqa2vfj11pz+vQXsrj4SpKvdY+WJOt/N9/M2NhzuXbt3ZSyZXdeAKDhtNkHOAC9mNman5/PrVtDSSa3GTGZ5eVkYWFh388FADSPgAawB7XWPPvsn2d19bV0ZrbWP/QqSb6W1dXXcvHi8/te7ri0tJShoUc2PP5mJUNDj+TGjRv7eh4AoJkENIA96NXM1ujoaNbW3ssnSxo3q1lb+2VOnTq1r+cBAJpJQAPYg17NbE1MTOTEiZrk8jYjLmdkJBkfH9/X8wAAzSSgAexBr2a2SimZnf1+hoefSvLmhuerSd7M8PBTmZl5SYMQADiiBDSAPejlzNbU1FQuXfpBxsaey/HjjyVJjh9/LGNjz+XSpR8cSLdIAKCZtNkH2KP1Lo6dRiGT6XzGtZbkcoaHnzrw8FRrzcLCQs6fP5+5ubmMj4+bOQOAI2CnNvsCGsBd+OQ6aMnKyi8O/DpoW+n+Eb8vjw0A9J6ABnCAej2zJaABwNHiQtUA+9BqtVJKuX0bGhrK+fPnk3T2pll2CAAcFDNoAHdh42xWr2a2zKABwNFiBg0AAOAQENAAAAAaYteAVko5V0p5tpTyQinl9VLK2R3GvlhKqaWUX5dSLu80FgAAgN+1Y0ArpYwk+VKtdbbW+r0kM9n+Kq1J8k+11lJr/UytdbLWevUgiwUYFBsbkyS5/X2r1epvYQDAfbVjk5BSyrkkr9daP9+9P5Lk10k+U2td3mL8s7XW2R2fUJMQ4BDrR5MQAOBouecmIbXWd5JMbjj0pSTLW4WzrpFSynQp5UJ3uePIdgX5NBg4TGqtmZubSxIzWwDAXdu8OmY7d9Vmv5TyepK/qbVe2ub8yHp4686+vVpr/eKmMWbQgEOl3W7n4sXns7xcsrLyDzl+/NGMjCQzMy9lamqq3+UBAIfMgbTZL6U8mx3CWZJsnFnrzr6d224WDeAwaLfbmZ5+Otevv5yVlZ8nSVZW3s316y9nevrptNvtPlfYX5sv4m1GEQD2Z08zaKWUC0lSa72yw5g7Zsy6s2Vl0zgzaMChUGvN6dNfyOLiK0kmM/H785n/xy9n4vd/mvl/nEhyOWNjz+XatXd3Xa4wCOzJA4C92dcMWjd4vb8ezkop0xvOnd0wQ3Y1nS6P6+cuJNl2tg2g6ebn53Pr1lD+6D/8NtdeOZPL/6WzJffyf5nMtVfO5I/+w2+zvJwsLCz0uVIA4Kg4ttPJ7nXM/i6d5h/rh6/mk+D1Yjpt92drrcullKvdpZBJ8vkkzxx8yQC9sbS0lD/8dyfy2jNP5qFPr94+/m+GV/Jvhldy6bkn89Sr/z43btzoY5UAwFGyY0DrXsfsMzucf3LT/W2XQAIcNqP/9t/m5W+8lYc+/dGW5x/69Gr+jz99O4ujoz2uDAA4qvbcJARg0Ez8fnJi+OMdx4z8648z/vkeFQQAHHkCGsA2yr/8v/n0vxreccynPz2c8i9LPaqomTZeI25ubk6jEADYBwENYDvDo/nUsZ3/TH7q2FAyfKpHBTVPu93OmTOPZXLy6STJ5OS3cubMYwN/+QEAuFcCGsB2Tk4knzqx85gHR5KT472pp2FcIw4ADp6ABrCdUpLx2eSBbZY5Dj2Y/I//a29raohaa5599s+zuvpakq8lWe/0W5J8Laurr+XixectdwSAuySgAezk4ankK5eSh8aSY8e7B0vnVo4l/893kh+fSRYHa7Zo/RpxyeQ2IyZdIw4A7oGABrCbh6eSJ64l/+5/7x6ondvHHyQfrSQfXE9+Mj1QIW1paSlDQ4/kk5mzzUqGhh5xjTgAuEsCGsBe/eKl7c99vJr87GIyIEv6RkdHs7b2XjphdSs1a2u/zKlTg9tABQDuhYAGsBc355Pf3tp5zIfLyc3BWNI3MTGREydqksvbjLickZFkfHwwG6gAwL0S0AB20Wq1Uv6HL6f8yUqSpHyjc2v9t80jh5LVwVjSV0rJ7Oz3Mzz8VJI388lMWk3yZoaHn8rMzEspZbslkADAVkqvO2yVUqquXsCh86u55O8nO3vOtnPsePLVK8nnJnpXV5+12+1cvPh8lpeTlZVf5PjxRzMykszMvJSpqal+lwcAjVRKSa11y08xBTSAvag1+dHpZHVx+zEPjXWaiQzYrFGtNQsLCzl//nzm5uYyPj5u5gwAdiCgARyExXanW+PHq3eee2C4047/4cGdNeq+2PS7DABoPAEN4KAstjvdGj9cTmcb71ry4Ejy+MxAh7NEQAOAvRLQAA5SrZ1ujas3kuFTycnxgVvWuBUBDQD2RkAD4L4T0ABgb3YKaNrsAwAANISABsA9a7VaKaXc7tq4/n2r1epvYQBwSFniCAAA0EOWOAIAABwCAhoAAEBDCGgAAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAAAA0hoAEAADSEgAYAANAQAhoAAEBDCGgAAAANIaABAAA0hIAGAADQEMd2G1BKOZfkS0lGkjye5Nu11qvbjD2bZDrJO0nOJZmttS4fXLkAAABH144BrZQykuRLtdbZ7v0LSS4n+fw2/+T1WusXu2PfSvJqkicPrlwAAICja7cljmeTfHvD/beSnO0Gt9/RnWl7f/1+d+bswkEUCQAAMAh2DGi11neSTG449KUky9ssWzybZPPx97vB7XeUUm7fWq3WXZYMAABwuLRardsZaCel1rrnBy2lvJ7kb2qtl7Y492ySyVrrkxuO/VOSJ7tBb/1YvZvnBAAAOEpKKam1bpnU9tzFsRvAtgxnXe+n00hko8/u9fEBAAAG3a5dHJPbzUGu1lqv7DDsarYIZBtnzwAAANjerjNo680/1sNZKWV6w7nbDUM2B7Fuy/2dAh0AAAAb7LgHrRuy3s7vLl28Wmv9fPf860kub2jDfy6dzo3bXgfNHjQAAGCQ7bQH7a6ahBxQMQIacG9qTW7OJ6tLyfBocnIi2aUTEgBA0+wU0Pa0Bw2g7xbbyc8uJh8up7M6ey15cCR5fCZ5eKrf1QEAHAgzaEDzLbaTn0wnH6/eee6B4eQrl4Q0AODQsMQROLxqTX50Olld3H7MQ2PJE9csdwQADoUDuQ4aQF/cnE9+e2vnMR8uJzcXelNPL9Wa/Gou+ec3Ol99uAUAR549aECzrS5l98+ShpLVG72opnfsuQOAgWQGDWi24dEka7sMWkuGT/Wimt5Y33P3wfXko5Xko990vn5wvXN8sd3vCgGA+0RAA5rt5ETyqRM7j3lwJDk53pt67rdak4Vnt26IknSO/+yi5Y490Gq1Ukq549ZqtfpdGgBHmCYhQPMtttN67ol89799dMep7/wvx9J6+cdHZ9nfr+aSv5/szJht59jx5KtXks9N9K6uAdfdzN3vMgA4IlwHDTjcHp5K6+Ufp/XNzp6s8icrqX9z/GjuyRrUPXcAQBIBDTgsHp5KTl3rdms835lBOjl+9FrrD+KeuwartWZ+fj5JMjc3l4mJiZSj9jsHQKPYgwYcGjXJ3D92lpnN/WPNkVxwNmh77hqs3W7nzJnHMjn5dJJkcvJbOXPmsbTbmrQAcP8IaMChMDBvlktJxmeTB4a3Pv/AcGdZp1mc+6rdbmd6+ulcv/5yVlZ+niRZWXk316+/nOnpp4/e7x0AjaFJCNB462+WV1dfSzKZ29cFy+UMDz+VS5d+kKmpI7QPLXEdtD6qteb06S9kcfGVJF/rHi3J7TnbNzM29lyuXXvXckcA7slOTUIENKDRBvrNcq2dPXerNzp7zo7inrsGmpuby+Tk092Zs/X/3ht/52qOH38sV678MBMTOmkCcPd0cQQOrfn5+dy6NZTOzNlWJrO8nCwsLBy9N8ulaKXfB0tLSxkaeiSfhLPNSoaGHsmNGzppAnDwBDT6o9bk5nynpfjwaKcxgpkBtuDNMr02OjqatbX30pkx2+r3rmZt7Zc5dUonTQAOnoBG79lbw13wZplem5iYyIkTNSsrl/PJstqNLmdkJBkf10kTgIOniyO9tdhOfjKdfHA9+Wgl+eg3na8fXO8cX9QZjd+1/mY5ubzNCG+WOVillMzOfj/Dw08leTMb954lb2Z4+KnMzLx09PY8AtAIAhq9U2uy8Gzy8erW5z9e7cysaSLDBt4s0w9TU1O5dOkHGRt7LsePP5YkOX78sYyNPXc0u4YC0Bi6ONI7v5pL/n6yM2O2nWPHk69e0RiBO7Tb7Vy8+HyWl5OVlV/k+PFHMzKSzMy85M0y902tNQsLCzl//nzm5uYyPj7uwwAA9k2bfZrhn99IfvpnnWWN2zn2e8mXf5ic/nqvquIQ8WaZfum+kPa7DACOCG32aYbh0XQuLryTtc71nmALpZTbrfSPXEt9AIDYg0YvnZxIPnVi5zEPjnQuxgvQZ61WK6WU27O069+3Wq3+FgbAkSag0TulJOOzyQPDW59/YLjTat+SNbbgzTK91mq1Umu94+Z3DoD7yR40es910AAAGGCahNA8tSY3F5LVG509ZyfHzZwBADAQBDQAAICG2Cmg2YMGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAAAA0hoAEAADSEgAYAANAQewpopZTLexjzYimlllJ+XUq5XEo5u//yAAAABsexnU6WUi4kOZvkwh4e659qreVAqgIAABhAO86g1Vqv1Fpne1UMAADAIDvIPWgjpZTpUsqF7nLHke0GllJu31qt1gGWAAAA0DytVut2BtpJqbXu+mCllLrb8sVSykitdbn7/bkkr9Zav7jNY+36nAAAAEdRKSXb5asDC2hb/Zskn1kPbZse624eCgAA4MjYKaAdyBLHUsq5Usrbm49vDmcAAABs754DWinl7IZ9ZleTzGw4dyHJpX3WBgAAMFB2a7N/Lt0W+6WUF5NcrrVe6Z5+McnlJLO11uVSytVSyrPdc59P8sx9qhkAAOBI2tMetAN9QnvQAACAAXbf96ABAACwfwIaAAyQjdfhcV1SgOaxxBEABkytNfPz8/nyl7+cn/70p5mYmNj1wqkAHBxLHAGAJEm73c6ZM49lcvLpJMnk5Ldy5sxjabfbfa4MgMQMGgAMjHa7nenpp7O6+lqSyXQ+p11LcjnDw0/l0qUfZGpqqr9FAgyAnWbQBDQAGAC11pw+/YUsLr6S5GvdoyXJ+mvymxkbey7Xrr1ruSPAfWaJIwAMuPn5+dy6NZTOzNlWJrO8nCwsLPSyLAA2EdAAYAAsLS1laOiRdGbNtlIyNPRIbty40cuyANhEQAOAATA6Opq1tffyyZLGzWrW1n6ZU6dO9bIsADYR0ADYn1qTX80l//xG56t9xo00MTGREydqksvbjLickZFkfHy8l2UBsMmxfhcAwCG22E5+djH5cDm3OwI+OJI8PpM8rBtgk5RSMjv7/UxPP7Whi2PSmVHrdHGcmfmBBiEAfaaLIwD3ZrGd/GQ6+Xj1znMPDCdfuSSkNVC73c7Fi89neTlZWflFjh9/NCMjyczMS1rsA/SINvsAHKxakx+dTlYXtx/z0FjyxLXEjEzj1FqzsLCQ8+fPZ25uLuPj42bOAHpIm30ADtbN+eS3t3Ye8+FyclPL9qZptVoZGhrK+fPnkyTnz5/P0NBQWq1WfwsDIIkZNNifWjtvVFeXkuHR5OSE2QIGwz+/kfz0z5KPfrP9mGO/l3z5h8npr/eqKgA4FHaaQdMkBO7VUWuOIGxyN4ZHk6ztMmgtGdayHQDuhhk0uBdHrTnCUQub3H/2oAHAPbMHDQ5SrcnCs1uHs6Rz/GcXD8+1oNbD5gfXk49WOkvWPlrp3P/JdOc8bFZKMj7b+UBiKw8MdwK+cAYAd0VAg7t1lJojHLWwSW89PNWZLX5oLDl2vLPn7Njxzv3DNosMAA1hDxrcrdWl7P7ZxlCyeqMX1ezP3YTNz030piYOl4enklPXOr8jqzc6e85Ojps5A4B7JKDB3TpKzRGOUtikf0oR4AHggFjiCHfr5ETyqRM7j3lwpDOL0HRHKWwCABwBAhrcraPUHOEohU0AgCNAQIN7cVSaIxylsAkAcAS4DhrsR61HozmC66ABAPTMTtdBE9CAjqMSNgEAGk5AAwAAaIidApo9aAAAAA0hoAEAADSEgAYAANAQAhoAAEBDCGgAAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADbGngFZKubyHMWdLKS+UUi50v47svzwAAIDBUWqt258s5UKSs0lmaq1lxwcq5e1a6xe7348kebXW+uQW4+pOzwkAAHCUlVKyXb7aMaBteIC6U0ArpZxL8mKtdXLDsV/XWj+zzWPtrXIAAIAjZqeAdlB70M4mWd507P1ucAMAAGAPDiqgffZuBpdSbt9ardYBlQAAANBMrVbrdgbayUEtcZxOcnHzEsckf1BrfWeLx9r1OYHeqrVmfn4+S0tLGR0dzcTExK5/QAAAuHs7LXE8dkDPcTVbzKJtDmdAM33jG9/IX//1X99x/E//9E/zV3/1V32oCABgMN3zEsduW/2R5M4gVko5m+TKPmsDeqDdbueNN/4uyX9PstY9upbkv+eNN/4u7Xa7f8UBAAyY3drsn0tyIcmLSb6X5HKt9Ur33Ovd+7Obxr6T5FyS2Vrr5sYhljhCg9Rac/r0F7K4+EqSr3WPliTr/x99M2Njz+XatXctdwQAOCD7brN/wMUIaNAQc3NzmZx8OisrP08nmCW/G9Bqjh9/LFeu/DATExP9KRL6qdbk5nyyupQMjyYnJxIfVgCwT73YgwYcQktLSxkaeiSfhLPNSoaGHsmNGzd6WRY0w2I7+dnF5MPldHYErCUPjiSPzyQPT/W7OgCOqINqsw8cQqOjo1lbey+fzJhtVrO29sucOnWql2VB/y22k59MJx9cTz5aST76TefrB9c7xxftzQTg/hDQYIBNTEzkxIma5PI2Iy5nZCQZHx/vZVnQX7UmC88mH69uff7j1c7MmuX6ANwHAhoMsFJKZme/n+Hhp5K8mY17z5I3Mzz8VGZmXtIghMFycz757a2dx3y4nNxc6E09AAwUAQ0G3NTUVL7+9T9I8of55E/CUJI/zNe//geZmrLXhgGzupTdXx6HklV7MwE4eLo4Akk6LfcXFhZy48aNnDp1KuPj42bOGEy/mkv+frKz52w7x44nX72SfE53UwDunjb7ALBXtSY/Op2sLm4/5qGx5IlrWu4DcE92CmiWOALARqUk47PJA8Nbn39guNNqXzgD4D4wgwYAW3EdNADuE0scAeBe1Nrp1rh6Ixk+lZwcN3MGwL4JaAAAAA1hDxoAAMAhIKABAAA0hIAGAADQEAIaAABAQxzrdwEAQI/VmtycT1aXkuHR5OSE7pQADSGgAcAgcX03gEbTZh8ABsViO/nJdPLx6p3nHhhOvnJJSAPoAddBA4BBV2vyo9PJ6uL2Yx4aS564ZrkjwH3mOmgAMOhuzie/vbXzmA+Xk5sLvakHgC0JaAAwCFaXsvvL/lCyeqMX1QCwDQENAAbB8GiStV0GrSXDp3pRDQDbENAAYBCcnEg+dWLnMQ+OJCfHe1MPAFsS0ABgEJSSjM92ujVu5YHhTqt9DUIA+koXRwAYJK6DBtB3ujgCAB0PT3Va6X/1SvLlH3a+PnFNOGuIVquVUsodt1ar1e/SgB451u8CAIAeKyX53ES/q2AL4+PjGRt7NMvLJSsr/5Djxx/NyEjnODAYLHEEAGiAdrud6emns7r6WpLJ3F6CmssZHn4qly79IFNTZjrhKNhpiaOABgDQZ7XWnD79hSwuvpLka92jJcn6e6Y3Mzb2XK5dezdFIxc49OxBAwBosPn5+dy6NZTOzNlWJrO8nCwsLPSyLKAPBDQAgD5bWlrK0NAj6cyabaVkaOiR3Lhxo5dlAX0goAEA9Nno6GjW1t7LJ0saN6tZW/tlTp061cuygD4Q0AAA+mxiYiInTtQkl7cZcVk3RxgQAhoAQJ+VUjI7+/0MDz+V5M18MpNWk7yZ4eGnMjPzkgYhMAB0cQQAaIh2u52LF5/P8nKysvKL29dBm5l5SYt9OEK02QcAOCRqrVlYWMj58+czNzeX8fFxM2dwxGizD/dJrTVzc3N54403Mjc3Fx8+ALAfrVYrQ0NDOX/+fJLk/PnzGRoaSqvV6m9hQM+YQYN79MkylE7r47W19yxDAQBgV/ta4lhKOZtkOsk7Sc4lma21Lm8z9sUkLyRZTvJWkou11qubxghoHHrtdjvT009ndfW1dC4qWtLZyH05w8NP5dKlHwhpAABsab8B7e1a6xe7348kebXW+uQ2Y5+ttc7u8ngCGqm1Zn5+PktLSxkdHc3ExMShWV9fa83p01/I4uLjSf5yixHfzNjYW7l27d1D8zMBANA797wHrZRyLsn76/e7M2cXDqKg9Zs11YOn3W7nzJnHMjn5dP7sz/7PTE5+K2fOPJZ2u93v0vZkfn4+t24NJXktnVmzja2Qa5LXsrycLCws9KtEAAAaptVq3c5AO9mtScjZdJYrbvR+N7htZaSUMl1KuVBKebE743aHWuvtm4A2WNaXBl6//nJWVn6e3/zmR1lZeTfXr7+c6emnD0VIW1paytDQI+ksa9xKZ0/ajRs3elkWAAAN1mq1bmegnRzb5XE+e5fPe3t/Winl/SR/l+SLd/kYHFG11jz77J9ndXUyyR/ecX519Zu5ePH5XLv2R41eGjg6Opq1tffSmS3bqs6atbVf5tSpUz2uDACAw263GbT3k2yeBds2tG1sHlJrfSfJue1m0Rg8R2Vp4MTERE6c6DQE2drljIwk4+PjvSwLAIAjYLeAdjVbBLJu+PodpZRzpZS3txi7ZcdHBs9RWRpYSsns7Pdz7Nh/TOdnWf95Ot8fO/YfMzPzUqNnAQEAaKYdA9rmINZtuX9l4/0NM2RXk8xsOHchyaWDK5XD7neXBm7l8CwNnJqayo9//F8zNvZojh9/NL/3e/9zjh9/NGNjj+bHP/6vWuwDAHBP9tJm/1w6nRvvuA5aKeX1JJfXW+t3Q9nZ7j/9fJK/2DyDps3+4PqkPf0rSb7WPbp+/bAkeTNjY88dqvb0tdYsLCzkxo0bOXXqVMbHxw9N7QA0VK3JzflkdSkZHk1OTiReW+BI2dd10O5DMQLaAGu323niif+Ujz76/+44d+zYvzb7BMBgW2wnP7uYfLiczkKnteTBkeTxmeRhr49wVAhoNEq73c7Fi89neTkZGnoka2u/zMhIMjPzknAGwOBabCc/mU4+Xr3z3APDyVcuCWlwRAhoNI6lgQCwQa3Jj04nq4vbj3loLHnimuWOcATsFNB2uw4a3BellExMTPS7DABohpvzyW9v7Tzmw+Xk5kLyOa+fcJTt1mYfAID7bXUpu78tG0pWm30pGmD/BDQAgH4bHk2ytsugtWS4+ZeiAfZHQAMA6LeTE8mnTuw85sGR5OR4b+oB+kZAAwDot1KS8dlOt8atPDDcabWvQQgcebo4AgA0heugwUDQZh8A4LCotdOy1UtgAAAFQklEQVStcfVGZ8/ZyXEzZ3DECGgAAAANsVNAswcNAACgIVyoGgCAA1Nrzfz8fJaWljI6OpqJiYkUSzRhzwQ0AAAORLvdzsWLz2d5uWRo6JGsrb2XkZFkZualTE1pcgJ7YQ8aAAD71m63Mz39dFZXX0symaQkqUkuZ3j4qVy69AMhDbo0CQEA4L6pteb06S9kcfHxJH+5xYhvZmzsrVy79q7ljhBNQgAAuI/m5+dz69ZQktfSmTVb/zB+/fvXsrycLCws9KtEODQENAAA9mVpaSlDQ4+ks6xxK509aTdu3OhlWXAoCWgAAOzL6Oho1tbeyyczZ5vVrK39MqdOneplWXAoCWgAAOzLxMRETpyoSZ5KZxZtfSZt/funMjKSjI+P96tEODS02QcAYF9KKZmd/X5e/e438so3P5eRh/4la3UoQ2Utyx/8q/znv/y/8sx3/kqDENgDXRwBANi/xXY+/r+/ngfy4R2nPs6DeeB/eiN5WJt9SLTZBwDgfqo1+dHpZHVx+zEPjSVPXEvMooE2+wAA3Ec355Pf3tp5zIfLyU1t9mE3AhoAAPuzupTd31YOJava7MNuBDQAAPZneDTJ2i6D1pJhbfZhNwIaAAD7c3Ii+dSJncc8OJKc1GYfdiOgAQCwP6Uk47PJA8Nbn39gOHl8RoMQ2ANdHAEAOBiL7eRnFzsNQTKUZK0zc/b4jBb7sIE2+wAA9EatnW6Nqzc6e85Ojps5g00ENAAAgIZwHTQAAIBDQEADAABoCAENAACgIQQ0AACAhhDQAAAAGkJAAwAAaIhdA1op5Wwp5YVSyoXu15GDGAvAnVqtVr9LAIAj5bC9tu5lBu31Wuv3aq1XkswmefWAxrJPf/u3f9vvEnrmMP6sTa25n3X1+rnv9/Pdj8f/7ne/e+CPydHX1L83h9Eg/bc8jD9rE2vud029fP7D+LqaHL7X1h0DWinlXJL31+/XWpeTXNjvWA5Gv/8g9NJh/FmbWrOAdngeH/bK7+LBGaT/lofxZ21izf2uSUA7ekqtdfuTpUwn+ZNa65Mbjv1Tkidrre/cy9hSyvZPCAAAMABqrWWr48d2+XefvYvn2NPY7QoBAAAYdLvtQXs/yeZGH9sFsbsZCwAAwCa7BbSr2SJkbV7eeA9jAQAA2GTHgLbFPrOzSa5svL/eSn+3sQDcu1LKSCllpt91AMBRUEp5u5RyuZTyQr9r2WzHJiHJ7e6MF5K8k+Rcktluh8aUUl5PcrnWOrvbWADuXbcR08Va62S/awGAw66UcqF7abDG2TWgAXCwSimXNwet7qqD6Wz9YdiFJG+lc61JAQ0ANrmH19YXuscna63f7nW9OxHQAHqkG7TOJpnZ3NG2lPJ2rfWL3e9Hkrxaa32y++KSWuvVrV58AGCQ3ctr66YxLyRZXl8R2AS7NQkB4IDUWq9s9QLQXR7+/oZxy+ksF086n/id7S5x/Gx3LACQe3ttLaVMd4PduvfTIAIaQP+dTbJ5v+77pZRztdZLG9bIb76UCQCwtW1fW9NtZLge0mqtl3pc2452u1A1APffrteM7L54NOoFBAAabNvX1u5s2vqHn41rFGIGDaD/3s+ds2O7hjYAYFuH9rVVQAPov6vZ4kVj8/UlAYA9O7SvrQIaQJ9tfrHodm5s3JILADgsDvNrqz1oAD3S3Zi83kHqxSSXNzQAeWbDNVnOJXmmP1UCwOFxFF9bXQcNAACgISxxBAAAaAgBDQAAoCEENAAAgIYQ0AAAABpCQAMAAGgIAQ0AAKAhBDQAAICGENAAAAAaQkADAABoiP8f+Dk0gxFG9YEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "plt.errorbar(eff_wavs, fluxes, yerr=fluxerrs, lw=1.0, linestyle=\" \", capsize=3, capthick=1, color=\"black\")\n",
    "plt.scatter(eff_wavs, fluxes, s=75, linewidth=1, facecolor=\"blue\", edgecolor=\"black\")\n",
    "plt.scatter(eff_wavs, model_photometry, s=75, facecolor=\"orange\")\n",
    "plt.xscale(\"log\")\n",
    "plt.xlim(1000.*(1 + redshift), 60000.*(1 + redshift))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whilst it may not be obvious if you're new to this, once you've spent a long time looking at these kind of plots you'll be able to tell that it looks like the redshift is pretty good, but we've underestimated the age and the mass. Let's try again by stepping up by 0.5 dex in mass and 5 in our age index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi-squared value: 164.328659021\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAE9CAYAAAB+ykFQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X9s3OedJ/b3QyneUMdETJQrSpnSAcqukSiHA6rEpALkusBBTHEstm4AOe3hAnfjRaT7qy4SwEEPe/AYKLD1Ag6wwR5Qyj3f1rj9IxshG2OBQc/yLlo0QEg6MVpck7v4LgJOa5KH3kqhsloz68h8+seMbJriL0nkzHfI1wsghjPzaL4fyjK/857neT7fUmsNAAAA/TfU7wIAAADoENAAAAAaQkADAABoCAENAACgIQQ0AACAhjiQAa2UMtrvGgAAANY7cAGtlHImyQ/XPXaulHKmlHKhT2UBAAAcvIBWa309yfKd+93Adqr7+NVSyqm+FQcAABxojQ9opZQrGzx2qpTydHfm6+kHWbLYDWZT3eOcqrVefZB6AQAA7tfhfhewmVLKuSSnkpzb4Olv11o/3R33gyQvJHm8e3+jZYo3aq2XNznOqSTfSnI1yXOllB90QxsAAEBPlVprbw9YyoVa66V1j40mObdRiCql1FprWXP/TJLnaq1Tax77Wa31I/dQww/XBLynk1yqtS5367hQa/3de//JAAAAHkw/ljj+oJTy3J073VD0QpJXd/jnT2XNHrKuG93gtq07M3Pd2yS5nOSL3fufSXJp0z8MAACwh3q+xLHW+nopJd2Q9jtJnkvylVrr+tC1mY8+4PFfTfKRNfevRigDAAAaoC970Loh7aPptLv/9D2EsyS5kWR9U5AHCm0AAABN0Jcujt1ljY8n+XqS//Ee//jVbBDINPYAAAAGXb/a7D+X5OvdpiDfWrsnbTvrg1i3C+NO968BAAA0Vj+6OD6X5HfWLmvsNvj4zNrujt3HzqUT5n43yZXu/rG1z72e5Ey6XRh791MAAADsvp4HNAAAADbWryWOAAAArNPzgFZKMWUHAAAcWFtlor602RfSAAAA7taXJY61Vl+78PWVr3yl7zX4WQev5n7W1etj7/Xx9uL1/Y70dT9fTf19M4hfB+nvchB/1ibW3O+aenn8QTyv1trMc+tW7EEbYL/xG7/R7xJ6ZhB/1qbW3M+6en3svT5eU/8bc/D4t7h7DtLf5SD+rE2sud819fL4zqu90Y82+7XXxwQYFKWUbT9ZAwB2ronn1m5NZaPnzKABNMgzzzzT7xIAYF8ZtHOrGTQAAIAeMoMGAAAwAAQ0AACAhhDQAAAAGkJAAwAAaAgBDQAAoCEENAAAgIYQ0AAaqNVqpZRy11er1ep3aQDAHnIdNICG614rpd9lAAC7xHXQAHaJmS0AYC+ZQQO4T72a2TKDBgD7ixk0AACAASCgAQAANISABgAA0BACGgAAQEMIaAANVWvN7OxskmR2dlajEAA4AAQ0gAZqt9s5efJ0pqaeTJJMTX05J0+eTrvd7nNlAMBeEtAA7tFez2y12+2cP/9k3nzz93Lr1o+SJLdu/Thvvvl7OX/+SSENAPYx10EDuAftdjsXL34ty8slt27964yMfCKjo8nMzPOZnp5+4NevtebEiU9mYeGbST7ffbQkufN785WMjz+Va9d+nFI2vHwKANBwroMGsAt6MbM1NzeXmzeHkkxtMmIqy8vJ/Pz8Ax8LAGgeAQ1gB2qtuXDhq1lZeSmdma07H3qVJJ/PyspLuXjxaw+83HFpaSlDQ4+sef31SoaGHsni4uIDHQcAaCYBDWAHejWzNTY2ltXVn+S9JY3r1ayuvpHjx48/0HEAgGYS0AB2oFczW5OTkzl6tCa5ssmIKxkdTSYmJh7oOABAMwloADvQq5mtUkouXfpGhoefSPLKmuPVJK9kePiJzMw8r0EIAOxTAhrADrx/Zqtm8ldn819/Jpn81dl0wtPuzWxNT0/n8uUXMz7+VEZGTidJRkZOZ3z8qVy+/OKudIsEAJpJm32AHWq323nh2X+Yb37pcEaP/CKr9VaGykiW3/pg/vt/cTtfeeYPdzU81VozPz+fs2fPZnZ2NhMTE2bOAGAf2KrNvoAGsFML7bzzf34hh/L2XU+9k4dy6Nf/OHl492e3ur/Ed/11AYD+ENAAHlStyXdPJCsLm485Mp48di3Z5VkuAQ0A9hcXqgZ4AK1WK2cfGcpf3tginCXJ28vJdReQBgDunxk0gJ249p3k+08k7/zV5mMOfzj57B8kJ76wq4c2gwYA+4sZNIAHsdBOXvtHW4ezJMlqMuwC0gDA/Tvc7wIAGm2hnXzvfPLOyvZjHxpNjrmANABw/7adQSulnCmlXCilPF1K+XYp5dQWY58rpdRSys9KKVe2GgvQeLUm8xd2Fs4ODSePzuxag5BWq5VSyrtt9e9832q1duX1AYBm2nIPWillNMkXa62XuvfPJZmptX58k/EX7ozd4jXtQQMGw1/MJn82ldy+tfW4X/lPkrP/fE9a7AMA+8+D7EE7leTra+7/IMmpbnB7oIJ8Ggw03spStvs1eesXyRf+5/8vrRd0bwQANrd+dcxmtu3iWEo5VWu92v3+XJJv11o/ssnYp5NcTbKcZCrJ79Ral9eNMYMGDIadzKAdHkn+3qvJxyZ7VxcAMNAeqIvjnXDWdTHJV7YYfqnWernW+mqSbyX503uqFKBJjk0mHzi69ZgD3hhk7aeBVkcAwIPb8XXQSikXktyotV7e8YuXUpN8ZO0smhk0YJC89vKz+dTPWjny0N3PvfXXyY8+2sqjjz3T+8IayPXaAGBntppB21FA6y5tTHdmbLMxZ5K8UGv99JrH6voDC2jAoKi15sSJT+bvfOxLmfmtmYweWc5qvZWhMpLlt0Zz8Z9dzL+6/oe5du3H264nPwgENADYmQcKaN3glVrr69375+/MonXb6N+otS5v0vHxYq318XWvJ6ABA2F2djZTU0/m1q0fJUkmPj6f4x85m8WfzWb+p51ljSMjp/Pqq3+Qyck92INWa3J9rtOsZHiss+SywUFQQAOAndkqoG15oepuAPvTJKNrPh2+muTOMsfnklxJZ+/ZcinlancpZJJ8PFvvVwNotKWlpQwNPZKk8/tv/qd3Qth7YWxo6JEsLi7u/sEX2slrF5O3l9PZLrza2e/26Ix2/gCwj20Z0LoNQjbs2Nh9/vF19zddAgkwaMbGxrK6+pMkNXdC2vvVrK6+kePHj+/ugRfayffO332B7Nu3Oo9/7rKQBgD71LZdHAEOqsnJyRw9WtNZKLCRKxkdTSYmdrGLY63J/IW7w9kd76x0ZtYsJQSAfUlAA9hEKSWXLn0jw8NPJHklnZm0dG9fyfDwE5mZeX53G4Rcn0t+eXPrMW8vJ9ebc2HsWmtmZ2eTdPbt2YcGAPdPQAPYwvT0dC5ffjHj409lZOR0kk5jkPHxp3L58ouZnt7lpYYrS9n+V/NQsrIH+97uQ7vdzsmTpzM19WSSZGrqyzl58nTa7XafKwOAwbTj66Dt2gF1cQQGUK018/PzOXv2bGZnZzMxMbE3rfX/Yjb5s6nOfrPNHB5J/t6rycf2oHPkPWi32zl//smsrLyUZCrvNjPJlQwPP7E3ARYA9oEHvg7aLhcjoAEDa89bydeafPdEsrKw+Zgj48lj1/racv/ONeIWFr6Z5PPdR0veWwb6SsbHn3KNOADYwFYBzRJHgCYpJZm4lBwa3vj5Q8OdVvt9Dj1zc3O5eXMonZmzjUxleTmZn2/OXjkAGAQCGkDTPDzdaaV/ZLyznPHwhzu3R8Yb02J//TXi7lb27hpxALCPbXkdNAD65OHp5Pi1TrfGlcVk+HhybKLvM2d3vP8accnkr85lbDRZWp7N3L/r7I3bk2vEAcA+Zw8awD3Y8z1oA+LOHrS/87EvZea3ZjJ6ZDmr9VaGykiW3xrNxX92Mf/q+h/agwYAG7AHDeABtFqtlFLeDRp3vm+1Wv0trI9KKfnjf/oPcvl/+Cc5cezNfGj4Vo4eST40fCsnjr2Zy0/9k3zn9/9b4QwA7pEZNADu3YB0mwSAJjKDBsDuuj6X/PLm1mPeXu7soQMAdkxAA+DerSxl+1PIUKfBCQCwYwIaAPdueCzJ6jaDVjvdJwGAHRPQALh3xyaTDxzdesxDo51LAwAAOyagAXDvSkkmLiWHhjd+/tBw8uiMBiEAcI8ENADuz8PTyecud7o1Hh7JzbeSHB7p3P/c5c7zAMA9EdAAuG+tF+ZTvvBmJv/xrfzmTDL5j2+lfOHNtF7QvREA7ofroAEAAPSQ66ABAAAMAAENAACgIQQ0AACAhhDQAAAAGkJAAwAAaAgBDQAAoCEENAAAgIYQ0AAAABpCQAMAAGgIAQ0AAKAhBDQAAICGENAAAAAa4nC/CwCAxqo1uT6XrCwlw2PJscmklH5XBcA+JqABwEYW2slrF5O3l9NZcLKaPDSaPDqTPDzd7+oA2KdKrbW3Byyl9vqYAHBPFtrJ984n76zc/dyh4eRzl4U0AO5bKSW11g2XZAhoALBWrcl3TyQrC5uPOTKePHbNckcA7stWAU2TEABY6/pc8subW495ezm5Pt+begA4UAQ0AFhrZSnbnx6HkpXFXlQDwAGzbZOQUsqZJJ9JMprk0SRfr7Ve3WTsqSTnk7ye5EySS7XW5d0rFwD22PBYktVtBq0mw8d7UQ0AB8yWAa2UMprkM7XWS93755JcSfLxTf7It2utn+6O/UGSF5I8vnvlAsAeOzaZfOBocvvW5mMeGk2OTfSuJgAOjO3WcJxK8vU193+Q5FQ3uL1Pd6btxp373Zmzc7tRJAD0TCnJxKVOt8aNHBrutNrXIASAPbBlQKu1vp5kas1Dn0myvMmyxVNJ1j9+oxvc3qeU8u5Xq9W6x5IBYI89PN1ppX9kPDk8khz+cOf2yLgW+wDcl1ar9W4G2so9tdkvpXw7ybdqrZc3eO5Ckqla6+NrHvtpkse7Qe/OY9rsAzAYau10a1xZ7Ow5OzZh5gyAB7ZVm/1tm4SseZEL2SScdd1Ip5HIWh/d6esDQOOUknxsst9VAHCA7CigdZuDXK21vrrFsKvZIJCtnT0DAABgc9teB+1O84874ayUcn7Nc+82DFkfxLot97cKdAAAAKyx5R60bsj6Yd6/dPFqrfXj3ee/neTKmjb8Z9Lp3LjpddDsQQMAAA6yrfag3VOTkF0qRkADAAAOrK0C2rZLHAEAAOgNAQ0AAKAhBDQAAICGENAAAAAaQkADAABoCAENAACgIQQ0ANhAq9VKKeWur1ar1e/SANjHXAcNGAitVivPPvvsXY8/88wz3jCz57rXq+l3GQDsEy5UDewr3izTa/7NAbCbXKgaAO5DrTWzs7NJktnZWSENgD0noAEDw5tleqndbufkydOZmnoySTI19eWcPHk67Xa7z5UBsJ8JaMBA8GaZXmq32zl//sm8+ebv5datHyVJbt36cd588/dy/vyT/t0BsGfsQQMa786b5ZWVl5JMpfPZ0mqSKxkefiKXL7+Y6enp/hbJvlFrzYkTn8zCwjeTfL77aEly59z1SsbHn8q1az9OKRtuHwCALdmDBgysWmsuXPhqN5x9Pp03yunefj4rKy/l4sWvWe7Irpmbm8vNm0PpfBiwkaksLyfz8/O9LAuAA0JAAxrNm2V6bWlpKUNDj+S9DwPWKxkaeiSLi4u9LAuAA0JAAxrNm2V6bWxsLKurP8l7SxrXq1ldfSPHjx/vZVkAHBACGtBo3izTa5OTkzl6tCa5ssmIKxkdTSYmJnpZFgAHhIAGNJo3y/RaKSWXLn0jw8NPJHkl7304UJO8kuHhJzIz87wGIQDsCQENaDRvlumH6enpXL78YsbHn8rIyOkkycjI6YyPP6VrKAB7Spt9YCC02+1cvPjV/NpHfpGjv/Lvc/Ov/1b+7c8+mJmZb3izzJ6ptWZ+fj5nz57N7OxsJiYmfBgAwAPbqs2+gAYMhoV26vzFrP7iem791UpG/sZwhj54LGViJnlYQGNvdU+k/S4DgH1CQAMG20I7+d755J2Vu587NJx87rKQxp4S0ADYTQIaMLhqTb57IllZ2HzMkfHksWuJpWfskX0X0GpNrs8lK0vJ8FhybNL/PwA9tFVA0yQEaLbrc8kvb2495u3l5LoLVbO7Wq1WSinv7jm7832r1epvYQ9qoZ28fDL5s6nk+7/ZuX35ZOdxAPrODBrQbH/+x503kbd/vvmYwx9OPvsHyYkv9KoqGEyWCwM0ghk0YHANjyVZ3WbQajLsQtWwpVqT+Qsbh7Ok8/hrFzvjAOgbAQ1otmOTyQeObj3modHkmAtVw5YsFwYYCAIa0GylJBOXOsuvNnJoOHl0RoMD2M7KUrY/7Q8lK4u9qAaATQhoQPM9PN3ZG3NkPDk80tlzdnikc9+eGdgZy4UBBoImIcDgqLWz/GplsfMm8tiEmTPYKZesAGgM10EDAHRxBGgIAQ0A6Fhod7o1vr2czk6H1U6jnUdnhDOAHhHQAID3WC4M0FcCGgAAQEO4UDUAAMAAENAAAAAaQkADAABoiB0FtFLKlR2Mea6UUkspPyulXCmlnHrw8gAAAA6Ow1s9WUo5l+RUknM7eK2fbrbRDQAAgO1tOYNWa3211nqpV8UAAAAcZLu5B220lHK+lHKuu9xxdBdfGwAAYN/bzYB2qdZ6udb6apJvJfnTzQaWUt79arVau1gCAABA87RarXcz0FZ2dKHq7sWl72l/WSmlJvlIrXV5g9e6l5cCAADYN/b8QtWllDOllB+uf3x9OAMAAGBz9x3QSimn1uwzu5pkZs1z55JcfsDaAAAADpTt2uyfSbfFfinluSRXunvMkuS5JFfS2Xu2XEq5Wkq50H3u40m+skc1AwAA7Es72oO2qwe0Bw0AADjA9nwPGgAAAA9OQAOAA2Rtm2eXvQFoHkscAeCA6i6x6XcZAAeOJY4AwLtqrZmdnU2SzM7OCmkADSKgAcAB0m63c/Lk6UxNPZkkmZr6ck6ePJ12u93nygBILHEEgAOj3W7n/Pkns7LyUpKpdD6nXU1yJcPDT+Ty5RczPT3d3yIBDoCtljgKaABwANRac+LEJ7Ow8M0kn+8+WpLcOSe/kvHxp3Lt2o9TyobvGQDYJfagAcABNzc3l5s3h9KZOdvIVJaXk/n5+V6WBcA6AhoAHABLS0sZGnoknVmzjZQMDT2SxcXFXpYFwDoCGgAcAGNjY1ld/UneW9K4Xs3q6hs5fvx4L8sCYB0BDQAOgMnJyRw9WpNc2WTElYyOJhMTE70sC4B1BDQAOABKKbl06RsZHn4iySt5byatJnklw8NPZGbmeQ1CAPpMF0cAOEDa7XYuXvxalpeTW7f+TUZGPpHR0WRm5nkt9gF6RJt9AOBdtdbMz8/n7NmzmZ2dzcTEhJkzgB7SZh8ASJK0Wq0MDQ3l7NmzSZKzZ89maGgorVarv4UBkMQMGgAAQE+ZQQMAABgAAhoAAEBDCGgAAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAAAA0hoAEAADSEgAYAANAQAhoAAEBDCGgAAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEPsKKCVUq7sYMypUsrTpZRz3dvRBy8PAADg4Ci11s2fLOVcklNJZmqtZcsXKuWHtdZPd78fTfJCrfXxDcbVrY4JAACwn5VSslm+2jKgrXmBulVAK6WcSfJcrXVqzWM/q7V+ZJPX2lnlAAAA+8xWAW239qCdSrK87rEb3eAGAADADuxWQPvovQwupbz71Wq1dqkEAACAZmq1Wu9moK3sVkC7kWR9U5BNQ1ut9d0vAQ36784vjLO/VvKFRzu3PkABANg9rVbr3Qy0lcO7dLyr2SCQ1Vpf36XXB/bQf/mflVz4p4fzoQ++k9XVmkNDJT//xaEsPLz1JzwAAOyu+55B67bVH03uDmKllFNJXn3A2oAeeO3lZ/Opn7VyfPR2PvTBmqNHkpEP1hwfvZ1P3WjltZef7XeJAAAHxnZt9s8kOZfkuSS/m+RKrfXV7nPf7t6/tG7s60nOJLlUa13fOEQXR2iQurqapf/lV3J89PamYxaXD2fsH/11ypDr2gP0RK3J9blkZSkZHkuOTSbb7FkBBssDt9nf5WIENGiI//f/+F/zt65eyIc+uPn/k3/5i5J///EX8rd//bd6WBnAAbXQTl67mLy9nM5Cp9XkodHk0Znk4el+Vwfskl602QcG0K2/+LdZrYe2HLNaD+XWf3yjRxUBHGAL7eR755O33kxu30pu/7xz+9abnccX2v2uEOgBAQ0OsJGP/VoOlXe2HDNU3snI33ykRxUBHFC1JvMXkndWNn7+nZXOzJpVSLDvCWhwgH3qP38yP//F1jNof/mLQ/nU3/1yjyoCOKCuzyW/vLn1mLeXk+vzvakH6BsBDQ6wMjSUhYd/O2+9vfHzb/11svDwb2sQArDXVpay/duyoWRlsRfVAH3kXRcccI8+9kx+///+9fz59eQvV5Kbb3Vu//x68vv/z6/n0cee6XeJAPvf8FiS1W0GrSbDx3tRDdBHujgCSTot93/0f/3z3PqPb2Tkbz6ST/3dL5s5A+iVWpPvnkhWFjYfc2Q8eeyalvuwD2izDwDQdHe6OG7UKOTQcPK5y1rtwz4hoAEADALXQYMDQUADABgUtXa6Na4sdvacHZuwrBH2GQENAACgIbYKaDoAAAAANISABgAA0BACGgAAQEMIaAAAAA0hoAEAADSEgAYAANAQAhoAAEBDCGgAAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAADdFqtVJKueur1Wr1uzSgR0qttbcHLKX2+pgAAIOi1pq5ubl89rOfzfe///1MTk6mlNLvsoBdVEpJrXXD/7HNoAEANES73c7Jk6czNfVkkmRq6ss5efJ02u12nysDesUMGgBAA7Tb7Zw//2RWVl5KMpXO5+irSa5kePiJXL78Yqanp/tbJLArtppBE9AAAPqs1poTJz6ZhYVvJvl899GS5M57plcyPv5Url37seWOsA9Y4ggA0GBzc3O5eXMonZmzjUxleTmZn5/vZVlAHwhoAAB9trS0lKGhR9KZNdtIydDQI1lcXOxlWUAfCGgAAH02NjaW1dWf5L0ljevVrK6+kePHj/eyLKAPBDQAgD6bnJzM0aM1yZVNRlzJ6GgyMTHRy7KAPhDQAAD6rJSSS5e+keHhJ5K8kvdm0mqSVzI8/ERmZp7XIAQOAF0cAQAaot1u5+LFr2V5Obl1699kZOQTGR1NZmae12If9hFt9gEABkStNfPz8zl79mxmZ2czMTFh5gz2GW32AQAGQKvVytDQUM6ePZskOXv2bIaGhtJqtfpbGNAzZtDgAdRaMzc3l6WlpYyNjWVyctKnnAAAbGmrGbTDO/jDp5KcT/J6kjNJLtValzcZ+1ySp5MsJ/lBkou11qv3Wzg02Xv7BDrXplld/Yl9AgAAPJBtZ9BKKT+stX66+/1okhdqrY9vMvZCrfXSNq9nBo2Bn3lqt9t57LEv5vbtv7rrucOH/0ZefvmPhDQAADZ033vQSilnkty4c787c3ZuNwq682VN9cHTbrdz8uTpTE09md/8zf8tU1NfzsmTp9Nut/td2o7UWnPhwldz+/Z30ml/vLYVcs3t29/JxYtfiw8iAAC4o9VqvZuBtrLlDFop5XyS/2btjFkp5adJHq+1vr7B+KeTXE1nieNUkt9ZvxzSDNrB1m63c/78k1lZeSmdfyIlnWBzJcPDT+Ty5RcbP/M0Ozubqaknc+vWj9KpP3nv50iSmpGR03n11T/I5ORkf4oEAKCxHmQP2kfv8Vjv7k8rpdxI8qdJPn2Pr8E+dWfmaWVlKsl/cdfzKytfysWLX8u1a3+/0csdl5aWMjT0SN4LZ+t19qQtLi72siwAAPaB7drs30gyuu6xTUPb2tmy7gzbme6+Ncjc3Fxu3hxK8lI2WhqYvJTl5WR+fr5fJe7I2NhYVld/kvfqX69mdfWNHD9+vJdlAQCwD2wX0K5mg0C2yfLGM6WUH24wdsOOjxw8+2XmaXJyMkeP1iRPpPOzrF3mWJI8kdHRZGJiol8lAkD/1Jr8xWzy53/cubW1Be7Jlksca62vr11q1m25/+q6+ze6Iexqkpk1z51Lcnm3C2ZwvX/maaOQNhgzT6WUXLr0je5euv89k7/64YyN/ocsLf+nmft3P8/w8H+XmZkXG71MEwD2xEI7ee1i8vZyOvMAq8lDo8mjM8nDzd5jDk2xkzb7Z9Lp3HjXddBKKd9OcuVOa/1uKDvV/aMfjyYhrFFrzYkTn8zCwjeTfL776NrmGq9kfPypXLv244EIN6+9/GweXvif8qEPvpPVeiiHyjv5+S8OZeHh386jjz3T7/IAoLcW2sn3zifvrNz93KHh5HOXhTTo2qpJyLYBbQ+KEdAOsH1z/TAnIQB4T63Jd08kKwubjzkynjx2LRmAD2Fhr933ddBgt01PT+fll/8o4+OfyMjIJ/LhD/9XGRn5RMbHPzE44azWZP7CxuEs6Tz+2kVr7gE4OK7PJb+8ufWYt5eT681uBAZNsF2bfdh109PTuXbt72d+fj6Li4s5fvx4JiYmBmJZY5J7Owl9zHXQADgAVpay/ef+Q8lKsxuBQRMIaPRFKWVwL+LsJAQA7zc8lmR1m0GryXCzG4FBE1jiCPfKSQgA3u/YZPKBo1uPeWg0OeYSNLAdAQ3ulZMQALxfKcnEpU6jrI0cGu602h+U7QzQRwIa3CsnIQC428PTnS7GR8aTwyPJ4Q93bo+M624M90CbfbhfLsYJAHertdMoa2Wxs9z/2IQPLWEd10GDveIkBADAPRLQAAAAGmKrgKbNPv1Ra+d6YitLna6IxybNPAEAcOAJaPSevVsAALAhSxzprYV28r3zyTsrdz93aFiXJwAA9j170GiGWpPvnkhWFjYfc2Q8eeya5Y4AMKBqrZmbm8vS0lLGxsYyOTmZ4rwO72MPGs1wfS5GDCSeAAAEPElEQVT55c2tx7y93OmK+LHJ3tQEAOyadrudixe/luXlkqGhR7K6+pOMjiYzM89netoKGdgJAY3eWVnK9tdGH+q0rAcABkq73c5jj30xt2//VfeRf50kuXUreeyxL+bll/9ISIMd2O7dMuye4bEkq9sMWu1cTwwAGBi11ly48NXcvv2dJLX7lXe/v337O7l48WuxzQW2J6DRO8cmkw8c3XrMQ6Odiz0DAANjbm4uN28OJZnaZMRUlpeT+fn5XpYFA0lAo3dKSSYudbo1buTQcKfVvo3EADBQlpaWMjT0SJLNzuGdPWmLi7YxwHYENHrr4elOK/0j48nhkeTwhzu3R8a12AeAATU2NpbV1Z/kvaWN69Wsrr6R48dtY4DtaLNPf9Ta6da4stjZc3ZswswZAAyoWmtOnPhkFhYeTfIvNhjxpYyP/yDXrv1Yy33I1m32zaDRH6V0Wumf+ELn1i9rABhYpZRcuvSNDA9fSfIv02kKVru3/zLDw1cyM/O8cAY7YAYNAIBd8d510NK9DtobroMGG9hqBk1AAwBg19RaMz8/n8XFxRw/fjwTExNmzmAdAQ0AAKAh7EEDAAAYAAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAN0mq1+l0CAOwrg3ZuFdAG2J/8yZ/0u4SeGcSftak197OuXh97r4+3F6//7LPP7vprsv819ffNIDpIf5eD+LM2seZ+19TL4w/ieTUZvHOrgDbA+v0LoZcG8Wdtas0C2uC8PuyUf4u75yD9XQ7iz9rEmvtdk4C2/5Raa28PWEpvDwgAANAwtday0eM9D2gAAABszBJHAACAhhDQAAAAGkJAAxgApZTRUspMv+sAgP2glPLDUsqVUsrT/a5lvcP9LgCAHTmX5FS/iwCAfeLrtdZX+13ERsygAfRYKeXKBo+dKqU8XUo5170dXfPcuSSNPIkAQBPc67k1yZnu48/1sMwd0cURoEe6QetUkpn1rXVLKT+stX66+/1okhdqrY+XUk4lSa31ainlSq11queFA0BD3c+5dd2Yp5Ms11ov9arm7ZhBA+iRWuurG50ASilnktxYM245nSWNSXImyalSyvkkH+2OBQByf+fWUsr5brC740YaREAD6L9TSZbXPXajlHKm1np5zRr50QAAO7HpuTXdbQN3Qlqt9XKPa9uSJiEA/ffR7QZ0Tx6NOoEAQINtem7tzqbd+fCzcXu8zaAB9N+N3D07tm1oAwA2NbDnVgENoP+uZoOTRq319T7UAgD7wcCeWwU0gD5bf7Lodm5s3JILABgUg3xutQcNoEe6G5PvdJB6LsmVNQ1AvtJt9ft6Op0bv9KfKgFgcOzHc6vroAEAADSEJY4AAAANIaABAAA0hIAGAADQEAIaAABAQwhoAAAADSGgAQAANISABgAA0BACGgAAQEMIaAAAAA3x/wM4G/jJDZcakwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mass = 10**10.5  # Solar masses\n",
    "redshift = 1.\n",
    "age_index = 30\n",
    "\n",
    "# Our photometry is in erg/s/cm^2/A per Solar mass of stars\n",
    "model_photometry = get_model_photometry(redshift, age_index)\n",
    "\n",
    "# Convert to erg/s/cm^2/A by multiplying by the total mass we want.\n",
    "model_photometry *= mass\n",
    "\n",
    "diffs = fluxes - model_photometry\n",
    "\n",
    "chisq = np.sum(diffs**2/fluxerrs**2)\n",
    "\n",
    "print \"Chi-squared value:\", chisq\n",
    "\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.errorbar(eff_wavs, fluxes, yerr=fluxerrs, lw=1.0, linestyle=\" \", capsize=3, capthick=1, color=\"black\")\n",
    "plt.scatter(eff_wavs, fluxes, s=75, linewidth=1, facecolor=\"blue\", edgecolor=\"black\")\n",
    "plt.scatter(eff_wavs, model_photometry, s=75, facecolor=\"orange\", zorder=10)\n",
    "plt.xscale(\"log\")\n",
    "plt.xlim(1000.*(1 + redshift), 60000.*(1 + redshift))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're already doing a lot better, but it's still going to take forever to get anywhere this way. Let's go to the next step and iterate over a grid of points in order to make things faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}