#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;

}