%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eksperimentel Fysik, Niels Bohr Institutet, Københavns Universitet
% Tirsdag d. 15. maj og torsdag d. 17. maj (uge 4), 2012
%
% Dette dokument indeholder løsninger til opgaven med et spektrogram. Disse 
% opgaver blev lavet af Kim Lefmann og lagt op som et PDF dokument på Absalon.
% 
% Pia Jensen, bozack@bozack.dk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all     %tømmer alle gemte variable
close all     %luk alle eksterne vinduer (figurer o.lign.)
clc           %rens kommando-vinduet


%% Opgave 1 og 2 ----------------------------------------------------------------
%  I opgave 1 skal der bare fittes en enkelt gang, mens der i opgave 2 skal gøres
%  det samme mange gange i træk, og bagefter skal der ses på statistikken

% Definerer parametre for den Gaussiske model
a       = 100;                      %signal-amplitude
b       = 10;                       %baggrunds-amplitude
x0      = 5000;                     %centrum for linie
sig     = 50;                       %bredde (sigma) for linie
t_meas  = 100;                      %samlet tælletid
N       = 100;                      %antal ønskede måleserier

% Initialiserer arrays til fit-parametre
a_fit   = zeros(1,N);               %fittet signal-amplitude
b_fit   = zeros(1,N);               %fittet baggrunds-amplitude
x0_fit  = zeros(1,N);               %fittet centrum for linie
sig_fit = zeros(1,N);               %fittet bredde (sigma) for linie

% Plan over eksperiment
x_min   = 4500;                     %minimum bølgelængde
x_max   = 5500;                     %maksimum bølgelængde
x_step  = 10;                       %bølgelængde step-størrelse
x_exp   = x_min : x_step : x_max;   %bølgelængder der skal måles ved
t       = t_meas / length(x_exp);   %tid der tælles for hver bølgelængde

% Lav model-datasættet
y_model = t * a * exp( - (( x_exp - x0 ).^2) ./ (2*sig^2) ) + t * b;

% Loop over eksperimenter
for j = 1:N   % Loop over forskelige fits, opgave 2 del
    
    % Tilføj støj til måleserien
    err_model = sqrt(y_model);  %modellens usikkerhed
    y_exp     = abs(round( y_model + err_model .* randn(size(err_model)) ));
    err_exp   = sqrt(y_exp);    %usikkerheden som man kender i eksperimentet
    
    % Plot datasættet
    figure(1)
    errorbar(x_exp,y_exp,err_exp,'b.','MarkerFaceColor', [0.4 0.4 0.8],...
          'MarkerEdgeColor',[0.4 0.4 0.8])
    axis([x_min-50 x_max+50 0 (a+b)*1.2*t])
    xlabel('Bølgelængde [Å]')
    ylabel('Intensitet [tællinger/s]')
    hold on

    % Lav fittet
    f1 = fit(x_exp',y_exp','a*exp(-((x-b).^2)/(2*c^2))+d',...
             'Startpoint', [a*t x0 sig b*t]);

    % Gem fit-parametre (tag højde for tælletid)
    a_fit(j)   = f1.a/t;
    x0_fit(j)  = f1.b;
    sig_fit(j) = f1.c;
    b_fit(j)   = f1.d/t;

    % Plot fit kurve
    x_fit = x_min:0.05:x_max;
    y_fit = t*a_fit(j)*exp(-((x_fit-x0_fit(j)).^2)./(2*sig_fit(j)^2)) + t*b_fit(j);
    plot(x_fit,y_fit,'g-')
    hold off
    
end   % ... loop over forskellige fits

figure(2)
hist(a_fit)
xlabel('a-parameter (signal amplitude)','FontSize',14)

figure(3)
hist(b_fit)
xlabel('b-parameter (baggrund amplitude)','FontSize',14)

figure(4)
hist(x0_fit)
xlabel('x0-parameter (linie centrum)','FontSize',14)

figure(5)
hist(sig_fit)
xlabel('sigma-parameter (linie bredde)','FontSize',14)

% Formatér figurer
figure(2)
ah1 = gca;
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',12)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.4 0.4 0.8],'EdgeColor','k')

figure(3)
ah1 = gca;
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',12)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.4 0.4 0.8],'EdgeColor','k')

figure(4)
ah1 = gca;
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',12)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.4 0.4 0.8],'EdgeColor','k')

figure(5)
ah1 = gca;
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',12)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.4 0.4 0.8],'EdgeColor','k')


