function [Va,Iv,Ia,tf,G] = voltage_feedforward(t,Ia,dt,varargin)
% Computes voltage feedforward for a given t,Ia trajectory (assumed
% piecewise linear;
% optionally take vessel into account

p = inputParser;
p.addParameter('sela',{'E','F','OH'},@(x) iscell(x))
p.addParameter('selu','v',@(x)ischar(x));
p.addParameter('nu',mdsvalue('static("n_v")'),@(x) isscalar(x)&&isnumeric(x) );
p.addParameter('doplot',false)
parse(p,varargin{:});
P=p.Results;
selu = P.selu; sela = P.sela; nu = P.nu;
%%

G = get_vacuum_static(sela,selu,nu);

tf = t(1):dt:t(end); % fine time grid
Va = zeros(G.na,numel(tf)); % init size, number of intervals *2
Iu = zeros(G.nu,numel(tf));
Ia = interp1(t,Ia',tf)';
% vessel currents
Q1 = (dt*diag(G.Ru)+G.Muu) \  G.Muu;
Q2 = (dt*diag(G.Ru)+G.Muu) \ -G.Mua;
for kt=2:numel(tf)
  Iu(:,kt) = Q1*Iu(:,kt-1) + Q2*(Ia(:,kt)-Ia(:,kt-1));
end

% inductive + resistive coil voltage
for kt=2:numel(tf)
  Va(:,kt)  = diag(G.Ra)*Ia(:,kt) + (G.Mau*(Iu(:,kt)-Iu(:,kt-1)) + G.Maa*(Ia(:,kt)-Ia(:,kt-1)))/dt;
end

Iv = G.Tvu*Iu;

if P.doplot
  subplot(311)
  plot(tf,Ia);
  title('Ia');
  subplot(312);
  plot(tf,Va);
  title('Va');
  subplot(313);
  plot(tf,sum(Iv));
  title('sum(Iv)')
end
