## Contents

- Eleven FBT equilibria for Anamak
- 1: My first FBT equilibrium
- 2: Force points to be exactly on the boundary
- 3: Constrain Flux value on one point
- 4: Constrain one current value
- 5: X point(s)
- 6: Add a strike point
- 7: Add magnetic field angle
- 8: Add flux expansion in radial direction
- 9: Zero hessian at X point (snowflake)
- 10: Vacuum field curvature
- 11: Apply current limits
- Bonus: Internal profiles constraints in FBT

## Eleven FBT equilibria for Anamak

This tutorial illustrates the use of constraints to specify equilibria in FBT.

The function of `fbt`, the meaning of various variables, as well as the setup of the constrained optimization problem are described by

    help fbthelp
    % Initialize L structure and weights for coil
    % We first initialize an 'empty' case with no constraints

    shot = 0; t=0; % loads empty case - no shot parameters defined
    [L,LX0] = fbt('ana',shot,t,...
      'izgrid',1,... % add extended grid to plot flux beyond limiter
      'iterq',50,... % switch on post processing with 50 iterations
      'pq',[],'npq',10);  % Override default and set number of rho points for contours

      FBTHELP - Help for FBT

      FBT solves the inverse static problem:
      * Given:
        * Cost function weights and constraints on flux, field, currents
        * Internal plasma constraints (e.g. bp, Ip, qA)
        * Inequality constraints on (combinations of) Ia.
      * Finds:
        * Ia and corresponding Equilibrium minimizing cost function while
          satisfying constraints

      Setup of the constrained optimization problem:
      The fitting parameters on the left hand side are Ia and Fb (coil currents
      and an offset flux).

      On the right side are the user-prescribed values and the plasma contribution Iy
      The latter is iterated as it depends on Ia

      equals(=) sign means equality constraint if weight=0 (1/weight is infinity)
      Otherwise the equation is solved with the others in the least-squares
      sense

      1./(gpid*gpie) * [ Ia                    = gpia                          ]  % Currents
      1./(gpdd*gpde) * [ Dda*(gpdw.*Ia)        = 0                             ]  % Dipoles
      1./(gpfd*gpfe) * [ MSa*Ia   - gpfb*Fb    = gpfa  - MSy*Iy                ]  % Flux
      1./(gpbd*gpbe) * [ BrSa*Ia               = gpbr  - BrSa*Iy               ]  % Br
      1./(gpbd*gpbe) * [ BzSa*Ia               = gpbz  - BzSa*Iy               ]  % Bz
      1./(gpbd*gpbe) * [ Gba*[BSra;BSza]*Ia    = 0     - Gba*[BSry;BSzy]*Iy    ]  % Field angle
      1./(gpcd*gpce) * [ MrrSa*Ia              = gpcr  - MrrSy*Iy              ]  % Frr constraint
      1./(gpcd*gpce) * [ MrzSa*Ia              = gpcz  - MrzSy*Iy              ]  % Frz constraint
      1./(gpcd*gpca) * [ Gca*[MrrSa;MrzSa]*Ia  = 0     - Gca*[MrrSy;MrzSy]*Iy  ]  % Hessian angle
      1./(gpvd*gpve) * [ BrrSa*Ia              = gpvrr                         ]  % Vacuum Brr constraint
      1./(gpcd*gpve) * [ BrzSa*Ia              = gpvrz                         ]  % Vacuum Brz constraint
      1./(gpcd*gpve) * [ BzzSa*Ia              = gpvzz                         ]  % Vacuum Bzz constraint

      Dda Matrix such that Dda*Ia gives current differences e.g. [Ia(2)-Ia(1);...;Ia(na)-Ia(na-1)];
      Gba = [sin(gpba),-cos(gpba)];
      Gca = [sin(gpca),-cos(gpca)];

      Dda is contained in G, while all other matrices are calculated in each time loop by fbtt.m

      If one of the Br/Bz or Mrr,Mrz components are constrained, then
      the corresponding angle constraint is automatically ignored.

      Control point constraints gpf* gpb* gpc* gpv* can be supplied using FBTGP()
      See FBTGP for details on the syntax

      NOTE: Fb and FB (boundary flux) are not necessarily the same, unless the
      solution has LCFS passing through a point with gpfb=1 and gpfa=0 and gpfe*gpfd=0.

      Inequality constraints:
      liml*limm < limc*Ia < limu*limm
      limm is a margin, limu,liml are upper/lower bounds, limc is a matrix.
      All are contained in P

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

We then set the global weights for the currents

    % % global weights
    LX0.gpie = ones(L.G.na,1); % pf currents, one per coil
    LX0.gpfd = 1e-4; % Global weights for flux points
    LX0.gpid = 1e+3; % Global weight for Current constraints
    LX0.gpbd = 1e-6; % Global weight for Field constraints
    LX0.gpcd = 1e-6; % Global weight for Hessian constraints

## 1: My first FBT equilibrium

Let's make our first equilibrium by specifying three points on which we want the flux value to be the same. We specify this as a cost function term. We use the auxiliary function `fbtgp` to populate these parameters.

    help fbtgp

      X = fbtgp(X,r,z,b,fa,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve)

      Appends entry to the set of gp* variables describing geometry constraints
      for FBT. For details of the equations, see FBTHELP

      X: previous structure
      r,z: r,z location of point(s) to add
      b  : Point is on main plasma boundary (used only for initial guess of current distribution)
      fa : Flux value
      fb : if 0, then fa is the absolute flux value; if 1, then fa is defined relative to an unknown offset common to all points with fb=1
      fe : Weight of flux constraint (multiplies fd)
      br : Value of radial magnetic field
      bz : Value of vertical magnetic field
      ba : Angle of magnetic field
      be : Weight of magnetic field constraint (multiplies bd)
      cr : Value of d2/dr2(Psi)
      cz : Value of d2/drdz(Psi)
      ca : [To be clarified] Angle of hessian of Psi
      ce : Weight of hessian of Psi constraints (multiplies cd)
      vrr: Value of d/dr(Br), vacuum field only
      vrz: Value of d/dr(Bz), vacuum field only
      vzz: Value of d/dz(Bz), vacuum field only
      ve : Weight of vacuum field derivative constraints (multiplies vd)

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

We specify a set of points on a circle, which will correspond to a set of cost function terms encouraging the flux to be equal on all these points

    th = (1:10)'/10*2*pi;
    a0 = 0.4;
    rr = 1 + a0.*cos(th);
    zz =   + a0.*sin(th);
    gpb  = 1; % the points are on the main plasma domain boundary
    gpfa = 0; % relative flux offset
    gpfb = 1; % Scale w.r.t. boundary flux
    gpfe = 1; % These are all cost function terms and not equality constraints

    % assign to P structure
    LX = LX0; % init
    LX = fbtgp(LX,rr,zz,gpb,gpfa,gpfb,gpfe,[],[],[],[],[],[],[],[],[],[],[],[]);

The resulting `LX` structure contains the cost function and constraint terms that will be used in the equilibrium calculation in `fbtt`. Its contents can be displayed using `fbtxdisp`

    fbtxdisp(L,LX);

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.324 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.600 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.324 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.400 z=-0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04

We now run the optimization and plot the result.

    LY  = fbtt(L,LX);
    clf; fbtplot(L,LX,LY);snapnow;

Note that the plasma boundary flux `FB` is not necessarily the same as the fitting parameter `Fb` because none of the control points ended up being on the LCFS

    disp([LY.Fb, LY.FB])

        0.0768    0.0606

Store this for later

    LX1 = LX;

## 2: Force points to be exactly on the boundary

To force a point to be on the boundary two conditions must be met: * The point is the one defining the separatrix, e.g. a limiter point or an x point * The constraints must be exact.

We set an equality constraint by `gpfe=0`.

    LX = LX1; % init
    LX = fbtgp(LX,0.55,0,1,0,1,0,[],[],[],[],[],[],[],[],[],[],[],[]);
    fbtxdisp(L,LX);
    LY  = fbtt(L,LX);

    clf;fbtplot(L,LX,LY);snapnow;

    disp([LY.Fb, LY.FB])

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.324 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.600 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.324 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.400 z=-0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.550 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00

        0.0654    0.0645

## 3: Constrain Flux value on one point

We may want to impose the absolute value of the flux at one or more points. We can do this by setting `gpfb=0` for that point. In this example we set the absolute value for the limiter point flux. In order to make sure the other points have the same flux, we also set a constraint with `gpfb=1` on the same r,z point. So the same point appears twice in the `fbtgp` call.

    FBB = 0.2; % target boundary flux value
    LX = LX1; % init
    %             r   ,z,b,fa ,fb,fe,br,bz,ba,be,cr,cz,ca,ce,vrr,vrz,vzz,ve
    LX = fbtgp(LX,0.55,0,1,FBB, 0, 0,[],[],[],[],[],[],[],[], [], [], [],[]);
    LX = fbtgp(LX,0.55,0,1,  0, 1, 0,[],[],[],[],[],[],[],[], [], [], [],[]);

    fbtxdisp(L,LX);
    LY  = fbtt(L,LX);

    clf;fbtplot(L,LX,LY);snapnow;
    disp([LY.Fb, LY.FB])
    LX3 = LX; % store for later

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.324 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.600 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.324 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.400 z=-0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.550 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =   0.2 +    0*Fb                 1/w=0.00e+00
                        B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00

        0.2000    0.1991

## 4: Constrain one current value

We can constrain one current value by setting `gpie=0` and imposing the value via `gpia`

    LX = LX1; % Start again from case 1.

    ia_set = 2; % index of coil current to set
    LX.gpie(ia_set) = 0;     % Set Ia(2)=gpia(2) to be an equality constraint
    LX.gpia(ia_set) = +0e3; % Assign target value

    fbtxdisp(L,LX);
    % note that 1/w for PF_002 is zero.

    LY  = fbtt(L,LX);
    % Check result
    fprintf('Ia(2)=%2.2g',LY.Ia(ia_set))

    clf;fbtplot(L,LX,LY);snapnow;

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=0.00e+00
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.324 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.600 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.324 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.400 z=-0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    Ia(2)= 0

## 5: X point(s)

We now generate a new point distribution, where two points are X points with Br=Bz=0. One of them is set to be on the same surface as the other LCFS point, while one is given an offset. The one with an offset is assigned gpb=0

    LX = LX0;

    % Elongated distribution of points
    a0 = 0.3; kap = 1.3; % minor radius, elongation
    th = (1:10)'/10*2*pi+0.2;
    rr = 1    +     a0.*cos(th);
    zz = 0.05 + kap*a0.*sin(th);

    gpb  = 1; % the points are on the main plasma domain boundary
    gpfa = 0; % relative flux offset
    gpfb = 1; % Scale w.r.t. boundary flux
    gpfe = 1; % These are all cost function terms and not equality constraints

    % assign to P structure
    % Elongated set of points
    %               r,    z,  b,  fa,  fb,  fe,  br,bz,ba,be,  cr,cz,ca,ce,vrr,vrz,vzz,ve)
    LX = fbtgp(LX, rr,   zz,gpb,gpfa,gpfb,gpfe,  [],[],[],[],  [],[],[],[], [], [], [],[]);
    % Primary X-point with exact flux constraint
    LX = fbtgp(LX,0.8,-0.35,  1,   0,   1,   0,   0, 0,[], 0,  [],[],[],[], [], [], [],[]);

    % Second x-point outside vacuum vessel, at specified flux offset w.r.t. primary x point
    %                r,  z, b, fa,fb,fe,  br,bz,ba,be,  cr,cz,ca,ce,vrr,vrz,vzz,ve)
    FFB= -0.02; % flux ofset value
    LX = fbtgp(LX,0.95,0.6, 0,FFB, 1, 0,   0, 0,[], 0,  [],[],[],[], [], [], [],[]);

    fbtxdisp(L,LX);

    LY  = fbtt(L,LX);
    LX5 = LX; % save for later
    clf;fbtplot(L,LX,LY);snapnow;

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.203 z=+0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.034 z=+0.437   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.852 z=+0.390   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.727 z=+0.212   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.706 z=-0.027   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.797 z=-0.237   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.966 z=-0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.148 z=-0.290   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.273 z=-0.112   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.294 z=+0.127   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.800 z=-0.350   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
    r=+0.950 z=+0.600     Psi:   MSa*Ia +  MSy*Iy = -0.02 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00

## 6: Add a strike point

A strike point has the same flux as the other LCFS points, but is not on the main plasma domain boundary

    % assign to P structure
    LX = LX5; % init
    LX = fbtgp(LX,0.82,-0.45,0,0,1,0,[],[],[],[],[],[],[],[],[],[],[],[]);
    fbtxdisp(L,LX);
    LY  = fbtt(L,LX);

    clf;fbtplot(L,LX,LY);snapnow;

    LX6 = LX; % save for later

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.203 z=+0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.034 z=+0.437   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.852 z=+0.390   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.727 z=+0.212   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.706 z=-0.027   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.797 z=-0.237   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.966 z=-0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.148 z=-0.290   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.273 z=-0.112   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.294 z=+0.127   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.800 z=-0.350   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
    r=+0.950 z=+0.600     Psi:   MSa*Ia +  MSy*Iy = -0.02 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
    r=+0.820 z=-0.450     Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00

## 7: Add magnetic field angle

We now force the angle of the magnetic field on the strike point To allow this constraint to be satisfied exctly we also release the equality constraint on the x point.

    LX = LX6; % start from example 6

    rs=LX.gpr(end); zs=LX.gpz(end); % strike point r,z
    % Add another constraint on the last point, to force the magnetic field angle value
    %              r, z, b,  fa,  fb,  fe,  br,  bz,        ba,    be,cr,cz,ca,ce,vrr,vrz,vzz,ve
    LX = fbtgp(LX,rs,zs, 0, NaN, NaN, NaN,  NaN,NaN,   0.45*pi,     0,[],[],[],[], [], [], [],[]);

    ix=(LX.gpbr==0) | LX.gpbz==0; % indices of x point field constriants
    LX.gpbe(ix)=1; % Set cost function term instead of equality constraint
    fbtxdisp(L,LX);
    LY  = fbtt(L,LX);

    clf;fbtplot(L,LX,LY);snapnow;

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.203 z=+0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.034 z=+0.437   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.852 z=+0.390   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.727 z=+0.212   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.706 z=-0.027   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.797 z=-0.237   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.966 z=-0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.148 z=-0.290   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.273 z=-0.112   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.294 z=+0.127   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.800 z=-0.350   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=1.00e-06
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=1.00e-06
    r=+0.950 z=+0.600     Psi:   MSa*Ia +  MSy*Iy = -0.02 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=1.00e-06
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=1.00e-06
    r=+0.820 z=-0.450     Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                           BN:  1.0*(BrSa*Ia+BrSy*Iy) +  0.2*(BzSa*Ia+BzSy*Iy) = 0  1/w=0.00e+00

## 8: Add flux expansion in radial direction

    LX=LX6;
    % Add constraint for psirr (d2/dr2 Psi) at the strike point point
    %             r,  z,b,fa,fb,fe, br, bz, ba, be,cr, cz, ca,ce,vrr,vrz,vzz,ve
    LX = fbtgp(LX,rs,zs,0, 0, 1, 0,NaN,NaN,NaN,NaN,-1.2,NaN,NaN, 0, [], [], [],[]);
    fbtxdisp(L,LX);
    LY  = fbtt(L,LX);
    clf;fbtplot(L,LX,LY);snapnow;

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.203 z=+0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.034 z=+0.437   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.852 z=+0.390   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.727 z=+0.212   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.706 z=-0.027   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.797 z=-0.237   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.966 z=-0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.148 z=-0.290   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.273 z=-0.112   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.294 z=+0.127   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.800 z=-0.350   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
    r=+0.950 z=+0.600     Psi:   MSa*Ia +  MSy*Iy = -0.02 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
    r=+0.820 z=-0.450     Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                          Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                      d2r psi: MrrSa*Ia +MrrSy*Iy =  -1.2                           1/w=0.00e+00

Note how this yields tighter flux surfaces around the divertor and removes the second x point under the vessel w.r.t. case 7.

## 9: Zero hessian at X point (snowflake)

    LX=LX0; % Start from scratch

    % Elongated distribution of points
    a0 = 0.3; kap = 1.3; % minor radius, elongation
    th = (1:10)'/10*2*pi+0.2;
    rr = 1    +     a0.*cos(th);
    zz = 0.05 + kap*a0.*sin(th);

    %             [r, z   b ,fa, fb,fe,       br, bz, ba, be, cr, cz,ca, ce ,vrr,vrz,vzz,ve]
    LX = fbtgp(LX,rr,zz,gpb,gpfa,gpfb,gpfe,   [],[],[],[],[],[],[],[],[],[],[],[]);

    % Primary monkey-saddle-point (B=gradB=0) with exact flux constraint
    % Add constraint for Br,Bz to be both zero (exact)
    % Add constraint for psirz and psirr to be both zero at that point (exact)
    %                r,    z,b,fa,fb,fe,   br,bz,ba,be,    cr,cz,ca,ce   ,vrr,vrz,vzz,ve
    LX = fbtgp(LX,0.84,-0.48,1, 0, 1, 0,    0 ,0,[], 0,     0, 0,[], 0   , [], [], [],[]);
    fbtxdisp(L,LX);
    LY  = fbtt(L,LX);
    clf;fbtplot(L,LX,LY);snapnow;

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.203 z=+0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.034 z=+0.437   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.852 z=+0.390   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.727 z=+0.212   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.706 z=-0.027   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.797 z=-0.237   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.966 z=-0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.148 z=-0.290   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.273 z=-0.112   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.294 z=+0.127   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.840 z=-0.480   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
                      d2r psi: MrrSa*Ia +MrrSy*Iy =     0                           1/w=0.00e+00
                     drdz psi: MrzSa*Ia +MrzSy*Iy =     0                           1/w=0.00e+00

Two observations:

- It is enough to constrain psirr,psirz to be 0 to ensure that gradB=0 because the grad-shafranov equation will give us the final constraint psizz=0 if we assume a zero current density at the X-point. This is ensured by this particular choice of basis functions, but may not hold in general.
- The resulting equilibria will possibly not have an exact monkey point because of the difference in the way we impose the contraints (expressing psi as a combination of green functions) and the way we actually solve for psi (by inverting the GS operator).

## 10: Vacuum field curvature

    LX = LX3; % start from limited equilibrium
    LX.gpvd = 1;
    %             r,z, b, fa, fb, fe, br, bz, ba, be, cr, cz, ca, ce,vrr,  vrz,vzz,ve
    LX = fbtgp(LX,1,0, 0,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,+0.05,NaN, 0);

    fbtxdisp(L,LX);

    LY  = fbtt(L,LX);

    clf;
    subplot(121)
    fbtplot(L,LX,LY);

    subplot(122)
    % calculate and plot resulting vacuum field
    meqplott(L,LY); hold on;
    L.G = meqg(L.G,L.P,'Brxa','Bzxa');
    LY.Br0x = reshape(L.G.Brxa*LY.Ia,L.nzx,L.nrx);
    LY.Bz0x = reshape(L.G.Bzxa*LY.Ia,L.nzx,L.nrx);
    LY.F0x  = reshape(L.G.Mxa *LY.Ia,L.nzx,L.nrx);
    meqplotfield(L,LY,'vacuum',true)
    snapnow

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
       none

    Magnetic geometry cost / equality constraints:
    r=+1.324 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=+0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=+0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.600 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.676 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.876 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.124 z=-0.380   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.324 z=-0.235   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.400 z=-0.000   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.550 z=+0.000   B Psi:   MSa*Ia +  MSy*Iy =   0.2 +    0*Fb                 1/w=0.00e+00
                        B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=0.00e+00
    r=+1.000 z=+0.000  dz B0r: BrzSa*Ia           =  0.05                           1/w=0.00e+00

## 11: Apply current limits

    LX = LX6;
    L.P.limu =  300e3;
    L.P.liml = -300e3;
    L.P.limm = 0.95*ones(L.G.na,1); % margin
    L.P.limc = eye(L.G.na);
    LX.gpfe(:) = 1; % make flux constraints non-exact
    fbtxdisp(L,LX);

    LY  = fbtt(L,LX);
    LX11 = LX; % save this

    clf;
    subplot(1,3,[1,2])
    fbtplot(L,LX,LY);

    subplot(1,3,3)
    barh(LY.Ia/1e3); hold on;
    plot(L.P.limm*L.P.limu*[1 1]/1e3,[1,L.G.na]);
    plot(L.P.limm*[1,1]*L.P.liml/1e3,[1,L.G.na]);
    set(gca,'xlim',500*[-1 1]); xlabel('kA');
    set(gca,'TickLabelInterpreter','none');
    aa=set(gca,'YTickLabel',L.G.dima);
    title('Coil currents and limits')
    snapnow;
    % Notice how this makes a less shaped plasma at the expense of less
    % accurate flux point fitting

    FBTX cost function / constraints display
      Tokamak: anamak, shot#0, t=0.000

    Coil current cost / equality constraints:
    PF_001             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_002             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_003             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_004             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_005             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_006             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_007             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03
    PF_008             Current:  Ia = 0.00e+00[A]                               1/w=1.00e+03

    Dipole cost / equality constraints:
       none

    Coil current limits:
    -2.85e+05 <=     +1*PF_001     <= 2.8e+05
    -2.85e+05 <=     +1*PF_002     <= 2.8e+05
    -2.85e+05 <=     +1*PF_003     <= 2.8e+05
    -2.85e+05 <=     +1*PF_004     <= 2.8e+05
    -2.85e+05 <=     +1*PF_005     <= 2.8e+05
    -2.85e+05 <=     +1*PF_006     <= 2.8e+05
    -2.85e+05 <=     +1*PF_007     <= 2.8e+05
    -2.85e+05 <=     +1*PF_008     <= 2.8e+05

    Magnetic geometry cost / equality constraints:
    r=+1.203 z=+0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.034 z=+0.437   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.852 z=+0.390   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.727 z=+0.212   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.706 z=-0.027   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.797 z=-0.237   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.966 z=-0.337   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.148 z=-0.290   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.273 z=-0.112   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+1.294 z=+0.127   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
    r=+0.800 z=-0.350   B Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
    r=+0.950 z=+0.600     Psi:   MSa*Ia +  MSy*Iy = -0.02 +    1*Fb                 1/w=1.00e-04
                           Br:  BrSa*Ia + BrSy*Iy =     0                           1/w=0.00e+00
                           Bz:  BzSa*Ia + BzSy*Iy =     0                           1/w=0.00e+00
    r=+0.820 z=-0.450     Psi:   MSa*Ia +  MSy*Iy =     0 +    1*Fb                 1/w=1.00e-04

## Bonus: Internal profiles constraints in FBT

By default, three constraints can be set: Ip, bp and qA This is handled through `agfitfct`.

    help(func2str(L.P.agfitfct))

      function ag = meqfit3(L,rA,dr2FA,dz2FA,drzFA,fPg,fTg,Tpg,ITpg,gAg,IgAg,rBt,Ip,bp,qA)
      find ag such that the profiles match Ip,qA,bp
      Parameter L.P.agfitp = true(3,1) (default) uses all 3 constraints
      set one to false to disable it (for less BFs)

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

Three degrees of freedom in the basis functions are assumed, by default the basis function `bfab` has 1 (linear in psiN) basis function for p' and 2 basis functions a (linear+quadratic in psiN) for TT'. These are reflected in `bfp`

    help(func2str(L.bfct))
    disp(L.bfp)

     BFABMEX  Polynomial base functions for p',TT'
      BFABMEX uses P=[nP nT] and produces nP+nT base functions for p',TT' as
      the first nP<=3,nT<=3 polynomials of {F-FB, (F-FB)*(F-FA),
      (F-FB)*(F-FA)*(2*F-FB-FA)}. See also BFHELP.

      For details see: [MEQ-redbook]

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

         1     2

We now scan qA and bp

    LX=LX5;

    ii=1;
    clf;
    Ip=150e3;
    for bp = [0 1 2]
      for qA = [0.5 1 2]
        subplot(3,3,ii)
        LX.IpD = Ip;
        LX.bpD = bp;
        LX.qA = qA;
        LY(ii) = fbtt(L,LX);
        fbtplot(L,LX,LY(ii));
        legend('off');
        title(sprintf('Ip=%3.0f[kA]\nqA=%2.1f, bp=%3.1f',LY(ii).Ip/1e3,LY(ii).qA,LY(ii).bp));
        ii=ii+1;
      end
    end
    snapnow;

For a few of these, let's also plot the internal profiles

    LYcell=num2cell(LY(7:9));
    clf;
    meqplotQ(L,LYcell{:})
    snapnow;