/* @(#)cstvals.c	17.1.1.1 (ESO-DMD) 01/25/02 17:40:33 */
/*===========================================================================
  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 Massachusetts Ave, Cambridge, 
  MA 02139, USA.
 
  Correspondence 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
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*

.COPYRIGHT: Copyright (c) 1988 European Southern Observatory,         
                                         all rights reserved         

.LANGUAGE: C

.AUTHOR: K.Banse                                                      

.IDENTIFICATION                                                       
  routine Cstvals		version 1.00	880912                  
  K. Banse			ESO _ Garching                          
  1.10	   890426	1.20	900704	1.30	910228          
  2.00	   920401	2.10	920703	2.15	931208

.KEYWORDS                                                             
  statistics, mean, min, max                                          

.PURPOSE                                                              
  calculate min, max and other statistical values                     

.ALGORITHM                                                            
  usual mathematics                                                   

.INPUT/OUTPUT                                                         
  call as                                                             
  stat = Cstvals(action,a,naxis,npix,sublo,subhi,cutvls,results,respix,nopix)

  input par:                                                          
  action:      c_string    = MIN for min, max calculation            
                           = MEAN for above + mean values, std deviation
                           = ALL for above + moments of 2. + 3. order 
                           = SPEC for mean absdev of mean, median
                           = XMEAN, XALL, XSPEC for repeated calls
                           = ZMEAN, ZALL, ZSPEC for last repeated call

                           for X/ZMEAN (X/ZALL) results[3...] and nopix also
			   serve as input pars!!
 
  naxis:       I*4         no. of axis of A (maxium of 3, currently)  
  npix:        I*4 array   no. of pixels per axis in A                
  sublo:       I*4 array   start pixels of subframe                   
  subhi:       I*4 array   end pixels of subframe                     
  cutvls:      R*4 array   user supplied cutvalues for subframe       

  in/output par:                                                         
  results:      R*4 array  depending on ACTION, this array contains   
                           min, max, mean, sigma, 2-moment, 3-moment, 
                           total intensity                           
  respix:      I*4 array   holds pixel no. of min + max of array     
  nopix:       I*4         total no. of pixels inside subframe       
                           with value in (cutvls[0],cutvls[1])     
.RETURNS
  return status            = 0 o.k., otherwise error             

.VERSIONS                                                        
000904			last update

----------------------------------------------------------------------*/

#include <math.h>
#include <midas_def.h>

/*

*/

#ifdef __STDC__
int Cstvals(char *action,float *a,int naxis,int *npix,int *sublo,int *subhi,
            float *cutvls,float *results,int *respix,int *nopix)
 
#else
int Cstvals(action,a,naxis,npix,sublo,subhi,cutvls,results,respix,nopix)
char       *action;
float      *a, cutvls[2], *results;
int   naxis, *npix, *sublo, *subhi, *respix, *nopix;
#endif

{
int   hiy, hiz, lowy, lowz;
int   npx, npxy;
int   io, optio, yoff1, zoff1, mopix;
int   binc, nlo, xinc, lowx, hix, maxpix;
register int  nr, ny, nz, minpix;

float    dfcuts;
float    *b, fmin, fmax;
register float   pix, fc0, fc1;

register double   dmean, dmedian, pixd, pix2;
static double     sum=0.0, sum2=0.0, sum3=0.0, sum4=0.0;
static double     discr=0.0, stadev=0.0, dsize=0.0;
static double     tsum=0.0, tsum2=0.0, tsum3=0.0, tsum4=0.0;




if (naxis >  3) return(1);                  /* max. 3 dimensions supported */


/* set up window in father frame                                        */

lowx = sublo[0];
hix = subhi[0];
xinc = hix - lowx;
npx = npix[0];
binc = lowx + npx - hix - 1;	      /* pointer inc between adjacent lines */

optio = 0;
if (action[0] >= 'X') 
   io = 1;
else
   {
   io = 0;
   *nopix = 0;
   }

if (action[io] == 'A') 
   optio = 1;
else if (action[io] == 'S')
   optio = 2;


mopix = hix - lowx + 1;
if (naxis >= 2)     
   {
   lowy = sublo[1];
   hiy = subhi[1];
   npxy = npx * npix[1];
   mopix *= (hiy - lowy + 1);
   }
else
   {
   lowy = 0;
   hiy = 0;
   npxy = npx;
   }

if (naxis == 3)     
   {
   lowz = sublo[2];
   hiz = subhi[2];
   mopix *= (hiz - lowz + 1);
   }
else
   {
   lowz = 0;
   hiz = 0;
   }
if (mopix <= 0) return(3);			/* check correct limits */


/*  init offsets + min, max */

zoff1  = lowz * npxy;
yoff1= lowy * npx;


dfcuts = cutvls[1] - cutvls[0];
if (dfcuts <= 0.)     
   {
   minpix = zoff1 + yoff1 + lowx;
   maxpix = minpix;
   if (tsum2 > 0.0)				/* it's not the first time */
      {
      fmin = results[0];
      fmax = results[1];
      }
   else
      {
      fmin = a[minpix];
      fmax = fmin;
      }
   }
else
   {
   mopix =  0;					/* reset to zero */
   fc0 = cutvls[0];   fc1 = cutvls[1];

   for (nz=lowz; nz<=hiz; nz++)
      {
      nlo = lowx + yoff1 + zoff1;
      b = a + nlo;
      for (ny=lowy; ny<=hiy; ny++)
         {
         for (minpix=nlo; minpix<=nlo+xinc; minpix++)
            {
            fmin = *b++;

            if ((fmin >= fc0) && (fmin <= fc1))
               {        
               fmax = fmin;
               maxpix = minpix;
               if (naxis == 3)
                  {
                  lowz = nz;			/* save position */
                  zoff1  = lowz * npxy;
                  }
               else if (naxis == 2)
                  {
                  lowy = ny;
                  yoff1= lowy * npx;
                  }
               else				/* naxis = 1 */
                  {
                  lowx = minpix;
                  xinc = hix - lowx;
                  }
               goto first_pix_found;
               }
            }
         nlo += npx;
         b += binc;
         }
      zoff1 += npxy;
      }


   /*  we come here, if no pixel of this chunk is in the valid interval */

   if (action[0] != 'Z') return(2);	/* not last chunk - so go on */


   /*  last chunk, did we get anything at all?  */

   if (*nopix == 0) return (2);		/* no valid pixel in any chunk */


   dsize = (double) *nopix;
   fmin = results[0];                   /* use min, max found before */
   fmax = results[1];

   if (action[io] == 'A') 
      goto sect_9800;			/* ALL */
   else if (action[io+1] == 'E')
      goto sect_9800;			/* MEAN */
   else if (action[io] == 'S')		/* SPEC */
      goto sect_9850;
   else
      goto sect_9900;			/* MINMAX */


  first_pix_found:
   if (*nopix > 0)				/* if not the first time */
      {						/* compare with previous */
      if (fmin > results[0]) fmin = results[0];	/* min, max values */
      if (fmax < results[1]) fmax = results[1];
      }
   }

if (optio == 1)
   goto sect_6000;				/* ALL */
else if (optio  == 2)
   goto sect_5000;				/* SPEC */

if (action[io+1] == 'E') goto sect_4000;		/* MEAN */


/* ******                                                               */
/*  here for MINMAX calculation                                         */
/* ******                                                               */

if (dfcuts > 0.) 
   {
   for (nz=lowz; nz<=hiz; nz++)
      {
      nlo = lowx + yoff1 + zoff1;
      b = a + nlo;
      for (ny=lowy; ny<=hiy; ny++)
         {
         for (nr=nlo; nr<=nlo+xinc; nr++)
            {
            pix = *b++;
   
            if ((pix >= fc0) && (pix <= fc1)) 
               {		/*  test valid pixel against old min + max  */
               if (pix < fmin)     
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax) 
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               mopix ++;
               }
            }
         nlo += npx;
         b += binc;
         }
      zoff1 += npxy;
      }
   }

else
   {				/*  here without user cuts  */
   if (binc > 0)
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pix = *b++;

               if (pix <  fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax)
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               }
            nlo += npx;
            b += binc;
            }
         zoff1 += npxy;
         }
      }

   else
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pix = *b++;

               if (pix <  fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax)
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               }
            nlo += npx;
            }
         zoff1 += npxy;
         }
      }
   }


*nopix += mopix;
goto sect_9900;


/* ******                                                               */
/*  here for MINMAX + MEAN calculation                                  */
/* ******                                                               */

sect_4000:


/*  we use cut values...                                                */

if (dfcuts > 0.) 
   {
   for (nz=lowz; nz<=hiz; nz++)
      {
      nlo = lowx + yoff1 + zoff1;
      b = a + nlo;
      for (ny=lowy; ny<=hiy; ny++)
         {
         for (nr=nlo; nr<=nlo+xinc; nr++)
            {
            pix = *b++;

            if ((pix >= fc0) && (pix <= fc1))
               {                  /*  test valid pixel against old min + max  */
               if (pix < fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax) 
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               pixd = (double) pix;
               sum += pixd;				/*sum of pixels*/
               sum2 += (pixd*pixd);			/*sum of squares*/
               mopix ++;
               }
            }
         nlo += npx;
         b += binc;
         }
      zoff1 += npxy;
      }
   }
else
   { 				/*  here without user cuts */
   if (binc > 0)
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pix = *b++;

               if (pix <  fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax)
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               pixd = (double) pix;
               sum += pixd;				/*sum of pixels*/
               sum2 += (pixd*pixd);			/*sum of squares*/
               }
            nlo += npx;
            b += binc;
            }
         zoff1 += npxy;
         }
      }

   else
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pix = *b++;

               if (pix <  fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax)
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               pixd = (double) pix;
               sum += pixd;				/*sum of pixels*/
               sum2 += (pixd*pixd);			/*sum of squares*/
               }
            nlo += npx;
            }
         zoff1 += npxy;
         }
      }
   }


*nopix += mopix;
dsize = (double) *nopix;
goto sect_9800;


/* ******                                                               */
/*  here for SPEC                                                       */
/* ******                                                               */

sect_5000:
dmean = results[0];
dmedian = results[1];

if (dfcuts > 0.)                        /*  we use cut values...  */
   {
   for (nz=lowz; nz<=hiz; nz++)
      {
      nlo = lowx + yoff1 + zoff1;
      b = a + nlo;
      for (ny=lowy; ny<=hiy; ny++)
         {
         for (nr=nlo; nr<=nlo+xinc; nr++)
            {
            pix = *b++;

            if ((pix >= fc0) && (pix <= fc1))
               {                  /*  test valid pixel against old min + max  */
               pixd = (double) pix;                      /* move to double */
               sum += fabs (pixd - dmean);	/* sum of pixels - mean */
               sum2 += fabs (pixd - dmedian);	/* sum of pixels - median */
               mopix ++;
               }
            }
         nlo += npx;
         b += binc;
         }
      zoff1 += npxy;
      }
   }
else
   {                                    /*  here without user cuts */
   if (binc > 0)
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pixd = (double) *b++;		/* move to double */
               sum += fabs (pixd - dmean);	/* sum of pixels - mean */
               sum2 += fabs (pixd - dmedian);	/* sum of pixels - median */
               }
            nlo += npx;
            b += binc;
            }
         zoff1 += npxy;
         }
      }

   else
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pixd = (double) *b++;		/* move to double */
               sum += fabs (pixd - dmean);	/* sum of pixels - mean */
               sum2 += fabs (pixd - dmedian);	/* sum of pixels - median */
               }
            nlo += npx;
            }
         zoff1 += npxy;
         }
      }
   }

*nopix += mopix;
dsize = (double) *nopix;
goto sect_9850;



/* ******                                                               */
/*  here for ALL                                                        */
/* ******                                                               */

sect_6000:
if (dfcuts > 0.)			/*  we use cut values...  */
   {
   for (nz=lowz; nz<=hiz; nz++)
      {
      nlo = lowx + yoff1 + zoff1;
      b = a + nlo;
      for (ny=lowy; ny<=hiy; ny++)
         {
         for (nr=nlo; nr<=nlo+xinc; nr++)
            {
            pix = *b++;
   
            if ((pix >= fc0) && (pix <= fc1))
               {                  /*  test valid pixel against old min + max  */
               if (pix < fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax) 
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               pixd = (double) pix;                      /* move to double */
               sum += pixd;                           /*sum of pixels*/
               pix2 = pixd * pixd;
               sum2 += pix2;                          /*sum of squares*/
               pix2 *= pixd;
               sum3 += pix2;                          /*sum of 3rd power*/
               sum4 += (pix2*pixd);                   /*sum of 4th power*/
               mopix ++;
               }
            }
         nlo += npx;
         b += binc;
         }
      zoff1 += npxy;
      }
   }
else
   { 					/*  here without user cuts */
   if (binc > 0)
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pix = *b++;

               if (pix <  fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax) 
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               pixd = (double) pix;                      /* move to double */
               sum += pixd;                           /*sum of pixels*/
               pix2 = pixd * pixd;
               sum2 += pix2;                          /*sum of squares*/
               pix2 *= pixd;
               sum3 += pix2;                          /*sum of 3rd power*/
               sum4 += (pix2*pixd);                   /*sum of 4th power*/
               }
            nlo += npx;
            b += binc;
            }
         zoff1 += npxy;
         }
      }

   else
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         nlo = lowx + yoff1 + zoff1;
         b = a + nlo;
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=nlo; nr<=nlo+xinc; nr++)
               {
               pix = *b++;

               if (pix <  fmin)
                  {
                  fmin = pix;
                  minpix = nr;
                  }
               else if (pix >  fmax)
                  {
                  fmax = pix;
                  maxpix = nr;
                  }
               pixd = (double) pix;	            /* move to double */
               sum += pixd;                     /*sum of pixels*/
               pix2 = pixd * pixd;
               sum2 += pix2;                    /*sum of squares*/
               pix2 *= pixd;
               sum3 += pix2;                    /*sum of 3rd power*/
               sum4 += (pix2*pixd);             /*sum of 4th power*/
               }
            nlo += npx;
            }
         zoff1 += npxy;
         }
      }
   }

*nopix += mopix;
dsize = (double) *nopix;



sect_9800:
if ((io == 0) || (action[0] == 'Z'))
   {
   sum += tsum;
   sum2 += tsum2;
   if ((fmax-fmin) < 10.e-30) 
      {
      results[2] = fmin;
      results[3] = 0.0;
      stadev = 0.0;
      }
   else
      {
      discr = sum2  - (sum*sum)/dsize;
      if  (discr <  0.0) discr = -discr;
      if (dsize > 1.0) discr /= (dsize - 1.0); 
      dmean = sum/dsize;					/* mean */
      results[2] = (float) dmean;
      stadev = sqrt(discr);				/* std. dev */
      results[3] = (float) stadev;
      }
   }


/*  calculate 3rd + 4th moment according to (Schaums series Statistics):

3rd moment = m3 = 1/size * SUM of (x - mean)**3
4th moment = m4 = 1/size * SUM of (x - mean)**4    

    and finally the dimensionless form:

a3 = m3/s**3			(with s = standard deviation)
a4 = m4/s**4							*/



sect_9850:
if (optio == 1)
   {
   if ((io == 0) || (action[0] == 'Z'))
      {
      double   sigma, dm2, dva, dvb, dvc, dres;

      sum3 += tsum3;
      sum4 += tsum4;

      if (stadev < 1.e-30)
         {
         results[4] = results[5] = 0.0;
         }
      else
         {
         dm2 = dmean * dmean;

         dva = dmean * sum2;
         dvb = dm2 * sum;
         dvc = dsize * (dm2*dmean);
         dres = sum3 - 3.0*(dva - dvb) - dvc;
         sigma = stadev * stadev * stadev * dsize;	/* size * std**3 */
         results[4] = (float) (dres / sigma);

         dva = sum3 * dmean;
         dvb = sum * dmean * dm2;
         dvc = 6.0 * sum2 * dm2;
         dres = sum4 - 4.0*(dva + dvb) + dvc + (dsize*dm2*dm2);
         sigma = discr * discr * dsize;			/* size * std**4 */
         results[5] = (float) (dres / sigma);
         }
      results[6] = (float) sum;
      }

   }
else if (optio == 2)
   {
   if ((io == 0) || (action[0] == 'Z'))
      {
      results[0] = (float) (sum / dsize);
      results[1] = (float) (sum2 / dsize);
      sum = 0.0;                   /* initialize all sums for next time*/
      sum2 = 0.0;
      }
   return (0);
   }


sect_9900:
results[0] = fmin;
results[1] = fmax;
respix[0] = minpix;
respix[1] = maxpix;

if ((io == 0) || (action[0] == 'Z'))
   {
   sum = sum2 = sum3 = sum4 = 0.0;
   tsum = tsum2 = tsum3 = tsum4 = 0.0;
   }
else
   {
   tsum += sum;
   tsum2 += sum2;
   tsum3 += sum3;
   tsum4 += sum4;
   sum = sum2 = sum3 = sum4 = 0.0;
   }

return(0);
}

/*

*/

void statfunc(hstart,bsize,nobins,histo,adbins,results)
int   nobins, adbins;
float      hstart, bsize, *results;
int  *histo;

{
int    n, total;
register int  npos, mr, nr;

float f, oldf, half, x0, vr;


/*  do not use the excess bins!  */

if (adbins > 0)
   {
   histo[0] = 0;              /*also note, that hstart points to histo[0]...*/
   histo[nobins-1] = 0;
   }


total = histo[0];
for (nr=1; nr<nobins; nr++) total += histo[nr];


/*  check for 0, 1 or 2 valued histograms */

if (total <= 0)
   {
   SCTPUT("histogram null except excess bins - modes, median set to 0.0 ...");
   *results++ = 0.0;
   *results++ = 0.0;
   *results = 0.0;
   return;
   }


/*  calculate smallest mode:
    xmod such that p(xmod) = maximum (first, if several)
    as well as real mode, i.e. max of hist bins          */

oldf = histo[0];
npos = 0;
for (nr=1; nr<nobins; nr++)
   {
   f = histo[nr];
   if (f > oldf)
      {
      oldf = f;
      npos = nr;
      }
   else if (f < oldf)				/* 1. maximum found */
      break;
   }
*results++ = hstart + (npos+0.5)*bsize;		/* take midpoint of bin */

for (mr=nr+1; mr<nobins; mr++)
   {
   f = histo[mr];
   if (f > oldf)				/* look for real maximumn */
      {
      oldf = f;
      npos = mr;
      }
   }
*results++ = hstart + (npos+0.5)*bsize;		/* take midpoint of bin */


/*  calculate median:
    xmed such that F(xmed) = 0.5, F the cumulative distribution function*/

half = total/2;				/* calculate no. of points/2 */
f = 0.0;				/* cumulative distribution function */
oldf = 0.0;

for (npos=0; npos<nobins; npos++)
   {
   f += histo[npos];
   if (f < half)
      oldf = f;                       /*save last function value */
   else
      break;
   }


/*  find xmed according to: ymed = y0 + ((y1-y0)/binsize) * (xmed-x0)  */

if (npos <= 1)
   x0 = hstart;
else
   x0 = hstart + (npos-1)*bsize;                    /*get low x*/
vr = f - oldf;
if (vr < 0.0 ) vr = -vr;
if (vr < 10.e-30)
   *results = bsize*(half-oldf) + x0;
else
   *results = bsize*((half-oldf)/(f-oldf)) + x0;

}