Repatch (from linux) of OSX IRAF

1 files changed, 3888 insertions, 0 deletions

diff --git a/vendor/cfitsio/quantize.c b/vendor/cfitsio/quantize.c
new file mode 100644
index 00000000..8fe6acca
--- /dev/null
+++ b/vendor/cfitsio/quantize.c
@@ -0,0 +1,3888 @@
+/*
+  The following code is based on algorithms written by Richard White at STScI and made
+  available for use in CFITSIO in July 1999 and updated in January 2008. 
+*/
+
+# include <stdio.h>
+# include <stdlib.h>
+# include <math.h>
+# include <limits.h>
+# include <float.h>
+
+#include "fitsio2.h"
+
+/* nearest integer function */
+# define NINT(x)  ((x >= 0.) ? (int) (x + 0.5) : (int) (x - 0.5))
+
+#define NULL_VALUE -2147483647 /* value used to represent undefined pixels */
+#define N_RESERVED_VALUES 10   /* number of reserved values, starting with */
+                               /* and including NULL_VALUE.  These values */
+                               /* may not be used to represent the quantized */
+                               /* and scaled floating point pixel values */
+			       /* If lossy Hcompression is used, and the */
+			       /* array contains null values, then it is also */
+			       /* possible for the compressed values to slightly */
+			       /* exceed the range of the actual (lossless) values */
+			       /* so we must reserve a little more space */
+			       
+/* more than this many standard deviations from the mean is an outlier */
+# define SIGMA_CLIP     5.
+# define NITER          3	/* number of sigma-clipping iterations */
+
+static int FnMeanSigma_short(short *array, long npix, int nullcheck, 
+  short nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);       
+static int FnMeanSigma_int(int *array, long npix, int nullcheck,
+  int nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);       
+static int FnMeanSigma_float(float *array, long npix, int nullcheck,
+  float nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);       
+static int FnMeanSigma_double(double *array, long npix, int nullcheck,
+  double nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);       
+
+static int FnNoise5_short(short *array, long nx, long ny, int nullcheck, 
+   short nullvalue, long *ngood, short *minval, short *maxval, 
+   double *n2, double *n3, double *n5, int *status);   
+static int FnNoise5_int(int *array, long nx, long ny, int nullcheck, 
+   int nullvalue, long *ngood, int *minval, int *maxval, 
+   double *n2, double *n3, double *n5, int *status);   
+static int FnNoise5_float(float *array, long nx, long ny, int nullcheck, 
+   float nullvalue, long *ngood, float *minval, float *maxval, 
+   double *n2, double *n3, double *n5, int *status);   
+static int FnNoise5_double(double *array, long nx, long ny, int nullcheck, 
+   double nullvalue, long *ngood, double *minval, double *maxval, 
+   double *n2, double *n3, double *n5, int *status);   
+
+static int FnNoise3_short(short *array, long nx, long ny, int nullcheck, 
+   short nullvalue, long *ngood, short *minval, short *maxval, double *noise, int *status);       
+static int FnNoise3_int(int *array, long nx, long ny, int nullcheck, 
+   int nullvalue, long *ngood, int *minval, int *maxval, double *noise, int *status);          
+static int FnNoise3_float(float *array, long nx, long ny, int nullcheck, 
+   float nullvalue, long *ngood, float *minval, float *maxval, double *noise, int *status);        
+static int FnNoise3_double(double *array, long nx, long ny, int nullcheck, 
+   double nullvalue, long *ngood, double *minval, double *maxval, double *noise, int *status);        
+
+static int FnNoise1_short(short *array, long nx, long ny, 
+   int nullcheck, short nullvalue, double *noise, int *status);       
+static int FnNoise1_int(int *array, long nx, long ny, 
+   int nullcheck, int nullvalue, double *noise, int *status);       
+static int FnNoise1_float(float *array, long nx, long ny, 
+   int nullcheck, float nullvalue, double *noise, int *status);       
+static int FnNoise1_double(double *array, long nx, long ny, 
+   int nullcheck, double nullvalue, double *noise, int *status);       
+
+static int FnCompare_short (const void *, const void *);
+static int FnCompare_int (const void *, const void *);
+static int FnCompare_float (const void *, const void *);
+static int FnCompare_double (const void *, const void *);
+static float quick_select_float(float arr[], int n);
+static short quick_select_short(short arr[], int n);
+static int quick_select_int(int arr[], int n);
+static LONGLONG quick_select_longlong(LONGLONG arr[], int n);
+static double quick_select_double(double arr[], int n);
+
+/*---------------------------------------------------------------------------*/
+int fits_quantize_float (long row, float fdata[], long nxpix, long nypix, int nullcheck, 
+	float in_null_value, float qlevel, int idata[], double *bscale,
+	double *bzero, int *iminval, int *imaxval) {
+
+/* arguments:
+long row            i: if positive, tile number = row number in the binary table
+                       (this is only used when dithering the quantized values)
+float fdata[]       i: array of image pixels to be compressed
+long nxpix          i: number of pixels in each row of fdata
+long nypix          i: number of rows in fdata
+nullcheck           i: check for nullvalues in fdata?
+float in_null_value i: value used to represent undefined pixels in fdata
+float qlevel        i: quantization level
+int idata[]         o: values of fdata after applying bzero and bscale
+double bscale       o: scale factor
+double bzero        o: zero offset
+int iminval         o: minimum quantized value that is returned
+int imaxval         o: maximum quantized value that is returned
+
+The function value will be one if the input fdata were copied to idata;
+in this case the parameters bscale and bzero can be used to convert back to
+nearly the original floating point values:  fdata ~= idata * bscale + bzero.
+If the function value is zero, the data were not copied to idata.
+*/
+
+	int status, anynulls = 0, iseed;
+	long i, nx, ngood = 0;
+	double stdev, noise2, noise3, noise5;	/* MAD 2nd, 3rd, and 5th order noise values */
+	float minval = 0., maxval = 0.;  /* min & max of fdata */
+	double delta;		/* bscale, 1 in idata = delta in fdata */
+	double zeropt;	        /* bzero */
+	double temp;
+        int nextrand = 0;
+	extern float *fits_rand_value; /* this is defined in imcompress.c */
+	LONGLONG iqfactor;
+
+	nx = nxpix * nypix;
+	if (nx <= 1) {
+	    *bscale = 1.;
+	    *bzero  = 0.;
+	    return (0);
+	}
+
+        if (qlevel >= 0.) {
+
+	    /* estimate background noise using MAD pixel differences */
+	    FnNoise5_float(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
+	        &minval, &maxval, &noise2, &noise3, &noise5, &status);      
+
+	    if (nullcheck && ngood == 0) {   /* special case of an image filled with Nulls */
+	        /* set parameters to dummy values, which are not used */
+		minval = 0.;
+		maxval = 1.;
+		stdev = 1;
+	    } else {
+
+	        /* use the minimum of noise2, noise3, and noise5 as the best noise value */
+	        stdev = noise3;
+	        if (noise2 != 0. && noise2 < stdev) stdev = noise2;
+	        if (noise5 != 0. && noise5 < stdev) stdev = noise5;
+            }
+
+	    if (qlevel == 0.)
+	        delta = stdev / 4.;  /* default quantization */
+	    else
+	        delta = stdev / qlevel;
+
+	    if (delta == 0.) 
+	        return (0);			/* don't quantize */
+
+	} else {
+	    /* negative value represents the absolute quantization level */
+	    delta = -qlevel;
+
+	    /* only nned to calculate the min and max values */
+	    FnNoise3_float(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
+	        &minval, &maxval, 0, &status);      
+ 	}
+
+        /* check that the range of quantized levels is not > range of int */
+	if ((maxval - minval) / delta > 2. * 2147483647. - N_RESERVED_VALUES )
+	    return (0);			/* don't quantize */
+
+        if (row > 0) { /* we need to dither the quantized values */
+            if (!fits_rand_value) 
+	        if (fits_init_randoms()) return(MEMORY_ALLOCATION);
+
+	    /* initialize the index to the next random number in the list */
+            iseed = (int) ((row - 1) % N_RANDOM);
+	    nextrand = (int) (fits_rand_value[iseed] * 500.);
+	}
+
+        if (ngood == nx) {   /* don't have to check for nulls */
+            /* return all positive values, if possible since some */
+            /* compression algorithms either only work for positive integers, */
+            /* or are more efficient.  */
+
+            if ((maxval - minval) / delta < 2147483647. - N_RESERVED_VALUES )
+            {
+                zeropt = minval;
+		/* fudge the zero point so it is an integer multiple of delta */
+		/* This helps to ensure the same scaling will be performed if the */
+		/* file undergoes multiple fpack/funpack cycles */
+		iqfactor = (LONGLONG) (zeropt/delta  + 0.5);
+		zeropt = iqfactor * delta;               
+            }
+            else
+            {
+                /* center the quantized levels around zero */
+                zeropt = (minval + maxval) / 2.;
+            }
+
+            if (row > 0) {  /* dither the values when quantizing */
+              for (i = 0;  i < nx;  i++) {
+	    
+		idata[i] =  NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
+
+                nextrand++;
+		if (nextrand == N_RANDOM) {
+                    iseed++;
+		    if (iseed == N_RANDOM) iseed = 0;
+	            nextrand = (int) (fits_rand_value[iseed] * 500);
+                }
+              }
+            } else {  /* do not dither the values */
+
+       	        for (i = 0;  i < nx;  i++) {
+	            idata[i] = NINT ((fdata[i] - zeropt) / delta);
+                }
+            } 
+        }
+        else {
+            /* data contains null values; shift the range to be */
+            /* close to the value used to represent null values */
+            zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
+
+            if (row > 0) {  /* dither the values */
+	      for (i = 0;  i < nx;  i++) {
+                if (fdata[i] != in_null_value) {
+		    idata[i] =  NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
+                } else {
+                    idata[i] = NULL_VALUE;
+                }
+
+                /* increment the random number index, regardless */
+                nextrand++;
+		if (nextrand == N_RANDOM) {
+                      iseed++;
+		      if (iseed == N_RANDOM) iseed = 0;
+	              nextrand = (int) (fits_rand_value[iseed] * 500);
+                }
+              }
+            } else {  /* do not dither the values */
+	       for (i = 0;  i < nx;  i++) {
+                 if (fdata[i] != in_null_value)
+		    idata[i] =  NINT((fdata[i] - zeropt) / delta);
+                 else 
+                    idata[i] = NULL_VALUE;
+               }
+            }
+	}
+
+        /* calc min and max values */
+        temp = (minval - zeropt) / delta;
+        *iminval =  NINT (temp);
+        temp = (maxval - zeropt) / delta;
+        *imaxval =  NINT (temp);
+
+	*bscale = delta;
+	*bzero = zeropt;
+	return (1);			/* yes, data have been quantized */
+}
+/*---------------------------------------------------------------------------*/
+int fits_quantize_double (long row, double fdata[], long nxpix, long nypix, int nullcheck, 
+	double in_null_value, float qlevel, int idata[], double *bscale,
+	double *bzero, int *iminval, int *imaxval) {
+
+/* arguments:
+long row            i: tile number = row number in the binary table
+double fdata[]      i: array of image pixels to be compressed
+long nxpix          i: number of pixels in each row of fdata
+long nypix          i: number of rows in fdata
+nullcheck           i: check for nullvalues in fdata?
+double in_null_value i: value used to represent undefined pixels in fdata
+int noise_bits      i: quantization level (number of bits)
+int idata[]         o: values of fdata after applying bzero and bscale
+double bscale       o: scale factor
+double bzero        o: zero offset
+int iminval         o: minimum quantized value that is returned
+int imaxval         o: maximum quantized value that is returned
+
+The function value will be one if the input fdata were copied to idata;
+in this case the parameters bscale and bzero can be used to convert back to
+nearly the original floating point values:  fdata ~= idata * bscale + bzero.
+If the function value is zero, the data were not copied to idata.
+*/
+
+	int status, anynulls = 0, iseed;
+	long i, nx, ngood = 0;
+	double stdev, noise2, noise3, noise5;	/* MAD 2nd, 3rd, and 5th order noise values */
+	double minval = 0., maxval = 0.;  /* min & max of fdata */
+	double delta;		/* bscale, 1 in idata = delta in fdata */
+	double zeropt;	        /* bzero */
+	double temp;
+        int nextrand = 0;
+	extern float *fits_rand_value;
+	LONGLONG iqfactor;
+
+	nx = nxpix * nypix;
+	if (nx <= 1) {
+	    *bscale = 1.;
+	    *bzero  = 0.;
+	    return (0);
+	}
+
+        if (qlevel >= 0.) {
+
+	    /* estimate background noise using MAD pixel differences */
+	    FnNoise5_double(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
+	        &minval, &maxval, &noise2, &noise3, &noise5, &status);      
+
+	    if (nullcheck && ngood == 0) {   /* special case of an image filled with Nulls */
+	        /* set parameters to dummy values, which are not used */
+		minval = 0.;
+		maxval = 1.;
+		stdev = 1;
+	    } else {
+
+	        /* use the minimum of noise2, noise3, and noise5 as the best noise value */
+	        stdev = noise3;
+	        if (noise2 != 0. && noise2 < stdev) stdev = noise2;
+	        if (noise5 != 0. && noise5 < stdev) stdev = noise5;
+            }
+
+	    if (qlevel == 0.)
+	        delta = stdev / 4.;  /* default quantization */
+	    else
+	        delta = stdev / qlevel;
+
+	    if (delta == 0.) 
+	        return (0);			/* don't quantize */
+
+	} else {
+	    /* negative value represents the absolute quantization level */
+	    delta = -qlevel;
+
+	    /* only nned to calculate the min and max values */
+	    FnNoise3_double(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
+	        &minval, &maxval, 0, &status);      
+ 	}
+
+        /* check that the range of quantized levels is not > range of int */
+	if ((maxval - minval) / delta > 2. * 2147483647. - N_RESERVED_VALUES )
+	    return (0);			/* don't quantize */
+
+        if (row > 0) { /* we need to dither the quantized values */
+            if (!fits_rand_value) 
+	       if (fits_init_randoms()) return(MEMORY_ALLOCATION);
+
+	    /* initialize the index to the next random number in the list */
+            iseed = (int) ((row - 1) % N_RANDOM);
+	    nextrand = (int) (fits_rand_value[iseed] * 500);
+	}
+
+        if (ngood == nx) {   /* don't have to check for nulls */
+            /* return all positive values, if possible since some */
+            /* compression algorithms either only work for positive integers, */
+            /* or are more efficient.  */
+            if ((maxval - minval) / delta < 2147483647. - N_RESERVED_VALUES )
+            {
+                zeropt = minval;
+		/* fudge the zero point so it is an integer multiple of delta */
+		/* This helps to ensure the same scaling will be performed if the */
+		/* file undergoes multiple fpack/funpack cycles */
+		iqfactor = (LONGLONG) (zeropt/delta  + 0.5);
+		zeropt = iqfactor * delta;               
+            }
+            else
+            {
+                /* center the quantized levels around zero */
+                zeropt = (minval + maxval) / 2.;
+            }
+
+            if (row > 0) {  /* dither the values when quantizing */
+       	      for (i = 0;  i < nx;  i++) {
+
+		idata[i] =  NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
+
+                nextrand++;
+		if (nextrand == N_RANDOM) {
+                    iseed++;
+	            nextrand = (int) (fits_rand_value[iseed] * 500);
+                }
+              }
+            } else {  /* do not dither the values */
+
+       	        for (i = 0;  i < nx;  i++) {
+	            idata[i] = NINT ((fdata[i] - zeropt) / delta);
+                }
+            } 
+        }
+        else {
+            /* data contains null values; shift the range to be */
+            /* close to the value used to represent null values */
+            zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
+
+            if (row > 0) {  /* dither the values */
+	      for (i = 0;  i < nx;  i++) {
+                if (fdata[i] != in_null_value) {
+		    idata[i] =  NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
+                } else {
+                    idata[i] = NULL_VALUE;
+                }
+
+                /* increment the random number index, regardless */
+                nextrand++;
+		if (nextrand == N_RANDOM) {
+                        iseed++;
+	                nextrand = (int) (fits_rand_value[iseed] * 500);
+                } 
+              }
+            } else {  /* do not dither the values */
+	       for (i = 0;  i < nx;  i++) {
+                 if (fdata[i] != in_null_value)
+		    idata[i] =  NINT((fdata[i] - zeropt) / delta);
+                 else 
+                    idata[i] = NULL_VALUE;
+               }
+            }
+	}
+
+        /* calc min and max values */
+        temp = (minval - zeropt) / delta;
+        *iminval =  NINT (temp);
+        temp = (maxval - zeropt) / delta;
+        *imaxval =  NINT (temp);
+
+	*bscale = delta;
+	*bzero = zeropt;
+
+	return (1);			/* yes, data have been quantized */
+}
+/*--------------------------------------------------------------------------*/
+int fits_img_stats_short(short *array, /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+	long ny,            /* number of rows in the image */
+	                    /* (if this is a 3D image, then ny should be the */
+			    /* product of the no. of rows times the no. of planes) */
+	int nullcheck,      /* check for null values, if true */
+	short nullvalue,    /* value of null pixels, if nullcheck is true */
+
+   /* returned parameters (if the pointer is not null)  */
+	long *ngoodpix,     /* number of non-null pixels in the image */
+	short *minvalue,    /* returned minimum non-null value in the array */
+	short *maxvalue,    /* returned maximum non-null value in the array */
+	double *mean,       /* returned mean value of all non-null pixels */
+	double *sigma,      /* returned R.M.S. value of all non-null pixels */
+	double *noise1,     /* 1st order estimate of noise in image background level */
+	double *noise2,     /* 2nd order estimate of noise in image background level */
+	double *noise3,     /* 3rd order estimate of noise in image background level */
+	double *noise5,     /* 5th order estimate of noise in image background level */
+	int *status)        /* error status */
+
+/*
+    Compute statistics of the input short integer image.
+*/
+{
+	long ngood;
+	short minval, maxval;
+	double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.;
+
+	/* need to calculate mean and/or sigma and/or limits? */
+	if (mean || sigma ) {
+		FnMeanSigma_short(array, nx * ny, nullcheck, nullvalue, 
+			&ngood, &xmean, &xsigma, status);
+
+	    if (ngoodpix) *ngoodpix = ngood;
+	    if (mean)     *mean = xmean;
+	    if (sigma)    *sigma = xsigma;
+	}
+
+	if (noise1) {
+		FnNoise1_short(array, nx, ny, nullcheck, nullvalue, 
+		  &xnoise, status);
+
+		*noise1  = xnoise;
+	}
+
+	if (minvalue || maxvalue || noise3) {
+		FnNoise5_short(array, nx, ny, nullcheck, nullvalue, 
+			&ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status);
+
+		if (ngoodpix) *ngoodpix = ngood;
+		if (minvalue) *minvalue= minval;
+		if (maxvalue) *maxvalue = maxval;
+		if (noise2) *noise2  = xnoise2;
+		if (noise3) *noise3  = xnoise3;
+		if (noise5) *noise5  = xnoise5;
+	}
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int fits_img_stats_int(int *array, /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+	long ny,            /* number of rows in the image */
+	                    /* (if this is a 3D image, then ny should be the */
+			    /* product of the no. of rows times the no. of planes) */
+	int nullcheck,      /* check for null values, if true */
+	int nullvalue,    /* value of null pixels, if nullcheck is true */
+
+   /* returned parameters (if the pointer is not null)  */
+	long *ngoodpix,     /* number of non-null pixels in the image */
+	int *minvalue,    /* returned minimum non-null value in the array */
+	int *maxvalue,    /* returned maximum non-null value in the array */
+	double *mean,       /* returned mean value of all non-null pixels */
+	double *sigma,      /* returned R.M.S. value of all non-null pixels */
+	double *noise1,     /* 1st order estimate of noise in image background level */
+	double *noise2,     /* 2nd order estimate of noise in image background level */
+	double *noise3,     /* 3rd order estimate of noise in image background level */
+	double *noise5,     /* 5th order estimate of noise in image background level */
+	int *status)        /* error status */
+
+/*
+    Compute statistics of the input integer image.
+*/
+{
+	long ngood;
+	int minval, maxval;
+	double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.;
+
+	/* need to calculate mean and/or sigma and/or limits? */
+	if (mean || sigma ) {
+		FnMeanSigma_int(array, nx * ny, nullcheck, nullvalue, 
+			&ngood, &xmean, &xsigma, status);
+
+	    if (ngoodpix) *ngoodpix = ngood;
+	    if (mean)     *mean = xmean;
+	    if (sigma)    *sigma = xsigma;
+	}
+
+	if (noise1) {
+		FnNoise1_int(array, nx, ny, nullcheck, nullvalue, 
+		  &xnoise, status);
+
+		*noise1  = xnoise;
+	}
+
+	if (minvalue || maxvalue || noise3) {
+		FnNoise5_int(array, nx, ny, nullcheck, nullvalue, 
+			&ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status);
+
+		if (ngoodpix) *ngoodpix = ngood;
+		if (minvalue) *minvalue= minval;
+		if (maxvalue) *maxvalue = maxval;
+		if (noise2) *noise2  = xnoise2;
+		if (noise3) *noise3  = xnoise3;
+		if (noise5) *noise5  = xnoise5;
+	}
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+int fits_img_stats_float(float *array, /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+	long ny,            /* number of rows in the image */
+	                    /* (if this is a 3D image, then ny should be the */
+			    /* product of the no. of rows times the no. of planes) */
+	int nullcheck,      /* check for null values, if true */
+	float nullvalue,    /* value of null pixels, if nullcheck is true */
+
+   /* returned parameters (if the pointer is not null)  */
+	long *ngoodpix,     /* number of non-null pixels in the image */
+	float *minvalue,    /* returned minimum non-null value in the array */
+	float *maxvalue,    /* returned maximum non-null value in the array */
+	double *mean,       /* returned mean value of all non-null pixels */
+	double *sigma,      /* returned R.M.S. value of all non-null pixels */
+	double *noise1,     /* 1st order estimate of noise in image background level */
+	double *noise2,     /* 2nd order estimate of noise in image background level */
+	double *noise3,     /* 3rd order estimate of noise in image background level */
+	double *noise5,     /* 5th order estimate of noise in image background level */
+	int *status)        /* error status */
+
+/*
+    Compute statistics of the input float image.
+*/
+{
+	long ngood;
+	float minval, maxval;
+	double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.;
+
+	/* need to calculate mean and/or sigma and/or limits? */
+	if (mean || sigma ) {
+		FnMeanSigma_float(array, nx * ny, nullcheck, nullvalue, 
+			&ngood, &xmean, &xsigma, status);
+
+	    if (ngoodpix) *ngoodpix = ngood;
+	    if (mean)     *mean = xmean;
+	    if (sigma)    *sigma = xsigma;
+	}
+
+	if (noise1) {
+		FnNoise1_float(array, nx, ny, nullcheck, nullvalue, 
+		  &xnoise, status);
+
+		*noise1  = xnoise;
+	}
+
+	if (minvalue || maxvalue || noise3) {
+		FnNoise5_float(array, nx, ny, nullcheck, nullvalue, 
+			&ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status);
+
+		if (ngoodpix) *ngoodpix = ngood;
+		if (minvalue) *minvalue= minval;
+		if (maxvalue) *maxvalue = maxval;
+		if (noise2) *noise2  = xnoise2;
+		if (noise3) *noise3  = xnoise3;
+		if (noise5) *noise5  = xnoise5;
+	}
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnMeanSigma_short
+       (short *array,       /*  2 dimensional array of image pixels */
+        long npix,          /* number of pixels in the image */
+	int nullcheck,      /* check for null values, if true */
+	short nullvalue,    /* value of null pixels, if nullcheck is true */
+
+   /* returned parameters */
+   
+	long *ngoodpix,     /* number of non-null pixels in the image */
+	double *mean,       /* returned mean value of all non-null pixels */
+	double *sigma,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Compute mean and RMS sigma of the non-null pixels in the input array.
+*/
+{
+	long ii, ngood = 0;
+	short *value;
+	double sum = 0., sum2 = 0., xtemp;
+
+	value = array;
+	    
+	if (nullcheck) {
+	        for (ii = 0; ii < npix; ii++, value++) {
+		    if (*value != nullvalue) {
+		        ngood++;
+		        xtemp = (double) *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		    }
+		}
+	} else {
+	        ngood = npix;
+	        for (ii = 0; ii < npix; ii++, value++) {
+		        xtemp = (double) *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		}
+	}
+
+	if (ngood > 1) {
+		if (ngoodpix) *ngoodpix = ngood;
+		xtemp = sum / ngood;
+		if (mean)     *mean = xtemp;
+		if (sigma)    *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
+	} else if (ngood == 1){
+		if (ngoodpix) *ngoodpix = 1;
+		if (mean)     *mean = sum;
+		if (sigma)    *sigma = 0.0;
+	} else {
+		if (ngoodpix) *ngoodpix = 0;
+	        if (mean)     *mean = 0.;
+		if (sigma)    *sigma = 0.;
+	}	    
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnMeanSigma_int
+       (int *array,       /*  2 dimensional array of image pixels */
+        long npix,          /* number of pixels in the image */
+	int nullcheck,      /* check for null values, if true */
+	int nullvalue,    /* value of null pixels, if nullcheck is true */
+
+   /* returned parameters */
+   
+	long *ngoodpix,     /* number of non-null pixels in the image */
+	double *mean,       /* returned mean value of all non-null pixels */
+	double *sigma,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Compute mean and RMS sigma of the non-null pixels in the input array.
+*/
+{
+	long ii, ngood = 0;
+	int *value;
+	double sum = 0., sum2 = 0., xtemp;
+
+	value = array;
+	    
+	if (nullcheck) {
+	        for (ii = 0; ii < npix; ii++, value++) {
+		    if (*value != nullvalue) {
+		        ngood++;
+		        xtemp = (double) *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		    }
+		}
+	} else {
+	        ngood = npix;
+	        for (ii = 0; ii < npix; ii++, value++) {
+		        xtemp = (double) *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		}
+	}
+
+	if (ngood > 1) {
+		if (ngoodpix) *ngoodpix = ngood;
+		xtemp = sum / ngood;
+		if (mean)     *mean = xtemp;
+		if (sigma)    *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
+	} else if (ngood == 1){
+		if (ngoodpix) *ngoodpix = 1;
+		if (mean)     *mean = sum;
+		if (sigma)    *sigma = 0.0;
+	} else {
+		if (ngoodpix) *ngoodpix = 0;
+	        if (mean)     *mean = 0.;
+		if (sigma)    *sigma = 0.;
+	}	    
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnMeanSigma_float
+       (float *array,       /*  2 dimensional array of image pixels */
+        long npix,          /* number of pixels in the image */
+	int nullcheck,      /* check for null values, if true */
+	float nullvalue,    /* value of null pixels, if nullcheck is true */
+
+   /* returned parameters */
+   
+	long *ngoodpix,     /* number of non-null pixels in the image */
+	double *mean,       /* returned mean value of all non-null pixels */
+	double *sigma,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Compute mean and RMS sigma of the non-null pixels in the input array.
+*/
+{
+	long ii, ngood = 0;
+	float *value;
+	double sum = 0., sum2 = 0., xtemp;
+
+	value = array;
+	    
+	if (nullcheck) {
+	        for (ii = 0; ii < npix; ii++, value++) {
+		    if (*value != nullvalue) {
+		        ngood++;
+		        xtemp = (double) *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		    }
+		}
+	} else {
+	        ngood = npix;
+	        for (ii = 0; ii < npix; ii++, value++) {
+		        xtemp = (double) *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		}
+	}
+
+	if (ngood > 1) {
+		if (ngoodpix) *ngoodpix = ngood;
+		xtemp = sum / ngood;
+		if (mean)     *mean = xtemp;
+		if (sigma)    *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
+	} else if (ngood == 1){
+		if (ngoodpix) *ngoodpix = 1;
+		if (mean)     *mean = sum;
+		if (sigma)    *sigma = 0.0;
+	} else {
+		if (ngoodpix) *ngoodpix = 0;
+	        if (mean)     *mean = 0.;
+		if (sigma)    *sigma = 0.;
+	}	    
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnMeanSigma_double
+       (double *array,       /*  2 dimensional array of image pixels */
+        long npix,          /* number of pixels in the image */
+	int nullcheck,      /* check for null values, if true */
+	double nullvalue,    /* value of null pixels, if nullcheck is true */
+
+   /* returned parameters */
+   
+	long *ngoodpix,     /* number of non-null pixels in the image */
+	double *mean,       /* returned mean value of all non-null pixels */
+	double *sigma,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Compute mean and RMS sigma of the non-null pixels in the input array.
+*/
+{
+	long ii, ngood = 0;
+	double *value;
+	double sum = 0., sum2 = 0., xtemp;
+
+	value = array;
+	    
+	if (nullcheck) {
+	        for (ii = 0; ii < npix; ii++, value++) {
+		    if (*value != nullvalue) {
+		        ngood++;
+		        xtemp = *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		    }
+		}
+	} else {
+	        ngood = npix;
+	        for (ii = 0; ii < npix; ii++, value++) {
+		        xtemp = *value;
+		        sum += xtemp;
+		        sum2 += (xtemp * xtemp);
+		}
+	}
+
+	if (ngood > 1) {
+		if (ngoodpix) *ngoodpix = ngood;
+		xtemp = sum / ngood;
+		if (mean)     *mean = xtemp;
+		if (sigma)    *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
+	} else if (ngood == 1){
+		if (ngoodpix) *ngoodpix = 1;
+		if (mean)     *mean = sum;
+		if (sigma)    *sigma = 0.0;
+	} else {
+		if (ngoodpix) *ngoodpix = 0;
+	        if (mean)     *mean = 0.;
+		if (sigma)    *sigma = 0.;
+	}	    
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise5_short
+       (short *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	short nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	short *minval,    /* minimum non-null value */
+	short *maxval,    /* maximum non-null value */
+	double *noise2,      /* returned 2nd order MAD of all non-null pixels */
+	double *noise3,      /* returned 3rd order MAD of all non-null pixels */
+	double *noise5,      /* returned 5th order MAD of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the median and background noise in the input image using 2nd, 3rd and 5th
+order Median Absolute Differences.
+
+The noise in the background of the image is calculated using the MAD algorithms 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+3rd order:  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
+	int *differences2, *differences3, *differences5;
+	short *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
+	short xminval = SHRT_MAX, xmaxval = SHRT_MIN;
+	int do_range = 0;
+	double *diffs2, *diffs3, *diffs5; 
+	double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
+	
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 9 pixels */
+	if (nx < 9) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise2) *noise2 = 0.;
+		if (noise3) *noise3 = 0.;
+		if (noise5) *noise5 = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences2 = calloc(nx, sizeof(int));
+	if (!differences2) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences3 = calloc(nx, sizeof(int));
+	if (!differences3) {
+		free(differences2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences5 = calloc(nx, sizeof(int));
+	if (!differences5) {
+		free(differences2);
+		free(differences3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs2 = calloc(ny, sizeof(double));
+	if (!diffs2) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs3 = calloc(ny, sizeof(double));
+	if (!diffs3) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs5 = calloc(ny, sizeof(double));
+	if (!diffs5) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+		free(diffs3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+			
+		/* find the 5th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v5 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		}
+				
+		/* find the 6th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v6 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v6 < xminval) xminval = v6;
+			if (v6 > xmaxval) xmaxval = v6;
+		}
+				
+		/* find the 7th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v7 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v7 < xminval) xminval = v7;
+			if (v7 > xmaxval) xmaxval = v7;
+		}
+				
+		/* find the 8th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v8 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v8 < xminval) xminval = v8;
+			if (v8 > xmaxval) xmaxval = v8;
+		}
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		nvals2 = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		    v9 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v9 < xminval) xminval = v9;
+			if (v9 > xmaxval) xmaxval = v9;
+		    }
+
+		    /* construct array of absolute differences */
+
+		    if (!(v5 == v6 && v6 == v7) ) {
+		        differences2[nvals2] =  abs((int) v5 - (int) v7);
+			nvals2++;
+		    }
+
+		    if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
+		        differences3[nvals] =  abs((2 * (int) v5) - (int) v3 - (int) v7);
+		        differences5[nvals] =  abs((6 * (int) v5) - (4 * (int) v3) - (4 * (int) v7) + (int) v1 + (int) v9);
+		        nvals++;  
+		    } else {
+		        /* ignore constant background regions */
+			ngoodpix++;
+		    }
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+		    v5 = v6;
+		    v6 = v7;
+		    v7 = v8;
+		    v8 = v9;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the median diffs */
+		/* Note that there are 8 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 8);
+
+		if (nvals == 0) {
+		    continue;  /* cannot compute medians on this row */
+		} else if (nvals == 1) {
+		    if (nvals2 == 1) {
+		        diffs2[nrows2] = differences2[0];
+			nrows2++;
+		    }
+		        
+		    diffs3[nrows] = differences3[0];
+		    diffs5[nrows] = differences5[0];
+		} else {
+                    /* quick_select returns the median MUCH faster than using qsort */
+		    if (nvals2 > 1) {
+                        diffs2[nrows2] = quick_select_int(differences2, nvals);
+			nrows2++;
+		    }
+
+                    diffs3[nrows] = quick_select_int(differences3, nvals);
+                    diffs5[nrows] = quick_select_int(differences5, nvals);
+		}
+
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise3 = 0;
+	       xnoise5 = 0;
+	} else if (nrows == 1) {
+	       xnoise3 = diffs3[0];
+	       xnoise5 = diffs5[0];
+	} else {	    
+	       qsort(diffs3, nrows, sizeof(double), FnCompare_double);
+	       qsort(diffs5, nrows, sizeof(double), FnCompare_double);
+	       xnoise3 =  (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
+	       xnoise5 =  (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
+	}
+
+	if (nrows2 == 0) { 
+	       xnoise2 = 0;
+	} else if (nrows2 == 1) {
+	       xnoise2 = diffs2[0];
+	} else {	    
+	       qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
+	       xnoise2 =  (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise2)  *noise2  = 1.0483579 * xnoise2;
+	if (noise3)  *noise3  = 0.6052697 * xnoise3;
+	if (noise5)  *noise5  = 0.1772048 * xnoise5;
+
+	free(diffs5);
+	free(diffs3);
+	free(diffs2);
+	free(differences5);
+	free(differences3);
+	free(differences2);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise5_int
+       (int *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	int nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	int *minval,    /* minimum non-null value */
+	int *maxval,    /* maximum non-null value */
+	double *noise2,      /* returned 2nd order MAD of all non-null pixels */
+	double *noise3,      /* returned 3rd order MAD of all non-null pixels */
+	double *noise5,      /* returned 5th order MAD of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the median and background noise in the input image using 2nd, 3rd and 5th
+order Median Absolute Differences.
+
+The noise in the background of the image is calculated using the MAD algorithms 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+3rd order:  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
+	LONGLONG *differences2, *differences3, *differences5, tdiff;
+	int *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
+	int xminval = INT_MAX, xmaxval = INT_MIN;
+	int do_range = 0;
+	double *diffs2, *diffs3, *diffs5; 
+	double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
+	
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 9 pixels */
+	if (nx < 9) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise2) *noise2 = 0.;
+		if (noise3) *noise3 = 0.;
+		if (noise5) *noise5 = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences2 = calloc(nx, sizeof(LONGLONG));
+	if (!differences2) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences3 = calloc(nx, sizeof(LONGLONG));
+	if (!differences3) {
+		free(differences2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences5 = calloc(nx, sizeof(LONGLONG));
+	if (!differences5) {
+		free(differences2);
+		free(differences3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs2 = calloc(ny, sizeof(double));
+	if (!diffs2) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs3 = calloc(ny, sizeof(double));
+	if (!diffs3) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs5 = calloc(ny, sizeof(double));
+	if (!diffs5) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+		free(diffs3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+			
+		/* find the 5th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v5 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		}
+				
+		/* find the 6th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v6 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v6 < xminval) xminval = v6;
+			if (v6 > xmaxval) xmaxval = v6;
+		}
+				
+		/* find the 7th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v7 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v7 < xminval) xminval = v7;
+			if (v7 > xmaxval) xmaxval = v7;
+		}
+				
+		/* find the 8th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v8 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v8 < xminval) xminval = v8;
+			if (v8 > xmaxval) xmaxval = v8;
+		}
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		nvals2 = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		    v9 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v9 < xminval) xminval = v9;
+			if (v9 > xmaxval) xmaxval = v9;
+		    }
+
+		    /* construct array of absolute differences */
+
+		    if (!(v5 == v6 && v6 == v7) ) {
+		        tdiff =  (LONGLONG) v5 - (LONGLONG) v7;
+			if (tdiff < 0)
+		            differences2[nvals2] =  -1 * tdiff;
+			else
+		            differences2[nvals2] =  tdiff;
+
+			nvals2++;
+		    }
+
+		    if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
+		        tdiff =  (2 * (LONGLONG) v5) - (LONGLONG) v3 - (LONGLONG) v7;
+			if (tdiff < 0)
+		            differences3[nvals] =  -1 * tdiff;
+			else
+		            differences3[nvals] =  tdiff;
+
+		        tdiff =  (6 * (LONGLONG) v5) - (4 * (LONGLONG) v3) - (4 * (LONGLONG) v7) + (LONGLONG) v1 + (LONGLONG) v9;
+			if (tdiff < 0)
+		            differences5[nvals] =  -1 * tdiff;
+			else
+		            differences5[nvals] =  tdiff;
+
+		        nvals++;  
+		    } else {
+		        /* ignore constant background regions */
+			ngoodpix++;
+		    }
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+		    v5 = v6;
+		    v6 = v7;
+		    v7 = v8;
+		    v8 = v9;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the median diffs */
+		/* Note that there are 8 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 8);
+
+		if (nvals == 0) {
+		    continue;  /* cannot compute medians on this row */
+		} else if (nvals == 1) {
+		    if (nvals2 == 1) {
+		        diffs2[nrows2] = (double) differences2[0];
+			nrows2++;
+		    }
+		        
+		    diffs3[nrows] = (double) differences3[0];
+		    diffs5[nrows] = (double) differences5[0];
+		} else {
+                    /* quick_select returns the median MUCH faster than using qsort */
+		    if (nvals2 > 1) {
+                        diffs2[nrows2] = (double) quick_select_longlong(differences2, nvals);
+			nrows2++;
+		    }
+
+                    diffs3[nrows] = (double) quick_select_longlong(differences3, nvals);
+                    diffs5[nrows] = (double) quick_select_longlong(differences5, nvals);
+		}
+
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise3 = 0;
+	       xnoise5 = 0;
+	} else if (nrows == 1) {
+	       xnoise3 = diffs3[0];
+	       xnoise5 = diffs5[0];
+	} else {	    
+	       qsort(diffs3, nrows, sizeof(double), FnCompare_double);
+	       qsort(diffs5, nrows, sizeof(double), FnCompare_double);
+	       xnoise3 =  (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
+	       xnoise5 =  (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
+	}
+
+	if (nrows2 == 0) { 
+	       xnoise2 = 0;
+	} else if (nrows2 == 1) {
+	       xnoise2 = diffs2[0];
+	} else {	    
+	       qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
+	       xnoise2 =  (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise2)  *noise2  = 1.0483579 * xnoise2;
+	if (noise3)  *noise3  = 0.6052697 * xnoise3;
+	if (noise5)  *noise5  = 0.1772048 * xnoise5;
+
+	free(diffs5);
+	free(diffs3);
+	free(diffs2);
+	free(differences5);
+	free(differences3);
+	free(differences2);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise5_float
+       (float *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	float nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	float *minval,    /* minimum non-null value */
+	float *maxval,    /* maximum non-null value */
+	double *noise2,      /* returned 2nd order MAD of all non-null pixels */
+	double *noise3,      /* returned 3rd order MAD of all non-null pixels */
+	double *noise5,      /* returned 5th order MAD of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the median and background noise in the input image using 2nd, 3rd and 5th
+order Median Absolute Differences.
+
+The noise in the background of the image is calculated using the MAD algorithms 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+3rd order:  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
+	float *differences2, *differences3, *differences5;
+	float *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
+	float xminval = FLT_MAX, xmaxval = -FLT_MAX;
+	int do_range = 0;
+	double *diffs2, *diffs3, *diffs5; 
+	double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
+	
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 9 pixels */
+	if (nx < 9) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise2) *noise2 = 0.;
+		if (noise3) *noise3 = 0.;
+		if (noise5) *noise5 = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences2 = calloc(nx, sizeof(float));
+	if (!differences2) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences3 = calloc(nx, sizeof(float));
+	if (!differences3) {
+		free(differences2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences5 = calloc(nx, sizeof(float));
+	if (!differences5) {
+		free(differences2);
+		free(differences3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs2 = calloc(ny, sizeof(double));
+	if (!diffs2) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs3 = calloc(ny, sizeof(double));
+	if (!diffs3) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs5 = calloc(ny, sizeof(double));
+	if (!diffs5) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+		free(diffs3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+			
+		/* find the 5th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v5 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		}
+				
+		/* find the 6th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v6 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v6 < xminval) xminval = v6;
+			if (v6 > xmaxval) xmaxval = v6;
+		}
+				
+		/* find the 7th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v7 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v7 < xminval) xminval = v7;
+			if (v7 > xmaxval) xmaxval = v7;
+		}
+				
+		/* find the 8th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v8 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v8 < xminval) xminval = v8;
+			if (v8 > xmaxval) xmaxval = v8;
+		}
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		nvals2 = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		    v9 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v9 < xminval) xminval = v9;
+			if (v9 > xmaxval) xmaxval = v9;
+		    }
+
+		    /* construct array of absolute differences */
+
+		    if (!(v5 == v6 && v6 == v7) ) {
+		        differences2[nvals2] = (float) fabs(v5 - v7);
+			nvals2++;
+		    }
+
+		    if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
+		        differences3[nvals] = (float) fabs((2 * v5) - v3 - v7);
+		        differences5[nvals] = (float) fabs((6 * v5) - (4 * v3) - (4 * v7) + v1 + v9);
+		        nvals++;  
+		    } else {
+		        /* ignore constant background regions */
+			ngoodpix++;
+		    }
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+		    v5 = v6;
+		    v6 = v7;
+		    v7 = v8;
+		    v8 = v9;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the median diffs */
+		/* Note that there are 8 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 8);
+
+		if (nvals == 0) {
+		    continue;  /* cannot compute medians on this row */
+		} else if (nvals == 1) {
+		    if (nvals2 == 1) {
+		        diffs2[nrows2] = differences2[0];
+			nrows2++;
+		    }
+		        
+		    diffs3[nrows] = differences3[0];
+		    diffs5[nrows] = differences5[0];
+		} else {
+                    /* quick_select returns the median MUCH faster than using qsort */
+		    if (nvals2 > 1) {
+                        diffs2[nrows2] = quick_select_float(differences2, nvals);
+			nrows2++;
+		    }
+
+                    diffs3[nrows] = quick_select_float(differences3, nvals);
+                    diffs5[nrows] = quick_select_float(differences5, nvals);
+		}
+
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise3 = 0;
+	       xnoise5 = 0;
+	} else if (nrows == 1) {
+	       xnoise3 = diffs3[0];
+	       xnoise5 = diffs5[0];
+	} else {	    
+	       qsort(diffs3, nrows, sizeof(double), FnCompare_double);
+	       qsort(diffs5, nrows, sizeof(double), FnCompare_double);
+	       xnoise3 =  (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
+	       xnoise5 =  (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
+	}
+
+	if (nrows2 == 0) { 
+	       xnoise2 = 0;
+	} else if (nrows2 == 1) {
+	       xnoise2 = diffs2[0];
+	} else {	    
+	       qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
+	       xnoise2 =  (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise2)  *noise2  = 1.0483579 * xnoise2;
+	if (noise3)  *noise3  = 0.6052697 * xnoise3;
+	if (noise5)  *noise5  = 0.1772048 * xnoise5;
+
+	free(diffs5);
+	free(diffs3);
+	free(diffs2);
+	free(differences5);
+	free(differences3);
+	free(differences2);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise5_double
+       (double *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	double nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	double *minval,    /* minimum non-null value */
+	double *maxval,    /* maximum non-null value */
+	double *noise2,      /* returned 2nd order MAD of all non-null pixels */
+	double *noise3,      /* returned 3rd order MAD of all non-null pixels */
+	double *noise5,      /* returned 5th order MAD of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the median and background noise in the input image using 2nd, 3rd and 5th
+order Median Absolute Differences.
+
+The noise in the background of the image is calculated using the MAD algorithms 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+3rd order:  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
+	double *differences2, *differences3, *differences5;
+	double *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
+	double xminval = DBL_MAX, xmaxval = -DBL_MAX;
+	int do_range = 0;
+	double *diffs2, *diffs3, *diffs5; 
+	double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
+	
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 9 pixels */
+	if (nx < 9) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise2) *noise2 = 0.;
+		if (noise3) *noise3 = 0.;
+		if (noise5) *noise5 = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences2 = calloc(nx, sizeof(double));
+	if (!differences2) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences3 = calloc(nx, sizeof(double));
+	if (!differences3) {
+		free(differences2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+	differences5 = calloc(nx, sizeof(double));
+	if (!differences5) {
+		free(differences2);
+		free(differences3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs2 = calloc(ny, sizeof(double));
+	if (!diffs2) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs3 = calloc(ny, sizeof(double));
+	if (!diffs3) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs5 = calloc(ny, sizeof(double));
+	if (!diffs5) {
+		free(differences2);
+		free(differences3);
+		free(differences5);
+		free(diffs2);
+		free(diffs3);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+			
+		/* find the 5th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v5 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		}
+				
+		/* find the 6th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v6 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v6 < xminval) xminval = v6;
+			if (v6 > xmaxval) xmaxval = v6;
+		}
+				
+		/* find the 7th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v7 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v7 < xminval) xminval = v7;
+			if (v7 > xmaxval) xmaxval = v7;
+		}
+				
+		/* find the 8th valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v8 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v8 < xminval) xminval = v8;
+			if (v8 > xmaxval) xmaxval = v8;
+		}
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		nvals2 = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		    v9 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v9 < xminval) xminval = v9;
+			if (v9 > xmaxval) xmaxval = v9;
+		    }
+
+		    /* construct array of absolute differences */
+
+		    if (!(v5 == v6 && v6 == v7) ) {
+		        differences2[nvals2] =  fabs(v5 - v7);
+			nvals2++;
+		    }
+
+		    if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
+		        differences3[nvals] =  fabs((2 * v5) - v3 - v7);
+		        differences5[nvals] =  fabs((6 * v5) - (4 * v3) - (4 * v7) + v1 + v9);
+		        nvals++;  
+		    } else {
+		        /* ignore constant background regions */
+			ngoodpix++;
+		    }
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+		    v5 = v6;
+		    v6 = v7;
+		    v7 = v8;
+		    v8 = v9;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the median diffs */
+		/* Note that there are 8 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 8);
+
+		if (nvals == 0) {
+		    continue;  /* cannot compute medians on this row */
+		} else if (nvals == 1) {
+		    if (nvals2 == 1) {
+		        diffs2[nrows2] = differences2[0];
+			nrows2++;
+		    }
+		        
+		    diffs3[nrows] = differences3[0];
+		    diffs5[nrows] = differences5[0];
+		} else {
+                    /* quick_select returns the median MUCH faster than using qsort */
+		    if (nvals2 > 1) {
+                        diffs2[nrows2] = quick_select_double(differences2, nvals);
+			nrows2++;
+		    }
+
+                    diffs3[nrows] = quick_select_double(differences3, nvals);
+                    diffs5[nrows] = quick_select_double(differences5, nvals);
+		}
+
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise3 = 0;
+	       xnoise5 = 0;
+	} else if (nrows == 1) {
+	       xnoise3 = diffs3[0];
+	       xnoise5 = diffs5[0];
+	} else {	    
+	       qsort(diffs3, nrows, sizeof(double), FnCompare_double);
+	       qsort(diffs5, nrows, sizeof(double), FnCompare_double);
+	       xnoise3 =  (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
+	       xnoise5 =  (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
+	}
+
+	if (nrows2 == 0) { 
+	       xnoise2 = 0;
+	} else if (nrows2 == 1) {
+	       xnoise2 = diffs2[0];
+	} else {	    
+	       qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
+	       xnoise2 =  (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise2)  *noise2  = 1.0483579 * xnoise2;
+	if (noise3)  *noise3  = 0.6052697 * xnoise3;
+	if (noise5)  *noise5  = 0.1772048 * xnoise5;
+
+	free(diffs5);
+	free(diffs3);
+	free(diffs2);
+	free(differences5);
+	free(differences3);
+	free(differences2);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise3_short
+       (short *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	short nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	short *minval,    /* minimum non-null value */
+	short *maxval,    /* maximum non-null value */
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the median and background noise in the input image using 3rd order differences.
+
+The noise in the background of the image is calculated using the 3rd order algorithm 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nvals, ngoodpix = 0;
+	short *differences, *rowpix, v1, v2, v3, v4, v5;
+	short xminval = SHRT_MAX, xmaxval = SHRT_MIN, do_range = 0;
+	double *diffs, xnoise = 0, sigma;
+
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 5 pixels */
+	if (nx < 5) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise) *noise = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences = calloc(nx, sizeof(short));
+	if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs = calloc(ny, sizeof(double));
+	if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+		
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		    v5 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		    }
+
+		    /* construct array of 3rd order absolute differences */
+		    if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
+		        differences[nvals] = abs((2 * v3) - v1 - v5);
+		        nvals++;  
+		    } else {
+		        /* ignore constant background regions */
+			ngoodpix++;
+		    }
+
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the 3rd order diffs */
+		/* Note that there are 4 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 4);
+
+		if (nvals == 0) {
+		    continue;  /* cannot compute medians on this row */
+		} else if (nvals == 1) {
+		    diffs[nrows] = differences[0];
+		} else {
+                    /* quick_select returns the median MUCH faster than using qsort */
+                    diffs[nrows] = quick_select_short(differences, nvals);
+		}
+
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise = 0;
+	} else if (nrows == 1) {
+	       xnoise = diffs[0];
+	} else {	    
+
+
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+
+              FnMeanSigma_double(diffs, nrows, 0, 0.0, 0, &xnoise, &sigma, status); 
+
+	      /* do a 4.5 sigma rejection of outliers */
+	      jj = 0;
+	      sigma = 4.5 * sigma;
+	      for (ii = 0; ii < nrows; ii++) {
+		if ( fabs(diffs[ii] - xnoise) <= sigma)	 {
+		   if (jj != ii)
+		       diffs[jj] = diffs[ii];
+		   jj++;
+	        } 
+	      }
+	      if (ii != jj)
+                FnMeanSigma_double(diffs, jj, 0, 0.0, 0, &xnoise, &sigma, status); 
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise)  *noise  = 0.6052697 * xnoise;
+
+	free(diffs);
+	free(differences);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise3_int
+       (int *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	int nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	int *minval,    /* minimum non-null value */
+	int *maxval,    /* maximum non-null value */
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the background noise in the input image using 3rd order differences.
+
+The noise in the background of the image is calculated using the 3rd order algorithm 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nvals, ngoodpix = 0;
+	int *differences, *rowpix, v1, v2, v3, v4, v5;
+	int xminval = INT_MAX, xmaxval = INT_MIN, do_range = 0;
+	double *diffs, xnoise = 0, sigma;
+	
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 5 pixels */
+	if (nx < 5) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise) *noise = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences = calloc(nx, sizeof(int));
+	if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs = calloc(ny, sizeof(double));
+	if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+		
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		    v5 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		    }
+
+		    /* construct array of 3rd order absolute differences */
+		    if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
+		        differences[nvals] = abs((2 * v3) - v1 - v5);
+		        nvals++;  
+		    } else {
+		        /* ignore constant background regions */
+			ngoodpix++;
+		    }
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the 3rd order diffs */
+		/* Note that there are 4 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 4);
+
+		if (nvals == 0) {
+		    continue;  /* cannot compute medians on this row */
+		} else if (nvals == 1) {
+		    diffs[nrows] = differences[0];
+		} else {
+                    /* quick_select returns the median MUCH faster than using qsort */
+                    diffs[nrows] = quick_select_int(differences, nvals);
+		}
+
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise = 0;
+	} else if (nrows == 1) {
+	       xnoise = diffs[0];
+	} else {	    
+
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+
+              FnMeanSigma_double(diffs, nrows, 0, 0.0, 0, &xnoise, &sigma, status); 
+
+	      /* do a 4.5 sigma rejection of outliers */
+	      jj = 0;
+	      sigma = 4.5 * sigma;
+	      for (ii = 0; ii < nrows; ii++) {
+		if ( fabs(diffs[ii] - xnoise) <= sigma)	 {
+		   if (jj != ii)
+		       diffs[jj] = diffs[ii];
+		   jj++;
+	        }
+	      }
+	      if (ii != jj)
+                FnMeanSigma_double(diffs, jj, 0, 0.0, 0, &xnoise, &sigma, status); 
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise)  *noise  = 0.6052697 * xnoise;
+
+	free(diffs);
+	free(differences);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise3_float
+       (float *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	float nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	float *minval,    /* minimum non-null value */
+	float *maxval,    /* maximum non-null value */
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the median and background noise in the input image using 3rd order differences.
+
+The noise in the background of the image is calculated using the 3rd order algorithm 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nvals, ngoodpix = 0;
+	float *differences, *rowpix, v1, v2, v3, v4, v5;
+	float xminval = FLT_MAX, xmaxval = -FLT_MAX;
+	int do_range = 0;
+	double *diffs, xnoise = 0;
+
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 5 pixels to calc noise, so just calc min, max, ngood */
+	if (nx < 5) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise) *noise = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	if (noise) {
+	    differences = calloc(nx, sizeof(float));
+	    if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	    }
+
+	    diffs = calloc(ny, sizeof(double));
+	    if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	    }
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+		
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) {
+			  ii++;
+		        }
+			
+		    if (ii == nx) break;  /* hit end of row */
+		    v5 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		    }
+
+		    /* construct array of 3rd order absolute differences */
+		    if (noise) {
+		        if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
+
+		            differences[nvals] = (float) fabs((2. * v3) - v1 - v5);
+		            nvals++;  
+		       } else {
+		            /* ignore constant background regions */
+			    ngoodpix++;
+		       }
+		    } else {
+		       /* just increment the number of non-null pixels */
+		       ngoodpix++;
+		    }
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the 3rd order diffs */
+		/* Note that there are 4 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 4);
+
+		if (noise) {
+		    if (nvals == 0) {
+		        continue;  /* cannot compute medians on this row */
+		    } else if (nvals == 1) {
+		        diffs[nrows] = differences[0];
+		    } else {
+                        /* quick_select returns the median MUCH faster than using qsort */
+                        diffs[nrows] = quick_select_float(differences, nvals);
+		    }
+		}
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (noise) {
+	    if (nrows == 0) { 
+	       xnoise = 0;
+	    } else if (nrows == 1) {
+	       xnoise = diffs[0];
+	    } else {	    
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+	    }
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise) {
+		*noise  = 0.6052697 * xnoise;
+		free(diffs);
+		free(differences);
+	}
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise3_double
+       (double *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	double nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	long *ngood,        /* number of good, non-null pixels? */
+	double *minval,    /* minimum non-null value */
+	double *maxval,    /* maximum non-null value */
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+
+/*
+Estimate the median and background noise in the input image using 3rd order differences.
+
+The noise in the background of the image is calculated using the 3rd order algorithm 
+developed for deriving the signal to noise ratio in spectra
+(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
+
+  noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
+
+The returned estimates are the median of the values that are computed for each 
+row of the image.
+*/
+{
+	long ii, jj, nrows = 0, nvals, ngoodpix = 0;
+	double *differences, *rowpix, v1, v2, v3, v4, v5;
+	double xminval = DBL_MAX, xmaxval = -DBL_MAX;
+	int do_range = 0;
+	double *diffs, xnoise = 0;
+	
+	if (nx < 5) {
+		/* treat entire array as an image with a single row */
+		nx = nx * ny;
+		ny = 1;
+	}
+
+	/* rows must have at least 5 pixels */
+	if (nx < 5) {
+
+		for (ii = 0; ii < nx; ii++) {
+		    if (nullcheck && array[ii] == nullvalue)
+		        continue;
+		    else {
+			if (array[ii] < xminval) xminval = array[ii];
+			if (array[ii] > xmaxval) xmaxval = array[ii];
+			ngoodpix++;
+		    }
+		}
+		if (minval) *minval = xminval;
+		if (maxval) *maxval = xmaxval;
+		if (ngood) *ngood = ngoodpix;
+		if (noise) *noise = 0.;
+		return(*status);
+	}
+
+	/* do we need to compute the min and max value? */
+	if (minval || maxval) do_range = 1;
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	if (noise) {
+	    differences = calloc(nx, sizeof(double));
+	    if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	    }
+
+	    diffs = calloc(ny, sizeof(double));
+	    if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	    }
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v1 < xminval) xminval = v1;
+			if (v1 > xmaxval) xmaxval = v1;
+		}
+
+		/***** find the 2nd valid pixel in row (which we will skip over) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v2 = rowpix[ii];  /* store the good pixel value */
+		
+		if (do_range) {
+			if (v2 < xminval) xminval = v2;
+			if (v2 > xmaxval) xmaxval = v2;
+		}
+
+		/***** find the 3rd valid pixel in row */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v3 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v3 < xminval) xminval = v3;
+			if (v3 > xmaxval) xmaxval = v3;
+		}
+				
+		/* find the 4nd valid pixel in row (to be skipped) */
+		ii++;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v4 = rowpix[ii];  /* store the good pixel value */
+
+		if (do_range) {
+			if (v4 < xminval) xminval = v4;
+			if (v4 > xmaxval) xmaxval = v4;
+		}
+		
+		/* now populate the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		    v5 = rowpix[ii];  /* store the good pixel value */
+
+		    if (do_range) {
+			if (v5 < xminval) xminval = v5;
+			if (v5 > xmaxval) xmaxval = v5;
+		    }
+
+		    /* construct array of 3rd order absolute differences */
+		    if (noise) {
+		        if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
+
+		            differences[nvals] = fabs((2. * v3) - v1 - v5);
+		            nvals++;  
+		        } else {
+		            /* ignore constant background regions */
+			    ngoodpix++;
+		        }
+		    } else {
+		       /* just increment the number of non-null pixels */
+		       ngoodpix++;
+		    }
+
+		    /* shift over 1 pixel */
+		    v1 = v2;
+		    v2 = v3;
+		    v3 = v4;
+		    v4 = v5;
+	        }  /* end of loop over pixels in the row */
+
+		/* compute the 3rd order diffs */
+		/* Note that there are 4 more pixel values than there are diffs values. */
+		ngoodpix += (nvals + 4);
+
+		if (noise) {
+		    if (nvals == 0) {
+		        continue;  /* cannot compute medians on this row */
+		    } else if (nvals == 1) {
+		        diffs[nrows] = differences[0];
+		    } else {
+                        /* quick_select returns the median MUCH faster than using qsort */
+                        diffs[nrows] = quick_select_double(differences, nvals);
+		    }
+		}
+		nrows++;
+	}  /* end of loop over rows */
+
+	    /* compute median of the values for each row */
+	if (noise) {
+	    if (nrows == 0) { 
+	       xnoise = 0;
+	    } else if (nrows == 1) {
+	       xnoise = diffs[0];
+	    } else {	    
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+	    }
+	}
+
+	if (ngood)  *ngood  = ngoodpix;
+	if (minval) *minval = xminval;
+	if (maxval) *maxval = xmaxval;
+	if (noise) {
+		*noise  = 0.6052697 * xnoise;
+		free(diffs);
+		free(differences);
+	}
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise1_short
+       (short *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	short nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+/*
+Estimate the background noise in the input image using sigma of 1st order differences.
+
+  noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
+
+The returned estimate is the median of the values that are computed for each 
+row of the image.
+*/
+{
+	int iter;
+	long ii, jj, kk, nrows = 0, nvals;
+	short *differences, *rowpix, v1;
+	double  *diffs, xnoise, mean, stdev;
+
+	/* rows must have at least 3 pixels to estimate noise */
+	if (nx < 3) {
+		*noise = 0;
+		return(*status);
+	}
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences = calloc(nx, sizeof(short));
+	if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs = calloc(ny, sizeof(double));
+	if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		/* now continue populating the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		
+		    /* construct array of 1st order differences */
+		    differences[nvals] = v1 - rowpix[ii];
+
+		    nvals++;  
+		    /* shift over 1 pixel */
+		    v1 = rowpix[ii];
+	        }  /* end of loop over pixels in the row */
+
+		if (nvals < 2)
+		   continue;
+		else {
+
+		    FnMeanSigma_short(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+
+		    if (stdev > 0.) {
+		        for (iter = 0;  iter < NITER;  iter++) {
+		            kk = 0;
+		            for (ii = 0;  ii < nvals;  ii++) {
+		                if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
+			            if (kk < ii)
+			                differences[kk] = differences[ii];
+			            kk++;
+		                }
+		            }
+		            if (kk == nvals) break;
+
+		            nvals = kk;
+		            FnMeanSigma_short(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+	              }
+		   }
+
+		   diffs[nrows] = stdev;
+		   nrows++;
+		}
+	}  /* end of loop over rows */
+
+	/* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise = 0;
+	} else if (nrows == 1) {
+	       xnoise = diffs[0];
+	} else {
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+	}
+
+	*noise = .70710678 * xnoise;
+
+	free(diffs);
+	free(differences);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise1_int
+       (int *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	int nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+/*
+Estimate the background noise in the input image using sigma of 1st order differences.
+
+  noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
+
+The returned estimate is the median of the values that are computed for each 
+row of the image.
+*/
+{
+	int iter;
+	long ii, jj, kk, nrows = 0, nvals;
+	int *differences, *rowpix, v1;
+	double  *diffs, xnoise, mean, stdev;
+
+	/* rows must have at least 3 pixels to estimate noise */
+	if (nx < 3) {
+		*noise = 0;
+		return(*status);
+	}
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences = calloc(nx, sizeof(int));
+	if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs = calloc(ny, sizeof(double));
+	if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		/* now continue populating the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		
+		    /* construct array of 1st order differences */
+		    differences[nvals] = v1 - rowpix[ii];
+
+		    nvals++;  
+		    /* shift over 1 pixel */
+		    v1 = rowpix[ii];
+	        }  /* end of loop over pixels in the row */
+
+		if (nvals < 2)
+		   continue;
+		else {
+
+		    FnMeanSigma_int(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+
+		    if (stdev > 0.) {
+		        for (iter = 0;  iter < NITER;  iter++) {
+		            kk = 0;
+		            for (ii = 0;  ii < nvals;  ii++) {
+		                if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
+			            if (kk < ii)
+			                differences[kk] = differences[ii];
+			            kk++;
+		                }
+		            }
+		            if (kk == nvals) break;
+
+		            nvals = kk;
+		            FnMeanSigma_int(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+	              }
+		   }
+
+		   diffs[nrows] = stdev;
+		   nrows++;
+		}
+	}  /* end of loop over rows */
+
+	/* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise = 0;
+	} else if (nrows == 1) {
+	       xnoise = diffs[0];
+	} else {
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+	}
+
+	*noise = .70710678 * xnoise;
+
+	free(diffs);
+	free(differences);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise1_float
+       (float *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	float nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+/*
+Estimate the background noise in the input image using sigma of 1st order differences.
+
+  noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
+
+The returned estimate is the median of the values that are computed for each 
+row of the image.
+*/
+{
+	int iter;
+	long ii, jj, kk, nrows = 0, nvals;
+	float *differences, *rowpix, v1;
+	double  *diffs, xnoise, mean, stdev;
+
+	/* rows must have at least 3 pixels to estimate noise */
+	if (nx < 3) {
+		*noise = 0;
+		return(*status);
+	}
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences = calloc(nx, sizeof(float));
+	if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs = calloc(ny, sizeof(double));
+	if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		/* now continue populating the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		
+		    /* construct array of 1st order differences */
+		    differences[nvals] = v1 - rowpix[ii];
+
+		    nvals++;  
+		    /* shift over 1 pixel */
+		    v1 = rowpix[ii];
+	        }  /* end of loop over pixels in the row */
+
+		if (nvals < 2)
+		   continue;
+		else {
+
+		    FnMeanSigma_float(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+
+		    if (stdev > 0.) {
+		        for (iter = 0;  iter < NITER;  iter++) {
+		            kk = 0;
+		            for (ii = 0;  ii < nvals;  ii++) {
+		                if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
+			            if (kk < ii)
+			                differences[kk] = differences[ii];
+			            kk++;
+		                }
+		            }
+		            if (kk == nvals) break;
+
+		            nvals = kk;
+		            FnMeanSigma_float(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+	              }
+		   }
+
+		   diffs[nrows] = stdev;
+		   nrows++;
+		}
+	}  /* end of loop over rows */
+
+	/* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise = 0;
+	} else if (nrows == 1) {
+	       xnoise = diffs[0];
+	} else {
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+	}
+
+	*noise = .70710678 * xnoise;
+
+	free(diffs);
+	free(differences);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnNoise1_double
+       (double *array,       /*  2 dimensional array of image pixels */
+        long nx,            /* number of pixels in each row of the image */
+        long ny,            /* number of rows in the image */
+	int nullcheck,      /* check for null values, if true */
+	double nullvalue,    /* value of null pixels, if nullcheck is true */
+   /* returned parameters */   
+	double *noise,      /* returned R.M.S. value of all non-null pixels */
+	int *status)        /* error status */
+/*
+Estimate the background noise in the input image using sigma of 1st order differences.
+
+  noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
+
+The returned estimate is the median of the values that are computed for each 
+row of the image.
+*/
+{
+	int iter;
+	long ii, jj, kk, nrows = 0, nvals;
+	double *differences, *rowpix, v1;
+	double  *diffs, xnoise, mean, stdev;
+
+	/* rows must have at least 3 pixels to estimate noise */
+	if (nx < 3) {
+		*noise = 0;
+		return(*status);
+	}
+	
+        /* allocate arrays used to compute the median and noise estimates */
+	differences = calloc(nx, sizeof(double));
+	if (!differences) {
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	diffs = calloc(ny, sizeof(double));
+	if (!diffs) {
+		free(differences);
+        	*status = MEMORY_ALLOCATION;
+		return(*status);
+	}
+
+	/* loop over each row of the image */
+	for (jj=0; jj < ny; jj++) {
+
+                rowpix = array + (jj * nx); /* point to first pixel in the row */
+
+		/***** find the first valid pixel in row */
+		ii = 0;
+		if (nullcheck)
+		    while (ii < nx && rowpix[ii] == nullvalue) ii++;
+
+		if (ii == nx) continue;  /* hit end of row */
+		v1 = rowpix[ii];  /* store the good pixel value */
+
+		/* now continue populating the differences arrays */
+		/* for the remaining pixels in the row */
+		nvals = 0;
+		for (ii++; ii < nx; ii++) {
+
+		    /* find the next valid pixel in row */
+                    if (nullcheck)
+		        while (ii < nx && rowpix[ii] == nullvalue) ii++;
+		     
+		    if (ii == nx) break;  /* hit end of row */
+		
+		    /* construct array of 1st order differences */
+		    differences[nvals] = v1 - rowpix[ii];
+
+		    nvals++;  
+		    /* shift over 1 pixel */
+		    v1 = rowpix[ii];
+	        }  /* end of loop over pixels in the row */
+
+		if (nvals < 2)
+		   continue;
+		else {
+
+		    FnMeanSigma_double(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+
+		    if (stdev > 0.) {
+		        for (iter = 0;  iter < NITER;  iter++) {
+		            kk = 0;
+		            for (ii = 0;  ii < nvals;  ii++) {
+		                if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
+			            if (kk < ii)
+			                differences[kk] = differences[ii];
+			            kk++;
+		                }
+		            }
+		            if (kk == nvals) break;
+
+		            nvals = kk;
+		            FnMeanSigma_double(differences, nvals, 0, 0, 0, &mean, &stdev, status);
+	              }
+		   }
+
+		   diffs[nrows] = stdev;
+		   nrows++;
+		}
+	}  /* end of loop over rows */
+
+	/* compute median of the values for each row */
+	if (nrows == 0) { 
+	       xnoise = 0;
+	} else if (nrows == 1) {
+	       xnoise = diffs[0];
+	} else {
+	       qsort(diffs, nrows, sizeof(double), FnCompare_double);
+	       xnoise =  (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
+	}
+
+	*noise = .70710678 * xnoise;
+
+	free(diffs);
+	free(differences);
+
+	return(*status);
+}
+/*--------------------------------------------------------------------------*/
+static int FnCompare_short(const void *v1, const void *v2)
+{
+   const short *i1 = v1;
+   const short *i2 = v2;
+   
+   if (*i1 < *i2)
+     return(-1);
+   else if (*i1 > *i2)
+     return(1);
+   else
+     return(0);
+}
+/*--------------------------------------------------------------------------*/
+static int FnCompare_int(const void *v1, const void *v2)
+{
+   const int *i1 = v1;
+   const int *i2 = v2;
+   
+   if (*i1 < *i2)
+     return(-1);
+   else if (*i1 > *i2)
+     return(1);
+   else
+     return(0);
+}
+/*--------------------------------------------------------------------------*/
+static int FnCompare_float(const void *v1, const void *v2)
+{
+   const float *i1 = v1;
+   const float *i2 = v2;
+   
+   if (*i1 < *i2)
+     return(-1);
+   else if (*i1 > *i2)
+     return(1);
+   else
+     return(0);
+}
+/*--------------------------------------------------------------------------*/
+static int FnCompare_double(const void *v1, const void *v2)
+{
+   const double *i1 = v1;
+   const double *i2 = v2;
+   
+   if (*i1 < *i2)
+     return(-1);
+   else if (*i1 > *i2)
+     return(1);
+   else
+     return(0);
+}
+/*--------------------------------------------------------------------------*/
+
+/*
+ *  These Quickselect routines are based on the algorithm described in
+ *  "Numerical recipes in C", Second Edition,
+ *  Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
+ *  This code by Nicolas Devillard - 1998. Public domain.
+ */
+
+/*--------------------------------------------------------------------------*/
+
+#define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
+
+static float quick_select_float(float arr[], int n) 
+{
+    int low, high ;
+    int median;
+    int middle, ll, hh;
+
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for (;;) {
+        if (high <= low) /* One element only */
+            return arr[median] ;
+
+        if (high == low + 1) {  /* Two elements only */
+            if (arr[low] > arr[high])
+                ELEM_SWAP(arr[low], arr[high]) ;
+            return arr[median] ;
+        }
+
+    /* Find median of low, middle and high items; swap into position low */
+    middle = (low + high) / 2;
+    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
+    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
+    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
+
+    /* Swap low item (now in position middle) into position (low+1) */
+    ELEM_SWAP(arr[middle], arr[low+1]) ;
+
+    /* Nibble from each end towards middle, swapping items when stuck */
+    ll = low + 1;
+    hh = high;
+    for (;;) {
+        do ll++; while (arr[low] > arr[ll]) ;
+        do hh--; while (arr[hh]  > arr[low]) ;
+
+        if (hh < ll)
+        break;
+
+        ELEM_SWAP(arr[ll], arr[hh]) ;
+    }
+
+    /* Swap middle item (in position low) back into correct position */
+    ELEM_SWAP(arr[low], arr[hh]) ;
+
+    /* Re-set active partition */
+    if (hh <= median)
+        low = ll;
+        if (hh >= median)
+        high = hh - 1;
+    }
+}
+
+#undef ELEM_SWAP
+
+/*--------------------------------------------------------------------------*/
+
+#define ELEM_SWAP(a,b) { register short t=(a);(a)=(b);(b)=t; }
+
+static short quick_select_short(short arr[], int n) 
+{
+    int low, high ;
+    int median;
+    int middle, ll, hh;
+
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for (;;) {
+        if (high <= low) /* One element only */
+            return arr[median] ;
+
+        if (high == low + 1) {  /* Two elements only */
+            if (arr[low] > arr[high])
+                ELEM_SWAP(arr[low], arr[high]) ;
+            return arr[median] ;
+        }
+
+    /* Find median of low, middle and high items; swap into position low */
+    middle = (low + high) / 2;
+    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
+    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
+    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
+
+    /* Swap low item (now in position middle) into position (low+1) */
+    ELEM_SWAP(arr[middle], arr[low+1]) ;
+
+    /* Nibble from each end towards middle, swapping items when stuck */
+    ll = low + 1;
+    hh = high;
+    for (;;) {
+        do ll++; while (arr[low] > arr[ll]) ;
+        do hh--; while (arr[hh]  > arr[low]) ;
+
+        if (hh < ll)
+        break;
+
+        ELEM_SWAP(arr[ll], arr[hh]) ;
+    }
+
+    /* Swap middle item (in position low) back into correct position */
+    ELEM_SWAP(arr[low], arr[hh]) ;
+
+    /* Re-set active partition */
+    if (hh <= median)
+        low = ll;
+        if (hh >= median)
+        high = hh - 1;
+    }
+}
+
+#undef ELEM_SWAP
+
+/*--------------------------------------------------------------------------*/
+
+#define ELEM_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; }
+
+static int quick_select_int(int arr[], int n) 
+{
+    int low, high ;
+    int median;
+    int middle, ll, hh;
+
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for (;;) {
+        if (high <= low) /* One element only */
+            return arr[median] ;
+
+        if (high == low + 1) {  /* Two elements only */
+            if (arr[low] > arr[high])
+                ELEM_SWAP(arr[low], arr[high]) ;
+            return arr[median] ;
+        }
+
+    /* Find median of low, middle and high items; swap into position low */
+    middle = (low + high) / 2;
+    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
+    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
+    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
+
+    /* Swap low item (now in position middle) into position (low+1) */
+    ELEM_SWAP(arr[middle], arr[low+1]) ;
+
+    /* Nibble from each end towards middle, swapping items when stuck */
+    ll = low + 1;
+    hh = high;
+    for (;;) {
+        do ll++; while (arr[low] > arr[ll]) ;
+        do hh--; while (arr[hh]  > arr[low]) ;
+
+        if (hh < ll)
+        break;
+
+        ELEM_SWAP(arr[ll], arr[hh]) ;
+    }
+
+    /* Swap middle item (in position low) back into correct position */
+    ELEM_SWAP(arr[low], arr[hh]) ;
+
+    /* Re-set active partition */
+    if (hh <= median)
+        low = ll;
+        if (hh >= median)
+        high = hh - 1;
+    }
+}
+
+#undef ELEM_SWAP
+
+/*--------------------------------------------------------------------------*/
+
+#define ELEM_SWAP(a,b) { register LONGLONG  t=(a);(a)=(b);(b)=t; }
+
+static LONGLONG quick_select_longlong(LONGLONG arr[], int n) 
+{
+    int low, high ;
+    int median;
+    int middle, ll, hh;
+
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for (;;) {
+        if (high <= low) /* One element only */
+            return arr[median] ;
+
+        if (high == low + 1) {  /* Two elements only */
+            if (arr[low] > arr[high])
+                ELEM_SWAP(arr[low], arr[high]) ;
+            return arr[median] ;
+        }
+
+    /* Find median of low, middle and high items; swap into position low */
+    middle = (low + high) / 2;
+    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
+    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
+    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
+
+    /* Swap low item (now in position middle) into position (low+1) */
+    ELEM_SWAP(arr[middle], arr[low+1]) ;
+
+    /* Nibble from each end towards middle, swapping items when stuck */
+    ll = low + 1;
+    hh = high;
+    for (;;) {
+        do ll++; while (arr[low] > arr[ll]) ;
+        do hh--; while (arr[hh]  > arr[low]) ;
+
+        if (hh < ll)
+        break;
+
+        ELEM_SWAP(arr[ll], arr[hh]) ;
+    }
+
+    /* Swap middle item (in position low) back into correct position */
+    ELEM_SWAP(arr[low], arr[hh]) ;
+
+    /* Re-set active partition */
+    if (hh <= median)
+        low = ll;
+        if (hh >= median)
+        high = hh - 1;
+    }
+}
+
+#undef ELEM_SWAP
+
+/*--------------------------------------------------------------------------*/
+
+#define ELEM_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
+
+static double quick_select_double(double arr[], int n) 
+{
+    int low, high ;
+    int median;
+    int middle, ll, hh;
+
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for (;;) {
+        if (high <= low) /* One element only */
+            return arr[median] ;
+
+        if (high == low + 1) {  /* Two elements only */
+            if (arr[low] > arr[high])
+                ELEM_SWAP(arr[low], arr[high]) ;
+            return arr[median] ;
+        }
+
+    /* Find median of low, middle and high items; swap into position low */
+    middle = (low + high) / 2;
+    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
+    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
+    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
+
+    /* Swap low item (now in position middle) into position (low+1) */
+    ELEM_SWAP(arr[middle], arr[low+1]) ;
+
+    /* Nibble from each end towards middle, swapping items when stuck */
+    ll = low + 1;
+    hh = high;
+    for (;;) {
+        do ll++; while (arr[low] > arr[ll]) ;
+        do hh--; while (arr[hh]  > arr[low]) ;
+
+        if (hh < ll)
+        break;
+
+        ELEM_SWAP(arr[ll], arr[hh]) ;
+    }
+
+    /* Swap middle item (in position low) back into correct position */
+    ELEM_SWAP(arr[low], arr[hh]) ;
+
+    /* Re-set active partition */
+    if (hh <= median)
+        low = ll;
+        if (hh >= median)
+        high = hh - 1;
+    }
+}
+
+#undef ELEM_SWAP
+
+