1                 : #include <stdlib.h>
       2                 : #include <math.h>
       3                 : #include "grib2.h"
       4                 : 
       5                 : 
       6               0 : void mkieee(g2float *a,g2int *rieee,g2int num)
       7                 : //$$$  SUBPROGRAM DOCUMENTATION BLOCK
       8                 : //                .      .    .                                       .
       9                 : // SUBPROGRAM:    mkieee 
      10                 : //   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2002-10-29
      11                 : //
      12                 : // ABSTRACT: This subroutine stores a list of real values in 
      13                 : //   32-bit IEEE floating point format.
      14                 : //
      15                 : // PROGRAM HISTORY LOG:
      16                 : // 2002-10-29  Gilbert
      17                 : //
      18                 : // USAGE:    mkieee(g2float *a,g2int *rieee,g2int num);
      19                 : //   INPUT ARGUMENT LIST:
      20                 : //     a        - Input array of floating point values.
      21                 : //     num      - Number of floating point values to convert.
      22                 : //
      23                 : //   OUTPUT ARGUMENT LIST:      
      24                 : //     rieee    - Output array of data values in 32-bit IEEE format
      25                 : //                stored in g2int integer array.  rieee must be allocated
      26                 : //                with at least 4*num bytes of memory before calling this
      27                 : //                function.
      28                 : //
      29                 : // REMARKS: None
      30                 : //
      31                 : // ATTRIBUTES:
      32                 : //   LANGUAGE: C
      33                 : //   MACHINE:  
      34                 : //
      35                 : //$$$
      36                 : {
      37                 : 
      38                 :       g2int  j,n,ieee,iexp,imant;
      39                 :       double  alog2,atemp;
      40                 : 
      41                 :       static double  two23,two126;
      42                 :       static g2int test=0;
      43                 :     //g2intu msk1=0x80000000;        // 10000000000000000000000000000000 binary
      44                 :     //g2int msk2=0x7F800000;         // 01111111100000000000000000000000 binary
      45                 :     //g2int msk3=0x007FFFFF;         // 00000000011111111111111111111111 binary
      46                 : 
      47               0 :       if ( test == 0 ) {
      48               0 :          two23=(double)int_power(2.0,23);
      49               0 :          two126=(double)int_power(2.0,126);
      50               0 :          test=1;
      51                 :       }
      52                 : 
      53               0 :       alog2=0.69314718;       //  ln(2.0)
      54                 : 
      55               0 :       for (j=0;j<num;j++) {
      56                 :       
      57               0 :         ieee=0;
      58                 : 
      59               0 :         if (a[j] == 0.0) {
      60               0 :           rieee[j]=ieee;
      61               0 :           continue;
      62                 :         }
      63                 :         
      64                 : //
      65                 : //  Set Sign bit (bit 31 - leftmost bit)
      66                 : //
      67               0 :         if (a[j] < 0.0) {
      68               0 :           ieee= 1 << 31;
      69               0 :           atemp=-1.0*a[j];
      70                 :         }
      71                 :         else {
      72               0 :           ieee= 0 << 31;
      73               0 :           atemp=a[j];
      74                 :         }
      75                 :         //printf("sign %ld %x \n",ieee,ieee);
      76                 : //
      77                 : //  Determine exponent n with base 2
      78                 : //
      79               0 :         if ( atemp >= 1.0 ) {
      80               0 :            n = 0;
      81               0 :            while ( int_power(2.0,n+1) <= atemp ) {
      82               0 :               n++;
      83                 :            }
      84                 :         }
      85                 :         else {
      86               0 :            n = -1;
      87               0 :            while ( int_power(2.0,n) > atemp ) {
      88               0 :               n--;
      89                 :            }
      90                 :         }
      91                 :         //n=(g2int)floor(log(atemp)/alog2);
      92               0 :         iexp=n+127;
      93               0 :         if (n >  127) iexp=255;     // overflow
      94               0 :         if (n < -127) iexp=0;
      95                 :         //printf("exp %ld %ld \n",iexp,n);
      96                 :         //      set exponent bits ( bits 30-23 )
      97               0 :         ieee = ieee | ( iexp << 23 );
      98                 : //
      99                 : //  Determine Mantissa
     100                 : // 
     101               0 :         if (iexp != 255) {
     102               0 :           if (iexp != 0) 
     103               0 :             atemp=(atemp/int_power(2.0,n))-1.0;
     104                 :           else
     105               0 :             atemp=atemp*two126;
     106               0 :           imant=(g2int)RINT(atemp*two23);
     107                 :         }
     108                 :         else {
     109               0 :           imant=0;
     110                 :         }
     111                 :         //printf("mant %ld %x \n",imant,imant);
     112                 :         //      set mantissa bits ( bits 22-0 )
     113               0 :         ieee = ieee | imant;
     114                 : //
     115                 : //  Transfer IEEE bit string to rieee array
     116                 : //
     117               0 :         rieee[j]=ieee;
     118                 : 
     119                 :       }
     120                 : 
     121                 :       return;
     122                 : 
     123                 : }