Line data    Source code

       1              : ! ***********************************************************************
       2              : !
       3              : !   Copyright (C) 2022  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 ion_offset
      21              :       use const_def, only: dp, amu, ev2erg
      22              :       use math_lib
      23              : 
      24              :       implicit none
      25              : 
      26              :       logical, parameter :: dbg = .false.
      27              : 
      28              :       private
      29              :       public :: compute_ion_offset
      30              : 
      31              :       contains
      32              : 
      33              :       !> Computes the energy required to fully ionize
      34              :       !! the specified composition. This is an offset
      35              :       !! between Skye/HELM and FreeEOS.
      36              :       !!
      37              :       !! @param species Number of species
      38              :       !! @param xa Mass fractions (vector)
      39              :       !! @param offset Energy in erg/g
      40        45520 :       real(dp) function compute_ion_offset(species, xa, chem_id) result(offset)
      41              :          use chem_def, only: chem_isos
      42              :          integer, intent(in) :: species
      43              :          real(dp), intent(in) :: xa(species)
      44              :          integer, pointer :: chem_id(:)
      45              : 
      46              :          ! First 28 ionization energies, same species as used in FreeEOS, obtained from
      47              :          ! NIST (https://physics.nist.gov/PhysRefData/ASD/ionEnergy.html)
      48              :          ! using the query 'H-Ni' asking for 'Total binding energy'.
      49              :          real(dp), parameter :: ionization_table(28) = [13.598434599702,79.0051545,203.4861711,399.14864,670.9809,&
      50              :                                                    1030.1085,1486.058,2043.8428,2715.89,3511.696,4420.0,5451.06,6604.95,&
      51              :                                                    7888.53,9305.8,10859.7,12556.4,14400.8,16382.0,18510.0,20788.0,&
      52              :                                                    23221.0,25820.0,28582.0,31514.0,34619.0,37899.0,41356.0]
      53              : 
      54        91040 :          integer :: k, Z(species)
      55      1518460 :          real(dp) :: A(species), ya(species), norm
      56              : 
      57              :          ! Get basic species info
      58              :          norm = 0d0
      59       759230 :          do k=1,species
      60       713710 :             A(k) = chem_isos% Z_plus_N(chem_id(k))  ! baryon number
      61       713710 :             Z(k) = chem_isos% Z(chem_id(k))  ! charge
      62       713710 :             ya(k) = xa(k) / A(k)            ! number fraction (not normalized)
      63       759230 :             norm = norm + ya(k)             ! accumulate number fraction normalization
      64              :          end do
      65              : 
      66              :          ! Normalize number fraction
      67       759230 :          do k=1,species
      68       759230 :             ya(k) = ya(k) / norm
      69              :          end do
      70              : 
      71              :          ! Compute offset in eV/baryon
      72              :          offset = 0d0
      73       759230 :          do k=1,species
      74       759230 :             if (Z(k) <= 28 .and. Z(k) >= 1) then
      75       686820 :                offset = offset + ionization_table(Z(k)) * ya(k)
      76              :             end if
      77              :          end do
      78              : 
      79              :          ! Translate that into erg/g
      80        45520 :          offset = offset * ev2erg / amu
      81              : 
      82        45520 :       end function compute_ion_offset
      83              : 
      84              : end module ion_offset