<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Particule dans une boîte uni-dimensionnelle
Germain Salvato Vallverdu germain.vallverdu@univ-pau.fr
Step2:  Description du système
On considère une particule de masse $m$, se déplaçant librement sur un segment de longuer $L$. Pour un tel système, la partie énergie potentielle étant nulle (on choisit $\hat{V} = Cste = 0$) l'équation de Schrödinger s'écrit 
Step3:  Vérification de la condition de normalisation 
Step4:  Représentation graphique des solutions
Step5:  Visualisation de la densité de probabilité de présence
Le module au carré de la fonction d'onde représente la densité de probabilité de présence de la particule. C'est la seule information accessible, la trajectoire n'étant pas définie en mécanique quantique du fait qu'on ne peut pas connaîter à la fois la vitesse et la position de la particule. Pour visualiser la densité de probabilité de présence, on fait l'expérience suivante 
Step9:  Sur la figure ci-dessus 
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.constants import h, m_e, elementary_charge
%matplotlib inline
Explanation: Particule dans une boîte uni-dimensionnelle
Germain Salvato Vallverdu germain.vallverdu@univ-pau.fr
End of explanation
def phi(x, p=1, L=1):
    return np.sqrt(2 / L) * np.sin(p * np.pi * x / L)
def epsilon(p, m=m_e, L=1):
     return energy in eV 
    return h**2 * p**2 / (8 * m * L**2) / elementary_charge
Explanation: Description du système
On considère une particule de masse $m$, se déplaçant librement sur un segment de longuer $L$. Pour un tel système, la partie énergie potentielle étant nulle (on choisit $\hat{V} = Cste = 0$) l'équation de Schrödinger s'écrit :
\begin{align}
    \hat{\mathcal{H}} \phi & = \varepsilon\phi \
    -\frac{\hbar^2}{2m}\frac{d^2\phi}{dx^2} & = \varepsilon \phi
\end{align}
Les conditions aux limites ou sur la fonction d'onde sont les suivantes :
\begin{equation}
    \begin{cases}
        \phi(0) & = 0 \
        \phi(L) & = 0 \
        \displaystyle\int_0^L \left\vert\phi(x)\right\vert^2 dx & = 1
    \end{cases}
\end{equation}
Lles fonctions propres et valeurs propres de l'opérateur Hamiltonien dépendent alors d'un nombre quantique $p$ et ont pour expressions :
\begin{equation}
    \begin{cases}\displaystyle
        \displaystyle \phi_p(x) & =\displaystyle \sqrt{\frac{2}{L}} \sin\left(\frac{p\pi x}{L}\right) \
        \displaystyle \varepsilon_p & = \displaystyle\frac{h^2p^2}{8 m L^2}
    \end{cases}
\end{equation}
End of explanation
L = 1
x = np.linspace(0, L, 100)
np.trapz(phi(x, p=2, L=L)**2, x)
Explanation: Vérification de la condition de normalisation :
End of explanation
fig, axes = plt.subplots(
    nrows=3, figsize=(7.75, 10.25), 
    sharex=True, sharey=True, 
    gridspec_kw=dict(wspace=.05, hspace=.0)
)
L = 1
x = np.linspace(0, L, 100)
for i in range(3):
    p = 3 - i
    axes[i].plot(x, phi(x, p, L), label="$\phi_p(x)$")
    axes[i].plot(x, phi(x, p, L)**2, color="C3")    
    axes[i].fill(x, phi(x, p, L)**2, color="C3", label=r"$\left\vert\phi_p(x)\right\vert^2$", alpha=.2)
    axes[i].plot((0, L), (0, 0), color="C7", linewidth=.5)
    axes[i].set_ylabel("p = %d" % p, fontsize=20)
    axes[i].axhline(0, color="C7")
    axes[i].grid(False)
    
# layout
axes[2].set_xlim((0, L))
axes[2].set_xlabel("$x$", fontsize=18)
#axes[2].legend(fontsize=18)
axes[2].set_xticks([0, L/4, L/2, 3*L/4, L])
text = axes[2].set_xticklabels(["0", "L/4", "L/2", "3L/4", "L"], fontsize=14)
axes[0].legend(fontsize=18, ncol=2, bbox_to_anchor=(0.5, 1), loc="lower center")
plt.tight_layout()
plt.savefig("puits_infini.pdf")
fig, axes = plt.subplots(
    nrows=3, figsize=(7.75, 10.25), 
    sharex=True, sharey=True, 
    gridspec_kw=dict(wspace=.05, hspace=.0)
)
L = 1
x = np.linspace(0, L, 1000)
ps = [50, 10, 1]
for i in range(3):
    p = ps[i]
    axes[i].plot(x, phi(x, p, L), label="$\phi_p(x)$")
    axes[i].plot(x, phi(x, p, L)**2, color="C3")    
    axes[i].fill(x, phi(x, p, L)**2, color="C3", label=r"$\left\vert\phi_p(x)\right\vert^2$", alpha=.2)
    axes[i].plot((0, L), (0, 0), color="C7", linewidth=.5)
    axes[i].set_ylabel("p = %d" % p, fontsize=20)
    axes[i].axhline(0, color="C7")
    axes[i].grid(False)
    
# layout
axes[2].set_xlim((0, L))
axes[2].set_xlabel("$x$", fontsize=18)
#axes[2].legend(fontsize=18)
axes[2].set_xticks([0, L/4, L/2, 3*L/4, L])
text = axes[2].set_xticklabels(["0", "L/4", "L/2", "3L/4", "L"], fontsize=14)
axes[0].legend(fontsize=18, ncol=2, bbox_to_anchor=(0.5, 1), loc="lower center")
plt.tight_layout()
plt.savefig("puis_infini_heigh.pdf")
Explanation: Représentation graphique des solutions
End of explanation
fig, axes = plt.subplots(
    nrows=2, ncols=3, figsize=(14, 8), 
    gridspec_kw=dict(wspace=0, hspace=.0, height_ratios=[2, 1])
)
[ax.grid(False) for ax in axes.flatten()]
# nombre de photos
ntry = 500
# dimension de la boite
L = 1
jitter = .5
x = np.linspace(0, L, 100)
for i in range(3):
    p = i + 1
    
    # sample the density
    pos = np.random.uniform(0, L, ntry)
    rho = phi(pos, p, L)**2
    rnd = np.random.uniform(0, 2/L, ntry)
    ix = np.where(rho > rnd)
    pos = pos[ix]
    # plot
    axes[0, i].plot(x, phi(x, p, L)**2, color="C3")
    axes[0, i].hist(pos, bins=20, density=True, alpha=.4, edgecolor="C7")
    axes[1, i].scatter(pos, np.random.normal(0, jitter, ix[0].shape), color="C7", s=30, alpha=.6)
    # layout
    axes[0, i].set_title("$p$ = %d" % p, fontsize=18, y=1.05)
    axes[0, i].yaxis.set_visible(False)
    axes[1, i].set_ylim(-5, 5)
    axes[1, i].yaxis.set_visible(False)
    axes[1, i].set_xlim((0, L))
    axes[1, i].set_xticks([0, L/4, L/2, 3*L/4, L])
    axes[1, i].set_xlabel("$x$", fontsize=18)
    if i > 0:
        axes[0, i].yaxis.set_visible(False)
        axes[1, i].set_xticks([L/4, L/2, 3*L/4, L])
        text = axes[1, i].set_xticklabels(["L/4", "L/2", "3L/4", "L"], fontsize=18)
    else:
        axes[1, i].set_xticks([0, L/4, L/2, 3*L/4, L])
        text = axes[1, i].set_xticklabels(["0", "L/4", "L/2", "3L/4", "L"], fontsize=18)
fig.tight_layout()
fig.savefig("puits_infini_density.pdf")
Explanation: Visualisation de la densité de probabilité de présence
Le module au carré de la fonction d'onde représente la densité de probabilité de présence de la particule. C'est la seule information accessible, la trajectoire n'étant pas définie en mécanique quantique du fait qu'on ne peut pas connaîter à la fois la vitesse et la position de la particule. Pour visualiser la densité de probabilité de présence, on fait l'expérience suivante : on enregistre un grand nombre de fois (ntry fois ici) une photo du système et on superpose l'ensemble de ces photos. Le nuage de points obtenu correspond aux positions accessibles par la particule. Plus les points sont resserrés plus la probabilité de trouver la particule à cet endroit est grande.
Ici on réalise numériquement cette expérience en tirant aléatoirement des points en accord avec la densité de probabilité de présence $\left|\phi_p(x)\right|^2$.
End of explanation
%matplotlib widget
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np
def phi(x, p=1, L=1):
    return np.sqrt(2 / L) * np.sin(p * np.pi * x / L)
plt.ion()
def sample(p=1, ntry=100, L=1):
    n = 0
    pos = list()
    while n < ntry:
        x = np.random.uniform(0, L)
        rho = phi(x, p, L)**2
        if rho > np.random.uniform(0, 2 / L):
            n += 1
            pos.append(x)
    return pos
# initialize
L = 1
jitter = .5
x = np.linspace(0, L, 100)
nbins = 30
ntry = 1
p = 1
pos = sample(p, ntry, L)
# set up figure
output = widgets.Output()
with output:
    fig, axes = plt.subplots(
        nrows=2, ncols=1, figsize=(10, 6), 
        sharex=True,
        gridspec_kw=dict(wspace=0, hspace=.0, height_ratios=[2, 1])
    )
    fig.canvas.toolbar_visible = False
    fig.canvas.header_visible = False
    fig.canvas.footer_visible = False
    fig.patch.set_facecolor("white")
axes[0].grid(False)
axes[0].set_ylim(0, 1.2 * 2 / L)
axes[0].set_title(f"p = {p}", fontsize=18, y=1.05)
axes[0].yaxis.set_visible(False)
axes[1].grid(False)
axes[1].set_xlim(0, L)
axes[1].yaxis.set_visible(False)
axes[1].set_ylim(-5, 5)
axes[1].set_xlabel("x", fontsize=18)
axes[1].set_xticks([0, L/4, L/2, 3*L/4, L])
text = axes[1].set_xticklabels(["0", "L/4", "L/2", "3L/4", "L"], fontsize=18)
density = axes[0].plot(x, phi(x, p, L)**2, color="C3")
hist, bins, patchs = axes[0].hist(pos, bins=nbins, density=True, alpha=.4, edgecolor="C7", range=(0, L))
scatter = axes[1].scatter(pos, np.random.normal(0, jitter, ntry), color="C7", s=30, alpha=.6)
p_slider = widgets.IntSlider(value=1, min=1, max=5, step=1, description="p")
ntry_slider = widgets.IntSlider(value=1, min=1, max=2000, step=10, description="N points")
replot_btn = widgets.Button(description="replot", button_style="", icon="redo")
def update(change):
     callback function to update the plot depending on the triggered control 
    if isinstance(change, widgets.Button):
        update_plot(p=p_slider.value, ntry=ntry_slider.value)
    elif change.owner == p_slider:
        update_plot(p=change.new, ntry=ntry_slider.value)
    elif change.owner == ntry_slider:
        update_plot(p=p_slider.value, ntry=change.new)
    
def update_plot(p=1, ntry=1):
     update the plot according to p and ntry 
    # sample the density
    pos = sample(p, ntry, L)
    
    axes[0].set_title(f"p = {p}", fontsize=18, y=1.05)
    density[0].set_data(x, phi(x, p, L)**2)
    hist, bins = np.histogram(pos, bins=nbins, density=True, range=(0, L))
    [patch.set_height(hist[i]) for i, patch in enumerate(patchs)]
    scatter.set_offsets(np.c_[pos, np.random.normal(0, jitter, ntry)])
    
p_slider.observe(update, "value")
ntry_slider.observe(update, "value")
replot_btn.on_click(update)
controls = widgets.HBox([p_slider, ntry_slider, replot_btn])
widgets.VBox([controls, output])
# set up plot
fig, axes = plt.subplots(
    nrows=2, ncols=1, figsize=(10.25, 7.75), 
    sharex=True,
    gridspec_kw=dict(wspace=0, hspace=.0, height_ratios=[2, 1])
)
fig.canvas.toolbar_visible = False
fig.canvas.header_visible = False
fig.canvas.footer_visible = False
axes[0].grid(False)
axes[1].grid(False)
axes[1].set_xlim(0, L)
axes[1].yaxis.set_visible(False)
axes[1].set_ylim(-5, 5)
axes[1].set_xlabel("x", fontsize=18)
axes[1].set_xticks([0, L/4, L/2, 3*L/4, L])
text = axes[1].set_xticklabels(["0", "L/4", "L/2", "3L/4", "L"], fontsize=18)
# data
L = 1
jitter = .5
x = np.linspace(0, L, 100)
nbins = 30
# initial plot
ntry = 100
p = 1
pos = sample(p, ntry, L)
density = axes[0].plot(x, phi(x, p, L)**2, color="C3")
hist, bins, patchs = axes[0].hist(pos, bins=nbins, density=True, alpha=.4, edgecolor="C7", range=(0, L))
scatter = axes[1].scatter(pos, np.random.normal(0, jitter, ntry), color="C7", s=10)
@widgets.interact(p=(1, 5, 1), ntry=(1, 1010, 10))
def update(p=1, ntry=100):
    Remove old lines from plot and plot new one
    
    # sample the density
    pos = sample(p, ntry, L)
    
    axes[0].set_title(f"p = {p}", fontsize=18, y=1.05)
    density[0].set_data(x, phi(x, p, L)**2)
    hist, bins = np.histogram(pos, bins=nbins, density=True, range=(0, L))
    [patch.set_height(hist[i]) for i, patch in enumerate(patchs)]
    scatter.set_offsets(np.c_[pos, np.random.normal(0, jitter, ntry)])
Explanation: Sur la figure ci-dessus :
La courbe rouge est la densité de probabilité de présence exacte : $\left|\phi_p(x)\right|^2$
L'histogramme est construit à partir de l'ensemble des positions de la particules. Les positions ont été choisies de sortes que cette histogramme correspondent à la densité de probabilité exacte.
Les points gris correspondent à la superposition des photos du système.
On observe que :
le nombre de points est important au voisinage des maxima de la densité de probabilité de présence
le nombre de points est nul au voisinages des points nodaux (points où la fonction d'onde est nulle)
End of explanation
<END_TASK>