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
* gyroradius (rg, rgp, rge)
* 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')
0.9579 rad/s
'''
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')
#0.10 rad/s
#"""
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')
#1758.8 rad/s
#"""
#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)
[docs]def gyroradius(mass, charge, bfield, vel):
r""" Gyroradius.
.. math::
r_g = \frac{v}{\omega_g} = \frac{mv}{qB}
"""
w = wg(mass, charge, bfield)
return vel / w
gyroradius_u = u.make_unitful_function(gyroradius, [u.kg, u.C, u.T, u.m / u.s], u.m)
rg = gyroradius
rg_u = gyroradius_u
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)
#r"""Proton gyro radius.
#
#: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],
u.rad / u.s)
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)
[docs]def escape_speed(mass_kg, radius_km):
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')))
r = radius_km * 1e3
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)