/* @(#)zstats.c	17.1.1.1 (ESO-DMD) 01/25/02 17:40:37 */
/*===========================================================================
  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
===========================================================================*/


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

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

#define MAXBIN 1024
#define NINT(a)      ((a) < 0 ? (int)((a) - 0.5 ) : (int)((a) + 0.5) )

/*

*/

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

.IDENTIFICATION
  routine Chistvals             version 1.00    920403
  K. Banse                      ESO - Garching

.KEYWORDS
  histogram

.PURPOSE
  calculate histogram

.ALGORITHM
  the NSLOT bins are filled in the following way:

a) no excess bins:
  all pixels in the half-open interval [a,b[
  with a = MIN+n*binsize and b = MIN+(n+1)*binsize (n=0, .., NSLOT-2)
  contribute to histogram slot #n (in total NSLOT-1 slots)
  all pixels equal to MAX contribute to slot #(NSLOT-1) to get finally
  NSLOT slots filled

  (MIN = minimum, MAX = maximum of data)

b) excess bins:
  all pixels in the half-open interval [-inf,MIN[ goto the low excess bin
  which is slot #0
  all pixels in the half-open interval [a,b[
  with a = MIN+n*binsize and b = MIN+(n+1)*binsize (n=0, .., NSLOT-3)
  contribute to histogram slot #(n+1) (in total NSLOT-2 slots)
  all pixels in the half-open interval [MAX,+inf] goto the high excess bin
  which is slot #(NSLOT-1)

  (MIN = low excess value, MAX = high excess value)


.VERSIONS
  1.00            from Fortran version of HSTVLS
 980622
-----------------------------------------------------------------------*/

#ifdef __STDC__
void  Chistvals(float *array, int naxis, int *npix, int *sublo, 
                int *subhi, float *cuts, float sltsiz, int nslot, int *slot)
#else
void Chistvals(array,naxis,npix,sublo,subhi,cuts,sltsiz,nslot,slot)

float  *array;		/* IN: data buffer */
int    naxis;		/* IN: no. of axes of data array */
int    *npix;		/* IN: NPIX of array */
int    *sublo;		/* IN: start pixels in each axis */
int    *subhi;		/* IN: end pixels in each axis */
float  *cuts;		/* IN: low + high cuts of data */
float  sltsiz;		/* IN: binsize */
int    nslot;		/* IN: no. of bins (without excess bins...) */
int    *slot;		/* OUT: vector of bins for histogram */
#endif

{
int   lowx, lowy, lowz, hix, hiy, hiz;
int   yoffa, zoff, npx, npxy;
int   binc, x, nslotm;
register int  nr, ny, nz;

float  *b;
register float  freg, f, fmax, fmin;

double df;



if (nslot < 1)
   {
   *slot = 0;
   return;
   }


/*  initialize  +  determine subarea */

lowx = sublo[0];
hix = subhi[0];
npx = npix[0];
binc = lowx + npx - hix - 1;		/* increment between adjacent lines */
nslotm = nslot - 1;			/* constant we need later on */

if (naxis >= 2)     
   {
   lowy = sublo[1];
   hiy = subhi[1];
   npxy = npx * npix[1];
   }
else
   {
   lowy = 0;
   hiy = 0;
   npxy = npx;
   }
if (naxis >= 3)     
   {
   lowz = sublo[2];
   hiz = subhi[2];
   }
else
   {
   lowz = 0;
   hiz = 0;
   }

zoff = lowz * npxy;
yoffa = (lowy * npx) + lowx;
fmax = cuts[1];
fmin = cuts[0];
df = 1.0 / (double) sltsiz;


/*  test, if we have excess bins */

if (fmax > fmin)
   {		/*  main loop over all pixels in given area with excess bins */
   for (nz=lowz; nz<=hiz; nz++)
      {
      b = array + yoffa + zoff;		 /* reset y offset for each z-loop */
      for (ny=lowy; ny<=hiy; ny++)
         {
         for (nr=lowx; nr<=hix; nr++)
            {
            freg = *b++;
            if (freg >= fmax)
               x = nslotm;		/* nslot-1 : high excess bin */
            else
               {
               if (freg < fmin)
                  x = 0; 			/* low excess bin */
               else
                  {
                  x = (int) floor(df*((double)(freg-fmin)));
                  x ++;				/* because of excess bin */
                  }
               }
            slot[x] ++;
            }
         b += binc;
         }
      zoff += npxy;
      }
   }
else
   {	    /* main loop over all pixels in given area without excess bins */
   if (binc > 0)
      {
      for (nz=lowz; nz<=hiz; nz++)
         {
         b = array + yoffa + zoff;	 /* reset y offset for each z-loop */
         for (ny=lowy; ny<=hiy; ny++)
            {
            for (nr=lowx; nr<=hix; nr++)
               {
               freg = *b++;
               x = (int) floor(df*((double)(freg-fmin)));

               if (x == nslot)
                  slot[nslotm] ++;		/* avoid high out bin */
               else
                  slot[x] ++;
               }
            b += binc;
            }
         zoff += npxy;
         }
      }
   else
      {
      register double  dr;

      if (fabs(fmin) < 10.e-33)
         {
         for (nz=lowz; nz<=hiz; nz++)
            {
            b = array + yoffa + zoff;       /* reset y offset for each z-loop */
            for (ny=lowy; ny<=hiy; ny++)
               {
               for (nr=lowx; nr<=hix; nr++)
                  {
                  dr = (double) *b++;
                  x = (int) floor(df*dr);

                  if (x == nslot) 
                     slot[nslotm] ++;		/* avoid high out bin */
                  else
                     slot[x] ++;
                  }
               }
            zoff += npxy;
            }
         }
      else
         {
         register double dmin;

         dmin = (double) fmin;
         for (nz=lowz; nz<=hiz; nz++)
            {
            b = array + yoffa + zoff;       /* reset y offset for each z-loop */
            for (ny=lowy; ny<=hiy; ny++)
               {
               for (nr=lowx; nr<=hix; nr++)
                  {
                  dr = (double) *b++;
                  dr -= dmin;
                  x = (int) floor(df*dr);

                  if (x == nslot) 
                     slot[nslotm] ++;		/* avoid high out bin */
                  else
                     slot[x] ++;
                  }
               }
            zoff += npxy;
            }
         }
      }
   }
}

/*

*/

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

.IDENTIFICATION
  routine Zstats
  K. Banse                      ESO - Garching

.KEYWORDS
  maximum, minimum,mean value, std. dev., histogram

.PURPOSE
  calculate statistics of given bulk data frame (opened for real data)

.ALGORITHM
  use routine Cstvals to calculate moments of up to order 4
  and standard deviation as well as min + max of given frame
  use HISTO to calculate the histogram

.INPUT/OUTPUT
  call as Zstats(imno,area,naxpix,zbins,formstr,defaul)

  input par:
  imno:         int     no. of input frame
  area:         char. string    subframe area [...,...]
                                if it is R[...] or W[...] we work
                                in row/column mode
  naxpix:       int array       number of axes, no. of pixels per axis
  zbins:        float array     holds nobins/binsize + limits for excess bins
                                (lowcut + hicut)
  formstr:      char. string    format string for min, max, mean, std_dev
  defaul:       char. string    [0] = N, subframe is given in AREA
                                    = M, subframes are merged into artificial
                                         frame in memory (no descrs)
                                    = Y, use complete frame
                                [1] = N, binsize is given in zbins[0]
                                    = Y, total no. of bins is in zbins[0]
                                [2] = F(ull), do complete statistics;
                                    = G, as F with exact median calculation;
                                    = V, as F + mean deviation of mean + median
                                    = W, as G + mean deviation of mean + median
                                    = R(educed), omit histogram stuff,
                                       also no median + smallest mode, since
                                       histogram needed for their calculation;
                                    = X, as R and also exact median;
                                    = S(hort), only min, max, mean + stdev;
                                    = M(inmax) only min + max;
                                    = H(istogram) only min + max + histogram;
                                [3] = F(ull) display, S(hort) display,
                                      X(short) very short display, N(o) display
                                [4] = DSCYES flag,
                                    = N, write no descriptors
                                    = Y, write descriptors
                                    = H, write only descriptors HIST_BINS
                                         and HISTOGRAM


  also the following keywords are written:
  OUTPUTR/R/1/11        will receive: Min,Max,Mean,Std,3rdMom,4thMom,Intensity
                        and Median,SmallMode,Nobins,Binsize
  OUTPUTI/I/1/7 (1) = no. valid of pixels in window
                        (2,3,4) = window start pixels, set to 0 if not appl.
                        (5,6,7) = window end pixels, set to 0 if not appl.

  and the following descriptors are written:
  STATISTICS/R/1/11             as OUTPUTR(1,...,11)
  HISTOGRAM/I/1/NOBINS          NOBINS = no. of bins
  HIST_BINS/R/1/4               no. of bins, binsize + excess bins
  WIN_DOW_FROM/I/1/NAXIS        start pixels of window in frame
  WINDOW_TO/I/1/NAXIS           end pixels of window
  LHCUTS/R/3/2                  min, max in case of full frame only

.VERSIONS
[1.00]  870626  from version 3.50 of STATIST as of 870612

 010412		last modif

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

#ifdef __STDC__
int Zstats(int imno, char *area, int *naxpix, float *zbins,
           char *formstr, char *defaul)
#else
int Zstats(imno,area,naxpix,zbins,formstr,defaul)
char    *area, *formstr, *defaul;
int  	imno;
int     *naxpix;
float   *zbins;
#endif

{
int   naxis, npix[3], bins[MAXBIN], maxbin, outpi[7];
int   ln1, ec, el, ed;
int   iav, nulo, wno;
int   lsublo[3], lsubhi[3];

int   nolines, mylines, mynpix[3], mylo[3], myhi[3];
int   noplanes, myplanes;
int   adbins, ival, size;
int   n, nobins, norang, npx;
int   stat, zstat, subdim, hacksize, nwork;
int   kio[4], loopy, savoff, myoff, pixoff, linecnt;
int   uni, xsublo[3], xsubhi[3];
int   retbuf[3], sublo[3], subhi[3], respix[2];

unsigned int allbytes;

register int  nr;

char   *pntr;
char   input[80], output[80], fulfram, dscyes, displevl, statlevl, hack, wmode;
char   mf, myform[12];

static char mess1[] = "max. no.of bins = MAXBIN, binsize modified...";
static char err1[] = "dynamic range of data = 0.0 ... ";

float rio[2], cuts[4], range, perc;
float binsiz, hstinf[5], outpr[14];
float cutvls[2], absbin, vr, vp;
float   *a;

void hack_up();


 
naxis = *naxpix;
if (naxis >3) naxis = 3;
for (nr=0; nr<3; nr++)				/* init variables */
   {
   sublo[nr] = 0;
   subhi[nr] = 0;
   if (nr >= naxis)
      mynpix[nr] = 1;
   else
      mynpix[nr] = naxpix[nr+1];
   npix[nr] = mynpix[nr];
   }
for (nr=0; nr<14; nr++)	outpr[nr] = 0.0;	

maxbin = MAXBIN;
fulfram = defaul[0];
statlevl = defaul[2];
displevl = defaul[3];
dscyes = defaul[4];
size = 1;
ival = 0;				/* because it's an I/O par.   */
adbins = 0;			/* init to no excess bins in histogram */

cutvls[0] = zbins[1];		/*  set excess limits */
cutvls[1] = zbins[2];


/* work on format string */

(void)strcpy(myform,"15.6e");			/* default format */
if ((*formstr != ' ') && (*formstr != '\0'))
   {
   mf = CGN_LOWER(*formstr);
   if ((mf >= 'a') && (mf <= 'z'))
      {					/* Fortran style format specs */
      n = CGN_COPY(myform,&formstr[1]);
      myform[n++] = mf;
      myform[n] = '\0';
      }
   else
      (void)strcpy(myform,formstr);
   }


/*  get boundaries of subframe  */

if (fulfram == 'N') 
   {
   (void)strcpy(input,area);
   if (input[0] != '[')     
      {
      n = 1;
      wmode = input[0];
      }
   else
      n = 0;
   stat = Convcoo(1,imno,&input[n],3,&subdim,sublo,subhi);
   if (stat != 0)
      SCETER(1,"invalid coordinate input...");

   if (naxis > 2)
      npx = sublo[2] * mynpix[0] * mynpix[1];
   else
      npx = 0;
   pixoff = npx + (sublo[1]*mynpix[0]) + 1;	/* 1.pixel of starting line */
   for (nr=0; nr<subdim; nr++)
      size *= (subhi[nr]-sublo[nr]+1);
   }

else
   {
   pixoff = 1;
   subdim = naxis;			/* use complete frame */
   for (nr=0; nr<subdim; nr++)
      {
      subhi[nr] = mynpix[nr] - 1;	/* sublo already initialized to 0 */
      size *= mynpix[nr];
      }
   }

myhi[0] = subhi[0];
mylo[0] = sublo[0];			/* save sublo[0], subhi[0] */
myhi[1] = subhi[1];
mylo[1] = sublo[1];			/* save sublo[1], subhi[1] */
subhi[1] -= sublo[1];			/* and shift towards 0  */
sublo[1] = 0;
nolines = subhi[1] + 1;			/* no. of lines we really need */

myhi[2] = subhi[2];   
mylo[2] = sublo[2];
subhi[2] = 0;   
sublo[2] = 0;	
noplanes = myhi[2] - mylo[2] + 1;	/* no. of planes we really need */

savoff = pixoff;
myoff = pixoff;				/* save original offset */
myplanes = 1;
  

/*  check, if we have to work with pieces of the frame */

hack = 'N';
hack_up(npix,D_R4_FORMAT,retbuf);	/* see, if we hack it up */
hacksize = retbuf[0];
if (noplanes > 1) 
   {
   mylines = (noplanes*npix[1]) - sublo[1] - (npix[1]-subhi[1]-1);
   if (hacksize < mylines)
      {
      hack = 'X';
      if (hacksize > nolines) hacksize = nolines;
      }
   else
      {
      sublo[2] = 0;			/* start with first plane */
      subhi[2] = myhi[2] - mylo[2];
      sublo[1] = mylo[1];
      subhi[1] = myhi[1];
      }
   }
else if (nolines > hacksize)
   hack = 'X';

if (hack == 'X')
   {
   zstat = -1;
   mylines = hacksize;
   linecnt = 0;
   loopy = 1;
   subhi[1] = mylines - 1;		/* update subhi of y */
   nwork = mylines * npix[0];		/* actual space to work with */
   mynpix[1] = mylines;
   n = 1;
   }
else
   {
   nwork = noplanes * nolines * npix[0];	/* actual space to work with */
   n = 0;
   }

for (nr=0; nr<3; nr++)
   {					/* save sublo, subhi */
   xsublo[nr] = sublo[nr];
   xsubhi[nr] = subhi[nr];
   }


allbytes = nwork * sizeof(float);
pntr = malloc(allbytes);
if (pntr == (char *) 0) SCETER(33,"Could not allocate memory...");
(void) SCFGET(imno,pixoff,nwork,&iav,pntr);
a = (float *) pntr;

if (statlevl == 'M')
   (void)strcpy(&output[n],"MIN");		/* MINMAX option */
else if ((statlevl == 'H') || (statlevl == 'L'))
   (void)strcpy(&output[n],"MIN");		/* MINMAX option */
else if (statlevl == 'S')
   (void)strcpy(&output[n],"MEAN");		/* MEAN option */
else
   (void)strcpy(&output[n],"ALL");		/* ALL option */


for (nr=0; nr<7; nr++)
   outpi[nr] = 0;
for (nr=0; nr<11; nr++)
   outpr[nr] = 0.0;


/*  now look for min + max + moments + std. dev. of file                */

if (hack == 'N')                     /* no hacking */
   {
   stat = Cstvals(output,a,subdim,mynpix,sublo,subhi,cutvls,outpr,
                  respix,&ival);
   if (stat != 0)
      {
      if (stat == 2) 
         {
         SCTPUT("No pixels found with data in given interval...");
         goto no_data;
         }
      else
         SCETER(2,"problems in calculation of statistics...");
      }
   respix[0] += pixoff;		 /* [0,npix-1] => [1,npix]  */
   respix[1] += pixoff;		 /* because pixoff begins with 1  */
   }

else
   {
hack1_loop:
   subhi[1] = mynpix[1] - 1;         /* update subhi of y */
   output[0] = hack;
   stat = Cstvals(output,a,subdim,mynpix,sublo,subhi,cutvls,outpr,
                  respix,&ival);
   if (stat == 0)
      {
      zstat = 0;			/* set overall status to `valid' */
      respix[0] ++;
      respix[1] ++;

      if (loopy == 1)			/* first time here */
         {
         rio[0] = outpr[0];		/* save min + max */
         rio[1] = outpr[1];
         kio[0] = respix[0];		/* and min, max pixels */
         kio[1] = pixoff-1;		/* with current pixel offset */
         kio[2] = respix[1];
         kio[3] = pixoff-1;
         loopy = 2;
         }

      else	/* while looping we update min + max with their pixels */
         {
         if (rio[0] > outpr[0]) 		/* new min smaller ?  */
            {
            rio[0] = outpr[0];
            kio[0] = respix[0];		 /* [0,npix-1] => [1,npix]  */
            kio[1] = pixoff-1;	
            }
         if (rio[1] < outpr[1]) 		/* new max larger ?  */
            {
            rio[1] = outpr[1];
            kio[2] = respix[1];		/* [0,npix-1] => [1,npix]  */
            kio[3] = pixoff-1;	
            }
         }
     }

   if (hack == 'Z')			/* was it last loop ?  */
      {
      if (zstat != 0)
         {
         if (zstat == 2)
            {
            SCTPUT("No pixels found with data in given interval...");
            goto no_data;
            }
         else
            SCETER(2,"problems in calculation of statistics...");
         }

      respix[0] = kio[0] + kio[1];		/* add offset */
      respix[1] = kio[2] + kio[3];
      outpr[0] = rio[0];		/* move back to `outpr' array */
      outpr[1] = rio[1];		/* move back to `outpr' array */
      }

   else
      {
      linecnt += mynpix[1];
      mynpix[1] = nolines - linecnt;

      if (mynpix[1] > 0)		/* last chunk in a plane */
         {
         if (mynpix[1] > hacksize)
            mynpix[1] = hacksize;
         else				/* will be last piece of this plane */
            {
            myplanes ++;
            if (myplanes > noplanes) hack = 'Z';
            }
         pixoff = (linecnt*npix[0]) + myoff;
         }

      else				/* move to next plane */
         {
         if (hacksize == nolines)	/* only here in 3-dim case */
            {				/* will be last piece !!  */
            myplanes ++;
            if (myplanes >= noplanes) hack = 'Z';
            }
         myoff += (npix[0]*npix[1]);  
         pixoff = myoff;
         linecnt = 0;                     /* reset linecount */
         mynpix[1] = mylines;             /* and work line count */
         }

      nwork = mynpix[1] * npix[0];
      stat = SCFGET(imno,pixoff,nwork,&iav,pntr);
      goto hack1_loop;
      }
   }


if (displevl != 'N')     
   {
   if (fulfram == 'N')     		/* not full frame */
      {
      if (input[0] != '[')     
         {
         n = CGN_INDEXC(input,'@') + 1;
         if ((wmode == 'R') || (wmode == 'r'))     
            {
            (void)strcpy(output,"row ");
            (void)strncpy(&output[4],&input[n],5);
            (void)strcpy(&output[9]," of frame");
            }
         else if ((wmode == 'C') || (wmode == 'c'))     
            {
            (void)strcpy(output,"column ");
            (void)strncpy(&output[7],&input[n],5);
            (void)strcpy(&output[12]," of frame");
            }
         else 			/* must be PLANE */
            {
            (void)strcpy(output,"plane ");
            (void)strncpy(&output[6],&input[n],5);
            (void)strcpy(&output[11]," of frame");
            }
         }
      else
         {
         n = CGN_INDEXC(input,' ');
         if (n > 0) 
            input[n] = '\0';
         (void)strcpy(output,"area ");
         (void)strcpy(&output[5],input);
         (void)strcat(output," of frame ");
         }
      }
   else
      {
      if (fulfram == 'Y')
         (void)strcpy(output,"complete area of frame ");
      else
         (void)strcpy(output,"merged subframes of frame ");
      }

   SCTPUT(output);			/* display area */

   (void)sprintf(input,"minimum, maximum:\t\t%%%s\t%%%s",myform,myform);
   (void)sprintf(output,input,outpr[0],outpr[1]);
   SCTPUT(output);
   }


/*  get the pixel no.s of min + max */

if (subdim == 1)     
   {
   outpi[1] = respix[0];
   outpi[2] = respix[1];
   (void)sprintf(output,"at pixel (%d),(%d)",outpi[1],outpi[2]);
   }

else if (subdim == 2)     
   {
   n = (respix[0]-1)/npix[0];
   outpi[2] = n + 1;				/*y-index*/
   outpi[1] = respix[0] - (n*npix[0]);		/*x-index*/

   n = (respix[1]-1)/npix[0];
   outpi[4] = n + 1;				/*y-index*/
   outpi[3] = respix[1] - (n*npix[0]);		/*x-index*/

   (void)sprintf(output,"at pixel (%d,%d),(%d,%d)",
                  outpi[1],outpi[2],outpi[3],outpi[4]);
   }

else
   {
   npx = npix[0] * npix[1];
   n = (respix[0]-1)/npx;
   outpi[3] = n + 1;				/*z-index*/
   respix[0] -= (n*npx);
   n = (respix[0]-1)/npix[0];
   outpi[2] = n + 1;				/*y-index*/
   outpi[1] = respix[0] - (n*npix[0]);		/*x-index*/

   n = (respix[1]-1)/npx;
   outpi[6] = n + 1;                            /*z-index*/
   respix[1] -= (n*npx);
   n = (respix[1]-1)/npix[0];
   outpi[5] = n + 1;                            /*y-index*/
   outpi[4] = respix[1] - (n*npix[0]);          /*x-index*/

   (void)sprintf(output,"at pixel (%d,%d,%d),(%d,%d,%d)",
           outpi[1],outpi[2],outpi[3],outpi[4],outpi[5],outpi[6]);
   }
if (displevl == 'F') SCTPUT(output);	/* if desired, display it */

norang = 1;			/* set to NO histogram */

if (statlevl == 'M')
   goto sect_4000;
else if (statlevl == 'H')
   goto hist_prep;


/* handle mean + standard deviation */

if (displevl != 'N')
   {
   (void)sprintf(input,"mean, standard_deviation:\t%%%s\t%%%s",myform,myform);
   (void)sprintf(output,input,outpr[2],outpr[3]);
   SCTPUT(output);
   }

if (statlevl == 'S') goto sect_4000;


/* handle 3. + 4. moment, total intensity  */

if ((displevl == 'F') || (displevl == 'S'))
   {
   (void)sprintf(output,"3rd + 4th moment: \t\t%15g%15g",outpr[4],outpr[5]);
   SCTPUT(output);
   (void)sprintf(output,"total intensity: \t\t%15g",outpr[6]);
   SCTPUT(output);
   }


/* now prepare the histogram calculation */

if ((statlevl != 'F') && (statlevl != 'G') &&
    (statlevl != 'V') && (statlevl != 'W')) 
   goto sect_4000;				/* skip histogram stuff */

hist_prep:
if (cutvls[0] >= cutvls[1])     
   {
   range = outpr[1] - outpr[0];     /*take physical min + max*/
   hstinf[2] = outpr[0];
   hstinf[3] = outpr[1];
   cutvls[0] = outpr[0];           /*cutvls[1] = min is needed for histvals*/
   cutvls[1] = cutvls[0];           /*indicate, that no excess bins*/
   }
else
   {
   range = cutvls[1] - cutvls[0];     /*take user min + max*/
   hstinf[2] = cutvls[0];
   hstinf[3] = cutvls[1];
   adbins = 2;                       /*provide 2 excess bins*/
   }
hstinf[4] = adbins;

norang = 0;
nobins = 0;
binsiz = 0.0;

if (range <= 10.e-30)     
   {
   SCTPUT(err1);
   norang = -1;
   goto sect_4000;			/* and skip histogram business...*/
   }
 

/*  determine BINSIZE + no. of bins  */

if (defaul[1] == 'Y')                  /*test, if binsiz or no. of bins...*/
   {
   nobins = NINT(zbins[0]);
   if (nobins < 2) 
      SCETER(3,"problems with histogram binsize...");

   if (nobins > (maxbin-adbins))
      {
      SCTPUT(mess1);
      nobins = maxbin - adbins;
      }
   binsiz = range/(nobins-1);
   }

else
   {
   binsiz = zbins[0];
   if (binsiz >= 0.)
      absbin = binsiz;
   else
      absbin = -binsiz;

   if (absbin > 10.e-30)
      nobins = (range / binsiz) + 1;

   if (nobins < 2)
      SCETER(3,"problems with histogram binsize...");

   if (nobins > (maxbin-adbins))     
      {
      SCTPUT(mess1);
      nobins = maxbin - adbins;
      binsiz = range / (nobins-1);
      }
   }
nobins += adbins;                 /* maybe add 2 excess bins... */


/*  use last stored data to begin with the actual histogram calculation  */

for (n=0; n<nobins; n++)
   bins[n] = 0;			/* init bins  */

Chistvals(a,subdim,mynpix,sublo,subhi,cutvls,binsiz,nobins,bins);


if (hack != 'N')
   {
   hack = 'X';			/* do last piece first */
   loopy = pixoff;		/* save the last pixoff */
   myplanes = 1;
   myoff = savoff;
   pixoff = myoff;
   mynpix[1] = mylines;
   linecnt = 0;

hack2_loop:
   nwork = mynpix[1] * npix[0];
   stat = SCFGET(imno,pixoff,nwork,&iav,pntr);
   subhi[1] = mynpix[1] - 1;         /* update subhi of y */
   Chistvals(a,subdim,mynpix,sublo,subhi,cutvls,binsiz,nobins,bins);

   if (hack != 'Z')
      {
      linecnt += mynpix[1];
      mynpix[1] = nolines - linecnt;
      if (mynpix[1] > 0)
         {
         if (mynpix[1] > hacksize)
            mynpix[1] = hacksize;
         else                           /* will be last piece of this plane */
            {
            myplanes ++;
            if (myplanes > noplanes) hack = 'Z';
            }
         pixoff = linecnt*npix[0] + myoff;
         }
      else				/* move to next plane */
         {
         if (hacksize == nolines)
            {                           /* will be last piece !!  */
            myplanes ++;
            if (myplanes >= noplanes) hack = 'Z';
            }
         myoff += (npix[0]*npix[1]);  
         pixoff = myoff;
         linecnt = 0;                     /* reset linecount */
         mynpix[1] = mylines;             /* and work line count */
         }

      if (loopy != pixoff) goto hack2_loop;
      }
   }


/*  fill the descriptors                                                */

hstinf[0] = nobins;
hstinf[1] = binsiz;
outpr[9] = nobins;
outpr[10] = binsiz;


/*  write descriptor histogram + clear end, if necessary                */

if (fulfram != 'M')
   {
   kio[0] = nobins;

   SCECNT("GET",&ec,&el,&ed);
   ln1 = 1; iav = 0;
   SCECNT("PUT",&ln1,&iav,&iav);
   stat = SCDRDR(imno,"HIST_BINS",1,1,&iav,&perc,&uni,&nulo);
   SCECNT("PUT",&ec,&el,&ed);
   if (stat == 0)     
      {
      kio[0] = perc;
      if (kio[0] > maxbin) kio[0] = maxbin;
      if (kio[0] > nobins)     
         {
         for (n=nobins; n<kio[0]; n++)
            bins[n] = -1;		/* array `bins' is patched with -1 */
         if (dscyes != 'N')
            SCTPUT("Unused, last bins of histogram set to -1 ...");
         }
      else
         kio[0] = nobins;
      }
   if (dscyes != 'N')
      {
      char  help[72];
   
      (void) strcpy(help,"nobins,binsiz,lo_excs,hi_excs,no_exc_bins");
      (void) SCDHWR(imno,"HIST_BINS",hstinf,1,5,&uni,help);
      (void) strcpy(help,"no_pixls in each hist_bin");
      (void) SCDHWI(imno,"HISTOGRAM",bins,1,kio[0],&uni,help);
      }
   }


/*  update min,max if complete frame was taken                          */

sect_4000:
if (fulfram == 'Y')
   {
   SCECNT("GET",&ec,&el,&ed);
   ln1 = 1;				/* only for full frame see, */
   SCECNT("PUT",&ln1,&iav,&iav);	/* if descr LHCUTS exists */
   stat = SCDRDR(imno,"LHCUTS",1,4,&iav,cuts,&uni,&nulo);  
   SCECNT("PUT",&ec,&el,&ed);

   if (stat != 0)
      {
      cuts[0] = 0.;
      cuts[1] = 0.;
      }
   cuts[2] = outpr[0];
   cuts[3] = outpr[1];
   stat = SCDWRR(imno,"LHCUTS",cuts,1,4,&uni);
   }



/*  calculate median + smallest mode  */

if (norang != 1)                             /* only, if we got some range */
   {
   int  exmedian;
   float  resu[3];

   if (norang == 0)
      statfunc(hstinf[2],binsiz,nobins,bins,adbins,resu);
   else
      {
      resu[0] = resu[1] = resu[2] = outpr[0];		/* use minimum */
      }

   outpr[7] = resu[2];			/* Median */
   outpr[8] = resu[0];			/* 1. Mode */
   outpr[11] = resu[1];			/* real Mode */


   /* if defaul[2] == 'G' or 'W'  we calculated exact median before */

   if ((statlevl == 'G') || (statlevl == 'W'))
      {
      (void) SCKRDR("OUTPUTR",8,1,&iav,rio,&uni,&nulo);
      outpr[7] = rio[0];
      exmedian = 1;
      }
   else
      exmedian = 0;

   if ((displevl != 'N') && (displevl != 'X'))     
      {
      if (exmedian == 1)
         (void)sprintf(output,
               "exact median, 1. mode, mode: \t%15.6e%15.6e%15.6e",
               outpr[7],outpr[8],outpr[11]);
      else
         (void)sprintf(output,
               "median, 1. mode, mode: \t\t%15.6e%15.6e%15.6e",
               outpr[7],outpr[8],outpr[11]);
      SCTPUT(output);
      (void)sprintf(output,"total no. of bins, binsize: \t\t%d   %15.6e",
              nobins,binsiz);
      SCTPUT(output);
      }
   }
else if (statlevl == 'X') 
   {
   (void) SCKRDR("OUTPUTR",8,1,&iav,rio,&uni,&nulo);
   outpr[7] = rio[0];
   if ((displevl != 'N') && (displevl != 'X'))     
      {
      (void)sprintf(output,"exact median: \t\t\t%15.6e",outpr[7]);
      SCTPUT(output);
      }
   }


/*  write descriptors STATISTIC, WINDOW_FROM + WINDOW_TO                */

sublo[1] = mylo[1];   subhi[1] = myhi[1];	 /* reset sublo[1], subhi[1] */
sublo[2] = mylo[2];   subhi[2] = myhi[2]; 
for (nr=0; nr<3; nr++)
   {
   sublo[nr] ++;
   subhi[nr] ++;
   }

if (dscyes == 'Y')
   {
   char  help[72];

   (void) strcpy(help,
          "min,max,mean,std,mom3,mom4,intens,median,1.mode,...");
   (void) SCDHWR(imno,"STATISTIC",outpr,1,12,&uni,help);
   for (nr=0; nr<naxis; nr++)
      {
      lsublo[nr] = sublo[nr];
      lsubhi[nr] = subhi[nr];
      }
   (void) strcpy(help,"start frame pixels of used image area");
   (void) SCDHWI(imno,"WINDOW_FROM",lsublo,1,naxis,&uni,help);
   (void) strcpy(help,"end frame pixels of used image area");
   (void) SCDHWI(imno,"WINDOW_TO",lsubhi,1,naxis,&uni,help);
   }

if (displevl == 'F')     			/* full display */
   {
   if (adbins > 0)
      {
      vr = ival;
      vp = size;
      perc = 100.0 * (vr/vp);
      (void)sprintf(output,
      "# of pixels used = %d or %6.2f%% of all possible pixels (= %d)",
      ival,perc,size);
      SCTPUT(output);
      }
   else
      (void)sprintf(input,"# of pixels used = %d ",ival);
      
   if (naxis == 1)     
      (void)sprintf (output,"from %d to %d (in pixels)",sublo[0],subhi[0]);
   else if (naxis == 2)     
      (void)sprintf (output,"from %d,%d to %d,%d (in pixels)",
               sublo[0],sublo[1],subhi[0],subhi[1]);
   else
      {
      (void)sprintf(output,"from %d,%d,%d to %d,%d,%d (in pixels)",
              sublo[0],sublo[1],sublo[2],subhi[0],subhi[1],subhi[2]);
      }
   if (adbins == 0)
      {
      (void)strcat(input,output);
      SCTPUT(input);
      }
   else
      SCTPUT(output);
   }


/*  save displayed info also in keywords OUTPUTR + OUTPUTI              */

outpi[0] = ival;
(void) SCKWRR("OUTPUTR",outpr,1,12,&uni);
(void) SCKWRI("OUTPUTI",outpi,1,7,&uni);

if ((statlevl == 'V') || (statlevl == 'W'))
   {
   outpr[0] = outpr[2];
   outpr[1] = outpr[7];			/* mean, median already calculated */

   pixoff = savoff;
   ival = 0;

   if (hack == 'N')				/* no hacking */
      {
      for (nr=0; nr<3; nr++)
         {                                    /* get back sublo, subhi */
         sublo[nr] = xsublo[nr];
         subhi[nr] = xsubhi[nr];
         }

      (void) SCFGET(imno,pixoff,nwork,&iav,pntr);
      stat = Cstvals("SPEC",a,subdim,mynpix,sublo,subhi,cutvls,outpr,
                     respix,&ival);
      }
   else
      {
      for (nr=0; nr<3; nr++)
         {                                    /* get back sublo, subhi */
         sublo[nr] = mylo[nr];
         subhi[nr] = myhi[nr];
         }

      hack = 'X';
      myplanes = 1;
      myoff = savoff;
      mynpix[1] = mylines;
      linecnt = 0;

      (void) strcpy(&output[1],"SPEC");
      nwork = mynpix[1] * npix[0];

hack9_loop:
      (void) SCFGET(imno,pixoff,nwork,&iav,pntr);
      subhi[1] = mynpix[1] - 1;			/* update subhi of y */
      output[0] = hack;
      stat = Cstvals(output,a,subdim,mynpix,sublo,subhi,cutvls,outpr,
                     respix,&ival);
      if (hack != 'Z')				/* was it last loop ?  */
         {
         linecnt += mynpix[1];
         mynpix[1] = nolines - linecnt;

         if (mynpix[1] > 0)			/* last chunk in a plane */
            {
            if (mynpix[1] > hacksize)
               mynpix[1] = hacksize;
            else			/* will be last piece of this plane */
               {
               myplanes ++;
               if (myplanes > noplanes) hack = 'Z';
               }
            pixoff = (linecnt*npix[0]) + myoff;
            }

         else                              /* move to next plane */
            {
            if (hacksize == nolines)       /* only here in 3-dim case */
               {                           /* will be last piece !!  */
               myplanes ++;
               if (myplanes >= noplanes) hack = 'Z';
               }
            myoff += (npix[0]*npix[1]);
            pixoff = myoff;
            linecnt = 0;                     /* reset linecount */
            mynpix[1] = mylines;             /* and work line count */
            }
   
         nwork = mynpix[1] * npix[0];
         stat = SCFGET(imno,pixoff,nwork,&iav,pntr);
         goto hack9_loop;
         }
      }

   if (displevl != 'N') 
      {
      (void)sprintf(output,
            "mean absdev of mean, median: \t%15.6e%15.6e%",outpr[0],outpr[1]);
      SCTPUT(output);
      }

   (void) SCKWRR("OUTPUTR",outpr,13,2,&uni);

   if (dscyes == 'Y')
      (void) SCDWRR(imno,"STATISTIC",outpr,13,2,&uni);
   }



(void)free(pntr); 		/*  release allocated memory */
return(0);



no_data:
for (nr=0; nr<7; nr++)
   {
   outpr[nr] = -1.0;
   outpi[nr] = -1;
   }

(void) SCKWRR("OUTPUTR",outpr,1,7,&uni);
(void) SCKWRI("OUTPUTI",outpi,1,7,&uni);

(void)free(pntr); 		/*  release allocated memory */
}

/*

*/

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

.IDENTIFICATION
  function hack_up              version 1.00    920625
  K. Banse                      ESO - Garching

.KEYWORDS
  memory size, swap area

.PURPOSE
  calculate max. no. of lines which could be safely mapped for
  a singel 32-bit pixel frame

.ALGORITHM
  keyword MONITPAR(20) holds the no. of pixels/lines of a square
  real frame which may be safely mapped in one go


.INPUT/OUTPUT
  call as  hack_up(npix,pix_format,rebuf)

  input par:
  npix:         int array 	no. of pixels for 2 axes
  pix_format:   int 		data format (e.g. D_I2_FORMAT)
  retbuf:       int array 	no. of lines per chunk,
				no. of full chunks,
				remaining lines if any - else = 0

.RETURNS
  no. of lines to read in

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

void hack_up(npix,pix_format,retbuf)
int  npix[2], pix_format, retbuf[3];

{
int  strip, uni, mapsize, iav, ylines, nulo;


(void) SCKRDI("MONITPAR",20,1,&iav,&mapsize,&uni,&nulo);   /* for real data */
mapsize *= mapsize;			/* square image */

switch(pix_format)		/* adapt to data format */
   {
   case D_I1_FORMAT:
     mapsize *= RR_SIZE;
     break;

   case D_I2_FORMAT:
   case D_UI2_FORMAT:
     iav = (II_SIZE >> 2);
     mapsize *= iav;
     break;

   case D_R8_FORMAT:
     iav = DD_SIZE / RR_SIZE;
     mapsize /= iav;
     break;

   default:
   ;
   }


ylines = npix[1];

strip = mapsize / *npix;		/* no. of possible lines */
if (strip <= 0)
   {
   if (ylines <= 1)
      {
      retbuf[0] = 1;			/* take full 1-dim frame */
      retbuf[1] = 1;
      retbuf[2] = 0;
      }
   else
      SCETER(66,"value in MONITPAR(20) too small or NPIX(1) too large...");
   }

else if (strip >= ylines)
   {
   retbuf[0] = ylines;
   retbuf[1] = 1;
   retbuf[2] = 0;
   }
else
   {
   retbuf[0] = strip;
   retbuf[1] = ylines / strip;
   retbuf[2] = ylines - (strip*retbuf[1]);
   }
}