Functions¶

Fast Ion Distributions¶

Bi-Maxwellian Distribution¶

functions.biMaxx(vpara, vperp, varargin)¶

Evaluates the bi-Maxwellian fast-ion velocity distribution on a (vpara,vperp)-grid.

Usage:

x = biMaxx(vpara,vperp)

[x, info] = biMaxx(vpara,vperp,varargin)

Inputs:

vpara: Meshgrid for X to be evaluated on.
vperp: Meshgrid for X to be evaluated on.

Optional inputs:

Tpara: Parallel ion temperature. Default : 20 keV
Tperp: Perpendicular ion temperature. Default : 20 keV
vparadrift: Parallel-drift velocity. Default : 5e5 m/s
Mi: Mass of ion. Default : 2*Mp
ne: Number of ions. Default : 1e19

Output:

x: The Tikhonov solution for a given alpha.
info: MATLAB struct containing all inputs (including optional inputs) above.

Isotropic Slowing-Down Distribution¶

functions.isoSDx(vpara, vperp, varargin)¶

Evaluates the isotropic slowing-down distribution on a (vpara,vperp)-grid.

Usage:

x = isoSDx(vpara,vperp)

[x, info] = isoSDx(vpara,vperp,varargin)

Inputs:

vpara: Meshgrid for X to be evaluated on.
vperp: Meshgrid for X to be evaluated on.

Optional inputs:

Ecrit: Critical energy. Default : 44*20000 eV
Ebirth: Birth energy. Default : 3.5e6 eV
Ebirthwidth: Birth energy width. Default : 6e4
Mi: Mass of ions. Default : 4*Mp
ne: Number of ions. Default : 1e19

Output:

x: The isotropic slowing-down distribution evaluated on the (vpara,vperp)-grid.
info: MATLAB struct containing all inputs (including optional inputs) above.

Analytic Velocity Distributions¶

Bi-Maxwellian Distribution¶

functions.biMaxg(u, phi, varargin)¶

Evaluates the analytic projected velocity distribution of the bi-Maxwellian velocity distribution.

Usage:

g = biMaxg(u)

g = biMaxg(u,phi)

[g, info] = biMaxg(u,phi,varargin)

Inputs:

u: Vector of projected velocities.
phi: Observation angles.

Optional inputs:

Tpara: Parallel ion temperature. Default : 20 keV
Tperp: Perpendicular ion temperature. Default : 20 keV
vparadrift: Parallel-drift velocity. Default : 5e5
Mi: Mass of ion. Default : 2*Mp
ne: Number of ions. Default : 1e19

Output:

g: Analytic projected velocity distribution.
info: MATLAB struct containing all inputs (including optional inputs) above.

Isotropic Slowing-Down Distribution¶

functions.isoSDg(u, phi, varargin)¶

Evaluates the analytic projected velocity distribution of the slowing down fast-ion velocity distribution.

Usage:

g = isoSDg(u)

g = isoSDg(u,phi)

[g, info] = isoSDg(u,phi,varargin)

Inputs:

u: Vector of projected velocities.
phi: Observation angles.

Optional inputs:

Ecrit: Critical energy. Default : 44*20000 eV
Ebirth: Birth energy. Default : 3.5e6 eV
Ebirthwidth: Birth energy width. Default : 6e4
Mi: Mass of ions. Default : 4*Mp
ne: Number of ions. Default : 1e19

Output:

g: Projected velocity distribution.
info: MATLAB struct containing all inputs (including optional inputs) above.

Forward Model¶

Transfer Matrix¶

functions.transferMatrix(vpara, vperp, phi, u)¶

Constructs the forward-model matrix A.

Usage:

A = transferMatrix(vpara, vperp, phi, u)

Inputs:

vpara: vpara-grid.
vperp: vperp-grid.
phi: Vector of observation angles.
u: Vector of projected velocity values.

Output:

A: The forward model matrix A.

Solvers¶

Nonnegative Tikhonov Solver¶

functions.TikhNN(A, b, alpha, L, varargin)¶

Solves the Tikhonov problem formulated as:

\(\mathrm{x}_{\alpha} = \underset{x}{\min} \frac{1}{2} \left\| \mathrm{A} \mathrm{x} - \mathrm{b} \right\|^2 + \frac{\alpha}{2} \|\mathrm{L}\mathrm{x}\|^2\)

or equivalently

\(\mathrm{x}_{\alpha} = \underset{x}{\min} \frac{1}{2} \left\| \begin{bmatrix} \mathrm{A} \\ \sqrt{\alpha} \, \mathrm{L} \end{bmatrix} \mathrm{x} - \begin{bmatrix} \mathrm{b} \\ 0 \end{bmatrix} \right\|^2\)

Usage:

[x, alpha] = TikhNN(A,b,alpha)

[x, alpha] = TikhNN(A,b,alpha,L)

[x, alpha, relerr] = TikhNN(A,b,alpha,L,'return_relerr',true,'xtrue',xtrue)

Inputs:

A: System matrix
b: Right hand side
alpha: Regularization parameter
L: Regularization matrix

Optional inputs:

solver: String specifying the used solver. Either lsqnonneg (default) , GPCG or \.
scaling: Whether or not A should be scaled.
dispwaitbar: Whether or not to display waitbar if several alphas are passed.
return_relerr: If true, then relerr is returned as the third output.
xtrue: Used to calculate the relative error.

Output:

x: The Tikhonov solution for a given alpha.
alpha: Corresponding value of alpha.
relerr: Relative error

Samplers¶

Nonnegative Hierachical Gibbs Sampler¶

functions.NNHGS(A, b, L, n, varargin)¶

Nonnegative Hierachical Gibbs Sampler

Usage:

[x, alpha] = NNHGS(A,b,L,n)

[x, alpha, delta, lambda, info] = NNHGS(A,b,L,n,varargin)

Inputs:

A: The system matrix A.
b: The right hand side.
L: The regularization matrix. [] is interpreted as speye(N)
n: Number of samples to be generated.

Optional inputs:

disp_waitbar: Whether or not a waitbar should be displayed.
welfordsalgorithm: If true, then Welford’s Online Algorithm is used to return mean and standard deviation of the samples.
nburnin: Samples to be generated before Welford’s Algorithm starts. Only used if welford=true.
keep: Which samples we should keep. Only used if welford=true. Default value is 1, meaning that every sample is kept.
solver: Which solver is used. lsqnonneg is default, but GPCG or \ can be set.
scaling: Whether or not A should be rescaled such that the largest element in A = 1. Default is true.

Output:

x: Samples drawn from the posterior
alpha: alpha-chain sampled from the posterior.
delta: delta-chain sampled from the posterior.
lambda: lambda-chain sampled from the posterior.
info: MATLAB struct containing all sorts of information.

Auxilary functions¶

functions.construct_vgrid(vparadim, vperpdim, varargin)¶

Function for constructing (vpara,vperp)-grid on which we can evaluate our distribution functions.

Usage:

[vpara, vperp] = construct_vgrid(vparadim, vperpdim)

[vpara, vperp, ginfo] = construct_vgrid(vparadim, vperpdim, varargin)

Inputs:

vparadim: The dimensions in the parallel dimension.
vperpdim: The dimensions in the perpendicular dimension.

Optional inputs:

vparamin: Minimum parallel velocity.
vparamax: Maximum parallel velocity.
vperpmin: Minimum perpendicular velocity.
vperpmax: Maximum perpendicular velocity.

Output:

vpara: A vpara-grid.
vperp: A vperp-grid.
ginfo: Information about the grid, used for plotting.

functions.construct_uvec(varargin)¶

Construct the u vector.

Usage:

u = construct_uvec()

u = construct_uvec(varargin)

Optional inputs:

Emin: Minimum projected energy.
Emax: Maximum projected energy.
umin: Minimum projected velocity u.
umax: Maximum projected velocity u.
du: Difference between elements in the vector.
udim: Dimension of the output vector.

Output:

u: u-vector.

functions.reguL(varargin)¶

Constructs the 1st order Tikhonov Regularization matrix without any assumptions about the border.

Usage:

L = reguL(vpara, vperp)

Inputs:

vpara: vpara-grid
vperp: vperp-grid

Output:

L: The regularization matrix L.

functions.showDistribution(X, gridinfo, caxis_limits)¶

Displays a fast-ion velocity distribution.

Usage:

showDistribution(X, ginfo)

showDistribution(X, ginfo, caxis_limits)

Inputs:

X: Velocity distribution X as a vector.
gridinfo: gridinfo structure.
caxis_limits: Color-axis limits.

functions.generate_noisy_b(A, x, varargin)¶

Generates a noisy right hand side.

Usage:

[b, e] = generate_noisy_b(A, x)

[b, e] = generate_noisy_b(A, x, varargin)

Inputs:

A: System matrix
x: True solution

Optional inputs:

noise_level: Noise level
background_level: Background level for the noise.

Output:

b: Noisy right hand side.
e: Noise used for normalization.

functions.cbounds(xsamp)¶

Calculates the width of 95% credibility bounds.

Usage:

c = cbounds(xsamp)

Inputs:

xsamp: An N x n matrix, where each column is a sample.

Output:

c: The width of the credibility bounds.

functions.check_mosek()¶

A wrapper for mosekdiag, checking if mosek works correctly.

Usage:: check_mosek()
Output:: 1 if working correctly, 0 otherwise.

functions.geweketest(v, burnin, trim)¶

Calculates the Geweke test statistic. Here, 0 means that the chain is stationary and 1 means that the chain is not.

Usage:

p = geweketest(v)

[p, pall] = geweketest(v, burnin, trim)

Inputs:

v: Matrix where each column is a chain, e.g. [alpha lambda delta].
burnin: Removes first burnin rows, default is 0.
trim: Keep every trim elements. Default is 1.

Output:

p: Geweke p-score for each chain.
pall: Geweke p-score for all chains.