#include "midas_def.h"

int exmed(actio,ino,zbins,npix,spix,epix,valmed,stat)
int      actio, ino, *npix, *spix, *epix, stat;
float    *zbins, *valmed;
{
  int      iav,n,chunka,chunkb,nn;
  int      na,nb,insiz,outsiz;
  int      uni,ibuf[8];
  int      imnox,imnoy;
  char     *pntrx, *pntry;
  float    rbuf[8], *fptrx, *fptry;

  if (actio == 5)                           /* plane option */
    {
      insiz = npix[1] * npix[2];
      outsiz = insiz;
      na = (spix[2]*insiz) + 1;
      stat = SCFCRE("middstat2",D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,outsiz,&imnoy);
      stat = SCFMAP(imnoy,F_X_MODE,1,outsiz,&iav,&pntry);
      stat = SCFGET(ino,na,insiz,&iav,pntry);
      fptry = (float*) pntry;
    }
  else if (actio != 0)                      /* subframe input */
    {
      nn = epix[1]-spix[1]+1;
      insiz = npix[1] * nn;                 /* size of full y-strip */
      chunka = npix[1];
      chunkb = epix[0] - spix[0] + 1;
      outsiz = chunkb * nn;
      na = spix[1]*npix[1] + 1;

      stat = SCFCRE("middstat1",D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,chunka,&imnox);
      stat = SCFMAP(imnox,F_X_MODE,1,chunka,&iav,&pntrx);
      stat = SCFCRE("middstat2",D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,outsiz,&imnoy);
      stat = SCFMAP(imnoy,F_X_MODE,1,outsiz,&iav,&pntry);
      fptrx = (float*) pntrx;
      fptry = (float*) pntry;
      nb = 1;
      for (n=1;n<=nn;n++)
	{
	  stat = SCFGET(ino,na,chunka,&iav,pntrx);
	  datmov(fptrx,nb-1,chunkb,spix[0],fptry);
	  na = na + chunka;
	  nb = nb + chunkb;
	}
      stat = SCFCLO(imnox);
    }
  else                                         /* full frame */
    {
      insiz = npix[1] * npix[2] * npix[3];
      outsiz = insiz;
      stat = SCFCRE("middstat2",D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,outsiz,&imnoy);
      stat = SCFMAP(imnoy,F_X_MODE,1,outsiz,&iav,&pntry);
      stat = SCFGET(ino,1,insiz,&iav,pntry);
      fptry = (float*) pntry;
    }

/* sort array + return median */
  na = (outsiz+1)/2;                             /* index of median */
  fptry = fptry - 1;

  stat = sortmed(fptry,zbins,outsiz,na,valmed);

  if (stat == -1)
    {
      SCTPUT("no pixels found with data in given interval...");
      for (n=1;n<8;n++)
	{
	  rbuf[n] = -1.0;
	  ibuf[n] = -1;
	}
      stat = SCKWRR("OUTPUTR",rbuf,1,7,&uni);
      stat = SCKWRI("OUTPUTI",ibuf,1,7,&uni);
      SCFCLO(imnoy);                       /* keep value of stat */
    }
  else
    {
      stat = SCKWRR("OUTPUTR",valmed,8,1,&uni);
      stat = SCFCLO(imnoy);
    }

  return(0);
}



int datmov(a,off,chunk,spix,y)
int      off,chunk,spix;
float    *a, *y;
{
  int   n, nn;

  for (n=spix,nn=off;n<spix+chunk;n++,nn++) y[nn] = a[n];
 
  return(0);
}