//
//  program to read a BYU SIR image file header and generate SIR images
//  of the lat/lon pixel locations
//
//   Written by DGL 16 Apr 2003 + based on csir_locmap.c
//
//   should be linked with: libsir.cpp sir_ez.c sir_geom.c sir_io.c
//   requires includes: libsir.h sir3.h
//

#include "libsir.h"  // cpp interface routines

#include <iostream.h>

void byte_file(char *fname,float smin,float smax, float *stval, 
	       struct sirheader *head);

int main(int argc, char **argv)
{
  cout << "BYU SIR sirlocmap c++ program";
  if(argc < 2) {
    cout << "\nusage: %s file <ref> \n\n",argv[0];
    cout << " input parameters:\n";
    cout << "   file  = input SIR file\n";
    cout << "   <ref> = (optional) pixel lat,lon reference location\n";
    cout << "            L=lower-left corner (default), C=pixel center\n";
    return 0;
  }

  // read SIR image file into memory

  float *stval;                 // pointer to image storage
  struct sirheader head;        // SIR file header structure

  stval=LoadSIR(argv[1], &head);
  if (stval == NULL) {
    cout << "*** Error reading SIR header from file" << endl;
    return -1;
  }

  // print out contents of sir file header
  print_sir_header(&head);

  /* By default, SIR geometric transformations refer to the lower-left corner
     of the image. However, make it possible to output the center coordinates
     of the pixels.*/

  float half=0.0;
  if (argc > 2)
    if (*(argv[2]) == 'C' || *(argv[2]) == 'c') half=0.5;
  if (half > 0.1)
    printf(" Using pixel center lat,lon coordinates\n");
  else
    printf(" Using standard lower-left pixel corner lat,lon coordinates\n");
 
  /* The geometry routines compute both latitude and longitude at the same time
     However, to avoid having to have to create two large image arrays, we
     compute the geometry twice, once for the latitude and once for the 
     longitude, storing the results in an image array and writing the 
     output file before computing the next image */

  /* first do latitude */
 
  int nsx=head.nsx;
  int nsy=head.nsy;

  int   i, ix, iy;
  float x, y, alon, alat;

  cout << "Computing latitude image\n";
  // note: x,y indexing is 1-based (rather than 0-based) to be
  // consistent with fortran and the documentation
  for (ix = 1; ix <= nsx; ix++) {
    x = ix+half;
    for (iy = 1; iy <= nsy; iy++) {
      i = SIRlex(ix, iy, &head);  // get image pixel index  0 ... nsx*nxy-1
      y = iy+half;
      pixtolatlon(&head, x, y, &alon, &alat);  // (ix,iy) -> (lat,lon)
      stval[i]=alat;
    }
  }
  
  head.ioff=-91;        // change header scaling options to represent
  head.iscale = 360;    // latitude values
  head.anodata = -91.0;
  head.v_min = -90.0;
  head.v_max = 90.0;
  head.itype = 31;     // set image type to latitude
  head.idatatype = 2;  // default two byte integer storage

  sprintf(head.crproc,"BYU sirlocmap");
  sprintf(head.type,"Lat. of %s",argv[1]);

  char outfname[180];
   // generate output file name
  sprintf(outfname,"%s_lat",argv[1]);

  cout << "Writing output SIR file" << endl;  
  if (PutSIR(outfname, &head, stval) < 0) {
    cout << "*** ERROR writing output Latitude file ***" << endl;
    return -1;
  }

  /* now longitude */

  cout << "Computing longitude image\n"; 
  for (ix = 1; ix <= nsx; ix++) {
    x = ix+half;
    for (iy = 1; iy <= nsy; iy++) {
      i = SIRlex(ix, iy, &head);  // get image pixel index  0 ... nsx*nxy-1
      y = iy+half;
      pixtolatlon(&head, x, y, &alon, &alat);  // (ix,iy) -> (lat,lon)
      stval[i]=alon;
    }
  }

  head.ioff=-181;       // change header scaling options to represent
  head.iscale = 180;    // longitude values
  head.anodata = -181.0;
  head.v_min = -180.0;
  head.v_max = 180.0;
  head.itype = 39;     // set image type to longitude
  head.idatatype = 2;  // default two byte integer storage

  sprintf(head.crproc,"BYU sirlocmap");
  sprintf(head.type,"Long. of %s",argv[1]);

   // generate output file name
  sprintf(outfname,"%s_lon",argv[1]);

  cout << "Writing output SIR file" << endl;  
  if (PutSIR(outfname, &head, stval) < 0) {
    cout << "*** ERROR writing output Latitude file ***" << endl;
    return -1;
  }

  return 0;
}