#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 jpcunpack(unsigned char *cpack,g2int len,g2int *idrstmpl,g2int ndpts,
                g2float **fld)
//$$$  SUBPROGRAM DOCUMENTATION BLOCK
//                .      .    .                                       .
// SUBPROGRAM:    jpcunpack
//   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2003-08-27
//
// ABSTRACT: This subroutine unpacks a data field that was packed into a
//   JPEG2000 code stream
//   using info from the GRIB2 Data Representation Template 5.40 or 5.40000.
//
// PROGRAM HISTORY LOG:
// 2003-08-27  Gilbert
//
// USAGE:    jpcunpack(unsigned char *cpack,g2int len,g2int *idrstmpl,g2int ndpts,
//                     g2float *fld)
//   INPUT ARGUMENT LIST:
//     cpack    - The packed data field (character*1 array)
//     len      - length of packed field cpack().
//     idrstmpl - Pointer to array of values for Data Representation
//                Template 5.40 or 5.40000
//     ndpts    - The number of data values to unpack
//
//   OUTPUT ARGUMENT LIST:
//     fld[]    - Contains the unpacked data values
//
// REMARKS: None
//
// ATTRIBUTES:
//   LANGUAGE: C
//   MACHINE:  IBM SP
//
//$$$
{

      g2int  *ifld;
      g2int  j,nbits, iret;
      g2float  ref,bscale,dscale;

      rdieee(idrstmpl+0,&ref,1);
      bscale = DoubleToFloatClamp(int_power(2.0,idrstmpl[1]));
      dscale = DoubleToFloatClamp(int_power(10.0,-idrstmpl[2]));
      nbits = idrstmpl[3];
//
//  if nbits equals 0, we have a constant field where the reference value
//  is the data value at each gridpoint
//
      *fld = 0;
      if (nbits != 0) {

         ifld = NULL;
         iret= (g2int) dec_jpeg2000(cpack,len,&ifld,ndpts);
         if( iret != 0 )
         {
             free(ifld);
             return -1;
         }
         *fld =(g2float *)calloc(ndpts,sizeof(g2float));
         if( *fld == 0 )
         {
             free(ifld);
             return -1;
         }
         for (j=0;j<ndpts;j++) {
           (*fld)[j]=(((g2float)ifld[j]*bscale)+ref)*dscale;
         }
         free(ifld);
      }
      else {
         // Limit to 2 GB
         if( ndpts > 500 * 1024 * 1024 )
         {
             fprintf(stderr, "jpcunpack: ndpts = %d > 500 * 1024 * 1024", ndpts );
             return -1;
         }
         *fld =(g2float *)calloc(ndpts,sizeof(g2float));
         if( *fld == 0 )
         {
             return -1;
         }
         for (j=0;j<ndpts;j++) (*fld)[j]=ref * dscale;
      }

      return(0);
}