BroadLineRegions.jl

A fast and flexible toolkit for modeling the broad-line region (BLR) in Julia.

Installation

Julia

using Pkg
Pkg.add("BroadLineRegions")

Or install from the GitHub repo directly:

using Pkg
Pkg.add(url="https://github.com/kirklong/BroadLineRegions.jl")

To run the unit tests and verify you have installed everything properly, you can then run:

Pkg.test("BroadLineRegions")

Python

To access BroadLineRegions from within your Python installation, first you must install JuliaCall. After successfully installing JuliaCall you can add BroadLineRegions to your new Julia installation in Python with:

from juliacall import Main as jl
from juliacall import Pkg as jlPkg
jlPkg.add("BroadLineRegions") #or use the github link

Note that if you are using BroadLineRegions.jl from within Python, make sure that you prepend the commands with the corresponding name that you imported Julia with, e.g. the first example command below becomes:

from juliacall import Main as jl
import numpy as np
jl.seval("using BroadLineRegions") #load BLR into Main

# Use ASCII keyword aliases: Python normalizes some Unicode identifiers (e.g. θₒ).
mClouds = jl.BLR.cloudModel(100_000,mu=500.,beta=1.,F=0.5,theta_o=30/180*np.pi,gamma=1.,xi=1.,i=0.,
             I=jl.BLR.IsotropicIntensity,v=jl.BLR.vCircularCloud,rescale=1e-5,tau=0.0)

Quickstart

While this code is designed to be very flexible and modular such that you can implement your own bespoke models of the BLR more easily, two popular models of the BLR are included as default models which one can play with.

To generate a "cloud" model similar to that of Pancoast+ 2011 and 2014, use syntax like:

using BroadLineRegions #exports itself as both BroadLineRegions and the shorter acronym BLR
mClouds = BLR.cloudModel(100_000,μ=500.,β=1.,F=0.5,θₒ=30/180*π,γ=1.,ξ=1.,i=0.,
        I=BLR.IsotropicIntensity,v=BLR.vCircularCloud,rescale=1e-5,τ=0.0)

To generate a "disk-wind" model similar to that of Chiang and Murray 1996 and 1997 following the prescription laid out in Long+ 2023 and 2025 use syntax like:

mDisk = BLR.DiskWindModel(500.,5.,1.,30/180*π,nr=24,nϕ=48,scale=:log,
        f1=1.0,f2=1.0,f3=1.0,f4=1.0,reflect=false,τ=5.)

In case you do not like Unicode symbols, spelled out keyword aliases (for example mu, theta_o, and tau) are also available for all parameters (and can be mixed and matched with the Unicode symbols), i.e. for the cloud model example above we could also write:

mClouds = BLR.cloudModel(100_000,mu=500.,beta=1.,F=0.5,theta_o=30/180*pi,gamma=1.,xi=1.,i=0.,
        I=BLR.IsotropicIntensity,v=BLR.vCircularCloud,rescale=1e-5,tau=0.0)

See each function's documentation (either online or in the REPL with e.g. ?BLR.cloudModel) for a list of all parameters, their defaults, and available spelled out ASCII keyword aliases. Note that the internal code still uses largely Unicode names, so to access struct fields only the Unicode versions are supported (i.e. only m.camera.α works, and m.camera.alpha will produce an error).

If you prefer to start from the geometric basis of one of these models but modify how some internal portion is calculated, you can pass your own custom functions to them to enable this functionality. For example if you wanted to generate a geometry similar to the DiskWindModel example above but with a custom intensity prescription you could write something like:

myCustomIntensityFx(;r,level=1.0,kwargs...) = level./r #custom intensity that scales ∝ 1/r with some "level"
mDiskCustom = BLR.DiskWindModel(500.,5.,1.,30/180*π,nr=24,nϕ=48,scale=:log,
        I=myCustomIntensityFx,level=2.) #pass extra new parameter level to custom intensity function

The code is designed in such a way that by default coordinates are always in reference to a virtual "camera" imaging the system, and by default the inner struct we use to hold these points is called a ring (which acts as a "ring" of points on an image, from which we "raytrace" back to physical space at construction). You can of course go beyond these default recipes and define fully custom models by passing your own "rings" and your own "camera" system to the model, such that you can configure the code quite flexibly while still being able to interface with existing models. Below we show an example of this general process for a very simple custom BLR model consisting of just a single BLR "cloud":

α,β = (0.,-50.) #start by defining (in camera coordinates) a "cloud"
myCustomCamera = BLR.camera([α],[β],false) #define custom BLR "camera" struct (α and β must be vectors, because usually you have more than one point in your BLR...)

i = 30/180*π #suppose we view system at physical inclination of 30 deg

#"raytrace" from camera coordinates to physical coordinates, calculate r and ϕ of our BLR "cloud" in physical space
xRing = -β/cos(i); yRing = α #coordinate convention for "ring" plane defined such that ϕ=0 points towards observer
r = √(xRing^2 + yRing^2)
ϕ = atan(yRing, xRing)

myCustomRings = Vector{BLR.ring}(undef, 1) #allocate empty vector of rings, in this case of size 1 to hold our singular BLR "cloud"
myCustomRings[1] = BLR.ring(r=r, i=i, ϕ=ϕ, v=0.0, I=1.0, ΔA=1.) #fill our vector using the calculated parameters above, using the bare minimum required parameters (see full examples/documentation for all possibilities)

#create our custom BLR, which then can be combined with other BLR models and profiles/images can be calculated from as in other examples
mCustom = BLR.model(myCustomRings, Dict{Symbol,BLR.profile}(), myCustomCamera, [1]) #note we also pass empty dictionary of profiles (since none have been pre-calculated) and pass vector of submodel indices [1] as this model doesn't have any subcomponents. 

Note that BroadLineRegions.jl attempts to make sure you have correctly passed keyword arguments for functions, and will issue a warning if it detects a keyword argument that it does not know how to check. This warning would be issued in the above case, as level is not a parameter of any existing function. For custom callbacks, you can thus optionally register extra keywords to enable stricter validation and suppress this warning:

myCustomIntensityFx(; r, level=1.0, kwargs...) = level ./ r #a simple custom intensity scaling like 1/r
BLR.acceptedKeywords(::typeof(myCustomIntensityFx)) = (:level,)

Models can be combined simply by writing mCombined = m1 + m2. After combining, you can run quick occlusion/optical depth raytracing calculations through raytrace!, but note that this implementation is by default exceedingly rudimentary (no complicated radiative transfer is solved, instead emission is just added up proportionally to $e^{-\tau}$ at each point along a ray until it reaches a specified optical depth stopping point and after that it is assumed all remaining emission along the ray is blocked, and you must supply the expected $\tau$ and emission values when creating the submodels).

Generate profiles (i.e. line, phase, delay, or whatever else your heart desires) for models with syntax like:

m = mDisk #or mClouds, mCustom, however you have defined your model
p = BLR.getProfile(m,:line) #generate line profile with default parameters
BLR.setProfile!(m,p) #optionally store the profile in model data structure

These profiles are just simple estimates using weighted sums (rectangular Riemann sum integral approximations only), so smooth profiles can depend highly on the number of points, especially in the case of discrete "cloud" models.

A few default visualization recipes exist as well:

profiles = BLR.profile(m) #plot all of the model's stored profiles 
img = BLR.image(m,:I) #make an image of the model intensity, can pass any other parameter as well to "image" them
geometry = BLR.plot3d(m) #visualize the geometry of the system in a 3d plot, can also color points according to any parameter

Note that no functions are exported by default into the global namespace to prevent overlap with other modules, and you must prepend the module name to all methods to access them. BroadLineRegions.jl exports itself as both BroadLineRegions and BLR, so both of these prefixes are equivalent, i.e. BroadLineRegions.model == BLR.model. If adding this prefix annoys you, you can avoid this by specifying the functions you want added to the global namespace in your import/using statement. For example:

using BroadLineRegions: DiskWindModel, cloudModel
mDisk = DiskWindModel(params...) #instead of having to type BLR.DiskWindModel

For more detailed examples, see the Usage and Examples page.

Full documentation is available in the API section.

Remember that if you are using BroadLineRegions from within Python you will need to modify all the code examples with the corresponding name that you imported Julia with, as demonstrated in the cloudModel example in the Python installation section at the top of this page.

Performance tips

Starting Julia with multiple threads (julia -t N) speeds up disk-wind model construction, especially for larger models, but note that "cloud" models use a serial construction path by default to preserve RNG reproducibility. GPU support is now included for Metal (Apple) and CUDA devices, which greatly speed up model generation and observable calculations, especially when raytracing multicomponent models. You may also consider starting Julia with a sysimage to improve performance (see pure Julia example here, and instructions for enabling the functionality in juliacall here).

Referencing

If you find this code useful in your work, please cite it as:

lorem ipsum

Contributing

If you would like to contribute to the package, please open a pull request. For bug reports and feature requests, please open an issue.