% This file solves the pendant drop problem for a fixed interfacial tension
function[parmp, solutionp]=pendantliquid(ns,deltap,K,zsign,parmguess,solguess)
%{
Input:
ns=number of steps along profile
deltap= vertical distance between the deformed membrane apex and the
     membrane fixed edge, normalized by Rm undeformed membrane
K = Rm, normalized by the capillary length
parmguess(1) = total contour length (units of Rm)
parmguess(2) = normalized pressure where 

Outputs:
volume: drop volume, normalized by Rm^3
area:  surface area of drop, normalized by Rm^2
shapeparm:  Neumann shape parameter
%}

options = bvpset('RelTol', 1e-5,'AbsTol',1e-8);
solinit.x=linspace(1/ns,1,ns);
solinit.y=solguess;
solinit.parameters=parmguess;

solb = bvp4c(@odebvp,@odebc,solinit,options); %apply boundary value solver

s = linspace(1/ns,1,ns);
solutionp = deval(solb,s);
parmp=solb.parameters;

%--------------------------------------------------------------
%define governing equations
    function dydx = odebvp(~,y,p)
        
        lp=p(1);
        p0p=p(2);
        
        a=y(1);
        r=y(2);
        z=y(3);
        
        % now define the set of relations in Eq. 12
        dydx = [lp*(p0p-sin(a)/r-zsign*K^2*z);
            lp*cos(a);
            lp*sin(a)];
    end
% set up the boundary conditions
    function res = odebc(ya,yb,p)
        
        lp=p(1);
        p0p=p(2);
        s1=1/ns;
        % now define the boundary conditions in Eq. 13
        res = [  ya(1)-0.5*(s1*lp*p0p);  % 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
    end

end