Dust

[julienguy]

I have committed an implementation of the Fitzpatrick reddening law.
I added dust.ebv_sfd_scaling = 0.86 and a function dust_transmission(wave,ebv_sfd).
I compare extinction_total_to_selective_ratio and synthetic mags of an extincted dummy spectrum using filters from speclite in a main function (that I can remove if some of you don't like it).

I do not obtain a very good agreement so it would be great if an expert has a look at it.

python py/desiutil/dust.py

R_G(N:BASS+MZLS)= 3.214 =? 0.86x3.729 = 3.207 , delta = 0.0068
R_R(N:BASS+MZLS)= 2.165 =? 0.86x2.471 = 2.125 , delta = 0.0400
R_Z(N:BASS+MZLS)= 1.211 =? 0.86x1.361 = 1.170 , delta = 0.0410
R_G(S:DECALS)= 3.283 =? 0.86x3.678 = 3.163 , delta = 0.1201
R_R(S:DECALS)= 2.200 =? 0.86x2.459 = 2.115 , delta = 0.0847
R_Z(S:DECALS)= 1.215 =? 0.86x1.376 = 1.184 , delta = 0.0313

[coveralls]

Coverage decreased (-1.3%) to 66.706% when pulling 8b6af6d on dust into 9dd77fe on master.

[schlafly]

Trying to track down what values I think should be right... this is what tractor uses:
https://github.com/dstndstn/tractor/blob/main/tractor/sfd.py#L23
which look like your BASS/MzLS values, though they're actually intended for use with DECam. They also match what we quote for the Legacy Surveys:
https://www.legacysurvey.org/dr9/catalogs/#galactic-extinction-coefficients
They are slightly different from what I have at the imaging standard bandpass page
https://desi.lbl.gov/trac/wiki/ImagingStandardBandpass
which probably reflects a very slightly different choice of filter (1 part in 1000).

I have not been able to track down the provenance of the numbers you are currently using for DECaLS.

So, focusing on the numbers I do understand (3.214, 2.165, 1.211), and comparing those with the numbers you compute for DECam (3.678, 2.459, 1.376), I get that these are different by a fixed scale factor of 0.88 (c.f. SF11 rescaling of SFD by 0.86). That's looks close, but it's better than that---the difference between 0.88 and 0.86 is because 0.86 is calculated for a galaxy source spectrum more appropriate to the original SFD calibration, rather than for a 7000 K F star. So I think what you're doing is basically right.

I recommend applying the SF11 calibration in the following way:
A(lambda, F99) = A(lambda, F99, rv)/A(1 micron, F99, rv) * SFD_EBV * 1.029. That's a definition that only uses monochromatic extinctions, so a lot of ambiguity in what extinction means goes away. That doesn't look like the 0.86 rescaling that we quote in abstract, but it's convenient because it uses only monochromatic extinctions

As part of preparing this comment, I tried to trace the history of the number 1.029 again. Here's more information than you want. You should probably just stop reading here. You have been warned. 1.029 is 0.78*1.319, where 0.78 is from the SF11 appendix and 1.319 is a magic number buried in dust_intfilter.pro, a bit of IDL used to compute the coefficients for SFD98. 1.319 is supposed to represent A(1 micron)/E(B-V) for the O'Donnell extinction curve. I can't reproduce that number myself. But it's not far from ext_odonnell(10000)/(ext_odonnell(4404)-ext_odonnell(5428)) = 1.326, where 4404 and 5428 are the effective wavelengths cited for the Landolt B and V filters in the SFD appendix. And if I just treat it as a magic constant I can reproduce the SFD appendix numbers. Presumably SFD was using some other definition of E(B-V), possibly non-monochromatic, or just using different wavelengths. In other words, in the parameterization of SF11, the recalibration is a factor of 0.78, but much of that is because O'Donnell and F99 extrapolate from B-V to 1 micron differently.

[weaverba137]

@julienguy, any chance of adding unit tests to cover these changes?

In general, what is the relationship between this PR and #166?

[julienguy]

• The values of R_g,R_r,R_z in the funcition extinction_total_to_selective_ratio have been modified to match the values of MW_TRANSMISSION_X in the targeting catalogs (DR8,DR9) . There was an error in the previous version.
• Use suggestions of Eddie for the nomalization of the extinction curve.
  Now the agreement between the R_g,R_r,R_z values and values obtained with the extinction curve and synthetic magnitudes are pretty good for DECALS. The residual difference is due to the choice of the fiducial spectrum to compute the mags. For BASS+MZLS, the differences are larger. This is because I have made the choice to keep the values of R_g,R_r,R_z considered in the imaging catalogs which are the same as for DECALS despite the known differences in the filter transmissions.

R_G(N:BASS+MZLS)= 3.214 =? 3.285, delta = -0.0714
R_R(N:BASS+MZLS)= 2.165 =? 2.177, delta = -0.0116
R_Z(N:BASS+MZLS)= 1.211 =? 1.198, delta = 0.0127
R_G(S:DECALS)= 3.214 =? 3.240, delta = -0.0259
R_R(S:DECALS)= 2.165 =? 2.167, delta = -0.0015
R_Z(S:DECALS)= 1.211 =? 1.212, delta = -0.0015

(note the differences have to be multiplied by E(B-V) which is typically <0.05 )

[julienguy]

I have also removed the scale factor dust.ebv_sfd_scaling which is not needed.

[julienguy]

The documentation generation is failing. Except for that it is fine.
This PR addresses issue #165 and should replace PR #166.

[schlafly]

FWIW, in the event we wanted to diagnose the DECam g band difference more, my guess is that this is due to my using an airmass of 1.4 here:
https://desi.lbl.gov/trac/wiki/ImagingStandardBandpass
while speclite uses an airmass of 1.3.

[weaverba137]

@julienguy are you actively working on the documentation problem, or are you asking someone else to do that?

[julienguy]

@julienguy are you actively working on the documentation problem, or are you asking someone else to do that?

I would appreciate if you have a look at it.

[weaverba137]

OK, documentation fixed.

[weaverba137]

@julienguy, it looks like you got tests to pass by commenting out the tests. I don't think that's a very good practice. Are you planning to restore the tests?

[julienguy]

Have you actually seen why I commented out the tests?

[weaverba137]

@julienguy, no, in fact there are no additional comments in test_dust.py that explain why the tests were commented out. I also don't see any discussion of that in this thread.

[julienguy]

After discussion with Stephen, the default is now the most accurate estimate of the total to selective ratios.
In order to retrieve the same result as in the Legacy Surveys catalogs, one would have to use the option match_legacy_surveys=True in the calls to extinction_total_to_selective_ratio or mwdust_transmission.
With the new default values, retrieved from https://desi.lbl.gov/trac/wiki/ImagingStandardBandpass,
the comparison with synthetic magnitudes is better. The residual difference in g is due to a different the choice of airmass
in the effective pass bands.

R_G(N:BASS+MZLS)= 3.258 =? 3.285, delta = -0.0274
R_R(N:BASS+MZLS)= 2.176 =? 2.177, delta = -0.0006
R_Z(N:BASS+MZLS)= 1.199 =? 1.198, delta = 0.0007
R_G(S:DECALS)= 3.212 =? 3.240, delta = -0.0279
R_R(S:DECALS)= 2.164 =? 2.167, delta = -0.0025
R_Z(S:DECALS)= 1.211 =? 1.212, delta = -0.0015

[julienguy]

New doc failures. I let you @weaverba137 fix this if you can.

[sbailey]

I have not evaluated the correctness of this PR, but some structural comments:

• let's directly support Gaia here instead of treating it as a special case on the desispec side
• While we're at it, do we already have coefficients for WISE and other DECam (or even SDSS and PanSTARRS?) bands that we could plug in here to make this more generically useful even if we're not using them yet? I don't want to add feature creep, but it could head of future proliferation of dust-like code (like you are cleaning up in Dust desispec#1246).
• since there are some subtleties here about the default behavior, let's be more proactive with the tests and check:
  • results match Legacy Surveys if photsys='S' regardless of whether match_legacy_surveys=True or False. Test with a few hardcoded (but documented provenance) examples taken from DR9 catalogs spanning small, medium, large E(B-V),
  • for photsys='N', results match Legacy Surveys if match_legacy_surveys=True but does not match LS if match_legacy_surveys=False
  • calling with PHOTSYS='G' (Gaia) works (or at least doesn't crash)

[sbailey]

Let's also support PHOTSYS='G' (Gaia) since that appears in the targeting catalogs as well, with allowable band values of "G", "BP", and "RP". This would move desispec.scripts.stdstars.unextinct_gaia_mags logic over here instead of having Gaia treated as a special case on the desispec side (which I think it still is there in desihub/desispec#1246).

For consideration: older targets/fiberassign files also have PHOTSYS='' meaning Gaia; we should track down whether that is recent enough for us to care (Dec 2020 or later), but I haven't done that yet.

[sbailey]

@segasai could you help with moving the Gaia extinction code out of desispec and into this PR branch so that desiutil.dust is our one-stop-shop for extinction corrections?

[segasai]

I'll see what I can do (on a timescale of few days)

[sbailey]

Is this main leftover from debugging, or is it still useful? I'm ok with sometimes hiding super-specialized scripts under if __name__ == "__main__": blocks under otherwise module code, but if you do want to retain this please add some brief comments about what it does and add a basic argparse so that python dust.py --help would be informative without actually running anything.