aboutsummaryrefslogtreecommitdiff
path: root/pkg/dataio/fits/fits_read.x
blob: f9b0e46c4811b94026a5f201529a57078960fff7 (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
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