spacepy.radbelt.RBmodel¶
- class spacepy.radbelt.RBmodel(grid='L', NL=91, const_kp=False)[source]¶
1-D Radial diffusion class
This module contains a class for performing and visualizing 1-D radial diffusion simulations of the radiation belts.
Here is an example using the default settings of the model. Each instance must be initialized with (assuming import radbelt as rb):
>>> rmod = rb.RBmodel()
Next, set the start time, end time, and the size of the timestep:
>>> import datetime >>> start = datetime.datetime(2003,10,14) >>> end = datetime.datetime(2003,12,26) >>> delta = datetime.timedelta(hours=1) >>> rmod.setup_ticks(start, end, delta, dtype='UTC')
Now, run the model over the entire time range using the evolve method:
>>> rmod.evolve()
Finally, visualize the results:
>>> rmod.plot_summary()
Gaussian source term added to radiation belt model.
add_Lmax
(Lmax_model)add last closed drift shell Lmax
add_Lpp
(Lpp_model)add last closed drift shell Lmax
add_PSD_obs
([time, PSD, Lstar, satlist])add PSD observations
add_PSD_twin
([dt, Lt])add observations from PSD database using the ticks list the arguments are the following:
add_omni
([keylist])add omni data to instance according to the tickrange in ticks
add_source
([source, A, mu, sigma])add source parameters A, mu, and sigma for the Gaussian source function
assimilate
([method, inflation])Assimilates data for the radiation belt model using the Ensemble Kalman Filter.
evolve
()calculate the diffusion in L at constant mu,K coordinates
get_DLL
(Lgrid, params[, DLL_model])Calculate DLL as a simple power law function (alpha*L**Bbta) using alpha/beta values from popular models found in the literature and chosen with the kwarg "DLL_model".
plot
([Lmax, Lpp, Kp, Dst, clims, title, values])Create a summary plot of the RadBelt object distribution function.
plot_obs
([Lmax, Lpp, Kp, Dst, clims, title, ...])Create a summary plot of the observations.
set_lgrid
([NL])Using NL grid points, create grid in L.
setup_ticks
(start, end, delta[, dtype])Add time information to the simulation by specifying a start and end time, timestep, and time type (optional).
- Gaussian_source()[source]¶
Gaussian source term added to radiation belt model. The source term is given by the equation:
S = A exp{-(L-mu)^2/(2*sigma^2)}
with A=10^(-8), mu=5.0, and sigma=0.5 as default values
- add_PSD_obs(time=None, PSD=None, Lstar=None, satlist=None)[source]¶
add PSD observations
- Parameters
- timeTicktock datetime array
array of observation times
- PSDlist of numpy arrays
PSD observational data for each time. Each entry in the list is a numpy array with the observations for the corresponding time
- Lstarlist of numpy arrays
Lstar location of each PSD observations. Each entry in the list is a numpy array with the location of the observations for the corresponding time
- satlistlist of satellite names
- Returns
- outlist of dicts
Information of the observational data, where each entry contains the observations and locations of observations for each time specified in the time array. Each list entry is a dictionary with the following information:
- TicksTicktock array
time of observations
- Lstarnumpy array
location of observations
- PSDnumpy array
PSD observation values
- satlist of strings
satellite names
- MUscalar value
Mu value for the observations
- Kscalar value
K value for the observations
- add_PSD_twin(dt=0, Lt=1)[source]¶
add observations from PSD database using the ticks list the arguments are the following:
dt = observation time delta in seconds Lt = observation space delta
- add_source(source=True, A=1e-08, mu=5.0, sigma=0.5)[source]¶
add source parameters A, mu, and sigma for the Gaussian source function
- assimilate(method='EnKF', inflation=0)[source]¶
Assimilates data for the radiation belt model using the Ensemble Kalman Filter. The algorithm used is the SVD method presented by Evensen in 2003 (Evensen, G., Ocean dynamics, 53, pp.343–367, 2003). To compensate for model errors, three inflation algorithms are implemented. The inflation methodology is specified by the ‘inflation’ argument, and the options are the following:
inflation == 0: Add model error (perturbation for the ensemble) around model state values only where observations are available (DEFAULT).
inflation == 1: Add model error (perturbation for the ensemble) around observation values only where observations are available.
inflation == 2: Inflate around ensemble average for EnKF.
Prior to assimilation, a set of data values has to be speficied by setting the start and end dates, and time step, using the setup_ticks funcion of the radiation belt model:
>>> import spacepy >>> import datetime >>> from spacepy import radbelt
>>> start = datetime.datetime(2002,10,23) >>> end = datetime.datetime(2002,11,4) >>> delta = datetime.timedelta(hours=0.5) >>> rmod.setup_ticks(start, end, delta, dtype='UTC')
Once the dates and time step are specified, the data is added using the add_PSD function:
>>> rmod.add_PSD()
The observations are averaged over the time windows, whose interval is give by the time step.
Once the dates and data are set, the assimiation is performed using the ‘assimilate’ function:
>>> rmod.assimilate(inflation=1)
This function will add the PSDa values, which are the analysis state of the radiation belt using the observations within the dates. To plot the analysis simply use the plot funtion:
>>> rmod.plot(values=rmod.PSDa,clims=[-10,-6],Lmax=False,Kp=False,Dst=False)
- get_DLL(Lgrid, params, DLL_model='BA2000')[source]¶
Calculate DLL as a simple power law function (alpha*L**Bbta) using alpha/beta values from popular models found in the literature and chosen with the kwarg “DLL_model”.
The calculated DLL is returned, as is the alpha and beta values used in the calculationp.
The output DLL is in units of units/day.
- plot(Lmax=True, Lpp=False, Kp=True, Dst=True, clims=[0, 10], title=None, values=None)[source]¶
Create a summary plot of the RadBelt object distribution function. For reference, the last closed drift shell, Dst, and Kp are all included. These can be disabled individually using the corresponding Boolean kwargs.
The clims kwarg can be used to manually set the color bar range. To use, set it equal to a two-element list containing minimum and maximum Log_10 value to plot. Default action is to use [0,10] as the log_10 of the color range. This is good enough for most applications.
The title of the top most plot defaults to ‘Summary Plot’ but can be customized using the title kwarg.
The figure object and all three axis objects (PSD axis, Dst axis, and Kp axis) are all returned to allow the user to further customize the plots as necessary. If any of the plots are excluded, None is returned in their stead.
Examples
>>> rb.plot(Lmax=False, Kp=False, clims=[2,10], title='Good work!')
This command would create the summary plot with a color bar range of 100 to 10^10. The Lmax line and Kp values would be excluded. The title of the topmost plot (phase space density) would be set to ‘Good work!’.
- plot_obs(Lmax=True, Lpp=False, Kp=True, Dst=True, clims=[0, 10], title=None, values=None)[source]¶
Create a summary plot of the observations. For reference, the last closed drift shell, Dst, and Kp are all included. These can be disabled individually using the corresponding boolean kwargs.
The clims kwarg can be used to manually set the color bar range. To use, set it equal to a two-element list containing minimum and maximum Log_10 value to plot. Default action is to use [0,10] as the log_10 of the color range. This is good enough for most applications.
The title of the top most plot defaults to ‘Summary Plot’ but can be customized using the title kwarg.
The figure object and all three axis objects (PSD axis, Dst axis, and Kp axis) are all returned to allow the user to further customize the plots as necessary. If any of the plots are excluded, None is returned in their stead.
Examples
>>> rb.plot_obs(Lmax=False, Kp=False, clims=[2,10], title='Observations Plot')
This command would create the summary plot with a color bar range of 100 to 10^10. The Lmax line and Kp values would be excluded. The title of the topmost plot (phase space density) would be set to ‘Good work!’.
- set_lgrid(NL=91)[source]¶
Using NL grid points, create grid in L. Default number of points is 91 (dL=0.1).
- setup_ticks(start, end, delta, dtype='ISO')[source]¶
Add time information to the simulation by specifying a start and end time, timestep, and time type (optional).
Examples
>>> start = datetime.datetime(2003,10,14) >>> end = datetime.datetime(2003,12,26) >>> delta = datetime.timedelta(hours=1) >>> rmod.setup_ticks(start, end, delta, dtype='UTC')