tcsSite.cc

Go to the documentation of this file.

#include <sstream>
#include <stdexcept>
#include "slalib.h"
#include "coordConv/physConst.h"
#include "tcc/telConst.h"
#include "tcc/tcsShim.h"
#include "tcc/tcsSite.h"

namespace {

    void nullRefSub (int mode,
        double refa, double refb,
        double hm, double tlat, double temp, double press,
        double humid, double tlr, double wavel,
        double xin, double yin, double zin,
        double *xout, double *yout, double *zout
    ) {
        *xout = xin;
        *yout = yin;
        *zout = zin;
    }

}

namespace tcc {

    TcsSite::TcsSite( ) :
        timestamp(0)
    {
        for (int i = 0, end = aux.size(); i < end; ++i) {
            aux[i] = 0;
        }
        // this should work, but clang rejects it
        // aux.fill(0);
    };
    
    void TcsSite::update(
            double tai,
            coordConv::Site const &site
    ) {
        int stat;
        stat = tcsItimeo(
            -site.utc_tai,                      // TAI-UTC (s)
            site.ut1_tai - site.utc_tai,        // UT1-UTC (s)
            site.poleX * coordConv::RadPerDeg,  // polar-motion x-component (radians)
            site.poleY * coordConv::RadPerDeg,  // polar-motion y-component (radians)
            &teo);                              // TCSpk earth orientation
        if (stat != 0) {
            std::ostringstream os;
            os << "tcsItimeo error: " << stat;
            throw std::runtime_error(os.str());
        }
        
        stat = tcsIsite(
            &teo,                                   // TCSpk earth orientation
            site.meanLong * coordConv::RadPerDeg,   // longitude (radians, E +ve)
            site.meanLat * coordConv::RadPerDeg,    // geodetic latitude (radians)
            site.elev,                              // height above sea level (m)
            site.wavelen / coordConv::AngstromsPerMicron,   // wavelength (micrometres)
            0,
//            &nullRefSub,                            // special refraction function; disable refraction correction
            &tsite);                                // TCSpk site location and conditions
        if (stat != 0) {
            std::ostringstream os;
            os << "tcsIsite error: " << stat;
            throw std::runtime_error(os.str());
        }

        tsite.press = 0;    // disable refraction correction

        // from tcsSlow (expensive to call and we already have the results)
        double taiDays = tai / coordConv::SecPerDay;
        tsite.t0 = taiDays;
        double ttDays = taiDays + teo.ttmtai;
        tsite.tt0 = ttDays;
        tsite.ttj = slaEpj(ttDays);
        
        double d = floor(tai);
        tsite.st0 = slaDranrm(slaGmsta(d, (tai-d) - teo.delat + teo.delut) + tsite.tlong + slaEqeqx(ttDays));
        tsite.refar = 0;    // disable refraction correction
        tsite.refbr = 0;    // disable refraction correction
        // tcsite.amprms intentionally left blank; it isn't needed and isn't readily available from coordConv

        timestamp = tai;
    }

    std::string TcsSite::__repr__() const {
        std::ostringstream os;
        os << *this;
        return os.str();
    }

    std::ostream &operator<<(std::ostream &os, TcsSite const &tcsSite) {
        os << "# TcsSite" << std::endl;
        
        os << "timestamp " << tcsSite.timestamp << " # timestamp (TAI, MJD seconds)" << std::endl;
        
        os << "tcsSite.teo: TCSpk time information" << std::endl;
        os << tcsSite.teo;
        
        os << "tcsSite.tsite: TCSpk site information" << std::endl;
        os << tcsSite.tsite;

        // start with only three readings; show more if we use more
        for (int i = 0; i < 3; ++i) {
            os << "aux[" << i << "] " << tcsSite.aux.at(i) << " # Auxiliary data readings" << std::endl;
        }

        return os;
    }

}