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

extern int swap4(int);
extern short swap2(short);
extern float swap4_float(float);
extern float get_float(char *, int);
extern short get_short(char *, int);
extern void utmc(double, double, double *, double *, int *);
extern void utmll(double, double, double *, double *, int);
extern int isbig();

#define nint(a)  ( ( (a) > (0.0) ) ? ((int)(a+0.5)) : ((int)(a-0.5)) )

int main (int argc, char **argv)
{
  FILE *infile, *insfile, *outfile, *roifile;
  char *inbuf, *insbuf, *out, *lineheader;
  char *ptr;
  char *block_ptr;
  int i, j, n, loc, endian;
  int ds_offset, offset, datasize, numrecs, swath, numrecsin;
  int dsr_length, ins_length, mph_size, sph_size;
  int days, secs, musecs, packet_size, num_blocks;
  int samps_per_block, block_id, samp, sampi, sampq, samp_index;
  int data_flag;
  int *zone, zzone;
  int track, frame, orbit;
  int acc, last, delta, outreclen;
  short pri, start_window, length_window, pulse_samples, start_window_min, start_window_max;
  short amb_number[7];
  float fs,freq,time_0, time_swath[7],range_swath[7];
  float pulse_duration[7], bandwidth, chirp_amp[7][4], chirp_pha[7][4];
  float i_lut[256][16], q_lut[256][16];
  float prf, plen, bw, chirpslope;
  float re, im;
  double xpos,ypos,zpos,xvel,yvel,zvel,vel;
  double lat_start, lat_stop, lat_ctr, long_start, long_stop, long_ctr, trackangle;
  double *northing, *easting, framee, framen, *framelat, *framelon, frame_range;
  double h, frame_dist, cosbeta, radius;

  /* set aside 100000 as maximum line length, with line headers */
  inbuf=(char *)malloc(100000);
  out=(char *)malloc(100000);
  lineheader=(char *)malloc(68);
  northing=(double *)malloc(sizeof(double));
  easting=(double *)malloc(sizeof(double));
  framelat=(double *)malloc(sizeof(double));
  framelon=(double *)malloc(sizeof(double));
  zone=(int *)malloc(sizeof(int));

  /* are we on a big endian machine ? */
  endian=isbig();
  if(endian==1)printf("Machine is BIG-endian\n");
  if(endian==0)printf("Machine is LITTLE-endian\n");

  printf("Read Envisat data file\n");

  if(argc < 3){ 
    printf("\nusage: %s Input_L0_file INS_file <Number_of_lines> \n\n",argv[0]);
    exit(-1);
  }
  numrecsin=0;
  if(argc > 3){
  numrecsin=atoi(argv[3]);
  }

  infile = fopen(argv[1],"r");
  if(infile == 0){
    printf("infile %i %s does not exist.\n",infile,argv[1]);
    exit(-1);
  }

  insfile = fopen(argv[2],"r");
  if(infile == 0){
    printf("insfile %i %s does not exist.\n",infile,argv[2]);
    exit(-1);
  }
  fseek(insfile, 0, SEEK_END);
  ins_length=ftell(insfile);
  insbuf=(char *)malloc(ins_length);
  printf("Size of characterization file: %i\n",ins_length);
  rewind(insfile);

  /*  open output data file */
  outfile = fopen("raw.dat","w");

  /* read beginning of file to get Main Product Header (MPH) */
  n=fread(inbuf, sizeof(unsigned char), 10000, infile);
  printf("Bytes read: %i \n",n);

  /* get offset into data start */
  ptr=strstr(inbuf,"DS_OFFSET=");
  ds_offset = atoi(ptr+10);
  printf("Offset to data: %i\n",ds_offset);

  /* data set size, number of records, record size */
  ptr=strstr(inbuf,"DS_SIZE=");
  datasize = atoi(ptr+8);
  printf("Data set size: %i \n",datasize);
  ptr=strstr(inbuf,"NUM_DSR=");
  numrecs=atoi(ptr+8);
  printf("Number of records: %i \n",numrecs);
  ptr=strstr(inbuf,"SWATH=");
  swath=atoi(ptr+9); /* note longer jump for "IS in string */
  printf("Swath number: %i\n",swath);

  /*  ptr=strstr(inbuf,"MPH_SIZE=");
      mph_size=atoi(ptr+9);
  */
  mph_size=1247;
  printf("Main Product Header MPH size: %i \n",mph_size);

  ptr=strstr(inbuf,"SPH_SIZE=");
  sph_size=atoi(ptr+9);
  printf("Special Product Header SPH size: %i \n",sph_size);

  ptr=strstr(inbuf,"X_POSITION=");
  sscanf(ptr+11,"%lf",&xpos);
  ptr=strstr(inbuf,"Y_POSITION=");
  sscanf(ptr+11,"%lf",&ypos);
  ptr=strstr(inbuf,"Z_POSITION=");
  sscanf(ptr+11,"%lf",&zpos);
  ptr=strstr(inbuf,"X_VELOCITY=");
  sscanf(ptr+11,"%lf",&xvel);
  ptr=strstr(inbuf,"Y_VELOCITY=");
  sscanf(ptr+11,"%lf",&yvel);
  ptr=strstr(inbuf,"Z_VELOCITY=");
  sscanf(ptr+11,"%lf",&zvel);
  vel=sqrt(xvel*xvel+yvel*yvel+zvel*zvel);
  printf("Spacecraft position: %lf %lf %lf\n",xpos,ypos,zpos);
  printf("Spacecraft velocity: %lf %lf %lf %lf\n",xvel,yvel,zvel,vel);

  /* location of the swath */
  ptr=strstr(inbuf, "START_LAT=");
  lat_start=atoi(ptr+10);
  lat_start=lat_start*1.e-6;
  ptr=strstr(inbuf, "START_LONG=");
  long_start=atoi(ptr+11);
  long_start=long_start*1.e-6;
  ptr=strstr(inbuf, "STOP_LAT=");
  lat_stop=atoi(ptr+9);
  lat_stop=lat_stop*1.e-6;
  ptr=strstr(inbuf, "STOP_LONG=");
  long_stop=atoi(ptr+10);
  long_stop=long_stop*1.e-6;
  lat_ctr=(lat_start+lat_stop)/2.;
  long_ctr=(long_start+long_stop)/2.;
  printf("Satellite latitude start, stop, center: %lf %lf %lf\n",lat_start,lat_stop,lat_ctr);
  printf("Satellite longitude start, stop, center: %lf %lf %lf\n",long_start,long_stop,long_ctr)
;
  ptr=strstr(inbuf, "SAT_TRACK=");
  sscanf(ptr+11,"%lf",&trackangle);
  printf("Satellite track angle: %lf\n",trackangle);

  /*  What are the orbit, track and frame for characterizing the data set? */
  if((trackangle >= 270.0) || (trackangle <= 90.0)){     /* ascending */ 
    if(lat_ctr >= 0.0)frame=nint(lat_ctr/0.05/18)*18+9;
    if(lat_ctr < 0.0)frame=nint((lat_ctr+360.0)/0.05/18)*18+9;
  }
  if((trackangle < 270.0) && (trackangle > 90.0)){     /* descending */ 
    if(lat_ctr >= 0.0)frame=nint((180.0-lat_ctr)/0.05/18)*18+9;
    if(lat_ctr < 0.0)frame=nint((180.0-lat_ctr)/0.05/18)*18+9;
  }
  ptr=strstr(inbuf,"REL_ORBIT=");
  sscanf(ptr+11,"%d",&track);
  ptr=strstr(inbuf,"ABS_ORBIT=");
  sscanf(ptr+11,"%d",&orbit);
  printf("Scene orbit, track, frame, swath: %d %d %d %d\n",orbit,track,frame,swath);

  /*Some radar setup parameters from the instrument characterization file */

  /* read instrument characterization file parameters */
  n=fread(insbuf, sizeof(unsigned char), ins_length, insfile);
  insbuf=insbuf+1625;  /*  skip mph and sph headers  */

  offset=param_offset(2);
  dsr_length=swap4(get_int(insbuf,offset));
  if(endian==1)dsr_length=(get_int(insbuf,offset));
  printf("dsr_length= %i\n",dsr_length);

  offset=param_offset(3);
  freq=swap4_float(get_float(insbuf,offset));
  if(endian==1)freq=(get_float(insbuf,offset));
  printf("Radar Frequency, wavelength: %f %f \n",freq,299792458./freq);
  
  offset=param_offset(4);
  fs=swap4_float(get_float(insbuf,offset));
  if(endian==1)fs=(get_float(insbuf,offset));
  printf("Sample Frequency fs: %f\n",fs);

  offset=param_offset(54); /* chirp parameters for each swath */
  for(i=0;i<7;i++){
    for(j=0;j<4;j++){
      chirp_amp[i][j]=swap4_float(get_float(insbuf,offset+i*36+4*j));
      chirp_pha[i][j]=swap4_float(get_float(insbuf,offset+i*36+4*j+16));
      if(endian==1)chirp_amp[i][j]=(get_float(insbuf,offset+i*36+4*j));
      if(endian==1)chirp_pha[i][j]=(get_float(insbuf,offset+i*36+4*j+16));
    }
    pulse_duration[i]=swap4_float(get_float(insbuf,offset+i*36+32));
    if(endian==1)pulse_duration[i]=(get_float(insbuf,offset+i*36+32));
  }
  printf("Chirp amp, phase coefficients:\n");
  for(i=0;i<7;i++)printf("%f %f %f %f %f %f %f %f %f\n",chirp_amp[i][0],chirp_amp[i][1],chirp_amp[i][2],chirp_amp[i][3],chirp_pha[i][0],chirp_pha[i][1],chirp_pha[i][2],chirp_pha[i][3],pulse_duration[i]);
  bandwidth=2*pulse_duration[swath-1]*chirp_pha[swath-1][2];
  printf("Chirp length, bandwidth, slope: %12.8f %f %f\n",pulse_duration[swath-1],bandwidth,2*chirp_pha[swath-1][2]);
  printf("Chirp length, bandwidth, slope: %12.8f %f %f\n",pulse_duration[swath-1],bandwidth,bandwidth/pulse_duration[swath-1]);

  /* get the adc lookup tables to undo fbaq quantization */
  offset=param_offset(74);
  for(i=0;i<256;i++){
    for(j=0;j<16;j++){
      n=(i+j*256)*4;
      i_lut[i][j]=swap4_float(get_float(insbuf,offset+n));
      if(endian==1)i_lut[i][j]=(get_float(insbuf,offset+n));
    }
  }
  offset=param_offset(77);
  for(i=0;i<256;i++){
    for(j=0;j<16;j++){
      n=(i+j*256)*4;
      q_lut[i][j]=swap4_float(get_float(insbuf,offset+n));
      if(endian==1)q_lut[i][j]=(get_float(insbuf,offset+n));
    }
  }
  /*
  printf("ADC lookup tables i, q: \n");
  for(i=0;i<256;i++){
  for(j=0;j<16;j++){printf("%.3f ",i_lut[i][j]);}
  printf("\n");
  }
  printf("\n");
  for(i=0;i<256;i++){
  for(j=0;j<16;j++){printf("%.3f ",q_lut[i][j]);}
  printf("\n");
  }
  */

  /*  loop over number of data records */

  /* first, loop over the records looking for data window position changes */

  rewind(infile);
  fread(inbuf, sizeof(unsigned char), mph_size+sph_size, infile);
  if(numrecsin > 0)numrecs=numrecsin;
  printf("Converting %i lines.\n",numrecs);

  for(n=0;n<numrecs;n++){ 
    fread(inbuf, sizeof(unsigned char), 68, infile); /* first 68 bytes header each line */
    days=swap4(get_int(inbuf,0));
    secs=swap4(get_int(inbuf,4));
    musecs=swap4(get_int(inbuf,8));
    if(endian==1)days=(get_int(inbuf,0));
    if(endian==1)secs=(get_int(inbuf,4));
    if(endian==1)musecs=(get_int(inbuf,8));
    if(n == 0)printf("Days since 1/1/2000, seconds of day, musecs: %i %i %i\n",days,secs,musecs);
    packet_size=swap2(get_short(inbuf,36));  /* get data packet size */
    if(endian==1)packet_size=(get_short(inbuf,36));  /* get data packet size */
    /* if(n%100==0)printf("%i Packet size %i bytes %i %i %i\n",n,packet_size,days,secs,musecs); */
    pri=swap2(get_short(inbuf,16+38));  /*  pri is repetition interval in clock cycles */
    start_window=swap2(get_short(inbuf,18+38)); /* window start in clocks */
    length_window=swap2(get_short(inbuf,20+38)); /* window length in clocks */
    if(endian==1)pri=(get_short(inbuf,16+38));  /*  pri is repetition interval in clock cycles */
    if(endian==1)start_window=(get_short(inbuf,18+38)); /* window start in clocks */
    if(endian==1)length_window=(get_short(inbuf,20+38)); /* window length in clocks */
    pulse_samples=(((unsigned short)inbuf[24+38]&0xff)<<2) | (((unsigned short)inbuf[25+38]&0xa0)>>6); /* samples in tx pulse */
    if(n==0)printf("Packet size, winstart, winlen, tx samples : %i %i %i %i\n",packet_size, start_window, length_window, pulse_samples);
    /*  if(n%1==0)printf("%i Packet size %i bytes %i %i %i %i %i\n",n,packet_size,days,pri,pulse_samples,start_window,length_window); */

    /* is record an echo data set record? */
    data_flag=((unsigned char)inbuf[38+14] & 0xf0) >> 4;    
    if((data_flag&0x08) > 0){

    /* convert to human units */
    prf=fs/pri;
    plen=pulse_samples/fs;
    bw=((unsigned char)inbuf[38+26])*1.6e7/255.;  /* chirp bandwidth */
    if(n==0)num_blocks=(packet_size-30+1)/64;  /* number of blocks to decode for fbaq */
    /*check for partial block*/
    if(64*num_blocks < (packet_size-30+1))num_blocks++;
    chirpslope=bw/plen;
    if(n==0)printf("Prf, plen, bw, chirp slope, num_blocks %f %12.8f %f %f %i\n",prf,plen,bw,chirpslope,num_blocks);
    if(n==0)start_window_min=start_window;
    if(n==0)start_window_max=start_window;
    if(start_window < start_window_min)start_window_min=start_window;
    if(start_window > start_window_max)start_window_max=start_window;
    }

    /*    if(n%1000==0)printf("Processed %i lines.\n",n); */
    fread(inbuf, sizeof(unsigned char), packet_size-30+1, infile); /* read in data */
  }
  printf("Window start minimum, maximum: %i %i\n",start_window_min,start_window_max);
  time_0=start_window_min/fs;
  printf("Ambiguous Delay to first sample, s: %f\n",time_0);

  /* range to first point in each swath */
  offset=param_offset(106);
  for(i=0;i<7;i++){
    amb_number[i]=swap2(get_short(insbuf,offset+28+i*2));
    if(endian==1)amb_number[i]=(get_short(insbuf,offset+28+i*2));
  }
  printf("Ambiguity number for each swath: %i %i %i %i %i %i %i\n",amb_number[0],amb_number[1],amb_number[2],amb_number[3],amb_number[4],amb_number[5],amb_number[6]);
  for(i=0;i<7;i++){
    time_swath[i]=amb_number[i]/prf+time_0;
    range_swath[i]=299799458.*time_swath[i]/2.;
    printf("Time, Range to each swath: %f %f\n",time_swath[i],range_swath[i]);
  }
  printf("Time, Range to selected swath: %f %f\n",time_swath[swath-1],range_swath[swath-1]);

    /*where is the frame? */

    utmc(lat_ctr, long_ctr, easting, northing, zone);
    h=sqrt(xpos*xpos+ypos*ypos+zpos*zpos);
    radius=6378000.;  /*  get actual earth radius */
    frame_range=range_swath[swath-1]+2500/fs*299799458./2;
    cosbeta=(radius*radius+h*h-pow(frame_range,2))/2/radius/h;
    frame_dist=radius*sqrt(1.-cosbeta*cosbeta);
    framee=*easting+frame_dist*cos(trackangle*3.14159265/180.);
    framen=*northing-frame_dist*sin(trackangle*3.14159265/180.);
    zzone=*zone;
    utmll(framee, framen, framelat, framelon, zzone); 
    printf("Frame latitude, longitude %f %f\n",*framelat, *framelon);
  

    /* now actually read the data and decode the fbaq */

    rewind(infile);
    fread(inbuf, sizeof(unsigned char), mph_size+sph_size, infile);

    last=0;
    for(n=0;n<numrecs;n++){ 
      fread(inbuf, sizeof(unsigned char), 68, infile); /* first 68 bytes header each line */
      days=swap4(get_int(inbuf,0));
      secs=swap4(get_int(inbuf,4));
      musecs=swap4(get_int(inbuf,8));
      packet_size=swap2(get_short(inbuf,36));  /* get data packet size */
      pri=swap2(get_short(inbuf,16+38));  /*  pri is repetition interval in clock cycles */
      start_window=swap2(get_short(inbuf,18+38)); /* window start in clocks */
      length_window=swap2(get_short(inbuf,20+38)); /* window length in clocks */
      if(endian==1){
      days=(get_int(inbuf,0));
      secs=(get_int(inbuf,4));
      musecs=(get_int(inbuf,8));
      packet_size=(get_short(inbuf,36));  /* get data packet size */
      pri=(get_short(inbuf,16+38));  /*  pri is repetition interval in clock cycles */
      start_window=(get_short(inbuf,18+38)); /* window start in clocks */
      length_window=(get_short(inbuf,20+38)); /* window length in clocks */
      }

      pulse_samples=((unsigned short)inbuf[24+38]<<2) | (((unsigned short)inbuf[25+38]&0xa0)>>6); /* samples in tx pulse */
      data_flag=((unsigned char)inbuf[38+14] & 0xf0) >> 4;    

      memcpy(lineheader, inbuf, 68);  /*  save header for output  */
      fread(inbuf, sizeof(unsigned char), packet_size-30+1, infile); /* read in data */
      /* is record an echo data set record? */
      if((data_flag&0x08) > 0){

      /* decode a line counter in the file */
      acc=0;
      for(j = 0;j<4;j++){
       acc = 256*acc + (unsigned char)lineheader[47+j];
      }
      delta=acc-last;
      if(delta != 1)printf("Line counter, last: %d %d %d\n",acc,last,delta);
      last=acc;


	for(i=0;i<num_blocks;i++){     /* begin fbaq decoding */
	  block_ptr=inbuf+i*64;
	  samps_per_block=63;  /*  default value */
	  if(i==num_blocks-1)samps_per_block=(packet_size-30+1)-num_blocks*64; /* last block */
	  block_id=block_ptr[0]; /* block id code */
	  /*      printf("%i\n",block_id);  */
	  /* watch for bad block pointers, seem to occur in corrupted lines */
	  if(block_id <= 0)block_id=0;
	  if(block_id >= 255)block_id=255;
	  /*      if(n>895)printf("%i %i\n",block_ptr,block_id);  */
	  for(j=0;j<samps_per_block;j++){
	    samp=i*63+j;
	    sampi=(block_ptr[j+1] & 0xf0)>>4;
	    sampq=(block_ptr[j+1] & 0x0f);
	    
	    if(sampi & 0x8) sampi=15-sampi;
	    else sampi = (sampi | 0x08);
	    if(sampq & 0x8) sampq=15-sampq;
	    else sampq = (sampq | 0x08);
	    
	    re=i_lut[block_id][sampi];
	    im=q_lut[block_id][sampq];
	
	    samp_index = samp;
	    out[samp_index*2]=nint(re*127.5+127.5);
	    out[samp_index*2+1]=nint(im*127.5+127.5);
	  }
	}
	fwrite(lineheader, sizeof(unsigned char), 68, outfile);
	fwrite(out, sizeof(unsigned char),(packet_size-30+1)*2,outfile);
	outreclen=(packet_size-30+1)*2+68;
	if(n==0)printf("Written lines of %i bytes, first 68 are header\n",(packet_size-30+1)*2+68); 
      }
      else{
	fwrite(lineheader, sizeof(unsigned char), 68, outfile);
	for(i=0; i<outreclen-68; i++){out[i]=0;};
	fwrite(out, sizeof(unsigned char),outreclen-68,outfile);
	printf("Written line of zeros %d, length %d\n",acc,outreclen);
      }
    }
    /*  create a roi input file  */
  roifile = fopen("roi.in","w");
  fprintf(roifile,"raw.dat\n");
  fprintf(roifile,"raw.dat\n");
  fprintf(roifile,"env.int\n");
  fprintf(roifile,"env.amp\n");
  fprintf(roifile,"8lk\n");
  fprintf(roifile,"0                               ! debug flag 1=trans1,2=trans2,4=rcrm,8=rcaz\n");
  fprintf(roifile,"%i                           ! bytes per line\n",(packet_size-30+1)*2+68);
  fprintf(roifile,"10960                           ! good bytes per line\n");
  fprintf(roifile,"1                               ! first line (from 0)\n");
  fprintf(roifile,"1                               ! number of patches 12\n");
  fprintf(roifile,"68                               ! i,q byte samples to skip (from 0)\n");
  fprintf(roifile,"3000                            ! valid samples per patch 2870\n");
  fprintf(roifile,"0                               ! deskew flag\n");
  fprintf(roifile,"0. 0.                           ! Squint angles 1,2\n");
  fprintf(roifile,"1   5120                        ! first range bin, number of bins in file 1\n");
  fprintf(roifile,"0 0                             ! azimuth, range delta for file 2\n");
  fprintf(roifile,"-0. 0. 0.                        ! quadratic Doppler coefficients (file 1)\n");
  fprintf(roifile,"-0. 0. 0.                        ! quadratic Doppler coefficients (file 2)\n");
  fprintf(roifile,"1                               ! 1 = use file 1 doppler, 2=file 2, 3 = avg\n");
  fprintf(roifile,"6378000.                        ! Earth radius (m)\n");
  fprintf(roifile,"%f %f                 ! Body fixed S/C velocity (m/s)\n",vel,vel);
  fprintf(roifile,"781000. 781000.                 ! ht 1,2\n");
  fprintf(roifile,"%f %f       ! r0 1,2\n",range_swath[swath-1],range_swath[swath-1]);
  fprintf(roifile,"%f %f            ! prf 1,2\n",prf,prf);
  fprintf(roifile,"127.5 127.5 127.5 127.5         ! i,q bias 1,2  15.31 15.48 15.25 15.29\n");
  fprintf(roifile,"n                               ! Flip I/Q ?\n");
  fprintf(roifile,"6.                              ! single look azimuth resolution\n");
  fprintf(roifile,"5                               ! azimuth looks 1 or >= 4\n");
  fprintf(roifile,"%f                       ! range sampling frequency (Hz)\n",fs);
  fprintf(roifile,"%f                      ! chirp slope (Hz/s) 6.98426888d11\n",chirpslope);
  fprintf(roifile,"%12.8f                        ! pulse length (s)\n",plen);
  fprintf(roifile,"0                               ! chirp extension\n");
  fprintf(roifile,"n                               ! secondary range compression?\n");
  fprintf(roifile,"%12.8f                       ! wavelength (m)\n",299799458./freq);
  fprintf(roifile,"1.                              ! range spectrum weight (1=none, 0.54=hamming)\n");
  fprintf(roifile,"0. 0.                           ! bandwidth fractional truncation (ra, az)\n");
  fprintf(roifile,"0 0 0 0\n");
  fprintf(roifile,"0 0 0 0\n");
}



int param_offset(int param){ 
/* offset parameter "param" in INS file */
/* Note: ESA input file is in BIG ENDIAN byte order */
/* table of sizes of elements in the Instrument Characterization File (ASA_INS_AX_GADS), see section 7.3 of the ASAR product handbook */

 int i, offset;
 int sz[] = {
 12,4,4,4,4,256,256,256,256,256,
 256,256,256,1792,1792,1792,1792,1792,1792,1792,
 1792,1792,1792,1792,1792,1792,1792,1792,1792,1792,
 1792,1792,1792,1792,1792,1792,1792,1280,1280,1280,
 1280,1280,1280,1280,1280,1280,1280,1280,1280,1280,
 1280,1280,1280,252,252,252,180,180,404,404,
 404,404,404,404,404,404,4,4,1020,1020,
 648,1024,1024,16384,8192,4096,16384,8192,4096,16384,
 8192,4096,64,64,32,32,32,32,32,84,
 84,84,84,84,2,2,28,28,28,28,
 28,2,2,2,64,56,56,56,56,56,
 2,44,4,4,4,4,4,4,4,4,
 64,512,512,512,512,512,512,512,512,5120,
 4,64,4,64,4,64,4,64,4,64,
 4,72};

 offset=0;
 for(i=0;i<param-1;i++)offset += sz[i];
 return(offset);
}

float get_float(char *buffer, int offset){
  float value;
  memcpy(&value, buffer+offset, 4);
  return(value);
}
int get_int(char *buffer, int offset){
  int value;
  memcpy(&value, buffer+offset, 4);
  return(value);
}
short get_short(char *buffer, int offset){
  short value;
  memcpy(&value, buffer+offset, 2);
  return(value);
}
int swap4(int in){
  unsigned char t[4],tt;
  int value;

  /*  printf("Input to swap: %i\n",in); */
  memcpy(t,&in,4);
  /*  printf("swap in %i %i %i %i\n",t[0],t[1],t[2],t[3]); */
  tt=t[0];
  t[0]=t[3];
  t[3]=tt;
  tt=t[1];
  t[1]=t[2];
  t[2]=tt;
  memcpy(&value,t,4);
  return(value);
}
short swap2(short in){
  unsigned char t[2],tt;
  short value;

  memcpy(t,&in,2);
  tt=t[0];
  t[0]=t[1];
  t[1]=tt;
  memcpy(&value,t,2);
  return(value);
}
float swap4_float(float in){
  unsigned char t[4],tt;
  float value;

  /*  printf("Input to swap float: %f\n",in); */
  memcpy(t,&in,4);
  /*  printf("swap in %i %i %i %i\n",t[0],t[1],t[2],t[3]); */
  tt=t[0];
  t[0]=t[3];
  t[3]=tt;
  tt=t[1];
  t[1]=t[2];
  t[2]=tt;
  memcpy(&value,t,4);
  return(value);
}

/*
c convert lat/lon to universal transverse mercator(utm)
c coordinates (easting, northing)
c
c	NOTE: modified UTM ! Delta y is always 1e+7 also on the
c	northern hemisphere. Therefore the northing is 1e+7 larger
c	than the true northing for latitudes larger than 0.
c	The advantage of this approach is that inversion is possible
c	without additional information on whether the location is on the
c	southern or northern hemisphere.
c 
c input parameters:
c       tlat- array of latitude coordinates
c       tlon- array of longitude coordinates
c     numpts- number of points to convert
c    eoffset- cross track error \ of the input image frame, use in
c    noffset- along track error / terrain correction processing only
c          a- semi-major (equatorial) axis of the earth (m)
c        esq- eccectricity square
c         k0- central scale factor
c
c output parameters:
c         tx- array of easting coordinates
c         ty- array of northing coordinates
c
c variables:
c        esq- eccentricity square
c       epsq- e-prime squared
c          a- semi major axis
c         k0- center scale factor
c          n- zone
c      long0- central meridian
c        dtr- degree to radian
c       xadj- easting adjustment
c       yadj- northing adjustment
c
c Reference: Synder, J. P. 1983, Map projection used by the US
c            geological survey(bulletin 1532)
*/

void utmc(double tlat, double tlon, double *tx, double *ty, int *zone){

      int    n, i;
      double epsq, a, lat, lon, long0, a1, a2, a3, rn;
      double t, b1, b2, rm, rm0, dtr, eoffset, noffset;
      double tanlat, c;
      double yadj, xadj;
      double esq, k0;
      double e4, e6, c1, c2, c3, c4;

      eoffset=0;
      noffset=0;
      a=6378137.0;
      esq=0.00669438000;
      k0=0.9996;

      dtr = 3.14159265358979323846264/180.0;
      epsq = esq/(1.0 - esq);

      /* calculate the zone no. and long0 */

      n = (180 + tlon)/6+1;
      *zone = n;
      long0 = (n-1)*6.0 + 3.0 - 180.0;
      lat = tlat * dtr;
      lon = tlon * dtr;
      rn = a/sqrt(1.0 - esq*sin(lat)*sin(lat));
      tanlat = sin(lat)/cos(lat);
      t = tanlat * tanlat;
      c = epsq * cos(lat)*cos(lat);

      a1 = cos(lat) * (tlon - long0) * dtr;
      a2 = (1.0 - t + c) * a1*a1*a1 / 6.0;
      a3 = (5.0 - 18.0*t + t*t + 72.0*c - 58.0*epsq) * pow(a1,5) / 120.0;

      xadj = 500000.0;
      *tx = k0 * rn * (a1 + a2 + a3) + xadj + eoffset;

      e4 = esq * esq;
      e6 = e4 * esq;
      c1 = 1.0 - esq/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0;
      c2 = 3.0*esq/8.0 + 3.0*e4/32.0 + 45.0*e6/1024.0;
      c3 = 15.0*e4/256.0 + 45.0*e6/1024.0;
      c4 = 35.0*e6/3072.0;
      rm = a*(c1*lat - c2*sin(2.0*lat) + c3*sin(4.0*lat) - c4*sin(6.0*lat));
      rm0 = 0.0;
      b1 = (a1*a1)/2.0 + (5.0 - t + 9.0*c + 4.0 * c*c) * pow(a1,4) / 24.0;
      b2 = (61.0 - 58.0*t + t*t + 600.0*c + 330.0*epsq) * pow(a1,6) / 720.0;
      yadj = 1.7;
      *ty = k0 * (rm - rm0 + rn*tanlat*(b1 + b2)) + yadj + noffset;
}

/*
c convert the utm coordinates (northing & easting) of
c a target to lat/lon
c
c input parameters:
c          x- array of easting coordinates (in meters)
c          y- array of northing coordinates (in meters)
c        esq- eccentricity squared
c          a- semi-major (equatorial) axis of earth (in meters)
c         k0- center scale factor
c     numpts- number of points to convert
c       zone- zone number
c
c output parameters:
c        lat- array of latitude coordinates (in degrees)
c       long- array of longitude coordinates (in degrees)
c
c variables:
c          n- zone no.
c        dtr- degree to radian
c       xadj- easting adjustment
c       yadj- northing adjustment
c
c Reference: Synder, J. P. 1983, Map projection used by the US
c            geological survey(bulletin 1532)
c
c Added esq, a, and k0 as input parameters.  Added variable
c rm0.  Redefined PI to have double-precision accuracy.  Added
c the variables xadj and yadj.  Made x, y, lat, and long arrays
c and added the input parameter numpts.
c                                                      jgs 4/90
*/

void utmll(double x, double y, double *lat, double *lon, int zone){

  int i;
  double esq, a, k0, u1, u2, u3, lat1, esqsin2, lat1d, long0, dtr;
  double rm, e1, u, epsq, t1, c1, rn1, r1, rm0;
  double xadj, yadj, tanlat1, d;

  a=6378137.0;
  esq=0.00669438000;
  k0=0.9996;

  dtr = 3.14159265358979323846264/180.0;
  rm0 = 0.0;
  yadj = 1.7;
  long0 = (zone)*6.0 - 183.0;

  rm = (y - yadj)/k0 + rm0;
  e1 = (1.0 - sqrt(1.0 - esq))/(1.0 + sqrt(1.0 - esq));
  u = rm/(a * (1.0 - esq/4.0 - (3.0 * esq * esq/64.0) - (5.0 * esq * esq * esq/256.0) ) );
  u1 = (3.0 * e1 / 2.0 - (27.0 * pow(e1,3))/32.0) * sin(2.0*u);
  u2 = (21.0 * e1*e1/16.0 - (55.0 * pow(e1,4))/32.0) * sin(4.0*u);
  u3 = (151.0 * e1*e1*e1 / 96.0) * sin(6.0*u);
  lat1 = u + u1 + u2 + u3;
  lat1d = lat1/dtr;

  esqsin2 = 1.0 - esq*pow(sin(lat1),2);
  epsq = esq/(1.0 - esq);
  c1 = epsq * pow(cos(lat1),2);
  tanlat1 = sin(lat1)/cos(lat1);
  t1 = tanlat1 * tanlat1;
  rn1 = a/sqrt(esqsin2);
  r1 = a*(1.0 - esq)/sqrt(esqsin2 * esqsin2 * esqsin2);
  xadj = 500000.0;
  d = (x - xadj)/(rn1 * k0);

  *lat = lat1d - ((rn1 * tanlat1/r1) * (d*d*0.50 - (5.0 + 3.0*t1 - 10.0*c1 + 4.0*c1*c1 - 9.0*epsq) * pow(d,4)/24.0 + (61.0 + 90.0*t1 + 298.0*c1 + 45.0*t1*t1 - 252.0*epsq - 3.0*c1*c1) *pow(d,6)/720.0) )/dtr;
  *lon = long0 + ((1.0/cos(lat1)) * (d - (1.0 + 2.0*t1 + c1) * pow(d,3)/6.0 + (5.0 - 2.0*c1 + 28.0*t1 - 3.0*c1*c1 + 8.0*epsq + 24.0*t1*t1) *pow(d,5)/120.0) )/dtr;
}

int isbig(){

  int bigendian, littleendian, test;
  unsigned char t[4];

  littleendian=256;
  bigendian=256*256;

  t[0]=0;
  t[1]=1;
  t[2]=0;
  t[3]=0;

  memcpy(&test, &t[0], 4);

  /* printf("test: %i\n",test); */
  if(test==bigendian)return(1);
  if(test==littleendian)return(0);
  printf("Error in endian test, test= %i ********\n",test);
}