DougBurke
commented
2 years ago I've stuck ni sme debugging statements so that the following
import numpy as np
from bhjet import BHJet
egrid = np.arange(0.1, 2.0, 0.1)
elo = egrid[:-1]
ehi = egrid[1:]
mdl = BHJet()
print(mdl)
y = mdl(elo, ehi)will create
bhjet
Param Type Value Min Max Units
----- ---- ----- --- --- -----
bhjet.mbh frozen 1e+06 3 3e+10 msun
bhjet.incl frozen 40 2 80 deg
bhjet.dkpc frozen 3000 1 1e+06 kpc
bhjet.redshift frozen 0 0 7
bhjet.jetrat thawed 0.005 1e-07 1 L_edd
bhjet.r_0 thawed 20 2 80 r_g
bhjet.z_diss thawed 1000 50 500000 r_g
bhjet.z_acc frozen 1000 300 500000 r_g
bhjet.z_max frozen 1e+07 100000 1e+09 r_g
bhjet.t_e thawed 100 10 2000 keV
bhjet.f_nth frozen 0.1 0.05 0.95
bhjet.f_pl frozen 0 0 10
bhjet.pspec thawed 2 1.5 3
bhjet.f_heat frozen 1 1 50
bhjet.f_b frozen 0.1 0.001 1
bhjet.f_sc thawed 1.5e-07 1e-09 0.1
bhjet.p_beta thawed 0.05 0.001 1
bhjet.sig_acc frozen 0.1 0.01 1
bhjet.l_disk thawed 0.01 0.0001 1 L_edd
bhjet.r_in thawed 1 1 200 r_g
bhjet.r_out frozen 1000 10 100000 r_g
bhjet.compar1 frozen 1000 3 1e+06
bhjet.compar2 frozen 0 0 1
bhjet.compar3 frozen 3e-10 0 1
bhjet.compsw frozen 0 -3.40282e+38 3.40282e+38
bhjet.velsw frozen 1 1 25
bhjet.infosw frozen 0 -3.40282e+38 3.40282e+38
bhjet.norm frozen 1 0 1e+24
# DBG: user energy grid ..
energy[0] = 0.1
energy[1] = 0.2
energy[2] = 0.3
energy[3] = 0.4
energy[4] = 0.5
...
energy[17] = 1.8
energy[18] = 1.9
# DBG: egrid for jetmain ..
egrid[0] = 1e-11
egrid[1] = 1.1749e-11
egrid[2] = 1.38038e-11
egrid[3] = 1.62181e-11
egrid[4] = 1.90546e-11
..
egrid[295] = 4.46684e+09
egrid[296] = 5.24807e+09
egrid[297] = 6.16595e+09
egrid[298] = 7.24436e+09
gsl: interp.c:83: ERROR: x values must be strictly increasing
Default GSL error handler invoked.
Aborted (core dumped)which shows the inputs to jetmain(egrid, NBINS, const_cast<double *>(parameter), energ, phot); where it fails with the GSL error. The first range (0.1 to 1.9, which is in keV) are the values of the energ array and the 299 or so bins refer to the egrid input.
matteolucchini1
The primary change is to allow the BHJet model to be used from Sherpa. There's also some support for building an XSPEC interface [*], but I have not got the time/energy to do the useful part of getting the build to work (I don't think it would be hard).
There are a few minor tweaks to the README to make it read better on GitHub when parsed as markdown, and the use of gsl-config rather than hard-coding the values when building the BHJet code.
The default parameter values and ranges are based on bhjet.sl, as pointed out by Sera, but all mistakes are my own.
Is there a
Input/ip.datfile that could be included or referred to let users check that they get sensible results.[*] the Sherpa interface is built on Sherpa's XSPEC interface, so the
wrap_xspec.cppfile that I add which should be useful for XSPEC is also useful for Sherpa.Sherpa installation
This is added to the
BHJet/README.mdfile, but installation requires CIAO 4.14, preferably installed with conda, and then you can sayThen you can load it into Sherpa and use it as a model with
Implementation note
So
wrap_xspec.cppprovides an "XSPEC-compatible" interface using the C (rather than FORTRAN or C++) API - see https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixLocal.htmlThis wrapper routine, wonderfully named
xspec_jetinterp, is meant to be close to the ISISbhjet_fitroutine - i.e. it does both the call tojetmainand then the call tojetinterp, rather than having a Python-level routine analogous tobhfit_jet. This is less open to user tweaks, but easier to handle (and I think required if you wanted to use the model in XSPEC; in Sherpa we could use a similar approach to ISIS).Unfortunately as I'm not an ISIS user then I appear to have made a mistake somewhere (see next section). In the S-Lang/ISIS call to
bhjet_fit(), are theloandhiarguments always in Angstroms?Current problems
Importing the Sherpa model "directly" (a "lower-level" interface than users are expected to use) shows that I have not correctly handled the call to both
jetmainand thenjetinterp, but I haven't got the time to track this down just yetFor the really low-level interface, first get the default parameter values as I'm too-lazy to type them in myself: note that we have 28 parameters because the last-one is the normalization.
Now we can access the "C++ code" (essentially this is equivalent to slirp creating a S-Lang functino that calls the C++ code):
This routtine takes three array/list arguments: the parameters, the low-edge, and the high-edge of the bins (it's done this way so I can just re-use the code we have for interfacing to X-SPEC models, but we could change it if need be):
So, am I calling it wrong? The lo/hi edges are assumed to be keV.