1 #include <string.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4
5 /*
6 Every program which uses the CFITSIO interface must include the
7 the fitsio.h header file. This contains the prototypes for all
8 the routines and defines the error status values and other symbolic
9 constants used in the interface.
10 */
11 #include "fitsio.h"
12
13 int main( void );
14 void writeimage( void );
15 void writeascii( void );
16 void writebintable( void );
17 void copyhdu( void );
18 void selectrows( void );
19 void readheader( void );
20 void readimage( void );
21 void readtable( void );
22 void printerror( int status);
23
main()24 int main()
25 {
26 /*************************************************************************
27 This is a simple main program that calls the following routines:
28
29 writeimage - write a FITS primary array image
30 writeascii - write a FITS ASCII table extension
31 writebintable - write a FITS binary table extension
32 copyhdu - copy a header/data unit from one FITS file to another
33 selectrows - copy selected row from one HDU to another
34 readheader - read and print the header keywords in every extension
35 readimage - read a FITS image and compute the min and max value
36 readtable - read columns of data from ASCII and binary tables
37
38 **************************************************************************/
39
40 writeimage();
41 writeascii();
42 writebintable();
43 copyhdu();
44 selectrows();
45 readheader();
46 readimage();
47 readtable();
48
49 printf("\nAll the cfitsio cookbook routines ran successfully.\n");
50 return(0);
51 }
52 /*--------------------------------------------------------------------------*/
writeimage(void)53 void writeimage( void )
54
55 /******************************************************/
56 /* Create a FITS primary array containing a 2-D image */
57 /******************************************************/
58 {
59 fitsfile *fptr; /* pointer to the FITS file, defined in fitsio.h */
60 int status, ii, jj;
61 long fpixel, nelements, exposure;
62 unsigned short *array[200];
63
64 /* initialize FITS image parameters */
65 char filename[] = "atestfil.fit"; /* name for new FITS file */
66 int bitpix = USHORT_IMG; /* 16-bit unsigned short pixel values */
67 long naxis = 2; /* 2-dimensional image */
68 long naxes[2] = { 300, 200 }; /* image is 300 pixels wide by 200 rows */
69
70 /* allocate memory for the whole image */
71 array[0] = (unsigned short *)malloc( naxes[0] * naxes[1]
72 * sizeof( unsigned short ) );
73
74 /* initialize pointers to the start of each row of the image */
75 for( ii=1; ii<naxes[1]; ii++ )
76 array[ii] = array[ii-1] + naxes[0];
77
78 remove(filename); /* Delete old file if it already exists */
79
80 status = 0; /* initialize status before calling fitsio routines */
81
82 if (fits_create_file(&fptr, filename, &status)) /* create new FITS file */
83 printerror( status ); /* call printerror if error occurs */
84
85 /* write the required keywords for the primary array image. */
86 /* Since bitpix = USHORT_IMG, this will cause cfitsio to create */
87 /* a FITS image with BITPIX = 16 (signed short integers) with */
88 /* BSCALE = 1.0 and BZERO = 32768. This is the convention that */
89 /* FITS uses to store unsigned integers. Note that the BSCALE */
90 /* and BZERO keywords will be automatically written by cfitsio */
91 /* in this case. */
92
93 if ( fits_create_img(fptr, bitpix, naxis, naxes, &status) )
94 printerror( status );
95
96 /* initialize the values in the image with a linear ramp function */
97 for (jj = 0; jj < naxes[1]; jj++)
98 { for (ii = 0; ii < naxes[0]; ii++)
99 {
100 array[jj][ii] = ii + jj;
101 }
102 }
103
104 fpixel = 1; /* first pixel to write */
105 nelements = naxes[0] * naxes[1]; /* number of pixels to write */
106
107 /* write the array of unsigned integers to the FITS file */
108 if ( fits_write_img(fptr, TUSHORT, fpixel, nelements, array[0], &status) )
109 printerror( status );
110
111 free( array[0] ); /* free previously allocated memory */
112
113 /* write another optional keyword to the header */
114 /* Note that the ADDRESS of the value is passed in the routine */
115 exposure = 1500;
116 if ( fits_update_key(fptr, TLONG, "EXPOSURE", &exposure,
117 "Total Exposure Time", &status) )
118 printerror( status );
119
120 if ( fits_close_file(fptr, &status) ) /* close the file */
121 printerror( status );
122
123 return;
124 }
125 /*--------------------------------------------------------------------------*/
writeascii(void)126 void writeascii ( void )
127
128 /*******************************************************************/
129 /* Create an ASCII table extension containing 3 columns and 6 rows */
130 /*******************************************************************/
131 {
132 fitsfile *fptr; /* pointer to the FITS file, defined in fitsio.h */
133 int status;
134 long firstrow, firstelem;
135
136 int tfields = 3; /* table will have 3 columns */
137 long nrows = 6; /* table will have 6 rows */
138
139 char filename[] = "atestfil.fit"; /* name for new FITS file */
140 char extname[] = "PLANETS_ASCII"; /* extension name */
141
142 /* define the name, datatype, and physical units for the 3 columns */
143 char *ttype[] = { "Planet", "Diameter", "Density" };
144 char *tform[] = { "a8", "I6", "F4.2" };
145 char *tunit[] = { "\0", "km", "g/cm" };
146
147 /* define the name diameter, and density of each planet */
148 char *planet[] = {"Mercury", "Venus", "Earth", "Mars","Jupiter","Saturn"};
149 long diameter[] = {4880, 12112, 12742, 6800, 143000, 121000};
150 float density[] = { 5.1f, 5.3f, 5.52f, 3.94f, 1.33f, 0.69f};
151
152 status=0;
153
154 /* open with write access the FITS file containing a primary array */
155 if ( fits_open_file(&fptr, filename, READWRITE, &status) )
156 printerror( status );
157
158 /* append a new empty ASCII table onto the FITS file */
159 if ( fits_create_tbl( fptr, ASCII_TBL, nrows, tfields, ttype, tform,
160 tunit, extname, &status) )
161 printerror( status );
162
163 firstrow = 1; /* first row in table to write */
164 firstelem = 1; /* first element in row (ignored in ASCII tables) */
165
166 /* write names to the first column (character strings) */
167 /* write diameters to the second column (longs) */
168 /* write density to the third column (floats) */
169
170 fits_write_col(fptr, TSTRING, 1, firstrow, firstelem, nrows, planet,
171 &status);
172 fits_write_col(fptr, TLONG, 2, firstrow, firstelem, nrows, diameter,
173 &status);
174 fits_write_col(fptr, TFLOAT, 3, firstrow, firstelem, nrows, density,
175 &status);
176
177 if ( fits_close_file(fptr, &status) ) /* close the FITS file */
178 printerror( status );
179 return;
180 }
181 /*--------------------------------------------------------------------------*/
writebintable(void)182 void writebintable ( void )
183
184 /*******************************************************************/
185 /* Create a binary table extension containing 3 columns and 6 rows */
186 /*******************************************************************/
187 {
188 fitsfile *fptr; /* pointer to the FITS file, defined in fitsio.h */
189 int status, hdutype;
190 long firstrow, firstelem;
191
192 int tfields = 3; /* table will have 3 columns */
193 long nrows = 6; /* table will have 6 rows */
194
195 char filename[] = "atestfil.fit"; /* name for new FITS file */
196 char extname[] = "PLANETS_Binary"; /* extension name */
197
198 /* define the name, datatype, and physical units for the 3 columns */
199 char *ttype[] = { "Planet", "Diameter", "Density" };
200 char *tform[] = { "8a", "1J", "1E" };
201 char *tunit[] = { "\0", "km", "g/cm" };
202
203 /* define the name diameter, and density of each planet */
204 char *planet[] = {"Mercury", "Venus", "Earth", "Mars","Jupiter","Saturn"};
205 long diameter[] = {4880, 12112, 12742, 6800, 143000, 121000};
206 float density[] = { 5.1f, 5.3f, 5.52f, 3.94f, 1.33f, 0.69f};
207
208 status=0;
209
210 /* open the FITS file containing a primary array and an ASCII table */
211 if ( fits_open_file(&fptr, filename, READWRITE, &status) )
212 printerror( status );
213
214 if ( fits_movabs_hdu(fptr, 2, &hdutype, &status) ) /* move to 2nd HDU */
215 printerror( status );
216
217 /* append a new empty binary table onto the FITS file */
218 if ( fits_create_tbl( fptr, BINARY_TBL, nrows, tfields, ttype, tform,
219 tunit, extname, &status) )
220 printerror( status );
221
222 firstrow = 1; /* first row in table to write */
223 firstelem = 1; /* first element in row (ignored in ASCII tables) */
224
225 /* write names to the first column (character strings) */
226 /* write diameters to the second column (longs) */
227 /* write density to the third column (floats) */
228
229 fits_write_col(fptr, TSTRING, 1, firstrow, firstelem, nrows, planet,
230 &status);
231 fits_write_col(fptr, TLONG, 2, firstrow, firstelem, nrows, diameter,
232 &status);
233 fits_write_col(fptr, TFLOAT, 3, firstrow, firstelem, nrows, density,
234 &status);
235
236 if ( fits_close_file(fptr, &status) ) /* close the FITS file */
237 printerror( status );
238 return;
239 }
240 /*--------------------------------------------------------------------------*/
copyhdu(void)241 void copyhdu( void)
242 {
243 /********************************************************************/
244 /* copy the 1st and 3rd HDUs from the input file to a new FITS file */
245 /********************************************************************/
246 fitsfile *infptr; /* pointer to the FITS file, defined in fitsio.h */
247 fitsfile *outfptr; /* pointer to the new FITS file */
248
249 char infilename[] = "atestfil.fit"; /* name for existing FITS file */
250 char outfilename[] = "btestfil.fit"; /* name for new FITS file */
251
252 int status, morekeys, hdutype;
253
254 status = 0;
255
256 remove(outfilename); /* Delete old file if it already exists */
257
258 /* open the existing FITS file */
259 if ( fits_open_file(&infptr, infilename, READONLY, &status) )
260 printerror( status );
261
262 if (fits_create_file(&outfptr, outfilename, &status)) /*create FITS file*/
263 printerror( status ); /* call printerror if error occurs */
264
265 /* copy the primary array from the input file to the output file */
266 morekeys = 0; /* don't reserve space for additional keywords */
267 if ( fits_copy_hdu(infptr, outfptr, morekeys, &status) )
268 printerror( status );
269
270 /* move to the 3rd HDU in the input file */
271 if ( fits_movabs_hdu(infptr, 3, &hdutype, &status) )
272 printerror( status );
273
274 /* copy 3rd HDU from the input file to the output file (to 2nd HDU) */
275 if ( fits_copy_hdu(infptr, outfptr, morekeys, &status) )
276 printerror( status );
277
278 if (fits_close_file(outfptr, &status) ||
279 fits_close_file(infptr, &status)) /* close files */
280 printerror( status );
281
282 return;
283 }
284 /*--------------------------------------------------------------------------*/
selectrows(void)285 void selectrows( void )
286
287 /*********************************************************************/
288 /* select rows from an input table and copy them to the output table */
289 /*********************************************************************/
290 {
291 fitsfile *infptr, *outfptr; /* pointer to input and output FITS files */
292 unsigned char *buffer;
293 char card[FLEN_CARD];
294 int status, hdutype, nkeys, keypos, nfound, colnum, anynulls, ii;
295 long naxes[2], frow, felem, noutrows, irow;
296 float nullval, density[6];
297
298 char infilename[] = "atestfil.fit"; /* name for existing FITS file */
299 char outfilename[] = "btestfil.fit"; /* name for new FITS file */
300
301 status = 0;
302
303 /* open the existing FITS files */
304 if ( fits_open_file(&infptr, infilename, READONLY, &status) ||
305 fits_open_file(&outfptr, outfilename, READWRITE, &status) )
306 printerror( status );
307
308 /* move to the 3rd HDU in the input file (a binary table in this case) */
309 if ( fits_movabs_hdu(infptr, 3, &hdutype, &status) )
310 printerror( status );
311
312 if (hdutype != BINARY_TBL) {
313 printf("Error: expected to find a binary table in this HDU\n");
314 return;
315 }
316 /* move to the last (2rd) HDU in the output file */
317 if ( fits_movabs_hdu(outfptr, 2, &hdutype, &status) )
318 printerror( status );
319
320 /* create new extension in the output file */
321 if ( fits_create_hdu(outfptr, &status) )
322 printerror( status );
323
324 /* get number of keywords */
325 if ( fits_get_hdrpos(infptr, &nkeys, &keypos, &status) )
326 printerror( status );
327
328 /* copy all the keywords from the input to the output extension */
329 for (ii = 1; ii <= nkeys; ii++) {
330 fits_read_record (infptr, ii, card, &status);
331 fits_write_record(outfptr, card, &status);
332 }
333
334 /* read the NAXIS1 and NAXIS2 keyword to get table size */
335 if (fits_read_keys_lng(infptr, "NAXIS", 1, 2, naxes, &nfound, &status) )
336 printerror( status );
337
338 /* find which column contains the DENSITY values */
339 if ( fits_get_colnum(infptr, CASEINSEN, "density", &colnum, &status) )
340 printerror( status );
341
342 /* read the DENSITY column values */
343 frow = 1;
344 felem = 1;
345 nullval = -99.;
346 if (fits_read_col(infptr, TFLOAT, colnum, frow, felem, naxes[1],
347 &nullval, density, &anynulls, &status) )
348 printerror( status );
349
350 /* allocate buffer large enough for 1 row of the table */
351 buffer = (unsigned char *) malloc(naxes[0]);
352
353 /* If the density is less than 3.0, copy the row to the output table */
354 for (noutrows = 0, irow = 1; irow <= naxes[1]; irow++) {
355 if (density[irow - 1] < 3.0) {
356 noutrows++;
357 fits_read_tblbytes( infptr, irow, 1, naxes[0], buffer, &status);
358 fits_write_tblbytes(outfptr, noutrows, 1, naxes[0], buffer, &status);
359 } }
360
361 /* update the NAXIS2 keyword with the correct number of rows */
362 if ( fits_update_key(outfptr, TLONG, "NAXIS2", &noutrows, 0, &status) )
363 printerror( status );
364
365 if (fits_close_file(outfptr, &status) || fits_close_file(infptr, &status))
366 printerror( status );
367
368 return;
369 }
370 /*--------------------------------------------------------------------------*/
readheader(void)371 void readheader ( void )
372
373 /**********************************************************************/
374 /* Print out all the header keywords in all extensions of a FITS file */
375 /**********************************************************************/
376 {
377 fitsfile *fptr; /* pointer to the FITS file, defined in fitsio.h */
378
379 int status, nkeys, keypos, hdutype, ii, jj;
380 char filename[] = "atestfil.fit"; /* name of existing FITS file */
381 char card[FLEN_CARD]; /* standard string lengths defined in fitsioc.h */
382
383 status = 0;
384
385 if ( fits_open_file(&fptr, filename, READONLY, &status) )
386 printerror( status );
387
388 /* attempt to move to next HDU, until we get an EOF error */
389 for (ii = 1; !(fits_movabs_hdu(fptr, ii, &hdutype, &status) ); ii++)
390 {
391 /* get no. of keywords */
392 if (fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
393 printerror( status );
394
395 printf("Header listing for HDU #%d:\n", ii);
396 for (jj = 1; jj <= nkeys; jj++) {
397 if ( fits_read_record(fptr, jj, card, &status) )
398 printerror( status );
399
400 printf("%s\n", card); /* print the keyword card */
401 }
402 printf("END\n\n"); /* terminate listing with END */
403 }
404
405 if (status == END_OF_FILE) /* status values are defined in fitsioc.h */
406 status = 0; /* got the expected EOF error; reset = 0 */
407 else
408 printerror( status ); /* got an unexpected error */
409
410 if ( fits_close_file(fptr, &status) )
411 printerror( status );
412
413 return;
414 }
415 /*--------------------------------------------------------------------------*/
readimage(void)416 void readimage( void )
417
418 /************************************************************************/
419 /* Read a FITS image and determine the minimum and maximum pixel values */
420 /************************************************************************/
421 {
422 fitsfile *fptr; /* pointer to the FITS file, defined in fitsio.h */
423 int status, nfound, anynull;
424 long naxes[2], fpixel, nbuffer, npixels, ii;
425
426 #define buffsize 1000
427 float datamin, datamax, nullval, buffer[buffsize];
428 char filename[] = "atestfil.fit"; /* name of existing FITS file */
429
430 status = 0;
431
432 if ( fits_open_file(&fptr, filename, READONLY, &status) )
433 printerror( status );
434
435 /* read the NAXIS1 and NAXIS2 keyword to get image size */
436 if ( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status) )
437 printerror( status );
438
439 npixels = naxes[0] * naxes[1]; /* number of pixels in the image */
440 fpixel = 1;
441 nullval = 0; /* don't check for null values in the image */
442 datamin = 1.0E30f;
443 datamax = -1.0E30f;
444
445 while (npixels > 0)
446 {
447 nbuffer = npixels;
448 if (npixels > buffsize)
449 nbuffer = buffsize; /* read as many pixels as will fit in buffer */
450
451 /* Note that even though the FITS images contains unsigned integer */
452 /* pixel values (or more accurately, signed integer pixels with */
453 /* a bias of 32768), this routine is reading the values into a */
454 /* float array. Cfitsio automatically performs the datatype */
455 /* conversion in cases like this. */
456
457 if ( fits_read_img(fptr, TFLOAT, fpixel, nbuffer, &nullval,
458 buffer, &anynull, &status) )
459 printerror( status );
460
461 for (ii = 0; ii < nbuffer; ii++) {
462 if ( buffer[ii] < datamin )
463 datamin = buffer[ii];
464
465 if ( buffer[ii] > datamax )
466 datamax = buffer[ii];
467 }
468 npixels -= nbuffer; /* increment remaining number of pixels */
469 fpixel += nbuffer; /* next pixel to be read in image */
470 }
471
472 printf("\nMin and max image pixels = %.0f, %.0f\n", datamin, datamax);
473
474 if ( fits_close_file(fptr, &status) )
475 printerror( status );
476
477 return;
478 }
479 /*--------------------------------------------------------------------------*/
readtable(void)480 void readtable( void )
481
482 /************************************************************/
483 /* read and print data values from an ASCII or binary table */
484 /************************************************************/
485 {
486 fitsfile *fptr; /* pointer to the FITS file, defined in fitsio.h */
487 int status, hdunum, hdutype, nfound, anynull, ii;
488 long frow, felem, nelem, longnull, dia[6];
489 float floatnull, den[6];
490 char strnull[10], *name[6], *ttype[3];
491
492 char filename[] = "atestfil.fit"; /* name of existing FITS file */
493
494 status = 0;
495
496 if ( fits_open_file(&fptr, filename, READONLY, &status) )
497 printerror( status );
498
499 for (ii = 0; ii < 3; ii++) /* allocate space for the column labels */
500 ttype[ii] = (char *) malloc(FLEN_VALUE); /* max label length = 69 */
501
502 for (ii = 0; ii < 6; ii++) /* allocate space for string column value */
503 name[ii] = (char *) malloc(10);
504
505 for (hdunum = 2; hdunum <= 3; hdunum++) /*read ASCII, then binary table */
506 {
507 /* move to the HDU */
508 if ( fits_movabs_hdu(fptr, hdunum, &hdutype, &status) )
509 printerror( status );
510
511 if (hdutype == ASCII_TBL)
512 printf("\nReading ASCII table in HDU %d:\n", hdunum);
513 else if (hdutype == BINARY_TBL)
514 printf("\nReading binary table in HDU %d:\n", hdunum);
515 else
516 {
517 printf("Error: this HDU is not an ASCII or binary table\n");
518 printerror( status );
519 }
520
521 /* read the column names from the TTYPEn keywords */
522 fits_read_keys_str(fptr, "TTYPE", 1, 3, ttype, &nfound, &status);
523
524 printf(" Row %10s %10s %10s\n", ttype[0], ttype[1], ttype[2]);
525
526 frow = 1;
527 felem = 1;
528 nelem = 6;
529 strcpy(strnull, " ");
530 longnull = 0;
531 floatnull = 0.;
532
533 /* read the columns */
534 fits_read_col(fptr, TSTRING, 1, frow, felem, nelem, strnull, name,
535 &anynull, &status);
536 fits_read_col(fptr, TLONG, 2, frow, felem, nelem, &longnull, dia,
537 &anynull, &status);
538 fits_read_col(fptr, TFLOAT, 3, frow, felem, nelem, &floatnull, den,
539 &anynull, &status);
540
541 for (ii = 0; ii < 6; ii++)
542 printf("%5d %10s %10ld %10.2f\n", ii + 1, name[ii], dia[ii], den[ii]);
543 }
544
545 for (ii = 0; ii < 3; ii++) /* free the memory for the column labels */
546 free( ttype[ii] );
547
548 for (ii = 0; ii < 6; ii++) /* free the memory for the string column */
549 free( name[ii] );
550
551 if ( fits_close_file(fptr, &status) )
552 printerror( status );
553
554 return;
555 }
556 /*--------------------------------------------------------------------------*/
printerror(int status)557 void printerror( int status)
558 {
559 /*****************************************************/
560 /* Print out cfitsio error messages and exit program */
561 /*****************************************************/
562
563
564 if (status)
565 {
566 fits_report_error(stderr, status); /* print error report */
567
568 exit( status ); /* terminate the program, returning error status */
569 }
570 return;
571 }
572