/* @(#)mossky.c	17.1.1.1 (ESO-IPG) 01/25/02 17:54:49 */
/* @(#)mossky.c       1.3 (ESO-IPG) 2/11/94 12:11:45 */
/* fit sky with median or polynomial                 */
/*                   Otmar Stahl, Nov.28, 1994       */

#include <tbldef.h>
#include <midas_def.h>
#include <ok.h>
#include <math.h>
#include <proto_mutil.h>
#include <proto_nrutil.h>

#define MEDIAN 0
#define POLYNOMIAL 1

#ifdef __STDC__
/*void median2 (float[], float[], int[], int,
  double[], double[], int, int, int, int, double[], double[]);
  void poly2 (float[], float[], int[], int,
  double[], double[], int, int, int, int, double[], double[]);*/
	 	
void median2 (float[], float[], int[],
	      double[], double[], int,  int, int, double[], double[]);
void poly2 (float[], float[], int[], int , double[], double[], int, int, int, double[], double[], double[], double[],float[]);
void copy_frame(float[],float[],int);
#else
void median2 ();
void poly2();
void copy_frame();
#endif

double *A;

main ()

{
  int parm[5],method;
  int i, iav, order;
  int nax, kun, knul, tid, imno1, imno2;
  int col, ncol, nslits, acol, arow, nsort, select;
  int npix[2], null[3], slit[100], upper[100], lower[100];
  int sky_upper, sky_lower;
  int icol[3],ii,npos,k;
  
  char  inframe[60], mos_table[60], win_table[60], outframe[60], mode[11];
  float *inimage, *outimage, cosmic[3];
  char unit[64], inst[72], line[80];

  float moscol[3];
  double *flux, *position, *tmp1, *tmp2;

  double dstart[2], dstep[2];

  /*         get into Midas and get parameters             */

  SCSPRO ("mossky");

  strcpy(unit,""), strcpy(inst,"");

  SCKGETC ("IN_A", 1, 60, &iav, inframe);
  SCKGETC ("IN_B", 1, 60, &iav, mos_table);
  SCKGETC ("IN_C", 1, 60, &iav, win_table);
  SCKGETC ("OUT_A", 1, 60, &iav, outframe);
  SCKGETC ("INPUTC", 1 , 10, &iav, mode);
  SCKRDI ("INPUTI", 1, 4, &iav, parm, &kun, &knul);
  SCKRDR ("INPUTR", 1, 3, &iav, cosmic, &kun, &knul);

  order = parm[0];

  A = dvector(1, order);
  
  SCTPUT ("\n ----------------------- ");
  sprintf (line, "Input image:         %s ", inframe);
  SCTPUT (line);
  sprintf (line, "Input table:         %s ", win_table);
  SCTPUT (line);
  sprintf (line, "Output image:        %s\n ", outframe);
  SCTPUT (line);
  SCTPUT ("input parameters:\n ");

  method = MEDIAN;
  sprintf(line,"Fitting method: median");

  if (!strncmp(mode,"POL",3) || !strncmp(mode,"pol",3))
    {
      method = POLYNOMIAL;
      sprintf(line,"Fitting method: polynomial");
    }	

  SCTPUT (line);
  sprintf (line, "order of fit:        %i", order);
  SCTPUT (line);

/*          map input file                              */

  SCIGET (inframe, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE,
	  2, &nax, npix, dstart, dstep, inst, unit, (char **)&inimage, &imno1);

  SCIPUT (outframe, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	  nax, npix, dstart, dstep, inst, unit, (char **)&outimage, &imno2);

  copy_frame(inimage,outimage,npix[0]*npix[1]);

  /*      arrays for intermediate data                                */

  position = (double *) osmmget (npix[1] * sizeof (double));
  flux     = (double *) osmmget (npix[1] * sizeof (double));
  tmp1     = (double *) osmmget (npix[1] * sizeof (double));
  tmp2     = (double *) osmmget (npix[1] * sizeof (double));


/*-------read table MOS------------------------------------------------*/

  TCTOPN (mos_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]);
  
/* 1.step-----------read table MOS------------------------------------*/
  /*copy selected slitlets to slit[0,..,nslits-1] , */
  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]-dstart[1])/dstep[1]) + 1;
        upper[nslits] =  (int) ((moscol[2]-dstart[1])/dstep[1]) + 1;
        nslits++;
      }
    }

  TCTCLO (tid);

/*           read table WINDOWS                           */

  TCTOPN (win_table, F_I_MODE, &tid);
  TCIGET (tid, &col, &ncol, &nsort, &acol, &arow);
 
  TCCSER (tid, ":Sky_Slit", &icol[0]);
  TCCSER (tid, ":Sky_Strt", &icol[1]);
  TCCSER (tid, ":Sky_End", &icol[2]);

/* 2.step-----------loop through slitlets, search sky windows-------------*/

  SCTPUT ("\n ----------------------- ");
  SCTPUT (" slit: rows for sky:");

  /* for all slitlets*/
  for (ii = 0; ii < nslits; ii++) {
    npos = 0;
    /*for all selected skypositions*/
    for (i = 1; i <= ncol; i++) {
      TCSGET(tid, i, &select);
      if (select) {
	TCRRDR (tid, i, 3, icol, moscol, null);
	if ((int) (moscol[0]) == slit[ii] &&  ! null[0] ) {
	  sky_lower =  (int) ((moscol[1]-dstart[1])/dstep[1]) + 1;
	  sky_upper =  (int) ((moscol[2]-dstart[1])/dstep[1]) + 1;

	  /* 3.step----------produce array of sky positions----------------*/
	  /* position[0..npos-1]*/
	  for (k=sky_lower; k<=sky_upper;k++) {
	    position[npos] =  (k-1)*dstep[1] + dstart[1];
	    npos++;
	  }
	}
      }
    }

    if (npos > 0) {
      sprintf(line,"%4i       %4i",slit[ii],npos);
      SCTPUT(line);
      switch (method)
	{
	case MEDIAN:
	  /* median2 (inimage, outimage, npix, order, position, flux, 
	     npos,ii,lower[ii-1],upper[ii-1],dstart,dstep); */
	  median2 (inimage, outimage, npix, position, flux, 
		   npos,lower[ii],upper[ii],dstart,dstep);
	  break;

	case POLYNOMIAL:
	  /* poly2 (inimage, outimage, npix, order, position, flux, 
	     npos,ii,lower[ii-1],upper[ii-1],dstart,dstep); */
	  poly2 (inimage, outimage, npix, order, position, flux, 
		 npos,lower[ii],upper[ii],dstart,dstep,tmp1,tmp2,cosmic);
	  break;
	}
    }
    else {
      sprintf(line,"%4i       no sky ",ii+1);
      SCTPUT(line);
    }
    
  }

/*          epilogue ...                 */


  TCTCLO (tid);
  osmmfree((char *) position);
  osmmfree((char *) flux);
  osmmfree((char *) tmp1);
  osmmfree((char *) tmp2);
  free_dvector(A,1,order);
  SCSEPI ();
}

/* #ifdef __STDC__
void
poly2 (float inframe[], float outframe[], int npix[],
	 int order,  double position[], double flux[], int npos, int slit,
	 int lower, int upper, double dstart[], double dstep[])
#else
void
poly2 (inframe, outframe, npix, order, position, flux, npos,slit,
       lower,upper,dstart,dstep)
  float inframe[],outframe[];
  int npix[],order,npos,slit,lower,upper;
  double dstart[], dstep[], position[], flux[];
#endif */

#ifdef __STDC__
void
poly2 (float inframe[], float outframe[], int npix[],
	 int order,  double position[], double flux[], int npos, 
	 int lower, int upper, double dstart[], double dstep[], double tmp1[],
         double tmp2[], float cosmic[])
#else
void
poly2 (inframe, outframe, npix, order, position, flux, npos,
       lower,upper,dstart,dstep,tmp1,tmp2,cosmic)
  float inframe[],outframe[],cosmic[];
  int npix[],order,npos,lower,upper;
  double dstart[], dstep[], position[], flux[], tmp1[], tmp2[];
#endif
     
{
  int i, j, k, jj,nsel;
  char line[80];
  double xout;
  float med,limit;

  /*for all pixel - x-axis*/
  for (k = 0; k < npix[0]; k++) 
    {
    /*for all pixel within sky-position --  y-axis */
    for (j = 0; j < npos; j++)
      {
	jj = ((int)((position[j]-dstart[1])/dstep[1]+0.5) * npix[0] + k);
	flux[j] = (double) inframe[jj];
	tmp1[j+1] = (double) inframe[jj];
      }

    if ( (int) cosmic[2] > 0) /* cosmic cleaning & polynom fit */
      {
      /*1              median of flux and           */
      if  (   (npos+1)/2 ==  npos/2 ){      /*even*/
	med =  (float) select_pos (npos/2 , npos, tmp1) / 2. +
	  (float) select_pos (npos/2+1 , npos, tmp1) / 2.  ;}
      else          /*odd*/
	med =  (float) select_pos ((npos+1)/2, npos, tmp1);
      /*2              detect cosmics and           */ 
      /*3              pik selected clean position  */
      nsel = 0;
      if (med > 0.)
	limit = (float) (cosmic[2]*(sqrt(med)/sqrt(cosmic[1]) + cosmic[0]/cosmic[1]));
      for (j = 0; j < npos; j++)
	{
	  if (med > 0.)
	    {
	      if ((float) flux[j] < med + limit && 
		  (float) flux[j] > med - limit) 
		{
		  tmp1[nsel] = flux[j];
		  tmp2[nsel] = position[j];  
		  nsel++;
		}
	    }
	}
      /*4              polynomial fit - after cleaning */
      if (nsel > order)
	lfit(tmp2,tmp1,nsel,A,order,dpoly);
      else
	lfit(position,flux,npos,A,order,dpoly);
      }
    else
      {
	/*              polynomial fit - no clean function */
	lfit(position,flux,npos,A,order,dpoly);
      }

    /*                  prepare output frame           */
    for (j = lower; j <= upper; j++)
      {
	xout = dstep[1] * (j-1) + dstart[1];
	jj = (j-1) * npix[0] + k;
	outframe[jj] = (float) eval_dpoly(xout,A,order);
      }
    }

}

/*-------------------------------------------------------------------------*/
/* #ifdef __STDC__
void
median2 (float inframe[], float outframe[], int npix[],
	 int order,  double position[], double flux[], int npos, int slit,
	 int lower, int upper, double dstart[], double dstep[])
#else
void
median2 (inframe, outframe, npix, order, position, flux, npos,slit,
	 lower,upper,dstart,dstep)
  float inframe[],outframe[];
  int npix[],order,npos,slit,lower,upper;
  double dstart[2],dstep[2],position[],flux[];   
#endif */

#ifdef __STDC__
void
median2 (float inframe[], float outframe[], int npix[],
	 double position[], double flux[], int npos,
	 int lower, int upper, double dstart[], double dstep[])
#else
void
median2 (inframe, outframe, npix, position, flux, npos,
	 lower,upper,dstart,dstep)
  float inframe[],outframe[];
  int npix[],npos,lower,upper;
  double dstart[],dstep[],position[],flux[];   
#endif

{
  float med;
  int j, k, jj;
  char line[80];

  /*for all pixel - x-axis*/
  for (k = 0; k < npix[0]; k++) {
    for (j = 0; j < npos; j++)
      {
	/*next integer pixel number*/
	jj = ((int)((position[j]-dstart[1])/dstep[1]+0.5) * npix[0] + k );
	/*NR starts with flux[1]*/
	flux[j+1] = (double) inframe[jj];
      }
    /*              fast median with select                 */
    if  (   (npos+1)/2 ==  npos/2 ){
      /*even*/
      med =  (float) select_pos (npos/2 , npos, flux) / 2. +
	(float) select_pos (npos/2+1 , npos, flux) / 2.  ;}
    else    
      /*odd*/
      med =  (float) select_pos ((npos+1)/2, npos, flux);
    /*output column*/
    for (j = lower; j <= upper; j++)
      {	
	jj = (j-1) * npix[0] + k;
	outframe[jj] =  med;
      }
  }
}

#ifdef __STDC__
void copy_frame(float inframe[], float outframe[], int npix)
#else
void copy_frame(inframe, outframe, npix)
float inframe[],outframe[];
int npix;
#endif
{
int j;

  /* init outframe to inframe                                */

  for (j = 0; j < npix; j++)
    {
	  outframe[j] = inframe[j];
    }
}