clear
close all

% Theoretical model of atomic resonance
x=-15:.1:15;
Lw=1.4;%  true value of linewidth
x0=3;% true value of central frequency
BG=.2% true value of background
Y=@(z)1/pi*(Lw./(Lw^2+(z-x0).^2))+BG;% Lorentzian with true 
plot(x,Y(x),'linewidth',2)

xlabel('Frequency (MHz)','fontsize',16)
ylabel('Intensity (Fluorescence)','fontsize',16)
set(gca,'fontsize',16)
title('Lorentzian distribution evaluated with true values')
%%Generate data for fitting using random Gaussian noise
N=10%% Number of experimental runs
xdata=-15:1.6:15;% Points in the experiment that the data was taken
Sigy=.04;% Noise level of our experiment
for ii=1:1:N%Generate the data for 10 independent runs
Ynoise=Sigy*randn(size(xdata));
Ymeas=Y(xdata)+Ynoise;
YY(:,ii)=Ymeas;%This is the data for all our "N runs"
end
figure

Yave=mean(YY,2);%% Average of the "N" experimental runs

%Deviations
for jj=1:1:N
    Dev(:,jj)=YY(:,jj)-Yave;
end
  %Population standard deviation
SigYave= sqrt(1/(N-1))*sqrt(sum((Dev).^2,2));

plot(x,Y(x))
hold on
errorbar(xdata,Yave,SigYave,'b.')

xlabel('Frequency (MHz)','fontsize',16)
ylabel('Intensity (Fluorescence)','fontsize',16)
set(gca,'fontsize',16)
title('Generated data with random noise. This is our data')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% FIT A LORENTZIAN TO OUR DATA
%%%FIT A FUNCTION
%% Find A and B by minimizing a function varying fit parameters "a"

% Function for our model in range "xx" and parameters a(1),a(2),a(3)
% We are assuming that our model is a Lorentzian. If you 
% want to fit another function, you need to change this model
Ymodel=@(xx,a)1/pi*(a(1)./(a(1)^2+(xx-a(2)).^2))+a(3);

% This is the function defining chi square.
Chi2 = @(a) sum((Yave-Ymodel(xdata,a)').^2./SigYave.^2);

%Initial guess parameters 
GuessParam=[2 1 .2];

% "fminsearch" is the minimization function. Minimizes Chi2 taking the
% initial guess values of the parameters a(1)=Lw and a(2)=x0 and a(3)=BG
[afit,fval] = fminsearch(Chi2,GuessParam)

% "afit" are the results of the fit

Yth=@(xx,a)1/pi*(afit(1)./(afit(1)^2+(xx-afit(2)).^2))+afit(3);
Ygeneral=@(xx,a)1/pi*(a(1)./(a(1)^2+(xx-a(2)).^2))+afit(3);
figure

plot(x,Ygeneral(x,GuessParam),'linewidth',2)
hold on
errorbar(xdata,Yave,SigYave,'r.','linewidth',2)
xlabel('Frequency (MHz)','fontsize',16)
ylabel('Intensity (Fluorescence)','fontsize',16)
set(gca,'fontsize',16)
legend('Guess parameters','Noisy data')

figure
plot(x,Y(x),'linewidth',2)
hold on
errorbar(xdata,Yave,SigYave,'r.','linewidth',2)
plot(x,Yth(x),'k','linewidth',2)

legend('True parent dist.','Noisy data','Fit results')
xlabel('Frequency (MHz)','fontsize',16)
ylabel('Intensity ((Fluorescence)','fontsize',16)
set(gca,'fontsize',16)
%%%%%%%%%%%5 HERE WE WILL FIND THE ERRORS IN THE PARAMETER ESTIMATION
%Find error parameter afit(1)
Xa1=(afit(1)-afit(1)/2):afit(1)/50:(afit(1)+afit(1)/2);

for kk=1:1:length(Xa1)
    
    Ya1chi2(kk)=Chi2([Xa1(kk) afit(2) afit(3)]);
end

figure
plot(Xa1,Ya1chi2)
c=min(Ya1chi2);
cpar=Xa1(Ya1chi2==min(Ya1chi2));

title(['Chi square varying parameter a1, Linewidth. Chi^2_{min}=',num2str(c),' at a1=',num2str(cpar)])
xlabel('Change parameter Fit param 1','fontsize',16)
ylabel('Chi square','fontsize',16)
set(gca,'fontsize',16)

%Find zeroes
f=@(x) Chi2([x afit(2) afit(3)])-(min(Ya1chi2)+1); 
z1(1) = fzero(f,.5*afit(1));
z1(2) = fzero(f,1.5*afit(1));
%Standard deviation for parameter 1. Note that may not be symmetric
afit(1)
Sig1=abs((afit(1)-z1))

%Find error parameter afit(2)
Xa2=(afit(2)-afit(2)/2):afit(2)/50:(afit(2)+afit(2)/2);

for kk=1:1:length(Xa2)
    
    Ya2chi2(kk)=Chi2([afit(1) Xa2(kk) afit(3)]);
end

figure
plot(Xa2,Ya2chi2)
c=min(Ya2chi2);
cpar=Xa2(Ya2chi2==min(Ya2chi2));

title(['Chi square varying parameter a2, Central frequency. Chi^2_{min}=',num2str(c),' at a2=',num2str(cpar)])
xlabel('Change parameter Fit param 2','fontsize',16)
ylabel('Chi square','fontsize',16)
set(gca,'fontsize',16)
%Find zeroes
f=@(x) Chi2([afit(1) x afit(3)])-(min(Ya2chi2)+1); 
z2(1) = fzero(f,.5*afit(2));
z2(2) = fzero(f,1.5*afit(2));
%Standard deviation for parameter 1. Note that may not be symmetric
afit(2)
Sig2=abs((afit(2)-z2))

%Find error parameter afit(3)
Xa3=(afit(3)-afit(3)/2):afit(3)/50:(afit(3)+afit(3)/2);

for kk=1:1:length(Xa3)
    
    Ya3chi2(kk)=Chi2([afit(1) afit(2) Xa3(kk)]);
end

figure
plot(Xa3,Ya3chi2)
c=min(Ya3chi2);
cpar=Xa3(Ya3chi2==min(Ya3chi2));
title(['Chi square varying parameter a3, background. Chi^2_{min}=',num2str(c),' at a3=',num2str(cpar)])
xlabel('Change parameter Fit param 3','fontsize',16)
ylabel('Chi square','fontsize',16)
set(gca,'fontsize',16)

%Find zeroes
f=@(x) Chi2([afit(1) afit(2) x])-(min(Ya3chi2)+1); 
z3(1) = fzero(f,.8*afit(3));
z3(2) = fzero(f,1.5*afit(3));
%Standard deviation for parameter 3. Note that may not be symmetric
afit(3)
Sig3=abs((afit(3)-z3))

%Calculate Reduced chi square
Chi2(afit)
Chi2(afit)/(length(xdata)-3)
    
    

