# Premessa: come installare la libreria che contiene la distribuzione di Maxwell-Boltzmann
# Andato qui:     https://cran.r-project.org/web/packages/VGAM/
# ti colleghi all'URL (dal sito o battendo direttamente):
# https://www.stat.auckland.ac.nz/~yee/VGAM
# Fai (usando il link) il   download VGAM  e poi da R operi sul menu
# "Packages" e quindi "Install package(s) from local zip files..." e
# scegli dalla cartella di Download il file   VGAM.zip
# Poi clicchi  "Packages", "Load package..." e scegli VGAM
#
# Un codice per generare velocità distribuite secondo Maxwell-Boltzmann
# Definiamo l'unità di massa atomica
u <- 1/6.02e26
# Definiamo la massa dell'atomo di elio
m <- 4*u
# Definiamo la costante di Boltzmann
k <- 1.38e-23
# Lavoriamo a temperatura ambiente (27°C)
T <- 300
# Facciamo la chiamata della libreria VGAM
library(VGAM)
# Qui sono presenti le definzioni sia di rmaxwell(N,R) che di dmaxwell(x,R)
# che, rispettivamente, genera a caso N valori distribuiti secondo Maxwell-Boltzmann
# e rappresenta la corrispondente funzione di densità, al variare del valore di R.
# Il valore di R è m/(k*T).
R <- m/(k*T)
# generiamo le velocità per 1 milione di molecole
v <- rmaxwell(1e6, R)
# L'istogramma di distribuzione:
hist(v, probability=TRUE)
# La distribuzione teorica:
g <- function(x) dmaxwell(x, R)
curve(g, add=TRUE, col="red", lwd=2)
# Il grafico di sovrappone "bene" all'istogramma

# Infine una verifica a posteriori: dalle velocità ricaviamo le energie cinetiche
EC <- 0.5*m*v^2
# L'energia cinetica media per molecola è 1.5*k*T, quindi stimiamo T
T <- mean(EC)/1.5/k; T
# Abbiamo ritrovato, circa, 300
# Volendo si può esplicitare la funzione di densità:
f <- function(x) sqrt(2/pi)*R^(3/2)*x^2*exp(-1/2*R*x^2)
# Si vede che il grafico si sovrappone al precedente:
curve(f, add=TRUE, lty=3, lwd=2)

#
# Una animazione, da poche centinaia ad un milione di valori generati