Module MOD_F_VAN_REGEMORTER
   !
   USE MPR09_FUNCTIONS, ONLY : INTERP
   !
   Implicit None
   !
   PRIVATE
   PUBLIC :: INTERP_P0_VAN_REGEMORTER, INTERP_P1_VAN_REGEMORTER
   PUBLIC :: GTEST_IN_P_VAN_REGEMORTER, G_VAN_REGEMORTER_ZIRKER
   PUBLIC :: GZ0_IN_P_VAN_REGEMORTER, GZ1_IN_P_VAN_REGEMORTER
   PUBLIC :: Q_VAN_REGEMORTER, OMEGA_VAN_REGEMORTER
   !
   Integer,          Public    :: Z
   Real,             Public    :: TEMP, X0, Fij, gi, gj
   !
   Double Precision, parameter :: c  = 2.99792458D+8		! m/s (Célérité de la lumiére)
   Double Precision, parameter :: a0 = 5.2917720859D-11 	! m   (Rayon de Bohr)
   Double Precision, parameter :: h = 6.62606896D-34		! Js  (Constante de planck)
   Double Precision, parameter :: k_B = 1.3806504D-23		! J/K (Constante de Boltzman) 
   Double Precision, parameter :: m_e = 9.10938215D-31		! kg  (Masse de l'électron)
   Double Precision, parameter :: Ryd = 10973731.568525D0	! 1/m (Energie de Rydberg)
   Double Precision, parameter :: gam = 0.57721566490153286D0
   Double Precision, parameter :: pi = 3.14159265359D0
   !
CONTAINS
   !-----------------------------------------------------------
   Real Function Q_VAN_REGEMORTER( X_PRIME )  
      ! choix = 1 : avec facteur de Gaunt effectif de Van Regemorter (1962)
      ! choix = 2 : avec ajustement de ZIRKER (1962) pour le facteur de Gaunt effectif
      Double precision, intent(in) :: X_PRIME
      !
      Integer, parameter :: choix = 1 ! 1 par défaut (2 NON TESTE)
      Real               :: C_Q
      !
      C_Q = (8 * pi)/SQRT(3.) * (Ryd * h * c / k_B / TEMP)**2 * (1/X0) *  Fij * 1. / (X_PRIME + X0)
      !
      If (Z == 1) Then
         ! ATOME NEUTRE
         Select Case(choix)
            Case(1)
               Q_VAN_REGEMORTER = C_Q * GZ0_IN_P_VAN_REGEMORTER(X_PRIME)
            Case(2) 
               Q_VAN_REGEMORTER = C_Q * G_VAN_REGEMORTER_ZIRKER(X_PRIME)
            Case Default
         End Select
      Else If ( Z > 1) Then
         ! ION POSITIF
         Select Case(choix)
            Case(1)
               Q_VAN_REGEMORTER = C_Q * GZ1_IN_P_VAN_REGEMORTER(X_PRIME)      !
            Case(2)
               Q_VAN_REGEMORTER = C_Q * 0.2D0
            Case Default
         End Select
      Else
         STOP "Q_VAN_REGEMORTER : PB-00" 
      End If
      !  
   End Function Q_VAN_REGEMORTER
   !-----------------------------------------------------------
   Double Precision Function OMEGA_VAN_REGEMORTER( X_PRIME )  
      ! Avec ajustement proposé par Zirker (1962)
      ! pour le facteur de Gaunt effectif
      Double precision, intent(in) :: X_PRIME
      !
      Integer, parameter :: choix = 1 ! (2 NONTESTE)
      Real               :: C_OMEGA
      !
      C_OMEGA = (8 * pi)/SQRT(3.) * (Ryd * h * c / k_B / TEMP) * (1/X0) * gi * Fij
      !
      If (Z == 1) Then
         Select Case(choix)
            Case(1)
               OMEGA_VAN_REGEMORTER = DBLE( C_OMEGA * GZ0_IN_P_VAN_REGEMORTER(X_PRIME) )
            Case(2) 
               OMEGA_VAN_REGEMORTER = DBLE( C_OMEGA * G_VAN_REGEMORTER_ZIRKER(X_PRIME) )
            Case Default
         End Select
      Else If ( Z > 1) Then
         Select Case(choix)
            Case(1)
               OMEGA_VAN_REGEMORTER = DBLE( C_OMEGA * GZ1_IN_P_VAN_REGEMORTER(X_PRIME) ) 
            Case(2)
               OMEGA_VAN_REGEMORTER = DBLE( C_OMEGA ) * 0.2D0
            Case Default
         End Select
      Else
         STOP "Q_VAN_REGEMORTER : PB-00" 
      End If
      ! 
   End Function OMEGA_VAN_REGEMORTER
   !-----------------------------------------------------------
   Real Function INTERP_P0_VAN_REGEMORTER(X0)
      ! Tabulation de la fonction P de Van Regemorter (1962) pour un atome neutre
      Real, intent(in) :: X0
      !
      Integer                        :: I_INF, I_SUP
      Real                           :: X0_INF, X0_SUP
      Real                           :: P0_INF, P0_SUP
      !
      Real, dimension(10), parameter :: VEC_X0 = (/ 0.01, 0.02, 0.04, 0.1, 0.2, 0.4, 1., 2., 4., 10./)
      Real, dimension(10), parameter :: VEC_P0 = (/ 1.160, 0.956, 0.758, 0.493, 0.331, 0.209, 0.100, 0.063, 0.040, 0.023 /)
      !
      Logical, parameter             :: TEST = .FALSE.
      !
      IF (X0 < 0.005 .OR. X0 > 10.) STOP "INTERP_P0_VAN_REGEMORTER : PB-00"
      !
      I_INF = MAXLOC(VEC_X0, dim = 1, mask = VEC_X0 <= X0)
      I_SUP = MINLOC(VEC_X0, dim = 1, mask = VEC_X0 >= X0)
      !
      X0_INF = VEC_X0( I_INF )
      X0_SUP = VEC_X0( I_SUP ) 
      P0_INF = VEC_P0( I_INF )
      P0_SUP = VEC_P0( I_SUP )
      !
      INTERP_P0_VAN_REGEMORTER = INTERP( X0_INF, X0_SUP, P0_INF, P0_SUP, X0)
      !
      IF(TEST) Write(*, *) X0_INF, X0_SUP, P0_INF, P0_SUP, X0, INTERP_P0_VAN_REGEMORTER
      !      
   End Function INTERP_P0_VAN_REGEMORTER
   !-----------------------------------------------------------
   Real Function INTERP_P1_VAN_REGEMORTER(X0)
      ! Tabulation de la fonction P de Van Regemorter (1962) pour un ion positif
      Real, intent(in) :: X0
      !
      Integer                        :: I_INF, I_SUP
      Real                           :: X0_INF, X0_SUP
      Real                           :: P1_INF, P1_SUP
      !
      Real, dimension(10), parameter :: VEC_X0 = (/ 0.01, 0.02, 0.04, 0.1, 0.2, 0.4, 1., 2., 4., 10./)
      Real, dimension(10), parameter :: VEC_P1 = (/ 1.160, 0.977, 0.788, 0.554, 0.403, 0.290, 0.214, 0.201, 0.200, 0.200 /)
      !
      Logical, parameter             :: TEST = .FALSE.
      !
      IF (X0 < 0.005 .OR. X0 > 10.) STOP "INTERP_P1_VAN_REGEMORTER : PB-00"
      !
      I_INF = MAXLOC(VEC_X0, dim = 1, mask = VEC_X0 <= X0)
      I_SUP = MINLOC(VEC_X0, dim = 1, mask = VEC_X0 >= X0)
      !
      X0_INF = VEC_X0( I_INF )
      X0_SUP = VEC_X0( I_SUP ) 
      P1_INF = VEC_P1( I_INF )
      P1_SUP = VEC_P1( I_SUP )
      !
      INTERP_P1_VAN_REGEMORTER = INTERP( X0_INF, X0_SUP, P1_INF, P1_SUP, X0)
      !
      IF(TEST) Write(*, *) X0_INF, X0_SUP, P1_INF, P1_SUP, X0, INTERP_P1_VAN_REGEMORTER
      !      
   End Function INTERP_P1_VAN_REGEMORTER
   !-----------------------------------------------------------
   Double Precision Function GTEST_IN_P_VAN_REGEMORTER(X_PRIME)
      ! Facteur de Gaunt effectif test (plus simple) pour Z > 1
      ! intervenant dans le calcul de P_VAN_REGEMORTER 
      !
      Double Precision, intent(in) :: X_PRIME
      !
      IF (X_PRIME < 0.D0) STOP "GTEST_IN_P_VAN_REGEMORTER : PB-00"
      !
      If ( ( X_PRIME / X0 ) < 2.066 ) Then
         GTEST_IN_P_VAN_REGEMORTER = 0.200D0
      Else
         GTEST_IN_P_VAN_REGEMORTER = SQRT(3.D0) /2.D0/ pi * LOG( X_PRIME / X0 )
      End If
      !Write(*, '(A,2ES10.2)') "X_PRIME / X0, GTEST = ", X_PRIME / X0, GTEST_IN_P_VAN_REGEMORTER
      !
   End Function GTEST_IN_P_VAN_REGEMORTER
   !-----------------------------------------------------------
   Double Precision Function GZ0_IN_P_VAN_REGEMORTER(X_PRIME)
      ! Facteur de Gaunt effectif pour Z = 1
      ! intervenant dans le calcul de P_VAN_REGEMORTER 
      ! FINIR DE TESTER CAR PB A L'INTEGRATION
      !
      Double Precision, intent(in) :: X_PRIME
      !
      Integer                         :: I_INF, I_SUP
      Real                            :: X_VR, X_VR_INF, X_VR_SUP
      Real                            :: GZ0_INF, GZ0_SUP
      !
      Real, dimension (10), parameter :: XZ0 = (/ 0.2, 0.4, 0.6, 0.8, 1., 2., 3., 4., 5., 6. /)
      Real, dimension (10), parameter :: GZ0 = (/ 0.015, 0.034, 0.057, 0.084, 0.124, 0.328, 0.561, 0.775, 0.922, 1.040 /)
      !
      Logical, parameter             :: TEST = .FALSE.
      !
      IF (X_PRIME < 0.D0) STOP "GZ0_IN_P_VAN_REGEMORTER : PB-00" 
      !
      X_VR = SQRT( X_PRIME / X0 )
      !
      If ( X_VR < 0.2 )  Then
         GZ0_IN_P_VAN_REGEMORTER = DBLE( 0.074 * X_VR * (1 + X_VR**2) )
      Else If ( X_VR > 6. ) Then
         GZ0_IN_P_VAN_REGEMORTER = DBLE( SQRT(3.)/2./pi * LOG(X_VR**2) )
      Else
         !
         I_INF = MAXLOC(XZ0, dim = 1, mask = XZ0 <= X_VR)
         I_SUP = MINLOC(XZ0, dim = 1, mask = XZ0 >= X_VR)
         !
         X_VR_INF = XZ0(I_INF) 
         X_VR_SUP = XZ0(I_SUP)
         !
         GZ0_INF = GZ0(I_INF)
         GZ0_SUP = GZ0(I_SUP)
         !  
         GZ0_IN_P_VAN_REGEMORTER = INTERP(X_VR_INF, X_VR_SUP, GZ0_INF, GZ0_SUP, X_VR) * 1.D0
         !
         IF(TEST) Write(*, *) X_VR_INF, X_VR_SUP, GZ0_INF, GZ0_SUP, X_VR, GZ0_IN_P_VAN_REGEMORTER 
         !
      End If
   End Function GZ0_IN_P_VAN_REGEMORTER
   !-----------------------------------------------------------
   Double Precision Function GZ1_IN_P_VAN_REGEMORTER(X_PRIME)
      ! Facteur de Gaunt effectif pour Z > 1
      ! intervenant dans le calcul de P_VAN_REGEMORTER 
      ! FINIR DE TESTER CAR PB A L'INTEGRATION
      !
      Double Precision, intent(in) :: X_PRIME
      !
      Integer                         :: I_INF, I_SUP
      Real                            :: X_VR, X_VR_INF, X_VR_SUP
      Real                            :: GZ1_INF, GZ1_SUP
      Real, dimension (10), parameter :: XZ1 = (/ 0.2, 0.4, 0.6, 0.8, 1., 2., 3., 4., 5., 6. /)
      Real, dimension (10), parameter :: GZ1 = (/ 0.2, 0.2, 0.2, 0.2, 0.2, 0.328, 0.561, 0.775, 0.922, 1.040 /)
      !
      Logical, parameter             :: TEST = .FALSE.
      !
      IF (X_PRIME < 0.D0) STOP "GZ1_IN_P_VAN_REGEMORTER : PB-00" 
      !
      !
      X_VR = SQRT( X_PRIME / X0 )
      !
      If (X_VR < 0.2) Then
         GZ1_IN_P_VAN_REGEMORTER = 0.200D0
      Else If ( X_VR > 6. ) Then
         GZ1_IN_P_VAN_REGEMORTER = DBLE( SQRT(3.)/2./pi * ALOG(X_VR**2) )
      Else
         !
         I_INF = MAXLOC(XZ1, dim = 1, mask = XZ1 <= X_VR)
         I_SUP = MINLOC(XZ1, dim = 1, mask = XZ1 >= X_VR)
         !
         X_VR_INF = XZ1(I_INF) 
         X_VR_SUP = XZ1(I_SUP)
         !
         GZ1_INF = GZ1(I_INF)
         GZ1_SUP = GZ1(I_SUP)
         !  
         GZ1_IN_P_VAN_REGEMORTER = INTERP(X_VR_INF, X_VR_SUP, GZ1_INF, GZ1_SUP, X_VR) * 1.D0
         !
         IF(TEST) Write(*, *) X_VR_INF, X_VR_SUP, GZ1_INF, GZ1_SUP, X_VR, GZ1_IN_P_VAN_REGEMORTER 
         !
      End If
   End Function GZ1_IN_P_VAN_REGEMORTER
   !-----------------------------------------------------------
   Double Precision Function G_VAN_REGEMORTER_ZIRKER(X_PRIME)
      ! POUR ATOME NEUTRE SEULEMENT !
      ! Ajustement proposé par Zirker (1962) pour le facteur de gaunt effectif
      ! vu dans "Spectral line formation", Jefferies, 1968
      ! NON TESTE !
      !
      Double Precision, intent(in) :: X_PRIME
      !
      G_VAN_REGEMORTER_ZIRKER = 0.12D0 * ( (X_PRIME + X0) / X0 - 1.D0 )**0.68D0 
      !
   End Function G_VAN_REGEMORTER_ZIRKER
   !
End Module MOD_F_VAN_REGEMORTER
!
!**************************************************************
!
Module VAN_REGEMORTER
   !
   USE MTD_STRUCT
   USE MPR09_FUNCTIONS
   USE MOD_F_VAN_REGEMORTER, ONLY : INTERP_P0_VAN_REGEMORTER, INTERP_P1_VAN_REGEMORTER, & 
                                    G_VAN_REGEMORTER_ZIRKER, GTEST_IN_P_VAN_REGEMORTER, &
                                    Q_VAN_REGEMORTER, OMEGA_VAN_REGEMORTER,             &
                                    GZ0_IN_P_VAN_REGEMORTER, GZ1_IN_P_VAN_REGEMORTER,   &
                                    TEMP, Fij, gi, gj, X0, Z
   !
   Implicit None
   !
   PRIVATE
   PUBLIC :: VAN_REGEMORTER_BB
   !
   Double Precision, parameter :: c  = 2.99792458D+8		! m/s (Célérité de la lumiére)
   Double Precision, parameter :: a0 = 5.2917720859D-11 	! m   (Rayon de Bohr)
   Double Precision, parameter :: h = 6.62606896D-34		! Js  (Constante de planck)
   Double Precision, parameter :: k_B = 1.3806504D-23		! J/K (Constante de Boltzman) 
   Double Precision, parameter :: m_e = 9.10938215D-31		! kg  (Masse de l'électron)
   Double Precision, parameter :: Ryd = 10973731.568525D0	! 1/m (Energie de Rydberg)
   Double Precision, parameter :: q = 1.60217653D-19            ! Charge élémentaire     A.s
   Double Precision, parameter :: pi = 3.14159265359D0
   Double Precision, parameter :: gam = 0.57721566490153286D0
   !
CONTAINS   
!-----------------------------------------------------------
Subroutine VAN_REGEMORTER_BB(K, TEMP_IN, S0, S1, S2, UPSILON_TEMP_IN)
   !
   Integer,                        intent(in) :: K
   Real,                           intent(in) :: TEMP_IN
   Type(ATOM)      ,               intent(in) :: S0
   Type(LEVELS)    , dimension(:), intent(in) :: S1  !N1
   Type(INDEXLINES), dimension(:), intent(in) :: S2  !N2
   Real,                          intent(out) :: UPSILON_TEMP_IN
   !
   Integer          :: N1
   Integer          :: Ni, Nj 
   Double Precision :: En_i_cm, En_j_cm, DeltaE_cm, En_Ion_cm
   Real             :: loggf, Fji
   Character(len = 50) :: Config_i, Config_j
   !
   Double Precision :: X_PRIME
   !
   Logical          :: BLABLA = .TRUE. 
   !
   Integer :: I
   !TEMPORAIRE
   Real, dimension(16), parameter             :: X_TEMPO =  &
   (/0.001, 0.002, 0.004, 0.01, 0.02, 0.04, 0.1, 0.2, 0.4, 1., 2., 4., 10., 20., 40., 100. /)
   Double Precision, dimension(12), parameter :: X_PRIMEtemp = &
   (/0.12,0.61,1.51,2.83,4.60,6.84,9.62,13.01,17.12,22.15,28.49,37.10/)
   Double Precision, dimension(20), parameter :: X_OVER_X0 = &
   (/0.0001, 0.01, 0.02, 0.04, 0.16, 0.36, 0.64, 1.00, 1.44, 4.00, 9.00, 16.00, 25.00, 36.00, 49.00, 64.00, 81.00, 100.00, 400.00, 900.00/) 
   !
   TEMP = TEMP_IN
   !
   N1 = SIZE(S1)
   Z  = DEGRE_IONISATION( S0%Ion )
   Ni = S2(K)%Ni
   Nj = S2(K)%Nj
   Config_i = S1(Ni)%config
   Config_j = S1(Nj)%config
   !
   En_i_cm   = S1(Ni)%En_cm
   En_j_cm   = S1(Nj)%En_cm    
   !EN_Ion_cm = S1(N1)%En_cm
   DeltaE_cm = En_j_cm - En_i_cm
   !
   gi = S1(Ni)%g
   gj = S1(Nj)%g
   loggf = S2(K)%loggf
   Fij   = 10**( loggf ) / gi
   Fji   = Fij * gi/gj
   !
   X0 = DeltaE_cm * h * c *100. / k_B / TEMP
   !
   If (BLABLA) Call BAVARD( )
   !
   UPSILON_TEMP_IN = UPSILON(TEMP, Z, X0, gi, Fij)
   !
   If (BLABLA) Write(*, *) "UPSILON( TEMP = ", INT(TEMP), " ) = ", UPSILON_TEMP_IN, " sans unité" 
   ! 
   !TEMPORAIRE
   !
   !Write(*, *) "GAUNT VAN REGEMORTER"
   !
   !Do I = 1, 20
   !   Z = 1
   !   Write(*, *) G_VAN_REGEMORTER_ZIRKER(X_OVER_X0(I))
   !End Do
   !
   !Write(*, *) "P_VAN_REGEMORTER"
   !
   !Do I = 1, 16
   !   X0 = X_TEMPO(I)
   !   Z = 1
   !   Write(*, *) P_VAN_REGEMORTER(Z, X0)
   !End Do
   !
   !Write(*, *) "Q_VAN_REGEMORTER"
   !
   !Do I = 1, 12 
   !   Write(*, *) Q_VAN_REGEMORTER(X_PRIMEtemp(I))
   !End Do 
   !
   !Write(*, *) "OMEGA_VAN_REGEMORTER"
   !
   !Do I = 1, 12 
   !   Write(*, *) OMEGA_VAN_REGEMORTER(X_PRIMEtemp(I))
   !End Do
   !  
   CONTAINS
   !-----------------------------------------------------------
   Subroutine BAVARD( )
      !
      Write(*, *) "Degré d'ionisation Z = ", Z
      Write(*, *) TRIM( Config_i ), " => ", TRIM( Config_j )
      Write(*,* ) "gi = ", gi
      Write(*, *) "Fij = ", Fij
      Write(*, *) "DeltaE_eV = ", X0 * k_B * TEMP / q 
      Write(*, *) "X0        = ", X0      
      !
   End Subroutine BAVARD
   !-----------------------------------------------------------
   Real Function P_VAN_REGEMORTER(Z, X0)
      ! ATTENTION : L'UTILISATION DU FACTEUR DE GAUNT EFFECTIF
      !             N'EST PAS ENCORE VALIDE ! 
      !
      USE INTEGRALE, ONLY : GAUSS_LAGUERRE_12PTS, GAUSS_LAGUERRE_8PTS
      !
      Integer, intent(in) :: Z
      Real, intent(in)    :: X0
      !
      Integer, parameter  :: choix  = 1 ! 1 par défaut
      !        
      ! ENERGIE DE LA TRANSITION POSITIVE !
      !
      IF (X0 < 0.) STOP "P_VAN_REGEMORTER : PB-00"
      !
      Select Case(choix)
         !
         Case(0)	! UTILISATION D'UN FACTEUR DE GAUNT EFFECTIF TEST
            !
            If (Z > 1) Then
               P_VAN_REGEMORTER = GAUSS_LAGUERRE_12PTS( GTEST_IN_P_VAN_REGEMORTER )
            Else
               STOP "P_VAN_REGEMORTER : PB-00"
            End If          
            !
         Case(1)	! UTILISATION DU FACTEUR DE GAUNT EFFECTIF
            !
            If ( Z == 1 ) Then 
               ! ATOME NEUTRE   
               P_VAN_REGEMORTER = GAUSS_LAGUERRE_12PTS( GZ0_IN_P_VAN_REGEMORTER )  
            Else If (Z > 1) Then
               ! ION POSITIF
               P_VAN_REGEMORTER = GAUSS_LAGUERRE_12PTS( GZ1_IN_P_VAN_REGEMORTER )
            Else
               Write (*,*) "Pas de méthode prévue par Van Regemorter pour Z < 0 !"
               STOP "P_VAN_REGEMORTER : PB-01"         
            End If
            !
         Case(2)	! INTERPOLATION DANS LA TABLE DE P0 ET P1  
            !
            If (X0 < 0.005) Then
               P_VAN_REGEMORTER = SQRT(3.) /(2. * pi) * ( -gam - ALOG(X0) )           
            Else If ( Z == 1 ) Then 
               ! ATOME NEUTRE 
               If (X0 > 10) Then
                  P_VAN_REGEMORTER = 0.066 / SQRT(X0)
               Else
                  ! INTERPOLATION DANS LA TABLE DE P0
                  P_VAN_REGEMORTER = INTERP_P0_VAN_REGEMORTER(X0)
               End If
            Else If ( Z > 1) Then
               ! ION POSITIF
               If (X0 > 4.) Then
                  P_VAN_REGEMORTER = 0.200
               Else
                  ! INTERPOLATION DANS LA TABLE DE P1
                  P_VAN_REGEMORTER = INTERP_P1_VAN_REGEMORTER(X0)  
               End If
            Else
               Write (*,*) "Pas de méthode prévue par Van Regemorter pour Z < 0 !"
               STOP "P_VAN_REGEMORTER : PB-02" 
            End If              
            !
         Case(3)	! UTILISATION DE LA FORMULE DE ZIRKER POUR FACTEUR DE GAUNT EFFECTIF
            !
            If (Z == 1) Then
               ! ATOME NEUTRE
               P_VAN_REGEMORTER = GAUSS_LAGUERRE_12PTS( G_VAN_REGEMORTER_ZIRKER )
            Else If (Z > 1) Then
               ! ION POSITIF
               P_VAN_REGEMORTER = 0.2
            Else             
               Write (*,*) "Pas de méthode prévue par Van Regemorter pour Z < 0 !"
               STOP "P_VAN_REGEMORTER : PB-03"    
            End If
         !
         Case Default
            Write (*,*) "Choisir une méthode implémenter !"
            STOP "P_VAN_REGEMORTER : PB-04"
      End Select
      !
   End Function P_VAN_REGEMORTER
   !-----------------------------------------------------------
   Real Function UPSILON( TEMP, Z, X0, gi, Fij )
      !
      USE INTEGRALE, ONLY : GAUSS_LAGUERRE_12PTS, GAUSS_LAGUERRE_8PTS
      !
      Integer, intent(in) :: Z
      Real,    intent(in) :: TEMP, X0, gi, Fij 
      !
      Integer, parameter  :: choix = 1
      Real :: C_UPSILON
      !
      C_UPSILON = (8 * pi)/SQRT(3.) * (Ryd * h * c / k_B / TEMP) * (1/X0) * gi * Fij
      !
      Select Case (choix)
         Case(1)
            UPSILON = C_UPSILON * P_VAN_REGEMORTER(Z, X0)
         Case(2)
            UPSILON = GAUSS_LAGUERRE_12PTS( OMEGA_VAN_REGEMORTER )
         Case Default
      End Select
      !
   End Function UPSILON

   !-----------------------------------------------------------
  
End Subroutine VAN_REGEMORTER_BB
   
End Module VAN_REGEMORTER
!
!**************************************************************
!
Program TEST
   !
   ! pathf95 -o TEST_VAN_REGEMORTER mtd_struct.f90 mpr10_external.f90 mpr09_functions.f90 VanRegemorter.f90  
   ! ./TEST_VAN_REGEMORTER
   !
   USE MTD_STRUCT
   USE MPR09_FUNCTIONS
   USE VAN_REGEMORTER
   !
   Implicit None
   !
   Integer                        :: K
   Type(ATOM)                     :: S0
   Type(LEVELS)    , dimension(3) :: S1
   Type(INDEXLINES), dimension(1) :: S2
   Real :: TEMP, UPSILON
   Integer :: choix
   !
   Do choix = 14, 14

   Select Case (choix)
      Case (1)
   S0 = ATOM(20, 'Ca', 'II')
   S1(1) = LEVELS('3p6.4s', 200, 0.D0, 0.D0, 0.D0, 2.0, 0)
   S1(2) = LEVELS('3p6.4p', 211, 25191.510D0, 2.2956193D-1, 3.1234106D0, 2.0, 0)
   S1(3) = LEVELS('Ca III', 999, 95751.870D0, 8.7255524D-1, 1.1871952D1, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.200, 1.336E+08, 8.152E+00, 396.8469)
      Case (2)
   S0 = ATOM(11, 'Na', 'I')
   S1(1) = LEVELS('3s',200, 0.D0, 0.D0, 0.D0, 2.0, 0)
   S1(2) = LEVELS('3p',211, 16967.636D0, 1.5462048D-1, 2.1037186D0, 6.0, 0)
   S1(3) = LEVELS('Na II', 999, 41449.44D0, 0.3777151D0, 5.139075D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, 0.301,0.,0., 588.)
      Case (3)
   S0 = ATOM(12, 'Mg', 'I')
   S1(1) = LEVELS('3p', 311, 21890.854D0, 0.D0, 0.D0, 9, 0)
   S1(2) = LEVELS('4s', 300, 41197.403D0, 0.D0, 0.D0, 3, 0)
   S1(3) = LEVELS('Mg II', 999, 61671.020D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, 0.021, 0., 0., 517.5)  ! 0.021 
      Case (4)
   S0 = ATOM(12, 'Mg', 'II')
   S1(1) = LEVELS('3s', 200, 0.D0, 0.D0, 0.D0, 2, 0)
   S1(2) = LEVELS('3p', 211, 35730.357D0, 0.D0, 0.D0, 6, 0)
   S1(3) = LEVELS('Mg III', 999, 121267.610D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, 0.277, 0., 0., 279.7)    ! 0.260
      Case(5)
   S0 = ATOM(12, 'Mg', 'II')
   S1(1) = LEVELS('4d', 220, 93310.798D0, 0.D0, 0.D0, 10, 0)
   S1(2) = LEVELS('4f', 231, 93799.699D0, 0.D0, 0.D0,14, 0)
   S1(3) = LEVELS('Mg III', 999, 121267.610D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.424, 0., 0., 20454.1) !  -0.242
      Case(6)
   S0 = ATOM(20, 'Ca', 'II')
   S1(1) = LEVELS('3d',220, 13650.190D0, 0.D0, 0.D0, 4, 0)
   S1(2) = LEVELS('4p', 211, 25191.510D0, 0.D0, 0.D0, 2, 0)
   S1(3) = LEVELS('Ca III', 999, 95751.870D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.723, 0., 0., 866.214)
      Case(7)
   S0 = ATOM(20, 'Ca', 'II')
   S1(1) = LEVELS('3d',220, 13650.190D0, 0.D0, 0.D0, 4, 0)
   S1(2) = LEVELS('4p', 211, 25414.400D0, 0.D0, 0.D0, 4, 0)
   S1(3) = LEVELS('Ca III', 999, 95751.870D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -1.416, 0., 0., 849.802)
      Case(8)
   S0 = ATOM(20, 'Ca', 'II')
   S1(1) = LEVELS('3d',220, 13710.880D0, 0.D0, 0.D0, 6, 0)
   S1(2) = LEVELS('4p', 211, 25414.400D0, 0.D0, 0.D0, 4, 0)
   S1(3) = LEVELS('Ca III', 999, 95751.870D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.463, 0., 0., 854.209)
      Case(9)
   S0 = ATOM(20, 'Ca', 'II')!multiplet
   S1(1) = LEVELS('3d',220, 13686.604D0, 0.D0, 0.D0, 10, 0)
   S1(2) = LEVELS('4p', 211, 25340.103D0, 0.D0, 0.D0, 6, 0)
   S1(3) = LEVELS('Ca III', 999, 95751.870D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.243, 0., 0., 857.875) !multiplet
      Case (10)
   S0 = ATOM(12, 'Mg', 'I')
   S1(1) = LEVELS('3p', 311, 21850.405D0, 0.D0, 0.D0, 1, 0)
   S1(2) = LEVELS('4s', 300, 41197.403D0, 0.D0, 0.D0, 3, 0)
   S1(3) = LEVELS('Mg II', 999, 61671.020D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.931, 0., 0., 516.7) ! Aldenius et al. (2007)   
      Case(11)
   S0 = ATOM(12, 'Mg', 'I')
   S1(1) = LEVELS('3p', 311, 21870.464, 0.D0, 0.D0, 3, 0)
   S1(2) = LEVELS('4s', 300, 41197.403D0, 0.D0, 0.D0, 3, 0)
   S1(3) = LEVELS('Mg II', 999, 61671.020D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.450, 0., 0., 517.2) ! Aldenius et al. (2007)   
      Case(12)
   S0 = ATOM(12, 'Mg', 'I')
   S1(1) = LEVELS('3p', 311, 21911.178D0, 0.D0, 0.D0, 5, 0)
   S1(2) = LEVELS('4s', 300, 41197.403D0, 0.D0, 0.D0, 3, 0)
   S1(3) = LEVELS('Mg II', 999, 61671.020D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.239, 0., 0., 518.3) ! Aldenius et al. (2007)
      Case(13)
   S0 = ATOM(12, 'Mg', 'I')
   S1(1) = LEVELS('3p', 311, 21890.854D0, 0.D0, 0.D0, 9, 0)
   S1(2) = LEVELS('4s', 300, 41197.403D0, 0.D0, 0.D0, 3, 0)
   S1(3) = LEVELS('Mg II', 999, 61671.020D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, 0.021, 0., 0., 517.5) !multiplet Aldenius et al. (2007)
      Case(14)
   S0 = ATOM(12, 'Mg', 'I')
   S1(1) = LEVELS('3d', 320, 47957.042D0, 0.D0, 0.D0, 15, 0)
   S1(2) = LEVELS('7f', 331, 59400.763D0, 0.D0, 0.D0, 21, 0)
   S1(3) = LEVELS('Mg II', 999, 61671.020D0, 0.D0, 0.D0, 999, 0)   
   S2(1) = INDEXLINES(2, 1, -0.351, 0., 0., 873.6)!multiplet de 9 composantes (NIST, Hirata, Kurucz) 
      Case Default 
   STOP "PROGRAM TEST : PB-00"
   End Select
   !
   K = 1
   TEMP = 100. !K
   Write(*, *) S0
   !Write(*, *) "T = ", TEMP
   !
   !Do TEMP = 2000, 10000, 2000
      Call VAN_REGEMORTER_BB(K, TEMP, S0, S1, S2, UPSILON)
      Write(*, *) "UPSILON(" , TEMP, " ) = ", UPSILON
   !End Do 
   !
   End Do  
End Program TEST