root/applications/dual/SIDP/SIDPS/matlab/rizeni.m @ 1185

function  [ztrata u y b P]=rizeni(H,C,system,apriori,seed)
%RAND
sigma=system.sigma;
horizont=system.horizont;

y=zeros(1,horizont);
b=zeros(1,horizont);
P=zeros(1,horizont);
u=zeros(1,horizont-1);
ztrata=0;

y(1)=apriori.y0;
b(1)=apriori.b0;
P(1)=apriori.P0;

    for i=1:horizont
        eta=abs(y(i))/sigma;
        beta=abs(b(i))/sqrt(P(i));
                
        %plot(H,'b.'); hold on; plot(index, beta,'r.'); hold off           
        
        ny=eta*beta/(1+beta^2)+interpol(beta,H,C);
        %ny=(eta*beta/(1+beta^2)+C(index)*abs(rand));                
        znamenko=sign(y(i)*((-1)^(b(i)<0)));
        u(i)=-ny*sigma/sqrt(P(i))*znamenko;
        
        y(i+1)=y(i)+system.b*u(i)+sigma*seed(i);
        K=u(i)*P(i)/(P(i)*u(i)^2+sigma^2);
        b(i+1)=b(i)+K*(y(i+1)-y(i)-b(i)*u(i));
        P(i+1)=(1-K*u(i))*P(i);   
        
        ztrata=ztrata+y(i+1)^2;
    end
    
 end