source: trunk/Cosy/catalog/StarCatalog.cc@ 19661

#include "StarCatalog.h"

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

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

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

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

ClassImp(StarCatalog);

using namespace std;

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

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

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

double StarCatalog::GetPixSize() const
{
    return fPixSize*3600*TMath::RadToDeg();
}

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 * TMath::DegToRad();
    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->Delete();
    fAstro->ReadCompressed("ppm9.bin");
    fAstro->SetLimMag(fLimitMag);
}

void StarCatalog::SetRaDec(const RaDec &radec)
{
    const RaDec rd = fRaDec*TMath::RadToDeg();

    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 * TMath::DegToRad();
    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;
        Int_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/2;
    CalcStars(list, 768/2, 576/2, 0, 0);
}
*/

void StarCatalog::CalcStars(MStarList &list, int xc, int yc, 
                            float offx, float offy) const
{
//    const Int_t offx = 768/2 + xo;
//    const Int_t offy = 576/2 + yo;

    // Allow catalog stars to be a bit outside [0.2deg] of the
    // monitored window. To get the std behaviour set offset=0
    const Int_t offset = TMath::Nint(0.2*TMath::DegToRad()/fPixSize);
    const Int_t boxx   = fBoxX+offset;
    const Int_t boxy   = fBoxY+offset;

    // CalcStars flips the picture in X defaultwise now
    // This defined the box in which stars are really returned
    const int x0 = TMath::Max((768-xc)-boxx,   -offset);
    const int y0 = TMath::Max(yc-boxy,         -offset);
    const int x1 = TMath::Min((768-xc)+boxx, fW+offset);
    const int y1 = TMath::Min(yc+boxy,       fH+offset);
/*
    const int x0 = TMath::Max((768-xc)-box,   -offset);
    const int x1 = TMath::Min((768-xc)+box, fW+offset);

    const int y0 = TMath::Max(yc-box-100,         -offset);
    const int y1 = TMath::Min(yc+box+100,       fH+offset);
 */

/*
    // 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);
 */

    TRotation rot;
    rot.RotateY(-(TMath::Pi()/2-fAltAz.Alt()));

    // 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));

        TVector3 v;
        //v.SetMagThetaPhi(1., TMath::Pi()/2-za.Alt(), za.Az()-fAltAz.Az());
        v.SetMagThetaPhi(1., za.Zd(), za.Az()-fAltAz.Az());
        v *= rot;

        if (v(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...)
        v *= 1./v(2);

        // Do unit conversion to pixels
        v *= 1./fPixSize;

        const Double_t x = -v.Y();
        const Double_t y =  v.X();

/*
        // 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 *= 1./fPixSize;

        const Double_t x = loc.X();
        const Double_t y = loc.Y();
*/
        // if (loc.Mod2()>fRadiusFOV*fRadiusFOV)
        //     continue;


        // Rotate by the rotation angle of the video camera
        // and add the offsets on both axis
        const Double_t xx = x*fCosAngle - y*fSinAngle + 768 - offx;
        const Double_t yy = x*fSinAngle + y*fCosAngle + offy;

        // Check if the resulting star is in the
        // search box for the real stars
        if (xx<x0 || xx>=x1 || yy<y0 || yy>=y1)
            continue;

        // Store pixel coordinates of star in list
        list.Add(xx, yy, -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());
}
*/

ZdAz StarCatalog::CalcDeltaZdAzFromPix(Double_t dpixx, Double_t dpixy) const
{
    double dx =  dpixx*fCosAngle + dpixy*fSinAngle;
    double dy = -dpixx*fSinAngle + dpixy*fCosAngle;

    TVector3 loc(dy, -dx, 1./fPixSize);

    loc.RotateY(TMath::Pi()/2-fAltAz.Alt());

    return ZdAz(loc.Theta()-TMath::Pi()/2+fAltAz.Alt(), -loc.Phi());

    /*
    // 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);


    TVector3 loc(dx, dy, 1./fPixSize);

    loc *= align.Inverse();

    cout << (TMath::Pi()/2-loc.Theta()-alt)*TMath::RadToDeg() << " " << (loc.Phi()-az)*TMath::RadToDeg() << endl;



    TVector3 loc(dx, -dy, 1./fPixSize);

    loc *= align.Inverse();

    return ZdAz(loc.Theta(), loc.Phi());*/
}