!-------------------------------------------------------------------
! calcul de la solution de l'équation de la chaleur axisymétrique 1d
! à l'instant t et au rayon r, dans un disque plein centré en 0, de rayon Rmax.
!
! Les variables D (coefficient de diffusion), Rmax, Tinit (donnée initiale constante),
! Tb (donnée de bord constante) sont supposées connues (variables globales)
!
! la solution est calculée avec 20 modes
!
! les zéros de la fonction de Bessel de première espèce J0 (stockées dans z) ont été obtenus
! avec maple, commande BesselJZeros(0, 1 .. 20)
!
!  VERSION
!  	15/10/2022 par Luc Mieussens
!-------------------------------------------------------------------

function T1d_axi(t,r) result(y)

  real(pr), intent(in) :: t, r
  real(pr) :: y

  integer :: n
  real(pr), dimension(1:20) :: z, c

  !--- zeros de J0
  z = (/ 2.404825558, 5.520078110, 8.653727913, 11.79153444, 14.93091771 &
       & , 18.07106397, 21.21163663, 24.35247153, 27.49347913, 30.63460647 &
       & , 33.77582021, 36.91709835, 40.05842576, 43.19979171, 46.34118837 &
       & , 49.48260990, 52.62405184, 55.76551076, 58.90698393, 62.04846919/)

  !--- coeff cn
  do n=1,20
     c(n) = 2 * (Tinit - Tb) / ( z(n) * Bessel_JN(1,z(n)) )
  end do

  !--- T exacte
  y  = Tb
  do n=1,20
     y = y + c(n) * exp(-z(n)**2 * D/Rmax**2 * t) * Bessel_JN(0,r*z(n)/Rmax)
  end do

end function T1d_axi


