# Path: DiscEvolution/constants.py
# G = 1.
# AU = 1.496e13


import numpy as np
from scipy import optimize
from .constants import AU, sig_H2, m_H, yr, Msun
        return self._FUV

    """Methods to determine global properties of a viscous accretion disc"""
    def Rout(self, thresh=1e-5):
        """Determine the outer radius via density threshold"""
        notempty = self.Sigma_G > thresh
        notempty_cells = self.R_edge[1:][notempty]
        if np.size(notempty_cells>0):
            R_outer = notempty_cells[-1]
        else:
            R_outer = 0.0

        return R_outer

    def RC(self):
        """Fit an LBP profile to the disc and return the scale radius"""
        not_empty = (self.R < self.Rout())
        popt,pcov = optimize.curve_fit(LBP_profile,self.R[not_empty],np.log(self.Sigma_G[not_empty]),p0=[100,0.01],maxfev=5000)
        return popt[0]

    def Mtot(self):
        """Determine the total mass"""
        Re = self.R_edge * AU
        dA = np.pi * (Re[1:] ** 2 - Re[:-1] ** 2)
        dM_tot = self.Sigma * dA
        M_tot = np.sum(dM_tot)
        return M_tot

    def Mdot(self,viscous_velocity):
        """Determine the viscous accretion rate"""

# Path: DiscEvolution/constants.py
# G = 1.
# AU = 1.496e13


# Path: DiscEvolution/disc.py
import numpy as np
from scipy import optimize
from .constants import AU, sig_H2, m_H, yr, Msun

        return self._FUV

    """Methods to determine global properties of a viscous accretion disc"""
    def Rout(self, thresh=1e-5):
        """Determine the outer radius via density threshold"""
        notempty = self.Sigma_G > thresh
        notempty_cells = self.R_edge[1:][notempty]
        if np.size(notempty_cells>0):
            R_outer = notempty_cells[-1]
        else:
            R_outer = 0.0

        return R_outer

    def RC(self):
        """Fit an LBP profile to the disc and return the scale radius"""
        not_empty = (self.R < self.Rout())
        popt,pcov = optimize.curve_fit(LBP_profile,self.R[not_empty],np.log(self.Sigma_G[not_empty]),p0=[100,0.01],maxfev=5000)
        return popt[0]

    def Mtot(self):
        """Determine the total mass"""
        Re = self.R_edge * AU
        dA = np.pi * (Re[1:] ** 2 - Re[:-1] ** 2)
        dM_tot = self.Sigma * dA
        M_tot = np.sum(dM_tot)
        return M_tot

    def Mdot(self,viscous_velocity):
        """Determine the viscous accretion rate"""

        M_visc_out = 2*np.pi * self.R[0] * self.Sigma[0] * viscous_velocity * (AU**2)