clf     // czyszczenie okna graficznego
clc     // czyszczenie okna polecen
R1n=1000;R2n=1000;R3n=1000;kn=10;   // wartości nominalne parametrów


function x = RweNIC(R1,R2,R3,k)
// funkcja układowa określejąca zależność na Rwe
// Rwe(R1,R2,R3,k)
// parametry wpływowe: R1, R2, R3, k
x=-R1*((R3*(k+1)+R2))/(R2*(k-1)-R3);
endfunction

Rwen=RweNIC(R1n,R2n,R3n,kn)         // wartość nominalna Rwe
disp(['Rwen=',string(Rwen)])


// współczynniki wrazliwosci:
dY=0.001                        // względne zaburzenie parametru wpływowego
SR1=(RweNIC(R1n*(1+dY),R2n,R3n,kn)-Rwen)/(dY*Rwen)  // wsp. SR1
SR2=(RweNIC(R1n,R2n*(1+dY),R3n,kn)-Rwen)/(dY*Rwen)  // wsp. SR2
SR3=(RweNIC(R1n,R2n,R3n*(1+dY),kn)-Rwen)/(dY*Rwen)  // wsp. SR3
Sk=(RweNIC(R1n,R2n,R3n,kn*(1+dY))-Rwen)/(dY*Rwen)   // wsp. Sk
disp(['SR1=',string(SR1)])
disp(['SR2=',string(SR2)])
disp(['SR3=',string(SR3)])
disp(['Sk=',string(Sk)])

// względna zmiana Rwe dla zadanych względnych odchyłek parametrów
// wpływowych:
dR1=0.02;dR2=-0.01;dR3=0.03;dk=-0.05;       // względne odchyłki parametrów
dRwe=dR1*SR1+dR2*SR2+dR3*SR3+dk*Sk          // względna zmiana Rwe
disp(['względna oochyłka rezystancji dRwe=',string(dRwe)])


// względna zmiana Rwe spowodowana zmianą spadkiem temp. o DT = -30:
TWR = 50*10^-6; TWk=-150*10^-6;             // współczynniki temperaturowe
DT = -30;                                   // bezwzględna zmiana temperatury
dR=TWR*DT; dk=TWk*DT;           // termiczne wzgledne odchyłki R i k
dRwet=dR*(SR1+SR2+SR3)+dk*Sk    // termiczna względna zmiana Rwe
disp(['względna oschyłka termiczna rezystancji przy spadku temp.  o 30st dRwet=',string(dRwet)])


// względna skrajna odchyłka Rwe metodą najgorszego przypadku (worst case):
dR1tl=0.05;dR2tl=0.02;dR3tl=0.01;dktl=0.1;   // tolerancje parametrów R i k
dRwenp=dR1tl*abs(SR1)+dR2tl*abs(SR2)+dR3tl*abs(SR3)+dktl*abs(Sk)
disp(['względna oschyłka rezystancji met najgorszego przypadku dRwenp=',string(dRwenp)])

// statystyczna analiza wrażliwości
m = 1000      // liczba losowań
dR1norm = dR1tl*rand(m,1,"normal")/3;    // seria dR1 (rozkład normalny)
dR2norm = dR2tl*rand(m,1,"normal")/3;    // seria dR2 (rozkład normalny)
dR3norm = dR3tl*rand(m,1,"normal")/3;    // seria dR3 (rozkład normalny)
dknorm = dktl*rand(m,1,"normal")/3;      // seria dk (rozkład normalny)
dRwenorm = dR1norm*SR1 + dR2norm*SR2 + dR3norm*SR3 + dknorm*Sk;   // rozkład normalny m odchylek dRwe
Rwenorm = Rwen*(1+dRwenorm);            // rozkład normalny m odchylonych wartości Rwe(j)
Rweavg = mean(Rwenorm)                  // wartość średnia
sigma = st_deviation(Rwenorm,1)                  // odchylenie standardowe
disp(['wartość średnia rozkladu dRweavg=',string(Rweavg)])
disp(['odchylenie sredniokwadratowe rozkladu sigma=',string(sigma)])
disp(['niepewnośc graniczna wzgledna z poziomem ufności 68,3% dRwe68',string(sigma/Rweavg)])
disp(['niepewnośc graniczna wzgledna z poziomem ufności 99,7% dRwe99',string(3*sigma/Rweavg)])
//hist(Rwenorm)