/********************************/
/*      skyhst			*/
/*      Find sky histogram      */
/********************************/
#include    "focas1.h"

int	lseed = 1;

skyhst (ar, pic, hstgrm, nhstgrm, step, start)
struct arblk   *ar;
struct image   *pic;
int     hstgrm[];
int     nhstgrm;
PIXEL step, start;

{
    float   sum, bin, mean, mean1, mean2;
    int     npix;
    int     i, j, l, xl, xr, xi, yi;
    PIXEL * cp;

 /* Histogram */

    npix = 0;
    for (i = 0; i < nhstgrm; i++)
	hstgrm[i] = 0.;

    for (i = pic -> dy - 1, l = 0; i >= 0; i--) {
	yi = ar -> ychst - pic -> y0 - i;
	cp = pic -> pic + (i + 1) * pic -> dx - 1;
	if ((yi >= 0) && (yi < ar -> nlines)) {
	    for (xi = pic -> dx - 1; l < ar -> nx[yi]; l++) {
		xl = ar -> xla[l] - pic -> x0;
		xr = ar -> xra[l] - pic -> x0;
		for (; xi > xr; xi--) {
		    if (xi >= pic -> dx - sp.skwdth) {
			j = ((*cp--) - start) / step;
			if ((j >= 0) && (j < nhstgrm)) {
			    npix++;
			    hstgrm[j]++;
			}
		    }
		    else
			cp--;
		}
		cp -= xi - xl + 1;
		xi = xl - 1;
	    }
	    for (; xi >= 0; xi--) {
		if (xi < sp.skwdth) {
		    j = ((*cp--) - start) / step;
		    if ((j >= 0) && (j < nhstgrm)) {
			npix++;
			hstgrm[j]++;
		    }
		}
		else
		    cp--;
	    }
	}
	else
	    if ((i < sp.skwdth) || (i >= pic -> dy - sp.skwdth)) {
		for (xi = pic -> dx - 1; xi >= 0; xi--) {
		    j = ((*cp--) - start) / step;
		    if ((j >= 0) && (j < nhstgrm)) {
			npix++;
			hstgrm[j]++;
		    }
		}
	    }
	    else {
		for (xi = pic -> dx - 1; xi >= 0; xi--) {
		    if ((xi < sp.skwdth) || (xi >= pic -> dx - sp.skwdth)) {
			j = ((*cp--) - start) / step;
			if ((j >= 0) && (j < nhstgrm)) {
			    npix++;
			    hstgrm[j]++;
			}
		    }
		    else
			cp--;
		}
	    }
    }

    return (npix);
}

/*Return a random value based on the histogram in the interval min to max*/
ransky (min, max, peak, hstgrm)
int     min, max, peak;
int     hstgrm[];
{
    float   focasrand();
    int     i, j;

    for (;;) {
	i = focasrand(&lseed) * (max - min) + min;
	j = focasrand(&lseed) * peak;
	if ((i < min) || (i > max))
	    continue;
	if (j <= hstgrm[i])
	    return (i);
    }
}