#include <math.h>
#include <curses.h>
#include <stdio.h>
#include "nrutil.h"
#include "string.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, struct image *mask)
{
    extern int (*pfunc) (const char *, ...);
    extern char device[];
    extern struct runflags flags;
    float estsky (struct image *data, struct image *mask, float *mean, 
	int flag, float skyrms);
    float mean = 0.;

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

    sigma = (struct image *) calloc (1, (size_t)(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, 0.01, 0.01); 

    /******************************************************************\ 
    *  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);

    if (flags.noestsky == 1) {
	skysig = 0.;
	pfunc ("\n");
        pfunc ("-- WARNING: Intrinsic sky is not estimated to create the internal sigma image.\n"); 
        pfunc ("            Zero pix values are set to (eff_rdnoise)/(eff_gain) in sigma image.\n");
    } else { 

	/***********************************************************\
	*  If the image is sky subtracted, the dark current, night  *
	*  sky and readnoise terms can be summed up into one term:  *
        *  the RMS of the sky background, skysig.                   *
	\***********************************************************/

	if (flags.noestsky == 2) 
	    mean = flags.skyped;

            skysig = estsky (d, mask, &mean, flags.noestsky, flags.skyrms);

	if (skysig > 0)
	    effrdnoise = 0.;
    }; 

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

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

            varp = fabs((sigma->z[j][i] - mean) * effgain);

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

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

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

        };
    };

    if (flags.outsig) {
        sprintf (sigma->name, "sigma.fits");
        mkimg (sigma->z, sigma->naxes, sigma->name, "sigma.fits");
    };

    return (sigma);
}

/***************************************************************************\
*  Estimate the mean and RMS of the sky pedestal.                           *
\***************************************************************************/

float estsky (struct image *data, struct image *mask, float *mean, int flag, 
    float skyrms)
{

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

    int i, j, npix=0, ngood=0;
    float stdsky = 0., stdsky2=0., skysig = 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++) {
	    if (mask->z[j][i] == 0) {
	        ngood++;
	        fz[ngood] = data->z[j][i];
	    };
	};
    };

    sort (ngood, fz);  

    if (flag == 0)
        median = fz [(int)(ngood/2)];
    else
	median = *mean;

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

    if (skyrms < 0) {
        for (i=(int)(ngood * 0.2); i <= (int)(ngood * 0.8); ++i) 
	    stdsky += (fz[i] - median) * (fz[i] - median);

        stdsky = sqrt (stdsky/(ngood*3/5.-1));
    } else
	stdsky = skyrms;

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

    free_vector (fz, 1, npix); 


    if (flag == 0) {

        /*********************************************\
        *  Print message regarding sky estimation:    *
        \*********************************************/

        pfunc ("\n");
        pfunc ("-- Estimating the sky mean and RMS to generate an internal sigma image.  This\n");
        pfunc ("   estimate won't work right if the product ADUxGAIN has a unit other than \n");
        pfunc ("   [electrons] (e.g. MJy/sr or e-/sec). To turn off this feature, start by:   \n");
        pfunc ("                 galfit -noskyest <filename>  \n");
        pfunc ("   For other options, do:  \n");
        pfunc ("                 galfit -help  \n\n");

        /******************************************************************\
        *  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) <= 5 * stdsky) {
		    if (mask->z[j][i] == 0) {
	                *mean += data->z[j][i];
		        npix++;
		    };
	        };
	    };
        };

        *mean = *mean / (npix - 1);

    } else if (flag == 2) {

        pfunc ("\n");
        pfunc ("-- Using the RMS and user provided sky mean to generate an internal sigma image.  \n");
        pfunc ("   This estimate won't work right if the product ADUxGAIN has a unit other \n");
        pfunc ("   [electrons] (e.g. microjanskies). To turn off this feature, start by:   \n");
        pfunc ("                 galfit -noskyest <filename>  \n\n");

    };
    
    if (skyrms < 0) {
        for (j=1; j <= data->naxes[2]; ++j) {
	    for (i=1; i <= data->naxes[1]; ++i) {
		if (mask->z[j][i] == 0 &&
	            fabs(data->z[j][i] - *mean) <= 5 * stdsky)
	                stdsky2 += (data->z[j][i] - *mean) * (data->z[j][i] - 
				    *mean);
	    };
        };

        skysig = sqrt (stdsky2/(npix-1));
    } else
	skysig = skyrms;

    pfunc ("   Est'd sky mean = %.2f, RMS = %.2f ADUs.  (Are these fairly reasonable?)\n", *mean, skysig);

    return (skysig);
     
}