/*      FILE:  vec_median.c
 *   PURPOSE:  
 *
 *  Determine the robust median and average deviation of the pixel vector pxv
 *  of size px_n using a temporary pixel vector tv of size px_n.  t_n is the
 *  number of pixels used in the determination of the robust median.
 *
 *    AUTHOR:  Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE:  ANSI C
 *      DATE:  2001JUL03
 * COPYRIGHT:  (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 */

#include "mx.h"
#include "qdphot.h"

#ifndef DBL_MAX
#define DBL_MAX 0.99E307
#endif

void qdphot5_VecdRobustMdn_vf6 (
  double *pxv,
  int px_n,
  double *tv,
  int *median_n_,
  double *median_,
  double *adev_
){

  char mxfunc[] = "qdphot5_VecdRobustMdn_vf6";

  double median;
  double adev;

  int t_n;
  int t_n_old;
  double t_mdn;
  double t_mdn_old;
  double t_adv;
  double t_adv_old;
  double t_limit;
  double t_px;
  double t_mean;
  double t_min;
  int j;
  int jj;
  int m;
  int k;

  t_n = px_n;
  t_n_old = t_n;

  t_mdn_old = pxv[0];
  t_adv_old = DBL_MAX;

  for (j=0; j<t_n; ++j) {
    tv[j] = pxv[j];
  }

  t_mdn = 0;
  t_adv = 0;
  t_limit = DBL_MAX;

  for (k=0; k<=5; ++k) {

    jj = t_n;
    t_n = 0;
    t_mean = 0 ;


    for (j=0; j<jj; ++j) {
      if ( fabs(tv[j] - t_mdn) < t_limit ) {
         t_px = tv[j];
         t_mean = t_mean + t_px;
         tv[t_n] = t_px;
         t_n = t_n + 1;
      }
    }

    if(t_n==0) goto done;

    t_mean = t_mean / t_n;

    /* sort tv */
    /* 
     * SPP version used the following SPP routine: 
     *
     *   call asrtr (tv, tv, t_n) 
     */
    /* CVOS: */
    c_asrtd (tv, tv, &t_n);


    /* minimum value is the 1s element */
    t_min = tv[0];

    /* determine the median value of tv */
    m = t_n / 2;
    if ((m*2)==t_n) {                 /* t_n is even */
      t_mdn = (tv[m] + tv[m-1])/2.0;  /* .. so average these values */
    } else {                          /* t_n is odd */
      t_mdn = tv[m];                  /* .. so take the middle value */
    }

    /* determine the average (mean) deviation */
    t_adv = 0;
    for (j=0; j<t_n; ++j) {
      t_adv = t_adv + fabs(tv[j] - t_mean);
    }
    t_adv = t_adv / t_n;

    /* set t_limit for the next iteration */
    t_limit = 3*t_adv;

    if ( 
      (t_n == t_n_old) 
      && (t_mdn == t_mdn_old) 
      && (t_adv == t_adv_old) 
    ) { 
      goto done;
    } else {
      t_n_old = t_n;
      t_mdn_old = t_mdn;
      t_adv_old = t_adv;
    }

  } /* for k */

done:
  median = t_mdn;
  adev = t_adv;

bye:
  *median_n_ = t_n;
  *median_ = median;
  *adev_ = adev;
  return;

}
/* end-of-file */