Line data    Source code

       1              : ! ***********************************************************************
       2              : !
       3              : !   Copyright (C) 2010-2021  The MESA Team
       4              : !
       5              : !   This program is free software: you can redistribute it and/or modify
       6              : !   it under the terms of the GNU Lesser General Public License
       7              : !   as published by the Free Software Foundation,
       8              : !   either version 3 of the License, or (at your option) any later version.
       9              : !
      10              : !   This program is distributed in the hope that it will be useful,
      11              : !   but WITHOUT ANY WARRANTY; without even the implied warranty of
      12              : !   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
      13              : !   See the GNU Lesser General Public License for more details.
      14              : !
      15              : !   You should have received a copy of the GNU Lesser General Public License
      16              : !   along with this program. If not, see <https://www.gnu.org/licenses/>.
      17              : !
      18              : ! ***********************************************************************
      19              : 
      20              :       module mod_typical_charge
      21              : 
      22              :       use const_def, only: one_third, two_thirds, dp
      23              :       use math_lib
      24              : 
      25              :       implicit none
      26              : 
      27              :       private
      28              :       public :: eval_typical_charge
      29              : 
      30              :       logical, parameter :: dbg = .false.
      31              : 
      32              :       contains
      33              : 
      34              :       ! average ionic charges based on
      35              :       ! C. Paquette, C. Pelletier, G. Fontaine, G. Michaud,
      36              :       ! "Diffusion in White Dwarfs: New Results and Comparative Study",
      37              :       ! ApJ Supp. Series, 61:197-217, 1986.
      38              : 
      39              :       ! Originally used eqn 21 of above paper for depression of the continuum,
      40              :       ! but that expression was missing a rho^1/3, so now corrected to match
      41              :       ! eqn 3 of
      42              :       ! J. Dupuis, G. Fontaine, C. Pelletier, F. Wesemael,
      43              :       ! "A Study of Metal Abundance Patterns in Cool White Dwarfs I.
      44              :       ! Time-dependent Calculations of Gravitational Settling",
      45              :       ! Apj Supp. Series, 82:505-521, 1992
      46              : 
      47              :       ! uses ionization potentials from
      48              :       ! Allen, C.W., 1973, "Astrophysical Quantities", 3rd edition, pg 37-38.
      49              : 
      50              :       ! tables go from H to Zn.
      51              : 
      52              :       ! the routine is written so that it doesn't ever return 0 as a typical charge.
      53              :       ! these are "typical" charges rather than average.  the values are whole numbers.
      54              : 
      55            0 :       real(dp) function eval_typical_charge( &
      56              :             cid, abar, free_e, T, log10_T, rho, log10_rho)
      57              :          integer, intent(in) :: cid  ! chem id such as ihe4.  defined in chem_def.
      58              :          real(dp), intent(in) :: abar  ! average atomic weight (from chem_lib)
      59              :          real(dp), intent(in) :: free_e
      60              :             ! mean number of free electrons per nucleon (from eos_lib)
      61              :             ! abar*free_e = (nucleons/particle)*(charge/nucleon) = charge/particle = z1
      62              :          real(dp), intent(in) :: T, log10_T, rho, log10_rho
      63              :          eval_typical_charge = get_typical_charge( &
      64            0 :             cid, abar, abar*free_e, T, log10_T, rho, log10_rho)
      65            0 :       end function eval_typical_charge
      66              : 
      67              : 
      68            0 :       subroutine chi_info(a1, z1, T, log_T, rho, log_rho, chi, c0, c1, c2)
      69              :          real(dp), intent(in) :: a1, z1, T, log_T, rho, log_rho
      70              :          real(dp), intent(out) :: chi, c0, c1, c2
      71            0 :          chi = 1.987d-4*T*(-8.392d0 - log_rho + 1.5d0*log_T - log10(z1/a1))  ! eqn 20
      72              :          ! coef's used in eqn 21
      73            0 :          c0 = 1.085d-4*rho*T/a1
      74            0 :          c1 = 1.617d4*sqrt(rho*(z1*z1 + z1)/(T*a1))
      75            0 :          c2 = 29.38d0*z1*pow(rho/a1,one_third)
      76              :          ! c2 had a typo in eqn 21, now corrected to match Dupuis et al. (1992) eqn 3
      77            0 :       end subroutine chi_info
      78              : 
      79            0 :       real(dp) function chi_effective(chi, c0, c1, c2, z1, z2)
      80              :          real(dp), intent(in) :: chi, c0, c1, c2, z1, z2
      81              :          chi_effective = chi + c0/(z2*z2*z2) + &
      82            0 :             min(c1*z2, c2*(pow(z2/z1,two_thirds) + 0.6d0))
      83            0 :       end function chi_effective
      84              : 
      85            0 :       real(dp) function get_typical_charge(cid, a1, z1, T, log_T, rho, log_rho)
      86              :          use ionization_potentials
      87              :          use chem_def
      88              :          integer, intent(in) :: cid
      89              :          real(dp), intent(in) :: a1, z1, T, log_T, rho, log_rho
      90            0 :          real(dp) :: chi, c0, c1, c2, z2, chi_eff
      91              :          integer :: i, izmax
      92              :          include 'formats'
      93            0 :          if (.not. ionization_tables_okay) then
      94            0 :             call set_ionization_potentials
      95            0 :             ionization_tables_okay = .true.
      96              :          end if
      97            0 :          izmax = int(chem_isos% Z(cid))
      98            0 :          get_typical_charge = dble(izmax)
      99            0 :          if (izmax > 30) return
     100            0 :          call chi_info(a1, z1, T, log_T, rho, log_rho, chi, c0, c1, c2)
     101            0 :          do i=1, izmax-1
     102            0 :             z2 = dble(i)
     103            0 :             chi_eff = chi_effective(chi, c0, c1, c2, z1, z2+1)
     104            0 :             if (chi_eff < ip(izmax,i+1)) then
     105            0 :                get_typical_charge = z2
     106            0 :                return
     107              :             end if
     108              :          end do
     109            0 :       end function get_typical_charge
     110              : 
     111              :       end module mod_typical_charge