diff --git a/mapext/_version.py b/mapext/_version.py index 662c870..4cba524 100644 --- a/mapext/_version.py +++ b/mapext/_version.py @@ -1,2 +1,2 @@ -__version__ = "2024.0.15" -__date__ = "20th February 2024" \ No newline at end of file +__version__ = "2024.0.16" +__date__ = "11th March 2024" \ No newline at end of file diff --git a/mapext/core/mapClass.py b/mapext/core/mapClass.py index 178d0b1..749c581 100644 --- a/mapext/core/mapClass.py +++ b/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,7 +68,7 @@ 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): @@ -74,7 +76,8 @@ def set_map_projection(self, **kwargs): """ 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): @@ -82,7 +85,8 @@ def set_map_resolution(self, **kwargs): """ 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): @@ -90,7 +94,8 @@ def set_map_units(self, **kwargs): """ 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 \ No newline at end of file + 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 + \ No newline at end of file diff --git a/mapext/emission/emission_models.py b/mapext/emission/emission_models.py index 735ece4..c160ec8 100644 --- a/mapext/emission/emission_models.py +++ b/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. diff --git a/pyproject.toml b/pyproject.toml index e034ed0..5e04bf4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [ project ] name = "mapext" -version = "2024.0.15" +version = "2024.0.16" authors = [ { name="Thomas J. Rennie" }, ] diff --git a/tests/test_sedAnalysis.py b/tests/test_sedAnalysis.py index 69ac441..54300f4 100644 --- a/tests/test_sedAnalysis.py +++ b/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__':