adding knot size following disk profile

skycatalogs/objects/galaxy_object.py

@@ -45,6 +45,22 @@ def _get_sed(self, component=None, resolution=None):
 
         return sed, magnorm
 
+    def _get_knot_size(self,z):
+        """
+        Return the angular knot size. Knots are modelled as the same physical size
+        """
+        # Deceleration paramameter
+        q  = -0.55
+        # Angular diameter scaling approximation in pc
+        dA = (3e9/q**2)*(z*q+(q-1)*(np.sqrt(2*q*z+1)-1))/(1+z)**2*(1.4-0.53*z)
+        # Using typical knot size 250pc, convert to sigma in arcmin
+        if z<0.6:
+            return 206264.8*250/dA/2.355
+        else:
+            # Above z=0.6, fractional contribution to post-convolved size 
+            # is <20% for smallest Roman PSF size, so can treat as point source
+            return 0
+
     def get_wl_params(self):
         """Return the weak lensing parameters, g1, g2, mu."""
         gamma1 = self.get_native_attribute('shear_1')
@@ -94,9 +110,13 @@ def get_gsobject_components(self, gsparams=None, rng=None):
             if component == 'knots':
                 npoints = self.get_native_attribute('n_knots')
                 assert npoints > 0
+                z = self.get_native_attribute('redshift')
+                size = self._get_knot_size(z) # get knot size
                 obj = galsim.RandomKnots(npoints=npoints,
-                                         half_light_radius=hlr, rng=rng,
+                                         profile=obj_dict['disk'], rng=rng,
                                          gsparams=gsparams)
+                obj = galsim.Convolve(obj, galsim.Gaussian(sigma=size))
+
             else:
                 n = self.get_native_attribute(f'sersic_{component}')
                 # Quantize the n values at 0.05 so that galsim can
@@ -113,9 +133,13 @@ def get_gsobject_components(self, gsparams=None, rng=None):
             # position angle in the old code. Newer input galaxy catalogs
             # most likely will not need this adjustment
             shear = galsim.Shear(g1=-e1, g2=-e2)
-            obj = obj._shear(shear)
-            g1, g2, mu = self.get_wl_params()
-            obj_dict[component] = obj._lens(g1, g2, mu)
+            obj_dict[component] = obj._shear(shear)
+
+        # Apply lensing
+        g1, g2, mu = self.get_wl_params()
+        for component in self.subcomponents:
+            obj_dict[component] = obj_dict[component]._lens(g1, g2, mu)
+
         return obj_dict
 
     def get_observer_sed_component(self, component, mjd=None, resolution=None):