/* @(#)mosnorm.c	17.1.1.1 (ESO-IPG) 01/25/02 17:54:48 */
/* @(#)mosnorm.c      1.2 (ESO-IPG) 2/11/94 12:11:44 */

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

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

#define max(A,B) ((A)  > (B) ? (A) : (B))
#define min(A,B) ((A)  < (B) ? (A) : (B))

#define MAXSLITS 100

#define POLYNOM 0
#define MEDIAN 1

int Method, Npix[2], Order;
double Start[2], Step[2];

#ifdef __STDC__
void normalize (float[], float[], float[], float[], int[], int[],
		int[], int);
void smooth (float[], float[], int, int);
#else
void normalize ();
void smooth ();
#endif

main ()

{
  int parm[5];
  int i, iav;
  int nax, kun, knul, tid, imno1, imno2;
  int col, ncol, nslits, acol, arow, nsort, select;
  int null[3];
  int slit[MAXSLITS], upper[MAXSLITS], lower[MAXSLITS];

  int icol[3];

  char inframe[60], table[60], outframe[60], method[20];
  float *inimage, *outimage;
  char unit[64], inst[72], line[80];

  float moscol[3];
  float *buff, *buff2;

/*         get into Midas and get parameters             */

  SCSPRO ("mosnorm");

  strcpy(unit,""), strcpy(inst,""); /* Purify */

  SCKGETC ("IN_A", 1, 60, &iav, inframe);
  SCKGETC ("IN_B", 1, 60, &iav, table);
  SCKGETC ("OUT_A", 1, 60, &iav, outframe);
  SCKGETC ("INPUTC", 1, 20, &iav, method);
  SCKRDI ("INPUTI", 1, 4, &iav, parm, &kun, &knul);

  Order = parm[0];

  SCTPUT ("\n ----------------------- ");
  sprintf (line, "Input image:         %s ", inframe);
  SCTPUT (line);
  sprintf (line, "Input table:         %s ", table);
  SCTPUT (line);
  sprintf (line, "Output image:        %s ", outframe);
  SCTPUT (line);
  SCTPUT ("\nInput parameters: \n");

  if (!strncmp (method, "MED", 3) || !strncmp (method, "med", 3))
    {
      Method = MEDIAN;
      sprintf (line, "Smoothing method: Median");
      SCTPUT (line);
      sprintf (line, "Window width:        %i\n", Order);
      SCTPUT (line);
    }
  else
    {
      Method = POLYNOM;
      sprintf (line, "Smoothing method: Polynomial");
      SCTPUT (line);
      sprintf (line, "Order of fit:        %i\n", Order);
      SCTPUT (line);
    }

/*          map input and output file                              */

  SCIGET (inframe, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE,
	  2, &nax, Npix, Start, Step, inst, unit, (char **)&inimage, &imno1);

  SCIPUT (outframe, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	  nax, Npix, Start, Step, inst, unit, (char **)&outimage, &imno2);

/*           read table MOS                            */

  TCTOPN (table, F_I_MODE, &tid);
  TCIGET (tid, &col, &ncol, &nsort, &acol, &arow);
  TCCSER (tid, ":slit",&icol[0]);
  TCCSER (tid, ":ystart",&icol[1]);
  TCCSER (tid, ":yend",&icol[2]);

  nslits = 0;
  for (i = 1; i <= ncol; i++)
    {
      TCSGET(tid, i, &select);
      if (select) {
	TCRRDR (tid, i, 3, icol, moscol, null);
	slit[nslits] = (int) moscol[0];
        lower[nslits] =  (int) ((moscol[1]-Start[1])/Step[1]) + 1;
        upper[nslits] =  (int) ((moscol[2]-Start[1])/Step[1]) + 1;
	nslits++;
      }
    }

  TCTCLO (tid);

/*              reserve space for arrays                */

  buff = (float *) osmmget (Npix[0] * sizeof (float));
  buff2 = (float *) osmmget (Npix[0] * sizeof (float));

  normalize (inimage, outimage, buff, buff2, slit, upper, lower, nslits);

  osmmfree((char *) buff);
  osmmfree((char *) buff2);

  /*          epilogue ...                 */

  SCSEPI ();

}

#ifdef __STDC__
void
normalize (float inframe[], float outframe[], float rval[],
	   float yval[], int slit[], int upper[], int lower[],
	   int nslits)
#else
void
normalize ( inframe, outframe, rval, yval, slit, upper, lower, nslits)
  float inframe[],outframe[],rval[],yval[];
  int   slit[],upper[],lower[],nslits;
#endif

{

  float nrow;
  int i, nslit, j, k, jj;
  char line[80];

  SCTPUT (" slit no. ");

/* init outframe  to zero                                  */

  for (j = 0; j < Npix[1]; j++)
    {
      for (i = 0; i < Npix[0]; i++)
	{
	  jj = j * Npix[0] + i;
	  outframe[jj] = 0.0;
	}
    }


/*          loop over all slitlets                         */

  for (nslit = 0; nslit < nslits ; nslit++) 
    {
      sprintf (line, "    %4i", slit[nslit]);
      SCTPUT (line);

/*        average all  rows in one slitlet                 */

      for (j = 0; j < Npix[0]; j++)
	rval[j] = 0.0;
      for (j = lower[nslit] - 1; j < upper[nslit]; j++)
	{
	  for (k = 0; k < Npix[0]; k++)
	    {
	      jj = j * Npix[0] + k;
	      rval[k] = rval[k] + inframe[jj];
	    }
	}

      nrow = upper[nslit] - lower[nslit] + 1.0;
      for (j = 0; j < Npix[0]; j++)
	rval[j] = rval[j] / nrow;

      switch (Method)
	{

	case MEDIAN:

	  smooth (rval, yval, Npix[0], Order);
	  break;

	case POLYNOM:

	  fit_poly (rval, yval, 1.0, 1.0, Npix[0], Order);
	  break;

	}

/*      now loop through the rows of one slitlet and normalize       */

      for (j = lower[nslit] -1 ; j < upper[nslit]; j++)
	{
	  for (k = 0; k < Npix[0]; k++)
	    {
	      jj = j * Npix[0] + k;
	      outframe[jj] = inframe[jj] / yval[k];
	    }
	}
    }

  SCTPUT (" ----------------------- ");
}

#ifdef __STDC__
void
smooth (float rval[], float yval[], int npix, int order)
#else
void
smooth ( rval, yval, npix, order )
  float rval[],yval[];
  int   npix,order;
#endif

/*      median filtering of one row: rval in, yval out  */

{

  int j, first, last;

  first = (order - 1) / 2;
  last = npix - first - 1;

  for (j = first; j <= last; j++)
    {
      yval[j] = pik_median (order, &rval[j - 2]);

    }

  for (j = 0; j < first; j++)
    yval[j] = yval[first];
  for (j = last + 1; j < npix; j++)
    yval[j] = yval[last];
}