MODULE AnalyticCooling2
   USE GlobalDeclarations
   USE PhysicsDeclarations
   REAL(KIND=qPrec) :: A_alpha,A_beta
   
CONTAINS

   !> Finalize initialization of analytic cooling source
   SUBROUTINE InitAnalyticCool2()
      !For power laws that have the same strength as brehm strahlung at a given shock speed (cs*Mach) 
      !use alpha=4.76e-20 (ergs*cm^3/s/K^.5) * (3d0/16d0*TempScale*(cs*Mach)**2)**(.5d0)  * (3d0/16d0 * (cs*Mach)**2) **(-coolingobj%beta)

      !Note if beta=.5 this is just pure Brehmstrahlung and
      ! alpha = 4.76e-20 (ergs*cm^3/s/K^.5) * (TempScale)**(.5d0)
      !See Imamura, Wolff, and Durisen (1983) although their alpha is our beta and their lambda is our alpha

      !alpha should have units of ergs*cm^3/s/K^beta
   END SUBROUTINE InitAnalyticCool2



  !>  Analytic cooling source term
  !! @details Analytic cooling source term, where the cooling
  !!          function is a single power law, alpha*T^beta
  !! @param q Fluid variable vector q
  !! @param dqdt Update to be applied to q
  !! @param coolingobj Cooling object
  SUBROUTINE Analytic_Cooling2(q,dqdt, Temp)
    ! Interface declarations
    REAL(KIND=qPrec) :: q(:), Temp
    REAL(KIND=qPrec),INTENT(INOUT) :: dqdt(:)
    ! Local declarations

    dqdt(iE) = dqdt(iE)-AnalyticCoolingStrength2(q, Temp)
 END SUBROUTINE Analytic_Cooling2
  
  FUNCTION AnalyticCoolingStrength2(q, Temp)
     REAL(KIND=qPREC) :: AnalyticCoolingStrength2
     REAL(KIND=qPrec) :: q(:)
     ! Local declarations
     REAL(KIND=qPrec) :: Temp, Temp0, T0

     Temp0 = 360000  !10^5 Kelvin
     T0 = Temp/Temp0
     AnalyticCoolingStrength2=(q(1)-q(9))**2 * A_alpha*T0**A_beta*ScaleCool

     !! only for ambient tracer > 1/3 which is the scaless value for initial ambient density
     !! Jet density is the initial jet density
     if( (q(1)-q(9)) <= 2d0 ) AnalyticCoolingStrength2=0d0
        
     
  END FUNCTION AnalyticCoolingStrength2

END MODULE AnalyticCooling2
