API Reference

To learn how to quickly get started with pyAMPS, we recommend that people start looking at Usage, and use this page once a comprehensive listing of all functionality is needed. Documentation for each function and its arguments can be found using the help command.

Core functionality (pyamps)

class pyamps.AMPS(v, By, Bz, tilt, f107, minlat=60, maxlat=89.99, height=110.0, dr=2, M0=4, resolution=100, coeff_fn='/home/docs/checkouts/readthedocs.org/user_builds/pyamps/checkouts/latest/pyamps/coefficients/SW_OPER_MIO_SHA_2E_00000000T000000_99999999T999999_0103.txt')[source]

Calculate and plot maps of the model Average Magnetic field and Polar current System (AMPS)

Parameters:
  • v (float) – solar wind velocity in km/s
  • By (float) – IMF GSM y component in nT
  • Bz (float) – IMF GSM z component in nT
  • tilt (float) – dipole tilt angle in degrees
  • f107 (float) – F10.7 index in s.f.u.
  • minlat (float, optional) – low latitude boundary of grids (default 60)
  • maxlat (float, optional) – low latitude boundary of grids (default 89.99)
  • height (float, optional) – altitude of the ionospheric currents in km (default 110)
  • dr (int, optional) – latitudinal spacing between equal area grid points (default 2 degrees)
  • M0 (int, optional) – number of grid points in the most poleward circle of equal area grid points (default 4)
  • resolution (int, optional) – resolution in both directions of the scalar field grids (default 100)
  • coeff_fn (str, optional) – file name of model coefficients - must be in format produced by model_vector_to_txt.py (default is latest version)

Examples

>>> # initialize by supplying a set of external conditions:
>>> m = AMPS(solar_wind_velocity_in_km_per_s,
             IMF_By_in_nT, IMF_Bz_in_nT,
             dipole_tilt_in_deg,
             F107_index)

>>> # make summary plot:
>>> m.plot_currents()

>>> # calculate field-aligned currents on a pre-defined grid
>>> Ju = m.get_upward_current()

>>> # Ju will be evaluated at the following coords:
>>> mlat, mlt = m.scalargrid

>>> # It is also possible to specify coordinates (can be slow with
>>> # repeated calls)
>>> Ju = m.get_upward_current(mlat = my_mlats, mlt = my_mlts)

>>> # get components of the total height-integrated horizontal current,
>>> # calculated on a pre-defined grid
>>> j_east, j_north = m.get_total_current()

>>> # j_east, and j_north will be evaluated at the following coords
>>> # (default grids are different with vector quantities)
>>> mlat, mlt = m.vectorgrid

>>> # update model vectors (tor_c, tor_s, etc.) without
>>> # recalculating the other matrices:
>>> m.update_model(new_v, new_By, new_Bz, new_tilt, new_f107)
tor_c

vector of cos term coefficents in the toroidal field expansion

Type:numpy.ndarray
tor_s

vector of sin term coefficents in the toroidal field expansion

Type:numpy.ndarray
pol_c

vector of cos term coefficents in the poloidal field expansion

Type:numpy.ndarray
pol_s

vector of sin term coefficents in the poloidal field expansion

Type:numpy.ndarray
keys_P

list of spherical harmonic wave number pairs (n,m) corresponding to elements of pol_c and pol_s

Type:list
keys_T

list of spherical harmonic wave number pairs (n,m) corresponding to elements of tor_c and tor_s

Type:list
vectorgrid

grid used to calculate and plot vector fields

Type:tuple
scalargrid

grid used to calculate and plot scalar fields

The grid formats are as follows (see also example below): (np.hstack((mlat_north, mlat_south)), np.hstack((mlt_north, mlt_south)))

The grids can be changed directly, but member function calculate_matrices() must then be called for the change to take effect.

Type:tuple
calculate_matrices()[source]

Calculate the matrices that are needed to calculate currents and potentials

Call this function if and only if the grid has been changed manually

get_AE_indices()[source]

Calculate model synthetic auroral electrojet (AE) indices: AL and AU. The unit is nT

Note

Here, AL and AU are defined as the lower/upper envelope curves for the northward component of the ground magnetic field perturbation that is equivalent with the divergence-free current, evaluated on scalargrid. Thus all the caveats for the get_ground_perturbation() function applies to these calculations as well. An additional caveat is that we have in principle perfect coverage with the model, while the true AE indices are derived using a small set of magnetometers in the auroral zone. The model values are also based on QD northward component, instead of the “H component”, which is used in the official measured AL index. It is possible to calculate model AE indices that are more directly comparable to the measured indices.

Returns:
  • AL_n (float) – Model AL index in the northerm hemisphere
  • AL_s (float) – Model AL index in the southern hemisphere
  • AU_n (float) – Model AU index in the northerm hemisphere
  • AU_s (float) – Model AU index in the southern hemisphere
get_curl_free_current(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the curl-free part of the horizontal current, in units of mA/m. The calculations refer to the height chosen upon initialization of the AMPS object (default 110 km).

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the current. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the vectorgrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the current. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the vectorgrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

  • jcf_eastward (numpy.ndarray, float) – eastward component of the curl-free current evalulated at the coordinates given by the vectorgrid attribute
  • jcf_northward (numpy.ndarray, float) – northward component of the curl-free current evalulated at the coordinates given by the vectorgrid attribute

See also

get_divergence_free_current()
Calculate divergence-free part of the horizontal current
get_total_current()
Calculate total horizontal current
get_curl_free_current_potential(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the curl-free current potential (unit is kA). The curl-free current potential is a scalar alpha which relates to the curl-free part of the horizontal current by J_{cf} = grad(alpha). The calculations refer to the height chosen upon initialization of the AMPS object (default 110 km).

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the current. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the current. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

alpha – Curl-free current potential evaulated at self.scalargrid, or, if specified, mlat/mlt
Return type:

numpy.ndarray
get_divergence_free_current(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the divergence-free part of the horizontal current, in units of mA/m. The calculations refer to the height chosen upon initialization of the AMPS object (default 110 km).

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the current. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the vectorgrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the current. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the vectorgrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

  • jdf_eastward (numpy.ndarray, float) – eastward component of the divergence-free current evalulated at the coordinates given by the vectorgrid attribute
  • jdf_northward (numpy.ndarray, float) – northward component of the divergence-free current evalulated at the coordinates given by the vectorgrid attribute

See also

get_curl_free_current()
Calculate curl-free part of the current
get_total_current()
Calculate total horizontal current
get_divergence_free_current_function(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the divergence-free current function

Isocontours of the divergence-free current function indicates the alignment of the divergence-free part of the horizontal current. Its direction is given by the cross product between a vertical vector and the gradient of the divergence-free current function. A fixed amount of current flows between isocontours. The calculations refer to the height chosen upon initialization of the AMPS object (default 110 km). Divergence-free current function unit is kA.

Note

The divergence-free current is similar to what is often termed the equivalent current, that is derived from ground based magnetic field measurements. However, the present divergence-free current is derived from measurements above the ionosphere, and thus it contains signal both from ionospheric currents below low Earth orbit, and from subsurface induced currents. See Laundal et al. (2016) [1] where this current is called Psi for more detail.

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the current. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the current. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

Psi – Divergence-free current function evaluated at self.scalargrid, or, if specified, mlat/mlt
Return type:

numpy.ndarray

References

[1]K. M. Laundal, C. C. Finlay, and N. Olsen, “Sunlight effects on the 3D polar current system determined from low Earth orbit measurements” Earth, Planets and Space, 2016, https://doi.org/10.1186/s40623-016-0518-x
get_ground_Beqd(height=0)[source]

Calculate ground magnetic field perturbations in the QD east direction, in units of nT.

Note

These calculations are made by assuming that the divergende-free current function calculated with the AMPS model correspond to the equivalent current function of an external magnetic potential, as described by Chapman & Bartels 1940 [2]_. Induced components are thus ignored. The height of the current function also becomes important when propagating the model values to the ground.

Also note that the output parameters will be QD components, and that they can be converted to geographic by use of QD base vectors [3]_

This function is not optimized for calculating long time series of model ground magnetic field perturbations, although it is possible to use for that.

Parameters:height (float, optional) – height, in km, where the magnetic field is evaluated. Must be less than self.height, which is the height of the current. Default is 0 (ground).
Returns:dB_east – Eastward component of the magnetic field disturbance on ground
Return type:numpy.ndarray

References

[2]
  1. Chapman & J. Bartels “Geomagnetism Vol 2” Oxford University Press 1940
[3]A. D. Richmond, “Ionospheric Electrodynamics Using Magnetic Apex Coordinates”, Journal of geomagnetism and geoelectricity Vol. 47, 1995, http://doi.org/10.5636/jgg.47.191

See also

get_ground_Bnqd()
Calculate ground perturbation in northward qd direction on scalargrid
get_ground_Buqd()
Calculate ground perturbation in upward qd direction on scalargrid
get_ground_perturbation()
Calculate ground perturbation in east/north qd direction
get_ground_Bnqd(height=0)[source]

Calculate ground magnetic field perturbations in the QD north direction, in units of nT.

Note

These calculations are made by assuming that the divergende-free current function calculated with the AMPS model correspond to the equivalent current function of an external magnetic potential, as described by Chapman & Bartels 1940 [2]_. Induced components are thus ignored. The height of the current function also becomes important when propagating the model values to the ground.

Also note that the output parameters will be QD components, and that they can be converted to geographic by use of QD base vectors [3]_

This function is not optimized for calculating long time series of model ground magnetic field perturbations, although it is possible to use for that.

Parameters:height (float, optional) – height, in km, where the magnetic field is evaluated. Must be less than self.height, which is the height of the current. Default is 0 (ground).
Returns:dB_north – Northward component of the magnetic field disturbance on ground
Return type:numpy.ndarray

References

[2]
  1. Chapman & J. Bartels “Geomagnetism Vol 2” Oxford University Press 1940
[3]A. D. Richmond, “Ionospheric Electrodynamics Using Magnetic Apex Coordinates”, Journal of geomagnetism and geoelectricity Vol. 47, 1995, http://doi.org/10.5636/jgg.47.191

See also

get_ground_Beqd()
Calculate ground perturbation in easthward qd direction on scalargrid
get_ground_Buqd()
Calculate ground perturbation in upward qd direction on scalargrid
get_ground_perturbation()
Calculate ground perturbation in east/north qd direction
get_ground_Buqd(height=0.0)[source]

Calculate ground magnetic field perturbations in the QD up direction, in units of nT.

Note

These calculations are made by assuming that the divergende-free current function calculated with the AMPS model correspond to the equivalent current function of an external magnetic potential, as described by Chapman & Bartels 1940 [2]_. Induced components are thus ignored. The height of the current function also becomes important when propagating the model values to the ground.

Also note that the output parameters will be QD components, and that they can be converted to geographic by use of QD base vectors [3]_

This function is not optimized for calculating long time series of model ground magnetic field perturbations, although it is possible to use for that.

Parameters:height (float, optional) – height, in km, where the magnetic field is evaluated. Must be less than self.height, which is the height of the current. Default is 0 (ground).
Returns:dB_up – Upward component of the magnetic field disturbance on ground
Return type:numpy.ndarray

References

[2]
  1. Chapman & J. Bartels “Geomagnetism Vol 2” Oxford University Press 1940
[3]A. D. Richmond, “Ionospheric Electrodynamics Using Magnetic Apex Coordinates”, Journal of geomagnetism and geoelectricity Vol. 47, 1995, http://doi.org/10.5636/jgg.47.191

See also

get_ground_Beqd()
Calculate ground perturbation in easthward qd direction on scalargrid
get_ground_Bnqd()
Calculate ground perturbation in northward qd direction on scalargrid
get_ground_perturbation()
Calculate ground perturbation in east/north qd direction
get_ground_perturbation(mlat=<object object>, mlt=<object object>, height=0)[source]

Calculate magnetic field perturbations on ground, in units of nT, that corresponds to the divergence-free current function.

Parameters:
  • mlat (numpy.ndarray, float, optional) – magnetic latitude of the output. The array shape will not be preserved, and the results will be returned as a 1-dimensional array. Default value is from self.vectorgrid
  • mlt (numpy.ndarray, float, optional) – magnetic local time of the output. The array shape will not be preserved, and the results will be returned as a 1-dimensional array. Default value is from self.vectorgrid
  • height (numpy.ndarray, optional) – geodetic height at which the field should be evalulated. Should be < current height set at initialization. Default 0 (ground)

Note

These calculations are made by assuming that the divergende-free current function calculated with the AMPS model correspond to the equivalent current function of an external magnetic potential, as described by Chapman & Bartels 1940 [2]_. Induced components are thus ignored. The height of the current function also becomes important when propagating the model values to the ground.

Also note that the output parameters will be QD components, and that they can be converted to geographic by use of QD base vectors [3]_

This function is not optimized for calculating long time series of model ground magnetic field perturbations, although it is possible to use for that.

Returns:
  • dB_east (numpy.ndarray) – Eastward component of the magnetic field disturbance on ground
  • dB_north (numpy.ndarray) – Northward component of the magnetic field disurubance on ground

References

[2]
  1. Chapman & J. Bartels “Geomagnetism Vol 2” Oxford University Press 1940
[3]A. D. Richmond, “Ionospheric Electrodynamics Using Magnetic Apex Coordinates”, Journal of geomagnetism and geoelectricity Vol. 47, 1995, http://doi.org/10.5636/jgg.47.191
get_integrated_upward_current()[source]

Calculate the integrated upward and downward current, poleward of minlat, in units of MA.

Note

This calculation uses the scalargrid attribute. By default this is a regular grid, with coordinates from north and south tiled side-by-side (equal number of coords, with north first). If scalargrid has been changed, and has a different structure, this function will return wrong values.

Returns:
  • J_up_n (float) – Total upward current in the northern hemisphere
  • J_down_n (float) – Total downward current in the northern hemisphere
  • J_up_s (float) – Total upward current in the southern hemisphere
  • J_down_s (float) – Total downward current in the southern hemisphere
get_poloidal_scalar(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the poloidal scalar potential values (unit is microTm).

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the poloidal scalar. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the poloidal scalar. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

V – Poloidal scalar evalulated at self.scalargrid, or, if specified, mlat/mlt
Return type:

numpy.ndarray
get_toroidal_scalar(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the toroidal scalar values (unit is nT).

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the toroidal scalar. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the toroidal scalar. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

T – Toroidal scalar evaluated at self.scalargrid, or, if specified, mlat/mlt
Return type:

numpy.ndarray
get_total_current(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the total horizontal current, in units of mA/m. This is calculated as the sum of the curl-free and divergence-free parts. The calculations refer to the height chosen upon initialization of the AMPS object (default 110 km).

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the current. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the vectorgrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the current. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the vectorgrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

  • j_eastward (numpy.ndarray, float) – eastward component of the horizontal current evalulated at the coordinates given by the vectorgrid attribute
  • j_northward (numpy.ndarray, float) – northward component of the horizontal current evalulated at the coordinates given by the vectorgrid attribute

See also

get_divergence_free_current()
Calculate divergence-free part of the horizontal current
get_curl_free_current()
Calculate curl-free part of the horizontal current
get_total_current_magnitude()[source]

Calculate the total horizontal current density magnitude, in units of mA/m. This is calculated as the sum of the curl-free and divergence-free parts. The calculations refer to the height chosen upon initialization of the AMPS object (default 110 km). The calculations are performed on the coordinates of self.scalargrid. This is useful for making contour plots of the horizontal current density magnitude, and faster than calculating the magnitude from the output of get_total_current

Returns:j – horizontal current density magnitude, evalulated at the coordinates given by the scalargrid attribute
Return type:numpy.ndarray, float

See also

get_total_current()
Calculate total current density vector components
get_upward_current(mlat=<object object>, mlt=<object object>, grid=False)[source]

Calculate the upward current (unit is microAmps per square meter). The calculations refer to the height chosen upon initialization of the AMPS object (default 110 km).

Parameters:
  • mlat (numpy.ndarray, optional) – array of mlats at which to calculate the current. Will be ignored if mlt is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • mlt (numpy.ndarray, optional) – array of mlts at which to calculate the current. Will be ignored if mlat is not also specified. If not specified, the calculations will be done using the coords of the scalargrid attribute.
  • grid (bool, optional, default False) – if True, mlat and mlt are interpreted as coordinates in a regular grid. They must be 1-dimensional, and the output will have dimensions len(mlat) x len(mlt). If mlat and mlt are not set, this keyword is ignored.
Returns:

Ju – Upward current evaulated at self.scalargrid, or, if specified, mlat/mlt
Return type:

numpy.ndarray
plot_currents(vector_scale=200)[source]

Create a summary plot of the current fields

Parameters:vector_scale (optional) – Current vector lengths will be shown relative to a template. This parameter determines the magnitude of that template, in mA/m. Default is 200 mA/m

Examples

>>> # initialize by supplying a set of external conditions:
>>> m = AMPS(300, # solar wind velocity in km/s
             -4, # IMF By in nT
             -3, # IMF Bz in nT
             20, # dipole tilt angle in degrees
             150) # F10.7 index in s.f.u.
>>> # make summary plot:
>>> m.plot_currents()
update_model(v, By, Bz, tilt, f107, coeff_fn=<object object>)[source]

Update the model vectors without updating all the other matrices. This leads to better performance than just making a new AMPS object.

Parameters:
  • v (float) – solar wind velocity in km/s
  • By (float) – IMF GSM y component in nT
  • Bz (float) – IMF GSM z component in nT
  • tilt (float) – dipole tilt angle in degrees
  • f107 (float) – F10.7 index in s.f.u.

Examples

If model currents shall be calculated on the same grid for a range of external conditions, it is faster to do this:

>>> m1 = AMPS(solar_wind_velocity_in_km_per_s, IMF_By_in_nT, IMF_Bz_in_nT, dipole_tilt_in_deg, F107_index)
>>> # ... current calculations ...
>>> m1.update_model(new_v, new_By, new_Bz, new_tilt, new_f107)
>>> # ... new current calcuations ...

than to make a new object:

>>> m2 = AMPS(new_v, new_By, new_Bz, new_tilt, new_f107)
>>> # ... new current calculations ...

Also note that the inputs are scalars in both cases. It is possible to optimize the calculations significantly by allowing the inputs to be arrays. That is not yet implemented.

pyamps.get_B_ground(qdlat, mlt, height, v, By, Bz, tilt, f107, current_height=110, epsilon_multiplier=1.0, chunksize=25000, coeff_fn='/home/docs/checkouts/readthedocs.org/user_builds/pyamps/checkouts/latest/pyamps/coefficients/SW_OPER_MIO_SHA_2E_00000000T000000_99999999T999999_0103.txt')[source]

Calculate model magnetic field on ground

This function uses dask to parallelize computations. That means that it is quite fast and that the memory consumption will not explode unless chunksize is too large.

Parameters:
  • qdlat (array_like or float) – quasi-dipole latitude, in degrees. Can be either a scalar (float), or an array with an equal number of elements as mlt
  • mlt (array_like) – array of magnetic local times (hours)
  • height (float) – geodetic height, in km (0 <= height <= current_height)
  • v (array_like) – array of solar wind velocities in GSM/GSE x direction (km/s)
  • By (array_like) – array of solar wind By values (nT)
  • Bz (array_like) – array of solar wind Bz values (nT)
  • tilt (array_like) – array of dipole tilt angles (degrees)
  • f107 (array_like) – array of F10.7 index values (SFU)
  • current_height (float, optional) – height (km) of the current sheet. Default is 110.
  • epsilon_multiplier (float, optional) – multiplier for the epsilon parameter. Default is 1.
  • chunksize (int) – the input arrays will be split in chunks in order to parallelize computations. Larger chunks consumes more memory, but might be faster. Default is 25000.
  • coeff_fn (str, optional) – file name of model coefficients - must be in format produced by model_vector_to_txt.py (default is latest version)
Returns:

  • Bqphi (array_like) – magnetic field in quasi-dipole eastward direction
  • Bqlambda (array_like) – magnetic field in quasi-dipole northward direction
  • Bqr (array_like) – magnetic field in upward direction. See notes

Note

We assume that there are no induced currents. The error in this assumption will be larger for the radial component than for the horizontal components

Array inputs should have the same dimensions.

pyamps.get_B_space(glat, glon, height, time, v, By, Bz, tilt, f107, epoch=2015.0, h_R=110.0, chunksize=15000, coeff_fn='/home/docs/checkouts/readthedocs.org/user_builds/pyamps/checkouts/latest/pyamps/coefficients/SW_OPER_MIO_SHA_2E_00000000T000000_99999999T999999_0103.txt')[source]

Calculate model magnetic field in space

This function uses dask to parallelize computations. That means that it is quite fast and that the memory consumption will not explode unless chunksize is too large.

Parameters:
  • glat (array_like) – array of geographic latitudes (degrees)
  • glon (array_like) – array of geographic longitudes (degrees)
  • height (array_like) – array of geodetic heights (km)
  • time (array_like) – list/array of datetimes, needed to calculate magnetic local time
  • v (array_like) – array of solar wind velocities in GSM/GSE x direction (km/s)
  • By (array_like) – array of solar wind By values (nT)
  • Bz (array_like) – array of solar wind Bz values (nT)
  • tilt (array_like) – array of dipole tilt angles (degrees)
  • f107 (array_like) – array of F10.7 index values (SFU)
  • epoch (float, optional) – epoch (year) used in conversion to magnetic coordinates with the IGRF. Default = 2015.
  • h_R (float, optional) – reference height (km) used when calculating modified apex coordinates. Default = 110.
  • chunksize (int, optional) – the input arrays will be split in chunks in order to parallelize computations. Larger chunks consumes more memory, but might be faster. Default is 15000.
  • coeff_fn (str, optional) – file name of model coefficients - must be in format produced by model_vector_to_txt.py (default is latest version)
Returns:

  • Be (array_like) – array of model magnetic field (nT) in geodetic eastward direction (same dimension as input)
  • Bn (array_like) – array of model magnetic field (nT) in geodetic northward direction (same dimension as input)
  • Bu (array_like) – array of model magnetic field (nT) in geodetic upward direction (same dimension as input)

Note

Array inputs should have the same dimensions.

pyamps.mlon_to_mlt(mlon, times, epoch)[source]

Calculate magnetic local time from magnetic longitude and time(s).

This is an implementation of the formula recommended in Laundal & Richmond, 2017 [4]_. It uses the subsolar point geomagnetic (centered dipole) longitude to define the noon meridian.

Parameters:
  • mlon (array_like) – array of magnetic longitudes
  • times (datetime or list of datetimes) – datetime object, or list of datetimes with equal number of elements as mlon
  • epoch (float) – the epoch (year, ) used for geo->mag conversion

References

[4]Laundal, K.M. & Richmond, A.D. Space Sci Rev (2017) 206: 27. https://doi.org/10.1007/s11214-016-0275-y

See also

mlt_to_mlon()
Magnetic longitude from magnetic local time and universal time

Magnetic local time (pyamps.mlt_utils)

pyamps.mlt_utils.car_to_sph(car, deg=True)[source]

Convert from spherical to cartesian coordinates

Parameters:
  • car (3 x N array) – 3 x N array, where the rows are, from top to bottom: x, y, z, in ECEF coordinates
  • deg (bool, optional) – set to True if output is wanted in degrees. False if radians
Returns:

sph – 3 x N array, where the rows are, from top to bottom: radius, colatitude, and longitude
Return type:

3 x N array
pyamps.mlt_utils.geo2mag(glat, glon, epoch, deg=True, inverse=False)[source]

Convert geographic (geocentric) to centered dipole coordinates

The conversion uses IGRF coefficients directly, interpolated to the provided epoch. The construction of the rotation matrix follows Laundal & Richmond (2017) [4]_ quite directly.

Parameters:
  • glat (array_like) – array of geographic latitudes
  • glon (array_like) – array of geographic longitudes
  • epoch (float) – epoch (year) for the dipole used in the conversion
  • deg (bool, optional) – True if input is in degrees, False otherwise
  • inverse (bool, optional) – set to True to convert from magnetic to geographic. Default is False
Returns:

  • cdlat (ndarray) – array of centered dipole latitudes [degrees]
  • cdlon (ndarray) – array of centered dipole longitudes [degrees]
pyamps.mlt_utils.is_leapyear(year)[source]

Return True if leapyear else False

Handles arrays and preserves shape

Parameters:year (array_like) – array of years
Returns:is_leapyear – True where input is leapyear, False elsewhere
Return type:ndarray of bools
pyamps.mlt_utils.mlon_to_mlt(mlon, times, epoch)[source]

Calculate magnetic local time from magnetic longitude and time(s).

This is an implementation of the formula recommended in Laundal & Richmond, 2017 [4]_. It uses the subsolar point geomagnetic (centered dipole) longitude to define the noon meridian.

Parameters:
  • mlon (array_like) – array of magnetic longitudes
  • times (datetime or list of datetimes) – datetime object, or list of datetimes with equal number of elements as mlon
  • epoch (float) – the epoch (year, ) used for geo->mag conversion

References

[4]Laundal, K.M. & Richmond, A.D. Space Sci Rev (2017) 206: 27. https://doi.org/10.1007/s11214-016-0275-y

See also

mlt_to_mlon()
Magnetic longitude from magnetic local time and universal time
pyamps.mlt_utils.mlt_to_mlon(mlt, times, epoch)[source]

Calculate quasi-dipole magnetic longitude from magnetic local time and universal time(s).

This is an implementation of the formula recommended in Laundal & Richmond, 2017 [4]_. It uses the subsolar point geomagnetic (centered dipole) longitude to define the noon meridian.

Parameters:
  • mlt (array_like) – array of magnetic local times
  • times (datetime or list of datetimes) – datetime object, or list of datetimes with equal number of elements as mlon
  • epoch (float) – the epoch (year, ) used for geo->mag conversion

References

[4]Laundal, K.M. & Richmond, A.D. Space Sci Rev (2017) 206: 27. https://doi.org/10.1007/s11214-016-0275-y

See also

mlon_to_mlt()
Magnetic local time from magnetic longitude and universal time
pyamps.mlt_utils.sph_to_car(sph, deg=True)[source]

Convert from spherical to cartesian coordinates

Parameters:
  • sph (3 x N array) – 3 x N array, where the rows are, from top to bottom: radius, colatitude, and longitude
  • deg (bool, optional) – set to True if input is given in degrees. False if radians
Returns:

car – 3 x N array, where the rows are, from top to bottom: x, y, z, in ECEF coordinates
Return type:

3 x N array
pyamps.mlt_utils.subsol(datetimes)[source]

calculate subsolar point at given datetime(s)

Parameters:datetimes (datetime or list of datetimes) – datetime or list (or other iterable) of datetimes
Returns:
  • subsol_lat (ndarray) – latitude(s) of the subsolar point
  • subsol_lon (ndarray) – longiutde(s) of the subsolar point
Raises:ValueError – if any datetime.year value provided is not within (1600,2100)

Note

The code is vectorized, so it should be fast.

After Fortran code by: 961026 A. D. Richmond, NCAR

Documentation from original code: Find subsolar geographic latitude and longitude from date and time. Based on formulas in Astronomical Almanac for the year 1996, p. C24. (U.S. Government Printing Office, 1994). Usable for years 1601-2100, inclusive. According to the Almanac, results are good to at least 0.01 degree latitude and 0.025 degree longitude between years 1950 and 2050. Accuracy for other years has not been tested. Every day is assumed to have exactly 86400 seconds; thus leap seconds that sometimes occur on December 31 are ignored: their effect is below the accuracy threshold of the algorithm.

Spherical harmonics (pyamps.sh_utils)

class pyamps.sh_utils.SHkeys(Nmax, Mmax)[source]

container for n and m in spherical harmonics

keys = SHkeys(Nmax, Mmax)

keys will behave as a tuple of tuples, more or less keys[‘n’] will return a list of the n’s keys[‘m’] will return a list of the m’s keys[3] will return the fourth n,m tuple

keys is also iterable

Mge(limit)[source]

set m >= limit

MleN()[source]

set m <= n

NminusMeven()[source]

remove keys if n - m is odd

NminusModd()[source]

remove keys if n - m is even

make_arrays()[source]

prepare arrays with shape ( 1, len(keys) ) these are used when making G matrices

negative_m()[source]

add negative m to the keys

setNmin(nmin)[source]

set minimum n

pyamps.sh_utils.getG0(glat, glon, height, time, epoch=2015.0, h_R=110.0, NT=65, MT=3, NV=45, MV=3)[source]

calculate the G matrix for the constant term in the AMPS model. The constant term is the term that depends only on the spherical harmonic coefficients that are not scaled by external parameters. This G matrix can be used to produce the full matrix.

The structure of the matrix is such that G0.dot(m), the product of the matrix with the first 1/19 of the model vector, will be model values of the eastward, northward, and upward components of the magnetic field, stacked in a column vector.

glat, glon, time, and height must all have the same number of elements (let’s call this N)

Parameters:
  • glat (array) – Geodetic latitude (degrees)
  • glon (array) – Geographic/geodetic longitude (degrees)
  • height (array) – Geodetic heights, in km
  • time (array) – Array of datetimes corresponding to each point. This is needed to calculate magnetic local time.
  • epoch (float, optional) – The epoch used for conversion to apex coordinates. Default 2015.
  • h_R (float, optional) – Reference height used in conversion to modified apex coordinates. Default 110 km.
  • MT, NV, MV (NT,) – Truncation level. Must match coefficient file
Returns:

G0 – an 3N by M matrix, where N is the number of elements in the input coordinates. There will be 3 times as many rows G0, since there are 3 components. Partionining G0 in thirds, from top to bottom, gives the parts that correspond to east, north, and up, respectively. M is the number of terms in the spherical harmonic expansion of B
Return type:

array
pyamps.sh_utils.get_ground_field_G0(qdlat, mlt, height, current_height, N=45, M=3)[source]

calculate the G matrix for the constant term in the AMPS model needed to calculate corresponding ground magnetic field perturbations. The constant term is the term that depends only on the spherical harmonic coefficients that are not scaled by external parameters. This G matrix can be used to produce the full matrix.

Parameters:
  • qdlat (array) – Quasi-dipole latitude (degrees)
  • MLT (array) – Magnetic local time (MLT)
  • height (float) – Geodetic height, in km (0 <= height <= current_height)
  • current_height (float) – height of the current
  • M (N,) –
Returns:

G0 – an 3N by M matrix, where N is the number of elements in the input coordinates. There will be 3 times as many rows G0, since there are 3 components. Partionining G0 in thirds, from top to bottom, gives the parts that correspond to east, north, and up, respectively. M is the number of terms in the spherical harmonic expansion of B
Return type:

array
pyamps.sh_utils.legendre(nmax, mmax, theta, schmidtnormalize=True, keys=None)[source]

Calculate associated Legendre function P and its derivative

Algorithm from “Spacecraft Attitude Determination and Control” by James Richard Wertz

Parameters:
  • nmax (int) – highest spherical harmonic degree
  • mmax (int) – hightest spherical harmonic order
  • theta (array, float) – colatitude in degrees (shape is not preserved)
  • schmidtnormalize (bool, optional) – True if Schmidth seminormalization is wanted, False otherwise. Default True
  • keys (SHkeys, optional) – If this parameter is set, an array will be returned instead of a dict. The array will be (N, 2M), where N is the number of elements in theta, and M is the number of keys. The first M columns represents a matrix of P values, and the last M columns represent values of dP/dtheta
Returns:

  • P (dict) – dictionary of Legendre function evalulated at theta. Dictionary keys are spherical harmonic wave number tuples (n, m), and values will have shape (N, 1), where N is number of elements in theta.
  • dP (dict) – dictionary of Legendre function derivatives evaluated at theta. Dictionary keys are spherical harmonic wave number tuples (n, m), and values will have shape (N, 1), where N is number of elements in theta.
  • PdP (array (only if keys != None)) – if keys != None, PdP is returned instaed of P and dP. PdP is an (N, 2M) array, where N is the number of elements in theta, and M is the number of keys. The first M columns represents a matrix of P values, and the last M columns represent values of dP/dtheta
pyamps.sh_utils.nterms(NT=0, MT=0, NVi=0, MVi=0, NVe=0, MVe=0)[source]

return number of coefficients in an expansion in real spherical harmonics of toroidal magnetic potential truncated at NT, MT poloidal magnetic potential truncated at NVi, MVi for internal sources poloidal magnetic potential truncated at NVe, MVe for external sources

AMPS model parameters (pyamps.model_utils)

pyamps.model_utils.get_m_matrix(coeff_fn='/home/docs/checkouts/readthedocs.org/user_builds/pyamps/checkouts/latest/pyamps/coefficients/SW_OPER_MIO_SHA_2E_00000000T000000_99999999T999999_0103.txt')[source]

make matrix of model coefficients - used in get_B_space for fast calculations of model field time series along trajectory with changing input

pyamps.model_utils.get_m_matrix_pol(coeff_fn='/home/docs/checkouts/readthedocs.org/user_builds/pyamps/checkouts/latest/pyamps/coefficients/SW_OPER_MIO_SHA_2E_00000000T000000_99999999T999999_0103.txt')[source]

make matrix of model coefficients - only poloidal part used in get_B_ground for fast calculations of model field time series on ground

pyamps.model_utils.get_model_vectors(v, By, Bz, tilt, f107, epsilon_multiplier=1.0, coeff_fn='/home/docs/checkouts/readthedocs.org/user_builds/pyamps/checkouts/latest/pyamps/coefficients/SW_OPER_MIO_SHA_2E_00000000T000000_99999999T999999_0103.txt')[source]

tor_c, tor_s, pol_c, pol_s = get_model_vectors(v, By, Bz, tilt, F107, epsilon_multiplier = 1., coeffs = coeffs)

returns column vectors ((K,1)-shaped) corresponding to the spherical harmonic coefficients of the toroidal and poloidal parts, with _c and _s denoting cos and sin terms, respectively.

This function is used by amps.AMPS class

pyamps.model_utils.get_truncation_levels(coeff_fn='/home/docs/checkouts/readthedocs.org/user_builds/pyamps/checkouts/latest/pyamps/coefficients/SW_OPER_MIO_SHA_2E_00000000T000000_99999999T999999_0103.txt')[source]

read model truncation levels from coefficient file returns NT, MT, NV, MV (spherical harmonic degree (N) and order (M) for toroidal (T) and poloidal (V) fields, respectively)

Polar coordinate plots (pyamps.plot_utils)

pyamps.plot_utils.equal_area_grid(dr=2, K=0, M0=8, N=20)[source]

function for calculating an equal area grid in polar coordinates

Parameters:
  • dr (float (optional)) – The latitudinal resolution of the grid. Default 2 degrees
  • K (int (optional)) – This number determines the colatitude of the inner radius of the post poleward ring of grid cells (the pole is not inlcluded!). It relates to this colatitude (r0) as r0/dr = (2K + 1)/2 => K = (2r0/dr - 1)/2. Default value is 0
  • M0 (int (optional)) – The number of sectors in the most poleward ring. Default is 8
  • N (int (optional)) – The number of rings to be included. This determiend how far equatorward the grid extends. Typically if dr is changed from 2 to 1, N should be doubled to reach the same latitude. Default is 20, which means that the equatorward edge of the grid is 89 - 20*2 = 49 degrees (the most poleward latitude is 89 with default values)
Returns:

  • mlat (array) – Array of latitudes for the equatorward west (“lower left”) corners of the grid cells.
  • mlt (array) – Array of magnetic local times for the equatorward west (“lower left”) corner of the grid cells.
  • mltres (array) – width, in magnetic local time, of the cells with lower left corners described by mlat and mlt. Notice that this width changes with latitude, while the latitudinal width is fixed, determined by the dr parameter

Indices and tables