function project()

%-------------------------------------------------------------------------
% ME 574 Course Project - Winter 2010
% Coded by Sukhoon Pyo, Junghwan Seo, Win Khwa, Nien-Tsu Huang
% Simulate The Evaporation and Rehydration Process of Aqueous Two Phase System (ATPS)
%-------------------------------------------------------------------------
clear all; close all; clc;

%----------------------------------------------------------------------
% Evaporation Process
%----------------------------------------------------------------------

% Set Simulation parameters
gamma     = 2.4;      % set coeficient of isotropic surface energy
dgamma    = 0;
g1        = 0.9;        % set coeficients of volume energy
g2        = 0.2;
m         = 1;        % set element mass
timesteps = 1800;     % set number of timesteps to calculate
dt        = 0.001;    % set timestep length

% Set dimensions of first droplet
TTThe=95;   %input initial contact angle  
n=15;        %input #of nodes of upper side
LL=3.02;       %bottom line
TTH=TTThe*pi/180;
RR=LL/2/sin(TTH);
PHII=2.*TTH/n;

% Set figure size
xmin = -LL/2 -1;               % set Animation xmin
xmax =  LL/2+1;              % set Animation xmax
ymin = -1;              % set Animation ymin
ymax = 3;              % set Animation ymax


% Movie Settings
filename    = 'test1'; % first filename to save as
filename2   = 'test2'; % second filename to save as
fps         = 40;               % set animation frames per second
compression = 'Cinepak';        % set 'none' on OsX, Linux & Unix, 'Indeo3' for WinXP, 'Cinepak' for Vista 
AVIquality  = 70;               % set animation quality

% Set initial node coordination
nodeCoordinates=[];
nodeCoordinates(1,:)=[-LL/2 0]; 
for i=1:n
    x1=RR*(-cos((2.*TTH/n)*i+0.5*pi-TTH)+cos(0.5*pi-TTH))-LL/2;
    y1=RR*(sin((2.*TTH/n)*i+0.5*pi-TTH)-sin(0.5*pi-TTH));
    nodeCoordinates(i+1,:)=[x1 y1]; 
end

x=nodeCoordinates(:,1);
y=nodeCoordinates(:,2);


% Solving iteration
for t = 1:timesteps;
        
        % store current frame for animation
        Frames(t) = plotgraph(x,y,xmin,xmax,ymin,ymax);
        %evolve nodes to next timestep
        [ x, y] = SOLVE1(x,y,n,gamma,dgamma,g1,g2,m,dt,t);
        CANG(t) = ContactAngle(x,y);
        vol(t) = Volume(x,y,n);
        SArea(t) = SurfaceArea(x,y,n);
        time(t) = t*dt;
       
end % end of the iteration 1

% store final frame for animation
%movie2avi(Frames,filename,'fps',fps,'compression',compression,'quality',AVIquality);

% Plot Contact angle corresponds to Time
figure(2)
plot(time,CANG);
xlabel('Time');
ylabel('Contact Angle');
hold;

% Plot Volume corresponds to Time
figure(3)
plot(time,vol);
xlabel('Time');
ylabel('Volume');
hold;


% Plot Surface Area corresponds to Time
figure(4)
plot(time,SArea);
xlabel('Time');
ylabel('Surface Area');
hold;

%----------------------------------------------------------------------
% Rehydration Process
%----------------------------------------------------------------------
% Set Simulation parameters

gamma     = 0.1;      % set coeficient of isotropic surface energy
dgamma    = 0;
g1        = 1;        % set coeficients of volume energy
g2        = 0.2;
m         = 1;        % set element mass
TimeSteps = 1300;     % set number of timesteps to calculate
dt        = 0.001;    % set timestep length

GAMM     = 2.4;      % set coeficient of isotropic surface energy
DGAMM    = 0;
GG1        = 0.9;        % set coeficients of volume energy
GG2        = 0.2;

% Set dimensions of second droplet
TTTheta=95;   %input initial contact angle
nn=14;        %input #of nodes of upper side
LLL=5;       %bottom line
TTHH=TTTheta*pi/180;
RRR=LLL/2/sin(TTHH);
PPHII=2.*TTHH/nn;

% Set figure size
xmin = -LLL/2 -1;               % set Animation xmin
xmax =  LLL/2+1;              % set Animation xmax
ymin = -1;              % set Animation ymin
ymax = 4;              % set Animation ymax


% Set initial node coordination
nodeCoord=[];
nodeCoord(1,:)=[-LLL/2 0]; 
for i=1:nn
    x1=RRR*(-cos((2.*TTHH/nn)*i+0.5*pi-TTHH)+cos(0.5*pi-TTHH))-LLL/2;
    y1=RRR*(sin((2.*TTHH/nn)*i+0.5*pi-TTHH)-sin(0.5*pi-TTHH));
    nodeCoord(i+1,:)=[x1 y1]; 
end

for i=1:nn+1
    j=i+n+1;
    x(j,1)=nodeCoord(i,1);
    y(j,1)=nodeCoord(i,2);
end

% Solving iteration
for t = 1:TimeSteps;
        
        % store current frame for animation
        Frames(t) = plotgraph(x,y,xmin,xmax,ymin,ymax);
        %evolve nodes to next timestep
        [ x, y] = SOLVE2(x,y,n,gamma,dgamma,g1,g2,m,dt,t);
        [ x, y] = SOLVE3(x,y,n,nn,GAMM,DGAMM,GG1,GG2,m,dt,t);
        CANG2(t) = ContactAngle(x,y);
        vol2(t) = Volume(x,y,nn);
        SArea2(t) = SurfaceArea(x,y,nn);
        Time(t) = t*dt+timesteps*dt;
       
end 

% Plot Contact angle corresponds to Time
figure(2)
plot(Time,CANG2);
xlabel('Time');
ylabel('Contact Angle');

% Plot Volume corresponds to Time
figure(3)
plot(Time,vol2);
xlabel('Time');
ylabel('Volume');

% Plot Surface Area corresponds to Time
figure(4)
plot(Time,SArea2);
xlabel('Time');
ylabel('Surface Area');


% store final frame for animation
%movie2avi(Frames,filename2,'fps',fps,'compression',compression,'quality',AVIquality);

end % end of iteration 2

function [Frame] = plotgraph(X1val, Y1val, xmin, xmax, ymin, ymax)
% Function grabplot saves frame of current state for animation
x0 = xmin:xmax;
y0 = zeros(1,length(x0));
figure(1);
plot(X1val,Y1val,'o-',x0,y0);
axis([xmin, xmax, ymin, ymax]);
Frame = getframe;
end


function [ x y ] = SOLVE1( x, y, n, gamma, dgamma, g1, g2, m, dt,timesteps)
%Function to formulate H&F matricies and updates x and y values to next
%timestep in structural evolution


% Create the global H and F matrices

% zero setting
H = zeros(2*n+2,2*n+2);
Hl = zeros(4);
f = zeros(2*n+2,1);
fl = zeros(4,1);
theta = zeros(n,1);
l = zeros(n,1);

%define elements
%calculate lenght and angle of nodes
for i = 1:n
    %calculate length
    l(i) = ((x(i+1)-x(i))^2+(y(i+1)-y(i))^2)^(1/2);
    %calculate angle
    theta(i) = asin((y(i+1)-y(i))/l(i));
    if x(i+1) < x(i)
        theta(i)=pi-theta(i);
    end
end

for j=1:n;

    %Calculate local H matrix elements
    Hl(1,1)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(1,2)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(1,3)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(1,4)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,1)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(2,2)=l(j)/(6*m)*2*(cos(theta(j)))^2;
    Hl(2,3)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,4)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(3,1)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(3,2)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(3,3)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(3,4)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,1)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(4,2)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(4,3)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,4)=l(j)/(6*m)*2*(cos(theta(j)))^2;

    %add local h matrix to global h matrix
    H(2*j-1:2*j+2,2*j-1:2*j+2) = H(2*j-1:2*j+2,2*j-1:2*j+2) + Hl;


    % calculate Force Vector
    fl(1)=gamma*cos(theta(j))+dgamma*(-sin(theta(j)))+l(j)*(g1-g2)/2*(-sin(theta(j)));
    fl(2)=gamma*sin(theta(j))+dgamma*(cos(theta(j)))+l(j)*(g1-g2)/2*(cos(theta(j)));
    fl(3)=gamma*-cos(theta(j))+dgamma*(sin(theta(j)))+l(j)*(g1-g2)/2*(-sin(theta(j)));
    fl(4)=gamma*-sin(theta(j))+dgamma*(-cos(theta(j)))+l(j)*(g1-g2)/2*(cos(theta(j)));

    %add local force vector to global force vector
    f(2*j-1:2*j+2)=f(2*j-1:2*j+2)+fl;


end

%Add Artificial Finite Mobility (small number to diagonal points in H matrix)
for m=1:(2*n+2);
    H(m,m)=0.000001+H(m,m);
end

%boundary conditions

H(end,end) = (H(2*n+2,2*n+2))*1E10;  %fix last point in y direction
H(2,2) = (1+H(2,2))*1E10;                  %fix first point in y direction
%H(end-1,end-1) = (H(2*n+1,2*n+1))*1E10;  %fix last point in x direction
%H(1,1) = (1+H(1,1))*1E10;                  %fix first point in x direction

% frictional forces (Viscous forces)
if (timesteps<700)
    FRF=timesteps*0.002;
elseif (timesteps>=700 & timesteps<1300)
    FRF=1.4+(timesteps-700)*0.001;
else
    FRF=2+(timesteps-1300)*0.0003;
end


f(1)=f(1)-0.988*FRF;
f(3)=f(3)-FRF/1500;
f(4)=f(4)+FRF/500;
f(5)=f(5)-FRF/3000;
f(6)=f(6)-FRF/3000;
f(end-1)=f(end-1)+0.988*FRF;
f(end-2)=f(end-2)+FRF/500;
f(end-3)=f(end-3)+FRF/1500;
f(end-5)=f(end-5)+FRF/3000;
f(end-4)=f(end-4)-FRF/3000;


%solve for velocities
vel=H\f;

%update positions using Euler method
for i=1:(n+1);
    x(i)=x(i)+vel(2*i-1)*dt;
    y(i)=y(i)+vel(2*i)*dt;
end
end % end of SOLVE1 function


function [ x y ] = SOLVE2( x, y, n, gamma, dgamma, g1, g2, m, dt,timesteps)
%Function to formulate H&F matricies and updates x and y values to next
%timestep in structural evolution


% Create the global H and F matrices

% zero setting
H = zeros(2*n+2,2*n+2);
Hl = zeros(4);
f = zeros(2*n+2,1);
fl = zeros(4,1);
theta = zeros(n,1);
l = zeros(n,1);

%define elements
%calculate lenght and angle of nodes
for i = 1:n
    %calculate length
    l(i) = ((x(i+1)-x(i))^2+(y(i+1)-y(i))^2)^(1/2);
    %calculate angle
    theta(i) = asin((y(i+1)-y(i))/l(i));
    if x(i+1) < x(i)
        theta(i)=pi-theta(i);
    end
end

for j=1:n;

    %Calculate local H matrix elements
    Hl(1,1)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(1,2)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(1,3)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(1,4)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,1)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(2,2)=l(j)/(6*m)*2*(cos(theta(j)))^2;
    Hl(2,3)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,4)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(3,1)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(3,2)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(3,3)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(3,4)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,1)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(4,2)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(4,3)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,4)=l(j)/(6*m)*2*(cos(theta(j)))^2;

    %add local h matrix to global h matrix
    H(2*j-1:2*j+2,2*j-1:2*j+2) = H(2*j-1:2*j+2,2*j-1:2*j+2) + Hl;


    % calculate Force Vector
    fl(1)=gamma*cos(theta(j))+dgamma*(-sin(theta(j)))+l(j)*(g1-g2)/2*(-sin(theta(j)));
    fl(2)=gamma*sin(theta(j))+dgamma*(cos(theta(j)))+l(j)*(g1-g2)/2*(cos(theta(j)));
    fl(3)=gamma*-cos(theta(j))+dgamma*(sin(theta(j)))+l(j)*(g1-g2)/2*(-sin(theta(j)));
    fl(4)=gamma*-sin(theta(j))+dgamma*(-cos(theta(j)))+l(j)*(g1-g2)/2*(cos(theta(j)));

    %add local force vector to global force vector
    f(2*j-1:2*j+2)=f(2*j-1:2*j+2)+fl;


end

%Add Artificial Finite Mobility (small number to diagonal points in H matrix)
for m=1:(2*n+2);
    H(m,m)=0.000001+H(m,m);
end

%boundary conditions
%fix end points comment and uncomment lines to turn off/on boundary
%conditions

H(end,end) = (H(2*n+2,2*n+2))*1E10;  %fix last point in y direction
H(2,2) = (1+H(2,2))*1E10;                  %fix first point in y direction
%H(end-1,end-1) = (H(2*n+1,2*n+1))*1E10;  %fix last point in y direction
%H(1,1) = (1+H(1,1))*1E10;                  %fix first point in y direction

% friction force
%if (timesteps<1100)
%    FRF=timesteps*0.0018;
%elseif (timesteps>=1100 & timesteps<1200)
%    FRF=1.9;
%else
%    FRF=1.98+timesteps*0.00008;
%end

%f(1)=f(1)-FRF;
%f(end-1)=f(end-1)+FRF;



%solve for velocities
vel=H\f;

%update positions using Euler method
for i=1:(n+1);
    x(i)=x(i)+vel(2*i-1)*dt;
    y(i)=y(i)+vel(2*i)*dt;
end
end % end of SOVLE2 function




function [ a b ] = SOLVE3(X, Y, n, nn, GAMM, DGAMM, GG1, GG2, m, dt, timesteps)
%Function to formulate H&F matricies and updates x and y values to next
%timestep in structural evolution

a = X;
b = Y;
for i = 1:(nn+1);
    x(i)=X(i+n+1);
    y(i)=Y(i+n+1);
end

% Create the global H and F matrices

% zero setting
H = zeros(2*nn+2,2*nn+2);
Hl = zeros(4);
f = zeros(2*nn+2,1);
fl = zeros(4,1);
theta = zeros(nn,1);
l = zeros(nn,1);

%define elements
%calculate lenght and angle of nodes
for i = 1:nn
    %calculate length
    l(i) = ((x(i+1)-x(i))^2+(y(i+1)-y(i))^2)^(1/2);
    %calculate angle
    theta(i) = asin((y(i+1)-y(i))/l(i));
    if x(i+1) < x(i)
        theta(i)=pi-theta(i);
    end
end

for j=1:nn;

    %Calculate local H matrix elements
    Hl(1,1)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(1,2)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(1,3)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(1,4)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,1)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(2,2)=l(j)/(6*m)*2*(cos(theta(j)))^2;
    Hl(2,3)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(2,4)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(3,1)=l(j)/(6*m)*(sin(theta(j)))^2;
    Hl(3,2)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(3,3)=l(j)/(6*m)*2*(sin(theta(j)))^2;
    Hl(3,4)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,1)=l(j)/(6*m)*-sin(theta(j))*cos(theta(j));
    Hl(4,2)=l(j)/(6*m)*(cos(theta(j)))^2;
    Hl(4,3)=l(j)/(6*m)*-2*sin(theta(j))*cos(theta(j));
    Hl(4,4)=l(j)/(6*m)*2*(cos(theta(j)))^2;

    %add local h matrix to global h matrix
    H(2*j-1:2*j+2,2*j-1:2*j+2) = H(2*j-1:2*j+2,2*j-1:2*j+2) + Hl;


    % calculate Force Vector
    fl(1)=GAMM*cos(theta(j))+DGAMM*(-sin(theta(j)))+l(j)*(GG1-GG2)/2*(-sin(theta(j)));
    fl(2)=GAMM*sin(theta(j))+DGAMM*(cos(theta(j)))+l(j)*(GG1-GG2)/2*(cos(theta(j)));
    fl(3)=GAMM*-cos(theta(j))+DGAMM*(sin(theta(j)))+l(j)*(GG1-GG2)/2*(-sin(theta(j)));
    fl(4)=GAMM*-sin(theta(j))+DGAMM*(-cos(theta(j)))+l(j)*(GG1-GG2)/2*(cos(theta(j)));

    %add local force vector to global force vector
    f(2*j-1:2*j+2)=f(2*j-1:2*j+2)+fl;


end

%Add Artificial Finite Mobility (small number to diagonal points in H matrix)
for m=1:(2*nn+2);
    H(m,m)=0.000001+H(m,m);
end

%boundary conditions

H(end,end) = (H(2*nn+2,2*nn+2))*1E10;  %fix last point in y direction
H(2,2) = (1+H(2,2))*1E10;                  %fix first point in y direction
%H(end-1,end-1) = (H(2*n+1,2*n+1))*1E10;  %fix last point in y direction
%H(1,1) = (1+H(1,1))*1E10;                  %fix first point in y direction

% friction force
if (timesteps<300)
    FRF=0;
else
    FRF=timesteps*0.000005;
end
%elseif (timesteps>=1100 & timesteps<1200)
%    FRF=1.9;
%else
%    FRF=1.98+timesteps*0.00008;
%end

f(1)=f(1)-0.988*FRF;
f(3)=f(3)-FRF/1500;
f(4)=f(4)+FRF/500;
f(5)=f(5)-FRF/3000;
f(6)=f(6)-FRF/3000;
f(end-1)=f(end-1)+0.988*FRF;
f(end-2)=f(end-2)+FRF/500;
f(end-3)=f(end-3)+FRF/1500;
f(end-5)=f(end-5)+FRF/3000;
f(end-4)=f(end-4)-FRF/3000;


%solve for velocities
vel=H\f;

%update positions using Euler method
for i=1:(nn+1);
    x(i)=x(i)+vel(2*i-1)*dt;
    y(i)=y(i)+vel(2*i)*dt;
end

for i = 1:(nn+1);
    a(i+n+1)=x(i);
    b(i+n+1)=y(i);
end

end % end of SOVLE3 function

%calculate Contact angle
function [TH] = ContactAngle(x,y)
TH = atan((y(2)-y(1))/(x(2)-x(1)))*180/pi;
end


%calculate surface area
function [sa] = SurfaceArea(x,y,n)
sa = 0;
for i = 1: n
    sa = sa + ((x(i+1)-x(i))^2+(y(i+1)-y(i))^2)^(1/2);
end
end

%calculate volume
function [v] = Volume(x,y,n)
v = 0;
for i = 1: n
    v = v + (x(i+1)-x(i))*(y(i+1)+y(i))/2;
end
end
