postpic package

POSTPIC
The open source particle-in-cell post processor.

Particle-in-cell simulations are a valuable tool for the simulation of non-equelibrium systems in plasma- or astrophysics. Such simulations usually produce a large amount of data consisting of electric and magnetic field data as well as particle positions and momenta. While there are various PIC codes freely available, the task of post-processing – essentially condensing the large amounts of data into small units suitable for plotting routines – is typically left to each user individually. As post-processing may be a time consuming and error-prone process, this python package has been developed.

Postpic can handle two different types of data:

Field data
which is data sampled on a predefined grid, such as electic and magnetic fields, particle- or charge densities, currents, etc. Fields are usually the data, which can be plotted directly. See postpic.Field.
Particle data
which is data of multiple particles and for each particle positions (x, y, z) and momenta (px, py, pz) are known. Particles usually also have weight, charge, time and a unique id. Postpic can transform particle data to field data using the same algorithm and particle shapes, which are used in most PIC Simulations. The particle-to-grid routines are written in C for maximum performance. See postpic.MultiSpecies.
class postpic.Field(matrix, name='', unit='', **kwargs)[source]

Bases: postpic._compat.mixins.NDArrayOperatorsMixin

The Field Object carries data in form of an numpy.ndarray together with as many Axis objects as the data’s dimensions. Additionaly the Field object provides any information that is necessary to plot _and_ annotate the plot.

Create a Field object from scratch. The only required argument is matrix which contains the actual data.

A name and a unit may be supplied.

The axis may be specified in different ways:

  • by passing a list of Axis object as axes
  • by passing arrays with the grid_nodes as xedges, yedges and zedges. This is intended to work with np.histogram.
  • by not passing anything, which will create default axes from 0 to 1.
T

Return the Field with the order of axes reversed. In 2D this is the usual matrix transpose operation.

__array__(dtype=None)[source]

will be called by numpy function in case an numpy array is needed.

__array_priority__ = 1
__array_ufunc__(ufunc, method, *inputs, **kwargs)[source]
__array_wrap__(array, context=None)[source]
__copy__()[source]

returns a shallow copy of the object. This method is called by copy.copy(obj). Just copy enough to create copies for operator overloading.

__getitem__(key)[source]
__setitem__(key, other)[source]
adjust_stagger_to(other)[source]
all(axis=None, out=None, keepdims=False)

Returns True if all elements evaluate to True.

Refer to numpy.all for full documentation.

See also

numpy.all()
equivalent function
angle
any(axis=None, out=None, keepdims=False)

Returns True if any of the elements of a evaluate to True.

Refer to numpy.any for full documentation.

See also

numpy.any()
equivalent function
atleast_nd(n)[source]

Make sure the field has at least ‘n’ dimensions

autocutout(axes=None, fractions=(0.001, 0.002))[source]

Automatically cuts out the main feature of the field by removing border regions that only contain small numbers.

This is done axis by axis. For each axis, the mean across all other axes is taken. The maximum max of the remaining 1d-array is taken and searched for the outermost boundaries a, d such that all values out of array[a:d] are smaller then fractions[0]*max. A second set of boundaries b, c is searched such that all values out of array[b:c] are smaller then fractions[1]*max. Because fractions[1] should be larger than fractions[0], array[b:c] should be contained completely in array[a:d].

A padding length x is chosen such that array[b-x:c+x] is entirely within array[a:d].

Then the corresponding axis of the field is sliced to [b-x:c+x] and multiplied with a tukey-window such that the region [b:c] is left untouched and the field in the padding region smoothly vanishes on the outer border.

This process is repeated for all axes in axes or for all axes if axes is None.

autoreduce(maxlen=4000)[source]

Reduces the Grid to a maximum length of maxlen per dimension by just executing half_resolution as often as necessary.

clip(a_min, a_max, out=None)[source]
conj()[source]
cutout(newextent)[source]

only keeps that part of the data, that belongs to newextent.

derivative(axis, staggered=False)[source]

Calculate the derivative of the field with respect to axis.

Uses np.gradient by default which outputs the second order accurate derivative on the same grid as the input field.

If staggered=True is passed, the method will instead calculate the second order accurate derivative at the points centered between the input grid points.

dimensions

returns only present dimensions. [] and [[]] are interpreted as -1 np.array(2) is interpreted as 0 np.array([1,2,3]) is interpreted as 1 and so on…

ensure_frequency_domain()[source]
ensure_positive_axes()[source]

ensures, that all axis are going from smaller to greater numbers, i.e. none of the axis is reversed.

ensure_spatial_domain()[source]
ensure_transform_state(transform_states)[source]

Makes sure that the field has the given transform_states. transform_states may be a single boolean, indicating the same desired transform_state for all axes. It may be a list of the desired transform states for all the axes or a dictionary indicating the desired transform states of specific axes.

evaluate(ex, local_dict=None, global_dict=None, **kwargs)[source]

Evaluates the expression ex using NumExpr and returns a field containing the result. This copies all metadata from self and just replaces the matrix.

This function is basically syntactic sugar simplifying

`field.replace_data(ne.evaluate(expr))`

to

`field.evaluate(expr)`

This method replicates some logic from NumExpr.necompiler.getArguments(), seems there is no way around it.

export(filename, **kwargs)[source]
Uses postpic.export_field to export this field to a file. All ``**kwargs` will be forwarded to this function. Format is recognized by the extension of the filename.

export Field object as a file. Format depends on the extention of the filename. Currently supported are: .npz:

uses numpy.savez.
.csv:
uses numpy.savetxt.
.vtk:
vtk export to paraview
extent

returns the extents in a linearized form, as required by “matplotlib.pyplot.imshow”.

fft(axes=None, exponential_signs='spatial', **kwargs)[source]

Performs Fourier transform on any number of axes.

The argument axis is either an integer indicating the axis to be transformed or a tuple giving the axes that should be transformed. Automatically determines forward/inverse transform. Transform is only applied if all mentioned axes are in the same transform state. If an axis is transformed twice, the origin of the axis is restored.

Parameters:
  • exponential_signs

    configures the sign convention of the exponential.

    • exponential_signs == ‘spatial’: fft using exp(-ikx), ifft using exp(ikx)
    • exponential_signs == ‘temporal’: fft using exp(iwt), ifft using exp(-iwt)
  • **kwargs – keyword-arguments are passed to the underlying fft implementation.
fft_autopad(axes=None, fft_padsize=<postpic.helper.FFTW_Pad object>)[source]

Automatically pad the array to a size such that computing its FFT using FFTW will be fast.

Parameters:fft_padsize (callable) –

The default for keyword argument fft_padsize is a callable, that is used to calculate the padded size for a given size.

By default, this uses fft_padsize=helper.fftw_padsize which finds the next larger “good” grid size according to what the FFTW documentation says.

However, the FFTW documentation also says: “(…) Transforms whose sizes are powers of 2 are especially fast.”

If you don’t worry about the extra padding, you can pass fft_padsize=helper.fft_padsize_power2 and this method will pad to the next power of 2.

flip(axis)[source]

functionality of numpy.flip.

field.flip(0) returns a postpic.Field object with the specified axis flipped. np.flip(field) returns only the numpy.ndarray and the axis information is lost.

grid
grid_nodes
half_resolution(axis)[source]

Halfs the resolution along the given axis by removing every second grid_node and averaging every second data point into one.

If there is an odd number of grid points, the last point will be ignored (that means, the extent will change by the size of the last grid cell).

Returns:the modified Field.
Return type:Field
imag
classmethod importfrom(filename, **kwargs)[source]

Construct a new field object from foreign data. Currently, the following file formats are supported explicitly: *.csv *.png

All other files will be opened using Pillow.

integrate(axes=None, method=<function simps>)[source]

Calculates the definite integral along the given axes.

Parameters:method (callable) –

Choose the method to use. Available options:

  • ’constant’
  • any function with the same signature as scipy.integrate.simps (default).
islinear()[source]
label
classmethod loadfrom(filename)[source]

Load a field object previously stored using the saveto method. These are .npz files with a specific metadata layout.

map_axis_grid(axis, transform, preserve_integral=True, jacobian_func=None)[source]

Transform the Field to new coordinates along one axis.

This function transforms the coordinates of one axis according to the function transform and applies the jacobian to the data.

Please note that no interpolation is applied to the data, instead a non-linear axis grid is produced. If you want to interpolate the data to a new (linear) grid, use the method map_coordinates() instead.

In contrast to map_coordinates(), the function transform is not used to pull the new data points from the old grid, but is directly applied to the axis. This reverses the direction of the transform. Therfore, in order to preserve the integral, it is necessary to divide by the Jacobian.

Parameters:
  • axis (int) – the index or name of the axis you want to apply transform to.
  • transform (callable) – the transformation function which takes the old coordinates as an input and returns the new grid
  • preserve_integral (bool) – Divide by the jacobian of transform, in order to preserve the integral.
  • jacobian_func (callable) – If given, this is expected to return the derivative of transform. If not given, the derivative is numerically approximated.
map_coordinates(newaxes, transform=None, complex_mode='polar', preserve_integral=True, jacobian_func=None, jacobian_determinant_func=None, **kwargs)[source]

Transform the Field to new coordinates.

Parameters:
  • newaxes (list) – The new axes of the new coordinates.
  • transform (callable) –

    a callable that takes the new coordinates as input and returns the old coordinates from where to sample the Field. It is basically the inverse of the transformation that you want to perform. If transform is not given, the identity will be used. This is suitable for simple interpolation to a new extent/shape. Example for cartesian -> polar:

    >>> def T(r, theta):
    >>>    x = r * np.cos(theta)
    >>>    y = r * np.sin(theta)
    >>>    return x, y
    

    Note that this function actually computes the cartesian coordinates from the polar coordinates, but stands for transforming a field in cartesian coordinates into a field in polar coordinates.

    However, in order to preserve the definite integral of the field, it is necessary to multiply with the Jacobian determinant of T.

    \[\tilde{U}(r, \theta) = U(T(r, \theta)) \cdot \det \frac{\partial (x, y)}{\partial (r, \theta)}\]

    such that

    \[\int_V \mathop{\mathrm{d}x} \mathop{\mathrm{d}y} U(x,y) = \int_{T^{-1}(V)} \mathop{\mathrm{d}r}\mathop{\mathrm{d}\theta} \tilde{U}(r,\theta)\,.\]
  • complex_mode

    The complex_mode specifies how to proceed with complex data.

    • complex_mode = ‘cartesian’ - interpolate real/imag part (fastest)
    • complex_mode = ‘polar’ - interpolate abs/phase If skimage.restoration is available, the phase will be unwrapped first (default)
    • complex_mode = ‘polar-no-unwrap’ - interpolate abs/phase Skip unwrapping the phase, even if skimage.restoration is available
  • preserve_integral (bool) –

    If True (the default), the data will be multiplied with the Jacobian determinant of the coordinate transformation such that the integral over the data will be preserved.

    In general, you will want to do this, because the physical unit of the new Field will correspond to the new axis of the Fields. Please note that Postpic, currently, does not automatically change the unit members of the Axis and Field objects, this you will have to do manually.

    There are, however, exceptions to this rule. Most prominently, if you are converting to polar coordinates, it depends on what you are going to do with the transformed Field. If you intend to do a Cartesian r-theta plot or are interested in a lineout for a single value of theta, you do want to apply the Jacobian determinant. If you had a density in e.g. J/m^2 than, in polar coordinates, you want to have a density in J/m/rad. If you intend, on the other hand, to do a polar plot, you do not want to apply the Jacobian. In a polar plot, the data points are plotted with variable density which visually takes care of the Jacobian automatically. A polar plot of the polar data should look like a Cartesian plot of the original data with just a peculiar coordinate grid drawn over it.

  • jacobian_determinant_func (callable) – A callable that returns the jacobian determinant of the transform. If given, this takes precedence over the following option.
  • jacobian_func (callable) – a callable that returns the jacobian of the transform. If this is not given, the jacobian is numerically approximated.
  • **kwargs – Additional keyword arguments are passed to scipy.ndimage.map_coordinates, see the documentation of that function.
matrix
max(axis=None, out=None)

Return the maximum along a given axis.

Refer to numpy.amax for full documentation.

See also

numpy.amax()
equivalent function
mean(axis=None, dtype=None, out=None, keepdims=False)

Returns the average of the array elements along given axis.

Refer to numpy.mean for full documentation.

See also

numpy.mean()
equivalent function
meshgrid(sparse=True)[source]
min(axis=None, out=None, keepdims=False)

Return the minimum along a given axis.

Refer to numpy.amin for full documentation.

See also

numpy.amin()
equivalent function
ndim
pad(pad_width, mode='constant', **kwargs)[source]

Pads the data using np.pad and takes care of the axes. See documentation of numpy.pad.

In contrast to np.pad, pad_width may be given as integers, which will be interpreted as pixels, or as floats, which will be interpreted as distance along the appropriate axis.

All other parameters are passed to np.pad unchanged.

phase(do_unwrap_phase=True)[source]

Returns the (unwrapped) phase of the complex Field.

do_unwrap_phase: True, if skimage.restoration.unwrap_phase should be applied to data

prod(axis=None, dtype=None, out=None, keepdims=False)

Return the product of the array elements over the given axis

Refer to numpy.prod for full documentation.

See also

numpy.prod()
equivalent function
ptp(axis=None, out=None)

Peak to peak (maximum - minimum) value along a given axis.

Refer to numpy.ptp for full documentation.

See also

numpy.ptp()
equivalent function
real
replace_data(other)[source]
rot90(k=1, axes=(0, 1))[source]

Rotates the field by 90 degrees, k times. Rotates the field in the plane given by axes.

saveto(filename)[source]

Save a Field object as a file. Use loadfrom() to load Field objects.

setaxisobj(axis, axisobj)[source]

replaces the current axisobject for axis axis by the new axisobj axisobj.

shape
shift_grid_by(dx, interpolation='fourier')[source]

Translate the Grid by dx. This is useful to remove the grid stagger of field components.

If all axis will be shifted, dx may be a list. Otherwise dx should be a mapping from axis to translation distance.

The keyword-argument interpolation indicates the method to be used and may be one of [‘linear’, ‘fourier’]. In case of interpolation = ‘fourier’ all axes must have same transform_state.

spacing

returns the grid spacings for all axis.

squeeze()[source]

removes axes that have length 1, reducing self.dimensions.

Note, that axis with length 0 will not be removed! numpy.squeeze also does not remove length=0 directions.

Same as numpy.squeeze.

std(axis=None, dtype=None, out=None, ddof=0, keepdims=False)

Returns the standard deviation of the array elements along given axis.

Refer to numpy.std for full documentation.

See also

numpy.std()
equivalent function
sum(axis=None, dtype=None, out=None, keepdims=False)

Return the sum of the array elements over the given axis.

Refer to numpy.sum for full documentation.

See also

numpy.sum()
equivalent function
swapaxes(axis1, axis2)[source]

Swaps the axes axis1 and axis2, equivalent to numpy.swapaxes.

topolar(extent=None, shape=None, angleoffset=0, **kwargs)[source]

Transform the Field to polar coordinates.

This is a convenience wrapper for map_coordinates() which will let you easily define the desired grid in polar coordinates.

Parameters:
  • extent – should be of the form extent=(phimin, phimax, rmin, rmax) or extent=(phimin, phimax)
  • shape – should be of the form shape=(N_phi, N_r),
  • angleoffset – can be any real number and will rotate the zero-point of the angular axis.
  • complex_mode

    The complex_mode specifies how to proceed with complex data.

    • complex_mode = ‘cartesian’ - interpolate real/imag part (fastest)
    • complex_mode = ‘polar’ - interpolate abs/phase If skimage.restoration is available, the phase will be unwrapped first (default)
    • complex_mode = ‘polar-no-unwrap’ - interpolate abs/phase Skip unwrapping the phase, even if skimage.restoration is available
  • preserve_integral (bool) –

    If True (the default), the data will be multiplied with the Jacobian determinant of the coordinate transformation such that the integral over the data will be preserved.

    In general, you will want to do this, because the physical unit of the new Field will correspond to the new axis of the Fields. Please note that Postpic, currently, does not automatically change the unit members of the Axis and Field objects, this you will have to do manually.

    There are, however, exceptions to this rule. Most prominently, if you are converting to polar coordinates, it depends on what you are going to do with the transformed Field. If you intend to do a Cartesian r-theta plot or are interested in a lineout for a single value of theta, you do want to apply the Jacobian determinant. If you had a density in e.g. J/m^2 than, in polar coordinates, you want to have a density in J/m/rad. If you intend, on the other hand, to do a polar plot, you do not want to apply the Jacobian. In a polar plot, the data points are plotted with variable density which visually takes care of the Jacobian automatically. A polar plot of the polar data should look like a Cartesian plot of the original data with just a peculiar coordinate grid drawn over it.

  • jacobian_determinant_func (callable) – A callable that returns the jacobian determinant of the transform. If given, this takes precedence over the following option.
  • jacobian_func (callable) – a callable that returns the jacobian of the transform. If this is not given, the jacobian is numerically approximated.
  • **kwargs – Additional keyword arguments are passed to scipy.ndimage.map_coordinates, see the documentation of that function.
transpose(*axes)[source]

transpose method equivalent to numpy.ndarray.transpose. If axes is empty, the order of the axes will be reversed. Otherwise axes[i] == j means that the i’th axis of the returned Field will be the j’th axis of the input Field.

var(axis=None, dtype=None, out=None, ddof=0, keepdims=False)

Returns the variance of the array elements, along given axis.

Refer to numpy.var for full documentation.

See also

numpy.var()
equivalent function
class postpic.Axis(name='', unit='', **kwargs)[source]

Bases: object

Axis handling for a single Axis.

Create an Axis object from scratch.

The least required arguments are any of:
  • grid
  • grid_node
  • extent _and_ n

The remaining fields will be deduced from the givens.

More arguments may be supplied, as long as they are compatible.

__eq__(other)[source]

equality test for axis

__getitem__(key)[source]

Returns an Axis which consists of a sub-part of this object defined by a slice containing floats or integers or a float or an integer

__getstate__()[source]

Excludes self._inv_map from the pickled state

extent
grid
grid_node
half_resolution()[source]

removes every second grid_node.

islinear(force=False)[source]

Checks if the axis has a linear grid.

isreversed
label
physical_length
reversed()[source]

returns an reversed Axis object

spacing
value_to_index(value)[source]

This funtion is used to map values to indices in an interpolating manner, this is mainly used by the map_coordinates method of the Field class.

In contrast to the _find_nearest_index method, this method does not return an integer but a fractional index that refers to a position between actual pixels.

In general the equality

ax._find_nearest_index(x) == np.round(ax.value_to_index(x))

should hold.

class postpic.PhysicalConstants[source]

Bases: object

gives you some constants.

c = 299792458.0
epsilon0 = 8.854187817620389e-12
mass_u = 1.67266490646e-27
me = 9.109383e-31
mu0 = 1.2566370614359173e-06
static ncrit(laslambda)[source]

Critical plasma density in particles per m^3 for a given wavelength laslambda in m.

static ncrit_um(lambda_um)[source]

Critical plasma density in particles per m^3 for a given wavelength lambda_um in microns.

qe = 1.602176565e-19
postpic.unstagger_fields(*fields, **kwargs)[source]

Unstagger a collection of fields.

This functions shifts the origins of the grids of the given fields such that they coincide. Since the choice of the common origin is somewhat arbitrary, it might be overriden by a keyword-argument origin, as may be the interpolation method. See Field.shift_grid_by for available methods.

postpic.kspace_epoch_like(component, fields, dt, extent=None, omega_func=<function omega_free>, align_to='B')[source]

Reconstruct the physical kspace of one polarization component See documentation of kspace

This function will use special care to make sure, that the implicit linear interpolation introduced by Epochs half-steps will not impede the accuracy of the reconstructed k-space. The frequency response of the linear interpolation is modelled and removed from the interpolated fields.

dt: time-step of the simulation, this is used to calculate the frequency response due to the linear interpolated half-steps

For the current version of EPOCH, v4.9, use the following: align_to == ‘B’ for intermediate dumps, align_to == “E” for final dumps

postpic.kspace(component, fields, extent=None, interpolation=None, omega_func=<function omega_free>)[source]

Reconstruct the physical kspace of one polarization component This function basically computes one component of

E = 0.5*(E - omega/k^2 * Cross[k, B])
or
B = 0.5*(B + 1/omega * Cross[k, E]).

component must be one of [“Ex”, “Ey”, “Ez”, “Bx”, “By”, “Bz”].

The necessary fields must be given in the dict fields with keys chosen from [“Ex”, “Ey”, “Ez”, “Bx”, “By”, “Bz”]. Which are needed depends on the chosen component and the dimensionality of the fields. In 3D the following fields are necessary:

Ex, By, Bz -> Ex Ey, Bx, Bz -> Ey Ez, Bx, By -> Ez

Bx, Ey, Ez -> Bx By, Ex, Ez -> By Bz, Ex, Ey -> Bz

In 2D, components which have “k_z” in front of them (see cross-product in equations above) are not needed. In 1D, components which have “k_y” or “k_z” in front of them (see cross-product in equations above) are not needed.

The keyword-argument extent may be a list of values [xmin, xmax, ymin, ymax, …] which denote a region of the Fields on which to execute the kspace reconstruction.

The keyword-argument interpolation indicates whether interpolation should be used to remove the grid stagger. If interpolation is None, this function works only for non-staggered grids. Other choices for interpolation are “linear” and “fourier”.

The keyword-argument omega_func may be used to pass a function that will calculate the dispersion relation of the simulation may be given. The function will receive one argument that contains the k mesh.

postpic.kspace_propagate(kspace, dt, nsteps=1, **kwargs)[source]

Evolve time on a field. This function checks the transform_state of the field and transforms first from spatial domain to frequency domain if necessary. In this case the inverse transform will also be applied to the result before returning it. This works, however, only correctly with fields that are the inverse transforms of a k-space reconstruction, i.e. with complex fields.

dt: time in seconds

This function will return an infinite generator that will do arbitrary many time steps.

If yield_zeroth_step is True, then the kspace will also be yielded after removing the antipropagating waves, but before the first actual step is done.

If a vector moving_window_vect is passed to this function, which is ideally identical to the mean propagation direction of the field in forward time direction, an additional linear phase is applied in order to keep the pulse inside of the box. This effectively enables propagation in a moving window. If dt is negative, the window will actually move the opposite direction of moving_window_vect. Additionally, all modes which propagate in the opposite direction of the moving window, i.e. all modes for which dot(moving_window_vect, k)<0, will be deleted.

The motion of the window can be inhibited by specifying move_window=False. If move_window is None, the moving window is automatically enabled if moving_window_vect is given.

The deletion of the antipropagating modes can be inhibited by specifying remove_antipropagating_waves=False. If remove_antipropagating_waves is None, the deletion of the antipropagating modes is automatically enabled if moving_window_vect is given.

nsteps: number of steps to take

If nsteps == 1, this function will just return the result. If nsteps > 1, this function will return a generator that will generate the results. If you want a list, just put list(…) around the return value.

postpic.time_profile_at_plane(kspace_or_complex_field, axis='x', value=None, dir=1, t_input=0.0, **kwargs)[source]

‘Measure’ the time-profile of the propagating complex_field while passing through a plane.

The arguments axis, value and dir specify the plane and main propagation direction.

axis specifies the axis perpendicular to the measurement plane.

dir=1 specifies propagation towards positive axis, dir=-1 specifies the opposite direction of propagation.

value specifies the position of the plane along axis. If value=None, a default is chosen, depending on dir.

If dir=-1, the starting point of the axis is used, which lies at the 0-component of the inverse transform.

If dir=1, the end point of the axis + one axis spacing is used, which, via periodic boundary conditions of the fft, also lies at the 0-component of the inverse transform.

If the given value differs from these defaults, an initial propagation with moving window will be performed, such that the desired plane lies in the default position.

t_input specifies the point in time at which the input field or kspace is given. This is used to specify the time axis of the output fields.

For example axis=’x’ and value=0.0 specifies the ‘x=0.0’ plane while dir=1 specifies propagation towards positive ‘x’ values. The ‘x’ axis starts at 2e-5 and ends at 6e-5 with a grid spacing of 1e-6. The default value for the measurement plane would have been 6.1e-5 so an initial backward propagation with dt = -6.1e-5/c is performed to move the pulse in front of the’x=0.0 plane.

Additional kwargs are passed to kspace_propagate if they are not overridden by this function.

class postpic.ScalarProperty(expr, name=None, unit=None, symbol=None)[source]

Bases: object

__iter__()[source]
evaluate(vars)[source]

vars must be a dictionary containing variables used within the expression “expr”.

expr
input_names

The list of variables used within this expression.

name
symbol
unit
class postpic.MultiSpecies(dumpreader, *speciess, **kwargs)[source]

Bases: object

The MultiSpecies class. Different MultiSpecies can be added together to create a combined collection.

ignore_missing_species = False
set to true to ignore missing species.

The MultiSpecies class will return a list of values for every particle property.

Ekin()[source]

Deprecated since version unknown: The function Ekin is deprecated. Use self(“Ekin”) instead.

Ekin_MeV()[source]

Deprecated since version unknown: The function Ekin_MeV is deprecated. Use self(“Ekin_MeV”) instead.

Ekin_MeV_amu()[source]

Deprecated since version unknown: The function Ekin_MeV_amu is deprecated. Use self(“Ekin_MeV_amu”) instead.

Ekin_MeV_qm()[source]

Deprecated since version unknown: The function Ekin_MeV_qm is deprecated. Use self(“Ekin_MeV_qm”) instead.

Ekin_keV()[source]

Deprecated since version unknown: The function Ekin_keV is deprecated. Use self(“Ekin_keV”) instead.

Ekin_keV_amu()[source]

Deprecated since version unknown: The function Ekin_keV_amu is deprecated. Use self(“Ekin_keV_amu”) instead.

Ekin_keV_qm()[source]

Deprecated since version unknown: The function Ekin_keV_qm is deprecated. Use self(“Ekin_keV_qm”) instead.

Eruhe()[source]

Deprecated since version unknown: The function Eruhe is deprecated. Use self(“Eruhe”) instead.

ID()[source]

Deprecated since version unknown: The function ID is deprecated. Use self(“id”) instead.

P()[source]

Deprecated since version unknown: The function P is deprecated. Use self(“p”) instead.

Px()[source]

Deprecated since version unknown: The function Px is deprecated. Use self(“px”) instead.

Py()[source]

Deprecated since version unknown: The function Py is deprecated. Use self(“py”) instead.

Pz()[source]

Deprecated since version unknown: The function Pz is deprecated. Use self(“pz”) instead.

V()[source]

Deprecated since version unknown: The function V is deprecated. Use self(“v”) instead.

Vx()[source]

Deprecated since version unknown: The function Vx is deprecated. Use self(“vx”) instead.

Vy()[source]

Deprecated since version unknown: The function Vy is deprecated. Use self(“vy”) instead.

Vz()[source]

Deprecated since version unknown: The function Vz is deprecated. Use self(“vz”) instead.

X()[source]

Deprecated since version unknown: The function X is deprecated. Use self(“x”) instead.

X_um()[source]

Deprecated since version unknown: The function X_um is deprecated. Use self(“x_um”) instead.

Y()[source]

Deprecated since version unknown: The function Y is deprecated. Use self(“y”) instead.

Y_um()[source]

Deprecated since version unknown: The function Y_um is deprecated. Use self(“Y_mu”) instead.

Z()[source]

Deprecated since version unknown: The function Z is deprecated. Use self(“z”) instead.

Z_um()[source]

Deprecated since version unknown: The function Z_um is deprecated. Use self(“z_um”) instead.

__add__(other)[source]

Operator overloading for the + operator.

Returns:A new MultiSpecies object containing the combined collection of particles.

Note

Particles present in self and other will be present twice in the returned object!

__call__(expr)[source]

Access to particle properties via the expression, which is used to calculate them.

This is only function to actually access the data. Every other function which allows data access must call this one internally!

Supported types
  • ScalarProperty
  • str: will be converted to a ScalarProperty by particle_scalars.__call__. Therefore known quantities will be recognized
  • callable, which acts on the MultiSpecies object. This will work, but maybe removed in a future release.

The list of known particle scalars can be accessed by postpic.particle_scalars.

Examples
  • self('x')
  • self('sqrt(px**2 + py**2 + pz**2)')
Variable resolution order
  1. try to find the value as a atomic particle property.
  2. try to find the value as a defined particle property in particle_scalars.
  3. if not found look for an equally named attribute in scipy.constants.
__copy__()[source]

returns a shallow copy of the object. This method is called by copy.copy(obj).

__iadd__(other)[source]

Operator overloading for the += operator. Adding MultiSpecies should give the feeling as if you were adding their particle lists. That is why there is no append function. Compare those outputs (numpy.array handles that differently!):

>>> a=[1,2,3]; a.append([4,5]); print a
>>> [1,2,3,[4,5]]
>>> a=[1,2,3]; a += [4,5]; print a
>>> [1,2,3,4,5]

This function modifies the current object. Same behaviour as

>>> a = [1,2]; b = a; b+=[7,8]; print(a,b)
>>> [1, 2, 7, 8] [1, 2, 7, 8]

See also

__add__()

__invert__()[source]

invert the selection of particles.

add(dumpreader, species, ignore_missing_species=False)[source]

adds a species to this MultiSpecies. This function modifies the current Object and always returns None.

species can be a single species name

or a reserved name for collection of species, such as ions adds all available particles that are ions nonions adds all available particles that are not ions ejected noejected all

Optional arguments
--------
ignore_missing_species = False

set to True to ignore if the species is missing.

angle_xaxis()[source]

Deprecated since version unknown: The function angle_xaxis is deprecated. Use self(“angle_xaxis”) instead.

angle_xy()[source]

Deprecated since version unknown: The function angle_xy is deprecated. Use self(“angle_xy”) instead.

angle_xz()[source]

Deprecated since version unknown: The function angle_xz is deprecated. Use self(“angle_xz”) instead.

angle_yx()[source]

Deprecated since version unknown: The function angle_yx is deprecated. Use self(“angle_yx”) instead.

angle_yz()[source]

Deprecated since version unknown: The function angle_yz is deprecated. Use self(“angle_yz”) instead.

angle_zx()[source]

Deprecated since version unknown: The function angle_zx is deprecated. Use self(“angle_zx”) instead.

angle_zy()[source]

Deprecated since version unknown: The function angle_zy is deprecated. Use self(“angle_zy”) instead.

beta()[source]

Deprecated since version unknown: The function beta is deprecated. Use self(“beta”) instead.

betax()[source]

Deprecated since version unknown: The function betax is deprecated. Use self(“betax”) instead.

betay()[source]

Deprecated since version unknown: The function betay is deprecated. Use self(“betay”) instead.

betaz()[source]

Deprecated since version unknown: The function betaz is deprecated. Use self(“betaz”) instead.

charge()[source]

Deprecated since version unknown: The function charge is deprecated. Use self(“charge”) instead.

charge_e()[source]

Deprecated since version unknown: The function charge_e is deprecated. Use self(“charge_e”) instead.

compress(condition, name='unknown condition')[source]

works like numpy.compress. Additionaly you can specify a name, that gets saved in the compresslog. Returns a new MultiSpecies instance. compress gives you a lot of control, but in most cases filter() will be sufficient and keeps your code more readable.

Parameters:
  • condition (1D numpy array) –

    Condition can be one out of two choices:

    • condition = [True, False, True, True, … , True, False]

      condition is a list of length N, specifing which particles to keep. The length of the list must be equal to the length of the MultiSpecies instance, otherwise a ValueError is raised. Example:

      >>> cfintospectrometer = lambda ms: ms('abs(angle_xaxis) < 30e-3')
      >>> ms2 = ms.compress(cfintospectrometer(ms), name='< 30mrad offaxis')
      

      Consider using filter() instead.

    • condtition = [7, 2000, 4, 5, 91, … , 765, 809]

      condition can be a list of arbitraty length. Only the particles with the ids listed here will be kept.

  • name (str, optional) – an optional name of the condition. This can later be reviewed by calling ‘self.compresslog()’
compressfn(conditionf, name='unknown condition')[source]

like compress(), but accepts a function.

Returns a new MultiSpecies instance.

createField(*sps, **kwargs)[source]

Creates an n-d Histogram enclosed in a Field object.

Parameters:
  • *sps – list of scalarfunctions/strings/scalar-properties, that will be evaluated to data for each axis. the number of args given determins the dimensionality of the field returned by this function (maximum 3)
  • name (string, optional) – addes a name. usually used for generating savenames. Defaults to “distfn”.
  • title (string, options) – overrides the title. Autocreated if title==None. Defaults to None.
  • simgrid (boolean, optional) – enforces the same grid as used in the simulation. Implies simextent=True. Defaults to False.
  • simextent (boolean, optional) – enforces, that the axis show the same extent as used in the simulation. Defaults to False.
  • weights (function, optional) – applies additional weights to the macroparticles, for example ‘gamma’ or ‘q’ to weight the particle by its charge. Defaults to ‘1’ (no additional weight).
  • rangex (list of two values, optional) – the xrange to include into the histogram. Defaults to None, determins the range by the range of scalars given.
  • rangey (list of two values, optional) – the yrange to include into the histogram. Defaults to None, determins the range by the range of scalars given.
  • rangez (list of two values, optional) – the zrange to include into the histogram. Defaults to None, determins the range by the range of scalars given.
  • bins (sequence or int) – The number of bins to use for each dimension
  • shape (int) – possible choices are: * 0 - use nearest grid point (NGP) * 1 - use tophat shape of width 1 bin * 2 - triangular shape (default) * 3 - spline 3 shape
dumpreader

returns the dumpreader if the dumpreader of all species are pointing to the same dump. This should be mostly the case.

Otherwise returns None.

filter(condition, name=None)[source]

like compress(), but takes a ScalarProperty object instead which is required to evalute to a boolean list. Example:

>>> ms2 = ms.filter('gamma > 12')
Parameters:condition (str) – A string, which can also be used at ms(condition) and evaluates
gamma()[source]

Deprecated since version unknown: The function gamma is deprecated. Use self(“gamma”) instead.

gamma_m1()[source]

Deprecated since version unknown: The function gamma_m1 is deprecated. Use self(“gamma_m1”) instead.

getcompresslog()[source]
initial_npart

Original number of particles (before the use of compression or filter).

mass()[source]

Deprecated since version unknown: The function mass is deprecated. Use self(“mass”) instead.

mass_u()[source]

Deprecated since version unknown: The function mass_u is deprecated. Use self(“mass_u”) instead.

mean(expr, weights=None)[source]

The mean of a value given by the expression expr. The particle weight of the individual particles will be automatically included in the calculation. An additional weight can be given using the keyword weights.

median(expr, weights=None)[source]

The median of a value given by the expression expr. The particle weight of the individual particles will be automatically included in the calculation. An additional weight can be given using the keyword weights.

name

an alias to self.species

npart

Number of Particles.

nspecies

Number of species.

quantile(expr, q, weights=None)[source]

The qth-quantile of the distribution of a value given by the expression expr. q can be a scalar or a list of quantiles to be calculated simultaneously. The particle weight of the individual particles will be automatically included in the calculation. An additional weight can be given using the keyword weights.

r_xy()[source]

Deprecated since version unknown: The function r_xy is deprecated. Use self(“r_xy”) instead.

r_xyz()[source]

Deprecated since version unknown: The function r_xyz is deprecated. Use self(“r_xyz”) instead.

r_yz()[source]

Deprecated since version unknown: The function r_yz is deprecated. Use self(“r_yz”) instead.

r_zx()[source]

Deprecated since version unknown: The function r_zx is deprecated. Use self(“r_zx”) instead.

simextent(axis)[source]

the combined simextent for all species and dumps included in this MultiSpecies object.

simgridpoints(axis)[source]

this function is for convenience only and is likely to be removed in the future. Particlarly it is impossible to define the grid of the simulation if the MultiSpecies object consists of multiple dumps from different simulations.

species

returns an string name for the species involved. Basically only returns unique names from all species (used for plotting and labeling purposes – not for completeness). May be overwritten.

speciess

a complete list of all species involved.

time()[source]

Deprecated since version unknown: The function time is deprecated. Use self(“time”) instead.

uncompress()[source]

Returns a new MultiSpecies instance, with all previous calls of compress() or filter() undone.

var(expr, weights='1')[source]

The variance of a value given by the expression expr. The particle weight of the individual particles will be automatically included in the calculation. An additional weight can be given using the keyword weights.

weight()[source]

Deprecated since version unknown: The function weight is deprecated. Use self(“weight”) instead.

class postpic.ParticleHistory(sr, speciess, ids=None)[source]

Bases: object

Represents a list of particles including their history that can be found in all the dumps defined by the simulation reader sr.

Parameters:
  • sr (iterable of datareader) – a collection of datareader to use. Usually a Simulationreader object
  • speciess (string or iterable of strings) – a species name or a list of species names. Those particles can be included into the history.
  • ids (iterable of int) – list of ids to use (default: None). If this is None all particles in speciess will be tracked. If a list of ids is given, these ids will be serached in speciess only.
__copy__()[source]

returns a shallow copy of the object. This method is called by copy.copy(obj).

collect(*scalarfs)[source]

Collects the given particle properties for all particles for all times.

*scalarfs: the scalarfunction(s) defining the particle property

numpy.ndarray holding the different particles in the same order as the list of self.ids, meaning the particle on position particle_idx has the ID self.ids[particle_idx]. every array element holds the history for a single particle. Indexorder of returned array: [particle_idx][scalarf_idx, collection_idx]

skip(n)[source]

takes only everth (n+1)-th particle

postpic.histogramdd(data, **kwargs)[source]

Creates a histogram of the data. This function has the similar signature and return values as numpy.histogramdd. In addition this function supports the shape keyword argument to choose the particle shape used. If used with shape=0 the results of this function and the numpy.histogramdd are identical, however, this function is approx. factor 2 or 3 faster.

Parameters:
  • data (sequence of ndarray or ndarray (1D or 2D)) –
    The input (particle) data for the histogram.
    • A 1D numpy array (for 1D histogram).
    • A sequence providing the data for the different axis, i.e. (datax, datay, dataz) (preferred).
    • A (N, D)-array, i.e. [[x1, y1, z1], [x2, y2, z2]] – must be a numpy array!
  • bins (sequence or int) – The number of bins to use for each dimension
  • range (sequence, optional) – A sequence of lower and upper bin edges to be used if the edges are not given explicitly in bins. Defaults to the minimum and maximum values along each dimension.
  • weights (1D numpy array) – The weights to be used for each data point
  • shape (int) –
    possible choices are:
    • 0 - use nearest grid point (NGP)
    • 1 - use tophat shape of width 1 bin
    • 2 - triangular shape (default)
    • 3 - spline 3 shape
Returns:

  • H (ndarray) – the final histogram
  • edges (list) – A list of D arrays describing the edges for each dimension

class postpic.SpeciesIdentifier[source]

Bases: postpic.helper.PhysicalConstants

This Class provides static methods for deriving particle properties from species Names. The only reason for this to be a class is that it can be used as a mixin.

classmethod identifyspecies(species)[source]

Returns a dictionary containing particle informations deduced from the species name. The following keys in the dictionary will always be present: name species name string mass kg (SI) charge C (SI) tracer boolean ejected boolean

Valid Examples: Periodic Table symbol + charge state: c6, F2, H1, C6b ionm#c# defining mass and charge: ionm12c2, ionc20m110 advanced examples: ejected_tracer_ionc5m20b, ejected_tracer_electronx, ejected_c6b, tracer_proton, protonb

static isejected(species)[source]
classmethod ision(species)[source]
postpic.chooseCode(code)[source]

Chooses appropriate reader for the given simulation code. After choosing a preset of the correct reader, the functions postpic.readDump() and postpic.readSim() are setup for this preset.

Parameters:code (string) –
Possible options are:
  • ”DUMMY”: dummy class creating fake data.
  • ”EPOCH”: .sdf files written by EPOCH1D, EPOCH2D or EPOCH3D.
  • ”openPMD”: .h5 files written in openPMD Standard
  • ”piconGPU”: same as “openPMD”
  • ”VSIM”: .hdf5 files written by VSim.
postpic.readDump(dumpidentifier, **kwargs)[source]

After using the fucntion postpic.chooseCode(), this function should be the main function for reading a dump into postpic.

Parameters:
  • dumpidentifier (str) – Identifies the dump. For most dumpreaders this is a string poining to the file or folder. See what the specific reader of your format expects.
  • **kwargs – will be forwarded to the dumpreader.
Returns:

the dumpreader for this specific data dump.

Return type:

Dumpreader

postpic.readSim(simidentifier, **kwargs)[source]

After using the function postpic.chooseCode(), this function should be the main function for reading a simulation into postpic. A simulation is equivalent to a series of dumps in a specific order (not neccessarily time order).

Parameters:
  • simidentifier (str) – Identifies the simulation. For EPOCH this should be a string pointing to a .visit file. Specifics depend on the current simreader class, as set by chooseCode.
  • **kwargs – will be forwarded to the simreader.
Returns:

the Simulationreader

postpic.export_field(filename, field, **kwargs)[source]

export Field object as a file. Format depends on the extention of the filename. Currently supported are: .npz:

uses numpy.savez.
.csv:
uses numpy.savetxt.
.vtk:
vtk export to paraview
postpic.load_field(filename)[source]

Load a field object previously stored using the saveto method. These are .npz files with a specific metadata layout.

postpic.export_scalar_vtk(filename, scalarfield, **kwargs)[source]

exports one 2D or 3D scalar field object to a VTK file which is suitable for viewing in ParaView. It is assumed that all fields are defined on the same grid.

postpic.export_scalars_vtk(filename, *fields, **kwargs)[source]

exports a set of scalar fields to a VTK file suitable for viewing in ParaView. Up to four fields may be given

postpic.export_vector_vtk(filename, *fields, **kwargs)[source]

exports a vector field to a VTK file suitable for viewing in ParaView. Three 3D fields are expected, which will form the X, Y and Z component of the vector field. If less than tree fields are given, the missing components will be assumed to be zero.

Submodules

postpic.datahandling module

The Core module for final data handling.

This module provides classes for dealing with axes, grid as well as the Field class – the final output of the postpic postprocessor.

Terminology

A data field with N numeric points has N ‘grid’ points, but N+1 ‘grid_nodes’ as depicted here:

+---+---+---+---+---+
|   |   |   |   |   |
+---+---+---+---+---+
|   |   |   |   |   |
+---+---+---+---+---+
|   |   |   |   |   |
+---+---+---+---+---+
  o   o   o   o   o     grid      (coordinates where data is sampled at)
o   o   o   o   o   o   grid_node (coordinates of grid cell boundaries)
|                   |   extent
class postpic.datahandling.Field(matrix, name='', unit='', **kwargs)[source]

Bases: postpic._compat.mixins.NDArrayOperatorsMixin

The Field Object carries data in form of an numpy.ndarray together with as many Axis objects as the data’s dimensions. Additionaly the Field object provides any information that is necessary to plot _and_ annotate the plot.

Create a Field object from scratch. The only required argument is matrix which contains the actual data.

A name and a unit may be supplied.

The axis may be specified in different ways:

  • by passing a list of Axis object as axes
  • by passing arrays with the grid_nodes as xedges, yedges and zedges. This is intended to work with np.histogram.
  • by not passing anything, which will create default axes from 0 to 1.
T

Return the Field with the order of axes reversed. In 2D this is the usual matrix transpose operation.

__array__(dtype=None)[source]

will be called by numpy function in case an numpy array is needed.

__array_priority__ = 1
__array_ufunc__(ufunc, method, *inputs, **kwargs)[source]
__array_wrap__(array, context=None)[source]
__copy__()[source]

returns a shallow copy of the object. This method is called by copy.copy(obj). Just copy enough to create copies for operator overloading.

__getitem__(key)[source]
__setitem__(key, other)[source]
adjust_stagger_to(other)[source]
all(axis=None, out=None, keepdims=False)

Returns True if all elements evaluate to True.

Refer to numpy.all for full documentation.

See also

numpy.all()
equivalent function
angle
any(axis=None, out=None, keepdims=False)

Returns True if any of the elements of a evaluate to True.

Refer to numpy.any for full documentation.

See also

numpy.any()
equivalent function
atleast_nd(n)[source]

Make sure the field has at least ‘n’ dimensions

autocutout(axes=None, fractions=(0.001, 0.002))[source]

Automatically cuts out the main feature of the field by removing border regions that only contain small numbers.

This is done axis by axis. For each axis, the mean across all other axes is taken. The maximum max of the remaining 1d-array is taken and searched for the outermost boundaries a, d such that all values out of array[a:d] are smaller then fractions[0]*max. A second set of boundaries b, c is searched such that all values out of array[b:c] are smaller then fractions[1]*max. Because fractions[1] should be larger than fractions[0], array[b:c] should be contained completely in array[a:d].

A padding length x is chosen such that array[b-x:c+x] is entirely within array[a:d].

Then the corresponding axis of the field is sliced to [b-x:c+x] and multiplied with a tukey-window such that the region [b:c] is left untouched and the field in the padding region smoothly vanishes on the outer border.

This process is repeated for all axes in axes or for all axes if axes is None.

autoreduce(maxlen=4000)[source]

Reduces the Grid to a maximum length of maxlen per dimension by just executing half_resolution as often as necessary.

clip(a_min, a_max, out=None)[source]
conj()[source]
cutout(newextent)[source]

only keeps that part of the data, that belongs to newextent.

derivative(axis, staggered=False)[source]

Calculate the derivative of the field with respect to axis.

Uses np.gradient by default which outputs the second order accurate derivative on the same grid as the input field.

If staggered=True is passed, the method will instead calculate the second order accurate derivative at the points centered between the input grid points.

dimensions

returns only present dimensions. [] and [[]] are interpreted as -1 np.array(2) is interpreted as 0 np.array([1,2,3]) is interpreted as 1 and so on…

ensure_frequency_domain()[source]
ensure_positive_axes()[source]

ensures, that all axis are going from smaller to greater numbers, i.e. none of the axis is reversed.

ensure_spatial_domain()[source]
ensure_transform_state(transform_states)[source]

Makes sure that the field has the given transform_states. transform_states may be a single boolean, indicating the same desired transform_state for all axes. It may be a list of the desired transform states for all the axes or a dictionary indicating the desired transform states of specific axes.

evaluate(ex, local_dict=None, global_dict=None, **kwargs)[source]

Evaluates the expression ex using NumExpr and returns a field containing the result. This copies all metadata from self and just replaces the matrix.

This function is basically syntactic sugar simplifying

`field.replace_data(ne.evaluate(expr))`

to

`field.evaluate(expr)`

This method replicates some logic from NumExpr.necompiler.getArguments(), seems there is no way around it.

export(filename, **kwargs)[source]
Uses postpic.export_field to export this field to a file. All ``**kwargs` will be forwarded to this function. Format is recognized by the extension of the filename.

export Field object as a file. Format depends on the extention of the filename. Currently supported are: .npz:

uses numpy.savez.
.csv:
uses numpy.savetxt.
.vtk:
vtk export to paraview
extent

returns the extents in a linearized form, as required by “matplotlib.pyplot.imshow”.

fft(axes=None, exponential_signs='spatial', **kwargs)[source]

Performs Fourier transform on any number of axes.

The argument axis is either an integer indicating the axis to be transformed or a tuple giving the axes that should be transformed. Automatically determines forward/inverse transform. Transform is only applied if all mentioned axes are in the same transform state. If an axis is transformed twice, the origin of the axis is restored.

Parameters:
  • exponential_signs

    configures the sign convention of the exponential.

    • exponential_signs == ‘spatial’: fft using exp(-ikx), ifft using exp(ikx)
    • exponential_signs == ‘temporal’: fft using exp(iwt), ifft using exp(-iwt)
  • **kwargs – keyword-arguments are passed to the underlying fft implementation.
fft_autopad(axes=None, fft_padsize=<postpic.helper.FFTW_Pad object>)[source]

Automatically pad the array to a size such that computing its FFT using FFTW will be fast.

Parameters:fft_padsize (callable) –

The default for keyword argument fft_padsize is a callable, that is used to calculate the padded size for a given size.

By default, this uses fft_padsize=helper.fftw_padsize which finds the next larger “good” grid size according to what the FFTW documentation says.

However, the FFTW documentation also says: “(…) Transforms whose sizes are powers of 2 are especially fast.”

If you don’t worry about the extra padding, you can pass fft_padsize=helper.fft_padsize_power2 and this method will pad to the next power of 2.

flip(axis)[source]

functionality of numpy.flip.

field.flip(0) returns a postpic.Field object with the specified axis flipped. np.flip(field) returns only the numpy.ndarray and the axis information is lost.

grid
grid_nodes
half_resolution(axis)[source]

Halfs the resolution along the given axis by removing every second grid_node and averaging every second data point into one.

If there is an odd number of grid points, the last point will be ignored (that means, the extent will change by the size of the last grid cell).

Returns:the modified Field.
Return type:Field
imag
classmethod importfrom(filename, **kwargs)[source]

Construct a new field object from foreign data. Currently, the following file formats are supported explicitly: *.csv *.png

All other files will be opened using Pillow.

integrate(axes=None, method=<function simps>)[source]

Calculates the definite integral along the given axes.

Parameters:method (callable) –

Choose the method to use. Available options:

  • ’constant’
  • any function with the same signature as scipy.integrate.simps (default).
islinear()[source]
label
classmethod loadfrom(filename)[source]

Load a field object previously stored using the saveto method. These are .npz files with a specific metadata layout.

map_axis_grid(axis, transform, preserve_integral=True, jacobian_func=None)[source]

Transform the Field to new coordinates along one axis.

This function transforms the coordinates of one axis according to the function transform and applies the jacobian to the data.

Please note that no interpolation is applied to the data, instead a non-linear axis grid is produced. If you want to interpolate the data to a new (linear) grid, use the method map_coordinates() instead.

In contrast to map_coordinates(), the function transform is not used to pull the new data points from the old grid, but is directly applied to the axis. This reverses the direction of the transform. Therfore, in order to preserve the integral, it is necessary to divide by the Jacobian.

Parameters:
  • axis (int) – the index or name of the axis you want to apply transform to.
  • transform (callable) – the transformation function which takes the old coordinates as an input and returns the new grid
  • preserve_integral (bool) – Divide by the jacobian of transform, in order to preserve the integral.
  • jacobian_func (callable) – If given, this is expected to return the derivative of transform. If not given, the derivative is numerically approximated.
map_coordinates(newaxes, transform=None, complex_mode='polar', preserve_integral=True, jacobian_func=None, jacobian_determinant_func=None, **kwargs)[source]

Transform the Field to new coordinates.

Parameters:
  • newaxes (list) – The new axes of the new coordinates.
  • transform (callable) –

    a callable that takes the new coordinates as input and returns the old coordinates from where to sample the Field. It is basically the inverse of the transformation that you want to perform. If transform is not given, the identity will be used. This is suitable for simple interpolation to a new extent/shape. Example for cartesian -> polar:

    >>> def T(r, theta):
    >>>    x = r * np.cos(theta)
    >>>    y = r * np.sin(theta)
    >>>    return x, y
    

    Note that this function actually computes the cartesian coordinates from the polar coordinates, but stands for transforming a field in cartesian coordinates into a field in polar coordinates.

    However, in order to preserve the definite integral of the field, it is necessary to multiply with the Jacobian determinant of T.

    \[\tilde{U}(r, \theta) = U(T(r, \theta)) \cdot \det \frac{\partial (x, y)}{\partial (r, \theta)}\]

    such that

    \[\int_V \mathop{\mathrm{d}x} \mathop{\mathrm{d}y} U(x,y) = \int_{T^{-1}(V)} \mathop{\mathrm{d}r}\mathop{\mathrm{d}\theta} \tilde{U}(r,\theta)\,.\]
  • complex_mode

    The complex_mode specifies how to proceed with complex data.

    • complex_mode = ‘cartesian’ - interpolate real/imag part (fastest)
    • complex_mode = ‘polar’ - interpolate abs/phase If skimage.restoration is available, the phase will be unwrapped first (default)
    • complex_mode = ‘polar-no-unwrap’ - interpolate abs/phase Skip unwrapping the phase, even if skimage.restoration is available
  • preserve_integral (bool) –

    If True (the default), the data will be multiplied with the Jacobian determinant of the coordinate transformation such that the integral over the data will be preserved.

    In general, you will want to do this, because the physical unit of the new Field will correspond to the new axis of the Fields. Please note that Postpic, currently, does not automatically change the unit members of the Axis and Field objects, this you will have to do manually.

    There are, however, exceptions to this rule. Most prominently, if you are converting to polar coordinates, it depends on what you are going to do with the transformed Field. If you intend to do a Cartesian r-theta plot or are interested in a lineout for a single value of theta, you do want to apply the Jacobian determinant. If you had a density in e.g. J/m^2 than, in polar coordinates, you want to have a density in J/m/rad. If you intend, on the other hand, to do a polar plot, you do not want to apply the Jacobian. In a polar plot, the data points are plotted with variable density which visually takes care of the Jacobian automatically. A polar plot of the polar data should look like a Cartesian plot of the original data with just a peculiar coordinate grid drawn over it.

  • jacobian_determinant_func (callable) – A callable that returns the jacobian determinant of the transform. If given, this takes precedence over the following option.
  • jacobian_func (callable) – a callable that returns the jacobian of the transform. If this is not given, the jacobian is numerically approximated.
  • **kwargs – Additional keyword arguments are passed to scipy.ndimage.map_coordinates, see the documentation of that function.
matrix
max(axis=None, out=None)

Return the maximum along a given axis.

Refer to numpy.amax for full documentation.

See also

numpy.amax()
equivalent function
mean(axis=None, dtype=None, out=None, keepdims=False)

Returns the average of the array elements along given axis.

Refer to numpy.mean for full documentation.

See also

numpy.mean()
equivalent function
meshgrid(sparse=True)[source]
min(axis=None, out=None, keepdims=False)

Return the minimum along a given axis.

Refer to numpy.amin for full documentation.

See also

numpy.amin()
equivalent function
ndim
pad(pad_width, mode='constant', **kwargs)[source]

Pads the data using np.pad and takes care of the axes. See documentation of numpy.pad.

In contrast to np.pad, pad_width may be given as integers, which will be interpreted as pixels, or as floats, which will be interpreted as distance along the appropriate axis.

All other parameters are passed to np.pad unchanged.

phase(do_unwrap_phase=True)[source]

Returns the (unwrapped) phase of the complex Field.

do_unwrap_phase: True, if skimage.restoration.unwrap_phase should be applied to data

prod(axis=None, dtype=None, out=None, keepdims=False)

Return the product of the array elements over the given axis

Refer to numpy.prod for full documentation.

See also

numpy.prod()
equivalent function
ptp(axis=None, out=None)

Peak to peak (maximum - minimum) value along a given axis.

Refer to numpy.ptp for full documentation.

See also

numpy.ptp()
equivalent function
real
replace_data(other)[source]
rot90(k=1, axes=(0, 1))[source]

Rotates the field by 90 degrees, k times. Rotates the field in the plane given by axes.

saveto(filename)[source]

Save a Field object as a file. Use loadfrom() to load Field objects.

setaxisobj(axis, axisobj)[source]

replaces the current axisobject for axis axis by the new axisobj axisobj.

shape
shift_grid_by(dx, interpolation='fourier')[source]

Translate the Grid by dx. This is useful to remove the grid stagger of field components.

If all axis will be shifted, dx may be a list. Otherwise dx should be a mapping from axis to translation distance.

The keyword-argument interpolation indicates the method to be used and may be one of [‘linear’, ‘fourier’]. In case of interpolation = ‘fourier’ all axes must have same transform_state.

spacing

returns the grid spacings for all axis.

squeeze()[source]

removes axes that have length 1, reducing self.dimensions.

Note, that axis with length 0 will not be removed! numpy.squeeze also does not remove length=0 directions.

Same as numpy.squeeze.

std(axis=None, dtype=None, out=None, ddof=0, keepdims=False)

Returns the standard deviation of the array elements along given axis.

Refer to numpy.std for full documentation.

See also

numpy.std()
equivalent function
sum(axis=None, dtype=None, out=None, keepdims=False)

Return the sum of the array elements over the given axis.

Refer to numpy.sum for full documentation.

See also

numpy.sum()
equivalent function
swapaxes(axis1, axis2)[source]

Swaps the axes axis1 and axis2, equivalent to numpy.swapaxes.

topolar(extent=None, shape=None, angleoffset=0, **kwargs)[source]

Transform the Field to polar coordinates.

This is a convenience wrapper for map_coordinates() which will let you easily define the desired grid in polar coordinates.

Parameters:
  • extent – should be of the form extent=(phimin, phimax, rmin, rmax) or extent=(phimin, phimax)
  • shape – should be of the form shape=(N_phi, N_r),
  • angleoffset – can be any real number and will rotate the zero-point of the angular axis.
  • complex_mode

    The complex_mode specifies how to proceed with complex data.

    • complex_mode = ‘cartesian’ - interpolate real/imag part (fastest)
    • complex_mode = ‘polar’ - interpolate abs/phase If skimage.restoration is available, the phase will be unwrapped first (default)
    • complex_mode = ‘polar-no-unwrap’ - interpolate abs/phase Skip unwrapping the phase, even if skimage.restoration is available
  • preserve_integral (bool) –

    If True (the default), the data will be multiplied with the Jacobian determinant of the coordinate transformation such that the integral over the data will be preserved.

    In general, you will want to do this, because the physical unit of the new Field will correspond to the new axis of the Fields. Please note that Postpic, currently, does not automatically change the unit members of the Axis and Field objects, this you will have to do manually.

    There are, however, exceptions to this rule. Most prominently, if you are converting to polar coordinates, it depends on what you are going to do with the transformed Field. If you intend to do a Cartesian r-theta plot or are interested in a lineout for a single value of theta, you do want to apply the Jacobian determinant. If you had a density in e.g. J/m^2 than, in polar coordinates, you want to have a density in J/m/rad. If you intend, on the other hand, to do a polar plot, you do not want to apply the Jacobian. In a polar plot, the data points are plotted with variable density which visually takes care of the Jacobian automatically. A polar plot of the polar data should look like a Cartesian plot of the original data with just a peculiar coordinate grid drawn over it.

  • jacobian_determinant_func (callable) – A callable that returns the jacobian determinant of the transform. If given, this takes precedence over the following option.
  • jacobian_func (callable) – a callable that returns the jacobian of the transform. If this is not given, the jacobian is numerically approximated.
  • **kwargs – Additional keyword arguments are passed to scipy.ndimage.map_coordinates, see the documentation of that function.
transpose(*axes)[source]

transpose method equivalent to numpy.ndarray.transpose. If axes is empty, the order of the axes will be reversed. Otherwise axes[i] == j means that the i’th axis of the returned Field will be the j’th axis of the input Field.

var(axis=None, dtype=None, out=None, ddof=0, keepdims=False)

Returns the variance of the array elements, along given axis.

Refer to numpy.var for full documentation.

See also

numpy.var()
equivalent function
class postpic.datahandling.Axis(name='', unit='', **kwargs)[source]

Bases: object

Axis handling for a single Axis.

Create an Axis object from scratch.

The least required arguments are any of:
  • grid
  • grid_node
  • extent _and_ n

The remaining fields will be deduced from the givens.

More arguments may be supplied, as long as they are compatible.

__eq__(other)[source]

equality test for axis

__getitem__(key)[source]

Returns an Axis which consists of a sub-part of this object defined by a slice containing floats or integers or a float or an integer

__getstate__()[source]

Excludes self._inv_map from the pickled state

extent
grid
grid_node
half_resolution()[source]

removes every second grid_node.

islinear(force=False)[source]

Checks if the axis has a linear grid.

isreversed
label
physical_length
reversed()[source]

returns an reversed Axis object

spacing
value_to_index(value)[source]

This funtion is used to map values to indices in an interpolating manner, this is mainly used by the map_coordinates method of the Field class.

In contrast to the _find_nearest_index method, this method does not return an integer but a fractional index that refers to a position between actual pixels.

In general the equality

ax._find_nearest_index(x) == np.round(ax.value_to_index(x))

should hold.

postpic.experimental module

Some experimental algorithms for your reference. Please note that these algorithms are not meant to be used as is and may need adjustment in order to be applicable to a wider range of cases.

postpic.experimental.kspace_propagate_adaptive(field_in, axis=0, t_final=None, **kwargs)[source]

An adaptive method to use Fourier propagation (provided by the function helper.kspace_propagate) to get far field data. The field is padded, propagated and automatically sliced in repeating steps.

Note that this method is highly experimental and should not be trusted as is: It is merely meant as a recipe so you don’t have to write your own function from scratch!

field_in: input field in either spatial or frequency domain axis: The direction in which to propagate. Currently only propagation parallel to the positive x, y or z direction is implemented. yield_zeroth_step: boolean that determines if the initial step is also output.

t_final: The time at which to stop the adaptive propagation.

postpic.helper module

Some global constants that are used in the code.

class postpic.helper.PhysicalConstants[source]

Bases: object

gives you some constants.

c = 299792458.0
epsilon0 = 8.854187817620389e-12
mass_u = 1.67266490646e-27
me = 9.109383e-31
mu0 = 1.2566370614359173e-06
static ncrit(laslambda)[source]

Critical plasma density in particles per m^3 for a given wavelength laslambda in m.

static ncrit_um(lambda_um)[source]

Critical plasma density in particles per m^3 for a given wavelength lambda_um in microns.

qe = 1.602176565e-19
postpic.helper.unstagger_fields(*fields, **kwargs)[source]

Unstagger a collection of fields.

This functions shifts the origins of the grids of the given fields such that they coincide. Since the choice of the common origin is somewhat arbitrary, it might be overriden by a keyword-argument origin, as may be the interpolation method. See Field.shift_grid_by for available methods.

postpic.helper.kspace_epoch_like(component, fields, dt, extent=None, omega_func=<function omega_free>, align_to='B')[source]

Reconstruct the physical kspace of one polarization component See documentation of kspace

This function will use special care to make sure, that the implicit linear interpolation introduced by Epochs half-steps will not impede the accuracy of the reconstructed k-space. The frequency response of the linear interpolation is modelled and removed from the interpolated fields.

dt: time-step of the simulation, this is used to calculate the frequency response due to the linear interpolated half-steps

For the current version of EPOCH, v4.9, use the following: align_to == ‘B’ for intermediate dumps, align_to == “E” for final dumps

postpic.helper.kspace(component, fields, extent=None, interpolation=None, omega_func=<function omega_free>)[source]

Reconstruct the physical kspace of one polarization component This function basically computes one component of

E = 0.5*(E - omega/k^2 * Cross[k, B])
or
B = 0.5*(B + 1/omega * Cross[k, E]).

component must be one of [“Ex”, “Ey”, “Ez”, “Bx”, “By”, “Bz”].

The necessary fields must be given in the dict fields with keys chosen from [“Ex”, “Ey”, “Ez”, “Bx”, “By”, “Bz”]. Which are needed depends on the chosen component and the dimensionality of the fields. In 3D the following fields are necessary:

Ex, By, Bz -> Ex Ey, Bx, Bz -> Ey Ez, Bx, By -> Ez

Bx, Ey, Ez -> Bx By, Ex, Ez -> By Bz, Ex, Ey -> Bz

In 2D, components which have “k_z” in front of them (see cross-product in equations above) are not needed. In 1D, components which have “k_y” or “k_z” in front of them (see cross-product in equations above) are not needed.

The keyword-argument extent may be a list of values [xmin, xmax, ymin, ymax, …] which denote a region of the Fields on which to execute the kspace reconstruction.

The keyword-argument interpolation indicates whether interpolation should be used to remove the grid stagger. If interpolation is None, this function works only for non-staggered grids. Other choices for interpolation are “linear” and “fourier”.

The keyword-argument omega_func may be used to pass a function that will calculate the dispersion relation of the simulation may be given. The function will receive one argument that contains the k mesh.

postpic.helper.kspace_propagate(kspace, dt, nsteps=1, **kwargs)[source]

Evolve time on a field. This function checks the transform_state of the field and transforms first from spatial domain to frequency domain if necessary. In this case the inverse transform will also be applied to the result before returning it. This works, however, only correctly with fields that are the inverse transforms of a k-space reconstruction, i.e. with complex fields.

dt: time in seconds

This function will return an infinite generator that will do arbitrary many time steps.

If yield_zeroth_step is True, then the kspace will also be yielded after removing the antipropagating waves, but before the first actual step is done.

If a vector moving_window_vect is passed to this function, which is ideally identical to the mean propagation direction of the field in forward time direction, an additional linear phase is applied in order to keep the pulse inside of the box. This effectively enables propagation in a moving window. If dt is negative, the window will actually move the opposite direction of moving_window_vect. Additionally, all modes which propagate in the opposite direction of the moving window, i.e. all modes for which dot(moving_window_vect, k)<0, will be deleted.

The motion of the window can be inhibited by specifying move_window=False. If move_window is None, the moving window is automatically enabled if moving_window_vect is given.

The deletion of the antipropagating modes can be inhibited by specifying remove_antipropagating_waves=False. If remove_antipropagating_waves is None, the deletion of the antipropagating modes is automatically enabled if moving_window_vect is given.

nsteps: number of steps to take

If nsteps == 1, this function will just return the result. If nsteps > 1, this function will return a generator that will generate the results. If you want a list, just put list(…) around the return value.

postpic.helper.time_profile_at_plane(kspace_or_complex_field, axis='x', value=None, dir=1, t_input=0.0, **kwargs)[source]

‘Measure’ the time-profile of the propagating complex_field while passing through a plane.

The arguments axis, value and dir specify the plane and main propagation direction.

axis specifies the axis perpendicular to the measurement plane.

dir=1 specifies propagation towards positive axis, dir=-1 specifies the opposite direction of propagation.

value specifies the position of the plane along axis. If value=None, a default is chosen, depending on dir.

If dir=-1, the starting point of the axis is used, which lies at the 0-component of the inverse transform.

If dir=1, the end point of the axis + one axis spacing is used, which, via periodic boundary conditions of the fft, also lies at the 0-component of the inverse transform.

If the given value differs from these defaults, an initial propagation with moving window will be performed, such that the desired plane lies in the default position.

t_input specifies the point in time at which the input field or kspace is given. This is used to specify the time axis of the output fields.

For example axis=’x’ and value=0.0 specifies the ‘x=0.0’ plane while dir=1 specifies propagation towards positive ‘x’ values. The ‘x’ axis starts at 2e-5 and ends at 6e-5 with a grid spacing of 1e-6. The default value for the measurement plane would have been 6.1e-5 so an initial backward propagation with dt = -6.1e-5/c is performed to move the pulse in front of the’x=0.0 plane.

Additional kwargs are passed to kspace_propagate if they are not overridden by this function.