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

function dA = flopFullMatrix(l,A1)
    A = real(A1);
    Kc = A(1); vc = A(2); Ks = A(3); vs = A(4); 
    y = A(5); Db = A(6); s = A(7); lam = A(8);
        
    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*(1-(2/3)*r1*r3*(2.75-Kc-Ks-y));
    dvc = -0.5*(Kc/vs)*(vc^2)*Db*(1+(2/3)*r1*r4*(2.75-Kc-Ks-y));
    dKs = -0.5*Ks^2*(y^2+Db*(1-(2/3)*r2*(r3s*(2.75-Kc-Ks-y)+4*y)));
    dvs = -0.5*Ks*vs*Db*(1+(2/3)*r1*r4s*(2.75-Kc-Ks-y));
    dy  = (2-2*Ks)*y - Db*(1-r2*(1+(1/3)*r4s*y*Ks*(2.75-Kc-Ks-y)+(2/3)*r3s*y*Ks^2-(4/3)*(Kc-Ks+y)));
    dDb = (2.75-Kc-Ks-y)*Db;
    dlam  = -lam*sqrt(s)*Db*((1-f)+(4/7)*f*(2-Kc-Ks-y)*(1-exp(-7*l/4)));
    ds  = -s*(1-dlam/lam);
    
    params = [l Kc Db y]

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