spacepy.empiricals¶
Module with some useful empirical models (plasmapause, magnetopause, Lmax)
Authors: Steve Morley, Josef Koller Institution: Los Alamos National Laboratory Contact: smorley@lanl.gov
Copyright 2010 Los Alamos National Security, LLC.
Functions
|
Calculates the Shue et al. (1997) subsolar magnetopause radius. |
|
Calculate the pressure-corrected Dst index, Dst* |
|
Return the expected solar wind temperature based on the bulk velocity |
|
calculate a simple empirical model for Lmax - last closed drift-shell |
|
Calculates the Shue et al. (1997) subsolar magnetopause radius. |
|
Calculates the Shue et al. (1997) position in equatorial plane. |
|
Plasmapause location model(s) |
|
Returns a SpaceData with energy and fluence spectra of solar particle events |
|
Calculates solar rotation number (Carrington or Bartels) for a given date/time |
Empirical lookup of power for sin^n pitch angle model from Vampola (1996) |
|
|
calculate a simple empirical model for Lmax - last closed drift-shell |
|
Plasmapause location model(s) |
|
Calculate omnidirectional flux [(s cm^2 kev)^-1] from directional flux [(s sr cm^2 keV)^-1] array |
|
Pitch angle model of sin^n form |
- spacepy.empiricals.ShueMP(ticks, dbase='QDhourly', alpha=[])¶
Calculates the Shue et al. (1997) subsolar magnetopause radius
Shue et al. (1997), A new functional form to study the solar wind control of the magnetopause size and shape, J. Geophys. Res., 102(A5), 9497–9511, doi:10.1029/97JA00196.
- Parameters:
- ticksspacepy.time.Ticktock
TickTock object of desired times (will be interpolated from hourly OMNI data) OR dictionary of form {‘P’: [], ‘Bz’: []} Where P is SW ram pressure [nPa] and Bz is IMF Bz (GSM) [nT]
- alphalist
Used as an optional return value to obtain the flaring angles. To use, assign an empty list and pass to this function through the keyword argument. The list will be modified in place, adding the flaring angles for each time step.
- Returns:
- outfloat
Magnetopause (sub-solar point) standoff distance [Re]
Examples
>>> import spacepy.time as spt >>> import spacepy.empiricals as emp >>> ticks = spt.tickrange('2002-01-01T12:00:00','2002-01-04T00:00:00',.25) >>> emp.getMPstandoff(ticks) array([ 10.57319537, 10.91327764, 10.75086873, 10.77577207, 9.78180261, 11.0374474 , 11.4065 , 11.27555451, 11.47988573, 11.8202582 , 11.23834814]) >>> data = {'P': [2,4], 'Bz': [-2.4, -2.4]} >>> emp.getMPstandoff(data) array([ 9.96096838, 8.96790412])
- spacepy.empiricals.getDststar(ticks, model='OBrien', dbase='QDhourly')[source]¶
Calculate the pressure-corrected Dst index, Dst*
We need to add in the references to the models here!
- Parameters:
- ticksspacepy.time.Ticktock
TickTock object of desired times (will be interpolated from hourly OMNI data) OR dictionary including ‘Pdyn’ and ‘Dst’ keys where data are lists or arrays and Dst is in [nT], and Pdyn is in [nPa]
- Returns:
- outfloat
Dst* - the pressure corrected Dst index from OMNI [nT]
Examples
Coefficients are applied to the standard formulation e.g. Burton et al., 1975 of Dst* = Dst - b*sqrt(Pdyn) + c The default is the O’Brien and McPherron model (2002). Other options are Burton et al. (1975) and Borovsky and Denton (2010)
>>> import spacepy.time as spt >>> import spacepy.omni as om >>> import spacepy.empiricals as emp >>> ticks = spt.tickrange('2000-10-16T00:00:00', '2000-10-31T12:00:00', 1/24.) >>> dststar = emp.getDststar(ticks) >>> dststar[0] -21.317220132108943
User-determined coefficients can also be supplied as a two-element list or tuple of the form (b,c), e.g.
>>> dststar = emp.getDststar(ticks, model=(2,11)) #b is extreme driving from O'Brien
We have chosen the OBrien model as the default here as this was rigorously determined from a very long data set and is pertinent to most conditions. It is, however, the most conservative correction. Additionally, Siscoe, McPherron and Jordanova (2005) argue that the pressure contribution to Dst diminishes during magnetic storms.
To show the relative differences, run the following example:
>>> import matplotlib.pyplot as plt >>> params = [('Burton','k-'), ('OBrien','r-'), ('Borovsky','b-')] >>> for model, col in params: dststar = getDststar(ticks, model=model) plt.plot(ticks.UTC, dststar, col)
- spacepy.empiricals.getExpectedSWTemp(velo, model='XB15', units='K')[source]¶
Return the expected solar wind temperature based on the bulk velocity
The formulations used by this function are those given by, L87 – Lopez, R.E., J. Geophys. Res., 92, 11189-11194, 1987 BS06 – Borovsky, J.E. and J.T. Steinberg, Geophysical Monograph Series 167, 59-76, 2006 XB15 – Xu, F. and J.E. Borovsky, J. Geophys. Res., 120, 70-100, 2015
- Parameters:
- veloarray-like
Array like of solar wind bulk velocity values [km/s]
- modelstr [optional]
Name of model to use. Valid choices are L87, BS06 and XB15. Default is XB15
- unitsstr [optional]
Units for output temperature, options are eV or K. Default is Kelvin [K]
- Returns:
- Texparray-like
The expected solar wind temperature given the bulk velocity [K] or [eV]
- spacepy.empiricals.getLmax(ticks, model='JKemp', dbase='QDhourly')[source]¶
calculate a simple empirical model for Lmax - last closed drift-shell
Uses the parametrized Lmax from: Koller and Morley (2010) ‘Magnetopause shadowing effects for radiation belt models during high-speed solar wind streams’ American Geophysical Union, Fall Meeting 2010, abstract #SM13A-1787
- Parameters:
- ticksspacepy.time.Ticktock
Ticktock object of desired times
- modelstring, optional
‘JKemp’ (default - empirical model of J. Koller)
- Returns:
- outnp.ndarray
Lmax - L* of last closed drift shell
See also
Examples
>>> from spacepy.empiricals import getLmax >>> import spacepy.time as st >>> import datetime >>> ticks = st.tickrange(datetime.datetime(2000, 1, 1), datetime.datetime(2000, 1, 3), deltadays=1) array([ 7.4928412, 8.3585632, 8.6463423])
- spacepy.empiricals.getMPstandoff(ticks, dbase='QDhourly', alpha=[])[source]¶
Calculates the Shue et al. (1997) subsolar magnetopause radius
Shue et al. (1997), A new functional form to study the solar wind control of the magnetopause size and shape, J. Geophys. Res., 102(A5), 9497–9511, doi:10.1029/97JA00196.
- Parameters:
- ticksspacepy.time.Ticktock
TickTock object of desired times (will be interpolated from hourly OMNI data) OR dictionary of form {‘P’: [], ‘Bz’: []} Where P is SW ram pressure [nPa] and Bz is IMF Bz (GSM) [nT]
- alphalist
Used as an optional return value to obtain the flaring angles. To use, assign an empty list and pass to this function through the keyword argument. The list will be modified in place, adding the flaring angles for each time step.
- Returns:
- outfloat
Magnetopause (sub-solar point) standoff distance [Re]
Examples
>>> import spacepy.time as spt >>> import spacepy.empiricals as emp >>> ticks = spt.tickrange('2002-01-01T12:00:00','2002-01-04T00:00:00',.25) >>> emp.getMPstandoff(ticks) array([ 10.57319537, 10.91327764, 10.75086873, 10.77577207, 9.78180261, 11.0374474 , 11.4065 , 11.27555451, 11.47988573, 11.8202582 , 11.23834814]) >>> data = {'P': [2,4], 'Bz': [-2.4, -2.4]} >>> emp.getMPstandoff(data) array([ 9.96096838, 8.96790412])
- spacepy.empiricals.getMagnetopause(ticks, LTs=None, dbase='QDhourly')[source]¶
Calculates the Shue et al. (1997) position in equatorial plane
Shue et al. (1997), A new functional form to study the solar wind control of the magnetopause size and shape, J. Geophys. Res., 102(A5), 9497–9511, doi:10.1029/97JA00196.
- Parameters:
- ticksspacepy.time.Ticktock
TickTock object of desired times (will be interpolated from hourly OMNI data) OR dictionary of form {‘P’: [], ‘Bz’: []} Where P is SW ram pressure [nPa] and Bz is IMF Bz (GSM) [nT]
- LTsarray-like
Array-like of local times for evaluating the magnetopause model. Default is 6 LT to 18 LT in steps of 20 minutes.
- Returns:
- outarray
NxMx2 array of magnetopause positions [Re] N is number of timesteps, M is number of local times. The 2 positions are the X_GSE and Y_GSE positions of the magnetopause
Examples
>>> import spacepy.time as spt >>> import spacepy.empiricals as emp >>> ticks = spt.Ticktock(['2002-01-01T12:00:00','2002-01-04T00:00:00']) >>> localtimes = [13,12,11] >>> emp.getMagnetopause(ticks, LTs=localtimes) array([[[ 10.27674331, -2.75364507], [ 10.52909163, 0. ], [ 10.27674331, 2.75364507]], [[ 10.91791834, -2.9254474 ], [ 11.18712131, 0. ], [ 10.91791834, 2.9254474 ]]]) >>> emp.getMPstandoff(ticks) #should give same result as getMagnetopause for 12LT array([ 10.52909163, 11.18712131])
To plot the magnetopause: >>> import numpy as np >>> import spacepy.plot as splot >>> import matplotlib.pyplot as plt >>> localtimes = np.arange(5, 19.1, 0.5) >>> mp_pos = emp.getMagnetopause(ticks, localtimes) >>> plt.plot(mp_pos[0,:,0], mp_pos[0,:,1]) >>> ax1 = plt.gca() >>> ax1.set_xlim([-5,20]) >>> ax1.set_xlabel(‘X$_{GSE}$ [R$_E$]’) >>> ax1.set_ylabel(‘Y$_{GSE}$ [R$_E$]’) >>> splot.dual_half_circle(ax=ax1) >>> ax1.axes.set_aspect(‘equal’)
- spacepy.empiricals.getPlasmaPause(ticks, model='M2002', LT='all', omnivals=None)[source]¶
Plasmapause location model(s)
CA1992 – Carpenter, D. L., and R. R. Anderson, An ISEE/whistler model of equatorial electron density in the magnetosphere, J. Geophys. Res., 97, 1097, 1992. M2002 – Moldwin, M. B., L. Downward, H. K. Rassoul, R. Amin, and R. R. Anderson, A new model of the location of the plasmapause: CRRES results, J. Geophys. Res., 107(A11), 1339, doi:10.1029/2001JA009211, 2002. RT1970 – Rycroft, M. J., and J. O. Thomas, The magnetospheric plasmapause and the electron density trough at the alouette i orbit, Planetary and Space Science, 18(1), 65-80, 1970
- Parameters:
- ticksspacepy.time.Ticktock
TickTock object of desired times
- Lpp_modelstring, optional
‘CA1992’ or ‘M2002’ (default) CA1992 returns the Carpenter and Anderson model, M2002 returns the Moldwin et al. model
- LTint, float
requested local time sector, ‘all’ is valid option
- omnivalsspacepy.datamodel.SpaceData, dict
dict-like containing UTC (datetimes) and Kp keys
- Returns:
- outfloat
Plasmapause radius in Earth radii
- Warns:
- RuntimeWarning
If the CA1992 model is called with LT as it is not implemented
Examples
>>> import spacepy.time as spt >>> import spacepy.empiricals as emp >>> ticks = spt.tickrange('2002-01-01T12:00:00','2002-01-04T00:00:00',.25) >>> emp.getPlasmaPause(ticks) array([ 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.26859998, 5.772 , 5.6574 , 5.6574 ])
- spacepy.empiricals.getSolarProtonSpectra(norm=32000000.0, gamma=-0.96, E0=15.0, Emin=0.1, Emax=600, nsteps=100)[source]¶
Returns a SpaceData with energy and fluence spectra of solar particle events
The formulation follows that of: Ellison and Ramaty ApJ 298: 400-408, 1985 dJ/dE = K^{-gamma}exp(-E/E0)
and the defualt values are the 10/16/2003 SEP event of: Mewaldt, R. A., et al. (2005), J. Geophys. Res., 110, A09S18, doi:10.1029/2005JA011038.
- Returns:
- datadm.SpaceData
SpaceData with the energy and fluence values
- Other Parameters:
- normfloat
Normilization factor for the intensity of the SEP event
- gammafloat
Power law index
- E0float
Expoential scaling factor
- Eminfloat
Minimum energy for fit
- Emaxfloat
Maximum energy for fit
- nstepsint
The number of log spaced energy steps to return
- spacepy.empiricals.getSolarRotation(ticks, rtype='carrington', fp=False, reverse=False)[source]¶
Calculates solar rotation number (Carrington or Bartels) for a given date/time
- Parameters:
- ticksspacepy.time.Ticktock or datetime.datetime
- Returns:
- rnumberinteger or array
Carrington (or Bartels) rotation number
- spacepy.empiricals.getVampolaOrder(L)[source]¶
Empirical lookup of power for sin^n pitch angle model from Vampola (1996)
Vampola, A.L. Outer zone energetic electron environment update, Final Report of ESA/ESTEC/WMA/P.O. 151351, ESA-ESTEC, Noordwijk, The Netherlands, 1996.
- Parameters:
- Larraylike or float
- Returns:
- orderarray
coefficient for sin^n model corresponding to McIlwain L (computed for OP77?)
Examples
Apply Vampola pitch angle model at L=[4, 6.6]
>>> from spacepy.empiricals import vampolaPA, getVampolaOrder >>> order = getVampolaOrder([4,6.6]) >>> order array([ 3.095 , 1.6402]) >>> vampolaPA([3000, 3000], alpha=[45, 90], order=order) (array([[ 140.08798878, 192.33572182], [ 409.49143136, 339.57417256]]), [45, 90])
- spacepy.empiricals.get_Lmax(ticks, model='JKemp', dbase='QDhourly')¶
calculate a simple empirical model for Lmax - last closed drift-shell
Uses the parametrized Lmax from: Koller and Morley (2010) ‘Magnetopause shadowing effects for radiation belt models during high-speed solar wind streams’ American Geophysical Union, Fall Meeting 2010, abstract #SM13A-1787
- Parameters:
- ticksspacepy.time.Ticktock
Ticktock object of desired times
- modelstring, optional
‘JKemp’ (default - empirical model of J. Koller)
- Returns:
- outnp.ndarray
Lmax - L* of last closed drift shell
See also
Examples
>>> from spacepy.empiricals import getLmax >>> import spacepy.time as st >>> import datetime >>> ticks = st.tickrange(datetime.datetime(2000, 1, 1), datetime.datetime(2000, 1, 3), deltadays=1) array([ 7.4928412, 8.3585632, 8.6463423])
- spacepy.empiricals.get_plasma_pause(ticks, model='M2002', LT='all', omnivals=None)¶
Plasmapause location model(s)
CA1992 – Carpenter, D. L., and R. R. Anderson, An ISEE/whistler model of equatorial electron density in the magnetosphere, J. Geophys. Res., 97, 1097, 1992. M2002 – Moldwin, M. B., L. Downward, H. K. Rassoul, R. Amin, and R. R. Anderson, A new model of the location of the plasmapause: CRRES results, J. Geophys. Res., 107(A11), 1339, doi:10.1029/2001JA009211, 2002. RT1970 – Rycroft, M. J., and J. O. Thomas, The magnetospheric plasmapause and the electron density trough at the alouette i orbit, Planetary and Space Science, 18(1), 65-80, 1970
- Parameters:
- ticksspacepy.time.Ticktock
TickTock object of desired times
- Lpp_modelstring, optional
‘CA1992’ or ‘M2002’ (default) CA1992 returns the Carpenter and Anderson model, M2002 returns the Moldwin et al. model
- LTint, float
requested local time sector, ‘all’ is valid option
- omnivalsspacepy.datamodel.SpaceData, dict
dict-like containing UTC (datetimes) and Kp keys
- Returns:
- outfloat
Plasmapause radius in Earth radii
- Warns:
- RuntimeWarning
If the CA1992 model is called with LT as it is not implemented
Examples
>>> import spacepy.time as spt >>> import spacepy.empiricals as emp >>> ticks = spt.tickrange('2002-01-01T12:00:00','2002-01-04T00:00:00',.25) >>> emp.getPlasmaPause(ticks) array([ 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.42140002, 6.26859998, 5.772 , 5.6574 , 5.6574 ])
- spacepy.empiricals.omniFromDirectionalFlux(fluxarr, alphas, norm=True)[source]¶
Calculate omnidirectional flux [(s cm^2 kev)^-1] from directional flux [(s sr cm^2 keV)^-1] array
J = 2.pi integ(j sin(a) da) If kwarg ‘norm’ is True (default), the omni flux is normalized by 4.pi to make it per steradian, in line with the PRBEM guidelines
- Parameters:
- fluxarrarraylike
Array of directional flux values
- alphasarraylike
Array of pitch angles corresponding to fluxarr
- Returns:
- omnifluxfloat
Omnidirectional flux value
Examples
Roundtrip from omni flux, to directional flux (Vampola model), integrate to get back to omni flux.
>>> from spacepy.empiricals import vampolaPA, omniFromDirectionalFlux >>> dir_flux, pa = vampolaPA(3000, alpha=range(0,181,2), order=4) >>> dir_flux[:10], pa[:10] (array([ 0.00000000e+00, 6.64032473e-04, 1.05986545e-02, 5.34380898e-02, 1.67932162e-01, 4.06999226e-01, 8.36427502e-01, 1.53325140e+00, 2.58383611e+00, 4.08170975e+00]), [0, 2, 4, 6, 8, 10, 12, 14, 16, 18]) >>> omniFromDirectionalFlux(dir_flux, pa, norm=False) 3000.0000008112293
Calculate “spin-averaged” flux, giving answer per steradian
>>> omniFromDirectionalFlux(dir_flux, pa, norm=True) 238.73241470239859
- spacepy.empiricals.vampolaPA(omniflux, **kwargs)[source]¶
Pitch angle model of sin^n form
- Parameters:
- omnifluxarraylike or float
omnidirectional number flux data
- orderinteger or float (optional)
order of sin^n functional form for distribution (default=2)
- alphasarraylike (optional)
pitch angles at which to evaluate the differential number flux (default is 5 to 90 degrees in 36 steps)
- Returns:
- dnfluxarray
differential number flux corresponding to pitch angles alphas
- alphasarray
pitch angles at which the differential number flux was evaluated
Notes
Directional number flux integrated over pitch angle from 0 to 90 degrees is a factor of 4*pi lower than omnidirectional number flux.
Examples
Omnidirectional number flux of [3000, 6000]
>>> from spacepy.empiricals import vampolaPA >>> vampolaPA(3000, alpha=[45, 90]) (array([ 179.04931098, 358.09862196]), [45, 90]) >>> data, pas = vampolaPA([3000, 6000], alpha=[45, 90]) >>> pas [45, 90] >>> data array([[ 179.04931098, 358.09862196], [ 358.09862196, 716.19724391]])
