Created
November 11, 2015 09:28
-
-
Save aasensio/af42a869f11ffc947cdd to your computer and use it in GitHub Desktop.
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
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as pl\n", | |
| "import seaborn as sn\n", | |
| "import scipy.special as sp\n", | |
| "import scipy.stats as st\n", | |
| "from ipdb import set_trace as stop\n", | |
| "import daft\n", | |
| "%matplotlib inline" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Some chromospheric lines display filamentary structures when observed in a filter very close to the core. Examples of these lines are the H$\\alpha$ or some of the lines of the Ca II infrared triplet. These filamentary structures are thought to be proxies of magnetic field lines. Recently, de la Cruz Rodríguez & Socas-Navarro (2011) used observations of linear polarization in the Ca II 8542 A line to infer, using a few points, that the linear polarization line seems to be aligned with the fibrils but some misalignment was found. The reason for this scarcity of results is that the observed signals are very weak and the inference of the field is not reliable at all.\n", | |
| "\n", | |
| "We plan in this work to make a statistical re-analysis of that work and tried to infer whether the alignment is supported by the data or not. The idea is then to compute the angle of the fibril from images on the core of the Ca II line and compare these angles with those inferred from the spectropolarimetric data. We do this study using a Bayesian approach to fully take into account the possible uncertainties in the inferred azimuth angles." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Obtaining azimuth angles from fibril images" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "We obtain the azimuth angle of the fibrils using the rolling Hough transform (Clark et al. 2014). This tranform is a generalization of the Hough transform and is obtained after the following steps:\n", | |
| "\n", | |
| "- A smoother version of the image is computed and then substracted from the original image. The smoothing is obtained using a top-hat kernel of a certain width $D_k$.\n", | |
| "- The resulting image is then thresholded and binarized to create a bitmask.\n", | |
| "- A disk of a certain diameter $D_w$ is extracted at each point in the image and the standard Hough transform is computed in each disk.\n", | |
| "- Finally, the Hough transform is again thresholded at a certain level $Z$ to make sure that only obvious fibrils are detected as so.\n", | |
| "\n", | |
| "The output of the rolling Hough transform is the array $R(\\theta,x,y)$ at each pixel position $(x,y)$ and angle $\\theta$. The dominant angle at each pixel position is obtained by computing the circular statistics average at each pixel:\n", | |
| "\n", | |
| "$$\n", | |
| "\\theta_\\mathrm{RHT} = \\frac{1}{2} \\arctan \\frac{\\int d\\theta \\sin 2\\theta R(\\theta,x,y)}{\\int d\\theta \\cos 2\\theta R(\\theta,x,y)}\n", | |
| "$$\n", | |
| "\n", | |
| "The following image shows the inferred azimuth angle of the fibril from the line core images.\n", | |
| "\n", | |
| "<img src=\"figs/fibrils.png\">" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Azimuth angle from polarimetry" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The inference is done using the Zeeman effect. In the weak field regime, it is possible to demonstrate that the linear polarization profile are given by the following expressions (where all the quantities are constant along the line-of-sight):\n", | |
| "\n", | |
| "$$\n", | |
| "Q(\\lambda) = \\beta B_\\perp^2 \\frac{\\partial I(\\lambda)}{\\partial \\lambda} \\frac{1}{\\lambda_0-\\lambda} \\cos 2 \\phi \\\\\n", | |
| "U(\\lambda) = \\beta B_\\perp^2 \\frac{\\partial I(\\lambda)}{\\partial \\lambda} \\frac{1}{\\lambda_0-\\lambda} \\sin 2 \\phi, \\\\\n", | |
| "$$\n", | |
| "\n", | |
| "that is particularized only for the wavelengths on the wings of the line. In the previous equations, \n", | |
| "$\\beta=4.09 \\times 10^{-26} \\bar{G} \\lambda^4$ is a constant that depends on the specific spectral line of interest. The magnetic field is the component perpendicular to the LOS. Finally, $\\phi$ is the field azimuth in the plane of the sky, $\\lambda$ is the wavelength and $I(\\lambda)$ is the intensity profile of the spectral line." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Bayesian modeling" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.axes._axes.Axes at 0x10950a390>" | |
| ] | |
| }, | |
| "execution_count": 2, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAALgAAAENCAYAAAC8Q6+1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmcFtWV97+nGojNIpuGGJYAal81hgANisQdSDIahxjE\nmHFJTITBLJPkjQiibxwTNaJJHLMZhTGiGA1oYnRifFVAo1FWNR01VDuCETEuNGBYBZ467x+3Gpq2\nl2epepbb9/v59Kf7qa4691TVr+5z6i7niqri8bhKUGoHPJ408QL3OI0XuMdpvMA9TuMF7nEaL3CP\n03iBe5zGC9zjNF7gHqfxAvc4jRe4x2m8wD1O4wXucRovcI/TeIF7nMYL3OM0XuAep/EC9ziNF7jH\naTqV2gFPsohIX2AUUAscCnwA2AWsB1YCq1T19dJ5WFzETzqufEQkAD4FfBU4AXgOWAWEwE6gCzAY\nK/rRQD3wC+BeVX2vBC4XDS/wCkdEhgO/AiLg58A9qrq9jf07AZ/BPgxHAFNV9eFi+FoKvMArlLjW\nvhz4BjAduENzvJkiMh6YCywCvq6qOxJ3tMR4gVcgsbjnAEcCk1V1fQG2esS2+gFnqOrWZLwsD7zA\nKxAR+SkwDDhNVbclYK8KK/IBwOmqurtQm+WCbyasMETks8CnsbVtweIGUNUMMCX+OD0Jm+WCr8Er\niLgJ8K/A2ar6VAr2BwLPAqeo6gtJ2y8FXuAVhIj8EKhW1a+lWMbF2G+H09Iqo5h4gVcIIlINvAYc\nq6prUiznAGAdMEZVX0mrnGLhY/DKYTKwPE1xA6jqTmy7+tQ0yykWvqu+chgH3N/WDsaYXtgHYSUw\nPgzDG4wxPcMwfDfHsu4HfpKfm+WFr8Erh1qscNtiQRiGc8IwfA5YY4yZgh2XshdjzLgsynoeOEpE\nPpCfq+WDF3gFEAvtMODF1vYxxpwFPNJk07vABGBFruXFXf1rsV35FY0XeGXQDdiuqrva2GcIVpSN\nKLAxDMN/gn0A4hodY8ykLMrcCByYp79lg4/BKwPJYp9bganGmBFAH6xAiT8fihX8K2EYLjbG/BK4\nrx17igMVoBd4ZbAN6CoinVvrRo9fJG9otnkagDHm82EYzjTGjItfRLOhN/DPvD0uE3w7eIUgIn8D\nzlHVv+R6rDFmCDAUGAlsDsNwTjtlVQMNQJ+42bBi8TV45bAK25KSs8DDMFwLrDXGbIxbWNpjGBBW\nurjBgRirA/E48K/5HmyMGZGluAHOiMureHyIUiGISHfg78BwVV2XYjldsEMCTlbV1WmVUyx8DV4h\nxBMR7gK+lXJR5wIvuiBu8DV4RSEiH8LG4Keranu9moXYP01VVyVtvxT4GryCUNU3gW8Dt4tI1yRt\nx9PgbgHmuCJu8DV4xSEiAtwGfAg4M4mWjtjmTdhmxHEupZLwNXiFEc+cn4LthPkfETmoEHvx+O85\nwLHAZ1wSN3iBVySqugf4N2yCnzoROTMfOyIyJrZxIDBeVTcn52V54EOUCkdEPoGdoLAem/jn923N\nio/DkVOwiX+OB76hqguL4Wsp8AJ3gLjt+kysaI8ElrIvddsObOq2Idie0GOBLdjUbfNVdUspfC4W\nXuCOISIfweYfbEy+ORl4GCv2VdhJE6tzzYJVqXiBO46IKPAxV9JA5Ip/yfQ4jRe4x2m8wD1O4wXu\ncRovcI/TeIF7nMYL3OM0ZdcOHrfbeiocVc0m1UXqlOWk43K5OC5Qio6ecqqkfIjicRovcI/TeIF7\nnMYL3OM0XuAep/EC9ziNF7jHabzAPU7jBe5xGi9wRxGRzk2yX3V1YUGpfPACd5czga3x38uAWSX0\npWR4gbvLKmzKCLBZsJzJN5gLXuDusoZ9i1d1wQvc4xJx3pOX4o/vAW+U0J2S4QXuNn+Kf9d1lEQ/\nzfECd5tl8e8nSupFCfECd5vGVSByXs7bFcpyyloaM3pEpAoYD5yIzdt3NNCDeMlrbBrhVcAjaSwP\nUgxEpCfwWeAY9uUm7IVNtrkOK/iVwP2q+o8U/UjlHuaD8wIXkQOBi7Gr/m4AHsIK+S/Au9iWhg9i\nVzcYBXwOeAebini+qmaS8iUtRMRglzb5PLAEeJJ92WV3Ap3Zl132OOxyhI8C/6WqT6fgT9kIHFUt\nqx/iBoCEbE3ALr13FzA6y2OqgNOxL2jPAEeU+pq04Wsn4DLsA/ld4JAsjzsQ+Bq2Vv8F0L1c72HB\nvpTagTQuDvbd4kexuD9ZgI2vxuL5UqmvSwv+fRCbB/xR4CN52uiFTZ6/BjiqnO5hYr6U2oGkL04s\nzF/FNXCvBPw5In5Qvlrqa9PEp37AauD7xGFmgfbOB94EhpXDPUzypyzTRhTI94EabM29rVBjqrpa\nRE4GnhSRN1T1/kJtFkI8aOoh4DeqemUSNlX1ThHZBfxRRGrVLlfoBE69ZIrIccBvsTXROwn7NRa4\nLw3bOfpxDfBR7BKCid682PZRwOcKsV1OL5nOCDxeDu954ApVvTd5z0BEfggMUNVz0rCfRfm12Nr7\n42nUsvG3w0rgWlW9uwA7XuCtUYDALwDOV9UJKbjVWEY1sBY4WUuwlruI3AcsUdWfpVjGycDN2JfO\nvMThBd4GBQj8GeAHqvpAS/83xgwBbgU2Ab+JN08AFoZhuCiHcq4Buqnqt3L1sRBEpD/wV2yLSWor\no8XLDP4Vu7zgkjxteIG3Rj4XR0SOBB4BBmsbHTPGmAXApWEYvhp/Hgo8EobhYTmUNQgbCh3cVllJ\nIyIzseK+uAhlfR0Yo6rn5Xl82QjclbEoxwFPZCG4kY3ijpkB5LQIqqq+BjRgmw+LyXHA4nwONMaM\nyPGQxcCYfMoqN1xpJqylnQH9xpiRwJr4ZvfBhiePhGF4Xx7lrYrLfDGPY/OlFtsdnw99ctw/BPqJ\nSC+t8OW9XanBDe2LbTzwyzAMnwvDcFEYhjOxYzf2N2Rj9fZ4iSLW4CLSHSvStXma6JXLzvE3YYi9\nrhWNKwLvCrTXqTMeeKzZtpHGmMHNt2VR3lagOjvXEuEAYEcurRrGmF7GmOnGmElAH2PMFGNMzxzK\n3EZxzzEVXBH4HtoPt4aGYfjPxg/GmPHApmYxebZ0BnbncVy+7MEOAsuFW8MwvAFYE4bhHGABMDuH\n4ztR3HNMBVdi8LeA/i39I661ZgEa12Zgv+5rgXF5ltcfeC3PY/NhC9BZRHpk00RojDkLWB5/bIy/\nD8U2kWZLf+DtnLwsQ1wReONL36+b/yMMw3exrSUzmv1rTgHl1QJFG5OiqhkRqQNGsG+eZZuH8P54\nfSZwaTbliUhfoDfwSi5+liOutINPAL6rqifkU2b8YtkYe0/ADkGlpRYWEemCnQE0oJgtDCLyM+BV\nVf1hNvsbY6Zjh8E21uCPZhuOici/ADNU9eQ8XC2rdvCSD2ds/kMeQy2xL0PvAEOL4N/Z2O7yYl+X\nCdhZSFkPj62pqRmXZ1l3Ad8q5j1M68eJl0xV3QHMA/69CMV9DTsLptgswj7Ix6VZiIh8EDgNez0r\nHidClPi4Q7EzXD6uqqkkuRGRU4HbgUNVtegtDCLyDex8yk9qSjdORG4EqlV1WgE28rqHaeBEDQ6g\nqq9ga9Zb4wFDiRJ3tvw3MK0U4o65Gfvy95U0jIvI8djOr8vTsF8KnKnB42O7YHOA3KKqiYUR8QNz\nOxCp6oVJ2c3Tl6OxM+dPUtWX2ts/B7sHYydZX6IFzloqpxrclWZCAFR1l4hMAp4QkS2qemehNkUk\nAOZi49JsuvFTRVVfEJHrgWUiMloTGJcuIgcBD2OnwZV0Sl7SOBOiNKKq/4ttcbhWRK6Ka/W8iNuD\n7wGOBQ7CDsktKfHQ4OuB7sASEflMgfZqgaewTaNXFO5heeGcwAHir+4x2EQ+y0VkVC7Hi0gQfxP8\nFZuVdTQ2V8pxIpJ4opwc/DqSfRljA+ALwE9EZF4cYuRiq5uIXA38EbhaVWem9eJaUkrdTplmGyo2\na9WXgNex2Z6+APRuY9/+wLewI+meB05ots9p2F7Cp0twXY6My1aatIVja/KbsN3wdwAnAB9oxUYV\n8PF4/wbszKaskgWV6h4W+uPUS2YbNjsDZ2Dbycdix1j8BcgAETZH4UisAB7BtsY8rS1cHBE5DfgD\n8Iyqjk3Szzb836/mbsWvPtiH+QJs2oy/AfXYkZZbgcFYcf8DO/DqFrWTN9LwN/F7mC8dQuDN7Fdh\nBfAxbA22FjuR4DlgXUviacFG0USejbhbOKYrVsyHAndiB5stA57VIgwv8AJvg2JeHBFRYIWqHpPH\nsamLPB9xt2BDgY+p6guJOtdOmeUicCdfMouBqj5Eii+eSYjb4wVeEGmJ3Is7ObzACyRpkXtxJ4sX\neAIkJXIv7uTxAk+IQkXuxZ0OXuAJkq/IvbjTwws8YXIVuRd3uniBp0C2IvfiTh8v8JRoT+Re3MXB\nCzxFWhO5F3fx8AJPmeYi9+IuLl7gRaCpyPHiLipe4MWjaaapZ7y4i4MXeBFoFpaUfGZQR8ILPGVa\niLlTHYXo2R8v8BRp7YXSi7x4eIGnRHutJV7kxcELPAWybQr0Ik+fDilwEakWkcbVikeLyMQEbefU\nzp3ipInDmpzjmXH+kw5Hh5yTGU883g40JgV6VlULFkAhnThJz/EUkU8BD2KXW8kAl6rqjwu1m2XZ\nfk5mKVG7itjLjR+BPxdqs9AeyhRq8lXYcwObNmJlAjYrjg4p8Jgn499bsWmX8yap7vckRa6qG4DG\nRbe6YtNidDg6ssCfwYo7oJ1FZNsi6bElCdfkjaJ+R1Nc376c6cgCX4WNT6vYF65khYgMiF9UUxk4\n1ZLIReSwPEw9gQ1TViThVyXSkQW+Gnv+q1U1yvagOJ3y/diYNrWBU81EvgN4SUROytHMSmzOxSeS\n9K2S6JACj1ObDccm5XxFRI6O8xdmwzewy3gfFX/umtbAqVjkN2FXOu4M3CMi3bI5VkT6se/+vici\ng9NY+aLc6TDNhCIyHJgKnAgMxSamBLuabw9gIPACNvnmraq6rgUbhwJ12Jc2sIk7p6lqIWtutuVz\nd2wSzX5Yge8A7lTV9y22JSKdgM9gk28eE/tYjz23LcAgbLPoSmAhcLeqbk/J77JpJix5etvmPySc\nehc4CXgauzLxFdgssl1a2K8HVvw/waYW/i12sanG/wdYcWSwQnsHuxhU2tejFzYJ/zZsPL0du3xJ\n4/+rgP8A1mET2X8R+wC/b7lB4BDgc8ADwAbgB9hvoLK+hwX5UmoH0ro4QDfgp9gw5GygUw7Hdgcu\niUX8jVjc38TW2NuxebgPLPJ1OS1+8DLYFMjdsFly/ww8DgzP0d5g7HqY9cDx5XgPE/Gl1A6kcXGA\nD2FDiTtpJeF9lnYOj2v/PwLvFavWbsOfxtpcgf+J/fk69iU3X5ufjR+YqeV0DxPzpdQOJH1xgIOx\ncet3W/qazsPeB4DHsDNy8n5YEr5GM+NvkhMSsncY8GpSIi8ngTv1khk34T0GLFfVmQn61AU7ruO5\nJO3m6ctAbBv+ZFVNrPkvbmd/Mrb7VIG28r6HSeOawL8GnIeNKTMJ+9UPu+zJRFVdlqTtHHwQbLj0\npKpek4L9zwI3YFeLzruFxQu8DfK9OCJyCHZVtOM1gbUjWynj88D/BYZpDp1DCZZ/HnaRrDGquiel\nMu4C/q6qswqw4QXeGgUI/LvYFcMubrrdGDOuycdeQJ8wDPNqt45r0GeBGapa1DUz47JXAbNU9WEA\nY8x4bLv9VGBhGIbvxtuuA64Nw/C3eZQzGNscOijfWrycBF7yl4AkXlCwnSCvY2vWvdtrampuqamp\nObXZtutqamrGFeDfRcDvS3BdjgVeoUmLSU1NzdCampqVzfetqan5XIFlPQh8uZj3MK0fV7rqxwBv\nqWpd4wZjzKXAy2EYLm627wpgcgFl3Q2Mz7bLPEEmA/N0/9BoPHaF4r3E31iLCizrV9i+g4rHFYGP\nwg5/BcAY0wuYGYbhD1vYt28hBanqNmwz5PBC7OTBfucYM4FmAgd6hWH4boFlLQVGuTB2xSWBr2r2\nubUhohMofM35VUDR5jjGzZ8jsfF/U8a18A1VMKr6BnaMzqCkbRcbVwQ+gP1To/Vs9hkAY8xQYHA+\nL1/NWItd9rtY9ABQ1YbGDfG5rGm6kzFmSPNtBVDsc0yFTqV2ICE6Y2ucRh4DLgMwxozEdm2vBS4F\nxsXbx2HHSo8ANgO1YRhOM8YswE7+PR24BRgZhuENzcrbxb4Jy8WgS1xmUza2sN/IMAzvg4LPD+z1\nLOY5poIrNfh72DHTAMQx6BRjzHRsKDEemByG4bQwDBvnKY6Mtz8WNxv2ibf/ADuhYWMYhouwy2E3\npzous1i8BxzQNCYOw3AzcIsxZroxZpIxZlKjuGMKOT+w17OY55gKrtTgrwCGJq0HYRg+Rzwn0Rgz\nKf49HtD4xt4b7zob22KwwhgzJQzDOfFLaltrutdgh6YWiy3YIbofBtY3bmynPT/v84sfpBrsda1o\nXKnB23vpexaYA0yKxQ1wFnbcdOML573sq81GASPimHaUMWZwM3u1FDBROVfUNi7n+mJbyPkNBbaq\n6tt5O10mONGTGWdtuhswmucJGWN6AkPjmr+tsg7CxvN9VbV5XJwaIvI9oFpVp+dzfLbnF5f1JeAM\nVZ2UT1nl1JPpSg3+LPal6MR8Do5j9fHZ3HzgK8DCYoo7Zj5wgYh8INcDczw/gCnYSR0VjxM1eHzc\n14ETVTW1Hrg45dv/YoeUFj1TlIg8CtyuqnelWMZw7JS2oZrngC5fg6fDHcAn8kitkAv/B1hTCnHH\nzAauEZED0zAedyjdBPwwX3GXG87U4PGxZwD/hR10tS1hv47ETggYrarv60QqFiIyF8hoCzPrE7D9\nH9gWl5O0gPH05VSDOyXw+PjbsZOGz0mqFhKRvtjkOT9X1ZuTsFmALz2B54HvqeqvErR7InAf8AlV\nrW9v/3ZslY3AXQpRGvl34EBgfg7JfFpFRD6IHdD0EPDLQu0Viqq+C3waG6pcmIRNETkZ24x4TqHi\nLjecE7iqvgdMxE4W/rOIHNXOIa0iIv+CHfz/IHaSQ1l83alqCJwCfFdEbsk3JheRziJyBbAA+Lyq\nFjrMtuxwTuAAqroDm+Dmv4EnROR7IvKhbI8XkaNEZB7wC+BCVb2yXMTdSCzy4dh7WCci52XbhCgi\nQfzwLgWOB2pVdUl63pYO52LwFux9BJiFfXn6f9hJuyuxSTcz8T4HAMOwPXyTsbkH5wA3aAWkHRaR\nccAM4OPYXDBPYns+1zc+mPF7RC12csgXsV31PwZ+nfTDW04xuPMCb2K3J3AucAL2Rg/CjjKMiLPM\nYkXxCHB/CTpyCkZEarDneAz2HLthRwTuxJ7nc9iHewGwIq1vJS/wNijWxYlr7R3YXtAxqrq7nUMq\ninjAVHfsKg9jsIIuSiaAchK4kzF4NqjqzvjPjGviBjtAq0l4ta1Y4i43OqzAPR0DL3CP03iBe5zG\nC9zjNF7gHqfxAvc4jRe4x2m8wD1O4wXucRovcI/TeIF7nMYL3OM0XuAep/EC9ziNF7jHabzAPU7j\nBe5xGi9wj9N0SIGLSHcRaVwNYbSI/KykDqWAiJwkIg/EH+fFqyR3ODqkwIFt2OU9GvloqRxJkZ3A\nJ+O/j8aB9XbyoUMKPE6X0Lho7B7g8dJ5kxp1QFX893vYdBEdjg4p8JgngAy2Nm9tTc2KJc7u9Xr8\n8QDgpRK6UzI6ssBXYMV9AEVcb6fILIt/v+JKvu9c6cgCX4lNjLNdVd8qtTMp8af4959L6kUJ6cgC\nfwPYis217SqN30wdVuAdKnVbnM7sZOBfgVoRGQFsUdUXsSK4T1X/mkbZxUJEugGfxeYQP0JERgEv\nxtloH8We4zsp+1A2qds6hMBFZISI3KiqJ4oIPXr0yAwYMKBTt27dUFU2b94cvfHGG9H27ds7icgu\nVb0XuERV/5GkH2kiIl8MguDKKIqGVFVVRX379o0+/OEPd+rSpQt79uzhrbfeyrz99tvs3r27KgiC\nDVEU3QjMLmSpkjZ88QJvjSQvjoj0DoJgURRFIwYMGJA5//zzq8aOHdvq/lEUcf/99/O73/0us3nz\n5ioR+bWqnlduucGbIiInBUHwANBj+PDhTJkyRQYMGNDq/ps3b+a2227jySefjDKZTEZVp6nqbQn7\n5AXeGkldHBH5kojM6d27t1x99dVVAwcOzOn4ZcuWcf3110e7d+/eoqrjVLWsWlrEcreqfn7UqFHR\nZZddFnTpkn1fThRFzJ07lwcffJAgCFZGUXR8vDpGEr55gbdGEhdHRGYB10ycOJGLLroobzu7du3i\nyiuvjF544QWACaq6uBC/kkJEJAiCFSIy8vLLL5fRo0fnbevVV19l5syZmR07djREUXRYEgn/vcDb\nIIFV1r4G/Oziiy/mtNNOS8SnH/3oR/r4448rcGwJ18jcS1VV1VNVVVXH/fznPw8OOeSQgu1t376d\nqVOnZrZs2fJ2FEUDC43Ly0ngTjUTxsuV/OTcc89NTNwA3/nOd+Too49GRB6LW2JKhohMV9WxN910\nUyLiBujatSu33nprVefOnfuJyG8SMVomOCXwIAiWDBgwQM8555zEbV911VVB586de4jI/MSNZ4mI\nHAJcd95550mu7xTt0bVrV2bNmhWo6iQROSVR4yXEGYGLyPmqOvjaa6+tan/v3OnSpQvTp08PVPXf\nRKRfGmW0h4gs7Nevn5599tmp2B85ciS1tbVREAR3p1JACXBG4EEQXDV8+HDt3bv3ftvr6+v3/tTV\n1bF06dK8yxgzZgw9e/bMAD8q0N2cEZFuqjp26tSpex/g+vp6LrnkEpYuXcqOHTv2brvxxhupq6tr\n1VZbfPOb3wyiKOoXdxBVPE4IXEQGR1E0ZMqUKfudz8KFCxERampqqKmpYdiwYWzYsIH6+vwX8z3z\nzDOrRGRywU7nzjXV1dXRMcccs3dD37596d+/P2PGjKG6uhqAmpoaxo0bx7Bhw/IqpHfv3gwaNCgj\nIjcm4nWJcULgwEVdu3bd0zQuXbx4MQcddBCHH374fjsOGjQo79oNYOLEiahqFxEp6iSJIAjOqK2t\n3S/8qq+vp6amhubbmp9zrnz605+u8jV4GSEiJwwaNGjvzd+xYweLFy/mlFPe/660bdu2Fm28/vrr\nNDQ0tFtWp06dqK6u3gNMyt/j3FHVAbW1tftta0ngO3fu3FubNyXb8wMYO3YsURQdICLvN1RhOCHw\nIAg+esQRR+xtvlu3bh2ttTK0JAqA6urqrGv2/v37B9gFZYtCHH93GTNmzH7bX3755axr61zOr2/f\nvgRBEAGfytXXcsMJgatqdb9++xo2du7cSd++fd+3X0NDA5s2bWoxPu3atWvW5fXq1SsQkT75eZsX\n/QC6d+++d0NDQ8P7zrGhoYE+fVp2K5fzA+jcubMC/XP0s+zoVGoHkqJp/8vhhx/OokWLAPvVDLZW\nWrJkCdOmTQOsGOrq6ujbty8HHHAAAwcOZN26ddTX17N+/XpOOeWUvS+jjZ+bliUixawc3ldWS4Jd\nv3793oe3kPNrQipNrsXECYGLyK7NmzfvvePV1dVMnjyZJUuWUF1dzY4dO3j99dc566yz9h7zhz/8\ngQsuuACAO+64g4EDBzJw4EBqamrYuHEjdXV1NDQ0sG3bNoYPH75fedu2bSOKooLHbORAA9ixMY0D\nqqqrqznuuONYsmTJ3pq86TdTIecHkMlkpLHcSsYJgUdRtGb16tUjgL3V+IABA2gcNtoYezbWWM1j\n8MY25EYaa8dGwTQVC8C6dev2UMSJyqq6KQiCzLPPPlvVNA5vHpO3Rq7nt337dvbs2RMADxfoeslx\nJQZ/eu3ata0OEOrfvz8LFiygrq5ur7hPP/10li5dSl1dHaeeeipga+b6+noaGhoYNmzY3lqu+QOx\ndevWKuCB5uWkSRAEb61cmf04r0LOb9myZYjIblWt+BrcidGEIjIGeOauu+7iwAMPTMkzy1NPPcXs\n2bMjoFMxJ0KIyD19+vQ5a968eanHxbNmzYpefPHF+kwmc2Q+x/vRhAmjqkuDINh4++23p17W/Pnz\nMyLypxLM8pmxcePGqrVr16ZayK5du3jhhRckiqL/TLWgIuGEwAGiKPrpkiVLol27dqVWxrp161i/\nfn2Vqn4rtUJaQVX/HgTByzfffHOUZjnz5s1DRHaoqhPDZp0ROPD9TCazPQ4fUuGKK67IBEHwnKr+\nJa0y2iKKoi/97W9/C5YtW9b+znnw5ptv8uCDD2oURTNSKaAEOCNwVc2o6sTly5cHK1Yk38Axd+5c\nNm3apFEUjUvceJao6tMi8uvrr78+2rlzZ6K2oyhi1qxZGRF5UVWdybbrjMABVHWxiMy7+uqrtZAR\ng8156KGH+P3vf4+qXqSqmxIznAeqet6ePXs2TJs2LZNkODZr1qzMhg0bMlEUnZqY0TLAKYGD/RoH\nHp4+fbomUZPffffd3HzzzQCXqeq8gg0WiKpqFEU1mzZt2vKVr3wls2lTYc/brl27+Pa3v5156aWX\nVFVHpZ0UqNg40UzYEkEQ3K6qXxw9enQ0c+bMnFIqgO3qvvzyyzPr168X4OuqenOhPiWJiPQIgqAO\n+MiXv/xlmThxYs42li9fzuzZs6Pdu3dvVdVjVXV1Qr6VTTOhswKPbZ0iIg9UVVV1Pfnkk4MLL7yw\n3XbydevWMWfOnOj5558XEVkTRdE4Vf17Ev6kgYh8D7iid+/e0aRJk6rOOOMMgqDtL+annnqK+fPn\nZ9avX18F3AOcl2SGKy/wNkj64ohIFXB5EATfjKKoT8+ePfcMGTKk6ogjjpBevXoRRRFvv/02q1ev\n1tdeey2zffv2TkEQrImi6D9V9c6k/EgTEekP/FhEzhSRTgcffHDm8MMP73TYYYdRXV3Nrl27ePXV\nV6mvr8+8+eabsnv3bhGRx1X122m0CHmBt0HKyTePAc4VkU8EQXCoqnaJt2+PouglVX0CuE1VX02j\n/LSJU1p8ATg9CIJjsMNsOwEZEdmUyWSexSbgnJtUFqtW/PACb41yujie/Cine+hcK4rH0xQvcI/T\neIF7nMYAprDBAAACgklEQVQL3OM0TgjcGLPAGDMlz2MvNcZszPd4T3njhMCBa7HNXzkThuH12BXX\nyqs5yZMITszJDMPQ5ZXSPAVQEQI3xkwFfgncG28aCswIw3CRMeZSYCbwgzAMbzDGzAamA48BM4DF\nwPIwDD9ljDkLu0Z9H2BFGIY3FPtcPMWlIgQehuGtsTiHhGE42hgzAlhljOkVhuH1xphRxCFGGIYz\njDE9gVHAZuCeMAwvNsYMBa4Lw/AwgDjufrSt2t8YMwT7QACsiX8PDcNwTjpn6kmaSovBfwMQhuFz\n8eemCSL39pyFYTgN6A08GobhxfHmswCMMdcZY67Dpn0YmkWZa4DaMAwXYRdWnVDQGXiKSkXU4E1o\nq/u3+Uvio8AUY8yI+IHoC2wOw3BmtoWFYbjWGDMNGx6Brc2X5+Kwp7RUUg0uxLWnMaYWK+iVTf63\nV/xxXL4QG5svjDc/CoyMww6MMWcZY8a1dHwzxjUJY84G5sQhkqcCqKQaXIE1xpgF2NBiQhiG/4xj\n83HAYGPMY8A5wCXYh2EzMNQYswIrzhnAQmPMRuCRMAzvjY+vBXoZY1Y2CX8aeaXJ32uwYVHJV1rz\nZEfFjCY0xjwCLAjDcG4J3PLkgB9NmB9thREeT4tUhMCbhBFTffzryYWKCVE8lUM53cOKqME9nnwp\ny1YUESmvrxVPxVJ2IYrHkyQ+RPE4jRe4x2m8wD1O4wXucRovcI/TeIF7nMYL3OM0XuAep/EC9ziN\nF7jHabzAPU7jBe5xGi9wj9N4gXucxgvc4zRe4B6n8QL3OI0XuMdpvMA9TuMF7nEaL3CP03iBe5zG\nC9zjNF7gHqfxAvc4jRe4x2m8wD1O4wXucZr/D5u0Mi1AdqpnAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x10950ba58>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "from matplotlib import rc\n", | |
| "rc(\"font\", family=\"serif\", size=12)\n", | |
| "rc(\"text\", usetex=True)\n", | |
| "\n", | |
| "pgm = daft.PGM([3, 4.5], origin=[0, 0])\n", | |
| "pgm.add_node(daft.Node(\"alpha\", r\"$\\alpha_\\phi$\", 2, 4))\n", | |
| "pgm.add_node(daft.Node(\"BT\", r\"$B_\\perp$\", 1, 3))\n", | |
| "pgm.add_node(daft.Node(\"phi\", r\"$\\phi$\", 2, 3))\n", | |
| "pgm.add_node(daft.Node(\"QSyn\", r\"$Q_\\mathrm{syn}$\", 1, 2))\n", | |
| "pgm.add_node(daft.Node(\"USyn\", r\"$U_\\mathrm{syn}$\", 2, 2))\n", | |
| "pgm.add_node(daft.Node(\"QObs\", r\"$Q_\\mathrm{obs}$\", 1, 1, observed=True))\n", | |
| "pgm.add_node(daft.Node(\"UObs\", r\"$U_\\mathrm{obs}$\", 2, 1, observed=True))\n", | |
| "pgm.add_plate(daft.Plate([0.5, 0.4, 2, 3.2], label=r\"pixel $n$\"))\n", | |
| "pgm.add_edge(\"BT\", \"QSyn\")\n", | |
| "pgm.add_edge(\"BT\", \"USyn\")\n", | |
| "pgm.add_edge(\"phi\", \"QSyn\")\n", | |
| "pgm.add_edge(\"phi\", \"USyn\")\n", | |
| "pgm.add_edge(\"alpha\", \"phi\")\n", | |
| "pgm.add_edge(\"QSyn\", \"QObs\")\n", | |
| "pgm.add_edge(\"USyn\", \"UObs\")\n", | |
| "pgm.render()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The Bayesian model that we implement is given by the previous graphical model. Our aim is to infer the full posterior distribution $p(\\mathbf{B_\\perp},\\phibold$\n", | |
| "\n", | |
| "\n", | |
| "We assume that the observed Stokes $Q$ and $U$ of each pixel can be explained using the weak-field approximation. Therefore, each pixel can potentially have a different value of $B_\\perp$ and azimuth $\\phi$. Our aim is to get information about the misalignment between the azimuth inferred from the fibrils and that inferred from the linear polarization. To this end, we first rotate the observed Stokes parameter in each pixel so that the Stokes $Q$ reference direction is aligned with the direction of the fibril. We do it using the standard rotation:\n", | |
| "\n", | |
| "$$\n", | |
| "Q' = Q \\cos 2 \\phi_\\mathrm{fib} + U \\sin 2 \\phi_\\mathrm{fib} \\\\\n", | |
| "U' = -Q \\sin 2 \\phi_\\mathrm{fib} + U \\cos 2 \\phi_\\mathrm{fib}.\n", | |
| "$$\n", | |
| "\n", | |
| "If the magnetic field is aligned with the direction of the fibril, then only Stokes $Q'$ should be nonzero, and the inferred azimuth should be zero. This might not be the case in, at least, two situations: first, when there is a misalignment between the magnetic field and the direction of the fibril and, second, when the azimuth inferred from the Stokes parameters has an uncertainty inherited from the uncertainty in the measurement of the Stokes parameters.\n", | |
| "\n", | |
| "When the spectral lines are observed, we come out with a discretized version of the generative model. We assume that we observe the Stokes parameters at a set of fixed $N$ wavelength points. Our generative model for the rotated Stokes parameters is the following:\n", | |
| "\n", | |
| "$$\n", | |
| "Q'_i = \\beta B_\\perp^2 \\hat{I}_i \\cos 2\\phi + e_i \\\\\n", | |
| "U'_i = \\beta B_\\perp^2 \\hat{I}_i \\sin 2\\phi + f_i,\n", | |
| "$$\n", | |
| "where $\\hat{I}=\\frac{\\partial I(\\lambda)}{\\partial \\lambda} \\frac{1}{\\lambda_0-\\lambda}$ has been used to simplify the notation. The terms $e_i$ and $f_i$ refer to the presence of measurement noise, that we assume to follow two Gaussian distributions with zero mean and variance $\\sigma_n^2$.\n", | |
| "\n", | |
| "The likelihood or sampling distribution for each pixel can be written as:\n", | |
| "\n", | |
| "$$\n", | |
| "p(D_i|B_\\perp,\\phi) = p(\\mathbf{Q}|B_\\perp,\\phi) p(\\mathbf{U}|B_\\perp,\\phi) = \\prod_{i=1}^N \\frac{1}{2\\pi\\sigma_n^2} \\exp \\left[ - \\frac{(Q_i-\\beta B_\\perp^2 \\hat{I}_i \\cos 2\\phi)^2}{2 \\sigma_n^2}\\right] \\exp \\left[ - \\frac{(U_i-\\beta B_\\perp^2 \\hat{I}_i \\sin 2\\phi)^2}{2 \\sigma_n^2}\\right]\n", | |
| "$$\n", | |
| "\n", | |
| "This likelihood can be simplified if the products are moved inside the exponential as sums:\n", | |
| "\n", | |
| "$$\n", | |
| "p(D|B_\\perp,\\phi) = \\left[ \\frac{1}{2\\pi\\sigma_n^2} \\right]^N \\exp \\left[ - (A_Q + B_Q B_\\perp^2 \\cos 2\\phi)^2}{2 \\sigma_n^2}\\right] \\exp \\left[ - \\frac{(U_i-\\beta B_\\perp^2 \\ddot{I}_i \\sin 2\\phi)^2}{2 \\sigma_n^2}\\right]\n", | |
| "$$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.4.3" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 0 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment