/*
   program to read BYU SIR image files and write out as .bmp file

   Written by DGL Aug. 8, 1997

   Written in ANSI C.  Define or undefine SWAP based on machine byte order

*/

#define SWAP 1                /* swap bytes (little endian: dec, pc) */
#undef SWAP                   /* no byte swapping (big endian: HP,SUN) */

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

#define REL_BEGIN 0           /* fseek relative to beginning of file */
#define SIR_HEADER_SIZE 512   /* SIR file header size in bytes */
#define BUFFER 6000           /* space in buffer */ 

int read_sir_header(FILE *imf, int *nhead, int *ndes, int *nhtype,
		    int *idatatype, int *nsx, int *nsy, 
		    float *xdeg, float *ydeg, 
		    float *ascale, float *bscale, float *a0, float *b0, 
		    int *ioff, int *iscale, int *iyear, 
		    int *isday, int *ismin, int *ieday, int *iemin, 
		    int *iregion, int *itype, int *iopt,
		    int *ispare1, int *ispare2, int *ispare3,
		    float *anodata, float *v_min, float *v_max,
		    char *sensor, char *title, char *type, char *tag,
		    char *crproc, char *crtime, int maxdes, 
		    char *descrip, int *ldes, int maxi, short *iaopt,
		    int *nia);

void print_head(int nhead, int ndes, int nhtype, int idatatype, int nsx,
		int nsy, float xdeg, float ydeg, float ascale,
		float bscale, float a0, float b0, 
		int ioff, int iscale, int iyear, int isday, int ismin,
		int ieday, int iemin, int iregion, int itype, int iopt,
		int ipol, int ifreqhm, int ispare1,
		float anodata, float v_min, float v_max,
		char *sensor, char *title, char *type, char *tag, char
		*crproc, char *crtime, 
		char *descrip, int ldes, short *iaopt, int nia);

void swap(short *, int);

#ifndef SWAP
void swapi(int *, int);
#endif

int writebmp(char *name_in, char *image, int *nw, int *nh,
	     char *rmap, char *gmap, char *bmap, int *numcols);


int main(int argc, char **argv)
{
  FILE *imf;
  int i, j, k;
  int ix, iy, ix1, iy1;
  float x, y, x1, y1, thelon, thelat, alon, alat, am, amag;
  char outfname[80];

  char rtab[256], gtab[256], btab[256];
	
  float smin, smax;
  float s, soff;
  float scale, scaleoffset;

/* other SIR file header information */

  float v_min, v_max, anodata;
  int nsx, nsy, ioff, iscale, iyear, isday, ismin, ieday, iemin;
  int iregion, itype, nhead, ndes, nhtype, idatatype, ldes, nia;
  int ipol, ifreqhm, ispare1;
  char title[101], sensor[41], crproc[101], type[139], tag[101], crtime[29];
  float xdeg, ydeg, ascale, bscale, a0, b0;
  int iopt;

#define MAXDES 1024
  char descrip[MAXDES+1];
#define MAXI 128
  short iaopt[MAXI];

  float *stval;

  int ierr;
  short *a;
  char *data;
  int nlines, width;

  union 
  {
    short i2[2];
    float f2;
  } un;

  fprintf(stdout,"BYU SIR to BMP conversionprogram\n") ;
  if(argc < 3) {
    fprintf(stdout,"\nusage: %s sir_in bmp_out <dmin> <dmax> <ctab>\n\n",argv[0]);
    fprintf(stdout," input parameters:\n");
    fprintf(stdout,"   sir_in  =  input SIR file\n");
    fprintf(stdout,"   bmp_out =  output bmp file\n");
    fprintf(stdout,"   dmin    =  min saturation value (optional,def=vmin)\n");
    fprintf(stdout,"   dmax    =  max saturation value (optional,def=vmax)\n");
    fprintf(stdout,"   ctab    =  color table file (optional,def=greyscale)\n"); 
    return(0);
  }

  imf = fopen(argv[1],"r"); 
  if (imf == NULL) {
     fprintf(stdout,"ERROR: cannot open input file: %s\n",argv[1]); 
     exit(-1);
  }

/* get SIR image header information */

  ierr = read_sir_header(imf, &nhead, &ndes, &nhtype,
		    &idatatype, &nsx, &nsy, 
		    &xdeg, &ydeg, &ascale, &bscale, &a0, &b0, 
		    &ioff, &iscale, &iyear, &isday, &ismin, &ieday, &iemin, 
		    &iregion, &itype, &iopt, &ipol, &ifreqhm, &ispare1,
		    &anodata, &v_min, &v_max,
		    sensor, title, type, tag, crproc, crtime, 
		    MAXDES, descrip, &ldes, MAXI, iaopt, &nia);  
  if (ierr < 0) {
     fprintf(stdout,"*** Error reading SIR header from file '%s'\n",argv[1]);
     exit(-1);
  }

/* write out SIR file header information */

  fprintf(stdout,"\nSIR file header: '%s'\n",argv[1]);
  print_head(nhead, ndes, nhtype, idatatype, nsx, nsy, 
	     xdeg, ydeg, ascale, bscale, a0, b0, 
	     ioff, iscale, iyear, isday, ismin, ieday, iemin, 
	     iregion, itype, iopt, ipol, ifreqhm, ispare1,
	     anodata, v_min, v_max,
	     sensor, title, type, tag, crproc, crtime, 
	     descrip, ldes, iaopt, nia);  

  rewind(imf);

  if (nhtype < 20) {  /* old style header */
    /*    fprintf(stdout," (Old style header) %d\n",nhtype); */
    v_min=-32.0;
    v_max=0;
    anodata=-32.0;
    idatatype=2;
  }

  smin = v_min;
  smax = v_max;

  if (argc > 3) sscanf(argv[3],"%f",&smin); 
  if (argc > 4) sscanf(argv[4],"%f",&smax); 
  if (smax - smin == 0.0) smax = smin + 1.0;

  width = nsx;
  nlines = nsy;

  a =(short *) malloc(sizeof(short)*width);
  data = (char *) malloc(sizeof(char)*width*nlines);
  stval = (float *) malloc(sizeof(float)*width*nlines);
  if (a == NULL || data == NULL || stval == NULL) {
     fprintf(stderr,"*** ERROR: memory allocation failure...\n");
     exit(-1);
  }

  s = 1.0 / ((float) iscale);
  soff = 32767.0/ ((float) iscale) + ((float) ioff);
  if (idatatype == 1) soff = 128.0/ ((float) iscale) + ((float) ioff);

  scale = (smax-smin);
  if (scale > 0.) scale = 255./ scale;
  scaleoffset = smin;

/* read file, converting two byte integers to floats them to byte array
   For later use, copy floating point valuses to an array */

  fseek(imf,SIR_HEADER_SIZE*nhead,REL_BEGIN);

  k = 0;
  for (i = 0; i < nlines; i++) {

    switch(idatatype) {
    case 1:     /* bytes */
      fread((char *)a, sizeof(char), width, imf) ;
      for (j = 0; j < width; j++){
	amag = ((float) a[j]) * s + soff; /* pixel value */
	*(stval+k*width+j) = amag;              /* floating point array */
      }
      break;
      
    case 4:     /* floats */
#ifdef SWAP
      fread((char *)a, sizeof(short), 2*width, imf) ;
      swap(a,width);
      for (j = 0; j < width; j++){
	un.i2[1]=a[j*2];
	un.i2[0]=a[j*2+1];
	amag = un.f2;
	*(stval+k*width+j) = amag;              /* floating point array */
      }
#else
      fread((char *) (stval+k*width), sizeof(float), width, imf);
#endif
      break;
      
    default:    /* 0 or 2: read two byte integers */
      fread((char *)a, sizeof(short), width, imf) ;
#ifdef SWAP	
      swap(a,width);
#endif
      for (j = 0; j < width; j++){
	amag = ((float) a[j]) * s + soff; /* pixel value */
	*(stval+k*width+j) = amag;              /* floating point array */
      }
    }
    
    for (j = 0; j < width; j++){      
      amag = *(stval+k*width+j);

      /* scale floating point to byte values */
      am = scale * (amag - scaleoffset);
      if(am > 255.) am = 255.;		/* check overflow */
      if(am < 0.) am = 0.;		/* check underflow */
      *(data+k*width+j) = (char)((int)(am));    /* byte array for output*/
    }
    k++;
  }
  fclose(imf);
  fprintf(stdout,"SIR input file successfully read\n");

  fprintf(stdout,"\nBMP Min, Max: %f , %f\n",smin,smax);

  /* default greyscale color table */

  j=255;
  for (i=0; i<256; i++) {
    if (i > 127)
      rtab[i]=i-256;
    else
      rtab[i]=i;
    gtab[i]=rtab[i];
    btab[i]=rtab[i];
  }
  
  if (argc > 5) { /* read user input color table file */
    fprintf(stdout,"Reading input color table file '%s'\n",argv[5]);
    imf = fopen(argv[5],"r"); 
    if (imf == NULL) {
      fprintf(stdout,"*** ERROR: cannot open color table file: %s\n",argv[5]);
      fprintf(stdout,"Using greyscale color table\n");
    } else {
      if (fread(rtab, sizeof(char), 256, imf) == 0) 
	fprintf(stdout," *** Error reading color table file\n");
      if (fread(gtab, sizeof(char), 256, imf) == 0) 
	fprintf(stdout," *** Error reading color table file\n");
      if (fread(btab, sizeof(char), 256, imf) == 0) 
	fprintf(stdout," *** Error reading color table file\n");
      fclose(imf);
    }
  } else
    fprintf(stdout,"Using greyscale color table\n");

  fprintf(stdout,"Writing output bmp file '%s'\n", argv[2]);

  ierr=writebmp(argv[2], data, &nsx, &nsy, rtab, gtab, btab, &j);
  if (ierr < 0) {
     fprintf(stderr,"*** ERROR writing output file ***\n");
     fflush(stderr);
  }
  else
    fprintf(stdout,"output file successfully written\n");

  return(0);
}

	    
int read_sir_header(FILE *imf, int *nhead, int *ndes, int *nhtype,
		    int *idatatype, int *nsx, int *nsy, 
		    float *xdeg, float *ydeg, 
		    float *ascale, float *bscale, float *a0, float *b0, 
		    int *ioff, int *iscale, int *iyear, 
		    int *isday, int *ismin, int *ieday, int *iemin, 
		    int *iregion, int *itype, int *iopt,
		    int *ipol, int *ifreqhm, int *ispare1,
		    float *anodata, float *v_min, float *v_max,
		    char *sensor, char *title, char *type, char *tag,
		    char *crproc, char *crtime, int maxdes, 
		    char *descrip, int *ldes, int maxi, short *iaopt,
		    int *nia)
{
   /* v2 read header of BYU SIR file --- written by D.Long March 1997 */

/*  Header consists of a variable number of 512 byte blocks.  The
    block contains scaling information and strings:
    note that character strings may not be null terminated in file

    The first block consists of 256 short integers.  Indexing them
    from 1..256 in an array temp, they contain:

	temp(1)	= nsx			! pixels in x direction
	temp(2) = nsy			! pixels in y direction
        temp(3) <= xdeg			! span of deg in x
	temp(4) <= ydeg			! span of deg in y
	temp(5) = nhtype                ! header type (old<15)
	temp(6) <= ascale               ! x scaling
	temp(7) <= bscale               ! y scaling
	temp(8) <= a0                   ! x (lon) origin
         note: longitudes should be in the range -180 to 180
	temp(9) <= b0                   ! y (lat) origin
	temp(10) = ioff			! offset to be added to scaled val
	temp(11) = iscale		! scale factor ival=(val-ioff)/iscale
	temp(12) = iyear		! year for data used
	temp(13) = isday		! starting JD
	temp(14) = ismin		! time of day for first data (in min)
	temp(15) = ieday		! ending JD
	temp(16) = iemin		! time of day for last data (in min)
	temp(17) = iopt			! projection type
					!  -1 = no projection, image only
					!   0 = rectalinear lat/lon
					!   1 = lambert equal area
					!   2 = lambert equal area (local rad)
					!   5 = polar stereographic
					!  11 = EASE north equal area grid
					!  12 = EASE south equal area grid
					!  13 = EASE cylindrical grid
	temp(18) = iregion		! region id code
	temp(19) = itype		! image type code
                                        ! standard values: 0=unknown or n/a
                                        ! 1 = scatterometer A (dB)
                                        ! 2 = scatterometer B (dB/deg)
                                        ! 3 = radiometer Tb (K)
                                        ! 9 = topography (m)
	temp(20)-temp(39) 40 chars of sensor
        temp(40) = 0 
        temp(41) = nhead                ! number of 512 byte header blocks
        temp(42) = ndes                 ! number of 512 byte blocks description
        temp(43) = ldes                 ! number of bytes of description
        temp(44) = nia                  ! number of optional integers
        temp(45) = ipol                 ! polarization (0=n/a,1=H,2=V)
        temp(46) = ifreqhm              ! frequency in 100's of MHz
        temp(47) = ispare1              ! spare
        temp(48) = idatatype            ! data type code 0,2=i*2,1=i*1,4=f

       the value of idata type determines how data is stored and how
       anodata, vmin, and vmax are stored.

        if idatatype = 1 data is stored as bytes and minv=128
        if idatatype = 2 data is stored as 2 byte integers and minv=32767
        if idatatype = 4 data is stored as IEEE floating point

        if idatatype = 1,2 anodata,vmin,vmax are stored as 2 byte integers
         in temp(49)..temp(51)  minv, ioff and iscal used to convert
         integers or bytes into floating point values
         nodata, vmin, and vmax must be representable with ioff and iscale
            temp(*) = (value-ioff)*iscale-minv
            value = float(temp(*)+minv)/float(iscale)+ioff
        idatatype=2 is considered the SIR standard format

        if idatatype = f anodata,vmin,vmax are stored as floating points
         in temp(52)..temp(57) and minv, ioff and iscale are ignored here
         and when reading the file.
         floating point numbers are NOT standard across platforms and
         is not recommended

        temp(49) <= anodata            ! value representing no data
        temp(50) <= vmin               ! minimum useful value from creator prg
        temp(51) <= vmax               ! maximum useful value from creator prg
        temp(52,53) = anodata          ! IEEE floating value of no data
        temp(54,55) = vmin             ! IEEE floating minimum useful value
        temp(56,57) = vmax             ! IEEE floating maximum useful value

	temp(58)-temp(126) 138 chars of type
        temp(127) = 0 
        temp(128) = 0 
	temp(129)-temp(168) 80 chars of title
        temp(169) = 0 
	temp(170)-temp(189) 40 chars of tag
        temp(190) = 0 
	temp(191)-temp(240) 100 chars of crproc
        temp(241) = 0 
	temp(242)-temp(255) 28 chars of crtime
        temp(256) = 0 

     optional header blocks:

	ndes header blocks of 512 bytes: chars of description
	nhead-ndes-1 header blocks of 512 bytes: values of iaopt
        by convention, first value iaopt is a code telling how to interpret
        the rest of the array if nia>0

        remainder of file is image data in a multiple of 512 byte blocks
        two byte integer and byte scaling is
          intval = (fvalue-ioff)*iscale-minv
          fvalue = float(intval+minv)/float(iscale)+ioff

*/

   short in, temp[256];
   int nch,i;
   float f,soff;

   union 
     {
       short i2[2];
       float f2;
   } un;

   rewind(imf);

/* read first header block */

   if (fread(temp, sizeof(short), 256, imf) == 0) return(-1);
#ifdef SWAP
   swap(temp, 256);
#endif

/* decode nsx */

   *nsx = temp[0];
   *nsy = temp[1];
   *xdeg = 0.01 * temp[2];
   *ydeg = 0.01 * temp[3];
   *nhtype = temp[4];
   if (*nhtype < 20) *nhtype=1;
   *ascale = 0.001 * temp[5];
   *bscale = 0.001 * temp[6];
   *a0 = 0.01 * temp[7];
   *b0 = 0.01 * temp[8];
   *ioff = temp[9];
   *iscale = temp[10];
   if (*iscale == 0) *iscale=1;
   *iyear = temp[11];
   *isday = temp[12];
   *ismin = temp[13];
   *ieday = temp[14];
   *iemin = temp[15];
   *iopt  = temp[16];
   *iregion = temp[17];
   *itype = temp[18];

   switch (*iopt){
    case -1:
      *xdeg = *xdeg * 10.0;
      *ydeg = *ydeg * 10.0;
      break;
    case 0:
      *xdeg = *xdeg + 100.0;
      break;
    case 1:
    case 2:
      *ascale = 1./ *ascale;
      *bscale = 1./ *bscale;
      *a0 = *a0 * 100.0;
      *b0 = *b0 * 100.0;
      break;
    case 5:
      *ascale = *ascale * 10.0;
      *bscale = *bscale * 10.0;
      *a0 = *a0 *100.0;
      *b0 = *b0 *100.0;
      *xdeg = *xdeg + 100.0;
      break;
    case 11:
    case 12:
    case 13:
      *ascale = (float) (((double) *ascale) * (double) 2.0 * 6371.228/25.067525);
      *bscale = (float) (((double) *bscale) * (double) 2.0 * 25.067525);
      *xdeg = *xdeg * 10.0;
      *ydeg = *ydeg * 10.0;
      *a0 = *a0 * 10.0;
      *b0 = *b0 * 10.0;
      break;
    default:
      fprintf(stderr,"\n *** Unrecognized SIR option in read_sir_header ***\n");
   }

   *nhead = temp[40];
   if (*nhead == 0) *nhead=1;
   *ndes = temp[41];
   *ldes = temp[42];
   *nia = temp[43];
   *ipol = temp[44];
   *ifreqhm = temp[45];
   *ispare1 = temp[46];
   *idatatype = temp[47];
   if (*idatatype == 0) *idatatype = 2;

   if (*iscale == 0) *iscale=1;
   soff = 32767.0/(float) *iscale;
   if (*idatatype == 1) soff = 128.0/(float) *iscale;

   *anodata = (float) temp[48]/(float) *iscale + *ioff + soff;
   *v_min = (float) temp[49]/(float) *iscale + *ioff + soff;
   *v_max = (float) temp[50]/(float) *iscale + *ioff + soff;

   if (*idatatype == 4) {
#ifdef SWAP
     un.i2[1]=temp[51];
     un.i2[0]=temp[52];
     *anodata = un.f2;
     un.i2[1]=temp[53];
     un.i2[0]=temp[54];
     *v_min = un.f2;
     un.i2[1]=temp[55];
     un.i2[0]=temp[56];
     *v_max = un.f2;
#else
     un.i2[0]=temp[51];
     un.i2[1]=temp[52];
     *anodata = un.f2;
     un.i2[0]=temp[53];
     un.i2[1]=temp[54];
     *v_min = un.f2;
     un.i2[0]=temp[55];
     un.i2[1]=temp[56];
     *v_max = un.f2;
#endif
   }
   
   for (in = 0; in < 20; in++) {
      sensor[2*in] = temp[19+in] % 256;
      sensor[2*in+1] = temp[19+in]/256; 
   }
   sensor[40] = '\0';

   for (in = 0; in < 69; in++) {
      type[2*in] = temp[57+in] % 256;
      type[2*in+1] = temp[57+in]/256;
   }
   type[138] = '\0';

   for (in = 0; in < 40; in++) {
      title[2*in] = temp[128+in] % 256;
      title[2*in+1] = temp[128+in]/256;
   }
   title[80] = '\0';

   for (in = 0; in < 20; in++) {
      tag[2*in] = temp[169+in] % 256;
      tag[2*in+1] = temp[169+in]/256;
   }
   tag[40] = '\0';

   for (in = 0; in < 50; in++) {
      crproc[2*in] = temp[190+in] % 256;
      crproc[2*in+1] = temp[190+in]/256;
   }
   crproc[100] = '\0';

   for (in = 0; in < 14; in++) {
      crtime[2*in] = temp[241+in] % 256;
      crtime[2*in+1] = temp[241+in]/256;
   }
   crtime[28] = '\0';

   if (*nhtype == 1) {  /* old style header */
     *nhead=1;
     *ndes=0;
     *ldes=0;
     *nia=0;
     type[0]='\0';
     tag[0]='\0';
     crproc[0]='\0';
     crtime[0]='\0';
   }

   if (*nhead > 1) {   /* read additional header blocks */
      if (*ndes > 0) {
	 fseek(imf, 512, REL_BEGIN);
	 nch=(maxdes < *ldes ? maxdes : *ldes);
	 nch=0;
	 for (i = 0; i < *ndes; i++) {
	   if (fread(temp, sizeof(short), 256, imf) == 0) return(-1);
#ifdef SWAP
	   swap(temp, 256);
#endif
	   for (in = 0; in < 256; in++) {
	     if (nch < maxdes) descrip[nch] = temp[in] % 256;
	     nch++;
	     if (nch < maxdes) descrip[nch] = temp[in]/256;
	     nch++;
	   }
	 }
      }
      if (maxdes > 0) descrip[maxdes-1] = '\0';

      if (*nhead-*ndes-1 > 0) {
	 fseek(imf, 512*(*ndes+1), REL_BEGIN);
	 nch=(maxi < *nia ? maxi : *nia);
	 if (fread(iaopt, sizeof(short), nch, imf) == 0) return(-1);
#ifdef SWAP
	 swap(iaopt, nch);
#endif
      }
   }

   return(0);
}


void print_head(int nhead, int ndes, int nhtype, int idatatype, int nsx,
		int nsy, float xdeg, float ydeg, float ascale,
		float bscale, float a0, float b0, 
		int ioff, int iscale, int iyear, int isday, int ismin,
		int ieday, int iemin, int iregion, int itype, int iopt,
		int ipol, int ifreqhm, int ispare1,
		float anodata, float v_min, float v_max,
		char *sensor, char *title, char *type, char *tag, char
		*crproc, char *crtime, 
		char *descrip, int ldes, short *iaopt, int nia)
{
  int i;

  if (nhtype < 20) fprintf(stdout,"  (Old style header) %d\n",nhtype);
  fprintf(stdout,"  Title:   '%s'\n",title);
  fprintf(stdout,"  Sensor:  '%s'\n",sensor);
  if (nhtype > 16) {
     fprintf(stdout,"  Type:    '%s'\n",type);
     fprintf(stdout,"  Tag:     '%s'\n",tag);
     fprintf(stdout,"  Creator: '%s'\n",crproc);
     fprintf(stdout,"  Created: '%s'\n",crtime);
  }
  fprintf(stdout,"  Size: %d x %d    Total:%d\n",nsx,nsy, nsx*nsy);
  fprintf(stdout,"  Offset: %d  Scale: %d\n",ioff,iscale);
  fprintf(stdout,"  Year: %d  JD range: %d-%d\n",iyear,isday,ieday);
  fprintf(stdout,"  Region Number: %d  Type: %d  Form: %d\n",iregion,itype,iopt);
  if (nhtype > 16) {
     fprintf(stdout,"  Polarization: %d  Frequency: %f MHz\n",ipol,ifreqhm*0.1);
     fprintf(stdout,"  Datatype: %d  Headers: %d  Ver:%d\n",idatatype,nhead,nhtype);
     fprintf(stdout,"  Nodata: %f\n   Vmin: %f  Vmax: %f\n",anodata,v_min,v_max);
     if (ldes > 0)
	fprintf(stdout,"  Description: (%d) '%s'\n",ldes,descrip);
     if (nia > 0) {
	fprintf(stdout,"  Extra Ints: %d\n",nia);
	for (i=0; i < nia; i++)
	   fprintf(stdout,"     %d %d",i,iaopt[i]);
     }	
  }

  switch(iopt) {
   case -1:
     fprintf(stdout,"  Rectangular image-only: \n");
     fprintf(stdout,"   Xspan,  Yspan:  %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"   Xscale, Yscale: %f , %f\n",ascale,bscale);
     fprintf(stdout,"   Xorg,   Yorg:   %f , %f\n",a0,b0);
     break;

   case 0:
     fprintf(stdout,"  Rectangular Lat/Long form: \n");
     fprintf(stdout,"   Size (deg):     %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"   Lon, Lat scale: %f , %f (pix/deg)\n",ascale,bscale);
     fprintf(stdout,"   Offsets:        %f , %f\n",a0,b0);
     break;

   case 2:
     fprintf(stdout,"  Lambert form: (local radius used)\n");
   case 1:
     if (iopt==1) fprintf(stdout,"  Lambert form:\n");
     fprintf(stdout,"   Center point:      %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"   Lon, Lat scale:    %f , %f (km/pix)\n",1./ascale,1./bscale);
     fprintf(stdout,"   Lower-Left Corner: %f , %f\n",a0,b0);
     break;

   case 5:
     fprintf(stdout,"  Polar sterographic form: \n");
     fprintf(stdout,"   Center Lon,Lat:    %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"   X,Y scales:        %f , %f (km/pix)\n",ascale,bscale);
     fprintf(stdout,"   Lower-Left Corner: %f , %f\n",a0,b0);
     break;

   case 11:
   case 12:
     fprintf(stdout,"  EASE polar azimuthal form: \n");
     fprintf(stdout,"   Map center (col,row): %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"   A,B scales:           %f , %f\n",ascale,bscale);
     fprintf(stdout,"   Map origin (col,row): %f , %f\n",a0,b0);
     break;

   case 13:
     fprintf(stdout,"  EASE cylindrical form: \n");
     fprintf(stdout,"   Map center (col,row): %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"   A,B scales:           %f , %f\n",ascale,bscale);
     fprintf(stdout,"   Map origin (col,row): %f , %f\n",a0,b0);
     break;

   default:
     fprintf(stdout,"  Unrecognized SIR file option: \n");
     fprintf(stdout,"   Xspan,  Yspan:  %f , %f\n",xdeg,ydeg);
     fprintf(stdout,"   Xscale, Yscale: %f , %f\n",ascale,bscale);
     fprintf(stdout,"   Xorg,   Yorg:   %f , %f\n",a0,b0);
     break;
  }
  
  fprintf(stdout,"  Image Min, Max: %f , %f\n\n",v_min,v_max);
  return;
}


void swap(short *i, int n)   /* byte swap short (2 byte) integers */
{
   char *s, t;
   int j;
   
   for (j = 0; j < n; j++) {
     s = (char *) (i+j);
     t = *s;
     *s = *(s+1);
     *(s+1) = t;
   }
   return;
}

#ifndef SWAP
void swapi(int *i, int n)  /* byte swap long (4 byte) integers */
{
   char *s, t, t1;
   int j;
   
   for (j = 0; j < n; j++) {
     s = (char *) (i+j);
     t = *s;
     t1 = *(s+1);
     *s = *(s+3);
     *(s+1) = *(s+2);
     *(s+2) = t1;
     *(s+3) = t;
   }
   return;
}
#endif

int writebmp(char *name_in, char *image, int *nw, int *nh,
	     char *rmap, char *gmap, char *bmap, int *numcols)
{
  FILE *imf;
  int len,i;
  short slen;
  char c;

  imf = fopen(name_in,"w"); 
  if (imf == NULL) {
     return(-1);
  }

  if (fwrite("BM",sizeof(char),2,imf) != 2) {
    return(-1);
    fclose(imf);
  }

  len=256*4+54+ *nh * *nw;  /* length of file */
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=0;
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=54+256*4;    /* offset of data in file */
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=40;    /* length of bitmap header */
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=*nw;
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=*nh;
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }
  
  slen=1;  /* planes */
#ifndef SWAP
  swap(&slen, 1);
#endif
  if (fwrite(&slen,sizeof(short),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  slen=8;  /* bits per pixel */
#ifndef SWAP
  swap(&slen, 1);
#endif
  if (fwrite(&slen,sizeof(short),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }
  
  len=0;  /* compression */
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }
  
  len=*nw * *nh;  /* size of compressed image */
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=2925;  /* horizontal pixels/meter */
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }
  len=2925;  /* horizontal pixels/meter */
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=256;  /* colors used*/
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  len=256;  /* important colors*/
#ifndef SWAP
  swapi(&len, 1);
#endif
  if (fwrite(&len,sizeof(int),1,imf) != 1) {
    return(-1);
    fclose(imf);
  }

  for (i=0; i<256; i++) {  /* color map */
    if (fwrite(&bmap[i],sizeof(char),1,imf) != 1) {
      return(-1);
      fclose(imf);
    }
    if (fwrite(&gmap[i],sizeof(char),1,imf) != 1) {
      return(-1);
      fclose(imf);
    }
    if (fwrite(&rmap[i],sizeof(char),1,imf) != 1) {
      return(-1);
      fclose(imf);
    }
    c='\0';
    if (fwrite(&c,sizeof(char),1,imf) != 1) {
      return(-1);
      fclose(imf);
    }
  }

  /* write image bytes */

  if (fwrite(image,sizeof(char),*nw * *nh,imf) != *nw * *nh) {
    return(-1);
    fclose(imf);
  }
  
  fclose(imf);
  
  return(0);
}