/* @(#)fitsrdm.c	17.1.1.1 (ES0-DMD) 01/25/02 17:39:09 */
/*===========================================================================
  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 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT  (c)  1994  European Southern Observatory
.IDENT      fitsrdm.c
.LAUGUAGE   C
.AUTHOR     P.Grosbol   ESO/IPG
.KEYWORDS   FITS data matrix, decode, read
.COMMENT    read the prime data matrix of FITS file
.VERSION    1.0  1988-Dec-10 : Creation,   PJG 
.VERSION    1.1  1989-Oct-08 : Save Group parm. in table,   PJG 
.VERSION    1.2  1989-Oct-24 : Check of data to store,   PJG 
.VERSION    1.3  1990-Feb-04 : Change call-seq. for cv-routine, PJG 
.VERSION    1.4  1991-Feb-14 : Correct pointer error, PJG
.VERSION    2.0  1991-Feb-24 : Change call-seq. structures, PJG 
.VERSION    2.1  1991-Mar-23 : Compute min/max if not given, PJG
.VERSION    2.2  1993-Sep-23 : Check overflow in cut values, PJG
.VERSION    2.3  1993-Dec-02 : Update to new SC + prototypes, PJG
.VERSION    2.4  1994-May-09 : Cast variable for ANSI-C, PJG
.VERSION    2.5  1994-Jun-28 : Include UI2 as bitpix=-16, PJG
.VERSION    2.6  1994-Sep-07 : Add explicit cast for unsigned, PJG
.VERSION    2.7  1997-Jan-08 : Only warning on wrong record size, PJG
---------------------------------------------------------------------*/
#include   <computer.h>            /* computer specific constants    */
#include   <fitsfmt.h>             /* general data definitions       */
#include   <fitsdef.h>             /* basic FITS definitions         */
#include   <midas_def.h>           /* MIDAS definitions              */

#define    MXFB              2880  /* max. size of scaling buffer    */

#ifdef __STDC__
int fitsrdm(int mfd , BFDEF * bfdef , int size , int mfdt , char fmt,
            int Midas_flag)
#else
int fitsrdm(mfd,bfdef,size,mfdt,fmt,Midas_flag)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE       Read prime data matrix in FITS
.RETURN        error code - 0: OK, -1: error
---------------------------------------------------------------------*/
int       mfd;                 /* IN: MIDAS file descriptor          */
BFDEF     *bfdef;              /* IN: pointer to FITS header parm's  */
int       size;                /* IN: size of data matrix (bytes)    */
int       mfdt;                /* IN: MIDAS file desc. Groups table  */
char      fmt;                 /* IN: data format Original/Float     */
int       Midas_flag;	       /* IN: = 0, `usual' mode
                                      = 1, `internal FITS access'  */

#endif
{
  register int   bf, k;
  unsigned char  uc, *pb;
  short          s, *ps;
  unsigned short u, *pu;
  int            fpf, dno, gno, pno, npix, n, fmm;
  int            dn, pn, ndata, dfmt, nb, pcnt, dsf, unit[4];
  int            felem, i, *pi, mnvi, mxvi;
  float          *pf, cuts[2];
  double         fac, zero, d, *pd, mnvd, mxvd;
  PDEF           *pp;
  union { 
	  float    f[2*MXFB];
	  double     d[MXFB];
        } buf;
  union { 
	  unsigned char   *b;
	  char            *c;
	  short           *s;
	  unsigned short  *u;
	  int             *i;
	  float           *f;
	  double          *d;
        } p;


  int    option;



  mnvd = mxvd = 0.0;

  pno = bfdef->pcount; pcnt = 0;             /* initiate counters     */
  pp = bfdef->parm;
  gno = 0;
  bf = bfdef->bflag;
  dfmt = bfdef->bitpix;
  nb = (dfmt<0) ? -dfmt/8 : dfmt/8;
  ndata = size/(nb*bfdef->gcount) - pno; dno = ndata;
  fac = bfdef->bscale; zero = bfdef->bzero;
  felem = 1;
  fmm = ((bfdef->mflag != 3) && (0<size));
  fpf = (bfdef->sflag || dfmt!=-32) ;
  dsf = ((bfdef->sflag && dfmt!=-64) || fmt=='F') ? -32 : dfmt;

  while (0<size) {                 /* read all data in prime matrix  */
    if ((n=dread(&p.c,FITSLR)) != FITSLR) 	/* read next data record  */
       {			
       if (size<=n) 
	  SCTPUT("Warning: incomplete FITS record read!");
       else 
          {
          int   unit;
          char  tbuf[80];

	  SCTPUT("Error: unexpected EOF");
          switch (dfmt)
             {
             case   8 : break;
             case -16 :
             case  16 : size /= 2; break;
             case  32 : 
             case -32 : size /= 4; break;
             case -64 : size /= 8; break;
             }
	  (void) sprintf(tbuf,"%d data values still missing",size);
          SCTPUT(tbuf);
	  if (0<=mfd) SCFCLO(mfd); mfd = -1;
          (void) SCKWRI("OUTPUTI",&size,16,1,&unit);	/* save value */
	  return NOFITS;
          }
       }
    if (size<n) n = size;
    size -= n;                     /* decrement remaining bytes      */
    switch (dfmt) {                /* convert to local data format   */
       case   8 : npix = n; break;
       case -16 :
       case  16 : npix = n/2; cvi2(p.s,npix,0); break;
       case  32 : npix = n/4; cvi4(p.i,npix,0); break;
       case -32 : npix = n/4; cvr4(p.f,npix,0); break;
       case -64 : npix = n/8; cvr8(p.d,npix,0); break;
    }
    do {                           /* scale all values if needed     */
       if (0<pno) {                /* decode groups parameters       */
	  pn = (pno<npix) ? pno : npix;
	  pno -= pn; npix -= pn;
	  while (pn--) {
	      switch (dfmt) {
	         case   8 : d = pp->pscal * (*(p.b++)) + pp->pzero;
			    break;
                 case -16 :
	         case  16 : d = pp->pscal * (*(p.s++)) + pp->pzero;
			    break;
	         case  32 : d = pp->pscal * (*(p.i++)) + pp->pzero;
			    break;
	         case -32 : d = pp->pscal * (*(p.f++)) + pp->pzero;
			    break;
	         case -64 : d = pp->pscal * (*(p.d++)) + pp->pzero;
			    break;
              }
	      pcnt++; pp++;
              if (0<=mfdt) TCEWRD(mfdt,gno+1,pcnt,&d);
	  }
       }
       if (!pno && 0<dno && npix) {    /* decode data values            */
	  dn = (dno<npix) ? dno : npix;
	  dno -= dn; npix -= dn;
	  if (bfdef->sflag) {
	     switch (dfmt) {
		case   8 : pf = buf.f; uc = bfdef->blank;
			   for (k=0; k<dn; k++, pf++, p.b++) {
                               if (bf && *p.b==uc) toNULLFLOAT(*pf);
                                  else *pf = fac*((int)(*p.b)) + zero;
                           }
	                   break;
		case -16 :
		case  16 : pf = buf.f; s = bfdef->blank;
			   for (k=0; k<dn; k++, pf++, p.s++) {
                               if (bf && *p.s==s) toNULLFLOAT(*pf);
                                  else *pf = fac*(*p.s) + zero;
                           }
	                   break;
		case  32 : pf = buf.f; i = bfdef->blank;
			   for (k=0; k<dn; k++, pf++, p.i++) {
                               if (bf && *p.i==i) toNULLFLOAT(*pf);
                                  else *pf = fac*(*p.i) + zero;
                           }
	                   break;
		case -32 : pf = buf.f;
			   for (k=0; k<dn; k++)
                               *pf++ = fac * (*(p.f++)) + zero;
	                   break;
		case -64 : pd = buf.d;
			   for (k=0; k<dn; k++)
                               *pd++ = fac * (*(p.d++)) + zero;
	                   break;
	     }
             if (Midas_flag != 0)
                {
	        if (dfmt!=-64)             /* store data in MIDAS file */
	          MyPut(-32,felem,dn,(char *) buf.f);
	        else  MyPut(dfmt,felem,dn,(char *) buf.d);
                }
             else
                {
	        if (dfmt!=-64)             /* store data in MIDAS file */
	          SCFPUT(mfd,felem,dn,(char *) buf.f);
	        else  SCFPUT(mfd,felem,dn,(char *) buf.d);
                }
	   }
          else if (fmt=='F') { 
	     switch (dfmt) {
		case   8 : pf = buf.f; uc = bfdef->blank;
			   for (k=0; k<dn; k++, pf++, p.b++) {
                               if (bf && *p.b==uc) toNULLFLOAT(*pf);
                                  else *pf = *p.b;
                           }
	                   break;
		case -16 :
		case  16 : pf = buf.f; s = bfdef->blank;
			   for (k=0; k<dn; k++, pf++, p.s++) {
                               if (bf && *p.s==s) toNULLFLOAT(*pf);
                                  else *pf = *p.s;
                           }
	                   break;
		case  32 : pf = buf.f; i = bfdef->blank;
			   for (k=0; k<dn; k++, pf++, p.i++) {
                               if (bf && *p.i==i) toNULLFLOAT(*pf);
                                  else *pf = *p.i;
                           }
	                   break;
		case -32 : if (Midas_flag != 0)
                              MyPut(dfmt,felem,dn,(char *) p.f);
                           else
                              SCFPUT(mfd,felem,dn,(char *) p.f);
		           p.f += dn;
                           break;
		case -64 : pf = buf.f;
			   for (k=0; k<dn; k++, pf++, p.d++) {
                               if (isNULLDOUBLE(*p.d)) toNULLFLOAT(*pf);
                                  else *pf = *p.d;
                           }
                           break;
	     }
             if (Midas_flag != 0)
                {
	        if (dfmt!=-32)             /* store data in MIDAS file */
	          MyPut(-32,felem,dn,(char *) buf.f);
                }
             else
                {
	        if (dfmt!=-32)             /* store data in MIDAS file */
	          SCFPUT(mfd,felem,dn,(char *) buf.f);
                }
	   }
          else {
	     switch (dfmt) {
	        case   8 : if (bf) { 
		              pb = p.b; uc = bfdef->blank;
			      for (k=0; k<dn; k++, pb++)
				if (*pb==uc) *pb = 0xFF;
			    }
                           if (Midas_flag != 0)
		              MyPut(dfmt,felem,dn,(char *) p.b);
                           else
		              SCFPUT(mfd,felem,dn,(char *) p.b);
		           p.b += dn;
                           break;
	        case -16 : pu = p.u; u = bfdef->blank; ps = p.s;
			   for (k=0; k<dn; k++, pu++, ps++) {
			      *pu = *ps + 32768;
			      if (bf && *pu==u) toNULLSHORT(*pu);
			    }
                           if (Midas_flag != 0)
		              MyPut(dfmt,felem,dn,(char *) p.u);
                           else
		              SCFPUT(mfd,felem,dn,(char *) p.u);
		           p.u += dn;
                           break;
	        case  16 : if (bf) {
		              ps = p.s; s = bfdef->blank;
			      for (k=0; k<dn; k++, ps++)
				if (*ps==s) toNULLSHORT(*ps);
			    }
                          if (Midas_flag != 0)
                              MyPut(dfmt,felem,dn,(char *) p.s);
                           else
		              SCFPUT(mfd,felem,dn,(char *) p.s);
		           p.s += dn;
                           break;
	        case  32 : if (bf) {
		              pi = p.i; i = bfdef->blank;
			      for (k=0; k<dn; k++, pi++)
				if (*pi==i) toNULLINT(*pi);
			    }
                          if (Midas_flag != 0)
                              MyPut(dfmt,felem,dn,(char *) p.i);
                           else
                              SCFPUT(mfd,felem,dn,(char *) p.i);
		           p.i += dn;
                           break;
	        case -32 : if (Midas_flag != 0)
                              MyPut(dfmt,felem,dn,(char *) p.f);
                           else
                              SCFPUT(mfd,felem,dn,(char *) p.f);
		           p.f += dn;
                           break;
	        case -64 : if (Midas_flag != 0)
                              MyPut(dfmt,felem,dn,(char *) p.d);
                           else
                              SCFPUT(mfd,felem,dn,(char *) p.d);
		           p.d += dn;
                           break;
	     }
	   }
	  if (fmm) {       /* find data min/max values if not known */
	    switch (dsf) {
	       case   8 : pb = p.b - dn;
                          if (felem==1) mnvi = mxvi = *pb;
			  for (k=0; k<dn; k++, pb++)
			    if (*pb!=0xFF) {
			      if ((int)(*pb)<mnvi) mnvi = *pb;
			      else if (mxvi<(int)(*pb)) mxvi = *pb;
			    }
                          break;
	       case -16 : pu = p.u - dn;
                          if (felem==1) mnvi = mxvi = *pu;
			  for (k=0; k<dn; k++, pu++)
			    if (!isNULLSHORT(*pu)) {
			      if ((double)(*pu)<mnvi) mnvi = *pu;
			      else if (mxvi<(double)(*pu)) mxvi = *pu;
			    }
                          break;
	       case  16 : ps = p.s - dn;
                          if (felem==1) mnvi = mxvi = *ps;
			  for (k=0; k<dn; k++, ps++)
			    if (!isNULLSHORT(*ps)) {
			      if (*ps<mnvi) mnvi = *ps;
			      else if (mxvi<*ps) mxvi = *ps;
			    }
                          break;
	       case  32 : pi = p.i - dn;
                          if (felem==1) mnvi = mxvi = *pi;
			  for (k=0; k<dn; k++, pi++)
			    if (!isNULLINT(*pi)) {
			      if (*pi<mnvi) mnvi = *pi;
			      else if (mxvi<*pi) mxvi = *pi;
			    }
                          break;
	       case -32 : pf = (fpf) ? buf.f : p.f - dn;
                          if (felem==1) mnvd = mxvd = *pf;
			  for (k=0; k<dn; k++, pf++)
			    if (!isNULLFLOAT(*pf)) {
			      if (*pf<mnvd) mnvd = *pf;
			      else if (mxvd<*pf) mxvd = *pf;
			    }
                          break;
	       case -64 : pd = (bfdef->sflag) ? buf.d : p.d - dn;
                          if (felem==1) mnvd = mxvd = *pd;
			  for (k=0; k<dn; k++, pd++)
			    if (!isNULLDOUBLE(*pd)) {
			      if (*pd<mnvd) mnvd = *pd;
			      else if (mxvd<*pd) mxvd = *pd;
			    }
                          break;
	    }
	  }
          felem += dn;
          if (!dno) {
             gno++;
             pno = bfdef->pcount; pcnt = 0;
             pp = bfdef->parm;
             dno = ndata;
          }
       }
    } while (npix && gno<bfdef->gcount);
  }
  if (fmm) {                         /* save computed min/max values  */
     if (dsf > -32) 
        {
        cuts[0] = mnvi; cuts[1] = mxvi;		/* handle also UI2 here */
        }
     else {
	if (MAXFLOAT<mnvd) mnvd = MAXFLOAT;
	  else if (mnvd<MINFLOAT) mnvd = MINFLOAT;
	if (MAXFLOAT<mxvd) mxvd = MAXFLOAT;
	  else if (mxvd<MINFLOAT) mxvd = MINFLOAT;
	cuts[0] = mnvd; cuts[1] = mxvd;
      }
    SCDWRR(mfd,"LHCUTS",cuts,3,2,unit);
  }


  if (Midas_flag == 0)
     {
     if (0<=mfd) SCFCLO(mfd);                       /* close data files  */
     if (0<=mfdt) { TCSINI(mfdt); TCTCLO(mfdt); }
     mfd = -1; mfdt = -1;
     }

  return 0;
}