Add option for Hestroffer (1997) limb-darkened disk.

imageoi/GreyImg.py

@@ -402,6 +402,42 @@ def add_gaussian(self, xpos, ypos, flux, fwhm):
                 rsq = (ix - xpos)**2 + (iy - ypos)**2
                 self.image[iy][ix] += peak * exp(-4 * log(2) * rsq / fwhm**2)
 
+    def add_hestroffer_disk(self, xpos, ypos, flux, diameter, alpha):
+        """Add a circular limb-darkened disk component to the current image.
+
+        Uses the ad-hoc model of Hestroffer (1997) A&A 327, 199.
+
+        Args:
+          xpos (float): Component position on first (fast) FITS axis /pix.
+          ypos (float): Component position on second (slow) FITS axis /pix.
+          flux (float): Integrated flux of component.
+          diameter (float): Disk diameter /pix.
+          alpha (float): Limb-darkening parameter.
+
+        Example:
+
+        >>> import numpy as np
+        >>> from imageoi.GreyImg import GreyImg
+        >>> img = GreyImg('test', 64, 64, 0.25)
+        >>> img.add_hestroffer_disk(12.0, 37.0, 0.5, 11, 0.5)
+        >>> np.all(np.argmax(img.image, axis=1) == 12)
+        True
+        >>> np.all(np.argmax(img.image, axis=0) == 37)
+        True
+        >>> img.add_hestroffer_disk(13.5, 42.8, 0.25, 11, 2.2)
+        >>> np.abs(np.sum(img.image) - 0.75) < 1e-2
+        True
+
+        """
+        peak = flux * (4 + 2*alpha) / (pi*diameter**2)
+        radius = diameter / 2
+        for iy in range(self.naxis2):
+            for ix in range(self.naxis1):
+                r = np.sqrt((ix - xpos)**2 + (iy - ypos)**2)
+                if r <= radius:
+                    mu = np.sqrt(1.0 - (2.0*r/diameter)**2)
+                    self.image[iy][ix] += peak * mu**alpha
+
 
 if __name__ == "__main__":
     import doctest

imageoi/imgen/__main__.py

@@ -24,6 +24,9 @@ def generate(args):
     elif args.modeltype == 'gaussian':
         img.add_gaussian(args.naxis1 / 2, args.naxis1 / 2, 1.0,
                          args.modelwidth / args.pixelsize)
+    elif args.modeltype == 'ld':
+        img.add_hestroffer_disk(args.naxis1 / 2, args.naxis1 / 2, 1.0,
+                                args.modelwidth / args.pixelsize, args.ldalpha)
     img.make_primary_hdu().writeto(args.imagefile, overwrite=args.overwrite)
 
 
@@ -53,10 +56,12 @@ def create_parser():
     parser.add_argument('pixelsize', type=float,
                         help='Pixel size /mas')
     parser.add_argument('-mt', '--modeltype', default='blank',
-                        choices=['blank', 'dirac', 'uniform', 'gaussian'],
+                        choices=['blank', 'dirac', 'uniform', 'gaussian', 'ld'],
                         help='Image model type')
     parser.add_argument('-mw', '--modelwidth', type=float, default=10.0,
                         help='Initial image model width /mas')
+    parser.add_argument('--ldalpha', type=float, default=1.0,
+                        help='Limb-darkening parameter for -mt ld')
     return parser
 
 

tests/test_imgen.py

@@ -49,6 +49,22 @@ def test_generate(self):
             self.assertAlmostEqual(imageHdu.header['CDELT2'],
                                    pixelsize * MAS_TO_DEG)
 
+    def test_ld(self):
+        """Test Hestroffer limb-darkened disk"""
+        naxis1 = 128
+        pixelsize = 0.5
+        args = self.parser.parse_args(['--overwrite', self.tempResult.name,
+                                       str(naxis1), str(pixelsize), '-mt=ld'])
+        generate(args)
+        args = self.parser.parse_args(['--overwrite', self.tempResult.name,
+                                       str(naxis1), str(pixelsize), '-mt=ld',
+                                       '-mw=10.5'])
+        generate(args)
+        args = self.parser.parse_args(['--overwrite', self.tempResult.name,
+                                       str(naxis1), str(pixelsize), '-mt=ld',
+                                       '-mw=13.2', '--ldalpha=2.0'])
+        generate(args)
+
 
 if __name__ == '__main__':
     unittest.main()