#include	"focas1.h"


/* TEMPLATE -- Make a PSF template from catalog objects */

mktemplate (catalog, psffile, nx, ny, rflag, verbose)
char	*catalog, *psffile;
int	*nx, *ny, *rflag, *verbose;
{
    struct objrec   ob;
    static struct image    p;
    FILE	*fout;
    int		i, j, k, l, n, xc, yc;
    float	ir1, lum, val;
    short	rct[4];
    PIXEL	*cp;

    if ((*nx > NXMAX) || (*ny > NYMAX))
	focaserr (1, "template: Desired template is too large", "");

    cfd = catopen (catalog, 0);
    afd = fndar (sp.arfl, 0);
    pfd = fndfld (sp.ptfl, 0);
    sp.nx = *nx;
    sp.ny = *ny;
    xc = (sp.nx - 1) / 2;
    yc = (sp.nx - 1) / 2;
 
    p.pic = 0;
    for (i = 0; i < sp.ny; i++)
	for (j = 0; j < sp.nx; j++)
	    template[i][j] = 0.;
    for (n = 0;; n++) {
	if (rdcatob (cfd, 0L, &ob))
	    break;
	rct[0] = lpcoords (ob.xc, 1, 2) +.5 - xc;
	rct[1] = lpcoords (ob.yc, 2, 2) +.5 - yc;
	rct[2] = rct[0] + sp.nx - 1;
	rct[3] = rct[1] + sp.ny - 1;
	if (rdimg (pfd, rct, &p) < 0) {
	    n--;
	    continue;
	}
	cp = p.pic;
	for (i = 0, lum = 0.; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++)
		lum += film (*cp++) - ob.sbr;
	cp = p.pic;
	for (i = 0; i < sp.ny; i++)
	    for (j = 0; j < sp.nx; j++)
		if (*rflag) {
		    k = 2 * yc - i;
		    l = 2 * xc - j;
		    if (i <= sp.nx && j <= sp.ny) {
		        val = (film (*cp++) - ob.sbr) / (8 * lum);
		        template[i][j] += val;
		        template[k][j] += val;
		        template[i][l] += val;
		        template[k][l] += val;
		        template[j][i] += val;
		        template[l][i] += val;
		        template[j][k] += val;
		        template[l][k] += val;
		    } else {
		        val = (film (*cp++) - ob.sbr) / (4 * lum);
		        template[i][j] += val;
		        template[k][j] += val;
		        template[i][l] += val;
		        template[k][l] += val;
		    }
		} else
		    template[i][j] += (film (*cp++) - ob.sbr) / lum;
    }
    fclose (afd);
    fclose (cfd);
    imclose (pfd);

    for (i = 0; i < sp.ny; i++)
	for (j = 0; j < sp.nx; j++)
	    template[i][j] /= n;

    if (strcmp (psffile, "STDOUT") == 0)
	fout = stdout;
    else {
	fout = fopen (psffile, "w");
	if (fout == NULL) {
	    c_delete (psffile);
	    fout = fopen (psffile, "w");
	    if (fout == NULL) {
		focaserr (1, "template: Can't create ", psffile);
		return;
	    }
	}
    }

    fprintf (fout, "%2d %2d\n", sp.nx, sp.ny);
    ir1 = lum = 0;
    for (i = 0; i < sp.ny; i++)
	for (j = 0; j < sp.nx; j++) {
	    fprintf (fout, "%g\n", template[i][j]);
	    ir1 += sqrt ((float) (i-yc)*(i-yc)+(j-xc)*(j-xc)) * template[i][j];
	    lum += template[i][j];
	}
    fclose (fout);

    if (verbose) {
	fprintf (stderr, "ir1 = %g\n", ir1 / lum);
	lum = template[xc][yc] / 100.;
	for (i = 0; i < sp.ny; i++) {
	    for (j = 0; j < sp.nx; j++)
		fprintf (stderr, "%4d", (int) (template[i][j] / lum + 0.5));
	    fprintf (stderr, "\n");
	}
    }
}