! Module permettant de partager une variable indiquant
! la convergence de l'algorithme de Newton.
module convergence_module
  implicit none

  logical :: convergence
end


! Algorithme de Newton
real function newton(f, x0, eps)
  use convergence_module
  implicit none

  real :: y, dy, x, dx, x0, eps
  integer :: imax, i

  ! La convergence étant rapide, il n'est pas nécessaire de prendre
  ! un imax très grand.
  ! Penser à bien réinitialiser la variable de convergence.
  x = x0
  imax = 1000
  convergence = .false.

  do i = 1, imax
    call f(x,y,dy)
    dx = - (y / dy)

    ! Si la variation de x est faible, on a atteint la convergence
    if (abs(y) <= eps) then
      convergence = .true.
    endif

    x = x + dx

    ! Sortir si l'algorithme a convergé
    if ( convergence ) exit
  enddo

  newton = x
end


! Fonction calculant la dérivée d'une fonction en un point de
! manière approchée.
real function approx_dy(f, x)
  implicit none

  real :: x, dx, a, b, f
  
  dx = 1e-4

  a = f(x-dx/2)
  b = f(x+dx/2)

  approx_dy = (b - a)/dx
end


! Fonction calculant un sinus et sa dérivée de manière précise
! (en utilisant sin' = cos)
subroutine sine_precise(x, y, dy)
  implicit none

  real :: x, y, dy

  y  = sin(x)
  dy = cos(x)
end


! Fonction calculant le sinus et sa dérivée en un point de 
! manière approchée en utilisant approx_dy
subroutine sine_approx(x, y, dy)
  implicit none

  real :: x, y, dy, approx_dy

  y  = sin(x)
  dy = approx_dy(sin, x)
end


! Fonction calculant exp(x) et sa dérivée (identique) de manière exacte
subroutine exp_precise(x, y, dy)
  implicit none

  real :: x, y, dy

  y  = exp(x)
  dy = exp(x)
end


! Fonction calculant 1/2 - tanh(x-1) et sa dérivée de manière exacte
subroutine tanh_precise(x, y, dy)
  implicit none
  
  real :: x, y, dy

  y  = 0.5 - tanh(x-1.0)
  dy = -1.0 / (cosh(x-1.0)**2.0)

end


! Module stockant les constantes nécessaires à la fonction calculée
! dans le TP, données dans l'énoncé
module const_module
  implicit none

  real :: const_g = 9.81
  real :: const_l = 0.1 
  real :: const_m = 0.006 
  real :: const_R = 8.32
  real :: const_S = 0.0001 
  real :: theta_0 = 1.0
  real :: const_P = 200
  real :: const_T = 300
  real :: const_q = 1.0

  ! Paramètres à déduire des autres
  real :: const_x, const_gamma, const_omega2, const_Tc

  ! Paramètres de la fonction theta_eq(T,E)
  real :: T, E

end


! Fonction calculant f(theta), la fonction calculée en question 5 du TP
real function ftheta(theta)
  use const_module
  implicit none

  real :: theta

  ! On sépare le calcul sur trois lignes pour éviter de dépasser 80
  ! caractères par ligne à l'impression...
  ftheta = (sin(theta) + (const_q * E) / (const_m * const_g))
  ftheta = ftheta * (theta_0 ** 2 - theta ** 2) / theta
  ftheta = ftheta - const_gamma * T / const_omega2
end


! Subroutine calculant f(theta) et sa dérivée de manière approchée
subroutine ftheta_approx(x, y, dy)
  implicit none

  real :: x, y, dy 
  real :: ftheta, approx_dy

  y  = ftheta(x)
  dy = approx_dy(ftheta, x)

end


! Programme principal
program main
  use convergence_module
  use const_module
  implicit none

  integer :: i
  real :: newton, res
  external :: ftheta_approx

  ! Calcul des constantes de l'équation
  const_x = (const_P * theta_0 * const_l * const_S) / (const_R * const_T)
  const_gamma = (2 * const_x * const_R) / (const_m * (const_l ** 2.0))
  const_omega2 = const_g / const_l
  const_Tc = (const_m * const_g * const_l * (theta_0 ** 2.0)) / (2.0 * const_x * const_R)

  ! Création de la courbe
  open(1, file='theta_E.res')

  T = const_Tc / 3.0
  ! On prend une bonne valeur initiale de x0
  res = -0.8

  do i = 0, 200
    E = real(i) / 1000.0 - 0.1
    
    ! Pour chaque iteration, on passe l'ancien résultat en valeur initiale de x
    ! Ceci permet de voir l'hystérèse
    res = newton(ftheta_approx, res, 1e-5)
    if(convergence) then
      write(1,*) E, ' ', res
    endif
  enddo

  ! On refait la même boucle mais décroissante cette fois
  do i = 0, 200
    E = 0.1 - real(i) / 1000.0 
    
    res = newton(ftheta_approx, res, 1e-5)
    if(convergence) then
      write(1,*) E, ' ', res
    endif
  enddo

  close(1)
end
