function [var_ome, var_th] = fpcrb2(model, control, amp, injom, Q, R, T, ref_profile, Q0)
    %function for computing PCRB of PMSM model
       

    %   machine parameters
    Rs = 0.28;
    Ls = 0.003465;
    psipm = 0.1989;
    B = 0;    
    kp = 1.5;
    pp = 4.0;
    J = 0.04;

    dt = 0.000125;

    Lq = 1.0*Ls;
    Ld = 0.9*Ls;

    %aliases
    kpp = kp*pp*pp;
    kppj = kpp/J;
    psi = psipm;

    a = 0.9898;
    b = 0.0072;
    c = 0.0361;
    d = 1.0;
    e = 0.0149;    

    %reference signal
    ref_ome = zeros(1, T);  
    for k = 1:T,
           index = floor(k*dt);
           if(index>0)
               lower = ref_profile(index);
           else
               lower = 0;
           end
           if(index<T*dt)
               upper = ref_profile(index+1);
           else
               upper = 0;
           end
           ref_ome(k) = lower + (upper-lower)*dt*(k-index/dt);
    end



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %   system variables

    %system state and control
    x_sys = zeros(4, T);     
    u_dq = zeros(2, T);

    %covariance alternatives
    Q4 = Q; 
    Q2 = Q4(3:4,3:4);
    iQ4 = inv(Q4);  
    iQ2 = inv(Q2);
    iR = inv(R);
    

    %PCRB J matrices
    if(model(1) == 1)
        iJn1 = zeros(2,T);
        Jj1 = inv(Q0);
    end
    if(model(2) == 1)
        iJn2 = zeros(2,T);
        Jj2 = inv(Q0);
    end
    if(model(3) == 1)
        iJn3 = zeros(2,T);
        Jj3 = inv(Q0(3:4,3:4));
    end
    if(model(4) == 1)
        iJn4 = zeros(2,T);
        Jj4 = inv(Q0(3:4,3:4));
    end
    
    

    %PI control
    sum_iq = 0;
    sum_ud = 0;
    sum_uq = 0;
    kon_pi = 30.0;
    kon_ii = 0.0;
    kon_pu = 20.0;
    kon_iu = 0.0;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %   init

    % measurement function derivatives
    C4 = zeros(2,4);
    C4(1,1) = 1.0;
    C4(2,2) = 1.0;

    % system function derivatives
        %alpha-beta, equal Ls
    A14 = zeros(4);
    A14(1,1) = a;
    A14(2,2) = a;
    A14(3,3) = d;
    A14(4,3) = dt;
    A14(4,4) = 1.0;       

    % new figure              
%     figure;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %   main loop

    for t = 1:T-1,
        %aliases
        ial = x_sys(1, t);
        ibe = x_sys(2, t);
        ome = x_sys(3, t);
        the = x_sys(4, t); 

        ia = ial;
        ib = ibe;

        id = ial*cos(the) + ibe*sin(the);
        iq = ibe*cos(the) - ial*sin(the);
        %PCRB models
        if(model(1) == 1) %al-be/Ls
            A14(1,3) = b*sin(the);
            A14(1,4) = b*ome*cos(the);
            A14(2,3) = -b*cos(the);
            A14(2,4) = b*ome*sin(the);
            A14(3,1) = -e*sin(the);
            A14(3,2) = e*cos(the);
            A14(3,4) = -e*(ial*cos(the)+ibe*sin(the));

            D11 = A14'*iQ4*A14;
            D12 = -A14'*iQ4;
            D21 = D12';
            D22 = iQ4 + C4'*iR*C4;

            Jj1 = D22 - D21/(Jj1 + D11)*D12;
            tmp = inv(Jj1);        
            iJn1(1,t) = tmp(3,3);
            iJn1(2,t) = tmp(4,4);
        end        
        if(model(2) == 1) %al-be/Ldq        
            A70 = [[ (Lq - Rs*dt*sin(the)^2)/Lq - (dt*ome*sin(the)*Lq^2*cos(the) + Rs*dt*Lq*cos(the)^2)/(Ld*Lq) + (Ld*dt*ome*cos(the)*sin(the))/Lq,                                       (dt*(Ld - Lq)*(- Lq*ome*cos(the)^2 + Rs*cos(the)*sin(the) + Ld*ome*sin(the)^2))/(Ld*Lq),    dt*cos(the)*(ia*sin(the) - ib*cos(the) + (Lq*(ib*cos(the) - ia*sin(the)))/Ld) + dt*sin(the)*(psi/Lq - ia*cos(the) - ib*sin(the) + (Ld*(ia*cos(the) + ib*sin(the)))/Lq),                                               (dt*(ome*psi*cos(the) + Rs*ib*cos(2*the) - Rs*ia*sin(2*the)))/Lq + (Ld*dt*(ia*ome*cos(2*the) + ib*ome*sin(2*the)))/Lq - (dt*(Lq^2*ia*ome*cos(2*the) + Lq^2*ib*ome*sin(2*the) + Lq*Rs*ib*cos(2*the) - Lq*Rs*ia*sin(2*the)))/(Ld*Lq)];...
                    [                                       (dt*(Ld - Lq)*(- Ld*ome*cos(the)^2 + Rs*cos(the)*sin(the) + Lq*ome*sin(the)^2))/(Ld*Lq), (Lq - Rs*dt*cos(the)^2)/Lq - (Lq*Rs*dt*sin(the)^2 - Lq^2*dt*ome*cos(the)*sin(the))/(Ld*Lq) - (Ld*dt*ome*cos(the)*sin(the))/Lq, (dt*(Lq*ia - psi*cos(the)))/Lq + (dt*((Lq^2*ia*cos(2*the))/2 - (Lq^2*ia)/2 + (Lq^2*ib*sin(2*the))/2))/(Ld*Lq) - (Ld*dt*(ia/2 + (ia*cos(2*the))/2 + (ib*sin(2*the))/2))/Lq, (dt*ome*psi*sin(the) - Rs*dt*ia*(2*sin(the)^2 - 1) + Rs*dt*ib*sin(2*the))/Lq + (Ld*(dt*ib*ome*(2*sin(the)^2 - 1) + dt*ia*ome*sin(2*the)))/Lq - (Lq*Rs*dt*ib*sin(2*the) + Lq^2*dt*ib*ome*(2*sin(the)^2 - 1) + Lq^2*dt*ia*ome*sin(2*the) - Lq*Rs*dt*ia*(2*sin(the)^2 - 1))/(Ld*Lq)];...
                    [     -dt*kppj*(psi*sin(the) - cos(the)*(Ld - Lq)*(ib*cos(the) - ia*sin(the)) + sin(the)*(Ld - Lq)*(ia*cos(the) + ib*sin(the))),      dt*kppj*(psi*cos(the) + cos(the)*(Ld - Lq)*(ia*cos(the) + ib*sin(the)) + sin(the)*(Ld - Lq)*(ib*cos(the) - ia*sin(the))),                                                                                                                                                                         1.0,                                                                                                                                                   -dt*kppj*(psi*(ia*cos(the) + ib*sin(the)) + (Ld - Lq)*(ia*cos(the) + ib*sin(the))^2 - (Ld - Lq)*(ib*cos(the) - ia*sin(the))^2)];...
                    [                                                                                                                             0.0,                                                                                                                             0.0,                                                                                                                                                                        dt,                                                                                                                                                                                                                                                                                1.0]];

            D11 = A70'*iQ4*A70;
            D12 = -A70'*iQ4;
            D21 = D12';
            D22 = iQ4 + C4'*iR*C4;

            Jj2 = D22 - D21/(Jj2 + D11)*D12;
            tmp = inv(Jj2);        
            iJn2(1,t) = tmp(3,3);
            iJn2(2,t) = tmp(4,4);        
        end
        if(model(3) == 1) %reduced al-be/Ls
            A3 = [d, -e*(ibe*sin(the)+ial*cos(the)); dt, 1];
            
            D11 = A3'*iQ2*A3;
            D12 = -A3'*iQ2;
            D21 = D12';
            D22 = iQ2 + [b^2, 0; 0, b^2*(ome^2 + Q4(3,3))]*iR;
%             D22 = iQ2 + [b^2, 0; 0, b^2*(ome^2)]*iR;
%             tmp = inv(Jj3);
%             D22 = iQ2 + [b^2, 0; 0, b^2*(ome^2 + tmp(1,1))]*iR;
% tady nahore to dela extremni rozdil
            
            
            Jj3 = D22 - D21/(Jj3 + D11)*D12;
            tmp = inv(Jj3);        
            iJn3(1,t) = tmp(1,1);
            iJn3(2,t) = tmp(2,2);
        end  
        if(model(4) == 1) %reduced al-be/Ldq
            A3 = [d, -dt*kppj*(psi*(ia*cos(the) + ib*sin(the)) + (Ld - Lq)*(ia*cos(the) + ib*sin(the))^2 - (Ld - Lq)*(ib*cos(the) - ia*sin(the))^2); dt, 1.0];
            
            D11 = A3'*iQ2*A3;
            D12 = -A3'*iQ2;
            D21 = D12';
            D22 = iQ2 + [b^2, 0; 0, b^2*(ome^2 + Q4(3,3))]*iR;
%             D22 = iQ2 + [b^2, 0; 0, b^2*(ome^2)]*iR;
%             tmp = inv(Jj3);
%             D22 = iQ2 + [b^2, 0; 0, b^2*(ome^2 + tmp(1,1))]*iR;
% tady nahore to dela extremni rozdil
            
            
            Jj4 = D22 - D21/(Jj4 + D11)*D12;
            tmp = inv(Jj4);        
            iJn4(1,t) = tmp(1,1);
            iJn4(2,t) = tmp(2,2);
        end 

        %control & simulation
        if(control == 1) %ome = ref_ome, the = integral(ome), ial=ibe=0
            x_sys(1, t+1) = 0;
            x_sys(2, t+1) = 0;
            x_sys(3, t+1) = ref_ome(t);
            x_sys(4, t+1) = the + dt*ome;
        else
            if(control == 6) %rnd. err. ref_ome
               ref_ome(t) = ref_ome(t) + amp*randn(); 
            end

            %PI (only PI for control == 2)        
            sum_iq = sum_iq + ref_ome(t) - ome;
            ref_iq = kon_pi*(ref_ome(t) - ome) + kon_ii*sum_iq;
            sum_ud = sum_ud - id;
            u_dq(1, t) = kon_pu*(-id) + kon_iu*sum_ud;
            sum_uq = sum_uq + ref_iq - iq;
            u_dq(2, t) = kon_pu*(ref_iq - iq) + kon_iu*sum_uq;
            u_dq(1, t) = u_dq(1, t) - Ls*ome*ref_iq;
            u_dq(2, t) = u_dq(2, t) + psipm*ome;        

            if(control == 3) %inj. sin -> dq
                u_dq(1,t) = u_dq(1,t) + amp*sin(injom*dt*t);
                u_dq(2,t) = u_dq(2,t) + amp*cos(injom*dt*t);
            elseif(control == 4) %inj. rect. -> dq
                u_dq(1,t) = u_dq(1,t) + amp*sign(sin(injom*dt*t));
                u_dq(2,t) = u_dq(2,t) + amp*sign(cos(injom*dt*t));
            elseif(control == 5) %inj. const. -> d
    %             u_dq(1,t) = 0.1*Ld/dt - (1.0*Ld/dt - Rs)*id - Lq*ome*iq;
                u_dq(1,t) = amp;            
            elseif(control == 9) %bic. sign. ome           
                u_dq(1,t) = u_dq(1,t) - amp*sign(ome);
                u_dq(2,t) = u_dq(2,t) + amp*sign(ome);
            end

            ual = u_dq(1,t)*cos(the) - u_dq(2,t)*sin(the);
            ube = u_dq(1,t)*sin(the) + u_dq(2,t)*cos(the);

            if(control == 7) %inj. sin -> al-be
                ual = ual + amp*cos(injom*dt*t);
                ube = ube + amp*sin(injom*dt*t);
            elseif(control == 8) %inj. rect. -> al-be
                ual = ual + amp*sign(cos(injom*dt*t));
                ube = ube + amp*sign(sin(injom*dt*t));        
            elseif((control == 10)&&(rand() > (1/5))) %bic rnd. 5
                duab = sign(rand(2,1)-0.5)*amp;
                ual = ual + duab(1);
                ube = ube + duab(2);
            end

            ud = ual*cos(the) + ube*sin(the);
            uq = ube*cos(the) - ual*sin(the);

            %simulation model Ldq
            idpl = (1.0 - Rs*dt/Ld)*id + Lq*dt/Ld*ome*iq + dt/Ld*ud;
            iqpl = (1.0 - Rs*dt/Lq)*iq - psipm*dt/Lq*ome - Ld*dt/Lq*ome*id + dt/Lq*uq;

            x_sys(1, t+1) = idpl*cos(the) - iqpl*sin(the);
            x_sys(2, t+1) = idpl*sin(the) + iqpl*cos(the);
            x_sys(3, t+1) = (1.0-B*dt/J)*ome + kp*pp*pp*dt/J*((Ld-Lq)*id*iq + psipm*iq);
            x_sys(4, t+1) = the + dt*ome;
        end

    end

    ssz = sum(model);
    var_ome = zeros(T,ssz);
    var_th = zeros(T,ssz);    
    ind = 1;
    if(model(1) == 1)
        var_ome(:,ind) = iJn1(1,:);
        var_th(:,ind) = iJn1(2,:);   
        ind = ind+1;
    end
    if(model(2) == 1)
        var_ome(:,ind) = iJn2(1,:);
        var_th(:,ind) = iJn2(2,:);
        ind = ind+1;
    end    
    if(model(3) == 1)
        var_ome(:,ind) = iJn3(1,:);
        var_th(:,ind) = iJn3(2,:);
        ind = ind+1;
    end
    if(model(4) == 1)
        var_ome(:,ind) = iJn4(1,:);
        var_th(:,ind) = iJn4(2,:);
        ind = ind+1;
    end
end