function [LX] = meqinit(L,LX)
% function [LX] = meqinit(L,LX)
% Compute Ia, Iv, Va consistent with initial equilibrium for forward codes, 
% also in the case of a time-varying equilibrium.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(numel(LX.t)==1,'only one time slice for LX in fgeinit');

if L.P.ssinit % stationary-state initial condition
  if any(L.icde) && LX.Ip~=0
    LX = meqinitcde(L,LX);
  else
    LX = meqinitIp(L,LX);
  end
else
  % non-stationary-state init, just pass default Va 
  % to sustain initial present Ia value.
  [LX.Vadot0,LX.Iadot0] = deal(zeros(L.G.na,1));
  LX.Va = L.G.Ra.*LX.Ia;
end
end

function LX = meqinitcde(L,LX)
% Plasma current needs to be induced by a combination of ramping coil currents.
% The nominal simulation features time-varying Ia=Ia0+Iadot0*(t-t(1)); Iv~=0;
% This function computes the additional compensating Ia to counteract Iv fields
% adding: 
%  * The corrected LX.Ia, and corresponding LX.Va.
%  * The non-zero LX.Iv, LX.Iu for this loop voltage 
%  * LX.Iadot0, LX.Vadot0 for the coil and voltage time derivatives

ia = 1:L.G.na; iu = L.G.na+(1:L.G.nu);
% Coil direction for inductive current drive with minimal change in fields.
Iaind = fgeIaind(L,LX);

% Scale the coil direction to obtain desired vloop
Ix = zeros(L.nzx,L.nrx); Ix(L.lxy) = LX.Iy;
Fpa = dot(Ix(:),(L.G.Mxa*Iaind))/sum(Ix(:)); % plasma current-weighted flux induced by Iaind

switch L.P.cde
 % Stationary condition depending on type of cde
 case 'cde0Dss'
  % Determine required loop voltage
  iIp = (abs(LX.IpD)>0); % domains with current
  Vloops = -(LX.IniD(iIp)-LX.IpD(iIp))./(LX.SQ(end,iIp))'; % mean vloop for all domains
  Vloop0 = mean(Vloops); % take mean Vloop - assumes doublets are symmetric
  Iadot0 = -(Vloop0/Fpa)*Iaind; % Ia ramp rate to sustain this Vloop
 case {'OhmTor', 'OhmTor_rigid'}
  % In ss condition Rp Ip = - Mpa Iadot0 , where by definition Mpa =  (Iy' Mya)/Ip
  Fpa = LX.Iy(:)'*(L.G.Mxa(L.lxy,:)*Iaind)/LX.IpD;
  Iadot0 = -(LX.Rp*LX.IpD/Fpa)*Iaind;
end
% Compute stationary vessel currents
Iu0 = - L.Mee(iu,ia)*Iadot0 ./ L.G.Ru;

% Compute compensation coil modification to compensate for this extra Iu;
dIe = zeros(L.ne,1); dIe(iu) = Iu0-LX.Iu(:,1);
dIacomp = fgeIacomp(L,dIe);
Ia0 = LX.Ia(:,1) + dIacomp;

%% Compute feed-forward Va(t) to achieve Ia(t) from coil circuit equation
Va0 = L.Mee(ia,ia)*Iadot0 + L.Re(ia).*Ia0;
Vadot0 = L.Re(ia).*Iadot0;

LX.Ia = Ia0;
LX.Iu = Iu0;
LX.Iv = L.G.Tvu*Iu0;
LX.Va = Va0;
LX.Iadot0 = Iadot0;
LX.Vadot0 = Vadot0;

if L.P.debug
  fprintf('\n FGE initial ss condition for CDE (fgeinit)\n')
  switch L.P.cde
    case 'cde0Dss'
      fprintf(' Ip=%2.2f[kA], Ini=%2.2f[kA], Vloop=%2.2f[V], sum(Iv)=%2.2f[kA] \n',...
        LX.Ip/1e3, LX.Ini/1e3, Vloop0, sum(LX.Iv)/1e3);
      if L.nD>1
        fprintf('Individual Vloop diff per node:\n');
        fprintf('%2.2e [V]\n',Vloops'-Vloop0);
      end
    case {'OhmTor', 'OhmTor_rigid'}
      fprintf(' Ip=%2.2f[kA], Ini=%2.2f[kA], Mpa dIadt0=%2.2f[V], sum(Iv)=%2.2f[kA] \n',...
        LX.Ip/1e3, LX.Ini/1e3, -LX.Rp*LX.IpD, sum(LX.Iv)/1e3);
  end
  fprintf(' Coil currents were modified to compensate for stray fields\n');
end
end

function LX = meqinitIp(L,LX)
% With constant Ip and no current diffusion, no vloop is necessary, so there are no vessel currents.
LX.Iu = zeros(L.G.nu, 1);
LX.Iv = zeros(L.G.nv,1);

if L.P.debug
  fprintf('\n FGE initial condition (fgeinit)\n')
  if any(LX.Iu) || any(LX.Iv)
    fprintf('No CDE and constant Ip - changed Iu=Iv=0 for stationary-state initialization\n');
  else
    fprintf('No CDE and constant Ip and Iv=0, no change needed\n');
  end
end
LX.Va = L.G.Ra.*LX.Ia;
LX.Vadot0 = zeros(L.G.na,1);
LX.Iadot0 = zeros(L.G.na,1);
end