#include <stdlib.h>
#include <cstdlib>
#include <iostream>
#include <iterator>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <list>
#include <random>

#include "VarStr.h"
#include "SpatialVector.h"
#include "SpatialInterface.h"
#include "gener.h"
#include "chealpix.h"

static const double pi=3.141592653589793238462643383279502884197;
static const double twopi=6.283185307179586476925286766559005768394;
static const double halfpi=1.570796326794896619231321691639751442099;
static const double inv_halfpi=0.6366197723675813430755350534900574;


using namespace std;

std::random_device rd;
std::default_random_engine generator(rd()); 

void split_line(string& line, string delim, string& values)
{
    size_t pos = 0;
    while ((pos = line.find(delim, (pos + 1))) != string::npos) {
	values = line.substr(0, pos);
	line = line.substr(pos + 1);
    }
}


void usage(char *name) {

  cout << "usage: " << endl
       << name << " [-depth depth] [-nside ns] [-file f] [-infile f_in] [-rapos rapos] [-decpos decpos]" << endl;
  cout << " [-depth  depth]  : htm depth" << endl      
       << " [-nside ns]      : healpix depth" << endl 
       << " [-file f]        : output file name" << endl
       << " [-infile f_in]     : input file name" << endl
       << " [-rapos rapos]     : position of RA (0 for first column) in file" << endl
       << " [-decpos decpos]   : position of DEC (0 for first column) in file" << endl;
       
  exit(0);
}



int
main(int argc, char *argv[]) {

//*******************************************************
//
// Initialization
//
//*******************************************************

  int args = 1;			// counter
  int arg = 4;			// number of required arguments
  double hp_theta=0.0;
  double hp_phi=0.0;
  uint64 htm_idx;
  hpint64 healpixring;
  hpint64 healpixnest;
  size_t depth=25;
  float64 ra, dec;
  VarStr varg;
  hpint64 nside=8192;
  htmInterface *htm;
  string value;
  string value_end;
  list<string> values;
  char filename[100], filename_in[100];
  int c_in=0, count=0;
  int rapos=0, decpos=1;
  argc--;
  
  if (argc<=0) usage(argv[0]); 
  if (argc>12) usage(argv[0]); 

  while(argc > 0) {
    if(strcmp(argv[args],"-depth")==0) {
      varg = argv[++args];
      if(!htmInterface::isInteger(varg))usage(argv[0]);
      depth = atoi(varg.data());
      argc--; 
    } else if(strcmp(argv[args],"-nside")==0) { 
      varg = argv[++args];
      if(!htmInterface::isInteger(varg))usage(argv[0]);
      nside = atol(varg.data());
      argc--; 
    } else if(strcmp(argv[args],"-file")==0) { 
      varg = argv[++args];
      sprintf(filename,varg.data());
      argc--; 
    } else if(strcmp(argv[args],"-infile")==0) { 
      varg = argv[++args];
      sprintf(filename_in,varg.data());
      argc--; 
    } else if(strcmp(argv[args],"-rapos")==0) { 
      varg = argv[++args];
      rapos= atoi(varg.data());
      argc--; 
    } else if(strcmp(argv[args],"-decpos")==0) { 
      varg = argv[++args];
      decpos= atoi(varg.data());
      argc--; 
    } else {
//      if( *argv[args] == '-' ) usage(argv[0]);
      usage(argv[0]);
      arg--;
    }
    args++;
    argc--;
  }

  htm = new htmInterface(depth);
  
  time_t t0=clock();
  std::ifstream  filein(filename_in);
  std::ofstream  fileout(filename);
  
  while(filein.good())
  {
    value=' ';
    getline(filein, value, '|'); 
    if (value.find('\n') != string::npos) {
        value_end=' ';
        split_line(value, "\n", value_end);
        fileout<<value_end<<"|";
        c_in=0;
        count++;
    } else {
        fileout<<value<<"|";
        c_in++;
    } 
    if (c_in==rapos) {
      ra=strtod(value.c_str(),NULL);
    } else if (c_in==decpos) {
      dec=strtod(value.c_str(),NULL);
    } 
    if ((c_in==0) && (count>0)) {
      htm_idx = htm->lookupID(ra,dec);
      hp_theta=(90.0-dec)/180.0*pi;
      hp_phi=ra/180.0*pi; 
      ang2pix_nest64(nside,hp_theta,hp_phi,&healpixnest);
      ang2pix_ring64(nside,hp_theta,hp_phi,&healpixring);
      fileout<<htm_idx<<"|"<<healpixnest<<"|"<<healpixring<<endl;
      value.erase(std::remove(value.begin(),value.end(),'\n'),value.end()); 
      if ((value!="") && (value!="\n")) { fileout<<value<<"|";} 
   }
  }
  time_t t1=clock();
  double nsec=(double) (t1-t0)/(double)CLOCKS_PER_SEC; 
  cout<<"Number of stars: "<<count<<" processed for "<<nsec<<" sec"<<endl; 
  return 0;
}