aboutsummaryrefslogtreecommitdiff
path: root/pkg/dataio/fits/fits_read.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/dataio/fits/fits_read.x')
-rw-r--r--pkg/dataio/fits/fits_read.x469
1 files changed, 469 insertions, 0 deletions
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
generated by cgit (git 2.34.1) at 2025-09-20 08:23:48 -0400