529.sagews

#####################################################################
##                                                                 ##
## Distribution Statistique des Distances entre Antisites Voisins  ##
##                                                                 ##
#####################################################################

#Variables : 
    #r : distance entre sites, 
    #N : 0.22 nm^-3
    #p : probabilite d'avoir la distance r
    #rmax : distance la plus probable
r,N,p,rmax= var('r N p rmax')

p(r,N)=4*pi*r^2*N*exp(-N*4*pi*r^3/3)
p(r,N)

\newcommand{\Bold}[1]{\mathbf{#1}}4 \, \pi N r^{2} e^{-4/3*pi*N*r^3}

graph=plot(p(r,0.22),r,0,3)
graph.show()

#Distribution des distances entre antisites voisins dans CLN (en nm)
#Maximum atteint à rmax=0.898 nm avec une largeurà mi-hauteur de 0.84507 nm (voir plus bas)

#Recherche des max de la courbe
solve(diff(p(r,N),r),r)

\newcommand{\Bold}[1]{\mathbf{#1}}[ r == 1/4*(I*sqrt(3) - 1)*2^(2/3)/(pi^(1/3)*N^(1/3)), r == 1/4*(-I*sqrt(3) - 1)*2^(2/3)/(pi^(1/3)*N^(1/3)), r == 1/2*2^(2/3)/(pi^(1/3)*N^(1/3)), r == 0 ]

#valeur rmax de la distance la plus probable retenue (réelle et positive)
rmax=1/(2^(1/3)*pi^(1/3)*N^(1/3))

#valeur numérique en nm
rmaxPhy=1/(2^(1/3)*pi^(1/3)*0.22^(1/3))
n(rmaxPhy)

\newcommand{\Bold}[1]{\mathbf{#1}}0.897702601202902

# Valeur de la probabilité à mi hauteur pour évaluer la largeur a mi hauteur
pdemi=p(rmaxPhy,0.22)/2
n(pdemi)

\newcommand{\Bold}[1]{\mathbf{#1}}0.571923394612675

# Recherche des distances pour lesquelles p(r)=pdemi
# Avec la méthode de Newton, on cherche les zéros de la fonction :

p(r,N)-pdemi

\newcommand{\Bold}[1]{\mathbf{#1}}4 \, \pi N r^{2} e^{-4/3*pi*N*r^3} - 0.603681073679769 \, \pi^{\left(\frac{1}{3}\right)} 2^{\left(\frac{1}{3}\right)} e^{-0.666666666666667}

diff(4*pi*r^2*N*e^(-(4*pi*r^3*N/3)) - 0.619880395326*pi^(1/3)/2^(2/3),r)

\newcommand{\Bold}[1]{\mathbf{#1}}-16 \, \pi^{2} N^{2} r^{4} e^{-4/3*pi*N*r^3} + 8 \, \pi N r e^{-4/3*pi*N*r^3}

def R2(R):
    return 4*pi*R^2*0.22*e^(-(4*pi*R^3*0.22/3)) - 0.619880395326*pi^(1/3)/2^(2/3)
def dR2(R):
    return n(8*pi*R*0.22*e^(-(4*pi*R^3*0.22/3)) - 16*pi^2*R^4*0.22^2*e^(-(4*pi*R^3*0.22/3)))
def Phi(x):
    return n(x-R2(x)/dR2(x))

def Newton(depart):
    R1=depart
    R3=0
    err=2
    tour=0
    while (err>0.0001 and tour<10000):
        R3=Phi(R1)
        err=abs(R3-R1)
        R1=R3
        tour=tour+1
    if tour==10000:
        return I
    else:
        return R1

#Calcul de la largeur à mi hauteur de p :
rdemi2=Newton(1)#premier point
rdemi1=Newton(0.2)#second point
sigma=abs(rdemi2-rdemi1)
sigma

\newcommand{\Bold}[1]{\mathbf{#1}}0.845078226562872