viscid.cotr module

Numpy implementation of space physics coordinate transformations

The abbreviations and equations in this module follow [Hapgood1992],

• gei: Geocentric equatorial inertial

  • X = First Point of Aries
  • Z = Geographic North Pole

• geo: Geographic

  • X = Intersection of Greenwich meridian and geographic equator
  • Z = Geographic North Pole

• gse: Geocentric solar ecliptic

  • X = Earth-Sun line (points toward Sun)
  • Z = Ecliptic North Pole

• gsm: Geocentric solar magnetospheric

  • X = Earth-Sun line (points toward Sun)
  • Z = Projection of dipole axis on GSE Y-Z plane

• sm: Solar magnetic

  • Y = Perpendicular to plane containing Earth-Sun line and dipole axis. Positive is opposite to Earth’s orbital motion
  • Z = Earth’s dipole axis

• mag: Geomagnetic

  • Y = Intersection between geographic equator and the geographic meridian 90deg East of the meditian containing the dipole axis
  • Z = Earth’s dipole axis

• mhd: OpenGGCM’s native coordinate system (based on GSE)

  • X = Sun-Earth line (points AWAY from Sun)
  • Z = Ecliptic North Pole

• hae: Heliocentric Aries ecliptic

  • X = First Point of Aries
  • Z = Ecliptic North Pole

• hee: Heliocentric Earth ecliptic

  • X = Sun-Earth line
  • Z = Ecliptic North Pole

• heeq: Heliocentric Earth equatorial

  • X = Intersection between solar equator and solar central meridian as seen from Earth
  • Z = North Pole of solar rotation

Dipole info between 2010-01-01 00:00 and 2010-06-20 23:59,

• Smallest dipole tilt angle: 0.00639 deg @ 2010-04-16 05:10
• Largest dipole tilt angle: 33.5 deg @ 2010-06-19 16:59
• Smallest angle between dip and GSE-Z: 13.3 deg @ 2010-01-03 15:54
• Largest angle between dip and GSE-Z: 33.6 deg @ 2010-01-13 03:19

Dipole info between 1952-01-01 00:00 and 1955-12-31 23:59,

• Smallest dipole tilt angle: -0.00145 deg @ 1952-09-01 01:35
• Largest dipole tilt angle: 35.2 deg @ 1954-06-21 16:38
• Smallest angle between dip and GSE-Z: 11.7 deg @ 1954-12-14 17:05
• Largest angle between dip and GSE-Z: 35.2 deg @ 1955-01-02 03:53

Note

IGRF coefficients are taken from [IGRF], and coefficients are assumed constant outside of the given range. This means every 5 years, the new coefficients should be added manually.

References

[IGRF]

<http://wdc.kugi.kyoto-u.ac.jp/igrf/coef/igrf12coeffs.txt>

[Hapgood1992]

Hapgood, M. A., 2011, “Space Physics Coordinate Transformations: A User Guide”, Planet. Space Sci. Vol 40, No. 5. pp. 711-717, 1992. <http://www.igpp.ucla.edu/public/vassilis/ESS261/Lecture03/Hapgood_sdarticle.pdf>

This module only depends on npdatetime, which is itself orthogonal to Viscid, so these two modules can be ripped out and used more generally. Please note that Viscid is MIT licensed, which requires attribution.

The MIT License (MIT) Copyright (c) 2017 Kristofor Maynard

viscid.cotr.as_nvec(arr, ndim=4, init_vals=(0, 0, 0, 1))

Extend vectors to ndim dimensions

The first dimension or arr that is the one that is extended. If the last dimension has shape >= 4, leave arr unchanged

Parameters:

• arr (sequence) – array with shape (3,) or (3, N) for N vectors
• ndim (int) – extend arr to this many dimensions
• init_vals (sequence) – initialized values for new dimensions, note that init_vals[i] is the initialization for the ith dimension of the result. (0, 0, 0, 1) is useful for 4x4 rotation+translation matrices

Returns:

(ndim,) or (ndim, N) depending on input

Return type:

ndarray

viscid.cotr.make_rotation(theta, axis='z', rdim=3)

Make rotation matrix of theta degrees around axis

Warning

The Hapgood paper uses the convention that rotations around y are right-handed, but rotations around x and z are left handed!

Parameters:

• theta (float) – angle (degrees)
• axis (int, str) – one of (0, ‘x’, 1, ‘y’, 2, ‘z’)
• rdim (int) – 3 to make 3x3 matrices, or 4 to make 4x4 where mat[:, 3] == mat[3, :] == [0, 0, 0, 1]. 4x4 matrices are useful for performing translations.

Raises:

ValueError – On invalid axis

viscid.cotr.make_translation(v)

Make a 4x4 translation matrix to act on [x, y, z, 1] vectors

Parameters:v (sequence) – Translation vector, must have at least 3 entries Returns:

mat will translate an [x, y, z, 1] vector, and

imat will do the inverse translation

Return type:(mat, imat)

viscid.cotr.get_rotation_wxyz(mat, check_inverse=True, quick=True)

Find rotation axis and angle of a rotation matrix

Parameters:

• mat (ndarray) – 3x3 rotation matrix with determinant +1
• check_inverse (bool) – check if mat.T == inverse(mat)
• quick (bool) – Try a quicker, less well tested method first

Returns:

[w (rotation angle), ux, uy, uz] where u is the

normalized axis vector

Return type:

ndarray

viscid.cotr.as_crd_system(thing, default=<class 'viscid.cotr._NOT_GIVEN'>)

Try to figure out a crd_system for thing

Parameters:

• thing – Either a string with a crd_system abbreviation, or something that that can do thing.find_info(‘crd_system’)
• default (str) – fallback value

Returns:

crd system abbreviation, stripped and lower case

Return type:

str

Raises:

TypeError – if no default supplied and thing can’t be turned into a crd_system

class viscid.cotr.Cotr(time='1967-01-01', dip_tilt=None, dip_gsm=None, rdim=3, notilt1967=True)

Bases: object

Transform geocentric vectors between crd systems @ a given UTC

GEOCENTRIC = ['gei', 'geo', 'gse', 'gsm', 'sm', 'mag', 'mhd']

HELIOCENTRIC = ['hea', 'hee', 'heeq']

Q_e

Q_g

XFORM_LOOKUP = {'gei->gei': (0,), 'geo->gei': (-1,), 'geo->geo': (0,), 'gse->gei': (-2,), 'gse->geo': (1, -2), 'gse->gse': (0,), 'gsm->gei': (-2, -3), 'gsm->geo': (1, -2, -3), 'gsm->gse': (-3,), 'gsm->gsm': (0,), 'hae->hae': 0, 'hee->hae': -11, 'hee->hee': 0, 'heeq->hae': -12, 'heeq->heeq': 0, 'mag->gei': (-1, -5), 'mag->geo': (-5,), 'mag->gse': (2, -1, -5), 'mag->gsm': (3, 2, -1, -5), 'mag->mag': (0,), 'mag->sm': (4, 3, 2, -1, -5), 'mhd->gei': (-2, 6), 'mhd->geo': (1, -2, 6), 'mhd->gse': (6,), 'mhd->gsm': (3, 6), 'mhd->mag': (5, 1, -2, 6), 'mhd->mhd': (0,), 'mhd->sm': (4, 3, 6), 'sm->gei': (-2, -3, -4), 'sm->geo': (1, -2, -3, -4), 'sm->gse': (-3, -4), 'sm->gsm': (-4,), 'sm->sm': (0,)}

dip_geolat

dip_geolat_rad

dip_geolon

dip_geolon_rad

dip_gsm

Angle psi between GSE-Z and dipole projected into GSE-YZ

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

dip_mu

Dipole tilt angle mu (angle between GSM-Z and dipole)

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

dip_psi

Angle psi between GSE-Z and dipole projected into GSE-YZ

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

dip_tilt

Dipole tilt angle mu (angle between GSM-Z and dipole)

Note

SM -> GSE, rotation is <-mu, Y> * <-psi, X> which is equivalent to <psi, X> * <mu, Y>

get_dipole_angles()

Get rotation angles between dipole and GSE

Note

psi refers to the angle between GSE and GSM, while mu is the dipole tilt angle.

Returns:

The SM -> GSE, rotation is <-mu, Y> * <-psi, X>

which is equivalent to <psi, X> * <mu, Y>

Return type:(tilt, gsm)

get_dipole_moment(crd_system='gse', strength=32707.5292732387)

Get Earth’s dipole moment

Parameters:

• crd_system (str) – crd_system of resulting moment
• strength (float) – magnitude of dipole moment

Returns:

3-vector of dipole tilt in crd_system

Return type:

ndarray

get_rotation_wxyz(from_system, to_system)

Get rotation axis and angle for a transformation

Parameters:

• from_system (str) – abbreviation of crd_system
• to_system (str) – abbreviation of crd_system

Returns:

[w, x, y, z] where w is the angle around the axis

vector [x, y, z]

Return type:

ndarray

get_transform_matrix(from_system, to_system)

Get a transformation matrix from one system to another

Parameters:

• from_system (str) – abbreviation of crd_system
• to_system (str) – abbreviation of crd_system

Returns:

self.rdim x self.rdim transformation matrix

Return type:

ndarray

sun_ecliptic_longitude

transform(from_system, to_system, vec)

Transform a vector from one system to another

Parameters:

• from_system (str) – abbreviation of crd_system
• to_system (str) – abbreviation of crd_system
• vec (ndarray) – should have shape (3, ) or (3, N) to transform N vectors at once

Returns:

vec in new crd system shaped either (3, ) or

(3, N) mirroring the shape of vec

Return type:

ndarray

viscid.cotr.as_cotr(thing=None, default=<class 'viscid.cotr._NOT_GIVEN'>)

Try to make cotr into a Cotr instance

Inputs understood are:

• None for North-South dipole
• anything with a __cotr__ property
• datetimes or similar, specifies the time for finding dip angles
• mapping (dict or similar), passed to Cotr constructor as kwargs

Parameters:

• thing (obj) – something to turn into a Cotr
• default (None, obj) – fallback value

Returns:

Cotr instance

Raises:

TypeError – if no default supplied and thing can’t be turned into a Cotr instance

viscid.cotr.cotr_transform(date_time, from_system, to_system, vec, notilt1967=True)

Transform a vector from one system to another

Parameters:

• date_time (str, datetime64) – datetime of transformation
• from_system (str) – abbreviation of crd_system
• to_system (str) – abbreviation of crd_system
• vec (ndarray) – should have shape (3, ) or (N, 3) to transform N vectors at once
• notilt1967 (bool) – is 1 Jan 1967 the special notilt time?

Returns:

vec in new crd system

Return type:

ndarray

viscid.cotr.get_dipole_moment(date_time, crd_system='gse', strength=32707.5292732387, notilt1967=True)

Get Earth’s dipole moment at datetime

Parameters:

• date_time (str, datetime64) – datetime to get dipole
• crd_system (str) – crd_system of resulting moment
• strength (float) – magnitude of dipole moment
• notilt1967 (bool) – is 1 Jan 1967 the special notilt time?

Returns:

3-vector of dipole tilt in crd_system

Return type:

ndarray

viscid.cotr.get_dipole_moment_ang(dip_tilt=0.0, dip_gsm=0.0, crd_system='gse', strength=32707.5292732387)

Get dipole moment from tilt and gsm angles

Parameters:

• dip_tilt (float) – if given, override the dipole tilt angle (mu degrees, positive points north pole sunward)
• dip_gsm (float) – if given, override the dipole gse->gsm angle (psi in degrees, positive points north pole duskward)
• crd_system (str) – coordinate system of the result
• strength (float) – magnitude of the dipole moment

Returns:

dipole moment [mx, my, mz]

Return type:

ndarray

viscid.cotr.get_dipole_angles(date_time, notilt1967=True)

Get rotation angles between dipole and GSE

Parameters:

• date_time (str, datetime64) – datetime to get dipole
• notilt1967 (bool) – is 1 Jan 1967 the special notilt time?

Note

psi refers to the angle between GSE and GSM, while mu is the dipole tilt angle.

Returns:

The SM -> GSE, rotation is <-mu, Y> * <-psi, X>

which is equivalent to <psi, X> * <mu, Y>

Return type:(tilt, gsm)

viscid.cotr.dipole_moment2cotr(m, crd_system='gse')

Turn dipole moment vector into a Cotr instance

Parameters:

• m (sequence) – dipole moment as sequence with length 3
• crd_system (str) – crd_system of m

Returns:

Cotr instance