{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nCold Magnetized Plasma Waves Tensor Elements (S, D, P in Stix's notation)\n=========================================================================\n\nThis example shows how to calculate the values of the cold plasma tensor\nelements for various electromagnetic wave frequencies.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# First, import some basics (and `PlasmaPy`!)\nimport numpy as np\nimport plasmapy as pp\nimport matplotlib.pyplot as plt\nfrom astropy import units as u"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's define some parameters, such as the magnetic field magnitude,\nthe plasma species and densities and the frequency band of interest\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "B = 2 * u.T\nspecies = ['e', 'D+']\nn = [1e18 * u.m ** -3, 1e18 * u.m ** -3]\n\nf = np.logspace(start=6, stop=11.3, num=3001)  # 1 MHz to 200 GHz\nomega_RF = f * (2 * np.pi) * (u.rad / u.s)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "help(pp.physics.cold_plasma_permittivity_SDP)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "S, D, P = pp.physics.cold_plasma_permittivity_SDP(B, species, n, omega_RF)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Filter positive and negative values, for display purposes only.\nStill for display purposes, replace 0 by NaN to NOT plot 0 values\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "S_pos = S * (S > 0)\nD_pos = D * (D > 0)\nP_pos = P * (P > 0)\nS_neg = S * (S < 0)\nD_neg = D * (D < 0)\nP_neg = P * (P < 0)\nS_pos[S_pos == 0] = np.NaN\nD_pos[D_pos == 0] = np.NaN\nP_pos[P_pos == 0] = np.NaN\nS_neg[S_neg == 0] = np.NaN\nD_neg[D_neg == 0] = np.NaN\nP_neg[P_neg == 0] = np.NaN"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(12, 6))\nplt.semilogx(f, abs(S_pos),\n             f, abs(D_pos),\n             f, abs(P_pos), lw=2)\nplt.semilogx(f, abs(S_neg), '#1f77b4',\n             f, abs(D_neg), '#ff7f0e',\n             f, abs(P_neg), '#2ca02c', lw=2, ls='--')\nplt.yscale('log')\nplt.grid(True, which='major')\nplt.grid(True, which='minor')\nplt.ylim(1e-4, 1e8)\nplt.xlim(1e6, 200e9)\nplt.legend(('S > 0', 'D > 0', 'P > 0', 'S < 0', 'D < 0', 'P < 0'),\n           fontsize=16, ncol=2)\nplt.xlabel('RF Frequency [Hz]', size=16)\nplt.ylabel('Absolute value', size=16)\nplt.tick_params(labelsize=14)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Cold Plasma tensor elements in the rotating basis\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "L, R, P = pp.physics.cold_plasma_permittivity_LRP(B, species, n, omega_RF)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "L_pos = L * (L > 0)\nR_pos = R * (R > 0)\nL_neg = L * (L < 0)\nR_neg = R * (R < 0)\nL_pos[L_pos == 0] = np.NaN\nR_pos[R_pos == 0] = np.NaN\nL_neg[L_neg == 0] = np.NaN\nR_neg[R_neg == 0] = np.NaN\n\nplt.figure(figsize=(12, 6))\nplt.semilogx(f, abs(L_pos),\n             f, abs(R_pos),\n             f, abs(P_pos), lw=2)\nplt.semilogx(f, abs(L_neg), '#1f77b4',\n             f, abs(R_neg), '#ff7f0e',\n             f, abs(P_neg), '#2ca02c', lw=2, ls='--')\nplt.yscale('log')\nplt.grid(True, which='major')\nplt.grid(True, which='minor')\nplt.xlim(1e6, 200e9)\nplt.legend(('L > 0', 'R > 0', 'P > 0', 'L < 0', 'R < 0', 'P < 0'),\n           fontsize=16, ncol=2)\nplt.xlabel('RF Frequency [Hz]', size=16)\nplt.ylabel('Absolute value', size=16)\nplt.tick_params(labelsize=14)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Checks if the values obtained are coherent. They should satisfy\nS = (R+L)/2 and D = (R-L)/2\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "try:\n    np.testing.assert_allclose(S, (R + L) / 2)\n    np.testing.assert_allclose(D, (R - L) / 2)\nexcept AssertionError as e:\n    print(e)\n# Checks for R=S+D and L=S-D\ntry:\n    np.testing.assert_allclose(R, S + D)\n    np.testing.assert_allclose(L, S - D)\nexcept AssertionError as e:\n    print(e)"
      ]
    }
  ],
  "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
}