function [A_k, C_k, pre_k, A_l] = assembDeriv4(ksi, uab, ksi0, Q, R, ref_ome, inddq)    
    a = 0.9898;
    b = 0.0072;
    c = 0.0361;    
    d = 1.0;
    e = 0.0149;
    dt = 0.000125;
    
    Rs = 0.28;
    Ls = 0.003465;
    psi = 0.1989;
    B = 0;    
    kp = 1.5;
    pp = 4.0;
    J = 0.04;
    Lq = 1.0*Ls;
    Ld = 0.9*Ls;
    kpp = kp*pp*pp;
    kppj = kpp/J;
    
    ia = ksi(1);
    ib = ksi(2);
    ome = ksi(3);
    the = ksi(4);
    P = [ksi(5), ksi(6), ksi(8), ksi(11);...
         ksi(6), ksi(7), ksi(9), ksi(12);...
         ksi(8), ksi(9), ksi(10), ksi(13);...
         ksi(11), ksi(12), ksi(13), ksi(14)];    
    A_k = zeros(14);
    A_l = zeros(15);
       
    
    %puvodni matice derivaci
    if(inddq == 0)
        %stejne indukcnosti
        A = [a, 0, b*sin(the), b*ome*cos(the);...
             0, a, -b*cos(the), b*ome*sin(the);...
             -e*sin(the), e*cos(the), d, -e*(ib*sin(the)+ia*cos(the));...
             0, 0, dt, 1.0];
    else
        %ruzne indukcnosti
        A = [[ (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]];
    end
    C = [1, 0, 0, 0;... 
         0, 1, 0, 0];
    
    %dalsi matice EKF
    Pp = A*P*A' + Q;
    S = C*Pp*C' + R;
    K = Pp*C'/S;
    Pnew = (eye(4) - K*C)*Pp;    
    
    %derivace zakladnich matic EKF stavu
    dAdia = [0, 0, 0, 0;...
             0, 0, 0, 0;...
             0, 0, 0, -e*cos(the);...
             0, 0, 0, 0];
    dAdib = [0, 0, 0, 0;...
             0, 0, 0, 0;...
             0, 0, 0, -e*sin(the);...
             0, 0, 0, 0];
    dAdom = [0, 0, 0, b*cos(the);...
             0, 0, 0, b*sin(the);...
             0, 0, 0, 0;...
             0, 0, 0, 0];
    dAdth = [0, 0, b*cos(the), -b*ome*sin(the);...
             0, 0, b*sin(the), b*ome*cos(the);...
             -e*cos(the), -e*sin(the), 0, -e*(ib*cos(the)-ia*sin(the));...
             0, 0, 0, 0];
%     dAdP = zeros(4);
%     dCdiaibomth = zeros(2,4);
%     dCdP = zeros(2,4);
%     dPdiaibomth = zeros(4);
    dPdP5 = zeros(4);
    dPdP5(1,1) = 1;
    dPdP6 = zeros(4);
    dPdP6(1,2) = 1;
    dPdP6(2,1) = 1;
    dPdP7 = zeros(4);
    dPdP7(2,2) = 1;
    dPdP8 = zeros(4);
    dPdP8(1,3) = 1;
    dPdP8(3,1) = 1;
    dPdP9 = zeros(4);
    dPdP9(2,3) = 1;
    dPdP9(3,2) = 1;
    dPdP10 = zeros(4);
    dPdP10(3,3) = 1;
    dPdP11 = zeros(4);
    dPdP11(1,4) = 1;
    dPdP11(4,1) = 1;
    dPdP12 = zeros(4);
    dPdP12(2,4) = 1;
    dPdP12(4,2) = 1;
    dPdP13 = zeros(4);
    dPdP13(3,4) = 1;
    dPdP13(4,3) = 1;
    dPdP14 = zeros(4);
    dPdP14(4,4) = 1;
    
    
    %derivace dalsich matic EKF stavu
    dPpdia = dAdia*P*A' + A*P*dAdia';
    dPpdib = dAdib*P*A' + A*P*dAdib';
    dPpdom = dAdom*P*A' + A*P*dAdom';
    dPpdth = dAdth*P*A' + A*P*dAdth';
    dPpdP5 = A*dPdP5*A';
    dPpdP6 = A*dPdP6*A';
    dPpdP7 = A*dPdP7*A';
    dPpdP8 = A*dPdP8*A';
    dPpdP9 = A*dPdP9*A';
    dPpdP10 = A*dPdP10*A';
    dPpdP11 = A*dPdP11*A';
    dPpdP12 = A*dPdP12*A';
    dPpdP13 = A*dPdP13*A';
    dPpdP14 = A*dPdP14*A';
    
    dSdia = C*dPpdia*C';
    dSdib = C*dPpdib*C';
    dSdom = C*dPpdom*C';
    dSdth = C*dPpdth*C';
    dSdP5 = C*dPpdP5*C';
    dSdP6 = C*dPpdP6*C';
    dSdP7 = C*dPpdP7*C';
    dSdP8 = C*dPpdP8*C';
    dSdP9 = C*dPpdP9*C';
    dSdP10 = C*dPpdP10*C';
    dSdP11 = C*dPpdP11*C';
    dSdP12 = C*dPpdP12*C';
    dSdP13 = C*dPpdP13*C';
    dSdP14 = C*dPpdP14*C';
    
    dKdia = dPpdia*C'/S - Pp*C'/S*dSdia/S;
    dKdib = dPpdib*C'/S - Pp*C'/S*dSdib/S;
    dKdom = dPpdom*C'/S - Pp*C'/S*dSdom/S;
    dKdth = dPpdth*C'/S - Pp*C'/S*dSdth/S;
    dKdP5 = dPpdP5*C'/S - Pp*C'/S*dSdP5/S;
    dKdP6 = dPpdP6*C'/S - Pp*C'/S*dSdP6/S;
    dKdP7 = dPpdP7*C'/S - Pp*C'/S*dSdP7/S;
    dKdP8 = dPpdP8*C'/S - Pp*C'/S*dSdP8/S;
    dKdP9 = dPpdP9*C'/S - Pp*C'/S*dSdP9/S;
    dKdP10 = dPpdP10*C'/S - Pp*C'/S*dSdP10/S;
    dKdP11 = dPpdP11*C'/S - Pp*C'/S*dSdP11/S;
    dKdP12 = dPpdP12*C'/S - Pp*C'/S*dSdP12/S;
    dKdP13 = dPpdP13*C'/S - Pp*C'/S*dSdP13/S;
    dKdP14 = dPpdP14*C'/S - Pp*C'/S*dSdP14/S;    
    
    dPnewdia = dPpdia - dKdia*C*Pp - K*C*dPpdia;
    dPnewdib = dPpdib - dKdib*C*Pp - K*C*dPpdib;
    dPnewdom = dPpdom - dKdom*C*Pp - K*C*dPpdom; 
    dPnewdth = dPpdth - dKdth*C*Pp - K*C*dPpdth; 
    dPnewdP5 = dPpdP5 - dKdP5*C*Pp - K*C*dPpdP5;
    dPnewdP6 = dPpdP6 - dKdP6*C*Pp - K*C*dPpdP6;
    dPnewdP7 = dPpdP7 - dKdP7*C*Pp - K*C*dPpdP7;
    dPnewdP8 = dPpdP8 - dKdP8*C*Pp - K*C*dPpdP8;
    dPnewdP9 = dPpdP9 - dKdP9*C*Pp - K*C*dPpdP9;
    dPnewdP10 = dPpdP10 - dKdP10*C*Pp - K*C*dPpdP10;
    dPnewdP11 = dPpdP11 - dKdP11*C*Pp - K*C*dPpdP11;
    dPnewdP12 = dPpdP12 - dKdP12*C*Pp - K*C*dPpdP12;
    dPnewdP13 = dPpdP13 - dKdP13*C*Pp - K*C*dPpdP13;
    dPnewdP14 = dPpdP14 - dKdP14*C*Pp - K*C*dPpdP14;
    
    
    A_k(1:4,1:4) = A;
    %A_k(1:4,5:14) = zeros(4,10);
    
    pttrn = [1 5 6 9 10 11 13 14 15 16]'; %pattern for selecting proper elements from P matrices
    A_k(5:14,1) =  dPnewdia(pttrn);
    A_k(5:14,2) =  dPnewdib(pttrn);
    A_k(5:14,3) =  dPnewdom(pttrn);
    A_k(5:14,4) =  dPnewdth(pttrn);
    A_k(5:14,5) =  dPnewdP5(pttrn);
    A_k(5:14,6) =  dPnewdP6(pttrn);
    A_k(5:14,7) =  dPnewdP7(pttrn);
    A_k(5:14,8) =  dPnewdP8(pttrn);
    A_k(5:14,9) =  dPnewdP9(pttrn);
    A_k(5:14,10) =  dPnewdP10(pttrn);
    A_k(5:14,11) =  dPnewdP11(pttrn);
    A_k(5:14,12) =  dPnewdP12(pttrn);
    A_k(5:14,13) =  dPnewdP13(pttrn);
    A_k(5:14,14) =  dPnewdP14(pttrn);    
    
    phi = zeros(14,1);
    phi(1) = a*ksi0(1) + b*ksi0(3)*sin(ksi0(4)) + c*uab(1);
    phi(2) = a*ksi0(2) - b*ksi0(3)*cos(ksi0(4)) + c*uab(2);
    phi(3) = d*ksi0(3) + e*(ksi0(2)*cos(ksi0(4)) - ksi0(1)*sin(ksi0(4)));
    phi(4) = ksi0(4) + dt*ksi0(3);    
    phi(5:14) = Pnew(pttrn);
       
    A_l(1:14,1:14) = A_k;
    A_l(1:14,15) = phi - A_k*ksi0;
    A_l(15,15) = 1.0;
        
    C_k = zeros(2,14);
    C_k(1:2,1:4) = C;
    
    pre_k = Pnew(pttrn);      
    %max x(14) = pi^2/3 ... variance of uniform -pi,pi
    if(pre_k(10) > pi^2/3)
        pre_k(10) = pi^2/3;
    end
end