
import numpy as np
from astropy import units as u
from scipy import special


def encmass(r, norm, rs, alpha, beta, gamma):
    
    """
    Enclosed mass profile of a generalised Hernquist model.
    
    INPUTS
      r     : radial variable (requires unit)
      norm  : normalisation (requires unit)
      rs    : scale radius of model (requires unit)
      alpha : sharpness of transition between inner and outer
      beta  : outer logarithmic slope
      gamma : inner logarithmic slope
    """
    
    a = (3.-gamma)/alpha
    b = (gamma-beta)/alpha
    y = (r/rs)**alpha
    
    fn = lambda x: x**a * special.hyp2f1(a, -b, 1+a, -x)/a
    
    encmass = (4*np.pi*norm*rs**3*fn(y.value)).to(u.Msun)
    
    return encmass


import numpy as np
from astropy import units as u
from scipy import special


def encmass(r, norm, rs, alpha, beta, gamma):
    
    """
    Enclosed mass profile of a generalised Hernquist model.
    
    INPUTS
      r     : radial variable (requires unit)
      norm  : normalisation (requires unit)
      rs    : scale radius of model (requires unit)
      alpha : sharpness of transition between inner and outer
      beta  : outer logarithmic slope
      gamma : inner logarithmic slope
    """
    
    a = (3.-gamma)/alpha
    b = (gamma-beta)/alpha
    y = ((r/rs).to(u.dimensionless_unscaled).value)**alpha
    
    fn = lambda x: x**a * special.hyp2f1(a, -b, 1+a, -x)/a
    
    encmass = (4*np.pi*norm*rs**3*fn(y)/alpha).to(u.Msun)
    
    return encmass

Fix error in enclosed mass calculation and fix unit handling.
Fix error in enclosed mass calculation and fix unit handling.

  
  import numpy as np
  from astropy import units as u
  from scipy import special
  
  
  def encmass(r, norm, rs, alpha, beta, gamma):
      
      """
      Enclosed mass profile of a generalised Hernquist model.
      
      INPUTS
        r     : radial variable (requires unit)
        norm  : normalisation (requires unit)
        rs    : scale radius of model (requires unit)
        alpha : sharpness of transition between inner and outer
        beta  : outer logarithmic slope
        gamma : inner logarithmic slope
      """
      
      a = (3.-gamma)/alpha
      b = (gamma-beta)/alpha
-     y = (r/rs)**alpha
+     y = ((r/rs).to(u.dimensionless_unscaled).value)**alpha
      
      fn = lambda x: x**a * special.hyp2f1(a, -b, 1+a, -x)/a
      
-     encmass = (4*np.pi*norm*rs**3*fn(y.value)).to(u.Msun)
+     encmass = (4*np.pi*norm*rs**3*fn(y)/alpha).to(u.Msun)
      
      return encmass
  
Fix error in enclosed mass calculation and fix unit handling.
## Code Before:

import numpy as np
from astropy import units as u
from scipy import special


def encmass(r, norm, rs, alpha, beta, gamma):
    
    """
    Enclosed mass profile of a generalised Hernquist model.
    
    INPUTS
      r     : radial variable (requires unit)
      norm  : normalisation (requires unit)
      rs    : scale radius of model (requires unit)
      alpha : sharpness of transition between inner and outer
      beta  : outer logarithmic slope
      gamma : inner logarithmic slope
    """
    
    a = (3.-gamma)/alpha
    b = (gamma-beta)/alpha
    y = (r/rs)**alpha
    
    fn = lambda x: x**a * special.hyp2f1(a, -b, 1+a, -x)/a
    
    encmass = (4*np.pi*norm*rs**3*fn(y.value)).to(u.Msun)
    
    return encmass

## Instruction:
Fix error in enclosed mass calculation and fix unit handling.
## Code After:

import numpy as np
from astropy import units as u
from scipy import special


def encmass(r, norm, rs, alpha, beta, gamma):
    
    """
    Enclosed mass profile of a generalised Hernquist model.
    
    INPUTS
      r     : radial variable (requires unit)
      norm  : normalisation (requires unit)
      rs    : scale radius of model (requires unit)
      alpha : sharpness of transition between inner and outer
      beta  : outer logarithmic slope
      gamma : inner logarithmic slope
    """
    
    a = (3.-gamma)/alpha
    b = (gamma-beta)/alpha
    y = ((r/rs).to(u.dimensionless_unscaled).value)**alpha
    
    fn = lambda x: x**a * special.hyp2f1(a, -b, 1+a, -x)/a
    
    encmass = (4*np.pi*norm*rs**3*fn(y)/alpha).to(u.Msun)
    
    return encmass
