function [Ia,Va,tVa,G] = voltage_stray_traces(varargin)
% function [Ia,Va,t] = voltage_stray_traces('param',value,...)%
% Returns current references, voltage feedforwards for PCS stray shot.
% Assumes all coils are in voltage mode.
%
% Ia: Currents in time [A]
% Va: Voltages in time [V]
%  t: Time base (non-equidistant) [s];
%
% ordering as mga.m: E1-8,F1-8,OH1-8
%
% Optional input parameters can be given, inspect this function for
% details.

p=inputParser;
p.addParameter('EFbias',2000, @isnumeric); % Bias current for EF
p.addParameter('EFamp', 1000, @isnumeric); % Excursion amplitude for EF
p.addParameter('OHbias',2000, @isnumeric); % Bias curent for OH
p.addParameter('OHamp', 1000, @isnumeric); % Excursion amplitude for OH
p.addParameter('ddt',   0.015,@isnumeric); % macro time grid per segment
p.addParameter('dt',    1e-3, @isnumeric); % fine time grid for solving Va;
p.addParameter('dtout', 1e-3, @isnumeric); % output resampling;
p.addParameter('dosave',false,@islogical); % save to local .mat file
p.addParameter('fname','voltage_stray_traces.mat',@ischar); % path to saved file
p.addParameter('tC',{[-0.3,3],[]},@(x) iscell(x) && numel(x)==2);
p.addParameter('selu','v',@ischar);
p.addParameter('nu',256,@isnumeric);
p.addParameter('doplot',false,@islogical);
p.addParameter('warnlimit',true,@islogical);
p.addParameter('scaling',1,@isnumeric);
p.addParameter('type','pulses',@ischar);
parse(p,varargin{:});

P = p.Results;

%% Time for each time segment
%% Order of actuation combinations
% E,F coils
EFcomb = eye(8);

% OH coils
OHcomb = eye(2);

% one trajectory
yo = [0  1  1  0]';
to = [1  2  3  4]*P.ddt;

nef = numel(yo)*size(EFcomb,1);
noh = numel(yo)*size(OHcomb,1);

% Repeat them
Etrajectory = []; tt=[];
for ii=1:size(EFcomb,1)
  Etrajectory = [Etrajectory;
    P.EFamp *yo*EFcomb(ii,:)]; %#ok<*AGROW>
  if isempty(tt)
    tt = to;
  else
    tt = [tt,tt(end)+to];
  end
end

Ftrajectory = [];
for ii=1:size(EFcomb,1)
  Ftrajectory = [Ftrajectory;
    P.EFamp *yo*EFcomb(ii,:)];
  if isempty(tt)
    tt = to;
  else
    tt = [tt,tt(end)+to];
  end
end

OHtrajectory = [];
for ii=1:size(OHcomb,1)
  OHtrajectory = [OHtrajectory;
    P.OHamp * yo*OHcomb(ii,:)];
  tt = [tt,tt(end)+to];
end

%% Construct full time sequence for Ia
% Coil indices from mga.m
na = 18; tmp = num2cell(1:na);
Ia = zeros(na,2*nef+noh);
[E1,E2,E3,E4,E5,E6,E7,E8,F1,F2,F3,F4,F5,F6,F7,F8,OH1,OH2] = deal(tmp{:}); %#ok<ASGLU>

Ia(E1:E8,1:nef) = Etrajectory';
Ia(F1:F8,(nef+1):2*nef) = Ftrajectory';
Ia(OH1:OH2,2*nef + (1:noh)) = OHtrajectory';

Ia(E1:F8,:)   = Ia( E1: F8,:) + P.EFbias;
Ia(OH1:OH2,:) = Ia(OH1:OH2,:) + P.OHbias;

%% Add stairs at start and ramp and end at 0
Ialevel = zeros(na,1); Ialevel(E1:F8) = P.EFbias; Ialevel(OH1:OH2) = P.OHbias;

Iastairs = Ialevel*[0 0 1 1  2  2  3  3  4  4  5  5 ]/5;
tstairs  =         (0:2:2*(size(Iastairs,2)-1))*2*P.ddt;

%Ia = [Iastairs,Ia,Ia(:,1),zeros(na,2)];
%tt = [ tstairs (nn+1)*P.ddt+tt tt(end)+(nn+2)*P.ddt, tt(end)+(nn+15+[1,2])*P.ddt];

%% Assemble
switch P.type
  case 'stairs'
    Ia = [Iastairs,fliplr(Iastairs),-Iastairs,-fliplr(Iastairs)];
    tl=(tstairs(end)+0.01);
    tt = [tstairs,tl+tstairs,2*tl+tstairs,3*tl+tstairs];
  case 'pulses'
    Ia = [Ia(:,1),Ia,Ia(:,1)]; tt = [0,tt+0.01,tt(end)+0.02];
    zz = zeros(na,1);
    Ia = [zz,zz,Ia,-Ia,zz,zz];
    tt = [0,0.05,tt+0.1,tt(end)+0.3+tt,2*tt(end)+[0.35,0.4]];
  otherwise
    error('unrecognized case %s',P.type);
end

%% Scaling
Ia = P.scaling*Ia;
%Ia = [Ia,-Ia];
%tt = [tt , tt(end)+0.05+tt];
%% Feedforward voltages for this current trajectory (no plasma)
sela = {'E','F','OH'};

[Va,Iv,Ia,tVa,G] = voltage_feedforward(tt,Ia,P.dt,'sela',sela,'selu',P.selu,'nu',P.nu);

% interpolate Ia,Va to final t grid
tout = tt(1):P.dtout:tt(end);
Ia = interp1(tVa,Ia',tout,'linear')';
Va = interp1(tVa,Va',tout,'linear')';

%% limits
if P.warnlimit
  iex = any(max(abs(Ia),[],2)>G.Ialim,2);
  vex = any(max(abs(Va),[],2)>G.Valim,2);
  if any(iex), warning('exceeded current limit on coil %s\n',G.dima{iex}); end
  if any(vex), warning('exceeded voltage limit on coil %s\n',G.dima{vex}); end
end

%save to .mat
if P.dosave
  toh = tVa'; %#ok<NASGU>
  efcur = Ia'; %#ok<NASGU>
  efvolt = Va'; %#ok<NASGU>
  tC = P.tC;
  save(P.fname,'toh','efcur','efvolt','tC');
  fprintf('saved %s\n',P.fname);
elseif nargout==0
  fprintf('trace not saved, call with ''dosave''=true to save\n');
end

%% optional plot
for c={'E','F','OH'}
  switch c{:}
  case 'E',  n=8; o= 0; j=1;
  case 'F',  n=8; o= 8; j=2;
  case 'OH', n=2; o=16; j=3;
  end
  subplot(3,2,2*(j-1)+1);
  plot(tout,Ia(o+(1:n),:));
  set(gca,'Ylim',G.Ialim(o+1)*[-1 1]);
  title(sprintf('%s current',c{:}))
  subplot(3,2,2*(j-1)+2);
  plot(tout,Va(o+(1:n),:));
  title(sprintf('%s voltage',c{:}))
  set(gca,'Ylim',G.Valim(o+1)*[-1 1]);
end
end
