Line data Source code
1 : #include <stdlib.h>
2 : #include <math.h>
3 : #include <limits.h>
4 : #include <float.h>
5 : #include "grib2.h"
6 :
7 4 : void misspack(g2float *fld,g2int ndpts,g2int idrsnum,g2int *idrstmpl,
8 : unsigned char *cpack, g2int *lcpack)
9 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
10 : // . . . .
11 : // SUBPROGRAM: misspack
12 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-06-21
13 : //
14 : // ABSTRACT: This subroutine packs up a data field using a complex
15 : // packing algorithm as defined in the GRIB2 documentation. It
16 : // supports GRIB2 complex packing templates with or without
17 : // spatial differences (i.e. DRTs 5.2 and 5.3).
18 : // It also fills in GRIB2 Data Representation Template 5.2 or 5.3
19 : // with the appropriate values.
20 : // This version assumes that Missing Value Management is being used and that
21 : // 1 or 2 missing values appear in the data.
22 : //
23 : // PROGRAM HISTORY LOG:
24 : // 2000-06-21 Gilbert
25 : //
26 : // USAGE: misspack(g2float *fld,g2int ndpts,g2int idrsnum,g2int *idrstmpl,
27 : // unsigned char *cpack, g2int *lcpack)
28 : // INPUT ARGUMENT LIST:
29 : // fld[] - Contains the data values to pack
30 : // ndpts - The number of data values in array fld[]
31 : // idrsnum - Data Representation Template number 5.N
32 : // Must equal 2 or 3.
33 : // idrstmpl - Contains the array of values for Data Representation
34 : // Template 5.2 or 5.3
35 : // [0] = Reference value - ignored on input
36 : // [1] = Binary Scale Factor
37 : // [2] = Decimal Scale Factor
38 : // .
39 : // .
40 : // [6] = Missing value management
41 : // [7] = Primary missing value
42 : // [8] = Secondary missing value
43 : // .
44 : // .
45 : // [16] = Order of Spatial Differencing ( 1 or 2 )
46 : // .
47 : // .
48 : //
49 : // OUTPUT ARGUMENT LIST:
50 : // idrstmpl - Contains the array of values for Data Representation
51 : // Template 5.3
52 : // [0] = Reference value - set by misspack routine.
53 : // [1] = Binary Scale Factor - unchanged from input
54 : // [2] = Decimal Scale Factor - unchanged from input
55 : // .
56 : // .
57 : // cpack - The packed data field (character*1 array)
58 : // *lcpack - length of packed field cpack(). (or -1 in case of error)
59 : //
60 : // REMARKS: None
61 : //
62 : // ATTRIBUTES:
63 : // LANGUAGE: C
64 : // MACHINE:
65 : //
66 : //$$$
67 : {
68 :
69 : g2int *ifld, *ifldmiss, *jfld;
70 : g2int *jmin, *jmax, *lbit;
71 4 : const g2int zero=0;
72 : g2int *gref, *gwidth, *glen;
73 : g2int glength, grpwidth;
74 4 : g2int i, n, iofst, ival1, ival2, isd, minsd, nbitsd = 0;
75 : g2int nbitsgref, left, iwmax, ngwidthref, nbitsgwidth, ilmax;
76 : g2int nglenref, nglenlast, nbitsglen /* , ij */;
77 : g2int j, missopt, nonmiss, itemp, maxorig, nbitorig, miss1, miss2;
78 : g2int ngroups, ng, num0, num1, num2;
79 4 : g2int imax, lg, mtemp, ier = 0, igmax;
80 : g2int kfildo, minpk, inc, maxgrps, ibit, jbit, kbit, novref, lbitref;
81 : g2float rmissp, rmisss, bscale, dscale, rmin, rmax, temp;
82 4 : const g2int simple_alg = 0;
83 4 : const g2float alog2=0.69314718f; // ln(2.0)
84 4 : const g2int one=1;
85 :
86 4 : bscale=(float)int_power(2.0,-idrstmpl[1]);
87 4 : dscale=(float)int_power(10.0,idrstmpl[2]);
88 4 : missopt=idrstmpl[6];
89 4 : if ( missopt != 1 && missopt != 2 ) {
90 0 : printf("misspack: Unrecognized option.\n");
91 0 : *lcpack=-1;
92 0 : return;
93 : }
94 : else { // Get missing values
95 4 : rdieee(idrstmpl+7,&rmissp,1);
96 4 : if (missopt == 2) rdieee(idrstmpl+8,&rmisss,1);
97 : }
98 : //
99 : // Find min value of non-missing values in the data,
100 : // AND set up missing value mapping of the field.
101 : //
102 4 : ifldmiss = calloc(ndpts,sizeof(g2int));
103 4 : if( ifldmiss == NULL )
104 : {
105 0 : *lcpack = -1;
106 0 : return;
107 : }
108 4 : rmin=FLT_MAX;
109 4 : rmax=-FLT_MAX;
110 4 : if ( missopt == 1 ) { // Primary missing value only
111 1084310 : for ( j=0; j<ndpts; j++) {
112 1084310 : if (fld[j] == rmissp) {
113 672513 : ifldmiss[j]=1;
114 : }
115 : else {
116 411793 : ifldmiss[j]=0;
117 411793 : if (fld[j] < rmin) rmin=fld[j];
118 411793 : if (fld[j] > rmax) rmax=fld[j];
119 : }
120 : }
121 : }
122 4 : if ( missopt == 2 ) { // Primary and secondary missing values
123 0 : for ( j=0; j<ndpts; j++ ) {
124 0 : if (fld[j] == rmissp) {
125 0 : ifldmiss[j]=1;
126 : }
127 0 : else if (fld[j] == rmisss) {
128 0 : ifldmiss[j]=2;
129 : }
130 : else {
131 0 : ifldmiss[j]=0;
132 0 : if (fld[j] < rmin) rmin=fld[j];
133 0 : if (fld[j] > rmax) rmax=fld[j];
134 : }
135 : }
136 : }
137 :
138 4 : if( !(floor((double)rmin*dscale) >= -FLT_MAX && floor((double)rmin*dscale) <= FLT_MAX) )
139 : {
140 0 : fprintf(stderr,
141 : "Scaled min value not representable on IEEE754 "
142 : "single precision float\n");
143 0 : *lcpack = -1;
144 0 : free(ifldmiss);
145 0 : return;
146 : }
147 4 : if (idrstmpl[1] == 0)
148 : {
149 3 : double max_diff = RINT(rmax*dscale) - RINT(rmin*dscale);
150 3 : if( !(max_diff >= 0 && max_diff <= INT_MAX) )
151 : {
152 0 : fprintf(stderr,
153 : "Difference of scaled max value with scaled min value "
154 : "not representable on int \n");
155 0 : *lcpack = -1;
156 0 : free(ifldmiss);
157 0 : return;
158 : }
159 : }
160 : else
161 : {
162 1 : double max_diff = RINT(rmax*dscale - rmin*dscale) * bscale;
163 1 : if( !(max_diff >= 0 && max_diff <= INT_MAX) )
164 : {
165 0 : fprintf(stderr,
166 : "Difference of scaled max value with scaled min value "
167 : "not representable on int \n");
168 0 : *lcpack = -1;
169 0 : free(ifldmiss);
170 0 : return;
171 : }
172 : }
173 : //
174 : // Allocate work arrays:
175 : // Note: -ifldmiss[j],j=0,ndpts-1 is a map of original field indicating
176 : // which of the original data values
177 : // are primary missing (1), secondary missing (2) or non-missing (0).
178 : // -jfld[j],j=0,nonmiss-1 is a subarray of just the non-missing values
179 : // from the original field.
180 : //
181 : //if (rmin != rmax) {
182 4 : iofst=0;
183 4 : ifld = calloc(ndpts,sizeof(g2int));
184 4 : jfld = calloc(ndpts,sizeof(g2int));
185 4 : gref = calloc(ndpts,sizeof(g2int));
186 4 : gwidth = calloc(ndpts,sizeof(g2int));
187 4 : glen = calloc(ndpts,sizeof(g2int));
188 4 : if( ifld == NULL || jfld == NULL || gref == NULL || gwidth == NULL ||
189 : glen == NULL )
190 : {
191 0 : free(ifld);
192 0 : free(jfld);
193 0 : free(ifldmiss);
194 0 : free(gref);
195 0 : free(gwidth);
196 0 : free(glen);
197 0 : *lcpack = -1;
198 0 : return;
199 : }
200 : //
201 : // Scale original data
202 : //
203 4 : nonmiss=0;
204 4 : if (idrstmpl[1] == 0) { // No binary scaling
205 3 : rmin=(g2float)RINT(rmin*dscale);
206 1061480 : for ( j=0; j<ndpts; j++) {
207 1061470 : if (ifldmiss[j] == 0) {
208 392716 : jfld[nonmiss]=(g2int)(RINT(fld[j]*dscale)-rmin);
209 392716 : nonmiss++;
210 : }
211 : }
212 : }
213 : else { // Use binary scaling factor
214 1 : rmin=rmin*dscale;
215 : //rmax=rmax*dscale;
216 22834 : for ( j=0; j<ndpts; j++ ) {
217 22833 : if (ifldmiss[j] == 0) {
218 19077 : jfld[nonmiss]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
219 19077 : nonmiss++;
220 : }
221 : }
222 : }
223 : //
224 : // Calculate Spatial differences, if using DRS Template 5.3
225 : //
226 4 : if (idrsnum == 3) { // spatial differences
227 0 : if (idrstmpl[16]!=1 && idrstmpl[16]!=2) idrstmpl[16]=2;
228 0 : if (idrstmpl[16] == 1) { // first order
229 0 : ival1=jfld[0];
230 0 : for ( j=nonmiss-1; j>0; j--)
231 0 : jfld[j]=jfld[j]-jfld[j-1];
232 0 : jfld[0]=0;
233 : }
234 0 : else if (idrstmpl[16] == 2) { // second order
235 0 : ival1=jfld[0];
236 0 : ival2=jfld[1];
237 0 : for ( j=nonmiss-1; j>1; j--)
238 0 : jfld[j]=jfld[j]-(2*jfld[j-1])+jfld[j-2];
239 0 : jfld[0]=0;
240 0 : jfld[1]=0;
241 : }
242 : //
243 : // subtract min value from spatial diff field
244 : //
245 0 : isd=idrstmpl[16];
246 0 : minsd=jfld[isd];
247 0 : for ( j=isd; j<nonmiss; j++ ) if ( jfld[j] < minsd ) minsd=jfld[j];
248 0 : for ( j=isd; j<nonmiss; j++ ) jfld[j]=jfld[j]-minsd;
249 : //
250 : // find num of bits need to store minsd and add 1 extra bit
251 : // to indicate sign
252 : //
253 0 : temp=(float)(log((double)(abs(minsd)+1))/alog2);
254 0 : nbitsd=(g2int)ceil(temp)+1;
255 : //
256 : // find num of bits need to store ifld[0] ( and ifld[1]
257 : // if using 2nd order differencing )
258 : //
259 0 : maxorig=ival1;
260 0 : if (idrstmpl[16]==2 && ival2>ival1) maxorig=ival2;
261 0 : temp=(float)(log((double)(maxorig+1))/alog2);
262 0 : nbitorig=(g2int)ceil(temp)+1;
263 0 : if (nbitorig > nbitsd) nbitsd=nbitorig;
264 : // increase number of bits to even multiple of 8 ( octet )
265 0 : if ( (nbitsd%8) != 0) nbitsd=nbitsd+(8-(nbitsd%8));
266 : //
267 : // Store extra spatial differencing info into the packed
268 : // data section.
269 : //
270 0 : if (nbitsd != 0) {
271 : // pack first original value
272 0 : if (ival1 >= 0) {
273 0 : sbit(cpack,&ival1,iofst,nbitsd);
274 0 : iofst=iofst+nbitsd;
275 : }
276 : else {
277 0 : sbit(cpack,&one,iofst,1);
278 0 : iofst=iofst+1;
279 0 : itemp=abs(ival1);
280 0 : sbit(cpack,&itemp,iofst,nbitsd-1);
281 0 : iofst=iofst+nbitsd-1;
282 : }
283 0 : if (idrstmpl[16] == 2) {
284 : // pack second original value
285 0 : if (ival2 >= 0) {
286 0 : sbit(cpack,&ival2,iofst,nbitsd);
287 0 : iofst=iofst+nbitsd;
288 : }
289 : else {
290 0 : sbit(cpack,&one,iofst,1);
291 0 : iofst=iofst+1;
292 0 : itemp=abs(ival2);
293 0 : sbit(cpack,&itemp,iofst,nbitsd-1);
294 0 : iofst=iofst+nbitsd-1;
295 : }
296 : }
297 : // pack overall min of spatial differences
298 0 : if (minsd >= 0) {
299 0 : sbit(cpack,&minsd,iofst,nbitsd);
300 0 : iofst=iofst+nbitsd;
301 : }
302 : else {
303 0 : sbit(cpack,&one,iofst,1);
304 0 : iofst=iofst+1;
305 0 : itemp=abs(minsd);
306 0 : sbit(cpack,&itemp,iofst,nbitsd-1);
307 0 : iofst=iofst+nbitsd-1;
308 : }
309 : }
310 : //print *,'SDp ',ival1,ival2,minsd,nbitsd
311 : } // end of spatial diff section
312 : //
313 : // Expand non-missing data values to original grid.
314 : //
315 4 : miss1=jfld[0];
316 411797 : for ( j=0; j<nonmiss; j++) if (jfld[j] < miss1) miss1 = jfld[j];
317 4 : if( miss1 <= INT_MIN + 1 )
318 : {
319 : // E. Rouault: no idea if this is correct, but avoids integer
320 : // wrap over
321 0 : miss1++;
322 0 : miss2 = miss1 + 1;
323 : }
324 : else
325 : {
326 4 : miss1--;
327 4 : miss2=miss1-1;
328 : }
329 4 : n=0;
330 1084310 : for ( j=0; j<ndpts; j++) {
331 1084310 : if ( ifldmiss[j] == 0 ) {
332 411793 : ifld[j]=jfld[n];
333 411793 : n++;
334 : }
335 672513 : else if ( ifldmiss[j] == 1 ) {
336 672513 : ifld[j]=miss1;
337 : }
338 0 : else if ( ifldmiss[j] == 2 ) {
339 0 : ifld[j]=miss2;
340 : }
341 : }
342 : //
343 : // Determine Groups to be used.
344 : //
345 4 : if ( simple_alg == 1 ) {
346 : // set group length to 10 : calculate number of groups
347 : // and length of last group
348 0 : ngroups=ndpts/10;
349 0 : for (j=0;j<ngroups;j++) glen[j]=10;
350 0 : itemp=ndpts%10;
351 0 : if (itemp != 0) {
352 0 : ngroups++;
353 0 : glen[ngroups-1]=itemp;
354 : }
355 : }
356 : else {
357 : // Use Dr. Glahn's algorithm for determining grouping.
358 : //
359 4 : kfildo=6;
360 4 : minpk=10;
361 4 : inc=1;
362 4 : maxgrps=(ndpts/minpk)+1;
363 4 : jmin = calloc(maxgrps,sizeof(g2int));
364 4 : jmax = calloc(maxgrps,sizeof(g2int));
365 4 : lbit = calloc(maxgrps,sizeof(g2int));
366 4 : pack_gp(&kfildo,ifld,&ndpts,&missopt,&minpk,&inc,&miss1,&miss2,
367 : jmin,jmax,lbit,glen,&maxgrps,&ngroups,&ibit,&jbit,
368 : &kbit,&novref,&lbitref,&ier);
369 : //printf("SAGier = %d %d %d %d %d %d\n",ier,ibit,jbit,kbit,novref,lbitref);
370 41281 : for ( ng=0; ng<ngroups; ng++) glen[ng]=glen[ng]+novref;
371 4 : free(jmin);
372 4 : free(jmax);
373 4 : free(lbit);
374 4 : if( ier != 0 )
375 : {
376 0 : free(ifld);
377 0 : free(jfld);
378 0 : free(ifldmiss);
379 0 : free(gref);
380 0 : free(gwidth);
381 0 : free(glen);
382 0 : *lcpack = -1;
383 0 : return;
384 : }
385 : }
386 : //
387 : // For each group, find the group's reference value (min)
388 : // and the number of bits needed to hold the remaining values
389 : //
390 4 : n=0;
391 41281 : for ( ng=0; ng<ngroups; ng++) {
392 : // how many of each type?
393 41277 : num0=num1=num2=0;
394 1125580 : for (j=n; j<n+glen[ng]; j++) {
395 1084310 : if (ifldmiss[j] == 0 ) num0++;
396 1084310 : if (ifldmiss[j] == 1 ) num1++;
397 1084310 : if (ifldmiss[j] == 2 ) num2++;
398 : }
399 41277 : if ( num0 == 0 ) { // all missing values
400 5341 : if ( num1 == 0 ) { // all secondary missing
401 0 : gref[ng]=-2;
402 0 : gwidth[ng]=0;
403 : }
404 5341 : else if ( num2 == 0 ) { // all primary missing
405 5341 : gref[ng]=-1;
406 5341 : gwidth[ng]=0;
407 : }
408 : else { // both primary and secondary
409 0 : gref[ng]=0;
410 0 : gwidth[ng]=1;
411 : }
412 : }
413 : else { // contains some non-missing data
414 : // find max and min values of group
415 35936 : gref[ng]=2147483647;
416 35936 : imax=-2147483647;
417 35936 : j=n;
418 490683 : for ( lg=0; lg<glen[ng]; lg++ ) {
419 454747 : if ( ifldmiss[j] == 0 ) {
420 411793 : if (ifld[j] < gref[ng]) gref[ng]=ifld[j];
421 411793 : if (ifld[j] > imax) imax=ifld[j];
422 : }
423 454747 : j++;
424 : }
425 35936 : if (missopt == 1) imax=imax+1;
426 35936 : if (missopt == 2) imax=imax+2;
427 : // calc num of bits needed to hold data
428 35936 : if ( gref[ng] != imax ) {
429 35936 : temp=(float)(log((double)(imax-gref[ng]+1))/alog2);
430 35936 : gwidth[ng]=(g2int)ceil(temp);
431 : }
432 : else {
433 0 : gwidth[ng]=0;
434 : }
435 : }
436 : // Subtract min from data
437 41277 : j=n;
438 41277 : mtemp=(g2int)int_power(2.,gwidth[ng]);
439 1125580 : for ( lg=0; lg<glen[ng]; lg++ ) {
440 1084310 : if (ifldmiss[j] == 0) // non-missing
441 411793 : ifld[j]=ifld[j]-gref[ng];
442 672513 : else if (ifldmiss[j] == 1) // primary missing
443 672513 : ifld[j]=mtemp-1;
444 0 : else if (ifldmiss[j] == 2) // secondary missing
445 0 : ifld[j]=mtemp-2;
446 :
447 1084310 : j++;
448 : }
449 : // increment fld array counter
450 41277 : n=n+glen[ng];
451 : }
452 : //
453 : // Find max of the group references and calc num of bits needed
454 : // to pack each groups reference value, then
455 : // pack up group reference values
456 : //
457 : //printf(" GREFS: ");
458 : //for (j=0;j<ngroups;j++) printf(" %d",gref[j]); printf("\n");
459 4 : igmax=gref[0];
460 41277 : for (j=1;j<ngroups;j++) if (gref[j] > igmax) igmax=gref[j];
461 4 : if (missopt == 1) igmax=igmax+1;
462 4 : if (missopt == 2) igmax=igmax+2;
463 4 : if (igmax != 0) {
464 4 : temp=(float)(log((double)(igmax+1))/alog2);
465 4 : nbitsgref=(g2int)ceil(temp);
466 4 : if( nbitsgref < 0 || nbitsgref >= 31 )
467 : {
468 0 : free(ifld);
469 0 : free(jfld);
470 0 : free(ifldmiss);
471 0 : free(gref);
472 0 : free(gwidth);
473 0 : free(glen);
474 0 : *lcpack = -1;
475 0 : return;
476 : }
477 : // reset the ref values of any "missing only" groups.
478 4 : mtemp=(g2int)int_power(2.,nbitsgref);
479 41281 : for ( j=0; j<ngroups; j++ ) {
480 41277 : if (gref[j] == -1) gref[j]=mtemp-1;
481 41277 : if (gref[j] == -2) gref[j]=mtemp-2;
482 : }
483 4 : sbits(cpack,gref,iofst,nbitsgref,0,ngroups);
484 4 : itemp=nbitsgref*ngroups;
485 4 : iofst=iofst+itemp;
486 : // Pad last octet with Zeros, if necessary,
487 4 : if ( (itemp%8) != 0) {
488 4 : left=8-(itemp%8);
489 4 : sbit(cpack,&zero,iofst,left);
490 4 : iofst=iofst+left;
491 : }
492 : }
493 : else {
494 0 : nbitsgref=0;
495 : }
496 : //
497 : // Find max/min of the group widths and calc num of bits needed
498 : // to pack each groups width value, then
499 : // pack up group width values
500 : //
501 : //write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups)
502 4 : iwmax=gwidth[0];
503 4 : ngwidthref=gwidth[0];
504 41277 : for (j=1;j<ngroups;j++) {
505 41273 : if (gwidth[j] > iwmax) iwmax=gwidth[j];
506 41273 : if (gwidth[j] < ngwidthref) ngwidthref=gwidth[j];
507 : }
508 4 : if (iwmax != ngwidthref) {
509 4 : temp=(float)(log((double)(iwmax-ngwidthref+1))/alog2);
510 4 : nbitsgwidth=(g2int)ceil(temp);
511 41281 : for ( i=0; i<ngroups; i++) gwidth[i]=gwidth[i]-ngwidthref;
512 4 : sbits(cpack,gwidth,iofst,nbitsgwidth,0,ngroups);
513 4 : itemp=nbitsgwidth*ngroups;
514 4 : iofst=iofst+itemp;
515 : // Pad last octet with Zeros, if necessary,
516 4 : if ( (itemp%8) != 0) {
517 3 : left=8-(itemp%8);
518 3 : sbit(cpack,&zero,iofst,left);
519 3 : iofst=iofst+left;
520 : }
521 : }
522 : else {
523 0 : nbitsgwidth=0;
524 0 : for (i=0;i<ngroups;i++) gwidth[i]=0;
525 : }
526 : //
527 : // Find max/min of the group lengths and calc num of bits needed
528 : // to pack each groups length value, then
529 : // pack up group length values
530 : //
531 : //printf(" GLENS: ");
532 : //for (j=0;j<ngroups;j++) printf(" %d",glen[j]); printf("\n");
533 4 : ilmax=glen[0];
534 4 : nglenref=glen[0];
535 41273 : for (j=1;j<ngroups-1;j++) {
536 41269 : if (glen[j] > ilmax) ilmax=glen[j];
537 41269 : if (glen[j] < nglenref) nglenref=glen[j];
538 : }
539 4 : nglenlast=glen[ngroups-1];
540 4 : if (ilmax != nglenref) {
541 4 : temp=(float)(log((double)(ilmax-nglenref+1))/alog2);
542 4 : nbitsglen=(g2int)ceil(temp);
543 41277 : for ( i=0; i<ngroups-1; i++) glen[i]=glen[i]-nglenref;
544 4 : sbits(cpack,glen,iofst,nbitsglen,0,ngroups);
545 4 : itemp=nbitsglen*ngroups;
546 4 : iofst=iofst+itemp;
547 : // Pad last octet with Zeros, if necessary,
548 4 : if ( (itemp%8) != 0) {
549 1 : left=8-(itemp%8);
550 1 : sbit(cpack,&zero,iofst,left);
551 1 : iofst=iofst+left;
552 : }
553 : }
554 : else {
555 0 : nbitsglen=0;
556 0 : for (i=0;i<ngroups;i++) glen[i]=0;
557 : }
558 : //
559 : // For each group, pack data values
560 : //
561 : //write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
562 4 : n=0;
563 : // ij=0;
564 41281 : for ( ng=0; ng<ngroups; ng++) {
565 41277 : glength=glen[ng]+nglenref;
566 41277 : if (ng == (ngroups-1) ) glength=nglenlast;
567 41277 : grpwidth=gwidth[ng]+ngwidthref;
568 : //write(77,*)'NGP ',ng,grpwidth,glength,gref(ng)
569 41277 : if ( grpwidth != 0 ) {
570 35936 : sbits(cpack,ifld+n,iofst,grpwidth,0,glength);
571 35936 : iofst=iofst+(grpwidth*glength);
572 : }
573 : // do kk=1,glength
574 : // ij=ij+1
575 : //write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
576 : // enddo
577 41277 : n=n+glength;
578 : }
579 : // Pad last octet with Zeros, if necessary,
580 4 : if ( (iofst%8) != 0) {
581 4 : left=8-(iofst%8);
582 4 : sbit(cpack,&zero,iofst,left);
583 4 : iofst=iofst+left;
584 : }
585 4 : *lcpack=iofst/8;
586 : //
587 4 : free(ifld);
588 4 : free(jfld);
589 4 : free(ifldmiss);
590 4 : free(gref);
591 4 : free(gwidth);
592 4 : free(glen);
593 : //}
594 : //else { // Constant field ( max = min )
595 : // nbits=0;
596 : // *lcpack=0;
597 : // nbitsgref=0;
598 : // ngroups=0;
599 : //}
600 :
601 : //
602 : // Fill in ref value and number of bits in Template 5.2
603 : //
604 4 : mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format
605 4 : idrstmpl[3]=nbitsgref;
606 4 : idrstmpl[4]=0; // original data were reals
607 4 : idrstmpl[5]=1; // general group splitting
608 4 : idrstmpl[9]=ngroups; // Number of groups
609 4 : idrstmpl[10]=ngwidthref; // reference for group widths
610 4 : idrstmpl[11]=nbitsgwidth; // num bits used for group widths
611 4 : idrstmpl[12]=nglenref; // Reference for group lengths
612 4 : idrstmpl[13]=1; // length increment for group lengths
613 4 : idrstmpl[14]=nglenlast; // True length of last group
614 4 : idrstmpl[15]=nbitsglen; // num bits used for group lengths
615 4 : if (idrsnum == 3) {
616 0 : idrstmpl[17]=nbitsd/8; // num bits used for extra spatial
617 : // differencing values
618 : }
619 :
620 : }
|
