Initial commit

1 files changed, 1399 insertions, 0 deletions

diff --git a/noao/artdata/mktemplates.xNEW b/noao/artdata/mktemplates.xNEW
new file mode 100644
index 00000000..048abf21
--- /dev/null
+++ b/noao/artdata/mktemplates.xNEW
@@ -0,0 +1,1399 @@
+include	<error.h>
+include	<imhdr.h>
+include	<math.h>
+include	<math/iminterp.h>
+
+# Template data structure
+define	LEN_MKT		18
+define	MKT_PROF	Memi[$1]	# Pointer to profile
+define	MKT_MSI		Memi[$1+1]	# MSI interpolation pointer
+define	MKT_NXM		Memi[$1+2]	# Number of X points in model
+define	MKT_NYM		Memi[$1+3]	# Number of Y points in model
+define	MKT_F		Memr[P2R($1+4)]	# Fraction of total flux in profile
+define	MKT_SCALE	Memr[P2R($1+5)]	# Radius scale
+
+define	MKT_NALLOC	Memi[$1+6]	# Allocated space for saved templates
+define	MKT_N		Memi[$1+7]	# Number of saved templates
+
+define	MKT_DATA	Memi[$1+8]	# Data pointer
+define	MKT_PTRS	Memi[$1+9]	# Data pointers
+define	MKT_NX		Memi[$1+10]	# Number of X pixels
+define	MKT_NY		Memi[$1+11]	# Number of Y pixels
+define	MKT_XC		Memi[$1+12]	# Subpixel X center
+define	MKT_YC		Memi[$1+13]	# Subpixel Y center
+define	MKT_FLUX	Memi[$1+14]	# Flux
+define	MKT_R		Memi[$1+15]	# Radius
+define	MKT_AR		Memi[$1+16]	# Axial ratio
+define	MKT_PA		Memi[$1+17]	# Position angle
+
+define	NALLOC	25		# Allocation block for saved templates
+define	NPROF	5001		# Profile length
+define	NY	11		# BINPROF binning parameter
+
+
+# MKT_INIT    -- Initialize template memory.
+# MKT_FREE    -- Free template memory.
+# MKT_SAVE    -- Save a template
+# MKT_GET     -- Get a template
+# MKT_STAR    -- Set star and seeing templates.
+# MKT_OBJECT  -- Set object profiles.
+# MKT_GOBJECT -- Get image raster.
+# MKT_BINPROF -- Bin intensity profile
+# MKT_GSTAR   -- Get the precomputed template with center nearest (x,y).
+# MKT_PROFILE -- Make template from profile.
+# MKT_MSI     -- Make template from image interpolation function.
+
+
+# MKT_INIT -- Initialize template memory.
+# The symbol table is used as a simple way to store the object types by name.
+
+procedure mkt_init ()
+
+int	clgeti()
+real	clgetr()
+
+pointer	stopen()
+include	"mktemplates.com"
+
+begin
+	nxc = clgeti ("nxc")
+	nyc = clgeti ("nyc")
+	nxssub = clgeti ("nxsub")
+	nyssub = clgeti ("nysub")
+	nxgsub = clgeti ("nxgsub")
+	nygsub = clgeti ("nygsub")
+	dynrange = clgetr ("dynrange")
+	psfrange = clgetr ("psfrange")
+	stp = stopen ("mkt", 10, 10, 10*SZ_FNAME)
+end
+
+
+# MKT_FREE -- Free template memory.
+
+procedure mkt_free ()
+
+int	i
+pointer sym, mkt, sthead(), stnext()
+include	"mktemplates.com"
+
+begin
+	# For each object type free the profile and image interpolator data,
+	# the last unsaved data buffer, all saved templates, and the object
+	# structure.  Finally free the symbol table.
+
+	for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {
+	    mkt = Memi[sym]
+	    if (mkt != NULL) {
+	        call mfree (MKT_PROF(mkt), TY_REAL)
+		if (MKT_MSI(mkt) != NULL)
+		   call msifree (MKT_MSI(mkt))
+	        call mfree (MKT_DATA(mkt), TY_REAL)
+		if (MKT_NALLOC(mkt) > 0) {
+		    do i = 0, MKT_N(mkt)-1
+		        call mfree (Memi[MKT_PTRS(mkt)+i], TY_REAL)
+		    call mfree (MKT_PTRS(mkt), TY_POINTER)
+		    call mfree (MKT_NX(mkt), TY_INT)
+		    call mfree (MKT_NY(mkt), TY_INT)
+		    call mfree (MKT_XC(mkt), TY_REAL)
+		    call mfree (MKT_YC(mkt), TY_REAL)
+		    call mfree (MKT_FLUX(mkt), TY_REAL)
+		    call mfree (MKT_R(mkt), TY_REAL)
+		    call mfree (MKT_AR(mkt), TY_REAL)
+		    call mfree (MKT_PA(mkt), TY_REAL)
+		}
+	        call mfree (mkt, TY_STRUCT)
+	    }
+	}
+	call stclose (stp)
+end
+
+
+# MKT_SAVE -- Save a template
+# If a template may be used more than once it may be saved upon direction of
+# the user in the object list.  Otherwise the last unsaved template is
+# freed for the next object.
+
+procedure mkt_save (mkt, data, nx, ny, xc, yc, flux, r, ar, pa, save)
+
+pointer	mkt		#I MKT pointer
+pointer	data		#I Pointer to template data
+int	nx, ny		#I Size of template
+real	xc, yc		#I Subpixel center
+real	flux		#I Flux
+real	r		#I Effective radius
+real	ar		#I Axial ratio
+real	pa		#I Position angle
+int	save		#I Save data?
+
+int	i
+
+begin
+	if (save == NO) {
+	    MKT_DATA(mkt) = data
+	    return
+	}
+
+	if (MKT_NALLOC(mkt) == 0) {
+	    i = NALLOC
+	    call malloc (MKT_PTRS(mkt), i, TY_POINTER)
+	    call malloc (MKT_NX(mkt), i, TY_INT)
+	    call malloc (MKT_NY(mkt), i, TY_INT)
+	    call malloc (MKT_XC(mkt), i, TY_REAL)
+	    call malloc (MKT_YC(mkt), i, TY_REAL)
+	    call malloc (MKT_FLUX(mkt), i, TY_REAL)
+	    call malloc (MKT_R(mkt), i, TY_REAL)
+	    call malloc (MKT_AR(mkt), i, TY_REAL)
+	    call malloc (MKT_PA(mkt), i, TY_REAL)
+	    MKT_NALLOC(mkt) = i
+	} else if (MKT_N(mkt) == MKT_NALLOC(mkt)) {
+	    i = MKT_NALLOC(mkt) + NALLOC
+	    call realloc (MKT_PTRS(mkt), i, TY_POINTER)
+	    call realloc (MKT_NX(mkt), i, TY_INT)
+	    call realloc (MKT_NY(mkt), i, TY_INT)
+	    call realloc (MKT_XC(mkt), i, TY_REAL)
+	    call realloc (MKT_YC(mkt), i, TY_REAL)
+	    call realloc (MKT_FLUX(mkt), i, TY_REAL)
+	    call realloc (MKT_R(mkt), i, TY_REAL)
+	    call realloc (MKT_AR(mkt), i, TY_REAL)
+	    call realloc (MKT_PA(mkt), i, TY_REAL)
+	    MKT_NALLOC(mkt) = i
+	}
+	i = MKT_N(mkt)
+	Memi[MKT_PTRS(mkt)+i] = data
+	Memi[MKT_NX(mkt)+i] = nx
+	Memi[MKT_NY(mkt)+i] = ny
+	Memr[MKT_XC(mkt)+i] = xc
+	Memr[MKT_YC(mkt)+i] = yc
+	Memr[MKT_FLUX(mkt)+i] = flux
+	Memr[MKT_R(mkt)+i] = r
+	Memr[MKT_AR(mkt)+i] = ar
+	Memr[MKT_PA(mkt)+i] = pa
+	MKT_N(mkt) = i + 1
+end
+
+
+# MKT_GET -- Get a template
+# If not a saved template just free last unsaved template.
+# If saved search for match with position, size, axial ratio, and pa.
+# Return null if not found.
+
+procedure mkt_get (mkt, data, nx, ny, xc, yc, flux, r, ar, pa, save)
+
+pointer	mkt		#I MKT pointer
+pointer	data		#O Pointer to template data
+int	nx, ny		#O Size of template
+real	xc, yc		#I Subpixel center
+real	flux		#I Flux
+real	r		#I Effective radius
+real	ar		#I Axial ratio
+real	pa		#I Position angle
+int	save		#I Get saved template?
+
+int	i
+real	f
+
+begin
+	data = NULL
+	call mfree (MKT_DATA(mkt), TY_REAL)
+	if (save == NO)
+	    return
+
+	for (i=0; i<MKT_N(mkt); i=i+1) {
+	    if (xc != Memr[MKT_XC(mkt)+i])
+		next
+	    if (yc != Memr[MKT_YC(mkt)+i])
+		next
+	    if (r != Memr[MKT_R(mkt)+i])
+		next
+	    if (ar != Memr[MKT_AR(mkt)+i])
+		next
+	    if (pa != Memr[MKT_PA(mkt)+i])
+		next
+	    data = Memi[MKT_PTRS(mkt)+i]
+	    nx = Memi[MKT_NX(mkt)+i]
+	    ny = Memi[MKT_NY(mkt)+i]
+	    f = Memr[MKT_FLUX(mkt)+i]
+	    if (f != flux) {
+	        call amulkr (Memr[data], flux/f, Memr[data], nx*ny)
+	        Memr[MKT_FLUX(mkt)+i] = flux
+	    }
+	    return
+	}
+end
+
+
+# MKT_STAR -- Define star and seeing templates.
+# The seeing template has a smaller range for efficiency.
+# THe star templates are determined once over a grid of centers and
+# then not evaluated again.
+
+pointer procedure mkt_star (name)
+
+char	name[ARB]	# Profile name or file
+
+# Star and seeing parameters obatined through CLIO.
+real	r		# Major axis sigma (pixels)
+real	ar		# Axial ratio (minor / major)
+real	pa		# Position angle (radians)
+
+int	i, j, nxm, nym, nx, ny, fd
+real	dr, flux, radius, seeing, beta, xc, yc, dxc, dyc, der[2]
+pointer	sym, mkt1, mkt2, prof, prof1, asi, msi, data, im
+
+bool	streq()
+real	clgetr(), asieval()
+int	open(), fscan(), nscan()
+pointer	immap(), imgs2r(), stfind(), stenter()
+errchk	immap, open, imgs2r, asifit, asieval, asider
+
+include	"mktemplates.com"
+
+begin
+	# Check if previously defined.
+	sym = stfind (stp, "star")
+	if (sym != NULL)
+	    return (Memi[sym])
+
+	# Select type of star profile and set intensity profile array.
+	# Compute the fraction of the total flux in the profile.
+	# Insure that the profile subsamples a subpixel.
+
+	star = NULL
+	see = NULL
+	prof = NULL
+	msi = NULL
+	if (streq (name, "gaussian")) {
+	    r = clgetr ("radius") / sqrt (log (2.))
+	    radius = sqrt (log (dynrange))
+	    seeing = sqrt (log (psfrange))
+	    nxm = max (NPROF, 1 + nint (radius*r*nxssub*nyssub))
+	    call malloc (prof, nxm, TY_REAL)
+
+	    dr = radius / (nxm - 1)
+	    do i = 0, nxm - 1
+		Memr[prof+i] = exp (-(i * dr) ** 2)
+	    flux = 1 - Memr[prof+nxm-1]
+
+	    r = sqrt (log (2.))
+	    radius = radius / r
+	    seeing = seeing / r
+	} else if (streq (name, "moffat")) {
+	    beta = clgetr ("beta")
+	    r = clgetr ("radius") / sqrt (2. ** (1/beta) - 1.)
+	    radius = sqrt ((dynrange) ** (1/beta) - 1)
+	    seeing = sqrt ((psfrange) ** (1/beta) - 1)
+	    nxm = max (NPROF, 1 + nint (radius*r*nxssub*nyssub))
+	    call malloc (prof, nxm, TY_REAL)
+
+	    dr = radius / (nxm - 1)
+	    flux = 0
+	    do i = 0, nxm - 1 {
+		r = i * dr
+		Memr[prof+i] =  1. / ((1 + r**2) ** beta)
+		flux = flux + r * Memr[prof+i]
+	    }
+
+	    # Compute the fraction of the total flux in the profile.
+	    # The last part of the total flux below is computed by expanding
+	    # (1+r**2) --> r**2 under the approximation that r >> 1.
+	    # Note that it is possible to explicitly compute the total
+	    # flux  F(total) = beta / (2 * beta - 2) (CRC 53rd edition)
+	    # I found errors in other versions of CRC for this integral!
+
+	    r = r + dr / 2
+	    xc = 2 * beta - 2
+	    flux = flux / (flux + 1. / (xc * r ** xc))
+
+	    r = sqrt (2. ** (1/beta) - 1.)
+	    radius = radius / r
+	    seeing = seeing / r
+	} else ifnoerr (im = immap (name, READ_ONLY, 0)) {
+	    iferr {
+		nxm = IM_LEN(im,1)
+		nym = IM_LEN(im,2)
+		data = imgs2r (im, 1, nxm, 1, nym)
+		call msiinit (msi, II_BILINEAR)
+		call msifit (msi, Memr[data], nxm, nym, nxm)
+	    } then
+		call erract (EA_WARN)
+	    call imunmap (im)
+
+	    flux = 1.
+	    radius = 1
+	    seeing = 0.8
+	} else ifnoerr (fd = open (name, READ_ONLY, TEXT_FILE)) {
+	    nxm = NPROF
+	    call malloc (prof1, nxm, TY_REAL)
+
+	    j = 0
+	    while (fscan (fd) != EOF) {
+	 	call gargr (flux)
+		if (nscan() < 1)
+		    next
+		if (j == nxm) {
+		    nxm = nxm + NPROF
+		    call realloc (prof1, nxm, TY_REAL)
+		}
+		Memr[prof1+j] = flux
+		j = j + 1
+	    }
+	    call close (fd)
+	    if (j == 0) {
+		call mfree (prof1, TY_REAL)
+		call error (1, "PSF template not found")
+	    }
+
+	    r = clgetr ("radius")
+	    nxm = max (NPROF, 1 + nint (r*nxssub*nyssub))
+	    call malloc (prof, nxm, TY_REAL)
+	    dr = 1. / (nxm - 1)
+	    j = j - 1
+
+	    call asiinit (asi, II_SPLINE3)
+	    call asifit (asi, Memr[prof1], j)
+	    xc = Memr[prof1]
+	    call mfree (prof1, TY_REAL)
+
+	    if (xc == 0.) {
+		flux = 0.
+		do i = 1, nxm - 1 {
+		    r = i * dr
+		    call asider (asi, 1+j*r, der, 2)
+		    Memr[prof+i] = max (0., der[2] / r)
+		    flux = flux + r * Memr[prof+i]
+		}
+		Memr[prof] = max (0., 2 * Memr[prof+1] - Memr[prof+2])
+	    } else {
+		flux = 0.
+		do i = 0, nxm - 1 {
+		    r = i * dr
+		    Memr[prof+i] = asieval (asi, 1+j*r)
+		    flux = flux + r * Memr[prof+i]
+		}
+	    }
+	    call asifree (asi)
+
+	    xc = 0.9 * flux
+	    flux = 0.
+	    for (i=1; i<nxm && flux<xc; i=i+1)
+		flux = flux + i * dr * Memr[prof+i]
+
+	    flux = 1.
+	    radius = 1.
+	    seeing = (i - 1.) * dr
+	} else
+	    call error (1, "PSF template not found")
+
+	# Set size and orientation parameters.
+	r = clgetr ("radius")
+	ar = clgetr ("ar")
+	pa = DEGTORAD (clgetr ("pa"))
+	radius = r * radius
+	seeing = r * seeing
+
+	# Compute templates with subsampling  over a grid of centers.
+	# Do this for the full star profile and a smaller region for
+	# convolving extended objects.
+
+	# Seeing kernel.
+	sym = stenter (stp, "seeing", 1)
+	call calloc (Memi[sym], LEN_MKT, TY_STRUCT)
+	mkt1 = Memi[sym]
+
+	nx = 2 * nint (seeing) + 1
+	ny = nx
+	dxc = 1. / nxc
+	dyc = 1. / nyc
+	if (prof != NULL) {
+	    nym = 1 + nint (radius * nxssub * nyssub)
+	    call malloc (prof1, nym, TY_REAL)
+	    call mkt_binprof (Memr[prof], nxm, Memr[prof1], nym, radius, nxssub)
+	    for (yc = -0.5+dyc/2; yc < 0.5; yc = yc+dyc) {
+	        for (xc = -0.5+dxc/2; xc < 0.5; xc = xc+dxc) {
+		    call malloc (data, nx*ny, TY_REAL)
+		    call mkt_profile (data, nx, ny, xc, yc, 1., Memr[prof1],
+		        nym, radius, ar, pa, nxssub, nyssub)
+		    call mkt_save (mkt1, data, nx, ny, xc, yc, 1., 0., 0., 0.,
+			YES)
+	        }
+	    }
+	} else {
+	    for (yc = -0.5+dyc/2; yc < 0.5; yc = yc+dyc) {
+	        for (xc = -0.5+dxc/2; xc < 0.5; xc = xc+dxc) {
+		    call malloc (data, nx*ny, TY_REAL)
+		    call mkt_msi (data, nx, ny, xc, yc, 1., msi, nxm, nym,
+		        radius, ar, pa, nxssub, nyssub)
+		    call mkt_save (mkt1, data, nx, ny, xc, yc, 1., 0., 0., 0.,
+			YES)
+	        }
+	    }
+	}
+
+	# Full star templates.
+	sym = stenter (stp, "star", 1)
+	call calloc (Memi[sym], LEN_MKT, TY_STRUCT)
+	mkt2 = Memi[sym]
+
+	nx = 2 * nint (radius) + 1
+	ny = nx
+	dxc = 1. / nxc
+	dyc = 1. / nyc
+	if (prof != NULL) {
+	    for (yc = 0.5+dyc/2; yc < 1.5; yc = yc+dyc) {
+	        for (xc = 0.5+dxc/2; xc < 1.5; xc = xc+dxc) {
+		    call malloc (data, nx*ny, TY_REAL)
+		    call mkt_profile (data, nx, ny, xc, yc, flux, Memr[prof1],
+		        nym, radius, ar, pa, nxssub, nyssub)
+		    call mkt_save (mkt2, data, nx, ny, xc, yc, 1., 0., 0., 0.,
+			YES)
+	        }
+	    }
+	    call mfree (prof, TY_REAL)
+	    call mfree (prof1, TY_REAL)
+	} else {
+	    for (yc = 0.5+dyc/2; yc < 1.5; yc = yc+dyc) {
+	        for (xc = 0.5+dxc/2; xc < 1.5; xc = xc+dxc) {
+		    call malloc (data, nx*ny, TY_REAL)
+		    call mkt_msi (data, nx, ny, xc, yc, flux, msi, nxm, nym,
+		        radius, ar, pa, nxssub, nyssub)
+		    call mkt_save (mkt2, data, nx, ny, xc, yc, 1., 0., 0., 0.,
+			YES)
+	        }
+	    }
+	    call msifree (msi)
+	}
+
+	see = mkt1
+	star = mkt2
+	return (star)
+end
+
+
+# MKT_OBJECT -- Set object profiles.
+
+pointer procedure mkt_object (name)
+
+char	name[ARB]	# Profile name or file
+
+int	i, j, n, nxm, nym, fd
+real	radius, r, dr, s, b, flux, der[2]
+pointer	sym, mkt, prof, asi, msi, buf, im
+
+real	asieval(), uigamma()
+int	open(), fscan(), nscan(), strncmp(), ctor()
+pointer	immap(), imgs2r(), stfind(), stenter()
+bool	streq()
+errchk	open, immap, asifit, asieval, asider
+
+include	"mktemplates.com"
+
+begin
+	# Check if previously defined.
+	sym = stfind (stp, name)
+	if (sym != NULL)
+	    return (Memi[sym])
+
+	# Select type of profile and set intensity profile array.
+	# Compute the fraction of the total flux in the profile.
+
+	prof = NULL
+	msi = NULL
+	if (strncmp (name, "sersic", 6) == 0) {
+	    i = 7
+	    if (ctor (name, i, s) == 0) {
+		call eprintf ("WARNING: Bad sersic profile syntax (%s).\n")
+		    call pargstr (name)
+		return (NULL)
+	    }
+	    n = nint (2 * (s + 0.01))
+	    s = n / 2
+call eprintf ("n = %.1f, 2n = %d\n")
+call pargr (s)
+call pargi (n)
+	    if (n < 1 || n > 20) {
+		call eprintf ("WARNING: Sersic index out of allowed range (%f).\n")
+		    call pargi (s)
+		return (NULL)
+	    }
+	    nxm = NPROF
+	    call malloc (prof, nxm, TY_REAL)
+call eprintf ("A\n")
+	    radius = log (dynrange) ** s
+call eprintf ("B\n")
+	    dr = radius / (nxm - 1)
+	    do i = 0, nxm - 1 {
+		r = (i * dr) ** (1/s)
+		Memr[prof+i] = exp (-r)
+	    }
+
+call eprintf ("C\n")
+	    flux = 1 - uigamma (n, r)
+call eprintf ("D\n")
+	    r = n - 1./3. + 4./(405.*s) + 46./(25515.*s*s)
+call eprintf ("E\n")
+	    radius = radius / r ** s
+call eprintf ("F\n")
+	} else if (streq (name, "expdisk")) {
+	    s = 1.
+	    n = nint (2 * s)
+
+	    nxm = NPROF
+	    call malloc (prof, nxm, TY_REAL)
+	    radius = log (dynrange) ** s
+	    dr = radius / (nxm - 1)
+	    do i = 0, nxm - 1 {
+		r = (i * dr) ** (1/s)
+		Memr[prof+i] = exp (-r)
+	    }
+
+	    flux = 1 - uigamma (n, r)
+	    r = n - 1./3. + 4./(405.*s) + 46./(25515.*s*s)
+	    radius = radius / r ** s
+	} else if (streq (name, "devauc")) {
+	    s = 4.
+	    n = nint (2 * s)
+	    nxm = NPROF
+	    call malloc (prof, nxm, TY_REAL)
+	    radius = log (dynrange) ** s
+	    dr = radius / (nxm - 1)
+	    do i = 0, nxm - 1 {
+		r = (i * dr) ** (1/s)
+		Memr[prof+i] = exp (-r)
+	    }
+
+	    flux = 1 - uigamma (n, r)
+	    r = n - 1./3. + 4./(405.*s) + 46./(25515.*s*s)
+	    radius = radius / r ** s
+	} else ifnoerr (im = immap (name, READ_ONLY, 0)) {
+	    iferr {
+		nxm = IM_LEN(im,1)
+		nym = IM_LEN(im,2)
+		buf = imgs2r (im, 1, nxm, 1, nym)
+		call msiinit (msi, II_BILINEAR)
+		call msifit (msi, Memr[buf], nxm, nym, nxm)
+	    } then
+		call erract (EA_WARN)
+	    call imunmap (im)
+
+	    flux = 1.
+	    radius = 1.
+	} else ifnoerr (fd = open (name, READ_ONLY, TEXT_FILE)) {
+	    nxm = NPROF
+	    call malloc (buf, nxm, TY_REAL)
+
+	    j = 0
+	    while (fscan (fd) != EOF) {
+	 	call gargr (flux)
+		if (nscan() < 1)
+		    next
+		if (j == nxm) {
+		    nxm = nxm + NPROF
+		    call realloc (buf, nxm, TY_REAL)
+		}
+		Memr[buf+j] = flux
+		j = j + 1
+	    }
+	    call close (fd)
+	    if (j == 0) {
+		call mfree (buf, TY_REAL)
+		nxm = 0
+		call error (1, "PSF template not found")
+	    }
+
+	    nxm = NPROF
+	    call malloc (prof, nxm, TY_REAL)
+	    dr = 1. / (nxm - 1)
+	    j = j - 1
+
+	    call asiinit (asi, II_SPLINE3)
+	    call asifit (asi, Memr[buf], j)
+	    s = Memr[buf]
+	    call mfree (buf, TY_REAL)
+
+	    if (s == 0.) {
+		do i = 1, nxm - 1 {
+		    r = i * dr
+		    call asider (asi, 1+j*r, der, 2)
+		    Memr[prof+i] = max (0., der[2] / r)
+		}
+		Memr[prof] = max (0., 2 * Memr[prof+1] - Memr[prof+2])
+	    } else {
+		do i = 0, nxm - 1 {
+		    r = i * dr
+		    Memr[prof+i] = asieval (asi, 1+j*r)
+		}
+	    }
+	    call asifree (asi)
+
+	    flux = 1.
+	    radius = 1.
+	} else {
+	    call eprintf ("WARNING: Object template %s not found.\n")
+		call pargstr (name)
+	    return (NULL)
+	}
+call eprintf ("G\n")
+
+	# Create the template structure if a model is defined..
+	if (prof == NULL && msi == NULL)
+	    mkt = NULL
+	else {
+	    call calloc (mkt, LEN_MKT, TY_STRUCT)
+	    MKT_PROF(mkt) = prof
+	    MKT_MSI(mkt) = msi
+	    MKT_NXM(mkt) = nxm
+	    MKT_NYM(mkt) = nym
+	    MKT_F(mkt) = flux
+	    MKT_SCALE(mkt) = radius
+call eprintf ("flux = %g, radius = %g\n")
+call pargr (MKT_F(mkt))
+call pargr (MKT_SCALE(mkt))
+do i = 0, nxm {
+call eprintf ("%d: %g\n")
+call pargi (i)
+call pargr (Memr[prof+i])
+}
+	}
+call eprintf ("H\n")
+
+	# Enter object model name in symbol table.
+	sym = stenter (stp, name, 1)
+	Memi[sym] = mkt
+	return (mkt)
+end
+
+
+# MKT_GOBJECT -- Get image raster.
+# The profile type is set by the template structure.
+
+procedure mkt_gobject (mkt, data, nx, ny, x, y, z, r, ar, pa, save)
+
+pointer	mkt		# Object template
+pointer	data		# Data
+int	nx, ny		# Size of returned data
+real	x, y		# Position of object
+real	z		# Flux of object
+real	r		# Major axis scale (pixels)
+real	ar		# Axial ratio (minor / major)
+real	pa		# Position angle (radians)
+int	save		# Use/save template?
+
+real	xc, yc, radius
+int	nprof
+pointer	prof
+
+include	"mktemplates.com"
+
+begin
+	data = NULL
+	if (mkt == NULL)
+	    return
+
+	# Stars are predefined.  Return the nearest template center.
+	# Other objects are computed with or without seeing.
+
+	if (mkt == star)
+	    call mkt_gstar (star, data, nx, ny, x, y, z)
+	else {
+	    xc = x + 0.5 - int (x + 0.5)
+	    yc = y + 0.5 - int (y + 0.5)
+	    call mkt_get (mkt, data, nx, ny, xc, yc, z, r, ar, pa, save)
+	    if (data != NULL)
+		return
+
+	    radius = r * MKT_SCALE(mkt)
+	    nx = max (3, 2 * nint (radius) + 1)
+	    ny = nx
+	    if (see != NULL) {
+	        nx = nx + Memi[MKT_NX(see)] / 2
+	        ny = ny + Memi[MKT_NY(see)] / 2
+	    }
+	    call malloc (data, nx * ny, TY_REAL)
+
+	    if (MKT_PROF(mkt) != 0) {
+	        nprof = 1 + nint (radius * nxgsub * nxgsub)
+	        call malloc (prof, nprof, TY_REAL)
+	        call mkt_binprof (Memr[MKT_PROF(mkt)], MKT_NXM(mkt),
+		    Memr[prof], nprof, radius, nxgsub)
+	        call mkt_profile (data, nx, ny, x, y, z*MKT_F(mkt), Memr[prof],
+		    nprof, radius, ar, pa, nxgsub, nygsub)
+	        call mfree (prof, TY_REAL)
+	    } else {
+	        call mkt_msi (data, nx, ny, x, y, z, MKT_MSI(mkt),
+		    MKT_NXM(mkt), MKT_NYM(mkt), radius, ar, pa,
+		    nxgsub, nygsub)
+	    }
+
+	    call mkt_save (mkt, data, nx, ny, xc, yc, z, r, ar, pa, save)
+	}
+end
+
+
+# MKT_BINPROF -- Bin intensity profile into subpixels
+
+procedure mkt_binprof (prof, nprof, prof1, nprof1, radius, nsub)
+
+real	prof[nprof]		# Input intensity profile
+int	nprof			# Number of input points
+real	prof1[nprof]		# Output binned intensity profile
+int	nprof1			# Number of output points
+real	radius			# Radius of profile
+int	nsub			# Maximum subsampling
+
+int	i, j, k, k1, k2, l, dx
+real	scale, dy, val
+
+int	debug, open()
+data	debug/0/
+
+begin
+	if (radius < 0.1) {
+	    call amovkr (1., prof1, nprof1)
+	    return
+	} else
+	    call aclrr (prof1, nprof1)
+
+	# Set binning parameters
+	scale = (nprof - 1.) / (nprof1 - 1.)
+	dx = nint ((nprof1 - 1.) / nsub / radius / 2.)
+	dy = dx / (NY - 1.)
+
+	# Bin central pixels
+	do i = -dx, 2*dx {
+	    k = abs (i)
+	    k1 = max (1, i - dx + 1) 
+	    k2 = i + dx + 1
+	    do j = 0, NY-1 {
+		if (j == 0)
+		    val = k
+		else if (k == 0)
+		    val = dy * j
+		else
+		    val = k * sqrt (1. + (dy * j / k) ** 2)
+		l = nint (scale * val + 1)
+		if (l > nprof)
+		    next
+		val = prof[l] / NY
+		do l = k1, k2
+		    prof1[l] = prof1[l] + val
+	    }
+	}
+
+	# Now bin remainder of pixels more crudely
+	do i = 2*dx+1, nprof1-1 {
+	    k1 = i - dx + 1
+	    k2 = min (nprof1, i + dx + 1)
+	    val = prof[nint (scale * i + 1)]
+	    do k = k1, k2
+		prof1[k] = prof1[k] + val
+	}
+
+	if (debug == YES) {
+	    j = open ("debug1.dat", APPEND, TEXT_FILE)
+	    do i = 1, nprof {
+		call fprintf (j, "%d %g\n")
+		    call pargi (i)
+		    call pargr (prof[i])
+	    }
+	    call close (j)
+
+	    j = open ("debug2.dat", APPEND, TEXT_FILE)
+	    do i = 1, nprof1 {
+		call fprintf (j, "%d %g\n")
+		    call pargi (i)
+		    call pargr (prof1[i])
+	    }
+	    call close (j)
+	}
+end
+
+
+# MKT_GSTAR -- Get the precomputed template with center nearest (x,y).
+
+procedure mkt_gstar (mkt, data, nx, ny, x, y, z)
+
+pointer	mkt		# Template
+pointer	data		# Data
+int	nx, ny		# Size of data
+real	x, y		# Position of object
+real	z		# Flux of object
+
+int	i, j
+real	f
+
+include	"mktemplates.com"
+
+begin
+	i = (x + 0.5 - int (x + 0.5)) * nxc
+	j = (y + 0.5 - int (y + 0.5)) * nyc
+	i = j * nxc + i
+	data = Memi[MKT_PTRS(mkt)+i]
+	nx = Memi[MKT_NX(mkt)+i]
+	ny = Memi[MKT_NY(mkt)+i]
+	f = Memr[MKT_FLUX(mkt)+i]
+	if (f != z) {
+	    call amulkr (Memr[data], z/f, Memr[data], nx*ny)
+	    Memr[MKT_FLUX(mkt)+i] = z
+	}
+end
+
+
+# MKT_PROFILE -- Make template from profile.
+
+procedure mkt_profile (data, nx, ny, xc, yc, flux, prof, npts, radius,
+	ar, pa, nxsub, nysub)
+
+pointer	data		# Pointer to data array
+int	nx, ny		# Size of template
+real	xc, yc		# Model center
+real	flux		# Model flux
+real	prof[npts]	# 1D profile
+int	npts		# Number of points in profile
+real	radius		# Major axis radius of profile (pixels)
+real	ar		# Axial ratio (minor / major)
+real	pa		# Position angle relative to major axis (radians)
+int	nxsub, nysub	# Number of subpixel samples
+
+int	i, n, nxs, nys, nxs2, nys2, xs1, xs2, ys1, ys2
+int	x, y, x2, y2
+real	a, b, c, r, s, t, z, sum, sum1, asumr()
+real	dx, dy, dsub, dsub2
+real	x1, y1, xc1, yc1
+pointer	ptr, ptr1, ptr2, see1, see2
+
+include	"mktemplates.com"
+
+define	see_	99
+
+begin
+	# Switch on the size of the seeing templates.
+	if (see != NULL)
+	    if (Memi[MKT_NX(see)] * Memi[MKT_NY(see)] > 1)
+	        goto see_
+
+# NO SEEING:
+
+	# If the radius is very small return delta function.
+	if (radius < 0.1) {
+	    call aclrr (Memr[data], nx*ny)
+	    Memr[data+(ny/2)*nx+(nx/2)] = flux
+	    return
+	}
+
+	# Compute elliptical scale factors for entry into profile array.
+	r = ((npts - 1) / radius) ** 2
+	t = ((npts - 1) / (ar * radius)) ** 2
+	c = cos (pa)
+	s = sin (pa)
+	a = r * c * c + t * s * s
+	b = 2 * (r - t) * c * s
+	c = r * s * s + t * c * c
+
+	# Subsample the profile and sum into each pixel.
+	n = nx * ny
+	xc1 = xc - int (xc + 0.5) + nx/2
+	yc1 = yc - int (yc + 0.5) + ny/2
+
+	sum1 = 0.
+	ptr = data
+	do y = 0, ny-1 {
+	    t = (y - yc1) ** 2
+	    do x = 0, nx-1 {
+		i = sqrt ((x - xc1) ** 2 + t)
+		dsub = 1. / max (1, nxsub - i)
+		sum = 0.
+		for (y1 = y-0.5+dsub/2; y1 < y+0.5; y1 = y1+dsub) {
+		    dy = (y1 - yc1)
+		    s = c * dy**2
+		    for (x1=x-0.5+dsub/2; x1<x+0.5; x1=x1+dsub) {
+			dx = (x1-xc1)
+			i = sqrt (a * dx**2 + b * dx * dy + s) + 1.5
+			if (i <= npts)
+			    sum = sum + prof[i]
+		    }
+		}
+		sum = sum * dsub ** 2
+		sum1 = sum1 + sum
+		Memr[ptr] = sum
+		ptr = ptr + 1
+	    }
+	}
+
+	# If the subsamples completely miss and signal is zero then return
+	# delta function otherwise scale to requested flux.
+
+	if (sum1 == 0.)
+	    Memr[data+(ny/2)*nx+(nx/2)] = flux
+	else
+	    call amulkr (Memr[data], flux/sum1, Memr[data], n)
+	return
+
+# WITH SEEING:
+
+see_	n = nx * ny
+	call aclrr (Memr[data], n)
+	sum = 0.
+
+	nxs = Memi[MKT_NX(see)]
+	nys = Memi[MKT_NY(see)]
+	nxs2 = nxs/2
+	nys2 = nys/2
+
+	# If the profile is very small return full star image rather than
+	# convolution with truncated seeing template.
+
+	if (radius > 0.01) {
+	    r = ((npts - 1) / radius) ** 2
+	    t = ((npts - 1) / (ar * radius)) ** 2
+	    c = cos (pa)
+	    s = sin (pa)
+	    a = r * c * c + t * s * s
+	    b = 2 * (r - t) * c * s
+	    c = r * s * s + t * c * c
+    
+	    xc1 = xc - int (xc + 0.5) + nx/2
+	    yc1 = yc - int (yc + 0.5) + ny/2
+	    ptr = data-nys2*nx-nxs2
+	    do y = 0, ny-1 {
+		t = (y - yc1) ** 2
+		ys1 = max (0, nys2 - y)
+		ys2 = min (nys-1, ny - y + nys2 - 1)
+		do x = 0, nx-1 {
+		    i = sqrt ((x - xc1) ** 2 + t)
+		    dsub = 1. / max (1, nxsub - i)
+		    dsub2 = dsub ** 2
+		    xs1 = max (0, nxs2 - x)
+		    xs2 = min (nxs-1, nx - x + nxs2 - 1)
+		    ptr1 = ptr + xs1
+		    for (y1=y-0.5+dsub/2; y1<y+0.5; y1=y1+dsub) {
+			dy = (y1 - yc1)
+			s = c * dy**2
+			for (x1=x-0.5+dsub/2; x1<x+0.5; x1=x1+dsub) {
+			    dx = (x1-xc1)
+			    i = sqrt (a * dx**2 + b * dx * dy + s) + 1.5
+			    if (i <= npts) {
+				z = prof[i] * dsub2
+				call mkt_gstar (see, see1, nxs, nys, x1, y1, 1.)
+				see1 = see1 + xs1
+				do y2 = ys1, ys2 {
+				   see2 = see1+y2*nxs
+				   ptr2 = ptr1+y2*nx
+				   do x2 = xs1, xs2 {
+				       Memr[ptr2] = Memr[ptr2] +
+					   z * Memr[see2]
+				       ptr2 = ptr2 + 1
+				       see2 = see2 + 1
+				    }
+				}
+			    }
+			}
+		    }
+		    ptr = ptr + 1
+		}
+	    }
+	    sum = asumr (Memr[data], n)
+	}
+
+	# If no flux is accumulated or radius is very small return star image.
+	# Otherwise scale to requested flux.
+
+	if (sum == 0.) {
+	    call mkt_gstar (star, see1, nxs, nys, xc, yc, flux)
+	    x = nx/2
+	    y = ny/2
+	    nxs2 = nxs / 2
+	    nys2 = nys / 2
+	    xs1 = max (0, nxs2 - x)
+	    xs2 = min (nxs-1, nx - x + nxs2 - 1)
+	    ys1 = max (0, nys2 - y)
+	    ys2 = min (nys-1, ny - y + nys2 - 1)
+	    ptr1 = data-nys2*nx-nxs2+(y*nx+x+xs1)
+	    see1 = see1 + xs1
+	    do y2 = ys1, ys2 {
+	       see2 = see1+y2*nxs
+	       ptr2 = ptr1+y2*nx
+	       do x2 = xs1, xs2 {
+		   Memr[ptr2] = Memr[ptr2] + Memr[see2]
+		   ptr2 = ptr2 + 1
+		   see2 = see2 + 1
+		}
+	    }
+	} else
+	    call amulkr (Memr[data], flux/sum, Memr[data], n)
+end
+
+
+# MKT_MSI -- Make template from image interpolation function.
+
+procedure mkt_msi (data, nx, ny, xc, yc, flux, model, nxm, nym, radius,
+	ar, pa, nxsub, nysub)
+
+pointer	data		# Pointer to data array
+int	nx, ny		# Size of template
+real	xc, yc		# Model center
+real	flux		# Model flux
+pointer	model		# Surface interpolation pointer for image template
+int	nxm, nym	# Number of points in image template
+real	radius		# Major axis radius of profile (pixels)
+real	ar		# Axial ratio (minor / major)
+real	pa		# Position angle relative to major axis (radians)
+int	nxsub, nysub	# Number of subpixel samples
+
+int	i, n, nxs, nys, nxs2, nys2, xs1, xs2, ys1, ys2, x, y, x2, y2
+real	a, b, c, s, xcm, ycm, x1, y1, xc1, yc1, dsub, sum, sum1
+real	ac, as, bc, bs, acdx1, acdx2, bsdx1, bsdx2, bcdy1, bcdy2, asdy1, asdy2
+real	val1, val2, val3, val4, minval, maxval, xm[5], ym[5]
+real	asumr(), msigrl(), msisqgrl()
+pointer	ptr, ptr1, ptr2, see1, see2
+
+include	"mktemplates.com"
+
+define	see_	99
+
+begin
+	# Switch on the size of the seeing templates.
+	if (see != NULL)
+	    if (Memi[MKT_NX(see)] * Memi[MKT_NY(see)] > 1)
+	        goto see_
+
+# NO SEEING:
+
+	# If the radius is very small return delta function.
+	if (radius < 0.1) {
+	    call aclrr (Memr[data], nx*ny)
+	    Memr[data+(ny/2)*nx+(nx/2)] = flux
+	    return
+	}
+
+	a = (nxm / 2.) / radius
+	b = (nym / 2.) / (ar * radius)
+	c = cos (pa)
+	s = sin (pa)
+	ac = a * c
+	as = a * s
+	bc = b * c
+	bs = b * s
+	a = nxm
+	b = nym
+	xcm = nxm / 2 + 1.
+	ycm = nym / 2 + 1.
+
+	# Subsample the profile and sum into each pixel.
+	n = nx * ny
+	xc1 = xc - int (xc + 0.5) + nx/2
+	yc1 = yc - int (yc + 0.5) + ny/2
+
+	sum1 = 0.
+	ptr = data
+	do y = 0, ny-1 {
+	    c = (y - yc1) ** 2
+	    do x = 0, nx-1 {
+		i = sqrt ((x - xc1) ** 2 + c)
+		dsub = 1. / max (1, nxsub - i)
+		sum = 0.
+		for (y1 = y-0.5; y1 < y+0.5-dsub/2; y1 = y1+dsub) {
+		    asdy1 = (y1 - yc1)
+		    asdy2 = asdy1 + dsub
+		    bcdy1 = bc * asdy1 + ycm
+		    bcdy2 = bc * asdy2 + ycm
+		    if (pa == 0.) {
+			val3 = max (1., bcdy1)
+			if (val3 >= nym)
+			    next
+			val4 = min (b, bcdy2)
+			if (val3 >= val4)
+			    next
+		    }
+		    asdy1 = as * asdy1
+		    asdy2 = as * asdy2
+		    for (x1=x-0.5; x1<x+0.5-dsub/2; x1=x1+dsub) {
+			bsdx1 = (x1-xc1)
+			bsdx2 = bsdx1 + dsub
+			acdx1 = ac * bsdx1 + xcm
+			acdx2 = ac * bsdx2 + xcm
+			if (pa == 0.) {
+			    val1 = max (1., acdx1)
+			    if (val1 >= nxm)
+				next
+			    val2 = min (a, acdx2)
+			    if (val1 >= val2)
+				next
+			    s = msisqgrl (model, val1, val2, val3, val4)
+			} else {
+			    bsdx1 = bs * bsdx1
+			    bsdx2 = bs * bsdx2
+
+			    val1 = acdx1 + asdy1
+			    val2 = acdx2 + asdy1
+			    val3 = acdx2 + asdy2
+			    val4 = acdx1 + asdy2
+			    minval = min (val1, val2, val3, val4)
+			    maxval = max (val1, val2, val3, val4)
+			    if (minval >= a || maxval <= 1.)
+				next
+			    xm[1] = max (1., min (a, val1))
+			    xm[2] = max (1., min (a, val2))
+			    xm[3] = max (1., min (a, val3))
+			    xm[4] = max (1., min (a, val4))
+			    xm[5] = xm[1]
+
+			    val1 = bcdy1 - bsdx1
+			    val2 = bcdy1 - bsdx2
+			    val3 = bcdy2 - bsdx2
+			    val4 = bcdy2 - bsdx1
+			    minval = min (val1, val2, val3, val4)
+			    maxval = max (val1, val2, val3, val4)
+			    if (minval >= b || maxval <= 1.)
+				next
+			    ym[1] = max (1., min (b, val1))
+			    ym[2] = max (1., min (b, val2))
+			    ym[3] = max (1., min (b, val3))
+			    ym[4] = max (1., min (b, val4))
+			    ym[5] = ym[1]
+			    y2 = 1
+			    do x2 = 2, 4
+				if (ym[x2] < ym[y2])
+				    y2 = x2
+			    switch (y2) {
+			    case 2:
+				xm[1] = xm[2]; ym[1] = ym[2]
+				xm[2] = xm[3]; ym[2] = ym[3]
+				xm[3] = xm[4]; ym[3] = ym[4]
+				xm[4] = xm[5]; ym[4] = ym[5]
+				xm[5] = xm[1]; ym[5] = ym[1]
+			    case 3:
+				xm[1] = xm[3]; ym[1] = ym[3]
+				xm[3] = xm[5]; ym[3] = ym[5]
+				xm[5] = xm[2]; ym[5] = ym[2]
+				xm[2] = xm[4]; ym[2] = ym[4]
+				xm[4] = xm[5]; ym[4] = ym[5]
+				xm[5] = xm[1]; ym[5] = ym[1]
+			    case 4:
+				xm[5] = xm[4]; ym[5] = ym[4]
+				xm[4] = xm[3]; ym[4] = ym[3]
+				xm[3] = xm[2]; ym[3] = ym[2]
+				xm[2] = xm[1]; ym[2] = ym[1]
+				xm[1] = xm[5]; ym[1] = ym[5]
+			    }
+
+			    s = msigrl (model, xm, ym, 5)
+			}
+			sum = sum + s
+		    }
+		}
+		sum1 = sum1 + sum
+		Memr[ptr] = sum
+		ptr = ptr + 1
+	    }
+	}
+
+	call amulkr (Memr[data], flux/sum1, Memr[data], n)
+	return
+
+# SEEING:
+
+see_	n = nx * ny
+	call aclrr (Memr[data], n)
+	sum = 0.
+
+	nxs = Memi[MKT_NX(see)]
+	nys = Memi[MKT_NY(see)]
+	nxs2 = nxs/2
+	nys2 = nys/2
+
+	# If the profile is very small return full star image rather than
+	# convolution with truncated seeing template.
+
+	if (radius > 0.01) {
+	    a = (nxm / 2.) / radius
+	    b = (nym / 2.) / (ar * radius)
+	    c = cos (pa)
+	    s = sin (pa)
+	    ac = a * c
+	    as = a * s
+	    bc = b * c
+	    bs = b * s
+	    a = nxm
+	    b = nym
+	    xcm = nxm / 2 + 1.
+	    ycm = nym / 2 + 1.
+    
+	    xc1 = xc - int (xc + 0.5) + nx/2
+	    yc1 = yc - int (yc + 0.5) + ny/2
+	    ptr = data-nys2*nx-nxs2
+	    do y = 0, ny-1 {
+		c = (y - yc1) ** 2
+		ys1 = max (0, nys2 - y)
+		ys2 = min (nys-1, ny - y + nys2 - 1)
+		do x = 0, nx-1 {
+		    i = sqrt ((x - xc1) ** 2 + c)
+		    dsub = 1. / max (1, nxsub - i)
+		    xs1 = max (0, nxs2 - x)
+		    xs2 = min (nxs-1, nx - x + nxs2 - 1)
+		    ptr1 = ptr + xs1
+		    for (y1=y-0.5; y1<y+0.5-dsub/2; y1=y1+dsub) {
+			asdy1 = (y1 - yc1)
+			asdy2 = asdy1 + dsub
+			bcdy1 = bc * asdy1 + ycm
+			bcdy2 = bc * asdy2 + ycm
+		        if (pa == 0.) {
+			    val3 = max (1., bcdy1)
+			    if (val3 >= nym)
+			        next
+			    val4 = min (b, bcdy2)
+			    if (val3 >= val4)
+			        next
+		        }
+			asdy1 = as * asdy1
+			asdy2 = as * asdy2
+		        for (x1=x-0.5; x1<x+0.5-dsub/2; x1=x1+dsub) {
+			    bsdx1 = (x1-xc1)
+			    bsdx2 = bsdx1 + dsub
+			    acdx1 = ac * bsdx1 + xcm
+			    acdx2 = ac * bsdx2 + xcm
+			    if (pa == 0.) {
+			        val1 = max (1., acdx1)
+			        if (val1 >= nxm)
+				    next
+			        val2 = min (a, acdx2)
+			        if (val1 >= val2)
+				    next
+			        sum = msisqgrl (model, val1, val2, val3, val4)
+			    } else {
+			        bsdx1 = bs * bsdx1
+			        bsdx2 = bs * bsdx2
+
+			        val1 = acdx1 + asdy1
+			        val2 = acdx2 + asdy1
+			        val3 = acdx2 + asdy2
+			        val4 = acdx1 + asdy2
+			        minval = min (val1, val2, val3, val4)
+			        maxval = max (val1, val2, val3, val4)
+			        if (minval >= a || maxval <= 1.)
+				    next
+			        xm[1] = max (1., min (a, val1))
+			        xm[2] = max (1., min (a, val2))
+			        xm[3] = max (1., min (a, val3))
+			        xm[4] = max (1., min (a, val4))
+			        xm[5] = xm[1]
+
+			        val1 = bcdy1 - bsdx1
+			        val2 = bcdy1 - bsdx2
+			        val3 = bcdy2 - bsdx2
+			        val4 = bcdy2 - bsdx1
+			        minval = min (val1, val2, val3, val4)
+			        maxval = max (val1, val2, val3, val4)
+			        if (minval >= b || maxval <= 1.)
+				    next
+			        ym[1] = max (1., min (b, val1))
+			        ym[2] = max (1., min (b, val2))
+			        ym[3] = max (1., min (b, val3))
+			        ym[4] = max (1., min (b, val4))
+			        ym[5] = ym[1]
+
+# The following is put in to work around a bug in msigrl in V2.8.  When
+# V2.8 is long gone we can take this stuff out since msigrl will do the
+# rotating of the coordinates itself.
+
+				minval = max (1., minval)
+				y2 = 1
+				do x2 = 2, 4
+				    if (ym[x2] < ym[y2])
+					y2 = x2
+				switch (y2) {
+				case 2:
+				    xm[1] = xm[2]; ym[1] = ym[2]
+				    xm[2] = xm[3]; ym[2] = ym[3]
+				    xm[3] = xm[4]; ym[3] = ym[4]
+				    xm[4] = xm[5]; ym[4] = ym[5]
+				    xm[5] = xm[1]; ym[5] = ym[1]
+				case 3:
+				    xm[1] = xm[3]; ym[1] = ym[3]
+				    xm[3] = xm[5]; ym[3] = ym[5]
+				    xm[5] = xm[2]; ym[5] = ym[2]
+				    xm[2] = xm[4]; ym[2] = ym[4]
+				    xm[4] = xm[5]; ym[4] = ym[5]
+				    xm[5] = xm[1]; ym[5] = ym[1]
+				case 4:
+				    xm[5] = xm[4]; ym[5] = ym[4]
+				    xm[4] = xm[3]; ym[4] = ym[3]
+				    xm[3] = xm[2]; ym[3] = ym[2]
+				    xm[2] = xm[1]; ym[2] = ym[1]
+				    xm[1] = xm[5]; ym[1] = ym[5]
+				}
+
+			        sum = msigrl (model, xm, ym, 5)
+			    }
+	    		    call mkt_gstar (see, see1, nxs, nys, x1, y1, 1.)
+			    see1 = see1 + xs1
+			    do y2 = ys1, ys2 {
+				see2 = see1+y2*nxs
+				ptr2 = ptr1+y2*nx
+				do x2 = xs1, xs2 {
+				    Memr[ptr2] = Memr[ptr2] + sum * Memr[see2]
+				    ptr2 = ptr2 + 1
+				    see2 = see2 + 1
+				}
+			    }
+			}
+		    }
+		    ptr = ptr + 1
+		}
+	    }
+	    sum = asumr (Memr[data], n)
+	}
+
+	# If no flux is accumulated or radius is very small return star image.
+	# Otherwise scale to requested flux.
+
+	if (sum == 0.) {
+	    call mkt_gstar (star, see1, nxs, nys, xc, yc, flux)
+	    x = nx/2
+	    y = ny/2
+	    nxs2 = nxs / 2
+	    nys2 = nys / 2
+	    xs1 = max (0, nxs2 - x)
+	    xs2 = min (nxs-1, nx - x + nxs2 - 1)
+	    ys1 = max (0, nys2 - y)
+	    ys2 = min (nys-1, ny - y + nys2 - 1)
+	    ptr1 = data-nys2*nx-nxs2+(y*nx+x+xs1)
+	    see1 = see1 + xs1
+	    do y2 = ys1, ys2 {
+	       see2 = see1+y2*nxs
+	       ptr2 = ptr1+y2*nx
+	       do x2 = xs1, xs2 {
+		   Memr[ptr2] = Memr[ptr2] + Memr[see2]
+		   ptr2 = ptr2 + 1
+		   see2 = see2 + 1
+		}
+	    }
+	} else
+	    call amulkr (Memr[data], flux/sum, Memr[data], n)
+end
+
+
+# UIGAMMA -- Upper Incomplete Gamma Function
+#
+#	uigamma(n,x) = (n-1)!e^(x) sum (x^k/k!) for k=0, n-1
+
+real procedure uigamma (n, x)
+
+int	n		#I argument
+real	x		#I argument
+
+int	i
+real	uigamma
+double	sum, numerator, denominator
+
+begin
+	numerator = exp(-x)
+	denominator = 1
+	sum = numerator / denominator
+	do i = 1, n {
+	    numerator = numerator * x
+	    denominator = denominator * i
+	    sum = sum + numerator / denominator
+	}
+	uigamma = denominator * sum
+	return (uigamma)
+end