From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/artdata/lists/gallist.key | 63 +++ noao/artdata/lists/mkpkg | 18 + noao/artdata/lists/starlist.h | 136 +++++ noao/artdata/lists/starlist.key | 47 ++ noao/artdata/lists/stcolon.x | 835 +++++++++++++++++++++++++++++ noao/artdata/lists/stdbio.x | 350 +++++++++++++ noao/artdata/lists/stlum.x | 342 ++++++++++++ noao/artdata/lists/stmix.x | 144 +++++ noao/artdata/lists/stplot.x | 1102 +++++++++++++++++++++++++++++++++++++++ noao/artdata/lists/stshow.x | 142 +++++ noao/artdata/lists/stspatial.x | 264 ++++++++++ noao/artdata/lists/t_gallist.x | 453 ++++++++++++++++ noao/artdata/lists/t_starlist.x | 303 +++++++++++ 13 files changed, 4199 insertions(+) create mode 100644 noao/artdata/lists/gallist.key create mode 100644 noao/artdata/lists/mkpkg create mode 100644 noao/artdata/lists/starlist.h create mode 100644 noao/artdata/lists/starlist.key create mode 100644 noao/artdata/lists/stcolon.x create mode 100644 noao/artdata/lists/stdbio.x create mode 100644 noao/artdata/lists/stlum.x create mode 100644 noao/artdata/lists/stmix.x create mode 100644 noao/artdata/lists/stplot.x create mode 100644 noao/artdata/lists/stshow.x create mode 100644 noao/artdata/lists/stspatial.x create mode 100644 noao/artdata/lists/t_gallist.x create mode 100644 noao/artdata/lists/t_starlist.x (limited to 'noao/artdata/lists') 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 + stdbio.x starlist.h + stlum.x + stmix.x starlist.h + stplot.x starlist.h + stshow.x starlist.h + stspatial.x + t_gallist.x starlist.h + t_starlist.x starlist.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 +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 +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 +include +include +include + + +# 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 +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 +include +include +include +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 +include +include + + +# 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 +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 +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 -- cgit