Hopefully this is a quick fix, but I've been working on this for 3 days to no avail, so perhaps some external perspective will help.
Basically, I'm building a MatLab module to calculate mean particle $\frac{dE}{dx}$ from the Bethe-Bloch equation. To verify that I'd implemented it correctly, I used the particle Data Group chapter on "Particle Interaction with matter" (see link at end for a pdf thereof). Unfortunately, as you can see in the following plot, my implementation does not match the data very well at all.
The PDG data was digitized from the widely available plot seen below:
To generate the blue data series in the first plot, I merely calculated the following quantities for a range of (total) proton energies $E_n$ in MeV:
$\gamma_n = \frac{E_n}{E_0}$$$\gamma_n = \frac{E_n}{E_0}$$
where $E_0$ is the rest energy of the incident particle (proton) in MeV.
$\beta_n = \sqrt{1-(\frac{1}{\gamma_n})^2}$$$\beta_n = \sqrt{1-\bigg(\frac{1}{\gamma_n}\bigg)^2}$$
$p = \sqrt{E_n^2 - (\frac{E_0}{c^2})^2} = \frac{E_0}{c^2} \beta \gamma$$$p = \sqrt{E_n^2 - \bigg(\frac{E_0}{c^2}\bigg)^2} = \frac{E_0}{c^2} \beta \gamma$$
$T_{max}=\frac{2 m_e c^2 \beta^2 \gamma^2}{1+\frac{2 \gamma m_e}{M}+(\frac{m_e}{M})^2}$$$T_{max}=\cfrac{2 m_e c^2 \beta^2 \gamma^2}{1+\cfrac{2 \gamma m_e}{M}+\bigg(\cfrac{m_e}{M}\bigg)^2}$$
Where $p$ is the particle momentum, $M$ is the rest mass of the incident proton in g, $m_e$ is the rest mass of an electron in g,and $c$ is, of course, the speed of light in m/s. I calculated the density effect, $\delta$ with the Sternheimer parameterization where:
$x = \log{\beta \gamma}$
$$x = \log{\beta \gamma}$$
$$
\delta(\beta \gamma) =
\begin{cases}
2(\ln 10)x -\bar{C} & \text{if }x\geq x_1; \\
2(\ln 10)x -\bar{C}+a(x_1-x)^k & \text{if }x_0\leq x < x_1; \\
0 & \text{if }x<x_0\ \text{(nonconductors)}; \\
\delta_010^{2(x-x_0)} & \text{if }x<x_0\ \text{(conductors)}; \\
\end{cases}
$$
Finally, the Bethe-Bloch equation I used (straight from the PDG pdf) was as follows:
$-\frac{dE}{dx} = K z^2 \frac{Z_{target}}{A_{target}} \frac{1}{\beta^2} [\frac{1}{2} \ln(\frac{2 m_e c^2 \beta^2 \gamma^2 Tmax}{I^2}-\beta^2 - \delta/2]$$$-\frac{dE}{dx} = K z^2 \frac{Z_{target}}{A_{target}} \frac{1}{\beta^2} \bigg[\frac{1}{2} \ln\frac{2 m_e c^2 \beta^2 \gamma^2 Tmax}{I^2}-\beta^2 - \delta/2\bigg]$$
For lead, the following material parameters were used: $Z_{target} = 82$ amu, $A_{target}=207.2$ g/mol, $I= 823$ eV. For the calculation of $\delta$, the following parameters were used (see link below for source): $C=6.2018$, $X_0=0.3776$, $X_1=3.8073$, $a=0.09359$, $k=3.1608$,and $\delta_0=0.14$. $K$ was taken to be constant at a value of 0.307075 $\frac{MeV mol}{g}$.
Unfortunately I still get the error evident in the first plot where the Fermi Plateau doesn't seem to be in evidence. Does anybody see anything wrong in the steps I took to calculate $\frac{dE}{dx}$?
For the record, this is NOT a problem with programming/MatLab, my trouble lies with the computation of $\frac{dE}{dx}$. The issue is simply that the Bethe Bloch equation, when computed, does not appear to reflect the behavior seen in the PDG plot from literature.
PDG Link: http://pdg.lbl.gov/2009/reviews/rpp2009-rev-passage-particles-matter.pdf
Sternheimer Parameterization Link: https://www.sciencedirect.com/science/article/pii/0092640X84900020
