# An introduction to MEQ

In this example we will use the 'anamak' tokamak, a tokamak with a simplified geometry description for educational and code testing purposes.

## Contents

- Codes In MEQ
- The MEQ structures
- The `P` structure
- The `G` structure
- Fluxes and fields
- Coordinate system
- `meq` data
- My first meq equilibrium

## Codes In MEQ

- FBTE (`fbt`) - Fit desired equilibrium description
- LIUQE (`liu`) - Equilibrium reconstruction from measurements
- FGS (`fgs`) - Forward Grad-Shafranov Static
- FGE (`fge`) - Forward Grad-Shafranov Evolutive

Each code is treated in more detail in later tutorials.

## The MEQ structures

MEQ uses several standardized structures that contain information. It is good to be familiar with these structures since they appear in various places in the code.

If we call any code of the `meq` suite with only one output argument, the code is not run but only the so-called `L` structure is calculated. This structure contains parameters, geometry, and pre-calculated data used to run the code. Each code produces a different `L` structure although some fields will be shared.

Let's first get an `L` structure by calling `fbt`. The first parameter is the tokamak, where we choose `'ana'` for the anamak tokamak. The second argument represents the shot number, which will be discussed later.

    L = fbt('ana',1);

We see that this structure contains two substructures, `P` and `G` as well as other quantities. Let's first look at those.

## The `P` structure

    L.P;

This contains parameters that will be used to run the code. These parameters are documented in: * `meqp.m` for general parameters * `meqp<tok>.m` for tokamak-specific parameters (e.g. `meqpana.m`) * `fbtp.m` for code-specific parameters general to all tokamaks * `fbtp<tok>.m` for code-specific and tokamak-specific parameters, for example `fbtpana.m`

    help meqp
    help meqpana
    help fbtp
    help fbtpana

     MEQP  MEQ default algorithm configuration parameters
      P = MEQP defines default algorithm configuration parameters common to
      MEQ.

      .tokamak      Full name of tokamak
      .tok          Short name of tokamak
      .shot         #shotnum
      .dasm         Minimal distance (grid points) between X and axis points [asxy]
      .lasy         Search axes and X-points also outside the limiter contour
      .xdoma        Angle of X-point polygon edges with perpendicular (>0 values to get a smaller polygon) [pdom]
      .isaddl       If ~isaddl, ignore X-points when searching LCFS
      .ilim         If 0, ignore limiter when searching LCFS. If 1, use limiter points, if 2, use quadratic interpolation, if 3: cubic splines (requires icsint)
      .idoublet     If idoublet, allow doublet equilibria
      .ihole        If ihole (default), allow axes with current densities opposite to Ip when multiple axes are found
      .r0           default major radius used for outputs/normalizations
      .b0           default toroidal field used for outputs/normalizations
      .pq           If not empty provides user-defined radial grid (excl LCFS and axis)
      .pqu          maximum rho-psi for radial grid                   [active only if pq is empty]
      .pql          minimum rho-psi for radial grid                   [active only if pq is empty]
      .npq          number of radial grid points (excl LCFS and axis) [active only if pq is empty]
      .noq          size of theta-grid mesh for contours
      .raS          r/a values of flux surfaces for which to return r,z, coordinates
      .ri       [m] Inner boundary of computational grid                       [active only if selx is empty for TCV]
      .ro       [m] Outer boundary of computational grid                       [active only if selx is empty for TCV]
      .zl       [m] Lower boundary of computational grid                       [active only if selx is empty for TCV]
      .zu       [m] Upper boundary of computational grid                       [active only if selx is empty for TCV]
      .nr           Set radial grid points, true grid will be of size nr+1     [active only if selx is empty for TCV]
      .nz           Set vertical grid points, true grid will be of size 2*nz+1 [active only if selx is empty for TCV]
      .selu         Generalized vessel description selector string e.g. 's','e','v';
      .nu           Number of vessel variables (eigenmodes, segments)

      .cde          CDE used (also includes static CDEs for fgs), empty for no
                    cde, otherwise string or function handle for cde type

      .ilackner     flag for boundary flux computation in meqfbp: 0: use Mby*Iy 1: use Lackners trick, 2: Use a precomputed Mby representing Lackners trick

      .nFW          Number of fluxes to track for wall gaps. >1 gives secondary separatrices
      .rn       [m] Radial coordinate for points for flux, field interpolations
      .zn       [m] Z coordinate for points for flux, field interpolations
      .bfct         Basis function handle
      .bfp          Basis function parameters

      .infct        Interpolation function handle
      .inp          Interpolation function parameters

      .interpmethod Interpolation method for data (default: nearest)
      .iterq        Iterations of contour algorithm, enables contour-dependent calculations
      .tolq         Tolerance for convergence of contour algorithm, normalized to the radial grid size
      .iqR          Value of 1/q for which to return r/a location (only if iterq>0)
      .naR          Number of occurrences of each 1/q surfaces to track (size of raR will be [nR,naR])
      .ifield       Calculate Br,Bz fields as additional outputs
      .ivacuum      Calculate vacuum fluxes and fields Br0x, Bz0x, F0x (excluding plasma contribution)
      .izgrid       Additional 'z' grid for post-processing flux, field calculations
      .izxoverlap   Force partial overlapping of z and x grids, allows re-use of Fx in Fz
      .relax        Relaxation parameter for SQP and Newton methods

      .Ipmin        Minimum plasma current for which GS solution will be attempted

      .icsint       If icsint, use cubic spline interpolation for domain id. and contouring (default: false)
      .tolcs        Tolerance on position of axes, saddle points, limiter points                        [active only if icsint is true]
      .itercs       Maximum number of Newton-Raphson iterations for asxycs/flcs                         [active only if icsint is true]
      .nthreadscs   Maximum number of threads to be used in bspsum MEX-file for cubic spline evaluation [active only if icsint is true]

      .mdsserver    MDS server hostname for parameters/measurements loading

      .gsxe         Switch for treatment of external currents inside computational grid
      .gsxetol      Tolerance for mapping of external currents inside computational grid to y grid (when gsxe=2)
      .gfile        If non-empty, attempt to load geometry data from this file (see meqgload/meqgsave)
      .LYall        If true, LY contains all time slices, converged or not

      .ohcoils      Names of coils used for driving Ohmic current (from G.dima)
      .compcoils    Names of coils used for compensation of stray fields from Ohmic transformer coils
      .ohequal     Force currents of coils for OH induction to be equal

      .debug        debug flag
      .debugplot    level of debug plotting
      .debugfail    trigger extra debugging upon encontering failed slices/equilibria
      .argout       Additional output arguments (syntax: '[file:]fld[=exp]')

      [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

      Default parameters for Anamak - an analytical test tokamak

      Anamak vessel, coils, probes, flux loops are distributed along
      concentric elliptical shapes.

      .ac       PF coil minor radius
      .af       Flux loop minor radius
      .al       Limiter minor radius
      .am       Magnetic probe minor radius
      .av       Vessel minor radius
      .cappav   Vessel elongation
      .nc       Number of PF coils
      .nf       Number of flux loops
      .nl       Number of limiter points
      .nm       Number of magnetic probes
      .nv       Number of vessel elements
      .nW       Number of wall gaps
      .wc       Coil width
      .fastcoil Internal fast coil (true/false)

      [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

     FBTP  FBTE default configuration parameters
      P = FBTP(P,'PAR',VAL,...) fills the structure P with default
      configuration parameters, optionally replacing or adding parameters with
      specified values.

      FBT-specific parameters:
      .fbtagcon    Internal profile constraints (see meqagcon for details)
      .limm        Margin factor for coil current limits, values lower than 1 are more conservative. Can be a scalar or a column vector with as many elements as limits.
      .mergex      If true, merge redundant X-point constraints in LX
      .circuit     If true, solve time-dependent problem with circuit equations (GSpulse algorithm)
      .itert       If not 0, re-use previous Iy from previous time as starting guess
      .iterfrz     Number of iterations before which plasma domain is frozen to ease convergence
      .interpx     Interpolation method for LX time base. Options: 'none', 'nearest','linear','previous','next'
      .MHyinterp   If true, use qintmex interpolation of Mxy from meqFxjac for computing control point sensitivity MHy
      .qpsolver    Solver for contrained QP problem. 'quadprog' (default if available), or 'ipm-optimized' (fallback) or 'ipm' (see HELP IPM).
      .useSQP      If true, use sequential quadratic programming
      .wsSQP       Number of picard iterations to perform to warm-start the SQP iterations
      .initineq    If true, inequality constraints will be included to compute the initial guess

      .algoNL      Plasma convergence method
      .xscal       If true, scale internal plasma state as in FGS
      .tolstep     Tolerance for step size
      .timeder     Type of stencil for time derivatives, 'centered' or 'backward'
      .jacobian_handle    Use function handle form of jacobian
      .jacobian_sparse    Use sparse-matrix form of jacobian (algoNL=Newton-GS)

      SQP Solver parameters
      .userowmask  Reduce Newton problem size when possible
      .algoGMRES   Algorithm variant for Block-GMRES method (direct_inversion,sim,giv,qr)
      .mkryl       Maximum krylov subspace dimension for Block-GMRES
      .epsilon_res Residual tolerance for Block-GMRES
      .usepreconditioner  If true, use Picard approx. preconditioner for Block-GMRES

      Quantities which will end up in LX
      .tol         Tolerance for residual
      .capaj       High capaj means stronger peaking of initial plasma current density guess
      .niter       Maximum number of iterations

      if bffbt basis function is used, additional parameters exist which are documented in BFPFBT

      [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

     FBTPANA  ANA FBT configuration parameters
      P = FBTPANA(SHOT,'PAR',VAL,...) returns a structure P with configuration
      parameters for ANA tokamak, optionally replacing or adding
      parameters with specified values. See also FBTP.

      [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

These parameters can be customized by entering parameter, value pairs in the call to the code (after the first three arguments being `tokamak`, `shot` and `time`), for example

    fprintf('iterq=%d\n',L.P.iterq)  % old L.P.iterq value
    L = fbt('ana',1,[],'iterq',1); % set another value
    fprintf('iterq=%d\n',L.P.iterq) % new L.P.iterq value

    iterq=0
    iterq=1

## The `G` structure

The `G` structure contains the geometry description of the tokamak, including many useful quantities such as mutual inductances, grids, coil locations and more.

    L.G;

We can plot the geometry using `meqgplot.m`

    clf; meqgplot(L.G)

Or we can plot quantites we want manually by hand:

    figure; clf; hold on;
    plot(L.G.rw,L.G.zw,'sr','MarkerFaceColor','r');  % coil windings
    plot(L.G.rv,L.G.zv,'ob'); % vessel
    plot(L.G.rl,L.G.zl,'.k'); % limiter
    plot(L.G.rm,L.G.zm,'sm','MarkerFaceColor','m'); % magnetic probes
    plot(L.G.rf,L.G.zf,'oc','MarkerFaceColor','c'); % flux loops
    axis equal;

and add the computational grid. Note that we use the derived quantities `L.nzx`,|L.nrx| stored in the `L` structure which are the grid sizes

    mesh(L.G.rx,L.G.zx,zeros(L.nzx,L.nrx),...
        'edgecolor','k','facecolor','none'); view(2); shg
    axis equal

## Fluxes and fields

We can use the `L` structure already to do very useful things. For example, we can compute the flux generated by each circuit on the computational grid, as follows

    Ia = zeros(L.G.na,1); % initialize current vector
    ia = 3; % index of circuit to plot
    Ia(ia) = 1000; % put 1kA of current
    Fx = reshape(L.G.Mxa*Ia,L.nzx,L.nrx); % compute flux, and reshape it into right size for plotting
    clf; hold on;
    plot(L.G.rl,L.G.zl,'.'); contour(L.rx,L.zx,Fx*1e3,21); axis equal; colorbar;
    title(sprintf('Flux [mWb] produced by a %2.2f kA current in coil %s',Ia(ia)/1e3,L.G.dima{ia}),...
        'Interpreter','none');

Note how `L.G.dima` contains a label for the circuit.

Similarly we can easily plot the radial, vertical and poloidal field, but first we need to add the relevant matrices to `G`. (They are not included by default since they are not always needed).

    L.G = meqg(L.G,L.P,'Brxa','Bzxa'); % get required matrices
    Br = reshape(L.G.Brxa*Ia,L.nzx,L.nrx);
    Bz = reshape(L.G.Bzxa*Ia,L.nzx,L.nrx);
    Bp = sqrt(Br.^2 + Bz.^2);

    clf; hold on;
    vars = {'Fx','Br','Bz','Bp'}, units = {'mWb','mT','mT','mT'};
    for ii=1:4
      subplot(2,2,ii); hold on;
      var = eval(vars{ii}); % get variable
      plot(L.G.rl,L.G.zl,'.'); contour(L.rx,L.zx,var*1e3,21); axis equal; colorbar;
      title(sprintf('%s [%s]',...
        vars{ii},units{ii}),...
        'Interpreter','none');
    end

    vars =

      1x4 cell array

        {'Fx'}    {'Br'}    {'Bz'}    {'Bp'}

## Coordinate system

Since we just plotted fluxes and fields, we should talk about coordinate systems. MEQ uses coordinate system as defined by COCOS=17 O. Sauter and S. Y. Medvedev, Comput. Phys. Commun., 184(2), 293–302, (2013)

Specifically it * Uses a right-handed coordinate system * Defines flux [in Wb] as . * Uses a right-handed toroidal coordinate system

## `meq` data

Calling a meq code with two arguments gives the input data structure `LX` required for the code to run. We separate the data from the structure `L`, where `LX` contains (time/shot-dependent) data while `L` is intended to contain fixed tokamak/parameter information. Let's use another code for this: the equilibrium reconstruction code `liuqe`

    [L,LX] = liu('ana',1,0);
    LX

    LX =

      struct with fields:

        tokamak: 'anamak'
           shot: 1
              t: 0
             Bm: [16x1 double]
             Ff: [16x1 double]
             Ft: 2.6832e-04
             Ia: [8x1 double]
             Ip: 2.0000e+05
             Is: [0x1 double]
             Iu: [200x1 double]
             Uf: [16x1 double]
             Xt: 2.6832e-04
             aW: []
             aq: []
            rBt: 1

The meaning of these quantities is documented in

    help liux

     LIUX  LIUQE diagnostic measurements
      LX = LIUX<TOK>(SHOT,T,L) returns a structure with diagnostic measurements
      at time T for SHOT

      .t     [s] (t)   Time base
      .Ff   [Wb] (:,t) Flux loop flux
      .Bm    [T} (:,t) Magnetic probes
      .Ia    [A} (:,t) Active coil currents
      .Uf    [V] (:,t) Flux loop voltage
      .Ft   [Wb] (1,t) Diamagnetic flux
      .Xt   [Wb] (1,t) Diamagnetic flux measurement (excludes Iy correction)
      .rBt [m T] (1,t) Vaccum toroidal field *r
      .Is    [A] (:,t) Vessel segment current
      .Ip    [A] (1,:) Plasma current from discrete Rogowski

      [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

In the case of `liuqe`, the data structure contains measurements that can be used to reconstruct the equilibrium, for example `Bm` are the fields measured at the probes and `Ff` are fluxes measured at flux loops. From `LX.t` you can see that only one time point is provided for this (ficticious) shot.

Let's plot some quantities

    figure; hold on;
    subplot(131);
    barh(1:L.G.nm,LX.Bm);
    set(gca,'YTick',1:L.G.nm,'YTickLabels',L.G.dimm,'TickLabelInterpreter','none');
    subplot(132);
    barh(1:L.G.nf,LX.Ff);
    set(gca,'YTick',1:L.G.nf,'YTickLabels',L.G.dimf,'TickLabelInterpreter','none');
    subplot(133);
    barh(1:L.G.na,LX.Ia);
    set(gca,'YTick',1:L.G.na,'YTickLabels',L.G.dima,'TickLabelInterpreter','none');

With our previous example, we can plot the vacuum flux generated by the measured coil currents, showing a dominantly vertical field.

    Fx = reshape(L.G.Mxa*LX.Ia,L.nzx,L.nrx);
    figure; hold on;
    plot([L.G.rl;L.G.rl(1)],[L.G.zl;L.G.zl(1)],'-k'); % plot closed limiter contour
    contour(L.rx,L.zx,Fx,21); axis equal;

## My first meq equilibrium

Let's now calculate our first meq equilibrium using the `FBT` code. As you know, this code determines the coil currents required to maintain a desired equilibrium. We now call the function with 3 output arguments to force the calculation.

    [L,~,LY] = fbt('ana',1);

We now get an output structure `LY` that contains the equilibrium. In this case only 1 time slice. The parameters contained in `LY` are documented in `meqt.m` for generic quantities, and |<code>t.m| for code-specific quantities, for example `fbtt.m`

    help meqt
    help fbtt

     MEQT helper function, lists meaning of output variables stored in LY
      Code-specific output quantities are documented in code-specific files,
      see  FBTT, LIUT, FGET, FGST

       Dimensions: Q: normalized sqrt poloidal flux including axis (L.pQ)
                   q: normalized sqrt poloidal flux excluding axis (L.pq)
                   rx,zx: computational grid L.rx,L.zx
                   ry,zy: inner computational grid excluding boundary: L.ry,L.zy
                   o: poloidal angle theta grid (L.oq)
                   g: basis functions
                   D: plasma domain (>1 for e.g. doublets)
                   W: wall gaps, nFW: number of flux values being tracked
                   R: rational 1/q surfaces, nR: number of 1/q surfaces being tracked
                   S: flux surfaces with requested outboard r/a values, nS: number of tracked values
                   t: time

      NB: In the code the q and Q suffixes can be used to designate 2D (rho,theta) quantities.
          For clarity they are broken down here into q/Q for rho and o for theta.

      .shot  Shot number
      .t     Time base                                                   (t)       [s]

      Main equilibrium and integral quantities:
      .ag    Basis function coefficients                                 (g,t)     [-]
      .Bm    Simulated magnetic probe measurements                       (*,t)     [T]
      .Um    Simulated magnetic field time derivative measurements       (*,t)     [T/s]
      .bp    Beta poloidal                                               (t)       [-]
      .bpli2 Betap + li/2                                                (t)       [-]
      .bt    Beta toroidal                                               (t)       [-]
      .Ff    Simulated flux loop poloidal flux                           (*,t)     [Wb]
      .Uf    Simulated flux loop voltage measurements (=d/dt(Ff))        (*,t)     [V]
      .FR    Ratio of normalized X-point/boundary flux                   (D,t)     [-]
      .Ft0   Vacuum contribution to the toroidal flux                    (t)       [Wb]
      .Ft    Plasma contribution to the toroidal flux                    (t)       [Wb]
      .Fx    Plasma poloidal flux map on                                 (rx,zx,t) [Wb]
      .Ia    Fitted poloidal field coil currents                         (*,t)     [A]
      .Ip    Total plasma current                                        (t)       [A]
      .Is    Vessel segment currents                                     (*,t)     [A]
      .Iu    Vessel currents in generalized description                  (*,t)     [A]
      .Iv    Vessel filament currents                                    (*,t)     [A]
      .Iy    Plasma current distribution                                 (ry,zy,t) [A]
      .Opy   Plasma domain index                                         (ry,zy,t) [-]
      .li    Normalised internal inductance                              (t)       [-]
      .mu    Normalised diamagnetism                                     (t)       [-]
      .nA    Number of magnetic axes found                               (t)       [-]
      .PA    Pressure on axis                                            (D,t)     [Pa]
      .qA    q on axis using small diamagnetic approximation             (D,t)     [-]
      .FA    Poloidal flux on magnetic axis                              (D,t)     [Wb]
      .rA    r position of magnetic axis                                 (D,t)     [m]
      .zA    z position of magnetic axis                                 (D,t)     [m]
      .dr2FA 2nd derivative wrt r of poloidal flux on axis               (D,t)     [Wb/m^2]
      .drzFA 2nd derivative wrt r,z of poloidal flux on axis             (D,t)     [Wb/m^2]
      .dz2FA 2nd derivative wrt z of poloidal flux on axis               (D,t)     [Wb/m^2]
      .nB    Number of plasma domain boundaries found                    (t)       [-]
      .lB    Boolean flag indicating if domain boundary has been found   (D,t)     [-]
      .FB    Poloidal flux at pt. defining plasma domain boundary        (D,t)     [Wb]
      .rB    r position of point defining plasma domain boundary         (D,t)     [m]
      .zB    z position of point defining plasma domain boundary         (D,t)     [m]
      .F0    Poloidal flux at domain center point                        (D,t)     [Wb]
      .F1    Poloidal flux at domain boundary point                      (D,t)     [Wb]
      .rBt   Vacuum toroidal field * r                                   (t)       [T.m]
      .rIp   current distribution weighted by r grid                     (t)       [A.m]
      .zIp   current distribution weighted by z grid                     (t)       [A.m]
      .rY    r position of centroid of current distribution (=rIp/Ip)    (t)       [m]
      .zY    z position of centroid of current distribution (=zIp/Ip)    (t)       [m]
      .lX    Boolean flag indicating if domain is diverted               (D,t)     [-]
      .nX    Number of X points found                                    (t)       [-]
      .FX    Flux at X points                                            (*,t)     [Wb]
      .rX    r position of X points                                      (*,t)     [m]
      .zX    z position of X points                                      (*,t)     [m]
      .dr2FX 2nd derivative wrt r of poloidal flux at X-points           (*,t)     [Wb/m^2]
      .drzFX 2nd derivative wrt r,z of poloidal flux at X-points         (*,t)     [Wb/m^2]
      .dz2FX 2nd derivative wrt z of poloidal flux at X-points           (*,t)     [Wb/m^2]
      .PQ    Pressure profile                                            (Q,D,t)   [Pa]
      .TQ    T = r*Btor profile                                          (Q,t)     [T.m]
      .iTQ   1/TQ (if P.smalldia: with small diagmagnetic approximation) (Q,t)     [1/(T.m)]
      .Wk    Plasma kinetic energy                                       (t)       [J]
      .WN    Normalisation energy                                        (t)       [J]
      .Wp    Plasma poloidal magnetic energy                             (t)       [J]
      .Wt0   Vacuum toroidal magnetic energy                             (t)       [J]
      .Wt    Plasma toroidal magnetic energy                             (t)       [J]
      .Vp    Plasma volume                                               (t)       [m^3]

      % Integral quantities per plasma domain
      .bpD   Beta poloidal per domain                                    (D,t)     [-]
      .btD   Beta toroidal per domain                                    (D,t)     [-]
      .FtD   Plasma toroidal flux per domain                             (D,t)     [Wb]
      .Ft0D  Vacuum toroidal flux per domain                             (D,t)     [Wb]
      .IpD   Plasma current per domain                                   (D,t)     [A]
      .liD   Normalised internal inductance per domain                   (D,t)     [-]
      .TpDg  Integral of basis function values per domain                (g,D,t)   [-]
      .ITpDg Integral of primitives of basis function values per domain  (g,D,t)   [-]
      .rYD   r position of centroid of current distribution (=rIpD./IpD) (D,t)     [m]
      .zYD   z position of centroid of current distribution (=zIpD./IpD) (D,t)     [m]
      .rIpD  current distribution weighted on r grid                     (D,t)     [A.m]
      .zIpD  current distribution weighted on z grid                     (D,t)     [A.m]
      .VpD   Plasma volume per domain                                    (D,t)     [m^3]
      .WkD   Plasma kinetic energy per domain                            (D,t)     [J]
      .WpD   Plasma poloidal magnetic energy per domain                  (D,t)     [J]
      .WtD   Plasma toroidal magnetic energy per domain                  (D,t)     [J]
      .Wt0D  Vacuum toroidal magnetic energy per domain                  (D,t)     [J]
      .WND   Normalisation energy per domain                             (D,t)     [J]

      % Integral quantities on rho grid per plasma domain
      .FtQ   Plasma toroidal flux contained in flux surface              (Q,D,t)   [Wb]
      .Ft0Q  Vacuum toroidal flux contained in flux surface              (Q,D,t)   [Wb]
      .FtPQ  Total toroidal flux contained in flux surface (FtQ+Ft0Q)    (Q,D,t)   [Wb]
      .IpQ   Plasma current contained in flux surface                    (Q,D,t)   [A]
      .OpQ   Number of grid points contained in flux surface             (Q,t)     [-]
      .PpQg  p' base functions multiplied by their coefficients          (Q,D,g,t) [Pa/Wb]
      .PpQ   p' profile                                                  (Q,D,t)   [Pa/Wb]
      .TTpQg TT' base functions multiplied by their coefficients         (Q,D,g,t) [T^2.m^2/Wb]
      .TTpQ  TT' profile                                                 (Q,D,t)   [T^2.m^2/Wb]
      .VpQ   Plasma volume contained in flux surface                     (Q,D,t)   [m^3]
      .WkQ   Plasma kinetic energy contained in flux surface             (Q,D,t)   [J]
      .WpQ   Plasma poloidal magnetic energy contained in flux surface   (Q,D,t)   [J]

      % Contour-related quantities on rho grid per plasma domain D, only available for iterq>0
      .aq      distance between point on flux surface and magnetic axis  (o,q,D,t) [m]
      .zq      z position of flux surface points                         (o,q,D,t) [m]
      .rq      r position of flux surface points                         (o,q,D,t) [m]
      .raQ     Outboard r/A value on pQ grid                             (Q,D,t)   [-]
      .iqQ     1/q profile where q is the safety factor                  (Q,D,t)   [-]
      .jtorQ   Toroidal current density defined as R0*<jphi/R>           (Q,D,t)   [A/m2]
      .raqmin  Outboard r/a value of location of minimum safety factor   (t)       [-]
      .qmin    Minimum safety factor                                     (t)       [-]
      .q95     safety factor at 95% normalized poloidal flux surface     (D,t)       [-]
      .Q0Q     <1/R>                                                     (Q,D,t)   [1/m]
      .Q1Q     -dpsi/dV                                                  (Q,D,t)   [T/m]
      .Q2Q     <1/R^2>                                                   (Q,D,t)   [1/m^2]
      .Q3Q     <|grad psi|^2/R^2>                                        (Q,D,t)   [T^2]
      .Q4Q     <|grad psi|^2>                                            (Q,D,t)   [T^2.m^2]
      .Q5Q     <|grad psi|/(2pi)>                                        (Q,D,t)   [T.m]
      .ItQ     Plasma current contained in flux surface                  (Q,D,t)   [A]
      .LpQ     Length of the poloidal contour: \oint dl                  (Q,D,t)   [m]
      .rbQ     Radial location of contour barycenter: (\oint R dl) / dl  (Q,D,t)   [m]
      .VQ      Plasma volume contained in flux surface                   (Q,D,t)   [m^3]
      .AQ      Cross-section area of flux-surface                        (Q,D,t)   [m^2]
      .SlQ     Lateral area of flux-surface                              (Q,D,t)   [m^2]

      .raR     Outboard r/a value of location of requested 1/q values    (R,t)     [-]
      .rS      r position of flux surfaces with requested outboard r/a   (S,t)     [-]
      .zS      z position of flux surfaces with requested outboard r/a   (S,t)

      % Shape related quantities, only available for iterq>0
      .aminor  Minor radius: (Rmax-Rmin)/2                               (Q,D,t)   [m]
      .epsilon Aspect ratio: aminor/rgeom                                (Q,D,t)   [-]
      .delta   Triangularity: (deltal+deltau)/2                          (Q,D,t)   [-]
      .deltal  Lower triangularity: (r(zmin )-rgeom)/aminor              (Q,D,t)   [-]
      .deltau  Upper triangularity: (r(zmmax)-rgeom)/aminor              (Q,D,t)   [-]
      .kappa   Elongation: (Zmax-Zmin)/(Rmax-Rmin)                       (Q,D,t)   [-]
      .lp      Length of the LCFS poloidal contour: \oint dl             (D,t)     [m]
      .rbary   Radial location of LCFS barycenter: (\oint R dl) / dl     (D,t)     [m]
      .rgeom   r position of geom. center of flux surface: (Rmax+Rmin)/2 (Q,D,t)   [m]
      .zgeom   z position of geom. center of flux surface: (Zmax+Zmin)/2 (Q,D,t)   [m]

      % Wall gap information available only for iterq>0
      .aW      Gap distance (defined by L.G.rW,.zW and angle .oW)        (W,nFW,t) [m]
      .FW      Flux values that are being tracked by gap algorithm       (1,nFW,t) [Wb]

      % Magnetic fields (ifield=true) and vacuum fluxes (ivacuum=true)
      % on x and z grids (z if izgrid=true)
      .Brx   Radial magnetic field                                       (rx,zx,t) [T]
      .Bzx   Vertical magnetic field                                     (rx,zx,t) [T]
      .Brz   Radial magnetic field on z grid                             (rz,zz,t) [T]
      .Bzz   Vertical magnetic field on z grid                           (rz,zz,t) [T]
      .Br0x  Radial magnetic field                                       (rx,zx,t) [T]
      .Bz0x  Vertical magnetic field                                     (rx,zx,t) [T]
      .Br0z  Radial magnetic field on zgrid                              (rz,zz,t) [T]
      .Bz0z  Vertical magnetic field on zgrid                            (rz,zz,t) [T]
      .Fx0   Vacuum flux                                                 (rx,zx,t) [Wb]
      .Fz0   Vacuum flux on z grid                                       (rz,zz,t) [Wb]

      % Quantities on interpolation points (when using 'infct')
      .Fn    Flux at interpolation points L.P.rn, L.P.zn                 (*,t)     [Wb]
      .Un    Loop voltage at interpolation points                        (*,t)     [V]
      .Brn   Br at interpolation points                                  (*,t)     [T]
      .Bzn   Bz at interpolation points                                  (*,t)     [T]
      .Brrn  dBr/dr at interpolation points                              (*,t)     [T/m]
      .Brzn  dBr/dz at interpolation points                              (*,t)     [T/m]
      .Bzzn  dBz/dz at interpolation points                              (*,t)     [T/m]
      .Bzrn  dBz/dr at interpolation points                              (*,t)     [T/m]

      % Quantities related to evolutive equations
      .Va       Voltage for active circuits                              (*,t)       [V]
      .IniD     Non-inductive plasma current                             (D,t)       [A]

      % Solver output:
      .err   sum of error flags (slx only):
                err=(|Ip|<L.P.Ipmin)+2*(nA==0)+4*(lB==0)+8*(PA<0) (t)[-]
      .isconverged Boolean flag indicating if code has converged         (t)       [-]
      .niter Number of iterations done                                   (t)       [-]
      .res   Solver residual                                             (t)       [-]

      .rese  Residual of circuit equations                               (t)       [-]
      .resy  Residual norm of Iy                                         (t)       [-]
      .resp  Residual norm of plasma current equation                    (t)       [-]
      .resC  Residual norm of ag constraints                             (t)       [-]
      .resFx Residual norm of Fx                                         (t)       [-]

      [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

     FBTT  FBT inverse equilibrium solver
      LY = FBTT(L,LX[,'PAR',VAL,...]) calculates equilibrium using the parameters
      in L as obtained by FBTC, optionally replacing or adding parameters with
      specified values. It returns a structure LY with the calculated
      equilibria.

      For help on output structure: see meqt.m
      FBT-specific outputs:

      .Sp    Sum of surfaces of all grid points containing plasma [m^2]
      .Ep    Plasma box elongation
      .chi   Root-mean-square of cost function residual
      .chie  Root-mean-square of constraint residual
      .werr  Weighted cost function error
      .uerr  Unweighted cost function error
      .ID    Equation id for each cost function term
      .Fb    Boundary flux (fitting parameter)

      The current version extends both the original FBT algorithm [Hofmann CPC
      1988] as well as the GSPD algorithm [J. Wai et al, arXiv 2023].
      The optimization algorithm is described in MEQOPT.

      See also FBTLSQINIT FBTINIT MEQOPT

      [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

we can plot the key quantities for illustration:

    figure;
    subplot(121); hold on;
    contourf(L.rx,L.zx,LY.Fx,21);
    contour(L.rx,L.zx,LY.Fx,LY.FB*[1 1],'linecolor','w','linewidth',2); % LCFS
    plot(LY.rA,LY.zA,'or'); % magnetic axis
    plot([L.G.rl;L.G.rl(1)],[L.G.zl;L.G.zl(1)],'-k','linewidth',2); % plot closed limiter contour
    axis equal; title('flux distribution')
    subplot(122); hold on; axis xy
    imagesc(L.ry,L.zy,LY.Iy); % current distribution (note different grid)
    contour(L.rx,L.zx,LY.Fx,LY.FB*[1 1],'linecolor','w','linewidth',2); % LCFS
    plot(LY.rIp/LY.Ip,LY.zIp/LY.Ip,'*r'); % current centroid
    plot([L.G.rl;L.G.rl(1)],[L.G.zl;L.G.zl(1)],'-k','linewidth',2); % plot closed limiter contour
    title('current distribution');
    axis equal;

As you can see this is a rather dull circular equilibrium. Let's plot something more exciting, and use the function `meqplott` that does all the plotting for us.

    [L,~,LY] = fbt('ana',2);
    figure; meqplott(L,LY);

If you want to see more details of the flux in the region around the coils, as well as the coil currents, use `meqplotfancy` function. This needs a specific parameter izgrid to be set first to compute the flux on a secondary grid.

    [L,~,LY] = fbt('ana',2,[],'izgrid',true);
    figure; meqplotfancy(L,LY);