-----Gibbs distribution 1-------
M=1000
me=runif(M)
hist(me)
N=20000
for (n in 1:N) {
  k=floor(M*runif(1))+1
  l=floor(M*runif(1))+1
  tot=me[k]+me[l]
  r1=runif(1)
  me[k]=tot*r1
  me[l]=tot-me[k]
  hist(me,breaks="FD")
}
hist(me)
----------------------------------

-----Gibbs distribution 2-------
M=1000
me <- rep(1, times = M)
N=20000
for (n in 1:N) {
  k=floor(M*runif(1))+1
  l=floor(M*runif(1))+1
  tot=me[k]+me[l]
  r1=runif(1)
  me[k]=tot*r1
  me[l]=tot-me[k]
  hist(me,breaks="FD")
}
hist(me)
--------------------------------

--Maxell-Boltmann distribution--
M=2000
ve=runif(M)
hist(ve)
N=50000
R=100
for (r in 1:R) {
for (n in 1:N) {
  k=floor(M*runif(1))+1
  l=floor(M*runif(1))+1
  t=pi*runif(1)/2.0
  f=pi*runif(1)/2.0
  tot1=sqrt((ve[k]*sin(t))^2+(ve[l]*cos(f))^2)
  tot2=sqrt((ve[k]*cos(t))^2+(ve[l]*sin(f))^2)
  ve[k]=tot1
  ve[l]=tot2
}
hist(ve)}

hist(ve,breaks="FD")
hist(ve)
----------------------------------