

#include "ast4vl.hpp"
#include "ast_frameset.hpp"
#include "io.hpp"
#include "math.h" // round() needed
#include "cutout_ostream.hpp" // coordinates needed

#include <ostream>
//#include <vector>

using namespace std;


std::ostream& operator << (std::ostream & o, const point2d &a)
{
         return o << "(" << a.lon << ", " << a.lat << ")" ;
}



std::ostream& operator<<( std::ostream &out, struct Bounds const& p)
{
      out << p.label <<" ["<< p.naxis <<"] " << "[" << p.unit  << "] " << p.low_str << " .. " << p.up_str;
         return out;
}



std::ostream& operator<<( std::ostream &out, struct uint_bounds const& p)
{
      out << "(" << p.pix1 << ", " << p.pix2 << ")";
         return out;
}

std::ostream& operator<<( std::ostream &out, overlap_ranges const& p)
{
   out << p.ov_code;
   for(double_xy r : p.pixel_ranges) out  <<" ("<< r.x << ", " << r.y  << ")";

   return out;
}



std::vector<point2d> calc_skyvertices(std::string header, std::string skysys)
{
   LOG_trace(__func__);

   ast::frameset frm_set(header);
   frm_set.set_skysystem(skysys);

   vector<point2d> skyvert = frm_set.sky_vertices();

   return skyvert;
}


std::vector<Bounds> calc_bounds(std::string header, std::string skysys, std::string specsys)
{
   LOG_trace(__func__);

   ast::frameset frm_set(header);
   frm_set.set_skysystem(skysys); 
   if(frm_set.has_specframe())
      frm_set.set_specsystem(specsys);

   vector<Bounds> bounds_vec = frm_set.bounds();

   return bounds_vec;
}



/*
 * ASTlib manual p.604: H.2 Changes Introduced in V1.2
 * ...
 * 6. When a FrameSet is created from a set of FITS header cards (by reading from a FitsChan
 *    using a “foreign” encoding), the base Frame of the resulting FrameSet now has its Domain
 *    attribute set to “GRID”. This reflects the fact that this Frame represents FITS data grid
 *    coordinates (equivalent to FITS pixel coordinates—see §7.13). Previously, this Domain
 *    value was not set.
 * ...
 *
 * p77: 7.13 Conventions for Domain Names
 * ...
 * GRID
 * Identifies the instantaneous data grid used to store and handle data, together
 * with an associated coordinate system. In this coordinate system, the first el-
 * ement stored in an array of data always has a coordinate value of unity at its
 * centre and all elements have unit extent. This applies to all dimensions.
 * ...
 * PIXEL
 * Identifies an array of pixels and an associated pixel-based coordinate system
 * which is related to the GRID coordinate system (above) simply by a shift of
 * origin along each axis. This shift may be integral, fractional, positive, negative
 * or zero. The data elements retain their unit extent along each axis.
 * ...
 * The GRID domain (which corresponds with the pixel-numbering convention
 * used by FITS) is a special case of the PIXEL domain and avoids this uncertainty.
 *
 *
 * FIXME add check : 'Domain ?= GRID' where assumption is made that first data point is centered at (1,1)
 */
std::vector<uint_bounds> calc_overlap(const std::string header, const coordinates coord, int& ov_code)
{
   LOG_trace(__func__);

   ast::frameset frm_set(header);

   LOG_STREAM << "INPUT coord: " << coord << endl;

   frm_set.set(coord.skysys);

   if(frm_set.has_specframe())
      frm_set.set(coord.specsys);

   if(frm_set.has_timeaxis())
      frm_set.set(coord.timesys);

   overlap_ranges pix_ranges = frm_set.overlap(coord);

   ov_code = pix_ranges.ov_code;

   LOG_STREAM << "ov-code & pix ranges[double]: " << pix_ranges << endl;

   /* convert to uint */

   vector<uint_bounds> uint_bounds_vec;

   LOG_STREAM << "pix ranges[uint]:";
   for(double_xy dbl_range : pix_ranges.pixel_ranges)
   {
      if(dbl_range.x < 0)
         throw out_of_range(string{__FILE__} + ":" + to_string(__LINE__)
               + " pixel axis from overlap x is negative " + to_string(dbl_range.x));

      if(dbl_range.y < 0)
         throw out_of_range(string{__FILE__} + ":" + to_string(__LINE__)
               + " pixel axis from overlap y is negative " + to_string(dbl_range.y));

      // FIXME review conversion double -> uint: result must be: 1 <= result <= NAXISn
      // because NAXISn start with 1 (FITS standard) which corresponds to ASTlib's GRID-domain
      // FitsChan uses GRID Domain for FITS-pixel coords
      if(dbl_range.x <= dbl_range.y)
      {
         uint_bounds ui_range{round(dbl_range.x), /*round*/(dbl_range.y)};
         uint_bounds_vec.push_back(ui_range);
         LOG_STREAM << " " << ui_range;
      }
      else
      {
         uint_bounds ui_range{round(dbl_range.y), /*round*/(dbl_range.x)};
         uint_bounds_vec.push_back(ui_range);
         LOG_STREAM << " " << ui_range;
      }

   }
   LOG_STREAM << endl;

   return uint_bounds_vec;
}


