File: simpack.c

package info (click to toggle)
g2clib 1.4.0-2
  • links: PTS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 608 kB
  • ctags: 226
  • sloc: ansic: 4,995; makefile: 194
file content (181 lines) | stat: -rwxr-xr-x 6,457 bytes parent folder | download | duplicates (33)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
#include <stdlib.h>
#include <math.h>
#include "grib2.h"


void simpack(g2float *fld,g2int ndpts,g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
//$$$  SUBPROGRAM DOCUMENTATION BLOCK
//                .      .    .                                       .
// SUBPROGRAM:    simpack
//   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2002-11-06
//
// ABSTRACT: This subroutine packs up a data field using the simple
//   packing algorithm as defined in the GRIB2 documention.  It
//   also fills in GRIB2 Data Representation Template 5.0 with the
//   appropriate values.
//
// PROGRAM HISTORY LOG:
// 2002-11-06  Gilbert
//
// USAGE:    CALL simpack(fld,ndpts,idrstmpl,cpack,lcpack)
//   INPUT ARGUMENT LIST:
//     fld[]    - Contains the data values to pack
//     ndpts    - The number of data values in array fld[]
//     idrstmpl - Contains the array of values for Data Representation
//                Template 5.0
//                [0] = Reference value - ignored on input
//                [1] = Binary Scale Factor
//                [2] = Decimal Scale Factor
//                [3] = Number of bits used to pack data, if value is
//                      > 0 and  <= 31.
//                      If this input value is 0 or outside above range
//                      then the num of bits is calculated based on given 
//                      data and scale factors.
//                [4] = Original field type - currently ignored on input
//                      Data values assumed to be reals.
//
//   OUTPUT ARGUMENT LIST: 
//     idrstmpl - Contains the array of values for Data Representation
//                Template 5.0
//                [0] = Reference value - set by simpack routine.
//                [1] = Binary Scale Factor - unchanged from input
//                [2] = Decimal Scale Factor - unchanged from input
//                [3] = Number of bits used to pack data, unchanged from 
//                      input if value is between 0 and 31.
//                      If this input value is 0 or outside above range
//                      then the num of bits is calculated based on given 
//                      data and scale factors.
//                [4] = Original field type - currently set = 0 on output.
//                      Data values assumed to be reals.
//     cpack    - The packed data field
//     lcpack   - length of packed field starting at cpack.
//
// REMARKS: None
//
// ATTRIBUTES:
//   LANGUAGE: C
//   MACHINE:  
//
//$$$
{

      static g2int zero=0;
      g2int  *ifld;
      g2int  j,nbits,imin,imax,maxdif,nbittot,left;
      g2float  bscale,dscale,rmax,rmin,temp;
      double maxnum;
      static g2float alog2=0.69314718;       //  ln(2.0)
      
      bscale=int_power(2.0,-idrstmpl[1]);
      dscale=int_power(10.0,idrstmpl[2]);
      if (idrstmpl[3] <= 0 || idrstmpl[3] > 31)
         nbits=0;
      else
         nbits=idrstmpl[3];
//
//  Find max and min values in the data
//
      rmax=fld[0];
      rmin=fld[0];
      for (j=1;j<ndpts;j++) {
        if (fld[j] > rmax) rmax=fld[j];
        if (fld[j] < rmin) rmin=fld[j];
      }
     
      ifld=calloc(ndpts,sizeof(g2int));
//
//  If max and min values are not equal, pack up field.
//  If they are equal, we have a constant field, and the reference
//  value (rmin) is the value for each point in the field and
//  set nbits to 0.
//
      if (rmin != rmax) {
        //
        //  Determine which algorithm to use based on user-supplied 
        //  binary scale factor and number of bits.
        //
        if (nbits==0 && idrstmpl[1]==0) {
           //
           //  No binary scaling and calculate minumum number of 
           //  bits in which the data will fit.
           //
           imin=(g2int)rint(rmin*dscale);
           imax=(g2int)rint(rmax*dscale);
           maxdif=imax-imin;
           temp=log((double)(maxdif+1))/alog2;
           nbits=(g2int)ceil(temp);
           rmin=(g2float)imin;
           //   scale data
           for(j=0;j<ndpts;j++)
             ifld[j]=(g2int)rint(fld[j]*dscale)-imin;
        }
        else if (nbits!=0 && idrstmpl[1]==0) {
           //
           //  Use minimum number of bits specified by user and
           //  adjust binary scaling factor to accomodate data.
           //
           rmin=rmin*dscale;
           rmax=rmax*dscale;
           maxnum=int_power(2.0,nbits)-1;
           temp=log(maxnum/(rmax-rmin))/alog2;
           idrstmpl[1]=(g2int)ceil(-1.0*temp);
           bscale=int_power(2.0,-idrstmpl[1]);
           //   scale data
           for (j=0;j<ndpts;j++)
             ifld[j]=(g2int)rint(((fld[j]*dscale)-rmin)*bscale);
        }
        else if (nbits==0 && idrstmpl[1]!=0) {
           //
           //  Use binary scaling factor and calculate minumum number of 
           //  bits in which the data will fit.
           //
           rmin=rmin*dscale;
           rmax=rmax*dscale;
           maxdif=(g2int)rint((rmax-rmin)*bscale);
           temp=log((double)(maxdif+1))/alog2;
           nbits=(g2int)ceil(temp);
           //   scale data
           for (j=0;j<ndpts;j++)
             ifld[j]=(g2int)rint(((fld[j]*dscale)-rmin)*bscale);
        }
        else if (nbits!=0 && idrstmpl[1]!=0) {
           //
           //  Use binary scaling factor and use minumum number of 
           //  bits specified by user.   Dangerous - may loose
           //  information if binary scale factor and nbits not set
           //  properly by user.
           //
           rmin=rmin*dscale;
           //   scale data
           for (j=0;j<ndpts;j++)
             ifld[j]=(g2int)rint(((fld[j]*dscale)-rmin)*bscale);
        }
        //
        //  Pack data, Pad last octet with Zeros, if necessary,
        //  and calculate the length of the packed data in bytes
        //
        sbits(cpack,ifld+0,0,nbits,0,ndpts);
        nbittot=nbits*ndpts;
        left=8-(nbittot%8);
        if (left != 8) {
          sbit(cpack,&zero,nbittot,left);   // Pad with zeros to fill Octet
          nbittot=nbittot+left;
        }
        *lcpack=nbittot/8;
      }
      else {
        nbits=0;
        *lcpack=0;
      }

//
//  Fill in ref value and number of bits in Template 5.0
//
      //printf("SAGmkieee %f\n",rmin);
      mkieee(&rmin,idrstmpl+0,1);   // ensure reference value is IEEE format
      //printf("SAGmkieee %ld\n",idrstmpl[0]);
      idrstmpl[3]=nbits;
      idrstmpl[4]=0;         // original data were reals

      free(ifld);
}