1 : #include <stdio.h>
2 : #include <stdlib.h>
3 : #include <math.h>
4 : #include "grib2.h"
5 :
6 :
7 0 : g2int specunpack(unsigned char *cpack,g2int *idrstmpl,g2int ndpts,g2int JJ,
8 : g2int KK, g2int MM, g2float *fld)
9 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
10 : // . . . .
11 : // SUBPROGRAM: specunpack
12 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2000-06-21
13 : //
14 : // ABSTRACT: This subroutine unpacks a spectral data field that was packed
15 : // using the complex packing algorithm for spherical harmonic data as
16 : // defined in the GRIB2 documention,
17 : // using info from the GRIB2 Data Representation Template 5.51.
18 : //
19 : // PROGRAM HISTORY LOG:
20 : // 2000-06-21 Gilbert
21 : //
22 : // USAGE: int specunpack(unsigned char *cpack,g2int *idrstmpl,
23 : // g2int ndpts,g2int JJ,g2int KK,g2int MM,g2float *fld)
24 : // INPUT ARGUMENT LIST:
25 : // cpack - pointer to the packed data field.
26 : // idrstmpl - pointer to the array of values for Data Representation
27 : // Template 5.51
28 : // ndpts - The number of data values to unpack (real and imaginary parts)
29 : // JJ - J - pentagonal resolution parameter
30 : // KK - K - pentagonal resolution parameter
31 : // MM - M - pentagonal resolution parameter
32 : //
33 : // OUTPUT ARGUMENT LIST:
34 : // fld() - Contains the unpacked data values. fld must be allocated
35 : // with at least ndpts*sizeof(g2float) bytes before
36 : // calling this routine.
37 : //
38 : // REMARKS: None
39 : //
40 : // ATTRIBUTES:
41 : // LANGUAGE: C
42 : // MACHINE:
43 : //
44 : //$$$
45 : {
46 :
47 : g2int *ifld,j,iofst,nbits;
48 : g2float ref,bscale,dscale,*unpk;
49 : g2float *pscale,tscale;
50 : g2int Js,Ks,Ms,Ts,Ns,Nm,n,m;
51 : g2int inc,incu,incp;
52 :
53 0 : rdieee(idrstmpl+0,&ref,1);
54 0 : bscale = int_power(2.0,idrstmpl[1]);
55 0 : dscale = int_power(10.0,-idrstmpl[2]);
56 0 : nbits = idrstmpl[3];
57 0 : Js=idrstmpl[5];
58 0 : Ks=idrstmpl[6];
59 0 : Ms=idrstmpl[7];
60 0 : Ts=idrstmpl[8];
61 :
62 0 : if (idrstmpl[9] == 1) { // unpacked floats are 32-bit IEEE
63 :
64 0 : unpk=(g2float *)malloc(ndpts*sizeof(g2float));
65 0 : ifld=(g2int *)malloc(ndpts*sizeof(g2int));
66 :
67 0 : gbits(cpack,ifld,0,32,0,Ts);
68 0 : iofst=32*Ts;
69 0 : rdieee(ifld,unpk,Ts); // read IEEE unpacked floats
70 0 : gbits(cpack,ifld,iofst,nbits,0,ndpts-Ts); // unpack scaled data
71 : //
72 : // Calculate Laplacian scaling factors for each possible wave number.
73 : //
74 0 : pscale=(g2float *)malloc((JJ+MM+1)*sizeof(g2float));
75 0 : tscale=(g2float)idrstmpl[4]*1E-6;
76 0 : for (n=Js;n<=JJ+MM;n++)
77 0 : pscale[n]=pow((g2float)(n*(n+1)),-tscale);
78 : //
79 : // Assemble spectral coeffs back to original order.
80 : //
81 0 : inc=0;
82 0 : incu=0;
83 0 : incp=0;
84 0 : for (m=0;m<=MM;m++) {
85 0 : Nm=JJ; // triangular or trapezoidal
86 0 : if ( KK == JJ+MM ) Nm=JJ+m; // rhombodial
87 0 : Ns=Js; // triangular or trapezoidal
88 0 : if ( Ks == Js+Ms ) Ns=Js+m; // rhombodial
89 0 : for (n=m;n<=Nm;n++) {
90 0 : if (n<=Ns && m<=Ms) { // grab unpacked value
91 0 : fld[inc++]=unpk[incu++]; // real part
92 0 : fld[inc++]=unpk[incu++]; // imaginary part
93 : }
94 : else { // Calc coeff from packed value
95 0 : fld[inc++]=(((g2float)ifld[incp++]*bscale)+ref)*
96 0 : dscale*pscale[n]; // real part
97 0 : fld[inc++]=(((g2float)ifld[incp++]*bscale)+ref)*
98 0 : dscale*pscale[n]; // imaginary part
99 : }
100 : }
101 : }
102 :
103 0 : free(pscale);
104 0 : free(unpk);
105 0 : free(ifld);
106 :
107 : }
108 : else {
109 0 : printf("specunpack: Cannot handle 64 or 128-bit floats.\n");
110 0 : for (j=0;j<ndpts;j++) fld[j]=0.0;
111 0 : return -3;
112 : }
113 :
114 0 : return 0;
115 : }
|
