function AllCoords_SmoothDensity(ScriptVar, AllCoords,...
    filename, ParticleTypesToBePlotted, frames_of_interest, ...
    NumBinsOrCellOfXEdgesYEdges, SlidingWindow, legend_numbers)
    % Output - density plots (along x-axis) depicting particle distribution 
    %    across different timepoints in the simulation run
    %
    % Inputs - 
    %    ScriptVar - Parameters container map
    %    AllCoordsCell - Cell containing particle coordinates,
    %    filename - string; don't append '.gif', '.fig' etc. here
    %    ParticleTypesToBePlotted - set to 1 in case of one particle type
    %       in the simulation
    %    frames_of_interest - 1xk array of frames for which density plots
    %       are to be plotted on the figure
    %    NumBinsOrCellOfXEdgesYEdges - set to 500 for smoothening along x
    %    SlidingWindow - width of the sliding window (set to 100 given 500
    %       bins along x)
    %
    % Usage - 
    %    AllCoords_SmoothDensity(ScriptVar, AllCoords, 'trytry.fig', 1,...
    %       [1 2 3 9 11], 500, 100)
    
    xLimit = ScriptVar('xLimit');
    yLimit = ScriptVar('yLimit');
    
    if ~exist('DensityOrHeatmap','var')
        DensityOrHeatmap = 'density';
    end
    
    %if length(ParticleTypesToBePlotted) == 1
    %    ColourArray = {'r','g','m'};
    %    disp('Using red-green-magenta colour order')
    %end
    
    N_Particles = ScriptVar('N_Particles');
    ParticleTypesToBePlotted = sort(ParticleTypesToBePlotted);
    AllCoords_new = []; %zeros(size(AllCoords));
    for j = 1:length(ParticleTypesToBePlotted)
        if ParticleTypesToBePlotted(j) == 1
            start_index = 1;
        else
            start_index = sum(N_Particles(1:ParticleTypesToBePlotted(j)-1))+1;
        end
        end_index = sum(N_Particles(1:ParticleTypesToBePlotted(j)));
        a = AllCoords(start_index:end_index,:,:);
        if isKey(ScriptVar, 'xyFRAPRegion')
            ind = a(:,3,:)==1;
            ind = cat(2,ind,ind,ind);
            a(ind) = NaN;
            %a = sum(~isnan(a),1);
            %a = squeeze(a(1,1,:));
            ThisParticleTypeCoords = a;  %a(a(:,3)==0,:,:);
        else
            ThisParticleTypeCoords = a;
        end
        AllCoords_new = cat(1,AllCoords_new,ThisParticleTypeCoords);
        %ThisFrameCoords = cat(1,ThisFrameCoords,ThisParticleTypeCoords);
    end
        
    if isnumeric(NumBinsOrCellOfXEdgesYEdges) == 1
        NumBins = NumBinsOrCellOfXEdgesYEdges;
        XEdges = 0:xLimit/NumBins:xLimit;
        %XEdges = XEdges(1:end-1);
        if strcmpi(DensityOrHeatmap, 'Heatmap')
            YEdges = 0:yLimit/NumBins:yLimit;
        end
    else
        XEdges = NumBinsOrCellOfXEdgesYEdges{1};
        NumBins = length(XEdges) - 1;
        if strcmpi(DensityOrHeatmap, 'Heatmap')
            YEdges = NumBinsOrCellOfXEdgesYEdges{2};
        end
    end
    
    N_Frames = ScriptVar('N_Frames'); 
    if strcmpi(DensityOrHeatmap, 'Heatmap')
        %NumBinsX = length(XEdges) - 1;
        %NumBinsY = length(YEdges) - 1;
        densitycell = zeros(length(XEdges)-1, length(YEdges)-1, N_Frames); % 7, 7, 5000
    else
        densitycell = zeros(length(XEdges)-1, N_Frames);
    end
    
    for i = 1:N_Frames
        k = AllCoords_new(:,1,i);
        if strcmpi(DensityOrHeatmap, 'Heatmap')
            k_y = AllCoords_new(:,2,i);
            [BinIndices, ~, ~] = histcounts2(k, k_y, XEdges, YEdges);
            densitycell(:,:,i) = rmmissing(rot90(BinIndices));
        else
            [BinIndices, ~] = histcounts(k,XEdges);
            densitycell(:,i) = rmmissing(BinIndices);
            densitycell(:,i) = densitycell(:,i)/N_Particles;
        end
    end
    
    figure();clf;set(gcf,'Position',[816   865   700   300]);
    movegui('center')
    if ~strcmp(frames_of_interest,'gif')
        if strcmpi(DensityOrHeatmap, 'Heatmap')
            h = heatmap(densitycell(:,:,frames_of_interest));
            title(['Frame = ' num2str(frames_of_interest)]);
            h.YDisplayLabels = flip(h.YDisplayLabels); 
        else
            [~, t] = size(frames_of_interest);
            for i = 1:t
                plot(linspace(0,1,NumBinsOrCellOfXEdgesYEdges),...
                    smoothdata(densitycell(:,frames_of_interest(i)),...
                    'gaussian',SlidingWindow),...
                    'LineWidth',2);
                hold on
            end
            %set(gca,'FontSize',15);
            legend_cell = cell(1,length(frames_of_interest));
            if ~exist('legend_numbers','var') || ~legend_numbers
                for ii = 1:length(frames_of_interest)
                    legend_cell{ii} = append("Frame = ", num2str(frames_of_interest(ii)));
                    %," (t=",num2str(time_points(frames_of_interest(1))), "s)");
                end
            elseif legend_numbers
                frames_legend = ScriptVar('CheckDensity');
                for ii = 1:length(frames_of_interest)
                    j = frames_of_interest(ii);
                    legend_cell{ii} = append("Frame = ", num2str(frames_legend(j)));
                    %," (t=",num2str(time_points(frames_of_interest(1))), "s)");
                end
            end
            legend(legend_cell,'FontSize',15,'Interpreter','latex')
            xlim([0 1]); 
            %ylim([0 1])
            %ylim([0 max(densitycell, [], 'all')+5]);
            xlabel('Bin');
            ylabel('Number of particles');
        end
    end
    %set(gca,'FontSize',15);
    set(findall(gcf,'-property','FontSize'),'FontSize',20)
    set(legend,'FontSize',15)
    xlabel('$x\ (\mu m)$', 'Interpreter', 'LaTex');
    ylabel('Particle frequency', 'Interpreter', 'LaTex');
    %xticks(1:NumBinsOrCellOfXEdgesYEdges)
    %yticks([0 0.25 0.5 0.75 1])
    set(gca,'TickLabelInterpreter', 'latex');
    savefig(strcat(filename,'.fig')) 
end