/****************************************/
/*	spotcalib	valdes 8 30 82	*/
/*					*/
/*	Generate plate dens to		*/
/*	ints conversion array		*/
/****************************************/
#include	"focas1.h"
#define		SPOTS	16

struct imgdef   phdr;
PIXEL spot[40000];
double  ints[SPOTS];
double  intsF[] = {
    0.000, -0.145, -0.339, -0.540, -0.728, -0.894, -1.082,
    -1.265, -1.476, -1.648, -1.841, -2.028, -2.242, -2.467, -2.636, -2.851
};
double  intsJ[] = {
    0.000, -0.162, -0.352, -0.558, -0.751, -0.964, -1.202,
    -1.438, -1.683, -1.883, -2.111, -2.344, -2.575, -2.792, -2.984, -3.203
};
double  va[4], bm[4][4], am[4][4], cf[4], vb[4];
double  b;
int     lnb[300];

spotcalib () {
    int     i, j, nfit, in[SPOTS], fin, stcmmt;
    char    spfl[12], spotstr[80], strg[80];
    int     nints, npixels, nspots, spdx, spdy, sp1;
    double  den[SPOTS], std[SPOTS], temp;

    int     k, x, y;
    double  f, g, ia, s, ib, d, dqe, q;
    double  ntsft (), gam ();


    printf ("  Enter spots filename: "); fflush (stdout);
    while (strlen (gets (spfl)) == 0);
    switch (sp.band[0]) {
	case 'J': 
	case 'j': 
	    nints = 16;
	    sprintf (comments[(sp.comment++) % ncmmnts],
		    "  Using Kitt Peak Sensitometer Wedge 129 values for band %s.",
		    sp.band);
	    printf ("  %s\n", comments[(sp.comment - 1) % ncmmnts]);
	    for (i = 0; i < nints; i++)
		ints[i] = intsJ[i];
	    break;
	case 'F': 
	case 'f': 
	    nints = 16;
	    sprintf (comments[(sp.comment++) % ncmmnts],
		    "  Using Kitt Peak Sensitometer Wedge 129 values for band %s.",
		    sp.band);
	    printf ("  %s\n", comments[(sp.comment - 1) % ncmmnts]);
	    for (i = 0; i < nints; i++)
		ints[i] = intsF[i];
	    break;
	default: 
	    printf ("  Don't know sensitometer values for band %s\n",
		    sp.band);
	    sprintf (comments[(sp.comment++) % ncmmnts],
		    "  Using user supplied sensitometer values.");
	    printf ("  Enter number of spot values: "); fflush (stdout);
	    while (strlen (gets (strg)) == 0);
	    sscanf (strg, "%d", &nints);
	    if (nints > SPOTS)
		nints = SPOTS;
	    printf ("  Enter %d exposure values.\n", nints);
	    for (i = 0; i < nints; i++) {
		printf ("    Spot %2d : ", i + 1); fflush (stdout);
		while (strlen (gets (strg)) == 0);
		sscanf (strg, "%hf", &ints[i]);
	    }
    }
 /* Sort darkest to lightest */
    for (i = 0; i < nints - 1; i++)
	for (j = i + 1; j < nints; j++) {
	    if (ints[j] > ints[i]) {
		temp = ints[i];
		ints[i] = ints[j];
		ints[j] = temp;
	    }
	}
    if (ints[nints - 1] < 0)
	for (i = 0; i < nints; i++)
	    ints[i] -= ints[nints - 1];

    fin = imopen (spfl, &phdr, 0);
    spdx = spdy = phdr.naxis1;
    sp1 = 0;

    npixels = spdx * spdy;
    if (npixels * phdr.bytepix > MAXBUF)
	npixels = MAXBUF / phdr.bytepix;
    nspots = phdr.naxis2 / spdy;
    if (nspots > SPOTS)
	nspots = SPOTS;
    for (i = 0; i < nspots; i++) {
	den[i] = std[i] = 0.;
	rdimage (fin, 0, i * spdx, npixels, &phdr, spot);
	for (j = 0; j < npixels; j++) {
	    den[i] += spot[j];
	    std[i] += spot[j] * spot[j];
	}
	den[i] /= npixels;
	std[i] = sqrt ((std[i] - npixels * den[i] * den[i]) / (npixels - 1));
	if (std[i] == 0.)
	    std[i] = 1.;
    }
    imclose (fin);
 /* Sort darkest to lightest */
    for (i = 0; i < nspots - 1; i++)
	for (j = i + 1; j < nspots; j++) {
	    if (den[j] > den[i]) {
		temp = den[i];
		den[i] = den[j];
		den[j] = temp;
		temp = std[i];
		std[i] = std[j];
		std[j] = temp;
	    }
	}
    if (nints < nspots)
	nspots = nints;
    printf ("      n  <den>  std   log(i)\n");
    for (i = 0; i < nspots; i++) {
	f = exp (2.303 * ints[sp1 + i]);
	g = (den[i] / std[i]) * (den[i] / std[i]) *.1 / f;
	printf ("   %2d  %6.2f  %6.4f  %6.4f\n",
		sp1 + i, den[i], std[i], ints[sp1 + i]);
    }
    sprintf (comments[(sp.comment++) % ncmmnts],
	    "Intensity calibration from sensitometer spots:");
    sprintf (comments[(sp.comment++) % ncmmnts],
	    "  Sensitometer spot file: %s", spfl);
    sprintf (comments[(sp.comment++) % ncmmnts],
	    "  Scan size of spots: %d x %d", spdx, spdy);
    sprintf (comments[(sp.comment++) % ncmmnts],
	    "  Spots scanned: %d - %d", sp1, sp1 + nspots - 1);
    stcmmt = sp.comment;
    for (;;) {			/* iterate to get eqn parameter 'b' */
	printf ("  Enter points for fit ( <cr> to use previously chosen points): ");
	fflush (stdout);
	gets (strg);
	if (strlen (strg) > 0) {
	    nfit = sscanf (strg, "%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d%d",
		    in, in + 1, in + 2, in + 3, in + 4,
		    in + 5, in + 6, in + 7, in + 8, in + 9, in + 10, in + 11, in + 12, in + 13, in + 14, in + 15);
	    if (nfit <= 0)
		break;
	    strcpy (spotstr, strg);
	}
	printf ("  Enter 'b' ( < cr > to exit): "); fflush (stdout);
	gets (strg);
	j = sscanf (strg, "%hf", &b);
	if (j <= 0)
	    break;
	sp.comment = stcmmt;
	sprintf (comments[(sp.comment++) % ncmmnts],
		"  n < den > log(i)    fit    gamma   DQE");
	printf ("  %s\n", comments[(sp.comment - 1) % ncmmnts]);
	for (i = 0; i < 4; i++) {
	    va[i] = 0.0;
	    for (j = 0; j < 4; j++)
		am[i][j] = 0.0;
	}
	for (i = 0; i < nspots; i++) {
	    for (j = 0; j < nfit; j++)
		if (i + sp1 == in[j])
		    break;
	    if (j == nfit)
		continue;
	    vb[0] = 1.0;
	    d = den[i] / sp.satr;
	    ia = ints[sp1 + i];
	    vb[1] = d;
	    g = exp (b * d);
	    vb[2] = log (g - 1.0);
	    vb[3] = g;
	    for (j = 0; j < 4; j++) {
		va[j] += ia * vb[j];
		for (k = 0; k < 4; k++)
		    am[j][k] += vb[j] * vb[k];
	    }
	}
	gauss3 (4, 0.0001, am, bm);
	vmul (4, 4, bm, va, cf);
	j = 0;
	for (i = 1; i < 135; i++) {
	    d = i / 50.0;
	    s = 150.0 * ntsft (d) + 40;
	    if ((s < 40.0) || (s > 600.0))
		continue;
	    lnb[j++] = s;
	    lnb[j++] = 120.0 * d + 40.0;
	}
	s = 0;
	for (i = 0; i < nspots; i++) {
	    for (j = 0; j < nfit; j++)
		if (i + sp1 == in[j])
		    break;
	    if (j == nfit)
		continue;
	    d = den[i] / sp.satr;
	    ia = ints[sp1 + i];
	    x = 150.0 * ia + 40.0;
	    y = 120.0 * d + 40.0;
	    lnb[0] = x - 3;
	    lnb[2] = x + 3;
	    lnb[1] = y - 3;
	    lnb[3] = y + 3;
	    ib = ntsft (d);
	    q = exp (ib * 2.303);
	    g = gam (d);
	    dqe = g * g / (std[i] * std[i]) / q;
	    sprintf (comments[(sp.comment++) % ncmmnts],
		    "%3d: %7.2f %7.4f %7.4f %7.4f %7.4f", sp1 + i, d * sp.satr, ia, ib, g, sp.satr * dqe);
	    printf ("  %s\n", comments[(sp.comment - 1) % ncmmnts]);
	    s += (ib - ia) * (ib - ia);
	}
	printf ("  Std of fit = %8.4f\n", sqrt (s / nspots));
	printf ("  Parameters: %2.0f %e %e %e %e\n", b, cf[0], cf[1], cf[2], cf[3]);
    }
    sprintf (comments[(sp.comment++) % ncmmnts],
	    "  Std of fit: %8.4f", sqrt (s / nspots));
    sprintf (comments[(sp.comment++) % ncmmnts],
	    "  Fit parameters: %2.0f %e %e %e %e", b, cf[0], cf[1], cf[2], cf[3]);
    setflm ();
}

double  ntsft (ds)
double  ds;
{
    double  a, f, g;
    ds = ds >.001 ? ds : 0.001;
    g = exp (ds * b);
    f = log (g - 1.0);
    a = cf[0] + ds * cf[1] + f * cf[2] + g * cf[3];
    a = a < 0.? 0.: a;
    a = a > 4.? 4.: a;
    return (a);
}

double  gam (ds)
double  ds;
{
    double  g;
    g = exp (ds * b);
    return (1./ (cf[1] + b * g * (cf[2] / (g - 1) + cf[3])));
}

setflm () {
    int     i;
    float   flm ();
    char    strg[20];

    sp.pflags &= ~LNR;
    sp.pflags &= ~INTP;
    sp.pflags &= ~DLUT;
    if (sp.satr == 0) {
	printf ("  Saturation value: "); fflush (stdout);
	gets (strg);
	sscanf (strg, "%d", &sp.satr);
    }
    if (sp.satr > MXDEN) {
	sp.pflags |= FUNC;
	sp.nden = 5;
	flmcv[0].i = b;
	flmcv[1].i = cf[0];
	flmcv[2].i = cf[1];
	flmcv[3].i = cf[2];
	flmcv[4].i = cf[3];
	return;
    }
    sp.pflags |= DLUT;
    sp.nden = sp.satr;

    for (i = 0; i < sp.nden; i++) {
	flmcv[i].d = i;
	flmcv[i].i = flm (i);
    }
    for (i = 0; i < sp.nden - 1; i++)
	flmcv[i].slope = (flmcv[i + 1].i - flmcv[i].i) / (flmcv[i + 1].d - flmcv[i].d);
}

float   flm (d)
int     d;
{
    float   i, c, f, g;

    if (b == 0.)
	return ((float) i);
    if (d == 0)
	return (0.);
    i = (float) d / sp.satr;
    i = i >.001 ? i : 0.001;
    g = exp (b * i);
    f = log (g - 1.0);
    c = cf[0] + i * cf[1] + f * cf[2] + g * cf[3];
    c = c < 0.? 0.: c;
    c = c > 4.? 4.: c;
    return (exp (2.303 * c));
}