Initial commit

1 files changed, 571 insertions, 0 deletions

diff --git a/vendor/cfitsio/cookbook.c b/vendor/cfitsio/cookbook.c
new file mode 100644
index 00000000..3b42ac61
--- /dev/null
+++ b/vendor/cfitsio/cookbook.c
@@ -0,0 +1,571 @@
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/*
+  Every program which uses the CFITSIO interface must include the
+  the fitsio.h header file.  This contains the prototypes for all
+  the routines and defines the error status values and other symbolic
+  constants used in the interface.  
+*/
+#include "fitsio.h"
+
+int main( void );
+void writeimage( void );
+void writeascii( void );
+void writebintable( void );
+void copyhdu( void );
+void selectrows( void );
+void readheader( void );
+void readimage( void );
+void readtable( void );
+void printerror( int status);
+
+int main()
+{
+/*************************************************************************
+   This is a simple main program that calls the following routines:
+
+    writeimage    - write a FITS primary array image
+    writeascii    - write a FITS ASCII table extension
+    writebintable - write a FITS binary table extension
+    copyhdu       - copy a header/data unit from one FITS file to another
+    selectrows    - copy selected row from one HDU to another
+    readheader    - read and print the header keywords in every extension
+    readimage     - read a FITS image and compute the min and max value
+    readtable     - read columns of data from ASCII and binary tables
+
+**************************************************************************/
+
+    writeimage();
+    writeascii();
+    writebintable();
+    copyhdu();
+    selectrows();
+    readheader();
+    readimage();
+    readtable();
+
+    printf("\nAll the cfitsio cookbook routines ran successfully.\n");
+    return(0);
+}
+/*--------------------------------------------------------------------------*/
+void writeimage( void )
+
+    /******************************************************/
+    /* Create a FITS primary array containing a 2-D image */
+    /******************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status, ii, jj;
+    long  fpixel, nelements, exposure;
+    unsigned short *array[200];
+
+    /* initialize FITS image parameters */
+    char filename[] = "atestfil.fit";             /* name for new FITS file */
+    int bitpix   =  USHORT_IMG; /* 16-bit unsigned short pixel values       */
+    long naxis    =   2;  /* 2-dimensional image                            */    
+    long naxes[2] = { 300, 200 };   /* image is 300 pixels wide by 200 rows */
+
+    /* allocate memory for the whole image */ 
+    array[0] = (unsigned short *)malloc( naxes[0] * naxes[1]
+                                        * sizeof( unsigned short ) );
+
+    /* initialize pointers to the start of each row of the image */
+    for( ii=1; ii<naxes[1]; ii++ )
+      array[ii] = array[ii-1] + naxes[0];
+
+    remove(filename);               /* Delete old file if it already exists */
+
+    status = 0;         /* initialize status before calling fitsio routines */
+
+    if (fits_create_file(&fptr, filename, &status)) /* create new FITS file */
+         printerror( status );           /* call printerror if error occurs */
+
+    /* write the required keywords for the primary array image.     */
+    /* Since bitpix = USHORT_IMG, this will cause cfitsio to create */
+    /* a FITS image with BITPIX = 16 (signed short integers) with   */
+    /* BSCALE = 1.0 and BZERO = 32768.  This is the convention that */
+    /* FITS uses to store unsigned integers.  Note that the BSCALE  */
+    /* and BZERO keywords will be automatically written by cfitsio  */
+    /* in this case.                                                */
+
+    if ( fits_create_img(fptr,  bitpix, naxis, naxes, &status) )
+         printerror( status );          
+
+    /* initialize the values in the image with a linear ramp function */
+    for (jj = 0; jj < naxes[1]; jj++)
+    {   for (ii = 0; ii < naxes[0]; ii++)
+        {
+            array[jj][ii] = ii + jj;
+        }
+    }
+
+    fpixel = 1;                               /* first pixel to write      */
+    nelements = naxes[0] * naxes[1];          /* number of pixels to write */
+
+    /* write the array of unsigned integers to the FITS file */
+    if ( fits_write_img(fptr, TUSHORT, fpixel, nelements, array[0], &status) )
+        printerror( status );
+      
+    free( array[0] );  /* free previously allocated memory */
+
+    /* write another optional keyword to the header */
+    /* Note that the ADDRESS of the value is passed in the routine */
+    exposure = 1500.;
+    if ( fits_update_key(fptr, TLONG, "EXPOSURE", &exposure,
+         "Total Exposure Time", &status) )
+         printerror( status );           
+
+    if ( fits_close_file(fptr, &status) )                /* close the file */
+         printerror( status );           
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void writeascii ( void )
+
+    /*******************************************************************/
+    /* Create an ASCII table extension containing 3 columns and 6 rows */
+    /*******************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status;
+    long firstrow, firstelem;
+
+    int tfields = 3;       /* table will have 3 columns */
+    long nrows  = 6;       /* table will have 6 rows    */
+
+    char filename[] = "atestfil.fit";           /* name for new FITS file */
+    char extname[] = "PLANETS_ASCII";             /* extension name */
+
+    /* define the name, datatype, and physical units for the 3 columns */
+    char *ttype[] = { "Planet", "Diameter", "Density" };
+    char *tform[] = { "a8",     "I6",       "F4.2"    };
+    char *tunit[] = { "\0",      "km",       "g/cm"    };
+
+    /* define the name diameter, and density of each planet */
+    char *planet[] = {"Mercury", "Venus", "Earth", "Mars","Jupiter","Saturn"};
+    long  diameter[] = {4880,    12112,    12742,   6800,  143000,   121000};
+    float density[]  = { 5.1,     5.3,      5.52,   3.94,    1.33,    0.69};
+
+    status=0;
+
+    /* open with write access the FITS file containing a primary array */
+    if ( fits_open_file(&fptr, filename, READWRITE, &status) ) 
+         printerror( status );
+
+    /* append a new empty ASCII table onto the FITS file */
+    if ( fits_create_tbl( fptr, ASCII_TBL, nrows, tfields, ttype, tform,
+                tunit, extname, &status) )
+         printerror( status );
+
+    firstrow  = 1;  /* first row in table to write   */
+    firstelem = 1;  /* first element in row  (ignored in ASCII tables) */
+
+    /* write names to the first column (character strings) */
+    /* write diameters to the second column (longs) */
+    /* write density to the third column (floats) */
+
+    fits_write_col(fptr, TSTRING, 1, firstrow, firstelem, nrows, planet,
+                   &status);
+    fits_write_col(fptr, TLONG, 2, firstrow, firstelem, nrows, diameter,
+                   &status);
+    fits_write_col(fptr, TFLOAT, 3, firstrow, firstelem, nrows, density,
+                   &status);
+
+    if ( fits_close_file(fptr, &status) )       /* close the FITS file */
+         printerror( status );
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void writebintable ( void )
+
+    /*******************************************************************/
+    /* Create a binary table extension containing 3 columns and 6 rows */
+    /*******************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status, hdutype;
+    long firstrow, firstelem;
+
+    int tfields   = 3;       /* table will have 3 columns */
+    long nrows    = 6;       /* table will have 6 rows    */
+
+    char filename[] = "atestfil.fit";           /* name for new FITS file */
+    char extname[] = "PLANETS_Binary";           /* extension name */
+
+    /* define the name, datatype, and physical units for the 3 columns */
+    char *ttype[] = { "Planet", "Diameter", "Density" };
+    char *tform[] = { "8a",     "1J",       "1E"    };
+    char *tunit[] = { "\0",      "km",       "g/cm"    };
+
+    /* define the name diameter, and density of each planet */
+    char *planet[] = {"Mercury", "Venus", "Earth", "Mars","Jupiter","Saturn"};
+    long  diameter[] = {4880,     12112,   12742,   6800,  143000,   121000};
+    float density[]  = { 5.1,      5.3,     5.52,   3.94,   1.33,     0.69};
+
+    status=0;
+
+    /* open the FITS file containing a primary array and an ASCII table */
+    if ( fits_open_file(&fptr, filename, READWRITE, &status) ) 
+         printerror( status );
+
+    if ( fits_movabs_hdu(fptr, 2, &hdutype, &status) ) /* move to 2nd HDU */
+         printerror( status );
+
+    /* append a new empty binary table onto the FITS file */
+    if ( fits_create_tbl( fptr, BINARY_TBL, nrows, tfields, ttype, tform,
+                tunit, extname, &status) )
+         printerror( status );
+
+    firstrow  = 1;  /* first row in table to write   */
+    firstelem = 1;  /* first element in row  (ignored in ASCII tables) */
+
+    /* write names to the first column (character strings) */
+    /* write diameters to the second column (longs) */
+    /* write density to the third column (floats) */
+
+    fits_write_col(fptr, TSTRING, 1, firstrow, firstelem, nrows, planet,
+                   &status);
+    fits_write_col(fptr, TLONG, 2, firstrow, firstelem, nrows, diameter,
+                   &status);
+    fits_write_col(fptr, TFLOAT, 3, firstrow, firstelem, nrows, density,
+                   &status);
+
+    if ( fits_close_file(fptr, &status) )       /* close the FITS file */
+         printerror( status );
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void copyhdu( void)
+{
+    /********************************************************************/
+    /* copy the 1st and 3rd HDUs from the input file to a new FITS file */
+    /********************************************************************/
+    fitsfile *infptr;      /* pointer to the FITS file, defined in fitsio.h */
+    fitsfile *outfptr;                 /* pointer to the new FITS file      */
+
+    char infilename[]  = "atestfil.fit";  /* name for existing FITS file   */
+    char outfilename[] = "btestfil.fit";  /* name for new FITS file        */
+
+    int status, morekeys, hdutype;
+
+    status = 0;
+
+    remove(outfilename);            /* Delete old file if it already exists */
+
+    /* open the existing FITS file */
+    if ( fits_open_file(&infptr, infilename, READONLY, &status) ) 
+         printerror( status );
+
+    if (fits_create_file(&outfptr, outfilename, &status)) /*create FITS file*/
+         printerror( status );           /* call printerror if error occurs */
+
+    /* copy the primary array from the input file to the output file */
+    morekeys = 0;     /* don't reserve space for additional keywords */
+    if ( fits_copy_hdu(infptr, outfptr, morekeys, &status) )
+         printerror( status );
+
+    /* move to the 3rd HDU in the input file */
+    if ( fits_movabs_hdu(infptr, 3, &hdutype, &status) )
+         printerror( status );
+
+    /* copy 3rd HDU from the input file to the output file (to 2nd HDU) */
+    if ( fits_copy_hdu(infptr, outfptr, morekeys, &status) )
+         printerror( status );
+
+    if (fits_close_file(outfptr, &status) ||
+        fits_close_file(infptr, &status)) /* close files */
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void selectrows( void )
+
+    /*********************************************************************/
+    /* select rows from an input table and copy them to the output table */
+    /*********************************************************************/
+{
+    fitsfile *infptr, *outfptr;  /* pointer to input and output FITS files */
+    unsigned char *buffer;
+    char card[FLEN_CARD];
+    int status, hdutype, nkeys, keypos, nfound, colnum, anynulls, ii;
+    long naxes[2], frow, felem, noutrows, irow;
+    float nullval, density[6];
+
+    char infilename[]  = "atestfil.fit";  /* name for existing FITS file   */
+    char outfilename[] = "btestfil.fit";  /* name for new FITS file        */
+
+    status = 0;
+
+    /* open the existing FITS files */
+    if ( fits_open_file(&infptr,  infilename,  READONLY,  &status) ||
+         fits_open_file(&outfptr, outfilename, READWRITE, &status) ) 
+         printerror( status );
+
+    /* move to the 3rd HDU in the input file (a binary table in this case) */
+    if ( fits_movabs_hdu(infptr, 3, &hdutype, &status) )
+         printerror( status );
+
+    if (hdutype != BINARY_TBL)  {
+        printf("Error: expected to find a binary table in this HDU\n");
+        return;
+    }
+    /* move to the last (2rd) HDU in the output file */
+    if ( fits_movabs_hdu(outfptr, 2, &hdutype, &status) )
+         printerror( status );
+
+    /* create new extension in the output file */
+    if ( fits_create_hdu(outfptr, &status) )
+         printerror( status );
+
+    /* get number of keywords */
+    if ( fits_get_hdrpos(infptr, &nkeys, &keypos, &status) ) 
+         printerror( status );
+
+    /* copy all the keywords from the input to the output extension */
+    for (ii = 1; ii <= nkeys; ii++)  {
+        fits_read_record (infptr, ii, card, &status); 
+        fits_write_record(outfptr,    card, &status); 
+    }
+
+    /* read the NAXIS1 and NAXIS2 keyword to get table size */
+    if (fits_read_keys_lng(infptr, "NAXIS", 1, 2, naxes, &nfound, &status) )
+         printerror( status );
+
+    /* find which column contains the DENSITY values */
+    if ( fits_get_colnum(infptr, CASEINSEN, "density", &colnum, &status) )
+         printerror( status );
+
+    /* read the DENSITY column values */
+    frow = 1;
+    felem = 1;
+    nullval = -99.;
+    if (fits_read_col(infptr, TFLOAT, colnum, frow, felem, naxes[1], 
+        &nullval, density, &anynulls, &status) )
+        printerror( status );
+
+    /* allocate buffer large enough for 1 row of the table */
+    buffer = (unsigned char *) malloc(naxes[0]);
+
+    /* If the density is less than 3.0, copy the row to the output table */
+    for (noutrows = 0, irow = 1; irow <= naxes[1]; irow++)  {
+      if (density[irow - 1] < 3.0)  {
+        noutrows++;
+        fits_read_tblbytes( infptr, irow,      1, naxes[0], buffer, &status); 
+        fits_write_tblbytes(outfptr, noutrows, 1, naxes[0], buffer, &status); 
+    } }
+
+    /* update the NAXIS2 keyword with the correct number of rows */
+    if ( fits_update_key(outfptr, TLONG, "NAXIS2", &noutrows, 0, &status) )
+         printerror( status );
+
+    if (fits_close_file(outfptr, &status) || fits_close_file(infptr, &status))
+        printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void readheader ( void )
+
+    /**********************************************************************/
+    /* Print out all the header keywords in all extensions of a FITS file */
+    /**********************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+
+    int status, nkeys, keypos, hdutype, ii, jj;
+    char filename[]  = "atestfil.fit";     /* name of existing FITS file   */
+    char card[FLEN_CARD];   /* standard string lengths defined in fitsioc.h */
+
+    status = 0;
+
+    if ( fits_open_file(&fptr, filename, READONLY, &status) ) 
+         printerror( status );
+
+    /* attempt to move to next HDU, until we get an EOF error */
+    for (ii = 1; !(fits_movabs_hdu(fptr, ii, &hdutype, &status) ); ii++) 
+    {
+        /* get no. of keywords */
+        if (fits_get_hdrpos(fptr, &nkeys, &keypos, &status) )
+            printerror( status );
+
+        printf("Header listing for HDU #%d:\n", ii);
+        for (jj = 1; jj <= nkeys; jj++)  {
+            if ( fits_read_record(fptr, jj, card, &status) )
+                 printerror( status );
+
+            printf("%s\n", card); /* print the keyword card */
+        }
+        printf("END\n\n");  /* terminate listing with END */
+    }
+
+    if (status == END_OF_FILE)   /* status values are defined in fitsioc.h */
+        status = 0;              /* got the expected EOF error; reset = 0  */
+    else
+       printerror( status );     /* got an unexpected error                */
+
+    if ( fits_close_file(fptr, &status) )
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void readimage( void )
+
+    /************************************************************************/
+    /* Read a FITS image and determine the minimum and maximum pixel values */
+    /************************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status,  nfound, anynull;
+    long naxes[2], fpixel, nbuffer, npixels, ii;
+
+#define buffsize 1000
+    float datamin, datamax, nullval, buffer[buffsize];
+    char filename[]  = "atestfil.fit";     /* name of existing FITS file   */
+
+    status = 0;
+
+    if ( fits_open_file(&fptr, filename, READONLY, &status) )
+         printerror( status );
+
+    /* read the NAXIS1 and NAXIS2 keyword to get image size */
+    if ( fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status) )
+         printerror( status );
+
+    npixels  = naxes[0] * naxes[1];         /* number of pixels in the image */
+    fpixel   = 1;
+    nullval  = 0;                /* don't check for null values in the image */
+    datamin  = 1.0E30;
+    datamax  = -1.0E30;
+
+    while (npixels > 0)
+    {
+      nbuffer = npixels;
+      if (npixels > buffsize)
+        nbuffer = buffsize;     /* read as many pixels as will fit in buffer */
+
+      /* Note that even though the FITS images contains unsigned integer */
+      /* pixel values (or more accurately, signed integer pixels with    */
+      /* a bias of 32768),  this routine is reading the values into a    */
+      /* float array.   Cfitsio automatically performs the datatype      */
+      /* conversion in cases like this.                                  */
+
+      if ( fits_read_img(fptr, TFLOAT, fpixel, nbuffer, &nullval,
+                  buffer, &anynull, &status) )
+           printerror( status );
+
+      for (ii = 0; ii < nbuffer; ii++)  {
+        if ( buffer[ii] < datamin )
+            datamin = buffer[ii];
+
+        if ( buffer[ii] > datamax )
+            datamax = buffer[ii];
+      }
+      npixels -= nbuffer;    /* increment remaining number of pixels */
+      fpixel  += nbuffer;    /* next pixel to be read in image */
+    }
+
+    printf("\nMin and max image pixels =  %.0f, %.0f\n", datamin, datamax);
+
+    if ( fits_close_file(fptr, &status) )
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void readtable( void )
+
+    /************************************************************/
+    /* read and print data values from an ASCII or binary table */
+    /************************************************************/
+{
+    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
+    int status, hdunum, hdutype,  nfound, anynull, ii;
+    long frow, felem, nelem, longnull, dia[6];
+    float floatnull, den[6];
+    char strnull[10], *name[6], *ttype[3]; 
+
+    char filename[]  = "atestfil.fit";     /* name of existing FITS file   */
+
+    status = 0;
+
+    if ( fits_open_file(&fptr, filename, READONLY, &status) )
+         printerror( status );
+
+    for (ii = 0; ii < 3; ii++)      /* allocate space for the column labels */
+        ttype[ii] = (char *) malloc(FLEN_VALUE);  /* max label length = 69 */
+
+    for (ii = 0; ii < 6; ii++)    /* allocate space for string column value */
+        name[ii] = (char *) malloc(10);   
+
+    for (hdunum = 2; hdunum <= 3; hdunum++) /*read ASCII, then binary table */
+    {
+      /* move to the HDU */
+      if ( fits_movabs_hdu(fptr, hdunum, &hdutype, &status) ) 
+           printerror( status );
+
+      if (hdutype == ASCII_TBL)
+          printf("\nReading ASCII table in HDU %d:\n",  hdunum);
+      else if (hdutype == BINARY_TBL)
+          printf("\nReading binary table in HDU %d:\n", hdunum);
+      else
+      {
+          printf("Error: this HDU is not an ASCII or binary table\n");
+          printerror( status );
+      }
+
+      /* read the column names from the TTYPEn keywords */
+      fits_read_keys_str(fptr, "TTYPE", 1, 3, ttype, &nfound, &status);
+
+      printf(" Row  %10s %10s %10s\n", ttype[0], ttype[1], ttype[2]);
+
+      frow      = 1;
+      felem     = 1;
+      nelem     = 6;
+      strcpy(strnull, " ");
+      longnull  = 0;
+      floatnull = 0.;
+
+      /*  read the columns */  
+      fits_read_col(fptr, TSTRING, 1, frow, felem, nelem, strnull,  name,
+                    &anynull, &status);
+      fits_read_col(fptr, TLONG, 2, frow, felem, nelem, &longnull,  dia,
+                    &anynull, &status);
+      fits_read_col(fptr, TFLOAT, 3, frow, felem, nelem, &floatnull, den,
+                    &anynull, &status);
+
+      for (ii = 0; ii < 6; ii++)
+        printf("%5d %10s %10ld %10.2f\n", ii + 1, name[ii], dia[ii], den[ii]);
+    }
+
+    for (ii = 0; ii < 3; ii++)      /* free the memory for the column labels */
+        free( ttype[ii] );
+
+    for (ii = 0; ii < 6; ii++)      /* free the memory for the string column */
+        free( name[ii] );
+
+    if ( fits_close_file(fptr, &status) ) 
+         printerror( status );
+
+    return;
+}
+/*--------------------------------------------------------------------------*/
+void printerror( int status)
+{
+    /*****************************************************/
+    /* Print out cfitsio error messages and exit program */
+    /*****************************************************/
+
+
+    if (status)
+    {
+       fits_report_error(stderr, status); /* print error report */
+
+       exit( status );    /* terminate the program, returning error status */
+    }
+    return;
+}