function readerplots(Sdata, species, ninst, nvars, nRuns, averages)

%close all

variance = var(Sdata);
std = sqrt(variance);

%Allows user to choose what species to display graphs for
OrderedSpecies = sort(species);
count=1;
specs=zeros(1,1);
fprintf('\n\n\nChoose species to plot: \n');
fprintf('0: ALL\n');
for n=1:length(species)
    fprintf('%d: %s\n', n, OrderedSpecies{n,1});
end
fprintf('99: DONE\n\n');
fprintf('Please input species one at a time')
fprintf('\n  (press enter after each entry).')
fprintf('\n  Enter "99" to finish.\n\n')

specs(1,count) = input('Species to plot: ');
while (specs(1,count) ~= 99) && (specs(1,count) ~= 0)
    count = count + 1;
    specs(1,count) = input('Species to plot: ');
end

m=1;
stop=0;
while specs(m) ~= 99 && stop ~=1
    if(specs(m)==0)
        SpecsToPlot = OrderedSpecies;
        for p=1:length(species)
            specs(p)=p;
        end
        stop=1;
    else
        SpecsToPlot{m,1} = OrderedSpecies{specs(m), 1};
    end
    m=m+1;
end

%Plotting Choices
fprintf('How to display?   \n')
fprintf(' 1 - One plot for each species (all four displays)  \n')
fprintf(' 2 - Four Species per plot (raw plots only)  \n')
fprintf(' 3 - Four Species per plot (error bar plots only)  \n')
fprintf(' 4 - Four Species per plot (color maps only)  \n')
fprintf(' 5 - Four Species per plot (histograms only)  \n')
fprintf(' 6 - All Plots  \n\n')
PlotChoice=input('How to display: ');
fprintf('\nShow all chosen species on a separate single plot?   \n')
PlotChoice2=input('1-Yes, 0-No:  ');

%Main Code
color = ['m' 'c' 'r' 'g' 'b' 'k'];
line = {'-' '-.' '--' ':'};
hist3d = [];                     
x = [];
paint = 1;
type = 1;
%plot results
plotcount=0;
figcount=1;
for ind=1:nvars
    for o=1:length(specs)
        if ind==specs(o)
            if plotcount == 4  %keeps track of number of plots/subplots
                figcount=figcount+1;
                figure
                plotcount=0;
            end
            plotcount=plotcount+1;
            
            if PlotChoice > 3 || PlotChoice == 1
                %sets up for histograms/colormaps
                Conc_max = max(averages(1,:,ind));
                clims = [0 nRuns];
                for t = 1:ninst
                    if Conc_max < 100
                        x=1:ninst;
                    else
                        %change step size to equal 1 if you want
                        %actual dimensions
                        step = (Conc_max)/ninst;
                        x = 1:step:Conc_max;
                    end
                    y = Sdata(1:nRuns,t,ind);
                    hist2d= hist(y,x);
                    [m,n] = size(hist2d);
                    for L = 1:n
                        hist3d(t,L) = hist2d(L);
                    end
                end
            end

            if PlotChoice==1  %plots runs, error bars, histogram, and color bars on one plot
                figure(ind)
                for runs=1:nRuns
                    hold on;
                    subplot(2,2,1), plot(Sdata(runs,:,ind))
                    title(sprintf('species %s', OrderedSpecies{ind,1}));
                end
                subplot(2,2,2), errorbar(averages(:,:,ind),std(:,:,ind))
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                subplot(2,2,3), bar3(rot90(hist3d),5,'detached')
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                subplot(2,2,4), imagesc(x,1:ninst,rot90(hist3d),clims)
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                colorbar();
            elseif PlotChoice==2  %plots runs, 4 to a plot
                figure(figcount)
                for runs=1:nRuns
                    hold on;
                    subplot(2,2,plotcount), plot(Sdata(runs,:,ind))
                    title(sprintf('species %s', OrderedSpecies{ind,1}));
                end                
            elseif PlotChoice==3  %plots errorbars, 4 to a plot
                figure(figcount)
                subplot(2,2,plotcount), errorbar(averages(:,:,ind),std(:,:,ind));
                title(sprintf('species %s', OrderedSpecies{ind,1}));
            elseif PlotChoice==4  %plots colormaps
                figure(figcount)
                subplot(2,2,plotcount), imagesc(x,1:ninst,rot90(hist3d),clims)
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                colorbar();
            elseif PlotChoice==5  %plots histograms
                figure(figcount)
                subplot(2,2,plotcount), bar3(rot90(hist3d),5,'detached')                
                title(sprintf('species %s', OrderedSpecies{ind,1}));
            elseif PlotChoice==6  %plots everything
                %plots runs
                figure(ind)
                plot(Sdata(runs,:,ind))
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                %plots error bars
                figure(ind+length(specs))
                errorbar(averages(:,:,ind),std(:,:,ind));   
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                %plots histograms
                figure(ind+2*length(specs))
                bar3(rot90(hist3d),5,'detached')
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                zlabel('Number of runs');
                %plots colormaps
                figure(ind+3*length(specs))
                imagesc(x,1:ninst,rot90(hist3d),clims)
                title(sprintf('species %s', OrderedSpecies{ind,1}));
                colorbar();
            end
            for subplots=1:4
                subplot(2,2,subplots)
                xlabel('time(min)');
                ylabel('Concentration(density)');
            end
            if PlotChoice2==1   %plots all species on the same plot
                %varies the colors and line types on the plots
                if paint > 6
                    paint = 1;
                    type = type + 1;
                end
                if type > 4
                    type = 1;
                end
                %plots the average if more then one sim is run
                figure(99)
                if nRuns < 2
                    plot(Sdata(runs,:,ind),[line{1,type} color(paint)]);
                else
                    plot(averages(:,:,ind),[line{1,type} color(paint)]);
                end
                hold on
                title('All Chosen Species')
                paint = paint + 1;
                legend(SpecsToPlot);
            end
            xlabel('time(min)');
            ylabel('Concentration(density)');
            grid on;
            hold on;
        end
    end
end
