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

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT:  Copyright (c) 1987 European Southern Observatory,
                                         all rights reserved
.IDENTIFIER  module CUTIL
.LANGUAGE    C
.AUTHOR      K. Banse			IPG-ESO Garching
.KEYWORDS    
.COMMENTS    holds  FRAMOU_C, BLANKO_C, OSCHAR_C, ReadASCI, Newsort

.ENVIRONment MIDAS
             #include <midas_def.h>      Prototypes for MIDAS interfaces
 
.VERSIONS    1.00       940411  converted to C from UTIL.FOR    RvH
	     1.10       940901  add GetWC routine KB

 010423		last modif
 
------------------------------------------------------------*/

/* Define _POSIX_SOURCE to indicate that this is a POSIX program */

#define  _POSIX_SOURCE 1

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

#include <osparms.h>
#include <midas_def.h>

#define MAXDIM 3
#define   SWAP(a,b)  temp=(a);(a)=(b);(b)=temp; 

/*

*/
 
/*++++++++++++++++++++++++++++++
.PURPOSE     display name of image and its data type
.ALGORITHM   use STFINF to get data type
.INPUT/OUTPUT
  call as    FRAMOU_C(frame)
    
  input:
          char *frame		name of data frame

.RETURNS     nothing
------------------------------*/

#ifdef __STDC__
void FRAMOU_C(char *frame)

#else
void FRAMOU_C(frame)
char *frame;
#endif

{
int  nn, ibuf[6];

char cbuff[4], file[80], output[100];


/* make sure that the frame name ends with a null character */

if ((int)strlen(frame) > 82)
   {
   (void) strncpy(output,frame,82);
   output[82] = '\0';
   (void) CGN_CUTOFF(output,file);
   }
else
   (void) CGN_CUTOFF(frame,file);

nn = (int)strlen(file);
if (nn > 66) 
   {
   file[nn++] = '\n';		/* too long: split the output */
   file[nn] = '\0';
   }


/* Get data format of frame, as well as file type */

(void) SCFINF(frame,7,ibuf);
if (ibuf[1] == D_R4_FORMAT)
   (void) strcpy(cbuff,"R4");
else if (ibuf[1] == D_I4_FORMAT)
   (void) strcpy(cbuff,"I4");
else if (ibuf[1] == D_I2_FORMAT)
   (void) strcpy(cbuff,"I2");
else if (ibuf[1] == D_I1_FORMAT)
   (void) strcpy(cbuff,"I1");
else if (ibuf[1] == D_UI2_FORMAT)
   (void) strcpy(cbuff,"UI2");
else if (ibuf[1] == D_R8_FORMAT)
   (void) strcpy(cbuff,"R8");
else
   {
   (void) sprintf(output,"frame: %s  unknown data type",file);
   SCTPUT(output);
   return;
   }

(void) sprintf(output,"frame: %s  (data = %s",file,cbuff);
if (ibuf[2] == 0)			/* Midas frame */
   {
   if (ibuf[5] == 2)
      (void) strcat(output,")");
   else if (ibuf[5] == 1)
      (void) strcat(output,") (desc = ZFormat)");
   else
      (void) strcat(output,") (desc = oFormat!!)");
   }
else   					/* FITS file */
   {
   if (ibuf[5] == 2)
      (void) strcat(output,", format = FITS)");
   else if (ibuf[5] == 1)
      (void) strcat(output,", format = FITS) (desc = ZFormat)");
   else
      (void) strcat(output,", format = FITS) (desc = oFormat!!)");
   }

SCTPUT(output);
}

/*

*/
 
/*++++++++++++++++++++++++++++++
.PURPOSE     remove spaces in input string
.INPUT/OUTPUT
  call as    BLANKO_C( string )
  in/output:
          char *string		input string

.RETURNS     nothing
------------------------------*/
#ifdef __STDC__
void BLANKO_C(char *string)

#else
void BLANKO_C(string)
char *string;
#endif

{
register char *pntr;

pntr = string;
while (*pntr != '\0')
   {
   if (*pntr != ' ') *string++ = *pntr;
   pntr++;
   }
*string = '\0';
}

/*+++++++++++++++++++++++++++++++
.PURPOSE     return the backslash character
.INPUT/OUTPUT
  call as    OSCHAR_C( cval )
   
  output:
          char *cval		backslash character

.RETURNS     nothing
-------------------------------*/

#ifdef __STDC__
void OSCHAR_C(char *cval)

#else
void OSCHAR_C(cval)
char *cval;
#endif

{
static char backsl = '\\';             
*cval = backsl;
}

/*

*/
 
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

.IDENTIFICATION
  function ReadASCI		version 1.00	951130
  K. Banse			ESO - Garching

.KEYWORDS
  ASCII data file 

.PURPOSE
  get data from ASCII file

.ALGORITHM
  read all records of an ASCII file, convert to real data and fill buffer

.INPUT/OUTPUT
  use as: void ReadASCI(infile,dattyp,total,a,b,minflg,fmin,fmax)

  input parameter:
  INFILE:	char *		name of ASCII data file
  DATTYP:	int		1 = integer, 2 = real,
				4 = double 
  TOTAL:	int		no. of pixels to fill
  MINFLG:	int		= 1 to also calculate MIN + MAX, else = 0

  output parameter:
  A:		float array	data buffer to be filled
  B:		int array	data buffer to be filled
  FMIN:		float *		minimum of data
  FMAX:		float *		maximum

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

#ifdef __STDC__
void ReadASCI(char *infile,int dattyp,int total,float *a,int *b,
              int minflg,float *fmin,float *fmax)

#else
void ReadASCI(infile,dattyp,total,a,b,minflg,fmin,fmax)
char  *infile;
int   dattyp, total, minflg, *b;
float *a, *fmin, *fmax;
#endif

{
int fb, no, ic, finis, uni;
int linecnt, reclen, inull, more;
int fid;
register int nr;

float rnull;

char *cptr, *ncptr, *oldrec, *newrec;

#define NOELM	2000
#define RECMAX	8000



(void) SCKRDR("NULL",2,1,&ic,&rnull,&uni,&no);     /* get user null value */

more = 0;
linecnt = 0;

fid = osaopen(infile,READ);
if (fid < 0)
   {
   char  errmess[80];

   (void) sprintf(errmess,"Problems opening data file %s",infile);
   SCETER(1,errmess);
   return;
   }

oldrec = malloc((unsigned int) RECMAX);		/* get space for file record */
newrec = malloc((unsigned int) RECMAX);

ic = 0;						/* init data counter */
finis = 0;

/*  ----------------------------- */
/*  reading loop for integer data */
/*  ----------------------------- */

if (dattyp == 1)     
   {
   int  *ispace, *iptr, bmin = 99999, bmax = -99999;
   register int  ir;
   float rdum = 0;
   double ddum = 0;

   inull = CGN_NINT(rnull);
   ispace = (int *) malloc((unsigned int) (NOELM*sizeof(int)));

sect_200:
   linecnt ++;
   reclen = osaread(fid,oldrec,RECMAX);
   if (reclen < 0) goto sect_930;		/* EOF encountered */
   if (reclen == 0) goto sect_200;		/* skip empty records */

   fb = 1;                      /* first-blank-flag is set in the beginning */
   ncptr = newrec;
   cptr = oldrec;
   for (nr=0; nr<reclen; nr++)
      {
      if (*cptr != ' ')     
         {
         *ncptr++ = *cptr;
         fb = 0;
         }
      else if (fb == 0)     
         {
         *ncptr++ = ',';
         fb = 1;
         }
      cptr++;
      }
   *ncptr = '\0';


   /*  convert to integer numbers */

   no = CGN_CNVT(newrec,1,NOELM,ispace,&rdum,&ddum);
   if (no < 1)
      {
      char  outbuf[80];

      (void) sprintf(outbuf,"no data in line no. %d",linecnt);
      SCTPUT(outbuf);
      goto sect_200;              /* skip comment lines */
      }

   iptr = ispace;
   if ((no+ic) > total)     
      {
      more = - (no + ic - total);       /* more goes negatively  */
      no = total - ic;                  /* cut off overflow data */
      finis = 1;
      }

   if (minflg == 1)     
      {
      for (nr=0; nr<no; nr++)
         {
         ir = *iptr++;
         b[ic++] = ir;
         if (ir > bmax)     
            bmax = ir;
         else if (ir < bmin)
            bmin =  ir;
         }
      }
   else
      {
      for (nr=0; nr<no; nr++)
         b[ic++] = *iptr++;
      }
   if (finis == 0) goto sect_200;            /* get next record */


   /* look if more data there + count */

sect_909:
   linecnt ++;
   reclen = osaread(fid,oldrec,RECMAX);
   if (reclen < 0) goto sect_930;               /* EOF encountered */
   if (reclen == 0) goto sect_909;              /* skip empty records */

   fb = 1;                      /* first-blank-flag is set in the beginning */
   cptr = oldrec;
   ncptr = newrec;
   for (nr=0; nr<reclen; nr++)
      {
      if (*cptr != ' ')
         {
         *ncptr++ = *cptr;
         fb = 0;
         }
      else if (fb == 0)
         {
         *ncptr++ = ',';
         fb = 1;
         }
      cptr++;
      }
   *ncptr = '\0';

   no = CGN_CNVT(newrec,1,NOELM,ispace,&rdum,&ddum);
   if (no < 1)
      {
      char  outbuf[80];

      (void) sprintf(outbuf,"no data in line no. %d",linecnt);
      SCTPUT(outbuf);
      goto sect_909;              /* skip comment lines */
      }

   more -= no;
   goto sect_909;


   /* end-of-file reached */

sect_930:
   (void) osaclose(fid);
   if (ic < total)		/* less data then required */
      {
      char outbuf[80];

      (void) sprintf(outbuf,
             "%d data values read in - %d values expected",ic,total);
      SCTPUT(outbuf);
      (void) sprintf(outbuf,"missing data set to %d",inull);
      SCTPUT(outbuf);
      for (nr=ic; nr<total; nr++)	/* fill with User Null */
         b[nr] = inull;
      }
   if (minflg == 1)                     /* return real min, max */
      {
      *fmin = bmin;
      *fmax = bmax;
      }

   (void) free(ispace);
   }

/*  ----------------------------- */
/*  reading loop for real data */
/*  ----------------------------- */

else if (dattyp == 2)
   {
   float  *rspace, *rptr, rmin = 99999.0, rmax = -99999.0;
   register float  rr;
   int    idum = 0;
   double ddum = 0;

   rspace = (float *) malloc((unsigned int) (NOELM*sizeof(float)));

sect_1000:
   linecnt ++;
   reclen = osaread(fid,oldrec,RECMAX);
   if (reclen < 0) goto sect_5000;               /* EOF encountered */
   if (reclen == 0) goto sect_1000;              /* skip empty records */

   fb = 1;                      /* first-blank-flag is set in the beginning */
   cptr = oldrec;
   ncptr = newrec;
   for (nr=0; nr<reclen; nr++)
      {
      if (*cptr != ' ')
         {
         *ncptr++ = *cptr;
         fb = 0;
         }
      else if (fb == 0)
         {
         *ncptr++ = ',';
         fb = 1;
         }
      cptr++;
      }
   *ncptr = '\0';


   /*  convert to real numbers */

   no = CGN_CNVT(newrec,2,NOELM,&idum,rspace,&ddum);
   if (no < 1) 
      {
      char  outbuf[80];

      (void) sprintf(outbuf,"no data in line no. %d",linecnt);
      SCTPUT(outbuf);
      goto sect_1000;              /* skip comment lines */
      }

   rptr = rspace;
   if ((no+ic) > total)
      {
      more = - (no + ic - total);       /* more goes negatively  */
      no = total - ic;                  /* cut off overflow data */
      finis = 1;
      }

   if (minflg == 1)
      {
      for (nr=0; nr<no; nr++)
         {
         rr = *rptr++;
         a[ic++] = rr;
         if (rr > rmax)
            rmax = rr;
         else if (rr < rmin)
            rmin = rr;
         }
      }
   else
      {
      for (nr=0; nr<no; nr++)
         a[ic++] = *rptr++;
      }
   if (finis == 0) goto sect_1000;            /* get next record */


   /* look if more data there + count  */

sect_3000:
   linecnt ++;
   reclen = osaread(fid,oldrec,RECMAX);
   if (reclen < 0) goto sect_5000;               /* EOF encountered */
   if (reclen == 0) goto sect_3000;              /* skip empty records */

   fb = 1;                      /* first-blank-flag is set in the beginning */
   cptr = oldrec;
   ncptr = newrec;
   for (nr=0; nr<reclen; nr++)
      {
      if (*cptr != ' ')
         {
         *ncptr++ = *cptr;
         fb = 0;
         }
      else if (fb == 0)
         {
         *ncptr++ = ',';
         fb = 1;
         }
      cptr++;
      }
   *ncptr = '\0';

   no = CGN_CNVT(newrec,2,NOELM,&idum,rspace,&ddum);
   if (no < 1)
      {
      char  outbuf[80];

      (void) sprintf(outbuf,"no data in line no. %d",linecnt);
      SCTPUT(outbuf);
      goto sect_3000;              /* skip comment lines */
      }

   rptr = rspace;

   more -= no;
   goto sect_3000;


   /* end-of-file reached */

sect_5000:
   (void) osaclose(fid);
   if (ic < total)		/* less data then required */
      {
      char outbuf[80];

      (void) sprintf(outbuf,
             "%d data values read in - %d values expected",ic,total);
      SCTPUT(outbuf);
      (void) sprintf(outbuf,"missing data set to %f",rnull);
      SCTPUT(outbuf);
      for (nr=ic; nr<total; nr++)		/* fill with user Null */
         a[nr] = rnull;
      }

   if (minflg == 1)                     /* return real min, max */
      {
      *fmin = rmin;
      *fmax = rmax;
      }

   (void) free(rspace);
   }

/*  ----------------------------- */
/*  reading loop for double data */
/*  ----------------------------- */

else
   {
   double  *dspace, *dptr, dmin = 99999.0, dmax = -99999.0;
   double  *ddat, dnull;
   register double  dd;
   float  rdum = 0;
   int    idum = 0;

   dspace = (double *) malloc((unsigned int) (NOELM*sizeof(double)));
   ddat = (double *) a;

sect_11000:
   linecnt ++;
   reclen = osaread(fid,oldrec,RECMAX);
   if (reclen < 0) goto sect_5000;               /* EOF encountered */
   if (reclen == 0) goto sect_1000;              /* skip empty records */

   fb = 1;                      /* first-blank-flag is set in the beginning */
   cptr = oldrec;
   ncptr = newrec;
   for (nr=0; nr<reclen; nr++)
      {
      if (*cptr != ' ')
         {
         *ncptr++ = *cptr;
         fb = 0;
         }
      else if (fb == 0)
         {
         *ncptr++ = ',';
         fb = 1;
         }
      cptr++;
      }
   *ncptr = '\0';


   /*  convert to double numbers */

   no = CGN_CNVT(newrec,4,NOELM,&idum,&rdum,dspace);
   if (no < 1)
      {
      char  outbuf[80];

      (void) sprintf(outbuf,"no data in line no. %d",linecnt);
      SCTPUT(outbuf);
      goto sect_11000;              /* skip comment lines */
      }

   dptr = dspace;
   if ((no+ic) > total)
      {
      more = - (no + ic - total);       /* more goes negatively  */
      no = total - ic;                  /* cut off overflow data */
      finis = 1;
      }

   if (minflg == 1)
      {
      for (nr=0; nr<no; nr++)
         {
         dd = *dptr++;
         ddat[ic++] = dd;
         if (dd > dmax)
            dmax = dd;
         else if (dd < dmin)
            dmin = dd;
         }
      }
   else
      {
      for (nr=0; nr<no; nr++)
         ddat[ic++] = *dptr++;
      }
   if (finis == 0) goto sect_11000;            /* get next record */


   /* look if more data there + count  */

sect_13000:
   linecnt ++;
   reclen = osaread(fid,oldrec,RECMAX);
   if (reclen < 0) goto sect_15000;               /* EOF encountered */
   if (reclen == 0) goto sect_13000;              /* skip empty records */

   fb = 1;                      /* first-blank-flag is set in the beginning */
   cptr = oldrec;
   ncptr = newrec;
   for (nr=0; nr<reclen; nr++)
      {
      if (*cptr != ' ')
         {
         *ncptr++ = *cptr;
         fb = 0;
         }
      else if (fb == 0)
         {
         *ncptr++ = ',';
         fb = 1;
         }
      cptr++;
      }
   *ncptr = '\0';

   no = CGN_CNVT(newrec,4,NOELM,&idum,&rdum,dspace);
   if (no < 1)
      {
      char  outbuf[80];

      (void) sprintf(outbuf,"no data in line no. %d",linecnt);
      SCTPUT(outbuf);
      goto sect_3000;              /* skip comment lines */
      }

   dptr = dspace;

   more -= no;
   goto sect_13000;


   /* end-of-file reached */

sect_15000:
   (void) osaclose(fid);
   if (ic < total)              /* less data then required */
      {
      char outbuf[80];

      (void) sprintf(outbuf,
             "%d data values read in - %d values expected",ic,total);
      SCTPUT(outbuf);
      (void) sprintf(outbuf,"missing data set to %f",rnull);
      SCTPUT(outbuf);
      dnull = (double) rnull;
      for (nr=ic; nr<total; nr++)               /* fill with user Null */
         ddat[nr] = dnull;
      }

   if (minflg == 1)                     /* return real min, max */
      {
      *fmin = (float) dmin;
      *fmax = (float) dmax;
      }

   (void) free(dspace);
   }


(void) free(oldrec);
(void) free(newrec);


/*  fill keyword NULL(1) */

rnull = more;
(void) SCKWRR("NULL",&rnull,1,1,&uni);
}

/*

*/
 
/* the following routine is needed for OSF1/DEC Unix on Alpha 
   which have a problem in FORTRAN with (-1) ** 2   !!!       */

#ifdef __STDC__
void MyPower(float *consta, float *constb, float *constc)

#else
void MyPower(consta,constb,constc)
float  *consta, *constb, *constc;
#endif

{
double da, db, dc;

da = (double) *consta;
db = (double) *constb;

dc = pow(da,db);
*constc = (float) dc;
}

/*

*/
 
#ifdef __STDC__
void Newsort(float *arr,int nsize,int nmed,float *fmedian)

#else
void Newsort(arr,nsize,nmed,fmedian)
float  *arr;		/* IN: pointer to float array - 1 !! */
int    nsize;		/* IN: size of arr */
int    nmed;		/* IN: indx of median - FORTRAN counting (nsize+1)/2 */
float  *fmedian;	/* OUT: value of median pixel */
#endif


{
register int  i, ir, j, l, ll, mid;

register float    a, temp;


l = 1;					/* arr actually points to 1.elem - 1 */
ir = nsize;

for (;;)
   {
   if (ir <= (l+1))
      {
      if ( (ir == (l+1)) && (arr[ir] < arr[l]) ) 
         {
         SWAP(arr[l],arr[ir])
         }
      *fmedian = arr[nmed];
      return;
      } 

   else 
      {
      mid = (l+ir) >> 1;
      SWAP(arr[mid],arr[l+1])
      if (arr[l] > arr[ir]) 
         {
         SWAP(arr[l],arr[ir])
         }
      if (arr[l+1] > arr[ir]) 
         {
         SWAP(arr[l+1],arr[ir])
         }
      if (arr[l] > arr[l+1]) 
         {
         SWAP(arr[l],arr[l+1])
         }
      i = ll = l + 1;
      j = ir;
      a = arr[ll];

      for (;;) 
         {
         do i++; while (arr[i] < a);
         do j--; while (arr[j] > a);
         if (j < i) break;

         SWAP(arr[i],arr[j])
         }

      arr[ll] = arr[j];
      arr[j] = a;
      if (j >= nmed) ir = j-1;
      if (j <= nmed) l = i;
      }
   }
}