function solution = planeShear_JR(param, xdisp)
% This function sloves the governing eqn derived from the full
% Jenkins & Richman theory in the dense limit for sheared flows between two
% parallel plates.
%
% It solves a set of ODEs:
%
%           dy/dx = f(x,y)
% where x is the distance between the parallel walls and y is a vector of
% unknowns:
%                   /x
%   y(1) = Phi(x) = |  phi(x')dx'
%                   /0
%
%   y(2) = d[y(1)]/dx = phi(x)
%
%   y(3) = u(x)
%   
%   y(4) = d[y(3)]/dx = du/dx
%
%   y(5) = w(x) = (T(x))^(1/2)
%
%   y(6) = d[y(5)]/dx = dw/dx
%
% The boundary conditions are given as:
%
% at x=0:   y(1) = 0, y(3) = u_bot, y(5) = w_bot;
% at x=1;   y(1) = phibar, y(3) = u_top, y(5) = w_top.
%
% syntax:
%   solution = planeShear_JR(param, xdisp)
%
% input:
%   param   --  A structure specifying the control parameters, including
%       param.e         --  coefficient of restitution
%       param.phibar    --  average volume fraction
%       u_bot           --  dimensionless velocity (T^(1/2))at the bottom (x=0)
%       u_top           --  dimensionless velocity at the top (x=1)
%       w_bot           --  dimensionless fluctuation velocity at x=0
%       w_top           --  dimensionless fluctuation velocity at x=1
%       Fr              --  streamwise Froude number (g*d/U^2)
%       H               --  normalized channel width (H/d)
%   xdisp   --  The points at which the solution will be evaluated 
%
% output:
%   solution    --  A structure contains the solution at points given by xdisp
%       solution.phi    --  volume fraction profile
%       solution.u      --  velocity profile
%       solution.w      --  fluctuation velocity (T^(1/2)) profile
%       solution.dudx   --  shear rate
%       solution.dwdx   --  gradient of T^(1/2)
%

e = param.e;            %   0.7;
phibar = param.phibar;  %   0.45;
u_bot = param.u_bot;    %   -0.2;
u_top = param.u_top;    %   0.2;
w_bot = param.w_bot;    %   0.4;
w_top = param.w_top;    %   0.4;
Fr = param.Fr;          %   0;
H = param.H;            %   8;         % dimensionless cell width

% some constants
E1 = 4/((1+e)*(3-e));
E2 = (2/5)*(1+e)*(3*e+1);
E3 = (8/25)*(1+e)^2*((3*e-1)+(12/pi)*(3-e));
E4 = 32/((1+e)*(49-33*e));
E5 = (3/5)*(1+e)*(1+e+2*e^2);
E6 = (2/25)*(1+e)^2*(9*(1+e)*(2*e-1)+(4/pi)*(49-33*e));
E7 = (24/25/pi)*H^2*(1+e)^2*(49-33*e)*(1-e);

% initial condition, given as simple shear flow
xin = [0:0.1:1];
Nx = length(xin);
solinit.x = xin;
solinit.y = [   phibar*xin
                phibar*ones(1,Nx)
                u_bot+(u_top-u_bot)*xin
                (u_top-u_bot)*ones(1,Nx)
                w_bot+(w_top-w_bot)*xin
                (w_top-w_bot)*ones(1,Nx)    ];


% solve the ODEs using bvp4c
sol = bvp4c(@planeShearODE, @planeShearBC, solinit);


% evaluate the solution at given points
Sx = deval(sol, xdisp);
solution.phi = (Sx(2,:))';
solution.u = (Sx(3,:))';
solution.w = (Sx(5,:))';
solution.dudx = (Sx(4,:))';
solution.dwdx = (Sx(6,:))';


% functions defining the ODEs
    function dydx = planeShearODE(x,y)
        dydx = [            y(2)
                    h1(y(2)).*y(6)./y(5)
                            y(4)
               h2(y(2)).*h1(y(2)).*y(6).*y(4)./y(5) - y(6).*y(4)./y(5) + Fr*h3(y(2))./y(5)
                            y(6)
               ((h4(y(2)).*h1(y(2)) - 1).*(y(6)).^2./y(5) + h5(y(2))*(y(4)).^2./y(5)) + h6(y(2))*y(5)     ];
    end

    function res = planeShearBC(ya, yb)
        res = [         ya(1)
                        yb(1) - phibar
                        ya(3) - u_bot
                        yb(3) - u_top
                        ya(5) - w_bot
                        yb(5) - w_top   ];
    end

% functions defined for the coefficients in the ODEs
    function f=h1(phi)
        fp = 1+2*(1+e)*phi.*chiCS(phi);
        f = -2*phi.*(1-phi).*(2-phi).*fp./ ...
            ((1-phi).*(2-phi)+(4-phi).*fp);
    end

    function f = h2(phi)
        chi = chiCS(phi);
        G = phi.*chi;
        f = (phi.*(5-2*phi)-E2*(1-phi).*(2-phi).*G-E3*(4-phi).*G.^2)./ ...
            (phi.*(1-phi).*(2-phi).*(1+E2*G+E3*G.^2));
    end

    function f = h3(phi)
        G = phi.*chiCS(phi);
        f = (96/(5*sqrt(pi)*E1))*G./(1+E2*G+E3*G.^2);
    end

    function f = h4(phi)
        G = phi.*chiCS(phi);
        f = (phi.*(5-2*phi)-E5*(1-phi).*(2-phi).*G-E6*(4-phi).*G.^2)./ ...
            (phi.*(1-phi).*(2-phi).*(1+E5*G+E6*G.^2));
    end

    function f = h5(phi)
        G = phi.*chiCS(phi);
        f = (4*E1/(15*E4))*(1+E2*G+E2*G.^2)./(1+E5*G+E6*G.^2);
    end

    function f = h6(phi)
        G = phi.*chiCS(phi);
        f = E7*(G.^2)./(1+E5*G+E6*G.^2);
    end

end