# Source code for irfpy.util.physics

r''' Collection of physical formulae

.. codeauthor:: Yoshifumi Futaana

A simple command line calculator for physical formulae.

.. seealso::

- :mod:irfpy.util.unitq: the unit handling
- :mod:irfpy.util.constant: the constant with and without unit

**Supported formula**

* temperature_2_pressure (vth, vthp, vthe)
* pressure_2_temperature
* dflux_2_vdf (j2f)
* vdf_2_dflux (f2j)
* kinetic_energy
* energy_2_velocity
* temperature_2_thermal_velocity
* gyrofrequency (wg, wgp, wge)
* gyrotime
* plasmafrequency (wp, wpp, wpe)
* plasmafrequency_ordinary (fp, fpp, fpe)
* inertial_length
* skin_depth
* alfven_speed
* ion_acoustic_speed
* debye_length
* escape_speed
* mean_free_path
* mean_free_path_kin
* mean_free_path_vis
* orbital_semimajor_axis
* orbital_period

Each functions has (in most of the cases) the associated functions with _u behiind.
This is for the unitful-version of each method using :mod:irfpy.util.unitq module.

**Example**

For example, you can get the kinetic energy of a mass of 5 kg particle with 3 m/s object as 22.5 J.

>>> import irfpy.util.physics as f
>>> print(f.kinetic_energy(5, 3))
22.5

Most of them has two variants.

- MKSA (or other pre-defined) unit system
- User-defined unit system at runtime using  :mod:irfpy.util.unitq system.

Usually, the second system has usually a suffix of _u.

For example, if you want to know the energy (in Joule) for a mass of 4 kg moving with 15 m/s.

>>> print(f.kinetic_energy(4, 15))   # Kinetic energy for 4kg mass point with 15m/s speed. 450J.
450.0

But the same works as follows.

>>> from irfpy.util import unitq as u
>>> print(f.kinetic_energy(u.nit('4 kg'), u.nit('15 m/s')).rescale(u.J))
450.0 J

Even, you can use another unit system.

>>> print(f.kinetic_energy(u.nit('4000 g'), u.nit('0.015 km/s')).rescale(u.J))
450.0 J

The :mod:irfpy.util.unitq module contains the definitions of various units and the constant.

>>> from irfpy.util.unitq import k
>>> print(f.kinetic_energy(k.proton_mass, u.speed_of_light))
0.5 m_p*c**2

Also note that the :mod:irfpy.util.constant module also contains the definitions of various constants.

'''
import numpy as _np
import functools as _functools

from sympy.physics.mechanics.functions import inertia

from irfpy.util import constant as k
from irfpy.util import unitq as u

[docs]def temperature_2_pressure(kelvin, per_cm3):
r""" Temperature to pressure

The temperature in Kelvin under the density per_cm3 are
converted to the pressure in nPa.

.. math::

p \mathrm{[nPa]}\times 10^{-9} = (n\mathrm{[cm}^{-3}\mathrm{]}\times 10^6)kT

where k is the Boltzman constant.

>>> loschmidt = 2.6868e19    # 2.69e19 /cm3 -- Loschmidt's constant
>>> temp_ideal = 273.15      # 0 degreeC
>>> pressure_at_surface = temperature_2_pressure(273.15, 2.6868e19)  # This should be the ground pressure, 1 bar ~= 1e5 Pa ~= 1e14 nPa
>>> print('{:.3e}'.format(pressure_at_surface))
1.013e+14

Also, unitful version.

>>> loschmidt = 2.6868e19 / (u.cm ** 3)
>>> temp_ideal = 273.15 * u.K
>>> pressure_at_surface = temperature_2_pressure_u(temp_ideal, loschmidt)
>>> print('{:.3e}'.format(float(pressure_at_surface.r(u.nPa))))
1.013e+14

"""
n = per_cm3 * 1e6    # /cm3 to /m3
kb = k.kB         # Boltzman constant
t = kelvin

p = n * kb * t    # Pressure in Pa
p /= 1e-9         # Pressure in nPa

return p

temperature_2_pressure_u = u.make_unitful_function(temperature_2_pressure, [u.K, u.cm**-3], u.nPa)
""" Unitful version of :meth:temperature_2_pressure
"""

[docs]def pressure_2_temperature(nPa, per_cm3):
r""" Pressure to temperature.

See also :meth:temperature_2_pressure.

>>> p = 1e5 * 1e9      # 1 bar as nPa
>>> loschmidt = 2.6868e19   # in /cm3
>>> temperature = pressure_2_temperature(p, loschmidt)
>>> print('{:.3f}'.format(temperature))   # Should be around 0degC
269.576

Unitful version also exists

>>> p = 1 * u.bar      # 1 bar
>>> loschmidt = 2.6868e19 / u.cm **3
>>> temperature = pressure_2_temperature_u(p, loschmidt)
>>> temperature = float(temperature.rescale(u.K))
>>> print('{:.3f}'.format(temperature))
269.576
"""
p = nPa * 1e-9   # Pressure in Pa
kb = k.kB
n = per_cm3 * 1e6   # /cm3 to /m3

t = p / n / kb

return t

pressure_2_temperature_u = u.make_unitful_function(pressure_2_temperature, [u.nPa, u.cm**-3], u.K)
""" Unitful version of :meth:temperature_2_pressure
"""

[docs]def dflux_2_vdf(j, mass, energy):
r''' Differnetial flux to velocity distribution function

This is MKSA version.  For plasma-unit version (but still unitless),
you may refer to :meth:dflux_2_vdf_pl.
For unitq version, you may refere to "math:dflux_2_vdf_u.

The relation is

.. math::

f = \frac{m^2}{2E} j

Let us consider the differential flux of j=1e-5 /cm2 eV sr s.
The particle is proton (for simplicity mass is assumed to be 1 amu),
and the energy is 500 eV.

>>> j = 1e+5   # /cm2 eV sr s
>>> m = 1      # amu
>>> e = 500    # eV
>>> print('%.2e' % dflux_2_vdf_pl(j, m, e))
1.07e-10

Now, we get the velocity distribution function of 1.1e-10 s^3/m^6.

Obviously, we also call via MKSA unit version.

>>> jmksa = j * 1e4 / k.qe
>>> print('%.2e' % jmksa)   # in /m2 J sr s
6.24e+27

>>> mmksa = m * k.amu   # kg
>>> emksa = e * k.qe   # J

>>> print('%.2e' % dflux_2_vdf(jmksa, mmksa, emksa))
1.07e-10

For interactive use, unitq module will be a significant help.

>>> j_u = 1e5 / u.cm**2 / u.eV / u.sr / u.s
>>> m_u = 1 * u.amu
>>> e_u = 500 * u.eV
>>> vdf_u = dflux_2_vdf_u(j_u, m_u, e_u)
>>> u.pprint(vdf_u, fmt='%.2e')
1.07e-10 s**3/m**6
'''

f = (mass ** 2) / (2 * energy) * j
return f

dflux_2_vdf_u = u.make_unitful_function(dflux_2_vdf, [u.m**-2 / u.J / u.sr / u.s, u.kg, u.J], u.s**3 / u.m**6)

j2f = dflux_2_vdf
''' Alias to :func:dflux_2_vdf
'''
j2f_u = dflux_2_vdf_u

[docs]def dflux_2_vdf_pl(j, mass, energy):
r''' Differnetial flux with plasma-frequently-used units

This provides the alias but different unit system to
:func:dflux_2_vdf.

The input unit is j as [/cm2 eV sr s],
mass for [amu], and energy for [eV].

The returned vdf is [s**3 / m**6].
'''
jmksa = j * 1e4 / k.qe
massmksa = mass * k.amu
energymksa = energy * k.qe
return dflux_2_vdf(jmksa, massmksa, energymksa)

j2f_pl = dflux_2_vdf_pl
''' Alias to :func:dflux_2_vdf_pl
'''

[docs]def vdf_2_dflux(f, mass, velocity):
r''' J from f. (MKSA)

.. math::

J = \frac{v^2}{m} f

>>> f = 1.07e-10    # [s^3 / m^6]
>>> j = vdf_2_dflux(f, 1.66e-27, 300e3)
>>> print('%.2e' % j)    # [/m^2 J sr s]
5.80e+27

The dflux should be converted back to the origical vdf.

>>> f2 = dflux_2_vdf(j, 1.66e-27, kinetic_energy(1.66e-27, 300e3))
>>> print('%.2e' % f2)
1.07e-10

The unitful version:

>>> f_u = 1.07e-10 * u.s ** 3 / u.m ** 6
>>> j_u = vdf_2_dflux_u(f_u, 1 * u.amu, 300 * u.km/u.s)
>>> u.pprint(j_u, fmt='%.2e')
5.80e+27 1/(m**2*s*sr*J)

>>> j_u.units = u.cu_diffflux
>>> u.pprint(j_u, fmt='%.2e')
9.29e+04 (1/cm**2 /sr /eV / s)
'''
return f * velocity ** 2 / mass

vdf_2_dflux_u = u.make_unitful_function(vdf_2_dflux, [u.s**3 / u.m**6, u.kg, u.m/u.s], 1 / u.m**2 / u.J / u.sr / u.s)

f2j = vdf_2_dflux
f2j_u = vdf_2_dflux_u

[docs]def kinetic_energy(mass, velocity):
r''' Return the kinetic energy.

.. math::

KE = \frac{1}{2}mv^2

>>> print(kinetic_energy(4, 15))   # Kinetic energy for 4kg mass point with 15m/s speed. 450J.
450.0
>>> print(kinetic_energy(u.nit('4 kg'), u.nit('15 m/s')).rescale(u.J))
450.0 J
'''
return mass * velocity ** 2 / 2.

kinetic_energy_u = u.make_unitful_function(kinetic_energy, [u.kg, u.m/u.s], u.J)
''' Unitful version of :func:kinetic_energy
'''

[docs]def kinetic_energy_proton(velocity_in_km_s):
""" Return the proton energy.

:param velocity_in_m_s: Velocity in km/s.
:return: The kinetic energy in eV.

Simple formulation: 1 eV energy proton is about 13.8 km/s

>>> print('{:.3f}'.format(kinetic_energy_proton(13.8)))
0.994

Solar wind (4e5 m/s) proton has energy of 835 eV.
The call below is the unitful version of the formulation.

>>> energy = kinetic_energy_proton_u(400000 * u.m / u.s)
>>> print('{:.3f}'.format(float(energy.nr(u.keV))))
0.835
"""
e_in_J = kinetic_energy(k.proton_mass, velocity_in_km_s * 1000.)   # in J
return e_in_J / k.e

kinetic_energy_proton_u = u.make_unitful_function(kinetic_energy_proton,
[u.km / u.s], u.eV)

[docs]def kinetic_energy_electron(velocity_in_km_s):
""" Return the electron energy

:param velocity_in_km_s: Velocity in km/s
:return: The kinetic energy in eV

Simple formulation: 1 eV energy electron is about 600 km/s

>>> print('{:.3f}'.format(kinetic_energy_electron(600)))
1.023

or one can use unitful version.

>>> energy = kinetic_energy_electron_u(600000 * u.m / u.s)
>>> print('{:.3e}'.format(float(energy.nr(u.J))))
1.640e-19

"""
e_in_J = kinetic_energy(k.electron_mass, velocity_in_km_s * 1000.)   # In J
return e_in_J / k.e

kinetic_energy_electron_u = u.make_unitful_function(kinetic_energy_electron,
[u.km / u.s], u.eV)

[docs]def energy_2_velocity(mass, energy):
r''' Return the velocity from the kinetic energy

Return the velocity from the kinetic energy.
This is the inversion function of :func:kinetic_energy.

>>> print(energy_2_velocity(4, 450.))   # 4kg mass with 450 J => v=15 m/s
15.0

You may also use unitful version.

>>> print(energy_2_velocity_u(u.nit('4 kg'), u.nit('450 J')))
15.0 m/s
'''
return _np.sqrt(2 * energy / mass)

energy_2_velocity_u = u.make_unitful_function(energy_2_velocity,
[u.kg, u.J], u.m / u.s)

[docs]def energy_2_velocity_proton(energy_eV):
""" Return the velocity (in km/s) of the proton with the given energy (eV)

:param energy_eV: Energy of proton in eV.
:return: Speed of proton in km/s.

>>> print('{:.2f}'.format(energy_2_velocity_proton(1)))    # 1 eV proton has 13.8 km/s
13.84

or use unit-version.

>>> speed = energy_2_velocity_proton_u(1 * u.eV)
>>> u.pprint(speed)
13.841 km/s
"""
energy_J = energy_eV * k.e
speed_m_s = energy_2_velocity(k.proton_mass, energy_J)
return speed_m_s / 1000.

energy_2_velocity_proton_u = u.make_unitful_function(energy_2_velocity_proton,
[u.eV], u.km / u.s)

[docs]def energy_2_velocity_electron(energy_eV):
""" Return the velocity (in km/s) of the electron with the given energy (eV)

:param energy_eV: Energy of electron in eV.
:return: Speed of electron in km/s.

>>> print('{:.2f}'.format(energy_2_velocity_electron(1)))    # 1 eV proton has 593 km/s
593.10

or use unit-version.

>>> speed = energy_2_velocity_electron_u(1 * u.eV)
>>> u.pprint(speed)
593.097 km/s
"""
energy_J = energy_eV * k.e
speed_m_s = energy_2_velocity(k.electron_mass, energy_J)
return speed_m_s / 1000.

energy_2_velocity_electron_u = u.make_unitful_function(energy_2_velocity_electron,
[u.eV], u.km / u.s)

[docs]def temperature_2_thermal_velocity(mass, temperature):
r''' Return the thermal velocity from temperature

.. math::

v_{th} = \sqrt{\frac{kT}{m}}

Thermal velocity is calculated from the temperature.

Roughly the proton with 1eV temperature has 9.8 km/s thermal velocity.

>>> print('%.3f' % temperature_2_thermal_velocity(1.67e-27, 11600))
9792.924

Unitful version also supported.

>>> print('%.3f' % temperature_2_thermal_velocity_u(k.mp_u,
...        1 * u.eVT).nr(u.nit('km/s')))
9.787
'''
return _np.sqrt(k.kB * temperature / mass)

temperature_2_thermal_velocity_u = u.make_unitful_function(
temperature_2_thermal_velocity, [u.kg, u.K], u.m / u.s)
''' Unitful version of :func:temperature_2_thermal_velocity.
'''

vth = temperature_2_thermal_velocity
vth_u = temperature_2_thermal_velocity_u

[docs]def gyrofrequency(mass, charge, magnetic_flux_density):
r''' Gyro frequency, :math:\omega_{g}.

.. math::

\omega_g = \frac{qB}{m}

>>> omega = gyrofrequency_u(k.mp_u, k.e_u, 10 * u.nT)
>>> u.pprint(omega, fmt='%.4f')
'''
return charge * magnetic_flux_density / mass

gyrofrequency_u = u.make_unitful_function(
gyrofrequency, [u.kg, u.C, u.T], u.rad / u.s)
_functools.update_wrapper(gyrofrequency_u, gyrofrequency)

wg = gyrofrequency
wg_u = gyrofrequency_u

wgp = _functools.partial(gyrofrequency, k.proton_mass, k.qe)
_functools.update_wrapper(wgp, gyrofrequency)
wgp_u = u.make_unitful_function(wgp, [u.T], u.rad / u.s)

#""" Return the proton gyro frequency
#
#:param bfield: Magnetic flux density (Tesla)
#:return: Proton gyro frequency in rad/s.
#
#>>> u.pprint(wgp_u(u.nit('1 nT')), fmt='%.2f')
#"""

wge = _functools.partial(wg, k.electron_mass, k.qe)
_functools.update_wrapper(wge, wg)
wge_u = u.make_unitful_function(wge, [u.T], u.rad / u.s)
#"""
#
#>>> u.pprint(wge_u(u.nit('10 nT')), fmt='%.1f')
#"""
#wge_u = u.make_unitful_function(wge, [u.T], u.rad / u.s)

[docs]def gyrotime(mass, charge, magnetic_flux_density):
r""" Gyro time

Return the gyration time.  :math:2\pi/\omega_g.
"""
return 2 * _np.pi / wg(mass, charge, magnetic_flux_density)

gyrotime_u = u.make_unitful_function(gyrotime, [u.kg, u.C, u.T], u.s)

.. math::

r_g = \frac{v}{\omega_g} = \frac{mv}{qB}
"""
w = wg(mass, charge, bfield)
return vel / w

rgp = _functools.partial(rg, k.proton_mass, k.qe)
_functools.update_wrapper(rgp, rg)
rgp_u = u.make_unitful_function(rgp, [u.T, u.m / u.s], u.m)
_functools.update_wrapper(rgp_u, rgp)
#
#:param b: magnetic field
#:param v: velocity
#"""

rge = _functools.partial(rg, k.electron_mass, k.qe)
rge_u = u.make_unitful_function(rge, [u.T, u.m / u.s], u.m)

[docs]def plasmafrequency(mass, charge, density):
r''' Plasma frequency, :math:\omega_{p}

.. math::

\omega_p = \sqrt{\frac{ne^2}{m\epsilon_0}}

Return the plasma frequency.

Unitless version is MKSA; you may also use the alias version wp.

>>> print('{:.1f}'.format(plasmafrequency(1.67e-27, 1.60e-19, 10e6)))
4160.9
>>> print('{:.1f}'.format(wp(1.67e-27, 1.60e-19, 10e6)))
4160.9
'''
return _np.sqrt(density * (charge ** 2) / (mass * k.eps0))

plasmafrequency_u = u.make_unitful_function(plasmafrequency,
[u.kg, u.C, u.m ** -3],

wp = plasmafrequency
wp_u = plasmafrequency_u

wpp = _functools.partial(wp, k.proton_mass, k.qe)
wpp_u = u.make_unitful_function(wpp, [u.m ** -3], u.rad / u.s)

wpe = _functools.partial(wp, k.me, k.qe)
wpe_u = u.make_unitful_function(wpe, [u.m ** -3], u.rad / u.s)

[docs]def plasmafrequency_ordinary(mass, charge, density):
r""" Plasma ordinary frequency, :math:f_p

.. math::

f_p = \frac{\omega_p}{2\pi}

For example, the electron ordinary plasma frequency is

.. math::

f_{pe} \sim 8980 \sqrt{n_e} \mathrm{Hz}

>>> n_e = 100 / (u.cm ** 3)
>>> f_pe =plasmafrequency_ordinary_u(k.me_u, k.qe_u, n_e).nr(u.kHz)
>>> print('{:.1f}'.format(float(f_pe)))
89.8
"""
return plasmafrequency(mass, charge, density) / (2 * _np.pi)

plasmafrequency_ordinary_u = u.make_unitful_function(plasmafrequency_ordinary,
[u.kg, u.C, u.m ** -3],
1 / u.s)

fp = plasmafrequency_ordinary
fp_u = plasmafrequency_ordinary_u
fpp = _functools.partial(fp, k.mp, k.qe)
fpp_u = u.make_unitful_function(fpp, [u.m ** -3], 1 / u.s)
fpe = _functools.partial(fp, k.me, k.qe)
fpe_u = u.make_unitful_function(fpe, [u.m ** -3], 1 / u.s)

[docs]def inertial_length(mass, charge, density):
"""

>>> d = inertial_length(1.67e-27, 1.60e-19, 10e6)
>>> print('{:.2f}'.format(d))
72049.89
"""
c = k.speed_of_light
return c / wp(mass, charge, density)

inertial_length_u = u.make_unitful_function(inertial_length,
[u.kg, u.C, u.m ** -3],
u.m)

[docs]def skin_depth(density):
return inertial_length(k.m_e, k.qe, density)

skin_depth_u = u.make_unitful_function(skin_depth,
[u.m ** -3], u.m)

[docs]def alfven_speed(mass, density, mag):
r""" Return Alfven speed

.. math::

v_a = \frac{B}{\sqrt{\mu0nm}}
"""
u0 = k.mu_0
dn = _np.sqrt(u0 * density * mass)
return mag / dn

alfven_speed_u = u.make_unitful_function(alfven_speed,
[u.kg, u.m ** -3, u.T], u.m / u.s)

[docs]def ion_acoustic_speed(mass, temp_i, temp_e, gamma_i, gamma_e):
r"""Ion acoustic speed

.. math::

v_c = \sqrt{\frac{\gamma_ekT_e+\gamma_ikT_i}{m_i}}

:math:\gamma_e=1 is frequently used.
"""
kb = k.kB
vc = _np.sqrt((kb * temp_e * gamma_e + kb * temp_i * gamma_i) / mass)
return vc

ion_acoustic_speed_u = u.make_unitful_function(ion_acoustic_speed,
[u.kg, u.K, u.K, u.dimensionless, u.dimensionless], u.m / u.s)

[docs]def debye_length(ne, temp_e, charge):
r""" debye length.

.. math::

\lambda_D = \sqrt{\frac{\epsilon_0 k T_e}{n_e q_e^2}}
"""
eps0 = k.eps0
kB = k.kB
inside =  eps0 * kB * temp_e / ne
return _np.sqrt(inside) / charge

debye_length_u = u.make_unitful_function(debye_length,
[u.m ** -3, u.K, u.C],
u.m)

r""" Return the escape speed in km/s.

.. math::

v_\mathrm{esc} = \sqrt{\frac{2GM}{R}}

For the Earth:

>>> me = 5.972e24   # Mass of Earth in kg
>>> re = 6371.      # Radius of Earth in km
>>> esc_e = escape_speed(me, re)    # Escape speed of Earth
>>> print('{:.1f}'.format(esc_e))
11.2

For the Sun:

>>> ms = 1.989e30 * u.kg
>>> rs = 695700 * u.km
>>> esc_s = escape_speed_u(ms, rs)
>>> print('%.1f' % esc_s.nr(u.km / u.s))
617.8
"""
gravcon = float(u.k.G.nr(u.nit('m**3/kg/s**2')))
return _np.sqrt(2 * gravcon * mass_kg / r) * 1e-3

escape_speed_u = u.make_unitful_function(escape_speed,
[u.kg, u.km],
u.km / u.s)

[docs]def mean_free_path(density_cm3, crosssection_cm2):
r""" Return the mean free path in cm.

A simple theory (scattering theory) gives the mean free path of

.. math::

\lambda = \frac{1}{n\sigma}

where :math:n is the number density and :math:\sigma is the
crosssection.

This is based on the scattering theory.

:param density_cm3: Density in cm^-3
:param crosssection_cm2: Cross section in cm^-2
:returns: Mean free path in the unit of cm.

>>> n = 2.5e19    # Approximate molecule number density at Earth surface.
>>> sigma = 1e-15   # Typical collision cross-section
>>> mfp = mean_free_path(n, sigma)
>>> print('{:.3e}'.format(mfp))
4.000e-05

Or derive with unit support.

>>> from irfpy.util import unitq as u
>>> n = 2.5e25 * u.nit('m^-3')    # 2.5 x 10^25 m^-3
>>> sigma = 1e-15 * u.nit('cm^2')  # 1 x 10^25 cm^2
>>> mfp = mean_free_path_u(n, sigma)
>>> mfp = mfp.r(u.m)    # The unit is rescaled to meter.
>>> print('{:.2e}'.format(mfp))     # The mean free path value as the unit of meter.
4.00e-07 m
"""
return 1 / (density_cm3 * 1e6) / (crosssection_cm2 * 1e-4) * 100

mean_free_path_u = u.make_unitful_function(mean_free_path,
[u.cm ** -3, u.cm ** 2],
u.cm)

[docs]def mean_free_path_kin(density_cm3, crosssection_cm2):
r""" Return the mean free path in cm.

Kinetic theory (including temperature) suggests a correction
by a factor of :math:\sqrt{2} compared to :func:mean_free_path.

.. math::

\lambda = \frac{1}{\sqrt{2}n\sigma}

where :math:n is the number density and :math:\sigma is the
crosssection.

:param density_cm3: Density in cm^-3
:param crosssection_cm2: Cross section in cm^-2
:returns: Mean free path in the unit of cm.

>>> n = 2.5e19    # Approximate molecule number density at Earth surface.
>>> sigma = 1e-15   # Typical collision cross-section
>>> mfp = mean_free_path_kin(n, sigma)
>>> print('{:.2e}'.format(mfp))  # Returned value is in cm.
2.83e-05

Or derive with unit support.

>>> from irfpy.util import unitq as u
>>> n = 2.5e25 * u.nit('m^-3')    # 2.5 x 10^25 m^-3
>>> sigma = 1e-15 * u.nit('cm^2')  # 1 x 10^25 cm^2
>>> mfp = mean_free_path_kin_u(n, sigma)
>>> mfp = mfp.r(u.m)    # The unit is rescaled to meter.
>>> print('{:.2e}'.format(mfp))     # The mean free path value as the unit of meter.
2.83e-07 m
"""
return mean_free_path(density_cm3, crosssection_cm2) / _np.sqrt(2)

mean_free_path_kin_u = u.make_unitful_function(mean_free_path_kin,
[u.cm ** -3, u.cm ** 2],
u.cm)

[docs]def mean_free_path_vis(mu_Pa_s, p_Pa, T_K, m_amu):
r""" Based on kinetic theory, the m.f.p is calcuated from viscosity.

.. math::

l = \frac{\mu}{p}\sqrt{\frac{\pi k_B T}{2m}}

Here, :math:\mu is the viscosity, p is the pressure,
:math:k_B is the Boltzmann constant, T is the temperature,
and m is the mass of the gas.

https://en.wikipedia.org/wiki/Mean_free_path#Kinetic_theory_of_gases

It is not easy to get the viscosity for specific gas.
For example of the air at 20 degreeC, the viscosity is ~1.8x10^-5 Pa s.

:param mu_Pa_s: Viscosity in the unit of Pa s. Just a note
1 Pa s = 10 P = 1000 cP.
:param p_Pa: Pressure in the unit of Pa.
:param T: Temperuter in Kelvin.
:param m_amu: Mass in atomic mass unit.
:returns: Mean free path in the unit of cm.

>>> mu = 1.8e-5    # Pa s
>>> p = 101100   # Pascal
>>> T = 20 + 273.15   # Kelvin
>>> m = 28    # amu

>>> l = mean_free_path_vis(mu, p, T, m)
>>> print('{:.2e}'.format(l))
6.58e-06

Unitful version:

>>> mu = 0.018 * u.nit('centipoise')    # Air viscosity
>>> p = 1013 * u.nit('hPa')    # Air pressure
>>> T = (20 + 273.15) * u.nit('K')    # Unit does not support the conversion of degC <-> degK !!
>>> m = 28 * u.amu

>>> l = mean_free_path_vis_u(mu, p, T, m)    # Calculate the mean free path
>>> l = l.r(u.angstrom)    # Rescale to the unit of Ångström
>>> print('{:.2e}'.format(l))
6.57e+02 angstrom
"""
mu = mu_Pa_s
p = p_Pa
pi = _np.pi
kb = k.kB    # Boltzman constant
T = T_K
m = m_amu * k.amu

l_m = (mu / p) * _np.sqrt(pi * kb * T / 2 / m)

return l_m * 100    # Return in cm

mean_free_path_vis_u = u.make_unitful_function(mean_free_path_vis,
[u.Pa * u.s,
u.Pa,
u.Kelvin,
u.amu],
u.cm)

[docs]def orbital_semimajor_axis(period, mass):
r""" Calculate the semi-major axis from the period of satellite.

:param period: Period in [s]
:param mass: Mass of the central body in [kg]
:returns: Semi-jmajor axis (the radius for a circular orbit) in [km]

Using the Kepler's law:

.. math::

\frac{T^2}{a^3} = \frac{4\pi^2}{GM}

# Lunar orbiter with period of 7000s.  The semi-major axis is ~1825 km (100 km alt).
>>> from irfpy.util import planets as pla
>>> a = orbital_semimajor_axis(7000, pla.moon['mass'])
>>> print(f'{a:.1f}')
1825.8

# Moon around the earth. 27.5 day period give the distance of ~60 Rm (384000 km)
>>> a_u = orbital_semimajor_axis_u(27.4 * u.day,
...                      pla.earth['mass'] * u.kg)
>>> a = a_u.nr(u.km)    # Represented by days in float.
>>> print(f'{a:.1f}')
383946.1
"""
gravity_constant = k.G     # Gravity constant in MKSA
gm = gravity_constant * mass

t2 = period * period

pi2 = _np.pi ** 2

a = _np.power(t2 * gm / 4 / pi2, 1/3.0)

a = a / 1000   # Convert to km.

return a

orbital_semimajor_axis_u = u.make_unitful_function(orbital_semimajor_axis,
[u.s, u.kg], u.km)

[docs]def orbital_period(a, m):
r""" Calculate the orbital period.

:param a: Semi-jmajor axis (the radius for a circular orbit) in [km]
:param m: Mass of the central body in [kg]
:returns: Period in [s]

Using the Kepler's law:

.. math::

\frac{T^2}{a^3} = \frac{4\pi^2}{GM}

# Lunar orbiter at 100 km.  It is about 2 hours (~7050s)
>>> from irfpy.util import planets as pla
>>> p = orbital_period(1834, pla.moon['mass'])
>>> print(f'{p:.1f}')
7046.9

# Moon around the earth. ~27 days
>>> p_u = orbital_period_u(pla.moon['semi-major-axis'] * u.km,
...                      pla.earth['mass'] * u.kg)
>>> p = p_u.nr(u.day)    # Represented by days in float.
>>> print(f'{p:.1f}')
27.4
"""
gravity_constant = k.G     # Gravity constant in MKSA
pi = _np.pi
a_m = a * 1000     # In meter

period = 2 * pi * a_m * _np.sqrt(a_m / gravity_constant / m)

return period

orbital_period_u = u.make_unitful_function(orbital_period,
[u.km, u.kg], u.s)
`