% The matrix of flow equations for the gaussian flow 
% the components are y(1)=Kc(l) y(2) = vc(l) y(3)=y(l) y(4)=eta(l) y(5)=lambda(l) 
% y(6) = D(k1) y(7) = D(k2) ...
function dA = gaussRedFlowMatrix(l,A1)
    A = real(A1);
    Kc = A(1); vc = A(2); y = A(3); eta = A(4); lambda = A(5);
    Db = A(6:length(A)); vf =1;
    N = length(A)-5;

    k = (2*pi/N).*(0:N-1)';
    r1 = (vf/vc)^(2-Kc);
    r2 = (vf/vc)^(-Kc);
    r3 = 1+(1/Kc^2);
    r4 = 1-(1/Kc^2);   
    r5 = 1-((k.^2*lambda^2)/(2*pi));
    r6 = 1-exp(-l.*r5);
    
    B = lambda;
    qKc = (B/N)*sum(Db.*(1-0.5*r1*r3*(1-Kc-(3/2)*y+r5)./r5));
    qvc = (B/N)*sum(Db.*(1+0.5*r1*r4*(1-Kc-(3/2)*y+r5)./r5));
    qy =  (B/N)*sum(Db.*(1-r2*(1-Kc-(3/2)*y+r5)./r5));
    qeta =(B/N)*sum(Db.*(k.^2).*(1+0.25*(r1*Kc+r2*y)).*(1+(r5-Kc-(3/2)*y).*r6./r5));
    
    dKc = -0.5*(vc/vf)*(Kc^2)*qKc;
    dvc = -0.5*(Kc/vf)*(vc^2)*qvc;
    dy  = -y^2 - qy;
    dDb = (1-Kc-(3/2)*y+r5).*Db;
    deta = -lambda^2*eta^2/(4*pi)*qeta;
    dlambda = 0.5*(lambda/eta)*deta;
    
    par = [Kc eta lambda Db(1) r5(end)]
    %if r5==0
    %    dlmwrite('r5Log.m',[l min(abs(r5)) max(abs(r5)) lambda Db(1) Kc],'-append','roffset',1,'delimiter','\t');
    %end
    
dA = vertcat([dKc; dvc; dy; deta; dlambda],dDb);
