sdcc package

Contents

sdcc package#

Submodules#

sdcc.analysis module#

sdcc.analysis.analyze_hyst_data(vs, steps, d, hels, plot=False, ax=None)#
sdcc.analysis.get_critical_sizes(energyLandscape, Blocking_Ts)#
sdcc.analysis.process_thellier_data(vs_list, routine, weights)#

Processes data from a Thellier experiment to create Arai and Zijderveld plot data.

Inputs#

vs_list: list Set of moment vectors at every time step from an experiment run with the SDCC.

routine: list Routine of treatment steps used for the experiment.

weights: list or array Relative contribution of each grain to the result - a higher weight would mean relatively more of those grains.

returns:
  • Zs_list (array)

  • Array of moment vectors for Zijderveld plot data

  • Is_list (array)

  • In field step moment vectors (vector subtracted from

  • zero field vectors)

  • Zs_mag, Is_mag (array)

  • Magnitudes of above vectors, normalized by the TRM

  • magnitude (Arai plot data).

sdcc.barriers module#

class sdcc.barriers.GEL(TMx, alignment, PRO, OBL, T_spacing=1)#

Bases: object

Class for storing energy barrier results at all temperatures for a given grain geometry and composition.

Todo: 1. This should inherit from some base class shared with Hysteresis. 2. The run through of temperatures should be parallelized for much quicker object creation.

Parameters:
  • TMx (float) – Titanomagnetite composition (0 - 100)

  • alignment (str) – Alignment of magnetocrystalline and shape axis. Either ‘hard’ or ‘easy’ magnetocrystalline always aligned with shape easy.

  • PRO (float) – Prolateness of ellipsoid (major / intermediate axis)

  • OBL (float) – Oblateness of ellipsoid (intermediae / minor axis)

  • T_spacing (float) – Spacing of temperature steps (Degrees C)

get_params(T)#

Gets directions and energies associated with LEM states and barriers for a grain as a function of temperature.

Parameters:

T (int, float or numpy array) – Temperature(s) (degrees C)

Returns:

params – Dictionary of arrays for directions and energies.

Return type:

dict

to_file(fname)#

Saves GEL object to file.

Parameters:

fname (string) – Filename

Return type:

None

class sdcc.barriers.HEL(TMx, alignment, PRO, OBL, rot_mat, B_max, B_spacing=0.001, T=20)#

Bases: object

Class for storing energy barrier results at all field strengths for a given grain geometry, composition and field direction.

Todo: 1. This should inherit from some base class shared with GEL. 2. The run through of temperatures should be parallelized for much quicker object creation.

Parameters:
  • TMx (float) – Titanomagnetite composition (0 - 100)

  • alignment (str) – Alignment of magnetocrystalline and shape axis. Either ‘hard’ or ‘easy’ magnetocrystalline always aligned with shape easy.

  • PRO (float) – Prolateness of ellipsoid (major / intermediate axis)

  • OBL (float) – Oblateness of ellipsoid (intermediae / minor axis)

  • B_dir (length 3 array) – Unit vector of field direction

  • B_max (float) – Maximum field (Tesla)

  • B_spacing (float) – Spacing of field steps (Tesla)

get_params(B)#

Gets directions and energies associated with LEM states and barriers for a grain as a function of field.

Parameters:

B (int, float or numpy array) – External Field (T)

Returns:

params – Dictionary of arrays for directions and energies.

Return type:

dict

to_file(fname)#

Saves GEL object to file.

Parameters:

fname (string) – Filename

Return type:

None

class sdcc.barriers.HELs(TMx, alignment, PRO, OBL, B_max, B_spacing=0.001, T=20, rot_mats=None, n_dirs=30, HEL_list=None)#

Bases: object

to_file(fname)#

Saves GEL object to file.

Parameters:

fname (string) – Filename

Return type:

None

sdcc.barriers.blocking_temperature(gel: GEL, d, i, j, block_t=100.0)#

Calculates the blocking temperature associated with an energy barrier. Involves a minimization to obtain the correct time.

Parameters:
  • gel (GEL object) – Energy landscape of particle to be considered

  • d (float) – Size of particle (nm)

  • i (int) – index of first state in barrier

  • j (int) – index of second state in barrier

  • block_t (float) – Timescale at which the barrier is considered ‘blocked’

Returns:

T – Blocking temperature

Return type:

float

sdcc.barriers.construct_energy_mat_fast(thetas, phis, energies, labels)#

Constructs an energy matrix from a set of labels and energies Does so in a rapid fashion by finding the minimum along an edge between two minima.

Parameters:
  • thetas (numpy array) – grid of theta angles

  • phis (numpy array) – grid of phi angles

  • energies (numpy array) – grid of energies corresponding to theta, phi angles

  • labels (numpy array) – grid of integer labels for each drainage basin

Returns:

  • theta_mat (numpy array) – matrix of energy barrier thetas between state i and j

  • phi_mat (numpy array) – matrix of energy barrier phis between state i and j

  • energy_mat (numpy array) – matrix of energy barriers between state i and state j

sdcc.barriers.delete_splits(new_indices, old_indices, high_energy_list, high_theta_list, high_phi_list, theta_lists, phi_lists)#

Function that deletes “new” states. The function differentiates between “new” states that have appeared before and then disappeared due to numerical noise, so called “lazarus” states, and “new” states that appear to occur from a splitting of existing states - called “spurious” states, as the number of states is not expected to increase with temperature.

Parameters:
  • new_indices (numpy array) – Array of “new” set of indices of states

  • old_indices (numpy array) – Array of “old” set of indices of states

  • high_energy_list (numpy array) – Energies at higher temperature or field.

Returns:

new_indices,old_indices – Arrays of new and old indices with “new” states deleted

Return type:

numpy arrays

sdcc.barriers.dijkstras(i, j, barriers)#

A modified version of dijkstras algorithm used to “prune” energy barriers between states which already have a smaller barrier between them.

Parameters:#

i, j: ints

Indices for barriers

barriers: numpy array

matrix of energy barriers

Returns:#

barriers: numpy array

barrier array with steps pruned

sdcc.barriers.direction_result(t, c, k, T)#

Calculates a direction from a temperature and a set of spline coefficients.

Parameters:
  • t (numpy array) – Cubic spline knot positions

  • c (numpy array) – B-Spline coefficients

  • k (numpy array) – Polynomial degrees

  • T (float) – Temperature (degrees C).

Returns:

direction – Direction at that temperature.

Return type:

numpy array or float

sdcc.barriers.direction_spline(x, y)#

Creates a piecewise 3D unit vector B-Spline from a set of x,y data.

Parameters:
  • x (numpy array) – Input x data for cubic spline (usually temperature)

  • y (numpy array of 3D unit vectors.) – Input y data for cubic spline (usually a direction)

Returns:

  • t (numpy array) – Cubic spline knot positions

  • c (numpy array) – B-Spline coefficients

  • k (numpy array) – Polynomial degrees

sdcc.barriers.energy_result(t, c, k, T)#

Calculates an energy from a temperature and a set of spline coefficients.

Parameters:
  • t (numpy array) – Cubic spline knot positions

  • c (numpy array) – B-Spline coefficients

  • k (numpy array) – Polynomial degrees

  • T (float) – Temperature (degrees C).

Returns:

result – Energy at that temperature.

Return type:

numpy array or float

sdcc.barriers.energy_spline(x, y)#

Creates a piecewise scalar B-Spline from a set of x,y data.

Parameters:
  • x (numpy array) – Input x data for cubic spline (usually temperature)

  • y (numpy array) – Input y data for cubic spline (usually energy or magnetization)

Returns:

  • t (numpy array) – Cubic spline knot positions

  • c (numpy array) – B-Spline coefficients

  • k (numpy array) – Polynomial degrees

sdcc.barriers.find_B_barriers(TMx, alignment, PRO, OBL, B_dir, B_max, B_spacing, T=20)#

Finds all LEM states and energy barriers for a grain composition and geometry over a range of fields in a specfic direction. Runs from zero-field to B_max. This function could also easily be parallelized using multiprocessing.

Parameters:
  • TMx (float) – Titanomagnetite composition (0 - 100)

  • alignment (str) – Alignment of magnetocrystalline and shape axis. Either ‘hard’ or ‘easy’ magnetocrystalline always aligned with shape easy.

  • PRO (float) – Prolateness of ellipsoid (major / intermediate axis)

  • OBL (float) – Oblateness of ellipsoid (intermediae / minor axis)

  • B_dir (length 3 array of floats) – Cartesian direction (unit vector) of field.

  • B_max (float) – Maximum field (T).

  • B_spacing (float) – Spacing of field steps (T)

  • T (float) – Temperature (Degrees C) - must be between 0 and 80.

Returns:

  • theta_lists, phi lists (numpy arrays) – Arrays of the minimum energy directions in the SD energy surface

  • min_energy_lists (numpy array) – Minimum energy at the minimum

  • theta_mats,phi_mats (numpy arrays) – Arrays of the directions associated with the energy barriers.

  • energy_mats (numpy array) – Matrix of energy barriers between state i and state j.

  • Ts (numpy array) – List of temperatures.

sdcc.barriers.find_T_barriers(TMx, alignment, PRO, OBL, T_spacing=1)#

Finds all LEM states and energy barriers for a grain composition and geometry at all temperatures. Runs from room temperature (20 C) to Tc. Todo: Write one of these for Hysteresis. This function could also easily be parallelized using multiprocessing.

Parameters:
  • TMx (float) – Titanomagnetite composition (0 - 100)

  • alignment (str) – Alignment of magnetocrystalline and shape axis. Either ‘hard’ or ‘easy’ magnetocrystalline always aligned with shape easy.

  • PRO (float) – Prolateness of ellipsoid (major / intermediate axis)

  • OBL (float) – Oblateness of ellipsoid (intermediae / minor axis)

  • T_spacing (float) – Spacing of temperature steps (Degrees C)

Returns:

  • theta_lists, phi lists (numpy arrays) – Arrays of the minimum energy directions in the SD energy surface

  • min_energy_lists (numpy array) – Minimum energy at the minimum

  • theta_mats,phi_mats (numpy arrays) – Arrays of the directions associated with the energy barriers.

  • energy_mats (numpy array) – Matrix of energy barriers between state i and state j.

  • Ts (numpy array) – List of temperatures.

sdcc.barriers.find_all_barriers(TMx, alignment, PRO, OBL, T=20, ext_field=array([0, 0, 0]), prune=True, trim=True, tol=2.0, return_labels=False)#

Finds all the minimum energy states in an SD grains and the barriers between them.

Parameters:
  • TMx (float) – Titanomagnetite composition (0 - 100)

  • alignment (str) – Alignment of magnetocrystalline and shape axis. Either ‘hard’ or ‘easy’ magnetocrystalline always aligned with shape easy.

  • PRO (float) – Prolateness of ellipsoid (major / intermediate axis)

  • OBL (float) – Oblateness of ellipsoid (intermediae / minor axis)

  • T (float) – Temperature (degrees C)

  • ext_field (numpy array) – Array of field_theta,field_phi,field_str where theta and phi are in radians and str is in Tesla.

  • prune (boolean) – if True, uses dijkstras algorithm to prune extraneous barriers from the result.

  • trim (boolean) – if True, removes minima with more barriers than the number of magnetocrystalline easy directions. Can cause problems if tol is set too low, because some close together minima may be chosen instead of distinct ones, removing legitimate minima.

  • tol (float) – Angular distance in degrees between minima below which they’re considered “the same” state. Used to fix numerical errors where multiple minima are found in extremely flat regions of the energy surface.

  • return_labels (boolean) – If True, returns the “labels” array. Default behavior is False.

Returns:

  • theta_list, phi list (numpy arrays) – Arrays of the minimum energy directions in the SD energy surface

  • min_energy_list (numpy array) – Minimum energy at the minimum

  • theta_mat,phi_mat (numpy arrays) – Arrays of the directions associated with the energy barriers.

  • barriers (numpy array) – Matrix of energy barriers between state i and state j.

sdcc.barriers.find_global_minimum(thetas, phis, energies, mask=None)#

Finds the global minimum of a masked region of a theta-phi map

Parameters:
  • thetas (numpy array) – grid of theta angles

  • phis (numpy array) – grid of phi angles

  • energies (numpy array) – grid of energies corresponding to theta, phi angles

  • mask (numpy array) – grid of True or False indicating the location within the grid which should be evaluated to find the minimum

Returns:

  • best_coords (numpy array) – Indices of minima

  • best_theta (float) – Theta value of minima

  • best_phi (float) – Phi value of minima

  • best_energy (float) – Energy value of minima

sdcc.barriers.fix_minima(min_coords, energies, max_minima)#

Reduces a set of minima on an energy surface to some nominal “maximum number” (usually set by the magnetocrystalline anisotropy) Trims out the lowest minima and makes sure things on the border wrap for good minima.

Parameters:
  • min_coords (list) – List of indices of LEM values on our grid

  • energies (numpy array) – SD energy surface as a function of direction

  • max_minima (int) – Nominal maximum number of minima

Returns:

new_markers – Corrected indices of energy minima.

Return type:

numpy array

sdcc.barriers.get_antipode_order(theta_list, phi_list)#

Calculates the nearest neighbor between a set of LEM states and the set of LEM states located in an antipodal orientation.

Parameters:
  • theta_list (numpy arrays) – Directions of minima on energy surface at first temperature.

  • phi_list (numpy arrays) – Directions of minima on energy surface at first temperature.

Returns:

antipode_indices – indices of numpy arrays

Return type:

numpy array

sdcc.barriers.get_min_regions(energies, markers=None)#

Partitions an SD energy surface into a set of “drainage basins” each of which contain an LEM state. This uses the watershed algorithm in scipy. We tile the maps so the connectivity constraints for this algorithm work. The minimum points along the drainage divides between two watersheds will later be used as the energy barriers.

Parameters:
  • energies (numpy array) – SD energy surface

  • markers (None or numpy array) – If supplied, uses “markers” as the starting minima to calculate the drainage divides

Returns:

  • min_coords (numpy array) – indices of the minima in each region (should correspond to markers if this was supplied.

  • labels (numpy array) – array of labelled drainage divides corresponding to the

sdcc.barriers.get_minima(thetas, phis, energies, labels)#

Gets the minimum energies for each label in an array with labelled regions.

Parameters:
  • thetas (numpy array) – grid of theta angles

  • phis (numpy array) – grid of phi angles

  • energies (numpy array) – grid of energies corresponding to theta, phi angles

  • labels (numpy array) – grid of integer labels for each drainage basin

Returns:

theta_coords, phi_coords – arrays of theta and phi coordinates of each minimum.

Return type:

arrays

sdcc.barriers.make_antipode_array(Bs, theta_lists, phi_lists, min_energy_lists, theta_mats, phi_mats, energy_mats)#

Takes a set of LEM states and energy barriers calculated a range of fields and calculates the barriers for a range of antipodal fields. This works because all other anisotropies are axisymmetric.

Parameters:
  • Bs (list) – List of field strengths

  • theta_lists (list) – List of state theta coordinates at each field

  • phi_lists (list) – List of state phi coordinates at each field

  • min_energy_lists (list) – List of state energies at each field

  • theta_mats (list) – List of barrier theta coordinates at each field

  • phi_mats (list) – List of barrier phi coordinates at each field

  • energy_mats (list) – List of barrier energies at each field

Returns:

  • Bs_new (list) – List of field strengths

  • theta_lists_new (list) – List of state theta coordinates at each field

  • phi_lists_new (list) – List of state phi coordinates at each field

  • min_energy_lists_new (list) – List of state energies at each field

  • theta_mats_new (list) – List of barrier theta coordinates at each field

  • phi_mats_new (list) – List of barrier phi coordinates at each field

  • energy_mats_new (list) – List of barrier energies at each field

sdcc.barriers.mat_to_mask(theta_mat, phi_mat)#

Creates a mask for use in plot_energy_barriers

Parameters:
  • theta_mat (numpy arrays) – Matrices of directions associated with energy barriers

  • phi_mat (numpy arrays) – Matrices of directions associated with energy barriers

Returns:

mask – Array of bools showing where the energy minima are.

Return type:

numpy array

sdcc.barriers.merge_similar_minima(min_coords, thetas, phis, energies, tol)#

Merges together close together LEM states assuming they’re the same state. Takes the state with lowest energy as the “true” minimum.

Parameters:
  • min_coords (list) – List of minima coordinate values.

  • thetas (numpy arrays) – Grid of theta,phi values on which energies are evaluated.

  • phis (numpy arrays) – Grid of theta,phi values on which energies are evaluated.

  • energies (numpy array) – Energies at theta, phi locations.

  • tol (float) – Angular distance in degrees, below which two states considered the same.

Returns:

  • min_coords (list) – List of minima index values:

  • markers (numpy array) – Locations of these minima in the energy matrix.

sdcc.barriers.prune_energies(barriers, thetas, phis)#

Uses modified dijkstras algorithm to “prune” extraneous barriers

Parameters:
  • barriers (numpy array) – Matrix of energy barriers

  • thetas (numpy array) – Theta locations of barriers

  • phis (numpy array) – Phi locations of barriers

Returns:

  • temp_barriers (numpy array) – “Pruned” barrier array

  • thetas (numpy array) – “Pruned” theta array

  • phis (numpy array) – “Pruned” phi array

sdcc.barriers.segment_region(mask)#

Splits an energy array into a segmented regions according to a mask Produces a set of labels for the array. Creates regions of “drainage basins” for the array.

Parameters:

mask (numpy array) – mask of regions to segment

Returns:

labels – array of drainage basin labels.

Return type:

numpy array

sdcc.barriers.smart_sort(low_theta_list, low_phi_list, low_energy_list, low_theta_mat, low_phi_mat, low_barriers, low_labels, high_theta_list, high_phi_list, high_energy_list, high_theta_mat, high_phi_mat, high_barriers, high_labels, theta_lists, phi_lists)#

Function for tracking LEM states or energy barriers as a function of temperature or field. Uses the final set of labels from the watershed algorithm to calculate what each state would minimize to.

Parameters:
  • low_theta_list (length n numpy array) – Array of theta coordinates for lower temperature or field states

  • low_phi_list (length n numpy array) – Ditto, for phi coordinates

  • low_energy_list (length n numpy array) – Ditto for energies

  • low_theta_mat (n x n array) – Ditto for energy barrier theta coordinate

  • low_phi_mat (n x n array) – Ditto for energy barrier phi coordinate

  • low_barriers (n x n array) – Ditto for energy barriers

  • low_labels (1001 x 1001 array) – Labels for low temperature/field energy barriers

  • high_theta_list (length m numpy array) – Array of theta coordinates for lower temperature or field states

  • high_phi_list (length m numpy array) – Ditto, for phi coordinates

  • high_energy_list (length m numpy array) – Ditto for energies

  • high_theta_mat (m x m array) – Ditto for energy barrier theta coordinate

  • high_phi_mat (m x m array) – Ditto for energy barrier phi coordinate

  • high_barriers (m x m array) – Ditto for energy barriers

  • high_labels (1001 x 1001 array) – Labels for high temperature/field energy barriers

  • theta_lists (list) – list of all previous LEM theta coordinates

  • phi_lists (list) – list of all previous LEM phi coordinates

Returns:

  • thetas_new (length m numpy array) – Reindexed high_theta_list

  • phis_new (length m numpy array) – Reindexed high_phi_list

  • energy_list_new (length m numpy array) – Reindexed high_energy_list

  • theta_mat_new (m x m array) – Reindexed high_theta_mat

  • phi_mat_new (m x m array) – Reindexed high_phi_mat

  • barriers_new (m x m array) – Reindexed high_barriers

  • labels_new (1001 x 1001 array) – Reindexed high_labels

sdcc.barriers.uniaxial_critical_size(K, T, t=100)#

Calculates the critical size of a particle assuming only a uniaxial transition time

Parameters:
  • K (float) – Energy barrier (J/m^3)

  • T (float) – Temperature energy barrier is calculated at

  • t (float) – Desired timescale to target

Returns:

d – Grain equivalent sphere volume diameter (nm)

Return type:

float

sdcc.barriers.uniaxial_relaxation_time(d, T, K)#

Calculates the Neel relaxation time for an energy barrier at a particular temperature

Parameters:
  • d (float) – Equivalent spherical volume diameter of particle

  • T (float) – Temperature

  • K (float) – Energy density of energy barrier

Returns:

t – Relaxation time

Return type:

float

sdcc.barriers.wrap_labels(labels)#

Wraps a set of labels assigned to an array so that the labels in the first and last columns are the same - this is useful as the arrays we use wrap at theta = 0 and 2*pi.

Parameters:

Labels (array) – Array of integer labels for the “drainage basins” for each LEM

Returns:

Labels – Array of wrapped labels.

Return type:

array

sdcc.energy module#

sdcc.energy.Ea(k1, k2, theta, phi, rot_mat)#

Calculates the CUBIC magnetocrystalline anisotropy energy as a function of theta, phi angles. Anisotropy energy field can be rotated using a rotation matrix.

Parameters:
  • k1 (floats) – Cubic magnetocrystalline anisotropy constants. If k2 is not calculated, leave as 0.

  • k2 (floats) – Cubic magnetocrystalline anisotropy constants. If k2 is not calculated, leave as 0.

  • theta (floats) – Angles in spherical coordinates along which the energy is calculated (in radians).

  • phi (floats) – Angles in spherical coordinates along which the energy is calculated (in radians).

  • rot_mat (numpy array) – Rotation matrix to rotate anisotropy energy field by.

Returns:

Ea – Magnetocrystalline anisotropy energy.

Return type:

float

sdcc.energy.Ed(LMN, theta, phi, Ms)#

Calculates the demagnetizing field energy for an ellipsoid as a function of theta, phi angles and saturation magnetization. Only applicable to single domain grains.

Parameters:
  • LMN (numpy array) – Demagnetizing factors of the ellipsoid along x, y, z directions.

  • theta (float) – horizontal angle (radians) between 0 and 2*pi

  • phi (float) – vertical angle (radians) between -pi and pi

  • Ms (float) – Saturation magnetization in A/m

Returns:

Ed – Demagnetizing field energy in J/m3

Return type:

float

sdcc.energy.Ez(theta, phi, field_theta, field_phi, field_str, Ms)#

Calculates the Zeeman energy as a function of theta, phi angles.

Parameters:
  • theta (field) – Angles along which energy is evaluated (radians).

  • phi (field) – Angles along which energy is evaluated (radians).

  • theta – Direction of field specified as an angle (radians).

  • phi – Direction of field specified as an angle (radians).

  • str (field) – Strength of field (in Tesla).

  • Ms (float) – Saturation magnetization of sample.

Returns:

Ez – Zeeman energy.

Return type:

float

sdcc.energy.angle2xyz(theta, phi)#

Converts from coordinates on a sphere surface (theta, phi) radians to cartesian coordinates (x,y,z) as a unit vector

Parameters:
  • theta (float) – horizontal angle (radians) between 0 and 2*pi

  • phi (float) – vertical angle (radians) between -pi and pi

Returns:

xyz – array of the x,y and z cartesian coordinates

Return type:

jax numpy array

sdcc.energy.calculate_anisotropies(TMx)#

Calculates the room temperature easy, intermediate and hard directions as a function of titanomagnetite composition. This function could probably be handled by materials.py in the future

Parameters:

TMx (float) – Titanomagnetite composition (0 to 100).

Returns:

sorted_axes – “Special Directions” array of [easy, intermediate, hard] axes. Direction is given as a space separated string e.g. ‘1 0 0’ for compatibility with MERRILL script.

Return type:

numpy array

sdcc.energy.demag_factors(PRO, OBL)#

Calculates the demagnetizing factors for a standard ellipsoid using the formulae of Osborn (1945). We assume that the major axis of the ellipsoid is always along x and the minor axis of the ellipsoid is always along z.

Parameters:
  • PRO (float) – Prolateness ratio (major axis/intermediate axis)

  • OBL (float) – Oblateness ratio (intermediate axis/minor axis)

Returns:

LMN – Array of L, M, N, corresponding to the demagnetizing factors

Return type:

numpy array

sdcc.energy.dir_to_rot_mat(x, x_prime)#

Creates a rotation matrix that rotates from one crystallographic direction to another.

Parameters:
  • x (string) – vector to be rotated from. Specified as a space separated string e.g. ‘1 0 0’ for compatibility with MERRILL script.

  • x_prime (string) – vector to be rotated to. Specified in the same way.

  • Returns

  • rot_mat (numpy array) – Rotation matrix.

sdcc.energy.energy_ang(angles, k1, k2, rot_mat, LMN, Ms, ext_field)#

Calculates the total energy (Zeeman + Anisotropy + Demag) as a function of theta, phi angle. Requires material parameters and demagnetizing factors.

Parameters:
  • angles (array) – Array of angles specified as [theta, phi]

  • k1 (floats) – Magnetocrystalline anisotropy constant

  • k2 (floats) – Magnetocrystalline anisotropy constant

  • rot_mat (numpy array) – Rotation matrix for magnetocrystalline anisotropy field.

  • LMN (numpy array) – Array of demagnetizing factors for Ellipsoid shape.

  • Ms (float) – Saturation magnetization of material

  • ext_field (numpy array) – Array of [theta, phi, B] for the external field, where theta and phi are expressed in DEGREES and B is in Tesla.

Returns:

Etot – Total energy as a function of theta, phi angle.

Return type:

float

sdcc.energy.energy_surface(k1, k2, rot_mat, Ms, LMN, ext_field, n_points=100, bounds=Array([[0., 6.28318531], [-1.57079633, 1.57079633]], dtype=float64))#

Calculates the total SD energy for theta, phi angles along an equirectangular grid. Note that this is not an area preserving grid which may have some numerical effects - perhaps a UV map would be better.

N.B. This has a lot of arguments, maybe specify a dictionary?

Parameters:
  • k1 (floats) – Magnetocrystalline anisotropy constants

  • k2 (floats) – Magnetocrystalline anisotropy constants

  • rot_mat (numpy array) – Rotation matrix for magnetocrystalline anisotropy field.

  • Ms (float) – Saturation magnetization.

  • LMN (numpy array) – Demagnetizing factors of Ellipsoid along x, y, z directions.

  • ext_field (numpy array) – Array of [theta, phi, B] for the external field, where theta and phi are expressed in DEGREES and B is in Tesla.

  • n_points (int) – Number of points in theta, phi direction to calculate energy.

  • bounds (numpy array) – Theta, phi bounds within which to calculate the energies. Specified as [[theta_min,theta_max],[phi_min,phi_max]]

Returns:

  • thetas (numpy array) – grid of theta angles

  • phis (numpy array) – grid of phi angles

  • energies (numpy array) – grid of energies associated with these thetas and phis.

sdcc.energy.energy_xyz(xyz, k1, k2, rot_mat, LMN, Ms, ext_field)#

Calculates the total energy (Zeeman + Anisotropy + Demag) as a function of x, y, z coordinate. Requires material parameters and demagnetizing factors.

Parameters:
  • xyz (array) – Array of coordinates specified as [x,y,z]

  • k1 (floats) – Magnetocrystalline anisotropy constant

  • k2 (floats) – Magnetocrystalline anisotropy constant

  • rot_mat (numpy array) – Rotation matrix for magnetocrystalline anisotropy field.

  • LMN (numpy array) – Array of demagnetizing factors for Ellipsoid shape.

  • Ms (float) – Saturation magnetization of material

  • ext_field (numpy array) – Array of [theta, phi, B] for the external field, where theta and phi are expressed in DEGREES and B is in Tesla.

Returns:

Etot – Total energy as a function of theta, phi angle.

Return type:

float

sdcc.energy.get_material_parms(TMx, alignment, T)#

Todo: Incorporate materials.py framework into material parameters Calculates the material parameters for titanomagnetites of different compositions.

Parameters:
  • TMx (float) – Titanomagnetite composition % (0 - 100)

  • alignment (string) – Either easy or hard. Specifies the magnetocrystalline direction that should be aligned with the x direction, which for our ellipsoids is the major (shape easy) axis.

  • T (float) – Temperature (degrees C).

Returns:

  • rot_mat – Rotation matrix associated with magnetocrystalline anisotropy energy.

  • K1, K2 (floats) – magnetocrystalline anisotropy constants.

  • Ms (float) – Saturation magnetization.

sdcc.energy.xyz2angle(xyz)#

Converts cartesian coordinates to coordinates on a sphere surface. Vector length is lost in conversion.

Parameters:
  • x (float) – cartesian x coordinate

  • y (float) – cartesian y coordinate

  • z (float) – cartesian z coordinate

Returns:

  • theta (float) – horizontal angle (radians) between 0 and 2*pi

  • phi (float) – vertical angle (radians) between -pi and pi

sdcc.materials module#

class sdcc.materials.Anisotropy(anisotropyDict, solidSolution=False)#

Bases: object

make_functions(param)#
class sdcc.materials.CubicAnisotropy(k1, k2)#

Bases: Anisotropy

class sdcc.materials.Material(materialDict, anisotropy)#

Bases: object

Aex(T)#
Kd(T)#
Ms(T)#
anisotropy_constants(T)#
exch_len(T)#
mesh_sizing(ESVD, Tmin, Tmax, availSizes)#
min_exch_len(Tmin, Tmax)#
class sdcc.materials.SolidSolution(solidSolutionDict, anisotropy=None)#

Bases: object

composition(comp)#
sdcc.materials.TM_k1(comp, T, Tc)#
sdcc.materials.TM_k2(comp, T, Tc)#

sdcc.plotting module#

sdcc.plotting._descent_update(i, j, energies)#

Computes the minimum energy in a 10x10 grid around energies[i,j], goes to that energy, like a fast gradient descent except not using gradients.

sdcc.plotting.dimap(D, I)#

Function to map directions to x,y pairs in equal area projection. This function is reproduced from the PmagPy Project (PmagPy/Pmagpy).

Copyright (c) 2023, PmagPy contributors All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

  • Redistributions of source code must retain the above copyright

notice, this list of conditions and the following disclaimer.

  • Redistributions in binary form must reproduce the above copyright

notice,this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

  • Neither the name of PmagPy nor the names of its contributors may

be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Parameters:
  • D (list or array of declinations (as float))

  • I (list or array or inclinations (as float))

Returns:

XY

Return type:

x, y values of directions for equal area projection [x,y]

sdcc.plotting.dimap_V(D, I)#

Maps declinations and inclinations into equal area projections. This function is a modified version of the dimap function from pmagpy (PmagPy/PmagPy).

Copyright (c) 2023, PmagPy contributors All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

  • Redistributions of source code must retain the above copyright

notice, this list of conditions and the following disclaimer.

  • Redistributions in binary form must reproduce the above copyright

notice,this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

  • Neither the name of PmagPy nor the names of its contributors may

be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Parameters:
  • D (numpy arrays)

  • I (numpy arrays)

Returns:

XY

Return type:

array of equal area projections

Examples

>>> dimap_V([35,60,20],[70,80,-10])
array([[0.140856382055789, 0.20116376126988 ],
   [0.106743548942519, 0.061628416716219],
   [0.310909633795401, 0.85421719834377 ]])
sdcc.plotting.energy_matrix_plot(barriers)#

Plots the energy barriers between states as a matrix pair-plot.

Parameters:

barriers (numpy array) – Matrix of energy barriers

Return type:

None

sdcc.plotting.fast_path_to_min(energies, i, j)#

Fast ‘gradient descent’ like function for getting energy barrier paths.

sdcc.plotting.gradient_descent(max_iterations, threshold, xyz_init, k1, k2, rot_mat, LMN, Ms, ext_field=array([0, 0, 0]), learning_rate=0.0001)#

Slow gradient descent function, more accurate than fast one

sdcc.plotting.plot_arai(Zs_mag, Is_mag, B_anc, B_lab, temp_steps, ax=None)#

Plots an Arai plot, and expected Arai plot line for a paleointensity experiment.

Parameters:
  • Zs_mag (length n arrays) – Zero-field and in field data for Arai plot

  • Is_mag (length n arrays) – Zero-field and in field data for Arai plot

  • B_anc (floats) – Ancient and laboratory fields used in paleointensity simulation

  • B_lab (floats) – Ancient and laboratory fields used in paleointensity simulation

Return type:

None

sdcc.plotting.plot_barriers(barrier_thetas, barrier_phis, projection='equirectangular', ax=None)#

Plots a set of barrier locations on an energy surface plot as numerals.

Inputs#

barrier_thetas,barrier_phis: numpy arrays Locations of the barriers on the surface.

projection: str Either “equirectangular” or “stereonet”.

ax: matplotlib axes Axes to plot on top of.

sdcc.plotting.plot_cubeocta(ax, rot_mat=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), kind='cube', proj='upper')#

Plots the edges of a cube or an octahedron on the equal area plot, along the edges of the magnetocrystalline directions

Parameters:
  • ax (matplotlib axis)

  • to (axis to be plotted)

  • rot_mat (3x3 numpy array)

  • directions (rotation matrix for magnetocrystalline)

  • kind (string)

  • edges ("octa" for octahedron)

  • edges

  • proj (string)

  • hemisphere ("lower" for lower)

  • hemisphere

sdcc.plotting.plot_energy_path(TM, alignment, PRO, OBL, mask, T=20, ext_field=array([0, 0, 0]), n_perturbations=10, n_saddles=5, projection='equirectangular', method='fast', **kwargs)#

Plots the path for an energy surface by first finding energy barriers, then doing a gradient descent from nearby regions. Efficiency scales with n_perturbations * n_saddles. It’s not a fast function either way.

Parameters:
  • TMx (float) – Titanomagnetite composition (0 - 100)

  • alignment (str) – Alignment of magnetocrystalline and shape axis. Either ‘hard’ or ‘easy’ magnetocrystalline always aligned with shape easy.

  • PRO (float) – Prolateness of ellipsoid (major / intermediate axis)

  • OBL (float) – Oblateness of ellipsoid (intermediae / minor axis)

  • mask (numpy array) – An array of True/False wherever the energy barriers are -> this could probably be taken out and the function could just calculate where they are.

  • T (float) – Temperature (degrees C)

  • ext_field (numpy array) – Array of field_theta,field_phi,field_str where theta and phi are in radians and str is in Tesla.

  • n_perturbations (int) – number of perturbations to run from each energy barrier

  • n_saddles (int) – number of barriers to find paths from.

  • projection (str) – Either “equirectangular” or “stereonet”.

  • method (str) – Either fast or slow

  • kwargs – arguments to pass to plot_energy_surface function

Return type:

None

sdcc.plotting.plot_energy_surface(TMx, alignment, PRO, OBL, T=20, ext_field=Array([0, 0, 0], dtype=int64), levels=10, n_points=100, projection='equirectangular', cubic_dirs=False, ax=None)#

Plots an energy density surface for a single domain grain.

Parameters:
  • TMx (float) – Titanomagnetite composition (0 - 100)

  • alignment (str) – Alignment of magnetocrystalline and shape axis. Either ‘hard’ or ‘easy’ magnetocrystalline always aligned with shape easy.

  • PRO (float) – Prolateness of ellipsoid (major / intermediate axis)

  • OBL (float) – Oblateness of ellipsoid (intermediae / minor axis)

  • T (float) – Temperature (degrees C)

  • ext_field (numpy array) – Array of field_theta,field_phi,field_str where theta and phi are in radians and str is in Tesla.

  • levels (int) – Number of contour levels

  • n_points (int) – Number of grid points to evaluate the energy surface at.

  • projection (string) – Either “equirectangular” or “stereonet”.

Return type:

None

sdcc.plotting.plot_minima(minima_thetas, minima_phis, projection='equirectangular', ax=None)#

Plots a set of LEM states on an energy surface plot as numerals.

Inputs#

minima_thetas,minima_phis: numpy arrays Locations of the minima on the surface.

projection: str Either “equirectangular” or “stereonet”.

sdcc.plotting.plot_net(ax=None)#

Draws circle and tick marks for equal area projection. Adapted from pmagpy (PmagPy/PmagPy).

Copyright (c) 2023, PmagPy contributors All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

  • Redistributions of source code must retain the above copyright

notice, this list of conditions and the following disclaimer.

  • Redistributions in binary form must reproduce the above copyright

notice,this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

  • Neither the name of PmagPy nor the names of its contributors may

be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Parameters:
  • ax (matplotlib axis or None)

  • None) (axis to draw net onto (creates new axis if)

Return type:

None

sdcc.plotting.plot_pullaiah_curves(gel, ds, i, j, ax=None, plot_size=True, color='k', add_ticks=True, **kwargs)#

Plots Pullaiah curves for a particular energy barrier in a grain as a function of size.

Parameters:
  • gel (sdcc.barriers.gel object) – Energy landscape that can calculate barriers as a function of temperature.

  • ds (float or array of floats) – Equivalent sphere volume diameters of grains to plot.

  • i (int) – Index of state we are switching from in energy barrier.

  • j (int) – Index of state we are switching to in energy barrier.

  • ax (None or matplotlib axis) – axis to plot data to. If no axis specified, creates a new figure.

  • plot_size (bool) – If True, plots size next to pullaiah curve as text.

  • color (string or array of floats) – Matplotlib color for plotting

  • add_ticks – Whether to add y ticks to plot.

Return type:

None

sdcc.plotting.plot_relax_experiment(vs, relax_routine, ax=None)#

Plots a VRM decay experiment for a mono-dispersion of particles

Parameters:
  • vs (list) – Result from mono-dispersion

  • relax_routine (list) – Relaxation time routine

  • ax (matplotlib axis) – Axis to plot to

Return type:

None

sdcc.plotting.plot_routine(steps)#

Plots a set of thermal steps for a particular experiment or routine

Parameters:

steps (list of treatment.TreatmentStep objects)

Return type:

None

sdcc.plotting.update(xyz, k1, k2, rot_mat, LMN, Ms, ext_field, lr)#

Update for gradient descent function

sdcc.simulation module#

sdcc.simulation.Q_matrix(params: dict, d, field_dir=array([1, 0, 0]), field_str=0.0)#

Constructs a Q matrix of rates of transition between LEM states.

Parameters:
  • params (dictionary) – Dictionary output from a barriers.GEL object - contains states and energy barriers.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • field_dir (numpy array) – Unit vector with external field direction.

  • field_str (numpy array) – Strength of field (uT)

Returns:

Q – Array of transition rates.

Return type:

numpy array

class sdcc.simulation.SDCCResult(sizes, thermal_steps, vs, ps)#

Bases: object

Class to store results from a set of grain results - able to be dumped to file.

to_file(fname)#
sdcc.simulation._update_p_vector(p_vec, Q, dt)#

Given an initial state vector, a Q matrix and a time, calculates a new state vector. This is very slow due to the high floating point precision which is required and could probably benefit from a C++ implementation. Additionally - this is very susceptible to floating point errors even with the high precision when dt gets large. Using mpmath’s Pade approximations is slower than Taylor series and doesn’t seem to help much. If there’s an algorithm that improves this it would be extremely helpful as we’re dealing with some large numbers (age of Solar System) here.

Parameters:
  • p_vec (numpy array) – Vector of relative proportions of grains in each state.

  • Q (numpy array) – Rate matrix of transition times between states.

  • dt (float) – Amount of time spent in these field conditions/temperature.

Returns:

p_vec_new – New state vector after treatment applied.

Return type:

numpy array

sdcc.simulation.calc_relax_time(start_p, d, relax_routine, energy_landscape, ts)#

Function for calculating the relaxation time of a mono-dispersion of grains. The relaxation time is calculated as the time it takes for the magnetization to decay to 1/e.

Parameters:
  • start_p (numpy array) – Starting state vector

  • relax_routine (list of treatment.TreatmentStep objects.) – Steps describing a relaxation time treatment (cooling infield, followed by hold at room temperature infield).

  • energy_landscape (barriers.GEL object) – Object describing LEM states and energy barriers as a function of temperature.

  • ts (numpy array) – Array of time steps to check relaxation time at N.B. this should be roughly the same as relax_routine[1].ts - relax_routine[0].ts[-1]

sdcc.simulation.critical_size(K)#

Calculates the critical size (nm) given an energy barrier simply using the Neel relaxation time equation and nothing else.

sdcc.simulation.eq_ps(params, field_str, field_dir, d)#

Get the probabilities of each state in a grain under a specific set of conditions after an infinite amount of time.

Parameters:
  • params (dictionary) – Dictionary output from a barriers.GEL object - contains states and energy barriers.

  • field_str (float) – Field strength (uT)

  • field_dirs (numpy array) – Unit vector of field direction

  • T (float) – Temperature (degrees C)

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • Ms (float) – Saturation magnetization of grain (A/m)

Returns:

ps – Equilibrium state vector.

Return type:

numpy array

sdcc.simulation.full_crit_size(TMx, PRO, OBL, alignment)#

Calculate critical SD size of a grain, taking some shortcuts by using Neel relaxation time in cases when there should be only one energy barrier.

Parameters:
  • TMx (float) – Titanomagnetite composition % (0 - 100)

  • PRO (float) – Prolateness ratio (major axis/intermediate axis)

  • OBL (float) – Oblateness ratio (intermediate axis/minor axis)

  • alignment (string) – Either easy or hard. Specifies the magnetocrystalline direction that should be aligned with the x direction, which for our ellipsoids is the major (shape easy) axis.

Returns:

d – Critical SD size in nm.

Return type:

int

sdcc.simulation.get_avg_vectors(ps, theta_lists, phi_lists, Ts, rot_mat, energy_landscape, d)#

Obtains the average magnetization vectors for a grain during a thermal experiment, given the probabilities and magnetization directions of said state.

Parameters:
  • ps (numpy array) – Array of state vectors at each time step.

  • theta_lists (numpy arrays) – Arrays of magnetization directions at each time step

  • phi_lists (numpy arrays) – Arrays of magnetization directions at each time step

  • Ts (numpy array) – Array of temperatures at each time step

  • rot_mat (numpy array) – Rotation matrix applied to field direction - the inverse of this is applied to the states

  • energy_landscape (barriers.GEL object) – Object describing energy barriers and LEM states for a particular grain geometry. Here it’s used to get Ms.

Returns:

vs – Array of average magnetization directions at each time step

Return type:

numpy array

sdcc.simulation.grain_hyst_vectors(start_t, start_p, Bs, ts, d, energy_landscape: HEL, rot_mat, eq=False)#

Gets the state vectors and average magnetization vectors at each time step in a hysteresis experiment for a single direction in a mono-dispersion of grains. This calculation is performed for a single treatment step - i.e. a single heating or cooling. See treatment.TreatmentStep for a full description of this.

Parameters

start_p: numpy array

Initial state vector

Bs: numpy array

Set of field strengths at the times corresponding to ts.

ts: numpy array

Time steps at which we calculate the state.

d: float

Equivalent volume spherical diameter of grain (nm).

energy_landscape: barriers.HEL object

Object describing energy barriers for a particular grain geometry.

rot_mat: 3x3 matrix

Orientation associated with this grain.

field_dir: numpy array

Direction of field relative to grain - will be rotated to 1,0,0.

eq: bool

If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • vs (numpy array) – Array of average magnetization vectors at each time step

  • ps (numpy array) – Array of state vectors at each time step

sdcc.simulation.grain_vectors(start_t, start_p, Ts, ts, d, energy_landscape: GEL, rot_mat, field_strs, field_dirs, eq=False)#

Gets the state vectors and average magnetization vectors at each time step in a thermal treatment for a single direction in a mono-dispersion of grains. This calculation is performed for a single treatment step - i.e. a single heating or cooling. See treatment.TreatmentStep for a full description of this.

Parameters:
  • start_t (float) – Time at which this experiment step starts

  • start_p (numpy array) – Initial state vector

  • Ts (numpy array) – Set of temperatures at the times corresponding to ts.

  • ts (numpy array) – Time steps at which we calculate the state.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • energy_landscape (barriers.GEL object) – Object describing energy barriers for a particular grain geometry.

  • rot_mat (3x3 matrix) – Orientation of grain

  • field_strs (numpy array) – Array of field strengths at each time step.

  • field_dirs (numpy array) – Array of field directions at each time step.

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • vs (numpy array) – Array of average magnetization vectors at each time step

  • ps (numpy array) – Array of state vectors at each time step

sdcc.simulation.hyst_mono_dispersion(d, steps, energy_landscape, eq=False)#

Gets the state vectors and average magnetization vectors at each time step in a high-field treatment for all directions in a mono-dispersion of grains. This calculation is performed for a set of treatment steps - see treatment.TreatmentStep for more details.

Parameters:
  • start_p (numpy array) – Initial state vector of grain.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • steps (list of treatment.TreatmentStep objects) – Set of steps that describe a hysteresis experiment.

  • energy_landscape (barriers.HELs object) – Object describing LEM states and energy barriers as a function of fields.

  • n_dirs (int) – Number of Fibonacci sphere directions to use for mono-dispersion

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • vs (numpy array) – List of arrays of average magnetization vectors at each time step, in each treatment step, for each mono-dispersion direction.

  • ps (numpy array) – List of arrays of state vectors at each time step, in each treatment step, for each mono-dispersion direction.

sdcc.simulation.hyst_treatment(start_t, start_p, Bs, ts, d, energy_landscape: HEL, eq=False)#

Function for calculating the probability of different LEM states in a grain during a hysteresis experiment.

Parameters:
  • start_t (float) – Time at which this experiment step starts

  • start_p (numpy array) – Initial state vector

  • Bs (numpy array) – Set of field strengths at the times corresponding to ts.

  • ts (numpy array) – Time steps at which we calculate the state.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • energy_landscape (barriers.HEL object) – Object describing energy barriers for a particular grain geometry.

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • ps (numpy array) – Array of state vectors at each time step

  • theta_lists (numpy array) – Magnetization directions at each time step

  • phi_lists (numpy array) – Magnetization magnitudes at each time step

sdcc.simulation.mono_direction(rot_mat, start_p, d, steps, energy_landscape: GEL, eq=[False])#

Gets the state vectors and average magnetization vectors at each time step in a thermal treatment for a single direction in a mono-dispersion of grains. This calculation is performed for a set of treatment steps - see treatment.TreatmentStep for more details.

Parameters:
  • start_p (numpy array) – Initial state vector of grain.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • steps (list of treatment.TreatmentStep objects) – Set of steps that describe a thermal experiment.

  • energy_landscape (barriers.GEL object) – Object describing LEM states and energy barriers as a function of temperature

  • rot_mat (numpy array) – Direction of this grain in the mono dispersion.

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • vs (lists) – List of arrays of average magnetization vectors at each time step, in each treatment step.

  • ps (list) – List of arrays of state vectors at each time step, in each treatment step.

sdcc.simulation.mono_dispersion(start_p, d, steps, energy_landscape: GEL, n_dirs=50, eq=False)#

Gets the state vectors and average magnetization vectors at each time step in a thermal treatment for all directions in a mono-dispersion of grains. This calculation is performed for a set of treatment steps - see treatment.TreatmentStep for more details. N.B. - Recommend using parallelized_mono_dispersion instead of this, it’s a lot faster.

Parameters:
  • start_p (numpy array) – Initial state vector of grain.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • steps (list of treatment.TreatmentStep objects) – Set of steps that describe a thermal experiment.

  • energy_landscape (barriers.GEL object) – Object describing LEM states and energy barriers as a function of temperature

  • n_dirs (int) – Number of Fibonacci sphere directions to use for mono-dispersion

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • vs (numpy array) – List of arrays of average magnetization vectors at each time step, in each treatment step, for each mono-dispersion direction.

  • ps (numpy array) – List of arrays of state vectors at each time step, in each treatment step, for each mono-dispersion direction.

sdcc.simulation.mono_hyst_direction(start_p, d, steps, energy_landscape: HEL, eq=[False])#

Gets the state vectors and average magnetization vectors at each time step in a thermal treatment for a single direction in a mono-dispersion of grains. This calculation is performed for a set of treatment steps - see treatment.TreatmentStep for more details.

Parameters:
  • start_p (numpy array) – Initial state vector of grain.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • steps (list of treatment.TreatmentStep objects) – Set of steps that describe a thermal experiment.

  • energy_landscape (barriers.GEL object) – Object describing LEM states and energy barriers as a function of temperature

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • vs (lists) – List of arrays of average magnetization vectors at each time step, in each treatment step.

  • ps (list) – List of arrays of state vectors at each time step, in each treatment step.

sdcc.simulation.parallelized_hyst_mono_dispersion(start_p, d, steps, energy_landscape, eq=False, cpu_count=0, ctx=None)#

Gets the state vectors and average magnetization vectors at each time step in a hysteresis treatment for all directions in a mono-dispersion of grains. This calculation is performed for a set of treatment steps - see treatment.TreatmentStep for more details.

Parameters:
  • start_p (numpy array) – Initial state vector of grain.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • steps (list of treatment.TreatmentStep objects) – Set of steps that describe a thermal experiment.

  • energy_landscape (barriers.GEL object) – Object describing LEM states and energy barriers as a function of temperature

  • n_dirs (int) – Number of Fibonacci sphere directions to use for mono-dispersion

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

  • cpu_count (int) – Number of cores to parallelize with, if 0 chooses available cores

  • ctx (None or str) – context used for multiprocessing. If working with None, then don’t change this. Otherwise, set to “spawn” and call this function in a python script inside an if __name__ == “__main__” statement. Note: if set to spawn, multiprocessing will not work in a jupyter notebook

Returns:

  • vs (numpy array) – List of arrays of average magnetization vectors at each time step, in each treatment step, for each mono-dispersion direction.

  • ps (numpy array) – List of arrays of state vectors at each time step, in each treatment step, for each mono-dispersion direction.

sdcc.simulation.parallelized_mono_dispersion(start_p, d, steps, energy_landscape: GEL, n_dirs=50, eq=False, cpu_count=0, ctx=None)#

Gets the state vectors and average magnetization vectors at each time step in a thermal treatment for all directions in a mono-dispersion of grains. This calculation is performed for a set of treatment steps - see treatment.TreatmentStep for more details.

Parameters:
  • start_p (numpy array) – Initial state vector of grain.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • steps (list of treatment.TreatmentStep objects) – Set of steps that describe a thermal experiment.

  • energy_landscape (barriers.GEL object) – Object describing LEM states and energy barriers as a function of temperature

  • n_dirs (int) – Number of Fibonacci sphere directions to use for mono-dispersion

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

  • cpu_count (int) – Number of cores to parallelize with, if 0 chooses available cores

  • ctx (None or str) – context used for multiprocessing. If working with None, then don’t change this. Otherwise, set to “spawn” and call this function in a python script inside an if __name__ == “__main__” statement. Note: if set to spawn, multiprocessing will not work in a jupyter notebook

Returns:

  • vs (numpy array) – List of arrays of average magnetization vectors at each time step, in each treatment step, for each mono-dispersion direction.

  • ps (numpy array) – List of arrays of state vectors at each time step, in each treatment step, for each mono-dispersion direction.

sdcc.simulation.relax_time_crit_size(relax_routine, energy_landscape, init_size=[5], size_incr=10)#

Finds the critical SP size of a grain.

Parameters:
  • relax_routine (list of treatment.TreatmentStep objects.) – Steps describing a relaxation time treatment (cooling infield, followed by hold at room temperature infield).

  • energy_landscape (barriers.GEL object) – Object describing LEM states and energy barriers as a function of temperature.

  • init_size (list of ints) – Initial grain sizes to try.

  • size_incr (int) – Amount to increment size by in certain situations.

Returns:

d – Critical SD size in nm.

Return type:

int

sdcc.simulation.result_to_file(energyLandscape, size, routine, moments, probabilities, file_ext: str, directory='./')#

Saves an SDCC result to file, along with information about the result

Parameters:
  • energyLandscape (GEL, HEL or HELs object) – The object that describes the particle energies

  • size (float) – The size of the particle (in nm)

  • routine (list)

  • steps (The routine of treatment)

  • moments (list) – The moments output by the SDCC

  • probabilities (list) – The probabilities of each state in each particle at each time-step.

  • file_ext (str) – The file extension - should relate to the experiment type.

Return type:

None

sdcc.simulation.thermal_treatment(start_t, start_p, Ts, ts, d, energy_landscape: GEL, field_strs, field_dirs, eq=False)#

Function for calculating the probability of different LEM states in a grain during a thermal experiment.

Parameters:
  • start_t (float) – Time at which this experiment step starts

  • start_p (numpy array) – Initial state vector

  • Ts (numpy array) – Set of temperatures at the times corresponding to ts.

  • ts (numpy array) – Time steps at which we calculate the state.

  • d (float) – Equivalent volume spherical diameter of grain (nm).

  • energy_landscape (barriers.GEL object) – Object describing energy barriers for a particular grain geometry.

  • field_strs (numpy array) – Array of field strengths at each time step.

  • field_dirs (numpy array) – Array of field directions at each time step.

  • eq (bool) – If True, ignore time steps and run magnetization to equilibrium.

Returns:

  • ps (numpy array) – Array of state vectors at each time step

  • theta_lists (numpy array) – Magnetization directions at each time step

  • phi_lists (numpy array) – Magnetization magnitudes at each time step

sdcc.treatment module#

class sdcc.treatment.CoolingStep(t_start, T_start, T_end, field_str, field_dir, char_time=1, max_temp=None, char_temp=None)#

Bases: TreatmentStep

Treatment Step for cooling a specimen in a known field. Assumes Newtonian cooling.

Parameters:
  • t_start (float) – Start time of step (in seconds)

  • T_start (float) – Start temperature for step (degrees C)

  • T_end (float) – Ambient temperature cooled to

  • field_str (float) – Magnetic field strength during step (in µT)

  • field_dir (array) – Magnetic field direction during step (cartesian unit vector)

  • char_time (float) – Characteristic cooling time.

  • max_temp (floats) – Cooling will take char_time to cool from max_temp to char_temp

  • char_temp (floats) – Cooling will take char_time to cool from max_temp to char_temp

class sdcc.treatment.HeatingStep(t_start, T_start, T_end, field_str, field_dir, lin_rate=0.3333333333333333)#

Bases: TreatmentStep

Treatment step for heating a specimen in a known field. Assumes linear heating ramp.

Parameters:
  • t_start (float) – Start time for step (s)

  • T_start (float) – Start temperature for ramp (degrees C)

  • T_end (float) – Peak temperature of ramp (degrees C)

  • field_str (float) – Strength of field (in µT) during step

  • field_dir (array) – Direction of field (Cartesian unit vector) during step

  • lin_rate (float) – Rate (degrees / s) of ramp.

class sdcc.treatment.HoldStep(t_start, T_start, field_str, field_dir, hold_steps=2, hold_time=1800.0)#

Bases: TreatmentStep

Treatment step for holding a specimen at constant field and temperature.

Parameters:
  • t_start (float) – start time of step (s)

  • T_start (float) – hold temperature (degrees C)

  • field_str (float) – strength of field during step (µT)

  • field_dir (array) – direction of field during step (cartesian unit vector)

  • hold_steps (int) – number of steps to break hold time into, unless you want measurements at a particular time, leave at 2.

  • hold_time (float) – length of hold time (in seconds)

class sdcc.treatment.HystBranch(t_start, B_start, B_end, B_step=1000.0, t_step=1, T=20, B_dir=[1, 0, 0])#

Bases: TreatmentStep

Treatment step for hysteresis loop branch.

Parameters:
  • t_start (float) – start time of step (s)

  • B_start (float) – start field of step (µT)

  • B_end (float) – end field of step (µT)

  • B_step (float) – field step at which measurements are made (µT)

  • t_step (float) – length of time to ramp field by B_step (s)

  • T (float) – temperature at which hysteresis loop measured (degrees C)

  • B_dir (array) – direction of field (cartesian unit vector)

class sdcc.treatment.SuscStep(t_start, t_dur, freq, field_peak, field_dir, n_points=100, T=20.0)#

Bases: TreatmentStep

Treatment step for susceptibility measurements. Field changes rather than temperature. Experimental and liable to change!

Parameters:
  • t_start (float) – start time of step (s)

  • t_dur (float) – length of susceptibility treatment (s)

  • freq (float) – frequency of AC susceptibility (hz)

  • field_peak (float) – maximum field amplitude (µT)

  • field_dir (array) – field direction (cartesian unit vector)

  • n_points (int) – number of points to evaluate susceptibility at

  • T (float) – temperature susceptibility experiment conducted at

class sdcc.treatment.TreatmentStep(ts, Ts, field_strs, field_dirs)#

Bases: object

Class for a thermal treatment step. Instead of having a lot of cases, I think it would make more sense to make each type of treatment step a separate class that inherits from Step - more flexible to create new stuff this way.

class sdcc.treatment.VRMStep(t_start, T_start, field_str, field_dir, hold_steps=100, hold_time=1e+17)#

Bases: TreatmentStep

Treatment step for holding a specimen in a field at constant temperature with logarithmically spaced time steps (meant for VRM acquisitions).

Parameters:
  • t_start (float) – start time of step (s)

  • T_start (float) – temperature of VRM acquisition (degrees C)

  • field_str (float) – strength of external field (µT)

  • field_dir (array) – direction of external field (cartesian unit vector)

  • hold_steps (int) – number of hold steps

  • hold_time (float) – time for VRM acquisition

sdcc.treatment.coe_experiment(temp_steps, B_anc, B_lab, B_ancdir, B_labdir)#

Creates a set of thermal treatment steps for an IZZI-Coe experiment

Parameters:
  • temp_steps (numpy array) – Set of temperatures for the Coe experiment

  • B_anc (float) – Ancient field strength (T)

  • B_lab (float) – Lab field strength (T)

  • B_ancdir (numpy array) – Unit vector ancient field direction

  • B_labdir (numpy array) – Unit vector lab field direction

Returns:

steps – Set of steps for coe experiment.

Return type:

list of treatment.ThermalStep objects

sdcc.treatment.hyst_loop(B_max, B_step=0.001, **kwargs)#

Creates a set of thermal treatment steps for a hysteresis loop

Parameters:
  • B_max (float) – Maximum field (in T)

  • B_min (float) – Field step value (in T)

  • **kwargs (keyword argument dictionary) – keyword arguments that are passed to sdcc.treatment.HystBranch

Returns:

steps – Set of steps for hysteresis loop

Return type:

list of treatment.HystBranch objects

sdcc.treatment.make_thellier_step(steps, T, T_max, T_min, B_str, B_dir)#

Adds a sequence of steps in a Thellier experiment (cooling, hold, heating).

Parameters:
  • steps (list of treatment.TreatmentStep objects)

  • far (steps in the experiment so)

  • T (float) – peak temperature of step

  • T_max (float) – T_max parameter of original NRM step

  • T_min (float) – T_min parameter of original NRM step

  • B_str (float) – Field strength (microT)

  • B_dir (numpy array) – Cartesian unit vector of field direction

Returns:

  • step (list of treatment.TreatmentStep objects)

  • List containing steps for thellier experiment.

sdcc.treatment.overprint(TRM_mag, TRM_dir, oPrint_mag, oPrint_dir, oPrint_T, oPrint_t, gel)#

Creates a set of thermal treatment steps for a TRM followed by a VRM overprint at room temperature.

Parameters:
  • TRM_mag (float) – magnitude of TRM field (µT)

  • TRM_dir (array) – direction of TRM field (cartesian unit vector)

  • oPrint_mag (float) – magnitude of overprinting field (µT)

  • oprint_dir (array) – direction of overprinting field (cartesian unit vector)

  • oPrint_T (float) – temperature overprint acquired at.

  • oPrint_t (float) – time over which overprint acquired

  • gel (SDCC.barriers.GEL object) – Energy landscape of particle.

Returns:

steps – steps for TRM and VRM acquisitions

Return type:

list of SDCC.treatment.ThermalStep objects

sdcc.treatment.relaxation_time(energy_landscape: GEL, B_dir, B)#

Creates a set of thermal treatment steps for a relaxation time estimate.

Parameters:
  • energy_landscape (barriers.GEL object) – Object describing energy barriers and LEM states for a particular grain geometry. Here it’s used to get Ms.

  • B_dir (numpy array) – Unit vector field direction

  • B (float) – Field strength (uT)

Returns:

steps – Set of steps for coe experiment.

Return type:

list of treatment.ThermalStep objects

sdcc.treatment.temp2time(T, t1, T0, T1, T_amb)#

Converts a temperature to a time assuming Newtonian cooling.

Parameters:
  • T (numpy array) – Array of temperature steps

  • t1 (float) – Characteristic cooling time (i.e. time to T1)

  • T0 (float) – Maximum temperature.

  • T1 (float) – Characteristic temperature at time T1

  • T_amb (float) – Ambient temperature (usually 20C)

Returns:

t – Array of times at temperatures T.

Return type:

numpy array

sdcc.treatment.thellier_experiment(temp_steps, B_anc, B_lab, B_ancdir, B_labdir, type='coe', ptrm_checks=0, **kwargs)#

Creates a set of thermal treatment steps for an IZZI-Coe experiment

Parameters:
  • temp_steps (numpy array) – Set of temperatures for the Coe experiment

  • B_anc (float) – Ancient field strength (T)

  • B_lab (float) – Lab field strength (T)

  • B_ancdir (numpy array) – Unit vector ancient field direction

  • B_labdir (numpy array) – Unit vector lab field direction

  • type (string) – Either “coe”, “aitken” or “izzi”. Step order for experiment

  • ptrm_checks (int) – If 0, no pTRM checks. If > 0, include pTRM checks (Coe and Aitken only). If 1, perform pTRM check between Z and I in ZI step (IZZI only). If 2, perform pTRM check after Z in IZ step (IZZI only).

Returns:

steps – Set of steps for coe experiment.

Return type:

list of treatment.ThermalStep objects

sdcc.treatment.thermal_demag(Ts)#

Creates a set of thermal treatment steps for a simple thermal demagnetization

Parameters:

Ts (array of floats) – temperatures of thermal demagnetization steps (degrees C)

Returns:

steps – steps for thermal demagnetization experient

Return type:

list of SDCC.treatment.ThermalStep objects

sdcc.treatment.time2temp(t, t1, T0, T1, T_amb)#

Converts a time to a temperature assuming Newtonian cooling.

Parameters:
  • t (numpy array) – Array of time steps

  • t1 (float) – Characteristic cooling time (i.e. time to T1)

  • T0 (float) – Maximum temperature.

  • T1 (float) – Characteristic temperature at time T1

  • T_amb (float) – Ambient temperature (usually 20C)

Returns:

T – Array of temperatures at times ts.

Return type:

numpy array

sdcc.utils module#

sdcc.utils.calc_d_min(TMx, alignment, PRO, OBL)#

Calculates the SD size limit (d_min) using data from Cych et al (Magnetic Domain States and Critical Sizes in the Titanomagnetite Series, 2024)

Parameters:
  • TMx (float) – Titanomagnetite Composition

  • alignment (str) – Magnetocrystalline axis aligned with shape easy axis, either ‘easy’,’intermediate’ or ‘hard’

  • PRO (float) – Prolateness (> 1)

  • OBL (float) – Oblateness (> 1)

Returns:

d_min – Critical size below which the SV state does not exist.

Return type:

float

sdcc.utils.fib_hypersphere(n)#

Algorithm for producing uniformly distributed orientations from a super Fibonacci spiral according to the method of Alexa (2022). This is required because the fibonacci sphere does not describe a uniformly destributed set of rotation matrices, but instead describes a uniformly distributed set of directions. This is important when rotation a magnetization which is not in the same orientation as the field direction, as those vectors may not be uniformly oriented.

Parameters:

n (int) – Number of directions

Returns:

rot_mats – Array of 3x3 rotation matrices of length n describing a set of uniformly distributed rotations.

Return type:

n x 3 x 3 array

sdcc.utils.fib_sphere(n=1000)#

Algorithm for producing directions from a Fibonacci sphere

Parameters:

n (int) – Number of directions

Returns:

xyz – Array of 3D cartesian directions

Return type:

array

sdcc.utils.smooth_surface(xx, yy, zs, grid_x, grid_y, kind='cubic')#

Interpolates critical size values on a grid to obtain a smoothed surface.

Parameters:
  • xx (numpy array) – Meshgrid of TM compositions

  • yy (numpy array) – Meshgrid of shapes

  • zs (numpy array) – Critical sizes at (xx,yy)

  • grid_x (float or array) – TMx value(s) to calculate smoothed surface for.

  • grid_y (float or array) – Shapes to calculate smoothed surface for.

  • kind (string) – What kind of smoothing to do (should always be ‘cubic’)

Returns:

grid_z – Interpolated critical size(s)

Return type:

numpy array

Module contents#