function [LY] = meqpost(L,LY,TpDg,ITpDg,ag,...
  Fx,FA,FB,rA,zA,dr2FA,dz2FA,drzFA,rB,zB,lB,lX,rX,zX,FX,...
  rBt,Ia,Iu,Iy,Opy)
%MEQPOST Post-processing for various meq functions
% LY = meqpost(L,LY,TpDg,ITpDg,ag,...
%  Fx,FA,FB,rA,zA,dr2FA,dz2FA,drzFA,rB,zB,lX,...
%  rBt,Ia,Iu,Iy,Opy)
% Computes integral quantities, estimated magnetic measurements,
% profiles, vessel currents. Given:
% MEQ structure L, 
% LY with shot, t, and containing optional initial aq,aW guesses.
% TpDg,ITpDg from bf(1,. call (see BFHELP),
% flux map Fx(Z,R), axis, boundary flux FA,FB, axis position and
% flux gradients rA,zA,dr2FA,dz2FA,drzFA.
% Boundary flux position rB,zB. Locicals lB, lX indicating whether boundary
% was found and whether the plasma is diverted.
% rX,zX,FX positions and fluxes of all X points.
% rBt the vacuum r*Bt.
% Ia,Iu,Iy the circuit, vessel and plasma currents.
% Opy the plasma current mask.
%
% if L.P.iterq>0, then initial guesses for aq and aW can be supplied in LY
%
%   See MEQINT, MEQPOSTQ, MEQPROF, BFHELP
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% bfct
fPg = L.fPg; fTg = L.fTg;

nA = int32(numel(FA)); % number of axes
nB = int32(numel(FB)); % number of bounded domains
nX = int32(numel(FX)); % number of X points

% When L.P.LYall is true, some slices can have more domains than is allowed
if nA>L.nD
  nA    = L.nD;
  FA    =    FA(1:L.nD);
  rA    =    rA(1:L.nD);
  zA    =    zA(1:L.nD);
  dr2FA = dr2FA(1:L.nD);
  dz2FA = dz2FA(1:L.nD);
  drzFA = drzFA(1:L.nD);
end
if nB>L.nD
  nB = L.nD;
  FB = FB(1:L.nD);
  rB = rB(1:L.nD);
  zB = zB(1:L.nD);
end


% Total integrals
[Ip,Wk,Wp,Wt,Wt0,WN,Vp,Ft,Ft0,bp,bt,mu,li,bpli2] = ...
  meqint(fPg,fTg,TpDg,ITpDg,ag,Fx,Opy,L,rBt);

rIp = sum(L.rry(:).*Iy(:));
zIp = sum(L.zzy(:).*Iy(:));

[F0,F1] = meqpdomlim(FA,FB,L.nD);
FR = calcFR(F0,F1,FX,lX); % X-point/boundary normalized flux ratio

% Integrals per domain
[IpD,WkD,WpD,WtD,Wt0D,WND,VpD,FtD,Ft0D,bpD,btD,liD,zIpD,rIpD] = deal(zeros(L.nD,1));
for iD=1:nB
  [IpD(iD),WkD(iD),WpD(iD),WtD(iD),Wt0D(iD),WND(iD),VpD(iD),...
    FtD(iD),Ft0D(iD),bpD(iD),btD(iD),~,liD(iD),~] = ...
    meqint(fPg,fTg,TpDg(iD,:),ITpDg(iD,:),ag,Fx,Opy==iD,L,rBt);
  
  % Centroids per domain
  rIpD(iD) = sum(L.rry(Opy(:)==iD).*Iy(Opy(:)==iD));
  zIpD(iD) = sum(L.zzy(Opy(:)==iD).*Iy(Opy(:)==iD));
end

% cumulative integrals per domain
[IpQ,~,~,VpQ,FtQ,Ft0Q,~] = ...
  meqintQ(L,F0,F1,rBt,ag,Fx,Opy,sign(Ip));
FtPQ = FtQ+Ft0Q;

% q axis using small diamagnetic approximation (LIUQE 2015 paper eq. 88)
qA = zeros(nA,1);
[gAg,IgAg] = L.bfct(5,L.bfp,[],F0,F1); IgAg = 2e-7*L.idsx*IgAg;
for iA=1:nA
  TyAg = gAg(iA,:).*(rA(iA)*fPg+fTg/rA(iA)).'; % TyAg: Axis current (jphi.dR.dZ) contribution from each basis function
  cqA = 4e-7*pi*rA(iA)^2*rBt*sqrt(dr2FA(iA)*dz2FA(iA)-drzFA(iA)^2)/(abs(dr2FA(iA)+dz2FA(iA))*L.dsx);
  qA(iA) = (rBt.^2+(IgAg(iA,:).*fTg.')*ag)/(cqA*TyAg*ag);
end

% Profiles
[PpQ,TTpQ,PQ,TQ,iTQ,PpQg,TTpQg] = ...
  meqprof(fPg,fTg,ag,L.pQ(:).^2,F0,F1,rBt,L.bfct,L.bfp,L.idsx,L.smalldia);
% Temporary fix, store only one domain for PpQg/TTpQg
PpQg = squeeze(PpQg(:,1,:));
TTpQg = squeeze(TTpQg(:,1,:));

% Fictional measurements
Bm  = L.G.Bma*Ia + L.G.Bmx(:,L.lxy)*Iy(:) + L.G.Bmu*Iu;
Ff  = L.G.Mfa*Ia + L.G.Mfx(:,L.lxy)*Iy(:) + L.G.Mfu*Iu;

% Vessel current from generic representation
Iv = L.G.Tvu*Iu;
Is = L.G.Tius*Iu;

% Store in LY
LY = meqlarg(LY,...
  Fx,rB,zB,FB,rA,zA,FA,lB,lX,dr2FA,dz2FA,drzFA,rX,zX,FX,nX,Opy,nA,nB,FR,...   % Flux map
  Iy,Ia,Iu,Iv,Is,Bm,Ff,ag,rIp,zIp,...                                      % Currents, measurements
  IpQ,VpQ,PpQg,TTpQg,PpQ,TTpQ,PQ,TQ,iTQ,FtPQ,...                        % integrals per Q
  IpD,zIpD,rIpD,WkD,WpD,WtD,Wt0D,WND,VpD,FtD,Ft0D,bpD,btD,liD,...       % integrals per domain
  Ip,zIp,rIp,rBt,Wk,Wp,Wt,Wt0,WN,Vp,Ft,Ft0,bp,bt,mu,li,bpli2,qA...      % integrals
  );

if L.P.iterq
  LY = meqpostq(L,LY);
end

% optional interpolation of flux, fields in desired points
if L.nn
  [Fn,Brn,Bzn,Brrn,Brzn,Bzrn,Bzzn] = L.P.infct(L.rx,L.zx,Fx,L.P.rn,L.P.zn,L.inM);
  LY = meqlarg(LY,Fn,Brn,Bzn,Brrn,Brzn,Bzrn,Bzzn);
end

% optional computation of Br,Bz fields on x grid
if L.P.ifield
  [Brx,Bzx] = meqBrBz(Fx,L.i2pirxdzx,L.i2pirxdrx,L.nzx,L.nrx);
  Btx = meqBt(L,Fx,Opy,ag,rBt,F0,F1,TQ);
  LY = meqlarg(LY,Brx,Bzx,Btx);
end

% optional computation of fields/fluxes on z grid
if L.P.izgrid
  F0z = reshape(L.G.Mza*LY.Ia + L.G.Mzu*LY.Iu,L.nzz,L.nrz);
  Fz  = F0z + reshape(L.G.Mzy*Iy(:),L.nzz,L.nrz);
  LY = meqlarg(LY,Fz);
  if L.P.ifield
    [Brz,Bzz] = meqBrBz(Fz,L.i2pirzdzz,L.i2pirzdrz,L.nzz,L.nrz);
    LY = meqlarg(LY,Brz,Bzz);
  end
end

% optional vacuum field/flux computation
if L.P.ivacuum
  if all(Iy(:)==0)
    % no plasma, vacuum values same as Br,Bz,Fx
    LY.F0x = Fx;
    if L.P.izgrid,                LY.F0z  = Fz;                 end
    if L.P.ifield,                LY.Br0x = Brx; LY.Bz0x = Bzx; end
    if L.P.ifield && L.P.izgrid,  LY.Br0z = Brz; LY.Bz0z = Bzz; end
  else
    LY.F0x = meqFx(L,zeros(L.nzy,L.nry),[Ia;Iu]); % recompute flux with Iy=0
    if L.P.izgrid,                LY.F0z = F0z;                 end
    if L.P.ifield,               [LY.Br0x,LY.Bz0x] = meqBrBz(LY.F0x,L.i2pirxdzx,L.i2pirxdrx,L.nzx,L.nrx); end
    if L.P.ifield && L.P.izgrid, [LY.Br0z,LY.Bz0z] = meqBrBz(LY.F0z,L.i2pirzdzz,L.i2pirzdrz,L.nzz,L.nrz); end
  end
end

end

function FR = calcFR(F0,F1,FX,lX)
%% Calculate ratio of normalized flux at x-point vs normalized flux at boudnary.
% For limited cases, FR>1 indicating the distance from the diverted case
% FR=1;
% For doublet cases, FR<1 for mantle domains, indicating the ratio of flux
% contained within the mantle w.r.t. maximum flux contained in the
% separatrix-enclosed domains.

sFBA = sign(F1(1)-F0(1)); % +1 is flux is increasing
nB = size(F1,1);
FR = zeros(nB,1);
for iB = 1:nB
  if F0(iB)==F1(iB),  FR(iB)=nan; % empty domain -> nan
  elseif  lX(iB),     FR(iB)=1;   % diverted -> 1
  elseif isempty(FX), FR(iB)=inf; % no X points-> inf
  else
    FB = F1(iB);
    if any(F0(iB)==FX)
      % mantle domain - choose minimum axis flux and inner limiting X point.
      FA = sFBA*min(F0*sFBA);
      FXi = F0(iB);
      FR(iB) = (FXi-FA)/(FB-FA);
    else % non-diverted axis domain
      FA = F0(iB); % axis flux for this domain
      [~,iX] = min(abs(FX-FB)); % closest X point in terms of flux
      FXi = FX(iX);
    end
    FR(iB) = (FXi-FA)/(FB-FA);
  end
end
end

function [Br,Bz] = meqBrBz(Fx,i2pirdz,i2pirdr,nz,nr)
% [Br,Bz] = meqBrBz(Fx,i2pirdz,i2pirdr,nz,nr)
% Compute Br,Bz fields
% General version that also accepts time-varying Fx

[Br,Bz] = deal(zeros(nz,nr,size(Fx,3))); % init

% Br = -1/(2*pi*R)* dF/dz
% Central differences dF/dz[i] =  F[i-1] - F[i+1]/(2*dz)
Br(2:end-1,:,:) = -i2pirdz.* (Fx(3:end,:,:) - Fx(1:end-2,:,:));
% At grid boundary i, use: dF/dz[i] = (-F(i-2) + 4*F(i-1) - 3*F(i))/(2*dz)
Br(end,:  ,:) = -i2pirdz          .* (+Fx(end-2,:,:) - 4*Fx(end-1,:,:) + 3*Fx(end,:,:));
Br(1  ,:  ,:) = -i2pirdz          .* (-Fx(    3,:,:) + 4*Fx(    2,:,:) - 3*Fx(  1,:,:));

% Bz = 1/(2*pi*R)* dF/dr
Bz(:,2:end-1,:) =  i2pirdr(2:end-1) .* (+Fx(:,  3:end,:) - Fx(:,1:end-2,:));
% Same as for Br
Bz(:,end    ,:) =  i2pirdr(end)     .* (+Fx(:,end-2,:) - 4*Fx(:,end-1,:) + 3*Fx(:,end,:));
Bz(:,1      ,:) =  i2pirdr(1)       .* (-Fx(:,    3,:) + 4*Fx(:,    2,:) - 3*Fx(:,  1,:));
end

function Btx = meqBt(L,Fx,Opy,ag,rBt,F0,F1,TQ)
% Btx = meqBt(L,Fx,Opy,ag,rBt,F0,F1,TQ)
% Computes toroidal field on x grid

Btx  = rBt*repmat(1./L.rx',L.nzx,1);
Bty = Btx(L.lxy);

nB = numel(F1);
for iB = 1:nB
  Opyi = (Opy==iB); % mask for this domain
  if isequal(func2str(L.bfct),'bf3pmex')
    Btyi = L.bfct(8,L.bfp, Fx,F0(iB),F1(iB),int8(Opyi),ag(:,iB),rBt,L.idsx,L.iry);
    Bty(Opyi) = Btyi(Opyi);
  else
    % Mode 8 not yet available for other bfs, to be added later
    FyN = (Fx(L.lxy)-F0(iB))/(F1(iB)-F0(iB)); % rhopol.^2
    Bty(Opyi) = interp1(L.pQ.^2,TQ(:,iB)',FyN(Opyi))./L.rry(Opyi);
  end
end
Btx(L.lxy) = Bty;
end
