# include <math.h>
# include <float.h>
# include <stdio.h>
# include <string.h>
# include <xtables.h>	/* For the IRAF_INDEF macros */
# include "pedsky.h"

# define  CGOLD	0.3819660	/* Num. Rec. golden section search ratio */
# define  ZEPS	1.0e-10		/* Num. Rec. zero min. fractional accuracy */
# define  SHFT(a,b,c,d)	(a)=(b);(b)=(c);(c)=(d);
# define  SIGN(a,b)	((b) >= 0.0 ? fabs(a) : -fabs(a))
# define  INDEF		(float)IRAF_INDEFR

/*   FINDSKY_ITER  --  Find best sky value for a NICMOS image.
**		  The input image is not modified. Only the optimal
**		  sky value is returned.
**
**	Based on Brent minimization routine from Numerical Recipes.
**
**
**
**   Revision history:
**   ----------------
**   H. Bushouse	04-May-1999	Implementation
**
*/

int findSky_Iter (TaskInfo *info, SingleNicmosGroup *input,
		  SingleNicmosGroup *Flat, Bool verbose) {

/* Arguments:
**	info	io: task info structure
**	input	 i: input image
**	Flat	 i: flatfield image
**	verbose	 i: verbose output switch
*/

	/* Local variables */
	int	iter;					/* loop index */
	float	mean, median, mode, stdv, min, max;	/* image statistics */
	float	a, b, d, e, etemp;
	float	fu, fv, fw, fx, p, q, r;
	float	tol1, tol2, u, v, w, x, xm;
	float	ax, bx, cx;				/* scaling params */
	SingleNicmosGroup testim;			/* temp image */

	/* Function declarations */
	int n_copyGroup (SingleNicmosGroup *, SingleNicmosGroup *);
	int n_iterstat (SingleNicmosGroup *, int, int, int, int, float, float,
			short, float, int, float *, float *, float *, float *,
			float *, float *);
	int transform (TaskInfo *, SingleNicmosGroup *, SingleNicmosGroup *,
		       Bool);

	/* Copy input image to working image */
	if (n_copyGroup (&testim, input))
	    return (1);

	/* Setup sky scaling params: if the user did not specify them,
	** then get statistics for the image and compute the scaling
	** params from them. */
	if (info->SkySmin == INDEF || info->SkySmax == INDEF) {
	    if (n_iterstat (&testim, info->statlim[0], info->statlim[1],
			    info->statlim[2], info->statlim[3], -FLT_MAX,
			    FLT_MAX, info->BitMask, 5.0, 10, &mean, &median,
			    &mode, &stdv, &min, &max))
		return (1);
	}

	if (info->SkySmin == INDEF) {
	    ax = median - 2.0*stdv;
	    ax = (ax < 0.0 ? ax : 0.0);
	} else
	    ax = info->SkySmin;

	if (info->SkySmax == INDEF)
	    cx = mean + stdv;
	else
	    cx = info->SkySmax;

	bx = (ax + cx) / 2.0;
	
	if (verbose) {
	    sprintf (MsgText, "  Sky scaling limits: %9.7g, %9.7g, %9.7g\n",
		     ax, bx, cx);
	    n_message (MsgText);
	}

	/* Initialize the minimization routine parameters */
	a = (ax < cx ? ax : cx);
	b = (ax > cx ? ax : cx);
	x = w = v = bx;
	e = 0.0;

	/* Subtract initial guess for sky value from working image */
	info->SkyValue = x;

	if (transform (info, &testim, Flat, verbose))
	    return (1);

	if (verbose) {
	    sprintf (MsgText, "  Testing sky value %9.7g\n", info->SkyValue);
	    n_message (MsgText);
	    sprintf (MsgText, "  RMS of corrected image = %9.7g\n", info->rms);
	    n_message (MsgText);
	}

	/* Iterate to find best sky value */
	fw = fv = fx = info->rms;

	for (iter = 1; iter <= info->MaxIter; iter++) {
	
	     xm = 0.5 * (a+b);
	     tol1 = info->Tolerance * fabs(x) + ZEPS;
	     tol2 = 2.0 * tol1;

	     /* Convergence check */
	     if (fabs(x-xm) <= tol2-0.5*(b-a))
		 break;

	     /* Reset params for next iteration */
	     if (fabs(e) > tol1) {
		 r = (x-w) * (fx-fv);
		 q = (x-v) * (fx-fw);
		 p = (x-v)*q - (x-w)*r;
		 q = 2.0 * (q-r);
		 if (q > 0.0)
		     p = -p;
		 q = fabs(q);
		 etemp = e;
		 e = d;


		 if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) ||
		     p >= q*(b-x)) {

		     e = ( x >= xm ? a-x : b-x);
		     d = CGOLD * e;

		 } else {

		     d = p/q;
		     u = x+d;
		     if (u-a < tol2 || b-u < tol2)
			 d = SIGN (tol1, xm-x);
		 }

	     } else {
		 e = (x >= xm ? a-x : b-x);
		 d = CGOLD * e;
	     }

	     u = (fabs(d) >= tol1 ? x+d : x + SIGN(tol1,d));
	     
	     /* Subtract new sky and pedestal values */
	     freeSingleNicmosGroup (&testim);
	     if (n_copyGroup (&testim, input))
		 return (1);
	     info->SkyValue = u;
	     if (transform (info, &testim, Flat, verbose))
		 return (1);
	     if (verbose) {
		 sprintf (MsgText, "  Iteration %d: Testing sky value %9.7g\n",
			  iter, info->SkyValue);
		 n_message (MsgText);
		 sprintf (MsgText, "  RMS of corrected image = %9.7g\n",
			  info->rms);
		 n_message (MsgText);
	     }

	     fu = info->rms;
	     if (fu <= fx) {
		 if (u >= x) 
		     a=x;
		 else
		     b=x;
		 SHFT(v,w,x,u)
		 SHFT(fv,fw,fx,fu)
	     } else {
		 if (u < x)
		     a=u;
		 else
		     b=u;
		 if (fu <= fw || w == x) {
		     v=w;
		     w=u;
		     fv=fw;
		     fw=fu;
		 } else if (fu <= fv || v == x || v == w) {
		     v=u;
		     fv=fu;
		 }
	     }
	} /* End of iteration loop */

	/* Save final sky value */
	info->SkyValue = x;

	if (verbose) {
	    sprintf (MsgText, "  Converged on sky value %9.7g\n",
		     info->SkyValue);
	    n_message (MsgText);
	}

	/* Find pedestal values */
	freeSingleNicmosGroup (&testim);
	if (n_copyGroup (&testim, input))
	    return (1);
	if (transform (info, &testim, Flat, verbose))
	    return (1);

	/* Free working image */
	freeSingleNicmosGroup (&testim);

	/* Successful return */
	return (0);
}

#undef CGOLD
#undef ZEPS
#undef SHFT
#undef SIGN