## Contents

- Eleven FBT equilibria for Anamak
- Set some global plasma quantities
- Cost function weights/errors
- 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
- Bonus 2: Custom basis functions

## 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 = 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 by default 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.
        * If P.circuit = true, additionally solves for Iu and Va satisfying
          the time-dependent circuit equations including the plasma.

      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 Ie.

      Each constraint group has a global absolute error gp?d (a scalar) and a
      relative error gp?e (a vector with as many elements as elements in the
      group), such that the expected error for each element is gp?d*gp?e.

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

      1./(gpid*gpie) * [ Ia                    = gpia                          ]  % ID=1  Active (coil) currents
      1./(gpod*gpoe) * [ Coa*Ia                = gpco                          ]  % ID=1  Active (coil) constraints
      1./(gpdd*gpde) * [ Dda*(gpdw.*Ia)        = gpda                          ]  % ID=2  Dipoles
      1./(gpud*gpue) * [ Iu                    = gpua                          ]  % ID=3  Passive currents
      1./(gpfd*gpfe) * [ MSe*Ie   - gpfb*Fb    = gpfa  - MSy*Iy                ]  % ID=4  Flux
      1./(gpbd*gpbe) * [ BrSe*Ie               = gpbr  - BrSy*Iy               ]  % ID=5  Br
      1./(gpbd*gpbe) * [ BzSe*Ie               = gpbz  - BzSy*Iy               ]  % ID=5  Bz
      1./(gpbd*gpbe) * [ Gba*[BrSe;BzSe]*Ie    = 0     - Gba*[BrSy;BzSy]*Iy    ]  % ID=5  Field angle
      1./(gpcd*gpce) * [ MrrSe*Ie              = gpcr  - MrrSy*Iy              ]  % ID=6  Frr constraint
      1./(gpcd*gpce) * [ MrzSe*Ie              = gpcz  - MrzSy*Iy              ]  % ID=6  Frz constraint
      1./(gpcd*gpca) * [ Gca*[MrrSe;MrzSe]*Ie  = 0     - Gca*[MrrSy;MrzSy]*Iy  ]  % ID=6  Hessian angle
      1./(gpvd*gpve) * [ BrrSe*Ie              = gpvrr                         ]  % ID=7  Vacuum Brr constraint
      1./(gpcd*gpve) * [ BrzSe*Ie              = gpvrz                         ]  % ID=7  Vacuum Brz constraint
      1./(gpcd*gpve) * [ BzzSe*Ie              = gpvzz                         ]  % ID=7  Vacuum Bzz constraint
      1./(gpad*gpae) * [ Va                    = gpaa                          ]  % ID=8  Active coil voltage

      Ie is the vector of conductor currents [Ia;Iu] and M??Se, B??Se matrices
      are Green's functions linking these currents to flux/field values on the control points.

      Control point constraints gpf* gpb* gpc* gpv* as well as locations gpr, gpz
      can be supplied using FBTGP(). See FBTGP for details on the syntax

      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.

      Passive current (Iu) model:
        If P.circuit = false, then Iu is assumed to be known
        It can be prescribed via LX.gpua (time-dependent)
        In this case, the corresponding M?Su*Iu, B?Su*Iu terms go to the
        right-hand side

        If P.circuit = true, Iu is solved-self-consistently following the
        circuit equation
          Mee*Iedot + Re*Iu + Mey*Iydot = [Va;0]

      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

      See also: FBT, FBTT, FBTGP, FBTX

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

## Set some global plasma quantities

    LX0.Ip = 200e3; % plasma current
    LX0.bp = 1;     % beta poloidal
    LX0.qA =  1;     % q on axis
    LX0.rBt = 1;     % rBt

## Cost function weights/errors

    % Each constraint group has a global absolute error |gp?d| (a scalar) and a
    % relative error |gp?e| (a vector with as many elements as elements in the
    % group), such that the expected error for each element is |gp?d*gp?e|.
    %
    % The global errors should be set such that a normalized error of one is
    % acceptable. For example on anamak the gpfd error is set to
    % |0.01*L.Fx0*Ip/L.Ip0| which is 1 percent of the typical poloidal flux
    % value for the machine at this Ip value.
    %
    % Hint: when running |fbt| with |debugplot=2| the bar plot on the right
    % side should show residuals of order 1.

    % We then set the global errors for the currents

    % % global errors
    LX0.gpfd = 1e-2*L.Fx0*sum(LX0.Ip,1)/L.Ip0;  % Global errors for flux points (~ 0.008 Wb)
    LX0.gpid = 5e-2*max(abs(L.Ia0));             % Global error for Current constraints (~ 80 kA)
    LX0.gpbd = LX0.gpfd/100;                     % Global error for Field constraints (~ 8e-5 T)
    LX0.gpcd = LX0.gpfd/100;                     % Global error for Hessian constraints (~ 8e-5 Wb/m^2)

## 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,timeder)

      Appends entry to the set of gp*/g1*/g2* variables describing geometry constraints
      for FBT and their time derivatives. 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)

      timeder: time derivative order of the constraint (default=0)
        td=0/1/2 populate gp*/g1*/g2* variables respectively

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

    % check and set defaults
    LX = fbtx(L,LX);

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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

    Magnetic geometry  cost / equality constraints:
    r=+1.324 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.600 z=+0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.324 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.400 z=-0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

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.0752    0.0604

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,[],[],[],[],[],[],[],[],[],[],[],[]);
    LX = fbtx(L,LX); % check and add defaults

    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

    Magnetic geometry  cost / equality constraints:
    r=+1.324 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.600 z=+0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.324 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.400 z=-0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.550 z=+0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=0.00e+00[Wb]

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

        0.0649    0.0651

## 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,[],[],[],[],[],[],[],[], [], [], [],[]);
    LX = fbtx(L,LX); % check and add defaults

    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

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

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

        0.2000    0.2002

## 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

    LX = fbtx(L,LX); % check and add defaults
    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=0.00e+00[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

    Magnetic geometry  cost / equality constraints:
    r=+1.324 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.600 z=+0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.324 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.400 z=-0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1
    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,  [],[],[],[], [], [], [],[]);

    LX = fbtx(L,LX); % check and add defaults
    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

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

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

## 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,[],[],[],[],[],[],[],[],[],[],[],[]);
    LX = fbtx(L,LX); % check and add defaults
    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

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

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

## 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, Inf,  NaN,NaN,   0.45*pi,     0,[],[],[],[], [], [], [],[]);
    LX = fbtx(L,LX); % check and add defaults

    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

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

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

## 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,[] ,[], [] ,[] ,[] ,[] ,[],-1.2,NaN,NaN, 0, [], [], [],[]);
    LX = fbtx(L,LX); % check and add defaults

    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

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

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

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.25; kap = 1.5; % minor radius, elongation
    nt = 10;
    th = (1:nt)'/nt*2*pi+0.2;
    rr = 1    +     a0.*cos(th);
    zz = 0.02 + 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   , [], [], [],[]);
    LX = fbtx(L,LX); % check and add defaults

    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

    Magnetic geometry  cost / equality constraints:
    r=+1.169 z=+0.296        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.028 z=+0.393        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.877 z=+0.347        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.773 z=+0.176        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.755 z=-0.055        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.831 z=-0.256        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.972 z=-0.353        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.123 z=-0.307        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.227 z=-0.136        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.245 z=+0.095        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.840 z=-0.480        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=0.00e+00[Wb]
                                Br:  BrSa*Ia + BrSy*Iy = +0                                err=0.00e+00[T]
                                Bz:  BzSa*Ia + BzSy*Iy = +0                                err=0.00e+00[T]
                           d2r psi: MrrSa*Ia +MrrSy*Iy = +0                                err=0.00e+00[Wb/m2]
                          drdz psi: MrzSa*Ia +MrzSy*Iy = +0                                err=0.00e+00[Wb/m2]

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

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, [], [], [], [], [], [], [], [], [], [], [], [],+0.05,[] , 0);
    LX = fbtx(L,LX); % check and add defaults

    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

    Coil current limits:
       none

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

    Magnetic geometry  cost / equality constraints:
    r=+1.324 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=+0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=+0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.600 z=+0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.676 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.876 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.124 z=-0.380        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.324 z=-0.235        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+1.400 z=-0.000        B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=7.60e-03[Wb]
    r=+0.550 z=+0.000        B Psi:   MSa*Ia +  MSy*Iy = +0.2    +    0*Fb                 err=0.00e+00[Wb]
                             B Psi:   MSa*Ia +  MSy*Iy = +0      +    1*Fb                 err=0.00e+00[Wb]
    r=+1.000 z=+0.000       dz B0r: BrzSa*Ia           = +0.05                    err=0.00e+00[T/m]

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

    ans =

      Axes (anamak#0 0.0000s/5, r,z=1.098,+0.000m+0.0mm:571,…) with properties:

                 XLim: [0.5209 1.4725]
                 YLim: [-0.5119 0.5260]
               XScale: 'linear'
               YScale: 'linear'
        GridLineStyle: '-'
             Position: [0.5703 0.1100 0.3347 0.8150]
                Units: 'normalized'

      Use GET to show all properties

## 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
    LX = fbtx(L,LX); % check and add defaults

    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]                               err=7.95e+04[A]
    PF_002             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_003             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_004             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_005             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_006             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_007             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]
    PF_008             Current:  Ia = +0.00e+00[A]                               err=7.95e+04[A]

    Coil current  combination cost / equality constraints:
       none

    Dipole  cost / equality constraints:
       none

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

    Passive currents cost / equality constraints:
       none

    Coil voltage  cost / equality constraints:
       none

    Coil voltage limits:
       none

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

    Plasma profile constraints:
    Ip         iD=1 LX.(Ip)(01)=2e+05
    Wk         iD=1 LX.(Wk)(01)=1.9e+04
    qA         iD=1 LX.(qA)(01)=    1

## Bonus: Internal profiles constraints in FBT

In FBT, the internal profiles are fully constrained by a set of (potentially non-linear) constraints as set by the `fbtagcon` parameter. This parameter is analogous to the `agcon` parameter for FGS and FGE. As many constraints as basis functions must be provided.

    % For example the default settings will use 3 basis functions and constrain
    % the values of Ip, Wk and qA.
    disp(L.ng);disp(L.P.fbtagcon)
    % Note that the value of Wk is computed using bp and Ip, so that the user
    % only needs to provide the values of Ip, bp and qA

         3

        {'Ip'}    {'Wk'}    {'qA'}

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.Ip  = Ip;
        LX.bp  = bp;
        LX.qA = qA;
        LX = rmfield(LX,{'Wk'}); % removing it causes it to be recomputed from `Ip`,`bp` in `fbtx`
        LX = fbtx(L,LX); % check and add defaults
        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;

## Bonus 2: Custom basis functions

We can use bf3imex to specify custom basis functions:

    n=21;
    GN = [linspace(1,0,n)',linspace(1,0,n)',zeros(n,1)];
    GN(end-2:end,1) = 1; % add a little pressure pedestal
    IGN = bfprmex(GN); % integrate to get IGN
    FP = [1;0;0]; FT = [0;1;0]; % first BF applies to p', second to TT', third to none.

    % define parameters
    bfct = @bf3imex;
    bfp = struct('gNg',GN,'IgNg',IGN,'fPg',FP,'fTg',FT)';

    fbtagcon = {'Ip','Wk','ag'}; % third one will be ag=0 for the third, unused basis function
    [L,LX,LY] = fbt('ana',1,0,'bfct',bfct,'bfp',bfp,'fbtagcon',fbtagcon);

    meqplotQ(L,LY);