/* anarw.c, standalone C functions for reading and writing f0 type data files */
#include<stdio.h>
 typedef unsigned char byte;
 int	ana_type_size[5]={1,2,4,4,8};
 static	int	runlengthflag=0, crunch_slice, crunch_bpp, crunch_bits;
 /*------------------------------------------------------------------------- */
 /* stand alone routines for reading and writing ana style files
 the format for the uncompressed is simple enough for an idl routine
 but the compressed is more complicated */
 struct fzhead {				/* first block for fz files */
  int     synch_pattern;
  byte    subf;
  byte    source;
  byte    nhb,datyp,ndim,free1;
  byte    cbytes[4];      /* can't do as int because of %4 rule*/
  byte    free[178];
  int     dim[16];
  char    txt[768];};
  /* note that this fzhead structure has room for 2*512 bytes to allow for
  long comments in some files (though rather rare) */
 /*------------------------------------------------------------------------- */
unsigned char *fzread(name, dtype, num_dim, dims, s)	/* fzread subroutine */
 /* read standard f0 files, compressed or not */
 char	*name, **s;
 int	*dtype, *num_dim, dims[8];
 /* name is a string containing file name, the  return value will be a
 pointer to the read in array, num_dim is the number of dimensions of the
 array, type is the data type (0=byte, 1=i*2, 2=i*4, 3=f*4, 4=f*8), and
 dims is an int array with the actual dimensions (up to 8, so make it 8 long)
 s is a pointer to the comment string
 note that dtype, num_dim, dims, s are returned parameters and should be
 called by ref; i.e., a sample call is:
 
 char	fname[100], *comments;
 byte	*data;
 int	dims[8], num_dim, type;
 data = fzread(fname, &type, &num_dim, dims, &comments);

 */
 {
 int	iq, n, nb, nc, nd, j, i, mq, nq, sbit, nelem, wwflag=0, type;
 char	*p, *q;
 struct	fzhead	*fh, fzheader;
 FILE *fopen(),*fin;
 
 struct compresshead {
		 int     tsize,nblocks,bsize;
		 byte    slice_size,type; } ch;
 
 union { int i;  byte b[4];} lmap;
						 /* try to open the file */
 if ((fin=fopen(name,"r")) == NULL) { perror("file open error");
 	return 0; }
 fh = &fzheader;  /* done mainly to keep code similar to ana module */
				 		/* read in header block */
 if (fread(fh,1,512,fin) != 512) { perror("fzread in header");
 	return 0; }
 type = fh->datyp;
#if LITTLEENDIAN
 if (fh->synch_pattern != 0x5555aaaa) {
  if (fh->synch_pattern == 0xaaaa5555) {
#else
 if (fh->synch_pattern != 0xaaaa5555) {
  if (fh->synch_pattern == 0x5555aaaa) {
#endif
  printf("reversed F0 synch pattern, assuming wwflag\n"); wwflag = 1; } else {
  printf("file does not have the F0 synch pattern %x\n",fh->synch_pattern); 
  fclose(fin);
  return 0; } }
#if LITTLEENDIAN
#else
 swapl(&fh->dim[0],(int) fh->ndim);	/* for BIGENDIAN machines (SGI) */
#endif
 /* compute size of array */
 nelem=1;
 for (i=0;i<fh->ndim;i++) { nelem *= fh->dim[i];}
 nb = nelem * ana_type_size[type];
 /* if the header is long, read in the rest now */
 if (fh->nhb > 1 ) { if (fh->nhb > 2) {
	 printf("header more than 2 blocks ! Cannot handle it!\n");
	 fclose(fin);
	 return 0; }
 fread((char *) fh + 512,1,512*(fh->nhb-1),fin); }

 p = (char *) fh->txt; q = p;
 n = 256 + 512*(fh->nhb-1); nc = 0;
   while (n--) { if ( *p++ == 0 ) break;  nc++; }	/* count to null */
			 		/* malloc space for the comment */
 *s = p = (char *) malloc(nc+1);
 if (p == NULL) {printf("could not malloc %d bytes for comment\n", (nc+1));
	 fclose(fin);
	 return 0; }
 while (nc--) *p++ = *q++;	*p = 0;			/*load it */
						 /* create the output array */
						 /* nb is size in bytes */
 q = (char *) malloc(nb);
 if (q == NULL) {printf("could not malloc %d bytes for array\n", nb);
	 fclose(fin);
	 return 0; }
				 /* big branch on compression flag */
 if (fh->subf%128 == 1 )      {                       /* compression case */
 /* read in the compression header */ 
 nq=fread(&ch,1,14,fin);
 if (nq != 14) { perror("error reading in compression header"); }
#if LITTLEENDIAN
#else
 swapl(&ch.tsize,1);     swapl(&ch.nblocks,1);   swapl(&ch.bsize,1);
#endif
 for (i=0;i<4;i++) lmap.b[i]=fh->cbytes[i];
#if LITTLEENDIAN
#else
 swapl(&lmap.i,1);
#endif
 mq = ch.tsize-14;
 p = (char *) malloc(mq);
 if (p == NULL) {printf("could not malloc %d bytes for compressed data\n", mq);
	 fclose(fin);
	 return 0; }
 nq=fread(p, 1, mq, fin);
 if (nq != mq) { perror("error reading in compressed data");
	 fclose(fin);
	 return 0; }
 sbit=ch.slice_size;
 /* fix a problem with ch.nblocks on some old files*/
 if ( ch.bsize * ch.nblocks > nelem ) {
  printf("warning, bad ch.nblocks = %d\n", ch.nblocks);
  ch.nblocks = nelem / ch.bsize;
  printf("correcting to %d, hope this is right!\n", ch.nblocks);
  }
 /* some consistency checks */
 if (ch.type%2 == type) { printf("inconsistent compression type\n");
	 fclose(fin);
	 return 0; }
 switch (ch.type) {
 case 0: iq = anadecrunch(p,q,sbit,ch.bsize,ch.nblocks); break;
 case 1: iq = anadecrunch8(p,q,sbit,ch.bsize,ch.nblocks); break;
 case 2: iq = anadecrunchrun(p,q,sbit,ch.bsize,ch.nblocks); break;
 case 3: iq = anadecrunchrun8(p,q,sbit,ch.bsize,ch.nblocks); break;
 /* added 32 bit decompression 2/12/96 based on files.c */
 case 4: iq = anadecrunch32(p,q,sbit,ch.bsize,ch.nblocks); break;
 default: printf("error in data type for compressed data, fh->datyp =%d\n",type);
 	iq = -1; break;
 }
 if (iq != 1)	 { fclose(fin);	return 0; }
 } else {                                        /* non-compression case */
						 /* just read into the array */
 /* the top bit of the byte fh->subf denotes endian type, 0 for little
 and 1 for big */
 if (fread(q, 1, nb, fin) != nb) {
	 perror("error reading compression header");  fclose(fin);  return 0; }
#if defined(LITTLEENDIAN)
 if (type == 1 && (fh->subf >= 128))  swapb((char *)q, nb);
 if (type == 2 && fh->subf >= 128)  swapl(q, nb/4);
 if (type == 3 && fh->subf >= 128)  swapl(q, nb/4);
 if (type == 4 && fh->subf >= 128)  swapd(q, nb/8);
#else
 if (type == 1 && (fh->subf < 128 || wwflag == 1) )  swapb((char *)q, nb);
 if (type == 2 && fh->subf < 128)  swapl(q, nb/4);
 if (type > 2 && fh->subf < 128 ) {
	printf("warning, might be unsupported floating point data from VMS\n");
	swapl(q, nb/4); }
#endif
 }
 /* common section again */
 /* set the output arguments describing the array */
 *dtype = type;
 *num_dim = fh->ndim;
 for (i=0;i<fh->ndim;i++) { dims[i] = fh->dim[i];}
 fclose(fin);
 return (unsigned char *) q;
 }
 /*------------------------------------------------------------------------- */
int fzwrite(array, name, type, num_dim, dims, s)/* fzwrite subroutine */
 /* write standard f0 files, uncompressed */
 char	*name, *s;
 char	*array;
 int	type, num_dim, dims[8];
 /* arguments similar to fzread except for the array (a pointer) added at
 the beginning and all arguments must be defined for the call and are passed
 by value
 name is a string containing file name,
 num_dims is the number of dimensions of the
 array, type is the data type (0=byte, 1=i*2, 2=i*4, 3=f*4, 4=f*8), and
 dims is an int array with the actual dimensions (up to 8, so make it 8 long)
 s is a pointer to the comment string, returns 1 if OK
 
 a sample call is:
 
 fzwrite(data, fname, type, num_dim, dims, comments);

 */
 {
 int	iq, n, nc, nd, j, i, mq, nq, sbit;
 char	*p, *q;
 struct	fzhead	*fh, fzheader;
 FILE *fopen(),*fout;
 nd = num_dim;
						 /* try to open the file */
 if ((fout=fopen(name,"w")) == NULL) { perror("file open error");
 	return 0; }
					 /* use scrat to setup header */
 fh = &fzheader;  /* done mainly to keep code similar to ana module */
 bzero( (char *) fh, 1024);		/*zero it first */
 					/* note 1024 size here */
			 /* these have to be readable by the Vax's */
#if defined(LITTLEENDIAN)
 fh->synch_pattern = 0x5555aaaa;	fh->subf = 0;	fh->source = 0;
#else		/* big endian option */
 fh->synch_pattern = 0xaaaa5555;	fh->subf = 128;	fh->source = 0;
#endif
 fh->nhb = 1;					/*may be changed later */
 fh->datyp = type;	fh->ndim = nd;
 n = ana_type_size[type];
 for(j=0;j<nd;j++) { n *= dims[j];  fh->dim[j] = dims[j]; }
#if LITTLEENDIAN
#else
 swapl(&fh->dim[0],(int) fh->ndim);	/* for BIGENDIAN machines (SGI) */
#endif
					 /* check for text to put in header */
   p = s;
   mq = (int) strlen(s) + 1;
   q = (char *) fh->txt;
   if ( mq > 256 ) {
   j = (mq-256)/512 + 2; if (j > 2) {
     printf("header more than 2 blocks ! Cannot handle it!\n");
     fclose(fout); return 0; }
   fh->nhb = j;  }
   while (mq--) *q++ = *p++;

 j=fwrite(fh, 1, 512 * fh->nhb, fout);		/*write header */
 if (j != 512 * fh->nhb) {
	 perror("error writing header");	fclose(fout);	return 0; }
 
 printf("n = %d\n",n);
 q = array;
 if (fwrite(q, 1, n, fout) != n) {
	 perror("error writing file");	fclose(fout);	return 0; }
 fclose(fout);
 return	1;
 }
 /*------------------------------------------------------------------------- */
int fcrunwrite(array, name, type, num_dim, dims, s)/* fcrunwrite subroutine */
 /* write standard f0 files, compressed format with additional
 run length encoding which helps when there are large blank areas*/
 char	*name, *s;
 char	*array;
 int	type, num_dim, dims[8];
 {
 /* we just set the runlengthflag to 1 and call fcwrite */
 int	iq;
 runlengthflag = 1;
 iq = fcwrite(array, name, type, num_dim, dims, s);
 runlengthflag = 0;	/* re-set to default */
 return iq;
 }
 /*------------------------------------------------------------------------- */
int fcwrite(array, name, type, num_dim, dims, s)/* fcwrite subroutine */
 /* write standard f0 files, compressed format */
 char	*name, *s;
 char	*array;
 int	type, num_dim, dims[8];
 {
 int	iq, n, nc, nd, j, i, mq, nq, sbit, nx, ny, limit;
 char	*p, *q;
 struct	fzhead	*fh, fzheader;
 union { int i;  byte b[4];} lmap;
 FILE *fopen(),*fout;
 nd = num_dim;
 if (type > 1) { printf("compression supported only for I*1 and I*2\n");
 	return 0; }
 /* reasonable defaults */
 if (type ==0 ) crunch_slice=2; else crunch_slice=5;
						 /* try to open the file */
 if ((fout=fopen(name,"w")) == NULL) { perror("file open error");
 	return 0; }
					 /* use scrat to setup header */
 fh = &fzheader;  /* done mainly to keep code similar to ana module */
 bzero( (char *) fh, 1024);		/*zero it first */
 					/* note 1024 size here */
			 /* these have to be readable by the Vax's */
#if LITTLEENDIAN
 fh->synch_pattern = 0x5555aaaa;	fh->subf = 1;	fh->source = 0;
#else
 fh->synch_pattern = 0xaaaa5555;	fh->subf = 129;	fh->source = 0;
#endif
 fh->nhb = 1;					/*may be changed later */
 fh->datyp = type;	fh->ndim = nd;
 n = 1;
 for(j=0;j<nd;j++) { n *= dims[j];  fh->dim[j] = dims[j]; }
 nx = fh->dim[0];	ny = n/nx;
 n = n * ana_type_size[type];
#if LITTLEENDIAN
#else
 swapl(&fh->dim[0],(int) fh->ndim);	/* for BIGENDIAN machines (SGI) */
#endif
					 /* check for text to put in header */
   p = s;
   mq = (int) strlen(s) + 1;
   q = (char *) fh->txt;
   if ( mq > 256 ) {
   j = (mq-256)/512 + 2; if (j > 2) {
     printf("header more than 2 blocks ! Cannot handle it!\n");
     fclose(fout); return 0; }
   fh->nhb = j;  }
   while (mq--) *q++ = *p++;

			 /* now compress the array, must be a byte or short */
 p = array;
 limit = 2*n;
 q = (char *) malloc(limit);	/* some room for mistakes */
 switch (type) {
 case 0:
 if (runlengthflag == 0)
 iq = anacrunch8( q, p, crunch_slice, nx, ny, limit); else
 iq = anacrunchrun8( q, p, crunch_slice, nx, ny, limit); break;
 case 1:
 if (runlengthflag == 0)
 iq = anacrunch( q, p, crunch_slice, nx, ny, limit); else
 iq = anacrunchrun( q, p, crunch_slice, nx, ny, limit);  break;
 }
 if ( iq < 0 ) {
 printf("not enough space allocated (%d bytes) for compressed array\n",limit);
 free(q);	fclose(fout); return 0; }
 /* nice to know these values */
 crunch_bpp = 8.0 * iq / (nx * ny);
 crunch_bits = 8 * iq;

 lmap.i = iq;
#if LITTLEENDIAN
#else
 swapl(&lmap.i,1);
#endif
 for (i=0;i<4;i++) fh->cbytes[i]=lmap.b[i];
 j=fwrite(fh, 1, 512 * fh->nhb, fout);		/*write header */
 if (j != 512 * fh->nhb) {	free(q);
	 perror("error writing header");	fclose(fout);	return 0; }
 if (fwrite(q, 1, iq, fout) != iq) {	free(q);
	 perror("error writing data");	fclose(fout);	return 0; }
 free(q);	fclose(fout);
 return	1;
 }
 /*------------------------------------------------------------------------- */
unsigned char *cread(fin, size)	/* cread subroutine */
 /* decompresses the compressed part of a compressed FITS file */
 FILE   *fin;
 int	*size;
 {
 int	iq, n, nb, nc, nd, j, i, mq, nq, sbit, nelem, type;
 char	*p, *q;
 
 struct compresshead {
		 int     tsize,nblocks,bsize;
		 byte    slice_size,type; } ch;
 
 nq=fread(&ch,1,14,fin);
 if (nq != 14) { perror("error reading in compression header"); }
#if LITTLEENDIAN
#else
 swapl(&ch.tsize,1);     swapl(&ch.nblocks,1);   swapl(&ch.bsize,1);
#endif
 mq = ch.tsize-14;
 p = (char *) malloc(mq);
 if (p == NULL) {printf("could not malloc %d bytes for compressed data\n", mq);
	 fclose(fin);
	 return 0; }
 nq=fread(p, 1, mq, fin);
 if (nq != mq) { perror("error reading in compressed data");
	 fclose(fin);
	 return 0; }
 /* now get memory for uncompressed data */
 switch (ch.type) {
 case 0: type = 1; break;
 case 1: type = 0; break;
 case 2: type = 1; break;
 case 3: type = 0; break;
 /* added 32 bit decompression 2/12/96 based on files.c */
 case 4: type = 2; break;
 }
 /* note that ana type and ch.type are not the same */
 
 mq =  ch.bsize * ch.nblocks * ana_type_size[type];
 *size = mq;
 q = (char *) malloc(mq);
 if (q == NULL) {printf("could not malloc %d bytes for uncompressed data\n", mq);
	 fclose(fin);
	 return 0; }
 sbit=ch.slice_size;
 switch (ch.type) {
 case 0: iq = anadecrunch(p,q,sbit,ch.bsize,ch.nblocks); break;
 case 1: iq = anadecrunch8(p,q,sbit,ch.bsize,ch.nblocks); break;
 case 2: iq = anadecrunchrun(p,q,sbit,ch.bsize,ch.nblocks); break;
 case 3: iq = anadecrunchrun8(p,q,sbit,ch.bsize,ch.nblocks); break;
 /* added 32 bit decompression 2/12/96 based on files.c */
 case 4: iq = anadecrunch32(p,q,sbit,ch.bsize,ch.nblocks); break;
 default: printf("error in data type for compressed data,  ch.type =%d\n",ch.type);
 	iq = -1; break;
 }
 if (iq != 1)	 { fclose(fin);	return 0; }

 /* common section again */
 return (unsigned char *) q;
 }
 /*------------------------------------------------------------------------- */