Extinction for chromatic objects

[rmandelb]

I've been thinking about the issue of extinction for chromatic objects. For example, for photo-z tests or whatever you might want to be able to draw an image through various bandpasses and also allow for a different extinction law and extinction zero-point.

I thought of a few options and am curious for opinions on what makes the most sense (pinging @jmeyers314 @barnabytprowe @rmjarvis ):

1. Make a method of the SED class called something like withExtinction that will create a new SED object from another one, with an extinction law applied (as specified via kwargs),

extincted_sed = sed.withExtinction(...)

2. Make some kwargs for the ChromaticObject.drawImage() method that applies the extinction at the time of drawing, without requiring a new SED to be created.

I can see a use for both of these. The first could be useful if you anticipate re-using the object with extinction and don't want to have to do the calculations each time you draw, but the second could be useful if you only plan to use the object once.

Any thoughts? Any interest? :)

[barnabytprowe]

I think that sounds very cool!

[jmeyers314]

Applying extinction to an SED should be as simple as multiplying it by a function of wavelength, I think, which is already supported. I'd be happy to see some methods added to make this more convenient, though.

We could steal this file from astropy/specutils and add it to the GalSim repository (which I believe is allowed by the specutils BSD license), and/or modify it if needed to conform to GalSim conventions.

A more complicated question is whether it makes sense to be able to add extinction to ChromaticObjects inside of drawImage. I think this is okay, but we would need to think through all the possible (or at least meaningful) combinations of sums and convolutions of profiles. Another thing to consider is whether to provide keywords for both Milky Way and host galaxy extinction in drawImage.

[rmandelb]

We could steal this file from astropy/specutils and add it to the GalSim repository (which I believe is
allowed by the specutils BSD license), and/or modify it if needed to conform to GalSim conventions.

Oh, that's handy. I think we should steal and modify, which I believe is permitted.

A more complicated question is whether it makes sense to be able to add extinction to
ChromaticObjects inside of drawImage.

Hmmmm. I originally thought "of course it does" because you could imagine getting an extincted SED and then using it for a bunch of objects. In what case do you think this could be problematic? Are you worried about people adding two objects corresponding to galaxy components, and forgetting to apply extinction to one of them?

If there are any problematic cases, we could simply have SEDs store whatever extinction law has been applied to them, so we can check for those problematic cases and throw an exception or warning when they happen.

Another thing to consider is whether to provide keywords for both Milky Way and host galaxy
extinction in drawImage.

I worry about making drawImage really cumbersome. For a single type of extinction, we presumably will have to require people to pass in E(B-V), and allow them to optionally pass in R_v and the extinction model (assuming we're borrowing the file you linked to), and maybe to use A_v instead of E(B-V). So that's something like 3-4 kwargs. If we want to permit MW and host galaxy extinction, that's 6-8 kwargs, which seems a bit much.

Another option would be as follows:

• drawImage would only allow kwargs to do something basic, like choose E(B-V) for the MW and for the host galaxy assuming some particular extinction model and R_v. That's just 2 extra kwargs.
• the SED.withExtinction() method would allow for much more freedom in choosing these parameters.

[rmjarvis]

Personally, I think it makes more sense to have withExtinction be a method of ChromaticObject rather than SED. Or both I guess, since applying to an SED also makes sense. But in particular, if a galaxy is a sum of bulge + disk, it doesn't (typically) have a single well-defined SED, so you would want to apply the extinction to the whole thing.

As Josh said, this is already possible by multiplying the object by a function:

gal_raw = chromatic_disk + chromatic_bulge
gal_obs = gal_raw * extinction_function
im = gal_obs.drawImage(bandpass, ...)

Having a withExtinction method would mostly just make it easier to use standard A_v or E(B-V) formulations of the extinctions so the user wouldn't have to look them up themselves to get the right function.

gal_obs = gal_raw.withExtinction(Ebv = 0.2, Rv = 'MW')

I think that's clearer than adding extinction parameters to the drawImage method. Especially since you can chain them together, so it's not much extra typing:

im = gal.withExtinction(0.2).drawImage(bandpass)

vs.

im = gal.drawImage(bandpass, Ebv = 0.2)

Another possible UI for this would be to have an extinction function (a bare function that is, not a method of anything) that would return an SED object corresponding to the correct extinction function. Then the syntax for the above would be:

gal_obs = gal_raw * galsim.extinctin(Ebv = 0.2, Rv = 'MW')

And the same syntax would work to apply it to an SED. (At least I think you can currently multiply two SEDs together.)

sed_obs = sed_raw * galsim.extinction(Ebv = 0.2, Rv = 'MW')

[rmandelb]

But in particular, if a galaxy is a sum of bulge + disk, it doesn't (typically) have a single well-defined
SED, so you would want to apply the extinction to the whole thing.

I would think it would be the case that the result should be the same regardless of whether you do (bulge+disk)*extinction or bulge*extinction + disk*extinction. What am I missing?

At least I think you can currently multiply two SEDs together.

You can multiply an SED by a callable function, and an SED is a callable function, so I would think this would work. Hopefully Josh will set me straight if I'm wrong.

Is there any reason not to have a bare extinction function and a withExtinction() method of the SED class? (with the latter calling the former) I realize it provides semi-redundant functionality, but there could be use cases where people will want to get the function themselves and do something with it, and others where they just want to apply it.

[rmjarvis]

I would think it would be the case that the result should be the same regardless of whether you do (bulge+disk)_extinction or bulge_extinction + disk*extinction. What am I missing?

It is. But if you already have the galaxy profile made, you don't want to have to remake it multiple times to view it at different galactic latitudes by applying the extinction separately to each component. Plus, we may someday have chromatic galaxies that are not separable this way, so it might not be possible to effect an extinction only via the SED.

Is there any reason not to have a bare extinction function and a withExtinction() method of the SED class?

That would be a bit unpythonic. The maxim is one and only one obvious way to do anything. Plus if obj.withExtinction(...) is equivalent to obj * extinction(...), there doesn't seem to be any advantage to the method version.

[rmandelb]

The maxim is one and only one obvious way to do anything.

Yes, I realize that, the reason i was asking is that there are some applications where people will want to actually check out / manipulate / plot the extinction they are using, so having an extinction function would facilitate that. Having both would facilitate that at the expense of making multiple ways to draw an object with extinction.

[rmjarvis]

Why not just use the function version then? I don't see the need for the withExtinction method if we have the extinction function. And, as you say, there are other advantages to having the function available.

[rmandelb]

Yes, I guess that makes sense.

I volunteer to do this (relatively small) bit of work sometime in the next ~month or so. If someone is super eager and wants to jump in and do it first, that’s fine with me ;) but I am happy to do it myself if that doesn’t happen!

[jmeyers314]

As Josh said, this is already possible by multiplying the object by a function:

gal_raw = chromatic_disk + chromatic_bulge
gal_obs = gal_raw * extinction_function
im = gal_obs.drawImage(bandpass, ...)

This does probably work right now but it's kind of a side-effect. I never intended for ChromaticObjects to be multipliable by SEDs or other functions wavelength. The reason this exists is due to ChromaticObject.dilate(), which works by expanding and then rescaling the flux. Since dilate() can accept a wavelength-dependent function as its argument, I needed to be able to apply wavelength-dependent flux rescalings. In retrospect, I probably should have done this by directly manipulating the object attributes, rather than allowing __mul__ to accept functions of wavelength. Note that __div__ doesn't work the same way: you can only divide ChromaticObjects by scalars.

I'm a bit torn about applying extinction directly to ChromaticObjects and instead of to their constituent SEDs. This is mostly because for internal extinction, you would want to apply a different curve to each component of a multi-component galaxy. If we want Milky Way extinction to be directly multipliable into ChromaticObjects, we should probably take a closer look at __mul__.

[rmandelb]

Just realized I never responded to this: if this wasn't an intended use (multiplying ChromaticObjects by SEDs) do you think the output of doing so is going to be actually wrong? If so, we should probably fix it so it works, or refactor how dilate works.

[jmeyers314]

Hi @rmandelb,
I think multiplying a ChromaticObject by a function of wavelength probably will accomplish what a user expects, it's just inconsistent with __div__, and the multiplier gets stored in a private attribute with a somewhat cryptic name (_fluxFactor) that should probably be made public and less cryptic if we want this to be a supported part of the API.

One other small point is that I'm not sure we actually want to allow multiplication by SEDs. Atomic ChromaticObjects should already have an SED associated with them, and I can't think of a reason why one would want to multiply together two SEDs. Multiplying by an extinction function would be okay, of course, but that's not an SED. Maybe that's too pedantic to actually enforce though, I'm not sure.

[rmandelb]

Sorry, I actually meant to say "multiplying by an extinction function" (that is a function of wavelength), not multiplying by an SED object per se.

[jmeyers314]

I figured :), but thought it was worth raising the question: should we actively try to catch and raise an exception if you try to multiply a ChromaticObject by an SED?

[rmjarvis]

I think this could be valid either way:

1. Let SEDs model either flux(lambda) or flux_ratio(lambda). Then extinction would also be an "SED" at least as far as the code is concerned and it would be appropriate to multiply them.
2. Have a new class for wavelength-dependent flux scaling, which is what an extinction would be. Then let all ChromaticObjects be multipliable by an Extinction, using the existing _fluxFactor implementation. Also, a GSObject times an Extinction should result in a ChromaticObject just like an SED does.

Either of these options seems fine by me.

[jmeyers314]

I'd vote for option 2 for two related reasons:

1. Dimensionally consistent. Although SED can accept input proportional to any of erg/nm, erg/Hz, or photons/nm, it always keeps an internal representation in photons/nm. I'm not sure there's a good way to convert dimensionless/nm to photons/nm though.
2. I've been thinking of adding an issue to add SED methods for computing AB and Vega magnitudes, and also chromatic quantities like the mean centroid shift due to DCR. I think these only make dimensional sense for option 2 though.

[rmandelb]

I also prefer option 2.

[barnabytprowe]

Also a vote for 2 from me, I think the arguments you put forward are persuasive Josh!