/*
 .IDENTIFICATION
   subroutine SORTMED
 .AUTHOR
   K. Banse                  	ESO - Garching
 
 .KEYWORDS
   sorting
 
 .PURPOSE
   sort an array to get the median value
 
 .ALGORITHM
   we only sort fully until index of median - then we insert
   for sorting we use the Heapsort algorithm from "Numerical Recipes", page 231
 
 .INPUT/OUTPUT
 
   call as STAT = SORTMED(RA,ZBINS,N,LL,VALU)
 
   input par:
   ZBINS:	R*4 array	low, high excess bins
 				if low >= high, no excess bins
   NUM:		I*4		size of array RA
   LL:		I*4		index of array element to be returned
 
   in/output par:
   RA:		R*4 array	array to be sorted
 
   output par:
   VALU:	R*4 		value = RA(LL) of sorted array RA
   STATUS:	I*4		return status, = 0 is o.k.
 
 .VERSIONS
   1.00		pulled out from `smosub3.for'
   1.50         transl.  FORTRAN to C, makes use of sortr() [swolf@eso.org]
*/

/* PROTOTYPE 
   void sortr(float *rka,int ndim);	does NOT work on HP */

void sortr();

int sortmed(ra,zbins,num,ll,valu)
int   num, ll;
float *ra, *zbins, *valu;
{
  int   ntrue, lltrue, l, ir, ir1, j, j1, i, i1, status, k;
  float rra, zl, zh;

  status = 0;

  /*  first check for excess bins */
  if (zbins[2] > zbins[1])
    {
      zl = zbins[1];
      zh = zbins[2];
      i = 0;
      for (l=0;l<num;l++)
	if ((ra[l] >= zl) && (ra[l] <= zh)) {
	  ra[i] = ra[l];
	  i++;
	}
      if (i > 3)
	{
	  ntrue = i ;
	  lltrue = (ntrue+1)/2;
	}
      else if (i < 1)
	{
	  status = -1;
	  return(status);
	}
      else if (i <= 2)
	{
	  *valu = ra[0];
	  return(status);
	}
      else                     /* i == 3 */
	{
	  *valu = ra[1];
	  return(status);
	}
    }
  else
    {
      ntrue = num;
      lltrue = ll;
    }
 
  sortr(ra,ntrue);
  *valu = ra[lltrue];
}



void sortr(rka,ndim)

/*
  ndim  = dimension of array rka for sorting
  but we pass the arrays with 1 element in front, so that the algorithm
  sorts from [1] -> [ndim]
  this is the Heapsort algorithm from "Numerical Recipes", page 231

  In case you an array has to be sorted starting at index 0 then you have
  to pass the reduced adress to rka:
                float *fp, farr;
		fp = farr-1;
		sortit(fp,ndim);
  
*/
  
float *rka;
int   ndim;

{
register int  m, j, i;
int k, ntrue=ndim;
float  zka;
  

m = (ndim>>1) + 1; 		     /* ndim/2 + 1 */
  
for (;;)
   {
   if (m > 1)			     /* still hiring */
      zka = rka[--m];
   else			             /* in retirement + promotion phase */
      {
   /* for (k=1;k<=ntrue;k++) printf("%f ",rka[k]); printf("\nir=%d\n",ndim); */
      zka = rka[ndim];	             /* clear a space at end of array rka */
      rka[ndim] = rka[1];            /* retire the top of the heap into it */
      if (--ndim == 1)	             /* done with last promotion */
         {
         rka[1] = zka;	             /* the least competent guy of all... */
         return;		     /* that's it folks  */
         }
      }
     
/* in hiring as well as in promotion phase */
/* here set up to shift down zka to its right level */

   i = m;	
   j = m << 1;                       /* in FORTRAN: j = m + m */
  

   while (j <= ndim)	
      {
      if (j < ndim)
         {
         if (rka[j] < rka[j+1]) j++;
         }
   
      if (zka < rka[j]) 	     /* demote zka */
         {
         rka[i] = rka[j];
         i = j;
         j = j << 1;		     /* in FORTRAN: j = j + j */
         }
      else
         j = ndim + 1;
 
      }
   rka[i] = zka;		     /* put zka into its slot */
   }
}