function [Forces, Density]=ForceBalancespline(M,it,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 2
        contourscale=0.6;
        norm=false;
        log=false;
        zlims=[-inf inf];
    case 3
        log=false;
        contourscale=0.6;
        zlims=[-inf inf];
    case 4
        contourscale=0.6;
        zlims=[-inf inf];
    case 5
        zlims=[-inf inf];
    case 6
        
    otherwise
        error("Invalid number of arguments")
end

contourcolor=[0 0 0];%[255 20 147]/255;
fluiduThet=M.fluidUTHET(:,:,it);
Density.fluiduThet=fluiduThet;

N=M.N(:,:,it);
[R,~]=meshgrid(M.rgrid,M.zgrid);
Rinv=1./R';
Rinv(isinf(Rinv))=0;
Density.N=N;
maxdens=max(N(:));

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


Forces.Eforcer=-M.qe*N.*M.Er(:,:,it);
Forces.Eforcez=-M.qe*N.*M.Ez(:,:,it);
Bforcer=zeros(size(N,1),size(N,2),size(N,3));
Bforcez=zeros(size(N,1),size(N,2),size(N,3));
inertforce=zeros(size(N,1),size(N,2),size(N,3));
Pressforcer=zeros(size(N,1),size(N,2),size(N,3));
Pressforcez=zeros(size(N,1),size(N,2),size(N,3));
parfor j=1:size(N,3)
    Bforcer(:,:,j)=-M.qe*N(:,:,j).*fluiduThet(:,:,j).*M.Bz';
    Bforcez(:,:,j)=M.qe*N(:,:,j).*fluiduThet(:,:,j).*M.Br';
    inertforce(:,:,j)=M.me*N(:,:,j).*fluiduThet(:,:,j).^2.*Rinv;
    Pressforcer(:,:,j)=-( squeeze(M.Presstens.der(1,:,:,it(j),[1 0]))...
        + squeeze(M.Presstens(1,:,:,it(j)) - M.Presstens(4,:,:,it(j))).*Rinv...
        + squeeze(M.Presstens.der(3,:,:,it(j),[0 1])) );
    
    Pressforcez(:,:,j)=-( squeeze(M.Presstens.der(3,:,:,it(j),[1 0]))...
        + squeeze(M.Presstens(3,:,:,it(j))).*Rinv...
        + squeeze(M.Presstens.der(6,:,:,it(j),[0 1])) );
end
Forces.Bforcer=Bforcer;
Forces.Bforcez=Bforcez;
Forces.inertforce=inertforce;
Forces.Pressforcer=Pressforcer;
Forces.Pressforcez=Pressforcez;

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(Eforcer~=0);
    Br=Forces.Bforcer(Bforcer~=0);
    inert=Forces.inertforce(inertforce~=0);
    pressr=Forces.Pressforcer(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}]',false);
% 
% %zlim=[Fmin Fmax];
% zlim=[-inf inf];
% 
% ax2=subplot(4,1,2);
% plotvalue(ax2,mean(Forces.Bforcer,3),'Magnetic force density','F_{Br} [Nm^{-3}]',false,zlim)
% 
% ax3=subplot(4,1,3);
% plotvalue(ax3,mean(Forces.inertforce,3),'Centrifugal force density','F_{i} [Nm^{-3}]',false,zlim)
% 
% ax4=subplot(4,1,4);
% plotvalue(ax4,mean(Forces.Pressforcer,3),'Pressure force density','F_{pr} [Nm^{-3}]',false,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')
% 
% %% Radial relative forces
% f=figure();
% ax1=subplot(4,1,1);
% reference=mean(Forces.Bforcer,3);
% plotvalue(ax1,mean(Forces.Eforcer,3)./reference,'Relative Electric force density', 'F_{Er}/F_{Br}',false);
% 
% %zlim=[Fmin Fmax];
% zlim=[-inf inf];
% 
% ax2=subplot(4,1,2);
% plotvalue(ax2,mean(Forces.Bforcer,3)./reference,'Relative Magnetic force density','F_{Br}/F_{Br}',false)
% 
% ax3=subplot(4,1,3);
% plotvalue(ax3,mean(Forces.inertforce,3)./reference,'Relative Centrifugal force density','F_i/F_{Br}',false)
% 
% ax4=subplot(4,1,4);
% plotvalue(ax4,mean(Forces.Pressforcer,3)./reference,'Relative Pressure force density','F_{pr}/F_{Br}',false)
% 
% 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')
% 
% %% Axial forces
% f=figure();
% ax1=subplot(3,1,1);
% plotvalue(ax1,mean(Forces.Eforcez,3),'Electric force density', 'F_{Ez} [Nm^{-3}]',false);
% 
% zlim=[-inf inf];
% 
% ax2=subplot(3,1,2);
% plotvalue(ax2,mean(Forces.Bforcez,3),'Magnetic force density','F_{Bz} [Nm^{-3}]',false,zlim)
% 
% ax3=subplot(3,1,3);
% plotvalue(ax3,mean(Forces.Pressforcez,3),'Pressure force density','F_{pz} [Nm^{-3}]',false,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')
% 
% %% 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./Forcemaxr,'Total radial force density','F [Nm^{-3}]',false,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./Forcemaxz,'Total axial force density','F [Nm^{-3}]',false,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('Radial 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')

    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(2:end),abs(matrix(2:end,:)));
            else
                h=surface(axis,M.zgrid,M.rgrid(2:end),matrix(2: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')
        end
        title(axis,titl)
        grid(axis, 'on')
        view(axis,2)
    end
end