% ALL COMMENTS IN CAPS ARE FOR GUIDING YOU TO WRITE A PROPER PROGRAM
% THE COMMENTS THAT ARE IN ALL CAPS ARE NOT EXPECTED WHEN YOU SUBMIT 
% YOUR PROGRAM.

% SEE HOW THERE IS AN INTRODUCTION PART AND IT HAS THESE
% COMPONENTS OF REVISED, PURPOSE, KEYWORD, AUTHOR
%% Introduction
% Revised:
%       May 7, 2020

% Purpose
%       This program will demonstrate the appropriate format to be
%       used for submissions. This is not a MATLAB tutorial. In addition,
%       it demonstrates finding the flow rate in a pipe when 
%       the velocity vs radius data is given.

% Keywords
%       Sample Format; Water Flow; Regression; Flow rate

% Author
%       Autar Kaw
%       Semester: Fall 2015

% SEE THE USE OF CLC, CLF AND CLEAR ALL
% SEE THE IDENTIFICATION OF THE AUTHOR AND PROJECT TITLE
clc
clf
clear all
disp ('Computational Methods')
disp ('EML3041')
disp ('Fall 2015')
disp ('Autar Kaw')
disp('Project Name: Flow rate in a pipe')
disp('*******************************************************************')
disp('  ')

%% Problem 1
% SEE HOW EACH PROBLEM IS IDENTIFIED EVEN IF IT IS NOT MATLAB RELATED
% DO NOT FORGET THE SPACE AFTER %%
% See attached sheet for typed input data
disp('Problem 1')
disp('See typed document in the report')
disp('*******************************************************************')
disp(' ')
% SEE HOW THE COMMENTS ARE WRITTEN AND THEY ARE SUFFICIENT BUT NOT
% REDUNDANT
% THEY ARE ALWAYS IN THEIR OWN LINE NOT APPENDED TO A MATLAB STATEMENT
%% Problem 2
% Attached is the data taken in the lab
% velocity (ft/s) vs radial location (ft) data
radial =   [0  0.083  0.17  0.25  0.33  0.42  0.5];
velocity= [10  9.72   8.88  7.5     5.6      3.1    0];
% Radius of pipe (inches)
Radius=6;
% APPRECIATE THE IDENTIFICATION OF PROBLEM NUMBER AND 
% USE OF DISP STATEMENTS
% ALL INPUT DATA IS DISPLAYED USING PROPER DISP AND FPRINTF STATEMENTS
disp ('Problem 2')
disp('_________________________________')
disp('   Radial    Velocity' )
disp('  Location     ' )
disp('    (ft)      (ft/s)' )
disp('_________________________________')
dataval=[radial;velocity]';
disp(dataval)
disp('_________________________________')
fprintf('The radius of the pipe is =%g inches',Radius)
disp('   ')
disp('*******************************************************************')
disp ('  ')

%% Problem 3
% Changing units of needed variables to USCS 
Radius_ft=Radius/12;

% SEE HOW LINES ARE SKIPPED SO THAT THE MFILE CAN BE READ NICELY AS WELL
%% Problem 4
% SEE HOW THE PROGRAM IS WRITTEN CLEARLY WITH VARIABLE NAMES THAT 
% MAKE SENSE, COMMENTS PRECEEDING NEW VARIABLES, ETC
% THIS PROBLEM IS AN EXAMPLE OF A CALCULATION PROBLEM
% Using regression formula to find the velocity profile
n=length(radial);
% Using the regression formula for velocity profile
C_Numer=sum(velocity.*(1-radial.^2/Radius_ft^2));
C_Denom=sum((1-radial.^2/Radius_ft^2).^2);
C=C_Numer/C_Denom;
disp('Problem 4')
fprintf('The velocity profile is %g(1-r^2/%g)',C,Radius_ft^2)
disp('   ')
disp('*******************************************************************')
disp ('  ')

%% Problem 5
% Plot of velocity profile as a function of radial location
% defining radial locaton, r as a symbolic variable 
syms r

% Velocity profile
vel_profile=C*(1-r^2/Radius_ft^2);

% Determing points for the plot
r_val=radial(1):(radial(n)-radial(1))/1000:radial(n);
v_val=subs(vel_profile,r,r_val);

% LOOK AT HOW THE PLOT IS LABELED WITH AXES AND LEGEND
% plotting radial location vs velocity
plot(radial,velocity,'o',r_val,v_val,'-');
xlabel('Radial Location, ft')
ylabel('Velocity, ft/s')
title ('Velocity vs radial location plot')
legend('Experimental data','Regression Curve')

disp('Problem 5')
disp('See Figure 1 for velocity profile')
disp('*******************************************************************')
disp ('  ')

%% Problem 6
% Finding the flow rate from the regression curve
% Integrating to find the flow rate from Equation (1)
flow_rate=vpaintegral(2*pi*r*vel_profile,r,0,Radius_ft);
flow_rate=double(flow_rate);
disp('Problem 6')
fprintf ('The flow rate from the regression curve is=  %g  ft^3/s',...
    flow_rate)
disp('  ')
disp('*******************************************************************')
disp ('  ')

%% Problem 7
% Finding the flow rate from the average velocity x Area method
% Average Velocity
avg_vel=mean(velocity);
% Area of pipe
Area=pi*Radius_ft^2;
flow_rate=avg_vel*Area;
flow_rate=double(flow_rate);
disp('Problem 7')
fprintf('The flow rate from the average velocity method is= %g ft^3/s',...
    flow_rate)
disp('  ')
disp('*******************************************************************')
disp ('  ')

%% Problem 8
% Using trapezoidal rule with unequal segments to find flow rate
% Using equation (1) to find te integrand
flow_rate_alt=0;
for i=1:1:n-1
    fun_up=2*pi*radial(i+1)*velocity(i+1);
    fun_low=2*pi*radial(i)*velocity(i);
    flow_annulus=(radial(i+1)-radial(i))/2*(fun_up+fun_low);
    flow_rate_alt=flow_rate_alt+flow_annulus;
end

disp('Problem 8')
disp('Using trapezoidal rule with unequal segments to find flow rate')
fprintf ('The flow rate from an alternative method is=  %g  ft^3/s',...
    flow_rate_alt)
disp('  ')
disp('********************************************************************')
disp ('  ')

%% Problem 9
% See typed document in the report
disp('Problem 9')
disp('See typed document in the report')
disp('*******************************************************************')
disp ('  ')