Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix some issues with SkyCalcTERCurve #245

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Conversation

hugobuddel
Copy link
Collaborator

This PR makes rescaling of a SkycalcTERCurve possible, but I'm not sure that is a reasonable thing to do. However, from a technical perspective, it should be possible to rescale a SkyCalcTERCurve, but there were two bugs preventing this. So this PR has two components:

  • @teutoburg (and perhaps @oczoske ) could you review my solution to the two technical problems with SkycalcTERCurve?
  • @oczoske (and perhaps @teutoburg ) could you please comment on whether it would make physical sense to rescale a SkycalcTERCurve?

Intro

Adding a rescale_emission argument to a SkycalcTERCurve is now technically possible, see the tests:

effects :
-   name : skycalc_atmosphere
    description : atmospheric spectra pulled from the skycalc server
    class : SkycalcTERCurve
    include : True
    kwargs :
        ...
        rescale_emission:
            filter_name: "!ATMO.background.filter_name"
            filename_format: "!INST.filter_file_format"
            value: "!ATMO.background.value"
            unit: "!ATMO.background.unit"

Bugs fixed

There where two problems, both related to the fact that synphot (by default) refuses to rescale a spectrum if the overlap between the spectrum and the filter it rescales to is not correct. That is, the spectrum returned by Skycalc, should overlap the desired filter, but is just too short:

  • The table returned by Skycalc does not actually contain the specified end point, e.g. in this case, 2.5 um. This is resolved by duplicating the last row of the table, but setting the wavelength to the specified endpoint.
  • The input to Skycalc needs to be in nm, but converting 2.5 um to nm makes it just short enough for synphot to complain: (2.5 * u.um).to(u.nm).value = 2499.9999999999995. This is fixed by 'improving' the unit conversion.

Now irrespective of any rescaling, it is my opinion that fixing these two problems is worthwhile anyway.

@teutoburg, could you turn on your WTFs/minute-meter and determine how stupid these solutions are?

Physics

Our documentation seems to imply that it is possible to rescale the background. In particular, the getting started notebook explicitly states:

For example, setting the atmospheric background level is achieved thusly:

opt.cmds["!ATMO.background.filter_name"] = "K"
opt.cmds["!ATMO.background.value"] = 13.6

I think this makes some sense for the AtmosphericTERCurve, which is what that notebook uses. However, then users go and proceed to rescale the atmosphere when using MICADO or METIS, but the SkycalcTERCurve that is used by default in Armazones.yaml does not allow the background to be scaled.

At least two people have complained about this, some with more detail than others, which led me to investigate this, and fix the bugs that prevented the background to be scaled.

However, I'm not sure it makes sense at all to rescale a SkycalcTERCurve background in this way:

  • Supposedly, the arguments given to Skycalc should capture everything so it should not be necessary to scale the background.
  • If scaling is required, then there is (probably?) a physical way to do it through Skycalc, e.g. by increasing the temperature.

And honestly, I don't understand what the rescaling does anyway. The unit is mag, but the background is not a point source. So it is not clear to me what it physically means to "rescale the background to magnitude 17 in the Ks band".

@oczoske could you perhaps give your thoughts on the background rescaling?

TODO later:

Whatever we decide (whether rescaling of an AtmosphericTERCurve should be allowed or not), we should update the documentation, because now it is confusing to multiple users.

Why is this even necessary?? The class should either have a wunit
or not, but not like have it half of the time in meta
@codecov
Copy link

codecov bot commented Jul 5, 2023

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.09 🎉

Comparison is base (8b5e08b) 75.16% compared to head (81bb307) 75.26%.

Additional details and impacted files
@@              Coverage Diff               @@
##           dev_master     #245      +/-   ##
==============================================
+ Coverage       75.16%   75.26%   +0.09%     
==============================================
  Files             152      152              
  Lines           15637    15676      +39     
==============================================
+ Hits            11754    11798      +44     
+ Misses           3883     3878       -5     
Impacted Files Coverage Δ
scopesim/effects/ter_curves.py 84.59% <100.00%> (+1.43%) ⬆️
...opesim/tests/tests_effects/test_SkycalcTERCurve.py 93.44% <100.00%> (+3.44%) ⬆️
scopesim/utils.py 68.16% <100.00%> (+0.61%) ⬆️

... and 1 file with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

Copy link
Contributor

@teutoburg teutoburg left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Other than my comment about the test raising an error:
My only "WTF/min" was in the table (row) indexing when the last row is duplicated. But that's probably because I haven't worked with astropy.table.Table a lot before, and not because there's anything wrong in your implementation 😅. And the new test(s) should make sure it's correct anyway.

sto1 = str(to1)
to2 = 1 / unit_to.to(unit_from)
sto2 = str(to2)
if len(sto2) < len(sto1):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if len(sto2) < len(sto1):
if len(str(to2)) < len(str(to1)):

Save two lines and two temp variables? It's not like they're used again, nor are the expressions long...


2499.99 nm is less than 2.5 um, so there is no overlap anymore between
the source spectrum and the to-be-scaled to spectrum, leading to an
error.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does it still "lead to an error"? Isn't that what was fixed here? Same for the comment at the bottom of this method.

@teutoburg
Copy link
Contributor

To play Devil's advocate (or advocaat 😋): Alternative solution to the unit conversion issue:
Limit the precision in the table, i.e. perform some kind of rounding after the unit conversion, which would also get rid of the .99999999999999999999999995
Not saying this would be better, but I guess we should have a good reason why this would not be better, to show that the proposed solution is better. If that makes sense 😂



def save_unit_to(unit_from: u.Unit, unit_to: u.Unit):
"""Like unit_from.to(unit_to), but u.um to u.nm == 1000, not 999.999.."""
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you state the purpose of the function more explicitely, please? "Like this or that" is a bit vague. Also, the docstring suggests to me that the function should be called safe_unit_to, or what is it that is being saved?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not very fond of that function sav/fe_unit_to - I feel it shouldn't be necessary (why does Unit.to do that anyway?). But I did a similar thing when I put in utils.seq as an improvement over numpy.arange, so I shouldn't say anything other than "use cautiously".

@oczoske
Copy link
Collaborator

oczoske commented Jul 5, 2023

I don't actually know much about skycalc and sky emission in general. In principle, skycalc should provide a "complete" physical model of the sky emission, so manual rescaling makes little sense if you know how to handle the physical parameters. It was introduced with AtmosphericTERCurve as that effect used a very limited number of emission tables provided with SimCADO, and rescaling those to a given sky brightness for some site was quite useful. It might work as a shortcut to specifying involved combinations of physical parameters, but might also easily lead to disastrous results. It might be a good idea to hear Wolfgang Kausch on this question.
The magnitude specified in the rescaling is for a patch of 1 arcsec^2, as is conventional for surface brightnesses. As far as I know, it is not possible to make that explicit in synphot (mag/arcsec^2 is not a proper unit because mag is not a proper unit, but while mags on their own are just about manageable, mag/arcsec^2 are not).

@hugobuddel
Copy link
Collaborator Author

Thanks @teutoburg and @oczoske . I'll think about this a bit more. I'll probably merge it after applying your suggestions, but maybe I'll change my mind, it is a bit ugly this.

@hugobuddel hugobuddel marked this pull request as draft August 23, 2023 12:44
Base automatically changed from dev_master to main April 18, 2024 09:55
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
Status: 📋 Backlog
Development

Successfully merging this pull request may close these issues.

3 participants