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

   Otmar Stahl
   
   fitfunc.c 

   Gauss-Hermite fitting routine
   determines Gauss-Hermite moments of line profiles

*/

/* system includes */

#include <math.h>
#include <stdio.h>

/* general Midas includes */


#include <midas_def.h>

/* FEROS specific includes */

#include <proto_nrutil.h>
#include <proto_fitnol.h>
#include <misc.h>
#include <fitfunc.h>
#include <dlfit.h>

/*           Fit Gauss-Hermite function                        */
/*           uses   dlfit, fit_gauss                           */
void 
#ifdef __STDC__
fit_gh (double *Xgaus, double *Ygaus, int npix, double *A, 
	int ncoef, double *coef)
#else
     fit_gh (Xgaus, Ygaus, npix, A, ncoef, coef)
     double *Xgaus, *Ygaus, *A, *coef; 
     int   npix, ncoef;
#endif
{
  
  double **covar, *sig, chisq;
  
  int i, imax, *ia;
  double max;

  imax = npix / 2;
  max = -1e99;

  for (i=1; i<=npix; i++)
    {
      if (Ygaus[i]>max)
	{
	  max = Ygaus[i];
	  imax = i;
	}
    }
  A[1] = Ygaus[imax];
  A[2] = Xgaus[imax];
  A[3] = fabs(Xgaus[2]-Xgaus[1])*2.;
      
  fit_gauss (Xgaus, Ygaus, npix, A, 3);
  
  for (i=1; i<=npix; i++)
    {
      Xgaus[i] = (Xgaus[i] - A[2]) / A[3] ;
      Ygaus[i] = Ygaus[i] / A[1];
    }

  covar = (double **) dmatrix (1,npix, 1, npix); 
  sig = (double *) dvector (1, npix); 
  ia = (int *) ivector (1, ncoef); 

  for (i = 1; i <= ncoef; i++) 
    ia[i] = 1; 
  for (i= 1; i <= npix; i++) 
    sig[i] = 1.0; 
  
  dlfit(Xgaus, Ygaus, sig, npix, coef, ia, ncoef, covar, 
	&chisq, dhermite); 

}

void 
#ifdef __STDC__
dhermite (double x, double p[], int np)
#else
     dhermite (x, p, np)
     double x, p[];
     int np;
#endif
{

  /* Gauss-Hermite fitting routine */

#ifdef __STDC__
  int facul[10] = {1,1,2,6,24,120,720,5040,40320,362880};
  int powerof2[10] = {1,2,4,8,16,32,64,128,256,512};
#else
  int facul[10], powerof2[10];
#endif
  int i;

  double exp1, h[10];

#ifndef __STDC__
  facul[1] = 1;
  facul[2] = 1;
  facul[3] = 2;
  facul[4] = 6;
  facul[5] = 24;
  facul[6] = 120;
  facul[7] = 720;
  facul[8] = 5040;
  facul[9] = 40320;
  facul[10] = 362880;
 
  powerof2[1] = 1;
  powerof2[2] = 2;
  powerof2[3] = 4;
  powerof2[4] = 8;
  powerof2[5] = 16;
  powerof2[6] = 32;
  powerof2[7] = 64;
  powerof2[8] = 128;
  powerof2[9] = 256;
  powerof2[10] = 512;
#endif

  /* compute Gauss-Hermite functions from recursion relation 
     Up to 8 fit parameters are allowed */

  h[0] = 1.0;
  h[1] = 2 * x;
  
  for (i=2; i<=np+1; i++)
    h[i] = 2 * x  * h[i-1] - 2 * (i-1) * h[i-2];

  /* h1 and h2 are zero for best fitting gaussian */

  exp1 = exp( -(x*x) / 2.0 );
  
  p[1] = h[0] * exp1;
  
  for (i=2; i<=np; i++)
    {
      p[i] = exp1 * h[i+1] /  sqrt( (double) ( facul[i+1] * powerof2[i+1]));
    }
}

double 
#ifdef __STDC__
fdhermite (double x, double p[], int np)
#else
     fdhermite (x, p, np)
     double x, p[];
     int np;
#endif
{

  /* Gauss-Hermite evaluation routine */

#ifdef __STDC__
  int facul[10] = {1,1,2,6,24,120,720,5040,40320,362880};
  int powerof2[10] = {1,2,4,8,16,32,64,128,256,512};
#else
  int facul[10], powerof2[10];
#endif

  double exp1, result, h[10];
  int i;

#ifndef __STDC__
  facul[1] = 1;
  facul[2] = 1;
  facul[3] = 2;
  facul[4] = 6;
  facul[5] = 24;
  facul[6] = 120;
  facul[7] = 720;
  facul[8] = 5040;
  facul[9] = 40320;
  facul[10] = 362880;
 
  powerof2[1] = 1;
  powerof2[2] = 2;
  powerof2[3] = 4;
  powerof2[4] = 8;
  powerof2[5] = 16;
  powerof2[6] = 32;
  powerof2[7] = 64;
  powerof2[8] = 128;
  powerof2[9] = 256;
  powerof2[10] = 512;
#endif

  h[0] = 1.0;
  h[1] = 2 * x;
  
  for (i=2; i<=np+1; i++)
    h[i] = 2 * x  * h[i-1] - 2 * (i-1) * h[i-2];

  exp1 = exp( -(x*x) / 2.0 );

  result = p[1] * h[0];

  for (i=2; i<=np; i++)
    {
      result += p[i] * h[i+1] /  sqrt( (double) ( facul[i+1] * powerof2[i+1]));
    }

  result *= exp1;
  return (float) result;
}