124 lines
3 KiB
C
124 lines
3 KiB
C
|
#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 double two23,two126;
|
||
|
|
static g2int test=0;
|
||
|
|
//g2intu msk1=0x80000000; // 10000000000000000000000000000000 binary
|
||
|
|
//g2int msk2=0x7F800000; // 01111111100000000000000000000000 binary
|
||
|
|
//g2int msk3=0x007FFFFF; // 00000000011111111111111111111111 binary
|
||
|
|
|
||
|
|
if ( test == 0 ) {
|
||
|
|
two23=(double)int_power(2.0,23);
|
||
|
|
two126=(double)int_power(2.0,126);
|
||
|
|
test=1;
|
||
|
|
}
|
||
|
|
|
||
|
|
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= 1 << 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;
|
||
|
|
|
||
|
|
}
|