1 files changed, 264 insertions, 0 deletions
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