Created
May 20, 2013 12:09
-
-
Save aasensio/5611869 to your computer and use it in GitHub Desktop.
Faraday rotation IPython Notebook
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"metadata": { | |
"name": "Faraday rotation" | |
}, | |
"nbformat": 3, | |
"nbformat_minor": 0, | |
"worksheets": [ | |
{ | |
"cells": [ | |
{ | |
"cell_type": "heading", | |
"level": 1, | |
"metadata": {}, | |
"source": [ | |
"Rotation measure" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"import numpy as np\n", | |
"import matplotlib.pylab as plt" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [], | |
"prompt_number": 12 | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The relation between the rotation angle and the rotation measure measure is linear:" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"$\\beta = \\mathrm{RM} \\lambda^2$" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Our aim is to compute the rotation measure from a measurement of different rotation angles at different wavelengths. We first define the set of wavelengths where we make the measurement. Let us assume that we observe at " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"cLight = 3e10\n", | |
"frequency = np.asarray([4.4, 4.9, 5.4, 6.7, 8.9, 11.0])*1e9\n", | |
"wavelength = cLight / frequency\n", | |
"ObservedBeta = 0.03 * wavelength**2 + 0.1*np.random.randn(6)" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [], | |
"prompt_number": 13 | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Let's do the plot of what we observe, noting that the points have some random noise added." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"plt.errorbar(wavelength**2, ObservedBeta,yerr=0.2,marker='.',color='k',linestyle='none')\n", | |
"plt.xlabel(r'$\\lambda^2$')\n", | |
"plt.ylabel(r'$\\beta$')" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"output_type": "pyout", | |
"prompt_number": 14, | |
"text": [ | |
"<matplotlib.text.Text at 0x9917f6c>" | |
] | |
}, | |
{ | |
"output_type": "display_data", | |
"png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAESCAYAAADTx4MfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFnlJREFUeJzt3X9M1Pfhx/HXKTYbop4mFcnBhh04geJxrZYljSnVUn+s\nZba6RTutUSzMpbouzbLZLCumnVFns7h16crSmjrjxJi1JSuYuU3UaJS5nd0Su6BVJqCSKV7qr07E\n9/cPv95A7op9C/f+AM9HYsLdvYWX75yfF5/P5z7vj88YYwQAwOc0xHUAAED/RIEAAKxQIAAAKxQI\nAMAKBQIAsEKBAACseKZAmpqa9OijjyovL0/333+/fvGLX8Qct3LlSmVnZysYDCocDic4JQDgliTX\nAW4ZNmyYfv7zn6ugoECXLl3Sgw8+qOLiYuXk5ETH1NTU6Pjx4zp27JgOHTqk5cuX6+DBgw5TA8Dg\n5Zk9kHHjxqmgoECSlJKSopycHJ0+fbrLmOrqai1evFiSVFhYqEgkotbW1oRnBQB4aA+ks8bGRoXD\nYRUWFnZ5vqWlRRkZGdHH6enpam5uVmpqavQ5n8+XsJwAMFDYLErimT2QWy5duqR58+Zp48aNSklJ\n6fb67f/IWIVhjPHUn5dfftl5BjINnExezUWm/pvJlqcKpL29XXPnztXChQs1Z86cbq8HAgE1NTVF\nHzc3NysQCCQyIgDg/3mmQIwxKi0tVW5url544YWYY0pKSrR582ZJ0sGDB+X3+7scvgIAJI5nzoHs\n379fW7Zs0aRJkxQKhSRJa9as0alTpyRJ5eXlmj17tmpqapSVlaXhw4dr06ZNLiPfsaKiItcRuiHT\nnfFiJsmbuch0Z7yYyZbP3M0BMA/y+Xx3dUwPAAYb2+2mZw5hAQD6FwoEAGCFAgEAWKFAAABWKBAA\ngBUKBABghQIBAFihQAAAVigQAIAVCgQAYIUCAQBYoUAAAFYoEACAFQoEAGCFAgEAWKFAAABWKBAA\ngBUKBABghQIBAFihQAAAVigQAIAVCgQAYIUCAQBYoUAAAFYoEACAlSTXAQAA/1NXV6e6ujpJ0o0b\nN3T9+nXdc889KioqUlFRkdNst/MZY4zrEL3J5/NpgP2TAAxSe/fu1Y9//GPt3bu3T3+O7XaTQ1gA\nACsUCADACgUCALBCgQAArFAgAAArFAgAwAoFAgCwQoEAAKx4pkCWLl2q1NRU5efnx3y9rq5Oo0aN\nUigUUigU0quvvprghACAzjyzlMmSJUu0YsUKPfvss3HHPPLII6qurk5gKgBAPJ7ZA5k6dapGjx79\nmWNYogQAvMMzeyA98fl8OnDggILBoAKBgDZs2KDc3NyYYysqKqJfe3EBMgBwqfOCjXfDU4spNjY2\n6sknn9Q///nPbq9dvHhRQ4cOVXJysmpra/W9731PDQ0N3caxmCKAgYLFFHvJiBEjlJycLEmaNWuW\n2tvb1dbW5jgVAAxe/aZAWltbow1ZX18vY4zGjBnjOBUADF6eOQeyYMEC7dmzR+fOnVNGRoZWr16t\n9vZ2SVJ5ebl27NihN954Q0lJSUpOTta2bdscJwaAvlNWVqa//vWv+ve//61IJCK/3+86UjeeOgfS\nGzgHAmAgKCoq0p49eyRJ3/zmN7V9+/Y++1kD/hwIAAwmt875pqSkqLKy0nGa2NgDAQAPikQieuqp\np/Tf//5XBw4c6NOfZbvdpEAAOBfvuoTBfh2X1z/G65mT6AAGr85FceTIEW3YsEFbtmxxGwo94hwI\nAE+5cuWKTpw44ToG7gAFAgCwQoEAAKxQIAAAKxQIAMAKBQIAsEKBAACsUCAAACsUCADACgUCALBC\ngQAArFAgAAArFAgAwAqr8QKAh3Re2v7cuXMaOnSoKioqPLm0PQUCAB7ixaKIhwIB4BllZWU6fPiw\nTp48qUgkIr/f7zoSPgPnQAB4RkNDg8LhsCKRiMrKylzHQQ8oEACekZycLEkaPny4KisrHadBT7gn\nOgDPiEQimjt3ri5evKj6+nrXcQYN2+0meyAAPMPv9+uVV15RUhKnZ/sDCgQAYIWaB/q5ztcNdNaf\nPg6K/okCAfq5zkWxe/duXbt2TTNmzHAbCoMCh7CAAeTAgQPau3ev6xgYJCgQAIAVCgQAYIUCAQBY\noUAAAFYoEACAFQoEAGCFAgEAWKFAAABWPFMgS5cuVWpqqvLz8+OOWblypbKzsxUMBhUOhxOYDgBw\nO88UyJIlS7Rz5864r9fU1Oj48eM6duyYKisrtXz58gSmAwDczjNrYU2dOlWNjY1xX6+urtbixYsl\nSYWFhYpEImptbVVqamqCEgLoK50XhDx9+rTa2tpUUVHBgpAe55kC6UlLS4syMjKij9PT09Xc3Byz\nQCoqKqJf8wYEvI//p4kVbwXnz6vfFIikbnfM8vl8Mcd1LhAAQFe3F/bq1autvo9nzoH0JBAIqKmp\nKfq4ublZgUDAYSIAGNz6TYGUlJRo8+bNkqSDBw/K7/dz/gPopKysTG+99ZaqqqoUiURcx8Eg4JlD\nWAsWLNCePXt07tw5ZWRkaPXq1Wpvb5cklZeXa/bs2aqpqVFWVpaGDx+uTZs2OU4MeEtDQ4NOnjwp\n6WaZbN++3XEiDHQ+c/uJhX7O5/N1O1cCDAazZ89WbW2t0tLSdPToUfn9fteR0E/Ybjc9swcC4O5s\n3bpVU6dO1eOPP055ICH6zTkQAJ/N7/dr/vz5+sIXvuA6CgYJCgQAYIUCAQBYoUAAAFYoEACAFQoE\nAGCFAgEAWKFAAABWKBAAgBUKBABghQIBAFihQAAAVigQAIAVCgQAYKXH5dx//etf69ChQ0pPT1dZ\nWZmqqqo0ZswYzZkzR2PGjElERgCAB/V4Q6na2lrNmjVLH3/8sVatWqXS0lJ9/PHH2rlzp37yk59o\n8uTJicp6R7ihFAaburo61dXVSZL27dun9vZ2TZs2TUVFRSoqKnKaDf1Dn91Q6urVq7px44a+8pWv\naNKkSZoxY4Ykafny5dq4caPnCgT9T+cNYGdsAO9M53lqbW3VjRs3lJaW5jYUBoUe90A+/fRTVVdX\na/z48ZoyZUqX137/+9/r6aef7tOAnxd7IP1bKBTSvn37lJKS4joKMGjYbjfv+J7oZ86cUTgcljFG\nHR0damxsVF5enqZPn/65f2hfokD6t5EjR6q5uVkjR450HQUYNPq8QG4XiUQUDod19epVDRkyRDNn\nzrT5Nr2OAunfKBAg8RJeIF5FgfRvFAiQeLbbTa4DAQBYoUAAAFYoEACAFQoEAGCFAgEAWKFAAABW\nKBAAgBUKBABghQIBAFihQAAAVnpczh19i6XMAfRXrIXlIS+99JLuu+8+LVu2zHUUZ1gLC0g81sIa\nAC5evKirV6+6juFMWVmZrly5onnz5ikSibiOA6AHFAg8o6GhQR0dHdq1a5fKyspcxwHQA08VyM6d\nOzVx4kRlZ2dr3bp13V6vq6vTqFGjFAqFFAqF9OqrrzpIib6SnJwsSXrggQdUWVnpOA2AnnjmJHpH\nR4eef/55/elPf1IgENCUKVNUUlKinJycLuMeeeQRVVdXO0qJvrR161bde++9evfdd+X3+13HAdAD\nz+yB1NfXKysrS5mZmRo2bJjmz5+v999/v9u4/nqCHD3z+/364he/SHkA/YRn9kBaWlqUkZERfZye\nnq5Dhw51GePz+XTgwAEFg0EFAgFt2LBBubm53b5XRUVF9Gs+DgsAXcW7fODz8kyB+Hy+Hsc88MAD\nampqUnJysmprazVnzhw1NDR0G9e5QAAAXd3+i/Xq1autvo9nDmEFAgE1NTVFHzc1NSk9Pb3LmBEj\nRkRPtM6aNUvt7e1qa2tLaE4AwE2eKZDJkyfr2LFjamxs1LVr11RVVaWSkpIuY1pbW6PnQOrr62WM\n0ZgxY1zEBYBBzzOHsJKSkvT6669rxowZ6ujoUGlpqXJycvTmm29KksrLy7Vjxw698cYbSkpKUnJy\nsrZt2+Y4NQAMXixl4iErVqzQhAkTtGLFCtdRnGEpEyDxWMoEAJBQFAgAwAoFAgCwQoEAAKx45lNY\ng11ZWZlqamo0cuRILVq0iOU8AHgeeyAe0dDQoJaWFn300UcsZQ6gX6BAPOLWFfZf+tKXWMocQL9A\ngXjE1q1blZWVpe9+97scvgLQL3AOpJPOK1R++OGHGj9+vEaOHJmQFX39fr9mzpwZ3RMZTDrPe3Jy\nstauXat77rmHlZQBj+NK9DgefvhhrV+/Xg8//HAvpLozXIkOwAWuRAcAJBQFAgCwQoEAAKxQIAAA\nKxQIAMAKBQIAsEKBAACsUCAAACsUCADACgUCALBCgQAArFAgAAArFAgAwArLuTvWeSnzf/3rXzpz\n5ozOnz/PUuYAPI/l3GMoKytTVVWVcnNzVVtbyw2eAAxoLOfeixoaGvTJJ5/o4MGD3J8cAOKgQGK4\ndVfAiRMncn9yAIiDQ1gxRCIRZWVlacuWLZo5c2YvJQMAb+IQVi/y+/366le/qhEjRriOAgCeRYEA\nAKxQIAAAKxQIAMAKBQIAsEKBAACsUCAAACueKZCdO3dq4sSJys7O1rp162KOWblypbKzsxUMBhUO\nhxOcEADQmScKpKOjQ88//7x27typo0eP6ne/+50++uijLmNqamp0/PhxHTt2TJWVlVq+fLmjtAAA\nySMFUl9fr6ysLGVmZmrYsGGaP3++3n///S5jqqurtXjxYklSYWGhIpGIWltbXcQFAMgjy7m3tLQo\nIyMj+jg9PV2HDh3qcUxzc7NSU1O7fb+Kioro1yyLDgBddb6NxN3wRIH4fL47Gnf7Wi3x/l7nAgEA\ndHX7L9arV6+2+j6eOIQVCATU1NQUfdzU1KT09PTPHNPc3KxAIJCwjACArjxRIJMnT9axY8fU2Nio\na9euqaqqSiUlJV3GlJSUaPPmzZKkgwcPyu/3xzx8BQBIDE8cwkpKStLrr7+uGTNmqKOjQ6WlpcrJ\nydGbb74pSSovL9fs2bNVU1OjrKwsDR8+XJs2bXKcGgAGN+4H0knnE0tvv/22iouLlZGRwYl4AAOa\n7XaTAonjww8/VGZmpkaNGtULqQDAuyiQ/9dbBQIAgwV3JAQAJBQFAgCwQoEAAKxQIAAAKxQIAMAK\nBQIAsEKBAACsUCAAACsUCADACgUCALBCgQAArFAgAAArFAgAwAoFAgCwQoEAAKxQIAAAKxQIAMAK\nBQIAsEKBAACsUCAAACsUCADACgUCALBCgQAArFAgAAArFAgAwAoFAgCwQoEAAKxQIAAAKxQIAMAK\nBQIAsEKBAACsUCAAACsUCADACgWSAHV1da4jdEOmO+PFTJI3c5Hpzngxky1PFEhbW5uKi4s1YcIE\nPf7444pEIjHHZWZmatKkSQqFQnrooYcSnNKeF98wZLozXswkeTMXme6MFzPZ8kSBrF27VsXFxWpo\naND06dO1du3amON8Pp/q6uoUDodVX1+f4JQAgM48USDV1dVavHixJGnx4sV677334o41xiQqFgDg\nM/iMB7bIo0eP1oULFyTdLIgxY8ZEH3d23333adSoURo6dKjKy8v13HPPdRvj8/n6PC8ADDQ2VZDU\nBzliKi4u1tmzZ7s9/9Of/rTLY5/PF7cE9u/fr7S0NP3nP/9RcXGxJk6cqKlTp3YZ44E+BIBBIWEF\nsmvXrrivpaam6uzZsxo3bpzOnDmjsWPHxhyXlpYmSbr33nv11FNPqb6+vluBAAASwxPnQEpKSvTO\nO+9Ikt555x3NmTOn25grV67o4sWLkqTLly/rj3/8o/Lz8xOaEwDwP544B9LW1qZvfetbOnXqlDIz\nM7V9+3b5/X6dPn1azz33nD744AOdOHFCTz/9tCTp+vXr+va3v61Vq1Y5Tg4Ag5gZQL785S+b/Px8\nU1BQYKZMmeIkw5IlS8zYsWPN/fffH33u/Pnz5rHHHjPZ2dmmuLjYXLhwwXmml19+2QQCAVNQUGAK\nCgpMbW1tQjOdOnXKFBUVmdzcXJOXl2c2btxojHE7V/EyuZyrq1evmoceesgEg0GTk5NjfvSjHxlj\n3L+n4uVy/b4yxpjr16+bgoIC88QTTxhj3M9VrEyu5ynWttJmngZUgWRmZprz5887zbB3717z97//\nvcvG+gc/+IFZt26dMcaYtWvXmh/+8IfOM1VUVJjXXnstoTk6O3PmjAmHw8YYYy5evGgmTJhgjh49\n6nSu4mVyPVeXL182xhjT3t5uCgsLzb59+5y/p+Llcj1Xxhjz2muvmWeeecY8+eSTxhj3//9iZXI9\nT7G2lTbz5IlzIL3JOD4iN3XqVI0ePbrLc5/nOpdEZZLcztW4ceNUUFAgSUpJSVFOTo5aWlqczlW8\nTJLbuUpOTpYkXbt2TR0dHRo9erTz91S8XJLbuWpublZNTY2WLVsWzeF6rmJlMjd/eU9ojtvd/vNt\n5mlAFYjP59Njjz2myZMn6ze/+Y3rOFGtra1KTU2VdPMTZ62trY4T3fTLX/5SwWBQpaWlcZePSYTG\nxkaFw2EVFhZ6Zq5uZfra174mye1c3bhxQwUFBUpNTdWjjz6qvLw8T8xTrFyS27n6/ve/r5/97Gca\nMuR/mzbXcxUrk8/nczpPsbaVNvM0oApk//79CofDqq2t1a9+9Svt27fPdaRuPus6l0Ravny5Tp48\nqSNHjigtLU0vvviikxyXLl3S3LlztXHjRo0YMaLLa67m6tKlS5o3b542btyolJQU53M1ZMgQHTly\nRM3Nzdq7d692797d5XVX83R7rrq6Oqdz9Yc//EFjx45VKBSK+9t9oucqXibX76metpV3Ok8DqkBi\nXSfiBbeuc5H0mde5JNLYsWOjb5Jly5Y5mav29nbNnTtXixYtin502/Vc3cq0cOHCaCYvzJUkjRo1\nSl//+tf1t7/9zfk8xcp1+PBhp3N14MABVVdXa/z48VqwYIH+8pe/aNGiRU7nKlamZ5991vl7Kta2\n0maeBkyBePk6kTu5ziXRzpw5E/363XffTfhcGWNUWlqq3NxcvfDCC9HnXc5VvEwu5+rcuXPRwxtX\nr17Vrl27FAqFnL+n4uXqvNpEoudqzZo1ampq0smTJ7Vt2zZNmzZNv/3tb53OVaxMmzdvdvqeiret\ntJqn3jqr79qJEydMMBg0wWDQ5OXlmTVr1jjJMX/+fJOWlmaGDRtm0tPTzdtvv23Onz9vpk+f7uxj\nhLdneuutt8yiRYtMfn6+mTRpkvnGN75hzp49m9BM+/btMz6fzwSDwS4fZXQ5V7Ey1dTUOJ2rf/zj\nHyYUCplgMGjy8/PN+vXrjTHG+XsqXi7X76tb6urqop94cj1Xt+zevTuaaeHChc7mKd620maePHEh\nIQCg/xkwh7AAAIlFgQAArFAgAAArFAgAwAoFAgCwQoEAAKxQIEACHT58WHv27NH69etdRwHuGgUC\nJNDhw4dVWFioc+fO6dKlS67jAHclYfdEByB95zvfUUdHh65fv66UlBTXcYC7wh4I0IcuXLigZ555\nRm1tbdHnqqqq9NJLL6m9vd1hMuDuUSBAHxo9erSmTZumHTt2SLq5SN2f//xnrVq1qsv9IYD+iLWw\ngD7W2tqqpUuX6oMPPnAdBehV/AoE9LHU1FRdvnxZn3zyiesoQK+iQIA+9umnnyolJYU9EAw4FAjQ\nhzo6OlRRUaFXXnlF7733nus4QK+iQIA+9OKLL2rRokUKhUI6deqUrl275joS0GsoEKCP7NixQw8+\n+KDy8vIkSU888YRqamocpwJ6D5/CAgBYYQ8EAGCFAgEAWKFAAABWKBAAgBUKBABghQIBAFihQAAA\nVigQAIAVCgQAYOX/AP/Sb/EzqXJbAAAAAElFTkSuQmCC\n" | |
} | |
], | |
"prompt_number": 14 | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"In order to get the value of the rotation measure, we need to minimize the following merit function:" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"\\begin{equation}\n", | |
"\\chi^2 = \\sum_{i=1}^N \\frac{\\left(\\beta_i - \\lambda_i^2 \\mathrm{RM} \\right)^2}{\\sigma^2}\n", | |
"\\end{equation}" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"which can be done analytically by computing the first derivative and equating it to zero. We obtain:" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"\\begin{equation}\n", | |
"\\mathrm{RM} = \\frac{\\sum_i \\lambda_i^2 \\beta_i}{\\sum_i \\lambda_i^4}\n", | |
"\\end{equation}" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"InferredRM = np.sum(wavelength**2 * ObservedBeta) / np.sum(wavelength**4)" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [], | |
"prompt_number": 15 | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The inferred value of the rotation measure is, then:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"InferredRM" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"output_type": "pyout", | |
"prompt_number": 16, | |
"text": [ | |
"0.030986361647827209" | |
] | |
} | |
], | |
"prompt_number": 16 | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"which is very close to the original one that we assumed. We can plot the resulting straight line:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"plt.errorbar(wavelength**2, ObservedBeta,yerr=0.2,marker='.',color='k',linestyle='none')\n", | |
"plt.plot(wavelength**2, InferredRM*wavelength**2)\n", | |
"plt.xlabel(r'$\\lambda^2$')\n", | |
"plt.ylabel(r'$\\beta$')" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"output_type": "pyout", | |
"prompt_number": 17, | |
"text": [ | |
"<matplotlib.text.Text at 0x9c23b4c>" | |
] | |
}, | |
{ | |
"output_type": "display_data", | |
"png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAESCAYAAADTx4MfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1UlvXhx/HPrWiJqOhJ0YClTkxEBXwIm/MXaWSWkVMr\nLZGphblTra11Vp1+C0/NWbNtNteDnVLUUsx8oEI2KxGmKZnYcz9IIwGVhsh8ToTr98c3CRQSL2/u\n6wLer3M8B7ivbj5+z9318Xr6fj2WZVkCAOACtXI6AACgaaJAAAC2UCAAAFsoEACALRQIAMAWCgQA\nYItrCqSwsFDXXnutIiIiNGDAAD377LN1bnf//fcrLCxMkZGRys3N9XFKAMAZfk4HOKNNmzb661//\nqqioKB09elRDhgxRXFycwsPDq7dJT0/XV199pfz8fG3fvl2zZ8/Wtm3bHEwNAC2Xa45Aunfvrqio\nKElSQECAwsPDtW/fvlrbpKWlKTExUZIUExOj8vJylZSU+DwrAMBFRyA1FRQUKDc3VzExMbV+Xlxc\nrNDQ0OrvQ0JCVFRUpKCgoOqfeTwen+UEgObCzqQkrjkCOePo0aOaNGmSFixYoICAgHNeP/svWVdh\nWJblqj+PP/644xnI1HwyuTUXmZpuJrtcVSAVFRWaOHGipk6dqvHjx5/zenBwsAoLC6u/LyoqUnBw\nsC8jAgC+55oCsSxLM2fOVP/+/fXAAw/UuU18fLyWLl0qSdq2bZsCAwNrnb4CAPiOa66BbNmyRcuX\nL9egQYMUHR0tSZo7d6727t0rSZo1a5ZuvPFGpaenq0+fPmrfvr0WL17sZOQGi42NdTrCOcjUMG7M\nJLkzF5kaxo2Z7PJYF3MCzIU8Hs9FndMDgJbG7n7TNaewAABNCwUCALCFAgEA2EKBAABsoUAAALZQ\nIAAAWygQAIAtFAgAwBYKBABgCwUCALCFAgEA2EKBAABsoUAAALZQIAAAWygQAIAtFAgAwBYKBABg\nCwUCALCFAgEA2EKBAABsoUAAALZQIAAAWygQAIAtFAgAwBYKBABgi5/TAQAAP8jMzFRmZqYkqaqq\nSqdPn1bbtm0VGxur2NhYR7OdzWNZluV0CG/yeDxqZn8lAC1UVlaWHnvsMWVlZTXq77G73+QUFgDA\nFgoEAGALBQIAsIUCAQDYQoEAAGyhQAAAtlAgAABbKBAAgC2uKZAZM2YoKChIAwcOrPP1zMxMderU\nSdHR0YqOjtaTTz7p44QAgJpcM5XJ9OnTdd9992natGn1bnPNNdcoLS3Nh6kAAPVxzRHIyJEj1blz\n5x/dhilKAMA9XHMEcj4ej0dbt25VZGSkgoODNX/+fPXv37/ObZOTk6u/duMEZADgpJoTNl4MV02m\nWFBQoJtvvlmffPLJOa8dOXJErVu3lr+/vzZs2KBf//rXysvLO2c7JlME0FwwmaKXdOjQQf7+/pKk\nsWPHqqKiQmVlZQ6nAoCWq8kUSElJSXVD5uTkyLIsdenSxeFUANByueYayJQpU7R582aVlpYqNDRU\nc+bMUUVFhSRp1qxZWr16tZ5//nn5+fnJ399fK1eudDgxADSepKQkffDBB/rmm29UXl6uwMBApyOd\nw1XXQLyBayAAmoPY2Fht3rxZUmvdeusErVq1qtF+l939pmuOQAAAxsmT0n//e72k36l160u1aNFQ\npyPVqclcAwGA5syypO3bpV/9SgoJkTp2fEj9+n2iIUP+5MrTVxJHIABcoL7nElrCc1xFRdLy5VJK\ninT6tJSYKH34oXTFFW2UlTVCjz22wemI9aJAADiuZlHs2rVL8+fP1/Lly50N1YiOH5fWrjWlsWOH\ndOut0ssvS1dfLXk8TqdrOAoEgKscP35ce/bscTqG11mW9O9/m9J44w1p+HBpxgxp/XqpXTun09lD\ngQBAI/r6a2npUvOnXTtziuqzz6TLL3c62cWjQADAy44ckVavNkcbn30mTZ4spaZKQ4Y0rVNU50OB\nAIAXVFVJmzZJS5ZIb74pxcZKv/61dNNNUtu2TqdrHBQIAFyEvDxzpLFsmXTZZeYU1V/+InXt6nSy\nxkeBAMAFOnTInJJKSTHXOKZOld5+W6pnQdVmiwIBgAY4fVr6179Mafzzn9L110uPPSaNGSP5tdA9\naQv9awNAw3zyiSmNV1+VevY0p6heeEE6zwKqLQIFAgBn+c9/pBUrTHF8+62UkCBlZkpXXul0Mneh\nQABA0qlTUnq6uYsqM1MaN0566inp2mul1q2dTudOFAiAFsuypJ07zZHGypVSeLg5RbV0qdSxo9Pp\n3I8CAdDi7N9vrmksWWLmpUpMlLZtk3r3djpZ00KBAGgRTp40806lpEjvvy9NmCA995z0859LrVjY\nwhYKBECzZVnmyGLJEjO1yJAh5mhj9WrJ39/pdHWrObV9aWmpWrdureTkZFdObU+BAGh29u41T4an\npJiji8REadcuKTTU6WTn58aiqA8FAsA1kpKStGPHDn399dcqLy+/oJX4jh2T1qwxRxu7dkm33WZK\n5KqrmtcEhm7CmT8ArpGXl6fc3FyVl5crKSnpvNtXVZlbbqdPN8vApqZK99wjFRdLzz8vxcRQHo2J\nIxAAruH//YWJ9u3ba9GiRfVut3v3D2tsdOgg/fKX0p/+JHXv7qOgkCR5LMuynA7hTR6PR83srwS0\nGOXl5Zo4caKOHDminJycWq/997/S66+b6xp5edKUKebaRlQURxkXy+5+k1NYAFwjMDBQTzzxhPy+\nn52wstJMYHjnndIVV5gnxR96SCoqkv72Nyk6mvJwEqewALjO8eNX6OGHpeXLpR49zJHGggVmvQ24\nBwUCNHE1nxuoqSndDipJZWVmOpF//GOg8vOf1fXXm2nTIyKcTob6cA0EaEY2bdqkU6dOacyYMU5H\naZCKCikjw1zXeOcdaexYadiwz/X667P0/vvZTsdrMbgGAkBbt25VVlaW0zHO66OPpN/+1tx6O2+e\nWZSpoMBMoT58eLk8nkqnI6IBOIUFwCe+/dZMYJiSYpaEnTZN+ve/pbAwp5PBLgoEQKP57jvprbdM\naWRlSbfcIv31r9I11zCBYXNAgQDwKsuSPvjAlEZqqjRokLmL6rXXpIAAp9PBmygQAF5RXGxuu01J\nMav7JSZKO3aYdcTRPFEgAGw7flxat86UxgcfSJMmSS+9JP3sZzzg1xJQIAAuiGVJW7aY0njjDTPb\n7fTppkjatXM6HXyJAgHQIAUFP0xgeMkl5hTVp59Kl1/udDI4hQIBUK+jR83qfSkppixuv908qzF0\nKKeo4KIHCWfMmKGgoCANHDiw3m3uv/9+hYWFKTIyUrm5uT5MB7QcVVXSe++Z5zRCQqS1a6X77jMT\nGC5cKA0bRnnAcE2BTJ8+XRkZGfW+np6erq+++kr5+flatGiRZs+e7cN0QPOXny899pjUq5f04IPS\n4MFm2vT166UJE8xpK6Am15zCGjlypAoKCup9PS0tTYmJiZKkmJgYlZeXq6SkREFBQT5KCDQ/5eXS\nqlVmGdg9e8y06W++aZ7d8KWaE0Lu27dPZWVlSk5ObnITQrY0rimQ8ykuLlZoaGj19yEhISoqKqqz\nQJKTk6u/5gMI1Hb6tLRxo7mukZEhxcVJjz5q5qNq08aZTPx/6lv1zeB8oZpMgUg6Z7ZITz0nYmsW\nCADj009Nabz6qvSTn5i7qJ57TurSxelk8LWzC3vOnDm23qfJFEhwcLAKCwurvy8qKlJwcLCDiQD3\nKy01d02lpEgHDkgJCeYCeb9+TidDc9BkCiQ+Pl4LFy7U5MmTtW3bNgUGBnL9A6ghKSlJ77zzjjye\nSxQR8Yhefz1AmzZJ48ZJf/qTNGqU1Lq10ynRnLhmQakpU6Zo8+bNKi0tVVBQkObMmaOKigpJ0qxZ\nsyRJ9957rzIyMtS+fXstXrxYgwcPPud9WFAKLZFlSUOH3q2dOwdImqLLLivVvHn9deutUseOTqeD\n29ndb7qmQLyFAkFLcuCAuaaxZIm0e/cBnTjxgrp2zVBeXoYCAwOdjocmghUJgRbi5Elz6+1NN0nh\n4dJnn5kH/IqKLtWAAW8oIWEE5QGfaDLXQICWzLKk7dvNkcbrr5uH/BITTZG0b39mq0BNnjxZx48f\ndzApWhIKBHCxwkJp2TJzF5VkSmPXLqnGI1GAYygQwGWOHZPWrDGlkZsr3Xab+Tomhjmo4C4UCOAC\nVVVSdrYpirVrzYJMSUlSfLx06aVOpwPqRoEADtq9+4c1NgICpF/+Upo7V+re3elkwPlRIICPHT5s\nLoSnpEhffilNmWJW9ouO5hQVmhYKBPCBykrp3XdNabz9tnTttWbK9LFjpbZtnU4H2EOBAI3oyy9N\naSxbZk5LJSZKCxZIl13mdDLg4lEggJeVlUkrV5riKCw0a2xkZEgDBjidDPAuCgTwgooK6Z//NKWx\ncaN0ww1ScrJZa8OP/8vQTPHRBi7CRx+Z0njtNal3b3MX1UsvScwkgpaAAgEu0LffmsJYssScrpo2\nTcrKkvr2dToZ4FsUCNAA330nvfWWOdrIyjIP+P3lL1JsrNSKKUnRQlEgQD0sS9qxw5TGypXSwIHm\nLqpXX5U6dHA6HeC88xbICy+8oO3btyskJERJSUlKTU1Vly5dNH78eHVhMWU0Q8XF0vLlpjhOnTKn\nqHbskHr2dDoZ4C7nXVBqw4YNGjt2rHbv3q1HHnlEM2fO1O7du5WRkaE//OEPGjp0qK+yNggLSsGO\nEyekdevMdY0PPpAmTjRHGyNGuP/p8MzMTGVmZkqSsrOzVVFRoVGjRik2NlaxsbGOZkPTYHe/ed4j\nkBMnTqiqqko//elPNWjQII0ZM0aSNHv2bC1YsMB1BYKmp+YOsKbG3gFalrRliznSeOMNadgwcxfV\n2rWSv3+j/VqvqzlOJSUlqqqqUo8ePZwNhRbhvEcgJ0+eVFpamnr16qVhw4bVem3NmjWaMGFCowa8\nUByBNG3R0dHKzs5WQEBAo/2OgoIfJjBs29YcaUydKgUHN9qvBFyt0ddE379/v3Jzc2VZliorK1VQ\nUKCIiAiNHj36gn9pY6JAmraOHTuqqKhIHTt29Or7Hj0qrV5tjjY++US6/XZztDF0qPtPUQGNrdEL\n5Gzl5eXKzc3ViRMn1KpVK91www123sbrKJCmzZsFUlUlbdpkSiMtTfqf/zFHG+PGSZdc4oWwQDPh\n8wJxKwqkafNGgeTn/zCBYefOpjTuuEMKCvJiUKAZabSL6EBTUF4urVpl7qLavdtMYJiWJkVGOp0M\naL4oEDRZp0+biQtTUsxst9ddJz36qDRmjNSmjdPpgOaPAkGT8+mnpjRefVUKDTWnqJ57TuK5VsC3\nKBA0CaWl0ooVpjgOHDC33b77rhQe7nQyoOWiQOBap05JGzaY6xqbNkk33STNnSuNHi21bu10OgAU\nCFzFsqRdu1rpjTfMEceVV5pTVEuWSJ06OZ0OQE0UCFzhwAFzTeP48fc1bZq/EhOl99+XfvpTp5MB\nqA8FAsecPGlutU1JkbZulcaPly655HfatStVgYHefRIdgPexFA58yrKkbduke+4xc08tWiRNniwV\nFUmLF0t+fltYoAloIjgCgU8UFponw1NSTIkkJkq5udJPfuJ0MgB2USAOc2oqc184dsxMjb5kiSmL\nW281Xw8fzgSGQHPAXFgu8uijj6p379666667nI5iW1WVlJ1tjjTWrpWuvtocbdxyi3Tppef/7xtr\nNl4A9bO73+Rss4scOXJEJ06ccDqGLXv2SMnJ5q6pe++V+veXPv9cSk83U6c3pDySkpJ0/PhxTZo0\nSeXl5Y2eGcDFoUBg2+HD0ssvm2nSY2KksjKzst/HH0u/+510oYvi5eXlqbKyUhs3blRSUlLjhAbg\nNa4qkIyMDPXr109hYWF66qmnznk9MzNTnTp1UnR0tKKjo/Xkk086kLJlq6w0Exjeeae5AP7WW9Jv\nfysVF0vPPisNHmz/+ob/9+vIDh48WIsWLfJiagCNwTUX0SsrK3XvvffqnXfeUXBwsIYNG6b4+HiF\nnzXZ0TXXXKO0tDSHUrZcX375wxobQUHmusbf/iZ17eq93/Haa6+pa9euWrt2rQIDA733xgAahWsK\nJCcnR3369FHPnj0lSZMnT9b69evPKZCmeoG8KSork1JTTXF8842ZwDAjQxowoHF+X2BgoNq1a0d5\nAE2EawqkuLhYoaGh1d+HhIRo+/bttbbxeDzaunWrIiMjFRwcrPnz56t///7nvFdycnL1183hdlhf\nOn3alERKijlVNWaM9Ic/SNdfL/m55tMC4GLU9/jAhXLNLsHTgBPngwcPVmFhofz9/bVhwwaNHz9e\neXl552xXs0DQMB9//MMaG717m1NUixaZJWEBNC9n/8N6zpw5tt7HNRfRg4ODVVhYWP19YWGhQkJC\nam3ToUOH6gutY8eOVUVFhcrKynyaszn59ltzHSM6Who3ztxqu3mzmZdq1izKA8CPc02BDB06VPn5\n+SooKNCpU6eUmpqq+Pj4WtuUlJRUXwPJycmRZVnqwjJ0F+S776Q1a8yDfX37Sjt3SvPnSwUF0h//\naKZPB4CGcM0pLD8/Py1cuFBjxoxRZWWlZs6cqfDwcL344ouSpFmzZmn16tV6/vnn5efnJ39/f61c\nudLh1E2DZUk7dphTVKmpUkSEOUW1fLnUoYPT6QA0VUxl4iL33Xef+vbtq/vuu88r77dvnymJlBQz\ndXpiopSQIPXq5ZW3bxRMZQL4nt39pmuOQOAdJ05I69aZ0sjJkSZMkF54Qfr5z5nAEIB3USDNgGWZ\nC98pKdLq1dKwYeZoY80a6ft7DgDA6yiQJuybb6SlS80fPz9TGh9/LJ118xoANAoKpIk5etRMWJiS\nYsri9tvNsxvDhnGKCoBvUSAukZSUpPT0dHXs2FEJCQm1pvOoqpIyM01prF8vjRwp/epX0s03S5dc\n4lxmAC2ba54Daeny8vJUXFysL774onoq8/x86X//19w19ZvfSFFR0v/9n/Tmm9KkSZQHAGdxBOIS\nZ56wDwkZoJ/9bIlGjJC++kq64w5z1BEV5XBAADgLBeIClZXSjBmrlJ29SaWlccrOvlS//700dqzU\npo3T6QCgbhRIDTVnqPzoo4/Uq1cvdezYsdFm9P3sM3NdY/lyKSQkQIMHn9DYscv08MN3e/13uVnN\ncff399e8efPUtm1bZlIGXI4n0esxYsQIPf300xoxYoQXUv3g4EFpxQpTHPv2mTU2EhPNGuLefhId\nABqCJ9FdrKJCSk83pfHee9KNN0pPPildd53UurXT6QDAHgqkkViWtGuXKY0VK6SwMHOksXix1KmT\n0+kA4OJRIF5WUmIe7FuyRDp8WJo2TdqyRerTx+lkAOBdFIgXnDxpns1ISTFlccst0oIF0jXXSK14\n0gZAM0WB2GRZZrbbM2tsREWZU1QrV0oBAU6nA4DGR4FcoKIiadkyUxyVlaY0du6UrrjC6WQA4FsU\nSAMcPy6tXWuua3z4oXTrrdIrr0hXX80EhgBaLgqkHuYuqo565RWzrsbw4dLMmVJamtSundPpAMB5\nFEg9vvhivp55prdmzzZPjF9+udOJAMBduEeoHmFhT2jp0l166CHKAwDqQoHUo02b/3J9AwB+BAUC\nALCFAgEA2MJFdIfVnMr8yy+/1P79+3Xw4EGmMgfgekznXoekpCSlpqaqf//+2rBhQ631yQGgubG7\n3+QUVh3y8vJ0+PBhbdu2rXp9cgBAbRRIHc6sT96vXz8tWrTI4TQA4E6cwqpDeXm5+vTpo+XLl+uG\nG27wUjIAcCdOYXlRYGCgrrzySnXo0MHpKADgWhQIAMAWCgQAYAsFAgCwhQIBANhCgQAAbKFAAAC2\nuKZAMjIy1K9fP4WFhempp56qc5v7779fYWFhioyMVG5uro8TAgBqckWBVFZW6t5771VGRoY+//xz\nrVixQl988UWtbdLT0/XVV18pPz9fixYt0uzZsx1KCwCQXFIgOTk56tOnj3r27Kk2bdpo8uTJWr9+\nfa1t0tLSlJiYKEmKiYlReXm5SkpKnIgLAJBLpnMvLi5WaGho9fchISHavn37ebcpKipSUFDQOe+X\nnJxc/TXTogNAbTWXkbgYrigQTwPXjj17rpb6/ruaBQIAqO3sf1jPmTPH1vu44hRWcHCwCgsLq78v\nLCxUSEjIj25TVFSk4OBgn2UEANTmigIZOnSo8vPzVVBQoFOnTik1NVXx8fG1tomPj9fSpUslSdu2\nbVNgYGCdp68AAL7hilNYfn5+WrhwocaMGaPKykrNnDlT4eHhevHFFyVJs2bN0o033qj09HT16dNH\n7du31+LFix1ODQAtG+uB1FDzwtIrr7yiuLg4hYaGciEeQLNmd79JgdTjo48+Us+ePdWpUycvpAIA\n96JAvuetAgGAloIVCQEAPkWBAABsoUAAALZQIAAAWygQAIAtFAgAwBYKBABgCwUCALCFAgEA2EKB\nAABsoUAAALZQIAAAWygQAIAtFAgAwBYKBABgCwUCALCFAgEA2EKBAABsoUAAALZQIAAAWygQAIAt\nFAgAwBYKBABgCwUCALCFAgEA2EKBAABsoUAAALZQIAAAWygQAIAtFAgAwBYKBABgCwUCALCFAgEA\n2EKB+EBmZqbTEc5BpoZxYybJnbnI1DBuzGSXKwqkrKxMcXFx6tu3r66//nqVl5fXuV3Pnj01aNAg\nRUdH66qrrvJxSvvc+IEhU8O4MZPkzlxkahg3ZrLLFQUyb948xcXFKS8vT6NHj9a8efPq3M7j8Sgz\nM1O5ubnKycnxcUoAQE2uKJC0tDQlJiZKkhITE7Vu3bp6t7Usy1exAAA/wmO5YI/cuXNnHTp0SJIp\niC5dulR/X1Pv3r3VqVMntW7dWrNmzdLdd999zjYej6fR8wJAc2OnCvwaIUed4uLidODAgXN+/sc/\n/rHW9x6Pp94S2LJli3r06KH//Oc/iouLU79+/TRy5Mha27igDwGgRfBZgWzcuLHe14KCgnTgwAF1\n795d+/fvV7du3ercrkePHpKkrl276he/+IVycnLOKRAAgG+44hpIfHy8UlJSJEkpKSkaP378Odsc\nP35cR44ckSQdO3ZM//rXvzRw4ECf5gQA/MAV10DKysp02223ae/everZs6dWrVqlwMBA7du3T3ff\nfbfefvtt7dmzRxMmTJAknT59WnfeeaceeeQRh5MDQAtmNSNXXHGFNXDgQCsqKsoaNmyYIxmmT59u\ndevWzRowYED1zw4ePGhdd911VlhYmBUXF2cdOnTI8UyPP/64FRwcbEVFRVlRUVHWhg0bfJpp7969\nVmxsrNW/f38rIiLCWrBggWVZzo5VfZmcHKsTJ05YV111lRUZGWmFh4dbDz/8sGVZzn+m6svl9OfK\nsizr9OnTVlRUlDVu3DjLspwfq7oyOT1Ode0r7YxTsyqQnj17WgcPHnQ0Q1ZWlrVz585aO+uHHnrI\neuqppyzLsqx58+ZZv//97x3PlJycbD3zzDM+zVHT/v37rdzcXMuyLOvIkSNW3759rc8//9zRsaov\nk9NjdezYMcuyLKuiosKKiYmxsrOzHf9M1ZfL6bGyLMt65plnrDvuuMO6+eabLcty/v+/ujI5PU51\n7SvtjJMrroF4k+XwGbmRI0eqc+fOtX52Ic+5+CqT5OxYde/eXVFRUZKkgIAAhYeHq7i42NGxqi+T\n5OxY+fv7S5JOnTqlyspKde7c2fHPVH25JGfHqqioSOnp6brrrruqczg9VnVlssw/3n2a42xn/347\n49SsCsTj8ei6667T0KFD9dJLLzkdp1pJSYmCgoIkmTvOSkpKHE5k/P3vf1dkZKRmzpxZ7/QxvlBQ\nUKDc3FzFxMS4ZqzOZBo+fLgkZ8eqqqpKUVFRCgoK0rXXXquIiAhXjFNduSRnx+o3v/mN/vznP6tV\nqx92bU6PVV2ZPB6Po+NU177Szjg1qwLZsmWLcnNztWHDBv3jH/9Qdna205HO8WPPufjS7Nmz9fXX\nX2vXrl3q0aOHHnzwQUdyHD16VBMnTtSCBQvUoUOHWq85NVZHjx7VpEmTtGDBAgUEBDg+Vq1atdKu\nXbtUVFSkrKwsbdq0qdbrTo3T2bkyMzMdHau33npL3bp1U3R0dL3/uvf1WNWXyenP1Pn2lQ0dp2ZV\nIHU9J+IGZ55zkfSjz7n4Urdu3ao/JHfddZcjY1VRUaGJEycqISGh+tZtp8fqTKapU6dWZ3LDWElS\np06ddNNNN+nDDz90fJzqyrVjxw5Hx2rr1q1KS0tTr169NGXKFL333ntKSEhwdKzqyjRt2jTHP1N1\n7SvtjFOzKRA3PyfSkOdcfG3//v3VX69du9bnY2VZlmbOnKn+/fvrgQceqP65k2NVXyYnx6q0tLT6\n9MaJEye0ceNGRUdHO/6Zqi9XzdkmfD1Wc+fOVWFhob7++mutXLlSo0aN0rJlyxwdq7oyLV261NHP\nVH37Slvj5K2r+k7bs2ePFRkZaUVGRloRERHW3LlzHckxefJkq0ePHlabNm2skJAQ65VXXrEOHjxo\njR492rHbCM/O9PLLL1sJCQnWwIEDrUGDBlm33HKLdeDAAZ9mys7OtjwejxUZGVnrVkYnx6quTOnp\n6Y6O1ccff2xFR0dbkZGR1sCBA62nn37asizL8c9Ufbmc/lydkZmZWX3Hk9NjdcamTZuqM02dOtWx\ncapvX2lnnFzxICEAoOlpNqewAAC+RYEAAGyhQAAAtlAgAABbKBAAgC0UCADAFgoE8KEdO3Zo8+bN\nevrpp52OAlw0CgTwoR07digmJkalpaU6evSo03GAi+KzNdEBSPfcc48qKyt1+vRpBQQEOB0HuCgc\ngQCN6NChQ7rjjjtUVlZW/bPU1FQ9+uijqqiocDAZcPEoEKARde7cWaNGjdLq1aslmUnq3n33XT3y\nyCO11ocAmiLmwgIaWUlJiWbMmKG3337b6SiAV/FPIKCRBQUF6dixYzp8+LDTUQCvokCARnby5EkF\nBARwBIJmhwIBGlFlZaWSk5P1xBNPaN26dU7HAbyKAgEa0YMPPqiEhARFR0dr7969OnXqlNORAK+h\nQIBGsnq228drAAAATklEQVT1ag0ZMkQRERGSpHHjxik9Pd3hVID3cBcWAMAWjkAAALZQIAAAWygQ\nAIAtFAgAwBYKBABgCwUCALCFAgEA2EKBAABsoUAAALb8P7OS1r/JRF8mAAAAAElFTkSuQmCC\n" | |
} | |
], | |
"prompt_number": 17 | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Note that the same cannot be done using standard linear regression routines, because the linear fit has to go through the (0,0) point. The result of such a fit gives a different slope to what we used at the beginning:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"(a,b) = polyfit(wavelength**2, ObservedBeta, 1)" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [], | |
"prompt_number": 18 | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [ | |
"a,b" | |
], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"output_type": "pyout", | |
"prompt_number": 19, | |
"text": [ | |
"(0.03312031780914055, -0.070886302133340476)" | |
] | |
} | |
], | |
"prompt_number": 19 | |
}, | |
{ | |
"cell_type": "code", | |
"collapsed": false, | |
"input": [], | |
"language": "python", | |
"metadata": {}, | |
"outputs": [], | |
"prompt_number": 19 | |
} | |
], | |
"metadata": {} | |
} | |
] | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment