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

  Thomas Szeifert, Anton Malina

  echrebin1d.c 
  
  rebin echelle spectra

*/

/* system includes */

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

/* general Midas includes */

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

/* FEROS specific includes */

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

int main ()
{
  int Method, NpixOut, Npix[2], Npix2[2], kun, knul;
  int dunit;
  int status, actvals, null;	
  int naxis, i, k, tid;
  int dnull, imno, imno2;
  int iparam[1], *Npixo, *ibuff;
  int fitdeg, norder, size;
  float lhcuts[4];
  float *yin, *yout;
  char scaling;
  char inframe[80], intable[80], outframe[80];
  char ident[80], cunit[80], output[80];
  char line[80];
  char cparam[1];
  double *x, *win, *wout;
  double Starti, Stepi;
  double *Starto, *LogStarto, Stepo, vcorr;
  double dparam[2];
  double Start2[2], Step2[2];
  double *dcoef, dbuff, *buffcoef;
  double tmp;

  SCSPRO ("echrebin1d");

  /* read the input parameters */

  status = SCKGETC ("IN_A", 1, 80, &actvals, inframe);
  if (status)
    SCTPUT ("Error while reading IN_A");
  status = SCKGETC ("IN_B", 1, 80, &actvals, intable);
  if (status)
    SCTPUT ("Error while reading IN_B");
  status = SCKGETC ("OUT_A", 1, 80, &actvals, outframe);
  if (status)
    SCTPUT ("Error while reading OUT_A");
  status = SCKRDI ("INPUTI", 1, 1, &actvals, iparam, &kun, &knul);
  status = SCKRDC ("INPUTC", 1, 1, 1, &actvals, cparam, &kun, &knul);
  status = SCKRDD ("INPUTD", 1, 2, &actvals, dparam, &kun, &knul);

  Method = iparam[0];

  Stepo = dparam[0];
  vcorr = 1.0+dparam[1]/LIGHTSPEED;

  scaling = toupper(cparam[0]);

  SCTPUT ("rebin spectrum");
  SCTPUT ("--------------\n");

  switch (scaling)
    {
    case 'I':
      {
	sprintf (line, "rebinning to constant steps in lambda");
	break;
      }
    case 'O':
      {
	sprintf (line, "rebinnig to constant steps in log lambda");
	break;
      }
    default:
      {
	sprintf (line, 
		 "invalid scaling method specified. using linear scaling.");
	scaling = 'I';
	break;
      }
    }
  SCTPUT (line);

  /* open "line-by-line" input frame and read descriptors */

  if (SCFOPN(inframe, D_R4_FORMAT, 0, F_IMA_TYPE, &imno) != 0)
    {
      sprintf(output, "Frame %s invalid...", inframe );
      SCTPUT(output);
      return(255);
    }

  SCDRDI(imno, "NPIX", 1, 2, &actvals, Npix, &dunit, &dnull);
  SCDRDI(imno, "NAXIS", 1, 1, &actvals, &naxis, &dunit, &dnull);
  if (naxis == 1) Npix[1] = 1;
  SCDRDD(imno, "START", 1, 1, &actvals, &Starti, &dunit, &dnull);
  SCDRDD(imno, "STEP", 1, 1, &actvals, &Stepi, &dunit, &dnull);
  SCDRDR(imno, "LHCUTS", 1, 4, &actvals, lhcuts, &dunit, &dnull);
  SCDGETC(imno, "IDENT", 1, 72, &actvals, ident);
  SCDGETC(imno, "CUNIT", 1, 72, &actvals, cunit);

  /* open WLC table to read the dispersion coeficients and the other
     related information */

  TCTOPN(intable, F_IO_MODE, &tid);
  SCDRDI(tid, "FITDEG", 1, 1, &actvals, &fitdeg, &dunit, &dnull);
  SCDRDI(tid, "NORDER", 1, 1, &actvals, &norder, &dunit, &dnull);
  dcoef = dvector (0, fitdeg * norder);
  buffcoef = dvector (1, fitdeg + 1);
  SCDRDD(tid, "DCOEF", 1, fitdeg * norder, &actvals, &dcoef[0], &dunit, &dnull);
  TCTCLO(tid);

  /* compute start and end wavelength and number of pixels of the output */ 

  Starto = dvector(0, norder);
  LogStarto = dvector(0, norder);
  Npixo  = ivector(0, norder);
  NpixOut = 0;
  for (k = 0; k < norder; k++)
    {
      for (i = 1; i <= fitdeg; i++)
	buffcoef[i] = dcoef[k * fitdeg + (i - 1)];

      /* start position (A) */

      Starto[k] = eval_dpoly(Starti, buffcoef, fitdeg); 

      /* correction for barycentric velocity */

      Starto[k] = Starto[k]*vcorr; 
      LogStarto[k] = log(Starto[k]);

      /* start each order at integer value of step */

      Starto[k] = Stepo * (int) (Starto[k]/Stepo) ; 
      LogStarto[k] = Stepo * (int) (LogStarto[k]/Stepo) ; 

      /* end position (A) */

      dbuff = eval_dpoly((double) 
			 (Starti + Npix[0] * Stepi), buffcoef, fitdeg); 

      /* correction for barycentric velocity */

      dbuff = dbuff*vcorr; 
      switch (scaling)
	{
	case 'I':
	  {
	    /* number of pixels */
	    Npixo[k] = (int) ((dbuff - Starto[k]) / Stepo); 
	    break;
	  }
	case 'O':
	  {
	    tmp = log(dbuff) - LogStarto[k];
	    Npixo[k] = (int) (tmp / Stepo);
	    break;
	  }
	}
      if (NpixOut < Npixo[k]) NpixOut = Npixo[k];
    }

  /* allocate memory for the input (x,win,yin) and the output (wout,yout) */

  x    = dvector(0, Npix[0]);
  win  = dvector(0, Npix[0]);
  yin  = vector(0, Npix[0]);

  wout = dvector(0, NpixOut);
  yout = vector(0, NpixOut);

  /* prepare output arrays of the descriptors of the output frame*/

  Npix2[0] = NpixOut;
  Npix2[1] = norder;

  Start2[0] = 1.0;
  Start2[1] = 1.0;

  Step2[0] = Stepo;
  Step2[1] = 1.0;

  /* open output frame with the SCF-interface and write some descriptors */

  size = (NpixOut * norder);
  SCFCRE (outframe, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE, size, &imno2);
  SCDWRI (imno2, "NAXIS", &naxis, 1, 1, &kun);
  SCDWRI (imno2, "NPIX",  Npix2, 1, 2, &kun);
  SCDWRD (imno2, "START", Start2, 1, 2, &kun);
  SCDWRD (imno2, "STEP", Step2, 1, 2, &kun);
  SCDWRC (imno2, "IDENT", 1, ident, 1, 72, &kun);
  SCDWRC (imno2, "CUNIT", 1, cunit, 1, 72, &kun);
  SCDWRR (imno2, "LHCUTS", lhcuts, 1, 4, &kun);

  /* do the real rebin work */

  for (i = 0; i < Npix[0]; i++)
    x[i]  = (double) (Starti + Stepi * i);

  for (k = 0; k < norder; k++)   /* loop over all orders */
    {
      for (i = 0; i < NpixOut; i++)
	{
	  switch (scaling)
	    {
	    case 'I':
	      {
		/* output wavelength */
		wout[i] = (double) (Starto[k] + Stepo * i); 
		break;
	      }
	    case 'O':
	      {
		wout[i] = exp(LogStarto[k] + (double) (Stepo * i));
		break;
	      }
	    }
	}
      for (i = 1; i <= fitdeg; i++)
	{
	  buffcoef[i] = dcoef[k * fitdeg + (i - 1)];
	}

      for (i = 0; i < Npix[0]; i++)
	{

	  /* input wavelength */

	  win[i] = eval_dpoly(x[i], buffcoef, fitdeg); 

	  /* correction for barycentric velocity */

	  win[i] = win[i]*vcorr; 

	}

      /* yin ... input flux */

      SCFGET(imno, k * Npix[0] + 1, Npix[0], &null, (char *) yin); 
      rebin(win, wout, yin, yout, Npix[0], NpixOut, Stepo, Method, scaling);
      SCFPUT(imno2, k * NpixOut + 1, NpixOut , (char *) yout);
    }

  /* write some more descriptors to the output frame to be consistent with
     old echelle context */

  if (scaling == 'I')
    SCDWRD(imno2, "WSTART", Starto, 1, norder, &kun);
  else
    SCDWRD(imno2, "WSTART", LogStarto, 1, norder, &kun);
  SCDWRI(imno2, "NPTOT", Npixo, 1, norder, &kun);
  ibuff = ivector(0, norder);
  for (k = 0; k < norder; k++)
    ibuff[k] = k + 1;

  SCDWRI(imno2, "NORDER", ibuff, 1, norder, &kun);

  /* Close everything and good bye */

  free_dvector(dcoef, 0, fitdeg * norder);
  free_dvector(buffcoef, 1, fitdeg);

  free_dvector(Starto, 0, norder);
  free_dvector(LogStarto, 0, norder);
  free_ivector(Npixo, 0, norder);

  free_dvector(x, 0, Npix[0]);
  free_dvector(win, 0, Npix[0]);
  free_vector(yin, 0, Npix[0]);

  free_dvector(wout, 0, NpixOut);
  free_vector(yout,0, NpixOut);

  SCFCLO(imno2);
  SCFCLO(imno);
  SCSEPI();
  return 0;
}