%
%   PARABOLIC_FIN1.M    Find the surface area of a parabolic fin
%
%   Assume that a parabolic fin has a shape described by
%       y = a1 + a2*r + a3*r^2
%   and it must satisfy three constraints
%      y(Ro) = H/2,   y(Ro+L) = h/2,   and   dy/dr at Ro+L = 0
%   These three constraints can be used to find the equation coefficients
%   for each value of L.
%
%   Once these are known, one can evaluate the following expression for the 
%   surface area of the fin for each L:
%      A = 2*top surface area + tip area
%      A = 4*pi* int (r*sqrt(1+(a2+2*a3*r)^2) dr  +  2*pi*(Ro+L)*h
%   The integral goes from Ro to Ro+L.  We will use the built-in integral routine 
%   to do these integrals.  An anonymous function, f(r), is used to define the
%   integrand for the integration.
%
%   Finally, a plot and summary table of results are generated to describe
%   the fin's geometry versus length.
%
%   File prepared by J. R. White, UMass-Lowell  (last update:  Nov. 2017)
%

      clear all,   close all
%
%   model constants
      Ro = 5;      % cylindrical tube's outside radius (cm)
      H = 6;       % height of fin at base of cylindrical tube (cm)
      h = 1;       % height of fin at the tip (cm)
%
%   loop over different values of fin length, L
      L = 2:2:10;   NL = length(L);   % set vector of fin lengths (cm)
      Sarea = zeros(1,NL);            % initialize storage for surface area
      Coeff = zeros(3,NL);            % initialize storage for coeff vector
      for i = 1:NL
%  
%   set up coeff matrices and solve resultant linear system
        RoL = Ro + L(i);
        A = [1   Ro   Ro^2; 1   RoL   RoL^2; 0   1  2*RoL];    b = [H/2  h/2  0]';
        x = A\b;    Coeff(:,i) = x;  % store solution vector for ith length
%
%   find surface area for given L
        a2 = x(2);   a3 = x(3);
        f = @(r) r.*sqrt(1 + (a2 + 2*a3*r).^2); % anonymous function for use in integral
        Sarea(i) = 4*pi*integral(f,Ro,RoL) + 2*pi*RoL*h; % compute integral to find area
      end
%
%   plot and tabulate results
      fprintf('\n  Parabolic_Fin1:  Summary Data for Geometry of Parabolic Fin \n\n');
      fprintf('   Length (cm)     a1            a2            a3     Surface Area (cm^2) \n');
      for i = 1:NL
        fprintf('  %6.2f    %12.3e  %12.3e  %12.3e    %13.4e  \n ',L(i),Coeff(:,i),Sarea(i));
      end
%
      style  = ['r- ';'b--';'g-.';'c: ';'m- '];  % 5 different plot styles
      for  i = 1:NL
        r = linspace(Ro,Ro+L(i),20);  % need a new r vector for each length, L
        y = Coeff(1,i) + Coeff(2,i)*r + Coeff(3,i)*r.^2;
        plot(r-Ro,y,style(i,:),'LineWidth',2),hold on          % plots +y(r)
        plot(r-Ro,-y,style(i,:),'LineWidth',2)                 % plots -y(r)
        plot([L(i) L(i)],[-h/2 h/2],style(i,:),'LineWidth',2)  % plots vertical end
      end
      axis equal; grid % use axis equal to display proper aspect ratio for real geometry
      title('Shape of Parabolic Fin for Several Fin Lengths')
      xlabel('position along fin (cm)'),ylabel('axial position (cm)')
      for i=1:NL,  text(L(i)-1.35,0.25,['L = ',num2str(L(i)),' cm']);   end
%
%  end of problem

