diff options
| author | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
|---|---|---|
| committer | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
| commit | fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch) | |
| tree | bdda434976bc09c864f2e4fa6f16ba1952b1e555 /pkg/utilities/nttools/stxtools/wcslab | |
| download | iraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz | |
Initial commit
Diffstat (limited to 'pkg/utilities/nttools/stxtools/wcslab')
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/mkpkg | 17 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/mkpkg.ori | 18 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/psiescape.h | 80 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/t_wcslab.x | 136 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wcs_desc.h | 219 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wcslab.h | 98 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wcslab.x | 935 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wlgrid.x | 448 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wllabel.x | 1100 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wllabel.x.ori | 1077 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wlsetup.x | 1000 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wlutil.x | 390 | ||||
| -rw-r--r-- | pkg/utilities/nttools/stxtools/wcslab/wlwcslab.x | 181 |
13 files changed, 5699 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/stxtools/wcslab/mkpkg b/pkg/utilities/nttools/stxtools/wcslab/mkpkg new file mode 100644 index 00000000..f5083925 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/mkpkg @@ -0,0 +1,17 @@ +# WCSLAB + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + wlutil.x <imio.h> <imhdr.h> <gset.h> <math.h> + wcslab.x <gset.h> <imhdr.h> <mwset.h> <math.h> "wcslab.h"\ + "wcs_desc.h" + wlwcslab.x <gio.h> <gset.h> "wcslab.h" "wcs_desc.h" + wlsetup.x <gset.h> <mach.h> <math.h> <math/curfit.h>\ + "wcslab.h" "wcs_desc.h" + wlgrid.x <gset.h> <math.h> "wcslab.h" "wcs_desc.h" + wllabel.x <gset.h> <math.h> "psiescape.h" "wcslab.h" "wcs_desc.h" + ; diff --git a/pkg/utilities/nttools/stxtools/wcslab/mkpkg.ori b/pkg/utilities/nttools/stxtools/wcslab/mkpkg.ori new file mode 100644 index 00000000..7108366c --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/mkpkg.ori @@ -0,0 +1,18 @@ +# WCSLAB + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + t_wcslab.x <gset.h> <imhdr.h> + wlutil.x <imio.h> <imhdr.h> <gset.h> <math.h> + wcslab.x <gset.h> <imhdr.h> <mwset.h> <math.h> "wcslab.h"\ + "wcs_desc.h" + wlwcslab.x <gio.h> <gset.h> "wcslab.h" "wcs_desc.h" + wlsetup.x <gset.h> <mach.h> <math.h> <math/curfit.h>\ + "wcslab.h" "wcs_desc.h" + wlgrid.x <gset.h> <math.h> "wcslab.h" "wcs_desc.h" + wllabel.x <gset.h> <math.h> "wcslab.h" "wcs_desc.h" + ; diff --git a/pkg/utilities/nttools/stxtools/wcslab/psiescape.h b/pkg/utilities/nttools/stxtools/wcslab/psiescape.h new file mode 100644 index 00000000..46da4ea2 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/psiescape.h @@ -0,0 +1,80 @@ +.help psiescape.h 1May92 plot +.ih +.NAME +psiescape.h -- Define the special GIO escape instructions. +.ih +DESCRIPTION +The following escape instructions are defined for the PostScript GKI +Interpreter: + + PS_CODE - Send the raw PostScript code to the output. + PS_IMAGE_RED_LUT - Download a new lookup table for the Red component + of an image. + PS_IMAGE_GREEN_LUT - Download a new lookup table for the Green component + of an image. + PS_IMAGE_BLUE_LUT - Download a new lookup table for the Blue component + of an image. + PS_GR_RED_LUT - Download a new lookup table for the Red component + for the graphics. + PS_GR_GREEN_LUT - Download a new lookup table for the Green component + for the graphics. + PS_GR_BLUE_LUT - Download a new lookup table for the Blue component + for the graphics. + PS_ROMAN_FONT - Specify a new font for the normal font. + PS_GREEK_FONT - Specify a new font for the greek font. + PS_ITALIC_FONT - Specify a new font for the italic font. + PS_BOLD_FONT - Specify a new font for the bold font. + PS_VARIABLE_SPACE - Change whether characters are mono-spaced or + variable-spaced. + PS_DASH, PS_DOT, + PS_SPACE - Change the sizes of a dash, a dot, and the space + between them. + PS_FILL_PATTERN - Add/change fill patterns. + +The size of the instruction array should have the following minimums: + - For the PS_CODE instruction, the array should be the length of the string + being passed. + - The size of each image LUT array is PS_IMAGE_LUT_SIZE. + - The size of each graphics LUT array is PS_GR_LUT_SIZE. + - The font arrays should be the length of the string containing the name + of the font to use. + - For the PS_VARIABLE_SPACE, the size is PS_VARIABLE_SPACE_SIZE. + - For PS_FILL_PATTERN, the size is PS_FILL_SIZE. + +.ih +SEE ALSO +t_psikern +.endhelp +#--------------------------------------------------------------------------- + +# Define the escape instructions. +define PS_CODE 1001 +define PS_IMAGE_RED_LUT 1002 +define PS_IMAGE_GREEN_LUT 1003 +define PS_IMAGE_BLUE_LUT 1004 +define PS_GR_RED_LUT 1005 +define PS_GR_GREEN_LUT 1006 +define PS_GR_BLUE_LUT 1007 +define PS_ROMAN_FONT 1008 +define PS_GREEK_FONT 1009 +define PS_ITALIC_FONT 1010 +define PS_BOLD_FONT 1011 +define PS_VARIABLE_SPACE 1012 +define PS_DOT 1013 +define PS_DASH 1014 +define PS_SPACE 1015 +define PS_FILL_PATTERN 1016 +define PS_IMAGE_LUT 1017 +define PS_GRAPHICS_LUT 1018 + +# Define the sizes of the instruction arrays. +define PS_MAX_LUT_VALUE 255 +define PS_IMAGE_LUT_SIZE 256 +define PS_GR_LUT_SIZE 16 +define PS_VARIABLE_SPACE_SIZE 1 +define PS_FILL_SIZE 9 + +# Define how to pack/unpack a LUT defined as reals from 0-1 into +# a short array from 0-PS_MAX_LUT_VALUE. +define PS_PACKLUT (int($1*PS_MAX_LUT_VALUE)) +define PS_UNPACKLUT ($1/real(PS_MAX_LUT_VALUE)) diff --git a/pkg/utilities/nttools/stxtools/wcslab/t_wcslab.x b/pkg/utilities/nttools/stxtools/wcslab/t_wcslab.x new file mode 100644 index 00000000..acafdf2f --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/t_wcslab.x @@ -0,0 +1,136 @@ +include <gset.h> +include <imhdr.h> + +# T_WCSLAB -- Procedure to draw labels and grids in sky projection coordinates. +# +# Description +# T_wcslab produces a labelling and grid based on the MWCS of a +# specified image. This is the task interface to the programmer interface +# wcslab. See wcslab.x for more information. +# +# Bugs +# Can only handle sky projections for Right Ascension/Declination. This +# should be able to deal with any of the projections for this system, but +# has only been tested with the Tangent projection. +# + +procedure t_wcslab() + +pointer image # I: name of the image +int frame # I: display frame containing the image +bool do_fill # I: true if the graph fills the specified viewport +int mode # I: the graphics stream mode +pointer device # I: the name of the graphics device +real vl, vr, vb, vt # I: the edges of the graphics viewport + +pointer sp, title, gp, im, mw +real c1, c2, l1, l2 +bool clgetb() +int clgeti(), strncmp() +pointer gopen(), immap(), mw_openim() +real clgetr() + +begin + # Get memory. + call smark (sp) + call salloc (device, SZ_FNAME, TY_CHAR) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (title, SZ_LINE, TY_CHAR) + + # Since all the MWCS information comes from an image open it. + call clgstr ("image", Memc[image], SZ_FNAME) + + if (Memc[image] != EOS) { + + # Open the image. + im = immap (Memc[image], READ_ONLY, 0) + + # Quit if the image is not 2-dimensional. + if (IM_NDIM(im) != 2) { + call eprintf ("Image: %s is not 2-dimensional\n") + call pargstr (Memc[image]) + call sfree (sp) + return + } + + # Set the default input image column and line limits. + c1 = 1.0 + c2 = real (IM_LEN(im,1)) + l1 = 1.0 + l2 = real (IM_LEN(im,2)) + + # Open the WCS. + mw = mw_openim (im) + + # Set up the default image title. + call strcpy (Memc[image], Memc[title], SZ_LINE) + call strcat (": ", Memc[title], SZ_LINE) + call strcat (IM_TITLE(im), Memc[title], SZ_LINE) + + } else { + + # Set the image information to undefined. All this will + # be determined in wcslab. + Memc[title] = EOS + im = NULL + mw = NULL + c1 = 0.0 + c2 = 1.0 + l1 = 0.0 + l2 = 1.0 + } + + # Set the graphics mode depending on whether we are appending to a plot + # or starting a new plot. + do_fill = clgetb ("fill") + if (clgetb ("append")) + mode = APPEND + else + mode = NEW_FILE + + # Open graphics. + call clgstr ("device", Memc[device], SZ_FNAME) + + # If we are appending, get the previous viewing parameters. + if (mode == APPEND) { + + gp = gopen (Memc[device], APPEND, STDGRAPH) + call ggview (gp, vl, vr, vb, vt) + do_fill = true + + # If drawing on the image display device try to match viewports. + } else if (strncmp (Memc[device], "imd", 3) == 0) { + + frame = clgeti ("frame") + vl = clgetr ("vl") + vr = clgetr ("vr") + vb = clgetr ("vb") + vt = clgetr ("vt") + call wl_imd_viewport (frame, im, c1, c2, l1, l2, vl, vr, vb, vt) + gp = gopen (Memc[device], NEW_FILE, STDGRAPH) + + # Otherwise set up a standard viewport. + } else { + vl = clgetr ("vl") + vr = clgetr ("vr") + vb = clgetr ("vb") + vt = clgetr ("vt") + gp = gopen (Memc[device], NEW_FILE, STDGRAPH) + } + + # Set the viewport. + call gseti (gp, G_WCS, 1) + call wl_map_viewport (gp, c1, c2, l1, l2, vl, vr, vb, vt, do_fill) + + # All reading from CL parameters is now done. Everything necessary to + # do the plotting is in the WCSLAB descriptor. Do it. + call wcslab (mw, c1, c2, l1, l2, gp, Memc[title]) + + # Release the memory. + call gclose (gp) + if (mw != NULL) + call mw_close (mw) + if (im != NULL) + call imunmap (im) + call sfree (sp) +end diff --git a/pkg/utilities/nttools/stxtools/wcslab/wcs_desc.h b/pkg/utilities/nttools/stxtools/wcslab/wcs_desc.h new file mode 100644 index 00000000..4f6b2a30 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wcs_desc.h @@ -0,0 +1,219 @@ +# WCS_DESC - The definition of the WCSLAB descriptor memory structure. +# +# Description +# This include file defines the memory structures and macros needed to +# access elements of a WCSLAB descriptor. The descriptor provides all +# the necessary elements for the routine wcslab to produce a labeled +# graph. +# +# History +# 9May91 - Created the descriptor. Jonathan D. Eisenhamer, STScI. +# 15May91 - Modified the descriptor to contain only pointers to arrays. +# Two routines, wcs_create and wcs_destroy are required to +# create the arrays that are pointed to in the descriptor. +# Also seperated the include file from the wcslab.h file. jde +# 12Jun91 - Rewrote some of the labelling parameters. jde +# 20Jun91 - Redesigned much of the parameters. jde +#--------------------------------------------------------------------------- + +# Value of opposite axis that polar labels should appear along. +define WL_POLAR_LABEL_POSITION Memd[P2D($1)] + +# The rotation between the Logical and World coordinate systems. +define WL_ROTA Memd[P2D($1+2)] + +# Size of the axis titles. +define WL_AXIS_TITLE_SIZE Memr[P2R($1+4)] + +# The offset required to properly calculate positions in the image display. +define WL_IMAGE_X_OFF Memr[P2R($1+5)] +define WL_IMAGE_Y_OFF Memr[P2R($1+6)] + +# Size of the grid labels. +define WL_LABEL_SIZE Memr[P2R($1+7)] + +# Major tick mark size. +define WL_MAJ_TICK_SIZE Memr[P2R($1+8)] + +# Minor tick mark size. +define WL_MIN_TICK_SIZE Memr[P2R($1+9)] + +# Magnification of the text size for the title. +define WL_TITLE_SIZE Memr[P2R($1+10)] + +# The side in polar/near-polar plots not to put Axis 1 labels. +define WL_BAD_LABEL_SIDE Memi[$1+11] + +# The type of graph that will be produced. The possible value are: +# +# UNKNOWN -> Graph type will be determined +# NORMAL -> Approximate a cartesian grid +# POLAR -> Graph center on a pole +# NEAR_POLAR -> Graph very close to a pole + +define WL_GRAPH_TYPE Memi[$1+12] + +# Number of segments each line should be broken into to plot it. +define WL_LINE_SEGMENTS Memi[$1+13] + +# The grid line type for major grids. The possible values are to standard +# IRAF GIO polyline types. +define WL_MAJ_LINE_TYPE Memi[$1+14] + +# The grid line type for minor grids. The possible values are to standard +# IRAF GIO polyline types. +define WL_MIN_LINE_TYPE Memi[$1+15] + +# The number of label points. +define WL_N_LABELS Memi[$1+16] + +# The graphic WCS that is set to NDC units. +define WL_NDC_WCS Memi[$1+17] + +# The graphic WCS used to plot the grid lines. +define WL_PLOT_WCS Memi[$1+18] + +# The direction of the latitude labelling on polar graphs. Possible values are: +# +# BOTTOM -> Towards the bottom of the graph. +# TOP -> Towards the top of the graph. +# RIGHT -> Towards the right of the graph. +# LEFT -> Towards the left of the graph. + +define WL_POLAR_LABEL_DIRECTION Memi[$1+19] + +# The possible axis types. The possible values are: +# +# RA_DEC_TAN - The tangential display in right ascension and declination. +# LINEAR - General linear systems. + +define WL_SYSTEM_TYPE Memi[$1+20] + +# Define which side of the graph will have the title. +define WL_TITLE_SIDE Memi[$1+21] + +# True if the axis mapping has reversed the order of the axis relative +# to the logical system. +define WL_AXIS_FLIP Memi[$1+22] + +# TRUE if the labels should always be printed in full form. +define WL_ALWAYS_FULL_LABEL Memi[$1+23] + +# TRUE if the grid labels should rotate with the grid lines. +define WL_LABEL_ROTATE Memi[$1+26] + +# True if coordinate labels are to be written. +define WL_LABON Memi[$1+27] + +# True if we are to write labels outside the window borders. Else, write +# them inside. +define WL_LABOUT Memi[$1+28] + +# True if we are to draw the major grid lines. +define WL_MAJ_GRIDON Memi[$1+29] + +# True if we are to draw the minor grid lines. +define WL_MIN_GRIDON Memi[$1+30] + +# True if the graph parameters should be written back out to the +# parameter file. +define WL_REMEMBER Memi[$1+31] + +# TRUE if tick marks should point into the graph. +define WL_TICK_IN Memi[$1+32] + +# Titles to label each axis. +define WL_AXIS_TITLE_PTR Memi[$1+33] +define WL_AXIS_TITLE Memc[WL_AXIS_TITLE_PTR($1)+(($2-1)*SZ_LINE)] + +# The sides the axis titles will appear. +define WL_AXIS_TITLE_SIDE_PTR Memi[$1+34] +define WL_AXIS_TITLE_SIDE Memi[WL_AXIS_TITLE_SIDE_PTR($1)+$2-1] + +# Beginning values to start labeling the axes. +define WL_BEGIN_PTR Memi[$1+35] +define WL_BEGIN Memd[WL_BEGIN_PTR($1)+$2-1] + +# The name of the graphics device. +#define WL_DEVICE_PTR Memi[$1+36] +#define WL_DEVICE Memc[WL_DEVICE_PTR($1)] + +# Value to stop labeling the axes. +define WL_END_PTR Memi[$1+37] +define WL_END Memd[WL_END_PTR($1)+$2-1] + +# The graphics descriptor. +define WL_GP Memi[$1+38] + +# The angle of text at this label point. +define WL_LABEL_ANGLE_PTR Memi[$1+40] +define WL_LABEL_ANGLE Memd[WL_LABEL_ANGLE_PTR($1)+$2-1] + +# Which axis the label represents. +define WL_LABEL_AXIS_PTR Memi[$1+41] +define WL_LABEL_AXIS Memi[WL_LABEL_AXIS_PTR($1)+$2-1] + +# The positions of tick mark/grid labels. +define WL_LABEL_POSITION_PTR Memi[$1+42] +define WL_LABEL_POSITION Memd[WL_LABEL_POSITION_PTR($1)+$2-1+(($3-1)*MAX_LABEL_POINTS)] +# +# NOTE: If the axis are transposed, the positions represented here are +# the corrected, transposed values. + +# The sides the labels for each axis should appear on. +define WL_LABEL_SIDE_PTR Memi[$1+43] +define WL_LABEL_SIDE Memb[WL_LABEL_SIDE_PTR($1)+$2-1+(($3-1)*N_SIDES)] + +# The value of the label. +define WL_LABEL_VALUE_PTR Memi[$1+44] +define WL_LABEL_VALUE Memd[WL_LABEL_VALUE_PTR($1)+$2-1] + +# The center of the transformations in the logical system. +define WL_LOGICAL_CENTER_PTR Memi[$1+45] +define WL_LOGICAL_CENTER Memd[WL_LOGICAL_CENTER_PTR($1)+$2-1] + +# The coordinate transformation from Logical to World. +define WL_LWCT Memi[$1+46] + +# Major grid intervals for the axis. +define WL_MAJ_I_PTR Memi[$1+47] +define WL_MAJOR_INTERVAL Memd[WL_MAJ_I_PTR($1)+$2-1] + +# The minor intervals for the axis. +define WL_MIN_I_PTR Memi[$1+48] +define WL_MINOR_INTERVAL Memi[WL_MIN_I_PTR($1)+$2-1] + +# Remember the extent of the labels around the plot box. +define WL_NV_PTR Memi[$1+49] +define WL_NEW_VIEW Memr[WL_NV_PTR($1)+$2-1] + +# The MWL structure. +define WL_MW Memi[$1+50] + +# The values of the sides of the screen. The indexes are defined as follows: +# +# TOP -> Y-axis value at the top of display. +# BOTTOM -> Y-axis value at bottom of display +# RIGHT -> X-axis value at right of display. +# LEFT -> X-axis value at left of display. +# +define WL_SCREEN_BOUNDARY_PTR Memi[$1+51] +define WL_SCREEN_BOUNDARY Memd[WL_SCREEN_BOUNDARY_PTR($1)+$2-1] + +# The title that will be placed on the plot. +define WL_TITLE_PTR Memi[$1+52] +define WL_TITLE Memc[WL_TITLE_PTR($1)] + +# The coordinate transformation from World to Logical. +define WL_WLCT Memi[$1+53] + +# The center of the transformations in the world system. +define WL_WORLD_CENTER_PTR Memi[$1+54] +define WL_WORLD_CENTER Memd[WL_WORLD_CENTER_PTR($1)+$2-1] + +# The length of this structure. +define WL_LEN 55+1 + +#--------------------------------------------------------------------------- +# End of wcs_desc +#--------------------------------------------------------------------------- diff --git a/pkg/utilities/nttools/stxtools/wcslab/wcslab.h b/pkg/utilities/nttools/stxtools/wcslab/wcslab.h new file mode 100644 index 00000000..cf088323 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wcslab.h @@ -0,0 +1,98 @@ +# Definitions file for WCSLAB + +# Define various important dimensions + +define MAX_DIM 10 # Maximum number of dimensions +define N_DIM 2 # Dimensionality of plotting space +define N_SIDES 4 # Number of sides to a window +define MAX_LABEL_POINTS 100 # The maximum number of possible label points +define N_EDGES 20 # Number of edges being examined from the window + +# Define the types of graphs possible. + +define GRAPHTYPES "|normal|polar|near_polar|" +define NORMAL 1 +define POLAR 2 +define NEAR_POLAR 3 + +# Define the graph sides. The ordering matches the calls to the GIO package. + +define GRAPHSIDES "|left|right|bottom|top|" +define LEFT 1 +define RIGHT 2 +define BOTTOM 3 +define TOP 4 + +# Define which index refers to the X-axis and which refers to the Y-axis. + +define X_DIM 1 +define Y_DIM 2 +define AXIS1 1 +define AXIS2 2 + +# Define which axis is longitude and which axis is latitude. + +define LONGITUDE 1 +define LATITUDE 2 + +# Define the available precisions for labelling + +define HOUR 1 +define DEGREE 1 +define MINUTE 2 +define SECOND 3 +define SUBSEC_LOW 4 +define SUBSEC_HIGH 5 + +# Define the possible MWCS transformation types. + +define RA_DEC_DICTIONARY "|tan|arc|sin|" +define LINEAR_DICTIONARY "|linear|" + +define NUMBER_OF_SUPPORTED_TYPES 2 +define RA_DEC 1 +define LINEAR 2 + +define AXIS 3B # transform all axes in any MWCS call + +# Some useful graphics definitions and defaults + +define NDC_WCS 0 # the base graphics WCS +define POLE_MARK_SHAPE 4 # the pole mark is a cross +define POLE_MARK_SIZE 3.0 # the half-size of the cross +define DISTANCE_TO_POLE 0.1 # % distance to pole for lines of longitude +define LINE_SIZE 1. # line width for lines and ticks +define MIN_ANGLE 10. # minimum angle for text rotation +define BOTTOM_LEFT .1 # default bottom left corner of viewport +define TOP_RIGHT .9 # default top right corner of viewport + +# Units conversion macros + +define RADTOST (240*RADTODEG($1)) # Radians to seconds of time +define RADTOSA (3600*RADTODEG($1)) # Radians to seconds of arc +define STTORAD (DEGTORAD(($1)/240)) # Seconds of time to radians +define SATORAD (DEGTORAD(($1)/3600)) # Seconds of arc to radians +define RADTOHRS (RADTODEG(($1)/15)) # Radians to hours +define HRSTORAD (DEGTORAD(15*($1))) # Hours to radians +define DEGTOST (240*($1)) # Degrees to seconds of time. +define STTODEG (($1)/240) # Seconds of time to degrees. +define DEGTOSA (3600*($1)) # Degrees to seconds of arc. +define SATODEG (($1)/3600) # Seconds of arc to degrees. +define HRSTODEG (($1)*15) # Hours to degrees. +define DEGTOHRS (($1)/15) # Degrees to hours. +define STPERDAY 86400 # Seconds per day + +# Other useful macros + +define INVERT ($1 < 45 || $1 > 315 || ( $1 > 135 && $1 < 225 )) + +# Define the latitudes of the north and south poles + +define NORTH_POLE_LATITUDE 90.0D0 +define SOUTH_POLE_LATITUDE -90.0D0 + +# Define sections of a circle + +define QUARTER_CIRCLE 90.0D0 +define HALF_CIRCLE 180.0D0 +define FULL_CIRCLE 360.0D0 diff --git a/pkg/utilities/nttools/stxtools/wcslab/wcslab.x b/pkg/utilities/nttools/stxtools/wcslab/wcslab.x new file mode 100644 index 00000000..2e6974c6 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wcslab.x @@ -0,0 +1,935 @@ +include <gset.h> +include <imhdr.h> +include <math.h> +include <mwset.h> +include "wcslab.h" +include "wcs_desc.h" +include <ctype.h> + + +# WCSLAB -- Procedure to draw labels and grids in sky projection coordinates. +# +# Description +# Wcslab produces a labelling and grid based on the MWCS of a +# specified image. +# +# The only things necessary to run this routine are: +# 1) Open an image and pass the image descriptor in im. +# 2) Open the graphics device and set the desired viewport (with a +# gsview call). +# 3) Make sure that the wlpars pset is available. +# +# Upon return, the graphics system will be in the state that it had been +# left in and a "virtual viewport" will be returned in the arguments +# left, right, bottom, top. This viewport defines the region where labels +# and/or titles were written. If any graphics is performed within this +# region, chances are that something will be overwritten. If any other +# graphics remain outside this region, then what was produced by this +# subroutine will remain untouched. +# +# Bugs +# Can only handle sky projections for Right Ascension/Declination. This +# should be able to deal with any of the projections for this system, but +# has only been tested with the Tangent projection. + +procedure wcslab (mw, log_x1, log_x2, log_y1, log_y2, gp, title) + +pointer mw # I: the wcs descriptor +real log_x1, log_x2 # I/O: the viewport +real log_y1, log_y2 # I/O: the viewport +pointer gp # I: the graphics descriptor +char title[ARB] # I: the image title + +pointer wd +real junkx1, junkx2, junky1, junky2 +bool clgetb() +pointer wl_create() +errchk clgstr + +begin + # Allocate the descriptor. + wd = wl_create() + + # Set the title name. + call strcpy (title, WL_TITLE(wd), SZ_LINE) + + # Set the WCS descriptor. If the descriptor is NULL or if + # the use_wcs parameter is yes, retrieve the parameter + # specified wcs. + if (mw == NULL) + call wl_wcs_params (mw, log_x1, log_x2, log_y1, log_y2) + else if (clgetb ("usewcs")) + call wl_wcs_params (mw, junkx1, junkx2, junky1, junky2) + WL_MW(wd) = mw + + # Determine axis types. + call wl_get_system_type (WL_MW(wd), WL_SYSTEM_TYPE(wd), + WL_LOGICAL_CENTER(wd,1), WL_WORLD_CENTER(wd,1), WL_AXIS_FLIP(wd)) + if (IS_INDEFI(WL_SYSTEM_TYPE(wd))) + call error (0, "WCSLAB: Image WCS is unsupported\n") + + # Get the parameters. + call wl_gr_inparams (wd) + + # Copy the graphics descriptor. + WL_GP(wd) = gp + + # Set the plot window in pixels (the logical space of the WCS). + WL_SCREEN_BOUNDARY(wd,LEFT) = log_x1 + WL_SCREEN_BOUNDARY(wd,RIGHT) = log_x2 + WL_SCREEN_BOUNDARY(wd,BOTTOM) = log_y1 + WL_SCREEN_BOUNDARY(wd,TOP) = log_y2 + + # Plot and label the coordinate grid. + call wl_wcslab (wd) + + # Return the possibly modified graphics descriptor and viewport. + gp = WL_GP(wd) + call gsview (gp, WL_NEW_VIEW(wd,LEFT), WL_NEW_VIEW(wd,RIGHT), + WL_NEW_VIEW(wd,BOTTOM), WL_NEW_VIEW(wd,TOP)) + + # Save the current parameters. + if (WL_REMEMBER(wd) == YES) + call wl_gr_remparams (wd) + + # Release the memory. + call wl_destroy (wd) +end + + +# WL_CREATE -- Create a WCSLAB descriptor and initialize it. +# +# Description +# This routine allocates the memory for the WCSLAB descriptor and +# subarrays and initializes values. +# +# Returns +# the pointer to the WCSLAB descriptor. + +pointer procedure wl_create() + +int i,j +pointer wd + +begin + # Allocate the descriptor memory. + call malloc (wd, WL_LEN, TY_STRUCT) + + # Allocate the subarrays. + call malloc (WL_AXIS_TITLE_PTR(wd), SZ_LINE * N_DIM, TY_CHAR) + call malloc (WL_AXIS_TITLE_SIDE_PTR(wd), N_SIDES * N_DIM, TY_BOOL) + call malloc (WL_BEGIN_PTR(wd), N_DIM, TY_DOUBLE) + call malloc (WL_END_PTR(wd), N_DIM, TY_DOUBLE) + call malloc (WL_LABEL_ANGLE_PTR(wd), MAX_LABEL_POINTS, TY_DOUBLE) + call malloc (WL_LABEL_AXIS_PTR(wd), MAX_LABEL_POINTS, TY_INT) + call malloc (WL_LABEL_POSITION_PTR(wd), N_DIM * MAX_LABEL_POINTS, + TY_DOUBLE) + call malloc (WL_LABEL_SIDE_PTR(wd), N_DIM * N_SIDES, TY_BOOL) + call malloc (WL_LABEL_VALUE_PTR(wd), MAX_LABEL_POINTS, TY_DOUBLE) + call malloc (WL_LOGICAL_CENTER_PTR(wd), N_DIM, TY_DOUBLE) + call malloc (WL_MAJ_I_PTR(wd), N_DIM, TY_DOUBLE) + call malloc (WL_MIN_I_PTR(wd), N_DIM, TY_INT) + call malloc (WL_NV_PTR(wd), N_SIDES, TY_REAL) + call malloc (WL_SCREEN_BOUNDARY_PTR(wd), N_SIDES, TY_DOUBLE) + call malloc (WL_TITLE_PTR(wd), SZ_LINE, TY_CHAR) + call malloc (WL_WORLD_CENTER_PTR(wd), N_DIM, TY_DOUBLE) + + # Initialize the simple values (should be the same as the parameter + # file). + WL_POLAR_LABEL_POSITION(wd) = INDEF + WL_AXIS_TITLE_SIZE(wd) = 1.5 + WL_LABEL_SIZE(wd) = 1.0 + WL_MAJ_TICK_SIZE(wd) = .03 + WL_MIN_TICK_SIZE(wd) = .01 + WL_TITLE_SIZE(wd) = 2.0 + WL_GRAPH_TYPE(wd) = INDEFI + WL_MAJ_LINE_TYPE(wd) = GL_SOLID + WL_MIN_LINE_TYPE(wd) = GL_DOTTED + WL_TITLE_SIDE(wd) = TOP + WL_ALWAYS_FULL_LABEL(wd) = NO + WL_LABEL_ROTATE(wd) = YES + WL_LABON(wd) = YES + WL_LABOUT(wd) = YES + WL_MAJ_GRIDON(wd) = YES + WL_MIN_GRIDON(wd) = NO + WL_REMEMBER(wd) = NO + WL_TICK_IN(wd) = YES + + # Initialize any strings. + call strcpy ("imtitle", WL_TITLE(wd), SZ_LINE) + + # Initialize the axis dependent values. + do i = 1, N_DIM { + WL_AXIS_TITLE(wd,i) = EOS + WL_AXIS_TITLE_SIDE(wd,i) = INDEFI + WL_BEGIN(wd,i) = INDEFD + WL_END(wd,i) = INDEFD + WL_MAJOR_INTERVAL(wd,i) = INDEFD + WL_MINOR_INTERVAL(wd,i) = 5 + do j = 1, N_SIDES + WL_LABEL_SIDE(wd,j,i) = false + } + + # Return the descriptor. + return (wd) +end + + +# WL_WCS_PARAMS -- Read the WCS descriptor from the parameters. +# +# Description +# This procedure returns the WCS descriptor created from task parameters +# and the logical space that will be graphed. +# +# Bugs +# This only deals with two axes. + +procedure wl_wcs_params (mw, log_x1, log_x2, log_y1, log_y2) + +pointer mw # O: The MWCS descriptor. +real log_x1, log_x2, # O: The extent of the logical space to graph. +real log_y1, log_y2 + +real cd[2,2], r[2], w[2] +pointer sp, input, pp +pointer clopset(), mw_open() +real clgpsetr() + +begin + call smark (sp) + call salloc (input, SZ_LINE, TY_CHAR) + + # Open the pset. + pp = clopset ("wcspars") + + # Create an MWCS descriptor. + mw = mw_open (NULL, 2) + + # Get the types. + call clgpset (pp, "ctype1", Memc[input], SZ_LINE) + call wl_decode_ctype (mw, Memc[input], 1) + call clgpset (pp, "ctype2", Memc[input], SZ_LINE) + call wl_decode_ctype (mw, Memc[input], 2) + + # Get the reference coordinates. + r[1] = clgpsetr (pp, "crpix1") + r[2] = clgpsetr (pp, "crpix2") + w[1] = clgpsetr (pp, "crval1") + w[2] = clgpsetr (pp, "crval2") + + # Get the CD matrix. + cd[1,1] = clgpsetr (pp, "cd1_1") + cd[1,2] = clgpsetr (pp, "cd1_2") + cd[2,1] = clgpsetr (pp, "cd2_1") + cd[2,2] = clgpsetr (pp, "cd2_2") + + # Set the Wterm. + call mw_swtermr (mw, r, w, cd, 2) + + # Get the extent of the logical space. + log_x1 = clgpsetr (pp, "log_x1") + log_x2 = clgpsetr (pp, "log_x2") + log_y1 = clgpsetr (pp, "log_y1") + log_y2 = clgpsetr (pp, "log_y2") + + # Close the pset. + call clcpset (pp) + + call sfree (sp) +end + + +# WL_DECODE_CTYPE -- Decode the ctype string into axis type and system type. +# +# Description +# The CTYPE is what is found in FITS keywords CTYPEn. The value may +# contain two pieces of information, always the system type and possibly +# an individual axis type. For systems such as plain old linear systems +# just a system type is defined. However, for celestial systems, both +# types are defined in the form "axistype-systemtype". There may be +# any number of '-' in between the values. + +procedure wl_decode_ctype (mw, input, axno) + +pointer mw # I: the MWCS descriptor +char input[ARB] # I: the string input +int axno # I: the axis being worked on + +int i, input_len +int strncmp(), strldx(), strlen() +string empty "" + +begin + + input_len = strlen (input) + + # Fix some characters. + do i = 1, input_len { + if (input[i] == ' ' || input[i] == '\'') + break + else if (IS_UPPER(input[i])) + input[i] = TO_LOWER(input[i]) + else if (input[i] == '_') + input[i] = '-' + } + + # Determine the type of function on this axis. + if (strncmp (input, "linear", 6) == 0) { + call mw_swtype (mw, 1, 2, "linear", empty) + + } else if (strncmp (input, "ra--", 4) == 0) { + i = strldx ("-", input) + 1 + call mw_swtype (mw, 1, 2, input[i], empty) + call mw_swattrs (mw, axno, "axtype", "ra") + + } else if (strncmp (input, "dec-", 4) == 0) { + i = strldx ("-", input) + 1 + call mw_swtype (mw, 1, 2, input[i], empty) + call mw_swattrs (mw, axno, "axtype", "dec") + + } else { + # Since we have to be able to read any FITS header, we have + # no control over the value of CTYPEi. If the value is + # something we don't know about, assume a LINEAR axis, using + # the given value of CTYPEi as the default axis label. + call mw_swtype (mw, 1, 2, "linear", empty) + call mw_swattrs (mw, axno, "label", input) + } + +end + + +# WL_GET_SYSTEM_TYPE -- Determine type of transformation the MWCS represents. +# +# Note +# For some systems, the axis mapping reverses the order to make +# the rest of the code tractable. The only problem is that when graphing, +# the graph routines need to "fix" this reversal. Also note that this +# occurs only for systems that have distinct axis types, such as RA and +# DEC. +# +# Bugs +# A potential problem: For a WCS that has more axes than necessary +# for the sky projections, those axis are set such that during +# transformations, the first index position is used. For the one +# example I have seen, the "third" axis is time and this interpretation +# works. But, I am sure something will fall apart because of this. + +procedure wl_get_system_type (mw, system_type, logical_center, world_center, + flip) + +pointer mw # I: the MWCS descriptor. +int system_type # O: the transformation type: + # RA_DEC -> tan, sin, or arc projection + # in right ascension and + # declination + # LINEAR -> any regular linear system + # INDEFI -> could not be determined +double logical_center[N_DIM] # O: the center point in the logical system. +double world_center[N_DIM] # O: the center point in the world system. +int flip # O: true if the order of the axes have been + # changed by axis mappins + +double tmp_logical[MAX_DIM], tmp_world[MAX_DIM] +int wcs_dim, axis, index_sys1, index_sys2, found_axis +int axno[MAX_DIM], axval[MAX_DIM], found_axis_list[N_DIM] +pointer sp, axtype, cd, cur_type +int mw_stati(), strncmp(), strdic() +errchk mw_gwattrs + +begin + # Get some memory. + call smark (sp) + call salloc (axtype, SZ_LINE, TY_CHAR) + call salloc (cur_type, SZ_LINE, TY_CHAR) + call salloc (cd, MAX_DIM, TY_DOUBLE) + + # Get the dimensionality of the WCS. + call mw_seti (mw, MW_USEAXMAP, NO) + wcs_dim = mw_stati (mw, MW_NDIM) + + # Initialize the two dimensions. + index_sys1 = INDEFI + index_sys2 = INDEFI + + # Look through the possible supported axis types. When a type has + # exactly N_DIM axes defined, that will be the one used. + + for (system_type = 1; system_type <= NUMBER_OF_SUPPORTED_TYPES; + system_type = system_type + 1) { + + # Determine the string that should be looked for. + switch (system_type) { + case RA_DEC: + call strcpy (RA_DEC_DICTIONARY, Memc[cur_type], SZ_LINE) + case LINEAR: + call strcpy (LINEAR_DICTIONARY, Memc[cur_type], SZ_LINE) + } + + # Initialize the number of found axis. + found_axis = 0 + + # Examine each axis to determine whether the current axis type is + # the one to use. + for (axis = 1; axis <= wcs_dim; axis = axis + 1) { + + # If the current physical axis is not mapped, ignore it. + # This statement is causing a problem in 2.10.3, not sure + # why but am removing it for now. + #if (axno[axis] == 0) + #next + + ifnoerr (call mw_gwattrs( mw, axis, "wtype", Memc[axtype], + SZ_LINE)) { + call strlwr (Memc[axtype]) + + # If this axis type matches the one being looked for, add + # it to the axis list. If there are too many axis of the + # current type found, don't add to the found axis list. + + if (strdic (Memc[axtype], Memc[axtype], SZ_LINE, + Memc[cur_type]) > 0) { + found_axis = found_axis + 1 + if (found_axis <= N_DIM) + found_axis_list[found_axis] = axis + } + } + } + + # Check to see whether we have the right number axes. + if (found_axis == N_DIM) + break + + } + + # If any axes were found, then further check axis types. + # Depending on the axis type, there may be need to distinguish + # between the two possible axis further. + + if (found_axis == N_DIM) + switch (system_type) { + case RA_DEC: + for (axis = 1; axis <= N_DIM; axis = axis + 1) + ifnoerr (call mw_gwattrs (mw, found_axis_list[axis], + "axtype", Memc[axtype], SZ_LINE)) { + call strlwr( Memc[axtype] ) + if (strncmp (Memc[axtype], "ra", 2) == 0) + index_sys1 = found_axis_list[axis] + else if (strncmp (Memc[axtype], "dec", 3) == 0) + index_sys2 = found_axis_list[axis] + } + + # The "default" seems to be the LINEAR case for MWCS. + # Since no other information is provided, this is all we know. + default: + index_sys1 = found_axis_list[1] + index_sys2 = found_axis_list[2] + } + + # If either axis is unknown, something is wrong. If the WCS has two + # axes defined, then make some grand assumptions. If not, then there + # is nothing more to be done. + + if (IS_INDEFI (index_sys1) || IS_INDEFI (index_sys2)) { + if (wcs_dim >= N_DIM) { + index_sys1 = 1 + index_sys2 = 2 + } else + call error (0, "Wcslab: Fewer than two defined axes") + } + + # Zero the axis values and set any "unknown" axis to always use the + # "first" position in that axis direction. This will more than likely + # be a problem, but no general solution comes to mind this second. + + call amovki (0, axno, wcs_dim) + call amovki (0, axval, wcs_dim) + + # Setup so that the desired axes are set as the X and Y axis. + axno[index_sys1] = X_DIM + axno[index_sys2] = Y_DIM + call mw_saxmap (mw, axno, axval, wcs_dim) + + # Recover the center points of the Logical and World systems. + call mw_gwtermd (mw, tmp_logical, tmp_world, Memd[cd], wcs_dim) + + logical_center[X_DIM] = tmp_logical[index_sys1] + logical_center[Y_DIM] = tmp_logical[index_sys2] + world_center[X_DIM] = tmp_world[index_sys1] + world_center[Y_DIM] = tmp_world[index_sys2] + + # Check for reversal of axes + if (index_sys1 > index_sys2) + flip = YES + else + flip = NO + + # Release the memory. + call sfree (sp) +end + + +# WL_GR_INPARAMS -- Read in the graphics parameters for wcslab. +# +# Description +# Read all the parameters in and make some decisions about what +# will be done. + +procedure wl_gr_inparams (wd) + +pointer wd # I: the WCSLAB descriptor + +pointer sp, aline, pp +bool clgpsetb(), streq() +double wl_string_to_internal() +int btoi(), strdic(), wl_line_type(), clgpseti() +pointer clopset() +real clgpsetr() + +begin + # Get some memory. + call smark (sp) + call salloc (aline, SZ_LINE, TY_CHAR) + + # Open the pset. + pp = clopset ("wlpars") + + # Get the title if other than the default. + call clgpset (pp, "title", Memc[aline], SZ_LINE) + if (! streq (Memc[aline], "imtitle")) + call strcpy (Memc[aline], WL_TITLE(wd), SZ_LINE) + + # Get the axis titles. + call clgpset (pp, "axis1_title", WL_AXIS_TITLE(wd,AXIS1), SZ_LINE) + call clgpset (pp, "axis2_title", WL_AXIS_TITLE(wd,AXIS2), SZ_LINE) + + # Get the parameters. + WL_ALWAYS_FULL_LABEL(wd) = btoi (clgpsetb (pp,"full_label")) + WL_AXIS_TITLE_SIZE(wd) = clgpsetr (pp, "axis_title_size") + WL_LABEL_ROTATE(wd) = btoi (clgpsetb (pp, "rotate")) + WL_LABEL_SIZE(wd) = clgpsetr (pp, "label_size") + WL_LABON(wd) = btoi (clgpsetb (pp, "dolabel")) + WL_LABOUT(wd) = btoi (clgpsetb (pp, "labout")) + WL_MAJ_GRIDON(wd) = btoi (clgpsetb (pp, "major_grid")) + WL_MAJ_TICK_SIZE(wd) = clgpsetr (pp, "major_tick") + WL_MIN_GRIDON(wd) = btoi (clgpsetb (pp, "minor_grid")) + WL_MINOR_INTERVAL(wd,AXIS1) = clgpseti (pp, "axis1_minor") + WL_MINOR_INTERVAL(wd,AXIS2) = clgpseti (pp, "axis2_minor") + WL_MIN_TICK_SIZE(wd) = clgpsetr (pp, "minor_tick") + WL_REMEMBER(wd) = btoi (clgpsetb (pp, "remember")) + WL_TICK_IN(wd) = btoi (clgpsetb (pp, "tick_in")) + WL_TITLE_SIZE(wd) = clgpsetr (pp, "title_size") + + # Set what type of graph will be plotted. + call clgpset (pp, "graph_type", Memc[aline], SZ_LINE) + call strlwr (Memc[aline]) + WL_GRAPH_TYPE(wd) = strdic (Memc[aline], Memc[aline], SZ_LINE, + GRAPHTYPES) + if (WL_GRAPH_TYPE(wd) <= 0) + WL_GRAPH_TYPE(wd) = INDEFI + + # Get which sides labels will appear on. + call clgpset (pp, "axis1_side", Memc[aline], SZ_LINE) + call strlwr (Memc[aline]) + call wl_label_side (Memc[aline], WL_LABEL_SIDE(wd,1,AXIS1)) + + call clgpset (pp, "axis2_side", Memc[aline], SZ_LINE) + call strlwr (Memc[aline]) + call wl_label_side (Memc[aline], WL_LABEL_SIDE(wd,1,AXIS2)) + + # Get the polar justification direction. + call clgpset (pp, "justify", Memc[aline], SZ_LINE) + call strlwr (Memc[aline]) + WL_POLAR_LABEL_DIRECTION(wd) = strdic (Memc[aline], Memc[aline], + SZ_LINE, GRAPHSIDES) + if (WL_POLAR_LABEL_DIRECTION(wd) <= 0) + WL_POLAR_LABEL_DIRECTION(wd) = INDEFI + + # Decode the graphing parameters. + call clgpset (pp, "axis1_int", Memc[aline], SZ_LINE) + WL_MAJOR_INTERVAL(wd,AXIS1) = wl_string_to_internal (Memc[aline], + WL_SYSTEM_TYPE(wd), AXIS1) + call clgpset (pp, "axis1_beg", Memc[aline], SZ_LINE) + WL_BEGIN(wd,AXIS1) = wl_string_to_internal (Memc[aline], + WL_SYSTEM_TYPE(wd), AXIS1) + call clgpset (pp, "axis1_end", Memc[aline], SZ_LINE) + WL_END(wd,AXIS1) = wl_string_to_internal (Memc[aline], + WL_SYSTEM_TYPE(wd), AXIS1) + + call clgpset (pp, "axis2_int", Memc[aline], SZ_LINE) + WL_MAJOR_INTERVAL(wd,AXIS2) = wl_string_to_internal (Memc[aline], + WL_SYSTEM_TYPE(wd), AXIS2) + call clgpset (pp, "axis2_beg", Memc[aline], SZ_LINE) + WL_BEGIN(wd,AXIS2) = wl_string_to_internal(Memc[aline], + WL_SYSTEM_TYPE(wd), AXIS2 ) + call clgpset (pp, "axis2_end", Memc[aline], SZ_LINE) + WL_END(wd,AXIS2) = wl_string_to_internal (Memc[aline], + WL_SYSTEM_TYPE(wd), AXIS2) + + # Get the polar label position. + call clgpset (pp, "axis2_dir", Memc[aline], SZ_LINE) + WL_POLAR_LABEL_POSITION(wd) = wl_string_to_internal( Memc[aline], + WL_SYSTEM_TYPE(wd), AXIS1) + + # Get the axis titles. + call clgpset (pp, "axis1_title_side", Memc[aline], SZ_LINE) + call strlwr (Memc[aline]) + WL_AXIS_TITLE_SIDE(wd,AXIS1) = strdic (Memc[aline], Memc[aline], + SZ_LINE, GRAPHSIDES) + if (WL_AXIS_TITLE_SIDE(wd,AXIS1) <= 0) + WL_AXIS_TITLE_SIDE(wd,AXIS1) = INDEFI + + call clgpset (pp, "axis2_title_side", Memc[aline], SZ_LINE) + call strlwr (Memc[aline]) + WL_AXIS_TITLE_SIDE(wd,AXIS2) = strdic (Memc[aline], Memc[aline], + SZ_LINE, GRAPHSIDES) + if (WL_AXIS_TITLE_SIDE(wd,AXIS2) <= 0) + WL_AXIS_TITLE_SIDE(wd,AXIS2) = INDEFI + + # Decode the grid line types. + call clgpset (pp, "major_line", Memc[aline], SZ_LINE) + WL_MAJ_LINE_TYPE(wd) = wl_line_type (Memc[aline]) + call clgpset (pp, "minor_line", Memc[aline], SZ_LINE) + WL_MIN_LINE_TYPE(wd) = wl_line_type (Memc[aline]) + + # Get the title side. + call clgpset (pp, "title_side", Memc[aline], SZ_LINE) + call strlwr (Memc[ aline]) + WL_TITLE_SIDE(wd) = strdic (Memc[aline], Memc[aline], SZ_LINE, + GRAPHSIDES) + + # Close the pset. + call clcpset (pp) + + # Free memory. + call sfree (sp) +end + + +# WL_GR_REMPARAMS -- Write out the graphing parameters. + +procedure wl_gr_remparams (wd) + +pointer wd # I: the WCSLAB descriptor. + +pointer sp, output, pp +pointer clopset() + +begin + # Get some memory. + call smark (sp) + call salloc (output, SZ_LINE, TY_CHAR) + + # Open the pset. + pp = clopset ("wlpars") + + # Set the graph type. + switch (WL_GRAPH_TYPE(wd)) { + case NORMAL: + call clppset (pp, "graph_type", "normal") + case POLAR: + call clppset (pp, "graph_type", "polar") + case NEAR_POLAR: + call clppset (pp, "graph_type", "near_polar") + default: + call clppset (pp, "graph_type", "default") + } + + # Write back the labelling parameters. + call wl_internal_to_string (WL_MAJOR_INTERVAL(wd,AXIS1), + WL_SYSTEM_TYPE(wd), AXIS1, Memc[output]) + call clppset (pp, "axis1_int", Memc[output]) + call wl_internal_to_string (WL_BEGIN(wd,AXIS1), WL_SYSTEM_TYPE(wd), + AXIS1, Memc[output]) + call clppset (pp, "axis1_beg", Memc[output]) + call wl_internal_to_string (WL_END(WD,AXIS1), WL_SYSTEM_TYPE(wd), + AXIS1, Memc[output]) + call clppset (pp, "axis1_end", Memc[output]) + call wl_internal_to_string (WL_MAJOR_INTERVAL(wd,AXIS2), + WL_SYSTEM_TYPE(wd), AXIS2, Memc[output]) + call clppset (pp, "axis2_int", Memc[output]) + call wl_internal_to_string (WL_BEGIN(wd,AXIS2), WL_SYSTEM_TYPE(wd), + AXIS2, Memc[output]) + call clppset (pp, "axis2_beg", Memc[output]) + call wl_internal_to_string (WL_END(wd,AXIS2), WL_SYSTEM_TYPE(wd), + AXIS2, Memc[output]) + call clppset (pp, "axis2_end", Memc[output]) + call wl_internal_to_string (WL_POLAR_LABEL_POSITION(wd), + WL_SYSTEM_TYPE(wd), AXIS1, Memc[output]) + call clppset (pp, "axis2_dir", Memc[output]) + + # Write back labelling justification. + call wl_side_to_string (WL_POLAR_LABEL_DIRECTION(wd), Memc[output], + SZ_LINE) + call clppset (pp, "justify", Memc[output]) + + # Put the axis title sides out. + call wl_side_to_string (WL_AXIS_TITLE_SIDE(wd,AXIS1), Memc[output], + SZ_LINE) + call clppset (pp, "axis1_title_side", Memc[output]) + call wl_side_to_string (WL_AXIS_TITLE_SIDE(wd,AXIS2), Memc[output], + SZ_LINE ) + call clppset (pp, "axis2_title_side", Memc[output]) + + # Put the label sides out. + call wl_put_label_sides (WL_LABEL_SIDE(wd,1,AXIS1), Memc[output], + SZ_LINE ) + call clppset (pp, "axis1_side", Memc[output]) + call wl_put_label_sides (WL_LABEL_SIDE(wd,1,AXIS2), Memc[output], + SZ_LINE) + call clppset (pp, "axis2_side", Memc[output]) + + # Close the pset. + call clcpset (pp) + + # Free memory. + call sfree (sp) +end + + +# WL_DESTROY -- Deallocate the WCSLAB descriptor. + +procedure wl_destroy (wd) + +pointer wd # I: the WCSLAB descriptor to be destroyed + +begin + # Deallocate all the subarrays. + call mfree (WL_WORLD_CENTER_PTR(wd), TY_DOUBLE) + call mfree (WL_TITLE_PTR(wd), TY_CHAR) + call mfree (WL_SCREEN_BOUNDARY_PTR(wd), TY_DOUBLE) + call mfree (WL_NV_PTR(wd), TY_REAL) + call mfree (WL_MIN_I_PTR(wd), TY_INT) + call mfree (WL_MAJ_I_PTR(wd), TY_DOUBLE) + call mfree (WL_LOGICAL_CENTER_PTR(wd), TY_DOUBLE) + call mfree (WL_LABEL_VALUE_PTR(wd), TY_DOUBLE) + call mfree (WL_LABEL_SIDE_PTR(wd), TY_BOOL) + call mfree (WL_LABEL_POSITION_PTR(wd), TY_DOUBLE) + call mfree (WL_LABEL_AXIS_PTR(wd), TY_INT) + call mfree (WL_LABEL_ANGLE_PTR(wd), TY_DOUBLE) + call mfree (WL_END_PTR(wd), TY_DOUBLE) + call mfree (WL_BEGIN_PTR(wd), TY_DOUBLE) + call mfree (WL_AXIS_TITLE_SIDE_PTR(wd), TY_BOOL) + call mfree (WL_AXIS_TITLE_PTR(wd), TY_CHAR) + + # Now deallocate the structure. + call mfree (wd, TY_STRUCT) +end + + +# WL_LABEL_SIDE -- Decode string into set of booleans sides. + +procedure wl_label_side (input, flag) + +char input[ARB] # I: string listing the sides to be labeled +bool flag[N_SIDES] # O: the flags indicating which sides wll be labeled + +int i +int strmatch() + +begin + # Initialize all the flags to false. + do i = 1, N_SIDES + flag[i] = false + + # Now set each side that is in the list. + if (strmatch (input, "right") != 0) + flag[RIGHT] = true + if (strmatch (input, "left") != 0) + flag[LEFT] = true + if (strmatch (input, "top") != 0) + flag[TOP] = true + if (strmatch (input, "bottom") != 0) + flag[BOTTOM] = true +end + + +# WL_STRING_TO_INTERVAL -- Convert from a string to a number. +# +# Description +# Since (ideally) the wcslab task should be able to handle any sky +# map transformation, there are a number of potential units that can be +# transformed from. The specification of coordinates in these systems +# are also quite varied. Thus, for input purposes, coordinates are entered +# as strings. This routine decodes the strings to a common unit (degrees) +# based on the type of system being graphed. +# +# Function Returns +# This returns the single coordinate value converted to a base system +# (degrees). + +double procedure wl_string_to_internal (input, axis_type, which_axis) + +char input[ARB] # I; the string containing the numerical value +int axis_type # I: the type of wcs +int which_axis # I: the axis number + +double value +int strlen(), nscan() + +begin + # It is possible that the value was not defined. + if (strlen (input) <= 0) + value = INDEFD + + # Decode based on the system. + else + switch (axis_type) { + + # The RA and DEC systems. + case RA_DEC: + + # Since SPP FMTIO can handle the HH:MM:SS format, just let it + # read in the value. However, there is no way to distinquish + # H:M:S from D:M:S. If the axis being read is RA, assume that + # it was H:M:S. + + call sscan (input) + call gargd (value) + + # If the axis is Longitude == RA, then convert the hours to + # degrees. + if (nscan() < 1) { + value = INDEFD + } else { + if (which_axis == AXIS1) + value = HRSTODEG (value) + } + + # Default- unknown system, just read the string as a double + # precision and return it. + default: + call sscan (input) + call gargd (value) + if (nscan() < 1) + value = INDEFD + } + + return (value) +end + + +# WL_LINE_TYPE -- Decode a string into an IRAF GIO polyline type. + +int procedure wl_line_type (line_type_string) + +char line_type_string[ARB] # I: the string specifying the line type + # "solid" -> GL_SOLID + # "dotted" -> GL_DOTTED + # "dashed" -> GL_DASHED + # "dotdash" -> GL_DOTDASH +int type +bool streq() + +begin + if (streq (line_type_string, "solid")) + type = GL_SOLID + else if (streq (line_type_string, "dotted")) + type = GL_DOTTED + else if (streq( line_type_string, "dashed")) + type = GL_DASHED + else if (streq (line_type_string, "dotdash")) + type = GL_DOTDASH + else { + call eprintf ("Pattern unknown, using 'solid'.\n") + type = GL_SOLID + } + + return (type) +end + + +# WL_INTERNAL_TO_STRING - Convert internal representation to a string. + +procedure wl_internal_to_string (value, system_type, which_axis, output) + +double value # I: the value to convert +int system_type # I: the wcs type +int which_axis # I: the axis +char output[ARB] # O: the output string + +begin + # If the value is undefined, write an empty string. + if (IS_INDEFD (value)) + output[1] = EOS + + # Else, convert the value depending on the axis types. + else + switch (system_type) { + + # Handle the RA, DEC + case RA_DEC: + + # If this is Axis1 == Right Ascension, then convert to hours. + if (which_axis == AXIS1) + value = value / 15.0D0 + + call sprintf (output, SZ_LINE, "%.6h") + call pargd (value) + + # Else, just write a value. + default: + call sprintf (output, SZ_LINE, "%.7g") + call pargd (value) + } + +end + + +# WL_SIDE_TO_STRING -- Convert a side to its string representation. + +procedure wl_side_to_string (side, output, max_len) + +int side # I: the side to convert +char output[max_len] # O: the string representation of the side +int max_len # I: the maximum length of the output string + +begin + switch (side) { + case RIGHT: + call strcpy ("right", output, max_len) + case LEFT: + call strcpy ("left", output, max_len) + case TOP: + call strcpy ("top", output, max_len) + case BOTTOM: + call strcpy ("bottom", output, max_len) + default: + call strcpy ("default", output, max_len) + } +end + + +# WL_PUT_LABEL_SIDES -- Create a string containing the sides specified. + +procedure wl_put_label_sides (side_flags, output, max_len) + +bool side_flags[N_SIDES] # I: the boolean array of sides +char output[ARB] # O: the output comma separated list of sides +int max_len # I: maximum length of the output string + +int i +pointer sp, side +int strlen() + +begin + # Get memory. + call smark (sp) + call salloc (side, max_len, TY_CHAR) + + # Build the list. + output[1] = EOS + do i = 1, N_SIDES + if (side_flags[i]) { + if (strlen (output) != 0) + call strcat (",", output, max_len) + call wl_side_to_string (i, Memc[side], max_len) + call strcat (Memc[side], output, max_len) + } + + if (strlen (output) == 0) + call strcat ("default", output, max_len) + + # Free memory. + call sfree (sp) +end diff --git a/pkg/utilities/nttools/stxtools/wcslab/wlgrid.x b/pkg/utilities/nttools/stxtools/wcslab/wlgrid.x new file mode 100644 index 00000000..4f457af4 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wlgrid.x @@ -0,0 +1,448 @@ +include <gset.h> +include <math.h> +include "wcslab.h" +include "wcs_desc.h" + + +# WL_GRID -- Put the grid lines/tick marks on the plot. +# +# Description +# Based on previously determined parameters., draw the grid lines and/or +# tick marks onto the graph. While in the process of doing this, create +# a list of possible label points for use by the label_grid routine. + +procedure wl_grid (wd) + +pointer wd # I: the WCSLAB descriptor + +double current, tmp_begin, tmp_end, tmp_minor_interval +int old_type, old_n_labels, min_counter +int gstati() + +begin + # Initialize the label counter. + WL_N_LABELS(wd) = 0 + + # Remember what line type is currently active. + old_type = gstati (WL_GP(wd), G_PLTYPE) + + # Determine integer range for axis 1. + tmp_minor_interval = WL_MAJOR_INTERVAL(wd,AXIS1) / + double (WL_MINOR_INTERVAL(wd,AXIS1)) + + # If near-polar, the lines should go all the way to the poles. + if (WL_GRAPH_TYPE(wd) == NEAR_POLAR) + if (abs (WL_BEGIN(wd,AXIS2)) < abs (WL_END(wd,AXIS2))) { + tmp_begin = WL_BEGIN(wd,AXIS2) + tmp_end = NORTH_POLE_LATITUDE + } else { + tmp_begin = SOUTH_POLE_LATITUDE + tmp_end = WL_END(wd,AXIS2) + } + else { + tmp_begin = WL_BEGIN(wd,AXIS2) + tmp_end = WL_END(wd,AXIS2) + } + + # Plot lines of constant value in axis 1. + current = WL_BEGIN(wd,AXIS1) + min_counter = 0 + repeat { + + if (mod (min_counter, WL_MINOR_INTERVAL(wd,AXIS1)) == 0) { + call gseti (WL_GP(wd), G_PLTYPE, WL_MAJ_LINE_TYPE(wd)) + call wl_graph_constant_axis1 (wd, current, tmp_begin, tmp_end, + WL_MAJ_GRIDON(wd), WL_LABON(wd), WL_MAJ_TICK_SIZE(wd)) + } else { + call gseti (WL_GP(wd), G_PLTYPE, WL_MIN_LINE_TYPE(wd)) + call wl_graph_constant_axis1 (wd, current, tmp_begin, tmp_end, + WL_MIN_GRIDON(wd), NO, WL_MIN_TICK_SIZE(wd)) + } + + min_counter = min_counter + 1 + current = WL_BEGIN(wd,AXIS1) + tmp_minor_interval * min_counter + + } until (real (current) > real (WL_END(wd,AXIS1))) + + # Determine the interval range for the second axis. + tmp_minor_interval = WL_MAJOR_INTERVAL(wd,AXIS2) / + double (WL_MINOR_INTERVAL(wd,AXIS2)) + + # Plot lines of constant value in axis 2. + if (WL_END(wd,AXIS2) < WL_BEGIN(wd,AXIS2)) { + current = WL_END(wd,AXIS2) + tmp_minor_interval = -tmp_minor_interval + tmp_end = WL_BEGIN(wd,AXIS2) + } else { + current = WL_BEGIN(wd,AXIS2) + tmp_end = WL_END(wd,AXIS2) + } + + min_counter = 0 + tmp_begin = current + repeat { + if (mod (min_counter, WL_MINOR_INTERVAL(wd,AXIS2)) == 0) { + + call gseti (WL_GP(wd), G_PLTYPE, WL_MAJ_LINE_TYPE(wd)) + old_n_labels = WL_N_LABELS(wd) + call wl_graph_constant_axis2 (wd, current, WL_BEGIN(wd,AXIS1), + WL_END(wd,AXIS1), WL_MAJ_GRIDON(wd), WL_LABON(wd), + WL_MAJ_TICK_SIZE(wd)) + + # If this is a polar or near_polar plot, the latitudes + # should be placed near the line, not where it crosses the + # window boundary. + + if (WL_GRAPH_TYPE(wd) == POLAR && + (WL_MAJ_GRIDON(wd) == YES) && (WL_LABON(wd) == YES)) { + WL_N_LABELS(wd) = old_n_labels + 1 + call wl_w2ld (WL_WLCT(wd), WL_AXIS_FLIP(wd), + WL_POLAR_LABEL_POSITION(wd), current, + WL_LABEL_POSITION(wd,WL_N_LABELS(wd),X_DIM), + WL_LABEL_POSITION(wd,WL_N_LABELS(wd),Y_DIM), 1) + WL_LABEL_VALUE(wd,WL_N_LABELS(wd)) = current + WL_LABEL_AXIS(wd,WL_N_LABELS(wd)) = AXIS2 + } + + } else { + call gseti (WL_GP(wd), G_PLTYPE, WL_MIN_LINE_TYPE(wd)) + call wl_graph_constant_axis2 (wd, current, WL_BEGIN(wd,AXIS1), + WL_END(wd,AXIS1), WL_MIN_GRIDON(wd), NO, + WL_MIN_TICK_SIZE(wd)) + } + + # Increment and continue + min_counter = min_counter + 1 + current = tmp_begin + tmp_minor_interval * min_counter + + } until (real (current) > real (tmp_end)) + + # Set the line type back to the way it was. + call gseti (WL_GP(wd), G_PLTYPE, old_type) +end + + +# WL_GRAPH_CONSTANT_AXIS1 - Graph lines of constant X-axis values. +# +# Description +# Because projections are rarely linear, the basic GIO interface to draw +# lines cannot be used. Instead, this routine handles the line drawing. +# Also, possible label points are found and added to a label list array. +# +# CLUDGE! Finding labels here is WRONG. Ideally, crossing points (where the +# line crosses a screen boundary) should be determined analytically. However, +# the MWCS interface lacks the required "cross-transformations". It can +# still be done, but requires a total bypassing of MWCS. Instead, this +# simplistic approach is used. + +procedure wl_graph_constant_axis1 (wd, x, ymin, ymax, gridon, label, tick_size) + +pointer wd # I: the WCSLAB descriptor +double x # I: X value to hold constant +double ymin, ymax # I: Y values to vary between +int gridon # I: true if gridding is on +int label # I: true if the points should be labelled +real tick_size # I: size of tick marks + +bool done +double lastx, lasty, lx, ly, y, yinc +real rlx, rly + +begin + # Determine the scale at which Y should be incremented. + yinc = (ymax - ymin) / WL_LINE_SEGMENTS(wd) + + # Now graph the line segments. + y = ymin + call wl_w2ld (WL_WLCT(wd), WL_AXIS_FLIP(wd), x, y, lastx, lasty, 1) + + rlx = lastx + rly = lasty + call gamove (WL_GP(wd), rlx, rly) + + repeat { + call wl_w2ld (WL_WLCT(wd), WL_AXIS_FLIP(wd), x, y, lx, ly, 1) + call wl_point_to_label (wd, lastx, lasty, lx, ly, AXIS1, x, gridon, + label, tick_size) + if (gridon == YES) { + rlx = lx + rly = ly + call gadraw (WL_GP(wd), rlx, rly) + } + if (yinc < 0.) + done = y < ymax + else + done = y > ymax + y = y + yinc + lastx = lx + lasty = ly + } until (done) +end + + +# WL_GRAPH_CONSTANT_AXIS2 -- Graph lines of constant Y-axis values. +# +# Description +# Because projections are rarely linear, the basic GIO interface to draw +# lines cannot be used. Instead, this routine handles the line drawing. +# Also, possible label points are found and added to an label list array. +# +# CLUDGE! Finding labels here is WRONG. Ideally, crossing points (where the +# line crosses a screen boundary) should be determined analytically. However, +# the MWCS interface lacks the required "cross-transformations". It can +# still be done, but requires a total bypassing of MWCS. Instead, this +# simplistic approach is used. + +procedure wl_graph_constant_axis2 (wd, y, xmin, xmax, gridon, label, tick_size) + +pointer wd # I: the WCSLAB descriptor +double y # I: Y value to hold constant +double xmin, xmax # I: X values to vary between +int gridon # I: true if gridding is on +int label # I: true if points should be labelled +real tick_size # I: tick mark size + +bool done +double lx, ly, lastx, lasty, x, xinc +real rlx, rly + +begin + # Determine the scale at which X should be incremented. + xinc = (xmax - xmin) / WL_LINE_SEGMENTS(wd) + + # Now graph the line segments. + x = xmin + call wl_w2ld (WL_WLCT(wd), WL_AXIS_FLIP(wd), x, y, lastx, lasty, 1) + + rlx = lastx + rly = lasty + call gamove (WL_GP(wd), rlx, rly) + + repeat { + call wl_w2ld (WL_WLCT(wd), WL_AXIS_FLIP(wd), x, y, lx, ly, 1) + call wl_point_to_label (wd, lastx, lasty, lx, ly, AXIS2, y, gridon, + label, tick_size) + if (gridon == YES) { + rlx = lx + rly = ly + call gadraw (WL_GP(wd), rlx, rly) + } + if (xinc < 0.) + done = x < xmax + else + done = x > xmax + lastx = lx + lasty = ly + x = x + xinc + } until (done) +end + + +# Define the inside and outside of the window. + +define OUT (($1<=WL_SCREEN_BOUNDARY(wd,LEFT))||($1>=WL_SCREEN_BOUNDARY(wd,RIGHT))||($2<=WL_SCREEN_BOUNDARY(wd,BOTTOM))||($2>=WL_SCREEN_BOUNDARY(wd,TOP))) + +define IN (($1>WL_SCREEN_BOUNDARY(wd,LEFT))&&($1<WL_SCREEN_BOUNDARY(wd,RIGHT))&&($2>WL_SCREEN_BOUNDARY(wd,BOTTOM))&&($2<WL_SCREEN_BOUNDARY(wd,TOP))) + + +# WL_POINT_TO_LABEL - Record a points position along a window boundary. +# +# Description +# Since the MWCS interface lacks "cross-transformations", i.e. If given +# RA and and X axis location, find DEC and Y axis, we need a different +# method of determining when lines of constant Axis 1/Axis 2 cross +# the window boundary. Since each line is drawn by small increments, each +# increment is watched to see if a window boundary has been crossed. This +# is what this routine does: Confirms that a boundary has been crossed, +# records this position and label value. Tick marks are also drawn here +# because all the necessary information is known at this point. +# +# NOTE: THIS WAY IS A CLUDGE ! A more formal method of finding +# cross-transformations is needed- most likely an iterative method. This +# way was just "convenient at the time". + +procedure wl_point_to_label (wd, x1, y1, x2, y2, axis, axis_value, gridon, + label, tick_size) + +pointer wd # I: the WCSLAB descriptor +double x1, y1, x2, y2 # I: the two possible points to label +int axis # I: which axis are we dealing with ? +double axis_value # I: the value of the axis at this point +int gridon # I: true if gridding is on +int label # I: true if this point should have a label +real tick_size # I: size of the tick mark + +double nx, ny, tick_x, tick_y +double wl_vector_angle() + +begin + # Determine whether the two points straddle a window boundary. If they + # do, then this is the point to label. + if (OUT (x1, y1) && IN (x2, y2)) { + + call wl_axis_on_line (x1, y1, x2, y2, WL_SCREEN_BOUNDARY(wd,1), + nx, ny) + + if (gridon == NO) { + call wl_mark_tick (WL_GP(wd), WL_NDC_WCS(wd), tick_size, + WL_TICK_IN(wd), x1, y1, x2, y2, nx, ny, tick_x, tick_y) + if (WL_TICK_IN(wd) != WL_LABOUT(wd)) { + nx = tick_x + ny = tick_y + } + } + + if ((label == YES) && WL_N_LABELS(wd) < MAX_LABEL_POINTS) { + WL_N_LABELS(wd) = WL_N_LABELS(wd) + 1 + WL_LABEL_POSITION(wd,WL_N_LABELS(wd),AXIS1) = nx + WL_LABEL_POSITION(wd,WL_N_LABELS(wd),AXIS2) = ny + WL_LABEL_VALUE(wd,WL_N_LABELS(wd)) = axis_value + WL_LABEL_AXIS(wd,WL_N_LABELS(wd)) = axis + WL_LABEL_ANGLE(wd,WL_N_LABELS(wd)) = + wl_vector_angle (WL_GP(wd), x1, y1, x2, y2) + } + } + + if (IN (x1, y1) && OUT (x2, y2)) { + + call wl_axis_on_line (x2, y2, x1, y1, WL_SCREEN_BOUNDARY(wd,1), + nx, ny) + + if (gridon == NO) { + call wl_mark_tick (WL_GP(wd), WL_NDC_WCS(wd), tick_size, + WL_TICK_IN(wd), x2, y2, x1, y1, nx, ny, tick_x, tick_y) + if (WL_TICK_IN(wd) != WL_LABOUT(wd)) { + nx = tick_x + ny = tick_y + } + } + + if ((label == YES) && WL_N_LABELS(wd) < MAX_LABEL_POINTS) { + WL_N_LABELS(wd) = WL_N_LABELS(wd) + 1 + WL_LABEL_POSITION(wd,WL_N_LABELS(wd),AXIS1) = nx + WL_LABEL_POSITION(wd,WL_N_LABELS(wd),AXIS2) = ny + WL_LABEL_VALUE(wd,WL_N_LABELS(wd)) = axis_value + WL_LABEL_AXIS(wd,WL_N_LABELS(wd)) = axis + WL_LABEL_ANGLE(wd,WL_N_LABELS(wd)) = + wl_vector_angle (WL_GP(wd), x1, y1, x2, y2) + } + } + +end + + +# WL_MARK_TICK - Draw the tick mark at the point. +# +# Description +# Draw a tick mark rooted at (sx,sy), whose direction is defined by +# the vector (x0,y0) to (x1,y1). The other end of the tick mark is +# returned in (tick_x,tick_y). + +procedure wl_mark_tick (gp, wcs, tick_size, in, x0, y0, x1, y1, sx, sy, + tick_x, tick_y) + +pointer gp # I: the graphics pointer +int wcs # I: the WCS to use to draw the tick marks +real tick_size # I: size of the tick mark +int in # I: true if ticks should be into the graph +double x0, y0, x1, y1 # I: the points defining the tick direction +double sx, sy # I: the root point of the tick mark +double tick_x, tick_y # O: the end point of the tick mark + +int old_line, old_wcs +real dx, dy, t, ndc_x0, ndc_y0, ndc_x1, ndc_y1, ndc_x2, ndc_y2 +real ndc_sx, ndc_sy +int gstati() +real wl_distancer() + +begin + # Change graphics coordinates to NDC. + old_wcs = gstati (gp, G_WCS) + old_line = gstati (gp, G_PLTYPE) + call gseti (gp, G_WCS, wcs) + call gseti (gp, G_PLTYPE, GL_SOLID) + + # Convert the points to NDC coordinates. + ndc_x2 = real (sx) + ndc_y2 = real (sy) + call gctran (gp, ndc_x2, ndc_y2, ndc_sx, ndc_sy, old_wcs, wcs) + ndc_x2 = real (x0) + ndc_y2 = real (y0) + call gctran (gp, ndc_x2, ndc_y2, ndc_x0, ndc_y0, old_wcs, wcs) + ndc_x2 = real (x1) + ndc_y2 = real (y1) + call gctran (gp, ndc_x2, ndc_y2, ndc_x1, ndc_y1, old_wcs, wcs) + + # Determine the parameterized line parameters. + dx = ndc_x1 - ndc_x0 + dy = ndc_y1 - ndc_y0 + + # Determine how large in "time" the tick mark is. + t = tick_size / wl_distancer (ndc_x0, ndc_y0, ndc_x1, ndc_y1) + + # If tick marks are to point out of the graph, reverse the sign of t. + # Also need to turn clipping off for the ticks appear. + if (in == NO) { + t = -t + call gseti (gp, G_CLIP, NO) + } + + # Determine the end point of the tick mark. + ndc_x2 = t * dx + ndc_sx + ndc_y2 = t * dy + ndc_sy + + # Now draw the tick mark. + call gamove (gp, ndc_sx, ndc_sy) + call gadraw (gp, ndc_x2, ndc_y2) + + # Restore clipping if necessary. + if (in == NO) + call gseti (gp, G_CLIP, YES) + + # Restore previous settings. + call gseti (gp, G_WCS, old_wcs) + call gseti (gp, G_PLTYPE, old_line) + + # Transform the end of the tick mark. + call gctran (gp, ndc_x2, ndc_y2, dx, dy, wcs, old_wcs) + tick_x = double (dx) + tick_y = double (dy) +end + + +# WL_VECTOR_ANGLE -- Return the angle represented by the given vector. +# +# Returns +# The angle of the given vector. + +double procedure wl_vector_angle (gp, x1, y1, x2, y2) + +pointer gp # I: the graphics descriptor +double x1, y1, x2, y2 # I: the end points of the vector + +double dangle +real angle, delx, dely, ndc_x1, ndc_x2, ndc_y1, ndc_y2 +bool fp_equalr() +int gstati() + +begin + # Translate the input points to NDC coordinates. + ndc_x1 = real (x1) + ndc_x2 = real (x2) + ndc_y1 = real (y1) + ndc_y2 = real (y2) + call gctran (gp, ndc_x1, ndc_y1, ndc_x1, ndc_y1, gstati (gp, G_WCS), + NDC_WCS) + call gctran (gp, ndc_x2, ndc_y2, ndc_x2, ndc_y2, gstati (gp, G_WCS), + NDC_WCS) + + dely = ndc_y2 - ndc_y1 + delx = ndc_x2 - ndc_x1 + if (fp_equalr (delx, 0.) && fp_equalr (dely, 0.)) + angle = 0.0 + else + angle = RADTODEG (atan2 (dely, delx)) + dangle = angle + + return (dangle) +end diff --git a/pkg/utilities/nttools/stxtools/wcslab/wllabel.x b/pkg/utilities/nttools/stxtools/wcslab/wllabel.x new file mode 100644 index 00000000..4578f89c --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wllabel.x @@ -0,0 +1,1100 @@ +include <gset.h> +include <math.h> +include "psiescape.h" +include "wcslab.h" +include "wcs_desc.h" + + +# Define the offset array. +define OFFSET Memr[$1+$2-1] + +# Define the postscript kernel +define PSIKERN "psikern" + +# WL_LABEL -- Place the labels on the grids. +# +# Description +# Format and write the labels for the grid/tick marks. Much of this +# is wading through conditions to decide whether a label should be +# written or not. + +procedure wl_label (wd) + +pointer wd # I: the WCSLAB descriptor + +bool no_side_axis1, no_side_axis2, streq() +int i, axis1_side, axis2_side +short flag +pointer kernel, sp, offset_ptr +real offset + +begin + # Get some memory. + call smark (sp) + call salloc (offset_ptr, N_SIDES, TY_REAL) + do i = 1, N_SIDES + OFFSET(offset_ptr,i) = 0. + call salloc (kernel, SZ_LINE, TY_CHAR ) + + # Decide whether any sides were specified for either axis. + no_side_axis1 = true + no_side_axis2 = true + do i = 1, N_SIDES { + if (WL_LABEL_SIDE(wd,i,AXIS1)) + no_side_axis1 = false + if (WL_LABEL_SIDE(wd,i,AXIS2)) + no_side_axis2 = false + } + + # If polar, then label the axis 2's next to their circles on the + # graph and allow the Axis 1s to be labeled on all sides of the graph. + + if (WL_GRAPH_TYPE(wd) == POLAR) { + + call wl_polar_label (wd) + + if (no_side_axis1) { + do i = 1, N_SIDES { + WL_LABEL_SIDE(wd,i,AXIS1) = true + } + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(WD,AXIS1))) + WL_AXIS_TITLE_SIDE(WD,AXIS1) = BOTTOM + } + + # If we are near-polar, label the Axis 2 as if polar, and label + # Axis1 on all sides except the side closest to the pole. + + } else if (WL_GRAPH_TYPE(wd) == NEAR_POLAR) { + + if (no_side_axis1) { + WL_LABEL_SIDE(wd,WL_BAD_LABEL_SIDE(wd),AXIS1) = true + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(wd,AXIS1))) + WL_AXIS_TITLE_SIDE(wd,AXIS1) = WL_BAD_LABEL_SIDE(wd) + } + + if (no_side_axis2) { + WL_LABEL_SIDE(wd,WL_POLAR_LABEL_DIRECTION(wd),AXIS2) = true + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(wd,AXIS2))) + WL_AXIS_TITLE_SIDE(wd,AXIS2) = WL_POLAR_LABEL_DIRECTION(wd) + } + + # Final case- adjacent sides should be labelled. + + } else { + + # Determine the best sides for labelling. + if (INVERT (WL_ROTA(wd))) { + axis1_side = LEFT + axis2_side = BOTTOM + } else { + axis1_side = BOTTOM + axis2_side = LEFT + } + + # If no sides were specified, use the calculated ones above. + if (no_side_axis1) + WL_LABEL_SIDE(wd,axis1_side,AXIS1) = true + if (no_side_axis2) + WL_LABEL_SIDE(wd,axis2_side,AXIS2) = true + } + + # Check to see if this is a psikern printer. If so, set text + # so that it is mono-spaced. The superscripting algorithm + # doesn't work too well in a proportional-spaced system. + call ggets (WL_GP(wd), "tn", Memc[kernel], SZ_LINE ) + if (streq (Memc[kernel], PSIKERN)) { + flag = NO + call gescape (WL_GP(wd), PS_VARIABLE_SPACE, flag, + PS_VARIABLE_SPACE_SIZE) + } + + # Now draw the labels for axis 1. + do i = 1, N_SIDES { + + if (WL_LABEL_SIDE(wd,i,AXIS1)) { + call wl_lab_edges (wd, AXIS1, i, offset) + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(WD,AXIS1))) + WL_AXIS_TITLE_SIDE(WD,AXIS1) = i + } else + offset = 0. + + # Modify the bounding box for the new viewport. + if (abs (offset) > abs (OFFSET(offset_ptr,i))) + OFFSET(offset_ptr,i) = offset + } + + # Draw the labels for axis 2. + if (WL_GRAPH_TYPE(wd) != POLAR) + do i = 1, N_SIDES { + + if (WL_LABEL_SIDE(wd,i,AXIS2)) { + call wl_lab_edges (wd, AXIS2, i, offset) + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(wd,AXIS2))) + WL_AXIS_TITLE_SIDE(wd,AXIS2) = i + } else + offset = 0. + + # Modify the bounding box for the new viewport. + if (abs (offset) > abs (OFFSET(offset_ptr,i))) + OFFSET(offset_ptr,i) = offset + } + + # Reset to variable spacing. + if (streq (Memc[kernel], PSIKERN)) { + flag = YES + call gescape (WL_GP(wd), PS_VARIABLE_SPACE, flag, + PS_VARIABLE_SPACE_SIZE) + } + + # Set the bounding box. + do i = 1, N_SIDES + WL_NEW_VIEW(wd,i) = WL_NEW_VIEW(wd,i) + OFFSET(offset_ptr,i) + + # Now write the graph title. + call wl_title (WL_GP(wd), WL_AXIS_TITLE(wd,AXIS1), + WL_AXIS_TITLE_SIDE(wd,AXIS1), WL_AXIS_TITLE_SIZE(wd), + WL_NEW_VIEW(wd,1)) + if (WL_GRAPH_TYPE(wd) != POLAR) + call wl_title (WL_GP(wd), WL_AXIS_TITLE(wd,AXIS2), + WL_AXIS_TITLE_SIDE(wd,AXIS2), WL_AXIS_TITLE_SIZE(WD), + WL_NEW_VIEW(wd,1)) + if (! IS_INDEFI (WL_TITLE_SIDE(wd))) + call wl_title (WL_GP(wd), WL_TITLE(wd), WL_TITLE_SIDE(wd), + WL_TITLE_SIZE(wd), WL_NEW_VIEW(wd,1)) + + # Release memory. + call sfree (sp) +end + + +# Define what is in the screen. + +define IN (($1>WL_SCREEN_BOUNDARY(wd,LEFT))&&($1<WL_SCREEN_BOUNDARY(wd,RIGHT))&&($2>WL_SCREEN_BOUNDARY(wd,BOTTOM))&&($2<WL_SCREEN_BOUNDARY(wd,TOP))) + +# WL_POLAR_LABEL -- Place Latitude labels next to Latitude circles. +# +# Description +# Since Lines of constant Latitude on a polar graph are usually circles +# around the pole, the lines may never cross edges. Instead, the labels +# are placed next to circles. The grid-drawing routines should setup +# the label position array such that each line has only one label point. + +procedure wl_polar_label (wd) + +pointer wd # I: the WCSLAB descriptor + +int i, prec +pointer sp, label, units, label_format, units_format +real char_height, char_width, ndc_textx, ndc_texty, old_text_size +real textx, texty +int wl_precision() +real gstatr(), ggetr() + +begin + # Get some memory. + call smark (sp) + call salloc (label, SZ_LINE, TY_CHAR) + call salloc (units, SZ_LINE, TY_CHAR) + call salloc (label_format, SZ_LINE, TY_CHAR) + call salloc (units_format, SZ_LINE, TY_CHAR) + + # Get the character height and width. This is used to ensure that we + # have moved the label strings off the border. + + char_height = ggetr (WL_GP(wd), "ch") * gstatr (WL_GP(wd), G_TXSIZE) / + 2. + char_width = ggetr (WL_GP(wd), "cw") * gstatr (WL_GP(wd), G_TXSIZE) / + 2. + + # Get the text size and cut it in half for on the plot labelling. + old_text_size = gstatr (WL_GP(wd), G_TXSIZE) + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size) + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size * 0.80) + + # Determine the precision of the output. + prec = wl_precision (wd, AXIS2) + + # Place the labels. + for (i = 1; i <= WL_N_LABELS(wd); i = i + 1) + if (WL_LABEL_AXIS(wd,i) == AXIS2) { + + # Decode the coordinate into a text string. + call wl_dms (WL_LABEL_VALUE(wd,i), Memc[label], Memc[units], + SZ_LINE, prec, true) + + # Convert text position from "unknown" coordinates to NDC. + call gctran (WL_GP(wd), real (WL_LABEL_POSITION(wd,i,AXIS1)), + real (WL_LABEL_POSITION(wd,i,AXIS2)), ndc_textx, ndc_texty, + WL_PLOT_WCS(wd), WL_NDC_WCS(wd)) + + # Determine the text justification. + switch (WL_POLAR_LABEL_DIRECTION(wd)) { + case BOTTOM: + call strcpy ("h=c;v=t", Memc[label_format], SZ_LINE) + call strcpy ("h=c;v=c", Memc[units_format], SZ_LINE) + ndc_texty = ndc_texty - char_height + case TOP: + call strcpy ("h=c;v=c", Memc[label_format], SZ_LINE) + call strcpy ("h=c;v=b", Memc[units_format], SZ_LINE) + ndc_texty = ndc_texty + char_height + case LEFT: + call strcpy ("h=r;v=c", Memc[label_format], SZ_LINE) + call strcpy ("h=r;v=b", Memc[units_format], SZ_LINE) + ndc_textx = ndc_textx - char_width + case RIGHT: + call strcpy ("h=l;v=c", Memc[label_format], SZ_LINE) + call strcpy ("h=l;v=b", Memc[units_format], SZ_LINE) + ndc_textx = ndc_textx + char_width + } + + # Convert the text position from NDC back to the "unknown" + # system. + call gctran (WL_GP(wd), ndc_textx, ndc_texty, textx, texty, + WL_NDC_WCS(wd), WL_PLOT_WCS(wd)) + + # Print the label. + if (IN (textx, texty)) { + call gtext (WL_GP(wd), textx, texty, Memc[label], + Memc[label_format]) + call gtext (WL_GP(wd), textx, texty, Memc[units], + Memc[units_format]) + } + + } + + # Set the text size back. + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size) + + # Release memory. + call sfree (sp) + +end + + +# Memory management for labels + +define LABEL_LIST Memi[labels+$1-1] + +# WL_LAB_EDGES -- Place labels along the edges of the window. +# +# Description +# Place labels on the specified side of the graph. + +procedure wl_lab_edges (wd, axis, side, offset) + +pointer wd # I: the WCSLAB descriptor +int axis # I: the type of axis being labeled +int side # I: the side to place the labels +real offset # O: offset in NDC units for titles + +bool do_full +double angle, tangle +int i, full_label, nlabels, old_wcs, prec +pointer sp, labels +real ndc_textx, ndc_texty, old_text_size, textx, texty + +int wl_full_label_position(), wl_find_side() +double wl_string_angle(), wl_angle() +int gstati(), wl_precision() +real gstatr() + +begin + call smark (sp) + + # All label placement is done in NDC coordinates. + old_wcs = gstati (WL_GP(wd), G_WCS) + call gseti (WL_GP(wd), G_WCS, WL_NDC_WCS(wd)) + + # Set text labelling size. + old_text_size = gstatr (WL_GP(wd), G_TXSIZE) + call gsetr (WL_GP(wd), G_TXSIZE, WL_LABEL_SIZE(wd)) + + # Get the precision of the axis interval. + prec = wl_precision (wd, axis) + + # Initialize string size. + offset = 0. + + # Build a list of possible labels for this side. The conditions are + # that the label should be for the current axis and that it lies on + # the current side. + + call salloc (labels, WL_N_LABELS(wd), TY_INT) + nlabels = 0 + for (i = 1; i <= WL_N_LABELS(wd); i = i + 1) + if (WL_LABEL_AXIS(wd,i) == axis && + wl_find_side (WL_LABEL_POSITION(wd,i,AXIS1), + WL_LABEL_POSITION(wd,i,AXIS2), + WL_SCREEN_BOUNDARY(wd,1)) == side) { + nlabels = nlabels + 1 + LABEL_LIST(nlabels) = i + } + + # If no labels found, then just forget it. If labels found, well + # write them out. + + if (nlabels != 0) { + + # Determine which label should be written out in full. + full_label = wl_full_label_position (wd, Memi[labels], nlabels, + axis, side, prec) + + # Determine the angle that all the labels will be written at. + if ((WL_LABOUT(wd) == NO) && (WL_GRAPH_TYPE(wd) != NORMAL) && + (WL_LABEL_ROTATE(wd) == YES)) + angle = INDEFR + else if ((WL_GRAPH_TYPE(wd) == NORMAL) && ((WL_LABEL_ROTATE(wd) == + YES) || ((WL_LABOUT(wd) == NO) && (WL_MAJ_GRIDON(wd) == YES)))) + angle = wl_angle (wd, Memi[labels], nlabels) + else + angle = 0.0 + + # Place the labels. + for (i = 1; i <= nlabels; i = i + 1) { + + # Save some pertinent information. + textx = real (WL_LABEL_POSITION(wd,LABEL_LIST(i),AXIS1)) + texty = real (WL_LABEL_POSITION(wd,LABEL_LIST(i),AXIS2)) + do_full = ((LABEL_LIST(i) == full_label) || + (WL_ALWAYS_FULL_LABEL(wd) == YES)) + + # Transform the "unknown" coordinate system to a known + # coordinate system, NDC, for text placement. + call gctran (WL_GP(wd), textx, texty, ndc_textx, ndc_texty, + old_wcs, WL_NDC_WCS(wd)) + + # If angle is undefined, determine the angle for each label. + if (IS_INDEFR(angle)) + tangle = wl_string_angle (WL_LABEL_ANGLE(wd, + LABEL_LIST(i)), WL_LABOUT(wd)) + else + tangle = angle + + # Format and write the label. + call wl_write_label (wd, WL_LABEL_VALUE(wd,LABEL_LIST(i)), + side, ndc_textx, ndc_texty, tangle, axis, prec, do_full, + offset) + } + } + + # Reset the graphics WCS. + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size) + call gseti (WL_GP(wd), G_WCS, old_wcs) + + call sfree (sp) +end + + +# WL_TITLE - Write the title of the graph. + +procedure wl_title (gp, title, side, size, viewport) + +pointer gp # I: the graphics descriptor +char title[ARB] # I: the title to write +int side # I: which side the title will go +real size # I: the character size to write the title +real viewport[N_SIDES] # I: the viewport in NDC to keep the title out of + +int old_wcs +real char_height, char_width, left, right, top, bottom, old_rotation +real old_text_size, x, y +int gstati(), strlen() +real ggetr(), gstatr() + +begin + # Make sure there is a title to write. If not, then punt. + if (strlen (title) <= 0) + return + + # Get/Set pertinent graphics info. + call ggview (gp, left, right, bottom, top) + + old_text_size = gstatr (gp, G_TXSIZE) + call gsetr (gp, G_TXSIZE, size) + old_rotation = gstatr (gp, G_TXUP) + + char_height = ggetr (gp, "ch") * size + char_width = ggetr (gp, "cw") * size + + old_wcs = gstati (gp, G_WCS) + call gseti (gp, G_WCS, NDC_WCS) + + # Depending on side, set text position and rotation. + switch (side) { + case TOP: + call gsetr (gp, G_TXUP, 90.) + x = (right + left) / 2. + y = viewport[TOP] + (2 * char_height) + viewport[TOP] = y + (char_height / 2.) + case BOTTOM: + call gsetr (gp, G_TXUP, 90.) + x = (right + left) / 2. + y = viewport[BOTTOM] - (2 * char_height) + viewport[BOTTOM] = y - (char_height / 2.) + case RIGHT: + call gsetr (gp, G_TXUP, 180.) + x = viewport[RIGHT] + (2 * char_width) + y = (top + bottom) / 2. + viewport[RIGHT] = x + (char_width / 2.) + case LEFT: + call gsetr (gp, G_TXUP, 180.) + x = viewport[LEFT] - (2 * char_width) + y = (top + bottom) / 2. + viewport[LEFT] = x - (char_width / 2.) + } + + # Write the puppy out. + call gtext (gp, x, y, title, "h=c;v=c") + + # Set the graphics state back. + call gseti (gp, G_WCS, old_wcs) + call gsetr (gp, G_TXSIZE, old_text_size) + call gsetr (gp, G_TXUP, old_rotation) +end + + +# WL_PRECISION -- Determine the precision of the interval. + +int procedure wl_precision (wd, axis) + +pointer wd # I: the WCSLAB descriptor +int axis # I: which axis is being examined ? + +int prec + +begin + # Handle the sky coordinates. + if (WL_SYSTEM_TYPE(wd) == RA_DEC) + + if (axis == AXIS1) { + if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (3600.0D0)) + prec = HOUR + else if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (60.0D0)) + prec = MINUTE + else if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (1.0D0)) + prec = SECOND + else if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (.01D0)) + prec = SUBSEC_LOW + else + prec = SUBSEC_HIGH + } else { + if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (3600.0D0)) + prec = DEGREE + else if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (60.0D0)) + prec = MINUTE + else if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (1.0D0)) + prec = SECOND + else if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (.01D0)) + prec = SUBSEC_LOW + else + prec = SUBSEC_HIGH + } + + # Handle other coordinate types. + else + prec = INDEFI + + return (prec) + +end + + +# Define some value constraints. + +define LOW_ACCURACY .01 +define HIGH_ACCURACY .0001 + +# WL_HMS -- Convert value to number in hours, minutes, and seconds. + +procedure wl_hms (rarad, hms, units, maxch, precision, all) + +double rarad # I: the value to format into a string (degrees) +char hms[ARB] # O: string containing formatted value +char units[ARB] # O: string containing formatted units +int maxch # I: the maximum number of characters allowed +int precision # I: how precise the output should be +bool all # I: true if all relevent fields should be formatted + +double accuracy, fraction +int sec, h, m, s +pointer sp, temp_hms, temp_units + +begin + # Get some memory. + call smark (sp) + call salloc (temp_hms, maxch, TY_CHAR) + call salloc (temp_units, maxch, TY_CHAR) + + units[1] = EOS + hms[1] = EOS + + # Define how close to zero is needed. + accuracy = LOW_ACCURACY + if (precision == SUBSEC_HIGH) + accuracy = HIGH_ACCURACY + + # Seconds of time. + fraction = double (abs(DEGTOST (rarad))) + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + sec = int (fraction) + fraction = fraction - double (sec) + } else { + sec = int (fraction + 0.5) + fraction = 0. + } + + # Range: 0 to 24 hours. + if (sec < 0) + sec = sec + STPERDAY + else if (sec >= STPERDAY) + sec = mod (sec, STPERDAY) + + # Separater fields. + s = mod (sec, 60) + m = mod (sec / 60, 60) + h = sec / 3600 + + # Format fields. + + # Subseconds. + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + fraction = s + fraction + if (precision == SUBSEC_LOW) { + call sprintf (hms, 6, "%05.2f") + call pargd (fraction) + call strcpy (" s ", units, maxch) + } else { + call sprintf (hms, 8, "%07.4f") + call pargd (fraction) + call strcpy (" s ", units, maxch) + } + if (!all) + all = (fraction < accuracy) + + # Seconds + } else if (precision == SECOND) { + + # NOTE: The all is not part of the if statement because if + # SUBSEC's have been printed, then seconds have already been + # dealt with. If SUBSEC's have not been dealt with, then this + # is the first field to be checked anyways. + + call sprintf (hms, 3, "%02d ") + call pargi (s) + call strcpy (" s", units, maxch) + if (! all) + all = (s == 0) + } + + # Minutes. + if (precision == MINUTE || (precision > MINUTE && all)) { + if (all) { + call strcpy (hms, Memc[temp_hms], maxch) + call strcpy (units, Memc[temp_units], maxch) + } + call sprintf (hms, 3, "%02d ") + call pargi (m) + call strcpy (" m", units, maxch) + if (all) { + call strcat (Memc[temp_hms], hms, maxch) + call strcat (Memc[temp_units], units, maxch) + } else + all = (m == 0) + } + + # Non-zero hours. + if (precision == HOUR || all) { + if (all) { + call strcpy (hms, Memc[temp_hms], maxch) + call strcpy (units, Memc[temp_units], maxch) + } + call sprintf (hms, 3, "%2.2d ") + call pargi (h) + call strcpy(" h", units, maxch) + if (all) { + call strcat (Memc[temp_hms], hms, maxch) + call strcat (Memc[temp_units], units, maxch) + } + } + + # Release memory + call sfree (sp) +end + + +# WL_DMS - Convert value to number in degrees, minutes, and seconds. + +procedure wl_dms (arcrad, dms, units, maxch, precision, all) + +double arcrad # I: the value to format into a string (degrees) +char dms[ARB] # O: string containing formatted value +char units[ARB] # O: string containing formatted units +int maxch # I: the maximum number of characters allowed +int precision # I: how precise the output should be ? +bool all # I: true if all relavent fields should be formatted + +double accuracy, fraction +int sec, h, m, s +pointer sp, temp_dms, temp_units +int strlen() + +begin + # Get some memory. + call smark (sp) + call salloc (temp_dms, maxch, TY_CHAR) + call salloc (temp_units, maxch, TY_CHAR) + + units[1] = EOS + dms[1] = EOS + + # Define how close to zero is needed. + accuracy = LOW_ACCURACY + if (precision == SUBSEC_HIGH) + accuracy = HIGH_ACCURACY + + # Seconds of time. + fraction = double (abs (DEGTOSA (arcrad))) + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + sec = int (fraction) + fraction = fraction - double (sec) + } else { + sec = nint (fraction) + fraction = 0. + } + + # Separater fields. + s = mod (abs(sec), 60) + m = mod (abs(sec) / 60, 60) + h = abs(sec) / 3600 + + # Format fields + # + # Subseconds. + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + + fraction = s + fraction + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + if (precision == SUBSEC_LOW) { + call sprintf (dms, 6, "%05.2f\"") + call pargd (fraction) + call strcpy (" ", units, maxch) + } else { + call sprintf (dms, 8, "%07.4f\"") + call pargd (fraction) + call strcpy (" ", units, maxch) + } + if (! all) + all = (fraction < accuracy) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + + # Seconds + } else if (precision == SECOND) { + + # NOTE: The all is not part of the if statement because if + # SUBSEC's have been printed, then seconds have already been + # dealt with. If SUBSEC's have not been dealt with, then this + # is the first field to be checked anyways. + + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + call sprintf (dms, 3, "%02d\"") + call pargi (s) + call strcpy (" ", units, maxch) + if (! all) + all = (s == 0) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + } + + # Minutes. + if (precision == MINUTE || (precision > MINUTE && all)) { + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + call sprintf (dms, 3, "%02d'") + call pargi (m) + call strcpy (" ", units, maxch) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + if (! all) + all = (m == 0) + } + + # Hours. + if (precision == DEGREE || all) { + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + if (sec + fraction < accuracy) + call strcpy (" 0 ", dms, maxch) + else if (arcrad < 0.) { + call sprintf (dms, 4, "-%d ") + call pargi (h) + } else { + call sprintf (dms, 4, "+%d ") + call pargi (h) + } + call sprintf(units, 4, "%*wo") + call pargi (strlen (dms) - 1) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + } + + # Release memory. + call sfree (sp) +end + + +# WL_FULL_LABEL_POSTION -- Find the position where the full label should be. +# +# Description +# This routine returns the index to the label that should be printed +# in its full form, regardless of its value. This is so there is always +# at least one labelled point with the full information. This point is +# choosen by examining which label is the closest to the passed point +# (usually one of the four corners of the display). +# +# Returns +# Index into the labell arrays of the label to be fully printed. +# If the return index is 0, then there are no labels for the given +# side. + +int procedure wl_full_label_position (wd, labels, nlabels, axis, side, + precision) + +pointer wd # I: the WCSLAB descriptor +int labels[nlabels] # I: array of indexes of labels to be printed +int nlabels # I: the number of labels in labels +int axis # I: the axis being dealt with +int side # I: the side being dealt with +int precision # I: precision of the label + +bool all +double cur_dist, dist +int i, cur_label, xside, yside +pointer sp, temp1 +double wl_distanced() + +begin + # Allocate some working space. + call smark (sp) + call salloc (temp1, SZ_LINE, TY_CHAR) + + # Initialize. + xside = INDEFI + yside = INDEFI + + # Determine which corner will have the full label. + if (side == TOP || side == BOTTOM) { + yside = side + if (axis == AXIS1) { + if (WL_LABEL_SIDE(wd,RIGHT,AXIS2)) + xside = RIGHT + if (WL_LABEL_SIDE(wd,LEFT,AXIS2)) + xside = LEFT + } else { + if (WL_LABEL_SIDE(wd,RIGHT,AXIS1)) + xside = RIGHT + if (WL_LABEL_SIDE(wd,LEFT,AXIS1)) + xside = LEFT + } + if (IS_INDEFI (xside)) + xside = LEFT + } else { + xside = side + if (axis == AXIS1) { + if (WL_LABEL_SIDE(wd,TOP,AXIS2)) + yside = TOP + if (WL_LABEL_SIDE(wd,BOTTOM,AXIS2)) + yside = BOTTOM + } else { + if (WL_LABEL_SIDE(wd,TOP,AXIS1)) + yside = TOP + if (WL_LABEL_SIDE(wd,BOTTOM,AXIS1)) + yside = BOTTOM + } + if (IS_INDEFI (yside)) + yside = BOTTOM + } + + # Find the full label. + cur_label = labels[1] + cur_dist = wl_distanced (WL_SCREEN_BOUNDARY(wd,xside), + WL_SCREEN_BOUNDARY(wd,yside), + WL_LABEL_POSITION(wd,cur_label,AXIS1), + WL_LABEL_POSITION(wd,cur_label,AXIS2)) + + # Now go through the rest of the labels to find a closer label. + for (i = 2; i <= nlabels; i = i + 1) { + + # Check to see if the label would be written in full anyways. + all = false + if (WL_SYSTEM_TYPE(wd) == RA_DEC) { + if (WL_LABEL_AXIS(wd, labels[i]) == LONGITUDE) + call wl_hms (WL_LABEL_VALUE(wd, labels[i]), + Memc[temp1], Memc[temp1], SZ_LINE, precision, all) + else + call wl_dms (WL_LABEL_VALUE(wd, labels[i]), + Memc[temp1], Memc[temp1], SZ_LINE, precision, all) + } + + # If so, don't figure out which label should be full, there + # will be one someplace. + if (all) { + cur_label = INDEFI + break + } + + dist = wl_distanced (WL_SCREEN_BOUNDARY(wd,xside), + WL_SCREEN_BOUNDARY(wd,yside), + WL_LABEL_POSITION(wd,labels[i],AXIS1), + WL_LABEL_POSITION(wd,labels[i],AXIS2)) + if (dist < cur_dist) { + cur_dist = dist + cur_label = labels[i] + } + } + + # Release memory. + call sfree (sp) + + # Return the label index. + return (cur_label) +end + + +# WL_WRITE_LABEL - Write the label in the format specified by the WCS type. + +procedure wl_write_label (wd, value, side, x, y, angle, axis, precision, + do_full, offset) + +pointer wd # I: the WCSLAB descriptor +double value # I: the value to use as the label +int side # I: the side the label is going on +real x, y # I: position of the label in NDC coordinates +double angle # I: the angle the text should be written at +int axis # I: which axis is being labelled +int precision # I: level of precision for labels +bool do_full # I: true if the full label should be printed +real offset # I/O: offset for titles in NDC units + +int tside +pointer sp, label, label_format, units, units_format +real char_height, char_width, in_off_x, in_off_y, length +real lx, ly, new_offset, rx, ry, text_angle +real unit_off_x, unit_off_y, ux, uy + +bool fp_equalr() +double wl_string_angle() +int wl_opposite_side(), strlen() +real ggetr(), gstatr() + +begin + # Get some memory. + call smark (sp) + call salloc (label, SZ_LINE, TY_CHAR) + call salloc (units, SZ_LINE, TY_CHAR) + call salloc (label_format, SZ_LINE, TY_CHAR) + call salloc (units_format, SZ_LINE, TY_CHAR) + + # Get character size. This info is used to move the character string + # by the appropriate amounts. + + char_height = ggetr (WL_GP(wd), "ch") * gstatr (WL_GP(wd), G_TXSIZE) + char_width = ggetr (WL_GP(wd), "cw") * gstatr (WL_GP(wd), G_TXSIZE) + + # Determine the "corrected" angle to write text in. + text_angle = wl_string_angle (angle, WL_LABOUT(wd)) + + # Determine the units offset. + call wl_rotate (0., char_height / 2., 1, text_angle - 90., unit_off_x, + unit_off_y) + + # If the labels are to appear inside the graph and the major grid lines + # have been drawn, then determine the necessary offset to get the label + # off the line. + + if ((WL_LABOUT(wd) == NO) && (WL_MAJ_GRIDON(wd) == YES)) + call wl_rotate (0., 0.75 * char_height, 1, text_angle - 90., + in_off_x, in_off_y) + else { + in_off_x = 0. + in_off_y = 0. + } + + # Decode the coordinate into a text string. + switch (WL_SYSTEM_TYPE(wd)) { + case RA_DEC: + if (axis == LONGITUDE) + call wl_hms (value, Memc[label], Memc[units], SZ_LINE, + precision, do_full) + else + call wl_dms (value, Memc[label], Memc[units], SZ_LINE, + precision, do_full) + default: + call sprintf (Memc[label], SZ_LINE, "%.2g") + call pargd (value) + } + + # Set the text justification. + call sprintf (Memc[label_format], SZ_LINE, "h=c;v=c;u=%f") + call pargr (text_angle) + call sprintf (Memc[units_format], SZ_LINE, "h=c;v=c;u=%f") + call pargr (text_angle) + + # Determine offset needed to rotate text about the point of placement. + # NOTE: The STDGRAPH kernel messes up rotate text placement. Try to + # accomodate with extra offset. + + length = .5 * char_width * (2 + strlen (Memc[label])) + call wl_rotate (length, 0., 1, text_angle - 90., rx, ry) + rx = abs (rx) + ry = abs (ry) + + # If labels are to appear inside the graph, then justification should + # appear as if it were done for the opposite side. + if (WL_LABOUT(wd) == YES) + tside = side + else + tside = wl_opposite_side (side) + + # Now add the offsets appropriately. + switch (tside) { + case TOP: + ly = y + ry + in_off_y + unit_off_y + if (fp_equalr (text_angle, 90.)) { + lx = x + ly = ly + unit_off_y + } else if (text_angle < 90.) + lx = x - rx + else + lx = x + rx + lx = lx + in_off_x + new_offset = ry + ry + + case BOTTOM: + ly = y - ry - in_off_y - unit_off_y + if (fp_equalr (text_angle, 90.)) { + lx = x + ly = ly - unit_off_y + } else if (text_angle < 90.) + lx = x + rx + else + lx = x - rx + lx = lx - in_off_x + new_offset = ry + ry + + case LEFT: + lx = x - rx - abs (unit_off_x) + if (text_angle < 90.) { + ly = y + ry - in_off_y + lx = lx - in_off_x + } else { + ly = y - ry + in_off_y + lx = lx + in_off_x + } + new_offset = rx + rx + abs (unit_off_x) + + case RIGHT: + lx = x + rx + abs (unit_off_x) + if (text_angle < 90.) { + ly = y - ry + in_off_y + lx = lx + in_off_x + } else { + ly = y + ry - in_off_y + lx = lx - in_off_x + } + new_offset = rx + rx + abs (unit_off_x) + } + + lx = lx - (unit_off_x / 2.) + ly = ly - (unit_off_y / 2.) + ux = lx + unit_off_x + uy = ly + unit_off_y + + # Print the label. + call gtext (WL_GP(wd), lx, ly, Memc[label], Memc[label_format]) + + # Print the units (if appropriate). + if (WL_SYSTEM_TYPE(wd) == RA_DEC) + call gtext (WL_GP(wd), ux, uy, Memc[units], Memc[units_format]) + + # Determine new maximum string size. + if ((WL_LABOUT(wd) == YES) && (abs (offset) < new_offset)) + if (side == LEFT || side == BOTTOM) + offset = -new_offset + else + offset = new_offset + + # Release memory. + call sfree (sp) +end + + +# WL_STRING_ANGLE -- Produce the angle that a label string should be written to. +# +# Description +# Fixes the input angle so that the output angle is in the range 0 to 180. +# +# Returns +# the angle that the label should be written as. + +double procedure wl_string_angle (angle, right_to_up) + +double angle # I: the input angle in degrees +int right_to_up # I: true if angle near horizontal/vertical are fixed + +double output_angle + +begin + # Try to ensure that the angle is "upright", i.e. the string will not + # be printed upside-down. + + output_angle = angle + if (output_angle > QUARTER_CIRCLE) + output_angle = output_angle - HALF_CIRCLE + if (output_angle < -QUARTER_CIRCLE) + output_angle = output_angle + HALF_CIRCLE + + # If the angle is close to parallel with one of the axis, then just + # print it normally. + + if ((right_to_up == YES) && ((mod (abs (output_angle), + QUARTER_CIRCLE) < MIN_ANGLE) || (QUARTER_CIRCLE - + mod (abs (output_angle), QUARTER_CIRCLE) < MIN_ANGLE))) + output_angle = 0. + + # Return the angle modified for the idiocincracy of GIO text angle + # specification. + + return (output_angle + QUARTER_CIRCLE) +end + + +# WL_ANGLE -- Return the average angle of the labels in the list. +# +# Returns +# Average angle +# +# Description +# So that labels on a side are uniform (in some sense), the average angle +# of all the labels is taken and is defined as the angle that all the labels +# will be printed at. + +double procedure wl_angle (wd, labels, nlabels) + +pointer wd # I: the WCSLAB descriptor +int labels[nlabels] # I: the indexes of the labels to be printed out +int nlabels # I: the number of indexes in the list + +double total, average +int i + +begin + total = 0.0 + for (i = 1; i <= nlabels; i = i + 1) + total = total + WL_LABEL_ANGLE(wd,labels[i]) + average = real (total / nlabels) + + return (average) +end diff --git a/pkg/utilities/nttools/stxtools/wcslab/wllabel.x.ori b/pkg/utilities/nttools/stxtools/wcslab/wllabel.x.ori new file mode 100644 index 00000000..33e86878 --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wllabel.x.ori @@ -0,0 +1,1077 @@ +include <gset.h> +include <math.h> +include "wcslab.h" +include "wcs_desc.h" + + +# Define the offset array. +define OFFSET Memr[$1+$2-1] + +# WL_LABEL -- Place the labels on the grids. +# +# Description +# Format and write the labels for the grid/tick marks. Much of this +# is wading through conditions to decide whether a label should be +# written or not. + +procedure wl_label (wd) + +pointer wd # I: the WCSLAB descriptor + +bool no_side_axis1, no_side_axis2 +int i, axis1_side, axis2_side +pointer sp, offset_ptr +real offset + +begin + # Get some memory. + call smark (sp) + call salloc (offset_ptr, N_SIDES, TY_REAL) + do i = 1, N_SIDES + OFFSET(offset_ptr,i) = 0. + + # Decide whether any sides were specified for either axis. + no_side_axis1 = true + no_side_axis2 = true + do i = 1, N_SIDES { + if (WL_LABEL_SIDE(wd,i,AXIS1)) + no_side_axis1 = false + if (WL_LABEL_SIDE(wd,i,AXIS2)) + no_side_axis2 = false + } + + # If polar, then label the axis 2's next to their circles on the + # graph and allow the Axis 1s to be labeled on all sides of the graph. + + if (WL_GRAPH_TYPE(wd) == POLAR) { + + call wl_polar_label (wd) + + if (no_side_axis1) { + do i = 1, N_SIDES { + WL_LABEL_SIDE(wd,i,AXIS1) = true + } + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(WD,AXIS1))) + WL_AXIS_TITLE_SIDE(WD,AXIS1) = BOTTOM + } + + # If we are near-polar, label the Axis 2 as if polar, and label + # Axis1 on all sides except the side closest to the pole. + + } else if (WL_GRAPH_TYPE(wd) == NEAR_POLAR) { + + if (no_side_axis1) { + WL_LABEL_SIDE(wd,WL_BAD_LABEL_SIDE(wd),AXIS1) = true + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(wd,AXIS1))) + WL_AXIS_TITLE_SIDE(wd,AXIS1) = WL_BAD_LABEL_SIDE(wd) + } + + if (no_side_axis2) { + WL_LABEL_SIDE(wd,WL_POLAR_LABEL_DIRECTION(wd),AXIS2) = true + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(wd,AXIS2))) + WL_AXIS_TITLE_SIDE(wd,AXIS2) = WL_POLAR_LABEL_DIRECTION(wd) + } + + # Final case- adjacent sides should be labelled. + + } else { + + # Determine the best sides for labelling. + if (INVERT (WL_ROTA(wd))) { + axis1_side = LEFT + axis2_side = BOTTOM + } else { + axis1_side = BOTTOM + axis2_side = LEFT + } + + # If no sides were specified, use the calculated ones above. + if (no_side_axis1) + WL_LABEL_SIDE(wd,axis1_side,AXIS1) = true + if (no_side_axis2) + WL_LABEL_SIDE(wd,axis2_side,AXIS2) = true + } + + # Now draw the labels for axis 1. + do i = 1, N_SIDES { + + if (WL_LABEL_SIDE(wd,i,AXIS1)) { + call wl_lab_edges (wd, AXIS1, i, offset) + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(WD,AXIS1))) + WL_AXIS_TITLE_SIDE(WD,AXIS1) = i + } else + offset = 0. + + # Modify the bounding box for the new viewport. + if (abs (offset) > abs (OFFSET(offset_ptr,i))) + OFFSET(offset_ptr,i) = offset + } + + # Draw the labels for axis 2. + if (WL_GRAPH_TYPE(wd) != POLAR) + do i = 1, N_SIDES { + + if (WL_LABEL_SIDE(wd,i,AXIS2)) { + call wl_lab_edges (wd, AXIS2, i, offset) + if (IS_INDEFI (WL_AXIS_TITLE_SIDE(wd,AXIS2))) + WL_AXIS_TITLE_SIDE(wd,AXIS2) = i + } else + offset = 0. + + # Modify the bounding box for the new viewport. + if (abs (offset) > abs (OFFSET(offset_ptr,i))) + OFFSET(offset_ptr,i) = offset + } + + # Set the bounding box. + do i = 1, N_SIDES + WL_NEW_VIEW(wd,i) = WL_NEW_VIEW(wd,i) + OFFSET(offset_ptr,i) + + # Now write the graph title. + call wl_title (WL_GP(wd), WL_AXIS_TITLE(wd,AXIS1), + WL_AXIS_TITLE_SIDE(wd,AXIS1), WL_AXIS_TITLE_SIZE(wd), + WL_NEW_VIEW(wd,1)) + if (WL_GRAPH_TYPE(wd) != POLAR) + call wl_title (WL_GP(wd), WL_AXIS_TITLE(wd,AXIS2), + WL_AXIS_TITLE_SIDE(wd,AXIS2), WL_AXIS_TITLE_SIZE(WD), + WL_NEW_VIEW(wd,1)) + if (! IS_INDEFI (WL_TITLE_SIDE(wd))) + call wl_title (WL_GP(wd), WL_TITLE(wd), WL_TITLE_SIDE(wd), + WL_TITLE_SIZE(wd), WL_NEW_VIEW(wd,1)) + + # Release memory. + call sfree (sp) +end + + +# Define what is in the screen. + +define IN (($1>WL_SCREEN_BOUNDARY(wd,LEFT))&&($1<WL_SCREEN_BOUNDARY(wd,RIGHT))&&($2>WL_SCREEN_BOUNDARY(wd,BOTTOM))&&($2<WL_SCREEN_BOUNDARY(wd,TOP))) + +# WL_POLAR_LABEL -- Place Latitude labels next to Latitude circles. +# +# Description +# Since Lines of constant Latitude on a polar graph are usually circles +# around the pole, the lines may never cross edges. Instead, the labels +# are placed next to circles. The grid-drawing routines should setup +# the label position array such that each line has only one label point. + +procedure wl_polar_label (wd) + +pointer wd # I: the WCSLAB descriptor + +int i, prec +pointer sp, label, units, label_format, units_format +real char_height, char_width, ndc_textx, ndc_texty, old_text_size +real textx, texty +int wl_precision() +real gstatr(), ggetr() + +begin + # Get some memory. + call smark (sp) + call salloc (label, SZ_LINE, TY_CHAR) + call salloc (units, SZ_LINE, TY_CHAR) + call salloc (label_format, SZ_LINE, TY_CHAR) + call salloc (units_format, SZ_LINE, TY_CHAR) + + # Get the character height and width. This is used to ensure that we + # have moved the label strings off the border. + + char_height = ggetr (WL_GP(wd), "ch") * gstatr (WL_GP(wd), G_TXSIZE) / + 2. + char_width = ggetr (WL_GP(wd), "cw") * gstatr (WL_GP(wd), G_TXSIZE) / + 2. + + # Get the text size and cut it in half for on the plot labelling. + old_text_size = gstatr (WL_GP(wd), G_TXSIZE) + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size) + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size * 0.80) + + # Determine the precision of the output. + prec = wl_precision (wd, AXIS2) + + # Place the labels. + for (i = 1; i <= WL_N_LABELS(wd); i = i + 1) + if (WL_LABEL_AXIS(wd,i) == AXIS2) { + + # Decode the coordinate into a text string. + call wl_dms (WL_LABEL_VALUE(wd,i), Memc[label], Memc[units], + SZ_LINE, prec, true) + + # Convert text position from "unknown" coordinates to NDC. + call gctran (WL_GP(wd), real (WL_LABEL_POSITION(wd,i,AXIS1)), + real (WL_LABEL_POSITION(wd,i,AXIS2)), ndc_textx, ndc_texty, + WL_PLOT_WCS(wd), WL_NDC_WCS(wd)) + + # Determine the text justification. + switch (WL_POLAR_LABEL_DIRECTION(wd)) { + case BOTTOM: + call strcpy ("h=c;v=t", Memc[label_format], SZ_LINE) + call strcpy ("h=c;v=c", Memc[units_format], SZ_LINE) + ndc_texty = ndc_texty - char_height + case TOP: + call strcpy ("h=c;v=c", Memc[label_format], SZ_LINE) + call strcpy ("h=c;v=b", Memc[units_format], SZ_LINE) + ndc_texty = ndc_texty + char_height + case LEFT: + call strcpy ("h=r;v=c", Memc[label_format], SZ_LINE) + call strcpy ("h=r;v=b", Memc[units_format], SZ_LINE) + ndc_textx = ndc_textx - char_width + case RIGHT: + call strcpy ("h=l;v=c", Memc[label_format], SZ_LINE) + call strcpy ("h=l;v=b", Memc[units_format], SZ_LINE) + ndc_textx = ndc_textx + char_width + } + + # Convert the text position from NDC back to the "unknown" + # system. + call gctran (WL_GP(wd), ndc_textx, ndc_texty, textx, texty, + WL_NDC_WCS(wd), WL_PLOT_WCS(wd)) + + # Print the label. + if (IN (textx, texty)) { + call gtext (WL_GP(wd), textx, texty, Memc[label], + Memc[label_format]) + call gtext (WL_GP(wd), textx, texty, Memc[units], + Memc[units_format]) + } + + } + + # Set the text size back. + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size) + + # Release memory. + call sfree (sp) + +end + + +# Memory management for labels + +define LABEL_LIST Memi[labels+$1-1] + +# WL_LAB_EDGES -- Place labels along the edges of the window. +# +# Description +# Place labels on the specified side of the graph. + +procedure wl_lab_edges (wd, axis, side, offset) + +pointer wd # I: the WCSLAB descriptor +int axis # I: the type of axis being labeled +int side # I: the side to place the labels +real offset # O: offset in NDC units for titles + +bool do_full +double angle, tangle +int i, full_label, nlabels, old_wcs, prec +pointer sp, labels +real ndc_textx, ndc_texty, old_text_size, textx, texty + +int wl_full_label_position(), wl_find_side() +double wl_string_angle(), wl_angle() +int gstati(), wl_precision() +real gstatr() + +begin + call smark (sp) + + # All label placement is done in NDC coordinates. + old_wcs = gstati (WL_GP(wd), G_WCS) + call gseti (WL_GP(wd), G_WCS, WL_NDC_WCS(wd)) + + # Set text labelling size. + old_text_size = gstatr (WL_GP(wd), G_TXSIZE) + call gsetr (WL_GP(wd), G_TXSIZE, WL_LABEL_SIZE(wd)) + + # Get the precision of the axis interval. + prec = wl_precision (wd, axis) + + # Initialize string size. + offset = 0. + + # Build a list of possible labels for this side. The conditions are + # that the label should be for the current axis and that it lies on + # the current side. + + call salloc (labels, WL_N_LABELS(wd), TY_INT) + nlabels = 0 + for (i = 1; i <= WL_N_LABELS(wd); i = i + 1) + if (WL_LABEL_AXIS(wd,i) == axis && + wl_find_side (WL_LABEL_POSITION(wd,i,AXIS1), + WL_LABEL_POSITION(wd,i,AXIS2), + WL_SCREEN_BOUNDARY(wd,1)) == side) { + nlabels = nlabels + 1 + LABEL_LIST(nlabels) = i + } + + # If no labels found, then just forget it. If labels found, well + # write them out. + + if (nlabels != 0) { + + # Determine which label should be written out in full. + full_label = wl_full_label_position (wd, Memi[labels], nlabels, + axis, side, prec) + + # Determine the angle that all the labels will be written at. + if ((WL_LABOUT(wd) == NO) && (WL_GRAPH_TYPE(wd) != NORMAL) && + (WL_LABEL_ROTATE(wd) == YES)) + angle = INDEFR + else if ((WL_GRAPH_TYPE(wd) == NORMAL) && ((WL_LABEL_ROTATE(wd) == + YES) || ((WL_LABOUT(wd) == NO) && (WL_MAJ_GRIDON(wd) == YES)))) + angle = wl_angle (wd, Memi[labels], nlabels) + else + angle = 0.0 + + # Place the labels. + for (i = 1; i <= nlabels; i = i + 1) { + + # Save some pertinent information. + textx = real (WL_LABEL_POSITION(wd,LABEL_LIST(i),AXIS1)) + texty = real (WL_LABEL_POSITION(wd,LABEL_LIST(i),AXIS2)) + do_full = ((LABEL_LIST(i) == full_label) || + (WL_ALWAYS_FULL_LABEL(wd) == YES)) + + # Transform the "unknown" coordinate system to a known + # coordinate system, NDC, for text placement. + call gctran (WL_GP(wd), textx, texty, ndc_textx, ndc_texty, + old_wcs, WL_NDC_WCS(wd)) + + # If angle is undefined, determine the angle for each label. + if (IS_INDEFR(angle)) + tangle = wl_string_angle (WL_LABEL_ANGLE(wd, + LABEL_LIST(i)), WL_LABOUT(wd)) + else + tangle = angle + + # Format and write the label. + call wl_write_label (wd, WL_LABEL_VALUE(wd,LABEL_LIST(i)), + side, ndc_textx, ndc_texty, tangle, axis, prec, do_full, + offset) + } + } + + # Reset the graphics WCS. + call gsetr (WL_GP(wd), G_TXSIZE, old_text_size) + call gseti (WL_GP(wd), G_WCS, old_wcs) + + call sfree (sp) +end + + +# WL_TITLE - Write the title of the graph. + +procedure wl_title (gp, title, side, size, viewport) + +pointer gp # I: the graphics descriptor +char title[ARB] # I: the title to write +int side # I: which side the title will go +real size # I: the character size to write the title +real viewport[N_SIDES] # I: the viewport in NDC to keep the title out of + +int old_wcs +real char_height, char_width, left, right, top, bottom, old_rotation +real old_text_size, x, y +int gstati(), strlen() +real ggetr(), gstatr() + +begin + # Make sure there is a title to write. If not, then punt. + if (strlen (title) <= 0) + return + + # Get/Set pertinent graphics info. + call ggview (gp, left, right, bottom, top) + + old_text_size = gstatr (gp, G_TXSIZE) + call gsetr (gp, G_TXSIZE, size) + old_rotation = gstatr (gp, G_TXUP) + + char_height = ggetr (gp, "ch") * size + char_width = ggetr (gp, "cw") * size + + old_wcs = gstati (gp, G_WCS) + call gseti (gp, G_WCS, NDC_WCS) + + # Depending on side, set text position and rotation. + switch (side) { + case TOP: + call gsetr (gp, G_TXUP, 90.) + x = (right + left) / 2. + y = viewport[TOP] + (2 * char_height) + viewport[TOP] = y + (char_height / 2.) + case BOTTOM: + call gsetr (gp, G_TXUP, 90.) + x = (right + left) / 2. + y = viewport[BOTTOM] - (2 * char_height) + viewport[BOTTOM] = y - (char_height / 2.) + case RIGHT: + call gsetr (gp, G_TXUP, 180.) + x = viewport[RIGHT] + (2 * char_width) + y = (top + bottom) / 2. + viewport[RIGHT] = x + (char_width / 2.) + case LEFT: + call gsetr (gp, G_TXUP, 180.) + x = viewport[LEFT] - (2 * char_width) + y = (top + bottom) / 2. + viewport[LEFT] = x - (char_width / 2.) + } + + # Write the puppy out. + call gtext (gp, x, y, title, "h=c;v=c") + + # Set the graphics state back. + call gseti (gp, G_WCS, old_wcs) + call gsetr (gp, G_TXSIZE, old_text_size) + call gsetr (gp, G_TXUP, old_rotation) +end + + +# WL_PRECISION -- Determine the precision of the interval. + +int procedure wl_precision (wd, axis) + +pointer wd # I: the WCSLAB descriptor +int axis # I: which axis is being examined ? + +int prec + +begin + # Handle the sky coordinates. + if (WL_SYSTEM_TYPE(wd) == RA_DEC) + + if (axis == AXIS1) { + if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (3600.0D0)) + prec = HOUR + else if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (60.0D0)) + prec = MINUTE + else if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (1.0D0)) + prec = SECOND + else if (WL_MAJOR_INTERVAL(wd,AXIS1) >= STTODEG (.01D0)) + prec = SUBSEC_LOW + else + prec = SUBSEC_HIGH + } else { + if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (3600.0D0)) + prec = DEGREE + else if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (60.0D0)) + prec = MINUTE + else if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (1.0D0)) + prec = SECOND + else if (WL_MAJOR_INTERVAL(wd,AXIS2) >= SATODEG (.01D0)) + prec = SUBSEC_LOW + else + prec = SUBSEC_HIGH + } + + # Handle other coordinate types. + else + prec = INDEFI + + return (prec) + +end + + +# Define some value constraints. + +define LOW_ACCURACY .01 +define HIGH_ACCURACY .0001 + +# WL_HMS -- Convert value to number in hours, minutes, and seconds. + +procedure wl_hms (rarad, hms, units, maxch, precision, all) + +double rarad # I: the value to format into a string (degrees) +char hms[ARB] # O: string containing formatted value +char units[ARB] # O: string containing formatted units +int maxch # I: the maximum number of characters allowed +int precision # I: how precise the output should be +bool all # I: true if all relevent fields should be formatted + +double accuracy, fraction +int sec, h, m, s +pointer sp, temp_hms, temp_units + +begin + # Get some memory. + call smark (sp) + call salloc (temp_hms, maxch, TY_CHAR) + call salloc (temp_units, maxch, TY_CHAR) + + units[1] = EOS + hms[1] = EOS + + # Define how close to zero is needed. + accuracy = LOW_ACCURACY + if (precision == SUBSEC_HIGH) + accuracy = HIGH_ACCURACY + + # Seconds of time. + fraction = double (abs(DEGTOST (rarad))) + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + sec = int (fraction) + fraction = fraction - double (sec) + } else { + sec = int (fraction + 0.5) + fraction = 0. + } + + # Range: 0 to 24 hours. + if (sec < 0) + sec = sec + STPERDAY + else if (sec >= STPERDAY) + sec = mod (sec, STPERDAY) + + # Separater fields. + s = mod (sec, 60) + m = mod (sec / 60, 60) + h = sec / 3600 + + # Format fields. + + # Subseconds. + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + fraction = s + fraction + if (precision == SUBSEC_LOW) { + call sprintf (hms, 6, "%05.2f") + call pargd (fraction) + call strcpy (" s ", units, maxch) + } else { + call sprintf (hms, 8, "%07.4f") + call pargd (fraction) + call strcpy (" s ", units, maxch) + } + if (!all) + all = (fraction < accuracy) + + # Seconds + } else if (precision == SECOND) { + + # NOTE: The all is not part of the if statement because if + # SUBSEC's have been printed, then seconds have already been + # dealt with. If SUBSEC's have not been dealt with, then this + # is the first field to be checked anyways. + + call sprintf (hms, 3, "%02d ") + call pargi (s) + call strcpy (" s", units, maxch) + if (! all) + all = (s == 0) + } + + # Minutes. + if (precision == MINUTE || (precision > MINUTE && all)) { + if (all) { + call strcpy (hms, Memc[temp_hms], maxch) + call strcpy (units, Memc[temp_units], maxch) + } + call sprintf (hms, 3, "%02d ") + call pargi (m) + call strcpy (" m", units, maxch) + if (all) { + call strcat (Memc[temp_hms], hms, maxch) + call strcat (Memc[temp_units], units, maxch) + } else + all = (m == 0) + } + + # Non-zero hours. + if (precision == HOUR || all) { + if (all) { + call strcpy (hms, Memc[temp_hms], maxch) + call strcpy (units, Memc[temp_units], maxch) + } + call sprintf (hms, 3, "%2.2d ") + call pargi (h) + call strcpy(" h", units, maxch) + if (all) { + call strcat (Memc[temp_hms], hms, maxch) + call strcat (Memc[temp_units], units, maxch) + } + } + + # Release memory + call sfree (sp) +end + + +# WL_DMS - Convert value to number in degrees, minutes, and seconds. + +procedure wl_dms (arcrad, dms, units, maxch, precision, all) + +double arcrad # I: the value to format into a string (degrees) +char dms[ARB] # O: string containing formatted value +char units[ARB] # O: string containing formatted units +int maxch # I: the maximum number of characters allowed +int precision # I: how precise the output should be ? +bool all # I: true if all relavent fields should be formatted + +double accuracy, fraction +int sec, h, m, s +pointer sp, temp_dms, temp_units +int strlen() + +begin + # Get some memory. + call smark (sp) + call salloc (temp_dms, maxch, TY_CHAR) + call salloc (temp_units, maxch, TY_CHAR) + + units[1] = EOS + dms[1] = EOS + + # Define how close to zero is needed. + accuracy = LOW_ACCURACY + if (precision == SUBSEC_HIGH) + accuracy = HIGH_ACCURACY + + # Seconds of time. + fraction = double (abs (DEGTOSA (arcrad))) + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + sec = int (fraction) + fraction = fraction - double (sec) + } else { + sec = nint (fraction) + fraction = 0. + } + + # Separater fields. + s = mod (abs(sec), 60) + m = mod (abs(sec) / 60, 60) + h = abs(sec) / 3600 + + # Format fields + # + # Subseconds. + if (precision == SUBSEC_LOW || precision == SUBSEC_HIGH) { + + fraction = s + fraction + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + if (precision == SUBSEC_LOW) { + call sprintf (dms, 6, "%05.2f\"") + call pargd (fraction) + call strcpy (" ", units, maxch) + } else { + call sprintf (dms, 8, "%07.4f\"") + call pargd (fraction) + call strcpy (" ", units, maxch) + } + if (! all) + all = (fraction < accuracy) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + + # Seconds + } else if (precision == SECOND) { + + # NOTE: The all is not part of the if statement because if + # SUBSEC's have been printed, then seconds have already been + # dealt with. If SUBSEC's have not been dealt with, then this + # is the first field to be checked anyways. + + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + call sprintf (dms, 3, "%02d\"") + call pargi (s) + call strcpy (" ", units, maxch) + if (! all) + all = (s == 0) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + } + + # Minutes. + if (precision == MINUTE || (precision > MINUTE && all)) { + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + call sprintf (dms, 3, "%02d'") + call pargi (m) + call strcpy (" ", units, maxch) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + if (! all) + all = (m == 0) + } + + # Hours. + if (precision == DEGREE || all) { + call strcpy (dms, Memc[temp_dms], maxch) + call strcpy (units, Memc[temp_units], maxch) + if (sec + fraction < accuracy) + call strcpy (" 0 ", dms, maxch) + else if (arcrad < 0.) { + call sprintf (dms, 4, "-%d ") + call pargi (h) + } else { + call sprintf (dms, 4, "+%d ") + call pargi (h) + } + call sprintf(units, 4, "%*wo") + call pargi (strlen (dms) - 1) + call strcat (Memc[temp_dms], dms, maxch) + call strcat (Memc[temp_units], units, maxch) + } + + # Release memory. + call sfree (sp) +end + + +# WL_FULL_LABEL_POSTION -- Find the position where the full label should be. +# +# Description +# This routine returns the index to the label that should be printed +# in its full form, regardless of its value. This is so there is always +# at least one labelled point with the full information. This point is +# choosen by examining which label is the closest to the passed point +# (usually one of the four corners of the display). +# +# Returns +# Index into the labell arrays of the label to be fully printed. +# If the return index is 0, then there are no labels for the given +# side. + +int procedure wl_full_label_position (wd, labels, nlabels, axis, side, + precision) + +pointer wd # I: the WCSLAB descriptor +int labels[nlabels] # I: array of indexes of labels to be printed +int nlabels # I: the number of labels in labels +int axis # I: the axis being dealt with +int side # I: the side being dealt with +int precision # I: precision of the label + +bool all +double cur_dist, dist +int i, cur_label, xside, yside +pointer sp, temp1 +double wl_distanced() + +begin + # Allocate some working space. + call smark (sp) + call salloc (temp1, SZ_LINE, TY_CHAR) + + # Initialize. + xside = INDEFI + yside = INDEFI + + # Determine which corner will have the full label. + if (side == TOP || side == BOTTOM) { + yside = side + if (axis == AXIS1) { + if (WL_LABEL_SIDE(wd,RIGHT,AXIS2)) + xside = RIGHT + if (WL_LABEL_SIDE(wd,LEFT,AXIS2)) + xside = LEFT + } else { + if (WL_LABEL_SIDE(wd,RIGHT,AXIS1)) + xside = RIGHT + if (WL_LABEL_SIDE(wd,LEFT,AXIS1)) + xside = LEFT + } + if (IS_INDEFI (xside)) + xside = LEFT + } else { + xside = side + if (axis == AXIS1) { + if (WL_LABEL_SIDE(wd,TOP,AXIS2)) + yside = TOP + if (WL_LABEL_SIDE(wd,BOTTOM,AXIS2)) + yside = BOTTOM + } else { + if (WL_LABEL_SIDE(wd,TOP,AXIS1)) + yside = TOP + if (WL_LABEL_SIDE(wd,BOTTOM,AXIS1)) + yside = BOTTOM + } + if (IS_INDEFI (yside)) + yside = BOTTOM + } + + # Find the full label. + cur_label = labels[1] + cur_dist = wl_distanced (WL_SCREEN_BOUNDARY(wd,xside), + WL_SCREEN_BOUNDARY(wd,yside), + WL_LABEL_POSITION(wd,cur_label,AXIS1), + WL_LABEL_POSITION(wd,cur_label,AXIS2)) + + # Now go through the rest of the labels to find a closer label. + for (i = 2; i <= nlabels; i = i + 1) { + + # Check to see if the label would be written in full anyways. + all = false + if (WL_SYSTEM_TYPE(wd) == RA_DEC) { + if (WL_LABEL_AXIS(wd, labels[i]) == LONGITUDE) + call wl_hms (WL_LABEL_VALUE(wd, labels[i]), + Memc[temp1], Memc[temp1], SZ_LINE, precision, all) + else + call wl_dms (WL_LABEL_VALUE(wd, labels[i]), + Memc[temp1], Memc[temp1], SZ_LINE, precision, all) + } + + # If so, don't figure out which label should be full, there + # will be one someplace. + if (all) { + cur_label = INDEFI + break + } + + dist = wl_distanced (WL_SCREEN_BOUNDARY(wd,xside), + WL_SCREEN_BOUNDARY(wd,yside), + WL_LABEL_POSITION(wd,labels[i],AXIS1), + WL_LABEL_POSITION(wd,labels[i],AXIS2)) + if (dist < cur_dist) { + cur_dist = dist + cur_label = labels[i] + } + } + + # Release memory. + call sfree (sp) + + # Return the label index. + return (cur_label) +end + + +# WL_WRITE_LABEL - Write the label in the format specified by the WCS type. + +procedure wl_write_label (wd, value, side, x, y, angle, axis, precision, + do_full, offset) + +pointer wd # I: the WCSLAB descriptor +double value # I: the value to use as the label +int side # I: the side the label is going on +real x, y # I: position of the label in NDC coordinates +double angle # I: the angle the text should be written at +int axis # I: which axis is being labelled +int precision # I: level of precision for labels +bool do_full # I: true if the full label should be printed +real offset # I/O: offset for titles in NDC units + +int tside +pointer sp, label, label_format, units, units_format +real char_height, char_width, in_off_x, in_off_y, length +real lx, ly, new_offset, rx, ry, text_angle +real unit_off_x, unit_off_y, ux, uy + +bool fp_equalr() +double wl_string_angle() +int wl_opposite_side(), strlen() +real ggetr(), gstatr() + +begin + # Get some memory. + call smark (sp) + call salloc (label, SZ_LINE, TY_CHAR) + call salloc (units, SZ_LINE, TY_CHAR) + call salloc (label_format, SZ_LINE, TY_CHAR) + call salloc (units_format, SZ_LINE, TY_CHAR) + + # Get character size. This info is used to move the character string + # by the appropriate amounts. + + char_height = ggetr (WL_GP(wd), "ch") * gstatr (WL_GP(wd), G_TXSIZE) + char_width = ggetr (WL_GP(wd), "cw") * gstatr (WL_GP(wd), G_TXSIZE) + + # Determine the "corrected" angle to write text in. + text_angle = wl_string_angle (angle, WL_LABOUT(wd)) + + # Determine the units offset. + call wl_rotate (0., char_height / 2., 1, text_angle - 90., unit_off_x, + unit_off_y) + + # If the labels are to appear inside the graph and the major grid lines + # have been drawn, then determine the necessary offset to get the label + # off the line. + + if ((WL_LABOUT(wd) == NO) && (WL_MAJ_GRIDON(wd) == YES)) + call wl_rotate (0., 0.75 * char_height, 1, text_angle - 90., + in_off_x, in_off_y) + else { + in_off_x = 0. + in_off_y = 0. + } + + # Decode the coordinate into a text string. + switch (WL_SYSTEM_TYPE(wd)) { + case RA_DEC: + if (axis == LONGITUDE) + call wl_hms (value, Memc[label], Memc[units], SZ_LINE, + precision, do_full) + else + call wl_dms (value, Memc[label], Memc[units], SZ_LINE, + precision, do_full) + default: + call sprintf (Memc[label], SZ_LINE, "%.2g") + call pargd (value) + } + + # Set the text justification. + call sprintf (Memc[label_format], SZ_LINE, "h=c;v=c;u=%f") + call pargr (text_angle) + call sprintf (Memc[units_format], SZ_LINE, "h=c;v=c;u=%f") + call pargr (text_angle) + + # Determine offset needed to rotate text about the point of placement. + # NOTE: The STDGRAPH kernel messes up rotate text placement. Try to + # accomodate with extra offset. + + length = .5 * char_width * (2 + strlen (Memc[label])) + call wl_rotate (length, 0., 1, text_angle - 90., rx, ry) + rx = abs (rx) + ry = abs (ry) + + # If labels are to appear inside the graph, then justification should + # appear as if it were done for the opposite side. + if (WL_LABOUT(wd) == YES) + tside = side + else + tside = wl_opposite_side (side) + + # Now add the offsets appropriately. + switch (tside) { + case TOP: + ly = y + ry + in_off_y + unit_off_y + if (fp_equalr (text_angle, 90.)) { + lx = x + ly = ly + unit_off_y + } else if (text_angle < 90.) + lx = x - rx + else + lx = x + rx + lx = lx + in_off_x + new_offset = ry + ry + + case BOTTOM: + ly = y - ry - in_off_y - unit_off_y + if (fp_equalr (text_angle, 90.)) { + lx = x + ly = ly - unit_off_y + } else if (text_angle < 90.) + lx = x + rx + else + lx = x - rx + lx = lx - in_off_x + new_offset = ry + ry + + case LEFT: + lx = x - rx - abs (unit_off_x) + if (text_angle < 90.) { + ly = y + ry - in_off_y + lx = lx - in_off_x + } else { + ly = y - ry + in_off_y + lx = lx + in_off_x + } + new_offset = rx + rx + abs (unit_off_x) + + case RIGHT: + lx = x + rx + abs (unit_off_x) + if (text_angle < 90.) { + ly = y - ry + in_off_y + lx = lx + in_off_x + } else { + ly = y + ry - in_off_y + lx = lx - in_off_x + } + new_offset = rx + rx + abs (unit_off_x) + } + + lx = lx - (unit_off_x / 2.) + ly = ly - (unit_off_y / 2.) + ux = lx + unit_off_x + uy = ly + unit_off_y + + # Print the label. + call gtext (WL_GP(wd), lx, ly, Memc[label], Memc[label_format]) + + # Print the units (if appropriate). + if (WL_SYSTEM_TYPE(wd) == RA_DEC) + call gtext (WL_GP(wd), ux, uy, Memc[units], Memc[units_format]) + + # Determine new maximum string size. + if ((WL_LABOUT(wd) == YES) && (abs (offset) < new_offset)) + if (side == LEFT || side == BOTTOM) + offset = -new_offset + else + offset = new_offset + + # Release memory. + call sfree (sp) +end + + +# WL_STRING_ANGLE -- Produce the angle that a label string should be written to. +# +# Description +# Fixes the input angle so that the output angle is in the range 0 to 180. +# +# Returns +# the angle that the label should be written as. + +double procedure wl_string_angle (angle, right_to_up) + +double angle # I: the input angle in degrees +int right_to_up # I: true if angle near horizontal/vertical are fixed + +double output_angle + +begin + # Try to ensure that the angle is "upright", i.e. the string will not + # be printed upside-down. + + output_angle = angle + if (output_angle > QUARTER_CIRCLE) + output_angle = output_angle - HALF_CIRCLE + if (output_angle < -QUARTER_CIRCLE) + output_angle = output_angle + HALF_CIRCLE + + # If the angle is close to parallel with one of the axis, then just + # print it normally. + + if ((right_to_up == YES) && ((mod (abs (output_angle), + QUARTER_CIRCLE) < MIN_ANGLE) || (QUARTER_CIRCLE - + mod (abs (output_angle), QUARTER_CIRCLE) < MIN_ANGLE))) + output_angle = 0. + + # Return the angle modified for the idiocincracy of GIO text angle + # specification. + + return (output_angle + QUARTER_CIRCLE) +end + + +# WL_ANGLE -- Return the average angle of the labels in the list. +# +# Returns +# Average angle +# +# Description +# So that labels on a side are uniform (in some sense), the average angle +# of all the labels is taken and is defined as the angle that all the labels +# will be printed at. + +double procedure wl_angle (wd, labels, nlabels) + +pointer wd # I: the WCSLAB descriptor +int labels[nlabels] # I: the indexes of the labels to be printed out +int nlabels # I: the number of indexes in the list + +double total, average +int i + +begin + total = 0.0 + for (i = 1; i <= nlabels; i = i + 1) + total = total + WL_LABEL_ANGLE(wd,labels[i]) + average = real (total / nlabels) + + return (average) +end diff --git a/pkg/utilities/nttools/stxtools/wcslab/wlsetup.x b/pkg/utilities/nttools/stxtools/wcslab/wlsetup.x new file mode 100644 index 00000000..c37e24ca --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wlsetup.x @@ -0,0 +1,1000 @@ +include <gset.h> +include <mach.h> +include <math.h> +include <math/curfit.h> +include "wcslab.h" +include "wcs_desc.h" + +# WL_SETUP -- Determine all the basic characteristics of the plot. +# +# Description +# Determine basic characteristics of the plot at hand. This involved +# "discovering" what part of the world system covers the screen, the +# orientation of the world to logical systems, what type of graph will +# be produced, etc. Many of the parameters determined here can be +# over-ridden by user-specified values. + +procedure wl_setup (wd) + +pointer wd # I: the WCSLAB descriptor + +bool north +double array[N_EDGES,N_DIM], max_value[N_DIM], min_value[N_DIM] +double range[N_DIM], pole_position[N_DIM], view_edge[N_EDGES,N_DIM] +double wl_coord_rotation() +pointer mw_sctran() +string logtran "logical" +string wrldtran "world" + +begin + # Calculate the transformations from the Logical (pixel space) system + # to the World (possibly anything) system and back. + WL_LWCT(wd) = mw_sctran (WL_MW(wd), logtran, wrldtran, AXIS) + WL_WLCT(wd) = mw_sctran (WL_MW(wd), wrldtran, logtran, AXIS) + + # Indicate whether the center of the transformation is north. + if (WL_SYSTEM_TYPE(wd) == RA_DEC) + north = (WL_WORLD_CENTER(wd,LATITUDE) > 0.0D0) + + # Determine the poles position. + if (WL_SYSTEM_TYPE(wd) == RA_DEC) + call wl_pole_position (WL_WLCT(wd), WL_AXIS_FLIP(wd), + WL_WORLD_CENTER(wd,LONGITUDE), north, WL_SYSTEM_TYPE(wd), + pole_position) + + # Determine graph type based on the system type. + call wl_determine_graph_type (WL_SYSTEM_TYPE(wd), pole_position, + WL_SCREEN_BOUNDARY(wd,1), WL_GRAPH_TYPE(wd)) + + # Now find the extent of the WCS the window views, by constructing + # x,y vectors containing evenly spaced points around the edges of + # the viewing window. + + call wl_construct_edge_vectors (WL_SCREEN_BOUNDARY(wd,1), + view_edge[1,X_DIM], view_edge[1,Y_DIM], N_EDGES) + + # Find the range of the axes over the graphics viewport. + call wl_l2wd (WL_LWCT(wd), WL_AXIS_FLIP(wd), view_edge[1,X_DIM], + view_edge[1,Y_DIM], array[1,AXIS1], array[1,AXIS2], N_EDGES) + call alimd (array[1,AXIS1], N_EDGES, min_value[AXIS1], max_value[AXIS1]) + call alimd (array[1,AXIS2], N_EDGES, min_value[AXIS2], max_value[AXIS2]) + range[AXIS1] = abs (max_value[AXIS1] - min_value[AXIS1]) + range[AXIS2] = abs (max_value[AXIS2] - min_value[AXIS2]) + + # The above isn't good enough for the sky projections. Deal with those. + if (WL_SYSTEM_TYPE(wd) == RA_DEC) + call wl_sky_extrema (wd, array[1,AXIS1], N_EDGES, pole_position, + north, min_value[AXIS1], max_value[AXIS1], range[AXIS1], + min_value[AXIS2], max_value[AXIS2], range[AXIS2]) + + # Determine the rotation between the systems. + WL_ROTA(wd) = wl_coord_rotation (WL_WLCT(wd), WL_AXIS_FLIP(wd), + WL_WORLD_CENTER(wd,AXIS1), max_value[AXIS2], + WL_WORLD_CENTER(wd,AXIS1), min_value[AXIS2]) + + # Round the intervals. This is done to make the labelling "nice" and + # to smooth edge effects. + if (IS_INDEFD (WL_MAJOR_INTERVAL(wd,AXIS1)) || + IS_INDEFD (WL_BEGIN(wd,AXIS1)) || IS_INDEFD (WL_END(wd,AXIS1))) + call wl_round_axis (wd, AXIS1, min_value[AXIS1], max_value[AXIS1], + range[AXIS1]) + + if (IS_INDEFD (WL_MAJOR_INTERVAL(wd,AXIS2)) || + IS_INDEFD (WL_BEGIN(wd,AXIS2)) || IS_INDEFD (WL_END(wd,AXIS2))) + call wl_round_axis (wd, AXIS2, min_value[AXIS2], max_value[AXIS2], + range[AXIS2]) +end + + +# WL_POLE_POSITION -- Determine logical coordinates of a pole. +# +# Description +# Calculate the pole's position in the Logical system. +# +# Bugs +# Can only deal with Right Ascension/Declination. + +procedure wl_pole_position (wlct, flip, longitude, north, system_type, + pole_position) + +pointer wlct # I: the world-to-logical transformation +int flip # I: true if the axes are transposed +double longitude # I: the longitude to determine latitude +bool north # I: true if the pole is in the north +int system_type # I: type of system being examined +double pole_position[N_DIM] # O: the pole's logical coordinates + +double sgn + +begin + switch (system_type) { + + # For Right Ascension/Declination, the pole is at any longitude but + # at only 90 degrees (north) or -90 degrees (south) latitude. + case RA_DEC: + if (north) + sgn = NORTH_POLE_LATITUDE + else + sgn = SOUTH_POLE_LATITUDE + call wl_w2ld (wlct, flip, longitude, sgn, pole_position[X_DIM], + pole_position[Y_DIM], 1) + } + + # Sanity check on the pole position. It is very likely that there is + # no valid position in pixel space for the pole. This is checked for + # by looking for extremely large numbers. + if (abs (pole_position[X_DIM]) > abs (double (MAX_INT))) + pole_position[X_DIM] = real (MAX_INT) + if (abs (pole_position[Y_DIM]) > abs (double (MAX_INT))) + pole_position[Y_DIM] = real (MAX_INT) +end + + +# How close can the pole be to the center of the screen to be near-polar. +define HOW_CLOSE 3. + +# WL_DETERMINE_GRAPH_TYPE -- Determine the actual graph type. + +procedure wl_determine_graph_type (system_type, pole_position, + screen_boundary, graph_type) + +int system_type # I: the type of WCS being dealt with +double pole_position[N_DIM] # I: the location of the pole +double screen_boundary[N_SIDES] # I: the edges of the display +int graph_type # O: the graph type + +double max_dist, pole_dist, xcen, ycen + +begin + # Determine graph type based on axis type. + switch (system_type) { + + # If the pole is on the graph then force a graph_type of polar. + case RA_DEC: + + xcen = (screen_boundary[LEFT] + screen_boundary[RIGHT]) / 2. + ycen = (screen_boundary[BOTTOM] + screen_boundary[TOP]) / 2. + max_dist = min ((screen_boundary[LEFT] - xcen) ** 2, + (screen_boundary[TOP] - ycen)**2) + pole_dist = (pole_position[X_DIM] - xcen) ** 2 + + (pole_position[Y_DIM] - ycen) ** 2 + + # Check to see whether the graph is "polar", "near_polar" + # or "normal". If the pole lies within middle part of the + # viewport, then the graph is "polar". If the pole is within + # a certain maximum distance then it is "near_polar". + # Otherwise it is normal. + + switch (graph_type) { + case NORMAL: + # do nothing + case POLAR: + # do nothing + case NEAR_POLAR: + # do nothing + default: + if (pole_dist < max_dist) + graph_type = POLAR + else if (pole_dist < HOW_CLOSE * max_dist) + graph_type = NEAR_POLAR + else + graph_type = NORMAL + } + + # For all other cases, explicitely set this to normal. + default: + graph_type = NORMAL + } +end + + +# WL_CONSTRUCT_EDGE_VECTORS -- Construct vectors of values along window's edge. +# +# Description +# This routines filles two arrays, with the x-values and y-values of +# evenly spaced points along the edges of the screen. This is used to +# make transformation of the logical edges into the world system +# more convenient. + +procedure wl_construct_edge_vectors (screen_boundary, x, y, vector_size) + +double screen_boundary[N_SIDES] # I: the side values +double x[vector_size], y[vector_size] # O: the edge vector points +int vector_size # I: the number of edge vector points + +double current, interval +int i, left_over, offset1, offset2, side_length + +begin + # Divide the vectors into equal amounts for each side. + side_length = vector_size / N_SIDES + left_over = mod (vector_size, N_SIDES) + + # Calculate the horizontal components. + interval = (screen_boundary[RIGHT] - screen_boundary[LEFT]) / + side_length + current = screen_boundary[LEFT] + offset1 = side_length + for (i = 1; i <= side_length; i = i + 1) { + x[i] = current + interval + y[i] = screen_boundary[BOTTOM] + x[i+offset1] = current + y[i+offset1] = screen_boundary[TOP] + current = current + interval + } + + # Calculate the verticle components. + interval = (screen_boundary[TOP] - screen_boundary[BOTTOM]) / + side_length + current = screen_boundary[BOTTOM] + offset1 = 2 * side_length + offset2 = 3 * side_length + for (i = 1; i <= side_length; i = i + 1) { + x[i+offset1] = screen_boundary[LEFT] + y[i+offset1] = current + x[i+offset2] = screen_boundary[RIGHT] + y[i+offset2] = current + interval + current = current + interval + } + + # Fill in the left over with a single point. + offset1 = 4 * side_length + for (i = 1; i <= left_over; i = i + 1) { + x[i+offset1] = screen_boundary[LEFT] + y[i+offset1] = screen_boundary[BOTTOM] + } + +end + + +# WL_SKY_EXTREMA -- Determine what range the view window covers in the sky. +# This routine is only called if the WCS RA,DEC. +# +# Description +# Because of the different graph types and the fact that axis 1 usually +# wraps, more work needs to be done to determine what part of the sky +# is covered by the viewing window. + +procedure wl_sky_extrema (wd, ax1_array, n_points, pole_position, north, + ax1min, ax1max, ax1ran, ax2min, ax2max, ax2ran) + +pointer wd # I: the WCSLAB descriptor +double ax1_array[n_points] # I: the axis 1 edge vector +int n_points # I: the length of the edge vector +double pole_position[N_DIM] # I: the pole position +bool north # I: is the pole in the north ? +double ax1min, ax1max, ax1ran # I/O: the minimum, maximum, range in axis 1 +double ax2min, ax2max, ax2ran # I/O: the minimum, maximum, range in axis 2 + +bool is_pole +double nx, ny, xcen, ycen +int wl_direction_from_axis1(), wl_find_side(), wl_opposite_side() + +begin + # Is the pole on the graph ? + if ((pole_position[X_DIM] < WL_SCREEN_BOUNDARY(wd,LEFT)) || + (pole_position[X_DIM] > WL_SCREEN_BOUNDARY(wd,RIGHT)) || + (pole_position[Y_DIM] < WL_SCREEN_BOUNDARY(wd,BOTTOM)) || + (pole_position[Y_DIM] > WL_SCREEN_BOUNDARY(wd,TOP))) + is_pole = false + else + is_pole = true + + # If so adjust the RA and DEC ranges appropriately. + if (is_pole) { + + # Set the RA range. + ax1min = 0.0D0 + ax1max = 359.9D0 + ax1ran = 360.0D0 + + # Set the dec range. + if (north) + ax2max = NORTH_POLE_LATITUDE - ((NORTH_POLE_LATITUDE - + ax2min) * DISTANCE_TO_POLE ) + else + ax2min = SOUTH_POLE_LATITUDE + ((NORTH_POLE_LATITUDE + + ax2max) * DISTANCE_TO_POLE) + ax2ran = abs (ax2max - ax2min) + + # Mark the pole. + call gmark (WL_GP(wd), real (pole_position[X_DIM]), + real (pole_position[Y_DIM]), POLE_MARK_SHAPE, POLE_MARK_SIZE, + POLE_MARK_SIZE) + + } else { + # Only the RA range needs adjusting. + call wl_ra_range (ax1_array, n_points, ax1min, ax1max, ax1ran) + } + + # Adjust the labelling characteristics appropritatley for various + # types of graphs. + + if (WL_GRAPH_TYPE(wd) == POLAR) { + + # Determine which direction the axis 2's will be labeled on polar + # graphs. + if (IS_INDEFD (WL_POLAR_LABEL_POSITION(wd))) { + call wl_get_axis2_label_direction (WL_LWCT(wd), + WL_AXIS_FLIP(wd), pole_position, WL_SCREEN_BOUNDARY(wd,1), + WL_POLAR_LABEL_POSITION(wd), WL_BAD_LABEL_SIDE(wd)) + } else { + WL_BAD_LABEL_SIDE(wd) = wl_direction_from_axis1 (WL_WLCT(wd), + WL_AXIS_FLIP(wd), pole_position, north, + WL_POLAR_LABEL_POSITION(wd), WL_BEGIN(wd,AXIS2), + WL_END(wd,AXIS2), WL_SCREEN_BOUNDARY(wd,1)) + if (IS_INDEFI (WL_BAD_LABEL_SIDE(wd))) + WL_BAD_LABEL_SIDE(wd) = BOTTOM + } + + # If the graph type is polar, then determine how to justify + # the labels. + + if (IS_INDEFI (WL_POLAR_LABEL_DIRECTION(wd))) + WL_POLAR_LABEL_DIRECTION(wd) = + wl_opposite_side (WL_BAD_LABEL_SIDE(wd)) + + # If the graph_type is near-polar, then handle the directions a bit + # differently. + } else if (WL_GRAPH_TYPE(wd) == NEAR_POLAR) { + + # Find the side that the pole is on. + xcen = (WL_SCREEN_BOUNDARY(wd,LEFT) + + WL_SCREEN_BOUNDARY(wd,RIGHT)) / 2. + ycen = (WL_SCREEN_BOUNDARY(wd,BOTTOM) + + WL_SCREEN_BOUNDARY(wd,TOP)) / 2. + call wl_axis_on_line (xcen, ycen, pole_position[X_DIM], + pole_position[Y_DIM], WL_SCREEN_BOUNDARY(wd,1), nx, ny) + + if (IS_INDEFD(nx) || IS_INDEFD(ny)) { + WL_BAD_LABEL_SIDE(wd) = BOTTOM + WL_POLAR_LABEL_DIRECTION(wd) = LEFT + } else { + WL_BAD_LABEL_SIDE(wd) = wl_find_side (nx, ny, + WL_SCREEN_BOUNDARY(wd,1)) + if (WL_BAD_LABEL_SIDE(wd) == LEFT || WL_BAD_LABEL_SIDE(wd) == + RIGHT) + if (abs (ny - WL_SCREEN_BOUNDARY(wd,BOTTOM)) < + abs (ny - WL_SCREEN_BOUNDARY(wd,TOP))) + WL_POLAR_LABEL_DIRECTION(wd) = BOTTOM + else + WL_POLAR_LABEL_DIRECTION(wd) = TOP + else + if (abs (nx - WL_SCREEN_BOUNDARY(wd,LEFT)) < + abs (nx - WL_SCREEN_BOUNDARY(wd,RIGHT))) + WL_POLAR_LABEL_DIRECTION(wd) = LEFT + else + WL_POLAR_LABEL_DIRECTION(wd) = RIGHT + } + + } +end + + +# WL_COORD_ROTATION -- Determine "rotation" between the coordinate systems. +# +# Description +# This routine takes the world-to-logical coordinate transformation and +# two points in the world system which should define the positive verticle +# axis in the world system. These points are translated into the logical +# system and the angle between the logical vector and its positive verticle +# vector is calculated and returned. The rotation angle is returned +# in degrees and is always positive. + +double procedure wl_coord_rotation (wlct, flip, wx1, wy1, wx2, wy2) + +pointer wlct # I: the world-to-logical transformation +int flip # I: true if the coordinates are transposed +double wx1, wy1, wx2, wy2 # I: points in world space to figure rotation from + +double delx, dely, rota, x1, y1, x2, y2 +bool fp_equald() + +begin + # Transform the points to the logical system. + call wl_w2ld (wlct, flip, wx1, wy1, x1, y1, 1) + call wl_w2ld (wlct, flip, wx2, wy2, x2, y2, 1) + + # Determine the rotation. + delx = x2 - x1 + dely = y2 - y1 + if (fp_equald (delx, 0.0D0) && fp_equald (dely, 0.0D0)) + rota = 0. + else + rota = RADTODEG (atan2 (dely, delx)) + + if (rota < 0.0D0) + rota = rota + FULL_CIRCLE + + return (rota) +end + + +# Define how many axis one should go for. + +define RA_NUM_TRY 6 +define DEC_NUM_TRY 6 +define DEC_POLAR_NUM_TRY 4 + +# WL_ROUND_AXIS - Round values for the axis. + +procedure wl_round_axis (wd, axis, minimum, maximum, range) + +pointer wd # I: the WCSLAB descriptor +int axis # I: the axis being worked on +double minimum, maximum, range # I: raw values to be rounded + +int num_try + +begin + # Depending on axis type, round the values. + switch (WL_SYSTEM_TYPE(wd)) { + case RA_DEC: + if (axis == LONGITUDE) + call wl_round_ra (minimum, maximum, range, RA_NUM_TRY, + WL_BEGIN(wd,LONGITUDE), WL_END(wd,LONGITUDE), + WL_MAJOR_INTERVAL(wd,LONGITUDE)) + else { + if (WL_GRAPH_TYPE(wd) == POLAR) + num_try = DEC_POLAR_NUM_TRY + else + num_try = DEC_NUM_TRY + call wl_round_dec (minimum, maximum, range, num_try, + WL_BEGIN(wd,LATITUDE), WL_END(wd,LATITUDE), + WL_MAJOR_INTERVAL(wd,LATITUDE)) + } + + default: + call wl_generic_round (minimum, maximum, range, WL_BEGIN(wd,axis), + WL_END(wd,axis), WL_MAJOR_INTERVAL(wd,axis)) + } + +end + + +# WL_GET_AXIS2_LABEL_DIRECTION -- Dertermine label direction for latitides. +# +# Description +# Determine from which edge of the graph the axis 2 labels are to +# appear. This (in general) is the opposite edge from which the pole +# is nearest to. Move the pole to the closest edges, determine which +# side it is, then chose the direction as the opposite. Also determines +# the Axis 1 at which the Axis 2 labels will appear. + +procedure wl_get_axis2_label_direction (lwct, flip, pole_position, + screen_boundary, pole_label_position, bad_label_side) + +pointer lwct # I: logical-to-world transformation +int flip # I: true if the axis are transposed +double pole_position[N_DIM] # I: the position of the pole +double screen_boundary[N_SIDES] # I: the edges of the screen +double pole_label_position # O: the axis 1 that axis 2 labels should + # appear for polar|near-polar graphs +int bad_label_side # O: side not to place axis 1 labels + +double dif, tdif, dummy + +begin + # Determine which direction, up or down, the axis 2's will be labelled. + dif = abs (screen_boundary[TOP] - pole_position[AXIS2]) + bad_label_side= TOP + tdif = abs (screen_boundary[BOTTOM] - pole_position[AXIS2]) + if (tdif < dif) { + dif = tdif + bad_label_side = BOTTOM + } + + # Determine at what value of Axis 1 the Axis 2 labels should appear. + switch (bad_label_side) { + case TOP: + call wl_l2wd (lwct, flip, pole_position[AXIS1], + screen_boundary[BOTTOM], pole_label_position, dummy, 1) + case BOTTOM: + call wl_l2wd (lwct, flip, pole_position[AXIS1], + screen_boundary[TOP], pole_label_position, dummy, 1) + case LEFT: + call wl_l2wd (lwct, flip, screen_boundary[RIGHT], + pole_position[AXIS2], pole_label_position, dummy, 1) + case RIGHT: + call wl_l2wd (lwct, flip, screen_boundary[LEFT], + pole_position[AXIS2], pole_label_position, dummy, 1) + } + +end + + +# WL_DIRECTION_FROM_AXIS1 -- Determine axis 2 label direction from axis 1. +# +# Function Returns +# This returns the side where Axis 1 should not be labelled. + +int procedure wl_direction_from_axis1 (wlct, flip, pole_position, north, + polar_label_position, lbegin, lend, screen_boundary) + +pointer wlct # I: world-to-logical transformation +int flip # I: true if the axes are transposed +double pole_position[N_DIM] # I: the pole position +bool north # I: true if the pole is the north pole +double polar_label_position # I: the axis 1 where axis 2 will be + # marked +double lbegin # I: low end of axis 2 +double lend # I: high end of axis 2 +double screen_boundary[N_SIDES] # I: the window boundary + +double nx, ny, cx, cy +int wl_find_side() + +begin + # Determine the point in logical space where the axis 1 and the + # minimum axis 2 meet. + + if (north) + call wl_w2ld (wlct, flip, polar_label_position, lbegin, nx, ny, 1) + else + call wl_w2ld (wlct, flip, polar_label_position, lend, nx, ny, 1) + + # This line should cross a window boundary. Find that point. + + call wl_axis_on_line (pole_position[X_DIM], pole_position[Y_DIM], + screen_boundary, nx, ny, cx, cy) + + # Get the side that the crossing point is. This is the axis 2 labelling + # direction. + + if (IS_INDEFD(cx) || IS_INDEFD(cy)) + return (INDEFI) + else + return (wl_find_side (cx, cy, screen_boundary)) +end + + +# WL_OPPOSITE_SIDE - Return the opposite of the given side. +# +# Returns +# The opposite side of the specified side as follows: +# RIGHT -> LEFT +# LEFT -> RIGHT +# TOP -> BOTTOM +# BOTTOM -> TOP + +int procedure wl_opposite_side (side) + +int side # I: the side to find the opposite of + +int new_side + +begin + switch (side) { + case LEFT: + new_side = RIGHT + case RIGHT: + new_side = LEFT + case TOP: + new_side = BOTTOM + case BOTTOM: + new_side = TOP + } + + return (new_side) +end + + +# Define whether things are on the screen boundary or on them. + +define IN (($1>=screen_boundary[LEFT])&&($1<=screen_boundary[RIGHT])&&($2>=screen_boundary[BOTTOM])&&($2<=screen_boundary[TOP])) + + +# WL_AXIS_ON_LINE - Determine intersection of line and a screen boundary. +# +# Description +# Return the point where the line defined by the two input points +# crosses a screen boundary. The boundary is choosen by determining +# which one is between the two points. + +procedure wl_axis_on_line (x0, y0, x1, y1, screen_boundary, nx, ny) + +double x0, y0, x1, y1 # I: random points in space +double screen_boundary[N_SIDES] # I: sides of the window +double nx, ny # O: the closest point on a window boundary + +double x_val[N_SIDES], y_val[N_SIDES], tx0, ty0, tx1, ty1, w[2] +int i +pointer cvx, cvy +double dcveval() + +begin + # Get the line parameters. + x_val[1] = x0 + x_val[2] = x1 + y_val[1] = y0 + y_val[2] = y1 + + iferr (call dcvinit (cvx, CHEBYSHEV, 2, min (x0, x1), max (x0, x1))) + cvx = NULL + else { + call dcvfit (cvx, x_val, y_val, w, 2, WTS_UNIFORM, i) + if (i != OK) + call error (i, "wlaxie: Error solving on X") + } + + iferr (call dcvinit (cvy, CHEBYSHEV, 2, min (y0, y1), max (y0, y1))) + cvy = NULL + else { + call dcvfit (cvy, y_val, x_val, w, 2, WTS_UNIFORM, i) + if (i != OK) + call error (i, "wlaxie: Error solving on Y") + } + + # Solve for each side. + x_val[LEFT] = screen_boundary[LEFT] + if (cvx == NULL) + y_val[LEFT] = screen_boundary[LEFT] + else + y_val[LEFT] = dcveval (cvx, x_val[LEFT]) + + x_val[RIGHT] = screen_boundary[RIGHT] + if (cvx == NULL ) + y_val[RIGHT] = screen_boundary[RIGHT] + else + y_val[RIGHT] = dcveval (cvx, x_val[RIGHT]) + + y_val[TOP] = screen_boundary[TOP] + if (cvy == NULL) + x_val[TOP] = screen_boundary[TOP] + else + x_val[TOP] = dcveval (cvy, y_val[TOP]) + + y_val[BOTTOM] = screen_boundary[BOTTOM] + if (cvy == NULL) + x_val[BOTTOM] = screen_boundary[BOTTOM] + else + x_val[BOTTOM] = dcveval (cvy, y_val[BOTTOM]) + + # Rearrange the input points to be in ascending order. + if (x0 < x1) { + tx0 = x0 + tx1 = x1 + } else { + tx0 = x1 + tx1 = x0 + } + + if (y0 < y1) { + ty0 = y0 + ty1 = y1 + } else { + ty0 = y1 + ty1 = y0 + } + + # Now find which point is between the two given points and is within + # the viewing area. + # NOTE: Conversion to real for the check- if two points are so close + # for double, any of them would serve as the correct answer. + + nx = INDEFD + ny = INDEFD + for (i = 1; i <= N_SIDES; i = i + 1) + if (real (tx0) <= real (x_val[i]) && + real (x_val[i]) <= real (tx1) && + real (ty0) <= real (y_val[i]) && + real (y_val[i]) <= real (ty1) && + IN (x_val[i], y_val[i]) ) { + nx = x_val[i] + ny = y_val[i] + } + + # Release the curve fit descriptors. + if (cvx != NULL) + call dcvfree (cvx) + if (cvy != NULL) + call dcvfree (cvy) +end + + +# WL_FIND_SIDE -- Return the side that the given point is lying on. +# +# Function Returns +# Return the side, TOP, BOTTOM, LEFT, or RIGHT, that the specified +# point is lying on. One of the coordinates must be VERY CLOSE to one of +# the sides or INDEFI will be returned. + +int procedure wl_find_side (x, y, screen_boundary) + +double x, y # I: the point to inquire about +double screen_boundary[N_SIDES] # I: the edges of the screen + +double dif, ndif +int side + +begin + dif = abs (x - screen_boundary[LEFT]) + side = LEFT + + ndif = abs (x - screen_boundary[RIGHT]) + if (ndif < dif) { + side = RIGHT + dif = ndif + } + + ndif = abs (y - screen_boundary[BOTTOM]) + if (ndif < dif) { + side = BOTTOM + dif = ndif + } + + ndif = abs (y - screen_boundary[TOP]) + if (ndif < dif) + side = TOP + + return (side) +end + + +# WL_RA_RANGE -- Determine the range in RA given a list of possible values. +# +# Description +# Determine the largest range in RA from the provided list of values. +# The problem here is that it is unknown which way the graph is oriented. +# To simplify the problem, it is assume that the graph range does not extend +# beyond a hemisphere and that all distances in RA is less than a hemisphere. +# This assumption is needed to decide when the 0 hour is on the graph. + +procedure wl_ra_range (ra, n_values, min, max, diff) + +double ra[ARB] # I: the possible RA values +int n_values # I: the number of possible RA values +double min # I/O: the minimum RA +double max # I/O: the maximum RA +double diff # I/O: the difference between minimum and maximum + +bool wrap +int i, j, n_diffs +pointer sp, max_array, min_array, ran_array +int wl_max_element_array() + +begin + call smark (sp) + call salloc (max_array, n_values * n_values, TY_DOUBLE) + call salloc (min_array, n_values * n_values, TY_DOUBLE) + call salloc (ran_array, n_values * n_values, TY_DOUBLE) + + # Check whether the RA is wrapped or not. + n_diffs = 0 + do i = 1, n_values { + if (ra[i] >= min && ra[i] <= max) + next + n_diffs = n_diffs + 1 + } + if (n_diffs > 0) + wrap = true + else + wrap = false + + n_diffs = 0 + for (i = 1; i <= n_values; i = i + 1) { + for (j = i + 1; j <= n_values; j = j + 1) { + n_diffs = n_diffs + 1 + call wl_getradif (ra[i], ra[j], Memd[min_array+n_diffs-1], + Memd[max_array+n_diffs-1], Memd[ran_array+n_diffs-1], + wrap) + } + } + + i = wl_max_element_array (Memd[ran_array], n_diffs) + min = Memd[min_array+i-1] + max = Memd[max_array+i-1] + diff = Memd[ran_array+i-1] + + call sfree (sp) +end + + +# WL_GETRADIFF -- Get differences in RA based on degrees. +# +# Description +# This procedure determines, given two values in degrees, the minimum, +# maximum, and difference of those values. The assumption is that no +# difference should be greater than half a circle. Based on this assumption, +# a difference is found and the minimum and maximum are determined. The +# maximum can be greater than 360 degrees. + +procedure wl_getradif (val1, val2, min, max, diff, wrap) + +double val1, val2 # I: the RA values +double min, max # O: the min RA and max RA (possibly > 360.0) +double diff # O: the min, max difference +bool wrap # I: is the ra wrapped ? + +begin + if (! wrap && (abs (val1 - val2) > HALF_CIRCLE)) + if (val1 < val2) { + min = val2 + max = val1 + FULL_CIRCLE + } else { + min = val1 + max = val2 + FULL_CIRCLE + } + else + if (val1 < val2) { + min = val1 + max = val2 + } else { + min = val2 + max = val1 + } + diff = max - min +end + + +define NRAGAP 26 + +# WL_ROUND_RA -- Modify the RA limits and calculate an interval to label. +# +# Description +# The RA limits determine by just the extremes of the window ususally do +# not fall on "reasonable" boundaries; i.e. essentially they are random +# numbers. However, for labelling purposes, it is nice to have grids and +# tick marks for "rounded" numbers- For RA, this means values close to +# whole hours, minutes, or seconds. For example, if the span across the +# plot is a few hours, the marks and labels should represent simply whole +# hours. This routine determines new RA limits based on this and some +# interval to produce marks between the newly revised limits. + +procedure wl_round_ra (longmin, longmax, longran, num_try, minimum, maximum, + major_interval) + +double longmin # I: longitude minimum +double longmax # I: longitude maximum +double longran # I: longitude range +int num_try # I: the number of intervals to try for +double minimum # O: the minimum RA value (in degrees) +double maximum # O: the maximum RA value (in degrees) +double major_interval # O: the appropriate interval (in degrees) for the + # major line marks. + +double ragap[NRAGAP] +double wl_check_arrayd(), wl_round_upd() +data ragap / 1.0D-4, 2.0D-4, 5.0D-4, 1.0D-3, 2.0D-3, 5.0D-3, + 0.01D0, 0.02D0, 0.05D0, 0.1D0, 0.2D0, 0.5D0, 1.0D0, + 2.0D0, 5.0D0, 10.0D0, 20.0D0, 30.0D0, 60.0D0, 120.0D0, + 300.0D0, 600.0D0, 1.2D3, 1.8D3, 3.6D3, 7.2D3 / + + +begin + major_interval = wl_check_arrayd (DEGTOST (longran) / num_try, + ragap, NRAGAP) + minimum = STTODEG (wl_round_upd (DEGTOST (longmin), major_interval) - + major_interval) + maximum = STTODEG (wl_round_upd (DEGTOST (longmax), major_interval)) + major_interval = STTODEG (major_interval) +end + + +define NDECGAP 28 + +# WL_ROUND_DEC -- Modify the DEC limits and calculate an interval to label. +# +# Description +# The DEC limits determine by just the extremes of the window ususally do +# not fall on "reasonable" boundaries; i.e. essentially they are random +# numbers. However, for labelling purposes, it is nice to have grids and +# tick marks for "rounded" numbers- For DEC, this means values close to +# whole degrees, minutes, or seconds. For example, if the span across the +# plot is a few degrees, the marks and labels should represent simply whole +# degrees. This routine determines new DEC limits based on this and some +# interval to produce marks between the newly revised limits. + +procedure wl_round_dec (latmin, latmax, latran, num_try, minimum, maximum, + major_interval) + +double latmin # I: the latitude minimum +double latmax # I: the latitude maximum +double latran # I: the latitude range +int num_try # I: number of intervals to try for +double minimum # O: the DEC minimum +double maximum # O: the DEC maximum +double major_interval # O: the labelling interval to use for major lines + +double decgap[NDECGAP] +double wl_check_arrayd(), wl_round_upd() +data decgap / 1.0D-4, 2.0D-4, 5.0D-4, 1.0D-3, 2.0D-3, 5.0D-3, + 0.01D0, 0.02D0, 0.05D0, 0.1D0, 0.2D0, 0.5D0, 1.0D0, + 2.0D0, 5.0D0, 10.0D0,20.0D0, 30.0D0, 60.0D0, 120.0d0, + 300.0D0, 600.0D0, 1.2D3, 1.8D3, 3.6D3, 7.2D3, 1.8D4, 3.6D4 / + +begin + major_interval = wl_check_arrayd (DEGTOSA (latran) / num_try, + decgap, NDECGAP) + minimum = SATODEG (wl_round_upd (DEGTOSA (latmin), major_interval) - + major_interval) + maximum = SATODEG (wl_round_upd (DEGTOSA (latmax), major_interval)) + major_interval = SATODEG (major_interval) + + # Make sure that the grid marking does not include the pole. + maximum = min (maximum, NORTH_POLE_LATITUDE - major_interval) + minimum = max (minimum, SOUTH_POLE_LATITUDE + major_interval) +end + + +# WL_GENERIC_ROUND -- Round the values (if possible). +# +# History +# 7Feb91 - Created by Jonathan D. Eisenhamer, STScI. + +procedure wl_generic_round (minimum, maximum, range, lbegin, lend, interval) + +double minimum, maximum, range # I: the raw input values +double lbegin, lend # O: the begin and end label points +double interval # O: the major label interval + +double amant, diff +int iexp, num +double wl_round_upd() + +begin + diff = log10 (abs (range) / 4.D0) + iexp = int (diff) + if (diff < 0) + iexp = iexp - 1 + + amant = diff - double (iexp) + if (amant < 0.15D0) + num = 1 + else if (amant < 0.50D0) + num = 2 + else if (amant < 0.85D0) + num = 5 + else + num = 10 + + interval = double (num) * 10.0D0 ** iexp + lbegin = wl_round_upd (minimum, interval) - interval + lend = wl_round_upd (maximum, interval) +end + + +# WL_ROUND_UPD -- Round X up to nearest whole multiple of Y. + +double procedure wl_round_upd (x, y) + +double x # I: value to be rounded +double y # I: multiple of X is to be rounded up in + +double z, r + +begin + if (x < 0.0D0) + z = 0.0D0 + else + z = y + r = y * double (int ((x + z) / y)) + + return (r) +end + + + +# WL_CHECK_ARRAYD -- Check proximity of array elements to each other. +# +# Description +# Returns the element of the array arr(n) which is closest to an exact +# value EX. + +double procedure wl_check_arrayd (ex, arr, n) + +double ex # I: the exact value +double arr[ARB] # I: the array of rounded values +int n # I: dimension of array of rounded values + +int j + +begin + for (j = 1; j < n && (ex - arr[j]) > 0.0D0; j = j + 1) + ; + if (j > 1 && j < n) + if (abs (ex - arr[j-1]) < abs (ex - arr[j])) + j = j - 1 + + return (arr[j]) +end diff --git a/pkg/utilities/nttools/stxtools/wcslab/wlutil.x b/pkg/utilities/nttools/stxtools/wcslab/wlutil.x new file mode 100644 index 00000000..c79b8f5e --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wlutil.x @@ -0,0 +1,390 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imio.h> +include <imhdr.h> +include <gset.h> +include <math.h> + +# WL_IMD_VIEWPORT -- Map the viewport and window of the image display. + +procedure wl_imd_viewport (frame, im, c1, c2, l1, l2, vl, vr, vb, vt) + +int frame # I: display frame to be overlayed +pointer im # I: pointer to the input image +real c1, c2, l1, l2 # I/O: input/output window +real vl, vr, vb, vt # I/O: input/output viewport + +int wcs_status, dim1, dim2, step1, step2 +pointer sp, frimage, frim, iw +real x1, x2, y1, y2, fx1, fx2, fy1, fy2, junkx, junky +real vx1, vx2, vy1, vy2, nx1, nx2, ny1, ny2 +pointer imd_mapframe(), iw_open() + + +begin + # If all of the viewport parameters were defined by the user + # use the default viewport and window. + if (! IS_INDEFR(vl) && ! IS_INDEFR(vr) && ! IS_INDEFR(vb) && + ! IS_INDEFR(vt)) + return + + # Allocate some memory. + call smark (sp) + call salloc (frimage, SZ_FNAME, TY_CHAR) + + # Open the requested display frame and get the loaded image name. + # If this name is blank, use the default viewport and window. + + frim = imd_mapframe (frame, READ_ONLY, YES) + iw = iw_open (frim, frame, Memc[frimage], SZ_FNAME, wcs_status) + if (Memc[frimage] == EOS || wcs_status == ERR) { + call iw_close (iw) + call imunmap (frim) + call sfree (sp) + return + } + + # Find the beginning and end points of the requested image section. + # We already know at this point that the input logical image is + # 2-dimensional. However this 2-dimensional section may be part of + # n-dimensional image. + + # X dimension. + dim1 = IM_VMAP(im,1) + step1 = IM_VSTEP(im,1) + if (step1 >= 0) { + x1 = IM_VOFF(im,dim1) + 1 + x2 = x1 + IM_LEN(im,1) - 1 + } else { + x1 = IM_VOFF(im,dim1) - 1 + x2 = x1 - IM_LEN(im,1) + 1 + } + + # Y dimension. + dim2 = IM_VMAP(im,2) + step2 = IM_VSTEP(im,2) + if (step2 >= 0) { + y1 = IM_VOFF(im,dim2) + 1 + y2 = y1 + IM_LEN(im,2) - 1 + } else { + y1 = IM_VOFF(im,dim2) - 1 + y2 = y1 - IM_LEN(im,2) + 1 + } + + # Get the frame buffer coordinates corresponding to the lower left + # and upper right corners of the image section. + + call iw_im2fb (iw, x1, y1, fx1, fy1) + call iw_im2fb (iw, x2, y2, fx2, fy2) + if (fx1 > fx2) { + junkx = fx1 + fx1 = fx2 + fx2 = junkx + } + if (fy1 > fy2) { + junky = fy1 + fy1 = fy2 + fy2 = junky + } + + # Check that some portion of the input image is in the display. + # If not select the default viewport and window coordinates. + if (fx1 > IM_LEN(frim,1) || fx2 < 1.0 || fy1 > IM_LEN(frim,2) || + fy2 < 1.0) { + call iw_close (iw) + call imunmap (frim) + call sfree (sp) + return + } + + # Compute a new viewport and window for X. + if (fx1 >= 1.0) { + vx1 = max (0.0, min (1.0, (fx1 - 0.5) / IM_LEN(frim,1))) + nx1 = 1.0 + } else { + vx1 = 0.0 + call iw_fb2im (iw, 1.0, 1.0, junkx, junky) + if (step1 >= 0) + nx1 = max (1.0, junkx - x1 + 1.0) + else + nx2 = max (1.0, junkx - x2 + 1.0) + } + if (fx2 <= IM_LEN(frim,1)) { + vx2 = max (0.0, min (1.0, (fx2 + 0.5) / IM_LEN(frim,1))) + nx2 = IM_LEN(im,1) + } else { + vx2 = 1.0 + call iw_fb2im (iw, real(IM_LEN(frim,1)), real (IM_LEN(frim,2)), + junkx, junky) + if (step1 >= 0) + nx2 = min (real (IM_LEN(im,1)), junkx - x1 + 1.0) + else + nx1 = min (real (IM_LEN(im,1)), junkx - x2 + 1.0) + } + + # Compute a new viewport and window for Y. + if (fy1 >= 1.0) { + vy1 = max (0.0, min (1.0, (fy1 - 0.5) / IM_LEN(frim,2))) + ny1 = 1.0 + } else { + vy1 = 0.0 + call iw_fb2im (iw, 1.0, 1.0, junkx, junky) + if (step2 >= 0) + ny1 = max (1.0, junky - y1 + 1) + else + ny2 = max (1.0, junky - y2 + 1) + } + if (fy2 <= IM_LEN(frim,2)) { + vy2 = max (0.0, min (1.0, (fy2 + 0.5) / IM_LEN(frim,2))) + ny2 = IM_LEN(im,2) + } else { + vy2 = 1.0 + call iw_fb2im (iw, real (IM_LEN(frim,1)), real (IM_LEN(frim,2)), + junkx, junky) + if (step2 >= 0) + ny2 = min (real (IM_LEN(im,2)), junky - y1 + 1.0) + else + ny1 = min (real (IM_LEN(im,2)), junky - y2 + 1.0) + } + + # Define a the new viewport and window. + if (IS_INDEFR(vl)) { + vl = vx1 + c1 = nx1 + } + if (IS_INDEFR(vr)) { + vr = vx2 + c2 = nx2 + } + if (IS_INDEFR(vb)) { + vb = vy1 + l1 = ny1 + } + if (IS_INDEFR(vt)) { + vt = vy2 + l2 = ny2 + } + + # Clean up. + call iw_close (iw) + call imunmap (frim) + call sfree (sp) +end + + +define EDGE1 0.1 +define EDGE2 0.9 +define EDGE3 0.12 +define EDGE4 0.92 + +# WL_MAP_VIEWPORT -- Set device viewport wcslab plots. If not specified by +# user, a default viewport centered on the device is used. + +procedure wl_map_viewport (gp, c1, c2, l1, l2, ux1, ux2, uy1, uy2, fill) + +pointer gp # I: pointer to graphics descriptor +real c1, c2, l1, l2 # I: the column and line limits +real ux1, ux2, uy1, uy2 # I/O: NDC coordinates of requested viewort +bool fill # I: fill viewport (vs preserve aspect ratio) + +int ncols, nlines +real xcen, ycen, ncolsr, nlinesr, ratio, aspect_ratio +real x1, x2, y1, y2, ext, xdis, ydis +bool fp_equalr() +real ggetr() +data ext /0.0625/ + +begin + ncols = nint (c2 - c1) + 1 + ncolsr = real (ncols) + nlines = nint (l2 - l1) + 1 + nlinesr = real (nlines) + + # Determine the standard window sizes. + if (fill) { + x1 = 0.0; x2 = 1.0 + y1 = 0.0; y2 = 1.0 + } else { + x1 = EDGE1; x2 = EDGE2 + y1 = EDGE3; y2 = EDGE4 + } + + # If any values were specified, then replace them here. + if (! IS_INDEFR(ux1)) + x1 = ux1 + if (! IS_INDEFR(ux2)) + x2 = ux2 + if (! IS_INDEFR(uy1)) + y1 = uy1 + if (! IS_INDEFR(uy2)) + y2 = uy2 + + # Calculate optimum viewport, as in NCAR's conrec, hafton. + if (! fill) { + ratio = min (ncolsr, nlinesr) / max (ncolsr, nlinesr) + if (ratio >= ext) { + if (ncols > nlines) + y2 = (y2 - y1) * nlinesr / ncolsr + y1 + else + x2 = (x2 - x1) * ncolsr / nlinesr + x1 + } + } + + xdis = x2 - x1 + ydis = y2 - y1 + xcen = (x2 + x1) / 2. + ycen = (y2 + y1) / 2. + + # So far, the viewport has been calculated so that equal numbers of + # image pixels map to equal distances in NDC space, regardless of + # the aspect ratio of the device. If the parameter "fill" has been + # set to no, the user wants to compensate for a non-unity aspect + # ratio and make equal numbers of image pixels map to into the same + # physical distance on the device, not the same NDC distance. + + if (! fill) { + aspect_ratio = ggetr (gp, "ar") + if (fp_equalr (aspect_ratio, 0.0)) + aspect_ratio = 1.0 + + if (aspect_ratio < 1.0) + # Landscape + xdis = xdis * aspect_ratio + else if (aspect_ratio > 1.0) + # Portrait + ydis = ydis / aspect_ratio + } + + ux1 = xcen - (xdis / 2.0) + ux2 = xcen + (xdis / 2.0) + uy1 = ycen - (ydis / 2.0) + uy2 = ycen + (ydis / 2.0) + + call gsview (gp, ux1, ux2, uy1, uy2) + call gswind (gp, c1, c2, l1, l2) +end + + +# WL_W2LD -- Transform world coordinates to logical coordinates. + +procedure wl_w2ld (wlct, flip, wx, wy, lx, ly, npts) + +pointer wlct # I: the MWCS coordinate transformation descriptor +int flip # I: true if the axes are transposed +double wx[npts], wy[npts] # I: the world coordinates +double lx[npts], ly[npts] # O: the logical coordinates +int npts # I: the number of points to translate + +begin + if (flip == YES) + call mw_v2trand (wlct, wx, wy, ly, lx, npts) + else + call mw_v2trand (wlct, wx, wy, lx, ly, npts) +end + + +# WL_L2WD -- Transform logical coordinates to world coordinates. + +procedure wl_l2wd (lwct, flip, lx, ly, wx, wy, npts) + +pointer lwct # I: the MWCS coordinate transformation descriptor +int flip # I: true if the axes are transposed +double lx[npts], ly[npts] # I: the logical coordinates +double wx[npts], wy[npts] # O: the world coordinates +int npts # I: the number of points to translate + +begin + if (flip == YES) + call mw_v2trand (lwct, ly, lx, wx, wy, npts) + else + call mw_v2trand (lwct, lx, ly, wx, wy, npts) +end + + +# WL_MAX_ELEMENT_ARRAY -- Return the index of the maximum array element. +# +# Description +# This function returns the index of the maximum value of the input array. + +int procedure wl_max_element_array (array, npts) + +double array[ARB] # I: the array to look through for the maximum +int npts # I: the number of points in the array + +int i, maximum + +begin + maximum = 1 + for (i = 2; i <= npts; i = i + 1) + if (array[i] > array[maximum]) + maximum = i + + return (maximum) +end + + +# WL_DISTANCED - Determine the distance between two points. + +double procedure wl_distanced (x1, y1, x2, y2) + +double x1, y1 # I: coordinates of point 1 +double x2, y2 # I: coordinates of point 2 + +double a, b + +begin + a = x1 - x2 + b = y1 - y2 + return (sqrt ((a * a) + (b * b))) +end + + +# WL_DISTANCER -- Determine the distance between two points. + +real procedure wl_distancer (x1, y1, x2, y2) + +real x1, y1 # I: coordinates of point 1 +real x2, y2 # I: coordinates of point 2 + +real a, b + +begin + a = x1 - x2 + b = y1 - y2 + return (sqrt ((a * a) + (b * b))) +end + + +# The dimensionality. +define N_DIM 2 + +# Define some memory management. +define ONER Memr[$1+$2-1] + +# WL_ROTATE -- Rotate a vector. + +procedure wl_rotate (x, y, npts, angle, nx, ny) + +real x[npts], y[npts] # I: the vectors to rotate +int npts # I: the number of points in the vectors +real angle # I: the angle to rotate (radians) +real nx[npts], ny[npts] # O: the transformed vectors + +pointer sp, center, mw +pointer mw_open(), mw_sctran() + +begin + # Get some memory. + call smark (sp) + call salloc (center, N_DIM, TY_REAL) + + mw = mw_open (NULL, N_DIM) + ONER(center,1) = 0. + ONER(center,2) = 0. + call mw_rotate (mw, -DEGTORAD( angle ), ONER(center,1), 3b) + call mw_v2tranr (mw_sctran (mw, "physical", "logical", 3b), + x, y, nx, ny, npts) + + call mw_close (mw) + call sfree (sp) +end diff --git a/pkg/utilities/nttools/stxtools/wcslab/wlwcslab.x b/pkg/utilities/nttools/stxtools/wcslab/wlwcslab.x new file mode 100644 index 00000000..156c9a8a --- /dev/null +++ b/pkg/utilities/nttools/stxtools/wcslab/wlwcslab.x @@ -0,0 +1,181 @@ +include <gio.h> +include <gset.h> +include "wcslab.h" +include "wcs_desc.h" + +# Define the memory structure for saving the graphics wcs. +define SAVE_BLOCK_SIZE 16 +define OLD_NDC_VIEW Memr[wcs_save_block-1+$1] +define OLD_NDC_WIND Memr[wcs_save_block+3+$1] +define OLD_PLT_VIEW Memr[wcs_save_block+7+$1] +define OLD_PLT_WIND Memr[wcs_save_block+11+$1] + +# WL_WCSLAB -- Label using a defined wcs. +# +# Description +# This routine uses the information in the WCSLAB descriptor to perform +# labelling. +# +# Before this routine can be called, several things must have already +# occured. They are as follows: +# 1 A call to wl_create must be made to create the WCSLAB descriptor. +# 2 The WCS_MW component must be set to the MWCS object of the +# desired transformations. +# 3 A call to wl_get_system_type must be made. +# 4 The graphics device must have been opened and the window defined. +# The WCS_GP component of the WCSLAB descriptor must be set to the +# graphics window descriptor. +# +# When done with this routine, the WL_GP and WL_MW components must be +# deallocated seperately. Then only wlab_destroy need be called to +# remove the WCSLAB descriptor. +# +#--------------------------------------------------------------------------- + +procedure wl_wcslab (wd) + +pointer wd # I: the WCSLAB descriptor + +int old_clip, old_pltype, old_txquality, old_wcs +pointer sp, wcs_save_block +real old_plwidth, old_txsize, old_txup +int gstati() +real gstatr() + +begin + # Allocate working space. + call smark(sp) + call salloc(wcs_save_block, SAVE_BLOCK_SIZE, TY_STRUCT) + + # Store certain graphics parameters. + old_plwidth = gstatr (WL_GP(wd), G_PLWIDTH) + old_txsize = gstatr (WL_GP(wd), G_TXSIZE) + old_txup = gstatr (WL_GP(wd), G_TXUP) + old_clip = gstati (WL_GP(wd), G_CLIP) + old_pltype = gstati (WL_GP(wd), G_PLTYPE) + old_txquality= gstati (WL_GP(wd), G_TXQUALITY) + old_wcs = gstati (WL_GP(wd), G_WCS) + + # Choose two other graphics wcs' for internal use. Save the wcs for + # later restoration. + if( old_wcs < MAX_WCS - 2 ) { + WL_NDC_WCS(wd) = old_wcs + 1 + WL_PLOT_WCS(wd) = WL_NDC_WCS(wd) + 1 + } else { + WL_NDC_WCS(wd) = old_wcs - 1 + WL_PLOT_WCS(wd) = WL_NDC_WCS(wd) - 1 + } + call gseti(WL_GP(wd), G_WCS, WL_NDC_WCS(wd)) + call ggview(WL_GP(wd), OLD_NDC_VIEW(LEFT), OLD_NDC_VIEW(RIGHT), + OLD_NDC_VIEW(BOTTOM), OLD_NDC_VIEW(TOP)) + call ggwind(WL_GP(wd), OLD_NDC_WIND(LEFT), OLD_NDC_WIND(RIGHT), + OLD_NDC_WIND(BOTTOM), OLD_NDC_WIND(TOP)) + call gseti(WL_GP(wd), G_WCS, WL_PLOT_WCS(wd)) + call ggview(WL_GP(wd), OLD_PLT_VIEW(LEFT), OLD_PLT_VIEW(RIGHT), + OLD_PLT_VIEW(BOTTOM), OLD_PLT_VIEW(TOP)) + call ggwind(WL_GP(wd), OLD_PLT_WIND(LEFT), OLD_PLT_WIND(RIGHT), + OLD_PLT_WIND(BOTTOM), OLD_PLT_WIND(TOP)) + + # Set the graphics device the way wcslab requires it. + call gseti (WL_GP(wd), G_WCS, old_wcs) + call wl_graphics (wd) + + # Determine basic characteristics of the plot. + call wl_setup (wd) + + # Plot the grid lines. + call wl_grid (wd) + + # Put the grid labels on the lines. + if (WL_LABON(wd) == YES) + call wl_label (wd) + + # Restore the original graphics wcs. + call gseti(WL_GP(wd), G_WCS, WL_NDC_WCS(wd)) + call gsview(WL_GP(wd), OLD_NDC_VIEW(LEFT), OLD_NDC_VIEW(RIGHT), + OLD_NDC_VIEW(BOTTOM), OLD_NDC_VIEW(TOP)) + call gswind(WL_GP(wd), OLD_NDC_WIND(LEFT), OLD_NDC_WIND(RIGHT), + OLD_NDC_WIND(BOTTOM), OLD_NDC_WIND(TOP)) + call gseti(WL_GP(wd), G_WCS, WL_PLOT_WCS(wd)) + call gsview(WL_GP(wd), OLD_PLT_VIEW(LEFT), OLD_PLT_VIEW(RIGHT), + OLD_PLT_VIEW(BOTTOM), OLD_PLT_VIEW(TOP)) + call gswind(WL_GP(wd), OLD_PLT_WIND(LEFT), OLD_PLT_WIND(RIGHT), + OLD_PLT_WIND(BOTTOM), OLD_PLT_WIND(TOP)) + + # Restore original graphics state. + call gsetr (WL_GP(wd), G_PLWIDTH, old_plwidth) + call gsetr (WL_GP(wd), G_TXSIZE, old_txsize) + call gsetr (WL_GP(wd), G_TXUP, old_txup) + call gseti (WL_GP(wd), G_CLIP, old_clip) + call gseti (WL_GP(wd), G_PLTYPE, old_pltype) + call gseti (WL_GP(wd), G_TXQUALITY, old_txquality) + call gseti (WL_GP(wd), G_WCS, old_wcs) + + call sfree (sp) +end + + +# WL_GRAPHICS -- Setup the graphics device appropriate for the occasion. + +procedure wl_graphics (wd) + +pointer wd # I: the WCSLAB descriptor + +real relative_size, vl, vr, vb, vt +real ggetr() + +begin + # Setup a graphics WCS that mimics the NDC coordinate WCS, + # but with clipping. + call ggview (WL_GP(wd), vl, vr, vb, vt) + call gseti (WL_GP(wd), G_WCS, WL_NDC_WCS(wd)) + call gsview (WL_GP(wd), vl, vr, vb, vt) + call gswind (WL_GP(wd), vl, vr, vb, vt) + call gseti (WL_GP(wd), G_CLIP, YES) + + # Setup the initial viewport. + WL_NEW_VIEW(wd,LEFT) = vl + WL_NEW_VIEW(wd,RIGHT) = vr + WL_NEW_VIEW(wd,BOTTOM) = vb + WL_NEW_VIEW(wd,TOP) = vt + + # Setup some parameters. + call gseti (WL_GP(wd), G_PLTYPE, GL_SOLID) + call gsetr (WL_GP(wd), G_PLWIDTH, LINE_SIZE) + + # Draw the edges of the viewport. + call gamove (WL_GP(wd), vl, vb) + call gadraw (WL_GP(wd), vr, vb) + call gadraw (WL_GP(wd), vr, vt) + call gadraw (WL_GP(wd), vl, vt) + call gadraw (WL_GP(wd), vl, vb) + + # Determine the tick mark size. + relative_size = max (abs (vr - vl), abs (vt - vb )) + WL_MAJ_TICK_SIZE(wd) = relative_size * WL_MAJ_TICK_SIZE(wd) + WL_MIN_TICK_SIZE(wd) = relative_size * WL_MIN_TICK_SIZE(wd) + + # Determine various character sizes. + WL_TITLE_SIZE(wd) = WL_TITLE_SIZE(wd) * relative_size + WL_AXIS_TITLE_SIZE(wd) = WL_AXIS_TITLE_SIZE(wd) * relative_size + WL_LABEL_SIZE(wd) = WL_LABEL_SIZE(wd) * relative_size + + # Now setup the general plotting WCS. + call gseti (WL_GP(wd), G_WCS, WL_PLOT_WCS(WD)) + call gsview (WL_GP(wd), vl, vr, vb, vt) + vl = real (WL_SCREEN_BOUNDARY(wd,LEFT)) + vr = real (WL_SCREEN_BOUNDARY(wd,RIGHT)) + vb = real (WL_SCREEN_BOUNDARY(wd,BOTTOM)) + vt = real (WL_SCREEN_BOUNDARY(wd,TOP)) + call gswind (WL_GP(wd), vl, vr, vb, vt) + call gseti (WL_GP(wd), G_CLIP, YES) + + # Set some characteristics of the graphics device. + call gseti (WL_GP(wd), G_TXQUALITY, GT_HIGH) + call gseti (WL_GP(wd), G_CLIP, YES) + call gsetr (WL_GP(wd), G_PLWIDTH, LINE_SIZE) + + # Determine the number of segments a "line" should consist of. + WL_LINE_SEGMENTS(wd) = int (min (ggetr (WL_GP(wd), "xr"), + ggetr (WL_GP(wd), "yr")) / 5) +end |