() = evalfile("bhjet.sl");
variable loval = _A(10^[-11.0:13.0:0.07]);
variable hival = make_hi_grid(loval);
cache_fun("bhjet",loval,hival);
model("constant(Isis_Active_Dataset)*tbabs(1)*reflect(1,bhjet_cache(1))+gaussian(1)");

() = xspec_abund("wilm"); %// Wilms abundances, needed for tbnew model
() = xspec_xsect("vern"); %// Verner cross-sections, needed for tbnew model

load_par("new_model/HH/hadronic_pspec17.dat");
%load_par("new_model/HS/hadronic_pspec22.dat");
%load_par("new_model/LH/leptonic_pspec17.dat");
%load_par("new_model/LS/leptonic_pspec22.dat");
set_par_fun("reflect(1).cosIncl","cos(bhjet_cache(1).incl*PI/180.)"); 

ec;