не совсем понял, или, если правильнее расставить слова, ...... Сходил по ссылкам, там файлы с разными расширениями, посмотрел как мне показалось все, - не увидел файла где есть на каждый спутник те переменные, что я использую в вычислениях. --- Сообщения объединены, 11 фев 2019, Оригинальное время сообщения: 11 фев 2019 --- Вот у меня такие результаты: файл SP3. 2019 2 11 2 15 PG10 -4236.001426 19163.988727 17856.678884 PG20 -11704.584987 21701.419721 9579.911955 * 2019 2 11 2 30 PG10 -5824.409922 20286.280184 16057.525893 PG20 -12607.680880 22183.682643 7001.452871 Мои данные время местное, примерно 9-30 Спутник 10 -4326, хххх 19910,хххххх ............... Спутник 20 -11501,хххх 22331,ххххх ................... дробные значения не печатаю здесь, но они есть. вот код: /////////////////////////////вызывается при совпадении наличия псевдодальности и эфемерид void calk_SV(int x) // расчет положения спутника по его эфемеридам { double Ttr; a_bp = pow(GP[x].square_big_axis,2); // вычмсляем большую полуось n_o = sqrt(gp_earth/(a_bp * a_bp * a_bp)); // n = n_o + GP[x].delta_n; // находим среднее движение Ttr = t - PR_sv[x] / c; // PR_sv[x] – псевдодальность вычисляемого спутника t_k = Ttr - GP[x].t_oe; // вычисляем промежуток времени от опорной эпохи эфемерид if ( t_k > 302400) // если t больше 302400 секунд, вычтите 604800 секунд из T. t_k = t_k - 604800; else if ( t_k < -302400) // Если T меньше -302400 секунд, добавьте 604800 секунд к T. t_k = t_k + 604800; M_sr = GP[x].Mo + (n * t_k); // определяем среднюю аномалию E1 = M_sr; // начальное условие для вычислений for (;;) // находим эксцентрическую аномалию через уравнение Кеплера { E1 = M_sr + (e * sin(GP[x].E_a)); // Kepler_equation(x); if(fabs(E1 - GP[x].E_a) > 1e-9) GP[x].E_a = E1; // E_a - это бывшая E(кампилятор не различает Е и е) else break; } v = (atan((sqrt((1 + e)/(1 - e))) * tan(GP[x].E_a/2)))*2; nash = v + GP[x].omega; // находим невозмущенный аргумент широты sigma_u = GP[x].Cus * sin(2 * nash) + GP[x].Cuc * cos(2 * nash); sigma_r = GP[x].Crs * sin(2 * nash) + GP[x].Crc * cos(2 * nash); sigma_i = GP[x].Cis * sin(2 * nash) + GP[x].Cic * cos(2 * nash); u_eph = nash + sigma_u; r_eph = a_bp * (1 - e * cos(GP[x].E_a)) + sigma_r; i_eph = GP[x].io + sigma_i + GP[x].IDOT * t_k; X1 = r_eph * cos(u_eph); Y1 = r_eph * sin(u_eph); Z1 = 0; OMEGA_r = GP[x].omega_0 + (GP[x].OMEGADOT - OMEGA_zemli)*t_k - OMEGA_zemli*GP[x].t_oe; X_wgs84 = X1 * cos(OMEGA_r)- Y1 * cos(i_eph) * sin(OMEGA_r); Y_wgs84 = X1 * sin(OMEGA_r)+ Y1 * cos(i_eph) * cos(OMEGA_r); Z_wgs84 = Y1 * sin(i_eph); //// временный кусочек кода tft_fill_screen(); // очистить экран out_b (x, 4 ,5); // отображение номера спутника для которого будут вычислены параметры tftout_shislo_1 (X_wgs84/1000); tftout_shislo_2 (Y_wgs84/1000); Delay_ms(60); } ///////////////////////////////////////////////////
спасибо. постараюсь воспользоваться........ но может вы по коду скажете что то умное? там правильно все..... или там не правильно все...... или еще варианты есть?
ftp://cddis.gsfc.nasa.gov/gps/data/daily/2019/042/19n/brdc0420.19n.Z описание Код: +----------------------------------------------------------------------------+ | TABLE A3 | | GPS NAVIGATION MESSAGE FILE - HEADER SECTION DESCRIPTION | +--------------------+------------------------------------------+------------+ | HEADER LABEL | DESCRIPTION | FORMAT | | (Columns 61-80) | | | +--------------------+------------------------------------------+------------+ |RINEX VERSION / TYPE| - Format version (2.10) | F9.2,11X, | | | - File type ('N' for Navigation data) | A1,19X | +--------------------+------------------------------------------+------------+ |PGM / RUN BY / DATE | - Name of program creating current file | A20, | | | - Name of agency creating current file | A20, | | | - Date of file creation | A20 | +--------------------+------------------------------------------+------------+ *|COMMENT | Comment line(s) | A60 |* +--------------------+------------------------------------------+------------+ *|ION ALPHA | Ionosphere parameters A0-A3 of almanac | 2X,4D12.4 |* | | (page 18 of subframe 4) | | +--------------------+------------------------------------------+------------+ *|ION BETA | Ionosphere parameters B0-B3 of almanac | 2X,4D12.4 |* +--------------------+------------------------------------------+------------+ *|DELTA-UTC: A0,A1,T,W| Almanac parameters to compute time in UTC| 3X,2D19.12,|* | | (page 18 of subframe 4) | 2I9 | | | A0,A1: terms of polynomial | | | | T : reference time for UTC data | *) | | | W : UTC reference week number. | | | | Continuous number, not mod(1024)! | | +--------------------+------------------------------------------+------------+ *|LEAP SECONDS | Delta time due to leap seconds | I6 |* +--------------------+------------------------------------------+------------+ |END OF HEADER | Last record in the header section. | 60X | +--------------------+------------------------------------------+------------+ Records marked with * are optional +----------------------------------------------------------------------------+ | TABLE A4 | | GPS NAVIGATION MESSAGE FILE - DATA RECORD DESCRIPTION | +--------------------+------------------------------------------+------------+ | OBS. RECORD | DESCRIPTION | FORMAT | +--------------------+------------------------------------------+------------+ |PRN / EPOCH / SV CLK| - Satellite PRN number | I2, | | | - Epoch: Toc - Time of Clock | | | | year (2 digits, padded with 0 | | | | if necessary) | 1X,I2.2, | | | month | 1X,I2, | | | day | 1X,I2, | | | hour | 1X,I2, | | | minute | 1X,I2, | | | second | F5.1, | | | - SV clock bias (seconds) | 3D19.12 | | | - SV clock drift (sec/sec) | | | | - SV clock drift rate (sec/sec2) | *) | +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 1| - IODE Issue of Data, Ephemeris | 3X,4D19.12 | | | - Crs (meters) | | | | - Delta n (radians/sec) | | | | - M0 (radians) | | +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 2| - Cuc (radians) | 3X,4D19.12 | | | - e Eccentricity | | | | - Cus (radians) | | | | - sqrt(A) (sqrt(m)) | | +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 3| - Toe Time of Ephemeris | 3X,4D19.12 | | | (sec of GPS week) | | | | - Cic (radians) | | | | - OMEGA (radians) | | | | - CIS (radians) | | +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 4| - i0 (radians) | 3X,4D19.12 | | | - Crc (meters) | | | | - omega (radians) | | | | - OMEGA DOT (radians/sec) | | +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 5| - IDOT (radians/sec) | 3X,4D19.12 | | | - Codes on L2 channel | | | | - GPS Week # (to go with TOE) | | | | Continuous number, not mod(1024)! | | | | - L2 P data flag | | +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 6| - SV accuracy (meters) | 3X,4D19.12 | | | - SV health (bits 17-22 w 3 sf 1) | | | | - TGD (seconds) | | | | - IODC Issue of Data, Clock | | +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 7| - Transmission time of message **) | 3X,4D19.12 | | | (sec of GPS week, derived e.g. | | | | from Z-count in Hand Over Word (HOW) | | | | - Fit interval (hours) | | | | (see ICD-GPS-200, 20.3.4.4) | | | | Zero if not known | | | | - spare | | | | - spare | | +--------------------+------------------------------------------+------------+ **) Adjust the Transmission time of message by -604800 to refer to the reported week, if necessary. *) In order to account for the various compilers, E,e,D, and d are allowed letters between the fraction and exponent of all floating point numbers in the navigation message files. Zero-padded two-digit exponents are required, however. t – глобальная переменная? Вычисление координат "верха" лучше не завязывать на псевдодальность. Т.е. программа вычисления координат спутника (любого) должна зависеть только от времени (по любой разумной шкале) или его аналога. --- Сообщения объединены, 11 фев 2019, Оригинальное время сообщения: 11 фев 2019 --- Так как сутки ещё не закончились, то навигационный файл периодически обновляется. --- Сообщения объединены, 11 фев 2019 --- В этих делах принято время давать "по гринвичу", местное время нафиг никому не нужно, если только вы не занимаетесь определением долгот. Да и само понятие "местное время" требует определённых уточнений. Чтобы не лезть в дебри, время лучше давать по шкале GPS или всемирному координированному времени UTC.
описание что то не понимается....... обведите в этом файле что нибудь. например Crs , что бы начал понимать , где искать, и как понимать формат.... --- Сообщения объединены, 11 фев 2019, Оригинальное время сообщения: 11 фев 2019 --- новосибирск, разница во времени как я понимаю - 7 часов. у меня 9-00, там 2-00.
по поводу разницы "19 секунд", - я может недопонимаю. но как мне показалось если "время жизни" эфемерид 15 минут, то я, считая координаты спутника, должен получить значения которые лежат в файле на это время. то есть если посчитаю в 9-10, или в 9-12 то значения должны быть одни и те же, и совпасть с файлом.
Сейчас точность эфемерид в навигационном сообщении лучше 1.5 метров, причем основная ошибка по трансверсали. A Long-term Analysis of the GPS Broadcast Orbit and Clock Error Variations
это значит, что в целых числах мои вычисления уж точно должны бы совпадать?.............. с файлом..... Где ошибка?
http://leapsecond.com/java/gpsclock.htm Одна из лучших страничек, если не самая лучшая, по шкалам времени https://www.ucolick.org/~sla/leapsecs/timescales.html 19 секунд – это разница с TAI. Описание (столбец FORMAT) – это фортрановская нотация. Например, F9.2,11X, – число простой точности, ширина поля - 9 символов, 2 знака после запятой, затем 11 пробелов. A20 – символьное поле шириной 20 символов. Код: +--------------------+------------------------------------------+------------+ | BROADCAST ORBIT - 1| - IODE Issue of Data, Ephemeris | 3X,4D19.12 | | | - Crs (meters) | | | | - Delta n (radians/sec) | | | | - M0 (radians) | | +--------------------+------------------------------------------+------------+ 3 пробела, затем 4 числа double шириной поля 19 и 12 знаков после запятой. Остальные строчки – аналогично. P.S. Заодно можно проверить правильность декодирования фреймов. --- Сообщения объединены, 11 фев 2019, Оригинальное время сообщения: 11 фев 2019 --- Да кто ж его знает, входных значений я не вижу, кода тоже. Спросил откуда переменная t в операторе вы молчите.
я вот про это: 2 NAVIGATION DATA RINEX VERSION / TYPE CCRINEXN V1.6.0 UX CDDIS 11-FEB-19 02:26 PGM / RUN BY / DATE IGS BROADCAST EPHEMERIS FILE COMMENT 0.8382D-08 -0.7451D-08 -0.5960D-07 0.5960D-07 ION ALPHA 0.8806D+05 -0.1638D+05 -0.1966D+06 0.6554D+05 ION BETA -0.186264514923D-08 0.266453525910D-14 319488 2040 DELTA-UTC: A0,A1,T,W 18 LEAP SECONDS END OF HEADER 1 19 2 11 0 0 0.0- 0.160989351571D-03-0.693489710102D-11 0.000000000000D+00 0.900000000000D+02 0.141250000000D+02 0.447661504041D-08-0.488271654406D+00 0.700354576111D-06 0.837390543893D-02 0.587105751038D-05 0.515365682030D+04 0.864000000000D+05 0.633299350738D-07-0.109884769581D+01 0.147148966789D-06 0.974178706545D+00 0.271656250000D+03 0.684905223726D+00-0.828427364468D-08 0.154649298907D-09 0.100000000000D+01 0.204000000000D+04 0.000000000000D+00 0.280000000000D+01 0.000000000000D+00 0.558793544769D-08 0.900000000000D+02 0.813480000000D+05 0.400000000000D+01 0.000000000000D+00 0.000000000000D+00 2 19 2 11 0 0 0.0-0.136601272970D-03-0.989075488178D-11 0.000000000000D+00 где в этих столбиках цифр найти нужную мне Crs или остальные......... --- Сообщения объединены, 11 фев 2019, Оригинальное время сообщения: 11 фев 2019 --- Время для расчета у меня берется из сообщения от навигатора, B5 62 02 10 - заголовок 40 01 - размер данных C0 7F 77 07 - Время недели измерения по локальному времени приемника. rcvTow \ в мс C0 7F 77 07 это и есть t t_oe - из файла эфемерид спутника.
Дано: сверхбыстрые эфемериды в формате SP3 * 2019 2 11 0 30 0.00000000 PG01 13692.480507 -20362.626145 9584.830824 -161.000753 2 4 8 235 P P PG02 -15261.494016 -899.262822 -21105.503768 -136.621949 8 11 8 225 P P PG03 19492.344528 -12952.351276 -12524.272973 182.333328 7 9 11 233 P P PG05 -25925.158766 1538.517184 -6171.278307 1.062307 11 6 8 142 P P PG06 -9939.063604 -13936.841330 -20259.450653 278.416006 6 9 6 250 P P PG07 5561.383441 -25386.833729 4331.836586 34.910982 6 5 11 234 P P PG08 16746.555142 -2583.630590 20549.144730 -132.280770 8 6 9 307 P P PG09 -2505.284424 -20424.154311 -16836.497751 460.255664 11 9 10 185 P P PG10 10885.726916 12521.137570 20833.643221 119.208861 12 5 10 204 P P навигационное сообщение для 10 спутника 10 19 2 11 0 0 0.0 0.119221396744D-03-0.738964445190D-11 0.000000000000D+00 0.530000000000D+02-0.560312500000D+02 0.448268672189D-08-0.251927565154D+01 -0.293180346489D-05 0.432240834925D-02 0.723637640476D-05 0.515366408348D+04 0.864000000000D+05-0.931322574616D-08-0.596143988819D-01-0.484287738800D-07 0.962769973131D+00 0.242750000000D+03-0.275661561130D+01-0.805533553706D-08 -0.129291099779D-09 0.100000000000D+01 0.204000000000D+04 0.000000000000D+00 0.200000000000D+01 0.000000000000D+00 0.186264514923D-08 0.530000000000D+02 0.792180000000D+05 0.400000000000D+01 0.000000000000D+00 0.000000000000D+00 грязный код Код: #include <stdio.h> #include "float.h" #include "math.h" const double mu = 3.986005E+14; // WGS 84 value of earth's univ. grav. par. const double We = 7.292115E-5; // WGS-84 earth rotation rate const double Wedot = 7.2921151467E-5; // WGS 84 value of earth's rotation rate const double pi = 4*atan(1.0); const double tol = 0.1*sqrt(DBL_EPSILON); int main(int argc, char **argv) { /* int prn = 10; int Year = 2019; int Month = 2; int Day = 11; int hour = 0; int minute = 0; double sec = 0; */ double Crs = -0.560312500000E+02; double dn = +0.448268672189E-08; double M0 = -0.251927565154E+01; double Cuc = -0.293180346489E-05; double ec = 0.432240834925E-02; double Cus = 0.723637640476E-05; double a = 0.515366408348E+04*0.515366408348E+04; double Toe = 0.864000000000E+05; double Cic = -0.931322574616E-08; double W0 = -0.596143988819E-01; double Cis = -0.484287738800E-07; double i0 = 0.962769973131E+00; double Crc = 0.242750000000E+03; double w = -0.275661561130E+01; double Wdot= -0.805533553706E-08; double idot= -0.129291099779E-09; double beta; double nu; double n0 = sqrt(mu/(a*a*a)); double n = n0 + dn; double T = 1800; double M = M0 + n*T; double E = M; // start value for E double deltaE, sinE, cosE, FK, dFKdE; int iter = 0; printf("\ntol = %10.2e",tol); do { sinE = sin(E); cosE = cos(E); // sincos(E, *sinE, *cosE); FK = E - M - ec*sinE; dFKdE = 1.0 - ec*cosE; deltaE = -FK/dFKdE; printf("\niter = %3d FKepler = %12.2e deltaE = %12.1e\n", ++iter, FK, deltaE); E += deltaE; } while (fabs(deltaE) > tol); sinE = sin(E); cosE = cos(E); beta = (1.0 + sqrt(1.0 + ec*ec))/ec; nu = atan(sinE/(beta - cosE)); nu += nu + E; double phi, du, dr, di, u, r, i, Xdash, Ydash, Wc, X, Y, Z; phi = nu + w; du = Cuc*cos(2*phi) + Cus*sin(2*phi); dr = Crc*cos(2*phi) + Crs*sin(2*phi); di = Cic*cos(2*phi) + Cis*sin(2*phi); u = phi + du; r = a*(1 - ec*cos(E)) + dr; i = i0 + idot*T + di; Xdash = r*cos(u); Ydash = r*sin(u); Wc = W0 + (Wdot - Wedot)*T - Wedot*Toe; X = Xdash*cos(Wc) - Ydash*cos(i)*sin(Wc); Y = Xdash*sin(Wc) + Ydash*cos(i)*cos(Wc); Z = Ydash*sin(i); printf("\nPG10 10885726.916 12521137.570 20833643.221"); printf("\nX = %15.3f Y = %15.3f Z = %15.3f\n\n", X, Y, Z); system("PAUSE"); return 0; } результат Код: PG10 10885726.916 12521137.570 20833643.221 X = 10885725.043 Y = 12521137.636 Z = 20833642.775
Добрый день. Спасибо. Посмотрел код, Посмотрел исходные данные, "разобрался" с файлом. 10 19 2 11 0 0 0.0 0.119221396744D-03 -0.738964445190D-11 0.000000000000D+00 0.530000000000D+02 --------- Crs ----------- ------------- dn------- ------------M0-------- --------- Cuc-------- ---------- ec-------- ------------Cus------- -------------a-------- -----------Toe ------------ ---------- Cic-------- ------------- W0------- ----------- Cis------- ---------- i0-------- -------- Crc--------- -------------- w-------- ---------- Wdot------ -----------idot --------- 0.100000000000D+01 0.204000000000D+04 0.000000000000D+00 0.200000000000D+01 0.000000000000D+00 0.186264514923D-08 0.530000000000D+02 0.792180000000D+05 0.400000000000D+01 0.000000000000D+00 0.000000000000D+00 сравнил значения, и понял что часть переменных у меня все таки не совпадают (декодируются неправильно) сейчас пытаюсь разобраться почему, но кажется вот увидел что то: В документе написано что то про "дополнение до двух", - я этого не делал. скорее всего - ошибка тут.
Добрый день. Кажется разобрался со своей ошибкой декодирования. Не внимательно читал, мои переменные в кадрах передаются в дополнительном коде, поэтому все положительные у меня преобразовывались правильно, - а отрицательные неверно. сейчас, надеюсь, я это исправил. Начал разбираться с Вашим кодом (спасибо за него), вот вопросы появились: вот заголовок файла, - и дата его: 12 февраля 10 спутник время 00 10 19 2 12 0 0 0.0 считаю время, которое Вы поставили, - чтобы получить 1800, пришлось выставить время 30 минут. почему так? и второе - пробую посчитать время для файла с заголовком 12 февраля 10 спутник время 6 часов 00 10 19 2 12 6 0 0.0 считаю, и по этому времени (23400) пробую считать координаты. Не получается правильно.
http://geodesist.ru/resources/gps-date-time-ver-1-0-0-0.19/ 30m – это разность момента на который нам даны "контрольные" координаты в формате SP3 и ближайшим моментом, на который известны "кеплеровы элементы орбиты" и их изменения (эфемериды) из навигационного сообщения (в нашем примере – из файла *.19n).
