clc;
clear all;
close all;

%definice hodnot
M=10e6;                %pocet opakovani
R0=99.989;%ohm         %hodnota odporu pri 0°C
Utro=0.29e-3;%ohm      %opakovatelnost referencniho teplomeru (RT)
Utrn=10.8e-3;%ohm      %nejistota RT
Utrh=2e-3;%ohm         %hystereze RT
Utrd=4.8e-3;%°C        %dlouhodoby drift RT
Utrs=4e-3;%°C          %samoohřev RT
Utko=0.27e-3;%ohm      %opakovatelnost kalibrovaneho teplomeru (KT)
Utkh=3e-3;%ohm         %hystereze KT
Utks=4e-3;%°C          %samoohrev KT
Umrn=1.99e-3;%ohm      %nejistota referencniho multimetru (RM)
Umrd=0.29e-3;%ohm      %drift RM
Umrr=0.01e-3;%ohm      %rozliseni RM
Umrt=20e-3;%°C         %parazitní napětí RM 
Umkn=15.95e-3;%ohm     %nejistota kalibrovaneho multimetru (KM)
Umkr=0.1e-3;%ohm       %rozliseni KM
Umkt=20e-3;%°C         %parazitní napětí KM 
Ulh=30e-3;%°C          %homogenita lazne
Uls=5e-3;%°C           %stabilita lazne
A=3.9083e-3;%°C-1      %koeficient A
B=-5.775e-7;%°C-2      %koeficient B
t = 50;%°C             %teplota 20 stupnu 

%generovani nahodnych hodnot pro zdroje nejistoty pusobici na R0[Ohm]
rUtro = random('norm',R0,Utro,1,M);       %opakovatelnost RT, normalni rozdeleni, stredni hodnota R0
rUtrn = random('norm',0,Utrn,1,M);        %nejistota RT, stredni hodnota 0 
rUtrh = 0 - Utrh + (2*Utrh).*rand(1,M);   %hystereze RT, rovnomerne rozdeleni
rUtko = random('norm',0,Utko,1,M);        %opakovatelnost KT
rUtkh = 0 - Utkh + (2*Utkh).*rand(1,M);   %hystereze KT
rUmrn = random('norm',0,Umrn,1,M);        %nejistota RM
rUmrr = 0 - Umrr + (2*Umrr).*rand(1,M);   %rozliseni RM
rUmrd = 0 - Umrd + (2*Umrd).*rand(1,M);   %drift RM
rUmkn = random('norm',0,Umkn,1,M);        %nejistota KM
rUmkr = 0 - Umkr + (2*Umkr).*rand(1,M);   %rozliseni KM

%soucet vygenerovanych nahodnych hodnot R0[Ohm]
soucetR0 = rUtro + rUtrn + rUtrh + rUtko + rUtkh + rUmrn + rUmrr + rUmrd + rUmkn + rUmkr;

%vykresleni histogramu pro R0
figure(1);                                                                    
hist(soucetR0,100);
xlabel('odpor R0 [Ohm]');
ylabel('cetnost [-]');
title('Histogram odporu R0');
grid on;

%generovani nahodnych hodnot pro zdroje nejistoty pusobici na t[°C]
rUtrd = 0 - Utrd + (2*Utrd).*rand(1,M);   %drift RT
rUtrs = 0 - Utrs + (2*Utrs).*rand(1,M);   %samoohrev RT
rUtks = 0 - Utks + (2*Utks).*rand(1,M);   %samoohrev KT
rUmrt = 0 - Umrt + (2*Umrt).*rand(1,M);   %parazitni napeti RM  
rUmkt = 0 - Umkt + (2*Umkt).*rand(1,M);   %parazitni napeti KM
rUlh = 0 - Ulh + (2*Ulh).*rand(1,M);      %homogenita lazne
rUls = 0 - Uls + (2*Uls).*rand(1,M);      %stabilita lazne

%soucet vygenerovanych nahodnych hodnot t[°C]
soucett = t + rUtrd + rUtrs + rUtks + rUmrt + rUmkt + rUlh + rUls;

%vykreslení histogramu pro t
figure(2);                                                                    
hist(soucett,100);
xlabel('teplota t[°C]');
ylabel('cetnost [-]');
title('Histogram teploty t');
grid on;

%model mereni
Rt = soucetR0.*(1+A*soucett+B*(soucett.^2));
%Rt = soucetR0.*(1+A*soucett);

%vykresleni histogramu pro Rt
figure(3);                                                                    
hist(Rt,100);
xlabel('odpor Rt [Ohm]');
ylabel('cetnost [-]');
title('Histogram výsledného odporu Rt');
grid on;

%prepocet na teplotu
D = A^2-4*B*(1-(Rt./R0));
tvys = (-A+D.^0.5)/(2*B);
%tvys = (Rt - R0)/(R0*A);

%prumerna hodnota potrebna k vypoctu vysledne nejistoty
prumtvys = mean(tvys);

%vykresleni histogramu pro tvys
figure(4);                                                                    
hist(tvys,100);
xlabel('teplota tvys [°C]');
ylabel('cetnost [-]');
title('Histogram výsledné teploty t');
grid on;

%serazeni do neklesajiciho poradi
tvys = sort(tvys);

%vypocet smerodatne odchylky
%uA=std(tvys);

%vypocet intervalu pokryti 95% pro tvys
p = 0.95;
q = p*M;
r = (M-q)/2;
f=r+q;
tmin = tvys(r);
tmax = tvys(r+q);

%vypocet vysledne nejistoty
U = tmax - prumtvys;

%vypocet mezi a Gaussova krivka z GUMu
Ugum = 0.07036;                     %hodnota kombinovane nejistoty vypocitana GUMem 
Uso=Ugum/2;
tmingum = t - Ugum;
tmaxgum = t + Ugum;
GUM = random('norm',t,Uso,1,M);     %generovani Gaussovy krivky

%histogram kumulativni distribucni funkce
figure(5);
[F,X]=ecdf(tvys);
ecdfhist(F,X,100);
title('Histogram kumulativní distribuční funkce výsledné teploty t');
xlabel('teplota tvys [°C]');
ylabel('cetnost [%]');
grid on;

%vykresleni jednotlivych intervalu 
line([tmin tmin],[0 12],'LineWidth',2,'Color','b');
hold on;
legend('','Monte Carlo');
line([tmaxgum tmaxgum],[0 12],'LineWidth',2,'Color','r','LineStyle','--');
legend('','Monte Carlo','GUM');
hold on;
line([tmax tmax],[0 12],'LineWidth',2,'Color','b');
hold on;
line([tmingum tmingum],[0 12],'LineWidth',2,'Color','r','LineStyle','--');
hold on;

% vykresleni prubehu Gaussovy krivky z GUMu do grafu 
[F,X]=ecdf(GUM);
[N,C]=ecdfhist(F,X,100);
[mu,var]=normfit(GUM);
t=49.8:0.001:50.2;
y=normpdf(t,mu,var);
plot(t,y,'r','LineWidth',2);
hold on;