function X = fbtgp(X,varargin)
% 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.

% initialize
if ~isfield(X,'t')
  X.t = 0; %default
end
nt = numel(X.t);
list = strcat('gp',{'r','z','b','fa','fb','fe','br','bz','ba','be','cr','cz','ca','ce','vrr','vrz','vzz','ve'});

assert(numel(varargin)==numel(list),'fbtgp expects need %d+1 arguments, found %d+1',numel(list),numel(varargin))

for var = list
  if ~isfield(X,var{1})
    X.(var{1}) = NaN(0,nt);
  end  
end

[m,n] = size(varargin{1});
for ii = 1:numel(list)
  myfield = list{ii};
  X.(myfield) = gpcat(X.(myfield),varargin{ii},m,n);
end

% add scalars if any are missing
list = strcat('gp',{'fd','bd','cd','id','dd','vd'});
for var = list
  if ~isfield(X,var{1})
    X.(var{1}) = NaN(1);
  end
end
end

function c = gpcat(a,b,m,n)
if isempty(b),      b = NaN(m,n);
elseif isscalar(b), b = repmat(b,m,n);
end
assert(isempty(a) || size(a,2)==size(b,2),'FBTE:Incompatible number of equilibria when converting geometrical parameters')
c = [a;b];
end
