Initial commit

18 files changed, 5108 insertions, 0 deletions

diff --git a/pkg/dataio/fits/fits_cards.x b/pkg/dataio/fits/fits_cards.x
new file mode 100644
index 00000000..0ddfa230
--- /dev/null
+++ b/pkg/dataio/fits/fits_cards.x
@@ -0,0 +1,292 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include "wfits.h"
+
+# WFT_STANDARD_CARD -- Procedure for fetching the minimum header
+# parameters required by fits. The end card is encoded separately.
+
+int procedure wft_standard_card (cardno, im, fits, axisno, card)
+
+int	cardno		# number of FITS standard card
+pointer	im		# pointer to the IRAF image
+pointer	fits		# pointer to the FITS structure
+int	axisno		# axis number
+char	card[ARB]	# FITS card image
+
+char	keyword[LEN_KEYWORD]
+int	len_object
+int	strlen()
+errchk	wft_encodeb, wft_encodei, wft_encodel, wft_encode_axis
+
+begin
+	# Get mandatory keywords.
+	switch (cardno) {
+	case FIRST_CARD:
+	    if (XTENSION(fits) == EXT_PRIMARY) {
+	        call wft_encodeb ("SIMPLE", YES, card, "FITS STANDARD")
+	    } else {
+		len_object = max (min (LEN_OBJECT, strlen ("IMAGE")),
+		    LEN_STRING)
+		call wft_encodec ("XTENSION", "IMAGE", len_object, card,
+		    "IMAGE EXTENSION")
+	    }
+	case SECOND_CARD:
+	    call wft_encodei ("BITPIX", FITS_BITPIX(fits), card,
+	    "FITS BITS/PIXEL")
+	case THIRD_CARD:
+	    call wft_encodei ("NAXIS", NAXIS(im), card, "NUMBER OF AXES")
+	default:
+	    call wft_encode_axis ("NAXIS", keyword, axisno)
+	    call wft_encodel (keyword, NAXISN(im, axisno), card, "")
+	    axisno = axisno + 1
+	}
+
+	return (YES)
+end
+
+
+# WFT_OPTION_CARD -- Procedure for fetching  optional FITS header parameters.
+# At present these are bscale, bzero, bunit, blank, object, origin, date,
+# irafmax, irafmin, iraf type and iraf bits per pixel. Blank is only encoded
+# if there are a nonzero number of blanks in the IRAF image. Bunit and object
+# are only encoded if the appropriate IRAF strings are defined. Bzero, bscale,
+# irafmax, irafmin, iraf type and iraf bits per pixel are only encoded if
+# there is a pixel file.
+
+int procedure wft_option_card (im, fits, optiono, card)
+
+pointer	im		# pointer to the IRAF image
+pointer	fits		# pointer to FITS structure
+int	optiono		# number of the option card
+char	card[ARB]	# FITS card image
+
+char	datestr[LEN_DATE]
+int	len_object, stat
+int	strlen()
+errchk	wft_encoded, wft_encodec, wft_encode_blank, wft_encoder, wft_encodei
+errchk	wft_encode_date 
+include "wfits.com"
+
+begin
+	 stat = YES
+
+	 # get optional keywords
+	 switch (optiono) {
+	 case KEY_EXTEND:
+	    if (XTENSION(fits) == EXT_IMAGE || wextensions == NO)
+		stat = NO
+	    else
+		call wft_encodeb ("EXTEND", YES, card,
+		    "STANDARD EXTENSIONS MAY BE PRESENT")
+	 case KEY_PCOUNT:
+	    if (XTENSION(fits) == EXT_PRIMARY)
+		stat = NO
+	    else
+		call wft_encodei ("PCOUNT", 0, card, "NO RANDOM PARAMETERS")
+	 case KEY_GCOUNT:
+	    if (XTENSION(fits) == EXT_PRIMARY)
+		stat = NO
+	    else
+		call wft_encodei ("GCOUNT", 1, card, "ONLY ONE GROUP")
+	 case KEY_BSCALE:
+	    if ((NAXIS(im) <= 0) || (FITS_BITPIX(fits) < 0))
+		stat = NO
+	    else {
+	        call wft_encoded ("BSCALE", BSCALE(fits), card,
+	        "REAL = TAPE*BSCALE + BZERO", NDEC_DOUBLE)
+	    }
+	 case KEY_BZERO:
+	    if ((NAXIS(im) <= 0) || (FITS_BITPIX(fits) < 0))
+		stat = NO
+	    else
+	        call wft_encoded ("BZERO", BZERO(fits), card, "", NDEC_DOUBLE)
+	 case KEY_BUNIT:
+	    stat = NO
+	 case KEY_BLANK:
+	    stat = NO
+	    #if (NBPIX(im) == 0)
+		#stat = NO
+	    #else
+		#call wft_encode_blank ("BLANK", BLANK_STRING(fits), card,
+		#"TAPE VALUE OF BLANK PIXEL")
+	case KEY_OBJECT:
+	    if (OBJECT(im) == EOS)
+		stat = NO
+	    else {
+		len_object = max (min (LEN_OBJECT, strlen (OBJECT(im))),
+		    LEN_STRING)
+		call wft_encodec ("OBJECT", OBJECT(im), len_object, card, "")
+	    }
+	case KEY_ORIGIN:
+	    call wft_encodec ("ORIGIN", "KPNO-IRAF", LEN_ORIGIN, card, "")
+	case KEY_DATE:
+	    call wft_encode_date (datestr, LEN_DATE)
+	    len_object = max (min (LEN_OBJECT, strlen (datestr)), LEN_STRING)
+	    call wft_encodec ("DATE", datestr, len_object, card, "")
+	case KEY_IRAFNAME:
+	    len_object = max (min (LEN_OBJECT, strlen (IRAFNAME(fits))),
+		LEN_STRING)
+	    call wft_encodec ("IRAFNAME", IRAFNAME(fits), len_object, card, 
+		"NAME OF IRAF IMAGE FILE")
+	case KEY_IRAFMAX:
+	    if (NAXIS(im) <= 0)
+		stat = NO
+	    else
+	        call wft_encoder ("IRAF-MAX", IRAFMAX(fits), card, "DATA MAX",
+				NDEC_REAL)
+	case KEY_IRAFMIN:
+	    if (NAXIS(im) <= 0)
+		stat = NO
+	    else
+	        call wft_encoder ("IRAF-MIN", IRAFMIN(fits), card, "DATA MIN",
+		NDEC_REAL)
+	case KEY_IRAFBP:
+	    if (NAXIS(im) <= 0)
+		stat = NO
+	    else
+	        call wft_encodei ("IRAF-BPX", DATA_BITPIX(fits), card,
+	    	    "DATA BITS/PIXEL")
+	case KEY_IRAFTYPE:
+	    if (NAXIS(im) <= 0)
+		stat = NO
+	    else
+	        call wft_encodec ("IRAFTYPE", TYPE_STRING(fits), LEN_STRING,
+	    		    card, "PIXEL TYPE")
+	default:
+	    stat = NO
+	}
+
+	optiono = optiono + 1
+
+	return (stat)
+end
+
+
+# WFT_HISTORY_CARD -- Procedure to fetch a single history line, trim newlines
+# and pad with blanks to size LEN_CARD in order to create a FITS HISTORY card.
+
+int procedure wft_history_card (im, hp, card)
+
+pointer	im		# pointer to the IRAF image
+int	hp		# pointer to first character to extract from string
+char	card[ARB]	# FITS card image
+
+char	cval
+char	chfetch()
+
+begin
+	if (chfetch (HISTORY(im), hp, cval) == EOS)
+	    return (NO)
+	else {
+	    hp = hp - 1
+	    call strcpy ("HISTORY ", card, LEN_KEYWORD)
+	    call wft_fits_card (HISTORY(im), hp, card, COL_VALUE - 2, LEN_CARD,
+		'\n')
+	    return (YES)
+	}
+end
+
+
+# WFT_UNKNOWN_CARD  -- Procedure to fetch a single unknown
+# "line", trim newlines and pad blanks to size LEN_CARD in order to
+# create an unknown keyword card. At present user area information is
+# assumed to be in the form of FITS card images, less then or equal to
+# 80 characters and delimited by a newline.
+
+int procedure wft_unknown_card (im, up, card)
+
+pointer	im		# pointer to the IRAF image
+int	up		# pointer to next character in the unknown string
+char	card[ARB]	# FITS card image
+
+char	cval
+int	stat, axis, index
+char	chfetch()
+int	strmatch(), ctoi()
+
+begin
+	if (chfetch (UNKNOWN(im), up, cval) == EOS)
+	    return (NO)
+	else {
+	    up = up - 1
+	    stat = NO
+	    while (stat == NO)  {
+	        call wft_fits_card (UNKNOWN(im), up, card, 1, LEN_CARD, '\n')
+		if (card[1] == EOS)
+		    break
+	        if (strmatch (card, "^GROUPS  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^SIMPLE  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^BITPIX  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^NAXIS   ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^NAXIS") != 0) {
+		    index = LEN_NAXIS_KYWRD + 1
+		    if (ctoi (card, index, axis) > 0)
+		        stat = NO
+		    else
+			stat = YES
+	        } else if (strmatch (card, "^GCOUNT  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^PCOUNT  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^PSIZE   ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^BSCALE  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^BZERO   ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^BLANK   ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^IRAF-MAX") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^IRAF-MIN") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^IRAFTYPE") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^IRAF-B/P") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^IRAF-BPX") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^FILENAME") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^IRAFNAME") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^EXTEND  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^EXTNAME ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^EXTVER  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^INHERIT ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^IRAF-TLM") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^OBJECT  ") != 0) {
+		    stat = NO
+	        } else if (strmatch (card, "^END     ") != 0) {
+		    stat = NO
+	        } else
+                    stat = YES
+	    }
+
+	    return (stat)
+	}
+end
+
+
+# WFT_LAST_CARD -- Procedure to encode the FITS end card.
+
+int procedure wft_last_card (card)
+
+char	card[ARB]	# FITS card image
+
+begin
+	call sprintf (card, LEN_CARD, "%-8.8s  %70w")
+	    call pargstr ("END")
+
+	return (YES)
+end
diff --git a/pkg/dataio/fits/fits_files.x b/pkg/dataio/fits/fits_files.x
new file mode 100644
index 00000000..ce2c553c
--- /dev/null
+++ b/pkg/dataio/fits/fits_files.x
@@ -0,0 +1,374 @@
+include <ctype.h>
+include <plset.h>
+
+define	DEF_MAXNCOLS	1000
+define	DEF_MAXNLINES	4096
+define	MAX_NRANGES	100
+
+
+# RFT_FLIST -- Decode a list of files and associated extensions into a pixel
+# list.
+
+pointer procedure rft_flist (file_list, first_file, last_file, nfiles) 
+
+char	file_list[ARB]		# the input file list string
+int	first_file		# the first file in the list
+int	last_file		# the last file in the list
+int	nfiles			# the number of files in the list
+
+int	i, j, maxncols, maxnlines, nrfiles, rp, rbegin, rend, rstep
+int	last_ext, ebegin, eend, estep, ep, nefiles
+pointer	sp, extensions, str, axes, pl
+
+bool	pl_linenotempty()
+int	rft_gfranges()
+pointer	pl_create()
+
+begin
+	# Allocate some working space.
+	call smark (sp)
+	call salloc (extensions, SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (axes, 2, TY_INT)
+
+	# Initialize the file list.
+	pl = NULL
+	maxncols = DEF_MAXNCOLS
+	maxnlines = DEF_MAXNLINES
+
+	repeat {
+
+	    # Initialize the file counter parameters.
+	    first_file = INDEFI
+	    last_file = INDEFI
+	    nfiles = 0
+	    rp = 1
+
+	    # Open the file list.
+	    Memi[axes] = maxncols
+	    Memi[axes+1] = maxnlines
+	    pl = pl_create (2, Memi[axes], 1)
+
+	    # Decode the file list.
+	    nrfiles = rft_gfranges (file_list, rp, YES, maxnlines, rbegin,
+	        rend, rstep, NO, Memc[extensions])
+	    while (nrfiles > 0) {
+
+		# Check the file number limits and terminate the loop if
+		# the current list size is exceeded.
+		if (IS_INDEFI(first_file))
+		    first_file = rbegin
+		else
+		    first_file = min (rbegin, first_file)
+		if (IS_INDEFI(last_file))
+		    last_file = rend
+		else
+		    last_file = max (last_file, rend)
+		if (last_file > maxnlines)
+		    break
+
+		# Initialize the extensions list decoding.
+		ep = 1
+		last_ext = INDEFI
+
+		# Decode the associated extensions files. If the extensions
+		# list is empty
+		nefiles = rft_gfranges (Memc[extensions], ep, YES, maxncols,
+		    ebegin, eend, estep, YES, Memc[str])
+		while (nefiles > 0) {
+
+		    # Check the  extensions number limits and quit if they
+		    # are exceeded.
+		    if (IS_INDEFI(last_ext))
+		        last_ext = eend
+		    else
+		        last_ext = max (last_ext, eend)
+		    if (last_ext > maxncols)
+		        break
+
+		    # Set the appropriate elements in the list.
+		    if (rstep == 1) {
+			if (estep == 1)
+		    	    call pl_box (pl, ebegin, rbegin, eend,
+			        rend, PIX_SET + PIX_VALUE(1))
+			else {
+			    do i = ebegin, eend, estep
+			        call pl_box (pl, i, rbegin, i, rend, 
+			            PIX_SET + PIX_VALUE(1))
+			}
+		    } else {
+		        do i = rbegin, rend, rstep {
+			    do j = ebegin, eend, estep
+		                call pl_point (pl, j, i, PIX_SET +
+			            PIX_VALUE(1))
+			}
+		    }
+
+		    nefiles = rft_gfranges (Memc[extensions], ep, NO, maxncols,
+		        ebegin, eend, estep, YES, Memc[str])
+		}
+
+		# Break if an extensions list decode error occurs.
+		if (nefiles == 0)
+		    break
+
+	        nrfiles = rft_gfranges (file_list, rp, NO, maxnlines, rbegin,
+		    rend, rstep, NO, Memc[extensions])
+	    }
+
+	    # Break if a file or extensions list decode error ocurred.
+	    if (nrfiles == 0 || nefiles == 0)
+		break
+
+	    # If the file or extensions list is larger than the current maximum
+	    # then free the list increase the default space and repeat the
+	    # procedure.
+
+	    if (!IS_INDEFI(last_file)) {
+	        if (last_file > maxnlines) {
+		    if (pl != NULL)
+		        call pl_close (pl)
+		    pl = NULL
+		    maxnlines = maxnlines + DEF_MAXNLINES
+	        } else {
+		    do i = first_file, last_file {
+		        Memi[axes] = 1
+		        Memi[axes+1] = i
+		        if (pl_linenotempty (pl, Memi[axes]))
+			    nfiles = nfiles + 1
+		    }
+		}
+	    }
+
+	    if (!IS_INDEFI(last_ext)) {
+		if (last_ext > maxncols) {
+		    if (pl != NULL)
+		        call pl_close (pl)
+		    pl = NULL
+		    maxncols = maxncols + DEF_MAXNCOLS
+		}
+	    }
+
+	} until (pl != NULL)
+
+	# Free space.
+	call sfree (sp)
+
+	# If the file list is empty or a decode error occurred return
+	# a NULL list, otherwise return the pointer to the list.
+
+	if (nfiles <= 0 || nrfiles == 0 || nefiles == 0) {
+	    nfiles = 0
+	    first_file = INDEFI
+	    last_file = INDEFI
+	    call pl_close (pl)
+	    return (NULL)
+	} else
+	    return (pl)
+end
+
+
+# RFT_GFRANGES -- Parse a string containing a list of integer numbers or
+# ranges, delimited by either spaces or commas.  Each range may have an
+# associated extensions lists delimited by square brackets. Return as output
+# each range in sequence. Range limits must be positive nonnegative integers.
+# EOF is returned if the  end of the file_list is reached.  0 is returned if a
+# conversion error takes place. Otherwise the number of elements in the
+# range is returned.
+
+int procedure rft_gfranges (range_string, ip, firstr, rmax, rbegin, rend,
+	rstep, zeroindex, extensions)
+
+char    range_string[ARB]       # range string to be decoded
+int	ip			# the range string pointer
+int	firstr			# first range to be returned
+int	rmax			# the maximum file number
+int	rbegin			# the begining of the range			
+int	rend			# the end of the range
+int	rstep			# the range step size
+int	zeroindex		# allow zero indexing ?
+char	extensions[ARB]		# the output extensions string
+
+int	ep, itemp
+int	ctoi()
+
+begin
+	# Initialize.
+	if (zeroindex == YES) {
+	    rbegin = 0
+            rend = rmax - 1
+	} else {
+	    rbegin = 1
+            rend = rmax
+	}
+        rstep = 1
+	extensions[1] = EOS
+
+	# Return default range if the range string is NULL.
+	if (range_string[ip] == EOS) {
+	    rbegin = 1
+	    rend = rmax
+	    if (firstr == YES)
+	        return (rend)
+	    else
+	        return (EOF)
+	}
+
+        # Skip delimiters
+        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+            ip = ip + 1
+
+        # Get first limit.
+        # Must be a number, '-', 'x', *, or EOS.
+        if (range_string[ip] == EOS) {                      # end of list
+	    rbegin = 1
+	    rend = rmax
+	    if (firstr == YES)
+	        return (rend)
+	    else
+		return (EOF)
+	} else if (range_string[ip] == '*') {
+	    ;
+        } else if (range_string[ip] == '-') {
+            ;
+        } else if (range_string[ip] == 'x') {
+            ;
+        } else if (IS_DIGIT(range_string[ip])) {              # ,n..
+            if (ctoi (range_string, ip, rbegin) == 0)
+                return (0)
+	    else if (zeroindex == NO) {
+	        if (rbegin <= 0)
+		    return (0)
+	    } else {
+	        if (rbegin < 0)
+		    return (0)
+	    }
+        } else
+            return (0)
+
+	# Extract extensions file list.
+	if (range_string[ip] == '[') {
+	    ip = ip + 1
+	    ep = 1
+	    while (range_string[ip] != EOS) {
+		if (range_string[ip] == ']') {
+		    ip = ip + 1
+		    break
+		}
+		extensions[ep] = range_string[ip]
+		ip = ip + 1
+		ep = ep + 1
+	    }
+	    extensions[ep] = EOS
+	    if (range_string[ip-1] != ']')
+		return (0)
+	}
+
+        # Skip delimiters
+        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+            ip = ip + 1
+
+        # Get last limit
+        # Must be '-', 'x', or '*' otherwise last = first.
+	if (range_string[ip] == '*')
+	    ;
+        else if (range_string[ip] == 'x')
+            ;
+        else if (range_string[ip] == '-') {
+            ip = ip + 1
+            while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+                ip = ip + 1
+            if (range_string[ip] == EOS)
+                ;
+            else if (IS_DIGIT(range_string[ip])) {
+                if (ctoi (range_string, ip, rend) == 0)
+                    return (0)
+		else if (zeroindex == NO) {
+		    if (rend <= 0)
+		        return (0)
+		} else {
+		    if (rend < 0)
+		        return (0)
+		}
+            } else if (range_string[ip] == 'x')
+                    ;
+            else
+                return (0)
+        } else
+                rend = rbegin
+
+	# Skip extensions files for now.
+	if (range_string[ip] == '[') {
+	    ip = ip + 1
+	    ep = 1
+	    while (range_string[ip] != EOS) {
+		if (range_string[ip] == ']') {
+		    ip = ip + 1
+		    break
+		}
+		extensions[ep] = range_string[ip]
+		ip = ip + 1
+		ep = ep + 1
+	    }
+	    extensions[ep] = EOS
+	    if (range_string[ip-1] != ']')
+	        return (0)
+	}
+
+        # Skip delimiters
+        while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+            ip = ip + 1
+
+        # Get step.
+        # Must be 'x' or assume default step.
+        if (range_string[ip] == '*')
+	    ip = ip + 1
+        else if (range_string[ip] == 'x') {
+            ip = ip + 1
+            while (IS_WHITE(range_string[ip]) || range_string[ip] == ',')
+                ip = ip + 1
+            if (range_string[ip] == EOS)
+                ;
+            else if (IS_DIGIT(range_string[ip])) {
+                if (ctoi (range_string, ip, rstep) == 0)
+                    ;
+	        else if (rstep <= 0)
+		    return (0)
+            } else if (range_string[ip] == '-')
+                ;
+            else
+                return (0)
+        }
+
+	# Skip extensions files for now.
+	if (range_string[ip] == '[') {
+	    ip = ip + 1
+	    ep = 1
+	    while (range_string[ip] != EOS) {
+		if (range_string[ip] == ']') {
+		    ip = ip + 1
+		    break
+		}
+		extensions[ep] = range_string[ip]
+		ip = ip + 1
+		ep = ep + 1
+	    }
+	    extensions[ep] = EOS
+	    if (range_string[ip-1] != ']')
+	        return (0)
+	}
+
+        # Output the range triple.
+	if (rend < rbegin) {
+	    itemp = rbegin
+	    rbegin = rend
+	    rend = itemp
+	}
+	if (zeroindex == YES) {
+	    rbegin = rbegin + 1
+	    rend = rend + 1
+	}
+        return (abs (rend - rbegin) / rstep + 1 )
+end
+
+
diff --git a/pkg/dataio/fits/fits_params.x b/pkg/dataio/fits/fits_params.x
new file mode 100644
index 00000000..6911a925
--- /dev/null
+++ b/pkg/dataio/fits/fits_params.x
@@ -0,0 +1,248 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <time.h>
+include "wfits.h"
+
+# WFT_ENCODEB -- Procedure to encode a boolean parameter into a FITS card.
+
+procedure wft_encodeb (keyword, param, card, comment)
+
+char	keyword[ARB]	# FITS keyword
+int	param		# integer parameter equal to YES/NO
+char	card[ARB]	# FITS card image
+char	comment[ARB]	# FITS comment string
+
+char	truth
+
+begin
+	if (param == YES)
+	    truth = 'T'
+	else
+	    truth = 'F'
+
+	call sprintf (card, LEN_CARD, "%-8.8s= %20c  /  %-45.45s")
+	    call pargstr (keyword)
+	    call pargc (truth)
+	    call pargstr (comment)
+end
+
+
+# WFT_ENCODEI -- Procedure to encode an integer parameter into a FITS card.
+
+procedure wft_encodei (keyword, param, card, comment)
+
+char	keyword[ARB]	# FITS keyword
+int	param		# integer parameter
+char	card[ARB]	# FITS card image
+char	comment[ARB]	# FITS comment string
+
+begin
+	call sprintf (card, LEN_CARD, "%-8.8s= %20d  /  %-45.45s")
+	    call pargstr (keyword)
+	    call pargi (param)
+	    call pargstr (comment)
+end
+
+
+# WFT_ENCODEL -- Procedure to encode a long parameter into a FITS card.
+
+procedure wft_encodel (keyword, param, card, comment)
+
+char	keyword[ARB]		# FITS keyword
+long	param			# long integer parameter
+char	card[ARB]		# FITS card image
+char	comment[ARB]		# FITS comment string
+
+begin
+	call sprintf (card, LEN_CARD, "%-8.8s= %20d  /  %-45.45s")
+	    call pargstr (keyword)
+	    call pargl (param)
+	    call pargstr (comment)
+end
+
+
+# WFT_ENCODER -- Procedure to encode a real parameter into a FITS card.
+
+procedure wft_encoder (keyword, param, card, comment, precision)
+
+char	keyword[ARB]		# FITS keyword
+real	param			# real parameter
+char	card[ARB]		# FITS card image
+char	comment[ARB]		# FITS comment card
+int	precision		# precision of real
+
+begin
+	call sprintf (card, LEN_CARD, "%-8.8s= %20.*e  /  %-45.45s")
+	    call pargstr (keyword)
+	    call pargi (precision)
+	    call pargr (param)
+	    call pargstr (comment)
+end
+
+
+# WFT_ENCODED -- Procedure to encode a double parameter into a FITS card.
+
+procedure wft_encoded (keyword, param, card, comment, precision)
+
+char	keyword[ARB]		# FITS keyword
+double	param			# double parameter
+char	card[ARB]		# FITS card image
+char	comment[ARB]		# FITS comment string
+int	precision		# FITS precision
+
+begin
+	call sprintf (card, LEN_CARD, "%-8.8s= %20.*e  /  %-45.45s")
+	    call pargstr (keyword)
+	    call pargi (precision)
+	    call pargd (param)
+	    call pargstr (comment)
+end
+
+
+# WFT_ENCODE_AXIS -- Procedure to add the axis number to axis dependent
+# keywords.
+
+procedure wft_encode_axis (root, keyword, axisno)
+
+char	root[ARB]		# FITS root keyword
+char	keyword[ARB]		# FITS keyword
+int	axisno			# FITS axis number
+
+begin
+	call strcpy (root, keyword, LEN_KEYWORD)
+	call sprintf (keyword, LEN_KEYWORD, "%-5.5s%-3.3s")
+	    call pargstr (root)
+	    call pargi (axisno)
+end
+
+
+# WFT_ENCODEC -- Procedure to encode an IRAF string parameter into a FITS card.
+
+procedure wft_encodec (keyword, param, maxch, card, comment)
+
+char	keyword[ARB]	# FITS keyword
+char	param[ARB]	# FITS string parameter
+int	maxch		# maximum number of characters in string parameter
+char	card[ARB]	# FITS card image
+char	comment[ARB]	# comment string
+
+char	strparam[LEN_ALIGN+2]
+int	maxchar, nblanks
+
+begin
+	maxchar = min (maxch, LEN_OBJECT)
+	if (maxchar <= LEN_ALIGN - 1) {
+	    strparam[1] = '\''
+	    call sprintf (strparam[2], maxchar, "%*.*s")
+		call pargi (-maxchar)
+		call pargi (maxchar)
+		call pargstr (param)
+	    strparam[maxchar+2] = '\''
+	    strparam[maxchar+3] = EOS
+            call sprintf (card, LEN_CARD, "%-8.8s= %-20.20s  /  %-45.45s")
+	        call pargstr (keyword)
+	        call pargstr (strparam)
+	        call pargstr (comment)
+	} else {
+	    nblanks = LEN_OBJECT - maxchar
+	    if (comment[1] == EOS)
+                call sprintf (card, LEN_CARD, "%-8.8s= '%*.*s'     %*.*s")
+	    else
+                call sprintf (card, LEN_CARD, "%-8.8s= '%*.*s'  /  %*.*s")
+	        call pargstr (keyword)
+	        call pargi (-maxchar)
+	        call pargi (maxchar)
+	        call pargstr (param)
+	        call pargi (-nblanks)
+	        call pargi (nblanks)
+	        call pargstr (comment)
+	}
+end
+
+
+# WFT_ENCODE_BLANK -- Procedure to encode the FITS blank parameter. Necessary
+# because the 32 bit blank value equals INDEFL.
+
+procedure wft_encode_blank (keyword, blank_str, card, comment)
+
+char	keyword[ARB]		# FITS keyword
+char	blank_str[ARB]		# string containing values of FITS blank integer
+char	card[ARB]		# FITS card image
+char	comment[ARB]		# FITS comment string
+
+begin
+    call sprintf (card, LEN_CARD, "%-8.8s= %20.20s  /  %-45.45s")
+	call pargstr (keyword)
+	call pargstr (blank_str)
+	call pargstr (comment)
+end
+
+
+# WFT_ENCODE_DATE -- Procedure to encode the date in the form dd-mm-yy.
+
+procedure wft_encode_date (datestr, szdate)
+
+char	datestr[ARB]	# string containing the date
+int	szdate		# number of chars in the date string
+
+long	ctime
+int	time[LEN_TMSTRUCT]
+long	clktime(), lsttogmt()
+
+begin
+	ctime = clktime (long (0))
+	ctime = lsttogmt (ctime)
+	call brktime (ctime, time)
+
+	if (TM_YEAR(time) >= NEW_CENTURY) {
+	    call sprintf (datestr, szdate, "%04d-%02d-%02dT%02d:%02d:%02d")
+	        call pargi (TM_YEAR(time))
+	        call pargi (TM_MONTH(time))
+	        call pargi (TM_MDAY(time))
+	        call pargi (TM_HOUR(time))
+	        call pargi (TM_MIN(time))
+	        call pargi (TM_SEC(time))
+	} else {
+	    call sprintf (datestr, szdate, "%02d-%02d-%02d")
+	        call pargi (TM_MDAY(time))
+	        call pargi (TM_MONTH(time))
+	        call pargi (mod (TM_YEAR(time), CENTURY))
+	}
+end
+
+
+# WFT_FITS_CARD --  Procedure to fetch a single line from a string parameter
+# padding it to a maximum of maxcols characters and trimmimg the delim
+# character.
+
+procedure wft_fits_card (instr, ip, card, col_out, maxcols, delim)
+
+char	instr[ARB]	# input string
+int	ip		# input string pointer, updated at each call
+char	card[ARB]	# FITS card image
+int	col_out		# pointer to column in card
+int	maxcols		# maximum columns in card
+int	delim		# 1 character string delimiter
+
+int	op
+
+begin
+	op = col_out
+
+	# Copy string
+	while (op <= maxcols && instr[ip] != EOS && instr[ip] != delim) {
+	    card[op] = instr[ip]
+	    ip = ip + 1
+	    op = op + 1
+	}
+
+	# Fill remainder of card with blanks
+	while (op <= maxcols ) {
+	    card[op] = ' '
+	    op = op + 1
+	}
+
+	if (instr[ip] == delim)
+	    ip = ip + 1
+
+end
diff --git a/pkg/dataio/fits/fits_read.x b/pkg/dataio/fits/fits_read.x
new file mode 100644
index 00000000..f9b0e46c
--- /dev/null
+++ b/pkg/dataio/fits/fits_read.x
@@ -0,0 +1,469 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include <error.h>
+include <imhdr.h>
+include <fset.h>
+include <plset.h>
+include	"rfits.h"
+
+define	MAX_RANGES	100		# the maximum number of ranges
+
+# RFT_READ_FITZ -- Convert a FITS file. An EOT is signalled by returning EOF.
+
+int procedure rft_read_fitz (fitsfile, iraffile, pl, file_number)
+
+char	fitsfile[ARB]		# FITS file name
+char	iraffile[ARB]		# root IRAF file name
+pointer	pl			# pointer to the file/extensions list
+int	file_number		# the current file number
+
+bool	strne()
+int	fits_fd, stat, min_lenuserarea, ip, len_elist, oshort_header
+int	olong_header, ext_count, ext_number, max_extensions, naxes
+pointer	im, gim, sp, fits, axes, extensions, imname, gimname, gfname, str
+pointer	himname
+int	rft_read_header(), mtopen(), immap(), strlen(), envfind(), ctoi()
+int	rft_ext_skip()
+real	asumi()
+errchk	smark, sfree, salloc, rft_read_header, rft_read_image, rft_find_eof()
+errchk	rft_scan_file, mtopen, immap, imdelete, close, imunmap
+
+include	"rfits.com"
+
+begin
+	# Open input FITS data.
+	fits_fd = mtopen (fitsfile, READ_ONLY, 0)
+
+	# Allocate memory for the FITS data structure and initialize the file
+	# dependent  components of that structure.
+	call smark (sp)
+	call salloc (fits, LEN_FITS, TY_STRUCT)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+	call salloc (himname, SZ_FNAME, TY_CHAR)
+	call salloc (gimname, SZ_FNAME, TY_CHAR)
+	call salloc (gfname, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Initialize.
+	SIMPLE(fits) = NO
+	EXTEND(fits) = NO
+	GLOBALHDR(fits) = NO
+	gim = NULL
+	Memc[gfname] = EOS
+
+	# Determine the length of the user area.
+	if (envfind ("min_lenuserarea", Memc[imname], SZ_FNAME) > 0) {
+	    ip = 1
+	    if (ctoi (Memc[imname], ip, min_lenuserarea) <= 0)
+		min_lenuserarea = LEN_USERAREA
+	    else
+		min_lenuserarea = max (LEN_USERAREA, min_lenuserarea)
+	} else
+	    min_lenuserarea = LEN_USERAREA
+
+	# Store the current values of the header printing options.
+	olong_header = long_header
+	oshort_header = short_header
+
+	# Get the extensions list for a given line and count the number of
+	# extensions files.
+	call salloc (axes, 2, TY_INT)
+	call pl_gsize (pl, naxes, Memi[axes], stat)
+	max_extensions = Memi[axes+1]
+	call salloc (extensions, max_extensions, TY_INT)
+	Memi[axes] = 1
+	Memi[axes+1] = file_number
+	call pl_glpi (pl, Memi[axes], Memi[extensions], 1, max_extensions,
+	    PIX_SRC)
+	len_elist = nint (asumi (Memi[extensions], max_extensions))
+
+	# Loop over the extensions.
+	ext_count = 1; stat = BOF
+	do ext_number = 1, max_extensions {
+
+	    if (stat == EOF)
+		break
+	    if (Memi[extensions+ext_number-1] == 0)
+		next
+
+	    # Locate the next extension to be read.
+	    while (ext_count <= ext_number) {
+
+	        # Create the IRAF image header. If only a header listing is
+	        # desired or the image extension is to be skipped then map
+		# the scratch image onto DEV$NULL (faster than a real file).
+		# If more than one extension is to be read then append the
+		# extension number to the input name.
+
+	        if (make_image == NO || ext_count != ext_number) {
+	            call strcpy ("dev$null", Memc[imname], SZ_FNAME)
+		} else if (len_elist > 1 && ext_count == ext_number) {
+		    call sprintf (Memc[imname], SZ_FNAME, "%s%04d")
+			call pargstr (iraffile)
+			call pargi (ext_number - 1)
+		} else
+		    call strcpy (iraffile, Memc[imname], SZ_FNAME) 
+	        im = immap (Memc[imname], NEW_IMAGE, min_lenuserarea)
+		call strcpy (IM_HDRFILE(im), Memc[himname], SZ_FNAME)
+
+	        # Skip any extensions the user does not want. In order to do
+	        # this we must read the header to see how big the data array
+	        # to be skipped is.
+		if (ext_count != ext_number) {
+
+		    # Turn off header printing.
+		    long_header = NO
+		    short_header = NO
+
+		    # Decode the header and skip the data.
+	            iferr {
+	                stat = rft_read_header (fits_fd, fits, im, gim)
+			if (stat != EOF)
+			    stat = rft_ext_skip (fits_fd, fits, im)
+			if (stat == EOF) {
+			    if (ext_count == 1) {
+		                call printf ("File:  %s\n")
+		                    call pargstr (fitsfile)
+			    } else if (asumi(Memi[extensions],
+			        ext_count - 1) < 1.0) {
+		                call printf ("File:  %s\n")
+		                    call pargstr (fitsfile)
+			    }
+			    if (ext_count > 1) {
+		                call printf ("Extension: %d End of data\n")
+		                    call pargi (ext_count - 1)
+			    } else
+	                        call printf ("    End of data\n")
+			} else if (EXTEND(fits) == NO) {
+		            call printf ("File: %s\n")
+		                call pargstr (fitsfile)
+	                    call printf ("Extension: 1 End of data\n")
+			}
+	            } then {
+	                call flush (STDOUT)
+	                call erract (EA_WARN)
+	            }
+
+		    # Restore the default header printing values.
+		    long_header = olong_header
+		    short_header = oshort_header
+
+		# Read the extension the user specified. If the extension
+		# is not the primary data or IMAGE skip the data and
+		# continue.
+		} else {
+
+	            # Set up for printing a long or a short header.
+	            if (long_header == YES || short_header == YES) {
+			if (long_header == YES) {
+			    if (ext_number == 1) {
+				if (make_image == YES) {
+		                    call printf ("File: %s  Image: %s")
+		                        call pargstr (fitsfile)
+					#call pargstr (Memc[imname])
+					call pargstr (Memc[himname])
+				} else {
+		                    call printf ("File: %s")
+		                        call pargstr (fitsfile)
+				}
+			    } else if (asumi (Memi[extensions],
+			        ext_number -1) < 1.0) {
+				if (make_image == YES) {
+		                    call printf (
+				        "File: %s\nExtension: %d  Image: %s")
+		                        call pargstr (fitsfile)
+		                        call pargi (ext_number - 1)
+					#call pargstr (Memc[imname])
+					call pargstr (Memc[himname])
+				} else {
+		                    call printf ("File: %s  Extension: %d")
+		                        call pargstr (fitsfile)
+		                        call pargi (ext_number - 1)
+				}
+			    } else {
+				if (make_image == YES) {
+		                    call printf ("Extension: %d  Image: %s")
+		                        call pargi (ext_number - 1)
+					#call pargstr (Memc[imname])
+					call pargstr (Memc[himname])
+				} else {
+		                    call printf ("File: %s  Extension: %d")
+					call pargstr (fitsfile)
+		                        call pargi (ext_number - 1)
+				}
+			    }
+			} else {
+			    if (ext_number == 1) {
+		                call printf ("File:  %s ")
+		                    call pargstr (fitsfile)
+			    } else if (asumi (Memi[extensions],
+			        ext_number - 1) < 1.0) {
+		                call printf ("File:  %s\nExtension: %d ")
+		                    call pargstr (fitsfile)
+		                    call pargi (ext_number - 1)
+			    } else {
+		                call printf ("Extension: %d ")
+		                    call pargi (ext_number - 1)
+			    }
+			}
+	                if (long_header == YES)
+		            call printf ("\n")
+	            }
+	            call flush (STDOUT)
+
+	            # Read header.  EOT is signalled by an EOF status from
+	            # fits_read_header. Create an IRAF image if desired.
+
+	            iferr {
+	                stat = rft_read_header (fits_fd, fits, im, gim)
+	                if (stat == EOF) {
+	                    call printf ("End of data\n")
+	                } else if (make_image == YES) {
+			    if (XTENSION(fits) == EXT_PRIMARY ||
+			        XTENSION(fits) == EXT_IMAGE) {
+	                        call rft_read_image (fits_fd, fits, im)
+			    } else if (EXTEND(fits) == YES) {
+				stat = rft_ext_skip (fits_fd, fits, im)
+				if (stat == EOF)
+	                    	    call printf ("End of data\n")
+			    } else if (EXTEND(fits) == NO && fe > 0.0) {
+		                call rft_find_eof (fits_fd)
+			    }
+			} else {
+			    if (EXTEND(fits) == YES) {
+				stat = rft_ext_skip (fits_fd, fits, im)
+				if (stat == EOF)
+	                    	    call printf ("End of data\n")
+			    } else if (EXTEND(fits) == NO && fe > 0.0)
+		                call rft_scan_file (fits_fd, fits, im, fe)
+			}
+	            } then {
+	                call flush (STDOUT)
+	                call erract (EA_WARN)
+	            }
+		}
+
+
+		# Deal with the global header issue. Save the global header
+		# file name for possible future use.
+		if (GLOBALHDR(fits) == YES) {
+		    if (gim == NULL && XTENSION(fits) == EXT_PRIMARY) {
+		        call mktemp ("tmp$", Memc[gimname], SZ_FNAME)
+		        gim = immap (Memc[gimname], NEW_COPY, im)
+		        call strcpy (IRAFNAME(fits), Memc[gfname], SZ_FNAME)
+		    } else if (IRAFNAME(fits) == EOS)
+		        call strcpy (Memc[gfname], IRAFNAME(fits), SZ_FNAME)
+
+		}
+
+	        # Close the output image.
+	        call imunmap (im)
+
+	        # Optionally restore the old IRAF name.
+	        if (stat == EOF) {
+	            call imdelete (Memc[imname])
+		    break
+	        } else if (make_image == NO || ext_number != ext_count) {
+	            call imdelete (Memc[imname])
+		} else if (XTENSION(fits) != EXT_PRIMARY && XTENSION(fits) !=
+		    EXT_IMAGE) {
+	            call imdelete (Memc[imname])
+		    if (XTENSION(fits) != EXT_SPECIAL && ext_count ==
+		        ext_number)
+		        call printf ("    Skipping non-image data\n")
+	        } else if (old_name == YES && strlen (IRAFNAME(fits)) != 0) {
+	            iferr {
+		        call imgimage (IRAFNAME(fits), IRAFNAME(fits), SZ_FNAME)
+	                call imrename (Memc[imname], IRAFNAME(fits))
+	            } then {
+			if (len_elist > 1) {
+			    call sprintf (Memc[str], SZ_FNAME, ".%d")
+			        call pargi (ext_number - 1)
+			    call strcat (Memc[str], IRAFNAME(fits), SZ_FNAME)
+	                    iferr (call imrename (Memc[imname],
+			        IRAFNAME(fits))) {
+		                call printf (
+				    "    Cannot rename image %s to %s\n")
+		                    #call pargstr (Memc[imname])
+		                    call pargstr (Memc[himname])
+		                    call pargstr (IRAFNAME(fits))
+	                        call flush (STDOUT)
+	                        call erract (EA_WARN)
+			    } else {
+	                	call printf ("    Image %s renamed to %s\n")
+		            	    #call pargstr (Memc[imname])
+		                call pargstr (Memc[himname])
+		                call pargstr (IRAFNAME(fits))
+			    }
+			} else {
+		            call printf ("    Cannot rename image %s to %s\n")
+		                #call pargstr (Memc[imname])
+		                call pargstr (Memc[himname])
+		                call pargstr (IRAFNAME(fits))
+	                    call flush (STDOUT)
+	                    call erract (EA_WARN)
+			}
+	            } else {
+	                call printf ("    Image %s renamed to %s\n")
+		            #call pargstr (Memc[imname])
+		            call pargstr (Memc[himname])
+		            call pargstr (IRAFNAME(fits))
+	            }
+	        } else if (EXTEND(fits) == NO && strne (Memc[imname],
+		    iraffile)) {
+	            iferr {
+	                call imrename (Memc[imname], iraffile)
+	            } then {
+		        call printf ("    Cannot rename image %s to %s\n")
+		            #call pargstr (Memc[imname])
+		            call pargstr (Memc[himname])
+		            call pargstr (iraffile)
+	                call flush (STDOUT)
+	                call erract (EA_WARN)
+	            } else {
+	                call printf (
+			    "    No FITS extensions Image renamed to %s\n")
+		            #call pargstr (Memc[imname])
+		            call pargstr (iraffile)
+		    }
+		}
+
+		if (EXTEND(fits) == YES && XTENSION(fits) == EXT_PRIMARY &&
+		    len_elist == 1 && ext_number == 1) {
+	            if (short_header == YES || long_header == YES) {
+			if (long_header == NO)
+			    call printf ("    ")
+			call printf (
+			    "Warning: FITS extensions may be present\n")
+		    }
+		}
+	        if (long_header == YES)
+	            call printf ("\n")
+
+		ext_count = ext_count + 1
+		if (EXTEND(fits) == NO || XTENSION(fits) == EXT_SPECIAL)
+		    break
+	    }
+
+	    if (EXTEND(fits) == NO || XTENSION(fits) == EXT_SPECIAL)
+		break
+	}
+
+	if (gim != NULL) {
+	    call imunmap (gim)
+	    call imdelete (Memc[gimname])
+	}
+	call close (fits_fd)
+	call sfree (sp)
+
+	if (ext_count == 1)
+	    return (EOF)
+	else
+	    return (OK)
+end
+
+
+# RFT_FIND_EOF -- Read the FITS data file until EOF is reached.
+
+procedure rft_find_eof (fd)
+
+int	fd			# the FITS file descriptor
+
+int	szbuf
+pointer	sp, buf
+int	fstati(), read()
+errchk	read
+
+begin
+	# Scan through the file.
+	szbuf = fstati (fd, F_BUFSIZE)
+	call smark (sp)
+	call salloc (buf, szbuf, TY_CHAR)
+	while (read (fd, Memc[buf], szbuf) != EOF)
+	    ;
+	call sfree (sp)
+end
+
+
+# RFT_SCAN_FILE -- Determine whether it is more efficient to read the
+# entire file or to skip forward to the next file if the parameter
+# make_image was set to no.
+
+procedure  rft_scan_file (fd, fits, im, fe)
+
+int	fd			# the FITS file descriptor
+pointer	fits			# pointer to the FITS descriptor
+pointer	im			# pointer to the output image
+real	fe			# maximum file size in Kb for scan mode
+
+int	i, szbuf
+pointer	sp, buf
+real	file_size
+int	fstati(), read()
+errchk	read
+
+begin
+	# Compute the file size in Kb and return if it is bigger than fe.
+	file_size = 1.0
+	do i = 1, IM_NDIM(im)
+	    file_size = file_size * IM_LEN(im,i)
+	if (IM_NDIM(im) <= 0)
+	    file_size = 0.0
+	else
+	    file_size = file_size * abs (BITPIX(fits)) / FITS_BYTE / 1.0e3
+	if (file_size >= fe)
+	    return
+
+	# Scan through the file.
+	szbuf = fstati (fd, F_BUFSIZE)
+	call smark (sp)
+	call salloc (buf, szbuf, TY_CHAR)
+	while (read (fd, Memc[buf], szbuf) != EOF)
+	    ;
+	call sfree (sp)
+end
+
+
+# RFT_EXT_SKIP -- Compute the size of the data extension to be skipped
+# and do the skipping.
+
+int procedure rft_ext_skip (fits_fd, fits, im)
+
+int	fits_fd			# fits file descriptor
+pointer	fits			# pointer to the fits structure
+pointer	im			# pointer to the output image
+
+int	i, nbits, nblocks, sz_rec, blksize, stat
+pointer	buf
+int	fstati(), rft_getbuf()
+
+begin
+	# Compute the number of blocks to skip.
+	nbits = NAXISN(im,1)
+	do i = 2, NAXIS(im)
+	    nbits = nbits * NAXISN(im,i)
+	nbits = nbits + PCOUNT(fits)
+	nbits = abs (BITPIX(fits)) * GCOUNT(fits) * nbits
+	nblocks = int ((nbits + 23039) / 23040)
+
+	sz_rec = FITS_RECORD / SZB_CHAR
+	call malloc (buf, sz_rec, TY_CHAR)
+	blksize = fstati (fits_fd, F_SZBBLK)
+        if (mod (blksize, FITS_RECORD) == 0)
+            blksize = blksize / FITS_RECORD
+        else
+            blksize = 1
+
+	# Skip the blocks.
+	do i = 1, nblocks {
+	    stat = rft_getbuf (fits_fd, Memc[buf], sz_rec, blksize,
+	        NRECORDS(fits))
+	    if (stat == EOF)
+		break
+	}
+
+	call mfree (buf, TY_CHAR)
+
+	return (stat)
+end
diff --git a/pkg/dataio/fits/fits_rheader.x b/pkg/dataio/fits/fits_rheader.x
new file mode 100644
index 00000000..e15a3559
--- /dev/null
+++ b/pkg/dataio/fits/fits_rheader.x
@@ -0,0 +1,888 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <ctype.h>
+include <imhdr.h>
+include <imio.h>
+include <mach.h>
+include	"rfits.h"
+
+define	NEPSILON 10.0d0    # number of machine epsilon
+
+# RFT_READ_HEADER -- Read a FITS header.
+# If BSCALE and BZERO are different from 1.0 and 0.0  scale is set to true
+# otherwise scale is false.
+# EOT is detected by an EOF on the first read and EOF is returned to the calling
+# routine.  Any error is passed to the calling routine.
+
+int procedure rft_read_header (fits_fd, fits, im, gim)
+
+int	fits_fd			# FITS file descriptor
+pointer	fits			# FITS data structure
+pointer	im			# IRAF image descriptor
+pointer	gim			# IRAF global header image descriptor
+
+int	i, stat, nread, max_lenuser, fd_usr, ndiscard
+char	card[LEN_CARD+1], type_str[LEN_TYPESTR]
+int	rft_decode_card(), rft_init_read_pixels(), rft_read_pixels(), strmatch()
+int	stropen()
+errchk	rft_decode_card, rft_init_read_pixels, rft_read_pixels
+errchk	stropen, close
+
+include "rfits.com"
+
+begin
+	# Initialization.
+	XTENSION(fits) = EXT_PRIMARY
+	BITPIX(fits) = INDEFI
+	NAXIS(im) = 0
+	do i = 1, IM_MAXDIM
+	    IM_LEN(im,i) = 0
+	PCOUNT(fits) = 0
+	GCOUNT(fits) = 1
+	SCALE(fits) = NO
+	FITS_BSCALE(fits) = 1.0d0
+	FITS_BZERO(fits) = 0.0d0
+	BLANKS(fits) = NO
+	BLANK_VALUE(fits) = INDEFL
+	NRECORDS(fits) = 0
+	IRAFNAME(fits) = EOS
+	INHERIT(fits) = NO
+	ndiscard = 0
+	OBJECT(im) = EOS
+	UNKNOWN(im) = EOS
+	max_lenuser = (LEN_IMDES + IM_LENHDRMEM(im) - IMU) * SZ_STRUCT - 1
+
+	# The FITS header is character data in FITS_BYTE form. Open the
+	# header for reading. Open the user area which is a character
+	# string as a file.
+
+	i = rft_init_read_pixels (len_record, FITS_BYTE, LSBF, TY_CHAR)
+	fd_usr = stropen (UNKNOWN(im), max_lenuser, NEW_FILE)
+
+	# Loop until the END card is encountered.
+	nread = 0
+	repeat {
+
+	    # Read the card.
+	    i = rft_read_pixels (fits_fd, card, LEN_CARD, NRECORDS(fits), 1)
+	    card[LEN_CARD + 1] = '\n'
+	    card[LEN_CARD + 2] = EOS
+
+	    # Decode the card images.
+	    if ((i == EOF) && (nread == 0)) {
+		call close (fd_usr)
+		return (EOF)
+	    } else if ((nread == 0) && SIMPLE(fits) == NO &&
+	        strmatch (card, "^SIMPLE  ") == 0) {
+		call flush (STDOUT)
+		call close (fd_usr)
+		call error (30,
+		    "RFT_READ_HEADER: Not a FITS file (no SIMPLE keyword)")
+	    } else if ((nread == 0) && EXTEND(fits) == YES &&
+	        strmatch (card, "^XTENSION") == 0) {
+		XTENSION(fits) = EXT_SPECIAL
+		call flush (STDOUT)
+		call close (fd_usr)
+		call error (30,
+		"RFT_READ_HEADER: Not a FITS extension (no XTENSION keyword)")
+	    } else if (i != LEN_CARD) {
+		call close (fd_usr)
+	        call error (2, "RFT_READ_HEADER: Error reading FITS header")
+	    } else
+	        nread = nread + 1
+
+	    # Remove contaminating control characters and replace with blanks.
+	    call rft_control_to_blank (card, card, LEN_CARD)
+
+	    # Print FITS card images if long_header option specified.
+	    if (long_header == YES) {
+		call printf ("%-80.80s\n")
+		    call pargstr (card)
+	    }
+
+	    # Stat = YES if FITS END card is encountered.
+	    stat = rft_decode_card (fits, im, fd_usr, card, ndiscard)
+
+	} until (stat == YES)
+
+	# Check for the possibility of a global header.
+	if (NAXIS(im) == 0 && XTENSION(fits) == EXT_PRIMARY)
+	    GLOBALHDR(fits) = YES
+
+	# Set the output image pixel type.
+	call rft_set_image_pixtype (fits, im, FITS_BSCALE(fits),
+	    FITS_BZERO(fits))
+
+	# Copy the global header title and user area into the output image.
+	if (GLOBALHDR(fits) == YES) {
+	    if (XTENSION(fits) == EXT_IMAGE && INHERIT(fits) == YES &&
+	        gim != NULL) {
+		if (OBJECT(im) == EOS)
+		    call strcpy (OBJECT(gim), OBJECT(im), SZ_OBJECT)
+		call close (fd_usr)
+		fd_usr = stropen (UNKNOWN(im), max_lenuser, APPEND)
+	        call rft_gheader (im, gim, fd_usr, card, LEN_CARD, ndiscard,
+		    long_header)
+	    }
+	}
+
+	# Print optional short header.
+	if (short_header == YES && long_header == NO) {
+	    call printf ("%s  ")
+	    switch (XTENSION(fits)) {
+	    case EXT_PRIMARY:
+		call pargstr ("")
+	    case EXT_IMAGE:
+		call pargstr ("IMAGE")
+	    case EXT_TABLE:
+		call pargstr ("TABLE")
+	    case EXT_BINTABLE:
+		call pargstr ("BINTABLE")
+	    case EXT_UNKNOWN:
+		call pargstr ("UNKNOWN")
+	    default:
+		call pargstr ("UNDEFINED")
+	    }
+	    if (make_image == NO) {
+		if (old_name == YES) {
+	            call printf ("-> %s  ")
+		        call pargstr (IRAFNAME(fits))
+		}
+	    } else {
+		call printf ("-> %s  ")
+		    call pargstr (IM_HDRFILE(im))
+	    }
+	    call printf ("%-20.20s  ")
+	        call pargstr (OBJECT(im))
+	    call printf ("size=")
+	    if (NAXIS(im) == 0)
+		call printf ("0")
+	    else {
+	        do i = 1, NAXIS(im) {
+		    if (i == 1) {
+		        call printf ("%d")
+		            call pargl (NAXISN(im,i))
+		    } else {
+		        call printf ("x%d")
+		            call pargl (NAXISN(im,i))
+		    }
+	        }
+	    }
+	    call printf ("\n")
+	    if (XTENSION(fits) == EXT_PRIMARY || XTENSION(fits) == EXT_IMAGE) {
+	        call printf ("    bitpix=%d")
+		    call pargi (BITPIX(fits))
+	        if (SCALE(fits) == NO) {
+		    call printf ("  scaling=none")
+	        } else {
+	            call printf ("  bscale=%.7g  bzero=%.7g")
+		        call pargd (FITS_BSCALE(fits))
+			call pargd (FITS_BZERO(fits))
+	        }
+		call rft_typestring (PIXTYPE(im), type_str, LEN_TYPESTR)
+		call strlwr (type_str)
+	        call printf ("  pixtype=%s")
+		    call pargstr (type_str)
+	        call printf ("\n")
+	    }
+	}
+
+	# Let the user know if there is not enough space in the user area.
+	if (ndiscard > 0) {
+	    if (short_header == YES || long_header == YES) {
+		if (long_header == NO)
+		    call printf ("    ")
+	        call printf (
+	            "Warning: User area too small %d card images discarded\n")
+	        call pargi (ndiscard)
+	    }
+	    call rft_last_user (UNKNOWN(im), max_lenuser)
+	}
+
+	call close (fd_usr)
+	return (OK)
+end
+
+
+# RFT_CONTROL_TO_BLANK -- Replace an ACSII control characters in the
+# FITS card image with blanks.
+
+procedure rft_control_to_blank (incard, outcard, len_card)
+
+char	incard[ARB]		# the input FITS card image
+char	outcard[ARB]		# the output FITS card image
+int	len_card		# the length of the FITS card image
+
+int	i
+
+begin
+	for (i = 1; i <= len_card; i = i + 1) {
+	    if (IS_PRINT(incard[i]))
+	        outcard[i] = incard[i]
+	    else
+		outcard[i] = ' '
+	}
+end
+
+
+# RFT_DECODE_CARD -- Decode a FITS card and return YES when the END
+# card is encountered.  The keywords understood are given in rfits.h.
+
+int procedure rft_decode_card (fits, im, fd_usr, card, ndiscard)
+
+pointer	fits		# FITS data structure
+pointer	im		# IRAF image descriptor
+int	fd_usr		# file descriptor of user area
+char	card[ARB]	# FITS card
+int	ndiscard	# Number of cards for which no space available
+
+char	cval
+double	dval
+int	nchar, i, j, k, len
+pointer	sp, str, comment
+
+bool	rft_equald()
+int	strmatch(), ctoi(), ctol(), ctod(), cctoc(), rft_hms()
+errchk	putline
+
+include	"rfits.com"
+
+begin
+	call smark (sp)
+	call salloc (str, LEN_CARD, TY_CHAR)
+	call salloc (comment, SZ_LINE, TY_CHAR)
+
+	i = COL_VALUE
+	if (strmatch (card, "^END     ") != 0) {
+	    call sfree (sp)
+	    return(YES)
+	} else if (strmatch (card, "^SIMPLE  ") != 0) {
+	    if (SIMPLE(fits) == YES) {
+		if (short_header == YES || long_header == YES) {
+		    if (long_header == NO)
+			call printf ("    ")
+		    call printf ("Warning: Duplicate SIMPLE keyword ignored\n")
+		}
+	    } else {
+	        nchar = cctoc (card, i, cval)
+	        if (cval != 'T')
+		    call error (13, "RFT_DECODE_CARD: Non-standard FITS format")
+		else
+		    SIMPLE(fits) = YES
+	    }
+	} else if (strmatch (card, "^XTENSION") != 0) {
+	    call rft_get_fits_string (card, Memc[str], LEN_CARD)
+	    if (strmatch (Memc[str], "^IMAGE") != 0)
+		XTENSION(fits) = EXT_IMAGE
+	    else if (strmatch (Memc[str], "^TABLE") != 0)
+		XTENSION(fits) = EXT_TABLE
+	    else if (strmatch (Memc[str], "^BINTABLE") != 0)
+		XTENSION(fits) = EXT_BINTABLE
+	    else
+		XTENSION(fits) = EXT_UNKNOWN
+	} else if (strmatch (card, "^BITPIX  ") != 0) {
+	    if (! IS_INDEFI(BITPIX(fits))) {
+		if (short_header == YES || long_header == YES) {
+		    if (long_header == NO)
+			call printf ("    ")
+		    call printf ("Warning: Duplicate BITPIX keyword ignored\n")
+		}
+	    } else
+	        nchar = ctoi (card, i, BITPIX(fits))
+	} else if (strmatch (card, "^NAXIS   ") != 0) {
+	    if (NAXIS(im) != 0) {
+		if (short_header == YES || long_header == YES) {
+		    if (long_header == NO)
+			call printf ("    ")
+		    call printf ("Warning: Duplicate NAXIS keyword ignored\n")
+		}
+	    } else
+	        nchar = ctoi (card, i, NAXIS(im))
+	    if (NAXIS(im) > IM_MAXDIM)
+		call error (5, "RFT_DECODE_CARD: FITS NAXIS too large")
+	} else if (strmatch (card, "^NAXIS") != 0) {
+	    k = strmatch (card, "^NAXIS")
+	    nchar = ctoi (card, k, j)
+	    if (NAXISN(im,j) != 0) {
+		if (short_header == YES || long_header == YES) {
+		    if (long_header == NO)
+			call printf ("    ")
+		    call printf ("Warning: Duplicate NAXIS%d keyword ignored\n")
+		        call pargi (j)
+		}
+	    } else
+	        nchar = ctol (card, i, NAXISN(im, j))
+	} else if (strmatch (card, "^GROUPS  ") != 0) {
+	    nchar = cctoc (card, i, cval)
+	    if (cval == 'T') {
+		NAXIS(im) = 0
+		call error (6, "RFT_DECODE_CARD: Group data not implemented")
+	    }
+	} else if (strmatch (card, "^EXTEND  ") != 0) {
+	    if (EXTEND(fits) == YES) {
+		if (short_header == YES || long_header == YES) {
+		    if (long_header == NO)
+			call printf ("    ")
+		    call printf ("Warning: Duplicate EXTEND keyword ignored\n")
+		}
+	    } else {
+	        nchar = cctoc (card, i, cval)
+	        if (cval == 'T')
+		    EXTEND(fits) = YES
+	    }
+	} else if (strmatch (card, "^INHERIT ") != 0) {
+	    if (INHERIT(fits) == YES) {
+		if (short_header == YES || long_header == YES) {
+		    if (long_header == NO)
+			call printf ("    ")
+		    call printf ("Warning: Duplicate INHERIT keyword ignored\n")
+		}
+	    } else {
+	        nchar = cctoc (card, i, cval)
+	        if (cval == 'T')
+		    INHERIT(fits) = YES
+	    }
+	} else if (strmatch (card, "^PCOUNT  ") != 0) {
+	    nchar = ctoi (card, i, PCOUNT(fits))
+	    if (nchar <= 0)
+		PCOUNT(fits) = 0
+	} else if (strmatch (card, "^GCOUNT  ") != 0) {
+	    nchar = ctoi (card, i, GCOUNT(fits))
+	    if (nchar <= 0)
+		GCOUNT(fits) = 1
+	#} else if (strmatch (card, "^TABLES  ") != 0) {
+	    #nchar = ctoi (card, i, ival)
+	    #if (ival > 0)
+		#call printf ("Warning: FITS special records not decoded\n")
+	} else if (strmatch (card, "^BSCALE  ") != 0) {
+	    nchar = ctod (card, i, dval)
+	    if (nchar > 0)
+	        FITS_BSCALE(fits) = dval
+	    else if (short_header == YES || long_header == YES) {
+		if (long_header == NO)
+		    call printf ("    ")
+		call printf ("Warning: Error decoding BSCALE, BSCALE=1.0\n")
+	    }
+	    if (! rft_equald (dval, 1.0d0) && (scale == YES))
+		SCALE(fits) = YES
+	} else if (strmatch (card, "^BZERO   ") != 0) {
+	    nchar = ctod (card, i, dval)
+	    if (nchar > 0)
+	        FITS_BZERO(fits) = dval
+	    else if (short_header == YES || long_header == YES) {
+		if (long_header == NO)
+		    call printf ("    ")
+		call printf ("Warning: Error decoding BZERO, BZERO=0.0\n")
+	    }
+	    if (! rft_equald (dval, 0.0d0) && (scale == YES))
+		SCALE(fits) = YES
+	} else if (strmatch (card, "^BLANK   ") != 0) {
+	    BLANKS(fits) = YES
+	    nchar = ctol (card, i, BLANK_VALUE(fits))
+	} else if (strmatch (card, "^OBJECT  ") != 0) {
+	    call rft_get_fits_string (card, OBJECT(im), SZ_OBJECT)
+	} else if (strmatch (card, "^IRAFNAME") != 0) {
+	    call rft_get_fits_string (card, IRAFNAME(fits), SZ_FNAME)
+	} else if (strmatch (card, "^FILENAME") != 0) {
+	    if (IRAFNAME(fits) == EOS)
+	        call rft_get_fits_string (card, IRAFNAME(fits), SZ_FNAME)
+	} else if (strmatch (card, "^EXTNAME ") != 0) {
+	    if (XTENSION(fits) != EXT_PRIMARY && XTENSION(fits) != EXT_IMAGE)
+	        call rft_get_fits_string (card, IRAFNAME(fits), SZ_FNAME)
+	} else if (strmatch (card, "^EXTVER  ") != 0) {
+	    # Filter this quantitity out and ignore it for now.
+	    ;
+	} else if (strmatch (card, "^ORIGIN  ") != 0) {
+	    call rft_trim_card (card, card, LEN_CARD)
+	    call strcat (card[i], HISTORY(im), SZ_HISTORY)
+	} else if (strmatch (card, "^DATE    ") !=  0) {
+	    call rft_trim_card (card, card, LEN_CARD)
+	    call strcat (card[i], HISTORY(im), SZ_HISTORY)
+	} else if (strmatch (card, "^IRAF-TLM") !=  0) {
+	    call rft_trim_card (card, card, LEN_CARD)
+	    call strcat (card[i], HISTORY(im), SZ_HISTORY)
+	#} else if (strmatch (card, "^HISTORY ") != 0) {
+	    #call rft_trim_card (card, card, LEN_CARD)
+	    #call strcat (card[i - 2], HISTORY(im), SZ_HISTORY)
+	} else if (strmatch (card, "^UT      ") != 0) {
+	    len = rft_hms (card, Memc[str], Memc[comment], LEN_CARD)
+	    if (len > 0) {
+		call wft_encodec ("UT", Memc[str], len, card, Memc[comment])
+		card[LEN_CARD+1] = '\n'
+		card[LEN_CARD+2] = EOS
+	    }
+	    if (ndiscard > 1)
+		ndiscard = ndiscard + 1
+	    else {
+	        iferr (call putline (fd_usr, card))
+		    ndiscard = ndiscard + 1
+	    }
+	} else if (strmatch (card, "^ZD      ") != 0) {
+	    len = rft_hms (card, Memc[str], Memc[comment], LEN_CARD)
+	    if (len > 0) {
+		call wft_encodec ("ZD", Memc[str], len, card, Memc[comment])
+		card[LEN_CARD+1] = '\n'
+		card[LEN_CARD+2] = EOS
+	    }
+	    if (ndiscard > 1)
+		ndiscard = ndiscard + 1
+	    else {
+	        iferr (call putline (fd_usr, card))
+		    ndiscard = ndiscard + 1
+	    }
+	} else if (strmatch (card, "^ST      ") != 0) {
+	    len = rft_hms (card, Memc[str], Memc[comment], LEN_CARD)
+	    if (len > 0) {
+		call wft_encodec ("ST", Memc[str], len, card, Memc[comment])
+		card[LEN_CARD+1] = '\n'
+		card[LEN_CARD+2] = EOS
+	    }
+	    if (ndiscard > 1)
+		ndiscard = ndiscard + 1
+	    else {
+	        iferr (call putline (fd_usr, card))
+		    ndiscard = ndiscard + 1
+	    }
+	} else if (strmatch (card, "^RA      ") != 0) {
+	    len = rft_hms (card, Memc[str], Memc[comment], LEN_CARD)
+	    if (len > 0) {
+		call wft_encodec ("RA", Memc[str], len, card, Memc[comment])
+		card[LEN_CARD+1] = '\n'
+		card[LEN_CARD+2] = EOS
+	    }
+	    if (ndiscard > 1)
+		ndiscard = ndiscard + 1
+	    else {
+	        iferr (call putline (fd_usr, card))
+		    ndiscard = ndiscard + 1
+	    }
+	} else if (strmatch (card, "^DEC     ") != 0) {
+	    len = rft_hms (card, Memc[str], Memc[comment], LEN_CARD)
+	    if (len > 0) {
+		call wft_encodec ("DEC", Memc[str], len, card, Memc[comment])
+		card[LEN_CARD+1] = '\n'
+		card[LEN_CARD+2] = EOS
+	    }
+	    if (ndiscard > 1)
+		ndiscard = ndiscard + 1
+	    else {
+	        iferr (call putline (fd_usr, card))
+		    ndiscard = ndiscard + 1
+	    }
+	} else {
+	    if (ndiscard > 1)
+		ndiscard = ndiscard + 1
+	    else {
+	        iferr (call putline (fd_usr, card))
+		    ndiscard = ndiscard + 1
+	    }
+	}
+
+	call sfree (sp)
+
+	return (NO)
+
+end
+
+
+# RFT_HMS -- Procedure to decode a FITS HMS card from the mountain.
+
+int procedure rft_hms (card, str, comment, maxch)
+
+char	card[ARB]		# FITS card
+char	str[ARB]		# string
+char	comment[ARB]		# comment string
+int	maxch			# maximum number of characters
+
+char	colon, minus
+int	ip, nchar, fst, lst, deg, min
+real	sec
+int	stridx(), strldx(), strlen(), ctoi(), ctor()
+
+begin
+	# Return if not a FITS string parameter.
+	if (card[COL_VALUE] != '\'')
+	    return (0)
+
+	# Set up key characters.
+	colon = ':'
+	minus = '-'
+
+	# Get the FITS string.
+	call rft_get_fits_string (card, str, maxch)
+
+	# Get the comment string.
+	call rft_get_comment (card, comment, maxch)
+
+	# Test for blank string and for 2 colon delimiters.
+	if (str[1] == EOS)
+	    return (0)
+	fst = stridx (colon, str)
+	if (fst == 0)
+	    return (0)
+	lst = strldx (colon, str)
+	if (lst == 0)
+	    return (0)
+	if (fst == lst)
+	    return (0)
+
+	# Decode the degrees field.
+	ip = 1
+	while (IS_WHITE(str[ip]))
+	    ip = ip + 1
+	if (str[ip] == '+' || str[ip] == '-')
+	    ip = ip + 1
+	nchar = ctoi (str, ip, deg)
+	if (nchar == 0)
+	    deg = 0
+
+	# Decode the minutes field.
+	ip = fst + 1
+	while (IS_WHITE(str[ip]))
+	    ip = ip + 1
+	if (str[ip] == '+' || str[ip] == '-')
+	    ip = ip + 1
+	nchar = ctoi (str, ip, min)
+	if (nchar == 0)
+	    min = 0
+
+	# Decode the seconds field.
+	ip = lst + 1
+	while (IS_WHITE(str[ip]))
+	    ip = ip + 1
+	if (str[ip] == '+' || str[ip] == '-')
+	    ip = ip + 1
+	nchar = ctor (str, ip, sec)
+	if (nchar == 0)
+	    sec = 0.0
+
+	# Reformat the HMS card.
+	if (stridx (minus, str) > 0 || deg < 0 || min < 0 || sec < 0.0) {
+	    call sprintf (str, maxch, "%c%d:%02d:%05.2f")
+		call pargc (minus)
+	        call pargi (abs (deg))
+	        call pargi (abs (min))
+	        call pargr (abs (sec))
+	} else {
+	    call sprintf (str, maxch, "%2d:%02d:%05.2f")
+	        call pargi (deg)
+	        call pargi (abs (min))
+	        call pargr (abs (sec))
+	}
+
+	return (strlen (str))
+end
+
+
+# RFT_GET_COMMENT -- Extract the comment field from a FITS card.
+
+procedure rft_get_comment (card, comment, maxch)
+
+char	card[ARB]		# FITS card
+char	comment[ARB]		# comment string
+int	maxch			# maximum number of characters
+
+int	istart, j
+
+begin
+	istart = 0
+	for (j = LEN_CARD; (j >= 1) && (card[j] != '\''); j = j - 1) {
+	    if (card[j] == '/') {
+		for (istart = j + 1; IS_WHITE(card[istart]) && istart <=
+		    LEN_CARD; istart = istart + 1)
+		    ;
+		break
+	    }
+	}
+
+	if (istart == 0)
+	    comment[1] = EOS
+	else
+	    call strcpy (card[istart], comment, LEN_CARD - istart + 1 )
+end
+
+
+# RFT_GET_FITS_STRING -- Extract a string from a FITS card and trim trailing
+# blanks. The EOS is marked by either ', /, or the end of the card.
+# There may be an optional opening ' (FITS standard).
+
+procedure rft_get_fits_string (card, str, maxchar)
+
+char	card[ARB]		# FITS card
+char	str[ARB]		# FITS string
+int	maxchar			# maximum number of characters
+
+int	j, istart, nchar
+
+begin
+	# Check for opening quote
+	for (istart = COL_VALUE; istart <= LEN_CARD && card[istart] != '\'';
+	    istart = istart + 1)
+	    ;
+	istart = istart + 1
+
+	# Check for closing quote.
+	for (j = istart; (j<LEN_CARD)&&(card[j]!='\''); j = j + 1)
+	    ;
+	for (j = j - 1; (j >= istart) && (card[j] == ' '); j = j - 1)
+	    ;
+	nchar = min (maxchar, j - istart + 1)
+
+	# Copy the string.
+	if (nchar <= 0)
+	    str[1] = EOS
+	else
+	    call strcpy (card[istart], str, nchar)
+end
+
+
+# RFT_EQUALD -- Procedure to compare two double precision numbers for equality
+# to within the machine precision for doubles.
+
+bool procedure rft_equald (x, y)
+
+double	x, y		# the two numbers to be compared for equality
+
+int	ex, ey
+double	x1, x2, normed_x, normed_y
+
+begin
+	if (x == y)
+	    return (true)
+
+	call rft_normd (x, normed_x, ex)
+	call rft_normd (y, normed_y, ey)
+
+	if (ex != ey)
+	    return (false)
+	else {
+	    x1 = 1.0d0 + abs (normed_x - normed_y)
+	    x2 = 1.0d0 + NEPSILON * EPSILOND
+	    return (x1 <= x2)
+	}
+end
+
+
+# RFT_NORMED -- Normalize a double precision number x to the value normed_x,
+# in the range [1-10]. Expon is returned such that x = normed_x *
+# (10.0d0 ** expon).
+
+procedure rft_normd (x, normed_x, expon)
+
+double	x			# number to be normailized
+double	normed_x		# normalized number
+int	expon			# exponent
+
+double	ax
+
+begin
+	ax = abs (x)
+	expon = 0
+
+	if (ax > 0) {
+	    while (ax < (1.0d0 - NEPSILON * EPSILOND)) {
+		ax = ax * 10.0d0
+		expon = expon - 1
+	    }
+
+	    while (ax >= (10.0d0 - NEPSILON * EPSILOND)) {
+		ax = ax / 10.0d0
+		expon = expon + 1
+	    }
+	}
+
+	if (x < 0)
+	    normed_x = -ax
+	else
+	    normed_x = ax
+end
+
+
+# RFT_TRIM_CARD -- Procedure to trim trailing whitespace from the card
+
+procedure rft_trim_card (incard, outcard, maxch)
+
+char	incard[ARB]		# input FITS card image 
+char	outcard[ARB]		# output FITS card
+int	maxch			# maximum size of card
+
+int	ip
+
+begin
+	ip = maxch
+	while (incard[ip] == ' ' || incard[ip] == '\t' || incard[ip] == '\0')
+	    ip = ip - 1
+	call amovc (incard, outcard, ip)
+	outcard[ip+1] = '\n'
+	outcard[ip+2] = EOS
+end
+
+
+# RFT_LAST_CARD -- Remove a partially written card from the data base
+
+procedure rft_last_user (user, maxch) 
+
+char	user[ARB]	# user area
+int	maxch		# maximum number of characters
+
+int	ip
+
+begin
+	ip = maxch
+	while (user[ip] != '\n')
+	    ip = ip - 1
+	user[ip+1] = EOS
+end
+
+
+# RFT_SET_IMAGE_PIXTYPE -- Set remaining header fields not set in
+# rft_read_header.
+
+procedure rft_set_image_pixtype (fits, im, bscale, bzero)
+
+pointer fits            # FITS data structure
+pointer im              # IRAF image pointer
+double	bscale		# FITS scaling parameter
+double	bzero		# FITS offset parameter
+
+bool	rft_equald()
+include "rfits.com"
+
+begin
+        # Determine data type from BITPIX if user data type not specified.
+
+        if (data_type == ERR) {
+            if (BITPIX(fits) < 0) {
+                if (abs (BITPIX(fits)) <= (SZ_REAL * SZB_CHAR * NBITS_BYTE))
+                    PIXTYPE(im) = TY_REAL
+                else
+                    PIXTYPE(im) = TY_DOUBLE
+            } else if (SCALE(fits) == YES) {
+		if (rft_equald (bscale, 1.0d0)) {
+		    if (rft_equald (bzero / 32768.0d0, 1.0d0))
+                        PIXTYPE(im) = TY_USHORT
+		    else
+                        PIXTYPE(im) = TY_REAL
+		} else
+                    PIXTYPE(im) = TY_REAL
+            } else {
+                if (BITPIX(fits) <= (SZ_SHORT * SZB_CHAR * NBITS_BYTE))
+                    PIXTYPE(im) = TY_SHORT
+                else
+                    PIXTYPE(im) = TY_LONG
+            }
+
+        } else
+            PIXTYPE(im) = data_type
+end
+
+
+# Copy the global header into the output image header.
+
+procedure rft_gheader (im, gim, fd_usr, card, len_card, ndiscard, long_header)
+
+pointer	im		# IRAF image header descriptor
+pointer	gim		# IRAF global image header descriptor
+int	fd_usr		# IRAF image header user area
+char	card[ARB]	# FITS card
+int	len_card	# length of FITS card
+int	ndiscard	# number of cards discarded
+int	long_header	# print the long header
+
+int	ngcards, gim_lenuser, ninherit, count
+pointer	sp, indices, idb_gim, grp, irp
+bool	streq()
+int	strlen(), idb_nextcard(), idb_find()
+pointer	idb_open()
+errchk	putline()
+
+begin
+	# Initialize.
+	call smark (sp)
+	ngcards = strlen (UNKNOWN(gim)) / (len_card + 1)
+	call salloc (indices, ngcards, TY_INT)
+
+	# Mark the global header cards which are to be inherited. These
+	# include all COMMENT, HISTORY, and BLANK cards, plus all those
+	# cards which do not already have values in the extension header.
+	count = 0
+	idb_gim = idb_open (gim, gim_lenuser)
+	while (idb_nextcard (idb_gim, grp) != EOF) {
+	    if (count >= ngcards)
+		break
+	    call strcpy (Memc[grp], card, 8)
+	    if (streq (card, "COMMENT "))
+		Memi[indices+count] = YES
+	    else if (streq (card, "HISTORY "))
+		Memi[indices+count] = YES
+	    else if (streq (card, "        "))
+		Memi[indices+count] = YES
+	    else if (idb_find (im, card, irp) > 0)
+		Memi[indices+count] = NO
+	    else 
+		Memi[indices+count] = YES
+	    count = count + 1
+	}
+	call idb_close (idb_gim)
+
+	# Open the global header image user area and loop through the cards.
+	ninherit = 0
+	count = 0
+	idb_gim = idb_open (gim, gim_lenuser)
+	while (idb_nextcard (idb_gim, grp) != EOF) {
+	    if (Memi[indices+count] == YES) {
+	        call strcpy (Memc[grp], card, len_card)
+	        card[len_card+1] = '\n'
+	        card[len_card+2] = EOS
+	        if (ndiscard > 1)
+		    ndiscard = ndiscard + 1
+	        else {
+	            iferr (call putline (fd_usr, card))
+		        ndiscard = ndiscard + 1
+		    else
+		        ninherit = ninherit + 1
+	        }
+	    }
+	    count = count + 1
+	}
+	call idb_close (idb_gim)
+
+	if (long_header == YES) {
+	    call printf ("%d global header keywords were inherited\n")
+		call pargi (ninherit)
+	}
+
+	call sfree (sp)
+end
+
+
+# RFT_TYPESTRING -- Procedure to set the iraf datatype keyword.
+
+procedure rft_typestring (data_type, type_str, maxch)
+
+int     data_type	# the IRAF data type
+char    type_str[ARB]	# the output IRAF type string
+int	maxch		# maximum size of the type string
+
+begin
+        switch (data_type) {
+        case TY_SHORT:
+            call strcpy ("SHORT", type_str, maxch)
+        case TY_USHORT:
+            call strcpy ("USHORT", type_str, maxch)
+        case TY_INT:
+            call strcpy ("INTEGER", type_str, maxch)
+        case TY_LONG:
+            call strcpy ("LONG", type_str, maxch)
+        case TY_REAL:
+            call strcpy ("REAL", type_str, maxch)
+        case TY_DOUBLE:
+            call strcpy ("DOUBLE", type_str, maxch)
+        case TY_COMPLEX:
+            call strcpy ("COMPLEX", type_str, maxch)
+        default:
+            call strcpy ("UNKNOWN", type_str, maxch)
+        }
+end
+
+
diff --git a/pkg/dataio/fits/fits_rimage.x b/pkg/dataio/fits/fits_rimage.x
new file mode 100644
index 00000000..9994c1d4
--- /dev/null
+++ b/pkg/dataio/fits/fits_rimage.x
@@ -0,0 +1,557 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include <mach.h>
+include <fset.h>
+include	"rfits.h"
+
+# RFT_READ_IMAGE --  Convert FITS image pixels to IRAF image pixels.
+
+procedure rft_read_image (fits_fd, fits, im)
+
+int	fits_fd		# FITS file descriptor
+pointer	fits		# FITS data structure
+pointer	im		# IRAF image descriptor
+
+int	i, npix, npix_record, blksize, ndummy
+long	v[IM_MAXDIM], nlines, il
+pointer	tempbuf, buf
+real	linemax, linemin, lirafmin, lirafmax
+double	dblank
+
+long	clktime()
+int	fstati(), rft_init_read_pixels(), rft_read_pixels()
+
+errchk	malloc, mfree, rft_init_read_pixels, rft_read_pixels, rft_lscale_pix
+errchk	rft_lchange_pix, rft_rchange_pix, rfit_dchange_pix, rft_put_image_line
+errchk	rft_pix_limits, rft_rscale_pix, rft_dscale_pix
+
+include	"rfits.com"
+
+begin
+	# No pixel file was created.
+	if (NAXIS(im) == 0) {
+	    if (short_header == YES || long_header == YES) {
+		if (long_header == NO)
+		    call printf ("    ")
+	        call printf ("Warning: No pixel file created\n")
+	    }
+	    return
+	}
+
+	# Compute the number of columns and lines in the image.
+	npix = NAXISN(im, 1)
+	nlines = 1
+	do i = 2, NAXIS(im)
+	    nlines = nlines * NAXISN(im, i)
+	lirafmax = -MAX_REAL
+	lirafmin = MAX_REAL
+
+	# Compute the number of pixels per record and the number of records
+	# per output block.
+
+	npix_record = len_record * FITS_BYTE / abs (BITPIX(fits))
+	blksize = fstati (fits_fd, F_SZBBLK)
+	if (mod (blksize, FITS_RECORD) == 0)
+	    blksize = blksize / FITS_RECORD
+	else
+	    blksize = 1
+
+	# FITS data is converted to type  LONG, REAL or DOUBLE. If BITPIX is
+	# not one of the MII types then rft_read_pixels returns an ERROR.
+
+	call amovkl (long(1), v, IM_MAXDIM)
+	switch (BITPIX(fits)) {
+	case FITS_REAL:
+
+	    # Allocate temporary space.
+	    call malloc (tempbuf, npix, TY_REAL)
+
+	    # Initialize the read.
+	    i = rft_init_read_pixels (npix_record, BITPIX(fits), LSBF, TY_REAL)
+
+	    # Turn on the ieee NaN mapping.
+	    call ieesnanr (blank)
+	    call ieemapr (YES, NO)
+	    #call ieezstatr ()
+	    NBPIX(im) = 0
+
+	    # Allocate the space for the output line, read in the image
+	    # line, convert from the ieee to native format, and compute the
+	    # minimum and maximum.
+
+	    do il = 1, nlines {
+	        call rft_put_image_line (im, buf, v, PIXTYPE(im))
+	        if (rft_read_pixels (fits_fd, Memr[tempbuf], npix,
+	            NRECORDS(fits), blksize) != npix)
+		    call printf ("Error reading FITS data\n")
+	        if (SCALE(fits) == YES)
+	            call rft_rscale_pix (Memr[tempbuf], buf, npix,
+		        FITS_BSCALE(fits), FITS_BZERO(fits), PIXTYPE(im))
+		else
+		    call rft_rchange_pix (Memr[tempbuf], buf, npix, PIXTYPE(im))
+	        call rft_pix_limits (buf, npix, PIXTYPE(im), linemin, linemax)
+	        lirafmax = max (lirafmax, linemax)
+	        lirafmin = min (lirafmin, linemin)
+	    }
+
+	    # Set the number of bad pixels.
+	    call ieestatr (NBPIX(im), ndummy)
+
+	    # Free space.
+	    call mfree (tempbuf, TY_REAL)
+
+	case FITS_DOUBLE:
+
+	    # Allocate temporary space.
+	    call malloc (tempbuf, npix, TY_DOUBLE)
+
+	    # Initialize the read.
+	    i = rft_init_read_pixels (npix_record, BITPIX(fits), LSBF,
+	        TY_DOUBLE)
+
+	    # Turn on the ieee NaN mapping.
+	    dblank = blank
+	    call ieesnand (dblank)
+	    call ieemapd (YES, NO)
+	    #call ieezstatd ()
+	    NBPIX(im) = 0
+
+	    # Allocate the space for the output line, read in the image
+	    # line, convert from the ieee to native format, and compute the
+	    # minimum and maximum.
+
+	    do il = 1, nlines {
+	        call rft_put_image_line (im, buf, v, PIXTYPE(im))
+	        if (rft_read_pixels (fits_fd, Memd[tempbuf], npix,
+	            NRECORDS(fits), blksize) != npix)
+		    call printf ("Error reading FITS data\n")
+	        if (SCALE(fits) == YES)
+	            call rft_dscale_pix (Memd[tempbuf], buf, npix,
+		        FITS_BSCALE(fits), FITS_BZERO(fits), PIXTYPE(im))
+		else
+		    call rft_dchange_pix (Memd[tempbuf], buf, npix, PIXTYPE(im))
+	        call rft_pix_limits (buf, npix, PIXTYPE(im), linemin, linemax)
+		if (IS_INDEFR(linemax))
+		    lirafmax = INDEFR
+		else
+	            lirafmax = max (lirafmax, linemax)
+		if (IS_INDEFR(linemin))
+		    lirafmin = INDEFR
+		else
+	            lirafmin = min (lirafmin, linemin)
+	    }
+
+	    # Set the number of bad pixels.
+	    call ieestatd (NBPIX(im), ndummy)
+
+	    # Free space.
+	    call mfree (tempbuf, TY_DOUBLE)
+
+	default:
+
+	    # Allocate the required space.
+	    call malloc (tempbuf, npix, TY_LONG)
+
+	    # Allocate the space for the output line, read in the image
+	    # line, convert from the ieee to native format, and compute the
+	    # minimum and maximum.
+
+	    i = rft_init_read_pixels (npix_record, BITPIX(fits), LSBF, TY_LONG)
+	    do il = 1, nlines {
+	        call rft_put_image_line (im, buf, v, PIXTYPE(im))
+	        if (rft_read_pixels (fits_fd, Meml[tempbuf], npix,
+	            NRECORDS(fits), blksize) != npix)
+		    call printf ("Error reading FITS data\n")
+	        if (SCALE(fits) == YES)
+	            call rft_lscale_pix (Meml[tempbuf], buf, npix,
+		        FITS_BSCALE(fits), FITS_BZERO(fits), PIXTYPE(im))
+	        else
+		    call rft_lchange_pix (Meml[tempbuf], buf, npix, PIXTYPE(im))
+	        if (BLANKS(fits) == YES)
+	            call rft_map_blanks (Meml[tempbuf], buf, npix, PIXTYPE(im),
+		        BLANK_VALUE(fits), blank, NBPIX(im))
+	        call rft_pix_limits (buf, npix, PIXTYPE(im), linemin, linemax)
+	        lirafmax = max (lirafmax, linemax)
+	        lirafmin = min (lirafmin, linemin)
+	    }
+
+	    # Free space.
+	    call mfree (tempbuf, TY_LONG)
+	}
+
+	IRAFMIN(im) = lirafmin
+	IRAFMAX(im) = lirafmax
+	LIMTIME(im) = clktime (long(0))
+
+	if (short_header == YES || long_header == YES) {
+	    if (NBPIX (im) != 0) {
+		if (long_header == NO)
+		    call printf ("    ")
+	        call printf ("Warning: %d bad pixels replaced in image\n")
+		    call pargl (NBPIX (im))
+	    }
+	    if (IS_INDEFR(lirafmax) || lirafmax > MAX_REAL) {
+		if (long_header == NO)
+		    call printf ("    ")
+	        call printf ("Warning: image contains pixel values > %g\n")
+		    call pargr (MAX_REAL)
+	    }
+	    if (IS_INDEFR(lirafmin) || lirafmin < -MAX_REAL) {
+		if (long_header == NO)
+		    call printf ("    ")
+	        call printf ("Warning: image contains pixel values < %g\n")
+		    call pargr (-MAX_REAL)
+	    }
+	}
+end
+
+
+# RFT_PUT_IMAGE_LINE -- Procedure to output an image line to and IRAF file.
+
+procedure rft_put_image_line (im, buf, v, data_type)
+
+pointer	im			# IRAF image descriptor
+pointer	buf			# Pointer to output image line
+long	v[ARB]			# imio pointer
+int	data_type		# output pixel type
+
+int	impnll(), impnlr(), impnld(), impnlx()
+errchk	impnll, impnlr, impnld, impnlx
+
+begin
+	switch (data_type) {
+	case TY_SHORT, TY_INT, TY_USHORT, TY_LONG:
+	    if (impnll (im, buf, v) == EOF)
+		call error (3, "RFT_PUT_IMAGE_LINE: Error writing FITS data")
+	case TY_REAL:
+	    if (impnlr (im, buf, v) == EOF)
+		call error (3, "RFT_PUT_IMAGE_LINE: Error writing FITS data")
+	case TY_DOUBLE:
+	    if (impnld (im, buf, v) == EOF)
+		call error (3, "RFT_PUT_IMAGE_LINE: Error writing FITS data")
+	case TY_COMPLEX:
+	    if (impnlx (im, buf, v) == EOF)
+		call error (3, "RFT_PUT_IMAGE_LINE: Error writing FITS data")
+	default:
+	    call error (10, "RFT_PUT_IMAGE_LINE: Unsupported IRAF image type")
+	}
+end
+
+
+# RFT_RSCALE_PIX -- Procedure to convert an IRAF image line from type real
+# to the requested output data type with optional scaling using the
+# FITS parameters BSCALE and BZERO.
+
+procedure rft_rscale_pix (inbuf, outbuf, npix, bscale, bzero, data_type)
+
+real	inbuf[ARB]		# buffer of FITS integers
+pointer	outbuf			# pointer to output image line
+int	npix			# number of pixels
+double	bscale, bzero		# FITS bscale and bzero
+int	data_type		# IRAF image pixel type
+
+errchk	altmdr, achtrl, amovr, achtrd, achtrx
+
+begin
+	switch (data_type) {
+	case TY_SHORT, TY_USHORT, TY_INT, TY_LONG:
+	    call altmdr (inbuf, inbuf, npix, bscale, bzero)
+	    call achtrl (inbuf, Meml[outbuf], npix)
+	case TY_REAL:
+	    call altmdr (inbuf, inbuf, npix, bscale, bzero)
+	    call amovr (inbuf, Memr[outbuf], npix)
+	case TY_DOUBLE:
+	    call altmdr (inbuf, inbuf, npix, bscale, bzero)
+	    call achtrd (inbuf, Memd[outbuf], npix)
+	case TY_COMPLEX:
+	    call altmdr (inbuf, inbuf, npix, bscale, bzero)
+	    call achtrx (inbuf, Memx[outbuf], npix)
+	default:
+	    call error (10, "RFT_SCALE_LINE: Illegal IRAF image type")
+	}
+end
+
+
+# RFT_DSCALE_PIX -- Procedure to convert an IRAF image line from type double
+# to the requested output data type with optional scaling using the
+# FITS parameters BSCALE and BZERO.
+
+procedure rft_dscale_pix (inbuf, outbuf, npix, bscale, bzero, data_type)
+
+double	inbuf[ARB]		# buffer of FITS integers
+pointer	outbuf			# pointer to output image line
+int	npix			# number of pixels
+double	bscale, bzero		# FITS bscale and bzero
+int	data_type		# IRAF image pixel type
+
+errchk	altmd, achtdl, amovd, achtdr, achtdx
+
+begin
+	switch (data_type) {
+	case TY_SHORT, TY_USHORT, TY_INT, TY_LONG:
+	    call altmd (inbuf, inbuf, npix, bscale, bzero)
+	    call achtdl (inbuf, Meml[outbuf], npix)
+	case TY_REAL:
+	    call altmd (inbuf, inbuf, npix, bscale, bzero)
+	    call achtdr (inbuf, Memr[outbuf], npix)
+	case TY_DOUBLE:
+	    call altmd (inbuf, inbuf, npix, bscale, bzero)
+	    call amovd (inbuf, Memd[outbuf], npix)
+	case TY_COMPLEX:
+	    call altmd (inbuf, inbuf, npix, bscale, bzero)
+	    call achtdx (inbuf, Memx[outbuf], npix)
+	default:
+	    call error (10, "RFT_SCALE_LINE: Illegal IRAF image type")
+	}
+end
+
+
+
+# RFT_LSCALE_PIX -- Procedure to convert an IRAF image line from type long
+# to the requested output data type with optional scaling using the
+# FITS parameters BSCALE and BZERO.
+
+procedure rft_lscale_pix (inbuf, outbuf, npix, bscale, bzero, data_type)
+
+long	inbuf[ARB]		# buffer of FITS integers
+pointer	outbuf			# pointer to output image line
+int	npix			# number of pixels
+double	bscale, bzero		# FITS bscale and bzero
+int	data_type		# IRAF image pixel type
+
+errchk	achtll, achtlr, achtld, achtlx
+errchk	altml, altmr, altmd, altmx
+
+begin
+	switch (data_type) {
+	case TY_SHORT, TY_USHORT, TY_INT, TY_LONG:
+	    call achtll (inbuf, Meml[outbuf], npix)
+	    call altml (Meml[outbuf], Meml[outbuf], npix, bscale, bzero)
+	case TY_REAL:
+	    call altmlr (inbuf, Memr[outbuf], npix, bscale, bzero)
+	case TY_DOUBLE:
+	    call achtld (inbuf, Memd[outbuf], npix)
+	    call altmd (Memd[outbuf], Memd[outbuf], npix, bscale, bzero)
+	case TY_COMPLEX:
+	    call achtlx (inbuf, Memx[outbuf], npix)
+	    call altmx (Memx[outbuf], Memx[outbuf], npix, real (bscale),
+	        real (bzero))
+	default:
+	    call error (10, "RFT_SCALE_LINE: Illegal IRAF image type")
+	}
+end
+
+
+# RFT_RCHANGE_PIX -- Procedure to change a line of real numbers to the
+# IRAF image type.
+
+procedure rft_rchange_pix (inbuf, outbuf, npix, data_type)
+
+real	inbuf[ARB]		# array of FITS integers
+pointer	outbuf			# pointer to IRAF image line
+int	npix			# number of pixels
+int	data_type		# IRAF pixel type
+
+errchk	achtrl, amovr, achtrd, achtrx
+
+begin
+	switch (data_type) {
+	case TY_SHORT, TY_INT, TY_USHORT, TY_LONG:
+	    call achtrl (inbuf, Meml[outbuf], npix)
+	case TY_REAL:
+	    call amovr (inbuf, Memr[outbuf], npix)
+	case TY_DOUBLE:
+	    call achtrd (inbuf, Memd[outbuf], npix)
+	case TY_COMPLEX:
+	    call achtrx (inbuf, Memx[outbuf], npix)
+	default:
+	    call error (10, "RFT_RCHANGE_LINE: Illegal IRAF image type")
+	}
+end
+
+
+# RFT_DCHANGE_PIX -- Procedure to change a line of double precision numbers
+# to the IRAF image type.
+
+procedure rft_dchange_pix (inbuf, outbuf, npix, data_type)
+
+double  inbuf[ARB]		# array of FITS integers
+pointer	outbuf			# pointer to IRAF image line
+int	npix			# number of pixels
+int	data_type		# IRAF pixel type
+
+errchk	achtdl, achtdr, amovd, achtdx
+
+begin
+	switch (data_type) {
+	case TY_SHORT, TY_INT, TY_USHORT, TY_LONG:
+	    call achtdl (inbuf, Meml[outbuf], npix)
+	case TY_REAL:
+	    call achtdr (inbuf, Memr[outbuf], npix)
+	case TY_DOUBLE:
+	    call amovd (inbuf, Memd[outbuf], npix)
+	case TY_COMPLEX:
+	    call achtdx (inbuf, Memx[outbuf], npix)
+	default:
+	    call error (10, "RFT_DCHANGE_LINE: Illegal IRAF image type")
+	}
+end
+
+
+
+# RFT_LCHANGE_PIX -- Procedure to change a line of long integers to the
+# IRAF image type.
+
+procedure rft_lchange_pix (inbuf, outbuf, npix, data_type)
+
+long	inbuf[ARB]		# array of FITS integers
+pointer	outbuf			# pointer to IRAF image line
+int	npix			# number of pixels
+int	data_type		# IRAF pixel type
+
+begin
+	switch (data_type) {
+	case TY_SHORT, TY_USHORT, TY_INT, TY_LONG:
+	    call achtll (inbuf, Meml[outbuf], npix)
+	case TY_REAL:
+	    call achtlr (inbuf, Memr[outbuf], npix)
+	case TY_DOUBLE:
+	    call achtld (inbuf, Memd[outbuf], npix)
+	case TY_COMPLEX:
+	    call achtlx (inbuf, Memx[outbuf], npix)
+	default:
+	    call error (10, "RFT_CHANGE_LINE: Illegal IRAF image type")
+	}
+end
+
+
+# RFT_MAP_BLANKS -- Map the blank pixels.  Currently only the number of blank
+# pixels is determined without an further mapping.
+
+procedure rft_map_blanks (a, buf, npts, pixtype, blank_value, blank, nbadpix)
+
+long	a[ARB]		# integer input buffer
+pointer	buf		# pointer to output image buffer
+int	npts		# number of points
+int	pixtype		# image data type
+long	blank_value	# FITS blank value
+real	blank		# user blank value
+long	nbadpix		# number of bad pixels
+
+int	i
+
+begin
+	# Do blank mapping here
+	switch (pixtype) {
+	case TY_SHORT, TY_INT, TY_USHORT, TY_LONG:
+	    do i = 1, npts {
+	        if (a[i] == blank_value) {
+		    nbadpix = nbadpix + 1
+		    Meml[buf+i-1] = blank
+		}
+	    }
+	case TY_REAL:
+	    do i = 1, npts {
+	        if (a[i] == blank_value) {
+		    nbadpix = nbadpix + 1
+		    Memr[buf+i-1] = blank
+		}
+	    }
+	case TY_DOUBLE:
+	    do i = 1, npts {
+	        if (a[i] == blank_value) {
+		    nbadpix = nbadpix + 1
+		    Memd[buf+i-1] = blank
+		}
+	    }
+	case TY_COMPLEX:
+	    do i = 1, npts {
+	        if (a[i] == blank_value) {
+		    nbadpix = nbadpix + 1
+		    Memx[buf+i-1] = blank
+		}
+	    }
+	}
+end
+
+
+# RFT_PIX_LIMITS -- Procedure to determine to maxmimum and minimum values in a
+# line. Note that double precision is somewhat of a special case because
+# MAX_DOUBLE is a lot less than the maximum permitted ieee numbers for iraf.
+
+procedure rft_pix_limits (buf, npix, pixtype, linemin, linemax)
+
+pointer	buf				# pointer to IRAF image line
+int	npix				# number of pixels
+int	pixtype				# output data type
+real	linemax, linemin		# min and max pixel values
+
+long	lmax, lmin
+real	rmax, rmin
+double	dmax, dmin
+complex	xmax, xmin
+
+begin
+	switch (pixtype) {
+	case TY_SHORT, TY_INT, TY_USHORT, TY_LONG:
+	    call aliml (Meml[buf], npix, lmin, lmax)
+	    linemax = lmax
+	    linemin = lmin
+	case TY_REAL:
+	    call alimr (Memr[buf], npix, rmin, rmax)
+	    linemax = rmax
+	    linemin = rmin
+	case TY_DOUBLE:
+	    call alimd (Memd[buf], npix, dmin, dmax)
+	    if (dmax > MAX_REAL)
+	        linemax = INDEFR
+	    else
+		linemax = dmax
+	    if (dmin < -MAX_REAL)
+		linemin = INDEFR
+	    else
+	        linemin = dmin
+	case TY_COMPLEX:
+	    call alimx (Memx[buf], npix, xmin, xmax)
+	    linemax = xmax
+	    linemin = xmin
+	default:
+	    call error (30, "RFT_PIX_LIMITS: Unknown IRAF type")
+	}
+end
+
+
+# ALTMDR -- procedure to scale a long vector into a real vector using
+# double precision constants to preserve accuracy
+
+procedure altmlr (a, b, npix, bscale, bzero)
+
+long	a[ARB]		# input array
+real	b[ARB]		# output array
+int	npix		# number of pixels
+double	bscale, bzero	# scaling parameters
+
+int	i
+
+begin
+	do i = 1, npix
+	    b[i] = a[i] * bscale + bzero
+end
+
+
+# ALTMDR -- procedure to scale a real vector with double precision constants.
+
+procedure altmdr (a, b, npix, bscale, bzero)
+
+real	a[ARB]		# input array
+real	b[ARB]		# output array
+int	npix		# number of pixels
+double	bscale, bzero	# scaling parameters
+
+int	i
+
+begin
+	do i = 1, npix
+	    b[i] = a[i] * bscale + bzero
+end
diff --git a/pkg/dataio/fits/fits_rpixels.x b/pkg/dataio/fits/fits_rpixels.x
new file mode 100644
index 00000000..6183120d
--- /dev/null
+++ b/pkg/dataio/fits/fits_rpixels.x
@@ -0,0 +1,154 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <fset.h>
+include <mii.h>
+include <mach.h>
+
+# RFT_INIT_READ_PIXELS and READ_PIXELS -- Read pixel data with record buffering
+# and data type conversion.  The input data must meet the MII standard
+# except for possibly in the case of integers having the least significant
+# byte first.
+#
+# Read data in records of len_record and convert to the specified IRAF
+# data type.  Successive calls of rft_read_pixels returns the next npix pixels.
+# Read_pixels returns EOF or the number of pixels converted.
+# Init_read_pixels must be called before read_pixels.
+#
+# Error conditions are:
+# 1. A short input record
+# 2. Error in converting the pixels by miiup.
+#
+# This routine is based on the MII unpack routine which is machine dependent.
+# The bitpix must correspond to an MII type.  If the lsbf (least significant
+# byte first) flag is YES then the pixels do not satisfy the MII standard.
+# In this case the bytes are first swapped into most significant byte first
+# before the MII unpack routine is called.
+
+int procedure rft_init_read_pixels (npix_record, bitpix, lsbf, spp_type)
+
+int	npix_record	# number of pixels per input record
+int	bitpix		# bits per pixel (must correspond to an MII type)
+int	lsbf		# byte swap?
+int	spp_type	# SPP data type to be returned
+
+# entry rft_read_pixels (fd, buffer, npix)
+
+int	rft_read_pixels
+int	fd		# input file descriptor
+char	buffer[1]	# output buffer
+int	npix		# number of pixels to read
+
+int	swap
+int	ty_mii, ty_spp, npix_rec, nch_rec, sz_rec, nchars, len_mii, recptr
+int	bufsize, i, n, ip, op
+pointer	mii, spp
+
+int	rft_getbuf(), sizeof(), miilen()
+errchk	miilen, mfree, malloc, rft_getbuf, miiupk
+data	mii/NULL/, spp/NULL/
+
+begin
+	ty_mii = bitpix
+	ty_spp = spp_type
+	swap = lsbf
+	npix_rec = npix_record
+	nch_rec = npix_rec * sizeof (ty_spp)
+
+	len_mii = miilen (npix_rec, ty_mii)
+	sz_rec = len_mii * SZ_INT32
+
+	if (mii != NULL)
+	    call mfree (mii, TY_INT)
+	call malloc (mii, len_mii, TY_INT)
+
+	if (spp != NULL)
+	    call mfree (spp, TY_CHAR)
+	call malloc (spp, nch_rec, TY_CHAR)
+
+	ip = nch_rec
+	return (OK)
+
+entry	rft_read_pixels (fd, buffer, npix, recptr, bufsize)
+
+	nchars = npix * sizeof (ty_spp)
+	op = 0
+
+	repeat {
+
+	    # If data is exhausted read the next record
+	    if (ip == nch_rec) {
+
+		i = rft_getbuf (fd, Memi[mii], sz_rec, bufsize, recptr)
+		if (i == EOF)
+		    return (EOF)
+
+		if (swap == YES)
+		    switch (ty_mii) {
+		    case MII_SHORT:
+			call bswap2 (Memi[mii], 1, Memi[mii], 1,
+				sz_rec * SZB_CHAR)
+		    case MII_LONG:
+			call bswap4 (Memi[mii], 1, Memi[mii], 1,
+				sz_rec * SZB_CHAR)
+		    }
+
+		call miiupk (Memi[mii], Memc[spp], npix_rec, ty_mii, ty_spp)
+
+		ip = 0
+		#recptr = recptr + 1
+	    }
+
+	    n = min (nch_rec - ip, nchars - op)
+	    call amovc (Memc[spp+ip], buffer[1+op], n)
+	    ip = ip + n
+	    op = op + n
+
+	} until (op == nchars)
+
+	return (npix)
+end
+
+
+# RFT_GETBUF -- Procedure to get the buffer.
+
+int procedure rft_getbuf (fd, buf, sz_rec, bufsize, recptr)
+
+int	fd		# file descriptor
+char	buf[ARB]	# buffer to be filled
+int	sz_rec		# size in chars of record to be read
+int	bufsize		# buffer size in records
+int	recptr		# last successful FITS record read
+
+int	i, nchars
+int	read(), fstati()
+errchk	read
+
+begin
+	nchars = 0
+	repeat {
+	    iferr {
+	        i = read (fd, buf[nchars+1], sz_rec - nchars)
+	    } then {
+	        call printf ("Error reading FITS record %d\n")
+	        if (mod (recptr + 1, bufsize) == 0)
+		    call pargi ((recptr + 1) / bufsize)
+	        else
+		    call pargi ((recptr + 1) / bufsize + 1)
+	        call fseti (fd, F_VALIDATE, fstati (fd, F_SZBBLK) / SZB_CHAR)
+	        i = read (fd, buf[nchars+1], sz_rec - nchars)
+	    }
+
+	    if (i == EOF)
+		break
+	    else
+	        nchars = nchars + i
+
+	} until (nchars >= sz_rec)
+
+	if ((i == EOF) && (nchars == 0))
+	    return (EOF)
+	else {
+	    recptr = recptr + 1
+	    return (nchars)
+	}
+end
diff --git a/pkg/dataio/fits/fits_wheader.x b/pkg/dataio/fits/fits_wheader.x
new file mode 100644
index 00000000..ec06f968
--- /dev/null
+++ b/pkg/dataio/fits/fits_wheader.x
@@ -0,0 +1,469 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include <mach.h>
+include <fset.h>
+include "wfits.h"
+
+# WFT_WRITE_HEADER -- Write the FITS headers. The FITS header
+# parameters are encoded one by one until the FITS END keyword is detected.
+# If the long_header switch is set the full FITS header is printed on the
+# standard output. If the short header parameter is specified only the image
+# title and dimensions are printed.
+
+procedure wft_write_header (im, fits, fits_fd)
+
+pointer	im		# pointer to the IRAF image
+pointer	fits		# pointer to the FITS structure
+int	fits_fd		# the FITS file descriptor
+
+char	card[LEN_CARD+1], trim_card[LEN_CARD+1]
+int	nrecords, recntr, cardptr, cardcnt, stat, cards_per_rec, i
+int	wft_card_encode(), wft_set_bitpix(), sizeof(), strncmp()
+int	wft_init_card_encode(), fstati()
+
+errchk	wft_set_bitpix, wft_get_iraf_typestring, wft_set_scale, wft_set_blank
+errchk	wft_fits_set_scale, wft_init_card_encode, wft_card_encode
+errchk	wft_init_write_pixels, wft_write_pixels, wft_write_last_record
+
+include "wfits.com"
+
+begin
+	# SET the data and FITS bits per pixel.
+
+	DATA_BITPIX(fits) = sizeof (PIXTYPE(im)) * SZB_CHAR * NBITS_BYTE
+	FITS_BITPIX(fits) = wft_set_bitpix (bitpix, PIXTYPE(im),
+	    DATA_BITPIX(fits))
+
+	# Calculate the FITS bscale and bzero parameters. Notice for the
+	# time being that scaling is turned off if IEEE floating point
+	# output is selected. May decide to change this later after
+	# checking the specifications.
+
+	if (FITS_BITPIX(fits) < 0) {
+
+	    IRAFMIN(fits) = IM_MIN(im)
+	    IRAFMAX(fits) = IM_MAX(im)
+	    SCALE(fits) = NO
+	    BZERO(fits) = 0.0d0
+	    BSCALE(fits) = 1.0d0
+
+	} else if (autoscale == YES) {
+
+	    call wft_get_tape_limits (FITS_BITPIX(fits), TAPEMIN(fits),
+	        TAPEMAX(fits))
+	    call wft_data_limits (im, IRAFMIN(fits), IRAFMAX(fits))
+	    call wft_fits_set_scale (im, DATA_BITPIX(fits), FITS_BITPIX(fits),
+	        IRAFMIN(fits), IRAFMAX(fits), TAPEMIN(fits), TAPEMAX(fits),
+		SCALE(fits), BSCALE(fits), BZERO(fits))
+
+	} else {
+
+	    IRAFMIN(fits) = IM_MIN(im)
+	    IRAFMAX(fits) = IM_MAX(im)
+	    SCALE(fits) = scale
+	    BZERO(fits) = bzero
+	    BSCALE(fits) = bscale
+	}
+
+	# If blanks in the image set the blank parameter. Currently information
+	# on blanks is not written out so this is effectively a null operation
+	# in IRAF.
+
+	if (NBPIX(im) > 0)
+	    call wft_set_blank (FITS_BITPIX(fits), BLANK(fits),
+	        BLANK_STRING(fits))
+
+	# Set the IRAF datatype parameter.
+	call wft_get_iraf_typestring (PIXTYPE(im), TYPE_STRING(fits))
+
+	# Initialize the card counters. These counters are used only for
+	# information printed to the standard output.
+
+	recntr = 1
+	cardptr = 1
+	cardcnt = 1
+	cards_per_rec = len_record / LEN_CARD
+
+	# Get set up to write the FITS header. Initialize for an ASCII write.
+	stat = wft_init_card_encode (im, fits)
+	if (make_image == YES)
+	    call wft_init_wrt_pixels (len_record, TY_CHAR, FITS_BYTE, blkfac)
+
+	# Print short header.
+	if (short_header == YES && long_header == NO) {
+
+	    call printf ("%-20.20s  ")
+		call pargstr (OBJECT(im))
+	    do i = 1, NAXIS(im) {
+		if (i == 1) {
+		    call printf ("Size = %d")
+			call pargl (NAXISN(im,i))
+		} else {
+		    call printf (" x %d")
+			call pargl (NAXISN(im,i))
+		}
+	    }
+	    call printf ("\n")
+
+	    call strlwr (TYPE_STRING(fits))
+	    call printf ("\tpixtype=%s bitpix=%d")
+		call pargstr (TYPE_STRING(fits))
+		call pargi (FITS_BITPIX(fits))
+
+	    if (fstati (fits_fd, F_BLKSIZE) == 0) {
+		call printf (" blkfac=%d")
+		    call pargi (blkfac)
+	    } else
+		call printf (" blkfac=fixed")
+
+	    if (SCALE(fits) == YES) {
+	        call printf (" bscale=%.7g bzero=%.7g\n")
+		    call pargd (BSCALE(fits))
+		    call pargd (BZERO(fits))
+	    } else
+		call printf (" scaling=none\n")
+	    call strupr (TYPE_STRING(fits))
+	}
+
+	# Write the cards to the FITS header.
+	repeat {
+
+	    # Encode the card.
+	    stat = wft_card_encode (im, fits, card)
+	    if (stat == NO)
+		next
+
+	    # Write the card to the output file if make_image is yes.
+	    if (make_image == YES)
+	        call wft_write_pixels (fits_fd, card, LEN_CARD)
+
+	    # Trim the card and write is to the standard output if 
+	    # long_header is yes.
+
+	    if (long_header == YES) {
+		call wft_trimstr (card, trim_card, LEN_CARD)
+	        call printf ("%2d/%2d:--  %s\n")
+		    call pargi (recntr)
+		    call pargi (cardptr)
+		    call pargstr (trim_card)
+	    }
+
+	    if (mod (cardcnt, cards_per_rec) == 0) {
+	        recntr = recntr + 1
+	        cardptr = 1
+	    } else
+		cardptr = cardptr + 1
+	    cardcnt = cardcnt + 1
+
+	} until (strncmp (card, "END     ", LEN_KEYWORD) == 0)
+
+	# Issue warning about possible precision loss. Comment this out
+	# for time being, since the short header was modified.
+	#if (SCALE(fits) == YES && bitpix != ERR) {
+	    #call printf (
+	    #"\tDefault bitpix overridden: maximum precision loss ~%.7g\n")
+		#call pargd (BSCALE(fits))
+	#}
+
+	# Write the last header records.
+	if (make_image == YES) {
+	    call wft_write_last_record (fits_fd, nrecords)
+	    if (short_header == YES || long_header == YES) {
+	        call printf ("\t%d Header  ")
+		    call pargi (nrecords)
+	    }
+	}
+end
+
+
+# WFT_SET_BITPIX -- This procedure sets the FITS bitpix for each image based on
+# either the user given value or the precision of the IRAF data. Notice that
+# the user must explicitly set the bitpix parameter to -16 or -32 to select
+# the IEEE output format.
+
+int procedure wft_set_bitpix (bitpix, datatype, data_bitpix)
+
+int	bitpix		# the user set bits per pixel, ERR or legal bitpix
+int	datatype	# the IRAF image data type
+int	data_bitpix	# the bits per pixel in the data
+
+begin
+	if (bitpix == ERR) {
+	    switch (datatype) {
+	    case TY_SHORT, TY_INT, TY_USHORT, TY_LONG:
+	        if (data_bitpix <= FITS_BYTE)
+		    return (FITS_BYTE)
+	        else if (data_bitpix <= FITS_SHORT) {
+		    #if (datatype == TY_USHORT)
+			#return (FITS_LONG)
+		    #else
+		    return (FITS_SHORT)
+	        } else
+		    return (FITS_LONG)
+	    case TY_REAL, TY_COMPLEX:
+		return (FITS_REAL)
+	    case TY_DOUBLE:
+		 return (FITS_DOUBLE)
+	    default:
+		call error (2, "SET_BITPIX: Unknown IRAF data type.")
+	    }
+	} else
+	    return (bitpix)
+end
+
+
+# WFT_GET_IRAF_TYPESTRING -- Procedure to set the iraf datatype keyword.
+
+procedure wft_get_iraf_typestring (datatype, type_str)
+
+int	datatype	# the IRAF data type
+char	type_str[ARB]	# the output IRAF type string
+
+begin
+	switch (datatype) {
+	case TY_SHORT:
+	    call strcpy ("SHORT", type_str, LEN_STRING) 
+	case TY_USHORT:
+	    call strcpy ("USHORT", type_str, LEN_STRING) 
+	case TY_INT: 
+	    call strcpy ("INTEGER", type_str, LEN_STRING) 
+	case TY_LONG: 
+	    call strcpy ("LONG", type_str, LEN_STRING) 
+	case TY_REAL:
+	    call strcpy ("REAL", type_str, LEN_STRING)
+	case TY_DOUBLE:
+	    call strcpy ("DOUBLE", type_str, LEN_STRING)
+	case TY_COMPLEX:
+	    call strcpy ("COMPLEX", type_str, LEN_STRING)
+	default:
+	    call error (3, "IRAF_TYPE: Unknown IRAF image type.")
+	}
+end
+
+
+# WFT_FITS_SET_SCALE -- Procedure to set the FITS scaling parameters if
+# autoscaling is enabled.
+
+procedure wft_fits_set_scale (im, data_bitpix, fits_bitpix, irafmin, irafmax,
+	tapemin, tapemax, scale, bscale, bzero )
+
+pointer	im		# pointer to IRAF image
+int	data_bitpix	# bits per pixel of data
+int	fits_bitpix	# fits bits per pixel
+real	irafmin		# minimum picture value
+real	irafmax		# maximum picture value
+double	tapemin		# minimum tape value
+double	tapemax		# maximum tape value
+int	scale		# scale data ?
+double	bscale		# FITS bscale
+double	bzero		# FITS bzero
+
+errchk	wft_set_scale
+
+begin
+	switch (PIXTYPE(im)) {
+	case TY_SHORT, TY_INT, TY_LONG:
+	    if (data_bitpix > fits_bitpix) {
+		scale = YES
+		call wft_set_scale (fits_bitpix, irafmin, irafmax, tapemin, 
+		    tapemax, bscale, bzero)
+	    } else {
+		scale = NO
+		bscale = 1.0d0
+		bzero = 0.0d0
+	    }
+	case TY_USHORT:
+	    if (data_bitpix > fits_bitpix) {
+		scale = YES
+		call wft_set_scale (fits_bitpix, irafmin, irafmax, tapemin, 
+		    tapemax, bscale, bzero)
+	    } else if (data_bitpix == fits_bitpix) {
+		scale = YES
+		bscale = 1.0d0
+		bzero = 3.2768d4
+	    } else {
+		scale = NO
+		bscale = 1.0d0
+		bzero = 0.0d0
+	    }
+	case TY_REAL, TY_DOUBLE, TY_COMPLEX:
+	    scale = YES
+	    call wft_set_scale (fits_bitpix, irafmin, irafmax, tapemin, tapemax,
+	        bscale, bzero)
+	default:
+	    call error (1, "WRT_HEADER: Unknown IRAF image type.")
+	}
+
+end
+
+
+# WFT_SET_SCALE -- This procedure calculates bscale and bzero for each frame
+# from the FITS bitpix and the maximum and minimum data values.
+
+procedure wft_set_scale (fits_bitpix, datamin, datamax, mintape, maxtape,
+	bscale, bzero)
+
+int	fits_bitpix		# the FITS integer bits per pixels
+real	datamax, datamin	# the IRAF image data minimum and maximum
+double	mintape, maxtape	# min and max FITS tape values
+double	bscale, bzero		# the calculated bscale and bzero values
+
+double	maxdata, mindata, num, denom
+bool	rft_equald()
+
+begin
+	# Calculate the maximum and minimum values in the data.
+	maxdata = datamax #+ abs ((datamax / (10.0 ** (NDIGITS_RP - 1))))
+	mindata = datamin #- abs ((datamin / (10.0 ** (NDIGITS_RP - 1))))
+	num = maxdata - mindata
+	denom =  (maxtape - mintape) * PREC_RATIO
+
+	# Check for constant image case.
+	#mindata = datamin
+	#maxdata = datamax
+	if (rft_equald (num, 0.0d0)) {
+	    bscale = 1.0d0
+	    bzero = maxdata
+	} else {
+	    bscale = num / denom
+	    #bzero = (maxtape / denom) * mindata - (mintape / denom) * maxdata
+	    bzero = (maxdata + mindata) / 2.0d0
+	}
+end
+
+
+# WFT_GET_TAPE_LIMITS -- Procedure for calculating the maximum and minimum FITS
+# integer values from the FITS bitpix.
+
+procedure wft_get_tape_limits (fits_bitpix, mintape, maxtape)
+
+int	fits_bitpix		# the bits per pixel of a FITS integer
+double	maxtape, mintape	# the maximun and minimum FITS tape integers
+
+begin
+	switch (fits_bitpix) {
+	case FITS_BYTE:
+	    maxtape = BYTE_MAX
+	    mintape = BYTE_MIN
+	case FITS_SHORT:
+	    maxtape = SHORT_MAX
+	    mintape = SHORT_MIN
+	case FITS_LONG:
+	    maxtape = LONG_MAX
+	    mintape = LONG_MIN
+	default:
+	    call error (4, "TAPE_LIMITS: Unknown FITS type.")
+	}
+end
+
+
+# WFT_SET_BLANK -- Determine the FITS integer value for a blank pixel from the
+# FITS bitpix. Notice that these are null ops for IEEE floating point format.
+
+procedure wft_set_blank (fits_bitpix, blank, blank_str)
+
+int	fits_bitpix	# the requested FITS bits per pixel
+long	blank		# the FITS integer value of a blank pixel
+char	blank_str[ARB]	# the encoded FITS integer value of a blank pixel
+
+begin
+	switch (fits_bitpix) {
+	case FITS_BYTE:
+	    blank = long (BYTE_BLANK)
+	    call strcpy ("0", blank_str, LEN_BLANK)
+	case FITS_SHORT:
+	    blank = long (SHORT_BLANK)
+	    call strcpy ("-32768", blank_str, LEN_BLANK)
+	case FITS_LONG:
+	    blank = long (LONG_BLANK)
+	    call strcpy ("-2147483648", blank_str, LEN_BLANK)
+	case FITS_REAL:
+	    blank = INDEFL
+	    call strcpy ("", blank_str, LEN_BLANK)
+	case FITS_DOUBLE:
+	    blank = INDEFL
+	    call strcpy ("", blank_str, LEN_BLANK)
+	default:
+	    call error (5, "SET_BLANK: Unknown FITS type.")
+	}
+end
+
+
+# WFT_INIT_CARD_ENCODE -- This procedure initializes the card encoding
+# procedure.  The cards counters are initialized and the number of history cards
+# calculated.
+
+int procedure wft_init_card_encode (im, fits)
+
+# both entry points
+pointer	im	# pointer to the IRAF image
+pointer	fits	# pointer to the WFITS structure
+
+# entry wft_card_encode
+int	wft_card_encode		# entry point
+char	card[LEN_CARD+1]	# string containing the card image
+
+int	cardno, axisno, optiono, hist_ptr, unknown_ptr
+int	nstandard, noptions, stat
+int	wft_standard_card(), wft_option_card(), wft_last_card()
+int	wft_history_card(), wft_unknown_card()
+errchk	wft_standard_card, wft_option_card, wft_history_card
+errchk	wft_unknown_card, wft_last_card
+
+begin
+	# Initialize the card pointers.
+	cardno = 1
+	axisno = 1
+	optiono = 1
+	unknown_ptr = 1
+	hist_ptr = 1
+
+	# Initilaize the card counters.
+	nstandard = 3 + NAXIS(im)
+	noptions = NOPTIONS + nstandard
+
+	return (YES)
+
+
+# WFT_CARD_ENCODE -- Procedure to encode the FITS header parameters into
+# FITS card images.
+
+entry	wft_card_encode (im, fits, card)
+
+	# Fetch the appropriate FITS header card image.
+	if (cardno <= nstandard) {
+	    stat = wft_standard_card (cardno, im, fits, axisno, card)
+	} else if (cardno <= noptions) {
+	    stat = wft_option_card (im, fits, optiono, card)
+	} else if (wft_unknown_card (im, unknown_ptr, card) == YES) {
+	    stat = YES
+	} else if (wft_history_card (im, hist_ptr, card) == YES) {
+	    stat = YES
+	} else {
+	    stat = wft_last_card (card)
+	}
+
+	cardno = cardno + 1
+
+	return (stat)
+end
+
+
+# WFT_TRIMSTR -- Procedure to trim trailing blanks from a fixed size string.
+
+procedure wft_trimstr (instr, outstr, nchars)
+
+char	instr[ARB]	# input string
+char	outstr[ARB]	# output string
+int	nchars		# last character of instr
+
+int	ip
+
+begin
+	call strcpy (instr, outstr, nchars)
+	ip = nchars
+	while (outstr[ip] == ' ')
+	    ip = ip - 1
+	outstr[ip+1] = EOS
+end
diff --git a/pkg/dataio/fits/fits_wimage.x b/pkg/dataio/fits/fits_wimage.x
new file mode 100644
index 00000000..7ed00372
--- /dev/null
+++ b/pkg/dataio/fits/fits_wimage.x
@@ -0,0 +1,497 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+include <imhdr.h>
+include "wfits.h"
+
+# WFT_WRITE_IMAGE -- Procedure to convert IRAF image data to FITS format line by
+# line.
+
+procedure wft_write_image (im, fits, fits_fd)
+
+pointer	im			# IRAF image descriptor
+pointer	fits			# FITS data structure
+int	fits_fd			# FITS file descriptor
+
+int	npix, nlines, npix_record, i, stat, nrecords
+long	v[IM_MAXDIM]
+pointer	tempbuf, buf
+
+int	wft_get_image_line()
+errchk	malloc, mfree, wft_get_image_line, wft_lscale_line, wft_long_line
+errchk	wft_init_write_pixels, wft_write_pixels, wft_write_last_record
+errchk	wft_rscale_line, wft_dscale_line
+
+include "wfits.com"
+
+begin
+	if (NAXIS(im) == 0) {
+	    if (short_header == YES || long_header == YES) {
+	        call printf ("0  Data logical (2880 byte) records written\n")
+	    }
+	    return
+	}
+
+	# Initialize.
+	npix = NAXISN(im,1)
+	nlines = 1
+	do i = 2, NAXIS(im)
+	    nlines = nlines * NAXISN(im, i)
+	npix_record = len_record * FITS_BYTE / abs (FITS_BITPIX(fits))
+
+	call amovkl (long(1), v, IM_MAXDIM)
+	switch (FITS_BITPIX(fits)) {
+	case FITS_REAL:
+
+	    # Allocate temporary space.
+	    call malloc (tempbuf, npix, TY_REAL)
+
+	    # Initialize the pixel write.
+	    call wft_init_write_pixels (npix_record, TY_REAL,
+	        FITS_BITPIX(fits), blkfac)
+
+	    # For the time being explicitly turn off ieee NaN mapping.
+	    call ieemapr (NO, NO)
+
+	    # Scale the lines, deal with the blanks via the ieee code which
+	    # is currently turned off, and write the output records.
+
+	    do i = 1, nlines {
+	        iferr (stat =  wft_get_image_line (im, buf, v, PIXTYPE(im))) {
+		    call erract (EA_WARN)
+		    call error (10, "WRT_IMAGE: Error writing IRAF image.")
+	        }
+	        if (stat == EOF )
+		    return
+	        if (stat != npix)
+		    call error (10, "WRT_IMAGE: Error writing IRAF image.")
+	        if (SCALE(fits) == YES)
+	            call wft_rscale_line (buf, Memr[tempbuf], npix,
+		        1. / BSCALE(fits), -BZERO(fits), PIXTYPE(im))
+		else
+		    call wft_real_line (buf, Memr[tempbuf], npix, PIXTYPE(im)) 
+	        call wft_write_pixels (fits_fd, Memr[tempbuf], npix)
+	    }
+
+	    # Free space.
+	    call mfree (tempbuf, TY_REAL)
+
+	case FITS_DOUBLE:
+
+	    # Allocate temporary space.
+	    call malloc (tempbuf, npix, TY_DOUBLE)
+
+	    # Initialize the pixel write.
+	    call wft_init_write_pixels (npix_record, TY_DOUBLE,
+	        FITS_BITPIX(fits), blkfac)
+
+	    # For the time being explicitly turn off ieee NaN mapping.
+	    call ieemapd (NO, NO)
+
+	    # Scale the lines, deal with the blanks via the ieee code which
+	    # is currently turned off, and write the output records.
+
+	    do i = 1, nlines {
+	        iferr (stat =  wft_get_image_line (im, buf, v, PIXTYPE(im))) {
+		    call erract (EA_WARN)
+		    call error (10, "WRT_IMAGE: Error writing IRAF image.")
+	        }
+	        if (stat == EOF )
+		    return
+	        if (stat != npix)
+		    call error (10, "WRT_IMAGE: Error writing IRAF image.")
+	        if (SCALE(fits) == YES)
+	            call wft_dscale_line (buf, Memd[tempbuf], npix,
+		        1. / BSCALE(fits), -BZERO(fits), PIXTYPE(im))
+		else
+		    call wft_double_line (buf, Memd[tempbuf], npix,
+		        PIXTYPE(im)) 
+	        call wft_write_pixels (fits_fd, Memd[tempbuf], npix)
+	    }
+
+	    # Free space.
+	    call mfree (tempbuf, TY_DOUBLE)
+
+	default:
+
+	    # Allocate temporary space.
+	    call malloc (tempbuf, npix, TY_LONG)
+
+	    # Scale the line, deal with the blanks, and write the output
+	    # record. At the moement blanks are not dealt with.
+
+	    call wft_init_write_pixels (npix_record, TY_LONG, FITS_BITPIX(fits),
+	         blkfac)
+	    do i = 1, nlines {
+	        iferr (stat =  wft_get_image_line (im, buf, v, PIXTYPE(im))) {
+		    call erract (EA_WARN)
+		    call error (10, "WRT_IMAGE: Error writing IRAF image.")
+	        }
+	        if (stat == EOF )
+		    return
+	        if (stat != npix)
+		    call error (10, "WRT_IMAGE: Error writing IRAF image.")
+	        if (SCALE(fits) == YES)
+	            call wft_lscale_line (buf, Meml[tempbuf], npix,
+		        1. / BSCALE(fits), -BZERO(fits), PIXTYPE(im))
+	        else
+		    call wft_long_line (buf, Meml[tempbuf], npix, PIXTYPE(im)) 
+	        # call map_blanks (im, Meml[tempbuf], blank)
+	        call wft_write_pixels (fits_fd, Meml[tempbuf], npix)
+	    }
+	    # Free space.
+	    call mfree (tempbuf, TY_LONG)
+	}
+
+	# Write the final record.
+	call wft_write_last_record (fits_fd, nrecords)
+	if (short_header == YES || long_header == YES) {
+	    call printf ("%d  Data logical (2880 byte) records written\n")
+	        call pargi (nrecords)
+	}
+end
+
+
+# WFT_GET_IMAGE_LINE -- Procedure to fetch the next image line.
+
+int procedure wft_get_image_line (im, buf, v, datatype)
+
+pointer	im			# IRAF image descriptor
+pointer	buf			# pointer to image line
+long	v[ARB]			# imio dimension descriptor
+int	datatype		# IRAF image data type
+
+int	npix
+int	imgnll(), imgnlr(), imgnld(), imgnlx()
+errchk	imgnll, imgnlr, imgnld, imgnlx
+
+begin
+	switch (datatype) {
+	case TY_SHORT, TY_INT, TY_USHORT, TY_LONG:
+	    npix = imgnll (im, buf, v)
+	case TY_REAL:
+	    npix = imgnlr (im, buf, v)
+	case TY_DOUBLE:
+	    npix = imgnld (im, buf, v)
+	case TY_COMPLEX:
+	    npix = imgnlx (im, buf, v)
+	default:
+	    call error (11, "GET_IMAGE_LINE: Unknown IRAF image type.")
+	}
+
+	return (npix)
+end
+
+
+# WFT_RSCALE_LINE -- This procedure converts the IRAF data to type real
+# and scales by the FITS parameters bscale and bzero.
+
+procedure wft_rscale_line (buf, outbuffer, npix, bscale, bzero, datatype)
+
+pointer	buf			# pointer to IRAF image line
+real	outbuffer[ARB]		# FITS integer buffer
+int	npix			# number of pixels
+double	bscale, bzero		# FITS bscale and bzero parameters
+int	datatype		# data type of image
+
+errchk	achtlr, altadr, amovr, achtdr, acthxr
+
+begin
+	switch (datatype) {
+	case TY_SHORT, TY_INT, TY_LONG, TY_USHORT:
+	    call achtlr (Meml[buf], outbuffer, npix)
+	    call altadr (outbuffer, outbuffer, npix, bzero, bscale)
+	case TY_REAL:
+	    call amovr (Memr[buf], outbuffer, npix)
+	    call altadr (outbuffer, outbuffer, npix, bzero, bscale)
+	case TY_DOUBLE:
+	    call achtdr (Memd[buf], outbuffer, npix)
+	    call altadr (outbuffer, outbuffer, npix, bzero, bscale)
+	case TY_COMPLEX:
+	    call achtxr (Memx[buf], outbuffer, npix)
+	    call altadr (outbuffer, outbuffer, npix, bzero, bscale)
+	default:
+	    call error (12, "WFT_RSCALE_LINE: Unknown IRAF image type.")
+	}
+end
+
+
+# WFT_DSCALE_LINE -- This procedure converts the IRAF data to type double with
+# after scaling by the FITS parameters bscale and bzero.
+
+procedure wft_dscale_line (buf, outbuffer, npix, bscale, bzero, datatype)
+
+pointer	buf			# pointer to IRAF image line
+double	outbuffer[ARB]		# FITS integer buffer
+int	npix			# number of pixels
+double	bscale, bzero		# FITS bscale and bzero parameters
+int	datatype		# data type of image
+
+errchk	achtld, altad, amovd, achtrd, achtxd
+
+begin
+	switch (datatype) {
+	case TY_SHORT, TY_INT, TY_LONG, TY_USHORT:
+	    call achtld (Meml[buf], outbuffer, npix)
+	    call altad (outbuffer, outbuffer, npix, bzero, bscale)
+	case TY_REAL:
+	    call achtrd (Memr[buf], outbuffer, npix)
+	    call altad (outbuffer, outbuffer, npix, bzero, bscale)
+	case TY_DOUBLE:
+	    call amovd (Memd[buf], outbuffer, npix)
+	    call altad (outbuffer, outbuffer, npix, bzero, bscale)
+	case TY_COMPLEX:
+	    call achtxd (Memx[buf], outbuffer, npix)
+	    call altad (outbuffer, outbuffer, npix, bzero, bscale)
+	default:
+	    call error (12, "WFT_DSCALE_LINE: Unknown IRAF image type.")
+	}
+end
+
+
+# WFT_LSCALE_LINE -- This procedure converts the IRAF data to type long with
+# after scaling by the FITS parameters bscale and bzero.
+
+procedure wft_lscale_line (buf, outbuffer, npix, bscale, bzero, datatype)
+
+pointer	buf			# pointer to IRAF image line
+long	outbuffer[ARB]		# FITS integer buffer
+int	npix			# number of pixels
+double	bscale, bzero		# FITS bscale and bzero parameters
+int	datatype		# data type of image
+
+errchk  altall, altarl, altadl, altaxl
+
+begin
+	switch (datatype) {
+	case TY_SHORT, TY_INT, TY_LONG, TY_USHORT:
+	    call altall (Meml[buf], outbuffer, npix, bzero, bscale)
+	case TY_REAL:
+	    call altarl (Memr[buf], outbuffer, npix, bzero, bscale)
+	case TY_DOUBLE:
+	    call altadl (Memd[buf], outbuffer, npix, bzero, bscale)
+	case TY_COMPLEX:
+	    call altaxl (Memx[buf], outbuffer, npix, bzero, bscale)
+	default:
+	    call error (12, "WFT_LSCALE_LINE: Unknown IRAF image type.")
+	}
+end
+
+
+# WFT_REAL_LINE -- This procedure converts the IRAF image line to type long with
+# no scaling.
+
+procedure wft_real_line (buf, outbuffer, npix, datatype)
+
+pointer	buf			# pointer to IRAF image line
+real	outbuffer[ARB]		# buffer of FITS integers
+int	npix			# number of pixels
+int	datatype		# IRAF image datatype
+
+errchk	achtlr, achtdr, amovr, achtxr
+
+begin
+	switch (datatype) {
+	case TY_SHORT, TY_INT, TY_LONG, TY_USHORT:
+	    call achtlr (Meml[buf], outbuffer, npix)
+	case TY_REAL:
+	    call amovr (Memr[buf], outbuffer, npix)
+	case TY_DOUBLE:
+	    call achtdr (Memd[buf], outbuffer, npix)
+	case TY_COMPLEX:
+	    call achtxr (Memx[buf], outbuffer, npix)
+	default:
+	    call error (13, "WFT_REAL_LINE: Unknown IRAF data type.")
+	}
+end
+
+
+# WFT_DOUBLE_LINE -- This procedure converts the IRAF image line to type long
+# with no scaling.
+
+procedure wft_double_line (buf, outbuffer, npix, datatype)
+
+pointer	buf			# pointer to IRAF image line
+double  outbuffer[ARB]		# buffer of FITS integers
+int	npix			# number of pixels
+int	datatype		# IRAF image datatype
+
+errchk	achtld, achtrd, amovd, achtxd
+
+begin
+	switch (datatype) {
+	case TY_SHORT, TY_INT, TY_LONG, TY_USHORT:
+	    call achtld (Meml[buf], outbuffer, npix)
+	case TY_REAL:
+	    call achtrd (Memr[buf], outbuffer, npix)
+	case TY_DOUBLE:
+	    call amovd (Memd[buf], outbuffer, npix)
+	case TY_COMPLEX:
+	    call achtxd (Memx[buf], outbuffer, npix)
+	default:
+	    call error (13, "WFT_DOUBLE_LINE: Unknown IRAF data type.")
+	}
+end
+
+
+# WFT_LONG_LINE -- This procedure converts the IRAF image line to type long with
+# no scaling.
+
+procedure wft_long_line (buf, outbuffer, npix, datatype)
+
+pointer	buf			# pointer to IRAF image line
+long	outbuffer[ARB]		# buffer of FITS integers
+int	npix			# number of pixels
+int	datatype		# IRAF image datatype
+
+errchk	amovl, achtrl, achtdl, achtxl
+
+begin
+	switch (datatype) {
+	case TY_SHORT, TY_INT, TY_LONG, TY_USHORT:
+	    call amovl (Meml[buf], outbuffer, npix)
+	case TY_REAL:
+	    call achtrl (Memr[buf], outbuffer, npix)
+	case TY_DOUBLE:
+	    call achtdl (Memd[buf], outbuffer, npix)
+	case TY_COMPLEX:
+	    call achtxl (Memx[buf], outbuffer, npix)
+	default:
+	    call error (13, "WFT_LONG_LINE: Unknown IRAF data type.")
+	}
+end
+
+
+# ALTALL -- Procedure to linearly scale a long vector into a long vector
+# using double precision constants to preserve precision.
+
+procedure altall (a, b, npix, k1, k2)
+
+long	a[ARB]		# input vector
+long	b[ARB]		# output vector
+int	npix		# number of pixels
+double	k1, k2		# scaling factors
+
+double	dtemp
+int	i
+
+begin
+	do i = 1, npix {
+	    dtemp = (a[i] + k1) * k2
+	    if (dtemp >= 0.0d0)
+		dtemp = dtemp + 0.5d0
+	    else
+		dtemp = dtemp - 0.5d0
+	    b[i] = dtemp
+	}
+end
+
+
+# ALTARL -- Procedure to linearly scale a real vector into a long vector
+# using double precision constants to preserve precision.
+
+procedure altarl (a, b, npix, k1, k2)
+
+real	a[ARB]		# input vector
+long	b[ARB]		# output vector
+int	npix		# number of pixels
+double	k1, k2		# scaling factors
+
+int	i
+double	dtemp
+
+begin
+	do i = 1, npix {
+	    dtemp = (a[i] + k1) * k2
+	    if (dtemp >= 0.0d0)
+		dtemp = dtemp + 0.5d0
+	    else
+		dtemp = dtemp - 0.5d0
+	    b[i] = dtemp
+	}
+end
+
+
+# ALTADL -- Procedure to linearly scale a double vector into a long vector
+# using double precision constants to preserve precision.
+
+procedure altadl (a, b, npix, k1, k2)
+
+double	a[ARB]		# input vector
+long	b[ARB]		# output vector
+int	npix		# number of pixels
+double	k1, k2		# scaling factors
+
+int	i
+double	dtemp
+
+begin
+	do i = 1, npix {
+	    dtemp = (a[i] + k1) * k2
+	    if (dtemp >= 0.0d0)
+		dtemp = dtemp + 0.5d0
+	    else
+		dtemp = dtemp - 0.5d0
+	    b[i] = dtemp
+	}
+end
+
+
+# ALTAXL -- Procedure to linearly scale a complex vector into a long vector
+# using double precision constants to preserve precision.
+
+procedure altaxl (a, b, npix, k1, k2)
+
+complex	a[ARB]		# input vector
+long	b[ARB]		# output vector
+int	npix		# number of pixels
+double	k1, k2		# scaling factors
+
+int	i
+double	dtemp
+
+begin
+	do i = 1, npix {
+	    dtemp = (a[i] + k1) * k2
+	    if (dtemp >= 0.0d0)
+		dtemp = dtemp + 0.5d0
+	    else
+		dtemp = dtemp - 0.5d0
+	    b[i] = dtemp
+	}
+end
+
+
+# ALTADR -- Procedure to linearly scale a real vector in double precision
+
+procedure altadr (a, b, npix, k1, k2)
+
+real	a[ARB]		# input vector
+real	b[ARB]		# output vector
+int	npix		# number of pixels
+double	k1, k2		# scaling factors
+
+int	i
+
+begin
+	do i = 1, npix
+	    b[i] = (a[i] + k1) * k2
+end
+
+
+# ALTADX -- Procedure to linearly scale a complex vector in double precision
+
+procedure altadx (a, b, npix, k1, k2)
+
+complex	a[ARB]		# input vector
+complex	b[ARB]		# output vector
+int	npix		# number of pixels
+double	k1, k2		# scaling factors
+
+int	i
+
+begin
+	do i = 1, npix
+	    b[i] = (a[i] + k1) * k2
+end
+
diff --git a/pkg/dataio/fits/fits_wpixels.x b/pkg/dataio/fits/fits_wpixels.x
new file mode 100644
index 00000000..7a9389ac
--- /dev/null
+++ b/pkg/dataio/fits/fits_wpixels.x
@@ -0,0 +1,162 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include <fset.h>
+include "wfits.h"
+
+# WFT_INIT_WRITE_PIXELS -- This procedure calculates the input and
+# output buffer sizes based in the spp and mii data types and allocates
+# the required space.
+
+procedure wft_init_write_pixels (npix_record, spp_type, bitpix, blkfac)
+
+int	npix_record		# number of data pixels per record
+int	spp_type		# pixel data type
+int	bitpix			# output bits per pixel
+int	blkfac			# fits blocking factor (0 for disk)
+
+# entry wft_write_pixels, wft_write_last_record
+
+int	fd			# output file descriptor
+char	buffer[1]		# input buffer
+int	npix			# number of pixels in the input buffer
+int	nrecords		# number of FITS records written
+
+char	blank, zero
+int	ty_mii, ty_spp, npix_rec, nch_rec, len_mii, sz_rec, nchars, n, nrec
+int	bf, szblk
+pointer	spp, mii, ip, op
+
+int	sizeof(), miilen(), fstati()
+long	note()
+errchk	malloc, mfree, write, miipak, amovc
+data	mii /NULL/, spp/NULL/
+
+begin
+	# Change input parameters into local variables.
+	ty_mii = bitpix
+	ty_spp = spp_type
+	npix_rec = npix_record
+	nch_rec = npix_rec * sizeof (ty_spp)
+	bf = blkfac
+	blank = ' '
+	zero = 0
+
+	# Compute the size of the mii buffer.
+	len_mii = miilen (npix_rec, ty_mii)
+	sz_rec = len_mii * SZ_INT32
+
+	# Allocate space for the buffers.
+	if (spp != NULL)
+	    call mfree (spp, TY_CHAR)
+	call malloc (spp, nch_rec, TY_CHAR)
+	if (mii != NULL)
+	    call mfree (mii, TY_INT)
+	call malloc (mii, len_mii, TY_INT)
+
+	op = 0
+	nrec = 0
+
+	return
+
+# WFT_WRITE_PIXELS -- Wft_wrt_pixels gets an image line and places it in the
+# output buffer. When the output buffer is full the data are packed by the mii
+# routines and written to the specified output.
+
+entry	wft_write_pixels (fd, buffer, npix)
+
+	nchars = npix * sizeof (ty_spp)
+	ip = 0
+
+	repeat {
+
+	    # Fill output buffer.
+	    n = min (nch_rec - op, nchars - ip)
+	    call amovc (buffer[1 + ip], Memc[spp + op], n)
+	    ip = ip + n
+	    op = op + n
+
+	    # Write output record.
+	    if (op == nch_rec) {
+		call miipak (Memc[spp], Memi[mii], npix_rec, ty_spp, ty_mii)
+		iferr (call write (fd, Memi[mii], sz_rec)) {
+	            if (ty_spp == TY_CHAR) {
+		        call printf (" File incomplete: %d logical header")
+			    call pargi (nrec)
+			call printf (" (2880 byte) records written\n")
+	                call error (18,
+			    "WRT_RECORD: Error writing header record.")
+	            } else {
+		        call printf (" File incomplete: %d logical data")
+			    call pargi (nrec)
+			call printf (" (2880 byte) records written\n")
+		        call error (19,
+			    "WRT_RECORD: Error writing data record.")
+		    }
+		}
+
+		nrec = nrec + 1
+		op = 0
+	    }
+
+	} until (ip == nchars)
+
+	return
+
+
+# WFT_WRITE_LAST_RECORD -- Procedure to write the last partially filled record
+# to tape. Fill with blanks if header record otherwise fill with zeros.
+
+entry	wft_write_last_record (fd, nrecords)
+
+	if (op != 0) {
+
+	    # Blank or zero fill the last record.
+	    n = nch_rec - op
+	    if (ty_spp == TY_CHAR)
+		call amovkc (blank, Memc[spp + op], n)
+	    else
+		call amovkc (zero, Memc[spp + op], n)
+
+	    # Write last record.
+	    call miipak (Memc[spp], Memi[mii], npix_rec, ty_spp, ty_mii)
+	    iferr (call write (fd, Memi[mii], sz_rec)) {
+	        if (ty_spp == TY_CHAR) {
+		    call printf ("File incomplete: %d logical header")
+			call pargi (nrec)
+		    call printf (" (2880 byte) records written\n")
+	            call error (18,
+			"WRT_LAST_RECORD: Error writing last header record.")
+	        } else {
+		    call printf ("File incomplete: %d logical data")
+			call pargi (nrec)
+		    call printf (" (2880 byte) records written\n")
+		    call error (19,
+			"WRT_LAST_RECORD: Error writing last data record.")
+		}
+	    }
+
+
+	    nrec = nrec + 1
+
+	    # Pad out the record if the blocking is non-standard.
+	    szblk = fstati (fd, F_BUFSIZE) * SZB_CHAR
+	    if ((bf > 0) && mod (szblk, FITS_RECORD) != 0 &&
+	        (ty_spp != TY_CHAR)) {
+		szblk = szblk / SZB_CHAR
+		n = note (fd) - 1
+		if (mod (n, szblk) == 0)
+		    n = 0
+		else
+		    n = szblk - mod (n, szblk)
+		for (op = 1; op <= n; op = op + nch_rec) {
+		    szblk = min (nch_rec, n - op + 1)
+		    call amovkc (zero, Memc[spp], szblk)
+		    #call write (fd, Memc[spp], szblk)
+		}
+	    }
+
+	}
+
+	nrecords = nrec
+end
diff --git a/pkg/dataio/fits/fits_write.x b/pkg/dataio/fits/fits_write.x
new file mode 100644
index 00000000..edfc9f83
--- /dev/null
+++ b/pkg/dataio/fits/fits_write.x
@@ -0,0 +1,246 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <fset.h>
+include	<error.h>
+include <mach.h>
+include <imhdr.h>
+include "wfits.h"
+
+# WFT_WRITE_FITZ -- Procedure to convert a single IRAF file to a FITS file.
+# If the make_image switch is set the header and pixel files are output
+# to the output destination. If the make_image switch is off the header
+# is printed to the standard output.
+
+procedure wft_write_fitz (iraf_file, fits_file, image_number, nimages)
+
+char	iraf_file[ARB]		# IRAF file name
+char	fits_file[ARB]		# FITS file name
+int	image_number		# the current image number
+int	nimages			# the number of images
+
+int	fits_fd, chars_rec, nchars, ip, min_lenuserarea
+pointer	im, sp, fits, envstr
+
+int	mtfile(), mtopen(), open(), fnldir(), envfind(), ctoi()
+pointer	immap()
+errchk	immap, imunmap, open, mtopen, close, smark, salloc, sfree
+errchk	delete, wft_write_header, wft_write_image, wft_data_limits
+
+include "wfits.com"
+
+begin
+	# Open the output file. Check whether the output file is a magtape
+	# device or a binary file. If the output file is magtape check
+	# for a legal blocking factor.
+
+	if (image_number == 1 || wextensions == NO) {
+	    if (make_image == NO)
+	        call strcpy ("dev$null", fits_file, SZ_FNAME)
+	    if (mtfile (fits_file) == YES) {
+	        chars_rec = (blkfac * len_record * FITS_BYTE) / (SZB_CHAR *
+	            NBITS_BYTE)
+	        fits_fd = mtopen (fits_file, WRITE_ONLY, chars_rec)
+	    } else
+	        fits_fd = open (fits_file, NEW_FILE, BINARY_FILE)
+	}
+
+	# Allocate memory for program data structure.
+	call smark (sp)
+	call salloc (fits, LEN_FITS, TY_STRUCT)
+	call salloc (envstr, SZ_FNAME, TY_CHAR)
+
+	# Set up the minimum length of the user area.
+	if (envfind ("min_lenuserarea", Memc[envstr], SZ_FNAME) > 0) {
+	    ip = 1
+	    if (ctoi (Memc[envstr], ip, min_lenuserarea) <= 0)
+		min_lenuserarea = LEN_USERAREA
+	    else
+		min_lenuserarea = max (LEN_USERAREA, min_lenuserarea)
+	} else
+	    min_lenuserarea = LEN_USERAREA
+
+	# Write the global header.
+	if (image_number == 1 && gheader == YES) {
+
+	    XTENSION(fits) = EXT_PRIMARY
+
+	    # Open a dummy image.
+	    im = immap ("dev$null", NEW_IMAGE, 0)
+	    NAXIS(im) = 0
+	    PIXTYPE(im) = TY_SHORT
+	    OBJECT(im) = EOS
+	    IRAFNAME(fits) = EOS
+
+            if (long_header == YES || short_header == YES) {
+                call printf ("Global header")
+		if (make_image == YES) {
+		    call printf (" -> %s[0] ")
+		        call pargstr (fits_file)
+		}
+		if (long_header == YES)
+		    call printf ("\n")
+		else if (short_header == YES)
+		    call printf (" ")
+            }
+	    call flush (STDOUT)
+
+	    iferr {
+		call wft_write_header (im, fits, fits_fd)
+	        if (make_image == YES)
+		    call wft_write_image (im, fits, fits_fd)
+	    } then {
+
+	        # Print the error message.
+	        call flush (STDOUT)
+	        call erract (EA_WARN)
+
+	        # Close files and cleanup.
+	        call imunmap (im)
+	        #if (image_number == nimages || wextensions == NO)
+	            call close (fits_fd)
+	        if (make_image == NO)
+	            call delete (fits_file)
+	        call sfree (sp)
+
+	        # Assert an error.
+	        call erract (EA_ERROR)
+
+	    } else {
+		call imunmap (im)
+	    }
+
+	    if (long_header == YES)
+	        call printf ("\n")
+	}
+
+	# Map the input image. Construct the old iraf name by removing
+	# the directory specification.
+	# Print the id string.
+        if (long_header == YES || short_header == YES) {
+            call printf ("Image %d: %s")
+                call pargi (image_number)
+                call pargstr (iraf_file)
+        }
+
+	# Define whether the image to be written is to be the FITS primary
+	# data image file or a FITS image extension file.
+	if (image_number == 1) {
+	    if (wextensions == YES && gheader == YES)
+	        XTENSION(fits) = EXT_IMAGE
+	    else 
+	        XTENSION(fits) = EXT_PRIMARY
+	} else {
+	    if (wextensions == YES)
+	        XTENSION(fits) = EXT_IMAGE
+	    else
+	        XTENSION(fits) = EXT_PRIMARY
+	}
+
+	im = immap (iraf_file, READ_ONLY, min_lenuserarea)
+	call imgcluster (iraf_file, IRAFNAME(fits), SZ_FNAME)
+	nchars = fnldir (IRAFNAME(fits), IRAFNAME(fits), SZ_FNAME)
+	call strcpy (iraf_file[nchars+1], IRAFNAME(fits), SZ_FNAME)
+
+	# Write header and image.
+	iferr {
+
+	    if (short_header == YES || long_header == YES) {
+	        if (make_image == YES) {
+		    if (wextensions == YES && nimages > 1) {
+		        call printf (" -> %s[%d] ")
+		            call pargstr (fits_file)
+			    if (gheader == YES)
+			        call pargi (image_number)
+			    else
+			        call pargi (image_number - 1)
+		    } else {
+		        call printf (" -> %s ")
+		            call pargstr (fits_file)
+		    }
+	        }
+		if (long_header == YES)
+		    call printf ("\n")
+		else if (short_header == YES)
+		    call printf (" ")
+	    }
+	    call flush (STDOUT)
+
+	    call wft_write_header (im, fits, fits_fd)
+	    if (make_image == YES)
+		call wft_write_image (im, fits, fits_fd)
+
+	    if (long_header == YES)
+	        call printf ("\n")
+
+	} then {
+
+	    # Print the error message.
+	    call flush (STDOUT)
+	    call erract (EA_WARN)
+
+	    # Close files and cleanup.
+	    call imunmap (im)
+	    #if (image_number == nimages || wextensions == NO)
+	        call close (fits_fd)
+	    if (make_image == NO)
+	        call delete (fits_file)
+	    call sfree (sp)
+
+	    # Assert an error.
+	    call erract (EA_ERROR)
+
+	} else {
+
+	    # Close files and cleanup.
+	    call imunmap (im)
+	    if (image_number == nimages || wextensions == NO)
+	        call close (fits_fd)
+	    if (make_image == NO)
+	        call delete (fits_file)
+	    call sfree (sp)
+	}
+
+end
+
+
+# WFT_DATA_LIMITS -- Procedure to calculate the maximum and minimum data values
+# in an IRAF image. Values are only calculated if the max and min are unknown
+# or the image has been modified since the last values were calculated.
+
+procedure wft_data_limits (im, irafmin, irafmax)
+
+pointer	im		# image pointer
+real	irafmin		# minimum picture value
+real	irafmax		# maximum picture value
+
+int	npix
+long	v[IM_MAXDIM]
+pointer	buf
+real	maxval, minval
+int	imgnlr()
+errchk	imgnlr
+
+begin
+	# Compute the data minimum and maximum if the image values
+	# are undefined out-of-date.
+
+	if (LIMTIME(im) < MTIME(im) && NAXIS(im) > 0) {
+
+	    irafmax = -MAX_REAL
+	    irafmin = MAX_REAL
+	    npix = NAXISN(im,1)
+
+	    call amovkl (long(1), v, IM_MAXDIM)
+	    while (imgnlr (im, buf, v) != EOF) {
+	        call alimr (Memr[buf], npix, minval, maxval)
+	        irafmin = min (irafmin, minval)
+	        irafmax = max (irafmax, maxval)
+	    }
+
+	} else {
+
+	    irafmax = IM_MAX(im)
+	    irafmin = IM_MIN(im)
+
+	}
+end
diff --git a/pkg/dataio/fits/mkpkg b/pkg/dataio/fits/mkpkg
new file mode 100644
index 00000000..ac5201d0
--- /dev/null
+++ b/pkg/dataio/fits/mkpkg
@@ -0,0 +1,24 @@
+# Make the  RFITS / WFITS Tasks
+
+$checkout	libpkg.a ../
+$update		libpkg.a
+$checkin	libpkg.a ../
+$exit
+
+libpkg.a:
+	fits_cards.x	wfits.com wfits.h <imhdr.h>
+	fits_params.x	wfits.h <time.h>
+	fits_read.x	rfits.com rfits.h <error.h> <fset.h> <imhdr.h>\
+			<mach.h> <plset.h>
+	fits_rheader.x	rfits.com rfits.h <ctype.h> <imhdr.h> <imio.h>\
+			<mach.h>
+	fits_rimage.x	rfits.com rfits.h <imhdr.h> <fset.h> <mach.h>
+	fits_rpixels.x	<fset.h> <mach.h> <mii.h>
+	fits_wheader.x	wfits.com wfits.h <fset.h> <imhdr.h> <mach.h>
+	fits_wimage.x	wfits.com wfits.h <error.h> <imhdr.h>
+	fits_wpixels.x	wfits.h <fset.h> <mach.h>
+	fits_write.x	<error.h> wfits.com wfits.h <fset.h> <imhdr.h> <mach.h>
+	fits_files.x	<ctype.h> <plset.h>
+	t_rfits.x	rfits.com rfits.h <error.h> <fset.h>
+	t_wfits.x	wfits.com wfits.h <error.h> <fset.h> <mach.h>
+	;
diff --git a/pkg/dataio/fits/rfits.com b/pkg/dataio/fits/rfits.com
new file mode 100644
index 00000000..08f44c0e
--- /dev/null
+++ b/pkg/dataio/fits/rfits.com
@@ -0,0 +1,18 @@
+
+# FITS reader common
+
+int	len_record		# Length of FITS records in bytes
+int	data_type		# Output data type
+real	blank			# Blank value
+real	fe			# Maximum size in megabytes for scan mode
+
+# Option flags
+
+int	make_image		# Create an IRAF image
+int	long_header		# Print a long header (FITS header cards)
+int	short_header		# Print a short header (Title and size)
+int	scale			# Scale the data
+int	old_name		# Use old IRAF name?
+
+common	/rfitscom/ len_record, data_type, blank, fe, make_image, long_header,
+	short_header, scale, old_name
diff --git a/pkg/dataio/fits/rfits.h b/pkg/dataio/fits/rfits.h
new file mode 100644
index 00000000..bab29220
--- /dev/null
+++ b/pkg/dataio/fits/rfits.h
@@ -0,0 +1,96 @@
+# FITS Definitions
+
+# The FITS standard readable by the FITS reader using these definitions:
+# 
+# 1.  8 bits / byte
+# 2.  ASCII character code
+# 3.  MII integer data format (i.e. 8 bit unsigned integers and 16 and 32
+#     bit signed twos complement integers with most significant bytes first.)
+# 4.  IEEE 32 and 64 bit floating point format
+#
+#
+# The following deviations from the FITS standard are allowed:
+# 
+# 1. The number of FITS bytes per record is normally 2880 or up to 10 times
+#    2880 bytes but may be arbitrarily specified  by the user.
+
+# Define the bits per pixel, precision and byte order of the basic FITS types
+
+define	FITS_RECORD	2880	# number of bytes in a standard FITS record
+
+define	FITS_BYTE	 8	# Bits in a FITS byte
+define	FITS_SHORT	16	# Bits in a FITS short
+define	FITS_LONG	32	# Bits in a FITS long
+define	FITS_REAL      -32	# Bits in a FITS real * -1
+define	FITS_DOUBLE    -64	# Bits in a FITS double * -1
+
+define	FITSB_PREC	 3	# Decimal digits of precision in a FITS byte
+define	FITSS_PREC	 5	# Decimal digits of precision in a FITS short
+define	FITSL_PREC	10	# Decimal digits of precision in a FITS long
+
+define	LSBF		NO	# Least significant byte first
+
+# Define the basic format of a FITS cardimage
+
+define	LEN_CARD	80	# Length of FITS card in characters
+define	COL_VALUE	11	# Starting column for parameter values
+
+
+# FITS standards not recognized currently by IRAF.
+#
+# 1. SIMPLE	SIMPLE = 'F' not implemented, file skipped
+# 2. GROUPS	Group data not currently implemented, file skippped
+
+# FITS extension currently recognised by IRAF
+
+define	EXT_PRIMARY	1		# recognized and read
+define	EXT_IMAGE	2		# recognized and read
+define	EXT_TABLE	3		# recognized and skipped
+define	EXT_BINTABLE	4		# recognized and skipped
+define	EXT_UNKNOWN	5		# unrecognized and skipped
+define	EXT_SPECIAL	6		# undefined
+
+
+# Values for the following quantities are stored in the structure below.
+
+define	LEN_FITS	(20 + SZ_FNAME + 1)
+
+define	FITS_BSCALE	Memd[P2D($1)]	  # FITS scaling parameter
+define	FITS_BZERO	Memd[P2D($1+2)]   # FITS zero parameter
+define	BLANK_VALUE	Meml[P2L($1+4)]   # Blank value
+define	BLANKS		Memi[$1+5]	  # YES if blank keyword in header
+define	BITPIX		Memi[$1+6]	  # Bits per pixel (Must be an MII type)
+define	SCALE		Memi[$1+7]	  # Scale the data ?
+define	SIMPLE		Memi[$1+8]	  # Standard FITS format 
+define	NRECORDS	Memi[$1+9]	  # Number of FITS logical records
+define	EXTEND		Memi[$1+10]	  # FITS extensions may be present
+define	XTENSION	Memi[$1+11]	  # FITS extension type
+define	PCOUNT		Memi[$1+12]	  # Number of random parameters
+define	GCOUNT		Memi[$1+13]	  # Number of groups
+define	GLOBALHDR	Memi[$1+14]	  # Global header may be present
+define	INHERIT		Memi[$1+15]	  # Inherit global header if present
+define	IRAFNAME	Memc[P2C($1+16)]  # Old IRAF name
+
+# Mapping of additional IRAF header parameters
+
+define	PIXTYPE		IM_PIXTYPE($1)
+define	NBPIX		IM_NBPIX($1)
+define	IRAFMAX		IM_MAX($1)
+define	IRAFMIN		IM_MIN($1)
+define	LIMTIME		IM_LIMTIME($1)
+define	LEN_USERAREA	28800
+
+# Mapping of FITS Keywords to IRAF image header
+
+define	NAXIS	 IM_NDIM($1)
+define	NAXISN	 IM_LEN($1,$2)
+define	OBJECT	 IM_TITLE($1)
+define	HISTORY	 IM_HISTORY($1)
+define	UNKNOWN	 Memc[($1+IMU-1)*SZ_STRUCT+1]	# All unrecognized keywords
+						# are stored here
+# Miscellaneous definitions.
+
+define	SZ_OBJECT	SZ_IMTITLE
+define	SZ_HISTORY	SZ_IMHIST
+define	SZ_FCTYPE	SZ_CTYPE
+define	LEN_TYPESTR	8
diff --git a/pkg/dataio/fits/t_rfits.x b/pkg/dataio/fits/t_rfits.x
new file mode 100644
index 00000000..06c55ec1
--- /dev/null
+++ b/pkg/dataio/fits/t_rfits.x
@@ -0,0 +1,216 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+include <fset.h>
+include "rfits.h"
+
+define	NTYPES		7		# the number of image data types
+
+# RFITS -- Read FITS format data.  Further documentation given in rfits.hlp
+
+procedure t_rfits()
+
+int	inlist, outlist, len_inlist, len_outlist
+int	file_number, offset, stat, first_file, last_file
+pointer	sp, infile, file_list, outfile, ext_list, in_fname, out_fname
+pointer	pl, axes
+
+bool	clgetb(), pl_linenotempty()
+#char	clgetc()
+int	rft_get_image_type(), clgeti(), mtfile(), strlen(), btoi(), fntlenb()
+int	rft_read_fitz(), fntgfnb(), fstati(), mtneedfileno(), fntrfnb()
+pointer	fntopnb(), rft_flist()
+real	clgetr(), rft_fe()
+
+include	"rfits.com"
+
+begin
+	# Set up the standard output to flush on a newline.
+	if (fstati (STDOUT, F_REDIR) == NO)
+	    call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (infile, SZ_FNAME, TY_CHAR)
+	call salloc (file_list, SZ_LINE, TY_CHAR)
+	call salloc (outfile, SZ_FNAME, TY_CHAR)
+	call salloc (ext_list, SZ_LINE, TY_CHAR)
+	call salloc (in_fname, SZ_FNAME, TY_CHAR)
+	call salloc (out_fname, SZ_FNAME, TY_CHAR)
+	call salloc (axes, 2, TY_INT)
+
+	# Get RFITS parameters.
+	call clgstr ("fits_file", Memc[infile], SZ_FNAME)
+	long_header = btoi (clgetb ("long_header"))
+	short_header = btoi (clgetb ("short_header"))
+	len_record = FITS_RECORD
+	old_name = btoi (clgetb ("oldirafname"))
+	make_image = btoi (clgetb ("make_image"))
+
+	# Open the input file list.
+	call clgstr ("file_list", Memc[ext_list], SZ_LINE)
+	if (mtfile (Memc[infile]) == YES) {
+	    inlist = NULL
+	    if (mtneedfileno (Memc[infile]) == YES) {
+	        call strcpy (Memc[ext_list], Memc[file_list], SZ_LINE)
+	    } else {
+		call sprintf (Memc[file_list], SZ_LINE, "1[%s]")
+		    call pargstr (Memc[ext_list])
+	    }
+	} else {
+	    inlist = fntopnb (Memc[infile], NO)
+	    len_inlist = fntlenb (inlist)
+	    if (len_inlist > 0) {
+		if (Memc[ext_list] == EOS) {
+	            call sprintf  (Memc[file_list], SZ_LINE, "1-%d[0]")
+		        call pargi (len_inlist)
+		        #call pargstr (Memc[ext_list])
+		} else {
+	            call sprintf  (Memc[file_list], SZ_LINE, "1-%d[%s]")
+		        call pargi (len_inlist)
+		        call pargstr (Memc[ext_list])
+		}
+	    } else {
+	        call sprintf  (Memc[file_list], SZ_LINE, "0[%s]")
+		    call pargstr (Memc[ext_list])
+	    }
+	}
+
+	# Decode the ranges string.
+	pl = rft_flist (Memc[file_list], first_file, last_file, len_inlist) 
+	if (pl == NULL || len_inlist <= 0)
+	    call error (1, "T_RFITS: Illegal file/extensions number list")
+
+	# Open the output file list.
+	if (make_image == YES) {
+	    call clgstr ("iraf_file", Memc[outfile], SZ_FNAME)
+	    if (Memc[outfile] == EOS) {
+		if (old_name == YES)
+		    call mktemp ("tmp$", Memc[outfile], SZ_FNAME) 
+		else
+		    call error (0, "T_RFITS: Undefined output file name")
+	    }
+	    outlist = fntopnb (Memc[outfile], NO)
+	    len_outlist = fntlenb (outlist)
+	    offset = clgeti ("offset")
+	} else {
+	    Memc[outfile] = EOS
+	    outlist = NULL
+	    len_outlist = 1
+	}
+	if ((len_outlist > 1) && (len_outlist != len_inlist))
+	    call error (0,
+	        "T_RFITS: Output and input lists have different lengths")
+
+	# Get the remaining parameters. Use the string in_fname as a
+	# temporary variable.
+	#data_type = rft_get_image_type (clgetc ("datatype"))
+	call clgstr ("datatype", Memc[in_fname], SZ_FNAME)
+	data_type = rft_get_image_type (Memc[in_fname])
+	scale = btoi (clgetb ("scale"))
+	blank = clgetr ("blank")
+
+	# Get the scan size parameter.
+	fe = rft_fe (Memc[infile])
+
+	# Read successive FITS files, convert and write into a numbered
+	# succession of output IRAF files.
+
+	do file_number = first_file, last_file {
+
+	    # Get the next file number.
+	    Memi[axes] = 1
+	    Memi[axes+1] = file_number
+	    if (! pl_linenotempty (pl, Memi[axes]))
+		next
+
+	    # Get the input file name.
+	    if (inlist != NULL) {
+		if (fntgfnb (inlist, Memc[in_fname], SZ_FNAME) == EOF)
+		    call error (0, "T_RFITS: Error reading input file name")
+	    } else {
+		if (mtneedfileno (Memc[infile]) == YES)
+		    call mtfname (Memc[infile], file_number, Memc[in_fname],
+		        SZ_FNAME)
+		else
+	            call strcpy (Memc[infile], Memc[in_fname], SZ_FNAME)
+	    }
+
+	    # Get the output file name.
+	    if (outlist == NULL) {
+		Memc[out_fname] = EOS
+	    } else if (len_inlist > len_outlist) {
+		if (fntrfnb (outlist, 1, Memc[out_fname], SZ_FNAME) == EOF)
+	            call strcpy (Memc[outfile], Memc[out_fname], SZ_FNAME)
+	        if (len_inlist > 1) {
+		    call sprintf (Memc[out_fname+strlen(Memc[out_fname])],
+		        SZ_FNAME, "%04d")
+		        call pargi (file_number + offset)
+	        }
+	    } else if (fntgfnb (outlist, Memc[out_fname], SZ_FNAME) == EOF)
+		call error (0, "T_RFITS: Error reading output file name")
+
+	    # Convert FITS file to the output IRAF file. If EOT is reached
+	    # then exit. If an error is detected then print a warning and
+	    # continue with the next file.
+
+	    iferr (stat = rft_read_fitz (Memc[in_fname], Memc[out_fname],
+	        pl, file_number))
+		call erract (EA_FATAL)
+	    if (stat == EOF)
+		break
+	}
+
+	if (inlist != NULL)
+	    call fntclsb (inlist) 
+	if (outlist != NULL)
+	    call fntclsb (outlist) 
+	if (pl != NULL)
+	    call pl_close (pl)
+
+	call sfree (sp)
+end
+
+
+# RFT_GET_IMAGE_TYPE -- Convert a character to and IRAF image type.
+
+int procedure rft_get_image_type (c)
+
+char	c
+
+int	type_codes[NTYPES], i
+string	types "usilrdx"	
+int	stridx()
+data	type_codes /TY_USHORT, TY_SHORT, TY_INT, TY_LONG, TY_REAL,
+		    TY_DOUBLE, TY_COMPLEX/
+begin
+	i = stridx (c, types)
+	if (i == 0)
+	    return (ERR)
+	else
+	    return (type_codes[stridx(c,types)])
+end
+
+
+# RFT_FE -- Fetch the maximum file size in MB for tape scanning mode.
+
+real procedure rft_fe (file)
+
+char	file[ARB]		# the input file name
+
+pointer	gty
+real	fe
+int	mtfile(), gtygeti()
+pointer	mtcap()
+errchk	gtygeti()
+
+begin
+	if (mtfile (file) == NO)
+	    return (0.0)
+	iferr (gty = mtcap (file))
+	    return (0.0)
+	iferr (fe = gtygeti (gty, "fe"))
+	    fe = 0.0
+	call gtyclose (gty)
+	return (fe)
+end
diff --git a/pkg/dataio/fits/t_wfits.x b/pkg/dataio/fits/t_wfits.x
new file mode 100644
index 00000000..256b9cde
--- /dev/null
+++ b/pkg/dataio/fits/t_wfits.x
@@ -0,0 +1,253 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include <error.h>
+include <fset.h>
+include "wfits.h"
+
+# T_WFITS -- This procedure converts a series of IRAF image files to
+# FITS image files.
+
+procedure t_wfits ()
+
+char	iraf_files[SZ_FNAME]	# list of IRAF images
+char	fits_files[SZ_FNAME]	# list of FITS files
+bool	newtape			# new or used tape ?
+char	in_fname[SZ_FNAME]	# input file name
+char	out_fname[SZ_FNAME]	# output file name
+char	fextn[SZ_FNAME]		# the fits extension
+
+char	ch
+int	imlist, flist, nimages, nfiles, file_number, addext, index
+bool	clgetb(), streq()
+double	clgetd()
+int	imtopen(), imtlen (), wft_get_bitpix(), clgeti(), imtgetim()
+int	mtfile(), btoi(), fstati(), fntlenb(), fntgfnb(), mtneedfileno()
+int	wft_blkfac(), fntrfnb(), strlen(), strldx()
+pointer	fntopnb()
+
+include "wfits.com"
+
+begin
+	# Flush on a newline if STDOUT has not been redirected.
+	if (fstati (STDOUT, F_REDIR) == NO)
+	    call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Open iraf_files template and determine number of files in list.
+	call clgstr ("iraf_files", iraf_files, SZ_FNAME)
+	imlist = imtopen (iraf_files)
+	nimages = imtlen (imlist)
+
+	# Get the wfits parameters.
+	if (nimages == 1)
+	    wextensions = NO
+	else
+	    wextensions = btoi (clgetb ("extensions"))
+	if (wextensions == NO)
+	    gheader = NO
+	else
+	    gheader = btoi (clgetb ("global_hdr"))
+	long_header = btoi (clgetb ("long_header"))
+	short_header = btoi (clgetb ("short_header"))
+	make_image = btoi (clgetb ("make_image"))
+
+	# Get the FITS bits per pixel and the FITS logical record size.
+	bitpix = wft_get_bitpix (clgeti ("bitpix"))
+	len_record = FITS_RECORD
+
+	# Get the scaling parameters.
+	scale = btoi (clgetb ("scale"))
+	if (scale == YES) {
+	    if (clgetb ("autoscale"))
+		autoscale = YES
+	    else {
+	        bscale = clgetd ("bscale")
+	        bzero = clgetd ("bzero")
+		autoscale = NO
+	    }
+	} else {
+	    autoscale = NO
+	    bscale = 1.0d0
+	    bzero = 0.0d0
+	}
+
+	# Get the output file name and type (tape or disk). If no tape file
+	# number is given for output, the user is asked if the tape is blank
+	# or contains data. If the tape is blank output begins at BOT,
+	# otherwise at EOT.
+
+	call clgstr ("fextn", fextn, SZ_FNAME)
+	ch = '.'
+	if (make_image == YES) {
+	    call clgstr ("fits_files", fits_files, SZ_FNAME)
+	    if (mtfile (fits_files) == YES) {
+		flist = NULL
+		if (mtneedfileno (fits_files) == YES) {
+	            newtape = clgetb ("newtape")
+		    if (newtape)
+			call mtfname (fits_files, 1, out_fname, SZ_FNAME)
+		    else
+			call mtfname (fits_files, EOT, out_fname, SZ_FNAME)
+	        } else {
+		    call strcpy (fits_files, out_fname, SZ_FNAME)
+		    newtape = false
+		}
+	    } else {
+		flist = fntopnb (fits_files, NO)
+		nfiles = fntlenb (flist)
+		if (wextensions == YES && nfiles > 1)
+		    call error (0,
+			"Only one output FITS extensions file can be written")
+		if ((nfiles > 1) && (nfiles != nimages))
+		    call error (0,
+		    "T_WFITS: Input and output lists are not the same length")
+	    }
+	} else {
+	    fits_files[1] = EOS
+	    flist = NULL
+	}
+
+	# Get the fits file blocking factor.
+	blkfac = wft_blkfac (fits_files, clgeti ("blocking_factor"))
+
+	# Loop through the list of input images files.
+
+	file_number = 1
+	while (imtgetim (imlist, in_fname, SZ_FNAME) != EOF) {
+
+	    # Get the output file name. If single file output to disk, use
+	    # name fits_file. If multiple file output to disk, the file number
+	    # is added to the output file name, if no output name list is
+	    # supplied. If an output name list is supplied then the names
+	    # are extracted one by one from that list.
+
+	    if (make_image == YES) {
+	        if (mtfile (fits_files) == YES) {
+		    if (wextensions == NO && file_number == 2)
+			call mtfname (out_fname, EOT, out_fname, SZ_FNAME)
+	        } else if (nfiles > 1) {
+		    if (fntgfnb (flist, out_fname, SZ_FNAME) == EOF)
+			 call error (0, "Error reading output file name")
+		    if (fextn[1] != EOS) {
+			addext = OK
+			index = strldx (ch, out_fname)
+			if (index > 0) {
+			    if (streq (fextn, out_fname[index+1]))
+				addext = ERR
+			    else
+				addext = OK
+			}
+			if (addext == OK){
+		            call strcat (".", out_fname, SZ_FNAME)
+		            call strcat (fextn, out_fname, SZ_FNAME)
+			}
+		    }
+		} else {
+		    if (fntrfnb (flist, 1, out_fname, SZ_FNAME) == EOF)
+		        call strcpy (fits_files, out_fname, SZ_FNAME)
+		    if (nimages > 1 && wextensions == NO) {
+		        call sprintf (out_fname[strlen(out_fname)+1], SZ_FNAME,
+			    "%04d")
+			    call pargi (file_number)
+		    }
+		    if (fextn[1] != EOS) {
+			addext = OK
+			index = strldx (ch, out_fname)
+			if (index > 0) {
+			    if (streq (fextn, out_fname[index+1]))
+				addext = ERR
+			    else
+				addext = OK
+			}
+			if (addext == OK){
+		            call strcat (".", out_fname, SZ_FNAME)
+		            call strcat (fextn, out_fname, SZ_FNAME)
+			}
+		    }
+	        }
+	    }
+
+	    # Write each output file.
+	    iferr (call wft_write_fitz (in_fname, out_fname, file_number,
+		nimages)) {
+		call printf ("Error writing file: %s\n")
+		    call pargstr (out_fname)
+		call erract (EA_WARN)
+		break
+	    } else
+		file_number = file_number + 1
+	}
+
+	# Close up the input and output lists.
+	call clpcls (imlist)
+	if (flist != NULL)
+	    call fntclsb (flist)
+end
+
+
+# WFT_GET_BITPIX -- This procedure fetches the user determined bitpix or ERR if
+# the bitpix is not one of the permitted FITS types.
+
+int procedure wft_get_bitpix (bitpix)
+
+int	bitpix
+
+begin
+	switch (bitpix) {
+	case FITS_BYTE, FITS_SHORT, FITS_LONG, FITS_REAL, FITS_DOUBLE:
+	    return (bitpix)
+	default:
+	    return (ERR)
+	}
+end
+
+
+# WFT_BLKFAC -- Get the fits tape blocking factor.
+
+int procedure wft_blkfac (file, ublkfac)
+
+char	file[ARB]		# the input file name
+int	ublkfac			# the user supplied blocking factor
+
+int	bs, fb, blkfac
+pointer	gty
+int	mtfile(), mtcap(), gtygeti()
+errchk	mtcap(), gtygeti()
+
+begin
+	# Return a blocking factor of 1 if the file is a disk file.
+	if (mtfile (file) == NO)
+	    return (0)
+
+	# Open the tapecap device entry for the given device, and get
+	# the device block size and default FITS blocking factor
+	# parameters.
+
+	iferr (gty = mtcap (file))
+	    return (max (ublkfac,1))
+	iferr (bs = gtygeti (gty, "bs")) {
+	    call gtyclose (gty)
+	    return (max (ublkfac,1))
+	}
+	iferr (fb = max (gtygeti (gty, "fb"), 1))
+	    fb = 1
+
+	# Determine whether the device is a fixed or variable blocked
+	# device. Set the fits blocking factor to the value of the fb
+	# parameter if device is fixed block or if the user has
+	# requested the default blocking factor. Set the blocking factor
+	# to the user requested value if the device supports variable
+	# blocking factors.
+
+	if (bs == 0) {
+	    if (ublkfac <= 0)
+	        blkfac = fb
+	    else
+		blkfac = ublkfac
+	} else
+	    blkfac = fb
+
+	call gtyclose (gty)
+
+	return (blkfac)
+end
diff --git a/pkg/dataio/fits/wfits.com b/pkg/dataio/fits/wfits.com
new file mode 100644
index 00000000..04779ef3
--- /dev/null
+++ b/pkg/dataio/fits/wfits.com
@@ -0,0 +1,17 @@
+# FITS common block
+
+double	bscale		# FITS scaling factor
+double	bzero		# FITS offset factor
+int	bitpix		# Output bits per pixel
+int	len_record	# Record length in FITS bytes
+int	long_header	# Print long header?
+int	short_header	# Print short header?
+int	make_image	# Make a FITS image?
+int	scale		# Scale the data with bzero and bscale?
+int	autoscale	# Allow program to calculate bscale and bzero?
+int	blkfac		# FITS tape blocking factor
+int	wextensions	# Write a FITS extensions file
+int	gheader		# Write a global FITS extensions file header
+
+common /wfitscom/ bscale, bzero, bitpix, len_record, long_header, short_header,
+		  make_image, scale, autoscale, blkfac, wextensions, gheader
diff --git a/pkg/dataio/fits/wfits.h b/pkg/dataio/fits/wfits.h
new file mode 100644
index 00000000..a36caa89
--- /dev/null
+++ b/pkg/dataio/fits/wfits.h
@@ -0,0 +1,128 @@
+# WFITS header file
+
+# The basic FITS data structure
+
+define	LEN_FITS	(44 + SZ_FNAME + 1)
+
+define	BSCALE		Memd[P2D($1)]	  # FITS bscale value
+define	BZERO		Memd[P2D($1+2)]	  # FITS bzero value
+define	TAPEMAX		Memd[P2D($1+4)]	  # IRAF tape max
+define	TAPEMIN		Memd[P2D($1+6)]	  # IRAF tape min
+define	IRAFMAX		Memr[P2R($1+8)]	  # IRAF image maximum
+define	IRAFMIN		Memr[P2R($1+9)]	  # IRAF image minimum
+define	BLANK		Meml[P2L($1+10)]  # FITS blank value
+define	FITS_BITPIX	Memi[$1+11]	  # FITS bits per pixel
+define	DATA_BITPIX	Memi[$1+12]	  # Data bits per pixel
+define	SCALE		Memi[$1+13]	  # Scale data?
+define	XTENSION	Memi[$1+14]	  # FITS extension type
+define	BLANK_STRING	Memc[P2C($1+19)]  # String containing FITS blank value
+define	TYPE_STRING	Memc[P2C($1+31)]  # String containing IRAF type
+define	IRAFNAME	Memc[P2C($1+41)]  # IRAF file name
+
+
+# Define the FITS record size
+
+define	FITS_RECORD	2880		# Size of standard FITS record (bytes)
+
+# Define the supported FITS extensions
+
+define	EXT_PRIMARY	1		# the primary data array
+define	EXT_IMAGE	2		# the image extension
+
+# Define the FITS data types
+
+define	FITS_BYTE	  8		# Number of bits in a FITS byte
+define	FITS_SHORT	 16		# Number of bits in a FITS short
+define	FITS_LONG	 32		# Number of bits in a FITS long
+define	FITS_REAL	-32		# Number of bits in a FITS real * -1
+define	FITS_DOUBLE	-64		# Number of bits in a FITS double * -1
+
+# Define the FITS precision in decimal digits
+
+define	BYTE_PREC	 3		# Precision of FITS byte
+define	SHORT_PREC	 5		# Precision of FITS short
+define	LONG_PREC	10		# Precision of FITS long
+
+# Define the FITS blank data values
+
+define	BYTE_BLANK	 0.0d0			# Blank value for a FITS byte
+define	SHORT_BLANK	-3.2768d4		# Blank value for a FITS short
+define	LONG_BLANK	-2.147483648d9		# Blank value for a FITS long
+#define	BYTE_BLANK	 0			# Blank value for a FITS byte
+#define	SHORT_BLANK	-32768			# Blank value for a FITS short
+#define	LONG_BLANK	-2147483648		# Blank value for a FITS long
+
+# Define the FITS integer max and min values
+
+define	BYTE_MAX		  2.55d2	# Max value for a FITS byte
+define	BYTE_MIN		   1.0d0	# Min value for a FITS byte
+define	SHORT_MAX	        3.2767d4	# Max value for a FITS short
+define	SHORT_MIN              -3.2767d4	# Min value for a FITS short
+define	LONG_MAX	   2.147483647d9	# Max value for a FITS long
+define	LONG_MIN	  -2.147483647d9	# Min value for a FITS long
+define	PREC_RATIO	     .99978637d0	# Tape span reduction factor
+
+# Define the FITS card image parameters
+
+define	LEN_CARD		80	# Length of FITS header card
+define	LEN_KEYWORD		8	# Length of FITS keyword
+define	LEN_NAXIS_KYWRD		5       # Length of the NAXIS keyword string
+define	COL_VALUE		11	# First column of value field
+
+# Mapping of FITS task keywords to IRAF image header keywords
+
+define	NAXIS		IM_NDIM($1)		# Number of dimensions
+define	NAXISN		IM_LEN($1, $2)		# Length of each dimension
+define	OBJECT		IM_TITLE($1)		# Image title
+define	HISTORY		IM_HISTORY($1)		# History
+define	UNKNOWN	Memc[($1+IMU-1)*SZ_STRUCT+1]	# IRAF user area
+
+define	PIXTYPE		IM_PIXTYPE($1)		# Image pixel type
+define	NBPIX		IM_NBPIX($1)		# Number of bad pixels
+define	LIMTIME		IM_LIMTIME($1)		# Last modify limits time
+define	MTIME		IM_MTIME($1)		# Last modify time
+define	CTIME		IM_CTIME($1)		# Create time
+
+define	LEN_USERAREA	28800			# Default user area size
+
+# Set up a structure for the WFITS parameters
+
+# Define the required keywords
+
+define	FIRST_CARD	1		# FITS simple/xtension parameter
+define	SECOND_CARD	2		# FITS bitpix parameter
+define	THIRD_CARD	3		# FITS naxis parameter
+
+# Define the optional FITS KEYWORD parameters
+
+define	NOPTIONS	15		# Number of optional keywords
+
+define	KEY_EXTEND	1		# FITS EXTEND keyword
+define	KEY_PCOUNT	2		# Number of random parameter
+define	KEY_GCOUNT	3		# Number of groups
+define	KEY_BSCALE	4		# FITS bscale parameter
+define	KEY_BZERO	5		# FITS bzero parameter
+define	KEY_BUNIT	6		# FITS physical units
+define	KEY_BLANK	7		# FITS value of blank pixel
+define	KEY_OBJECT	8		# FITS title string
+define	KEY_ORIGIN	9		# Origin of FITS tape
+define	KEY_DATE	10		# Date the tape was written
+define	KEY_IRAFNAME	11		# Root name of IRAF image
+define	KEY_IRAFMAX	12		# Maximum value of IRAF image
+define	KEY_IRAFMIN	13		# Minimum value of IRAF image
+define	KEY_IRAFBP	14		# Bits per pixel in IRAF image
+define	KEY_IRAFTYPE	15		# IRAF image data type
+
+define	LEN_STRING	8		# Minimum length of a string parameter
+define	LEN_DATE	19		# Length of the new date string
+define	LEN_OBJECT	63		# Maximum length of string parameter
+define	LEN_ALIGN	18		# Maximum length for aligning parameter
+define	LEN_ORIGIN	9		# Length of origin string
+define	LEN_BLANK	11		# Length of the blank string
+define	NDEC_REAL	7		# Precision of real data
+define	NDEC_DOUBLE	11		# Precision of double precision data
+
+# Miscellaneous
+
+define	CENTURY		1900
+define	NEW_CENTURY	2000