module MPR09_FUNCTIONS

   use MTD_STRUCT

   implicit none
   private
   public :: TERME, 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 = .TRUE.
   include 'cstes.inc'

   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) :: Ryd = 109678.758
      !
      NSTAR = Z * (Ryd/(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 (scan(Ltemp,'SPDFGHIK')==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(1) = 0
      pos_L(1) = 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:)),'spdfghk')==0) exit B3
            pos_L(J+1) = pos_L(J) + scan(trim(mot(pos_L(J)+1:)),'spdfghk')
            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

!--------------------------------------------------------------------------------------------------
! 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                :: i,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) Result(y)
    ! 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

      ! Valeur retournée si pas modifiée
      y = -2

      conf   = Trim(x)
      taille = Len_trim(conf)

      !Write(*, *)

      pos = Scan(conf, '.', back = .true.)

      If (Index(conf, 'II') /= 0) Then
         y = -1
         If (wod) Then
            Write(*, *) conf       
            Write(*, *) "Détection d'une ionisation"
            Read(*, *)
         End If
         RETURN
      End If

      If (Scan(conf, '()<>/') /= 0) Then
         y = -1
         If (wod) Then         
            Write(*, *) conf  
            Write(*, *) "Attention, problème pour la détection du NQP."
         End If
         RETURN
      End If

      !Write(*, *) 'Config : ',trim(conf),' Taille : ',taille,' Position du point : ', pos

      ! Cas où config :3p2
      If (pos == 0) Read(conf(1 : 1),'(I2)') y      

      If (Scan(conf(taille-1 : taille-1), '0123456789') == 0) Then
      ! cas où Config : 3p6.3d2 
         Read(conf(pos+1 : taille-2),'(I2)') y
      else
      ! Cas général où  Config : 3p6.4s.10g, 3p6.4s.5g, ... 
         Read(conf(pos+1:taille-1),'(I2)') y
      end if

   end function NQP   

!--------------------------------------------------------------------------------------------------
! FONCTION DONNANT LE NOMBRE QUANTIQUE SECONDAIRE DE L'ETAT DE L'ELECTRON EXCITE
!--------------------------------------------------------------------------------------------------


  Integer Function NQS(x) Result(y)
    ! 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

      ! Valeur retournée si pas modifiée
      y = -2

      conf = trim(x)
      conf = trim(conf)
      taille = len_trim(conf)

      pos1 = Scan(conf, '.', back = .true.)
      pos2 = Scan(conf(pos1 : taille), 'spdfghik')

      If (Index(conf, 'II') /= 0) Then
         y = -1
         RETURN
      End If

      If (Scan(conf, '()<>/') /= 0) Then
         y = -1
         If (wod) Then
            Write(*, *) "Attention, problème pour la détection du NQS."
            Read(*, *)
         End If 
         RETURN
      End If

     ! Write(*, *) 'Config : ', trim(conf), ' Taille : ', taille, ' Position de l orbitale : ', pos2

      ! Cas où Config : 3p2
      IF (pos1 == 0) Read(conf(2 : 2),'(A)') cy

      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
         !Write(*, *) cy
      Else
      ! Cas général où  Config : 3p6.4s.10g, 3p6.4s.5g, ... 
         Read(conf(pos1+pos2-1 : pos1+pos2-1), '(A)') cy
         !Write(*, *) cy
      End If

      If (cy == 's') y = 0
      If (cy == 'p') y = 1
      If (cy == 'd') y = 2
      If (cy == 'f') y = 3
      If (cy == 'g') y = 4
      If (cy == 'h') y = 5
      If (cy == 'i') y = 6
      If (cy == 'k') y = 7

   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 = 100 * 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 !!! Continuer ? (o/n)"
         Read(*, '(A)') choix

         Select Case (choix)
            Case ('o', 'O', 'y', 'Y')
                 y = 0.
            Case Default
                 STOP
         End Select

      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 "l doit être inférieur à n !!!"
      IF (Z < 0 .OR. nu < 0 .OR. n < 0 .OR. l < 0) STOP "Z ou nu ou n ou l est négatif !!!"

      E      = nu / (100*c)
      etha   = (Z**2 * Ryd / 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 = (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) )

      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 "Section efficace négative ! (1)"

      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 "Section efficace négative ! (2)" 
   
      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

      Allocate(b(Ncoeff), stat = er)
      IF (er /= 0) STOP "Problème d'allocation de b."
      
      b(1) = 1.0_P
      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) STOP "Factorielle d'un nombre négatif !"

      Do I = N, 1, -1
         Y = Y * I  
      End Do
    
   End Function FACT

!--------------------------------------------------------------------------------------------------
! FONCTION CALCULANT LA SECTION EFFICACE 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
       

      ! Mihalas (Stellar Atmospheres, 2nd edition, p 99, 1978)
       y = 2.815D29/(n**5 * nu**3) * g_bf * Z**4

       Write(*, *) "Mihalas     : ", y*1D18, ' Mbarn'

      ! Karzas et Latter (ApJS, 1961)
      ! a     = q**2 / (me * c * nu)
      ! E     = nu / c
      ! rho   = (Z**2 * Ryd / 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

   End Function CROSSECT_BF_THRES

end module MPR09_FUNCTIONS