function [ header output_data ] = Plot_Atomic_Contribs_V3( folder )


%MERGE FRAGMETNS SEEMS TO NOT WORK AT THE MO. But only if merge_spd=0

%Makes a plot of the atomic contribs to a gelius spectrum when given a
%folder containing output files of Gen_Gelius_Spectrum script
%   Detailed explanation goes here


%HISTORY:
%V1 = Made sometime in 2016, not sure when.
%29/07/2016: V2 done. Only change is that it now outputs two arguments.
%           Wasnt massively worth calling it a different function over really....
%           The output data is just the stuff that gets written to the txt
%           file normally. 
%03/09/2016: Fixed a bug whereby if numb_frags>2, the atomic labels for
%            frag 2 were used to frag 2-->end. (current_column = final_col+1; changed
%            to curr_column = final_col+1;)
%04/09/2016: V3 finally (probably) fixed it so all options work.
%(merge_spd=0,merge_fragments=1 wasnt previously working. Still needs error checking tho.
%27/08/2017: Made it work on windows AND mac/linux by using fullfile function to join
%	     folder/file strings

%OUTPUT:
%header = a comma delim string with the names of each data series
%output_data = nxm matrix of data(duh). Find the label for each column in
%              header

dbstop if error
%WARNING:STILL NEEDS TESTING WITH >2 fragments being used
%BUG: If merge_spd=0 and merge_fragments=0, it doesnt merge the fragments
%IMPORTANT GLOBAL VARIABLES GO HERE.
merge_spd = 0;  %set to 1 to sum all the N2s and N2p (for example), so you just get ATOMIC contributions,not AO contributions
merge_fragments = 1;  %set to 1 to merge all the fragments (e.g cation/anion). So you'd get Nitrogen contributions to spectrum, but not cationic nitrogen/anionic nitrogen)
save_txt = 1;   %set to 1 to save .txt file with final labels/data in it
save_fig = 1;   %set to 1 to automatically save figure in Folder
lwidth=2;
font_size=10;
font_type='arial';








%Imports all important imformation from the relevant output file
%frag_name_list = cell array (1xn) with names of main fragments (e.g
%                   cation, anion)
%sub_frag_numb_list = cell array (1xn) with numbers of sub fragments (e.g N2s) in
%                     each fragment(e.g cation). eg = {7,3,4}
%atom_label_list = cell array (1xn) containing ALL sub-fragments (e.g N1s)
%energies = mx1 column matrix with energies for the plot we'll make
%           eventually
%intensities = mxn matrix, each column is one sub-fragments contribution to
%              the calculated spectrum
    function [frag_name_list,sub_frag_numb_list,atom_label_list,energies,intensities] = Extract_Data(folder); 
        
        %Gets the full path of the relevant file
        dir_struct = dir(folder)
        for i=1:size(dir_struct)
            if strfind(dir_struct(i).name,'ALL_Contribs')~=-1
		fname = fullfile(folder,dir_struct(i).name)
%                fname = strcat(folder,'\',dir_struct(i).name);
                break;
            end
        end
        
        %Extracting the useful info from the file
        raw_data = dlmread(fname,',',3);
        energies = raw_data(:,1);
        intensities = raw_data(:,2:end);
        clear raw_data;
        
        %Extracts the file header
        fobject1=fopen(fname);
        header = fileread(fname);
        header = strsplit(header,'\n')';
        header = header(1:3);	%top 3 lines of input file as 3x1 cell array
        fclose(fobject1);
        
        %Organising header data in a better way
        header=strrep(header,'#','');
        header{2,1} = strrep(header{2,1},'Sub fragments in each=,','');
        header{3,1} = strrep(header{3,1},'Energy(eV),','');
        frag_name_list = strsplit(header{1},',');
        sub_frag_numb_list = strsplit(header{2},',');
        atom_label_list = strsplit(header{3},',');
    end



%Generates intensities and labels for atomic contribs of single fragment
%(possibly with s/p/d intensities summed togerher
%intensities = matrix of intensities for THIS FRAGMENT ONLY (NOT the full
%               matrix)
%sub_fragment_list = list of N2S/S3s labels FOR THIS FRAGMENT ONLY
%fragment name = name of the fragment, to be put in front of all the labels
%merge_spd_opt = 1 if you want the s/p/d of individual elements summed.
%fragment_intensities = list of intensities for the gaussian peaks from this
%                   fragment. Same as input intensities if merge_spd_opt=0
%all_frag_labels = 1xn cell array with labels of each fragment (1 per col in
%               intensities)
    function [fragment_intensities,all_frag_labels] = Gen_Fragment_Intensities(intensities,sub_fragment_list,fragment_name,merge_spd_opt)

        
        if merge_spd_opt==1;
            %Removing s/p/d labels off the sub-fragments
            for i=1:size(sub_fragment_list,2)
                numb_start = regexp(sub_fragment_list{i},'[0-9]');
                tstring = sub_fragment_list{i}
                tstring = tstring(1:numb_start-1)
                sub_fragment_list{i} = tstring
            end
        end
            %Grouping all atomic contribs together
            all_frag_labels{1} = sub_fragment_list{1} %holds list of each type of element
            counter=2;		%tracks which cell of all_frag_labels we're on
            for i=1:size(sub_fragment_list,2)
                if ~any(strcmp(sub_fragment_list{i},all_frag_labels))
                    all_frag_labels{counter} = sub_fragment_list{i};
                    counter=counter+1;
                end
            end
            
            
            
            clear counter;
            %Now summing together the relevant columns + creating new intensities matrix
            fragment_intensities = zeros(size(intensities,1),size(all_frag_labels,2));
            counter=1;
            for i=1:size(all_frag_labels,2);
                temp_string = all_frag_labels{i};
                for j=1:size(sub_fragment_list,2);
                    if strcmp(temp_string,sub_fragment_list{j});
                        fragment_intensities(:,i) = fragment_intensities(:,i) + intensities (:,j);
                    end
                end
            end
          
        
%         if merge_spd_opt~=1
%            fragment_intensities = intensities;
%            all_frag_labels = sub_fragment_list;
%         end
        
        %finally comibining the fragment labels and atomic labels
        for i=1:size(all_frag_labels,2)
            all_frag_labels{i} = strcat(fragment_name,'(',all_frag_labels{i},')');
        end
    end


 



    %Just creates a plot. Plots summed intensities as part of it
    function [header output_data] = Create_Plot(energies,intensities,label_list,txt_save,fig_save,save_folder)
        
        %Assing summed_intensities into the mix
        summed_intensities = zeros(size(intensities,1),1);
        for i=1:size(summed_intensities,1);
           summed_intensities(i,1) = sum(intensities(i,:)); 
        end
        intensities = [summed_intensities intensities];
        label_list = ['Summed Intensities' label_list];
        
        xlim = [min(energies(:,1)) max(energies(:,1))];
        ylim = [0 max(intensities(:,1))+(0.1*max(intensities(:,1)))];
        
        
        figure1=figure;
        ax1=axes('parent',figure1);
        hold all;
        for i=1:size(intensities,2)
            plot(energies(:,1),intensities(:,i),'linewidth',lwidth);
        end
        
        leg1 = legend(label_list);
        
        
        set(ax1,'ylim',ylim);
        set(ax1,'xlim',xlim);
        set(get(ax1,'xlabel'),'string','Binding Energy / eV','fontsize',font_size,'fontname',font_type);
        set(get(ax1,'ylabel'),'string','Intensity / arb.units','fontsize',font_size,'fontname',font_type);
        set(ax1,'xdir','reverse');
        set(ax1,'fontsize',font_size,'fontname',font_type);
        set(leg1,'fontsize',font_size,'fontname',font_type);
        
        
        %File saving:
        txt_name =strcat(save_folder,'\Plotted_Atomic_Contribs_spd',num2str(merge_spd),'_frag',num2str(merge_fragments),'.txt');
        fig_name =strcat(save_folder,'\Plotted_Atomic_Contribs_spd',num2str(merge_spd),'_frag',num2str(merge_fragments),'.fig');
        
        
        if fig_save==1;
            savefig(fig_name);
        end
        
        if txt_save==1;
           %writes header
           fobject1=fopen(txt_name,'w');
           header = strcat('Energies / eV,',strjoin(label_list,','));
           fprintf(fobject1,'%s\n',header);
           fclose(fobject1);
           %writes the rest
           output_data = [energies intensities];
           dlmwrite(txt_name,output_data,'-append');
        end

        
        
        
        
        end



    function [] = Main()
        [frag_name_list,sub_frag_numb_list,atom_label_list,energies,intensities_init] = Extract_Data(folder); 
        
        
        if merge_fragments~=1
        end_col = str2num(sub_frag_numb_list{1});
        [all_fragment_intensities,all_frag_labels] = Gen_Fragment_Intensities(intensities_init(:,1:end_col),atom_label_list(1:end_col),frag_name_list{1},merge_spd);
         
        curr_column = end_col+1;   %tracks the current column we're on in intensities matrix
        for fragment=2:size(frag_name_list,2);
             final_col = curr_column + str2num(sub_frag_numb_list{fragment})-1;  %column this fragments intensities finish on
             [curr_fragment_intensities,curr_frag_labels] = Gen_Fragment_Intensities(intensities_init(:,curr_column:final_col),atom_label_list(curr_column:final_col),frag_name_list{fragment},merge_spd);
             all_fragment_intensities = [all_fragment_intensities curr_fragment_intensities];
             all_frag_labels =[all_frag_labels curr_frag_labels];
             curr_column = final_col+1;    %where next fragment starts.
         end
        end
        
        
        %NOT really TESTED BECAUSE I CBA.(probably works tho)
        if merge_fragments==1;
        end_col = str2num(sub_frag_numb_list{1});
        [all_fragment_intensities,all_frag_labels] = Gen_Fragment_Intensities(intensities_init,atom_label_list,frag_name_list{1},merge_spd);
        all_frag_labels = strrep(all_frag_labels,frag_name_list{1},''); %Just leaving elemental symbol for the legend
        all_frag_labels = strrep(strrep(all_frag_labels,'(',''),')','');
        end

        [header output_data] = Create_Plot(energies,all_fragment_intensities,all_frag_labels,save_txt,save_fig,folder);

        
end


Main;


end