function drupHiBo(inputfile)                  % S.W. Rienstra 8-11-2022, based on "The shape of a sessile drop for small and large surface tension", J. Engin. Math. 1990, 24, 193-202
   %clear all;
   opengl('save', 'software')
   warning off MATLAB:MKDIR:DirectoryExists
   format long
   format compact
   %% start input parameters
   if nargin == 1
      %% in case the input is provided by an inputfile
      display(['inputfile = ', inputfile]);
      fid = fopen(inputfile,'r');
      eval(fscanf(fid,'%c'));
      fclose(fid);
   else
    %% input parameters
    Bond     = [4,8,16,32,64,128,256,512,1024];             % Bond number. Results acceptable for Bo>=20
   %Bond     = 64;
   %Bond     = [20,40,80,160,320,640,1280];
    alph     = 165*pi/180;           % contact angle
    %
    L        = 1;                    %   L in mm
    Nt       =  1000;                % number of t-steps minus 1
   end
   %% end input parameters
   tstart = tic;
   %%%%%%%%%%%%%%
   eps  = 1./sqrt(Bond);
   lbd  = sqrt(2)-1;
   tnu  = 1 + eps.*(2-sqrt(2) + eps*(1-2*log(1+lbd^2))/3);
   bt   = 4*lbd*sqrt(2*pi*tnu)./eps./sqrt(eps).*exp(-1./eps-2+sqrt(2)+eps*(sqrt(2)/6-5/24))
   for j=1:length(Bond)
     tend = fsolve(@(t) psit(t,eps(j),tnu(j))-alph,1.3);
     t    = tend*[0:Nt]/Nt;
     psi(j,:)  = psit(t,eps(j),tnu(j))*180/pi;
     xsi(j,:)  = xsit(t,eps(j),tnu(j));
     eta(j,:)  = etat(t,eps(j),tnu(j),bt(j));
     pa(j)   = psi(end);
     xa(j)   = xsi(end);
     ea(j)   = eta(end);
   end
   V       = 2*pi./Bond.*xa/2.*((Bond.*ea+bt).*xa - 2*sin(pa));
   Va      = 2*pi.*eps*sin(alph/2);                                  % large Bond-limit
   if length(Bond) == 1
    fprintf(' s_a   = %3.4f mm     \n',tend*L)                                                % s_a   in mm
    fprintf(' psi_a = %3.4f degrees\n',pa       )                                             % psi_a in degrees
    fprintf(' r_a   = %1.6f mm     \n',xa*L     )                                             % r_a in mm
    fprintf(' z_a   = %1.6f mm     \n',ea*L     )                                             % z_a in mm
    fprintf(' V     = %1.6f mm^3   \n',V*L^3    )                                             % V in mm^3
    fprintf(' V (ap)= %1.6f mm^3   \n',Va*L^3   )                                             % V in mm^3
   end
   %% make plots of psi, (r,z)-shape, log(beta), and volume
   figure(1)
   clf
   for j=1:length(Bond)
    plot(t,psi(j,:),'-k','LineWidth',1)
    hold on
   end
   xlabel('t')
   ylabel('\psi')
   hYLabel = get(gca,'YLabel');
   set(hYLabel,'rotation',0,'VerticalAlignment','middle')
   hold off
   figure(2)
   clf
   for j=1:length(Bond)
     plot(xsi(j,:),ea(j)-eta(j,:),'-k','LineWidth',1)
     hold on
     plot([0,xa(j)],[0,0],'-k','LineWidth',1)
   end
   hold off
   axis on
   axis equal
   axis([0,1.02,0,1.10*max(ea)])
   xlabel('r')
   ylabel('-z')
   hYLabel = get(gca,'YLabel');
   set(hYLabel,'rotation',0,'VerticalAlignment','middle')
   figure(3)
   clf
   plot(eps,log(bt),'-k','LineWidth',1)
   xlabel('\epsilon')
   ylabel('log(\beta)')
   hYLabel = get(gca,'YLabel');
   set(hYLabel,'rotation',0,'VerticalAlignment','middle')
   figure(4)
   clf
   plot(eps,V,'-k','LineWidth',1)
   xlabel('\epsilon')
   ylabel('V')
   hYLabel = get(gca,'YLabel');
   set(hYLabel,'rotation',0,'VerticalAlignment','middle')
   hold on
   plot(eps,Va,'-','LineWidth',1,'Color','#EDB120')
   hold off
   %
   fprintf(' time  = %2.4f sec\n',toc(tstart) );
   %%  functions to define psi, xi, eta
   function p = psit(t,ep,tu)
       tau  = (t-tu)/ep;
       Q    = 1/12*lbd^2*(exp(2*tau)-1)-tau/2+2/3*log((1+lbd^2*exp(2*tau))/(1+lbd^2));
       p0   = 4*atan(lbd*exp(tau));
       p1   = Q.*(4*lbd*exp(tau)./(1+lbd^2*exp(2*tau)));
       p    = p0 + ep*p1; % Note the spill-over error in the composite expansion from the boundary layer at t=tnu to the outer solution at t=0
   end
   function x = xsit(t,ep,tu);
       tau  = (t-tu)/ep;
       Q    = 1/12*lbd^2*(exp(2*tau)-1)-tau/2+2/3*log((1+lbd^2*exp(2*tau))/(1+lbd^2));
       x0   = tau+2-sqrt(2)-4*lbd^2*exp(2*tau)./(1+lbd^2*exp(2*tau));
       x1   = 1/3+2/3*log((1+lbd^2*exp(2*tau))/(1+lbd^2))-1/3*exp(2*tau).*((1+lbd^2)./(1+lbd^2*exp(2*tau))).^2;
       x1   = x1-1/2*Q.*((4*lbd*exp(tau)./(1+lbd^2*exp(2*tau))).^2);
       xs   = 1 + ep*(x0 + ep*x1);
       x    = (xs-xs(1))/(1-xs(1));
   end
   function y = etat(t,ep,tu,b);
       tau  = (t-tu)/ep;
       Q    = 1/12*lbd^2*(exp(2*tau)-1)-tau/2+2/3*log((1+lbd^2*exp(2*tau))/(1+lbd^2));
       et0  = 4*lbd*exp(tau)./(1+lbd^2*exp(2*tau));
       et1  = ((1-lbd^2*exp(2*tau)).*Q+lbd^4*exp(4*tau)/6+1/2).*et0./(1+lbd^2*exp(2*tau));
       et   = ep^2*b*(et0 +ep*et1 ).*besseli(0,t/ep)./(4*lbd*exp((1-ep/2)*tau+ep*(1/2-lbd^2/12-2/3*log(1+lbd^2))));
       y    = ((ep^2*b/et(1)).^(1-t/tu)).*et - ep^2*b;
   end
end  % end main program
