#include <math.h>
#include <curses.h>
#include <stdio.h>
#include "nrutil.h"
#include "structs.h"

#define NCUTOFF 10.

void smooth (float **in, float **out, long naxes[], float sigma, float radius);
void mkimg (float **array, long naxes[], char *outname, char *tname);
void sort(unsigned long n, float arr[]);

struct image *sigma_img (struct image *d)
{
    extern int (*pfunc) (const char *, ...);
    extern char device[];
    float estsky (struct image *data, float effgain, float effrdnoise);

    int i, j, npix;
    float var, varp, varn, effgain, effrdnoise, skyped=0.;
    struct image *sigma;

    sigma = (struct image *) calloc (1, sizeof (struct image));
    if (sigma == (struct image *) 0) {
        pfunc ("Error in allocating memory for mini-image....");
        if (strncmp(device, "curses", 6)== 0)
            refresh();
        exit (1);
    };

    npix = d->naxes[1] * d->naxes[2]; 
    sigma->naxes[1] = d->naxes[1];
    sigma->naxes[2] = d->naxes[2];

    sigma->z = matrix (1, sigma->naxes[2], 1, sigma->naxes[1]);
    effgain = FMAX(d->gain * d->ncombine, 1.);

    smooth (d->z, sigma->z, sigma->naxes, 2., 3.); 

/*    smooth (d->z, sigma->z, sigma->naxes, 10, 20.);  */

    /******************************************************************\ 
    *  Check to see if we need to add a sky pedestal, i.e., if the     *
    *  value is less than 1, then return a best guess sky pedestal.    *
    *  Otherwise, return 0.                                            *
    \******************************************************************/

    effrdnoise = d->rdnoise * sqrt (d->ncombine) / sqrt (d->nsamp+d->nread);
    skyped = estsky (d, effgain, effrdnoise);

    /*********************************\
    *  Now create the variance array  *
    \*********************************/

    for (j=1; j <= sigma->naxes[2]; ++j) {
        for (i=1; i <= sigma->naxes[1]; ++i) {

/*	    if (d->z[j][i] == 0.)
		sigma->z[j][i] = 1.e12;
	    else {
*/              varp = (sigma->z[j][i] + skyped) * effgain;

	        if (varp <= -3 * effrdnoise) 
	            varp = 1.e30;

                varn = pow(effrdnoise, 2.);
	        var = varp + varn;
	        if (var <= 1.)
	            var = 1.;

 	        sigma->z[j][i] = sqrt(var) / effgain; 
/*	    }; */
        };
    };

    return (sigma);

}


/*
    sprintf (sigma->name, "sigma.fits");
    mkimg (sigma->z, sigma->naxes, sigma->name, "sigma.fits");
*/

/***************************************************************************\
*  Estimate the sky pedestal based on the RMS fluctuation of the sky        *
\***************************************************************************/

float estsky (struct image *data, float effgain, float effrdnoise)
{

    extern int (*pfunc) (const char *, ...);

    int i, j, npix=0, n=0;
    float stdsky = 0., stdsky2 = 0., mean = 0., skyped=0., median;
    float *fz;

    npix = data->naxes[2] * data->naxes[1];
    fz = vector (1, npix);

    for (i=1; i <= data->naxes[1]; i++) {
        for (j=1; j <= data->naxes[2]; j++) {
	    n++;
	    fz[n] = data->z[j][i];
	};
    };

    sort (npix, fz);  
    median = fz [(int)(npix/2)];

    /* Exclude the top 20 and bottom 15% when calculating stddev  */

    for (i=(int)(npix * 0.15); i <= (int)(npix * 0.8); ++i) 
	   stdsky += (fz[i] - median) * (fz[i] - median);

    stdsky = sqrt (stdsky/(npix*0.65-1));

/*    pfunc ("median = %.2f, stdsky = %.2f, npix=%d \n", median, stdsky, npix); */

    free_vector (fz, 1, npix); 

    /*********************************************************************\
    *  Either print warning message that we're generating internal sigma  * 
    *  image or quit immediately (if the sky is present)                  *
    \*********************************************************************/

    if (0.5*pow(stdsky*effgain, 2) < median*effgain + pow(effrdnoise,2.) &&
	2 * pow(stdsky*effgain, 2) > median*effgain + pow(effrdnoise,2.) )

        return (0.);    
    else {
        pfunc ("\n");
        pfunc ("   WARNING: The input data has been sky subtracted.  Estimating original\n");
        pfunc ("   sky to generate an internal sigma image.  This estimate won't work right\n");
        pfunc ("   if the ADU has a unit other than counts (e.g. counts/sec, microjanskies).\n\n");
/*      pfunc ("   Median sky value = %f\n\n", median); */
    };

    /******************************************************************\
    *  Iterate on the mean and stddev one more time, rejecting +/- 3   *
    *  data pixels as objects.                                         *
    \******************************************************************/

    npix = 0;
    for (j=1; j <= data->naxes[2]; ++j) {
	for (i=1; i <= data->naxes[1]; ++i) {
	    if (fabs(data->z[j][i] - median) <= 3 * stdsky) {
	        mean += data->z[j][i];
		npix++;
	    };
	};
    };

    mean = mean / npix;
    
    for (j=1; j <= data->naxes[2]; ++j) {
	for (i=1; i <= data->naxes[1]; ++i) {
	    if (fabs(data->z[j][i] - mean) <= 3 * stdsky)
	        stdsky2 += (data->z[j][i] - mean) * (data->z[j][i] - mean);
	};
    };

    stdsky2 = sqrt (stdsky2/(npix-1));
    skyped = sqrt (pow (stdsky2*effgain, 2) - pow(effrdnoise,2))/effgain;

    pfunc ("   Est'd sky mean = %.2f, stddev = %.1f ADUs.  (Are these reasonable?)\n", mean, stdsky2);
    pfunc ("   Est'd original sky = %.1f ADUs.\n\n", skyped); 

    if (skyped <= 0.) {
	pfunc ("   WARNING: Bad sky estimate.  The GAIN and/or READNOISE may be incorrect.");
	pfunc ("            Assuming READNOISE is 0.");

        skyped = stdsky2/effgain;
  	pfunc ("Assume a sky value of %.2f ADUs.\n", skyped);
    };

    return (skyped - mean);
     
}