function [psstate,Vact] = meqps(G,psstate,dt,it,Vreq)
% MEQPS - General power supply model with voltage saturations;
% alarms, and delays.
% [psstate,Vact,Ialarm] = meqps(mode,G,psstate,dt,it,Vreq,Ia)
%
% [psstate] = meqps('init',G,[],dt,[],Vreq); % initialize state
% [psstate,Vact,Ialarm] = meqps('step',G,psstate,dt,it,Vreq,Ia); % initialize state
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if ~isfield(G,'Vadelay')
  Vadelay = zeros(G.na,1); % default if none defined
else
  Vadelay = G.Vadelay;
end

if dt>0
  idelays = round(Vadelay/dt);
else
  % case called with dt=0 (one time step only, delay is not relevant)
  idelays = zeros(size(Vadelay));
end

if it==0 % init
    if nargin<5
      Vreq = zeros(G.na,1); % default init
    end
    % initialize circular delay buffer as power supply state
    psstate = repmat(Vreq,1,max(idelays)+1); % initialize buffer with initial Vreq value
    Vact = Vreq; % init
else 
    %% Saturate input voltage
    Vsat = max(min(Vreq,G.Valim),-G.Valim);
    
    %% Delays
    [Vact,psstate] = delay_buffer(idelays,Vsat,it,psstate);
    
end

end
%%

function [skd,buffer] = delay_buffer(idelays,sk,it,buffer)

nu = size(sk,1);
nbuf = size(buffer,2);
skd = zeros(nu,1);

it_now = mod(it-1,nbuf)+1; % current index in buffer
buffer(:,it_now) = sk; % fill buffer with latest input

it_delays = mod(it-idelays-1,nbuf)+1; % indices for signals in buffer
for iu = 1:numel(sk)
  skd(iu) = buffer(iu,it_delays(iu)); % delayed output from buffer
end
end

