HinodeXRT / xrtpy

https://xrtpy.readthedocs.io
BSD 2-Clause "Simplified" License
11 stars 8 forks source link

Error Raised by `TemperatureResponseFundamental.temperature_response` #74

Closed DanRyanIrish closed 1 year ago

DanRyanIrish commented 2 years ago

Hi XRTpy team. I'm excited to see that this package is being developed. When trying it out myself I encountered a ValueError error when executing the temperature_response function for the Be-thick filter. Below is the code I wrote that raised the error as well as the traceback.

>>> from xrtpy.response import TemperatureResponseFundamental
>>> be_thick = TemperatureResponseFundamental("Be-thick", "2021-05-07")
>>> be_thick.temperature_response()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [26], in <cell line: 1>()
----> 1 be_thick_temp_resp = be_thick.temperature_response()

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/astropy/units/decorators.py:304, in QuantityInput.__call__.<locals>.wrapper(*func_args, **func_kwargs)
    302 # Call the original function with any equivalencies in force.
    303 with add_enabled_equivalencies(self.equivalencies):
--> 304     return_ = wrapped_function(*func_args, **func_kwargs)
    306 # Return
    307 ra = wrapped_signature.return_annotation

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/temperature_response.py:183, in TemperatureResponseFundamental.temperature_response(self)
    180 @u.quantity_input
    181 def temperature_response(self) -> u.DN * u.cm**5 / (u.s * u.pix):
    182     """Apply gain value to the Temperature Response in units of DN cm\ :sup:`5` s\ :sup:`-1` pix\ :sup:`-1`."""
--> 183     return self.integration() / self.ccd_gain_right

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/astropy/units/decorators.py:304, in QuantityInput.__call__.<locals>.wrapper(*func_args, **func_kwargs)
    302 # Call the original function with any equivalencies in force.
    303 with add_enabled_equivalencies(self.equivalencies):
--> 304     return_ = wrapped_function(*func_args, **func_kwargs)
    306 # Return
    307 ra = wrapped_signature.return_annotation

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/temperature_response.py:162, in TemperatureResponseFundamental.integration(self)
    160 constants = (_c_Å_per_s * _h_eV_s / self.channel_wavelength).value
    161 factors = (self.solid_angle_per_pixel / self.ev_per_electron).value
--> 162 effective_area = (self.effective_area()).value
    164 temp_resp_w_u_c = [
    165     integrate.simpson(
    166         self.spectra()[i] * effective_area * constants * factors,
   (...)
    169     for i in range(61)
    170 ]
    172 return temp_resp_w_u_c * (u.electron * u.cm**5 * (1 / u.s) * (1 / u.pix))

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/astropy/units/decorators.py:304, in QuantityInput.__call__.<locals>.wrapper(*func_args, **func_kwargs)
    302 # Call the original function with any equivalencies in force.
    303 with add_enabled_equivalencies(self.equivalencies):
--> 304     return_ = wrapped_function(*func_args, **func_kwargs)
    306 # Return
    307 ra = wrapped_signature.return_annotation

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/temperature_response.py:154, in TemperatureResponseFundamental.effective_area(self)
    152 @u.quantity_input
    153 def effective_area(self) -> u.cm**2:
--> 154     return effective_area(self.name, self.observation_date)

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/effective_area.py:449, in effective_area(filter_name, observation_date)
    447 def effective_area(filter_name, observation_date):
    448     EAP = EffectiveAreaFundamental(filter_name, observation_date)
--> 449     return EAP.effective_area()

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/astropy/units/decorators.py:304, in QuantityInput.__call__.<locals>.wrapper(*func_args, **func_kwargs)
    302 # Call the original function with any equivalencies in force.
    303 with add_enabled_equivalencies(self.equivalencies):
--> 304     return_ = wrapped_function(*func_args, **func_kwargs)
    306 # Return
    307 ra = wrapped_signature.return_annotation

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/effective_area.py:443, in EffectiveAreaFundamental.effective_area(self)
    436 @u.quantity_input
    437 def effective_area(self) -> u.cm**2:
    438     """Calculation of the Effective Area."""
    439     return (
    440         self.channel_geometry_aperture_area
    441         * self.channel_transmission
    442         * self.interpolated_CCD_contamination_transmission
--> 443         * self.interpolated_filter_contamination_transmission
    444     )

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/effective_area.py:432, in EffectiveAreaFundamental.interpolated_filter_contamination_transmission(self)
    428 @property
    429 def interpolated_filter_contamination_transmission(self):
    430     """Interpolate filter contam transmission to the wavelength."""
    431     Filter_contam_transmission = interpolate.interp1d(
--> 432         self.n_DEHP_wavelength, self.filter_contamination_transmission
    433     )
    434     return Filter_contam_transmission(self.channel_wavelength)

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/functools.py:981, in cached_property.__get__(self, instance, owner)
    979 val = cache.get(self.attrname, _NOT_FOUND)
    980 if val is _NOT_FOUND:
--> 981     val = self.func(instance)
    982     try:
    983         cache[self.attrname] = val

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/effective_area.py:400, in EffectiveAreaFundamental.filter_contamination_transmission(self)
    397 i_i = complex(0, 1)  # Define complex number
    399 # Define transfer matrix
--> 400 M = [
    401     [
    402         [
    403             np.cos(self.filterwheel_angular_wavenumber[i]),
    404             (-i_i * np.sin(self.filterwheel_angular_wavenumber[i])) / index[i],
    405         ],
    406         [
    407             -i_i * np.sin(self.filterwheel_angular_wavenumber[i]) * index[i],
    408             np.cos(self.filterwheel_angular_wavenumber[i]),
    409         ],
    410     ]
    411     for i in range(4000)
    412 ]
    414 transmittance = [
    415     2
    416     * n_o
   (...)
    423     for i in range(4000)
    424 ]
    426 return [abs(transmittance[i] ** 2) for i in range(4000)]

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/effective_area.py:403, in <listcomp>(.0)
    397 i_i = complex(0, 1)  # Define complex number
    399 # Define transfer matrix
    400 M = [
    401     [
    402         [
--> 403             np.cos(self.filterwheel_angular_wavenumber[i]),
    404             (-i_i * np.sin(self.filterwheel_angular_wavenumber[i])) / index[i],
    405         ],
    406         [
    407             -i_i * np.sin(self.filterwheel_angular_wavenumber[i]) * index[i],
    408             np.cos(self.filterwheel_angular_wavenumber[i]),
    409         ],
    410     ]
    411     for i in range(4000)
    412 ]
    414 transmittance = [
    415     2
    416     * n_o
   (...)
    423     for i in range(4000)
    424 ]
    426 return [abs(transmittance[i] ** 2) for i in range(4000)]

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/functools.py:981, in cached_property.__get__(self, instance, owner)
    979 val = cache.get(self.attrname, _NOT_FOUND)
    980 if val is _NOT_FOUND:
--> 981     val = self.func(instance)
    982     try:
    983         cache[self.attrname] = val

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/effective_area.py:321, in EffectiveAreaFundamental.filterwheel_angular_wavenumber(self)
    313 angular_wavenumber = np.array(
    314     [
    315         (2.0 * math.pi * index[i] * cos_a) / self.n_DEHP_wavelength[i]
    316         for i in range(4000)
    317     ]
    318 )
    320 # Multiply by thickness
--> 321 angular_wavenumber_thickness = angular_wavenumber * self.contamination_on_filter
    323 real_angular_wavenumber = angular_wavenumber_thickness.real
    324 imaginary_angular_wavenumber = angular_wavenumber_thickness.imag

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/xrtpy/response/effective_area.py:193, in EffectiveAreaFundamental.contamination_on_filter(self)
    187 """
    188 Thickness of the contamination layer on a filter."""
    190 interpolater = scipy.interpolate.interp1d(
    191     self.filter_data_dates_to_seconds, self.filter_data, kind="linear"
    192 )
--> 193 return interpolater(self.filter_observation_date_to_seconds)

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/scipy/interpolate/_polyint.py:78, in _Interpolator1D.__call__(self, x)
     57 """
     58 Evaluate the interpolant
     59 
   (...)
     75 
     76 """
     77 x, x_shape = self._prepare_x(x)
---> 78 y = self._evaluate(x)
     79 return self._finish_y(y, x_shape)

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/scipy/interpolate/_interpolate.py:695, in interp1d._evaluate(self, x_new)
    693 y_new = self._call(self, x_new)
    694 if not self._extrapolate:
--> 695     below_bounds, above_bounds = self._check_bounds(x_new)
    696     if len(y_new) > 0:
    697         # Note fill_value must be broadcast up to the proper size
    698         # and flattened to work here
    699         y_new[below_bounds] = self._fill_value_below

File ~/miniconda3/envs/xray-stereoscopy/lib/python3.10/site-packages/scipy/interpolate/_interpolate.py:727, in interp1d._check_bounds(self, x_new)
    724     raise ValueError("A value in x_new is below the interpolation "
    725                      "range.")
    726 if self.bounds_error and above_bounds.any():
--> 727     raise ValueError("A value in x_new is above the interpolation "
    728                      "range.")
    730 # !! Should we emit a warning if some values are out of bounds?
    731 # !! matlab does not.
    732 return below_bounds, above_bounds

ValueError: A value in x_new is above the interpolation range.

Package Versions

xrtpy: 0.1.0 astropy: 5.0.1 scipy: 1.8.0

wtbarnes commented 1 year ago

I've just run into this as well using a slightly earlier date, 2020-11-09 17:59:59.120.

Naively, I would've expected the contamination data to extend as far as there were observations available (I'm pulling this date from a Be-thin file).

joyvelasquez commented 1 year ago

Close by #77

joyvelasquez commented 1 year ago

This issue has been fixed in release v0.1.1.