aboutsummaryrefslogtreecommitdiff
path: root/noao/obsutil/src/starfocus
diff options
context:
space:
mode:
authorJoe Hunkeler <jhunkeler@gmail.com>2015-08-11 16:51:37 -0400
committerJoe Hunkeler <jhunkeler@gmail.com>2015-08-11 16:51:37 -0400
commit40e5a5811c6ffce9b0974e93cdd927cbcf60c157 (patch)
tree4464880c571602d54f6ae114729bf62a89518057 /noao/obsutil/src/starfocus
downloadiraf-osx-40e5a5811c6ffce9b0974e93cdd927cbcf60c157.tar.gz
Repatch (from linux) of OSX IRAF
Diffstat (limited to 'noao/obsutil/src/starfocus')
-rw-r--r--noao/obsutil/src/starfocus/Revisions162
-rw-r--r--noao/obsutil/src/starfocus/mkpkg22
-rw-r--r--noao/obsutil/src/starfocus/psfhelp.key60
-rw-r--r--noao/obsutil/src/starfocus/psfmeasure.par24
-rw-r--r--noao/obsutil/src/starfocus/starfocus.h140
-rw-r--r--noao/obsutil/src/starfocus/starfocus.key15
-rw-r--r--noao/obsutil/src/starfocus/starfocus.par32
-rw-r--r--noao/obsutil/src/starfocus/stfgraph.x2682
-rw-r--r--noao/obsutil/src/starfocus/stfhelp.key63
-rw-r--r--noao/obsutil/src/starfocus/stfmeasure.x134
-rw-r--r--noao/obsutil/src/starfocus/stfprofile.x1189
-rw-r--r--noao/obsutil/src/starfocus/t_starfocus.x1240
-rw-r--r--noao/obsutil/src/starfocus/x_starfocus.x2
13 files changed, 5765 insertions, 0 deletions
diff --git a/noao/obsutil/src/starfocus/Revisions b/noao/obsutil/src/starfocus/Revisions
new file mode 100644
index 00000000..d3f93e90
--- /dev/null
+++ b/noao/obsutil/src/starfocus/Revisions
@@ -0,0 +1,162 @@
+.help revisions Nov01 obsutil
+.nf
+t_starfocus.x
+starfocus.par
+psfmeasure.par
+starfocus.hlp
+psfmeasure.hlp
+ Added a "wcs" parameter to allow specifying the coordinate system type
+ for cursor input. (10/7/98, Valdes)
+
+stfprofile.x
+ The logic in STF_FIT for determining the points to fit and the point
+ to use for the initial width estimate was faulty allowing some bad
+ cases to get through. (7/31/98, Valdes)
+
+t_starfocus.x
+ There was an error in the changes to allow focus sequences to go
+ up or down such that simple PSF measurement failed.
+ (6/11/97, Valdes)
+
+stfmeasure.x
+stfprofile.x
+ Changes were need to better deal with error conditions.
+ (6/11/97, Valdes)
+
+stfmeasure.x
+ Changes to deal with error conditions.
+ (6/11/97, Valdes)
+
+t_starfocus.x
+ Fixed a typo error in an earlier change which left out the pointer
+ variable in a macro. (5/8/97, Valdes)
+
+t_starfocus.x
+ 1. Added a minimum radius of 3 to the input and to interative setting.
+ 2. The estimate for the next focus position is not based on the previous
+ center plus the step.
+ 3. Added commented out code to allow setting the selected star to always
+ be the top or bottom, etc.
+ (2/7/97, Valdes)
+
+stfprofile.x
+ 1. Added a minimum radius of 3 to stf_find.
+ 2. Added a error return if sum2 <= 0 in stf_find
+ 3. Added errchks in stf_widths for called procedures.
+ 4. Added error returns in stf_fit for bad profile or parameters.
+ (2/6/97, Valdes)
+
+t_starfocus.x
+ The parabola fitting routine would fail if the independent variable
+ (the focus values) became too large because a) use of real and b)
+ unscaled variables. The routine was revised for both these problems.
+ (10/24/96, Valdes)
+
+stfprofile.x
+stfmeasure.x
+ Fixed bug in evaluation of enclosed flux profile in which the scaled
+ radius was used for the gaussian subtraction stage instead of pixels.
+ This does not currently affect IMEXAM because the scale is fixed
+ at 1. (8/29/96, Valdes)
+
+stfmeasure.x
+stfprofile.x
+ Fixed incorrect datatype declaration "real np" -> "int np" in various
+ related places. (4/9/96, Valdes)
+
+stfmeasure.x -
+ 1. Restricted the peak pixel to be within the specified radius.
+ 2. Renamed file and procedures to avoid library conflict.
+ (3/14/96, Valdes)
+
+stfprofile.x
+ Minor bug fix that does not affect anything. (3/14/96, Valdes)
+
+kpnofocus.cl
+starfocus.par
+psfmeasure.par
+ Changed the defaults for radius=5, sbuffer=5, swidth=5 to try and make
+ IMEXAMINE and the PSF tasks measure the same thing by default.
+ (3/13/96, Valdes)
+
+stfprofile.x
+ Changed the centering routine to only use the data in the specified
+ radius rather than the extended region including the sky. (11/7/95, Valdes)
+
+stfmeasure.x
+mkpkg
+ Added a routine which can be called separately to compute the
+ enclosed flux width or radius and the direct FWHM. (11/6/95, Valdes)
+
+stfprofile.x
+starfocus.h
+ Added a computation of the direct FWHM. (11/6/95, Valdes)
+
+t_starfocus.x
+ Make a change to avoid reloading the image display when the specified
+ image name and the image name in the display are same apart from
+ the image extension. (6/26/95, Valdes)
+
+stfprofile.h
+t_starfocus.x
+stfprofile.x
+starfocus.par
+psfmeasure.par
+kpnofocus.cl
+starfocus.par
+psfmeasure.par
+kpnofocus.cl
+ Added saturation and ignore_sat parameters to allow flagging measurements
+ with saturated pixels and either ignoring them or including them but
+ showing them in the output log. (10/28/94, Valdes)
+
+t_starfocus.x
+ Fixed bug in interpreting the focus parameter. (9/15/94, Valdes)
+
+starfocus.h
+t_starfocus.x
+stfprofile.x
+stfgraph.x
+psfmeasure.par
+starfocus.par
+kpnofocus.cl
+psfmeasure.hlp
+starfocus.hlp
+kpnofocus.hlp
+ Added an "iteration" parameter that, if greater than 1, uses the
+ previous FWHM estimate to adjust the "radius" parameter. (8/5/94, Valdes)
+
+psfmeasure.hlp
+starfocus.hlp
+ Clarified the description of the imagecur parameter. (2/1/94, Valdes)
+
+t_starfocus.x
+starfocus.h
+starfocus.par
+psfmeasure.par
+kpnofocus.cl +
+starfocus.hlp
+psfmeasure.hlp
+kpnofocus.hlp
+ 1. A new parameter, fstep, was added to STARFOCUS to allow specifying
+ the focus sequence as a starting value, by the focus parameters,
+ and a step.
+ 2. STARFOCUS was modified to allow specification of header keywords
+ for the focus (this was true previously), the focus step,
+ the number of exposures, and the multiple exposure shift.
+ This allows multiple exposure images to be completely header
+ drivien if the appropriate keywords are present.
+ 3. A new task KPNOFOCUS was added. This is a script calling
+ STARFOCUS with parameters set specifically for Kitt Peak
+ headers containing the focus parameters in the header.
+ Many of the parameters are fixed in the script and the task
+ parameters are then simpler.
+ 4. STARFOCUS/PSFMEASURE were modified to search all display frames
+ for the requested image rather than just frame 1.
+ 5. A new parameters, frame, was added to STARFOCUS/PSFMEASURE
+ to specify the display frame to load if necessary. Previously
+ it was always frame 1.
+ 7. In STARFOCUS the default number of steps was changed to 7 and
+ the default step size to 30.
+ (11/13/93, Valdes)
+.endhelp
diff --git a/noao/obsutil/src/starfocus/mkpkg b/noao/obsutil/src/starfocus/mkpkg
new file mode 100644
index 00000000..3feffc76
--- /dev/null
+++ b/noao/obsutil/src/starfocus/mkpkg
@@ -0,0 +1,22 @@
+# Make the STARFOCUS tasks.
+
+$checkout libpkg.a ../
+$update libpkg.a
+$checkin libpkg.a ../
+$exit
+
+standalone:
+ $update libpkg.a
+ $omake x_starfocus.x
+ $link x_starfocus.o libpkg.a -lxtools -lnlfit -liminterp -lgsurfit\
+ -lds -o xx_starfocus.e
+ ;
+
+libpkg.a:
+ stfgraph.x starfocus.h <error.h> <gset.h> <mach.h>
+ stfmeasure.x starfocus.h <error.h> <imhdr.h> <imset.h>\
+ <math/nlfit.h>
+ stfprofile.x starfocus.h <imhdr.h> <mach.h> <math.h>\
+ <math/iminterp.h> <math/nlfit.h>
+ t_starfocus.x starfocus.h <error.h> <imhdr.h> <imset.h> <mach.h>
+ ;
diff --git a/noao/obsutil/src/starfocus/psfhelp.key b/noao/obsutil/src/starfocus/psfhelp.key
new file mode 100644
index 00000000..9e617a62
--- /dev/null
+++ b/noao/obsutil/src/starfocus/psfhelp.key
@@ -0,0 +1,60 @@
+ PSFMEASURE COMMAND OPTIONS
+
+ SUMMARY
+
+? Help m Magnitude r Redraw z Zoom
+a Spatial n Normalize s Mag symbols <space> Next
+d Delete o Offset t Field radius
+e Enclosed flux p Radial profile u Undelete
+i Info q Quit x Delete
+
+:level <val> :radius <val> :show <file> :xcenter <val>
+:overplot <y|n> :scale <val> :size <type> :ycenter <val>
+
+
+ CURSOR COMMANDS
+
+All plots may not be available depending on the number of stars.
+
+? Page this help information
+a Spatial plot
+d Delete star nearest to cursor
+e Enclosed flux for all stars
+i Information about star nearest the cursor
+m Size and ellipticity vs relative magnitude
+n Normalize enclosed flux at x cursor position
+o Offset enclosed flux to x,y cursor position by adjusting background
+p Radial profiles for all stars
+ Profiles are determined from the derivatives of the enclosed flux
+q Quit
+r Redraw
+s Toggle magnitude symbols in spatial plot
+t Size and ellipticity vs radius from field center
+u Undelete all deleted points
+x Delete nearest point or star (selected by query)
+z Zoom to a single measurement showing enclosed flux and radial profile
+<space> Step through different stars in some plots
+
+
+ COLON COMMANDS
+
+A command without a value generally shows the current value of the
+parameter while with a value it sets the value of the parameter.
+
+:level <val> Level at which the size parameter is evaluated
+:overplot <y|n> Overplot the profiles from the narrowest profile?
+:radius <val> Change profile radius(*)
+:show <file> Page all information for the current set of objects
+:size <type> Size type (Radius|FWHM|GFWHM|MFWHM) (**)
+:scale <val> Pixel scale for size values
+:xcenter <val> X field center for radius from field center plots
+:ycenter <val> Y field center for radius from field center plots
+
+(*) The profile radius may not exceed the initial value set by the task
+ parameter.
+
+(**)
+Radius = radius enclosing the fraction of the flux specified by "level"
+FWHM = Full-width at half-maximum based on the radially smoothed profile
+GFWHM = Full-width at half-maximum of Gaussian function fit to enclosed flux
+MFWHM = Full-width at half-maximum of Moffat function fit to enclosed flux
diff --git a/noao/obsutil/src/starfocus/psfmeasure.par b/noao/obsutil/src/starfocus/psfmeasure.par
new file mode 100644
index 00000000..16dea3e9
--- /dev/null
+++ b/noao/obsutil/src/starfocus/psfmeasure.par
@@ -0,0 +1,24 @@
+# PSFMEASURE -- Measure PSF.
+
+images,s,a,,,,"List of images"
+coords,s,h,"markall","center|mark1|markall",,"Object coordinates"
+wcs,s,h,"logical","logical|physical|world",,"Coordinate system"
+display,b,h,yes,,,"Display images?"
+frame,i,h,1,1,,"Display frame to use
+"
+level,r,h,0.5,,,"Measurement level (fraction or percent)"
+size,s,h,"FWHM","Radius|FWHM|GFWHM|MFWHM",,"Size to display"
+beta,r,h,INDEF,2.1,,Moffat beta parameter
+scale,r,h,1.,,,"Pixel scale"
+radius,r,h,5.,,,"Measurement radius (pixels)"
+sbuffer,r,h,5.,,,"Sky buffer (pixels)"
+swidth,r,h,5.,,,"Sky width (pixels)"
+saturation,r,h,INDEF,,,"Saturation level"
+ignore_sat,b,h,no,,,"Ignore objects with saturated pixels?"
+iterations,i,h,2,1,,"Number of radius adjustment iterations"
+xcenter,r,h,INDEF,,,X field center (pixels)
+ycenter,r,h,INDEF,,,X field center (pixels)
+logfile,s,h,"logfile",,,"Logfile
+"
+imagecur,*imcur,h,"",,,"Image cursor input"
+graphcur,*gcur,h,"",,,"Graphics cursor input"
diff --git a/noao/obsutil/src/starfocus/starfocus.h b/noao/obsutil/src/starfocus/starfocus.h
new file mode 100644
index 00000000..871e8c98
--- /dev/null
+++ b/noao/obsutil/src/starfocus/starfocus.h
@@ -0,0 +1,140 @@
+# STARFOCUS
+
+# Types of coordinates
+define SF_TYPES "|center|mark1|markall|"
+define SF_CENTER 1 # Star at center of image
+define SF_MARK1 2 # Mark stars in first image
+define SF_MARKALL 3 # Mark stars in all images
+
+# Task type
+define STARFOCUS 1
+define PSFMEASURE 2
+
+# Radius types
+define SF_WTYPES "|Radius|FWHM|GFWHM|MFWHM|"
+
+define SF_RMIN 16 # Minimum centering search radius
+define MAX_FRAMES 8 # Maximum number of display frames
+
+# Data structures for STARFOCUS
+
+define NBNDRYPIX 100 # Number of boundary pixels
+define TYBNDRY BT_REFLECT # Type of boundary extension
+define SAMPLE .2 # Subpixel sampling size
+define SF_SZFNAME 79 # Length of file names
+define SF_SZWTYPE 7 # Length of width type string
+
+# Main data structure
+define SF 40
+define SF_TASK Memi[$1] # Task type
+define SF_WTYPE Memc[P2C($1+1)] # Width type string
+define SF_WCODE Memi[$1+5] # Width code
+define SF_BETA Memr[P2R($1+6)] # Moffat beta
+define SF_SCALE Memr[P2R($1+7)] # Pixel scale
+define SF_LEVEL Memr[P2R($1+8)] # Profile measurement level
+define SF_RADIUS Memr[P2R($1+9)] # Profile radius
+define SF_SBUF Memr[P2R($1+10)]# Sky region buffer
+define SF_SWIDTH Memr[P2R($1+11)]# Sky region width
+define SF_SAT Memr[P2R($1+12)]# Saturation
+define SF_NIT Memi[$1+13] # Number of iterations for radius
+define SF_OVRPLT Memi[$1+14] # Overplot the best profile?
+define SF_NCOLS Memi[$1+15] # Number of image columns
+define SF_NLINES Memi[$1+16] # Number of image lines
+define SF_XF Memr[P2R($1+17)]# X field center
+define SF_YF Memr[P2R($1+18)]# Y field center
+define SF_GP Memi[$1+19] # GIO pointer
+define SF_F Memr[P2R($1+20)]# Best focus
+define SF_W Memr[P2R($1+21)]# Width at best focus
+define SF_M Memr[P2R($1+22)]# Brightest star magnitude
+define SF_XP1 Memr[P2R($1+23)]# First derivative point to plot
+define SF_XP2 Memr[P2R($1+24)]# Last derivative point to plot
+define SF_YP1 Memr[P2R($1+25)]# Minimum of derivative profile
+define SF_YP2 Memr[P2R($1+26)]# Maximum of derivative profile
+define SF_N Memi[$1+27] # Number of points not deleted
+define SF_NSFD Memi[$1+28] # Number of data points
+define SF_SFDS Memi[$1+29] # Pointer to data structures
+define SF_NS Memi[$1+30] # Number of stars not deleted
+define SF_NSTARS Memi[$1+31] # Number of stars
+define SF_STARS Memi[$1+32] # Pointer to star groups
+define SF_NF Memi[$1+33] # Number of focuses not deleted
+define SF_NFOCUS Memi[$1+34] # Number of different focus values
+define SF_FOCUS Memi[$1+35] # Pointer to focus groups
+define SF_NI Memi[$1+36] # Number of images not deleted
+define SF_NIMAGES Memi[$1+37] # Number of images
+define SF_IMAGES Memi[$1+38] # Pointer to image groups
+define SF_BEST Memi[$1+39] # Pointer to best focus star
+
+define SF_SFD Memi[SF_SFDS($1)+$2-1]
+define SF_SFS Memi[SF_STARS($1)+$2-1]
+define SF_SFF Memi[SF_FOCUS($1)+$2-1]
+define SF_SFI Memi[SF_IMAGES($1)+$2-1]
+
+# Basic data structure.
+define SFD 94
+define SFD_IMAGE Memc[P2C($1)] # Image name
+define SFD_DATA Memi[$1+40] # Pointer to real image raster
+define SFD_RADIUS Memr[P2R($1+41)]# Profile radius
+define SFD_NP Memi[$1+42] # Number of profile points
+define SFD_NPMAX Memi[$1+43] # Maximum number of profile points
+define SFD_X1 Memi[$1+44] # Image raster limits
+define SFD_X2 Memi[$1+45]
+define SFD_Y1 Memi[$1+46]
+define SFD_Y2 Memi[$1+47]
+define SFD_ID Memi[$1+48] # Star ID
+define SFD_X Memr[P2R($1+49)]# Star X position
+define SFD_Y Memr[P2R($1+50)]# Star Y position
+define SFD_F Memr[P2R($1+51)]# Focus
+define SFD_W Memr[P2R($1+52)]# Width to use
+define SFD_M Memr[P2R($1+53)]# Magnitude
+define SFD_E Memr[P2R($1+54)]# Ellipticity
+define SFD_PA Memr[P2R($1+55)]# Position angle
+define SFD_R Memr[P2R($1+56)]# Radius at given level
+define SFD_DFWHM Memr[P2R($1+57)]# Direct FWHM
+define SFD_GFWHM Memr[P2R($1+58)]# Gaussian FWHM
+define SFD_MFWHM Memr[P2R($1+59)]# Moffat FWHM
+define SFD_ASI1 Memi[$1+60] # Pointer to enclosed flux profile
+define SFD_ASI2 Memi[$1+61] # Pointer to derivative profile
+define SFD_YP1 Memr[P2R($1+62)]# Minimum of derivative profile
+define SFD_YP2 Memr[P2R($1+63)]# Maximum of derivative profile
+define SFD_FWHM Memr[P2R($1+$2+63)] # FWHM vs level=0.5*i (i=1-19)
+define SFD_BKGD Memr[P2R($1+83)]# Background value
+define SFD_BKGD1 Memr[P2R($1+84)]# Original background value
+define SFD_MISO Memr[P2R($1+85)]# Moment isophote
+define SFD_SIGMA Memr[P2R($1+86)]# Moffat alpha
+define SFD_ALPHA Memr[P2R($1+87)]# Moffat alpha
+define SFD_BETA Memr[P2R($1+88)]# Moffat beta
+define SFD_STATUS Memi[$1+89] # Status
+define SFD_NSAT Memi[$1+90] # Number of saturated pixels
+define SFD_SFS Memi[$1+91] # Pointer to star group
+define SFD_SFF Memi[$1+92] # Pointer to focus group
+define SFD_SFI Memi[$1+93] # Pointer to image group
+
+
+# Structure grouping data by star.
+define SFS ($1+7)
+define SFS_ID Memi[$1] # Star ID
+define SFS_F Memr[P2R($1+1)] # Best focus
+define SFS_W Memr[P2R($1+2)] # Best width
+define SFS_M Memr[P2R($1+3)] # Average magnitude
+define SFS_N Memi[$1+4] # Number of points used
+define SFS_NF Memi[$1+5] # Number of focuses
+define SFS_NSFD Memi[$1+6] # Number of data points
+define SFS_SFD Memi[$1+$2+6] # Array of data structures
+
+
+# Structure grouping stars by focus values.
+define SFF ($1+5)
+define SFF_F Memr[P2R($1)] # Focus
+define SFF_W Memr[P2R($1+1)] # Average width
+define SFF_N Memi[$1+2] # Number in average
+define SFF_NI Memi[$1+3] # Number of images
+define SFF_NSFD Memi[$1+4] # Number of data points
+define SFF_SFD Memi[$1+$2+4] # Array of data structures
+
+
+# Structure grouping stars by image.
+define SFI ($1+42)
+define SFI_IMAGE Memc[P2C($1)] # Image
+define SFI_N Memi[$1+40] # Number in imagE
+define SFI_NSFD Memi[$1+41] # Number of data points
+define SFI_SFD Memi[$1+$2+41] # Array of data structures
diff --git a/noao/obsutil/src/starfocus/starfocus.key b/noao/obsutil/src/starfocus/starfocus.key
new file mode 100644
index 00000000..119e74a6
--- /dev/null
+++ b/noao/obsutil/src/starfocus/starfocus.key
@@ -0,0 +1,15 @@
+ COMMAND OPTIONS
+
+? Page this help information.
+g Measure object and graph the results.
+m Measure object.
+q Quit object marking and go to next image.
+ At the end of all images go to analysis of all measurements.
+
+:show Show the current results.
+
+When an object is measured the center and enclosed flux profile is determined.
+When using the "mark1" option typing 'q' will measure all remaining images
+at the same cursor positions while the "markall" option goes to the
+next image until the image list is finished. If measuring only one
+object with the 'g' key then a 'q' will exit the program.
diff --git a/noao/obsutil/src/starfocus/starfocus.par b/noao/obsutil/src/starfocus/starfocus.par
new file mode 100644
index 00000000..b1b1d5b0
--- /dev/null
+++ b/noao/obsutil/src/starfocus/starfocus.par
@@ -0,0 +1,32 @@
+# STARFOCUS -- Determine best star focus.
+
+images,s,a,,,,"List of images"
+focus,s,h,"1x1",,,"Focus values"
+fstep,s,h,"",,,"Focus step
+"
+nexposures,s,h,"1",,,"Number of exposures"
+step,s,h,"30",,,"Step in pixels"
+direction,s,h,"-line","-line|+line|-column|+column",,"Step direction"
+gap,s,h,"end","none|beginning|end",,"Double step gap
+"
+coords,s,h,"mark1","center|mark1|markall",,"Object coordinates"
+wcs,s,h,"logical","logical|physical|world",,"Coordinate system"
+display,b,h,yes,,,"Display images?"
+frame,i,h,1,1,,"Display frame to use
+"
+level,r,h,0.5,,,"Measurement level (fraction or percent)"
+size,s,h,"FWHM","Radius|FWHM|GFWHM|MFWHM",,"Size to display"
+beta,r,h,INDEF,2.1,,Moffat beta parameter
+scale,r,h,1.,,,"Pixel scale"
+radius,r,h,5.,,,"Measurement radius (pixels)"
+sbuffer,r,h,5,,,"Sky buffer (pixels)"
+swidth,r,h,5.,,,"Sky width (pixels)"
+saturation,r,h,INDEF,,,"Saturation level"
+ignore_sat,b,h,no,,,"Ignore objects with saturated pixels?"
+iterations,i,h,2,1,,"Number of radius adjustment iterations"
+xcenter,r,h,INDEF,,,X field center (pixels)
+ycenter,r,h,INDEF,,,Y field center (pixels)
+logfile,s,h,"logfile",,,"Logfile
+"
+imagecur,*imcur,h,"",,,"Image cursor input"
+graphcur,*gcur,h,"",,,"Graphics cursor input"
diff --git a/noao/obsutil/src/starfocus/stfgraph.x b/noao/obsutil/src/starfocus/stfgraph.x
new file mode 100644
index 00000000..890dfc81
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfgraph.x
@@ -0,0 +1,2682 @@
+include <error.h>
+include <gset.h>
+include <mach.h>
+include "starfocus.h"
+
+# Interactive help files. There is one for STARFOCUS and one for PSFMEASUE.
+define STFHELP "starfocus$stfhelp.key"
+define PSFHELP "starfocus$psfhelp.key"
+define PROMPT "Options"
+
+# View ports for all plots.
+define VX1 .15 # Minimum X viewport for left graph
+define VX2 .47 # Maximum X viewport for left graph
+define VX3 .63 # Minimum X viewport for right graph
+define VX4 .95 # Maximum X viewport for right graph
+define VY1 .10 # Minimum Y viewport for bottom graph
+define VY2 .44 # Minimum Y viewport for bottom graph
+define VY3 .54 # Minimum Y viewport for top graph
+define VY4 .88 # Maximum Y viewport for top graph
+
+# Miscellaneous graphics parameters.
+define NMAX 5 # Maximum number of samples for labeling
+define HLCOLOR 2 # Highlight color
+define HLWIDTH 4. # Highlight width
+define GM_MARK GM_CROSS # Point marker
+define GM_MAG GM_PLUS+GM_CROSS # Magnitude marker
+
+
+# STF_GRAPH -- Interactive graphing of results.
+
+procedure stf_graph (sf)
+
+pointer sf #I Starfocus structure
+
+real wx, wy, x, y, r2, r2min, fa[8]
+int i, j, ix, iy, nx, ny, wcs, key, pkey, skey, redraw, clgcur()
+pointer sp, sysidstr, title, cmd, gp, gopen()
+pointer sfd, sfs, sff, current, nearest
+
+data fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+ call smark (sp)
+ call salloc (sysidstr, SZ_LINE, TY_CHAR)
+ call salloc (title, SZ_LINE, TY_CHAR)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ # Set system id label
+ call sysid (Memc[sysidstr], SZ_LINE)
+
+ # Open graphics and enter interactive graphics loop
+ SF_GP(sf) = gopen ("stdgraph", NEW_FILE, STDGRAPH)
+ gp = SF_GP(sf)
+ wcs = 0
+ if (SF_NF(sf) > 1)
+ key = 'f'
+ else if (SF_NS(sf) > 1)
+ key = 'a'
+ else
+ key = 'z'
+ pkey = 0
+ skey = 1
+ current = SF_BEST(sf)
+ repeat {
+ switch (key) {
+ case 'q': # Quit
+ break
+ case '?': # Help
+ if (SF_TASK(sf) == PSFMEASURE)
+ call gpagefile (gp, PSFHELP, PROMPT)
+ else
+ call gpagefile (gp, STFHELP, PROMPT)
+ next
+ case ':': # Colon commands
+ iferr (call stf_colon (sf, Memc[cmd], redraw))
+ redraw = NO
+ if (redraw == NO)
+ next
+ case 'a', 'b', 'e', 'f', 'g', 'm', 'p', 't', 'z': # Plots
+ # When there is not enough data for the requested plot
+ # map the key to another one. This is done mostly to
+ # avoid redrawing the same graph when different keys
+ # map to the same pkey. The 'e', 'g', and 'p' key may
+ # select a different object so the check for the same
+ # plot is deferred.
+
+ if (SF_NS(sf) > 1 && SF_NF(sf) > 1) {
+ ;
+ } else if (SF_NS(sf) > 1) {
+ if (key == 'b')
+ key = 'a'
+ if (key == 'f')
+ key = 'm'
+ } else if (SF_NF(sf) > 1) {
+ if (key == 'a' || key == 'b' || key == 'm' || key == 't')
+ key = 'f'
+ } else {
+ key = 'z'
+ }
+
+ switch (key) {
+ case 'e', 'g', 'p':
+ ;
+ default:
+ if (key == pkey)
+ next
+ }
+ case 's': # Toggle plotting of magnitude symbols
+ if (pkey != 'a' && pkey != 'b')
+ next
+ skey = mod (skey+1, 2)
+ case 'u': # Undelete all
+ j = 0
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ if (SFD_STATUS(sfd) != 0) {
+ SFD_STATUS(sfd) = 0
+ j = j + 1
+ }
+ }
+ if (j == 0)
+ next
+ call stf_fitfocus (sf)
+ case 'd', 'n', 'o', 'r', 'x', 'i', ' ': # Misc
+ ;
+ default: # Unknown
+ call printf ("\007")
+ next
+ }
+
+ # Find the nearest or next object if needed.
+ switch (key) {
+ case 'r', 's', 'u', ':': # Redraw last graph
+ pkey = pkey
+ nearest = current
+ case 'n', 'o': # Renormalize enclosed flux profile
+ if (wcs != 7 || pkey == 'p')
+ next
+ pkey = pkey
+ nearest = current
+ if (key == 'n')
+ call stf_norm (sf, nearest, wx, INDEF)
+ else
+ call stf_norm (sf, nearest, wx, wy)
+ call stf_widths (sf, nearest)
+ call stf_fwhms (sf, nearest)
+ call stf_fitfocus (sf)
+ case ' ': # Select next focus or star
+ switch (pkey) {
+ case 'a', 'm', 't':
+ sff = SFD_SFF(current)
+ for (i=1; SF_SFF(sf,i)!=sff; i=i+1)
+ ;
+ j = SF_NFOCUS(sf)
+ i = mod (i, j) + 1
+ for (; SFF_N(SF_SFF(sf,i))==0; i=mod(i,j)+1)
+ ;
+ if (SF_SFF(sf,i) == sff)
+ next
+ sff = SF_SFF(sf,i)
+ do i = 1, SFF_NSFD(sff) {
+ nearest = SFF_SFD(sff,i)
+ if (SFD_STATUS(nearest) == 0)
+ break
+ }
+ case 'e', 'g', 'p', 'z':
+ switch (wcs) {
+ case 7, 8, 11:
+ for (i=1; SF_SFD(sf,i)!=current; i=i+1)
+ ;
+ j = SF_NSFD(sf)
+ i = mod (i, j) + 1
+ for (; SFD_STATUS(SF_SFD(sf,i))!=0; i=mod(i,j)+1)
+ ;
+ nearest = SF_SFD(sf,i)
+ case 9:
+ sfs = SFD_SFS(current)
+ for (i=1; SFS_SFD(sfs,i)!=current; i=i+1)
+ ;
+ j = SFS_NSFD(sfs)
+ i = mod (i, j) + 1
+ for (; SFD_STATUS(SFS_SFD(sfs,i))!=0; i=mod(i,j)+1)
+ ;
+ nearest = SFS_SFD(sfs,i)
+ if (nearest == current)
+ next
+ case 10:
+ sff = SFD_SFF(current)
+ for (i=1; SFF_SFD(sff,i)!=current; i=i+1)
+ ;
+ j = SFF_NSFD(sff)
+ i = mod (i, j) + 1
+ for (; SFD_STATUS(SFF_SFD(sff,i))!=0; i=mod(i,j)+1)
+ ;
+ nearest = SFF_SFD(sff,i)
+ if (nearest == current)
+ next
+ }
+ default:
+ next
+ }
+ default: # Select nearest to cursor
+ switch (pkey) {
+ case 'a':
+ r2min = MAX_REAL
+ call gctran (gp, wx, wy, wx, wy, wcs, 0)
+ sff = SFD_SFF(current)
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ switch (wcs) {
+ case 1:
+ x = SFD_X(sfd)
+ y = SFD_Y(sfd)
+ case 2:
+ x = SFD_X(sfd)
+ y = SFD_W(sfd)
+ case 3:
+ x = SFD_W(sfd)
+ y = SFD_Y(sfd)
+ case 4:
+ x = SFD_X(sfd)
+ y = SFD_E(sfd)
+ case 5:
+ x = SFD_E(sfd)
+ y = SFD_Y(sfd)
+ }
+ call gctran (gp, x, y, x, y, wcs, 0)
+ r2 = (x-wx)**2 + (y-wy)**2
+ if (r2 < r2min) {
+ r2min = r2
+ nearest = sfd
+ }
+ }
+ case 'b':
+ r2min = MAX_REAL
+ call gctran (gp, wx, wy, wx, wy, wcs, 0)
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_N(sfs) == 0)
+ next
+ switch (wcs) {
+ case 1:
+ x = SFD_X(SFS_SFD(sfs,1))
+ y = SFD_Y(SFS_SFD(sfs,1))
+ case 2:
+ x = SFD_X(SFS_SFD(sfs,1))
+ y = SFS_W(sfs)
+ case 3:
+ x = SFS_W(sfs)
+ y = SFD_Y(SFS_SFD(sfs,1))
+ case 4:
+ x = SFD_X(SFS_SFD(sfs,1))
+ y = SFS_F(sfs)
+ case 5:
+ x = SFS_F(sfs)
+ y = SFD_Y(SFS_SFD(sfs,1))
+ }
+ call gctran (gp, x, y, x, y, wcs, 0)
+ r2 = (x-wx)**2 + (y-wy)**2
+ if (r2 < r2min) {
+ r2min = r2
+ nearest = sfs
+ }
+ }
+ sfs = nearest
+ r2min = MAX_REAL
+ do i = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ r2 = SFD_W(sfd)
+ if (r2 < r2min) {
+ r2min = r2
+ nearest = sfd
+ }
+ }
+ case 'e', 'g', 'p':
+ switch (wcs) {
+ case 9:
+ sfs = SFD_SFS(current)
+ i = SFS_N(sfs)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+ ix = max (1, min (nx, nint(wx)))
+ iy = max (1, min (ny, nint(wy)))
+
+ j = 0
+ do i = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ if (ix == 1 + mod (j, nx) && iy == 1 + j / nx) {
+ nearest = sfd
+ break
+ }
+ j = j + 1
+ }
+ case 10:
+ sff = SFD_SFF(current)
+ i = SFF_N(sff)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+ ix = max (1, min (nx, nint(wx)))
+ iy = max (1, min (ny, nint(wy)))
+
+ j = 0
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ if (ix == 1 + mod (j, nx) && iy == 1 + j / nx) {
+ nearest = sfd
+ break
+ }
+ j = j + 1
+ }
+ }
+ if (key == pkey && nearest == current)
+ next
+ default:
+ switch (wcs) {
+ case 1, 2:
+ r2min = MAX_REAL
+ call gctran (gp, wx, wy, wx, wy, wcs, 0)
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ switch (wcs) {
+ case 1:
+ x = SFD_F(sfd)
+ y = SFD_W(sfd)
+ case 2:
+ x = SFD_F(sfd)
+ y = SFD_E(sfd)
+ }
+ call gctran (gp, x, y, x, y, wcs, 0)
+ r2 = (x-wx)**2 + (y-wy)**2
+ if (r2 < r2min) {
+ r2min = r2
+ nearest = sfd
+ }
+ }
+ case 3, 4, 5, 6:
+ r2min = MAX_REAL
+ call gctran (gp, wx, wy, wx, wy, wcs, 0)
+ sff = SFD_SFF(current)
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ switch (wcs) {
+ case 3:
+ x = -2.5 * log10 (SFS_M(SFD_SFS(sfd))/SF_M(sf))
+ y = SFD_W(sfd)
+ case 4:
+ x = -2.5 * log10 (SFS_M(SFD_SFS(sfd))/SF_M(sf))
+ y = SFD_E(sfd)
+ case 5:
+ x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+ (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+ y = SFD_W(sfd)
+ case 6:
+ x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+ (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+ y = SFD_E(sfd)
+ }
+ call gctran (gp, x, y, x, y, wcs, 0)
+ r2 = (x-wx)**2 + (y-wy)**2
+ if (r2 < r2min) {
+ r2min = r2
+ nearest = sfd
+ }
+ }
+ default:
+ nearest = current
+ }
+ }
+
+ # Act on selection for delete or info.
+ switch (key) {
+ case 'd':
+ if (SF_NS(sf) > 1) {
+ sfs = SFD_SFS(nearest)
+ do i = 1, SFS_NSFD(sfs)
+ SFD_STATUS(SFS_SFD(sfs,i)) = 1
+ } else
+ SFD_STATUS(nearest) = 1
+ call stf_fitfocus (sf)
+ case 'x':
+ repeat {
+ switch (key) {
+ case 'f':
+ sff = SFD_SFF(nearest)
+ do i = 1, SFF_NSFD(sff)
+ SFD_STATUS(SFF_SFD(sff,i)) = 1
+ case 'i':
+ sfd = SFD_SFI(nearest)
+ do i = 1, SFI_NSFD(sfd)
+ SFD_STATUS(SFI_SFD(sfd,i)) = 1
+ case 'p':
+ SFD_STATUS(nearest) = 1
+ case 's':
+ sfs = SFD_SFS(nearest)
+ do i = 1, SFS_NSFD(sfs)
+ SFD_STATUS(SFS_SFD(sfs,i)) = 1
+ default:
+ call printf (
+ "Delete image, star, focus, or point? (i|s|f|p)")
+ next
+ }
+ call stf_fitfocus (sf)
+ break
+ } until (clgcur ("graphcur",
+ wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+ case 'i':
+ switch (pkey) {
+ case 'b':
+ sfs = SFD_SFS(nearest)
+ call stf_title (sf, NULL, sfs, NULL,
+ Memc[title], SZ_LINE)
+ default:
+ call stf_title (sf, nearest, NULL, NULL,
+ Memc[title], SZ_LINE)
+ }
+ call printf ("%s\n")
+ call pargstr (Memc[title])
+ next
+ default:
+ pkey = key
+ }
+ }
+
+ # If current object has been deleted select another.
+ if (SFD_STATUS(nearest) == 0)
+ current = nearest
+ else
+ current = SF_BEST(sf)
+
+ # Make the graphs. The graph depends on the number of stars
+ # and number of focus values. Note that the pkey has already
+ # been mapped but all the keys are shown for clarity.
+
+ call gclear (gp)
+ call gseti (gp, G_FACOLOR, 0)
+
+ if (SF_NS(sf) > 1 && SF_NF(sf) > 1) {
+ switch (pkey) {
+ case 'a':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 1)
+ call gsview (gp, VX1, VX4, VY1, VY4)
+ call stf_g11 (sf, current, skey, Memc[title])
+ case 'b':
+ call sprintf (Memc[title], SZ_LINE,
+ "Best focus estimates for each star")
+ call gseti (gp, G_WCS, 1)
+ call gsview (gp, VX1, VX4, VY1, VY4)
+ call stf_g12 (sf, current, skey, Memc[title])
+ case 'e':
+ sfs = SFD_SFS(current)
+ call sprintf (Memc[title], SZ_LINE,
+ "Star: x=%.2f, y=%.2f, m=%.2f")
+ call pargr (SFD_X(current))
+ call pargr (SFD_Y(current))
+ call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+ call gseti (gp, G_WCS, 9)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g2 (sf, current, Memc[title])
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 10)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call stf_g3 (sf, current, Memc[title])
+ case 'f':
+ call gseti (gp, G_WCS, 1)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g1 (sf, current, 'f', 'r', "", "", SF_WTYPE(sf))
+ call gseti (gp, G_WCS, 2)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g1 (sf, current, 'f', 'e', "", "Focus",
+ "Ellipticity")
+ case 'g':
+ sfs = SFD_SFS(current)
+ call sprintf (Memc[title], SZ_LINE,
+ "Star: x=%.2f, y=%.2f, m=%.2f")
+ call pargr (SFD_X(current))
+ call pargr (SFD_Y(current))
+ call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+ call gseti (gp, G_WCS, 9)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g9 (sf, current, Memc[title])
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 10)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call stf_g10 (sf, current, Memc[title])
+ case 'm':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 3)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g1 (sf, current, 'm', 'r', Memc[title],
+ "", SF_WTYPE(sf))
+ call gseti (gp, G_WCS, 4)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g1 (sf, current, 'm', 'e', "", "Magnitude",
+ "Ellipticity")
+ case 'p':
+ sfs = SFD_SFS(current)
+ call sprintf (Memc[title], SZ_LINE,
+ "Star: x=%.2f, y=%.2f, m=%.2f")
+ call pargr (SFD_X(current))
+ call pargr (SFD_Y(current))
+ call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+ call gseti (gp, G_WCS, 9)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g4 (sf, current, Memc[title])
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 10)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call stf_g5 (sf, current, Memc[title])
+ case 't':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 5)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g1 (sf, current, 't', 'r', Memc[title],
+ "", SF_WTYPE(sf))
+ call gseti (gp, G_WCS, 6)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g1 (sf, current, 't', 'e', "", "Field radius",
+ "Ellipticity")
+ case 'z':
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX2, VY3, VY4)
+ call stf_g6 (sf, current, "", "", "Enclosed flux")
+ call gseti (gp, G_WCS, 8)
+ call gsview (gp, VX1, VX2, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g7 (sf, current, "", "Radius", "Profile")
+ call gseti (gp, G_WCS, 11)
+ call gsview (gp, VX3, VX4, VY3, VY4)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 0)
+ call gsetr (gp, G_PLWIDTH, 2.0)
+ call gline (gp, 0., 0., 0., 0.)
+ call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+ }
+ } else if (SF_NS(sf) > 1) {
+ switch (pkey) {
+ case 'a', 'b':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 1)
+ call gsview (gp, VX1, VX4, VY1, VY4)
+ call stf_g11 (sf, current, skey, Memc[title])
+ case 'e':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 10)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g3 (sf, current, Memc[title])
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g6 (sf, current, Memc[title], "Radius",
+ "Enclosed flux")
+ case 'f', 'm':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 3)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g1 (sf, current, 'm', 'r', Memc[title], "",
+ SF_WTYPE(sf))
+ call gseti (gp, G_WCS, 4)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g1 (sf, current, 'm', 'e', "", "Magnitude",
+ "Ellipticity")
+ case 'g':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 10)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g10 (sf, current, Memc[title])
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 11)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g8 (sf, current, Memc[title], "Enclosed flux",
+ "FWHM")
+ case 'p':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 10)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g5 (sf, current, Memc[title])
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g7 (sf, current, Memc[title], "Radius", "Profile")
+ case 't':
+ sff = SFD_SFF(current)
+ call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 5)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g1 (sf, current, 't', 'r', Memc[title], "",
+ SF_WTYPE(sf))
+ call gseti (gp, G_WCS, 6)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g1 (sf, current, 't', 'e', "", "Field radius",
+ "Ellipticity")
+ case 'z':
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX2, VY3, VY4)
+ call stf_g6 (sf, current, "", "", "Enclosed flux")
+ call gseti (gp, G_WCS, 8)
+ call gsview (gp, VX1, VX2, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g7 (sf, current, "", "Radius", "Profile")
+ call gseti (gp, G_WCS, 11)
+ call gsview (gp, VX3, VX4, VY3, VY4)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 0)
+ call gsetr (gp, G_PLWIDTH, 2.0)
+ call gline (gp, 0., 0., 0., 0.)
+ call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+ }
+ } else if (SF_NF(sf) > 1) {
+ switch (pkey) {
+ case 'a', 'b', 'f', 'm', 't':
+ call gseti (gp, G_WCS, 1)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g1 (sf, current, 'f', 'r', "", "", SF_WTYPE(sf))
+ call gseti (gp, G_WCS, 2)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g1 (sf, current, 'f', 'e', "", "Focus",
+ "Ellipticity")
+ case 'e':
+ sfs = SFD_SFS(current)
+ call sprintf (Memc[title], SZ_LINE,
+ "Star: x=%.2f, y=%.2f, m=%.2f")
+ call pargr (SFD_X(current))
+ call pargr (SFD_Y(current))
+ call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+ call gseti (gp, G_WCS, 9)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g2 (sf, current, Memc[title])
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g6 (sf, current, Memc[title], "Radius",
+ "Enclosed flux")
+ case 'g':
+ sfs = SFD_SFS(current)
+ call sprintf (Memc[title], SZ_LINE,
+ "Star: x=%.2f, y=%.2f, m=%.2f")
+ call pargr (SFD_X(current))
+ call pargr (SFD_Y(current))
+ call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+ call gseti (gp, G_WCS, 9)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g9 (sf, current, Memc[title])
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 11)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g8 (sf, current, Memc[title], "Enclosed flux",
+ "FWHM")
+ case 'p':
+ sfs = SFD_SFS(current)
+ call sprintf (Memc[title], SZ_LINE,
+ "Star: x=%.2f, y=%.2f, m=%.2f")
+ call pargr (SFD_X(current))
+ call pargr (SFD_Y(current))
+ call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+ call gseti (gp, G_WCS, 9)
+ call gsview (gp, VX1, VX4, VY3, VY4)
+ call stf_g4 (sf, current, Memc[title])
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX4, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g7 (sf, current, Memc[title], "Radius",
+ "profile")
+ case 'z':
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX2, VY3, VY4)
+ call stf_g6 (sf, current, "", "", "Enclosed flux")
+ call gseti (gp, G_WCS, 8)
+ call gsview (gp, VX1, VX2, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g7 (sf, current, "", "Radius", "Profile")
+ call gseti (gp, G_WCS, 11)
+ call gsview (gp, VX3, VX4, VY3, VY4)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 0)
+ call gsetr (gp, G_PLWIDTH, 2.0)
+ call gline (gp, 0., 0., 0., 0.)
+ call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+ }
+ } else {
+ switch (pkey) {
+ case 'a', 'b', 'f', 'm', 'p', 'z', 'e', 't':
+ call gseti (gp, G_WCS, 7)
+ call gsview (gp, VX1, VX2, VY3, VY4)
+ call stf_g6 (sf, current, "", "", "Enclosed flux")
+ call gseti (gp, G_WCS, 8)
+ call gsview (gp, VX1, VX2, VY1, VY2)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g7 (sf, current, "", "Radius", "Profile")
+ call gseti (gp, G_WCS, 11)
+ call gsview (gp, VX3, VX4, VY3, VY4)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+ call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+ call stf_title (sf, current, NULL, NULL, Memc[title],
+ SZ_LINE)
+ call gseti (gp, G_WCS, 0)
+ call gsetr (gp, G_PLWIDTH, 2.0)
+ call gline (gp, 0., 0., 0., 0.)
+ call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+ }
+ }
+
+ # Add banner title.
+ call stf_title (sf, NULL, NULL, NULL, Memc[title], SZ_LINE)
+ call gseti (gp, G_WCS, 0)
+ call gsetr (gp, G_PLWIDTH, 2.0)
+ call gline (gp, 0., 0., 0., 0.)
+ call gtext (gp, 0.5, 0.99, Memc[sysidstr], "h=c,v=t")
+ call gtext (gp, 0.5, 0.96, Memc[title], "h=c,v=t")
+
+ if (SF_NSFD(sf) == 1)
+ break
+
+ } until (clgcur ("graphcur", wx, wy, wcs, key, Memc[cmd], SZ_LINE)==EOF)
+
+ call gclose (gp)
+ call sfree (sp)
+end
+
+
+# List of colon commands.
+define CMDS "|show|level|size|scale|radius|xcenter|ycenter\
+ |overplot|beta|"
+define SHOW 1 # Show current results
+define LEVEL 2 # Measurement level
+define SIZE 3 # Size type
+define SCALE 4 # Pixel scale
+define RADIUS 5 # Maximum radius
+define XCENTER 6 # X field center
+define YCENTER 7 # Y field center
+define OVERPLOT 8 # Overplot best profile
+define BETA 9 # Beta value for Moffat function
+
+# STF_COLON -- Respond to colon command.
+
+procedure stf_colon (sf, cmd, redraw)
+
+pointer sf #I Starfocus pointer
+char cmd[ARB] #I Colon command
+int redraw #O Redraw?
+
+bool bval
+real rval, stf_r2i()
+int i, j, ncmd, nscan(), strdic(), open(), btoi()
+pointer sp, str, sfd
+errchk open, delete, stf_log, stf_norm, stf_radius, stf_fitfocus
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Scan the command string and get the first word.
+ call sscan (cmd)
+ call gargwrd (Memc[str], SZ_FNAME)
+ ncmd = strdic (Memc[str], Memc[str], SZ_FNAME, CMDS)
+
+ switch (ncmd) {
+ case SHOW:
+ call gargwrd (Memc[str], SZ_FNAME)
+ iferr {
+ if (nscan() == 1) {
+ call mktemp ("tmp$iraf", Memc[str], SZ_FNAME)
+ i = open (Memc[str], APPEND, TEXT_FILE)
+ call stf_log (sf, i)
+ call close (i)
+ call gpagefile (SF_GP(sf), Memc[str], "starfocus")
+ call delete (Memc[str])
+ } else {
+ i = open (Memc[str], APPEND, TEXT_FILE)
+ call stf_log (sf, i)
+ call close (i)
+ }
+ } then
+ call erract (EA_WARN)
+ redraw = NO
+ case LEVEL:
+ call gargr (rval)
+ if (nscan() == 2) {
+ if (rval > 1.)
+ rval = rval / 100.
+ SF_LEVEL(sf) = max (0.05, min (0.95, rval))
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ call stf_radius (sf, sfd, SF_LEVEL(sf), SFD_R(sfd))
+ }
+ if (SF_WCODE(sf) == 1)
+ call stf_fitfocus (sf)
+ redraw = YES
+ } else {
+ call printf ("level %g\n")
+ call pargr (SF_LEVEL(sf))
+ redraw = NO
+ }
+ case SIZE:
+ call gargwrd (Memc[str], SZ_FNAME)
+ if (nscan() == 2) {
+ ncmd = strdic (Memc[str], Memc[str], SZ_FNAME, SF_WTYPES)
+ if (ncmd == 0) {
+ call eprintf ("Invalid size type\n")
+ redraw = NO
+ } else {
+ call strcpy (Memc[str], SF_WTYPE(sf), SF_SZWTYPE)
+ SF_WCODE(sf) = ncmd
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ switch (SF_WCODE(sf)) {
+ case 1:
+ SFD_W(sfd) = SFD_R(sfd)
+ case 2:
+ SFD_W(sfd) = SFD_DFWHM(sfd)
+ case 3:
+ SFD_W(sfd) = SFD_GFWHM(sfd)
+ case 4:
+ SFD_W(sfd) = SFD_MFWHM(sfd)
+ }
+ call stf_fwhms (sf, sfd)
+ }
+ call stf_fitfocus (sf)
+ redraw = YES
+ }
+ } else {
+ call printf ("size %s\n")
+ call pargstr (SF_WTYPE(sf))
+ redraw = NO
+ }
+ case SCALE:
+ call gargr (rval)
+ if (nscan() == 2) {
+ rval = rval / SF_SCALE(sf)
+ SF_SCALE(sf) = SF_SCALE(sf) * rval
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ switch (SF_WCODE(sf)) {
+ case 1:
+ SFD_R(sfd) = SFD_R(sfd) * rval
+ SFD_W(sfd) = SFD_R(sfd)
+ case 2:
+ SFD_DFWHM(sfd) = SFD_DFWHM(sfd) * rval
+ SFD_W(sfd) = SFD_DFWHM(sfd)
+ case 3:
+ SFD_SIGMA(sfd) = SFD_SIGMA(sfd) * rval
+ SFD_GFWHM(sfd) = SFD_GFWHM(sfd) * rval
+ SFD_W(sfd) = SFD_GFWHM(sfd)
+ case 4:
+ SFD_ALPHA(sfd) = SFD_ALPHA(sfd) * rval
+ SFD_MFWHM(sfd) = SFD_MFWHM(sfd) * rval
+ SFD_W(sfd) = SFD_MFWHM(sfd)
+ }
+ do j = 1, 19
+ SFD_FWHM(sfd,j) = SFD_FWHM(sfd,j) * rval
+ }
+ do i = 1, SF_NSTARS(sf) {
+ sfd = SF_SFS(sf,i)
+ SFS_W(sfd) = SFS_W(sfd) * rval
+ }
+ do i = 1, SF_NFOCUS(sf) {
+ sfd = SF_SFF(sf,i)
+ SFF_W(sfd) = SFF_W(sfd) * rval
+ }
+ SF_W(sf) = SF_W(sf) * rval
+ redraw = YES
+ } else {
+ call printf ("scale %g\n")
+ call pargr (SF_SCALE(sf))
+ redraw = NO
+ }
+ case RADIUS:
+ call gargr (rval)
+ if (nscan() == 2) {
+ j = stf_r2i (rval) + 1
+ SF_RADIUS(sf) = rval
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ if (j > SFD_NPMAX(sfd))
+ next
+ SFD_NP(sfd) = j
+ SFD_RADIUS(sf) = SF_RADIUS(sf)
+ call stf_norm (sf, sfd, INDEF, INDEF)
+ call stf_widths (sf, sfd)
+ call stf_fwhms (sf, sfd)
+ }
+ call stf_fitfocus (sf)
+ redraw = YES
+ } else {
+ call printf ("radius %g\n")
+ call pargr (SF_RADIUS(sf))
+ redraw = NO
+ }
+ case XCENTER:
+ call gargr (rval)
+ if (nscan() == 2) {
+ if (IS_INDEF(rval))
+ SF_XF(sf) = (SF_NCOLS(sf) + 1) / 2.
+ else
+ SF_XF(sf) = rval
+ redraw = NO
+ } else {
+ call printf ("xcenter %g\n")
+ call pargr (SF_XF(sf))
+ redraw = NO
+ }
+ case YCENTER:
+ call gargr (rval)
+ if (nscan() == 2) {
+ if (IS_INDEF(rval))
+ SF_YF(sf) = (SF_NLINES(sf) + 1) / 2.
+ else
+ SF_YF(sf) = rval
+ redraw = NO
+ } else {
+ call printf ("ycenter %g\n")
+ call pargr (SF_YF(sf))
+ redraw = NO
+ }
+ case OVERPLOT:
+ call gargb (bval)
+ if (nscan() == 2) {
+ SF_OVRPLT(sf) = btoi (bval)
+ redraw = YES
+ } else {
+ call printf ("overplot %b\n")
+ call pargi (SF_OVRPLT(sf))
+ redraw = NO
+ }
+ case BETA:
+ call gargr (rval)
+ if (nscan() == 2) {
+ SF_BETA(sf) = rval
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ call stf_widths (sf, sfd)
+ switch (SF_WCODE(sf)) {
+ case 1:
+ SFD_W(sfd) = SFD_R(sfd)
+ case 2:
+ SFD_W(sfd) = SFD_DFWHM(sfd)
+ case 3:
+ SFD_W(sfd) = SFD_GFWHM(sfd)
+ case 4:
+ SFD_W(sfd) = SFD_MFWHM(sfd)
+ }
+ call stf_fwhms (sf, sfd)
+ }
+ call stf_fitfocus (sf)
+ redraw = YES
+ } else {
+ call printf ("beta %g\n")
+ call pargr (SF_BETA(sf))
+ redraw = NO
+ }
+ default:
+ call printf ("Unrecognized or ambiguous command\007")
+ redraw = NO
+ }
+
+ call sfree (sp)
+end
+
+
+# STF_G1 -- Plot of size/ellip vs. focus/mag/radius.
+
+procedure stf_g1 (sf, current, xkey, ykey, title, xlabel, ylabel)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+int xkey #I X axis key
+int ykey #I Y axis key
+char title[ARB] #I Title
+char xlabel[ARB] #I X label
+char ylabel[ARB] #I Y label
+
+int i, j
+bool hl
+real x, x1, x2, dx, y, y1, y2, dy
+pointer gp, sff, sfd
+
+begin
+ # Determine data range
+ x1 = MAX_REAL
+ x2 = -MAX_REAL
+ switch (ykey) {
+ case 'r':
+ y1 = SF_W(sf)
+ y2 = 1.5 * SF_W(sf)
+ case 'e':
+ y1 = 0
+ y2 = 1
+ }
+ do j = 1, SF_NFOCUS(sf) {
+ sff = SF_SFF(sf,j)
+ if (xkey != 'f' && sff != SFD_SFF(current))
+ next
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) == 0) {
+ switch (xkey) {
+ case 'f':
+ x = SFD_F(sfd)
+ case 'm':
+ x = -2.5 * log10 (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+ case 't':
+ x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+ (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+ }
+ switch (ykey) {
+ case 'r':
+ y = SFD_W(sfd)
+ case 'e':
+ y = SFD_E(sfd)
+ }
+ x1 = min (x1, x)
+ x2 = max (x2, x)
+ y1 = min (y1, y)
+ y2 = max (y2, y)
+ }
+ }
+ }
+
+ dx = (x2 - x1)
+ dy = (y2 - y1)
+ x1 = x1 - dx * 0.05
+ x2 = x2 + dx * 0.05
+ y1 = y1 - dy * 0.05
+ y2 = y2 + dy * 0.05
+ gp = SF_GP (sf)
+ call gswind (gp, x1, x2, y1, y2)
+ call glabax (gp, title, xlabel, ylabel)
+
+ do j = 1, SF_NFOCUS(sf) {
+ sff = SF_SFF(sf,j)
+ if (xkey != 'f' && sff != SFD_SFF(current))
+ next
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) == 0) {
+ hl = false
+ switch (xkey) {
+ case 'f':
+ x = SFD_F(sfd)
+ #hl = (SFD_SFS(sfd) == SFD_SFS(current))
+ case 'm':
+ x = -2.5 * log10 (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+ #hl = (SFD_SFF(sfd) != SFD_SFF(current))
+ case 't':
+ x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+ (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+ #hl = (SFD_SFF(sfd) != SFD_SFF(current))
+ }
+ switch (ykey) {
+ case 'r':
+ y = SFD_W(sfd)
+ case 'e':
+ y = SFD_E(sfd)
+ }
+ if (hl) {
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ if (sfd == current)
+ call gmark (gp, x, y, GM_BOX, 3., 3.)
+ call gmark (gp, x, y, GM_PLUS, 3., 3.)
+ call gseti (gp, G_PLCOLOR, 1)
+ } else
+ call gmark (gp, x, y, GM_MARK, 2., 2.)
+ }
+ }
+ }
+
+ call gseti (gp, G_PLTYPE, 2)
+ if (xkey == 'f')
+ call gline (gp, SF_F(sf), y1, SF_F(sf), y2)
+ if (ykey == 'r')
+ call gline (gp, x1, SF_W(sf), x2, SF_W(sf))
+ call gseti (gp, G_PLTYPE, 1)
+end
+
+
+# STF_G2 -- Enclosed flux profiles for a given star.
+
+procedure stf_g2 (sf, current, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+char title[ARB] #I Title
+
+int i, j, np, np1, nx, ny, ix, iy
+real vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, fa[10]
+pointer sp, str, gp, sfs, sfd, asi
+real stf_i2r(), stf_r2i(), asieval()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ gp = SF_GP(sf)
+ sfs = SFD_SFS(current)
+ np = SFD_NP(current)
+
+ # Set grid layout
+ i = SFS_N(sfs)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+
+ # Set subview port parameters
+ call ggview (gp, vx, dvx, vy, dvy)
+ dvx = (dvx - vx) / nx
+ dvy = (dvy - vy) / ny
+
+ # Set data window parameters
+ x1 = -0.05
+ x2 = 1.05
+ y1 = -0.15
+ y2 = 1.05
+ call gswind (gp, x1, x2, y1, y2)
+
+ # Set fill area
+ fa[1] = x1; fa[6] = y1
+ fa[2] = x2; fa[7] = y1
+ fa[3] = x2; fa[8] = y2
+ fa[4] = x1; fa[9] = y2
+ fa[5] = x1; fa[10] = y1
+
+ # Draw profiles.
+ j = 0
+ do i = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ np1 = SFD_NP(sfd)
+ ix = 1 + mod (j, nx)
+ iy = 1 + j / nx
+ j = j + 1
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ call gfill (gp, fa, fa[6], 4, GF_SOLID)
+ call gseti (gp, G_DRAWTICKS, NO)
+ call glabax (gp, "", "", "")
+ if (sfd == current) {
+ call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+ vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+ call gsetr (gp, G_PLWIDTH, HLWIDTH)
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ call gpline (gp, fa, fa[6], 5)
+ call gsetr (gp, G_PLWIDTH, 1.)
+ call gseti (gp, G_PLCOLOR, 1)
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ }
+
+ asi = SFD_ASI1(sfd)
+ r2 = stf_i2r (real(np))
+ call gamove (gp, 0., 0.)
+ for (z = 1.; z <= np1; z = z + 0.1)
+ call gadraw (gp, stf_i2r(z)/r2, asieval(asi,z))
+ if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ np1 = SFD_NP(SF_BEST(sf))
+ asi = SFD_ASI1(SF_BEST(sf))
+ r1 = stf_i2r (1.)
+ r2 = stf_i2r (real(np))
+ dr = 0.05 * (r2 - r1)
+ for (r = r1; r <= r2; r = r + dr) {
+ z = stf_r2i (r)
+ z1 = stf_r2i (r+0.7*dr)
+ if (z > 1. && z1 <= np1)
+ call gline (gp, r/r2, asieval(asi,z),
+ (r+0.7*dr)/r2, asieval(asi,z1))
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call sprintf (Memc[str], SZ_LINE, "%.3g")
+ call pargr (SFD_W(sfd))
+ call gtext (gp, 0.95, -0.1, Memc[str], "h=r;v=b")
+ if (nx < NMAX && ny < NMAX) {
+ call sprintf (Memc[str], SZ_LINE, "%.4g")
+ call pargr (SFD_F(sfd))
+ call gtext (gp, 0.05, -0.1, Memc[str], "h=l;v=b")
+ }
+ }
+
+ call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+ call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+ call gamove (gp, 1., 1.)
+
+ # Draw label
+ call gseti (gp, G_DRAWAXES, 0)
+ call glabax (gp, title, "", "")
+ call gseti (gp, G_DRAWAXES, 3)
+
+ call sfree (sp)
+end
+
+
+# STF_G3 -- Enclosed flux profiles for a given focus.
+
+procedure stf_g3 (sf, current, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+char title[ARB] #I Title
+
+int i, j, np, np1, nx, ny, ix, iy
+real vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, fa[10]
+pointer sp, str, gp, sff, sfd, asi
+real stf_i2r(), stf_r2i(), asieval()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ gp = SF_GP(sf)
+ sff = SFD_SFF(current)
+ np = SFD_NP(current)
+
+ # Set grid layout
+ i = SFF_N(sff)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+
+ # Set subview port parameters
+ call ggview (gp, vx, dvx, vy, dvy)
+ dvx = (dvx - vx) / nx
+ dvy = (dvy - vy) / ny
+
+ # Set data window parameters
+ x1 = -0.05
+ x2 = 1.05
+ y1 = -0.2
+ y2 = 1.05
+ call gswind (gp, x1, x2, y1, y2)
+
+ # Set fill area
+ fa[1] = x1; fa[6] = y1
+ fa[2] = x2; fa[7] = y1
+ fa[3] = x2; fa[8] = y2
+ fa[4] = x1; fa[9] = y2
+ fa[5] = x1; fa[10] = y1
+
+ # Draw profiles.
+ j = 0
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ np1 = SFD_NP(sfd)
+ ix = 1 + mod (j, nx)
+ iy = 1 + j / nx
+ j = j + 1
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ call gfill (gp, fa, fa[6], 4, GF_SOLID)
+ call gseti (gp, G_DRAWTICKS, NO)
+ call glabax (gp, "", "", "")
+ if (sfd == current) {
+ call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+ vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+ call gsetr (gp, G_PLWIDTH, HLWIDTH)
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ call gpline (gp, fa, fa[6], 5)
+ call gsetr (gp, G_PLWIDTH, 1.)
+ call gseti (gp, G_PLCOLOR, 1)
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ }
+
+ asi = SFD_ASI1(sfd)
+ r2 = stf_i2r (real(np))
+ call gamove (gp, 0., 0.)
+ for (z = 1.; z <= np1; z = z + 0.1)
+ call gadraw (gp, stf_i2r(z)/r2, asieval(asi,z))
+ if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ np1 = SFD_NP(SF_BEST(sf))
+ asi = SFD_ASI1(SF_BEST(sf))
+ r1 = stf_i2r (1.)
+ r2 = stf_i2r (real(np))
+ dr = 0.05 * (r2 - r1)
+ for (r = r1; r <= r2; r = r + dr) {
+ z = stf_r2i (r)
+ z1 = stf_r2i (r+0.7*dr)
+ if (z > 1. && z1 <= np1)
+ call gline (gp, r/r2, asieval(asi,z),
+ (r+0.7*dr)/r2, asieval(asi,z1))
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call sprintf (Memc[str], SZ_LINE, "%.3g")
+ call pargr (SFD_W(sfd))
+ call gtext (gp, 0.95, -.1, Memc[str], "h=r;v=b")
+ if (nx < NMAX && ny < NMAX) {
+ call sprintf (Memc[str], SZ_LINE, "%d %d")
+ call pargr (SFD_X(sfd))
+ call pargr (SFD_Y(sfd))
+ call gtext (gp, 0.05, -.1, Memc[str], "h=l;v=b")
+ }
+ }
+
+ call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+ call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+ call gamove (gp, 1., 1.)
+
+ # Draw label
+ call gseti (gp, G_DRAWAXES, 0)
+ call glabax (gp, title, "", "")
+ call gseti (gp, G_DRAWAXES, 3)
+
+ call sfree (sp)
+end
+
+
+# STF_G4 -- Radial profiles (derivative of enclosed flux) for a given star.
+
+procedure stf_g4 (sf, current, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+char title[ARB] #I Title
+
+int i, j, np, np1, nx, ny, ix, iy
+real vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, rmax, fa[10]
+pointer sp, str, gp, sfs, sfd, asi
+real stf_i2r(), stf_r2i(), asieval()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ gp = SF_GP(sf)
+ sfs = SFD_SFS(current)
+ np = SFD_NP(current)
+
+ # Set grid layout
+ i = SFS_N(sfs)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+
+ # Set subview port parameters
+ call ggview (gp, vx, dvx, vy, dvy)
+ dvx = (dvx - vx) / nx
+ dvy = (dvy - vy) / ny
+
+ # Set data window parameters
+ x1 = -0.05
+ x2 = 1.05
+ z = SF_YP2(sf) - SF_YP1(sf)
+ y1 = SF_YP1(sf) - 0.05 * z
+ y2 = SF_YP2(sf) + 0.15 * z
+
+ # Set fill area
+ fa[1] = x1; fa[6] = y1
+ fa[2] = x2; fa[7] = y1
+ fa[3] = x2; fa[8] = y2
+ fa[4] = x1; fa[9] = y2
+ fa[5] = x1; fa[10] = y1
+
+ # Draw profiles.
+ j = 0
+ do i = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ np1 = SFD_NP(sfd)
+ ix = 1 + mod (j, nx)
+ iy = 1 + j / nx
+ j = j + 1
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ call gswind (gp, x1, x2, y1, y2)
+ call gfill (gp, fa, fa[6], 4, GF_SOLID)
+ call gseti (gp, G_DRAWTICKS, NO)
+ call glabax (gp, "", "", "")
+ if (sfd == current) {
+ call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+ vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+ call gsetr (gp, G_PLWIDTH, HLWIDTH)
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ call gpline (gp, fa, fa[6], 5)
+ call gsetr (gp, G_PLWIDTH, 1.)
+ call gseti (gp, G_PLCOLOR, 1)
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ }
+
+ asi = SFD_ASI2(sfd)
+ rmax = stf_i2r (real(np))
+ z = SF_XP1(sf)
+ call gamove (gp, stf_i2r(z)/rmax, asieval(asi,z))
+ for (; z <= SF_XP2(sf); z = z + 0.1)
+ call gadraw (gp, stf_i2r(z)/rmax, asieval(asi,z))
+ if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ np1 = SFD_NP(SF_BEST(sf))
+ asi = SFD_ASI2(SF_BEST(sf))
+ rmax = stf_i2r (real(np))
+ r1 = stf_i2r (SF_XP1(sf))
+ r2 = stf_i2r (SF_XP2(sf))
+ dr = 0.05 * (rmax - stf_i2r(1.))
+ for (r = r1; r <= r2; r = r + dr) {
+ z = stf_r2i (r)
+ z1 = stf_r2i (r+0.7*dr)
+ if (z > 1. && z1 <= np1)
+ call gline (gp, r/rmax, asieval(asi,z),
+ (r+0.7*dr)/rmax, asieval(asi,z1))
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gswind (gp, 0., 1., 0., 1.)
+ call sprintf (Memc[str], SZ_LINE, "%.3g")
+ call pargr (SFD_W(sfd))
+ call gtext (gp, 0.95, 0.98, Memc[str], "h=r;v=t")
+ if (nx < NMAX && ny < NMAX) {
+ call sprintf (Memc[str], SZ_LINE, "%.4g")
+ call pargr (SFD_F(sfd))
+ call gtext (gp, 0.05, 0.98, Memc[str], "h=l;v=t")
+ }
+ }
+
+ call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+ call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+ call gamove (gp, 1., 1.)
+
+ # Draw label
+ call gseti (gp, G_DRAWAXES, 0)
+ call glabax (gp, title, "", "")
+ call gseti (gp, G_DRAWAXES, 3)
+
+ call sfree (sp)
+end
+
+
+# STF_G5 -- Radial profiles (derivative of enclosed flux) for a given focus.
+
+procedure stf_g5 (sf, current, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+char title[ARB] #I Title
+
+int i, j, np, np1, nx, ny, ix, iy
+real vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, rmax, fa[10]
+pointer sp, str, gp, sff, sfd, asi
+real stf_i2r(), stf_r2i(), asieval()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ gp = SF_GP(sf)
+ sff = SFD_SFF(current)
+ np = SFD_NP(current)
+
+ # Set grid layout
+ i = SFF_N(sff)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+
+ # Set subview port parameters
+ call ggview (gp, vx, dvx, vy, dvy)
+ dvx = (dvx - vx) / nx
+ dvy = (dvy - vy) / ny
+
+ # Set data window parameters
+ x1 = -0.05
+ x2 = 1.05
+ z = SF_YP2(sf) - SF_YP1(sf)
+ y1 = SF_YP1(sf) - 0.05 * z
+ y2 = SF_YP2(sf) + 0.15 * z
+
+ # Set fill area
+ fa[1] = x1; fa[6] = y1
+ fa[2] = x2; fa[7] = y1
+ fa[3] = x2; fa[8] = y2
+ fa[4] = x1; fa[9] = y2
+ fa[5] = x1; fa[10] = y1
+
+ # Draw profiles.
+ j = 0
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ np1 = SFD_NP(sfd)
+ ix = 1 + mod (j, nx)
+ iy = 1 + j / nx
+ j = j + 1
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ call gswind (gp, x1, x2, y1, y2)
+ call gfill (gp, fa, fa[6], 4, GF_SOLID)
+ call gseti (gp, G_DRAWTICKS, NO)
+ call glabax (gp, "", "", "")
+ if (sfd == current) {
+ call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+ vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+ call gsetr (gp, G_PLWIDTH, HLWIDTH)
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ call gpline (gp, fa, fa[6], 5)
+ call gsetr (gp, G_PLWIDTH, 1.)
+ call gseti (gp, G_PLCOLOR, 1)
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ }
+
+ asi = SFD_ASI2(sfd)
+ rmax = stf_i2r (real(np))
+ z = SF_XP1(sf)
+ call gamove (gp, stf_i2r(z)/rmax, asieval(asi,z))
+ for (; z <= SF_XP2(sf); z = z + 0.1)
+ call gadraw (gp, stf_i2r(z)/rmax, asieval(asi,z))
+ if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ np1 = SFD_NP(SF_BEST(sf))
+ asi = SFD_ASI2(SF_BEST(sf))
+ rmax = stf_i2r (real(np))
+ r1 = stf_i2r (SF_XP1(sf))
+ r2 = stf_i2r (SF_XP2(sf))
+ dr = 0.05 * (rmax - stf_i2r (1.))
+ for (r = r1; r <= r2; r = r + dr) {
+ z = stf_r2i (r)
+ z1 = stf_r2i (r+0.7*dr)
+ if (z > 1. && z1 <= np1)
+ call gline (gp, r/rmax, asieval(asi,z),
+ (r+0.7*dr)/rmax, asieval(asi,z1))
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gswind (gp, 0., 1., 0., 1.)
+ call sprintf (Memc[str], SZ_LINE, "%.3g")
+ call pargr (SFD_W(sfd))
+ call gtext (gp, 0.95, 0.98, Memc[str], "h=r;v=t")
+ if (nx < NMAX && ny < NMAX) {
+ call sprintf (Memc[str], SZ_LINE, "%d %d")
+ call pargr (SFD_X(sfd))
+ call pargr (SFD_Y(sfd))
+ call gtext (gp, 0.05, 0.98, Memc[str], "h=l;v=t")
+ }
+ }
+
+ call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+ call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+ call gamove (gp, 1., 1.)
+
+ # Draw label
+ call gseti (gp, G_DRAWAXES, 0)
+ call glabax (gp, title, "", "")
+ call gseti (gp, G_DRAWAXES, 3)
+
+ call sfree (sp)
+end
+
+
+# STF_G6 -- Enclosed flux profile of a star.
+
+procedure stf_g6 (sf, current, title, xlabel, ylabel)
+
+pointer sf #I Starfocus pointer
+pointer current #I Star pointer
+char title[ARB] #I Title
+char xlabel[ARB] #I X label
+char ylabel[ARB] #I Y label
+
+int np, np1
+real scale, level, radius, flux, profile
+pointer gp, asi
+
+real x1, x2, y1, y2, z, z1, r, r1, r2, dr
+real stf_i2r(), stf_r2i(), asieval()
+
+begin
+ gp = SF_GP(sf)
+ level = SF_LEVEL(sf)
+ scale = SF_SCALE(sf)
+ np = SFD_NP(current)
+ asi = SFD_ASI1(current)
+
+ x1 = -0.5 * scale
+ x2 = (stf_i2r (real(np)) + 0.5) * scale
+ y1 = -0.05
+ y2 = 1.05
+ call gswind (gp, x1, x2, y1, y2)
+
+ call gseti (gp, G_DRAWTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call glabax (gp, title, xlabel, ylabel)
+
+ # Draw profiles.
+ if (SFD_STATUS(current) == 0) {
+ call gseti (gp, G_PLCOLOR, 1)
+ for (z = 1.; z <= np; z = z + 1)
+ call gmark (gp, stf_i2r(z)*scale, asieval(asi,z),
+ GM_PLUS, 2., 2.)
+ call gamove (gp, 0., 0.)
+ for (z = 1.; z <= np; z = z + 0.1)
+ call gadraw (gp, stf_i2r(z)*scale, asieval(asi,z))
+ switch (SF_WCODE(sf)) {
+ case 1:
+ radius = SFD_W(current)
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, x1, level, radius, level)
+ call gline (gp, radius, level, radius, y1)
+ call gseti (gp, G_PLTYPE, 1)
+ default:
+ radius = SFD_W(current) / 2.
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, radius, y1, radius, y2)
+ call gseti (gp, G_PLTYPE, 1)
+ }
+
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ call stf_model (sf, current, 0., profile, flux)
+ call gamove (gp, 0., flux)
+ for (z = 1.; z <= np; z = z + 0.1) {
+ r = stf_i2r(z) * scale
+ call stf_model (sf, current, r, profile, flux)
+ call gadraw (gp, r, flux)
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ if (SF_OVRPLT(sf) == YES && current != SF_BEST(sf)) {
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ np1 = SFD_NP(SF_BEST(sf))
+ asi = SFD_ASI1(SF_BEST(sf))
+ r1 = stf_i2r(1.)
+ r2 = stf_i2r (real(np))
+ dr = 0.05 * (r2 - r1)
+ for (r = r1; r <= r2; r = r + dr) {
+ z = stf_r2i (r)
+ z1 = stf_r2i (r+0.7*dr)
+ if (z > 1. && z1 <= np1)
+ call gline (gp, r*scale, asieval(asi,z),
+ (r+0.7*dr)*scale, asieval(asi,z1))
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+ }
+end
+
+
+# STF_G7 -- Radial profile (derivative of enclosed flux) for a star.
+
+procedure stf_g7 (sf, current, title, xlabel, ylabel)
+
+pointer sf #I Starfocus pointer
+pointer current #I Star pointer
+char title[ARB] #I Title
+char xlabel[ARB] #I X label
+char ylabel[ARB] #I Y label
+
+int np, np1
+real scale, level, radius, profile, flux
+pointer gp, asi
+
+real x1, x2, y1, y2, z, z1, r, r1, r2, dr
+real stf_i2r(), stf_r2i(), asieval()
+
+begin
+ gp = SF_GP(sf)
+ level = SF_LEVEL(sf)
+ scale = SF_SCALE(sf)
+ np = SFD_NP(current)
+ asi = SFD_ASI2(current)
+
+ x1 = -0.5 * scale
+ x2 = (stf_i2r (real(np)) + 0.5) * scale
+ z = SFD_YP2(current) - SFD_YP1(current)
+ y1 = SFD_YP1(current) - 0.05 * z
+ y2 = SFD_YP2(current) + 0.05 * z
+ call gswind (gp, x1, x2, y1, y2)
+
+ call gseti (gp, G_XDRAWTICKS, YES)
+ call gseti (gp, G_YDRAWTICKS, NO)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call glabax (gp, title, xlabel, ylabel)
+
+ # Draw profile
+ call gseti (gp, G_PLCOLOR, 1)
+ for (z = SF_XP1(sf); z <= SF_XP2(sf); z = z + 1)
+ call gmark (gp, stf_i2r(z)*scale, asieval(asi,z),
+ GM_PLUS, 2., 2.)
+ z = SF_XP1(sf)
+ call gamove (gp, stf_i2r(z)*scale, asieval(asi,z))
+ for (; z <= SF_XP2(sf); z = z + 0.1)
+ call gadraw (gp, stf_i2r (z)*scale, asieval(asi,z))
+
+ switch (SF_WCODE(sf)) {
+ case 1:
+ radius = SFD_W(current)
+ default:
+ radius = SFD_W(current) / 2.
+ }
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, radius, y1, radius, y2)
+ call gseti (gp, G_PLTYPE, 1)
+
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ z = SF_XP1(sf)
+ r = stf_i2r(z) * scale
+ call stf_model (sf, current, r, profile, flux)
+ call gamove (gp, r, profile)
+ for (; z <= np; z = z + 0.1) {
+ r = stf_i2r(z) * scale
+ call stf_model (sf, current, r, profile, flux)
+ call gadraw (gp, r, profile)
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ if (SF_OVRPLT(sf) == YES && current != SF_BEST(sf)) {
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ np1 = SFD_NP(SF_BEST(sf))
+ asi = SFD_ASI2(SF_BEST(sf))
+ r1 = stf_i2r (SF_XP1(sf))
+ r2 = stf_i2r (SF_XP2(sf))
+ dr = 0.05 * (r2 - r1)
+ for (r = r1; r <= r2; r = r + dr) {
+ z = stf_r2i (r)
+ z1 = stf_r2i (r+0.7*dr)
+ if (z > 1. && z1 <= np1)
+ call gline (gp, r*scale, asieval(asi,z),
+ (r+0.7*dr)*scale, asieval(asi,z1))
+ }
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+end
+
+
+# STF_G8 -- FWHM vs level.
+
+procedure stf_g8 (sf, current, title, xlabel, ylabel)
+
+pointer sf #I Starfocus pointer
+pointer current #I Star pointer
+char title[ARB] #I Title
+char xlabel[ARB] #I X label
+char ylabel[ARB] #I Y label
+
+real y1, y2, level, fwhm
+pointer gp
+
+begin
+ level = SF_LEVEL(sf)
+ if (SF_WCODE(sf) == 1)
+ fwhm = SFD_MFWHM(current)
+ else
+ fwhm = SFD_W(current)
+
+ call alimr (SFD_FWHM(current,2), 17, y1, y2)
+ y2 = y2 - y1
+ y1 = y1 - 0.05 * y2
+ y2 = y1 + 1.10 * y2
+ y1 = min (y1, 0.9 * fwhm)
+ y2 = max (y2, 1.1 * fwhm)
+
+ gp = SF_GP(sf)
+ call gseti (gp, G_DRAWTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call gswind (gp, 0., 1., y1, y2)
+ call glabax (gp, title, xlabel, ylabel)
+ call gvline (gp, SFD_FWHM(current,2), 17, 0.1, 0.9)
+ call gvmark (gp, SFD_FWHM(current,2), 17, 0.1, 0.9, GM_PLUS, 2., 2.)
+
+ switch (SF_WCODE(sf)) {
+ case 1:
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, 0., fwhm, level, fwhm)
+ call gline (gp, level, y1, level, fwhm)
+ call gseti (gp, G_PLTYPE, 1)
+ default:
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, 0., fwhm, 1., fwhm)
+ call gseti (gp, G_PLTYPE, 1)
+ }
+end
+
+
+# STF_G9 -- FWHM vs level for a given star.
+
+procedure stf_g9 (sf, current, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+char title[ARB] #I Title
+
+int i, j, nx, ny, ix, iy
+real level, fwhm, vx, dvx, vy, dvy, x1, x2, y1, y2, fa[10]
+pointer sp, str, gp, sfs, sfd
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ gp = SF_GP(sf)
+ sfs = SFD_SFS(current)
+ level = SF_LEVEL(sf)
+ if (SF_WCODE(sf) == 1)
+ fwhm = SFD_MFWHM(current)
+ else
+ fwhm = SFD_W(current)
+
+ # Set grid layout
+ i = SFS_N(sfs)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+
+ # Set subview port parameters
+ call ggview (gp, vx, dvx, vy, dvy)
+ dvx = (dvx - vx) / nx
+ dvy = (dvy - vy) / ny
+
+ # Set data window parameters
+ y1 = 0.9 * fwhm
+ y2 = 1.1 * fwhm
+ do i = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ call alimr (SFD_FWHM(sfd,2), 17, x1, x2)
+ x2 = x2 - x1
+ x1 = x1 - 0.05 * x2
+ x2 = x1 + 1.10 * x2
+ y1 = min (x1, y1)
+ y2 = max (x2, y2)
+ }
+ x2 = y2 - y1
+ y1 = min (y1, fwhm - 0.2 * x2)
+ y2 = max (y2, fwhm + 0.2 * x2)
+ x1 = 0.
+ x2 = 1.
+ call gswind (gp, x1, x2, y1, y2)
+
+ # Set fill area
+ fa[1] = x1; fa[6] = y1
+ fa[2] = x2; fa[7] = y1
+ fa[3] = x2; fa[8] = y2
+ fa[4] = x1; fa[9] = y2
+ fa[5] = x1; fa[10] = y1
+
+ # Draw profiles.
+ j = 0
+ do i = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ ix = 1 + mod (j, nx)
+ iy = 1 + j / nx
+ j = j + 1
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ call gfill (gp, fa, fa[6], 4, GF_SOLID)
+ call gseti (gp, G_DRAWTICKS, NO)
+ call glabax (gp, "", "", "")
+ if (sfd == current) {
+ call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+ vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+ call gsetr (gp, G_PLWIDTH, HLWIDTH)
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ call gpline (gp, fa, fa[6], 5)
+ call gsetr (gp, G_PLWIDTH, 1.)
+ call gseti (gp, G_PLCOLOR, 1)
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ }
+
+ call gvline (gp, SFD_FWHM(sfd,2), 17, 0.1, 0.9)
+ #call gseti (gp, G_PLTYPE, 2)
+ #call gline (gp, x1, fwhm, x2, fwhm)
+ #call gseti (gp, G_PLTYPE, 1)
+
+ call sprintf (Memc[str], SZ_LINE, "%.3g")
+ call pargr (SFD_W(sfd))
+ call gtext (gp, 0.95, 0.95*y2+0.05*y1, Memc[str], "h=r;v=t")
+ if (nx < NMAX && ny < NMAX) {
+ call sprintf (Memc[str], SZ_LINE, "%.4g")
+ call pargr (SFD_F(sfd))
+ call gtext (gp, 0.05, 0.95*y2+0.05*y1, Memc[str], "h=l;v=t")
+ }
+ }
+
+ call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+ call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+ call gamove (gp, 1., 1.)
+
+ # Draw label
+ call gseti (gp, G_DRAWAXES, 0)
+ call glabax (gp, title, "", "")
+ call gseti (gp, G_DRAWAXES, 3)
+
+ call sfree (sp)
+end
+
+
+# STF_G10 -- FWHM vs level for a given focus.
+
+procedure stf_g10 (sf, current, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+char title[ARB] #I Title
+
+int i, j, nx, ny, ix, iy
+real level, fwhm, vx, dvx, vy, dvy, x1, x2, y1, y2, fa[10]
+pointer sp, str, gp, sff, sfd
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ gp = SF_GP(sf)
+ sff = SFD_SFF(current)
+ level = SF_LEVEL(sf)
+ if (SF_WCODE(sf) == 1)
+ fwhm = SFD_MFWHM(current)
+ else
+ fwhm = SFD_W(current)
+
+ # Set grid layout
+ i = SFF_N(sff)
+ if (i < 4) {
+ nx = i
+ ny = 1
+ } else {
+ nx = nint (sqrt (real (i)))
+ if (mod (i-1, nx+1) >= mod (i-1, nx))
+ nx = nx + 1
+ ny = (i - 1) / nx + 1
+ }
+
+ # Set subview port parameters
+ call ggview (gp, vx, dvx, vy, dvy)
+ dvx = (dvx - vx) / nx
+ dvy = (dvy - vy) / ny
+
+ # Set data window parameters
+ y1 = 0.9 * fwhm
+ y2 = 1.1 * fwhm
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ call alimr (SFD_FWHM(sfd,2), 17, x1, x2)
+ x2 = x2 - x1
+ x1 = x1 - 0.05 * x2
+ x2 = x1 + 1.10 * x2
+ y1 = min (x1, y1)
+ y2 = max (x2, y2)
+ }
+ x2 = y2 - y1
+ y1 = min (y1, fwhm - 0.2 * x2)
+ y2 = max (y2, fwhm + 0.2 * x2)
+ x1 = 0.
+ x2 = 1.
+ call gswind (gp, x1, x2, y1, y2)
+
+ # Set fill area
+ fa[1] = x1; fa[6] = y1
+ fa[2] = x2; fa[7] = y1
+ fa[3] = x2; fa[8] = y2
+ fa[4] = x1; fa[9] = y2
+ fa[5] = x1; fa[10] = y1
+
+ # Draw plots.
+ j = 0
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff, i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ ix = 1 + mod (j, nx)
+ iy = 1 + j / nx
+ j = j + 1
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ call gfill (gp, fa, fa[6], 4, GF_SOLID)
+ call gseti (gp, G_DRAWTICKS, NO)
+ call glabax (gp, "", "", "")
+ if (sfd == current) {
+ call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+ vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+ call gsetr (gp, G_PLWIDTH, HLWIDTH)
+ call gseti (gp, G_PLCOLOR, HLCOLOR)
+ call gpline (gp, fa, fa[6], 5)
+ call gsetr (gp, G_PLWIDTH, 1.)
+ call gseti (gp, G_PLCOLOR, 1)
+ call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+ vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+ }
+
+ call gvline (gp, SFD_FWHM(sfd,2), 17, 0.1, 0.9)
+ #call gseti (gp, G_PLTYPE, 2)
+ #call gline (gp, x1, fwhm, x2, fwhm)
+ #call gseti (gp, G_PLTYPE, 1)
+
+ call sprintf (Memc[str], SZ_LINE, "%.3g")
+ call pargr (SFD_W(sfd))
+ call gtext (gp, 0.95, 0.95*y2+0.05*y1, Memc[str], "h=r;v=t")
+ if (nx < NMAX && ny < NMAX) {
+ call sprintf (Memc[str], SZ_LINE, "%d %d")
+ call pargr (SFD_X(sfd))
+ call pargr (SFD_Y(sfd))
+ call gtext (gp, 0.05, 0.95*y2+0.05*y1, Memc[str], "h=l;v=t")
+ }
+ }
+
+ call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+ call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+ call gamove (gp, 1., 1.)
+
+ # Draw label
+ call gseti (gp, G_DRAWAXES, 0)
+ call glabax (gp, title, "", "")
+ call gseti (gp, G_DRAWAXES, 3)
+
+ call sfree (sp)
+end
+
+
+# STF_G11 -- Spatial plot at one focus.
+
+procedure stf_g11 (sf, current, key, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+int key #I Plot magnitude symbol?
+char title[ARB] #I Title
+
+int i
+real x, y, z, x1, x2, y1, y2, rbest, rmin, rmax, emin, emax
+real vx[3,2], vy[3,2], dvx, dvy, fa[8]
+pointer gp, sfd, sff
+
+data fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+ gp = SF_GP(sf)
+ sff = SFD_SFF(current)
+
+ # Range of X, Y, R, E.
+ x1 = 1.
+ y1 = 1.
+ x2 = SF_NCOLS(sf)
+ y2 = SF_NLINES(sf)
+
+ rbest = SFD_W(SF_BEST(sf))
+ rmin = SF_W(sf)
+ rmax = 1.5 * SF_W(sf)
+ emin = 0
+ emax = 1
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ rmin = min (rmin, SFD_W(sfd))
+ rmax = max (rmax, SFD_W(sfd))
+ emin = min (emin, SFD_E(sfd))
+ emax = max (emax, SFD_E(sfd))
+ }
+ z = rmax - rmin
+ rmin = rmin - 0.1 * z
+ rmax = rmax + 0.1 * z
+
+ # Set view ports
+ call ggview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+ dvx = vx[3,2] - vx[1,1]
+ dvy = vy[3,2] - vy[1,1]
+ vx[1,1] = vx[1,1] + 0.00 * dvx
+ vx[1,2] = vx[1,1] + 0.20 * dvx
+ vx[2,1] = vx[1,1] + 0.25 * dvx
+ vx[2,2] = vx[1,1] + 0.75 * dvx
+ vx[3,1] = vx[1,1] + 0.80 * dvx
+ vx[3,2] = vx[1,1] + 1.00 * dvx
+ vy[1,1] = vy[1,1] + 0.00 * dvy
+ vy[1,2] = vy[1,1] + 0.20 * dvy
+ vy[2,1] = vy[1,1] + 0.25 * dvy
+ vy[2,2] = vy[1,1] + 0.75 * dvy
+ vy[3,1] = vy[1,1] + 0.80 * dvy
+ vy[3,2] = vy[1,1] + 1.00 * dvy
+
+ # (X,R)
+ call gseti (gp, G_WCS, 2)
+ call gseti (gp, G_DRAWAXES, 3)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gseti (gp, G_YLABELTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 4)
+ call gseti (gp, G_YNMINOR, 0)
+ call gsview (gp, vx[2,1], vx[2,2], vy[1,1], vy[1,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, x1, x2, rmin, rmax)
+ call glabax (gp, "", "Column", "")
+
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ x = SFD_X(sfd)
+ y = SFD_W(sfd)
+ if (key == 1) {
+ z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFD_W(sfd) < SF_W(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFD_W(sfd) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, x1, SF_W(sf), x2, SF_W(sf))
+ call gseti (gp, G_PLTYPE, 1)
+
+ # (R,Y)
+ call gseti (gp, G_WCS, 3)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gseti (gp, G_YLABELTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 4)
+ call gseti (gp, G_XNMINOR, 0)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call gsview (gp, vx[1,1], vx[1,2], vy[2,1], vy[2,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, rmin, rmax, y1, y2)
+ call glabax (gp, "", SF_WTYPE(sf), "Line")
+
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ x = SFD_W(sfd)
+ y = SFD_Y(sfd)
+ if (key == 1) {
+ z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFD_W(sfd) < SF_W(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFD_W(sfd) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, SF_W(sf), y1, SF_W(sf), y2)
+ call gseti (gp, G_PLTYPE, 1)
+
+ # (E,R)
+ call gseti (gp, G_WCS, 4)
+ call gseti (gp, G_DRAWAXES, 3)
+ call gseti (gp, G_XLABELTICKS, NO)
+ call gseti (gp, G_YLABELTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 4)
+ call gseti (gp, G_YNMINOR, 0)
+ call gsview (gp, vx[2,1], vx[2,2], vy[3,1], vy[3,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, x1, x2, emin, emax)
+ call glabax (gp, "", "", "Ellip")
+
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ x = SFD_X(sfd)
+ y = SFD_E(sfd)
+ if (key == 1) {
+ z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFD_W(sfd) < SF_W(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFD_W(sfd) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ # (E,Y)
+ call gseti (gp, G_WCS, 5)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gseti (gp, G_YLABELTICKS, NO)
+ call gseti (gp, G_XNMAJOR, 4)
+ call gseti (gp, G_XNMINOR, 0)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call gsview (gp, vx[3,1], vx[3,2], vy[2,1], vy[2,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, emin, emax, y1, y2)
+ call glabax (gp, "", "Ellip", "")
+
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ x = SFD_E(sfd)
+ y = SFD_Y(sfd)
+ if (key == 1) {
+ z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFD_W(sfd) < SF_W(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFD_W(sfd) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ # Label window.
+ call gseti (gp, G_WCS, 1)
+ call gseti (gp, G_DRAWAXES, 0)
+ call gsview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+ call glabax (gp, title, "", "")
+
+ # (X,Y)
+ call gseti (gp, G_DRAWAXES, 3)
+ call gseti (gp, G_LABELTICKS, NO)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call gsview (gp, vx[2,1], vx[2,2], vy[2,1], vy[2,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, x1, x2, y1, y2)
+ call glabax (gp, "", "", "")
+
+ do i = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,i)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ x = SFD_X(sfd)
+ y = SFD_Y(sfd)
+ if (key == 1) {
+ z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFD_W(sfd) < SF_W(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFD_W(sfd) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+end
+
+
+# STF_G9 -- Spatial plots at best focus.
+
+procedure stf_g12 (sf, current, key, title)
+
+pointer sf #I Starfocus pointer
+pointer current #I Current sfd pointer
+int key #I Plot magnitude symbol?
+char title[ARB] #I Title
+
+int i
+real x, y, z, x1, x2, y1, y2, fmin, fmax, rbest, rmin, rmax
+real vx[3,2], vy[3,2], dvx, dvy, fa[8]
+pointer gp, sfs, sfd
+
+data fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+ gp = SF_GP(sf)
+
+ # Range of X, Y, R, F.
+ x1 = 1.
+ y1 = 1.
+ x2 = SF_NCOLS(sf)
+ y2 = SF_NLINES(sf)
+
+ rbest = SFD_W(SF_BEST(sf))
+ fmin = MAX_REAL
+ fmax = -MAX_REAL
+ rmin = SF_W(sf)
+ rmax = 1.5 * SF_W(sf)
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_N(sfs) == 0)
+ next
+ fmin = min (fmin, SFS_F(sfs))
+ fmax = max (fmax, SFS_F(sfs))
+ rmin = min (rmin, SFS_W(sfs))
+ rmax = max (rmax, SFS_W(sfs))
+ }
+ z = fmax - fmin
+ fmin = fmin - 0.1 * z
+ fmax = fmax + 0.1 * z
+ z = rmax - rmin
+ rmin = rmin - 0.1 * z
+ rmax = rmax + 0.1 * z
+
+ # Set view ports
+ call ggview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+ dvx = vx[3,2] - vx[1,1]
+ dvy = vy[3,2] - vy[1,1]
+ vx[1,1] = vx[1,1] + 0.00 * dvx
+ vx[1,2] = vx[1,1] + 0.20 * dvx
+ vx[2,1] = vx[1,1] + 0.25 * dvx
+ if (SF_NF(sf) > 1) {
+ vx[2,2] = vx[1,1] + 0.75 * dvx
+ vx[3,1] = vx[1,1] + 0.80 * dvx
+ vx[3,2] = vx[1,1] + 1.00 * dvx
+ } else {
+ vx[2,2] = vx[1,1] + 1.00 * dvx
+ vx[3,1] = vx[1,1] + 1.00 * dvx
+ vx[3,2] = vx[1,1] + 1.00 * dvx
+ }
+ vy[1,1] = vy[1,1] + 0.00 * dvy
+ vy[1,2] = vy[1,1] + 0.20 * dvy
+ vy[2,1] = vy[1,1] + 0.25 * dvy
+ if (SF_NF(sf) > 1) {
+ vy[2,2] = vy[1,1] + 0.75 * dvy
+ vy[3,1] = vy[1,1] + 0.80 * dvy
+ vy[3,2] = vy[1,1] + 1.00 * dvy
+ } else {
+ vy[2,2] = vy[1,1] + 1.00 * dvy
+ vy[3,1] = vy[1,1] + 1.00 * dvy
+ vy[3,2] = vy[1,1] + 1.00 * dvy
+ }
+
+ dvx = vx[2,1] - vx[2,2]
+ dvy = vy[1,1] - vy[1,2]
+ if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+ # (X,R)
+ call gseti (gp, G_WCS, 2)
+ call gseti (gp, G_DRAWAXES, 3)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gseti (gp, G_YLABELTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 4)
+ call gseti (gp, G_YNMINOR, 0)
+ call gsview (gp, vx[2,1], vx[2,2], vy[1,1], vy[1,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, x1, x2, rmin, rmax)
+ call glabax (gp, "", "Column", "")
+
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_N(sfs) == 0)
+ next
+ x = SFD_X(SFS_SFD(sfs,1))
+ y = SFS_W(sfs)
+ if (key == 1) {
+ z = sqrt (SFS_M(sfs) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFS_F(sfs) < SF_F(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFS_W(sfs) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, x1, SF_W(sf), x2, SF_W(sf))
+ call gseti (gp, G_PLTYPE, 1)
+ }
+
+ dvx = vx[1,1] - vx[1,2]
+ dvy = vy[2,1] - vy[2,2]
+ if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+ # (R,Y)
+ call gseti (gp, G_WCS, 3)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gseti (gp, G_YLABELTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 4)
+ call gseti (gp, G_XNMINOR, 0)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call gsview (gp, vx[1,1], vx[1,2], vy[2,1], vy[2,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, rmin, rmax, y1, y2)
+ call glabax (gp, "", SF_WTYPE(sf), "Line")
+
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_N(sfs) == 0)
+ next
+ x = SFS_W(sfs)
+ y = SFD_Y(SFS_SFD(sfs,1))
+ if (key == 1) {
+ z = sqrt (SFS_M(sfs) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFS_F(sfs) < SF_F(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFS_W(sfs) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, SF_W(sf), y1, SF_W(sf), y2)
+ call gseti (gp, G_PLTYPE, 1)
+ }
+
+ dvx = vx[2,1] - vx[2,2]
+ dvy = vy[3,1] - vy[3,2]
+ if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+ # (X,F)
+ call gseti (gp, G_WCS, 4)
+ call gseti (gp, G_XLABELTICKS, NO)
+ call gseti (gp, G_YLABELTICKS, YES)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 4)
+ call gseti (gp, G_YNMINOR, 0)
+ call gsview (gp, vx[2,1], vx[2,2], vy[3,1], vy[3,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, x1, x2, fmin, fmax)
+ call glabax (gp, "", "", "Focus")
+
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_N(sfs) == 0)
+ next
+ x = SFD_X(SFS_SFD(sfs,1))
+ y = SFS_F(sfs)
+ if (key == 1) {
+ z = sqrt (SFS_M(sfs) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFS_F(sfs) < SF_F(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFS_W(sfs) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, x1, SF_F(sf), x2, SF_F(sf))
+ call gseti (gp, G_PLTYPE, 1)
+ }
+
+ dvx = vx[3,1] - vx[3,2]
+ dvy = vy[2,1] - vy[2,2]
+ if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+ # (F,Y)
+ call gseti (gp, G_WCS, 5)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gseti (gp, G_YLABELTICKS, NO)
+ call gseti (gp, G_XNMAJOR, 4)
+ call gseti (gp, G_XNMINOR, 0)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call gsview (gp, vx[3,1], vx[3,2], vy[2,1], vy[2,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, fmin, fmax, y1, y2)
+ call glabax (gp, "", "Focus", "")
+
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_N(sfs) == 0)
+ next
+ x = SFS_F(sfs)
+ y = SFD_Y(SFS_SFD(sfs,1))
+ if (key == 1) {
+ z = sqrt (SFS_M(sfs) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFS_F(sfs) < SF_F(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFS_W(sfs) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, SF_F(sf), y1, SF_F(sf), y2)
+ call gseti (gp, G_PLTYPE, 1)
+ }
+
+ # Label window.
+ call gseti (gp, G_WCS, 1)
+ call gseti (gp, G_DRAWAXES, 0)
+ call gsview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+ call glabax (gp, title, "", "")
+
+ dvx = vx[2,1] - vx[2,2]
+ dvy = vy[2,1] - vy[2,2]
+ if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+ # (X,Y)
+ call gseti (gp, G_DRAWAXES, 3)
+ call gseti (gp, G_LABELTICKS, NO)
+ call gseti (gp, G_XNMAJOR, 6)
+ call gseti (gp, G_XNMINOR, 4)
+ call gseti (gp, G_YNMAJOR, 6)
+ call gseti (gp, G_YNMINOR, 4)
+ call gsview (gp, vx[2,1], vx[2,2], vy[2,1], vy[2,2])
+ call gswind (gp, 0., 1., 0., 1.)
+ call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+ call gswind (gp, x1, x2, y1, y2)
+ call glabax (gp, "", "", "")
+
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_N(sfs) == 0)
+ next
+ sfd = SFS_SFD(sfs,1)
+ x = SFD_X(sfd)
+ y = SFD_Y(sfd)
+ if (key == 1) {
+ z = sqrt (SFS_M(sfs) / SF_M(sf))
+ z = max (0.005, 0.03 * z)
+ call gmark (gp, x, y, GM_MAG, z, z)
+ }
+ if (SFS_F(sfs) < SF_F(sf))
+ call gseti (gp, G_PLCOLOR, 2)
+ else
+ call gseti (gp, G_PLCOLOR, 3)
+ z = min (2., SFS_W(sfs) / rbest)
+ z = 0.010 * (1 + (z - 1) * 5)
+ call gmark (gp, x, y, GM_CIRCLE, z, z)
+ call gseti (gp, G_PLCOLOR, 1)
+ }
+ }
+end
diff --git a/noao/obsutil/src/starfocus/stfhelp.key b/noao/obsutil/src/starfocus/stfhelp.key
new file mode 100644
index 00000000..948bb9b7
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfhelp.key
@@ -0,0 +1,63 @@
+ STARFOCUS COMMAND OPTIONS
+
+ SUMMARY
+
+? Help f Focus p Radial profile u Undelete
+a Spatial i Info q Quit x Delete
+b Best m Magnitude r Redraw z Zoom
+d Delete n Normalize s Mag symbols <space> Next
+e Enclosed flux o Offset t Field radius
+
+:level <val> :radius <val> :show <file> :xcenter <val>
+:overplot <y|n> :scale <val> :size <type> :ycenter <val>
+
+
+ CURSOR COMMANDS
+
+All plots may not be available depending on the number of focus values and
+the number of stars.
+
+? Page this help information
+a Spatial plot at a single focus
+b Spatial plot of best focus values
+d Delete star nearest to cursor
+e Enclosed flux for all stars at one focus and all focus for one star
+f Size and ellipticity vs focus for all data
+i Information about point nearest the cursor
+m Size and ellipticity vs relative magnitude at one focus
+n Normalize enclosed flux at x cursor position
+o Offset enclosed flux to x,y cursor position by adjusting background
+p Radial profiles for all stars at one focus and all focus for one star
+ The profiles are determined from the derivatives of the enclosed flux
+q Quit
+r Redraw
+s Toggle magnitude symbols in spatial plots
+t Size and ellipticity vs radius from field center at one focus
+u Undelete all deleted points
+x Delete nearest point, star, or focus (selected by query)
+z Zoom to a single measurement showing enclosed flux and radial profile
+<space> Step through different focus or stars in current plot type
+
+
+ COLON COMMANDS
+
+A command without a value generally shows the current value of the
+parameter while with a value it sets the value of the parameter.
+
+:level <val> Level at which the size parameter is evaluated
+:overplot <y|n> Overplot the profiles from the narrowest profile?
+:radius <val> Change profile radius(*)
+:show <file> Page all information for the current set of objects
+:size <type> Size type (Radius|FWHM|GFWHM|MFWHM) (**)
+:scale <val> Pixel scale for size values
+:xcenter <val> X field center for radius from field center plots
+:ycenter <val> Y field center for radius from field center plots
+
+(*) The profile radius may not exceed the initial value set by the task
+ parameter.
+
+(**)
+Radius = radius enclosing the fraction of the flux specified by "level"
+FWHM = Full-width at half-maximum based on the radially smoothed profile
+GFWHM = Full-width at half-maximum of Gaussian function fit to enclosed flux
+MFWHM = Full-width at half-maximum of Moffat function fit to enclosed flux
diff --git a/noao/obsutil/src/starfocus/stfmeasure.x b/noao/obsutil/src/starfocus/stfmeasure.x
new file mode 100644
index 00000000..808b6058
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfmeasure.x
@@ -0,0 +1,134 @@
+include <error.h>
+include <imhdr.h>
+include <imset.h>
+include <math/nlfit.h>
+include "starfocus.h"
+
+
+# STF_MEASURE -- PSF measuring routine.
+# This is a stand-alone routine that can be called to return the FWHM.
+# It is a greatly abbreviated version of starfocus.
+
+procedure stf_measure (im, xc, yc, beta, level, radius, nit,
+ sbuffer, swidth, saturation, gp, logfd,
+ bkg, renclosed, dfwhm, gfwhm, mfwhm)
+
+pointer im #I Image pointer
+real xc #I Initial X center
+real yc #I Initial Y center
+real beta #I Moffat beta
+real level #I Measurement level
+real radius #U Profile radius
+int nit #I Number of iterations on radius
+real sbuffer #I Sky buffer (pixels)
+real swidth #I Sky width (pixels)
+real saturation #I Saturation
+pointer gp #I Graphics output if not NULL
+int logfd #I Log output if not NULL
+real bkg #O Background used
+real renclosed #O Enclosed flux radius
+real dfwhm #O Direct FWHM
+real gfwhm #O Gaussian FWHM
+real mfwhm #O Moffat FWHM
+
+int i
+bool ignore_sat
+pointer sp, str, sf, sfd, sfds
+
+int strdic()
+real stf_r2i()
+errchk stf_find, stf_bkgd, stf_profile, stf_widths, stf_fwhms, stf_organize
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call salloc (sf, SF, TY_STRUCT)
+ call salloc (sfd, SFD, TY_STRUCT)
+ call salloc (sfds, 1, TY_POINTER)
+ call aclri (Memi[sf], SF)
+ call aclri (Memi[sfd], SFD)
+ Memi[sfds] = sfd
+
+ # Initialize parameters.
+ SF_TASK(sf) = PSFMEASURE
+ SF_WCODE(sf) = strdic ("FWHM", SF_WTYPE(sf), SF_SZWTYPE, SF_WTYPES)
+ SF_SCALE(sf) = 1.
+ SF_LEVEL(sf) = level
+ SF_BETA(sf) = beta
+ SF_RADIUS(sf) = radius
+ SF_SBUF(sf) = sbuffer
+ SF_SWIDTH(sf) = swidth
+ SF_SAT(sf) = saturation
+ SF_NIT(sf) = nit
+ SF_OVRPLT(sf) = NO
+ SF_NCOLS(sf) = IM_LEN(im,1)
+ SF_NLINES(sf) = IM_LEN(im,2)
+ SF_XF(sf) = (IM_LEN(im,1) + 1) / 2.
+ SF_YF(sf) = (IM_LEN(im,2) + 1) / 2.
+ ignore_sat = false
+
+ call imstats (im, IM_IMAGENAME, SFD_IMAGE(sfd), SF_SZFNAME)
+ SFD_ID(sfd) = 1
+ SFD_X(sfd) = xc
+ SFD_Y(sfd) = yc
+ SFD_F(sfd) = INDEF
+ SFD_STATUS(sfd) = 0
+ SFD_SFS(sfd) = NULL
+ SFD_SFF(sfd) = NULL
+ SFD_SFI(sfd) = NULL
+
+ if (SF_LEVEL(sf) > 1.)
+ SF_LEVEL(sf) = SF_LEVEL(sf) / 100.
+ SF_LEVEL(sf) = max (0.05, min (0.95, SF_LEVEL(sf)))
+
+ # Evaluate PSF data.
+ iferr {
+ do i = 1, SF_NIT(sf) {
+ if (i == 1)
+ SFD_RADIUS(sfd) = SF_RADIUS(sf)
+ else
+ SFD_RADIUS(sfd) = 3. * SFD_DFWHM(sfd)
+ SFD_NPMAX(sfd) = stf_r2i (SFD_RADIUS(sfd)) + 1
+ SFD_NP(sfd) = SFD_NPMAX(sfd)
+ call stf_find (sf, sfd, im)
+ call stf_bkgd (sf, sfd)
+ if (SFD_NSAT(sfd) > 0 && i == 1) {
+ if (ignore_sat)
+ call error (0,
+ "Saturated pixels found - ignoring object")
+ else
+ call eprintf (
+ "WARNING: Saturated pixels found.\n")
+ }
+ call stf_profile (sf, sfd)
+ call stf_widths (sf, sfd)
+ call stf_fwhms (sf, sfd)
+ }
+
+ # Set output results.
+ radius = SFD_RADIUS(sfd)
+ bkg = SFD_BKGD(sfd)
+ renclosed = SFD_R(sfd)
+ dfwhm = SFD_DFWHM(sfd)
+ mfwhm = SFD_MFWHM(sfd)
+ gfwhm = SFD_GFWHM(sfd)
+
+ # Optional graph and log output. Note that the gp pointer is only
+ # used to indicate whether to make a graph. The stf_graph
+ # procedure opens its own graphics stream.
+
+ call stf_organize (sf, sfds, 1)
+ if (gp != NULL)
+ call stf_graph (sf)
+ if (logfd != NULL)
+ call stf_log (sf, logfd)
+
+ call asifree (SFD_ASI1(sfd))
+ call asifree (SFD_ASI2(sfd))
+ } then
+ call erract (EA_WARN)
+
+ # Finish up
+ call stf_free (sf)
+ call sfree (sp)
+end
diff --git a/noao/obsutil/src/starfocus/stfprofile.x b/noao/obsutil/src/starfocus/stfprofile.x
new file mode 100644
index 00000000..d26c085d
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfprofile.x
@@ -0,0 +1,1189 @@
+include <imhdr.h>
+include <mach.h>
+include <math.h>
+include <math/iminterp.h>
+include <math/nlfit.h>
+include "starfocus.h"
+
+
+# STF_FIND -- Find the object and return the data raster and object center.
+# STF_BKGD -- Compute the background.
+# STF_PROFILE -- Compute enclosed flux profile, derivative, and moments.
+# STF_NORM -- Renormalized enclosed flux profile
+# STF_WIDTHS -- Set widths.
+# STF_I2R -- Radius from sample index.
+# STF_R2I -- Sample index from radius.
+# STF_R2N -- Number of subsamples from radius.
+# STF_MODEL -- Return model values.
+# STF_DFWHM -- Direct FWHM from profile.
+# STF_FWHMS -- Measure FWHM vs level.
+# STF_RADIUS -- Measure the radius at the specified level.
+# STF_FIT -- Fit model.
+# STF_GAUSS1 -- Gaussian function used in NLFIT.
+# STF_GAUSS2 -- Gaussian function and derivatives used in NLFIT.
+# STF_MOFFAT1 -- Moffat function used in NLFIT.
+# STF_MOFFAT2 -- Moffat function and derivatives used in NLFIT.
+
+
+# STF_FIND -- Find the object and return the data raster and object center.
+# Centering uses centroid of marginal distributions of data above the mean.
+
+procedure stf_find (sf, sfd, im)
+
+pointer sf #I Starfocus pointer
+pointer sfd #I Object pointer
+pointer im #I Image pointer
+
+long lseed
+int i, j, k, x1, x2, y1, y2, nx, ny, npts
+real radius, buffer, width, xc, yc, xlast, ylast, r1, r2
+real mean, sum, sum1, sum2, sum3, asumr(), urand()
+pointer data, ptr, imgs2r()
+errchk imgs2r
+
+begin
+ radius = max (3., SFD_RADIUS(sfd))
+ buffer = SF_SBUF(sf)
+ width = SF_SWIDTH(sf)
+
+ xc = SFD_X(sfd)
+ yc = SFD_Y(sfd)
+ r1 = radius + buffer + width
+ r2 = radius
+
+ # Iterate on the center finding.
+ do k = 1, 3 {
+
+ # Extract region around current center.
+ xlast = xc
+ ylast = yc
+
+ x1 = max (1-NBNDRYPIX, nint (xc - r2))
+ x2 = min (IM_LEN(im,1)+NBNDRYPIX, nint (xc + r2))
+ nx = x2 - x1 + 1
+ y1 = max (1-NBNDRYPIX, nint (yc - r2))
+ y2 = min (IM_LEN(im,2)+NBNDRYPIX, nint (yc + r2))
+ ny = y2 - y1 + 1
+ npts = nx * ny
+ data = imgs2r (im, x1, x2, y1, y2)
+
+ # Find center of gravity of marginal distributions above mean.
+ npts = nx * ny
+ sum = asumr (Memr[data], npts)
+ mean = sum / nx
+ sum1 = 0.
+ sum2 = 0.
+
+ do i = x1, x2 {
+ ptr = data + i - x1
+ sum3 = 0.
+ do j = y1, y2 {
+ sum3 = sum3 + Memr[ptr]
+ ptr = ptr + nx
+ }
+ sum3 = sum3 - mean
+ if (sum3 > 0.) {
+ sum1 = sum1 + i * sum3
+ sum2 = sum2 + sum3
+ }
+ }
+ if (sum2 <= 0)
+ call error (1, "Centering failed to converge")
+ xc = sum1 / sum2
+ if (xlast - xc > 0.2 * nx)
+ xc = xlast - 0.2 * nx
+ if (xc - xlast > 0.2 * nx)
+ xc = xlast + 0.2 * nx
+
+ ptr = data
+ mean = sum / ny
+ sum1 = 0.
+ sum2 = 0.
+ do j = y1, y2 {
+ sum3 = 0.
+ do i = x1, x2 {
+ sum3 = sum3 + Memr[ptr]
+ ptr = ptr + 1
+ }
+ sum3 = sum3 - mean
+ if (sum3 > 0.) {
+ sum1 = sum1 + j * sum3
+ sum2 = sum2 + sum3
+ }
+ }
+ if (sum2 <= 0)
+ call error (1, "Centering failed to converge")
+ yc = sum1 / sum2
+ if (ylast - yc > 0.2 * ny)
+ yc = ylast - 0.2 * ny
+ if (yc - ylast > 0.2 * ny)
+ yc = ylast + 0.2 * ny
+
+ if (nint(xc) == nint(xlast) && nint(yc) == nint(ylast))
+ break
+ }
+
+ # Get a new centered raster if necessary.
+ if (nint(xc) != nint(xlast) || nint(yc) != nint(ylast) || r2 < r1) {
+ x1 = max (1-NBNDRYPIX, nint (xc - r1))
+ x2 = min (IM_LEN(im,1)+NBNDRYPIX, nint (xc + r1))
+ nx = x2 - x1 + 1
+ y1 = max (1-NBNDRYPIX, nint (yc - r1))
+ y2 = min (IM_LEN(im,2)+NBNDRYPIX, nint (yc + r1))
+ ny = y2 - y1 + 1
+ npts = nx * ny
+ data = imgs2r (im, x1, x2, y1, y2)
+ }
+
+ # Add a dither for integer data. The random numbers are always
+ # the same to provide reproducibility.
+
+ i = IM_PIXTYPE(im)
+ if (i == TY_SHORT || i == TY_INT || i == TY_LONG) {
+ lseed = 1
+ do i = 0, npts-1
+ Memr[data+i] = Memr[data+i] + urand(lseed) - 0.5
+ }
+
+ SFD_DATA(sfd) = data
+ SFD_X1(sfd) = x1
+ SFD_X2(sfd) = x2
+ SFD_Y1(sfd) = y1
+ SFD_Y2(sfd) = y2
+ SFD_X(sfd) = xc
+ SFD_Y(sfd) = yc
+end
+
+
+# STF_BKGD -- Compute the background.
+# A mode is estimated from the minimum slope in the sorted background pixels
+# with a bin width of 5%.
+
+procedure stf_bkgd (sf, sfd)
+
+pointer sf #I Parameter structure
+pointer sfd #I Star structure
+
+int i, j, x1, x2, y1, y2, xc, yc, nx, ny, npts, ns, nsat
+real sat, bkgd, miso
+real r, r1, r2, r3, dx, dy, dz
+pointer sp, data, bdata, ptr
+
+begin
+ data = SFD_DATA(sfd)
+ x1 = SFD_X1(sfd)
+ x2 = SFD_X2(sfd)
+ y1 = SFD_Y1(sfd)
+ y2 = SFD_Y2(sfd)
+ xc = SFD_X(sfd)
+ yc = SFD_Y(sfd)
+
+ nx = x2 - x1 + 1
+ ny = y2 - y1 + 1
+ npts = nx * ny
+
+ ns = 0
+ nsat = 0
+ r1 = SFD_RADIUS(sfd) ** 2
+ r2 = (SFD_RADIUS(sfd) + SF_SBUF(sf)) ** 2
+ r3 = (SFD_RADIUS(sfd) + SF_SBUF(sf) + SF_SWIDTH(sf)) ** 2
+ sat = SF_SAT(sf)
+ if (IS_INDEF(sat))
+ sat = MAX_REAL
+
+ call smark (sp)
+ call salloc (bdata, npts, TY_REAL)
+
+ ptr = data
+ do j = y1, y2 {
+ dy = (yc - j) ** 2
+ do i = x1, x2 {
+ dx = (xc - i) ** 2
+ r = dx + dy
+ if (r <= r1) {
+ if (Memr[ptr] >= sat)
+ nsat = nsat + 1
+ } else if (r >= r2 && r <= r3) {
+ Memr[bdata+ns] = Memr[ptr]
+ ns = ns + 1
+ }
+ ptr = ptr + 1
+ }
+ }
+
+ if (ns > 9) {
+ call asrtr (Memr[bdata], Memr[bdata], ns)
+ r = Memr[bdata+ns-1] - Memr[bdata]
+ bkgd = Memr[bdata] + r / 2
+ miso = r / 2
+
+ j = 1 + 0.50 * ns
+ do i = 0, ns - j {
+ dz = Memr[bdata+i+j-1] - Memr[bdata+i]
+ if (dz < r) {
+ r = dz
+ bkgd = Memr[bdata+i] + dz / 2
+ miso = dz / 2
+ }
+ }
+ } else {
+ bkgd = 0.
+ miso = 0.
+ }
+
+ SFD_BKGD1(sfd) = bkgd
+ SFD_BKGD(sfd) = bkgd
+ SFD_MISO(sfd) = miso
+ SFD_NSAT(sfd) = nsat
+
+ call sfree (sp)
+end
+
+
+# STF_PROFILE -- Compute enclosed flux profile, derivative, direct FWHM, and
+# profile moments..
+# 1. The flux profile is normalized at the maximum value.
+# 2. The radial profile is computed from the numerical derivative of the
+# enclose flux profile.
+
+procedure stf_profile (sf, sfd)
+
+pointer sf #I Parameter structure
+pointer sfd #I Star structure
+
+int np
+real radius, xc, yc
+
+int i, j, k, l, m, ns, nx, ny, x1, x2, y1, y2
+real bkgd, miso, sigma, peak
+real r, r1, r2, r3, dx, dy, dx1, dx2, dy1, dy2, dz, xx, yy, xy, ds, da
+pointer sp, data, profile, ptr, asi, msi, gs
+int stf_r2n()
+real asieval(), msieval(), gseval(), stf_i2r(), stf_r2i()
+errchk asiinit, asifit, msiinit, msifit, gsrestore
+
+real gsdata[24]
+data gsdata/ 1., 4., 4., 1., 0., 0.6726812, 1., 2., 1.630641, 0.088787,
+ 0.00389378, -0.001457133, 0.3932125, -0.1267456, -0.004864541,
+ 0.00249941, 0.03078612, 0.02731274, -4.875850E-4, 2.307464E-4,
+ -0.002134843, 0.007603908, -0.002552385, -8.010564E-4/
+
+begin
+ data = SFD_DATA(sfd)
+ x1 = SFD_X1(sfd)
+ x2 = SFD_X2(sfd)
+ y1 = SFD_Y1(sfd)
+ y2 = SFD_Y2(sfd)
+ xc = SFD_X(sfd)
+ yc = SFD_Y(sfd)
+ bkgd = SFD_BKGD(sfd)
+ miso = SFD_MISO(sfd)
+ radius = SFD_RADIUS(sfd)
+ np = SFD_NP(sfd)
+
+ nx = x2 - x1 + 1
+ ny = y2 - y1 + 1
+
+ # Use an image interpolator fit to the data.
+ call msiinit (msi, II_BISPLINE3)
+ call msifit (msi, Memr[data], nx, ny, nx)
+
+ # To avoid trying to interpolate outside the center of the
+ # edge pixels, a requirement of the interpolator functions,
+ # we reset the data limits.
+ x1 = x1 + 1
+ x2 = x2 - 1
+ y1 = y1 + 1
+ y2 = y2 - 1
+
+ # Compute the enclosed flux profile, its derivative, and moments.
+ call smark (sp)
+ call salloc (profile, np, TY_REAL)
+ call aclrr (Memr[profile], np)
+
+ xx = 0.
+ yy = 0.
+ xy = 0.
+ do j = y1, y2 {
+ ptr = data + (j-y1+1)*nx + 1
+ dy = j - yc
+ do i = x1, x2 {
+ dx = i - xc
+
+ # Set the subpixel sampling which may be a function of radius.
+ r = sqrt (dx * dx + dy * dy)
+ ns = stf_r2n (r)
+ ds = 1. / ns
+ da = ds * ds
+ dz = 0.5 + 0.5 * ds
+
+ # Sum the interpolator values over the subpixels and compute
+ # an offset to give the correct total for the pixel.
+
+ r2 = 0.
+ dy1 = dy - dz
+ do l = 1, ns {
+ dy1 = dy1 + ds
+ dy2 = dy1 * dy1
+ dx1 = dx - dz
+ do k = 1, ns {
+ dx1 = dx1 + ds
+ dx2 = dx1 * dx1
+ r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2)
+ r2 = r2 + r1
+ }
+ }
+
+ r1 = Memr[ptr] - bkgd
+ ptr = ptr + 1
+ r2 = r1 - r2 * da
+
+ # Accumulate the enclosed flux over the sub pixels.
+ dy1 = dy - dz
+ do l = 1, ns {
+ dy1 = dy1 + ds
+ dy2 = dy1 * dy1
+ dx1 = dx - dz
+ do k = 1, ns {
+ dx1 = dx1 + ds
+ dx2 = dx1 * dx1
+ r = max (0., sqrt (dx2 + dy2) - ds / 2)
+ if (r < radius) {
+ r1 = da * (msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2) +
+ r2)
+
+ # Use approximation for fractions of a subpixel.
+ for (m=stf_r2i(r)+1; m<=np; m=m+1) {
+ r3 = (stf_i2r (real(m)) - r) / ds
+ if (r3 >= 1.)
+ break
+ Memr[profile+m-1] = Memr[profile+m-1] + r3 * r1
+ }
+
+ # The subpixel is completely within these radii.
+ for (; m<=np; m=m+1)
+ Memr[profile+m-1] = Memr[profile+m-1] + r1
+
+ # Accumulate the moments above an isophote.
+ if (r1 > miso) {
+ xx = xx + dx2 * r1
+ yy = yy + dy2 * r1
+ xy = xy + dx1 * dy1 * r1
+ }
+ }
+ }
+ }
+ }
+ }
+
+ call msifree (msi)
+
+ # Compute the ellipticity and position angle from the moments.
+ r = (xx + yy)
+ if (r > 0.) {
+ r1 = (xx - yy) / r
+ r2 = 2 * xy / r
+ SFD_E(sfd) = sqrt (r1**2 + r2**2)
+ SFD_PA(sfd) = RADTODEG (atan2 (r2, r1) / 2.)
+ } else {
+ SFD_E(sfd) = 0.
+ SFD_PA(sfd) = 0.
+ }
+
+ # The magnitude and profile normalization is from the max enclosed flux.
+ call alimr (Memr[profile], np, r, SFD_M(sfd))
+ if (SFD_M(sfd) <= 0.)
+ call error (1, "Invalid flux profile")
+ call adivkr (Memr[profile], SFD_M(sfd), Memr[profile], np)
+
+ # Fit interpolator to the enclosed flux profile.
+ call asiinit (asi, II_SPLINE3)
+ call asifit (asi, Memr[profile], np)
+ SFD_ASI1(sfd) = asi
+
+ # Estimate a gaussian sigma (actually sqrt(2)*sigma) and if it is
+ # it is small subtract the gaussian so that the image interpolator
+ # can more accurately estimate subpixel values.
+
+ #call stf_radius (sf, sfd, SF_LEVEL(sf), r)
+ #sigma = r / sqrt (log (1/(1-SF_LEVEL(sf))))
+ call stf_radius (sf, sfd, 0.8, r)
+ r = r / SF_SCALE(sf)
+ sigma = 2 * r * sqrt (log(2.) / log (1/(1-0.8)))
+ if (sigma < 5.) {
+ if (sigma <= 2.) {
+ call gsrestore (gs, gsdata)
+ dx = xc - nint (xc)
+ dy = yc - nint (yc)
+ r = sqrt (dx * dx + dy * dy)
+ dx = 1.
+ ds = abs (sigma - gseval (gs, r, dx))
+ for (da = 1.; da <= 2.; da = da + .01) {
+ dz = abs (sigma - gseval (gs, r, da))
+ if (dz < ds) {
+ ds = dz
+ dx = da
+ }
+ }
+ sigma = dx
+ call gsfree (gs)
+ }
+
+ sigma = sigma / (2 * sqrt (log(2.)))
+ sigma = sigma * sigma
+
+ # Compute the peak that gives the correct central pixel value.
+ i = nint (xc)
+ j = nint (yc)
+ dx = i - xc
+ dy = j - yc
+ r = sqrt (dx * dx + dy * dy)
+ ns = stf_r2n (r)
+ ds = 1. / ns
+ da = ds * ds
+ dz = 0.5 + 0.5 * ds
+
+ r1 = 0.
+ dy1 = dy - dz
+ do l = 1, ns {
+ dy1 = dy1 + ds
+ dy2 = dy1 * dy1
+ dx1 = dx - dz
+ do k = 1, ns {
+ dx1 = dx1 + ds
+ dx2 = dx1 * dx1
+ r2 = (dx2 + dy2) / sigma
+ if (r2 < 25.)
+ r1 = r1 + exp (-r2)
+ }
+ }
+ ptr = data + (j - y1 + 1) * nx + (i - x1 + 1)
+ peak = (Memr[ptr] - bkgd) / (r1 * da)
+
+ # Subtract the gaussian from the data.
+ do j = y1, y2 {
+ ptr = data + (j - y1 + 1) * nx + 1
+ dy = j - yc
+ do i = x1, x2 {
+ dx = i - xc
+ r = sqrt (dx * dx + dy * dy)
+ ns = stf_r2n (r)
+ ds = 1. / ns
+ da = ds * ds
+ dz = 0.5 + 0.5 * ds
+
+ r1 = 0.
+ dy1 = dy - dz
+ do l = 1, ns {
+ dy1 = dy1 + ds
+ dy2 = dy1 * dy1
+ dx1 = dx - dz
+ do k = 1, ns {
+ dx1 = dx1 + ds
+ dx2 = dx1 * dx1
+ r2 = (dx2 + dy2) / sigma
+ if (r2 < 25.)
+ r1 = r1 + peak * exp (-r2)
+ }
+ }
+ Memr[ptr] = Memr[ptr] - r1 * da
+ ptr = ptr + 1
+ }
+ }
+
+ # Fit the image interpolator to the residual data.
+ call msiinit (msi, II_BISPLINE3)
+ call msifit (msi, Memr[data], nx, ny, nx)
+
+ # Recompute the enclosed flux profile and moments
+ # using the gaussian plus image interpolator fit to the residuals.
+
+ call aclrr (Memr[profile], np)
+
+ xx = 0.
+ yy = 0.
+ xy = 0.
+ do j = y1, y2 {
+ ptr = data + (j - y1 + 1) * nx + 1
+ dy = j - yc
+ do i = x1, x2 {
+ dx = i - xc
+ r = sqrt (dx * dx + dy * dy)
+ ns = stf_r2n (r)
+ ds = 1. / ns
+ da = ds * ds
+ dz = 0.5 + 0.5 * ds
+
+ # Compute interpolator correction.
+ r2 = 0.
+ dy1 = dy - dz
+ do l = 1, ns {
+ dy1 = dy1 + ds
+ dx1 = dx - dz
+ do k = 1, ns {
+ dx1 = dx1 + ds
+ r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2)
+ r2 = r2 + r1
+ }
+ }
+
+ r1 = Memr[ptr] - bkgd
+ ptr = ptr + 1
+ r2 = r1 - r2 * da
+
+ # Accumulate the enclosed flux and moments.
+ dy1 = dy - dz
+ do l = 1, ns {
+ dy1 = dy1 + ds
+ dy2 = dy1 * dy1
+ dx1 = dx - dz
+ do k = 1, ns {
+ dx1 = dx1 + ds
+ dx2 = dx1 * dx1
+ r3 = (dx2 + dy2) / sigma
+ if (r3 < 25.)
+ r3 = peak * exp (-r3)
+ else
+ r3 = 0.
+ r = max (0., sqrt (dx2 + dy2) - ds / 2)
+ if (r < radius) {
+ r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2)
+ r1 = da * (r1 + r2 + r3)
+
+ for (m=stf_r2i(r)+1; m<=np; m=m+1) {
+ r3 = (stf_i2r (real(m)) - r) / ds
+ if (r3 >= 1.)
+ break
+ Memr[profile+m-1] = Memr[profile+m-1] +
+ r3 * r1
+ }
+ for (; m<=np; m=m+1)
+ Memr[profile+m-1] = Memr[profile+m-1] + r1
+
+ if (r1 > miso) {
+ xx = xx + dx2 * r1
+ yy = yy + dy2 * r1
+ xy = xy + dx1 * dy1 * r1
+ }
+ }
+ }
+ }
+ }
+ }
+
+ call msifree (msi)
+
+ # Recompute the moments, magnitude, normalized flux, and interp.
+ r = (xx + yy)
+ if (r > 0.) {
+ r1 = (xx - yy) / r
+ r2 = 2 * xy / r
+ SFD_E(sfd) = sqrt (r1**2 + r2**2)
+ SFD_PA(sfd) = RADTODEG (atan2 (r2, r1) / 2.)
+ } else {
+ SFD_E(sfd) = 0.
+ SFD_PA(sfd) = 0.
+ }
+
+ call alimr (Memr[profile], np, r, SFD_M(sfd))
+ if (SFD_M(sfd) <= 0.)
+ call error (1, "Invalid flux profile")
+ call adivkr (Memr[profile], SFD_M(sfd), Memr[profile], np)
+
+ call asifit (asi, Memr[profile], np)
+ SFD_ASI1(sfd) = asi
+ }
+
+ # Compute derivative of enclosed flux profile and fit an image
+ # interpolator.
+
+ dx = 0.25
+ Memr[profile] = 0.
+ ns = 0
+ do i = 1, np {
+ r = stf_i2r (real(i))
+ r2 = stf_r2i (r + dx)
+ if (r2 > np) {
+ k = i
+ break
+ }
+ r1 = stf_r2i (r - dx)
+ if (r1 < 1) {
+ if (i > 1) {
+ dy = asieval (asi, real(i)) / r**2
+ Memr[profile] = (ns * Memr[profile] + dy) / (ns + 1)
+ ns = ns + 1
+ }
+ j = i
+ } else {
+ dy = (asieval (asi, r2) - asieval (asi, r1)) /
+ (4 * r * dx)
+ Memr[profile+i-1] = dy
+ }
+ }
+ do i = 2, j
+ Memr[profile+i-1] = (Memr[profile+j] - Memr[profile]) / j *
+ (i - 1) + Memr[profile]
+ do i = k, np
+ Memr[profile+i-1] = Memr[profile+k-2]
+
+ call adivkr (Memr[profile], SF_SCALE(sf)**2, Memr[profile], np)
+ call alimr (Memr[profile], np, SFD_YP1(sfd), SFD_YP2(sfd))
+ call asiinit (asi, II_SPLINE3)
+ call asifit (asi, Memr[profile], np)
+ SFD_ASI2(sfd) = asi
+ #SF_XP1(sf) = j+1
+ SF_XP1(sf) = 1
+ SF_XP2(sf) = k-1
+
+ call sfree (sp)
+end
+
+
+# STF_NORM -- Renormalize the enclosed flux profile.
+
+procedure stf_norm (sf, sfd, x, y)
+
+pointer sf #I Parameter structure
+pointer sfd #I Star structure
+real x #I Radius
+real y #I Flux
+
+int npmax, np
+pointer asi
+
+int i, j, k
+real r, r1, r2, dx, dy
+pointer sp, profile
+real asieval(), stf_i2r(), stf_r2i()
+errchk asifit
+
+begin
+ npmax = SFD_NPMAX(sfd)
+ np = SFD_NP(sfd)
+ asi = SFD_ASI1(sfd)
+
+ call smark (sp)
+ call salloc (profile, npmax, TY_REAL)
+
+ # Renormalize the enclosed flux profile.
+ if (IS_INDEF(x) || x <= 0.) {
+ dy = SFD_BKGD(sfd) - SFD_BKGD1(sfd)
+ SFD_BKGD(sfd) = SFD_BKGD(sfd) - dy
+ do i = 1, npmax
+ Memr[profile+i-1] = asieval (asi, real(i)) +
+ dy * stf_i2r(real(i)) ** 2
+ call alimr (Memr[profile], np, r1, r2)
+ call adivkr (Memr[profile], r2, Memr[profile], npmax)
+ } else if (IS_INDEF(y)) {
+ r = max (1., min (real(np), stf_r2i (x)))
+ r2 = asieval (asi, r)
+ if (r2 <= 0.)
+ return
+ do i = 1, npmax
+ Memr[profile+i-1] = asieval (asi, real(i))
+ call adivkr (Memr[profile], r2, Memr[profile], npmax)
+ } else {
+ r = max (1., min (real(np), stf_r2i (x)))
+ r1 = asieval (asi, r)
+ dy = (y - r1) / x ** 2
+ SFD_BKGD(sfd) = SFD_BKGD(sfd) - dy
+ do i = 1, npmax
+ Memr[profile+i-1] = asieval (asi, real(i)) +
+ dy * stf_i2r(real(i)) ** 2
+ }
+
+ call asifit (asi, Memr[profile], npmax)
+ SFD_ASI1(sfd) = asi
+
+ # Compute derivative of enclosed flux profile and fit an image
+ # interpolator.
+
+ dx = 0.25
+ do i = 1, npmax {
+ r = stf_i2r (real(i))
+ r2 = stf_r2i (r + dx)
+ if (r2 > np) {
+ k = i
+ break
+ }
+ r1 = stf_r2i (r - dx)
+ if (r1 < 1) {
+ if (i > 1) {
+ dy = asieval (asi, real(i)) / r**2
+ Memr[profile] = dy
+ }
+ j = i
+ } else {
+ dy = (asieval (asi, r2) - asieval (asi, r1)) /
+ (4 * r * dx)
+ Memr[profile+i-1] = dy
+ }
+ }
+ do i = 2, j
+ Memr[profile+i-1] = (Memr[profile+j] - Memr[profile]) / j *
+ (i - 1) + Memr[profile]
+ do i = k, npmax
+ Memr[profile+i-1] = Memr[profile+k-2]
+
+ call adivkr (Memr[profile], SF_SCALE(sf)**2, Memr[profile], np)
+ call alimr (Memr[profile], np, SFD_YP1(sfd), SFD_YP2(sfd))
+ asi = SFD_ASI2(sfd)
+ call asifit (asi, Memr[profile], np)
+ SFD_ASI2(sfd) = asi
+ #SF_XP1(sf) = min (j+1, np)
+ SF_XP1(sf) = 1
+ SF_XP2(sf) = min (k-1, np)
+
+ call sfree (sp)
+end
+
+
+# STF_WIDTHS -- Set the widhts.
+
+procedure stf_widths (sf, sfd)
+
+pointer sf #I Main data structure
+pointer sfd #I Star data structure
+
+errchk stf_radius, stf_dfwhm, stf_fit
+
+begin
+ call stf_radius (sf, sfd, SF_LEVEL(sf), SFD_R(sfd))
+ call stf_dfwhm (sf, sfd)
+ call stf_fit (sf, sfd)
+
+ switch (SF_WCODE(sf)) {
+ case 1:
+ SFD_W(sfd) = SFD_R(sfd)
+ case 2:
+ SFD_W(sfd) = SFD_DFWHM(sfd)
+ case 3:
+ SFD_W(sfd) = SFD_GFWHM(sfd)
+ case 4:
+ SFD_W(sfd) = SFD_MFWHM(sfd)
+ }
+end
+
+
+# STF_I2R -- Compute radius from sample index.
+
+real procedure stf_i2r (i)
+
+real i #I Index
+real r #O Radius
+
+begin
+ if (i < 20)
+ r = 0.05 * i
+ else if (i < 30)
+ r = 0.1 * i - 1
+ else if (i < 40)
+ r = 0.2 * i - 4
+ else if (i < 50)
+ r = 0.5 * i - 16
+ else
+ r = i - 41
+ return (r)
+end
+
+
+# STF_R2I -- Compute sample index from radius.
+
+real procedure stf_r2i (r)
+
+real r #I Radius
+real i #O Index
+
+begin
+ if (r < 1)
+ i = 20 * r
+ else if (r < 2)
+ i = 10 * (r + 1)
+ else if (r < 4)
+ i = 5 * (r + 4)
+ else if (r < 9)
+ i = 2 * (r + 16)
+ else
+ i = r + 41
+ return (i)
+end
+
+
+# STF_R2N -- Compute number of subsamples from radius.
+
+int procedure stf_r2n (r)
+
+real r #I Radius
+int n #O Number of subsamples
+
+begin
+ if (r < 1)
+ n = 20
+ else if (r < 2)
+ n = 10
+ else if (r < 4)
+ n = 5
+ else if (r < 9)
+ n = 2
+ else
+ n = 1
+ return (n)
+end
+
+
+# STF_MODEL -- Return model value.
+
+procedure stf_model (sf, sfd, r, profile, flux)
+
+pointer sf #I Main data structure
+pointer sfd #I Star data structure
+real r #I Radius at level
+real profile #I Profile value
+real flux #I Enclosed flux value
+
+real x, x1, x2, r1, r2, dr
+
+begin
+ dr = 0.25 * SF_SCALE(sf)
+ r1 = r - dr
+ r2 = r + dr
+ if (r1 < 0.) {
+ r1 = dr
+ r2 = r1 + dr
+ }
+
+ switch (SF_WCODE(sf)) {
+ case 3:
+ x = r**2 / (2. * SFD_SIGMA(sfd)**2)
+ if (x < 20.)
+ flux = 1 - exp (-x)
+ else
+ flux = 0.
+
+ x1 = r1**2 / (2. * SFD_SIGMA(sfd)**2)
+ x2 = r2**2 / (2. * SFD_SIGMA(sfd)**2)
+ if (x2 < 20.) {
+ x1 = 1 - exp (-x1)
+ x2 = 1 - exp (-x2)
+ } else {
+ x1 = 1.
+ x2 = 1.
+ }
+ if (r <= dr) {
+ x1 = x1 / dr ** 2
+ x2 = x2 / (4 * dr ** 2)
+ profile = (x2 - x1) / dr * r + x1
+ } else {
+ profile = (x2 - x1) / (4 * r * dr)
+ }
+ default:
+ x = 1 + (r / SFD_ALPHA(sfd)) ** 2
+ flux = 1 - x ** (1 - SFD_BETA(sfd))
+
+ x1 = 1 + (r1 / SFD_ALPHA(sfd)) ** 2
+ x2 = 1 + (r2 / SFD_ALPHA(sfd)) ** 2
+ x1 = 1 - x1 ** (1 - SFD_BETA(sfd))
+ x2 = 1 - x2 ** (1 - SFD_BETA(sfd))
+ if (r <= dr) {
+ x1 = x1 / dr ** 2
+ x2 = x2 / (4 * dr ** 2)
+ profile = (x2 - x1) / dr * r + x1
+ } else {
+ profile = (x2 - x1) / (4 * r * dr)
+ }
+ }
+end
+
+
+# STF_DFWHM -- Direct FWHM from profile.
+
+procedure stf_dfwhm (sf, sfd)
+
+pointer sf #I Main data structure
+pointer sfd #I Star data structure
+
+int np
+real r, rpeak, profile, peak, asieval(), stf_i2r()
+pointer asi
+
+begin
+ asi = SFD_ASI2(sfd)
+ np = SFD_NP(sfd)
+
+ rpeak = 1.
+ peak = 0.
+ for (r=1.; r <= np; r = r + 0.01) {
+ profile = asieval (asi, r)
+ if (profile > peak) {
+ rpeak = r
+ peak = profile
+ }
+ }
+
+ peak = peak / 2.
+ for (r=rpeak; r <= np && asieval (asi, r) > peak; r = r + 0.01)
+ ;
+
+ SFD_DFWHM(sfd) = 2 * stf_i2r (r) * SF_SCALE(sf)
+end
+
+
+# STF_FWHMS -- Measure FWHM vs level.
+
+procedure stf_fwhms (sf, sfd)
+
+pointer sf #I Main data structure
+pointer sfd #I Star data structure
+
+int i
+real level, r
+
+begin
+ do i = 1, 19 {
+ level = i * 0.05
+ call stf_radius (sf, sfd, level, r)
+ switch (SF_WCODE(sf)) {
+ case 3:
+ SFD_FWHM(sfd,i) = 2 * r * sqrt (log (2.) / log (1/(1-level)))
+ default:
+ r = r / sqrt ((1.-level)**(1./(1.-SFD_BETA(sfd))) - 1.)
+ SFD_FWHM(sfd,i) = 2 * r * sqrt (2.**(1./SFD_BETA(sfd))-1.)
+ }
+ }
+end
+
+
+# STF_RADIUS -- Measure the radius at the specified level.
+
+procedure stf_radius (sf, sfd, level, r)
+
+pointer sf #I Main data structure
+pointer sfd #I Star data structure
+real level #I Level to measure
+real r #O Radius
+
+int np
+pointer asi
+real f, fmax, rmax, asieval(), stf_i2r()
+
+begin
+ np = SFD_NP(sfd)
+ asi = SFD_ASI1(sfd)
+
+ for (r=1; r <= np && asieval (asi, r) < level; r = r + 0.01)
+ ;
+ if (r > np) {
+ fmax = 0.
+ rmax = 0.
+ for (r=1; r <= np; r = r + 0.01) {
+ f = asieval (asi, r)
+ if (f > fmax) {
+ fmax = f
+ rmax = r
+ }
+ }
+ r = rmax
+ }
+ r = stf_i2r (r) * SF_SCALE(sf)
+end
+
+
+# STF_FIT -- Fit models to enclosed flux.
+
+procedure stf_fit (sf, sfd)
+
+pointer sf #I Main data structure
+pointer sfd #I Star data structure
+
+int i, j, n, np, pfit[2]
+real beta, z, params[3]
+pointer asi, nl
+pointer sp, x, y, w
+
+int locpr()
+real asieval(), stf_i2r()
+extern stf_gauss1(), stf_gauss2(), stf_moffat1(), stf_moffat2()
+errchk nlinitr, nlfitr
+
+data pfit/2,3/
+
+begin
+ np = SFD_NP(sfd)
+ asi = SFD_ASI1(sfd)
+
+ call smark (sp)
+ call salloc (x, np, TY_REAL)
+ call salloc (y, np, TY_REAL)
+ call salloc (w, np, TY_REAL)
+
+ n = 0
+ j = 0
+ do i = 1, np {
+ z = 1. - max (0., asieval (asi, real(i)))
+ if (n > np/3 && z < 0.5)
+ break
+ if ((n < np/3 && z > 0.01) || z > 0.5)
+ j = n
+
+ Memr[x+n] = stf_i2r (real(i)) * SF_SCALE(sf)
+ Memr[y+n] = z
+ Memr[w+n] = 1.
+ n = n + 1
+ }
+
+ # Gaussian.
+ np = 1
+ params[2] = Memr[x+j] / sqrt (2. * log (1./min(0.99,Memr[y+j])))
+ params[1] = 1
+ call nlinitr (nl, locpr (stf_gauss1), locpr (stf_gauss2),
+ params, params, 2, pfit, np, .001, 100)
+ call nlfitr (nl, Memr[x], Memr[y], Memr[w], n, 1, WTS_USER, i)
+ if (i != SINGULAR && i != NO_DEG_FREEDOM) {
+ call nlpgetr (nl, params, i)
+ if (params[2] < 0.)
+ params[2] = Memr[x+j] / sqrt (2. * log (1./min(0.99,Memr[y+j])))
+ }
+ SFD_SIGMA(sfd) = params[2]
+ SFD_GFWHM(sfd) = 2 * SFD_SIGMA(sfd) * sqrt (2. * log (2.))
+
+ # Moffat.
+ if (SF_BETA(sf) < 1.1) {
+ call nlfreer (nl)
+ call sfree (sp)
+ call error (1, "Cannot measure FWHM - Moffat beta too small")
+ }
+
+ beta = SF_BETA(sf)
+ if (IS_INDEFR(beta)) {
+ beta = 2.5
+ np = 2
+ } else {
+ np = 1
+ }
+ params[3] = 1 - beta
+ params[2] = Memr[x+j] / sqrt (min(0.99,Memr[y+j])**(1./params[3]) - 1.)
+ params[1] = 1
+ call nlinitr (nl, locpr (stf_moffat1), locpr (stf_moffat2),
+ params, params, 3, pfit, np, .001, 100)
+ call nlfitr (nl, Memr[x], Memr[y], Memr[w], n, 1, WTS_USER, i)
+ if (i != SINGULAR && i != NO_DEG_FREEDOM) {
+ call nlpgetr (nl, params, i)
+ if (params[2] < 0.) {
+ params[3] = 1. - beta
+ params[2] = Memr[x+j] /
+ sqrt (min(0.99,Memr[y+j])**(1./params[3]) - 1.)
+ }
+ }
+ SFD_ALPHA(sfd) = params[2]
+ SFD_BETA(sfd) = 1 - params[3]
+ SFD_MFWHM(sfd) = 2 * SFD_ALPHA(sfd) * sqrt (2.**(1./SFD_BETA(sfd))-1.)
+
+ call nlfreer (nl)
+ call sfree (sp)
+end
+
+
+# STF_GAUSS1 -- Gaussian function used in NLFIT. The parameters are the
+# amplitude and sigma and the input variable is the radius.
+
+procedure stf_gauss1 (x, nvars, p, np, z)
+
+real x[nvars] #I Input variables
+int nvars #I Number of variables
+real p[np] #I Parameter vector
+int np #I Number of parameters
+real z #O Function return
+
+real r2
+
+begin
+ r2 = x[1]**2 / (2 * p[2]**2)
+ if (abs (r2) > 20.)
+ z = 0.
+ else
+ z = p[1] * exp (-r2)
+end
+
+
+# STF_GAUSS2 -- Gaussian function and derivatives used in NLFIT. The parameters
+# are the amplitude and sigma and the input variable is the radius.
+
+procedure stf_gauss2 (x, nvars, p, dp, np, z, der)
+
+real x[nvars] #I Input variables
+int nvars #I Number of variables
+real p[np] #I Parameter vector
+real dp[np] #I Dummy array of parameters increments
+int np #I Number of parameters
+real z #O Function return
+real der[np] #O Derivatives
+
+real r2
+
+begin
+ r2 = x[1]**2 / (2 * p[2]**2)
+ if (abs (r2) > 20.) {
+ z = 0.
+ der[1] = 0.
+ der[2] = 0.
+ } else {
+ der[1] = exp (-r2)
+ z = p[1] * der[1]
+ der[2] = z * 2 * r2 / p[2]
+ }
+end
+
+
+# STF_MOFFAT1 -- Moffat function used in NLFIT. The parameters are the
+# amplitude, alpha squared, and beta and the input variable is the radius.
+
+procedure stf_moffat1 (x, nvars, p, np, z)
+
+real x[nvars] #I Input variables
+int nvars #I Number of variables
+real p[np] #I Parameter vector
+int np #I Number of parameters
+real z #O Function return
+
+real y
+
+begin
+ y = 1 + (x[1] / p[2]) ** 2
+ if (abs (y) > 20.)
+ z = 0.
+ else
+ z = p[1] * y ** p[3]
+end
+
+
+# STF_MOFFAT2 -- Moffat function and derivatives used in NLFIT. The
+# parameters are the amplitude, alpha squared, and beta and the input
+# variable is the radius.
+
+procedure stf_moffat2 (x, nvars, p, dp, np, z, der)
+
+real x[nvars] #I Input variables
+int nvars #I Number of variables
+real p[np] #I Parameter vector
+real dp[np] #I Dummy array of parameters increments
+int np #I Number of parameters
+real z #O Function return
+real der[np] #O Derivatives
+
+real y
+
+begin
+ y = 1 + (x[1] / p[2]) ** 2
+ if (abs (y) > 20.) {
+ z = 0.
+ der[1] = 0.
+ der[2] = 0.
+ der[3] = 0.
+ } else {
+ der[1] = y ** p[3]
+ z = p[1] * der[1]
+ der[2] = -2 * z / y * p[3] / p[2] * (x[1] / p[2]) ** 2
+ der[3] = z * log (y)
+ }
+end
diff --git a/noao/obsutil/src/starfocus/t_starfocus.x b/noao/obsutil/src/starfocus/t_starfocus.x
new file mode 100644
index 00000000..2c3213e9
--- /dev/null
+++ b/noao/obsutil/src/starfocus/t_starfocus.x
@@ -0,0 +1,1240 @@
+include <error.h>
+include <imhdr.h>
+include <imset.h>
+include <mach.h>
+include "starfocus.h"
+
+define HELP "nmisc$src/starfocus.key"
+define PROMPT "Options"
+
+
+# T_STARFOCUS -- Stellar focusing task.
+
+procedure t_starfocus ()
+
+begin
+ call starfocus (STARFOCUS)
+end
+
+
+# T_PSFMEASURE -- PSF measuring task.
+
+procedure t_psfmeasure ()
+
+begin
+ call starfocus (PSFMEASURE)
+end
+
+
+# STARFOCUS -- Stellar focusing and PSF measuring main routine.
+
+procedure starfocus (type)
+
+int type #I Task type
+
+int list # List of images
+pointer fvals # Focus values
+pointer fstep # Focus step
+pointer nexposures # Number of exposures
+pointer step # step in pixels
+int direction # Step direction
+int gap # Double step gap
+int coords # Type of image data
+bool display # Display images?
+int frame # Display frame
+int logfd # Log file descriptor
+bool ignore_sat # Ignore saturation?
+
+real wx, wy, f, df, xstep, ystep
+int i, i1, i2, i3, j, k, l, ip, wcs, key, id, ncols, nlines
+int nexp, nsfd, nimages, nstars, ngraph, nmark
+pointer sp, sf, image, system, cmd, rg, mark, im, mw, ct
+pointer sfds, sfd
+
+bool clgetb(), streq()
+real clgetr(), imgetr(), stf_r2i()
+int clgeti(), clgwrd(), clgcur(), imtopenp(), imtgetim(), imgeti()
+int nowhite(), open(), rng_index(), strdic(), ctoi(), ctor()
+pointer rng_open(), immap(), mw_openim(), mw_sctran()
+errchk immap, open, imgetr, imgeti, mw_openim, mw_sctran
+errchk stf_find, stf_bkgd, stf_profile, stf_widths, stf_fwhms, stf_radius
+errchk stf_organize, stf_graph, stf_display
+
+begin
+ call smark (sp)
+ call salloc (sf, SF, TY_STRUCT)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (fvals, SZ_LINE, TY_CHAR)
+ call salloc (fstep, SZ_LINE, TY_CHAR)
+ call salloc (nexposures, SZ_LINE, TY_CHAR)
+ call salloc (step, SZ_LINE, TY_CHAR)
+ call salloc (system, SZ_LINE, TY_CHAR)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ call aclri (Memi[sf], SF)
+ SF_TASK(sf) = type
+
+ # Set task parameters.
+ switch (type) {
+ case STARFOCUS:
+ call clgstr ("focus", Memc[fvals], SZ_LINE)
+ call clgstr ("fstep", Memc[fstep], SZ_LINE)
+ call clgstr ("nexposures", Memc[nexposures], SZ_LINE)
+ call clgstr ("step", Memc[step], SZ_LINE)
+
+ direction = clgwrd ("direction", Memc[cmd], SZ_LINE,
+ "|-line|+line+|-column|+column|")
+ gap = clgwrd ("gap", Memc[cmd], SZ_LINE, "|none|beginning|end|")
+
+ if (nowhite (Memc[fvals], Memc[fvals], SZ_LINE) != 0) {
+ iferr (rg = rng_open (Memc[fvals], -MAX_REAL, MAX_REAL, 1.))
+ rg = NULL
+ } else
+ rg = NULL
+ case PSFMEASURE:
+ Memc[fvals] = EOS
+ rg = NULL
+ nexp = 1
+ }
+
+ list = imtopenp ("images")
+ display = clgetb ("display")
+ frame = clgeti ("frame")
+ coords = clgwrd ("coords", Memc[cmd], SZ_LINE, SF_TYPES)
+ call clgstr ("wcs", Memc[system], SZ_LINE)
+
+ SF_XF(sf) = clgetr ("xcenter")
+ SF_YF(sf) = clgetr ("ycenter")
+ SF_LEVEL(sf) = clgetr ("level")
+ SF_WCODE(sf) = clgwrd ("size", SF_WTYPE(sf), SF_SZWTYPE, SF_WTYPES)
+ SF_BETA(sf) = clgetr ("beta")
+ SF_SCALE(sf) = clgetr ("scale")
+ SF_RADIUS(sf) = max (3., clgetr ("radius"))
+ SF_NIT(sf) = clgeti ("iterations")
+ SF_SBUF(sf) = clgetr ("sbuffer")
+ SF_SWIDTH(sf) = clgetr ("swidth")
+ SF_SAT(sf) = clgetr ("saturation")
+ ignore_sat = clgetb ("ignore_sat")
+ SF_OVRPLT(sf) = NO
+
+ if (SF_LEVEL(sf) > 1.)
+ SF_LEVEL(sf) = SF_LEVEL(sf) / 100.
+ SF_LEVEL(sf) = max (0.05, min (0.95, SF_LEVEL(sf)))
+
+ # Accumulate the psf/focus data.
+ key = 'm'
+ mark = NULL
+ nstars = 0
+ nmark = 0
+ ngraph = 0
+ nimages = 0
+ nsfd = 0
+ while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
+ im = immap (Memc[image], READ_ONLY, 0)
+ call imseti (im, IM_TYBNDRY, TYBNDRY)
+ call imseti (im, IM_NBNDRYPIX, NBNDRYPIX)
+ if (streq (Memc[system], "logical")) {
+ mw = NULL
+ ct = NULL
+ } else {
+ mw = mw_openim (im)
+ ct = mw_sctran (mw, Memc[system], "logical", 03)
+ }
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+ nimages = nimages + 1
+ if (nimages == 1) {
+ SF_NCOLS(sf) = ncols
+ SF_NLINES(sf) = nlines
+ if (IS_INDEF(SF_XF(sf)))
+ SF_XF(sf) = (SF_NCOLS(sf) + 1) / 2.
+ if (IS_INDEF(SF_YF(sf)))
+ SF_YF(sf) = (SF_NLINES(sf) + 1) / 2.
+ } else if (ncols!=SF_NCOLS(sf)||nlines!=SF_NLINES(sf))
+ call eprintf ("WARNING: Images have different sizes\n")
+
+ # Display the image if needed.
+ if (display) {
+ switch (coords) {
+ case SF_MARK1:
+ if (nimages == 1)
+ call stf_display (Memc[image], frame)
+ case SF_MARKALL:
+ call stf_display (Memc[image], frame)
+ }
+ if (nimages == 1) {
+ call printf (
+ "** Select stars to measure with 'm' and finish with 'q'.\n")
+ call printf (
+ "** Additional options are '?', 'g', and :show.\n")
+ call flush (STDOUT)
+ }
+ }
+
+ # Accumulate objects.
+ repeat {
+ switch (coords) {
+ case SF_CENTER:
+ if (nstars == nimages)
+ break
+ if (rg == NULL && Memc[fvals] == EOS)
+ id = nstars
+ else
+ id = 0
+ wx = 1 + (ncols - 1) / 2.
+ wy = 1 + (nlines - 1) / 2.
+ key = 0
+ case SF_MARK1:
+ if (nimages == 1) {
+ if (clgcur ("imagecur", wx, wy, wcs, key,
+ Memc[cmd], SZ_LINE) == EOF)
+ break
+ switch (key) {
+ case '?':
+ call pagefile (HELP, PROMPT)
+ next
+ case ':':
+ if (strdic (Memc[cmd], Memc[cmd], SZ_LINE,
+ "|show|") == 1) {
+ if (nsfd > 0) {
+ call stf_organize (sf, sfds, nsfd)
+ call mktemp ("tmp$iraf", Memc[cmd], SZ_LINE)
+ logfd = open (Memc[cmd], APPEND, TEXT_FILE)
+ call stf_log (sf, logfd)
+ call close (logfd)
+ call pagefile (Memc[cmd], "starfocus")
+ call delete (Memc[cmd])
+ }
+ }
+ next
+ case 'q':
+ break
+ }
+ id = nstars
+ if (mark == NULL)
+ call malloc (mark, 3*10, TY_REAL)
+ else if (mod (nmark, 10) == 0)
+ call realloc (mark, 3*(nmark+10), TY_REAL)
+ Memr[mark+3*nmark] = id
+ Memr[mark+3*nmark+1] = wx
+ Memr[mark+3*nmark+2] = wy
+ nmark = nmark+1
+ } else {
+ if (nmark == 0)
+ break
+ if (nstars / nmark == nimages)
+ break
+ i = mod (nstars, nmark)
+ id = Memr[mark+3*i]
+ wx = Memr[mark+3*i+1]
+ wy = Memr[mark+3*i+2]
+ key = 0
+ }
+ if (ct != NULL)
+ call mw_c2tranr (ct, wx, wy, wx, wy)
+ case SF_MARKALL:
+ if (clgcur ("imagecur", wx, wy, wcs, key,
+ Memc[cmd], SZ_LINE) == EOF)
+ break
+ switch (key) {
+ case '?':
+ call pagefile (HELP, PROMPT)
+ next
+ case ':':
+ if (strdic(Memc[cmd],Memc[cmd],SZ_LINE,"|show|")==1) {
+ if (nsfd > 0) {
+ call stf_organize (sf, sfds, nsfd)
+ call mktemp ("tmp$iraf", Memc[cmd], SZ_LINE)
+ logfd = open (Memc[cmd], APPEND, TEXT_FILE)
+ call stf_log (sf, logfd)
+ call close (logfd)
+ call pagefile (Memc[cmd], "starfocus")
+ call delete (Memc[cmd])
+ }
+ }
+ next
+ case 'q':
+ break
+ }
+ id = nstars
+ if (ct != NULL)
+ call mw_c2tranr (ct, wx, wy, wx, wy)
+ }
+
+ if (type == STARFOCUS) {
+ ip = 1
+ if (ctoi (Memc[nexposures], ip, nexp) == 0)
+ nexp = imgeti (im, Memc[nexposures])
+ ip = 1
+ if (ctor (Memc[step], ip, f) == 0)
+ f = imgetr (im, Memc[step])
+
+ xstep = 0.
+ ystep = 0.
+ switch (direction) {
+ case 1:
+ ystep = -f
+ case 2:
+ ystep = f
+ case 3:
+ xstep = -f
+ case 4:
+ xstep = f
+ }
+
+ # Set the steps and order of evaluation.
+ # This depends on which star in the sequence is marked.
+ # Below we assume the minimum x or maximum y is marked.
+
+ i1 = 1; i2 = nexp; i3 = 1
+# if (xstep < 0.) {
+# i1 = nexp; i2 = 1; i3 = -1; xstep = -xstep
+# }
+# if (ystep > 0.) {
+# i1 = nexp; i2 = 1; i3 = -1; ystep = -ystep
+# }
+ } else {
+ i1 = 1; i2 = 1; i3 = 1
+ }
+
+ k = nsfd
+ do i = i1, i2, i3 {
+ if (i > 1) {
+ wx = wx + xstep
+ wy = wy + ystep
+ switch (gap) {
+ case 2:
+ if ((i==2 && i3==1) || (i==1 && i3==-1)) {
+ wx = wx + xstep
+ wy = wy + ystep
+ }
+ case 3:
+ if ((i==nexp && i3==1) || (i==nexp-1 && i3==-1)) {
+ wx = wx + xstep
+ wy = wy + ystep
+ }
+ }
+ }
+
+ if (wx < SF_RADIUS(sf)-NBNDRYPIX ||
+ wx > IM_LEN(im,1)-SF_RADIUS(sf)+NBNDRYPIX ||
+ wy < SF_RADIUS(sf)-NBNDRYPIX ||
+ wy > IM_LEN(im,2)-SF_RADIUS(sf)+NBNDRYPIX)
+ next
+ if (nexp == 1)
+ j = nimages
+ else
+ j = i
+ if (nsfd == 0)
+ call malloc (sfds, 10, TY_POINTER)
+ else if (mod (nsfd, 10) == 0)
+ call realloc (sfds, nsfd+10, TY_POINTER)
+ call malloc (sfd, SFD, TY_STRUCT)
+ call strcpy (Memc[image], SFD_IMAGE(sfd), SF_SZFNAME)
+ SFD_ID(sfd) = id
+ SFD_X(sfd) = wx
+ SFD_Y(sfd) = wy
+ if (Memc[fvals] == EOS)
+ f = INDEF
+ else if (Memc[fstep] != EOS) {
+ ip = 1
+ if (ctor (Memc[fvals], ip, f) == 0)
+ f = imgetr (im, Memc[fvals])
+ ip = 1
+ if (ctor (Memc[fstep], ip, df) == 0)
+ df = imgetr (im, Memc[fstep])
+ f = f + (i - 1) * df
+ } else if (rg != NULL) {
+ if (rng_index (rg, j, f) == EOF)
+ call error (1, "Focus list ended prematurely")
+ } else
+ f = imgetr (im, Memc[fvals])
+ SFD_F(sfd) = f
+ SFD_STATUS(sfd) = 0
+ SFD_SFS(sfd) = NULL
+ SFD_SFF(sfd) = NULL
+ SFD_SFI(sfd) = NULL
+
+ iferr {
+ do l = 1, SF_NIT(sf) {
+ if (l == 1)
+ SFD_RADIUS(sfd) = max (3., SF_RADIUS(sf))
+ else
+ SFD_RADIUS(sfd) = max (3., 3. * SFD_DFWHM(sfd))
+ SFD_NPMAX(sfd) = stf_r2i (SFD_RADIUS(sfd)) + 1
+ SFD_NP(sfd) = SFD_NPMAX(sfd)
+ call stf_find (sf, sfd, im)
+ call stf_bkgd (sf, sfd)
+ if (SFD_NSAT(sfd) > 0 && l == 1) {
+ if (ignore_sat)
+ call error (0,
+ "Saturated pixels found - ignoring object")
+ else
+ call eprintf (
+ "WARNING: Saturated pixels found.\n")
+ }
+ call stf_profile (sf, sfd)
+ call stf_widths (sf, sfd)
+ call stf_fwhms (sf, sfd)
+ }
+ Memi[sfds+nsfd] = sfd
+ nsfd = nsfd + 1
+ wx = SFD_X(sfd)
+ wy = SFD_Y(sfd)
+ } then {
+ call erract (EA_WARN)
+ call mfree (sfd, TY_STRUCT)
+ }
+ }
+ if (nsfd > k) {
+ nstars = nstars + 1
+ if (key == 'g') {
+ if (nsfd - k > 0) {
+ call stf_organize (sf, sfds+k, nsfd-k)
+ call stf_graph (sf)
+ ngraph = ngraph + 1
+ }
+ }
+ }
+ }
+ if (mw != NULL)
+ call mw_close (mw)
+ call imunmap (im)
+ }
+
+ if (nsfd == 0)
+ call error (1, "No input data")
+
+ # Organize the objects, graph the data, and log the results.
+ if (nstars > 1 || ngraph != nstars) {
+ call stf_organize (sf, sfds, nsfd)
+ call stf_graph (sf)
+ }
+ call stf_log (sf, STDOUT)
+ call clgstr ("logfile", Memc[image], SZ_FNAME)
+ ifnoerr (logfd = open (Memc[image], APPEND, TEXT_FILE)) {
+ call stf_log (sf, logfd)
+ call close (logfd)
+ }
+
+ # Finish up
+ call rng_close (rg)
+ call imtclose (list)
+ call stf_free (sf)
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ call asifree (SFD_ASI1(sfd))
+ call asifree (SFD_ASI2(sfd))
+ call mfree (sfd, TY_STRUCT)
+ }
+ call mfree (SF_SFDS(sf), TY_POINTER)
+ call mfree (mark, TY_REAL)
+ call sfree (sp)
+end
+
+
+# STF_FREE -- Free the starfocus data structures.
+
+procedure stf_free (sf)
+
+pointer sf #I Starfocus structure
+int i
+
+begin
+ do i = 1, SF_NSTARS(sf)
+ call mfree (SF_SFS(sf,i), TY_STRUCT)
+ do i = 1, SF_NFOCUS(sf)
+ call mfree (SF_SFF(sf,i), TY_STRUCT)
+ do i = 1, SF_NIMAGES(sf)
+ call mfree (SF_SFI(sf,i), TY_STRUCT)
+ call mfree (SF_STARS(sf), TY_POINTER)
+ call mfree (SF_FOCUS(sf), TY_POINTER)
+ call mfree (SF_IMAGES(sf), TY_POINTER)
+ SF_NSTARS(sf) = 0
+ SF_NFOCUS(sf) = 0
+ SF_NIMAGES(sf) = 0
+end
+
+
+# STF_ORGANIZE -- Organize the individual object structures by star, focus,
+# and image. Compute focus, radius, and magnitude by group and over all
+# data.
+
+procedure stf_organize (sf, sfds, nsfd)
+
+pointer sf #I Starfocus structure
+pointer sfds #I Pointer to array of object structures
+int nsfd #I Number of object structures
+
+int i, j, k, nstars, nfocus, nimages, key
+real f
+pointer stars, focus, images, sfd, sfs, sff, sfi
+pointer sp, image
+bool streq()
+errchk malloc
+
+int stf_focsort(), stf_magsort()
+extern stf_focsort, stf_magsort
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+
+ # Free previous structures.
+ call stf_free (sf)
+
+ # Organize sfds by star.
+ nstars = 0
+ for (i = 0; i < nsfd; i = i + 1) {
+ key = SFD_ID(Memi[sfds+i])
+ for (j = 0; SFD_ID(Memi[sfds+j]) != key; j = j + 1)
+ ;
+ if (j == i)
+ nstars = nstars + 1
+ }
+ call malloc (stars, nstars, TY_POINTER)
+
+ nstars = 0
+ for (i = 0; i < nsfd; i = i + 1) {
+ key = SFD_ID(Memi[sfds+i])
+ for (j = 0; j < nstars; j = j + 1)
+ if (SFS_ID(Memi[stars+j]) == key)
+ break
+ if (j == nstars) {
+ k = 0
+ for (j = i; j < nsfd; j = j + 1)
+ if (SFD_ID(Memi[sfds+j]) == key)
+ k = k + 1
+ call malloc (sfs, SFS(k), TY_STRUCT)
+ SFS_ID(sfs) = key
+ SFS_NSFD(sfs) = k
+ k = 0
+ for (j = i; j < nsfd; j = j + 1) {
+ sfd = Memi[sfds+j]
+ if (SFD_ID(sfd) == key) {
+ k = k + 1
+ SFD_SFS(sfd) = sfs
+ SFS_SFD(sfs,k) = sfd
+ }
+ }
+ Memi[stars+nstars] = sfs
+ nstars = nstars + 1
+ }
+ }
+
+ # Organize sfds by focus values. Sort by magnitude.
+ nfocus = 0
+ for (i = 0; i < nsfd; i = i + 1) {
+ f = SFD_F(Memi[sfds+i])
+ for (j = 0; SFD_F(Memi[sfds+j]) != f; j = j + 1)
+ ;
+ if (j == i)
+ nfocus = nfocus + 1
+ }
+ call malloc (focus, nfocus, TY_POINTER)
+
+ nfocus = 0
+ for (i = 0; i < nsfd; i = i + 1) {
+ f = SFD_F(Memi[sfds+i])
+ for (j = 0; j < nfocus; j = j + 1)
+ if (SFF_F(Memi[focus+j]) == f)
+ break
+ if (j == nfocus) {
+ k = 0
+ for (j = i; j < nsfd; j = j + 1)
+ if (SFD_F(Memi[sfds+j]) == f)
+ k = k + 1
+ call malloc (sff, SFF(k), TY_STRUCT)
+ SFF_F(sff) = f
+ SFF_NSFD(sff) = k
+ k = 0
+ for (j = i; j < nsfd; j = j + 1) {
+ sfd = Memi[sfds+j]
+ if (SFD_F(sfd) == f) {
+ k = k + 1
+ SFD_SFF(sfd) = sff
+ SFF_SFD(sff,k) = sfd
+ }
+ }
+ Memi[focus+nfocus] = sff
+ nfocus = nfocus + 1
+ }
+ }
+
+ # Organize sfds by image.
+ nimages = 0
+ for (i = 0; i < nsfd; i = i + 1) {
+ call strcpy (SFD_IMAGE(Memi[sfds+i]), Memc[image], SZ_FNAME)
+ for (j = 0; !streq (SFD_IMAGE(Memi[sfds+j]), Memc[image]); j = j+1)
+ ;
+ if (j == i)
+ nimages = nimages + 1
+ }
+ call malloc (images, nimages, TY_POINTER)
+
+ nimages = 0
+ for (i = 0; i < nsfd; i = i + 1) {
+ call strcpy (SFD_IMAGE(Memi[sfds+i]), Memc[image], SZ_FNAME)
+ for (j = 0; j < nimages; j = j + 1)
+ if (streq (SFI_IMAGE(Memi[images+j]), Memc[image]))
+ break
+ if (j == nimages) {
+ k = 0
+ for (j = i; j < nsfd; j = j + 1)
+ if (streq (SFD_IMAGE(Memi[sfds+j]), Memc[image]))
+ k = k + 1
+ call malloc (sfi, SFI(k), TY_STRUCT)
+ call strcpy (Memc[image], SFI_IMAGE(sfi), SF_SZFNAME)
+ SFI_NSFD(sfi) = k
+ k = 0
+ for (j = i; j < nsfd; j = j + 1) {
+ sfd = Memi[sfds+j]
+ if (streq (SFD_IMAGE(sfd), Memc[image])) {
+ k = k + 1
+ SFD_SFI(sfd) = sfi
+ SFI_SFD(sfi,k) = sfd
+ }
+ }
+ Memi[images+nimages] = sfi
+ nimages = nimages + 1
+ }
+ }
+
+ SF_NSFD(sf) = nsfd
+ SF_SFDS(sf) = sfds
+ SF_NSTARS(sf) = nstars
+ SF_STARS(sf) = stars
+ SF_NFOCUS(sf) = nfocus
+ SF_FOCUS(sf) = focus
+ SF_NIMAGES(sf) = nimages
+ SF_IMAGES(sf) = images
+
+ # Find the average and best focus values. Sort the focus groups
+ # by magnitude and the star groups by focus.
+
+ call stf_fitfocus (sf)
+ do i = 1, SF_NFOCUS(sf) {
+ sff = SF_SFF(sf,i)
+ call qsort (SFF_SFD(sff,1), SFF_NSFD(sff), stf_magsort)
+ }
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ call qsort (SFS_SFD(sfs,1), SFS_NSFD(sfs), stf_focsort)
+ }
+
+ call sfree (sp)
+end
+
+
+# STF_LOG -- Print log of results
+
+procedure stf_log (sf, fd)
+
+pointer sf #I Main data structure
+int fd #I File descriptor
+
+int i, j, n
+pointer sp, str, sfd, sfs, sff, sfi
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Print banner and overall result.
+ call sysid (Memc[str], SZ_LINE)
+ call fprintf (fd, "%s\n\n")
+ call pargstr (Memc[str])
+
+ # Print each individual object organized by image.
+ call fprintf (fd, "%15.15s %7s %7s %7s")
+ call pargstr ("Image")
+ call pargstr ("Column")
+ call pargstr ("Line")
+ call pargstr ("Mag")
+ if (IS_INDEF(SF_F(sf))) {
+ call fprintf (fd, " %7s")
+ call pargstr (SF_WTYPE(sf))
+ } else {
+ call fprintf (fd, " %7s %7s")
+ call pargstr ("Focus")
+ call pargstr (SF_WTYPE(sf))
+ }
+ if (SF_WCODE(sf) == 4) {
+ call fprintf (fd, " %4s")
+ call pargstr ("Beta")
+ }
+ call fprintf (fd, " %7s %7s %3s\n")
+ call pargstr ("Ellip")
+ call pargstr ("PA")
+ call pargstr ("SAT")
+
+ do i = 1, SF_NIMAGES(sf) {
+ sfi = SF_SFI(sf,i)
+ n = 0
+ do j = 1, SFI_NSFD(sfi) {
+ sfd = SFI_SFD(sfi,j)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ if (n == 0) {
+ call fprintf (fd, "%15.15s")
+ call pargstr (SFD_IMAGE(sfd))
+ } else
+ call fprintf (fd, "%15w")
+ call fprintf (fd, " %7.2f %7.2f %7.2f")
+ call pargr (SFD_X(sfd))
+ call pargr (SFD_Y(sfd))
+ call pargr (-2.5*log10 (SFD_M(sfd) / SF_M(sf)))
+ if (IS_INDEF(SFD_F(sfd))) {
+ call fprintf (fd,
+ " %7.3f")
+ call pargr (SFD_W(sfd))
+ } else {
+ call fprintf (fd, " %7.6g %7.3f")
+ call pargr (SFD_F(sfd))
+ call pargr (SFD_W(sfd))
+ }
+ if (SF_WCODE(sf) == 4) {
+ call fprintf (fd, " %4.1f")
+ call pargr (SFD_BETA(sfd))
+ }
+ call fprintf (fd, " %7.2f %7d")
+ call pargr (SFD_E(sfd))
+ call pargr (SFD_PA(sfd))
+ if (SFD_NSAT(sfd) == 0)
+ call fprintf (fd, "\n")
+ else
+ call fprintf (fd, " *\n")
+ n = n + 1
+ }
+ }
+ if (n > 0)
+ call fprintf (fd, "\n")
+
+ # Print best values by star.
+ if (SF_NS(sf) > 1) {
+ n = 0
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ if (SFS_NF(sfs) > 1 || SFS_N(sfs) > 1) {
+ call stf_title (sf, NULL, sfs, NULL, Memc[str], SZ_LINE)
+ call fprintf (fd, " %s\n")
+ call pargstr (Memc[str])
+ n = n + 1
+ }
+ }
+ if (n > 0)
+ call fprintf (fd, "\n")
+ }
+
+ # Print averages at each focus.
+ if (SF_NF(sf) > 1) {
+ n = 0
+ do i = 1, SF_NFOCUS(sf) {
+ sff = SF_SFF(sf,i)
+ if (SFF_N(sff) > 1) {
+ call stf_title (sf, NULL, NULL, sff, Memc[str], SZ_LINE)
+ call fprintf (fd, " %s\n")
+ call pargstr (Memc[str])
+ n = n + 1
+ }
+ }
+ if (n > 0)
+ call fprintf (fd, "\n")
+ }
+
+ call stf_title (sf, NULL, NULL, NULL, Memc[str], SZ_LINE)
+ call fprintf (fd, "%s\n")
+ call pargstr (Memc[str])
+end
+
+
+# STF_TITLE -- Return result title string.
+# The title is dependent on whether an overall title, a title for a star
+# group, for a focus group, or for an indivdual object is desired.
+# The title also is adjusted for the select size type and the number
+# of objects in a group.
+
+procedure stf_title (sf, sfd, sfs, sff, title, sz_title)
+
+pointer sf #I Starfocus pointer
+pointer sfd #I Data pointer
+pointer sfs #I Star pointer
+pointer sff #I Focus pointer
+char title[sz_title] #O Title string
+int sz_title #I Size of title string
+
+pointer ptr
+int i, fd, stropen()
+errchk stropen
+
+begin
+ fd = stropen (title, sz_title, WRITE_ONLY)
+
+ if (sfd != NULL) {
+ call fprintf (fd, "%s @ (%.2f, %.2f):")
+ call pargstr (SFD_IMAGE(sfd))
+ call pargr (SFD_X(sfd))
+ call pargr (SFD_Y(sfd))
+ switch (SF_WCODE(sf)) {
+ case 4:
+ call fprintf (fd, " %s=%.2f (%3.1f), e=%.2f, pa=%d")
+ call pargstr (SF_WTYPE(sf))
+ call pargr (SFD_W(sfd))
+ call pargr (SFD_BETA(sfd))
+ call pargr (SFD_E(sfd))
+ call pargr (SFD_PA(sfd))
+ default:
+ call fprintf (fd, " %s=%.2f, e=%.2f, pa=%d")
+ call pargstr (SF_WTYPE(sf))
+ call pargr (SFD_W(sfd))
+ call pargr (SFD_E(sfd))
+ call pargr (SFD_PA(sfd))
+ }
+ if (SFD_SFS(sfd) != NULL) {
+ if (SFS_M(SFD_SFS(sfd)) != SF_M(sf)) {
+ call fprintf (fd, " , m=%.2f")
+ call pargr (-2.5*log10 (SFS_M(SFD_SFS(sfd)) / SF_M(sf)))
+ }
+ }
+ if (!IS_INDEF(SFD_F(sfd))) {
+ call fprintf (fd, ", f=%.4g")
+ call pargr (SFD_F(sfd))
+ }
+ } else if (sfs != NULL) {
+ ptr = SFS_SFD(sfs,1)
+ call fprintf (fd, "%s")
+ if (SFS_NF(sfs) > 1)
+ call pargstr ("Best focus estimate")
+ else if (SFS_N(sfs) > 1)
+ call pargstr ("Average star")
+ else {
+ for (i=1; SFD_STATUS(SFS_SFD(sfs,i))!=0; i=i+1)
+ ;
+ call pargstr (SFD_IMAGE(SFS_SFD(sfs,i)))
+ }
+ call fprintf (fd, " @ (%.2f, %.2f): %s=%.2f")
+ call pargr (SFD_X(ptr))
+ call pargr (SFD_Y(ptr))
+ call pargstr (SF_WTYPE(sf))
+ call pargr (SFS_W(sfs))
+ #if (SFS_M(sfs) != SF_M(sf)) {
+ call fprintf (fd, ", m=%.2f")
+ call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+ #}
+ if (!IS_INDEF(SFS_F(sfs))) {
+ call fprintf (fd, ", f=%.4g")
+ call pargr (SFS_F(sfs))
+ }
+ } else if (sff != NULL) {
+ if (SFF_NI(sff) == 1) {
+ for (i=1; SFD_STATUS(SFF_SFD(sff,i))!=0; i=i+1)
+ ;
+ call fprintf (fd, "%s")
+ call pargstr (SFD_IMAGE(SFF_SFD(sff,i)))
+ if (!IS_INDEF(SFF_F(sff))) {
+ call fprintf (fd, " at focus %.4g")
+ call pargr (SFF_F(sff))
+ }
+ call fprintf (fd, " with average")
+ } else {
+ if (IS_INDEF(SFF_F(sff)))
+ call fprintf (fd, "Average")
+ else {
+ call fprintf (fd, "Focus %.4g with average")
+ call pargr (SFF_F(sff))
+ }
+ }
+ call fprintf (fd, " %s of %.2f")
+ call pargstr (SF_WTYPE(sf))
+ call pargr (SFF_W(sff))
+ } else {
+ if (IS_INDEF(SF_F(sf))) {
+ if (SF_WCODE(sf) == 1) {
+ call fprintf (fd, " %s%d%% enclosed flux radius of ")
+ if (SF_N(sf) > 1)
+ call pargstr ("Average ")
+ else
+ call pargstr ("")
+ call pargr (100 * SF_LEVEL(sf))
+ } else {
+ if (SF_N(sf) > 1)
+ call fprintf (fd,
+ " Average full width at half maximum (%s) of ")
+ else
+ call fprintf (fd,
+ " Full width at half maximum (%s) of ")
+ call pargstr (SF_WTYPE(sf))
+ }
+ call fprintf (fd, "%.4f")
+ call pargr (SF_W(sf))
+ } else {
+ call fprintf (fd, " %s of %.6g with ")
+ if (SF_NS(sf) > 1) {
+ if (SF_NF(sf) > 1)
+ call pargstr ("Average best focus")
+ else
+ call pargstr ("Average focus")
+ } else {
+ if (SF_NF(sf) > 1)
+ call pargstr ("Best focus")
+ else
+ call pargstr ("Focus")
+ }
+ call pargr (SF_F(sf))
+ if (SF_WCODE(sf) == 1) {
+ call fprintf (fd, "%d%% enclosed flux radius of ")
+ call pargr (100 * SF_LEVEL(sf))
+ } else {
+ call fprintf (fd, "%s of ")
+ call pargstr (SF_WTYPE(sf))
+ }
+ call fprintf (fd, "%.2f")
+ call pargr (SF_W(sf))
+ }
+ }
+
+ call strclose (fd)
+end
+
+
+# STF_FITFOCUS -- Find the best focus.
+
+procedure stf_fitfocus (sf)
+
+pointer sf #I Starfocus pointer
+
+int i, j, k, n, jmin
+pointer x, y, sfd, sfs, sff, sfi
+real f, r, m, wr, wf
+bool fp_equalr()
+
+begin
+ # Set number of valid points, stars, focuses, images.
+ SF_N(sf) = 0
+ SF_YP1(sf) = 0
+ SF_YP2(sf) = 0
+ do i = 1, SF_NSFD(sf) {
+ sfd = SF_SFD(sf,i)
+ if (SFD_STATUS(sfd) == 0) {
+ SF_N(sf) = SF_N(sf) + 1
+ SF_YP1(sf) = min (SF_YP1(sf), SFD_YP1(sfd))
+ SF_YP2(sf) = max (SF_YP2(sf), SFD_YP2(sfd))
+ }
+ }
+ SF_NS(sf) = 0
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ SFS_N(sfs) = 0
+ SFS_M(sfs) = 0.
+ SFS_NF(sfs) = 0
+ do j = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs,j)
+ if (SFD_STATUS(SFS_SFD(sfs,j)) != 0)
+ next
+ SFS_N(sfs) = SFS_N(sfs) + 1
+ SFS_M(sfs) = SFS_M(sfs) + SFD_M(sfd)
+ sff = SFD_SFF(sfd)
+ for (k = 1; SFD_SFF(SFS_SFD(sfs,k)) != sff; k = k + 1)
+ ;
+ if (k == j)
+ SFS_NF(sfs) = SFS_NF(sfs) + 1
+ }
+ if (SFS_N(sfs) > 0) {
+ SFS_M(sfs) = SFS_M(sfs) / SFS_N(sfs)
+ SF_NS(sf) = SF_NS(sf) + 1
+ }
+ }
+ SF_NF(sf) = 0
+ do i = 1, SF_NFOCUS(sf) {
+ sff = SF_SFF(sf,i)
+ SFF_W(sff) = 0.
+ SFF_N(sff) = 0
+ SFF_NI(sff) = 0
+ wr = 0
+ do j = 1, SFF_NSFD(sff) {
+ sfd = SFF_SFD(sff,j)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ m = SFS_M(SFD_SFS(sfd))
+ wr = wr + m
+ SFF_W(sff) = SFF_W(sff) + m * SFD_W(sfd)
+ SFF_N(sff) = SFF_N(sff) + 1
+ sfi = SFD_SFI(sfd)
+ for (k = 1; SFD_SFI(SFF_SFD(sff,k)) != sfi; k = k + 1)
+ ;
+ if (k == j)
+ SFF_NI(sff) = SFF_NI(sff) + 1
+ }
+ if (SFF_N(sff) > 0) {
+ SFF_W(sff) = SFF_W(sff) / wr
+ SF_NF(sf) = SF_NF(sf) + 1
+ }
+ }
+ SF_NI(sf) = 0
+ do i = 1, SF_NIMAGES(sf) {
+ sfi = SF_SFI(sf,i)
+ SFI_N(sfi) = 0
+ do j = 1, SFI_NSFD(sfi)
+ if (SFD_STATUS(SFI_SFD(sfi,j)) == 0)
+ SFI_N(sfi) = SFI_N(sfi) + 1
+ if (SFI_N(sfi) > 0)
+ SF_NI(sf) = SF_NI(sf) + 1
+ }
+
+ # Find the average magnitude, best focus, and radius for each star.
+ # Find the brightest magnitude and average best focus and radius
+ # over all stars.
+
+ SF_BEST(sf) = SF_SFD(sf,1)
+ SF_F(sf) = 0.
+ SF_W(sf) = 0.
+ SF_M(sf) = 0.
+ SF_NS(sf) = 0
+ wr = 0.
+ wf = 0.
+ do i = 1, SF_NSTARS(sf) {
+ sfs = SF_SFS(sf,i)
+ call malloc (x, SFS_NSFD(sfs), TY_REAL)
+ call malloc (y, SFS_NSFD(sfs), TY_REAL)
+ k = 0
+ n = 0
+ do j = 1, SFS_NSFD(sfs) {
+ sfd = SFS_SFD(sfs,j)
+ if (SFD_STATUS(sfd) != 0)
+ next
+ r = SFD_W(sfd)
+ f = SFD_F(sfd)
+ if (!IS_INDEF(f))
+ k = k + 1
+ Memr[x+n] = f
+ Memr[y+n] = r
+ n = n + 1
+ if (r < SFD_W(SF_BEST(sf)))
+ SF_BEST(sf) = sfd
+ }
+
+ # Find the best focus and radius.
+ if (n == 0) {
+ SFS_F(sfs) = INDEF
+ SFS_W(sfs) = INDEF
+ SFS_M(sfs) = INDEF
+ SFS_N(sfs) = 0
+ } else if (k == 0) {
+ call alimr (Memr[y], n, f, r)
+ f = INDEF
+ m = SFS_M(sfs)
+ wr = wr + m
+ SFS_F(sfs) = f
+ SFS_W(sfs) = r
+ SFS_M(sfs) = m
+ SFS_N(sfs) = n
+ SF_W(sf) = SF_W(sf) + m * r
+ SF_M(sf) = max (SF_M(sf), m)
+ SF_NS(sf) = SF_NS(sf) + 1
+ } else {
+ SFS_N(sfs) = n
+ if (k < n) {
+ k = 0
+ do j = 0, n-1 {
+ if (!IS_INDEF(Memr[x+j])) {
+ Memr[x+k] = Memr[x+j]
+ Memr[y+k] = Memr[y+j]
+ k = k + 1
+ }
+ }
+ }
+ call xt_sort2 (Memr[x], Memr[y], k)
+ n = 0
+ do j = 1, k-1 {
+ if (fp_equalr (Memr[x+j], Memr[x+n])) {
+ if (Memr[y+j] < Memr[y+n])
+ Memr[y+n] = Memr[y+j]
+ } else {
+ n = n + 1
+ Memr[x+n] = Memr[x+j]
+ Memr[y+n] = Memr[y+j]
+ }
+ }
+ n = n + 1
+
+ # Find the minimum radius
+ jmin = 0
+ do j = 0, n-1
+ if (Memr[y+j] < Memr[y+jmin])
+ jmin = j
+
+ # Use parabolic interpolation to find the best focus
+ if (jmin == 0 || jmin == n-1) {
+ f = Memr[x+jmin]
+ r = Memr[y+jmin]
+ } else
+ call stf_parab (Memr[x+jmin-1], Memr[y+jmin-1], f, r)
+
+ m = SFS_M(sfs)
+ wr = wr + m
+ wf = wf + m
+ SFS_F(sfs) = f
+ SFS_W(sfs) = r
+ SFS_M(sfs) = m
+ SF_F(sf) = SF_F(sf) + m * f
+ SF_W(sf) = SF_W(sf) + m * r
+ SF_M(sf) = max (SF_M(sf), m)
+ SF_NS(sf) = SF_NS(sf) + 1
+ }
+ call mfree (x, TY_REAL)
+ call mfree (y, TY_REAL)
+ }
+
+ if (wr > 0.)
+ SF_W(sf) = SF_W(sf) / wr
+ else {
+ SF_W(sf) = INDEF
+ SF_M(sf) = INDEF
+ }
+ if (wf > 0.)
+ SF_F(sf) = SF_F(sf) / wf
+ else
+ SF_F(sf) = INDEF
+end
+
+
+# STF_PARAB -- Find the minimum of a parabolic fit to three points.
+
+procedure stf_parab (x, y, xmin, ymin)
+
+real x[3]
+real y[3]
+real xmin
+real ymin
+
+double x12, x13, x23, x213, x223, y13, y23, a, b, c
+
+begin
+ x12 = x[1] - x[2]
+ x13 = x[1] - x[3]
+ x23 = x[2] - x[3]
+ x213 = x13 * x13
+ x223 = x23 * x23
+ y13 = y[1] - y[3]
+ y23 = y[2] - y[3]
+ c = (y13 - y23 * x13 / x23) / (x213 - x223 * x13 / x23)
+ b = (y23 - c * x223) / x23
+ a = y[3]
+ xmin = -b / (2 * c)
+ ymin = a + b * xmin + c * xmin * xmin
+ xmin = xmin + x[3]
+end
+
+
+# STF_MAGSORT -- Compare two star structures by average magnitude.
+
+int procedure stf_magsort (sfd1, sfd2)
+
+pointer sfd1, sfd2 # Structures to compare
+pointer sfs1, sfs2 # Star structures for magnitudes
+
+begin
+ sfs1 = SFD_SFS(sfd1)
+ sfs2 = SFD_SFS(sfd2)
+ if (SFS_M(sfs1) > SFS_M(sfs2))
+ return (-1)
+ else if (SFS_M(sfs1) < SFS_M(sfs2))
+ return (1)
+ else
+ return (0)
+end
+
+
+# STF_FOCSORT -- Compare two star structures by focus.
+
+int procedure stf_focsort (sfd1, sfd2)
+
+pointer sfd1, sfd2 # Structures to compare
+
+begin
+ if (SFD_F(sfd1) < SFD_F(sfd2))
+ return (-1)
+ else if (SFD_F(sfd1) > SFD_F(sfd2))
+ return (1)
+ else
+ return (0)
+end
+
+
+# STF_DISPLAY -- Display image if necessary.
+# The user is required to display the first image separately.
+
+procedure stf_display (image, frame)
+
+char image[ARB] #I Image to display
+int frame #I Display frame to use
+
+int i, status
+pointer sp, dname, ds, iw, imd_mapframe(), iw_open()
+bool xt_imnameeq()
+errchk clcmdw
+
+begin
+ call smark (sp)
+ call salloc (dname, SZ_LINE, TY_CHAR)
+
+ ds = imd_mapframe (1, READ_WRITE, NO)
+ do i = 1, MAX_FRAMES {
+ iferr (iw = iw_open (ds, i, Memc[dname], SZ_LINE, status))
+ next
+ call iw_close (iw)
+ if (xt_imnameeq (image, Memc[dname]))
+ break
+ }
+ call imunmap (ds)
+
+ if (!xt_imnameeq (image, Memc[dname])) {
+ call sprintf (Memc[dname], SZ_LINE, "display %s frame=%d fill=yes")
+ call pargstr (image)
+ call pargi (frame)
+ call clcmdw (Memc[dname])
+ }
+
+ call sfree (sp)
+end
+
+
+# STFCUR -- Debugging routine.
+# Replace calls to clgcur with stfcur so that an text file containing the
+# cursor coordinates may be specified when running standalone (such as
+# under a debugger).
+
+int procedure stfcur (cur, wx, wy, wcs, key, cmd, sz_cmd)
+
+char cur[ARB] # Cursor name
+real wx, wy # Cursor coordinate
+int wcs # WCS
+int key # Key
+char cmd[sz_cmd] # Command
+int sz_cmd # Size of command
+
+int fd, stat, open(), fscan()
+pointer fname
+errchk open
+
+begin
+ if (fd == NULL) {
+ call malloc (fname, SZ_FNAME, TY_CHAR)
+ call clgstr (cur, Memc[fname], SZ_FNAME)
+ fd = open (Memc[fname], READ_ONLY, TEXT_FILE)
+ call mfree (fname, TY_CHAR)
+ }
+
+ stat = fscan (fd)
+ if (stat == EOF) {
+ call close (fd)
+ return (stat)
+ }
+
+ call gargr (wx)
+ call gargr (wy)
+ call gargi (wcs)
+ call gargi (key)
+
+ return (stat)
+end
diff --git a/noao/obsutil/src/starfocus/x_starfocus.x b/noao/obsutil/src/starfocus/x_starfocus.x
new file mode 100644
index 00000000..a3642bae
--- /dev/null
+++ b/noao/obsutil/src/starfocus/x_starfocus.x
@@ -0,0 +1,2 @@
+task psfmeasure = t_psfmeasure,
+ starfocus = t_starfocus
generated by cgit (git 2.34.1) at 2025-09-20 11:27:34 -0400