% Mass assignment parameters
m1 = 48; % mass of lighter ion used for mass assignment
t1 = 9721; % bin number of center of peak m1
m2 = 256;% mass of heavier ion used for mass assignment
t2 = 22690; % bin number of center of peak m2
prominence = 40; %HOW sensitive

% turn different procedures on and off
AutoMass = 0; % Assign masses automatically by varying m1 and m2 for the given peaks t1 and t2
Generate_EPS = 1; % set to 1 to generate EPS file containing TPD profiles
VUV_Correction = 0; % set to 1 to generate VUV normalization constant for each point in TPD profile
Waterfall_Bin = 1; % Binning along mass axis for waterfall plots, set to 0 to turn off

% TPD profile parameters input
MinimumMass = 1; % Minimum m/z for which TPD profile/waterfall will be generated
MaximumMass = 400; % Maximum m/z for which TPD profile/waterfall will be generated
TempStart = 0; % Temperature at which TPD started
TempStep = 1; % TPD rate (K/min)
VUV_TStart = 19; %1 Temperature Tek Acquire was turned on
VUV_TEnd = 341; % Temperature Tek Acquire was turned off, usually 320 K
MassResolutionTPD = 1; % Width of integration channel for TPD profiles
MassResolutionBin = 0.1; % Width of integration channel for waterfall and contour plots
MassSearch = 3; % Search for correct mass assignment for given m1 and m2 by varying each by +/- MassSearch when AutoMass = 1
range = 30016; % From header of .mp file ir from older software will be in .888 or .889
Cycles = 325; % From header of .mp file ir from older software will be in .888 or .889

% Request input .mp and import relevant data
[FileName,PathName] = uigetfile('*.mp','Multiselect','on'); % Change *.mp to *.asc to analyze .asc data files
fileposition = strcat(PathName,FileName);
data_import = importdata(fileposition,' ',54);
RawData = reshape(data_import.data,range,Cycles);
MSData = sum(RawData,2);
ToF = (1:range); % All ToF data is analyzed in units of bin number. Conversion to units of seconds is unnescessary
clear FileName PathName fileposition

% Mass assignment
A = (m1 + m2 - 2 * sqrt(m1 * m2))/(t1 - t2)^2;
B = -t2 + sqrt(m2/(m1 + m2 - 2 * sqrt(m1 * m2))) * abs(t1 - t2);
Mass = A*(ToF+B).^2;
clear A B

% Automated mass assignment, not for use with blank/low intensity data
SmoothData = smoothdata(MSData,'gaussian',10); % smooth data so peak centers are more prominent
if AutoMass == 1 
    [~,Peaks] = findpeaks(SmoothData,'MinPeakDistance',35,'MinPeakProminence',prominence); % identify peaks, reduce MinPeakProminence for low intensity data or just use a prior mass assignment and do not use the automatic procedure 
    t1 = interp1(Peaks,Peaks,t1,"nearest"); % Round t1 to the nearest peak identified in prior line
    t2 = interp1(Peaks,Peaks,t2,"nearest"); % same as above for t2
    MassTest = zeros((MassSearch*2+1)^2,5);
    for i = 1:MassSearch*2+1 % This loop enters test values of m1 and m2 into column 1 and 2 of MassTest, respectively.
        for j = 1:MassSearch*2+1
            MassTest(i+(j-1)*(MassSearch*2+1),1) = m1-MassSearch+i-1;
            MassTest(i+(j-1)*(MassSearch*2+1),2) = m2-MassSearch+j-1;
        end
    end
    MassTest(MassTest(:,1)<1,:) = []; % these three lines discard unphysical values of m1 and m2, i.e. negative values of either mass and all cases of m1 > or equal to m2
    MassTest(MassTest(:,2)<1,:) = [];
    MassTest(MassTest(:,1)>=MassTest(:,2),:) = [];
    for i = 1:length(MassTest) % This loop calculates value of constants A and B for mass assignment using the test values enumerated in the previous loop and stores them in columns 3 and 4 of MassTest.
        m1 = MassTest(i,1);
        m2 = MassTest(i,2);
        MassTest(i,3) = (m1 + m2 - 2 * sqrt(m1 * m2))/(t1 - t2)^2;
        MassTest(i,4) = -t2 + sqrt(m2/(m1 + m2 - 2 * sqrt(m1 * m2))) * abs(t1 - t2);
    end
    for i = 1:length(MassTest) % This loop calculates the mean absolute error for each peak position to the nearest integer given a mass assignment, and lists the error in the 5th column of MassTest.
        A = MassTest(i,3);
        B = MassTest(i,4);
        PeakMass = A*(Peaks+B).^2;
        MassTest(i,5) = mean((PeakMass - round(PeakMass)).^2);
    end
    [~,idx] = min(MassTest(:,5)); % Finds the row of MassTest containing the most accurate mass assignment and uses these to assign all masses
    m1 = MassTest(idx,1);
    m2 = MassTest(idx,2);
    A = MassTest(idx,3);
    B = MassTest(idx,4);
    Mass = A*(ToF+B).^2; % Applies the optimized mass assignment to the data.
    clear AutoMass A B idx PeakMass MassTest i j Peaks
end

% Plot mass spectra with x-axes in mass and again with smoothed data in bins
figure(1)
Export_MS = [Mass' MSData];
plot(Mass,MSData);
figure(2)
plot(ToF,SmoothData);
findpeaks(SmoothData,'MinPeakDistance',100,'MinPeakProminence',prominence);

% Temperture scale calculation
temp = (TempStart+TempStep:TempStep*2:Cycles*TempStep*2+TempStart-TempStep);
clear TempStart TempStep

% VUV power correction calculations
if VUV_Correction == 1
    [FileName2,PathName2] = uigetfile('*.txt');
    VUV = importdata(strcat(PathName2,FileName2));
    VUV = VUV(:,2);
    T_VUV = linspace(VUV_TStart,VUV_TEnd,length(VUV));
    T_VUV = interp1(temp,temp,T_VUV,'nearest');
    Export_VUV = zeros(length(temp),1);
    for i = (1:length(temp))
      TIndex = T_VUV == temp(i);
      VUV_Intensity = mean(VUV(TIndex,:));
      Export_VUV(i) = VUV_Intensity;
    end
    Export_VUV = Export_VUV/1000;
    Export_VUV(isnan(Export_VUV)) = 1;
    clear VUV_TEnd VUV_TStart VUV_Correction TIndex T_VUV VUV FileName2 PathName2 VUV_Intensity
    else
    Export_VUV = ones(length(temp),1);
    clear VUV_Correction
end

% TPD profile calculation and .eps
if Generate_EPS == 1
    Export_TPD = zeros(Cycles,MaximumMass-MinimumMass+1);
    j = MinimumMass:MassResolutionTPD:MaximumMass;
    for i = MinimumMass:MaximumMass
        RoundedMass = interp1(j,j,Mass,"nearest");
        MassIndex = RoundedMass == i;
        TPD_Profile = sum(RawData(MassIndex,:),1);
        figure(3)
        plot(temp,TPD_Profile(1:length(temp)));
        Export_TPD(:,i) = TPD_Profile()';
        xlabel('Temperature (K)');
		    ylabel('Counts');
		    title(strcat({num2str(i)},{' m/z'}, {' '}),'Interpreter','none');
            YMax = max(TPD_Profile); % Calculates the maximum y-value in each TPD for scaling purposes, can be modified to TPD_Profile(n,end) to ignore regions of detector saturation.
		    if max(TPD_Profile) > 10 
			    YMax = YMax + 0.1 * YMax; %set ymax of plot
		    else
			    YMax = 10; % if nothing is there, set a ymax anyway to not be confused with noise
            end
		    axis([1 325 0 YMax]);
		    set(gcf, 'PaperUnits', 'inches');
		    set(gcf, 'PaperType', 'usletter');
		    set(gcf, 'PaperPosition', [0.5 0.5 8.5 4.5]);
		    set(gca,'FontSize',20,'LineWidth',1.25)
	    if i == MinimumMass
		    print ('-dpsc', strcat('TPD Profiles.eps'));
	    else
		    print ('-dpsc', strcat('TPD Profiles.eps'), '-append');
        end
    end
    Export_TPD = [temp',Export_TPD];
    clear Generate_EPS i MassIndex YMax
end

% Waterfall plot data processing. This bins data along the mass axis into bins of width MassResolutionBin so Origin can handle these large datasets more easily
if Waterfall_Bin == 1
    Export_Waterfall = [];
    j = MinimumMass:MassResolutionBin:MaximumMass;
    for i = j
        RoundedMass = interp1(j,j,Mass,"nearest");
        MassIndex = RoundedMass == i;
        TPD_Profile = sum(RawData(MassIndex,:),1);
        Export_Waterfall = [Export_Waterfall,TPD_Profile()'];
    end
    Export_Waterfall = [[MinimumMass:MassResolutionBin:MaximumMass]',Export_Waterfall'];
    clear TPD_Profile MassIndex RoundedMass i j Waterfall_Bin
end

clear Mass temp ToF Waterfall_Bin VUV_Correction Generate_EPS AutoMass MinimumMass MaximumMass TempStart TempStep VUV_TStart VUV_TEnd MassResolutionTPD MassResolutionBin MassSearch range Cycles sweeps


