source: trunk/MagicSoft/Cosy/catalog/StarCatalog.cc@ 5842

#include "StarCatalog.h"

#include <iomanip.h>   // cout
#include <iostream.h>  // cout

#include <TSystem.h>
#include <TRotation.h>

#include "slalib.h"
#include "slamac.h"
#include "File.h"

#include "MStarList.h"
#include "MAstroCatalog.h"

ClassImp(StarCatalog);

StarCatalog::StarCatalog(MObservatory::LocationName_t key) : SlaStars(key), fW(768), fH(576), fAstro(0), /*fSao(NULL), fSrt(NULL), fEntries(0),*/ fSinAngle(0), fCosAngle(1), fBox(768)
{
    fAstro = new MAstroCatalog;
    fAstro->SetObservatory(*this);
    fAstro->SetPlainScreen();
}

StarCatalog::~StarCatalog()
{
    delete fAstro;
}

void StarCatalog::SetPixSize(const double pixsize)
{
    // pixsize [arcsec/pixel]
    fPixSize = D2PI/360.0*pixsize/3600; // [rad / (deg*pixel)]
    fAstro->SetRadiusFOV(pixsize, 768, 576);
}

double StarCatalog::GetPixSize() const
{
    return fPixSize*3600*360.0/D2PI;
}

void StarCatalog::SetLimitMag(const float mag)
{
    fLimitMag = mag;
    fAstro->SetLimMag(mag);
}

void StarCatalog::SetMjd(double mjd)
{
    SlaStars::SetMjd(mjd);
    fAstro->SetTime(MTime(mjd));
}

void StarCatalog::SetAltAz(const AltAz &altaz)
{
    fAltAz = altaz * D2PI/360.0;
    fRaDec = CalcRaDec(fAltAz);

    fAstro->SetRaDec(fRaDec.Ra(), fRaDec.Dec());
}

void StarCatalog::Reload()
{
    fAstro->SetLimMag(99);
    //fAstro->ReadBSC("bsc5.dat");
    //fAstro->ReadHeasarcPPM("heasarc_ppm.tdat");
    fAstro->ReadCompressed("ppm9.bin");
    fAstro->SetLimMag(fLimitMag);
}

void StarCatalog::SetRaDec(const RaDec &radec)
{
    const RaDec rd = fRaDec*360.0/D2PI;;

    const Bool_t same =
        rd.Ra() >radec.Ra() -1e-5 && rd.Ra() <radec.Ra() +1e-5 &&
        rd.Dec()>radec.Dec()-1e-5 && rd.Dec()<radec.Dec()+1e-5;

    fRaDec = radec * D2PI/360.0;
    fAltAz = CalcAltAz(fRaDec);

    fAstro->SetRaDec(fRaDec.Ra(), fRaDec.Dec());
    if (!same)
        Reload();
}

void StarCatalog::DrawCross(byte *img, const int x, const int y)
{
    for (int dx=-4; dx<5; dx++)
        if (dx) img[y*768+x+dx] = 0xff;

    for (int dy=-4; dy<5; dy++)
        if (dy) img[(y+dy)*768+x] = 0xff;
}

void StarCatalog::GetImg(byte *img, byte *cimg, MStarList &list) const
{
    memset(cimg, 0, 768*576);

    DrawStars(list, cimg);
    DrawCross(img, 768/2, 576/2);
}

void StarCatalog::DrawCircle(int color, byte *img, int xx, int yy, int size)
{
    for (int x=xx-size; x<xx+size+1; x++)
    {
        if (x<0 || x>767)
            continue;

        const float p = xx+size-x;
        const float q = 2*size - p;
        const int h = (int)sqrt(p*q);

        const int y1 = yy-h;
        if (y1>=0 && y1<576)
            img[y1*768+x] = color;

        const int y2 = yy+h;
        if (y2>=0 && y2<576)
            img[y2*768+x] = color;
    }
}

void StarCatalog::PaintImg(unsigned char *buf, int w, int h)
{
    fAstro->PaintImg(buf, w, h);
}

void StarCatalog::DrawStars(MStarList &list, byte *img)
{
    MStarListIter Next(&list);

    MStar *star;
    while ((star=Next()))
    {
        const int mag = (10 - (star->GetMag()>1 ? (int)star->GetMag() : 1))/2;
        Double_t color = 0xf0; //0x0f;       
        // DrawStars flips the picture in X defaultwise now
        DrawCircle(color, img, 768-(int)star->GetX(), (int)star->GetY(), mag);
    }
}

void StarCatalog::CalcStars(MStarList &list)
{
    // full FOV
    fBox=768;
    CalcStars(list, 0, 576, 0, 0);
}

void StarCatalog::CalcStars(MStarList &list, int xc, int yc, 
                            int xo, int yo) const
{

    // CalcStars flips the picture in X defaultwise now
    int x0 = (768-xc)-fBox;
    int x1 = (768-xc)+fBox;
    int y0 = yc-fBox;
    int y1 = yc+fBox;

    if (x0<0) x0=0;
    if (y0<0) y0=0;
    if (x1>fW) x1=fW;
    if (y1>fH) y1=fH;

    // Align stars into telescope system
    // (Move the telescope to pointing position)
    TRotation align;
    align.RotateZ(-fAltAz.Az());
    align.RotateY(-(TMath::Pi()/2-fAltAz.Alt()));
    align.RotateZ(TMath::Pi()/2);

    // For an apropriate unit conversion to pixels
    const Double_t scale = TMath::RadToDeg()*sqrt(768*768 + 576*576)/(fAstro->GetRadiusFOV()*2);

    // Get List of stars from catalog
    TIter Next(fAstro->GetList());
    TVector3 *star=0;

    const Double_t limmag = pow(10, -fLimitMag/2.5);

    while ((star=(TVector3*)Next()))
    {
        // Check for limiting magnitude
        const Double_t mag = star->Mag();
        if (mag < limmag)
            continue;

        // Get star position and do an apropiate
        // conversion to local coordinates
        const RaDec rd(star->Phi(), TMath::Pi()/2-star->Theta());
        const ZdAz  za(CalcZdAz(rd));

        // Virtually move telescope to pointing position
        TVector3 loc;
        loc.SetMagThetaPhi(1, za.Zd(), za.Az());
        loc *= align;

        // Sanity check
        if (loc(2)<0)
            continue;

        // Stretch such, that the Z-component is alwas the same. Now
        // X and Y contains the intersection point between the star-light
        // and the plain of a virtual plain screen (ccd...)
        loc *= 1./loc(2);

        // Do an apropriate unit conversion to pixels
        loc *= scale;

        // if (loc.Mod2()>fRadiusFOV*fRadiusFOV)
        //     continue;

        // Rotate by the rotation angle of the video camera
        Float_t xx = loc.X()*fCosAngle - loc.Y()*fSinAngle;
        Float_t yy = loc.X()*fSinAngle + loc.Y()*fCosAngle;

        if (xx<(x0-768/2) || xx >=(x1-768/2) 
            || yy<(y0-(576/2+yo)) || yy>=(y1-(576/2+yo)))
            continue;

        // Store pixel coordinates of star in list
        list.Add(xx+768/2+xo, yy+576/2+yo, -2.5*log10(mag));

    }
}

AltAz StarCatalog::CalcAltAzFromPix(Double_t pixx, Double_t pixy) const
{
    double dx =  (pixx-576/2)*fCosAngle + (pixy-768/2)*fSinAngle;
    double dy = -(pixx-576/2)*fSinAngle + (pixy-768/2)*fCosAngle;

    dx *= fPixSize;
    dy *= fPixSize;

    //const double dx = (pixx-768.0)*fPixSize + fWidth+DPI;
    //const double dy =         pixy*fPixSize - fHeight;

    double ha, dec;
    slaDh2e(dx, dy, DPI/2-fAltAz.Alt(), &ha, &dec);

    return AltAz(-dec, ha+fAltAz.Az());
}