Module MOD_F_SEATON ! Implicit None ! PRIVATE PUBLIC :: F_SEATON, Q_SEATON, QX_SEATON, OMEGA_SEATON PUBLIC :: FCT_Q0, FCT_Q1, OMEGA_0, OMEGA_1 PUBLIC :: SEATON_PHI, SEATON_ZETA, BETA_0, BETA_1 ! Real, Public :: C_SEATON Real, Public :: TEMP, R0, Fij, gi, gj Double Precision, Public :: X0 Integer, Save, Public :: ON_DEPASSE = 0 ! 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 :: pi = 3.14159265359D0 Double Precision, parameter :: q = 1.60217653D-19 ! Charge élémentaire A.s ! CONTAINS !---------------------------------------------------------------------- Double Precision Function F_SEATON ( BETA1 ) Result(y) ! USE BESSEL, BESSEL_K0 => dbesk0, BESSEL_K1 => dbesk1 ! Double Precision, intent(in) :: BETA1 ! y = BESSEL_K0( BETA1 )**2 + BESSEL_K1( BETA1 )**2 - C_SEATON ! End Function F_SEATON !---------------------------------------------------------------------- Double Precision Function OMEGA_SEATON ( X_PRIME ) ! Double Precision, intent (in) :: X_PRIME ! Ej/kT ! Real :: OM_0, OM_1 ! Logical :: BLABLA = .FALSE. ! OM_0 = OMEGA_0(TEMP, X_PRIME, X0, R0, Fij) ! Tel que définit dans Proc. Phys. Soc., Seaton, 1962 OM_1 = OMEGA_1(TEMP, X_PRIME, X0, Fij) ! Tel que définit dans Proc. Phys. Soc., Seaton, 1962 If (BLABLA) Write(*, *) OM_0, OM_1 ! If ( OM_0 >= OM_1 ) Then OMEGA_SEATON = OM_0 !IF (BLABLA) Write(*, *) "OM_0(", X_PRIME, ") = ", OM_0 Else OMEGA_SEATON = OM_1 !IF (BLABLA) Write(*, *) "OM_1(", X_PRIME, ") = ", OM_1 End If ! End Function OMEGA_SEATON !---------------------------------------------------------------------- Double Precision Function Q_SEATON ( X_PRIME ) ! Double Precision, intent (in) :: X_PRIME ! Ej/kT ! Real :: Q0, Q1 ! Logical :: BLABLA = .FALSE. ! Q0 = FCT_Q0(TEMP, X_PRIME, X0, R0, Fij) ! Tel que définit dans Proc. Phys. Soc., Seaton, 1962 Q1 = FCT_Q1(TEMP, X_PRIME, X0, Fij) ! Tel que définit dans Proc. Phys. Soc., Seaton, 1962 If (BLABLA) Write(*, *) Q0, Q1 ! If ( Q0 <= Q1 ) Then Q_SEATON = Q0 !IF (BLABLA) Write(*, *) "Q0(", X_PRIME, ") = ", Q0 Else Q_SEATON = Q1 !IF (BLABLA) Write(*, *) "Q1(", X_PRIME, ") = ", Q1 End If ! End Function Q_SEATON !---------------------------------------------------------------------- Double Precision Function QX_SEATON ( X_PRIME ) ! Fonction à utiliser pour l'intégration en Gauss-Laguerre ! Double Precision, intent (in) :: X_PRIME ! Ej/kT ! QX_SEATON = (X_PRIME + X0) * Q_SEATON( X_PRIME ) ! End Function QX_SEATON !---------------------------------------------------------------------- Real Function OMEGA_0(TEMP, X_PRIME, X0, R0, Fij) ! Real, intent(in) :: TEMP, R0, Fij Double Precision, intent(in) :: X_PRIME, X0 ! Real :: BETA0, CSTE ! OMEGA_0 = 0. ! IF ( X_PRIME < 0.D0 ) RETURN ! BETA0 = BETA_0(TEMP, X_PRIME, X0, R0) CSTE = 8. * (Ryd* h * c / k_B / TEMP) * gi * Fij / X0 OMEGA_0 = CSTE * SEATON_PHI( DBLE(BETA0) ) ! End Function OMEGA_0 !---------------------------------------------------------------------- Real Function FCT_Q0(TEMP, X_PRIME, X0, R0, Fij) Result(y) ! Real, intent(in) :: TEMP, R0, Fij Double Precision, intent(in) :: X_PRIME, X0 ! Real :: BETA0, C_Q0 ! y = 0. ! IF ( X_PRIME < 0.D0 ) RETURN ! BETA0 = BETA_0(TEMP, X_PRIME, X0, R0) C_Q0 = 8. * (Ryd* h * c / k_B / TEMP)**2 * Fij / X0 y = C_Q0 * SEATON_PHI( DBLE(BETA0) ) / ( X_PRIME + X0) ! End Function FCT_Q0 !---------------------------------------------------------------------- Real Function OMEGA_1(TEMP, X_PRIME, X0, Fij) !Attention il faut que X_PRIME, Fij, gmin représente l'état de !plus faible poids statistique ! Real, intent(in) :: TEMP, Fij Double Precision, intent(in) :: X_PRIME, X0 ! Real :: BETA1, Fji, CSTE, g_min, Fo ! OMEGA_1 = 0. !gedit IF (X_PRIME < 0.D0) RETURN ! Fji = Fij * gi / gj ! If (gi <= gj) Then g_min = gi Fo = Fij Else g_min = gj Fo = Fji End If ! BETA1 = BETA_1(TEMP, X_PRIME, X0, Fo) CSTE = 8. * (Ryd* h * c / k_B / TEMP) * g_min * Fo / X0 OMEGA_1 = CSTE * ( 0.5 * SEATON_ZETA( DBLE(BETA1) ) + SEATON_PHI( DBLE(BETA1) ) ) / ( X_PRIME + X0 ) ! End Function OMEGA_1 !---------------------------------------------------------------------- Real Function FCT_Q1(TEMP, X_PRIME, X0, Fij) Result (y) ! Real, intent(in) :: TEMP, Fij Double Precision, intent(in) :: X_PRIME, X0 ! Real :: BETA1, C_Q1, Fji, g_min, Fo ! y = 0. ! IF (X_PRIME < 0.D0) RETURN ! Fji = Fij * gi / gj ! If (gi <= gj) Then g_min = gi Fo = Fij Else g_min = gj Fo = Fji End If ! BETA1 = BETA_1(TEMP, X_PRIME, X0, Fo) C_Q1 = 8. * (Ryd* h * c / k_B / TEMP)**2 * Fo / X0 y = C_Q1 * ( 0.5 * SEATON_ZETA( DBLE(BETA1) ) + SEATON_PHI( DBLE(BETA1) ) ) / ( X_PRIME + X0 ) ! End Function FCT_Q1 !---------------------------------------------------------------------- Real Function BETA_0(TEMP, X_PRIME, X0, R0) Result(BETA) Real, intent(in) :: TEMP, R0 Double Precision, intent(in) :: X_PRIME, X0 ! X_PRIME = Ej/kT ! IF (X_PRIME < 0.) STOP "FUNCTION BETA_0 : X_PRIME < 0 !!!" BETA = R0 * (k_B *TEMP * ( X_PRIME + X0 ) / Ryd / h / c)**0.5 * X0 / (2.*X_PRIME + X0) ! End Function BETA_0 !---------------------------------------------------------------------- Real Function BETA_1(TEMP, X_PRIME, X0, Fo) ! USE ZERO, ONLY : BRENT ! Real, intent(in) :: TEMP, Fo Double Precision, intent(in) :: X_PRIME, X0 ! Double Precision, parameter :: prec = 1.0d-12 Double Precision :: borne_inf, borne_sup, BETA1 ! !borne_inf = SQRT( TINY(borne_inf) ) !borne_sup = SQRT( HUGE(borne_sup) ) borne_inf = 1.D-10 borne_sup = 1000.D0 ! IF (X_PRIME < 0.) STOP "FUNCTION BETA_1 : X_PRIME < 0 !!!" ! C_SEATON = k_B * TEMP * X0 / (8. * Ryd * h * c * Fo) * ( (2*X_PRIME + X0) / X0 )**2 ! Write(*,*) " C_SEATON = ", C_SEATON ! ! Méthode BRENT(dichotomie, sécante et interpolation quadratique inverse) ! de résolution de l'équation F_SEATON(beta1) = 0 BETA1 = BRENT( borne_inf , borne_sup , EPSILON(1.D0), prec, F_SEATON ) ! !Write(*, *) "BETA1 = ", BETA1 If (ABS( F_SEATON( BETA1 ) ) > 1.0D4 * prec * C_SEATON ) Then !Write(*, *) "F_SEATON( BETA ) / C_SEATON , C_SEATON == ", F_SEATON( BETA1 ) / C_SEATON , C_SEATON ON_DEPASSE = ON_DEPASSE + 1 End If ! BETA_1 = REAL(BETA1) ! End Function BETA_1 !---------------------------------------------------------------------- Real Function SEATON_PHI(x) Result(y) ! USE BESSEL, BESSEL_K0 => dbesk0, BESSEL_K1 => dbesk1 ! Double Precision, intent(in) ::x ! y = x * BESSEL_K0(x) * BESSEL_K1(x) ! End Function SEATON_PHI !---------------------------------------------------------------------- Real Function SEATON_ZETA(x) Result(y) ! USE BESSEL, BESSEL_K0 => dbesk0, BESSEL_K1 => dbesk1 ! Double Precision, intent(in) :: x ! y = x**2 * (BESSEL_K0(x)**2 + BESSEL_K1(x)**2) ! End Function SEATON_ZETA !---------------------------------------------------------------------- End Module MOD_F_SEATON ! !********************************************************************** ! Module SEATON ! USE MTD_STRUCT USE MPR09_FUNCTIONS USE MOD_F_SEATON, ONLY : F_SEATON, Q_SEATON, QX_SEATON, OMEGA_SEATON, FCT_Q0, FCT_Q1, & OMEGA_0, OMEGA_1, BETA_0, BETA_1, SEATON_PHI, SEATON_ZETA, & C_SEATON, TEMP, R0, Fij, gi, gj, X0 ! Implicit None ! PRIVATE PUBLIC :: SEATON_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 :: pi = 3.14159265359D0 Double Precision, parameter :: q = 1.60217653D-19 ! Charge élémentaire A.s ! CONTAINS !---------------------------------------------------------------------- Subroutine SEATON_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, Z 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 = .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 R0 = WHICH_MEAN_RADIUS( ) ! If (BLABLA) Call BAVARD( ) ! UPSILON_TEMP_IN = UPSILON(TEMP) ! If (BLABLA) Write(*, *) "UPSILON( TEMP = ", INT(TEMP), " ) = ", UPSILON_TEMP_IN, " sans unité" ! CONTAINS !---------------------------------------------------------------------- Subroutine BAVARD( ) ! Real :: BETA0, BETA1 Logical :: TEMPERATURE = .FALSE. ! Write(*, *) "Degré d'ionisation Z = ", Z Write(*, *) TRIM( Config_i ), " => ", TRIM( Config_j ) Write(*, * ) "gi = ", gi, " => gj = ", gj Write(*, *) "DeltaE_eV = ", X0 * k_B * TEMP / q Write(*, *) "Fij = ", Fij Write(*, *) "NQP_EFFECTIF = ", NQP_EFFECTIF(En_i_cm, En_Ion_cm, Z) Write(*, *) "RADIUS_0 = ", R0 Write(*, *) "X0 = ", X0 ! If (TEMPERATURE) Then X_PRIME = 2.D0 BETA0 = BETA_0(TEMP, X_PRIME, X0, R0 ) BETA1 = BETA_1(TEMP, X_PRIME, X0, Fij) Write(*, *) "DeltaE_eV, Ej_eV = ", X0 * k_B * TEMP / q , X_PRIME * k_B * TEMP Write(*, *) "BETA_0(X_PRIME = ", X_PRIME, ") = ", BETA0 Write(*, *) "BETA_1(X_PRIME = ", X_PRIME, ") = ", BETA1 Write(*, *) "C_SEATON = ", C_SEATON Write(*, *) "F_SEATON(BETA1) = ", F_SEATON(DBLE(BETA1)) Write(*, *) "SEATON_PHI (", BETA0, ") = ", SEATON_PHI(DBLE(BETA0)) Write(*, *) "SEATON_ZETA(", BETA1, ") = ", SEATON_ZETA(DBLE(BETA1)) Write(*, *) "FCT_Q0(", X_PRIME, ") = ", FCT_Q0(TEMP, X_PRIME, X0, R0, Fij) Write(*, *) "FCT_Q1(", X_PRIME, ") = ", FCT_Q1(TEMP, X_PRIME, X0, Fij) Write(*, *) "OMEGA_0(", X_PRIME, ") = ", OMEGA_0(TEMP, X_PRIME, X0, R0, Fij) Write(*, *) "OMEGA_1(", X_PRIME, ") = ", OMEGA_1(TEMP, X_PRIME, X0, Fij) Write(*, *) "OMEGA_SEATON(", X_PRIME, ") = ", OMEGA_SEATON(X_PRIME) Write(*, *) "Q_SEATON(", X_PRIME, ") = ", Q_SEATON(X_PRIME) Write(*, *) "QX_SEATON(", X_PRIME, ") = ", QX_SEATON(X_PRIME) Write(*, *) "FCT_Q0/OMEGA_0(", X_PRIME, ") = ", FCT_Q0(TEMP, X_PRIME, X0, R0, Fij) / OMEGA_0(TEMP, X_PRIME, X0, R0, Fij) Write(*, *) "Q_SEATON/OMEGA_SEATON(", X_PRIME, ") = ", Q_SEATON(X_PRIME) / OMEGA_SEATON(X_PRIME) Write(*, *) "Ryd * h * c / (gi * k_B * TEMP * (X_PRIME + X0)) = ", Ryd * h * c / (gi * k_B * TEMP * (X_PRIME + X0)) Write(*, *) "QX_SEATON/OMEGA_SEATON(", X_PRIME, ") = ", QX_SEATON(X_PRIME) / OMEGA_SEATON(X_PRIME) Write(*, *) "Ryd * h * c / (gi * k_B * TEMP) = ", Ryd * h * c / (gi * k_B * TEMP) Write(*, *) "" End If ! End Subroutine BAVARD !---------------------------------------------------------------------- Real Function UPSILON(TEMP) ! USE INTEGRALE, ONLY : GAUSS_LAGUERRE_12PTS, GAUSS_LAGUERRE_8PTS ! Real, intent(in) :: TEMP ! Integer, parameter :: choix = 2 ! 2 par défaut (1 est douteux) Real :: C_UPSILON ! Select Case(choix) Case(1) If (BLABLA) Then Write(*, *) "CALCUL DE UPSILON AVEC OMEGA_SEATON" Write(*, *) " OM_0 OM_1" End If UPSILON = GAUSS_LAGUERRE_12PTS( OMEGA_SEATON ) Case(2) If (BLABLA) Then Write(*, *) "CALCUL DE UPSILON AVEC Q_SEATON" Write(*, *) " Q_0 Q_1" End If C_UPSILON = gi * k_B * TEMP / (Ryd * h * c) UPSILON = C_UPSILON * GAUSS_LAGUERRE_12PTS( QX_SEATON ) Case Default STOP "SEATON : FUNCTION UPSILON : FAIRE LE BON CHOIX." End Select ! End Function UPSILON !---------------------------------------------------------------------- Real Function RADIUS_0( Z, Config_i, En_i_cm, En_Ion_cm, DeltaE_cm ) Result(y) ! exprimé en unité de a0 Integer, intent(in) :: Z Double Precision, intent(in) :: En_i_cm, En_Ion_cm, DeltaE_cm Character(len = *), intent(in) :: Config_i ! Integer :: L Real :: NSTAR, mean_r Real :: K0, R0 ! Logical, parameter :: choix = .TRUE. ! si vraie utilise mean_r sinon R0 ! If (choix) Then NSTAR = NQP_EFFECTIF(En_i_cm, En_Ion_cm, Z) L = NQS(Config_i) mean_r = (1/2./Z) * (3 * NSTAR**2 - L*(L + 1)) y = mean_r Else K0 = (m_e * h * c * DeltaE_cm * 100. )**0.5 / (h / (2.D0 * pi)) ! 1/m R0 = 1.1229 / K0 / a0 y = R0 End If ! End Function RADIUS_0 !---------------------------------------------------------------------- Real Function WHICH_MEAN_RADIUS( ) Result (R0_min) ! Real :: R0_i, R0_j ! R0_min = -1. ! R0_i = RADIUS_0( Z, Config_i, En_i_cm, En_Ion_cm, DeltaE_cm ) R0_j = RADIUS_0( Z, Config_j, En_j_cm, En_Ion_cm, DeltaE_cm ) ! If (BLABLA) Then Write(*, *) "R0_i = ", R0_i, " a0" Write(*, *) "R0_j = ", R0_j, " a0" End If ! If (R0_i <= R0_j) Then R0_min = R0_i Else R0_min = R0_j End If ! End Function WHICH_MEAN_RADIUS !---------------------------------------------------------------------- End Subroutine SEATON_BB ! End Module SEATON