%Jacob Dresher
%jd352@njit.edu
%New Jersey Institute of Technology
%Simple Thin Film Fluid
%Open Source Code
clear;clc;close all;

%intalize our grid and time
a = 0; b = 1;
x = linspace(a,b,100);
t = linspace(0,10,100);

%initalize our units
%[p g mew sigma lamda1 theta alpha kapa delsig ~ dx dt]
coe = [1 1 1 1 1 0 1 1 1];tic
[~,decay] = Navier(x,t,b,coe);toc

fprintf('Theoretical Decay: %f\nExperimental Decay: %f\nReletive Error: %f\n',decay(1),decay(2),decay(3));

%--------------------------------------------------------------------------------------------------------------
function [ h,decay ] = Navier( x,t,b,coe )
%solves the fluid height over time with periodic boundary conditions

%save the coeiffents we will use
eps = .001;
h0 = .1;

%set up our grid and step sizes
dx = x(3) - x(2);
dt = t(2) - t(1);

%[p g mew sigma lamda1 theta alpha kapa delsig ~ dx dt]
coe = [coe dx dt]; 

%our intial conditioneps
h = h0*ones(length(x),1) + eps*sin(((2*pi)/b)*x)';
indexHight = find(h==max(h),1);

%ploting set up 
figure; hold on;
n = floor(length(t)/10);
saveHight = zeros(1,length(t));

%implcit Euler solver
for i=1 : length(t)
   saveHight(i) = h(indexHight);
   if(n==floor(length(t)/10))
      plot(x,h'); 
      n = 1;
   else
      n = n + 1; 
   end
   %root solver
   h = NewtonRoots(h,coe);
end
%plot last time step
plot(x,h'); hold off;
axis([0 x(end) min(h)-(10*eps) max(h)+(10*eps)]);
xlabel('X');
ylabel('Hight');
title('Fluid Hight vs Time');

%finds numerical omega
omega = real(log(saveHight - h0) - log(eps*h0));
%cases when to take the slope depending on lamda
if (coe(6) == 0)
    omega =  ((omega(end) - omega(1))/(t(end) - t(1)));
else
    omega =  ((omega(2) - omega(1))/(t(2) - t(1)));
end

%finds the theoretical decay
p = coe(1);
gr = coe(2);
mew = coe(3);
sigma = coe(4);
lam1 = coe(5);
theta = coe(6);
alpha = coe(7);
kapa = coe(8);
delsig = coe(9);

%get our theoretical omega
w = @(k) (-(h0^3*sigma/3)*k.^4 - (cos(theta)*p*gr/3)*h0^3*k.^2 + (h0^2)*(alpha*delsig/(2*kapa))*k.^2)/(mew + (k.^4)*(sigma*lam1/3)*h0^3);
decayT = w(2*pi/b);
decayE = omega;
%reletive Error
RE = abs(decayT - decayE)/abs(decayT);

decay = [decayT decayE RE];
end

function [ x0 ] = NewtonRoots( h,coe )
%newton root finding method

%inialize our guess and tolerance
x0 = h;
eps = .0000001;
MAXITTER = 10000;

%itterative root finding method
for i=1 : MAXITTER
   Jacob = myJacobian(x0,h,coe);
   F = functionVector(x0,h,coe);
   xP = x0 - Jacob\F;
   if (norm(xP - x0) < eps)
      break; 
   end
   x0 = xP;
end
x0 = xP;
end

function [ Jacob ] = myJacobian( g,hP,coe )
%Generates the Jacobian for the root finding step

%initalize our coeiffents
p = coe(1);
gr = coe(2);
mew = coe(3);
sigma = coe(4);
lam1 = coe(5);
theta = coe(6);
alpha = coe(7);
kapa = coe(8);
delsig = coe(9);
dx = coe(10);
dt = coe(11);

%initalize function diffs
fm2 = @(hm2,hm1,h,hp1,hp2,hIm2,hIm1,hI,hIp1,hIp2) -(lam1*sigma*(h + hm1)^2*(h - hI - hIm1 + hm1))/(16*dt*dx^4)...
    + (sigma*(h + hm1)^3)/(24*dx^4)...
    + (lam1*sigma*(h + hm1)^3)/(24*dt*dx^4);

fm1 = @(hm2,hm1,h,hp1,hp2,hIm2,hIm1,hI,hIp1,hIp2) -(lam1*sigma*((h + hm1)^2*(3*h - 3*hm1 + hm2 - hp1) - 3*(h + hm1)^2*(h - hI - hIm1 + hm1) - (h + hp1)^2*(h - hI - hIp1 + hp1) + (2*h + 2*hm1)*(3*h - 3*hm1 + hm2 - hp1)*(h - hI - hIm1 + hm1)))/(16*dt*dx^4)...
    + -(gr*p*(3*(h + hm1)^2*(sin(theta) - (cos(theta)*(h - hm1))/dx) + (cos(theta)*(h + hm1)^3)/dx))/(8*dx)...
    + -(sigma*(3*(h + hm1)^3 - 3*(h + hm1)^2*(3*h - 3*hm1 + hm2 - hp1) + (h + hp1)^3))/(24*dx^4)...
    + -(lam1*sigma*(3*(h + hm1)^3 - 3*(h + hm1)^2*(3*h - 3*hm1 + hm2 - hp1) + (h + hp1)^3))/(24*dt*dx^4)...
    + (alpha*delsig*((h + hm1)^2 - (h - hm1)*(2*h + 2*hm1)))/(8*dx^2*kapa);

f = @(hm2,hm1,h,hp1,hp2,hIm2,hIm1,hI,hIp1,hIp2) mew/dt + -(lam1*sigma*((h + hm1)^2*(3*h - 3*hm1 + hm2 - hp1) + (h + hp1)^2*(3*h - hm1 - 3*hp1 + hp2) + 3*(h + hm1)^2*(h - hI - hIm1 + hm1) + 3*(h + hp1)^2*(h - hI - hIp1 + hp1) + (2*h + 2*hm1)*(3*h - 3*hm1 + hm2 - hp1)*(h - hI - hIm1 + hm1) + (2*h + 2*hp1)*(3*h - hm1 - 3*hp1 + hp2)*(h - hI - hIp1 + hp1)))/(16*dt*dx^4)...
    + (gr*p*(3*(h + hp1)^2*(sin(theta) + (cos(theta)*(h - hp1))/dx) - 3*(h + hm1)^2*(sin(theta) - (cos(theta)*(h - hm1))/dx) + (cos(theta)*(h + hm1)^3)/dx + (cos(theta)*(h + hp1)^3)/dx))/(8*dx)...
    + (sigma*(3*(h + hm1)^2*(3*h - 3*hm1 + hm2 - hp1) + 3*(h + hp1)^2*(3*h - hm1 - 3*hp1 + hp2) + 3*(h + hm1)^3 + 3*(h + hp1)^3))/(24*dx^4)...
    + (lam1*sigma*(3*(h + hm1)^2*(3*h - 3*hm1 + hm2 - hp1) + 3*(h + hp1)^2*(3*h - hm1 - 3*hp1 + hp2) + 3*(h + hm1)^3 + 3*(h + hp1)^3))/(24*dt*dx^4)...
    + -(alpha*delsig*((h + hm1)^2 + (h + hp1)^2 + (h - hm1)*(2*h + 2*hm1) + (h - hp1)*(2*h + 2*hp1)))/(8*dx^2*kapa);

fp1 = @(hm2,hm1,h,hp1,hp2,hIm2,hIm1,hI,hIp1,hIp2) -(lam1*sigma*((h + hp1)^2*(3*h - hm1 - 3*hp1 + hp2) - (h + hm1)^2*(h - hI - hIm1 + hm1) - 3*(h + hp1)^2*(h - hI - hIp1 + hp1) + (2*h + 2*hp1)*(3*h - hm1 - 3*hp1 + hp2)*(h - hI - hIp1 + hp1)))/(16*dt*dx^4)...
    + (gr*p*(3*(h + hp1)^2*(sin(theta) + (cos(theta)*(h - hp1))/dx) - (cos(theta)*(h + hp1)^3)/dx))/(8*dx)...
    + -(sigma*((h + hm1)^3 - 3*(h + hp1)^2*(3*h - hm1 - 3*hp1 + hp2) + 3*(h + hp1)^3))/(24*dx^4)...
    + -(lam1*sigma*((h + hm1)^3 - 3*(h + hp1)^2*(3*h - hm1 - 3*hp1 + hp2) + 3*(h + hp1)^3))/(24*dt*dx^4)...
    + (alpha*delsig*((h + hp1)^2 - (h - hp1)*(2*h + 2*hp1)))/(8*dx^2*kapa);

fp2 = @(hm2,hm1,h,hp1,hp2,hIm2,hIm1,hI,hIp1,hIp2) -(lam1*sigma*(h + hp1)^2*(h - hI - hIp1 + hp1))/(16*dt*dx^4)...
    + (sigma*(h + hp1)^3)/(24*dx^4)...
    + (lam1*sigma*(h + hp1)^3)/(24*dt*dx^4);
    
%initalize the Jacobian
Jacob = zeros(length(g),length(g));

%plug in boundary functions
%i = 1
Jacob(1,1) = f(g(end-2),g(end-1),g(1),g(2),g(3),hP(end-2),hP(end-1),hP(1),hP(2),hP(3));
Jacob(1,2) = fp1(g(end-2),g(end-1),g(1),g(2),g(3),hP(end-2),hP(end-1),hP(1),hP(2),hP(3));
Jacob(1,3) = fp2(g(end-2),g(end-1),g(1),g(2),g(3),hP(end-2),hP(end-1),hP(1),hP(2),hP(3));
%i = 2
Jacob(2,1) = fm1(g(end-1),g(1),g(2),g(3),g(4),hP(end-1),hP(1),hP(2),hP(3),hP(4));
Jacob(2,2) = f(g(end-1),g(1),g(2),g(3),g(4),hP(end-1),hP(1),hP(2),hP(3),hP(4));
Jacob(2,3) = fp1(g(end-1),g(1),g(2),g(3),g(4),hP(end-1),hP(1),hP(2),hP(3),hP(4));
Jacob(2,4) = fp2(g(end-1),g(1),g(2),g(3),g(4),hP(end-1),hP(1),hP(2),hP(3),hP(4));
%i = N-1
Jacob(end-1,end) = fp1(g(end-3),g(end-2),g(end-1),g(end),g(2),hP(end-3),hP(end-2),hP(end-1),hP(end),hP(2));
Jacob(end-1,end-1) = f(g(end-3),g(end-2),g(end-1),g(end),g(2),hP(end-3),hP(end-2),hP(end-1),hP(end),hP(2));
Jacob(end-1,end-2) = fm1(g(end-3),g(end-2),g(end-1),g(end),g(2),hP(end-3),hP(end-2),hP(end-1),hP(end),hP(2));
Jacob(end-1,end-3) = fm2(g(end-3),g(end-2),g(end-1),g(end),g(2),hP(end-3),hP(end-2),hP(end-1),hP(end),hP(2));
%i = N
Jacob(end,end) = f(g(end-2),g(end-1),g(end),g(2),g(3),hP(end-2),hP(end-1),hP(end),hP(2),hP(3));
Jacob(end,end-1) = fm1(g(end-2),g(end-1),g(end),g(2),g(3),hP(end-2),hP(end-1),hP(end),hP(2),hP(3));
Jacob(end,end-2) = fm2(g(end-2),g(end-1),g(end),g(2),g(3),hP(end-2),hP(end-1),hP(end),hP(2),hP(3));

%build the internal part of the matrix
for i=3 : length(g)-2
   Jacob(i,i-2) = fm2(g(i-2),g(i-1),g(i),g(i+1),g(i+2),hP(i-2),hP(i-1),hP(i),hP(i+1),hP(i+2));
   Jacob(i,i-1) = fm1(g(i-2),g(i-1),g(i),g(i+1),g(i+2),hP(i-2),hP(i-1),hP(i),hP(i+1),hP(i+2));
   Jacob(i,i) = f(g(i-2),g(i-1),g(i),g(i+1),g(i+2),hP(i-2),hP(i-1),hP(i),hP(i+1),hP(i+2));
   Jacob(i,i+1) = fp1(g(i-2),g(i-1),g(i),g(i+1),g(i+2),hP(i-2),hP(i-1),hP(i),hP(i+1),hP(i+2));
   Jacob(i,i+2) = fp2(g(i-2),g(i-1),g(i),g(i+1),g(i+2),hP(i-2),hP(i-1),hP(i),hP(i+1),hP(i+2));
end
end

function [ F ] = functionVector( g,hP,coe )
%generates our function vector

%initalize our coeiffents
p = coe(1);
gr = coe(2);
mew = coe(3);
sigma = coe(4);
lam1 = coe(5);
theta = coe(6);
alpha = coe(7);
kapa = coe(8);
delsig = coe(9);
dx = coe(10);
dt = coe(11);

%initalize the function we will use
f = @(hm2,hm1,h,hp1,hp2,hIm2,hIm1,hI,hIp1,hIp2) (sigma*((h + hm1)^3*(3*h - 3*hm1 + hm2 - hp1)...
    + (h + hp1)^3*(3*h - hm1 - 3*hp1 + hp2)))/(24*dx^4) ...
    + (mew*(h - hI))/dt - (gr*p*((h + hm1)^3*(sin(theta) ...
    - (cos(theta)*(h - hm1))/dx) - (h + hp1)^3*(sin(theta) ...
    + (cos(theta)*(h - hp1))/dx)))/(8*dx) ...
    + (lam1*sigma*((h + hm1)^3*(3*h - 3*hm1 + hm2 - hp1) ...
    - (hI + hIm1)^3*(3*hI - 3*hIm1 + hIm2 - hIp1) ...
    - (hI + hIp1)^3*(3*hI - hIm1 - 3*hIp1 + hIp2) ...
    + (h + hp1)^3*(3*h - hm1 - 3*hp1 + hp2)))/(24*dt*dx^4) ...
    - (lam1*sigma*((h + hm1)^2*(3*h - 3*hm1 + hm2 - hp1)*(h - hI - hIm1 + hm1) ...
    + (h + hp1)^2*(3*h - hm1 - 3*hp1 + hp2)*(h - hI - hIp1 + hp1)))/(16*dt*dx^4)...
    + -(alpha*delsig*((h + hm1)^2*(h - hm1) + (h + hp1)^2*(h - hp1)))/(8*dx^2*kapa);

%initalize our function vector
F = zeros(length(g),1);

%plug in boundary functions
%i = 1
F(1) = f(g(end-2),g(end-1),g(1),g(2),g(3),hP(end-2),hP(end-1),hP(1),hP(2),hP(3));
%i = 2
F(2) = f(g(end-1),g(1),g(2),g(3),g(4),hP(end-1),hP(1),hP(2),hP(3),hP(4));
%i = N-1
F(end-1) = f(g(end-3),g(end-2),g(end-1),g(end),g(2),hP(end-3),hP(end-2),hP(end-1),hP(end),hP(2));
%i = N
F(end) = f(g(end-2),g(end-1),g(end),g(2),g(3),hP(end-2),hP(end-1),hP(end),hP(2),hP(3));

%internal functions
for i=3 : length(F)-2
   F(i) =  f(g(i-2),g(i-1),g(i),g(i+1),g(i+2),hP(i-2),hP(i-1),hP(i),hP(i+1),hP(i+2));
end
end


