specunpack.c

Go to the documentation of this file.

 
#include "grib2_int.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
 
g2int
specunpack(unsigned char *cpack, g2int *idrstmpl, g2int ndpts, g2int JJ,
           g2int KK, g2int MM, float *fld)
{
    g2int *ifld, j, iofst, nbits;
    float ref, bscale, dscale, *unpk;
    float *pscale, tscale;
    g2int Js, Ks, Ms, Ts, Ns, Nm, n, m;
    g2int inc, incu, incp;
 
    rdieee(idrstmpl + 0, &ref, 1);
    bscale = int_power(2.0, idrstmpl[1]);
    dscale = int_power(10.0, -idrstmpl[2]);
    nbits = idrstmpl[3];
    Js = idrstmpl[5];
    Ks = idrstmpl[6];
    Ms = idrstmpl[7];
    Ts = idrstmpl[8];
 
    if (idrstmpl[9] == 1)
    { /* unpacked floats are 32-bit IEEE */
 
        unpk = malloc(ndpts * sizeof(float));
        ifld = malloc(ndpts * sizeof(g2int));
 
        gbits(cpack, ifld, 0, 32, 0, Ts);
        iofst = 32 * Ts;
        rdieee(ifld, unpk, Ts);                          /* read IEEE unpacked floats */
        gbits(cpack, ifld, iofst, nbits, 0, ndpts - Ts); /* unpack scaled data */
 
        /* Calculate Laplacian scaling factors for each possible wave
         * number. */
        pscale = malloc((JJ + MM + 1) * sizeof(float));
        tscale = idrstmpl[4] * 1E-6;
        for (n = Js; n <= JJ + MM; n++)
            pscale[n] = pow((float)(n * (n + 1)), -tscale);
 
        /* Assemble spectral coeffs back to original order. */
        inc = 0;
        incu = 0;
        incp = 0;
        for (m = 0; m <= MM; m++)
        {
            Nm = JJ; /* triangular or trapezoidal */
            if (KK == JJ + MM)
                Nm = JJ + m; /* rhombodial */
            Ns = Js;         /* triangular or trapezoidal */
            if (Ks == Js + Ms)
                Ns = Js + m; /* rhombodial */
            for (n = m; n <= Nm; n++)
            {
                if (n <= Ns && m <= Ms)
                {                              /* grab unpacked value */
                    fld[inc++] = unpk[incu++]; /* real part */
                    fld[inc++] = unpk[incu++]; /* imaginary part */
                }
                else
                { /* Calc coeff from packed value */
                    fld[inc++] = (((float)ifld[incp++] * bscale) + ref) *
                                 dscale * pscale[n]; /* real part */
                    fld[inc++] = (((float)ifld[incp++] * bscale) + ref) *
                                 dscale * pscale[n]; /* imaginary part */
                }
            }
        }
 
        free(pscale);
        free(unpk);
        free(ifld);
    }
    else
    {
        printf("specunpack: Cannot handle 64 or 128-bit floats.\n");
        for (j = 0; j < ndpts; j++)
            fld[j] = 0.0;
        return G2_SPECUNPACK_TYPE;
    }
 
    return G2_NO_ERROR;
}