/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz, 03/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2002-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MAstroCatalog
//
//  THIS IMPLEMENTATION IS PRELIMINARY AND WILL BE MERGED WITH
//  SOME PARTS OF THE DRIVE SOFTWARE SOON!
//
//////////////////////////////////////////////////////////////////////////////
#include "MAstroCatalog.h"

#include <fstream>
#include <stdlib.h>

#include <KeySymbols.h>

#include <TPad.h> // TPad::GetMaxPickDistance
#include <TLine.h>
#include <TMarker.h>
#include <TCanvas.h>
#include <TArrayI.h>
#include <TGToolTip.h>
#include <TRotation.h>
#include <TStopwatch.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MTime.h"
#include "MAstro.h"
#include "MObservatory.h"

ClassImp(MVector3);
ClassImp(MAstroCatalog);

using namespace std;

MVector3 MVector3::GetZdAz(const MObservatory &obs, Double_t gmst) const
{
    if (!fType==kIsRaDec)
        return MVector3();

    const Double_t alpha = gmst + obs.GetElong();

    MVector3 zdaz;
    zdaz.SetZdAz(Theta(), alpha-Phi(), Mag());
    zdaz.RotateY(obs.GetPhi()-TMath::Pi()/2);

    return zdaz;

    /*
     // ------ The same using slalib, tested in the drive system -------
     const Double_t alpha = slaGmst(mjd) + obs.GetElong();
     Double_t el;
     slaDe2h(fAlpha-ra, dec, obs.GetPhi(), &az, &el);
     zd = TMath::Pi()/2-el;
     return;
     */
}

MVector3 MVector3::GetZdAz(const MTime &time, MObservatory &obs) const
{
    return GetZdAz(obs, time.GetGmst());
}

MVector3 MVector3::GetRaDec(const MObservatory &obs, Double_t gmst) const
{
    if (!fType==kIsZdAz)
        return MVector3();

    const Double_t alpha = gmst + obs.GetElong();

    MVector3 v(*this);
    v.RotateY(TMath::Pi()/2-obs.GetPhi());

    MVector3 rd;
    rd.SetRaDec(alpha-v.Phi(), TMath::Pi()/2-v.Theta(), Mag());
    return rd;

    /*
     // ------ The same using slalib, tested in the drive system -------
     const Double_t alpha = slaGmst(mjd) + obs.GetElong();
     Double_t el;
     slaDe2h(fAlpha-ra, dec, obs.GetPhi(), &az, &el);
     zd = TMath::Pi()/2-el;
     return;
     */
}

MVector3 MVector3::GetRaDec(const MTime &time, MObservatory &obs) const
{
    return GetRaDec(obs, time.GetGmst());
}

MAstroCatalog::MAstroCatalog() : fLimMag(99), fRadiusFOV(99), fToolTip(0), fObservatory(0), fTime(0)
{
    fList.SetOwner();
    fToolTip = new TGToolTip(0, "", 0);
}

MAstroCatalog::~MAstroCatalog()
{
    if (fTime)
        delete fTime;
    if (fObservatory)
        delete fObservatory;

    fToolTip->Hide();
    delete fToolTip;

    DeleteMap();

    // FIXME: There must be an easier way!
    TIter Next(gROOT->GetListOfCanvases());
    TCanvas *c;
    while ((c=(TCanvas*)Next()))
        c->Disconnect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", this,
                      "EventInfo(Int_t,Int_t,Int_t,TObject*)");

}

TString MAstroCatalog::FindToken(TString &line, Char_t tok)
{
    Ssiz_t token = line.First(tok);
    if (token<0)
    {
        const TString copy(line);
        line = "";
        return copy;
    }

    const TString res = line(0, token);
    line.Remove(0, token+1);
    return res;
}

Int_t MAstroCatalog::atoi(const TSubString &sub)
{
    return atoi(TString(sub));
}

Float_t MAstroCatalog::atof(const TSubString &sub)
{
    return atof(TString(sub));
}

Int_t MAstroCatalog::atoi(const TString &s)
{
    return std::atoi(s);
}

Float_t MAstroCatalog::atof(const TString &s)
{
    return std::atof(s);
}

Int_t MAstroCatalog::ReadXephem(TString catalog)
{
    SetBit(kHasChanged);

    gLog << inf << "Reading Xephem catalog: " << catalog << endl;

    ifstream fin(catalog);

    Int_t add =0;
    Int_t cnt =0;
    Int_t line=0;

    Double_t maxmag=0;

    while (1)
    {
        TString row;
        row.ReadLine(fin);
        if (!fin)
            break;

        line++;

        if (row[0]=='#')
            continue;

        TString line(row);

        TString name  = FindToken(line);
        TString dummy = FindToken(line);
        TString r     = FindToken(line);
        TString d     = FindToken(line);
        TString m     = FindToken(line);
        TString epoch = FindToken(line);

        if (name.IsNull() || r.IsNull() || d.IsNull() || m.IsNull() || epoch.IsNull())
        {
            gLog << warn << "Invalid Entry Line #" << line << ": " << row << endl;
            continue;
        }

        cnt++;

        const Double_t mag = atof(m);

        maxmag = TMath::Max(maxmag, mag);

        if (mag>fLimMag)
            continue;

        if (epoch!="2000")
        {
            gLog << warn << "Epoch != 2000... skipped." << endl;
            continue;
        }

        Double_t ra0, dec0;
        MAstro::Coordinate2Angle(r, ra0);
        MAstro::Coordinate2Angle(d, dec0);

        ra0  *= TMath::Pi()/12;
        dec0 *= TMath::Pi()/180;

        MVector3 *star=new MVector3;
        star->SetRaDec(ra0, dec0, mag);
        star->SetName(name);
        if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV)
        {
            delete star;
            continue;
        }

        fList.Add(star);
        add++;
    }
    gLog << inf << "Read " << add << " out of " << cnt << " (Total max mag=" << maxmag << ")" << endl;

    return add;
}

Int_t MAstroCatalog::ReadNGC2000(TString catalog)
{
    SetBit(kHasChanged);

    gLog << inf << "Reading NGC2000 catalog: " << catalog << endl;

    ifstream fin(catalog);

    Int_t add=0;
    Int_t cnt=0;
    Int_t n  =0;

    Double_t maxmag=0;

    while (1)
    {
        TString row;
        row.ReadLine(fin);
        if (!fin)
            break;

        cnt++;

        const Int_t   rah  = atoi(row(13, 2));
        const Float_t ram  = atof(row(16, 4));
        const Char_t  decs = row(22);
        const Int_t   decd = atoi(row(23, 2));
        const Int_t   decm = atoi(row(26, 2));
        const TString m    = row(43, 4);

        if (m.Strip().IsNull())
            continue;

        n++;

        const Double_t mag = atof(m);

        maxmag = TMath::Max(maxmag, mag);

        if (mag>fLimMag)
            continue;

        const Double_t ra  = MAstro::Hms2Rad(rah,  (int)ram, fmod(ram, 1)*60);
        const Double_t dec = MAstro::Dms2Rad(decd, decm, 0, decs);

        MVector3 *star=new MVector3;
        star->SetRaDec(ra, dec, mag);
        star->SetName(row(0, 8));
        if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV)
        {
            delete star;
            continue;
        }

        fList.Add(star);
        add++;
    }

    gLog << inf << "Read " << add << " out of " << n << " (Total max mag=" << maxmag << ")" << endl;

    return add;
}

Int_t MAstroCatalog::ReadBSC(TString catalog)
{
    SetBit(kHasChanged);

    gLog << inf << "Reading Bright Star Catalog (BSC5) catalog: " << catalog << endl;

    ifstream fin(catalog);

    Int_t add=0;
    Int_t cnt=0;
    Int_t n  =0;

    Double_t maxmag=0;

    while (1)
    {
        TString row;
        row.ReadLine(fin);
        if (!fin)
            break;

        cnt++;

        const Int_t   rah    = atoi(row(75, 2));
        const Int_t   ram    = atoi(row(77, 2));
        const Float_t ras    = atof(row(79, 4));
        const Char_t  decsgn = row(83);
        const Int_t   decd   = atoi(row(84, 2));
        const Int_t   decm   = atoi(row(86, 2));
        const Int_t   decs   = atoi(row(88, 2));
        const TString m      = row(102, 5);

        if (m.Strip().IsNull())
            continue;

        n++;

        const Double_t mag = atof(m.Data());

        maxmag = TMath::Max(maxmag, mag);

        if (mag>fLimMag)
            continue;

        const Double_t ra  = MAstro::Hms2Rad(rah, ram, ras);
        const Double_t dec = MAstro::Dms2Rad(decd, decm, decs, decsgn);

        MVector3 *star=new MVector3;
        star->SetRaDec(ra, dec, mag);
        star->SetName(row(4,9));
        if (star->Angle(fRaDec)*TMath::RadToDeg()>fRadiusFOV)
        {
            delete star;
            continue;
        }

        fList.Add(star);
        add++;
    }

    gLog << inf << "Read " << add << " out of " << n << " (Total max mag=" << maxmag << ")" << endl;

    return add;
}

Int_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v, Int_t type)
{
    MVector3 w;

    switch (type)
    {
    case 1:
        w.SetRaDec(v.X()-fRaDec.Phi(), v.Y(), 1);
        break;

    case 2:
        if (!fTime || !fObservatory)
            return kFALSE;
        w.SetZdAz(v.Y(), v.X(), 1);
        w = w.GetRaDec(*fTime, *fObservatory);
        w.RotateZ(-fRaDec.Phi());
        break;

    default:
        return kFALSE;
    }

    w *= rot;

    // Stretch such, that the X-component is alwas the same. Now
    // Y and Z contains the crossing point between the star-light
    // and the plain of a virtual screen (ccd...)
    if (TestBit(kPlainScreen))
        w *= TMath::RadToDeg()/w.X();
    else
        w.SetMag(TMath::RadToDeg());

    v.Set(w(1), w(2));

    if (w(0)<0)
        return kERROR;

    return w(1)>gPad->GetX1() && w(2)>gPad->GetY1() &&
           w(1)<gPad->GetX2() && w(2)<gPad->GetY2();
}

Bool_t MAstroCatalog::DrawLine(const TVector2 &v, Double_t dx, Double_t dy, const TRotation &rot, Int_t type)
{
    const TVector2 add(dx*TMath::DegToRad(), dy*TMath::DegToRad());

    TVector2 v0 = v;
    TVector2 v1 = v+add;

    const Int_t rc0 = Convert(rot, v0, type);
    const Int_t rc1 = Convert(rot, v1, type);
    // Both are kFALSE or both are kERROR
    if ((rc0|rc1)==kFALSE || (rc0&rc1)==kERROR)
        return kFALSE;

    TLine *line = new TLine(v0.X(), v0.Y(), v1.X(), v1.Y());
    line->SetBit(kCanDelete);
    line->SetLineStyle(kDashDotted); //kDashed, kDotted, kDashDotted
    line->SetLineColor(kWhite+type*2);
    line->SetBit(kCannotPick);
    fMapG.Add((Long_t)line, 0);

    const TVector2 deg = v*TMath::RadToDeg();
    TString txt = type==1 ?
        Form("Ra=%.1fh  Dec=%.1fd", fmod(deg.X()/15+48, 24),  fmod(deg.Y()+270,180)-90) :
        Form("Zd=%.1fd  Az=%.1fd",  fmod(deg.Y()+270,180)-90, fmod(deg.X()+720, 360));

    TMarker *tip=new TMarker(v0.X(), v0.Y(), kDot);
    tip->SetBit(kCanDelete);
    tip->SetBit(kCannotPick);
    tip->SetMarkerColor(kWhite+type*2);
    fMapG.Add((Long_t)tip, (Long_t)new TString(txt));

    return kTRUE;
}

void MAstroCatalog::Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type)
{
    const TVector2 v1 = v0 + TVector2(dx[0]*TMath::DegToRad(), dy[0]*TMath::DegToRad());

    //    if (TMath::Abs(v1.Y())>TMath::Pi()/2)
    //        return;

    const Int_t v0x = (int)(v0.X()*TMath::RadToDeg());
    const Int_t v0y = (int)(v0.Y()*TMath::RadToDeg());

    Int_t idx[] = {1, 1, 1, 1};

    Int_t dirs[4][2] = { {0, stepy}, {stepx, 0}, {0, -stepy}, {-stepx, 0} };

    for (int i=0; i<dx.GetSize(); i++)
    {
        for (int j=0; j<4; j++)
        {
            const Bool_t rcx0 = (dx[i]+720)%360==(dx[0]+dirs[j][0]+720)%360;
            const Bool_t rcy0 = (dy[i]+360)%180==(dy[0]+dirs[j][1]+360)%180;
            if (rcx0&&rcy0)
                idx[j] = 0;
        }
    }

    // Enhance size of array by 1, copy current
    // position as last entry
    dx.Set(dx.GetSize()+1);
    dy.Set(dy.GetSize()+1);

    dx[dx.GetSize()-1] = dx[0];
    dy[dy.GetSize()-1] = dy[0];

    // Store current positon
    const Int_t d[2] = { dx[0], dy[0] };

    for (int i=0; i<4; i++)
        if (idx[i])
        {
            // Calculate new position
            //dx[0] = (d[0]+dirs[i][0]+540-v0x)%360-180+v0x;
            //dy[0] = (d[1]+dirs[i][1]+270-v0y)%180- 90+v0y;
            dx[0] = d[0]+dirs[i][0];
            dy[0] = d[1]+dirs[i][1];

            // Draw corresponding line and iterate through grid
            if (DrawLine(v1, dirs[i][0], dirs[i][1], rot, type))
                Draw(v0, rot, dx, dy, stepx, stepy, type);

            dx[0]=d[0]; dy[0]=d[1];
        }
}

void MAstroCatalog::DrawNet(const TVector2 &v0, const TRotation &rot, Int_t type)
{
    //const Double_t step = TMath::DegToRad();

    TArrayI dx(1);
    TArrayI dy(1);

    // align to 1deg boundary
    TVector2 v = v0*TMath::RadToDeg();
    v.Set((Float_t)TMath::Nint(v.X()), (Float_t)TMath::Nint(v.Y()));

    // calculate stepsizes based on visible FOV
    Int_t stepx=1;
    if (fabs(v.Y())>90-fRadiusFOV || fabs(v.Y())<fRadiusFOV)
        stepx = 180/10;
    else
    {
        // This is a rough estimate how many degrees are visible
        const Float_t m = log(fRadiusFOV/180.)/log(90./fRadiusFOV-1);
        const Float_t t = log(180.)-m*log(fRadiusFOV);
        const Int_t n = (Int_t)(exp(m*log(90-fabs(v.Y()))+t)+0.5);
        stepx = n<4 ? 1 : n/4;
    }

    const Int_t n = (Int_t)(fRadiusFOV+1);
    Int_t stepy = n<4 ? 1 : n/4;

    // align stepsizes to be devisor or 180 and 90
    while (180%stepx)
        stepx++;
    while (90%stepy)
        stepy++;

    // align to step-size boundary
    while ((int)(v.X())%stepx)
        v.Set(v.X()+1, v.Y());

    while ((int)(v.Y())%stepy)
        v.Set(v.X(), v.Y()+1);

    // draw...
    v *= TMath::DegToRad();
    Draw(v, rot, dx, dy, stepx, stepy, type);
}

void MAstroCatalog::Paint(Option_t *o)
{
    if (!gPad->IsBatch())
        gVirtualX->ClearWindow();

    if (TestBit(kHasChanged))
        DrawPrimitives(o);
}

void MAstroCatalog::DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t transparent)
{
    const Double_t ra  = v.Phi()*TMath::RadToDeg()/15;
    const Double_t dec = (TMath::Pi()/2-v.Theta())*TMath::RadToDeg();

    TString str = v.GetName();
    str += Form(":  Ra=%.1fh", ra);
    str += Form("  Dec=%.1fd", dec);
    str += Form("  Mag=%.1f", -2.5*log10(v.Mag()));

    // draw star on the camera display
    TMarker *tip=new TMarker(x, y, transparent ? kDot : kFullDotLarge);;
    tip->SetMarkerColor(kBlack);
    tip->SetBit(kCanDelete);
    tip->SetBit(kCannotPick);
    fMapG.Add((Long_t)tip, (Long_t)new TString(str));
}

void MAstroCatalog::Update()
{
    if (gPad && TestBit(kMustCleanup))
    {
        SetBit(kHasChanged);
        gPad->Modified();
    }
}

void MAstroCatalog::SetTime(const MTime &time)
{
    if (fTime)
        delete fTime;
    fTime=(MTime*)time.Clone();
}

void MAstroCatalog::SetObservatory(const MObservatory &obs)
{
    if (fObservatory)
        delete fObservatory;
    fObservatory=(MObservatory*)obs.Clone();
}

void MAstroCatalog::AddPrimitives(Option_t *o)
{
    const Double_t ra  = fRaDec.Phi();
    const Double_t dec = TMath::Pi()/2-fRaDec.Theta();

    // Precalc Sin/Cos...
    TRotation trans;
    trans.Rotate(dec, TVector3(0, 1, 0));

    const TVector3 zdaz0 = fRaDec.GetZdAz(*fTime, *fObservatory);
    TVector2 zdaz(zdaz0.Phi(), zdaz0.Theta());
    MAstroCatalog::DrawNet(zdaz, trans, 2);

    TVector2 radec(ra, dec);
    MAstroCatalog::DrawNet(radec, trans, 1);

    TIter Next(&fList);
    TVector3 *v=0;
    while ((v=(TVector3*)Next()))
    {
        // FIXME: Check Magnitude!
        TVector2 s(v->Phi(), TMath::Pi()/2-v->Theta());
        if (Convert(trans, s, 1)==kTRUE)
            DrawStar(s.X(), s.Y(), *v, kFALSE);
    }
}

void MAstroCatalog::SetRangePad()
{
     const Double_t edge = fRadiusFOV/TMath::Sqrt(2.);
     gPad->Range(-edge, -edge, edge, edge);

     cout << gPad->GetWw() << " " << gPad->GetWh() << endl;

     const Float_t w = gPad->GetWw();
     const Float_t h = gPad->GetWh();

     if (w<h)
         gPad->Range(-edge, -edge*h/w, edge, edge*h/w);
     else
         gPad->Range(-edge*w/h, -edge, edge*w/h, edge);
}

void MAstroCatalog::DrawPrimitives(Option_t *o)
{
    DeleteMap();

    SetRangePad();

    TStopwatch clk;
    clk.Start();
    AddPrimitives(o);
    clk.Stop();
    clk.Print();

    // Append all objects to pad
    DrawMap();

    ResetBit(kHasChanged);
}

void MAstroCatalog::Draw(Option_t *o)
{
    TObject::Draw(o);

    if (!TestBit(kHasChanged))
        DrawPrimitives(o);

    // Connect all TCanvas::ProcessedEvent to this->EventInfo
    // This means, that after TCanvas has processed an event
    // EventInfo of this class is called, see TCanvas::HandleInput
    gPad->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)",
                               "MAstroCatalog", this,
                               "EventInfo(Int_t,Int_t,Int_t,TObject*)");

    // Do this instead of fListG.Draw, because
    // TCollection overwrites Draw
    // Would be nice, but doesn't work because the single
    // graphical object are not handled by TPad anymore...
    //    fListG.AppendPad();
}

// --------------------------------------------------------------------------
//
// This function was connected to all created canvases. It is used
// to redirect GetObjectInfo into our own status bar.
//
// The 'connection' is done in AddTab
//
void MAstroCatalog::EventInfo(Int_t event, Int_t px, Int_t py, TObject *selected)
{
    // Try to find a corresponding object with kCannotPick set and
    // an available TString (for a tool tip)
    if (!selected || selected==this)
    {
        Long_t key, val;
        TExMapIter map(&fMapG);
        while (map.Next(key, val))
        {
            if (!val)
                continue;

            TObject *o=(TObject*)key;
            if (o->DistancetoPrimitive(px, py)>TPad::GetMaxPickDistance())
                continue;

            selected = o;
            break;
        }
    }

    if (!selected)
        return;

    TCanvas *c = (TCanvas*)gTQSender;

    TVirtualPad* save=gPad;

    gPad = c ? c->GetSelectedPad() : NULL;
    if (!gPad)
        return;

    switch (event)
    {
    case kMouseMotion:
        {
            TString *s = (TString*)fMapG.GetValue((Long_t)selected);
            if (!fToolTip->IsMapped() && s)
                ShowToolTip(px, py, *s);
        }
        break;

    case kMouseLeave:
        if (fToolTip->IsMapped())
            fToolTip->Hide();
        break;

    case kKeyPress: // Unresolved keyboard event
        /*
        switch (px)
        {
        case kKey_Plus:
        case kKey_Minus:*/
            ExecuteEvent(kKeyPress, px, py);/*
            break;
        }
        break;*/
    }

    gPad=save;
}

void MAstroCatalog::ExecuteEventKbd(Int_t keycode, Int_t keysym)
{
    Double_t dra =0;
    Double_t ddec=0;

    switch (keysym)
    {
    case kKey_Left:
        dra = -TMath::DegToRad();
        break;
    case kKey_Right:
        dra = +TMath::DegToRad();
        break;
    case kKey_Up:
        ddec = +TMath::DegToRad();
        break;
    case kKey_Down:
        ddec = -TMath::DegToRad();
        break;
    case kKey_Plus:
        SetRadiusFOV(fRadiusFOV+1);
        break;
    case kKey_Minus:
        SetRadiusFOV(fRadiusFOV-1);
        break;

    default:
        return;
    }

    const Double_t r = fRaDec.Phi();
    const Double_t d = TMath::Pi()/2-fRaDec.Theta();

    SetRaDec(r+dra, d+ddec);

    gPad->Update();
}

// ------------------------------------------------------------------------
//
// Execute a gui event on the camera
//
void MAstroCatalog::ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2)
{
    if (!TestBit(kGuiActive))
        return;

    if (event==kKeyPress)
        ExecuteEventKbd(mp1, mp2);

}

Int_t MAstroCatalog::DistancetoPrimitive(Int_t px, Int_t py)
{
    Int_t min = INT_MAX;

    Long_t key, val;
    TExMapIter map(&fMapG);
    while (map.Next(key, val))
    {
        TObject *o=(TObject*)key;

        const Int_t d = o->DistancetoPrimitive(px, py);

        if (d<TPad::GetMaxPickDistance())
            return 0;

        if (d<min)
            min=d;
    }

    return min;
}

void MAstroCatalog::ShowToolTip(Int_t px, Int_t py, const char *txt)
{
    Int_t x=0;
    Int_t y=0;

    const Window_t id1 = gVirtualX->GetWindowID(gPad->GetCanvasID());
    const Window_t id2 = fToolTip->GetParent()->GetId();

    Window_t id3;
    gVirtualX->TranslateCoordinates(id1, id2, px, py, x, y, id3);

    // Show tool tip
    fToolTip->SetText(txt);
    fToolTip->Show(x+4, y+4);
}

