! ***********************************************************************
!
!   Copyright (C) 2012  Bill Paxton
!
!   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_binary_extras 

      use star_lib
      use star_def
      use const_def
      use chem_def
      use num_lib
      use binary_def
      use crlibm_lib
      
      
      implicit none
      
      contains
      
      subroutine extras_binary_controls(binary_id, ierr)
         integer :: binary_id
         integer, intent(out) :: ierr
         type (binary_info), pointer :: b
         ierr = 0

         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if

         ! Set these function pinters to point to the functions you wish to use in
         ! your run_binary_extras. Any which are not set, default to a null_ version
         ! which does nothing.
         b% how_many_extra_binary_history_columns => how_many_extra_binary_history_columns
         b% data_for_extra_binary_history_columns => data_for_extra_binary_history_columns

         b% extras_binary_startup=> extras_binary_startup
         b% extras_binary_check_model=> extras_binary_check_model
         b% extras_binary_finish_step => extras_binary_finish_step
         b% extras_binary_after_evolve=> extras_binary_after_evolve

         ! Once you have set the function pointers you want, then uncomment this (or set it in your star_job inlist)
         ! to disable the printed warning message,
          b% warn_binary_extra =.false.
          b% other_jdot_ml => chen_jdot_ml
          b% other_jdot_gr => chen_jdot_gr
      end subroutine extras_binary_controls

      integer function how_many_extra_binary_history_columns(binary_id)
         use binary_def, only: binary_info
         integer, intent(in) :: binary_id
         how_many_extra_binary_history_columns = 0
      end function how_many_extra_binary_history_columns

      !Added model program, By Chen   
      subroutine chen_jdot_gr(binary_id, ierr)
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr
         type (binary_info), pointer :: b
         real(dp) :: bsep4, clight5, cgrav3
         ierr = 0
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if
         bsep4 = pow4(b% separation)
         clight5 = pow5(clight)
         cgrav3 = b% s_donor% cgrav(1)*b% s_donor% cgrav(1)*b% s_donor% cgrav(1)
         b% jdot_gr = -32d0 * cgrav3 * b% m(b% a_i) * b% m(b% d_i) * (b% m(b% a_i) + b% m(b% d_i)) / &
             (5d0 * clight5 * bsep4) * b% angular_momentum_j

      end subroutine chen_jdot_gr
    
      subroutine chen_jdot_ml(binary_id, ierr)
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr
         integer :: model_number_last
         integer, parameter :: flag_ml = 3
         !flag_ml = 1, 2, 3 for mass lost from vicinity of donor, accretor, or circumbinary coplanar toroid

         common /binary/ model_number_last
         common/wind/ irradiation_wind
         type (binary_info), pointer :: b
         real(dp) :: irradiation_wind, L_edd, Lx, clight2
         real, parameter :: eff_wind = 0.0
         !eff_wind, from 0 to 1, wind-driving effiency
         real, parameter :: dotM_cr_gs = 6.3077d18
         !critical accretion rate in unit of g/s
         common /luninosity/ Lx, L_edd
         real :: eta_acc = 0.1
         !accretion-luninosity effiency, from 0 to 1, 0.0001 for WD, 0.1 for NS, 0.1-0.4 for BH ... 
         
         ierr = 0
         call binary_ptr(binary_id, b, ierr)

         clight2 = pow2(clight)

         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if

         if (flag_ml == 1) then
         !mass lost from vicinity of donor
            b% jdot_ml = (b% mdot_system_transfer(b% d_i) + b% mdot_system_wind(b% d_i))*&
               (b% m(b% a_i)/(b% m(b% a_i)+b% m(b% d_i))*b% separation)**2*2*pi/b% period *&
               sqrt(1 - b% eccentricity**2)
         else if (flag_ml == 2) then
         !mass lost from vicinity of accretor
            b% jdot_ml = b% jdot_ml + (b% mdot_system_transfer(b% a_i) + b% mdot_system_wind(b% a_i))*&
               (b% m(b% d_i)/(b% m(b% a_i)+b% m(b% d_i))*b% separation)**2*2*pi/b% period *&
               sqrt(1 - b% eccentricity**2)
         else if (flag_ml == 3) then
         !mass lost from circumbinary coplanar toroid
            b% jdot_ml = b% jdot_ml + b% mdot_system_cct * b% mass_transfer_gamma * &
               sqrt(b% s_donor% cgrav(1) * (b% m(1) + b% m(2)) * b% separation)
         end if 

         L_edd = (b% m(b% a_i) * 1.3d38)/msol
        
         if (abs(b% mtransfer_rate)*b% xfer_fraction.gt.dotM_cr_gs .and. &
            b% mdot_edd.gt.abs(b% mtransfer_rate)*b% xfer_fraction) then
            Lx = eta_acc * (abs(b% mtransfer_rate)*b% xfer_fraction)*clight2
         end if


         if (dotM_cr_gs.gt.(abs(b% mtransfer_rate)*b% xfer_fraction)) then
            Lx = eta_acc * (abs(b% mtransfer_rate)*b% xfer_fraction)**2*clight2/dotM_cr_gs 
         end if
        
         irradiation_wind = eff_wind*Lx*min(b% r(b% d_i), b% rl(b% d_i))**3/ &
            (4*standard_cgrav*b% m(b% d_i)*b% separation**2.0)
         
         if (b% s_donor% model_number.gt.model_number_last) then
            write(15,*) b% s_donor% model_number, b% m(b% a_i), b% m(b% d_i), b% separation, & 
               b% period, b% binary_age, Lx, b% mtransfer_rate, irradiation_wind, b% xfer_fraction  
         end if
    
         model_number_last = b% s_donor% model_number

      end subroutine chen_jdot_ml
      !end added model program

      subroutine data_for_extra_binary_history_columns(binary_id, n, names, vals, ierr)
         use const_def, only: dp
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id
         integer, intent(in) :: n
         character (len=maxlen_binary_history_column_name) :: names(n)
         real(dp) :: vals(n)
         integer, intent(out) :: ierr
         real(dp) :: beta
         ierr = 0
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then
            write(*,*) 'failed in binary_ptr'
            return
         end if
         
      end subroutine data_for_extra_binary_history_columns
      
      
      integer function extras_binary_startup(binary_id,restart,ierr)
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr
         logical, intent(in) :: restart
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if
         
!          b% s1% job% warn_run_star_extras = .false.
          extras_binary_startup = keep_going
      end function  extras_binary_startup
      
      !Return either rety,backup,keep_going or terminate
      integer function extras_binary_check_model(binary_id)
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id
         integer :: ierr
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if  
         extras_binary_check_model = keep_going
        
      end function extras_binary_check_model
      
      
      ! returns either keep_going or terminate.
      ! note: cannot request retry or backup; extras_check_model can do that.
      integer function extras_binary_finish_step(binary_id)
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id
         integer :: ierr
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if  
         extras_binary_finish_step = keep_going
         
      end function extras_binary_finish_step
      
      subroutine extras_binary_after_evolve(binary_id, ierr)
         type (binary_info), pointer :: b
         integer, intent(in) :: binary_id
         integer, intent(out) :: ierr
         call binary_ptr(binary_id, b, ierr)
         if (ierr /= 0) then ! failure in  binary_ptr
            return
         end if      
         
 
      end subroutine extras_binary_after_evolve     
      
      end module run_binary_extras
