clear all;

global X1 elcon;

for m=4:4
    nnodes=2^(m+2)+1
%*********** Geometry ***********************************
PHIL = [1];
PHIR = [nnodes];
UL = [1];
UR = [nnodes];
PL = [1];
PR = [nnodes];

%*********** Physical Parameters ***********************
%nnodes = 9;
left = 0;
right = 500/1000;
X1 = linspace(left,right,nnodes);
B = 50/1000;
C = 100/1000;
n = 2;
L = right-left;
A = B*sin(n*pi*X1/L) + C;
Uin = 5;
rho_2 = 960/2;
rho_4 = 960/4;
po = (200)*ones(nnodes,1);
zero = 0;
one = 1;
[Xnodes Ynodes] = size(X1);
phi = ones(Xnodes.*Ynodes,1)*0;
u = ones(Xnodes.*Ynodes,1)*0;
p = ones(Xnodes.*Ynodes,1)*0;

%********************* Apply Boundary Conditions to Initial Condition for Set 0 ********************
for ii=1:1:length(PHIR)
phi(PHIR(ii)) = +zero;
end

%********************* Equations for Set 0*****************
phi0 = +asres1D(1,[],A,-1,'A3011L',phi);

%Apply Neumann boundary condition to boundary PHIL.
phi0(1) = phi0(1) -A(1).*Uin;

%Apply Dirichlet boundary condition to boundary PHIR.
for ii=1:1:length(PHIR)
phi0(PHIR(ii)) = 0;
end

%***************** Jacobians for Set 0**********
%Equation phi by phi
phi_jac0 = +asjac1D(1,[],A,-1,'A3011L',[]);
DIFFPHI=phi_jac0;

%Apply Dirichlet boundary condition to boundary PHIR.
for ii=1:1:length(PHIR)
phi_jac0(PHIR(ii),:) = 0;
phi_jac0(:,PHIR(ii)) = 0;
phi_jac0(PHIR(ii),PHIR(ii)) = 1;
end

%***************** Solver ********************************************
dphi = phi_jac0\(-phi0);

%************* Solve Iterate *****************************************
phi = phi + dphi

%********************* Apply Boundary Conditions to Initial Condition for Set 1 ********************

%********************* Equations for Set 1*****************
u1 = +asres1D(1,[],[],1,'A200L',u)...
+asres1D(1,[],[],0,'A201L',phi);

%Apply Neumann boundary condition to boundary UL.
u1(1) = u1(1) +zero;

%Apply Neumann boundary condition to boundary UR.
u1(nnodes) = u1(nnodes) +zero;

%***************** Jacobians for Set 1**********
%Equation u by u
u_jac1 = +asjac1D(one,[],[],1,'A200L',[]);
DIFFU=u_jac1;

%***************** Solver ********************************************
du = u_jac1\(-u1);

%************* Solve Iterate *****************************************
u = u + du

%********************* Apply Boundary Conditions to Initial Condition for Set 2 ********************

%********************* Equations for Set 2*****************
%p2 = +asres1D(1,[],[],1,'A200L',p)...
%-asres1D(1,[],[],1,'A200L',[p])...
%+asres1D(rho_2,[],u,1,'A3000L',u);

p2 = +asres1D(1,[],[],1,'A200L',p)...
-asres1D(1,[],[],1,'A200L',po)...
+asres1D(rho_4,[],phi,-1,'A3101L',phi);

%Apply Neumann boundary condition to boundary PL.
p2(1) = p2(1) +zero;

%Apply Neumann boundary condition to boundary PR.
p2(nnodes) = p2(nnodes) +zero;

%***************** Jacobians for Set 2**********
%Equation p by p
p_jac2 = +asjac1D(one,[],[],1,'A200L',[]);
DIFFP=p_jac2;

%***************** Solver ********************************************
dp = p_jac2\(-p2);

%************* Solve Iterate *****************************************
p = p + dp
if m==1
    pm8=p
elseif m==2
    pm16=p
elseif m==3
    pm32=p
elseif m==4
    pm64=p
elseif m==5
    pm128=p
end
%***Computing enorm****
enormPHI=0.5*phi'*DIFFPHI*phi
enormP=0.5*p'*DIFFP*p
enormU=0.5*u'*DIFFU*u

figure(1)
hold on
plot(X1,u)
xlabel('X1','FontSize',16),ylabel('u','FontSize',16)
title('Centerline Velocity','Fontsize',16)
figure(2)
hold on
plot(X1,p,'r')
xlabel('X1','FontSize',16),ylabel('p','FontSize',16)
title('Centerline Pressure','FontSize',16)
figure(3)
hold on
plot(X1,phi,'m')
xlabel('X1','FontSize',16),ylabel('phi','FontSize',16)
title('Centerline Potential','FontSize',16)
end

%END 
