function [Tstate,LY,alarm] = meqlim(L,LY,dt,it,Tstate)
% function [Tstate,LY,alarm] = meqlim(L,LY,dt,it,Tstate)
% Checks coil system operational limits
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if it==0
  Tstate = zeros(L.G.na,1); % init
end

%% Check coil limits

% Coil current limit
[Iarel,Ialarm] = meqIalim(L.G,LY.Ia);

% Thermal limit
[Tarel,Talarm,Tstate] = meqtherm(L.G,LY.Ia,dt,Tstate);

% Coil combination equation limit
if isfield(L.P,'eqfct')
  % optional equation or coil current limits
  [Parel,Palarm] = L.P.eqfct(L.G,LY.Ia);
else
  Parel = zeros(0,1); Palarm = false;
end

%% Alarms
if any(Ialarm) && (L.P.debug || L.P.stoponalarm)
  fprintf('\n*** Current limit exceeded on coil %s',L.G.dima{Ialarm});
  fprintf('\n');
end
if any(Palarm) && (L.P.debug || L.P.stoponalarm)
  fprintf('\n*** Protection limit equation %d violated',find(Palarm));
  fprintf('\n');
end
if any(Talarm) && (L.P.debug || L.P.stoponalarm)
  fprintf('\n*** Thermal limit on coil %d violated',L.G.dima{Talarm});
  fprintf('\n');
end

alarm = any(Talarm) || any(Ialarm) || any(Palarm);

%% Outputs
LY.Tarel = Tarel;
LY.Iarel = Iarel;
LY.Parel = Parel;
end

function [Iarel,Ialarm] = meqIalim(G,Ia)
% [Iarel,Ialarm] = meqIalim(G,Ia)
% Iarel: Relative violation of coil current limits
% Ialarm: Current limits violation

% Check current limits
Ialarm = abs(Ia) >  G.Ialim;
Iarel = Ia./G.Ialim;

end

function [Tarel,Talarm,Tstate] = meqtherm(G,Ia,dt,Tstate)
% Thermal protection implementation for coil currents
% Assume dissipated power: R*I^2.
% Thermal limit implemented as time integral of I.^2.

Tstate = Tstate + dt*Ia.^2; % time integral of current squared
Tarel = Tstate./G.Talim;
Talarm = (Tarel > 1);

end