%FBTT  FBTE equilibrium calculation
% LY = FBTT(L,LX[,'PAR',VAL,...]) calculates equilibrium using the parameters
% in L as obtained by FBTC, optionally replacing or adding parameters with
% specified values. It returns a structure LY with the calculated
% equilibria.
%
% For help on output structure: see meqt.m
% FBT-specific outputs:
% 
% .Sp    Sum of surfaces of all grid points containing plasma [m^2] 
% .Ep    Plasma box elongation    
% .chi   Cost function residual
% .chie  Equality constraints global residual
% .Fb    Boundary flux (fitting parameter)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

function LY = fbtt(L,LX,varargin)
 
 %% Arguments
 for k = 1:2:length(varargin), L.P.(varargin{k}) = varargin{k+1}; end
 
 %% Checks
 % True if bfct is bffbt (this covers also cases where bfgenD is automatically selected)
 isbffbt = strcmp(func2str(L.P.bfct),'bffbt');

 % Check sign congruence of qA,rBt,Ip
 assert(all(~any(LX.qA.*LX.rBt.*LX.IpD(1:size(LX.qA,1),:) < 0) ),'signs of qA,rBt,Ip are not congruent, product of signs must be positive')

 ctol = 1e-6; % Constraint tolerance for post-LSQ checks
 
 %% Add PFCs
 if isfield(L.P,'addpfc')
  [r,z,T] = feval(L.P.addpfc,L.P.addpfcpar);
  L.G.rw = [L.G.rw;r]; L.G.zw = [L.G.zw;z]; L.G.Twa = blkdiag(L.G.Twa,T);
 elseif isfield(L.P,'addpfcpar')
  L.G.rw = [L.G.rw;L.P.addpfcpar.r]; L.G.zw = [L.G.zw;L.P.addpfcpar.z]; L.G.Twa = blkdiag(L.G.Twa,L.P.addpfcpar.T);
 end
 na = size(L.G.Twa,2);

 %% Debug
 switch L.P.debug
   case {0,1,2}, display = 'off'  ;
   case 3,       display = 'final';
   otherwise,    display = 'iter' ;
 end
 
 assert(~isempty(which('quadprog')), 'quadprog missing, need optimization toolbox for fbtt') 
 if ~isempty(which('optimoptions')) % optimoptions not available in Octave
   osel = optimoptions('quadprog','Display',display,'Algorithm','interior-point-convex');
 else
   osel = [];
 end
 
 %% Base function parameters
 if isbffbt
   bfpt = LX.bfp;
 end

 %% Equilibrium loop
 for kt = 1:numel(LX.t)
  
  shot = LX.shot(min(kt,end));
  t    = LX.t(kt);
  rBt  = LX.rBt(kt);
  IpD  = LX.IpD(:,kt);
  qA   = LX.qA(:,kt);
  bpD  = LX.bpD(:,kt);
  capaj = LX.capaj(:,kt);
  tol = LX.tol(:,kt);
  idoublet = LX.idoublet(kt);
  isaddl = LX.isaddl(kt);
  issym = LX.issym(kt);
  niter = LX.niter(kt);
  sIp  = sign(sum(IpD));
  
  % Check input symmetry for symmetric doublets
  symdoublet = idoublet && issym; % indices with symmetric doublets
  checkfields = {'IpD','bpD','qA'};
  if any(symdoublet)
    for field = checkfields
      assert(all(~symdoublet | diff(LX.(field{:})(1:2,kt))==0),'%s must be equal if issym=1',field{:});
    end
  end
  
  %% Geometry
  rH  = LX.gpr(:,kt);
  zH  = LX.gpz(:,kt);
  kH  = ~isnan(rH) & ~isnan(zH);
  rH  = rH(kH);
  zH  = zH(kH);
  nH  = length(rH);
  bH  = logical(LX.gpb(kH,kt));
  FaH = LX.gpfa(kH,kt);
  FbH = LX.gpfb(kH,kt);
  FeH = LX.gpfe(kH,kt)*LX.gpfd(kt);
  FdH = LX.gpfd(kt);
  BrH = LX.gpbr(kH,kt);
  BzH = LX.gpbz(kH,kt);
  BaH = LX.gpba(kH,kt);
  BnH = zeros(nH,1); BnH(~isnan(BrH) & ~isnan(BzH) | isnan(BaH)) = NaN;
  BaH(isnan(BaH)) = 0;
  BeH = LX.gpbe(kH,kt)*LX.gpbd(kt);
  BdH = LX.gpbd(kt);
  CrH = LX.gpcr(kH,kt);
  CzH = LX.gpcz(kH,kt);
  CaH = LX.gpca(kH,kt);
  CnH = zeros(nH,1); CnH(~isnan(CrH) & ~isnan(CzH) | isnan(CaH)) = NaN;
  CaH(isnan(CaH)) = 0;
  CeH = LX.gpce(kH,kt)*LX.gpcd(kt);
  CdH = LX.gpcd(kt);
  VrrH = LX.gpvrr(kH,kt);
  VrzH = LX.gpvrz(kH,kt);
  VzzH = LX.gpvzz(kH,kt);
  VeH  = LX.gpve(kH,kt)*LX.gpvd(kt);
  VdH  = LX.gpvd(kt);
  
  %% PFC currents
  IaH = LX.gpia(:,kt);
  IeH = LX.gpie(:,kt)*LX.gpid(kt);
  IdH = LX.gpid(kt);
  
  %% PFC dipoles
  DaH = zeros(na-1,1);
  DeH = LX.gpde(:,kt)*LX.gpdd(kt);
  DdH = LX.gpdd(kt);

  Dax = L.G.Dda.*LX.gpdw(:,kt)'; % Assign the correct weights for indiv. coil currents
  
  %% Bounding rectangle
  if ~idoublet
    bb  =  [min(rH(bH)) max(rH(bH)) min(zH(bH))      max(zH(bH))     ]*0.5;
  else
    % assume two regions above and below axis, generate only upper box
    bb  = [[min(rH(bH)) max(rH(bH)) min(zH(bH&zH>0)) max(zH(bH&zH>0))]*0.5;
           [min(rH(bH)) max(rH(bH)) min(zH(bH&zH<0)) max(zH(bH&zH<0))]*0.5];
  end
  rh = bb(:,1)+bb(:,2); zh = bb(:,3)+bb(:,4); wh = bb(:,2)-bb(:,1); hh = bb(:,4)-bb(:,3);
  
  %% Initial plasma current distribution
  Iy = zeros(L.nzy,L.nry);
  for iD=1:numel(rh)
    IyD = max(1-((L.rry-rh(iD))/wh(iD)).^2-((L.zzy-zh(iD))/hh(iD)).^2,0).^capaj;
    IyD = IpD(iD)/sum(IyD(:))*IyD;
    Iy = Iy + IyD;
  end
  
  if issym
     Iy = (Iy+flipud(Iy))*0.5;     % up-down symmetric init
  end
  
  %% Base function parameters for this slice
  if isbffbt
    if L.nD>1, for ii = 1:size(L.bfp,1), L.bfp{ii,2} = bfpt{ii,2}(:,kt); end
    else,                                L.bfp       = bfpt(:,kt); end
  end
  %% Equilibrium dependent Green's functions
  if     any(~isnan(VrrH) | ~isnan(VrzH) | ~isnan(VzzH) & ~isnan(VeH))
   gsel = {'mut' 'br' 'bz' 'dr1r1mut' 'dr1z1mut' 'dbrdr' 'dbrdz' 'dbzdz'};
  elseif any(~isnan(CrH) | ~isnan(CzH) | ~isnan(CaH) & ~isnan(CeH))
   gsel = {'mut' 'br' 'bz' 'dr1r1mut' 'dr1z1mut' '0' '0' '0'};
  elseif any(~isnan(BrH) | ~isnan(BzH) | ~isnan(BaH) & ~isnan(BeH))
   gsel = {'mut' 'br' 'bz' '0' '0' '0' '0' '0'};
  else
   gsel = {'mut' '0' '0' '0' '0' '0' '0' '0'};
  end
  [MHy,BrHy,BzHy,CrHy,CzHy,~    ,~    ,~    ] = greenem(gsel,rH,zH,L.rry ,L.zzy         );
  [MHa,BrHa,BzHa,CrHa,CzHa,VrrHa,VrzHa,VzzHa] = greenem(gsel,rH,zH,L.G.rw,L.G.zw,L.G.Twa);
  BnHy = sin(BaH).*BrHy      + cos(BaH).*BzHy;
  BnHa = sin(BaH).*BrHa      + cos(BaH).*BzHa;
  CnHy = atan2(CzHy,CrHy);
  CnHa = atan2(CzHa,CrHa);
  CnHy = sin(CaH).*cos(CnHy) + cos(CaH).*sin(CnHy);
  CnHa = sin(CaH).*cos(CnHa) + cos(CaH).*sin(CnHa);
  
  %% Constraints and cost function
  % Cx*x = d-Cy*Iy, w.*(Ax*x ~= b-Ay*Iy), x = [Ia;Fb] (Note: x will later be rescaled).
  %     PFC current     PFC dipole           Flux            Br                  Bz                  Bn                  Cr (flux curvature) Cz                  Cn                  Vrr (vacuum dBr/dr)  Vrz                  Vzz                  
  ID = [ones(na,1)    ; repmat(5,na-1,1)   ; repmat(2,nH,1); repmat(3,nH,1)    ; repmat(3,nH,1)    ; repmat(3,nH,1)    ; repmat(4,nH,1)    ; repmat(4,nH,1)    ; repmat(4,nH,1)    ; repmat(6,nH,1)     ; repmat(6,nH,1)     ; repmat(6,nH,1)     ];
  w  = [IeH           ; DeH                ; FeH           ; BeH               ; BeH               ; BeH               ; CeH               ; CeH               ; CeH               ; VeH                ; VeH                ; VeH                ];
  Ax = [eye(na,na+1)  ; [Dax zeros(na-1,1)]; [MHa -FbH]    ; [BrHa zeros(nH,1)]; [BzHa zeros(nH,1)]; [BnHa zeros(nH,1)]; [CrHa zeros(nH,1)]; [CzHa zeros(nH,1)]; [CnHa zeros(nH,1)]; [VrrHa zeros(nH,1)]; [VrzHa zeros(nH,1)]; [VzzHa zeros(nH,1)]];
  b  = [IaH           ; DaH                ; FaH           ; BrH               ; BzH               ; BnH               ; CrH               ; CzH               ; CnH               ; VrrH               ; VrzH               ; VzzH               ];
  Ay = [zeros(na,L.ny); zeros(na-1,L.ny)   ; MHy           ; BrHy              ; BzHy              ; BnHy              ; CrHy              ; CzHy              ; CnHy              ; zeros(nH,L.ny)     ; zeros(nH,L.ny)     ; zeros(nH,L.ny)     ];
  wC = [IdH*ones(na,1); DdH*ones(na-1,1)   ; FdH*ones(nH,1); BdH*ones(nH,1)    ; BdH*ones(nH,1)    ; BdH*ones(nH,1)    ; CdH*ones(nH,1)    ; CdH*ones(nH,1)    ; CdH*ones(nH,1)    ; VdH*ones(nH,1)     ; VdH*ones(nH,1)     ; VdH*ones(nH,1)     ];
  % Remove lines with zero or NaN weights or NaN values
  k = ~isfinite(w) | isnan(b); w(k) = []; Ax(k,:) = []; b(k) = []; Ay(k,:) = []; wC(k) = []; ID(k) = [];
  % Lines with infinite weights (w == 0 now) are constraints
  k = w == 0; Cx = Ax(k,:);  d = b(k) ; Cy = Ay(k,:); wC = wC(k);
  w(k) = [];  Ax(k,:) = [];  b(k) = []; Ay(k,:) = []; ID(k) = [];
  w  = 1./w ; Ax = w .*Ax; b = w .*b; Ay = w .*Ay;
  wC = 1./wC; Cx = wC.*Cx; d = wC.*d; Cy = wC.*Cy;

  % Inequality constraints:
  % Starting from liml < limc*Ia < limu
  limu = L.P.limm.*L.P.limu; liml = L.P.limm.*L.P.liml; % Add margin factor to limits
  % limits that involve only one coil can be translated into c1 < x < c2
  l = sum(L.P.limc ~= 0,2) == 1;                  % limit on just one Ia value
  C = L.P.limc(l,:); c1 = liml(l); c2 = limu(l);
  [k,~] = find(C); cscal = 1./diag(C(k,:)); % scale to constrain x directly
  c1 = [cscal.*c1(k);-Inf]; c2 = [cscal.*c2(k); Inf];
  % limits that involve more coils are translated into C*x < c3
  C  = [[L.P.limc(~l,:); -L.P.limc(~l,:)] zeros(2*sum(~l),1)];
  c3 =  [    limu(~l  ); -    liml(~l  )];
  k = all(C==0,2) | any(isnan(C),2) | ~isfinite(c3);
  C (k,:) = [];
  c3(k  ) = [];
  
  % Now rescale x such that x = [Ia;Fb/L.ffb]
  Ax(:,end) = Ax(:,end) * L.ffb;
  Cx(:,end) = Cx(:,end) * L.ffb;
  C (:,end) = C (:,end) * L.ffb;
  c1(end)   = c1(end)   * L.ffb;
  c2(end)   = c2(end)   * L.ffb;
  Hx = (Ax'*Ax); % Hessian for quadprog
  
  %% 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')));
  Ia                           =      zeros(meqsize(L,  'Ia'));
  Iu                           =      zeros(meqsize(L,  'Iu'));
  Opy                          =      zeros(meqsize(L, 'Opy'),'int8');
  lB = false;
  lX = false;

  %% Picard iterations
  kit = 1; kip = 1; conv = false; 
  while 1
   
   %% Flux and fields from plasma current distribution
   if issym
    Iy = (Iy+flipud(Iy))*0.5; % up-down symmetric
   end
   
   %% Find x = [Ia;Fb/L.ffb]
   by = b-Ay*Iy(:);
   dy = d-Cy*Iy(:);
   fx = -Ax'*by;
   
   % Temporary rescaling aiming at x, cost function and constraints all being of order 1
   xscale = 10.^round(log10( sum(abs(fx))/sum(abs(Hx(:))) )); % unknowns scale
   fscale = 10.^round(log10( sum(abs(fx))                 )); % cost function scale
   escale = 10.^round(log10( sum(abs(Cx),2)               )); % exact constraints scale (vector)
   iscale = 10.^round(log10( sum(abs(C ),2)               )); % inequality constraints scale (vector)
   
   % catch case that all target values are zero
   xscale(xscale==0 | ~isfinite(xscale))=1;
   fscale(fscale==0 | ~isfinite(fscale))=1;
   escale(escale==0 | ~isfinite(escale))=1;
   iscale(iscale==0 | ~isfinite(iscale))=1;
   % 
   H_   = Hx./fscale.*xscale;
   f_   = fx./fscale;
   A_   = C ./iscale;
   b_   = c3./iscale./xscale;
   Aeq_ = Cx./escale;
   beq_ = dy./escale./xscale;
   lb_  = c1        ./xscale;
   ub_  = c2        ./xscale;
   %
   if kip == 1
     % Disable inequality constraints
     A_ = []; b_ = []; lb_ = []; ub_ = [];
   end
   [x_,~,stat] = quadprog(H_,f_,A_,b_,Aeq_,beq_,lb_,ub_,[],osel);
   % If problem is infeasible, quadprog can return an empty solution
   if stat <= 0 && isempty(x_), x_ = NaN(na+1,1);end
   x = x_*xscale;
   %
   ex = Ax*x-by; % residual
   
   if (stat <= 0)
     meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'Ia fit by quadprog failed','lsqlinFail');break;
   end
   
   if (norm(Cx*x-dy,Inf) > ctol)
     meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'Ia fit by quadprog does not satisfy exact constraints','lsqlinEqFail');break;
   end
   
   if kip>1 && any( (c1-ctol) > x )
     meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'Ia fit by quadprog violates lower bound','lsqlinLBFail'); break
   end
   
   if kip>1 && any( x > (c2+ctol)  )
     meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'Ia fit by quadprog violates lower bound','lsqlinUBFail'); break
   end
      
   if kip>1 && any( C*x > (c3+ctol)  )
     meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'Ia fit by quadprog violates inequality constraint','lsqlinIneqFail'); break
   end

   Ia = x(1:na);
   
   %% Update flux distribution
   Fx = meqFx(L,Iy,[Ia;Iu]);
   
   Fy = Fx(2:end-1,2:end-1); % this is 30% faster than Fx(L.lxy) apparently

   %% Plasma domain and current distribution
   [rA,zA,FA,dr2FA,dz2FA,drzFA,rX,zX,FX,~,~,~,~,~,...
     rB,zB,FB,lB,lX,Opy] = meqpdom(Fx,sIp,isaddl,L);
   nA=numel(rA);
   if ~lB,                        meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'could not find LCFS'                   ,'NoLCFS'   ); break; end
   if nA==0,                      meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'No magnetic axis'                      ,'NoMagAx'  ); break; end
   if ~idoublet&&nA>1,            meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'Multiple magnetic axes but not doublet','MultMagAx'); break; end
   if any(dr2FA.*dz2FA<drzFA.^2), meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'Tilted Hessian at magnetic axis'       ,'MagAxTilt'); break; end
   
   %% plasma current from basis function coefficients
   nB = numel(FB); % number of domains
   [F0,F1] = meqpdomlim(FA,FB,L.nD); % limiting fluxes on various domains
   
   if nB > L.nD
     meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,...
       sprintf('number of found domains (%d) exceeds maximum number (%d)',nB,L.nD),'TooManyDomains');
     break
   end
   % Fit assumes all domain fits are independent
   ag = zeros(L.ng,1); % init
   [Tyg,TpDg,ITpDg] = L.bfct(1,L.bfp,Fx,F0,F1,Opy,L.ry,L.iry);
   [g0g,Ig0g      ] = L.bfct(5,L.bfp,[],F0,F1);
   for iD=1:nA % fit ag for axis domains
     ig = logical(L.TDg(iD,:));
     ag(ig) = L.P.agfitfct(L,rA(iD),dr2FA(iD),dz2FA(iD),drzFA(iD),...
       L.fPg(ig),L.fTg(ig),TpDg(iD,ig),ITpDg(iD,ig),g0g(iD,ig),Ig0g(iD,ig),rBt,IpD(iD),bpD(iD),qA(iD)); % ag fit per domain
   end
   for iD = (nA+1):nB % possible mantle domains
     if ~any(Opy==iD), continue; end % if domain has no Opy points, skip
     ig = logical(L.TDg(iD,:));
     ag(ig) = L.P.agfitfct(L,[],[],[],[],...
       L.fPg(ig),L.fTg(ig),TpDg(iD,ig),ITpDg(iD,ig),g0g(iD,ig),Ig0g(iD,ig),rBt,IpD(iD),bpD(iD),[]); % ag fit per domain
   end
   Iy(:) = Tyg*ag;

   %% Flux residual and tolerance
   if kit>1; res= max(abs(Fy(:)-zFy(:)))/abs(FA(1)-FB(1)); else, res=NaN; end

   %% Debug
   if L.P.debug
     % display information
     nc = fprintf('%s#%d %6.4fs/%02d/%02d Ip=%.0fkA,nA=%d,res=%8.2e,tol=%8.2e,',...
       L.P.tokamak,shot,t,kit,kip,sum(IpD)*1e-3,nA,res,tol);
     for iA = 1:nA
       if iA > 1, fprintf('%*s',nc,'');end % Add spaces for alignment
       fprintf('r,z=%5.3f,%+6.3fm\n',rA(iA),zA(iA));
     end
   end
   
   if L.P.debug>2
     lim = L.P.limc*Ia;
     k = lim > 0;
     lim( k) =  lim( k)./limu( k);
     lim(~k) = -lim(~k)./liml(~k);
     % Same contour selection as in liut
     if ~lB % if meqpdom could not find a valid LCFS point
       F = 21;
     else % plot specific contour lines depending on FB,FB-FA
       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),0,sum(Opy(:) == iA));
     end
  
     clf;
     meqvis([{sprintf('%d/%d.%d',kt,kip,kit)}...
       rztext...
       {sprintf('Ip=%.0fkA',sum(IpD)*1e-3)}],...
       L.G.rx,L.G.zx,Fx,F,FB,...
       {L.G.rl([1:end 1]) rH rA rX rh+wh*[-1 1 1 -1 -1] L.ry},...
       {L.G.zl([1:end 1]) zH zA zX zh+hh*[-1 -1 1 1 -1] Iy(iround(L.zy,zA),:)/max(abs(Iy(:)))*.2+zA},...
       {'-m' 'om' '+m' 'xm' '-m' '-w'},...
       {lim(1:na) lim(na+1:end) ex(ID==1) ex(ID==2) ex(ID==3) ex(ID==4) ex(ID==5)},'ycrbkmg',...
       {'Limits: Ia, C*Ia' 'Cost fct: Ia \psi B H_{\psi} Da'},...
       '')
     drawnow;
   end
   
   %% Convergence test
   if kit>1 && res<tol
    % Skip second round with inequality constraints if they are already satisfied
    if kip == 2 || ( ~any((c1-ctol) > x) && ~any(x > (c2+ctol)) && ~any(C*x > (c3+ctol)))
     conv = true; break % Success!
    else
     kit = 1;
     kip = 2;
    end
   elseif kit == niter
    meqmsge('w',mfilename,L.P.tokamak,t,[kip,kit],shot,'iterations exceeded','nIterMax');break
   else
    kit = kit + 1;
   end
   zFy = Fy;
   
  end % while 1 (iterations over kip=[1,2], kit=[1:niter])
  
  %% Post-processing
  aq = []; aW = [];  % no initial guess as time slices are independent
  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);
  
  % fbt specific
  isconverged = conv; niter = kit;
  Sp   = sum(Opy(:)) * L.dsx;
  bb   = bboxmex(logical(Opy),L.zy,L.ry);
  Ep   = (bb(3) - bb(1)) / (bb(4) - bb(2));
  Fb   = x(end)*L.ffb;                                                % Value of the "boundary" flux (free fitting parameter)
  chi  = sqrt(sum((Ax*x + Ay*Iy(:) - b).^2)/numel(b));                % From approximate FBT constraints
  chie = sqrt(sum((Cx*x + Cy*Iy(:) - d).^2)/numel(d));                % From exact FBT constraints
  res  = chi+chie;                                                    % composite residual
  LYt  = meqlarg(LYt,Sp,Ep,isconverged,niter,chi,chie,Fb,res); % FBT specific

  if isconverged || L.P.LYall
   LY(kt) = LYt;%#ok<AGROW>
  end
  
  if L.P.debug>1
   cla; meqplott(L,LYt); drawnow;
  end
 end % kt

 if exist('LY','var') % There is at least one converged time slice
  LY = meqlpack(LY,{'FW'});
 elseif nargout > 0
  LY = struct([]);
 end
end
