Line data    Source code

       1             : #include <stdlib.h>
       2             : #include <math.h>
       3             : #include "grib2.h"
       4             : 
       5             : 
       6          77 : 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 const double two23 = 8388608.0;              // pow(2,23)
      42             :       static const double two126 = 8.507059173023462e+37; // pow(2,126)
      43             : 
      44             :     //g2intu msk1=0x80000000;        // 10000000000000000000000000000000 binary
      45             :     //g2int msk2=0x7F800000;         // 01111111100000000000000000000000 binary
      46             :     //g2int msk3=0x007FFFFF;         // 00000000011111111111111111111111 binary
      47             : 
      48             :       // alog2=0.69314718;       //  ln(2.0)
      49             : 
      50         154 :       for (j=0;j<num;j++) {
      51             : 
      52          77 :         ieee=0;
      53             : 
      54          77 :         if (a[j] == 0.0) {
      55          12 :           rieee[j]=ieee;
      56          12 :           continue;
      57             :         }
      58             : 
      59             : //
      60             : //  Set Sign bit (bit 31 - leftmost bit)
      61             : //
      62          65 :         if (a[j] < 0.0) {
      63           1 :           ieee= (g2int) (1U << 31);
      64           1 :           atemp=-1.0*a[j];
      65             :         }
      66             :         else {
      67          64 :           ieee= 0 << 31;
      68          64 :           atemp=a[j];
      69             :         }
      70             :         //printf("sign %ld %x \n",ieee,ieee);
      71             : //
      72             : //  Determine exponent n with base 2
      73             : //
      74          65 :         if ( atemp >= 1.0 ) {
      75          65 :            n = 0;
      76         491 :            while ( int_power(2.0,n+1) <= atemp ) {
      77         426 :               n++;
      78             :            }
      79             :         }
      80             :         else {
      81           0 :            n = -1;
      82           0 :            while ( int_power(2.0,n) > atemp ) {
      83           0 :               n--;
      84             :            }
      85             :         }
      86             :         //n=(g2int)floor(log(atemp)/alog2);
      87          65 :         iexp=n+127;
      88          65 :         if (n >  127) iexp=255;     // overflow
      89          65 :         if (n < -127) iexp=0;
      90             :         //printf("exp %ld %ld \n",iexp,n);
      91             :         //      set exponent bits ( bits 30-23 )
      92          65 :         ieee = ieee | ( iexp << 23 );
      93             : //
      94             : //  Determine Mantissa
      95             : //
      96          65 :         if (iexp != 255) {
      97          65 :           if (iexp != 0)
      98          65 :             atemp=(atemp/int_power(2.0,n))-1.0;
      99             :           else
     100           0 :             atemp=atemp*two126;
     101          65 :           imant=(g2int)RINT(atemp*two23);
     102             :         }
     103             :         else {
     104           0 :           imant=0;
     105             :         }
     106             :         //printf("mant %ld %x \n",imant,imant);
     107             :         //      set mantissa bits ( bits 22-0 )
     108          65 :         ieee = ieee | imant;
     109             : //
     110             : //  Transfer IEEE bit string to rieee array
     111             : //
     112          65 :         rieee[j]=ieee;
     113             : 
     114             :       }
     115             : 
     116          77 :       return;
     117             : 
     118             : }