% Matlab version of the Python script "NormalEquation.py"
% Create a set of vectors G representing the spectra of the 3 components
NumPoints=1000;
X=linspace(200,800,NumPoints)'; % Creates NumPoints x 1 vector X
center1=450;
G1=gaussian(X,center1,1,100); % Gaussian function defined in lines 38-43
center2=500;
G2=gaussian(X,center2,1,110);
center3=550;
G3=gaussian(X,center3,1,120);

E=[G1,G2,G3]; % compose the separate G vectors into a matrix E
noiseamp=.01;
% Specify the amplitudes of the components in the mixture
amp1=1; % Amplitude (height) of peak 1
amp2=.5; % Amplitude (height) of peak 2
amp3=.3; % Amplitude (height) of peak 2
% Compute the simulated noisy mixture spectrum
A=amp1.*G1+amp2.*G2+amp3.*G3+noiseamp.*randn(length(X),1) ;

tic % Start clock for execution timimg
% Apply normal equation to compute concentration vector C
C = inv(E'*E)*E'*A;

plot(X,A,X,C(1).*G1,X,C(2).*G2,X,C(3).*G3)
PercentDifference1=100.*(amp1-C(1))/amp1;
PercentDifference2=100*(amp2-C(2))/amp2;
PercentDifference3=100*(amp3-C(3))/amp2;
text(center1,C(1),num2str(PercentDifference1),'FontSize',14);
text(center2,C(2),num2str(PercentDifference2),'FontSize',14);
text(center3,C(3),num2str(PercentDifference3),'FontSize',14);
title('Matlab: Percent errors of peak heights measured by Classical Least Squares')
toc % stop clock and print elapsed time

disp(C) % Print measured amplitudes

% Define a Gaussian peak function
function g=gaussian(x,center,amplitude,width)
% This generates a single peak with a gaussian lineshape,
% with a specific center, amplitude, and width.
% width is the FWHM (full width at half maximum)
g=amplitude*exp(-((x-center)/(0.60056120439323*width)).^2);
end