/* @(#)imaqual.c	17.1.1.1 (ESO-IPG) 01/25/02 17:52:20 */
/* 

   Andreas Kaufer

   imaqual.c

   determine image quality by fitting Gauss-hermite functions 
   
*/

/* system includes */

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

/* general Midas includes */

#include <tbldef.h>
#include <midas_def.h>

/* FEROS specific includes */

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

int main()
{

  /* general */

  int i, j, m, actvals, kun, knul, selected;
  char line[80];

  /* io-related */

  int tid, icol[11];
  int acol, arow;
  int mcol, mrow, msort;
  int nax, imno, Npix[2];
  char unit[64], inst[72];
  float *inima;
  double Start[2], Step[2];

  int inull;
  float fnull;
  double dnull;
  
  char in_table[60], in_image[60];
  float xpos, ypos;

  double *Xgaus, *Ygaus;
  double *A, *coef;
  int ncoef, npixwin;
  float winsize;

  SCSPRO("imaqual");
  TCMNUL(&inull, &fnull, &dnull);
  SCKGETC("IN_A",   1, 60, &actvals, in_image);
  SCKGETC("IN_B",   1, 60, &actvals, in_table);
  SCKRDI("INPUTI",  1, 1,  &actvals, &ncoef, &kun, &knul);
  SCKRDR("INPUTR",  1, 1,  &actvals, &winsize , &kun, &knul);

  if (ncoef >= 4) ncoef = 4;

  /* READ THE COMPLETE INPUT TABLE of observed line centers */

  TCTOPN (in_table, F_IO_MODE, &tid);
  TCIGET(tid, &mcol, &mrow, &msort, &acol, &arow);
  TCCSER(tid, "X", &icol[0]);
  TCCSER(tid, "Y", &icol[1]);

  TCCSER(tid, "GPOS",  &icol[2]);
  TCCSER(tid, "GINT",  &icol[3]);
  TCCSER(tid, "GFWHM", &icol[4]);

  TCCSER(tid, "GH1", &icol[5]);
  TCCSER(tid, "GH2", &icol[6]);
  TCCSER(tid, "GH3", &icol[7]);
  TCCSER(tid, "GH4", &icol[8]);

  
  if (icol[0] == -1 || icol[1] == -1)
  {
    sprintf(line, "Error: missing input column in table %s", in_table);
    SCTPUT(line);
    SCETER(9, "Exiting...");
  }
  if (icol[2] == -1)
    TCCINI (tid,D_R8_FORMAT, 1, "F15.7", "World Coord", ":GPOS", &icol[2]);
  if (icol[3] == -1)
    TCCINI (tid,D_R8_FORMAT, 1, "F15.7", "DN", ":GINT", &icol[3]);
  if (icol[4] == -1)
    TCCINI (tid,D_R8_FORMAT, 1, "F15.7", "World Coord", ":GFWHM", &icol[4]);

  if (icol[5] == -1)
    TCCINI (tid,D_R8_FORMAT, 1, "F15.7", " ", ":GH1", &icol[5]);
  if (icol[6] == -1)
    TCCINI (tid,D_R8_FORMAT, 1, "F15.7", " ", ":GH2", &icol[6]);
  if (icol[7] == -1)
    TCCINI (tid,D_R8_FORMAT, 1, "F15.7", " ", ":GH3", &icol[7]);
  if (icol[8] == -1)
    TCCINI (tid,D_R8_FORMAT, 1, "F15.7", " ", ":GH4", &icol[8]);


  /* map input file                              */

  SCIGET (in_image, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE,
	  2, &nax, Npix, Start, Step, inst, unit, (char **)&inima, &imno);

  npixwin = (int)(winsize/Step[0]);
  Xgaus = dvector (1, 2*npixwin+1);
  Ygaus = dvector (1, 2*npixwin+1);

  /* loop over found lines and fit gauss-hermite function in pix-order frame*/

  A = dvector (1, 3);
  coef = (double *) dvector (1, ncoef); 

  for (m = 1; m <= mrow; m++)
    {
      TCSGET(tid, m, &selected);
      if (selected) {
      TCERDR(tid, m, icol[0], &xpos, &knul); /* x-position of line */
      TCERDR(tid, m, icol[1], &ypos, &knul); /* y-position = echelle order of line */

      for (i=-npixwin; i<=npixwin; i++)
	{
	  j = (int)((ypos-Start[1])/Step[1])*Npix[0] + (int)((xpos-Start[0])/Step[0]) + i;
	  Xgaus[i+npixwin+1] =  (double)(xpos + i*Step[0]);
	  Ygaus[i+npixwin+1] =  (double)(inima[j]);
	}
      
      fit_gh (Xgaus, Ygaus, 2*npixwin+1, A, ncoef, coef); 
      
      /* convert width of Gaussian to FWHM */
      A[3] = A[3]*2.354;

      TCEWRD(tid, m, icol[2], &A[2]);
      TCEWRD(tid, m, icol[3], &A[1]);
      TCEWRD(tid, m, icol[4], &A[3]);

      for (i=1;i<=ncoef;i++) 
	{

	  TCEWRD(tid, m, icol[i+4], &coef[i]);
	}
    }
    }
  sprintf(line, "%d input lines fitted with Gauss-Hermite polynomials", mrow);
  SCTPUT(line);


  /* deallocate memory */
  free_dvector(coef, 1, ncoef);
  free_dvector(A, 1, 3);
  free_dvector (Xgaus, 1, 2*npixwin+1);
  free_dvector (Ygaus, 1, 2*npixwin+1);

  TCTCLO(tid);  
  SCSEPI();
  return 0;
} /* end of main */