
% a code to compute the "2D drop" with contact angle, surface tension, and weight
% using the Newton iterations
%
% here the pressure increase between atmosphere and bottom is taken as a variable
% and the left contact angle is imposed
%
% there is a loop on L, the drop length, to see the change from drop to "flaque"...

clear all; clf

g=9.81 % gravity
ro=1000 % density
sig=72e-3 % surface tension
N=200 % number of points
h=0    % left height
ang=pi/5 % left contact angle

testvec=linspace(0.003,0.05,5)
%testvec=0.05
for ind=1:length(testvec) % loop on the drop length
    L=testvec(ind);
    
    %%%% finite differences
    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);
    
    %%%% state constraints
    bou=[1;N]; 
    C=[[I(bou,:); D(1,:)],zeros(3,1)]; % constraints
    loc=[bou;N+1]; % location where to put the constraints
    val=[h;0;tan(ang)]; % values of the constraints
    
    %%%% initial guess
    F=(-(tan(ang)*L+h)/L^2)*x.^2+tan(ang)*x+h;
    FF=[F;0];
    
    quit=0; count=1;
    while ~quit % Loop of Newton
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        F=FF(1:N); dp=FF(end);
        
        %%%% 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)-dp;0];
        q(loc)=C*FF-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
        J=ro*g*I-diag(3*sig*(Fxx.*Fx./(1-Fx.^2).^2.5))*D-diag(sig./(1-Fx.^2).^1.5)*DD;
        A=[J,-ones(N,1); zeros(1,N+1)]; 
        A(loc,:)=C; % state constraints
        
        %%%% iterate Newton
        FF=FF-A\q;
        count=count+1;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    plot(x,F(1:N)); axis equal; hold on; 
    drawnow
end

