Repatch (from linux) of OSX IRAF

13 files changed, 4199 insertions, 0 deletions

diff --git a/noao/artdata/lists/gallist.key b/noao/artdata/lists/gallist.key
new file mode 100644
index 00000000..aa9d1822
--- /dev/null
+++ b/noao/artdata/lists/gallist.key
@@ -0,0 +1,63 @@
+	Gallist Keystroke Commands
+
+?	Print options
+f	Fit one or more of following 
+	    Spatial density function (SDF)
+            Luminosity  function (LF)
+	    Distribution of morphological type
+	    Diameter distribution
+	    Roundness distribution
+	    Position angle distribution 
+x	Plot the x-y spatial density function
+r	Plot the histogram of the radial density function
+m	Plot the histogram of the luminosity function
+d	Plot the histogram of the diameter values
+e	Plot the histogram of the roundness values 
+p	Plot the histogram of the position angle values
+:	Colon escape commands (see below)
+q	Exit program
+
+
+	Gallist Colon Commands
+
+:show			Show gallist parameters
+:ngal       [value]	Number of galaxies
+
+:spatial    [string]	Spatial density function (SDF) (uniform|hubble|file) 
+:xmin       [value]	Minimum X value
+:xmax       [value]	Maximum X value
+:ymin       [value]	Minimum Y value
+:ymax       [value]	Maximum Y value
+:xcenter    [value]	X center for radial density function
+:ycenter    [value]	Y center for radial density function
+:core       [value]	Core radius for Hubble density function
+:base       [value]	Background density for Hubble density function
+
+:luminosity [string]	Luminosity function (LF)
+			(uniform|powlaw|shcecter|file)
+:minmag     [value]	Minimum magnitude
+:maxmag     [value]	Maximum magnitude
+:power      [value]     Power law coefficient
+:mzero	    [value]     Zero point for Schecter luminosity function
+:alpha      [value]	Schecter parameter
+:mstar      [value]	Characteristic mag for Schecter function
+
+:egalmix    [value]	Elliptical/Spiral galaxy ratio
+:ar         [value]     Minimum elliptical galaxy axial ratio
+:eradius    [value]     Maximum elliptical half flux radius
+:sradius    [value]     Spiral/elliptical radius at same magnitude
+:absorption [value]	Absorption correction for edge-on spirals
+:z          [value]     Minimum redshift
+
+:lfile      [string]    Name of the luminosity function file
+:sfile	    [string]    Name of the spatial density function file
+:nlsample   [value]	Number of LF sample points 
+:lorder	    [value]	Order of spline approximation to the integrated LF
+:nssample   [value]	Number of SDF sample points
+:sorder	    [value]	Order of spline approximation to the integrated SDF
+
+:rbinsize   [value]	Resolution of radial SDF histogram in pixels
+:mbinsize   [value]	Resolution of magnitude histogram in magnitudes
+:dbinsize   [value]	Resolution of diameter histogram in pixels
+:ebinsize   [value]	Resolution of roundness histogram in pixels
+:pbinsize   [value]     Resolution of position angle histogram in degrees
diff --git a/noao/artdata/lists/mkpkg b/noao/artdata/lists/mkpkg
new file mode 100644
index 00000000..d8281097
--- /dev/null
+++ b/noao/artdata/lists/mkpkg
@@ -0,0 +1,18 @@
+# ARTDATA STARLIST and GALLIST task libraries
+
+$checkout libpkg.a ../
+$update libpkg.a
+$checkin libpkg.a ../
+$exit
+
+libpkg.a:
+	stcolon.x	starlist.h <error.h>
+	stdbio.x	starlist.h <pkg/dttext.h>
+	stlum.x		<mach.h> <math.h> <math/curfit.h> <math/iminterp.h>
+	stmix.x		<math.h> starlist.h
+	stplot.x	starlist.h <gset.h> <mach.h> <math.h> <pkg/gtools.h>
+	stshow.x	starlist.h
+	stspatial.x	<math.h> <math/curfit.h> <math/iminterp.h>
+	t_gallist.x	starlist.h <fset.h>
+	t_starlist.x	starlist.h <fset.h>
+	;
diff --git a/noao/artdata/lists/starlist.h b/noao/artdata/lists/starlist.h
new file mode 100644
index 00000000..01498222
--- /dev/null
+++ b/noao/artdata/lists/starlist.h
@@ -0,0 +1,136 @@
+# STARLIST/GALLIST task definitions file
+
+define	ST_STARS		1		# Make star list
+define	ST_GALAXIES		2		# Make galaxies list
+
+# Spatial distribution functions
+
+define	ST_UNIFORM		1		# Uniform spatial distribution
+define	ST_HUBBLE		2		# Hubble law
+define	ST_SPFILE		3		# User input
+
+# Luminosity distribution function
+
+define	ST_UNIFORM		1		# Uniform luminosity function
+define	ST_SALPETER		2		# Salpeter luminosity function
+define	ST_BANDS		3		# Bahcall and Soneira
+define	ST_LFFILE		4		# User input
+define	ST_POWLAW		5		# Power law
+define	ST_SCHECTER		6		# Schecter luminosity function
+
+# Galaxies types
+
+define	ST_DEVAUC		1		# Ellipticals
+define	ST_EXP			2		# Spirals
+
+define LEN_STSTRUCT		(45 + 4 * SZ_FNAME + 4)
+
+define	ST_TYPE			Memi[$1]	# Stars or galaxies
+define	ST_SPATIAL		Memi[$1+1]	# Spatial function
+define	ST_XC			Memr[P2R($1+2)]	# X center
+define	ST_YC			Memr[P2R($1+3)]	# Y center
+define	ST_CORE			Memr[P2R($1+4)]	# Hubble core radius
+define	ST_BASE			Memr[P2R($1+5)]	# Hubble baseline probability
+define	ST_XMIN			Memr[P2R($1+6)]	# Minimum x value
+define	ST_XMAX			Memr[P2R($1+7)]	# Maximum x value
+define	ST_YMIN			Memr[P2R($1+8)]	# Minimum y value
+define	ST_YMAX			Memr[P2R($1+9)]	# Maximum y value
+
+define	ST_LUMINOSITY		Memi[$1+10]	# Luminosity function
+define	ST_POWER		Memr[P2R($1+11)]# Power law
+define	ST_MZERO		Memr[P2R($1+12)]# Zero point of magnitudes
+define	ST_ALPHA		Memr[P2R($1+13)]# Bands function alpha
+define	ST_BETA			Memr[P2R($1+14)]# Bands function beta
+define	ST_DELTA		Memr[P2R($1+15)]# Bands function delta
+define	ST_MSTAR		Memr[P2R($1+16)]# Bands function mstar
+define	ST_MINMAG		Memr[P2R($1+17)]# Minimum magnitude
+define	ST_MAXMAG		Memr[P2R($1+18)]# Maximum magnitude
+
+define	ST_Z			Memr[P2R($1+19)]# Minimum redshift
+define	ST_AR			Memr[P2R($1+20)]# Minimum roundness
+define	ST_ERADIUS		Memr[P2R($1+21)]# Maximum elliptical radius
+define	ST_SRADIUS		Memr[P2R($1+22)]# Maximum spiral radius
+define	ST_EGALMIX		Memr[P2R($1+23)]# Egal fraction
+define	ST_ABSORPTION		Memr[P2R($1+24)]# Absorption
+
+define	ST_SSEED		Meml[$1+25]	# Spatial function seed
+define	ST_LSEED		Meml[$1+26]	# Luminosity function seed
+define	ST_NSSAMPLE		Memi[$1+27]	# Spatial function sampling
+define	ST_NLSAMPLE		Memi[$1+28]	# Luminosity function sampling
+define	ST_SORDER		Memi[$1+29]	# Spatial spline order
+define	ST_LORDER		Memi[$1+30]	# Luminosity spline order
+define	ST_NSTARS		Memi[$1+31]	# Number of stars
+
+define	ST_RBINSIZE		Memr[P2R($1+32)]# Radial histogram resolution
+define	ST_MBINSIZE		Memr[P2R($1+33)]# Magnitude histogram resolution
+define	ST_DBINSIZE		Memr[P2R($1+34)]# Diameter histogram resolution
+define	ST_EBINSIZE		Memr[P2R($1+35)]# Roundness histogram resolution
+define	ST_PBINSIZE		Memr[P2R($1+36)]# Posang histogram resolution
+
+define	ST_SPSTRING		Memc[P2C($1+37)]
+define	ST_LFSTRING		Memc[P2C($1+37+SZ_FNAME+1)]
+define	ST_SFILE		Memc[P2C($1+37+2*SZ_FNAME+2)]
+define	ST_LFILE		Memc[P2C($1+37+3*SZ_FNAME+3)]
+
+define	STCMDS "|show|nstars|spatial|xcenter|ycenter|core|base|xmin|xmax|\
+ymin|ymax|luminosity|power|alpha|beta|delta|mstar|minmag|maxmag|||nssample|\
+nlsample|sorder|lorder|sfile|lfile|rbinsize|mbinsize|ar|z|eradius|sradius|\
+egalmix|dbinsize|ebinsize|pbinsize|ngals|mzero|absorption|"
+
+define	SPFUNCS		"|uniform|hubble|file|"
+define	LUMFUNCS	"|uniform|salpeter|bands|file|powlaw|"
+define	GLUMFUNCS	"|uniform|||file|powlaw|schecter|"
+
+define	STCMD_SHOW		1
+define	STCMD_NSTARS		2
+define	STCMD_SPATIAL		3
+define	STCMD_XCENTER		4
+define	STCMD_YCENTER		5
+define	STCMD_CORE		6
+define	STCMD_BASE		7
+define	STCMD_XMIN		8
+define	STCMD_XMAX		9
+define	STCMD_YMIN		10
+define	STCMD_YMAX		11
+define	STCMD_LUMINOSITY	12
+define	STCMD_POWER		13
+define	STCMD_ALPHA		14
+define	STCMD_BETA		15
+define	STCMD_DELTA		16
+define	STCMD_MSTAR		17
+define	STCMD_MINMAG		18
+define	STCMD_MAXMAG		19
+define	STCMD_SSEED		20
+define	STCMD_LSEED		21
+define	STCMD_NSSAMPLE		22
+define	STCMD_NLSAMPLE		23
+define	STCMD_SORDER		24
+define	STCMD_LORDER		25
+define	STCMD_SFILE		26
+define	STCMD_LFILE		27
+define	STCMD_RBINSIZE		28
+define	STCMD_MBINSIZE		29
+define	STCMD_AR		30
+define	STCMD_Z			31
+define	STCMD_ERADIUS		32
+define	STCMD_SRADIUS		33
+define	STCMD_EGALMIX		34
+define	STCMD_DBINSIZE		35
+define	STCMD_EBINSIZE		36
+define	STCMD_PBINSIZE		37
+define	STCMD_NGALS		38
+define	STCMD_MZERO		39
+define	STCMD_ABSORPTION	40
+
+# Miscellaneous default values
+
+define	DEF_CORE		20.0
+define	DEF_BASE		0.00
+
+define	DEF_ALPHA		0.74
+define	DEF_BETA		0.04
+define	DEF_DELTA		0.294
+define	DEF_MSTAR		1.28
+
+define	DEF_GMSTAR		-20.6
+define	DEF_GALPHA		-1.25
diff --git a/noao/artdata/lists/starlist.key b/noao/artdata/lists/starlist.key
new file mode 100644
index 00000000..45f6c1fc
--- /dev/null
+++ b/noao/artdata/lists/starlist.key
@@ -0,0 +1,47 @@
+	Starlist Keystroke Commands
+
+?	Print options
+f	Fit  one or more of the following
+	    Spatial density function (SDF)
+	    Luminosity functions (LF)
+x	Plot the x-y spatial density function
+r	Plot the histogram of the radial density function
+m	Plot the histogram of the luminosity function
+:	Colon escape commands (see below)
+q	Exit program
+
+
+	Starlist Colon Commands
+
+:show			Show starlist parameters
+:nstars     [value]	Number of stars
+
+:spatial    [string]	Spatial density function (uniform|hubble|file) 
+:xmin       [value]	Minimum X value
+:xmax       [value]	Maximum X value
+:ymin       [value]	Minimum Y value
+:ymax       [value]	Maximum Y value
+:xcenter    [value]	X center for radial density function
+:ycenter    [value]	Y center for radial density function
+:core       [value]	Core radius for Hubble density function
+:base       [value]	Background density for Hubble density function
+
+:luminosity [string]	Luminosity function (uniform|powlaw|salpeter|bands|file)
+:minmag     [value]	Minimum magnitude
+:maxmag     [value]	Maximum magnitude
+:power	    [value]     Coefficent for powlaw magnitude distribution
+:mzero	    [value]	Magnitude zero point for salpeter and bands functions
+:alpha      [value]	Alpha parameter for bands luminosity function
+:beta       [value]	Beta parameter for bands luminosity function
+:delta      [value]	Delta parameter for bands luminosity function
+:mstar      [value]	Mstar parameter for bands luminosity function
+
+:sfile	    [string]    File containing the user spatial function
+:nssample   [value]	Number of SDF sample points
+:sorder	    [value]	Order of spline fit to spatial function
+:lfile	    [string]    File containing the user luminosity function
+:nlsample   [value]	Number of LF sample points 
+:lorder	    [value]	Order of spline fit to the integrated LF
+
+:rbinsize   [value]	Resolution of radial profile histogram (pixels)
+:mbinsize   [value]	Resolution of magnitude histogram (mag)
diff --git a/noao/artdata/lists/stcolon.x b/noao/artdata/lists/stcolon.x
new file mode 100644
index 00000000..9aadd578
--- /dev/null
+++ b/noao/artdata/lists/stcolon.x
@@ -0,0 +1,835 @@
+include <error.h>
+include "starlist.h"
+
+# ST_COLON -- Process colon commands for the STARLIST task.
+
+procedure st_colon (gd, st, sf, lf, cmdstr, newspace, newlum, newplot)
+
+pointer	gd		# pointer to the graphics stream
+pointer	st		# pointer to starlist structure
+int	sf		# spatial density function file descriptor
+int	lf		# luminosity function file descriptor
+char	cmdstr[ARB]	# input command string
+int	newspace	# new spatial distribution function
+int	newlum		# new luminosity function
+int	newplot		# new plot
+
+int	ival, ncmd, stat
+pointer	sp, cmd, str
+long	lval
+real	rval
+bool	streq()
+int	strdic(), nscan(), open()
+
+string	lumfuncs LUMFUNCS
+
+begin
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the command.
+	call sscan (cmdstr)
+	call gargwrd (Memc[cmd], SZ_LINE)
+	if (Memc[cmd] == EOS)
+	    return
+
+	# Process the command.
+	ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, STCMDS)
+	switch (ncmd) {
+
+	case STCMD_SHOW:
+	    call gdeactivate (gd, 0)
+	    call st_show (st)
+	    call greactivate (gd, 0)
+
+	case STCMD_SPATIAL:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (nscan () == 1) {
+	        call printf ("spatial = %s\n")
+		    call pargstr (ST_SPSTRING(st))
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, SPFUNCS)
+		if (stat > 0) {
+		    ST_SPATIAL(st) = stat
+		    call strcpy (Memc[cmd], ST_SPSTRING(st), SZ_FNAME)
+		}
+		newspace = YES
+	    }
+
+	case STCMD_XCENTER:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("xcenter = %g pixels\n")
+		    call pargr (ST_XC(st))
+	    } else {
+		ST_XC(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_YCENTER:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ycenter = %g pixels\n")
+		    call pargr (ST_YC(st))
+	    } else {
+		ST_YC(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_CORE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("core = %g pixels\n")
+		    call pargr (ST_CORE(st))
+	    } else {
+		ST_CORE(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_BASE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("base = %g pixels\n")
+		    call pargr (ST_BASE(st))
+	    } else {
+		ST_BASE(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_XMIN:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("xmin = %g pixels\n")
+		    call pargr (ST_XMIN(st))
+	    } else {
+		ST_XMIN(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_XMAX:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("xmax = %g pixels\n")
+		    call pargr (ST_XMAX(st))
+	    } else {
+		ST_XMAX(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_YMIN:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ymin = %g pixels\n")
+		    call pargr (ST_YMIN(st))
+	    } else {
+		ST_YMIN(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_YMAX:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ymax = %g pixels\n")
+		    call pargr (ST_YMAX(st))
+	    } else {
+		ST_YMAX(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_SFILE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], ST_SFILE(st))) {
+		call printf ("sfile: %s\n")
+		    call pargstr (ST_SFILE(st))
+	    } else {
+		if (sf != NULL) {
+		    call close (sf)
+		    sf = NULL
+		}
+		iferr {
+		    sf = open (Memc[cmd], READ_ONLY, TEXT_FILE)
+		} then {
+		    sf = NULL
+		    call erract (EA_WARN)
+		    call strcpy ("", ST_SFILE(st), SZ_FNAME)
+		    call printf (
+		        "Spatial distribution function file is undefined.\n")
+		} else {
+		    call strcpy (Memc[cmd], ST_SFILE(st), SZ_FNAME)
+		    newspace = YES
+		}
+	    }
+
+	case STCMD_LUMINOSITY:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (nscan () == 1) {
+	        call printf ("luminosity = %s\n")
+		    call pargstr (ST_LFSTRING(st))
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, LUMFUNCS)
+		if (stat > 0) {
+		    ST_LUMINOSITY (st) = stat
+		    call strcpy (Memc[cmd], ST_LFSTRING(st), SZ_FNAME)
+		}
+		newlum = YES
+	    }
+
+	case  STCMD_POWER:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("power = %g\n")
+		    call pargr (ST_POWER(st))
+	    } else {
+		ST_POWER(st) = rval
+		newlum = YES
+	    }
+
+	case  STCMD_MZERO:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("mzero = %g\n")
+		    call pargr (ST_MZERO(st))
+	    } else {
+		ST_MZERO(st) = rval
+		newlum = YES
+	    }
+
+	case  STCMD_ALPHA:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("alpha = %g\n")
+		    call pargr (ST_ALPHA(st))
+	    } else {
+		ST_ALPHA(st) = rval
+		newlum = YES
+	    }
+
+	case  STCMD_BETA:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+	        call printf ("beta = %g\n")
+		    call pargr (ST_BETA(st))
+	    } else {
+		ST_BETA(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_DELTA:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("delta = %g\n")
+		    call pargr (ST_DELTA(st))
+	    } else {
+		ST_DELTA(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_MSTAR:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("mstar = %g\n")
+		    call pargr (ST_MSTAR(st))
+	    } else {
+		ST_MSTAR(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_MINMAG:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("minmag = %g\n")
+		    call pargr (ST_MINMAG(st))
+	    } else {
+		ST_MINMAG(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_MAXMAG:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("maxmag = %g\n")
+		    call pargr (ST_MAXMAG(st))
+	    } else {
+		ST_MAXMAG(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_LFILE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], ST_LFILE(st))) {
+		call printf ("lfile: %s\n")
+		    call pargstr (ST_LFILE(st))
+	    } else {
+		if (lf != NULL) {
+		    call close (lf)
+		    lf = NULL
+		}
+		iferr {
+		    lf = open (Memc[cmd], READ_ONLY, TEXT_FILE)
+		} then {
+		    lf = NULL
+		    call erract (EA_WARN)
+		    call strcpy ("", ST_LFILE(st), SZ_FNAME)
+		    call printf (
+		        "Luminosity function file is undefined.\n")
+		} else {
+		    call strcpy (Memc[cmd], ST_LFILE(st), SZ_FNAME)
+		    newlum = YES
+		}
+	    }
+
+	case STCMD_NSTARS:
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("nstars = %d\n")
+		    call pargi (ST_NSTARS(st))
+	    } else {
+		ST_NSTARS(st) = ival
+		newlum = YES
+		newspace = YES
+	    }
+
+	case STCMD_NSSAMPLE:
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("nssample = %d\n")
+		    call pargi (ST_NSSAMPLE(st))
+	    } else {
+		ST_NSSAMPLE(st) = ival
+		newspace = YES
+	    }
+
+	case STCMD_SORDER:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("sorder = %d\n")
+		    call pargl (ST_SORDER(st))
+	    } else {
+		ST_SORDER(st) = lval
+		newspace = YES
+	    }
+
+	case STCMD_SSEED:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("sseed = %d\n")
+		    call pargl (ST_SSEED(st))
+	    } else {
+		ST_SSEED(st) = lval
+		newspace = YES
+	    }
+
+	case STCMD_NLSAMPLE:
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("nlsample = %d\n")
+		    call pargi (ST_NLSAMPLE(st))
+	    } else {
+		ST_NLSAMPLE(st) = ival
+		newlum = YES
+	    }
+
+	case STCMD_LORDER:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("lorder = %d\n")
+		    call pargl (ST_LORDER(st))
+	    } else {
+		ST_LORDER(st) = lval
+		newlum = YES
+	    }
+
+	case STCMD_LSEED:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("lseed = %d\n")
+		    call pargl (ST_LSEED(st))
+	    } else {
+		ST_LSEED(st) = lval
+		newlum = YES
+	    }
+
+	case STCMD_RBINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("rbinsize = %g\n")
+		    call pargr (ST_RBINSIZE(st))
+	    } else {
+		ST_RBINSIZE(st) = rval
+		newplot = YES
+	    }
+
+	case STCMD_MBINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("mbinsize = %g\n")
+		    call pargr (ST_MBINSIZE(st))
+	    } else {
+		ST_MBINSIZE(st) = rval
+		newplot = YES
+	    }
+
+	default:
+	    call printf ("\7\n")
+	}
+
+	call sfree (sp)
+end
+
+
+# ST_GCOLON -- Process colon commands for the GALAXIES task.
+
+procedure st_gcolon (gd, st, sf, lf, cmdstr, newspace, newlum, newmix, newaxis,
+	newround, newphi, newplot)
+
+pointer	gd		# pointer to the graphics stream
+pointer	st		# pointer to starlist structure
+int	sf		# spatial distribution function file descriptor
+int	lf		# luminosity function file descriptor
+char	cmdstr[ARB]	# command string
+int	newspace	# new spatial distribution function
+int	newlum		# new luminosity function
+int	newmix		# new E/S galaxy mixture
+int	newaxis		# new axis parameter
+int	newround	# compute new roundness parameters
+int	newphi		# compute new position angles
+int	newplot		# make a newplot
+
+int	ival, ncmd, stat
+pointer	sp, cmd, str
+long	lval
+real	rval
+bool	streq()
+int	strdic(), nscan(), open()
+
+string	lumfuncs LUMFUNCS
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the command.
+	call sscan (cmdstr)
+	call gargwrd (Memc[cmd], SZ_LINE)
+	if (Memc[cmd] == EOS)
+	    return
+
+	# Process the command.
+	ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, STCMDS)
+	switch (ncmd) {
+	case STCMD_SHOW:
+	    call gdeactivate (gd, 0)
+	    call st_gshow (st)
+	    call greactivate (gd, 0)
+
+	case STCMD_SPATIAL:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (nscan () == 1) {
+	        call printf ("spatial = %s\n")
+		    call pargstr (ST_SPSTRING(st))
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, SPFUNCS)
+		if (stat > 0) {
+		    ST_SPATIAL(st) = stat
+		    call strcpy (Memc[cmd], ST_SPSTRING(st), SZ_FNAME)
+		}
+		newspace = YES
+	    }
+
+	case STCMD_XCENTER:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("xcenter = %g pixels\n")
+		    call pargr (ST_XC(st))
+	    } else {
+		ST_XC(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_YCENTER:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ycenter = %g pixels\n")
+		    call pargr (ST_YC(st))
+	    } else {
+		ST_YC(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_CORE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("core = %g pixels\n")
+		    call pargr (ST_CORE(st))
+	    } else {
+		ST_CORE(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_BASE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("base = %g pixels\n")
+		    call pargr (ST_BASE(st))
+	    } else {
+		ST_BASE(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_XMIN:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("xmin = %g pixels\n")
+		    call pargr (ST_XMIN(st))
+	    } else {
+		ST_XMIN(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_XMAX:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("xmax = %g pixels\n")
+		    call pargr (ST_XMAX(st))
+	    } else {
+		ST_XMAX(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_YMIN:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ymin = %g pixels\n")
+		    call pargr (ST_YMIN(st))
+	    } else {
+		ST_YMIN(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_YMAX:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ymax = %g pixels\n")
+		    call pargr (ST_YMAX(st))
+	    } else {
+		ST_YMAX(st) = rval
+		newspace = YES
+	    }
+
+	case STCMD_SFILE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], ST_SFILE(st))) {
+		call printf ("sfile: %s\n")
+		    call pargstr (ST_SFILE(st))
+	    } else {
+		if (sf != NULL) {
+		    call close (sf)
+		    sf = NULL
+		}
+		iferr {
+		    sf = open (Memc[cmd], READ_ONLY, TEXT_FILE)
+		} then {
+		    sf = NULL
+		    call erract (EA_WARN)
+		    call strcpy ("", ST_SFILE(st), SZ_FNAME)
+		    call printf (
+		        "Spatial distribution function file is undefined.\n")
+		} else {
+		    call strcpy (Memc[cmd], ST_SFILE(st), SZ_FNAME)
+		    newspace = YES
+		}
+	    }
+
+	case STCMD_LUMINOSITY:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (nscan () == 1) {
+	        call printf ("luminosity = %s\n")
+		    call pargstr (ST_LFSTRING(st))
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, GLUMFUNCS)
+		if (stat > 0) {
+		    ST_LUMINOSITY (st) = stat
+		    call strcpy (Memc[cmd], ST_LFSTRING(st), SZ_FNAME)
+		}
+		newlum = YES
+	    }
+
+	case  STCMD_POWER:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("power = %g\n")
+		    call pargr (ST_POWER(st))
+	    } else {
+		ST_POWER(st) = rval
+		newlum = YES
+	    }
+
+	case  STCMD_MZERO:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("mzero = %g\n")
+		    call pargr (ST_MZERO(st))
+	    } else {
+		ST_MZERO(st) = rval
+		newlum = YES
+	    }
+
+	case  STCMD_ALPHA:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("alpha = %g\n")
+		    call pargr (ST_ALPHA(st))
+	    } else {
+		ST_ALPHA(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_MSTAR:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("mstar = %g\n")
+		    call pargr (ST_MSTAR(st))
+	    } else {
+		ST_MSTAR(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_MINMAG:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("minmag = %g\n")
+		    call pargr (ST_MINMAG(st))
+	    } else {
+		ST_MINMAG(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_MAXMAG:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("maxmag = %g\n")
+		    call pargr (ST_MAXMAG(st))
+	    } else {
+		ST_MAXMAG(st) = rval
+		newlum = YES
+	    }
+
+	case STCMD_LFILE:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS || streq (Memc[cmd], ST_LFILE(st))) {
+		call printf ("lfile: %s\n")
+		    call pargstr (ST_LFILE(st))
+	    } else {
+		if (lf != NULL) {
+		    call close (lf)
+		    lf = NULL
+		}
+		iferr {
+		    lf = open (Memc[cmd], READ_ONLY, TEXT_FILE)
+		} then {
+		    lf = NULL
+		    call erract (EA_WARN)
+		    call strcpy ("", ST_LFILE(st), SZ_FNAME)
+		    call printf (
+		        "Luminosity function file is undefined.\n")
+		} else {
+		    call strcpy (Memc[cmd], ST_LFILE(st), SZ_FNAME)
+		    newlum = YES
+		}
+	    }
+
+
+	case STCMD_EGALMIX:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("egalmix = %g\n")
+		    call pargr (ST_EGALMIX(st))
+	    } else {
+		ST_EGALMIX(st) = rval
+		newmix = YES
+	    }
+
+	case STCMD_ERADIUS:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("eradius = %g\n")
+		    call pargr (ST_ERADIUS(st))
+	    } else {
+		ST_ERADIUS(st) = rval
+		newaxis = YES
+	    }
+
+	case STCMD_SRADIUS:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("sradius = %g\n")
+		    call pargr (ST_SRADIUS(st))
+	    } else {
+		ST_SRADIUS(st) = rval
+		newaxis = YES
+	    }
+
+	case STCMD_AR:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ar = %g\n")
+		    call pargr (ST_AR(st))
+	    } else {
+		ST_AR(st) = rval
+		newround = YES
+	    }
+
+	case STCMD_Z:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("z = %g\n")
+		    call pargr (ST_Z(st))
+	    } else {
+		ST_Z(st) = rval
+		newaxis = YES
+	    }
+
+	case STCMD_ABSORPTION:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("absorption = %g\n")
+		    call pargr (ST_ABSORPTION(st))
+	    } else {
+		ST_ABSORPTION(st) = rval
+		newround = YES
+	    }
+
+	case STCMD_NGALS:
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("ngals = %d\n")
+		    call pargi (ST_NSTARS(st))
+	    } else {
+		ST_NSTARS(st) = ival
+		newlum = YES
+		newspace = YES
+		newmix = YES
+		newaxis = YES
+		newround = YES
+		newphi = YES
+	    }
+
+	case STCMD_NSSAMPLE:
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("nssample = %d\n")
+		    call pargi (ST_NSSAMPLE(st))
+	    } else {
+		ST_NSSAMPLE(st) = ival
+		newspace = YES
+	    }
+
+	case STCMD_SORDER:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("sorder = %d\n")
+		    call pargl (ST_SORDER(st))
+	    } else {
+		ST_SORDER(st) = lval
+		newspace = YES
+	    }
+
+	case STCMD_SSEED:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("sseed = %d\n")
+		    call pargl (ST_SSEED(st))
+	    } else {
+		ST_SSEED(st) = lval
+		newspace = YES
+	    }
+
+	case STCMD_NLSAMPLE:
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("nlsample = %d\n")
+		    call pargi (ST_NLSAMPLE(st))
+	    } else {
+		ST_NLSAMPLE(st) = ival
+		newlum = YES
+	    }
+
+	case STCMD_LORDER:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("lorder = %d\n")
+		    call pargl (ST_LORDER(st))
+	    } else {
+		ST_LORDER(st) = lval
+		newlum = YES
+	    }
+
+	case STCMD_LSEED:
+	    call gargl (lval)
+	    if (nscan() == 1) {
+		call printf ("lseed = %d\n")
+		    call pargl (ST_LSEED(st))
+	    } else {
+		ST_LSEED(st) = lval
+		newlum = YES
+	    }
+
+	case STCMD_RBINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("rbinsize = %g\n")
+		    call pargr (ST_RBINSIZE(st))
+	    } else {
+		ST_RBINSIZE(st) = rval
+		newplot = YES
+	    }
+
+	case STCMD_MBINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("mbinsize = %g\n")
+		    call pargr (ST_MBINSIZE(st))
+	    } else {
+		ST_MBINSIZE(st) = rval
+		newplot = YES
+	    }
+
+	case STCMD_DBINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("dbinsize = %g\n")
+		    call pargr (ST_DBINSIZE(st))
+	    } else {
+		ST_DBINSIZE(st) = rval
+		newplot = YES
+	    }
+
+	case STCMD_EBINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("ebinsize = %g\n")
+		    call pargr (ST_EBINSIZE(st))
+	    } else {
+		ST_EBINSIZE(st) = rval
+		newplot = YES
+	    }
+
+	case STCMD_PBINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("pbinsize = %g\n")
+		    call pargr (ST_PBINSIZE(st))
+	    } else {
+		ST_PBINSIZE(st) = rval
+		newplot = YES
+	    }
+
+	default:
+	    call printf ("\7\n")
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/artdata/lists/stdbio.x b/noao/artdata/lists/stdbio.x
new file mode 100644
index 00000000..1ec411b8
--- /dev/null
+++ b/noao/artdata/lists/stdbio.x
@@ -0,0 +1,350 @@
+include <pkg/dttext.h>
+include "starlist.h"
+
+# ST_DTINIT -- Write the header to the database.
+
+procedure st_dtinit (dt, st, starlist, sseed, lseed)
+
+pointer	dt				# pointer to output database
+pointer	st				# pointer to structure
+char	starlist[ARB]			# name of output text file
+long	sseed				# spatial function seed
+long	lseed				# luminsosity function seed
+
+begin
+	call dtptime (dt)
+	call dtput (dt, "# begin\t%s\n")
+	    call pargstr (starlist)
+
+	# Write out the spatial density function parameters.
+	call dtput (dt, "#\tspatial\t\t%s\n")
+	switch (ST_SPATIAL(st)) {
+	case ST_UNIFORM:
+	    call pargstr ("uniform")
+	case ST_HUBBLE:
+	    call pargstr ("hubble")
+	    call dtput (dt, "#\txcenter\t\t%g\n")
+		call pargr (ST_XC(st))
+	    call dtput (dt, "#\tycenter\t\t%g\n")
+		call pargr (ST_YC(st))
+	    call dtput (dt, "#\tcoreradius\t%g\n")
+		call pargr (ST_CORE(st))
+	    call dtput (dt, "#\tbaseline\t%g\n")
+		call pargr (ST_BASE(st))
+	case ST_SPFILE:
+	    call pargstr (ST_SFILE(st))
+	    call dtput (dt, "#\txcenter\t\t%g\n")
+		call pargr (ST_XC(st))
+	    call dtput (dt, "#\tycenter\t\t%g\n")
+		call pargr (ST_YC(st))
+	}
+	call dtput (dt, "#\txmin\t\t%g\n")
+	    call pargr (ST_XMIN(st))
+	call dtput (dt, "#\txmax\t\t%g\n")
+	    call pargr (ST_XMAX(st))
+	call dtput (dt, "#\tymin\t\t%g\n")
+	    call pargr (ST_YMIN(st))
+	call dtput (dt, "#\tymax\t\t%g\n")
+	    call pargr (ST_YMAX(st))
+
+	# Write out the luminsosity function parameters.
+	call dtput (dt, "#\tluminosity\t%s\n")
+	switch (ST_LUMINOSITY(st)) {
+	case ST_UNIFORM:
+	    call pargstr ("uniform")
+	case ST_POWLAW:
+	    call pargstr ("powlaw")
+	    call dtput (dt, "#\tpower\t\t%g\n")
+	        call pargr (ST_POWER(st))
+	case ST_SALPETER:
+	    call pargstr ("salpeter")
+	    call dtput (dt, "#\tmzero\t\t%g\n")
+	        call pargr (ST_MZERO(st))
+	case ST_BANDS:
+	    call pargstr ("bands")
+	    call dtput (dt, "#\tmzero\t\t%g\n")
+	        call pargr (ST_MZERO(st))
+	    call dtput (dt, "#\talpha\t\t%g\n")
+	        call pargr (ST_ALPHA(st))
+	    call dtput (dt, "#\tbeta\t\t%g\n")
+	        call pargr (ST_BETA(st))
+	    call dtput (dt, "#\tdelta\t\t%g\n")
+	        call pargr (ST_DELTA(st))
+	    call dtput (dt, "#\tmstar\t\t%g\n")
+	        call pargr (ST_MSTAR(st))
+	case ST_LFFILE:
+	    call pargstr (ST_LFILE(st))
+	}
+	call dtput (dt, "#\tminmag\t\t%g\n")
+	    call pargr (ST_MINMAG(st))
+	call dtput (dt, "#\tmaxmag\t\t%g\n")
+	    call pargr (ST_MAXMAG(st))
+
+	# Save the spatial density function fitting parameters.
+	call dtput (dt, "#\tnssample\t%d\n")
+	    call pargi (ST_NSSAMPLE(st))
+	call dtput (dt, "#\tsorder\t\t%d\n")
+	    call pargi (ST_SORDER(st))
+	call dtput (dt, "#\tsseed\t\t%d\n")
+	    call pargl (sseed)
+
+	# Save the luminosity function fitting parameters.
+	call dtput (dt, "#\tnlsample\t%d\n")
+	    call pargi (ST_NLSAMPLE(st))
+	call dtput (dt, "#\tlorder\t\t%d\n")
+	    call pargi (ST_LORDER(st))
+	call dtput (dt, "#\tlseed\t\t%d\n")
+	    call pargl (lseed)
+
+	# Save the number of stars.
+	call dtput (dt, "#\tnstars\t\t%d\n")
+	    call pargi (ST_NSTARS(st))
+end
+
+
+# ST_DTGINIT -- Write the GALLIST header to the database.
+
+procedure st_dtginit (dt, st, galaxies, sseed, lseed)
+
+pointer	dt				# pointer to database
+pointer	st				# pointer to starlist structure
+char	galaxies[ARB]			# name of output text file
+long	sseed				# spatial function seed
+long	lseed				# luminsosity function seed
+
+begin
+	call dtptime (dt)
+	call dtput (dt, "# begin\t%s\n")
+	    call pargstr (galaxies)
+
+	# Save the spatial distribution function parameters.
+	call dtput (dt, "#\tspatial\t\t%s\n")
+	switch (ST_SPATIAL(st)) {
+	case ST_UNIFORM:
+	    call pargstr ("uniform")
+	case ST_HUBBLE:
+	    call pargstr ("hubble")
+	    call dtput (dt, "#\txcenter\t\t%g\n")
+		call pargr (ST_XC(st))
+	    call dtput (dt, "#\tycenter\t\t%g\n")
+		call pargr (ST_YC(st))
+	    call dtput (dt, "#\tcoreradius\t%g\n")
+		call pargr (ST_CORE(st))
+	    call dtput (dt, "#\tbaseline\t%g\n")
+		call pargr (ST_BASE(st))
+	case ST_SPFILE:
+	    call pargstr (ST_SFILE(st))
+	    call dtput (dt, "#\txcenter\t\t%g\n")
+		call pargr (ST_XC(st))
+	    call dtput (dt, "#\tycenter\t\t%g\n")
+		call pargr (ST_YC(st))
+	}
+	call dtput (dt, "#\txmin\t\t%g\n")
+	    call pargr (ST_XMIN(st))
+	call dtput (dt, "#\txmax\t\t%g\n")
+	    call pargr (ST_XMAX(st))
+	call dtput (dt, "#\tymin\t\t%g\n")
+	    call pargr (ST_YMIN(st))
+	call dtput (dt, "#\tymax\t\t%g\n")
+	    call pargr (ST_YMAX(st))
+
+	# Save the luminsosity function parameters.
+	call dtput (dt, "#\tluminosity\t%s\n")
+	switch (ST_LUMINOSITY(st)) {
+	case ST_UNIFORM:
+	    call pargstr ("uniform")
+	case ST_POWLAW:
+	    call pargstr ("powlaw")
+	    call dtput (dt, "#\tpower\t\t%g\n")
+		call pargr (ST_POWER(st))
+	case ST_SCHECTER:
+	    call pargstr ("shechter")
+	    call dtput (dt, "#\tmzero\t\t%g\n")
+	        call pargr (ST_MZERO(st))
+	    call dtput (dt, "#\talpha\t\t%g\n")
+		call pargr (ST_ALPHA(st))
+	    call dtput (dt, "#\tmstar\t\t%g\n")
+		call pargr (ST_MSTAR(st))
+	case ST_LFFILE:
+	    call pargstr (ST_LFILE(st))
+	}
+	call dtput (dt, "#\tminmag\t\t%g\n")
+	    call pargr (ST_MINMAG(st))
+	call dtput (dt, "#\tmaxmag\t\t%g\n")
+	    call pargr (ST_MAXMAG(st))
+	call dtput (dt, "#\teradius\t\t%g\n")
+	    call pargr (ST_ERADIUS(st))
+	call dtput (dt, "#\tsradius\t\t%g\n")
+	    call pargr (ST_SRADIUS(st))
+
+	call dtput (dt, "#\tegalmix\t\t%g\n")
+	    call pargr (ST_EGALMIX(st))
+	call dtput (dt, "#\tar\t\t%g\n")
+	    call pargr (ST_AR(st))
+	call dtput (dt, "#\tabsorption\t%g\n")
+	    call pargr (ST_ABSORPTION(st))
+	call dtput (dt, "#\tz\t\t%g\n")
+	    call pargr (ST_Z(st))
+
+	# Save the spatial distribution fitting parameters.
+	call dtput (dt, "#\tnssample\t%d\n")
+	    call pargi (ST_NSSAMPLE(st))
+	call dtput (dt, "#\tsorder\t\t%d\n")
+	    call pargi (ST_SORDER(st))
+	call dtput (dt, "#\tsseed\t\t%d\n")
+	    call pargl (sseed)
+
+	# Save the spatial function fitting parameters.
+	call dtput (dt, "#\tnlsample\t%d\n")
+	    call pargi (ST_NLSAMPLE(st))
+	call dtput (dt, "#\tlorder\t\t%d\n")
+	    call pargi (ST_LORDER(st))
+	call dtput (dt, "#\tlseed\t\t%d\n")
+	    call pargl (lseed)
+
+	# Save the number of stars.
+	call dtput (dt, "#\tngals\t\t%d\n")
+	    call pargi (ST_NSTARS(st))
+end
+
+
+# ST_DTWRITE -- Write the starlist to the database.
+
+procedure st_dtwrite (dt, x, y, mag, nstars)
+
+pointer	dt		# pointer to the output database
+real	x[ARB]		# array of x coordinates
+real	y[ARB]		# array of y coordinates
+real	mag[ARB]	# array of magnitude values
+int	nstars		# number of stars
+
+int	i, j
+pointer	sp, index
+
+begin
+	call smark (sp)
+	call salloc (index, nstars, TY_INT)
+	call st_qsort (y, Memi[index], Memi[index], nstars)
+
+	do i = 1, nstars {
+	    j = Memi[index+i-1]
+	    call dtput (dt, "\t%8.3f  %8.3f  %7.3f\n")
+		call pargr (x[j])
+		call pargr (y[j])
+		call pargr (mag[j])
+	}
+
+	call sfree (sp)
+end
+
+
+# ST_DTGWRITE -- Procedure to write the galaxy list to the database
+
+procedure st_dtgwrite (dt, x, y, mag, egal, axis, round, phi, nstars)
+
+pointer	dt		# pointer to database
+real	x[ARB]		# x values
+real	y[ARB]		# y values
+real	mag[ARB]	# magnitude values
+int	egal[ARB]	# galaxy types
+real	axis[ARB]	# galaxy diameters
+real	round[ARB]	# galaxy roundness
+real	phi[ARB]	# galaxy position angles
+int	nstars		# number of stars
+
+int	i, j
+pointer	sp, index
+
+begin
+	call smark (sp)
+	call salloc (index, nstars, TY_INT)
+	call st_qsort (y, Memi[index], Memi[index], nstars)
+
+	do i = 1, nstars {
+	    j = Memi[index+i-1]
+	    call dtput (dt,
+	        "\t%8.3f  %8.3f  %7.3f  %7s  %7.2f  %5.3f  %5.1f\n")
+		call pargr (x[j])
+		call pargr (y[j])
+		call pargr (mag[j])
+		if (egal[j] == ST_DEVAUC)
+		    call pargstr ("devauc")
+		else
+		    call pargstr ("expdisk")
+		call pargr (axis[j])
+		call pargr (round[j])
+		call pargr (phi[j])
+	}
+
+	call sfree (sp)
+end
+
+
+# ST_QSORT -- Vector Quicksort. In this version the index array is
+# sorted.
+
+define	LOGPTR		20			# log2(maxpts) (1e6)
+
+procedure st_qsort (data, a, b, npix)
+
+real	data[ARB]		# data array
+int	a[ARB], b[ARB]		# index array
+int	npix			# number of pixels
+
+int	i, j, lv[LOGPTR], p, uv[LOGPTR], temp
+real	pivot
+
+begin
+	# Initialize the indices for an inplace sort.
+	do i = 1, npix
+	    a[i] = i
+	call amovi (a, b, npix)
+
+	# Initialize.
+	p = 1
+	lv[1] = 1
+	uv[1] = npix
+
+	# Sort.
+	while (p > 0) {
+
+	    # If only one elem in subset pop stack otherwise pivot line.
+	    if (lv[p] >= uv[p])
+		p = p - 1
+	    else {
+		i = lv[p] - 1
+		j = uv[p]
+		pivot = data[b[j]]
+
+		while (i < j) {
+		    for (i=i+1;  data[b[i]] < pivot;  i=i+1)
+			;
+		    for (j=j-1;  j > i;  j=j-1)
+			if (data[b[j]] <= pivot)
+			    break
+		    if (i < j) {		# Out of order pair
+			temp = b[j]		# Interchange elements
+			b[j] = b[i]
+			b[i] = temp
+		    }
+		}
+
+		j = uv[p]			# Move pivot to position i
+		temp = b[j]			# Interchange elements
+		b[j] = b[i]
+		b[i] = temp
+
+		if (i-lv[p] < uv[p] - i) {	# Stack so shorter done first
+		    lv[p+1] = lv[p]
+		    uv[p+1] = i - 1
+		    lv[p] = i + 1
+		} else {
+		    lv[p+1] = i + 1
+		    uv[p+1] = uv[p]
+		    uv[p] = i - 1
+		}
+
+		p = p + 1			# Push onto stack
+	    }
+	}
+end
diff --git a/noao/artdata/lists/stlum.x b/noao/artdata/lists/stlum.x
new file mode 100644
index 00000000..0abe9e1e
--- /dev/null
+++ b/noao/artdata/lists/stlum.x
@@ -0,0 +1,342 @@
+include <mach.h>
+include <math.h>
+include <math/curfit.h>
+include <math/iminterp.h>
+
+
+# ST_MAGUNIFORM -- Compute a set of magnitude values which are uniformly
+# distributed between minmag and maxmag.
+
+procedure st_maguniform (mag, nstars, minmag, maxmag, seed)
+
+real	mag[ARB]	# output array of magnitudes
+int	nstars		# number of stars
+real	minmag, maxmag	# minimum and maximum magnitude values
+long	seed		# seed for random number generator
+
+int	i
+real	urand()
+
+begin
+	# Get values between 0 and 1.
+	do i = 1, nstars
+	    mag[i] = urand (seed)
+
+	# Map values into data range.
+	call amapr (mag, mag, nstars, 0.0, 1.0, minmag, maxmag)
+end
+
+
+# ST_POWER -- Compute a set of magnitude values which are power law
+# distributed between minmag and maxmag.
+
+procedure st_power (mag, nstars, power, minmag, maxmag, seed)
+
+real	mag[ARB]	# output array of magnitudes
+int	nstars		# number of stars
+real	power		# power law exponent
+real	minmag, maxmag	# minimum and maximum magnitude values
+long	seed		# seed for random number generator
+
+int	i
+real	a, urand()
+
+begin
+	# Get values between 0 and 1.
+	a = 10. ** (power * (maxmag - minmag)) - 1
+	do i = 1, nstars
+	    mag[i] = minmag + log10 (a * urand (seed) + 1) / power
+end
+
+
+define	MIN_BANDS	-6.
+define	MAX_BANDS	19.
+define	MID_BANDS	15.
+
+# ST_BANDS -- Compute the Bahcall and Soneira luminosity function.
+
+procedure st_bands (mag, nstars, alpha, beta, delta, mstar, minmag, maxmag,
+	mzero, nsample, order, seed)
+
+real	mag[ARB]		# array of output magnitudes
+int	nstars			# number of stars
+real	alpha, beta		# Bahcall and Soneira parameters
+real	delta, mstar		# Bahcall and Soneira parameters
+real	minmag, maxmag		# minimum and maximum magnitude values
+real	mzero			# zero point between relative and absolute mags.
+int	nsample			# number of points in sampling function
+int	order			# order of the spline fit
+long	seed			# value of the seed
+
+int	i, ier
+pointer	sp, m, iprob, w, cv, asi
+real	dmag, magval, mtemp, temp1, temp2, imin, imax
+real	cveval(), asigrl(), urand()
+
+begin
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (m, nsample, TY_REAL)
+	call salloc (iprob, nsample, TY_REAL)
+	call salloc (w, nsample, TY_REAL)
+
+	# Compute the probability function.
+	magval = max (minmag - mzero, MIN_BANDS)
+	dmag = (min (maxmag - mzero, MAX_BANDS) - magval) / (nsample - 1)
+	do i = 1, nsample {
+	    Memr[m+i-1] = magval
+	    if (magval > MID_BANDS)
+		mtemp = MID_BANDS
+	    else
+	        mtemp = magval - mstar
+	    temp1 = 10.0 ** (beta * mtemp)
+	    temp2 = (1.0 + 10.0 ** ((beta - alpha) * delta * mtemp)) **
+	            (1.0 / delta)
+	    Memr[iprob+i-1] = temp1 / temp2
+	    magval = magval + dmag
+	}
+
+	# Integrate the probablity function.
+	call asiinit (asi, II_SPLINE3)
+	call asifit (asi, Memr[iprob], nsample)
+	Memr[iprob] = 0.0
+	do i = 2, nsample
+	    Memr[iprob+i-1] = Memr[iprob+i-2] + asigrl (asi, real (i-1),
+	        real (i))
+	call alimr (Memr[iprob], nsample, imin, imax)
+	call amapr (Memr[iprob], Memr[iprob], nsample, imin, imax, 0.0, 1.0)
+	call asifree (asi)
+
+	# Fit the inverse of the integral of the probability function.
+	call cvinit (cv, SPLINE3, order, 0.0, 1.0)
+	call cvfit (cv, Memr[iprob], Memr[m], Memr[w], nsample, WTS_UNIFORM,
+	    ier)
+
+	# Compute the magnitudes.
+	if (ier == OK) {
+	    do i = 1, nstars
+	        mag[i] = cveval (cv, urand (seed)) + mzero
+	} else {
+	    call printf ("Error computing the bands luminosity function.\n")
+	    call amovkr ((minmag + maxmag) / 2.0, mag, nstars)
+	}
+	call cvfree (cv)
+
+	# Free space.
+	call sfree (sp)
+end
+
+
+define	MIN_SALPETER	-4.0
+define	MAX_SALPETER	16.0
+
+# ST_SALPETER  -- Compute the Salpter luminosity function.
+
+procedure st_salpeter (mag, nstars, minmag, maxmag, mzero, nsample, order, seed)
+
+real	mag[ARB]		# array of output magnitudes
+int	nstars			# number of stars
+real	minmag, maxmag		# minimum and maximum magnitude values
+real	mzero			# zero point between relative and absolute mags.
+int	nsample			# number of points in sampling function
+int	order			# order of the spline fit
+long	seed			# value of the seed
+
+int	i, ier
+pointer	sp, m, iprob, w, cv, asi
+real	dmag, magval, imin, imax
+real	cveval(), asigrl(), urand()
+
+begin
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (m, nsample, TY_REAL)
+	call salloc (iprob, nsample, TY_REAL)
+	call salloc (w, nsample, TY_REAL)
+
+	# Compute the probability function.
+	magval = max (minmag - mzero, MIN_SALPETER)
+	dmag = (min (maxmag - mzero, MAX_SALPETER) - magval) / (nsample - 1)
+	do i = 1, nsample {
+	    Memr[m+i-1] = magval
+	    Memr[iprob+i-1] = 10.0 ** (-3.158375 + 1.550629e-1 * magval -
+	        5.19388e-3 * magval ** 2)
+	    magval = magval + dmag
+	}
+
+	# Integrate the probablity function.
+	call asiinit (asi, II_SPLINE3)
+	call asifit (asi, Memr[iprob], nsample)
+	Memr[iprob] = 0.0
+	do i = 2, nsample
+	    Memr[iprob+i-1] = Memr[iprob+i-2] + asigrl (asi, real (i-1),
+	        real (i))
+	call alimr (Memr[iprob], nsample, imin, imax)
+	call amapr (Memr[iprob], Memr[iprob], nsample, imin, imax, 0.0, 1.0)
+	call asifree (asi)
+
+	# Fit the inverse of the integral of the probability function.
+	call cvinit (cv, SPLINE3, order, 0.0, 1.0)
+	call cvfit (cv, Memr[iprob], Memr[m], Memr[w], nsample, WTS_UNIFORM,
+	    ier)
+
+	# Compute the magnitudes.
+	if (ier == OK) {
+	    do i = 1, nstars
+	        mag[i] = cveval (cv, urand (seed)) + mzero
+	} else {
+	    call printf ("Error computing the Salpeter luminosity function.\n")
+	    call amovkr ((minmag + maxmag) / 2.0, mag, nstars)
+	}
+	call cvfree (cv)
+
+	# Free space.
+	call sfree (sp)
+end
+
+
+# ST_SCHECTER -- Compute the Schecter luminosity function.
+
+procedure st_schecter (mag, nstars, alpha, mstar, minmag, maxmag, mzero,
+	nsample, order, seed)
+
+real	mag[ARB]		# array of output magnitudes
+int	nstars			# number of stars
+real	alpha, mstar		# Schecter luminosity function parameters
+real	minmag, maxmag		# minimum and maximum magnitude values
+real	mzero			# zero point between relative and absolute mags.
+int	nsample			# number of points in the sampling function
+int	order			# order of the spline fit
+long	seed			# value of the seed
+
+int	i, ier
+pointer	sp, m, iprob, w, cv, asi
+real	dmag, magval, temp, imin, imax
+real	cveval(), asigrl(), urand()
+
+begin
+	# Allocate space for fitting.
+	call smark (sp)
+	call salloc (m, nsample, TY_REAL)
+	call salloc (iprob, nsample, TY_REAL)
+	call salloc (w, nsample, TY_REAL)
+
+	# Sample the luminosity function.
+	magval = minmag - mzero
+	dmag = (maxmag - minmag) / (nsample - 1)
+	do i = 1, nsample {
+	    Memr[m+i-1] = magval
+	    temp = 0.4 * (mstar - magval)
+	    Memr[iprob+i-1] = 10.0 ** ((alpha + 1) * temp) *
+	        exp (- 10.0 ** temp)
+	    magval = magval + dmag
+	}
+
+	# Integrate the sampling function.
+	call asiinit (asi, II_SPLINE3)
+	call asifit (asi, Memr[iprob], nsample)
+	Memr[iprob] = 0.0
+	do i = 2, nsample
+	    Memr[iprob+i-1] = Memr[iprob+i-2] + asigrl (asi, real (i-1),
+	        real(i))
+	call alimr (Memr[iprob],nsample, imin, imax)
+	call amapr (Memr[iprob], Memr[iprob], nsample, imin, imax, 0.0, 1.0)
+	call asifree (asi)
+
+	# Fit the inverse of the integral of the probability function.
+	call cvinit (cv, SPLINE3, order, 0.0, 1.0)
+	call cvfit (cv, Memr[iprob], Memr[m], Memr[w], nsample, WTS_UNIFORM,
+	    ier)
+	if (ier == OK) {
+	    do i = 1, nstars
+	        mag[i] = cveval (cv, urand (seed)) + mzero
+	} else {
+	    call printf ("Error fitting the Schecter luminosity function.\n")
+	    call amovkr ((minmag + maxmag) / 2.0, mag, nstars)
+	}
+	call cvfree (cv)
+
+	# Free space.
+	call sfree (sp)
+end
+
+
+# ST_LFSAMPLE -- Compute the luminosity function using a user supplied
+# function.
+
+procedure st_lfsample (smag, mprob, nlf, mag, nstars, minmag, maxmag, nsample,
+	order, seed)
+
+real	smag[ARB]		# input array of magnitudes
+real	mprob[ARB]		# input array of relative probabilities
+int	nlf			# number of input points
+real	mag[ARB]		# output magnitude array
+int	nstars			# number of stars
+real	minmag, maxmag		# minimum and maximum magnitude values
+int	nsample			# number of sample points
+int	order			# order of the spline fit
+long	seed			# value of the seed
+
+int	npts, i, ier
+pointer	sp, m, w, iprob, cv, asi
+real	mval, dm, sfmin, sfmax, imin, imax
+real	cveval(), asigrl(), urand()
+
+begin
+	# Allocate space for fitting.
+	npts = max (nlf, nsample)
+	call smark (sp)
+	call salloc (m, nsample, TY_REAL)
+	call salloc (iprob, nsample, TY_REAL)
+	call salloc (w, npts, TY_REAL)
+
+	# Smooth the relative probability function.
+	call alimr (smag, nlf, sfmin, sfmax)
+	call cvinit (cv, SPLINE3, max (1, nlf / 4), sfmin, sfmax)
+	call cvfit (cv, smag, mprob, Memr[w], nlf, WTS_UNIFORM, ier) 
+
+	# Evaluate the smoothed function at equal intervals in r.
+	if (ier == OK) {
+	    mval = max (minmag, sfmin)
+	    dm = (min (maxmag, sfmax) - mval) / (nsample - 1)
+	    do i = 1, nsample {
+		Memr[m+i-1] = mval
+	        Memr[iprob+i-1] = cveval (cv, mval)
+	        mval = mval + dm
+	    }
+	    call cvfree (cv)
+	} else {
+	    call printf ("Error smoothing user luminosity function.\n")
+	    call amovkr ((minmag + maxmag) / 2.0, mag, nstars)
+	    call cvfree (cv)
+	    call sfree (sp)
+	}
+
+	# Evaluate the integral.
+	call asiinit (asi, II_SPLINE3)
+	call asifit (asi, Memr[iprob], nsample)
+	Memr[iprob] = 0.0
+	do i = 2, nsample
+	    Memr[iprob+i-1] = Memr[iprob+i-2] + asigrl (asi, real(i-1), real(i))
+	call alimr (Memr[iprob], nsample, imin, imax)
+	call amapr (Memr[iprob], Memr[iprob], nsample, imin, imax, 0.0, 1.0)
+	call asifree (asi)
+
+	# Fit the inverse of the integral of the probability function.
+	call cvinit (cv, SPLINE3, order, 0.0, 1.0)
+	call cvfit (cv, Memr[iprob], Memr[m], Memr[w], nsample, WTS_UNIFORM,
+	    ier)
+
+	# Sample the computed function.
+	if (ier == OK) {
+	    do i = 1, nstars
+	        mag[i] = cveval (cv, urand (seed))
+	} else {
+	    call printf ("Error computing the user luminosity function.\n")
+	    call amovkr ((minmag + maxmag) / 2.0, mag, nstars)
+	}
+	call cvfree (cv)
+
+	# Free space.
+	call sfree (sp)
+end
diff --git a/noao/artdata/lists/stmix.x b/noao/artdata/lists/stmix.x
new file mode 100644
index 00000000..4254e8eb
--- /dev/null
+++ b/noao/artdata/lists/stmix.x
@@ -0,0 +1,144 @@
+include	<math.h>
+include "starlist.h"
+
+# ST_ESMIX -- Compute the percentage of elliptical galaxies.
+
+procedure st_esmix (egal, nstar, esmix, seed)
+
+int	egal[ARB]		# array of types
+int	nstar			# number of objects
+real	esmix			# fraction of elliptical galaxies
+long	seed			# seed for random number generator
+
+int	i
+real	fraction
+real	urand()
+
+begin
+	do i = 1, nstar {
+	    fraction = urand (seed)
+	    if (fraction <= esmix)
+	        egal[i] = ST_DEVAUC
+	    else
+	        egal[i] = ST_EXP
+	}
+end
+
+
+# ST_ROUND -- Compute an array of roundness parameters.
+# For ellipical models use a uniform distribution of axial ratios.
+# For spiral models use a uniform inclination distribution.  Compute
+# the axial ratio with a .1 flatness parameter from Holmberg (Stars
+# and Stellar Systems).  Then add cosecant absorption factor to the
+# the magnitudes with a limit of 10 in the cosecant.
+
+procedure st_round (egal, mag, round, nstars, ar, alpha, seed)
+
+int	egal[ARB]		# array of galaxy types
+real	mag[ARB]		# magnitudes
+real	round[ARB]		# array of roundness values
+int	nstars			# number of stars
+real	ar			# minimum roundness value
+real	alpha			# absorption coefficent
+long	seed			# seed for the random number generator
+
+int	i
+real	dr, s, urand()
+
+begin
+	dr = (1. - ar)
+	do i = 1, nstars {
+	    if (egal[i] == ST_DEVAUC)
+	        round[i] = ar + dr * urand (seed)
+	    else {
+		s = sin (HALFPI * urand (seed))
+		round[i] = sqrt (s**2 * .99 + .01)
+		mag[i] = mag[i] + alpha * (1 / max (0.1, s) - 1) / 9.
+	    }
+	}
+
+end
+
+
+# ST_PHI -- Compute an array of position angles.
+
+procedure st_phi (phi, nstars, seed)
+
+real	phi[ARB]		# array of position angles
+int	nstars			# number of stars
+long	seed			# seed for random number generator
+
+int	i
+real	urand()
+
+begin
+	do i = 1, nstars
+	    phi[i] = urand (seed)
+	call amapr (phi, phi, nstars, 0.0, 1.0, 0.0, 360.0)
+end
+
+
+# ST_DIAMETERS -- Compute the effective diameters of the galaxies based
+# on their magnitudes.  The relation used is from Holmberg (Stars and
+# Stellar Systems).  A uniform dispersion of 20% is added.
+
+procedure st_diameters (mag, egal, axis, nstars, minmag, maxmag, eradius,
+	sradius, seed)
+
+real	mag[ARB]		# input array of magnitudes
+int	egal[ARB]		# array of galaxy types
+real	axis[ARB]		# output array of diameters
+int	nstars			# number of stars
+real	minmag			# minimum magnitude
+real	maxmag			# maximum magnitude
+real	eradius			# maximum elliptical radius
+real	sradius			# maximum spiral radius
+long	seed			# seed for random number generator
+
+int	i
+real	urand()
+
+begin
+	do i = 1, nstars {
+	    if (egal[i] == ST_DEVAUC)
+	        axis[i] = eradius * 10.0 ** ((minmag - mag[i]) / 6.)
+	    else
+	        axis[i] = sradius * eradius * 10.0 ** ((minmag - mag[i]) / 6.)
+	    axis[i] = axis[i] * (0.8 + 0.4 * urand (seed))
+	}
+end
+
+
+# ST_ZDIAMETERS -- Compute the effective diameters of the galaxies based
+# on their magnitudes and redshift. The relation used is that the redshift
+# is proportional to the luminousity and the diameters include the
+# factor of (1+z)**2.  A uniform dispersion of 50% is added.
+
+procedure st_zdiameters (mag, egal, axis, nstars, minmag, maxmag,
+	z, eradius, sradius, seed)
+
+real	mag[ARB]		# input array of magnitudes
+int	egal[ARB]		# array of galaxy types
+real	axis[ARB]		# output array of diameters
+int	nstars			# number of stars
+real	minmag			# minimum magnitude
+real	maxmag			# maximum magnitude
+real	z			# minumum redshift
+real	eradius			# maximum elliptical radius
+real	sradius			# maximum spiral radius
+long	seed			# seed for random number generator
+
+int	i
+real	z0, z1, urand()
+
+begin
+	z0 = z / (1 + z) ** 2
+	do i = 1, nstars {
+	    z1 = z * 10.0 ** (-0.2 * (minmag - mag[i]))
+	    if (egal[i] == ST_DEVAUC)
+	        axis[i] = eradius * z0 * (1 + z1) ** 2 / z1
+	    else
+	        axis[i] = sradius * eradius * z0 * (1 + z1) ** 2 / z1
+	    axis[i] = axis[i] * (0.5 + urand (seed))
+	}
+end
diff --git a/noao/artdata/lists/stplot.x b/noao/artdata/lists/stplot.x
new file mode 100644
index 00000000..2e79b6c1
--- /dev/null
+++ b/noao/artdata/lists/stplot.x
@@ -0,0 +1,1102 @@
+include <mach.h>
+include <math.h>
+include <gset.h>
+include <pkg/gtools.h>
+include "starlist.h"
+
+define	HELPFILE1  "artdata$lists/starlist.key"
+
+# ST_PLOTS -- Interactively examine the spatial density and luminosity
+# functions.
+
+procedure st_plots (sf, lf, gd, st, x, y, mag)
+
+int	sf			# spatial density file descriptor
+int	lf			# luminsosity function file descriptor
+pointer	gd			# graphics stream pointer
+pointer	st			# pointer to starlist structure
+pointer	x			# pointer to x array
+pointer	y			# pointer to y array
+pointer mag			# pointer to mag array
+
+int	wcs, key, plottype, newplot, newspace, newlum
+pointer	sp, cmd, gt1, gt2, gt3
+real	wx, wy
+int	gt_gcur()
+pointer	gt_init()
+
+begin
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Intialize the plots.
+	gt1 = gt_init()
+	gt2 = gt_init()
+	gt3 = gt_init()
+
+	newspace = NO
+	newlum = NO
+	newplot = NO
+	plottype = 1
+
+	# Draw the first plot.
+	call st_pfield (gd, gt1, st, Memr[x], Memr[y], Memr[mag], ST_NSTARS(st))
+
+	while (gt_gcur ("cursor", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+	    EOF) {
+	    switch (key) {
+
+	    case '?':
+		call gpagefile (gd, HELPFILE1, "")
+
+	    case ':':
+		call st_colon (gd, st, sf, lf, Memc[cmd], newspace, newlum,
+		    newplot)
+		switch (plottype) {
+		case 1:
+		    call gt_colon (Memc[cmd], gd, gt1, newplot)
+		case 2:
+		    call gt_colon (Memc[cmd], gd, gt2, newplot)
+		case 3:
+		    call gt_colon (Memc[cmd], gd, gt3, newplot)
+		}
+
+	    case 'q':
+		break
+
+	    case 'f':
+
+		if (newspace == YES) {
+		    call st_mkspatial (sf, st, x, y, mag)
+		    newspace = NO
+		    newplot = YES
+		}
+
+		if (newlum == YES) {
+		    call st_mklum (lf, st, x, y, mag)
+		    newlum = NO
+		    newplot = YES
+		}
+
+		if (newplot == YES) {
+		    switch (plottype) {
+		    case 1:
+			call st_pfield (gd, gt1, st, Memr[x], Memr[y],
+			    Memr[mag], ST_NSTARS(st))
+		    case 2:
+		        call st_prhist (gd, gt2, st, Memr[x], Memr[y],
+		            ST_NSTARS(st))
+		    case 3:
+			call st_pmhist (gd, gt3, st, Memr[mag], ST_NSTARS(st))
+		    default:
+			call st_pfield (gd, gt1, st, Memr[x], Memr[y],
+			    Memr[mag], ST_NSTARS(st))
+		    }
+		    newplot = NO
+		}
+
+	    case 'x':
+		if (newspace == YES || newlum == YES)
+		    call printf ("Type the f key to remake the star list\n")
+		else if (plottype != 1 || newplot == YES) {
+		    call st_pfield (gd, gt1, st, Memr[x], Memr[y], Memr[mag],
+		        ST_NSTARS(st))
+		    plottype = 1
+		    newplot = NO
+		}
+
+	    case 'r':
+		if (newspace == YES)
+		    call printf ("Type the f key to remake the star list\n")
+		else if (plottype != 2 || newplot == YES) {
+		    call st_prhist (gd, gt2, st, Memr[x], Memr[y],
+		        ST_NSTARS(st))
+		    plottype = 2
+		    newplot = NO
+		}
+
+	    case 'm':
+		if (newlum == YES)
+		    call printf ("Type the f key to remake the star list\n")
+		else if (plottype != 3 || newplot == YES) {
+		    call st_pmhist (gd, gt3, st, Memr[mag], ST_NSTARS(st))
+		    plottype = 3
+		    newplot = NO
+		}
+
+	    default:
+	    }
+	}
+
+	# Recompute the results if necessary.
+	if (newspace == YES)
+	    call st_mkspatial (sf, st, x, y, mag)
+	if (newlum == YES)
+	    call st_mklum (lf, st, x, y, mag)
+
+	# Free space for the plots.
+	call gt_free (gt1)
+	call gt_free (gt2)
+	call gt_free (gt3)
+
+	call sfree (sp)
+end
+
+
+define	HELPFILE2  "artdata$lists/gallist.key"
+
+# ST_GPLOTS -- Fit the luminosity function and make plots.
+
+procedure st_gplots (sf, lf, gd, st, x, y, mag, egal, axis, round, phi)
+
+int	sf			# spatial distribution file descriptor
+int	lf			# luminsosity distribution file descriptor
+pointer	gd			# graphics stream pointer
+pointer	st			# pointer to starlist structure
+pointer	x			# pointer to x array
+pointer	y			# pointer to y array
+pointer mag			# pointer to mag array
+pointer	egal			# pointer to type array
+pointer axis			# pointer to the diameters array
+pointer	round			# pointer to the roundness array
+pointer	phi			# pointer to the position angle array
+
+int	wcs, key, plottype
+int	newplot, newspace, newlum, newmix, newaxis, newround, newphi
+pointer	sp, cmd, gt1, gt2, gt3, gt4, gt5, gt6
+real	wx, wy
+int	gt_gcur()
+pointer	gt_init()
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Intialize the plots.
+	gt1 = gt_init()
+	gt2 = gt_init()
+	gt3 = gt_init()
+	gt4 = gt_init()
+	gt5 = gt_init()
+	gt6 = gt_init()
+
+	newspace = NO
+	newlum = NO
+	newmix = NO
+	newaxis = NO
+	newround = NO
+	newplot = NO
+	newphi = NO
+	plottype = 1
+
+	# Draw the first plot.
+	call st_pgfield (gd, gt1, st, Memr[x], Memr[y], Memr[mag], Memi[egal],
+	    Memr[axis], Memr[round], ST_NSTARS(st))
+
+	while (gt_gcur ("cursor", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+	    EOF) {
+	    switch (key) {
+
+	    case '?':
+		call gpagefile (gd, HELPFILE2, "")
+
+	    case ':':
+		call st_gcolon (gd, st, sf, lf, Memc[cmd], newspace, newlum,
+		    newmix, newaxis, newround, newphi, newplot)
+		switch (plottype) {
+		case 1:
+		    call gt_colon (Memc[cmd], gd, gt1, newplot)
+		case 2:
+		    call gt_colon (Memc[cmd], gd, gt2, newplot)
+		case 3:
+		    call gt_colon (Memc[cmd], gd, gt3, newplot)
+		case 4:
+		    call gt_colon (Memc[cmd], gd, gt4, newplot)
+		case 5:
+		    call gt_colon (Memc[cmd], gd, gt5, newplot)
+		case 6:
+		    call gt_colon (Memc[cmd], gd, gt6, newplot)
+		}
+
+	    case 'q':
+		break
+
+	    case 'f':
+
+
+		if (newphi == YES) {
+		    call st_gmkphi (st, x, y, mag, egal, axis, round, phi)
+		    newphi = NO
+		    newplot = YES
+		}
+
+		if (newspace == YES) {
+		    call st_gmkspatial (sf, st, x, y, mag, egal, axis, round,
+			phi)
+		    newspace = NO
+		    newplot = YES
+		}
+
+		if (newmix == YES) {
+		    call st_gmkmix (st, x, y, mag, egal, axis, round, phi)
+		    newmix = NO
+		    newplot = YES
+		}
+
+		if (newlum == YES) {
+		    call st_gmklum (lf, st, x, y, mag, egal, axis, round, phi)
+		    call st_gmkround (st, x, y, mag, egal, axis, round, phi)
+		    call st_gmkaxis (st, x, y, mag, egal, axis, round, phi)
+		    newlum = NO
+		    newaxis = NO
+		    newround = NO
+		    newplot = YES
+		}
+
+		if (newround == YES) {
+		    call st_gmkround (st, x, y, mag, egal, axis, round, phi)
+		    newround = NO
+		    newplot = YES
+		}
+
+		if (newaxis == YES) {
+		    call st_gmkaxis (st, x, y, mag, egal, axis, round, phi)
+		    newaxis = NO
+		    newplot = YES
+		}
+
+		if (newplot == YES) {
+		    switch (plottype) {
+		    case 1:
+			call st_pgfield (gd, gt1, st, Memr[x], Memr[y],
+			    Memr[mag], Memi[egal], Memr[axis], Memr[round],
+			    ST_NSTARS(st))
+		    case 2:
+		        call st_prhist (gd, gt2, st, Memr[x], Memr[y],
+		            ST_NSTARS(st))
+		    case 3:
+			call st_pmhist (gd, gt3, st, Memr[mag], ST_NSTARS(st))
+		    case 4:
+			call st_pdhist (gd, gt4, st, Memr[axis], ST_NSTARS(st))
+		    case 5:
+			call st_pehist (gd, gt5, st, Memr[round], ST_NSTARS(st))
+		    case 6:
+			call st_pphist (gd, gt6, st, Memr[phi], ST_NSTARS(st))
+		    default:
+			call st_pgfield (gd, gt1, st, Memr[x], Memr[y],
+			    Memr[mag], Memi[egal], Memr[axis], Memr[round],
+			    ST_NSTARS(st))
+		    }
+		    newplot = NO
+		}
+
+	    case 'x':
+		if (newspace == YES || newlum == YES || newmix == YES ||
+		    newround == YES)
+		    call printf ("Type the f key to refit the galaxies list\n")
+		else if (plottype != 1 || newplot == YES) {
+		    call st_pgfield (gd, gt1, st, Memr[x], Memr[y], Memr[mag],
+		        Memi[egal], Memr[axis], Memr[round], ST_NSTARS(st))
+		    plottype = 1
+		    newplot = NO
+		}
+
+	    case 'r':
+		if (newspace == YES)
+		    call printf ("Type the f key to refit the galaxies list\n")
+		else if (plottype != 2 || newplot == YES) {
+		    call st_prhist (gd, gt2, st, Memr[x], Memr[y],
+		        ST_NSTARS(st))
+		    plottype = 2
+		    newplot = NO
+		}
+
+	    case 'm':
+		if (newlum == YES)
+		    call printf ("Type the f key to refit the galaxies list\n")
+		else if (plottype != 3 || newplot == YES) {
+		    call st_pmhist (gd, gt3, st, Memr[mag], ST_NSTARS(st))
+		    plottype = 3
+		    newplot = NO
+		}
+
+	    case 'd':
+		if (newlum == YES)
+		    call printf ("Type the f key to refit the galaxies list\n")
+		else if (plottype != 4 || newplot == YES) {
+		    call st_pdhist (gd, gt4, st, Memr[axis], ST_NSTARS(st))
+		    plottype = 4
+		    newplot = NO
+		}
+
+	    case 'e':
+		if (newround == YES)
+		    call printf ("Type the f key to refit the galaxies list\n")
+		else if (plottype != 5 || newplot == YES) {
+		    call st_pehist (gd, gt5, st, Memr[round], ST_NSTARS(st))
+		    plottype = 5
+		    newplot = NO
+		}
+
+	    case 'p':
+		if (newphi == YES)
+		    call printf ("Type the f key to refit the galaxies list\n")
+		else if (plottype != 6 || newplot == YES) {
+		    call st_pphist (gd, gt6, st, Memr[phi], ST_NSTARS(st))
+		    plottype = 6
+		    newplot = NO
+		}
+
+	    default:
+	    }
+	}
+
+	# Recompute the functions if necessary.
+	if (newphi == YES)
+	    call st_gmkphi (st, x, y, mag, egal, axis, round, phi)
+	if (newspace == YES)
+	    call st_gmkspatial (sf, st, x, y, mag, egal, axis, round, phi)
+	if (newmix == YES)
+	    call st_gmkmix (st, x, y, mag, egal, axis, round, phi)
+	if (newlum == YES) {
+	    call st_gmklum (lf, st, x, y, mag, egal, axis, round, phi)
+	    call st_gmkaxis (st, x, y, mag, egal, axis, round, phi)
+	    call st_gmkround (st, x, y, mag, egal, axis, round, phi)
+	} else {
+	    if (newaxis == YES)
+	        call st_gmkaxis (st, x, y, mag, egal, axis, round, phi)
+	    if (newround == YES)
+	        call st_gmkround (st, x, y, mag, egal, axis, round, phi)
+	}
+
+	# Free space for the plots.
+	call gt_free (gt1)
+	call gt_free (gt2)
+	call gt_free (gt3)
+	call gt_free (gt4)
+	call gt_free (gt5)
+	call gt_free (gt6)
+
+	call sfree (sp)
+end
+
+
+# ST_PFIELD -- Plot distribution of stars in the x-y plane.
+
+procedure st_pfield (gd, gt, st, x, y, mag, npts)
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# pointer to graphics descriptor
+pointer	st		# pointer to starlist structure
+real	x[ARB]		# x coords
+real	y[ARB]		# y coords
+real	mag[ARB]	# magnitudes
+int	npts		# number of points
+
+int	i
+pointer	sp, sizes, par1, par2, title
+bool	fp_equalr()
+
+begin
+	call smark (sp)
+	call salloc (sizes, npts, TY_REAL)
+	call salloc (par1, SZ_LINE, TY_CHAR)
+	call salloc (par2, SZ_LINE, TY_CHAR)
+	call salloc (title, 3 * SZ_LINE, TY_CHAR)
+
+	# Clear screen.
+	call gclear (gd)
+
+	# Set the labels.
+	call gt_sets (gt, GTXLABEL, "X coordinate")
+	call gt_sets (gt, GTYLABEL, "Y coordinate")
+
+	# Set the title.
+	call sprintf (Memc[par1], SZ_LINE,
+	    "Nstars: %d  Spatial model: %s  Luminosity model: %s\n")
+	    call pargi (ST_NSTARS(st))
+	    if (ST_SPATIAL(st) == ST_SPFILE)
+		call pargstr (ST_SFILE(st))
+	    else
+		call pargstr (ST_SPSTRING(st))
+	    if (ST_LUMINOSITY(st) == ST_LFFILE)
+		call pargstr (ST_LFILE(st))
+	    else
+		call pargstr (ST_LFSTRING(st))
+	call sprintf (Memc[par2], SZ_LINE,
+	    "X range: %g to %g  Yrange: %g to %g  Mag range: %g to %g\n") 
+	    call pargr (ST_XMIN(st))
+	    call pargr (ST_XMAX(st))
+	    call pargr (ST_YMIN(st))
+	    call pargr (ST_YMAX(st))
+	    call pargr (ST_MINMAG(st))
+	    call pargr (ST_MAXMAG(st))
+	call sprintf (Memc[title], 3 * SZ_LINE, "%s%s%s")
+	    call pargstr ("MAP OF STARLIST\n")
+	    call pargstr (Memc[par1])
+	    call pargstr (Memc[par2])
+	call gt_sets (gt, GTTITLE, Memc[title])
+
+	# Set the plot axes min and max values.
+	call gt_setr (gt, GTXMIN, ST_XMIN(st))
+	call gt_setr (gt, GTXMAX, ST_XMAX(st))
+	call gt_setr (gt, GTYMIN, ST_YMIN(st))
+	call gt_setr (gt, GTYMAX, ST_YMAX(st))
+
+	# Set the window and labels.
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	# Plot.
+	if (fp_equalr (ST_MAXMAG(st), ST_MINMAG(st)))
+	    call amovkr (2.0, Memr[sizes], npts)
+	else
+	    call amapr (mag, Memr[sizes], npts, ST_MAXMAG(st), ST_MINMAG(st),
+		1.0, 4.0)
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_sets (gt, GTMARK, "box")
+	do i = 1, npts
+	    call gmark (gd, x[i], y[i], GM_BOX, Memr[sizes+i-1],
+		Memr[sizes+i-1])
+
+	call sfree (sp)
+end
+
+
+# ST_PGFIELD -- Plot distribution of stars in the x-y plane.
+
+procedure st_pgfield (gd, gt, st, x, y, mag, egal, axis, round, npts)
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# pointer to plot descriptor
+pointer	st		# pointer to starlist structure
+real	x[ARB]		# array of x coordinates
+real	y[ARB]		# array of y coordinates
+real	mag[ARB]	# array of magnitudes
+int	egal[ARB]	# array of galaxy types
+real	axis[ARB]	# array of diameters
+real	round[ARB]	# array of roundness values
+int	npts		# number of points
+
+int	i
+pointer	sp, par1, par2, title, xsizes
+real	amin, amax
+bool	fp_equalr()
+
+begin
+	call smark (sp)
+	call salloc (par1, SZ_LINE, TY_CHAR)
+	call salloc (par2, SZ_LINE, TY_CHAR)
+	call salloc (title, 3 * SZ_LINE, TY_CHAR)
+	call salloc (xsizes, npts, TY_REAL)
+
+	# Clear screen.
+	call gclear (gd)
+
+	# Set the labels.
+	call gt_sets (gt, GTXLABEL, "X coordinate")
+	call gt_sets (gt, GTYLABEL, "Y coordinate")
+
+	# Set the title.
+	call sprintf (Memc[par1], SZ_LINE,
+	    "Ngals: %d  Spatial model: %s  Luminosity model: %s\n")
+	    call pargi (ST_NSTARS(st))
+	    if (ST_SPATIAL(st) == ST_SPFILE)
+		call pargstr (ST_SFILE(st))
+	    else
+		call pargstr (ST_SPSTRING(st))
+	    if (ST_LUMINOSITY(st) == ST_LFFILE)
+		call pargstr (ST_LFILE(st))
+	    else
+		call pargstr (ST_LFSTRING(st))
+	call sprintf (Memc[par2], SZ_LINE,
+	    "X range: %g to %g  Yrange: %g to %g  Mag range: %g to %g\n") 
+	    call pargr (ST_XMIN(st))
+	    call pargr (ST_XMAX(st))
+	    call pargr (ST_YMIN(st))
+	    call pargr (ST_YMAX(st))
+	    call pargr (ST_MINMAG(st))
+	    call pargr (ST_MAXMAG(st))
+	call sprintf (Memc[title], 3 * SZ_LINE, "%s%s%s")
+	    call pargstr ("MAP OF GALLIST\n")
+	    call pargstr (Memc[par1])
+	    call pargstr (Memc[par2])
+	call gt_sets (gt, GTTITLE, Memc[title])
+
+	# Set the x and y axis minimums and maximums.
+	call gt_setr (gt, GTXMIN, ST_XMIN(st))
+	call gt_setr (gt, GTXMAX, ST_XMAX(st))
+	call gt_setr (gt, GTYMIN, ST_YMIN(st))
+	call gt_setr (gt, GTYMAX, ST_YMAX(st))
+
+	# Set the window and labels.
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	# Compute the marksizes.
+	call alimr (axis, npts, amin, amax)
+	if (fp_equalr (amin, amax))
+	    call amovkr (2.0, Memr[xsizes], npts)
+	else
+	    call amapr (axis, Memr[xsizes], npts, amin, amax, 1.0, 4.0)
+
+	#call amulr (axis, round, Memr[ysizes], npts)
+	#call alimr (Memr[ysizes], npts, amin, amax)
+	#call amapr (Memr[ysizes], Memr[ysizes], npts, amin, amax, 1.0, 4.0)
+
+	# Plot.
+	call gt_sets (gt, GTTYPE, "mark")
+	do i = 1, npts {
+	    if (egal[i] == ST_DEVAUC)
+	        call gmark (gd, x[i], y[i], GM_CIRCLE, Memr[xsizes+i-1],
+		    Memr[xsizes+i-1])
+	    else
+	        call gmark (gd, x[i], y[i], GM_DIAMOND, Memr[xsizes+i-1],
+		    Memr[xsizes+i-1])
+	}
+
+	call sfree (sp)
+end
+
+
+# ST_PMHIST -- Plot luminosity function of the stars or galaxies.
+
+procedure st_pmhist (gd, gt, st, mag, npts)
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# pointer to the plot descriptor
+pointer	st		# pointer to starlist structure
+real	mag[ARB]	# array of magnitudes
+int	npts		# number of points
+
+int	i, nbins
+pointer	sp, par, title, hx, hgm
+real	mval, mmin, mmax, dm, hmin, hmax
+
+begin
+	# Compute the parameters of the histogram to be plotted.
+	mmin = ST_MINMAG(st)
+	if (ST_MBINSIZE(st) <= 0.0) {
+	    nbins = 1
+	    mmax = ST_MAXMAG(st)
+	    dm = mmax - mmin
+	} else {
+	    dm = ST_MBINSIZE(st)
+	    mval = (ST_MAXMAG(st) - mmin) / dm
+	    i =  int (mval)
+	    if (abs (mval -i ) <= 100.0 * EPSILONR) {
+		nbins = i + 1
+		mmax = ST_MAXMAG(st)
+	    } else {
+		mmax = mmin + (i + 1) * dm
+		nbins = i + 2
+	    }
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (par, SZ_LINE, TY_CHAR)
+	call salloc (title, 2 * SZ_LINE, TY_CHAR)
+	call salloc (hx, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+
+	# Make the histogram.
+	call aclrr (Memr[hgm], nbins)
+	call st_hgmr (mag, npts, Memr[hgm], nbins, mmin, mmax)
+	call alimr (Memr[hgm], nbins, hmin, hmax)
+
+	# Make the histogram x scale.
+	mval = mmin + dm / 2.0
+	do i = 1, nbins {
+	    Memr[hx+i-1] = mval
+	    mval = mval + dm
+	}
+	mval = mmin + nbins * dm
+
+	# Clear the screen.
+	call gclear (gd)
+
+	# Set the labels.
+	call gt_sets (gt, GTXLABEL, "Magnitude")
+	call gt_sets (gt, GTYLABEL, "N(M)")
+
+	# Set the title.
+	call sprintf (Memc[par], SZ_LINE,
+	    "%s: %d  Model: %s\nMag: %g to %g in steps of %g\n")
+	    if (ST_TYPE(st) == ST_STARS)
+		call pargstr ("Nstars")
+	    else
+		call pargstr ("Ngals")
+	    call pargi (ST_NSTARS(st))
+	    if (ST_LUMINOSITY(st) == ST_LFFILE)
+		call pargstr (ST_LFILE(st))
+	    else
+		call pargstr (ST_LFSTRING(st))
+	    call pargr (mmin)
+	    call pargr (mmax)
+	    call pargr (dm)
+	call sprintf (Memc[title], 2 * SZ_LINE, "%s%s")
+	    call pargstr ("LUMINOSITY FUNCTION\n")
+	    call pargstr (Memc[par])
+	call gt_sets (gt, GTTITLE, Memc[title])
+
+	# Set the mins and maxs.
+	call gt_setr (gt, GTXMIN, mmin)
+	call gt_setr (gt, GTXMAX, mval)
+	call gt_setr (gt, GTYMIN, 0.0)
+	call gt_setr (gt, GTYMAX, hmax)
+
+	# Set the window.
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	# Plot.
+	call gt_sets (gt, GTTYPE, "histogram")
+	call gt_plot (gd, gt, Memr[hx], Memr[hgm], nbins)
+	call gline (gd, mmin, Memr[hgm], Memr[hx], Memr[hgm])
+	call gline (gd, Memr[hx+nbins-1], Memr[hgm+nbins-1], mval,
+	    Memr[hgm+nbins-1])
+
+	call sfree (sp)
+end
+
+
+# ST_PRHIST -- Plot radial density distribution of the stars.
+
+procedure st_prhist (gd, gt, st, x, y, npts)
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# pointer to plot descriptor
+pointer	st		# pointer to starlist structure
+real	x[ARB]		# array of x values
+real	y[ARB]		# array of y values
+int	npts		# number of points
+
+int	i, nbins
+pointer	sp, par, title, r, hx, hgm
+real	rval, rmin, rmax, dr, hmin, hmax
+
+begin
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (r, npts, TY_REAL)
+	call salloc (par, SZ_LINE, TY_CHAR)
+	call salloc (title, 2 * SZ_LINE, TY_CHAR)
+
+	# Compute the radial coordinate values.
+	do i = 1, npts {
+	    Memr[r+i-1]  = sqrt ((x[i] - ST_XC(st)) ** 2 +
+	        (y[i] - ST_YC(st)) ** 2)
+	}
+	call alimr (Memr[r], npts, rmin, rmax)
+
+	# Compute the size of the histogram.
+	rmin = 0.0
+	if (ST_RBINSIZE(st) <= 0) {
+	    nbins = 1
+	    dr = rmax - rmin
+	} else {
+	    dr = ST_RBINSIZE(st)
+	    rval = (rmax - rmin) / dr
+	    i = int (rval)
+	    if (abs (rval - i) <= 100.0 * EPSILONR)
+		nbins = i + 1
+	    else {
+		rmax = rmin + (i + 1) * dr
+		nbins = i + 2
+	    }
+	}
+
+	# Make the histogram and normalize by area.
+	call salloc (hx, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+	call aclrr (Memr[hgm], nbins)
+	call st_hgmr (Memr[r], npts, Memr[hgm], nbins, rmin, rmax)
+	rval = rmin
+	do i = 1, nbins {
+	    Memr[hgm+i-1] = Memr[hgm+i-1] / ( PI * dr * (2.0 * rval + dr))
+	    rval = rval + dr
+	}
+	call alimr (Memr[hgm], nbins, hmin, hmax)
+
+	# Make the histogram x scale.
+	rval = rmin + dr / 2.0
+	do i = 1, nbins {
+	    Memr[hx+i-1] = rval
+	    rval = rval + dr
+	}
+	rval = rmin + nbins * dr
+
+	# Clear the screen.
+	call gclear (gd)
+
+	# Set the labels.
+	call gt_sets (gt, GTXLABEL, "Radial Distance")
+	call gt_sets (gt, GTYLABEL, "N(R) / Area")
+
+	# Set the title.
+	call sprintf (Memc[par], SZ_LINE,
+	    "Model: %s  %s: %d\nRadius: %g to %g in steps of %g\n")
+	    if (ST_SPATIAL(st) == ST_SPFILE)
+		call pargstr (ST_SFILE(st))
+	    else
+		call pargstr (ST_SPSTRING(st))
+	    if (ST_TYPE(st) == ST_STARS)
+		call pargstr ("Nstars")
+	    else
+		call pargstr ("Ngals")
+	    call pargi (ST_NSTARS(st))
+	    call pargr (rmin)
+	    call pargr (rmax)
+	    call pargr (dr)
+	call sprintf (Memc[title], 2 * SZ_LINE, "%s%s")
+	    call pargstr ("RADIAL DENSITY FUNCTION\n")
+	    call pargstr (Memc[par])
+	call gt_sets (gt, GTTITLE, Memc[title])
+
+	# Set the x and y axes minimum and maximum values.
+	call gt_setr (gt, GTXMIN, rmin) 
+	call gt_setr (gt, GTXMAX, rval)
+	call gt_setr (gt, GTYMIN, 0.0)
+	call gt_setr (gt, GTYMAX, hmax)
+
+	# Set the window.
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	# Plot.
+	call gt_sets (gt, GTTYPE, "histogram")
+	call gt_plot (gd, gt, Memr[hx], Memr[hgm], nbins)
+	call gline (gd, rmin, Memr[hgm], Memr[hx], Memr[hgm])
+	call gline (gd, Memr[hx+nbins-1], Memr[hgm+nbins-1], rval,
+	    Memr[hgm+nbins-1])
+
+	call sfree (sp)
+end
+
+
+# ST_PDHIST -- Plot the distribution of galaxy diameters.
+
+procedure st_pdhist (gd, gt, st, axis, npts)
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# pointer to plot descriptor
+pointer	st		# pointer to starlist structure
+real	axis[ARB]	# array of diameters
+int	npts		# number of points
+
+int	i, nbins
+pointer	sp, par, title, hx, hgm
+real	aval, amin, amax, da, hmin, hmax
+
+begin
+
+	# Allocate space for the histogram.
+	call alimr (axis, npts, amin, amax)
+	amax = max (ST_ERADIUS(st), ST_SRADIUS(st))
+	if (ST_DBINSIZE(st) <= 0) {
+	    nbins = 1
+	    da = amax - amin
+	} else {
+	    da = ST_DBINSIZE(st)
+	    aval = (amax - amin) / da
+	    i = int (aval)
+	    if (abs (aval - i) <= 100.0 * EPSILONR)
+		nbins = i + 1
+	    else {
+		amin = amax - (i + 1) * da
+		nbins = i + 2
+	    }
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (par, SZ_LINE, TY_CHAR)
+	call salloc (title, 2 * SZ_LINE, TY_CHAR)
+	call salloc (hx, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+
+	# Make the histogram.
+	call aclrr (Memr[hgm], nbins)
+	call st_hgmr (axis, npts, Memr[hgm], nbins, amin, amax)
+	call alimr (Memr[hgm], nbins, hmin, hmax)
+
+	# Make the histogram x scale.
+	aval = amin + da / 2.0
+	do i = 1, nbins {
+	    Memr[hx+i-1] = aval
+	    aval = aval + da
+	}
+	aval = amin + nbins * da
+
+	# Clear the screen.
+	call gclear (gd)
+
+	# Set the labels.
+	call gt_sets (gt, GTXLABEL, "Half-Flux Radius")
+	call gt_sets (gt, GTYLABEL, "N(radius)")
+
+	# Set the title.
+	call sprintf (Memc[par], SZ_LINE,
+     "Luminosity function: %s  Ngals: %d\nDiameter: %g to %g in steps of %g\n")
+	    if (ST_LUMINOSITY(st) == ST_LFFILE)
+		call pargstr (ST_LFILE(st))
+	    else
+		call pargstr (ST_LFSTRING(st))
+	    call pargi (ST_NSTARS(st))
+	    call pargr (amin)
+	    call pargr (amax)
+	    call pargr (da)
+	call sprintf (Memc[title], 2 * SZ_LINE, "%s%s")
+	    call pargstr ("HALF-FLUX RADIUS DISTRIBUTION\n")
+	    call pargstr (Memc[par])
+	call gt_sets (gt, GTTITLE, Memc[title])
+
+	# Set the mins and maxs.
+	call gt_setr (gt, GTXMIN, amin)
+	call gt_setr (gt, GTXMAX, aval)
+	call gt_setr (gt, GTYMIN, 0.0)
+	call gt_setr (gt, GTYMAX, hmax)
+
+	# Set the window.
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	# Plot.
+	call gt_sets (gt, GTTYPE, "histogram")
+	call gt_plot (gd, gt, Memr[hx], Memr[hgm], nbins)
+	call gline (gd, amin, Memr[hgm], Memr[hx], Memr[hgm])
+	call gline (gd, Memr[hx+nbins-1], Memr[hgm+nbins-1], aval,
+	    Memr[hgm+nbins-1])
+
+	call sfree (sp)
+end
+
+
+# ST_PEHIST -- Plot the distribution of galaxy diameters.
+
+procedure st_pehist (gd, gt, st, round, npts)
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# pointer to plot descriptor
+pointer	st		# pointer to starlist structure
+real	round[ARB]	# array of roundness values
+int	npts		# number of points
+
+int	i, nbins
+pointer	sp, par, title, hx, hgm
+real	eval, emin, emax, de, hmin, hmax
+
+begin
+	# Compute the size of the histogram.
+	emin = .1
+	if (ST_EBINSIZE(st) <= 0) {
+	    nbins = 1
+	    emax = 1.
+	    de = emax - emin
+	} else {
+	    de = ST_EBINSIZE(st)
+	    eval = (1. - emin) / de
+	    i = int (eval)
+	    if (abs (eval - i) <= 100.0 * EPSILONR) {
+		nbins = i + 1
+		emax = 1.
+	    } else {
+		nbins = i + 2
+		emax = emin + (i + 1) * de
+	    }
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (par, SZ_LINE, TY_CHAR)
+	call salloc (title, 2 * SZ_LINE, TY_CHAR)
+	call salloc (hx, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+
+	# Make the histogram and normalize by area.
+	call aclrr (Memr[hgm], nbins)
+	call st_hgmr (round, npts, Memr[hgm], nbins, emin, emax)
+	call alimr (Memr[hgm], nbins, hmin, hmax)
+
+	# Make the histogram x scale.
+	eval = emin + de / 2.0
+	do i = 1, nbins {
+	    Memr[hx+i-1] = eval
+	    eval = eval + de
+	}
+	eval = emin + nbins * de
+
+	# Clear the screen.
+	call gclear (gd)
+
+	# Set the labels.
+	call gt_sets (gt, GTXLABEL, "Axial Ratio")
+	call gt_sets (gt, GTYLABEL, "N(Axial Ratio)")
+
+	# Set the title.
+	call sprintf (Memc[par], SZ_LINE,
+	    "Ngals: %d\nRoundness: %g to %g in steps of %g\n")
+	    call pargi (ST_NSTARS(st))
+	    call pargr (emin)
+	    call pargr (emax)
+	    call pargr (de)
+	call sprintf (Memc[title], 2 * SZ_LINE, "%s%s")
+	    call pargstr ("AXIAL RATIO DISTRIBUTION\n")
+	    call pargstr (Memc[par])
+	call gt_sets (gt, GTTITLE, Memc[title])
+
+	# Set the mins and maxs.
+	call gt_setr (gt, GTXMIN, emin)
+	call gt_setr (gt, GTXMAX, eval)
+	call gt_setr (gt, GTYMIN, 0.0)
+	call gt_setr (gt, GTYMAX, hmax)
+
+	# Set the window.
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	# Plot.
+	call gt_sets (gt, GTTYPE, "histogram")
+	call gt_plot (gd, gt, Memr[hx], Memr[hgm], nbins)
+	call gline (gd, emin, Memr[hgm], Memr[hx], Memr[hgm])
+	call gline (gd, Memr[hx+nbins-1], Memr[hgm+nbins-1], eval,
+	    Memr[hgm+nbins-1])
+
+	call sfree (sp)
+end
+
+
+define	MIN_ANGLE	0.0
+define	MAX_ANGLE	360.0
+
+# ST_PPHIST -- Plot the distribution of galaxy diameters.
+
+procedure st_pphist (gd, gt, st, phi, npts)
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# pointer to plot descriptor
+pointer	st		# pointer to starlist structure
+real	phi[ARB]	# array of position angle values
+int	npts		# number of points
+
+int	i, nbins
+pointer	sp, par, title, hx, hgm
+real	pval, pmin, pmax, dp, hmin, hmax
+
+begin
+	# Compute the size of the histogram.
+	pmin = MIN_ANGLE
+	if (ST_PBINSIZE(st) <= 0) {
+	    nbins = 1
+	    pmax = MAX_ANGLE
+	    dp = pmax - pmin
+	} else {
+	    dp = ST_PBINSIZE(st)
+	    pval = (MAX_ANGLE - pmin) / dp
+	    i = int (pval)
+	    if (abs (pval - i) <= 100.0 * EPSILONR) {
+		nbins = i + 1
+		pmax = MAX_ANGLE
+	    } else {
+	        pmax = MIN_ANGLE + (i + 1) * dp
+		nbins = i + 2
+	    }
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (par, SZ_LINE, TY_CHAR)
+	call salloc (title, 2 * SZ_LINE, TY_CHAR)
+	call salloc (hx, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+
+	# Make the histogram.
+	call aclrr (Memr[hgm], nbins)
+	call st_hgmr (phi, npts, Memr[hgm], nbins, pmin, pmax)
+	call alimr (Memr[hgm], nbins, hmin, hmax)
+
+
+	# Make the histogram x scale.
+	pval = dp / 2.0
+	do i = 1, nbins {
+	    Memr[hx+i-1] = pval
+	    pval = pval + dp
+	}
+	pval = pmin + nbins * dp
+
+	# Clear the screen.
+	call gclear (gd)
+
+	# Set the axis labels.
+	call gt_sets (gt, GTXLABEL, "Position Angle")
+	call gt_sets (gt, GTYLABEL, "N(Position Angle)")
+
+	# Set the title.
+	call sprintf (Memc[par], SZ_LINE,
+	    "Ngals: %d\nPosition angle: %g to %g in steps of %g\n")
+	    call pargi (ST_NSTARS(st))
+	    call pargr (pmin)
+	    call pargr (pmax)
+	    call pargr (dp)
+	call sprintf (Memc[title], 2 * SZ_LINE, "%s%s")
+	    call pargstr ("POSITION ANGLE DISTRIBUTION\n")
+	    call pargstr (Memc[par])
+	call gt_sets (gt, GTTITLE, Memc[title])
+
+	# Set the axis mins and maxs.
+	call gt_setr (gt, GTXMIN, pmin)
+	call gt_setr (gt, GTXMAX, pval)
+	call gt_setr (gt, GTYMIN, 0.0)
+	call gt_setr (gt, GTYMAX, hmax)
+
+	# Set the window.
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	# Plot.
+	call gt_sets (gt, GTTYPE, "histogram")
+	call gt_plot (gd, gt, Memr[hx], Memr[hgm], nbins)
+	call gline (gd, pmin, Memr[hgm], Memr[hx], Memr[hgm])
+	call gline (gd, Memr[hx+nbins-1], Memr[hgm+nbins-1], pval,
+	    Memr[hgm+nbins-1])
+
+	call sfree (sp)
+end
+
+
+# ST_HGMR -- Accumulate the histogram of the input vector.  The output vector
+# hmg (the histogram) should be cleared prior to the first call.
+
+procedure st_hgmr (data, npix, hgm, nbins, z1, z2)
+
+real 	data[ARB]		# data vector
+int	npix			# number of pixels
+real	hgm[ARB]		# output histogram
+int	nbins			# number of bins in histogram
+real	z1, z2			# greyscale values of first and last bins
+
+int	bin, i
+real	z, dz
+
+begin
+	dz = real (nbins - 1) / real (z2 - z1)
+	if (abs (dz - 1.0) < (EPSILONR * 2.0)) {
+	    do i = 1, npix {
+		z = data[i]
+	        if (z < z1 || z > z2)
+		    next
+	        bin = int (z - z1) + 1
+	        hgm[bin] = hgm[bin] + 1.0
+	    }
+	} else {
+	    do i = 1, npix {
+		z = data[i]
+	        if (z < z1 || z > z2)
+		    next
+	        bin = int ((z - z1) * dz) + 1
+	        hgm[bin] = hgm[bin] + 1.0
+	    }
+	}
+end
diff --git a/noao/artdata/lists/stshow.x b/noao/artdata/lists/stshow.x
new file mode 100644
index 00000000..4aa8d853
--- /dev/null
+++ b/noao/artdata/lists/stshow.x
@@ -0,0 +1,142 @@
+include "starlist.h"
+
+procedure st_show (st)
+
+pointer st		# pointer to starlist structure
+
+begin
+	call printf ("nstars = %d\n")
+	    call pargi (ST_NSTARS(st))
+
+	call printf ("spatial: %s\n")
+	    call pargstr (ST_SPSTRING(st))
+	call printf ("    xmin = %g  xmax = %g\n")
+	    call pargr (ST_XMIN(st))
+	    call pargr (ST_XMAX(st))
+	call printf ("    ymin = %g  ymax = %g\n")
+	    call pargr (ST_YMIN(st))
+	    call pargr (ST_YMAX(st))
+	call printf ("    xc = %g  yc = %g\n")
+	    call pargr (ST_XC(st))
+	    call pargr (ST_YC(st))
+	if (ST_SPATIAL(st) == ST_HUBBLE) {
+	    call printf ("    core = %g  base = %g\n")
+	        call pargr (ST_CORE(st))
+	        call pargr (ST_BASE(st))
+	} else if (ST_SPATIAL(st) == ST_SPFILE) {
+	    call printf ("    file: %s\n")
+	        call pargstr (ST_SFILE(st))
+	}
+	call printf ("    nssample = %d  sorder = %d\n")
+	    call pargi (ST_NSSAMPLE(st))
+	    call pargi (ST_SORDER(st))
+	call printf ("    rbinsize = %d\n")
+	    call pargr (ST_RBINSIZE(st))
+
+	call printf ("luminosity = %s\n")
+	    call pargstr (ST_LFSTRING(st))
+	call printf ("    minmag = %g  maxmag = %g\n")
+	    call pargr (ST_MINMAG(st))
+	    call pargr (ST_MAXMAG(st))
+	if (ST_LUMINOSITY(st) == ST_BANDS) {
+	    call printf ( "    bands: alpha = %g  beta = %g\n") 
+	        call pargr (ST_ALPHA(st))
+	        call pargr (ST_BETA(st))
+	    call printf ("    delta = %g  mstar = %g  mzero = %g\n")
+	        call pargr (ST_DELTA(st))
+	        call pargr (ST_MSTAR(st))
+		call pargr (ST_MZERO(st))
+	} else if (ST_LUMINOSITY(st) == ST_SALPETER) {
+	    call printf ("    mzero = %g\n")
+		call pargr (ST_MZERO(st))
+	} else if (ST_LUMINOSITY(st) == ST_POWLAW) {
+	    call printf ("    power = %g\n")
+		call pargr (ST_POWER(st))
+	} else if (ST_LUMINOSITY(st) == ST_LFFILE) {
+	    call printf ("    file: %s\n")
+	        call pargstr (ST_LFILE(st))
+	}
+	call printf ( "    nlsample = %d  lorder = %d\n")
+	    call pargl (ST_NLSAMPLE(st))
+	    call pargi (ST_LORDER(st))
+	call printf ("    mbinsize = %g\n")
+	    call pargr (ST_MBINSIZE(st))
+end
+
+
+# ST_GSHOW -- Display the GALAXIES parameters.
+
+procedure st_gshow (st)
+
+pointer st		# pointer to starfield structure
+
+begin
+	call printf ("ngals = %d\n")
+	    call pargi (ST_NSTARS(st))
+
+	call printf ("spatial: %s\n")
+	    call pargstr (ST_SPSTRING(st))
+	call printf ("    xmin = %g  xmax = %g\n")
+	    call pargr (ST_XMIN(st))
+	    call pargr (ST_XMAX(st))
+	call printf ("    ymin = %g  ymax = %g\n")
+	    call pargr (ST_YMIN(st))
+	    call pargr (ST_YMAX(st))
+	call printf ("    xc = %g  yc = %g\n")
+	    call pargr (ST_XC(st))
+	    call pargr (ST_YC(st))
+	if (ST_SPATIAL(st) == ST_HUBBLE) {
+	    call printf ("    core = %g  base = %g\n")
+	        call pargr (ST_CORE(st))
+	        call pargr (ST_BASE(st))
+	} else if (ST_SPATIAL(st) == ST_SPFILE) {
+	    call printf ("    file: %s\n")
+	        call pargstr (ST_SFILE(st))
+	}
+	call printf ("    nssample = %d  sorder = %d\n")
+	    call pargi (ST_NSSAMPLE(st))
+	    call pargi (ST_SORDER(st))
+	call printf ("    rbinsize = %d\n")
+	    call pargr (ST_RBINSIZE(st))
+
+	call printf ("luminosity = %s\n")
+	    call pargstr (ST_LFSTRING(st))
+	call printf ("    minmag = %g  maxmag = %g\n")
+	    call pargr (ST_MINMAG(st))
+	    call pargr (ST_MAXMAG(st))
+	if (ST_LUMINOSITY(st) == ST_POWLAW) {
+	    call printf ( "    powlaw: power = %g\n") 
+	        call pargr (ST_POWER(st))
+	} else if (ST_LUMINOSITY(st) == ST_SCHECTER) {
+	    call printf ( "    schecter: alpha = %g  mstar = %g  mzero = %g\n") 
+	        call pargr (ST_ALPHA(st))
+	        call pargr (ST_MSTAR(st))
+		call pargr (ST_MZERO(st))
+	} else if (ST_LUMINOSITY(st) == ST_LFFILE) {
+	    call printf ("    file: %s\n")
+	        call pargstr (ST_LFILE(st))
+	}
+	call printf ("    nlsample = %d  lorder = %d\n")
+	    call pargi (ST_NLSAMPLE(st))
+	    call pargi (ST_LORDER(st))
+	call printf ("    mbinsize = %d\n")
+	    call pargr (ST_MBINSIZE(st))
+	call printf ("    eradius = %g  sradius = %g  dbinsize = %g\n")
+	    call pargr (ST_ERADIUS(st))
+	    call pargr (ST_SRADIUS(st))
+	    call pargr (ST_DBINSIZE(st))
+	call printf ("    z = %g\n")
+	    call pargr (ST_Z(st))
+	call printf ("    absorption = %g\n")
+	    call pargr (ST_ABSORPTION(st))
+
+	call eprintf ("egalmix = %g\n")
+	    call pargr (ST_EGALMIX(st))
+	call printf ("ar = %g\n")
+	    call pargr (ST_AR(st))
+	call printf ("    ebinsize = %g\n")
+	    call pargr (ST_EBINSIZE(st))
+	call printf ("posmin = 0.0  posmax = 360.0\n")
+	call printf ("    pbinsize = %g\n")
+	    call pargr (ST_PBINSIZE(st))
+end
diff --git a/noao/artdata/lists/stspatial.x b/noao/artdata/lists/stspatial.x
new file mode 100644
index 00000000..1171dfda
--- /dev/null
+++ b/noao/artdata/lists/stspatial.x
@@ -0,0 +1,264 @@
+include <math.h>
+include <math/curfit.h>
+include <math/iminterp.h>
+
+
+# ST_XYUNIFORM -- Compute a set of x and y values uniformly distributed in x and
+# y between xmin, xmax, ymin and ymax.
+
+procedure st_xyuniform (x, y, nstars, xmin, xmax, ymin, ymax, seed)
+
+real	x[ARB]		# output array of x values
+real	y[ARB]		# output array of y values
+int	nstars		# number of stars
+real	xmin, xmax	# x coordinate limits
+real	ymin, ymax	# y coordinate limits
+long	seed		# seed for random number generator
+
+int	i
+real	urand()
+
+begin
+	# Compute x and y values between 0 and 1.
+	do i = 1, nstars {
+	    x[i] = urand (seed)
+	    y[i] = urand (seed)
+	}
+
+	# Map these values into the data range.
+	call amapr (x, x, nstars, 0.0, 1.0, xmin, xmax)
+	call amapr (y, y, nstars, 0.0, 1.0, ymin, ymax)
+end
+
+
+# ST_HBSAMPLE -- Compute a set of x and y values with a Hubble density
+# distribution.
+
+procedure st_hbsample (x, y, nstars, core, base, xc, yc, xmin, xmax, ymin, ymax,
+        nsample, order, seed)
+
+real	x[ARB]			# output array of x values
+real	y[ARB]			# output array of y values
+int	nstars			# number of stars
+real	core			# Hubble core radius
+real	base			# baseline density
+real	xc, yc			# x and y center coordinates
+real	xmin, xmax		# x range
+real	ymin, ymax		# y range
+int	nsample			# number of sample points
+int	order			# order of spline fit
+long	seed			# seed for random number generator
+
+int	i, ier
+pointer	sp, rad, prob, w, cv
+real	r1, r2, r3, r4, rmin, rmax, rval, dr, theta
+real	urand(), cveval()
+
+begin
+	# Allocate space for the fits.
+	call smark (sp)
+	call salloc (rad, nsample, TY_REAL)
+	call salloc (prob, nsample, TY_REAL)
+	call salloc (w, nsample, TY_REAL)
+
+	# Compute the maximum radial distance from the center and
+	# the sampling interval.
+
+	r1 = (xmin - xc) ** 2 + (ymin - yc) ** 2
+	r2 = (xmax - xc) ** 2 + (ymin - yc) ** 2
+	r3 = (xmax - xc) ** 2 + (ymax - yc) ** 2
+	r4 = (xmin - xc) ** 2 + (ymax - yc) ** 2
+	if (xc >= xmin && xc <= xmax && yc >= ymin && yc <= ymax)
+	    rmin = 0.0
+	else if (yc >= ymin && yc <= ymax)
+	    rmin = min (abs (xmin - xc), abs (xmax - xc))
+	else if (xc >= xmin && xc <= xmax)
+	    rmin = min (abs (ymin - yc), abs (ymax - yc))
+	else
+	    rmin = sqrt (min (r1, r2, r3, r4))
+	rmax = sqrt (max (r1, r2, r3, r4))
+	dr = (rmax - rmin) / (nsample - 1)
+
+	# Compute the integral of the sampling function.
+	r1 = core ** 2
+	rval = rmin
+	do i = 1, nsample {
+	    Memr[rad+i-1] = rval
+	    r2 = (core + rval) / core
+	    Memr[prob+i-1] = r1 * (log (r2) + 1.0 / r2 - 1.0) +
+	        base * rval ** 2 / 2.0
+	    rval = rval + dr
+	}
+
+	# Normalize the probability function.
+	call alimr (Memr[prob], nsample, rmin, rmax)
+	call amapr (Memr[prob], Memr[prob], nsample, rmin, rmax, 0.0, 1.0)
+
+	# Fit the inverse of the integral of the probability function
+	call cvinit (cv, SPLINE3, order, 0.0, 1.0)
+	call cvfit (cv, Memr[prob], Memr[rad], Memr[w], nsample, WTS_UNIFORM,
+	    ier)
+
+	# Sample the computed function.
+	if (ier == OK) {
+	    i = 0
+	    repeat {
+	        rval = cveval (cv, urand (seed))
+	        theta = DEGTORAD (360.0 * urand (seed))
+	        x[i+1] = rval * cos (theta) + xc
+	        y[i+1] = rval * sin (theta) + yc
+	        if (x[i+1] >= xmin && x[i+1] <= xmax && y[i+1] >= ymin &&
+	            y[i+1] <= ymax)
+		    i = i + 1
+	    } until (i >= nstars)
+	} else {
+	    call amovkr ((xmin + xmax) / 2.0, x, nstars)
+	    call amovkr ((ymin + ymax) / 2.0, y, nstars)
+	    call printf ("Error computing the spatial probability function.\n")
+	}
+
+	# Free up the space.
+	call cvfree (cv)
+	call sfree (sp)
+end
+
+
+# ST_SFSAMPLE -- Compute a sample of x and y coordinate values based
+# on a user supplied spatial density function.
+
+procedure st_sfsample (r, rprob, nsf, x, y, nstars, nsample, order, xc, yc,
+	xmin, xmax, ymin, ymax, seed)
+
+real	r[ARB]			# input array of radii
+real	rprob[ARB]		# input array of relative probabilities
+int	nsf			# number of input points
+real	x[ARB]			# output x coordinate array
+real	y[ARB]			# output y coordinate array
+int	nstars			# number of stars
+int	nsample			# number of sample points
+int	order			# order of the spline fit
+real	xc, yc			# x and y center coordiantes
+real	xmin, xmax		# min and max x values
+real	ymin, ymax		# min and max y values
+long	seed			# value of the seed
+
+int	itemp, i, ier
+pointer	sp, w, rad, iprob, cv, asi
+real	rfmin, rfmax, dr, rval, theta, imin, imax
+real	cveval(), asigrl(), urand()
+
+begin
+	# Allocate space for fitting.
+	itemp = max (nsf, nsample)
+	call smark (sp)
+	call salloc (rad, nsample, TY_REAL)
+	call salloc (iprob, nsample, TY_REAL)
+	call salloc (w, itemp, TY_REAL)
+
+	# Smooth the relative probability function function.
+	call alimr (r, nsf, rfmin, rfmax)
+	itemp = min (order, max (1, nsf / 4))
+	call cvinit (cv, SPLINE3, itemp, rfmin, rfmax)
+	call cvfit (cv, r, rprob, Memr[w], nsf, WTS_UNIFORM, ier) 
+
+	# Evaluate the smoothed function at equal intervals in r,
+	# Multiplying by r to prepare for the area integration.
+	if (ier == OK) {
+	    rval = rfmin
+	    dr = (rfmax - rfmin) / (nsample - 1)
+	    do i = 1, nsample {
+		Memr[rad+i-1] = rval
+	        Memr[iprob+i-1] = rval * cveval (cv, rval)
+	        rval = rval + dr
+	    }
+	    call cvfree (cv)
+	} else {
+	    call printf ("Error smoothing the user spatial density function.\n")
+	    call amovkr ((xmin + xmax) / 2.0, x, nstars)
+	    call amovkr ((ymin + ymax) / 2.0, y, nstars)
+	    call cvfree (cv)
+	    call sfree (sp)
+	    return
+	}
+
+
+	# Evaluate the integral.
+	call asiinit (asi, II_SPLINE3)
+	call asifit (asi, Memr[iprob], nsample)
+	Memr[iprob] = 0.0
+	do i = 2, nsample
+	    Memr[iprob+i-1] = Memr[iprob+i-2] + asigrl (asi, real (i - 1),
+		real (i))
+	call alimr (Memr[iprob], nsample, imin, imax)
+	call amapr (Memr[iprob], Memr[iprob], nsample, imin, imax, 0.0, 1.0)
+	call asifree (asi)
+
+	# Fit the inverse of the integral of the probability function.
+	call cvinit (cv, SPLINE3, order, 0.0, 1.0)
+	call cvfit (cv, Memr[iprob], Memr[rad], Memr[w], nsample, WTS_UNIFORM,
+	    ier)
+
+	# Sample the computed function.
+	if (ier == OK) {
+	    i = 0
+	    repeat {
+	        rval = cveval (cv, urand (seed))
+	        theta = DEGTORAD (360.0 * urand (seed))
+	        x[i+1] = rval * cos (theta) + xc
+	        y[i+1] = rval * sin (theta) + yc
+	        if (x[i+1] >= xmin && x[i+1] <= xmax && y[i+1] >= ymin &&
+	            y[i+1] <= ymax)
+		    i = i + 1
+	    } until (i >= nstars)
+	} else {
+	    call printf (
+	    "Error fitting the spatial probability function.\n")
+	    call amovkr ((xmin + xmax) / 2.0, x, nstars)
+	    call amovkr ((ymin + ymax) / 2.0, y, nstars)
+	}
+	call cvfree (cv)
+
+	# Free space.
+	call sfree (sp)
+end
+
+
+define	BUFSIZE		200
+
+# ST_GFETCHXY -- Fetch two real values from a text file.
+
+int procedure st_gfetchxy (sf, x, y)
+
+int	sf			# input text file descriptor
+pointer	x			# pointer to the x array
+pointer	y			# pointer to the y array
+
+int	bufsize, npts
+int	fscan(), nscan()
+
+begin
+	call seek (sf, BOF)
+
+	call malloc (x, BUFSIZE, TY_REAL)
+	call malloc (y, BUFSIZE, TY_REAL)
+	bufsize = BUFSIZE
+
+	npts = 0
+	while (fscan (sf) != EOF) {
+	    call gargr (Memr[x+npts])
+	    call gargr (Memr[y+npts])
+	    if (nscan () != 2)
+		next
+	    npts = npts + 1
+	    if (npts < bufsize)
+		next
+	    bufsize = bufsize + BUFSIZE
+	    call realloc (x, bufsize, TY_REAL)
+	    call realloc (y, bufsize, TY_REAL)
+	}
+
+	call realloc (x, npts, TY_REAL)
+	call realloc (y, npts, TY_REAL)
+
+	return (npts)
+end
diff --git a/noao/artdata/lists/t_gallist.x b/noao/artdata/lists/t_gallist.x
new file mode 100644
index 00000000..48f1f333
--- /dev/null
+++ b/noao/artdata/lists/t_gallist.x
@@ -0,0 +1,453 @@
+include <fset.h>
+include "starlist.h"
+
+procedure t_gallist()
+
+pointer	galaxies		# pointer to the name of the output file
+pointer	graphics		# poionter to graphics device name
+
+int	sf, lf
+long	seed, sseed, lseed
+long	sseed1, lseed1
+pointer	sp, str, x, y, mag, egal, axis, round, phi, dt, st, gd
+
+bool	clgetb()
+int	clgeti(), clgwrd(), open()
+long	clgetl(), clktime()
+pointer	dtmap(), gopen()
+real	clgetr()
+
+begin
+	call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (galaxies, SZ_FNAME, TY_CHAR)
+	call salloc (graphics, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Allocate the starlist / galaxies structure.
+	call malloc (st, LEN_STSTRUCT, TY_STRUCT)
+	ST_TYPE(st) = ST_GALAXIES
+
+	# Get the parameters.
+	call clgstr ("gallist", Memc[galaxies], SZ_FNAME)
+	ST_NSTARS(st) = clgeti ("ngals")
+
+	# Get the parameters of the spatial density function.
+	ST_SPATIAL(st) = clgwrd ("spatial", ST_SPSTRING(st), SZ_FNAME, SPFUNCS)
+	ST_XMIN(st) = clgetr ("xmin")
+	ST_XMAX(st) = clgetr ("xmax")
+	ST_YMIN(st) = clgetr ("ymin")
+	ST_YMAX(st) = clgetr ("ymax")
+	ST_SFILE(st) = EOS
+	switch (ST_SPATIAL(st)) {
+	case ST_UNIFORM:
+	    sf = NULL
+	    ST_XC(st) = (ST_XMAX(st) + ST_XMIN(st)) / 2.0
+	    ST_YC(st) = (ST_YMAX(st) + ST_YMIN(st)) / 2.0
+	case ST_HUBBLE:
+	    sf = NULL
+	    ST_XC(st) = clgetr ("xcenter")
+	    if (IS_INDEFR(ST_XC(st)))
+	        ST_XC(st) = (ST_XMAX(st) + ST_XMIN(st)) / 2.0
+	    ST_YC(st) = clgetr ("ycenter")
+	    if (IS_INDEFR(ST_YC(st)))
+	        ST_YC(st) = (ST_YMAX(st) + ST_YMIN(st)) / 2.0
+	case ST_SPFILE:
+	    call clgstr ("sfile", ST_SFILE(st), SZ_FNAME)
+	    sf = open (ST_SFILE(st), READ_ONLY, TEXT_FILE)
+	    ST_XC(st) = clgetr ("xcenter")
+	    if (IS_INDEFR(ST_XC(st)))
+	        ST_XC(st) = (ST_XMAX(st) + ST_XMIN(st)) / 2.0
+	    ST_YC(st) = clgetr ("ycenter")
+	    if (IS_INDEFR(ST_YC(st)))
+	        ST_YC(st) = (ST_YMAX(st) + ST_YMIN(st)) / 2.0
+	}
+	ST_CORE(st) = clgetr ("core_radius")
+	ST_BASE(st) = clgetr ("base")
+
+
+	# Get the parameters of the luminosity function.
+	ST_LUMINOSITY(st) = clgwrd ("luminosity", ST_LFSTRING(st), SZ_FNAME,
+	    GLUMFUNCS)
+	ST_MINMAG(st) = clgetr ("minmag")
+	ST_MAXMAG(st) = clgetr ("maxmag")
+	ST_LFILE(st) = EOS
+	switch (ST_LUMINOSITY(st)) {
+	case ST_UNIFORM, ST_POWLAW, ST_SCHECTER:
+	    lf = NULL
+	case ST_LFFILE:
+	    call clgstr ("lfile", ST_LFILE(st), SZ_FNAME)
+	    lf = open (ST_LFILE(st), READ_ONLY, TEXT_FILE)
+	}
+	ST_POWER(st) = clgetr ("power")
+	ST_MZERO(st) = clgetr ("mzero")
+	ST_ALPHA(st) = clgetr ("alpha")
+	ST_MSTAR(st) = clgetr ("mstar")
+
+	# Get the remaining parameters.
+	ST_Z(st) = clgetr ("z")
+	ST_AR(st) = clgetr ("ar")
+	ST_ERADIUS(st) = clgetr ("eradius")
+	ST_SRADIUS(st) = clgetr ("sradius")
+	ST_EGALMIX(st) = clgetr ("egalmix")
+	ST_ABSORPTION(st) = clgetr ("absorption")
+
+	# Get the spatial density and luminosity function sampling parameters.
+	seed = clktime (long (0))
+	sseed1 = clgetl ("sseed")
+	if (IS_INDEFL(sseed1))
+	    sseed = sseed + seed
+	else
+	    sseed = sseed1
+	ST_SSEED(st) = sseed
+	lseed1 = clgetl ("lseed")
+	if (IS_INDEFL(lseed1))
+	    lseed = lseed + seed + 1
+	else
+	    lseed = lseed1
+	ST_LSEED(st) = lseed
+	ST_NSSAMPLE(st) = clgeti ("nssample")
+	ST_NLSAMPLE(st) = clgeti ("nlsample")
+	ST_SORDER(st) = clgeti ("sorder")
+	ST_LORDER(st) = clgeti ("lorder")
+	ST_RBINSIZE(st) = clgetr ("rbinsize")
+	ST_MBINSIZE(st) = clgetr ("mbinsize")
+	ST_DBINSIZE(st) = clgetr ("dbinsize")
+	ST_EBINSIZE(st) = clgetr ("ebinsize")
+	ST_PBINSIZE(st) = clgetr ("pbinsize")
+
+	x = NULL
+	y = NULL
+	mag = NULL
+	egal = NULL
+	axis = NULL
+	round = NULL
+	phi = NULL
+
+	#  Compute the spatial and luminosity functions.
+	call st_gmkspatial (sf, st, x, y, mag, egal, axis, round, phi)
+	call st_gmklum (lf, st, x, y, mag, egal, axis, round, phi)
+	call st_gmkmix (st, x, y, mag, egal, axis, round, phi)
+	call st_gmkaxis (st, x, y, mag, egal, axis, round, phi)
+	call st_gmkround (st, x, y, mag, egal, axis, round, phi)
+	call st_gmkphi (st, x, y, mag, egal, axis, round, phi)
+
+	# Plot the results.
+	if (clgetb ("interactive")) {
+	    call clgstr ("graphics", Memc[graphics], SZ_FNAME)
+	    if (Memc[graphics] != EOS) {
+	        gd = gopen (Memc[graphics], NEW_FILE, STDGRAPH)
+	        call st_gplots (sf, lf, gd, st, x, y, mag, egal, axis, round,
+		    phi)
+	        call gclose (gd)
+	    }
+	}
+
+	# Write the database.
+	dt = dtmap (Memc[galaxies], APPEND)
+	call st_dtginit (dt, st, Memc[galaxies], sseed, lseed)
+	call st_dtgwrite (dt, Memr[x], Memr[y], Memr[mag], Memi[egal],
+	    Memr[axis], Memr[round], Memr[phi], ST_NSTARS(st))
+	call dtunmap (dt)
+
+	# Free up memory.
+	if (x != NULL)
+	    call mfree (x, TY_REAL)
+	if (y != NULL)
+	    call mfree (y, TY_REAL)
+	if (mag != NULL)
+	    call mfree (mag, TY_REAL)
+	if (egal != NULL)
+	    call mfree (egal, TY_INT)
+	if (axis != NULL)
+	    call mfree (axis, TY_REAL)
+	if (round != NULL)
+	    call mfree (round, TY_REAL)
+	if (phi != NULL)
+	    call mfree (phi, TY_REAL)
+	call mfree (st, TY_STRUCT)
+
+	# Close files.
+	if (sf != NULL)
+	    call close (sf)
+	if (lf != NULL)
+	    call close (lf)
+
+	call sfree (sp)
+end
+
+
+# ST_GMKSPATIAL -- Compute the galactic spatial density function.
+
+procedure st_gmkspatial (sf, st, x, y, mag, egal, axis, round, phi)
+
+int	sf			# spatial density function file descriptor
+pointer	st			# pointer to the starlist strucuture
+pointer	x			# pointer to the x array
+pointer	y			# pointer to the y array
+pointer	mag			# pointer to the magnitude array
+pointer	egal			# pointer to the galaxy type array
+pointer	axis			# pointer to half-power diameter array
+pointer	round			# pointer to roundness array
+pointer	phi			# pointer to an array of position angles
+
+int	nsf
+pointer	r, rprob 
+int	st_gfetchxy()
+
+begin
+	# Check the sizes of the arrays.
+	call st_gmalloc (st, x, y, mag, egal, axis, round, phi)
+
+	# Compute the x and y values.
+	switch (ST_SPATIAL(st)) {
+	case ST_UNIFORM:
+	    call st_xyuniform (Memr[x], Memr[y], ST_NSTARS(st), ST_XMIN(st),
+	        ST_XMAX(st), ST_YMIN(st), ST_YMAX(st), ST_SSEED(st))
+
+	case ST_HUBBLE:
+	    call st_hbsample (Memr[x], Memr[y], ST_NSTARS(st), ST_CORE(st),
+	        ST_BASE(st), ST_XC(st), ST_YC(st), ST_XMIN(st), ST_XMAX(st),
+		ST_YMIN(st), ST_YMAX(st), ST_NSSAMPLE(st), ST_SORDER(st),
+		ST_SSEED(st)) 
+
+	case ST_SPFILE:
+	    if (sf == NULL) {
+		call printf ("The spatial density file is not open.\n")
+		call amovkr ((ST_XMIN(st) + ST_XMAX(st)) / 2.0, Memr[x],
+		    ST_NSTARS(st))
+		call amovkr ((ST_YMIN(st) + ST_YMAX(st)) / 2.0, Memr[y],
+		    ST_NSTARS(st))
+	    } else {
+	        nsf = st_gfetchxy (sf, r, rprob)
+	        if (nsf > 0) {
+	            call st_sfsample (Memr[r], Memr[rprob], nsf, Memr[x],
+		        Memr[y], ST_NSTARS(st), ST_NSSAMPLE(st), ST_SORDER(st),
+			ST_XC(st), ST_YC(st), ST_XMIN(st), ST_XMAX(st),
+			ST_YMIN(st), ST_YMAX(st), ST_SSEED(st))
+	        } else {
+		    call printf (
+		        "The spatial density function file is empty.\n")
+		    call amovkr ((ST_XMIN(st) + ST_XMAX(st)) / 2.0, Memr[x],
+		        ST_NSTARS(st))
+		    call amovkr ((ST_YMIN(st) + ST_YMAX(st)) / 2.0, Memr[y],
+		        ST_NSTARS(st))
+	        }
+	        call mfree (r, TY_REAL)
+	        call mfree (rprob, TY_REAL)
+	    }
+
+	default:
+	    call printf ("Unknown spatial density function.\n")
+	}
+end
+
+
+# ST_GMKLUM -- Compute the luminosity function and the diameter distribution
+# function.
+
+procedure st_gmklum (lf, st, x, y, mag, egal, axis, round, phi)
+
+int	lf		# luminsosity function file descriptor
+pointer	st		# pointer to starlist structure
+pointer	x		# pointer to the x array
+pointer	y		# pointer to the y array
+pointer	mag		# pointer to magnitude array
+pointer	egal		# pointer to the galaxy type array
+pointer	axis		# pointer to half-power diameter array
+pointer	round		# pointer to roundness array
+pointer	phi		# pointer to an array of position angles
+
+int	nlf
+pointer	m, mprob
+int	st_gfetchxy()
+
+begin
+	# Check the sizes of the arrays.
+	call st_gmalloc (st, x, y, mag, egal, axis, round, phi)
+
+	# Compute the luminosity  function.
+	switch (ST_LUMINOSITY(st)) {
+	case ST_UNIFORM:
+	    call st_maguniform (Memr[mag], ST_NSTARS(st), ST_MINMAG(st),
+	        ST_MAXMAG(st), ST_LSEED(st))
+
+	case ST_POWLAW:
+	    call st_power (Memr[mag], ST_NSTARS(st), ST_POWER(st),
+	        ST_MINMAG(st), ST_MAXMAG(st), ST_LSEED(st))
+
+	case ST_SCHECTER:
+	    call st_schecter (Memr[mag], ST_NSTARS(st), ST_ALPHA(st),
+	        ST_MSTAR(st), ST_MINMAG(st), ST_MAXMAG(st), ST_MZERO(st),
+		ST_NLSAMPLE(st), ST_LORDER(st), ST_LSEED(st))
+
+	case ST_LFFILE:
+	    if (lf == NULL) {
+		call printf ("The luminosity function file is not open.\n")
+		call amovkr ((ST_MINMAG(st) + ST_MAXMAG(st)) / 2.0, Memr[mag],
+		    ST_NSTARS(st))
+	    } else {
+	        nlf = st_gfetchxy (lf, m, mprob)
+	        if (nlf > 0) {
+	            call st_lfsample (Memr[m], Memr[mprob], nlf, Memr[mag],
+	                ST_NSTARS(st), ST_MINMAG(st), ST_MAXMAG(st),
+		        ST_NLSAMPLE(st), ST_LORDER(st), ST_LSEED(st))
+	        } else {
+		    call printf (
+		    "The luminosity function file is empty.\n")
+		    call amovkr ((ST_MINMAG(st) + ST_MAXMAG(st)) / 2.0,
+		        Memr[mag], ST_NSTARS(st))
+	        }
+	        call mfree (m, TY_REAL)
+	        call mfree (mprob, TY_REAL)
+	    }
+
+	default:
+	    call printf ("The luminosity function is unknown.\n")
+	}
+end
+
+
+# ST_GMKMIX -- Compute the percentage of elliptical versus spiral galaxies.
+
+procedure st_gmkmix (st, x, y, mag, egal, axis, round, phi)
+
+pointer	st		# pointer to the starlist structure
+pointer	x		# pointer to the x array
+pointer	y		# pointer to the y array
+pointer	mag		# pointer to magnitude array
+pointer	egal		# pointer to the galaxy type array
+pointer	axis		# pointer to half-power diameter array
+pointer	round		# pointer to roundness array
+pointer	phi		# pointer to an array of position angles
+
+begin
+	# Check the sizes of the arrays.
+	call st_gmalloc (st, x, y, mag, egal, axis, round, phi)
+
+	# Compute the elliptical / spiral galaxy mix.
+	call st_esmix (Memi[egal], ST_NSTARS(st), ST_EGALMIX(st), ST_SSEED(st))
+end
+
+
+# ST_GMKROUND -- Compute the roundness values.
+
+procedure st_gmkround (st, x, y, mag, egal, axis, round, phi)
+
+pointer	st		# pointer to the starlist structure
+pointer	x		# pointer to the x array
+pointer	y		# pointer to the y array
+pointer	mag		# pointer to magnitude array
+pointer	egal		# pointer to galaxy type array
+pointer	axis		# pointer to half-power diameter array
+pointer	round		# pointer to roundness array
+pointer	phi		# pointer to an array of position angles
+
+begin
+	# Check the sizes of the arrays.
+	call st_gmalloc (st, x, y, mag, egal, axis, round, phi)
+
+	# Compute the roundness values.
+	call st_round (Memi[egal], Memr[mag], Memr[round], ST_NSTARS(st),
+	    ST_AR(st), ST_ABSORPTION(st),  ST_SSEED(st))
+end
+
+
+# ST_GMKPHI -- Compute the position angles.
+
+procedure st_gmkphi (st, x, y, mag, egal, axis, round, phi)
+
+pointer	st		# pointer to the starlist structure
+pointer	x		# pointer to the x array
+pointer	y		# pointer to the y array
+pointer	mag		# pointer to mag array
+pointer	egal		# pointer to the  galaxy type array
+pointer	axis		# pointer to half-power diameter array
+pointer	round		# pointer to roundness array
+pointer	phi		# pointer to an array of position angles
+
+begin
+	# Check the sizes of the arrays.
+	call st_gmalloc (st, x, y, mag, egal, axis, round, phi)
+
+	# Compute the position angles.
+	call st_phi (Memr[phi], ST_NSTARS(st), ST_SSEED(st))
+end
+
+
+# ST_GMKAXIS -- Compute the semi-major axes.
+
+procedure st_gmkaxis (st, x, y, mag, egal, axis, round, phi)
+
+pointer	st		# pointer to the starlist structure
+pointer	x		# pointer to the x array
+pointer	y		# pointer to the y array
+pointer	mag		# pointer to magnitude array
+pointer	egal		# pointer to E or S galaxy array
+pointer	axis		# pointer to half-power diameter array
+pointer	round		# pointer to roundness array
+pointer	phi		# pointer to an array of position angles
+
+begin
+	# Check the sizes of the arrays.
+	call st_gmalloc (st, x, y, mag, egal, axis, round, phi)
+
+	# Make the effective diameters array.
+	switch (ST_LUMINOSITY(st)) {
+	case ST_POWLAW:
+	    call st_zdiameters (Memr[mag], Memi[egal], Memr[axis],
+		ST_NSTARS(st), ST_MINMAG(st), ST_MAXMAG(st), ST_Z(st),
+		ST_ERADIUS(st), ST_SRADIUS(st), ST_LSEED(st))
+	default:
+	    call st_diameters (Memr[mag], Memi[egal], Memr[axis],
+		ST_NSTARS(st), ST_MINMAG(st), ST_MAXMAG(st),
+		ST_ERADIUS(st), ST_SRADIUS(st), ST_LSEED(st))
+	}
+end
+
+
+# ST_GMALLOC -- Allocate array space for the gallist task.
+
+procedure st_gmalloc (st, x, y, mag, egal, axis, round, phi)
+
+pointer	st		# pointer to the starlist structure
+pointer	x		# pointer to the x array
+pointer	y		# pointer to the y array
+pointer	mag		# pointer to magnitude array
+pointer	egal		# pointer to galaxy type array
+pointer	axis		# pointer to half-power diameter array
+pointer	round		# pointer to roundness array
+pointer	phi		# pointer to an array of position angles
+
+begin
+	if (x == NULL)
+	    call malloc (x, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (x, ST_NSTARS(st), TY_REAL)
+	if (y == NULL)
+	    call malloc (y, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (y, ST_NSTARS(st), TY_REAL)
+	if (mag == NULL)
+	    call malloc (mag, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (mag, ST_NSTARS(st), TY_REAL)
+	if (egal == NULL)
+	    call malloc (egal, ST_NSTARS(st), TY_INT)
+	else
+	    call realloc (egal, ST_NSTARS(st), TY_INT)
+	if (axis == NULL)
+	    call malloc (axis, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (axis, ST_NSTARS(st), TY_REAL)
+	if (round == NULL)
+	    call malloc (round, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (round, ST_NSTARS(st), TY_REAL)
+	if (phi == NULL)
+	    call malloc (phi, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (phi, ST_NSTARS(st), TY_REAL)
+end
diff --git a/noao/artdata/lists/t_starlist.x b/noao/artdata/lists/t_starlist.x
new file mode 100644
index 00000000..f110a27c
--- /dev/null
+++ b/noao/artdata/lists/t_starlist.x
@@ -0,0 +1,303 @@
+include <fset.h>
+include "starlist.h"
+
+procedure t_starlist ()
+
+pointer	starlist		# pointer to the name of the output file
+pointer	graphics		# pointer to the graphics device name
+
+int	sf, lf
+long	seed, sseed, lseed
+long	sseed1, lseed1
+pointer	sp, str, x, y, mag, dt, st, gd
+
+bool	clgetb()
+int	clgeti(), clgwrd(), open()
+long	clgetl(), clktime()
+pointer	dtmap(), gopen()
+real	clgetr()
+
+begin
+	call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (starlist, SZ_FNAME, TY_CHAR)
+	call salloc (graphics, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Allocate the starlist structure.
+	call malloc (st, LEN_STSTRUCT, TY_STRUCT)
+	ST_TYPE(st) = ST_STARS
+
+	# Get the principal parameters.
+	call clgstr ("starlist", Memc[starlist], SZ_FNAME)
+	ST_NSTARS(st) = clgeti ("nstars")
+
+	# Get the charactersitics of the spatial density function.
+	ST_SPATIAL(st) = clgwrd ("spatial", ST_SPSTRING(st), SZ_FNAME, SPFUNCS)
+	ST_XMIN(st) = clgetr ("xmin")
+	ST_XMAX(st) = clgetr ("xmax")
+	ST_YMIN(st) = clgetr ("ymin")
+	ST_YMAX(st) = clgetr ("ymax")
+	ST_SFILE(st) = EOS
+	switch (ST_SPATIAL(st)) {
+	case ST_UNIFORM:
+	    sf = NULL
+	    ST_XC(st) = (ST_XMAX(st) + ST_XMIN(st)) / 2.0
+	    ST_YC(st) = (ST_YMAX(st) + ST_YMIN(st)) / 2.0
+	case ST_HUBBLE:
+	    sf = NULL
+	    ST_XC(st) = clgetr ("xcenter")
+	    if (IS_INDEFR(ST_XC(st)))
+	        ST_XC(st) = (ST_XMAX(st) + ST_XMIN(st)) / 2.0
+	    ST_YC(st) = clgetr ("ycenter")
+	    if (IS_INDEFR(ST_YC(st)))
+	        ST_YC(st) = (ST_YMAX(st) + ST_YMIN(st)) / 2.0
+	case ST_SPFILE:
+	    call clgstr ("sfile", ST_SFILE(st), SZ_FNAME)
+	    sf = open (ST_SFILE(st), READ_ONLY, TEXT_FILE)
+	    ST_XC(st) = clgetr ("xcenter")
+	    if (IS_INDEFR(ST_XC(st)))
+	        ST_XC(st) = (ST_XMAX(st) + ST_XMIN(st)) / 2.0
+	    ST_YC(st) = clgetr ("ycenter")
+	    if (IS_INDEFR(ST_YC(st)))
+	        ST_YC(st) = (ST_YMAX(st) + ST_YMIN(st)) / 2.0
+	}
+	ST_CORE(st) = clgetr ("core_radius")
+	ST_BASE(st) = clgetr ("base")
+
+
+	# Get the luminosity function parameters.
+	ST_LUMINOSITY(st) = clgwrd ("luminosity", ST_LFSTRING(st), SZ_FNAME,
+	    LUMFUNCS)
+	ST_MINMAG(st) = clgetr ("minmag")
+	ST_MAXMAG(st) = clgetr ("maxmag")
+	ST_LFILE(st) = EOS
+	switch (ST_LUMINOSITY(st)) {
+	case ST_UNIFORM, ST_SALPETER, ST_BANDS, ST_POWLAW:
+	    lf = NULL
+	case ST_LFFILE:
+	    call clgstr ("lfile", ST_LFILE(st), SZ_FNAME)
+	    lf = open (ST_LFILE(st), READ_ONLY, TEXT_FILE)
+	}
+	ST_POWER(st) = clgetr ("power")
+	ST_MZERO(st) = clgetr ("mzero")
+	ST_ALPHA(st) = clgetr ("alpha")
+	ST_BETA(st) = clgetr ("beta")
+	ST_DELTA(st) = clgetr ("delta")
+	ST_MSTAR(st) = clgetr ("mstar")
+
+	# Get the function sampling parameters.
+	seed = clktime (long (0))
+	sseed1 = clgetl ("sseed")
+	if (IS_INDEFL(sseed1))
+	    sseed = sseed + seed
+	else
+	    sseed = sseed1
+	ST_SSEED(st) = sseed
+	lseed1 = clgetl ("lseed")
+	if (IS_INDEFL(lseed1))
+	    lseed = lseed + seed + 1
+	else
+	    lseed = lseed1
+	ST_LSEED(st) = lseed
+	ST_NSSAMPLE(st) = clgeti ("nssample")
+	ST_NLSAMPLE(st) = clgeti ("nlsample")
+	ST_SORDER(st) = clgeti ("sorder")
+	ST_LORDER(st) = clgeti ("lorder")
+
+
+	#  Compute the initial spatial and luminosity functions.
+	x = NULL
+	y = NULL
+	mag = NULL
+	call st_mkspatial (sf, st, x, y, mag)
+	call st_mklum (lf, st, x, y, mag)
+
+	# Plot the results.
+	if (clgetb ("interactive")) {
+	    call clgstr ("graphics", Memc[graphics], SZ_FNAME)
+	    ST_RBINSIZE(st) = clgetr ("rbinsize")
+	    ST_MBINSIZE(st) = clgetr ("mbinsize")
+	    if (Memc[graphics] != EOS) {
+	        gd = gopen (Memc[graphics], NEW_FILE, STDGRAPH)
+	        call st_plots (sf, lf, gd, st, x, y, mag)
+	        call gclose (gd)
+	    }
+	}
+
+	# Write the database.
+	dt = dtmap (Memc[starlist], APPEND)
+	call st_dtinit (dt, st, Memc[starlist], sseed, lseed)
+	call st_dtwrite (dt, Memr[x], Memr[y], Memr[mag], ST_NSTARS(st))
+	call dtunmap (dt)
+
+	# Free up memory.
+	if (x != NULL)
+	    call mfree (x, TY_REAL)
+	if (y != NULL)
+	    call mfree (y, TY_REAL)
+	if (mag != NULL)
+	    call mfree (mag, TY_REAL)
+	call mfree (st, TY_STRUCT)
+
+	# Close files.
+	if (sf != NULL)
+	    call close (sf)
+	if (lf != NULL)
+	    call close (lf)
+
+	call sfree (sp)
+end
+
+
+# ST_MKSPATIAL -- Compute the spatial density function.
+
+procedure st_mkspatial (sf, st, x, y, mag)
+
+int	sf			# spatial function file descriptor
+pointer	st			# pointer to the starlist structure
+pointer	x			# pointer to the x coordinate array
+pointer	y			# pointer to the y coordinate array
+pointer	mag			# pointer to the magnitude array
+
+int	nsf
+pointer	r, rprob 
+int	st_gfetchxy()
+
+begin
+	# Dynamically reallocate the x, y and magnitude arrays.
+	if (x == NULL)
+	    call malloc (x, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (x, ST_NSTARS(st), TY_REAL)
+	if (y == NULL)
+	    call malloc (y, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (y, ST_NSTARS(st), TY_REAL)
+	if (mag == NULL)
+	    call malloc (mag, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (mag, ST_NSTARS(st), TY_REAL)
+
+	# Compute the x and y values.
+	switch (ST_SPATIAL(st)) {
+	case ST_UNIFORM:
+	    call st_xyuniform (Memr[x], Memr[y], ST_NSTARS(st), ST_XMIN(st),
+	        ST_XMAX(st), ST_YMIN(st), ST_YMAX(st), ST_SSEED(st))
+
+	case ST_HUBBLE:
+	    call st_hbsample (Memr[x], Memr[y], ST_NSTARS(st), ST_CORE(st),
+	        ST_BASE(st), ST_XC(st), ST_YC(st), ST_XMIN(st), ST_XMAX(st),
+		ST_YMIN(st), ST_YMAX(st), ST_NSSAMPLE(st), ST_SORDER(st),
+		ST_SSEED(st)) 
+
+	case ST_SPFILE:
+	    if (sf == NULL) {
+		call printf ("The spatial density file is not open.\n")
+		call amovkr ((ST_XMIN(st) + ST_XMAX(st)) / 2.0, Memr[x],
+		    ST_NSTARS(st))
+		call amovkr ((ST_YMIN(st) + ST_YMAX(st)) / 2.0, Memr[y],
+		    ST_NSTARS(st))
+	    } else {
+	        nsf = st_gfetchxy (sf, r, rprob)
+	        if (nsf > 0) {
+	            call st_sfsample (Memr[r], Memr[rprob], nsf, Memr[x],
+		        Memr[y], ST_NSTARS(st), ST_NSSAMPLE(st), ST_SORDER(st),
+			ST_XC(st), ST_YC(st), ST_XMIN(st), ST_XMAX(st),
+			ST_YMIN(st), ST_YMAX(st), ST_SSEED(st))
+	        } else {
+		    call printf (
+		        "The spatial density function file is empty.\n")
+		    call amovkr ((ST_XMIN(st) + ST_XMAX(st)) / 2.0, Memr[x],
+		        ST_NSTARS(st))
+		    call amovkr ((ST_YMIN(st) + ST_YMAX(st)) / 2.0, Memr[y],
+		        ST_NSTARS(st))
+	        }
+	        call mfree (r, TY_REAL)
+	        call mfree (rprob, TY_REAL)
+	    }
+
+	default:
+	    call printf ("Unknown spatial density function.\n")
+	}
+end
+
+
+# ST_MKLUM -- Compute the luminosity function.
+
+procedure st_mklum (lf, st, x, y, mag)
+
+int	lf		# luminsosity function file descriptor
+pointer	st		# pointer to starlist structure
+pointer	x		# pointer to the x coordinate array
+pointer	y		# pointer to the y coordinate array
+pointer	mag		# pointer to the magnitude array
+
+int	nlf
+pointer	m, mprob
+int	st_gfetchxy()
+
+begin
+	# Dynamically reallocate the array space.
+	if (x == NULL)
+	    call malloc (x, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (x, ST_NSTARS(st), TY_REAL)
+	if (y == NULL)
+	    call malloc (y, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (y, ST_NSTARS(st), TY_REAL)
+	if (mag == NULL)
+	    call malloc (mag, ST_NSTARS(st), TY_REAL)
+	else
+	    call realloc (mag, ST_NSTARS(st), TY_REAL)
+
+	# Compute the magnitudes.
+	switch (ST_LUMINOSITY(st)) {
+	case ST_UNIFORM:
+	    call st_maguniform (Memr[mag], ST_NSTARS(st), ST_MINMAG(st),
+	        ST_MAXMAG(st), ST_LSEED(st))
+
+	case ST_POWLAW:
+	    call st_power (Memr[mag], ST_NSTARS(st), ST_POWER(st),
+	        ST_MINMAG(st), ST_MAXMAG(st), ST_LSEED(st))
+
+	case ST_SALPETER:
+	    call st_salpeter (Memr[mag], ST_NSTARS(st),
+	        ST_MINMAG(st), ST_MAXMAG(st), ST_MZERO(st),
+		ST_NLSAMPLE(st), ST_LORDER(st), ST_LSEED(st))
+
+	case ST_BANDS:
+	    call st_bands (Memr[mag], ST_NSTARS(st), ST_ALPHA(st),
+	        ST_BETA(st), ST_DELTA(st), ST_MSTAR(st),
+	        ST_MINMAG(st), ST_MAXMAG(st), ST_MZERO(st),
+		ST_NLSAMPLE(st), ST_LORDER(st), ST_LSEED(st))
+
+	case ST_LFFILE:
+	    if (lf == NULL) {
+		call printf ("The luminosity function file is not open.\n")
+		call amovkr ((ST_MINMAG(st) + ST_MAXMAG(st)) / 2.0, Memr[mag],
+		    ST_NSTARS(st))
+	    } else {
+	        nlf = st_gfetchxy (lf, m, mprob)
+	        if (nlf > 0) {
+	            call st_lfsample (Memr[m], Memr[mprob], nlf, Memr[mag],
+	                ST_NSTARS(st),
+			ST_MINMAG(st), ST_MAXMAG(st),
+		        ST_NLSAMPLE(st), ST_LORDER(st), ST_LSEED(st))
+	        } else {
+		    call printf (
+		    "The luminosity function file is empty.\n")
+		    call amovkr ((ST_MINMAG(st) + ST_MAXMAG(st)) / 2.0,
+		        Memr[mag], ST_NSTARS(st))
+	        }
+	        call mfree (m, TY_REAL)
+	        call mfree (mprob, TY_REAL)
+	    }
+
+	default:
+	    call printf ("The luminosity function is unknown.\n")
+	}
+end