kill(all)$

Maxwell - Boltzmann distribution .(appendix 1),version 1.0
peter.vlasschaert@gmail.com,23/02/2019

load packages
=============

load(lrats)$

approximation :  log (x!)  = x*log(x) - x , 'stirling approximation'
===============================
rem : log(x)   base :e 

integrate(log(x), x);

W_1 => n_1
W_2 => n_2
---------------
W_i  => n_i
==================================
W_i : energy @  level i   with n_i particles , i=1,2,..k levels
total energy :  W_i for  'k' levels  = W

total number particles = N  for  'k' levels
=================================

p1:'N = sum(n[i], i, 1, k);

total energy = W  for 'k' levels
========================

p2:'W =  sum(n[i]*W[i], i, 1, k);

How to fill energy levels with particles.
================================

  N     ways choose  ( = first particle)
  N-1  ways choose  (= second particle)   
---------------------------------------------------------------------
  N-i  ways choose (= i th particle)      
======================================
boleaan operation :  And  
--------------------------------------------------------------------           
result : N (and )  N-1 (and)  N-2  .... (and) N-i        
--------------------------------------------------------------------
math : N*(N-1)*(N-2) *......*(N-i)   =  N!/(N-i+1) !

by example : first three particle a,b,c ( with energy W_1)
==================================
i=2 : for three particles  ( = n_1)

p3: N!/(N-3)!;

We have :  3 !  permutations
=======================
a,b,c
a,c,b
b,a,c
b,c,a
c,a,b
c,b,a
------------------------------------------
3! = 6 ways to select

p4:p3/3!;

general fill up level W_1 with n_1 particles
===================================

p5a:  N!/((n[1])! *(N-n[1])!);

general fill up level W_2 with n_2 particles
===================================
N  replaced by  N - n_1

p5b: subst([n[1] = n[2]],p5a)$
p5b1:subst([N=N-n[1]],p5b);

general fill up level W_3 with n_3 particles
===================================
N-n_1  replaced by  N - n_1-n_2

p5c1:subst([N=N-n[2]],num(p5b1));
p5c2:subst([n[2]=n[3]],denom(p5b1))$
p5c3:subst([N = N-n[2]],p5c2);
p5c4:p5c1/p5c3;

probability ( when filling particles 'is the same' for 'every energy level')
===========================================================
example : three energy levels  with  W_1  => n_1
                                                                   W_2  => n_2
                                                                   W_3  => n_3

p6:P =p5a*p5b1*p5c4;

rem   from p1  (  use  k =3 ,three   levels)

we make function : 'p1' 
===================

particles(klevel):= sum(n[i], i, 1, klevel);

assume(n_1>0,n_2>0,n_3>0);

p6a:lhs(p1) = particles(3);

p6b:lratsubst([lhs(p6a)=rhs(p6a)],p6);

p6c:subst([rhs(p6a)=lhs(p6a)],num(rhs(p6b)))$
p6c1:lhs(p6b) = p6c/denom(rhs(p6b));

probability ( when filling particles 'is the same' for 'every energy level')
===========================================================
example : three energy levels  with  W_1  => n_1
                                                                   W_2  => n_2
                                                                   W_3  => n_3
=> condition removed   ' same '  = >   'different'
============================================================
reason :  first level can be have higher probability than the higher levels
-----------------------------------------------------------------------------------------------------------
probability one   particle    level W_i   = >  (g_i)^1
                      two   particles level W_i   =>   (g_i)*(g_i) = (g_i)^2
i = the level for particles
-------------------------------------------------------------------------------------------------------------
i=1  =>  'one particle '  g_1   = >  'n_1 particles'   => (g_1)^(n_1)
i=2  =>  'one particle '  g_2   = >  'n_2 particles'   => (g_2)^(n_2)
------------------------------------------------------------------------------------------------------------
i=k  =>  'one particle '  g_k   = >  'n_k particles'   => (g_k)^(n_k)

correction ( for  probability P)
=========================

p7:product((g[i])^(n[i]) ,i, 1, k);

total probability ( after use correction 'p7') , example ' three energy levels'
===================================

p7a:subst([k=3],p7);

p7b:lhs(p6c1) = rhs(p6c1)*p7a;

p7c:subst([p6a],p7b);

last correction : remove permutation from above 
-----------------------------------------------------------------------------------
because : all particle or the same in structure and
no interaction between each other .( perfect = ideal gas).

p7d:lhs(p7b) = rhs(p7b)/N!;

how to find the equilibrium ' P' : Max(P)
-----------------------------------------------------------
constraint : from 'p1' and 'p2'
-----------------------------------------------------------
solve : lagrange multiplicator :α,β

result : 'particles , i th level with energy W_i'
=====

p8: n[i] = g[i]*exp(-α-β*W[i]);

p9:sum(p8,i,1,k);

example : three levels  k=3
==================

p9a:subst([k=3],p9),expand;

p9b:N = expand(sum(g[i]*%e^(-W[i]*β-α),i,1,3));

p9c:rhs(p9b)*exp(α)$

p9d:lhs(p9b) = expand(p9c)*exp(-α);

p9e:lhs(p9d)=subst([part(p9d,2,2) = Z],rhs(p9d));

p9e1:Z= part(p9d,2,2) ;

p9f:p9e/Z;

New equation : 'p8'
================

p10:lratsubst([rhs(p9f) = lhs(p9f)],p8);

because from mathematics :  exp(x)
==============================
===========
x = number :
===========

==============================
=>  W_i (energy = SI units 'Joule')
result : 1/β = energy 
==============================

p11:lhs(p10) = n0[i]*exp(-W[i]/(kb*T));

legend 
======

-------------------------------------------

T     = absolute temperature
kb   = Boltzmann constante
W_i = energy @ level i
n0_i= bulk factor
i       = i the level.

-------------------------------------------

How to find :n0_i  : ' from  p11' ( = bulk factor)
======================================

p11a:n0[i] =subst([β = 0],rhs(p10));

assume(kb>0);

p12:'limit(rhs(p11),T,inf) =limit(rhs(p11),T,inf); 

p12a:lhs(p11) = rhs(p12);

rem : T → infinity : kinetic energy , leaves the energy levels.(bulk factor)