Initial commit

1 files changed, 342 insertions, 0 deletions

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