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

/* ++++++++++++++++++++++++++++++++++++++++++++++++++
 
.IDENTIFICATION: program ccdcomb
 
.KEYWORDS:       bulk data frames, average
 
.PURPOSE:        average a no. of image frames, the result will either be a 
                 frame with size equal to common area of all frames involved 
                 or size equal to union of all frame areas
 
.ALGORITHM:      extract input frames either from given parameter string or 
                 from catalogue, add up all frames + divide in the end
 
.INPUT/OUTPUT:   the following keys are used:
                 ACTION/C/1/2  (1:1) = M or N, for merging option
 		               (2:2) = I or W, for AVERAGE/IMAGE, /WEIGHTS
                 OUT_A/C/1/60  result frame
                 P3/C/1/80  list of frames to be averaged
 
.VERSIONS:       930311 RHW Created; original version from Klaus (G.) Banse
.VERSION:        980303 SW  Include mean-median method
                            same as in AVER/IMA: av_option = median,low,high
 010802		last modif

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

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

#define  MAX(x,y) (((x) < (y)) ? (y) : (x))

#define  MAXIMS  80
#define  MAXIMSONE  MAXIMS+1


main()

{

int          uni, ec, ed, el, ln0, ln1;
int          imnoc, imnox, imnos, imnol[MAXIMS];
int          felm, sizec, chunk[MAXIMS], chunkc;
int          naxis[MAXIMS], npix[MAXIMS][3], naxisc, npixc[3];
int          iseq, iav, nulo, zpix[3];

int          apix[3][2], cpix[3];
int          iaux[10];
int          ncombine[MAXIMS], nout;
int          stat, nolin, begin;
int          stripe, sig, optio, exist;
int          nrcol, nrrow, nsort, allcol, allrow;
int          tid, nrow, sel, colnum, icount;
int          frmcnt, nulcnt, ll, lll, lu, inp, m, n, ibuf;
int          weight, mmmscale, mmmoffset, expscale;
int          scale, statist, verbose;
int          inul;

int          nr, iact, snpix[4], spix[3], epix[3];
 
short int    *xpntr;

char         *wpntr, *tmpntr;
char         line[84], action[4], mexp[3];
char         expname[21], exptype[21], expfram[21];
char         frame[MAXIMS][64], framec[64], framesig[64], framecnt[64];
char         prefix[64], frname[64];
char         catfil[64], tblfil[64], colnam[68];
char         cunit[64], ident[72], output[64];
char         mod_area[40], im_area[40], otime[20];
char         defaul1[5],  null[2];
char         meth[30];

static float minundef = -99999;
static float maxundef =  99999;

static char  error1[] = "*** FATAL: operands do not match in stepsize... ";
static char  error2[] = "*** FATAL: operands do not overlap... ";
static char  error3[] = "*** FATAL: step sizes of different sign... ";
static char  error4[] = "*** FATAL: catalogue empty... ";
static char  error5[] = "*** FATAL: Unknown combining option... ";
static char  frmstr[] = "\0"; 

double       step[MAXIMS][3], start[MAXIMS][3];
double       stepc[3], ostep[3], endc[3];
double       startc[3], ostart[3], oldend[3];
double       exptime[MAXIMS];
double       aostep;
double       sstart[3], sstep[3];

 
float        Zeroes[13];
float        zbins[3];
float        *pntr[MAXIMS], *pntrc, *zpntr, *pntrs, *pntrm;
float        w[MAXIMS];
float        dif, eps[3], cutc[4], cuts[4], cutn[4];
float        exposure;
float        statics[12];
float        zeros[MAXIMS];
float        scales[MAXIMS];
float        mmm[MAXIMS];
float        wts[MAXIMS];
float        usrnul;
float        faux[5];
float        sumz, sums, sumw;
float        meanz, means;
float        low, high;
float        lowsigma, highsigma;
float        clip[2]; 
float        medx;
 
/* set up MIDAS environment + enable automatic error abort */

SCSPRO("ccdcomb");

for (n=0; n<3; n++)
   {
   apix[n][0] = 0;
   apix[n][1] = 0;
   cpix[n] = 0;
   startc[n] = 0.0;
   stepc[n] = 1.0;
   endc[n] = 0.0;
   npixc[n] = 1;
   zpix[n] = 1;
   ostart[n] = 0.0;
   oldend[n] = 0.;
   ostep[n] = 1.0;
   npix[0][n] = 1;
   start[0][n] = 0.0;
   step[0][n] = 1.0;
   spix[n] = 0;                            /* added by swolf@eso.org */
   epix[n] = 0;                            /* used for exmed() */
   }
nulcnt = 0;
icount = 0; 
zbins[0] = 256; zbins[1] = zbins[2] = 0;

/* get result frame, input frame list */
stat = SCKGETC("INCAT",1,80,&iav,line);

stat = SCKGETC("OUTFRM",1,60,&iav,framec);     /* overwritten if table input */

stat = SCKGETC("EXP",1,2,&iav,mexp);
CGN_UPSTR(mexp); 		                    /* convert to upper case */

stat = SCKGETC("EXP_DESC",1,20,&iav,expname);        /* get exp. descr. name */
CGN_UPSTR(expname);		                    /* convert to upper case */

stat = SCKGETC("EXPTYP",1,20,&iav,exptype);          /* get allowed exp type */
CGN_UPSTR(exptype);		                    /* convert to upper case */

stat = SCKGETC("ACTION",1,2,&iav,action);
CGN_UPSTR(action);

/* get the verbose option */
stat = SCKGETC("VERBOSE",1,3,&iav,output);
output[4] = '\0';
CGN_UPSTR(output);		                    /* convert to upper case */
if (strcmp(output,"YES") == 0)		                         /* set flag */
   {
   verbose = 1;
   strcpy(defaul1,"NYGXY");  /* subframe, binsize in zbins[0], exact median, */
   }                         /* short output, update descr. STATISTIC */
else
   {
   verbose = 0;
   strcpy(defaul1,"NYGNY");  /* subframe, binsize in zbins[0], exact median, */
   }                         /* NO output, update descr. STATISTIC */


stat = SCKGETC("IM_SEC",1,40,&iav,im_area);

/* ---------------------------------------------------------------------------
the auxiliary arrays iaux + faux contain the following:

iaux[0]   = 1, input data is inside result space
            0, only initialize counter pixels if necessary
iaux[1]   = used for mean-median-method
iaux[2]   = 0, take all data
          = 1, take only data in interval [ faux[2],faux[3] ]
iaux[3,4] = used to take index interval in mean-median-method (2.3.98)
iaux[5]   = 0, `nomerge' option
          = 1, `merge' option
iaux[6]   = frame count
iaux[7]   = no. of current frame
iaux[8]   = 0, use user BLANK  value
          = 1, use last pixel as BLANK  value

faux[0,1] = average limits if iaux[1] == 1 (2.3.98, swolf@eso.org)
faux[2,3] = valid_data_interval if iaux[2] = 1
faux[4]   = weighting factor
---------------------------------------------------------------------------- */

for (n=0; n<10; n++)
  {
  iaux[n] = 0;
  Zeroes[n] = 0.0;
  }
Zeroes[10] = 0.0;
Zeroes[11] = 0.0;
Zeroes[12] = 0.0;

 
/* get average option */
optio = 2;
stat = SCKGETC("METHOD",1,60,&iav,output);
CGN_UPCOPY(meth,output,10);	                      /* upper case -> frame */

if (meth[0] == 'S')
   if (meth[1] == 'I')                                     /* sigma clipping */
      {                
      optio = 7;
      strcpy(meth,"sigma clipping");
      }
   else                                                    /* Sum data points*/
      {
      optio = 1;
      strcpy(meth,"sum");
      }

else if (meth[0] == 'A')                                  /* Average options */
   {
   if (meth[2] == 'S')                             /* AVerage Sigma clipping */
      {
      optio = 8;
      strcpy(meth,"average sigma clipping");
      }
   else                                                  /* straight average */
      {
      optio = 2;                                       
      strcpy(meth,"average");
      }
   }

else if (meth[0] == 'M')                                        /* M options */
   {
   if (meth[1] == 'A')                                      /* MAx rejection */
      {
      optio = 5;
      strcpy(meth,"maximum rejection");
      }

   else if (meth[1] == 'I')
      { 
      if (meth[3] == 'M')                                /* MINMAX rejection */
	 {
         optio = 6;
         strcpy(meth,"minmax rejection");
         }
      else
         { 
         optio = 4;                                     /* MINimum rejection */
         strcpy(meth,"minimum rejection");
         } 
      }
   else if (meth[1] == 'M')                         /* MMedian (mean median) */
      {
      optio = 9;
      strcpy(meth,"mean median");
      }

   else 
      {
      optio = 3;                                                   /* MEDIAN */
      strcpy(meth,"median");
      }
   }
    
else
   SCETER(4,error5);                             /* unknown combining option */

/* get averaging options (only for optio > 1) */
if (optio > 1) 
   {
/* get the exposure time scaling option */
   stat = SCKGETC("EXPOSU",1,3,&iav,output);
   output[4] = '\0';
   CGN_UPSTR(output);		                    /* convert to upper case */
   if (strcmp(output,"YES") == 0)		                 /* set flag */
      expscale = 1;
   else
      expscale = 0;

/* get the sigma flag */
   stat = SCKGETC("SIGMA",1,3,&iav,output);
   output[4] = '\0';
   CGN_UPSTR(output);		                    /* convert to upper case */
   if (strcmp(output,"YES") == 0)		                 /* set flag */
      sig = 1;
   else
      sig = 0;

/* get scaling by weight */
   stat = SCKGETC("WEIGHT",1,3,&iav,output);
   output[4] = '\0';
   CGN_UPSTR(output);		                    /* convert to upper case */
   if (strcmp(output,"YES") == 0)		                 /* set flag */
      weight = 1;
   else
      weight = 0;
   }

/* get statistics by mean/median/mode */
   stat = SCKGETC("STATIST",1,6,&iav,output);
   output[4] = '\0';
   CGN_UPSTR(output);		                    /* convert to upper case */
   if (strncmp(output,"MEA",3) == 0)		                 /* set flag */
      statist = 0;
   else if (strncmp(output,"MED",3) == 0)
      statist = 1;
   else if (strncmp(output,"MOD",3) == 0)
      statist = 2;
   else
      {
      SCTPUT(
      "*** INFO: Illegal statistical option (_STA): mode statistics assumed");
      statist = 2;
      }

/* get scaling by mean/median/mode */
   stat = SCKGETC("SCALE",1,3,&iav,output);
   output[4] = '\0';
   CGN_UPSTR(output);		                    /* convert to upper case */
   if (strcmp(output,"YES") == 0)		                 /* set flag */
      mmmscale = 1;
   else
      mmmscale = 0;

/* get offset by mean/median/mode */
   stat = SCKGETC("OFFSET",1,3,&iav,output);
   output[4] = '\0';
   CGN_UPSTR(output);		                    /* convert to upper case */
   if (strcmp(output,"YES") == 0)		                 /* set flag */
      mmmoffset = 1;
   else
      mmmoffset = 0;


/* get the range of valid pixels */
stat = SCKRDR("RANGE",1,2,&iav,&faux[2],&uni,&nulo);
if ((faux[2] != minundef) || (faux[3] != maxundef))
   {
   if (faux[2] > faux[3])
      SCETER(11,"*** FATAL: invalid valid_data_interval ...");
   else
      {
      low  = faux[2];
      high = faux[3];
      iaux[2] = 1;
      }
   }
else 
  {
  low = 0.0;
  high = 0.0; 
  }

if ((optio == 7) || (optio == 8) || (optio == 9))
   {
   stat = SCKRDR("CLIP",1,2,&iav,clip,&uni,&nulo);
   if (optio == 9)
     {
     if ( (clip[0] < 0.0) && (clip[1] < 0.0) )             /* INDEX interval */
       {
       iaux[1] = 0;
       iaux[3] = -CGN_NINT(clip[0]);
       iaux[4] = -CGN_NINT(clip[1]);
       }
     else                                                  /* DATA interval */
       {
       iaux[1] = 1;
       faux[0] = clip[0];
       faux[1] = clip[1];
       }
     }
   else
     {
     lowsigma  = clip[0];                              /* sigma clip factors */
     highsigma = clip[1];
     }
   }

/* get Null value */
stat = SCKGETC("BLANK",1,20,&iav,output);
if ((output[0] == '+') && (output[1] == '\0'))
   iaux[8] = 1;                                 /*  use `last' value as Null */
else
  {
  iav = CGN_CNVT(output,2,1,npixc,&usrnul,&aostep);
  if (iav < 1)
     SCETER(19,"*** FATAL: Invalid Null value ... ");
  }

/* test, if we have list of frames or catalog */
exist = 0;
inp = 0;                                          /* default, list of frames */
ll = CGN_INDEXS(line,".cat");
lu = CGN_INDEXS(line,".CAT");
if ((ll > 0) || (lu > 0)) inp = 1;                    /* input was catalogue */
ll = CGN_INDEXS(line,".tbl");
lu = CGN_INDEXS(line,".TBL");
if ((ll > 0) || (lu > 0)) inp = 2;                        /* input was table */

/* handle input from a catalog - get name of first image file in catalog */
if (inp == 1)                                             /* catalogue input */
   {
   if ((int) strlen(line) > 63)
      SCETER(3,"*** INFO: Catalogue name too long...");
   else
      strcpy(catfil,line);

   iseq = 0;
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)	       /* max. MAXIMS frames */
      {
      stat = SCCGET(catfil,0,frame[frmcnt],output,&iseq);
      if (frame[frmcnt][0] == ' ') break;	   /*  indicates the end ... */
      stat = SCFOPN(frame[frmcnt],D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
      }
   if (frmcnt <= 0) SCETER(4,error4);     /* there must be at least 1 image  */
   }

else if (inp == 2)                                            /* table input */
   { 
   if ((int) strlen(line) > 63)
      SCETER(3,"*** INFO: Table name too long...");
   else
      {
      strcpy(tblfil,line);
      stat = SCKGETC("ASSCOL", 1, 60, &iav, colnam);
      (void) TCTOPN(tblfil, F_I_MODE, &tid );     /* open table and read info*/
      (void) TCIGET(tid, &nrcol, &nrrow, &nsort, &allcol, &allrow); 
      if ( nrrow == 0 ) 
         SCETER(1, "*** INFO: No rows in this table (selected)");
      stat = TCCSER(tid, colnam, &colnum );         /* get the column number */
      if ( stat < 0 ) 
         SCETER(1, "*** FATAL: Column with calibration data not found");
      }
   }
else                                               /* input was frame string */
   {
   begin = 0;
   m = strlen(line);
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)
      {
      ll = CGN_EXTRSS(line,m,',',&begin,output,60);
      if (ll <= 0) break;
       
      CGN_FRAME(output,(int)F_IMA_TYPE,frame[frmcnt],0);  /* convert frames */
      stat = SCFOPN(frame[frmcnt],D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
      }
   if (frmcnt <= 0) SCETER(4,error4);     /* there must be at least 1 image  */
   }

sprintf(output,"Method=%s, low=%6.2f, high=%6.2f",meth,low,high);
SCTPUT(output);
if ((optio == 7) || (optio == 8))
   {
   sprintf(line,"lowclip=%6.2f, highclip=%6.2f",lowsigma,highsigma);
   SCTPUT(line);
   }
if (optio == 9)
   {
   if (iaux[1] == 1)
     sprintf(line,"    data  interval: low=%6.2g, high=%6.2g (rel. to median) to take mean",faux[0],faux[1]);
   else
     sprintf(line,"    index interval: low=%d, high=%d (rel. to median) to take mean",iaux[3],iaux[4]);
   SCTPUT(line);
   }
SCTPUT( 
       "------------------------------------------------------------------------------");
if (verbose == 0) 
   SCTPUT(
      "Output frame              Ncomb  Exp_time   Mean    Median     Mode   N_undef");

/* 
This part is special for the MIDAS Data Organizer Association Table 
*/ 
nrow = 1;                                     /* start loop trough the rows */
sect_tbl_start:
if (inp == 2)                                /* only in the case of a table */
   {
   stat = TCSGET(tid, nrow, &sel);                   /* check the selection */
   if (sel)                                                 /* row selected */
      {
      stat = TCERDC(tid, nrow, colnum, framec, &inul);   /* read the column */

      if (strcmp(&framec[0],"*") == 0) 
	 {
         stat = TCERDC(tid, nrow, 1, framec, &inul);     /* read sci column */
         strcat(framec,"_");
         strcat(framec,colnam);
         sprintf(line, "%s: '*' Option not implemented yet; skipped", framec);
         SCTPUT(line);
         goto sect_tbl_end;
         }

      begin = 0;
      m = strlen(framec);
      ll = CGN_EXTRSS(framec,m,'_',&begin,prefix,60);         /* name prefix */
      if (ll <= 0) 
	 {
         sprintf(line,"Output frame not or incorrectly defined; skipped");
         SCTPUT(line);
         goto sect_tbl_end;
         }

      stat = SCFINF(framec, 99, &ibuf);             /* check if input exists */
      if (stat == 0)
         {
         sprintf(line, "%s: Output frame already exists; skipped", framec);
         SCTPUT(line);
         exist++;
         goto sect_tbl_end;                          /* output already exist */
         }
      for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)
         {
         strcpy(frname,prefix);
         ll = CGN_EXTRSS(framec, m, '_', &begin, output, 60);
         if (ll <= 0) break; 
         lll = strlen(output);
         switch (lll)
	    { case 1:
                 strcat(frname,"000");
                 strcat(frname,output);                 
                 break;
              case 2:
                 strcat(frname,"00");
                 strcat(frname,output);                 
                 break;
              case 3:
                 strcat(frname,"0");
                 strcat(frname,output);                 
	         break;
	      case 4:
                 strcat(frname,output);                 
                 break;
	      default:
                 sprintf(line, 
                         "%s: Error in name of calibration frame; skipped", 
                         frname); 
                 SCTPUT(line);
                 goto sect_tbl_end;
            }                              

         CGN_FRAME(frname,(int)F_IMA_TYPE,frame[frmcnt],0);       /* convert */
         stat = SCFOPN(frame[frmcnt],D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
         if (stat < 0) 
	    {
            sprintf(line, "%s: Error in name of calibration frame; skipped", 
                    frname); 
            SCTPUT(line);
            goto sect_tbl_end;
	    }
         }
 
      if (frmcnt <= 0)                    /* there must be at least 1 image  */
	 {
         if (verbose == 1)  
	    {
            SCTPUT("No calibration frame found; skipped");
            goto sect_tbl_end;
	    }
         }
      }
   
   else
      goto sect_tbl_end;
   }
iaux[6] = frmcnt;

/* Check the exposure type */
if ((exptype[0] != '*') &&  (exptype[0] != '?') && (exptype[0] != '\0')) 
   {
   for (n=0; n< frmcnt; n++)
      {
      stat = SCDRDC(imnol[n],expname,1,1,20,&iav,expfram,&uni,&nulo);
      CGN_UPSTR(expfram);		         /* convert to upper case */
      if (stat != 0)
         {
         if (inp == 2) 
            sprintf(line,
                    "*** INFO: frame %s: Exp. type descriptor %s absent",
                    frame[n], expname);
         else
            sprintf(line,
                    "*** FATAL: frame %s: Exp. type descriptor %s absent",
                    frame[n], expname);
         SCTPUT(line);
         goto sect_tbl_end;

         }
      else if (strncmp(exptype,expfram,2) != 0) 
         { 
         if (inp == 2) 
	    {
            sprintf(line,
                    "*** WARNING: frame %s: Exposure type %s is wrong",
                    frame[n], expfram);
            SCTPUT(line);
	    }

         else
	    {
            sprintf(line,
                    "*** FATAL: frame %s: Exposure type %s is wrong",
                    frame[n], expfram);
            SCTPUT(line);
	    }
         goto sect_tbl_end;
         }
      }
   }

/*  --------------------------------------------*/
/*   get initial area from 1. frame             */
/*  --------------------------------------------*/
stat = SCDRDI(imnol[0],"NAXIS",1,1,&iav,&naxis[0],&uni,&nulo);
stat = SCDRDI(imnol[0],"NPIX",1,3,&iav,&npix[0][0],&uni,&nulo);
stat = SCDRDD(imnol[0],"START",1,3,&iav,&start[0][0],&uni,&nulo);
stat = SCDRDD(imnol[0],"STEP",1,3,&iav,&step[0][0],&uni,&nulo);
stat = SCDGETC(imnol[0],"CUNIT",1,64,&iav,cunit);

stat = SCDRDR(imnol[0],"LHCUTS",1,4,&iav,cutc,&uni,&nulo);
cutc[0] = cutc[3];
cutc[1] = cutc[2];

for (m=0; m<3; m++)
   {
   ostart[m] = start[0][m];
   ostep[m] = step[0][m];
   aostep = ostep[m];
   if (aostep < 0.0) aostep = - aostep;
   eps[m] = 0.0001 * aostep;	   		   /* take 0.01% of stepsize */
   oldend[m] = ostart[m] + (npix[0][m]-1)*ostep[m];
   }
 
/* now loop through the other input frames  */
for (n=1; n< frmcnt; n++)
   {
   for (m=0; m<3; m++)                                     /* init values... */
      {
      start[n][m] = 0.0;
      step[n][m] = 1.0;
      npix[n][m] = 1;
      }
  
   stat = SCDRDI(imnol[n],"NAXIS",1,1,&iav,&naxis[n],&uni,&nulo);
   stat = SCDRDI(imnol[n],"NPIX",1,3,&iav,&npix[n][0],&uni,&nulo);
   stat = SCDRDD(imnol[n],"START",1,3,&iav,&start[n][0],&uni,&nulo);
   stat = SCDRDD(imnol[n],"STEP",1,3,&iav,&step[n][0],&uni,&nulo);
   
/* stepsizes should have same sign and not differ too much...  */
   for (m=0; m<3; m++)
      {
      if ((ostep[m]*step[n][m]) <= 0.) SCETER(1,error3);
      aostep = step[n][m] - ostep[m];
      if (aostep < 0.) aostep = -aostep;
      
      if (aostep > eps[m]) SCETER(5,error1);
      }
 
/* get intersection or union of image areas  */
   if (action[0] != 'M') 
      {
      for (m=0; m<3; m++)
         {
         dif = start[n][m] + (npix[n][m]-1)*step[n][m];
         if (ostep[m] > 0.) 
            {
            if (ostart[m] < start[n][m]) ostart[m] = start[n][m];     /* MAX */
            if (oldend[m] > dif) oldend[m] = dif;		      /* MIN */
            }
         else
            {
            if (ostart[m] > start[n][m]) ostart[m] = start[n][m];     /* MIN */
            if (oldend[m] < dif) oldend[m] = dif;                     /* MAX */
            }
         }
      }
   else
      {
      for (m=0; m<3; m++)
         {
         dif = start[n][m] + (npix[n][m]-1)*step[n][m];
         if (ostep[m] < 0.) 
            {
            if (ostart[m] < start[n][m]) ostart[m] = start[n][m];     /* MAX */
            if (oldend[m] > dif) oldend[m] = dif;		      /* MIN */
            }
         else
            {
            if (ostart[m] > start[n][m]) ostart[m] = start[n][m];     /* MIN */
            if (oldend[m] < dif) oldend[m] = dif;                     /* MAX */
            }
         }
      }
 
   iaux[5] = 1;
   }
 
/* test, if something is left...  */
for (m=0; m<3; m++)
  if (ostep[m]*(oldend[m]-ostart[m]) < 0.) SCETER(2,error2);


 
/* read the the exposure time and number of combined images */
stat = SCECNT("GET",&ec,&el,&ed);
ln1 = 1; ln0 = 0;
stat = SCECNT("PUT",&ln1,&ln0,&ln0);

stat = SCKGETC("O_DESC",1,20,&iav,otime);
for (n=0; n< frmcnt; n++)
   {
   stat = SCDRDD(imnol[n],"O_TIME",7,1,&iav,&exptime[n],
                 &uni,&nulo);                               
   if ((stat != ERR_NORMAL) || (exptime[n] < 0.0))
      {
      exptime[n] = 0.0;
/*    
      if (verbose == 1)  
	 {
         sprintf(line,"frame %s: descr. %s missing: exp_time = 0.0 ",
                 frame[n],otime);
         SCTPUT(line);
         }
*/
      }

/* if weighting option, loop + look for descriptor NCOMBINE  */
   stat = SCDRDI(imnol[n],"NCOMBINE",1,1,&iav,&ncombine[n],
                 &uni,&nulo);                                     /* weights */
   if ((stat != ERR_NORMAL) || (ncombine[n] < 1)) 
      ncombine[n] = 1;

/* determine image section: set indices spix, epix, snpix, sstart, sstep */
   stat = SCKGETC("SECTION",1,40,&iav,mod_area);
   stat = Convcoo(1,imnol[n],mod_area,3,&nr,spix,epix);
   if (stat != 0)
     SCETER(1,"invalid coordinate input...");

   snpix[0] = (int)naxis[n];
   for (nr = 0; nr < snpix[0]; nr++)
     {
     snpix[1 + nr]  = (int)npix[n][nr];
     sstart[1 + nr] = (double)start[n][nr];
     sstep[1 + nr]  = (double)step[n][nr];
     }
   if ( (epix[0] - spix[0] + 1 != snpix[1])
     || (epix[1] - spix[1] + 1 != snpix[2]) )
     iact = 4;
   else
     iact = 5;                                            /* full plane */

   mmm[n]  = usrnul;
   wts[n]    = 1.0;
   scales[n] = 1.0;
   zeros[n]  = 0.0;
   w[n]      = 1.0;            /* dummy array; use in original average/image */
   }


/* now calculate the combined weigthing factors */
if (mmmscale == 1)                                            /* mmm scaling */
   {
   for (n=0; n< frmcnt; n++)                                  /* get the mmm */
      {
      if (verbose == 1)
	 {
         sprintf(line," ");
         SCTPUT(line);
         sprintf(line,"Statistics of input frame %s: ",frame[n]);
         SCTPUT(line);
         }
      stat = SCKWRR("OUTPUTR",Zeroes,1,12,&uni);            /* reset OUTPUTR */

      exmed(iact,imnol[n],zbins,snpix,spix,epix,&medx,stat);   /* exact median */
      Zstats(imnol[n], mod_area, snpix, zbins, frmstr, defaul1);
      stat = SCKRDR("OUTPUTR",1,12,&iav,statics,&uni,&nulo);/* statistics */
      stat = SCDWRR(imnol[n], "MEAN", &statics[2], 1 , 1, &uni);
      stat = SCDWRR(imnol[n], "MEDIAN", &statics[7], 1 , 1, &uni);
      stat = SCDWRR(imnol[n], "MODE", &statics[11], 1 , 1, &uni);
      stat = SCDWRR(imnol[n], "STATISTIC", statics, 1 , 12, &uni);
      if (statist == 0)
         stat = SCKRDR("OUTPUTR",3,1,&iav,&mmm[n],&uni,&nulo);       /* mode */
      else if (statist == 1)
         stat = SCKRDR("OUTPUTR",8,1,&iav,&mmm[n],&uni,&nulo);     /* median */
      else
         stat = SCKRDR("OUTPUTR",12,1,&iav,&mmm[n],&uni,&nulo);      /* mode */

      scales[n] = mmm[n];
      if (scales[n] <= 0.0)
	 {
         SCTPUT(
         "*** INFO: Mean/Median/Mode must be positive for scaling; no output");
         goto sect_tbl_end;
         }

      if (weight == 1)
         wts[n] = (float) (sqrt((float) (ncombine[n]*scales[n])));
      }
   }

else                                                       /* no mmm scaling */
   { 
   if (expscale == 1)                                    /* exposure scaling */
      {
	for (n=0; n< frmcnt; n++)
	  {
	  scales[n] = MAX(0.001,exptime[n]);                       /* scales */
	  if (weight == 1)                                     /* mmm offset */
	    wts[n] = (float)(sqrt((float) (ncombine[n]*scales[n])));
	                                                          /* weights */
	  }
      }
   if (mmmoffset == 1)
      {
      for (n=0; n< frmcnt; n++)
	 {
         if (verbose == 1)
            {
            sprintf(line," ");
            SCTPUT(line);
            sprintf(line,"Statistics of input frame %s: ",frame[n]);
            SCTPUT(line);
            }
	 stat = SCKWRR("OUTPUTR",Zeroes,1,12,&uni);         /* reset OUTPUTR */
	 exmed(iact,imnol[n],zbins,snpix,spix,epix,&medx,stat); /* exact median */
         Zstats(imnol[n], mod_area, snpix, sstart, frmstr, defaul1);
         stat = SCKRDR("OUTPUTR",1,12,&iav,statics,&uni,&nulo);/* statistics */
	 stat = SCDWRR(imnol[n], "MEAN", &statics[2], 1 , 1, &uni);
	 stat = SCDWRR(imnol[n], "MEDIAN", &statics[7], 1 , 1, &uni);
	 stat = SCDWRR(imnol[n], "MODE", &statics[11], 1 , 1, &uni);
	 stat = SCDWRR(imnol[n], "STATISTIC", statics, 1 , 12, &uni); 
         if (statist == 0)
            stat = SCKRDR("OUTPUTR",3,1,&iav,&mmm[n],&uni,&nulo);    /* mode */
         else if (statist == 1)
            stat = SCKRDR("OUTPUTR",8,1,&iav,&mmm[n],&uni,&nulo);  /* median */
         else
           stat = SCKRDR("OUTPUTR",12,1,&iav,&mmm[n],&uni,&nulo);    /* mode */
 
         zeros[n] = mmm[n]/scales[n];
         if (weight == 1)
            {
            if (zeros[n] <= 0.0)
	       {
               SCTPUT(
         "*** INFO: Mean/Median/Mode must be positive for scaling; no output");
               goto sect_tbl_end;
	       }
            else
               wts[n] = (float) (sqrt((float) 
                        (ncombine[n]*scales[n]/zeros[n])));
            }
         }
      } 
   }

/* change to relative scaling factors */
sumz = 0.0;
sums = 0.0;
sumw = 0.0;
for (n=0; n< frmcnt; n++)
  {
  sumz = sumz + zeros[n];
  sums = sums + scales[n];
  sumw = sumw + wts[n];
  }
meanz = sumz/(float) (frmcnt);
means = sums/(float) (frmcnt);

for (n=0; n< frmcnt; n++)
   {
   zeros[n]  = zeros[n] - meanz;
   scales[n] = scales[n]/means;
   zeros[n]  = zeros[n]*means;
   wts[n]    = wts[n]/sumw;
   }

/* check for the actual zeros and weights; if all equal don't scale */
scale = 0;
for (n=1; n<frmcnt; n++) 
  if ((zeros[n] != zeros[0]) || (wts[n] != wts[0]) || (scales[n] != scales[0]))
     scale = 1;

/* now take care of the output image */
nout     = 0;
exposure = 0.0;
for (n=0; n< frmcnt; n++)
   {
   nout     = nout + ncombine[n];
   exposure = exposure + wts[n]*exptime[n]/scales[n];
   }


/* new result frame with dimension as intersection or union of input frame */
naxisc = naxis[0];
if (action[0] == 'M') 
   {
   for (n=1; n<frmcnt; n++)
      {
      if (naxisc < naxis[n]) naxisc = naxis[n];			      /* MAX */
      }
   }
else
   {
   for (n=1; n<frmcnt; n++)
      {
      if (naxisc > naxis[n]) naxisc = naxis[n];                       /* MIN */
      }
   }
 

/* currently we only work on max 2-dim. frames...   */
if (naxisc > 2) 
   {
   SCTPUT("currently only 1-dim and 2-dim frames supported...");
   naxisc = 2;
   }


/* set up standard stuff for result frame        */
sizec = 1;
for (m=0; m<naxisc; m++)
   {
   startc[m] = ostart[m];
   stepc[m]  = ostep[m];
   npixc[m]  = CGN_NINT((oldend[m]-ostart[m]) / stepc[m]) + 1;
   sizec     = sizec * npixc[m];
   }

stat = SCFCRE(framec,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnoc);
stat = SCDWRI(imnoc,"NAXIS",&naxisc,1,1,&uni);
stat = SCDWRI(imnoc,"NPIX",npixc,1,naxisc,&uni);
stat = SCDWRD(imnoc,"START",startc,1,naxisc,&uni);
stat = SCDWRD(imnoc,"STEP",stepc,1,naxisc,&uni);
if (strcmp(mexp,"BS") == 0)	
   strcpy(ident,"CCD combined bias frame");
else if (strcmp(mexp,"DK") == 0)	                         /* set flag */
   strcpy(ident,"CCD combined dark frame");
else if (strcmp(mexp,"FF") == 0)	                         /* set flag */
   strcpy(ident,"CCD combined flat frame");
else if (strcmp(mexp,"SK") == 0)	                         /* set flag */
   strcpy(ident,"CCD combined sky frame");
else
   strcpy(ident,"CCD combined other frame");
stat = SCDWRC(imnoc,"IDENT",1,ident,1,72,&uni);
ln1  = (naxisc+1) * 16;
stat = SCDWRC(imnoc,"CUNIT",1,cunit,1,ln1,&uni);
  
/* in case of debugging save the counts in MIDAS image averdumy.dum  */
if (sig == 1)
   {
   strcpy(framecnt,framec);
   strcat(framecnt,"_cnt.bdf");
   strcpy(ident,"valid pixel counts");
   stat = SCFCRE(framecnt,D_I2_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnox);
   stat = SCDWRI(imnox,"NAXIS",&naxisc,1,1,&uni);
   stat = SCDWRI(imnox,"NPIX",npixc,1,naxisc,&uni);
   stat = SCDWRD(imnox,"START",startc,1,naxisc,&uni);
   stat = SCDWRD(imnox,"STEP",stepc,1,naxisc,&uni);
   stat = SCDWRC(imnox,"IDENT",1,ident,1,72,&uni);
   stat = SCDWRC(imnox,"CUNIT",1,cunit,1,ln1,&uni);   

   if (scale == 1)
      { 
      strcpy(framesig,framec);
      strcat(framesig,"_wtsig.bdf");
      strcpy(ident,"weighted sigma frame");
      }
   else
      { 
      strcpy(framesig,framec);
      strcat(framesig,"_sig.bdf");
      strcpy(ident,"sigma frame");
      }
   stat = SCFCRE(framesig,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnos);
   stat = SCDWRI(imnos,"NAXIS",&naxisc,1,1,&uni);
   stat = SCDWRI(imnos,"NPIX",npixc,1,naxisc,&uni);
   stat = SCDWRD(imnos,"START",startc,1,naxisc,&uni);
   stat = SCDWRD(imnos,"STEP",stepc,1,naxisc,&uni);
   stat = SCDWRC(imnos,"IDENT",1,ident,1,72,&uni);
   stat = SCDWRC(imnos,"CUNIT",1,cunit,1,ln1,&uni);
   }

/* see, if keyword MONITPAR(20) holds amount of virtual memory available */
stat = SCKRDI("MONITPAR",20,1,&iav,&chunkc,&uni,&nulo);
ln1 = (chunkc * chunkc)/5;     /* `chunkc' = x-dim square image to be mapped */

if (ln1 < npixc[0])
   stripe = 1;			                   /* at least a single line */
else
   stripe = ln1 / npixc[0];


/*  get virtual memory buffers */
if (npixc[1] < stripe) stripe = npixc[1];
chunkc = stripe * npixc[0];                        /* size of 1 result strip */

ll = chunkc * sizeof(float);                      /* here for the output map */
wpntr = osmmget((unsigned int)ll);
if (wpntr == (char *) 0)
   SCETER(66,"*** FATAL: Could not allocate virtual memory...");  
else
   pntrc = (float *) wpntr;


if ((optio == 7) || (optio == 8) || (sig == 1))        /* here for sigma map */
   { 
   wpntr = osmmget((unsigned int)ll);
   if (wpntr == (char *) 0)
      SCETER(66,"** FATAL: Could not allocate virtual memory...");  
   else
      pntrs = (float *) wpntr;                       /* here for sigma image */

   wpntr = osmmget((unsigned int)ll);
   if (wpntr == (char *) 0)
      SCETER(66,"*** FATAL: Could not allocate virtual memory...");  
   else
      pntrm = (float *) wpntr;                        /* here for mean image */
   }


for (n=0; n<frmcnt; n++)
   {
   if (npix[n][1] < stripe) 
      iav = npix[n][1];				             /* iav = y-dim  */
   else
      iav = stripe;

   chunk[n] = iav * npix[n][0];			     /* always stripe lines  */
   ll = chunk[n] * sizeof(float);
   wpntr = osmmget((unsigned int)ll);
   if (wpntr == (char *) 0)
      SCETER(66,"*** FATAL: Could not allocate virtual memory...");  
   else
      pntr[n] = (float *) wpntr;
   }
 

/*  now map chunk for z-direction and count buffer  */
n = frmcnt * chunkc;
ll = n * sizeof(float);
wpntr = osmmget((unsigned int)ll);
if (wpntr == (char *) 0)
   SCETER(66,"*** FATAL: Could not allocate virtual memory...");  
else
   zpntr = (float *) wpntr;

ll = chunkc * sizeof(short int);          /* here for the valid counts frame */
wpntr = osmmget((unsigned int)ll);
if (wpntr == (char *) 0)
   SCETER(66,"*** FATAL: Could not allocate virtual memory...");  
else
   xpntr = (short int *) wpntr;

 
/*  here the main loops over all chunks first fill the cube chunk,  
    work on it + store result                                                */

endc[0] = startc[0] + (npixc[0]-1)*stepc[0];
endc[1] = startc[1] + (stripe-1)*stepc[1];
zpix[0] = npixc[0];
  

/* print output for uses */
if (verbose == 1) 
  {
    SCTPUT(" " );
    SCTPUT( 
	   "------------------------------------------------------------------------------");
    if (statist == 0)
      sprintf(output,"Input frame                Ncomb   Exp_time     Mean    Scale   Offset   Weight");
    else if (statist == 1)
      sprintf(output,"Input frame                Ncomb   Exp_time     Median  Scale   Offset   Weight");
    else
      sprintf(output,"Input frame                Ncomb   Exp_time     Mode    Scale   Offset   Weight");

    SCTPUT(output);
    for (n=0; n< frmcnt; n++)
      {
	sprintf(line,"%-25s %6d  %8g %8g %8g %8g %8g",
		frame[n], ncombine[n], exptime[n], mmm[n], 
		1./scales[n], -zeros[n], wts[n]);
	SCTPUT(line);
      }
  }

for (nolin=0; nolin<npixc[1]; nolin+=stripe)
   {
   if ((nolin+stripe) > npixc[1])                       /* adjust chunk size */
      {
      stripe = npixc[1] - nolin;
      chunkc = stripe * npixc[0];
      for (n=0; n<frmcnt; n++)
         chunk[n] = stripe * npix[n][0];
      }
   zpix[1] = stripe;
 
   for (n=0; n<frmcnt; n++)
      {
 
      /*  convert start + end of overlap region into pixel no.'s   */

      for (m=0; m<naxisc; m++)
         {
         dif = (start[n][m]-startc[m]) / stepc[m];
         if (dif < 0.1) 
            cpix[m] = 0;
         else
            {
            cpix[m] = CGN_NINT(dif);			 /* offset in output */
            if (cpix[m] >= zpix[m]) 
               {
               iaux[0] = 0;
               goto sect_5360;
               }
            }

         dif = (startc[m]-start[n][m]) / step[n][m];
         if (dif < 0.1)					  /* offset in input */
            apix[m][0] = 0;
         else
            {
            apix[m][0] = CGN_NINT(dif);
            if (apix[m][0] >= npix[n][m]) 
               {
               iaux[0] = 0;
               goto sect_5360;
               }
            }
         
         dif = (endc[m]-start[n][m]) / step[n][m];
         apix[m][1] = CGN_NINT(dif);
         if (apix[m][1] >= npix[n][m]) apix[m][1] = npix[n][m] - 1;
         }
  
      iaux[0] = 1;
				  /* remember, SCFGET/SCFPUT begins with 1!  */
      felm = (apix[1][0] * npix[n][0]) + apix[0][0] + 1;
      tmpntr = (char *) pntr[n];
      stat = SCFGET(imnol[n],felm,chunk[n],&iav,tmpntr);
 

sect_5360:						    /* fill z-buffer */
      iaux[7] = n;
      faux[4] = w[n];
      fill(iaux,faux,pntr[n],xpntr,zpntr,apix,cpix,npix[n][0],zpix);
      }
  
   /*   now do the calculus  */
   if (optio == 1)                                             /* simple sum */
      ssum(iaux,faux,xpntr,zpntr,pntrc,usrnul,cutc,zpix,&ll);             

   else if (optio == 2)                                    /* simple average */
      if (scale == 1)    
         wtaver(iaux,faux,xpntr,zpntr,pntrc,scales,zeros,wts,
                usrnul,cutc,zpix,&ll);
      else
         aver(iaux,faux,xpntr,zpntr,pntrc,usrnul,cutc,zpix,&ll);

   else if (optio == 3)                                            /* median */
      if (scale == 1)                                            /* weighted */
         scmedian(iaux,faux,xpntr,zpntr,pntrc,scales,zeros, 
                usrnul,cutc,zpix,&ll);                           
      else                                                   /* not weighted */
         median(iaux,faux,xpntr,zpntr,pntrc,usrnul,cutc,zpix,&ll);

   else if (optio == 4)                                     /* min rejection */
      if (scale == 1)                                            /* weighted */
         wtminrej(iaux,faux,xpntr,zpntr,pntrc,scales,zeros,wts,    
                usrnul,cutc,zpix,&ll);                           
      else                                                   /* not weighted */
         minrej(iaux,faux,xpntr,zpntr,pntrc,usrnul,cutc,zpix,&ll);

   else if (optio == 5)                                     /* max rejection */
      if (scale == 1)                                            /* weighted */
         wtmaxrej(iaux,faux,xpntr,zpntr,pntrc,scales,zeros,wts,
                usrnul,cutc,zpix,&ll);                           
      else                                                   /* not weighted */
         maxrej(iaux,faux,xpntr,zpntr,pntrc,usrnul,cutc,zpix,&ll); 

   else if (optio == 6)                                 /* min max rejection */
      if (scale == 1)                                            /* weighted */
         wtmmrej(iaux,faux,xpntr,zpntr,pntrc,scales,zeros,wts,
                usrnul,cutc,zpix,&ll);                           
      else                                                   /* not weighted */
         mmrej(iaux,faux,xpntr,zpntr,pntrc,usrnul,cutc,zpix,&ll);
  
   else if (optio == 7)                                    /* sigma clipping */
      if (scale == 1)                                            /* weighted */
         wtsigclip(iaux,faux,xpntr,zpntr,pntrc,pntrs,pntrm,
                   scales,zeros,wts,lowsigma,highsigma,usrnul,cutc,zpix,&ll);
      else                                                   /* not weighted */
         sigclip(iaux,faux,xpntr,zpntr,pntrc,pntrs,pntrm,
                 lowsigma,highsigma,usrnul,cutc,zpix,&ll);                  

   else if (optio == 8)                            /* average sigma clipping */
      if (scale == 1)                                            /* weighted */
         wtavsigclip(iaux,faux,xpntr,zpntr,pntrc,pntrs,pntrm,
                     scales,zeros,wts,lowsigma,highsigma,usrnul,cutc,zpix,&ll);
      else                                                   /* not weighted */
         avsigclip(iaux,faux,xpntr,zpntr,pntrc,pntrs,pntrm,
                   lowsigma,highsigma,usrnul,cutc,zpix,&ll);
   else if (optio == 9)                                       /* mean median */
      if (scale == 1)                                            /* weighted */
         scmedian(iaux,faux,xpntr,zpntr,pntrc,scales,zeros, 
                usrnul,cutc,zpix,&ll);                           
      else                                                   /* not weighted */
         median(iaux,faux,xpntr,zpntr,pntrc,usrnul,cutc,zpix,&ll);



   /*  and write results to disk */
   felm = nolin*npixc[0] + 1;
   tmpntr = (char *) pntrc;
   stat = SCFPUT(imnoc,felm,chunkc,tmpntr);
   nulcnt += ll;				        /* update null count */


/* create debug image if wanted */
   if (sig == 1)
      {
      felm   = nolin*npixc[0] + 1;
      tmpntr = (char *) xpntr;
      stat   = SCFPUT(imnox,felm,chunkc,tmpntr);

      if (scale == 1)
         wtsigma(iaux,faux,xpntr,zpntr,pntrc,pntrs,
                 scales,zeros,wts,usrnul,cuts,zpix,&ll);
      else
         sigma(iaux,faux,xpntr,zpntr,pntrc,pntrs,
               usrnul,cuts,zpix,&ll);

      felm   = nolin*npixc[0] + 1;
      tmpntr = (char *) pntrs;
      stat   = SCFPUT(imnos,felm,chunkc,tmpntr);
      }

/* follow chunks with start + end value  */
   startc[1] = startc[1] + stripe*stepc[1];
   endc[1]   = endc[1] + stripe*stepc[1];
   }

/*  update descriptors + cuts of result frame */
null[0] = ' ';
null[1] = '\0';
CGN_DSCUPD(imnol[0],imnoc,null);
cutc[2] = cutc[0];
cutc[3] = cutc[1];
stat = SCDWRR(imnoc,"LHCUTS",cutc,1,4,&uni);
stat = SCDWRI(imnoc,"NCOMBINE",&nout,1,1,&uni);
stat = SCDWRR(imnoc,"O_TIME",&exposure,7,1,&uni);

/*  update descriptors + cuts of error and count frame */
if (sig == 1)
   {
   cutn[0] = cuts[2];
   cutn[1] = cuts[3];
   cutn[2] = cuts[2];
   cutn[3] = cuts[3];
   cuts[2] = cuts[0];
   cuts[3] = cuts[1];

   null[0] = ' ';
   null[1] = '\0';
   CGN_DSCUPD(imnol[0],imnos,null);
   stat = SCDWRR(imnos,"LHCUTS",cuts,1,4,&uni);

   null[0] = ' ';
   null[1] = '\0';
   CGN_DSCUPD(imnol[0],imnox,null);
   stat = SCDWRR(imnox,"LHCUTS",cutn,1,4,&uni);
   }

cutc[0] = nulcnt;			                  /* update key NULL */
if (iaux[8] == 0)
   {
   cutc[1] = usrnul;
   ln1 = 2;
   }
else
   ln1 = 1;
stat = SCKWRR("NULL",cutc,1,ln1,&uni);

snpix[0] = naxisc;
for (nr = 0; nr < snpix[0]; nr++)
    snpix[1 + nr] = npixc[nr];

stat = SCKWRR("OUTPUTR",Zeroes,1,12,&uni);               /* reset OUTPUTR */
if (verbose == 1)
   {
   SCTPUT(" ");
   sprintf(line,"Statistics of output frame %s",framec);
   SCTPUT(line);
   exmed(iact,imnoc,zbins,snpix,spix,epix,&medx,stat);   /* exact median */
   Zstats(imnoc, im_area, snpix, zbins, frmstr, defaul1);
   }
else{
   exmed(iact,imnoc,zbins,snpix,spix,epix,&medx,stat);   /* exact median */
   Zstats(imnoc, im_area, snpix, zbins, frmstr, defaul1);
}
stat = SCKRDR("OUTPUTR",1,12,&iav,statics,&uni,&nulo);   /* statistics */
if (verbose == 1) 
   {
   SCTPUT(" ");
   SCTPUT(
     "Output frame              Ncomb  Exp_time   Mean    Median    Mode   N_undef");

   sprintf(line,"%-25s %3d  %8g %8g %8g %8g %8d", 
           framec, nout, exposure, statics[2], statics[7], statics[11], nulcnt);
   SCTPUT(line);
   SCTPUT( 
       "------------------------------------------------------------------------------");
   }
else
   {
   sprintf(line,"%-25s %3d  %8g %8g %8g %8g %8d", 
           framec, nout, exposure, statics[2], statics[7], statics[11], nulcnt);
   SCTPUT(line);
   }

/* Write the statistics in the output frame */
stat = SCKRDR("OUTPUTR",1,9,&iav,&statics[0],&uni,&nulo);   /* statistics */
stat = SCDWRR(imnoc, "MEAN", &statics[2], 1 , 1, &uni);
stat = SCDWRR(imnoc, "MEDIAN", &statics[7], 1 , 1, &uni);
stat = SCDWRR(imnoc, "MODE", &statics[11], 1 , 1, &uni);
for (n=0; n< frmcnt; n++){
  SCDWRC(imnol[n],"CCDMETHOD",strlen(meth),meth,1,1,&uni);
  SCDWRR(imnol[n],"CCDMMM",&mmm[n],1,1,&uni);
  SCDWRR(imnol[n],"CCDSCALE",&scales[n],1,1,&uni);
  SCDWRR(imnol[n],"CCDOFFSET",&zeros[n],1,1,&uni);
  SCDWRR(imnol[n],"CCDWEIGHT",&wts[n],1,1,&uni);
}

icount++;                           /* increase the number of output frames */

sect_tbl_end:
if (inp == 2)
   {
   nrow ++;
   if (nrow <= nrrow) goto sect_tbl_start;
   }

if (icount > 0) 
   {
   SCTPUT(
       "------------------------------------------------------------------------------");
   if (iaux[8] == 0)
      sprintf(
        output,"Undefined pixels, set to `null value' (= %12.6f)",usrnul);
   else
      sprintf(output,"Undefined pixels, set to `previous pixel'");
   SCTPUT(output);
   if (exist > 0) 
      SCTPUT(
   "*** INFO: One of more output frames already existed; no new ones created");
   }

else
   {
   if (exist > 0) 
      SCTPUT(
      "*** INFO: Your output frame(s) already existed; no new one(s) created");
   else      
      SCTPUT(
         "*** INFO: No output frame(s) are created; check your input specs");
   }
SCSEPI();
}



/*

*/

fill(iaux,faux,a,x,z,apix,cpix,npixa,npixc)
 
int         *iaux, apix[3][2], *cpix;
int         npixa, *npixc;
 
short int   *x;

float       *faux, *a, *z;

{
 
int         frmcnt, count;
int         n, nn, nin, nini, nout, jump, indx, cntdx;
int         nx, ny, xfirst, yfirst, xend, yend;

register float  rr, ww;

 

/* ---------------------------
the intermediate z-space is structured as follows:
   
pix(1) pix(2) up to  pix(FRMCNT)      1              )
pix(1) pix(2) up to  pix(FRMCNT)      2              )
. . .                                                ) = 1 line
pix(1) pix(2) up to  pix(FRMCNT)      NPIXC(1)       )
  
STRIP lines of above
--------------------------- */
 

/*   if here for the first time, init the count `pixels'  */
ww     = faux[4];		                     /* get weighting factor */
frmcnt = iaux[6];
count  = iaux[7];

if (count == 0)
   {
   nn = npixc[0] * npixc[1];
   if ((iaux[5] == 0) && (iaux[2] == 0))
      nini = frmcnt;			     /* for `nomerge' and `all data' */
   else
      nini = 0;

   for (n=0; n<nn; n++)
      x[n] = nini;
   }

if (iaux[0] == 0) return;		   /* if iaux[0] = 0, not much to do */

nini = 0;			   /* input frame begins with 1. valid pixel */
if (iaux[5] == 0)
   {

/* --------------------------------- */
/*  here for the non-merging option  */
/* --------------------------------- */

   if (iaux[2] == 0)		                          /* take all pixels */
      {
      indx = count;
      for (nn=0; nn<npixc[1]; nn++)
         {
         nin = nini;
         for (n=0; n<npixc[0]; n++)
            {
            z[indx] = ww * a[nin++];
            indx += frmcnt;
            }
         nini += npixa;     	                  /* follow with input index */
         }
      }
   else			               /* take only pixels in valid interval */
      {
      cntdx = 0;
      indx = 0;
      for (nn=0; nn<npixc[1]; nn++)
         {
         nin = nini;
         for (n=0; n<npixc[0]; n++)
            {
            rr = a[nin++];
            if ((rr >= faux[2]) && (rr <= faux[3]))
               {
               nout = indx + x[cntdx];
               z[nout] = ww * rr;
               x[cntdx] ++;			          /* increment count */
               }
            indx += frmcnt;
            cntdx ++;
            }
         nini += npixa;     	                  /* follow with input index */
         }
      }
   }
else
   {


/* ----------------------------- */
/*  here for the merging option  */
/* ----------------------------- */

   nx = apix[0][1] - apix[0][0];		  /*   get overlapping part  */
   ny = apix[1][1] - apix[1][0];
   xfirst = cpix[0];
   yfirst = cpix[1];
   xend = xfirst + nx;
   yend = yfirst + ny;
   jump = frmcnt * npixc[0];

   cntdx = 0;
   indx = 0;

   if (iaux[2] == 0)		                          /* take all pixels */
      {
      for (nn=0; nn<npixc[1]; nn++)
         {
         if ((nn >= yfirst) && (nn <= yend)) 
            {
            nin = nini;
            for (n=0; n<npixc[0]; n++)
               {
               if ((n >= xfirst) && (n <= xend)) 
                  {
                  nout = indx + x[cntdx];
                  z[nout] = ww * a[nin++];
                  x[cntdx] ++;			          /* increment count */
                  }
               indx += frmcnt;
               cntdx ++;
               }

            nini += npixa;			  /* follow with input index */
            }
         else
            {
            indx += jump;			       /* jump a whole line; */
            cntdx += npixc[0];
            }
         }
      }
   else			/* take only pixels in valid interval */
      {
      for (nn=0; nn<npixc[1]; nn++)
         {
         if ((nn >= yfirst) && (nn <= yend))
            {
            nin = nini;
            for (n=0; n<npixc[0]; n++)
               {
               if ((n >= xfirst) && (n <= xend))
                  {
                  rr = a[nin++];
                  if ((rr >= faux[2]) && (rr <= faux[3]))
                     {
                     nout = indx + x[cntdx];
                     z[nout] = ww * rr;
                     x[cntdx] ++;                         /* increment count */
                     }
                  }
               indx += frmcnt;
               cntdx ++;
               }

            nini += npixa;                        /* follow with input index */
            }
         else
            {
            indx += jump;                              /* jump a whole line; */
            cntdx += npixc[0];
            }
         }
      }
   }
}




/*

*/
ssum(iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
/*
Here the simple summing of all deata values 
*/
      
int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;

{
int        frmcnt, n, nn, mcc, indx, cntdx, size;
register float       rval;
double     sum;

static float old = 0.0;


/*  init  */
mcc = 0;
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];
  
for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn = x[cntdx];
  
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 
   else
      {
      sum = 0.0;
      for (n=indx; n<indx+nn; n++)
         sum += z[n];
         rval = sum;
      }

   c[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}



/*

*/
wtaver(iaux,faux,x,z,c,scales,zeros,wts,usrnul,cuts,npixc,nc)
/*
Here the averaging with weighting
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;
float      *zeros, *wts, *scales;

{
int        frmcnt, n, nn;
int        mcc, indx, cntdx, size;
int        idx;
register float       rval;
double     sum;

static float old = 0.0;


/*  init  */
mcc = 0;
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];
  
for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn = x[cntdx];
  
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 
   else
      {
      sum = 0.0;
      for (n=indx; n<indx+nn; n++) 
	 {
         idx = n-indx;
         sum += wts[idx]*(z[n]/scales[idx]-zeros[idx]);
         }
      rval = sum;
      }

   c[cntdx] = rval;
   old = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}

         
/*

*/
aver(iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
/*
Here the simple averaging without weighting
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;

{
int        frmcnt, n, nn, mcc, indx, cntdx, size;
register float       rval;
double     sum;

static float old = 0.0;


/*  init  */
mcc = 0;
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];
  
for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn = x[cntdx];
  
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 
   else
      {
      sum = 0.0;
      for (n=indx; n<indx+nn; n++)
         sum += z[n];
         rval = sum/nn;
      }

   c[cntdx] = rval;
   old = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}




/*

*/
scmedian(iaux,faux,x,z,c, scales,zeros,usrnul,cuts,npixc,nc)
/*
Here the median with scaling
*/
      
int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;
float      *zeros, *scales;

{
int        frmcnt, n, nn, nh, k, mcc, ma, mb, idx, indx, cntdx, size;
register float       rval, va, vb;
float      rbuf[MAXIMSONE];			       /* buffer for sorting */
static float old = 0.0;
double  sum;

/*  init  */
mcc = 0;
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];
      
/*
............................................................


here comes the code for handling the median
if iaux[1] = 0, use index interval [ iaux[3],iaux[4] ]
           = 1, use data interval [ faux[0],faux[1] ]

............................................................
								*/

if (iaux[1] == 0)
   {
   k = iaux[3] + iaux[4];         /* k==0 if both [3] == [4] = 0 */
   if (k > 0)                     /* take average over index interval  */
      {                           /* [ MEDIAN-iaux[3] , MEDIAN+iaux[4] ] */
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];
   
         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                         /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx]/scales[0]-zeros[0];
            else if (nn == 2)
	      {
	      rbuf[0] = z[indx]/scales[0]-zeros[0];
	      rval    = z[indx+1]/scales[1]-zeros[1];
	      if (rval > rbuf[0]) rval = rbuf[0];        /* take min of them */
	      }
            else                		/* here we have: nn > 2 */
               {
               k = 1;				/* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
		 {
		 idx = n-indx;
		 rbuf[k++] = z[n]/scales[idx]-zeros[idx];
		 }
               sortr(rbuf,nn);			/* sort array */
               nh = (nn+1) / 2;	     /* index of median (starting at 1 ...) */
               ma = nh - iaux[3];
               if (ma < 1) ma = 1;
               mb = nh + iaux[4];
               if (mb > nn) mb = nn;
               nn = mb - ma + 1;
               sum = rbuf[ma];
               for (k=ma+1; k<=mb; k++)
                   sum += rbuf[k];
               rval = sum / nn;
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }
   else                       /* DEFAULT: no averaging, take only the MEDIAN */
      {  
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];

         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                      /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx]/scales[0]-zeros[0];
            else if (nn == 2)
	      {
	      rbuf[0] = z[indx]/scales[0]-zeros[0];
	      rval    = z[indx+1]/scales[1]-zeros[1];
	      if (rval > rbuf[0]) rval = rbuf[0];       /* take min of them */
	      }
            else		                /* here we have: nn > 2 */
               {
               k = 1;                           /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
		 {
		 idx = n-indx;
		 rbuf[k++] = z[n]/scales[idx]-zeros[idx];
		 }
               sortr(rbuf,nn);                 /* sort array */
               nh = (nn+1) / 2;      /* index of median (starting at 1 ...) */
               rval = rbuf[nh];
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }
   }

else                                /* do the averaging in a data interval */
   {			
   rval = faux[0] + faux[1];
   if (rval > 0.0)                     /* take average over data interval  */
      {                           /* [ MEDIAN-faux[0] , MEDIAN+faux[4] ] */
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];

         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                         /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx];
            else if (nn == 2)
	      {
	      rbuf[0] = z[indx]/scales[0]-zeros[0];
	      rval    = z[indx+1]/scales[1]-zeros[1];
	      if (rval > rbuf[0]) rval = rbuf[0];       /* take min of them */
	      }
            else                                /* here we have: nn > 2 */
               {
               k = 1;                           /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
		 {
		   idx = n-indx;
		   rbuf[k++] = z[n]/scales[idx]-zeros[idx];
		 }
               sortr(rbuf,nn);                 /* sort array */
               nh = (nn+1) / 2;      /* index of median (starting at 1 ...) */
               va = rbuf[nh] - faux[0];
               vb = rbuf[nh] + faux[1];
               sum = 0.0;
               nh = 0;
               for (k=1; k<=nn; k++)
                  {
                  rval = rbuf[k];
                  if (rval > vb) break;		/* array is sorted ... */
                  if (rval >= va) 
                     {
                     sum += rval;
                     nh ++;
                     }
                  }
               rval = sum / nh;
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }

   else					/* just use median */
      {
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];

         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                         /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx]/scales[0]-zeros[0];
            else if (nn == 2)
	      {
	      rbuf[0] = z[indx]/scales[0]-zeros[0];
	      rval    = z[indx+1]/scales[1]-zeros[1];
	      if (rval > rbuf[0]) rval = rbuf[0];       /* take min of them */
	      }
            else                                /* here we have: nn > 2 */
               {
               k = 1;                           /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
		 {
		 idx = n-indx;
		 rbuf[k++] = z[n]/scales[idx]-zeros[idx];
		 }
               sortr(rbuf,nn);                 /* sort array */
               nh = (nn+1) / 2;      /* index of median (starting at 1 ...) */
               rval = rbuf[nh];
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }
   }

*nc = mcc;
return;
}

/*

*/
oldscmedian(iaux,faux,x,z,c, scales,zeros,usrnul,cuts,npixc,nc)
/*
Here the median with scaling
*/
      
int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;
float      *zeros, *scales;

{ 
int        frmcnt, n, nn, nh, k, mcc, idx, indx, cntdx, size;
register float       rval;
float      rbuf[MAXIMSONE];			/* buffer for sorting */
float      sum;

static float old = 0.0;

/*  init  */
mcc    = 0;
indx   = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];
      
for (cntdx=0; cntdx<size; cntdx++)              /* loop through count buffer */
   {
   nn = x[cntdx];

   if (nn == 0)
      {
      if (iaux[8] == 1)
         rval = old;                                       /* use last value */
      else
         rval = usrnul;
         mcc ++;                                    /* increment null count  */
      }
   else
      {
      if (nn == 1)
         rval = z[indx]/scales[0]-zeros[0];
      else if (nn == 2)
         {
         rbuf[0] = z[indx]/scales[0]-zeros[0];
         rval    = z[indx+1]/scales[1]-zeros[1];
         if (rval > rbuf[0]) rval = rbuf[0];	         /* take min of them */
         }
      else		                             /* here we have: nn > 2 */
         {
         k   = 1;                                    /* start at index 1 ... */
         for (n=indx; n<indx+nn; n++)
  	    {
            idx = n-indx;
            rbuf[k++] = z[n]/scales[idx]-zeros[idx];
	    }
         sortr(rbuf,nn);                                      /* sort array */
         nh = (nn+1) / 2; 
         rval = rbuf[nh];
         }
      }

   c[cntdx] = rval;
   old   = rval;
   if (cuts[0] > rval) cuts[0] = rval;            /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}
     


/*

*/
median(iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
/*
Here the median without scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;

{
int        frmcnt, n, nn, nh, k, mcc, ma, mb, indx, cntdx, size;
register float       rval, va, vb;
float      rbuf[MAXIMSONE];			       /* buffer for sorting */
static float old = 0.0;
double  sum;

/*  init  */
mcc = 0;
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];
      
/*
............................................................


here comes the code for handling the median
if iaux[1] = 0, use index interval [ iaux[3],iaux[4] ]
           = 1, use data interval [ faux[0],faux[1] ]

............................................................
								*/

if (iaux[1] == 0)
   {
   k = iaux[3] + iaux[4];         /* k==0 if both [3] == [4] = 0 */
   if (k > 0)                     /* take average over index interval  */
      {                           /* [ MEDIAN-iaux[3] , MEDIAN+iaux[4] ] */
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];
   
         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                         /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx];
            else if (nn == 2)
               rval = (z[indx] + z[indx+1]) / 2;
            else                		/* here we have: nn > 2 */
               {
               k = 1;				/* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
                  rbuf[k++] = z[n];
               sortr(rbuf,nn);			/* sort array */
               nh = (nn+1) / 2;	     /* index of median (starting at 1 ...) */
               ma = nh - iaux[3];
               if (ma < 1) ma = 1;
               mb = nh + iaux[4];
               if (mb > nn) mb = nn;
               nn = mb - ma + 1;
               sum = rbuf[ma];
               for (k=ma+1; k<=mb; k++)
                   sum += rbuf[k];
               rval = sum / nn;
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }
   else                       /* DEFAULT: no averaging, take only the MEDIAN */
      {  
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];

         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                      /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx];
            else if (nn == 2)
               {
               rbuf[0] = z[indx];
               rval = z[indx+1];
               if (rval > rbuf[0]) rval = rbuf[0];	/* take min of them */
               }
            else		                /* here we have: nn > 2 */
               {
               k = 1;                           /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
                  rbuf[k++] = z[n];
               sortr(rbuf,nn);                 /* sort array */
               nh = (nn+1) / 2;      /* index of median (starting at 1 ...) */
               rval = rbuf[nh];
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }
   }

else                                /* do the averaging in a data interval */
   {			
   rval = faux[0] + faux[1];
   if (rval > 0.0)                     /* take average over data interval  */
      {                           /* [ MEDIAN-faux[0] , MEDIAN+faux[4] ] */
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];

         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                         /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx];
            else if (nn == 2)
               rval = (z[indx] + z[indx+1]) / 2;
            else                                /* here we have: nn > 2 */
               {
               k = 1;                           /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
                  rbuf[k++] = z[n];
               sortr(rbuf,nn);                 /* sort array */
               nh = (nn+1) / 2;      /* index of median (starting at 1 ...) */
               va = rbuf[nh] - faux[0];
               vb = rbuf[nh] + faux[1];
               sum = 0.0;
               nh = 0;
               for (k=1; k<=nn; k++)
                  {
                  rval = rbuf[k];
                  if (rval > vb) break;		/* array is sorted ... */
                  if (rval >= va) 
                     {
                     sum += rval;
                     nh ++;
                     }
                  }
               rval = sum / nh;
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }

   else					/* just use median */
      {
      for (cntdx=0; cntdx<size; cntdx++)        /* loop through count buffer */
         {
         nn = x[cntdx];

         if (nn == 0)
            {
            if (iaux[8] == 1)
               rval = old;                         /* use last value */
            else
               rval = usrnul;
            mcc ++;                             /* increment null count  */
            }
         else
            {
            if (nn == 1)
               rval = z[indx];
            else if (nn == 2)
	      {
              rbuf[0] = z[indx];
	      rval = z[indx+1];
	      if (rval > rbuf[0]) rval = rbuf[0];	/* take min of them */
	      }
            else                                /* here we have: nn > 2 */
               {
               k = 1;                           /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
                  rbuf[k++] = z[n];
               sortr(rbuf,nn);                 /* sort array */
               nh = (nn+1) / 2;      /* index of median (starting at 1 ...) */
               rval = rbuf[nh];
               }
            }

         c[cntdx] = rval;
         old = rval;
         if (cuts[0] > rval) cuts[0] = rval;       /* update min + max */
         if (cuts[1] < rval) cuts[1] = rval;
         indx += frmcnt;
         }
      }
   }

*nc = mcc;
return;
}

/*

*/
oldmedian(iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
/*
Here the median without scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;

{
int        frmcnt, n, nn, nh, k, mcc, indx, cntdx, size;
register float       rval;
float      rbuf[MAXIMSONE];			       /* buffer for sorting */

static float old = 0.0;

/*  init  */
mcc = 0;
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];
      
for (cntdx=0; cntdx<size; cntdx++)              /* loop through count buffer */
   {
   nn = x[cntdx];

   if (nn == 0)
      {
      if (iaux[8] == 1)
         rval = old;                                       /* use last value */
      else
         rval = usrnul;
         mcc ++;                                    /* increment null count  */
      }
   else
      {
      if (nn == 1)
         rval = z[indx];
      else if (nn == 2)
         {
         rbuf[0] = z[indx];
         rval = z[indx+1];
         if (rval > rbuf[0]) rval = rbuf[0];	         /* take min of them */
         }
      else		                             /* here we have: nn > 2 */
         {
         k = 1;                                      /* start at index 1 ... */
         for (n=indx; n<indx+nn; n++)
            rbuf[k++] = z[n];
         sortr(rbuf,nn);                                      /* sort array */
         nh = (nn+1) / 2;             /* index of median (starting at 1 ...) */
         rval = rbuf[nh];
         }
      }

   c[cntdx] = rval;
   old   = rval;
   if (cuts[0] > rval) cuts[0] = rval;            /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}
      



/*

*/
wtminrej(iaux,faux,x,z,c,scales,zeros,wts,usrnul,cuts,npixc,nc)
/*
Here the minimum rejection with scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;
float      *zeros, *wts, *scales;

{
int        frmcnt, k, n, nn, mcc, idx, indx, cntdx, size;
register float       rval;
float      mean, val, minval, minwt, wt;

static float old = 0.0;


/*  init  */
mcc    = 0;
indx   = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];
  
for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      {
      mean   = 0.0;
      minval = z[indx]/scales[0] - zeros[0];
      minwt  = wts[0];
      k      = indx;

      for (n=indx+1; n<indx+nn; n++) 
         {
         idx = n-indx;
         val = z[n]/scales[idx]-zeros[idx];
         wt  = wts[idx];

         if (val < minval)
	    {
            mean  += minwt*minval;
            minval = val;
            minwt  = wt;
            k      = n;
            }
         else
            mean  += wt*val;
         }

      rval = mean/(1.- minwt);
      z[k] = usrnul;
      }

   c[cntdx] = rval;
   old = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}



/*

*/
minrej(iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
/*
Here the minimum rejection without scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;

{
int        frmcnt, k, n, nn, mcc, indx, cntdx, size;
int        nims;
float      minval, val, mean;
register float       rval;

static float old = 0.0;


/*  init  */
mcc    = 0;
indx   = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];
nims   = frmcnt - 1;

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 
   else
      {
      mean   = 0.0;
      minval = z[indx];
      k      = indx;

      for (n=indx+1; n<indx+nn; n++) 
         {
         val  = z[n];
         if (val < minval)
	    {
            mean  += minval;
            minval = val;
            k      = n;
            }
         else
            mean  += val;

         }
      rval = mean/nims;
      z[k] = usrnul;
      }

   c[cntdx] = rval;
   old = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}



/*

*/
wtmaxrej(iaux,faux,x,z,c,scales,zeros,wts,usrnul,cuts,npixc,nc)
/*
Here the maximum rejection with scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;
float      *zeros, *wts, *scales;

{
int        frmcnt, k, n, nn, mcc, idx, indx, cntdx, size;
register float       rval;
float      mean, val, maxval, maxwt, wt;

static float old = 0.0;


/*  init  */
mcc    = 0;
indx   = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];
  
for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      {
      mean   = 0.0;
      maxval = z[indx]/scales[0] - zeros[0];
      maxwt  = wts[0];
      k      = indx;

      for (n=indx+1; n<indx+nn; n++) 
         {
         idx = n-indx;
         val = z[n]/scales[idx]-zeros[idx];
         wt  = wts[idx];

         if (val > maxval)
	    {
            mean  += maxwt*maxval;
            maxval = val;
            maxwt  = wt;
            k      = n;

            }
         else
            mean += wt*val;

         }
      rval = mean/(1.- maxwt);
      z[k] = usrnul;
      }

   c[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx    += frmcnt;
   }

*nc = mcc;
return;
}



/*

*/
maxrej(iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
/*
Here the maximum rejection without scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;

{
int        frmcnt, k, n, nn, mcc, indx, cntdx, size;
int        nims;
float      maxval, val, mean;
register float       rval;

static float old = 0.0;


/*  init  */
mcc    = 0;
indx   = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];
nims   = frmcnt - 1;

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      {
      mean   = 0.0;
      maxval = z[indx];
      k      = indx;

      for (n=indx+1; n<indx+nn; n++) 
         {
         val  = z[n];
         if (val > maxval)
	    {
            mean  += maxval;
            maxval = val;
            k      = n;
            }
         else
            mean  += val;

         }
      rval = mean/nims;
      z[k] = usrnul;
      }

   c[cntdx] = rval;
   old = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}

         
/*

*/
wtmmrej(iaux,faux,x,z,c,scales,zeros,wts,usrnul,cuts,npixc,nc)
/*
Here the minimum, maximum rejection with scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;
float      *zeros, *wts, *scales;

{
int        frmcnt, k, l, n, nn, mcc, idx, indx, cntdx, size;
register float       rval;
float      mean, val, minval, maxval;
float      wt, minwt, maxwt;

static float old = 0.0;


/*  init  */
mcc    = 0;
indx   = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];
  
for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      {
      mean   = 0.0;
      minval = z[indx]/scales[0] - zeros[0];
      maxval = z[indx+1]/scales[1] - zeros[1];

      if (minval < maxval)
 	 {
         minwt  = wts[0];
         maxwt  = wts[1];
         k      = indx;
         l      = indx+1;
         }

      else
	 {
         val    = minval;
         minval = maxval;
         maxval = val;
         minwt  = wts[1];
         maxwt  = wts[0];
         k      = indx+1;
         l      = indx;
         }

      for (n=indx+2; n<indx+nn; n++) 
         {
         idx = n-indx;
         val  = z[n]/scales[idx]-zeros[idx];
         wt   = wts[idx];

         if (val < minval)
            {
            mean  += minwt*minval;
            minval = val;
            minwt  = wt;
            k      = n;
            }

         else if (val > maxval)
            {
            mean  += maxwt*maxval;
            maxval = val;
            maxwt  = wt;
            l      = n;
            }

         else
            mean += wt*val;

         }
      rval = mean/(1.- maxwt - minwt);
      z[k] = usrnul;
      z[l] = usrnul;
      }

   c[cntdx] = rval;
   old = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}


/*

*/
mmrej(iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
/*
Here the minimum, maximum rejection without scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, usrnul, *cuts;

{
int        frmcnt, k, l, n, nn, mcc, indx, cntdx, size;
int        nims;
float      minval, maxval, val, mean;
register float       rval;

static float old = 0.0;


/*  init  */
mcc    = 0;
indx   = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];
nims   = frmcnt - 2;

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn = x[cntdx];

   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      {
      mean   = 0.0;
      minval = z[indx];
      maxval = z[indx+1];

      if (minval < maxval)
 	 {
         k      = indx;
         l      = indx+1;
         }
      else
	 {
         val    = minval;
         minval = maxval;
         maxval = val;
         k      = indx+1;
         l      = indx;
         }

      for (n=indx+2; n<indx+nn; n++) 
         { 
         val = z[n];
         if (val < minval)
	    {
            mean  += minval;
            minval = val;
            k      = n;
            }

         else if (val > maxval)
 	    {
            mean  += maxval;
            maxval = val;
            l      = n;
	    }

         else
            mean += val;

         }
      rval = mean/nims;
      z[k] = usrnul;
      z[l] = usrnul;
      }

   c[cntdx] = rval;
   old = rval;
   if (cuts[0] > rval) cuts[0] = rval;	                 /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}




/*

*/
wtsigclip(iaux,faux,x,z,c,s,m,scales,zeros,wts,
          lowsigma,highsigma,usrnul,cuts,npixc,nc)
/*
Here the sigma clipping with scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, *s, *m, usrnul, *cuts;
float      lowsigma, highsigma;
float      *zeros, *wts, *scales;

{
int        frmcnt, n, nn, mcc, k, l, idx, indx, cntdx, size;
register float       rval;
double     sum;
float      low, high, minwt, maxwt;
float      val, wt;
float      sig, resid, resid2;
float      maxresid, maxresid2;

static float old = 0.0;

/*  init  */
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn > 0) 
      {
      sum    = 0.0;
      low    = z[indx]/scales[0]-zeros[0];
      high   = z[indx+1]/scales[1]-zeros[1];
      if (low < high)
         {
         minwt  = wts[0];
         maxwt  = wts[1];
         }

      else
	 {
         minwt  = low;
         low    = high;
         high   = minwt;
         minwt  = wts[1];
         maxwt  = wts[0];
         }      

      for (n=indx+2; n<indx+nn; n++) 
         {
         idx  = n-indx;
         val  = z[n]/scales[idx]-zeros[idx];
         wt   = wts[idx];
         if (val < low) 
	    {
            sum  += minwt*low;
            low   = val;
            minwt = wt;
	    } 

         else if (val > high)
	    {
            sum  += maxwt*high;
            high  = val;
            maxwt = wt;
            }

         else
            sum  +=wt*val;
         }

      m[cntdx]    = sum/(1-maxwt-minwt);
      sum        += minwt*low + maxwt*high;
      c[cntdx]    = sum;     
      }
   indx       += frmcnt;
   }

/* Now compute sigma at each point. Correct individual residuals for the
   the image scaling */

indx = 0;

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];
   
   if (nn > 0) 
      {
      val = m[cntdx];
      sig = 0.0;
      idx = n -indx;
      for (n=indx; n<indx+nn; n++)
         {
         idx   = n-indx;
         resid = (z[n]/scales[idx]-zeros[idx]-val) * wts[idx];
         sig  += resid*resid;
         }
      s[cntdx] = sqrt(sig/(nn-2));
      }
   indx    += frmcnt;
   }


/* Now reject pixels */
indx   = 0;
mcc    = 0;
for (cntdx=0; cntdx<size; cntdx++)
   {
   nn     = x[cntdx];
   
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      { 
      rval = c[cntdx];
      val  = m[cntdx];
      sig  = s[cntdx];
      low  = -lowsigma*sig;
      high = highsigma*sig;
      
      maxresid  = (z[indx]/scales[0]-zeros[0]-val)*wts[0];
      maxresid2 = maxresid*maxresid;
      k         = indx;
      l         = 0;

      for (n=indx+1; n<indx+nn; n++) 
         {
         idx    = n-indx;
         resid  = (z[n]/scales[idx]-zeros[idx]-val)*wts[idx];
         resid2 = resid*resid;
         if (resid2 > maxresid2)
            {
            maxresid = resid;
            maxresid2 = resid2;
            k = n;
            l = n-indx;
            }
         }
    
      if ((maxresid > high) || (maxresid < low))
	 {
         rval = (rval - (z[k]/scales[l]-zeros[l])*wts[l])/(1-wts[l]);
         z[k] = usrnul;
         }
      }

   c[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;	         /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}




/*

*/
sigclip(iaux,faux,x,z,c,s,m,lowsigma,highsigma,usrnul,cuts,npixc,nc)
/*
Here the sigma clipping without scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, *s, *m, usrnul, *cuts;
float      lowsigma, highsigma;

{
int        frmcnt, n, nn, k, mcc, indx, cntdx, size;
register float       rval;
double     sum;
float      resid, resid2;
float      maxresid, maxresid2;
float      low, high;
float      val, sig;
static float old = 0.0;

/*  init  */
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn > 0)
      {
      sum    = 0.0;
      low    = z[indx];
      high   = z[indx+1];
      if (low > high)
	 {
         val    = low;
         low    = high;
         high   = val;
         }      

      for (n=indx+2; n<indx+nn; n++) 
         {
         val  = z[n];

         if (val < low) 
	    {
            sum  += low;
            low   = val;
	    } 

         else if (val > high)
	    {
            sum  += high;
            high  = val;
            }
         else
            sum  += val;
         }

      m[cntdx]    = sum/(nn-2);
      sum        += low + high;
      c[cntdx]    = sum/nn;     
      }
   indx       += frmcnt;
   }

/* Now compute sigm at each point. Correct individual residuals for the
   the image scaling */

indx   = 0;

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];
   
   if (nn > 0) 
      {
      val = m[cntdx];
      sig = 0.0;
      for (n=indx; n<indx+nn; n++)
         {
         resid = z[n] - val;
         sig  += resid*resid;
         }
      s[cntdx] = sqrt(sig/(nn-1));
      }
   indx    += frmcnt;
   }


/* Now reject pixels */
indx   = 0;
mcc    = 0;
for (cntdx=0; cntdx<size; cntdx++)
   {
   nn     = x[cntdx];
   
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      { 
      rval = c[cntdx];
      val  = m[cntdx];
      sig  = s[cntdx];
      low  = -lowsigma*sig;
      high = highsigma*sig;
      
      maxresid  = z[indx] - val;
      maxresid2 = maxresid*resid;
      k         = indx;

      for (n=indx+1; n<indx+nn; n++) 
         {
         resid  = z[n] - val;
         resid2 = resid*resid;
         if (resid2 > maxresid2)
            {
            maxresid = resid;
            maxresid2 = resid2;
            k = n;
            }
         }
    
      if ((maxresid > high) || (maxresid < low))
	 {
         rval = (nn*rval - z[k])/(nn - 1);
         z[k] = usrnul;
         }
      }

   c[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;	         /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}


/*

*/
wtavsigclip(iaux,faux,x,z,c,s,m,scales,zeros,wts,
            lowsigma,highsigma,usrnul,cuts,npixc,nc)
/*
Here the average sigma clipping with scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, *s, *m, usrnul, *cuts;
float      *zeros, *wts, *scales;
float      lowsigma, highsigma;

{
int        frmcnt, n, nn, mcc, k, l, idx, indx, cntdx, size;
register float       rval;
double     sum;
float      low, high, minwt, maxwt;
float      val, wt;
float      sigall, sig, resid, resid2;
float      maxresid, maxresid2;

static float old = 0.0;

/*  init  */
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn > 0)
      {
      sum    = 0.0;
      low    = z[indx]/scales[0]-zeros[0];
      high   = z[indx+1]/scales[1]-zeros[1];
      if (low < high)
         {
         minwt  = wts[0];
         maxwt  = wts[1];
         }

      else
	 {
         minwt  = low;
         low    = high;
         high   = minwt;
         minwt  = wts[1];
         maxwt  = wts[0];
         }      

      for (n=indx+2; n<indx+nn; n++) 
         {
         idx  = n-indx;
         val  = z[n]/scales[idx]-zeros[idx];
         wt   = wts[idx];
         if (val < low) 
	    {
            sum  += minwt*low;
            low   = val;
            minwt = wt;
	    } 

         else if (val > high)
	    {
            sum  += maxwt*high;
            high  = val;
            maxwt = wt;
            }

         else
            sum  +=wt*val;
         }

      m[cntdx]    = sum/(1-maxwt-minwt);
      sum        += minwt*low + maxwt*high;
      c[cntdx]    = sum;     
      }
   indx       += frmcnt;
  }

/* Now compute sigm at each point. Correct individual residulas for the
   the image scaling */

indx = 0;
sigall = 0.0;

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];
   
   if (nn > 0) 
      {
      val = m[cntdx];
      sig = 0.0;
      idx = n -indx;
      for (n=indx; n<indx+nn; n++)
         {
         idx   = n-indx;
         resid = (z[n]/scales[idx]-zeros[idx]-val) * wts[idx];
         sig  += resid*resid;
         }

      if (val > 0.0) 
         s[cntdx] = sqrt(val);
      else
         s[cntdx] = 1.;

      sigall += sqrt(sig)/s[cntdx];
      }
   indx    += frmcnt;
   }

sigall = sigall/sqrt((float)(nn-1))/size;
for (cntdx=0; cntdx<size; cntdx++)
    s[cntdx] = sigall * s[cntdx];


/* Now reject pixels */
indx   = 0;
mcc    = 0;
for (cntdx=0; cntdx<size; cntdx++)
   {
   nn     = x[cntdx];
   
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      { 
      rval = c[cntdx];
      val  = m[cntdx];
      sig  = s[cntdx];
      low  = -lowsigma*sig;
      high = highsigma*sig;
      
      maxresid  = (z[indx]/scales[0]-zeros[0]-val)*wts[0];
      maxresid2 = maxresid*maxresid;
      k         = indx;
      l         = 0;

      for (n=indx+1; n<indx+nn; n++) 
         {
         idx    = n-indx;
         resid  = (z[n]/scales[idx]-zeros[idx]-val)*wts[idx];
         resid2 = resid*resid;
         if (resid2 > maxresid2)
            {
            maxresid = resid;
            maxresid2 = resid2;
            k = n;
            l = n-indx;
            }
         }
    
      if ((maxresid > high) || (maxresid < low))
	 {
         rval = (rval - (z[k]/scales[l]-zeros[l])*wts[l])/(1-wts[l]);
         z[k] = usrnul;
         }
      }

   c[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;	         /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}




/*

*/
avsigclip(iaux,faux,x,z,c,s,m,lowsigma,highsigma,usrnul,cuts,npixc,nc)
/*
Here the average sigma clipping without scaling
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, *s, *m, usrnul, *cuts;
float      lowsigma, highsigma;

{
int        frmcnt, n, nn, mcc, k, indx, cntdx, size;
register float       rval;
double     sum;
float      low, high;
float      val;
float      sigall, sig, resid, resid2;
float      maxresid, maxresid2;

static float old = 0.0;

/*  init  */
indx = 0;
frmcnt = iaux[6];
size = npixc[0] * npixc[1];

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn > 0)
      {
      sum    = 0.0;
      low    = z[indx];
      high   = z[indx+1];
      if (low > high)
	 {
         val    = low;
         low    = high;
         high   = val;
         }      

      for (n=indx+2; n<indx+nn; n++) 
         {
         val  = z[n];
         if (val < low) 
	    {
            sum  += low;
            low   = val;
	    } 

         else if (val > high)
	    {
            sum  += high;
            high  = val;
            }

         else
            sum  += val;
         }

      m[cntdx]    = sum/(nn-2);
      sum        += low + high;
      c[cntdx]    = sum/nn;     
      }
   indx       += frmcnt;
 }

/* Now compute sigm at each point. Correct individual residulas for the
   the image scaling */

indx = 0;
sigall = 0.0;

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];
   
   if (nn > 0) 
      {
      val = m[cntdx];
      sig = 0.0;
      for (n=indx; n<indx+nn; n++)
         {
         resid = z[n]-val;
         sig  += resid*resid;
         }

      if (val > 0.0) 
         s[cntdx] = sqrt(val);
      else
         s[cntdx] = 1.;

      sigall += sqrt(sig)/s[cntdx];
      }
   indx    += frmcnt;
 }

sigall = sigall/sqrt((float)(nn-1))/size;
for (cntdx=0; cntdx<size; cntdx++)
    s[cntdx] = sigall * s[cntdx];


/* Now reject pixels */
indx   = 0;
mcc    = 0;
for (cntdx=0; cntdx<size; cntdx++)
   {
   nn     = x[cntdx];
   
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      mcc ++;		                            /* increment null count  */
      } 

   else
      { 
      rval = c[cntdx];
      val  = m[cntdx];
      sig  = s[cntdx];
      low  = -lowsigma*sig;
      high = highsigma*sig;
      
      maxresid  = z[indx]-val;
      maxresid2 = maxresid*maxresid;
      k         = indx;

      for (n=indx+1; n<indx+nn; n++) 
         {
         resid  = z[n] - val;
         resid2 = resid*resid;
         if (resid2 > maxresid2)
            {
            maxresid = resid;
            maxresid2 = resid2;
            k = n;
            }
         }

      if ((maxresid > high) || (maxresid < low))
	 {
         rval = (nn*rval - z[k])/(nn-1);
         z[k] = usrnul;
         }
      }

   c[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;	         /* update min + max */
   if (cuts[1] < rval) cuts[1] = rval;
   indx += frmcnt;
   }

*nc = mcc;
return;
}

/*

*/
wtsigma(iaux,faux,x,z,c,s,scales,zeros,wts,usrnul,cuts,npixc,nc)
/*
Compute the sigma frame with weigthing
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, *s, usrnul, *cuts;
float      *scales,*zeros,*wts;

{
int        frmcnt, n, nn, indx, idx, cntdx, size;
register float       rval;
float      sig, val, pixval, sumwts;
int        np;

static float old = 0.0;


/*  init  */
indx   = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];
   if (nn == 0) 
      {
      if (iaux[8] == 1)
         s[cntdx] = old;	                           /* use last value */
      else
         s[cntdx] = usrnul;
      }
   else
      {
      val    = c[cntdx];
      sig    = 0.0;
      sumwts = 0.0;
      np     = 0;

      for (n=indx; n<indx+nn; n++) 
         {
         idx    = n -indx;
         pixval = z[n];
         if (pixval != usrnul)
	    {
            np++;
            sig     =  wts[idx] *
                      ((pixval/scales[idx]-zeros[idx]-val) * 
                       (pixval/scales[idx]-zeros[idx]-val));
            sumwts += wts[idx];
            }
         }
 
      if (np > 0)
         rval = sqrt(sig/sumwts * np/(np-1));
      else
         rval = 0.0;
      }

   s[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;
   if (cuts[1] < rval) cuts[1] = rval;
   if (cuts[2] > nn) cuts[2] = nn;
   if (cuts[3] < nn) cuts[3] = nn;
   indx += frmcnt;
   }

return;
}


/*

*/
sigma(iaux,faux,x,z,c,s,usrnul,cuts,npixc,nc)
/*
Computer the sigma frame without weighting
*/

int        *npixc;
int        *iaux, *nc;
short int  *x;
float      *faux, *z, *c, *s, usrnul, *cuts;

{
int        frmcnt, n, nn, indx, cntdx, size;
register float       rval;
float      sig, val, pixval;
int        np;

static float old = 0.0;


/*  init  */
indx = 0;
frmcnt = iaux[6];
size   = npixc[0] * npixc[1];

for (cntdx=0; cntdx<size; cntdx++)		/* loop through count buffer */
   {
   nn     = x[cntdx];

   if (nn == 0) 
      {
      if (iaux[8] == 1)
         rval = old;		                           /* use last value */
      else
         rval = usrnul;
      }
   else
      {
      val  = c[cntdx];
      sig  = 0.0;
      np   = 0;

      for (n=indx; n<indx+nn; n++) 
         {
         pixval  = z[n];
         if (pixval != usrnul)
	    {
            sig += (pixval-val)*(pixval-val);
            np++;
            } 
         }

      if (np > 0)
         rval = sqrt(sig/np);
      else
         rval = 0.0;
      }

   s[cntdx] = rval;
   old      = rval;
   if (cuts[0] > rval) cuts[0] = rval;
   if (cuts[1] < rval) cuts[1] = rval;
   if (cuts[2] > np) cuts[2] = nn;
   if (cuts[3] < np) cuts[3] = nn;
   indx += frmcnt;
   }

return;
}