"Fossies" - the Fresh Open Source Software Archive

Member "cfitsio-4.0.0/cookbook.c" (20 May 2021, 21858 Bytes) of package /linux/misc/cfitsio-4.0.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "cookbook.c" see the Fossies "Dox" file reference documentation.

    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 
   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 /*--------------------------------------------------------------------------*/
   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 /*--------------------------------------------------------------------------*/
  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 /*--------------------------------------------------------------------------*/
  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 /*--------------------------------------------------------------------------*/
  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 /*--------------------------------------------------------------------------*/
  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 /*--------------------------------------------------------------------------*/
  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 /*--------------------------------------------------------------------------*/
  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 /*--------------------------------------------------------------------------*/
  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 /*--------------------------------------------------------------------------*/
  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 }