{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n1D Maxwellian distribution function\n===================================\n\nWe import the usual and the hero of this notebook, the Maxwellian 1D distribution:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom astropy import units as u\nimport plasmapy\nimport matplotlib.pyplot as plt\nfrom plasmapy.constants import (m_p, m_e, c, mu0, k_B, e, eps0, pi, e)\nfrom plasmapy.physics.distribution import Maxwellian_1D"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Given we'll be plotting:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from astropy.visualization import quantity_support\nquantity_support()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's get the probability density of finding an electron at 1 m/s if we\nhave a plasma at 30 000 K:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Maxwellian_1D(\n    v=1 * u.m / u.s,\n    T=30000 * u.K,\n    particle='e',\n    V_drift=0 * u.m / u.s)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note the units! Integrated over velocities, this will give us a\nprobability. Let's test that for a bunch of particles:\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "With a vector (from Numpy)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "start = -500000\nstop = -start\nv = np.arange(start, stop) * 10 * u.m / u.s\ndv = v[1] - v[0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Test the normalization to 1 (the particle has to be somewhere)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for particle in ['p', 'e']:\n    pdf = Maxwellian_1D(v, T=30000 * u.K, particle=particle)\n    integral = (pdf).sum() * dv\n    print(f\"Integral value for {particle}: {integral}\")\n    plt.plot(v, pdf, label=particle)\nplt.legend()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The standard deviation of this distribution should give us back our\ntemperature:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "T = 30000 * u.K\nstd = np.sqrt((Maxwellian_1D(v, T=T, particle='e') * v**2 * dv).sum())\nT_theo = (std**2 / k_B * m_e).to(u.K)\n\nprint(T_theo / T)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And the center of the distribution is, as can be seen below:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "T_e = 30000 * u.K\nV_drift = 10 * u.km / u.s\n\nstart = -5000\nstop = - start\ndv = 10000 * u.m / u.s\n\nv_vect = np.arange(start, stop, dtype='float64') * dv\n\nprint(v_vect[Maxwellian_1D(v_vect, T=T_e,\n                           particle='e', V_drift=V_drift).argmax()])"
      ]
    }
  ],
  "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.6.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}