# 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
      .debug        debug flag
      .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)

      .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
      .iqR          Value of 1/q for which to return r/a location (only if iterq>0)
      .naR          Number of ocurrences 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

      .icsint       If icsint, use cubic spline interpolation for domain id. and contouring (default: false)
      .tolcs        Tolerance on position of axes, saddle points, limiter points and flux-surface contours [active only if icsint is true]
      .itercs       Maximum number of Newton-Raphson iterations                                            [active only if icsint is true]

      .mdsserver    MDS server hostname for parameters/measurements loading

      .gsxe         Switch for treatment of external currents inside computational grid
      .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

      [+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:
      .agfitfct    Fitting function handle for basis function coefficients
      .agfitp      Parameters for fitting function
      .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.
      .dissi       High dissi promotes low coil currents   (in gpid)
      .dipol       High dipol promotes low dipole currents (in gpdd)
      .mergex      If true, merge redundant X-point constraints in LX

      Quantities which will end up in LX
      .tol         Default tol
      .issym       If true, force up-down summetric current distribution
      .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 =

      1×4 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: [16×1 double]
             Ff: [16×1 double]
             Ft: 2.6795e-04
             Ia: [8×1 double]
             Ip: 2.0000e+05
            IpD: 2.0000e+05
             Iu: [200×1 double]
             aW: []
             ag: [3×1 double]
             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
      .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
                   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,D,t)   [-]
      .Bm    Fitted magnetic probe measurements                          (*,t)     [T]
      .bp    Beta poloidal                                               (t)       [-]
      .bpli2 Betap + li/2                                                (t)       [-]
      .bt    Beta toroidal                                               (t)       [-]
      .Ff    Fitted flux loop poloidal flux                              (*,t)     [Wb]
      .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]
      .rBt   Vacuum toroidal field * r                                   (t)       [T.m]
      .rIp   R of centroid of current distribution                       (t)       [m]
      .zIp   z position of centroid of current distribution              (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]
      .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)   [-]
      .rIpD  r position of centroid of current distribution              (D,t)     [m]
      .zIpD  z position of centroid of current distribution              (D,t)     [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
              for first domain                                           (Q,g,t)   [Pa/Wb]
      .PpQ   p' profile                                                  (Q,D,t)   [Pa/Wb]
      .TTpQg TT' base functions multiplied by their coefficients
              for first domain                                           (Q,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  (q,o,D,t) [m]
      .zq      z position of flux surface points                         (q,o,D,t) [m]
      .rq      r position of flux surface points                         (q,o,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     (t)       [-]
      .Q0Q     <1/R>                                                     (Q,t)     [1/m]
      .Q1Q     -dpsi/dV                                                  (Q,t)     [T/m]
      .Q2Q     <1/R^2>                                                   (Q,t)     [1/m^2]
      .Q3Q     <(grad psi)^2/R^2>                                        (Q,t)     [T^2]
      .Q4Q     <(grad psi)^2>                                            (Q,t)     [T^2.m^2]

      % Shape related quantities, only availble 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]
      .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]

      % 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)       [-]

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

     FBTT  FBTE equilibrium calculation
      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   Cost function residual
      .chie  Equality constraints global residual
      .Fb    Boundary flux (fitting parameter)

      [+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],'color','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],'color','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);