/*	initsky			*/
/*				*/
/*	Establish initial sky	*/
/********************************/
#include	"focas1.h"

float  *dsky, th0, th1, wts, wts2;
int     xfltsz, yfltsz;
PIXEL * lnptrs[NYFLT];
float  *smthln;

intsky (minsky, maxsky, skyblk)
float   minsky, maxsky;
int	skyblk;
{
    int     i, j, l;
    float   nn, ns, ssum1;
    float   ssum2, ssum3, ssum4;
    float   a0, a1, b, c, sigma, min, max;
    PIXEL * cpr;

 /* Find line average */
    ssum1 = 0;
    ssum2 = ssum3 = ssum4 = 0.;
    for (i = 0; i < sp.fy; i++) {
	cpr = lnptrs[i];
	for (j = 0; j < sp.kpts; j++) {
	    c = *cpr++;
	    if ((c < minsky) || (c > maxsky))
		continue;
	    ssum1++;
	    ssum2 += c;
	    ssum3 += c * c;
	}
    }
    if (ssum1 > 0)
	c = ssum2 * wts / ssum1;
    else
	c = 0.;
    sigma = ssum3 * ssum1 - ssum2 * ssum2;
    if (sigma <= 0.)
	sigma = 0.;
    else
	sigma = sqrt (wts2) * sqrt (sigma) / ssum1;
    if (sigma < sp.skyhw)
	sigma = sp.skyhw;
    sigma *= sqrt (wts2);

 /* printf("intsky: Initial sky average = %g\n", c / wts); */

    for (i = 0; i < (1 + sp.kpts / skyblk); i++)
	dsky[i] = c;

 /* iterate to get initial sky values */
    max = (maxsky * wts - dsky[0]) / 100.;
    min = (minsky * wts - dsky[0]) / 100.;

    pixel (lnptrs);

    i = 0;
    for (ns = nn = 0.; i < 101; i++) {
	ssum1 = 0;
	ssum2 = ssum3 = ssum4 = 0.;
	for (j = xfltsz; j < (sp.kpts - xfltsz); j++) {
	    l = j / skyblk;
	    if (((j-xfltsz) % skyblk) == 0) {
		if (sp.thrln0 < 0.)
		    a0 = th0;
		else
		    a0 = dsky[l] + th0 + max * (100 - i);
		if (sp.thrln1 < 0.)
		    a1 = th1;
		else
		    a1 = dsky[l] - th1 + min * (100 - i);
	    }

	/* Detection filter value */
	    b = smthln[j];

	    if ((b < a0) && (b > a1)) {
		c = b - dsky[l];
		nn++;
		ns += b;
		ssum1++;
		ssum2 += c;
		ssum3 += c * c;
		ssum4 += b;
	    }
	    if (((j-xfltsz) % skyblk) == 0) {
		l = l > 0 ? l : 1;
		if (nn > sp.nupdte)
		    dsky[l] = dsky[l] + sp.sa * (ns / nn - dsky[l]) +
			2 * sp.sb * ((dsky[l + 1] + dsky[l - 1]) / 2 - dsky[l]);
		else
		    dsky[l] = dsky[l] + 2 * sp.sb * ((dsky[l + 1] + dsky[l - 1]) / 2 - dsky[l]);
		nn = ns = 0.;
	    }
	}
	dsky[l + 1] = dsky[l];
	dsky[0] = dsky[1];
	sigma = ssum3 * ssum1 - ssum2 * ssum2;
	if (sigma <= 0.)
	    sigma = 0.;
	else
	    sigma = sqrt (sigma) / ssum1;
	if (sigma < sp.skyhw * sqrt (wts2))
	    sigma = sp.skyhw * sqrt (wts2);
	if (sp.thrln0 >= 0.)
	    th0 = sigma * sp.thrln0;
	if (sp.thrln1 >= 0.)
	    th1 = sigma * sp.thrln1;

    /* printf("%3d: %4.0f  %g\n", i, ssum1, th0); */

    }
    if (ssum1 > 0)
	ssum4 /= ssum1 * wts;
    if (sp.skyhw == 0.)
	sp.skyhw = sigma / sqrt (wts2);
    else
	sigma = sp.skyhw * sqrt (wts2);
    if (sp.thrln0 < 0.)
	a0 = -sp.thrln0;
    else {
	th0 = sigma * sp.thrln0;
	a0 = th0 / sqrt (wts2);
    }
    if (sp.thrln1 < 0.)
	a1 = -sp.thrln1;
    else {
	th1 = sigma * sp.thrln1;
	a1 = th1 / sqrt (wts2);
    }
    sprintf (comments[sp.comment % ncmmnts],
	    "sky: pts = %4.0f avg = %g sigma = %g\nthreshold above = %g threshold below = %g\n",
	    ssum1, ssum4, sigma / sqrt (wts2), a0, a1);

 /* printf("%s\n", comments[sp.comment%ncmmnts]); */

    sp.comment++;
}