! -------------------------------------------------------------------------
! Evaluate effective core mass on the main sequence, as a function of
! stellar (ZAMS) mass and metallity. Useful for calculating the lifetime
! of stellar collision products.
!
! See Glebbeek & Pols, 2008 A&A 488 for details.
!
! Terms of use:
!
! THIS CODE IS PROVIDED ON AN 'AS-IS' BASIS, WITHOUT ANY CLAIMS ABOUT ITS
! CORRECTNESS OR USEFULNESS. IN PARTICULAR WE WILL NOT ACCEPT
! RESPONSIBILITY IF THIS CODE CAUSES YOUR COMPUTER TO CRASH, BLOW UP, SING
! A SONG OR MISBEHAVE IN ANY OTHER WAY. USE AT YOUR OWN RISK!
!
! That said, this code has worked fine and without problems for me.
! If you use this code and find it at all useful I would be very interested
! in hearing from you!
! If you use this code for published work let me know. I may be interested
! in collaborating on your project! At the very least, you are REQUIRED to
! cite the above paper.
! If you cannot accept these terms, you should delete this code from your
! computer.
!
! Evert Glebbeek           March 13 2009
! -------------------------------------------------------------------------



! -------------------------------------------------------------------------
!     get_hydrogen_burning_fraction:
!
!     Returns the mass ratio qc of the total amount of hydrogen burned on
!     the main sequence to the total stellar mass.
!
!     Input parameters:
!      * mass  - ZAMS mass of the star
!      * Z     - metallicity
! -------------------------------------------------------------------------
      function get_hydrogen_burning_fraction(mass, Z)
      implicit none
      double precision, intent(in) :: mass, Z
      double precision :: get_hydrogen_burning_fraction
!     Fitting coefficients, after Glebbeek & Pols, for the mass range
!     0.4-75Msun. For masses outside this range we just use the boundary
!     values. Metallicities below 0.001 and above 0.02 are likewise
!     truncated at the available metallicties.
      double precision, parameter :: qc_coeff_02(8) =
     &         (/ -0.685213d0, 0.289269d0, 0.0123223d0, 3.33570d-6,
     &             2.709640d0, 1.449630d0, 0.6291210d0, 4.93117d-4 /)
      double precision, parameter :: qc_coeff_001(8) =
     &         (/ -0.682671d0, 0.185036d0, 0.00658021d0, 7.5450d-7,
     &             2.714370d0, 0.334367d0, 0.34388900d0, 1.8516d-4 /)

      double precision :: x, ox
      double precision :: m, m2, m3, m4, m5, m6, m7
      double precision :: numerator, denominator
      double precision :: c(8)

!     Determine mass for calculation
      m = max(0.4, min(mass, 75.0))
!     Determine fitting coefficients
      if (Z >= 0.02d0) then
         c(:) = qc_coeff_02(:)
      elseif (Z <= 0.001d0) then
         c(:) = qc_coeff_001(:)
      else
         x = (0.02d0 - Z) / (0.02d0 - 0.001d0)
         ox = 1.0d0 - x
         c(:) = x * qc_coeff_001(:) + ox * qc_coeff_02(:)
      end if

!     Powers of m
      m2 = m*m
      m3 = m2*m
      m4 = m2*m2
      m5 = m3*m2
      m6 = m3*m3
      m7 = m5*m2
!     Rational function fit to qc(M,Z)
      numerator   =   1  + c(1)*m  + c(2)*m3 + c(3)*m5 + c(4)*m7
      denominator = c(5) + c(6)*m2 + c(7)*m4 + c(8)*m6
      get_hydrogen_burning_fraction = numerator / denominator
      end function



! -------------------------------------------------------------------------
!     get_mixing_parameter:
!
!     Returns the mixing parameter ("alpha") for collisions between main
!     sequence stars.
!     See Glebbeek & Pols, 2008 A&A 488 for details.
!
!     Input parameters:
!      * Z     - metallicity
! -------------------------------------------------------------------------
      function get_mixing_parameter(Z)
      implicit none
      double precision, intent(in) :: Z
      double precision :: get_mixing_parameter
      double precision :: x, ox
      double precision, parameter :: a02 = 1.67, a001 = 1.43

      if (Z >= 0.02d0) then
         get_mixing_parameter = a02
      elseif (Z <= 0.001d0) then
         get_mixing_parameter = a001
      else
         x = (0.02d0 - Z) / (0.02d0 - 0.001d0)
         ox = 1.0d0 - x
         get_mixing_parameter = x * a001 + ox * a02
      end if

      end function