1 /*
2 * Copyright 2005-2018 ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence Version 2.0
5 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
8 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
9 */
10
11 #include "grib_api_internal.h"
12 #include <math.h>
13 /*
14 This is used by make_class.pl
15
16 START_CLASS_DEF
17 CLASS = accessor
18 SUPER = grib_accessor_class_data_simple_packing
19 IMPLEMENTS = unpack_double
20 IMPLEMENTS = pack_double
21 IMPLEMENTS = value_count
22 IMPLEMENTS = init
23 MEMBERS= const char* GRIBEX_sh_bug_present
24 MEMBERS= const char* ieee_floats
25 MEMBERS= const char* laplacianOperatorIsSet
26 MEMBERS= const char* laplacianOperator
27 MEMBERS= const char* sub_j
28 MEMBERS= const char* sub_k
29 MEMBERS= const char* sub_m
30 MEMBERS= const char* pen_j
31 MEMBERS= const char* pen_k
32 MEMBERS= const char* pen_m
33 END_CLASS_DEF
34
35 */
36
37
38 /* START_CLASS_IMP */
39
40 /*
41
42 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
43 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
44 or edit "accessor.class" and rerun ./make_class.pl
45
46 */
47
48 static int pack_double(grib_accessor*, const double* val,size_t *len);
49 static int unpack_double(grib_accessor*, double* val,size_t *len);
50 static int value_count(grib_accessor*,long*);
51 static void init(grib_accessor*,const long, grib_arguments* );
52 static void init_class(grib_accessor_class*);
53
54 typedef struct grib_accessor_data_complex_packing {
55 grib_accessor att;
56 /* Members defined in gen */
57 /* Members defined in values */
58 int carg;
59 const char* seclen;
60 const char* offsetdata;
61 const char* offsetsection;
62 int dirty;
63 /* Members defined in data_simple_packing */
64 int edition;
65 const char* units_factor;
66 const char* units_bias;
67 const char* changing_precision;
68 const char* number_of_values;
69 const char* bits_per_value;
70 const char* reference_value;
71 const char* binary_scale_factor;
72 const char* decimal_scale_factor;
73 /* Members defined in data_complex_packing */
74 const char* GRIBEX_sh_bug_present;
75 const char* ieee_floats;
76 const char* laplacianOperatorIsSet;
77 const char* laplacianOperator;
78 const char* sub_j;
79 const char* sub_k;
80 const char* sub_m;
81 const char* pen_j;
82 const char* pen_k;
83 const char* pen_m;
84 } grib_accessor_data_complex_packing;
85
86 extern grib_accessor_class* grib_accessor_class_data_simple_packing;
87
88 static grib_accessor_class _grib_accessor_class_data_complex_packing = {
89 &grib_accessor_class_data_simple_packing, /* super */
90 "data_complex_packing", /* name */
91 sizeof(grib_accessor_data_complex_packing), /* size */
92 0, /* inited */
93 &init_class, /* init_class */
94 &init, /* init */
95 0, /* post_init */
96 0, /* free mem */
97 0, /* describes himself */
98 0, /* get length of section */
99 0, /* get length of string */
100 &value_count, /* get number of values */
101 0, /* get number of bytes */
102 0, /* get offset to bytes */
103 0, /* get native type */
104 0, /* get sub_section */
105 0, /* grib_pack procedures long */
106 0, /* grib_pack procedures long */
107 0, /* grib_pack procedures long */
108 0, /* grib_unpack procedures long */
109 &pack_double, /* grib_pack procedures double */
110 &unpack_double, /* grib_unpack procedures double */
111 0, /* grib_pack procedures string */
112 0, /* grib_unpack procedures string */
113 0, /* grib_pack procedures bytes */
114 0, /* grib_unpack procedures bytes */
115 0, /* pack_expression */
116 0, /* notify_change */
117 0, /* update_size */
118 0, /* preferred_size */
119 0, /* resize */
120 0, /* nearest_smaller_value */
121 0, /* next accessor */
122 0, /* compare vs. another accessor */
123 0, /* unpack only ith value */
124 0, /* unpack a subarray */
125 0, /* clear */
126 };
127
128
129 grib_accessor_class* grib_accessor_class_data_complex_packing = &_grib_accessor_class_data_complex_packing;
130
131
init_class(grib_accessor_class * c)132 static void init_class(grib_accessor_class* c)
133 {
134 c->dump = (*(c->super))->dump;
135 c->next_offset = (*(c->super))->next_offset;
136 c->string_length = (*(c->super))->string_length;
137 c->byte_count = (*(c->super))->byte_count;
138 c->byte_offset = (*(c->super))->byte_offset;
139 c->get_native_type = (*(c->super))->get_native_type;
140 c->sub_section = (*(c->super))->sub_section;
141 c->pack_missing = (*(c->super))->pack_missing;
142 c->is_missing = (*(c->super))->is_missing;
143 c->pack_long = (*(c->super))->pack_long;
144 c->unpack_long = (*(c->super))->unpack_long;
145 c->pack_string = (*(c->super))->pack_string;
146 c->unpack_string = (*(c->super))->unpack_string;
147 c->pack_bytes = (*(c->super))->pack_bytes;
148 c->unpack_bytes = (*(c->super))->unpack_bytes;
149 c->pack_expression = (*(c->super))->pack_expression;
150 c->notify_change = (*(c->super))->notify_change;
151 c->update_size = (*(c->super))->update_size;
152 c->preferred_size = (*(c->super))->preferred_size;
153 c->resize = (*(c->super))->resize;
154 c->nearest_smaller_value = (*(c->super))->nearest_smaller_value;
155 c->next = (*(c->super))->next;
156 c->compare = (*(c->super))->compare;
157 c->unpack_double_element = (*(c->super))->unpack_double_element;
158 c->unpack_double_subarray = (*(c->super))->unpack_double_subarray;
159 c->clear = (*(c->super))->clear;
160 }
161
162 /* END_CLASS_IMP */
163
164 typedef unsigned long (*encode_float_proc)(double);
165 typedef double (*decode_float_proc)(unsigned long);
166
init(grib_accessor * a,const long v,grib_arguments * args)167 static void init(grib_accessor* a,const long v, grib_arguments* args)
168 {
169 grib_accessor_data_complex_packing *self =(grib_accessor_data_complex_packing*)a;
170
171 self->GRIBEX_sh_bug_present = grib_arguments_get_name(a->parent->h,args,self->carg++);
172 self->ieee_floats = grib_arguments_get_name(a->parent->h,args,self->carg++);
173 self->laplacianOperatorIsSet = grib_arguments_get_name(a->parent->h,args,self->carg++);
174 self->laplacianOperator = grib_arguments_get_name(a->parent->h,args,self->carg++);
175 self->sub_j = grib_arguments_get_name(a->parent->h,args,self->carg++);
176 self->sub_k = grib_arguments_get_name(a->parent->h,args,self->carg++);
177 self->sub_m = grib_arguments_get_name(a->parent->h,args,self->carg++);
178 self->pen_j = grib_arguments_get_name(a->parent->h,args,self->carg++);
179 self->pen_k = grib_arguments_get_name(a->parent->h,args,self->carg++);
180 self->pen_m = grib_arguments_get_name(a->parent->h,args,self->carg++);
181
182 a->flags |= GRIB_ACCESSOR_FLAG_DATA;
183 }
184
value_count(grib_accessor * a,long * count)185 static int value_count(grib_accessor* a,long* count)
186 {
187 grib_accessor_data_complex_packing *self =(grib_accessor_data_complex_packing*)a;
188 int ret = 0;
189
190 long pen_j= 0;
191 long pen_k= 0;
192 long pen_m= 0;
193
194 *count=0;
195
196 if(a->length == 0) return 0;
197
198 if((ret = grib_get_long_internal(a->parent->h,self->pen_j,&pen_j)) != GRIB_SUCCESS)
199 return ret;
200 if((ret = grib_get_long_internal(a->parent->h,self->pen_k,&pen_k)) != GRIB_SUCCESS)
201 return ret;
202 if((ret = grib_get_long_internal(a->parent->h,self->pen_m,&pen_m)) != GRIB_SUCCESS)
203 return ret;
204
205 if (pen_j != pen_k || pen_j!=pen_m ) {
206 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,"pen_j=%ld, pen_k=%ld, pen_m=%ld\n",pen_j,pen_k,pen_m);
207 Assert ((pen_j == pen_k) && (pen_j == pen_m));
208 }
209 *count=(pen_j+1)*(pen_j+2);
210
211 return ret;
212 }
213
214
unpack_double(grib_accessor * a,double * val,size_t * len)215 static int unpack_double(grib_accessor* a, double* val, size_t *len)
216 {
217 grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a;
218
219 size_t i = 0;
220 int ret = GRIB_SUCCESS;
221 long hcount = 0;
222 long lcount = 0;
223 long hpos = 0;
224 long lup = 0;
225 long mmax = 0;
226 long n_vals = 0;
227 double *scals = NULL;
228 double *pscals=NULL,*pval=NULL;
229
230 double s = 0;
231 double d = 0;
232 double laplacianOperator = 0;
233 unsigned char* buf = NULL;
234 unsigned char* hres = NULL;
235 unsigned char* lres = NULL;
236 unsigned long packed_offset;
237 long lpos = 0;
238
239 long maxv = 0;
240 long GRIBEX_sh_bug_present =0;
241 long ieee_floats = 0;
242
243 long offsetdata = 0;
244 long bits_per_value = 0;
245 double reference_value = 0;
246 long binary_scale_factor = 0;
247 long decimal_scale_factor = 0;
248
249 long sub_j= 0;
250 long sub_k= 0;
251 long sub_m= 0;
252 long pen_j= 0;
253 long pen_k= 0;
254 long pen_m= 0;
255
256 double operat= 0;
257 int bytes;
258 int err=0;
259
260 decode_float_proc decode_float = NULL;
261
262 err=grib_value_count(a,&n_vals);
263 if (err) return err;
264
265 if(*len < n_vals){
266 *len = n_vals;
267 return GRIB_ARRAY_TOO_SMALL;
268 }
269
270 if((ret = grib_get_long_internal(a->parent->h,self->offsetdata,&offsetdata))
271 != GRIB_SUCCESS) return ret;
272 if((ret = grib_get_long_internal(a->parent->h,self->bits_per_value,&bits_per_value))
273 != GRIB_SUCCESS) return ret;
274 if((ret = grib_get_double_internal(a->parent->h,self->reference_value,&reference_value))
275 != GRIB_SUCCESS) return ret;
276 if((ret = grib_get_long_internal(a->parent->h,self->binary_scale_factor,&binary_scale_factor))
277 != GRIB_SUCCESS) return ret;
278
279 if((ret = grib_get_long_internal(a->parent->h,self->decimal_scale_factor,&decimal_scale_factor))
280 != GRIB_SUCCESS)
281 return ret;
282
283 if((ret = grib_get_long_internal(a->parent->h,self->GRIBEX_sh_bug_present,&GRIBEX_sh_bug_present))
284 != GRIB_SUCCESS)
285 return ret;
286
287 if((ret = grib_get_long_internal(a->parent->h,self->ieee_floats,&ieee_floats)) != GRIB_SUCCESS)
288 return ret;
289
290 if((ret = grib_get_double_internal(a->parent->h,self->laplacianOperator,&laplacianOperator))
291 != GRIB_SUCCESS)
292 return ret;
293 if((ret = grib_get_long_internal(a->parent->h,self->sub_j,&sub_j)) != GRIB_SUCCESS)
294 return ret;
295 if((ret = grib_get_long_internal(a->parent->h,self->sub_k,&sub_k)) != GRIB_SUCCESS)
296 return ret;
297 if((ret = grib_get_long_internal(a->parent->h,self->sub_m,&sub_m)) != GRIB_SUCCESS)
298 return ret;
299 if((ret = grib_get_long_internal(a->parent->h,self->pen_j,&pen_j)) != GRIB_SUCCESS)
300 return ret;
301 if((ret = grib_get_long_internal(a->parent->h,self->pen_k,&pen_k)) != GRIB_SUCCESS)
302 return ret;
303 if((ret = grib_get_long_internal(a->parent->h,self->pen_m,&pen_m)) != GRIB_SUCCESS)
304 return ret;
305
306 self->dirty=0;
307
308 switch (ieee_floats) {
309 case 0:
310 decode_float=grib_long_to_ibm;
311 bytes=4;
312 break;
313 case 1:
314 decode_float=grib_long_to_ieee;
315 bytes=4;
316 break;
317 case 2:
318 decode_float=grib_long_to_ieee64;
319 bytes=8;
320 break;
321 default:
322 return GRIB_NOT_IMPLEMENTED;
323 }
324
325 Assert (sub_j == sub_k);
326 Assert (sub_j == sub_m);
327 Assert (pen_j == pen_k);
328 Assert (pen_j == pen_m);
329
330 buf = (unsigned char*)a->parent->h->buffer->data;
331
332 maxv = pen_j+1;
333
334 buf += grib_byte_offset(a);
335 hres = buf;
336 lres = buf;
337
338 if (pen_j == sub_j) {
339 n_vals = (pen_j+1)*(pen_j+2);
340 d = grib_power(-decimal_scale_factor,10) ;
341 grib_ieee_decode_array(a->parent->h->context,buf,n_vals,bytes,val);
342 if (d) {
343 for (i=0;i<n_vals;i++) val[i]*=d;
344 }
345 return 0;
346 }
347
348 packed_offset = grib_byte_offset(a) + 4*(sub_k+1)*(sub_k+2);
349
350 lpos = 8*(packed_offset-offsetdata);
351
352 s = grib_power(binary_scale_factor,2);
353 d = grib_power(-decimal_scale_factor,10) ;
354
355 scals = (double*)grib_context_malloc(a->parent->h->context,maxv*sizeof(double));
356 Assert(scals);
357
358 scals[0] = 0;
359 for(i=1;i<maxv;i++){
360 operat = pow(i*(i+1),laplacianOperator);
361 if(operat != 0)
362 scals[i] = (1.0/operat);
363 else{
364 grib_context_log(a->parent->h->context,GRIB_LOG_WARNING,
365 "COMPLEX_PACKING : problem with operator div by zero at index %d of %d \n",
366 i , maxv);
367 scals[i] = 0;
368 }
369 }
370
371 /*
372 printf("UNPACKING LAPLACE=%.20f\n",laplacianOperator);
373
374 printf("packed offset=%ld\n",packed_offset);
375 for(i=0;i<maxv;i++)
376 printf("scals[%d]=%g\n",i,scals[i]);*/
377
378 i=0;
379
380 while(maxv>0)
381 {
382 lup=mmax;
383 if(sub_k>=0)
384 {
385 for(hcount=0;hcount<sub_k+1;hcount++)
386 {
387 val[i++] = decode_float(grib_decode_unsigned_long(hres,&hpos,32))*d;
388 val[i++] = decode_float(grib_decode_unsigned_long(hres,&hpos,32))*d;
389
390 if (GRIBEX_sh_bug_present && hcount==sub_k){
391 /* bug in ecmwf data, last row (K+1)is scaled but should not */
392 val[i-2] *= scals[lup];
393 val[i-1] *= scals[lup];
394 }
395 lup++;
396 }
397 sub_k--;
398 }
399
400 pscals=scals+lup;
401 pval=val+i;
402 #if FAST_BIG_ENDIAN
403 grib_decode_double_array_complex(lres,
404 &lpos,bits_per_value,
405 reference_value,s,pscals,(maxv-hcount)*2,pval);
406 i+=(maxv-hcount)*2;
407 #else
408 (void)pscals; /* suppress gcc warning */
409 (void)pval; /* suppress gcc warning */
410 for(lcount=hcount; lcount < maxv ; lcount++)
411 {
412 val[i++] = (double) ((grib_decode_unsigned_long(lres, &lpos,
413 bits_per_value)*s)+reference_value)*scals[lup];
414 val[i++] = (double) ((grib_decode_unsigned_long(lres, &lpos,
415 bits_per_value)*s)+reference_value)*scals[lup];
416 lup++;
417 }
418 #endif
419
420 maxv--;
421 hcount=0;
422 mmax++;
423 }
424
425 Assert(*len >= i);
426 *len = i;
427
428 if(d != 1) {
429 for(i=0;i<*len;i++)
430 val[i++] *= d;
431 }
432
433 grib_context_free(a->parent->h->context,scals);
434
435 return ret;
436 }
437
438
439 #define MAXVAL(a,b) a>b?a:b
440
calculate_pfactor(grib_context * ctx,const double * spectralField,long fieldTruncation,long subsetTruncation)441 static double calculate_pfactor(grib_context *ctx,const double* spectralField, long fieldTruncation, long subsetTruncation)
442 {
443 /*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/
444 long loop, index, m, n = 0;
445 double pFactor, zeps = 1.0e-15;
446 long ismin = (subsetTruncation+1), ismax = (fieldTruncation+1);
447 double* weights, range, * norms;
448 double weightedSumOverX = 0.0, weightedSumOverY = 0.0, sumOfWeights = 0.0, x, y;
449 double numerator = 0.0, denominator = 0.0, slope;
450
451 /* Catch corner case. See GRIB-172 */
452 if (ismax-ismin <= 1) {
453 return 1; /* any value will do! we cannot do linear fit on a single point! */
454 }
455
456 /*
457 * Setup the weights
458 */
459 range = (double) (ismax - ismin +1);
460 weights = (double*) grib_context_malloc(ctx,(ismax+1)*sizeof(double));
461 for( loop = ismin; loop <= ismax; loop++ )
462 weights[loop] = range / (double) (loop-ismin+1);
463
464 /*
465 * Compute norms
466 * Handle values 2 at a time (real and imaginary parts).
467 */
468 norms = (double*) grib_context_malloc(ctx,(ismax+1)*sizeof(double));
469 for( loop = 0; loop < ismax+1; loop++ ) norms[loop] = 0.0;
470
471 /*
472 * Form norms for the rows which contain part of the unscaled subset.
473 */
474
475 index = -2;
476 for( m = 0; m < subsetTruncation; m++ ) {
477 for( n = m; n <= fieldTruncation; n++ ) {
478 index += 2;
479 if( n >= subsetTruncation ) {
480 norms[n] = MAXVAL(norms[n],fabs(spectralField[index]));
481 norms[n] = MAXVAL(norms[n],fabs(spectralField[index+1]));
482 }
483 }
484 }
485
486 /*
487 * Form norms for the rows which do not contain part of the unscaled subset.
488 */
489
490 for( m = subsetTruncation; m <= fieldTruncation; m++ ) {
491 for( n = m; n <= fieldTruncation; n++ ) {
492 index += 2;
493 norms[n] = MAXVAL(norms[n],fabs(spectralField[index]));
494 norms[n] = MAXVAL(norms[n],fabs(spectralField[index+1]));
495 }
496 }
497
498 /*
499 * Ensure the norms have a value which is not too small in case of
500 * problems with math functions (e.g. LOG).
501 */
502 for( loop = ismin; loop <= ismax; loop++ ) {
503 norms[loop] = MAXVAL(norms[loop],zeps);
504 if( norms[loop] == zeps ) weights[loop] = 100.0 * zeps;
505 }
506
507 /*
508 * Do linear fit to find the slope
509 */
510 for( loop = ismin; loop <= ismax; loop++ ) {
511 x = log( (double) (loop*(loop+1)) );
512 Assert( norms[loop] > 0 );
513 y = log( norms[loop] );
514 weightedSumOverX = weightedSumOverX + x * weights[loop];
515 weightedSumOverY = weightedSumOverY + y * weights[loop];
516 sumOfWeights = sumOfWeights + weights[loop];
517 }
518 weightedSumOverX = weightedSumOverX / sumOfWeights;
519 weightedSumOverY = weightedSumOverY / sumOfWeights;
520
521 /*
522 * Perform a least square fit for the equation
523 */
524
525 for( loop = ismin; loop <= ismax; loop++ ) {
526 x = log( (double)(loop*(loop+1)) );
527 y = log( norms[loop] );
528 numerator =
529 numerator + weights[loop] * (y-weightedSumOverY) * (x-weightedSumOverX);
530 denominator =
531 denominator + weights[loop] * ((x-weightedSumOverX) * (x-weightedSumOverX));
532 }
533 slope = numerator / denominator;
534
535 grib_context_free(ctx,weights);
536 grib_context_free(ctx,norms);
537
538 pFactor = -slope;
539 if( pFactor < -9999.9 ) pFactor = -9999.9;
540 if( pFactor > 9999.9 ) pFactor = 9999.9;
541 return pFactor;
542 }
543
544
pack_double(grib_accessor * a,const double * val,size_t * len)545 static int pack_double(grib_accessor* a, const double* val, size_t *len) {
546
547 grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a;
548
549 size_t i = 0;
550 int ret = GRIB_SUCCESS;
551 long hcount = 0;
552 long lcount = 0;
553 long hpos = 0;
554 long lup = 0;
555 long mmax = 0;
556 long n_vals = 0;
557 double *scals = NULL;
558
559 double s = 0;
560 double d = 0;
561
562 unsigned char* buf = NULL;
563
564 size_t buflen = 0;
565 size_t hsize = 0;
566 size_t lsize = 0;
567
568 unsigned char* hres = NULL;
569 unsigned char* lres = NULL;
570
571 long lpos = 0;
572 long maxv = 0;
573
574 long offsetdata = 0;
575 long bits_per_value = 0;
576 double reference_value = 0;
577 long binary_scale_factor = 0;
578 long decimal_scale_factor = 0;
579 long laplacianOperatorIsSet = 0;
580
581 double laplacianOperator = 0;
582 long sub_j= 0;
583 long sub_k= 0;
584 long sub_m= 0;
585 long pen_j= 0;
586 long pen_k= 0;
587 long pen_m= 0;
588 long GRIBEX_sh_bug_present =0;
589 long ieee_floats =0;
590 double min = 0;
591 double max = 0;
592 double current_val = 0;
593 short mixmax_unset = 0;
594 int bytes;
595
596 encode_float_proc encode_float = NULL;
597
598 if (*len ==0) return GRIB_NO_VALUES;
599
600 if((ret = grib_get_long_internal(a->parent->h,self->offsetdata,&offsetdata)) != GRIB_SUCCESS)
601 return ret;
602 if((ret = grib_get_long_internal(a->parent->h,self->bits_per_value,&bits_per_value)) != GRIB_SUCCESS)
603 return ret;
604
605 if((ret = grib_get_long_internal(a->parent->h,self->decimal_scale_factor,&decimal_scale_factor)) != GRIB_SUCCESS)
606 return ret;
607
608 if((ret = grib_get_long_internal(a->parent->h,self->GRIBEX_sh_bug_present,&GRIBEX_sh_bug_present)) != GRIB_SUCCESS)
609 return ret;
610
611 if((ret = grib_get_long_internal(a->parent->h,self->ieee_floats,&ieee_floats)) != GRIB_SUCCESS)
612 return ret;
613
614 if((ret = grib_get_long_internal(a->parent->h,self->laplacianOperatorIsSet,&laplacianOperatorIsSet)) != GRIB_SUCCESS)
615 return ret;
616 if((ret = grib_get_double_internal(a->parent->h,self->laplacianOperator,&laplacianOperator)) != GRIB_SUCCESS)
617 return ret;
618
619 if((ret = grib_get_long_internal(a->parent->h,self->sub_j,&sub_j)) != GRIB_SUCCESS)
620 return ret;
621 if((ret = grib_get_long_internal(a->parent->h,self->sub_k,&sub_k)) != GRIB_SUCCESS)
622 return ret;
623 if((ret = grib_get_long_internal(a->parent->h,self->sub_m,&sub_m)) != GRIB_SUCCESS)
624 return ret;
625 if((ret = grib_get_long_internal(a->parent->h,self->pen_j,&pen_j)) != GRIB_SUCCESS)
626 return ret;
627 if((ret = grib_get_long_internal(a->parent->h,self->pen_k,&pen_k)) != GRIB_SUCCESS)
628 return ret;
629 if((ret = grib_get_long_internal(a->parent->h,self->pen_m,&pen_m)) != GRIB_SUCCESS)
630 return ret;
631
632 self->dirty=1;
633
634
635 switch (ieee_floats) {
636 case 0:
637 encode_float =grib_ibm_to_long;
638 bytes=4;
639 break;
640 case 1:
641 encode_float =grib_ieee_to_long;
642 bytes=4;
643 break;
644 case 2:
645 encode_float =grib_ieee64_to_long;
646 bytes=8;
647 break;
648 default:
649 return GRIB_NOT_IMPLEMENTED;
650 }
651
652 Assert (sub_j == sub_k); Assert( sub_j == sub_m);
653 Assert (pen_j == pen_k); Assert( pen_j == pen_m);
654
655 n_vals = (pen_j+1)*(pen_j+2);
656 d = grib_power(decimal_scale_factor,10) ;
657
658 if(*len != n_vals){
659 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,"COMPLEX_PACKING : wrong number of values, expected %d - got %d",n_vals,*len);
660 return GRIB_INTERNAL_ERROR;
661 }
662
663 if (pen_j == sub_j) {
664 double* values;
665 if (d) {
666 values=(double*)grib_context_malloc_clear(a->parent->h->context,sizeof(double)*n_vals);
667 for (i=0;i<n_vals;i++) values[i]=val[i]*d;
668 } else {
669 values=(double*)val;
670 }
671 buflen=n_vals*bytes;
672 buf = (unsigned char*)grib_context_malloc_clear(a->parent->h->context,buflen);
673 grib_ieee_encode_array(a->parent->h->context,values,n_vals,bytes,buf);
674 if (d) grib_context_free(a->parent->h->context,values);
675 grib_buffer_replace(a, buf, buflen,1,1);
676 grib_context_free(a->parent->h->context,buf);
677 return 0;
678 }
679
680 if(!laplacianOperatorIsSet) {
681 laplacianOperator = calculate_pfactor(a->parent->h->context,val,pen_j,sub_j);
682 if((ret = grib_set_double_internal(a->parent->h,self->laplacianOperator,laplacianOperator))
683 != GRIB_SUCCESS) return ret;
684 grib_get_double_internal(a->parent->h,self->laplacianOperator,&laplacianOperator);
685 }
686
687 /*
688 printf("PACKING LAPLACE set=%ld value=%.20f\n",laplacianOperatorIsSet,laplacianOperator);
689 */
690
691 hsize = 4*(sub_k+1)*(sub_k+2);
692 lsize = ((n_vals - ((sub_k+1)*(sub_k+2)))*bits_per_value)/8;
693
694 buflen = hsize+lsize;
695
696 buf = (unsigned char*)grib_context_malloc(a->parent->h->context,buflen);
697 hres = buf;
698 lres = buf+hsize;
699
700 maxv = pen_j+1;
701
702 lpos = 0;
703 hpos = 0;
704
705 scals = (double*) grib_context_malloc(a->parent->h->context,maxv*sizeof(double));
706 Assert(scals);
707
708 scals[0] =0;
709 for(i=1;i<maxv;i++)
710 scals[i] = ((double)pow(i*(i+1),laplacianOperator));
711
712
713 i=0;
714
715 mmax = 0;
716 maxv = pen_j+1;
717 i=0;
718 lcount=0;
719 hcount=0;
720 sub_k = sub_j;
721
722 while(maxv>0)
723 {
724 lup=mmax;
725
726 if(sub_k>=0)
727 {
728 i += 2*(sub_k+1);
729 lup += sub_k+1 ;
730 hcount += sub_k+1 ;
731 sub_k--;
732 }
733
734 for(lcount=hcount; lcount < maxv ; lcount++)
735 {
736 current_val = ((val[i++]*d) * scals[lup]);
737 if(mixmax_unset == 0){
738 max = current_val;
739 min = current_val;
740 mixmax_unset = 1;
741 }
742
743 if(current_val > max) max = current_val;
744 if(current_val < min) min = current_val;
745
746 current_val = ((val[i++]*d) * scals[lup]);
747 if(current_val > max) max = current_val;
748 if(current_val < min) min = current_val;
749
750 lup++;
751 }
752 maxv--;
753 hcount=0;
754 mmax++;
755 }
756
757 if (grib_get_nearest_smaller_value(a->parent->h,self->reference_value,min,&reference_value)
758 !=GRIB_SUCCESS) {
759 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,
760 "unable to find nearest_smaller_value of %g for %s",min,self->reference_value);
761 exit(GRIB_INTERNAL_ERROR);
762 }
763 binary_scale_factor = grib_get_binary_scale_fact(max,reference_value,bits_per_value,&ret);
764
765 if (ret==GRIB_UNDERFLOW) {
766 d=0;
767 binary_scale_factor = 0;
768 reference_value=0;
769
770 }
771 s = grib_power(-binary_scale_factor,2);
772
773 /* printf("D : %.30f\n",d); */
774
775 i=0;
776
777 mmax = 0;
778 maxv = pen_j+1;
779 i=0;
780 lcount=0;
781 hcount=0;
782 sub_k = sub_j;
783
784 while(maxv>0)
785 {
786 lup=mmax;
787
788 if(sub_k>=0)
789 {
790 for(hcount=0;hcount<sub_k+1;hcount++)
791 {
792 if ( GRIBEX_sh_bug_present && hcount==sub_k ) {
793 /* _test(val[i]*d*scals[lup],1); */
794 grib_encode_unsigned_long(hres, encode_float((val[i++]*d)*scals[lup]) , &hpos, 32);
795 /* _test(val[i]*d*scals[lup],1); */
796 grib_encode_unsigned_long(hres, encode_float((val[i++]*d)*scals[lup]) , &hpos, 32);
797 }else{
798
799 /* _test(val[i]*d,0); */
800
801 grib_encode_unsigned_long(hres, encode_float(val[i++]*d) , &hpos, 32);
802 /* _test(val[i]*d,0); */
803 grib_encode_unsigned_long(hres, encode_float(val[i++]*d) , &hpos, 32);
804 }
805 lup++;
806 }
807 sub_k--;
808 }
809
810 #if FAST_BIG_ENDIAN
811 grib_encode_double_array_complex((maxv-hcount)*2,&(val[i]),bits_per_value,reference_value,&(scals[lup]),d,s,lres,&lpos);
812 i+=(maxv-hcount)*2;
813 #else
814 if (bits_per_value % 8) {
815 for(lcount=hcount; lcount < maxv ; lcount++)
816 {
817 current_val = (((((val[i++]*d) * scals[lup])-reference_value)*s)+0.5);
818 if(current_val < 0)
819 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,
820 "COMPLEX_PACKING : negative coput before packing (%g)", current_val);
821 grib_encode_unsigned_longb(lres, current_val, &lpos, bits_per_value);
822
823 current_val = (((((val[i++]*d) * scals[lup])-reference_value)*s)+0.5);
824 if(current_val < 0)
825 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,
826 "COMPLEX_PACKING : negative coput before packing (%g)", current_val);
827 grib_encode_unsigned_longb(lres, current_val, &lpos, bits_per_value);
828 lup++;
829 }
830 } else {
831 for(lcount=hcount; lcount < maxv ; lcount++)
832 {
833 current_val = (((((val[i++]*d) * scals[lup])-reference_value)*s)+0.5);
834 if(current_val < 0)
835 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,
836 "COMPLEX_PACKING : negative coput before packing (%g)", current_val);
837 grib_encode_unsigned_long(lres, current_val, &lpos, bits_per_value);
838
839 current_val = (((((val[i++]*d) * scals[lup])-reference_value)*s)+0.5);
840 if(current_val < 0)
841 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,
842 "COMPLEX_PACKING : negative coput before packing (%g)", current_val);
843 grib_encode_unsigned_long(lres, current_val, &lpos, bits_per_value);
844 lup++;
845 }
846 }
847 #endif
848
849 maxv--;
850 hcount=0;
851 mmax++;
852 }
853
854 if(((hpos/8) != hsize) &&((lpos/8) != lsize))
855 {
856 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,
857 "COMPLEX_PACKING : Mismatch in packing between high resolution and low resolution part");
858 grib_context_free(a->parent->h->context,buf);
859 grib_context_free(a->parent->h->context,scals);
860 return GRIB_INTERNAL_ERROR;
861 }
862
863 buflen = ((hpos + lpos)/8);
864
865 if((ret = grib_set_double_internal(a->parent->h,self->reference_value, reference_value)) != GRIB_SUCCESS)
866 return ret;
867 {
868 /* Make sure we can decode it again */
869 double ref = 1e-100;
870 grib_get_double_internal(a->parent->h,self->reference_value,&ref);
871 Assert(ref == reference_value);
872 }
873
874 if((ret = grib_set_long_internal(a->parent->h,self->binary_scale_factor, binary_scale_factor)) != GRIB_SUCCESS)
875 return ret;
876
877 grib_buffer_replace(a, buf, buflen,1,1);
878 grib_context_free(a->parent->h->context,buf);
879 grib_context_free(a->parent->h->context,scals);
880
881 return ret;
882
883 }
884