Skip to content

Commit

Permalink
1D maxwellian distribution
Browse files Browse the repository at this point in the history
created distribution.py
created test_distribution.py
change `part` argument to `particle` and minor docstrings
Added is_electron


Update formula docstring
  • Loading branch information
Antoine T committed Oct 6, 2017
1 parent 1b1db10 commit a6fbf18
Show file tree
Hide file tree
Showing 4 changed files with 392 additions and 0 deletions.
252 changes: 252 additions & 0 deletions plasmapy/examples/Distribution_exemple.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,252 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-04T10:38:17.117221Z",
"start_time": "2017-10-04T10:38:16.455734Z"
}
},
"outputs": [],
"source": [
"import importlib\n",
"import sys\n",
"sys.path.insert(0, '../../../PlasmaPy')\n",
"\n",
"import numpy as np\n",
"from astropy import units as u\n",
"import plasmapy\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from plasmapy.constants import (m_p, m_e, c, mu0, k_B, e, eps0, pi, e)\n",
"from plasmapy.physics.distribution import Maxwellian_1D\n",
"\n",
"from imp import reload\n",
"reload(plasmapy.physics.distribution)\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-04T10:38:24.834239Z",
"start_time": "2017-10-04T10:38:24.818209Z"
}
},
"outputs": [
{
"data": {
"text/latex": [
"$5.9163297 \\times 10^{-7} \\; \\mathrm{\\frac{s}{m}}$"
],
"text/plain": [
"<Quantity 5.916329687405701e-07 s / m>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Maxwellian_1D(v=1*u.m/u.s, T= 30000*u.K, particle='e',V_drift=0*u.m/u.s)\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-04T10:38:34.361995Z",
"start_time": "2017-10-04T10:38:34.187756Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.9999999999999999\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f101b8146a0>]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEQCAYAAACQip4+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFkFJREFUeJzt3W1sZGd5xvHrmhd7s5uEbIgLKdmwEa1AKBACLpQX0RJa\nGkIEVQsoFJAoqVZIFAUVhKB8aj9VqgRUbQGtKG0REBQgaSvEawUUkBKKF5Iob6UoBNgUtE6WwL57\nzpy7H2bGa6/H9tg+Z4bnmf9PirJrnz2+rThXntzPfZ7jiBAAIB2NSRcAANgaghsAEkNwA0BiCG4A\nSAzBDQCJIbgBIDG1Bbftj9o+YvueCu71Ett3rvjrtO0/rKJOAEiN65rjtv1iScclfSwirqzwvhdL\n+oGkyyLiZFX3BYBU1LbijohvSDq68mO2n2L7i7YP2f6m7adt49avlvQFQhvAtBp3j/ugpLdFxHMk\nvVPSB7dxjxsk3VxpVQCQkNa4vpDt8yW9QNKnbQ8+PNv/3B9J+ushf+zhiPiDFfe4VNIzJH2p3moB\n4FfX2IJbvdX9YxHxrHM/ERG3Srp1hHu8VtJtEdGpujgASMXYWiUR8UtJP7T9Gklyz1VbvM3rRJsE\nwJSrcxzwZkm3S3qq7cO2b5T0ekk32r5L0r2SXrWF++2XtE/Sf1VfLQCko7ZxQABAPXhyEgASU8vm\n5CWXXBL79++v49YAkKVDhw49EhFzo1xbS3Dv379fCwsLddwaALJk+0ejXkurBAASQ3ADQGJGCm7b\nF9n+jO0HbN9v+/l1FwYAGG7UHvffSfpiRLza9oyk3TXWBADYwKbBbftxkl4s6U2SFBFLkpbqLQsA\nsJ5RWiVXSFqU9M+2v2f7I7b3nHuR7QO2F2wvLC4uVl4oAKBnlOBuSXq2pA9FxNWSTkh697kXRcTB\niJiPiPm5uZFGEQEA2zBKcB+WdDgivt3//WfUC3IgKQsPHdWhHx3d/ELgV9ymwR0RP5P0E9tP7X/o\npZLuq7UqoAav/vDt+uMP3T7pMoAdG3Wq5G2SPtGfKHlQ0p/WVxJQPQ5TQ05GCu6IuFPSfM21ALU5\nsdQ9++szhfbMjvMdIkC1eHISU+Ho8bMTrMdOFxOsBNg5ghtT4VRnxYp7ieBG2ghuTIWVwX1qRdsE\nSBHBjalw6pweN5AyghtT4fSKFfdJVtxIHMGNqbCyVXKcFTcSR3BjKqxslZwpyglWAuwcwY2pcLo4\nG9xLBDcSR3BjKqxccS8V9LiRNoIbU2Hl5iStEqSO4MZUWLk5SasEqSO4MRVOd0rtnmmqYWmpS3Aj\nbZy0g6lQdEu1mw2VEbRKkDyCG1OhU4baTSuiQasEySO4MRWKbqlWoyG3rTNMlSBxBDemQtENtZpW\nM0yrBMkjuDEVeq2ShiymSpA+ghtTodcqsZoNVtxIH8GNqdDphlrNhmaaZsWN5DHHjalQlKXaTWu2\n1WRzEskjuDEVim6o2bBaTavo8sZ3pI3gxlQoylLtRkOtZkOdkuBG2ghuTIXBOGC7YRU88o7EEdyY\nCp2ytzlJqwQ5ILgxFYpuqXbDajcb6pSsuJG2kcYBbT8k6ZikrqQiIubrLAqo2nKrpNlQh1YJEreV\nOe6XRMQjtVUC1KhTlr1WSYNWCdJHqwRToeiG2g33pkoIbiRu1OAOSV+2fcj2gWEX2D5ge8H2wuLi\nYnUVAhUour0Vd7tpFfS4kbhRg/tFEfFsSS+X9FbbLz73gog4GBHzETE/NzdXaZHAThVlqNWwWo0G\nrRIkb6TgjoiH+38/Iuk2Sc+tsyigakXZ35xsmc1JJG/T4La9x/YFg19Lepmke+ouDKhSp/8ihXaj\noYInJ5G4UaZKniDpNtuD6z8ZEV+stSqgYkW39+qyVtPqlqGyDDUannRZwLZsGtwR8aCkq8ZQC1Cb\nohxsTvb+J7NTlpptNCdcFbA9jAMiexGhzmAcsL/KZoMSKSO4kb1uv6fdO6uk9yNPcCNlBDeyVywH\nt9Vu9lbcnFeClBHcyN5g/K/dONvjZsWNlBHcyN6gVdJc0eNmlhspI7iRvcHZJO3+6YCSmOVG0ghu\nZG9wNsngRQqSeAsOkkZwI3uDfvbgrBJJWiK4kTCCG9lb3pzsnw4osTmJtBHcyN7qccBBj5sVN9JF\ncCN7gxV3q3G2x83LFJAyghvZK4ZNlRDcSBjBjeytmipp8OQk0kdwI3srp0pYcSMHBDeyt7w52TBz\n3MgCwY3sLW9OrjiPmzlupIzgRvZWbU42aJUgfQQ3sre8OdloqNlvlXQ5qwQJI7iRvVWHTDFVggwQ\n3Mje6kOmaJUgfQQ3stfprp0q4TxupIzgRvaW57hXbk7S40bCCG5kr1uuPauEOW6kjOBG9lZuTp59\ndRkrbqRr5OC23bT9Pdufq7MgoGorNydtq9kw44BI2lZW3DdJur+uQoC6rNycHPydcUCkbKTgtn2Z\npFdI+ki95QDVO/vkZGP574wDImWjrrg/IOldktZdptg+YHvB9sLi4mIlxQFVKMpSttQcrLibZnMS\nSds0uG1fL+lIRBza6LqIOBgR8xExPzc3V1mBwE51urE8Bij1pks69LiRsFFW3C+U9ErbD0n6lKRr\nbH+81qqAChXdcnkMUOpNl7DiRso2De6IeE9EXBYR+yXdIOmrEfGG2isDKlKUsdwmkQatElbcSBdz\n3MheUZbLG5MSrRKkr7WViyPi65K+XkslQE2KbiyPAkq9ccAu44BIGCtuZK/TjdUr7maDJyeRNIIb\n2StKNieRF4Ib2RvWKuF0QKSM4Eb2Ot1ySKuEFTfSRXAje0UZQ1olrLiRLoIb2et0S7VWPDnZZBwQ\niSO4kb1ze9xtxgGROIIb2eue0yrhyUmkjuBG9jrnPjnJ5iQSR3Aje8NaJYwDImUEN7LX6ZZqnbPi\nplWClBHcyF5RhtrnjAPSKkHKCG5kr1gzDkirBGkjuJG9TvecqZIGm5NIG8GN7BVluerVZe2m1WXF\njYQR3Mhe0Q01V81xszmJtBHcyF5RhtrnjAN2eHISCSO4kb1iyDhghGiXIFkEN7LXGfLIuyQ2KJEs\nghvZK7qrNycHT1EyEohUEdzIWlmGytCacUBJvL4MySK4kbXBJuTKQ6bay60SVtxIE8GNrA3G/la9\nc7If4mxOIlUEN7K2HNzNtT1uNieRKoIbWRu0SlYd69oPcTYnkapNg9v2Ltv/bfsu2/fa/qtxFAZU\nYdAOGTYOyOYkUtUa4Zozkq6JiOO225K+ZfsLEXFHzbUBOzZohwwbB2RzEqnaNLgjIiQd7/+23f+L\nn3gk4WyPe8g4II+9I1Ej9bhtN23fKemIpK9ExLeHXHPA9oLthcXFxarrBLZlEM6rH3lnxY20jRTc\nEdGNiGdJukzSc21fOeSagxExHxHzc3NzVdcJbMsgnNtDNicZB0SqtjRVEhGPSfqapGvrKQeo1kbj\ngGxOIlWjTJXM2b6o/+vzJP2+pAfqLgyowvI4YHPtAzgdVtxI1ChTJZdK+lfbTfWC/paI+Fy9ZQHV\nKJZbJWsfeWfFjVSNMlVyt6Srx1ALULnB5mRzRY+7yTggEseTk8ja8oq7OezJSVbcSBPBjawNHQdc\n3pxkxY00EdzIWmfI6YCDFTeHTCFVBDeyNgjnmdbaB3CY40aqCG5kbeh53A3GAZE2ghtZWz5kasgb\ncBgHRKoIbmStGHKsa5PNSSSO4EbWBqvq1qoHcAatElbcSBPBjax1hsxxMw6I1BHcyNqwOe4mh0wh\ncQQ3sjZsjtu22k3zzkkki+BG1s4+8r76R73VaBDcSBbBjax1uqUaXn3IlNSbMuHJSaSK4EbWOmW5\nqr890GqYzUkki+BG1opurHpt2UCr2eB0QCSL4EbWiu7wFXe7Yc7jRrIIbmStU8aqGe6BVrPBOCCS\nRXAja0W3XPXU5ECraQ6ZQrIIbmSt6Maqc0oG2o2GurRKkCiCG1nrtUqGr7jZnESqCG5krVOUw3vc\nbE4iYQQ3slaU6/W4GQdEughuZK3TXWeqhBU3EkZwI2vFOk9OthkHRMI2DW7b+2x/zfZ9tu+1fdM4\nCgOq0OnGqpMBB1qcDoiEtUa4ppD0joj4ru0LJB2y/ZWIuK/m2oAdK7qlds+s/TFvNRqcVYJkbbri\njoifRsR3+78+Jul+SU+quzCgCkW5zhw344BI2JZ63Lb3S7pa0reHfO6A7QXbC4uLi9VUB+xQr1Wy\n9se8yemASNjIwW37fEmflfT2iPjluZ+PiIMRMR8R83Nzc1XWCGxb0R0+x91uNnhZMJI1UnDbbqsX\n2p+IiFvrLQmoTqdbDn9ykhU3EjbKVIkl/ZOk+yPiffWXBFSns85ZJa1mgzluJGuUFfcLJb1R0jW2\n7+z/dV3NdQGVKMpS7SE9bjYnkbJNxwEj4luS1i5ZgASsdzog44BIGU9OImvr9bhZcSNlBDeyVpTD\nn5xkHBApI7iRtV6rZL3TAUMRhDfSQ3AjWxGhTrnOHHd/Fc55JUgRwY1sdctQhNY9j1sS7RIkieBG\ntgZz2jOt4ZuTknh6EkkiuJGtpaIXysOCe7BhyYobKSK4ka0z3a6kdYJ70CphxY0EEdzI1nKrZJ1X\nl0msuJEmghvZ2rBVwuYkEkZwI1vLwd1srvkcm5NIGcGNbG28OcmKG+kiuJGtpQ03J/srbt70jgQR\n3MjWmeVWyfpz3Dw5iRQR3MjW2Qdw1k6VzLZ6fe9BOwVICcGNbG20OTnbb5+cKbpjrQmoAsGNbG20\nOTlYcZ/psOJGeghuZGujzcmZ5RU3wY30ENzI1sYrblolSBfBjWwtbTBVMtturLoGSAnBjWxtNA64\n3OMmuJEgghvZ2ug8blolSBnBjWyN1ONmqgQJIriRraVuV82G1RzylvdWs6Fmw7RKkKRNg9v2R20f\nsX3POAoCqrJUlEP72wOzrQatEiRplBX3v0i6tuY6gMotFeXQNslAL7hZcSM9mwZ3RHxD0tEx1AJU\naqlbqr3hirtJjxtJoseNbJ3ulNo9s/ackoHZdkNLHOuKBFUW3LYP2F6wvbC4uFjVbYFtO7XU1a42\nPW7kp7LgjoiDETEfEfNzc3NV3RbYtlOdrs5rr7/inmk1aJUgSbRKkK1Tna52bRDcs60mm5NI0ijj\ngDdLul3SU20ftn1j/WUBO3em09V5G/W4aZUgUa3NLoiI142jEKBqpzpdXdraOLiPnynGWBFQDVol\nyNapTVfcjAMiTQQ3snVqqdywx72r3dCpDq0SpIfgRrbObDJVsnu2pZNLBDfSQ3AjW72pkvV/xPfM\nNHVyiR430kNwI0udbqmijI1X3DO9FXdZxhgrA3aO4EaWBr3rjTYnz5/tDVWdpM+NxBDcyNLpfu96\no83J3bO9z51kJBCJIbiRpcGm40aHTO2Z6a24T7BBicQQ3MjSsdO9VfQFu9rrXjMI9ROsuJEYghtZ\nOnamI+lsH3uYwecIbqSG4EaWji+vuNcP7t2DzUlaJUgMwY0sDc4g2WjFvWfQKmGWG4khuJGlY1tY\ncdMqQWoIbmRpecW9QXAPVuODkAdSQXAjS8dOF5ppNjS7wbGuF+5qqdmwHjvZGWNlwM4R3MjSL051\ndOF5Gx83b1t7d7d19OTSmKoCqkFwI0uPHj+jx++Z3fS6vbtn9PMTBDfSQnAjS0dPLOnx589set3e\nPTM6SnAjMQQ3svToiSVdvGfz4L5494x+TqsEiSG4kaVHjp/RJeeP0CphxY0EEdzIzqmlro6dLnTJ\nCK2SJ1w4q0dPLPG2dySF4EZ2Dv/8pCRp38W7N712397dipAe/vmpussCKkNwIzs/PtoL7stHCO7L\nH7971Z8BUkBwIzs/fOSEpBGDu3/NQ/0/A6SA4EZ27j78C/3643bp8SNsTv7aBbOau2BWdx/+xRgq\nA6oxUnDbvtb2/9j+ge13110UsF3dMnTHg4/q6sv3jnS9bT3n8r2648FHeWkwkrFpcNtuSvpHSS+X\n9HRJr7P99LoLA7bjs4cO68ixM7ruGZeO/Geue+al+r9fnNa/3flwjZUB1dn4MIee50r6QUQ8KEm2\nPyXpVZLuq7qY6//+mzrdKVd9LGLtKmjoumjIB4ddN+r9hlymGHLl0OtGXLiNpZYR7zfsytHvt4Pv\no8J/vmWETi519Zwn79W1Vz5x2F2GevmVT9RV+y7SX9xyl/7mCw9oz2zv8KlxGM9Xwbjs3T2jW97y\n/Nq/zijB/SRJP1nx+8OSnnfuRbYPSDogSZdffvm2ivmNufPV6Q75N3TIT/ewH3h77UeHX1ft/YbX\nN+TPjvx1d3C/EQscSy1D7zdaVG336+7be55e+1v7thS87WZDn/yz5+lT3/mJHvjpL3Wq0x35P747\nMew/vkjbhRu847RKowT3SCLioKSDkjQ/P7+tn8gP3HB1VeUAW7JntqUbX3TFpMsARjLK5uTDkvat\n+P1l/Y8BACZglOD+jqTftH2F7RlJN0j6j3rLAgCsZ9NWSUQUtv9c0pckNSV9NCLurb0yAMBQI/W4\nI+Lzkj5fcy0AgBHw5CQAJIbgBoDEENwAkBiCGwAS42GPHO/4pvaipB9VfuN6XSLpkUkXMWZ8z9OB\n7zkNT46IuVEurCW4U2R7ISLmJ13HOPE9Twe+5/zQKgGAxBDcAJAYgvusg5MuYAL4nqcD33Nm6HED\nQGJYcQNAYghuAEgMwT2E7XfYDtuXTLqWutn+W9sP2L7b9m22L5p0TXWYthde295n+2u277N9r+2b\nJl3TuNhu2v6e7c9Nupa6ENznsL1P0ssk/XjStYzJVyRdGRHPlPR9Se+ZcD2Vm9IXXheS3hERT5f0\n25LeOgXf88BNku6fdBF1IrjXer+kd2mdd9bmJiK+HBFF/7d3qPeGo9wsv/A6IpYkDV54na2I+GlE\nfLf/62PqBdmTJltV/WxfJukVkj4y6VrqRHCvYPtVkh6OiLsmXcuEvFnSFyZdRA2GvfA6+xAbsL1f\n0tWSvj3ZSsbiA+otvMpJF1Knyl4WnArb/ynpiUM+9V5Jf6lemyQrG33PEfHv/Wveq97/Xn9inLWh\nXrbPl/RZSW+PiF9Oup462b5e0pGIOGT7dyddT52mLrgj4veGfdz2MyRdIeku21KvZfBd28+NiJ+N\nscTKrfc9D9h+k6TrJb008hzsn8oXXttuqxfan4iIWyddzxi8UNIrbV8naZekC21/PCLeMOG6KscD\nOOuw/ZCk+YhI7YSxLbF9raT3SfqdiFicdD11sN1Sb+P1peoF9nck/UnO7051b/Xxr5KORsTbJ13P\nuPVX3O+MiOsnXUsd6HHjHyRdIOkrtu+0/eFJF1S1/ubr4IXX90u6JefQ7nuhpDdKuqb/z/XO/koU\nGWDFDQCJYcUNAIkhuAEgMQQ3ACSG4AaAxBDcALBDtj9q+4jte0a49v0rJn2+b/uxLX89pkoAYGds\nv1jScUkfi4grt/Dn3ibp6oh481a+HituANihiPiGpKMrP2b7Kba/aPuQ7W/aftqQP/o6STdv9etN\n3SPvADAmByW9JSL+1/bzJH1Q0jWDT9p+snrHbHx1qzcmuAGgYv3DvV4g6dP9s48kafacy26Q9JmI\n6G71/gQ3AFSvIemxiHjWBtfcIOmt2705AKBC/SN0f2j7NVLv0C/bVw0+3+9375V0+3buT3ADwA7Z\nvlm9EH6q7cO2b5T0ekk32r5L0r1a/dalGyR9arvHKDMOCACJYcUNAIkhuAEgMQQ3ACSG4AaAxBDc\nAJAYghsAEkNwA0Bi/h99lO882H4frAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f101d8817b8>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# With a vertor (from Numpy)\n",
"start = -5000;\n",
"stop = - start;\n",
"dv = 10000 * u.m/u.s;\n",
"\n",
"v = np.arange(start,stop,dtype='float64')* dv;\n",
"\n",
"#Test the normationation to 1\n",
"integral = (Maxwellian_1D(v,T= 30000*u.K, particle='e')).sum()*dv\n",
"\n",
"print(integral)\n",
"\n",
"plt.plot(v,(Maxwellian_1D(v,T= 30000*u.K, particle='e')))\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-04T10:38:47.008876Z",
"start_time": "2017-10-04T10:38:46.994885Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.0\n"
]
}
],
"source": [
"#test std \n",
"T = 30000*u.K\n",
"std = np.sqrt((Maxwellian_1D(v,T=T, particle='e')*v**2*dv).sum())\n",
"T_theo = (std**2/k_B*m_e).to(u.K)\n",
"\n",
"print(T_theo/T)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-04T10:39:00.264312Z",
"start_time": "2017-10-04T10:39:00.225670Z"
}
},
"outputs": [],
"source": [
"T_e = 30000*u.K;\n",
"V_drift = 10*u.m/u.s;\n",
"\n",
"start = -5000;\n",
"stop = - start;\n",
"dv = 10000 * u.m/u.s;\n",
"\n",
"v_vect = np.arange(start,stop,dtype='float64')* dv;"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-04T10:39:00.477616Z",
"start_time": "2017-10-04T10:39:00.467525Z"
}
},
"outputs": [
{
"data": {
"text/latex": [
"$0 \\; \\mathrm{\\frac{m}{s}}$"
],
"text/plain": [
"<Quantity 0.0 m / s>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v_vect[Maxwellian_1D(v_vect,T=T_e, particle='e',V_drift = 0*u.m/u.s ).argmax()]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-04T08:41:33.617076Z",
"start_time": "2017-10-04T08:41:33.603795Z"
}
},
"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.6.2"
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}
2 changes: 2 additions & 0 deletions plasmapy/physics/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,3 +21,5 @@
from .relativity import Lorentz_factor

from .transport import Coulomb_logarithm

from .distribution import (Maxwellian_1D)
94 changes: 94 additions & 0 deletions plasmapy/physics/distribution.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
"""Functions to deal with distribution : generate, fit, calculate"""

from astropy import units

from astropy.units import (UnitConversionError, UnitsError, quantity_input,
Quantity)

from ..constants import (m_p, m_e, c, mu0, k_B, e, eps0, pi, e)
from ..atomic import (ion_mass, charge_state)
from ..atomic.atomic import _is_electron as is_electron
import numpy as np

from ..utils import _check_quantity, check_relativistic, check_quantity


@units.quantity_input
def Maxwellian_1D(v: units.m/units.s,
T: units.K, particle="e",
V_drift=0*units.m/units.s):
r"""Return the probability at the velocity `v` in m/s
to find a particle `particle` in a plasma of temperature `T`
following the Maxwellian distribution function
Parameters
----------
v: Quantity
The velocity in units convertible to m/s
T: Quantity
The temperature in Kelvin
particle: string, optional
Representation of the particle species(e.g., 'p' for protons, 'D+'
for deuterium, or 'He-4 +1' for $He_4^{+1}$ : singly ionized helium-4),
which defaults to electrons.
Returns
-------
f : Quantity
probability in Velocity^-1, normized so that:
$\int_{- \infty}^{\infty} f(v) dv = 1}
Raises
------
TypeError
The parameter arguments are not Quantities and
cannot be converted into Quantities.
UnitConversionError
If the parameters is not in appropriate units.
ValueError
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In one dimension, the Maxwellian distribution function for a particle of
mass m, velocity v, a drift velocity V and with temperature T is:
.. math::
f = \sqrt{\frac{m}{2 \pi k_B T}} \exp(-\frac{m}{2 k_B T} (v-V)^2)
Example
-------
>>> from plasmapy.physics import Maxwellian_1D
>>> from astropy import units as u
>>> v=1*u.m/u.s
>>> Maxwellian_1D(v=v, T= 30000*u.K, particle='e',V_drift=0*u.m/u.s)
<Quantity 5.916329687405701e-07 s / m>
"""

# pass to Kelvin

if T.unit.is_equivalent(units.K):
T = T.to(units.K, equivalencies=units.temperature_energy())

# Get mass

if is_electron(particle):
m_s = m_e

else:
try:
m_s = ion_mass(particle)
except Exception:
raise ValueError("Unable to find ion mass in Maxwellian_1D")

T = k_B*T

f = np.sqrt((m_s/(2*pi*T)))*np.exp(-m_s/(2*T)*(v-V_drift)**2)

return f.to(units.s/units.m)
44 changes: 44 additions & 0 deletions plasmapy/physics/tests/test_distribution.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
"""Tests for functions that uses Distribution functions."""

import numpy as np
import pytest
from astropy import units as u


from ...constants import (m_p, m_e, c, mu0, k_B, e, eps0, pi, e)
from ...atomic import (ion_mass, charge_state)
from ..distribution import (Maxwellian_1D)

T_e = 30000*u.K
V_drift = 1000000*u.m/u.s

start = -5000
stop = - start
dv = 10000 * u.m/u.s

v_vect = np.arange(start, stop, dtype='float64') * dv


def test_Maxwellian_1D():
"""test the 1D maxwellian distribution function"""

max_index = Maxwellian_1D(v_vect, T=T_e, particle='e', V_drift=0*u.m/u.s
).argmax()
assert np.isclose(v_vect[max_index].value, 0.0)

max_index = Maxwellian_1D(v_vect, T=T_e, particle='e', V_drift=V_drift
).argmax()
assert np.isclose(v_vect[max_index].value, V_drift.value)

# integral of the distribution over v_vect
integral = (Maxwellian_1D(v_vect, T=30000*u.K, particle='e')).sum()*dv
assert np.isclose(integral, 1.0)

std = (Maxwellian_1D(v_vect, T=T_e, particle='e')*v_vect**2*dv).sum()
std = np.sqrt(std)
T_distri = (std**2/k_B*m_e).to(u.K)

assert np.isclose(T_distri.value, T_e.value)

with pytest.raises(ValueError):
Maxwellian_1D(1*u.m/u.s, T=1*u.K, particle='XXX')

0 comments on commit a6fbf18

Please sign in to comment.