
%%Uses massSpectrumSimulator to generate a mass spectrum
addpath('.\src\');
newPos = '';
%set seed
rng(1)
totalAtomCount = 1000000;
totalNoiseCount = 0;
%Ion Input from excel sheet Model0.6nJ
%               Name   at%  wth  +   ++   +++ ++++
%ionInput = {  'Ni' ,'bal', 1, 0.7, 0.3,   0, 0; ...
ionInput = { 'Zr','bal', 1, 0.004, 0.995, 0.0005, 0; ...
             'ZrH', 17.18, 1, 0.08, 0.92, 0, 0; ...
             'Nb', 0.4, 1, 0.33, 0.66, 0.02, 0;...
             'H', 2, 1, 1, 0, 0, 0;...
             'H2', 0.027, 1, 1, 0, 0, 0;...
             'Sn', 0.90, 1, 0, 1, 0, 0;...
             'Fe', 0.006, 1, 1, 0, 0, 0;...
             'Al', 0.01, 1, 1, 0, 0, 0;};...
peakShape = 'norm';
minimumAbundance = 0.02;
confidenceAlpha = 0.1; 

%Use massSpecSimulator with the input above to generate a mass spectrum
[masses, peakCounts, ionNames, eleTable] = massSpecSimulator(totalAtomCount, totalNoiseCount, ionInput, peakShape, minimumAbundance);
eleTable
% calc simulated composition
if 0
    T = array2table(peakCounts);
    [itab,e]=ions2ionTable(myIonsList);
    peakCountIonSum = varfun(@sum,T,'InputVariables','peakCounts2','GroupingVariables','peakCounts4');
    simulatedElementalCounts = peakCountIonSum.sum_peakCounts2'*itab;
    simulatedElementalCounts./sum(simulatedElementalCounts)
end

posData = masses;
%generating random coordinates -50<x,y<50 and 0<Z<100
x = (20+20).*rand((totalAtomCount + totalNoiseCount),1)-20;
y = (20+20).*rand((totalAtomCount + totalNoiseCount),1)-20;
z = (70-0).*rand((totalAtomCount + totalNoiseCount),1)+0;

m = masses;
%Saves Posfile with name given at top
fileName = newPos;
savepos2(x',y',z',m',fileName);

%takes ionNames list and puts it in the ion cell format {
[myIonsList]=ionStr2ions(ionNames);

%% Defines Overlap Problem and solved the overlaps
overlapProblem = overlapProblemBuilder('ions',myIonsList);

[overlapProblem] = overlapProblem_getRangeCounts(overlapProblem, posData);

OP = [overlapProblem];

[comp_ele,rangeOutput,comp_ionic,ionNames,count_ele_APL,asRangedCounts,rangeTableDecompEle,rangeTableDecompIonic,overlapIonicRangeComp,ionAmounts]=massOverlapSolver(OP,confidenceAlpha);
%% Taken from demoscriptBulkDecomp

if confidenceAlpha>0
    [{'Element','LowerCI','Comp','UpperCI'}; OP.rangeElements num2cell(comp_ele)]
else
    [{'Element','Comp'}; OP.rangeElements num2cell(comp_ele)]
end

graphvizApp = '"C:\Program Files (x86)\Graphviz2.38\bin\dot.exe"';
allRangeCounts=0;
for o = 1:length(overlapIonicRangeComp)
    % (2)a
    % Create a new blank figure to draw on
    figure
    
    % Get overlap fit details
    s = OP.overlapTable{o};
    ionAmounts1 = ionAmounts{o}(1:end-1);
    ionicRangeComp = overlapIonicRangeComp{o}(:,2:end);
    rangeCounts1 = overlapIonicRangeComp{o}(:,1);
    allRangeCounts=allRangeCounts + sum(rangeCounts1);
    % The amount of each ion in each range:
    barValues = bsxfun(@times,s(:,2:end),ionAmounts1);
    
    % Plot a bar graph of the computed peak heights
    b = bar(s(:,1),barValues,1,'stacked');
      %James F addition to change colour map
    %b = bar(s(:,1),barValues,1,'stacked','FaceColor','flat');
    %for k = 1:size(barValues,2)
    %b(k).CData = k;
        %colormap lines %or could come up with my own colourmap%
    %end
    hold on
    % Actual measured range counts as x's
    plot(s(:,1),rangeCounts1,'kx');
    hold off
    
    % legend values
    ionHeadings = ionNames(OP.ionIndex{o});
    assert(size(barValues,2) == length(ionNames(OP.ionIndex{o})),'Error: result and heading size mismatch');
    legend(ionHeadings); % place legend of ion names
    fprintf('==========================\n')
    
    
    ions = (OP.overlapIonNames{o,1})'
     %table of solved overlaps
    fprintf('=== Summary (mass %s) ===\n','ions')
    format('long','g')
    disp([s(:,1) barValues])
    format('short')
    fprintf('==========================\n')
  
    
   %residual counts in peak
   counts = [sum(barValues,2), rangeCounts1];
   residualPerPeak = [(counts(:,1)-counts(:,2)).^2];
   ResidualPerOverlapGroup = sum(residualPerPeak);
   disp(ResidualPerOverlapGroup)

    % (2)b plot dot diagram
    overlaps2graphviz_colour2(OP.overlapIons(o),OP.overlapTable(o),[num2str(o) '.txt'],[num2str(o) '.png'],[],0.2,0.01,0.01,graphvizApp);
    try
        winopen([num2str(o) '.png']);
    end
end


