/* ======================================================================== *\
!
! *
! * 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 <TLine.h>
#include <TMarker.h>
#include <TArrayI.h>
#include <TRotation.h>
#include <TStopwatch.h>
#include <TVirtualPad.h>

#include <stdlib.h>

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

#include "MAstro.h"
#include "MTime.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());
}

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)
{
    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);
        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)
{
    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);
        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)
{
    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);
        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;
}

Bool_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v, Int_t type)
{
    MVector3 w;
    w.SetRaDec(v.Y(), v.X()-fRaDec.Phi(), v.Y());
    w *= rot;

    v.Set(w.Phi(), w.Theta());

    return w.Angle(TVector3(1, 0, 0))*TMath::RadToDeg()<=fRadiusFOV;
}

Bool_t MAstroCatalog::PaintLine(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 Bool_t rc0 = Convert(rot, v0, type);
    const Bool_t rc1 = Convert(rot, v1, type);
    if (!rc0 && !rc1)
        return kFALSE;

    TLine line;
    line.SetLineColor(kGreen);
    line.PaintLine(v0.X(), TMath::Pi()/2-v0.Y(),
                   v1.X(), TMath::Pi()/2-v1.Y());

    return kTRUE;
}

void MAstroCatalog::Paint(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Byte_t type)
{
    Int_t idx[] = {1, 1, 1, 1};

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

    for (int i=0; i<dx.GetSize(); i++)
    {
        for (int j=0; j<4; j++)
            if (dx[i]==dx[0]+dirs[j][0] && dy[i]==dy[0]+dirs[j][1])
                idx[j] = 0;
    }

    TVector2 v1 = v0 + TVector2(dx[0]*TMath::DegToRad(), dy[0]*TMath::DegToRad());

    for (int i=0; i<4; i++)
        if (idx[i])
        {
            dx.Set(dx.GetSize()+1);
            dy.Set(dy.GetSize()+1);

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

            dx[0] += dirs[i][0];
            dy[0] += dirs[i][1];

            if (PaintLine(v1, dirs[i][0], dirs[i][1], rot, type))
                Paint(v0, rot, dx, dy, type);

            dx[0] -= dirs[i][0];
            dy[0] -= dirs[i][1];
        }
}

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

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

    TVector2 v = v0*TMath::RadToDeg();
    v.Set(TMath::Floor(v.X()), TMath::Floor(v.Y()));
    v *= TMath::DegToRad();

    Paint(v, rot, dx, dy, type);
}

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

    TIter Next(&fList);
    TVector3 *v;

    Double_t minra=360, maxra=0, mindec=360, maxdec=0;

    while ((v=(TVector3*)Next()))
    {
        minra = TMath::Min(minra, v->Phi());
        maxra = TMath::Max(maxra, v->Phi());

        mindec = TMath::Min(mindec, TMath::Pi()/2-v->Theta());
        maxdec = TMath::Max(maxdec, TMath::Pi()/2-v->Theta());
    }

    cout << gPad << endl;

    cout << "Minra: " << (minra-ra)*TMath::RadToDeg() << endl;
    cout << "Maxra: " << (maxra-ra)*TMath::RadToDeg() << endl;

    cout << "Mindec: " << (mindec-dec)*TMath::RadToDeg() << endl;
    cout << "Maxdec: " << (maxdec-dec)*TMath::RadToDeg() << endl;

    //gPad->Range(minra-ra, mindec-dec, maxra-ra, maxdec-dec);
    gPad->Range(-fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.), -fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.),
                fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.),  fRadiusFOV*TMath::DegToRad()/TMath::Sqrt(2.));

    Next.Reset();

    cout << "Dec: " << dec*TMath::RadToDeg() << endl;
    cout << "Ra:  " << ra*TMath::RadToDeg() << endl;

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

    TStopwatch clk;
    clk.Start();

    TMarker mark;
    mark.SetMarkerColor(kBlack);
    mark.SetMarkerStyle(kCircle);
    while ((v=(TVector3*)Next()))
    {
        MVector3 v0;
        v0.SetRaDec(v->Phi()-ra, TMath::Pi()/2-v->Theta(), 1);
        v0 *= trans;

        mark.SetMarkerSize((fLimMag-TMath::Log(v->Mag()))/4);
        mark.PaintMarker(v0.Phi(), TMath::Pi()/2-v0.Theta());
    }

    TMarker m;
    m.SetMarkerStyle(kCross);
    m.PaintMarker(0, 0);

    m.SetMarkerColor(kRed);
    m.SetMarkerStyle(kFullDotSmall);

    TVector2 v0(ra, dec);
    PaintNet(v0, trans);

    clk.Stop();
    clk.Print();
}
