function [Forces, Continuity]=FBspline(M,it,Forces, Continuity, Density,norm,log,contourscale, zlims)
%ForceBalance Show the radial force balance
%   Plot the three radial forces for the given time-step it or averaged
%   over the range of time steps defined in it

if it==':'
    it=floor(0.95*size(M.t2d)):size(M.t2d);
end
switch nargin
    case 5
        contourscale=0.6;
        norm=false;
        log=false;
        zlims=[-inf inf];
    case 6
        log=false;
        contourscale=0.6;
        zlims=[-inf inf];
    case 7
        contourscale=0.6;
        zlims=[-inf inf];
    case 8
        zlims=[-inf inf];
    case 9
        
    otherwise
        error("Invalid number of arguments")
end

contourcolor=[0 0 0];%[255 20 147]/255;
N=Density.N;

maxdens=max(N(:));

densitycontour=contour(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens);

Fmaxr=max([max(Forces.Eforcer(:)),max(Forces.Bforcer(:)),max(Forces.inertforce(:)),max(Forces.Pressforcer(:))]);
Fmaxz=max([max(Forces.Eforcez(:)),max(Forces.Bforcez(:)),max(Forces.Pressforcez(:))]);
if log
    Er=Forces.Eforcer(Forces.Eforcer~=0);
    Br=Forces.Bforcer(Forces.Bforcer~=0);
    inert=Forces.inertforce(Forces.inertforce~=0);
    pressr=Forces.Pressforcer(Forces.Pressforcer~=0);
    Fminr=min([min(abs(Er(:))),min(abs(Br(:))),min(abs(inert(:))),min(abs(pressr(:)))]);
else
    Fminr=min([min(Forces.Eforcer(:)),min(Forces.Bforcer(:)),min(Forces.inertforce(:)),min(Forces.Pressforcer(:))]);
end
rgridmax=sum(M.nnr(1:2));

%% Radial forces
f=figure();

ax1=subplot(4,1,1);
plotvalue(ax1,mean(Forces.Eforcer,3),'Electric force density', 'F_{Er} [Nm^{-3}]',true);

%zlim=[Fmin Fmax];
zlim=[-inf inf];

ax2=subplot(4,1,2);
plotvalue(ax2,mean(Forces.Bforcer,3),'Magnetic force density','F_{Br} [Nm^{-3}]',true,zlim)

ax3=subplot(4,1,3);
plotvalue(ax3,mean(Forces.inertforce,3),'Centrifugal force density','F_{i} [Nm^{-3}]',true,zlim)

ax4=subplot(4,1,4);
plotvalue(ax4,mean(Forces.Pressforcer,3),'Pressure force density','F_{pr} [Nm^{-3}]',true,zlim)

linkaxes([ax1 ax2 ax3 ax4],'xy')
sgtitle(sprintf('Radial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))

f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%sforce_balancer',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sforce_balancer',name))

%% Radial relative forces
f=figure();
ax1=subplot(4,1,1);
color_limits=[-1 1];
reference=Forces.Bforcer;
plotvalue(ax1,mean(Forces.Eforcer./reference,3),'Relative Electric force density', 'F_{Er}/F_{Br}',true,color_limits);

%zlim=[Fmin Fmax];
zlim=[-inf inf];

ax2=subplot(4,1,2);
plotvalue(ax2,mean(Forces.Bforcer./reference,3),'Relative Magnetic force density','F_{Br}/F_{Br}',true,color_limits)

ax3=subplot(4,1,3);
plotvalue(ax3,mean(Forces.inertforce./reference,3),'Relative Centrifugal force density','F_i/F_{Br}',true,color_limits)

ax4=subplot(4,1,4);
plotvalue(ax4,mean(Forces.Pressforcer./reference,3),'Relative Pressure force density','F_{pr}/F_{Br}',true,color_limits)

linkaxes([ax1 ax2 ax3 ax4],'xy')
sgtitle(sprintf('Relative radial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))

f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%s_relatforce_balance',name),'-dpdf','-fillpage')
savefig(f,sprintf('%s_relatforce_balance',name))

%% Axial forces
f=figure();
ax1=subplot(3,1,1);
plotvalue(ax1,mean(Forces.Eforcez,3),'Electric force density', 'F_{Ez} [Nm^{-3}]',true);

zlim=[-inf inf];

ax2=subplot(3,1,2);
plotvalue(ax2,mean(Forces.Bforcez,3),'Magnetic force density','F_{Bz} [Nm^{-3}]',true,zlim)

ax3=subplot(3,1,3);
plotvalue(ax3,mean(Forces.Pressforcez,3),'Pressure force density','F_{pz} [Nm^{-3}]',true,zlim)

linkaxes([ax1 ax2 ax3],'xy')
sgtitle(sprintf('Axial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))

f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%sforce_balancez',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sforce_balancez',name))

%% Forces sums
f2=figure();
ax4=subplot(3,1,1);
totforcer=mean(Forces.Eforcer+Forces.Bforcer+Forces.inertforce+Forces.Pressforcer,3);
Forcemaxr=cat(3,mean(Forces.Eforcer,3),mean(Forces.Bforcer,3),mean(Forces.inertforce,3),mean(Forces.Pressforcer,3));
Forcemaxr=max(Forcemaxr,[],3);
plotvalue(ax4,totforcer,'Total radial force density','F [Nm^{-3}]',true,zlims)

ax5=subplot(3,1,2);
totforcez=mean(Forces.Eforcez+Forces.Bforcez+Forces.Pressforcez,3);
Forcemaxz=cat(3,mean(Forces.Eforcez,3),mean(Forces.Bforcez,3),mean(Forces.Pressforcez,3));
Forcemaxz=max(Forcemaxz,[],3);
plotvalue(ax5,totforcez,'Total axial force density','F [Nm^{-3}]',true,zlims)

ax6=subplot(3,1,3);
surface(ax6,M.zgrid,M.rgrid,mean(N,3),'edgecolor','none');
hold on;
zlevel=interp2(M.zgrid, M.rgrid, mean(N,3), densitycontour(1,2:end), densitycontour(2,2:end));
plot3(ax6,densitycontour(1,2:end),densitycontour(2,2:end),zlevel,'-.','linewidth',2,'color',contourcolor)
xlim(ax6,[M.zgrid(1) M.zgrid(end)])
ylim(ax6,[M.rgrid(1) M.rgrid(rgridmax)])
colormap(ax6,'parula')
xlabel(ax6,'z [m]')
ylabel(ax6,'r [m]')
title(ax6,'Density')
c = colorbar(ax6);
c.Label.String= 'n[m^{-3}]';
view(ax6,2)

sgtitle(sprintf('Forces sum t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f2.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f2.PaperUnits='centimeters';
papsize=[20 29];
f2.PaperSize=papsize;
print(f2,sprintf('%sforce_dens',name),'-dpdf','-fillpage')
savefig(f2,sprintf('%sforce_dens',name))


%% Continuity
if( ~isempty(Continuity))
f3=figure();
t=M.t2d(it);

ax7=subplot(4,1,1);
plotvalue(ax7,mean(Continuity.dndt,3),'\partialn/\partialt','\partialn/\partialt [m^{-3}s^{-1}]',true)


zlim=[max([min(min(mean(Continuity.ndivu(2:end,:,:),3))),min(min(mean(Continuity.ugradn(2:end,:,:),3)))]) min([max(max(mean(Continuity.ndivu(2:end,:,:),3))),max(max(mean(Continuity.ugradn(2:end,:,:),3)))])];

ax8=subplot(4,1,2);
plotvalue(ax8,mean(Continuity.ndivu,3),'n\nabla(u)','n\nabla(u)[m^{-3}s^{-1}]',true,zlim)

ax9=subplot(4,1,3);
plotvalue(ax9,mean(Continuity.ugradn,3),'u\nabla(n)','u\nabla(n)[m^{-3}s^{-1}]',true,zlim)

ax10=subplot(4,1,4);
plotvalue(ax10,mean(Continuity.ndivu+Continuity.ugradn+Continuity.dndt,3),'Total','total[m^{-3}s^{-1}]',true,zlim)

sgtitle(sprintf('Continuity t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
f3.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f3.PaperUnits='centimeters';
papsize=[20 19];
f3.PaperSize=papsize;
print(f3,sprintf('%scontinuity',name),'-dpdf','-fillpage')
savefig(f,sprintf('%scontinuity',name))
end

    function plotvalue(axis,matrix,titl, clab, rm1strow, zlims)
        if nargin < 5
            rm1strow=false;
        end
        if nargin>5
            zmin=zlims(1);
            zmax=zlims(2);
        else
            zmin=-inf;
            zmax=inf;
        end
        if rm1strow
            if log
                h=surface(axis,M.zgrid,M.rgrid(12:end),abs(matrix(12:end,:)));
            else
                h=surface(axis,M.zgrid,M.rgrid(12:end),matrix(12:end,:));
            end
        else
            if log
                h=surface(axis,M.zgrid,M.rgrid,abs(matrix));
            else
                h=surface(axis,M.zgrid,M.rgrid,matrix);
            end
        end
        set(h,'edgecolor','none');
        hold on;
        if log
            zpos=interp2(M.zgrid, M.rgrid, abs(matrix), densitycontour(1,2:end), densitycontour(2,2:end));
        else
            zpos=interp2(M.zgrid, M.rgrid, matrix, densitycontour(1,2:end), densitycontour(2,2:end));
        end
        plot3(axis,densitycontour(1,2:end),densitycontour(2,2:end),zpos,'-.','linewidth',2,'color',contourcolor)
        xlabel(axis,'z [m]')
        ylabel(axis,'r [m]')
        xlim(axis,[M.zgrid(1) M.zgrid(end)])
        ylim(axis,[M.rgrid(1) M.rgrid(rgridmax)])
        colormap(axis,'jet')
        c = colorbar(axis);
        c.Label.String=clab;
        caxis(axis,[zmin zmax])
        if norm
            caxis(axis,[Fminr Fmaxr]);
        end
        if log
            set(axis,'colorscale','log')
            caxis(axis,'auto')
        end
        title(axis,titl)
        grid(axis, 'on')
        view(axis,2)
    end
end