%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eksperimentel Fysik, Niels Bohr Institutet, Københavns Universitet
% Opgaver torsdag d. 26. april (uge 1), 2012
%
% Dette dokument indeholder løsninger til opgaverne fra bogen
%   J.R. Taylor, An Introduction to Error Analysis
% efter at opgaverne er blevet omskrevet til at blive lavet i MATLAB i
% stedet for i hånden
% 
% Opgaverne er uploadet i 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 (Taylor opgave 3.4) -------------------------------------------------

% Definerer datapunkterne (usikkerheden på tælletal er Poisson-fordelt)
t = 0:0.5:2;
nu = [214 134 101 61 54];
nu_usikk = sqrt(nu);

% Plotter data med errorbars
errorbar(t,nu,nu_usikk,'ro')
hold on  %denne gør at jeg bliver ved med at plotte oven i figuren senere

% Sætter akser, plot-titel og aksenavne
axis([0 2 40 240])
title('Opgave 1 - Radioaktivt henfald')
xlabel('Tid i timer')
ylabel('Tælletal')

% Laver teori-linien
tt = -0.25:0.01:2.25;
nu_0 = 200; r = 0.693;
nu_exp = nu_0 * exp(-r*tt);

% Plotter teori-linien oven i data-grafen
plot(tt,nu_exp,'b:')

% Sætter en forklaring på for at vise hvad der er hvad
legend('Datapunkter','Teori-linie')
% hold off  %denne "slipper" figuren igen


%% Opgave 2 (Taylor opgave 3.25) ------------------------------------------------

% Konklusionen i den fejlagtige udledning er forkert fordi man har brugt 
% udtrykket i formel uden at overholde dennes antagelser - nemlig at de
% variable man ganger sammen er uafhængige. Og nu er x jo fuldkommen afhængig 
% af sig selv - hvorfor konklusionen altså kun kan være forkert.

% Definerer liste med normalfordelte tal
x = 10 + randn(1,500);

% Ser på relativ usikkerhed for denne og for denne i anden
std(x)/mean(x)
std(x.^2)/mean(x.^2)     % = 0.2029

% Og man finder altså at nummer to usikkerhed ligner at man ganger den første 
% med 2, ikke sqrt(2)
2*std(x)/mean(x)         % = 0.2052
sqrt(2)*std(x)/mean(x)   % = 0.1451


%% Opgave 3 (Taylor opgave 4.6) -------------------------------------------------

% Definerer datasættet
nu = [10 13 8 15 8 13 14 13 19 8 13 13 7 8 6 8 11 12 8 7];

% Udregner gennemsnittet og SD
gn = sum(nu)/length(nu)   % = 10.7000
std(nu)                   % = 3.4043

% Udregner kvadratroden af gennemsnittet
sqrt(gn)                  % = 3.2711
% standardafvigelsen ovenfor er faktisk et over-estimat, men det passer ok


%% Opgave 4 (Taylor opgave 4.28) ------------------------------------------------

% Definerer datasættet
l = [51.2 59.7 68.2 79.7 88.3];
T = [1.448 1.566 1.669 1.804 1.896];

% Udregner hans værdier for g
g = 4*pi^2*l./T.^2

% Finder gennemsnittet, SD og SDOM af g-værdierne
g_mean     = mean(g)                 % = 965.6387
sigma      = std(g)                  % = 3.2557
sigma_mean = std(g)/sqrt(length(g))  % = 1.4560
% hans endelige resultat findes til at være  mean \pm SDOM = 966 \pm 1

% Definerer den accepterede værdi for g og finder diskrepansen
g_right = 979.6;
discrepancy = g_mean - g_right       % = -13.9613
% diskrepansen er næsten 14!

% Finder hvor stor fejlen skulle være for at den rigtige værdi er indenfor SDOM
x = 0.6:0.02:1.2;
for i=1:length(x)
    g_new = 4*pi^2*(l+x(i))./T.^2;
    discr_new(i) = mean(g_new) + std(g_new)/sqrt(length(g_new)) - g_right;
%   discr_new(i) = mean(g_new) - g_right;
end
plot(x,discr_new,'.b')
axis([0.6 1.2 -5 5])
xlabel('x')
ylabel('g teori - g fundet ved l+x')
% man kan se at den skærer 0 ca. omkring 1

% Retter den systematiske fejl med de 2.00 cm og finder den nye g-værdi
l = l + 1;
g = 4*pi^2*l./T.^2;
g_mean     = mean(g)                 % = 980.0811
sigma_mean = std(g)/sqrt(length(g))  % = 0.9585

g_mean - sigma_mean                  % = 979.1226


%% Opgave 5 (Taylor opgave 5.11) ------------------------------------------------

% Definerer min x-akse
x = -2:0.05:6;

% Definerer mine givne parametre
X1     = 2;
sigma1 = 1;
X2     = 3;
sigma2 = 0.3;

% Udregner værdierne der skal plottes
y1 = 1/( sigma1*sqrt(2*pi) ) * exp(-(x-X1).^2/(2*sigma1^2));
y2 = 1/( sigma2*sqrt(2*pi) ) * exp(-(x-X2).^2/(2*sigma2^2));

% Alternativt kan man definere en funktion
% MinGauss = @(x,a,b) 1/( b*sqrt(2*pi) ) * exp(-(x-a).^2/(2*b^2));
% og så bare skrive
% y1 = MinGauss(x,X1,sigma1);
% y2 = MinGauss(x,X2,sigma2);

% Plotter det hele sammen
plot(x,y1,'r-',x,y2,'b-')
xlabel('x')
ylabel('y')


%% Opgave 6 (Taylor opgave 5.31) ------------------------------------------------

% Definerer datasættet
t = [63 58 74 78 70 74 75 82 68 69 76 62 72 88 65 81 79 77 66 76 86 72 79 ...
     77 60 70 65 69 73 77 72 79 65 66 70 74 84 76 80 69];

% Standardafvigelsen findes til at være 7.0376 med
std(t)     % = 7.0376

% Udregner gennemsnit og SD for bins af størrelse 4
ncol = 4;
gn   = zeros(1,10);
sd   = zeros(1,10);

for n = 0:9
  gn(n+1) = sum(t(n*ncol+1:n*ncol + ncol))/ncol;
  sd(n+1) = std(t(n*ncol+1:n*ncol + ncol));
end

gn
sd
% Ud fra del (a) ville man forvente at SD roder lidt rundt om den SD man fandt 
% tidligere, og det er også det der sker

% For at plotte histogrammer der kan sammenlignes skal man bruge histc, hvor 
% man kan vælge bin-grænser manuelt. Først bestemmes egenskaberne for bins'ene
nbins   = 5;    %antal bins
min     = 58;   %position af første bin grænse
dist    = 6;    %bredde af bins
max     = min+nbins*dist;
edges   = min:dist:max;
centers = min+dist*0.5:dist:max+dist*0.5;

% Her udregnes histogrammerne
hist40  = histc(t,edges);
hist10  = histc(gn,edges);

% Plotter histogrammerne sammen
figure
bar(centers,hist40,0.5,'y')
hold on
bar(centers,hist10,0.3,'r')
xlabel('bin nummer')
ylabel('antal i bin')
legend('hist40','hist10')
% man kan tydeligt se at hist10 er mere centreret, som forventet

axis([58 94 0 13])