tcsShim.cc

Go to the documentation of this file.

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

namespace tcc {
    namespace tcspk {

        bool tcsFromPVT(double &pos, double &vel, coordConv::PVT const &pvt, double tai, bool isAz) {
            bool isfinite = pvt.isfinite();
            if (isfinite) {
                pos = pvt.getPos(tai) * coordConv::RadPerDeg;
                vel = pvt.vel * coordConv::RadPerDeg * coordConv::SecPerDay;
            } else {
                pos = 0;
                vel = 0;
            }
            if (isAz) {
                pos = coordConv::Pi - pos;
                vel = -vel;
            }
            return isfinite;
        }
    
        void pvtFromTCS(coordConv::PVT &pvt, double posArr[2], double tai, double deltaTAI, bool isAz) {
            // slaDrange wraps an angle into the range [-pi, pi]
            double locPosArr[2];
            for (int i = 0; i < 2; ++i) {
                if (isAz) {
                    locPosArr[i] = coordConv::Pi - posArr[i];
                } else {
                    locPosArr[i] = posArr[i];
                }
            }
            pvt.pos = locPosArr[0] / coordConv::RadPerDeg;
            pvt.vel = slaDrange(locPosArr[1] - locPosArr[0]) / (coordConv::RadPerDeg * deltaTAI);
            pvt.t = tai;
        }
    }

    std::ostream &operator<<(std::ostream &os, ASTROM const &astrom) {
        os << "astrom.cosys  " << astrom.cosys << " # Type of sky reference system; AZEL_TOPO=" << AZEL_TOPO << std::endl;
        os << "astrom.eqx    " << astrom.eqx  << " # Equinox (for mean RA,Dec systems)" << std::endl;
        os << "astrom.wavel  " << astrom.wavel  << " # Wavelength (micrometers)" << std::endl;
        os << "astrom.refa   " << astrom.refa  << " # Refraction coefficient A" << std::endl;
        os << "astrom.refb   " << astrom.refb  << " # Refraction coefficient B" << std::endl;
        os << "astrom.sth    " << astrom.sth  << " # Sine of intermediate z-rotation" << std::endl;
        os << "astrom.cth    " << astrom.cth  << " # Cosine of intermediate z-rotation" << std::endl;
        os << "astrom.spm1   # Sky-patch matrix #1" << std::endl;
        os << formatCMat<double, 3, 3>(astrom.spm1) << std::endl;
        os << "astrom.spm1_i   # Inverse sky-patch matrix #1" << std::endl;
        os << formatCMat<double, 3, 3>(astrom.spm1_i) << std::endl;
        os << "astrom.spm2   # Sky-patch matrix #2" << std::endl;
        os << formatCMat<double, 3, 3>(astrom.spm2) << std::endl;
        os << "astrom.spm2_i   # Inverse sky-patch matrix #2" << std::endl;
        os << formatCMat<double, 3, 3>(astrom.spm2_i) << std::endl;
        return os;
    }

    std::ostream &operator<<(std::ostream &os, FLDOR const &fldor) {
        os << "fldor.sia cia " << fldor.sia << "  " << fldor.cia << " # Sine and cosine of IAA" << std::endl;
        os << "fldor.pai " << fldor.pai << " # Instrument Position-Angle" << std::endl;
        os << "fldor.jf " << fldor.jf << " # Optimization: SLITO = " << SLITO << " = 1D (slit), FIELDO = " << FIELDO << " or else 2D (field)" << std::endl;
        return os;
    }

    std::ostream &operator<<(std::ostream &os, PORIG const &porig) {
        os << "porig.p0  " << formatCArr(porig.p0, 2) << " # (x,y) position (user units)" << std::endl;
        os << "porig.ob  # Offset from base (user units)" << std::endl;
        os << formatCMat<double, 3, 2>(porig.ob) << std::endl;
        os << "porig.p   " << formatCArr(porig.p, 2) << " # Current pointing-origin including offsets (user units)" << std::endl;
        return os;
    }
    
    std::ostream &operator<<(std::ostream &os, SITE const &site) {
        os << "site.tlongm, tlatm " << site.tlongm << ", " << site.tlatm << " # Telescope longitude and latitude (radians)" << std::endl;
        os << "site.tlong,  tlat  " << site.tlong << ", " << site.tlat << " # Telescope longitude and latitude (polar motion adjusted, radians)" << std::endl;
        os << "site.slat, clat " << site.slat << ", " << site.clat << " # Functions of latitude" << std::endl;
        os << "site.hm     " << site.hm << " # Site elevation, metres above sea-level" << std::endl;
        os << "site.diurab " << site.diurab << " # Diurnal aberration amplitude (radians)" << std::endl;
        os << "site.daz    " << site.daz << " # Azimuth correction (terrestrial-celestial, radians)" << std::endl;
        os << "site.temp   " << site.temp << " # Ambient temperature (K)" << std::endl;
        os << "site.press  " << site.press << " # Pressure (mB=hPa)" << std::endl;
        os << "site.humid  " << site.humid << " # Relative humidity (0-1)" << std::endl;
        os << "site.tlr " << site.tlr << " # Tropospheric lapse rate (K per metre)" << std::endl;
        os << "site.wavelr " << site.wavelr << " # Reference wavelength (micron) or minus frequency (GHz)" << std::endl;
        os << "site.refar, refbr " << site.refar << ", " << site.refbr << " # Refraction coefficients for wavelr (radians)" << std::endl;
        os << "site.rfun   " << site.rfun << " # Optional refraction function" << std::endl;
        os << "site.t0     " << site.t0 << " # Raw clock time at reference time" << std::endl;
        os << "site.st0    " << site.st0 << " # LAST at reference time (radians)" << std::endl;
        os << "site.tt0    " << site.tt0 << " # TT at reference time (MJD)" << std::endl;
        os << "site.ttj    " << site.ttj << " # TT at reference time (Julian Epoch)" << std::endl;
        return os;
    }

    std::ostream &operator<<(std::ostream &os, TARG const &targ) {
        os << "targ.p0    " << formatCArr(targ.p0, 2) << " # Position at reference time (spherical coordinates, radians)" << std::endl;
        os << "targ.dt    " << formatCArr(targ.dt, 2) << " # Differential rates (radians/day)" << std::endl;
        os << "targ.t0    " << targ.t0 << " # Reference time (TAI MJD days)" << std::endl;
        os << "targ.ob # Offsets from base (radians, additive)" << std::endl;
        os << formatCMat<double, 3, 2>(targ.ob) << std::endl;
        os << "targ.op0   " << formatCArr(targ.op0, 2) << " # Position at reference time including offsets from base" << std::endl;
        os << "targ.op    " << formatCArr(targ.p, 2) << " # Current position including offsets and differential tracking" << std::endl;
        return os;
    }

    std::ostream &operator<<(std::ostream &os, TIMEO const &timeo) {
        os << "timeo.ttmtai " << timeo.ttmtai << " # TT-TAI (day)" << std::endl;
        os << "timeo.delat  " << timeo.delat  << " # TAI-UTC (day)" << std::endl;
        os << "timeo.delut  " << timeo.delut  << " # UT1-UTC (day)" << std::endl;
        os << "timeo.xpmr, ypmr " << timeo.xpmr << "  " << timeo.ypmr << " # Polar motion angle in x, y (radians)" << std::endl;
        return os;
    }

    std::ostream &operator<<(std::ostream &os, TSCOPE const &tscope) {
        os << "tscope.fl " << tscope.fl << " # Telescope focal length (user units)" << std::endl;
        os << "tscope.mount " << tscope.mount << " # Mount type; ALTAZ=" << ALTAZ << std::endl;
        os << "tscope.rotl " << tscope.rotl << " # Location of instrument rotator:"
            << " OTA=" << OTA
            << " NASMYTH_L=" << NASMYTH_L
            << " NASMYTH_R=" << NASMYTH_R
            << std::endl;
        os << "tscope.rnogo " << tscope.rnogo << " # Mount pole avoidance distance (radians)" << std::endl;
        os << "tscope.ae2nm  # Rotation matrix, [Az,El] to nominal mount" << std::endl;
        os << formatCMat<double, 3, 3>(tscope.ae2nm) << std::endl;
        os << "tscope.ae2mt  # Rotation matrix, [Az,El] to actual mount" << std::endl;
        os << formatCMat<double, 3, 3>(tscope.ae2mt) << std::endl;
        os << "tscope.ia " << tscope.ia << " # Roll zero point (radians)" << std::endl;
        os << "tscope.ib " << tscope.ib << " # Pitch zero point (radians)" << std::endl;
        os << "tscope.np " << tscope.np << " # Mount axes nonperpendicularity (radians)" << std::endl;
        os << "tscope.xt " << tscope.xt << " # Telescope vector, x" << std::endl;
        os << "tscope.yt " << tscope.yt << " # Telescope vector, y" << std::endl;
        os << "tscope.zt " << tscope.zt << " # Telescope vector, z" << std::endl;
        return os;
    }

    template <typename T>
    std::string formatCArr(T const*arr, int n) {
        std::ostringstream os;
        os << std::setprecision(12);
        for (int i = 0; i < n; ++i) {
            if (i > 0) {
                os << "  ";
            }
            os << arr[i];
        }
        return os.str();
    }
    
    template <typename T, int ROWS, int COLS>
    std::string formatCMat(T const mat[ROWS][COLS]) {
        std::ostringstream os;
        for (int r = 0; r < ROWS; ++r) {
            os << "  " << formatCArr(mat[r], COLS);
            if (r < ROWS-1) {
                os << std::endl;
            }
        }
        return os.str();
    }

}