function frontend_getconcsweeplastsnap(varargin)
% This code plots images from last snap at various cutoff F-actin
% concentrations  from .mat files generated by readMEDYANtraj.
%% Input for the function
% varargin{1} - matfileloadpath - path to find .mat files in
% varargin{2} - filenamecell - cell consisting of filename strings that will help the user
% to load .mat files.
% For example: {'Arp1nM','Arp5nM','Arp10nM','Arp25nM','Arp50nM'};
%% Outputs
% 1. Plots require the use of tight_subplot. Please download from
% https://www.mathworks.com/matlabcentral/fileexchange/27991-tight_subplot-nh-nw-gap-marg_h-marg_w
matloadfilepath = varargin{1};
filenamecell=varargin{2};
%cutoff concentrations in micromolar
cutoff = [0 10 20 40 100];
% Dimensions of a pixel
pixelsize = [100 100 100];
disp(['Pixel size in nm ',num2str(pixelsize)]);
pixelVinL = pixelsize(1)*pixelsize(2)*pixelsize(3)*1e-27*1e3;
pixelvolumeinmum3 = pixelsize(1)*pixelsize(2)*pixelsize(3)*1e-9;
% Dimensions of reaction volume used in this study
Rxnvoldim = [2000 2000 7500];
% Number of pixels generated per image
Npixels = ceil(Rxnvoldim./pixelsize);
% Npixels(1)*Npixels(2)*Npixels(3)*64/1e9
% 100nm -> 204Kb
% 50nm -> 15.4 Mb
% 25nm -> 9.8 Mb
% 10nm -> 1.92 Gb
% 5nm -> 15.36 Gb
% 1nm -> 1920 Gb

fig = figure('units','inch','position',[1,1,20,6],'Color','w');
% [yspacing xspacing], [bottommargin, topmargin], [leftmargin, rightmargin]
[ha, pos] = tight_subplot(numel(filenamecell),numel(cutoff),[.2 .05],[.07 .1],[.07 .07]);
%


for fileid=1:numel(filenamecell)
    loadmatfilename = loadmatfilenamecell{fid};
    
    clear p;
    %Constructor for instances of class pixel
    for i = 1:Npixels(1)*Npixels(2)*Npixels(3)
        p(i) = pixel;
    end
    %Set the total number of pixels
    p.setgetNpixels(Npixels);
    disp('pixels created');
    runid = 1;
    
    %Fetch filenames of all replicates of the snapshot file
    filelist = ls([filepath,'/',loadmatfilename,'_S*.mat']);
    %If Windows
    filelisttemp={};
    if ispc
        for fid = 1:size(filelist,1)
            if(numel(filelist(fid,:))==0)
                %if filename is empty, skip.
                continue;
            end
            filelisttemp = [filelisttemp, [filepath,filelist(fid,:)]];
        end
        filelist = filelisttemp;
    end
    %Just use the first replicate run to generate the images
    for fid = 1
        %Find last snapshot. Ignore snaps with imperfect data
        for snap = 1:size(r(1).s,2)
            filcoord = r(1).s(snap).f.coord_cell1;
            if(isempty(filcoord))
                break;
            end
        end
        lastsnap = snap;
        %
        count = 1;
        for snap = lastsnap
            %Get the filament coordinate cell corresponding to the snapshot
            filcoordcell{count} = r(1).s(snap).f.coord_cell1;
            count = count + 1;
        end
        %
        resetneeded = true;
        for iter = 1:numel(snaplist)
            %Get the filament coordinate cell
            filcoord = filcoordcell{iter};
            disp(['pixellating snap ',num2str(snaplist(iter))]);
            % reset pixels
            if(resetneeded)
                %Go through all pixels, reset pixels
                for pidx = 1:numel(p)
                    %Set them as unoccupied, and unaccounted pixels, with no F-actin
                    p(pidx).occupiedstatus = false;
                    p(pidx).accountedstatus = false;
                    p(pidx).count = 0;
                    p(pidx).meanconc = 0;
                end
            end
            % Mark cubes outside the boundary as occupied.
            p = pixelcylinderinitialize(p,pixelsize);
            disp('pixel initialized');
            crosssection = Npixels(1)*Npixels(2);
            % Go through filaments in this snap
            countoccupied = 0;
            for fpos = 1:size(filcoord,1)
                % reshape the coordinate to Nx3 matrix where N is number of
                % cylinder ends
                fil = reshape(filcoord{fpos}',3,[])';
                %Interpolate each cylinder in a filament to get coordinates of
                %monomers
                Coord_cat = interpolateallmonomersinfilemant(fil);
                %Convert them to coordinates on the bin axis
                binvec = ceil(Coord_cat./pixelsize);
                %Get the 1 dimensional ID corresponding to each monomer
                onedpos = (binvec(:,3)-1).*crosssection + (binvec(:,2)-1).*Npixels(1) + binvec(:,1);
                onedpos = onedpos';
                %Go through the list of one dimensional IDs
                for x=onedpos
                    if(~p(x).occupiedstatus)
                        %Count pixels that are occupied
                        countoccupied = countoccupied + 1;
                    end
                    p(x).occupiedstatus = true;
                    p(x).accountedstatus = false;
                    %Count the number of monomers
                    p(x).count = p(x).count + 1;
                end
            end
            disp('pixel occupancy updated');
            % Each entry is 1 monomer long.
            VinL=pi*1e3*1e3*7.5e3*1e-27*1e3;
            monomerspercount = 1;
            %Factor to convert count to micromoles
            factor = monomerspercount*1e6/(6.023e23*VinL);
            %Total number of monomers
            Nmontotal = sum([p(:).count]);
            %Total concentration of F-actin
            Capprox = factor*Nmontotal;
            factorpixel = monomerspercount*1e6/(6.023e23*pixelVinL);
            % Update mean concentration of F-actin in each pixel
            for i = 1:size(p,2)
                pvecids = [p.getneighbors(i),i];
                p(i).meanconc = mean([p(pvecids).count])*factorpixel;
            end
        end
        %% Plotter section of the code
        %Get mean concentration of all pixels
        C=[p(:).meanconc];
        %Get the accounted status of all pixels
        Accounted = [p(:).accountedstatus];
        withincylinder = find(Accounted==false);
        for j = 1:numel(cutoff)
            %Find pixels with actin larger than the cutoff
            nonzeroelems = find(C>cutoff(j));
            %remove out of boundary pixels
            nonzeroelems = intersect(withincylinder,nonzeroelems);
            Coords=zeros(numel(nonzeroelems),3);
            count = 1;
            %Get 3D cooordinates of occupied pixels
            for i = nonzeroelems
                Coords(count,:) = pixelsize.*pixel.getidx3dfrom1d(i);
                count = count + 1;
            end
            Coords = Coords./1e3;
            % Plot
            
            axes(ha(plotid));plot3(Coords(:,3),Coords(:,2),Coords(:,1),'rs','LineWidth',0.1,'MarkerFaceColor','r','MarkerSize',6);
            title(num2str(cutoff(j)));
            xlim([0,7.5]);
            ylim([0,2]);
            zlim([0,2]);
            view([-10,12]);
            set(gca,'LineWidth',2);
            set(gca,'FontSize',20);
            % Plot cylinder
            plotcylinder();
        end
    end
end
%save figure
savefig(fig,'Concsweep_lastsnap.fig');
print(fig,'Concsweep_lastsnap','-dpng','-r300')
end
%% END of main function
%
%
%% ADDITIONAL AUXILLARY FUNCTIONS
%Function initializes cylinders by marking those that are out of boundary
%appropriately
function p = pixelcylinderinitialize(p,pixelsize)
cylindercenter = [1000 1000]./pixelsize(1:2);
%Radius of 1000nm
Radius = 1000./pixelsize(1);
count = 1;
for pid = 1:size(p,2)
    offsetintegers = p.getidx3dfrom1d(pid);
    if(norm(offsetintegers(1:2)-cylindercenter-[0.5 0.5]) > Radius)
        p(pid).occupiedstatus = false;
        p(pid).accountedstatus = true;
        count = count + 1;
    end
end
disp(['Number of pixels outside boundary=',num2str(count)]);
end
%Function plots the circular ends of a cylinder
function plotcylinder()
hold on;
[X,Y,Z]=cylinder(1,50);
Z=Z*0.25;
surf(Z,Y+1,X+1,'FaceAlpha',0.4,'FaceColor','y','EdgeColor','k');
surf(Z+7.5-0.25,Y+1,X+1,'FaceAlpha',0.4,'FaceColor','y','EdgeColor','k');
hold off;
end