% The matrix of flow equations for the rigid flow 
% the components are [Kc vc y Db Df w M L]

function dA = rigdFullMatrix(l,A1)
    A = real(A1);
    Kc = A(1); vc = A(2); Ks = A(3); vs = A(4); 
    y = A(5); Db = A(6); w = A(7); lam = A(8); L = A(9);
        
    F = 1./sqrt((1-exp(-w)));
    g = 1 + (2./w)*log(1+1./F);
    I1 = Iw1(l,w);
    I2 = IM1(l,w);
        
    r1 = (vs/vc)^(2-Kc);
    r2 = (vs/vc)^(-Kc);
    r3 = 1+(1/Kc^2);
    r4 = 1-(1/Kc^2);
    r3s = 1+(1/Ks^2);
    r4s = 1-(1/Ks^2);
    f  = -0.25*(r1*r4*Kc+r2*r4s*Ks);
    
    dKc = -0.5*(vc/vs)*(Kc^2)*Db*(F*(1-0.5*r1*r3)-g*0.5*r1*r3*(2-Kc-Ks-y));
    dvc = -0.5*(Kc/vs)*(vc^2)*Db*(F*(1+0.5*r1*r4)+g*0.5*r1*r4*(2-Kc-Ks-y));
    dKs = -0.5*Ks^2*(y^2+Db*F*(1-0.5*r2*r3s)-Db*g*r2*(2*y+0.5*r3s*(2-Kc-Ks-y)));
    dvs = -0.5*Ks*vs*Db*(F*(1+0.5*r2*r4s)+g*0.5*r2*r4s*(2-Kc-Ks-y));
    dy  = (2-2*Ks)*y - Db*F*(1-r2*(1+0.25*r4s*y*Ks)) + Db*g*r2*(0.5*r3s*y*Ks^2+0.25*r4s*y*Ks-Kc+Ks-y);
    dDb = (3-Kc-Ks-y)*Db;
    dw  = w + L*Db*((1-cosh(w))*(1-f)*F^3-(2-Kc-Ks-y)*f*w*I1);
    dlam  = -lam*L*Db*((1-cosh(w)+w*sinh(w))*(1-f)/w-(2-Kc-Ks-y)*f*I2);
    dL  = -L;
    
    params = [l Kc Db y]

dA = [dKc; dvc; dKs; dvs; dy; dDb; dw; dlam; dL];