/*
   xmunipack - coordinates

   Copyright © 1997-2013, 2017, 2019-2022 F.Hroch (hroch@physics.muni.cz)

   This file is part of Munipack.

   Munipack is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   Munipack is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with Munipack.  If not, see <http://www.gnu.org/licenses/>.

*/

#include "fits.h"
#include <wx/wx.h>
#include <math.h>

#ifdef __WXDEBUG__
#include <wx/debug.h>
#endif

#define RAD 57.295779513082322865


// -----------------------------------------------------------------

FitsProjection::FitsProjection():
  acen(0.0),dcen(0.0),xcen(0.0),ycen(0.0),scale(1.0),angle(0.0),reflex(1.0) {}


FitsProjection::FitsProjection(const wxString& t, double a,double d,double x,
			       double y,double s,double r,double z):
  type(t.Upper()),acen(a/RAD),dcen(d/RAD),xcen(x),ycen(y),scale(s*RAD),
  angle(r/RAD), reflex(z) {}

FitsProjection::FitsProjection(const wxString& t, double a,double d, double x,
			       double y, double cd11, double cd12, double cd21,
			       double cd22):
  type(t.Upper()),acen(a/RAD),dcen(d/RAD),xcen(x),ycen(y),scale(1.0/RAD),
  angle(0.0), reflex(1.0)
{
  reflex = -cd11*cd22 > 0 ? 1.0 : -1.0;

  if( fabs(cd12) > 0 || fabs(cd22) > 0 )
    angle = atan2(reflex*cd12,cd22);
  else
    type = "";

  double q = fabs(cd11*cd22) + fabs(cd12*cd21);
  if( q > 0 )
    scale = RAD/sqrt(q);
  else
    type = "";

}

bool FitsProjection::IsOk() const { return type != "";  }

double FitsProjection::GetScale() const { return scale; }

void FitsProjection::ad2xy(double a,double d,double& x,double& y) const
{
  if( ! IsOk() ) {
    x = 0;
    y = 0;
    return;
  }

  double aa = cos(angle)*scale;
  double bb = sin(angle)*scale;

  double u,v;

  gnomon(a/RAD,d/RAD,u,v);
  x = xcen +  (aa*u + bb*v)*reflex;
  y = ycen + (-bb*u + aa*v);
}

void FitsProjection::xy2ad(double x,double y,double& a,double& d) const
{
  if( ! IsOk() ) {
    a = 0;
    d = 0;
    return;
  }

  double aa = cos(angle)/scale;
  double bb = sin(angle)/scale;

  double u,v,xx,yy;

  u = x - xcen;
  v = y - ycen;
  xx = (aa*u - bb*v)*reflex;
  yy =  bb*u + aa*v;

  ignomon(xx,yy,a,d);
  a = RAD*a;
  d = RAD*d;
}


void FitsProjection::gnomon(double a,double d,double& x,double& y) const
{
  double c,p,q,r,v,w,s;

  c = cos(d);
  p = sin(d);
  q = c*sin(a - acen);
  r = c*cos(a - acen);
  v = sin(dcen);
  w = cos(dcen);
  s = p*v + r*w;

  x = -q/s;
  y = (w*p - v*r)/s;
}

void FitsProjection::ignomon(double x,double y,double& a,double& d) const
{
  double p,q,r,v,w,t;

  v = sin(dcen);
  w = cos(dcen);
  t = sqrt(1.0 + x*x + y*y);
  p = (v + w*y)/t;
  q = -x/t;
  r = (w - v*y)/t;

  d = asin(p);
  a = atan2(q,r) + acen;
}

double FitsProjection::xy2dist(double x1,double y1,double x2,double y2) const
{
  // https://en.wikipedia.org/wiki/Haversine_formula
  double r1,d1,r2,d2;
  xy2ad(x1,y1,r1,d1);
  xy2ad(x2,y2,r2,d2);

  double sind = sin((d2 - d1)/2/RAD);
  double sinr = sin((r2 - r1)/2/RAD);
  double h = sind*sind + cos(d1/RAD)*cos(d2/RAD)*sinr*sinr;
  return h > 0 ? 2*RAD*asin(sqrt(h)) : 0;
}

// -----------------------------------------------------------------


FitsCoo::FitsCoo(): FitsArray(),type(COO_PIXEL),haswcs(false),digits(4) {}

FitsCoo::FitsCoo(const FitsArray& a):
  FitsArray(a),type(COO_PIXEL),haswcs(false),digits(4)
{
  // decode WCS

  if( a.GetKey("CTYPE1").Find("TAN") != wxNOT_FOUND &&
      a.GetKey("CTYPE2").Find("TAN") != wxNOT_FOUND ) {

    double crpix1 = a.GetKeyDouble("CRPIX1");
    double crpix2 = a.GetKeyDouble("CRPIX2");
    double crval1 = a.GetKeyDouble("CRVAL1");
    double crval2 = a.GetKeyDouble("CRVAL2");
    double cd1_1 = a.GetKeyDouble("CD1_1");
    double cd1_2 = a.GetKeyDouble("CD1_2");
    double cd2_1 = a.GetKeyDouble("CD2_1");
    double cd2_2 = a.GetKeyDouble("CD2_2");
    proj = FitsProjection("GNONOMIC",crval1,crval2,crpix1,crpix2,
			  cd1_1,cd1_2,cd2_1,cd2_2);

    if( proj.IsOk() )
      // significant digits for the print format are determined by scale
      digits = wxMax(int(log10(proj.GetScale()/RAD) + 1),0);

    haswcs = proj.IsOk();
  }

}

bool FitsCoo::HasWCS() const
{
  return haswcs;
}

void FitsCoo::SetType(int t)
{
  type = static_cast<coords_type>(t);
}

void FitsCoo::SetType(const wxString& a)
{
  for(int i = COO_FIRST+1; i < COO_LAST; i++)
    if( a == Label_str(i) ) {
      SetType(i);
      return;
    }
}

coords_type FitsCoo::GetType() const
{
  return type;
}

void FitsCoo::GetEq(int xcen, int ycen, double& alpha, double& delta) const
{
  proj.xy2ad(xcen,ycen,alpha,delta);
}


void FitsCoo::RaSix(double alpha, int& h, int& min, double& sec) const
{
  double x = alpha/15.0;
  h = int(x);
  min = int(60.0*(x - h));
  sec = 3600.0*(x - (h + min/60.0));
}

void FitsCoo::DecSix(double delta, char& sign, int& d, int& m, double& s) const
{
  sign = delta > 0.0 ? '+' : '-';

  double y = fabs(delta);
  d = int(y);
  m = int(60.0*(y - d));
  s = 3600.0*(y - (d + m/60.0));
}


void FitsCoo::GetPix(int i, int j, wxString& a, wxString& b) const
{
  if( 0 <= i && i < Width() && 0 <= j && j < Height() ) {
    a.Printf("%d",i+1); // C-style indexing should be considered very funny!
    b.Printf("%d",j+1);
  }
  else {
    a.Clear();
    b.Clear();
  }
}

void FitsCoo::GetCoo(int i, int j, wxString& a, wxString& b) const
{
  if( haswcs && 0 <= i && i < Width() && 0 <= j && j < Height() ) {

    int h,min,d,m;
    double x,y, sec,s;
    char sg;

    GetEq(i,j,x,y);

    if( type == COO_EQDEG ) {
      a.Printf(L"%.*f°",digits,x);
      b.Printf(L"%.*f°",digits,y);
    }

    else if( type == COO_EQSIX ) {
      int n = wxMax(digits-4,0);
      RaSix(x,h,min,sec);
      DecSix(y,sg,d,m,s);
      a.Printf(L"%02dh%02dm%04.*fs",h,min,n+1,sec);
      b.Printf(L"%c%02d°%02d'%02.*f\"",sg,d,m,n,s);
    }
  }
  else {
    a.Clear();
    b.Clear();
  }
}

wxString FitsCoo::Label_str(int n)
{
  switch(n){
  case COO_PIXEL: return  "pixels";
  case COO_EQDEG: return L"Equatorial (°)";
  case COO_EQSIX: return  "Equatorial (h:m:s)";
  default:        return  "";
  }
}

wxArrayString FitsCoo::Label_str()
{
  wxArrayString a;

  for(int i = COO_FIRST+1; i < COO_LAST; i++)
    a.Add(Label_str(i));

  return a;
}