AGN Model SEDs don't work & general issues with "Models" code logic

[jhmatthews]

This is a bug related to #294 but with some issues of larger significance.

The base level problem is that when one tries to run an AGN model with "models" for the spectral type, it doesn't work. I'm running Nick's SEDs documented in issues #249. If you enter "1" in answer to the question "Rad_type_for_agn()_to_make_wind", the code will try and create a model, but it doesn't actually run and fails with

model: Wtot must be greater than 0 or something is badly wrong

which is a result of not having any "weight" in the models routine. I'm not entirely clear why this happens, but I note in passing that the code is passing around "geo.alpha_agn" as the variable to use to work out which model to use, but this is never obtained from the user and so is zero.

The code logic in the general get_models machinery is incredibly confusing to me - we seem to use spectype in very convoluted ways in get_spectype() for example. My main questions/issues are as follows:

• I/we want to be able to read in arbritrary SEDs for the AGN continuum source. These are normally going to be one model, so it seems silly to supply a grid of identical models. However, we also want to avoid code duplication so we want to use as much of the current models machinery as possible. For example, the code should use the same machinery to generate a photon from a model AGN SED as from a disk annulus. How should I go about this in the clearest way possible?

The code logic for spectypes needs improving. I note that string I/O parameters with a map to the relevant #define quantities from python.h would get around a lot of the confusion in get_spectype()

[kslong]

One of the problems I noticed when I looked at this previously was that it appeared T and g were the parameters that Python used to interpolate models and of course this is irrelevant for AGN. I am not sure this was fixed.

In any event agree this is messed up conceptually and we need to think how to deal with this in a general way.

Note that I would have thought that one could read in a single model without causing anything to fail.

[jhmatthews]

Right Knox. I think one can provide a fudgey fix by ensuring that geo.agn_alpha is set to a non-zero value, but I'm not certain here. The user would then have to supply a list of four models for python to read in and do a weird dummy interpolation. However...

Clearly this is not the best if you just want one SED, but it might take a bit of thinking to work out how best to do this. Currently, the models are actually read in by get_models(), which is called by get_spectype() with a hardwired number of params (2) in the call.

We could allow get_spectype to know if the user required a grid of models to be interpolated on or one individual model. Then the code could use get_one_model() to get the individual model component. One would need to design this in a robust enough way so that it allowed to grid of models option to still work in tandem with this approach. It's a little tricky.

I note also that the temperature and log_g that are passed to one_continuum are not actually used, except to check if they have changed from last time.

[jhmatthews]

Again, one could consider a more "object oriented approach", where one specifies a number of radiation sources, and each radiation source has a geometry (disk, point, sphere), spectype (model, BB, power law) and corresponding set of parameters. It depends how much we care about preserving the "defaults" approach python has historically had, and how much work this would be (probably a lot).

[jhmatthews]

@kslong are you able to clarify a little as to what the motivation for the Models code logic is. I am trying to figure out a solution that avoids code duplication and is readable and simple, but I don't know what the philosophy behind the "comps" structure and get_models() function really is.

My feeling is that I need to include an additional "type" of some kind here, so that the code knows whether the model in question is something that needs to be interpolated on, because it has parameters, or whether the model is just one single file. Bearing in mind that my solution should be:

Rad_type_for_agn()_to_make_wind    1
SED_file         my_spectrum.dat 

where my_spectrum.dat not a list of models but one individual model.

[jhmatthews]

Note also, that currently spectype is not set in the case of reading in a model. It defaults to 0, so it manages to get through the code logic in the case of one component that uses models, but in principle the code allows you to have multiple model components, so this needs fixing.

(all the other spectypes are set to negatives via #define statements)

[kslong]

I promise to look at this today or tomorrow, hopefully today.

[jhmatthews]

Right. Multiple issues here. Here is my best attempt to create what is probably non-exhaustive list of problems I've found. This is complicated, so sorry if they aren't as clear as they could be, it's almost impossible to understand and explain all the issues.

Mostly specific to "AGN"

• i) currently, the "correct" dummy variables for interpolating on a dummy grid (two parameters) aren't read in properly or passed to agn_init properly, so when it calls agn_init it goes wrong. This produces the first issue documented, I think.
• ii) emittance_continuum returns an emittance (i.e. a radiative flux through a radiating surface) which needs to be turned into a luminosity by multiplying by the area of said surface. In the AGN case, this isn't done and it's not clear what the emittance . Furthermore, emittance expects to read in an Eddington flux, I think, which is the emittance over 4PI.
• iii) in the ionization cycles, emittance_continuum computes an emittance in a given spectral band. That's fine, but we currently make sure the photon generation CDF is normalised properly by integrating over this band and comparing to the total luminosity. This normalisation isn't done correctly, because we don't know what lum_agn to compare it to, and we also don't know what fraction of the central source luminosity comes out in that band. Similar issues may occur in the spectral cycles.
• iv) AGN models need to be done by interpolating on a pseuo-grid of models, which is messy. Ideally one might be able to read one single file.

General problems

• v) spectype is not currently set to the component in question. In theory, this means one could read in four components (star + disk for each of ionziation and spectral cycles).
• vi) if one sets spectype to be incremented -- so that one could in theory have different model files for different cycles, say -- the code can segfault,
• vii) even if we read in a different spectype for the spectral cycles, we don't actually read in the parameters pertaining to that spectral type. This is OK in the case of, say, stars, because we always have to have a temperature and mass, but in AGN one could change the "BB central source" in ion cycles to a "power law central source" in spec cycles and the code would just use garbage parameters.

I believe ii) and iii) can be fixed by ensuring that the model file for AGN is in luminosity units and having a variable passed to emittance_continuum which specifies the units of the models. i) can be fixed by passing some dummy variables to emittance_continuum which gets the models, but that's a kluge.

[kslong]

There are a lot of issues being co-mingled here, which perhaps is necessary, and we need to figure out how to break them apart. Probably, @jhmatthews you and I at least should chat. I do think we need to get the "use_cases" down that we are trying to solve. That said....

1. I have no objection to creating a new type of radiative model that reads in a single model, but we should make it a possibility for all sources of radiation, not just AGN of course. Presumably, to keep consistence with the other models this should be lambda, flambda. According to the notes in get_models.c that is the units for the Kurucz models, etc, which also conforms to my memory of this. Ideally, we would use the Model, ModSum structures in some way. If we create the option for reading in a single model it should be generic to all radiation sources, though physically, we would not normally use such a model for something like a disk. (Aside, even if we were to do as Nick is doing for his X-ray binary case, we may want to preserve the possibility of re-radiation)

If we have a single model, I would suggest we want to read in (or calculate a luminosity) somehow, to determine how many photons to create for this component of the system. I assume then to deal with bands, one would assume that the model integrated over all wavelength ranges equals the luminosity, and figure out the relative luminosity from the fraction in a particular band.

2. The normalizations, in continuum.c are indeed set up for what both the Kurucz models and Tlusty/Synspec produce, which is the Eddington flux.

This is from my old notes

f Rwd2

• = 4pi ----- H d2

with f the observed flux, H the model flux, Rwd the white dwarf radius and d the distance. Ivan's models are in Eddington flux (H) as are the Kurucz models.

I did a lot of work to verify that Python could did actually reproduce this in the past, for the WD, and I also did a good bit of work to assure that disks produced spectra of the right luminosity, though this could be re-investigated. The notes I found were more about whether BB disks and Kurucz model disks produced the same luminosity, than whether either was correct.

3. It's not obvious to me why reading in a sequence of models for AGN does not work, except that as I mentioned earlier, the current code thinks the two dimensions should be T and g, and we don't read T and g in for AGN. If this were fixed, I believe it should work properly, but the real (If by AGN we really mean central source this may not be necessary)

(Aside, the current program assumes that the first parameter is a T, and extrapolates assuming the spectrum is BB-like when you get out of the T range of the models. If a T-like parameter is not what one wants, then obviously there are problems.)

4. As I recall, and in looking at the code, the idea of comps was that one might have multiple sets of models that one needed to read in, e.g one for the star and another for the disk. The way this is used in practice is that one often reads in Kurucz models for the ionization stage, and finer Hubeny models for the detailed spectrum stage. All of the models are read into the program at the beginning so you need multiple components. The general structure resembles the sorts of things we do for the atomic data, where the various structures are linked together, with comps being the more general structure which looks at the detailed models. I think we should generally continue with this approach.

5. It seems to me we need to distinguish clearly between the situation where we want to have a single spectrum, in which one can do almost anything, and a distributed source where one wants a run of models with radius, and some simple way to say how the model parameters scale with radius.

6. This history of all of this, of course, was that Python was originally designed for CVs, and thermal sources. AGN have been an add-on; we have never gone back and re-engineered the code for this in a systematic way for sources that were not in some sense thermal. So I guess we should not be surprised there are issues. We should just make sure we don't create yet another band-aid.

[Higginbottom]

Looking at this today - I can see that the way we have tried to do this in the past is to send dummy model variables to the AGN part of the code, that then allows us to pretend to interpolate.

There is an XXX comment in agn.c referring to this - and the way forward as I can see it is to actually code up a new models code, where one sends a simple flux vs frequency set of points that the code will use to make a spectrum. It would be used when one doesn't want to use the whole stellar atmospheres part of the code.

[kslong]

So I take it what you all would like is for a single model to be read in (presumably using the code we have to read in each model and not some duplicate code) that would be renormlized by a luminosity, and then to bypass the interpolation. We'd like to store that directly in comp[spectype].mod, and somehow make sure that when one tries to create a photon from the model, we make sure that the routine does not try to interpolate.

This seems reasonable to me, my only concern being that this option be melded as much as possible with our current approach.

[Higginbottom]

Thinking g a bit more....

If we want to keep things as closely allied to the current modes as possible, we should probably keep the wavelength vs flux scheme (even though this will be a bit annoying for X-ray work!!).

We can use the Model structure, the name contains the name of the file we are reading in, the par field might be superfluous, or it could be used to hold the luminosity of this spectrum between some limits, in order to scale. This would make things a bit more flexible.

We also need to think about how we are going to scale this - for the power law spectrum, Bremstrahlung, and broken power law we use the 2-10keV luminosity. I guess we could keep this - but if we wanted this bit of code to be a bit more general, then 2-10keV is not ideal for lower frequency SEDs. I guess that the simplest way of doing things would be to make the luminosity - by default - run over the entire range of the imported model, or alternatively we could have the limits also a user entered number, or we could default to 13.6eV to infinity. So a few choices:

Definitions of luminosity for the single SED file: 2-10keV 13.6eV to infty User defined (with extra inputs) Latched to the end points of the supplied model file

I'm going to start off with 2-10keV - and we can mod later.

I'll also start off by leaving the par field empty at read in - but storing the computed 2-10keV luminosity there...

[Higginbottom]

I'm working with the idea that if one only puts one model in, then that is the model that is used with no interpolation. Just wrangling with the very odd (simplistic) way in which the emittance of a model is calculated. Essentially if there are no points in the model that lie in a band (i.e. I was trying to make a very simple broken power law with just a few points defining bands) then emittance_continuum returns a zero. So one needs a super fine lambda vs flux file. Should really code up a better calculation, but at the moment I'm going to just make a super fine file to use....

[Higginbottom]

I’ve tried to make modifications to python to allow one to properly use a model SED with an AGN type source.

Now, there were some strictures I was working under, primarily that I was to use the current code rather than making new code up. So thats what I did, and this is why there are a couple of fiddles….

The code is in a branch of the main python repo, called agn_model

The way one applies a model SED is as follows:

1: Make yourself an SED file - this must be a two column file, with lambda in angtroms and L_lambda in ergs/s/angstrom. This file must: a: Be sufficiently finely sampled that you get lots of points in any band b: Have at least some power in the 2-10keV band

The reason for a is that the way the code computes the emissivity of a model is rather simplistic - lambda F_lambda is hard wired, and there is no interpolation if bands don’t line up nicely with points in the model file. And we need to have some power in the 2-10keV band so we can use the agn_luminoisty to scale the model in the absence of a temperature.

2: Make a file that ‘points’ to this SED file. It should contain just the name of the file (assuming it will be in the execution directory) and some number - I use -1 just as a marker. This allows the existing model machinery to read in the file and populate the model structure.

3: In the .pf file, just set the source type to model and point the code to the file you made in (2) and away you go. The code should use the model you have provided, and scale it so the 2-10keV luminosity is as requested in lum_agn.

So the relevant bits of the .pf file are:

Parameters for BL or AGN

QSO_BH_radiation(yes,no) yes Rad_type_for_agn(0=bb,1=models,3=power_law,4=cloudy_table,5=bremsstrahlung)_to_make_wind 1 Input_spectra.model_file AGN_spec.ls lum_agn(ergs/s) 1e44

I attach my AGN_spec.ls file (the pointer file) and my AGN_bpow.txt which is a model for a cloudy style broken power law.

AGN_spec.ls txt AGN_bpow.txt

The code changes are very minimal, and essentially just say “If there is only one model file - don’t try and interpolate, just use that one” and a few tests to tell the user that if you want an AGN to use a model, you should only supply one model file.

Please feel free to try and break this - in particular with any favourite runs that use the original models in disks or stars. The examples file runs fine, but I’m quite sure there is a lot of scope for terrible errors.

This is a bit of a kludge - but to do anything more would be a big job and would require lots of new modes/switches/parameters - particularly if one wanted to be able to use nu Lnu data for example..

I’m going to take a little look at the emittance_continuum code because I don’t like the way it is unable to interpolate between points if a band edge falls between two points. It means (as I said above) that you need a very fine model to make things work properly, and I really don’t see why is shouldn’t be easy to make it so one could supply a broken power law model with essentially just a start and end point and a point at each break….

[Higginbottom]

I have further modified emittance_continuum to use a numerical integration to compute the total emittance between bands rather than the rather 'ad hoc' scheme.

[Higginbottom]

From James:

I tried this out with Gordon Richards' mean SED which I interpolated to get it high res. It looks to be working well at least in terms of the shape, although I haven't checked the normalisation yet (note that the blue python curve is higher than the orange original one because it's been scaled up to the luminosity in the parameter file). so thanks!

I did wonder about 2 things:

• is there a reason why we can't have the user just specify the filename of the spectrum in the parameter file? It might take a bit of extra code logic, but better user experience?
• Probably I would have thought it would be better to define the lum_agn in a non 2-10keV way, or be flexible.

[Higginbottom]

I agree with both your points. However as you say both will require significant extra coding because of the way spectral types are input. The code that reads in models has no knowledge of wether it is reading in a model for a star or an agn. So one would need a new mode - say agn model. Then one would need to disable the stellar model for an agn type source or we are back to square one. Same with having a more complex luminosity command. You would either need a special luminosity for a model input, or one would need to recode all the other bits of code that currently assume 2-10 Kev luminosity in the case of an agn. My remit was to minimally change code and reuse as much as possible, so that’s the reason I did it as I did. We should all have a talk about it before making more significant changes. I don’t want to do a load of work only to have it chucked out!

[kslong]

I fully agree. You have done what was asked and we should talk about this before expanding the issue to make it more user friendly/general. We really need a discussion at the same time of a couple of other issues #301 and #307 and the general issues that the AGN models inputs are treated specially (and as a result) it is not possible to specify all spectrum types for all types of models (e.g. cv’s, stars)

[kslong]

@Higginbottom I have compared cv_kur which is in the regression suite in dev and in agn_models. There are clearly detailed differences in the results of the spectra that were produced. Were you expecting the standard models to have been affected at all by the changes you made to allow for a single model. I regard this as a concern that should be understood before merging.

[Higginbottom]

There could well be slight differences now I think about it. One of the other changes I made to the code was to use qromb to integrate a model over a band. This was because the original scheme was not very good if you had a coarse grid of the model - I was getting glitches in my output spectrum which I thought was due to band edges interfering with points in the model. These should be very slight differences I'd have thought/hoped??

[kslong]

OK, we may be able to check this. It looks like you isolated this update in a specific commit. So “we” should compare a cv_kur.pf before this commit, to see if it agrees with dev, and after it to see whether things have changed.

From: Higginbottom notifications@github.com Reply-To: agnwinds/python reply@reply.github.com Date: Tuesday, November 27, 2018 at 8:51 AM To: agnwinds/python python@noreply.github.com Cc: Knox Long long@stsci.edu, Mention mention@noreply.github.com Subject: Re: [agnwinds/python] AGN Model SEDs don't work & general issues with "Models" code logic (#301)

There could well be slight differences now I think about it. One of the other changes I made to the code was to use qromb to integrate a model over a band. This was because the original scheme was not very good if you had a coarse grid of the model - I was getting glitches in my output spectrum which I thought was due to band edges interfering with points in the model. These should be very slight differences I'd have thought/hoped??

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/agnwinds/python/issues/301#issuecomment-442085703, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ACaeVVjNXtBsK9HV0KCPddy8171_8Nx0ks5uzVFrgaJpZM4P1Shc.

[Higginbottom]

Leave it with me!

[kslong]

Thanks, Nick. I think the only thing else need is a comparison of the spectra from cv_kur before and after the change. To be sure we don’t have a problem I would increase the number of photons by an order o magnitude, leaving the number of Cycles the same. That will eliminated questions of statistics. Obviously it is the disk and star spectra that are created that matter here.

Knox

From: Higginbottom notifications@github.com Reply-To: agnwinds/python reply@reply.github.com Date: Wednesday, November 28, 2018 at 9:07 AM To: agnwinds/python python@noreply.github.com Cc: Knox Long long@stsci.edu, Mention mention@noreply.github.com Subject: Re: [agnwinds/python] AGN Model SEDs don't work & general issues with "Models" code logic (#301)

I have run the model with python commits

98b23a0https://github.com/agnwinds/python/commit/98b23a0ccf02e8df859d8c3692a924cbd0eef253 - current dev branch f01d779https://github.com/agnwinds/python/commit/f01d7793546760717af4f6a1c73b39e30432a6b3 - agn_model branch before the qromb change

the cv_kur.spec files are numerically identical.

As noted by KSL, a run with commit

f01d779https://github.com/agnwinds/python/commit/f01d7793546760717af4f6a1c73b39e30432a6b3 - the latest commit in the agn_model branch produces different spectra.

The only question is - can we tolerate the changes?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/agnwinds/python/issues/301#issuecomment-442479597, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ACaeVcYOay9cVAvKJKCZfKUB5VcHvfsYks5uzqbBgaJpZM4P1Shc.

[Higginbottom]

UPDATED corrected.. I have run the model with python commits

98b23a0ccf02e8df859d8c3692a924cbd0eef253 - current dev branch 4b4d6985dbbc68200a210c58aa86c289b0fd8615 - agn_model branch before the qromb change

the cv_kur.spec files are numerically identical.

As noted by KSL, a run with commit

6956df4ccc3bb511f883fa6e4e591e9a1114e5b6 - the commit in the agn_model directly after qromb was used branch produces different spectra.

The only question is - can we tolerate the changes?

[Higginbottom]

OK, here is a comparison plot with 10x the number of photons, for 62 degrees..

It is very similar! Here is the same thing - slightly smoothed..

[kslong]

That looks good enough to me. I’d be happy for you to go ahead and merge.

Knox

From: Higginbottom notifications@github.com Reply-To: agnwinds/python reply@reply.github.com Date: Wednesday, November 28, 2018 at 10:56 AM To: agnwinds/python python@noreply.github.com Cc: Knox Long long@stsci.edu, Mention mention@noreply.github.com Subject: Re: [agnwinds/python] AGN Model SEDs don't work & general issues with "Models" code logic (#301)

OK, here is a comparison plot with 10x the number of photons, for 62 degrees..

[a62p0 50spectrum_summary]https://user-images.githubusercontent.com/3329213/49167948-65fe2600-f32e-11e8-9b9f-79df7f8b76ba.png

It is very similar! Here is the same thing - slightly smoothed..

[a62p0 50spectrum_summary]https://user-images.githubusercontent.com/3329213/49168004-862de500-f32e-11e8-85b5-13ce508d6aca.png

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://github.com/agnwinds/python/issues/301#issuecomment-442521810, or mute the threadhttps://github.com/notifications/unsubscribe-auth/ACaeVXt98vEGD-rNVaYamlVISaQLmnVBks5uzsBYgaJpZM4P1Shc.

[jhmatthews]

Likewise - just to check, is the -0.5 just an offset? and what are the units of the differences?

[Higginbottom]

Yeah, the -0.5 is an offset so you can have multiple files to compare.

The calculation shown as the difference is actually the ratio -0.5. It is really intended to show where there are differences rather than quantify them…

I can make more plots if necessary..

N

On 28 Nov 2018, at 17:01, James Matthews notifications@github.com wrote:

Likewise - just to check, is the -0.5 just an offset? and what are the units of the differences?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/agnwinds/python/issues/301#issuecomment-442523772, or mute the thread https://github.com/notifications/unsubscribe-auth/ADLMvWmwTUXXFI2-3_n7OYdI-TV-LAzhks5uzsGBgaJpZM4P1Shc.

[jhmatthews]

I would say go ahead!

[Higginbottom]

I think this issue is all fixed up now. Closing