! ***********************************************************************
!
!   Copyright (C) 2010-2019  Bill Paxton & The MESA Team
!
!   this file is part of mesa.
!
!   mesa is free software; you can redistribute it and/or modify
!   it under the terms of the gnu general library public license as published
!   by the free software foundation; either version 2 of the license, or
!   (at your option) any later version.
!
!   mesa is distributed in the hope that it will be useful, 
!   but without any warranty; without even the implied warranty of
!   merchantability or fitness for a particular purpose.  see the
!   gnu library general public license for more details.
!
!   you should have received a copy of the gnu library general public license
!   along with this software; if not, write to the free software
!   foundation, inc., 59 temple place, suite 330, boston, ma 02111-1307 usa
!
! ***********************************************************************
 
      module run_star_extras

      use star_lib
      use star_def
      use const_def
      use math_lib
      
      implicit none
      
      ! these routines are called by the standard run_star check_model
      contains
      
     
      subroutine extras_controls(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         ! this is the place to set any procedure pointers you want to change
         ! e.g., other_wind, other_mixing, other_energy  (see star_data.inc)


         ! the extras functions in this file will not be called
         ! unless you set their function pointers as done below.
         ! otherwise we use a null_ version which does nothing (except warn).

         s% extras_startup => extras_startup
         s% extras_start_step => extras_start_step
         s% extras_check_model => extras_check_model
         s% extras_finish_step => extras_finish_step
         s% extras_after_evolve => extras_after_evolve
         s% how_many_extra_history_columns => how_many_extra_history_columns
         s% data_for_extra_history_columns => data_for_extra_history_columns
         s% how_many_extra_profile_columns => how_many_extra_profile_columns
         s% data_for_extra_profile_columns => data_for_extra_profile_columns  

         s% how_many_extra_history_header_items => how_many_extra_history_header_items
         s% data_for_extra_history_header_items => data_for_extra_history_header_items
         s% how_many_extra_profile_header_items => how_many_extra_profile_header_items
         s% data_for_extra_profile_header_items => data_for_extra_profile_header_items

         s% other_wind => eval_Leuven_wind

      end subroutine extras_controls
      
      
      subroutine extras_startup(id, restart, ierr)
         integer, intent(in) :: id
         logical, intent(in) :: restart
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      end subroutine extras_startup
      

      integer function extras_start_step(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_start_step = 0
      end function extras_start_step


      ! returns either keep_going, retry, backup, or terminate.
      integer function extras_check_model(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_check_model = keep_going         
         if (s% center_he4 < 1.d-3 .and. s% center_c12 < 1.d-3) then
            ! stop when star hydrogen mass drops to specified level
            extras_check_model = terminate
            write(*, *) 'He4 and C12 depleted simultaneously'
            return
         end if


         ! if you want to check multiple conditions, it can be useful
         ! to set a different termination code depending on which
         ! condition was triggered.  MESA provides 9 customizeable
         ! termination codes, named t_xtra1 .. t_xtra9.  You can
         ! customize the messages that will be printed upon exit by
         ! setting the corresponding termination_code_str value.
         ! termination_code_str(t_xtra1) = 'my termination condition'

         ! by default, indicate where (in the code) MESA terminated
         if (extras_check_model == terminate) s% termination_code = t_extras_check_model
      end function extras_check_model


      integer function how_many_extra_history_columns(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_columns = 0
      end function how_many_extra_history_columns
      
      
      subroutine data_for_extra_history_columns(id, n, names, vals, ierr)
         integer, intent(in) :: id, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         ! note: do NOT add the extras names to history_columns.list
         ! the history_columns.list is only for the built-in history column options.
         ! it must not include the new column names you are adding here.
         

      end subroutine data_for_extra_history_columns

      
      integer function how_many_extra_profile_columns(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_profile_columns = 0
      end function how_many_extra_profile_columns
      
      
      subroutine data_for_extra_profile_columns(id, n, nz, names, vals, ierr)
         integer, intent(in) :: id, n, nz
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(nz,n)
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         integer :: k
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         
         ! note: do NOT add the extra names to profile_columns.list
         ! the profile_columns.list is only for the built-in profile column options.
         ! it must not include the new column names you are adding here.

         ! here is an example for adding a profile column
         !if (n /= 1) stop 'data_for_extra_profile_columns'
         !names(1) = 'beta'
         !do k = 1, nz
         !   vals(k,1) = s% Pgas(k)/s% P(k)
         !end do
         
      end subroutine data_for_extra_profile_columns


      integer function how_many_extra_history_header_items(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_history_header_items = 0
      end function how_many_extra_history_header_items


      subroutine data_for_extra_history_header_items(id, n, names, vals, ierr)
         integer, intent(in) :: id, n
         character (len=maxlen_history_column_name) :: names(n)
         real(dp) :: vals(n)
         type(star_info), pointer :: s
         integer, intent(out) :: ierr
         ierr = 0
         call star_ptr(id,s,ierr)
         if(ierr/=0) return

         ! here is an example for adding an extra history header item
         ! also set how_many_extra_history_header_items
         ! names(1) = 'mixing_length_alpha'
         ! vals(1) = s% mixing_length_alpha

      end subroutine data_for_extra_history_header_items


      integer function how_many_extra_profile_header_items(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         how_many_extra_profile_header_items = 0
      end function how_many_extra_profile_header_items


      subroutine data_for_extra_profile_header_items(id, n, names, vals, ierr)
         integer, intent(in) :: id, n
         character (len=maxlen_profile_column_name) :: names(n)
         real(dp) :: vals(n)
         type(star_info), pointer :: s
         integer, intent(out) :: ierr
         ierr = 0
         call star_ptr(id,s,ierr)
         if(ierr/=0) return

         ! here is an example for adding an extra profile header item
         ! also set how_many_extra_profile_header_items
         ! names(1) = 'mixing_length_alpha'
         ! vals(1) = s% mixing_length_alpha

      end subroutine data_for_extra_profile_header_items


      ! returns either keep_going or terminate.
      ! note: cannot request retry or backup; extras_check_model can do that.
      integer function extras_finish_step(id)
         integer, intent(in) :: id
         integer :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
         extras_finish_step = keep_going

         ! to save a profile, 
            ! s% need_to_save_profiles_now = .true.
         ! to update the star log,
            ! s% need_to_update_history_now = .true.

         ! see extras_check_model for information about custom termination codes
         ! by default, indicate where (in the code) MESA terminated
         if (extras_finish_step == terminate) s% termination_code = t_extras_finish_step
      end function extras_finish_step
      
      
      subroutine extras_after_evolve(id, ierr)
         integer, intent(in) :: id
         integer, intent(out) :: ierr
         type (star_info), pointer :: s
         ierr = 0
         call star_ptr(id, s, ierr)
         if (ierr /= 0) return
      end subroutine extras_after_evolve


      subroutine eval_Leuven_wind(id, Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z, w, ierr)

         use star_def

         integer, intent(in) :: id
         real(dp), intent(in) :: Lsurf, Msurf, Rsurf, Tsurf, X, Y, Z
         real(dp) :: kappa, vturb
         real(dp) :: cs, Gamma, Rpmod, HoR, rhoo, rho, mdotAn
         real(dp) :: alfa, Teff_jump, dT, vinf_div_vesc
         real(dp) :: wLh, wVh, w1, w2, wNLh, wh, wdJ, wLc, wc
         real(dp) :: switch1, switch2
         real(dp) :: log10wLh, log10wVh, log10wdJ, log10wLc
	 real(dp) :: Zsolar = 0.019d0
         real(dp), intent(out) :: w
         integer, intent(out) :: ierr

         include 'formats'

	 ! Leuven Hot Wind (Bjorklund et al. 2020)

         log10wLh = -5.55 + 0.79*log10(Z/Zsolar) &
                    +(2.16 - 0.32*log10(Z/Zsolar))*log10(Lsurf/Lsun/1.d6)
	 wLh = exp10(log10wLh)


         ! 1/3 Vink Hot Wind as Teff < 28 kK extrapolation of Bjorklund (Vink et al. 2001, Bjorklund priv. comm.)

            if (Tsurf < 22500d0) then
               alfa = 0
            else ! use Vink et al 2001, eqns 14 and 15 to set "jump" temperature
               Teff_jump = 1d3*(61.2d0 + 2.59d0*(-13.636d0 + 0.889d0*log10(Z/Zsolar)))
               dT = 100d0
               if (Tsurf > Teff_jump + dT) then
                  alfa = 1
               else if (Tsurf < Teff_jump - dT) then
                  alfa = 0
               else
                  alfa = (Tsurf - (Teff_jump - dT)) / (2*dT)
               end if
            end if

            if (alfa > 0) then ! eval hot side wind (eqn 24)
               vinf_div_vesc = 2.6d0 ! this is the hot side galactic value
               vinf_div_vesc = vinf_div_vesc*pow(Z/Zsolar,0.13d0) ! corrected for Z
               log10wVh = &
                  - 6.697d0 &
                  + 2.194d0*log10(Lsurf/Lsun/1d5) &
                  - 1.313d0*log10(Msurf/Msun/30) &
                  - 1.226d0*log10(vinf_div_vesc/2d0) &
                  + 0.933d0*log10(Tsurf/4d4) &
                  - 10.92d0*pow2(log10(Tsurf/4d4)) &
                  + 0.85d0*log10(Z/Zsolar)
               w1 = exp10(log10wVh)
            else
               w1 = 0
            end if

            if (alfa < 1) then ! eval cool side wind (eqn 25)
               vinf_div_vesc = 1.3d0 ! this is the cool side galactic value
               vinf_div_vesc = vinf_div_vesc*pow(Z/Zsolar,0.13d0) ! corrected for Z
               log10wVh = &
                  - 6.688d0 &
                  + 2.210d0*log10(Lsurf/Lsun/1d5) &
                  - 1.339d0*log10(Msurf/Msun/30) &
                  - 1.601d0*log10(vinf_div_vesc/2d0) &
                  + 1.07d0*log10(Tsurf/2d4) &
                  + 0.85d0*log10(Z/Zsolar)
               w2 = exp10(log10wVh)
            else
               w2 = 0
            end if

            wVh = (alfa*w1 + (1 - alfa)*w2)/3.d0

         !  helium rich Wolf-Rayet star wind (Nugis & Lamers 2000)

         wNLh = 1d-11 * pow(Lsurf/Lsun,1.29d0) * pow(Y,1.7d0) * sqrt(Z) ! helium rich Wolf-Rayet star: Nugis & Lamers

         
         ! Hot Wind Selection Criteria

         if (X < 0.4d0) then 
            wh = wNLh
            write(*,*) 'Nugis & Lamers'
         else if (Tsurf < 2.8e4) then
            wh = wVh
         else
            wh = wLh
         endif


         ! Leuven Cool Wind (Kee et al. 2020)

         kappa = 0.01e0
         vturb = 18.2e5 

         cs = sqrt(1.38e-16 * Tsurf/1.6e-24)
         Gamma = kappa * Lsurf/(4.0 * 3.14159 * 6.67e-8 * Msurf * 3e10)
         Rpmod = 6.67e-8 * Msurf * (1.0 - Gamma)/(2.0 * (cs*cs + vturb*vturb))
         HoR = 0.5 * Rsurf/Rpmod
         rhoo = 2.e0/3.e0 * 1.e0/(kappa * HoR * Rsurf)/(1 - exp(-1/HoR))
         rho = rhoo * exp(-1.e0/HoR*(1-Rsurf/Rpmod)) * 0.606531
         mdotAn =  4.e0 * 3.14159 * rho * sqrt(vturb*vturb + cs*cs) * Rpmod*Rpmod/Msun*3600*24*365.25
    
         wLc = mdotAn * ((vturb/17e5)/(sqrt(2.e0 * 6.67e-8 * Msurf/Rsurf)/60e5))**1.30


         ! de Jager cool wind (de Jager et al. 1988)

         log10wdJ = -8.158d0 &
                  +1.769d0*log10(Lsurf/Lsun) &
                  -1.676d0*log10(Tsurf)
         wdJ = exp10(log10wdJ)


        ! Cool Wind Selection Criteria

        if (Tsurf > 5000.d0) then
          wc = max(wdJ,wLc)
          if (wdJ > wLc) then
             write(*,*) "de Jager:", wdJ
          else
             write(*,*) "Leuven:", wLc
          endif
	  ierr = 0
        else
	  wc = wLc
	  ierr = 0
        endif


         ! Interpolation between Hot and Cool

         switch1 = 10000.d0
         switch2 = 8000.d0

         if (Tsurf > switch1) then
             w = wh
             ierr = 0
         else if (Tsurf > switch2) then
             w = (wh*(Tsurf - switch2) + wc*(switch1-Tsurf))/(switch1-switch2)
             ierr = 0
         else
             w = wc
             ierr = 0
         end if

      end subroutine eval_Leuven_wind

      end module run_star_extras
      
