#include <math.h>
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include "grib2.h"

static float DoubleToFloatClamp(double val) {
   if (val >= FLT_MAX) return FLT_MAX;
   if (val <= -FLT_MAX) return -FLT_MAX;
   return (float)val;
}

g2int specunpack(unsigned char *cpack,g2int *idrstmpl,g2int ndpts,g2int JJ,
               g2int KK, g2int MM, g2float *fld)
//$$$  SUBPROGRAM DOCUMENTATION BLOCK
//                .      .    .                                       .
// SUBPROGRAM:    specunpack
//   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2000-06-21
//
// ABSTRACT: This subroutine unpacks a spectral data field that was packed
//   using the complex packing algorithm for spherical harmonic data as
//   defined in the GRIB2 documentation,
//   using info from the GRIB2 Data Representation Template 5.51.
//
// PROGRAM HISTORY LOG:
// 2000-06-21  Gilbert
//
// USAGE:    int specunpack(unsigned char *cpack,g2int *idrstmpl,
//                          g2int ndpts,g2int JJ,g2int KK,g2int MM,g2float *fld)
//   INPUT ARGUMENT LIST:
//     cpack    - pointer to the packed data field.
//     idrstmpl - pointer to the array of values for Data Representation
//                Template 5.51
//     ndpts    - The number of data values to unpack (real and imaginary parts)
//     JJ       - J - pentagonal resolution parameter
//     KK       - K - pentagonal resolution parameter
//     MM       - M - pentagonal resolution parameter
//
//   OUTPUT ARGUMENT LIST:
//     fld()    - Contains the unpacked data values.   fld must be allocated
//                with at least ndpts*sizeof(g2float) bytes before
//                calling this routine.
//
// REMARKS: None
//
// ATTRIBUTES:
//   LANGUAGE: C
//   MACHINE:
//
//$$$
{

      g2int  *ifld,j,iofst,nbits;
      g2float  ref,bscale,dscale,*unpk;
      g2float  *pscale,tscale;
      g2int   Js,Ks,Ms,Ts,Ns,Nm,n,m;
      g2int   inc,incu,incp;

      rdieee(idrstmpl+0,&ref,1);
      bscale = DoubleToFloatClamp(int_power(2.0,idrstmpl[1]));
      dscale = DoubleToFloatClamp(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=(g2float *)malloc(ndpts*sizeof(g2float));
         ifld=(g2int *)malloc(ndpts*sizeof(g2int));

         gbits(cpack,G2_UNKNOWN_SIZE,ifld,0,32,0,Ts);
         iofst=32*Ts;
         rdieee(ifld,unpk,Ts);          // read IEEE unpacked floats
         gbits(cpack,G2_UNKNOWN_SIZE,ifld,iofst,nbits,0,ndpts-Ts);  // unpack scaled data
//
//   Calculate Laplacian scaling factors for each possible wave number.
//
         pscale=(g2float *)calloc((JJ+MM+1),sizeof(g2float));
         tscale=(g2float)(idrstmpl[4]*1E-6);
         for (n=Js;n<=JJ+MM;n++)
              pscale[n]=(float)pow((g2float)(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++]=(((g2float)ifld[incp++]*bscale)+ref)*
                            dscale*pscale[n];          // real part
                  fld[inc++]=(((g2float)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 -3;
      }

      return 0;
}