#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h> /* for ldexp(3) */
#include <limits.h> /* for CHAR_BIT */

unsigned long
get2u(unsigned char *buf)
{
  return (buf[0] << 8) | buf[1];
}

long
get2i(unsigned char *buf)
{
  unsigned long u;
  u = ((0x7F & buf[0]) << 8) | buf[1];
  return (buf[0] & 0x80) ? -u : u;
}

unsigned long
get3u(unsigned char *buf)
{
  return (buf[0] << 16) | (buf[1] << 8) | buf[2];
}

long
get3i(unsigned char *buf)
{
  unsigned long u;
  u = ((0x7F & buf[0]) << 16) | (buf[1] << 8) | buf[2];
  return (buf[0] & 0x80) ? -u : u;
}

float
get4f(unsigned char *buf)
{
  double r;
  r = ldexp((float)get3u(buf + 1),
    ((0x7F & *buf) - 64) * 4 - 24);
  if (0x80 & *buf) {
    r = -r;
  }
  return r;
}

unsigned long
getbits(unsigned char *bits, unsigned width, unsigned pos)
{
  unsigned long r = 0;
  unsigned ihave, iused, ineed;
  ineed = width;
  bits += (width * pos) / CHAR_BIT;
  iused = (width * pos) % CHAR_BIT;
  ihave = CHAR_BIT - iused;
  while (ineed) {
    unsigned val, mask;
    mask = (1 << (CHAR_BIT - iused)) - 1;
    val = bits[0] & mask;
    if (ineed > ihave) {
      r <<= ihave;
      r |= val;
      ineed -= ihave;
      bits++;
      ihave = CHAR_BIT;
      iused = 0;
    } else {
      r |= (val >> (ihave - ineed));
      break;
    }
  }
  return r;
}

int
msgcheck(unsigned char *msgbuf, unsigned long msgl)
{
  unsigned char *sec1;
  unsigned long sec1len;
  unsigned char *sec2;
  unsigned long sec2len;
  unsigned nv, pl;
  unsigned  nx, ny, nxmax;
  long la1, lo1, la2, lo2;
  unsigned flags;
  unsigned long ncell, icell;

  unsigned char *sec4;
  unsigned char *bits;
  unsigned w;
  float r;
  int d, e;

  sec1 = msgbuf + 8;
  sec1len = get3u(sec1);
  flags = sec1[7];
  if ((flags & 0x80) == 0) {
    fprintf(stderr, "sec2 missing\n"); return -1;
  }
  if (flags & 0x40) {
    fprintf(stderr, "sec3 present\n"); return -1;
  }
  sec2 = sec1 + sec1len;
  nv = sec2[3];
  pl = sec2[4] + 4 * nv;
  if (sec2[5] != 0) {
    fprintf(stderr, "sec2[5]=%u: only 0 (lat-long) supported\n", sec2[5]);
    return -1;
  }
  nx = get2u(sec2 + 6);
  if (nx == 0xFFFF) nx = 0;
  nxmax = nx;
  ny = get2u(sec2 + 8);
  if (ny == 0xFFFF) {
    fprintf(stderr, "variable ny\n"); return -1;
  }
  la1 = get3i(sec2 + 10);
  lo1 = get3i(sec2 + 13);
  la2 = get3i(sec2 + 17);
  lo2 = get3i(sec2 + 20);
  if (nx != 0) {
    ncell = nx * ny;
  } else {
    unsigned char *nxlist = sec2 + pl - 1;
    unsigned inx, iy;
    ncell = 0;
    for (iy = 0; iy < ny; iy++) {
      inx = get2u(nxlist + 2 * iy);
      ncell += inx;
      if (inx > nxmax) nxmax = inx;
    }
  }
  printf("c%u p%u v%u/%04lu t%u-%02u-%u r%04u%02u%02uT%02u%02u",
    sec1[7], sec1[8], sec1[9], get2u(sec1 + 10),
    sec1[17], sec1[18], sec1[19],
    (sec1[24] - 1) * 100 + sec1[12],
    sec1[13], sec1[14], sec1[15], sec1[16]);
  printf(" %ux%u=%lu %.3g:%.3g %.3g:%.3g",
    nxmax, ny, ncell, la1 * 0.001, la2 * 0.001,
    lo1 * 0.001, lo2 * 0.001);
  putchar('\n');

  sec4 = sec2 + get3u(sec2);
  if ((sec4[3] & 0xD0) != 0) {
    fprintf(stderr, "sec2[3]=%02x: bits 0xD0 are to be cleared\n", sec4[3]);
    return -1;
  }
  d = get2i(sec1 + 26);
  e = get2i(sec4 + 4);
  r = get4f(sec4 + 6);
  w = sec4[10];
  printf(" d%+03d e%+03d r%-.3f w%u\n", d, e, r, w);
  bits = sec4 + 11;
  for (icell = 0; icell < ncell; icell++) {
    unsigned long y = getbits(bits, w, icell);
    printf(" %lu", y);
  }
  putchar('\n');

  return 0;
}

int
scanfile(FILE *gfp)
{
  unsigned long msgl;
  unsigned char *msgbuf;
  int state, c, r;
  r = 0;
  state = 0;
  while ((c = fgetc(gfp)) != EOF) {
    if (state == 0) {
      if (c == 'G') { state = 'G'; }
    } else if (state == 'G') {
      if (c == 'R') { state = 'R'; } else { state = 0; }
    } else if (state == 'R') {
      if (c == 'I') { state = 'I'; } else { state = 0; }
    } else if (state == 'I') {
      if (c == 'B') { state = 'B'; } else { state = 0; }
    } else if (state == 'B') {
      msgl = c;
      state = 'C';
    } else if (state == 'C') {
      msgl = ((msgl << 8) | c);
      state = 'D';
    } else if (state == 'D') {
      msgl = ((msgl << 8) | c);
      state = 0;
      msgbuf = malloc(msgl);
      if (msgbuf == NULL) {
        perror("malloc");
	return 1;
      }
      memcpy(msgbuf, "GRIB", 4);
      msgbuf[4] = (0xFF & (msgl >> 16));
      msgbuf[5] = (0xFF & (msgl >> 8));
      msgbuf[6] = (0xFF & msgl);
      if (fread(msgbuf + 7, 1, msgl - 7, gfp) < (msgl - 7)) {
        perror("fread");
	return 1;
      }
      if (memcmp(msgbuf + msgl - 4, "7777", 4) != 0) {
        fprintf(stderr, "GRIB 7777 check failed\n");
	r = 1;
      } else {
        msgcheck(msgbuf, msgl);
      }
      free(msgbuf);
    }
  }
  return r;
}

int
main(int argc, char **argv)
{
  char *infile;
  FILE *gfp;
  infile = getenv("GRIBIN");
  if (!infile) {
    infile = "GRIBIN";
  }
  gfp = fopen(infile, "rb");
  if (gfp == NULL) {
    perror(infile);
    return 1;
  }
  scanfile(gfp);
  fclose(gfp);
  return 0;
}