%% Demonstration of boundary current regularization
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

clear all;
clc;
clf;

[L,LX,LY0] = liuqe(61400,1,'bfct',@bfef,'bfp',[2 4],'wreg',5e-8*[1 1 1]);
assert(~isempty(LY0),'0 failed');

%%
[L,LX,LY1] = liuqe(61400,1,'bfct',@bfef,'bfp',[2 4],'wreg',5e-8*[1,1,1e-5]);
assert(~isempty(LY1),'1 failed');

%%
[L,LX,LY2] = liuqe(61400,1,'bfct',@bfef,'bfp',[2 4],'wreg',5e-8*[1,1,1e+4]);
assert(~isempty(LY2),'2 failed');
%%
LY0.jB = -2*pi*(LY0.PpQ*LY0.rA+LY0.TTpQ/(4e-7*pi)/LY0.rA); % jB ~ -2*pi*(rA*p'(1) + TT'(1)/mu0/rA)
LY1.jB = -2*pi*(LY1.PpQ*LY1.rA+LY1.TTpQ/(4e-7*pi)/LY1.rA); % 
LY2.jB = -2*pi*(LY2.PpQ*LY2.rA+LY2.TTpQ/(4e-7*pi)/LY2.rA); % 

%%
subplot(131)
contour(L.rx,L.zx,LY0.Fx-LY0.FB,21,'b'); hold on;
contour(L.rx,L.zx,LY1.Fx-LY1.FB,21,'r')
contour(L.rx,L.zx,LY2.Fx-LY2.FB,21,'k')

legstr{1} = ['j_B = ',sprintf('%2.2e',LY0.jB(end))];
legstr{2} = ['j_B = ',sprintf('%2.2e',LY1.jB(end))];
legstr{3} = ['j_B = ',sprintf('%2.2e',LY2.jB(end))];

legend(legstr)


subplot(232);
plot(L.pQ.^2,LY0.PpQ,'b',...
  L.pQ.^2,LY1.PpQ,'r',...
  L.pQ.^2,LY2.PpQ,'k');
title('p''')
subplot(233);
plot(L.pQ.^2,LY0.TTpQ,'b',...
  L.pQ.^2,LY1.TTpQ,'r',...
   L.pQ.^2,LY2.TTpQ,'k');
title('TT''')

subplot(235);
plot(L.pQ.^2,LY0.PQ,'b',...
  L.pQ.^2,LY1.PQ,'r',...
  L.pQ.^2,LY2.PQ,'k');
title('p')
xlabel('\psiN')

subplot(236);
plot(L.pQ.^2,LY0.TQ,'b',...
  L.pQ.^2,LY1.TQ,'r',...
   L.pQ.^2,LY2.TQ,'k');
title('T')
xlabel('\psiN')

shg
