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

#include "StarCatalog.h"

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

#include <TSystem.h>

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

#include "MStarList.h"

ClassImp(StarCatalog);

StarCatalog::StarCatalog(MObservatory::LocationName_t key) : SlaStars(key), fSao(NULL), fSrt(NULL), fEntries(0), fSinAngle(0), fCosAngle(1)
{
    // p = pointer to MainFrame (not owner)

    //
    // read index file
    //
    File idx("sao/sao-sort.idx", "r");
    if (!idx)
        return;

    while (!idx.Eof())
    {
        idx.Newline();
        fEntries++;
    }

    idx.Reset();

    fSrt = new sort_t[fEntries];

    for (int i=0; i<fEntries; i++)
    {
        fSrt[i].ra  = idx.Geti(4);
        fSrt[i].dec = idx.Geti(4);
        fSrt[i].nr  = idx.Geti(10);
        idx.Newline();
    }

    //
    // open catalog
    //
    fSao = new SaoFile("sao/sao-sort.cmp");
}

StarCatalog::~StarCatalog()
{
    if (fSrt)
        delete fSrt;
    if (fSao)
        delete fSao;
}

void StarCatalog::SetPixSize(const double pixsize)
{
    fPixSize = D2PI/360.0 * pixsize;

    fWidth  = fPixSize * 768/2;
    fHeight = fPixSize * 576/2;
}

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


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

    cout << "Set --> Alt: " << 360.0/D2PI*fAltAz.Alt();
    cout << "  Az: " << fAltAz.Az() << endl;

    fRaDec = CalcRaDec(fAltAz);

    cout << "Ra: " << 360.0/D2PI*fRaDec.Ra();
    cout << "  Dec: " << 360.0/D2PI*fRaDec.Dec() << endl;

    CalcAltAzRange();
    CalcRaDecRange();
}

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

    fAltAz = CalcAltAz(fRaDec);

    //cout << "Alt: " << 360.0/D2PI*fAltAz.Alt() << "  ";
    //cout << "Az: "  << 360.0/D2PI*fAltAz.Az()  << endl;

    CalcRaDecRange();
    CalcAltAzRange();
}

void StarCatalog::CalcAltAzRange()
{
    byte fAlt0[180];

    for (int h=0; h<180; h++)
        fAlt0[h] = kFALSE;

    for (int h=0; h<360; h++)
        fAz0[h] = kFALSE;

    double az0, alt0;
    double az1, alt1;
    //
    // scan horizontal border
    //
    for (int x=-768/2; x<768/2+1; x++)
    {
        slaDh2e(DPI+x*fPixSize, -fHeight, DPI/2-fAltAz.Alt(), &az0, &alt0);
        slaDh2e(DPI+x*fPixSize, +fHeight, DPI/2-fAltAz.Alt(), &az1, &alt1);

        const int z0 = ((int)(360.0/D2PI*(az0+fAltAz.Az()))+360)%360;
        const int t0 = (int)(360.0/D2PI*alt0);

        fAz0[z0] = kTRUE;

        if (-89<=t0 && t0<=90)
            fAlt0[90-t0] = kTRUE;

        const int z1 = ((int)(360.0/D2PI*(az1+fAltAz.Az()))+360)%360;
        const int t1 = (int)(360.0/D2PI*alt1);

        fAz0[z1] = kTRUE;

        if (-89<=t1 && t1<=90)
            fAlt0[90-t1] = kTRUE;
    }

    //
    // scan vertical border
    //
    for (int y=-576/2; y<576/2+1; y++)
    {
        slaDh2e(DPI-fWidth, y*fPixSize, DPI/2-fAltAz.Alt(), &az0, &alt0);
        slaDh2e(DPI+fWidth, y*fPixSize, DPI/2-fAltAz.Alt(), &az1, &alt1);

        const int z0 = ((int)(360.0/D2PI*(az0+fAltAz.Az()))+360)%360;
        const int t0 = (int)(360.0/D2PI*alt0);

        fAz0[z0] = kTRUE;

        if (-89<=t0 && t0<=90)
            fAlt0[90-t0] = kTRUE;

        const int z1 = ((int)(360.0/D2PI*(az1+fAltAz.Az()))+360)%360;
        const int t1 = (int)(360.0/D2PI*alt1);

        fAz0[z1] = kTRUE;

        if (-89<=t1 && t1<=90)
            fAlt0[90-t1] = kTRUE;
    }

    //
    // count degrees of azimut
    //
    fAzCnt=0;
    for (int x=0; x<360; x++)
        if (fAz0[x])
            fAzCnt++;

    //cout << "fAzCnt: " << setw(3) << fAzCnt << "  " << flush;

    //
    // calculate min and max of altitude
    //
    fAltMin=0;
    fAltMax=0;
    for (int y=0; y<180; y++)
    {
        if (fAlt0[y])
            fAltMax = y;

        if (fAlt0[179-y])
            fAltMin = 179-y;
    }

    fAltMin -= 90;
    fAltMax -= 90;

    //
    // check whether altaz north- or south-pole is in the visible region
    //
    byte img[768*576];
    if (DrawAltAz(0, img, 90, 0))
    {
        fAltMax=89;
        cout << "Alt Az Pole1 Inside!" << endl;
    }
    if (DrawAltAz(0, img, -90, 0))
    {
        fAltMin=-90;
        cout << "Alt Az Pole2 Inside!" << endl;
    }

    //cout << "fAltMin: " << setw(3) << fAltMin << "  ";
    //cout << "fAltMax: " << setw(3) << fAltMax << endl;
}

void StarCatalog::CalcRaDecRange()
{
    //
    // calculate range to search in
    //
    byte fDec[180];

    for (int h=0; h<180; h++)
        fDec[h] = kFALSE;

    for (int h=0; h<360; h++)
        fRa0[h] = kFALSE;

    double ha0, ha1;
    double de0, de1;

    const double phi   = GetPhi();
    const double alpha = GetAlpha();
    //
    // scan horizontal border
    //
    for (int x=-768/2; x<768/2+1; x++)
    {
        double dx, dy;
        slaDh2e(DPI-x*fPixSize, -fHeight, DPI/2-fAltAz.Alt(), &dx, &dy);
        slaDh2e(fAltAz.Az()+dx, -dy, phi, &ha0, &de0);

        slaDh2e(DPI-x*fPixSize, +fHeight, DPI/2-fAltAz.Alt(), &dx, &dy);
        slaDh2e(fAltAz.Az()+dx, -dy, phi, &ha1, &de1);

        const int h0 = ((int)(360.0/D2PI*(alpha-ha0))+360)%360;
        const int d0 = (int)(360.0/D2PI*de0);

        fRa0[h0] = kTRUE;

        if (-90<=d0 && d0<=89)
            fDec[d0+90] = kTRUE;

        const int h1 = ((int)(360.0/D2PI*(alpha-ha1))+360)%360;
        const int d1 = (int)(360.0/D2PI*de1);

        fRa0[h1] = kTRUE;

        if (-90<=d1 && d1<=89)
            fDec[d1+90] = kTRUE;
    }

    //
    // scan vertical border
    //
    for (int y=-576/2; y<576/2+1; y++)
    {
        double dx, dy;
        slaDh2e(DPI-fWidth, -y*fPixSize, DPI/2-fAltAz.Alt(), &dx, &dy);
        slaDh2e(fAltAz.Az()+dx, -dy, phi, &ha0, &de0);

        slaDh2e(DPI+fWidth, -y*fPixSize, DPI/2-fAltAz.Alt(), &dx, &dy);
        slaDh2e(fAltAz.Az()+dx, -dy, phi, &ha1, &de1);

        const int h0 = ((int)(360.0/D2PI*(alpha-ha0))+360)%360;
        const int d0 = (int)(360.0/D2PI*de0);

        fRa0[h0] = kTRUE;

        if (-90<=d0 && d0<=89)
            fDec[d0+90] = kTRUE;

        const int h1 = ((int)(360.0/D2PI*(alpha-ha1))+360)%360;
        const int d1 = (int)(360.0/D2PI*de1);

        fRa0[h1] = kTRUE;

        if (-90<=d1 && d1<=89)
            fDec[d1+90] = kTRUE;
    }

    //
    // count degrees of right ascension
    //
    fRaCnt=0;
    for (int x=0; x<360; x++)
        if (fRa0[x])
            fRaCnt++;
    //cout << "fRaCnt: " << setw(3) << fRaCnt << "  " << flush;

    //
    // calculate min and max of declination
    //
    for (int y=0; y<180; y++)
    {
        if (fDec[y])
            fDecMax = y;

        if (fDec[179-y])
            fDecMin = 179-y;
    }

    fDecMin -= 90;
    fDecMax -= 90;

    //
    // check whether radec north- or south-pole is in the visible region
    //
    byte img[768*576];
    if (DrawRaDec(0, img, 0, 90))
    {
        fDecMax=89;
        cout << "Ra Dec Pole1 Inside!" << endl;
    }
    if (DrawRaDec(0, img, 0, -90))
    {
        fDecMin=-90;
        cout << "Ra Dec Pole1 Inside!" << endl;
    }

    //cout << "fDecMin: " << setw(3) << fDecMin << "  ";
    //cout << "fDecMax: " << setw(3) << fDecMax << endl;
}

void StarCatalog::DrawSCAltAz(byte *img, const int color) const
{
    //
    // ------------ draw az lines ---------------
    //
    for (int az=0; az<360; az++)
    {
        if (!fAz0[az])
            continue;

        for (double alt=fAltMin-1; alt<fAltMax+1; alt+=0.006*(fAltMax-fAltMin))
        {
            if ((alt>88 && az%5) || alt>89.5)
                continue;

            DrawAltAz(color, img, alt, az);
        }
    }

    //
    // ------------ draw alt lines ---------------
    //
    for (int alt=fAltMin; alt<fAltMax+1; alt++)
    {
        for (double az=0; az<360; az+=0.004*fAzCnt)
        {
            if (!fAz0[(int)az] && !fAz0[(int)(az+359)%360] && !fAz0[(int)(az+1)%360])
                continue;

            DrawAltAz(color, img, alt, az);
        }
    }
}

void StarCatalog::DrawSCRaDec(byte *img, const int color) const
{
    //
    // ------------ draw ra lines ---------------
    //
    for (int ra=0; ra<360; ra++)
    {
        if (!fRa0[ra])
            continue;

        for (double dec=fDecMin-1; dec<fDecMax+1; dec+=0.005*(fDecMax-fDecMin))
        {
            if ((dec>88 && ra%5) || dec>89.5)
                continue;

            DrawRaDec(color, img, ra, dec, ra==0||ra==90);
        }
    }

    //
    // ------------ draw dec lines ---------------
    //
    for (int dec=fDecMin; dec<fDecMax+1; dec++)
    {
        for (double ra=0; ra<360; ra+=0.003*fRaCnt)
        {
            if (!fRa0[(int)ra])
                continue;

            DrawRaDec(color, img, ra, dec, dec==89);
        }
    }
}

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

    DrawSCAltAz(cimg, 2<<4);
    DrawSCRaDec(cimg, 2);

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

void StarCatalog::GetImg(byte *img, byte *cimg, const double utc,
                         const RaDec &radec)
{
    MStarList list;
    GetStars(list, utc, radec);
    GetImg(img, cimg, list);
    /*
     // memset(img,  0, 768*576);
     SetMjd(utc);
     //fAlpha = sla.GetAlpha();
     SetRaDec(radec);
     //CalcImg(cimg);
     */
}

void StarCatalog::GetImg(byte *img, byte *cimg, const double utc,
                         const AltAz &altaz)
{
    MStarList list;
    GetStars(list, utc, altaz);
    GetImg(img, cimg, list);
    /*
     // memset(img,  0, 768*576);

     SetMjd(utc);
     //fAlpha = sla.GetAlpha();
     SetAltAz(altaz);

     CalcRaDecRange();

     //CalcImg(img);
     */

}

void StarCatalog::GetStars(MStarList &list, const double utc, const RaDec &radec)
{
    SetMjd(utc);
    SetRaDec(radec);

    CalcStars(list);
}

void StarCatalog::GetStars(MStarList &list, const double utc, const AltAz &altaz)
{
    SetMjd(utc);
    SetAltAz(altaz);

    CalcRaDecRange();
    CalcStars(list);
}

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;
    }
}

Bool_t StarCatalog::DrawAltAz(const int color, byte *img, double alt, double az, int size) const
{
    //
    // alt/az[deg] -> alt/az[rad]
    //
    alt *= D2PI/360.0;
    az  *= D2PI/360.0;

    //
    // alt/az[rad] -> alt/az[pix]
    //
    double dx, dy;
    slaDe2h(az-fAltAz.Az(), -alt, DPI/2-fAltAz.Alt(), &dx, &dy);

    //
    // Align alt/az[pix]
    //
    const int xx = (int)(((dx-DPI)*fCosAngle - dy*fSinAngle + fWidth)/fPixSize);
    const int yy = (int)(((dx-DPI)*fSinAngle + dy*fCosAngle + fHeight)/fPixSize);
    //const int xx = 767-(int)((fWidth-dx+DPI)/fPixSize);
    //const int yy =     (int)((fHeight+dy)/fPixSize);

    //
    // Range Check
    //
    if (!(0<=xx && xx<768 && 0<=yy && yy<576))
        return kFALSE;

    //
    // Draw
    //
    DrawCircle(color, img, xx, yy, size);

    return kTRUE;
}

Bool_t StarCatalog::Draw(const int color, byte *img, const AltAz &altaz)
{
    return DrawAltAz(color, img, altaz.Alt(), altaz.Az());
}

/*
Bool_t StarCatalog::Draw(const int color, byte *img, const SaoFile *sao)
{
    if (sao->MagV() > fLimitMag)
        return kFALSE;

    //
    // ---- mean to observed ---
    //
    AltAz altaz=CalcAltAz(sao->GetRaDec()) * 360.0/D2PI;

    const int mag = (10 - (sao->MagV()>1 ? (int)sao->MagV() : 1))/2;

    //
    // ---- imaging -----
    //
    return DrawAltAz(color, img, altaz.Alt(), altaz.Az(), mag);
}
*/

Bool_t StarCatalog::DrawRaDec(const int color, byte *img, double ra, double dec, int size) const
{
    //
    // radec[deg] -> radec[rad]
    //
    ra  *= D2PI/360.0;
    dec *= D2PI/360.0;

    //
    // radec[rad] -> hadec[rad]
    //
    const double ha = GetAlpha()-ra;

    //
    // hadec[rad] -> altaz[rad]
    //
    double alt, az;
    slaDe2h(ha, dec, GetPhi(), &az, &alt);

    //
    // altaz[rad] -> altaz[deg]
    //
    alt *= 360.0/D2PI;
    az  *= 360.0/D2PI;

    return DrawAltAz(color, img, alt, az, size);
}

Bool_t StarCatalog::Draw(const int color, byte *img, const RaDec &radec)
{
    return DrawRaDec(color, img, radec.Ra(), radec.Dec());
}
/*
void StarCatalog::CalcImg(byte *img)
{

    //
    // --------- search for stars in catalog ----------
    //
    int count   = 0;
    int deleted = 0;

    int idx     = 0;

    while (fSrt[idx].dec<fDecMin)
        idx++;

    idx--;
    while (++idx<fEntries && fSrt[idx].dec<fDecMax+1)
    {
        const int ra = fSrt[idx].ra;

        if (!fRa0[ra])
            continue;

        int nr = fSrt[idx].nr;
        do
        {
            //
            // Get entry from catalog
            //
            fSao->GetEntry(nr++);

            //
            // Try to draw star into the image
            //  white = 0xff
            //
            if (!Draw(0x0f, img, fSao))
                deleted++;

            count++;
        }
        while ((int)(360.0/D2PI*fSao->Ra())==ra);
    }

    cout << " " << count << "-" << deleted << "=" << count-deleted << " " << flush;
}
*/
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 = 0x0f;

        DrawCircle(color, img, (int)star->GetX(), (int)star->GetY(), mag);
    }
}

void StarCatalog::CalcStars(MStarList &list) const
{
    //
    // --------- search for stars in catalog ----------
    //
    if (fEntries==0)
        return;

    int count   = 0;
    int deleted = 0;

    int idx     = 0;

    while (fSrt[idx].dec<fDecMin)
        idx++;

    idx--;
    while (++idx<fEntries && fSrt[idx].dec<fDecMax+1)
    {
        const int ra = fSrt[idx].ra;

        if (!fRa0[ra])
            continue;

        int nr = fSrt[idx].nr;
        do
        {
            //
            // Get entry from catalog
            //
            fSao->GetEntry(nr++);

            if (fSao->MagV() > fLimitMag)
                continue;

            //
            // ---- mean to observed ---
            //
            AltAz altaz=CalcAltAz(fSao->GetRaDec());

            //
            // alt/az[rad] -> alt/az[pix]
            //
            double dx, dy;
            slaDe2h(altaz.Az()-fAltAz.Az(), -altaz.Alt(),
                    DPI/2-fAltAz.Alt(), &dx, &dy);

            //
            // Align and rotate alt/az[pix]
            //
            float xx = ((dx-DPI)*fCosAngle - dy*fSinAngle + fWidth)/fPixSize;
            float yy = ((dx-DPI)*fSinAngle + dy*fCosAngle + fHeight)/fPixSize;

            //
            // Range Check, add stars to the list
            //
            if (!(0<=xx && xx<768 && 0<=yy && yy<576))
            {
                deleted++;
                continue;
            }

            list.Add(xx, yy, fSao->MagV());
            count++;
        }
        while ((int)(360.0/D2PI*fSao->Ra())==ra);
    }

    cout << "Showing " << count+deleted << "-" << deleted << "=" << count << " stars." << endl;
}

AltAz StarCatalog::CalcAltAzFromPix(Double_t pixx, Double_t pixy) const
{
    pixx *= fPixSize;
    pixy *= fPixSize;

    const double dx =  (pixx-fWidth)*fCosAngle + (pixy-fHeight)*fSinAngle;
    const double dy = -(pixx-fWidth)*fSinAngle + (pixy-fHeight)*fCosAngle;

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