%TESTCREATE CREATE(ITER) LIUQE tests
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function tests = testcreate
 tests = functiontests(localfunctions);
end

%% CREATE DATABASE2
function testdatabase2(test)
 [~,equil,cnl] = createf(10000:999999,0);
 L = grt047(equil(1),cnl(1),'wreg',1e-5,'nz',128,'Iaerr',100);
 LX = liuxcreate(equil,cnl,L);
 nt = numel(LX.t);
 for k = nt:-1:1
  L.Wq(:) = 0;
  ly = liut(L,meqxk(LX,k));
  if isempty(ly) || ly.Wk<0
   L.Wq(:) = 1e-5;
   ly = liut(L,meqxk(LX,k));
  end
  LY(k) = ly;
 end
 LY = meqlpack(LY,{'ppQg' 'TTpQg'});
 save(mymatfile)
end

%% CREATE STREAMING SAMPLE FILE
% function teststreamingsample(test)
%  [filename,equil,cnl] = createf('../data/Case001_bellshape_curdriv_t625d000_t629d140_CNL4LIUQE.mat');
%  L = grt047(filename,cnl,'wreg',1e-5,'nz',128,'itert',1,'debug', 0,'bfct',@bf3pmex);
%  LX = liuxcreate(filename,cnl,L);
%  LY = liut(L,LX);
%  save(mymatfile)
% end

%% CREATE GRT047/S0 grid scan
function testgrt047s01(test)
 [~,equil,cnl] = createf(1:999,0);
 p = { {'nr',256,'nz',128} {'nz',128} {'ri',3.5,'ro',9,'zl',-5,'zu',5.4,'nr',128,'nz',64} {} };
 parfor k=1:length(p)
  [L{k},LX{k},LY{k},D{k}] = grt047(equil,cnl,p{k}{:});
 end
 save(mymatfile)
 close all
 for k=1:length(p)
  grt047plot(L{k},LX{k},LY{k},D{k},sprintf(' %dx%d GRT047/S0',L{k}.nrx,L{k}.nzx))
 end
 myprint
end

%% CREATE GRT047/S1 Iaerr scan
function testgrt047s02(test)
 [~,equil,cnl] = createf(1000:9999,0);
 p = { {'Iaerr',1} {'Iaerr',10} {'Iaerr',100} {'Iaerr',1000} };
 parfor k=1:length(p)
  [L{k},LX{k},LY{k},D{k}] = grt047(equil,cnl,'nz',128,p{k}{:});
 end
 save(mymatfile)
 
 close all
 for k=1:length(p)
  grt047plot(L{k},LX{k},LY{k},D{k},sprintf(' %gA GRT047/S2',L{k}.P.Iaerr))
 end
 myprint
end

%% CREATE GRT047 DZ
function testdzinstab(test)
 p = { {'bfct',@bfabmex} {} };
 close all
 for k = 1:length(p)
  figure
  [L{k},LX,LY{k}] = grt047(436,0,'debug',1,'nz',128,'psichco',2.5e-4,p{k}{:});
 end
 nN = 251;
 g1 = interp1(LX.FNQ,LX.PpQ * L{1}.dsx     ,linspace(0,1,nN)');
 g2 = interp1(LX.FNQ,LX.TTpQ*(L{1}.dsx/mu0),linspace(0,1,nN)');
 g3 = zeros(size(g1));
 [Ig1,Ig2,Ig3] = bfprmex(g1,g2,g3);
 bfp = struct('gNg',[g1 g2 g3],'IgNg',[Ig1 Ig2 Ig3],'fPg',[1;0;0],'fTg',[0;1;0]);
 L{3} = liuc(L{1}.P,L{1}.G,'bfct',@bf3imex,'bfp',bfp);
 figure
 warning off
 LY{3} = liut(L{3},LX,'debug',1);
 
 e = { 'LY{k}.Fx' 'Fx0' 'Fx1' 'Fx0+Fx1' 'Fxe' 'LY{k}.Fx' };
 t = { '\psi'     '\psi_{p0}' '\psi_{p1}' '\psi_p' '\psi_e' '\psi_{tot}' };
 er = { 'LY{k}.rA' 'r0' '[]' 'rp' '[]' 'rt' };
 ez = { 'LY{k}.zA' 'z0' '[]' 'zp' '[]' 'zt' };
 for k = 1:length(LY)
  Iy0 = (LY{k}.Iy + interp1(L{k}.zy,LY{k}.Iy,2*LY{k}.zA-L{k}.zy,'linear',0)) * 0.5;
  Iy1 = LY{k}.Iy - Iy0;
  Fx0 = gszrmex(meqfbp(Iy0,L{k}),Iy0,L{k}.cx,L{k}.cq,L{k}.cr,L{k}.cs,L{k}.ci,L{k}.co,0);
  Fx1 = gszrmex(meqfbp(Iy1,L{k}),Iy1,L{k}.cx,L{k}.cq,L{k}.cr,L{k}.cs,L{k}.ci,L{k}.co,0);
  Fxe = gszrmex(L{k}.Mbe*[LY{k}.Ia;LY{k}.Is],zeros(size(Iy0)),L{k}.cx,L{k}.cq,L{k}.cr,L{k}.cs,L{k}.ci,L{k}.co,0);
  [r0,z0] = meqpdom(Fx0,LY{k}.Ip,L{k}.P.isaddl,L{k});
  [rp,zp] = meqpdom(Fx0+Fx1,LY{k}.Ip,L{k}.P.isaddl,L{k});
  [rt,zt] = meqpdom(Fx0+Fx1+Fxe,LY{k}.Ip,L{k}.P.isaddl,L{k});
  
  figure
  for l = 1:length(e)
   x = eval(e{l}); r = eval(er{l}); z = eval(ez{l});
   subplot(4,3,l)
   contour(L{k}.rx,L{k}.zx,x,20), axis equal
   title(t{l})
   subplot(4,3,l+6)
   contour(L{k}.rx,L{k}.zx,x,200,'-k'), axis equal
   hold on, plot(r,z,'om',[0 100],LY{k}.zA([1 1]),'-m'), hold off
   axis equal, xlim([LY{k}.rA-.25,LY{k}.rA+.75]), ylim([LY{k}.zA-.5,LY{k}.zA+.5])
   title(t{l})
  end
 end
 
 save(mymatfile)
 myprint
end

%% CREATE GRT047/S0 DML problem
function testgrt047dml(test)
 [~,equil,cnl] = createf(1:999,0);
 [L,LX,LY] = grt047(equil,cnl,'nz',128,'iterq',0);
 save(mymatfile)
 
 close all
 [~,kl,kx] = intersect(LY.shot,LX.shot,'stable');
 subplot(3,1,1)
 histogram(LY.Ft(kl)-LX.Ft(kx)')
 title('Diamagnetic flux mismatch [Wb]')
 subplot(3,1,2)
 histogram(100-100*LY.Ft(kl)./LX.Ft(kx)')
 title('Diamagnetic flux mismatch [%]')
 subplot(3,1,3)
 histogram(LY.Ft(kl)./LY.Ft0(kl))
 title('Diamagnetic/vacuum flux')
 myprint
end

function [L,LX,LY,D] = grt047(equil,cnl,varargin)
 [L,LX,LY] = liu('create',equil,cnl,'iterq',100,'psichco',1e-6,'itera',200,...
  'Iaerr',.01,'idml',0,'Iperr',10,'selfm','55\.A[ABDEL].*',...
  'bfct',@bfef,'bfp',[3 3],'wreg',0,...
  'zfct',@liudz,'stab',1,'stabint',.1,...
  varargin{:});
 if isfield(LY,'aW')
  LY.aW(2,:) = L.G.aW(2) - LY.aW(2,:);
  LY.aW(1:2,~LY.lX) = NaN;
 end
 if nargout > 3
  D = creategaps(L.G.dimW);
 end
end

function grt047plot(L,LX,LY,D,t)
 figure %Gaps
 [~,kl,kd] = intersect(LY.shot,D.shot,'stable');
 subplot(4,2,1)
 boxplot(LY.aW(:,kl)'-D.aW(:,kd)',L.G.dimW,'datalim',[-.02 .02],'extrememode','clip','labelorientation','inline')
 title(['Gap errors [m] ' t])
 subplot(4,2,2)
 semilogx(LY.chi(kl),LY.dz(kl),'.')
 title('dz[m] vs \chi')
 kW = find(startsWith(L.G.dimW,'c'));
 for k = 1:length(kW)
  subplot(4,2,k+2)
  semilogx(LY.chi(kl),LY.aW(kW(k),kl)-D.aW(kW(k),kd),'.'), ylim([-.02 .02])
  ylim([-.02 .02]), title([L.G.dimW{kW(k)} ' error [m]'])
 end
end

%% CREATE bf3i
function testbf3i(test)
 [~,equil,cnl] = createf;
 k = equil > 0 & equil < 1000 & cnl == 0;
 equil = equil(k); cnl = cnl(k);
 LX = liuxcreate(equil,cnl);
 nN = 1001;
 bfp = struct('gNg',zeros(nN,3),'IgNg',zeros(nN,3),'fPg',[1;0;0],'fTg',[0;1;0]);
 L = liu('create',equil(1),cnl(1),'selfm','55\.A[ABDEL].*','iterq',200,'itera',500,'nz',128,...
  'bfct',@bf3imex,'bfp',bfp,...
  'Iaerr',.01,'idml',0,'Iperr',10,...
  'zfct',@liudz,'stab',1,'stabint',.1);
 L = liuc(L.P,L.G);
 parfor k = 1:nt
  g1 = interp1(LX.FNQ(:,k),LX.PpQ (:,k)* L.dsx     ,linspace(0,1,nN)');
  g2 = interp1(LX.FNQ(:,k),LX.TTpQ(:,k)*(L.dsx/mu0),linspace(0,1,nN)');
  g3 = zeros(size(g1));
  [Ig1,Ig2,Ig3] = bfprmex(g1,g2,g3);
  bfp = struct('gNg',[g1 g2 g3],'IgNg',[Ig1 Ig2 Ig3],'fPg',[1;0;0],'fTg',[0;1;0]);
  warning off
  ly{k} = liut(L,meqxk(LX,k),'bfp',bfp);
  warning on
 end
 l = 0;
 for k = 1:nt
  if ~isempty(ly{k})
   l = l+1;
   LY(l)=ly{k};
  end
 end
 LY = meqlpack(LY,{'ppQg' 'TTpQg'});
 g = [1:2 7:25];
 gtxt = [num2str((1:2)','cg%d') ; num2str((1:19)','g%d')];
 gtxt([17 14 12 5],:) = num2str((3:6)','cg%d');
 LY.aW(2,:) = L.G.aW(2) - LY.aW(2,:);
 LY.aW = LY.aW(g,:);

 D = creategaps;
 D.aW = D.aW(g,:);
 
 clf
 [~,kl,kd] = intersect(LY.shot,D.shot,'stable');
 [~,kl,kx] = intersect(LY.shot,LX.shot,'stable');
 subplot(3,1,1)
 boxplot(LY.aW(:,kl)'-D.aW(:,kd)',gtxt,'datalim',[-.02 .02],'extrememode','clip')
 title('Gap errors [m]')
 subplot(3,1,2)
 boxplot(LY.ag(1:2,:)','datalim',[.98 1.02],'extrememode','clip')
 subplot(3,2,5)
 histogram(LY.dz(kl),'binlimits',[-.001 .001])
 title('Vertical shift [m]')
 subplot(3,2,6)
 histogram(LY.Ft(kl)-LX.Ft(kx)')
 title('Diamagnetic flux mismatch [Wb]')

 
 save /tmp/tmp.mat
end

%% CREATE constraint triangular FE for initial guess
function testipmh(test)
 [~,equil,cnl] = createf(1:10000,0);
 L = liu('create',equil(1),cnl(1),'iterq',0,'selfm','55\.A[ABDEL].*');
 LX = liuxcreate(equil,cnl,L);
 
 tolh = [1e-3 1e-2 1e-1 Inf];
 parfor k = 1:length(tolh)
  LY{k} = liut(L,LX,'tolh',tolh(k),'argout',{'hfinfo' 'chih'});
 end
 save(mymatfile,'tolh','L','LY')
 
 close all
 for k = 1:length(tolh)
  subplot(length(tolh),2,2*k-1)
  semilogy(LY{k}.hfinfo), xlim([1 size(LY{k}.hfinfo,1)]), ylim([1e-3 Inf]), xlim([1 10])
  title(sprintf('tolh=%.1e #=%d/%d',tolh(k),length(LY{k}.shot),length(LX.shot)))
  subplot(length(tolh),2,2*k)
  set(histogram(log10(LY{k}.chih(LY{k}.shot<1000)),-1:.1:1.2),'facecolor','g'), hold on
  set(histogram(log10(LY{k}.chih(LY{k}.shot>1000)),-1:.1:1.2),'facecolor','r'), hold off
  set(gca,'xlim',[-1 1.2]), xlabel('log_{10}\chi_h')
  legend({'S1', 'S2'})
 end
 myprint
 
end

%% Gap comparison for GRT47 equilibria
function testgrt47(test)
 
 %% Setup
 [~,equil,cnl] = createf;
 k = equil > 0 & cnl == 0;
 equil = equil(k); cnl = cnl(k);
 LX = liuxcreate(equil,cnl);
 parEFIT = {'ri',3.5,'ro',9,'zl',-5,'zu',5.4,'nr',128};
 parz = {'zbfc',@liuz0,'stabz',0,'stab',.8,'stabcint',.15,'itera',1000};
 LL{1} = liu('create',equil(1),cnl(1),'nz',128,parz{:},'selfm','55\.A[ABDEL].*');
 LL{2} = liu('create',equil(1),cnl(1),'nz', 64,parz{:},'selfm','55\.A[ABDEL].*',parEFIT{:});
 LL{3} = liu('create',equil(1),cnl(1),'nz', 32,parz{:},'selfm','55\.A[ABDEL].*');
 
 %% Calculation
 parfor k = 1:length(LL)
  L = LL{k};
  LY{k} = liut(LL{k},LX);
 end
 L = LL;
 
 %  save /Lac8_D/CREATE0883/SPCoutput/efit47.mat L LX LY
 save(mymatfile,'L','LX','LY')
 
 %% Plot
 g = [1:2 7:25];
 gtxt = [num2str((1:2)','cg%d') ; num2str((1:19)','g%d')];
 gtxt([17 14 12 5],:) = num2str((3:6)','cg%d');
 for kr = 1:length(LY)
  LY{kr}.aW(2,:) = L{kr}.G.aW(2) - LY{kr}.aW(2,:);
  LY{kr}.aW = LY{kr}.aW(g,:);
 end
 close all, kfig = 0;
 
 kfig = kfig+1; figure(kfig)
 l = LY{1}.shot < 1000;
 subplot(2,2,1)
 semilogy(LY{1}.dz(l),LY{1}.chi(l),'.')
 xlim([min(LY{1}.dz ) max(LY{1}.dz )]), xlabel('dz [m]')
 ylim([min(LY{1}.chi) max(LY{1}.chi)]), ylabel('\chi')
 grid on, title('S01')
 subplot(2,2,2)
 semilogy(LY{1}.dz(~l),LY{1}.chi(~l),'.')
 xlim([min(LY{1}.dz ) max(LY{1}.dz )]), xlabel('dz [m]')
 ylim([min(LY{1}.chi) max(LY{1}.chi)]), ylabel('\chi')
 grid on, title('S02')
 
 kfig = kfig+1; figure(kfig)
 equil = 435;
 lx = liuxcreate(equil,0);
 g1 = lx.PpQ*mu0;
 g2 = lx.TTpQ;
 g3 = zeros(size(g1));
 [Ig1,Ig2,Ig3] = bfprmex(g1,g2,g3);
 bfp = struct('gNg',[g1 g2 g3],'IgNg',[Ig1 Ig2 Ig3],'fPg',[1;0;0],'fTg',[0;1;0]);
 [~,~,ly] = liu('create',equil,0,'bfct',@bf3imex,'bfp',bfp);
 subplot(2,3,1)
 jpp  = lx.PpQ' *lx.rA    ;
 jttp = lx.TTpQ'/lx.rA/mu0;
 fp   = linspace(0,1,length(jpp));
 plot(sqrt(fp),[jpp;jttp;jpp+jttp]/1e6)
 ylim([-4 6]),xlabel('\rho'), ylabel('j_\phi [MA/m2]')
 title(sprintf('CREATE #%d',equil))
 subplot(2,3,2)
 l = find(LY{1}.shot == equil);
 jpp  = LY{1}.PpQ(:,l)' *LY{1}.rA(l)    *2*pi;
 jttp = LY{1}.TTpQ(:,l)'/LY{1}.rA(l)/mu0*2*pi;
 fp   = L{1}.pQ;
 plot(fp,[jpp;jttp;jpp+jttp]/1e6)
 ylim([-4 6]),xlabel('\rho'), ylabel('j_\phi [MA/m2]')
 title(sprintf('X=%3.1f dz=%+4.1fmm',LY{1}.chi(l),LY{1}.dz(l)*1000))
 subplot(2,3,3)
 jpp  = ly.PpQ' *ly.rA    *2*pi;
 jttp = ly.TTpQ'/ly.rA/mu0*2*pi;
 plot(fp,[jpp;jttp;jpp+jttp]/1e6)
 ylim([-4 6]),xlabel('\rho'), ylabel('j_\phi [MA/m2]')
 title(sprintf('X=%3.1f dz=%+4.1fmm',ly.chi,ly.dz*1000))
 legend({'p''' 'TT''' 'tot'})
 
 kfig = kfig+1; figure(kfig)
 ly = liut(L{3},meqxk(LX,LX.shot==5231),'debug',1);
 
 for S = {'S01' 'S02'}
  T = readtable(['/Lac8_D/CREATE0883/equil/Dati4SPC/Scenario' S{1} '.xlsx'],'Sheet','CREATEdata');
  D.shot = T.Number;
  D.aW = zeros(25,length(D.shot));
  for k = 1:6
   D.aW(k,:) = T.(num2str(k,'cg%d_m_'));
  end
  for k = 1:19
   D.aW(k+6,:) = T.(num2str(k,'g%d_m_'));
  end
  
  D.aW = D.aW(g,:);
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   LY{kr}.aW(:,kY) = LY{kr}.aW(:,kY) .* (D.aW(:,kD) ~= 0);
  end
  
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   kfig = kfig+1; figure(kfig)
   for k = 1:length(g)
    subplot(7,3,k)
    plot(LY{kr}.dz(kY),LY{kr}.aW(k,kY)-D.aW(k,kD),'.')
    set(gca,'xlimspec','tight','xgrid','on','ylim',[-.2 .2])
    ylabel(gtxt(k,:))
   end
   subplot(7,3,2)
   title(sprintf('%s %dx%d',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   kfig = kfig+1; figure(kfig)
   for k = 1:length(g)
    subplot(7,3,k)
    semilogx(LY{kr}.chi(kY),LY{kr}.aW(k,kY)-D.aW(k,kD),'.')
    set(gca,'xlimspec','tight','xtick',[.01 .1 1],'xgrid','on','ylim',[-.2 .2])
    ylabel(gtxt(k,:))
   end
   subplot(7,3,2)
   title(sprintf('%s %dx%d',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  bin = -.2:.01:.2;
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   kfig = kfig+1; figure(kfig)
   l = LY{kr}.chi(kY) > 1;
   for k = 1:length(g)
    subplot(7,3,k)
    set(histogram(LY{kr}.aW(k,kY(l))-D.aW(k,kD(l)),bin),'facecolor','r')
    hold on
    set(histogram(LY{kr}.aW(k,kY(~l))-D.aW(k,kD(~l)),bin),'facecolor','g')
    hold off
    set(gca,'xgrid','on')
    ylabel(gtxt(k,:))
   end
   subplot(7,3,2)
   title(sprintf('%s %dx%d',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  kfig = kfig+1; figure(kfig)
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   l = abs(LY{kr}.dz(kY)) < 0.005;
   subplot(3,1,kr)
   boxplot(LY{kr}.aW(:,kY(l))'-D.aW(:,kD(l))',gtxt,'datalim',[-.05 .05],'extrememode','clip')
   title(sprintf('%s %dx%d dz<5mm',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  kfig = kfig+1; figure(kfig)
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   l = LY{kr}.chi(kY) < 1;
   subplot(3,1,kr)
   boxplot(LY{kr}.aW(:,kY(l))'-D.aW(:,kD(l))',gtxt,'datalim',[-.05 .05],'extrememode','clip')
   title(sprintf('%s %dx%d X<1',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  kfig = kfig+1; figure(kfig)
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   l = abs(LY{kr}.dz(kY)) < Inf;
   subplot(3,1,kr)
   boxplot(LY{kr}.aW(:,kY(l))'-D.aW(:,kD(l))',gtxt,'datalim',[-.05 .05],'extrememode','clip')
   title(sprintf('%s %dx%d All',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  
 end
 myprint
 
 %% End
 verifyEqual(test,true,true)
end

%% Gap comparison for GRT47 equilibria
function testgrt47dz(test)
 
 %% Setup
 [~,equil,cnl] = createf;
 k = equil > 0 & equil < 1000 & cnl == 0;
 equil = equil(k); cnl = cnl(k);
 LX = liuxcreate(equil,cnl);
 LL{1} = liu('create',equil(1),cnl(1),'selfm','55\.A[ABDEL].*');
 LL{2} = liu('create',equil(1),cnl(1),'zfct',@liuz0,'stabz',0,'stab',.8,'stabcint',.15,'stabint',0,'stabder',0,'itera',1000,'selfm','55\.A[ABDEL].*');
 
 %% Calculation
 for k = 1:length(LL)
  L = LL{k};
  LY{k} = liut(LL{k},LX);
  for kt = 1:length(LY{k}.t)
   rA = LY{k}.rA(kt);
   zA = LY{k}.zA(kt);
   FA = LY{k}.FA(kt);
   dr2FA = LY{k}.dr2FA(kt);
   dz2FA = LY{k}.dz2FA(kt);
   drzFA = LY{k}.drzFA(kt);
   FB = LY{k}.FB(kt);
   rB = LY{k}.rB(kt);
   lX = LY{k}.lX(kt);
   Ip = LY{k}.Ip(kt);
   rBt = LY{k}.rBt(kt);
   FBA = FB-FA;
   irBt = 1/rBt;
   Fx = LY{k}.Fx(:,:,kt);
   ag = LY{k}.ag(:,kt);
   [fPg,fTg     ] = L.P.bfct(0,L.P.bfp                        );
   [ gQg, IgQg          ] = L.P.bfct(2,L.P.bfp,L.pQ(:).^2,FA,FB);
   [aPpg,aTTpg,aPg,ahqTg] = L.P.bfct(3,L.P.bfp,        ag,FA,FB,fPg,fTg,L.idsx);
   PpQg = gQg.* aPpg';  PpQ =  sum( PpQg,2);   PQ = IgQg*  aPg;
   TTpQg = gQg.*aTTpg'; TTpQ =  sum(TTpQg,2); hqTQ = IgQg*ahqTg;
   TQ =    rBt +      irBt*hqTQ; % small diamagnetic approximation
   iTQ = 2*irBt - irBt*irBt*  TQ;
   aq = repmat(min([rA-L.G.rx(1) L.G.rx(end)-rA zA-L.G.zx(1) L.G.zx(end)-zA])*L.pinit*L.pq,L.noq,1);
   aW = zeros(length(L.G.aW),1);
   for ki = 1:200
    [aq,rq,irq] = rtcimex(aq,rA,(rA-L.G.rx(1))*L.idrx,(zA-L.G.zx(1))*L.idzx,Fx,FB-L.fq*FBA,L.idrcosq,L.idzsinq,L.cosq,false,L.drx);
    aW          = rtci   (aW,[],L.erW                ,L.ezW                ,Fx,FB         ,L.idrcosW,L.idzsinW,[]    ,false,L.drx);
    aq = abs(aq);
   end
   qaq = aq.^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)
   [Q0Q,Q1Q,Q2Q,Q3Q,Q4Q,iqQ] = fsgimex(M1q,M2q,rq,irq,rA,FA,FB,sign(Ip)*sqrt(dr2FA*dz2FA-drzFA.^2),lX,rB,iTQ,1/L.doq);
   LY{k}.aq(:,:,kt) = aq;
   LY{k}.aW(:,kt) = aW;
   LY{k}.iqQ(:,kt) = iqQ;
  end
 end
 L = LL;
 
 %  save /Lac8_D/CREATE0883/SPCoutput/efit47.mat L LX LY
 save(mymatfile,'L','LX','LY')
 
 %% Plot
 g = [1:2 7:25];
 gtxt = [num2str((1:2)','cg%d') ; num2str((1:19)','g%d')];
 gtxt([17 14 12 5],:) = num2str((3:6)','cg%d');
 for kr = 1:length(LY)
  LY{kr}.aW(2,:) = L{kr}.G.aW(2) - LY{kr}.aW(2,:);
  LY{kr}.aW = LY{kr}.aW(g,:);
 end
 close all, kfig = 0;
 
 kfig = kfig+1; figure(kfig)
 l = LY{1}.shot < 1000;
 subplot(2,2,1)
 semilogy(LY{1}.dz(l),LY{1}.chi(l),'.')
 xlim([min(LY{1}.dz ) max(LY{1}.dz )]), xlabel('dz [m]')
 ylim([min(LY{1}.chi) max(LY{1}.chi)]), ylabel('\chi')
 grid on, title('S01')
  l = LY{2}.shot < 1000;
subplot(2,2,2)
 semilogy(LY{2}.dz(l),LY{2}.chi(l),'.')
 xlim([min(LY{2}.dz ) max(LY{2}.dz )]), xlabel('dz [m]')
 ylim([min(LY{2}.chi) max(LY{2}.chi)]), ylabel('\chi')
 grid on, title('S01')
 
 kfig = kfig+1; figure(kfig)
 equil = 435;
 lx = liuxcreate(equil,0);
 g1 = lx.PpQ*mu0;
 g2 = lx.TTpQ;
 g3 = zeros(size(g1));
 [Ig1,Ig2,Ig3] = bfprmex(g1,g2,g3);
 bfp = struct('gNg',[g1 g2 g3],'IgNg',[Ig1 Ig2 Ig3],'fPg',[1;0;0],'fTg',[0;1;0]);
 [~,~,ly] = liu('create',equil,0,'bfct',@bf3imex,'bfp',bfp);
 subplot(2,3,1)
 jpp  = lx.PpQ' *lx.rA    ;
 jttp = lx.TTpQ'/lx.rA/mu0;
 fp   = linspace(0,1,length(jpp));
 plot(sqrt(fp),[jpp;jttp;jpp+jttp]/1e6)
 ylim([-4 6]),xlabel('\rho'), ylabel('j_\phi [MA/m2]')
 title(sprintf('CREATE #%d',equil))
 subplot(2,3,2)
 l = find(LY{1}.shot == equil);
 jpp  = LY{1}.PpQ(:,l)' *LY{1}.rA(l)    *2*pi;
 jttp = LY{1}.TTpQ(:,l)'/LY{1}.rA(l)/mu0*2*pi;
 fp   = L{1}.pQ;
 plot(fp,[jpp;jttp;jpp+jttp]/1e6)
 ylim([-4 6]),xlabel('\rho'), ylabel('j_\phi [MA/m2]')
 title(sprintf('X=%3.1f dz=%+4.1fmm',LY{1}.chi(l),LY{1}.dz(l)*1000))
 subplot(2,3,3)
 jpp  = ly.PpQ' *ly.rA    *2*pi;
 jttp = ly.TTpQ'/ly.rA/mu0*2*pi;
 plot(fp,[jpp;jttp;jpp+jttp]/1e6)
 ylim([-4 6]),xlabel('\rho'), ylabel('j_\phi [MA/m2]')
 title(sprintf('X=%3.1f dz=%+4.1fmm',ly.chi,ly.dz*1000))
 legend({'p''' 'TT''' 'tot'})
 
 kfig = kfig+1; figure(kfig)
 ly = liut(L{3},meqxk(LX,LX.shot==5231),'debug',1);
 
 for S = {'S01' 'S02'}
  T = readtable(['/Lac8_D/CREATE0883/equil/Dati4SPC/Scenario' S{1} '.xlsx'],'Sheet','CREATEdata');
  D.shot = T.Number;
  D.aW = zeros(25,length(D.shot));
  for k = 1:6
   D.aW(k,:) = T.(num2str(k,'cg%d_m_'));
  end
  for k = 1:19
   D.aW(k+6,:) = T.(num2str(k,'g%d_m_'));
  end
  
  D.aW = D.aW(g,:);
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   LY{kr}.aW(:,kY) = LY{kr}.aW(:,kY) .* (D.aW(:,kD) ~= 0);
  end
  
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   kfig = kfig+1; figure(kfig)
   for k = 1:length(g)
    subplot(7,3,k)
    plot(LY{kr}.dz(kY),LY{kr}.aW(k,kY)-D.aW(k,kD),'.')
    set(gca,'xlimspec','tight','xgrid','on','ylim',[-.2 .2])
    ylabel(gtxt(k,:))
   end
   subplot(7,3,2)
   title(sprintf('%s %dx%d',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   kfig = kfig+1; figure(kfig)
   for k = 1:length(g)
    subplot(7,3,k)
    semilogx(LY{kr}.chi(kY),LY{kr}.aW(k,kY)-D.aW(k,kD),'.')
    set(gca,'xlimspec','tight','xtick',[.01 .1 1],'xgrid','on','ylim',[-.2 .2])
    ylabel(gtxt(k,:))
   end
   subplot(7,3,2)
   title(sprintf('%s %dx%d',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  bin = -.2:.01:.2;
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   kfig = kfig+1; figure(kfig)
   l = LY{kr}.chi(kY) > 1;
   for k = 1:length(g)
    subplot(7,3,k)
    set(histogram(LY{kr}.aW(k,kY(l))-D.aW(k,kD(l)),bin),'facecolor','r')
    hold on
    set(histogram(LY{kr}.aW(k,kY(~l))-D.aW(k,kD(~l)),bin),'facecolor','g')
    hold off
    set(gca,'xgrid','on')
    ylabel(gtxt(k,:))
   end
   subplot(7,3,2)
   title(sprintf('%s %dx%d',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  kfig = kfig+1; figure(kfig)
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   l = abs(LY{kr}.dz(kY)) < 0.005;
   subplot(3,1,kr)
   boxplot(LY{kr}.aW(:,kY(l))'-D.aW(:,kD(l))',gtxt,'datalim',[-.05 .05],'extrememode','clip')
   title(sprintf('%s %dx%d dz<5mm',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  kfig = kfig+1; figure(kfig)
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   l = LY{kr}.chi(kY) < 1;
   subplot(3,1,kr)
   boxplot(LY{kr}.aW(:,kY(l))'-D.aW(:,kD(l))',gtxt,'datalim',[-.05 .05],'extrememode','clip')
   title(sprintf('%s %dx%d X<1',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  kfig = kfig+1; figure(kfig)
  for kr = 1:length(LY)
   [~,kD,kY] = intersect(D.shot,LY{kr}.shot);
   l = abs(LY{kr}.dz(kY)) < Inf;
   subplot(3,1,kr)
   boxplot(LY{kr}.aW(:,kY(l))'-D.aW(:,kD(l))',gtxt,'datalim',[-.05 .05],'extrememode','clip')
   title(sprintf('%s %dx%d All',S{1},L{kr}.nrx,L{kr}.nzx))
  end
  
 end
 myprint
 
 %% End
 verifyEqual(test,true,true)
end

%% All GRT47 equilibria
function grt47test(test)
 [~,equil,cnl] = createf;
 k = equil > 0;
 equil = equil(k); cnl = cnl(k);
 L = liu('create',equil(1),cnl(1),'selfm','55\.A[ABDEL].*');
 LX = liuxcreate(equil,cnl);
 LY = liut(L,LX);
 verifyEqual(test,isempty(LY),false)
 
 save testcreate/grt47test L LX LY
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(nt,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

%% Gaps convergence
function testgapconvergence(test)
 eq = 1;
 L = liu('create',eq,0,'bfct',@bfabmex,'elomin',1.2,'itert',20);
 LX = liuxcreate(eq,0);
 aW = [];
 for kiter = 20:-1:1
  LY = liut(L,LX,'itert',kiter);
  aW(:,kiter) = LY.aW;
 end
 aWW = aW;
 aWW(:,10:end) = 10000*aWW(:,10:end) - 9999*aWW(:,end);
 aWW(:,1:9) = NaN;
 close all, figure(1)
 LY = liut(L,LX,'itert',20,'debug',1);
 subplot(1,2,2)
 plot(1:20,aWW(1:6,:),'-k',1:20,aW(1:6,:)','-o'), xlabel('# iterations')
 title('Gaps cg1-6'), legend('x1e4','location','best')
 myprint
 verifyEqual(test,0,0)
end

%% FEM definition
function testfemdef(test)
 eq = 1;
 L = liu('create',eq,0);
 LX = liuxcreate(eq,0);
 close all, figure(1)
 LY = liut(L,LX,'itera',1,'debug',1)
 subplot(1,2,2)
 for k = 1:L.nh
  contour(L.ry,L.zy,reshape(L.Tyh(:,k),L.nzy,[]),linspace(0,1,11))
  hold on
 end
 plot(L.G.rl,L.G.zl,'.m')
 axis equal off
 myprint
 verifyEqual(test,0,0)
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,0,'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;
 
 clf
 subplot(221)
 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,0,'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)
 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)
 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)
 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')
 myprint
 
 verifyEqual(test,0,0)
end

%% Impact of rz grid on magnetic axis and X point location
function testasxy(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,0);
 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
 
 clf
 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')
 myprint
 
 verifyEqual(test,0,0)
end

%% Impact of rz grid on LCFS location
function testrtci(test)
 eq = 1;
 L = liu('create',eq,0);
 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',eq,0,'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));
 
 clf
 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')
 myprint
 
 verifyEqual(test,0,0)
end

%% CREATE dz feedback
function testdz(test)
 [~,equil,cnl] = createf;
 k = equil>0 & cnl==0;
 cnl = cnl(k); equil = equil(k);
 LX = liuxcreate(equil,cnl);
 L = liu('create',equil(1),cnl(1));
 PI = [.5 0; 1 0; 1.5 0; 1 .5; 1 1; 1 1.5];
 LY = cell(0);
 n = size(PI,1);
 for k = n:-1:1
  LY{k} = liut(L,LX,'stab',PI(k,1),'stabint',PI(k,2));
  equil = intersect(equil,LY{k}.shot);
 end
 [zA,dz] = deal(zeros(n,length(equil)));
 for k = 1:n
  [~,l] = intersect(LY{k}.shot,equil);
  zA(k,:) = LY{k}.zA(l);
  dz(k,:) = LY{k}.dz(l);
 end
 [z0,alpha] = deal(zeros(length(equil),1));
 for k = 1:length(equil)
  x = [zA(:,k) ones(n,1)]\dz(:,k);
  z0(k) = -x(2)/x(1);
  alpha(k) = x(1);
 end
 beta = (zA(2,:)'-z0)\dz(2,:)';
 
 save testcreate/testdz PI L LY
 
 close all
 subplot(2,2,1)
 plot(zA,dz,'-',zA(2,:),dz(2,:),'.'), xlabel('zA'), ylabel('dz')
 subplot(2,2,2)
 plot(zA,dz,'-'), xlabel('zA'), ylabel('dz'), ylim([-.001 .001])
 subplot(2,2,3)
 plot(z0,alpha,'.'), xlabel('z0'), ylabel('slope')
 subplot(2,2,4)
 plot(zA-z0',dz,'-',zA(2,:)-z0',dz(2,:),'.'), xlabel('zA-z0'), ylabel('dz')
 title(sprintf('P=1,I=0,dz=%4.2f(zA-z0)',beta))
 
 
 myprint
 
end

function y=rhox(rho)
 y = 1-(-log(1-rho.^2)).^-.5;
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 myprint
 f = dbstack; f = [f(2).file(1:end-2) '/' f(2).name];
 hfig = get(0,'children');
 for k = 1:length(hfig)
  fk = sprintf('%s-%02d',f,hfig(k).Number);
  set(hfig(k),'paperorientation','portrait','paperunits','centimeters','papersize',[17 25],'paperposition',[0 0 17 25])
  text(axes(hfig(k),'position',[0 0 1 1],'nextplot','add','visible','off'),0,0,fk,'verticalalignment','bottom')
  print(hfig(k).Number,'-dpdf',fk)
 end
end

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

function D = creategaps(dim)
 T = [readtable('/Lac8_D/CREATE0883/equil/Dati4SPC/ScenarioS01.xlsx','Sheet','CREATEdata');
      readtable('/Lac8_D/CREATE0883/equil/Dati4SPC/ScenarioS02.xlsx','Sheet','CREATEdata')];
 D.shot = T.Number;
 D.aW = zeros(length(dim),length(D.shot));
 for k = 1:length(dim)
  D.aW(k,:) = T.([deblank(dim{k}) '_m_']);
 end
 k = D.aW(1,:) == 0 | D.aW(2,:) == 0;
 D.aW([1 2],k) = NaN;
 D.aW(D.aW == 0) = NaN;
end


