Reverberation Mapping Worked Example

Whilst there are quick breakdowns of how the code works on the previous page and the bindata/tfpy wiki, here's an example fully worked-through from the very start to the very end. This example is for a response map of the Hβ line of a QSO.

.pf file

This is the base .pf file for the 100% ionising luminosity version. This is run as a standard python model, and two others templated off of it with the lum_agn and disk.mdot values increased/decreased by 10% are also run.

System_type(0=star,1=binary,2=agn)	2
Number.of.wind.components		1
Wind_type(0=SV,1=Sphere,2=Previous,3=Proga,4=Corona,5=knigge,6=thierry,7=yso,8=elvis,9=shell)	0
Coord.system(0=spherical,1=cylindrical,2=spherical_polar,3=cyl_var)	1
Wind.dim.in.x_or_r.direction	200
Wind.dim.in.z_or_theta.direction	200
disk.type(0=no.disk,1=standard.flat.disk,2=vertically.extended.disk)	1
disk.atmosphere(0=no,1=yes) 	0
Atomic_data				data/h10_hetop_lohe1_standard78
photons_per_cycle	 		1e8
Ionization_cycles			25
spectrum_cycles			40
Wind_ionization(0=on.the.spot,1=LTE,2=fixed,3=recalc_bb,5=recalc_pow,6=pairwise_bb,7=pairwise_pow)	9
Line_transfer(0=pure.abs,1=pure.scat,2=sing.scat,3=escape.prob,6=macro_atoms,7=macro_atoms+aniso.scattering)	7
Thermal_balance_options(0=everything.on,1=no.adiabatic)	0
Disk_radiation(y=1)			1
Wind_radiation(y=1)			1
QSO_BH_radiation(y=1)		1
Rad_type_for_disk(0=bb,1=models)_to_make_wind	0
Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_to_make_wind	3
mstar(msol)				1e9
rstar(cm)				8.85667e+14
disk.mdot(msol/yr)			5.0
Disk.illumination.treatment(0=no.rerad,1=high.albedo,2=thermalized.rerad,3=analytic)	0
Disk.temperature.profile(0=standard;1=readin)	0
disk.radmax(cm)			1e17
lum_agn(ergs/s)			1e45
agn_power_law_index			-0.9
geometry_for_pl_source(0=sphere,1=lamp_post)	0
Torus(0=no,1=yes)			0
wind.radmax(cm)			1e20
wind.t.init				1e5
wind.mdot(msol/yr)			5
sv.diskmin(wd_rad)			50
sv.diskmax(wd_rad)			100
sv.thetamin(deg)			70
sv.thetamax(deg)			82
sv.mdot_r_exponent			0.0
sv.v_infinity(in_units_of_vescape	1.0
sv.acceleration_length(cm)		1e19
sv.acceleration_exponent		0.6
filling_factor(1=smooth,<1=clumped)	0.01
Rad_type_for_disk(0=bb,1=models,2=uniform)_in_final_spectrum	0
Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table)_in_final_spectrum	3
spectrum_wavemin			3500
spectrum_wavemax			5500
no_observers			1
angle(0=pole)			30
live.or.die(0).or.extract(anything_else)	1
Select_specific_no_of_scatters_in_spectra(y/n)	n
Select_photons_by_position(y/n)	n
spec.type(flambda(1),fnu(2),basic(other)	1
Extra.diagnostics(0=no)		0
Use.standard.care.factors(1=yes)	1
reverb.type(0=None,1=Photon,2=Wind)	2
reverb.disk_type			2
reverb.path_bins			1000
reverb.visualisation		0
reverb.filter_lines			1
reverb.filter_line			49
Photon.sampling.approach(0=T,1=(f1,f2),2=cv,3=yso,4=user_defined)	7

Running the .pf file

The .pf file is run as you would run a standard python input file, and will output a roughly GB-sized .delay_dump file along with the standard python outputs. You need at least twice the expected final file size in order to successfully produce this file- each thread writes to its own file, then they're all stapled together. If you're left with files like .delay_dump2, then you've potentially got a broken run.

Any non-photon (i.e. reverb.type > 0) run cannot be resumed successfully. This is because I've not had time to write the code to write all the arrays to disk, and they're likely to be quite large (reverb.path_bins * Wind.dim.in.x * Wind.dim.in.z * (1+reverb.matom_lines) doubles).

.py file

Running the .py file

Install all the required python libraries (Astropy, Numpy, SQLalchemy). Then just run like a normal python file!