program Kalman_Filter;

{$APPTYPE CONSOLE}

uses
  SysUtils, math, ShellAPI;

// predpoklad start k 21.6, max. delka dne 21.6 00 h
// omega hodinovy uhel, perioda 24 hod, nulovy 21.6 12 h
// Vyska - vyjadruje vysku Slunce nad obzorem v rozmezi -1 az +1, pokud je pod 0, je insolace nulova
// Neuvazovana obllacnost, clenitost, propustnost atmosfery, nadmorska vyska, advekce, srazky, latentni teplo

var P,Q,R,K,a,b,X,Ti,Tj,Tk,sigma,phi,dec,omega,omegas,Vyska,den: real;
var time: integer;

begin

Ti:=300; //pocetecni teplota v K
P:=0.6; // Chyba predpovedi
Q:=0.5; // Chyba modelu
R:=0.2; // Chyba mereni
a:=1;
X:=300; // Radiacni bilance Insolace - Vyzarovani W/m2
//b:=0.003*X; // X je intenzita zmeny teploty v zavislosi
time:=0;
sigma:=5.6704E-8;
phi:=0.89011792; // 51 st.
writeln('Den',';','Teplota',';','sin elevace',';','dec sun',';','sinus omega',';','zareni SV');

for time:=0 to 26299 do begin  // beh 26299 hodin, 3 roky

dec:=23.5*cos(6.2831853071796*time/8766.24);   // deklinace slunce
omega:=(6.2831853071796*(time-12)/24);  // hodinovy uhel, bezi od zacatku
Vyska:=sin(phi)*sin(dec)+cos(phi)*cos(dec)*cos(omega); // vyska jako sinus elevace alpha=(90-theta), theta je zenitovy uhel
if Vyska>0 then X:=(1368*Vyska)-(sigma*Tj*Tj*Tj*Tj) else X:=-(sigma*Tj*Tj*Tj*Tj); // bilance = insolace - vyzareni, pokud vyska<0 tak insolace = 0
b:=0.001*X;   // model
Tj:=a*Ti+b;  // model
P:=P+Q;      // Chyba predpoveni nova = Chyba predpovedi P + chyba modelu Q
K:=P/(P+R);  // zmena vahy K (podle chyby mereni R a chyby predpovedi P)
Tk:=Tj+K*(Ti-Tj); // nova hodnota
P:=(1-K)*P;   // oprava chyby predpovedi
Ti:=Tk;
den:=time/24;
omegas:=sin(omega);
writeln(den,';',Tk,';',Vyska,';',dec,';',omegas,';',X);

end;
end.

