%---
 close all
 clear all
%---

%=========================
% meniscus on an inclined plate
%=========================

%---
% data
%---

gac = 1.0;    % acceleration of gravity
rhop = 1.0;    % density of the pool
rhoa = 0.0;   % density of the ambient fluid
gamma = 1.0;  % surface tension

beta = 0.01*pi;   % plate inclination angle

alpha = 0.99*pi;  % contact angle
alpha = 0.01*pi;
alpha = 0.90*pi;

ndiv = 2*64;

%---
% prepare
%---

 drho = rhop-rhoa;

 if(drho<0)
  disp 'The density of the pool must be higher than'
  disp 'the density of the ambient fluid'
  return
 end

 capls = gamma/(gac*abs(drho));
 capl  = sqrt(capls);

 dbeta = 0.05;
 Irepeat = 1;

%========
while (Irepeat==1)  % loop over inclination angles
%========

%--------------------------------
% angle theta at the contact line
%--------------------------------

  thcl = alpha+beta;
  cst  = cos(thcl);

%-----------------------------------------
% compute the magnitude of meniscus rise h
%-----------------------------------------

  if((thcl>0.5*pi)&(thcl<1.5*pi))
   h = sqrt(2.0*capls*(1.0-abs(cst)));
  else
   h = sqrt(2.0*capls*(1.0+abs(cst)));
  end

%---
% sign of meniscus rise
%---

  if(thcl<pi)
     h = abs(h);   % meniscus goes up
  else
     h =-abs(h);   % meniscus goes down
  end

%---
% plate position at y=0
%---

  xpl = -h/tan(beta);

%------------------------------------------
% integrate meniscus equation dx/df = G(f)
% from f=h to f=0
% using the modified Euler method with constant step
%------------------------------------------

    df  = h/ndiv;  % note that df can be positive or negative
    dfh = 0.5*df;

    y(1) = h;
    x(1) = 0.0;

     for i=1:ndiv-1
       fred  = y(i)/capl;
       freds = fred*fred;
       xp    = (2.0-freds)/sqrt(4.0-freds)/fred;
       xsv  = x(i);         % save
       xpsv = xp;
       y(i+1) = y(i)-df;
       x(i+1) = x(i)+xp*df;
       fred  = y(i+1)/capl;
       freds = fred*fred;
       xp    = (2.0-freds)/(fred*sqrt(4.0-freds));
       x(i+1) = xsv+(xpsv+xp)*dfh;
     end

%---
% plotting
%---

 plot(x,y)
 hold on
 patch([xpl, 0, x, x(ndiv), (-2-h)/tan(beta)],[0, h, y,-2, -2],'y');
 plot([-1 3 3 -1 -1],[-2 -2 2 2 -2]);
 plot([-10 10],[-10*tan(beta)+h, 10*tan(beta)+h],'r','LineWidth',3);
 plot([xpl, 10],[0, 0],'c--','LineWidth',1);

 hold off

 axis equal
 set(gca,'fontsize',15)
 xlabel('x','fontsize',15)
 ylabel('y','fontsize',15)
 axis([-1 3 -2 2])
 pause(0.1)

%---
% tilt the plate
%---

 beta = beta+dbeta;

 if(beta>0.99*pi|beta<0.01*pi)
  dbeta=-dbeta;
 end

%====
end
%====