{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resolution Convolution (4D) Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from neutronpy import Instrument, Sample\n",
    "from neutronpy.instrument.tools import _modvec, _star\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matplotlib version:  2.0.0\n",
      "numpy version:  1.12.1\n",
      "neutronpy version:  1.0.3\n"
     ]
    }
   ],
   "source": [
    "import matplotlib as mpl\n",
    "import neutronpy as npy\n",
    "\n",
    "print('matplotlib version: ', mpl.__version__)\n",
    "print('numpy version: ', np.__version__)\n",
    "print('neutronpy version: ', npy.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def angle2(x, y, z, h, k, l, lattice):\n",
    "    V, Vstar, latticestar = _star(lattice)\n",
    "\n",
    "    return np.arccos(2 * np.pi * (h * x + k * y + l * z) / _modvec([x, y, z], lattice) / _modvec([h, k, l], latticestar))\n",
    "\n",
    "def sqw(H, K, L, W, p=None):\n",
    "    '''S(Q, w) for a gapped excitation in a 1D antiferromagnet'''\n",
    "    \n",
    "    dx, dy, dz, cc, gamma = p[:5]\n",
    "    \n",
    "    omega_x = np.sqrt(cc ** 2 * np.sin(2 * np.pi * H) ** 2 + dx ** 2)\n",
    "    omega_y = np.sqrt(cc ** 2 * np.sin(2 * np.pi * H) ** 2 + dy ** 2)\n",
    "    omega_z = np.sqrt(cc ** 2 * np.sin(2 * np.pi * H) ** 2 + dz ** 2)\n",
    "    \n",
    "    lor_x = 1 / np.pi * gamma / ((W - omega_x) ** 2 + gamma ** 2)\n",
    "    lor_y = 1 / np.pi * gamma / ((W - omega_y) ** 2 + gamma ** 2)    \n",
    "    lor_z = 1 / np.pi * gamma / ((W - omega_z) ** 2 + gamma ** 2)\n",
    "    \n",
    "    sqw = np.array([lor_x * (1 - np.cos(np.pi * H)) / omega_x / 2, \n",
    "                    lor_y * (1 - np.cos(np.pi * H)) / omega_y / 2, \n",
    "                    lor_z * (1 - np.cos(np.pi * H)) / omega_z / 2])\n",
    "\n",
    "    return sqw\n",
    "    \n",
    "def pref(H, K, L, W, EXP, p=None):\n",
    "    '''More complicated prefactor'''\n",
    "    I, bgr = p[5:]\n",
    "    \n",
    "    \n",
    "    sample, rsample = EXP.get_lattice()\n",
    "    q2 = _modvec([H, K, L], rsample) ** 2\n",
    "    sd = q2 / (16 * np.pi ** 2)\n",
    "    ff = 0.0163 * np.exp(-35.883 * sd) + 0.3916 * np.exp(-13.223 * sd) + 0.6052 * np.exp(-4.339 * sd) - 0.0133\n",
    "    \n",
    "    \n",
    "    # Calculate the polarization factors for transverse excitations\n",
    "    alphay = angle2(0, 1, 0, H, K, L, sample)\n",
    "    alphaz = angle2(0, 0, 1, H, K, L, sample)\n",
    "    alphax = angle2(1, 0, 0, H, K, L, sample)\n",
    "\n",
    "    # Polarization factors for each of the three modes.\n",
    "    polx = np.sin(alphax) ** 2\n",
    "    poly = np.sin(alphay) ** 2\n",
    "    polz = np.sin(alphaz) ** 2\n",
    "\n",
    "    prefactor = np.array([ff ** 2 * polx * I, ff ** 2 * poly * I, ff ** 2 * polz * I])\n",
    "\n",
    "    # Constant Background\n",
    "    bgr = np.ones(H.shape) * bgr\n",
    "    \n",
    "    return np.ones(H.shape)[np.newaxis, ...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sample = Sample(4., 4., 4., 90, 90, 90, mosaic=60., u=[1, 0, 0], v=[0, 1, 0])\n",
    "\n",
    "instr = Instrument(efixed=14.7, sample=sample, hcol=[50, 80, 50, 120], ana='PG(002)', mono='PG(002)', \n",
    "                   moncor=1, mono_mosaic=35., ana_mosaic=35.)\n",
    "\n",
    "instr.mono.dir = 1\n",
    "instr.sample.dir = -1\n",
    "instr.ana.dir = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "H1=1.5 \n",
    "K1=0\n",
    "L1=0.35\n",
    "W1=np.arange(20, 0, -0.5)\n",
    "\n",
    "q = [H1, K1, L1, W1]  # q = [2, -0.18, 0, eValues]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = [3, 3, 3, 30, 0.4, 6e4, 40]\n",
    "\n",
    "output_fix = instr.resolution_convolution(sqw, pref, 1, q, METHOD='fix', ACCURACY=[5,5], p=p)  # Fixed sample method\n",
    "\n",
    "output_mc = instr.resolution_convolution(sqw, pref, 1, q, METHOD='mc', ACCURACY=[5], p=p)  # Monte Carlo Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11596f470>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+UXWV97/H3xyQkmsjokEQ0kDshE+kNJUaZoEbtaqBg\noK6MbWgJVcvV3MVqGxS1eIGry6aodznVC1INeFNAMcWGHylmBORHIV69opAJhECCqQeIMAFMwtBo\ntJAffO8fe09yMnPmnJPZs8+ZM/N5rTVrznmeZ+/zPXtmznf28+z9PIoIzMzMBus19Q7AzMwamxOJ\nmZll4kRiZmaZOJGYmVkmTiRmZpaJE4mZmWXiRGJmZpk4kZiZWSZOJGZmlsnYegdQC5MnT46WlpZ6\nh2Fm1jA2bNiwKyKmVNN2VCSSlpYWurq66h2GmVnDkPTLatu6a8vMzDJxIjEzs0ycSMzMLBMnEjMz\ny8SJxMzMMhkVV22NNGs3bmfFugKFHXtonTqJZQtaaZ87rep6M7Oh5ETSYNZu3M5X79lKx+I5zGtp\nZv22Hi5ZswmA9rnTKtabmQ01J5IGs2JdgY7Fc5g/czIA82dOpmPxHJZ3bqZ97jTed8dptO97AVYl\n7ecDPwZ67jgW5m6tW9xmNnI5kTSYwo49zGtpPqxsXkszhR17AGje9wIs391vu+blTTWJz8xGHw+2\nN5jWqZNYv63nsLL123ponTqpThGZ2WjnM5IG07n/r5mwavthZfOBzonTgC11icnMRjcnkgYz4bfb\nS3ZdTXDXlZnVibu2zMwsE5+RjDRN06HU2UnT9NrHYmajghPJSPOpx+odgZmNMu7aMjOzTHxG0mjc\ndWVmw4wTSaNx15WZDTPu2jIzs0ycSMzMLBMnEjMzy8SJxMzMMsk1kUhaKGmrpIKkS0vUj5d0U1r/\noKSWorrL0vKtkt5fVP4GSbdK+rmkJyS9O8/3YGZm5eWWSCSNAVYAZwGzgfMkze7TbCnwUkS0AlcC\nHem2s4ElwEnAQuDqdH8AVwF3RcTvAW8DnsjrPZiZWWV5npGcChQi4qmI2AusBtr7tGkHbkgf3wqc\nLklp+eqIeCUingYKwKmSjgb+ALgOICL2RsR/5PgezMysgjwTyTTg2aLn3WlZyTYRsR/YDRxTZtsT\ngJ3AtyQ9IulaSRNLvbikCyR1SerauXPnULwfMzMrIc9EohJlUWWbgcrHAu8AromItwO/BfqNvQBE\nxMqIaIuItilTplQftZmZHZE8E0k3cHzR8+OA5wZqI2ks0AT0lNm2G+iOiAfT8ltJEouZmdVJnolk\nPTBL0gxJR5EMnnf2adMJnJ8+Pge4PyIiLV+SXtU1A5gFPBQRLwDPSjox3eZ0vCygmVld5TbXVkTs\nl3QhcDcwBrg+IjZLuhzoiohOkkHzVZIKJGciS9JtN0u6mSRJ7AeWRcSBdNcfB25Mk9NTwEfzeg9m\nZlaZkhOAka2trS26urrqHYaZWcOQtCEi2qpp6zvbzcwsEycSMzPLxInEzMwycSIxM7NMnEjMzCwT\nJxIzM8vEicTMzDJxIjEzs0ycSMzMLJPcpkixYejKk2H3M/3Lm6bDpx6rfTxmNiI4kYwmu5+B5bv7\nly9vqn0sZjZiuGvLzMwycSIxM7NMnEjMzCwTJxIzM8vEg+2jSdP00gPrTdNrH4uZjRhOJKOJL/E1\nsxw4kYwyazduZ8W6AoUde2idOollC1ppnzut3mGZWQNzIhlF1m7czlfv2UrH4jnMa2lm/bYeLlmz\nCcDJxMwGzYPto8iKdQU6Fs9h/szJjBvzGubPnEzH4jmsWFeod2hm1sB8RjLc5DiNSWHHHua1NB9W\nNq+lmcKOPZn2a2ajW66JRNJC4CpgDHBtRHy5T/144DvAKcCLwLkRsS2tuwxYChwAPhERd6fl24Df\npOX7I6Itz/dQczlOY9I6dRLrt/Uwf+bkg2Xrt/XQOnVS5n2b2eiVW9eWpDHACuAsYDZwnqTZfZot\nBV6KiFbgSqAj3XY2sAQ4CVgIXJ3ur9eCiJg74pJIzpYtaOWSNZt44Mld7DvwKg88uYtL1mxi2YLW\neodmZg0szzOSU4FCRDwFIGk10A5sKWrTDixPH98KfEOS0vLVEfEK8LSkQrq/n+YY74jXO6C+vHPz\nwau2Lj7zRA+0m1kmeSaSacCzRc+7gXcO1CYi9kvaDRyTlv+sz7a9n3YB3CMpgP8TEStziH3Eap87\nzYnDzIZUnolEJcqiyjbltn1PRDwnaSpwr6SfR8SP+r24dAFwAcD06b5z28wsL3kmkm7g+KLnxwHP\nDdCmW9JYoAnoKbdtRPR+3yHpNpIur36JJD1TWQnQ1tbWN4ENWy9PnMaEEgPrL0+cxoQ6xGNmVkme\n95GsB2ZJmiHpKJLB884+bTqB89PH5wD3R0Sk5UskjZc0A5gFPCRpoqTXA0iaCJwJPJ7je6i5CZ/Z\nwtoPbuHMpk5OeOW7nNnUydoPbmHCZ7ZU3tjMrA5yOyNJxzwuBO4mufz3+ojYLOlyoCsiOoHrgFXp\nYHoPSbIhbXczycD8fmBZRByQ9CbgtmQ8nrHAdyPirrzeQ714HMPMGomSE4CRra2tLbq6uuodhplZ\nw5C0odpbLDxFipmZZeJEYmZmmTiRmJlZJk4kZmaWiROJmZll4kRiZmaZOJGYmVkmTiRmZpaJE4mZ\nmWXiRGJmZpk4kZiZWSZOJGZmlokTiZmZZeJEYmZmmTiRmJlZJk4kZmaWiROJmZll4kRiZmaZOJGY\nmVkmTiRmZpaJE4mZmWWSayKRtFDSVkkFSZeWqB8v6aa0/kFJLUV1l6XlWyW9v892YyQ9Iun2POM3\nM7PKckskksYAK4CzgNnAeZJm92m2FHgpIlqBK4GOdNvZwBLgJGAhcHW6v14XAU/kFbuZmVUvzzOS\nU4FCRDwVEXuB1UB7nzbtwA3p41uB0yUpLV8dEa9ExNNAId0fko4D/hi4NsfYzcysSnkmkmnAs0XP\nu9Oykm0iYj+wGzimwrZfA/4H8OrQh2xmZkcqz0SiEmVRZZuS5ZI+AOyIiA0VX1y6QFKXpK6dO3dW\njtbMzAalqkQiqXkQ++4Gji96fhzw3EBtJI0FmoCeMtu+B1gkaRtJV9lpkv651ItHxMqIaIuItilT\npgwifDMzq0a1ZyQPSrpF0tnpGEY11gOzJM2QdBTJ4HlnnzadwPnp43OA+yMi0vIl6VVdM4BZwEMR\ncVlEHBcRLen+7o+ID1cZj5mZ5aDaRPJWYCXwEaAg6X9Jemu5DdIxjwuBu0musLo5IjZLulzSorTZ\ndcAxkgrAp4FL0203AzcDW4C7gGURceDI3pqZmdWCkhOAI9hAWgD8MzAReBS4NCJ+mkNsQ6atrS26\nurrqHYaZWcOQtCEi2qppO7bKHR4DfJjkjORXwMdJup/mArcAMwYXqpmZNbqqEgnwU2AV8MGI6C4q\n75L0zaEPy8zMGkW1YySfi4gvFCcRSX8GEBEduURmZmYNodpE0m+eLOCyoQzEzMwaU9muLUlnAWcD\n0yT9Y1HV0cD+PAMzM7PGUGmM5DmgC1gEFN9N/hvgU3kFZfWzduN2VqwrUNixh9apk1i2oJX2uX1n\ntjEzO6RsIomIR4FHJd2Y3hdiI9jajdv56j1b6Vg8h3ktzazf1sMlazYBOJmY2YAqdW3dHBF/Djwi\nqfiGEwEREXNyjc5q6n13nEb7vheS6/OA+cCPgZ47joW5W+sZmpkNY5W6ti5Kv38g70Cs/pr3vQDL\nd/cvX95Uh2jMrFGUvWorIp5PH+4Cno2IXwLjgbfRfwJGMzMbhaq9/PdHwARJ04D7gI8C384rKDMz\naxzVJhJFxO+APwW+HhF/QrJ8rpmZjXJVJxJJ7wY+BNyRllU7vYqZmY1g1SaDi0juZL8tnQr+BGBd\nfmFZXTRNh1ID603Tax+LmTWMI55GvhF5GnkzsyOTxzTybwUuBlqKt4mI0wYToJmZjRzVdm3dAnwT\nuBbwSoVmZnZQtYlkf0Rck2skZmbWkKq9auv7kv5G0pslNfd+5RqZmZk1hGrPSM5Pv3+mqCyAE4Y2\nnFHgypNh9zP9y5umw6ceq308ZmYZVZVIIsJrsg+V3c+UnM+q5GW3ZmYNoKquLUmvk/Q5SSvT57Mk\nVZzIUdJCSVslFST1W2VR0nhJN6X1D0pqKaq7LC3fKun9adkESQ9JelTSZkl/X+0bNTOzfFQ7RvIt\nYC/JzOIA3cAXy20gaQywAjiLZDqV8yT1nVZlKfBSRLQCVwId6bazgSXAScBC4Op0f68Ap0XE24C5\nwEJJ76ryPZiZWQ6qTSQzI+IfgH0AEfGfJGuSlHMqUIiIpyJiL7AaaO/Tph24IX18K3C6JKXlqyPi\nlYh4GigAp0ZiT9p+XPo18u+oNDMbxqpNJHslvZb0Q1vSTJKzg3KmAc8WPe9Oy0q2SVdg3A0cU25b\nSWMkbQR2APdGxINVvgczM8tBtVdtLQfuAo6XdCPwHpKp5MspdcbS9+xhoDYDbhsRB4C5kt4A3Cbp\n9yPi8X4vLl0AXAAwffrwmStq55ipTCkxsL5zzFSm1CEeM7Osqr1q6x5JG4B3kXzIXxQRuyps1g0c\nX/T8OPovhtXbplvSWKAJ6Klm24j4D0k/JBlD6ZdIImIlsBKSubYqxFozH5p0LcsXncT8mZMPlj3w\n5C6Wd27mnjrGZWY2WNVetXVfRLwYEXdExO0RsUvSfRU2Ww/MkjRD0lEkg+edfdp0cugelXOA+yOZ\nRbITWJJe1TUDmAU8JGlKeiZC2tX2R8DPq3kPw8WyBa1csmYTDzy5i30HXuWBJ3dxyZpNLFvQWu/Q\nzMwGpewZiaQJwOuAyZLeyKEup6OBt5TbNiL2S7oQuBsYA1yfTkF/OdAVEZ3AdcAqSQWSM5El6bab\nJd0MbAH2A8si4oCkNwM3pFdwvQa4OSJuH9Q7r5P2uckw0fLOzRR27KF16iQuPvPEg+VmZo2m7DTy\nki4CPkmSNLZzKJH8GviniPhG7hEOAU8jb2Z2ZIZsGvmIuAq4StLHI+LrQxKdmZmNKNUOtn9d0nz6\nr0fynZziMjOzBlHtwlargJnARg6tRxKAE4mZ2ShX7X0kbcDsGA3r8pqZ2RGp9s72x4Fj8wzEzMwa\nU7VnJJOBLZIeomhqlIhYlEtUZmbWMI5kihQzM7N+qr1q6//mHYiZmTWmSne2/4bS07QLiIg4Opeo\nzMysYVS6IfH1tQrEzMwaU7VXbZmZmZXkRGJmZpk4kZiZWSZOJGZmlkm195GYAbB243ZWrCscXEtl\n2YJWr6ViNso5kVjV1m7czlfv2UrH4jnMa2lm/bYeLlmzCcDJxGwUcyKxqr3vjtNo3/cCrEqezwd+\nDPTccSzM3VrP0MysjpxIrGrN+16A5bv7ly9vqkM0ZjZceLDdzMwycSIxM7NMnEjMzCwTJxIzM8sk\n18F2SQuBq4AxwLUR8eU+9eNJ1n0/BXgRODcitqV1lwFLSdaI/0RE3C3p+LT9scCrwMqIuCrP92BF\nmqZDqYH1pum1j8XMho3cEomkMcAK4AygG1gvqTMithQ1Wwq8FBGtkpYAHcC5kmYDS4CTgLcA/ybp\nrcB+4G8j4mFJrwc2SLq3zz4tL596rN4RmNkwlGfX1qlAISKeioi9wGqgvU+bduCG9PGtwOmSlJav\njohXIuJpoACcGhHPR8TDABHxG+AJwHfCmZnVUZ6JZBrwbNHzbvp/6B9sExH7gd3AMdVsK6kFeDvw\nYKkXl3SBpC5JXTt37hz0mzAzs/LyTCQqUdZ3tcWB2pTdVtIkYA3wyYj4dakXj4iVEdEWEW1Tpkyp\nMmQzMztSeSaSbuD4oufHAc8N1EbSWKAJ6Cm3raRxJEnkxoj411wiNzOzquWZSNYDsyTNkHQUyeB5\nZ582ncD56eNzgPsjItLyJZLGS5oBzAIeSsdPrgOeiIgrcozdzMyqlNtVWxGxX9KFwN0kl/9eHxGb\nJV0OdEVEJ0lSWCWpQHImsiTddrOkm4EtJFdqLYuIA5LeC3wEeEzSxvSl/mdE3JnX+zAzs/KUnACM\nbG1tbdHV1VXvMMzMGoakDRHRVk1b39luZmaZOJGYmVkmTiRmZpaJF7bKgdc1N7PRxIlkiHldczMb\nbdy1NcRWrCvQsXgO82dOZtyY1zB/5mQ6Fs9hxbpCvUMzM8uFz0iGWGHHHua1NB9WNq+lmcKOPXWK\nqIauPBl2P9O/vGm6Zw42G8GcSIZY69RJrN/Ww/yZkw+Wrd/WQ+vUSXWMqkZ2PwPLd/cvL7WGiZmN\nGO7aGmLLFrRyyZpNPPDkLvYdeJUHntzFJWs2sWxBa71DMzPLhc9IhljvgPryzs0Hr9q6+MwTPdBu\nZiOWE8lQu/Jk2nc/k6zgNZ5khZXvAes8TmBmI5MTyVDzOIGZjTJOJDZkesYdS3OJhNkz7liaS7Q3\ns5HBicSGTNueK9j6xbMYN+bQNRz7DrxK2+d+wFN1jMvM8uWrtmzI9F76XGzUXPpsNoo5kdiQ8aXP\nZqOTu7aGWtP00gPrTdNrH0uN+dJns9HJKySamVk/XiHRzMxqxonEzMwy8RiJ1Y5nBzYbkXJNJJIW\nAlcBY4BrI+LLferHA98BTgFeBM6NiG1p3WXAUuAA8ImIuDstvx74ALAjIn4/z/htiPmuf7MRKbeu\nLUljgBXAWcBs4DxJs/s0Wwq8FBGtwJVAR7rtbGAJcBKwELg63R/At9MyMzMbBvIcIzkVKETEUxGx\nF1gNyVyGRdqBG9LHtwKnS1JavjoiXomIp4FCuj8i4kdAD2ZmNizkmUimAc8WPe9Oy0q2iYj9JHPl\nHlPltmZmNgzkmUhUoqzvTSsDtalm2/IvLl0gqUtS186dO49kUzMzOwJ5DrZ3A8cXPT8OeG6ANt2S\nxgJNJN1W1WxbVkSsBFZCckPiEUVuuXh54jQmlBhYf3niNCbUIR4zGxp5npGsB2ZJmiHpKJLB884+\nbTqB89PH5wD3R3KrfSewRNJ4STOAWcBDOcZqNbBo7DU88JEnkyu30q8HPvIki8ZeU+/QzCyD3BJJ\nOuZxIXA38ARwc0RslnS5pEVps+uAYyQVgE8Dl6bbbgZuBrYAdwHLIuIAgKR/AX4KnCipW9LSvN6D\nDa3Cjj3Mazl8ZZJ5Lc0UduypU0RmNhRyvY8kIu4E7uxT9vmixy8DfzbAtl8CvlSi/LwhDnNQ1m7c\nzop1hYOTEy5b0OrJCSvonWZ+/szJB8s8zbxZ4/Od7YPw8ldm0/7b7f3WZX/53mlM+MyW+gY3jPVO\nM9+xeA7zWppZv62HS9Zs4uIzT6x3aGaWgRPJIEz47faSd2iXGki2Q8pOM+/pU8walhOJ1VT73Gml\nuwA9fYpZw3IiscbgMxazYcuJxBqDz1jMhi2vR2JmZpn4jGQwRvG67GZmfTmRDIb75IfczjFTmVIi\nOe8cM5UpdYjHzKrnRGLDwocmXcvyRScddrPiA0/uYnnnZu6pZgcejDerGyeSUvyhVHMVb1as1J3o\nwXizunEiKcUfSjVX9mZFYO2Cu/jqPVv7J5oFJ/ZbLa0k/3NglhsnEhs2BrxZEVixrkDH4jkHu77m\nz5xMx+I5LO/cXN0cZ/7nwCw3vvzXGoJnDjYbvnxGYg0h95mD3fVlNmhOJNYQKg3G94w7luYS3VQ9\n446luV9pCRW6vrxsgNnAnEhK8Q2Hw06lwfi2PVew9YtnMW7Mod7afQdepe1zP+CpjK9dzbIBlRJN\nnvX1fO2hqLfG50RSirsyhqVyg/GVur6y3PBYadmASommUv3ajdtLXpHW+57L1QOD3na41A/XJOcE\negQiYsR/nXLKKWEj2/ce6Y73dtwXPynsjL37D8RPCjvjvR33xfce6Y6IiDOu+GH8pLDzsG1+UtgZ\nZ1zxw+TJ3x1desd/d3T5ukrbVlFfKbZy9Vm2rXd9pZ9ZPevzfu3eNmdc8cOYcentB49HsXL1Wbat\nFtAVVX7G1v1DvhZfTiSjQ6U/vLJ/2HVMJDsubz30OkVfOy5vrVhfadsZl94ee/cfOOxl9+4/EDMu\nvb3u9cM5yeX92vVMctU6kkTiri0bMcp1fVUaY3l54rSSK1y+PHFa0rWVoykHdpTsOuvtiqtUX66u\na9KnGfeFvzisbhzQNelY4I8z11fqUixXX+mS7nrX57nvSvdFlasHBr1tXl1rTiQ2apRLNL1jFSX7\ntBv4psXmfS+UTDS9V7hlre/c/9dMWHV4op0PdE6cBmwpW/+7SQdyTXJZ6pe8buWgE2Q19fVMcnnI\nNZFIWghcBYwBro2IL/epHw98BzgFeBE4NyK2pXWXAUuBA8AnIuLuavZpNlgDJhpfxTegShcilKuf\nALkmuSz1WRJkNfXZkiCZEmgeckskksYAK4AzgG5gvaTOiNhS1Gwp8FJEtEpaAnQA50qaDSwBTgLe\nAvybpLem21Tap9nQqnQVX6VE40TUcLIkyGrqsybBTNvmIM8zklOBQkQ8BSBpNdAOFH/otwPL08e3\nAt+QpLR8dUS8AjwtqZDujyr2aVZblRKNLye3ES7PRDINeLboeTfwzoHaRMR+SbuBY9Lyn/XZtrfP\nodI+zRpL1jMan+1YneWZSFSiLKpsM1B5qUkm++4z2bF0AXABwPTp/qOyYSzPM5qsScrdclaFPBNJ\nN3B80fPjgOcGaNMtaSzQBPRU2LbSPgGIiJXASoC2traSycZsxMuapPIeH8ozyWWpLzWBpw0oz0Sy\nHpglaQawnWTw/C/6tOkEzgd+CpwD3B8RIakT+K6kK0gG22cBD5GcqVTap5nVSp5nU3knwXL1V548\nfJMcDLuzSCU3MOa0c+ls4Gskl+peHxFfknQ5yR2TnZImAKuAt5OciSwpGkj/LPAxYD/wyYj4wUD7\nrBRHW1tbdHV1Df0bNDMboSRtiIi2qtrmmUiGCycSM7MjcySJxCskmplZJk4kZmaWiROJmZll4kRi\nZmaZjIrBdkk7gV8OUD0Z2FXDcI6EYxscxzY4jm1wRmps/yUiKi0gCoySRFKOpK5qr0yoNcc2OI5t\ncBzb4Dg2d22ZmVlGTiRmZpaJE0k6H9cw5dgGx7ENjmMbnFEf26gfIzEzs2x8RmJmZpmMmkQiaaGk\nrZIKki4tUT9e0k1p/YOSWmoU1/GS1kl6QtJmSReVaPOHknZL2ph+fb4WsaWvvU3SY+nr9puwTIl/\nTI/bJknvqFFcJxYdj42Sfi3pk33a1Oy4Sbpe0g5JjxeVNUu6V9Iv0u9vHGDb89M2v5B0fo1i+4qk\nn6c/s9skvWGAbcv+/HOKbbmk7UU/t7MH2Lbs33ROsd1UFNc2SRsH2Dbv41byc6Nuv3MRMeK/SGYK\nfhI4ATgKeBSY3afN3wDfTB8vAW6qUWxvBt6RPn498O8lYvtD4PY6HbttwOQy9WcDPyCZ4v9dwIN1\n+vm+QHLde12OG/AHwDuAx4vK/gG4NH18KdBRYrtm4Kn0+xvTx2+sQWxnAmPTxx2lYqvm559TbMuB\ni6v4mZf9m84jtj71/xv4fJ2OW8nPjXr9zo2WM5KD68dHxF6gd633Yu3ADenjW4HTJZVaqXFIRcTz\nEfFw+vg3wBMcWla4EbQD34nEz4A3SHpzjWM4HXgyIga66TR3EfEjkqUQihX/Tt0AfLDEpu8H7o2I\nnoh4CbgXWJh3bBFxT0TsT5/+jGSRuJob4LhVo5q/6dxiSz8b/hz4l6F8zWqV+dyoy+/caEkkpdaP\n7/thfdj68UDv+vE1k3anvR14sET1uyU9KukHkk6qYVgB3CNpg5Lli/uq5tjmbQkD/0HX67gBvCki\nnofkDx+YWqLNcDh+HyM5qyyl0s8/Lxem3W7XD9A9U+/j9j7gVxHxiwHqa3bc+nxu1OV3brQkkizr\nx9eEpEnAGpJFvH7dp/phkm6btwFfB75Xq7iA90TEO4CzgGWS/qBPfb2P21HAIuCWEtX1PG7Vqvfx\n+yzJ4nE3DtCk0s8/D9cAM4G5wPMkXUh91fW4AedR/mykJsetwufGgJuVKMt07EZLIjmS9ePR4evH\n507SOJJfhhsj4l/71kfEryNiT/r4TmCcpMm1iC0inku/7wBuI+lSKFbNsc3TWcDDEfGrvhX1PG6p\nX/V286Xfd5RoU7fjlw6yfgD4UKSd531V8fMfchHxq4g4EBGvAv80wGvW87iNBf4UuGmgNrU4bgN8\nbtTld260JJKD68en/8EuIVkvvljv+vFQtH583oGlfa3XAU9ExBUDtDm2d7xG0qkkP7cXaxDbREmv\n731MMkD7eJ9mncBfKvEuYHfvqXWNDPifYb2OW5Hi36nzgbUl2twNnCnpjWkXzplpWa4kLQQuARZF\nxO8GaFPNzz+P2IrH2P5kgNes5m86L38E/DwiuktV1uK4lfncqM/vXF5XFQy3L5Kri/6d5EqPz6Zl\nl5P8IQFMIOkeKQAPASfUKK73kpxWbgI2pl9nA38F/FXa5kJgM8mVKT8D5tcothPS13w0ff3e41Yc\nm4AV6XF9DGir4c/0dSSJoamorC7HjSSZPQ/sI/mPbynJGNt9wC/S781p2zbg2qJtP5b+3hWAj9Yo\ntgJJP3nv71zvFYtvAe4s9/OvQWyr0t+lTSQfjG/uG1v6vN/fdN6xpeXf7v0dK2pb6+M20OdGXX7n\nfGe7mZllMlq6tszMLCdOJGZmlokTiZmZZeJEYmZmmTiRmJlZJmPrHYBZo5J0gOQy1V6rI+LL9YrH\nrF58+a/ZIEnaExGThnifY+PQZIpmDcFdW2ZDLF2L4u8lPZyuSfF7afnEdBLC9ZIekdSelv83SbdI\n+j7JRH+vkXR1us7E7ZLulHSOpNMl3Vb0OmdI6jeljlmtOZGYDd5rdfjiWucW1e2KZNK+a4CL07LP\nkky9Mw9YAHwlnUID4N3A+RFxGsk8Ti3AycB/T+sA7gf+q6Qp6fOPAt/K6b2ZVc1jJGaD958RMXeA\nut4zhQ0kiQGSOY0WSepNLBOA6enjeyOid5LQ9wK3RDJp4QuS1gFEREhaBXxY0rdIEsxfDt3bMRsc\nJxKzfLySfj/Aob8zAYsjYmtxQ0nvBH5bXFRmv98Cvg+8TJJsPJ5ideeuLbPauRv4eNGMxG8foN3/\nAxanYyWhyJyBAAAAs0lEQVRvIlkyGDg4PflzwOdIJg80qzufkZgN3mslbSx6fldEXFqm/ReArwGb\n0mSyjWQ9kL7WkCwf/DjJ7LYPkqzY2etGYEpEbMkQu9mQ8eW/ZsOQpEkRsUfSMSTLGrwnIl5I674B\nPBIR19U1SLOUz0jMhqfbJb0BOAr4QlES2UAynvK39QzOrJjPSMzMLBMPtpuZWSZOJGZmlokTiZmZ\nZeJEYmZmmTiRmJlZJk4kZmaWyf8HLfwpreLtKDsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x115921668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.plot(W1, output_fix, 'o', mfc='w')\n",
    "ax.plot(W1, output_mc, 's', mfc='w')\n",
    "\n",
    "ax.set_xlabel('Energy')\n",
    "ax.set_ylabel('Intensity')"
   ]
  }
 ],
 "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}