% 
% Calculates isotopic distributions including isotopic fine structure 
% of molecules using FFT and various scaling 'tricks'. Easily adopted
% to molecules of any elemental composition (by altering MAX_ELEMENTS 
% and the nuclide matrix A). To simulate spectra, convolute with peak
% shape using FFT.
%
% (C) 1999 by Magnus Palmblad, Division of Ion Physics, Uppsala Univ.
% Acknowledgements:
% Lars Larsson-Cohn, Dept. of Mathematical Statistics, Uppsala Univ.,
% for help on theory of convolutions and FFT. 
% Jan Axelsson, Div. of Ion Physics, Uppsala Univ. for comments and ideas
%
% Contact Magnus Palmblad at magnus.palmblad@gmail.com if you should
% have any questions or comments.
%

clear; clf; format long;

MAX_ELEMENTS=5+1;  % add 1 due to mass correction 'element'
MAX_ISOTOPES=4;    % maxiumum # of isotopes for one element
CUTOFF=1e-8;       % relative intensity cutoff for plotting

WINDOW_SIZE=input('Window size (in Da) ---> ');

RESOLUTION=input('Resolution (in Da) ----> ');  % mass unit used in vectors
if RESOLUTION < 0.00001  % minimal mass step allowed 
  RESOLUTION = 0.00001;
elseif RESOLUTION > 0.5  % maximal mass step allowed
  RESOLUTION = 0.5;
end

R=0.00001/RESOLUTION;  % R is used to scale nuclide masses (see below)

WINDOW_SIZE=WINDOW_SIZE/RESOLUTION;   % convert window size to new mass units
WINDOW_SIZE=2^nextpow2(WINDOW_SIZE);  % fast radix-2 fast-Fourier transform algorithm

if WINDOW_SIZE < round(496708*R)+1
  WINDOW_SIZE = 2^nextpow2(round(496708*R)+1);  % just to make sure window is big enough
end

disp(sprintf('Vector size: 1x%d',WINDOW_SIZE));

M=[378 254 65 75 6 0];  % empiric formula, e.g. bovine insulin 

% isotopic abundances stored in matrix A (one row for each element)
A=zeros(MAX_ELEMENTS,MAX_ISOTOPES,2);

A(1,1,:)=[100783 0.9998443];                 % 1H     
A(1,2,:)=[201410 0.0001557];                 % 2H
A(2,1,:)=[100000 0.98889];                   % 12C
A(2,2,:)=[200336 0.01111];                   % 13C
A(3,1,:)=[100307 0.99634];                   % 14N
A(3,2,:)=[200011 0.00366];                   % 15N
A(4,1,:)=[099492 0.997628];                  % 16O
A(4,2,:)=[199913 0.000372];                  % 17O
A(4,3,:)=[299916 0.002000];                  % 18O
A(5,1,:)=[097207 0.95018];                   % 32S
A(5,2,:)=[197146 0.00750];                   % 33S
A(5,3,:)=[296787 0.04215];                   % 34S
A(5,4,:)=[496708 0.00017];                   % 36S
A(6,1,:)=[100000 1.00000];                   % for shifting mass so that Mmi is
                                             % near left limit of window

Mmi=[round(100783*R) round(100000*R) round(100307*R) round(099492*R) ...
      round(097207*R) 0]*M';  % (Virtual) monoisotopic mass in new units

% mass shift so Mmi is in left limit of window:

FOLDED=floor(Mmi/(WINDOW_SIZE-1))+1;  % folded FOLDED times (always one folding due to shift below)

% shift distribution to 1 Da from lower window limit:
M(MAX_ELEMENTS)=ceil(((WINDOW_SIZE-1)-mod(Mmi,WINDOW_SIZE-1)+round(100000*R))*RESOLUTION);  

MASS_REMOVED=[0 11 13 15 31 -1]*M';  % correction for 'virtual' elements and mass shift

ptA=ones(WINDOW_SIZE,1);
t_fft=0; t_mult=0;

for i=1:MAX_ELEMENTS
  tA=zeros(WINDOW_SIZE,1);
  for j=1:MAX_ISOTOPES
    if A(i,j,1) ~= 0
      tA(round(A(i,j,1)*R)+1)=A(i,j,2);  % put isotopic distribution in tA
    end
  end
  disp('Calculate FFT...')
  t0=clock;
  tA=fft(tA) ;  % FFT along elements isotopic distribution  O(nlogn)
  t_fft=t_fft+etime(clock,t0);
  disp('Multiply vectors...');
  t0=clock;
  tA=tA.^M(i);  % O(n)
  ptA=ptA.*tA;  % O(n)
  t_mult=t_mult+etime(clock,t0);
end

clear tA;

disp(sprintf('Time for FFT: %4.2f s',t_fft))
disp(sprintf('Time for multiplications: %4.2f s',t_mult))

disp('Calculate IFFT...')
t0=clock;
ptA=real(ifft(ptA));  % O(nlogn)

disp(sprintf('Time for IFFT: %4.2f s',etime(clock,t0)))

disp('Plotting...')
t0=clock;
clf;
MA=linspace((FOLDED*(WINDOW_SIZE-1)+1)*RESOLUTION+MASS_REMOVED,(FOLDED+1)*(WINDOW_SIZE-1)*RESOLUTION+MASS_REMOVED, WINDOW_SIZE-1)';
axis([MA(1) MA(WINDOW_SIZE-1) 0 1]);

ind=find(ptA>CUTOFF)';
P=[MA(ind) ptA(ind)];

H=bar(P(:,1),P(:,2),0);
set(H,'FaceColor','k');
grid
zoom on;

disp(sprintf('Time for plotting: %4.2f s',etime(clock,t0)))