%LIUT  LIUQE reconstruction
% LY = LIUT(L,LX[,'PAR',VAL,...]) performs the LIUQE reconstuction using the
% parameters in L as obtained by LIUC and the measurements in LX as obtained by
% LIUX<TOK>, optionally replacing or adding parameters with specified values. It
% returns a structure LY with the reconstructed equilibria, with fields:
%
% LIUT(...,'rst',false,...) can be used to insert a custom initial guess.
% Note that this is only allowed when LX contains a single time
% slice.Values for IyD,Ie,dz should be also provided as name-value pairs.
% Values for aq,aW,z0,Fnp can also be provided when desired. Combined with
% itert>0, one can then use LIUT as an iterator. 
%
% For a listing of generic LY output fields, see help meqt.m
%
% LIU-specific output fields:
% .dz     Vertical shift
% .dzg    Fitted vertical shift
% .chi    Fit residual for liu least-squares problem
% .chih   Fit residual for FE least-squares fit of currents
% .rst    Reset flag
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function LY = liut(L,LX,varargin)

 %% Run on a single thread
 if ~isempty(which('maxNumCompThreads'))
  nthreads = maxNumCompThreads(1);
  cleanupObj = onCleanup(@() maxNumCompThreads(nthreads));
 end

 %% Arguments
 [IyD_init,Ie_init,dz_init,aq_init,aW_init,z0_init,Fnp_init,rst_init] = deal(true);
 for k = 1:2:length(varargin)
  switch varargin{k}
   % Initialise iterator state via arguments
   case 'IyD', IyD = varargin{k+1};IyD_init = false;
   case 'Ie' , Ie  = varargin{k+1}; Ie_init = false;
   case 'dz' , dz  = varargin{k+1}; dz_init = false;
   case 'aq' , aq  = varargin{k+1}; aq_init = false;
   case 'aW' , aW  = varargin{k+1}; aW_init = false;
   case 'z0' , z0  = varargin{k+1}; z0_init = false;
   case 'Fnp', Fnp = varargin{k+1};Fnp_init = false;
   case 'rst', rst = varargin{k+1};rst_init = false;
   % Additional parameters
   otherwise, L.P.(varargin{k}) = varargin{k+1};
  end
 end
 
 % Allocate variables not provided as inputs
 if IyD_init, IyD=zeros(L.nzy,L.nry,L.nD);end
 if  Ie_init, Ie =zeros(L.ne,1);end
 if  dz_init, dz =zeros(L.ndz,1);end
 if  aq_init, aq =[];end
 if  aW_init, aW =[];end
 if  z0_init, z0 =0;end
 if Fnp_init, Fnp=[];end
 if rst_init, rst=true;end
 
 Iy = sum(IyD,3);

 %% Basic checks
 assert(~L.P.slx,'Can not run liut.m with slx=1, call `liutsim.m` if you want to run Simulink version')
 assert(rst || numel(LX.t) == 1,'Providing a custom initial condition can only work with a single LX slice')

 %% Additional output arguments
 if isfield(L.P,'argout') % liut(...,'argout',{'var[=expr]' ...})
  L.P.argout = regexp(cellstr(L.P.argout),'(?<fld>\w+)=*(?<exp>.+)*','names');
 else
  L.P.argout = [];
 end

 %% Debug options
 % 0: No output
 % 1: Text output only
 % 2: Only plot at end of time iteration or on error,
 % 3: Plot all iterations
 if L.P.debug>2
  dplt = true; dpli = true; dtxt = true; % plot and text every iteration
 elseif L.P.debug>1
  dplt = true; dpli =false; dtxt = true; % plot each time slice, text every iteration
 elseif L.P.debug==1
  dplt =false; dpli =false; dtxt = true; % no plots, text every iteration
 elseif L.P.debug<=0
  dplt =false; dpli =false; dtxt =false; % no plots or text
 else
  error('invalid debug %d\n',L.P.debug);
 end

 %% ag/dz constraints (temporary ... or not)
 if ~isfield(LX,'ag'), LX.ag = zeros(L.ng ,numel(LX.t));end
 if ~isfield(LX,'dz'), LX.dz = zeros(L.ndz,numel(LX.t));end

 
 %% Allocate variables
 [Yq,Zq]    = deal(zeros(L.nq,1)); %#ok<ASGLU>
 Wij        = zeros(L.ni,L.nj);
 Wqj        = zeros(L.nq,L.nj);
 Qcj        = zeros(L.nc,L.nj);
 
 %% Time loop
 niter = int32(1);
 nt = numel(LX.t);
 if L.P.itert, nitmax = L.P.itert;
 else,         nitmax = L.P.itera; end

 for kt = 1:nt
  if L.P.itert, tstart = tic; end
  
  %% Measurements
  shot = LX.shot(min(kt,end)); t = LX.t(kt);
  Ip = LX.Ip(kt); rBt = LX.rBt(kt);
  sIp = sign(Ip);
  Xr = [LX.Ff(:,kt); LX.Bm(:,kt)]; % magnetic measurements
  Xe = [LX.Ia(:,kt); LX.Iu(:,kt)]; % external currents
  Xi = [LX.Ft(kt)*rBt ; Ip];       % integral measurements
  Xj = [LX.ag(:,kt); LX.dz(:,kt)*L.P.fdz];
  Yr = L.Wr.*Xr;
  Ye = L.We.*Xe;
  Yi = L.Wi.*Xi;
  Yj = L.Wj.*Xj;                   % 0 for exact constraints (see liuc)
  Yd = [Yr ; Ye ; Yi];
  
  % Initial values
  aj = Xj;
  a = [Xe;aj];
  
  %% Init variables to enable post-processing of failed slices
  [TpDg,ITpDg]                 = deal(zeros(meqsize(L,'TpDg')));
  ag                           =      zeros(meqsize(L,  'ag'));
  Fx                           =      zeros(meqsize(L,  'Fx'));
  [FA,rA,zA,dr2FA,dz2FA,drzFA] = deal(zeros(meqsize(L,  'FA')));
  [FB,rB,zB]                   = deal(zeros(meqsize(L,  'FB')));
  [FX,rX,zX]                   = deal(zeros(meqsize(L,  'FX')));
  Opy                          =      zeros(meqsize(L, 'Opy'),'int8');
  lB = false;
  lX = false;
  dzg = 0;
  Mrj = zeros(L.nr,L.nj);
  if L.P.debug>0
   nA = int32(0);      %#ok<NASGU>
   ie0 = 0;            %#ok<NASGU>
   [rl,zl] = deal([]); %#ok<ASGLU>
  end
  chih = NaN;
  Zd = zeros(L.nd,1);
  
  %% contours
  if rst, aq = []; aW = []; end
  
  %% Prepare for aj constraints
  if L.nje == 0
   ZYr = Yr;
  end
  
  %% FE fit
  Jh = repmat(Ip/L.nh,L.nh,1);
  if rst
   if L.P.iterh
    [Ie,Jh,~] = ipmhmex(L.Ahd,L.Aed,L.Ahe,L.Aeh,L.Ahh,L.uAhh,Yd,Xe,Jh,sIp,L.P.iterh,max(L.P.Ipmin/L.nh,abs(Jh(1)))*L.P.tolh);
    %if ~st, meqmsge('w',mfilename,L.P.tokamak,t,0,shot,'FE fit (iterh) did not converge','iterhFail'), end
   else
    Ie = L.Aed*Yd;
    Jh = L.Ahd*Yd;
   end
   Zd = L.Wdh*Jh+L.Wde*Ie; Zd(L.kdt) = Yd(L.kdt);
   chih = sqrt(sum((Zd-Yd).^2)/L.nd); % chi of current fit
   Iy = reshape(L.Tyh*Jh,L.nzy,L.nry);
   IyD(:,:,1) = Iy;
   dz = zeros(max(1,L.ndz),1); 
   Opy = int8(logical(Iy));
  end
  % Set z0 (only used when zfct is liuz0)
  if isfield(L.P,'z0'), z0 = L.P.z0;
  elseif rst,           z0 = (L.Ipzh*Jh)/(L.Iph*Jh);end
  
  %% Picard iterations
  res = 0;
  for kit = 1:nitmax
   if rst, niter = int32(1); else, niter = niter + 1; end
    
   %% Poisson solver
   Fx = meqFx(L,IyD,Ie,dz,rst,Jh);
   
   %% Axis, LCFS flux, elongation estimator
   [rA,zA,FA,dr2FA,dz2FA,drzFA,rX,zX,FX,~,~,~,rl,zl,...
    rB,zB,FB,lB,lX,Opy] = meqpdom(Fx,Ip,L.P.isaddl,L); %#ok<ASGLU>
  
   %% No plasma
   if kit == 1 && abs(Ip) < L.P.Ipmin
    meqmsge('w',mfilename,L.P.tokamak,t,0,shot,sprintf('Ip too low: %3.2e',Ip),'LowIp')
    rst = true; break;
   end
   
   if L.P.itert == 0 && kit > L.P.iterfrz
    if kit == L.P.iterfrz + 1, meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,sprintf('No convergence after %d iterations - freezing plasma domain',L.P.iterfrz),'fixOpy');end
    Opy = zOpy; % fix the plasma domain map after some iterations to avoid oscillations
   else
    zOpy = Opy;
   end
   
   ie0 = zeros(L.nD,1);
   b0 = zeros(4,L.nD);
      
   nA = int32(numel(FA)); nB = int32(numel(FB));
   if nA==0,                      meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,'No magnetic axis'                ,'NoMagAx'      ); rst = true; debug(dplt,dtxt) , break, end
   if (nA>1&&~L.P.idoublet),      meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,'Multiple magnetic axes'          ,'MultMagAx'    ); rst = true; debug(dplt,dtxt),  break, end
   if ~lB,                        meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,'LCFS not found'                  ,'NoLCFS'       ); rst = true; debug(dplt,dtxt),  break, end
   if any(dr2FA.*dz2FA<drzFA.^2), meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,'Tilted Hessian at magnetic axis' ,'MagAxTilt'    ); rst = true; debug(dplt,dtxt),  break; end
   if (nB>L.nD||nA>L.nD),         meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,'More boundaries/axis than allowed','nB>nD||nA>nD'); rst = true; debug(dplt,dtxt),  break; end
   
   for iD=1:nA
    b0(:,iD) = bboxmex(Opy==iD,L.zy,L.ry);
    ie0(iD) = (b0(4,iD) - b0(2,iD)) / (b0(3,iD) - b0(1,iD));
   end
   
   %% Debug reset
   if rst, debug(dpli,dtxt,false);  end
   
   %% Convergence test
   if L.P.itert == 0
    if kit > 1
     zFx = zFx - Fx;
     res = max(reshape(abs(zFx(2:end-1,2:end-1)),[],1))/abs(FB(1) - FA(1));
    end
    if kit > 4 && res < L.P.psichco  % successful convergence
     debug(dpli,dtxt,false); break
    elseif kit == nitmax
     meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,'iterations exceeded','nIterMax')
     rst = true; debug(dplt,dtxt,false); break
    end
    zFx = Fx;
   end
   
   %% Base functions, regularization and constraint matrix
   [F0,F1] = meqpdomlim(FA,FB,L.nD);
   r0 = [rA;zeros(L.nD-nA,1)]; ir0 = [1./rA;zeros(L.nD-nA,1)];
   [Tyg,TpDg,ITpDg] = L.bfct(1,L.bfp,Fx,F0,F1,Opy,L.ry,L.iry);
   [Qqg,Xq        ] = L.bfct(6,L.bfp,[],F0,F1,r0,ir0,L.idsx );
   [Qcg,Xc        ] = L.bfct(7,L.bfp,[],F0,F1               );
   
   %% Response matrices
   % Adapt current distribution's initial guess for dz fitting
   if L.ndz
    if rst && L.P.itert == 0
     if nB > 1,     IyD(:) = 0;
     else,          IyD(:,:,1) = reshape(bfabmex(1,[0 1],Fx,F0(1),F1(1),Opy,L.ry,L.iry),L.nzy,L.nry);
                    IyD = IyD*(Ip/sum(IyD(:)));
     end
    elseif rst
                    IyD(:,:, 1) = Iy.*(Opy <=  1);    % Do not modify Iy when only one domain (ensures backw. compatibility)
     for iD=2:L.nD, IyD(:,:,iD) = Iy.*(Opy == iD);end
    end
             Mrj = respmex(L.Mry,Tyg,L.Mdzry,reshape(IyD,L.ny,L.nD));
   else,     Mrj =         L.Mry*Tyg                                ;
   end
   Wrj = L.Wr.*Mrj;
   Wij(L.kit,1:L.ng) = L.Wt*(2e-7*sum(ITpDg,1).*L.fTg.'+L.Mty*Tyg*rBt); % Ft
   Wij(L.kip,1:L.ng) = L.Wp*sum(TpDg,1);                                % Ip
   
   Wq = L.Wq;
   % Adapt regularisation with elongation - domain specific
   if L.nD == 1 % Deactivate adaptive regularization for multiple domains (temp.)
    Wq = Wq.*min(exp(-(1./reshape(repmat(ie0, 1, L.nq/L.nD)',L.nq,1) -L.P.elomin).*L.P.wregadapt),1);
   end
   Wqj(:,1:L.ng) = Wq.*Qqg;
   Yq = Wq.*Xq;
   Qcj(:,1:L.ng) = Qcg;
   
   % Update of Yk for regularization
   Yk  = [Yi ;Yq ;  Yj ];
   Wkj = [Wij;Wqj;L.Wjj];
   
   %% aj constraints (so far only of the form aj=Xj with weights, no general linear constraints)
   kj = L.kj;
   if L.nje > 0 || (rst && L.P.itert == 0) || nB < L.nD % L.nje = number of exact constraints
    aj(kj) = 0; % This also sets coefficients from unused domains to 0
    if (rst && L.P.itert == 0 && nB>1), kj(L.ng+1:end) = false; end % All dzs removed
    if (nB < L.nD),                     kj(~any(L.TDg(1:nB,:),1)) = false;
                                        kj(L.ng+nB+1:end) = false; end % Remove all fitting coefficients for missing domains
    % Adjust A and b for exact constraints
    ZYr = Yr - Wrj*aj; ZWrj = Wrj(:,kj);
    ZYk = Yk - Wkj*aj; ZWkj = Wkj(:,kj);
    ZXc = Xc - Qcj*aj; ZQcj = Qcj(:,kj);
   else
    ZWrj = Wrj;
    ZWkj = Wkj;
    ZQcj = Qcj;
    ZYk  = Yk;
    ZXc  = Xc;
   end
   nja = sum(kj);
   
   %% Fit sources to measurements
   if L.P.iters
    if L.P.ipm
     [Ie,Zaj,~] = ipmj(ZWrj,ZWkj,L.Wre,L.Aer,L.Aee,ZYr,Ye,ZYk,Ie,[],ZQcj,ZXc,[],[],1e-6,L.P.iters);
    else
     if L.P.itert
      ZiAjj = iatamex(ZWrj,ZWkj);
     else
      ZiAjj = [ZWrj ; ZWkj]\([ZWrj ; ZWkj]'\eye(nja));
     end
     ZAjr = ZiAjj * ZWrj';
     [Zaj,Ie] = bslvmex(ZAjr*ZYr+ZiAjj*(ZWkj'*ZYk), L.Aer*ZYr+L.Aee*Ye, Xe, ZAjr*L.Wre,L.Aer*ZWrj, L.P.iters);
    end
   else
    % direct non-iterative solver
    A = [L.Wre            ZWrj            ; ...
         L.Wee            zeros(L.ne,nja) ; ...
         zeros(L.nk,L.ne) ZWkj            ];
    % if isfield(LX,'ro')
    %  AdG = [AdG ; L.wo.*[L.Moa L.Mos L.Moy*Tyg L.dzMoy*Iy(:)]];
    %  if isempty(LX.po)
    %   Yd = [Yd ; zeros(L.no,1)];
    %  else
    %   Yd = [Yd ; L.wo.*(FA + LX.po.^2*FBA)];
    %  end
    % end
    if L.P.robust
     a = robustfit(A,[ZYr ; Ye ; ZYk],[],[],'off');
    elseif L.P.ipm
     [a,~,~,~,k] = ipm(A'*A,-A'*[ZYr ; Ye ; ZYk],[zeros(L.nc,L.ne) ZQcj],ZXc,[],[],a,[],[],[],2);
     if k > 50, meqmsge('w',mfilaname,L.P.tokamak,t,it,shot,'IPM no convergence','IPMFail'), end
    else
     a = A\[ZYr ; Ye ; ZYk];
    end
    Ie = a(1:L.ne);
    Zaj = a(L.ne+1:end);
   end
   aj(kj) = Zaj;
   Zd = [L.Wre*Ie+Wrj*aj ; L.Wee*Ie ; Wij*aj];
   Zq =  Wqj*aj; %#ok<NASGU>
   ag = aj(1:L.ng);
   if L.ndz, dzg = L.P.fdz*aj(L.ng+1:end);
   else,     dzg = 0; end
   IyD = reshape(Tyg*(ag.*L.TDg.'),L.nzy,L.nry,L.nD);
   Iy  = sum(IyD,3);
   
   %% Vertical position feedback
   dz = L.P.zfct(L.P,rst,z0,zA,dzg);
   
   %% Check for -ve p (or T2)
   if L.P.itert && strcmp(func2str(L.bfct),'bf3pmex') && ag(1) < 0
    meqmsge('w',mfilename,L.P.tokamak,t,kit,shot,'negative pressure','NegPress')
    rst = true; debug(dpli,dtxt); break
   end
   
   %% Debug every iteration if dpli
   if kit~=nitmax, debug(dpli,dtxt,false); end
   rst = false;
  end % for kit
  
  if ~rst || L.P.LYall
   %% Fitted signals and residual
   Ia = Ie(1:L.G.na);
   Iu = Ie(L.G.na+1:end);
   chi= sqrt(sum((Yd-Zd).^2)/L.nd);
   if any(L.Wr==0)
    % Some sensors were deselected: cannot use Zd for their reconstructed value since it would yield NaNs.
    % Note: Computing all sensors is actually quicker than all indexing operations required otherwise
    Zr = L.Mre*Ie+Mrj*aj;
    Ff = Zr(L.kdf);
    Bm = Zr(L.kdm);
   else
    Zd = Zd./L.Wd;
    Ff = Zd(L.kdf);
    Bm = Zd(L.kdm);
   end
   
   %% Post-processing
   LYt = struct('shot',shot,'t',t,'aq',aq,'aW',aW);
   [LYt] = meqpost(L,LYt,TpDg,ITpDg,ag,...
    Fx,FA,FB,rA,zA,dr2FA,dz2FA,drzFA,rB,zB,lB,lX,rX,zX,FX,...
    rBt,Ia,Iu,Iy,Opy);
   
   % override measurement fits with LIUQE versions which include dz term
   LYt.Ff = Ff; LYt.Bm = Bm;
   
   % Loop voltage at interpolation points from finite differences
   if L.nn
    if kt==1 || rst || isempty(Fnp), Vn = zeros(L.nn,1);
    else, Vn = (LYt.Fn-Fnp)/(LX.t(kt)-LX.t(kt-1));
    end
    Fnp = LYt.Fn;
    LYt.Vn = Vn;
   end
   
   % Additional liuqe-specific outputs
   isconverged = ~rst;
   LYt = meqlarg(LYt,...
    dz,dzg,niter,rst,isconverged);
   
   if L.P.itert == 0
    LYt = meqlarg(LYt,chi,chih);
   else
    LYt.cycle = toc(tstart)*1e6; % cycle time in us
   end
   for k = 1:length(L.P.argout)
    % additional output arguments can be requested via L.P.argout{:}
    % each cell contains .exp and .fld such that LY.(fld) = eval(exp)
    if isempty(L.P.argout{k}.exp), L.P.argout{k}.exp = L.P.argout{k}.fld; end
    LYt.(L.P.argout{k}.fld) = eval(L.P.argout{k}.exp);
   end
   % Update aq, aW
   aq = LYt.aq; aW = LYt.aW;
   % Add to LY array 
   LY(kt) = LYt; %#ok<AGROW>
  end
  % debug plot including post-processing
  dpost = true; debug(dplt,dtxt,dpost);
  % Prepare FE fit for next iteration if itert=0
  rst = rst || L.P.itert == 0;
 end % for kt
 
 if (~exist('LY','var')), LY = struct([]);end
 LY = meqlpack(LY,{'ppQg','TTpQg','FW'});

 end

 %% Debug
 function debug(doplot,dotext,dopost)
 if nargin==2, dopost=false; end
 if ~(doplot || dotext), return; end
 % List of necessary variables from liut workspace, avoids shared variables and static workspaces
 vars = evalin('caller','{L,LX,rX,zX,rB,zB,rA,zA,nA,rl,zl,rst,Iy,Fx,Opy,rBt,sIp,Yd,Zd,lB,FB,FA,shot,t,niter,dz,Yq,Zq,res,ie0}');
 [L,LX,rX,zX,rB,zB,rA,zA,nA,rl,zl,rst,Iy,Fx,Opy,rBt,sIp,Yd,Zd,lB,FB,FA,shot,t,niter,dz,Yq,Zq,res,ie0] = deal(vars{:});
 
 % Optional variables
 if dopost && L.P.iterq && ~rst
  LYt = evalin('caller','LYt');
  rq = LYt.rq; zq = LYt.zq; 
  if L.nW
   aW=LYt.aW;
  end
 end
 
 if ~rst
  vars = evalin('caller','{TpDg,ITpDg,ag,Tyg}');
  [TpDg,ITpDg,ag,Tyg] = deal(vars{:});
 end
 
 nA = min(nA,L.nD);
 
 % Symbols-lines to plot on the cross section
 mvis = {'ok' '-m' '-w' '-oc' 'xm' 'dm' 'om' '.c'   'oc' '-y'              'oy'};
 rvis = {[]   []   []   []    rX   rB   rA   L.G.rl rl   L.bh([1 3 3 1 1]) L.rh};
 zvis = {[]   []   []   []    zX   zB   zA   L.G.zl zl   L.bh([2 2 4 4 2]) L.zh};
 if isfield(LX,'ro') % Optical measurements
  rvis{1} = LX.ro;
  zvis{1} = LX.zo;
 end
 % Flux contours, gaps
 if dopost && L.P.iterq && ~rst
  rvis{2} = rq([1:end 1],:);
  zvis{2} = zq([1:end 1],:);
  if L.nW
   rvis{4} = [L.G.rW L.G.rW-aW.*cos(L.G.oW)]';
   zvis{4} = [L.G.zW L.G.zW+aW.*sin(L.G.oW)]';
  end
 end
 if rst
  % on first iteration plot the contour lines of Iy
  ip = sum(Iy(:)); bp = 0; li = 0;
  [rvis{3},zvis{3}] = mycontourc(L.ry,L.zy,abs(Iy),linspace(.1,.95,5)*max(abs(Iy(:))));
 else
  [ip,~,~,~,~,~,~,~,~,bp,~,~,li] = meqint(L.fPg,L.fTg,TpDg,ITpDg,ag,Fx,logical(Opy),L,rBt);
  rvis{3} = L.ry;
  zvis{3} = NaN(L.nry,(1+L.ng)*nA);
  for iA = 1:nA
   kzA = round((zA(iA)-L.zy(1))*L.idzx+1);
   ic = (1+L.ng)*(double(iA)-1)+(1:1+L.ng);
   ig = find(L.TDg(iA,:)); % Only ag from selected domain (avoids lines of 0s)
   zvis{3}(:,ic([1;1+ig(:)])) = [Tyg(kzA+(0:L.nry-1)*L.nzy,ig)*ag(ig) Tyg(kzA+(0:L.nry-1)*L.nzy,ig).*ag(ig)']; % profiles and current values on z=zA
   zvis{3}(:,ic) = zvis{3}(:,ic)*(0.2*(L.G.zx(end)-L.G.zx(1))/max(sIp*zvis{3}(:,ic(1))))+zA(iA);
  end
 end
 % Measurement residuals
 Ed = Yd - Zd;
 chi = sqrt(sum(Ed.^2)/L.nd);
 
 % display information
 dz(end+1:nA) = 0;
 nc = fprintf('%s#%d %6.4fs/%02d Chi=%.1e res=%.1e Ip=%.0fkA bp=%.2f li=%.2f ',...
  L.P.tokamak,shot,t,niter,chi,res,ip*1e-3,bp,li);
 for iA = 1:nA
  if iA > 1, fprintf('%*s',nc,'');end % Add spaces for alignment
  fprintf('r,z=%5.3f,%+6.3fm%+05.1fmm k=%.2f\n',rA(iA),zA(iA),dz(iA)*1000,1./ie0(iA));
 end
 
 % optional plots
 if doplot
  if ~lB % if meqpdom could not find a valid LCFS point
   F = 21;
  else % plot specific contour lines depending on FB,FBA
   FA_ = sIp*max(sIp*FA);
   FB_ = sIp*min(sIp*FB);
   FBA = (FB_ - FA_);
   F = FB_-[linspace(-1,0,11), L.fq(end:-1:1),linspace(1,2.5,16) ]*FBA;
  end
  
  rztext = cell(1,nA);
  for iA = 1:nA
   rztext{iA} = sprintf('r,z=%5.3f,%+6.3fm%+4.1fmm:%d',rA(iA),zA(iA),dz(iA)*1000,sum(Opy(:) == iA));
  end
  
  clf;
  meqvis([{sprintf('%s#%d %6.4fs/%d',L.P.tokamak,shot,t,niter)} ...
           rztext  ...
           {sprintf('Ip=%.0fkA bp=%.2f li=%.2f',ip*1e-3,bp,li)}], ...
         L.G.rx,L.G.zx,Fx,F,FB,...
         rvis,zvis,mvis,...
         {Ed(L.kdf) Ed(L.kdm) Ed(L.kda) Ed(L.kdu) Ed(L.kdt) Ed(L.kdp) Yq-Zq},...
         'rbrbrbc',...
         {'Residuals (fmastpq)' sprintf('X=%.1e',chi) ...
          sprintf('Convergence=%.1e',res)},'')
  drawnow;
 end
 end

