
% a code to compute the nonlinear meniscus
% using Newton iteration
% with analytic jacobian

clear all; clf

g=9.81 % gravity
ro=1000 % density
sig=72e-3 % surface tension
L=0.05 % domain length
N=200 % number of grid points
ang=pi/3 % contact angle

%%%% finite differences matrices
x=linspace(0,L,N)'; dx=x(2)-x(1);
D=zeros(N,N);
D(1,1:3)=[-1.5,2,-0.5]/dx; 
for ind=2:N-1;    
    D(ind,ind-1:ind+1)=[-1,0,1]/(2*dx); 
end
D(N,N-2:N)=[0.5,-2,1.5]/dx;
DD=D^2; % second order derivative
I=eye(N);

testvec=linspace(pi/2,pi/3.5,5);
for ind=1:length(testvec);
    ang=testvec(ind)
    %%%% state constraints
loc=[1;N]; % place to substitute
C=[D(1,:); I(N,:)];   % constraints
val=[-1/tan(ang);0]; % values of the constraints

%%%% initial guess
F=0*x;

quit=0; count=1;
while ~quit % Loop of Newton
    
    %%%% derivatives of present state
    Fx=D*F;
    Fxx=DD*F;
    
    %%%% compute nonlinear function for present iterate
    q=ro*g*F-sig*Fxx./((1-Fx.^2).^1.5);
    q(loc)=C*F-val; % state constraints
    
    %%%% compute residual, plot and test for convergence
    res=max(abs(q)); 
   
    disp(['it: ' num2str(count) ' res: ' num2str(res)]);
    if res<1e-8; disp('converged'); quit=1; continue; 
    elseif res>1e4|count>20; disp('no convergence...'); break; 
    end ; 
    
    %%%% build jacobian
    A=ro*g*I-diag(3*sig*(Fxx.*Fx./(1-Fx.^2).^2.5))*D-diag(sig./(1-Fx.^2).^1.5)*DD;
    A(loc,:)=C; % state constraints
    
    %%%% iterate Newton
    F=F-A\q;
    
    count=count+1;
end
plot(x,F,'b');  axis equal; drawnow; hold on
end


