/*****************************************************************************/
/* util.c                                  Jan Grinder-Pedersen - EMI,DTU,DK */
/*****************************************************************************/
#define VERSION 2.03

#ifdef apollo
#define APOLLO  1
#else
#define APOLLO  0
#endif

#ifdef __hp9000s700
#define HP700   1
#else
#define HP700   0
#endif

#ifdef __sgi
#define SGI	1
#else 	
#define	SGI	0
#endif

#if (APOLLO || HP700 || SGI) /* If Apollo,  HP700 or SGI  Platform. */
#define _UNIX_ 
#endif /* (APOLLO || HP700 || SGI) */

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

#ifdef __cplusplus 
  #include  <malloc.h>
#endif

#define CHAR  char                 /*  8 bit character.           */
#define SBYTE   signed char        /*  8 bit signed integer.      */
#define UBYTE unsigned char        /*  8 bit unsigned integer.    */
#define SWORD   signed short int   /* 16 bit signed integer.      */
#define UWORD unsigned short int   /* 16 bit unsigned integer.    */

#ifdef _UNIX_
#define SINT    signed int         /* 32 bit signed integer.      */
#define UINT  unsigned int         /* 32 bit unsigned integer.    */
#else
#define SINT    signed long int    /* 32 bit signed integer.      */
#define UINT  unsigned long int    /* 32 bit unsigned integer.    */
#endif /* _UNIX_ */

#define SREAL float                /* 32 bit floating point real. */
#define DREAL double               /* 64 bit floating point real. */

#define LINT  long int             /* Long integer used for io.   */
                   
#define PI              3.14159265358979323846

#define ABS(a)     (((a) > 0)   ? (a) : (-(a)))
#define MAX(a,b)   (((a) > (b)) ? (a) : (b))
#define MIN(a,b)   (((a) < (b)) ? (a) : (b))
#define BRACKET(a,x,b) ( ((b)<(x)) ? (b) : (((x)<(a)) ? (a) : (x)) )
#define SQUARE(a)  ((a)*(a))
#define ROUND(a)   (((a) > 0)   ? (SINT)((a)+0.5) : (SINT)((a)-0.5) )
#define ATAN2(a,b) ((((a)==0)&&((b)==0)) ? 0.0 : atan2((a),(b)))

#define CHCK_INIT(a) CHAR *pid=a; SWORD rval=0;
#define EXIT(a,b)    {ErrMes(pid,a,b);rval=-1; goto nice_exit;}
#define EXIT2(a,b)    {ErrMes(pid,a,b);rval=-1;exit(EXIT_FAILURE);}
#define WARN(a,b)    {WarMes(pid,a,b);}

#ifdef _UNIX_ 
#define SSCANF_D(a,b)  sscanf((a),"%d",(b));
#else
#define SSCANF_D(a,b)  sscanf((a),"%ld",(b));
#endif /* _UNIX_ */

/* Union used to extract or set the first 16 bit of a 32 bit float. */
typedef union {
  SREAL r;
  SWORD s[2];
} RS_U;

typedef union {
  UWORD s;
  UBYTE b[2];
} SB_U;

/* Error types. */
typedef enum {
  E_FOPE,          /* File open error.        */
  E_FCLO,          /* File closing error.     */
  E_FPOS,          /* File positioning error. */
  E_FREA,          /* File reading error.     */
  E_FWRI,          /* File writting error.    */
  E_ALLO,          /* Allocation error.       */
  E_DIVZ,          /* Division by zero.       */
  E_TRAC,          /* Trace back.             */
  E_HEAD,          /* Wrong Header.           */
  E_LARG,          /* Arg. too large.         */
  E_SMAL,          /* Arg. too small.         */
  E_NEXI,          /* Non-exist item.         */
  E_NIMP,          /* Non-impl. func.         */
  E_SCAN,          /* Scan input error.       */
  E_SPEC,          /* Special error.          */
  E_TYPE           /* Unknown datatype        */
} ERR_T;

typedef enum {
  SINC,
  RCOS,
  IDMP
} FILTER_TYPE;

typedef enum {
  LOG10,
  POWER,
  LIN
} TRANSFORMATION_TYPE;

enum data_type {
  ubyte, sbyte, uword, sword, uInt, sInt, half_sreal, sreal, dreal
};

enum parm_type {
  amp, pwr, rad, deg
};
  
enum Format {
  Dos, Unix
};

enum Format machine_type;

DREAL round(DREAL inp);
void set_fixpix(DREAL *fixPix, enum data_type otype);  
void log_transform(DREAL *arr, SINT len, enum parm_type oPrmType, DREAL minPix); 
void root_transform(DREAL *arr, SINT len, DREAL pivot, DREAL root, DREAL minPix); 
void linear_transform(DREAL *arr, SINT len, DREAL minPix, DREAL maxPix, 
                      DREAL type_min, DREAL type_max);
void set_min_max_out(DREAL *omin, DREAL *omax, enum data_type dtype);
void find_type_min_max(DREAL *typeMin, DREAL *typeMax, enum data_type dtype);
void f_length(FILTER_TYPE type, DREAL scale, SINT *iflen, SINT *oflen);
void f_init(DREAL *filter, DREAL scale, SINT over_samp, SINT klen,
            FILTER_TYPE filter_type);
void float_array_power_real(DREAL *arr, SINT len);
void float_array_power(DREAL *iArr, DREAL *oArr, SINT len);
void float_array_sqrt(DREAL *arr, SINT len);
void float16_to_float64(SWORD *sbuf, DREAL *rbuf, SINT length, SINT step, 
                        SINT icr, DREAL inputScale);
void float64_to_float16(DREAL *rbuf, SWORD *sbuf, SINT length, SINT step, 
                        SINT icr);
SINT read_line(DREAL *iarr, UBYTE *ibufp, SINT step, 
               enum data_type itype, SINT size, 
               SINT nSam, SINT icr, enum Format format,
               DREAL inputScale);
SINT write_line(DREAL *oarr, UBYTE *obufp, SINT step, 
                enum data_type otype, SINT size, SINT nSam, SINT icr, 
                enum Format oFormat);
void nice_saturation(DREAL *arr, enum data_type dtype, SINT length, 
                     SINT *uSat, SINT *lSat);
void ErrMes(const CHAR *pid, ERR_T error, CHAR *name);
void WarMes(const CHAR *pid, ERR_T error, CHAR *name);
void update_menu(char *iNam, enum data_type itype, SINT icr, 
                 SINT nSam, SINT nLin,
                 enum Format iFormat, DREAL inputScale,
                 char *oNam, enum data_type otype, SINT ocr, 
                 enum Format oFormat, SINT win, SINT detect, 
                 SINT average, SINT pix, enum parm_type iPrmType);
void update_type_menu(void);
void type_menu(enum data_type *itype, SINT *icr, SINT *isize);
void window_menu(SINT *win, SINT *uSam, SINT *uLin, SINT *saOf, SINT *liOf,
                 SINT iSam, SINT iLin);
void update_window_menu(SINT win, SINT uSam, SINT uLin, SINT saOf, SINT liOf);
void detect_menu(SINT *detect, enum parm_type *oPrmType);
void average_menu(SINT *average, DREAL *scaleSam, DREAL *scaleLin,
                  FILTER_TYPE *filter_type);
void filter_menu(FILTER_TYPE *filter_type);
void update_pix_otype_menu(void);
void pix_otype_menu(enum data_type *otype, SINT *isize,DREAL *fixPix, 
                    DREAL *minOut, DREAL *maxOut);
void pix_menu(SINT *pix, TRANSFORMATION_TYPE *transType, enum data_type *otype, 
              SINT *osize, SINT *ocr, DREAL *minPix, DREAL *maxPix, 
              DREAL *powerPix, DREAL *fixPix, enum parm_type oPrmType,
              DREAL *minOut, DREAL *maxOut);
void pix_phase_menu(SINT *pix, TRANSFORMATION_TYPE *transType, enum data_type *otype, 
                    SINT *osize, SINT *ocr, DREAL *minPix, DREAL *maxPix,
                    DREAL *minOut, DREAL *maxOut, DREAL *fixPix);
char *parm_type_str(enum parm_type type);

/* Averaging definitions */
#define   OVER_SAMP  128
#define   MIN_F_LEN    8 

/*****************************************************************************/
/* Begin MAIN                                                                */
/*****************************************************************************/
SINT main (int argc, char *argv[])
{
  CHCK_INIT("main")
  
  SB_U test;

  /* I/O buffers */ 
  UBYTE  *ibuf = NULL;   
  UBYTE  *obuf = NULL;

  /* File names and pointers. */
  char   iNam[128], oNam[128];
  SINT   fPos;
  FILE   *iFp = NULL, *oFp = NULL;

  /* Common variables. */
  SINT   iline, oline, iSam = 1024 , iLin = 1024, i,
         icr = 2, ocr = 2, isize = 2, osize = 2, lSat = 0, uSat = 0,
         oSam = 1024, oLin = 1024;
  DREAL  *iarr = NULL, *oarr = NULL, inputScale = 1;
  enum data_type itype = half_sreal, otype = half_sreal;
  enum Format    iFormat = Unix; 
  #ifdef _UNIX_
  enum Format oFormat = Unix;
  #else
  enum Format oFormat = Dos;
  #endif /* _UNIX_ */
  
  /* Window variables. */
  SINT   win = 0, uSam = 1024, uLin = 1024, liOf = 0, saOf = 0;  
  
  /* Detection variables. */
  SINT detect = 0; 
  enum parm_type iPrmType = amp;
  enum parm_type oPrmType = amp;
  
  /* Averaging variables. */
  SINT   average = 0, saklen = 0, liklen = 0, isaflen = 0, osaflen = 0, 
         iliflen = 0, oliflen = 0, sa_int, li_int, k, n, samp, icrAv,
         saAv, liAv;
  DREAL  scaleSam = 1, scaleLin = 1;
  DREAL  *iarrp  = NULL, *arr    = NULL, *arrp  = NULL,
         *safilt = NULL, *lifilt = NULL, *filtp = NULL;
  DREAL  sa_frac, li_frac, li_offset, sa_offset;
  FILTER_TYPE filter_type = SINC;
  
  /* Pixel variables. */
  SINT   pix = 0; 
  TRANSFORMATION_TYPE transType = POWER;
  DREAL  minPix = -36, maxPix = 0, powerPix = 0.8, fixPix = 128; 
  DREAL  minOut, maxOut;
  
  /* Menu variables. */
  SINT   menu_ok = 0, choice = 0, max_choice = 15;
  char   chline[16],ctmp;

  /***************************************************************************/
  /* The code that does something starts here.                               */
  /***************************************************************************/

  fprintf(stderr, "util v. %4.2f \n", VERSION);

  test.s = 1;
  if (test.b[1]) machine_type = Unix;
  else           machine_type = Dos;

  strcpy(iNam, "");
  strcpy(oNam, "");
  if (argc > 1) strcat(iNam, *(argv+1));
  if (argc > 2) strcat(oNam, *(argv+2));

  /***************************************************************************/
  /* Mainloop - (Do forever loop)                                            */
  /***************************************************************************/
  do {
 
    /*************************************************************************/
    /* Section 1: Menu Loop                                                  */
    /*************************************************************************/
    do { 

      /*---------------------------------------------------------------------*/
      /* Write main menu, and prompt for choice.                             */
      /*---------------------------------------------------------------------*/
      update_menu(iNam, itype, icr, iSam, iLin, iFormat, inputScale, oNam, 
                 otype, ocr, oFormat, win, detect, average, pix, iPrmType);
      printf("\n");
      do {
        printf("Choose (1) to (%i) for changing, (%i) for execute or (0) for exit: ", 
             (int)max_choice-1, (int)max_choice);
        if (scanf("%s",chline) == EOF) goto final_exit;
        choice = -1;
        SSCANF_D(chline, &choice);
      } while ((choice < 0) || (choice > max_choice));

      /*---------------------------------------------------------------------*/
      /* Decode choice.                                                      */
      /*---------------------------------------------------------------------*/
      switch(choice) {
      case 1:  
        printf("\nInput file name : ");
        if (scanf("%s", iNam) == EOF) EXIT2(E_SCAN, "iNam"); 
        break;
      case 2:
        type_menu(&itype, &icr, &isize);
        if (icr == 2 && ocr == 1 && !detect) {
          detect = 1;
          printf("\n\nWARNING: Detection turned on due to complex input and real output \n");
        };
        break;
      case 3: 
        if (iPrmType==pwr) iPrmType = amp; 
        else               iPrmType = pwr;
        pix = 0;
        break;
      case 4: 
        printf("\nNumber of samples per line : ");
        if (scanf("%i", &iSam) == EOF) EXIT2(E_SCAN, "iSam");
        if (iSam <= 0) iSam = 1024;
        win = 0;
        uSam = iSam;
        saOf = 0;
        break;
      case 5: 
        printf("\nNumber of lines : ");
        if (scanf("%i", &iLin) == EOF) EXIT2(E_SCAN, "iLin"); 
        if (iLin <= 0) iLin = 1024;
        win = 0;
        uLin = iLin;
        liOf = 0;
        break;
      case 6: 
        if (iFormat == Dos) iFormat = Unix; 
        else                iFormat = Dos; 
        break;
      case 7:
        printf("\nOutput file name : ");
        if (scanf("%s", oNam) == EOF) EXIT2(E_SCAN, "oNam"); 
        break;
      case 8:
        type_menu(&otype, &ocr, &osize);
        /*if (ocr == 2 || otype > 5)*/ pix = 0;
        if (ocr == 1 && otype < 6) set_fixpix(&fixPix, otype); 
        if (ocr == 1 && otype < 6) set_min_max_out(&minOut,&maxOut, otype); 
        if (icr == 2 && ocr == 1 && !detect) {
          detect = 1;
          printf("\n\nWARNING: Detection turned on due to complex input and real output \n");
        };
        break;
      case 9: 
        if (oFormat == Dos) oFormat = Unix; 
        else                oFormat = Dos; 
        break;
      case 10: 
        printf("\nInput Scale factor : ");
        if (scanf("%lf", &inputScale) == EOF) EXIT2(E_SCAN, "inputScale"); 
        if (inputScale == 0) WARN(E_SPEC, "Input data will be multiplied by zero");
        break;
      case 11: 
        window_menu(&win, &uSam, &uLin, &saOf, &liOf, iSam, iLin); 
        break;
      case 12:
        detect_menu(&detect, &oPrmType);
        if (!detect) {
          if ((icr==2) && (ocr==1)) ocr = 2; 
          if (icr==2) pix = 0;
        } else {
          pix = 0;
        }
        break;
      case 13:
        average_menu(&average, &scaleSam, &scaleLin, &filter_type); 
        break;
      case 14:
        if(detect && (oPrmType == deg || oPrmType == rad)) {
          minPix = 0;
          maxPix = (oPrmType == deg) ? 360 : 2*PI;
          pix_phase_menu(&pix, &transType, &otype, &osize, &ocr, 
                         &minPix, &maxPix, &minOut, &maxOut, &fixPix);
        } else {
          if(!detect) oPrmType = iPrmType;
          pix_menu(&pix, &transType, &otype, &osize, &ocr, &minPix, 
                   &maxPix, &powerPix, &fixPix, oPrmType,
                   &minOut, &maxOut);
        }
        if (pix && icr == 2) detect = 1; 
        break;
      case 15:
        menu_ok = 1; 
        printf("\n"); 
        break;
      case 0:
        printf("\n"); 
        goto final_exit;
      default: 
        ;
      } 
    } while (!menu_ok); /* Menu loop. */
    menu_ok = 0;

    /*************************************************************************/
    /* Section 2: Prepare for processing.                                    */
    /*************************************************************************/

    /*-----------------------------------------------------------------------*/
    /* Prepare files.                                                        */
    /*-----------------------------------------------------------------------*/
  
    if (!strcmp(iNam,oNam)) { 
      EXIT(E_SPEC, "Input and output files must differ");
    } else {
      if ((iFp=fopen(iNam,"rb")) == NULL) EXIT(E_FOPE,iNam);
      if ((oFp=fopen(oNam,"wb")) == NULL) EXIT(E_FOPE,oNam);
    }
  
    /*-----------------------------------------------------------------------*/
    /* Prepare window, types etc.                                            */
    /*-----------------------------------------------------------------------*/
  
    if (!win) {
      uSam = iSam;
      uLin = iLin;
      saOf = 0;
      liOf = 0;
    }

    if (!detect) oPrmType = iPrmType;
    
    if (detect && iPrmType==amp && oPrmType!=rad && oPrmType!=deg) 
      icrAv = 1;
    else 
      icrAv = icr; 

    /*-----------------------------------------------------------------------*/
    /* Prepare for averaging (and set window in case of no average).         */
    /*-----------------------------------------------------------------------*/

    if (average) { 
      oLin = (SINT)(uLin*scaleLin); 
      oSam = (SINT)(uSam*scaleSam);
 
      /* Filter lengths. */
      f_length(filter_type,scaleSam,&isaflen,&osaflen);
      f_length(filter_type,scaleLin,&iliflen,&oliflen);
     
      saklen = OVER_SAMP*isaflen;
      liklen = OVER_SAMP*iliflen;
    
      if ((arr =(DREAL *)malloc(sizeof(DREAL)*icrAv*oSam*iliflen)) == NULL) 
        EXIT(E_ALLO,"arr");
      if ((safilt =(DREAL *)malloc(sizeof(DREAL)*saklen)) == NULL)
        EXIT(E_ALLO,"safilt");
      if ((lifilt =(DREAL *)malloc(sizeof(DREAL)*liklen)) == NULL) 
        EXIT(E_ALLO,"lifilt");
    
      /* Initialize filter look-up tables. */
      switch(filter_type) {
      case SINC: 
        f_init(safilt, scaleSam, OVER_SAMP, saklen, filter_type);
        f_init(lifilt, scaleLin, OVER_SAMP, liklen, filter_type);
        break;
      case RCOS: 
        f_init(safilt, 2.0/isaflen, OVER_SAMP, saklen, filter_type);
        f_init(lifilt, 2.0/iliflen, OVER_SAMP, liklen, filter_type);
        break;
      case IDMP: 
        saAv = isaflen/2;
        liAv = iliflen/2;
        break;  
      }
    } /* average */

    if (!average) {
      oLin = uLin; 
      oSam = uSam; 
      scaleSam = 1; 
      scaleLin = 1;
    }

    /*-----------------------------------------------------------------------*/
    /* Allocate memory.                                                      */
    /*-----------------------------------------------------------------------*/
  
    if ((ibuf=(UBYTE *)malloc(isize*icr*uSam)) == NULL) EXIT(E_ALLO,"ibuf");
    if ((obuf=(UBYTE *)malloc(osize*ocr*oSam)) == NULL) EXIT(E_ALLO,"obuf");
    if ((iarr=(DREAL *)malloc(2*sizeof(DREAL)*uSam))==NULL) EXIT(E_ALLO,"iarr");
    if ((oarr=(DREAL *)malloc(2*sizeof(DREAL)*oSam))==NULL) EXIT(E_ALLO,"oarr");

    /*************************************************************************/
    /* Section 3: Here comes the processing loop.                            */
    /*************************************************************************/

    iline = -1;
    for (oline = 0; oline < oLin; oline++) {

      /*---------------------------------------------------------------------*/
      /* While necessary: Read, detect & average (sample direction) a line.  */
      /*---------------------------------------------------------------------*/
      while ( (li_offset = oline/scaleLin) >= iline-iliflen/2+1 
              && iline+1 < uLin) {

        iline++;

        fPos = ((liOf+iline)*iSam+saOf)*isize*icr;
        if (fseek(iFp, fPos, SEEK_SET) == -1) EXIT(E_FPOS,"prior to read");
        if (fread(ibuf, isize, icr*uSam, iFp) != icr*uSam) EXIT(E_FREA,iNam);
        read_line(iarr, ibuf, 1, itype, isize, uSam, icr, iFormat, inputScale);

        if (detect && iPrmType==amp && oPrmType!=rad && oPrmType!=deg) {
          if (icr==2) float_array_power(iarr, iarr, uSam); 
          else        float_array_power_real(iarr, uSam);
          if (oPrmType == amp) float_array_sqrt(iarr, uSam);
        }
      
        if (average) {
          arrp = arr+icrAv*oSam*(iline%iliflen);
          for (samp = 0; samp < (icrAv*oSam); samp++) arrp[samp] = 0.0;  

          if (filter_type==SINC || filter_type==RCOS) {
            #pragma parallel local(samp,sa_offset,sa_int,sa_frac,iarrp,filtp) 
            #pragma byvalue(osaflen,oSam,scaleSam,iarr,icrAv,isaflen,safilt,saklen)
            #pragma shared(arrp) 
            {
              #pragma pfor iterate (samp = osaflen/2; oSam - osaflen; 1)
              for (samp   = osaflen/2; samp < oSam - osaflen/2; samp++) {
                sa_offset = (DREAL)samp/scaleSam;
                sa_int    = (SINT)sa_offset;
                sa_frac   = sa_offset-sa_int;
                iarrp     = iarr+icrAv*(sa_int+isaflen/2);
                if (icrAv == 1) {
                  for (filtp = safilt+(SINT)(sa_frac*OVER_SAMP+0.5);  
                       filtp < safilt+saklen; filtp += OVER_SAMP) 
                    arrp[samp]  += *filtp*(*iarrp--);  
                } else { 
                  iarrp++;
                  for (filtp = safilt+(SINT)(sa_frac*OVER_SAMP+0.5);        
                       filtp < safilt+saklen; filtp += OVER_SAMP) {
                    arrp[2*samp+1] += *filtp*(*iarrp--); 
                    arrp[2*samp  ] += *filtp*(*iarrp--);    
                  }
                }
              }
            }
          } else { /* IDMP */
            if (icrAv==1) {
              #pragma parallel local(samp,n) byvalue(oSam,saAv) shared(iarr,arrp)
              {
                #pragma pfor iterate (samp=0; oSam; 1)
                for (samp=0; samp<oSam; samp++)
                  for (n=0; n<saAv; n++)
                    arrp[samp] += iarr[samp*saAv+n];
              }
            } else {
              #pragma parallel local(samp,n) byvalue(oSam,saAv) shared(iarr,arrp)
              {
                #pragma pfor iterate (samp=0; oSam; 1)
                for (samp=0; samp<oSam; samp++)
                  for (n=0; n<saAv; n++) {
                    arrp[2*samp  ] += iarr[2*samp*saAv+2*n  ];
                    arrp[2*samp+1] += iarr[2*samp*saAv+2*n+1];
                  }
              }
            }
          } /* filter_type */
        } /* average */
      } /* end while */ 
      
      /*---------------------------------------------------------------------*/
      /* Average (line direction) or copy                                    */
      /*---------------------------------------------------------------------*/

      if (average) {
        for (samp = 0; samp < (icrAv*oSam); samp++) oarr[samp] = 0.0;   

        if (filter_type==SINC || filter_type==RCOS) {
          if (oline >= oliflen/2 && oline < oLin-oliflen/2) {  
            k = iline;
            li_int    = (SINT)li_offset;
            li_frac   = li_offset-li_int;  
            for (filtp = lifilt+(SINT)(li_frac*OVER_SAMP+0.5);  
                 filtp < lifilt+liklen; filtp += OVER_SAMP)  {
              arrp = arr+icrAv*oSam*(k%iliflen);
              #pragma parallel local(samp) byvalue(arrp, filtp, icrAv, oSam)
              #pragma shared(oarr)
              {
                #pragma pfor iterate (samp = 0; icrAv*oSam; 1)
                for (samp = 0; samp < icrAv*oSam; samp++) 
                  oarr[samp] += *filtp*arrp[samp];
              }
              k += iliflen-1;  /* decrease k by one modulo iliflen */
            }
          }
        } else { /* IDMP */
          for (n=0; n<liAv; n++) {
            arrp = arr+icrAv*oSam*((oline*liAv+n)%iliflen);
            #pragma parallel local(samp) byvalue(icrAv,oSam) shared(oarr,arrp)
            {
              #pragma pfor iterate (samp=0; icrAv*oSam; 1)
              for (samp=0; samp<icrAv*oSam; samp++)
                oarr[samp] += arrp[samp];
            }
            for (samp=0; samp<icrAv*oSam; samp++) oarr[samp] *= 1.0/(saAv*liAv);
          }
        } /* filter_type */
      } /* average */

      if (!average) { 
        for (i=0; i<icrAv*oSam; i++) oarr[i] = iarr[i];
      }

      /*---------------------------------------------------------------------*/
      /* Remaining transformations (limit, detect etc.).                     */
      /*---------------------------------------------------------------------*/
      if (average && icrAv==1 && (iPrmType==pwr || iPrmType==amp)) { 
        for (i=0; i<oSam; i++) oarr[i] = MAX(0,oarr[i]);
      }

      if (detect) {
        if (oPrmType==rad || oPrmType==deg) {
          for (i=0; i<oSam; i++) {
            oarr[i] = ATAN2(oarr[2*i+1],oarr[2*i]);
            if (oarr[i] < 0) oarr[i] += 2*PI;
          }
          if (oPrmType==deg) for (i=0; i<oSam; i++) oarr[i] *= 180/PI;
        } else if (iPrmType==pwr) {
          if (icrAv == 2) {
            float_array_power(oarr, oarr, oSam);
            float_array_sqrt(oarr, oSam);
          }
          if (oPrmType==amp)float_array_sqrt(oarr, oSam);
        }
      }
    
      if (pix) {
        switch(transType) {
        case LOG10:
          log_transform(oarr,oSam,oPrmType,minPix);
          linear_transform(oarr,oSam,minPix,maxPix,minOut,maxOut);
          break;                   
        case POWER:
          root_transform(oarr,oSam,maxPix*fixPix/maxOut,powerPix,minPix);
          linear_transform(oarr,oSam,minPix,maxPix,minOut,maxOut);
          break;
        case LIN:
          linear_transform(oarr,oSam,minPix,maxPix,minOut,maxOut);
          break;
        default: 
          ;
        }                    
      }
      
      if (ocr == 2 && icrAv == 1) {
        for (i = oSam-1; i >= 0; i--) {
          oarr[i*ocr  ] = oarr[i];
          oarr[i*ocr+1] = 0;
        }
      }
                    
      /*---------------------------------------------------------------------*/
      /* Report progress and write to file.                                  */
      /*---------------------------------------------------------------------*/
      if (oline%100==0) fprintf(stderr,"Writing line %5li\r", (LINT)oline);

      if (!pix) nice_saturation(oarr,otype,oSam*ocr,&uSat,&lSat);
      write_line(oarr, obuf, 1, otype, osize, oSam, ocr, oFormat);
      fPos = oline*oSam*osize*ocr;     
      if (fseek(oFp, fPos, SEEK_SET) == -1) EXIT(E_FPOS,"prior to write");
      if (fwrite(obuf, osize, ocr*oSam, oFp) != ocr*oSam)
        EXIT(E_FWRI,oNam);
    }

    /*************************************************************************/
    /* Section 4: Processing of image completed - write status and cleanup.  */
    /*************************************************************************/
 
    printf("\nThere are %li samples per line and %li lines in output file\n", 
           (LINT)oSam, (LINT)oLin);
  
    if (uSat > 0 && !pix)
      printf("\n%10.5lf percent of values are saturated upwards\n", 
             100*((DREAL)uSat/(oSam*oLin*ocr)), uSat);
    if (lSat > 0 && !pix) 
      printf("\n%10.5lf percent of values are saturated downwards\n", 
             100*((DREAL)lSat/(oSam*oLin*ocr)), lSat);
  
nice_exit:
    uSat = 0;
    lSat = 0;
  
    if (ibuf != NULL)  free((char *)ibuf);                         
    if (obuf != NULL)  free((char *)obuf);                         
    if (iarr != NULL)  free((char *)iarr);                         
    if (oarr != NULL)  free((char *)oarr);                         
  
    if (arr != NULL)     free((char *)arr);                         
    if (safilt != NULL)  free((char *)safilt); 
    if (lifilt != NULL)  free((char *)lifilt); 
                          
    fclose(iFp);
    fclose(oFp);

    if (rval == -1) {
      printf("\n  -- Press return to continue -- ");
      fflush(stdin);
      scanf("%c", &ctmp);
      fflush(stdin);
      rval = 0;
    }
  } while(1);

  /***************************************************************************/
  /* Exiting from main menu loop. Final cleanup.                             */
  /***************************************************************************/
  
final_exit:
  if (ibuf != NULL)  free((char *)ibuf);                         
  if (obuf != NULL)  free((char *)obuf);                         
  if (iarr != NULL)  free((char *)iarr);                         
  if (oarr != NULL)  free((char *)oarr);                         
  
  if (arr != NULL)     free((char *)arr);                         
  if (safilt != NULL)  free((char *)safilt); 
  if (lifilt != NULL)  free((char *)lifilt); 
                          
  fclose(iFp);
  fclose(oFp);

  printf("\n");

  return(EXIT_SUCCESS);

} /* end main */
 
/*****************************************************************************/
/* End MAIN.                                                                 */
/*****************************************************************************/
void set_fixpix(DREAL *fixPix, enum data_type otype)
{
  switch(otype) {
    case ubyte: *fixPix = 128; break;  
    case sbyte: *fixPix = 64; break;
    case uword: *fixPix = 32768; break;
    case sword: *fixPix = 16384; break;
    case uInt:  *fixPix = 2147483648; break;
    case sInt:  *fixPix = 1073741824; break;
  }
}    
/******************************************************************************/
void log_transform(DREAL *arr, SINT len, enum parm_type oPrmType, DREAL minPix) 
{
  SINT  count;
  SINT a = 20;

  if(oPrmType == pwr) a = 10;

  #pragma parallel local(count) byvalue(len, a) shared(arr)
  {
    #pragma pfor iterate (count = 0; len; 1)
    for (count = 0; count < len; count++) {
      if (arr[count] > 0)
        arr[count] = a*log10(arr[count]); 
      else 
        arr[count] = minPix;
    }
  }                                                
}
/******************************************************************************/
void root_transform(DREAL *arr, SINT len, DREAL pivot, DREAL root, DREAL minPix) 
{
  SINT  count;
  
  #pragma parallel local(count) byvalue(len, pivot, root) shared(arr)
  {
    #pragma pfor iterate (count = 0; len; 1)
    for (count = 0; count < len; count++) {
      if (arr[count] > 0)
        arr[count] = pow( (arr[count])/pivot, root)*pivot; 
      else 
        arr[count] = minPix;
    }                                                
  }
}
/******************************************************************************/
void linear_transform(DREAL *arr, SINT len, DREAL minPix, DREAL maxPix, 
                      DREAL minOut, DREAL maxOut) 
{
  SINT  count;
  DREAL a, b;

  a = (maxOut - minOut) / (maxPix - minPix);
  b = minOut - a * minPix;

  #pragma parallel local(count) byvalue(len, a, b, maxOut, minOut) shared(arr)
  {
    #pragma pfor iterate (count = 0; len; 1)
    for (count = 0; count < len; count++) {
      arr[count] = arr[count] * a + b; 
      if (arr[count] > maxOut) arr[count] = maxOut;
      if (arr[count] < minOut) arr[count] = minOut;
    }
  }
}
/******************************************************************************/
void set_min_max_out(DREAL *omin, DREAL *omax, enum data_type dtype)
{
  switch(dtype) {
    case ubyte: *omax = 255.0; *omin = 0.0; break;
    case sbyte: *omax = 127.0; *omin = 0.0; break;
    case uword: *omax = 65535.0; *omin = 0.0; break;
    case sword: *omax = 32767.0; *omin = 0.0; break;
    case uInt:  *omax = 4294967295.0; *omin = 0.0; break;
    case sInt:  *omax = 2147483647.0; *omin = 0.0; break;
  } 
}     
/******************************************************************************/
void find_type_min_max(DREAL *typeMin, DREAL *typeMax, enum data_type dtype)
{
  switch(dtype) {
    case ubyte: *typeMax = 255.0; *typeMin = 0.0; break;
    case sbyte: *typeMax = 127.0; *typeMin = -128.0; break;
    case uword: *typeMax = 65535.0; *typeMin = 0.0; break;
    case sword: *typeMax = 32767.0; *typeMin = -32768.0; break;
    case uInt:  *typeMax = 4294967295.0; *typeMin = 0.0; break;
    case sInt:  *typeMax = 2147483647.0; *typeMin = -2147483648.0; break;
  }
}
/******************************************************************************/
/*Can not be parallelized because of data dependencies because
  iarr and oarr might be the same array
DREAL float_array_power(DREAL *iArr, DREAL *oArr, SINT len)
{
  DREAL maxAmp=0.0;
  SINT  count;

  #pragma parallel local(count) byvalue(len, iArr) shared(oArr)
  {
    #pragma pfor iterate (count = 0; len; 1)
    for (count = 0; count < len; count++)
      oArr[count] = SQUARE(iArr[count*2])+SQUARE(iArr[count*2+1]); 
  }
  return(maxAmp);
}
*/
void float_array_power(DREAL *iArr, DREAL *oArr, SINT len)
{
  SINT  count;

  for (count = 0; count < len; count++) {
    *oArr = SQUARE(*iArr)+SQUARE(*(iArr+1));
    iArr += 2;
    oArr++;
  }
}
/******************************************************************************/
void f_length(FILTER_TYPE type, DREAL scale, SINT *iflen, SINT *oflen)
{
  if (scale < 1.0) {
    switch(type) {
    case SINC:
      *iflen = 2*(SINT)(MIN_F_LEN/(2*scale));
      *oflen = 2*(SINT)(MIN_F_LEN/2); 
      break;
    case RCOS:
      *iflen = 2*(SINT)(1.0/scale);
      *oflen = 2;
      break;
    case IDMP:
      *iflen = 2*ROUND(1.0/scale);
      *oflen = 2;
      break;
    }
  } else {
    switch(type) {
    case SINC:
      *iflen = 2*(SINT)(MIN_F_LEN/2);
      *oflen = 2*(SINT)(MIN_F_LEN*scale/2);
      break;
    case RCOS:
      *iflen = 2;
      *oflen = 2*(SINT)(scale);
      break;
    case IDMP: 
      ; /* Expansion not implemented. */
    }
  }
}
/******************************************************************************/
void f_init(DREAL *filter, DREAL scale, SINT over_samp, SINT klen, 
            FILTER_TYPE filter_type)
{
  SINT  samp;
  DREAL *kernel, aux;
                      
  if (scale > 1) scale = 1;
  kernel = filter;
  for (samp = -klen/2; samp < klen/2; samp++) {
    if (filter_type == SINC) {
      aux = PI*scale*samp/over_samp;
      *kernel    = ( samp != 0 )            ? (sin(aux)/aux) *
                   (1 + cos(2*PI*samp/klen))/2 : 1.0;/**/
      *kernel++ *= scale;
    } else if (filter_type == RCOS) {
     *kernel    = (1 + cos(2*PI*samp/klen))/2;
     *kernel++ *= scale;
    }
  }
}
/******************************************************************************/
void float_array_sqrt(DREAL *arr, SINT len)
{
  SINT  count;
   
  #pragma parallel local(count) byvalue(len) shared(arr)
  {
    #pragma pfor iterate (count = 0; len; 1)
    for (count = 0; count < len; count++) 
      arr[count] = sqrt(arr[count]);
  }
}
/******************************************************************************/
void float_array_power_real(DREAL *arr, SINT len)
{
  SINT  count;

  #pragma parallel local(count) byvalue(len) shared(arr)
  {
    #pragma pfor iterate (count = 0; len; 1)
    for (count = 0; count < len; count++)
      arr[count] *= arr[count];
  }                                                
}
/******************************************************************************/
char *parm_type_str(enum parm_type type)
{
  switch(type) {
    case amp: return("AMPLITUDE");
    case pwr: return("POWER");
    case rad: return("PHASE (rad)");
    case deg: return("PHASE (deg)");
    default : return("No such parameter type");
  }
}
/******************************************************************************/
char *on_off(SINT boolean)
{
  if (boolean) return("ON");
  else         return("OFF");
}
/******************************************************************************/
void update_pix_otype_menu(void)
{
  printf("\n\n       ---------- Data type menu ---------- \n\n");
  printf ("\n (1) 8  bit unsigned integer\n");
  printf ("\n (2) 8  bit signed integer\n");
  printf ("\n (3) 16 bit unsigned integer\n");
  printf ("\n (4) 16 bit signed integer\n");
  printf ("\n (5) 32 bit unsigned integer\n");
  printf ("\n (6) 32 bit signed integer\n\n");
}
/******************************************************************************/
void pix_otype_menu(enum data_type *otype, SINT *isize, DREAL *fixPix, 
                    DREAL *minOut, DREAL *maxOut)
{
  CHCK_INIT("pix_otype_menu")
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 6;
  char chline[16];
  
  do {
    update_pix_otype_menu();
    printf("\n");
    do {
      printf("Choose (1) to (%i) for setting output data type: ", 
             (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));

    switch(choice) {
      case 1: menu_ok = 1; *otype = ubyte; *isize = 1; break;
      case 2: menu_ok = 1; *otype = sbyte; *isize = 1; break;
      case 3: menu_ok = 1; *otype = uword; *isize = 2; break;
      case 4: menu_ok = 1; *otype = sword; *isize = 2; break;
      case 5: menu_ok = 1; *otype = uInt;  *isize = 4; break;
      case 6: menu_ok = 1; *otype = sInt;  *isize = 4; break;
    }
  } while (!menu_ok);
  
  set_min_max_out(minOut,maxOut,*otype); 
  set_fixpix(fixPix, *otype); 
}
/******************************************************************************/
void log_menu(DREAL *minPix, DREAL *maxPix, DREAL *minOut, DREAL *maxOut,
              enum data_type otype)
{
  CHCK_INIT("log_menu")
  DREAL inMax = *maxPix, dynRan = *maxPix - *minPix;
  DREAL typeMin, typeMax;
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 4;
  char chline[16];
  
  find_type_min_max(&typeMin, &typeMax, otype);

  do {
    printf("\n\n       ---------- Log transform menu ---------- \n\n");
    printf("\n (1) Maximum value in dB  : %-14.3lf\n", inMax);
    printf("\n (2) Dynamic range in dB  : %-14.3lf\n", dynRan);
    printf("\n (3) Minimum output value : %-14.0lf\n", *minOut);
    printf("\n (4) Maximum output value : %-14.0lf\n", *maxOut);
    printf("\n (0) Exit\n\n");
    
    do {
      printf("Choose (1) to (%i) for changing or (0) for exit: ", 
             (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
    case 1: 
      printf("\nMaximum value in dB : ");
      if (scanf("%lf", &inMax) == EOF) EXIT2(E_SCAN, "inMax"); 
      break;
    case 2: 
      printf("\nDynamic range in dB : ");
      if (scanf("%lf", &dynRan) == EOF) EXIT2(E_SCAN, "dynRan"); 
      if (dynRan <= 0) dynRan = 36;
      break;
    case 3: 
      printf("\nMinimum output value : ");
      if (scanf("%lf", minOut) == EOF) EXIT2(E_SCAN, "minOut"); 
      if (*minOut > *maxOut) *minOut = *maxOut;
      if (*minOut < typeMin) *minOut = typeMin;
      *minOut = round(*minOut);
      break;
    case 4: 
      printf("\nMaximum output value : ");
      if (scanf("%lf", maxOut) == EOF) EXIT2(E_SCAN, "maxOut"); 
      if (*maxOut < *minOut) *maxOut = *minOut;
      if (*maxOut > typeMax) *maxOut = typeMax;
      *maxOut = round(*maxOut);
      break;
    case 0:
      menu_ok = 1;
      break;
    }
  } while (!menu_ok);
  *maxPix = inMax;
  *minPix = inMax - dynRan;
}
/******************************************************************************/
void power_menu(DREAL *minPix, DREAL *maxPix, DREAL *powerPix, DREAL *fixPix,
                enum parm_type oPrmType, DREAL *minOut, DREAL *maxOut,
                enum data_type otype)
{
  CHCK_INIT("power_menu")
  DREAL inMax;
  DREAL typeMin, typeMax;
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 5;
  char chline[16];
  SINT a = 20;

  find_type_min_max(&typeMin, &typeMax, otype);

  if (oPrmType == pwr) a = 10;
  
  if (*maxPix <= 0) *maxPix = 1;
  inMax = a*log10(*maxPix);
  
  *minOut = 0;
  do {
    printf("\n\n       ---------- Power transform menu ---------- \n\n");
    printf("\n (1) Maximum value in dB  : %-14.3lf\n", inMax);
    printf("\n (2) Power                : %-14.3lf\n", *powerPix);
    printf("\n (3) Fixpoint             : %-14.3lf\n", *fixPix);
    printf("\n (4) Minimum output value : %-14.0lf\n", *minOut);
    printf("\n (5) Maximum output value : %-14.0lf\n", *maxOut);
    printf("\n (0) Exit\n\n");

    do {
      printf("Choose (1) to (%i) for changing or (0) for exit: ", 
             (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
    case 1: 
      printf("\nMaximum value in dB : ");
      if (scanf("%lf", &inMax) == EOF) EXIT2(E_SCAN, "inMax");
      break;
    case 2: 
      printf("\nPower : ");
      if (scanf("%lf", powerPix) == EOF) EXIT2(E_SCAN, "powerPix"); 
      if (*powerPix <= 0) *powerPix = 1;
      break;
    case 3: 
      printf("\nFixpoint : ");
      if (scanf("%lf", fixPix) == EOF) EXIT2(E_SCAN, "fixPix"); 
      if (*fixPix <= 0) set_fixpix(fixPix, otype);
      break;
    case 4: 
      printf("\nMinimum output value : ");
      if (scanf("%lf", minOut) == EOF) EXIT2(E_SCAN, "minOut"); 
      WARN(E_SPEC, "Minimum output value must be equal to 0 for this transformation ");
      *minOut = 0;
      break;
    case 5: 
      printf("\nMaximum output value : ");
      if (scanf("%lf", maxOut) == EOF) EXIT2(E_SCAN, "maxOut"); 
      if (*maxOut < *minOut) *maxOut = *minOut;
      if (*maxOut > typeMax) *maxOut = typeMax;
      *maxOut = round(*maxOut);
      break;
    case 0:
      menu_ok = 1; 
      break;
    }
  } while (!menu_ok);
  *maxPix = pow(10, inMax/a);
  *minPix = 0;
}  
/******************************************************************************/
void linear_menu(DREAL *minPix, DREAL *maxPix, DREAL *minOut, DREAL *maxOut,
                 enum data_type otype)
{
  CHCK_INIT("linear_menu")
  DREAL typeMin, typeMax;
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 4;
  char chline[16];

  find_type_min_max(&typeMin, &typeMax, otype);

  do {
    printf("\n\n       ---------- Linear transform menu ---------- \n\n");
    printf("\n (1) Minimum value        : %-14.3lf\n", *minPix);
    printf("\n (2) Maximum value        : %-14.3lf\n", *maxPix);
    printf("\n (3) Minimum output value : %-14.0lf\n", *minOut);
    printf("\n (4) Maximum output value : %-14.0lf\n", *maxOut);
    printf("\n (0) Exit\n\n");

    do {
      printf("Choose (1) to (%i) for changing or (0) for exit: ", 
             (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice"); 
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
    case 1: 
      printf("\nMinimum value : ");
      if (scanf("%lf", minPix) == EOF) EXIT2(E_SCAN, "minPix"); 
      break;
    case 2: 
      printf("\nMaximum value : ");
      if (scanf("%lf", maxPix) == EOF) EXIT2(E_SCAN, "maxPix"); 
      break;
    case 3: 
      printf("\nMinimum output value : ");
      if (scanf("%lf", minOut) == EOF) EXIT2(E_SCAN, "minOut"); 
      if (*minOut > *maxOut) *minOut = *maxOut;
      if (*minOut < typeMin) *minOut = typeMin;
      *minOut = round(*minOut);
      break;
    case 4: 
      printf("\nMaximum output value : ");
      if (scanf("%lf", maxOut) == EOF) EXIT2(E_SCAN, "maxOut"); 
      if (*maxOut < *minOut) *maxOut = *minOut;
      if (*maxOut > typeMax) *maxOut = typeMax;
      *maxOut = round(*maxOut);
      break;
    case 0: 
      menu_ok = 1; 
      break;
    }
  } while (!menu_ok);
}
/******************************************************************************/
void pix_phase_menu(SINT *pix, TRANSFORMATION_TYPE *transType, enum data_type *otype, 
                    SINT *osize, SINT *ocr, DREAL *minPix, DREAL *maxPix,
                    DREAL *minOut, DREAL *maxOut, DREAL *fixPix)
{
  CHCK_INIT("pix_phase_menu")
  SINT choice;
  SINT max_choice = 2;
  char chline[16];
  
  *pix = 1;
  
  printf("\n\n       ---------- Pixel transform menu ---------- \n\n");
  printf("\n (1)  Transform to pixels : %s\n", on_off(*pix));
  printf("\n (2)  linear\n\n");

  do {
    printf("choose (1) to (%i) for setting transformation type: ", 
           max_choice); 
    if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
    choice = -1;
    SSCANF_D(chline,&choice);
  } while ((choice < 0) || (choice > max_choice));
    
  switch(choice) {
  case 1: 
    if (*pix) *pix = 0; 
    else      *pix = 1; 
    break;
  case 2: 
    if (*otype > 5 || *ocr == 2) { 
      pix_otype_menu(otype, osize, fixPix, minOut, maxOut); 
      *ocr = 1; 
    }
    *pix = 1;
    *transType = LIN;
    linear_menu(minPix, maxPix, minOut, maxOut, *otype); 
    break;
  }
}
/******************************************************************************/
void pix_menu(SINT *pix, TRANSFORMATION_TYPE *transType, enum data_type *otype, 
              SINT *osize, SINT *ocr, DREAL *minPix, DREAL *maxPix, 
              DREAL *powerPix, DREAL *fixPix, enum parm_type oPrmType,
              DREAL *minOut, DREAL *maxOut)
{
  CHCK_INIT("pix_menu")
  SINT choice;
  SINT max_choice = 4;
  char chline[16];
  
  *pix = 1;
  
  printf("\n\n       ---------- Pixel transform menu ---------- \n\n");
  printf("\n (1)  Transform to pixels : %s\n", on_off(*pix));
  printf("\n (2)  log10\n");
  printf("\n (3)  power\n");
  printf("\n (4)  linear\n\n");

  do {
    printf("choose (2) to (%i) for setting transformation type: ", 
           max_choice); 
    if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
    choice = -1;
    SSCANF_D(chline,&choice);
  } while ((choice < 0) || (choice > max_choice));
    
  switch(choice) {
  case 1: 
    if (*pix) *pix = 0; 
    else      *pix = 1; 
    break;
  case 2: 
    if (*otype > 5 || *ocr == 2) { 
      pix_otype_menu(otype, osize, fixPix, minOut, maxOut); 
      *ocr = 1;
    }
    *pix = 1;
    if (*transType != LOG10) {
      *minPix = -36; 
      *maxPix = 0;
    }
    *transType = LOG10;
    log_menu(minPix, maxPix, minOut, maxOut, *otype); 
    break;
  case 3: 
    if (*otype > 5 || *ocr == 2) { 
      pix_otype_menu(otype, osize, fixPix, minOut, maxOut); 
      *ocr = 1; 
    }
    *pix = 1;
    if (*transType != POWER) *maxPix = 1; 
    *transType = POWER;
    power_menu(minPix, maxPix, powerPix, fixPix, oPrmType, minOut, maxOut, *otype); 
    break;
  case 4: 
    if (*otype > 5 || *ocr == 2) { 
      pix_otype_menu(otype, osize, fixPix, minOut, maxOut); 
      *ocr = 1; 
    }
    *pix = 1;
    *transType = LIN;
    linear_menu(minPix, maxPix, minOut, maxOut, *otype); break;
  }
}
/******************************************************************************/
void out_parm_menu(enum parm_type *type)
{
  CHCK_INIT("out_parm_menu")
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 4;
  char chline[16];
  
  do {
    printf("\n\n       ---------- Output parameter type menu ---------- \n\n");
    printf ("\n (1)  AMPLITUDE");
    printf ("\n (2)  POWER");
    printf ("\n (3)  PHASE (rad)");
    printf ("\n (4)  PHASE (deg)");
    printf ("\n\n");
    
    do {
      printf("Choose (1) to (%i) for setting data type: ", (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
      case 1: menu_ok = 1; *type = amp; break;
      case 2: menu_ok = 1; *type = pwr; break;
      case 3: menu_ok = 1; *type = rad; break;
      case 4: menu_ok = 1; *type = deg; break;
    }
  } while(!menu_ok);
}
/******************************************************************************/
void detect_menu(SINT *detect, enum parm_type *oPrmType)
{
  CHCK_INIT("detect_menu")
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 2;
  char chline[16];
  
  *detect = 1;
  do {
    printf("\n\n       ---------- Detection menu ---------- \n\n");
    printf("\n (1)  Detection             : %s\n", on_off(*detect));
    printf("\n (2)  Output parameter type : %s\n", parm_type_str(*oPrmType));
    printf("\n (0)  Exit\n\n");

    do {
      printf("Choose (1) to (%i) for changing or (0) for exit: ", 
             (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
    case 1: 
      if (*detect) *detect = 0; 
      else         *detect = 1; 
      break;
    case 2: 
      out_parm_menu(oPrmType); 
      break;
    case 0:
      menu_ok = 1; 
      break;
    }
  } while (!menu_ok);
}
/******************************************************************************/
CHAR *filter_name(FILTER_TYPE filter_type)
{
  switch(filter_type) {
    case SINC : return("WEIGHTED SINC");
    case RCOS : return("RAISED COSINE");
    case IDMP : return("INTEGRATE DUMP");
    default : return("No such filter type");
  }
}
/******************************************************************************/
void average_menu(SINT *average, DREAL *scaleSam, DREAL *scaleLin, 
                  FILTER_TYPE *filter_type)
{
  CHCK_INIT("average_menu")
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 4;
  char chline[16];
  
  *average = 1;
  do {
    printf("\n\n       ---------- Average menu ---------- \n\n");
    printf("\n (1)  Averaging          : %s\n", on_off(*average));
    printf("\n (2)  Sample scale factor: %-14.3lf\n", *scaleSam);
    printf("\n (3)  Line scale factor  : %-14.3lf\n", *scaleLin);
    printf("\n (4)  Filter type        : %s\n", filter_name(*filter_type));
    printf("\n (0)  Exit\n\n");

    do {
      printf("Choose (1) to (%i) for changing or (0) for exit: ", 
             (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
    case 1: 
      if (*average) *average = 0; 
      else          *average = 1; 
      break;
    case 2: 
      printf("\nSample scale factor: ");
      if (scanf("%lf", scaleSam) == EOF) EXIT2(E_SCAN, "scaleSam");
      if (*scaleSam <= 0) *scaleSam = 1; 
      if (*filter_type==IDMP)
        *scaleSam = 1.0/MAX(1,ROUND(1.0/(*scaleSam)));
      break;
    case 3: 
      printf("\nLine scale factor: ");
      if (scanf("%lf", scaleLin) == EOF) EXIT2(E_SCAN, "scaleLin"); 
      if (*scaleLin <= 0) *scaleLin = 1; 
      if (*filter_type==IDMP)
        *scaleLin = 1.0/MAX(1,ROUND(1.0/(*scaleLin)));
      break;
    case 4: 
      filter_menu(filter_type);
      if (*filter_type==IDMP) {
        *scaleSam = 1.0/MAX(1,ROUND(1.0/(*scaleSam)));
        *scaleLin = 1.0/MAX(1,ROUND(1.0/(*scaleLin)));
      }
      break;
    case 0: 
      menu_ok = 1; 
      break;
    }
  } while (!menu_ok);
}  
/******************************************************************************/
void filter_menu(FILTER_TYPE *filter_type)
{
  CHCK_INIT("filter_menu")
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 3;
  char chline[16];
  
  do {
    printf("\n\n       ---------- Filter type menu ---------- \n\n");
    printf("\n (1)  Weighted sinc \n");
    printf("\n (2)  Raised cosine \n");
    printf("\n (3)  Integrate & dump\n");
    printf("\n (0)  Exit\n\n");

    do {
      printf("Choose (1) to (%i) for changing og (0) for exit: ", 
              (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
    case 0: 
      menu_ok = 1;
      break;
    case 1: 
      *filter_type = SINC; 
      menu_ok = 1;
      break;
    case 2:
      *filter_type = RCOS;
      menu_ok = 1;
      break;
    case 3: 
      *filter_type = IDMP;
      menu_ok = 1;
      break;
    }
  } while (!menu_ok);
}  
/******************************************************************************/
char *name_type(enum data_type type, SINT icr)
{
  if (icr == 1) {
    switch(type) {
      case ubyte: return("8 bit unsigned integer");
      case sbyte: return("8 bit signed integer");
      case uword: return("16 bit unsigned integer"); 
      case sword: return("16 bit signed integer"); 
      case uInt: return("32 bit unsigned integer"); 
      case sInt: return("32 bit signed integer"); 
      case half_sreal: return("16 bit float");
      case sreal: return("32 bit float");
      case dreal: return("64 bit float"); 
    }
  } else {
    switch(type) {
      case ubyte: return("complex 8 bit unsigned integer");
      case sbyte: return("complex 8 bit signed integer");
      case uword: return("complex 16 bit unsigned integer"); 
      case sword: return("complex 16 bit signed integer"); 
      case uInt: return("complex 32 bit unsigned integer"); 
      case sInt: return("complex 32 bit signed integer"); 
      case half_sreal: return("complex 16 bit float");
      case sreal: return("complex 32 bit float");
      case dreal: return("complex 64 bit float"); 
    }
  }
  return("No such datatype");
}  
/******************************************************************************/
void update_menu(char *iNam, enum data_type itype, SINT icr, 
                 SINT iSam, SINT iLin, enum Format iFormat, 
                 DREAL inputScale, char *oNam, enum data_type otype, 
                 SINT ocr, enum Format oFormat, SINT win, SINT detect, 
                 SINT average, SINT pix, enum parm_type iPrmType)
{
  printf("\n       ---------- Main menu ---------- \n");
  printf("\n (1)  Input file name            : %s", iNam);
  printf("\n (2)  Input data type            : %s", name_type(itype, icr));  
  printf("\n (3)  Input parameter type       : %s", parm_type_str(iPrmType));
  printf("\n (4)  Number of samples per line : %li", (LINT)iSam);
  printf("\n (5)  Number of lines            : %li", (LINT)iLin);
  printf("\n (6)  Input image storing format : %s", 
         iFormat == Dos ? "DOS" : "UNIX");
  printf("\n\n (7)  Output file name           : %s", oNam);
  printf("\n (8)  Output data type           : %s" , name_type(otype, ocr));
  printf("\n (9)  Output image storing format: %s", 
         oFormat == Dos ? "DOS" : "UNIX");
  printf("\n\n (10) Input scale factor         : %e", inputScale);
  printf("\n (11) Windowing                  : %s", on_off(win));  
  printf("\n (12) Detection                  : %s", on_off(detect));  
  printf("\n (13) Averaging                  : %s", on_off(average));  
  printf("\n (14) Transform to pixels        : %s", on_off(pix));
  printf("\n (15) Go");
  printf("\n (0)  Exit\n");
}  
/******************************************************************************/
void update_type_menu(void)
{
  printf("\n\n       ---------- Data type menu ---------- \n\n");
  printf ("\n (1)  8  bit unsigned integer");
  printf ("\n (2)  8  bit signed integer");
  printf ("\n (3)  16 bit unsigned integer");
  printf ("\n (4)  16 bit signed integer");
  printf ("\n (5)  32 bit unsigned integer");
  printf ("\n (6)  32 bit signed integer");
  printf ("\n (7)  16 bit float");
  printf ("\n (8)  32 bit float");
  printf ("\n (9)  64 bit float");
  printf ("\n (10) complex 8  bit unsigned integer");
  printf ("\n (11) complex 8  bit signed integer");
  printf ("\n (12) complex 16 bit unsigned integer");
  printf ("\n (13) complex 16 bit signed integer");
  printf ("\n (14) complex 32 bit unsigned integer");
  printf ("\n (15) complex 32 bit signed integer");
  printf ("\n (16) complex 16 bit float");
  printf ("\n (17) complex 32 bit float");
  printf ("\n (18) complex 64 bit float\n\n");
  }
/******************************************************************************/
void type_menu(enum data_type *itype, SINT *icr, SINT *isize)
{
  CHCK_INIT("type_menu")
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 18;
  char chline[16];
  
  do {
    update_type_menu();

    do {
      printf("Choose (1) to (%i) for setting data type: ", (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
      case 1:  menu_ok = 1; *itype = ubyte; *icr = 1; *isize = 1; break;
      case 2:  menu_ok = 1; *itype = sbyte; *icr = 1; *isize = 1; break;
      case 3:  menu_ok = 1; *itype = uword; *icr = 1; *isize = 2; break;
      case 4:  menu_ok = 1; *itype = sword; *icr = 1; *isize = 2; break;
      case 5:  menu_ok = 1; *itype = uInt;  *icr = 1; *isize = 4; break;
      case 6:  menu_ok = 1; *itype = sInt;  *icr = 1; *isize = 4; break;
      case 7:  menu_ok = 1; *itype = half_sreal; *icr = 1; *isize = 2; break;
      case 8:  menu_ok = 1; *itype = sreal; *icr = 1; *isize = 4; break;
      case 9:  menu_ok = 1; *itype = dreal; *icr = 1; *isize = 8; break;
      case 10: menu_ok = 1; *itype = ubyte; *icr = 2; *isize = 1; break;
      case 11: menu_ok = 1; *itype = sbyte; *icr = 2; *isize = 1; break;
      case 12: menu_ok = 1; *itype = uword; *icr = 2; *isize = 2; break;
      case 13: menu_ok = 1; *itype = sword; *icr = 2; *isize = 2; break;
      case 14: menu_ok = 1; *itype = uInt;  *icr = 2; *isize = 4; break;
      case 15: menu_ok = 1; *itype = sInt;  *icr = 2; *isize = 4; break;
      case 16: menu_ok = 1; *itype = half_sreal; *icr = 2; *isize = 2; break;
      case 17: menu_ok = 1; *itype = sreal; *icr = 2; *isize = 4; break;
      case 18: menu_ok = 1; *itype = dreal; *icr = 2; *isize = 8; break;
    }
  } while (!menu_ok);
}
/******************************************************************************/
void update_window_menu(SINT win, SINT uSam, SINT uLin, SINT saOf, SINT liOf)
{
  printf("\n\n       ---------- Window menu ---------- \n\n");
  printf("\n (1)  Windowing                       : %s\n", on_off(win));
  printf("\n (2)  Number of used samples per line : %li\n", (LINT)uSam);
  printf("\n (3)  Number of used lines            : %li\n", (LINT)uLin);
  printf("\n (4)  Sample offset                   : %li\n", (LINT)saOf);
  printf("\n (5)  Line offset                     : %li\n", (LINT)liOf);
  printf("\n (0)  Exit\n\n");
}  
/******************************************************************************/
void window_menu(SINT *win, SINT *uSam, SINT *uLin, SINT *saOf, SINT *liOf, 
                 SINT iSam, SINT iLin)
{
  CHCK_INIT("window_menu")
  SINT menu_ok = 0;
  SINT choice;
  SINT max_choice = 5;
  char chline[16];
  
  *win = 1;
  do {
    update_window_menu(*win, *uSam, *uLin, *saOf, *liOf);

    do {
      printf("Choose (1) to (%i) for changing or (0) for exit: ", 
             (int)max_choice);
      if (scanf("%s",chline) == EOF) EXIT2(E_SCAN, "menu choice");
      choice = -1;
      SSCANF_D(chline,&choice);
    } while ((choice < 0) || (choice > max_choice));
    
    switch(choice) {
    case 1: 
      if (*win) *win = 0;
      else      *win = 1;
      break;
    case 2: 
      printf("\nNumber of used samples per line : ");
      if (scanf("%i", uSam) == EOF) EXIT2(E_SCAN, "uSam"); 
      if (*uSam + *saOf > iSam) *uSam = iSam - *saOf;
      if (*uSam <=0) *uSam = iSam - *saOf; 
      break;
    case 3: 
      printf("\nNumber of used lines            : ");
      if (scanf("%i", uLin) == EOF) EXIT2(E_SCAN, "uLin"); 
      if (*uLin + *liOf > iLin) *uLin = iLin - *liOf;
      if (*uLin <=0) *uLin = iLin - *liOf;
      break;
    case 4: 
      printf("\n Sample offset                   : ");
      if (scanf("%i", saOf) == EOF) EXIT2(E_SCAN, "saOf"); 
      if (*uSam + *saOf > iSam) *saOf = 0;
      if (*saOf < 0) *saOf = 0; 
      break;
    case 5: 
      printf("\n Line offset                   : ");
      if (scanf("%i", liOf) == EOF) EXIT2(E_SCAN, "liOf"); 
      if (*uLin + *liOf > iLin) *liOf = 0;
      if (*liOf < 0) *liOf = 0;
      break;
    case 0: 
      menu_ok = 1; 
      break;
    }
  } while (!menu_ok);
}
/******************************************************************************/
void float16_to_float64(SWORD *sbuf, DREAL *rbuf, SINT length,
                        SINT step, SINT icr, DREAL inputScale)
{
  int  i, first, second;
  RS_U rs;
 
  if(machine_type == Unix) { first = 0; second = 1; }
  else { first = 1; second = 0; }

  #pragma parallel local(i,rs) byvalue(length, icr, inputScale, sbuf, step,first,second) shared(rbuf)
  {
    #pragma pfor iterate (i = 0; length; 1)
    for (i = 0; i < length; i++) {
      rs.s[first] = sbuf[i*step];
      rs.s[second] = 0;
      rbuf[icr*i] = (DREAL) rs.r * inputScale;
      if(icr == 2) {
        rs.s[first] = sbuf[i*step+1];
        rs.s[second] = 0;
        rbuf[icr*i+1] = (DREAL) rs.r * inputScale;
      }
    }
  }
}
/******************************************************************************/
void float64_to_float16(DREAL *rbuf, SWORD *sbuf, SINT length,
                        SINT step, SINT icr)
{
  int  i, first;
  RS_U rs;

  if(machine_type == Unix) { first = 0; }
  else { first = 1; }

  #pragma parallel local(i,rs) byvalue(length, icr,step,first,rbuf) shared(sbuf)
  {
    #pragma pfor iterate (i = 0; length; 1)
    for (i = 0; i < length; i++) {
      rs.r = (SREAL) rbuf[i*icr];
      sbuf[i*step] = rs.s[first];
      if(icr == 2) {
        rs.r = (SREAL) rbuf[i*icr+1];
        sbuf[i*step+1] = rs.s[first];
      }
    }
  }
}
/******************************************************************************/
void swap(UBYTE *ibufp, SINT nSam, SINT icr, SINT step, SINT size)
{
  UBYTE tmp;
  SINT i,j,k;
  for (i = 0; i < nSam; i++, ibufp += step*size) {
    for (k = 0; k < icr; k++) {
      for (j = 0; j < size/2; j++) {
        tmp = *(ibufp+j+k*size);
        *(ibufp+j+k*size) = *(ibufp+size-j-1+k*size);
        *(ibufp+size-j-1+k*size) = tmp;
      }
    }
  }
}
/******************************************************************************/
SINT read_line(DREAL *iarr, UBYTE *ibufp, SINT step, 
               enum data_type itype, SINT size, SINT nSam, 
               SINT icr, enum Format format, DREAL inputScale)
{
  CHCK_INIT("read_line")
  int i;
  step *= icr;

  if ((format != machine_type) && (size > 1)) 
    swap(ibufp, nSam, icr, step, size);
    
  switch(itype) {
  case ubyte: 
    for(i = 0; i < nSam; i++, ibufp += step) {
      *iarr++ = (DREAL)*ibufp * inputScale;
      if(icr == 2) *iarr++ = (DREAL) *(ibufp+1) * inputScale;
    }
    break;
  case sbyte: 
    {
      SBYTE *ptr;
      ptr = (SBYTE*)ibufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *iarr++ = (DREAL)*ptr * inputScale;
        if(icr == 2) *iarr++ = (DREAL) *(ptr+1) * inputScale;
      }
    }
    break;
  case uword: 
    {
      UWORD *ptr;
      ptr = (UWORD*)ibufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *iarr++ = (DREAL)*ptr * inputScale;
        if(icr == 2) *iarr++ = (DREAL) *(ptr+1) * inputScale;
      }
    }
    break;
  case sword: 
    {
      SWORD *ptr;
      ptr = (SWORD*)ibufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *iarr++ = (DREAL)*ptr * inputScale;
        if(icr == 2) *iarr++ = (DREAL) *(ptr+1) * inputScale;
      }
    }
    break;
  case uInt: 
    {
      UINT *ptr;
      ptr = (UINT*)ibufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *iarr++ = (DREAL)*ptr * inputScale;
        if(icr == 2) *iarr++ = (DREAL) *(ptr+1) * inputScale;
      }
    }
    break;
  case sInt:  
    {
      SINT *ptr;
      ptr = (SINT*)ibufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *iarr++ = (DREAL)*ptr * inputScale;
        if(icr == 2) *iarr++ = (DREAL) *(ptr+1) * inputScale;
      }
    }
    break;
  case half_sreal:
    { 
      SWORD *ptr;
      ptr = (SWORD*)ibufp;
      float16_to_float64(ptr, iarr, nSam, step, 
                         icr, inputScale);
    } /* */
    break;
  case sreal:
    {
      SREAL *ptr;
      ptr = (SREAL*)ibufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *iarr++ = (DREAL)*ptr * inputScale;
        if(icr == 2) *iarr++ = (DREAL) *(ptr+1) * inputScale;
      }
    }
    break;
  case dreal:
    {
      DREAL *ptr;
      ptr = (DREAL*)ibufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *iarr++ = (DREAL)*ptr * inputScale;
        if(icr == 2) *iarr++ = (DREAL) *(ptr+1) * inputScale;
      }
    }
    break;
  default:    
    EXIT(E_TYPE, "itype");
  } 
  nice_exit: 
  return(rval);           
}
/******************************************************************************/
SINT write_line(DREAL *oarr, UBYTE *obufp, SINT step, 
                      enum data_type otype, SINT size, SINT nSam, SINT icr,
                      enum Format format)
{
  CHCK_INIT("write_line")
  int i;
  step *= icr;
  switch(otype) {
  case ubyte:
    for(i = 0; i < nSam; i++, obufp += step) {
      *obufp = (UBYTE) (0.5 + *oarr++);
      if(icr == 2) *(obufp+1) = (UBYTE) (0.5 + *oarr++);
    }
    break;
  case sbyte: 
    {
      SBYTE *ptr;
      ptr = (SBYTE*)obufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        if (*oarr >= 0) *ptr = (SBYTE)  (0.5 + *oarr++);
        else            *ptr = (SBYTE) (-0.5 + *oarr++);
        if(icr == 2) { 
          if (*oarr >= 0) *(ptr+1) = (SBYTE)  (0.5 + *oarr++);
          else            *(ptr+1) = (SBYTE) (-0.5 + *oarr++);
        } 
      }
    }
    break;
  case uword:
    {
      UWORD *ptr;
      ptr = (UWORD*)obufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *ptr = (UWORD) (0.5 + *oarr++);
        if(icr == 2) *(ptr+1) = (UWORD) (0.5 + *oarr++);
      }
    }
    break;
  case sword: 
    {
      SWORD *ptr;
      ptr = (SWORD*)obufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        if (*oarr >= 0) *ptr = (SWORD)  (0.5 + *oarr++);
        else            *ptr = (SWORD) (-0.5 + *oarr++);
        if(icr == 2) {
          if (*oarr >= 0) *(ptr+1) = (SWORD)  (0.5 + *oarr++);
          else            *(ptr+1) = (SWORD) (-0.5 + *oarr++);
        }
      }
    }
    break;
  case uInt: 
    {
      UINT *ptr;
      ptr = (UINT*)obufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *ptr = (UINT) (0.5 + *oarr++);
        if(icr == 2) *(ptr+1) = (UINT) (0.5 + *oarr++);
        /*printf("%20.15e, %u \n", *(oarr-1), *ptr);*/
      }
    }
    break;
  case sInt: 
    {
      SINT *ptr;
      ptr = (SINT*)obufp;
      for(i = 0; i < nSam; i++, ptr += step) { 
        if (*oarr >= 0) *ptr = (SINT)  (0.5 + *oarr++);
        else            *ptr = (SINT) (-0.5 + *oarr++);
        if (icr == 2) {
          if (*oarr >= 0) *(ptr+1) = (SINT)  (0.5 + *oarr++);
          else            *(ptr+1) = (SINT) (-0.5 + *oarr++);
        }
        /*printf("%20.15e, %d \n", *(oarr-1), *ptr);*/
      }
    }
    break;
  case half_sreal:
    {
      SWORD *ptr;
      ptr = (SWORD*)obufp;
      float64_to_float16(oarr, ptr, nSam, step, icr);
    } /* */ 
    break;
  case sreal:
    {
      SREAL *ptr;
      ptr = (SREAL*)obufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *ptr = (SREAL) *oarr++;
        if(icr == 2) *(ptr+1) = (SREAL) *oarr++;
      }
    }
    break;
  case dreal: 
    {
      DREAL *ptr;
      ptr = (DREAL*)obufp;
      for(i = 0; i < nSam; i++, ptr += step) {
        *ptr = (DREAL) *oarr++;
        if(icr == 2) *(ptr+1) = (DREAL) *oarr++;
      }
    }
    break;
  default:  
    EXIT(E_TYPE, "otype");
  } 

  if ((format != machine_type) && (size > 1)) 
    swap(obufp, nSam, icr, step, size);

  nice_exit: 
    return(rval);           
}    
/******************************************************************************/
void nice_saturation(DREAL *arr, enum data_type dtype, SINT length, 
                     SINT *uSat, SINT *lSat)
{
  /*CHCK_INIT("nice_saturation")*/
  
  SINT  count;
  DREAL uLimit, lLimit;
  
  switch(dtype) {
    case ubyte:      lLimit = 0; uLimit = 255.0; break;
    case sbyte:      lLimit = -128.0; uLimit = 127.0; break;
    case uword:      lLimit = 0; uLimit = 65535.0; break;
    case sword:      lLimit = -32768.0; uLimit = 32767.0; break;
    case uInt:       lLimit = 0; uLimit = 4294967295; break;
    case sInt:       lLimit = -2147483648.0; uLimit = 2147483647.0; break;
    case half_sreal: lLimit = -FLT_MAX; uLimit = FLT_MAX; break;
    case sreal:      lLimit = -FLT_MAX; uLimit = FLT_MAX; break;
    case dreal:      length = 0; break;
  } 

  for (count = 0; count < length; count++) {
    if (arr[count] < lLimit) {
      arr[count] = lLimit;
      *lSat = *lSat + 1;
    } else if (arr[count] > uLimit) {
      if (arr[count]> uLimit) arr[count] = uLimit;
      *uSat = *uSat + 1;
    }
  }
}
/******************************************************************************/
DREAL round(DREAL inp)
{
  DREAL tmp;
  if (inp < 0) tmp = (SINT)(inp - 0.5);
  else tmp = (UINT)(inp + 0.5);
  return(tmp);
}
/******************************************************************************/
static void ErrTxt(ERR_T error, CHAR *str)
{
  switch (error) {
    case E_FOPE : strcpy(str, "File open error");                break;
    case E_FCLO : strcpy(str, "File closing error");             break;
    case E_FPOS : strcpy(str, "File positioning error");         break;
    case E_FREA : strcpy(str, "File reading error");             break;
    case E_FWRI : strcpy(str, "File writting error");            break;
    case E_ALLO : strcpy(str, "Allocation error");               break;
    case E_DIVZ : strcpy(str, "Division by zero");               break;
    case E_TRAC : strcpy(str, "Trace back from previous error"); break;
    case E_HEAD : strcpy(str, "Wrong header");                   break;
    case E_LARG : strcpy(str, "Arg. to large");                  break;
    case E_SMAL : strcpy(str, "Arg. to small");                  break;
    case E_NEXI : strcpy(str, "Non-existent item");              break;
    case E_NIMP : strcpy(str, "Non-implented function");         break;
    case E_TYPE : strcpy(str, "Unknown datatype");               break;
    case E_SCAN : strcpy(str, "Stdin error: end of file");       break;
    case E_SPEC : strcpy(str, "Special error");                  break;
    default     : strcpy(str, "Unknown error discriptor");
  }
}
/******************************************************************************/
void ErrMes(const CHAR *pid, ERR_T error, CHAR *name)
{
  CHAR str[80];
  
  ErrTxt(error, str);
  printf("\n!!! ERROR !!!\n");   
  printf("%s\n  In function   : %s\n  File/variable : %s\n", str, pid, name);
}
/******************************************************************************/
void WarMes(const CHAR *pid, ERR_T error, CHAR *name)
{
  CHAR str[80];
  
  ErrTxt(error, str);
  printf("\n!!! WARNING !!!\n");   
  printf("%s\n  In function   : %s\n  File/variable : %s\n", str, pid, name);
}
/******************************************************************************/