function [Forces, Continuity]=FBspline(M,Forces, Continuity, Density,norm,log,parper,contourscale, zlims, rlims)
%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

it=Forces.it;
switch nargin
    case 4
        parper=false;
        contourscale=0.6;
        norm=false;
        log=false;
        zlims=[-inf inf];
    case 5
        parper=false;
        log=false;
        contourscale=0.6;
        zlims=[-inf inf];
    case 6
        parper=false;
        contourscale=0.6;
        zlims=[-inf inf];
    case 7
        contourscale=0.6;
        zlims=[-inf inf];
    case 8
        zlims=[-inf inf];
    case 9
        rlims=[];
    case 10
end
if nargin<10
    rlims=[];
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.inertforcer(:)),max(Forces.Pressforcer(:))]);
Fmaxz=max([max(Forces.Eforcez(:)),max(Forces.Bforcez(:)),max(Forces.Pressforcez(:)),max(Forces.inertforcez(:))]);
FEr=mean(Forces.Eforcer,3);
FBr=mean(Forces.Bforcer,3);
FPr=mean(Forces.Pressforcer,3);
FIr=mean(Forces.inertforcer,3);
Fdr=mean(Forces.Dragforcer,3);

FBthet=mean(Forces.Bforcethet,3);
FPthet=mean(Forces.Pressforcethet,3);
FIthet=mean(Forces.inertforcethet,3);
Fdthet=mean(Forces.Dragforcethet,3);

FEz=mean(Forces.Eforcez,3);
FBz=mean(Forces.Bforcez,3);
FPz=mean(Forces.Pressforcez,3);
FIz=mean(Forces.inertforcez,3);
Fdz=mean(Forces.Dragforcez,3);

durdt=mean(Forces.durdt,3);
duthetdt=mean(Forces.duthetdt,3);
duzdt=mean(Forces.duzdt,3);
totforcer=(FEr+FBr+FPr+FIr+Fdr-durdt);
totforcethet=(FBthet+FPthet+FIthet+Fdthet-duthetdt);
totforcez=(FEz+FBz+FPz+FIz+Fdz-duzdt);
if log
    Er=Forces.Eforcer(Forces.Eforcer~=0);
    Br=Forces.Bforcer(Forces.Bforcer~=0);
    inertr=Forces.inertforcer(Forces.inertforcer~=0);
    pressr=Forces.Pressforcer(Forces.Pressforcer~=0);
    Fminr=min([min(abs(Er(:))),min(abs(Br(:))),min(abs(inertr(:))),min(abs(pressr(:)))]);
else
    Fminr=min([min(Forces.Eforcer(:)),min(Forces.Bforcer(:)),min(Forces.inertforcer(:)),min(Forces.Pressforcer(:))]);
end
rgridmax=sum(M.nnr(1:2));

if parper
    costhet=(M.Br./M.B)';
    sinthet=(M.Bz./M.B)';
    FErd=FEr.*sinthet-FEz.*costhet;
    FBrd=FBr.*sinthet-FBz.*costhet;
    FPrd=FPr.*sinthet-FPz.*costhet;
    FIrd=FIr.*sinthet-FIz.*costhet;
    Fdrd=Fdr.*sinthet-Fdz.*costhet;
    durd=durdt.*sinthet-duzdt.*costhet;
    
    FBthetd=FBthet;
    FPthetd=FPthet;
    FIthetd=FIthet;
    Fdthetd=Fdthet;
    duthetd=duthetdt;
    
    FEzd=FEr.*costhet+FEz.*sinthet;
    FBzd=FBr.*costhet+FBz.*sinthet;
    FPzd=FPr.*costhet+FPz.*sinthet;
    FIzd=FIr.*costhet+FIz.*sinthet;
    Fdzd=Fdr.*costhet+Fdz.*sinthet;
    duzd=durdt.*costhet+duzdt.*sinthet;
    
    totforcerd=totforcer.*sinthet-totforcez.*costhet;
    totforcethetd=totforcethet;
    totforcezd=totforcer.*costhet+totforcez.*sinthet;
else
    FErd=FEr;
    FBrd=FBr;
    FPrd=FPr;
    FIrd=FIr;
    Fdrd=Fdr;
    durd=durdt;
    
    FBthetd=FBthet;
    FPthetd=FPthet;
    FIthetd=FIthet;
    Fdthetd=Fdthet;
    duthetd=duthetdt;
    
    FEzd=FEz;
    FBzd=FBz;
    FPzd=FPz;
    FIzd=FIz;
    Fdzd=Fdz;
    duzd=duzdt;
    
    totforcerd=totforcer;
    totforcethetd=totforcethet;
    totforcezd=totforcez;
end

    

%% Radial forces
f=figure('Name','radial');

clims=[-inf inf];
ax0=subplot(6,1,1);
plotvalue(ax0,durd,'Mass flux time derivative', 'mn\partialu_{r}/\partialt [Nm^{-3}]',true,clims,zlims,rlims);

ax1=subplot(6,1,2);
plotvalue(ax1,FErd,'Electric force density', 'F_{E,r} [Nm^{-3}]',true,clims,zlims,rlims);

ax2=subplot(6,1,3);
plotvalue(ax2,FBrd,'Magnetic force density','F_{B,r} [Nm^{-3}]',true,clims,zlims,rlims);

ax3=subplot(6,1,4);
plotvalue(ax3,FIrd,'Inertial force density','F_{i,r} [Nm^{-3}]',true,clims,zlims,rlims);

ax4=subplot(6,1,5);
plotvalue(ax4,FPrd,'Pressure force density','F_{p,r} [Nm^{-3}]',true,clims,zlims,rlims);

ax5=subplot(6,1,6);
plotvalue(ax5,Fdrd,'Drag force density','F_{d,r} [Nm^{-3}]',true,clims,zlims,rlims);



linkaxes([ax0 ax1 ax2 ax3 ax4 ax5],'xy')
if parper
    sgtitle(f,sprintf('Perp forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
sgtitle(f,sprintf('Radial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
end
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('Name','radial relative');
ax1=subplot(5,1,1);
color_limits=[-1 1]*10;
referencer=FBrd;
plotvalue(ax1,FErd./referencer,'Relative Electric force density', 'F_{Er}/F_{Br}',true,color_limits,zlims,rlims);

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

ax2=subplot(5,1,2);
plotvalue(ax2,FBrd./referencer,'Relative Magnetic force density','F_{Br}/F_{Br}',true,color_limits,zlims,rlims);

ax3=subplot(5,1,3);
plotvalue(ax3,FIrd./referencer,'Relative Inertial force density','F_i/F_{Br}',true,color_limits,zlims,rlims);

ax4=subplot(5,1,4);
plotvalue(ax4,FPrd./referencer,'Relative Pressure force density','F_{pr}/F_{Br}',true,color_limits,zlims,rlims);

ax5=subplot(5,1,5);
plotvalue(ax5,Fdrd./referencer,'Relative Drag force density','F_{dr}/F_{Br}',true,color_limits,zlims,rlims);

linkaxes([ax1 ax2 ax3 ax4 ax5],'xy')
if parper
    sgtitle(sprintf('Relative perp forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
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)))
end
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%s_relatforce_balancer',name),'-dpdf','-fillpage')
savefig(f,sprintf('%s_relatforce_balancer',name))


%% Azimuthal forces
f=figure('Name','azimuthal');

ax1=subplot(5,1,1);
plotvalue(ax1,duthetdt,'Mass flux time derivartive','mn\partialu_\theta/\partialt [Nm^{-3}]',true,clims,zlims,rlims);

ax2=subplot(5,1,2);
plotvalue(ax2,FBthetd,'Magnetic force density','F_{B,\theta} [Nm^{-3}]',true,clims,zlims,rlims);

ax3=subplot(5,1,3);
plotvalue(ax3,FIthetd,'Inertial force density','F_{i,\theta} [Nm^{-3}]',true,clims,zlims,rlims);

ax4=subplot(5,1,4);
plotvalue(ax4,FPthetd,'Pressure force density','F_{p,\theta} [Nm^{-3}]',true,clims,zlims,rlims);

ax5=subplot(5,1,5);
plotvalue(ax5,Fdthetd,'Drag force density','F_{d,\theta} [Nm^{-3}]',true,clims,zlims,rlims);

linkaxes([ax1 ax2 ax3 ax4 ax5],'xy')

sgtitle(sprintf('Azimuthal 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_balancethet',name),'-dpdf','-fillpage')
savefig(f,sprintf('%sforce_balancethet',name))

%% Axial forces
f=figure('Name','axial');
%zlim=[Fmin Fmax];
clims=[-inf inf];
%zlim=[-15 15];

ax0=subplot(6,1,1);
plotvalue(ax0,duzd,'Mass flux time derivartive','mn\partialu_z/\partialt [Nm^{-3}]',true,clims,zlims,rlims);

ax1=subplot(6,1,2);
plotvalue(ax1,FEzd,'Electric force density', 'F_{E,z} [Nm^{-3}]',true,clims,zlims,rlims);

ax2=subplot(6,1,3);
plotvalue(ax2,FBzd,'Magnetic force density','F_{B,z} [Nm^{-3}]',true,clims,zlims,rlims);

ax3=subplot(6,1,4);
plotvalue(ax3,FIzd,'Inertial force density','F_{i,z} [Nm^{-3}]',true,clims,zlims,rlims);

ax4=subplot(6,1,5);
plotvalue(ax4,FPzd,'Pressure force density','F_{p,z} [Nm^{-3}]',true,clims,zlims,rlims);

ax5=subplot(6,1,6);
plotvalue(ax5,Fdzd,'Drag force density','F_{d,z} [Nm^{-3}]',true,clims,zlims,rlims);

linkaxes([ax0 ax1 ax2 ax3 ax4 ax5],'xy')

if parper 
    sgtitle(sprintf('Parallel forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
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)))
end

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))

%% Axial relative forces
f=figure('Name','axial relative');
ax1=subplot(5,1,1);
color_limits=[-1 1]*10;
referencez=FPzd;
plotvalue(ax1,FEzd./referencez,'Relative Electric force density', 'F_{Ez}/F_{Ez}',true,color_limits,zlims,rlims);

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

ax2=subplot(5,1,2);
plotvalue(ax2,FBzd./referencez,'Relative Magnetic force density','F_{Bz}/F_{Ez}',true,color_limits,zlims,rlims);

ax3=subplot(5,1,3);
plotvalue(ax3,FIzd./referencez,'Relative Inertial force density','F_{Iz}/F_{Ez}',true,color_limits,zlims,rlims);

ax4=subplot(5,1,4);
plotvalue(ax4,FPzd./referencez,'Relative Pressure force density','F_{Pz}/F_{Ez}',true,color_limits,zlims,rlims);

ax5=subplot(5,1,5);
plotvalue(ax5,Fdzd./referencez,'Relative Drag force density','F_{dz}/F_{Ez}',true,color_limits,zlims,rlims);

linkaxes([ax1 ax2 ax3 ax4 ax5],'xy')
if parper
    sgtitle(sprintf('Relative Parallel forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
else
sgtitle(sprintf('Relative axial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens)))
end
f.PaperOrientation='portrait';
[~, name, ~] = fileparts(M.file);
f.PaperUnits='centimeters';
papsize=[20 29];
f.PaperSize=papsize;
print(f,sprintf('%s_relatforce_balancez',name),'-dpdf','-fillpage')
savefig(f,sprintf('%s_relatforce_balancez',name))

%% Forces sums
f2=figure('Name','total');
ax4=subplot(4,1,1);
maxforcerd=1;%cat(3,FErd,FBrd,FPrd,FIrd);
%maxforcerd=max(maxforcerd,[],3);
if parper
plotvalue(ax4,totforcerd./maxforcerd,'Total Perp force density','F [Nm^{-3}]',true,clims,zlims,rlims);
else
    plotvalue(ax4,totforcerd./maxforcerd,'Total radial force density','F [Nm^{-3}]',true,clims,zlims,rlims);
end

ax5=subplot(4,1,2);
maxforcethetd=1;%cat(3,FEzd,FBzd,FPzd,FIzd);
%maxforcezd=max(maxforcezd,[],3);
%zlims=[-15 15];
plotvalue(ax5,totforcethetd./maxforcethetd,'Total azimuthal force density','F [Nm^{-3}]',true,clims,zlims,rlims);


ax6=subplot(4,1,3);
maxforcezd=1;%cat(3,FEzd,FBzd,FPzd,FIzd);
%maxforcezd=max(maxforcezd,[],3);
%zlims=[-15 15];
if parper
    plotvalue(ax6,totforcezd./maxforcezd,'Total par force density','F [Nm^{-3}]',true,clims,zlims,rlims);
else
    plotvalue(ax6,totforcezd./maxforcezd,'Total axial force density','F [Nm^{-3}]',true,clims,zlims,rlims);
end
ax7=subplot(4,1,4);
surface(ax7,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(ax7,densitycontour(1,2:end),densitycontour(2,2:end),zlevel,'-.','linewidth',2,'color',contourcolor)

colormap(ax7,'parula')
xlabel(ax7,'z [m]')
ylabel(ax7,'r [m]')
title(ax7,'Density')
c = colorbar(ax7);
c.Label.String= 'n[m^{-3}]';
view(ax7,2)
linkaxes([ax4 ax5 ax6 ax7])

if (nargin>8 && ~isempty(zlims))
    xlim(ax7,zlims)
else
    xlim(ax7,[M.zgrid(1) M.zgrid(end)])
end
if (nargin>9 && ~isempty(rlims))
    ylim(ax7,rlims)
else
    if M.conformgeom
        ylim(ax7,[M.rgrid(1) M.rgrid(rgridmax)])
    else
        ylim(ax7,[M.rgrid(1) M.rgrid(end)])
    end
end

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('Name','Continuity');
it=Continuity.it;
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(sum(Continuity.ndivu,4),3),'n\nabla(u)','n\nabla(u)[m^{-3}s^{-1}]',true,zlim)

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

ax10=subplot(4,1,4);
plotvalue(ax10,mean(sum(Continuity.ndivu,4)+sum(Continuity.ugradn,4)+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, clims, zlims, rlims)
        if nargin < 5
            rm1strow=false;
        end
        if nargin>5
            cmin=clims(1);
            cmax=clims(2);
        else
            cmin=-inf;
            cmax=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]')
        if (nargin>6 && ~isempty(zlims))
            xlim(axis,zlims)
        else
            xlim(axis,[M.zgrid(1) M.zgrid(end)])
        end
        if (nargin>7 && ~isempty(rlims))
            ylim(axis,rlims)
        else
            if M.conformgeom
                ylim(axis,[M.rgrid(1) M.rgrid(rgridmax)])
            else
                ylim(axis,[M.rgrid(1) M.rgrid(end)])
            end
        end
        
        colormap(axis,'jet')
        c = colorbar(axis);
        c.Label.String=clab;
        caxis(axis,[cmin cmax])
        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