/* @(#)mgauss.c	17.1.1.1 (ES0-DMD) 01/25/02 17:26:59 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* @(#)mgauss.c	17.1.1.1  (ESO)  01/25/02  17:26:59 */

#include <math.h>
#include <midas_def.h>
#include <alice_util.h>
#include <alice_global.h>

#define ALICETABLE "TMPalice.tbl"
#define EPSILON 1e-15

void            fgauss();
float           fit_cont();
void            fpoly();

int             tid;
int             col_ident, col_fgauss, col_fline, col_xcen, col_xfwhm,
                col_cont, col_eqwt;
int             col_error, col_ystart, col_yend, currline;
float           line_error, xint1, xint2, yint1, yint2;
double          fitContError;

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  mgauss() : execute a multi-gaussian fitting and plot the solution. 
---------------------------------------------------------------------*/
mgauss()
{
  float          *xgauss, *ygauss;
  double         *dyda, *a, *sigma, **covar, **alpha, chisq, alamda;
  double          ochisq, oalamda;
  double          contin;
  int             ndata, ma, *lista, mfit;
  int             i, j;
  char            out[20];
  double          getFitValue();

  ma = gaussNumOfSol * 3;
  mfit = gaussNumOfSol * 3;
  ndata = specNpix[0];
  covar = dmatrix(1, ma, 1, ma);
  alpha = dmatrix(1, ma, 1, ma);
  lista = ivector(1, ma);
  a = dvector(1, ma);
  dyda = dvector(1, ma);
  xgauss = vector(1, ndata);
  ygauss = vector(1, ndata);
  sigma = dvector(1, ndata);
  alamda = -1;
  for (i = 1; i <= ma; i++)
    a[i] = (double) gaussFitValues[i - 1];
  for (i = 1; i <= ma; i++)
    a[i] = getFitValue(a, i, ma);
  for (i = 1, mfit = 1; i <= ma; i++)
    if (!gaussFixOpt[i - 1])
      lista[mfit++] = i;
  mfit--;
  ochisq = 1e20;
  i = 0;
  while (specX[i++] < specXcen - specDx & specX[i] < specXmax);
  i--;
  j = 1;
  while (specX[i + j] < specXcen + specDx & specX[i + j] < specXmax)
  {
    xgauss[j] = specX[i + j];
    contin = fit_cont(xgauss[j]);
    ygauss[j] = specY[i + j] - contin;
    j++;
  }
  ndata = j - 1;
  for (i = 1; i <= ndata; i++)
    sigma[i] = 1.0;
  mrqmin(xgauss, ygauss, sigma, ndata, a, ma, lista, mfit, covar, alpha,
         &chisq, fgauss, &alamda);
  for (j = 1; j <= ma; j++)
    a[j] = getFitValue(a, j, ma);
  i = 1;
  while (i++ < gaussMinIterations)
  {
    if (i % 4 == 0)
      ChangeCurs();
    mrqmin(xgauss, ygauss, sigma, ndata, a, ma, lista, mfit, covar, alpha,
           &chisq, fgauss, &alamda);
    for (j = 1; j <= ma; j++)
      a[j] = getFitValue(a, j, ma);
  }
  while (i++ < gaussMaxIterations & ochisq != chisq)
  {
    ochisq = chisq;
    mrqmin(xgauss, ygauss, sigma, ndata, a, ma, lista, mfit, covar, alpha,
           &chisq, fgauss, &alamda);
    for (j = 1; j <= ma; j++)
      a[j] = getFitValue(a, j, ma);
    if (i % 4 == 0)
      ChangeCurs();
  }
  alamda = 0.0;
  mrqmin(xgauss, ygauss, sigma, ndata, a, ma, lista, mfit, covar, alpha,
         &chisq, fgauss, &alamda);
  put_iterations(i - 1);
  DefaultCurs();
  for (i = 1; i <= ma; i++)
  {
    gaussFitValues[i - 1] = (float) a[i];
    gaussErrors[i - 1] = covar[i][i] == 0.0 ? 0.0 : sqrt(covar[i][i]);
  }
  fitRms = chisq;
  out_errors();
  ndata = specNpix[0];
  free_dvector(sigma, 1, ndata);
  free_vector(ygauss, 1, ndata);
  free_vector(xgauss, 1, ndata);
  free_dvector(dyda, 1, ma);
  free_dvector(a, 1, ma);
  free_ivector(lista, 1, ma);
  free_dmatrix(alpha, 1, ma, 1, ma);
  free_dmatrix(covar, 1, ma, 1, ma);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 sgauss(ycen, xcen, xstep) : execute a single gaussian fitting and
                             plot the solution.
             xcen , ycen   : x-cent and amplitude.
                    xstep  : sigma.
---------------------------------------------------------------------*/
sgauss(ycen, xcen, xstep)
  float           xcen, ycen, xstep;
{
  float          *xgauss, *ygauss;
  double         *dyda, *a, *sigma, **covar, **alpha, chisq, alamda;
  double          ochisq, oalamda;
  double          contin;
  double          fwhm, area, eqwt, conta2;
  int             ndata, ma, *lista, mfit;
  int             i, j;
  char            out[80];

  ma = 3;
  mfit = 3;
  ndata = specNpix[0];
  covar = dmatrix(1, ma, 1, ma);
  alpha = dmatrix(1, ma, 1, ma);
  lista = ivector(1, ma);
  a = dvector(1, ma);
  dyda = dvector(1, ma);
  xgauss = vector(1, ndata);
  ygauss = vector(1, ndata);
  sigma = dvector(1, ndata);
  alamda = -1;
  a[1] = (double) ycen - fit_cont(xcen);
  a[2] = (double) xcen;
  a[3] = (double) xstep;
  for (i = 1; i <= ma; i++)
    lista[i] = i;
  ochisq = 1e20;
  i = 0;
  while (specX[i++] < specXcen - specDx & specX[i] < specXmax);
  i--;
  j = 1;
  while (specX[i + j] < specXcen + specDx & specX[i + j] < specXmax)
  {
    xgauss[j] = specX[i + j];
    contin = fit_cont(xgauss[j]);
    ygauss[j] = specY[i + j] - contin;
    j++;
  }
  ndata = j - 1;
  for (i = 1; i <= ndata; i++)
    sigma[i] = 1.0;
  mrqmin(xgauss, ygauss, sigma, ndata, a, ma, lista, mfit, covar, alpha,
         &chisq, fgauss, &alamda);
  /* while( ochisq != chisq) */
  while ((ochisq - chisq) / chisq > EPSILON)
  {
    ochisq = chisq;
    mrqmin(xgauss, ygauss, sigma, ndata, a, ma, lista, mfit, covar, alpha,
           &chisq, fgauss, &alamda);
  }
  alamda = 0.0;
  mrqmin(xgauss, ygauss, sigma, ndata, a, ma, lista, mfit, covar, alpha,
         &chisq, fgauss, &alamda);
  draw_sgauss(a[1], a[2], a[3], 4);
  fwhm = 2.35482 * a[3];
  conta2 = fit_cont(a[2]);
  if (conta2 == 0.0)
    specFluxReal = 0.0;
  eqwt = specFluxReal / conta2;
  area = 2.50663 * a[1] * a[3];
  line_error = fitContError * 
               sqrt(fabs(specStep * (eqwt + 2 * (xint2 - xint1))));
  sprintf(out, "%8.5f   %8.5f   %6.0f   %9.5g   %8.5f   %8.5f", a[2],
    fwhm, conta2, specFluxReal, eqwt, line_error);
  SCTPUT(out);
  put_table_values(specFrameIdent, area, specFluxReal, a[2], fwhm,
       conta2, eqwt, line_error, specLineNum, specLineNum + specLineStep - 1);
  currline++;
  ndata = specNpix[0];
  free_dvector(sigma, 1, ndata);
  free_vector(ygauss, 1, ndata);
  free_vector(xgauss, 1, ndata);
  free_dvector(dyda, 1, ma);
  free_dvector(a, 1, ma);
  free_ivector(lista, 1, ma);
  free_dmatrix(alpha, 1, ma, 1, ma);
  free_dmatrix(covar, 1, ma, 1, ma);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  draw_gauss() : draw the multi-gaussian fitting solution.
---------------------------------------------------------------------*/
draw_gauss()
{
  double          wl, wr, contin;
  double         *a, *dyda;
  float           xline[1000], yline[1000], yres;
  float           fit_cont();
  int             i;

  a = dvector(1, 3 * gaussNumOfSol);
  dyda = dvector(1, 3 * gaussNumOfSol);
  for (i = 1; i <= 3 * gaussNumOfSol; i++)
    a[i] = (double) gaussFitValues[i - 1];
  wl = specXcen - specDx;
  wr = specXcen + specDx;
  for (i = 0; i < 1000; i++)
  {
    xline[i] = (float) i *(wr - wl) / 1000.0 + wl;

    fgauss(xline[i], a, &yres, dyda, 3 * gaussNumOfSol);
    contin = fit_cont(xline[i]);
    yline[i] = yres + contin;
  }
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_MOPN(PLOTAPP);
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  AG_SSET("lstyle=0;lwidt=2;color=2");
  AG_GPLL(xline, yline, 1000);
  AG_VUPD();
  AG_MCLS();
  AG_CLS();
  free_dvector(dyda, 1, 3 * gaussNumOfSol);
  free_dvector(a, 1, 3 * gaussNumOfSol);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
draw_error() : draw the residual.
---------------------------------------------------------------------*/
draw_error()
{
  double          wl, wr, contin;
  double         *a, *dyda;
  float           xline[MAXVALUES], yline[MAXVALUES], yres;
  float           fit_cont();
  int             i;

  a = dvector(1, 3 * gaussNumOfSol);
  dyda = dvector(1, 3 * gaussNumOfSol);
  for (i = 1; i <= 3 * gaussNumOfSol; i++)
    a[i] = (double) gaussFitValues[i - 1];
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_MOPN(PLOTAPP);
  AG_SSET("lwidth=0;lstyle=0");
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  xline[0] = specX[0];
  xline[1] = specX[specNpix[0] - 1];
  yline[0] = specYcen + specDy / 1.5;
  yline[1] = specYcen + specDy / 1.5;
  AG_GPLL(xline, yline, 2);
  for (i = 0; i < specNpix[0]; i++)
  {
    xline[i] = specX[i];
    fgauss(xline[i], a, &yres, dyda, 3 * gaussNumOfSol);
    contin = fit_cont(xline[i]);
    yline[i] = yres + contin - specY[i] + specYcen + specDy / 1.5;
  }
  AG_SSET("color = 2");
  AG_GPLL(xline, yline, specNpix[0]);
  AG_VUPD();
  AG_MCLS();
  AG_CLS();
  free_dvector(dyda, 1, 3 * gaussNumOfSol);
  free_dvector(a, 1, 3 * gaussNumOfSol);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 draw_sgauss(a1, a2, a3, color) : draw a single gaussian curve.
       a1  : x-cent.
       a2  : Amplitude.
       a3  : sigma.
---------------------------------------------------------------------*/
draw_sgauss(a1, a2, a3, color)
  float           a1, a2, a3;
  int             color;
{
  double          wl, wr, contin;
  double         *a, *dyda;
  float           xline[100], yline[100], yres;
  float           fit_cont();
  int             i;
  char            sset[80];

  a = dvector(1, 3);
  dyda = dvector(1, 3);
  sprintf(sset, "lstyle=2;lwidth=0;color=%d", color);
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_MOPN(PLOTAPP);
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  AG_SSET(sset);
  AG_MCLS();
  AG_SSET("lstyle=0");
  AG_MOPN(PLOTAPP);
  a[1] = (double) a1;
  a[2] = (double) a2;
  a[3] = (double) a3;
  wl = a[2] - 3 * a[3];
  wr = a[2] + 3 * a[3];
  for (i = 0; i < 100; i++)
  {
    xline[i] = (float) i *(wr - wl) / 100.0 + wl;

    fgauss(xline[i], a, &yres, dyda, 3);
    contin = fit_cont(xline[i]);
    yline[i] = yres + contin;
  }
  AG_GPLL(xline, yline, 100);
  AG_VUPD();
  AG_MCLS();
  /* AG_SSET("lstyle=0"); AG_GPLL(xline,yline,100); */
  AG_CLS();
  free_dvector(dyda, 1, 3 * gaussNumOfSol);
  free_dvector(a, 1, 3 * gaussNumOfSol);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
draw_number_comp(comp, sol, color) : put the number of the gaussian 
                 component over the solution.
---------------------------------------------------------------------*/
draw_number_comp(comp, sol, color)
  int             comp, sol, color;
{
  char            out[80];

  sprintf(out, "color=%d;chdi=1.2,1.2;changl=0.0", color);
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_MOPN(PLOTAPP);
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  AG_SSET(out);
  sprintf(out, "%d", sol + 1);
  AG_GTXT(gaussFitValues[comp * 3 + 1],
    gaussFitValues[comp * 3] + fit_cont(gaussFitValues[comp * 3 + 1]), out, 12);
  AG_VUPD();
  AG_MCLS();
  AG_CLS();
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 plot_fit(xv, yv, npt, n, color)
 xv,yv  : x and y values.
 npt    : number of data.
 n      : fit degree.
 color  : plot color.
---------------------------------------------------------------------*/
void 
plot_fit(xv, yv, npt, n, color)
  float           xv[], yv[];
  int             npt, n, color;
{
  int            *lista;
  int             i, k;
  double          chisq, *a, *sig, **covar, deltax, *powx;
  float           xl[2], yl[2], fit_cont();
  float           SQR();
  char            sset[40];

  lista = ivector(1, n);
  a = dvector(1, n);
  sig = dvector(1, npt);
  covar = dmatrix(1, n, 1, n);
  powx = dvector(1, n);
  for (i = 1; i <= npt; i++)
    sig[i] = 1.0;
  for (i = 1; i <= n; i++)
    lista[i] = i;
  lfit(xv, yv, sig, npt, a, n, lista, n, covar, &chisq, fpoly);
  for (i = 0; i <= fitDegree; i++)
    fitPolyValues[i] = a[i + 1];
  for (; i < 20; i++)
    fitPolyValues[i] = 0.0;
  deltax = 2 * specDx / 100.0;
  fpoly(specXcen - specDx, powx, n);
  xl[1] = specXcen - specDx;
  yl[1] = a[1] * powx[1];
  for (i = 2; i <= n; i++)
    yl[1] += a[i] * powx[i];
  sprintf(sset, "lstyle=1;lwidth=0;color=%d", color);
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_MOPN(PLOTAPP);
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  AG_SSET(sset);
  AG_MCLS();
  AG_SSET("lstyle=0");
  AG_MOPN(PLOTAPP);
  for (k = 1; k <= 99; k++)
  {
    xl[0] = xl[1];
    yl[0] = yl[1];
    xl[1] += deltax;
    fpoly(xl[1], powx, n);
    yl[1] = a[1] * powx[1];
    for (i = 2; i <= n; i++)
      yl[1] += a[i] * powx[i];
    AG_GPLL(xl, yl, 2);
  }
  AG_MCLS();
  AG_CLS();
  fitContError = 0.0;
  for (i = 1; i <= npt; i++)
    fitContError += SQR(fit_cont(xv[i]) - yv[i]);
  fitContError = (float) sqrt(fitContError / npt);
  free_dvector(powx, 1, n);
  free_dmatrix(covar, 1, n, 1, n);
  free_dvector(sig, 1, npt);
  free_ivector(lista, 1, n);
  free_dvector(a, 1, n);
  save_cont("TMPcont.bdf");
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 calc_fit(xv, yv, npt, n) : calculate the polynomial fit of degree n
                            and set the global fit array.
---------------------------------------------------------------------*/
calc_fit(xv, yv, npt, n)
  float           xv[], yv[];
  int             npt, n;
{
  int            *lista, i;
  double          chisq, *a, *sig, **covar, deltax, *powx;

  lista = ivector(1, n);
  a = dvector(1, n);
  sig = dvector(1, npt);
  covar = dmatrix(1, n, 1, n);
  powx = dvector(1, n);
  for (i = 1; i <= npt; i++)
    sig[i] = 1.0;
  for (i = 1; i <= n; i++)
    lista[i] = i;
  lfit(xv, yv, sig, npt, a, n, lista, n, covar, &chisq, fpoly);
  for (i = 0; i <= fitDegree; i++)
    fitPolyValues[i] = a[i + 1];
  for (; i < 20; i++)
    fitPolyValues[i] = 0.0;
  free_dvector(powx, 1, n);
  free_dmatrix(covar, 1, n, 1, n);
  free_dvector(sig, 1, npt);
  free_ivector(lista, 1, n);
  free_dvector(a, 1, n);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  plot_spline(n, color) : plot a line using splaine interpolation.
---------------------------------------------------------------------*/
void 
plot_spline(n, color)
  int             n, color;
{
  int             i, k;
  float           xl[2], yl[2], fit_cont(), deltax;
  char            sset[40];

  deltax = 2 * specDx / 100.0;
  sprintf(sset, "lstyle=1;lwidth=0;color=%d", color);
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_MOPN(PLOTAPP);
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  AG_SSET(sset);
  AG_MCLS();
  AG_SSET("lstyle=0");
  AG_MOPN(PLOTAPP);
  xl[1] = specXcen - specDx;
  yl[1] = fit_cont(xl[1]);
  for (k = 1; k <= 99; k++)
  {
    xl[0] = xl[1];
    yl[0] = yl[1];
    xl[1] += deltax;
    yl[1] = fit_cont(xl[1]);
    AG_GPLL(xl, yl, 2);
  }
  AG_MCLS();
  AG_CLS();
  save_cont("TMPcont.bdf");
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 fit_cont(xcont) : return the interpolated y value using polynomial or
                   spline fitting.
---------------------------------------------------------------------*/
float 
fit_cont(xcont)
  float           xcont;
{
  int            *lista;
  int             i, k;
  float           yret, dy;
  double         *powx;

  if (fitMode == POLYFITMODE)
  {
    powx = dvector(1, fitDegree + 1);
    fpoly(xcont, powx, fitDegree + 1);
    yret = fitPolyValues[0] * powx[1];
    for (i = 1; i <= fitDegree; i++)
      yret += fitPolyValues[i] * powx[i + 1];
    free_dvector(powx, 1, fitDegree + 1);
    return (yret);
  } else if (fitMode == SPLINEFITMODE)
  {
    ratint(specXaux, specYaux, gaussNumOfFitData, xcont, &yret, &dy);
    return (yret);
  } else
  {
    printf("\rContinuum error : 0.0 returned\n");
    return (0);
  }
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  add_fit(color) : add new x values to the fitting X array.
          color  : arrow color.
---------------------------------------------------------------------*/
add_fit(color)
  int             color;
{
  int             key, pix;
  float           xfit1, xfit2, yfit1, yfit2, aux;
  char            sset[20];

  key = 1;
  fitAddFit = 1;
  sprintf(sset, "COLOR=%d", color);
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  AG_SSET(sset);
  AG_SSET("CURSOR = 2");
  xfit1 = specXcen;
  yfit1 = aux = specYcen;
  if (fitMode == POLYFITMODE)
  {
    AG_SSET("SCALE = 1.5");
    while (key == 1)
    {
      xfit1 = xfit2;
      yfit1 = aux;
      AG_VLOC(&xfit1, &yfit1, &key, &pix);
      if (key == 1)
      {
        for (i = 0; specX[i] < xfit1; i++);
        AG_GTXT(xfit1, specY[i], "\\downarro", 2);
        xfit2 = xfit1;
        yfit2 = yfit1;
        AG_VLOC(&xfit2, &yfit2, &key, &pix);
        if (key == 1)
        {
          for (i = 0; specX[i] < xfit2; i++);
          AG_GTXT(xfit2, specY[i], "\\downarro", 2);
          if (xfit1 > xfit2)
          {
            aux = xfit1;
            xfit1 = xfit2;
            xfit2 = aux;
          }
          aux = yfit1;
          for (i = 0; specX[i] < xfit1; i++);
          yfit1 = specY[i];
          fitXminPair[fitPairNum] = xfit1;
          fitXmaxPair[fitPairNum] = xfit2;
          fitPairNum++;
          for (; specX[i] < xfit2; i++, gaussNumOfFitData++)
          {
            specXaux[gaussNumOfFitData] = specX[i];
            specYaux[gaussNumOfFitData] = specY[i];
          }
          yfit2 = specY[i];
        }
      }
    }
    AG_SSET("SCALE = 1.0");
  } else
    while (key == 1)
    {
      xfit1 = xfit2;
      yfit1 = aux;
      AG_VLOC(&xfit1, &yfit1, &key, &pix);
      if (key == 1)
      {
        gaussNumOfFitData++;
        specXaux[gaussNumOfFitData] = xfit1;
        specYaux[gaussNumOfFitData] = yfit1;
        AG_GPLM(specXaux + gaussNumOfFitData, specYaux + gaussNumOfFitData,
                  1, 5);
      }
    }
  AG_VUPD();
  AG_CLS();
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 auto_fit(spec) : calculate the continuum for the mentionated spectrum.
          spec  : 0 = current spectrum; 1-15 overplot spectrum.
---------------------------------------------------------------------*/
auto_fit(spec)
  int             spec;
{
  int             i, j = -1;
  float          *yvalues, *xvalues;
  extern          ifit();

  xvalues = spec == 0 ? specX : OverX[spec - 1];
  yvalues = spec == 0 ? specY : OverY[spec - 1];

  for (i = 0; i < fitPairNum; i++)
  {
    for (j = 0; xvalues[j] < fitXminPair[i]; j++);
    for (; xvalues[j] < fitXmaxPair[i]; j++, gaussNumOfFitData++)
    {
      specXaux[gaussNumOfFitData] = xvalues[j];
      specYaux[gaussNumOfFitData] = yvalues[j];
    }
  }

}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  integrate() : called for Integrate push button.
---------------------------------------------------------------------*/
integrate()
{
  int             key, pix;
  float           aux, xcen, ycen, xstep;
  char            sset[20];

  key = 1;
  create_table();
  SCTPUT("\n");
  SCTPUT("Center          FWHM    Contin       Flux       EQWT        Error");
  SCTPUT("--------------------------------------------------------------------");
  sprintf(sset, "COLOR=4");
  AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
  AG_WDEF(specXcen - specDx, specXcen + specDx,
    specYcen - specDy, specYcen + specDy);
  AG_SSET(sset);
  AG_SSET("CURSOR = 2");
  xint1 = specXcen;
  yint1 = aux = specYcen;
  AG_SSET("SCALE = 1.5");
  while (key == 1)
  {
    xint1 = xint2;
    yint1 = aux;
    AG_VLOC(&xint1, &yint1, &key, &pix);
    if (key == 1)
    {
      for (i = 0; specX[i] < xint1; i++);
      AG_GTXT(xint1, specY[i], "\\downarro", 2);
      xint2 = xint1;
      yint2 = yint1;
      AG_VLOC(&xint2, &yint2, &key, &pix);
      if (key == 1)
      {
        for (i = 0; specX[i] < xint2; i++);
        AG_GTXT(xint2, specY[i], "\\downarro", 2);
        if (xint1 > xint2)
        {
          aux = xint1;
          xint1 = xint2;
          xint2 = aux;
        }
        aux = yint1;
        for (i = 0; specX[i] < xint1; i++);
        ycen = specY[i] - fit_cont(specX[i]);
        xcen = specX[i];
        specFluxReal = 0.0;
        for (; specX[i] < xint2; i++)
        {
          aux = specY[i] - fit_cont(specX[i]);
          specFluxReal += specStep * aux;
          if (aux * aux > ycen * ycen)
          {
            ycen = aux;
            xcen = specX[i];
          }
        }
        AG_CLS();
        sgauss(ycen, xcen, specStep);
        AG_VDEF("graph_wnd0/n:", 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
        AG_CDEF(specClip[0], specClip[1], specClip[2], specClip[3]);
        AG_WDEF(specXcen - specDx, specXcen + specDx,
              specYcen - specDy, specYcen + specDy);
        AG_SSET(sset);
        AG_SSET("CURSOR = 2");
        AG_SSET("SCALE = 1.5");
      }
    }
  }
  AG_SSET("SCALE = 1.0");
  AG_VUPD();
  AG_CLS();
  TCTCLO(tid);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  print_plot(mode, printer) : print the current ploted frame.
        mode   : LANDSCAPE or PORTRAIT mode
---------------------------------------------------------------------*/
print_plot(mode, printer)
  int             mode;
  char           *printer;
{
  char            sys[60];
  char            dname[20];
  int             iac;
            
  system("rm -f pscrplot.*");
  
  (void) strcpy(dname, printer);
  if (mode == NORMALPRINT)  
    (void) strcat( dname, ":");
  else
    (void) strcat( dname, ".p:");
  AG_VDEF(dname, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
  AG_MRDW(PLOTNAME);
  AG_MRDW("alicel.plt");
  AG_CLS();

  (void) SCKGETC( "SYSCOMS", 1, 20, &iac, sys );
  (void) strcat( sys, printer );
  (void) strcat( sys, " pscrplot.0");
  system(sys);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 print_statistics() : print the current statistics on the Statistics
                      interface.
---------------------------------------------------------------------*/
void 
print_statistics()
{
  FILE           *stat;
  char            line[80];
  double          area, total = 0.0;

  if ((stat = fopen("TMPalice.stat", "w")) != NULL)
  {
    fprintf(stat, "Alice output file\n\n");
    fprintf(stat, "Input file : %s  ", specImageName);
    if (specLineStep > 1)
      fprintf(stat, "Lines %d to %d\n\n", 
              specLineNum, specLineNum + specLineStep - 1);
    else
      fprintf(stat, "Line %d\n\n", specLineNum);
    fprintf(stat,
            "Frame limits \n\txmin: %f  xmax: %f\n\tymin: %f  ymax %f\n\n",
      specXcen - specDx, specXcen + specDx,
      specYcen - specDy, specYcen + specDy);
    fprintf(stat, "Gaussian Values:\n");
    fprintf(stat, "\t\t\tInitial Guess\n\n Component\tAmplitude\t\tPosition\t\tWidth\n");
    read_init_guess();
    for (i = 0; i < gaussNumOfComp; i++)
    {
      sprintf(line, "    %d\t\t%f\t\t%f\t\t%f\n", i + 1, gaussFitValues[i * 3],
        gaussFitValues[i * 3 + 1], gaussFitValues[i * 3 + 2]);
      fprintf(stat, "%s", line);
    }
    fprintf(stat, "\n");
    sprintf(line, "\t\t\tSolution\n\n Component\tAmplitude\t\tPosition\t\tWidth\n");
    fprintf(stat, "%s", line);
    read_fit_values();
    for (i = 0; i < gaussNumOfComp; i++)
    {
      sprintf(line, "    %d\t\t%f\t\t%f\t\t%f\n", i + 1, gaussFitValues[i * 3],
        gaussFitValues[i * 3 + 1], gaussFitValues[i * 3 + 2]);
      fprintf(stat, "%s", line);
    }
    fprintf(stat, "\n\n Statistics\n");
    fprintf(stat, "\t\tGaussian Standard Deviation\n\n");
    fprintf(stat, " Component\tAmplitude\t\tPosition\t\tWidth\n");
    for (i = 0; i < gaussNumOfComp; i++)
    {
      sprintf(line, "    %d\t\t%f\t\t%f\t\t%f\n", i + 1, gaussErrors[i * 3],
          gaussErrors[i * 3 + 1], gaussErrors[i * 3 + 2]);
      fprintf(stat, "%s", line);
    }
    fprintf(stat, "\n RMS : %f\n", fitRms);
    fprintf(stat, "\n\t\tIntegration:\n");
    fprintf(stat, " Component\tflux\t\t\tfwhm\n");
    for (i = 0; i < gaussNumOfComp; i++)
    {
      area = 2.50663 * gaussFitValues[i * 3] * gaussFitValues[i * 3 + 2];
      total += area;
      fprintf(stat, "    %d\t\t%f\t\t%f\n", i + 1, area,
        2.35482 * gaussFitValues[i * 3 + 2]);
    }
    fprintf(stat, "Total flux: \t%f\n", total);
    fclose(stat);
  } else
    printf("Can't open tmp file\n");
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 create_table() : create a table for the integration output.
---------------------------------------------------------------------*/
create_table()
{
  TCTINI(ALICETABLE, F_TRANS, F_O_MODE, 20, 300, &tid);
  TCCINI(tid, D_C_FORMAT, 32, "A32", "", "IDENT", &col_ident);
  TCCINI(tid, D_R4_FORMAT, 1, "F10.3", "", "FLUX_GAUSS", &col_fgauss);
  TCCINI(tid, D_R4_FORMAT, 1, "F10.3", "", "FLUX_LINE", &col_fline);
  TCCINI(tid, D_R4_FORMAT, 1, "F10.3", "", "XCEN", &col_xcen);
  TCCINI(tid, D_R4_FORMAT, 1, "F10.3", "", "XFWHM", &col_xfwhm);
  TCCINI(tid, D_R4_FORMAT, 1, "F10.3", "", "CONT", &col_cont);
  TCCINI(tid, D_R4_FORMAT, 1, "F10.3", "", "EQWT", &col_eqwt);
  TCCINI(tid, D_R4_FORMAT, 1, "F10.3", "", "ERROR", &col_error);
  TCCINI(tid, D_I4_FORMAT, 1, "I8", "", "YSTART", &col_ystart);
  TCCINI(tid, D_I4_FORMAT, 1, "I8", "", "YEND", &col_yend);
  currline = 1;
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
put_table_values() 
---------------------------------------------------------------------*/
put_table_values(ident, fluxg, fluxl, xcen, xfwhm, 
                 cont, eqwt, error, ystart, yend)
  char           *ident;
  float           fluxg, fluxl, xcen, xfwhm, cont, eqwt, error;
  int             ystart, yend;
{
  float           fvalues[20];
  int             ivalues[20];
  int             colnum[20];
  int             column, row, nsort, allcol, allrow;

  TCIGET(tid, &column, &row, &nsort, &allcol, &allrow);
  colnum[0] = col_fgauss;
  colnum[1] = col_fline;
  colnum[2] = col_xcen;
  colnum[3] = col_xfwhm;
  colnum[4] = col_cont;
  colnum[5] = col_eqwt;
  colnum[6] = col_error;
  fvalues[0] = fluxg;
  fvalues[1] = fluxl;
  fvalues[2] = xcen;
  fvalues[3] = xfwhm;
  fvalues[4] = cont;
  fvalues[5] = eqwt;
  fvalues[6] = error;
  TCRWRR(tid, currline, 7, colnum, fvalues);
  colnum[0] = col_ystart;
  colnum[1] = col_yend;
  ivalues[0] = ystart;
  ivalues[1] = yend;
  TCRWRI(tid, currline, 2, colnum, ivalues);
  TCEWRC(tid, currline, col_ident, ident);
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 out_integration() : print the integration output.
---------------------------------------------------------------------*/
out_integration()
{
  int             i;
  char            out[80];
  float           fit_cont();
  double          area, fwhm, eqwt, conta2, xl, xr, aux, total = 0.0;

  read_fit_values();
  SCTPUT("\n");
  SCTPUT("------------------------------------------------------");
  SCTPUT("  Comp\tArea\t\tFWHM\t\tContin");
  xl = gaussFitValues[1] - 3 * gaussFitValues[2];
  xr = gaussFitValues[1] + 3 * gaussFitValues[2];
  for (i = 1; i < gaussNumOfSol; i++)
  {
    if ((aux = gaussFitValues[3 * i + 1] + 3 * gaussFitValues[3 * i + 2]) > xr)
      xr = aux;
    if ((aux = gaussFitValues[3 * i + 1] - 3 * gaussFitValues[3 * i + 2]) < xl)
      xl = aux;
  }
  for (i = 0; specX[i] < xl; i++);
  specFluxReal = 0.0;
  for (; specX[i] < xr; i++)
    specFluxReal += specStep * (specY[i] - fit_cont(specX[i]));
  for (i = 0; i < gaussNumOfSol; i++)
  {
    fwhm = 2.35482 * gaussFitValues[i * 3 + 2];
    conta2 = fit_cont(gaussFitValues[i * 3 + 1]);
    area = 2.50663 * gaussFitValues[i * 3] * gaussFitValues[i * 3 + 2];
    total += area;
    sprintf(out, "     %d\t%8.5f\t%8.5f\t%6.0f", i + 1, area, fwhm, conta2);
    SCTPUT(out);
  }
  sprintf(out, "\n Total area : \t%f", total);
  SCTPUT(out);
  sprintf(out, " Real flux : \t%f\t(x: %g - %g)", specFluxReal, xl, xr);
  SCTPUT(out);
  SCTPUT("------------------------------------------------------\n");
}

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 SQR(sq) 
---------------------------------------------------------------------*/
float 
SQR(sq)
  float           sq;
{
  return (sq * sq);
}