Repatch (from linux) of OSX IRAF

1 files changed, 2221 insertions, 0 deletions

diff --git a/vendor/cfitsio/histo.c b/vendor/cfitsio/histo.c
new file mode 100644
index 00000000..6af856f8
--- /dev/null
+++ b/vendor/cfitsio/histo.c
@@ -0,0 +1,2221 @@
+/*   Globally defined histogram parameters */
+#include <string.h>
+#include <ctype.h>
+#include <math.h>
+#include <stdlib.h>
+#include "fitsio2.h"
+
+typedef struct {  /*  Structure holding all the histogramming information   */
+   union {        /*  the iterator work functions (ffwritehist, ffcalchist) */
+      char   *b;  /*  need to do their job... passed via *userPointer.      */
+      short  *i;
+      int    *j;
+      float  *r;
+      double *d;
+   } hist;
+
+   fitsfile *tblptr;
+
+   int   haxis, hcolnum[4], himagetype;
+   long  haxis1, haxis2, haxis3, haxis4;
+   float amin1, amin2, amin3, amin4;
+   float maxbin1, maxbin2, maxbin3, maxbin4;
+   float binsize1, binsize2, binsize3, binsize4;
+   int   wtrecip, wtcolnum;
+   float weight;
+   char  *rowselector;
+
+} histType;
+
+/*--------------------------------------------------------------------------*/
+int ffbins(char *binspec,   /* I - binning specification */
+                   int *imagetype,      /* O - image type, TINT or TSHORT */
+                   int *histaxis,       /* O - no. of axes in the histogram */
+                   char colname[4][FLEN_VALUE],  /* column name for axis */
+                   double *minin,        /* minimum value for each axis */
+                   double *maxin,        /* maximum value for each axis */
+                   double *binsizein,    /* size of bins on each axis */
+                   char minname[4][FLEN_VALUE],  /* keyword name for min */
+                   char maxname[4][FLEN_VALUE],  /* keyword name for max */
+                   char binname[4][FLEN_VALUE],  /* keyword name for binsize */
+                   double *wt,          /* weighting factor          */
+                   char *wtname,        /* keyword or column name for weight */
+                   int *recip,          /* the reciprocal of the weight? */
+                   int *status)
+{
+/*
+   Parse the input binning specification string, returning the binning
+   parameters.  Supports up to 4 dimensions.  The binspec string has
+   one of these forms:
+
+   bin binsize                  - 2D histogram with binsize on each axis
+   bin xcol                     - 1D histogram on column xcol
+   bin (xcol, ycol) = binsize   - 2D histogram with binsize on each axis
+   bin x=min:max:size, y=min:max:size, z..., t... 
+   bin x=:max, y=::size
+   bin x=size, y=min::size
+
+   most other reasonable combinations are supported.        
+*/
+    int ii, slen, defaulttype;
+    char *ptr, tmpname[30], *file_expr = NULL;
+    double  dummy;
+
+    if (*status > 0)
+         return(*status);
+
+    /* set the default values */
+    *histaxis = 2;
+    *imagetype = TINT;
+    defaulttype = 1;
+    *wt = 1.;
+    *recip = 0;
+    *wtname = '\0';
+
+    /* set default values */
+    for (ii = 0; ii < 4; ii++)
+    {
+        *colname[ii] = '\0';
+        *minname[ii] = '\0';
+        *maxname[ii] = '\0';
+        *binname[ii] = '\0';
+        minin[ii] = DOUBLENULLVALUE;  /* undefined values */
+        maxin[ii] = DOUBLENULLVALUE;
+        binsizein[ii] = DOUBLENULLVALUE;
+    }
+
+    ptr = binspec + 3;  /* skip over 'bin' */
+
+    if (*ptr == 'i' )  /* bini */
+    {
+        *imagetype = TSHORT;
+        defaulttype = 0;
+        ptr++;
+    }
+    else if (*ptr == 'j' )  /* binj; same as default */
+    {
+        defaulttype = 0;
+        ptr ++;
+    }
+    else if (*ptr == 'r' )  /* binr */
+    {
+        *imagetype = TFLOAT;
+        defaulttype = 0;
+        ptr ++;
+    }
+    else if (*ptr == 'd' )  /* bind */
+    {
+        *imagetype = TDOUBLE;
+        defaulttype = 0;
+        ptr ++;
+    }
+    else if (*ptr == 'b' )  /* binb */
+    {
+        *imagetype = TBYTE;
+        defaulttype = 0;
+        ptr ++;
+    }
+
+    if (*ptr == '\0')  /* use all defaults for other parameters */
+        return(*status);
+    else if (*ptr != ' ')  /* must be at least one blank */
+    {
+        ffpmsg("binning specification syntax error:");
+        ffpmsg(binspec);
+        return(*status = URL_PARSE_ERROR);
+    }
+
+    while (*ptr == ' ')  /* skip over blanks */
+           ptr++;
+
+    if (*ptr == '\0')   /* no other parameters; use defaults */
+        return(*status);
+
+    /* Check if need to import expression from a file */
+
+    if( *ptr=='@' ) {
+       if( ffimport_file( ptr+1, &file_expr, status ) ) return(*status);
+       ptr = file_expr;
+       while (*ptr == ' ')
+               ptr++;       /* skip leading white space... again */
+    }
+
+    if (*ptr == '(' )
+    {
+        /* this must be the opening parenthesis around a list of column */
+        /* names, optionally followed by a '=' and the binning spec. */
+
+        for (ii = 0; ii < 4; ii++)
+        {
+            ptr++;               /* skip over the '(', ',', or ' ') */
+            while (*ptr == ' ')  /* skip over blanks */
+                ptr++;
+
+            slen = strcspn(ptr, " ,)");
+            strncat(colname[ii], ptr, slen); /* copy 1st column name */
+
+            ptr += slen;
+            while (*ptr == ' ')  /* skip over blanks */
+                ptr++;
+
+            if (*ptr == ')' )   /* end of the list of names */
+            {
+                *histaxis = ii + 1;
+                break;
+            }
+        }
+
+        if (ii == 4)   /* too many names in the list , or missing ')'  */
+        {
+            ffpmsg(
+ "binning specification has too many column names or is missing closing ')':");
+            ffpmsg(binspec);
+	    if( file_expr ) free( file_expr );
+            return(*status = URL_PARSE_ERROR);
+        }
+
+        ptr++;  /* skip over the closing parenthesis */
+        while (*ptr == ' ')  /* skip over blanks */
+            ptr++;
+
+        if (*ptr == '\0') {
+	    if( file_expr ) free( file_expr );
+            return(*status);  /* parsed the entire string */
+	}
+
+        else if (*ptr != '=')  /* must be an equals sign now*/
+        {
+            ffpmsg("illegal binning specification in URL:");
+            ffpmsg(" an equals sign '=' must follow the column names");
+            ffpmsg(binspec);
+	    if( file_expr ) free( file_expr );
+            return(*status = URL_PARSE_ERROR);
+        }
+
+        ptr++;  /* skip over the equals sign */
+        while (*ptr == ' ')  /* skip over blanks */
+            ptr++;
+
+        /* get the single range specification for all the columns */
+        ffbinr(&ptr, tmpname, minin,
+                                     maxin, binsizein, minname[0],
+                                     maxname[0], binname[0], status);
+        if (*status > 0)
+        {
+            ffpmsg("illegal binning specification in URL:");
+            ffpmsg(binspec);
+	    if( file_expr ) free( file_expr );
+            return(*status);
+        }
+
+        for (ii = 1; ii < *histaxis; ii++)
+        {
+            minin[ii] = minin[0];
+            maxin[ii] = maxin[0];
+            binsizein[ii] = binsizein[0];
+            strcpy(minname[ii], minname[0]);
+            strcpy(maxname[ii], maxname[0]);
+            strcpy(binname[ii], binname[0]);
+        }
+
+        while (*ptr == ' ')  /* skip over blanks */
+            ptr++;
+
+        if (*ptr == ';')
+            goto getweight;   /* a weighting factor is specified */
+
+        if (*ptr != '\0')  /* must have reached end of string */
+        {
+            ffpmsg("illegal syntax after binning range specification in URL:");
+            ffpmsg(binspec);
+	    if( file_expr ) free( file_expr );
+            return(*status = URL_PARSE_ERROR);
+        }
+
+        return(*status);
+    }             /* end of case with list of column names in ( )  */
+
+    /* if we've reached this point, then the binning specification */
+    /* must be of the form: XCOL = min:max:binsize, YCOL = ...     */
+    /* where the column name followed by '=' are optional.         */
+    /* If the column name is not specified, then use the default name */
+
+    for (ii = 0; ii < 4; ii++) /* allow up to 4 histogram dimensions */
+    {
+        ffbinr(&ptr, colname[ii], &minin[ii],
+                                     &maxin[ii], &binsizein[ii], minname[ii],
+                                     maxname[ii], binname[ii], status);
+
+        if (*status > 0)
+        {
+            ffpmsg("illegal syntax in binning range specification in URL:");
+            ffpmsg(binspec);
+	    if( file_expr ) free( file_expr );
+            return(*status);
+        }
+
+        if (*ptr == '\0' || *ptr == ';')
+            break;        /* reached the end of the string */
+
+        if (*ptr == ' ')
+        {
+            while (*ptr == ' ')  /* skip over blanks */
+                ptr++;
+
+            if (*ptr == '\0' || *ptr == ';')
+                break;        /* reached the end of the string */
+
+            if (*ptr == ',')
+                ptr++;  /* comma separates the next column specification */
+        }
+        else if (*ptr == ',')
+        {          
+            ptr++;  /* comma separates the next column specification */
+        }
+        else
+        {
+            ffpmsg("illegal characters following binning specification in URL:");
+            ffpmsg(binspec);
+	    if( file_expr ) free( file_expr );
+            return(*status = URL_PARSE_ERROR);
+        }
+    }
+
+    if (ii == 4)
+    {
+        /* there are yet more characters in the string */
+        ffpmsg("illegal binning specification in URL:");
+        ffpmsg("apparently greater than 4 histogram dimensions");
+        ffpmsg(binspec);
+        return(*status = URL_PARSE_ERROR);
+    }
+    else
+        *histaxis = ii + 1;
+
+    /* special case: if a single number was entered it should be      */
+    /* interpreted as the binning factor for the default X and Y axes */
+
+    if (*histaxis == 1 && *colname[0] == '\0' && 
+         minin[0] == DOUBLENULLVALUE && maxin[0] == DOUBLENULLVALUE)
+    {
+        *histaxis = 2;
+        binsizein[1] = binsizein[0];
+    }
+
+getweight:
+    if (*ptr == ';')  /* looks like a weighting factor is given */
+    {
+        ptr++;
+       
+        while (*ptr == ' ')  /* skip over blanks */
+            ptr++;
+
+        recip = 0;
+        if (*ptr == '/')
+        {
+            *recip = 1;  /* the reciprocal of the weight is entered */
+            ptr++;
+
+            while (*ptr == ' ')  /* skip over blanks */
+                ptr++;
+        }
+
+        /* parse the weight as though it were a binrange. */
+        /* either a column name or a numerical value will be returned */
+
+        ffbinr(&ptr, wtname, &dummy, &dummy, wt, tmpname,
+                                     tmpname, tmpname, status);
+
+        if (*status > 0)
+        {
+            ffpmsg("illegal binning weight specification in URL:");
+            ffpmsg(binspec);
+	    if( file_expr ) free( file_expr );
+            return(*status);
+        }
+
+        /* creat a float datatype histogram by default, if weight */
+        /* factor is not = 1.0  */
+
+        if ( (defaulttype && *wt != 1.0) || (defaulttype && *wtname) )
+            *imagetype = TFLOAT;
+    }
+
+    while (*ptr == ' ')  /* skip over blanks */
+         ptr++;
+
+    if (*ptr != '\0')  /* should have reached the end of string */
+    {
+        ffpmsg("illegal syntax after binning weight specification in URL:");
+        ffpmsg(binspec);
+        *status = URL_PARSE_ERROR;
+    }
+
+    if( file_expr ) free( file_expr );
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int ffbinr(char **ptr, 
+                   char *colname, 
+                   double *minin,
+                   double *maxin, 
+                   double *binsizein,
+                   char *minname,
+                   char *maxname,
+                   char *binname,
+                   int *status)
+/*
+   Parse the input binning range specification string, returning 
+   the column name, histogram min and max values, and bin size.
+*/
+{
+    int slen, isanumber;
+    char token[FLEN_VALUE];
+
+    if (*status > 0)
+        return(*status);
+
+    slen = fits_get_token(ptr, " ,=:;", token, &isanumber); /* get 1st token */
+
+    if (slen == 0 && (**ptr == '\0' || **ptr == ',' || **ptr == ';') )
+        return(*status);   /* a null range string */
+
+    if (!isanumber && **ptr != ':')
+    {
+        /* this looks like the column name */
+
+        if (token[0] == '#' && isdigit((int) token[1]) )
+        {
+            /* omit the leading '#' in the column number */
+            strcpy(colname, token+1);
+        }
+        else
+            strcpy(colname, token);
+
+        while (**ptr == ' ')  /* skip over blanks */
+             (*ptr)++;
+
+        if (**ptr != '=')
+            return(*status);  /* reached the end */
+
+        (*ptr)++;   /* skip over the = sign */
+
+        while (**ptr == ' ')  /* skip over blanks */
+             (*ptr)++;
+
+        slen = fits_get_token(ptr, " ,:;", token, &isanumber); /* get token */
+    }
+
+    if (**ptr != ':')
+    {
+        /* this is the first token, and since it is not followed by */
+        /* a ':' this must be the binsize token */
+        if (!isanumber)
+            strcpy(binname, token);
+        else
+            *binsizein =  strtod(token, NULL);
+
+        return(*status);  /* reached the end */
+    }
+    else
+    {
+        /* the token contains the min value */
+        if (slen)
+        {
+            if (!isanumber)
+                strcpy(minname, token);
+            else
+                *minin = strtod(token, NULL);
+        }
+    }
+
+    (*ptr)++;  /* skip the colon between the min and max values */
+    slen = fits_get_token(ptr, " ,:;", token, &isanumber); /* get token */
+
+    /* the token contains the max value */
+    if (slen)
+    {
+        if (!isanumber)
+            strcpy(maxname, token);
+        else
+            *maxin = strtod(token, NULL);
+    }
+
+    if (**ptr != ':')
+        return(*status);  /* reached the end; no binsize token */
+
+    (*ptr)++;  /* skip the colon between the max and binsize values */
+    slen = fits_get_token(ptr, " ,:;", token, &isanumber); /* get token */
+
+    /* the token contains the binsize value */
+    if (slen)
+    {
+        if (!isanumber)
+            strcpy(binname, token);
+        else
+            *binsizein = strtod(token, NULL);
+    }
+
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int ffhist2(fitsfile **fptr,  /* IO - pointer to table with X and Y cols;    */
+                             /*     on output, points to histogram image    */
+           char *outfile,    /* I - name for the output histogram file      */
+           int imagetype,    /* I - datatype for image: TINT, TSHORT, etc   */
+           int naxis,        /* I - number of axes in the histogram image   */
+           char colname[4][FLEN_VALUE],   /* I - column names               */
+           double *minin,     /* I - minimum histogram value, for each axis */
+           double *maxin,     /* I - maximum histogram value, for each axis */
+           double *binsizein, /* I - bin size along each axis               */
+           char minname[4][FLEN_VALUE], /* I - optional keywords for min    */
+           char maxname[4][FLEN_VALUE], /* I - optional keywords for max    */
+           char binname[4][FLEN_VALUE], /* I - optional keywords for binsize */
+           double weightin,        /* I - binning weighting factor          */
+           char wtcol[FLEN_VALUE], /* I - optional keyword or col for weight*/
+           int recip,              /* I - use reciprocal of the weight?     */
+           char *selectrow,        /* I - optional array (length = no. of   */
+                             /* rows in the table).  If the element is true */
+                             /* then the corresponding row of the table will*/
+                             /* be included in the histogram, otherwise the */
+                             /* row will be skipped.  Ingnored if *selectrow*/
+                             /* is equal to NULL.                           */
+           int *status)
+{
+    fitsfile *histptr;
+    int   bitpix, colnum[4], wtcolnum;
+    long haxes[4];
+    float amin[4], amax[4], binsize[4],  weight;
+
+    if (*status > 0)
+        return(*status);
+
+    if (naxis > 4)
+    {
+        ffpmsg("histogram has more than 4 dimensions");
+        return(*status = BAD_DIMEN);
+    }
+
+    /* reset position to the correct HDU if necessary */
+    if ((*fptr)->HDUposition != ((*fptr)->Fptr)->curhdu)
+        ffmahd(*fptr, ((*fptr)->HDUposition) + 1, NULL, status);
+
+    if (imagetype == TBYTE)
+        bitpix = BYTE_IMG;
+    else if (imagetype == TSHORT)
+        bitpix = SHORT_IMG;
+    else if (imagetype == TINT)
+        bitpix = LONG_IMG;
+    else if (imagetype == TFLOAT)
+        bitpix = FLOAT_IMG;
+    else if (imagetype == TDOUBLE)
+        bitpix = DOUBLE_IMG;
+    else
+        return(*status = BAD_DATATYPE);
+
+    
+    /*    Calculate the binning parameters:    */
+    /*   columm numbers, axes length, min values,  max values, and binsizes.  */
+
+    if (fits_calc_binning(
+      *fptr, naxis, colname, minin, maxin, binsizein, minname, maxname, binname,
+      colnum,  haxes, amin, amax, binsize, status) > 0)
+    {
+        ffpmsg("failed to determine binning parameters");
+        return(*status);
+    }
+ 
+    /* get the histogramming weighting factor, if any */
+    if (*wtcol)
+    {
+        /* first, look for a keyword with the weight value */
+        if (ffgky(*fptr, TFLOAT, wtcol, &weight, NULL, status) )
+        {
+            /* not a keyword, so look for column with this name */
+            *status = 0;
+
+            /* get the column number in the table */
+            if (ffgcno(*fptr, CASEINSEN, wtcol, &wtcolnum, status) > 0)
+            {
+               ffpmsg(
+               "keyword or column for histogram weights doesn't exist: ");
+               ffpmsg(wtcol);
+               return(*status);
+            }
+
+            weight = FLOATNULLVALUE;
+        }
+    }
+    else
+        weight = (float) weightin;
+
+    if (weight <= 0. && weight != FLOATNULLVALUE)
+    {
+        ffpmsg("Illegal histogramming weighting factor <= 0.");
+        return(*status = URL_PARSE_ERROR);
+    }
+
+    if (recip && weight != FLOATNULLVALUE)
+       /* take reciprocal of weight */
+       weight = (float) (1.0 / weight);
+
+    /* size of histogram is now known, so create temp output file */
+    if (fits_create_file(&histptr, outfile, status) > 0)
+    {
+        ffpmsg("failed to create temp output file for histogram");
+        return(*status);
+    }
+
+    /* create output FITS image HDU */
+    if (ffcrim(histptr, bitpix, naxis, haxes, status) > 0)
+    {
+        ffpmsg("failed to create output histogram FITS image");
+        return(*status);
+    }
+
+    /* copy header keywords, converting pixel list WCS keywords to image WCS form */
+    if (fits_copy_pixlist2image(*fptr, histptr, 9, naxis, colnum, status) > 0)
+    {
+        ffpmsg("failed to copy pixel list keywords to new histogram header");
+        return(*status);
+    }
+
+    /* if the table columns have no WCS keywords, then write default keywords */
+    fits_write_keys_histo(*fptr, histptr, naxis, colnum, status);
+    
+    /* update the WCS keywords for the ref. pixel location, and pixel size */
+    fits_rebin_wcs(histptr, naxis, amin, binsize,  status);      
+    
+    /* now compute the output image by binning the column values */
+    if (fits_make_hist(*fptr, histptr, bitpix, naxis, haxes, colnum, amin, amax,
+        binsize, weight, wtcolnum, recip, selectrow, status) > 0)
+    {
+        ffpmsg("failed to calculate new histogram values");
+        return(*status);
+    }
+              
+    /* finally, close the original file and return ptr to the new image */
+    ffclos(*fptr, status);
+    *fptr = histptr;
+
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int ffhist(fitsfile **fptr,  /* IO - pointer to table with X and Y cols;    */
+                             /*     on output, points to histogram image    */
+           char *outfile,    /* I - name for the output histogram file      */
+           int imagetype,    /* I - datatype for image: TINT, TSHORT, etc   */
+           int naxis,        /* I - number of axes in the histogram image   */
+           char colname[4][FLEN_VALUE],   /* I - column names               */
+           double *minin,     /* I - minimum histogram value, for each axis */
+           double *maxin,     /* I - maximum histogram value, for each axis */
+           double *binsizein, /* I - bin size along each axis               */
+           char minname[4][FLEN_VALUE], /* I - optional keywords for min    */
+           char maxname[4][FLEN_VALUE], /* I - optional keywords for max    */
+           char binname[4][FLEN_VALUE], /* I - optional keywords for binsize */
+           double weightin,        /* I - binning weighting factor          */
+           char wtcol[FLEN_VALUE], /* I - optional keyword or col for weight*/
+           int recip,              /* I - use reciprocal of the weight?     */
+           char *selectrow,        /* I - optional array (length = no. of   */
+                             /* rows in the table).  If the element is true */
+                             /* then the corresponding row of the table will*/
+                             /* be included in the histogram, otherwise the */
+                             /* row will be skipped.  Ingnored if *selectrow*/
+                             /* is equal to NULL.                           */
+           int *status)
+{
+    int ii, datatype, repeat, imin, imax, ibin, bitpix, tstatus, use_datamax = 0;
+    long haxes[4];
+    fitsfile *histptr;
+    char errmsg[FLEN_ERRMSG], keyname[FLEN_KEYWORD], card[FLEN_CARD];
+    tcolumn *colptr;
+    iteratorCol imagepars[1];
+    int n_cols = 1, nkeys;
+    long  offset = 0;
+    long n_per_loop = -1;  /* force whole array to be passed at one time */
+    histType histData;    /* Structure holding histogram info for iterator */
+    
+    float amin[4], amax[4], binsize[4], maxbin[4];
+    float datamin = FLOATNULLVALUE, datamax = FLOATNULLVALUE;
+    char svalue[FLEN_VALUE];
+    double dvalue;
+    char cpref[4][FLEN_VALUE];
+    char *cptr;
+
+    if (*status > 0)
+        return(*status);
+
+    if (naxis > 4)
+    {
+        ffpmsg("histogram has more than 4 dimensions");
+        return(*status = BAD_DIMEN);
+    }
+
+    /* reset position to the correct HDU if necessary */
+    if ((*fptr)->HDUposition != ((*fptr)->Fptr)->curhdu)
+        ffmahd(*fptr, ((*fptr)->HDUposition) + 1, NULL, status);
+
+    histData.tblptr     = *fptr;
+    histData.himagetype = imagetype;
+    histData.haxis      = naxis;
+    histData.rowselector = selectrow;
+
+    if (imagetype == TBYTE)
+        bitpix = BYTE_IMG;
+    else if (imagetype == TSHORT)
+        bitpix = SHORT_IMG;
+    else if (imagetype == TINT)
+        bitpix = LONG_IMG;
+    else if (imagetype == TFLOAT)
+        bitpix = FLOAT_IMG;
+    else if (imagetype == TDOUBLE)
+        bitpix = DOUBLE_IMG;
+    else
+        return(*status = BAD_DATATYPE);
+
+    /* The CPREF keyword, if it exists, gives the preferred columns. */
+    /* Otherwise, assume "X", "Y", "Z", and "T"  */
+
+    tstatus = 0;
+    ffgky(*fptr, TSTRING, "CPREF", cpref[0], NULL, &tstatus);
+
+    if (!tstatus)
+    {
+        /* Preferred column names are given;  separate them */
+        cptr = cpref[0];
+
+        /* the first preferred axis... */
+        while (*cptr != ',' && *cptr != '\0')
+           cptr++;
+
+        if (*cptr != '\0')
+        {
+           *cptr = '\0';
+           cptr++;
+           while (*cptr == ' ')
+               cptr++;
+
+           strcpy(cpref[1], cptr);
+           cptr = cpref[1];
+
+          /* the second preferred axis... */
+          while (*cptr != ',' && *cptr != '\0')
+             cptr++;
+
+          if (*cptr != '\0')
+          {
+             *cptr = '\0';
+             cptr++;
+             while (*cptr == ' ')
+                 cptr++;
+
+             strcpy(cpref[2], cptr);
+             cptr = cpref[2];
+
+            /* the third preferred axis... */
+            while (*cptr != ',' && *cptr != '\0')
+               cptr++;
+
+            if (*cptr != '\0')
+            {
+               *cptr = '\0';
+               cptr++;
+               while (*cptr == ' ')
+                   cptr++;
+
+               strcpy(cpref[3], cptr);
+
+            }
+          }
+        }
+    }
+
+    for (ii = 0; ii < naxis; ii++)
+    {
+
+      /* get the min, max, and binsize values from keywords, if specified */
+
+      if (*minname[ii])
+      {
+         if (ffgky(*fptr, TDOUBLE, minname[ii], &minin[ii], NULL, status) )
+         {
+             ffpmsg("error reading histogramming minimum keyword");
+             ffpmsg(minname[ii]);
+             return(*status);
+         }
+      }
+
+      if (*maxname[ii])
+      {
+         if (ffgky(*fptr, TDOUBLE, maxname[ii], &maxin[ii], NULL, status) )
+         {
+             ffpmsg("error reading histogramming maximum keyword");
+             ffpmsg(maxname[ii]);
+             return(*status);
+         }
+      }
+
+      if (*binname[ii])
+      {
+         if (ffgky(*fptr, TDOUBLE, binname[ii], &binsizein[ii], NULL, status) )
+         {
+             ffpmsg("error reading histogramming binsize keyword");
+             ffpmsg(binname[ii]);
+             return(*status);
+         }
+      }
+
+      if (binsizein[ii] == 0.)
+      {
+        ffpmsg("error: histogram binsize = 0");
+        return(*status = ZERO_SCALE);
+      }
+
+      if (*colname[ii] == '\0')
+      {
+         strcpy(colname[ii], cpref[ii]); /* try using the preferred column */
+         if (*colname[ii] == '\0')
+         {
+           if (ii == 0)
+              strcpy(colname[ii], "X");
+           else if (ii == 1)
+              strcpy(colname[ii], "Y");
+           else if (ii == 2)
+              strcpy(colname[ii], "Z");
+           else if (ii == 3)
+              strcpy(colname[ii], "T");
+         }
+      }
+
+      /* get the column number in the table */
+      if (ffgcno(*fptr, CASEINSEN, colname[ii], histData.hcolnum+ii, status)
+              > 0)
+      {
+        strcpy(errmsg, "column for histogram axis doesn't exist: ");
+        strcat(errmsg, colname[ii]);
+        ffpmsg(errmsg);
+        return(*status);
+      }
+
+      colptr = ((*fptr)->Fptr)->tableptr;
+      colptr += (histData.hcolnum[ii] - 1);
+
+      repeat = (int) colptr->trepeat;  /* vector repeat factor of the column */
+      if (repeat > 1)
+      {
+        strcpy(errmsg, "Can't bin a vector column: ");
+        strcat(errmsg, colname[ii]);
+        ffpmsg(errmsg);
+        return(*status = BAD_DATATYPE);
+      }
+
+      /* get the datatype of the column */
+      fits_get_coltype(*fptr, histData.hcolnum[ii], &datatype,
+         NULL, NULL, status);
+
+      if (datatype < 0 || datatype == TSTRING)
+      {
+        strcpy(errmsg, "Inappropriate datatype; can't bin this column: ");
+        strcat(errmsg, colname[ii]);
+        ffpmsg(errmsg);
+        return(*status = BAD_DATATYPE);
+      }
+
+      /* use TLMINn and TLMAXn keyword values if min and max were not given */
+      /* else use actual data min and max if TLMINn and TLMAXn don't exist */
+ 
+      if (minin[ii] == DOUBLENULLVALUE)
+      {
+        ffkeyn("TLMIN", histData.hcolnum[ii], keyname, status);
+        if (ffgky(*fptr, TFLOAT, keyname, amin+ii, NULL, status) > 0)
+        {
+            /* use actual data minimum value for the histogram minimum */
+            *status = 0;
+            if (fits_get_col_minmax(*fptr, histData.hcolnum[ii], amin+ii, &datamax, status) > 0)
+            {
+                strcpy(errmsg, "Error calculating datamin and datamax for column: ");
+                strcat(errmsg, colname[ii]);
+                ffpmsg(errmsg);
+                return(*status);
+            }
+         }
+      }
+      else
+      {
+        amin[ii] = (float) minin[ii];
+      }
+
+      if (maxin[ii] == DOUBLENULLVALUE)
+      {
+        ffkeyn("TLMAX", histData.hcolnum[ii], keyname, status);
+        if (ffgky(*fptr, TFLOAT, keyname, &amax[ii], NULL, status) > 0)
+        {
+          *status = 0;
+          if(datamax != FLOATNULLVALUE)  /* already computed max value */
+          {
+             amax[ii] = datamax;
+          }
+          else
+          {
+             /* use actual data maximum value for the histogram maximum */
+             if (fits_get_col_minmax(*fptr, histData.hcolnum[ii], &datamin, &amax[ii], status) > 0)
+             {
+                 strcpy(errmsg, "Error calculating datamin and datamax for column: ");
+                 strcat(errmsg, colname[ii]);
+                 ffpmsg(errmsg);
+                 return(*status);
+             }
+          }
+        }
+        use_datamax = 1;  /* flag that the max was determined by the data values */
+                          /* and not specifically set by the calling program */
+      }
+      else
+      {
+        amax[ii] = (float) maxin[ii];
+      }
+
+      /* use TDBINn keyword or else 1 if bin size is not given */
+      if (binsizein[ii] == DOUBLENULLVALUE)
+      {
+         tstatus = 0;
+         ffkeyn("TDBIN", histData.hcolnum[ii], keyname, &tstatus);
+
+         if (ffgky(*fptr, TDOUBLE, keyname, binsizein + ii, NULL, &tstatus) > 0)
+         {
+	    /* make at least 10 bins */
+            binsizein[ii] = (amax[ii] - amin[ii]) / 10. ;
+            if (binsizein[ii] > 1.)
+                binsizein[ii] = 1.;  /* use default bin size */
+         }
+      }
+
+      if ( (amin[ii] > amax[ii] && binsizein[ii] > 0. ) ||
+           (amin[ii] < amax[ii] && binsizein[ii] < 0. ) )
+          binsize[ii] = (float) -binsizein[ii];  /* reverse the sign of binsize */
+      else
+          binsize[ii] =  (float) binsizein[ii];  /* binsize has the correct sign */
+
+      ibin = (int) binsize[ii];
+      imin = (int) amin[ii];
+      imax = (int) amax[ii];
+
+      /* Determine the range and number of bins in the histogram. This  */
+      /* depends on whether the input columns are integer or floats, so */
+      /* treat each case separately.                                    */
+
+      if (datatype <= TLONG && (float) imin == amin[ii] &&
+                               (float) imax == amax[ii] &&
+                               (float) ibin == binsize[ii] )
+      {
+        /* This is an integer column and integer limits were entered. */
+        /* Shift the lower and upper histogramming limits by 0.5, so that */
+        /* the values fall in the center of the bin, not on the edge. */
+
+        haxes[ii] = (imax - imin) / ibin + 1;  /* last bin may only */
+                                               /* be partially full */
+        maxbin[ii] = (float) (haxes[ii] + 1.);  /* add 1. instead of .5 to avoid roundoff */
+
+        if (amin[ii] < amax[ii])
+        {
+          amin[ii] = (float) (amin[ii] - 0.5);
+          amax[ii] = (float) (amax[ii] + 0.5);
+        }
+        else
+        {
+          amin[ii] = (float) (amin[ii] + 0.5);
+          amax[ii] = (float) (amax[ii] - 0.5);
+        }
+      }
+      else if (use_datamax)  
+      {
+        /* Either the column datatype and/or the limits are floating point, */
+        /* and the histogram limits are being defined by the min and max */
+        /* values of the array.  Add 1 to the number of histogram bins to */
+        /* make sure that pixels that are equal to the maximum or are */
+        /* in the last partial bin are included.  */
+
+        maxbin[ii] = (amax[ii] - amin[ii]) / binsize[ii]; 
+        haxes[ii] = (long) (maxbin[ii] + 1);
+      }
+      else  
+      {
+        /*  float datatype column and/or limits, and the maximum value to */
+        /*  include in the histogram is specified by the calling program. */
+        /*  The lower limit is inclusive, but upper limit is exclusive    */
+        maxbin[ii] = (amax[ii] - amin[ii]) / binsize[ii];
+        haxes[ii] = (long) maxbin[ii];
+
+        if (amin[ii] < amax[ii])
+        {
+          if (amin[ii] + (haxes[ii] * binsize[ii]) < amax[ii])
+            haxes[ii]++;   /* need to include another partial bin */
+        }
+        else
+        {
+          if (amin[ii] + (haxes[ii] * binsize[ii]) > amax[ii])
+            haxes[ii]++;   /* need to include another partial bin */
+        }
+      }
+    }
+
+       /* get the histogramming weighting factor */
+    if (*wtcol)
+    {
+        /* first, look for a keyword with the weight value */
+        if (ffgky(*fptr, TFLOAT, wtcol, &histData.weight, NULL, status) )
+        {
+            /* not a keyword, so look for column with this name */
+            *status = 0;
+
+            /* get the column number in the table */
+            if (ffgcno(*fptr, CASEINSEN, wtcol, &histData.wtcolnum, status) > 0)
+            {
+               ffpmsg(
+               "keyword or column for histogram weights doesn't exist: ");
+               ffpmsg(wtcol);
+               return(*status);
+            }
+
+            histData.weight = FLOATNULLVALUE;
+        }
+    }
+    else
+        histData.weight = (float) weightin;
+
+    if (histData.weight <= 0. && histData.weight != FLOATNULLVALUE)
+    {
+        ffpmsg("Illegal histogramming weighting factor <= 0.");
+        return(*status = URL_PARSE_ERROR);
+    }
+
+    if (recip && histData.weight != FLOATNULLVALUE)
+       /* take reciprocal of weight */
+       histData.weight = (float) (1.0 / histData.weight);
+
+    histData.wtrecip = recip;
+        
+    /* size of histogram is now known, so create temp output file */
+    if (ffinit(&histptr, outfile, status) > 0)
+    {
+        ffpmsg("failed to create temp output file for histogram");
+        return(*status);
+    }
+
+    if (ffcrim(histptr, bitpix, histData.haxis, haxes, status) > 0)
+    {
+        ffpmsg("failed to create primary array histogram in temp file");
+        ffclos(histptr, status);
+        return(*status);
+    }
+
+    /* copy all non-structural keywords from the table to the image */
+    fits_get_hdrspace(*fptr, &nkeys, NULL, status);
+    for (ii = 1; ii <= nkeys; ii++)
+    {
+       fits_read_record(*fptr, ii, card, status);
+       if (fits_get_keyclass(card) >= 120)
+           fits_write_record(histptr, card, status);
+    }           
+
+    /* Set global variables with histogram parameter values.    */
+    /* Use separate scalar variables rather than arrays because */
+    /* it is more efficient when computing the histogram.       */
+
+    histData.amin1 = amin[0];
+    histData.maxbin1 = maxbin[0];
+    histData.binsize1 = binsize[0];
+    histData.haxis1 = haxes[0];
+
+    if (histData.haxis > 1)
+    {
+      histData.amin2 = amin[1];
+      histData.maxbin2 = maxbin[1];
+      histData.binsize2 = binsize[1];
+      histData.haxis2 = haxes[1];
+
+      if (histData.haxis > 2)
+      {
+        histData.amin3 = amin[2];
+        histData.maxbin3 = maxbin[2];
+        histData.binsize3 = binsize[2];
+        histData.haxis3 = haxes[2];
+
+        if (histData.haxis > 3)
+        {
+          histData.amin4 = amin[3];
+          histData.maxbin4 = maxbin[3];
+          histData.binsize4 = binsize[3];
+          histData.haxis4 = haxes[3];
+        }
+      }
+    }
+
+    /* define parameters of image for the iterator function */
+    fits_iter_set_file(imagepars, histptr);        /* pointer to image */
+    fits_iter_set_datatype(imagepars, imagetype);  /* image datatype   */
+    fits_iter_set_iotype(imagepars, OutputCol);    /* image is output  */
+
+    /* call the iterator function to write out the histogram image */
+    if (fits_iterate_data(n_cols, imagepars, offset, n_per_loop,
+                          ffwritehisto, (void*)&histData, status) )
+         return(*status);
+
+    /* write the World Coordinate System (WCS) keywords */
+    /* create default values if WCS keywords are not present in the table */
+    for (ii = 0; ii < histData.haxis; ii++)
+    {
+     /*  CTYPEn  */
+       tstatus = 0;
+       ffkeyn("TCTYP", histData.hcolnum[ii], keyname, &tstatus);
+       ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus);
+       if (tstatus)
+       {               /* just use column name as the type */
+          tstatus = 0;
+          ffkeyn("TTYPE", histData.hcolnum[ii], keyname, &tstatus);
+          ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus);
+       }
+
+       if (!tstatus)
+       {
+        ffkeyn("CTYPE", ii + 1, keyname, &tstatus);
+        ffpky(histptr, TSTRING, keyname, svalue, "Coordinate Type", &tstatus);
+       }
+       else
+          tstatus = 0;
+
+     /*  CUNITn  */
+       ffkeyn("TCUNI", histData.hcolnum[ii], keyname, &tstatus);
+       ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus);
+       if (tstatus)
+       {         /* use the column units */
+          tstatus = 0;
+          ffkeyn("TUNIT", histData.hcolnum[ii], keyname, &tstatus);
+          ffgky(*fptr, TSTRING, keyname, svalue, NULL, &tstatus);
+       }
+
+       if (!tstatus)
+       {
+        ffkeyn("CUNIT", ii + 1, keyname, &tstatus);
+        ffpky(histptr, TSTRING, keyname, svalue, "Coordinate Units", &tstatus);
+       }
+       else
+         tstatus = 0;
+
+     /*  CRPIXn  - Reference Pixel  */
+       ffkeyn("TCRPX", histData.hcolnum[ii], keyname, &tstatus);
+       ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus);
+       if (tstatus)
+       {
+         dvalue = 1.0; /* choose first pixel in new image as ref. pix. */
+         tstatus = 0;
+       }
+       else
+       {
+           /* calculate locate of the ref. pix. in the new image */
+           dvalue = (dvalue - amin[ii]) / binsize[ii] + .5;
+       }
+
+       ffkeyn("CRPIX", ii + 1, keyname, &tstatus);
+       ffpky(histptr, TDOUBLE, keyname, &dvalue, "Reference Pixel", &tstatus);
+
+     /*  CRVALn - Value at the location of the reference pixel */
+       ffkeyn("TCRVL", histData.hcolnum[ii], keyname, &tstatus);
+       ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus);
+       if (tstatus)
+       {
+         /* calculate value at ref. pix. location (at center of 1st pixel) */
+         dvalue = amin[ii] + binsize[ii]/2.;
+         tstatus = 0;
+       }
+
+       ffkeyn("CRVAL", ii + 1, keyname, &tstatus);
+       ffpky(histptr, TDOUBLE, keyname, &dvalue, "Reference Value", &tstatus);
+
+     /*  CDELTn - unit size of pixels  */
+       ffkeyn("TCDLT", histData.hcolnum[ii], keyname, &tstatus);
+       ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus);
+       if (tstatus)
+       {
+         dvalue = 1.0;  /* use default pixel size */
+         tstatus = 0;
+       }
+
+       dvalue = dvalue * binsize[ii];
+       ffkeyn("CDELT", ii + 1, keyname, &tstatus);
+       ffpky(histptr, TDOUBLE, keyname, &dvalue, "Pixel size", &tstatus);
+
+     /*  CROTAn - Rotation angle (degrees CCW)  */
+     /*  There should only be a CROTA2 keyword, and only for 2+ D images */
+       if (ii == 1)
+       {
+         ffkeyn("TCROT", histData.hcolnum[ii], keyname, &tstatus);
+         ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus);
+         if (!tstatus && dvalue != 0.)  /* only write keyword if angle != 0 */
+         {
+           ffkeyn("CROTA", ii + 1, keyname, &tstatus);
+           ffpky(histptr, TDOUBLE, keyname, &dvalue,
+                 "Rotation angle", &tstatus);
+         }
+         else
+         {
+            /* didn't find CROTA for the 2nd axis, so look for one */
+            /* on the first axis */
+           tstatus = 0;
+           ffkeyn("TCROT", histData.hcolnum[0], keyname, &tstatus);
+           ffgky(*fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus);
+           if (!tstatus && dvalue != 0.)  /* only write keyword if angle != 0 */
+           {
+             dvalue *= -1.;   /* negate the value, because mirror image */
+             ffkeyn("CROTA", ii + 1, keyname, &tstatus);
+             ffpky(histptr, TDOUBLE, keyname, &dvalue,
+                   "Rotation angle", &tstatus);
+           }
+         }
+       }
+    }
+
+    /* convert any TPn_k keywords to PCi_j; the value remains unchanged */
+    /* also convert any TCn_k to CDi_j; the value is modified by n binning size */
+    /* This is a bit of a kludge, and only works for 2D WCS */
+
+    if (histData.haxis == 2) {
+
+      /* PC1_1 */
+      tstatus = 0;
+      ffkeyn("TP", histData.hcolnum[0], card, &tstatus);
+      strcat(card,"_");
+      ffkeyn(card, histData.hcolnum[0], keyname, &tstatus);
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) 
+         ffpky(histptr, TDOUBLE, "PC1_1", &dvalue, card, &tstatus);
+
+      tstatus = 0;
+      keyname[1] = 'C';
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) {
+         dvalue *=  binsize[0];
+         ffpky(histptr, TDOUBLE, "CD1_1", &dvalue, card, &tstatus);
+      }
+
+      /* PC1_2 */
+      tstatus = 0;
+      ffkeyn("TP", histData.hcolnum[0], card, &tstatus);
+      strcat(card,"_");
+      ffkeyn(card, histData.hcolnum[1], keyname, &tstatus);
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) 
+         ffpky(histptr, TDOUBLE, "PC1_2", &dvalue, card, &tstatus);
+ 
+      tstatus = 0;
+      keyname[1] = 'C';
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) {
+        dvalue *=  binsize[0];
+        ffpky(histptr, TDOUBLE, "CD1_2", &dvalue, card, &tstatus);
+      }
+       
+      /* PC2_1 */
+      tstatus = 0;
+      ffkeyn("TP", histData.hcolnum[1], card, &tstatus);
+      strcat(card,"_");
+      ffkeyn(card, histData.hcolnum[0], keyname, &tstatus);
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) 
+         ffpky(histptr, TDOUBLE, "PC2_1", &dvalue, card, &tstatus);
+ 
+      tstatus = 0;
+      keyname[1] = 'C';
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) {
+         dvalue *=  binsize[1];
+         ffpky(histptr, TDOUBLE, "CD2_1", &dvalue, card, &tstatus);
+      }
+       
+       /* PC2_2 */
+      tstatus = 0;
+      ffkeyn("TP", histData.hcolnum[1], card, &tstatus);
+      strcat(card,"_");
+      ffkeyn(card, histData.hcolnum[1], keyname, &tstatus);
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) 
+         ffpky(histptr, TDOUBLE, "PC2_2", &dvalue, card, &tstatus);
+        
+      tstatus = 0;
+      keyname[1] = 'C';
+      ffgky(*fptr, TDOUBLE, keyname, &dvalue, card, &tstatus);
+      if (!tstatus) {
+         dvalue *=  binsize[1];
+         ffpky(histptr, TDOUBLE, "CD2_2", &dvalue, card, &tstatus);
+      }
+    }   
+       
+    /* finally, close the original file and return ptr to the new image */
+    ffclos(*fptr, status);
+    *fptr = histptr;
+
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int fits_calc_binning(
+      fitsfile *fptr,  /* IO - pointer to table to be binned      ;       */
+      int naxis,       /* I - number of axes/columns in the binned image  */
+      char colname[4][FLEN_VALUE],   /* I - optional column names         */
+      double *minin,     /* I - optional lower bound value for each axis  */
+      double *maxin,     /* I - optional upper bound value, for each axis */
+      double *binsizein, /* I - optional bin size along each axis         */
+      char minname[4][FLEN_VALUE], /* I - optional keywords for min       */
+      char maxname[4][FLEN_VALUE], /* I - optional keywords for max       */
+      char binname[4][FLEN_VALUE], /* I - optional keywords for binsize   */
+
+    /* The returned parameters for each axis of the n-dimensional histogram are */
+
+      int *colnum,     /* O - column numbers, to be binned */
+      long *haxes,     /* O - number of bins in each histogram axis */
+      float *amin,     /* O - lower bound of the histogram axes */
+      float *amax,     /* O - upper bound of the histogram axes */
+      float *binsize,  /* O - width of histogram bins/pixels on each axis */
+      int *status)
+/*_
+    Calculate the actual binning parameters, based on various user input
+    options.
+*/
+{
+    tcolumn *colptr;
+    char *cptr, cpref[4][FLEN_VALUE];
+    char errmsg[FLEN_ERRMSG], keyname[FLEN_KEYWORD];
+    int tstatus, ii;
+    int datatype, repeat, imin, imax, ibin,  use_datamax = 0;
+    float datamin, datamax;
+
+    /* check inputs */
+    
+    if (*status > 0)
+        return(*status);
+
+    if (naxis > 4)
+    {
+        ffpmsg("histograms with more than 4 dimensions are not supported");
+        return(*status = BAD_DIMEN);
+    }
+
+    /* reset position to the correct HDU if necessary */
+    if ((fptr)->HDUposition != ((fptr)->Fptr)->curhdu)
+        ffmahd(fptr, ((fptr)->HDUposition) + 1, NULL, status);
+    
+    /* ============================================================= */
+    /* The CPREF keyword, if it exists, gives the preferred columns. */
+    /* Otherwise, assume "X", "Y", "Z", and "T"  */
+
+    *cpref[0] = '\0';
+    *cpref[1] = '\0';
+    *cpref[2] = '\0';
+    *cpref[3] = '\0';
+
+    tstatus = 0;
+    ffgky(fptr, TSTRING, "CPREF", cpref[0], NULL, &tstatus);
+
+    if (!tstatus)
+    {
+        /* Preferred column names are given;  separate them */
+        cptr = cpref[0];
+
+        /* the first preferred axis... */
+        while (*cptr != ',' && *cptr != '\0')
+           cptr++;
+
+        if (*cptr != '\0')
+        {
+           *cptr = '\0';
+           cptr++;
+           while (*cptr == ' ')
+               cptr++;
+
+           strcpy(cpref[1], cptr);
+           cptr = cpref[1];
+
+          /* the second preferred axis... */
+          while (*cptr != ',' && *cptr != '\0')
+             cptr++;
+
+          if (*cptr != '\0')
+          {
+             *cptr = '\0';
+             cptr++;
+             while (*cptr == ' ')
+                 cptr++;
+
+             strcpy(cpref[2], cptr);
+             cptr = cpref[2];
+
+            /* the third preferred axis... */
+            while (*cptr != ',' && *cptr != '\0')
+               cptr++;
+
+            if (*cptr != '\0')
+            {
+               *cptr = '\0';
+               cptr++;
+               while (*cptr == ' ')
+                   cptr++;
+
+               strcpy(cpref[3], cptr);
+
+            }
+          }
+        }
+    }
+
+    /* ============================================================= */
+    /* Main Loop for calculating parameters for each column          */
+
+    for (ii = 0; ii < naxis; ii++)
+    {
+
+      /* =========================================================== */
+      /* Determine column Number, based on, in order of priority,
+         1  input column name, or
+	 2  name given by CPREF keyword, or
+	 3  assume X, Y, Z and T for the name
+      */
+	  
+      if (*colname[ii] == '\0')
+      {
+         strcpy(colname[ii], cpref[ii]); /* try using the preferred column */
+         if (*colname[ii] == '\0')
+         {
+           if (ii == 0)
+              strcpy(colname[ii], "X");
+           else if (ii == 1)
+              strcpy(colname[ii], "Y");
+           else if (ii == 2)
+              strcpy(colname[ii], "Z");
+           else if (ii == 3)
+              strcpy(colname[ii], "T");
+         }
+      }
+
+      /* get the column number in the table */
+      if (ffgcno(fptr, CASEINSEN, colname[ii], colnum+ii, status)
+              > 0)
+      {
+          strcpy(errmsg, "column for histogram axis doesn't exist: ");
+          strcat(errmsg, colname[ii]);
+          ffpmsg(errmsg);
+          return(*status);
+      }
+
+      /* ================================================================ */
+      /* check tha column is not a vector or a string                     */
+
+      colptr = ((fptr)->Fptr)->tableptr;
+      colptr += (colnum[ii] - 1);
+
+      repeat = (int) colptr->trepeat;  /* vector repeat factor of the column */
+      if (repeat > 1)
+      {
+        strcpy(errmsg, "Can't bin a vector column: ");
+        strcat(errmsg, colname[ii]);
+        ffpmsg(errmsg);
+        return(*status = BAD_DATATYPE);
+      }
+
+      /* get the datatype of the column */
+      fits_get_coltype(fptr, colnum[ii], &datatype,
+         NULL, NULL, status);
+
+      if (datatype < 0 || datatype == TSTRING)
+      {
+        strcpy(errmsg, "Inappropriate datatype; can't bin this column: ");
+        strcat(errmsg, colname[ii]);
+        ffpmsg(errmsg);
+        return(*status = BAD_DATATYPE);
+      }
+
+      /* ================================================================ */
+      /* get the minimum value */
+
+      datamin = FLOATNULLVALUE;
+      datamax = FLOATNULLVALUE;
+      
+      if (*minname[ii])
+      {
+         if (ffgky(fptr, TDOUBLE, minname[ii], &minin[ii], NULL, status) )
+         {
+             ffpmsg("error reading histogramming minimum keyword");
+             ffpmsg(minname[ii]);
+             return(*status);
+         }
+      }
+
+      if (minin[ii] != DOUBLENULLVALUE)
+      {
+        amin[ii] = (float) minin[ii];
+      }
+      else
+      {
+        ffkeyn("TLMIN", colnum[ii], keyname, status);
+        if (ffgky(fptr, TFLOAT, keyname, amin+ii, NULL, status) > 0)
+        {
+            /* use actual data minimum value for the histogram minimum */
+            *status = 0;
+            if (fits_get_col_minmax(fptr, colnum[ii], amin+ii, &datamax, status) > 0)
+            {
+                strcpy(errmsg, "Error calculating datamin and datamax for column: ");
+                strcat(errmsg, colname[ii]);
+                ffpmsg(errmsg);
+                return(*status);
+            }
+         }
+      }
+
+      /* ================================================================ */
+      /* get the maximum value */
+
+      if (*maxname[ii])
+      {
+         if (ffgky(fptr, TDOUBLE, maxname[ii], &maxin[ii], NULL, status) )
+         {
+             ffpmsg("error reading histogramming maximum keyword");
+             ffpmsg(maxname[ii]);
+             return(*status);
+         }
+      }
+
+      if (maxin[ii] != DOUBLENULLVALUE)
+      {
+        amax[ii] = (float) maxin[ii];
+      }
+      else
+      {
+        ffkeyn("TLMAX", colnum[ii], keyname, status);
+        if (ffgky(fptr, TFLOAT, keyname, &amax[ii], NULL, status) > 0)
+        {
+          *status = 0;
+          if(datamax != FLOATNULLVALUE)  /* already computed max value */
+          {
+             amax[ii] = datamax;
+          }
+          else
+          {
+             /* use actual data maximum value for the histogram maximum */
+             if (fits_get_col_minmax(fptr, colnum[ii], &datamin, &amax[ii], status) > 0)
+             {
+                 strcpy(errmsg, "Error calculating datamin and datamax for column: ");
+                 strcat(errmsg, colname[ii]);
+                 ffpmsg(errmsg);
+                 return(*status);
+             }
+          }
+        }
+        use_datamax = 1;  /* flag that the max was determined by the data values */
+                          /* and not specifically set by the calling program */
+      }
+
+
+      /* ================================================================ */
+      /* determine binning size and range                                 */
+
+      if (*binname[ii])
+      {
+         if (ffgky(fptr, TDOUBLE, binname[ii], &binsizein[ii], NULL, status) )
+         {
+             ffpmsg("error reading histogramming binsize keyword");
+             ffpmsg(binname[ii]);
+             return(*status);
+         }
+      }
+
+      if (binsizein[ii] == 0.)
+      {
+        ffpmsg("error: histogram binsize = 0");
+        return(*status = ZERO_SCALE);
+      }
+
+      /* use TDBINn keyword or else 1 if bin size is not given */
+      if (binsizein[ii] != DOUBLENULLVALUE)
+      { 
+         binsize[ii] = (float) binsizein[ii];
+      }
+      else
+      {
+         tstatus = 0;
+         ffkeyn("TDBIN", colnum[ii], keyname, &tstatus);
+
+         if (ffgky(fptr, TDOUBLE, keyname, binsizein + ii, NULL, &tstatus) > 0)
+         {
+	    /* make at least 10 bins */
+            binsize[ii] = (amax[ii] - amin[ii]) / 10.F ;
+            if (binsize[ii] > 1.)
+                binsize[ii] = 1.;  /* use default bin size */
+         }
+      }
+
+      /* ================================================================ */
+      /* if the min is greater than the max, make the binsize negative */
+      if ( (amin[ii] > amax[ii] && binsize[ii] > 0. ) ||
+           (amin[ii] < amax[ii] && binsize[ii] < 0. ) )
+          binsize[ii] =  -binsize[ii];  /* reverse the sign of binsize */
+
+
+      ibin = (int) binsize[ii];
+      imin = (int) amin[ii];
+      imax = (int) amax[ii];
+
+      /* Determine the range and number of bins in the histogram. This  */
+      /* depends on whether the input columns are integer or floats, so */
+      /* treat each case separately.                                    */
+
+      if (datatype <= TLONG && (float) imin == amin[ii] &&
+                               (float) imax == amax[ii] &&
+                               (float) ibin == binsize[ii] )
+      {
+        /* This is an integer column and integer limits were entered. */
+        /* Shift the lower and upper histogramming limits by 0.5, so that */
+        /* the values fall in the center of the bin, not on the edge. */
+
+        haxes[ii] = (imax - imin) / ibin + 1;  /* last bin may only */
+                                               /* be partially full */
+        if (amin[ii] < amax[ii])
+        {
+          amin[ii] = (float) (amin[ii] - 0.5);
+          amax[ii] = (float) (amax[ii] + 0.5);
+        }
+        else
+        {
+          amin[ii] = (float) (amin[ii] + 0.5);
+          amax[ii] = (float) (amax[ii] - 0.5);
+        }
+      }
+      else if (use_datamax)  
+      {
+        /* Either the column datatype and/or the limits are floating point, */
+        /* and the histogram limits are being defined by the min and max */
+        /* values of the array.  Add 1 to the number of histogram bins to */
+        /* make sure that pixels that are equal to the maximum or are */
+        /* in the last partial bin are included.  */
+
+        haxes[ii] = (long) (((amax[ii] - amin[ii]) / binsize[ii]) + 1.); 
+      }
+      else  
+      {
+        /*  float datatype column and/or limits, and the maximum value to */
+        /*  include in the histogram is specified by the calling program. */
+        /*  The lower limit is inclusive, but upper limit is exclusive    */
+        haxes[ii] = (long) ((amax[ii] - amin[ii]) / binsize[ii]);
+
+        if (amin[ii] < amax[ii])
+        {
+          if (amin[ii] + (haxes[ii] * binsize[ii]) < amax[ii])
+            haxes[ii]++;   /* need to include another partial bin */
+        }
+        else
+        {
+          if (amin[ii] + (haxes[ii] * binsize[ii]) > amax[ii])
+            haxes[ii]++;   /* need to include another partial bin */
+        }
+      }
+    }
+
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int fits_write_keys_histo(
+      fitsfile *fptr,   /* I - pointer to table to be binned              */
+      fitsfile *histptr,  /* I - pointer to output histogram image HDU      */
+      int naxis,        /* I - number of axes in the histogram image      */
+      int *colnum,      /* I - column numbers (array length = naxis)      */
+      int *status)     
+{      
+   /*  Write default WCS keywords in the output histogram image header */
+   /*  if the keywords do not already exist.   */
+
+    int ii, tstatus;
+    char keyname[FLEN_KEYWORD], svalue[FLEN_VALUE];
+    double dvalue;
+    
+    if (*status > 0)
+        return(*status);
+
+    for (ii = 0; ii < naxis; ii++)
+    {
+     /*  CTYPEn  */
+       tstatus = 0;
+       ffkeyn("CTYPE", ii+1, keyname, &tstatus);
+       ffgky(histptr, TSTRING, keyname, svalue, NULL, &tstatus);
+       
+       if (!tstatus) continue;  /* keyword already exists, so skip to next axis */
+       
+       /* use column name as the axis name */
+       tstatus = 0;
+       ffkeyn("TTYPE", colnum[ii], keyname, &tstatus);
+       ffgky(fptr, TSTRING, keyname, svalue, NULL, &tstatus);
+
+       if (!tstatus)
+       {
+         ffkeyn("CTYPE", ii + 1, keyname, &tstatus);
+         ffpky(histptr, TSTRING, keyname, svalue, "Coordinate Type", &tstatus);
+       }
+
+       /*  CUNITn,  use the column units */
+       tstatus = 0;
+       ffkeyn("TUNIT", colnum[ii], keyname, &tstatus);
+       ffgky(fptr, TSTRING, keyname, svalue, NULL, &tstatus);
+
+       if (!tstatus)
+       {
+         ffkeyn("CUNIT", ii + 1, keyname, &tstatus);
+         ffpky(histptr, TSTRING, keyname, svalue, "Coordinate Units", &tstatus);
+       }
+
+       /*  CRPIXn  - Reference Pixel choose first pixel in new image as ref. pix. */
+       dvalue = 1.0;
+       tstatus = 0;
+       ffkeyn("CRPIX", ii + 1, keyname, &tstatus);
+       ffpky(histptr, TDOUBLE, keyname, &dvalue, "Reference Pixel", &tstatus);
+
+       /*  CRVALn - Value at the location of the reference pixel */
+       dvalue = 1.0;
+       tstatus = 0;
+       ffkeyn("CRVAL", ii + 1, keyname, &tstatus);
+       ffpky(histptr, TDOUBLE, keyname, &dvalue, "Reference Value", &tstatus);
+
+       /*  CDELTn - unit size of pixels  */
+       dvalue = 1.0;  
+       tstatus = 0;
+       dvalue = 1.;
+       ffkeyn("CDELT", ii + 1, keyname, &tstatus);
+       ffpky(histptr, TDOUBLE, keyname, &dvalue, "Pixel size", &tstatus);
+
+    }
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int fits_rebin_wcs(
+      fitsfile *fptr,   /* I - pointer to table to be binned           */
+      int naxis,        /* I - number of axes in the histogram image   */
+      float *amin,        /* I - first pixel include in each axis        */
+      float *binsize,     /* I - binning factor for each axis            */
+      int *status)      
+{      
+   /*  Update the  WCS keywords that define the location of the reference */
+   /*  pixel, and the pixel size, along each axis.   */
+
+    int ii, jj, tstatus, reset ;
+    char keyname[FLEN_KEYWORD], svalue[FLEN_VALUE];
+    double dvalue;
+    
+    if (*status > 0)
+        return(*status);
+  
+    for (ii = 0; ii < naxis; ii++)
+    {
+       reset = 0;  /* flag to reset the reference pixel */
+       tstatus = 0;
+       ffkeyn("CRVAL", ii + 1, keyname, &tstatus);
+       /* get previous (pre-binning) value */
+       ffgky(fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); 
+       if (!tstatus && dvalue == 1.0)
+           reset = 1;
+
+       tstatus = 0;
+       /*  CRPIXn - update location of the ref. pix. in the binned image */
+       ffkeyn("CRPIX", ii + 1, keyname, &tstatus);
+
+       /* get previous (pre-binning) value */
+       ffgky(fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); 
+
+       if (!tstatus)
+       {
+           if (dvalue != 1.0)
+	      reset = 0;
+
+           /* updated value to give pixel location after binning */
+           dvalue = (dvalue - amin[ii]) / ((double) binsize[ii]) + .5;  
+
+           fits_modify_key_dbl(fptr, keyname, dvalue, -14, NULL, &tstatus);
+       } else {
+          reset = 0;
+       }
+
+       /*  CDELTn - update unit size of pixels  */
+       tstatus = 0;
+       ffkeyn("CDELT", ii + 1, keyname, &tstatus);
+
+       /* get previous (pre-binning) value */
+       ffgky(fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); 
+
+       if (!tstatus)
+       {
+           if (dvalue != 1.0)
+	      reset = 0;
+
+           /* updated to give post-binning value */
+           dvalue = dvalue * binsize[ii];  
+
+           fits_modify_key_dbl(fptr, keyname, dvalue, -14, NULL, &tstatus);
+       }
+       else
+       {   /* no CDELTn keyword, so look for a CDij keywords */
+          reset = 0;
+
+          for (jj = 0; jj < naxis; jj++)
+	  {
+             tstatus = 0;
+             ffkeyn("CD", jj + 1, svalue, &tstatus);
+	     strcat(svalue,"_");
+	     ffkeyn(svalue, ii + 1, keyname, &tstatus);
+
+             /* get previous (pre-binning) value */
+             ffgky(fptr, TDOUBLE, keyname, &dvalue, NULL, &tstatus); 
+
+             if (!tstatus)
+             {
+                /* updated to give post-binning value */
+               dvalue = dvalue * binsize[ii];  
+
+               fits_modify_key_dbl(fptr, keyname, dvalue, -14, NULL, &tstatus);
+             }
+	  }
+       }
+
+       if (reset) {
+          /* the original CRPIX, CRVAL, and CDELT keywords were all = 1.0 */
+	  /* In this special case, reset the reference pixel to be the */
+	  /* first pixel in the array (instead of possibly far off the array) */
+ 
+           dvalue = 1.0;
+           ffkeyn("CRPIX", ii + 1, keyname, &tstatus);
+           fits_modify_key_dbl(fptr, keyname, dvalue, -14, NULL, &tstatus);
+
+           ffkeyn("CRVAL", ii + 1, keyname, &tstatus);
+	   dvalue = amin[ii] + (binsize[ii] / 2.0);	  
+           fits_modify_key_dbl(fptr, keyname, dvalue, -14, NULL, &tstatus);
+	}
+
+    }
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+
+int fits_make_hist(fitsfile *fptr, /* IO - pointer to table with X and Y cols; */
+    fitsfile *histptr, /* I - pointer to output FITS image      */
+    int bitpix,       /* I - datatype for image: 16, 32, -32, etc    */
+    int naxis,        /* I - number of axes in the histogram image   */
+    long *naxes,      /* I - size of axes in the histogram image   */
+    int *colnum,    /* I - column numbers (array length = naxis)   */
+    float *amin,     /* I - minimum histogram value, for each axis */
+    float *amax,     /* I - maximum histogram value, for each axis */
+    float *binsize, /* I - bin size along each axis               */
+    float weight,        /* I - binning weighting factor          */
+    int wtcolnum, /* I - optional keyword or col for weight*/
+    int recip,              /* I - use reciprocal of the weight?     */
+    char *selectrow,        /* I - optional array (length = no. of   */
+                             /* rows in the table).  If the element is true */
+                             /* then the corresponding row of the table will*/
+                             /* be included in the histogram, otherwise the */
+                             /* row will be skipped.  Ingnored if *selectrow*/
+                             /* is equal to NULL.                           */
+    int *status)
+{		  
+    int ii, imagetype, datatype;
+    int n_cols = 1;
+    long imin, imax, ibin;
+    long  offset = 0;
+    long n_per_loop = -1;  /* force whole array to be passed at one time */
+    float taxes[4], tmin[4], tmax[4], tbin[4], maxbin[4];
+    histType histData;    /* Structure holding histogram info for iterator */
+    iteratorCol imagepars[1];
+
+    /* check inputs */
+    
+    if (*status > 0)
+        return(*status);
+
+    if (naxis > 4)
+    {
+        ffpmsg("histogram has more than 4 dimensions");
+        return(*status = BAD_DIMEN);
+    }
+
+    if   (bitpix == BYTE_IMG)
+         imagetype = TBYTE;
+    else if (bitpix == SHORT_IMG)
+         imagetype = TSHORT;
+    else if (bitpix == LONG_IMG)
+         imagetype = TINT;    
+    else if (bitpix == FLOAT_IMG)
+         imagetype = TFLOAT;    
+    else if (bitpix == DOUBLE_IMG)
+         imagetype = TDOUBLE;    
+    else
+        return(*status = BAD_DATATYPE);
+
+    /* reset position to the correct HDU if necessary */
+    if ((fptr)->HDUposition != ((fptr)->Fptr)->curhdu)
+        ffmahd(fptr, ((fptr)->HDUposition) + 1, NULL, status);
+
+    histData.weight     = weight;
+    histData.wtcolnum   = wtcolnum;
+    histData.wtrecip    = recip;
+    histData.tblptr     = fptr;
+    histData.himagetype = imagetype;
+    histData.haxis      = naxis;
+    histData.rowselector = selectrow;
+
+    for (ii = 0; ii < naxis; ii++)
+    {
+      taxes[ii] = (float) naxes[ii];
+      tmin[ii] = amin[ii];
+      tmax[ii] = amax[ii];
+      if ( (amin[ii] > amax[ii] && binsize[ii] > 0. ) ||
+           (amin[ii] < amax[ii] && binsize[ii] < 0. ) )
+          tbin[ii] =  -binsize[ii];  /* reverse the sign of binsize */
+      else
+          tbin[ii] =   binsize[ii];  /* binsize has the correct sign */
+          
+      imin = (long) tmin[ii];
+      imax = (long) tmax[ii];
+      ibin = (long) tbin[ii];
+    
+      /* get the datatype of the column */
+      fits_get_coltype(fptr, colnum[ii], &datatype, NULL, NULL, status);
+
+      if (datatype <= TLONG && (float) imin == tmin[ii] &&
+                               (float) imax == tmax[ii] &&
+                               (float) ibin == tbin[ii] )
+      {
+        /* This is an integer column and integer limits were entered. */
+        /* Shift the lower and upper histogramming limits by 0.5, so that */
+        /* the values fall in the center of the bin, not on the edge. */
+
+        maxbin[ii] = (taxes[ii] + 1.F);  /* add 1. instead of .5 to avoid roundoff */
+
+        if (tmin[ii] < tmax[ii])
+        {
+          tmin[ii] = tmin[ii] - 0.5F;
+          tmax[ii] = tmax[ii] + 0.5F;
+        }
+        else
+        {
+          tmin[ii] = tmin[ii] + 0.5F;
+          tmax[ii] = tmax[ii] - 0.5F;
+        }
+      } else {  /* not an integer column with integer limits */
+          maxbin[ii] = (tmax[ii] - tmin[ii]) / tbin[ii]; 
+      }
+    }
+
+    /* Set global variables with histogram parameter values.    */
+    /* Use separate scalar variables rather than arrays because */
+    /* it is more efficient when computing the histogram.       */
+
+    histData.hcolnum[0]  = colnum[0];
+    histData.amin1 = tmin[0];
+    histData.maxbin1 = maxbin[0];
+    histData.binsize1 = tbin[0];
+    histData.haxis1 = (long) taxes[0];
+
+    if (histData.haxis > 1)
+    {
+      histData.hcolnum[1]  = colnum[1];
+      histData.amin2 = tmin[1];
+      histData.maxbin2 = maxbin[1];
+      histData.binsize2 = tbin[1];
+      histData.haxis2 = (long) taxes[1];
+
+      if (histData.haxis > 2)
+      {
+        histData.hcolnum[2]  = colnum[2];
+        histData.amin3 = tmin[2];
+        histData.maxbin3 = maxbin[2];
+        histData.binsize3 = tbin[2];
+        histData.haxis3 = (long) taxes[2];
+
+        if (histData.haxis > 3)
+        {
+          histData.hcolnum[3]  = colnum[3];
+          histData.amin4 = tmin[3];
+          histData.maxbin4 = maxbin[3];
+          histData.binsize4 = tbin[3];
+          histData.haxis4 = (long) taxes[3];
+        }
+      }
+    }
+
+    /* define parameters of image for the iterator function */
+    fits_iter_set_file(imagepars, histptr);        /* pointer to image */
+    fits_iter_set_datatype(imagepars, imagetype);  /* image datatype   */
+    fits_iter_set_iotype(imagepars, OutputCol);    /* image is output  */
+
+    /* call the iterator function to write out the histogram image */
+    fits_iterate_data(n_cols, imagepars, offset, n_per_loop,
+                          ffwritehisto, (void*)&histData, status);
+       
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int fits_get_col_minmax(fitsfile *fptr, int colnum, float *datamin, 
+                     float *datamax, int *status)
+/* 
+   Simple utility routine to compute the min and max value in a column
+*/
+{
+    int anynul;
+    long nrows, ntodo, firstrow, ii;
+    float array[1000], nulval;
+
+    ffgky(fptr, TLONG, "NAXIS2", &nrows, NULL, status); /* no. of rows */
+
+    firstrow = 1;
+    nulval = FLOATNULLVALUE;
+    *datamin =  9.0E36F;
+    *datamax = -9.0E36F;
+
+    while(nrows)
+    {
+        ntodo = minvalue(nrows, 100);
+        ffgcv(fptr, TFLOAT, colnum, firstrow, 1, ntodo, &nulval, array,
+              &anynul, status);
+
+        for (ii = 0; ii < ntodo; ii++)
+        {
+            if (array[ii] != nulval)
+            {
+                *datamin = minvalue(*datamin, array[ii]);
+                *datamax = maxvalue(*datamax, array[ii]);
+            }
+        }
+
+        nrows -= ntodo;
+        firstrow += ntodo;
+    }
+    return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int ffwritehisto(long totaln, long pixoffset, long firstn, long nvalues,
+             int narrays, iteratorCol *imagepars, void *userPointer)
+/*
+   Interator work function that writes out the histogram.
+   The histogram values are calculated by another work function, ffcalchisto.
+   This work function only gets called once, and totaln = nvalues.
+*/
+{
+    iteratorCol colpars[5];
+    int ii, status = 0, ncols;
+    long rows_per_loop = 0, offset = 0;
+    histType *histData;
+
+    histData = (histType *)userPointer;
+
+    /* store pointer to the histogram array, and initialize to zero */
+
+    switch( histData->himagetype ) {
+    case TBYTE:
+       histData->hist.b = (char *  ) fits_iter_get_array(imagepars);
+       break;
+    case TSHORT:
+       histData->hist.i = (short * ) fits_iter_get_array(imagepars);
+       break;
+    case TINT:
+       histData->hist.j = (int *   ) fits_iter_get_array(imagepars);
+       break;
+    case TFLOAT:
+       histData->hist.r = (float * ) fits_iter_get_array(imagepars);
+       break;
+    case TDOUBLE:
+       histData->hist.d = (double *) fits_iter_get_array(imagepars);
+       break;
+    }
+
+    /* set the column parameters for the iterator function */
+    for (ii = 0; ii < histData->haxis; ii++)
+    {
+      fits_iter_set_by_num(&colpars[ii], histData->tblptr,
+			   histData->hcolnum[ii], TFLOAT, InputCol);
+    }
+    ncols = histData->haxis;
+
+    if (histData->weight == FLOATNULLVALUE)
+    {
+      fits_iter_set_by_num(&colpars[histData->haxis], histData->tblptr,
+			   histData->wtcolnum, TFLOAT, InputCol);
+      ncols = histData->haxis + 1;
+    }
+
+    /* call iterator function to calc the histogram pixel values */
+
+    /* must lock this call in multithreaded environoments because */
+    /* the ffcalchist work routine uses static vaiables that would */
+    /* get clobbered if multiple threads were running at the same time */
+    FFLOCK;
+    fits_iterate_data(ncols, colpars, offset, rows_per_loop,
+                          ffcalchist, (void*)histData, &status);
+    FFUNLOCK;
+
+    return(status);
+}
+/*--------------------------------------------------------------------------*/
+int ffcalchist(long totalrows, long offset, long firstrow, long nrows,
+             int ncols, iteratorCol *colpars, void *userPointer)
+/*
+   Interator work function that calculates values for the 2D histogram.
+*/
+{
+    long ii, ipix, iaxisbin;
+    float pix, axisbin;
+    static float *col1, *col2, *col3, *col4; /* static to preserve values */
+    static float *wtcol;
+    static long incr2, incr3, incr4;
+    static histType histData;
+    static char *rowselect;
+
+    /*  Initialization procedures: execute on the first call  */
+    if (firstrow == 1)
+    {
+
+      /*  Copy input histogram data to static local variable so we */
+      /*  don't have to constantly dereference it.                 */
+
+      histData = *(histType*)userPointer;
+      rowselect = histData.rowselector;
+
+      /* assign the input array pointers to local pointers */
+      col1 = (float *) fits_iter_get_array(&colpars[0]);
+      if (histData.haxis > 1)
+      {
+        col2 = (float *) fits_iter_get_array(&colpars[1]);
+        incr2 = histData.haxis1;
+
+        if (histData.haxis > 2)
+        {
+          col3 = (float *) fits_iter_get_array(&colpars[2]);
+          incr3 = incr2 * histData.haxis2;
+
+          if (histData.haxis > 3)
+          {
+            col4 = (float *) fits_iter_get_array(&colpars[3]);
+            incr4 = incr3 * histData.haxis3;
+          }
+        }
+      }
+
+      if (ncols > histData.haxis)  /* then weights are give in a column */
+      {
+        wtcol = (float *) fits_iter_get_array(&colpars[histData.haxis]);
+      }
+    }   /* end of Initialization procedures */
+
+    /*  Main loop: increment the histogram at position of each event */
+    for (ii = 1; ii <= nrows; ii++) 
+    {
+        if (rowselect)     /* if a row selector array is supplied... */
+        {
+           if (*rowselect)
+           {
+               rowselect++;   /* this row is included in the histogram */
+           }
+           else
+           {
+               rowselect++;   /* this row is excluded from the histogram */
+               continue;
+           }
+        }
+
+        if (col1[ii] == FLOATNULLVALUE)  /* test for null value */
+            continue;
+
+        pix = (col1[ii] - histData.amin1) / histData.binsize1;
+        ipix = (long) (pix + 1.); /* add 1 because the 1st pixel is the null value */
+
+	/* test if bin is within range */
+        if (ipix < 1 || ipix > histData.haxis1 || pix > histData.maxbin1)
+            continue;
+
+        if (histData.haxis > 1)
+        {
+          if (col2[ii] == FLOATNULLVALUE)
+              continue;
+
+          axisbin = (col2[ii] - histData.amin2) / histData.binsize2;
+          iaxisbin = (long) axisbin;
+
+          if (axisbin < 0. || iaxisbin >= histData.haxis2 || axisbin > histData.maxbin2)
+              continue;
+
+          ipix += (iaxisbin * incr2);
+
+          if (histData.haxis > 2)
+          {
+            if (col3[ii] == FLOATNULLVALUE)
+                continue;
+
+            axisbin = (col3[ii] - histData.amin3) / histData.binsize3;
+            iaxisbin = (long) axisbin;
+            if (axisbin < 0. || iaxisbin >= histData.haxis3 || axisbin > histData.maxbin3)
+                continue;
+
+            ipix += (iaxisbin * incr3);
+ 
+            if (histData.haxis > 3)
+            {
+              if (col4[ii] == FLOATNULLVALUE)
+                  continue;
+
+              axisbin = (col4[ii] - histData.amin4) / histData.binsize4;
+              iaxisbin = (long) axisbin;
+              if (axisbin < 0. || iaxisbin >= histData.haxis4 || axisbin > histData.maxbin4)
+                  continue;
+
+              ipix += (iaxisbin * incr4);
+
+            }  /* end of haxis > 3 case */
+          }    /* end of haxis > 2 case */
+        }      /* end of haxis > 1 case */
+
+        /* increment the histogram pixel */
+        if (histData.weight != FLOATNULLVALUE) /* constant weight factor */
+        {
+            if (histData.himagetype == TINT)
+              histData.hist.j[ipix] += (int) histData.weight;
+            else if (histData.himagetype == TSHORT)
+              histData.hist.i[ipix] += (short) histData.weight;
+            else if (histData.himagetype == TFLOAT)
+              histData.hist.r[ipix] += histData.weight;
+            else if (histData.himagetype == TDOUBLE)
+              histData.hist.d[ipix] += histData.weight;
+            else if (histData.himagetype == TBYTE)
+              histData.hist.b[ipix] += (char) histData.weight;
+        }
+        else if (histData.wtrecip) /* use reciprocal of the weight */
+        {
+            if (histData.himagetype == TINT)
+              histData.hist.j[ipix] += (int) (1./wtcol[ii]);
+            else if (histData.himagetype == TSHORT)
+              histData.hist.i[ipix] += (short) (1./wtcol[ii]);
+            else if (histData.himagetype == TFLOAT)
+              histData.hist.r[ipix] += (float) (1./wtcol[ii]);
+            else if (histData.himagetype == TDOUBLE)
+              histData.hist.d[ipix] += 1./wtcol[ii];
+            else if (histData.himagetype == TBYTE)
+              histData.hist.b[ipix] += (char) (1./wtcol[ii]);
+        }
+        else   /* no weights */
+        {
+            if (histData.himagetype == TINT)
+              histData.hist.j[ipix] += (int) wtcol[ii];
+            else if (histData.himagetype == TSHORT)
+              histData.hist.i[ipix] += (short) wtcol[ii];
+            else if (histData.himagetype == TFLOAT)
+              histData.hist.r[ipix] += wtcol[ii];
+            else if (histData.himagetype == TDOUBLE)
+              histData.hist.d[ipix] += wtcol[ii];
+            else if (histData.himagetype == TBYTE)
+              histData.hist.b[ipix] += (char) wtcol[ii];
+        }
+
+    }  /* end of main loop over all rows */
+
+    return(0);
+}
+