% This file solves the pendant drop problem for a fixed interfacial tension
function[parmsol, solution]=pendantsolid(ns,Kp,zsign,l0,kappa,deltap,parmguess,solguess)
%{
Input:
ns = number of points along arc
Kp = Rm, normalized by the capillary length for undeformed case
l0 = extension ratio at drop apex
kappa = normalized dilational modulus
parmguess = 4 element vector with guesss for input parameters to be returned by subroutie (see outputs)
solguess = 7xnspoint matrix containing guessed vales of a, z, lxi, lphi

Outputs:
parmsol = 3 element array with the following quantities:
  
solution = 4xnspoint matrix containing the following quantities as a function of sp,
the position along the normalized undeformed arc
1) ap:  angle 
2) rp:  radial distance
3) zp:  vertical distance
4) 
5)
6)
7)
ap,rp,zp, a, z, lxi, lphi
%}
s1=1/ns; 
sp = linspace(s1,1,ns);
options = bvpset('RelTol', 1e-5,'AbsTol',1e-8);
solinit.x=sp;
solinit.y=solguess;
solinit.parameters=parmguess;

solb = bvp4c(@odebvp,@odebc,solinit,options); %apply boundary value solver
solution = deval(solb,sp);
parmsol=solb.parameters;

%--------------------------------------------------------------
%define governing equations
    function dydx = odebvp(~,y,p)
        lp=p(1);
        p0p=p(2);
        p0=p(3);

        % x here is the sp - the fractional arc position in the undeformed
        % case
        ap=y(1); % alpha in undeformed state
        rp=y(2); % radius in undeformed state
        zp=y(3); % height in undeformed state
        a=y(4);  % angle
        z=y(5);  % z coordinate
        lxi=y(6);  % longitudinal extenion ratio (psi direction)
        lphi=y(7);  % transverse extension ratio (phi direction)

        % Note:  all tensions are normalized here by kappa
        txi=1+(kappa/3)*(lxi/lphi-1/(lxi^3*lphi^3));  % tension in psi direction
        tphi=1+(kappa/3)*(lphi/lxi-1/(lxi^3*lphi^3));  % tension in phi direction
        txixi=(kappa/(3*lxi))*(lxi/lphi+3/(lxi^3*lphi^3));  % derivative of txi wrt lxi
        txiphi=(kappa/(3*lxi))*(-lxi^2/lphi^2+3/(lxi^2*lphi^4)); %% derivative of txi wrt lphi       
        dydx = [ p0p*lp-zsign*Kp^2*zp*lp-lp*sin(ap)/rp;  %  ap
             lp*cos(ap);    %rp
             lp*sin(ap);  %zp
             (lxi*lp/txi)*(p0-zsign*z*Kp^2-tphi*sin(a)/(lphi*rp));  % a
                 lxi*lp*sin(a);   % z
             (lp*lxi*cos(a)*(tphi-txi)-lp*txiphi*lphi*(lxi*cos(a)-lphi*cos(ap)))/(txixi*rp*lphi);  % lxi
             (lp*lxi*cos(a)-lp*lphi*cos(ap))/rp;];  %  lphi
            
    end

% set up the boundary conditions
    function res = odebc(ya,yb,p)
        lp=p(1);
        p0p=p(2);
        p0=p(3);
        res = [ya(1)-0.5*p0p*s1*lp;                    % ap(s1)
            ya(2)-s1*lp;                            % rp(s1)
            ya(3)-0.25*s1^2*lp^2*p0p;               % zp(s1)
            yb(2)-1;                                % rp(1)
            yb(3)-deltap;                          % zp(1)
            ya(4)-0.5*p0*s1*lp;                  % a(s1)
            ya(5)-0.25*s1^2*lp^2*p0;                % z(s1)
            ya(6)-l0;                               % lxi(s1)
            ya(7)-l0;                               % lphi(s1)
            yb(7)-1;];                              % lphi(1)
     end
end