#include "SlaStars.h"

#include "slalib.h"

ClassImp(SlaStars);

SlaStars::SlaStars() : Slalib()
{
}

SlaStars::~SlaStars()
{
}

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

void SlaStars::Set(const ZdAz &zdaz)
{
    fAltAz = AltAz(kPiDiv2-zdaz.Zd(), zdaz.Az()) * kDeg2Rad;
    fRaDec = CalcRaDec(fAltAz);
}

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

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

    //
    // ----- calculate star independent parameters ----------
    //
    slaMappa(2000.0, mjd, fAmprms);
    slaAoppa(mjd, 0,                    // mjd, UT1-UTC
             GetElong(), GetPhi(), 148, // g”ttingen long, lat, height
             0, 0,                      // polar motion x, y-coordinate (radians)
             273.155, 1013.25, 0.5,     // temp, pressure, humidity
             0.2, 0.0065,               // wavelength, tropo lapse rate
             fAoprms);
}

RaDec SlaStars::CalcRaDec(const AltAz &altaz) const
{
    return CalcRaDec(ZdAz(kPiDiv2-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);
}

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)
              &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(kPiDiv2-zdaz.Zd(), zdaz.Az());
}

ZdAz SlaStars::GetApproxVel(const RaDec &radec) const // [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);
}
