/* @(#)fit_line.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:58 */
/* 

   Otmar Stahl

   fit_line.c 

   determine center of profile

*/

/* FEROS specific includes */

#include <echelle.h>
#include <proto_nrutil.h>
#include <proto_fitnol.h>
#include <proto_mutil.h>

int fit_line
#ifdef __STDC__
(
 float rval[], int index, float obj[],
 double start[], double step[], int method, 
 int win, int width, float rnull
 )
#else
     (
      rval, index, obj, start, step, method, 
      win, width, rnull
      )
     float rval[], obj[], rnull;
     double start[], step[];
     int method, win, width;
#endif

{
  int k;
  float factor, shift, a, b, aleft, aright;
  double x_guess;
  double *A, *xgaus, *ygaus;

  A = dvector(1, 3);
  xgaus = dvector(1, 2 * win + 1);
  ygaus = dvector(1, 2 * win + 1);  


  switch (method)
    {  
    case GRAVITY:    /* center of gravity */  

      aleft = rval[index - 1];
      aright = rval[index + 1];
      factor = 1;
      if (aleft >= aright)
	{
	  aleft = rval[index + 1];
	  aright = rval[index - 1];
	  factor = -1;
	}
      a = rval[index] - aleft;
      b = aright - aleft;
      shift = (a + b == 0.0) ? 0.0 : step[0] * b / (a + b);
      obj[0] = start[0] + step[0] * index + factor * shift;
      obj[2] = rval[index];
      break;

    case GAUSSIAN: /* gaussian fit */

      A[1] = rval[index];                /* A */
      A[2] = start[0] + step[0] * index; /* x0 */
      x_guess = A[2];
      A[3] = step[0];                    /* sigma */

      for(k = -win; k <= win; k++)
	{
	  xgaus[k + win + 1] = start[0] + step[0] * (index + k);
	  ygaus[k + win + 1] = rval[index + k];
	}

      fit_gauss(xgaus, ygaus, width, A, 3); /* writes width, A */
      obj[0] = A[2];
      obj[1] = A[3];
      obj[2] = A[1];

      /* check if result is out of limits */

      if (myabs (A[2] - x_guess)  > win) 
      {
	obj[0] = rnull;
	obj[1] = rnull;
	obj[2] = rnull;
      }
      break;
    }
  free_dvector(A, 1, 3);
  free_dvector(xgaus, 1, 2 * win + 1);
  free_dvector(ygaus, 1, 2 * win + 1);
  return 0;
}