#include "SlaStars.h"

#include "slalib.h"

#include "MPointing.h"

ClassImp(SlaStars);

SlaStars::SlaStars(MObservatory::LocationName_t key) : Slalib(key)
{
}

SlaStars::~SlaStars()
{
}

void SlaStars::Set(const AltAz &altaz)
{
    fAltAz = altaz * TMath::DegToRad();
    fRaDec = CalcRaDec(fAltAz);
}

void SlaStars::Set(const ZdAz &zdaz)
{
    fAltAz = AltAz(TMath::Pi()/2-zdaz.Zd(), zdaz.Az()) * TMath::DegToRad();
    fRaDec = CalcRaDec(fAltAz);
}

void SlaStars::Set(const RaDec &radec)
{
    fRaDec = radec * TMath::DegToRad();
    fAltAz = CalcAltAz(fRaDec);
}

void SlaStars::SetMjd(double mjd)
{
    Slalib::SetMjd(mjd);

    //
    // ----- calculate star independent parameters ----------
    //

    // prepare calculation: Mean Place to geocentric apperent
    slaMappa(2000.0, GetMjd(), fAmprms);   // Epoche, TDB

    // prepare: Apperent to observed place
    slaAoppa(GetMjd(), 0,                  // mjd, Delta UT=UT1-UTC
             GetElong(), GetPhi(), GetHeight(), // long, lat, height
             0, 0,                         // polar motion x, y-coordinate (radians)
             // 273.155, 1013.25, 0.5,     // temp, pressure, humidity
             273.155+20, 1013.25, 0.5,     // temp, pressure, humidity
             // 0.2, 0.0065,               // wavelength, tropo lapse rate
             0.55, 0.0065,                 // wavelength, tropo lapse rate
             fAoprms);
}

RaDec SlaStars::CalcRaDec(const AltAz &altaz) const
{
    return CalcRaDec(ZdAz(TMath::Pi()/2-altaz.Alt(), altaz.Az()));
}

RaDec SlaStars::CalcRaDec(const ZdAz &zdaz) const
{
    //
    // -- observed to apparent --
    //  Workaraound for slalib: discard const
    //
    double r=0, d=0;
    slaOapqk((char*)"A", zdaz.Az(), zdaz.Zd(), (double*)fAoprms, &r, &d);

    //
    // -- apparent to mean --
    //  Workaraound for slalib: discard const
    //
    double ra, dec;
    slaAmpqk(r, d, (double*)fAmprms, &ra, &dec);

    return RaDec(ra, dec);
}

RaDec SlaStars::CalcRaDecFast(const AltAz &altaz) const
{
    //
    // This function does a coordinate system transformation only.
    // This is very fast compared to a correct tranformation with all
    // effects, but much less accurate.
    //
    // It transforms Altitude/Azimuth [rad] into
    // Right Ascension/Declination [rad] 
    //
    double ha, dec;
    slaDh2e(altaz.Az(), altaz.Alt(), GetPhi(), &ha, &dec);
    return RaDec(GetAlpha()-ha, dec);
}

RaDec SlaStars::CalcRaDecFast(const ZdAz &zdaz) const
{
    //
    // This function does a coordinate system transformation only.
    // This is very fast compared to a correct tranformation with all
    // effects, but much less accurate.
    //
    // It transforms Zenith Distance/Azimuth [rad] into
    // Right Ascension/Declination [rad] 
    //
    return CalcRaDecFast(AltAz(TMath::Pi()/2-zdaz.Zd(), zdaz.Az()));
}

ZdAz SlaStars::CalcZdAzFast(const RaDec &radec) const
{
    //
    // This function does a coordinate system transformation only.
    // This is very fast compared to a correct tranformation with all
    // effects, but much less accurate.
    //
    // It transforms Right Ascension/Declination [rad] into
    // zenith Distance/Azimuth [rad]
    //
    AltAz altaz = CalcAltAzFast(radec);
    return ZdAz(TMath::Pi()/2-altaz.Alt(), altaz.Az());
}

AltAz SlaStars::CalcAltAzFast(const RaDec &radec) const
{
    //
    // This function does a coordinate system transformation only.
    // This is very fast compared to a correct tranformation with all
    // effects, but much less accurate.
    //
    // It transforms Right Ascension/Declination [rad] into
    // Altitude/Azimuth [rad]
    //
    double az, el;
    slaDe2h(GetAlpha()-radec.Ra(), radec.Dec(), GetPhi(), &az, &el);
    return AltAz(el, az);
}

ZdAz SlaStars::CalcZdAz(const RaDec &radec) const
{
    //
    // ---- Mean to apparent ----
    //

    double r=0, d=0;
    slaMapqkz(radec.Ra(), radec.Dec(), (double*)fAmprms, &r, &d);
    //
    // Doesn't work - don't know why
    //
    //    slaMapqk (radec.Ra(), radec.Dec(), rdpm.Ra(), rdpm.Dec(),
    //              0, 0, (double*)fAmprms, &r, &d);
    //

    //
    // -- apparent to observed --
    //
    double r1=0;  // ra
    double d1=0;  // dec
    double h0=0;  // ha

    double zd;
    double az;
    slaAopqk (r, d, (double*)fAoprms,
              &az,    // observed azimuth (radians: N=0,E=90) [-pi, pi]
              &zd,    // observed zenith distance (radians)   [-pi/2, pi/2]
              &h0,    // observed hour angle (radians)
              &d1,    // observed declination (radians)
              &r1);   // observed right ascension (radians)

    return ZdAz(zd, az);
}
AltAz SlaStars::CalcAltAz(const RaDec &radec) const
{
    ZdAz zdaz = CalcZdAz(radec);
    return AltAz(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
}

ZdAz SlaStars::GetApproxVel(const RaDec &radec) const // [rad/rad]
{
    // radec        [rad]
    // GetApproxVel [rad/rad]
    double az, vaz, aaz;
    double el, vel, ael;
    double pa, vpa, apa;
    slaAltaz(GetAlpha()-radec.Ra(), radec.Dec(), GetPhi(),
             &az, &vaz, &aaz,
             &el, &vel, &ael,
             &pa, &vpa, &apa);

    return ZdAz(-vel, vaz);
}