clear all;clc;
% FEm.PSE template, steady conduction with variable conductivity
% Tapered and uniform cross-section fin
% Convection applied on the left boundary instead of radiation
% convection applied as an element data instead of the global constant

% load the FE matrix library
load femlib;

global X1;        % array for node coordinates
XL = [1];
for step = 3:1:7
%*********** Physical Parameters ***********************
nnodes = 2^step+1;
XR = (2^(step-2)+1):nnodes; % nodes on the fin
M=nnodes-1
CYLelmts = floor((nnodes-1)/4);


% constant data
sigma = 5.67e-08;   % (W/m2K4)
epsilon = 0.8;     % (emissivity)
HC = 113.56;       % heat transfer coefficient (W/m2K), = 20Btu/hrft2F
HF=HC;             % (W/m2K)
HF=ones(M,1)*HF;
TRC= 1268;         % in Kelvin, = 1500F
TRF=294.26;        % in Kelvin, = 70F
RC=1.5/39.3708;      % inner radius of the cylinder(in mts) = 2in. 
Rright=2.5/39.3708;  % outer radius(in mts) = 4 in.
COND=186;% (W/mK)= 107.5 Btu/hrftF; (Al alloy 2024-T6 at 400K, melting point 775K)
% Note:- The fin material should have a high thermal conductivity to minimize
% temperature variations from its base to tip. Hence Aluminium alloy is used 
% as a fin and cylinder material(Al has relative cost and weight advantage 
% over copper). If we can have infinite thermal conductivity
% then entire fin would have been at the same temperature.

X1 = linspace(RC,Rright,nnodes);
% distributed data at the nodes
 R = X1' ;
 
 % Computation for uniform cross-section fin
 CYLelmts = floor((nnodes-1)/4); % 1/4th of the total elmts are in the cylinder
 FINelmts = (nnodes-1)-CYLelmts;
 CYLwidth = 0.33/39.3708; % (in mts.)
 FINthknss = 0.15/39.3708; % (in mts.)
 cylinder = ones(CYLelmts,1)*CYLwidth;
 fin = ones(FINelmts,1)*FINthknss;
 DELTAZ = cat(1,cylinder,fin);

%Computation for tapered-fin starts
%CYLwidth = 0.33/39.3708;
%FINthknss=0.15/39.3708;
%delR=(Rright-RC)/M;
%DELTAZ(1:M/4,1)=CYLwidth ;
%for ii=(M/4+1):M
%DELTAZ(ii,1)=2*(FINthknss/2)*(nnodes-ii)*delR/((3/4)*(Rright-RC));
%end
%% Note that the tapered fin relation is obtained using
%% similarity triangle concept

fcv=ones(M,1);
for i=1:(M/4)
    fcv(i,1)=fcv(i,1)*0;
end
HF=HF.*fcv; %%or use
%HF(1:CYLelmts)=0;

%  initial condition
%Q = (linspace(RC,Rright,nnodes))';   
Q = zeros(nnodes,1);

%  set max iteration count and convergence criteria
maxit = 10;  %  max iterations
eps = 1e-8; % convergence tolerance

%  iteration loop
for i = 1:maxit;
    
    res= asres1D(COND,DELTAZ,R,-1,'A3011L',Q);   
    restemp = asres1D(2,HF,R,1,'A3000L',Q-TRF);
        res=res+restemp;
    res(1,1)=res(1,1)+DELTAZ(1).*sigma.*epsilon.*RC.*(Q(1).^4-TRC^4);  
    % res(1,1)=res(1,1)+(HC*RC*DELTAZ(1)*(Q(1)-TRC));    % for convection
    % BC on the left side
    jac = asjac1D(COND,DELTAZ,R,-1,'A3011L',[]);
         jactemp = asjac1D(2,HF,R,1,'A3000L',[]);
    jac=jac+jactemp;    
    DIFFHBC=jac;
    
    jac(1,1)=jac(1,1)+4.*DELTAZ(1).*sigma.*epsilon.*RC.*(Q(1)).^3;
   
   % for convection BC on the left side
    %jac(1,1)=jac(1,1)+(HC*RC*DELTAZ(1)); 
    dQ=jac\(-1*(res));
        
    %************* Solve Iterate *****************************************
Q = Q + dQ;
dQmax = max(abs(dQ));
if dQmax < eps
    iteration =i
	break
end
end
% plot the results with X1 in inches
X1_P=X1*39.3708;
figure (1)
hold on
if step == 3
    plot(X1_P,Q)
else 
end
% compute Energy Norm using [DIFF + HBC] (without Dirichlet BC)
format long; % display Enorm in long format
MaxQ=max(Q)
Enorm = 0.5*Q'*DIFFHBC*Q
end
plot(X1_P,Q,'r')
xlabel('radial distance (in.)','Fontsize',14); ylabel('Temperature (K)','Fontsize',14)
title('Temperature distribution', 'Fontsize',16)
