irfpy.util.fields

Scalar field or 3-D vector field in 2-D or 3-D spaces

Code author: Yoshifumi Futaana

In physics, any value (scalar or vector) assigned at a specific position is called as the field. Mathematically, this module represents mapping from \(R^3 oR\) or \(R^3 o R^3\). Scalar field and 3D vector field are implemented in this module.

You can find a tutorial here.

Using a scalar field

A scalar field, \(f(x,y,z)=p\) are represented by ScalarField3D class and its derivatives. For example, you may find a FunctionScalarField3D class, with which the value p is obtained from a given function. A typical example is to represent gravity potential field.

Another example use case is a field represented by a gridded system. For example, you can use MHD simulation results with this. A cartesian field is implemented as GridScalarField3D.

Using a vector field

Sample vector field implementations are

• UniformVector3dField

• DipoleField

• More will come / Please contribute.

You can easily create the field that you need.

For example, if you have gridded field data, you may use the following classes.

• GridScalarField3D

• GridVector3dField3D

If you have functions that represent the field instead, you may use the following classes.

• FunctionScalarField3D

• FunctionVector3dField3D — Not yet implemented

Base classes of scalar field in 3-D space in ScalarField3D and 3-D vector field in 3-D space in :class`Vector3dField3D`.

class irfpy.util.fields.ScalarField3D(extrapolation='raise')

Bases: object

Baseclass of scalar field.

This is also callable. See __call__().

__call__(vector)

Return the value.

Equivalent to get(), but the argument is a array.

abstract get_field(x, y, z)

Override this method for implementation.

get(x, y, z)

User interface to get the field.

User interface to get the field. Domain check is done, and then call and return get_field().

Parameters:

• x – x value. Scalar.

• y – y value. Scalar.

• z – z value. Scalar.

Raises:

ValueError – if the given coordinate is not in domain.

indomain(x, y, z)

Return true if the given x, y, z is in the domain

This is the user interface.

abstract xdomain()

Return (x0, x1), xdomain.

abstract ydomain()

Return (y0, y1), ydomain.

abstract zdomain()

Return (z0, z1), zdomain.

set_extrapolation_method(method)

Set the extrapolation method.

Parameters:

method – An extrapolation method. Either “raise”, “nan”, or “flat”.

extrapolate(x, y, z)

class irfpy.util.fields.ScalarField3DInfty(sf1)

Bases: ScalarField3D

By wrapping other scalar field, the field is extended to infinity domain.

Warning

This class is deprecated. Use extrapolate="flat" keyword when creating scalar field.

The finite domain field is extended to the infitity domain, by taking the nearest point at the edge of the original domain.

Suppose the following 3D gridded field as a baseline.

\[f(x, y, z) = 1/|x - 0.5| + 1/(y - 0.7)^2 + 1/(z + 0.3)^4\]

>>> xlist = np.array([-2, 1, 0, 1, 2])
>>> ylist = np.array([-2, 1, 0, 1, 2])
>>> zlist = np.array([-2, 1, 0, 1, 2])
>>> xlist3d, ylist3d, zlist3d = np.meshgrid(xlist, ylist, zlist, indexing='ij')
>>> vallist = 1 / np.abs(xlist3d - 0.5) + 1 / ((ylist3d - 0.7) ** 2) + 1 / ((zlist3d + 0.3) ** 4)
>>> gsf = GridScalarField3D(xlist, ylist, zlist, vallist)

The domain of gsf is within the grid definition, namely, -2 to 2 in each axis.

>>> print(gsf.zdomain())
(-2, 2)

If you wants the value at (x, y, z) = (0, 0, 3), an exception is occurred, since the value is out of the domain defined.

By wrapping the gsf by ScalarField3DInfty, the domain is extended to the infinity.

>>> gsf_inf = ScalarField3DInfty(gsf)

gsf_inf can return the value at the point outside of the domain of gsf.

>>> print(gsf_inf.get(0, 0, 3))
4.076550904380177

This is the same value as the point at (0, 0, 2) for gsf, i.e. the nearest point in the doimain of gsf.

>>> print(gsf.get(0, 0, 2))
4.076550904380177

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

get_field(x, y, z)

Override this method for implementation.

class irfpy.util.fields.ScalarField3DInftyNan(scalarfield3d)

Bases: ScalarField3D

By wrapping other scalar field, the new field is extended to infinity domain.

Warning

This class is deprecated. Use extrapolate="flat" keyword when creating scalar field.

The finate domain field is extended to the infinity domain, making the user’s program not crash by ValueError.

The difference with ScalarField3DInfty is that this class returns np.nan if the field value is not defined.

Example

Let us make a scalar field defined by a grid value.

>>> xg = [0, 1, 2, 3]
>>> yg = [0, 3, 6]
>>> zg = [0, 10]
>>> field_value = np.arange(24).reshape([4, 3, 2])    # 4x3x2 array to represent a scalar field
>>> scalar_field = GridScalarField3D(xg, yg, zg, field_value)

In this case, the field of (1, 3, 5) should return 8.5.

>>> print(scalar_field([1, 3, 5]))
8.5

If you give the out of range place, a ValueError is raised.

>>> print(scalar_field([-1, -1, -1]))   
Traceback (most recent call last):
...
ValueError: x, y, z = (-1, -1, -1) is not in domain

To avoid ValueError to stop, you may use the following statements.

>>> scalar_field = ScalarField3DInftyNan(scalar_field)

Then, nan is returned.

>>> print(scalar_field([-1, -1, -1]))
nan

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

get_field(x, y, z)

Override this method for implementation.

irfpy.util.fields.add(sf1, sf2)

Return the ScalarField3D added.

Return type:

ScalarField3D

irfpy.util.fields.neg(sf1)

Return the negated ScalarField3D.

Return type:

ScalarField3D

irfpy.util.fields.sub(sf1, sf2)

Return sf1-sf2

Return type:

ScalarField3D

irfpy.util.fields.mul(k_or_sf, sf1)

Multiply by scalar or field.

Parameters:

• k_or_sf – Scalar or ScalarField3D.

• sf1 – ScalarField3D.

Return type:

ScalarField3D

irfpy.util.fields.power(sf1, powerindex)

Returns sf1 ** powerindex

Parameters:

sf1 – ScalarField3D.

Return type:

ScalarField3D

irfpy.util.fields.sqrt(sf1)

Take square root from the given scalar field.

Parameters:

sf1 – Scalar field.

Returns:

sqrt(sf1)

Return type:

ScalarField3D

Sample follows

>>> sf = GridScalarField3D([-2, 2], [-2, 2], [-1, 1], np.arange(8.).reshape([2, 2, 2]))
>>> print(sf([0, 0, 0]))
3.5
>>> sqrt_sf = sqrt(sf)
>>> print('{:.3f}'.format(sqrt_sf([0, 0, 0])))
1.871

irfpy.util.fields.div(sf1, k_or_sf)

Divide by scalar or field.

See mul().

Return type:

ScalarField3D

class irfpy.util.fields.IntegrateScalarField(scalarfield)

Bases: object

An integrator of scalar field along a segment. (Object-oriented version)

This class is to integrate a scalar value along a segment. Use can use integrate_trapz() function as well.

USAGE

You have to first make object of ScalarField3D.

>>> func = lambda x, y, z: z + np.sqrt(x*x + y*y)
>>> scalarfield = FunctionScalarField3D(func)

Then, instance IntegrateScalarField object.

>>> integfunc = IntegrateScalarField(scalarfield)

Below is the segment you want to calculate the integral along. This should be (3, N) array.

>>> segment = np.array([[0., 0, 0], [0, 0.5, 0], [0, 0.7, 0.2], [0.2, 0.7, 0.4],
...                 [0.2, 0.7, 0.5], [0.4, 0.7, 0.7], [1.0, 0.7, 0.7], [1.0, 1.0, 1.0]]).T

Then, to run trapezoid integral call trapz().

>>> print('%.4f' % integfunc.trapz(segment))
3.0619

You can use Simpson’s integral. Call simps().

>>> print('%.4f' % integfunc.simps(segment))
3.0539

trapz(segment)

Return the trapezoid integral.

Parameters:

segment – np.array of vortex with (3, N) shape

simps(segment)

Return the simpson’s integral.

Parameters:

segment – np.array of vortex with (3, N) shape

irfpy.util.fields.integrate_trapz(sf, segment)

Integrate the given scalar field along the given segment.

\(\int_{L} f(x(l), y(l), z(l)) dl\)

Parameters:

• sf – ScalarField3D

• segment – A np.array with dimension of (3, N).

Integral is made using scipy.integrate.trapz. The simplest algorithm of trapezoidal integration provides the most robust solution. For each given vertex in segment argument, the field value is calculated as y, and considers the L2 distance as x.

See IntegrateScalarField and its method of IntegrateScalarField.trapz().

A very simple example follows. Set up uniform field, and integrate over a certain segment.

>>> xlist = np.array([0, 4])
>>> ylist = np.array([0, 4])
>>> zlist = np.array([0, 4])
>>> vallist = np.ones([2, 2, 2]) * 1.5   # Uniform field with domain 0<4.
>>> scalarfield = GridScalarField3D(xlist, ylist, zlist, vallist)
>>> segment = np.array([[0., 0, 0], [0, 1, 0], [0, 1, 1], [0, 1, 2], [1, 1, 2],
...             [1, 1, 3], [1, 2, 3], [1, 3, 3], [2, 3, 3], [2, 4, 3], [2, 4, 4],]).T
>>> print(integrate_trapz(scalarfield, segment))
15.0

Another example is to use function type field.

>>> func = lambda x, y, z: z + np.sqrt(x*x + y*y)
>>> scalarfield = FunctionScalarField3D(func)

Below shows (non-practical) example of 1 segment integral. func(0, 0, 0) is 0 and func(0.1, 0.1, 0.1) is 0.2414. The distance is 0 (by definition) for the first point and 0.1732 for (0.1, 0.1, 0.1). Therefore the integral should be 0.2414 * 0.1732 / 2 = 0.0209.

>>> segment = np.array([[0., 0, 0], [0.1, 0.1, 0.1]]).T
>>> print('%.4f' % integrate_trapz(scalarfield, segment))
0.0209

More practical sample follows.

>>> segment = np.array([[0., 0, 0], [0, 0.5, 0], [0, 0.7, 0.2], [0.2, 0.7, 0.4],
...                 [0.2, 0.7, 0.5], [0.4, 0.7, 0.7], [1.0, 0.7, 0.7], [1.0, 1.0, 1.0]]).T
>>> print('%.4f' % integrate_trapz(scalarfield, segment))
3.0619

class irfpy.util.fields.GridScalarField3D(xlist, ylist, zlist, vallist, masked_invalid=False, interpolation='linear', extrapolation='raise')

Bases: ScalarField3D

Scalar field with gridded data.

An implementation of scalar field in 3D. The cartesian grid (xlist, ylist, and zlist) and the corresponding scalar value (:attr:`vallist) is treated as a field.

The get() method will return the interpolated data by trilinear algorithm.

>>> xlist = np.array([0, 1, 2, 3, 4])
>>> ylist = np.array([0, 1, 2, 3])
>>> zlist = np.array([0, 1, 2])
>>> vallist = np.arange(60).reshape([5, 4, 3])
>>> gsf = GridScalarField3D(xlist, ylist, zlist, vallist)
>>> print(gsf.get(1, 2, 0))
18.0

bin_index() returns the indexes.

>>> gsf.bin_index(1.5, 3, 1)
(1, 2, 1)

Note that the y index iy=2, not 3. This ensures the ylist[iy+1] can always be referred to.

Parameters:

• xlist – List of x coordinates. All elements should be increasing, i.e. xlist[i] < xlist[i+1] for all i.

• ylist – List of y coordinates.

• zlist – List of z coordinates.

• vallist – List of value for the field. np.array of (Nx, Ny, Nz) shape.

• masked_invalid – If True, the vallist is initialized by np.ma.maked_invalid, instead of np.array. Not applicable to the coordinates.

• interpolation – Either linear or nearest. Refer to set_interpolation_method().

• extrapolation – Either raise, nan or flat

set_interpolation_method(method)

Set the interpolation method.

Method:

Choose the interpolation method. Currently, “linear” and “nearest” is supported.

• Linear returns the linearly interpolated values (interpolate_trilinear())

• Nearest returns the nearest neighbor values (interpolate_nearest_neighbor())

get_field(x, y, z)

Override this method for implementation.

interpolate_trilinear(x, y, z)

Return the interpolated value. No domain check is done.

Parameters:

• x – X value

• y – Y value

• z – Z value

Returns:

the interpolated field value

interpolate_nearest_neighbor(x, y, z)

Return the interpolated value using nearest neighbor method.

Parameters:

• x – X value

• y – Y value

• z – Z value

Returns:

the interpolated field value

>>> xlist = np.array([0, 1, 2, 3, 4])
>>> ylist = np.array([0, 1, 2, 3])
>>> zlist = np.array([0, 1, 2])
>>> vallist = np.arange(60.).reshape([5, 4, 3])
>>> gsf = GridScalarField3D(xlist, ylist, zlist, vallist)
>>> print(gsf.interpolate_nearest_neighbor(1, 2, 0))
18.0
>>> print(gsf.interpolate_nearest_neighbor(1, 2, 0.4))
18.0
>>> print(gsf.interpolate_nearest_neighbor(1, 2, 0.6))
19.0
>>> print(gsf.interpolate_nearest_neighbor(1, 1.6, 0))
18.0
>>> print(gsf.interpolate_nearest_neighbor(1, 1.4, 0))
15.0
>>> print(gsf.interpolate_nearest_neighbor(1.4, 2, 0))
18.0
>>> print(gsf.interpolate_nearest_neighbor(1.6, 2, 0))
30.0

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

bin_index(x, y, z)

Return the (left) index of the given x, y, z

For each component, i value satisfying xlist[i] <= x < xlist[i+1] is returned.

If x is the upper most value, the index i = len(xlist)-2 is returned. This is to ensure that the returned index is “left” bin so that i+1 can be always used to refer to for further calculation.

class irfpy.util.fields.GridScalarFieldCacheFiled(xlist, ylist, zlist, vallist, cache_filename, keep_cache=False, interpolation='linear', extrapolation='raise')

Bases: GridScalarField3D

Cache file support of GridScalarField3D

Parameters:

• xlist – List of x coordinates. All elements should be increasing, i.e. xlist[i] < xlist[i+1] for all i.

• ylist – List of y coordinates.

• zlist – List of z coordinates.

• vallist – List of value for the field. np.array of (Nx, Ny, Nz) shape. If the file cache_filename exists, this argument is ignored.

• cache_filename – The vallist saved by numpy.save command.

• keep_cache – If the cache file is to be saved, this option should be set True.

unload()

load()

get_field(x, y, z)

Override this method for implementation.

interpolate_trilinear(x, y, z)

Return the interpolated value. No domain check is done.

Parameters:

• x – X value

• y – Y value

• z – Z value

Returns:

the interpolated field value

irfpy.util.fields.grid_gradient(gridscalarfield)

Gradient of the scalar field (del-del-x, del-del-y, del-del-z) is calculated.

The gradient of the scalar field is calculated, using numpy.gradient() for GridScalarField3D object. According to the np.gradient document:

The gradient is computed using second order accurate central differences
in the interior and either first differences or second order accurate
one-sides (forward or backwards) differences at the boundaries.

Parameters:

gridscalarfield – A GridScalarField3D object, g(x, y, z)

Returns:

Three new GridVector3dField3D representing \(\nabla g(x, y, z)\)

Return type:

tuple(GridVector3dField3D, GridVector3dField3D, GridVector3dField3D)

An original GridScalarField3D is first instanced. The field is \(g(x,y,z)=x^4\cdot y^3 \cdot z^2\).

>>> xlist = np.array([0, 1.9, 2, 2.1])
>>> ylist = np.array([0, 1, 3.9, 4, 4.1])
>>> zlist = np.array([0, 1, 2, 7.9, 8, 8.1])
>>> xlist3d, ylist3d, zlist3d = np.meshgrid(xlist, ylist, zlist, indexing='ij')
>>> val = (xlist3d ** 4) * (ylist3d ** 3) * (zlist3d ** 2)
>>> gsf = GridScalarField3D(xlist, ylist, zlist, val)

The gradient of g is taken by grid_gradient() function, creating GridVector3dField3D object. Mathematically, the gradient is

\[\begin{split}gx(x, y, z) = 4x^3 \cdot y^3 \cdot z^2, \\ gy(x, y, z) = x^4 \cdot 3y^2 \cdot z^2, \\ gz(x, y, z) = x^4 \cdot y^3 \cdot 2z.\end{split}\]

>>> grad_gsf = grid_gradient(gsf)
>>> vec = grad_gsf.get(2, 4, 8)
>>> print(vec)      
[131399.68   49162.24   16384.  ]

class irfpy.util.fields.FunctionScalarField3D(valfunc, xdomain=None, ydomain=None, zdomain=None, extrapolation='raise')

Bases: ScalarField3D

Scalar field defined by a function

>>> func = lambda x, y, z: x + np.cos(y) + 1. / (np.exp(z) ** 2)
>>> scf = FunctionScalarField3D(func)
>>> print(scf.get(0, 0, 0))
2.0

These two statements are identical.

>>> print('%.3f' % scf.get(-10, np.pi, 1.5))
-10.950
>>> print('%.3f' % scf([-10, np.pi, 1.5]))
-10.950

This class is a kind of “wrapper” of scalar field function to behave as ScalarField3D abstract class.

Initializer

Parameters:

• valfunc (A callable.) – A function (callable) that eats 3 float to return 1 float. v=f(x, y, z) where x, y, z is the position and v is the value of the field.

• xdomain – Domain defined by the valfunc. [x0, x1]. If not given or None is given, [-np.inf, +np.inf] is given. Both are inclusive.

• ydomain – Domain defined by the valfunc. [y0, y1].

• zdomain – Domain defined by the valfunc. [z0, z1].

get_field(x, y, z)

Override this method for implementation.

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

irfpy.util.fields.rebin_scalarfield(xlist, ylist, zlist, scalarfield, interpolation='linear')

Given scalar field is re-binned and a new scalar field (GridScalarField3D) is created.

Parameters:

• xlist – List of x coordinates. All elements should be increasing, i.e. xlist[i] < xlist[i+1] for all i.

• ylist – List of y coordinates.

• zlist – List of z coordinates.

• scalarfield – Scalar field

Returns:

A corresponding gridded scalar field.

Return type:

GridScalarField3D

The given scalar field is converted to GridScalarField3D using the given coordinate system. Below, a usual work flow is shown. First, prepare the original GridScalarField3D object, here we took a FunctionScalarField3D as an example. This object is converted to the GridScalarField3D by using this function.

As an example, one may first get object of GravityPotential. This is a subclass of FunctionScalarField3D.

>>> gPot = GravityPotential(1e12, center=(1e-8, 0, 0))   # Gravity potential produced by 1e12 kg mass at the center (slightly off to avoid infinity.).

Note that here an enormously high mass of 1e12 kg is selected to understand the gravity field around 1 m. To avoid the “infinity gravity at the point of mass, Very small (1e-8) offset was introduced.

Here, we see some potentials close-by.

>>> print('{:.2e}'.format(gPot.get(0.5, 0, 0)))   # At 0.5 m distance
-1.33e+02
>>> print('{:.2e}'.format(gPot.get(1, 0, 0)))   # At 1 m distance
-6.67e+01
>>> print('{:.2e}'.format(gPot.get(1.5, 0, 0)))   # At 1.5 m distance
-4.45e+01
>>> print('{:.2e}'.format(gPot.get(2, 0, 0)))   # At 2 m distance
-3.34e+01

Now, we get the GridScalarField3D.

>>> xlist = [-1, 0, 1, 2, 3]
>>> ylist = [-1, 0, 1, 2]
>>> zlist = [-1, 0, 1]
>>> gPotGrid = rebin_scalarfield(xlist, ylist, zlist, gPot)

>>> isinstance(gPotGrid, GridScalarField3D)
True

Now we get it! See the values on grid.

>>> print('{:.2e}'.format(gPotGrid.get(1, 0, 0)))   # At 1 m distance
-6.67e+01
>>> print('{:.2e}'.format(gPotGrid.get(2, 0, 0)))   # At 2 m distance
-3.34e+01

The value not on the exact grid should be just an interpolation >>> print(‘{:.2e}’.format(gPotGrid.get(1.5, 0, 0))) # At 1.5 m distance. The original function has returned -4.45e+01, but… -5.01e+01

class irfpy.util.fields.GravityPotential(mass, center=None)

Bases: FunctionScalarField3D

Implementation of gravity potential

Gravity potential implementated as FunctionScalarField3D.

\[\phi(x, y, z) = - GM / r\]

Parameters:

mass – Mass in kg.

Gravity field created from the Sun.

>>> mass_sun = 1.989e30
>>> gravity_potential_sun = GravityPotential(mass_sun)
>>> au = 1.49597871e11     # 1 AU in meter
>>> print('{:.2e}'.format(gravity_potential_sun.get(au, 0, 0)))
-8.87e+08
>>> print('{:.2e}'.format(gravity_potential_sun.get(0, au, 0)))
-8.87e+08

Gravity field created by the Earth

>>> mass_earth = 5.972e24
>>> gravity_potential_earth = GravityPotential(mass_earth, center=[au, 0, 0])   # earth at (1Au, 0, 0)
>>> print('{:.2e}'.format(gravity_potential_earth.get(0, 0, 0)))   # Earth's gravity potential at the Sun.
-2.66e+03

Gravity field created by the Sun & Earth

>>> gravity_potential_sun_earth = add(gravity_potential_sun, gravity_potential_earth)
>>> lagrange1 = au - 1.5e9
>>> print('{:.2e}'.format(gravity_potential_sun_earth.get(lagrange1, 0, 0)))
-8.97e+08

G = 6.67408e-11

irfpy.util.fields.xyz2rtp(x, y, z)

Convert from Cartesian to the polar frame.

Polar frame is defined as follows. r is radial, \(\theta\) is the angle from z-axis, and \(\phi\) is the angle in x-y plane from x-axis. The expression of the angle is in radians. It is very natural.

The formuation is

\[x = r \sin\theta \cos\phi y = r \sin\theta \sin\phi z = r \cos\theta\]

Parameters:

• x – x coordinates

• y – y coordinates

• z – z coordinates

Returns:

A tuple (r, \(\theta\), \(\phi\)).

>>> r, t, p = xyz2rtp(10, 10, 10)
>>> print('{:.2f}'.format(r))
17.32
>>> print('{:.2f}'.format(np.rad2deg(t)))
54.74
>>> print('{:.2f}'.format(np.rad2deg(p)))
45.00

class irfpy.util.fields.Vector3dField3D(extrapolation='raise')

Bases: object

Abstract base class of 3D vector field in 3-D space.

This is also callable.

__call__(vector)

Return the value.

Equivalent to get(), but the argument is a array.

abstract get_field(x, y, z)

Override this method for implementation.

get(x, y, z)

User interface to get the field.

User interface to get the field. Domain check is done, and then call and return get_field().

Parameters:

• x – x value. Scalar.

• y – y value. Scalar.

• z – z value. Scalar.

Returns:

(vx, vy, vz), the vector of the field. (3,) shaped numpy array

Raises ValueError, if the given coordinate is not in domain.

set_extrapolation_method(method)

Set the extrapolation method.

Parameters:

method – An extrapolation method. Either “raise”, “nan”, or “flat”.

indomain(x, y, z)

Return true if the given x, y, z is in the domain

This is the user interface.

abstract xdomain()

Return (x0, x1), xdomain.

abstract ydomain()

Return (y0, y1), ydomain.

abstract zdomain()

Return (z0, z1), zdomain.

extrapolate(x, y, z)

get_rtp(x, y, z)

Return the vector in RTP coordinates.

In the RTP coordinates, the vector is expressed by the r:radial, t:southward, p:eastward coordinates. It is a local coodinates. Mathematically it is expressed as follows.

\[ \begin{align}\begin{aligned}\begin{split}\hat{r} = \left(\begin{array}{c} \sin\theta\cos\phi \\ \sin\theta\sin\phi \\ \cos\theta \end{array}\right) \\\end{split}\\\begin{split} \hat{\theta} = \left(\begin{array}{c} \cos\theta\cos\phi \\ \cos\theta\sin\phi \\ -\sin\theta \end{array}\right) \\\end{split}\\\begin{split}\hat{\phi} = \left(\begin{array}{c} -\sin\phi \\ \cos\phi \\ 0 \end{array}\right)\end{split}\end{aligned}\end{align} \]

Parameters:

• x – X value

• y – Y value

• z – Z value

Returns:

(vr, vt, vp) in the RTP coordinates.

class irfpy.util.fields.Vector3dField3DInftyNan(vectorfield3d)

Bases: Vector3dField3D

By wrapping other vector field, the new vector field is extended to infinity domain.

The finate domain field is extended to the infinity domain, making the user’s program not crash by ValueError.

This provides the similar workflow as ScalarField3DInftyNan, but for a 3-D vector field.

Example

Without details, we make a vector field from a gridded data.

>>> xlist = [-1, 0, 1, 2, 3]
>>> ylist = [0, 1.5, 3, 4.5]
>>> zlist = [0.1, 0.2, 0.3]

>>> vx = np.ones([5, 4, 3])   # Should be (5, 4, 3) shape
>>> vy = np.arange(60).reshape([5, 4, 3])   # Should be (5, 4, 3) shape
>>> vz = np.arange(30, 90).reshape([5, 4, 3])   # Should be (5, 4, 3) shape

>>> gvf = GridVector3dField3D(xlist, ylist, zlist, vx, vy, vz)

The vector field has domain [-1, 3], [0, 4.5], [0.1, 0.3], but with wrapping by

Vector3dField3DInftyNan, np.nan is returned for the points outside the domain.

>>> gvf = Vector3dField3DInftyNan(gvf)
>>> print(gvf([-100, -100, -100]))    # The given point is out of the domain.
nan

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

get_field(x, y, z)

Override this method for implementation.

irfpy.util.fields.abs_vf(vf1)

Return the scalar field produced from a vector field.

Parameters:

vf1 (Vector3dField3D) – A vector field.

Returns:

A scalar filed.

Return type:

ScalarField3D

As an example, the following grid field is produced.

>>> x = y = z = np.arange(5)
>>> vx, vy, vz = np.meshgrid(x, y, z, indexing='ij')

This is 5x5x5 grid, with vx[i, *, *] = i, vy[*, j, *] = j, and vz[*, *, k] = k.

Thus, the absolute value of this field is \(|v[i, j, k]| = \sqrt(i^2 + j^2 + k^2)\).

>>> v = GridVector3dField3D(x, y, z, vx, vy, vz)
>>> vabs = abs_vf(v)
>>> print(vabs([0, 0, 0]))
0.0
>>> print('{:.3f}'.format(vabs([1, 2, 3])))   # sqrt(14)
3.742

class irfpy.util.fields.FunctionVector3dField3D

Bases: Vector3dField3D

class irfpy.util.fields.GridVector3dField3D(xlist, ylist, zlist, vxlist, vylist, vzlist, interpolation='linear', extrapolation='raise')

Bases: Vector3dField3D

Grid version of 3D vector field.

>>> xlist = [-1, 0, 1, 2, 3]
>>> ylist = [0, 1.5, 3, 4.5]
>>> zlist = [0.1, 0.2, 0.3]

>>> vx = np.ones([5, 4, 3])   # Should be (5, 4, 3) shape
>>> vy = np.arange(60).reshape([5, 4, 3])   # Should be (5, 4, 3) shape
>>> vz = np.arange(30, 90).reshape([5, 4, 3])   # Should be (5, 4, 3) shape

>>> gvf = GridVector3dField3D(xlist, ylist, zlist, vx, vy, vz)
>>> vec = gvf.get(3, 3, 0.2)
>>> print(vec[0], vec[1], vec[2])
1.0 55.0 85.0

You can also get interpolated value.

>>> vec = gvf.get(3, 3., 0.225)
>>> print(vec[0], vec[1], vec[2])
1.0 55.25 85.25

Parameters:

• xlist – X coordinates. Shape of (Nx,). See also GridScalarField3D.

• ylist – Y coordinates. Shape of (Ny,)

• zlist – Z coordinates. Shape of (Nz,)

• vxlist – x component of vector field. Shape of (Nx, Ny, Nz).

• vylist – y component of vector field. Shape of (Nx, Ny, Nz).

• vzlist – z component of vector field. Shape of (Nx, Ny, Nz).

get_field(x, y, z)

Override this method for implementation.

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

set_interpolation_method(method)

irfpy.util.fields.grid_divergence(grid_vector3dfield)

Return a numerically calculated divergence of the given vector field.

Parameters:

grid_vector3dfield – A GridVector3dField3D object

Returns:

A GridScalarField3D object

Return type:

GridScalarField3D

First, create a vector field as example.

>>> xlist = [9.8, 9.9, 10.0, 10.1, 10.2]    # Coordinate of grid
>>> ylist = [6.8, 6.9, 7.0, 7.1, 7.2]    # Coordinate of grid
>>> zlist = [-10.2, -10.1, -10, -9.9, -9.8]   # Coordinate of grid
>>> xlist3d, ylist3d, zlist3d = np.meshgrid(xlist, ylist, zlist, indexing='ij')   # Make corrdinates to 3D
>>> valx = 3 * xlist3d * ylist3d * zlist3d      # fx(x, y, z)
>>> valy = xlist3d ** 2 * ylist3d     # fy(x, y, z)
>>> valz = 2 * ylist3d ** 2 * zlist3d     # fz(x, y, z), where f=(fx, fy, fz)
>>> vectorfield = GridVector3dField3D(xlist, ylist, zlist, valx, valy, valz)

The original field is

\[\vec{F(x,y,z)} = (3xyz, x^2 y, 2 y^2 z)\]

The divergence field is

\[\nabla\cdot\vec{F(x, y, z)} = 3yz + x^2 + 2 y^2\]

>>> div = grid_divergence(vectorfield)
>>> print('{:.3f}'.format(div.get(10, 7, -10)))     
-12.000

irfpy.util.fields.grid_rotation(grid_vector3dfield)

Return a numerically calculated rotation of the given vector field.

param grid_vector3dfield:

A GridVector3dField3D object

return:

A GridVector3dField3D object.

rtype:

GridVector3dField3D

First, create a vector field as example.

>>> xlist = [9.8, 9.9, 10.0, 10.1, 10.2]    # Coordinate of grid
>>> ylist = [6.8, 6.9, 7.0, 7.1, 7.2]    # Coordinate of grid
>>> zlist = [-10.2, -10.1, -10, -9.9, -9.8]   # Coordinate of grid
>>> xlist3d, ylist3d, zlist3d = np.meshgrid(xlist, ylist, zlist, indexing='ij')   # Make corrdinates to 3D
>>> valx = 3 * xlist3d * ylist3d * zlist3d      # fx(x, y, z)
>>> valy = xlist3d ** 2 * ylist3d     # fy(x, y, z)
>>> valz = 2 * ylist3d ** 2 * zlist3d     # fz(x, y, z), where f=(fx, fy, fz)
>>> vectorfield = GridVector3dField3D(xlist, ylist, zlist, valx, valy, valz)

The original field is

\[\]

ec{F(x,y,z)} = (3xyz, x^2 y, 2 y^2 z)

The divergence field is

\[\]

abla imes ec{F(x, y, z)} = (4yz, 3xy, 2xy - 3xz)

>>> rot = grid_rotation(vectorfield)
>>> print(rot.get(10, 7, -10))
[-280.  210.  440.]

class irfpy.util.fields.UniformVector3dField(xval, yval, zval, xdomain=None, ydomain=None, zdomain=None)

Bases: Vector3dField3D

An implementation of the uniform field

Parameters:

• xval – X value of the field.

• yval – Y value of the field.

• zval – Z value of the field.

• xdomain – If one wants to specify the domain, use (xyz)domain. List or tuple with two elements, (lower, upper). Default is infinity, namely, (-np.inf, np.inf).

>>> bfield = UniformVector3dField(1e-9, -1e-9, 0)
>>> print(bfield.get(0, 0, 0))    # At the origin.   
[ 1.e-09 -1.e-09 0.e+00]
>>> print(bfield.get(1e30, 1e30, -1e30))    # At the large distance   
[ 1.e-09 -1.e-09 0.e+00]

get_field(x, y, z)

Override this method for implementation.

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

irfpy.util.fields.dot(vf1, vf2)

Produce a ScalarField3D from two Vector3dField3D by taking the dot (inner product)

Parameters:

• vf1 – A Vector3dField3D

• vf2 – A Vector3dField3D

Returns:

A ScalarField3D object :math:`vf1 cdot vf2.

Return type:

ScalarField3D

A simple example is to obtain the magnitude of a vector field. Assume the magnetic field of (60, 80, 100) nT:

>>> bfield = UniformVector3dField(60e-9, 80e-9, 100e-9)
>>> print(bfield([5, 3., 0]))  # Value at any position   
[6.e-08 8.e-08 1.e-07]

The square magnitude of the field, |B|^2, is obtained by the dot product:

>>> squ_mag = dot(bfield, bfield)
>>> print('{:.4e}'.format(squ_mag([5, 3., 0])))   # Value at any position
2.0000e-14

irfpy.util.fields.cross(vf1, vf2)

Produce a Vector3dField3D from two Vector3dField3D by taking the cross (outer product)

Parameters:

• vf1 – A Vector3dField3D

• vf2 – A Vector3dField3D

Returns:

A Vector3dField3D object :math:`vf1 times vf2.

Return type:

Vector3dField3D

A simple example is to obtain the convection electric field. Assume the magnetic field as uniform

>>> bfield = UniformVector3dField(0, 0, 100e-9)
>>> print(bfield([0, 0, 0]))  # Value at origin   
[0.e+00 0.e+00 1.e-07]

The velocity vector is also assumed to be uniform.

>>> vfield = UniformVector3dField(-400e3, 0, 0)
>>> print(vfield([0, 0, 0]))  # Value at origin
[-400000.       0.       0.]

Then, the convection electric field is \(E=-v\timesB\). (To make the code simpler, use \(E=B\timsv\).) In this case, the convection field must be always (0, -0.04, 0).

>>> efield = cross(bfield, vfield)
>>> print(efield([0, 0, 0]))  # Value at origin
[ 0.   -0.04  0.  ]

class irfpy.util.fields.DipoleField(moment_vector, coefficient=1.0)

Bases: Vector3dField3D

Implementation of the dipole field.

A dipole field is produced.

Parameters:

• moment_vector – A moment vector, like \(\vec{m}\).

• coefficient – Coefficient. For magnetic moment, MagneticCoefficient (\(\frac{\mu_0}{4\pi}\)) should be used. For electric moment, ElectricCoefficient (\(\frac{1}{4\pi\epsiron_0}\)) should be used. Default is 1.

For example, geomagnetic dipole moment 8e22 (A/m^2) is instanced as (no transpose, no tilt considered here)

>>> geomag = DipoleField([0, 0, -8e22], coefficient=DipoleField.MagneticCoefficient)
>>> print(geomag([0, 0, 6378e3]))  # At the (magnetic) pole     
[-0.00000000e+00  -0.00000000e+00  -6.16689335e-05]
>>> print(geomag([6378e3, 0, 0]))  # At the equator   
[0.00000000e+00   0.00000000e+00   3.08344668e-05]
>>> print(geomag([4000e3, 0, 4000e3]))  # At 45 deg north   
[-6.62912607e-05  -0.00000000e+00  -2.20970869e-05]

You get ~60000 nT at the pole and ~30000 nT at the equator.

You can also get the geomagnetic field in the local (RTP) coordinates.

>>> print('{b[0]:.2f} {b[1]:.2e} {b[2]:.2f}'.format(b=geomag.get_rtp(6378e3, 0, 0)))   # At the equator
0.00 -3.08e-05 0.00

>>> print('{b[0]:.2e} {b[1]:.2e} {b[2]:.2e}'.format(b=geomag.get_rtp(4000e3, 0, 4000e3)))   # At the 45 deg north
-6.25e-05 -3.13e-05 0.00e+00

MagneticCoefficient = 1e-07

Coefficient for the magnetic dipole.

ElectricCoefficient = 8987551787.997911

Coefficient for the electric dipole.

get_field(x, y, z)

Override this method for implementation.

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

zdomain()

Return (z0, z1), zdomain.

class irfpy.util.fields.ScalarField2D(extrapolation='raise')

Bases: object

Baseclass of scalar field.

This is also callable. See __call__().

ScalarField3D.__call__(vector)

Return the value.

Equivalent to get(), but the argument is a array.

abstract get_field(x, y)

Override this method for implementation.

get(x, y)

User interface to get the field.

User interface to get the field. Domain check is done, and then call and return get_field().

Parameters:

• x – x value. Scalar.

• y – y value. Scalar.

Raises:

ValueError – if the given coordinate is not in domain.

indomain(x, y)

Return true if the given x, y, z is in the domain

This is the user interface.

abstract xdomain()

Return (x0, x1), xdomain.

abstract ydomain()

Return (y0, y1), ydomain.

set_extrapolate_method(method)

Set the extrapolation method.

Parameters:

method – An extrapolation method. Either “riase”, “nan”, or “flat”

extrapolate(x, y)

class irfpy.util.fields.GridScalarField2D(xlist, ylist, vallist, interpolation='linear', extrapolation='raise')

Bases: ScalarField2D

Scalar field with gridded data.

The scalar field is defined in the 2-D gridded data. The get() method will return the interpolated data by trilinear algorithm.

>>> xlist = np.array([0, 1, 2, 3, 4])
>>> ylist = np.array([0, 1, 2, 3])
>>> vallist = np.arange(20).reshape([5, 4])
>>> gsf = GridScalarField2D(xlist, ylist, vallist)
>>> print(gsf.get(1, 2))
6.0

bin_index() returns the indexes.

>>> gsf.bin_index(1.5, 3)
(1, 2)

Note that the y index iy=2, not 3. This ensures the ylist[iy+1] can always refered to.

Parameters:

• xlist – List of x coordinates. All elements should be increasing, i.e. xlist[i] < xlist[i+1] for all i.

• ylist – List of y coordinates.

• vallist – List of value for the field. np.array of (Nx, Ny) shape.

set_intepolation_method(method)

get_field(x, y)

Override this method for implementation.

interpolate_trilinear(x, y)

Return the interpolated value. No domain check is done.

xdomain()

Return (x0, x1), xdomain.

ydomain()

Return (y0, y1), ydomain.

bin_index(x, y)

Return the (left) index of the given x, y, z

For each component, i value satisfying xlist[i] <= x < xlist[i+1] is returned.

If x is the upper most value, the index i = len(xlist)-2 is returned. This is to ensure that the returned index is “left” bin so that i+1 can be always used to refer to for further calculation.

pcolor(ax=None, *args, **kwds)

It draws a pseud-color map. See pcolor() function for details.

irfpy.util.fields.slice(scalar_field3d, neworigin, newxnorm, newynorm, newxgrid, newygrid)

Slice the 3D field, return the 2D field.

Parameters:

• scalar_field3d – A 3-D scalar field.

• neworigin – The 3-D coorinates corresponding to the origin for new 2D field.

• newxnorm – (3,)-shape array (i.e. vector) of the x-axis of new 2D field.

• newynorm – (3,)-shape array (i.e. vector) of the y-axis of new 2D field. It is not necessarily perpendicular to the newxnorm.

• newxgrid – (Nx,)-shape array of coordinates of the x-axis of new 2D field.

• newygrid – (Ny,)-shape array of coordinates of the y-axis of new 2D field.

Returns:

A sliced 2D field.

Return type:

GridScalarField2D

A 3D scalar field is sliced to make 2D scalar field. It is then easy to visualize.

A simple example follows.

First, a gravity potential field is created

>>> gp = GravityPotential(1/6.67408e-11)     # Gravity potential at the origin with a mass of 1/G

Then, slice it in the x-y plane.

• The origin is the same, [0, 0, 0].

• The new x-axis is the same as the original x-axis, i.e., give [1, 0, 0].

• The new y-axis is the same as the original y-axis, i.e., give [0, 1, 0].

• The new x-grid is defined by the domain between -10 and 10 with 30 steps, i.e., give np.linspace(-10, 10, 30).

• The new y-grid is defined by the domain between -5 and 5 with 30 steps, i.e., give np.linspace(-5, 5, 30).

Overall, the command of slice is

>>> gp_z = slice(gp, [0, 0, 0], [1, 0, 0], [0, 1, 0], np.linspace(-10, 10, 30), np.linspace(-5, 5, 30))

irfpy.util.fields.pcolor(grid_scalar_field2d, ax=None, *args, **kwds)

Plot the 2D field

Plot the given 2d field by using pcolor method of matplotlib.

Parameters:

• grid_scalar_field2d – Scalar field.

• ax – Axis object of matplotlib . If None, a new axis is produced.

• args – Not used

• kwds – Keywords that is fed into matplotlib.pcolormesh method.

Returns:

The same as matplotlib.pyplot.pcolormesh.

Note

It is almost always good to use mpl.pcolormesh instead of mpl.pcolor for large array. For this method, mpl.pcolormesh is called internally.

Give keywords that allows the pcolormesh method;

• Give cmap keyword to change the color map

• User vmin and vmax keywords for changing the linear color scale

• Use norm=matplotlib.colors.Lognorm(vmin=0.1, vmax=100) to use log scale color

An example workflow follows.

>>> xlist = np.linspace(-5, 5, 11)
>>> ylist = np.linspace(-3, 3, 7)
>>> _xfull, _yfull = np.meshgrid(xlist, ylist, indexing='ij')
>>> vallist = 1 / (_xfull ** 2 + _yfull ** 2)    # A field of 1 / r^2
>>> field2d = GridScalarField2D(xlist, ylist, vallist)

Here comes how to make a plot.

>>> pc = irfpy.util.fields.pcolor(field2d)     
>>> import matplotlib.pyplot as plt
>>> plt.gca().colorbar(pc)     

Other examples

To specify the color map and upper/lower limits

>>> field2d.pcolor(cmap='plasma', vmin=-3, vmax=3)    

To plot in log scale in value

>>> field2d.pcolor(norm=matplotlib.colors.LogNorm(vmin=0.01, vmax=1))