calculations

This is a module that includes most, if not all, calculations within the analysator. Please type “.” and press tab to see which functions this module includes.

Example: import analysator as pt

pt.calculations. #press tab

#Output: list of functions

pt.calculations.fourier? #Press Enter

class analysator.calculations.AMRInterpolator(reader, method='linear', cellids=array([1, 2, 3, 4, 5]))

Wrapper class for interpolators, esp. at refinement interfaces. Supported methods:

linear

(nearly) C0 continuous, regular-grid trilinear interpolant extended to collapsed hexahedral cells.
Non-parametric
Exact handling of multiply-degenerate hexahedra is missing, with the enabling hack causing errors in trilinear coordinate on the order of 1m

Radial Basis Functions, rbf

Accurate, slow-ish, but hard to make properly continuous and number of neighbors on which to base the interpolant is not trivial to find.
Not continuous with regular-grid trilinear interpolants, needs to be used in the entire interpolation domain.
kword options: “neighbors” for number of neighbors (64)
basis function uses SciPy default. A free parameter.

delaunay (not recommended)

Choice of triangulation is not unique with regular grids, including refinement interfaces.
kword options as in qhull; “qhull_options” : “QJ” (breaks degeneracies)

analysator.calculations.LMN_null_lines_FOTE(LMNs, jacobs, Bs, dxs, coords)

Function that uses a linear approximation of B from B and its Jacobian in the LMN coordinates to get an approximate distance to the neutral line. inputs: * LMN basis vectors stack * Jacobian of B in 9-element vector stack * vg_b_vol

return a measure for closeness to the cell center for a neutral line and other data. The return value is the minimum of the signed distance function between the neutral line

class analysator.calculations.VariableInfo(data_array, name='', units='', latex='', latexunits='')

A class/struct for holding variable info. This includes the variable data in array form, the name and the units. Variable info is in: data, name, units.

LIST VARIABLES IN VARIBLE INFO:
data              Data of the variable (array list)
name              Name of the variable
units             Units of the variable
latex             Name of the variable in LaTeX
latexunits        Units of the variable in LaTeX

get_scaled_units(vscale=None, env='EarthSpace', manualDict=None)

Return scaling metadata

Parameters:

env – A string to choose the scaling dictionary [default: EarthSpace]
manualDict – a dictionary of {units : {scalingparams}}; used to update the included dictionary
vscale – float, factor to scale the variable with

Returns:

(norming factor, scaledUnits, scaledLatexUnits)

get_scaled_var(vscale=None, data=None, env='EarthSpace', manualDict=None)

Automatically scales the variableinfo data and adjusts the units correspondingly with the default dictionaries.

Parameters:

data – in case you wish to provide new data array (why, though?)
env – A string to choose the scaling dictionary [default: EarthSpace]
manualDict – a dictionary of {units : {scalingparams}}; used to update the included dictionary

Returns:

self, with scaled units with pre-formatted units included in the varinfo.

get_variable(index)

Creates a new variable with identical data but only certain index is included

Parameters:: index – Vector index
Returns:: an edited version of the variable

Note

E.g. if the data is an array of 3d vectors, get_variable(0) would return the variable with data[:,0] as the data

analysator.calculations.cell_time_evolution(vlsvReader_list, variables, cellids, units='')

Returns variable data from a time evolution of some certain cell ids

Parameters:

vlsvReader_list (vlsvfile.VlsvReader) – List containing VlsvReaders with a file open
variables – Name of the variables
cellids – List of cell ids
units – List of units for the variables (OPTIONAL)

Returns:

an array containing the data for the time evolution for every cell id

import analysator as pt; import matplotlib.pyplot as pl
# Example of usage:
time_data = pt.calculations.cell_time_evolution( vlsvReader_list=[VlsvReader("bulk.000.vlsv"), VlsvReader("bulk.001.vlsv"), VlsvReader("bulk.002.vlsv")], variables=["rho", "Pressure", "B"], cellids=[2,4], units=["N", "Pascal", "T"] )

# Check output
logging.info time_data

# Now plot the results:
time = time_data[0]
rho = time_data[3]
pt.plot.plot_variables(time, rho)
pl.show()

# Do post processing:
rho_data = rho.data
non_existing_example_function(rho_data)

analysator.calculations.cut3d(vlsvReader, xmin, xmax, ymin, ymax, zmin, zmax, variable, operator='pass', trim_array=False)

Retrieves variables for the given 3d cut

Parameters:

vlsvReader (vlsvfile.VlsvReader) – Some VlsvReader with a file open
xmin – The minimum x coordinate of the 2d cut
xmax – The maximum x coordinate of the 2d cut
ymin – The minimum y coordinate of the 2d cut
ymax – The maximum y coordinate of the 2d cut
zmin – The minimum z coordinate of the 2d cut
zmax – The maximum z coordinate of the 2d cut
variable – Some variable to read from the vlsv file
operator – The variable operator
trim_array – If true, shapes the array into an array with minimum amount of dimensions, e.g. if the cut is zmax-zmin=0 then this will return a 2d array

Example:
import analysator as pt
f = pt.vlsvfile.VlsvReader('example.vlsv')
three_cut = pt.calculations.cut3d( vlsvReader=f, xmin=1e6, xmax=4e6, ymin=1e6, xmax=4e6, zmin=0, zmax=0, variable="rho" )
import numpy as np
# Now three_cut is a three-dimensional array (x,y,z), but we can transform it into a 2-d array (x,y) with:
dimensions = np.shape( three_cut )
two_cut = np.reshape( three_cut, dimensions[0:2] )

analysator.calculations.cut_through(vlsvReader, point1, point2)

Returns cell ids and distances from point 1 for every cell in a line between given point1 and point2

Parameters:

vlsvReader (vlsvfile.VlsvReader) – Some open VlsvReader
point1 – The starting point of a cut-through line
point2 – The ending point of a cut-through line

Returns:

an array containing cell ids, coordinates and distances in the following format: [cell ids, distances, coordinates]. NB. Last cellid is a duplicate.

Example:
vlsvReader = VlsvReader("testfile.vlsv")
cut_through = cut_through(vlsvReader, [0,0,0], [2,5e6,0])
cellids = cut_through[0]
distances = cut_through[1]
logging.info "Cell ids: " + str(cellids)
logging.info "Distance from point 1 for every cell: " + str(distances)

analysator.calculations.cut_through_curve(vlsvReader, curve)

Returns cell ids and distances from point 1 for every cell in a line between given point1 and point2

Parameters:

vlsvReader (vlsvfile.VlsvReader) – Some open VlsvReader
point1 – The starting point of a cut-through line

Returns:

an array containing cell ids, edge distances, and the coordinates of edges in the following format: [cell ids, distances, coordinates]. NB. Last cellid is a duplicate.

Example:
vlsvReader = VlsvReader("testfile.vlsv")
cut_through_curve = cut_through(vlsvReader, [[0,0,0], [2,5e6,0]])
cellids = cut_through[0]
distances = cut_through[1]
logging.info "Cell ids: " + str(cellids)
logging.info "Distance from point 1 for every cell: " + str(distances)

analysator.calculations.cut_through_step(vlsvReader, point1, point2)

Returns cell ids and distances from point 1 to point 2, returning not every cell in a line but rather the amount of cells which corresponds with the largest axis-aligned component of the line.

Parameters:

vlsvReader (vlsvfile.VlsvReader) – Some open VlsvReader
point1 – The starting point of a cut-through line
point2 – The ending point of a cut-through line

Returns:

an array containing cell ids, coordinates and distances in the following format: [cell ids, distances, coordinates]

Example:
vlsvReader = VlsvReader("testfile.vlsv")
cut_through = cut_through_step(vlsvReader, [0,0,0], [2,5e6,0])
cellids = cut_through[0]
distances = cut_through[1]
logging.info "Cell ids: " + str(cellids)
logging.info "Distance from point 1 for every cell: " + str(distances)

analysator.calculations.dynamic_field_tracer(vlsvReader_list, x0, max_iterations, dx)

Field tracer in a dynamic time frame

Parameters:

vlsvReader_list – List of vlsv readers
x0 – The starting point for the streamlines

analysator.calculations.epsilon_M(f, cell, pop='proton', m=1.67262192595e-27, bulk=None, B=None, model='bimaxwellian', normorder=1, norm=2, threshold=0, dummy=None)

Calculates the ‘non-maxwellianity’ parameter for a distribution function f_i and its corresponding Maxwellian g_M.

Parameters:

f – VlsvReader containing VDF data
cell – CellID for the queried cell

Kword pop:

Population to calculate the parameter for

Kword m:

Species mass (default: m_p)

Kword bulk:

Bulk file name to use for moments (if available and needed)

Kword B:

Optional, user-given magnetic field vector for a bimaxwellian model distribution

Kword model:

VDF model to be used. Available models “maxwellian”, “bimaxwellian” (default)

Kword normorder:

Norm used for model-data distance measure (default: 1)

Kword norm:

Constant norm (default 2, see below)

Kword threshold:

Disregard vspace cells under this threshold [0]

Kword dummy:

If not None, generate dummy data for e.g. integration.

Returns:

scalar, non-Maxwellianity parameter for given model and norm

The definition is given by Graham et al. (2021): (1/2n) * integral(|f_i - g_M|) d^3v

valued between 0 (bi-Maxwellian) and 1 (complete deviation from a bi-Maxwellian).

NOTE that this is different to the definition by Greco et al. (2012): (1/n) * sqrt(integral((f_i - g_M)^2) d^3v)

examples comparing these two definitions can be found in Settino et al. (2021).

analysator.calculations.fourier(t, y, kaiserwindowparameter=0)

Function for returning fourier series and frequencies of some given arrays t and y

Parameters:

t – Time
y – Some variable data

Returns:

the frequencies, new time variables and frequencies

Note

return format: [FFT, frequencies, t, y]

Note

t must have a constant time stepping

Example usage:
fourier_data = fourier( t=np.arange(0,500,0.5), y=rho_data, kaiserwindowparameter=14 )

analysator.calculations.gyrophase_angles_from_file(vlsvReader, cellid)

Calculates the gyrophase angle angle distribution for a given cell with a given file :param vlsvReader: Some VlsvReader class with a file open :type vlsvReader: vlsvfile.VlsvReader :param cellid: The cell id whose gyrophase angle the user wants NOTE: The cell id must have a velocity distribution! :returns: gyrophase angles and avgs [gyro_angles, avgs]

# Example usage: vlsvReader = VlsvReader(“fullf.0001.vlsv”) result = gyrophase_angles_from_file( vlsvReader=vlsvReader, cellid=1924) # Plot the data import matplotlib.pyplot as plt plt.hist(result[0].data, weights=result[1].data, bins=100, log=False)

analysator.calculations.lineout(vlsvReader, point1, point2, variable, operator='pass', interpolation_order=1, points=100)

Returns a line cut-through from a given VLSV file for distance, coordinates and variable values. The main difference between this and cut_through is that this function interpolates a given variable.

Parameters:

vlsvReader (vlsvfile.VlsvReader) – Some open VlsvReader
point1 – The starting point of a cut-through line
point2 – The ending point of a cut-through line
variable – Variable to return
operator – The operator for the variable, for example “x” for x-component or “magnitude” for magnitude
interpolation_order – Order of interpolation (0 or 1), defaults to 1
points – Number of points to return

Returns:

A tuple with output: (distance, coordinates, variable_values)

# Example:
import analysator as pt # import analysator

vlsvReader = pt.vlsvfile.VlsvReader("testfile.vlsv") # Open a vlsv file
lineout_rho = pt.calculations.lineout( vlsvReader=vlsvReader, point1=[1.0e5, 1.0e6, 0], point2=[2.0e5, 2.0e6, 0], variable="rho", interpolation_order=1, points=100 )
distance = lineout_rho[0]
coordinates = lineout_rho[1]
values = lineout_rho[2]

analysator.calculations.pitch_angles(vlsvReader, cellid, nbins=10, filename=None, filedir=None, step=None, outputdir=None, outputfile=None, cosine=False, plasmaframe=False, vcut=None, vcutmax=None, pop='proton')

Calculates the pitch angle distribution for a given cell

Parameters:

vlsvReader – Some VlsvReader class with a file open. Can be overriden by keywords.
cellid – The cell id whose pitch angle the user wants NOTE: The cell id must have a velocity distribution!

Kword filename:

path to .vlsv file to use for input.

Kword filedir:

Optionally provide directory where files are located and use step for bulk file name

Kword step:

output step index, used for constructing output (and possibly input) filename

Kword nbins:

How many bins to use for the distribution

Kword cosine:

True if returning the pitch angles as a cosine(alpha) plot [-1,1]. If false, returns as pitch angle in degrees [0,180].

Kword plasmaframe:

True if the user wants to get the pitch angle distribution in the plasma frame (for this population). If set to a string, will try to use the string as a variable for the frame to transform into. If set to a 3-element vector, will use that frame instead.

Kword vcut:

Set to True to ignore velocity cells below 2x the thermal speed. If set to a number, will use that velocity in m/s instead.

Kword vcutmax:

Set to True to ignore velocity cells above 2x the thermal speed. If set to a number, will use that velocity in m/s instead.

Kword outputdir:

Optional (recommended) method to save results to a file in the given directory. If directory does not exist, it will be created. Filenames within directory are generated automatically.

Kword outputfile:

Provide exact output file name (including complete path)

Kword pop:

Active population, defaults to proton (avgs)

Returns:

pitch angles and avgs [pitch_angles, avgs]

# Example usage:
vlsvReader = VlsvReader("restart.0000798.vlsv")
result = pitch_angles( vlsvReader=vlsvReader, 1924, cosine=True,
                     plasmaframe=True, outputdir="/wrk/username/pitchangledirectory/" )

analysator.calculations.point_time_evolution(vlsvReader_list, variables, coordinates, units='', method='nearest')

Returns variable data from a time evolution of some certain cell ids

Parameters:

vlsvReader_list (vlsvfile.VlsvReader) – List containing VlsvReaders with a file open
variables – Name of the variables
coordinates – List of coordinates [n,3]
units – List of units for the variables (OPTIONAL)
method – name of interpolation method [‘nearest’,’linear’]

Returns:

an array containing the data for the time evolution for every coordinate

import pytools as pt; import pylab as pl
# Example of usage:
time_data = pt.calculations.point_time_evolution( vlsvReader_list=[VlsvReader("bulk.000.vlsv"), VlsvReader("bulk.001.vlsv"), VlsvReader("bulk.002.vlsv")], variables=["rho", "Pressure", "B"], coordinates=[[1e8,0,0],[1.2e8,0,0]], units=["N", "Pascal", "T"] )

# Check output
logging.info time_data

# Now plot the results:
time = time_data[0]
rho = time_data[3]
pt.plot.plot_variables(time, rho)
pl.show()

# Do post processing:
rho_data = rho.data
non_existing_example_function(rho_data)

analysator.calculations.static_field_tracer(vlsvReader, x0, max_iterations, dx, direction='+', bvar='B', centering='default', boundary_inner=-1)

Field tracer in a static frame

Parameters:

vlsvReader – An open vlsv file
x – Starting point for the field trace
max_iterations – The maximum amount of iteractions before the algorithm stops
dx – One iteration step length
direction – ‘+’ or ‘-’ or ‘+-’ Follow field in the plus direction or minus direction
bvar – String, variable name to trace [default ‘B’]
centering – String, variable centering: ‘face’, ‘volume’, ‘node’ [defaults to ‘face’]
boundary_inner – Float, stop propagation if closer to origin than this value [default -1]

Returns:

List of coordinates

analysator.calculations.static_field_tracer_3d(vlsvReader, seed_coords, max_iterations, dx, direction='+', grid_var='vg_b_vol', stop_condition=<function default_stopping_condition>, centering=None)

static_field_tracer_3d() integrates along the (static) field-grid vector field to calculate a final position. Code uses forward Euler method to conduct the tracing. Based on Analysator’s static_field_tracer() :Inputs:

param vlsvReader: A vlsvReader object (~an open .vlsv file)

param coord_list: a list of 3-element array-like initial coordinates [ [x1,y1,z1], [x2,y2,z2], … ] ::
if considering just a single starting point, the code accepts a 3-element array-like object [x1,y1,z1]

param max_iterations: The maximum number of iterations (int) before the algorithm stops. Total traced length is dx*max_iterations

param dx: One iteration step length [meters] (ex. dx=1e4 for typical applications)

keyword direction: ‘+’ or ‘-’ or ‘+-’ Follow field in the plus direction, minus direction, or both

keyword grid_var: Variable to be traced (A string)

options include:

grid_var = some string
grid_var = ‘vg_b_vol’: AMR-grid B-field grid_var =’fg_b’: fieldsolver grid B-field, grid_var=’fg_e’: fieldsolver grid E-field static_field_tracer_3d() will load the appropriate variable via the vlsvReader object NOTE: volumetric variables, with ‘_vol’ suffix, may not work as intended. Use face-centered values: ‘fg_b’, ‘fg_e’ etc.

grid_var = some fieldsolver-grid (“fg”) array. dimensions [dimx,dimy,dimz,3]
ex. fg = vlsvobj.read_variable(‘fg_b’) field grid data is already loaded externally using read_variable() method (see vlsvreader.py). If fg keyword is set this way, the input vlsvReader is only referred to for metadata (esp. grid dimensions)

keyword stop_condition: Boolean array (seed_coords.shape[0],)

Determine when the iteration stop, for the vg trace only If not specified, it will always be True for each seed points. eg. def my_stop(points):

distances = np.linalg.norm(points[:,:],axis = 1) return (distances <= lower_bound) | (distances >= upper_bound)

Returns:: points_traced — a 3d numpy array [len(seed_coords),max_iterations (or 2*max_iterations-1 for ‘+-‘),3] Non-traced sections (e.g. iterations after reaching outer boundaries) filled with np.nan

keyword centering: Set to a string (‘face’ or ‘edge’) indicating whether the fg variable is face- or edge-centered: If keyword fg == ‘fg_b’, then centering = ‘face’ (overriding input) If keyword fg == ‘fg_e’, then centering = ‘edge’ (overriding input)
EXAMPLE: vlsvobj = analysator.vlsvfile.VlsvReader(vlsvfile): fg_b = vlsvobj.read_variable(‘fg_b’) traces = static_field_tracer_3d( vlsvobj, [[5e7,0,0], [0,0,5e7]], 10, 1e5, direction=’+’, fg = fg_b )

analysator.calculations.themis_observation_from_file(vlsvReader, cellid, matrix=array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), countrates=True, interpolate=True, binOffset=[0.0, 0.0], pop='proton')

Calculates artificial THEMIS EMS observation from the given cell :param vlsvReader: Some VlsvReader class with a file open :type vlsvReader: vlsvfile.VlsvReader :param cellid: The cell id where the distribution is supposet to be sampled NOTE: The cell id must have a velocity distribution! :param matrix: Matrix to transform velocities from simulation space into detector space (use simulation_to_spacecraft_frame helper function) :param countrates: Transform phase space densities into count rates? :param interpolate: interpolate into detector bins? :param binOffset: offset bin allocation in (angle, energy) by this amount (used to align polar plots) :returns: detector bins angles, energies and counts [bin_thetas,energies,min,max,counts]

# Example usage: vlsvReader = VlsvReader(“fullf.0001.vlsv”) angles, energies, vmin, vmax, values = themis_observation_from_file( vlsvReader=self.vlsvReader, cellid=cellid) # plot: grid_r, grid_theta = np.meshgrid(energies,angles) fig,ax=pl.subplots(subplot_kw=dict(projection=”polar”),figsize=(12,10)) ax.set_title(“Detector view at cell “ + str(cellid)) cax = ax.pcolormesh(grid_theta,grid_r,values, norm=matplotlib.colors.LogNorm(vmin=vmin,vmax=vmax)) pl.show()

analysator.calculations.themis_plot_detector(vlsvReader, cellID, detector_axis=array([0, 1, 0]), pop='proton'): Plots a view of the detector countrates using matplotlib :param vlsvReader: Some VlsvReader class with a file open :type vlsvReader: vlsvfile.VlsvReader :param cellid: The cell id where the distribution is supposet to be sampled NOTE: The cell id must have a velocity distribution! :param detector_axis: detector axis direction (note: this is not spacecraft spin axis!)

analysator.calculations.themis_plot_phasespace_contour(vlsvReader, cellID, plane_x=array([1., 0., 0.]), plane_y=array([0., 0., 1.]), smooth=False, xlabel='Vx', ylabel='Vy', pop='proton'): Plots a contour view of phasespace, as seen by a themis detector, at the given cellID :param vlsvReader: Some VlsvReader class with a file open :type vlsvReader: vlsvfile.VlsvReader :param cellid: The cell id where the distribution is supposet to be sampled NOTE: The cell id must have a velocity distribution! :param plane_x and plane_y: x and y direction of the resulting plot plane

analysator.calculations.themis_plot_phasespace_helistyle(vlsvReader, cellID, plane_x=array([1., 0., 0.]), plane_y=array([0., 0., 1.]), smooth=True, xlabel='Vx', ylabel='Vy'): Plots a view of phasespace, as seen by a themis detector, at the given cellID, in the style that heli likes. :param vlsvReader: Some VlsvReader class with a file open :type vlsvReader: vlsvfile.VlsvReader :param cellid: The cell id where the distribution is supposet to be sampled NOTE: The cell id must have a velocity distribution! :param smooth: Smooth re-gridded phasespace before plotting :param plane_x and plane_y: x and y direction of the resulting plot plane

analysator.calculations.vlsv_intpol_file(file_vlsv, file_orbit, varlist, file_output)

Writes interpolated values of variables in ascii file :param file_vlsv: VLSV file :param file_orbit: Orbit file (columns: x,y,z or t,x,y,z) :param varlist: Variable list :param file_output: Output ascii file (columns: x,y,z,cellid,var1,var2,var3,…) :returns: none

# Example:
import analysator as pt
pt.calculations.vlsv_intpol_file("state00040000.vlsv","orbit.dat",["cellB","n_H+sw_ave"],"output.dat")

analysator.calculations.vlsv_intpol_points(vlsvReader, points, varlist, operator='pass', interpolation_order=1)

Returns interpolated values of variables at given points :param vlsvReader: Some open VlsvReader :type vlsvReader: vlsvfile.VlsvReader :param points: Coordinate points :param varlist: Variable list, if empty all variables are got :param operator: The operator for the variable, for example “x” for x-component or “magnitude” for magnitude :param interpolation_order: Order of interpolation (0 or 1), defaults to 1 :returns: A tuple with output: (coordinates,variable_values,header_string)

# Example:
import analysator as pt
import numpy as np
f=pt.vlsvfile.VlsvReader(file_name="state00040000.vlsv")
mesh_limits = f.get_spatial_mesh_extent()
x = np.linspace(mesh_limits[0],mesh_limits[3],100) # points along x-axis
y = np.zeros(len(x))
z = np.zeros(len(x))
points=(np.array([x,y,z])).transpose()
varlist = ["cellB","n_H+sw_ave"]
[crd,cellids,params,hstr]=pt.calculations.vlsv_intpol_points(f,points,varlist)
x=crd[:,0]
y=crd[:,1]
z=crd[:,2]
Bx=params[:,0]
By=params[:,1]
Bz=params[:,2]
n=params[:,3]