%
%   Fuel_Temp_FD2_mx.M   Function file to compute the coefficient matrices (A & B) 
%                          for the cylindrical fuel pin heat transfer problem.  
%
%   This routine is called from Fuel_Temp_FD2.m.  Several variables are passed 
%   into this routine via Matlab's global statement.  The position vector, r, and 
%   the thermal conductivity vector, K, are passed into the subroutine via input
%   arguments and the resultant A and B matrices are returned.
%
%   File prepared by J. R. White, UMass-Lowell (Oct. 2003)
%

      function [a,b] = abmatrix(r,K)
      global  qqq Ts Ks h h2
%
%   initialize variables
      N = length(K);   a = zeros(N,N);   b = zeros(N,1);
%
%   matrix eqns. --  interior points
      for i = 2:N-1 
        rm = r(i)-h/2;   rp =r(i)+h/2;                                                      
        a(i,i-1) = rm*(K(i) + K(i-1))/2;    
        a(i,i) = -(rm*K(i-1) + 2*r(i)*K(i) + rp*K(i+1))/2;                                                     
        a(i,i+1) = rp*(K(i) + K(i+1))/2;     
        b(i) = -qqq*r(i)*h2;     
      end                                                                       
%
%   matrix eqns.  --  left boundary  
      a(1,1) = -(K(1) + K(2))/2;   a(1,2) = (K(1) + K(2))/2;   
      b(1) = -qqq*h2/4;    
%
%   matrix eqns.  --  right boundary                                                            
      rm = r(N)-h/2;   rp =r(N)+h/2;                                                      
      a(N,N-1) = rm*(K(N) + K(N-1))/2; 
      a(N,N) = -(rm*K(N-1) + 2*r(N)*K(N) + rp*Ks)/2;                                                  
      b(N) = -qqq*r(N)*h2 - Ts*rp*(K(N)+Ks)/2;
%
%   end of routine
%
