/* @(#)sort.c	17.1.1.1 (ES0-DMD) 01/25/02 17:49:08 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* @(#)sort.c	17.1.1.1 (ESO-DAG) 01/25/02 17:49:08
++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENTIFICATION: subroutine sort
.KEYWORDS:       sorting heap sort
.PURPOSE:        to sort an array 
.ALGORITHM:      Heapsort algorithm from "Numerical Recipes", page 231
.VERSION:        940221 RHW Modification for mosacing
--------------------------------------------------*/
/*
 * Define _POSIX_SOURCE to indicate
 * that this is a POSIX program
 */
#define  _POSIX_SOURCE 1

/*
 * definition of the used functions
 */
#include <midas_def.h>
#include <ccd_def.h>

/*
 * here start the code of the function
 */

void sorti( n, ra )
int   n;
int   *ra;

{
register int indx;

if ( n < 2 ) return;

indx = n - 1;
do { register int jj, ii = -1;

     for ( jj = 0; jj < indx; jj++ )
         { if ( ra[jj] > ra[jj+1] )
              { int itemp = ra[jj];

                ra[jj]   = ra[jj+1];
                ra[jj+1] = itemp;
                ii = jj;
              }
         }
     indx = ii;
   } 
while ( indx != -1 );

return;
}

/*

*/
void sortr( n, ra )
int   n;
float   *ra;

{
register int indx;

if ( n < 2 ) return;

indx = n - 1;
do { register int jj, ii = -1;

     for ( jj = 0; jj < indx; jj++ )
         { if ( ra[jj] > ra[jj+1] )
              { float rtemp = ra[jj];

                ra[jj]   = ra[jj+1];
                ra[jj+1] = rtemp;
                ii = jj;
              }
         }
     indx = ii;
   } 
while ( indx != -1 );

return;
}

/*

*/
void sortd( n, da )
int     n;
double  *da;

{
register int indx;

if ( n < 2 ) return;

indx = n - 1;
do { register int jj, ii = -1;

     for ( jj = 0; jj < indx; jj++ )
         { if ( da[jj] > da[jj+1] )
              { double dtemp = da[jj];

                da[jj]   = da[jj+1];
                da[jj+1] = dtemp;
                ii = jj;
              }
         }
     indx = ii;
   } 
while ( indx != -1 );
return;
}

/*

*/
sortii( n, xval, yval )
int   n;
int   *xval, *yval;

{
register int indx;

if ( n < 2 ) return;

indx = n - 1;
do { register int jj, ii = -1;

     for ( jj = 0; jj < indx; jj++ )
         { if ( xval[jj] > xval[jj+1] )
              { float xtemp = xval[jj],
                      ytemp = yval[jj];

                xval[jj]   = xval[jj+1];
                yval[jj]   = yval[jj+1];
                xval[jj+1] = xtemp;
                yval[jj+1] = ytemp;
                ii = jj;
              }
         }
     indx = ii;
   } 
while ( indx != -1 );
return;
}

/*

*/
sorti2(ndim,rka,ikb)
/*
  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
*/
  
int       ndim;
int       *rka;
int       *ikb;

{
int       m, ir, j, i;
int       zka;
int       jkb;

m = ndim/2 + 1;
ir = ndim;

for (;;)
   {
   if (m > 1)                   /* still hiring */
      {
      zka = rka[--m];
      jkb = ikb[--m];
      }
   else                 /* in retirement + promotion phase */
      {
      zka = rka[ir];    /* clear a space at end of array rka */
      jkb = ikb[ir];
      rka[ir] = rka[1]; /* retire the top of the heap into it */
      ikb[ir] = ikb[1];
      if (--ir == 1)    /* done with last promotion */
         {
         rka[1] = zka;  /* the least competent guy of all... */
         ikb[1] = jkb;
         return;                /* that's it folks  */
         }
      }
     
/* in hiring as well as in promotion phase */
/* here set up to sift down zka to its right level */

   i = m;
   j = m << 1;                          /* in FORTRAN: j = m + m */
  
   while (j <= ir)
      {
      if (j < ir)
         {
         if (rka[j] < rka[j+1]) j++;
         }
   
      if (zka < rka[j])                 /* demote zka */
         {
         rka[i] = rka[j];
         ikb[i] = ikb[j];
         i = j;
         j = j << 1;                    /* in FORTRAN: j = j + j */
         }
      else
         j = ir + 1;
 
      }

   rka[i] = zka;                /* put zka into its slot */
   ikb[i] = jkb;
   }

}

/*

*/
void sortri( n, xval, yval )
int   n;
float *xval;
int   *yval;
{
float *xsav;
float xtemp;
int   ytemp;
int   i;

register int indx;

if ( n < 2 ) return;

indx = n - 1;

xsav = (float *) osmmget(n * sizeof(float));
for (i = 0; i < n; i++) xsav[i] = xval[i];

do { register int jj, ii = -1;

     for ( jj = 0; jj < indx; jj++ )
         { if ( xsav[jj] > xsav[jj+1] )
              { xtemp      = xsav[jj],
                ytemp      = yval[jj];

                xsav[jj]   = xsav[jj+1];
                yval[jj]   = yval[jj+1];
                xsav[jj+1] = xtemp;
                yval[jj+1] = ytemp;
                ii = jj;
              }
         }
     indx = ii;
   } 
while ( indx != -1 );
osmmfree( (char *) xsav );
}