clear all
close all

% Mass assignment parameters
m1 = 1; 
bin1 = 1235; 
m2 = 83;
bin2 = 12861;

%On/Off switches
Waterfall_Bin = 0; % Binning along mass axis for waterfall plots, set to 0 to turn off
VUV_Correction = 1; % set to 1 to generate VUV normalization constant for each point in TPD profile
Generate_EPS = 1; % set to 1 to generate EPS file containing TPD profiles

% TPD profile parameters input
MinimumMass = 1; % Minimum m/z for which TPD profile will be generated
MaximumMass = 200; % Maximum m/z for which TPD profile will be generated
TempStart = 5; % Temperature at which TPD started
TempStep = 1; % TPD rate (K/min)
smoothing = 1; % # of points used for Boxcar averaging of TPD profiles
VUV_TStart = 160; % Temperature Tek Acquire was turned on
VUV_TEnd = 361; % Temperature Tek Acquire was turned off, usually 320 K

% Request input *.mp and import relevant data
[FileName,PathName] = uigetfile('*.mp','Multiselect','on');
fileposition = strcat(PathName,FileName);
data_import = importdata(fileposition,' ',54);
range = str2double(data_import.textdata{2}(7:end));
cycles = str2double(data_import.textdata{12}(8:end));
BinWidth = 0.8*2^str2double(data_import.textdata{43}(10));
sweeps = str2num(data_import.textdata{7}(10:end));
RawData = reshape(data_import.data,range,cycles);
Data = sum(RawData,2);
ToF = (1:range)*BinWidth;

% Mass assignment
t1 = bin1 * BinWidth;
t2 = bin2 * BinWidth;
A = (m1+m2)/(t1-t2)^2-2*sqrt(m1*m2)/abs(t1-t2)^2;
B = sqrt(m2/A)-t2;
Mass = A*(ToF+B).^2;
% t3=14961*BinWidth;
% mass_t3=A*(17581*BinWidth+B).^2


% Temperture scale calculation
temp = (TempStart+TempStep*sweeps/1800/2:TempStep*sweeps/1800:cycles);

% 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);
Export_VUV_e = zeros(length(temp),1);%% added error
for i = (1:length(temp))
  TIndex = T_VUV == temp(i);
  VUV_Intensity = mean(VUV(TIndex,:));
  Export_VUV(i) = VUV_Intensity;
  VUV_Intensity_e = std(VUV(TIndex,:));
  Export_VUV_e(i) = VUV_Intensity_e;
end
Export_VUV = Export_VUV/1000;
Export_VUV_e = Export_VUV_e/1000;
Export_VUV(isnan(Export_VUV)) = 1;
Export_VUV_e(isnan(Export_VUV)) = 1;
else
Export_VUV = ones(length(temp),1);
Export_VUV_e = ones(length(temp),1);
end

% TPD profile calculation and .eps
if Generate_EPS == 0
Export_TPD = [];
for i = MinimumMass:MaximumMass
    MassIndex = round(Mass) == i;
    TPD_Profile = smooth(sum(RawData(MassIndex,:),1),smoothing);
    plot(temp,TPD_Profile(1:length(temp)));
    Export_TPD = [Export_TPD,TPD_Profile()];
    xlabel('Temperature (K)');
		ylabel('Counts');
		title(strcat({num2str(i)},{' m/z'}, {' '}),'Interpreter','none');
        YMax = max(TPD_Profile);
		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
end

Export_MS = [Mass' Data];
plot(Mass,Data);

% Waterfall plot data processing
if Waterfall_Bin ~= 0
idx = (Mass > 20 & Mass < 200);
idx2 = (diff(Mass)>0);
j=1;
for i=1:length(Mass)-1
if idx(i)==1 && idx2(i)==1
mass3(j,1)=Mass(i);
counts2(j,:)=RawData(i,:);
j=j+1;
end
end
for i=1:floor(length(mass3)/Waterfall_Bin)
    mass4(i)=mean(mass3((i-1)*Waterfall_Bin+1:i*Waterfall_Bin+1));
    counts3(i,:)=sum(counts2((i-1)*Waterfall_Bin+1:i*Waterfall_Bin+1,:));
end
Export_Waterfall = [mass4',counts3];
end


figure (111)
plot(1:range,Data);