Line data    Source code

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