diff --git a/CHANGES.rst b/CHANGES.rst index 919f9b70aa..6c028d132d 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -6,6 +6,9 @@ New Features - Profile viewers now support plotting with profiles "as steps". [#1595] +- Line flux in the Line Analysis plugin are reported in W/m2 if Spectral Flux is given + in Jy [#1564] + Cubeviz ^^^^^^^ diff --git a/docs/specviz/plugins.rst b/docs/specviz/plugins.rst index aba060a158..7431456cc8 100644 --- a/docs/specviz/plugins.rst +++ b/docs/specviz/plugins.rst @@ -234,6 +234,8 @@ secondary region can be created and selected as the region to fit the linear con The statistics returned include the line centroid, gaussian sigma width, gaussian FWHM, total flux, and equivalent width. +The line flux results are automatically converted to Watts/meter^2, when appropriate. + Redshift from Centroid ---------------------- diff --git a/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py b/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py index 0414e88278..79872b0f8f 100644 --- a/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py +++ b/jdaviz/configs/specviz/plugins/line_analysis/line_analysis.py @@ -7,7 +7,7 @@ from glue_jupyter.common.toolbar_vuetify import read_icon from traitlets import Bool, List, Float, Unicode, observe from astropy import units as u -from specutils import analysis +from specutils import analysis, Spectrum1D from specutils.manipulation import extract_region from jdaviz.core.custom_traitlets import FloatHandleEmpty @@ -332,7 +332,56 @@ def _uncertainty(result): for function in FUNCTIONS: # TODO: update specutils to allow ALL analysis to take regions and continuum so we # don't need these if statements - if function == "Equivalent Width": + if function == "Line Flux": + flux_unit = spec_subtracted.flux.unit + # If the flux unit is equivalent to Jy, or Jy per spaxel for Cubeviz, + # enforce integration in frequency space + if (flux_unit.is_equivalent(u.Jy) or + flux_unit.is_equivalent(u.Jy/u.sr)): + # Perform integration in frequency space + freq_spec = Spectrum1D( + spectral_axis=spec_subtracted.spectral_axis.to(u.Hz, + equivalencies=u.spectral()), + flux=spec_subtracted.flux, + uncertainty=spec_subtracted.uncertainty) + + raw_result = analysis.line_flux(freq_spec) + # When flux is equivalent to Jy, lineflux result should be shown in W/m2 + if flux_unit.is_equivalent(u.Jy/u.sr): + final_unit = u.Unit('W/m2/sr') + else: + final_unit = u.Unit('W/m2') + + temp_result = raw_result.to(final_unit) + if getattr(raw_result, 'uncertainty', None) is not None: + temp_result.uncertainty = raw_result.uncertainty.to(final_unit) + + # If the flux unit is instead equivalent to power density + # (Jy, but defined in wavelength), enforce integration in wavelength space + elif (flux_unit.is_equivalent(u.Unit('W/m2/m')) or + flux_unit.is_equivalent(u.Unit('W/m2/m/sr'))): + # Perform integration in wavelength space using MKS unit (meters) + wave_spec = Spectrum1D( + spectral_axis=spec_subtracted.spectral_axis.to(u.m, + equivalencies=u.spectral()), + flux=spec_subtracted.flux, + uncertainty=spec_subtracted.uncertainty) + raw_result = analysis.line_flux(wave_spec) + # When flux is equivalent to Jy, lineflux result should be shown in W/m2 + if flux_unit.is_equivalent(u.Unit('W/m2/m'/u.sr)): + final_unit = u.Unit('W/m2/sr') + else: + final_unit = u.Unit('W/m2') + + temp_result = raw_result.to(final_unit) + if getattr(raw_result, 'uncertainty', None) is not None: + temp_result.uncertainty = raw_result.uncertainty.to(final_unit) + + # Otherwise, just rely on the default specutils line_flux result + else: + temp_result = analysis.line_flux(spec_subtracted) + + elif function == "Equivalent Width": if np.any(continuum <= 0): temp_results.append({'function': function, 'result': '', diff --git a/jdaviz/configs/specviz/plugins/line_analysis/tests/test_line_analysis.py b/jdaviz/configs/specviz/plugins/line_analysis/tests/test_line_analysis.py index 1caa5fd92f..c659c44e8d 100644 --- a/jdaviz/configs/specviz/plugins/line_analysis/tests/test_line_analysis.py +++ b/jdaviz/configs/specviz/plugins/line_analysis/tests/test_line_analysis.py @@ -136,7 +136,7 @@ def test_continuum_surrounding_spectral_subset(specviz_helper, spectrum1d): plugin.width = 3 # Values have not yet been validated - np.testing.assert_allclose(float(plugin.results[0]['result']), 350.89288537581467, atol=1e-5) + np.testing.assert_allclose(float(plugin.results[0]['result']), 2.153181e-13, atol=1e-15) def test_continuum_spectral_same_value(specviz_helper, spectrum1d): @@ -215,7 +215,7 @@ def test_continuum_subset_spectral_entire(specviz_helper, spectrum1d): plugin.width = 3 # Values have not yet been validated - np.testing.assert_allclose(float(plugin.results[0]['result']), -467.6854635447396, atol=1e-5) + np.testing.assert_allclose(float(plugin.results[0]['result']), -2.79572e-13, atol=1e-15) def test_continuum_subset_spectral_subset2(specviz_helper, spectrum1d): @@ -248,7 +248,7 @@ def test_continuum_subset_spectral_subset2(specviz_helper, spectrum1d): plugin.width = 3 # Values have not yet been validated - np.testing.assert_allclose(float(plugin.results[0]['result']), 32.520521, atol=1e-5) + np.testing.assert_allclose(float(plugin.results[0]['result']), 1.482418e-14, atol=1e-16) def test_continuum_surrounding_no_right(specviz_helper, spectrum1d): @@ -276,7 +276,7 @@ def test_continuum_surrounding_no_right(specviz_helper, spectrum1d): plugin.width = 3 # Values have not yet been validated - np.testing.assert_allclose(float(plugin.results[0]['result']), 39.76685499263615, atol=1e-5) + np.testing.assert_allclose(float(plugin.results[0]['result']), 4.204513e-14, atol=1e-16) def test_continuum_surrounding_no_left(specviz_helper, spectrum1d): @@ -304,7 +304,7 @@ def test_continuum_surrounding_no_left(specviz_helper, spectrum1d): plugin.width = 3 # Values have not yet been validated - np.testing.assert_allclose(float(plugin.results[0]['result']), 146.67186446784513, atol=1e-5) + np.testing.assert_allclose(float(plugin.results[0]['result']), 7.570859e-14, atol=1e-16) def test_subset_changed(specviz_helper, spectrum1d): @@ -335,4 +335,4 @@ def test_subset_changed(specviz_helper, spectrum1d): specviz_helper.app.state.drawer = True # Values have not yet been validated - np.testing.assert_allclose(float(plugin.results[0]['result']), 350.89288537581467, atol=1e-5) + np.testing.assert_allclose(float(plugin.results[0]['result']), 2.153181e-13, atol=1e-15) diff --git a/jdaviz/configs/specviz/plugins/line_analysis/tests/test_lineflux.py b/jdaviz/configs/specviz/plugins/line_analysis/tests/test_lineflux.py new file mode 100644 index 0000000000..44256ecb6c --- /dev/null +++ b/jdaviz/configs/specviz/plugins/line_analysis/tests/test_lineflux.py @@ -0,0 +1,123 @@ +from astropy.tests.helper import assert_quantity_allclose +import astropy.units as u +from glue.viewers.profile.state import FUNCTIONS as COLLAPSE_FUNCTIONS +import numpy as np +import pytest +from specutils import Spectrum1D + +from jdaviz import Cubeviz + +import warnings +warnings.filterwarnings('ignore') + +# Maps input spectrum flux unit to expected line analysis flux unit +expected_lineflux_results = { + u.Jy/u.sr: u.Unit('W/(m2*sr)'), + u.Jy: u.Unit('W/m2'), + u.Unit('W/m2/m')/u.sr: u.Unit('W/(m2*sr)'), + u.Unit('W/m2/m'): u.Unit('W/m2') +} + +test_cases = list(expected_lineflux_results.keys()) + + +# Test cases for unit gaussian tests +def _gauss_with_unity_area(x, mean, sigma): + ''' Gaussian function with area = 1 unit ''' + dx = x - mean + n = 1/(sigma*np.sqrt(2*np.pi)) + g = n * np.exp(-dx**2/(2*sigma**2)) + return g + + +mn = 1.5 +sig = 0.01 + +# 5 cases: replace fl_wave with fnu_freq, fnu_wave, flam_wave, flam_freq, or fl_wave +# to try the different cases here +unit_flux_gaussian_test_cases = [] + +# unit-flux gaussian in frequency space +freq = np.arange(1, 2, 0.001)*u.Hz +flux_freq = _gauss_with_unity_area(freq.value, mn, sig)*1.0E26*u.Jy +fnu_freq = Spectrum1D(spectral_axis=freq, flux=flux_freq) +unit_flux_gaussian_test_cases.append(fnu_freq) +fnu_wave = Spectrum1D(spectral_axis=fnu_freq.wavelength, flux=flux_freq) +unit_flux_gaussian_test_cases.append(fnu_wave) + +# unit-flux gaussian in wavelength space +lam = np.arange(1, 2, 0.001)*u.m +flux_wave = _gauss_with_unity_area(lam.value, mn, sig)*1.0*u.W/u.m**2/u.m +flam_wave = Spectrum1D(spectral_axis=lam, flux=flux_wave) +unit_flux_gaussian_test_cases.append(flam_wave) +flam_freq = Spectrum1D(spectral_axis=flam_wave.frequency, flux=flux_wave) +unit_flux_gaussian_test_cases.append(flam_freq) + + +def _calculate_line_flux(viz_helper): + ''' + Returns the line flux calculation by opening the Line Analysis Plugin + Assumes the plugin hasn't been opened yet + ''' + # Open the plugin and force the calculation + viz_helper.app.state.drawer = True + line_analysis_plugin = viz_helper.app.get_tray_item_from_name('specviz-line-analysis') + line_analysis_plugin.open_in_tray() + # Retrieve Results + for result in line_analysis_plugin.results: + if result['function'] == 'Line Flux': + return result + + +@pytest.mark.filterwarnings('ignore') +@pytest.mark.parametrize('spectra_fluxunit', test_cases) +def test_cubeviz_collapse_fluxunits(spectrum1d_cube_custom_fluxunit, spectra_fluxunit): + ''' Calculates line flux and checks the units for each collapse function ''' + data = spectrum1d_cube_custom_fluxunit(spectra_fluxunit) + for function in COLLAPSE_FUNCTIONS: + # Initialize Cubeviz with specific data and collapse function + cubeviz_helper = Cubeviz() + data_label = "Test Cube" + cubeviz_helper.load_data(data, data_label=data_label) + cubeviz_helper.app.get_viewer('spectrum-viewer').state.function = function + + lineflux_result = _calculate_line_flux(cubeviz_helper) + + autocollapsed_spectrum_unit = cubeviz_helper.specviz.get_spectra()[data_label + "[FLUX]" + ].flux.unit + # Futureproofing: Eventually Cubeviz autocollapse will change the flux units of the + # spectra depending on whether the spectrum was collapsed OVER SPAXELS or not. Only + # do the assertion check if we KNOW what the expected lineflux results should be + if autocollapsed_spectrum_unit in expected_lineflux_results.keys(): + assert u.Unit(lineflux_result['unit']) == expected_lineflux_results[spectra_fluxunit] + + +@pytest.mark.filterwarnings('ignore') +@pytest.mark.parametrize('test_case', unit_flux_gaussian_test_cases) +def test_unit_gaussian(specviz_helper, test_case): + ''' + Test an Area 1 Gaussian and ensure the result returns in W/m2 + Test provided by Patrick Ogle + ''' + specviz_helper.load_data(test_case) + + lineflux_result = _calculate_line_flux(specviz_helper) + assert_quantity_allclose(float(lineflux_result['result']) * u.Unit(lineflux_result['unit']), + 1*u.Unit('W/m2')) + + +@pytest.mark.filterwarnings('ignore') +def test_unit_gaussian_mixed_units_per_steradian(specviz_helper): + ''' + A special unit test of Area 1 with mixed units. Should return W/m2sr + Test provided by Patrick Ogle + ''' + # unit-flux gaussian in wavelength space, mixed units, per steradian + lam_a = np.arange(1, 2, 0.001)*u.Angstrom + flx_wave = _gauss_with_unity_area(lam_a.value, mn, sig)*1E3*u.erg/u.s/u.cm**2/u.Angstrom/u.sr + fl_wave = Spectrum1D(spectral_axis=lam_a, flux=flx_wave) + + specviz_helper.load_data(fl_wave) + lineflux_result = _calculate_line_flux(specviz_helper) + assert_quantity_allclose(float(lineflux_result['result']) * u.Unit(lineflux_result['unit']), + 1*u.Unit('W/(m2sr)')) diff --git a/jdaviz/conftest.py b/jdaviz/conftest.py index 217d3c5f1e..28f2dd1e53 100644 --- a/jdaviz/conftest.py +++ b/jdaviz/conftest.py @@ -156,9 +156,9 @@ def spectrum_collection(spectrum1d): return result -@pytest.fixture -def spectrum1d_cube(): - flux = np.arange(16).reshape((2, 2, 4)) * u.Jy +def _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy): + + flux = np.arange(16).reshape((2, 2, 4)) * fluxunit wcs_dict = {"CTYPE1": "RA---TAN", "CTYPE2": "DEC--TAN", "CTYPE3": "WAVE-LOG", "CRVAL1": 205, "CRVAL2": 27, "CRVAL3": 4.622e-7, "CDELT1": -0.0001, "CDELT2": 0.0001, "CDELT3": 8e-11, @@ -168,6 +168,16 @@ def spectrum1d_cube(): return Spectrum1D(flux=flux, wcs=w) +@pytest.fixture +def spectrum1d_cube(): + return _create_spectrum1d_cube_with_fluxunit(u.Jy) + + +@pytest.fixture +def spectrum1d_cube_custom_fluxunit(): + return _create_spectrum1d_cube_with_fluxunit + + @pytest.fixture def mos_spectrum1d(): '''