#include <stdlib.h>
#include <math.h>
#include "grib2.h"


void mkieee(g2float *a,g2int *rieee,g2int num)
//$$$  SUBPROGRAM DOCUMENTATION BLOCK
//                .      .    .                                       .
// SUBPROGRAM:    mkieee
//   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2002-10-29
//
// ABSTRACT: This subroutine stores a list of real values in
//   32-bit IEEE floating point format.
//
// PROGRAM HISTORY LOG:
// 2002-10-29  Gilbert
//
// USAGE:    mkieee(g2float *a,g2int *rieee,g2int num);
//   INPUT ARGUMENT LIST:
//     a        - Input array of floating point values.
//     num      - Number of floating point values to convert.
//
//   OUTPUT ARGUMENT LIST:
//     rieee    - Output array of data values in 32-bit IEEE format
//                stored in g2int integer array.  rieee must be allocated
//                with at least 4*num bytes of memory before calling this
//                function.
//
// REMARKS: None
//
// ATTRIBUTES:
//   LANGUAGE: C
//   MACHINE:
//
//$$$
{

      g2int  j,n,ieee,iexp,imant;
      double  /* alog2, */ atemp;

      static const double two23 = 8388608.0;              // pow(2,23)
      static const double two126 = 8.507059173023462e+37; // pow(2,126)

    //g2intu msk1=0x80000000;        // 10000000000000000000000000000000 binary
    //g2int msk2=0x7F800000;         // 01111111100000000000000000000000 binary
    //g2int msk3=0x007FFFFF;         // 00000000011111111111111111111111 binary

      // alog2=0.69314718;       //  ln(2.0)

      for (j=0;j<num;j++) {

        ieee=0;

        if (a[j] == 0.0) {
          rieee[j]=ieee;
          continue;
        }

//
//  Set Sign bit (bit 31 - leftmost bit)
//
        if (a[j] < 0.0) {
          ieee= (g2int) (1U << 31);
          atemp=-1.0*a[j];
        }
        else {
          ieee= 0 << 31;
          atemp=a[j];
        }
        //printf("sign %ld %x \n",ieee,ieee);
//
//  Determine exponent n with base 2
//
        if ( atemp >= 1.0 ) {
           n = 0;
           while ( int_power(2.0,n+1) <= atemp ) {
              n++;
           }
        }
        else {
           n = -1;
           while ( int_power(2.0,n) > atemp ) {
              n--;
           }
        }
        //n=(g2int)floor(log(atemp)/alog2);
        iexp=n+127;
        if (n >  127) iexp=255;     // overflow
        if (n < -127) iexp=0;
        //printf("exp %ld %ld \n",iexp,n);
        //      set exponent bits ( bits 30-23 )
        ieee = ieee | ( iexp << 23 );
//
//  Determine Mantissa
//
        if (iexp != 255) {
          if (iexp != 0)
            atemp=(atemp/int_power(2.0,n))-1.0;
          else
            atemp=atemp*two126;
          imant=(g2int)RINT(atemp*two23);
        }
        else {
          imant=0;
        }
        //printf("mant %ld %x \n",imant,imant);
        //      set mantissa bits ( bits 22-0 )
        ieee = ieee | imant;
//
//  Transfer IEEE bit string to rieee array
//
        rieee[j]=ieee;

      }

      return;

}