function X = radial_heat_conduction2(N)

% This function can be used to solve Part II, Problem 2.  The number of
% roots is now specified as an input into the function rather than being
% specified by editing the function.  To solve the given problem, we simply
% ran the function repeatedly, increasing the number of roots with each
% iteration until the difference temperature between successive iterations
% was less than 10^-4 K.  The answer is 3 roots for t=5 sec and 2 roots for
% t=30 sec.  For the time-dependent calculations, 11 roots were required
% until the change at all thermocouples became less than 10^-4 K.  To be 
% safe, we will use 15 roots when obtaining theoretical data for the
% compound.

% That is, to obtain theoretical data for II.3, all of our physical data
% was entered into radial_heat_conduction.m and the function was then run
% for N=5.

save 'matlab.mat' N % need to prevent N from being cleared
clear;
load 'matlab.mat' N 

k=125; %Thermal conductivity (W/m-K)
rho=8500; %Density of brass (kg/m^3)
Cp=380; %Heat capacity of brass (J/kg-K)
alpha=k/(rho*Cp); %Thermal diffusivity (m^2/s)
I=18.11; %Current (A)
V=2.9; %Voltage (V)
Q=I*V; %Heater power (W)
a=0.007; %Inner radius of disk (m)
b=0.055; %Outer radius of disk (m)
L=0.0032; %Thickness of disk (m)
A=2*pi*a*L;
Q_dp=Q/A; %Q_double prime (W/m^2)
U=9803.161; %Overall heat transfer rate (W/m2-K)
% N=1; %Number of roots desired is now an input of the function
xn=find_roots(N,a,b,U,k);

%For T versus r
t=[5 30]; %Times (s)
rstep=0.001; %Radius step (m)
rinit=0.007; %Initial radius (m)
rfin=0.007; %Final radius (m)

i=1;
j=1;

for j=1:length(t);
    R=rinit;
    i=1;
    for i=1:((rfin-rinit)/rstep+1);         
    sum=0;
        for n=1:length(xn);  
        Co=besselj(0,R*xn(n))*(xn(n)*bessely(1,a*xn(n)))-bessely(0,R*xn(n))*(xn(n)*besselj(1,a*xn(n)));
        F=(xn(n)^2+(U/k)^2)*(xn(n)*besselj(1,a*xn(n)))^2-(xn(n))^2*(xn(n)*besselj(1,b*xn(n))-U/k*besselj(0,b*xn(n)))^2;
        rest=(xn(n)*besselj(1,b*xn(n))-U/k*besselj(0,b*xn(n)))^2;
        sum=sum+exp(-alpha*xn(n)^2*t(j))*(Co/F*rest);
        end;
    Theta_r(i,1)=R;
    Theta_r(i,j+1)=a*Q_dp/b/U+a*Q_dp/k*log(b/R)+(pi*Q_dp/k)*sum;
    R=R+rstep;
    end;
end;

Theta_r

%For T versus t
R=[.007 .010 .020 .030 .040 .050]; %Locations of T1-T6 (m)
tstep=1; %Time step (s)
ttot=500; %Total time (s)

i=1;
j=1;

for j=1:6;
    t=0;
    i=1;
    for i=1:(ttot/tstep);         
    sum=0;
        for n=1:length(xn);  
            Co=besselj(0,R(j)*xn(n))*(xn(n)*bessely(1,a*xn(n)))-bessely(0,R(j)*xn(n))*(xn(n)*besselj(1,a*xn(n)));
            F=(xn(n)^2+(U/k)^2)*(xn(n)*besselj(1,a*xn(n)))^2-(xn(n))^2*(xn(n)*besselj(1,b*xn(n))-U/k*besselj(0,b*xn(n)))^2;
            rest=(xn(n)*besselj(1,b*xn(n))-U/k*besselj(0,b*xn(n)))^2;
            sum=sum+exp(-alpha*xn(n)^2*t)*(Co/F*rest);
        end;
    Theta_t(i,1)=t;
    Theta_t(i,j+1)=a*Q_dp/b/U+a*Q_dp/k*log(b/R(j))+(pi*Q_dp/k)*sum;
    t=t+tstep;
    end;
end;

Theta_t(1,:)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Xn = find_roots(n,a,b,U,k)

Xo(1)=50;
for j=2:n
    Xo(j)=Xo(j-1)+65;
end

for i=1:length(Xo)
    xo=Xo(i);
    X=fsolve(@roots,xo,[],a,b,U,k);
    Roots(i)=X;
end

Xn=Roots;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F = roots(x,a,b,U,k)

F = (x*besselj(1,a*x))*(x*bessely(1,b*x)-U/k*bessely(0,b*x))-(x*bessely(1,a*x))*(x*besselj(1,b*x)-U/k*besselj(0,b*x));