/* File idl_crunch.c */
/* A port of ANA's crunch.c to IDL linkimage format. */
/* Authors: Richard A. Shine (original) */
/*          Louis H. Strous (this port) */
/*          Lockheed-Martin Advanced Technology Center, Palo Alto, CA */

#include <stdio.h>
#include "/local/rsi/idl_5/external/export.h"
#include <malloc.h>
#undef putc
#undef putchar

#define IDL_TYP_INTEGER (IDL_TYP_MASK(IDL_TYP_BYTE) | IDL_TYP_MASK(IDL_TYP_INT) | IDL_TYP_MASK(IDL_TYP_LONG))
#define IDL_TYP_NUMERICAL (IDL_TYP_B_ALL & ~IDL_TYP_MASK(IDL_TYP_STRING) & ~IDL_TYP_MASK(IDL_TYP_STRUCT))

static union { short int i; unsigned char b[2]; } endian_temp = { 1 };
#define little_endian	(endian_temp.b[0])

typedef unsigned char	byte;
byte bits[8] = { 1, 2, 4, 8, 16, 32, 64, 128 };
int	ana_type_size[] = { 
  sizeof(char), sizeof(short int), sizeof(int), sizeof(float), sizeof(double)
};

struct fzHeadStruct {
  int	synch_pattern;
  byte	subf, source, nhb, datyp, ndim, free1, cbytes[4], free[178];
  int	dim[16];
  char	txt[256];
  char	*extratxt;
};
typedef struct fzHeadStruct	fzHead;

void	endian(void *, int, int);
#define swapl(p,n)	endian(p,n*sizeof(int),2)

IDL_VPTR ana_crunch(int argc, IDL_VPTR argv[])
/* call: result = CRUNCH(array, byte_slice, out [, crunch_run_flag]) */
/* <array>: numerical array to be compressed */
/* <byte_slice>: number of bits that get passed on without compression */
/* <out>: predefined output array */
/* <crunch_run_flag>: if non-zero (as an integer), then apply run-length */
/* encoding, too */
/* <result> will be a 2-element float array; element 0 contains */
/* the number of bits in the compressed data, and element 1 the number */
/* of bits per original element. */
{
  static IDL_EZ_ARG arg_struct[] = {
    { IDL_EZ_DIM_ARRAY, IDL_TYP_MASK(IDL_TYP_BYTE) | IDL_TYP_MASK(IDL_TYP_INT),
      IDL_EZ_ACCESS_R, IDL_TYP_UNDEF, 0, 0 }, /* <array>: integer array */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_NUMERICAL, IDL_EZ_ACCESS_R, IDL_TYP_INT,
      0, 0 },			/* <byte_slice>: long integer scalar */
    { IDL_EZ_DIM_ARRAY, IDL_TYP_NUMERICAL, IDL_EZ_ACCESS_RW, IDL_TYP_UNDEF,
      0, 0 }, /* <out> */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_NUMERICAL, IDL_EZ_ACCESS_R, IDL_TYP_LONG,
      0, 0 },			/* <crunch_run_flag>: long integer scalar */
  };
  void	*q1, *q2;
  IDL_VPTR	result;
  int	slice, nx, outer, i, limit, crunch_run_flag, type, ctype, iq;
  IDL_LONG	two = 2;
  float	*info;
  int anacrunch(byte *x, short array[], int slice, int nx, int ny, int limit),
    anacrunch8(byte *x, byte array[], int slice, int nx, int ny, int limit),
    anacrunchrun(byte *x, short array[], int slice, int nx, int ny, int limit),
    anacrunchrun8(byte *x, byte array[], int slice, int nx, int ny, int limit);

  if (argc < 3)
    IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP, "Too few arguments");

  IDL_EzCall(argc, argv, arg_struct); /* check arguments */

  /* <array> info */
  q1 = argv[0]->value.arr->data;/* pointer to its data */
  nx = argv[0]->value.arr->dim[0]; /* # elements in first dimension */
  outer = 1;			/* # elements in other dimensions: */
  for (i = 1; i < argv[0]->value.arr->n_dim; i++)
    outer *= argv[0]->value.arr->dim[i];
  type = argv[0]->type;

  /* <byte_slice> value */
  slice = (int) argv[1]->value.i;
  if (slice < 0) {
    IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_RET,
		"<slice> is negative: changing to 0");
    slice = 0;
  } else if (slice >= argv[0]->value.arr->elt_len*8) {
    IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_RET,
		"<slice> is too great: changing to 0");
    slice = 0;
  }
    
  /* <out> info */
  q2 = argv[2]->value.arr->data;
  limit = argv[2]->value.arr->arr_len; /* # bytes in <out> */

  /* <crunch_run_flag> value */
  crunch_run_flag = (argc > 3)? argv[3]->value.l: 0;
  
  switch (type) {
    case IDL_TYP_BYTE:
      ctype = 1;
      break;
    case IDL_TYP_INT:
      ctype = 0;
      break;
  }
  if (crunch_run_flag)
    ctype += 2;

  info = (float *) IDL_MakeTempArray(IDL_TYP_FLOAT, 1, &two,
				     IDL_BARR_INI_NOP, &result);
  
  switch (ctype) {
    case 0:
      iq = anacrunch(q2, q1, slice, nx, outer, limit);
      break;
    case 1:
      iq = anacrunch8(q2, q1, slice, nx, outer, limit);
      break;
    case 2:
      iq = anacrunchrun(q2, q1, slice, nx, outer, limit);
      break;
    case 3:
      iq = anacrunchrun8(q2, q1, slice, nx, outer, limit);
      break;
  }
  if (iq < 0) {
    IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_RET,
		"Not enough space allocated for compressed array");
    info[0] = info[1] = 0;
  } else {
    info[0] = 8*iq;
    info[1] = 8.0*iq/(nx*outer);
  }
  IDL_EzCallCleanup(argc, argv, arg_struct);
  return result;
}
/*--------------------------------------------------------------------*/
void ana_decrunch(int argc, IDL_VPTR argv[])
/* call: decrunch, in, out */
/* <int>: compressed data */
/* <out>: decompressed data */
{
  static IDL_EZ_ARG arg_struct[] = {
    { IDL_EZ_DIM_ARRAY, IDL_TYP_NUMERICAL, IDL_EZ_ACCESS_R, IDL_TYP_UNDEF,
      0, 0 }, /* <in>: numerical array */
    { IDL_EZ_DIM_ANY, IDL_TYP_B_ALL, IDL_EZ_ACCESS_W, IDL_TYP_UNDEF,
      0, 0 } /* <out>: anything */
  };
  union { byte *b; int *l; }	q1;
  int	bsize, outer, nx, slice, ctype, type, nd;
  IDL_VPTR	result;
  IDL_LONG	dim[2];
  void	*q2;

  IDL_EzCall(argc, argv, arg_struct); /* check arguments */

  q1.l = (int *) argv[0]->value.arr->data; /* pointer to compressed data */
  bsize = *q1.l++;
  outer = *q1.l++;
  nx = *q1.l++;
  slice = *q1.b++;
  ctype = *q1.b++;
  if (!little_endian) {
    swapl(&bsize,1);
    swapl(&nx,1);
    swapl(&outer,1); 
  }

  switch (ctype % 2) {
    case 0:
      type = IDL_TYP_INT;
      break;
    case 1:
      type = IDL_TYP_BYTE;
      break;
  }
  nd = outer? 2: 1;
  dim[0] = nx;
  dim[1] = outer;
  q2 = IDL_MakeTempArray(type, nd, dim, IDL_BARR_INI_NOP, &result);
  if (!q2)
    IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_RET,
		"Problem setting up result for DECRUNCH");
  else
    switch (ctype) {
      case 0:
	anadecrunch(q1.b, q2, slice, nx, outer);
	break;
      case 1:
	anadecrunch8(q1.b, q2, slice, nx, outer);
	break;
      case 2:
	anadecrunchrun(q1.b, q2, slice, nx, outer);
	break;
      case 3:
	anadecrunchrun8(q1.b, q2, slice, nx, outer);
	break;
    }

  IDL_VarCopy(result, argv[1]);
  IDL_EzCallCleanup(argc, argv, arg_struct);
}
/*--------------------------------------------------------------------*/
void endian(void *pp, int n, int type)
/* swap bytes according to data type, starting at p, spanning n bytes.
   goes from bigendian to little_endian or back.
   byte -> do nothing
   word -> swap 1 2 to 2 1
   long, float -> swap 1 2 3 4 to 4 3 2 1
   double -> swap 1 2 3 4 5 6 7 8 to 8 7 6 5 4 3 2 1
   Works between Ultrix and Irix.   LS 10/12/92	*/
{
 int	size, i, n2, n3;
 byte	temp, *p2, *p;

 p = (byte *) pp;
 size = ana_type_size[type];
 if (size <= 1)
   return;
 n /= size;			/* number of elements */
 n2 = size/2;
 n3 = size + n2;
 p2 = p + size;
 while (n--) {
   for (i = 0; i < n2; i++) {
     temp = *p;
     *p++ = *--p2;
     *p2 = temp;
   }
   p += n2;
   p2 += n3;
 }
}
/*--------------------------------------------------------------------*/
int anacrunch(byte *x, short array[], int slice, int nx, int ny, int limit)
/* compress 16 bit array into x (a byte array) using ny blocks each of size
   nx, bit slice size slice, returns # of bytes in x */
{
  struct compresshead {
    int     tsize,nblocks,bsize;
    byte    slice_size,type; } *ch;
  unsigned nb,ixa,ixb;
  unsigned register i,j,r1,in;
  int r0,r2,r3,mask,fac;
  int i2,k,iy;
  union { int i; short w; unsigned char b[4]; } y;

  /* begin execution */
  limit = limit - 10;	/* some margin since we don't check all times */
  mask=1; for (i=0;i<slice;i++) mask=2*mask;
  fac=mask;       mask=mask-1; /* no inline expon. in C */
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
  else { if (slice < 10) nb=2; else nb=3;    }};
  y.i=0;
  /* do the compression header */
  ch = (struct compresshead *) x;
  /* important note - can't use the sizeof(struct compresshead) because it
     is 14 on some machines and rounded up to 16 on others */
  /*x = x + sizeof(struct compresshead);*/
  x = x + 14;
  ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 0;
  i=0;    r1=0;   in=0;
  for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
    /* load the first value, reverse bytes (VAX style)*/
    if (little_endian) {
      y.w=array[in]   ;x[i]=y.b[0]    ;x[i+1]=y.b[1];
    } else {
      y.w=array[in]   ;x[i]=y.b[1]    ;x[i+1]=y.b[0];
    }
    r1=r1+16;
    r2=0;
    ixa=1+iy*nx;    ixb=(iy+1)*nx;
    for (in=ixa; in<ixb; in++)      { /* start of ix (inner) loop */
      /* first the fixed slice portion */
      y.i=array[in]-array[in-1];
      r3=(y.i>>slice);
      i=r1>>3;
      j=r1%8;
      if ( i > limit ) return -1;		/* bad news, went too far */
      /* now load nb bytes into x */
      /*low order byte of y.i is first in stream */
      if (j == 0) {
	y.i=(y.i & mask);	x[i]= (byte) y.i;
	/* since we started at bit 0, spillover to the next byte is
	   determined as follows (and is unlikely since slice gt 8 is unusual */
	if (slice > 8) x[i+1]= (byte) (y.i >> 8);
      } else {
	y.i=(y.i & mask)<<j;	x[i]=x[i] | (byte) y.i;
	/* spillover more likely here */
	if (nb>1) { x[i+1]= (byte) (y.i >> 8);
	if (nb>2) x[i+2]=(byte) (y.i >> 16); }
      }
      
      r1=r1+slice;       /* bump r1 pass the fixed part */
      i=r1>>3;                j=r1%8;
      /* note that r3 is the # of bits required minus 1 */
      if (r3==0) { if (j ==0 ) {x[i]= *bits;} else {x[i]=x[i]| bits[j];}
      r1++;} else    {
	r3=2*r3;        if (r3<0) r3 = -r3-1;
	if (r3<31)  {
	  r0=j+r3;		/* this is the bit that needs setting offset from x[i] */
	  if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
	  /* note, discovered on 2/3/96 that the j==0 case not done above for the
	     sunbow version, was OK in the umbra version, may have happened while cleaning
	     up code?, caused extra bits to be set if x[i] wasn't zero */
	  else  {
	    if (j == 0) x[i]=0;
	    j=r0%8;
	    if (r0 < 16) x[i+1]=bits[j]; else {
	      i2=i+r0/8;
	      for (k=i+1;k<i2;k++) x[k]=0;
	      x[i2]=bits[j]; }
	  }
	  r1+=1+r3;
	} else {        /* big one exception, should be rare */
	  /* does not need to be efficient, used rarely */
	  /*printf("big one \n");*/
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  r0=j+31;        j=r0%8;         i2=i+r0/8;
	  for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
	  /* recompute the difference and load 17 bits (always 3 bytes) */
	  r1=r1+32;
	  i=r1/8;
	  j=r1%8;
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  y.i=((array[in]-array[in-1])& 0x1ffff) << j;
	  if (little_endian) {
	    x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];      x[i+2]=y.b[2];
	  } else {
	    x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];      x[i+2]=y.b[1];
	  }
	  r1=r1+17;       } /* end of big one exception */
      }       /* end of (r3==0) conditional */
      /*in=in+1; */                   }       /* end of ix loop */
    /* some checks here */
    /* bump to next byte boundary */
    i=(r1+7)/8;     r1=8*i;                 }       /* end of iy loop */
  ch->tsize = i = i + 14;
  /* we have to put these in a form readable by the Vax (these may be used
     by fcwrite) */
  if (!little_endian) {
    swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
  }
  return  i;      /*return # of bytes used */
}       /* end of routine */
/*--------------------------------------------------------------------------*/
int anacrunch8(byte *x, byte array[], int slice, int nx, int ny, int limit)
/* compress 8 bit array into x (a byte array) using ny blocks each of size
   nx, bit slice size slice, returns # of bytes in x */
{
  struct compresshead {
    int     tsize,nblocks,bsize;
    byte    slice_size,type; } *ch;
  unsigned nb,ixa,ixb;
  unsigned register i,j,r1,in;
  int r0,r2,r3,mask,fac;
  int i2,k,iy;
  union { int i; short w; unsigned char b[4]; } y;

  /* begin execution */
  limit = limit - 10;	/* some margin since we don't check all times */
  mask=1; for (i=0;i<slice;i++) mask=2*mask;
  fac=mask;       mask=mask-1; /* no inline expon. in C */
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (slice > 8) slice = 8;
  if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
  else { if (slice < 10) nb=2; else nb=3;    }};
  y.i=0;
  /* do the compression header */
  ch = (struct compresshead *) x;
  /* important note - can't use the sizeof(struct compresshead) because it
     is 14 on some machines and rounded up to 16 on others */
  /*x = x + sizeof(struct compresshead);*/
  x = x + 14;
  ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 1;
  i=0;    r1=0;   in=0;
  for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
    /* load the first value */
    x[i] = array[in];
    r1=r1+8;
    r2=0;
    ixa=1+iy*nx;    ixb=(iy+1)*nx;
    for (in=ixa; in<ixb; in++)      {               /* start of ix (inner) loop */
      /* first the fixed slice portion */
      y.i= (int) array[in]- (int) array[in-1];
      r3=(y.i>>slice);
      i=r1>>3;
      j=r1%8;
      if ( i > limit ) return -1;		/* bad news, went too far */
      /* now load nb bytes into x */
      /*low order byte of y.i is first in stream */
      if (little_endian) {
	if (j == 0) {y.i=(y.i & mask); x[i]=y.b[0];}
	else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[0];}
	if (nb>1) { x[i+1]=y.b[1]; }
      } else {
	if (j == 0) {y.i=(y.i & mask); x[i]=y.b[3];}
	else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[3];}
	if (nb>1) { x[i+1]=y.b[2]; }
      }
      r1=r1+slice;       /* bump r1 pass the fixed part */
      i=r1>>3;                j=r1%8;
      /* note that r3 is the # of bits required minus 1 */
      if (r3==0) { if (j ==0 ) {x[i]=bits[j];} else {x[i]=x[i]|bits[j];}
      r1+=1;} else    {
	r3=2*r3;        if (r3<0) r3 = -r3-1;
	if (r3<31)  {
	  r0=j+r3;        /* this is the bit that needs setting offset from x[i] */
	  if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
	  else {if (j == 0) x[i]=0;  j=r0%8; if (r0 < 16) x[i+1]=bits[j];
	  else { i2=i+r0/8; for (k=i+1;k<i2;k++) x[k]=0;  x[i2]=bits[j]; }
	  }
	  r1+=1+r3;
	} else {        /* big one exception, should be rare */
	  /* does not need to be efficient, used rarely */
	  /* printf("big one \n"); */
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  r0=j+31;        j=r0%8;         i2=i+r0/8;
	  for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
	  /* recompute the difference and load 9 bits (always 2 bytes) */
	  r1=r1+32;
	  i=r1/8;
	  j=r1%8;
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  y.i=((array[in]-array[in-1])& 0x1ff) << j;
	  if (little_endian) {
	    x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];
	  } else {
	    x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];
	  }
	  r1=r1+9;       } /* end of big one exception */
      }       /* end of (r3==0) conditional */
    }       /* end of ix loop */
    /* some checks here */
    /* bump to next byte boundary */
    i=(r1+7)/8;     r1=8*i;                 }       /* end of iy loop */
  ch->tsize = i = i + 14;
  /* we have to put these in a form readable by the Vax (these may be used
     by fcwrite) */
  if (!little_endian) {
    swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
  }
  return  i;      /*return # of bytes used */
}       /* end of routine */
/*--------------------------------------------------------------------------*/
int anacrunchrun(byte *x, short array[], int slice, int nx, int ny, int limit)
/* compress 16 bit array into x (a byte array) using ny blocks each of size
   nx, bit slice size slice, returns # of bytes in x */
{
  struct compresshead {
    int     tsize,nblocks,bsize;
    byte    slice_size,type; } *ch;
  short	*p;
  unsigned nb;
  unsigned register i,j,r1;
  int	r0,r2,r3,mask,fac, nrun, lrun, ic;
  int	*dif, *d, nc, zq, yq, *dd;
  /* enum	state { RUN, LITERAL }; */
  int	i2,k,iy;
  union { int i; short w; unsigned char b[4]; } y;
  IDL_VPTR	tempmem;

  /* begin execution */
  limit = limit - 10;	/* some margin since we don't check all times */
  mask=1; for (i=0;i<slice;i++) mask=2*mask;
  fac=mask;       mask=mask-1; /* no inline expon. in C */
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
  else { if (slice < 10) nb=2; else nb=3;    }};
  y.i=0;
  /* do the compression header */
  ch = (struct compresshead *) x;
  /* important note - can't use the sizeof(struct compresshead) because it
     is 14 on some machines and rounded up to 16 on others */
  x = x + 14;
  ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 2;
  i=0;    r1=0;
  dif = (int *) malloc(nx*4);			/* line buffer */
  for (iy=0;iy<ny;iy++) {                 	/* start of iy (outer) loop */
    /* load the first value, reverse bytes (VAX style)*/
    if (little_endian) { 
      y.w=array[iy*nx]   ;x[i++]=y.b[0]    ;x[i++]=y.b[1];
    } else {
      y.w=array[iy*nx]   ;x[i++]=y.b[1]    ;x[i++]=y.b[0];
    }
    /* compute and store the first differences for this line */
    p = (array+nx*iy);	nc=nx-1;
    d=dif; yq=(int) *p++;	zq=(int) *p++;
    while (nc--) { *d++ = zq - yq; yq = zq; zq = (int) *p++; }
    r1=r1+16;
    r2=0;
    p = (array+nx*iy);	nc=nx-1;
    d=dif;
    ic = i++;			/* count position */
    r1=r1+8;		/* to cover first count */
    lrun = 0;		/* literal count) */
    while (1) {
      /* look for a run */
      /* printf("*d, *(d+1) = %d %d, nc = %d\n",*d, *(d+1), nc );*/
      y.i = *d++;
      if (nc > 1) {
	while ( y.i == *d ) {	/* at least a run of 2 */
	  dd = d+1;	nrun = 2;
	  while ( nc-- > 2 && y.i == *dd) {
	    /* printf("run!, y.i, *dd = %d %d, nc = %d\n", y.i, *dd, nc ); */
	    nrun++;  dd++; }
	  /* short runs are not worth it, make the legal limit 4 */
	  /* printf("nrun = %d, nc = %d\n", nrun,nc);*/
	  if ( nrun >= 4 ) {	/* code the run */
	    /* a previous literal ? */
	    if (lrun != 0) {
	      /* printf("previous was literal, ic, i = %d %d\n", ic,i);*/
	      x[ic] = lrun;	i = (r1+7)/8;	lrun = 0;
	      /* printf("now, i = %d\n",i );*/
	    } else i=ic;
	    while (nrun  > 128 )	{ /* a big one, multiple runs */
	      /* need only 2 bytes to represent run, runs can't be 17 bits */
	      if (nrun == 129)	/* beware the dreaded 129 */
	      { x[i++] = 0x82; nrun -= 127;} else { x[i++] = 0x81; nrun -= 128; }
	      if (little_endian) {
		x[i++]=y.b[0];	x[i++]=y.b[1];
	      } else {
		x[i++]=y.b[3];	x[i++]=y.b[2];
	      }
	    }
	    /* printf("encoding run, nrun = %d, i=%d, iy = %d\n",nrun,i,iy); */
	    if (little_endian) {
	      x[i++] = -(nrun-1);	x[i++]=y.b[0];	x[i++]=y.b[1];
	    } else {
	      x[i++] = -(nrun-1);	x[i++]=y.b[3];	x[i++]=y.b[2];
	    }
	    /* prepare for a literal and for next run check */
	    nc--;
	    if (nc <= 0) goto ended_on_run;
	    lrun = 0;	ic = i++;	r1 = 8*i;	d = dd;
	    y.i = *d++;
	    /* printf("after a run, new y.i = %d, new *d = %d\n", y.i, *d);*/
	  } else { nc = nc + nrun -1;	break; }
	}	/* not a run, do next literal, assume setup for literals */
      } else if (nc <= 0) break;
      nc--;
      /* increment the literal count */
      /* printf("literal, lrun = %d, nc = %d, ic,i = %d %d\n", lrun, nc, ic,i);*/
      if (++lrun > 127)	{		/* need a new literal run count */
	x[ic] = 127;
	/* printf("ic = %d, i,r1 = %d %d\n", ic,i,r1); */
	/* bump to next byte boundary */
	i = (r1+7)/8;	ic = i++;     r1 = 8*i;		lrun = 1;
	/* printf("next ic = %d\n", ic); */
      }
      /* first the fixed slice portion */
      /* printf("y.i = %d\n", y.i);*/
      r3=(y.i>>slice);
      i=r1>>3;	/* byte number */
      j=r1 & 7;		/* bit number */
      if ( i > limit ) return -1;			/* bad news, went too far */
      /* now load nb bytes into x */
      /*low order byte of y.i is first in stream */
      if (little_endian) {
	if (j == 0) {y.i=(y.i & mask); x[i]=y.b[0];}
	else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[0];}
	if (nb>1) { x[i+1]=y.b[1]; if (nb>2) x[i+2]=y.b[2]; }
      } else {
	if (j == 0) {y.i=(y.i & mask); x[i]=y.b[3];}
	else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[3];}
	if (nb>1) { x[i+1]=y.b[2]; if (nb>2) x[i+2]=y.b[1]; }
      }
      r1=r1+slice;		/* bump r1 pass the fixed part */
      i=r1>>3;                j=r1 & 7;
      /* note that r3 is the # of bits required minus 1 */
      /* printf("r3 = %d\n", r3);*/
      if (r3==0) { if (j ==0 ) {x[i]=bits[j];} else {x[i]=x[i]|bits[j];}
      r1+=1;} else    {
	r3=2*r3;        if (r3<0) r3 = -r3-1;
	if (r3<31)  {
	  r0=j+r3;		/* this is the bit that needs setting offset from x[i] */
	  if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
	  else {if (j == 0) x[i]=0;  j=r0%8; if (r0 < 16) x[i+1]=bits[j];
	  else { i2=i+r0/8; for (k=i+1;k<i2;k++) x[k]=0;  x[i2]=bits[j]; }
	  }
	  r1+=1+r3;
	} else {        /* big one exception, should be rare */
	  /* does not need to be efficient, used rarely */
	  /* printf("big one \n");*/
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  r0=j+31;        j=r0%8;         i2=i+r0/8;
	  for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
	  /* recompute the difference and load 17 bits (always 3 bytes) */
	  r1=r1+32;
	  i=r1/8;
	  j=r1%8;
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  y.i=((*(d-1))& 0x1ffff) << j;
	  if (little_endian) {
	    x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];      x[i+2]=y.b[2];
	  } else {
	    x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];      x[i+2]=y.b[1];
	  }
	  r1=r1+17;       } /* end of big one exception */
      }       /* end of (r3==0) conditional */
      /* printf("end of literal, i,r1 = %d %d\n", i,r1);*/
      /* printf(" x[r1/8] = %d\n",  x[r1/8]);*/
    }       /* end of ix loop */
    /* some checks here */
    /* bump to next byte boundary */
    /* a final literal ? */
    /* printf("final lrun = %d, ic = %d\n", lrun, ic);*/
    if (lrun != 0) { x[ic] = lrun;	lrun = 0; }
    i=(r1+7)/8;
    ended_on_run:
    r1=8*i;
  }       /* end of iy loop */
  ch->tsize = i = i + 14;
  /* we have to put these in a form readable by the Vax (these may be used
     by fcwrite) */
  if (!little_endian) {
    swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
  }
  free(dif);
  return  i;      /*return # of bytes used */
}       /* end of routine */
/*--------------------------------------------------------------------------*/
int anacrunchrun8(byte *x, byte array[], int slice, int nx, int ny, int limit)
/* compress 8 bit array into x (a byte array) using ny blocks each of size
   nx, bit slice size slice, returns # of bytes in x */
{
  struct compresshead {
    int     tsize,nblocks,bsize;
    byte    slice_size,type; } *ch;
  byte	*p;
  unsigned nb;
  unsigned register i,j,r1;
  int	r0,r2,r3,mask,fac, nrun, lrun, ic;
  int	*dif, *d, nc, zq, yq, *dd;
  int	i2,k,iy;
  union { int i; short w; unsigned char b[4]; } y;

  /* begin execution */
  limit = limit - 10;	/* some margin since we don't check all times */
  mask=1; for (i=0;i<slice;i++) mask=2*mask;
  fac=mask;       mask=mask-1; /* no inline expon. in C */
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (slice == 0) nb=0; else { if (slice < 2 ) nb=1;
  else { if (slice < 10) nb=2; else nb=3;    }};
  y.i=0;
  /* do the compression header */
  ch = (struct compresshead *) x;
  /* important note - can't use the sizeof(struct compresshead) because it
     is 14 on some machines and rounded up to 16 on others */
  x = x + 14;
  ch->bsize = nx;  ch->nblocks = ny;  ch->slice_size = slice;  ch->type = 3;
  i=0;    r1=0;
  dif = (int *) malloc(nx*4);			/* line buffer */
  for (iy=0;iy<ny;iy++) {                 	/* start of iy (outer) loop */
    /* load the first value */
    x[i++] = array[iy*nx];
    
    /* compute and store the first differences for this line */
    p = (array+nx*iy);	nc=nx-1;
    d=dif; yq=(int) *p++;	zq=(int) *p++;
    while (nc--) { *d++ = zq - yq; yq = zq; zq = (int) *p++; }
    r1=r1+8;
    r2=0;
    p = (array+nx*iy);	nc=nx-1;
    d=dif;
    ic = i++;			/* count position */
    r1=r1+8;		/* to cover first count */
    lrun = 0;		/* literal count) */
    while (1) {
      /* look for a run */
      /* printf("*d, *(d+1) = %d %d, nc = %d\n",*d, *(d+1), nc );*/
      y.i = *d++;
      if (nc > 1) {
	while ( y.i == *d ) {	/* at least a run of 2 */
	  dd = d+1;	nrun = 2;
	  while ( nc-- > 2 && y.i == *dd) {
	    /* printf("run!, y.i, *dd = %d %d, nc = %d\n", y.i, *dd, nc ); */
	    nrun++;  dd++; }
	  /* short runs are not worth it, make the legal limit 4 */
	  /* printf("nrun = %d, nc = %d\n", nrun,nc);*/
	  if ( nrun >= 4 ) {	/* code the run */
	    /* a previous literal ? */
	    if (lrun != 0) {
	      /* printf("previous was literal, ic, i = %d %d\n", ic,i);*/
	      x[ic] = lrun;	i = (r1+7)/8;	lrun = 0;
	      /* printf("now, i = %d\n",i );*/
	    } else i=ic;
	    while (nrun  > 128 )	{ /* a big one, multiple runs */
	      /* printf("big run, nrun = %d\n", nrun); */
	      /* need only 2 bytes to represent run, runs can't be 17 bits */
	      if (nrun == 129)	/* beware the dreaded 129 */
	      { x[i++] = 0x82; nrun -= 127;} else { x[i++] = 0x81;	nrun -= 128; }
	      if (little_endian) {
		x[i++]=y.b[0];	x[i++]=y.b[1];
	      } else {
		x[i++]=y.b[3];	x[i++]=y.b[2];
	      }
	    }
	    /* printf("encoding run, nrun = %d, i=%d, iy = %d\n",nrun,i,iy); */
	    if (little_endian) {
	      x[i++] = -(nrun-1);	x[i++]=y.b[0];	x[i++]=y.b[1];
	    } else {
	      x[i++] = -(nrun-1);	x[i++]=y.b[3];	x[i++]=y.b[2];
	    }
	    /* prepare for a literal and for next run check */
	    nc--;
	    if (nc <= 0) goto ended_on_run;
	    lrun = 0;	ic = i++;	r1 = 8*i;	d = dd;
	    y.i = *d++;
	    /* printf("after a run, new y.i = %d, new *d = %d\n", y.i, *d);*/
	  } else { nc = nc + nrun -1;	break; }
	}	/* not a run, do next literal, assume setup for literals */
      } else if (nc <= 0) break;
      nc--;
      /* increment the literal count */
      /* printf("literal, lrun = %d, nc = %d, ic,i = %d %d\n", lrun, nc, ic,i);*/
      if (++lrun > 127)	{		/* need a new literal run count */
	x[ic] = 127;
	/* printf("ic = %d, i,r1 = %d %d\n", ic,i,r1); */
	/* bump to next byte boundary */
	i = (r1+7)/8;	ic = i++;     r1 = 8*i;		lrun = 1;
	/* printf("next ic = %d\n", ic); */
      }
      /* first the fixed slice portion */
      /* printf("y.i = %d\n", y.i);*/
      r3=(y.i>>slice);
      i=r1>>3;	/* byte number */
      j=r1 & 7;		/* bit number */
      if ( i > limit ) return -1;			/* bad news, went too far */
      /* now load nb bytes into x */
      /*low order byte of y.i is first in stream */
      if (little_endian) {
	if (j == 0) {y.i=(y.i & mask); x[i]=y.b[0];}
	else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[0];}
	if (nb>1) { x[i+1]=y.b[1]; if (nb>2) x[i+2]=y.b[2]; }
      } else {
	if (j == 0) {y.i=(y.i & mask); x[i]=y.b[3];}
	else { y.i=(y.i & mask)<<j; x[i]=x[i] | y.b[3];}
	if (nb>1) { x[i+1]=y.b[2]; if (nb>2) x[i+2]=y.b[1]; }
      }
      
      r1=r1+slice;		/* bump r1 pass the fixed part */
      i=r1>>3;                j=r1 & 7;
      /* note that r3 is the # of bits required minus 1 */
      if (r3==0) { if (j ==0 ) {x[i]=bits[j];} else {x[i]=x[i]|bits[j];}
      r1+=1;} else    {
	r3=2*r3;        if (r3<0) r3 = -r3-1;
	if (r3<31)  {
	  r0=j+r3;        /* this is the bit that needs setting offset from x[i] */
	  if (r0 < 8) { if (j == 0) x[i]=bits[r0]; else x[i]=x[i]|bits[r0];}
	  else {if (j == 0) x[i]=0;  j=r0%8; if (r0 < 16) x[i+1]=bits[j];
	  else { i2=i+r0/8; for (k=i+1;k<i2;k++) x[k]=0;  x[i2]=bits[j]; }
	  }
	  r1+=1+r3;
	} else {        /* big one exception, should be rare */
	  /* does not need to be efficient, used rarely */
	  /* printf("big one \n");*/
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  r0=j+31;        j=r0%8;         i2=i+r0/8;
	  for (k=i+1;k<i2;k++) x[k]=0;    x[i2]=bits[j];
	  /* recompute the difference and load 9 bits (always 2 bytes) */
	  r1=r1+32;
	  i=r1/8;
	  j=r1%8;
	  if (j == 0) x[i]=0;     /* gotta zero the virgins */
	  y.i=((*(d-1))& 0x1ff) << j;
	  if (little_endian) {
	    x[i]=x[i] | y.b[0]; x[i+1]=y.b[1];
	  } else {
	    x[i]=x[i] | y.b[3]; x[i+1]=y.b[2];
	  }
	  r1=r1+9;       } /* end of big one exception */
      }       /* end of (r3==0) conditional */
      /* printf("end of literal, i,r1 = %d %d\n", i,r1);*/
      /* printf(" x[r1/8] = %d\n",  x[r1/8]);*/
    }       /* end of ix loop */
    /* some checks here */
    /* bump to next byte boundary */
    /* a final literal ? */
    /* printf("final lrun = %d, ic = %d\n", lrun, ic);*/
    if (lrun != 0) { x[ic] = lrun;	lrun = 0; }
    i=(r1+7)/8;
    ended_on_run:
    r1=8*i;
  }       /* end of iy loop */
  ch->tsize = i = i + 14;
  /* we have to put these in a form readable by the Vax (these may be used
     by fcwrite) */
  if (!little_endian) {
    swapl(&(ch->tsize),1); swapl(&(ch->bsize),1); swapl(&(ch->nblocks),1); 
  }
  free(dif);
  return  i;      /*return # of bytes used */
}
/*--------------------------------------------------------------------------*/
int anadecrunch(byte *x, short array[], int r9, int nx, int ny)
/* decompress a bit stream in x; result is n I*2 elements, put in array;
   bit slice size r9 */
{
  short iq;
  int r0,r1,r2,r4,nb,mask;
  int j,in,i,k,ix,iy;
  unsigned char xq;
  union { int i; short w; unsigned char b[4]; } y;

  /* begin execution */
  mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
  /*printf("slice width = %d\n",r9);*/
  /*printf ("mask = %x, %d\n",mask,mask);*/
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
  else { if (r9 < 10) nb=2; else nb=3;    }};
  y.i=0;
  i=0;    r1=0;   in=0;
  for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
    /* get the first value */
    if (little_endian) {
      y.b[0]=x[i];  y.b[1]=x[i+1];    iq=y.w; array[in++]=iq;
    } else {
      y.b[0]=x[i+1];  y.b[1]=x[i];    iq=y.w; array[in++]=iq;
    }
    /*printf("first value = %d 0x%x\n",iq, iq);*/
    r1=r1+16;
    r2=0;
    for (ix=1; ix<nx; ix++) {
      /* first the fixed slice portion */
      i=r1/8;         j=r1%8;
      if (little_endian) {
	y.b[0]=x[i];
      } else {
	y.b[3]=x[i];
      }
      /* test effect on timing */
      if (little_endian) {
	if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
      } else {
	if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
      }
      /* shift and mask out the bit slice */
      r2=(y.i>>j) & mask;
      /*printf("r2 = %x, %d\n",r2,r2);*/
      /* the variable bit portion, find the first set bit */
      r1=r1+r9;       /* bump r1 pass the fixed part */
      i=r1/8;         j=r1%8;
      if ((xq=x[i]>>j) != 0) {
	/* caught it on first byte, find the bit */
	if ((xq&1) != 0) r0=1; else {
	  if ((xq&2) != 0) r0=2; else {
	    if ((xq&4) != 0) r0=3; else {
	      if ((xq&8) != 0) r0=4; else {
		if ((xq&16) != 0) r0=5; else {
		  if ((xq&32) != 0) r0=6; else {
		    if ((xq&64) != 0) r0=7; else {
		      if ((xq&128) != 0) r0=8; }}}}}}}}       else {
			/* not in first byte (or part of one) checked, carry on, first count bits in
			   that first byte */
			r0=8-j;
			/* check up to 4 more bytes, if not found than an error */
			for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
			  /* caught it here, find the bit and then jump from loop */
			  if ((xq&1) != 0) r0+=1; else {
			    if ((xq&2) != 0) r0+=2; else {
			      if ((xq&4) != 0) r0+=3; else {
				if ((xq&8) != 0) r0+=4; else {
				  if ((xq&16) != 0) r0+=5; else {
				    if ((xq&32) != 0) r0+=6; else {
				      if ((xq&64) != 0) r0+=7; else {
					if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
					/* add 8 bits for each all zero byte */
					if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue\n");
					printf("i = %d, r1 = %d, ix= %d, iy = %d\n",i,r1,ix,iy);
					return -1; }       }       }       }
      r1=r1+r0;       /* update pointer */
      /* r0 even or odd determines sign of difference */
      if ((r0&1) != 0) { 
	/* positive case */
	r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
      } else
      { if (r0 == 32) { 
	/* a long one, yank out the next 17 bits and use as difference */
	i=r1/8;         j=r1%8;
	if (little_endian) {
	  y.b[0]=x[i];
	  y.b[1]=x[i+1]; y.b[2]=x[i+2];
	} else {
	  y.b[3]=x[i];
	  y.b[2]=x[i+1]; y.b[1]=x[i+2];
	}
	/* shift and mask out the 17 bit slice */
	r2=(y.i>>j) & 0x1ffff;
	r1=r1+17;
	/* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
	if ( (r2& 0x10000) != 0 ) r2=r2 | 0xffff0000;
	r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
      } else {
	/* minus case (normal) */
	r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
      }}
      in=in+1; }	/* end of ix loop */
    i=(r1+7)/8;     r1=8*i; }	/* end of iy loop */
  return 1;
}  						     /* end of routine */
/*--------------------------------------------------------------------------*/
int anadecrunch8(byte *x, byte array[], int r9, int nx, int ny)
/* decompress a bit stream in x; result is n I*1 elements, put in array;
   bit slice size r9 */
/* byte version, modified from I*2 version
   r. shine 6/5/91 */
{
  byte iq;
  int r0,r1,r2,r4,nb,mask;
  int j,in,i,k,ix,iy;
  byte xq;
  union { int i; short w; unsigned char b[4]; } y;

  /* begin execution */
  mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
  else { if (r9 < 10) nb=2; else nb=3;    }};
  y.i=0;
  i=0;    r1=0;   in=0;
  for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
    /* get the first value */
    iq=x[i];        array[in++]=iq;
    r1=r1+8;
    r2=0;
    for (ix=1; ix<nx; ix++) {
      /* first the fixed slice portion */
      i=r1/8;         j=r1%8;
      if (little_endian) {
	y.b[0]=x[i];
	if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
      } else {
	y.b[3]=x[i];
	if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
      }
      /* shift and mask out the bit slice */
      r2=(y.i>>j) & mask;
      /* the variable bit portion, find the first set bit */
      r1=r1+r9;			/* bump r1 pass the fixed part */
      i=r1/8;         j=r1%8;
      if ((xq=x[i]>>j) != 0) {
	/* caught it on first byte, find the bit */
	if ((xq&1) != 0) r0=1; else {
	  if ((xq&2) != 0) r0=2; else {
	    if ((xq&4) != 0) r0=3; else {
	      if ((xq&8) != 0) r0=4; else {
		if ((xq&16) != 0) r0=5; else {
		  if ((xq&32) != 0) r0=6; else {
		    if ((xq&64) != 0) r0=7; else {
		      if ((xq&128) != 0) r0=8; }}}}}}}}       else {
			/* not in first byte (or part of one) checked, carry on, first count bits in
			   that first byte */
			r0=8-j;
			/* check up to 4 more bytes, if not found than an error */
			for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
			  /* caught it here, find the bit and then jump from loop */
			  if ((xq&1) != 0) r0+=1; else {
			    if ((xq&2) != 0) r0+=2; else {
			      if ((xq&4) != 0) r0+=3; else {
				if ((xq&8) != 0) r0+=4; else {
				  if ((xq&16) != 0) r0+=5; else {
				    if ((xq&32) != 0) r0+=6; else {
				      if ((xq&64) != 0) r0+=7; else {
					if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
					/* add 8 bits for each all zero byte */
					if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue");
					return -1; }       }       }       }
      r1=r1+r0;       /* update pointer */
      /* r0 even or odd determines sign of difference */
      if ((r0&1) != 0) { 
	/* positive case */
	r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
      } else
      { if (r0 == 32) { 
	/* a long one, yank out the next 9 bits and use as difference */
	i=r1/8;         j=r1%8;
	if (little_endian) {
	  y.b[0]=x[i];
	  y.b[1]=x[i+1];
	} else {
	  y.b[3]=x[i];
	  y.b[2]=x[i+1];
	}
	/* shift and mask out the 9 bit slice */
	r2=(y.i>>j) & 0x1ff;
	r1=r1+9;
	/* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
	if ( (r2& 0x100) != 0 ) r2=r2 | 0xffffff00;
	r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
      } else {
	/* minus case (normal) */
	r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
      } }
      in=in+1;                                }       /* end of ix loop */
    i=(r1+7)/8;     r1=8*i;                 }       /* end of iy loop */
  return 1;
}       /* end of routine */
/*--------------------------------------------------------------------------*/
int anadecrunchrun(byte *x, short array[], int r9, int nx, int ny)
/* decompress a bit stream in x; result is n I*2 elements, put in array;
   bit slice size r9 */
/* this version handles the run length encoding used in anacrunchrun */
{
  short iq;
  int r0,r1,r2,r4,nb,mask,nrun,n,nc;
  int j,in,i,k,ix,iy;
  unsigned char xq;
  union { int i; short w; unsigned char b[4]; } y;

  /* begin execution */
  mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
  /*printf("slice width = %d\n",r9);*/
  /*printf ("mask = %x, %d\n",mask,mask);*/
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
  else { if (r9 < 10) nb=2; else nb=3;    }};
  y.i=0;
  i=0;    r1=0;   in=0;
  for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
    /* get the first value, assume bytes reversed */
    if (little_endian) {
      y.b[0]=x[i++];	y.b[1]=x[i++];    iq=y.w; array[in++]=iq;
    } else {
      y.b[1]=x[i++];	y.b[0]=x[i++];    iq=y.w; array[in++]=iq;
    }
    /* printf("first value = %d 0x%x\n",iq, iq); */
    r1=r1+16;
    r2=0;
    nc=nx-1;
    while (nc>0) {
      /* look at the next run length code */
      /* printf("i = %d\n", i); */
      nrun = (int) x[i++];
      /* printf("nrun = %d\n", nrun); */
      if (nrun > 127) {	/* a run of a constant difference */
	n = 255 - nrun + 2;
	nc = nc - n;
	if (little_endian) {
	  y.b[0]=x[i++];	y.b[1]=x[i++];
	} else {
	  y.b[1]=x[i++];	y.b[0]=x[i++];
	}
	/* printf("increment (run) = %d\n", y.w); */
	while (n--) {
	  array[in] = array[in-1] + y.w;	in++;
	}
	iq = array[in-1];	r1=8*i;
      } else {	/* a literal */
	r1 = 8 * i;
	nc = nc - nrun;
	while(nrun--) {
	  /* first the fixed slice portion */
	  i=r1/8;         j=r1%8;
	  /* printf("start literal, i,r1 = %d %d\n", i,r1); */
	  if (little_endian) {
	    y.b[0]=x[i];
	    /* test effect on timing */
	    if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
	  } else {
	    y.b[3]=x[i];
	    /* test effect on timing */
	    if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
	  }
	  /* shift and mask out the bit slice */
	  r2=(y.i>>j) & mask;
	  /* printf("r2 = %x, %d\n",r2,r2);*/
	  /* the variable bit portion, find the first set bit */
	  r1=r1+r9;       /* bump r1 pass the fixed part */
	  i=r1/8;         j=r1%8;
	  if ((xq=x[i]>>j) != 0) {
	    /* caught it on first byte, find the bit */
	    if ((xq&1) != 0) r0=1; else {
	      if ((xq&2) != 0) r0=2; else {
		if ((xq&4) != 0) r0=3; else {
		  if ((xq&8) != 0) r0=4; else {
		    if ((xq&16) != 0) r0=5; else {
		      if ((xq&32) != 0) r0=6; else {
			if ((xq&64) != 0) r0=7; else {
			  if ((xq&128) != 0) r0=8; }}}}}}}}       else {
			    /* not in first byte (or part of one) checked, carry on, first count bits in
			       that first byte */
			    r0=8-j;
			    /* check up to 4 more bytes, if not found than an error */
			    for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
			      /* caught it here, find the bit and then jump from loop */
			      if ((xq&1) != 0) r0+=1; else {
				if ((xq&2) != 0) r0+=2; else {
				  if ((xq&4) != 0) r0+=3; else {
				    if ((xq&8) != 0) r0+=4; else {
				      if ((xq&16) != 0) r0+=5; else {
					if ((xq&32) != 0) r0+=6; else {
					  if ((xq&64) != 0) r0+=7; else {
					    if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
					    /* add 8 bits for each all zero byte */
					    if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue\n");
					    printf("i = %d, r1 = %d, ix= %d, iy = %d\n",i,r1,ix,iy);
					    return -1; }       }       }       }
	  r1=r1+r0;       /* update pointer */
	  /* r0 even or odd determines sign of difference */
	  if ((r0&1) != 0) { 
	    /* positive case */
	    r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
	    /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq);*/
	  } else
	  { if (r0 == 32) { 
	    /* a long one, yank out the next 17 bits and use as difference */
	    i=r1/8;         j=r1%8;
	    if (little_endian) {
	      y.b[0]=x[i];
	      y.b[1]=x[i+1]; y.b[2]=x[i+2];
	    } else {
	      y.b[3]=x[i];
	      y.b[2]=x[i+1]; y.b[1]=x[i+2];
	    }
	    /* shift and mask out the 17 bit slice */
	    r2=(y.i>>j) & 0x1ffff;
	    r1=r1+17;
	    /* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
	    if ( (r2& 0x10000) != 0 ) r2=r2 | 0xffff0000;
	    /* printf("big one, r2 = %d, array[in-1] = %d\n", r2, array[in-1]);*/
	    r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
	  } else {
	    /* minus case (normal) */
	    r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
	    /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq); */
	  }}
	  /* printf("end literal, i,r1 = %d %d\n", i,r1); */
	  in=in+1;
	}
	/* printf("r1, i after nrun exhausted = %d %d\n",r1,i); */
	i=(r1+7)/8;     r1=8*i;
	/* printf("new i = %d\n",i); */
      }
    }   	    /* end of ix loop */
    if (nc < 0) {
      printf("bad loop in decrunchrun, nc=%d, iy=%d, in= %d\n",nc,iy,in);  return -1; }
    
    i=(r1+7)/8;     r1=8*i;                 }   	    /* end of iy loop */
  return 1;
}  						     /* end of routine */
/*--------------------------------------------------------------------------*/
int anadecrunchrun8(byte x[], byte array[], int r9, int nx, int ny)
/* decompress a bit stream in x; result is n I*2 elements, put in array;
   bit slice size r9 */
/* this version handles the run length encoding used in anacrunchrun */
{
  byte iq;
  int r0,r1,r2,r4,nb,mask,nrun,n,nc;
  int j,in,i,k,ix,iy;
  unsigned char xq;
  union { int i; short w; unsigned char b[4]; } y;

  /* begin execution */
  mask=1; for (i=0;i<r9;i++) mask=2*mask; mask=mask-1;
  /* determine the # of bytes to transfer to 32 bit int for fixed portion */
  if (r9 == 0) nb=0; else { if (r9 < 2 ) nb=1;
  else { if (r9 < 10) nb=2; else nb=3;    }};
  y.i=0;
  i=0;    r1=0;   in=0;
  for (iy=0;iy<ny;iy++) {                 /* start of iy (outer) loop */
    /* get the first value */
    iq=x[i++];	array[in++]=iq;
    /* printf("first value = %d 0x%x\n",iq, iq); */
    r1=r1+16;
    r2=0;
    nc=nx-1;
    /* printf("nc = %d\n", nc); */
    while (nc>0) {
      /* look at the next run length code */
      /* printf("i = %d\n", i); */
      nrun = (int) x[i++];
      /* printf("nrun = %d\n", nrun); */
      if (nrun > 127) {	/* a run of a constant difference */
	n = 255 - nrun + 2;
	nc = nc - n;
	if (little_endian) {
	  y.b[0]=x[i++];	y.b[1]=x[i++];
	} else {
	  y.b[1]=x[i++];	y.b[0]=x[i++];
	}
	/* printf("increment (run) = %d of length %d\n", y.w,n); */
	while (n--) {
	  array[in] = array[in-1] + y.w;	in++;
	}
	iq = array[in-1];	r1=8*i;
      } else {	/* a literal */
	r1 = 8 * i;
	nc = nc - nrun;
	while(nrun--) {
	  /* first the fixed slice portion */
	  i=r1/8;         j=r1%8;
	  /* printf("start literal, i,r1 = %d %d\n", i,r1); */
	  if (little_endian) {
	    y.b[0]=x[i];
	    /* test effect on timing */
	    if (nb>1) { y.b[1]=x[i+1]; if (nb>2) y.b[2]=x[i+2]; }
	  } else {
	    y.b[3]=x[i];
	    /* test effect on timing */
	    if (nb>1) { y.b[2]=x[i+1]; if (nb>2) y.b[1]=x[i+2]; }
	  }
	  /* shift and mask out the bit slice */
	  r2=(y.i>>j) & mask;
	  /* the variable bit portion, find the first set bit */
	  r1=r1+r9;       /* bump r1 pass the fixed part */
	  i=r1/8;         j=r1%8;
	  if ((xq=x[i]>>j) != 0) {
	    /* caught it on first byte, find the bit */
	    if ((xq&1) != 0) r0=1; else {
	      if ((xq&2) != 0) r0=2; else {
		if ((xq&4) != 0) r0=3; else {
		  if ((xq&8) != 0) r0=4; else {
		    if ((xq&16) != 0) r0=5; else {
		      if ((xq&32) != 0) r0=6; else {
			if ((xq&64) != 0) r0=7; else {
			  if ((xq&128) != 0) r0=8; }}}}}}}}       else {
			    /* not in first byte (or part of one) checked, carry on, first count bits in
			       that first byte */
			    r0=8-j;
			    /* check up to 4 more bytes, if not found than an error */
			    for (k=i+1;k<i+5;k++) { if ( (xq=x[k]) != 0 ) {
			      /* caught it here, find the bit and then jump from loop */
			      if ((xq&1) != 0) r0+=1; else {
				if ((xq&2) != 0) r0+=2; else {
				  if ((xq&4) != 0) r0+=3; else {
				    if ((xq&8) != 0) r0+=4; else {
				      if ((xq&16) != 0) r0+=5; else {
					if ((xq&32) != 0) r0+=6; else {
					  if ((xq&64) != 0) r0+=7; else {
					    if ((xq&128) != 0) r0+=8; }}}}}}} break; } else { r0=r0+8; 
					    /* add 8 bits for each all zero byte */
					    if (r0 > 32) { printf("DECRUNCH -- bad bit sequence, cannot continue\n");
					    printf("i = %d, r1 = %d, ix= %d, iy = %d\n",i,r1,ix,iy);
					    return -1; }       }       }       }
	  r1=r1+r0;       /* update pointer */
	  /* r0 even or odd determines sign of difference */
	  if ((r0&1) != 0) { 
	    /* positive case */
	    r0=(r0/2)<<r9;  iq=iq+r2;       iq=iq+r0;       array[in]=iq;
	    /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq);*/
	  } else
	  { if (r0 == 32) { 
	    /* a long one, yank out the next 9 bits and use as difference */
	    i=r1/8;         j=r1%8;
	    if (little_endian) {
	      y.b[0]=x[i];
	      y.b[1]=x[i+1];
	    } else {
	      y.b[3]=x[i];
	      y.b[2]=x[i+1];
	    }
	    /* shift and mask out the 9 bit slice */
	    r2=(y.i>>j) & 0x1ff;
	    r1=r1+9;
	    /* if the top bit was set, do a sign extend, note that 32 bit arithmetic used*/
	    if ( (r2& 0x100) != 0 ) r2=r2 | 0xffffff00;
	    /* printf("long one decoded, r2 = %d, array[in-1]=%d\n", r2, array[in-1]); */
	    r4=array[in-1]; r4=r4+r2; array[in]=r4; iq=r4;
	  } else {
	    /* minus case (normal) */
	    r0=(-r0/2)<<r9; iq=iq+r2;       iq=iq+r0;       array[in]=iq;
	    /* printf("r0,r2,iq = %d %d %d\n", r0,r2,iq); */
	  }}
	  /* printf("end literal, i,r1 = %d %d\n", i,r1); */
	  in=in+1;
	}
	/* printf("r1, i after literal nrun exhausted = %d %d\n",r1,i); */
	i=(r1+7)/8;     r1=8*i;
	/* printf("new i = %d\n",i); */
      }
    }   	    /* end of ix loop */
    if (nc < 0) {
      printf("bad loop in decrunchrun8, nc=%d, iy=%d, in= %d\n",nc,iy,in);
      return -1; }
    
    i=(r1+7)/8; r1=8*i; }	/* end of iy loop */
  return 1;
} /* end of routine */
/*--------------------------------------------------------------------------*/
void ana_fzread(int argc, IDL_VPTR argv[])
/* call: FZREAD, target, filename [, header] */
/* <target>: target variable */
/* <filename>: name of the FZ file (compressed or uncompressed) */
/* <header>: variable in which the header string will be returned */
{
  static IDL_EZ_ARG arg_struct[] = {
    { IDL_EZ_DIM_ANY, IDL_TYP_B_ALL, IDL_EZ_ACCESS_W, IDL_TYP_UNDEF,
      0, 0 },			/* <target>: anything */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_MASK(IDL_TYP_STRING), IDL_EZ_ACCESS_R,
      IDL_TYP_UNDEF, 0, 0 },	/* <filename>: name of the file */
    { IDL_EZ_DIM_ANY, IDL_TYP_B_ALL, IDL_EZ_ACCESS_W, IDL_TYP_UNDEF,
      0, 0 }			/* <header>: anything */
  };
  IDL_VPTR	arg[2], result;
  IDL_VARIABLE	lun;
  fzHead	fh;
  IDL_FILE_STAT	stat_blk;
  FILE	*fin;
  char	*p, *p2;
  int	wwflag = 0, nelem, type, nb, n, i, flag = 0, nq, mq, sbit, iq;
  IDL_LONG	dim[IDL_MAX_ARRAY_DIM];
  struct compresshead {
    int     tsize, nblocks, bsize;
    byte    slice_size, type; } ch;
  union { int i;  byte b[4];} lmap;
  void	*q;
  
  IDL_EzCall(argc, argv, arg_struct); /* check arguments */
  
  arg[0] = &lun;
  lun.type = IDL_TYP_LONG;
  lun.flags = 0;
  IDL_FileGetUnit(1, arg);	/* get a logical unit */
  arg[1] = argv[1];		/* the file name */
  if (IDL_FileOpen(2, arg, NULL, IDL_OPEN_R, 0, 0)) {
    IDL_FileStat(arg[0]->value.l, &stat_blk);
    fin = stat_blk.fptr;
    if (ck_synch_hd(fin, &fh, &p, &wwflag) < 0)
      IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_RET, "Invalid file header");
    else {
      type = fh.datyp;
      nelem = 1;		/* compute size of array */
      for (i = 0; i < fh.ndim; i++)
	nelem *= fh.dim[i];
      nb = nelem*ana_type_size[type];
      if (argc > 2)
	IDL_VarCopy(IDL_StrToSTRING(p), argv[2]); /* store in <header> */
      else {			/* print the header */
	p2 = p;
	while (*p2) {
	  if (isprint(*p2))
	    putchar(*p2);
	  else if (*p2 == 13)
	    putchar('\n');
	  else
	    putchar(' ');
	  p2++;
	}
	putchar('\n');
      }
      
      /* prepare output */
      for (i = 0; i < fh.ndim; i++)
	dim[i] = fh.dim[i];
      switch (type) {
	case 0:
	  type = IDL_TYP_BYTE;
	  break;
	case 1:
	  type = IDL_TYP_INT;
	  break;
	case 2:
	  type = IDL_TYP_LONG;
	  break;
	case 3:
	  type = IDL_TYP_FLOAT;
	  break;
	case 4:
	  type = IDL_TYP_DOUBLE;
	  break;
      }
      q = IDL_MakeTempArray(type, fh.ndim, dim, IDL_BARR_INI_NOP, &result);
      
      /* read the stuff */
      if ((fh.subf & 1) == 1) {	/* compression case */
	/* read in the compression header */ 
	nq = fread(&ch, 1, 14, fin);
	if (nq != 14 && !flag) 
	  perror("error reading in compression header");
	if (!little_endian) {
	  swapl(&ch.tsize, 1);
	  swapl(&ch.nblocks, 1);
	  swapl(&ch.bsize, 1);
	}
	for (i = 0; i < 4; i++)
	  lmap.b[i] = fh.cbytes[i];
	if (!little_endian)
	  swapl(&lmap.i, 1);
	mq = ch.tsize - 14;
	p = malloc(mq);
	nq = fread(p, 1, mq, fin);
	if (nq != mq) {
	  perror("error reading in compressed data");
	  printf("expected %1d bytes; found %1d\n", mq, nq);
	}
	sbit = ch.slice_size;
	/* fix a problem with ch.nblocks */
	if (ch.bsize * ch.nblocks > nelem) {
	  if (!flag)
	    printf("warning, bad ch.nblocks = %d\n", ch.nblocks);
	  ch.nblocks = nelem / ch.bsize;
	  if (!flag)
	    printf("correcting to %d, hope this is right!\n", ch.nblocks);
	}
	/* some consistency checks */
	if (ch.type % 2 == type)
	  puts("inconsisent compression type");
	switch (ch.type) {
	  case 0:
	    iq = anadecrunch((byte *) p, q, sbit, ch.bsize, ch.nblocks);
	    break;
	  case 1:
	    iq = anadecrunch8((byte *) p, q, sbit, ch.bsize, ch.nblocks);
	    break;
	  case 2:
	    iq = anadecrunchrun((byte *) p, q, sbit, ch.bsize, ch.nblocks);
	    break;
	  case 3:
	    iq = anadecrunchrun8((byte *) p, q, sbit, ch.bsize, ch.nblocks);
	    break;
	  default:
	    if (!flag)
	      printf("error in data type for compressed data, type = %d\n",
		     type);
	    iq = -1;
	    break;
	}
      } else {			/* non-compression case */
	nq = (((fh.subf >> 7) & 1) ^ !little_endian);
	/* just read into the array */
	iq = fread(q, 1, nb, fin);
	if (iq != nb)
	  puts("Error reading from FZ file");
	if (nq)
	  endian(q, nb, type);	/* to host machine native format */
      }
      if (p == fh.extratxt)
	free(p);
      IDL_VarCopy(result, argv[0]);
    }
    IDL_FileFreeUnit(1, arg);
  }
      
  IDL_EzCallCleanup(argc, argv, arg_struct);
} 
/*------------------------------------------------------------------------- */
int ck_synch_hd(FILE *fin, fzHead *fh, char **p, int *wwflag)
/* internal utility */
/* reads the start of an fz file, checks for synch and a reasonable header,
   returns -1 if something amiss, used by several fz readers.
   NOTE: it is assumed that fh is at the start of a block of memory of
   at least 8192 bytes great so that long headers can be stored! */
{
  int	i;
  int	synch_ok, synch_reverse;

  if (fread(fh, 1, 512, fin) != 512)
    perror("fzread in header");

  if (little_endian) {
    synch_ok = 0x5555aaaa;
    synch_reverse = 0xaaaa5555;
  } else {
    synch_ok = 0xaaaa5555;
    synch_reverse = 0x5555aaaa;
  }
  if (fh->synch_pattern != synch_ok) {
    if (fh->synch_pattern == synch_reverse) {
      printf("Reversed F0 magic number -- keep your fingers crossed.\n");
      *wwflag = 1;
    } else {
      printf("File does not have F0 magic number, but %x\n",
	     fh->synch_pattern);
      return -1;
    }
  } else
    *wwflag = 0;		/* sameendian */
  /* if the header is long, read in the rest now */
  i = fh->nhb;
  if (i > 1) {
    if (i > 16) {
      printf("header too long (up to %1d bytes) -- truncated to 7936 bytes\n",
	     512*(i - 1) + 256);
      i = 16;
    }
    fh->extratxt = malloc(512*(i - 1) + 256);
    if (!fh->extratxt)
      puts("Memory allocation error in ch_synch_hd()");
    else {
      memcpy(fh->extratxt, fh->txt, 256);
      fread(fh->extratxt + 256, 1, 512*(i - 1), fin);
      if (i == 16)
	fh->extratxt[i*512+255] = '\0'; /* ensure proper end */
    }
    *p = fh->extratxt;
  } else {
    fh->extratxt = NULL;
    *p = fh->txt;
  }
  if (!little_endian)
    endian(&fh->dim[0], fh->ndim*sizeof(int), 2);
  return 1;
}
/*------------------------------------------------------------------------- */
void ana_fzwrite(int argc, IDL_VPTR argv[])
/* call: FZWRITE, array, filename [, header] */
/* <array>: data to be stored */
/* <filename>: name of the uncompressed FZ file */
/* <header>: optional file header */
{
  static IDL_EZ_ARG arg_struct[] = {
    { IDL_EZ_DIM_ARRAY, IDL_TYP_NUMERICAL, IDL_EZ_ACCESS_R, IDL_TYP_UNDEF,
      0, 0 },			/* <array>: numerical array */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_MASK(IDL_TYP_STRING), IDL_EZ_ACCESS_R,
      IDL_TYP_UNDEF, 0, 0 },	/* <filename>: name of the file */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_MASK(IDL_TYP_STRING), IDL_EZ_ACCESS_R,
      IDL_TYP_UNDEF, 0, 0 }	/* <header>: string */
  };
  IDL_VPTR	arg[2];
  IDL_VARIABLE	lun;
  fzHead	fh;
  int	idl_type, type, n, i, mq, j;
  IDL_LONG	*dim;
  char	*p, *q;
  FILE	*fout;
  IDL_FILE_STAT	stat_blk;

  IDL_EzCall(argc, argv, arg_struct); /* check arguments */
  
  arg[0] = &lun;
  lun.type = IDL_TYP_LONG;
  lun.flags = 0;
  IDL_FileGetUnit(1, arg);	/* get a logical unit */
  arg[1] = argv[1];		/* the file name */
  if (IDL_FileOpen(2, arg, NULL, IDL_OPEN_W, 0, 0)) {
    IDL_FileStat(arg[0]->value.l, &stat_blk);
    fout = stat_blk.fptr;
    fh.synch_pattern = little_endian? 0x5555aaaa: 0xaaaa5555; /* magic # */
    fh.subf = (!little_endian) << 7; /* no compression, indicate endian */
    fh.source = 0;
    fh.nhb = 1;		/* may be changed later */
    idl_type = argv[0]->type;
    switch (idl_type) {
      case IDL_TYP_BYTE:
	type = 0;
	break;
      case IDL_TYP_INT:
	type = 1;
	break;
      case IDL_TYP_LONG:
	type = 2;
	break;
      case IDL_TYP_FLOAT:
	type = 3;
	break;
      case IDL_TYP_DOUBLE:
	type = 4;
	break;
    }
    fh.datyp = type;
    fh.ndim = argv[0]->value.arr->n_dim;
    n = argv[0]->value.arr->arr_len;
    dim = argv[0]->value.arr->dim;
    for (i = 0; i < fh.ndim; i++)
      fh.dim[i] = dim[i];
    if (!little_endian)
      endian(fh.dim, fh.ndim*sizeof(int), 2);
    
    /* check header */
    if (argc > 2) {		/* we have a header */
      p = argv[2]->value.str.s;
      mq = argv[2]->value.str.slen;	/* add one for terminating null */
      if (mq > 256) {
	j = (mq - 256)/512 + 2;	/* number of header blocks */
	if (j > 16) {
	  puts("Header for FZ file is more than 7935 bytes long; truncating");
	  j = 16;
	  mq = 7935;
	}
	fh.nhb = j;
      }
      if (fh.nhb == 1) {	/* only the single block */
	memcpy(fh.txt, p, mq);
	j = fwrite(&fh, 1, 512, fout);
      } else {			/* more blocks */
	/* we want to write a number of complete 512-byte blocks; */
	/* we take the first 256 from <fh>, and then as many complete */
	/* 512-byte blocks as fit in <p>.  Then we pad with zeros to fill */
	/* up the last 512-byte block. */
        j = fwrite(&fh, 1, 256, fout);/* data info block */
	j += fwrite(p, 1, mq, fout); /* header text */
	i = mq - (256 + (fh.nhb - 2)*512) + 1; /* # chars for last block */
					       /* (add 1 for terminating \0) */
	while (i--)
	  j += (putc(0, fout) != EOF);
	if (j != mq + 1)
	  printf("Wrote %d bytes of the header into an FZ file; wanted %d\n",
		 j, mq + 1);
      }
    } else {			/* no explicit header */
      for (i = 0; i < 256; i++)
	fh.txt[i] = '\0';
      j = fwrite(&fh, 1, 512, fout);
      if (j != 512)
	printf("Wrote %d bytes of the header into an FZ file; wanted %d\n",
	       j, 512);
    }
    /* write the data */
    j = fwrite(argv[0]->value.arr->data, 1, n, fout);
    if (j != n)
      printf("Wrote %d bytes of data into an FZ file; wanted %d\n",
	     j, n);
    IDL_FileClose(1, arg, NULL);
    IDL_FileFreeUnit(1, arg);    
  }

  IDL_EzCallCleanup(argc, argv, arg_struct);
}
/*------------------------------------------------------------------------- */
void ana_fcwrite(int argc, IDL_VPTR argv[])
/* call: FCWRITE, array, filename [, header] */
/* <array>: data to be stored */
/* <filename>: name of the compressed FZ file */
/* <header>: optional file header */
{
  static IDL_EZ_ARG arg_struct[] = {
    { IDL_EZ_DIM_ARRAY, IDL_TYP_MASK(IDL_TYP_BYTE)
      | IDL_TYP_MASK(IDL_TYP_INT), IDL_EZ_ACCESS_R, IDL_TYP_UNDEF,
      0, 0 },			/* <array>: numerical array */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_MASK(IDL_TYP_STRING), IDL_EZ_ACCESS_R,
      IDL_TYP_UNDEF, 0, 0 },	/* <filename>: name of the file */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_MASK(IDL_TYP_STRING), IDL_EZ_ACCESS_R,
      IDL_TYP_UNDEF, 0, 0 },	/* <header>: string */
    { IDL_EZ_DIM_MASK(0), IDL_TYP_B_ALL, IDL_EZ_ACCESS_R,
      IDL_TYP_LONG, 0, 0 }	/* <crunch_slice>: scalar */
  };
  IDL_VPTR	arg[2];
  IDL_VARIABLE	lun;
  fzHead	fh;
  int	idl_type, type, n, i, mq, j, limit, nx, ny, iq;
  static int	crunch_slice = 5;
  IDL_LONG	*dim;
  char	*p, *q;
  FILE	*fout;
  IDL_FILE_STAT	stat_blk;
  void	*q1, *q2;
  union { int i;  byte b[4];} lmap;

  IDL_EzCall(argc, argv, arg_struct); /* check arguments */
  
  if (argc > 3)
    crunch_slice = argv[3]->value.l;

  arg[0] = &lun;
  lun.type = IDL_TYP_LONG;
  lun.flags = 0;
  IDL_FileGetUnit(1, arg);	/* get a logical unit */
  arg[1] = argv[1];		/* the file name */
  if (IDL_FileOpen(2, arg, NULL, IDL_OPEN_W, 0, 0)) {
    IDL_FileStat(arg[0]->value.l, &stat_blk);
    fout = stat_blk.fptr;
    p = (char *) &fh;
    i = 512;
    while (i--)
      *p++ = '\0';
    fh.synch_pattern = little_endian? 0x5555aaaa: 0xaaaa5555; /* magic # */
    fh.subf = ((!little_endian) << 7) | 1; /* compression, indicate endian */
    fh.source = 0;
    fh.nhb = 1;		/* may be changed later */
    idl_type = argv[0]->type;
    switch (idl_type) {
      case IDL_TYP_BYTE:
	type = 0;
	break;
      case IDL_TYP_INT:
	type = 1;
	break;
      case IDL_TYP_LONG:
	type = 2;
	break;
      case IDL_TYP_FLOAT:
	type = 3;
	break;
      case IDL_TYP_DOUBLE:
	type = 4;
	break;
    }
    fh.datyp = type;
    fh.ndim = argv[0]->value.arr->n_dim;
    dim = argv[0]->value.arr->dim;
    for (i = 0; i < fh.ndim; i++)
      fh.dim[i] = dim[i];
    nx = fh.dim[0];
    endian(fh.dim, fh.ndim*sizeof(int), 2);
    n = argv[0]->value.arr->arr_len;
    ny = argv[0]->value.arr->n_elts/nx;
    q1 = argv[0]->value.arr->data;
    
    limit = 2*n;
    q2 = malloc(limit);
    switch (type) {
      case 0:
	iq = anacrunch8(q2, q1, crunch_slice, nx, ny, limit);
	break;
      case 1:
	iq = anacrunch(q2, q1, crunch_slice, nx, ny, limit);
	break;
    }
    if (iq < 0)
      IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_RET,
		  "Insufficient compression");
    else {
      lmap.i = iq;
      if (!little_endian)
	endian(&lmap.i, sizeof(int), 2);
      memcpy(&fh.cbytes, &lmap.i, sizeof(int));

      /* check header */
      if (argc > 2) {		/* we have a header */
	p = argv[2]->value.str.s;
	mq = argv[2]->value.str.slen;	/* add one for terminating null */
	if (mq > 256) {
	  j = (mq - 256)/512 + 2;	/* number of header blocks */
	  if (j > 16) {
	    puts("Header for FZ file is more than 7935 bytes long; truncating");
	    j = 16;
	    mq = 7935;
	  }
	  fh.nhb = j;
	}
	if (fh.nhb == 1) {	/* only the single block */
	  memcpy(fh.txt, p, mq);
	  j = fwrite(&fh, 1, 512, fout);
	} else {			/* more blocks */
	  /* we want to write a number of complete 512-byte blocks; */
	  /* we take the first 256 from <fh>, and then as many complete */
	  /* 512-byte blocks as fit in <p>.  Then we pad with zeros to fill */
	  /* up the last 512-byte block. */
	  j = fwrite(&fh, 1, 256, fout);/* data info block */
	  j += fwrite(p, 1, mq, fout); /* header text */
	  i = mq - (256 + (fh.nhb - 2)*512) + 1; /* # chars for last block */
	  /* (add 1 for terminating \0) */
	  while (i--)
	    j += (putc(0, fout) != EOF);
	  if (j != mq + 1)
	    printf("Wrote %d bytes of the header into an FZ file; wanted %d\n",
		   j, mq + 1);
	}
      } else {			/* no explicit header */
	for (i = 0; i < 256; i++)
	  fh.txt[i] = '\0';
	j = fwrite(&fh, 1, 512, fout);
	if (j != 512)
	  printf("Wrote %d bytes of the header into an FZ file; wanted %d\n",
		 j, 512);
      }
      if (fwrite(q2, 1, iq, fout) != iq)
	IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_RET,
		    "Error writing compressed data into FZ file");
      free(q2);
    }

    IDL_FileClose(1, arg, NULL);
    IDL_FileFreeUnit(1, arg);    
  }

  IDL_EzCallCleanup(argc, argv, arg_struct);
}
/*------------------------------------------------------------------------- */