function [LY] = meqpostq(L,LY)
%MEQPOSTQ Updates LY with quantities related to contour integrals and wall gaps
% function [LY] = meqpostq(L,LY)
% Restarts from initial guess of LY.aq, LY.aW (from previous iteration)
% 
% For details of gap/contour finding algorithm, see [MEQ-redbook]
% 
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if isempty(LY), return; end

aq = LY.aq; aW = LY.aW;
rA = LY.rA; zA = LY.zA; rB = LY.rB; zB = LY.zB; lX = LY.lX;
FA = LY.FA; FB = LY.FB; Fx = LY.Fx; Opy = LY.Opy;
Ip = LY.Ip; IpD = LY.IpD; sq = sign(Ip)*sign(LY.rBt);
FX = LY.FX; iTQ=LY.iTQ;
dr2FA = LY.dr2FA; dz2FA = LY.dz2FA; drzFA = LY.drzFA;

[rS,zS] = deal(zeros(L.noq,L.nS,L.nD));
[lp,rbary] = deal(zeros(L.nD,1));
[rgeom,zgeom,aminor,epsilon,kappa,delta,deltal,deltau] = deal(zeros(L.nQ,L.nD));
[rrmax,zrmax,rrmin,zrmin,rzmax,zzmax,rzmin,zzmin]      = deal(zeros(L.nQ,L.nD));
% initial guess
if isempty(aq)
  for iD = 1:LY.nA
    aq(:,:,iD) = repmat(min([rA(iD)-L.G.rx(1) L.G.rx(end)-rA(iD) zA(iD)-L.G.zx(1) L.G.zx(end)-zA(iD)])*...
      L.pinit*L.pq,L.noq,1);
  end
end
if isempty(aW)
  aW  = zeros(L.G.nW,L.P.nFW);
end

% Contour locations
for iD = 1:LY.nA
  F = FB(iD)-L.fq*(FB(iD)-FA(iD)); % flux surfaces to be tracked
  if L.P.icsint
    aqD = rtcics(Fx,F,rA(iD),zA(iD),L.crq,L.czq,L);
  else
    ero = (rA(iD)-L.G.rx(1))*L.idrx;
    ezo = (zA(iD)-L.G.zx(1))*L.idzx;
    aqD = aq(:,:,iD);
    for k = 1:L.P.iterq
      [aqD,daqD] = rtcimex(aqD,ero,ezo,Fx,F,L.cdrq,L.cdzq,Opy,FA(iD),int8(iD),L.drx);
      if ~L.liurtemu
        if all(daqD < L.drx*1e-3), break; % convergence tests only for non-itert
        elseif k==L.P.iterq
          meqmsge('w',mfilename,L.P.tokamak,LY.t,k,LY.shot,'meqpostq contour finding did not converge','iterq');
        end
      end
    end
  end
  aq(:,:,iD) = aqD;
end
rq  = reshape(rA,[1,1,LY.nA]) + L.crq.*aq;
zq  = reshape(zA,[1,1,LY.nA]) + L.czq.*aq;
irq = 1./rq;

% Wall gaps
% Flux values to track
FW = meqFW(FB(LY.nB),FX,lX(LY.nB),LY.nX,L.P.nFW);
if L.nW
  if L.P.icsint
    aW = rtcics(Fx,FW,L.G.rW,L.G.zW,-cos(L.G.oW),sin(L.G.oW),L);
  else
    % Values at origin points
    FoW = bintmex(Fx,int32(L.kxW-1),L.cWx);
    % default gaps, then secondary x-point gaps
    for k=1:L.P.iterq
      [aW,daW ] = rtcimex(aW,L.erW,L.ezW,Fx,FW,L.cdrW,L.cdzW,Opy,FoW,int8(0),L.drx);
      aW = min(aW,L.G.aW); % clamp aW
      if ~L.liurtemu
        if max(daW) < L.drx*1e-3, break, end % convergence tests only for non-itert
      end
    end
  end
end

% FS metrics
[Q0Q,Q1Q,Q2Q,Q3Q,Q4Q,iqQ] = deal(zeros(L.nQ,L.nD));
for iD=1:LY.nA
  qaq = aq(:,:,iD).^2;
  M1q = qaq*L.M1q; % -da^2/qdq
  M2q = ((L.M2q*qaq).^2./qaq+qaq)./M1q; % ((da^2/2do)^2/a^2+a^2)/(-da^2/qdq)
  
  % contour integrals
  [Q0Q(:,iD),Q1Q(:,iD),Q2Q(:,iD),Q3Q(:,iD),Q4Q(:,iD),iqQ(:,iD)] = ...
    fsgimex(M1q,M2q,rq(:,:,iD),irq(:,:,iD),rA(iD),FA(iD),FB(iD),sign(IpD(iD))*sqrt(dr2FA(iD)*dz2FA(iD)-drzFA(iD).^2),...
    lX(iD),rB(iD),iTQ(:,iD),1/L.doq);
  
  % shape profiles
  [rgeom(:,iD),zgeom(:,iD),aminor(:,iD),epsilon(:,iD),kappa(:,iD),delta(:,iD),deltal(:,iD),deltau(:,iD),...
   rrmax(:,iD),zrmax(:,iD),rrmin(:,iD),zrmin(:,iD),rzmax(:,iD),zzmax(:,iD),rzmin(:,iD),zzmin(:,iD)] = ...
    shapmex(rq(:,:,iD),zq(:,:,iD),rB(iD),zB(iD),rA(iD),zA(iD),dr2FA(iD),dz2FA(iD),L.drx,L.dzx);
  
  % temporary calculations of other flux-surface-dependent quantities pending their implementation in fsgimex
  dqD = sqrt(diff([rq(:,end,iD);rq(1,end,iD)]).^2+diff([zq(:,end,iD);zq(1,end,iD)]).^2);
  drq = 0.5*(rq(:,end,iD)+[rq(2:end,end,iD);rq(1,end,iD)]).*dqD;
  lp(:,iD) = sum(dqD)'; % circumference
  rbary(:,iD) = sum(drq)'./lp(:,iD); % contour barycenter
  
  % outboard radius and q surface location
  raQ = [0;aq(L.crq==1,:,iD)']; raQ=raQ/raQ(end); % outboard minor radius
  [raqmin,iqmin] = minQmex(raQ,iqQ(:,iD),-sq,1); qmin = 1./iqmin;
  if L.nR
    raR = locRmex(raQ,iqQ,sq*L.P.iqR,L.P.naR,L.raN); % radius of desired 1/q surfaces
  end
  
  % coordinates of flux surfaces specified as r/a
  if L.nS
    aS = locSmex([zeros(L.noq,1),aq],raQ,L.P.raS,L.raN); % interpolated flux values at this r/a location
    rS(:,:,iD) = L.crq.*aS + rA(iD); zS(:,:,iD) = L.czq.*aS + zA(iD);
  end
end

% q95
q95 = 1./(L.c95*iqQ(L.i95,:))';

% jtorQ = R0*<jphi/R>
jtorQ = 2*pi*L.P.r0 * (LY.PpQ +( LY.TTpQ.*Q2Q)/(4e-7*pi));

% collect new additions in LY
LY = meqlarg(LY,iqQ,jtorQ,q95,lp,rbary,...
  rgeom,zgeom,aminor,epsilon,kappa,delta,deltal,deltau,...
  rrmax,zrmax,rrmin,zrmin,rzmax,zzmax,rzmin,zzmin,...
  Q0Q,Q1Q,Q2Q,Q3Q,Q4Q,...
  aq,rq,zq,aW,FW,...
  raqmin,qmin,raQ);

if L.nR
  LY = meqlarg(LY,raR);
end

if L.nS
  LY = meqlarg(LY,rS,zS);
end

end
