Скрытый Ra, De c для Alt, Az - результаты немного - PullRequest
0 голосов
/ 22 марта 2020

Я пишу код для преобразования значений прямого склонения (Ra) (De c) в высоту (Alt) и азимут (Az) на основе этого веб-сайта: http://www.stargazing.net/kepler/altaz.html Я заинтересован на текущей позиции астрообъекта. (Не интересно, где это было 25 лет go). Код написан на c ++ (все еще выполняется):

#include <time.h>
#include <iostream>
#include <chrono>


#define _USE_MATH_DEFINES
#include <math.h>

void double_to_nice_out(double value) {
    double degree; double min; double sec;
    min = modf(value, &degree);
    min = min * 60;
    sec = modf(min, &min);
    sec = sec * 60;

    printf("%d%c %d' %f\"\t(%f)\n", (int)degree, 248, (int)min, sec, value);
}

void degree_to_rad(double *_value) {
    *_value = *_value * M_PI / 180;
}

void rad_to_degree(double* _value) {
    *_value = *_value * 180 / M_PI;
}

void conv_ra(bool is_pos, int _rah, int _ram, double _ras, double *_RA ) {
    *_RA = ( (double)_rah + ((double)_ram / 60) + (_ras / 3600) ) * 15; //*15 because hour to degree conversion
    if (!is_pos) {
        *_RA = (-1) * (*_RA);
    }
}
void conv_dec(bool is_pos, int _decd, int _decm, double _decs, double *_DEC) {
    *_DEC = (double)_decd + ((double)_decm / 60) + (_decs / 3600);
    if (!is_pos) {
        *_DEC = (-1) * (*_DEC);
    }
}
void conv_lat(int _latd, double _latm, std::string _latcd, double *_LAT) {
    //?south negatív?
    *_LAT = _latd + _latm / 60;
    if (_latcd == "S") {
        *_LAT = *_LAT * (-1);
    }
    else if (_latcd == "N") {
        //do nothing
    }
    else if (_latcd == "") {
        //do nothing
    }
    else {
        *_LAT = NULL;
    }
}
void conv_long(int _longd, double _longm, std::string _longcd, double* _LONG) {
    // is south negativ?
    *_LONG = _longd + _longm / 60;
    if (_longcd == "W") {
        *_LONG = *_LONG * (-1);
    }
    else if (_longcd == "E") {
        //do nothing
    }
    else if (_longcd == "") {
        //do nothing
    }
    else {
        *_LONG = NULL;
    }
}

void datumdiff(int _year, int _mon, int _mday, int _hour, int _min, int _sec, int _isdst, bool _isutc, double *_datediff) {
    struct tm t;
    time_t t_of_day;

    t.tm_year = _year - 1900;
    t.tm_mon = _mon;
    t.tm_mday = _mday;          // Day of the month
    t.tm_hour = _hour;     //UTC: add +2h, if you live in Hungary as I do
    t.tm_min = _min;
    t.tm_sec = _sec;
    t.tm_isdst = _isdst;        // Is DST on? 1 = yes, 0 = no, -1 = unknown
    t_of_day = mktime(&t);

    long J2000 = 946728000; //J2000: 2000.01.01 12:00 UT // interesting that it is not midnight


    if (_isutc == 1) {
        *_datediff = (((double)t_of_day + 7200) - (double)J2000) / 60 / 60 / 24;
    }
    else
    {
        *_datediff = (((double)t_of_day) - (double)J2000) / 60 / 60 / 24;
    }
}

void calc_lst(double _datediff, int _hour, int _min, double _LONG, double *_lst) {
    double ut = (double)_hour + (double)_min / 60;
    *_lst = 100.46 + 0.985647 * _datediff + _LONG + 15 * ut;
    //printf("_lst: %f\n", *_lst);
    for (int i = 0; i < 100; i++) { //this is clearly not ok but works for now
        if (*_lst < 0) {
            *_lst = *_lst + 360;
            //printf("_lst: %f\n", *_lst);
        }
        else if (*_lst > 360) {
            *_lst = *_lst - 360;
            //printf("_lst: %f\n", *_lst);
        }
        else {
            //do nothing
        }
    }
}

void calc_lst_now(double _datediff, int _hour, int _min, int _sec, int _millis, double _LONG, double* _lst) {
    //_sec = _sec - 45; //add leapsec
    double ut = (double)_hour + (double)_min / 60 + (double)_sec / 3600 + (double)_millis / 3600000;
    printf("Difference between [Date - J2000]: %f [days]\n", _datediff);

    *_lst = 100.46 + 0.985647 * _datediff + _LONG + 15 * ut;
    for (int i = 0; i < 100; i++) { //this is clearly not ok but works for now
        if (*_lst < 0) {
            *_lst = *_lst + 360;
            //printf("_lst: %f\n", *_lst);
        }
        else if (*_lst > 360) {
            *_lst = *_lst - 360;
            //printf("_lst: %f\n", *_lst);
        }
        else {
            //do nothing
        }
    }
}

void calc_hour_angle(double _lst, double _ra, double *_ha) {
    *_ha = _lst - _ra;
    if (*_ha < 0) {
        *_ha = *_ha + 360;
    }
}

void ha_dec_to_alt_az(double _ha, double _DEC, double _LAT, double* _alt, double* _az) {
    degree_to_rad(&_ha);
    degree_to_rad(&_DEC);
    degree_to_rad(&_LAT);

    *_alt = asin(sin(_DEC) * sin(_LAT) + cos(_DEC) * cos(_LAT) * cos(_ha));
    *_az = acos((sin(_DEC) - sin(*_alt) * sin(_LAT)) / (cos(*_alt) * cos(_LAT)));

    if (sin(_ha) < 0) {
        //do nothing
        rad_to_degree(_az);
    }
    else {
        rad_to_degree(_az);
        *_az = 360 - *_az;
    }

    rad_to_degree(_alt);
}

void radec_to_altaz_almanac(
    bool RA_is_positiv, int RAh, int RAm, double RAs,
    bool DEC_is_positiv, int DECd, int DECm, double DECs,
    int LATd, double LATm, std::string LATcd,
    int LONGd, double LONGm, std::string LONGcd,
    int year, int mon, int mday, int hour, int min, int sec, int isdst, bool isutc,
    double *alt, double *az
) {

    //sub result:
    double RA;          //[degree]
    double DEC;         //[degree]
    double LAT;         //[degree]
    double LONG;        //[degree]
    double datediff;    //[hour]
    double lst;         //[degree]
    double ha;          //[degree]

    conv_ra(RA_is_positiv, RAh, RAm, RAs, &RA);
    printf("RA: %f\n", RA );

    conv_dec(DEC_is_positiv, DECd, DECm, DECs, &DEC);
    printf("DEC: %f\n", DEC);

    conv_lat(LATd, LATm, LATcd, &LAT);
    printf("LAT: %f\n", LAT);

    conv_long(LONGd, LONGm, LONGcd, &LONG);
    printf("LONG: %f\n---------------------\n", LONG);

    datumdiff(year, mon, mday, hour, min, sec, isdst, isutc, &datediff);
    printf("Difference between [Date - J2000]: %f [days]\n", datediff);

    calc_lst(datediff, hour, min, LONG, &lst);
    printf("LST: %f\n", lst);

    calc_hour_angle(lst, RA, &ha);
    printf("Hour Angle: %f degrees\n", ha);

    ha_dec_to_alt_az(ha, DEC, LAT, alt, az);

}

void radec_to_altaz_now(
    bool RA_is_positiv, int RAh, int RAm, double RAs,
    bool DEC_is_positiv, int DECd, int DECm, double DECs,
    int LATd, double LATm, std::string LATcd,
    int LONGd, double LONGm, std::string LONGcd,
    unsigned __int64 now, int hour, int min, int sec, int millis,
    double* alt, double* az
) {

    //sub result:
    double RA;          //[degree]
    double DEC;         //[degree]
    double LAT;         //[degree]
    double LONG;        //[degree]
    double datediff;    //[hour]
    double lst;         //[degree]
    double ha;          //[degree]

    conv_ra(RA_is_positiv, RAh, RAm, RAs, &RA);
    printf("RA: %f\n", RA);

    conv_dec(DEC_is_positiv, DECd, DECm, DECs, &DEC);
    printf("DEC: %f\n", DEC);

    conv_lat(LATd, LATm, LATcd, &LAT);
    printf("LAT: %f\n", LAT);

    conv_long(LONGd, LONGm, LONGcd, &LONG);
    printf("LONG: %f\n---------------------\n", LONG);

    datediff = ((double)now / 1000 - 946728000 - 0) / 3600 / 24; //37 leap secunds until 2020 // 946728000 = J2000

    calc_lst_now(datediff, hour, min, sec, millis, LONG, &lst);
    printf("LST: %f\n", lst);

    calc_hour_angle(lst, RA, &ha);
    printf("Hour Angle: %f degrees\n", ha);

    ha_dec_to_alt_az(ha, DEC, LAT, alt, az);

} 


int main()
{
    //input data:
    bool RA_is_positiv = 1; //0 = false = negativ number
    int RAh = 6; //[hour]
    int RAm = 45; //[min]
    double RAs = 8.05; //[sec]

    bool DEC_is_positiv = 0; //0 = false = negativ number
    int DECd = 16; //[degree]
    int DECm = 43; //[min]
    double DECs = 23.9; //[sec]

    int LATd = 46; //[degree]
    double LATm = 19.11; //[min]
    std::string LATcd = "N";

    int LONGd = 17; //[degree]
    double LONGm = 45.09; //[min]
    std::string LONGcd = "E";

    int year = 1998;  // Year - 1900
    int mon = 7;      // Month, where 0 = jan
    int mday = 10;    // Day of the month
    int hour = 23;    //UTC: add +2h, if you live in Hungary as I do
    int min = 10;
    int sec = 0;
    int isdst = -1;   // Is DST on? 1 = yes, 0 = no, -1 = unknown
    bool isutc = 1;   // UTC: 1=yes, 0=no local time

    //################################ 
    //For current time, ra / dec values 
    struct tm newtime;
    __int64 ltime;
    _time64(&ltime);
    _gmtime64_s(&newtime, &ltime); // Obtain coordinated universal time hour, min, sec

    unsigned __int64 most = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
    int millis = most % 1000;       //millisecunds


    //final result:
    double alt, az;   //[degree] result

    printf("Almanac:\n");
    radec_to_altaz_almanac(RA_is_positiv, RAh, RAm, RAs, DEC_is_positiv, DECd, DECm, DECs, LATd, LATm, LATcd, LONGd, LONGm, LONGcd, year, mon, mday, hour, min, sec, isdst, isutc, &alt, &az);
    printf("AZ: ");
    double_to_nice_out(az);
    printf("ALT: ");
    double_to_nice_out(alt);

    printf("\n###########################\n\n");

    printf("NOW:\n");
    radec_to_altaz_now(RA_is_positiv, RAh, RAm, RAs, DEC_is_positiv, DECd, DECm, DECs, LATd, LATm, LATcd, LONGd, LONGm, LONGcd, most, newtime.tm_hour, newtime.tm_min, newtime.tm_sec, millis, &alt, &az);
    printf("AZ: ");
    double_to_nice_out(az);
    printf("ALT: "); 
    double_to_nice_out(alt);



    return 0;
} 

Кажется, что код работает, но результат немного отключен (Функция : radec_to_altaz_now). В случае Сириуса, около 0 ° 10 минут) по сравнению со Stellarium или Skyfield. Я знаю, что этот код не включает важные вещи, такие как дополнительные секунды (работа над ним) или скорость объекта и другие вещи. Интересно посмотреть, когда я вычитаю 50 секунд из текущего времени, мои результаты, кажется, становятся более точными. (Я думал, что из-за високосных секунд мне нужно добавить 37 se c ко времени [так как текущий год - 2020])

Вопрос: Как я могу сделать этот код более точным без особых усилий? Я не хочу делать из этого научное исследование, это мой хобби-проект. Но, возможно, я не учитываю то, что может быть легко закодировано. Или это отклонение не важно для хобби проекта? Как вы думаете, откуда исходит отклонение?

Спасибо! С уважением,

Роберт ПС: Я не знаю ни одного форума, где я мог бы публиковать вопросы об астрономическом кодировании.

...