#include    "focas1.h"


/* SKY -- Evaluate sky around objects. */

sky (catalog, skytype, sbr, buffer, width)

char	*catalog;
int	*skytype;
float	*sbr;
int	*buffer;
int	*width;

{
    int     i, n, nsbr, entnum;
    short   rct[4];
    long    szrec, szhdr, szobjrec(), szcathdr();
    float   ssbr, r1, prob, delta;
    int     hstgrm[512];
    PIXEL   step, start;

    struct objrec   ob;
    static struct areas    area;
    static struct image    pic;

    cfd = catopen (catalog, 2);
    cmmnt ("sky");
    wtcathdr (cfd, 0);

    if (*width > 0)
	sp.skwdth = *width;
    if (*buffer <= 0)
	*buffer = sp.buf - sp.skwdth;
    *buffer -= sp.buf - sp.skwdth;
    sp.skybw = 0.02;
    pfd = fndfld (sp.ptfl, 0);
    sp.xs = pthdr.crpix1;
    sp.ys = pthdr.crpix2;
    afd = fndar (sp.arfl, 0);
    szhdr = szcathdr();
    szrec = szobjrec();

    step = sp.skybw * sp.skyhw;
    if (step < 1)
	step = 1;

    delta = 0.;
    n = 0;
    for (;;) {
	if (rdcatob (cfd, 0L, &ob))
	    break;
	if (ffilter (&ob) != 1)
	    continue;

	switch (*skytype) {
	    case 2: 
		ob.sbr = *sbr;
		ob.ssbr = 0.;
		ob.nsbr = 1000;
		ob.r1 = 0.;
		ob.prob = 0.;
		break;
	    case 3: 
		ob.sbr += *sbr;
		break;
	    default: 
		if (rdarea (afd, ob.entnum, ob.subent, ob.arpos, &area))
		    continue;

		rct[0] = max (area.rct[0] - *buffer, 0);
		rct[1] = max (area.rct[1] - *buffer, 0);
		rct[2] = min (area.rct[2] + *buffer, pthdr.naxis1-1);
		rct[3] = min (area.rct[3] + *buffer, pthdr.naxis2-1);

		if (rdimg (pfd, rct, &pic) < 0)
		    continue;

		if (ob.subent == 0) {
		    evlsky (&ob, &area, &pic, hstgrm, 512, step, &start,
			*skytype);
		    delta += ob.r1 - ob.sbr;
		    n++;
		} else {
		    ob.sbr = *sbr;
		    ob.ssbr = ssbr;
		    ob.nsbr = nsbr;
		    ob.r1 = r1;
		    ob.prob = prob;
		}

		*sbr = ob.sbr;
		ssbr = ob.ssbr;
		nsbr = ob.nsbr;
		r1 = ob.r1;
		prob = ob.prob;
		entnum = ob.entnum;
	}

	fseek (cfd, -szrec, 1);
	wtcat (&ob, cfd);
    }

    if (*skytype == 4 && n > 0) {
	delta /= n;
	printf ("%s: <Mode - Mean> = %g\n", catalog, delta);

	fseek (cfd, szhdr, 0);
	for (;;) {
	    if (rdcatob (cfd, 0L, &ob))
		break;
	    if (ffilter (&ob) != 1)
		continue;
	    ob.sbr += delta;
	    fseek (cfd, -szrec, 1);
	    wtcat (&ob, cfd);
	}
    }

    fclose (afd);
    fclose (cfd);
}