function meqrzpconvergence
%% Function to test convergence of rzp calculations
% The convergence studies ensure that all the different ways to write the
% rzp equations are effectively equivalent and converge to the same value.
% To make the convergence study, the grid of tcv has been increased.
% Running the function produces convergence plots and displays the errors
% between the present grid mesh and the most fine grid computed for each
% possible way to rewrite rzp equations.
%
% You can change the parameters below:
%   tok: tokamak type (tcv or anamak)
%   shot: number of the discharge on which to perform the convergence study.
%   t: single time. 
%   nstep: number of steps. At each step both the number of vertical and 
%     horizontal grid points are changed. nz can be only a power of 2,
%     therefore doubles each time. nr maintains the intial ratio between nz
%     (28/32). For computing time reasons, it is not suggested to go beyond
%     5 steps.
%   y_on: select if displaying the convergence study with the use of the y
%     grid (only plasma domain).
%   x_on: select if displaying the convergence study with the use of the x
%     grid (plasma domain + border).
%   nr_on: selct if plotting convergence study wrt 1./nr, with nr number of
%     radial grid points as input in meqp (see help meqp).
%   nz_on: selct if plotting convergence study wrt 1./nz, with nz number of
%     vertical grid points as input in meqp (see help meqp).
%   plotting parameters: a series of plotting options for all the graphs.
%
% To use the function, move it outside scripts (for example copy and paste
% and remove it once used) and run it
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% Parameters
% Changes these lines to customize the convergence study
tok = 'tcv'; % tcv or ana;
% Suggested shots to use. Uncomment the lines you want to use, comment the
% others. If you prefer, change shot and t on your will;
% Diverted plasma. 
shot = 61400; 
t = 0.6; 
% Limited plasma. 
% shot = 63783;
% t = 0.8;
nstep = 4; % number of convergence points. Modify this if you want more or less convergence points;
y_on = 1; x_on = 1;   % bool for selecting if you want to display the convergence studies with the x (x_on) and/or y (y_on) grid;
nr_on = 1; nz_on = 1; % bool for selecting if you want to display the convergence studies wrt 1./nr (nr_on) and/or 1./nz (nz_on);

%% Plotting parameters
% Change these parameters to change general plotting features. Will apply 
% to all plots and figures.
plotMarkerSize = 14;  % marker size;
plotMarkerWidth = 3;  % marker width;
titleFontSize = 25;         % title font size;
titleFontWeight = 'bold';   % title font weight;
titleInterpreter = 'tex';   % title text interpreter;
legendLocation = 'best';        % legend location;
legendOrientation = 'vertical'; % legend orientations;
legendFontSize = 15;            % legend font size;
legendNumColumns = 1;           % legend number of columns;
legendInterpreter = 'tex';      % legend text interpreter;
xlabelFontSize = 22;         % xlabel font size;
xlabelFontWeight = 'normal'; % xlabel font weight;
xlabelInterpreter = 'tex';   % xlabel text interpreter;

PlotStyle = {'plotMarkerSize',plotMarkerSize,'plotMarkerWidth',plotMarkerWidth,...
             'titleFontSize',titleFontSize,'titleFontWeight',titleFontWeight,'titleInterpreter',titleInterpreter,...
             'legendLocation',legendLocation,'legendOrientation',legendOrientation,'legendFontSize',legendFontSize,'legendNumColumns',legendNumColumns,'legendInterpreter',legendInterpreter,...
             'xlabelFontSize',xlabelFontSize,'xlabelFontWeight',xlabelFontWeight,'xlabelInterpreter',xlabelInterpreter};

%% Convergence study
addpath scripts
nz = 2.^(log2(16)+(1:nstep)-1);
nr = round(nz.*27./32);
iterative_cell = {repmat({'nr'},1,nstep), nr, repmat({'nz'},1,nstep), nz};

convergence_function_list = {@sum_rrzz,@sum_IxdBxe,@sum_IxdB,@sum_drdzMxe,@sum_dFdag,@sum_Brxe,@sum_dFdIe,@sum_dIxd};%,@sum_rzpF};
nfunc = numel(convergence_function_list);
structure_function = @rzp_convergence_study;
fix_cell = {tok,shot,t,'selx',''};

convergence_structure = convergence_study(convergence_function_list,structure_function,iterative_cell,fix_cell);

%% Plotting convergence study
plot_flag.sum_rrzz    = logical([0*x_on   1*y_on]);
plot_flag.sum_IxdBxe  = logical([[1 1]*x_on   [1 1]*y_on]);
plot_flag.sum_IxdB    = logical([[1 1]*x_on   [1 1]*y_on]);
plot_flag.sum_drdzMxe = logical([1 1]);
plot_flag.sum_dFdag   = logical([[1 1 1]*x_on [1 1 1]*y_on]);
plot_flag.sum_Brxe    = logical([[1 1]*x_on [1 1]*y_on]);
plot_flag.sum_dFdIe   = logical([[1 1 1]*x_on [1 1 1]*y_on]);
plot_flag.sum_dIxd    = logical([1*x_on 1*y_on]);
% plot_flag.sum_rzpF    = logical([1 1]);

titles.sum_rrzz    = {'rrx/rry','zzx/zzy'}; 
titles.sum_IxdBxe  = {'sum_IxdBrdrxe','sum_IxdBrdzxe','sum_IxdBzdrxe','sum_IxdBzdzxe'};
titles.sum_IxdB    = {'sum_IxdBrdr','sum_IxdBrdz','sum_IxdBzdr','sum_IxdBzdz'}; 
titles.sum_drdzMxe = {'sum_drdzMxe'};
titles.sum_dFdag   = {'dFzdr','dFzdz','dFrdr','dFrdz'};
titles.sum_Brxe    = {'sum_(2*pi*R*Bzxe)','sum_(2*pi*R*Brxe)'};
titles.sum_dFdIe   = {'dFrdIe','dFzdIe'};
titles.sum_dIxd    = {'Ix','dIxdr','dIxdz'};
% titles.sum_rzpF    = {'dFdIy','dFdag','dFdIe'};

legends.sum_rrzz    = {'rrx','rry';...
                       'zzx','zzy'}';  
legends.sum_IxdBxe  = {' Ix*drBrxe','-dIxdr*Brxe',' Iy*drBrye','-dIydr*Brye';...
                       ' Ix*dzBrxe','-dIxdz*Brxe',' Iy*dzBrye','-dIydz*Brye';...
                       ' Ix*drBzxe','-dIxdr*Bzxe',' Iy*drBzye','-dIydr*Bzye';...
                       ' Ix*dzBzxe','-dIxdz*Bzxe',' Iy*dzBzye','-dIydz*Bzye'}';  
legends.sum_IxdB    = {' Ix*drBr','-dIxdr*Br',' Iy*drBr','-dIydr*Br';...
                       ' Ix*dzBr','-dIxdz*Br',' Iy*dzBr','-dIydz*Br';...
                       ' Ix*drBz','-dIxdr*Bz',' Iy*drBz','-dIydr*Bz';...
                       ' Ix*dzBz','-dIxdz*Bz',' Iy*dzBz','-dIydz*Bz'}';
legends.sum_drdzMxe = {'drdzMxe','dzdrMxe'};
legends.sum_dFdag   = {'-2*pi*((R.*Ix)*drBr)',' 2*pi*((R.*dIxdr).*Br)','(Ix*drdzMx)','-2*pi*((R.*Iy)*drBr)',' 2*pi*((R.*dIydr).*Br)','(Iy*drdzMy)';...
                       '-2*pi*((R.*Ix)*dzBr)',' 2*pi*((R.*dIxdz).*Br)','(Ix*dzdzMx)','-2*pi*((R.*Iy)*dzBr)',' 2*pi*((R.*dIydz).*Br)','(Iy*dzdzMy)';...
                       ' 2*pi*((R.*Ix)*drBz)','-2*pi*((R.*dIxdr).*Bz)','(Ix*drdrMx)',' 2*pi*((R.*Iy)*drBz)','-2*pi*((R.*dIydr).*Bz)','(Iy*drdrMy)';...
                       ' 2*pi*((R.*Ix)*dzBz)','-2*pi*((R.*dIxdz).*Bz)','(Ix*dzdrMx)',' 2*pi*((R.*Iy)*dzBz)','-2*pi*((R.*dIydz).*Bz)','(Iy*dzdrMy)'}';
legends.sum_Brxe    = {' 2*pi*R*Bzxe','drMxe',' 2*pi*R*Bzye','drMye';...
                       '-2*pi*R*Brxe','dzMxe','-2*pi*R*Brye','dzMye'}';
legends.sum_dFdIe   = {' 2*pi*((R.*Ix)*Bzxe)','-Mxe''*dIxdr','drMex*Ix',' 2*pi*((R.*Iy)*Bzye)','-Mye''*dIydr','drMey*Iy';...
                       ' 2*pi*((R.*Ix)*Brxe)','-Mxe''*dIxdz','drMex*Ix',' 2*pi*((R.*Iy)*Brye)','-Mye''*dIydz','dzMey*Iy'}';
legends.sum_dIxd    = {'   Iy','   Ix';...
                       'dIydr','dIxdr';...
                       'dIydz','dIxdz'}';
% legends.sum_rzpF    = {'dFydag_slow ','dFdIy_fast';...
%                        'dFrzdag_slow','dFdag_fast';...
%                        'dFrzdIe_slow','dFdIe_fast'}'; 

xlabels = {'1./nr','1./nz'};
              
for jj = 1:nfunc
  func_plotted = func2str(convergence_function_list{jj});
  if nr_on
    plot_convergence_study(1./nr,convergence_structure.(func_plotted),'titles',titles.(func_plotted),...
      'legends',legends.(func_plotted),'plot_flag',plot_flag.(func_plotted),'xlabels',xlabels{1},'isfigure',1,PlotStyle{:});
  end
  if nz_on
    plot_convergence_study(1./nz,convergence_structure.(func_plotted),'titles',titles.(func_plotted),...
      'legends',legends.(func_plotted),'plot_flag',plot_flag.(func_plotted),'xlabels',xlabels{2},'isfigure',1,PlotStyle{:});
  end
end

%% 
[L,LX] = rzp_convergence_study(tok,shot,t);

sum_dFdag_present = sum_dFdag(L,LX);
sum_dFdIe_present = sum_dFdIe(L,LX);
% sum_rzpF_present  = sum_rzpF(L,LX) ;

rel_err_x_dFdag = calculate_rel_error(sum_dFdag_present(1:3,:),squeeze(convergence_structure.sum_dFdag(end,3,:))');
rel_err_y_dFdag = calculate_rel_error(sum_dFdag_present(4:6,:),squeeze(convergence_structure.sum_dFdag(end,6,:))');
rel_err_x_dFdIe = calculate_rel_error(sum_dFdIe_present(1:3,:),squeeze(convergence_structure.sum_dFdIe(end,3,:))');
rel_err_y_dFdIe = calculate_rel_error(sum_dFdIe_present(4:6,:),squeeze(convergence_structure.sum_dFdIe(end,6,:))');
% rel_err_rzpF    = calculate_rel_error(sum_rzpF_present        ,squeeze(convergence_structure.sum_rzpF(end,1,:))' );

fprintf(sprintf('\n dFrdr errors for x grid:\n      2*pi*sum(rrx(:).*Ix(:).*drBz): %.5f\n    -2*pi*sum(rrx(:).*dIxdr(:).*Bz): %.5f\n    sum(Ix(:).*drdrMx(:)): %.5f\n', rel_err_x_dFdag(1,3), rel_err_x_dFdag(2,3), rel_err_x_dFdag(3,3)));
fprintf(sprintf('\n dFrdz errors for x grid:\n      2*pi*sum(rrx(:).*Ix(:).*dzBz): %.5f\n    -2*pi*sum(rrx(:).*dIxdz(:).*Bz): %.5f\n    sum(Ix(:).*dzdrMx(:)): %.5f\n', rel_err_x_dFdag(1,4), rel_err_x_dFdag(2,4), rel_err_x_dFdag(3,4)));
fprintf(sprintf('\n dFzdr errors for x grid:\n     -2*pi*sum(rrx(:).*Ix(:).*drBr): %.5f\n     2*pi*sum(rrx(:).*dIxdr(:).*Br): %.5f\n    sum(Ix(:).*drdzMx(:)): %.5f\n', rel_err_x_dFdag(1,1), rel_err_x_dFdag(2,1), rel_err_x_dFdag(3,1)));
fprintf(sprintf('\n dFzdz errors for x grid:\n     -2*pi*sum(rrx(:).*Ix(:).*dzBr): %.5f\n     2*pi*sum(rrx(:).*dIxdz(:).*Br): %.5f\n    sum(Ix(:).*dzdzMx(:)): %.5f\n', rel_err_x_dFdag(1,2), rel_err_x_dFdag(2,2), rel_err_x_dFdag(3,2)));

fprintf(sprintf('\n dFrdr errors for y grid:\n      2*pi*sum(rry(:).*Iy(:).*drBz(lxy(:))): %.5f\n    -2*pi*sum(rry(:).*dIydr(:).*Bz(lxy(:))): %.5f\n    sum(Iy(:).*drdrMx(lxy(:))): %.5f\n', rel_err_y_dFdag(1,3), rel_err_y_dFdag(2,3), rel_err_y_dFdag(3,3)));
fprintf(sprintf('\n dFrdz errors for y grid:\n      2*pi*sum(rry(:).*Iy(:).*dzBz(lxy(:))): %.5f\n    -2*pi*sum(rry(:).*dIydz(:).*Bz(lxy(:))): %.5f\n    sum(Iy(:).*dzdrMx(lxy(:))): %.5f\n', rel_err_y_dFdag(1,4), rel_err_y_dFdag(2,4), rel_err_y_dFdag(3,4)));
fprintf(sprintf('\n dFzdr errors for y grid:\n     -2*pi*sum(rry(:).*Iy(:).*drBr(lxy(:))): %.5f\n     2*pi*sum(rry(:).*dIydr(:).*Br(lxy(:))): %.5f\n    sum(Iy(:).*drdzMx(lxy(:))): %.5f\n', rel_err_y_dFdag(1,1), rel_err_y_dFdag(2,1), rel_err_y_dFdag(3,1)));
fprintf(sprintf('\n dFzdz errors for y grid:\n     -2*pi*sum(rry(:).*Iy(:).*dzBr(lxy(:))): %.5f\n     2*pi*sum(rry(:).*dIydz(:).*Br(lxy(:))): %.5f\n    sum(Iy(:).*dzdzMx(lxy(:))): %.5f\n', rel_err_y_dFdag(1,2), rel_err_y_dFdag(2,2), rel_err_y_dFdag(3,2)));

fprintf(sprintf('\n dFrdIe errors for x grid:\n     2*pi*sum((rrx(:).*Ix(:))''*Bzxe): %.5f\n    -sum(dIxdr*Mxe): %.5f\n    sum(Ix*drMxe): %.5f\n', rel_err_x_dFdIe(1,1), rel_err_x_dFdIe(2,1), rel_err_x_dFdIe(3,1)));
fprintf(sprintf('\n dFrdIe errors for x grid:\n    -2*pi*sum((rrx(:).*Ix(:))''*Brxe): %.5f\n    -sum(dIxdz*Mxe): %.5f\n    sum(Ix*dzMxe): %.5f\n', rel_err_x_dFdIe(1,2), rel_err_x_dFdIe(2,2), rel_err_x_dFdIe(3,2)));

fprintf(sprintf('\n dFzdIe errors for y grid:\n     2*pi*sum((rry(:).*Iy(:))''*Bzxe(lxy(:))): %.5f\n    -sum(dIydr*Mxe(lxy(:))): %.5f\n    sum(Iy*drMxe(lxy(:))): %.5f\n', rel_err_y_dFdIe(1,1), rel_err_y_dFdIe(2,1), rel_err_y_dFdIe(3,1))); 
fprintf(sprintf('\n dFzdIe errors for y grid:\n    -2*pi*sum((rry(:).*Iy(:))''*Brxe(lxy(:))): %.5f\n    -sum(dIydz*Mxe(lxy(:))): %.5f\n    sum(Iy*dzMxe(lxy(:))): %.5f\n', rel_err_y_dFdIe(1,2), rel_err_y_dFdIe(2,2), rel_err_y_dFdIe(3,2))); 

% fprintf(sprintf('\n dFdIy errors:\n     slow version: %.5f\n     fast version: %.5f\n', rel_err_rzpF(1,1), rel_err_rzpF(2,1))); 
% fprintf(sprintf('\n dFdag errors:\n     slow version: %.5f\n     fast version: %.5f\n', rel_err_rzpF(1,2), rel_err_rzpF(2,2))); 
% fprintf(sprintf('\n dFdIe errors:\n     slow version: %.5f\n     fast version: %.5f\n', rel_err_rzpF(1,3), rel_err_rzpF(2,3))); 
end

function [L,LX] = rzp_convergence_study(tok,shot,t,varargin)
  [L,LX]  = rzp(tok,shot,t,varargin{:});
  
  % Add additional fields (required by meqrzpconvergence)
  xargs = {'drBrxa','dzBrxa','drBzxa','dzBzxa',...
           'drBrxu','dzBrxu','drBzxu','dzBzxu',...
           'dzdrMxa','dzdrMxu','drMxa','drMxu','dzMxa','dzMxu'};
         
  L.G  = meqg(L.G ,L.P ,xargs{:}); 
  L = rzpc(L.P,L.G);
end

function rel_err = calculate_rel_error(var1,var2)
  rel_err = zeros(size(var1));
  for ii = 1:size(var1,1)
    rel_err = abs(var1./var2 - 1);
  end
end

%% List of function used in the convergence study

function x_grid = sum_rrzz(L,LX) %#ok<INUSD>
  rry = L.rry; zzy = L.zzy;
  rrx = L.rrx; zzx = L.zzx;
  dS = L.drx*L.dzx;
  
  x_grid = zeros(2,2);
  
  x_grid(1,1) = sum(rrx(:))*dS;
  x_grid(2,1) = sum(rry(:))*dS;
  x_grid(1,2) = sum(zzx(:))*dS;
  x_grid(2,2) = sum(zzy(:))*dS;
end

function x_IxdBxe = sum_IxdBxe(L,LX)
  
  [~,~,~,~,Iy,Ix,dIydr,dIydz,dIxdr,dIxdz,lxy] = init_convergence(L,LX);
  
  Bzxe    = [L.G.Bzxa   L.G.Bzxu]  ;
  Brxe    = [L.G.Brxa   L.G.Brxu]  ;
  drBrxe  = [L.G.drBrxa L.G.drBrxu];
  drBzxe  = [L.G.drBzxa L.G.drBzxu];
  dzBrxe  = [L.G.dzBrxa L.G.dzBrxu];
  dzBzxe  = [L.G.dzBzxa L.G.dzBzxu];
  
  x_IxdBxe = zeros(4,4);

  x_IxdBxe(1,1) =  sum(sum(Ix(:).*drBrxe));
  x_IxdBxe(2,1) = -sum(sum(dIxdr(:).*Brxe));
  x_IxdBxe(3,1) =  sum(sum(Iy(:).*drBrxe(lxy(:),:)));
  x_IxdBxe(4,1) = -sum(sum(dIydr(:).*Brxe(lxy(:),:)));
  x_IxdBxe(1,2) =  sum(sum(Ix(:).*dzBrxe));
  x_IxdBxe(2,2) = -sum(sum(dIxdz(:).*Brxe));
  x_IxdBxe(3,2) =  sum(sum(Iy(:).*dzBrxe(lxy(:),:)));
  x_IxdBxe(4,2) = -sum(sum(dIydz(:).*Brxe(lxy(:),:)));
  x_IxdBxe(1,3) =  sum(sum(Ix(:).*drBzxe));
  x_IxdBxe(2,3) = -sum(sum(dIxdr(:).*Bzxe));
  x_IxdBxe(3,3) =  sum(sum(Iy(:).*drBzxe(lxy(:),:)));
  x_IxdBxe(4,3) = -sum(sum(dIydr(:).*Bzxe(lxy(:),:)));
  x_IxdBxe(1,4) =  sum(sum(Ix(:).*dzBzxe));
  x_IxdBxe(2,4) = -sum(sum(dIxdz(:).*Bzxe));
  x_IxdBxe(3,4) =  sum(sum(Iy(:).*dzBzxe(lxy(:),:)));
  x_IxdBxe(4,4) = -sum(sum(dIydz(:).*Bzxe(lxy(:),:)));
end

function x_IxdB = sum_IxdB(L,LX)
  
  [~,~,~,~,Iy,Ix,dIydr,dIydz,dIxdr,dIxdz,lxy,~,Bz,Br,...
  drBr,drBz,dzBr,dzBz] = init_convergence(L,LX);

  x_IxdB(1,1) =  sum(Ix(:).*drBr );
  x_IxdB(2,1) = -sum(dIxdr(:).*Br);
  x_IxdB(3,1) =  sum(Iy(:).*drBr(lxy(:)) );
  x_IxdB(4,1) = -sum(dIydr(:).*Br(lxy(:)));
  x_IxdB(1,2) =  sum(Ix(:).*dzBr );
  x_IxdB(2,2) = -sum(dIxdz(:).*Br);
  x_IxdB(3,2) =  sum(Iy(:).*dzBr(lxy(:)) );
  x_IxdB(4,2) = -sum(dIydz(:).*Br(lxy(:)));  
  x_IxdB(1,3) =  sum(Ix(:).*drBz );
  x_IxdB(2,3) = -sum(dIxdr(:).*Bz);
  x_IxdB(3,3) =  sum(Iy(:).*drBz(lxy(:)) );
  x_IxdB(4,3) = -sum(dIydr(:).*Bz(lxy(:)));
  x_IxdB(1,4) =  sum(Ix(:).*dzBz );
  x_IxdB(2,4) = -sum(dIxdz(:).*Bz); 
  x_IxdB(3,4) =  sum(Iy(:).*dzBz(lxy(:)) );
  x_IxdB(4,4) = -sum(dIydz(:).*Bz(lxy(:)));
end

function x_drdzMxe = sum_drdzMxe(L,LX) %#ok<INUSD>
  dzdrMxe = [L.G.dzdrMxa L.G.dzdrMxu];
  drdzMxe = [L.G.drdzMxa L.G.drdzMxu];

  x_drdzMxe(1,1) = sum(drdzMxe(:));
  x_drdzMxe(2,1) = sum(dzdrMxe(:));
end

function x_dFdag = sum_dFdag(L,LX)

  [rry,~,rrx,~,Iy,Ix,dIydr,dIydz,dIxdr,dIxdz,lxy,~,Bz,Br,...
  drBr,drBz,dzBr,dzBz,drdrMx,drdzMx,dzdrMx,dzdzMx] = init_convergence(L,LX);

  x_dFdag(1,1) =  -2*pi*sum(rrx(:).*Ix(:).*drBr );
  x_dFdag(2,1) =   2*pi*sum(rrx(:).*dIxdr(:).*Br);
  x_dFdag(3,1) =        sum(Ix(:).*drdzMx(:));
  x_dFdag(4,1) =  -2*pi*sum(rry(:).*Iy(:).*drBr(lxy(:)) );
  x_dFdag(5,1) =   2*pi*sum(rry(:).*dIydr(:).*Br(lxy(:)));
  x_dFdag(6,1) =        sum(Iy(:).*drdzMx(lxy(:)));
  x_dFdag(1,2) =  -2*pi*sum(rrx(:).*Ix(:).*dzBr );
  x_dFdag(2,2) =   2*pi*sum(rrx(:).*dIxdz(:).*Br);
  x_dFdag(3,2) =        sum(Ix(:).*dzdzMx(:));
  x_dFdag(4,2) =  -2*pi*sum(rry(:).*Iy(:).*dzBr(lxy(:)) );
  x_dFdag(5,2) =   2*pi*sum(rry(:).*dIydz(:).*Br(lxy(:)));
  x_dFdag(6,2) =        sum(Iy(:).*dzdzMx(lxy(:)));
  x_dFdag(1,3) =   2*pi*sum(rrx(:).*Ix(:).*drBz );
  x_dFdag(2,3) =  -2*pi*sum(rrx(:).*dIxdr(:).*Bz);
  x_dFdag(3,3) =        sum(Ix(:).*drdrMx(:));
  x_dFdag(4,3) =   2*pi*sum(rry(:).*Iy(:).*drBz(lxy(:)) );
  x_dFdag(5,3) =  -2*pi*sum(rry(:).*dIydr(:).*Bz(lxy(:)));
  x_dFdag(6,3) =        sum(Iy(:).*drdrMx(lxy(:)));
  x_dFdag(1,4) =   2*pi*sum(rrx(:).*Ix(:).*dzBz );
  x_dFdag(2,4) =  -2*pi*sum(rrx(:).*dIxdz(:).*Bz); 
  x_dFdag(3,4) =        sum(Ix(:).*dzdrMx(:));
  x_dFdag(4,4) =   2*pi*sum(rry(:).*Iy(:).*dzBz(lxy(:)) );
  x_dFdag(5,4) =  -2*pi*sum(rry(:).*dIydz(:).*Bz(lxy(:)));
  x_dFdag(6,4) =        sum(Iy(:).*dzdrMx(lxy(:)));
end

function x_Brxe = sum_Brxe(L,LX)
  [rry,~,rrx,~,~,~,~,~,~,~,lxy] = init_convergence(L,LX);
  
  Bzxe    = [L.G.Bzxa   L.G.Bzxu]  ;
  Brxe    = [L.G.Brxa   L.G.Brxu]  ;
  dzMxe   = [L.G.dzMxa L.G.dzMxu];
  drMxe   = [L.G.drMxa L.G.drMxu];
  
  x_Brxe(1,1) =   2*pi*sum((rrx(:))'*Bzxe);
  x_Brxe(2,1) =   sum(drMxe(:));
  x_Brxe(3,1) =   2*pi*sum((rry(:))'*Bzxe(lxy(:),:));
  x_Brxe(4,1) =   sum(sum(drMxe(lxy(:),:)));
  x_Brxe(1,2) =  -2*pi*sum((rrx(:))'*Brxe);
  x_Brxe(2,2) =   sum(dzMxe(:));
  x_Brxe(3,2) =  -2*pi*sum((rry(:))'*Brxe(lxy(:),:));
  x_Brxe(4,2) =   sum(sum(dzMxe(lxy(:),:)));
end

function x_dFdIe = sum_dFdIe(L,LX)
  [rry,~,rrx,~,Iy,Ix,dIydr,dIydz,dIxdr,dIxdz,lxy] = init_convergence(L,LX);

  Bzxe    = [L.G.Bzxa   L.G.Bzxu]  ;
  Brxe    = [L.G.Brxa   L.G.Brxu]  ;
  Mxe     = [L.G.Mxa    L.G.Mxu ]  ;
  dzMxe   = [L.G.dzMxa L.G.dzMxu];
  drMxe   = [L.G.drMxa L.G.drMxu];

  x_dFdIe(1,1) =  2*pi*sum((rrx(:).*Ix(:))'*Bzxe);
  x_dFdIe(2,1) = -sum(Mxe' * dIxdr(:));
  x_dFdIe(3,1) =  sum(drMxe' * Ix(:));
  x_dFdIe(4,1) =  2*pi*sum((rry(:).*Iy(:))'*Bzxe(lxy(:),:));
  x_dFdIe(5,1) = -sum(L.Mey * dIydr(:));
  x_dFdIe(6,1) =  sum(drMxe(lxy(:),:)' * Iy(:));
  x_dFdIe(1,2) = -2*pi*sum((rrx(:).*Ix(:))'*Brxe);
  x_dFdIe(2,2) = -sum(Mxe' * dIxdz(:));
  x_dFdIe(3,2) =  sum(dzMxe' * Ix(:));
  x_dFdIe(4,2) = -2*pi*sum((rry(:).*Iy(:))'*Brxe(lxy(:),:));
  x_dFdIe(5,2) = -sum((L.Mey * dIydz(:)));
  x_dFdIe(6,2) =  sum(dzMxe(lxy(:),:)' * Iy(:));
end

function x_dIxd = sum_dIxd(L,LX)
  [~,~,~,~,Iy,Ix,dIydr,dIydz,dIxdr,dIxdz,~] = init_convergence(L,LX);
  
  x_dIxd(1,1) = sum(Iy(:));
  x_dIxd(2,1) = sum(Ix(:));
  x_dIxd(1,2) = sum(dIydr(:));
  x_dIxd(2,2) = sum(dIxdr(:));
  x_dIxd(1,3) = sum(dIydz(:));
  x_dIxd(2,3) = sum(dIxdz(:));
  
end

% function x_rzpF = sum_rzpF(L,LX)
%   [~,dFdag_slow,~,dFdIe_slow,~] = rzpFlin(L,LX);
%   [dFdIy_fast,dFdag_fast,~,dFdIe_fast] = rzpFlin(L,LX);
%   
%   x_rzpF(1,1) = sum(sum([dFdag_slow(1:L.ny,1) dFdag_slow(1:L.ny,2:end)]));
%   x_rzpF(2,1) = sum(sum([dFdIy_fast(:,1) dFdIy_fast(:,2:end)]));
% %   [dIydr, dIydz] = meqIyfd(LX.Iy,L.drx,L.dzx);
% %   temp_slow = reshape(dFdag_slow(1:L.ny,2),[L.nzy,L.nry]);
% %   temp_fast = reshape(dFdIy_fast(:,2),[L.nzy,L.nry]);
% %   x_rzpF(1,1) = sum(sum(dIydr(2:end-1,2:end-1)));
% %   x_rzpF(2,1) = sum(sum(temp_fast(2:end-1,2:end-1)));
%   x_rzpF(1,2) = sum(sum(dFdag_slow(L.ny+2:L.ny+L.ng,:)));
%   x_rzpF(2,2) = sum(dFdag_fast(:));
%   x_rzpF(1,3) = sum(sum(dFdIe_slow(L.ny+2:L.ny+L.ng,:)));
%   x_rzpF(2,3) = sum(dFdIe_fast(:));
% end

function [rry,zzy,rrx,zzx,Iy,Ix,dIydr,dIydz,dIxdr,dIxdz,lxy,Ie,Bz,Br,...
  drBr,drBz,dzBr,dzBz,drdrMx,drdzMx,dzdrMx,dzdzMx] = init_convergence(L,LX)
  rry = L.rry; zzy = L.zzy;
  rrx = L.rrx; zzx = L.zzx;

  Iy = LX.Iy;
  Ix = contour_Ix(Iy,rrx,zzx);
  [dIydr, dIydz] = meqIyfd(Iy,L.drx,L.dzx);
  [dIxdr, dIxdz] = meqIyfd(Ix,L.drx,L.dzx);
  
  lxy = L.lxy;

  Ie = [LX.Ia;LX.Iu];

  Bz     = [L.G.Bzxa   L.G.Bzxu  ]*Ie;
  Br     = [L.G.Brxa   L.G.Brxu  ]*Ie;
  drBr   = [L.G.drBrxa L.G.drBrxu]*Ie;
  drBz   = [L.G.drBzxa L.G.drBzxu]*Ie;
  dzBr   = [L.G.dzBrxa L.G.dzBrxu]*Ie;
  dzBz   = [L.G.dzBzxa L.G.dzBzxu]*Ie;
  drdrMx = [L.G.drdrMxa L.G.drdrMxu]*Ie;
  drdzMx = [L.G.drdzMxa L.G.drdzMxu]*Ie;
  dzdrMx = [L.G.dzdrMxa L.G.dzdrMxu]*Ie;
  dzdzMx = [L.G.dzdzMxa L.G.dzdzMxu]*Ie;
end

function Ix = contour_Ix(Iy,rx,zx)
  %grid size
  nr = size(rx,2); 
  nz = size(zx,1);
  
  Ix = zeros(nz,nr);
  
  %adapt the Iy matrix to use x grid instead of y
  if ~all(size(Iy) == size(rx))
    Ix(2:nz-1,2:nr-1) = Iy;
  else
    Ix = Iy;
  end
end