# Plotting Maxwellian

**URL:** <https://discourse.julialang.org/t/plotting-maxwellian/76385>\
**Category:** New to Julia\
**Tags:** plots\
**Created:** [February 14, 2022, 3:32am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385 "2022-02-14T03:32:18Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![CentraCep](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/centracep/32/32454_2.png) [@CentraCep](https://discourse.julialang.org/u/CentraCep)\
**Post date:** [February 14, 2022, 3:32am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/1 "2022-02-14T03:32:18Z")

</div>

I have a list of positive and negative values and a single temperature. I am trying to plot the [Maxwell-Boltzmann Distribution][1] using the equation for particles moving in only one direction.

```julia
m_e = 9.11E-28 # electron mass [g]
k = 1.38E-16 # boltzmann constant [erg*K^-1]
v = range(1e10, -1e10, step=-1e8) # velocity [cm/s]

T_M = 1e6 # temperature of Maxwellian [K]

function Maxwellian(v_Max, T_Max)
    normal = (m_e/(2*pi*k*T_Max))^1.5
    exp_term = exp(-((m_e).*v_Max.*v_Max)/(3*k*T_Max))
    return normal*exp_term
end

# Initially comparing chosen distribution f_s to Maxwellian F_s
plot(v, Maxwellian.(v, T_M), label= L"F_s" * " (Maxwellian)")
xlabel!("velocity (cm/s)")
ylabel!("probability density")

```

However, when, plotting this, my whole function is 0:

[![enter image description here][2]][2]

I tested out if I wrote my function correctly by replacing `return normal*exp_term` with `return exp_term` (i.e. ignoring any normalization constants) and this seems to produce the distinct of the bell curve:

[![enter image description here][3]][3]

Yet, without the normalization constant, this will not preserve the area under the curve. I was wondering what may I be doing incorrectly with setting up my Maxwellian function and the constant in front of the exponential.  
[1]: [Maxwell–Boltzmann distribution - Wikipedia](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution#Derivation_and_related_distributions)  
[2]: [https://i.stack.imgur.com/tYogq.png](https://i.stack.imgur.com/tYogq.png)  
[3]: [https://i.stack.imgur.com/wyKeE.png](https://i.stack.imgur.com/wyKeE.png)

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 14, 2022, 5:38am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/2 "2022-02-14T05:38:16Z")

</div>

The data extrema are:

```julia
(2.9269998814999053e-123, 1.0769341115495682e-27)

```

Did you try manually setting the `plot()` keyword `ylim` accordingly?

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [February 14, 2022, 6:16am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/3 "2022-02-14T06:16:44Z")

</div>

Maybe there’s a choice of units that gives more reasonable size numbers? If so, `Unitful` and `UnitfulRecipes` could be of use.

---

<div class="post-metadata">

**Author:** ![digital\_carver](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/digital_carver/32/33818_2.png) [@digital\_carver](https://discourse.julialang.org/u/digital_carver)\
**Post date:** [February 14, 2022, 6:30am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/4 "2022-02-14T06:30:48Z")

</div>

Settings `ylims` does show the correct plot, as [mentioned in my StackOverflow response](https://stackoverflow.com/questions/71106799/plotting-maxwellian-distribution-in-julia) to this.

Since we have the topic already, I might as well ask this here (I was planning to make a separate post): why are the default `ylims` of Plots.jl off by several orders of magnitude? Plotting this directly with `julia> GR.plot(v, Maxwellian.(v, T_M))` automatically uses appropriate `ylims`, that lead to a visible plot.

![Screenshot_20220214_115901](https://global.discourse-cdn.com/julialang/original/3X/4/d/4d8e8029c3646c154a84d7138a4a40ac5050438e.png)

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [February 14, 2022, 7:04am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/5 "2022-02-14T07:04:00Z")

</div>

I think Plots has a lower threshold where it assumes whatever you’re plotting is “supposed” to be zeros. Note that `eps(1.0)` is about `1e-16`, which could be the (somewhat flawed) reason for the specific choice. Either way, I would probably change the units to give a more intuitive scale.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [February 14, 2022, 7:37am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/6 "2022-02-14T07:37:16Z")

</div>

Linking a related [Plots.jl issue](https://github.com/JuliaPlots/Plots.jl/issues/3830#issuecomment-923051925).

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [February 14, 2022, 9:18am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/7 "2022-02-14T09:18:41Z")

</div>

```julia
using Plots, Unitful, UnitfulRecipes

m_e = 9.11E-28u"g"
k = 1.38E-16u"erg/K" # boltzmann constant [erg*K^-1]
v = range(1e10, -1e10, step=-1e8) .* 1u"cm/s" # velocity [cm/s]

T_M = 1e6u"K" # temperature of Maxwellian [K]

function Maxwellian(v_Max, T_Max)
    normal = (m_e/(2*pi*k*T_Max))^1.5
    exp_term = exp(-((m_e).*v_Max.*v_Max)/(3*k*T_Max))
    return normal*exp_term |> upreferred
end

```

results in a maximum at `1.077e-21 s^3 m^-3` which is not much more reasonable size number than `1e-26s^3 cm^-3`. But in any case `Unitful` is your friend in doing any physical calculations.

P.S. s3/m3 appears wrong units

> for particles moving in only one direction"

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [February 14, 2022, 9:28am UTC](https://discourse.julialang.org/t/plotting-maxwellian/76385/8 "2022-02-14T09:28:58Z")

</div>

Good thing about SI units: `1e-21 s^3/m^3 = 1000 μs^3/m^3`