mfouesneau / pyphot

suite to deal with passband photometry
https://mfouesneau.github.io/pyphot/
MIT License
57 stars 18 forks source link

Not using unit aware classes, loading the included dat files with Filter.from_ascii can result in different zero points than the HDF5 file #46

Closed christophra closed 1 year ago

christophra commented 2 years ago

I'm sorry the title is a bit vague, I was trying to cover two possibilities because I'm not sure about the expected behaviour. The docs say

the transmission is unitless. Its scaling factor will not affect the flux/magnitude calculations as these are implicitly normalizing the passband (see the equations) but will affect the number of photon photons/s/cm2

So if I understand that right, the issue is either

I understand this doesn't have high priority because _zero_flux and _zero_mag are the ones used for flux/magnitude calculations. But maybe it's still helpful to report.

Steps to reproduce:

  1. Load all the filters in libs/ascii_sources with Ascii_Library
  2. Use get_library to load libs/new_filters.hdf5
  3. Compare two ZTF filters from either.

Expected behaviour: Their zero points (etc.) are identical within rounding error.

Observed behaviour: AB_zero_Jy, ST_zero_Jy, Vega_zero_Jy are reduced by a factor 100

Minimal example

import pyphot
import os
import numpy as np

# Using provided library for comparison
lib = pyphot.get_library()
# Reloading and storing in dictionary
lib_reload = {}
pyphot_dir = os.path.dirname(pyphot.__file__)
pyphot_filter_dir = os.path.join(pyphot_dir, 'libs/ascii_sources/')
for _f in pyphot.phot.Ascii_Library(pyphot_filter_dir).load_all_filters():
    lib_reload[_f.name] = _f

# e.g. ZTF filter dat files are in percent
for filter_name in[ 'ZTF_g', 'ZTF_r', 'ZTF_i']:
    print(f'Filter: {filter_name}')

    filters = [_lib[filter_name] for _lib in [lib, lib_reload]]

    zeropoint_keys = []
    for _q in ['AB', 'ST', 'Vega']:
        for _u in ['mag', 'Jy', 'flux']:
            _k = f'{_q}_zero_{_u}'
            zeropoint_keys.append(_k)
    zeropoint_keys.append('Vega_zero_photons')

    for _k in zeropoint_keys:
        # get the zero point from both the original and reloaded filter
        quantities = [_f.__getattribute__(_k) for _f in filters]

        # accommodate both those with units and the unitless i.e. magnitudes
        # to help with printing
        if isinstance(quantities[0], pyphot.unit.Quantity):
            values = tuple([_q.magnitude for _q in quantities])
            unit = quantities[0].unit
        else:
            values = tuple(quantities)
            unit = ''

        # NumPy uses tolerances approximately to the effect of 
        # "numbers are close both in relative and absolute difference"
        match = np.isclose(*values, rtol=1e-10, atol=1e-10)

        print(f'{_k} is',
              [f'different: \n{values} {unit}', 'identical'][match],
             )
christophra commented 2 years ago

I now found out about pyphot.sandbox and its parallel Unit classes. The bug seems to be fixed in UnitFilter.from_ascii, i.e. when using pyphot.sandbox.UnitAscii_Library in the example above.

mfouesneau commented 2 years ago

Yes, I need to deprecate the old classes.