/* @(#)fitscdm.c	17.1.1.1 (ES0-DMD) 01/25/02 17:39:08 */
/*===========================================================================
  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      fitscdm.c
.LAUGUAGE   C
.AUTHOR     P.Grosbol   ESO/IPG
.KEYWORDS   FITS data matrix, decode, copy, buffer
.COMMENT    copy prime data matrix of FITS file to user buffer
.VERSION    1.0  1993-Mar-25 : Creation,   PJG 
.VERSION    1.1  1993-Oct-28 : Update to new SC + prototypes,   PJG 
.VERSION    1.2  1994-May-09 : Cast variables for ANSI-C,   PJG 
---------------------------------------------------------------------*/
#include   <computer.h>            /* computer specific constants    */
#include   <fitsfmt.h>             /* general data definitions       */
#include   <fitsdef.h>             /* basic FITS definitions         */
#include   <midas_def.h>

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

static     int                pno; /* Parameter count                */
static     int                gno; /* Group count                    */
static     int                dno; /* Data array size                */
static     int               dfmt; /* Data format i.e. BITPIX        */
static     int              sflag; /* Scale flag                     */
static     int              bflag; /* Blank flag                     */
static     int              blank; /* Blank value                    */
static     double             fac; /* Scale factor for pixel values  */
static     double            zero; /* Zero offset for pixel values   */

static     long int          size; /* Total no. of pixels in matrix  */
static     int                nvp; /* No. of valid pixels in buffer  */
static     float       fbuf[MXBF]; /* Internal data buffer           */

int fitsidm(bfdef)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE  Initiate copying of prime data matrix to user buffer
.RETURN   error code - 0: OK, -1: error
---------------------------------------------------------------------*/
BFDEF     *bfdef;              /* IN: pointer to FITS header parm's  */
{
  int     n;

  nvp = 0;                     /* Initiate static variables          */
  pno = bfdef->pcount;
  gno = bfdef->gcount;
  fac = bfdef->bscale;
  zero = bfdef->bzero;
  dfmt = bfdef->bitpix;
  sflag = bfdef->sflag;
  bflag = bfdef->bflag;
  blank = bfdef->blank;

  dno = 1;                     /* Compute size of data martix        */
  for (n=0; n<bfdef->naxis; n++)
    dno *= bfdef->data[n].naxis;
  size = gno * (pno + dno);

  if (pno!=0 || gno!=1) return -1;

  return (dno) ? 0 : -1;
}

int fitscdm(pbuf, no)
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.PURPOSE  Copy prime data matrix to user buffer
.RETURN   No of values copied; -1: error
---------------------------------------------------------------------*/
float     *pbuf;               /* IN: Pointer to user buffer         */
int       no;                  /* IN: No. of pixels to copy          */
{
  register int  bf;
  unsigned char uc;
  short         s;
  int           npix, k, n;
  int           i, dn, pn, npt;
  float         *pf;
  double        d;
  union { 
	  unsigned char   *b;
	  char            *c;
	  short           *s;
	  int             *i;
	  float           *f;
	  double          *d;
        } p;

  if (pno!=0 || gno!=1) return -1;   /* check if simple data matrix  */
  if (dno<no) no = dno;

  bf = bflag;
  npt = 0;
  npix = 0;
  if (nvp && 0<no) {               /* if values in internal buffer   */
     n = (nvp<no) ? nvp : no;
     pf = fbuf;
     k = n;
     while (k--) *pbuf++ = *pf++;
     nvp -= n;
     for (k=0; k<nvp; k++) fbuf[k] = *pf++;
     no -= n;
     npt += n;
   }

  while (0<no) {                              /* read remaining data */
    if ((n=dread(&p.c,FITSLR))!=FITSLR) break;   /* get next record  */

    switch (dfmt) {                /* convert to local data format   */
       case   8 : npix = n; break;
       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;
    }

    dn = (npix<no) ? npix : no;
    if (sflag)
      switch (dfmt) {
	 case   8 : uc = blank;
	            for (k=0; k<dn; k++, pbuf++, p.b++) {
		       if (bf && *p.b==uc) toNULLFLOAT(*pbuf);
		       else *pbuf = fac*((int)(*p.b)) + zero;
		     }
	            break;
	 case  16 : s = blank;
	            for (k=0; k<dn; k++, pbuf++, p.s++) {
		       if (bf && *p.s==s) toNULLFLOAT(*pbuf);
		       else *pbuf = fac*(*p.s) + zero;
		     }
	            break;
	 case  32 : i = blank;
	            for (k=0; k<dn; k++, pbuf++, p.i++) {
		       if (bf && *p.i==i) toNULLFLOAT(*pbuf);	
		       else *pbuf = fac*(*p.i) + zero;
		     }
	            break;
	 case -32 : for (k=0; k<dn; k++)
	               *pbuf++ = fac * (*(p.f++)) + zero;
	            break;
	 case -64 : for (k=0; k<dn; k++)
	               *pbuf++ = fac * (*(p.d++)) + zero;
	            break;
       }
    else
      switch (dfmt) {
	 case   8 : uc = blank;
	            for (k=0; k<dn; k++, pbuf++, p.b++) {
		       if (bf && *p.b==uc) toNULLFLOAT(*pbuf);
		       else *pbuf = *p.b;
		     }
	            break;
	 case  16 : s = blank;
	            for (k=0; k<dn; k++, pbuf++, p.s++) {
		       if (bf && *p.s==s) toNULLFLOAT(*pbuf);
		       else *pbuf = *p.s;
		     }
	            break;
	 case  32 : i = blank;
	            for (k=0; k<dn; k++, pbuf++, p.i++) {
		       if (bf && *p.i==i) toNULLFLOAT(*pbuf);
		       else *pbuf = *p.i;
		     }
	            break;
	 case -32 : for (k=0; k<dn; k++) *pbuf++ = *p.f++;
	            break;
	 case -64 : for (k=0; k<dn; k++, pbuf++, p.d++) {
	               if (isNULLDOUBLE(*p.d)) toNULLFLOAT(*pbuf);
		       else *pbuf = *p.d;
		     }
	            break;
       }
    dno -= dn;
    npix -= dn;
    no -= dn;
    npt +=dn;
  }

  if (npix) {                         /* something left in buffer   */
     pf = fbuf;
     nvp = npix;
     switch (dfmt) {
	case   8 : uc = blank;
	           for (k=0; k<npix; k++, pf++, p.b++) {
		      if (bf && *p.b==uc) toNULLFLOAT(*pf);
		      else *pf = fac*((int)(*p.b)) + zero;
		    }
	           break;
	case  16 : s = blank;
	           for (k=0; k<npix; k++, pf++, p.s++) {
		      if (bf && *p.s==s) toNULLFLOAT(*pf);
		      else *pf = fac*(*p.s) + zero;
		    }
	           break;
	case  32 : i = blank;
	           for (k=0; k<npix; k++, pf++, p.i++) {
		      if (bf && *p.i==i) toNULLFLOAT(*pf);	
		      else *pf = fac*(*p.i) + zero;
		    }
	           break;
	case -32 : for (k=0; k<npix; k++, pf++, p.f++) {
	              if (isNULLFLOAT(*p.f)) toNULLFLOAT(*pf);
	              else *pf = fac*(*p.f) + zero;
		    }
	           break;
	case -64 : for (k=0; k<npix; k++, pf++, p.d++) {
	              if (isNULLDOUBLE(*p.d)) toNULLFLOAT(*pf);
	              else *pf = fac*(*p.d) + zero;
		    }
	           break;
       }
   }

  return npt;
}