2024.0.16

* Added general free-free model (without T_e constraints)
* Bug fixes since 2024.0.15

mapext/_version.py

@@ -1,2 +1,2 @@
-__version__ = "2024.0.15"
-__date__ = "20th February 2024"
+__version__ = "2024.0.16"
+__date__ = "11th March 2024"

mapext/core/mapClass.py

@@ -7,6 +7,7 @@
 from astropy.wcs.utils import pixel_to_skycoord
 from astropy_healpix import HEALPix
 from reproject import reproject_interp,reproject_from_healpix
+from mapext import _version
 
 __all__ = ['astroMap', 'mapObj']
 
@@ -55,6 +56,7 @@ def __init__(self, load_maps='IQU', **kwargs):
         # INITIALISE MAP METADATA
         metadata = kwargs['Metadata']
         self.name = metadata.get('name', 'mapobject')
+        self.ID = metadata.get('id', self.name)
         self.freq = get_astropy_quantity(metadata.get('freq', None))
         if self.freq == None:
             raise ValueError('Key `freq` must be specified in Metadata kwarg dictionary.')
@@ -66,31 +68,34 @@ def __init__(self, load_maps='IQU', **kwargs):
             if (map_params is None) or (stokes_parameter not in load_maps):
                 setattr(self, stokes_parameter, None)
             else:
-                setattr(self, stokes_parameter, mapObj(name=f'{self.name}_{stokes_parameter}', freq=self.freq, **map_params))
+                setattr(self, stokes_parameter, mapObj(name=f'{self.ID}_{stokes_parameter}', freq=self.freq, **map_params))
         return
 
     def set_map_projection(self, **kwargs):
         """ Method to set map projection for all available maps in this object. For details see :func:`mapext.core.mapClass.mapObj.reproject`.
         """
         for comp_stokes in self.stokes:
             comp_obj = getattr(self, comp_stokes)
-            comp_obj.reproject(**kwargs)
+            if comp_obj != None:
+                comp_obj.reproject(**kwargs)
         return
 
     def set_map_resolution(self, **kwargs):
         """ Method to set map resolution for all available maps in this object. For details see :func:`mapext.core.mapClass.mapObj.smooth_to`.
         """
         for comp_stokes in self.stokes:
             comp_obj = getattr(self, comp_stokes)
-            comp_obj.smooth_to(**kwargs)
+            if comp_obj != None:
+                comp_obj.smooth_to(**kwargs)
         return
 
     def set_map_units(self, **kwargs):
         """ Method to set map units for all available maps in this object. For details see :func:`mapext.core.mapClass.mapObj.convert_units`.
         """
         for comp_stokes in self.stokes:
             comp_obj = getattr(self, comp_stokes)
-            comp_obj.smooth_to(**kwargs)
+            if comp_obj != None:
+                comp_obj.smooth_to(**kwargs)
         return
 
 
@@ -125,13 +130,12 @@ def __init__(self, name=None, freq=None, **kwargs):
         self.freq = get_astropy_quantity(freq)
         self.reso = get_astropy_quantity(kwargs.get('resolution', None))
         self.bwid = get_astropy_quantity(kwargs.get('beamwidth', str(self.reso)))
-        self.cal = kwargs.get('cal', 8)
         self.unit = get_astropy_unit(kwargs.get('unit', None))
         if self.unit in [u.K, u.mK, u.uK, u.nK]:
             self.temp = kwargs.get('temp', 'rj')
         else:
             self.temp = None
-        self.uncert_cal = kwargs.get('calibration_uncertainty', 8)
+        self.cal = kwargs.get('calibration', 8)
         self.tfield = kwargs.get('tfield', None)
         self.zaxis = kwargs.get('zaxis', None)
         self.fname = str(kwargs.get('fname', None))
@@ -183,10 +187,14 @@ def read_hpx_map(self):
         """
         hdu = fits.open(self.fname)[self.hdrno]
         coordsys = None
-        if hdu.header["COORDSYS"].lower() in coorddict:
-            coordsys = coorddict[hdu.header["COORDSYS"].lower()]
+        if "COORDSYS" in hdu.header.keys():
+            if hdu.header["COORDSYS"].lower() in coorddict:
+                coordsys = coorddict[hdu.header["COORDSYS"].lower()]
+            else:
+                coordsys = hdu.header["COORDSYS"].lower()
+
         else:
-            coordsys = hdu.header["COORDSYS"].lower()
+            coordsys = 'galactic'
         proj = HEALPix(nside=hdu.header['NSIDE'], order=hdu.header['ORDERING'], frame=coordsys)
         if self.tfield != None:
             data = np.array(hdu.data[self.tfield])
@@ -213,7 +221,7 @@ def reproject(self, wcs=None, shape_out=None, order='nearest_neighbour'):
             self.data, self.proj = new_map, wcs
         elif type(self.proj) == HEALPix:
             new_map, new_footprint = reproject_from_healpix((self.data, str(self.proj.frame.name).lower()), wcs, shape_out=shape_out, nested=True)
-            new_map[new_footprint==0] = np.nan
+            # new_map[new_footprint==0] = np.nan
             self.data, self.proj = new_map, wcs
         return
 
@@ -311,4 +319,41 @@ def get_beam_area(self):
         :return: Beam area.
         :rtype: astropy.unit.Quantity
         """
-        return np.pi * (self.bwid/2.355)**2
+        return np.pi * (self.bwid/2.355)**2
+
+    def meta_to_hdr(self):
+        unit = ''
+        if self.unit == u.K:
+            unit = f'{self.unit}_{self.temp.upper()}'
+        else:
+            unit = f'{self.unit}'
+        include = [
+            f'Resolution : {self.reso}',
+            f'Beamwidth : {self.bwid}',
+            f'Frequency : {self.freq}',
+            f'Units : {unit}',
+            f'calUncertainty : {self.uncert_cal}']
+        return include
+
+    def save_out(self, filename=None, dir = '.'):
+        """Method to save out the current state of a map.
+        """
+        hdu = fits.PrimaryHDU(self.data, header=self.proj.to_header())
+        hdul = fits.HDUList(hdu)
+
+        hdul[0].header['COMMENT'] = f'Processed with MAPEXT {_version.__version__}'
+        hdul[0].header['COMMENT'] = '-'*80
+        for _ in self.meta_to_hdr():
+            hdul[0].header['COMMENT'] = _
+
+        if filename is not None:
+            outfilename = filename
+        else:
+            outfilename = '{}-{}.fits'.format(self.name, self.reso)
+        try:
+            hdul.writeto(f'{dir}/{outfilename}', overwrite=True)
+        except:
+            print('File already exists - not overwritten')
+        del hdul
+        return
+

mapext/emission/emission_models.py

@@ -3,7 +3,7 @@
 
 __all__ = [
     'synchrotron_1comp',
-    'freeFree_7000k',
+    'freeFree','freeFree_7000k',
     'ame_lognormal',
     'thermalDust',
 ]
@@ -18,8 +18,9 @@ class FittableEmissionModel(FittableModel):
     n_inputs = 2
     n_outputs = 1
 
-    def fit_deriv(self, *args):
-        return self.deriv(self,*args)
+    @staticmethod
+    def fit_deriv(*args):
+        return self.deriv(self, *args)
 
 # Synchrotron emission model
 class synchrotron_1comp(FittableEmissionModel):
@@ -46,7 +47,7 @@ class synchrotron_1comp(FittableEmissionModel):
                           description='Synchrotron spectral index')
 
     @staticmethod
-    def evaluate(nu, area, synch_S1:float, synch_alp:float):
+    def evaluate(nu, area, synch_S1, synch_alp):
         """Evaluate the emission model.
 
         :param nu: Frequency in GHz
@@ -147,6 +148,91 @@ def deriv(self,nu, area, ff_em):
         d_ff_em = 2.0 * phys_const['k'] * area * np.power(np.multiply(nu, 1e9) / phys_const['c'], 2) * T_ff_derivative_ff_em * 1e26
         return [d_ff_em]
 
+# Free-free emission model (Draine 2011)
+class freeFree(FittableEmissionModel):
+    """Emission model for free-free UCHii emission 
+    
+    :param ff_em: Free-free emission measure
+    :type ff_em: float
+    :param ff_Te: Free-free electron temperature
+    :type ff_Te: float
+
+    :Formulism:
+        .. math::
+            g^\\mathrm{ff}_\\nu = \\ln\\left(\\exp\\left\\{5.90 - \\frac{\\sqrt{3}}{\\pi}\\ln\\left(\\left[ \\frac{\\nu}{\\mathrm{GHz}} \\right] \\left[\\frac{T_e}{10^4\\,\\mathrm{K}}\\right] ^\\frac{3}{2}\\right)\\right\\} + 2.71828\\right)
+        .. math::
+            \\tau^\\mathrm{ff}_\\nu = 5.468\\times 10^{-2} \\cdot T_e^{-\\frac{3}{2}} \\left[ \\frac{\\nu}{\\mathrm{GHz}} \\right]^{-2} \\left[\\frac{EM}{\\mathrm{pc\\,cm}^-6}\\right]  g^\\mathrm{ff}_\\nu
+        .. math::
+            T^\\mathrm{ff}_\\nu = T_e \\left(1-e^{-\\tau^\\mathrm{ff}_\\nu}\\right)
+        .. math::
+            S^\\mathrm{ff}_\\nu = \\frac{2k_B\\Omega\\nu^2}{c^2} T^\\mathrm{ff}_\\nu
+        For details see `Draine (2011)`_ .
+    
+    :note: Either all or none of input ``nu``, ``area`` and ``ff_em`` must be provided consistently with compatible units or as unitless numbers.
+    
+    .. _Draine (2011): https://ui.adsabs.harvard.edu/abs/2011piim.book.....D/abstract
+    """
+    ff_em = Parameter(default=100,
+                      min=0,
+                      description="Free-free emission measure")
+    ff_Te = Parameter(default=7000,
+                      min=0,
+                      description="Free-free electron temperature")
+
+    @staticmethod
+    def evaluate(nu, area, ff_em, ff_Te):
+        """Evaluate the emission model.
+
+        :param nu: Frequency in GHz
+        :type nu: float or numpy.ndarray
+        :param area: Beam area in steradians
+        :type name: float or numpy.ndarray
+        :param ff_em: Free-free emission measure
+        :type ff_em: float
+        :param ff_Te: Free-free electron temperature
+        :type ff_Te: float
+
+        :return: Evaluated function
+        :rtype: float or numpy.ndarray
+        """
+        g = np.log( np.exp( 5.90 - ( np.sqrt(3)/np.pi * np.log( nu * ((ff_Te/1e4)**1.5) ) ) ) + 2.71828 )
+        tau = 5.468e-2 * (ff_Te**-1.5) * (nu**-2) * ff_em * g
+        T_ff = ff_Te * (1 - np.exp(-1 * tau))
+        S = 2. * phys_const['k'] * area * np.power(np.multiply(nu,1e9),2)  / phys_const['c']**2 * T_ff * 1e26
+        return S
+
+    def deriv(self, nu, area, ff_em, ff_Te):
+        """Evaluate the first derivitives of emission model with respect to input parameters.
+
+        :param nu: Frequency in GHz
+        :type nu: float or numpy.ndarray
+        :param area: Beam area in steradians
+        :type name: float or numpy.ndarray
+        :param ff_em: Free-free emission measure
+        :type ff_em: float
+        :param ff_Te: Free-free electron temperature
+        :type ff_Te: float
+
+        :return: Evaluated function
+        :rtype: float or numpy.ndarray
+        """
+        g = np.log( np.exp( 5.90 - ( np.sqrt(3)/np.pi * np.log( nu * ((ff_Te/1e4)**1.5) ) ) ) + 2.71828 )
+        tau = 5.468e-2 * (ff_Te**-1.5) * (nu**-2) * ff_em * g
+        T_ff = ff_Te * (1 - np.exp(-1 * tau))
+        S = 2. * phys_const['k'] * area * np.power(np.multiply(nu,1e9),2)  / phys_const['c']**2 * T_ff * 1e26
+
+        dtau_dem = tau/ff_em
+        dTff_dem = T_ff*dtau_dem/(np.exp(tau) - 1)
+        dS_dem = S*dTff_dem/T_ff
+
+        dg_dTe = 225693 / ((ff_Te * ((nu * (ff_Te**1.5))**(np.sqrt(3)/np.pi))) + (272908*ff_Te))
+        dtau_dTe = (tau*dg_dTe/g) - (tau*1.5/ff_Te)
+        dTff_dTe = (T_ff*dtau_dTe/(np.exp(tau) - 1)) + T_ff/ff_Te
+        dS_dTe = S*dTff_dTe/T_ff
+
+        return [dS_dem, dS_dTe]
+
+
 # AME lognormal
 class ame_lognormal(FittableEmissionModel):
     """Emission model for an AME lognormal source.

pyproject.toml

@@ -1,6 +1,6 @@
 [ project ]
 name = "mapext"
-version = "2024.0.15"
+version = "2024.0.16"
 authors = [
   { name="Thomas J. Rennie" },
 ]

tests/test_sedAnalysis.py

@@ -56,8 +56,6 @@ def test_sedfit(self):
         sigmas = np.sqrt(src.flux['Sv_e']**2 + S_std**2)
         deltas_scaled = deltas/sigmas
 
-        print(deltas_scaled)
-
         self.assertTrue(np.all(np.abs(deltas_scaled)<1,axis=0))
 
 if __name__ == '__main__':