%FLUID_FLOW : Fluid flow between one or two moving plates. 
% The program solves the one-dimensional unsteady-state equation 
% of motion with a forcing function:
%
%       dv/dt=(mu/rho)*(d^2v/dy^2) + (1/rho)*f(y,t)
%
% using the function PARABOLIC1D.M.

clear; clc
disp(' Solution of the one-dimensional unsteady state nonhomogeneous')
disp(' parabolic partial differential equation of motion that simulates')
disp(' the flow of fluid between one or two moving plates.')
bcdialog = [' Lower y boundary condition:'
            ' Upper y boundary condition:'];
redo = 1;
while redo
   disp(' ')
   number_of_plates= input(' How many plates (1 or 2)? = ');
   disp(' ')
   h = input(' Height of fluid (m) = ');
   tmax = input(' Maximum time (s) = ');
   p = input(' Number of divisions in y-direction = ');
   q = input(' Number of divisions in t-direction = ');
   disp(' ')
   mu = input(' Viscosity of fluid (kg/m.s) = ');
   rho = input(' Density of fluid (kg/m^3) = ');
   disp(' ')
   disp(' 1 - No forcing function')
   disp(' 2 - Yes, there is a forcing function')
   forcing = input(' Enter your choice : ') - 1;
   if forcing
      disp(' ')
      f = input(' Name of the file containing the forcing function = ');
   end
   disp(' ')
   v0 = input(' Initial velocity of fluid (m/s) = ');
   disp(' ')
   disp(' Boundary conditions:')
   for k=1:2
       disp(' ')
       disp(bcdialog(k,:))
       disp('   1 - Dirichlet')
       disp('   2 - Neumann')
       disp('   3 - Robbins')
       bc(k,1) = input(' Enter your choice : ');
      if bc(k,1) < 3
         bc(k,2) = input(' Value = ');
      else
         disp(' u'' = (beta) + (gamma)*u')
         bc(k,2) = input(' Constant    (beta)  = ');
         bc(k,3) = input(' Coefficient (gamma) = ');
      end
   end
   u0 = v0*ones(p+1,1);
% Calculating velocity profile  
   if forcing
        [y,t,v] = parabolic1D(p,q,h/p,tmax/q,mu/rho,u0,bc,f,bc(2,2),bc(1,2),mu,rho,h);
   else
        [y,t,v] = parabolic1D(p,q,h/p,tmax/q,mu/rho,u0,bc);
   end
% Plotting the velocity profile
   figure(1);
   new_v(:,1)=v(:,1);
   count=fix((length(t)-1)/10);
   for i=1:10
       new_v(:,i+1)=v(:,1+i*count);
   end
   plot(new_v,y)
   xlabel('Velocity v, m/s')
   ylabel('Position y, m')
   switch number_of_plates
     case 1
       title('Velocity profiles above moving plate')
     case 2    
       title('Velocity profiles between two plates')
   end
   v_average=simpson(y,v(:,length(v)));
   fprintf('\n Average velocity at last t, < v > = %f \n',v_average)    
   disp(' ')
   redo = input(' Repeat calculations (0/1) ? ');
end