viscid.calculator package

Module contents

Calculate using fields

This package has some general purpose calculator stuff, like interpolation, streamline tracing, Div, Curl, calculate flux function, etc.

viscid.calculator.convective_deriv(a, b=None, bnd=True)[source]

Compute (a dot nabla) b for vector fields a and b

viscid.calculator.project_vector(a, b)[source]

Calculates the vector (a dot b_hat) * b_hat

viscid.calculator.project_along_line(line, fld, interp_kind='trilin')[source]

Project a Vector Field Parallel to a streamline

Parameters:
  • line (ndarray) – 3xN of points along the
  • fld (VectorField) – Field to interpolate and project onto the line
  • interp_kind (str) – which interpolation to use, const or trilin
viscid.calculator.resample_lines(lines, factor=1, kind='linear')[source]

Resample a bunch of lines

Parameters:
  • lines (sequence) – a list of 3xN ndarrays that are lines
  • factor (float) – how to resample. 2 will give a point at the midpoint of all points, 3 will give two points, etc. 0.5 will leave half the number of points, etc. This should be > 0.0.
  • kind (str) – How to interpolate; must be a value recognized by scipy.interpolate.interp1d()
Raises:

ValueError – List of resampled 3xN ndarrays
viscid.calculator.integrate_along_line(line, fld, reduction='dot', mask_func=None, interp_kind='trilin')[source]

Integrate the value of fld along a line

Parameters:
  • line (ndarray) – 3xN ndarray
  • fld (Field) – Field to interpolate / integrate
  • reduction (str) – If fld is a vector field, what quantity to integrate. Can be “dot” to dot the vectors with ds along the line, or “mag” to integrate the magnitude.
  • interp_kind (str) – which interpolation to use, const or trilin
Returns:

a scalar with the same dtype as fld
viscid.calculator.integrate_along_lines(lines, fld, reduction='dot', mask_func=None, interp_kind='trilin')[source]

Integrate the value of fld along a list of lines

Parameters:
  • lines (list) – list of 3xN ndarrays, N needs not be the same for all lines
  • fld (Field) – Field to interpolate / integrate
  • reduction (str) – If fld is a vector field, what quantity to integrate. Can be “dot” to dot the vectors with ds along the line, or “mag” to integrate the magnitude.
  • interp_kind (str) – which interpolation to use, const or trilin
Returns:

ndarray with shape (len(lines), )
viscid.calculator.jacobian_at_point(B, x, y, z, dx=None, dy=None, dz=None)[source]

Get the Jacobian at a point

If dx|dy|dz == None, then their set to the grid spacing at the point of interest

Returns:The jacobian as a 3x3 ndarray. The result is in xyz order, in other words:
[ [d_x Bx, d_y Bx, d_z Bx],
  [d_x By, d_y By, d_z By],
  [d_x Bz, d_y Bz, d_z Bz] ]
viscid.calculator.jacobian_at_ind(B, ix, iy, iz)[source]

Get the Jacobian at index

Returns:The jacobian as a 3x3 ndarray. The result is in xyz order, in other words:
[ [d_x Bx, d_y Bx, d_z Bx],
  [d_x By, d_y By, d_z By],
  [d_x Bz, d_y Bz, d_z Bz] ]
viscid.calculator.jacobian_eig_at_point(B, x, y, z, dx=None, dy=None, dz=None)[source]

Get the eigen vals/vecs of the jacobian

Returns: evals, evecs (3x3 ndarray)
The evec[:, i] corresponds to evals[i]. Eigen vectors are returned in xyz order, aka evec[:, 0] is [x, y, z] for the 0th eigen vector
viscid.calculator.jacobian_eig_at_ind(B, ix, iy, iz)[source]

Get the eigen vals/vecs of the jacobian

Returns: evals, evecs (3x3 ndarray)
The evec[:, i] corresponds to evals[i]. Eigen vectors are returned in xyz order, aka evec[:, 0] is [x, y, z] for the 0th eigen vector
viscid.calculator.div_at_point(A, x, y, z, dx=None, dy=None, dz=None)[source]

Returns divergence at a point

viscid.calculator.curl_at_point(A, x, y, z, dx=None, dy=None, dz=None)[source]

Returns curl at point as ndarray with shape (3,) xyz

viscid.calculator.extend_boundaries(fld, nl=1, nh=1, axes='all', nr_comp=None, order=1, crd_order=1, crds=None, invarient_dx=0.0, crd_invarient_dx=1.0)[source]

Extend and pad boundaries of field (leaves new corners @ 0.0)

Parameters:
  • fld (Field) – Field to extend
  • nl (int) – extend this many cells in lower direction
  • nh (int) – extend this many cells in upper direction
  • axes (str, list) – something like ‘xyz’, [‘x’, ‘theta’], etc.
  • nr_comp (int, None) – index of shape that corresponds to vector component dimension
  • order (int) – extrapolation order for data; 0 for repeating boundary values or 1 for linear extrapolation
  • crd_order (int) – extrapolation order for crds; 0 for repeating boundary values or 1 for linear extrapolation
  • crds (Coordinates) – Use these coordinates, no extrapolate them
  • invarient_dx (float) – if fld has a single cell in a given dimension, and order == 1, the field is extended with constant dx, using 0 here is synonymous with 0th order
  • invarient_dx – if fld has a single cell in a given dimension, and crd_order == 1, the coordinates in that dimension are extended with constant dx, using 0 here is synonymous with crd_order = 0.
Returns:

A new extended / padded ndarray
Return type:

ndarray

Warning

When using crd_order=0 (0 order hold), coordinates will be extended by repeating boundary values, so min_dx will be 0. Keep this in mind if dividing by dx.

viscid.calculator.extend_boundaries_ndarr(arr, nl=1, nh=1, axes='all', nr_comp=None, order=1, invarient_dx=0.0)[source]

Extend and pad boundaries of ndarray (leaves new corners @ 0.0)

Parameters:
  • arr (ndarray) – Array to extend
  • nl (int) – extend this many cells in lower direction
  • nh (int) – extend this many cells in upper direction
  • axes (list) – list of ints that correspond to axes that should be extended. If nr_comp != None, then axes > nr_comp will be shifted by 1.
  • nr_comp (int, None) – index of shape that corresponds to vector component dimension
  • order (int) – 0 for repeating boundary values or 1 for linear extrapolation
  • invarient_dx (float) – if len(arr) == 1 and order == 1, then the array is extended with constant dx, using 0 here is synonymous with 0th order.
Returns:

A new extended / padded ndarray
Return type:

ndarray
viscid.calculator.distance_to_clusters(point, clusters, alt=())[source]

L2 distance between point and clusters

viscid.calculator.find_clusters(indx, indy, x, y, multiple=True, periodic=(False, False), diagonals=True)[source]

Cluster and average groups of neighboring points

TODO: If absolutely necessary, could do some K-means clustering
here by calling into scikit-learn.
Parameters:
  • indx (sequence) – list of x indices
  • indy (sequence) – list of y indices
  • x (sequence) – x coordinate array
  • y (sequence) – y coordinate array
  • multiple (bool) – If False, average all points as a single cluster
  • periodic (sequence) – indicate whether that direction is periodic, and if so, whether the coordinate arrays are overlapped or not. Values can be True, False, or ‘+’. ‘+’ indicates that x[0] and x[-1] are not colocated, so assume they’re dx apart where dx = x[-1] - x[-2].
  • diagonals (bool) – if true, then diagonal points are considered neighbors
Returns:

2xN for N clusters
Return type:

ndarray
viscid.calculator.make_ecfc_field_leading(fc, trim_leading=True)[source]

Standardize staggering on edge / face centered fields

viscid.calculator.fc2cc(fc, force_numpy=False, bnd=True)[source]

Average a face centered field to cell centers

viscid.calculator.ec2cc(ec, force_numpy=False, bnd=True)[source]

Average an edge centered field to cell centers

viscid.calculator.div_fc(fc, force_numpy=False, bnd=True)[source]

Calculate cell centered divergence of face centered field

viscid.calculator.evaluate(grid, result_name, eqn, try_numexpr=True, slc=Ellipsis)[source]

Evaluate an equation on a grid

Examples

To use this function directly

>>> evaluator.enabled = True
>>> f = viscid.load_file("...")
>>> evaluator.evaluate(f.get_grid(),
                       "sqrt(vx**2+vy**2+vz**2)",
                       "speed")
<viscid.field.ScalarField object at ...>

Or, for short, you can as a grid to evaluate implicitly,

>> evaluator.enabled = True >> f = viscid.load_file(“…”) >> speed = f[“speed=sqrt(vx**2+vy**2+vz**2)”] <viscid.field.ScalarField object at …>
Parameters:
  • grid – a grid instance where the fields live
  • result_name (str) – Used for the name and pretty_name of the resulting field
  • eqn (str) – the equation, if a symbol exists in the numpy namespace, then that’s how it is interpreted, otherwise, the symbol will be looked up in the grid
Returns:

Field instance
viscid.calculator.minvar(B, p1, p2, n=40)[source]

Find minimum variance eigenvectors of B between p1 and p2

Do minimum variance analysis (MVA) of vector field B between two points, i.e., p1 and p2 are like two spacecraft locations on either side of the boundary for the MVA.

Parameters:
  • B (viscid.field.VectorField) – Vector field for MVA
  • p1 (sequence) – xyz point on one side of the boundary
  • p2 (sequence) – xyz point on the other side of the boundary
  • n (int) – number of points to sample B for the MVA
Returns:

(evals, evecs) All are ordered toward increasing eigenvalue magnitude. evecs is a 2d array where the vectors are columns (2nd axis), i.e., the min variance eigenvector is evecs[0, :]
viscid.calculator.minvar_series(arr, warn_finite=True)[source]

Find minimum variance eigenvectors of vector series

Parameters:arr (ndarray) – vector time series with shape (nsamples, 3)
Returns:(evals, evecs) All are ordered toward increasing eigenvalue magnitude. evecs is a 2d array where the vectors are columns (2nd axis), i.e., the min variance eigenvector is evecs[0, :]
viscid.calculator.minvar_around(B, p0, l=1.0, path_dir=(1, 0, 0), n=40)[source]

Find minimum variance eigenvectors of B around point p0

Do minimum variance analysis (MVA) of vector field B a distance l/2 around point p0 in the direction of path_dir, i.e., path_dir is like the spacecraft travel direction.

Parameters:
  • B (viscid.field.VectorField) – Vector field for MVA
  • p0 (sequence) – xyz center point, within -l/2 and +l/2 of the boundary in the path_dir directon
  • l (float) – distance on either side of center point
  • path_dir (tuple) – spacecraft travel direction
  • n (int) – number of points to sample B for the MVA
Returns:

(evals, evecs) All are ordered toward increasing eigenvalue magnitude. evecs is a 2d array where the vectors are columns (2nd axis), i.e., the min variance eigenvector is evecs[0, :]

See also

minvar()

viscid.calculator.find_minvar_lmn(B, p1, p2, n=40, l_basis=(0, 0, 1))[source]

Find rotation matrix for going to boundary normal (lmn) crds

Parameters:
  • B (viscid.field.VectorField) – Vector field for MVA
  • p1 (sequence) – xyz point on one side of the boundary
  • p2 (sequence) – xyz point on the other side of the boundary
  • n (int) – number of points to sample B for the MVA
  • l_basis (sequence) – vector used to define the l-direction, i.e., the l direction will be l_basis projected perpendicular to the boundary normal directon. If None, then l_basis will be set to evec_max.
Returns:

3x3 array of l,m,n basis column (2nd axis) vectors
Return type:

ndarray
viscid.calculator.find_minvar_lmn_around(B, p0, l=1.0, path_dir=(1, 0, 0), n=40, l_basis=(0, 0, 1))[source]

Find rotation matrix for going to boundary normal (lmn) crds

Finds boundary normal from minimum variance analysis (MVA) of vector field B a distance l/2 around point p0 in the direction of path_dir, i.e., path_dir is like the spacecraft travel direction.

Parameters:
  • B (viscid.field.VectorField) – Vector field for MVA
  • p0 (sequence) – xyz center point, within -l/2 and +l/2 of the boundary in the path_dir directon
  • l (float) – distance on either side of center point
  • path_dir (tuple) – spacecraft travel direction
  • n (int) – number of points to sample B for the MVA
  • l_basis (sequence) – vector used to define the l-direction, i.e., the l direction will be l_basis projected perpendicular to the boundary normal directon. If None, then l_basis will be set to evec_max.
Returns:

3x3 array of l,m,n basis column (2nd axis) vectors
Return type:

ndarray

See also

find_nlm()

viscid.calculator.paraboloid(y, z, x0, y0, z0, ax, ay, az)[source]

Generic paraboloid function

viscid.calculator.paraboloid_normal(y, z, x0, y0, z0, ax, ay, az, normalize=True)[source]

Normal vector of a generic paraboloid

viscid.calculator.fit_paraboloid(fld, p0=(9.0, 0.0, 0.0, 1.0, -1.0, -1.0), tolerance=0.0)[source]

Fit paraboloid it GSE coordinates x ~ y**2 + z**2

Parameters:
  • fld (viscid.field.ScalarField) – field of x values
  • p0 (sequence) – initial guess for parabaloid (x0, y0, z0, ax, ay, az), where (x0, y0, z0) is the nose location (should be subsolar for 0 dipole tilt), and the ay, az, and az coefficients determine the curvature
Returns:

record array of parameters with length 2; the

1st value is the fit value, and the 2nd is one sigma of the fit
Return type:

numpy.recarray
viscid.calculator.get_mp_info(pp, b, j, e, cache=True, cache_dir=None, slc='x=5.5j:11.0j, y=-4.0j:4.0j, z=-3.6j:3.6j', fit='mp_xloc', fit_p0=(9.0, 0.0, 0.0, 1.0, -1.0, -1.0))[source]

Get info about m-pause as flattened fields

Notes

The first thing this function does is mask locations where the GSE-y current density < 1e-4. This masks out the bow shock and current free regions. This works for southward IMF, but it is not very general.

Parameters:
  • pp (ScalarcField) – pressure
  • b (VectorField) – magnetic field
  • j (VectorField) – current density
  • e (VectorField, None) – electric field (same centering as b). If None, then the info that requires E will be filled with NaN
  • cache (bool, str) – Save to and load from cache, if “force”, then don’t load from cache if it exists, but do save a cache at the end
  • cache_dir (str) – Directory for cache, if None, same directory as that file to which the grid belongs
  • slc (str) – slice that gives a box that contains the m-pause
  • fit (str) – to which resulting field should the paraboloid be fit, defaults to mp_xloc, but pp_max_xloc might be useful in some circumstances
  • fit_p0 (tuple) – Initial guess vector for paraboloid fit
Returns:

Unless otherwise noted, the entiries are 2D (y-z) fields

  • mp_xloc location of minimum abs(Bz), this works better than max of J^2 for FTEs
  • mp_sheath_edge location where Jy > 0.1 * Jy when coming in from the sheath side
  • mp_sphere_edge location where Jy > 0.1 * Jy when coming in from the sphere side
  • mp_width difference between m-sheath edge and msphere edge
  • mp_shear magnetic shear taken 6 grid points into the m-sheath / m-sphere
  • pp_max max pp
  • pp_max_xloc location of max pp
  • epar_max max e parallel
  • epar_max_xloc location of max e parallel
  • paraboloid numpy.recarray of paraboloid fit. The parameters are given in the 0th element, and the 1st element contains the 1-sigma values for the fit
Return type:

dict
Raises:

RuntimeError – if using MHD crds instead of GSE crds
viscid.calculator.find_mp_edges(j_block, msphere_thresh=0.1, sheath_thresh=0.1, maskval=0.0001)[source]

Find x location of msphere and msheath edges using current (J)

Note

GSE coordinates only please

Parameters:
  • j_block (VectorField) – Current density containing the whole magnetopause
  • msphere_thresh (float) – thereshold of current on the magnetosphere side as a fraction of the maximum current density, i.e., 0.1 is 10% of the max
  • sheath_thresh (float) – thereshold of current on the magnetosheath side as a fraction of the maximum current density, i.e., 0.1 is 10% of the max
  • maskval (float, None) – if not None, then mask out J values less than maskval; useful for masking out bowshock, and current free regions
Returns:

sheath and sphere fields / values

  • sheath_edge: float or 2D ScalarField of x values
  • msphere_edge: float or 2D ScalarField of x values
  • mp_width: sheath_edge - msphere_edge
  • sheath_ind: index of sheath_edge x location
  • sphere_ind: index of msphere_edge x location
Return type:

tuple
viscid.calculator.calc_psi(B, rev=False)[source]

Calc Flux function (only valid in 2d)

Parameters:
  • B (VectorField) – magnetic field, should only have two spatial dimensions so we can infer the symmetry dimension
  • rev (bool) – since this integration doesn’t like going through undefined regions (like within 1 earth radius of the origin for openggcm), you can use this to start integrating from the opposite corner.
Returns:

2-D scalar flux function
Return type:

ScalarField
Raises:

ValueError – If B has <> 2 spatial dimensions
viscid.calculator.calc_beta(pp, B, scale=1.0)[source]

Calc plasma beta (2.0 * p / B^2)

Parameters:
  • pp (ScalarField or ndarray) – pressure
  • B (VectorField) – magnetic field
  • scale (float, optional) – overall scale factor
Returns:

Plasma beta
Return type:

ScalarField

Note

For OpenGGCM, where pp is in pPa and B is in nT, scale should be 40.0.

viscid.calculator.trace_separator(grid, b_slcstr='x=-25j:15j, y=-30j:30j, z=-15j:15j', r=1.0, plot=False, trace_opts=None, cache=True, cache_dir=None)[source]

Trace a separator line from most dawnward null

Still in testing Uses the bisection algorithm.

Parameters:
  • grid (Grid) – A grid that has a “b” field
  • b_slcstr (str) – Some valid slice for B field
  • r (float) – spatial step of separator line
  • plot (bool) – make debugging plots
  • trace_opts (dict) – passed to streamline function
  • cache (bool, str) – Save to and load from cache, if “force”, then don’t load from cache if it exists, but do save a cache at the end
  • cache_dir (str) – Directory for cache, if None, same directory as that file to which the grid belongs
Raises:

IOError – Description
Returns:

(separator_lines, nulls)

  • separator_lines (list): list of M 3xN ndarrays that represent M separator lines with N points
  • nulls (ndarray): 3xN array of N null points
Return type:

tuple
viscid.calculator.topology_bitor_clusters(fld, min_depth=1, max_depth=10, multiple=True, plot=False, sep_val=15, mask_limit=15, periodic='00', pt_bnds=())[source]

Find separator as intersection of all global topologies

Neighbors are bitwise ORed until at least one value matches sep_val which is presumably (Close | Open N | Open S | SW). This happens between min_depth and max_depth times, where the resolution of each iteration is reduced by a factor of two, ie, worst case 2**(max_depth).

Parameters:
  • fld (Field) – Topology (bitmask) as a field
  • min_depth (int) – Iterate at least this many times
  • max_depth (int) – Iterate at most this many times
  • multiple (bool) – passed to viscid.cluster()
  • sep_val (int) – Value of bitmask that indicates a separator
  • plot (bool) – Make a 2D plot of Fld and the sep candidates
  • mask_limit (int) – if > 0, then bitmask fld with mask_limit, i.e., fld = fld & mask_limit (bitwise and)
  • periodic (sequence) – indicate whether that direction is periodic, and if so, whether the coordinate arrays are overlapped or not. Values can be True, False, or ‘+’. ‘+’ indicates that x[0] and x[-1] are not colocated, so assume they’re dx apart where dx = x[-1] - x[-2].
  • pt_bnd (sequence) – Boundaries that come to a point, i.e., all values along that boundary are neighbors such as the poles of a sphere. Specified like “0-” for lower boundary of dimension 0 or “1+” for the upper boundary of dimension 1.
Returns:

2xN for N clusters of separator points in the same coordinates as fld
Return type:

ndarray
viscid.calculator.get_sep_pts_bitor(fld, seed, trace_opts=None, make_3d=True, **kwargs)[source]

bitor topologies to find separator points in uv map from seed

Parameters:
  • fld (VectorField) – Magnetic Field
  • seed (viscid.seed.SeedGen) – Any Seed generator with a 2d local representation
  • trace_opts (dict) – kwargs for calc_streamlines
  • make_3d (bool) – convert result from uv to 3d space
  • **kwargs – passed to topology_bitor_clusters()
Returns:

3xN ndarray of N separator points in uv space or 3d space depending on the make_3d kwarg
viscid.calculator.get_sep_pts_bisect(fld, seed, trace_opts=None, min_depth=3, max_depth=7, plot=False, perimeter_check=<function perimeter_check_bitwise_or>, make_3d=True)[source]

bisect uv map of seed to find separator points

Parameters:
  • fld (VectorField) – Magnetic Field
  • seed (viscid.seed.SeedGen) – Any Seed generator with a 2d local representation
  • trace_opts (dict) – kwargs for calc_streamlines
  • min_depth (int) – Min allowable bisection depth
  • max_depth (int) – Max bisection depth
  • plot (bool) – Useful for debugging the algorithm
  • perimeter_check (func) – Some func that returns a bool with the same signature as perimeter_check_bitwise_or()
  • make_3d (bool) – convert result from uv to 3d space
Returns:

3xN ndarray of N separator points in uv space or 3d space depending on the make_3d kwarg
viscid.calculator.topology2color(topology, topo_style='msphere', bad_color=None)[source]

Determine RGB from topology value

Parameters:
  • topology (int, list, ndarray) – some value in calculator.streamline.TOPOLOGY_*
  • topo_style (string) – msphere, or a dict with its own mapping
  • bad_color (tuple) – rgb color for invalid topologies
Returns:

Nx3 array of rgb data or (R, G, B) tuple if topology is a single value