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

#include "StarCatalog.h"

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

#include <TSystem.h>

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

StarCatalog::StarCatalog() : fEntries(0)
{
    // p = pointer to MainFrame (not owner)

    //
    // calculate observers location (goe)
    //
    int status;
    slaDaf2r(51, 38, 48.0, &fPhi,   &status);
    slaDaf2r( 9, 56, 36.0, &fElong, &status);

    cout << "fPhi:   51ø38'48.0\" = " <<  360.0/D2PI*fPhi << endl;
    cout << "fElong:  9ø56'36.0\" = " <<  360.0/D2PI*fElong << endl;

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

    if (!idx)
        exit(0);

    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()
{
    delete fSrt;
    delete fSao;
}

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

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

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 = AltAz2RaDec(fAltAz);

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

    CalcAltAzRange();
    CalcRaDecRange();
}

RaDec StarCatalog::AltAz2RaDec(const AltAz &altaz) const
{
    //
    // -- observed to apparent --
    //  Workaraound for slalib: discard const
    //
    double r=0, d=0;
    slaOapqk ("A", altaz.Az(), DPI/2-altaz.Alt(), (double*)fAoprms, &r, &d);

    //
    // -- apparent to mean --
    //  Workaraound for slalib: discard const
    //
    double ra, dec;
    slaAmpqk(r, d, (double*)fAmprms, &ra, &dec);

    return RaDec(ra, dec);
}

AltAz StarCatalog::RaDec2AltAz(const RaDec &radec, const RaDec &rdpm) const
{
    //
    // ---- Mean to apparent ----
    //

    double r=0, d=0;
    slaMapqkz(radec.Ra(), radec.Dec(), (double*)fAmprms, &r, &d);
    //
    // Doesn't work - don't know why
    //
    //    slaMapqk (radec.Ra(), radec.Dec(), rdpm.Ra(), rdpm.Dec(),
    //              0, 0, (double*)fAmprms, &r, &d);
    //

    //
    // -- apparent to observed --
    //
    double r1=0;  // ra
    double d1=0;  // dec
    double h0=0;  // ha

    double zd;
    double az;
    slaAopqk (r, d, (double*)fAoprms,
              &az,    // observed azimuth (radians: N=0,E=90)
              &zd,    // observed zenith distance (radians) [-pi/2, pi/2]
              &h0,    // observed hour angle (radians)
              &d1,    // observed declination (radians)
              &r1);   // observed right ascension (radians)

    return AltAz(DPI/2-zd, az);
}

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

    RaDec pm = rdpm * D2PI/360.0;

    fAltAz = RaDec2AltAz(fRaDec, pm);

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

    CalcRaDecRange();
    CalcAltAzRange();
}

void StarCatalog::CalcAltAzRange()
{
    char 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
    //
    char 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
    //
    char 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;

    //
    // 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, fPhi, &ha0, &de0);

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

        const int h0 = ((int)(360.0/D2PI*(fAlpha-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*(fAlpha-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, fPhi, &ha0, &de0);

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

        const int h0 = ((int)(360.0/D2PI*(fAlpha-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*(fAlpha-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
    //
    char 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(char *img, const int color)
{
    //
    // ------------ 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(char *img, const int color)
{
    //
    // ------------ 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);
        }
    }

//    DrawRaDec(color, img, 37.953, 89.2641, 10);
}


void StarCatalog::SetMjd(const double mjd)
{
    fMjd   = mjd;
    fAlpha = slaGmst(mjd) + fElong;

    cout << "UTC: " << mjd << endl;

    //
    // ----- calculate star independent parameters ----------
    //
    slaMappa(2000.0, fMjd, fAmprms);
    slaAoppa(fMjd, 0,                // mjd, UT1-UTC
             fElong, fPhi, 148,      // g”ttingen long, lat, height
             0, 0,                   // polar motion x, y-coordinate (radians)
             273.155, 1013.25, 0.5,  // temp, pressure, humidity
             0.2, 0.0065,            // wavelength, tropo lapse rate
             fAoprms);

}

void StarCatalog::DrawCross(char *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(char *img, char *cimg, const double utc,
                         const RaDec &radec)
{
    // memset(img,  0, 768*576);
    memset(cimg, 0, 768*576);

    SetMjd(utc);

    SetRaDec(radec, RaDec());

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

    CalcImg(img);

    DrawCross(img, 768/2, 576/2);
}

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

    SetMjd(utc);

    SetAltAz(altaz);

    CalcRaDecRange();
    DrawSCRaDec(cimg, 2);
    DrawSCAltAz(cimg, 2<<4);
    CalcImg(img);

    DrawCross(img, 768/2, 576/2);
}

void StarCatalog::DrawCircle(int color, char *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, char *img, double alt, double az, int size)
{
    //
    // 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 = 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, char *img, const AltAz &altaz)
{
    return DrawAltAz(color, img, altaz.Alt(), altaz.Az());
}

Bool_t StarCatalog::Draw(const int color, char *img, const SaoFile *sao)
{
    //
    // ---- mean to observed ---
    //
    AltAz altaz=RaDec2AltAz(sao->GetRaDec(), sao->GetRaDecPm()) * 360.0/D2PI;

    if (sao->MagV() > fLimitMag)
        return kFALSE;

    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, char *img, double ra, double dec, int size)
{
    //
    // radec[deg] -> radec[rad]
    //
    ra  *= D2PI/360.0;
    dec *= D2PI/360.0;

    //
    // radec[rad] -> hadec[rad]
    //
    const double ha = fAlpha-ra;

    //
    // hadec[rad] -> altaz[rad]
    //
    double alt, az;
    slaDe2h(ha, dec, fPhi, &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, char *img, const RaDec &radec)
{
    return DrawRaDec(color, img, radec.Ra(), radec.Dec());
}

void StarCatalog::CalcImg(char *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)
            //
            if (!Draw(0xff, img, fSao))
                deleted++;

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

    cout << " " << count << "-" << deleted << "=" << count-deleted << " " << flush;
}