%
%   CYLFINSH.M   Heat Transfer in a Cylindrical Fin (Shooting Method)
%
%   This file solves the cylindrical fin heat transfer problem using the 
%   Shooting Method. The base problem is defined via the following equation:
%
%     x^2*u'' + x*u' - ALF2*x^2*u = 0    where   ALF2 = 2*h*rs^2/[k*thk]
%
%     with B.C.   u(a) = 1    and    u'(b) = 0
%      and        a = rw/rs   and    b = rs/rs = 1
%
%     where   u = normalized temp = [T(r) - Tinf]/[Tw - Tinf]
%             x = normalized distance = r/rs
%
%     with   rw, rs = inside and outside radius of fin, respectively
%            thk = thickness of fin
%     and    h, k, Tw, and Tinf are all given quantities (fixed)
%
%   From the normalized solution we can construct absolute profiles (if desired):
%       T(r)  =  Tinf + u(x)[Tw - Tinf]
%       T'(r) =  u'(x)[Tw - Tinf]/rs
%
%   The above 2nd order ODE is first converted into a set of two first order
%   ODEs in standard format (see CYLFINSHF.M) and then solved as an IVP as part 
%   of the Shooting Method algorithm.  The basic Shooting method has been implemented
%   within the BVP2SH.M function file (requires linear BCs).
%
%   The above development is given as part of the course notes in the Math
%   Methods course (10/24.539).
%
%   File prepared by J. R. White, UMass-Lowell (Aug. 2003)
%

%
%   getting started
      clear all, close all, nfig = 0;
      global alpha
%
%   basic data for the problem
      rw = 1/12;    rs = 1.5/12;    % inside and outside radius (ft)
      thk = .0625/12;               % thickness of fin (ft)
      Tw = 200;    Tinf = 70;       % inside wall and ambient temps (F)
      h = 20;                       % heat transfer coeff (BTU/hr-ft^2-F)
      k = 75;                       % thermal conductivity (BTU/hr-ft-F)
%
%   derived constants
      a = rw/rs;   b = rs/rs;
      alf2 = (2*h*rs*rs)/(k*thk);    alpha = sqrt(alf2);
      Qideal = 2*pi*h*(rs*rs-rw*rw)*(Tw-Tinf);
%
%   write base data to output file
      fid = fopen('cylfinsh.out','w');
      fprintf(fid,'\n *** CYLFINSH.OUT ***  Data and Results from CYLFINSH.M \n');
      fprintf(fid,'\n \nBASIC DATA FOR PROBLEM \n');
      fprintf(fid,'Inside and outside radius:  rw =  %6.2f ft \t rs =   %5.2f ft \n',rw,rs);
      fprintf(fid,'Thickness of fin:           thk = %6.2f ft \n',thk);
      fprintf(fid,'Inside/ambient temps:       Tw =  %6.2f F  \t Tinf = %5.2f F \n',Tw,Tinf);
      fprintf(fid,'Heat transfer coeff:        h =   %6.2f Btu/hr-ft^2-F \n',h);
      fprintf(fid,'Thermal conductivity:       k =   %6.2f Btu/hr-ft-F \n',k);
%
%   initialize zbc array and call BVP2SH.M to solve problem
      zbc = [1 0 1.0 ;   % u(a) = 1
             0 1 0.0];   % u'(b) = 0
      [x,z] = bvp2sh('cylfinshf',[a b],zbc);
      [nr,nc] = size(z);
%
%   now have desired solution (do auxiliary calcs and plotting)
%   calc heat transferred and fin efficiency (edit key parameters)
      Qactual = -k*(2*pi*a*thk)*(Tw-Tinf)*z(1,2);
      Ttip = Tinf + (Tw-Tinf)*z(nr,1);
      fineff = Qactual/Qideal;
      fprintf(fid,'\n \nFINAL SHOOTING METHOD RESULTS \n');
      fprintf(fid,'Wall Temp (F) =     %8.3f \n',Tw);
      fprintf(fid,'Ambient Temp (F) =  %8.3f \n',Tinf);
      fprintf(fid,'Tip Temp (F) =      %8.3f \n',Ttip);
      fprintf(fid,'Qideal (BTU/hr) =   %8.3f \n',Qideal);
      fprintf(fid,'Qactual (BTU/hr) =  %8.3f \n',Qactual);
      fprintf(fid,'Fin Eff =           %8.3f \n',fineff);
%
%   plot normalized profiles
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(x,z(:,1),'LineWidth',2)
      title('CylFinSH:  Normalized Temp Profile for Cylindrical Fin (Shooting Method)')
      grid,ylabel('temperature')
      subplot(2,1,2),plot(x,z(:,2),'LineWidth',2)
      title('CylFinSH:  Normalized Temp Gradient for Cylindrical Fin (Shooting Method)')
      grid,xlabel('normalized distance'),ylabel('temp gradient')
%
%   close output file
      fclose(fid);
%
%   end simulation