/** This is a library for equations about thermodynamics. */

// We need this to obtain the erf function in maxwellDistributionProbability.
// THINK ABOUT:  Dynamically load this using eval["use statistics.frink"] in
// that function?
use statistics.frink

/** The following are the Maxwell-Boltzmann equations for particle speeds in
    an ideal gas.

    See:
    https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution

    Many of these equations require specifying the mass of a molecule of the
    gas.  This is most easily obtained with the elements.frink library.  For
    example, to obtain the mass of a nitrogen molecule (N₂), you could use:

    use elements.frink
    N2mass = 2 Element.N.atomicMass
*/


/** Calculate the mean particle speed of an ideal gas at temperature T with
    a particle mass of m.  This equation is equivalent to:

    sqrt[8 k T / (pi m)]

    Note that this is the *mean* particle speed in 3 dimensions, not the
    *most likely* particle speed as calculated by mostLikelyParticleSpeed.
*/
meanParticleSpeed[T is temperature, m is mass] :=
{
   a = maxwellDistributionParameter[T, m]
   return 2 a sqrt[2/pi]
}


/** Calculate the most likely particle speed in an ideal gas at temperature T
    with a particle mass of m.  This equation is equivalent to:

    sqrt[R T / M]

    Where R is the gas constant (Frink calls this R or gasconstant) and M is
    the *molar mass* of the substance, that is the particle mass multiplied
    by "avogadro" and has units of mass/mol.

    Note that this is different from the mean speed calculated by
    meanParticleSpeed.  The most likely particle speed is amaller than
    the mean speed by a factor of sqrt[pi]/2 or ~0.88622 or, in other words,
    mostLikelyParticleSpeed = sqrt[pi]/2 meanParticleSpeed 
*/
mostLikelyParticleSpeed[T is temperature, m is mass] :=
{
   return sqrt[2 k T / m]
}


/** The Maxwell distribution parameter (often called "a") for an ideal gas
    of temperature T and particle mass m.  This is used by many other
    equations.
*/
maxwellDistributionParameter[T is temperature, m is mass] :=
{
   return sqrt[k T / m]
}


/** Returns the probability density function (PDF) of a temperature distribution
    for a given target velocity targetV.

    The mean of the distribution is the most likely particle speed (as
    calculated by mostLikelyParticleSpeed above, and *not* the mean speed.)
    TODO:  Verify this.  And derive the peak.

    Note that this function is not particularly useful other than to draw the
    shape of a pretty (skewed) bell curve.

    The results have the units of inverse velocity, that is seconds/meter
*/
maxwellDistributionDensity[T is temperature, m is mass, targetV is velocity] :=
{
   a = maxwellDistributionParameter[T, m]
   return sqrt[2/pi] targetV^2 / a^3 exp[-targetV^2 / (2 a^2)]
}


/** Returns the cumulative density function (CDF) of a temperature distribution
    for a given target velocity targetV.  This gives the probability that the
    velocity of a particle is targetV *or less*.  The result is a dimensionless
    number between 0 and 1.

    You can use this function to find the probability that a particle velocity
    is between v1 (lower) and v2 (upper) by:

    maxwellDistributionProbability[T, m, v2] - maxwellDistributionProbability[T, m, v1]

    This distribution returns 0.5 at about targetV = 
    1.0876520317581671916 mostLikelyParticleSpeed[T, m]
    or
    1.0876520317581671916 sqrt[pi] meanParticleSpeed[T, m] / 2 

    (the constant 1.087... is a constant numerical solution to
      erf[x] − 2 x exp⁡[−x^2] / sqrt[π] = 1/2
     and probably does not have a closed-form solution.
    )

    THINK ABOUT:  Do we want to make a "medianParticleSpeed" function that
    gives the velocity to which this function returns 1/2 as given above?
*/
maxwellDistributionProbability[T is temperature, m is mass, targetV is velocity] :=
{
   a = maxwellDistributionParameter[T, m]
   return erf[targetV / (sqrt[2] a)] - sqrt[2/pi] targetV/a exp[-targetV^2 / (2 a^2)]
}


/** Calculate the mean kinetic energy of a molecule at temperature T in 3
    dimensions.  Note that this is independent of the mass of the molecule!
*/
meanParticleKineticEnergy[T is temperature] :=
{
   return 3/2 k T
}


/** Calculate the mean momentum of a molecule at temperature T.

    This is equivalent to:
    p = sqrt[8 k T m / pi]
*/
meanParticleMomentum[T is temperature, m is mass] :=
{
   return m meanParticleSpeed[T, m]
}