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 ! 1 par défaut (2 NON TESTE) 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 ! 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 ! 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 Real :: loggf, Fji Character(len = 50) :: Config_i, Config_j ! Logical :: BLABLA = .FALSE. ! 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é" ! 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) ! 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