diff options
Diffstat (limited to 'noao/digiphot/daophot/daoedit')
| -rw-r--r-- | noao/digiphot/daophot/daoedit/daoedit.h | 123 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/daoedit.key | 26 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/dpecolon.x | 441 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/dpeconfirm.x | 433 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/dpemark.x | 734 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/dpeomark.x | 68 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/dperplot.x | 172 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/dperprofile.x | 590 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/idaoedit.key | 22 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/mkpkg | 16 | ||||
| -rw-r--r-- | noao/digiphot/daophot/daoedit/t_daoedit.x | 317 |
11 files changed, 2942 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/daoedit/daoedit.h b/noao/digiphot/daophot/daoedit/daoedit.h new file mode 100644 index 00000000..c373f7c8 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/daoedit.h @@ -0,0 +1,123 @@ +# The definitions file for the DAOEDIT task +# At some point this should become paert of daophotdef.h + +define PSET_LIST "|datapars|centerpars|fitskypars|photpars|daopars|\ +findpars|" + +define PSET_DATAPARS 1 +define PSET_CENTERPARS 2 +define PSET_FITSKYPARS 3 +define PSET_PHOTPARS 4 +define PSET_DAOPARS 5 +define PSET_FINDPARS 6 + +define CMD_LIST "|lpar|epar|unlearn|scale|fwhmpsf|emission|sigma|datamin|\ +datamax|noise|ccdread|gain|readnoise|epadu|exposure|airmass|filter|obstime|\ +itime|xairmass|ifilter|otime|calgorithm|cbox|cthreshold|maxshift|minsnratio|\ +cmaxiter|clean|rclean|rclip|kclean|mkcenter|salgorithm|annulus|\ +dannulus|skyvalue|smaxiter|sloclip|shiclip|snreject|sloreject|shireject|\ +khist|binsize|smooth|rgrow|mksky|weighting|apertures|zmag|mkapert|\ +function|varorder|nclean|saturated|matchrad|psfrad|fitrad|recenter|fitsky|\ +groupsky|sannulus|wsannulus|flaterr|proferr|maxiter|clipexp|cliprange|\ +critoverlap|maxnstar|maxgroup|threshold|nsigma|ratio|theta|sharplo|sharphi|\ +roundlo|roundhi|mkdetections|" + +define CMD_LPAR 1 +define CMD_EPAR 2 +define CMD_UNLEARN 3 + +define CMD_SCALE 4 +define CMD_FWHMPSF 5 +define CMD_EMISSION 6 +define CMD_SIGMA 7 +define CMD_DATAMIN 8 +define CMD_DATAMAX 9 +define CMD_NOISE 10 +define CMD_CCDREAD 11 +define CMD_GAIN 12 +define CMD_READNOISE 13 +define CMD_EPADU 14 +define CMD_EXPOSURE 15 +define CMD_AIRMASS 16 +define CMD_FILTER 17 +define CMD_OBSTIME 18 +define CMD_ITIME 19 +define CMD_XAIRMASS 20 +define CMD_IFILTER 21 +define CMD_OTIME 22 + +define CMD_CALGORITHM 23 +define CMD_CBOX 24 +define CMD_CTHRESHOLD 25 +define CMD_MAXSHIFT 26 +define CMD_MINSNRATIO 27 +define CMD_CMAXITER 28 +define CMD_CLEAN 29 +define CMD_RCLEAN 30 +define CMD_RCLIP 31 +define CMD_KCLEAN 32 +define CMD_MKCENTER 33 + +define CMD_SALGORITHM 34 +define CMD_ANNULUS 35 +define CMD_DANNULUS 36 +define CMD_SKYVALUE 37 +define CMD_SMAXITER 38 +define CMD_SLOCLIP 39 +define CMD_SHICLIP 40 +define CMD_SNREJECT 41 +define CMD_SLOREJECT 42 +define CMD_SHIREJECT 43 +define CMD_KHIST 44 +define CMD_BINSIZE 45 +define CMD_SMOOTH 46 +define CMD_RGROW 47 +define CMD_MKSKY 48 + +define CMD_WEIGHTING 49 +define CMD_APERTURES 50 +define CMD_ZMAG 51 +define CMD_MKAPERT 52 + +define CMD_FUNCTION 53 +define CMD_VARORDER 54 +define CMD_NCLEAN 55 +define CMD_SATURATED 56 +define CMD_MATCHRAD 57 +define CMD_PSFRAD 58 +define CMD_FITRAD 59 +define CMD_RECENTER 60 +define CMD_FITSKY 61 +define CMD_GROUPSKY 62 +define CMD_SANNULUS 63 +define CMD_WSANNULUS 64 +define CMD_FLATERR 65 +define CMD_PROFERR 66 +define CMD_MAXITER 67 +define CMD_CLIPEXP 68 +define CMD_CLIPRANGE 69 +define CMD_CRITOVERLAP 70 +define CMD_MAXNSTAR 71 +define CMD_MAXGROUP 72 + +define CMD_THRESHOLD 73 +define CMD_NSIGMA 74 +define CMD_RATIO 75 +define CMD_THETA 76 +define CMD_SHARPLO 77 +define CMD_SHARPHI 78 +define CMD_ROUNDLO 79 +define CMD_ROUNDHI 80 +define CMD_MKDETECTIONS 81 + +# define the different WCSS for the radial profile plots + +define WCS_XPIX 1 +define WCS_XSCALE 2 + +define WCS_YCOUNT 1 +define WCS_YNORM 2 + +# miscellaneous definitions + +define MAX_NAPERTS 100 diff --git a/noao/digiphot/daophot/daoedit/daoedit.key b/noao/digiphot/daophot/daoedit/daoedit.key new file mode 100644 index 00000000..bd2f6b95 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/daoedit.key @@ -0,0 +1,26 @@ + Daoedit Commands + +? Print help +: Colon commands +i Interactive parameter setup menu +a Estimate center, sky, skysigma, fwhmpsf and magnitude of a star +r Plot the radial profile of an object and its integral +g Toggle between image and graphics cursor +x Toggle radial profile plot between pixel and scale units +y Toggle radial profile plot between counts and normal units +q Quit task + + Colon Commands + +:lparam/eparam/unlearn pset List/edit/unlearn the named pset +:parameter [value] List or set an individual pset parameter + + + Named Parameter Sets + +datapars The data dependent parameters +findpars The daofind task object detection parameters +centerpars The phot task centering algorithm parameters +fitskypars The phot task sky fitting algorithm parameters +photpars The phot task photometry algroithm parameters +daopars The psf fitting algorithm parameters diff --git a/noao/digiphot/daophot/daoedit/dpecolon.x b/noao/digiphot/daophot/daoedit/dpecolon.x new file mode 100644 index 00000000..36bccdf9 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/dpecolon.x @@ -0,0 +1,441 @@ +include <error.h> +include "daoedit.h" + +# DP_ECOLON -- Show/set the DAOPHOT algorithm parameters. + +procedure dp_ecolon (cmdstr, gd, redraw) + +char cmdstr[ARB] # input colon command +int gd # pointer to the graphics stream +int redraw # redraw the radial profile plot + +bool bval +int cmd, pset, ival +pointer sp, incmd, param, value +real rval +bool clgetb() +int strdic(), nscan(), clgeti() +real clgetr() +string datapars "datapars" +string centerpars "centerpars" +string fitskypars "fitskypars" +string photpars "photpars" +string daopars "daopars" +string findpars "findpars" + +begin + call smark (sp) + call salloc (incmd, SZ_LINE, TY_CHAR) + call salloc (param, SZ_LINE, TY_CHAR) + call salloc (value, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[incmd], SZ_LINE) + if (Memc[incmd] == EOS) { + call sfree (sp) + return + } + + cmd = strdic (Memc[incmd], Memc[incmd], SZ_LINE, CMD_LIST) + switch (cmd) { + case CMD_LPAR: + call gargwrd (Memc[incmd], SZ_LINE) + pset = strdic (Memc[incmd], Memc[incmd], SZ_LINE, PSET_LIST) + if (pset != 0) + call gdeactivate (gd, 0) + switch (pset) { + case PSET_DATAPARS: + call clcmdw ("lparam datapars") + case PSET_CENTERPARS: + call clcmdw ("lparam centerpars") + case PSET_FITSKYPARS: + call clcmdw ("lparam fitskypars") + case PSET_PHOTPARS: + call clcmdw ("lparam photpars") + case PSET_DAOPARS: + call clcmdw ("lparam daopars") + case PSET_FINDPARS: + call clcmdw ("lparam findpars") + default: + call printf ("Unknown parameter set\7\n") + } + + case CMD_EPAR: + call gargwrd (Memc[incmd], SZ_LINE) + pset = strdic (Memc[incmd], Memc[incmd], SZ_LINE, PSET_LIST) + if (pset != 0) + call gdeactivate (gd, 0) + switch (pset) { + case PSET_DATAPARS: + call clcmdw ("eparam datapars") + case PSET_CENTERPARS: + call clcmdw ("eparam centerpars") + case PSET_FITSKYPARS: + call clcmdw ("eparam fitskypars") + case PSET_PHOTPARS: + call clcmdw ("eparam photpars") + case PSET_DAOPARS: + call clcmdw ("eparam daopars") + case PSET_FINDPARS: + call clcmdw ("eparam findpars") + default: + call printf ("Unknown parameter set\7\n") + } + if (pset != 0) + redraw = YES + + case CMD_UNLEARN: + call gargwrd (Memc[incmd], SZ_LINE) + pset = strdic (Memc[incmd], Memc[incmd], SZ_LINE, PSET_LIST) + switch (pset) { + case PSET_DATAPARS: + call clcmdw ("unlearn datapars") + case PSET_CENTERPARS: + call clcmdw ("unlearn centerpars") + case PSET_FITSKYPARS: + call clcmdw ("unlearn fitskypars") + case PSET_PHOTPARS: + call clcmdw ("unlearn photpars") + case PSET_DAOPARS: + call clcmdw ("unlearn daopars") + case PSET_FINDPARS: + call clcmdw ("unlearn findpars") + default: + call printf ("Unknown parameter set\7\n") + } + if (pset != 0) + redraw = YES + + case CMD_SCALE: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (datapars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else { + call clputr (Memc[param], rval) + redraw = YES + } + + case CMD_FWHMPSF, CMD_SIGMA, CMD_DATAMIN, CMD_DATAMAX, + CMD_READNOISE, CMD_EPADU, CMD_ITIME, CMD_XAIRMASS: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (datapars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else + call clputr (Memc[param], rval) + + case CMD_EMISSION: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (datapars) + call pargstr (Memc[incmd]) + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (Memc[incmd]) + call pargb (clgetb (Memc[param])) + } else + call clputb (Memc[param], bval) + + case CMD_NOISE: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (datapars) + call pargstr (Memc[incmd]) + call clgstr (Memc[param], Memc[value], SZ_LINE) + call printf ("%s = %s\n") + call pargstr (Memc[incmd]) + call pargstr (Memc[value]) + + case CMD_CCDREAD, CMD_GAIN, CMD_EXPOSURE, CMD_AIRMASS, CMD_FILTER, + CMD_OBSTIME, CMD_IFILTER, CMD_OTIME: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (datapars) + call pargstr (Memc[incmd]) + call gargwrd (Memc[value], SZ_LINE) + if (nscan() == 1) { + call clgstr (Memc[param], Memc[value], SZ_LINE) + call printf ("%s = \"%s\"\n") + call pargstr (Memc[incmd]) + call pargstr (Memc[value]) + } else + call clpstr (Memc[param], Memc[value]) + + case CMD_CALGORITHM: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (centerpars) + call pargstr (Memc[incmd]) + call gargwrd (Memc[value], SZ_LINE) + if (nscan() == 1) { + call clgstr (Memc[param], Memc[value], SZ_LINE) + call printf ("%s = \"%s\"\n") + call pargstr (Memc[incmd]) + call pargstr (Memc[value]) + } else { + call clpstr (Memc[param], Memc[value]) + } + + case CMD_CBOX: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (centerpars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else { + call clputr (Memc[param], rval) + redraw = YES + } + + case CMD_CTHRESHOLD, CMD_MAXSHIFT, CMD_MINSNRATIO, CMD_RCLEAN, + CMD_RCLIP, CMD_KCLEAN: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (centerpars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else + call clputr (Memc[param], rval) + + case CMD_CMAXITER: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (centerpars) + call pargstr (Memc[incmd]) + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (Memc[incmd]) + call pargi (clgeti (Memc[param])) + } else + call clputi (Memc[param], ival) + + case CMD_CLEAN, CMD_MKCENTER: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (centerpars) + call pargstr (Memc[incmd]) + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (Memc[incmd]) + call pargb (clgetb (Memc[param])) + } else + call clputb (Memc[param], bval) + + case CMD_SALGORITHM: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (fitskypars) + call pargstr (Memc[incmd]) + call gargwrd (Memc[value], SZ_LINE) + if (nscan() == 1) { + call clgstr (Memc[param], Memc[value], SZ_LINE) + call printf ("%s = \"%s\"\n") + call pargstr (Memc[incmd]) + call pargstr (Memc[value]) + } else { + call clpstr (Memc[param], Memc[value]) + call clgstr (Memc[param], Memc[value], SZ_LINE) + } + + case CMD_ANNULUS, CMD_DANNULUS: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (fitskypars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else { + call clputr (Memc[param], rval) + redraw = YES + } + + case CMD_SLOCLIP, CMD_SHICLIP, CMD_SLOREJECT, CMD_SHIREJECT, + CMD_KHIST, CMD_BINSIZE, CMD_SKYVALUE, CMD_RGROW: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (fitskypars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else + call clputr (Memc[param], rval) + + case CMD_SMAXITER, CMD_SNREJECT: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (fitskypars) + call pargstr (Memc[incmd]) + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (Memc[incmd]) + call pargi (clgeti (Memc[param])) + } else + call clputi (Memc[param], ival) + + case CMD_SMOOTH, CMD_MKSKY: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (fitskypars) + call pargstr (Memc[incmd]) + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (Memc[incmd]) + call pargb (clgetb (Memc[param])) + } else + call clputb (Memc[param], bval) + + case CMD_WEIGHTING: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (photpars) + call pargstr (Memc[incmd]) + call clgstr (Memc[param], Memc[value], SZ_LINE) + call printf ("%s = %s\n") + call pargstr (Memc[incmd]) + call pargstr (Memc[value]) + + case CMD_APERTURES: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (photpars) + call pargstr (Memc[incmd]) + call gargwrd (Memc[value], SZ_LINE) + if (nscan() == 1) { + call clgstr (Memc[param], Memc[value], SZ_LINE) + call printf ("%s = \"%s\"\n") + call pargstr (Memc[incmd]) + call pargstr (Memc[value]) + } else { + call clpstr (Memc[param], Memc[value]) + redraw = YES + } + + case CMD_ZMAG: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (photpars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else { + call clputr (Memc[param], rval) + } + + case CMD_MKAPERT: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (photpars) + call pargstr (Memc[incmd]) + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (Memc[incmd]) + call pargb (clgetb (Memc[param])) + } else + call clputb (Memc[param], bval) + + case CMD_FUNCTION: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (daopars) + call pargstr (Memc[incmd]) + call clgstr (Memc[param], Memc[value], SZ_LINE) + call printf ("%s = %s\n") + call pargstr (Memc[incmd]) + call pargstr (Memc[value]) + + case CMD_RECENTER, CMD_FITSKY, CMD_GROUPSKY, CMD_SATURATED: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (daopars) + call pargstr (Memc[incmd]) + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (Memc[incmd]) + call pargb (clgetb (Memc[param])) + } else + call clputb (Memc[param], bval) + + case CMD_PSFRAD, CMD_FITRAD: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (daopars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else + call clputr (Memc[param], rval) + + case CMD_MATCHRAD, CMD_SANNULUS, CMD_WSANNULUS, CMD_FLATERR, + CMD_PROFERR, CMD_CRITOVERLAP, CMD_CLIPRANGE: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (daopars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else + call clputr (Memc[param], rval) + + case CMD_VARORDER, CMD_NCLEAN, CMD_MAXITER, CMD_MAXGROUP, + CMD_MAXNSTAR, CMD_CLIPEXP: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (daopars) + call pargstr (Memc[incmd]) + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (Memc[incmd]) + call pargi (clgeti (Memc[param])) + } else + call clputi (Memc[param], ival) + + case CMD_THRESHOLD, CMD_NSIGMA, CMD_RATIO, CMD_THETA, CMD_SHARPLO, + CMD_SHARPHI, CMD_ROUNDLO, CMD_ROUNDHI: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (findpars) + call pargstr (Memc[incmd]) + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (Memc[incmd]) + call pargr (clgetr (Memc[param])) + } else + call clputr (Memc[param], rval) + + case CMD_MKDETECTIONS: + call sprintf (Memc[param], SZ_LINE, "%s.%s") + call pargstr (findpars) + call pargstr (Memc[incmd]) + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (Memc[incmd]) + call pargb (clgetb (Memc[param])) + } else + call clputb (Memc[param], bval) + + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daoedit/dpeconfirm.x b/noao/digiphot/daophot/daoedit/dpeconfirm.x new file mode 100644 index 00000000..31aa359f --- /dev/null +++ b/noao/digiphot/daophot/daoedit/dpeconfirm.x @@ -0,0 +1,433 @@ +include "daoedit.h" + +# DP_CBANNER -- Pritnt the confirm banner. + +procedure dp_cbanner () + +begin + call printf ("\nVERIFY THE NEW VALUES\n") +end + + +# DP_CFWHMPSF -- Confirm the new value of fwhmpsf. + +procedure dp_cfwhmpsf() + +real fwhmpsf, scale, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + fwhmpsf = clgetr ("datapars.fwhmpsf") + + # Confirm the fwhmpsf. + call printf ("FWHM of features (%g scale units) (CR or value): ") + call pargr (fwhmpsf) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan() == 1) + fwhmpsf = rval + } + call printf ("\tNew FWHM of features: %g scale units %g pixels\n") + call pargr (fwhmpsf) + call pargr (scale * fwhmpsf) + + # Store the new fwhmpsf. + call clputr ("datapars.fwhmpsf", fwhmpsf) +end + + +# DP_CSIGMA -- Confirm the sigma parameters. + +procedure dp_csigma() + +real sigma, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current value. + sigma = clgetr ("datapars.sigma") + + # Confirm the sky sigma. + call printf ( + "Standard deviation of background (%g counts) (CR or value): ") + call pargr (sigma) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + sigma = rval + } + call printf ("\tNew standard deviation of background: %g counts\n") + call pargr (sigma) + + # Store the new sky sigma. + call clputr ("datapars.sigma", sigma) +end + + +# DP_CDMIN -- Confirm the good data minimum value. + +procedure dp_cdmin () + +real datamin, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current value. + datamin = clgetr ("datapars.datamin") + + # Confirm the new minimum good data value. + call printf ( + "Minimum good data value (%g counts) (CR or value): ") + call pargr (datamin) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + datamin = rval + } + call printf ("\tNew good data minimum: %g counts\n") + call pargr (datamin) + + # Store the new good data minimum. + call clputr ("datapars.datamin", datamin) +end + + +# DP_CDMAX -- Confirm the good data maximum value. + +procedure dp_cdmax () + +real datamax, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current value. + datamax = clgetr ("datapars.datamax") + + # Confirm the new maximum good data value. + call printf ( + "Maximum good data value (%g counts) (CR or value): ") + call pargr (datamax) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + datamax = rval + } + call printf ("\tNew good data maximum: %g counts\n") + call pargr (datamax) + + # Store the new maximum good data value. + call clputr ("datapars.datamax", datamax) +end + + +# DP_CCBOX -- Confirm the cbox parameter. + +procedure dp_ccbox() + +real scale, cbox, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + cbox = clgetr ("centerpars.cbox") + + # Confirm the centering box value. + call printf ("Centering box width (%g scale units) (CR or value): ") + call pargr (cbox) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + cbox = rval + } + call printf ("\tNew centering box width: %g scale units %g pixels\n") + call pargr (cbox) + call pargr (scale * cbox) + + # Store the new centering box. + call clputr ("centerpars.cbox", cbox) +end + + +# DP_CRCLEAN -- Confirm the cleaning radius. + +procedure dp_crclean() + +real scale, rclean, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + rclean = clgetr ("centerpars.rclean") + + # Confirm the cleaning radius.. + call printf ("Cleaning radius (%g scale units) (CR or value): ") + call pargr (rclean) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + rclean = rval + } + call printf ("\tNew cleaning radius: %g scale units %g pixels\n") + call pargr (rclean) + call pargr (scale * rclean) + + # Store the new cleaning radius. + call clputr ("centerpars.rclean", rclean) +end + + +# DP_CRCLIP -- Confirm the clipping radius. + +procedure dp_crclip() + +real scale, rclip, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + rclip = clgetr ("centerpars.rclip") + + # Confirm the cleaning radius.. + call printf ("Clipping radius (%g scale units) (CR or value): ") + call pargr (rclip) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + rclip = rval + } + call printf ("\tNew clipping radius: %g scale units %g pixels\n") + call pargr (rclip) + call pargr (scale * rclip) + + # Store the new clipping radius. + call clputr ("centerpars.rclip", rclip) +end + + +# DP_CANNULUS -- Confirm the inner radius of the sky annulus. + +procedure dp_cannulus() + +real scale, annulus, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + annulus = clgetr ("fitskypars.annulus") + + # Confirm the sky annulus. + call printf ( + "Inner radius of sky annulus (%g scale units) (CR or value): ") + call pargr (annulus) + call flush (STDOUT) + if (scan () != EOF) { + call gargr (rval) + if (nscan () == 1) + annulus = rval + } + call printf ( + "\tNew inner radius of sky annulus: %g scale units %g pixels\n") + call pargr (annulus) + call pargr (scale * annulus) + + # Store the new sky annulus. + call clputr ("fitskypars.annulus", annulus) +end + + +# DP_CDANNULUS -- Confirm the annulus width parameter. + +procedure dp_cdannulus() + +real scale, dannulus, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + dannulus = clgetr ("fitskypars.dannulus") + + # Confirm the sky annulus width. + call printf ( + "Width of the sky annulus (%g scale units) (CR or value): ") + call pargr (dannulus) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + dannulus = rval + } + call printf ( + "\tNew width of the sky annulus: %g scale units %g pixels\n") + call pargr (dannulus) + call pargr (scale * dannulus) + + # Save the new sky annulus width. + call clputr ("fitskypars.dannulus", dannulus) +end + + +# DP_CRGROW -- Confirm the region growing radius. + +procedure dp_crgrow() + +real scale, rgrow, rval +int scan(), nscan() +real clgetr() + +begin + # Mark the region growing radius. + scale = 1.0 / clgetr ("datapars.scale") + rgrow = clgetr ("fitskypars.rgrow") + + # Confirm the new region growing radius. + call printf ( + "Region growing radius (%g scale units) (CR or value): ") + call pargr (rgrow) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + rgrow = rval + } + call printf ( + "\tNew region growing radius: %g scale units %g pixels\n") + call pargr (rgrow) + call pargr (scale * rgrow) + + # Save the new region growing radius. + call clputr ("fitskypars.rgrow", rgrow) +end + + +# DP_CAPER -- Confirm the aperture string. + +procedure dp_caper() + +int i, naperts +pointer sp, apstr, newapstr, aperts +real scale +int scan(), nscan(), dp_gaperts() +real clgetr() + +begin + call smark (sp) + call salloc (apstr, SZ_LINE, TY_CHAR) + call salloc (newapstr, SZ_LINE, TY_CHAR) + call salloc (aperts, MAX_NAPERTS, TY_REAL) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + call clgstr ("photpars.apertures", Memc[apstr], SZ_LINE) + + # Confirm the aperture string. + call printf ( + "File/list of aperture radii (%s scale units) (CR or value): ") + call pargstr (Memc[apstr]) + call flush (STDOUT) + + # Get the new apertures. + if (scan() != EOF) { + call gargwrd (Memc[newapstr], SZ_LINE) + if (nscan () == 1) + call strcpy (Memc[newapstr], Memc[apstr], SZ_LINE) + } + + # Print the new apertures. + naperts = dp_gaperts (Memc[apstr], Memr[aperts], MAX_NAPERTS) + do i = 1, naperts { + call printf ("\tAperture radius %d: %g scale units %g pixels\n") + call pargi (i) + call pargr (Memr[aperts+i-1]) + call pargr (scale * Memr[aperts+i-1]) + } + + # Save the new aperture string. + call clpstr ("photpars.apertures", Memc[apstr]) + + call sfree (sp) +end + + +# DP_CPSFRAD -- Confirm the psf radius. + +procedure dp_cpsfrad() + +real scale, psfrad, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + psfrad = clgetr ("daopars.psfrad") + + # Confirm the new PSF radius. + call printf ("PSF radius (%g scale units) (CR or value): ") + call pargr (psfrad) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + psfrad = rval + } + call printf ("\tNew PSF radius: %g scale units %g pixels\n") + call pargr (psfrad) + call pargr (scale * psfrad) + + # Store the new PSF radius. + call clputr ("daopars.psfrad", psfrad) +end + + +# DP_CFITRAD -- Confirm the fitting radius. + +procedure dp_cfitrad () + +real scale, fitrad, rval +int scan(), nscan() +real clgetr() + +begin + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + fitrad = clgetr ("daopars.fitrad") + + # Confirm the new fitting radius. + call printf ("Fitting radius (%g scale units) (CR or value): ") + call pargr (fitrad) + call flush (STDOUT) + if (scan() != EOF) { + call gargr (rval) + if (nscan () == 1) + fitrad = rval + } + call printf ("\tNew fitting radius: %g scale units %g pixels\n") + call pargr (fitrad) + call pargr (scale * fitrad) + + # Store the new fitting radius. + call clputr ("daopars.fitrad", fitrad) +end diff --git a/noao/digiphot/daophot/daoedit/dpemark.x b/noao/digiphot/daophot/daoedit/dpemark.x new file mode 100644 index 00000000..de72fc2c --- /dev/null +++ b/noao/digiphot/daophot/daoedit/dpemark.x @@ -0,0 +1,734 @@ +include <ctype.h> +include "daoedit.h" + +# DP_MFWHMPSF -- Mark the fwhmpsf on the radial profile plot and confirm. + +procedure dp_mfwhmpsf (gd) + +pointer gd # pointer to the graphics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, fwhmpsf, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the x and y limits of the current plot. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current parameters. + scale = 1.0 / clgetr ("datapars.scale") + fwhmpsf = clgetr ("datapars.fwhmpsf") / 2.0 + + # Mark the FWHM of the PSF on the radial profile plot. + call printf ("Mark HWHM of the psf (%g pixels):") + call pargr (fwhmpsf * scale) + call gscur (gd, fwhmpsf * scale, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx <= 0.0 || wx > rmax) + ; + else + fwhmpsf = wx / scale + + # Store the new fwhmpsf. + call clputr ("datapars.fwhmpsf", 2.0 * fwhmpsf) + + call sfree (sp) +end + + +# DP_MSIGMA -- Mark the sky sigma on the radial profile plot and confirm. + +procedure dp_msigma (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, mean, sigma, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the range of the plot + call ggwind (gd, rmin, rmax, imin, imax) + + # Mark the mean sky on the radial profile plot. + call printf ("Mark mean sky background level:") + call gscur (gd, (rmin + rmax) / 2.0, imin) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wy < imin || wy > imax) + mean = imin + else + mean = wy + + # Get the current value. + sigma = clgetr ("datapars.sigma") + if (! IS_INDEFR (sigma)) + sigma = 3.0 * sigma + + # Mark the sky sigma on the radial profile plot. + call printf ("Mark 3 sigma sky level (%g counts):") + call pargr (sigma) + if (IS_INDEFR(sigma)) + call gscur (gd, (rmin + rmax) / 2.0, imax) + else + call gscur (gd, (rmin + rmax) / 2.0, mean + sigma) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wy < imin || wy > imax) + ; + else + sigma = abs (wy - mean) / 3.0 + + # Store the new sky sigma. + call clputr ("datapars.sigma", sigma) + + call sfree (sp) +end + + +# DP_MDMIN -- Mark the minimum good data value on the radial profile plot +# and confirm. + +procedure dp_mdmin (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, datamin, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the limits of the plot. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current value. + datamin = clgetr ("datapars.datamin") + + # Mark the threshold on the radial profile plot. + call printf ("Mark the minimum good data level (%g counts):") + call pargr (datamin) + + if (IS_INDEFR(datamin) || datamin < imin) + call gscur (gd, (rmin + rmax) / 2.0, imin) + else + call gscur (gd, (rmin + rmax) / 2.0, datamin) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wy < imin || wy > imax) + ; + else + datamin = wy + + # Store the new good data minimum. + call clputr ("datapars.datamin", datamin) + + call sfree (sp) +end + + +# DP_MDMAX -- Mark the maximum good data value on the radial profile plot +# and confirm. + +procedure dp_mdmax (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, datamax, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the limits of the plot. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current value. + datamax = clgetr ("datapars.datamax") + + # Mark the threshold on the radial profile plot. + call printf ("Mark the maximum good data level (%g counts):") + call pargr (datamax) + + if (IS_INDEFR(datamax) || datamax > imax) + call gscur (gd, (rmin + rmax) / 2.0, imax) + else + call gscur (gd, (rmin + rmax) / 2.0, datamax) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wy < imin || wy > imax) + ; + else + datamax = wy + + # Store the new maximum good data value. + call clputr ("datapars.datamax", datamax) + + call sfree (sp) +end + + +# DP_MCBOX -- Mark the centering aperture on the radial profile plot and +# confirm. + +procedure dp_mcbox (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, capert, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the x and y limits of the current plot. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + capert = clgetr ("centerpars.cbox") / 2.0 + + # Mark the centering aperture on the radial profile plot. + call printf ("Mark centering box half width (%g pixels):") + call pargr (capert * scale) + call gscur (gd, capert * scale, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx <= 0.0 || wx > rmax) + ; + else + capert = wx / scale + + # Store the new centering box. + call clputr ("centerpars.cbox", 2.0 * capert) + + call sfree (sp) +end + + +# DP_MRCLEAN -- Mark the cleaning radius on the radial profile plot and +# confirm. + +procedure dp_mrclean (gd) + +pointer gd # pointer to the graphics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, rclean, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Get the current plot window. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + rclean = clgetr ("centerpars.rclean") + + # Mark the cleaning radius on the plot. + call printf ( + "Mark the centering algorithm cleaning radius (%g pixels):") + call pargr (scale * rclean) + call gscur (gd, scale * rclean, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx <= 0.0 || wx > rmax) + ; + else + rclean = wx / scale + + # Store the new cleaning radius. + call clputr ("centerpars.rclean", rclean) + + call sfree (sp) +end + + +# DP_MRCLIP -- Mark the clipping radius on the radial profile plot and. +# confirm. + +procedure dp_mrclip (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, rclip, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Get the current plot window. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the clipping radius values. + scale = 1.0 / clgetr ("datapars.scale") + rclip = clgetr ("centerpars.rclip") + + # Mark clipping radius on the plot. + call printf ( + "Mark the centering algorithm clipping radius (%g pixels):") + call pargr (scale * rclip) + call gscur (gd, scale * rclip, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx <= 0.0 || wx > rmax) + ; + else + rclip = wx / scale + + # Store the new clipping radius. + call clputr ("centerpars.rclip", rclip) + + call sfree (sp) +end + + +# DP_MANNULUS -- Mark the sky annulus on the radial profile plot and confirm. + +procedure dp_mannulus (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, annulus, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Get the current plot window. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + annulus = clgetr ("fitskypars.annulus") + + # Mark the inner sky radius. + call printf ("Mark inner sky radius (%g pixels):") + call pargr (annulus * scale) + call gscur (gd, annulus * scale, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx < 0.0 || wx > rmax) + ; + else + annulus = wx / scale + + # Store the new sky annulus. + call clputr ("fitskypars.annulus", annulus) + + call sfree (sp) + +end + + +# DP_MDANNULUS -- Mark the sky annulus width on the radial profile plot and +# confirm. + +procedure dp_mdannulus (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, annulus, dannulus, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Get the current plot window. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + annulus = clgetr ("fitskypars.annulus") + dannulus = clgetr ("fitskypars.dannulus") + + # Mark the outer sky radius. + call printf ("Mark outer sky radius (%g pixels):") + call pargr (scale * (annulus + dannulus)) + call gscur (gd, scale * (annulus + dannulus), (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || (wx / scale < annulus) || wx > rmax) + ; + else + dannulus = (wx / scale - annulus) + + # Save the new sky annulus width. + call clputr ("fitskypars.dannulus", dannulus) + + call sfree (sp) +end + + +# DP_MRGROW -- Mark the regions growing radius the radial profile plot. + +procedure dp_mrgrow (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, rgrow, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Get the current plot window. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + rgrow = clgetr ("fitskypars.rgrow") + + # Mark the inner sky radius. + call printf ("Mark region growing radius (%g pixels):") + call pargr (rgrow * scale) + call gscur (gd, rgrow * scale, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx < 0.0 || wx > rmax) + ; + else + rgrow = wx / scale + + # Store the new sky annulus. + call clputr ("fitskypars.rgrow", rgrow) + + call sfree (sp) + +end + + +# DP_MAPER -- Mark the photometry apertures on the radial profile plot and +# confirm. + +procedure dp_maper (gd) + +pointer gd # pointer to the grapics stream + +int wcs, key, naperts +pointer sp, oapstr, aperts, tapstr, apstr, cmd +real rmin, rmax, imin, imax, scale, wx, wy +int dp_gaperts(), clgcur(), strlen() +real clgetr() + +begin + call smark (sp) + call salloc (oapstr, SZ_LINE, TY_CHAR) + call salloc (aperts, MAX_NAPERTS, TY_REAL) + call salloc (apstr, SZ_LINE, TY_CHAR) + call salloc (tapstr, SZ_LINE, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the current plot window. + call ggwind (gd, rmin, rmax, imin, imax) + + # Decode the apertures. + scale = 1.0 / clgetr ("datapars.scale") + call clgstr ("photpars.apertures", Memc[oapstr], SZ_LINE) + naperts = dp_gaperts (Memc[oapstr], Memr[aperts], MAX_NAPERTS) + + # Type prompt string. + call printf ("Mark apertures (%s pixels) [q=quit]:") + call pargstr (Memc[oapstr]) + call gscur (gd, Memr[aperts] * scale, (imin + imax) / 2.0) + + # Mark the apertures. + Memc[apstr] = EOS + Memc[tapstr] = EOS + while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + if (key == 'q') + break + if (wx <= 0.0 || wx > rmax) + next + call sprintf (Memc[apstr+strlen(Memc[apstr])], SZ_FNAME,"%.2f,") + call pargr (wx / scale) + call sprintf (Memc[tapstr+strlen(Memc[tapstr])], SZ_FNAME,"%.2f,") + call pargr (wx) + call printf ("Mark apertures (%s pixels) [q=quit]:") + call pargstr (Memc[tapstr]) + } + Memc[apstr+strlen(Memc[apstr])-1] = EOS + + # Save the new aperture string. + call clpstr ("photpars.apertures", Memc[apstr]) + + call sfree (sp) +end + + +# DP_MPSFRAD -- Mark the psf radius on the radial profile plot and confirm. + +procedure dp_mpsfrad (gd) + +pointer gd # pointer to the graphics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, psfrad, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the x and y limits of the current plot. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + psfrad = clgetr ("daopars.psfrad") + + # Mark the FWHM of the PSF on the radial profile plot. + call printf ("Mark the PSF radius (%g pixels):") + call pargr (psfrad * scale) + call gscur (gd, psfrad * scale, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx <= 0.0 || wx > rmax) + ; + else + psfrad = wx / scale + + # Store the new PSF radius. + call clputr ("daopars.psfrad", psfrad) + + call sfree (sp) +end + + +# DP_MFITRAD -- Mark the fitting radius on the radial profile plot and confirm. + +procedure dp_mfitrad (gd) + +pointer gd # pointer to the graphics stream + +int wcs, key, stat +pointer sp, cmd +real rmin, rmax, imin, imax, scale, fitrad, wx, wy +int clgcur() +real clgetr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Determine the x and y limits of the current plot. + call ggwind (gd, rmin, rmax, imin, imax) + + # Get the current values. + scale = 1.0 / clgetr ("datapars.scale") + fitrad = clgetr ("daopars.fitrad") + + # Mark the FWHM of the PSF on the radial profile plot. + call printf ("Mark the fitting radius (%g pixels):") + call pargr (fitrad * scale) + call gscur (gd, fitrad * scale, (imin + imax) / 2.0) + stat = clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) + if (stat == EOF || wx <= 0.0 || wx > rmax) + ; + else + fitrad = wx / scale + + # Store the new fitting radius. + call clputr ("daopars.fitrad", fitrad) + + call sfree (sp) +end + + +# DP_GAPERTS -- Decode the aperture string. + +int procedure dp_gaperts (str, aperts, max_naperts) + +char str[ARB] # aperture string +real aperts[ARB] # aperture array +int max_naperts # maximum number of apertures + +int naperts, ip, op, ndecode, nap +pointer sp, outstr +real apstart, apend, apstep +bool fp_equalr() +int dp_gctor() + +begin + call smark (sp) + call salloc (outstr, SZ_LINE, TY_CHAR) + + naperts = 0 + for (ip = 1; str[ip] != EOS && naperts < max_naperts;) { + + # Initialize. + apstart = 0.0 + apend = 0.0 + apstep = 0.0 + ndecode = 0 + + # Skip past white space and commas. + while (IS_WHITE(str[ip])) + ip = ip + 1 + if (str[ip] == ',') + ip = ip + 1 + + # Get the starting aperture number. + op = 1 + while (IS_DIGIT(str[ip]) || str[ip] == '.') { + Memc[outstr+op-1] = str[ip] + ip = ip + 1 + op = op + 1 + } + Memc[outstr+op-1] = EOS + + # Decode the starting aperture. + op = 1 + if (dp_gctor (Memc[outstr], op, apstart) > 0) { + apend = apstart + ndecode = 1 + } else + apstart = 0.0 + + # Skip past white space and commas. + while (IS_WHITE(str[ip])) + ip = ip + 1 + if (str[ip] == ',') + ip = ip + 1 + + # Search for the ending aperture. + if (str[ip] == ':') { + ip = ip + 1 + + # Get the ending aperture. + op = 1 + while (IS_DIGIT(str[ip]) || str[ip] == '.') { + Memc[outstr+op-1] = str[ip] + ip = ip + 1 + op = op + 1 + } + Memc[outstr+op-1] = EOS + + # Decode the ending aperture. + op = 1 + if (dp_gctor (Memc[outstr], op, apend) > 0) { + ndecode = 2 + apstep = apend - apstart + } + } + + # Skip past the white space. + while (IS_WHITE(str[ip])) + ip = ip + 1 + + # Skip past the commas. + if (str[ip] == ',') + ip = ip + 1 + + # Get the step size. + if (str[ip] == ':') { + ip = ip + 1 + + # Get the step size. + op = 1 + while (IS_DIGIT(str[ip]) || str[ip] == '.') { + Memc[outstr+op-1] = str[ip] + ip = ip + 1 + op = op + 1 + } + Memc[outstr+op-1] = EOS + + # Decode the step size. + op = 1 + if (dp_gctor (Memc[outstr], op, apstep) > 0) { + if (fp_equalr (apstep, 0.0)) + apstep = apend - apstart + else + ndecode = (apend - apstart) / apstep + 1 + if (ndecode < 0) { + ndecode = -ndecode + apstep = - apstep + } + } + } + + # Negative apertures are not permitted. + if (apstart <= 0.0 || apend <= 0.0) + break + + # Fill in the apertures. + if (ndecode == 0) { + ; + } else if (ndecode == 1) { + naperts = naperts + 1 + aperts[naperts] = apstart + } else if (ndecode == 2) { + naperts = naperts + 1 + aperts[naperts] = apstart + if (naperts >= max_naperts) + break + naperts = naperts + 1 + aperts[naperts] = apend + } else { + for (nap = 1; nap <= ndecode && naperts < max_naperts; + nap = nap + 1) { + naperts = naperts + 1 + aperts[naperts] = apstart + (nap - 1) * apstep + } + } + } + + call sfree (sp) + + return (naperts) +end + + +# DP_GCTOR -- Procedure to convert a character variable to a real number. +# This routine is just an interface routine to the IRAF procedure gctod. + +int procedure dp_gctor (str, ip, rval) + +char str[ARB] # string to be converted +int ip # pointer to the string +real rval # real value + +double dval +int nchars +int gctod() + +begin + nchars = gctod (str, ip, dval) + rval = dval + return (nchars) +end diff --git a/noao/digiphot/daophot/daoedit/dpeomark.x b/noao/digiphot/daophot/daoedit/dpeomark.x new file mode 100644 index 00000000..d635c31a --- /dev/null +++ b/noao/digiphot/daophot/daoedit/dpeomark.x @@ -0,0 +1,68 @@ +include <gset.h> + +# DP_EOMARK -- Procedure to mark center, fitsky and phot parameters on the +# display. + +procedure dp_eomark (id, xcenter, ycenter, srin, srout, apr, mkcenter, mksky, + mkapert) + +pointer id # pointer to image display stream +real xcenter # the center x coordinate +real ycenter # the center y coordinate +real srin # the inner radius of the sky annulus +real srout # the outer radius of the sky annulus +real apr # the maximum photometry aperture radius +int mkcenter # mark the computed center +int mksky # mark the sky annulus +int mkapert # mark the aperture(s) + +real rad +int marktype +int gstati() + +errchk greactivate, gdeactivate + +begin + if (id == NULL) + return + if (mkcenter == NO && mksky == NO && mkapert == NO) + return + iferr { + call greactivate (id, 0) + } then { + return + } + + # Save old mark type. + marktype = gstati (id, G_PMLTYPE) + + # Mark the center and shift on the display. + if (mkcenter == YES) { + call gseti (id, G_PMLTYPE, GL_SOLID) + call gmark (id, xcenter, xcenter, GM_PLUS, -2.0, -2.0) + } + + # Draw the sky annuli on the display. + if (mksky == YES) { + call gseti (id, G_PMLTYPE, GL_DASHED) + rad = 2.0 * srin + call gmark (id, xcenter, ycenter, GM_CIRCLE, -rad, -rad) + rad = 2.0 * srout + call gmark (id, xcenter, ycenter, GM_CIRCLE, -rad, -rad) + } + + # Draw the apertures on the display. + if (mkapert == YES) { + call gseti (id, G_PMLTYPE, GL_DASHED) + rad = 2.0 * apr + call gmark (id, xcenter, ycenter, GM_CIRCLE, -rad, -rad) + } + + # Restore the mark type. + call gseti (id, G_PMLTYPE, marktype) + + iferr { + call gdeactivate (id, 0) + } then + return +end diff --git a/noao/digiphot/daophot/daoedit/dperplot.x b/noao/digiphot/daophot/daoedit/dperplot.x new file mode 100644 index 00000000..95a78848 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/dperplot.x @@ -0,0 +1,172 @@ +include <gset.h> +include "daoedit.h" + +define FRACTION 0.10 + +# DP_ERPLOT -- Plot the radial profile. + +procedure dp_erplot (gd, title, xwcs, ywcs, radius, intensity, npts, + rcentroid, pmean, integral, nbins, rmin, rmax, iannulus, oannulus, + apradius, skyval, skysigma, rscale, pnorm) + +pointer gd # pointer to the graphics stream +char title[ARB] # the plot title +int xwcs # the x wcs type of the final plot +int ywcs # the y wcs type of the final plot +real radius[ARB] # the radius vector +real intensity[ARB] # the intensity vector +int npts # number of points in the profile +real rcentroid[ARB] # the radius centroid vector +real pmean[ARB] # the mean intensity vector +real integral[ARB] # the integral of the profile +int nbins # the number of bins +real rmin, rmax # min and max radius +real iannulus # the inner radius of the sky annulus +real oannulus # the outer radius of the sky annulus +real apradius # the aperture radius +real skyval # the sky value +real skysigma # the sigma of the sky value +real rscale # the image scale +real pnorm # the profile normalization factor + +int i, j +pointer sp, pxlabel, pylabel, sxlabel, sylabel +real r1, r2, rp1, rp2, rs1, rs2, i1, i2, ip1, ip2, is1, is2, dr, fraction + +begin + # Get space for the x and y labels. + call smark (sp) + call salloc (pxlabel, SZ_FNAME, TY_CHAR) + call salloc (pylabel, SZ_FNAME, TY_CHAR) + call salloc (sxlabel, SZ_FNAME, TY_CHAR) + call salloc (sylabel, SZ_FNAME, TY_CHAR) + + # Clear the plot. + call gclear (gd) + + # Determine the data range of the x and y axes. + + r1 = rmin - FRACTION * (rmax - rmin) + r2 = rmax + FRACTION * (rmax - rmin) + switch (xwcs) { + case WCS_XPIX: + rp1 = r1 + rp2 = r2 + rs1 = r1 / rscale + rs2 = r2 / rscale + case WCS_XSCALE: + rp1 = r1 / rscale + rp2 = r2 / rscale + rs1 = r1 + rs2 = r2 + default: + rp1 = r1 + rp2 = r2 + rs1 = r1 / rscale + rs2 = r2 / rscale + } + fraction = max (FRACTION, 5.0 * skysigma / pnorm) + i1 = -pnorm * fraction + i2 = pnorm * (1.0 + fraction) + switch (ywcs) { + case WCS_YNORM: + ip1 = i1 / pnorm + ip2 = i2 / pnorm + is1 = i1 + skyval + is2 = i2 + skyval + case WCS_YCOUNT: + ip1 = i1 + skyval + ip2 = i2 + skyval + is1 = i1 / pnorm + is2 = i2 / pnorm + default: + ip1 = i1 / pnorm + ip2 = i2 / pnorm + is1 = i1 + skyval + is2 = i2 + skyval + } + + # Draw the axes and axis labels. + call gsetr (gd, G_ASPECT, 0.0) + switch (xwcs) { + case WCS_XPIX: + call strcpy ("Radius (b=pixels, t=scale units)", Memc[pxlabel], + SZ_FNAME) + call strcpy ("", Memc[sxlabel], SZ_FNAME) + case WCS_XSCALE: + call strcpy ("Radius (b=scale units, t=pixels)", Memc[pxlabel], + SZ_FNAME) + call strcpy ("", Memc[sxlabel], SZ_FNAME) + default: + call strcpy ("Radius (b=pixels, t=scale units)", Memc[pxlabel], + SZ_FNAME) + call strcpy ("", Memc[sxlabel], SZ_FNAME) + } + switch (ywcs) { + case WCS_YCOUNT: + call strcpy ("Counts", Memc[pylabel], SZ_FNAME) + call strcpy ("Norm Counts", Memc[sylabel], SZ_FNAME) + case WCS_YNORM: + call strcpy ("Norm Counts", Memc[pylabel], SZ_FNAME) + call strcpy ("Counts", Memc[sylabel], SZ_FNAME) + default: + call strcpy ("Counts", Memc[pylabel], SZ_FNAME) + call strcpy ("Norm Counts", Memc[sylabel], SZ_FNAME) + } + + call gseti (gd, G_XDRAWAXES, 1) + call gseti (gd, G_YDRAWAXES, 1) + call gswind (gd, rp1, rp2, ip1, ip2) + call glabax (gd, title, Memc[pxlabel], Memc[pylabel]) + call gseti (gd, G_XDRAWAXES, 2) + call gseti (gd, G_YDRAWAXES, 2) + call gswind (gd, rs1, rs2, is1, is2) + call glabax (gd, title, Memc[sxlabel], Memc[sylabel]) + + # Plot the data points. + call gswind (gd, r1, r2, i1, i2) + call gpmark (gd, radius, intensity, npts, GM_PLUS, 1.0, 1.0) + + # Plot the smoothed radial profile skipping any points with no data. + call gswind (gd, r1, r2, -fraction, 1.0 + fraction) + for (i = 1; i <= nbins && rcentroid[i] <= 0.0; i = i + 1) + ; + call gamove (gd, rcentroid[i], pmean[i]) + do j = i, nbins { + if (pmean[j] == 0.0) + next + call gadraw (gd, rcentroid[j], pmean[j]) + } + + # Plot the integral. + call gswind (gd, r1, r2, -fraction, 1.0 + fraction) + call gamove (gd, 0.0, 0.0) + dr = (rmax - rmin) / real (nbins - 1) + do j = 2, nbins { + call gadraw (gd, real (j - 1) * dr, integral[j]) + } + + # Plot the sky annuli and the aperture radius. + call gamove (gd, iannulus, -fraction) + call gadraw (gd, iannulus, 1.0 + fraction) + call gamove (gd, oannulus, -fraction) + call gadraw (gd, oannulus, 1.0 + fraction) + call gamove (gd, apradius, -fraction) + call gadraw (gd, apradius, 1.0 + fraction) + + # Plot the zero level. + call gswind (gd, rp1, rp2, ip1, ip2) + switch (ywcs) { + case WCS_YCOUNT: + call gamove (gd, rp1, skyval) + call gadraw (gd, rp2, skyval) + case WCS_YNORM: + call gamove (gd, rp1, 0.0) + call gadraw (gd, rp2, 0.0) + default: + call gamove (gd, rp1, 0.0) + call gadraw (gd, rp2, 0.0) + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/daoedit/dperprofile.x b/noao/digiphot/daophot/daoedit/dperprofile.x new file mode 100644 index 00000000..37cccf97 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/dperprofile.x @@ -0,0 +1,590 @@ +include <mach.h> +include <imhdr.h> +include "daoedit.h" + + +# DP_ERPROFILE -- Compute and optionally plot the radial profile of the +# selected star. + +procedure dp_erprofile (gd, id, banner, xwcs, ywcs, im, xinit, yinit) + +pointer gd # pointer to the graphics descriptor +pointer id # pointer to the display descriptor +int banner # print banner for photometry results +int xwcs # the x wcs type +int ywcs # the y wcs type +pointer im # pointer to the input image +real xinit # initial x position +real yinit # initial y position + +int cboxsize, rboxsize, npts, naperts, nbins, nr +pointer sp, radius, intensity, apstr, aperts, rcentroid, pmean, integral, title +real scale, pradius, iannulus, oannulus, mean, median, sigma, pnorm, inorm +real fwhmpsf, xcenter, ycenter, apradius, counts +real clgetr(), dp_aprcounts() +int dp_rprofile(), strlen(), dp_gaperts(), btoi() +bool clgetb() + +begin + # Return if the initial x and y centers are undefined. + if (IS_INDEFR(xinit) || IS_INDEFR(yinit)) + return + + # Get the scale of the image from the datapars pset. + scale = clgetr ("datapars.scale") + if (scale <= 0.0) + scale = 1.0 + else + scale = 1.0 / scale + + # Estimate the center of the star using the datapars parameter + # scale and the centerpars parameter cbox and an initial center. + # Force the centering box to be an odd number of pixels. + cboxsize = scale * clgetr ("centerpars.cbox") + 1.0 + if (mod (cboxsize, 2) == 0) + cboxsize = cboxsize + 1 + call dp_rcenter (im, xinit, yinit, cboxsize, xcenter, ycenter) + + # Allocate working memory for the radial profile data. + pradius = scale * (clgetr ("fitskypars.annulus") + + clgetr ("fitskypars.dannulus") + 1.0) + rboxsize = 2 * nint (pradius + 1.0) + 1 + call smark (sp) + call salloc (radius, rboxsize * rboxsize, TY_REAL) + call salloc (intensity, rboxsize * rboxsize, TY_REAL) + + # Compute the radial profile. + npts = dp_rprofile (im, xcenter, ycenter, rboxsize, pradius, + Memr[radius], Memr[intensity]) + if (npts <= 0) { + call printf ("No data for computing radial profile\n") + call sfree (sp) + return + } + + # Estimate the sky level and subtract it from the profile. + iannulus = scale * clgetr ("fitskypars.annulus") + oannulus = iannulus + scale * clgetr ("fitskypars.dannulus") + call dp_rskyval (Memr[radius], Memr[intensity], npts, + iannulus, oannulus, mean, median, sigma) + call asubkr (Memr[intensity], median, Memr[intensity], npts) + + # Get the last photometry aperture and sum up the pixels inside it. + call salloc (apstr, SZ_LINE, TY_CHAR) + call salloc (aperts, MAX_NAPERTS, TY_REAL) + call clgstr ("photpars.apertures", Memc[apstr], SZ_LINE) + naperts = dp_gaperts (Memc[apstr], Memr[aperts], MAX_NAPERTS) + apradius = scale * Memr[aperts+naperts-1] + counts = dp_aprcounts (Memr[radius], Memr[intensity], npts, apradius) + + # Compute the normalized mean intensity as a function of radius + # and the normalized integral of the intensity as a function of + # radius. + nbins = int (pradius / 0.5 + 0.5) + call salloc (nr, nbins, TY_INT) + call salloc (rcentroid, nbins, TY_REAL) + call salloc (pmean, nbins, TY_REAL) + call salloc (integral, nbins, TY_REAL) + call dp_crprofile (Memr[radius], Memr[intensity], npts, Memi[nr], + Memr[rcentroid], Memr[pmean], Memr[integral], nbins, 0.0, + pradius, real ((cboxsize-1)/ 2), iannulus, fwhmpsf, pnorm, inorm) + + # Mark the objects on the display. + if (id != NULL) { + call dp_eomark (id, xcenter, ycenter, iannulus, oannulus, + apradius, btoi (clgetb ("centerpars.mkcenter")), + btoi (clgetb ("fitskypars.mksky")), + btoi (clgetb ("photpars.mkapert"))) + if (gd == id) + call gflush (gd) + else + call gframe (id) + } + # Convert the center coordinates if appropriate. + call dp_ltov (im, xcenter, ycenter, xcenter, ycenter, 1) + + # Draw the plot. + if (gd != NULL) { + call salloc (title, 3 * SZ_LINE, TY_CHAR) + call sprintf (Memc[title], 3 * SZ_LINE, + "%s: Radial profile at %0.2f %0.2f\n") + call pargstr (IM_HDRFILE(im)) + call pargr (xcenter) + call pargr (ycenter) + call sprintf (Memc[title+strlen(Memc[title])], 3 * SZ_LINE, + "Sky: mean %g median %g sigma %g\n") + call pargr (mean) + call pargr (median) + call pargr (sigma) + call sprintf (Memc[title+strlen(Memc[title])], 3 * SZ_LINE, + "Fwhmpsf: %0.2f Counts: %g Mag: %0.3f\n\n") + call pargr (fwhmpsf) + call pargr (counts) + if (counts <= 0.0) + call pargr (INDEFR) + else + call pargr (-2.5 * log10 (counts)) + call dp_erplot (gd, Memc[title], xwcs, ywcs, Memr[radius], + Memr[intensity], npts, Memr[rcentroid], Memr[pmean], + Memr[integral], nbins, 0.0, pradius, iannulus, oannulus, + apradius, median, sigma, scale, pnorm) + } + + # Print the results. + if (gd == NULL) { + if (banner == YES) + call dp_bprint (STDOUT) + call dp_aprint (STDOUT, xcenter, ycenter, median, sigma, fwhmpsf, + counts) + } + + # Free memory. + call sfree (sp) +end + + +# DP_RCENTER -- Compute the star center using a simple 1D centroiding +# algorithm on the x and y marginals, after thresholding at the mean. + +procedure dp_rcenter (im, xstart, ystart, boxsize, xcntr, ycntr) + +pointer im # pointer to the input image +real xstart, ystart # starting coordinates +int boxsize # width of the centering box +real xcntr, ycntr # centered coordinates + +int half_box, x1, x2, y1, y2, ncols, nrows, nx, ny, try +pointer bufptr, sp, x_vect, y_vect +real xinit, yinit +pointer imgs2r() + +begin + # Initialize. + half_box = (boxsize - 1) / 2 + xinit = xstart + yinit = ystart + ncols = IM_LEN (im,1) + nrows = IM_LEN (im,2) + + try = 0 + repeat { + + # Compute the parameters of the extraction region. + x1 = max (xinit - half_box, 1.0) + 0.5 + x2 = min (xinit + half_box, real (ncols)) + 0.5 + y1 = max (yinit - half_box, 1.0) + 0.5 + y2 = min (yinit + half_box, real (nrows)) + 0.5 + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + + # Get the data. + bufptr = imgs2r (im, x1, x2, y1, y2) + + # Allocate space for the marginals. + call smark (sp) + call salloc (x_vect, nx, TY_REAL) + call salloc (y_vect, ny, TY_REAL) + + # Compute the marginals. + call aclrr (Memr[x_vect], nx) + call aclrr (Memr[y_vect], ny) + call dp_rowsum (Memr[bufptr], Memr[x_vect], nx, ny) + call dp_colsum (Memr[bufptr], Memr[y_vect], nx, ny) + + # Compute the centers. + call dp_vcentroid (Memr[x_vect], nx, xcntr) + call dp_vcentroid (Memr[y_vect], ny, ycntr) + + # Add in offsets to image coordinate system. + xcntr = xcntr + x1 + ycntr = ycntr + y1 + + call sfree (sp) + + # If the shifts are greater than 1 pixel in either direction + # do 1 more iteration. + try = try + 1 + if (try == 1) { + if ((abs (xcntr - xinit) > 1.0) || (abs (ycntr - yinit) > + 1.0)) { + xinit = xcntr + yinit = ycntr + } + } else + break + } +end + + +# DP_RPROFILE -- Get the data and compute the radius and intensity vectors. + +int procedure dp_rprofile (im, xcntr, ycntr, rboxsize, pradius, radius, + intensity) + +pointer im # pointer to the input image +real xcntr, ycntr # the center of the extraction box +int rboxsize # the width of the extraction box +real pradius # the plotting radius +real radius[ARB] # the output radius vector +real intensity[ARB] # the output intensity vector + +int half_box, ncols, nrows, x1, x2, y1, y2, nx, ny, npts +pointer bufptr +real xinit, yinit +int dp_rivectors() +pointer imgs2r() + +begin + # Initialize. + half_box = (rboxsize - 1) / 2 + xinit = xcntr + yinit = ycntr + ncols = IM_LEN(im,1) + nrows = IM_LEN(im,2) + + # Get the data. + x1 = max (xinit - half_box, 1.0) + 0.5 + x2 = min (xinit + half_box, real (ncols)) + 0.5 + y1 = max (yinit - half_box, 1.0) + 0.5 + y2 = min (yinit + half_box, real (nrows)) + 0.5 + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + bufptr = imgs2r (im, x1, x2, y1, y2) + + # Compute the radius and intensity vectors. + npts = dp_rivectors (Memr[bufptr], nx, ny, x1, y1, xcntr, ycntr, + pradius, radius, intensity) + + return (npts) +end + + +# DP_RSKYVAL -- Compute the mean, median and sigma of the pixels in the +# sky annulus. + +procedure dp_rskyval (radius, intensity, npts, iannulus, oannulus, mean, + median, sigma) + +real radius[ARB] # the output radius vector +real intensity[ARB] # the output intensity vector +int npts # number of points in the profile +real iannulus # inner radius of sky annulus +real oannulus # outer radius of sky annulus +real mean # mean sky value +real median # median sky value +real sigma # standard deviation of sky values + +int i, nsky, il, ih, ngood +pointer skypix +real slocut, shicut + +begin + call malloc (skypix, npts, TY_REAL) + + nsky = 0 + do i = 1, npts { + if (radius[i] < iannulus || radius[i] > oannulus) + next + nsky = nsky + 1 + Memr[skypix+nsky-1] = intensity[i] + } + call asrtr (Memr[skypix], Memr[skypix], nsky) + + if (nsky == 0) { + mean = 0.0 + median = 0.0 + sigma = 0.0 + } else { + call aavgr (Memr[skypix], nsky, mean, sigma) + if (mod (nsky, 2) == 0) + median = 0.5 * (Memr[skypix+(nsky+1)/2-1] + + Memr[skypix+(nsky+1)/2]) + else + median = Memr[skypix+(nsky+1)/2-1] + } + + # Detect pixels to be rejected. + slocut = median - min (median - Memr[skypix], Memr[skypix+nsky-1] - + median, 3.0 * sigma) + shicut = median + min (median - Memr[skypix], Memr[skypix+nsky-1] - + median, 3.0 * sigma) + for (il = 1; il <= nsky; il = il + 1) { + if (Memr[skypix+il-1] >= slocut) + break + } + for (ih = nsky; ih >= 1; ih = ih - 1) { + if (Memr[skypix+ih-1] <= shicut) + break + } + ngood = ih - il + 1 + if (ngood < nsky) { + if (ngood == 0) { + mean = 0.0 + median = 0.0 + sigma = 0.0 + } else { + call aavgr (Memr[skypix+il-1], ngood, mean, sigma) + if (mod (ngood, 2) == 0) + median = 0.5 * (Memr[skypix+il-1+(ngood+1)/2-1] + + Memr[skypix+il-1+(ngood+1)/2]) + else + median = Memr[skypix+il-1+(ngood+1)/2-1] + } + } + + + call mfree (skypix, TY_REAL) +end + + +# DP_APRCOUNTS -- Sumup the counts inside the aperture. + +real procedure dp_aprcounts (radius, intensity, npts, apradius) + +real radius[ARB] # the output radius vector +real intensity[ARB] # the output intensity vector +int npts # number of points in the profile +real apradius # the aperture radius + +int i +real counts + +begin + counts = 0.0 + do i = 1, npts { + if (radius[i] > apradius) + next + counts = counts + intensity[i] + } + return (counts) +end + + +# DP_CRPROFILE -- Compute the smoothed radial profile and its integral at half +# pixel intervals. + +procedure dp_crprofile (radius, intensity, npts, nr, rcentroid, pmean, + integral, nbins, rmin, rmax, prad, irad, fwhmpsf, pnorm, inorm) + +real radius[ARB] # the output radius vector +real intensity[ARB] # the output intensity vector +int npts # number of points in the profile +int nr[ARB] # the number of points in each interval +real rcentroid[ARB] # the centroid of the radius values +real pmean[ARB] # the mean of the intensity values +real integral[ARB] # the integral of the curve +int nbins # the number of radius bins +real rmin, rmax # the radius min and max values +real prad # the normalization radius for the profile +real irad # the normalization radius for the integral +real fwhmpsf # computed full width halfmax psf +real pnorm # the maximum profile value +real inorm # the maximum count value + +int i, bin +real dr, r + +begin + # Initialize the arrays. + call aclri (nr, nbins) + call aclrr (rcentroid, nbins) + call aclrr (pmean, nbins) + call aclrr (integral, nbins) + + # Accumulate the data. + dr = real (nbins - 1) / real (rmax - rmin) + do i = 1, npts { + r = radius[i] + if (r < rmin || r > rmax) + next + bin = int ((r - rmin) * dr) + 1 + nr[bin] = nr[bin] + 1 + pmean[bin] = pmean[bin] + intensity[i] + rcentroid[bin] = rcentroid[bin] + r + } + + # Compute the integral of the radial profile and normalize it + # to 1.0 at the radius irad. + do i = 2, nbins + integral[i] = integral[i-1] + pmean[i-1] + bin = int ((irad - rmin) * dr) + 1 + inorm = integral[min (bin, nbins)] + call adivkr (integral, inorm, integral, nbins) + + # Compute the smoothed radial profile and normalize to the + # maximum data point inside the radius prad. + pnorm = -MAX_REAL + do i = 1, npts { + if (radius[i] > prad) + next + if (intensity[i] > pnorm) + pnorm = intensity[i] + } + do i = 1, nbins { + if (nr[i] <= 0) { + rcentroid[i] = (i - 1.0 + 0.5) / dr + if (rcentroid[i] < rmin || rcentroid[i] > rmax) + pmean[i] = 0.0 + else + pmean[i] = pnorm + } else { + rcentroid[i] = rcentroid[i] / nr[i] + pmean[i] = pmean[i] / nr[i] + } + } + call adivkr (pmean, pnorm, pmean, nbins) + + # Estimate the full width half max of the psf in pixels. + do i = 1, nbins { + if (pmean[i] == 0.0) + next + if (pmean[i] < 0.5) + break + } + if (i == 1) + fwhmpsf = 2.0 * rcentroid[1] + else if (i == nbins && pmean[nbins] >= 0.5) + fwhmpsf = 2.0 * rcentroid[nbins] + else if (pmean[i-1] == pmean[i]) + fwhmpsf = rcentroid[i-1] + rcentroid[i] + else + fwhmpsf = 2.0 * ((rcentroid[i] * (0.5 - pmean[i-1]) + + rcentroid[i-1] * (pmean[i] - 0.5)) / (pmean[i] - pmean[i-1])) +end + + +# DP_ROWSUM -- Sum all the rows in a raster. + +procedure dp_rowsum (raster, row, nx, ny) + +real raster[nx,ny] # the input subraster +real row[ARB] # the output summed row +int nx, ny # the dimensions of the input subraster + +int i, j + +begin + do i = 1, ny + do j = 1, nx + row[j] = row[j] + raster[j,i] +end + + +# DP_COLSUM -- Sum all the columns in a raster. + +procedure dp_colsum (raster, col, nx, ny) + +real raster[nx,ny] # the input subraster +real col[ARB] # the output summed column +int nx, ny # the dimensions of the input subraster + +int i, j + +begin + do i = 1, ny + do j = 1, nx + col[j] = col[j] + raster[i,j] +end + + +# DP_VCENTROID -- Compute the centroid of a vector. + +procedure dp_vcentroid (vector, nv, vc) + +real vector[ARB] # the input vector +int nv # length of the input array +real vc # the output centroid + +int i +real sum1, sum2, sigma, cont + +begin + sum1 = 0.0 + sum2 = 0.0 + + call aavgr (vector, nv, cont, sigma) + do i = 1, nv + if (vector[i] > cont) { + sum1 = sum1 + (i - 1) * (vector[i] - cont) + sum2 = sum2 + (vector[i] - cont) + } + + vc = sum1 / sum2 +end + + +# DP_RIVECTORS -- Compute the radius and intensity vectors. + +int procedure dp_rivectors (a, nx, ny, x1, y1, xcntr, ycntr, pradius, + radius, intensity) + +real a[nx,ny] # the input data array +int nx, ny # dimensions of the input array +int x1, y1 # lower left corner of input array +real xcntr, ycntr # coordinates of center pixel +real pradius # the plotting radius +real radius[ARB] # the output radius vector +real intensity[ARB] # the output intensity vector + +int i, j, npts +real pr2, r2, dy2 + +begin + pr2 = pradius * pradius + + npts = 0 + do i = 1, ny { + dy2 = (ycntr - y1 + 1 - i) ** 2 + do j = 1, nx { + r2 = (xcntr - x1 + 1 - j) ** 2 + dy2 + if (r2 > pr2) + next + npts = npts + 1 + radius[npts] = sqrt (r2) + intensity[npts] = a[j,i] + } + } + + return (npts) +end + + +# DP_BPRINT -- Print the photometry banner. + +procedure dp_bprint (fd) + +int fd # output file descriptor + +string banner "# XCENTER YCENTER SKY SKYSIGMA FWHM COUNTS MAG" + +begin + call fprintf (fd, "\n%s\n") + call pargstr (banner) +end + + +# DP_APRINT -- Print the photometry results. + +procedure dp_aprint (fd, xcenter, ycenter, skyval, sigma, fwhmpsf, inorm) + +int fd # output file descriptor +real xcenter # x coordinate of profile +real ycenter # y coordinate of profile +real skyval # the sky value +real sigma # the standard deviation of the sky pixels +real fwhmpsf # the estimated fwhmpsf +real inorm # the total counts + +begin + call fprintf (fd, " %7.2f %7.2f %8.1f %8.2f %6.2f %8.1f %7.3f\n") + call pargr (xcenter) + call pargr (ycenter) + call pargr (skyval) + call pargr (sigma) + call pargr (fwhmpsf) + call pargr (inorm) + if (inorm <= 0.0) + call pargr (INDEFR) + else + call pargr (-2.5 * log10 (inorm)) +end diff --git a/noao/digiphot/daophot/daoedit/idaoedit.key b/noao/digiphot/daophot/daoedit/idaoedit.key new file mode 100644 index 00000000..5967d1e8 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/idaoedit.key @@ -0,0 +1,22 @@ + Interactive Daoedit Parameter Setup Menu + +? Print help +spbar Mark/verify critical parameters (f, t, a, d, r, w, b) +q Quit + +f Mark/verify the fwhm of the psf on the radial profile plot +s Mark/verify the sky sigma on the radial profile plot +l Mark/verify the minimum good data value on the radial profile plot +u Mark/verify the maximum good data value on the radial profile plot + +c Mark/verify the centering box on the radial profile plot +n Mark/verify the cleaning radius on the radial profile plot +p Mark/verify the clipping radius on the radial profile plot + +a Mark/verify the inner sky annulus radius on the radial profile plot +d Mark/verify the width of the sky annulus on the radial profile plot +g Mark/verify the sky region growing radius on the radial profile plot + +r Mark/verify the photometry aperture(s) on the radial profile plot +w Mark/verify the psf function radius on the radial profile plot +b Mark/verify the psf fitting radius on the radial profile plot diff --git a/noao/digiphot/daophot/daoedit/mkpkg b/noao/digiphot/daophot/daoedit/mkpkg new file mode 100644 index 00000000..0fc962d1 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/mkpkg @@ -0,0 +1,16 @@ +# DAOEDIT task + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + t_daoedit.x <imhdr.h> <gset.h> ../lib/daophotdef.h daoedit.h + dpecolon.x <error.h> daoedit.h + dpemark.x <ctype.h> daoedit.h + dpeomark.x <gset.h> + dpeconfirm.x daoedit.h + dperprofile.x <mach.h> <imhdr.h> daoedit.h + dperplot.x <gset.h> daoedit.h + ; diff --git a/noao/digiphot/daophot/daoedit/t_daoedit.x b/noao/digiphot/daophot/daoedit/t_daoedit.x new file mode 100644 index 00000000..8c92e969 --- /dev/null +++ b/noao/digiphot/daophot/daoedit/t_daoedit.x @@ -0,0 +1,317 @@ +include <imhdr.h> +include <gset.h> +include "../lib/daophotdef.h" +include "daoedit.h" + +define HELPFILE "daophot$daoedit/daoedit.key" +define IHELPFILE "daophot$daoedit/daoedit.key" + +# T_DAOEDIT -- Edit the DAOPHOT parameters interactively using the image +# display and a radial profile plot. + +procedure t_daoedit () + +pointer image # the name of the input image +int cache # cache the image pixels +pointer graphics # the graphics device +pointer display # the image display + +real wx, wy, xlast, ylast +pointer sp, cmd, im, gd, id +int wcs, key, redraw, gcurtype, curtype, xwcs, ywcs, lastkey +int req_size, memstat, old_size, buf_size + +pointer immap(), gopen() +int dp_gcur(), btoi(), sizeof(), dp_memstat() +bool clgetb(), streq() +errchk gopen() + +data gcurtype /'g'/ + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (display, SZ_FNAME, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call clgstr ("image", Memc[image], SZ_FNAME) + im = immap (Memc[image], READ_ONLY, 0) + + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + iferr (gd = gopen (Memc[graphics], AW_DEFER+NEW_FILE, STDGRAPH)) + gd = NULL + call clgstr ("display", Memc[display], SZ_FNAME) + if (Memc[display] == EOS) + id = NULL + else if (streq (Memc[graphics], Memc[display])) + id = gd + else { + iferr { + id = gopen (Memc[display], APPEND, STDIMAGE) + } then { + call eprintf ( + "Warning: Graphics overlay not available for display device.\n") + id = NULL + } + } + cache = btoi (clgetb ("cache")) + + # Set up the display coordinate system. + if ((id != NULL) && (id != gd)) + call dp_gswv (id, Memc[image], im, 4) + + # Cache the input image pixels. + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = dp_memstat (cache, req_size, old_size) + if (memstat == YES) + call dp_pcache (im, INDEFI, buf_size) + + xlast = INDEFR + ylast = INDEFR + lastkey = INDEFI + xwcs = WCS_XPIX + ywcs = WCS_YCOUNT + curtype = 'i' + redraw = NO + + while (dp_gcur (curtype, wx, wy, wcs, key, Memc[cmd], SZ_LINE) != EOF) { + + # Convert the cursor coordinates if necessary. + if (curtype == 'i') + call dp_vtol (im, wx, wy, wx, wy, 1) + + switch (key) { + + # Print help page. + case '?': + if (curtype == 'i') { + call pagefile (HELPFILE, "[space=cmhelp,q=quit,?=help]") + } else if (curtype == 'g') { + #call greactivate (gd, 0) + call gpagefile (gd, HELPFILE, "") + } + + + # Quit the program. + case 'q': + break + + # Toggle between the image display and graphics cursor. + case 'g': + if (curtype == 'i') { + call greactivate (gd, 0) + curtype = 'g' + } else { + call gdeactivate (gd, 0) + curtype = 'i' + } + + # Toggle between the pixel and scale units plot scale in x. + case 'x': + if (xwcs == WCS_XPIX) + xwcs = WCS_XSCALE + else if (xwcs == WCS_XSCALE) + xwcs = WCS_XPIX + redraw = YES + + # Toggle between the counts and norm plot scale in y. + case 'y': + if (ywcs == WCS_YCOUNT) + ywcs = WCS_YNORM + else if (ywcs == WCS_YNORM) + ywcs = WCS_YCOUNT + redraw = YES + + case 'a': + if (curtype == 'i') { + if (lastkey == 'a') + call dp_erprofile (NULL, id, NO, xwcs, ywcs, im, wx, wy) + else + call dp_erprofile (NULL, id, YES, xwcs, ywcs, im, + wx, wy) + } else { + call dp_erprofile (NULL, id, NO, xwcs, ywcs, im, + xlast, ylast) + } + + case 'r': + if (curtype == 'i') { + xlast = wx + ylast = wy + } + redraw = YES + + case 'i': + if (curtype == 'i') { + xlast = wx + ylast = wy + } + xwcs = WCS_XPIX + ywcs = WCS_YCOUNT + redraw = YES + + case ':': + call dp_ecolon (Memc[cmd], gd, redraw) + + default: + call printf ("Unknown or ambiguous keystroke command\7\n") + } + + # Draw the plot. + if (redraw == YES) { + call dp_erprofile (gd, id, NO, xwcs, ywcs, im, xlast, ylast) + if (key == 'i') + call dp_isetup (gcurtype, gd) + redraw = NO + if (curtype == 'i') + call gdeactivate (gd, 0) + } + + lastkey = key + } + + if (id == gd && gd != NULL) { + call gclose (gd) + } else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + + call imunmap (im) + + # Uncache memory. + call fixmem (old_size) + + call sfree (sp) +end + + +# DP_GCUR --Get a SETPARS cursor value. + +int procedure dp_gcur (curtype, x, y, wcs, key, strval, maxch) + +int curtype # cursor type +real x, y # cursor position +int wcs # cursor wcs +int key # keystroke value of cursor event +char strval[ARB] # string value, if any +int maxch # max chars out + +int nitems +int clgcur() + +begin + # Initialize. + strval[1] = EOS + + # Get a cursor values from the desired cursor parameter. + switch (curtype) { + case 'i': + nitems = clgcur ("icommands", x, y, wcs, key, strval, maxch) + case 'g': + nitems = clgcur ("gcommands", x, y, wcs, key, strval, maxch) + } + + return (nitems) +end + + +# DP_ISETUP -- Setup the daophot parameters interactively using a +# radial profile plot. + +procedure dp_isetup (curtype, gd) + +int curtype # the cursor type graphics or display +pointer gd # pointer to the graphics stream + +int wcs, key +pointer sp, cmd +real wx, wy +int dp_gcur() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call printf ( + "Interactive setup menu (?=help, spbar=default, q=quit):\n") + while (dp_gcur (curtype, wx, wy, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + switch (key) { + case '?': + call gpagefile (gd, IHELPFILE, "") + case 'q': + break + case 'f': + call dp_mfwhmpsf(gd) + call dp_cfwhmpsf() + case 's': + call dp_msigma(gd) + call dp_csigma() + case 'u': + call dp_mdmax(gd) + call dp_cdmax() + case 'l': + call dp_mdmin(gd) + call dp_cdmin() + case 'c': + call dp_mcbox(gd) + call dp_ccbox() + case 'n': + call dp_mrclean(gd) + call dp_crclean() + case 'p': + call dp_mrclip(gd) + call dp_crclip() + case 'a': + call dp_mannulus(gd) + call dp_cannulus() + case 'd': + call dp_mdannulus(gd) + call dp_cdannulus() + case 'g': + call dp_mrgrow(gd) + call dp_crgrow() + case 'r': + call dp_maper(gd) + call dp_caper() + case 'w': + call dp_mpsfrad(gd) + call dp_cpsfrad() + case 'b': + call dp_mfitrad(gd) + call dp_cfitrad() + case ' ': + call dp_mfwhmpsf(gd) + call dp_mcbox(gd) + call dp_mannulus(gd) + call dp_mdannulus(gd) + call dp_maper(gd) + call dp_mpsfrad(gd) + call dp_mfitrad(gd) + call gdeactivate(gd, 0) + call dp_cbanner() + call dp_cfwhmpsf() + call dp_ccbox() + call dp_cannulus() + call dp_cdannulus() + call dp_caper() + call dp_cpsfrad() + call dp_cfitrad() + call greactivate(gd, 0) + default: + call printf ( + "Interactive setup menu (?=help, spbar=default, q=quit):\n") + } + call printf ( + "Interactive setup menu (?=help, spbar=default, q=quit):\n") + } + + if (curtype == 'i') + call gdeactivate (gd, 0) + call sfree (sp) +end |