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


void specpack(g2float *fld,g2int ndpts,g2int JJ,g2int KK,g2int MM,
              g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
//$$$  SUBPROGRAM DOCUMENTATION BLOCK
//                .      .    .                                       .
// SUBPROGRAM:    specpack
//   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2002-12-19
//
// ABSTRACT: This subroutine packs a spectral data field using the complex
//   packing algorithm for spherical harmonic data as 
//   defined in the GRIB2 Data Representation Template 5.51.
//
// PROGRAM HISTORY LOG:
// 2002-12-19  Gilbert
//
// USAGE:    void specpack(g2float *fld,g2int ndpts,g2int JJ,g2int KK,g2int MM,
//                        g2int *idrstmpl,insigned char *cpack,g2int *lcpack)
//   INPUT ARGUMENT LIST:
//     fld[]    - Contains the packed data values
//     ndpts    - The number of data values to pack
//     JJ       - J - pentagonal resolution parameter
//     KK       - K - pentagonal resolution parameter
//     MM       - M - pentagonal resolution parameter
//     idrstmpl - Contains the array of values for Data Representation
//                Template 5.51
//
//   OUTPUT ARGUMENT LIST:
//     cpack    - The packed data field (character*1 array)
//     lcpack   - length of packed field cpack().
//
// REMARKS: None
//
// ATTRIBUTES:
//   LANGUAGE: C
//   MACHINE:  IBM SP
//
//$$$
{

      g2int    *ifld,tmplsim[5];
      g2float  bscale,dscale,*unpk,*tfld;
      g2float  *pscale,tscale;
      g2int    Js,Ks,Ms,Ts,Ns,inc,incu,incp,n,Nm,m,ipos;

      bscale = int_power(2.0,-idrstmpl[1]);
      dscale = int_power(10.0,idrstmpl[2]);
      Js=idrstmpl[5];
      Ks=idrstmpl[6];
      Ms=idrstmpl[7];
      Ts=idrstmpl[8];

//
//   Calculate Laplacian scaling factors for each possible wave number.
//
      pscale=(g2float *)malloc((JJ+MM)*sizeof(g2float));
      tscale=(g2float)idrstmpl[4]*1E-6;
      for (n=Js;n<=JJ+MM;n++)
           pscale[n]=pow((g2float)(n*(n+1)),tscale);
//
//   Separate spectral coeffs into two lists; one to contain unpacked
//   values within the sub-spectrum Js, Ks, Ms, and the other with values 
//   outside of the sub-spectrum to be packed.
//
      tfld=(g2float *)malloc(ndpts*sizeof(g2float));
      unpk=(g2float *)malloc(ndpts*sizeof(g2float));
      ifld=(g2int *)malloc(ndpts*sizeof(g2int));
      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) {       // save unpacked value
               unpk[incu++]=fld[inc++];         // real part
               unpk[incu++]=fld[inc++];     // imaginary part
            }
            else {                       // Save value to be packed and scale
                                         // Laplacian scale factor
               tfld[incp++]=fld[inc++]*pscale[n];      // real part
               tfld[incp++]=fld[inc++]*pscale[n];      // imaginary part
            }
         }
      }

      free(pscale);

      if (incu != Ts) {
         printf("specpack: Incorrect number of unpacked values %d given:\n",(int)Ts);
         printf("specpack: Resetting idrstmpl[8] to %d\n",(int)incu);
         Ts=incu;
      }
//
//  Add unpacked values to the packed data array in 32-bit IEEE format
//
      mkieee(unpk,(g2int *)cpack,Ts);
      ipos=4*Ts;
//
//  Scale and pack the rest of the coefficients
// 
      tmplsim[1]=idrstmpl[1];
      tmplsim[2]=idrstmpl[2];
      tmplsim[3]=idrstmpl[3];
      simpack(tfld,ndpts-Ts,tmplsim,cpack+ipos,lcpack);
      *lcpack=(*lcpack)+ipos;
//
//  Fill in Template 5.51
//
      idrstmpl[0]=tmplsim[0];
      idrstmpl[1]=tmplsim[1];
      idrstmpl[2]=tmplsim[2];
      idrstmpl[3]=tmplsim[3];
      idrstmpl[8]=Ts;
      idrstmpl[9]=1;         // Unpacked spectral data is 32-bit IEEE

      free(tfld);
      free(unpk);
      free(ifld);

      return;
}