/* @(#)opt_ext.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:59 */
/* 

   Otmar Stahl, Anton Malina

   opt_ext.c  

   optimal extraction of spectra 
   based on Horne algorithm

*/

/* system includes */

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

/* general Midas includes */

#include <midas_def.h>

/* FEROS specific includes */

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

float opt_ext 
#ifdef __STDC__ 
( 
 float profile[], float image[], float variance[], int number, int niter, 
 float v0, float gain, float thres,
 float ea, int DoOptext, float *mask_row,
 int masked_counter[], int total_masked
) 
#else 
     ( 
      profile, image, variance, number, niter, v0, gain, thres,
      ea, DoOptext, mask_row, masked_counter, total_masked
      ) 
     float profile[], image[], variance[]; 
     int niter, number, DoOptext; 
     float v0, gain, thres,  ea; 
     float *mask_row; 
     int masked_counter[], total_masked;
#endif

/* extract weighted mean, spatial profile given, clip cosmics */
{
  
  int i, j, kk;
  int something_masked;
  float spec, var, diff, diff2, diffbyimg, ratio;
  float rmax, sum1, sum2, sum3, mp, xmax, spec_opt;
  float *mask;
  
  /* input parameters */

  spec = 0.0;
  var = 0.0;
  mask = vector(0, number);
  xmax = ea;

  for (j = 0; j < number; j++)
    {
      /* initial variance estimate */
      variance[j] = v0 + fabs(image[j]) / gain; 
      /* extract standard spectrum */
      spec = spec + image[j];                   
      /* variance of standard spectrum */ 
      var = var + variance[j];                  
      /* mask = 1 for all pixels */
      mask[j] = 1.0;                            
    } 

  if(niter <= 0) /* return with standard extraction */
    {
      free_vector(mask, 0, number);
      return spec; 
    }

  /* otherwise do optimum extraction */

  i = 0;

  do
    {
      ++i;
      for (j = 0; j < number; j++)
	
	/* revise variance estimates */

	variance[j] = v0 + fabs(spec * profile[j]) / gain; 
      
      rmax = 1.0;
      
      kk = -1;
      something_masked = FALSE;
      for (j = 0; j < number; j++)
	{
	  diff = image[j] - spec * profile[j];
	  diffbyimg = fabs(diff / image[j] * mask[j]);
	  
	  diff2 = SQUARE(diff) * mask[j];
	  ratio = (diff2 / (thres * variance[j]));
	  
	  /* && (diff > 0.0) cosmics always make higher values*/
	  
	  if ((ratio > rmax) )
	    {
	      rmax = ratio;
	      if (diffbyimg > xmax)
		{
		  kk = j;
		  something_masked = TRUE;
		}
	    }
	}
      
      if (kk >= 0)
	{
	  mask[kk] = 0.0;
	  masked_counter[kk]++;
	  total_masked++;
	}
      
      sum1 = sum2 = sum3 = 0.0;
      for (j = 0; j < number; j++)
	{
	  mp = mask[j] * profile[j];
	  sum1 += mp * image[j] / variance[j];
	  sum2 += mp * profile[j] / variance[j];
	  sum3 += mp;
	}
      /* optimal extracted spectrum */
      
      spec = sum1 / sum2; 
      
      /* variance of optimal extracted spectrum */
      
      var = (float) sqrt( (double) (sum3 / sum2)); 
  
    } while (something_masked && (i < niter));
  
  /* mask cosmics only */

  if (!DoOptext) 
    {
      spec_opt = spec;
      
      spec  = var = 0.0;
      for (j = 0; j < number; j++)
	{
	  if (mask[j] !=0 )
	    {
	      spec += image[j];
	    }
	  else 
	    {
	      spec += profile[j]*spec_opt;
	    }
	}
    } 
  
  for(j = 0; j < number; j++)
    {
      mask_row[j] = mask[j];
    }
  
  free_vector(mask, 0, number);
  return spec;
}