Merging the interpolated temperature function onto the main branch (#692)

* Adding the shell_interpolator following Das2025+

* temp_interp_everywhere: unify interpolation toggle and bug fixes

- standardize interpolation activation on `use_interpolated_temperature` across `HotRegion`, `Elsewhere`, and `Everywhere`
- keep `mycoolgrid` as a legacy alias with conflict checks to avoid ambiguous toggle state
- fix ceding-region assignment bug in `HotRegion` so interpolated temperatures update `_cede_cellParamVecs` (not superseding cells)
- preserve `first_spot`/`second_spot` compatibility and pass them through the shell interpolator flow
- align callsites with shell interpolator API (`read_regrid()` / `read_average_shell(...)`) using object attributes for shared config
- add robust NaN handling in interpolation (nearest fallback + finite clamp + `T_everywhere` floor)

* Cleaning up and adding comments...

* Phi-grid definition changed to match temperature map interpolation assumptions.

* The error message for using the old temperature grid argument name added also to Elsewhere and Everywhere.

* Changelog items added.

* Cleaning of an extra line.

* Mentioning "Adding docs" in the todo page.

---------

Co-authored-by: Pushpita Das <pushpita@Pushpitas-MacBook-Pro-9.local>
Co-authored-by: Tuomo Salmi <tuomo.salmi@helsinki.fi>

Author: 3 people

changelog.d/692.added.rst

@@ -0,0 +1 @@
+Added interpolated temperature functionality when using a pre-calculated temperature map.

changelog.d/692.attribution.rst

@@ -0,0 +1 @@
+Pushpita Das

changelog.d/692.changed.rst

@@ -0,0 +1 @@
+Phi-grid definition changed in cellmesh/global_mesh.pyx to match temperature map interpolation assumptions.

docs/source/todo.rst

@@ -13,4 +13,5 @@ We are currently working on the following upgrades:
 * Having the option of specifying EOS model-informed mass-radius priors.
 * Testing the robustness of universal relations used in the code.
 * Improving plotting and graphics routines to visualize results.
-* Defining and implementing code standards
+* Defining and implementing code standards.
+* Providing documentation for the interpolated temperature map functionality implemented in v3.3.0.

xpsi/Elsewhere.py

@@ -6,6 +6,7 @@
 
 from xpsi.Parameter import Parameter
 from xpsi.ParameterSubspace import ParameterSubspace
+from xpsi.shell_interpolator import Temp_Interpolator_shells
 
 class AtmosError(xpsiError):
     """ Raised if the numerical atmosphere data were not preloaded. """
@@ -91,8 +92,39 @@ class Elsewhere(ParameterSubspace):
         on a surface (accounting for the surface tilt due to rotation), then
         the iteration over image orders terminates.
 
+    :param bool use_interpolated_temperature:
+        If ``True``, replace the uniform elsewhere temperature field with
+        temperatures interpolated from ``filename`` onto the X-PSI elsewhere
+        mesh.
+
+    :param bool myelsewhere:
+        Region selector passed to the interpolation helper for elsewhere-only
+        shell data handling.
+
+    :param float T_everywhere:
+        Lower floor applied to interpolated log10 temperatures after
+        interpolation.
+
+    :param int coderes:
+        Resolution of the input shell grid. The file is assumed to represent
+        a ``coderes x coderes`` angular mesh.
+
+    :param str filename:
+        Path to the shell data file used when
+        ``use_interpolated_temperature=True``.
+
+    :param bool bhac_data:
+        If ``True``, derive temperatures from BHAC shell fluxes; otherwise
+        read temperatures directly from the file.
+
     """
     required_names = ['elsewhere_temperature (if no custom specification)']
+    optional_names = ['use_interpolated_temperature',
+                      'myelsewhere',
+                      'T_everywhere',
+                      'coderes',
+                      'filename',
+                      'bhac_data']
 
     def __init__(self,
                  sqrt_num_cells = 64,
@@ -101,16 +133,39 @@ def __init__(self,
                  values = None,
                  atm_ext="BB",
                  custom = None,
-                 image_order_limit = None):
+                 image_order_limit = None,
+                 use_interpolated_temperature = None,
+                 myelsewhere = False,
+                 T_everywhere = 5.5,
+                 coderes = 512,
+                 filename = False,
+                 bhac_data = True,
+                 **kwargs):
+
+        if 'mycoolgrid' in kwargs:
+            raise TypeError("The 'mycoolgrid' argument has been removed; "
+                            "use 'use_interpolated_temperature' instead.")
 
         self.sqrt_num_cells = sqrt_num_cells
         self.num_rays = num_rays
         self.image_order_limit = image_order_limit
         self.atm_ext = atm_ext
 
+        requested_interp = use_interpolated_temperature
+
+        self.use_interpolated_temperature = bool(requested_interp)
+        self.myelsewhere = myelsewhere
+        self.T_everywhere = T_everywhere
+        self.coderes = coderes
+        self.filename = filename
+        self.bhac_data = bhac_data
+
         if bounds is None: bounds = {}
         if values is None: values = {}
 
+        if self.use_interpolated_temperature and self.filename is False:
+            raise ValueError('Filename not given for interpolation')
+
         if not custom: # setup default temperature parameter
             T = Parameter('elsewhere_temperature',
                           strict_bounds = (3.0, 7.6), # very cold --> very hot
@@ -259,29 +314,74 @@ def _compute_cellParamVecs(self, *args):
             An *ndarray[n,n]* of mesh-point colatitudes.
 
         """
-        if args: # hot region mesh shape information
-            cellParamVecs = _np.ones((args[0].shape[0],
-                                      args[0].shape[1],
-                                      len(self.vector)+1),
-                                     dtype=_np.double)
+        if not self.use_interpolated_temperature:
+            if args: # hot region mesh shape information
+                cellParamVecs = _np.ones((args[0].shape[0],
+                                          args[0].shape[1],
+                                          len(self.vector)+1),
+                                         dtype=_np.double)
 
-            # get self.vector because there may be fixed variables
-            # that also need to be directed to the integrators
-            # for intensity evaluation
-            cellParamVecs[...,:-1] *= _np.array(self.vector)
+                # get self.vector because there may be fixed variables
+                # that also need to be directed to the integrators
+                # for intensity evaluation
+                cellParamVecs[...,:-1] *= _np.array(self.vector)
 
-            return cellParamVecs
+                return cellParamVecs
 
-        else:
-            self._cellParamVecs = _np.ones((self._theta.shape[0],
-                                            self._theta.shape[1],
-                                            len(self.vector)+1),
-                                           dtype=_np.double)
+            else:
+                self._cellParamVecs = _np.ones((self._theta.shape[0],
+                                                self._theta.shape[1],
+                                                len(self.vector)+1),
+                                               dtype=_np.double)
+
+                self._cellParamVecs[...,:-1] *= _np.array(self.vector)
 
-            self._cellParamVecs[...,:-1] *= _np.array(self.vector)
+                for i in range(self._cellParamVecs.shape[1]):
+                    self._cellParamVecs[:,i,-1] *= self._effGrav
 
-            for i in range(self._cellParamVecs.shape[1]):
-                self._cellParamVecs[:,i,-1] *= self._effGrav
+        else:
+            # Replace the uniform elsewhere temperature field with values read
+            # from the shell file and interpolated onto the X-PSI mesh.
+            shell_interp = Temp_Interpolator_shells()
+            shell_interp.coderes = self.coderes
+            shell_interp.filename = self.filename
+            shell_interp.bhac_shell_avg = self.bhac_data
+            shell_interp.elsewhere_xpsi = self.myelsewhere
+            data_avg = shell_interp.read_average_shell(bhac_shell_avg=self.bhac_data)
+
+            if args: # hot region mesh shape information
+                shell_interp.xpsi_theta = args[0]
+                shell_interp.xpsi_phi = self._phi
+
+                cellParamVecs = _np.ones((args[0].shape[0],
+                                          args[0].shape[1],
+                                          len(self.vector)+1),
+                                         dtype=_np.double)
+
+                temperature_interpolated = shell_interp.temp_interpolation_avg(
+                    thetacode=data_avg[0],
+                    phicode=data_avg[1],
+                    Tempcode=data_avg[2],
+                    T_everywhere=self.T_everywhere)
+                cellParamVecs[:,:,0] = temperature_interpolated[:,:]
+                return cellParamVecs
+
+            else:
+                shell_interp.xpsi_theta = self._theta
+                shell_interp.xpsi_phi = self._phi
+                self._cellParamVecs = _np.ones((self._theta.shape[0],
+                                                self._theta.shape[1],
+                                                len(self.vector)+1),
+                                               dtype=_np.double)
+
+                temperature_interpolated = shell_interp.temp_interpolation_avg(
+                    thetacode=data_avg[0],
+                    phicode=data_avg[1],
+                    Tempcode=data_avg[2],
+                    T_everywhere=self.T_everywhere)
+                self._cellParamVecs[:,:,0] = temperature_interpolated[:,:]
+                for i in range(self._cellParamVecs.shape[1]):
+                    self._cellParamVecs[:,i,-1] *= self._effGrav
 
     def embed(self, spacetime, threads):
         """ Embed the photosphere elsewhere. """

xpsi/Everywhere.py

@@ -6,6 +6,7 @@
 
 from xpsi.Parameter import Parameter
 from xpsi.ParameterSubspace import ParameterSubspace
+from xpsi.shell_interpolator import Temp_Interpolator_shells
 
 class AtmosError(xpsiError):
     """ Raised if the numerical atmosphere data were not preloaded. """
@@ -134,9 +135,41 @@ class Everywhere(ParameterSubspace):
         default global mesh (i.e., no subclassing) results in a time-invariant
         signal.
 
+    :param bool use_interpolated_temperature:
+        If ``True``, replace the uniform global temperature field with
+        temperatures interpolated from ``filename`` onto the X-PSI mesh.
+
+    :param bool myeverywhere:
+        Legacy alias for enabling interpolated temperatures.
+
+    :param bool everywhere_flag:
+        Legacy alias for enabling interpolated temperatures.
+
+    :param float T_everywhere:
+        Lower floor applied to interpolated log10 temperatures after
+        interpolation.
+
+    :param int coderes:
+        Resolution of the input shell grid. The file is assumed to represent
+        a ``coderes x coderes`` angular mesh.
+
+    :param str filename:
+        Path to the shell data file used when
+        ``use_interpolated_temperature=True``.
+
+    :param bool bhac_data:
+        If ``True``, derive temperatures from BHAC shell fluxes; otherwise
+        read temperatures directly from the file.
+
     """
     required_names = ['temperature (if no custom specification)']
-
+    optional_names = ['use_interpolated_temperature',
+                      'myeverywhere (legacy)',
+                      'everywhere_flag (legacy)',
+                      'T_everywhere',
+                      'coderes',
+                      'filename',
+                      'bhac_data']
     def __init__(self,
                  time_invariant,
                  bounds = None,
@@ -150,7 +183,19 @@ def __init__(self,
                  beam_opt=0,
                  custom = None,
                  image_order_limit = None,
-                 _integrator_toggle = False):
+                 _integrator_toggle = False,
+                 use_interpolated_temperature = None,
+                 myeverywhere = None,
+                 everywhere_flag = None,
+                 T_everywhere = 5.5,
+                 coderes = 512,
+                 filename=False,
+                 bhac_data = True,
+                 **kwargs):
+
+        if 'mycoolgrid' in kwargs:
+            raise TypeError("The 'mycoolgrid' argument has been removed; "
+                            "use 'use_interpolated_temperature' instead.")
 
         self.num_rays = num_rays
 
@@ -194,6 +239,34 @@ def __init__(self,
         self.time_invariant = time_invariant
         self._integrator_toggle = _integrator_toggle
 
+        requested_interp = use_interpolated_temperature
+        legacy_flags = [v for v in (myeverywhere, everywhere_flag) if v is not None]
+        if legacy_flags:
+            if any(bool(v) != bool(legacy_flags[0]) for v in legacy_flags[1:]):
+                raise ValueError('Conflicting legacy interpolation flags: '
+                                 'myeverywhere and everywhere_flag.')
+
+            legacy_flag = bool(legacy_flags[0])
+            if requested_interp is not None and legacy_flag != bool(requested_interp):
+                raise ValueError('Conflicting interpolation flags: '
+                                 'use_interpolated_temperature and legacy '
+                                 'flags (myeverywhere/everywhere_flag).')
+        else:
+            legacy_flag = None
+
+        if requested_interp is None:
+            requested_interp = legacy_flag if legacy_flag is not None else False
+
+        self.use_interpolated_temperature = bool(requested_interp)
+        self.myeverywhere = self.use_interpolated_temperature
+        self.T_everywhere = T_everywhere
+        self.coderes = coderes
+        self.filename = filename
+        self.bhac_data = bhac_data
+        
+        if self.use_interpolated_temperature and self.filename is False:
+            raise ValueError('Filename not given for interpolation')
+        
     @property
     def atm_ext(self):
         """ ... """
@@ -458,8 +531,25 @@ def _compute_cellParamVecs(self):
                                         2),
                                        dtype=_np.double)
 
-        self._cellParamVecs[...,:-1] *= _np.array([self.vector[0]])
-
+        if not self.use_interpolated_temperature:
+            self._cellParamVecs[...,:-1] *= _np.array([self.vector[0]])
+        else:
+            # Replace the uniform everywhere temperature field with values
+            # read from the shell file and interpolated onto the X-PSI mesh.
+            shell_interp = Temp_Interpolator_shells()
+            shell_interp.bhac_shell_avg = self.bhac_data
+            shell_interp.coderes= self.coderes
+            shell_interp.filename = self.filename
+            shell_interp.everywhere_xpsi= self.myeverywhere
+            shell_interp.xpsi_theta = self._theta
+            shell_interp.xpsi_phi = self._phi
+            data_avg = shell_interp.read_average_shell(bhac_shell_avg=self.bhac_data)
+
+            Temperature_interpolated = shell_interp.temp_interpolation_avg(thetacode=data_avg[0],phicode=data_avg[1],
+                                                                        Tempcode=data_avg[2], T_everywhere = self.T_everywhere)
+
+            self._cellParamVecs[:,:,0] = Temperature_interpolated[:,:]
+            
         for i in range(self._cellParamVecs.shape[1]):
             self._cellParamVecs[:,i,-1] *= self._effGrav