%% 
   clear all %clear all variables
   clc % clear output screen

   %Constants
   
   mu	  = 0.1;  %i.e. Re=100
   nx	  = 21+1;
   ny	  = nx;  %different nx & ny gives non-square
   itmax  = 10;
   
   w	= zeros(nx,ny); % form a matrix of zeros (1..nx, 1..ny) 
   p = w; % w = vorticty, p = streamfunction
   u = v = vt = ut = w; % wt = dw/dt
   
   %Initialization
   h = 1/(nx-2); %h = dx (=dy); so x-size is 1, y-size is ny/nx
   t = 0;
   r = 2/(1+ 0.8/(nx-1)); %SOR parameter
   dt = 0.2*h*h/mu; 
   
   for i= 1:nx
      for j=1:ny
	      p(i,j) = 0;
	      u(i,j) = 0;
              v(i,j) = 0;
      end
   end
   
   %% Main program starts here

   while(t<3) %time loop starts here...
   
           for i =2:nx-2
              u(i,1) = - u(i,2);
              u(i,ny)= 2*sin(pi*(i-1)*h)*sin(pi*(i-1)*h) - u(i,ny-1);
              for j = 2:ny-1
	              ut(i,j) = - u(i,j)*(u(i+1,j)-u(i-1,j))/(2*h)...
	                     - (v(i,j)+v(i+1,j)+v(i,j-1)+v(i+1,j-1))*(u(i,j+1)-u(i,j-1))/(8*h)...
	                      + mu*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-4*u(i,j))/(h*h);
              end
           end
   
           for j = 2:ny-2
              v(1,j) = - v(2,j);
              v(nx,j)= - v(nx-1,j);
              for i = 2:nx-1
	              vt(i,j) = - (u(i,j)+u(i,j+1)+u(i-1,j)+u(i-1,j+1))*(v(i+1,j)-v(i-1,j))/(8*h)...
	                      - v(i,j)*(v(i,j+1)-v(i,j-1))/(2*h)...
	                    + mu*(v(i+1,j)+v(i-1,j)+v(i,j+1)+v(i,j-1)-4*v(i,j))/(h*h);
              end
           end
   
           for i=2:nx-1
              for j=2:ny
	              u(i,j)= u(i,j) + dt*ut(i,j);
              end
           end
           for i= 2:nx
              for j= 2:ny-1 
	              v(i,j)= v(i,j) + dt*vt(i,j);
              end
           end

        %div u
           for i = 2:nx-1 
              for j = 2:ny-1 
	              w(i,j) = (u(i,j) - u(i-1,j) + v(i,j) - v(i,j-1))/(dt*h);  
              end
           end 
 
%----------% Main Poisson solver starts here: solve nabla^2 p = - w
           for it = 1:round(nx/2) % iterations start here

              %BC dp/dn=mu*d2un/dn2
              for i = 2:nx-1
                p(i,1) = p(i,2)- mu*(-v(i,4)+4*v(i,3)-5*v(i,2)+2*v(i,1))/h;
	              p(i,ny)= p(i,ny-1)...
	                      + mu*(-v(i,ny-4)+4*v(i,ny-3)-5*v(i,ny-2)+2*v(i,ny-1))/h;     
              end              
              for j = 1:ny
                p(1,j)= p(2,j)...
	                      - mu*(-u(4,j)+4*u(3,j)-5*u(2,j)+2*u(1,j))/h;
	              p(nx,j)= p(nx-1,j)...
	                      + mu*(-u(nx-4,j)+4*u(nx-3,j)-5*u(nx-2,j)+2*u(nx-1,j))/h; 
              end              
           
              % Field 
              for i = 2:nx-1
	              for j = 2:ny-1
	                p(i,j) =(1-r)*p(i,j) + ...
	                        r*0.25*(p(i+1,j) + p(i-1,j) + p(i,j+1) + p(i,j-1)...
		                      - w(i,j)*h*h );
                end
              end
           end    % iterations end here
           q = 0;
           for i = 1:nx
              for j = 1:ny
        	      q = q + p(i,j);
              end  
           end
           q = -q/((nx)*(ny));
           for i = 1:nx
              for j = 1:ny
	              p(i,j) = p(i,j) + q;
              end
           end
%-------%---- Poisson solver Ends here

           for i =2:nx-2
              for j = 2:ny-1
	              u(i,j) = u(i,j) - dt*(p(i+1,j) - p(i,j))/h;
              end
           end
           for j = 2:ny-2
              for i = 2:nx-1
	              v(i,j) = v(i,j) - dt*(p(i,j+1) - p(i,j))/h;
              end
           end
       
           t = t + dt;
      
     endwhile % time loop ends here
       
   
%       % Force calculation 
%       
%       % top, uses u = 0 at ends to save correction to int  
%       q = 0;
%       for i = 1:nx-1 
%          q = q + mu*(u(i,ny)-u(i,ny-1));          
%       end
%       
%       % bottom
%       qb = 0;
%       for i = 1:nx-1
%           qb = qb - mu*(u(i,2)-u(i,1));
%       end
%       
%       % sides
%       
%       qs = 0;
%       for j = 2:ny-1
%           qs = qs + h*0.5*(p(1,j)+p(1,j) - p(nx-1,j)-p(nx,j));
%       end
%
%          writeln(nx:4,1/(nx-1):12:4,1/((nx-1)*(nx-1)):12:8,
%           q:12:4,qs:12:4,qb:12:4,-qs-qb:12:4);
%       
%       formatSpecA='nx=%4d; 1/(nx-1)=%10.5f; 1/((nx-1)*nx-1))=%12.8f;q=%12.4f;qs=%12.4f; qb=%12.4f; -qs-qb=%12.4f\n ';
%       fprintf(formatSpecA, nx,1/(nx-1),1/((nx-1)*(nx-1)),q, qs, qb, -qs-qb);
%       % End of Force calculation
       
 
    
    %Plot 'u' 
    figure(1)   
    tx = linspace (0, (nx-1)*h, nx)';
    ty = linspace (0, (ny-1-0.5)*h, ny)';
    
    [xx, yy] = meshgrid (tx, ty);
    mesh (tx, ty, u');
    xlabel ("x");
    ylabel ("y");
    zlabel ("u(x,y)");
    title ("u plot");
    % End of plot
    
    %Plot 'v'
    figure(2)
    ty = linspace (0, (nx-1)*h, nx)';
    tx = linspace (0, (ny-1-0.5)*h, ny)';
    
    [xx, yy] = meshgrid (tx, ty);
    mesh (tx, ty, v');
    xlabel ("x");
    ylabel ("y");
    zlabel ("v(x,y)");
    title ("v plot");
    % End of plot
    
    %Plot 'p'
    figure(3)
    tx = linspace (0, (nx-1-0.5)*h, nx)';
    ty = linspace (0, (ny-1-0.5)*h, ny)';
    
    [xx, yy] = meshgrid (tx, ty);
    mesh (tx, ty, p');
    xlabel ("x");
    ylabel ("y");
    zlabel ("p(x,y)");
    title ("p plot");
    % End of plot
    
    % Final output
    for j = 1:ny
      formatSpec3='y=%10.5f; u(0.5,y)=%10.5f\n ';
      fprintf(formatSpec3, h*(j-1-0.5), u(round((nx-1)/2),j) );
    end
    
    
      
   
