расчет координат точки на местности используя исходные данные модуля навигатора

А зачем ему навигационный файл ГЛОНАСС, если алгоритм для GPS?

 

не совсем понял, или, если правильнее расставить слова, ...... Сходил по ссылкам, там файлы с разными расширениями, посмотрел как мне показалось все, - не увидел файла где есть на каждый спутник те переменные, что я использую в вычислениях.

--- Сообщения объединены, 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);
}
///////////////////////////////////////////////////

 

для кода есть специальный тэг - CODE

 

как этим пользоваться???

 

там кнопочки есть - для вставки смайликов, ссылок, картинок, а есть кнопочка "код"

 

спасибо. постараюсь воспользоваться........ но может вы по коду скажете что то умное? там правильно все..... или там не правильно все...... или еще варианты есть?

 

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.

 

разница 19 секунд?

 

описание что то не понимается....... обведите в этом файле что нибудь. например Crs , что бы начал понимать , где искать, и как понимать формат....

--- Сообщения объединены, 11 фев 2019, Оригинальное время сообщения: 11 фев 2019 ---

новосибирск, разница во времени как я понимаю - 7 часов. у меня 9-00, там 2-00.

 

по поводу разницы "19 секунд", - я может недопонимаю. но как мне показалось если "время жизни" эфемерид 15 минут, то я, считая координаты спутника, должен получить значения которые лежат в файле на это время. то есть если посчитаю в 9-10, или в 9-12 то значения должны быть одни и те же, и совпасть с файлом.

 

GPS живёт в своём времени

 

Сейчас точность эфемерид в навигационном сообщении лучше 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/

30^m – это разность момента на который нам даны "контрольные" координаты в формате SP3 и ближайшим моментом, на который известны "кеплеровы элементы орбиты" и их изменения (эфемериды) из навигационного сообщения (в нашем примере – из файла *.19n).