%  clear all  % clear all previously defined variables
clear all
deltaptarget=8; % normalized drop drop displacement before deformation
ns=1000; % number of steps in arc length
K=0.1; % shape parameter for droplet before deformation
l0target=0.9;  % stretch ratio after deformation at droplet apex
kappa=3;  % small strain dilational modulus, from Neo-Hookean strain energy function
zsign=1;  % 1 for pendant drop, -1 for pendant bubble
% the next two parameters can be adjusted to integral numbers larger than 1
% to make sure that sensible initial guesses are always obtained
ndeltapoints=1;  % number of steps used to reach target delta
nl0points=3;  % number of steps used to reach target deformation
plotsize=280;  % plotsize in points
%
deldelta=deltaptarget/ndeltapoints;
dell0=(l0target-1)/nl0points;
%
if zsign==1
    droptype='drop';
elseif zsign==-1
    droptype='bubble';
else
    pause 'not a valid value for zsign'
end
% parmpsol(1) is the overall countour length, normalized by Rm
% parmpsol(2) is the normalized pressure
% our initial guess, used only for the lowest value of delta, is a
% spherical cap
deltaiter=1;
deltap=deldelta;
s = linspace(1/ns,1,ns);  % this sets the original mesh of points equally spaced along the undeformed arc length
% in the following initial guess we assume that the membrane is deformed to a constant
% radius of curvature, which is determined by deltap
amax=2*atan(deltap);  % this is the angle at s=1 (R=Rm)
ap=s*amax; 
rp=s.*sinc(ap/pi);
zp=tan(ap/2);
solp(1,:)=ap;  % these are the angles along the whole membrane profile for the assumed cirular cross section
solp(2,:)=rp;  % these are the r values
solp(3,:)=zp;  % these are the z values
parmpsol(1)=1/sinc(amax/pi);  % initial guess for lp (normalized arc length of membrane)
parmpsol(2)=4*deltap/(1+deltap^2);  % inital guess for normalized pressure
% in the following loop we solve the pendant drop problem for a constant
% interfacial tension
while deltaiter<=ndeltapoints;
    parmpguess=parmpsol;
    solpguess=solp;
        [parmpsol,solp]=...
            pendantliquid(ns,deltap,K,zsign,parmpguess,solpguess);         
    deltaiter=deltaiter+1;
    deltap=deltap+deldelta;
end
deltap=deltap-deldelta;
lp=parmpsol(1);  % this is the total arc length, dived by Rm
p0p=parmpsol(2); % this is the normalized pressure at the apex
fprintf(1,'liquid case completed \n') %#ok<PRTCAL>
% the solution to the undeformed case is alrady known:
ap(:)=solp(1,:);  % ap, rp and zp are the solutions to the liquid problem
rp(:)=solp(2,:);
zp(:)=solp(3,:);
sol(1,:)=ap;  % the solutions to the liquid problem are needed
sol(2,:)=rp;  % for the full elastic problem - the program is set up 
sol(3,:)=zp;  % to solve 7 simultaneous first order differential equations,
              % although we already know the solutions to the first three
parmsol(1)=lp; 
parmsol(2)=p0p;
% For the first guess for the deformed profile, we assume a guess
% consistent with the inflation of a flat membrane to a constant radius of
% curvature
l0=1;
l0iter=1;
amaxinc=2*atan(dell0);  % this is the angle at s=1 (R=Rm)
a=ap+s*amaxinc; % these are the initial guesses
z=zp+tan(s*amaxinc/2);  % for a z and z for the eastic case
lxi(1:ns)=1+dell0;  % assume constant lxi for initial guess
lphi(1:ns)=1+dell0.*(1-s); % assume linear variation in lphi for initial guess
sol(4,:)=a(:); 
sol(5,:)=z(:);  
sol(6,:)=lxi(:);  
sol(7,:)=lphi(:);  
parmsol(3)=parmpsol(2)+4*dell0/(1+dell0^2);
% in the following loop we solve the full elastic problem
while l0iter<=nl0points;
    l0=l0+dell0;
    parmguess=parmsol;
    solguess=sol;   
        [parmsol, sol]=pendantsolid(ns,K,zsign,l0,kappa,deltap,parmguess,solguess);
    fprintf(1,'l0iter= %4.2u \n',l0iter);
    l0iter=l0iter+1;
end
fprintf(1,'solid case completed \n') %#ok<PRTCAL>
ap=[0,sol(1,:)];
rp=[0,sol(2,:)];
zp=[0,sol(3,:)];
a=[0,sol(4,:)];
z=[0,sol(5,:)];
lxi=[l0,sol(6,:)];
lphi=[l0,sol(7,:)]; 
r=lphi.*rp;
s=[0,s];
p0=parmsol(3); % this is the normalized pressure
%now make the plots
if ~exist('data', 'dir')
  mkdir('data');  % make the data directory to story the plots
end
set(0, 'DefaultAxesFontSize', 16);
set(0, 'DefaultLineLineWidth', 2);
% add the points at S=0 to the data
% now we generate the plots for the underformed and deformed drops
% Note that for the plots, we redefine z=0 at the the interface with the
% capillary
zrange=max(zp(ns+1),z(ns+1));
rrange=max(max(rp),max(r));
hold off
fig(1)=plot(rp,zp-zp(ns),'--b');
hold on
plot(-rp,zp-zp(ns),'--b');
fig(2)=plot(r,z-z(ns),'-r');
plot(-r,z-z(ns),'-r');
xlabel('$\overline{r}$','interpreter','latex','fontsize',24)
ylabel('$\overline{z}$','interpreter','latex','fontsize',28, 'Rotation',90)
space='\hspace{6pt}'; % spacing between title elements
title(['$K''','=',num2str(K),';',space,'\overline{\delta ''}=',num2str(deltap), ...
    ';', space, '\overline{\kappa}=',num2str(kappa),';',space,droptype,'$'],...
    'interpreter','latex','fontsize',16)
axis([-1.1*rrange,1.1*rrange,-1.1*zrange,0.1*zrange])
rsize=2.2*rrange;
zsize=1.2*zrange;
set(gca,'units','points')
set(gcf,'units','points')
if 2*rsize>zsize
    raxis=plotsize;
    zaxis=plotsize*zsize/rsize;
else
    zaxis=plotsize;
    raxis=plotsize*rsize/zsize;
end
set(gcf,'Position',[10, 10, raxis+100, zaxis+120]) 
set(gca,'Position',[60, 50, raxis, zaxis])
legend(fig,'$\lambda_{0}=1$',['$\lambda_{0}=\hspace{3pt}$',num2str(l0)])
h = legend;
set(h, 'interpreter', 'latex') % enable interpretation of Latex symbols in legends
set(h, 'location', 'best')
figname=['data/K',num2str(K),'d',num2str(deltap),'k',num2str(kappa),...
    'l',num2str(l0)];
xmark=[-1 1];
ymark=[0 0];
plot(xmark,ymark,'o','markerfacecolor','blue')
text(0.5,0.9,'(a)','Units','normalized','fontsize',24,'HorizontalAlignment','center')
saveas(gcf,'fig13a.fig');
% now plot tensions
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
hold off
plot(s,txi,'--b')
hold on
plot(s,tphi,'-r')
xlabel('$\sigma''$','fontsize',24,'interpreter','latex')
ylabel(['$\overline{T_{\xi}},',space,'\overline{T_{\phi}}$'], 'interpreter', 'latex','rotation',90,'fontsize',20);
% title(['$K''','=',num2str(K),';',space,'\overline{\delta ''}=',num2str(deltap), ...
%   ';', space, '\overline{\kappa}=',num2str(kappa),';', space, '\lambda_{0}=',num2str(l0),...
%   ';',space,droptype,'$'], 'interpreter','latex','fontsize',16)
set(gcf,'Position',[10, 10, plotsize+120, plotsize+120])
set(gca,'Position',[80, 50, plotsize, plotsize])
legend('$\overline{T_{\xi}}$', '$\overline{T_{\phi}}$')
h = legend;
set(h, 'interpreter', 'latex') % enable interpretation of Latex symbols in legends
set(h, 'location', 'northwest')
text(0.5,0.9,'(b)','Units','normalized','fontsize',24,'HorizontalAlignment','center')
saveas(gcf,'fig13b.fig');
hold off
% Now calculate some of the characteristic gemoetric parameters from the solution
Vp=trapz(zp,pi.*rp.^2);   % this is the normzlized volume of the undeformed drop
Ap=lp*(trapz(s,2*pi.*rp));   % this is the normalized area of the deformed drop
V=trapz(z,pi.*r.^2);   % this is the normzlized volume of the undeformed drop
A=lp*(trapz(s,2*pi.*r));   % this is the normalized area of the deformed drop