## Author: M Turhan ÇOBAN
## Created: 2021-10-28

## HT_cylinder_1
  h=500.0;
  kc=63.9;
  alpha=1.2e-5;
  b=0.04;
  t=[360 420 480];
  k=length(t);
  Tinf=60;
  Ti=-20;
  T=Ti-Tinf;
  Bi=h/kc*b

  f=@(x)x*besselj(1,x)-Bi*besselj(0,x);
  a1=0.0;
  b1=30.0;
  r=0;
  m=100;
  dx=b/(m-1);
  ksi=roots_bisection(f,a1,b1)
  n=length(ksi);
  TT=zeros(k,m);
  for i=1:n
    J1=besselj(1,ksi(i));
    J0=besselj(0,ksi(i));
	  Cn(i)=2.0*(Ti-Tinf)/ksi(i)*J1/(J0*J0+J1*J1);
  endfor
  Cn'
  for l=1:k 
      Fo(l)=alpha*t(l)/(b*b);
  endfor
  for l=1:k 
  theta=0.0;
  for j=1:m;
     theta=0.0;
     r=(j-1)*dx;
     for i=1:n
       theta+=Cn(i)*exp(-ksi(i)*ksi(i)*Fo(l))*besselj(0,ksi(i)*r/b);
     endfor
     TT(l,j)=theta+Tinf;
     xx(j)=r;
   endfor
   endfor
  plot(xx,TT(1,:));
   title('Solid cylinder Temperature profile t=360,420,480 s');
   xlabel('r') 
   ylabel('T* degree C') 
   hold on
   for i=2:k
     plot(xx,TT(i,:));
   endfor
   hold off