function Ia = fgeIaind(L,LX)
% function Ia = fgeIaind(L[,LX]);
% Compute an active coil direction to drive inductive loop voltage without
% perturbing the plasma current distribution.
% This is used in fgeinit() to compute initial conditions and Va trajectories
% for maintaining a stationary equilibrium while solving the CDE.
% It also has other uses e.g. for Ip feedback control design.
%
% Current vector Ia is determined so that it generates a constant flux of 1WB
% the region of interest, using a preferred set of coils specified via
% a specified set of OH coils. Another (possibly different) set of
% compensation coils is used to compensate the magnetic field variations.
%
% The region of interest is chosen via P.compregion and can be 'Iy' 
% (plasma current region), 'Oasx' (whole limiter region + surrounding points)
% or 'Oasx-sym' (Oasx enforcing up-down symmetry).
% The latter choice does not require knowledge of the plasma location.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

I = eye(L.G.na);

SOH = I(L.ioha,:);   % Selector for OH coils
SCO = I(L.icompa,:); % Selector for COmpensating coils
S0  = I(~(L.ioha|L.icompa),:); % Selector for unused coils

switch L.P.compregion
  % grid region in which compensation should be attempted
  case 'Iy'
    Ocx = zeros(L.nzx,L.nrx);
    assert(nargin==2,...
      'must call fgeIaind with L,LX inputs when compregion = ''Iy''')
    Ocx(L.lxy) = LX.Opy(:);
  case 'Oasx'
    Ocx = L.Oasx; 
  case 'Oasx-sym'
    Ocx = L.Oasx; 
    Ocx = (Ocx + flipud(Ocx))/2; % up-down symmetrize
  otherwise
    if isnumeric(L.P.compregion) && numel(L.P.compregion)==L.nx
      Ocx = logical(L.P.compregion); % optional manual user entry
    else
      error('unknown option %s for compregion',L.P.compregion);
    end
end
Ocx = Ocx(:)/sum(Ocx(:)); % normalize so sums to 1

% Cost function J=(Cx-d)
C = [Ocx.*L.G.Mxa;  % Make flux of all coils in region of interest equal to 1.
  L.P.wcompa*SCO;   % Penalize compensating currents
  L.P.wioha *SOH];  % Penalize OH currents

d = [Ocx.*ones(L.nx       ,1)
  zeros(size(SCO,1),1);...
  zeros(size(SOH,1),1)];

% Equality constraints Aeq*x-beq = 0
% Ocx'*Mxa(,:)*SOH*Ia = 1; % mean flux in Ocx region is equal to 1Wb.
%
% S0 *Ia = 0; % force non-use of S0 coils
% Optionally force OH coil currents to be equal;
if L.P.ohequal && (size(SOH,1)>1)
  DOH = SOH(1,:)-SOH(2:end,:); % matrix giving differences between OH coil currents
  bdoh = zeros(size(DOH,1),1);
else
  DOH = zeros(0,L.G.na); bdoh = [];
end

Aeq = [Ocx'*(L.G.Mxa(:,L.ioha)*SOH); % constrain OH coil generated flux to be equal to 1
       S0;
       DOH];
beq = [1;
       zeros(size(S0,1),1);
       bdoh];

% Solve least-squares problem with equality constraints to find complete
% coil direction.
% Parametrize ICO = w + N*v where w is a solution of the equality constraint and
% N is the null space of Aeq. Hence equality constraints are satisfied for all v.
N = null(Aeq); w = Aeq\beq;
ok = (norm(beq,Inf)==0 && norm(w,Inf)==0) ||norm(Aeq*w-beq) < (sqrt(eps)*norm(beq,Inf));
assert(ok,'Aeq*w ~= beq. Are there too many constraints?') 
Ia = w + N*( (C*N)\(d - C*w));
end
