""" Martian exosphere models.
Several density models for the Martian upper atmosphere.
- :class:`ISSI2008_SolarMax`: Solar max neutral density model for Mars
- :class:`ISSI2008_SolarMin`: Solar min neutral density model for Mars
- :class:`ZetaChamberlain1963` Synthesis of three implementations of the partition function.
Sample scripts are also available at :mod:`snippet_mars.exosphere_max`, :mod:`snippet_mars.exosphere_min`, :mod:`snippet_mars.exosphere_demo`
"""
import numpy as _np
[docs]class ISSI2008_SolarMax:
r""" The neutral upper atmosphere profiles of CO2, O, and H used in the ISSI SWIM model challenge (solar max)
.. codeauthor:: X.-D. Wang
Height profiles of five species (CO2, hot_O, cold_O, hot_H, cold_H) are obtained.
TODO: include the partition function into the solar maximum profiles too.
>>> from irfpy.mars.exosphere import ISSI2008_SolarMax as exosph
>>> co2_250 = exosph.CO2(250) # CO2 density at 250 km in /cm3.
>>> print('{:.3e}'.format(co2_250))
1.089e+07
>>> oc_250 = exosph.O_cold(250) # Cold O denisty at 250 km.
>>> print('{:.3e}'.format(oc_250))
1.655e+07
>>> oh_250 = exosph.O_hot(250) # Hot O density at 250 km.
>>> print('{:.3e}'.format(oh_250))
1.761e+04
>>> hc_250 = exosph.H_cold(250) # Cold H density at 250 km.
>>> print('{:.3e}'.format(hc_250))
2.923e+01
>>> hh_250 = exosph.H_hot(250) # Hot H density at 250 km.
>>> print('{:.3e}'.format(hh_250))
2.835e+04
The original models are fed by the MTGCM model.
This is an analytical copy obtained by manual scaling the plots in Brain's presentation.
Following are the original texts on http://www.issibern.ch/teams/martianplasma/
Neutral Atmosphere:
All information will be based on results from Steve Bougher's MTGCM,
and published results from Jean-Yves Chaufray for hydrogen densities /?? Reference preferable ??/.
The MTGCM was run for Ls=270, F10.7=105 at Mars.
.. Below Needs another module
.. Gif images show the MTGCM output in red,
.. and the fit/extrapolation in white. z is referenced from the Martian surface, in km.
.. Density has units of cm^-3.
CO2 number density (:meth:`CO2`):
5.88e18 * exp( -z / 7.00 ) + 3.55e13 * exp( -z / 16.67 )
'Cold' O number density (:meth:`O_cold`):
2.33e13 * exp( -z / 12.27 ) + 2.84e09 * exp( -z / 48.57 )
'Cold' H number density (:meth:`H_cold`):
1e3 * exp[ 9.25e5 * ( 1 / (z+3393.5) - 1/3593.5) ]
Neutral Temperature (?? Not implemented ??):
-161.13 * exp[ -0.5 * ( (z-112.6) / 25.25 )^2 ] + 291.78
Following are the description of the extra hot corona components of hydrogen and oxygen.
Hot corona / Extended exosphere:
'Hot' Neutral Hydrogen (:meth:`H_hot`):
3e4 * exp[ 1.48e4 * ( 1 / (z+3393.5) - 1/3593.5) ]
'Hot' Neutral Oxygen (:meth:`O_hot`):
1.56e4 * exp( -z / 696.9) + 2.92e3 * exp( -z / 2891.) + 5.01e4 * exp( -z / 99.19)
"""
[docs] @staticmethod
def CO2(z):
""" CO2 density model for solar max at Mars.
:param z: Altitude in km.
:return: CO2 density in /cm3.
"""
n = 5.88e18 * _np.exp(-z / 7.00) + 3.55e13 * _np.exp(-z / 16.67)
return n
[docs] @staticmethod
def O_cold(z):
""" Cold O density model for solar max at Mars.
:param z: Altitude in km.
:return: Cold O density in /cm3.
"""
n = 2.33e13 * _np.exp(-z / 12.27) + 2.84e09 * _np.exp(-z / 48.57)
return n
[docs] @staticmethod
def H_cold(z):
""" Cold H density model for solar max at Mars.
:param z: Altitude in km.
:return: Cold H density in /cm3.
"""
n = 1e3 * _np.exp(9.25e5 * (1 / (z + 3393.5) - 1 / 3593.5))
return n
[docs] @staticmethod
def H_hot(z):
""" Hot H density model for solar max at Mars.
:param z: Altitude in km.
:return: Hot H density in /cm3.
"""
n = 3e4 * _np.exp(1.48e4 * (1 / (z + 3393.5) - 1 / 3593.5))
return n
[docs] @staticmethod
def O_hot(z):
""" Hot O density model for solar max at Mars.
:param z: Altitude in km.
:return: Hot O density in /cm3.
"""
n = 1.56e4 * _np.exp(-z / 696.9) + 2.92e3 * _np.exp(-z / 2891.) + 5.01e4 * _np.exp(-z / 99.19)
return n
[docs]class ISSI2008_SolarMin:
''' The neutral upper atmosphere profiles of CO2, O, and H used in the ISSI SWIM model challenge (solar min).
.. codeauthor:: X.-D. Wang
Height profiles of the number densities for five species (CO2, hot_O, cold_O, hot_H, cold_H) are obtained.
>>> from irfpy.mars.exosphere import ISSI2008_SolarMin as exosph
>>> co2_250 = exosph.CO2(250) # CO2 density at 250 km in /cm3.
>>> print('{:.3e}'.format(co2_250))
5.950e+05
>>> oc_250 = exosph.O_cold(250) # Cold O denisty at 250 km.
>>> print('{:.3e}'.format(oc_250))
4.472e+06
>>> oh_250 = exosph.O_hot(250) # Hot O density at 250 km.
>>> print('{:.3e}'.format(oh_250))
6.599e+03
>>> hc_250 = exosph.H_cold(250) # Cold H density at 250 km.
>>> print('{:.3e}'.format(hc_250))
1.358e+05
>>> hh_250 = exosph.H_hot(250) # Hot H density at 250 km.
>>> print('{:.3e}'.format(hh_250))
1.826e+04
The original models are fed by the MTGCM model.
This is an analytical copy obtained by manual scaling the plots in Brain's presentation.
Following are the original texts on http://www.issibern.ch/teams/martianplasma/
Neutral Atmosphere:
All information will be based on results from Steve Bougher's MTGCM,
and published results from Jean-Yves Chaufray for hydrogen densities.
The MTGCM was run for Ls=270, F10.7 = 34 at Mars.
z is referenced from the Martian surface, in km.
Density has units of cm^-3.
CO2 number density:
6.04e18 * exp( -z / 6.98 ) + 1.67e15 * exp( -z / 11.49 )
O number density:
5.85e13 * exp( -z / 10.56 ) + 7.02e09 * exp( -z / 33.97 )
H number density:
1.5e5 * exp[ 25965 * ( 1 / (z+3393.5) - 1/3593.5) ]
Neutral Temperature:
-64.56 * exp[ -0.5 * ( (z-115.7) / 20.14 )^2 ] + 196.95
Following are the description of the extra hot corona components of hydrogen and oxygen.
Neutral Hydrogen (Cold and Hot):
Use the same cold hydrogen profile as before, but now add a hot component.
hot component:
1.9e4 * exp[ 10365 * ( 1 / (z+3393.5) - 1 / 3593.5) ]
Neutral Oxygen (Cold and Hot):
Use the same cold oxygen profile as before, but now add a hot component, supplied by Arnaud Valeille.
hot component:
5.23e3 * exp( -z / 626.2) + 9.76e2 * exp( -z / 2790.) + 3.71e4 * exp( -z / 88.47) ]
'''
R_m = 3376000.
[docs] @staticmethod
def part_func(r, lam=None, Tc=200, partfun=0):
''' Polynomial implementation of Chamberlain's partition function zeta. Obselete due to large error for high altitudes.
:param r: radius from Mars center, in meters
:param zc: exobase altitude, in meters
:param hc: scale height at exobase, in meters
:param partfun: flag to determine what partition function to use. 1 for polynomial, 2 for Chamberlain1963.
Default values (zc=200e3, hc=1176e3) are for the hot hydrogen component: Tc = 500 K, nc = 1.5e4 cm-3 [chaufray2008]
:return: (?? Any description ??)
'''
# values for hydrogen
# h0 = 164e3
# hs = 706e3 # 706 at solar min
# transform to planetocentric distances
if partfun == 0:
zeta = r * 0 + 1
elif partfun == 1:
lam, lam_c, z_bal, z_sat, z_esc = ZetaChamberlain1963.zeta_chamberlain1963(r, Tc=Tc, lam=lam)
zeta = z_bal + z_esc
else:
zeta = 0
return zeta
[docs] @staticmethod
def CO2(z):
""" CO2 density model for solar min at Mars.
:param z: Altitude in km.
:return: CO2 density in /cm3.
"""
n = 6.04e18 * _np.exp(-z / 6.98) + 1.67e15 * _np.exp(-z / 11.49)
return n
[docs] @staticmethod
def O_cold(z):
""" Cold O density model for solar min at Mars.
:param z: Altitude in km.
:return: Cold O density in /cm3.
"""
n = 5.85e13 * _np.exp(-z / 10.56) + 7.02e09 * _np.exp(-z / 33.97)
return n
[docs] @staticmethod
def H_cold(z, part=0):
""" Cold H density model for solar min at Mars.
:param z: Altitude in km.
:keyword part: /?? Any description ??/
:return: Cold H density in /cm3.
"""
zeta = ISSI2008_SolarMin.part_func((z + 3393.5) * 1000, partfun=part, lam=25965 / (z + 3393.5))
n = 1.5e5 * _np.exp(25965 * (1 / (z + 3393.5) - 1 / 3593.5)) * zeta
return n
[docs] @staticmethod
def H_hot(z, part=0):
""" Hot H density model for solar min at Mars.
:param z: Altitude in km.
:keyword part: /?? Any description ??/
:return: Hot H density in /cm3.
"""
if part == 0:
zeta = 1
elif part == 1:
zeta = ISSI2008_SolarMin.part_func((z + 3393.5) * 1000, partfun=part, lam=10365 / (z + 3393.5))
# elif part == 2:
# zeta = ISSI2008_SolarMin.part_func((z + 3393.5) * 1000,
# lam=10365 / (z + 3393.5)) # directly specify lambda in the partition function
else:
zeta = 0
n = 1.9e4 * _np.exp(10365 * (1 / (z + 3393.5) - 1 / 3593.5)) * zeta
return n
[docs] @staticmethod
def O_hot(z):
""" Hot O density model for solar min at Mars.
:param z: Altitude in km.
:return: Hot O density in /cm3.
"""
n = 5.23e3 * _np.exp(-z / 626.2) + 9.76e2 * _np.exp(-z / 2790.) + 3.71e4 * _np.exp(-z / 88.47)
return n
class _EF_const:
"""
.. codeauthor:: X.-D. Wang
"""
G = 6.67408e-11 # gravitational constant, m3 kg-1 s-2
R_m = 3376000. # meter
mp = 1.67e-27 # kg
eV = 1.6e-19 # joule
kb = 1.38e-23 # Boltzmann constatn in SI unit. m**2 kg s**-2 K**-1
Mars_mass = 6.417e23 # kg
from scipy.special import erfc, gamma, gammainc
[docs]class ZetaChamberlain1963:
""" Implementation of Chamberlain's partition function for ballistic, satellite and escape populations.
.. codeauthor:: X.-D. Wang
<?? More detailed comments if this class is intended to be used by users ??>
"""
G = _EF_const.G
Mars_mass = _EF_const.Mars_mass
mp = _EF_const.mp
kb = _EF_const.kb
R_m = _EF_const.R_m
lam_c_tab = _np.array([1, 2, 3, 4, 5, 7.5, 10, 15]) # lam_c = G M m/(kb Tc rc), G gravitational constant, M mass of Mars, m mass of hydrogen atom,
# Tc [144, 725, 483, 362, 290, 193, 145, 97] # kb boltzamnn constant, Tc exobase temerature, rc exobase altitude
lams_full = _np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,
2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0,
8.5, 9.0, 9.5, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0])
beta = G * Mars_mass * mp / kb
Tcs = beta / lam_c_tab / (200e3 + R_m)
#print(lam_c_tab, Tcs)
[docs] @staticmethod
def gam(a, x):
'''
UNNORMALIZED Incomplete Gamma function. = scipy.gamma * scipy.gammainc
:param a:
:param x:
:return:
'''
return gamma(a) * gammainc(a, x)
[docs] @classmethod
def zeta_chamberlain1963_tab(self, lam_c_orig):
'''
Tabulated Chamberlain's partition function. Copied from Chamberlain1963, used to verify the analytical implementation.
:param lam_c_orig: Original lambda value.
:return: zeta values as a function of lambda.
'''
zeta_bal = { # ballistic population
1:_np.array([0.0031541 , 0.0158253, 0.0388161, 0.0711837, 0.1115230, 0.1584466, 0.2108664, 0.2683348, 0.3321255, 0.4275932]),
2:_np.array([0.0016478, 0.0086060, 0.0218891, 0.0414704, 0.0668653, 0.0973691, 0.1321862, 0.1705066, 0.2115528, 0.2546090,
0.2990377, 0.3442906, 0.3899154, 0.4355640, 0.4810061, 0.5261605, 0.5711717, 0.6166366, 0.6645299, 0.7385358]),
3: _np.array([0.0011154, 0.0059134, 0.0152638, 0.0293366, 0.0479653, 0.0707944, 0.0973625, 0.1271557, 0.1596422, 0.1942969,
0.2306169, 0.2681325, 0.3064139, 0.3450743, 0.3837724, 0.4222115, 0.4601389, 0.4973445, 0.5336582, 0.5689490,
0.6031222, 0.6361199, 0.6679210, 0.6985453, 0.7280631, 0.7566176, 0.7844794, 0.8122052, 0.8412825, 0.8883897]),
4: _np.array([0.0008430, 0.0045046, 0.0117200, 0.0227057, 0.0374208, 0.0556711, 0.0771688, 0.1015702, 0.1285021, 0.1575804,
0.1884238, 0.2206628, 0.2539470, 0.2879492, 0.3223688, 0.3569330, 0.3913976, 0.4255469, 0.4591931, 0.4921747,
0.5243559, 0.5556242, 0.5858892, 0.6150808, 0.6431474, 0.6700545, 0.6957824, 0.7203260, 0.7436928, 0.7659024,
0.8611442, 0.9539883]),
5: _np.array([0.0006775, 0.0036380, 0.0095123, 0.0185221, 0.0306831, 0.0458852, 0.0639381, 0.0845997, 0.1075971, 0.1326409,
0.1594359, 0.1876891, 0.2171163, 0.2474455, 0.2784207, 0.3098037, 0.3413754, 0.3729367, 0.4043083, 0.4353311,
0.4658647, 0.4957879, 0.5249966, 0.5534036, 0.5809370, 0.6075394, 0.6331664, 0.6577855, 0.6813752, 0.7039237,
0.8011886, 0.8746162, 0.9289265, 0.9814339]),
7.5: _np.array([0.0004545, 0.0024567, 0.0064675, 0.0126828, 0.0211637, 0.0318873, 0.0447743, 0.0597070, 0.0765414, 0.0951166,
0.1152607, 0.1367970, 0.1595472, 0.1833347, 0.2079873, 0.2333387, 0.2592300, 0.2855109, 0.3120403, 0.3386868,
0.3653290, 0.3918555, 0.4181652, 0.4441668, 0.4697785, 0.4949283, 0.5195526, 0.5435969, 0.5670143, 0.5897658,
0.6926384, 0.7765107, 0.8421618, 0.8918194, 0.9282875, 0.9544006, 0.9727474, 0.9856534, 0.9981837]),
10: _np.array([0.0003420, 0.0018545, 0.0048995, 0.0096434, 0.0161540, 0.0244367, 0.0344551, 0.0461430, 0.0594135, 0.0741651,
0.0902866, 0.1076602, 0.1261649, 0.1456785, 0.1660789, 0.1872463, 0.2090636, 0.2314178, 0.2542002, 0.2773072,
0.3006408, 0.3241087, 0.3476242, 0.3711069, 0.3944822, 0.4176819, 0.4406431, 0.4633092, 0.4856288, 0.5075563,
0.6100652, 0.6986411, 0.7722627, 0.8314909, 0.8778145, 0.9131511, 0.9395052, 0.9587566, 0.9725498, 0.9822550,
0.9889730, 0.9935690, 0.9967392, 0.9998306]),
15: _np.array([0.0002287, 0.0012444, 0.0032996, 0.0065191, 0.0109642, 0.0166556, 0.0235868, 0.0317320, 0.0410510, 0.0514939,
0.0630029, 0.0755155, 0.0889653, 0.1032836, 0.1184008, 0.1342465, 0.1507507, 0.1678442, 0.1854591, 0.2035290,
0.2219893, 0.2407776, 0.2598337, 0.2791000, 0.2985212, 0.3180449, 0.3376212, 0.3572031, 0.3767460, 0.3962084,
0.4910461, 0.5791505, 0.6581367, 0.7268622, 0.7851346, 0.8334304, 0.8726492, 0.9039130, 0.9284158, 0.9473205,
0.9616937, 0.9724720, 0.9804492, 0.9862797, 0.9934941, 0.9970901, 0.9988011, 0.9995809, 0.9999990])}
zeta_esc = { # escaping population
1:_np.array([0.0054313, 0.0182501, 0.0355262, 0.0556409, 0.0777011, 0.1013701, 0.1269113, 0.1555973, 0.1916387, 0.2862033]),
2:_np.array([0.0021737, 0.0071668, 0.0136826, 0.0209618, 0.0285147, 0.0360195, 0.0432664, 0.0501251, 0.0565225, 0.0624283,
0.0678453, 0.0728034, 0.0773571, 0.0815874, 0.0856088, 0.0895902, 0.0938043, 0.0987735, 0.1058662, 0.1307320]),
3:_np.array([0.0013295, 0.0043529, 0.0082590, 0.0125751, 0.0169971, 0.0213247, 0.0254273, 0.0292223, 0.0326615, 0.0357213,
0.0383953, 0.0406901, 0.0426207, 0.0442084, 0.0454782, 0.0464580, 0.0471770, 0.0476651, 0.0479527, 0.0480704,
0.0480491, 0.0479206, 0.0477188, 0.0474811, 0.0472530, 0.0470948, 0.0470996, 0.0474424, 0.0485820, 0.0558051]),
4:_np.array([0.0009523, 0.0031064, 0.0058754, 0.0089191, 0.0120195, 0.0150339, 0.0178698, 0.0204692, 0.0227990, 0.0248434,
0.0265991, 0.0280718, 0.0292729, 0.0302181, 0.0309256, 0.0314152, 0.0317071, 0.0318214, 0.0317782, 0.0315964,
0.0312943, 0.0308891, 0.0303966, 0.0298317, 0.0292081, 0.0285385, 0.0278346, 0.0271070, 0.0263659, 0.0256205,
0.0221587, 0.0230058]),
5:_np.array([0.0007404, 0.0024094, 0.0045483, 0.0068921, 0.0092718, 0.0115769, 0.0137363, 0.0157060, 0.0174609, 0.0189898,
0.0202909, 0.0213694, 0.0222348, 0.0228999, 0.0233795, 0.0236892, 0.0238452, 0.0238633, 0.0237593, 0.0235480,
0.0232437, 0.0228594, 0.0224075, 0.0218992, 0.0213449, 0.0207537, 0.0201343, 0.0194940, 0.0188396, 0.0181771,
0.0149050, 0.0119878, 0.0096776, 0.0092830]),
7.5:_np.array([0.0004745, 0.0015391, 0.0028978, 0.0043806, 0.0058796, 0.0073248, 0.0086716, 0.0098926, 0.0109728, 0.0119056,
0.0126909, 0.0133327, 0.0138376, 0.0142145, 0.0144733, 0.0146245, 0.0146786, 0.0146463, 0.0145377, 0.0143625,
0.0141298, 0.0138484, 0.0135259, 0.0131696, 0.0127860, 0.0123810, 0.0119600, 0.0115276, 0.0110881, 0.0106449,
0.0084732, 0.0065348, 0.0049265, 0.0036538, 0.0026801, 0.0019546, 0.0014270, 0.0010573, 0.0009081]),
10:_np.array([0.0003487, 0.0011292, 0.0021232, 0.0032057, 0.0042978, 0.0053483, 0.0063248, 0.0072076, 0.0379859, 0.0086555,
0.0092163, 0.0096716, 0.0100266, 0.0102880, 0.0104631, 0.0105599, 0.0105862, 0.0105498, 0.0104584, 0.0103191,
0.0101386, 0.0099233, 0.0096789, 0.0094107, 0.0091234, 0.0088214, 0.0085085, 0.0081881, 0.0078633, 0.0075367,
0.0059457, 0.0045380, 0.0033792, 0.0024689, 0.0017773, 0.0012647, 0.0008918, 0.0006246, 0.0004354, 0.0003028,
0.0002106, 0.0001472, 0.0001046, 0.0000846]),
15:_np.array([0.0002277, 0.0007361, 0.0013821, 0.0020842, 0.0027911, 0.0034695, 0.0040986, 0.0046658, 0.0051642, 0.0055914,
0.0059475, 0.0062348, 0.0064569, 0.0066182, 0.0067237, 0.0067785, 0.0067880, 0.0067573, 0.0066912, 0.0065946,
0.0064718, 0.0063270, 0.0061639, 0.0059859, 0.0057961, 0.0055973, 0.0053920, 0.0051823, 0.0049703, 0.0047576,
0.0037273, 0.0028233, 0.0020849, 0.0015092, 0.0010751, 0.0007561, 0.0005260, 0.0003627, 0.0002482, 0.0001687,
0.0001140, 0.0000767, 0.0000513, 0.0000342, 0.0000150, 0.0000065, 0.0000027, 0.0000011, 0.0000005])}
if lam_c_orig not in self.lam_c_tab: # if input lam_c is not in the given lam_c list.
# pick the nearest one
diff = _np.abs(self.lam_c_tab - lam_c_orig)
minind = _np.argmin(diff)
lam_c = self.lam_c_tab[minind]
print('No matching lambda_c found, using the closest one: \lambda_c = {:}'.format(lam_c))
else:
lam_c = lam_c_orig
picflag = (self.lams_full <= lam_c)
currlams = self.lams_full[picflag]
Tc = self.beta / lam_c / (200e3 + self.R_m)
curralt = (self.beta / Tc / currlams - self.R_m) / 1000 # altitude of different lambda values for given Tc, in km
return curralt, Tc, zeta_bal[lam_c], zeta_esc[lam_c]
[docs] @classmethod
def zeta_bal_shen1963(self, r, T):
'''
Analytical format for the ballistic orbits from Shen 1963, cited by Chamberlain 1963.
:param r:
:param T:
:return:
'''
R_m = self.R_m
G = self.G
Mars_mass = self.Mars_mass
mp = self.mp
kb = self.kb
rc = R_m + 200e3 # exobase radius
alpha = rc/r
E = G * Mars_mass * mp/(rc * T * kb)
H_bottom = _np.pi ** 0.5 * _np.exp(E) / alpha / E ** 0.5 * erfc(E ** 0.5)
H_E = 3/(2 * E) + 1/(1 + H_bottom)
rho_t = _np.exp(-(1 - alpha) * E) - (1 - alpha ** 2) ** 0.5 * _np.exp(-E / (1 + alpha))
term2_top = alpha ** 2 * (1+E) * _np.exp(-E)
term2_bottom = (2 * _np.pi ** 0.5) * E ** 0.5 * (H_E - 1 + alpha * (1 - 2 / 3 * alpha * H_E)) ** 0.5
rho_s = term2_top / term2_bottom
zeta_bal = rho_t - rho_s
return zeta_bal, rho_t, rho_s, alpha
[docs] @classmethod
def zeta_chamberlain1963(self, r, Tc=200, planetmass=6.417e23, rc=3597000., particlemass=1.67e-27, lam=None):
'''
Partition function from Chamberlain 1963, analytical, with 3 types of orbits.
:param r: Radius of the point of interest, in meter.
:param Tc: Exobase temperature, in Kelvin. If lam kwarg is set,
At least ONE of r and Tc must be scalar.
:param planetmass: Mass of the planet, in kg. By default Mars.
:param rc: Radius of the exobase, in meter. By default the Mars exobase at R_m + 200 km.
:param particlemass: Mass of the interested particles, in kg. By default hydrogen atom, 1 amu.
:param lam: Input lambda values instead of radius. If lam is input, it must have the same dimension as r.
:return: All returned values have same dimension as the array in r or Tc.
lam, lambda values for r with exobase temperature Tc: lambda = G * M * m/(kb * Tc * r)
lam_c, lambda value for exobase with temperature Tc: lambda_c = G * M * m/(kb * Tc * rc)
z_bal, Ballistic orbit zeta function value(s) at radius r with exobase tmperature Tc. Range [0, 1]
z_sat, Satellite orbit zeta function.
z_esc, Escape orbit zeta function.
'''
if lam is None:
lam = self.G * planetmass * particlemass/(self.kb * Tc * r)
lam_c = lam * r/rc
psi1 = lam**2 / (lam + lam_c)
z_bal = 2 / _np.pi ** 0.5 * (
self.gam(1.5, lam) - (lam_c**2 - lam**2) ** 0.5 / lam_c * _np.exp(-psi1) * self.gam(1.5, lam - psi1)
)
z_sat = 2 / _np.pi ** 0.5 * self.gam(1.5, lam) - z_bal
z_esc = (
gamma(1.5) - self.gam(1.5, lam) -
(lam_c**2 - lam**2) ** 0.5 / lam_c * _np.exp(-psi1) * (gamma(1.5) - self.gam(1.5, lam - psi1))
) \
/ _np.pi ** 0.5
return lam, lam_c, z_bal, z_sat, z_esc