%IPMJ  Interior point method for base function fitting, under development
% [IE,AJ,ST] = IPMJ(WRJ,WIJ,WRE,AER,AEE,YR,YE,YI,IE0,AJ0,QCJ,XC,S,Z,TOL,NITS)
% solves for IE,AJ in a LS sense and imposing QCJ*AJ>XC
%       YR = [WRE WRJ]*[IE;AJ]
%  YD = YE = [WEE   0]*[IE;AJ]
%       YI = [  0 WIJ]*[IE;AJ]
% by iteratively solving for AJ and IE respectively
%  YR-WRE*IE = WRJ*AJ <=> min AJ'*AJJ*AJ/2+(-WRJ'*YR-WIJ'*YI+WRJ'*WRE*IE)'*AJ 
%  YI        = WIJ*AJ <=> QCJ*AJ>XC
% and
%  YR-WRJ*AJ = WRE*IE <=> IE=AER*YR+AEE*YE-AER*WRJ*AJ
%  YE        = WEE*IE
% AJJ=WRJ'*WRJ+WIJ'*WIJ, [AER;AEE]=[WRE'*WRE+WEE'*WEE]\[WRE' WEE']. IE0,AJ0
% are initial guesses. NITS,TOL control the iteration numbers and the
% convergence criterion. ST=1 if success.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function [Ie,aj,st] = ipmj(Wrj,Wkj,Wre,Aer,Aee,Yr,Ye,Yk,Ie,aj,Qcj,Xc,s,z,tol,iters)
 uAjj = uatamex(Wrj) + uatamex(Wkj);
 Aje = Wrj'*Wre;
 Aej = -Aer*Wrj;
 Ie0 = Aer*Yr+Aee*Ye;
 aj0 = -Wrj'*Yr-Wkj'*Yk;
 if isempty(aj), aj = double(any(Qcj,1)).'; end
 if isempty(s ), [s,z] = deal(repmat(tol,size(Xc))); end
 [Ie,aj,st] = ipmjmex(uAjj,Aje,Aej,Ie0,aj0,Ie,aj,Qcj,Xc,s,z,tol,iters);
end