%TESTLIUTCV  TCV LIUQE tests
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function tests = testliutcv
 tests = functiontests(localfunctions);
end

%% TCV constraint FEM for initial guess
function testipmh(test)
 LX = load('tcvdb/liutcvlx.mat'); LX = LX.LX;
 L = liuqe(50000,1,'rip',NaN,'rop',NaN,'zlp',NaN,'zup',NaN,'nelem',10,'iterh',10,'bfp',[0,1]);
 tolh = [1e-4 1e-3 1e-2 1e-1 Inf];
 parfor k = 1:length(tolh)
  LY{k} = liut(L,LX,'tolh',tolh(k),'argout',{'chih'});
 end
 save(mymatfile,'tolh','L','LY')
 close all
 for k = 1:length(tolh)
  subplot(length(tolh),2,2*k-1)
  title(sprintf('tolh=%.1e #equil=%d',tolh(k),length(LY{k}.shot)))
  subplot(length(tolh),2,2*k)
  set(histogram(log10(LY{k}.chih),-1:.1:1.2),'facecolor','g')
  set(gca,'xlim',[-1 1.2]), xlabel('log_{10}\chi_h')
 end
end

%% Tune regularisation
function wregcreatetest(test)
 D = grt47;
 equil = D.Number;
 L = liu('create',equil(1),1,'iters',0);
 LX = liuxcreate(equil,1);
 wq = logspace(-9,-2,15);
 clear li LY
 for k = 1:length(wq)
  L.Wq = wq(k).*ones(L.nq,1);
  LY(k) = liut(L,LX,'wreg',wq(k));
  li(LY(k).shot,k) = LY(k).li;
 end
 li = li(equil,:); li(li == 0) = NaN;
 li3 = D.li3;
 subplot(211)
 semilogx(wq,li,'-r'), xlim([wq(1) wq(end)]), title('Current profile regularisation l_i(w_q)'), xlabel('w_q')
 subplot(212)
 plot(li3,li(:,3:2:end-2),'.','markersize',12), grid, xlabel('l_{i3}')
 legend(num2str(wq(3:2:end-2)','wq=%1.0e'),'location','best')
 orient tall
 verifyEqual(test,isempty(LY),false)
end

%% Tune regularisation
function wregtcvtest(test)
 D = load('tcvdb/liutcvsamp.mat');
 L = liuqe(50000,1,'rip',NaN,'rop',NaN,'zlp',NaN,'zup',NaN,'nelem',10,'ipmi',repmat(20,20,1)); LX = D.LX;
 wq = logspace(-9,-2,15);
 clear LY
 [li,kappa] = deal(NaN(numel(LX.t),length(wq)));
 for k = 1:length(wq)
  L.Wq = wq(k).*ones(L.nq,1);
  LY(k) = liut(L,LX,'wreg',wq(k),'argout','ie0');
  [~,l] = intersect([LX.shot' LX.t'],[LY(k).shot LY(k).t],'rows');
  li(l,k) = LY(k).li;
  kappa(l,k) = LY(k).ie0;
 end
 subplot(211)
 semilogx(wq,li,'-r'), xlim([wq(1) wq(end)]), ylim([0 2])
 title('Current profile regularisation l_i(w_q)'), xlabel('w_q')
 subplot(212)
 kappam = repmat(kappa(:,end),length(wq)); kappam(~isnan(kappa)) = NaN;
 semilogx(wq,1./kappa,'.r',wq,1./kappam,'.k'), xlim([wq(1) wq(end)])
 orient tall
 verifyEqual(test,isempty(LY),false)
end

%% Tune weights
function tunewgttest(test)
 eq = 50;
 LX = liuxcreate(eq,1);
 for Iaerr = [1 10 100 1000]
  L = liu('create',eq,1,'wIa',[ones(12,1);1e6;1e6],'Iaerr',Iaerr);
  for k = 1:6
   LY = liut(L,LX,'argout',{'Yd','Zd'},'debug',1);
   L.P.Fferr = L.P.Fferr*sqrt(mean((LY.Yd(L.kdf)-LY.Zd(L.kdf).*L.Wd(L.kdf)).^2));
   L.P.Bmerr = L.P.Bmerr*sqrt(mean((LY.Yd(L.kdm)-LY.Zd(L.kdm).*L.Wd(L.kdm)).^2));
   L.P.Fterr = L.P.Fterr*sqrt(mean((LY.Yd(L.kdt)-LY.Zd(L.kdt).*L.Wd(L.kdt)).^2));
   L = liuc(L.P,L.G);
  end
  disp([L.P.Iaerr L.P.Fferr L.P.Bmerr L.P.Fterr])
 end
 verifyEqual(test,isempty(LY),false)
end

%% Rational q surfaces
function qproftest(test)
 eq = 10;
 iqR = 1./(1:.5:5);
 [~,~,LY] = liu('create',eq,1,'itert',20,'debug',1);
 aR = liur(LY.iqQ,iqR,[0 LY.aq(end/2+1,:)]);
 subplot(1,2,2)
 plot([0 LY.aq(end/2+1,:)],1./LY.iqQ,'-+',aR,1./iqR,'o')
 verifyEqual(test,isempty(LY),false)
end

%% rho grid
function rhogridtest(test)
 eq = 1;
 np = [101  51 21 13];
 no = [256 128 64 32];
 
 pP = [(0:70)/100 linspace(.7,sqrt(2)-sqrt(1-.95.^2),51)];
 pP(72:end) = sqrt(1-(sqrt(2)-pP(72:end)).^2); pP(72) = [];
 pP = [pP .96 .97 .98 .99 .999 .9999 1];
 oP = linspace(0,2*pi,257)';
 
 pL = (0:100)/100;
 oL = (0:no(1)-1)'/no(1)*2*pi;
 [L,LX,LY] = liu('create',eq,1,'itert',150,'rip',3.5,'rop',9,'zlp',-5,'zup',5.4,'nr',256,'nz',128,'pq',pP,'noq',256);
 fsd = psitbxp2p(meq2psi(L,LY(1)),'FS',psitbxgrid('Flux','Grid',{pP,oP,NaN}));
 fsg = psitbxfsg(fsd);
 [V,D] = eig([LY.dr2FA LY.drzFA ; LY.drzFA LY.dz2FA]); o0 = atan2(-V(1,2),V(1,1));
 aQ = sqrt(2*(LY.FB-LY.FA)./(D(1,1)*cos(oP+o0).^2+D(2,2)*sin(oP+o0).^2)).*pP;
 
 subplot(221), cla
 plot(-LY.aq.*cos(oL),LY.aq.*sin(oL),'-k',-aQ.*cos(oP),aQ.*sin(oP),'-b',-fsd.x'.*cos(oP),fsd.x'.*sin(oP),'-r')
 xlim([-.2 .2]), ylim([-.2 .2]), axis square, title('FS \Delta\rho=0.01'),
 
 BAL = sign(LY.FA-LY.FB)*sqrt(LY.dr2FA  *LY.dz2FA  -LY.drzFA  ^2)/(2*pi);
 BA2L = -(LY.FA-LY.FB)*(LY.dr2FA  + LY.dz2FA );
 rAL  = LY.rA;
 BAP = sign(fsd.psimag )*sqrt(fsd.dr2mag*fsd.dz2mag-fsd.drzmag^2)/(2*pi);
 BA2P = -fsd.psimag   *(fsd.dr2mag+fsd.dz2mag);
 q0P = fsg.ngpsi0r1.x*fsd.rmag;
 q1P = 2*pi*fsd.rmag/BAP./fsg.ugpsiir1.x;
 q2P = fsg.ngpsi0r2.x*fsd.rmag^2;
 q3P = fsg.ngpsi2r2.x*fsd.rmag^2/BA2P./pP'.^2;
 q4P = fsg.ngpsi2r0.x/BA2P./pP'.^2;
 
 for k = 1:length(np)
  pL1 = linspace(.2,.7,(np(k)+1)/2);
  pL2 = linspace(.7,sqrt(2)-sqrt(1-.95.^2),(np(k)+1)/2); pL2(1) = [];
  pL2 = sqrt(1-(sqrt(2)-pL2).^2);
  pL = [0 pL1 pL2 1];
  [L,LX,LY(k)] = liu('create',eq,1,'itert',150,'rip',3.5,'rop',9,'zlp',-5,'zup',5.4,'nr',256,'nz',128,'pq',pL,'noq',no(k));
  LQ(k).pQ  = pL';
  LQ(k).q0L = LY(k).Q0Q*rAL;
  LQ(k).q1L = LY(k).Q1Q*2*pi*rAL/BAL;
  LQ(k).q2L = LY(k).Q2Q*rAL^2;
  LQ(k).q3L = LY(k).Q3Q*rAL^2/BA2L./pL'.^2;
  LQ(k).q4L = LY(k).Q4Q/BA2L./pL'.^2;
 end
 
 subplot(222), cla
 plot(LQ(1).pQ.^2,[LQ(1).q0L LQ(1).q1L LQ(1).q2L LQ(1).q3L LQ(1).q4L],'.',...
  fsg.ngpsi0r1.grid.x{1}.^2,[q0P q1P q2P q3P q4P],'-k')
 set(gca,'xlim',[0 .09],'xtick',[0 .01 .04 .09],'xticklabel',{'\rho=0' '' '0.2' '0.3'})
 legend({'Q_0' 'Q_1' 'Q_2' 'Q_3' 'Q_4' 'bicubic'},'location','best')
 title('FSA'), xlabel('\rho^2'), grid
 
 subplot(223), cla
 xt = .8:.01:1; xa = string(num2str(xt')); xa([2:10 12:end-2])=""; xa(1) = strcat("\rho=",xa(1));
 mk = 'rgb';
 plot(rhox(LQ(1).pQ),[LQ(1).q0L LQ(1).q1L LQ(1).q2L LQ(1).q3L LQ(1).q4L],'.',rhox(pP),[q0P q1P q2P q3P q4P],'-k')
 set(gca,'xlim',[0 1],'xtick',rhox(xt),'xticklabel',xa), grid
 title('FSA'), xlabel('1-(-ln(1-\rho^2))^{-1/2}'), grid on
 
 subplot(224), cla
 mk = 'rgbm';
 xa = {};
 for k = 1:length(mk)
  plot(NaN,NaN,['.' mk(k)])
  hold on
  xa{k} = sprintf('%3dx%3d',np(k),no(k));
 end
 for k = 1:length(np)
  e0 = LQ(k).q0L - resamp(q0P,pP',LQ(k).pQ);
  e1 = LQ(k).q1L - resamp(q1P,pP',LQ(k).pQ);
  e2 = LQ(k).q2L - resamp(q2P,pP',LQ(k).pQ);
  e3 = LQ(k).q3L - resamp(q3P,pP',LQ(k).pQ);
  e4 = LQ(k).q4L - resamp(q4P,pP',LQ(k).pQ);
  plot(LQ(k).pQ,[e0 e1 e2 e3 e4],['.' mk(k)])
  hold on
 end
 hold off
 legend(xa,'location','best')
 set(gca,'xlim',[.2 .95]), xlabel('\rho')
 title('FSA error')
 orient tall
 
 verifyEqual(test,isempty(LY),false)
end

%% rz grid
function asxytest(test)
 eq = 1;
 nz = [32  64 128  32  64 128];
 nr = [31  61 121  46  90 180];
 n = length(nz);
 ka = [3 4  9 10];
 Bz = 1e-2;
 L = liu('create',eq,1);
 r1 = (L.G.rx(end)+L.G.rx(1))/2;
 z1 = (L.G.zx(end)+L.G.zx(1))/2;
 Bra = greenem('br',[r1 r1+.1],[z1 z1],L.G.rw,L.G.zw) * L.G.Twa;
 Bza = greenem('bz',[r1 r1+.1],[z1 z1],L.G.rw,L.G.zw) * L.G.Twa;
 Ia = zeros(L.G.na,1);
 Ia(ka) = [Bra(:,ka);Bza(:,ka)]\[0;0;0;Bz];
 Hrr = greenem('dr1r1m',r1,z1,L.G.rw,L.G.zw) * L.G.Twa * Ia;
 Hzz = greenem('dz1z1m',r1,z1,L.G.rw,L.G.zw) * L.G.Twa * Ia;
 Hrz = greenem('dr1z1m',r1,z1,L.G.rw,L.G.zw) * L.G.Twa * Ia;
 H0 = sqrt(Hrz.*Hrz-Hrr.*Hzz)';
 
 for l = 1:n
  L = liu('create',eq,1,'nz',nz(l),'nr',nr(l));
  Mba = L.G.Mxa(~L.lxy,:);
  txt{l} = sprintf('%dx%d',L.nzx,L.nrx);
  Fb = Mba*Ia;
  Fx0 = gszrmex(Fb,zeros(L.nzx-2,L.nrx-2),L.cx,L.cq,L.cr,L.cs,L.ci,L.co,0);
  [zA,rA,FA,dz2FA,dr2FA,drzFA,zX,rX,FX,dz2FX,dr2FX,drzFX] = asxymex(Fx0,L.G.zx,L.G.rx,L.P.dasm,L.dzx,L.drx,L.idzx,L.idrx,true(size(Fx0)),10);
  eX(l,1:2) = sqrt(((r1-rX).^2+(z1-zX)^2)./[1 L.drx^2+L.dzx^2]);
  eX(l,3) = abs(sqrt(drzFX*drzFX-dr2FX*dz2FX)./H0-1);
  for k = 1:L.G.na
   Fb = Mba(:,k);
   Fx = gszrmex(Fb,zeros(L.nzx-2,L.nrx-2),L.cx,L.cq,L.cr,L.cs,L.ci,L.co,0);
   eFx = Fx./reshape(L.G.Mxa(:,k),size(Fx))-1;
   eF(k,l) = max(abs(eFx(:)));
  end
 end
 subplot(211)
 bar(eF(1:end-1,:)), set(gca,'yscale','log'), legend(txt,'location','best')
 title('Flux error')
 subplot(223)
 contour(L.G.rx,L.G.zx,Fx0,50,'-r')
 hold on, plot(r1,z1,'xk'), hold off, axis equal
 title('X point')
 subplot(224)
 bar(eX')
 set(gca,'yscale','log','xticklabel',{'[m]' '[\Delta]' 'H'})
 title('Error')
 orient tall
 
 verifyEqual(test,0,0)
end

%% rz grid
function rtcitest(test)
 L = liu('create',1,1);
 rB = 2; zB = 4;
 nz = [32  64 128];
 n = length(nz);
 aB = repmat(rB/2,L.noq,n);
 for l = 1:n
  L = liu('create',1,1,'nz',nz(l));
  txt{l} = sprintf('%dx%d',L.nzx,L.nrx);
  rA = (L.rx(end)+L.rx(1))*0.5; zA = 0;
  Fx = ((L.rx'-rA)/rB).^2 + ((L.zx-zA)/zB).^2;
  for ki = 1:200
   aB(:,l) = rtcimex(aB(:,l),rA,(rA-L.G.rx(1))*L.idrx,(zA-L.G.zx(1))*L.idzx,Fx,1,L.idrcosq,L.idzsinq,L.cosq,false,L.drx);
  end
 end
 aB0 = 1./sqrt(((cos(L.oq)/rB).^2+(sin(L.oq)/zB).^2));
 subplot(221)
 plot(rA,zA,'ok')
 hold on
 contour(L.rx,L.zx,Fx,linspace(0,3,31),'-r')
 contour(L.rx,L.zx,Fx,[1 1],'-k')
 hold off, axis equal
 subplot(222)
 semilogy(L.oq,aB0-aB), xlabel('\theta'), xlim([0 2*pi]), title('LCFS position error')
 legend(txt,'location','best')
 orient tall
 
 verifyEqual(test,0,0)
end

function y=rhox(rho)
 y = 1-(-log(1-rho.^2)).^-.5;
end

function par = createp(L)
 wBm = L.P.wBm;
 wFf = L.P.wFf;
 wBm(~contains(extractBetween(L.G.dimm,4,5),{'AA' 'AB' 'AL' 'AD' 'AE'})) = 0;
 wFf(~contains(extractBetween(L.G.dimf,4,5),{'AA' 'AB' 'AL' 'AD' 'AE'})) = 0;
 par = {'nelem',10,'iterh',repmat(20,20,1),'tolh',1e-3,'zlp',-3.7,...
        'iters',0,'psichco',1e-6,...
        'wBm',wBm,'wFf',wFf};
end

function [equil,t] = grt47
 [~,p] = system('find /Lac8_D/CREATE0883/equil -regextype grep -regex  ".*Equil_[0-9]\{5\}_C[N]\?L4E.mat"');
 p = regexp(split(p(1:end-1),newline),'.*Equil_([0-9]{5})_C([N]?)L4E.mat','tokens');
 n = length(p);
 [equil,t] = deal(zeros(n,1));
 for k = 1:n
  equil(k) = str2double(p{k}{1}{1});
  t    (k) =  strcmp(p{k}{1}{2},'N')+1;
 end
 [~,k] = sortrows([equil t]);
 equil = equil(k);
 t     = t    (k);
end

function f = mymatfile
 f = dbstack; f = [f(2).file(1:end-2) '/' f(2).name];
end


