/* @(#)necregr.c	17.1.1.1 (ES0-DMD) 01/25/02 17:51:31 */
/*===========================================================================
  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
===========================================================================*/

/* Program  : echregr.c                          */
/* Author   : P. Ballester  -  ESO Garching      */
/* Date     : 17.03.92                           */
/*                                               */
/* Purpose  :                                    */
/*   Performs a linear fit on independent orders */
/*   and writes downs the rms in a table         */
/*                                               */
/* Input    :                                    */
/*      IN _A     Name of order  table           */
/*      INPUTI(1) Number of orders               */
/*      OUT_A     Name of output table           */
/*                                               */

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

main()

{
      int       kunit;
      char      outtab[TEXT_LEN], ordtab[TEXT_LEN];

      int  order, nb_order, present_order, order_nb, row;
      int  rmscol, ordcol;
      int  order_col, xcol, ycol;
      int  ncol, nrow, nsort, allcol, allrow;

      int  select, null;
      int  actvals, tid, tidord;

      double    x=0., y=0., cnt, sx, sy, sx2, sxy, sy2;
      double    det, a, b, rms;

      SCSPRO("echregr");

      SCKGETC ("IN_A", 1, 60,  &actvals, ordtab);
      SCKGETC ("OUT_A", 1, 60, &actvals, outtab);

      SCKRDI ("INPUTI", 1, 1, &actvals, &nb_order, &kunit, &null);

      TCTOPN (ordtab, F_I_MODE, &tidord);
      TCIGET (tidord, &ncol, &nrow, &nsort, &allcol, &allrow);
      TCCSER (tidord, "ORDER", &order_col);
      TCCSER (tidord, "X",     &xcol);
      TCCSER (tidord, "Y",     &ycol);

      TCTINI (outtab, F_TRANS, F_IO_MODE, 3, 100, &tid);

      TCCINI (tid, D_R4_FORMAT, 1, "I6",   "  ", "ORDER", &ordcol);
      TCCINI (tid, D_R4_FORMAT, 1, "F12.4", "  ", "RMS",   &rmscol);

      row = 1;

      for (order=1; order<=nb_order; order++) {

          cnt=0., sx=0., sy=0., sx2=0., sxy=0., sy2 = 0.;
          /* Reads present_order to start the while ( */
          TCERDI(tidord, row,   order_col, &order_nb, &null);
          present_order = order_nb;

          while (present_order == order_nb) {

            TCSGET(tidord, row, &select);

            if (select) {
            TCERDD(tidord, row,   xcol,  &x,  &null);
            TCERDD(tidord, row,   ycol,  &y,  &null);
            cnt += 1., sx += x, sy += y, sx2 += x*x, sy2 += y*y, sxy += x*y;
  	    }

            if (row >= nrow) break;
            TCERDI(tidord, ++row,   order_col, &present_order, &null);

	    }

          if (cnt >= 3) {
            det = cnt*sx2 - sx*sx;
            a   = (sy*sx2 - sx*sxy)/det;
            b   = (cnt*sxy - sx*sy)/det;
            rms = (sy2 - a*a*cnt - 2.*b*a*sx - b*b*sx2)/cnt;
            if (rms < 0. && rms > -0.05) rms = 0.;
            rms = sqrt(rms);
	  }
          else    rms = 99999.;

          /* printf("a %f  b %f  rms %f\n",a,b,rms); */

          TCEWRI(tid, order, ordcol, &order_nb);
          TCEWRD(tid, order, rmscol, &rms);

      }


      /* Close everything and good bye. */
      TCTCLO(tidord), TCTCLO(tid), SCSEPI();

    }