function THsolver_simple()

global Af Ac;
global Tin Pressure w;
global nr nz dz length;
global Rf0 Rf1 Rc0 Rc1;

close all;

% --- System parameters

Tin = 538; % Inlet temperature (K)

w = 1.0*7100*6*1000/60/60/313/126; % Flow rate (kg/s) based on total flow
			       % rate
%w = 1500e6/(35*5300)/313/126;  % Flow rate (kg/s) based on heat
			       % transfer

w = 0.5*w;

Pressure = 12.3e6; % System pressure (Pa)

% --- Rod geometry

Rf0 = 0.06e-2; % Fuel inner radius (m)
Rf1 = 0.38e-2; % Fuel outer radius (m)
Rc0 = 0.3865e-2; % Cladding inner radius (m)
Rc1 = 0.4535e-2; % Cladding outer radius (m)

length = 248e-2; % Total length (m)

pitch = 1.23e-2; % Pin pitch (m)

% --- Fuel and coolant area

%Ac = pitch^2-Rc1^2*pi;
Af = Rf1^2*pi;

% --- Number of axial and radial zones in fuel rod

nz = 100;
nr = 50;

dz = length/nz;

Z = linspace(-length/2, length/2, nz+1)';
R = [0 linspace(Rf0, Rf1, nr-3+1) Rc0 Rc1]';

% --- Calculate the area of each annulus

A = pi*(R(2:end).*R(2:end) - R(1:end-1).*R(1:end-1));

% --- Get power density

Pden = GetPowerDensity(Z, R);

% --- Solve the coolant temperature distribution

Tc = SolveCoolantTemperature(Z, Pden, A);

% --- Solve the fuel temperature distribution

Tf = SolveFuelTemperature(R, Pden, A, Tc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

display(GetFuelK(1000))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

display(GetGasK(700))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

display(GetCladK(560))

figure
%plot(0.5*(Z(1:end-1)+Z(2:end)), Pden,'rx')
hold on
plot(Z, Tc,'b-')
plot(0.5*(Z(1:end-1)+Z(2:end)), Tf(:,end),'k-')
hold off
grid on

xlabel('Axial coordinate')
ylabel('Temperature (K)')
legend('Coolant', 'Cladding','location','northwest')

figure

plot(0.5*(Z(1:end-1)+Z(2:end)), sum(Pden.*repmat(A',nz,1),2)*dz,'k.')
hold on
plot(0.5*(Z(1:end-1)+Z(2:end)), Tf(:,1),'r-')
hold off
grid on
xlabel('Axial coordinate')
ylabel('Fuel centerline temperature / power level')


figure

%plot(0.5*(R(1:end-1)+R(2:end)),Pden(ceil(nz/2),:),'rx')
hold on
plot(R*100,Tf(ceil(nz/2),:),'bx')
hold off
grid on

xlabel('Radius (cm)')
ylabel('Fuel temperature (K)')

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Tf = SolveFuelTemperature(R, P, A, Tc)
global nz nr Tin w dz Rf0 Rf1 Rc0 Rc1;

% --- R is a vertical vector containing the r-limits of the annular
%     zones
% --- P is a matrix containing the power density in each of the
%     fuel rod annular zones
% --- Tc is a vector containing the coolant temperatures at each of
%     the axial elevations (nz + 1)

% --- We'll start with an initial guess for the fuel temperature
%     distribtion

Tf = zeros(nz,nr+1)+100;

% --- Loop over axial zones

for i = 1:1:nz

    % --- Get bulk coolant temperature for this axial zone

    Tb = 0.5*(Tc(i) + Tc(i+1));

    % --- We will start at the cladding outer surface

    % --- Calculate the heat flux through the cladding surface

    Hflux = sum(P(i,1:end).*A'*dz)/(dz*2*pi*R(end));

    % --- Get heat transfer coefficient between cladding and coolant
    %     (now constant)

    hclad = GetCladH(Tb);

    % --- Calculate and store cladding outer temperature

    Tf(i, nr + 1) = Tb + Hflux/hclad;

    Tout = Tf(i, nr + 1);

    for j = nr:-1:1

      % --- Store zone inner and outer radius

      Ri0 = R(j);
      Ri1 = R(j+1);

      % --- Calculate the heat flux through the zone inner surface

      if (j > 2)
	Hflux = sum(P(i,1:j-1).*A(1:j-1)'*dz)/(dz*2*pi*Ri0);
      else
	Hflux = 0;
      end

      % --- Calculate the power produced in this zone

      Pdeni = P(i,j);

      % --- Calculate increase in temperature over this segment

      epsilon = 1.0;
      Told = Tf(i,j);

      while (epsilon > 0.01)

	% --- Get the heat conductivity of this zone

	Teff = 0.5*(Tf(i,j+1) + Tf(i,j));

	if (Ri0 >= Rc0)
	  k = GetCladK(Teff);
	elseif ((Ri0 >= Rf1) || (Ri1 <= Rf0))
	  k = GetGasK(Teff);
	else
	  k = GetFuelK(Teff);
	end

	% --- Calculate the constants C0, C1 and C2

	C0 = -Pdeni/(4*k);

	if (j > 1)
	  C1 = Pdeni*Ri0^2/(2*k) - Hflux*Ri0/k;
	else
	  C1 = 0;
	end

	C2 = Tout + Pdeni/(4*k)*Ri1^2 - C1*log(Ri1);

	% --- Calculate and store zone inner temperature

	if (j > 1)
	  Tf(i, j) = C0*Ri0^2 + C1*log(Ri0) + C2;
	else
	  Tf(i, j) = C0*Ri0^2 + C2;
	end

	epsilon = (Tf(i,j) - Told)/Told;

	Told = Tf(i,j);
      end

      % --- This inner temperature will be the outer temperature for
      %     the next zone

      Tout = Tf(i, j);

    end


    % --- Calculate heat transfer in each of the fuel rings
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Tc = SolveCoolantTemperature(Z, Pden, A)
global nz Tin w dz Af;

% --- Z is a vertical vector containing the z-limits of the axial
%     nodes
% --- Pden is a matrix containing the power density in each of the
%     fuel rod annular zones
% --- A is a matrix containing the cross sectional area for each of the
%     fuel rod annular zones

% --- We'll start with an initial guess for the coolant temperature
%     distribtion

Tc = ones(nz + 1,1)*Tin;

% --- Loop over the axial length of the rod and solve the temperature
%     of each axial zone

for i = 2:1:nz+1

    % --- Calculate amount of power flowing into this segment

    Pin = sum(Pden(i-1,:).*A'*dz);

    % --- Calculate increase in temperature over this segment

    epsilon = 1.0;
    Told = Tc(i);

    while (epsilon > 0.01)
      % --- Get specific heat capacity based on current temperature

      cp = GetSpecificHeat(0.5*(Tc(i)+Tc(i-1)));

      % --- Calculate new upper temperature

      Tc(i) = Tc(i-1) + Pin/(w*cp);

      epsilon = (Tc(i) - Told)/Told;
      Told = Tc(i);
    end
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cp = GetSpecificHeat(T)

cptbl = [4926 5061 5229 5443 5729 6136 6781 6976 6976]; % J/(kg K)
Ttbl = [538 548 558 568 578 588 598 608 1000]; % K

% --- Interpolate data to current temperature

cp = interp1(Ttbl, cptbl, T);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function k = GetFuelK(T)

  if (1 == 2)
    % --- Use constant heat conductivity

    k = 3.5; % W / (m K)

  else
    % --- Use modified NFI model (FRAPCON/FRAPTRAN)
    %     For fresh fuel, 0% gadolinium, UO2

    A = 0.0452;
    B = 2.46e-4;
    E = 3.5e9;
    F = 16361;

    k = 1/(A + B*T) + E/T^2*exp(-F/T);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function k = GetGasK(T)

  if (1 == 2)
    % --- Use constant heat conductivity

    k = 0.3; % W / (m K)

  else
    % --- Use power law (FRAPCON/FRAPTRAN)

    k = 2.531e-3*T^0.7146;

  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function k = GetCladK(T)

  if (1 == 2)
    % --- Use constant heat conductivity

    k = 17.0; % W / (m K)

  else
    % --- Use polynomial relation (FRAPCON/FRAPTRAN)

    k = 7.51 + 2.09e-2*T - 1.45e-5*T^2 + 7.67e-9*T^3;

  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function h = GetCladH(Tb)

% --- Use constant heat transfer coefficient

h = 2e6/35.0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Pden = GetPowerDensity(Z,R)
global nz nr Rf0 Rf1 dz length Af;

PRod = 1500e6/313/126;
C = PRod/(2*Af);

Pden = zeros(nz, nr);

for i = 1:nz
    for j = 1:nr
	if (R(j) >= Rf0) && (R(j+1) <= Rf1)
	  Pden(i,j) = C*cos(0.5*(Z(i) + Z(i+1))/(length/2)*(pi/2));
	else
	  Pden(i,j) = 0.0;
	end
    end
end



end