class NormBeta1D(ArithmeticModel):

    def __init__(self, name='normbeta1d'):
        self.pos = Parameter(name, 'pos', 0)
        self.width = Parameter(name, 'width', 1, tinyval, hard_min=tinyval)
        self.index = Parameter(name, 'index', 2.5, 0.5, 1000, 0.5)
        self.ampl = Parameter(name, 'ampl', 1, 0)
        ArithmeticModel.__init__(self, name,
                                 (self.pos, self.width, self.index, self.ampl))

    def get_center(self):
        return (self.pos.val,)

    def set_center(self, pos, *args, **kwargs):
        self.pos.set(pos)

    def guess(self, dep, *args, **kwargs):
        ampl = guess_amplitude(dep, *args)
        pos = get_position(dep, *args)
        fwhm = guess_fwhm(dep, *args)
        param_apply_limits(pos, self.pos, **kwargs)
        norm = (fwhm['val']*numpy.sqrt(numpy.pi)*
                numpy.exp(lgam(self.index.val-0.5)-lgam(self.index.val)))
        for key in ampl.keys():
            ampl[key] *= norm
        param_apply_limits(ampl, self.ampl, **kwargs)

    @modelCacher1d
    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.nbeta1d(*args, **kwargs)

class Lorentz1D(ArithmeticModel):

    def __init__(self, name='lorentz1d'):
        self.fwhm = Parameter(name, 'fwhm', 10, 0, hard_min=0)
        self.pos = Parameter(name, 'pos', 1)
        self.ampl = Parameter(name, 'ampl', 1)
        ArithmeticModel.__init__(self, name,
                                 (self.fwhm, self.pos, self.ampl))

    def get_center(self):
        return (self.pos.val,)

    def set_center(self, pos, *args, **kwargs):
        self.pos.set(pos)

    def guess(self, dep, *args, **kwargs):
        pos = get_position(dep, *args)
        fwhm = guess_fwhm(dep, *args)
        param_apply_limits(pos, self.pos, **kwargs)
        param_apply_limits(fwhm, self.fwhm, **kwargs)

        norm = guess_amplitude(dep, *args)
        if fwhm != 10:
            aprime = norm['val']*self.fwhm.val*numpy.pi/2.
            ampl = {'val': aprime,
                    'min': aprime/_guess_ampl_scale,
                    'max': aprime*_guess_ampl_scale}
            param_apply_limits(ampl, self.ampl, **kwargs)
        else:
            param_apply_limits(norm, self.ampl, **kwargs)

    @modelCacher1d
    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.lorentz1d(*args, **kwargs)

class Lorentz2D(ArithmeticModel):

    def __init__(self, name='lorentz2d'):
        self.fwhm = Parameter(name, 'fwhm', 10, tinyval, hard_min=tinyval)
        self.xpos = Parameter(name, 'xpos', 0)
        self.ypos = Parameter(name, 'ypos', 0)
        self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999,
                               frozen=True)
        self.theta = Parameter(name, 'theta', 0, 0, 2*numpy.pi, -2*numpy.pi,
                               4*numpy.pi, 'radians',frozen=True)
        self.ampl = Parameter(name, 'ampl', 1)
        ArithmeticModel.__init__(self, name,
                                 (self.fwhm, self.xpos, self.ypos, self.ellip,
                                  self.theta, self.ampl))
        self.cache = 0

    def get_center(self):
        return (self.xpos.val, self.ypos.val)

    def set_center(self, xpos, ypos, *args, **kwargs):
        self.xpos.set(xpos)
        self.ypos.set(ypos)

    def guess(self, dep, *args, **kwargs):
        xpos, ypos = guess_position(dep, *args)
        norm = guess_amplitude2d(dep, *args)
        param_apply_limits(xpos, self.xpos, **kwargs)
        param_apply_limits(ypos, self.ypos, **kwargs)
        param_apply_limits(norm, self.ampl, **kwargs)


    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.lorentz2d(*args, **kwargs)

class HubbleReynolds(ArithmeticModel):

    def __init__(self, name='hubblereynolds'):
        self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
        self.xpos = Parameter(name, 'xpos', 0)
        self.ypos = Parameter(name, 'ypos', 0)
        self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
        self.theta = Parameter(name, 'theta', 0, 0, 2*numpy.pi, -2*numpy.pi,
                               4*numpy.pi, 'radians')
        self.ampl = Parameter(name, 'ampl', 1)
        ArithmeticModel.__init__(self, name,
                                 (self.r0, self.xpos, self.ypos, self.ellip,
                                  self.theta, self.ampl))
        self.cache = 0

    def get_center(self):
        return (self.xpos.val, self.ypos.val)

    def set_center(self, xpos, ypos, *args, **kwargs):
        self.xpos.set(xpos)
        self.ypos.set(ypos)

    def guess(self, dep, *args, **kwargs):
        xpos, ypos = guess_position(dep, *args)
        norm = guess_amplitude2d(dep, *args)
        rad = guess_radius(*args)
        param_apply_limits(xpos, self.xpos, **kwargs)
        param_apply_limits(ypos, self.ypos, **kwargs)
        param_apply_limits(norm, self.ampl, **kwargs)
        param_apply_limits(rad, self.r0, **kwargs)


    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.hr(*args, **kwargs)

class Beta1D(ArithmeticModel):

    def __init__(self, name='beta1d'):
        self.r0 = Parameter(name, 'r0', 1, tinyval, hard_min=tinyval)
        self.beta = Parameter(name, 'beta', 1, 1e-05, 10, 1e-05, 10)
        self.xpos = Parameter(name, 'xpos', 0, 0, frozen=True)
        self.ampl = Parameter(name, 'ampl', 1, 0)
        ArithmeticModel.__init__(self, name,
                                 (self.r0, self.beta, self.xpos, self.ampl))

    def get_center(self):
        return (self.xpos.val,)

    def set_center(self, xpos, *args, **kwargs):
        self.xpos.set(xpos)

    def guess(self, dep, *args, **kwargs):
        pos = get_position(dep, *args)
        param_apply_limits(pos, self.xpos, **kwargs)

        ref = guess_reference(self.r0.min, self.r0.max, *args)
        param_apply_limits(ref, self.r0, **kwargs)

        norm = guess_amplitude_at_ref(self.r0.val, dep, *args)
        param_apply_limits(norm, self.ampl, **kwargs)


    @modelCacher1d
    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.beta1d(*args, **kwargs)

 def __init__(self, name='hubblereynolds'):
     self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
     self.xpos = Parameter(name, 'xpos', 0)
     self.ypos = Parameter(name, 'ypos', 0)
     self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
     self.theta = Parameter(name, 'theta', 0, 0, 2*numpy.pi, -2*numpy.pi,
                            4*numpy.pi, 'radians')
     self.ampl = Parameter(name, 'ampl', 1)
     ArithmeticModel.__init__(self, name,
                              (self.r0, self.xpos, self.ypos, self.ellip,
                               self.theta, self.ampl))
     self.cache = 0

class Synchrotron(ArithmeticModel):
    def __init__(self,name='IC'):
        self.index   = Parameter(name, 'index', 2.0, min=-10, max=10)
        self.ref     = Parameter(name, 'ref', 20, min=0, frozen=True, units='TeV')
        self.ampl    = Parameter(name, 'ampl', 1, min=0, max=1e60, hard_max=1e100, units='1e30/eV')
        self.cutoff  = Parameter(name, 'cutoff', 0.0, min=0,frozen=True, units='TeV')
        self.beta    = Parameter(name, 'beta', 1, min=0, max=10, frozen=True)
        self.B       = Parameter(name, 'B', 1, min=0, max=10, frozen=True, units='G')
        self.verbose = Parameter(name, 'verbose', 0, min=0, frozen=True)
        ArithmeticModel.__init__(self,name,(self.index,self.ref,self.ampl,self.cutoff,self.beta,self.B,self.verbose))
        self._use_caching = True
        self.cache = 10

    def guess(self,dep,*args,**kwargs):
        # guess normalization from total flux
        xlo,xhi=args
        model=self.calc([p.val for p in self.pars],xlo,xhi)
        modflux=trapz_loglog(model,xlo)
        obsflux=trapz_loglog(dep*(xhi-xlo),xlo)
        self.ampl.set(self.ampl.val*obsflux/modflux)

    @modelCacher1d
    def calc(self,p,x,xhi=None):

        index,ref,ampl,cutoff,beta,B,verbose = p

        # Sherpa provides xlo, xhi in KeV, we merge into a single array if bins required
        if xhi is None:
            outspec = x * u.keV
        else:
            outspec = _mergex(x,xhi) * u.keV

        if cutoff == 0.0:
            pdist = models.PowerLaw(ampl * 1e30 * u.Unit('1/eV'), ref * u.TeV, index)
        else:
            pdist = models.ExponentialCutoffPowerLaw(ampl * 1e30 * u.Unit('1/eV'),
                    ref * u.TeV, index, cutoff * u.TeV, beta=beta)

        sy = models.Synchrotron(pdist, B=B*u.G,
                log10gmin=5, log10gmax=10, ngamd=50)

        model = sy.flux(outspec, distance=1*u.kpc).to('1/(s cm2 keV)')

        # Do a trapz integration to obtain the photons per bin
        if xhi is None:
            photons = (model * outspec).to('1/(s cm2)').value
        else:
            photons = trapz_loglog(model,outspec,intervals=True).to('1/(s cm2)').value

        if verbose:
            print(self.thawedpars, trapz_loglog(outspec*model,outspec).to('erg/(s cm2)'))

        return photons

 def __init__(self,name='pp'):
     self.index   = Parameter(name , 'index'   , 2.1 , min=-10 , max=10)
     self.ref     = Parameter(name , 'ref'     , 60  , min=0   , frozen=True  , units='TeV')
     self.ampl    = Parameter(name , 'ampl'    , 100 , min=0   , max=1e60     , hard_max=1e100 , units='1e30/eV')
     self.cutoff  = Parameter(name , 'cutoff'  , 0   , min=0   , frozen=True  , units='TeV')
     self.beta    = Parameter(name , 'beta'    , 1   , min=0   , max=10       , frozen=True)
     self.nh      = Parameter(name , 'nH'      , 1   , min=0   , frozen=True  , units='1/cm3')
     self.verbose = Parameter(name , 'verbose' , 0   , min=0   , frozen=True)
     ArithmeticModel.__init__(self,name,(self.index,self.ref,self.ampl,self.cutoff,self.beta,self.nh,self.verbose))
     self._use_caching = True
     self.cache = 10

class HubbleReynolds(ArithmeticModel):
    """Two-dimensional Hubble-Reynolds model.

    Attributes
    ----------
    r0
        The core radius.
    xpos
        The center of the model on the x0 axis.
    ypos
        The center of the model on the x1 axis.
    ellip
        The ellipticity of the model.
    theta
        The angle of the major axis. It is in radians, measured
        counter-clockwise from the X0 axis (i.e. the line X1=0).
    ampl
        The amplitude refers to the maximum peak of the model.

    See Also
    --------
    Beta2D, DeVaucouleurs2D, Lorentz2D, Sersic2D

    Notes
    -----
    The functional form of the model for points is::

        f(x0,x1) = ampl / (1 + r(x0,x1))^2

        r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
                     -------------------------------------------
                                  r0^2 * (1-ellip)^2

        xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)

        yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)

    The grid version is evaluated by adaptive multidimensional
    integration scheme on hypercubes using cubature rules, based
    on code from HIntLib ([1]_) and GSL ([2]_).

    References
    ----------

    .. [1] HIntLib - High-dimensional Integration Library
           http://mint.sbg.ac.at/HIntLib/

    .. [2] GSL - GNU Scientific Library
           http://www.gnu.org/software/gsl/

    """


    def __init__(self, name='hubblereynolds'):
        self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
        self.xpos = Parameter(name, 'xpos', 0)
        self.ypos = Parameter(name, 'ypos', 0)
        self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
        self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
                               4*numpy.pi, 'radians')
        self.ampl = Parameter(name, 'ampl', 1)
        ArithmeticModel.__init__(self, name,
                                 (self.r0, self.xpos, self.ypos, self.ellip,
                                  self.theta, self.ampl))
        self.cache = 0

    def get_center(self):
        return (self.xpos.val, self.ypos.val)

    def set_center(self, xpos, ypos, *args, **kwargs):
        self.xpos.set(xpos)
        self.ypos.set(ypos)

    def guess(self, dep, *args, **kwargs):
        xpos, ypos = guess_position(dep, *args)
        norm = guess_amplitude2d(dep, *args)
        rad = guess_radius(*args)
        param_apply_limits(xpos, self.xpos, **kwargs)
        param_apply_limits(ypos, self.ypos, **kwargs)
        param_apply_limits(norm, self.ampl, **kwargs)
        param_apply_limits(rad, self.r0, **kwargs)


    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.hr(*args, **kwargs)

class test_parameter(SherpaTestCase):

    def setUp(self):
        self.p = Parameter('model', 'name', 0, -10, 10, -100, 100, 'units')
        self.afp = Parameter('model', 'name', 0, alwaysfrozen=True)

    def test_name(self):
        self.assertEqual(self.p.modelname, 'model')
        self.assertEqual(self.p.name, 'name')
        self.assertEqual(self.p.fullname, 'model.name')

    def test_alwaysfrozen(self):
        self.assertTrue(self.afp.frozen)
        self.afp.frozen = True
        self.assertTrue(self.afp.frozen)
        self.afp.freeze()
        self.assertTrue(self.afp.frozen)
        self.assertRaises(ParameterErr, self.afp.thaw)
        self.assertRaises(ParameterErr, setattr, self.afp, 'frozen', 0)

    def test_readonly_attributes(self):
        self.assertEqual(self.p.alwaysfrozen, False)
        self.assertRaises(AttributeError, setattr, self.p, 'alwaysfrozen', 1)
        self.assertEqual(self.p.hard_min, -100.0)
        self.assertRaises(AttributeError, setattr, self.p, 'hard_min', -1000)
        self.assertEqual(self.p.hard_max, 100.0)
        self.assertRaises(AttributeError, setattr, self.p, 'hard_max', 1000)

    def test_val(self):
        self.p.val = -7
        self.assertEqual(self.p.val, -7)
        self.assertTrue(type(self.p.val) is SherpaFloat)
        self.assertRaises(ValueError, setattr, self.p, 'val', 'ham')
        self.assertRaises(ParameterErr, setattr, self.p, 'val', -101)
        self.assertRaises(ParameterErr, setattr, self.p, 'val', 101)

    def test_min_max(self):
        for attr, sign in (('min', -1), ('max', 1)):
            setattr(self.p, attr, sign * 99)
            val = getattr(self.p, attr)
            self.assertEqual(val, sign * 99)
            self.assertTrue(type(val) is SherpaFloat)
            self.assertRaises(ValueError, setattr, self.p, attr, 'ham')
            self.assertRaises(ParameterErr, setattr, self.p, attr, -101)
            self.assertRaises(ParameterErr, setattr, self.p, attr, 101)

    def test_frozen(self):
        self.p.frozen = 1.0
        self.assertTrue(self.p.frozen is True)
        self.p.frozen = []
        self.assertTrue(self.p.frozen is False)
        self.assertRaises(TypeError, setattr, self.p.frozen, arange(10))
        self.p.link = self.afp
        self.assertTrue(self.p.frozen is True)
        self.p.link = None
        self.p.freeze()
        self.assertTrue(self.p.frozen is True)
        self.p.thaw()
        self.assertTrue(self.p.frozen is False)

    def test_link(self):
        self.p.link = None
        self.assertTrue(self.p.link is None)
        self.assertNotEqual(self.p.val, 17.3)
        self.afp.val = 17.3
        self.p.link = self.afp
        self.assertEqual(self.p.val, 17.3)
        self.p.unlink()
        self.assertTrue(self.p.link is None)
        self.assertRaises(ParameterErr, setattr, self.afp, 'link', self.p)
        self.assertRaises(ParameterErr, setattr, self.p, 'link', 3)
        self.assertRaises(ParameterErr, setattr, self.p, 'link',
                          3 * self.p + 2)

    def test_iter(self):
        for part in self.p:
            self.assertTrue(part is self.p)

class Beta2D(RegriddableModel2D):
    """Two-dimensional beta model function.

    The beta model is a Lorentz model with a varying power law.

    Attributes
    ----------
    r0
        The core radius.
    xpos
        X0 axis coordinate of the model center (position of the peak).
    ypos
        X1 axis coordinate of the model center (position of the peak).
    ellip
        The ellipticity of the model.
    theta
        The angle of the major axis. It is in radians, measured
        counter-clockwise from the X0 axis (i.e. the line X1=0).
    ampl
        The model value at the peak position (xpos, ypos).
    alpha
        The power-law slope of the profile at large radii.

    See Also
    --------
    Beta1D, DeVaucouleurs2D, HubbleReynolds, Lorentz2D, Sersic2D

    Notes
    -----
    The functional form of the model for points is::

        f(x0,x1) = ampl * (1 + r(x0,x1)^2)^(-alpha)

        r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
                     -------------------------------------------
                                  r0^2 * (1-ellip)^2

        xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)

        yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)

    The grid version is evaluated by adaptive multidimensional
    integration scheme on hypercubes using cubature rules, based
    on code from HIntLib ([1]_) and GSL ([2]_).

    References
    ----------

    .. [1] HIntLib - High-dimensional Integration Library
           http://mint.sbg.ac.at/HIntLib/

    .. [2] GSL - GNU Scientific Library
           http://www.gnu.org/software/gsl/

    """

    def __init__(self, name='beta2d'):
        self.r0 = Parameter(name, 'r0', 10, tinyval, hard_min=tinyval)
        self.xpos = Parameter(name, 'xpos', 0)
        self.ypos = Parameter(name, 'ypos', 0)
        self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999,
                               frozen=True)
        self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
                               4*numpy.pi, 'radians', True)
        self.ampl = Parameter(name, 'ampl', 1)
        self.alpha = Parameter(name, 'alpha', 1, -10, 10)
        ArithmeticModel.__init__(self, name,
                                 (self.r0, self.xpos, self.ypos, self.ellip,
                                  self.theta, self.ampl, self.alpha))
        self.cache = 0

    def get_center(self):
        return (self.xpos.val, self.ypos.val)

    def set_center(self, xpos, ypos, *args, **kwargs):
        self.xpos.set(xpos)
        self.ypos.set(ypos)

    def guess(self, dep, *args, **kwargs):
        xpos, ypos = guess_position(dep, *args)
        norm = guess_amplitude2d(dep, *args)
        rad = guess_radius(*args)
        param_apply_limits(xpos, self.xpos, **kwargs)
        param_apply_limits(ypos, self.ypos, **kwargs)
        param_apply_limits(norm, self.ampl, **kwargs)
        param_apply_limits(rad, self.r0, **kwargs)

    def calc(self, *args, **kwargs):
        kwargs['integrate'] = bool_cast(self.integrate)
        return _modelfcts.beta2d(*args, **kwargs)

 def setUp(self):
     self.p = Parameter('model', 'name', 0, -10, 10, -100, 100, 'units')
     self.afp = Parameter('model', 'name', 0, alwaysfrozen=True)

class Sersic2D(ArithmeticModel):
    """Two-dimensional Sersic model.

    This is a generalization of the ``DeVaucouleurs2D`` model,
    in which the exponent ``n`` can vary ([1]_, [2]_, and [3]_).

    Attributes
    ----------
    r0
        The core radius.
    xpos
        The center of the model on the x0 axis.
    ypos
        The center of the model on the x1 axis.
    ellip
        The ellipticity of the model.
    theta
        The angle of the major axis. It is in radians, measured
        counter-clockwise from the X0 axis (i.e. the line X1=0).
    ampl
        The amplitude refers to the maximum peak of the model.
    n
        The Sersic index (n=4 replicates the ``DeVaucouleurs2D``
        model).

    See Also
    --------
    Beta2D, DeVaucouleurs2D, HubbleReynolds, Lorentz2D

    Notes
    -----
    The functional form of the model for points is can be
    expressed as the following::

        f(x0,x1) = ampl * exp(-b(n) * (r(x0,x1)^(1/n) - 1))

            b(n) = 2 * n - 1 / 3 + 4 / (405 * n) + 46 / (25515 * n^2)

        r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
                     -------------------------------------------
                                  r0^2 * (1-ellip)^2

        xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)

        yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)

    The grid version is evaluated by adaptive multidimensional
    integration scheme on hypercubes using cubature rules, based
    on code from HIntLib ([4]_) and GSL ([5]_).

    References
    ----------

    .. [1] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html

    .. [2] Graham, A. & Driver, S., 2005, PASA, 22, 118
           http://adsabs.harvard.edu/abs/2005PASA...22..118G

    .. [3] Ciotti, L. & Bertin, G., A&A, 1999, 352, 447-451
           http://adsabs.harvard.edu/abs/1999A%26A...352..447C

    .. [4] HIntLib - High-dimensional Integration Library
           http://mint.sbg.ac.at/HIntLib/

    .. [5] GSL - GNU Scientific Library
           http://www.gnu.org/software/gsl/

    """

    def __init__(self, name='sersic2d'):
        self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
        self.xpos = Parameter(name, 'xpos', 0)
        self.ypos = Parameter(name, 'ypos', 0)
        self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
        self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
                               4*numpy.pi, 'radians')
        self.ampl = Parameter(name, 'ampl', 1)
        self.n = Parameter(name,'n', 1, .1, 10, 0.01, 100, frozen=True )
        ArithmeticModel.__init__(self, name,
                                 (self.r0, self.xpos, self.ypos, self.ellip,
                                  self.theta, self.ampl, self.n))
        self.cache = 0

    def get_center(self):
        return (self.xpos.val, self.ypos.val)

    def set_center(self, xpos, ypos, *args, **kwargs):
        self.xpos.set(xpos)
        self.ypos.set(ypos)

    def guess(self, dep, *args, **kwargs):
        xpos, ypos = guess_position(dep, *args)
        norm = guess_amplitude2d(dep, *args)
        rad = guess_radius(*args)
        param_apply_limits(xpos, self.xpos, **kwargs)
        param_apply_limits(ypos, self.ypos, **kwargs)
        param_apply_limits(norm, self.ampl, **kwargs)
        param_apply_limits(rad, self.r0, **kwargs)


#.........这里部分代码省略.........

galabso.nH.freeze()
galabso.nH.val = props['nhgal'] / 1e22

print('freezing background params')
for p in get_bkg_model(id).pars: 
	p.freeze()
print(get_model(id))
if id2:
	for p in get_bkg_model(id2).pars: 
		p.freeze()
	print(get_model(id2))

print('creating prior functions...')
srclevel = Parameter('src', 'level', numpy.log10(sphere.norm.val), -8, 3, -8, 3)
srclevel2 = Parameter('src2', 'level', numpy.log10(sphere2.norm.val), -8, 3, -8, 3)
srcnh = Parameter('src', 'nh', numpy.log10(sphere.nh.val)+22, 20, 26, 20, 26)
galnh = galabso.nH.val

sphere.norm = 10**srclevel
sphere2.norm = 10**srclevel2
sphere.nh = 10**(srcnh - 22)
sphere2.nh = 10**(srcnh - 22)
galabso.nH = 10**(galnh - 22)

priors = []
parameters = [srclevel, sphere.phoindex, srcnh]

import bxa.sherpa as bxa
priors += [bxa.create_uniform_prior_for(srclevel)]
priors += [bxa.create_gaussian_prior_for(sphere.phoindex, 1.95, 0.15)]
priors += [bxa.create_uniform_prior_for(srcnh)]

class Beta1D(ArithmeticModel):
    """One-dimensional beta model function.

    The beta model is a Lorentz model with a varying power law.

    Attributes
    ----------
    r0
        The core radius.
    beta
        This parameter controls the slope of the profile at large
        radii.
    xpos
        The reference point of the profile. This is frozen by default.
    ampl
        The amplitude refers to the maximum value of the model, at
        x = xpos.

    See Also
    --------
    Beta2D, Lorentz1D, NormBeta1D

    Notes
    -----
    The functional form of the model for points is::

        f(x) = ampl * (1 + ((x - xpos) / r0)^2)^(0.5 - 3 * beta)

    The grid version is evaluated by numerically intgerating the
    function over each bin using a non-adaptive Gauss-Kronrod scheme
    suited for smooth functions [1]_, falling over to a simple
    trapezoid scheme if this fails.

    References
    ----------

    .. [1] https://www.gnu.org/software/gsl/manual/html_node/QNG-non_002dadaptive-Gauss_002dKronrod-integration.html

    """

    def __init__(self, name='beta1d'):
        self.r0 = Parameter(name, 'r0', 1, tinyval, hard_min=tinyval)
        self.beta = Parameter(name, 'beta', 1, 1e-05, 10, 1e-05, 10)
        self.xpos = Parameter(name, 'xpos', 0, 0, frozen=True)
        self.ampl = Parameter(name, 'ampl', 1, 0)
        ArithmeticModel.__init__(self, name,
                                 (self.r0, self.beta, self.xpos, self.ampl))

    def get_center(self):
        return (self.xpos.val,)

    def set_center(self, xpos, *args, **kwargs):
        self.xpos.set(xpos)

    def guess(self, dep, *args, **kwargs):
        pos = get_position(dep, *args)
        param_apply_limits(pos, self.xpos, **kwargs)

        ref = guess_reference(self.r0.min, self.r0.max, *args)
        param_apply_limits(ref, self.r0, **kwargs)

        norm = guess_amplitude_at_ref(self.r0.val, dep, *args)
        param_apply_limits(norm, self.ampl, **kwargs)


    @modelCacher1d
    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.beta1d(*args, **kwargs)

class Lorentz1D(ArithmeticModel):
    """One-dimensional normalized Lorentz model function.

    Attributes
    ----------
    fwhm
        The full-width half maximum of the line.
    pos
        The center of the line.
    ampl
        The amplitude refers to the integral of the model.

    See Also
    --------
    Beta1D, NormBeta1D

    Notes
    -----
    The functional form of the model for points is::

        f(x) =                A * fwhm
               --------------------------------------
               2 * pi * (0.25 * fwhm^2 + (x - pos)^2)

           A = ampl / integral f(x) dx

    and for an integrated grid it is the integral of this over
    the bin.
    """

    def __init__(self, name='lorentz1d'):
        self.fwhm = Parameter(name, 'fwhm', 10, 0, hard_min=0)
        self.pos = Parameter(name, 'pos', 1)
        self.ampl = Parameter(name, 'ampl', 1)
        ArithmeticModel.__init__(self, name,
                                 (self.fwhm, self.pos, self.ampl))

    def get_center(self):
        return (self.pos.val,)

    def set_center(self, pos, *args, **kwargs):
        self.pos.set(pos)

    def guess(self, dep, *args, **kwargs):
        pos = get_position(dep, *args)
        fwhm = guess_fwhm(dep, *args)
        param_apply_limits(pos, self.pos, **kwargs)
        param_apply_limits(fwhm, self.fwhm, **kwargs)

        norm = guess_amplitude(dep, *args)
        if fwhm != 10:
            aprime = norm['val']*self.fwhm.val*numpy.pi/2.
            ampl = {'val': aprime,
                    'min': aprime/_guess_ampl_scale,
                    'max': aprime*_guess_ampl_scale}
            param_apply_limits(ampl, self.ampl, **kwargs)
        else:
            param_apply_limits(norm, self.ampl, **kwargs)

    @modelCacher1d
    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.lorentz1d(*args, **kwargs)

class NormBeta1D(ArithmeticModel):
    """One-dimensional normalized beta model function.

    This is the same model as the ``Beta1D`` model but with a
    different slope parameter and normalisation.

    Attributes
    ----------
    pos
        The center of the line.
    w
        The line width.
    alpha
        The slope of the profile at large radii.
    ampl
        The amplitude refers to the integral of the model.

    See Also
    --------
    Beta1D, Lorentz1D

    Notes
    -----
    The functional form of the model for points is::

        f(x) = A * (1 + ((x - pos) / w)^2)^(-alpha)

           A = ampl / integral f(x) dx

    The grid version is evaluated by numerically intgerating the
    function over each bin using a non-adaptive Gauss-Kronrod scheme
    suited for smooth functions [1]_, falling over to a simple
    trapezoid scheme if this fails.

    References
    ----------

    .. [1] https://www.gnu.org/software/gsl/manual/html_node/QNG-non_002dadaptive-Gauss_002dKronrod-integration.html

    """

    def __init__(self, name='normbeta1d'):
        self.pos = Parameter(name, 'pos', 0)
        self.width = Parameter(name, 'width', 1, tinyval, hard_min=tinyval)
        self.index = Parameter(name, 'index', 2.5, 0.5, 1000, 0.5)
        self.ampl = Parameter(name, 'ampl', 1, 0)
        ArithmeticModel.__init__(self, name,
                                 (self.pos, self.width, self.index, self.ampl))

    def get_center(self):
        return (self.pos.val,)

    def set_center(self, pos, *args, **kwargs):
        self.pos.set(pos)

    def guess(self, dep, *args, **kwargs):
        ampl = guess_amplitude(dep, *args)
        pos = get_position(dep, *args)
        fwhm = guess_fwhm(dep, *args)
        param_apply_limits(pos, self.pos, **kwargs)
        norm = (fwhm['val']*numpy.sqrt(numpy.pi)*
                numpy.exp(lgam(self.index.val-0.5)-lgam(self.index.val)))
        for key in ampl.keys():
            ampl[key] *= norm
        param_apply_limits(ampl, self.ampl, **kwargs)

    @modelCacher1d
    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.nbeta1d(*args, **kwargs)

class DeVaucouleurs2D(ArithmeticModel):
    """Two-dimensional de Vaucouleurs model.

    This is a formulation of the R^(1/4) law introduced by [1]_. It
    is a special case of the ``Sersic2D`` model with ``n=4``,
    as described in [2]_, [3]_, and [4]_.

    Attributes
    ----------
    r0
        The core radius.
    xpos
        The center of the model on the x0 axis.
    ypos
        The center of the model on the x1 axis.
    ellip
        The ellipticity of the model.
    theta
        The angle of the major axis. It is in radians, measured
        counter-clockwise from the X0 axis (i.e. the line X1=0).
    ampl
        The amplitude refers to the maximum peak of the model.

    See Also
    --------
    Beta2D, HubbleReynolds, Lorentz2D, Sersic2D

    Notes
    -----
    The model used is the same as the ``Sersic2D`` model with ``n=4``.

    References
    ----------

    .. [1] de Vaucouleurs G., 1948, Ann. d’Astroph. 11, 247
           http://adsabs.harvard.edu/abs/1948AnAp...11..247D

    .. [2] http://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html

    .. [3] Graham, A. & Driver, S., 2005, PASA, 22, 118
           http://adsabs.harvard.edu/abs/2005PASA...22..118G

    .. [4] Ciotti, L. & Bertin, G., A&A, 1999, 352, 447-451
           http://adsabs.harvard.edu/abs/1999A%26A...352..447C

    """

    def __init__(self, name='devaucouleurs2d'):
        self.r0 = Parameter(name, 'r0', 10, 0, hard_min=0)
        self.xpos = Parameter(name, 'xpos', 0)
        self.ypos = Parameter(name, 'ypos', 0)
        self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999)
        self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
                               4*numpy.pi, 'radians')
        self.ampl = Parameter(name, 'ampl', 1)
        ArithmeticModel.__init__(self, name,
                                 (self.r0, self.xpos, self.ypos, self.ellip,
                                  self.theta, self.ampl))
        self.cache = 0

    def get_center(self):
        return (self.xpos.val, self.ypos.val)

    def set_center(self, xpos, ypos, *args, **kwargs):
        self.xpos.set(xpos)
        self.ypos.set(ypos)

    def guess(self, dep, *args, **kwargs):
        xpos, ypos = guess_position(dep, *args)
        norm = guess_amplitude2d(dep, *args)
        rad = guess_radius(*args)
        param_apply_limits(xpos, self.xpos, **kwargs)
        param_apply_limits(ypos, self.ypos, **kwargs)
        param_apply_limits(norm, self.ampl, **kwargs)
        param_apply_limits(rad, self.r0, **kwargs)


    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.devau(*args, **kwargs)

class Lorentz2D(ArithmeticModel):
    """Two-dimensional un-normalised Lorentz function.

    Attributes
    ----------
    fwhm
        The full-width half maximum.
    xpos
        The center of the model on the x0 axis.
    ypos
        The center of the model on the x1 axis.
    ellip
        The ellipticity of the model.
    theta
        The angle of the major axis. It is in radians, measured
        counter-clockwise from the X0 axis (i.e. the line X1=0).
    ampl
        The amplitude refers to the maximum peak of the model.

    See Also
    --------
    Beta1D, DeVaucouleurs2D, HubbleReynolds, Lorentz1D, Sersic2D

    Notes
    -----
    The functional form of the model for points is::

        f(x0,x1) = ampl / (1 + 4 * r(x0,x1)^2)

        r(x0,x1)^2 = xoff(x0,x1)^2 * (1-ellip)^2 + yoff(x0,x1)^2
                     -------------------------------------------
                                 fwhm^2 * (1-ellip)^2

        xoff(x0,x1) = (x0 - xpos) * cos(theta) + (x1 - ypos) * sin(theta)

        yoff(x0,x1) = (x1 - ypos) * cos(theta) - (x0 - xpos) * sin(theta)

    and for an integrated grid it is the integral of this over
    the bin.
    """

    def __init__(self, name='lorentz2d'):
        self.fwhm = Parameter(name, 'fwhm', 10, tinyval, hard_min=tinyval)
        self.xpos = Parameter(name, 'xpos', 0)
        self.ypos = Parameter(name, 'ypos', 0)
        self.ellip = Parameter(name, 'ellip', 0, 0, 0.999, 0, 0.9999,
                               frozen=True)
        self.theta = Parameter(name, 'theta', 0, -2*numpy.pi, 2*numpy.pi, -2*numpy.pi,
                               4*numpy.pi, 'radians',frozen=True)
        self.ampl = Parameter(name, 'ampl', 1)
        ArithmeticModel.__init__(self, name,
                                 (self.fwhm, self.xpos, self.ypos, self.ellip,
                                  self.theta, self.ampl))
        self.cache = 0

    def get_center(self):
        return (self.xpos.val, self.ypos.val)

    def set_center(self, xpos, ypos, *args, **kwargs):
        self.xpos.set(xpos)
        self.ypos.set(ypos)

    def guess(self, dep, *args, **kwargs):
        xpos, ypos = guess_position(dep, *args)
        norm = guess_amplitude2d(dep, *args)
        param_apply_limits(xpos, self.xpos, **kwargs)
        param_apply_limits(ypos, self.ypos, **kwargs)
        param_apply_limits(norm, self.ampl, **kwargs)


    def calc(self, *args, **kwargs):
        kwargs['integrate']=bool_cast(self.integrate)
        return _modelfcts.lorentz2d(*args, **kwargs)

class InverseCompton(ArithmeticModel):
    def __init__(self,name='IC'):
        self.index   = Parameter(name, 'index', 2.0, min=-10, max=10)
        self.ref     = Parameter(name, 'ref', 20, min=0, frozen=True, units='TeV')
        self.ampl    = Parameter(name, 'ampl', 1, min=0, max=1e60, hard_max=1e100, units='1e30/eV')
        self.cutoff  = Parameter(name, 'cutoff', 0.0, min=0,frozen=True, units='TeV')
        self.beta    = Parameter(name, 'beta', 1, min=0, max=10, frozen=True)
        self.TFIR    = Parameter(name, 'TFIR', 70, min=0, frozen=True, units='K')
        self.uFIR    = Parameter(name, 'uFIR', 0.0, min=0, frozen=True, units='eV/cm3') # 0.2eV/cm3 typical in outer disk
        self.TNIR    = Parameter(name, 'TNIR', 3800, min=0, frozen=True, units='K')
        self.uNIR    = Parameter(name, 'uNIR', 0.0, min=0, frozen=True, units='eV/cm3') # 0.2eV/cm3 typical in outer disk
        self.verbose = Parameter(name, 'verbose', 0, min=0, frozen=True)
        ArithmeticModel.__init__(self,name,(self.index,self.ref,self.ampl,self.cutoff,self.beta,
            self.TFIR, self.uFIR, self.TNIR, self.uNIR, self.verbose))
        self._use_caching = True
        self.cache = 10

    def guess(self,dep,*args,**kwargs):
        # guess normalization from total flux
        xlo,xhi=args
        model=self.calc([p.val for p in self.pars],xlo,xhi)
        modflux=trapz_loglog(model,xlo)
        obsflux=trapz_loglog(dep*(xhi-xlo),xlo)
        self.ampl.set(se