module MPR09_FUNCTIONS

   use MTD_STRUCT
   use MOD_CSTES, ONLY: c => cste_c, me => cste_m_e, q => cste_q, pi => cste_pi, Ryd => cste_Ryd

   implicit none
   private
   public :: TERME, CONFIG, CHANGE_CONFIG, WINDOW, NQP, NQS, DEGRE_IONISATION, NQP_EFFECTIF
   public :: LAMBDA_AIR
   public :: INTERP, CALCLEVELS, GAUNT_MP, GAUNT_KL, CROSSECT_BF_THRES

   ! PRECISION
   Integer, parameter :: SP = Kind(1.0), DP = Kind(1.d0), P = DP  
   Integer :: I, er, ios

   Integer :: it1 = 100 ! longueur de la chaine de caractère de la configuration électronique
   Logical :: wod = .FALSE.
   !
   Contains
!--------------------------------------------------------------------------------------------------
! CONVERSION DE LONGUEUR D'ONDE DU VIDE A L'AIR
!--------------------------------------------------------------------------------------------------
   Double Precision Function LAMBDA_AIR(x)
   !Conversion des longueurs d'onde du vide dans l'air en A 
   !dans les conditions standards de température et de pression
   !IAU Morton (1991, ApJS, 77, 119)
   !http://www.sdss.org/DR6/products/spectra/vacwavelength.html
   
      Double Precision, intent(in) :: x !Longueur d'onde dans le vide en Angstroem

      LAMBDA_AIR = x / (1.0D0 + 2.735182D-4 + 131.4182D0 / x**2 + 2.76249D8 / x**4)

   End Function LAMBDA_AIR
!--------------------------------------------------------------------------------------------------
! CALCUL DU NOMBRE QUANTIQUE PRINCIPAL EFFECTIF N*
!--------------------------------------------------------------------------------------------------
   Real(SP) Function NQP_EFFECTIF(En_cm, En_Ion_cm, Z) Result(NSTAR)
      Real(DP) :: En_cm, En_Ion_cm
      Integer  :: Z
      Real(DP) :: R = 109678.758
      !
      NSTAR = Z * (R/(En_Ion_cm - En_cm))**0.5
      !
   End Function NQP_EFFECTIF

!--------------------------------------------------------------------------------------------------
! CONVERSION DU DEGRE D'IONISATION EN ENTIER ('I' = 1, 'II' = 2, etc.)
!--------------------------------------------------------------------------------------------------
   Integer Function DEGRE_IONISATION(cIon) Result(IIon)
      Character(len = *), intent(in) :: cIon
      !
      IIon = 0
      !
      If (cIon == 'I') Then
         IIon = 1
      ElseIf (cIon == 'II') Then
         IIon = 2
      ElseIf (cIon == 'III') Then
         IIon = 3  
      Else
         STOP "Degré d'ionisation pas encore implémenté (> 3)."
      End If
      !
   End Function DEGRE_IONISATION
!--------------------------------------------------------------------------------------------------
! CONVERSION DU TERME SPECTROSCOPIQUE (SLP) EN ENTIER = 100*S + 10*L + P (CODE DE LA TOPBASE)
!--------------------------------------------------------------------------------------------------
   function TERME(mot)

      integer          :: TERME
      integer          :: S,L,P
      character(len=*) :: mot
      character(len=1) :: Ltemp

      S = 8                  ! Valeur par défaut retournée si pas modifiée
      L = 8                  ! Valeur par défaut retournée si pas modifiée
      P = 0

      mot = trim(adjustl(mot))

      if (mot=='Limit') then                       ! Ionisation
         TERME = 999
         return
      end if
      if (scan(mot,'()[],/.')/=0.or.mot=='*') then ! Bizzareries
         TERME = 888 ! Bizarrerie 
         return
      endif
      if (len_trim(mot)==1) then
         !Cas 'L'
         Ltemp = mot
      end if
      if (len_trim(mot)==2) then
         !Cas VALD / KURUCZ  'SL'
         if(scan(mot(1:1),'0123456789')/=0) then
            read(mot(1:1),'(I1.1)') S
            Ltemp = mot(2:2)
         else!Autre cas (considéré bizarre)
            S = 8
            Ltemp='Z'
            P = 8
         end if
      end if
      if (len_trim(mot)==3) then
         if (scan(mot(1:1),'0123456789')/=0) then
            !Cas TOPBASE 'SLP'
            read(mot(1:1),'(I1.1)') S
            Ltemp = mot(2:2)
            if(mot(3:3)=='*') P = 1
            !Cas VALD 'SL?'
            if(mot(3:3)=='?') P = 0
         else
            !CAS KURUCZ  'aSL'
            read(mot(2:2),'(I1.1)') S
            Ltemp = mot(3:3)
         end if
      end if
      if (len_trim(mot)>3) then
         if (mot(2:2)==' ') then
            !Cas NIST       'a SL' ou 'a SLP'
            read(mot(3:3),'(I1.1)') S
            Ltemp = mot(4:4)
            if (len_trim(mot)==5) then
               if (mot(5:5)=='*') P = 1
            end if
         else
            !Autres cas (considéré bizarre)
            S = 8
            Ltemp  = 'Z'
            P = 8
         end if
      end if
      if (Ltemp == 'S') L = 0
      if (Ltemp == 'P') L = 1
      if (Ltemp == 'D') L = 2
      if (Ltemp == 'F') L = 3
      if (Ltemp == 'G') L = 4
      if (Ltemp == 'H') L = 5
      if (Ltemp == 'I') L = 6
      if (Ltemp == 'K') L = 7
      if (Ltemp == 'L') L = 8
      if (Ltemp == 'M') L = 9      
      if (scan(Ltemp,'SPDFGHIKLM')==0) L = 8  ! 04/28/2009 Rajout de K (sans doute un oubli)...

      !Calcul du TERME : TERME = S*100 + L*10 + P
      TERME = S*100 + L*10 + P

   end function TERME

!--------------------------------------------------------------------------------------------------
! CONVERSION DE LA CONFIGURATION ELECTRONIQUE TOPBASE  =>  CONFIGURATION ELECTRONIQUE NIST
!--------------------------------------------------------------------------------------------------
!Exemple : 3p6 4s10h  ==>   3p6.4s.10h

   function CONFIG(config_temp)

      integer                          :: I,J,NB,NL,long
      integer,dimension(20)            :: pos_B,pos_L
      character(len=*)                 :: config_temp
      character(len=100)               :: mot,CONFIG
      character(len=100),dimension(20) :: mot_temp,mot_fmt

      !Initialisation
      pos_B = 0
      pos_L = 0
      I = 1
      J = 1
      mot_fmt = ''
      mot_temp = ''
      CONFIG = ''
      config_temp = adjustl(config_temp)

      !Détection du nombre de blancs
   B1:do
         pos_B(I+1) = pos_B(I)+scan(trim(config_temp(pos_B(I)+1:)),' ')
         if ( scan(trim(config_temp(pos_B(I)+1:)),' ')==0) exit
         I = I+1
         if(I>=20) exit B1
      end do B1
      NB = I-1  ! Nombre de blancs dans la config
      pos_B(NB+2) = len_trim(config_temp)+1
      !Remplacement des blancs par des points
   B2:do I = 0,NB
         J = 1
         mot = config_temp(pos_B(I+1)+1:pos_B(I+2)-1)
         !Détermination du nombre de lettres dans un mot
      B3:do
            if (scan(trim(mot(pos_L(J)+1:)),'spdfghiklm')==0) exit B3
            pos_L(J+1) = pos_L(J) + scan(trim(mot(pos_L(J)+1:)),'spdfghiklm')
            J=J+1
            if (pos_L(J+1)==len_trim(mot)) exit B3
            if(J>=20) exit B3
         end do B3
         NL = J-1  ! Nombre de lettre dans le mot
         pos_L(NL+1) = len_trim(mot)
         if (NB == 0 .and. NL == 1) then
            CONFIG = mot
            return
         end if
         !Espacement des lettres par des points
      B4:do J = 1,NL
               mot_temp(J) = trim(mot(pos_L(J)+1:pos_L(J+1)))//'.'
               if (pos_L(J+2)-pos_L(J+1)==4) then  !Cas plus rare   (ex : 2p610f => 2p6.)
                    mot_temp(J) = trim(mot(pos_L(J)+1:pos_L(J+1)+1))//'.'
               end if
               if (pos_L(J+1)-pos_L(J)==4) then     !Cas plus rare (ex : 2p610f => 10f.)
                    mot_temp(J) = trim(mot(pos_L(J)+2:pos_L(J+1)))//'.'
               end if
            mot_fmt(I+1) = trim(mot_fmt(I+1))//trim(mot_temp(J))
         end do B4
      end do B2

      !Ecriture de la nouvelle configuration
      do I = 0,NB
         CONFIG = trim(CONFIG) // trim(mot_fmt(I+1))
      end do
      !Suppression du point de la fin
      long = len_trim(CONFIG)
      if (CONFIG(long:long)=='.') CONFIG = CONFIG(:long-1)

   end function CONFIG
!--------------------------------------------------------------------------------------------------
! ARRANGEMENT DE LA CONFIGURATION ELECTRONIQUE  (POUR IDENTIFICATION VALD => NIST OU TOPBASE)
!--------------------------------------------------------------------------------------------------
   Character(len = 100) Function CHANGE_CONFIG(x, Z) Result(y)
      !
      Character(len = *), intent(in) :: x
      Integer, optional,  intent(in) :: Z
      !
      Integer :: pos1, pos2, N
      Character(len = 100) :: x_temp
      Character(len = 10)  :: Base
      !
      !Exemple : 3d.4s      ==> 4s.3d
      !        : s.3p       ==> s.3p 
      !        : 3s2        ==> 3s2   
      !
      ! Valeur retourné si non modifié
      !
      y = x
      !
      ! Soustraction de la base de la config
      !
      Base = ''
      !
      If (Present(Z)) Then
         If (Z == 20) Then
            Base = '3p6'
         Else If (Z == 12 .OR. Z == 11) Then
            Base = '2p6' 
         Else If (Z == 13) Then
            Base = '3s2'
	    Else 
            STOP "In MPR09_FUNCTIONS, FUNCTION CHANGE_CONFIG: PB-00."
         End If
      End If
      ! 
      N = Len_trim(Base)
      ! 
      If (Index(x, Trim(Base)) /= 0) Then       
         x_temp = x(N + 2 : )
      Else
         x_temp = x
      End If
      ! 
      pos1 = Scan(x_temp, '.')
      !
      If ( pos1 /= 0) Then
         pos2 = Scan(x_temp, ' ')
         If (Scan(x_temp(1 : pos1-1), '12345678890') /= 0) Then
            y = x_temp(pos1 + 1 : pos2 - 1) // '.' // x_temp(1 : pos1 - 1)
         End If
      End If
      !
      IF (Trim(Base) /= '') y = Trim(Base) // '.' // y
      !
    ENd Function CHANGE_CONFIG
!--------------------------------------------------------------------------------------------------
! SELECTIONNE LES SECTIONS EFFICACES DE PHOTOIONISTATION
!--------------------------------------------------------------------------------------------------
   function WINDOW(seff,lambda,N,Nmin,M)
   ! seff   : sections efficaces de photoionisation d'un niveau d'énergie (en cm²)
   ! lambda : longueur d'onde (en Angstroems)
   ! N      : taille du vecteur seff et lambda
   ! Nmin   : nombre de valeurs de seff à ne pas modifier
   ! M      : taille de la fenêtre glissante

      implicit none
      integer                :: j,k,N,M,Nmin
      real,dimension(2000)   :: seff,lambda
      real,dimension(2000)   :: s
      real,dimension(2000,2) :: WINDOW
      real,dimension(:),allocatable :: h,e,r

         allocate(h(M))
         allocate(e(M))
         allocate(r(M))
      !   h = (/(0.5+0.5*cos(j*pi/M),j=1,M)/) ! Fenêtre de Hanning 
      !   h = 1/(sigma*sqrt(2*pi))*(/(exp(-j**2/(2*sigma**2)),j=1,M)/) ! Gaussienne
      !   print*,sum(h)

         do j=1,N
            if (j>=Nmin.and.j<=N-int(M/2)) then
               ! Plage des seff sur lequel on moyenne
               e = seff(j-int(M/2):j+int(M/2))
               ! Pondération des seff à cause de la non régularité des données
               do k = 1,M
                  if (abs(lambda(j)-lambda(j+k-int(M/2)-1))>10e-35) then
                     r(k) = 1./abs(lambda(j)-lambda(j+k-int(M/2)-1))
                  else
                     r(k) = 1 ! valeur au centre ou l'écart est nul.
                  end if
               end do
               !  r = 1./abs(lambda(j)-lambda(j-int(M/2):j+int(M/2)))
               ! Pour un écart nul de lambda, pondération à 1
               ! r((M+1)/2) = 0.5
               r = r/sum(r)  ! Normalisation
               ! Fenêtrage
               s(j) = sum(e*r)
             !  print'(2(ES11.3))',(e(i),r(i),i=1,M)
             !  print*,'Valeur moyenne',s(j),sum(r)
            else
               s(j) = seff(j) ! Aux bords, on garde les valeurs initiales
            end if
         end do
         deallocate(h,e,r)
    !  end do

      WINDOW(:,1) = lambda
      WINDOW(:,2) = s

   end function WINDOW
!--------------------------------------------------------------------------------------------------
! FONCTION D'INTERPOLATION LINEAIRE
!--------------------------------------------------------------------------------------------------
   Function INTERP(x1, x2, y1, y2, x)
     !
     Real, intent(in) :: x1, x2, x
     Real, intent(in) :: y1, y2
     !
     Real :: a, b
     Real :: INTERP
     !
     !Formule d'interpolation linéaire y = a*x + b
     !
     If (x1 /= x2) Then       
        a = (y1 - y2) / (x1 - x2)
     Else
        a = 0.
     End If
     !
     b = y1 - a*x1
     INTERP = a*x + b
     !
   End Function INTERP

!--------------------------------------------------------------------------------------------------
! FONCTION DONNANT LE NOMBRE QUANTIQUE PRINCIPAL DE L'ETAT DE L'ELECTRON EXCITE
!--------------------------------------------------------------------------------------------------
   Integer Function NQP(x)
      ! x     : chaine de caractère contenant la configuration électronique
      ! NQP   : Nombre quantique Principal
      !
      Character(len = *), intent(in) :: x 
      Character(len = it1)           :: conf   ! chaine de caractère contenant la configuration électronique
      Integer                        :: pos    ! Position du dernier point
      Integer                        :: taille ! Longueur de la configuration électronique
      Logical                        :: BLABLA = .FALSE.
      !
      ! Valeur retournée si pas modifiée
      !
      NQP = -2
      !
      conf   = Trim(x)
      taille = Len_trim(conf)
      !
      pos = Scan(conf, '.', back = .true.)
      !
      If (Index(conf, 'II') /= 0) Then
         NQP = -1
         If (wod) Then
            Write(*, *) conf       
            Write(*, *) "Détection d'une ionisation"
            Read(*, *)
         End If
         RETURN
      End If
      !
      If (Scan(conf, '()<>/') /= 0) Then
         NQP = -1
         If (BLABLA) Then         
            Write(*, *) conf  
            Write(*, *) "Attention, problème pour la détection du NQP."
         End If
         RETURN
      End If
      !
      IF(BLABLA) Write(*, *) 'Config : ',trim(conf),' Taille : ',taille,' Position du point : ', pos
      !
      ! Cas où config :3p2
      If (pos == 0) Read(conf(1 : 1),'(I2)') NQP      
      !
      If (Scan(conf(taille-1 : taille-1), '0123456789') == 0) Then
         ! Cas où Config : 3p6.3d2 
         Read(conf(pos+1 : taille-2),'(I2)') NQP
      else
         ! Cas général où  Config : 3p6.4s.10g, 3p6.4s.5g, ... 
         Read(conf(pos+1:taille-1),'(I2)') NQP
      end if
      !
      IF(BLABLA) Write(*, *) "NQP = ", NQP
      !
   End Function NQP   
   !
   !--------------------------------------------------------------------------------------------------
   ! FONCTION DONNANT LE NOMBRE QUANTIQUE SECONDAIRE DE L'ETAT DE L'ELECTRON EXCITE
   !--------------------------------------------------------------------------------------------------
   !
   Integer Function NQS(x)
      ! conf  : chaine de caractère contenant la configuration électronique
      ! NQS   : Nombre Quantique Secondaire
      !
      Character(len = *), intent(in) :: x 
      Character(len = it1)           :: conf
      Character(len = 1)             :: cy
      Integer                        :: pos1   ! Position du dernier point dans x
      Integer                        :: pos2   ! Position de l'orbitale l dans x
      Integer                        :: taille ! Longueur de la configuration électronique
      Logical                        :: BLABLA = .FALSE.
      !
      ! Valeur retournée si pas modifiée
      !
      NQS = -2
      !
      conf = trim(x)
      taille = len_trim(conf)
      !
      pos1 = Scan(conf, '.', back = .true.)
      !Write(*,*)"conf, pos1, pos2 :", conf, pos1, pos2
      !
      ! Cas où Config : 3p2
      If (pos1 == 0) Then
         Read(conf(2 : 2),'(A)') cy
         pos1 = 1
      End If
      !
      pos2 = Scan(conf(pos1 : taille), 'spdfghiklm')
      !
      If (Index(conf, 'II') /= 0) Then
         NQS = -1
         RETURN
      End If
      !
      If (Scan(conf, '()<>/') /= 0) Then
         NQS = -1
         If (BLABLA) Then
            Write(*, *) "Attention, problème pour la détection du NQS."
            Read(*, *)
         End If 
         RETURN
      End If
      !
      IF(BLABLA) Write(*, *) 'Config : ', Trim(conf), ' Taille : ', taille, ' Position de l orbitale : ', pos2
      !
      If (Scan(conf(taille-1 : taille-1), '0123456789') == 0) Then
         ! Cas où Config : 3p6.3d2 
         Read(conf(pos1+pos2-1 : pos1+pos2-1), '(A)') cy
      Else
         ! Cas général où  Config : 3p6.4s.10g, 3p6.4s.5g, ... 
         Read(conf(pos1+pos2-1 : pos1+pos2-1), '(A)') cy

      End If
      !
      IF(BLABLA) Write(*, *) "cy = ", cy
      !
      IF (cy == 's') NQS = 0
      IF (cy == 'p') NQS = 1
      IF (cy == 'd') NQS = 2
      IF (cy == 'f') NQS = 3
      IF (cy == 'g') NQS = 4
      IF (cy == 'h') NQS = 5
      IF (cy == 'i') NQS = 6
      IF (cy == 'k') NQS = 7
      IF (cy == 'l') NQS = 8
      IF (cy == 'm') NQS = 9
      !
      IF(BLABLA) Write(*, *) "NQS = ", NQS
      !
   End Function NQS
   !
!--------------------------------------------------------------------------------------------------
! FONCTION DONNANT LA SECTION EFFICACE AU SEUIL (MENZEL ET PICKERIS 1935)
!--------------------------------------------------------------------------------------------------
!   function SEFFSEUIL(lambda,n,Z)
!   ! lambda : énergie d'ionisation du niveau exprimée en A
!   ! n      : nombre quantique principal du terme spectro
!   ! Z      : numéro atomique
! 
!      real    :: lambda
!      integer :: n, Z
!      real    :: SEFFSEUIL
!      real    :: Gaunt
!      double precision :: a, nu, Ryd_Hz
!
!      Ryd_Hz = Ryd*100*c
!      print'(A,ES17.11)','Ryd_Hz = ',Ryd_Hz
!
!      nu = c/(lambda*1D-10)  
!      a = nu/(Ryd_Hz * Z**2)
!     !   Gaunt = 1 - 0.1728 * a**(1/3) * (2/(a*n**2) - 1)
!      Gaunt = 0.0 ! TEMPORAIRE
!      SEFFSEUIL = 7.91 * n * Gaunt * 1E-18 ! cm2
!      print*,"Fréquence d'ionisation : ",nu
!      print*,'a = ',a
!      print'(A,F6.3)',"Gaunt = ",Gaunt
!   end function SEFFSEUIL

!--------------------------------------------------------------------------------------------------
! DETERMINATION DES COEFFICIENTS DE DE-EXCITATION RADIATIF DES NIVEAUX SELECTIONNES
!--------------------------------------------------------------------------------------------------
   subroutine CALCLEVELS(N1,N2,l5,S1)
   ! S1     : structure des transitions indexées sélectionnées
   ! N1     : nombre de niveaux d'énergie sélectionnés
   ! N2     : nombre de transitions sélectionnés

      integer,                      intent(in)    :: N1,N2
      logical,                      intent(inout) :: l5
      type(INDEXLINES),dimension(:),intent(inout) :: S1

      integer               :: I,J,L
      real,dimension(N1,N1) :: Aij
      real                  :: gammarad
      character             :: c

      !Initialisation de la matrice
      Aij = 0.

      !Affectation de la matrice des coefficients de de-excitation radiatif Aij
      do L = 1,N2
         I = S1(L)%Ni
         J = S1(L)%Nj
         Aij(I,J) = S1(L)%Aij
      end do

      if (l5) then
          print*,'Calcul de la probabilité de dé-éxcitation radiative Gamma rad :'
          print '(/,A,/)','-------------------------------------------------------------------'
      end if

      !Calcul du gamma radiatif
      do L = 1,N2
         I = S1(L)%Ni
         J = S1(L)%Nj

         gammarad = log10(sum(Aij(1:I-1,I)) + sum(Aij(1:J-1,J)))

         if (l5)then
            print*,'Théo : ',S1(L)%Gr,' Calculé : ',gammarad,'  Encore (o/n) ? '
            read(*,*) c
            if (c=='n') l5 = .false.
         end if

         !Affectation du gamma radiatf dans la structure des niveaux indexés
         S1(L)%Gr = gammarad
      end do

   end subroutine CALCLEVELS
   !
   !--------------------------------------------------------------------
   ! FONCTION CALCULANT LE FACTEUR DE GAUNT BOUND-FREE  
   ! (MENZEL ET PICKERIS 1935)
   !--------------------------------------------------------------------
   !
   Real(SP) Function GAUNT_MP(Z, nu, n) Result(y)
      !
      Integer, intent(in) :: Z  ! Dégré d'ionisation
      Real(P), intent(in) :: nu ! Fréquence d'ionisation au seuil (Hz)
      Integer, intent(in) :: n  ! Nombre quantique principal
      !
      Real(P)             :: Ryd_Hz
      Real(P)             :: a
      Character(len = 1)  :: choix
      !
      Ryd_Hz = Ryd * c
      !
      a = nu / (Ryd_Hz * Z**2)
      y = 1 - 0.1728 * a**(1/3) * ( 2 / (n**2 * a) - 1)
      !
      Write(*, *) Z, c, Ryd_Hz, a, y
      !
      If (y < 0) then
         Write(*, *) "Facteur de Gaunt négatif avec le modèle de Menzel et Pickeris !!! &
                     &Continuer (utilisation du modèle de Karzas & Latter) ? (o/n)"
         !
B1 :     Do 
            Read(*, '(A)') choix
            Select Case (choix)
               Case ('o', 'O', 'y', 'Y')
                    y = 0.
                    EXIT B1
               Case ('n', 'N')
                    STOP 
               Case Default
                 CYCLE B1
            End Select
         End Do B1
         !
      End if
      !
   End Function GAUNT_MP     
   !
   !--------------------------------------------------------------------
   ! FONCTION CALCULANT LE FACTEUR DE GAUNT BOUND-FREE 
   ! (KARZAS et LATTER ApJS, 1961)
   !--------------------------------------------------------------------
   !
   Real(SP) Function GAUNT_KL(Z, nu, n, l) Result(y)
      !
      Integer, intent(in) :: Z  ! Dégré d'ionisation
      Real(P), intent(in) :: nu ! Fréquence d'ionisation au seuil
      Integer, intent(in) :: n  ! Nombre quantique principal
      Integer, intent(in) :: l  ! Nombre quantique secondaire
      !
      Real(P)            :: E      ! énergie d'ionisation cm-1
      Real(P)            :: rho, etha
      Real(P)            :: G11, G12, G21, G22 ! 
      Real(P)            :: Bra1, Bra2
      Real(P)            :: Sigma1, Sigma2, SigmaG ! Section efficace de Gaunt
      Real(P)            :: SigmaK                 ! Section efficace de Kramer
      Real(P)            :: a       ! constante de section efficace
      Real(P)            :: s1, s11, s111, s2, s22, s222
      !
      ! Sigma1 = f(a, s1(Bra1), s11, s111(G11, G12))
      ! Sigma2 = f(a, s2(Bra2), s22, s222(G21, G22))
      ! SigmaG = Sigma1 + Sigma2 ! Section efficace de Gaunt
      ! y = SigmaG / SigmaK      ! Facteur de Gaunt
      ! 
      IF (l >= n) STOP "IN MPR09_FUNCTIONS, FUNCTION GAUNT_KL: l MUST BE STRICTLY LESS THAN n !!!"
      IF ( Z < 0   ) STOP "IN MPR09_FUNCTIONS, FUNCTION GAUNT_KL: PB-00."
      IF (nu < 0.d0) STOP "IN MPR09_FUNCTIONS, FUNCTION GAUNT_KL: PB-01."
      IF (n  < 0   ) STOP "IN MPR09_FUNCTIONS, FUNCTION GAUNT_KL: PB-02."
      IF (l  < 0   ) STOP "IN MPR09_FUNCTIONS, FUNCTION GAUNT_KL: PB-03."
      !
      E      = nu / (100*c)
      etha   = (Z**2 * Ryd/100.D0 / E)**0.5       
      rho    = etha / n 
      a      = q**2 / (me * c * nu)
      !
      ! Calcul de Sigma(nl => E, l-1)
      !
      G11  = POLY_G(-(l + 1 - n), l, etha, rho)
      G12  = POLY_G(-(l - 1 - n), l, etha, rho)
      Bra1 = BRACES(l-1, etha)
      !
      s1 = 0.0
      !
      If ( l > 0 ) Then
         s1 = (2.0_P**(4*l)/3.0_P) * l**2 * fact(n+l) * Bra1 / &
              ( fact(2*l+1) * fact(2*l-1) * fact(n-l-1) )
      End If
      !
      s11  = dexp(-4*etha * (pi/2.0_P - atan(rho))) / (1 - dexp(-2.0_P*pi*etha)) * &
             rho**(2*l+2) / (1 + rho**2)**(2*n-2)
      !
      s111 = (G11 - (1 + rho**2)**(-2) * G12)**2 
      !
      IF (s1 < 0. .OR. s11 < 0. .OR. s111 < 0.) STOP "IN MPR09_FUNCTIONS, FUNCTION GAUNT_KL: PB-04."
      !
      Sigma1 = s1 * s11 * s111 * pi * a
      !
      ! Calcul de Sigma(nl => E, l-1)
      ! 
      G21 = POLY_G(-(l + 1 - n), l+1, etha, rho)
      G22 = POLY_G(-(l-n), l+1, etha, rho)
      Bra2 = BRACES(l+1, etha)
      !
      s2   = (2.0_P**(4*l+6)/3.0_P) * (l+1)**2 * fact(n+l) * Bra2 / & 
             ((2*l+1) * fact(2*l+1) * fact(2*l+2) * fact(n-l-1) * ((l+1)**2 + etha**2)**2)
      !
      s22  = dexp(-4*etha *(pi/2.0_P - atan(rho))) / (1 - dexp(-2*pi*etha)) * &
             rho**(2*l+4) * etha**2 / (1 + rho**2)**(2*n)
      !
      s222 = ( (l+1-n) * G21 + (l+1+n)/(1 + rho**2) * G22 )**2
      !
      IF (s2 < 0. .OR. s22 < 0. .OR. s222 < 0.) STOP "IN MPR09_FUNCTIONS, FUNCTION GAUNT_KL: PB-05." 
      !
      Sigma2 = s2 * s22 * s222 * pi * a  
      !
      ! Calcul de Sigma Gaunt
      !
      SigmaG = Sigma1 + Sigma2
      !
      ! Calcul de Sigma Kramer
      !   
      SigmaK = 2**4/(3.0_P * 3**0.5) * (1.0_P/n) * (rho**2/(1+rho**2))**2 * a 
      !
      ! Calcul du facteur de Gaunt lié - libre
      !
      y = SigmaG / SigmaK
      !
   End Function GAUNT_KL
   !
   !--------------------------------------------------------------------
   ! FONCTION CALCULANT LA SOLUTION DE 
   ! L'EQUATION DIFFERENTIELLE HYPERGEOMETRIQUE (POUR GAUNT_KL)
   !-------------------------------------------------------------------- 
   !
   Real(P) Function POLY_G(m, ll, etha, rho) Result(G)
      Integer, Intent(in)  :: m, ll
      Real(P), Intent(in)  :: etha, rho
      !
      Integer                            :: Ncoeff    ! Nombre de coefficient du polynome G
      Integer                            :: s         ! Indice de comptage
      Real(P), Dimension(:), Allocatable :: b         ! Coefficient du polynome G
      !  
      G = 0.0_P
      Ncoeff = 2*m + 1
      !
      IF (ll == 0) RETURN
      !      
      Allocate(b(Ncoeff), stat = er)
      IF (er /= 0) STOP "IN MPR09_FUNCTIONS, FUNCTION POLY_G: PB-00."
      !
      b    = 0
      b(1) = 1.0_P
      If (m /= 0) b(2) = 2 * m * etha / ll
      !    
      Do I = 3, Ncoeff
         s = I - 1
         b(I) = -(1.0_P/(s * (s + 2*ll -1))) * &
                 ( 4.0_P*etha * (s - 1 - m) * b(s) + (2*m + 2 - s) * (2*m + 2*ll + 1 - s) * b(s-1) ) 
      End Do
      !
      Do I = 1, Ncoeff
         G = G + b(I) * rho**(I-1)
      End Do 
      !
      Deallocate(b) 
      !
   End Function POLY_G  
   !
   !--------------------------------------------------------------------
   ! FONCTION CALCULANT UNE PARTIE DE LA SECTION EFFICACE LIE LIBRE 
   ! (POUR GAUNT_KL)
   !--------------------------------------------------------------------
   !
   Real(P) Function BRACES(lll, etha) Result(y)
      !
      Integer, Intent(in)  :: lll
      Real(P), Intent(in)  :: etha
      !
      y = 1.0_P
      !
      IF (lll == 0) RETURN
      !     
      Do I = 1, lll
         y = y * (I**2 + etha**2)
      End Do    
      !
   End Function BRACES    
   !
   !--------------------------------------------------------------------
   ! FONCTION FACTORIELLE
   !--------------------------------------------------------------------
   !
   Real(P) Function FACT(N) Result(Y)
      !
      Integer, Intent(in)  :: N
      Integer              :: I
      !
      Y = 1
      !
      IF (N == 0) RETURN
      !
      IF (N < 0) Then
       STOP "IN MPR09_FUNCTIONS, FUNCTION FACT: PB-00"
      END IF
      !
      Do I = N, 1, -1
         Y = Y * I  
      End Do
      !
   End Function FACT
   !
   !--------------------------------------------------------------------
   ! FONCTION CALCULANT LA SECTION EFF. DE PHOTOIONISATION AU SEUIL NU
   !--------------------------------------------------------------------
   !
   Real(SP) Function CROSSECT_BF_THRES(g_bf, Z, nu, n) Result(y)
      !
      Real(SP), intent(in) :: g_bf ! Facteur de Gaunt bound-free
      Integer, intent(in) :: Z     ! Degré d'ionisation de l'atome (1 pr neutre, 2 ionisé 1 fois...)
      Real(P), intent(in) :: nu    ! Fréquence d'ionisation au seuil
      Integer, intent(in) :: n     ! Nombre quantique principal
      !
      Real(P)             :: E, rho ! Energies
      Real(P)             :: a      ! constante de section efficace
      Real(P)             :: SigmaK ! Section efficace de Kramer
      !
      Integer, parameter  :: choix = 1       
      !
      Select Case (choix)
        Case (1)
          ! Mihalas (Stellar Atmospheres, 2nd edition, p 99, 1978)
          y = 2.815D29/(n**5 * nu**3) * g_bf * Z**4
          Write(*, *) "Mihalas     : ", y*1D18, ' Mbarn'
        Case (2)
          ! Karzas et Latter (ApJS, 1961)
          a     = q**2 / (me * c * nu)
          E     = nu / c
          rho   = (Z**2 * Ryd / 100.D0 / E)**0.5 / n       
          SigmaK = 2**4/(3.0_P * 3**0.5) * (1.0_P/n) * (rho**2/(1+rho**2))**2 * a 
          y =  g_bf * SigmaK
          Write(*, *) "Karzas & Latter : ", y*1D18, ' Mbarn'
        Case Default
      End Select
      !   
   End Function CROSSECT_BF_THRES
   !
end module MPR09_FUNCTIONS