/* @(#)averag.c	17.1.1.1 (ESO-DMD) 01/25/02 17:40:55 */
/*===========================================================================
  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 AVERAG			version 1.00	890120
  K. Banse                  		ESO - Garching
  1.10	   890918	1.20	900719		1.30	920313
 
.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
 000425

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

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

#define  MAXIMS  64
#define  MAXIMSONE  MAXIMS+1


/*

*/

void sortit(rka,ndim)

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

{
register int  m, j, i;

float  zka;
  

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

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

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

/*

*/

void 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    nn, nin, nini, nout, jump, indx, cntdx;
int    nx, ny, xfirst, yfirst, xend, yend;
register int nr;

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 (nr=0; nr<nn; nr++)
      x[nr] = 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 (nr=0; nr<npixc[0]; nr++)
            {
            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 (nr=0; nr<npixc[0]; nr++)
            {
            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 (nr=0; nr<npixc[0]; nr++)
               {
               if ((nr >= xfirst) && (nr <= 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 (nr=0; nr<npixc[0]; nr++)
               {
               if ((nr >= xfirst) && (nr <= 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];
            }
         }
      }
   }
}

/*

*/

void add(flag,iaux,faux,x,z,c,usrnul,cuts,npixc,nc)
int   *npixc;
int   flag, *iaux, *nc;
short int   *x;
float       *faux, *z, *c, usrnul, *cuts;

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



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

/*   branch according to FLAG */

if (flag == 4) goto median;			/* that is the `worst' part */


if (flag == 1)
  
   /* -------------------------------*/
   /*   here for the normal average  */
   /* -------------------------------*/
   
   {
   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) 
            {
            sum = z[indx];
            for (n=indx+1; n<indx+nn; n++)
                sum += z[n];
            rval = sum / nn;
            }
         else
            rval = z[indx];
         }

      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 if (flag == 2)
  
   /* -----------------------*/
   /*  here for the minimum  */
   /* -----------------------*/
 
   {
   if (iaux[3] > 0)			/* take average of MIN and the next */
      {					/* `iaux[3]' higher values */
      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)
               {
               k = 1;                              /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
                  rbuf[k++] = z[n];
               sortit(rbuf,nn);                    /* sort array */
               k = iaux[3] + 1;
               if (nn > k)
                  nh = k;
               else
                  nh = nn;
               sum = rbuf[1];		  /* avrage over max. iaux[3] values */
               for (k=2; k<=nh; k++)
                   sum += rbuf[k];
               rval = sum / nh;
               }
            else
               rval = z[indx];
            }

         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
      {
      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
            {
            rval = z[indx];
            if (nn > 1)
               {
               for (n=indx+1; n<indx+nn; n++)
                  if (rval > z[n]) rval = z[n];
               }
            }

         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 if (flag == 3)

   /* -----------------------*/
   /*  here for the maximum  */
   /* -----------------------*/

   {
   if (iaux[3] > 0)                     /* take average of MAX and the next */
      {                                 /* `iaux[3]' lower values */
      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)
               {
               k = 1;                              /* start at index 1 ... */
               for (n=indx; n<indx+nn; n++)
                  rbuf[k++] = z[n];
               sortit(rbuf,nn);                    /* sort array */
               k = iaux[3] + 1;
               if (nn > k)
                  nh = k;
               else
                  nh = nn;
               sum = rbuf[nn];            /* avrage over max. iaux[3] values */
               for (k=nn-1; k>nn-nh; k--)
                   sum += rbuf[k];
               rval = sum / nh;
               }
            else
               rval = z[indx];
            }

         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
      {
      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
            {
            rval = z[indx];
            if (nn > 1)
               {
               for (n=indx+1; n<indx+nn; n++)
                  if (rval < z[n]) rval = z[n];
               }
            }

         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:								/*
-------

............................................................


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];
   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];
               sortit(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
      {  
      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];
               sortit(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
   {			
   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];
               sortit(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)
               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];
               sortit(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;
}

/*

*/

main()

{
int    uni, ec, ed, el, ln0, ln1;
int    imnoc, imnox, 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    stat, nolin, begin;
int    stripe, optio;
int    frmcnt, nulcnt, kk, m, dim3, floff;
int    debufl=0, planesize=0;
short int   *xpntr;
register int nr;

char     *wpntr, *tmpntr;
char     line[84], action[4];
char     frame[84], framec[84], catfil[84];
char     cunit[64], ident[72], output[84];
static char  error1[] = "operands do not match in stepsize... ";
static char  error2[] = "operands do not overlap... ";
static char  error3[] = "stepsizes of different sign... ";
static char  error4[] = "catalog empty... ";
static char  mesalloc[] = "could not allocate virtual memory...";

double   step[MAXIMS][3], start[MAXIMS][3];
double   stepc[3], ostep[3], endc[3];
double   startc[3], ostart[3], oldend[3];
double   aostep;
 
float    *pntr[MAXIMS], *pntrc, *zpntr;
float    dif, eps[3], cuts[4], w[MAXIMS];
float    usrnul;
float    faux[5];
 
 
/* set up MIDAS environment + enable automatic error abort */

SCSPRO("averag");

for (nr=0; nr<3; nr++)
   {
   apix[nr][0] = 0;
   apix[nr][1] = 0;
   cpix[nr] = 0;
   startc[nr] = 0.0;
   stepc[nr] = 1.0;
   endc[nr] = 0.0;
   npixc[nr] = 1;
   zpix[nr] = 1;
   ostart[nr] = 0.0;
   oldend[nr] = 0.0;
   ostep[nr] = 1.0;
   npix[0][nr] = 1;
   start[0][nr] = 0.0;
   step[0][nr] = 1.0;
   }
 
 
/* get result frame, input frame list, action flag */

(void) SCKGETC("OUT_A",1,80,&iav,framec);
(void) SCKGETC("P3",1,80,&iav,line);
(void) SCKGETC("ACTION",1,2,&iav,action);
CGN_UPSTR(action);			/* convert to upper case */
nulcnt = 0;
(void) SCKGETC("MID$SPEC",1,5,&iav,output);
output[5] = '\0';
if ( strcmp(output,"DEBUG") == 0) debufl = 1;		/* set debug flag */


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

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] = INDEX/DATA flag - not used here
iaux[2] = 0, take all data
        = 1, take only data in interval [ faux[2],faux[3] ]
iaux[3,4] = average limits - not used here
iaux[5] = 0, `nomerge' option
        = 1, `merge' option
iaux[6] = frame count
iaux[7] = no. of current frame
iaux[8] = 0, use NULL value
        = 1, use last pixel as NULL value

faux[0,1] = average limits if iaux[1] == 1 - not used here
faux[2,3] = valid_data_interval if iaux[2] = 1
faux[4] = weighting factor

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

for (nr=0; nr<10; nr++) iaux[nr] = 0;

 
/* get Null value, average option+limits, valid_data_interval */

(void) SCKGETC("P5",1,40,&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,"invalid `null value' ...");
   }

optio = 1;
(void) SCKGETC("P6",1,60,&iav,output);
CGN_UPCOPY(frame,output,2);		/* upper case -> frame */

if (frame[0] == 'M')
   {
   if (frame[1] == 'I')
      optio = 2;
   else if (frame[1] == 'A')
      optio = 3;
   else if (frame[1] == 'E')
      optio = 4;
   }

if (optio != 1)
   {	
   begin = CGN_INDEXC(output,',')+1;
   if ((begin > 1) && (output[begin] != '\0'))
      {
      (void) strcpy(output,&output[begin]);
      if (optio < 4)			/* must be MIN,n or MAX,n */
         {
         kk = CGN_CNVT(output,1,1,npixc,&dif,&aostep);
         if (kk <= 0)
            SCETER(8,"invalid syntax in option_string ...");
         iaux[3] =  npixc[0];
         }
      else
         {						/* MED,r,q,DATA */
         iav = CGN_CNVT(output,2,2,npixc,faux,&aostep);
         if (iav < 1)
            SCETER(8,"invalid syntax in `option_string' ...");
         if (iav == 1) faux[1] = faux[0];
         iaux[3] = faux[0];
         iaux[4] = faux[1];
         					/* skip max. two numbers */
         begin = 0;
         m = (int) strlen(output);			/*          or: MAX,n */
         kk = CGN_EXTRSS(output,m,',',&begin,cunit,40);
         if (iav == 2)
            kk = CGN_EXTRSS(output,m,',',&begin,cunit,40);

         if ((output[begin] == 'D') || (output[begin] == 'd'))
            iaux[1] = 1;		   /* show that it's DATA method */
         }
      if ((iaux[3] < 0) || (iaux[4] < 0))
         SCETER(9,"negative `average limits' ...");
      }
   }
 

/*  now let's look if we have a valid data interval */

(void) SCKGETC("P7",1,60,&iav,output);
if (output[0] != '+')
   {
   kk = CGN_CNVT(output,2,2,npixc,&faux[2],&aostep);
   if (kk < 2)
      SCETER(10,"invalid syntax in `valid_data_interval' ...");
   if (faux[2] > faux[3])
      SCETER(11,"invalid `valid_data_interval' ...");
   iaux[2] = 1;			/* indicate we have a valid_data_interval */
   }


/*  test, if we have list of frames or catalog */

kk = CGN_INDEXS(line,".cat");
if (kk <= 0) kk = CGN_INDEXS(line,".CAT");
 

/*   here we handle input from a catalog - get name of first image file 
     in catalog                                                         */
 
if (kk > 0) 
   {
   if ((int)strlen(line) > 63)
      SCETER(3,"catalog name too long...");
   else
      (void) strcpy(catfil,line);

   iseq = 0;
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)	    /* max. MAXIMS frames */
      {
      (void) SCCGET(catfil,0,frame,output,&iseq);
      if (frame[0] == ' ') break;		/*  indicates the end ... */
    
      imnol[frmcnt] = 0;
      (void) SCFOPN(frame,D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
      }
   }

 
/*  here we handle input from single input line - pull out operands */

else
   {
   begin = 0;
   m = (int) strlen(line);
   for (frmcnt=0; frmcnt<MAXIMS; frmcnt++)
      {
      kk = CGN_EXTRSS(line,m,',',&begin,output,60);
      if (kk <= 0) break;
       
      CGN_FRAME(output,F_IMA_TYPE,frame,0);	/* convert frames  */
      imnol[frmcnt] = 0;
      stat = SCFOPN(frame,D_R4_FORMAT,0,F_IMA_TYPE,&imnol[frmcnt]);
      }
   }

if (frmcnt <= 0) SCETER(4,error4);       /* there must be at least 1 image  */

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


for (nr=0; nr<3; nr++)
   {
   ostart[nr] = start[0][nr];
   ostep[nr] = step[0][nr];
   aostep = ostep[nr];
   if (aostep < 0.0) aostep = - aostep;
   eps[nr] = 0.0001 * aostep;			/* take 0.01% of stepsize */
   oldend[nr] = ostart[nr] + (npix[0][nr]-1)*ostep[nr];
   }
cuts[0] = 999999.0;					/* cuts[0] > cuts[1] */
cuts[1] = -999999.0;
 

/*   if weighting option, loop + look for descriptor WEIGHT  */

for (nr=0; nr< frmcnt; nr++)		/* initialize weights to 1.0 */
   w[nr] = 1.0;

if (action[1] == 'W') 
   {
   stat = SCECNT("GET",&ec,&el,&ed);
   ln1 = 1; ln0 = 0;
   stat = SCECNT("PUT",&ln1,&ln0,&ln0);

   for (nr=0; nr< frmcnt; nr++)
      {
      stat = SCDRDR(imnol[nr],"WEIGHT",1,1,&iav,&w[nr],&uni,&nulo);
      if (stat != ERR_NORMAL) 
         {
         (void) sprintf(line,"frame no. %ld: descr. WEIGHT missing ",imnol[nr]);
         SCTPUT(line);
         SCTPUT("weight defaulted to 1.0 ...");
         }
      }

   stat = SCECNT("PUT",&ec,&el,&ed);
   }
 

/*   now loop through the other input frames  */
  
for (nr=1; nr< frmcnt; nr++)
   {
   for (m=0; m<3; m++)				/* init values... */
      {
      start[nr][m] = 0.0;
      step[nr][m] = 1.0;
      npix[nr][m] = 1;
      }
  
   (void) SCDRDI(imnol[nr],"NAXIS",1,1,&iav,&naxis[nr],&uni,&nulo);
   (void) SCDRDI(imnol[nr],"NPIX",1,3,&iav,&npix[nr][0],&uni,&nulo);
   (void) SCDRDD(imnol[nr],"START",1,3,&iav,&start[nr][0],&uni,&nulo);
   (void) SCDRDD(imnol[nr],"STEP",1,3,&iav,&step[nr][0],&uni,&nulo);
 

/*  stepsizes should have same sign and not differ too much...  */

   for (m=0; m<3; m++)
      {
      if ((ostep[m]*step[nr][m]) <= 0.) SCETER(1,error3);
      aostep = step[nr][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[nr][m] + (npix[nr][m]-1)*step[nr][m];
         if (ostep[m] > 0.) 
            {
            if (ostart[m] < start[nr][m]) ostart[m] = start[nr][m];   /* MAX */
            if (oldend[m] > dif) oldend[m] = dif;		      /* MIN */
            }
         else
            {
            if (ostart[m] > start[nr][m]) ostart[m] = start[nr][m];   /* MIN */
            if (oldend[m] < dif) oldend[m] = dif;                     /* MAX */
            }
         }
      }
   else
      {
      for (m=0; m<3; m++)
         {
         dif = start[nr][m] + (npix[nr][m]-1)*step[nr][m];
         if (ostep[m] < 0.) 
            {
            if (ostart[m] < start[nr][m]) ostart[m] = start[nr][m];   /* MAX */
            if (oldend[m] > dif) oldend[m] = dif;		      /* MIN */
            }
         else
            {
            if (ostart[m] > start[nr][m]) ostart[m] = start[nr][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);
   }

 
/*   create new result frame with dimension as intersection 
     or union of input frames                                   */

naxisc = naxis[0];
if (action[0] == 'M') 
   {
   for (nr=1; nr<frmcnt; nr++)
      {
      if (naxisc < naxis[nr]) naxisc = naxis[nr];		/* MAX */
      }
   }
else
   {
   for (nr=1; nr<frmcnt; nr++)
      {
      if (naxisc > naxis[nr]) naxisc = naxis[nr];		/* MIN */
      }
   }
 

/*  check, if input from a cube */

if (naxisc > 2) 
   {
   dim3 = npix[0][2];
   naxisc = 2;
   }
else
   dim3 = 0;


/*  set up standard stuff for result frame        */

sizec = 1;
for (nr=0; nr<naxisc; nr++)
   {
   startc[nr] = ostart[nr];
   stepc[nr] = ostep[nr];
   npixc[nr] = CGN_NINT((oldend[nr]-ostart[nr]) / stepc[nr]) + 1;
   sizec = sizec * npixc[nr];
   }

imnoc = 0;
(void) SCFCRE(framec,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnoc);
(void) SCDWRI(imnoc,"NAXIS",&naxisc,1,1,&uni);
(void) SCDWRI(imnoc,"NPIX",npixc,1,naxisc,&uni);
(void) SCDWRD(imnoc,"START",startc,1,naxisc,&uni);
(void) SCDWRD(imnoc,"STEP",stepc,1,naxisc,&uni);
(void) strcpy(ident,"average frame ");
(void) SCDWRC(imnoc,"IDENT",1,ident,1,72,&uni);
ln1 = (naxisc+1) * 16;
(void) SCDWRC(imnoc,"CUNIT",1,cunit,1,ln1,&uni);

  
/* in case of debugging save the counts in MIDAS image averdumy.dum  */

if (debufl == 1)
   {
   imnox = 0;
   (void) strcpy(frame,"averdumy.dum");
   (void) SCFCRE(frame,D_I2_FORMAT,F_O_MODE,F_IMA_TYPE,sizec,&imnox);
   (void) SCDWRI(imnox,"NAXIS",&naxisc,1,1,&uni);
   (void) SCDWRI(imnox,"NPIX",npixc,1,naxisc,&uni);
   (void) SCDWRD(imnox,"START",startc,1,naxisc,&uni);
   (void) SCDWRD(imnox,"STEP",stepc,1,naxisc,&uni);
   (void) strcpy(ident,"valid pixel counts ");
   (void) SCDWRC(imnox,"IDENT",1,ident,1,72,&uni);
   (void) SCDWRC(imnox,"CUNIT",1,cunit,1,ln1,&uni);
   }


/* see, if keyword MONITPAR(20) holds amount of virtual memory available */

(void) SCKRDI("MONITPAR",20,1,&iav,&chunkc,&uni,&nulo);
ln1 = (chunkc * chunkc) / 5;   		 /* `chunkc' = x-dim of square image
					    which can be fully 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 */
kk = chunkc * sizeof(float);
wpntr = malloc((unsigned int)kk);
if (wpntr == (char *) 0)
   SCETER(66,mesalloc);
else
   pntrc = (float *) wpntr;

if (dim3 > 1)
   {
   if (npix[0][1] < stripe) 
      iav = npix[0][1];				/* iav = y-dim  */
   else
      iav = stripe;

   chunk[0] = iav * npix[0][0];		/* always stripe lines  */
   for (nr=0; nr<dim3; nr++)
      {					/* all chunks like chunk[0] */
      chunk[nr] = chunk[0];
      kk = chunk[nr] * sizeof(float);
      wpntr = malloc((unsigned int)kk);
      if (wpntr == (char *) 0)
         SCETER(66,mesalloc);
      else
         pntr[nr] = (float *) wpntr;
      w[nr] = w[0];
      imnol[nr] = imnol[0];
      for (m=0; m<naxisc; m++)
         {
         npix[nr][m] = npix[0][m];
         start[nr][m] = start[0][m];
         step[nr][m] = step[0][m];
         }
      }
   }
else
   {
   for (nr=0; nr<frmcnt; nr++)
      {
      if (npix[nr][1] < stripe) 
         iav = npix[nr][1];				/* iav = y-dim  */
      else
         iav = stripe;

      chunk[nr] = iav * npix[nr][0];		/* always stripe lines */
      kk = chunk[nr] * sizeof(float);
      wpntr = malloc((unsigned int)kk);
      if (wpntr == (char *) 0)
         SCETER(66,mesalloc);
      else
         pntr[nr] = (float *) wpntr;
      }
   }
 

/*  now map chunk for z-direction and count buffer  */

if (dim3 > 1)
   {
   planesize = npixc[0] * npixc[1];
   frmcnt = dim3;
   }

iaux[6] = frmcnt;
m = frmcnt * chunkc;
kk = m * sizeof(float);
wpntr = malloc((unsigned int)kk);
if (wpntr == (char *) 0)
   SCETER(66,mesalloc);
else
   zpntr = (float *) wpntr;

kk = chunkc * sizeof(short int);
wpntr = malloc((unsigned int)kk);
if (wpntr == (char *) 0)
   SCETER(66,mesalloc);
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];
  

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 (nr=0; nr<frmcnt; nr++)
         chunk[nr] = stripe * npix[nr][0];
      }
   zpix[1] = stripe;
 
   floff = 0;				/* for data cube input */
   for (nr=0; nr<frmcnt; nr++)
      {
 
      /*  convert start + end of overlap region into pixel no.'s   */

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

sect_5360:						/* fill z-buffer */
      iaux[7] = nr;
      faux[4] = w[nr];
      fill(iaux,faux,pntr[nr],xpntr,zpntr,apix,cpix,npix[nr][0],zpix);
      floff += planesize;
      }
  

   /*   now do the calculus  */

   add(optio,iaux,faux,xpntr,zpntr,pntrc,usrnul,cuts,zpix,&kk);

  
   /*  and write results to disk */

   felm = nolin*npixc[0] + 1;
   tmpntr = (char *) pntrc;
   (void) SCFPUT(imnoc,felm,chunkc,tmpntr);
   nulcnt += kk;				/* update null count */
   if (debufl == 1)
      {
      tmpntr = (char *) xpntr;
      (void) SCFPUT(imnox,felm,chunkc,tmpntr);	/* if debug, update on disk */
      }
  

   /*  follow chunks with start + end value  */

   aostep = stripe*stepc[1];
   startc[1] += aostep;
   endc[1] += aostep;
   }

  
/*  update descriptors + cuts of result frame  */

frame[0] = ' ';
frame[1] = '\0';
CGN_DSCUPD(imnol[0],imnoc,frame);

cuts[2] = cuts[0];
cuts[3] = cuts[1];
(void) SCDWRR(imnoc,"LHCUTS",cuts,1,4,&uni);
 
cuts[0] = nulcnt;			/* update key NULL */
if (iaux[8] == 0)
   {
   cuts[1] = usrnul;
   ln1 = 2;
   }
else
   ln1 = 1;
(void) SCKWRR("NULL",cuts,1,ln1,&uni);
if (nulcnt > 0) 
   {
   if (iaux[8] == 0)
      (void) sprintf(output,
      "%d undefined pixels, set to `null value' (= %12.6f)",nulcnt,usrnul);
   else
      (void) sprintf(output,
      "%d undefined pixels, set to `previous pixel'",nulcnt);
   SCTPUT(output);
   }
  
SCSEPI();
}