Active dev

pylanetary/navigation/core.py

@@ -256,7 +256,7 @@
     https://www.astro.keele.ac.uk/jkt/codes/jktld.html 
     '''
     # Input check 
-    if np.any(mu <0) or np.any(mu>1): 
+    if np.any(mu <0.0) or np.any(mu>1.0): 
         raise ValueError('Cosine of emission angle range [0,1]')
 
     mu = np.array(mu) #necessary so when floats are passed in, the line mu[bad] = 0 doesnt complain about indexing a float
@@ -388,6 +388,10 @@
         self.sun_n = surface_normal(self.lat_g, self.lon_w, self.sun_lon)
         self.mu0 = emission_angle(self.ob_lat, self.sun_n)
 
+        # solve small numerical issue where mu = 1.0 + epsilon
+        self.mu[(self.mu > 1.0)*(self.mu < 1.00001)] = 1.0
+        self.mu0[(self.mu0 > 1.0)*(self.mu0 < 1.00001)] = 1.0
+
         avg_circumference = 2*np.pi*((self.req + self.rpol)/2.0)
         self.deg_per_px = self.pixscale_km * (1/avg_circumference) * 360
 
@@ -403,7 +407,7 @@
         return f'ModelEllipsoid instance; req={self.req}, rpol={self.rpol}'
 
 
-    def ldmodel(self, tb, a, beam=None, law='exp', mu0=None):
+    def ldmodel(self, tb, a, beam=None, law='exp', mu0=None, psf_mode='gaussian'):
         '''
         Make a limb-darkened model disk convolved with the beam
         See docstring of limb_darkening() for options
@@ -424,6 +428,8 @@
             has no effect unless law=="minnaert".
             if self.mu0 undefined, law=="minnaert", and mu0=None, then
             code will fail
+        psf_mode : str, optional, default "gaussian"
+            mode to use for convolve_with_beam(), options "airy", "gaussian"
         '''
         if mu0 is None:
             mu0=getattr(self, "mu0", None)
@@ -436,7 +442,7 @@
         if beam is None:
             return ldmodel
         else:
-            return convolve_with_beam(ldmodel, beam)   
+            return convolve_with_beam(ldmodel, beam, mode=psf_mode)   
 
 
     def zonalmodel(lats, lons, tbs, a=0.0):
@@ -573,6 +579,8 @@
                         type of limb darkening model to use
                     :beam: float, tuple, or np.array, optional. default None.
                         units pixels. see utils.convolve_with_beam
+                    :psf_mode: str, optional. default "gaussian"
+                        what beam shape to use, options "airy", "gaussian"
                     :low_thresh: float, required.
                         see documentation of skimage.feature.canny
                     :high_thresh: float, required.
@@ -594,6 +602,8 @@
                         type of limb darkening model to use
                     :beam: float, tuple, or np.array, optional. default None.
                         units pixels. see utils.convolve_with_beam
+                    :psf_mode: str, optional. default "gaussian"
+                        what beam shape to use, options "airy", "gaussian"
                     :err: float
                         per-pixel error in input image
 
@@ -628,17 +638,17 @@
         -----
         * sometimes dxerr, dyerr give unrealistic or undefined behavior
         '''
-        defaultKwargs={'err':None,'beam':None,'law':'exp'}
+        defaultKwargs={'err':None,'beam':None,'law':'exp', 'psf_mode':'gaussian'}
         kwargs = { **defaultKwargs, **kwargs }
 
         if (mode == 'convolution') or (mode == 'disk'):
             model = self.ldmodel(kwargs['tb'], kwargs['a'], law=kwargs['law'])
-            model = convolve_with_beam(model, kwargs['beam'])
+            model = convolve_with_beam(model, kwargs['beam'], mode=kwargs['psf_mode'])
             data_to_compare = self.data 
         elif mode == 'canny':
             model_planet = self.ldmodel(kwargs['tb'], kwargs['a'], law=kwargs['law'])
             if kwargs['beam'] is not None:
-                model_planet = convolve_with_beam(model_planet, kwargs['beam'])
+                model_planet = convolve_with_beam(model_planet, kwargs['beam'], mode=kwargs['psf_mode'])
 
             edges = feature.canny(self.data, sigma=kwargs['sigma'], low_threshold = kwargs['low_thresh'], high_threshold = kwargs['high_thresh'])
             model = feature.canny(model_planet, sigma=kwargs['sigma'], low_threshold = kwargs['low_thresh'], high_threshold = kwargs['high_thresh'])
@@ -734,8 +744,31 @@
         '''
 
         hdu0 = fits.PrimaryHDU()
-        hdu0.header = header
 
+        # turn ephem into OPAL-like header
+        # allow to be overwritten by input header dict
+        header_default = {
+            'SIMPLE':True,
+            'BITPIX':-32,
+            'NAXIS':0,
+            'EXTEND':True,
+            'NEXTEND':6,
+            'TRG_ROT':self.ephem['NPole_ang'],
+            'TRG_RA':self.ephem['RA'],
+            'TRG_DEC':self.ephem['DEC'],
+            'TRG_R_A':self.req, 
+            'TRG_R_B':self.rpol,
+            'TRG_LON':self.ephem['PDObsLon'],
+            'TRG_LAT':self.ephem['PDObsLat'],
+            'SUN_LAT':self.ephem['PDSunLat'],
+            'TRG_PHAS':self.ephem['alpha'],
+            'TRG_R':self.ephem['r'],
+            'TRG_D':self.ephem['delta'],
+        }
+        header_out = {**header_default, **header} #gives priority to input header dict
+        hdu0.header = fits.Header(header_out)
+
+        # set up metadata for hdul[1-5]
         hdulist = [hdu0]
         data_list = [self.data, 
                     self.lat_g, 
@@ -753,6 +786,7 @@
         date.format = 'iso'
         date = date.iso[:10]
 
+        # generate hdul[1-5] headers and data
         for i in range(len(data_list)):
 
             data = data_list[i]
@@ -769,6 +803,7 @@
             hdu = fits.ImageHDU(data=data, header=hdr, name=extnames_list[i])
             hdulist.append(hdu)
 
+        # write it
         hdul = fits.HDUList(hdulist)
         hdul.writeto(outstem, overwrite=True)
 

pylanetary/navigation/tests/test_nav_remote.py

@@ -142,11 +142,11 @@ def test_write(datadir):
     nav.mu0 = np.cos(np.deg2rad(hdul[5].data))
 
     # write it, then reload it
-    nav.write(os.path.join(datadir, 'tmp.fits'), header=hdul[0].header, flux_unit='I/F')
+    nav.write(os.path.join(datadir, 'tmp.fits'), header={}, flux_unit='I/F')
     hdul_new = fits.open(os.path.join(datadir, 'tmp.fits'))
 
     # assert header info propagated
-    assert hdul_new[0].header == hdul[0].header
+    assert np.isclose(hdul_new[0].header['TRG_LON'], hdul[0].header['TRG_LON'], rtol=1e-3)
 
     # assert data are in correct places
     for i in range(1,6):

pylanetary/utils/core.py

@@ -1,16 +1,18 @@
 # Licensed under a ??? style license - see LICENSE.rst
 import numpy as np
+from datetime import timedelta
+import importlib, yaml, importlib.resources
+from scipy.interpolate import interp1d
+from scipy import fftpack
+
 from astropy import convolution
 from astropy.coordinates import Angle
 import astropy.units as u
 from astropy.units import Quantity
-import importlib, yaml, importlib.resources
-from scipy.interpolate import interp1d
-from scipy import fftpack
-import pylanetary.utils.data as body_info
 from astropy.time import Time
 from astroquery.jplhorizons import Horizons
-from datetime import timedelta
+
+import pylanetary.utils.data as body_info
 
 '''
 To do:
@@ -363,42 +365,65 @@
     return (1/z**2)*arr.reshape(shape).mean(-1).mean(1)
 
 
-def convolve_with_beam(data, beam):
+def convolve_with_beam(data, beam, mode='gaussian'):
     '''
     Convolves input 2-D image with a Gaussian beam or an input PSF image
 
     Parameters
     ----------
     data : np.array, required.
     beam : float/int, 3-element array-like, np.array, or None, required.
-        if float/int, circular Gaussian beam assumed, and this sets the fwhm 
-        [pixels].
-        if 3-element array-like, those are (fwhm_x, fwhm_y, theta_deg) for a 2-D Gaussian
-        [pixels, pixels, degrees].
+        If float/int, this sets the fwhm [pixels] of either Airy disk or
+        circular Gaussian beam, depending on "mode" parameter.
+        If 3-element array-like, those are (fwhm_x, fwhm_y, theta_deg) for a 2-D Gaussian
+        [pixels, pixels, degrees]. In this case, "Airy" mode not supported.
         if np.array of size > 3, assumes input PSF image
-        if None, returns original data
+        if None or 0.0, simply returns original data
+    mode : str, optional, default "gaussian"
+        options "gaussian", "airy"; what beam shape to use. Case insensitive.
+        "airy" only accepts beam of dtype float 
+        or 1-element array-like (i.e., beam must be circular).
+        this parameter has no effect if beam is a 2-D array.
 
     Returns
     -------
     np.array
         beam-convolved data
     '''
-    if beam is None:
+    # allow pass through for beam of zero size
+    if (beam is None):
         return data
-    if np.array(beam).size == 1:
+
+    # check inputs  
+    mode = mode.lower()
+    if mode not in ['gaussian', 'airy']:
+        raise ValueError(f'mode {mode} not recognized; supported options are "gaussian", "airy"')
+    if (mode == 'airy') and (np.array(beam).size != 1):
+        raise ValueError(f'mode "airy" only accepts a single value for the "beam" parameter (distance to first null)')
+
+    if mode == 'airy':
+        if beam == 0.0:
+            return data
+        # convert FWHM to first-null distance
+        null = 0.5 * (2.44/1.02) * beam
+        psf = convolution.AiryDisk2DKernel(radius=null)
+    elif (mode == 'gaussian') and (np.array(beam).size == 1):
+        if beam == 0.0:
+            return data
         fwhm_x = beam
         fwhm_y = beam
         theta = 0.0
         psf = convolution.Gaussian2DKernel(fwhm_x / 2.35482004503,
-                                            fwhm_y / 2.35482004503,
-                                            Angle(theta, unit=u.deg))
-    elif np.array(beam).size == 3:
+                                                fwhm_y / 2.35482004503,
+                                                Angle(theta, unit=u.deg))
+    elif (mode == 'gaussian') and (np.array(beam).size == 3):
         (fwhm_x, fwhm_y, theta) = beam
         psf = convolution.Gaussian2DKernel(fwhm_x / 2.35482004503,
-                                            fwhm_y / 2.35482004503,
-                                            Angle(theta, unit=u.deg))
+                                                fwhm_y / 2.35482004503,
+                                                Angle(theta, unit=u.deg))
     else:
         psf = beam
+
     return convolution.convolve_fft(data, psf)
 
 

pylanetary/utils/data/Io.yaml

@@ -2,8 +2,8 @@ body:
   name: Io
   jpl_hor_id: 501
   mass: 89.3e21  # kg
-  req: 1821.  # km equatorial radius at 1 bar level
-  rpol: 1821.  # km polar radius
+  req: 1829.4  # km equatorial radius at 1 bar level
+  rpol: 1815.7  # km polar radius
   rvol: 1821.  # km volumetric mean radius
   accel_g: 1.80  # m/s^2 Gravity at 1 bar level
   num_moons: None