Функция gln_a
function [ ] = gln_a(A_PZ90_KM, ns, n0, ti_current, time_s0, alm_gln, timeUTC_leap, fid);
%Имя функции:gln_a
%Функция формирует альманах ГЛОНАСС в формате YUMA. Данные альманаха ГЛОНАСС соот- ветствуют
%данным альманаха ГЛОНАСС, записанным приемником СН 4701 и преобразованным в строгом со- ответствии
%интерфейсному контрольному документу ГЛОНАСС
%Постоянные:
I_MID = 1.0995574287564; %(double)(PI * 63.0 / 180.0). I_MID - Mean value of an inclination of a plane of orbit of a satellite
T_DR_MID = 43200.0; %Mean value a dragon of cycle time of a satellite
MU = 398600.44; % (Km^3/cek^2) гравитационная постоянная (constant of a gravitational)
OMEGA_Z = 0.7292115e-4; %(рад/cek)- Угловая скорость вращения Земли (Angular speed of rotation of the Earth )
SEC_IN_RAD = 7.2722052166430e-005; %(PI / 43200.0 ) Number radian in second of time
2*PI/(24*3600)=PI/43200
HALF_PI = pi / 2;
D7_3 = 7.0 / 3.0;
D7_4 = 7.0 / 4.0;
D7_6 = 7.0 / 6.0;
D7_24 = 7.0 / 24.0;
D49_72 = 49.0 / 72.0; nn = fix(ti_current / 86400); ti = ti_current - nn * 86400;%время в секундах i_incl = I_MID + alm_gln(ns).deltai; %наклонение орбиты спутника ГЛОНАСС с номером ns
156
ecc2_1 = 1.0 - alm_gln(ns).ecc * alm_gln(ns).ecc;
% t_dr = T_DR_MID + alm_gln(ns).Tdr; % GG24 t_dr = alm_gln(ns).Tdr; %драконический период спутника ГЛОНАСС (данные приемника СН 4701) n_dr = pi * 2 / t_dr;
%1. Вычисление полуоси a_n: a_n = semi_axis( A_PZ90_KM, t_dr, i_incl, alm_gln(ns).ecc, alm_gln(ns).omegan); alm_gln(ns).a = a_n;
%2. Вычисление t_lambda_k - время пересечения восходящего узла орбиты спутника ГЛОНАСС sin_i = sin(i_incl); cos_i = cos(i_incl); sin_i2 = sin_i * sin_i; cos_i2 = cos_i * cos_i; v = - alm_gln(ns).omegan; % omega_n - angle of a perigee
J = -0.00162393855 ; % J = 3/2 *C20; C20=1082.6257 * 10-6; ae_a = A_PZ90_KM / a_n; ae_a2 = ae_a * ae_a; j_ae_a2 = J * ae_a2; omega1 = j_ae_a2 * n_dr * cos_i / (ecc2_1 * ecc2_1); n0_na = (n0 - alm_gln(ns).Na) ; tz = ti - alm_gln(ns).TLambdaN + 86400.0 * n0_na; wk = tz / t_dr; wi = fix(wk); wi2 = wi * wi; t_lambda_kk = alm_gln(ns).TLambdaN + t_dr * wi + alm_gln(ns).dTdr * wi2; t_lambda_k = t_lambda_kk - n0_na * 86400.0 ; lambda_k = alm_gln(ns).LambdaN + (omega1 - OMEGA_Z) * (wi * t_dr + alm_gln(ns).dTdr * wi2) ;
%time_s0 - a true sidereal time to Greenwich midnight of date n0 to which time ti concerns
%time_s0 - истинное звездное время в текущий момент обсервации time_s0 = 0; time_s = time_s0 * SEC_IN_RAD + OMEGA_Z * t_lambda_k ; omega_0 = lambda_k + time_s;
% Auxiliary values: d1 = 1.0 - 1.5 * sin_i2; j_ae_a2_d1 = j_ae_a2 * d1; j_ae_a2_d = j_ae_a2 * (1.0 - 1.5 * sin_i); j_ae_a2_sin_i = j_ae_a2 * sin_i; j_ae_a2_sin_i2 = j_ae_a2 * sin_i2; j_ae_a2_cos_i2 = j_ae_a2 * cos_i2;
% 3. Calculation of constants: tau = 0;
157
l = cos(alm_gln(ns).omegan) * alm_gln(ns).ecc; h = alm_gln(ns).ecc * sin(alm_gln(ns).omegan); dop_x = (sqrt(1.0 - alm_gln(ns).ecc) * tan( v / 2.0)); dop_y = sqrt(1.0 + alm_gln(ns).ecc); ee = 2.0 * atan2(dop_x, dop_y ); ma = ee - alm_gln(ns).ecc * sin(ee); ll = ma + alm_gln(ns).omegan;
% for (j = 0; j <= 1; j++) for ( j= 1 : 2) sin_l = sin(ll); sin_2l = sin(2 * ll); sin_3l = sin(3 * ll); sin_4l = sin(4 * ll); cos_l = cos(ll); cos_2l = cos(2 * ll); cos_3l = cos(3 * ll); cos_4l = cos(4 * ll); l_cosl = l * cos_l; delta_a(j) = 2.0 * a_n * j_ae_a2_d1 * l_cosl +... j_ae_a2 * sin_i * (0.5 * h * sin_l - 0.5 * l_cosl +... cos(2.0 * alm_gln(ns).LambdaN) + 3.5 * l * cos_3l + 3.5 * h * sin_3l); delta_h(j) = j_ae_a2_d1 * (l * n_dr * tau + sin_3l +...
1.5 * l * sin_2l - 1.5 * h * cos_2l) -...
0.25 * j_ae_a2_sin_i2 * (sin_l - D7_3 * sin_3l +...
5.0 * l * sin_2l - 8.5 * l * sin_4l +...
8.5 * h * cos_4l + h * cos_2l) + j_ae_a2_cos_i2 *...
(tau * l * n_dr - 0.5 * l * sin_2l); delta_l(j) = j_ae_a2_d1 * ( - tau * h * n_dr + cos_l +...
1.5 * l * cos_2l + 1.5 * h * sin_2l) -...
0.25 * j_ae_a2_sin_i2 * ( - cos_l - D7_3 * cos_3l -...
5.0 * h * sin_2l - 8.5 * l * cos_4l -...
8.5 * h * sin_4l + l * cos_2l) + j_ae_a2_cos_i2 *...
( - tau * h * n_dr + 0.5 * h * sin_2l); d_omega(j) = j_ae_a2 * cos_i * (tau * n_dr + 3.5 * l * sin_l -...
3.5 * h * cos_l - 0.5 * sin_2l - D7_6 * sin_3l + ...
D7_6 * h * cos_3l); delta_i(j) = 0.5 * j_ae_a2_sin_i * cos_i * ( - l * cos_l +... h * sin_l + cos_2l + D7_3 * l * cos_3l + D7_3 * h * sin_3l); delta_ll(j) = 2.0 * j_ae_a2_d *... %
2.0 * j_ae_a2_d1 *
(tau * n_dr + D7_4 * l * sin_l -...
158
D7_4 * h * cos_l) + 3 * j_ae_a2_sin_i *...
( - D7_24 * h * cos_l - D7_24 * l * sin_l -...
D49_72 * h * cos_3l + D49_72 * l * sin_3l +...
0.25 * sin_2l) + j_ae_a2 * cos_i *...
(tau * n_dr + 3.5 * l * sin_l - 2.5 * h * cos_l -...
0.5 * sin_2l - D7_6 * l * sin_3l + D7_6 * h * cos_3l); tau = ti - t_lambda_k; ll = ma + alm_gln(ns).omegan + n_dr * tau; end;
% 4. Corrections to orbit elements of a satellite by second zone harmonic J20
% influence at the moments of time ti dda = delta_a(2) - delta_a(1); ddh = delta_h(2) - delta_h(1); ddl = delta_l(2) - delta_l(1); dd_omega = d_omega(2) - d_omega(1); ddi = delta_i(2) - delta_i(1); dd_ll = delta_ll(2) - delta_ll(1);
% 5. calculation of influenced elements of orbits at the moment of ti time ai = a_n + dda; hi = h + ddh; li = l + ddl;
% ecc_i = alm_gln[ns].ecc; ecc_i = sqrt(hi * hi + li * li); if (ecc_i == 0) w_i = 0; else if (li == ecc_i) w_i = HALF_PI; else if (li == - ecc_i) w_i = - HALF_PI; else if (li = 0) w_i = atan2(hi,li); else w_i = HALF_PI; end; end; end; end;
159
omega_i = omega_0 + dd_omega; ii_incl = i_incl + ddi; ll_z = ma + alm_gln(ns).omegan + n_dr * tau + dd_ll;
% ll_z = ma + alm_gln(ns).omegan + n_dr * (ti - t_lambda_k) + dd_ll; mm_i = ll_z - w_i;
%--------------------------------------------------
%timeUTC_leap
Na = alm_gln(ns).Na;
[wGPS,weekday] = week_GLONAS_gps(Na, timeUTC_leap);
Yuma(ns).week = wGPS;% - 1024 ;%варианты записи недели
Yuma(ns).id = alm_gln(ns).ID + 37;
%Запись альманаха ГЛОНАСС в формате YUMA
%Заголовок альманаха ГЛОНАСС. В заголовок вынесен номер (литера) частоты- Hn fprintf(fid,'GLO*Week %i almanac for prn-%i NAUHn-%i***** \n',Yuma(ns).week,
Yuma(ns).id,alm_gln(ns).Hn);
%Номер спутника ГЛОНАСС 1...24. В заголовке номера спутников ГЛОНАСС 38...61. fprintf(fid,'ID: %i\n',Yuma(ns).id-37);% для Planning
%fprintf(fid,'ID: %i\n',Yuma(ns).id);% вариант записи номера спутника
%Здоровье спутников ГЛОНАСС: 000- здоров
Yuma(ns).Health = alm_gln (ns).Health ; if (Yuma(ns).Health > 100) fprintf(fid,'Health: %i\n', Yuma(ns).Health); else if (Yuma(ns).Health > 10) fprintf(fid,'Health: 0%i\n', Yuma(ns).Health); else fprintf(fid,'Health: 00%i\n', Yuma(ns).Health); end; end;
%Ексцентриситет орбиты спутника strdop = e_norm( alm_gln(ns).ecc, 9); fprintf(fid,'Eccentricity: %s \n', strdop);
%Время, к которому относятся данные орбита спутника ГЛОНАСС, приведенное к началу недели
GPS.
%При необходимости нужно учесть поправку смещения времени GPS от UTC (в 2006 году- 14 се- кунд) strdop = e_norm((alm_gln(ns).TLambdaN + 86400 * (weekday)),9); fprintf(fid,'Time of Applicability(s): %s \n',strdop);
%Наклонение орбиты спутника ГЛОНАСС
Yuma(ns).inc = alm_gln(ns).deltai;
160
fprintf(fid,'Orbital Incluation(rad): %11.9e \n',Yuma(ns).inc);
%Скорость изменения восходящего узла орбиты спутника ГЛОНАСС strdop = e_norm(omega1,9); fprintf(fid,'Rate of Right Ascen(r/s): %s\n', strdop);
%Корень квадратный из большой полуоси орбиты спутника ГЛОНАСС
Yuma(ns).sqrt_a = sqrt(a_n * 1000); fprintf(fid,'SQRT(A) (m^1/2): %11.6f \n',Yuma(ns).sqrt_a);
%Долгота восходящего узла орбиты спутника ГЛОНАСС strdop = e_norm(omega_i,9); fprintf(fid,'Right Ascen at Week(rad): %s\n', strdop);
%Аргумент перигея орбиты спутника ГЛОНАСС strdop = e_norm(w_i,9); fprintf(fid,'Argument of Perigee(rad): %s\n',strdop);
%Средняя аномалия спутника ГЛОНАСС strdop = e_norm(mm_i,9); fprintf(fid,'Mean Anom(rad): %s\n', strdop);
%Коэффициенты полинома для учета поправок времени strdop = e_norm(alm_gln(ns).tau_n,9); fprintf(fid,'Af0(s): %s\n', strdop);
Yuma(ns).Af1 = 0.0; fprintf(fid,'Af1(s/s): %11.9e\n', Yuma(ns).Af1);
%Номер недели fprintf(fid,'week: %i \n\n',Yuma(ns).week);
%Номер дня от начала недели, к которому относятся данные альманаха
Yuma(ns).day = weekday;
Функция week_GLONAS_gps
function [wGPS,weekday] = week_GLONAS_gps(Na, timeUTC)
% Имя функции: week_GLONAS_gps
% Функция вычисляет :
% Выходные данные: wGPS, weekday.
%wGPS -номер GPS- недели,
% weekday - день недели (0 - воскресенье, 1 - понедельник, 2 - вторник, 3 - среда,
%4 - четверг, 5 - пятница, 6 - суббота.
% Входные данные:
%Na- номер суток, к которым относятся данные альманаха спутника ГЛОНАСС, отсчитываемые от
%ближайшего високосного года
% timesUTC -время UTC c начала текущей недели (секунда),
%Используемые константы:
% сдвиг времени между UTC и системным временем GPS на начало 2004 года
%diff_UTC_GPS = 13;(до конца 2005 года) diff_UTC_GPS = 14;%(с 2006 года)
161
%количество дней в месяцах
DnMon = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
% юлианский день начала отсчета недели GPS на ноль часов UTC c 5 на 6
% января 1980 года: wGPS_0 = 2444244.5;
% Расчет номеров юлианского дня и дня года, от которого отсчитываются дни ГЛОНАСС
[JD, day_year] = JD_data(timeUTC);
%Сдвиг на ноль часов UTC на гривчском меридиане
JD = JD - 0.5;
%Расчет номера GPS неделе стандартной функцией MatLab "floor" wGPS = floor((JD + Na - wGPS_0) / 7);
%Расчет дня недели стандартной функцией MatLab "mod" weekday = mod((JD+Na - wGPS_0) , 7);
Функция semi_axis
function [a_npp] = semi_axis(A_PZ90_KM, t_dr, i_incl, ecc, omega_n);
%Имя функции:semi_axis
%Функция вычисляет радиус орбиты спутника ГЛОНАСС в соответствии с интерфейсным контроль- ным
%документом ГЛОНАСС
% A_WGS84 = 6378.1370 ; %km
%A_PZ90_KM = A_WGS84;
% A_PZ90_KM = 6378.136; %(Km) Equatorial radius of the Earth
J20 = -1082625.7e-9; % Factor at the second zone harmonicm
% B_PZ90_KM = 6356.75136174571344; % AP_LAND (Km) Polar radius of the Earth */
%a_npp = 1;
MU = 398600.44; % (Km^3/cek^2) constant of a gravitational epsilon = 1.0e-3; sin_i = sin(i_incl); sin_i2 = sin_i * sin_i; v = -omega_n;% omega_n - angle of a perigee ecc2_1 = 1.0 - ecc * ecc; b1 = 2.0 - 2.5 * sin_i2; b2 = sqrt(ecc2_1 * ecc2_1 * ecc2_1 ); b3 = 1.0 + ecc * cos(omega_n); b3 = b3 * b3; b4 = 1.0 + ecc * cos(v); b5 = b4 * b4 * b4 / ecc2_1; b = b1 * b2 / b3 + b5; t_ock = t_dr;
162
tock_2pi = t_ock /(pi * 2); p1_3 = 1.0 / 3.0; a_dop = MU * tock_2pi * tock_2pi; a_n = a_dop^(1/3);
% a_n = pow(a_dop, p1_3); nn = 0; dda = epsilon + 1; while ( (dda > epsilon) & (nn < 50) ) p = a_n * ecc2_1; % Focal parameter ae_p = A_PZ90_KM / p; b0 = 1.0 + (1.5 * J20 * ae_p * ae_p) * b; t_ock = t_dr / b0; tock_2pi = t_ock / (pi * 2); a_dop = MU * tock_2pi * tock_2pi; a_npp = a_dop^(1/3); dda = abs(a_n - a_npp); a_n = a_npp; nn = nn + 1; end;
Функция JD_epohi
function jden = JD_epohi(epoha)
%Имя: JD_epohi
%Фукция JD_epohi(epoha) рассчитывает номер юлианского дня
%до начала года (epoha) на 12h, 0 день, январь.
%Входные данные: epoha, размерность-год
%Выходные данные:
% jden- номер юлианского дня на 12h, 0 день, январь ( размерность -дни) rk = mod(epoha,4); if ( rk == 0 ) rk = 1.0; else rk = 2.0 - rk * 0.25; end; n100 = floor(epoha / 100); n400 = floor(epoha / 400); jden = (4712 + epoha) * 365.25 + n400 - n100 + rk;
% fprintf('epoha=%d rk=%f jden=%6.2f \n', epoha, rk, jden);
163
Функция JD_data
function [JD, day_year] = JD_data(timeUTC)
%Имя:JD_data
% Функция JD_data(timeUTC) вычисляет :
%JD - номер юлианского дня, day_year - номер дня года.
%Входные данные:
%timeUTC.year - год,
% timeUTC.mon - месяц,
% timeUTC.day - день.
%Выходные данные:
%JD - юлианский день, day_year- день от начала года.
%Оригинальные функции: function jd0 = JD_epohi(year),
%(смотри комментарий).
% число дней от начала периода до 12ч. всемирного времени (полдень)
% указанной даты ( по Гринвичу)
%Входные данные для контрольного примера:
%year = 2000; mon = 1; day = 1;
%количество дней в месяцах
DnMon = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
%Вычисление номера юлианского дня опорной эпохи jd0 = JD_epohi(timeUTC.year);
%Учет високосного года nfebr = 0; if mod(timeUTC.year,4) == 0 nfebr = 1; end;
%Расчетномера дня года k = 0; for i = 2 : timeUTC.mon k = k + DnMon(i - 1); if (i == 2) k = k + nfebr; end; end; day_year = k + timeUTC.day;
%Расчет номера юлианского дня
JD = jd0 + day_year;
%2451545 - номер юлианского дня; 1- первый день года.
164
5.2 Декодирование данных альманаха спутников GPS
5.2.1 Краткие сведения из теории
Спутниковые навигационные приемники обрабатывают данные навигационных спутников и вырабатывают информации в определенных форматах, как правило, в форма- тах приемников характерных для фирм-производителей. Рассмотрим программу преобра- зования данных альманаха, получаемых с одного из широко распространенных приемни- ков фирмы Novatel, ProPak-G2, выполненного на базе платы OEM 4 [10, 11]..
Информация, вырабатываемая приемником, записывается в виде файла в бинарной форме и в форме кода ASCII. В программе декодирования альманаха используется бинар- ная форма. Декодируется сырой альманах (RAWALM, Message ID: 74) [10, 11]. Формат данных из оригинала [10, 11].имеет вид
Таблица
Формат данных альманаха RAWALM
Field
Field type
Data Descrsption
Format
Binary
Bytes
Binary
Offset
1 header
Log header
H
0 2 ref week
Almanac reference week number
Ulong
4
H
3 ref secs
Almanac reference time (seconds.)
Ulong
4
H+4 4 subframes
Number of subframes to follow
Ulong
4
H+8 5 svid
SV ID (satellite vehicle ID)a
UShort
2
H+12 6 data
Subframe page data.
Hex
30
H+14 7...
Next subframe offset = H + 12 + (subframe x 32) variable xxxx
32-bit CRC (ASCII and Binary only)
Hex
4
H + 12 +
(32 x subframes) variable [CR][LF]
Sentence terminator (ASCII only)
-
-
-
На CD-диске программа записана в папке ALM_ProPakG2
Цель лабораторной работы изучение процедур и правил декодирования данных навигационного приемника.
5.2.2 Лабораторная работа 5. 2 «Декодирование данных альманаха навигационных приемников на базе плат OEM- 4»
Пакет программ для выполнения лабораторной работы расположен в папке
ALM_ProPakG2 Для выполнения работы в качестве входных данных потребуется альма- нах спутников GPS в виде бинарного файла, который можно получить с навигационных приемников на базе плат OEM 4. Для тестирования программы такой альманах приводит- ся в приложении.
165
Рекомендуется следующий порядок (не обязательный) выполнения лабораторной работы.
1. Создайте папки
ALM_ProPakG2_My, RAW_ALM_BIN_DAT,
RAW_ALM_BIN_ALM.
2. Скопируйте в папку
ALM_ProPakG2_My
из ALM_ProPakG2 папку
RAW_ALM_PRG и RAW_ALM_BIN_DAT, RAW_ALM_BIN_ALM.
3. Запишите в папку RAW_ALM_BIN_DAT бинарный файл альманаха, который требу- ется декодировать.
4. Откроете папку RAW_ALM_PRG и изучите файлы и функции по комментария и про- граммным процедурам.
5. В папке RAW_ALM_PRG откройтефайл Raw_Alm.m и в соответствующие строки запишите имя файла альманаха, который требуется декодировать и имя файла альма- наха, под которым запишется декодированный альманах в формате YUMA.
6. Задание. Сравните данные декодированного альманаха с содержанием альманаха
YUMA на сайте
7. Выполните файл Raw_Alm.m и в папке RAW_ALM_BIN_ALM прочитайте декодиро- ванный альманах.
5.2.3 Задание для самоподготовки
В руководстве пользователя на приемники, выполненные на базе плат OEM 4 [10,
11]. Изучите соответствующие разделы.
5.2.4 Функции и файлы из папки RAW_ALM_PRG
Файл Raw_Alm.m
%Имя файла:Raw_Alm.m
%Программа преобразует данные бинарного файла альманаха спутников GPS в
%формат YUMA
%Входные данные:
%Бинарный файл "сырого" альманаха навигационных приемников ProPak G2, ProPak G2 plus, DL 4
%для приведенного примера имя файла (2006_05_09_16.31_G2_RAW).bin
%Выходные данные:
%Файл альманаха GPS в формате YUMA,
%для приведенного примера имя файла (2006_05_09_16.31_G2_RAW_OK).alm
%начало цикла для 63 спутников (в том числе и спутников-нагрузок) for i=1:63
%структура согласно формата альманаха YUMA и ICD GPS-200 C alm(i)= struct('SV_ID',0, 'week',0, 'health_0',0, 'health',0,'ecc',0,...
't0a',0, 'incl_angle',0, 'omega_dot',0,...
ЧЧЧЧЧЧЧЧЧЧЧЧЧЧЧЧЧЧЧЧЧ.
166
'sqrtA',0, 'omega0',0, 'omega', 0, 'M0',0,'af0',0, 'af1',0); end;
%диск и последовательность папок, из которых считывается файл с записью альманаха в бинарном
%формате. При необходимости можно указать любой другой путь к файлу name_in = 'E:\MATLAB7\work\RAW_ALM_BIN_DAT\(2006_05_09_16.31_G2_RAW).bin';
%открытие данных для чтения и записи fw = fopen('E:\MATLAB7\work\RAW_ALM_BIN_ALM\(2006_05_09_16.31_G2_RAW_OK).alm','Wt');
%формирование CRC-кода для проверки считываемых данных fid_crc = fopen(name_in,'rb'); fseek(fid_crc, 0, 'eof'); size1 = ftell(fid_crc)- 4; fseek(fid_crc, size1, 'bof');
CRC = fread(fid_crc,1,'uint32'); fprintf('CRC_hex = %X \n',CRC); fseek(fid_crc, 0, 'bof');
CRC_own = CalculateBlockCRC32(fid_crc, size1); fprintf('CRC_own_hex= %X \n',CRC_own); fclose(fid_crc);
%Чтение "сырых" данных альманаха fid = fopen(name_in,'rb');
% Заголовок данных- Header Structure (Binary Message)
Syn_c = fread(fid,3,'char'); %В OEM-4 'Char'
Hear_d= fread(fid,1,'uchar');%В OEM-4 'Uchar'
Message_ID = fread(fid,1,'uint16');%В OEM-4 'Long'
Message_Type=fread(fid,1,'char'); %В OEM-4 'Char'
Port_Address = fread(fid,1,'char'); %В OEM-4 'Char'
Message_Length = fread(fid,1,'uint16'); %В OEM-4 'Ushort'
Sequenc_e = fread(fid,1,'uint16'); %В OEM-4 'Ushort'
Idle_Time = fread(fid,1,'char'); %В OEM-4 'Char'
Time_Status = fread(fid,1,'char'); %В OEM-4 'Enum' week_1 = fread(fid,1,'uint16'); %В OEM-4 'Ushort'
Time = fread(fid,1,'ulong'); %В OEM-4 'Double'
N_ul = fread(fid,1,'uint'); d_fcc=fread(fid,1,'uint16'); d_fcc = dec2hex(d_fcc); e1 = fread(fid,1,'uint16'); week =fread(fid,1,'uint');% неделя GPS
T0a= fread(fid,1,'ulong');% время привязки альманаха kol_str= fread(fid,1,'ulong'); %количество строк
%Чтение данных по каждому спутнику for i=1:kol_str
167
Nsv= fread(fid,1,'ushort'); % номер спутника согласно табл.20-V ICD GPS-200 C
TLM= fread(fid,3,'char');
TL_M = dec2hex(TLM,2);
HOW= fread(fid,3,'char');
HO_W = dec2hex(HOW,2); e6= fread(fid,1,'char'); maska = bin2dec('11000000');
ID_dop = bitand(e6, maska); data_ID = bitshift(ID_dop, -6); maska = bin2dec('00111111');
SV_ID = bitand(e6, maska); ecc1 = fread(fid,1,'char'); ecc2 = fread(fid,1,'char'); ecc = (bitshift(ecc1, 8) + ecc2)*2^(-21);%ексцентриситет t0a= fread(fid,1,'char')*2^12 ;% время привязки альманаха d1= fread(fid,1,'bit8'); d2=fread(fid,1,'char'); incl_angle= ((bitshift(d1, 8) + d2)*2^(-19) + 0.3) * pi;%наклонение орбиты спутника d1= fread(fid,1,'bit8'); d2=fread(fid,1,'char'); omega_dot= (bitshift(d1, 8) + d2)*2^(-38) * pi;%скорость изменения восходящего узла maska = bin2dec('00111111');
%здоровье спутника GPS health_1 = fread(fid,1,'char'); if health_1 > 255 health_1 = health_1; end; health = bitand(health_1,maska); if health > 63 health = health; end; maska = bin2dec('11000000'); health_0 = bitand(health_1,maska); health_0 = bitshift(health_0, -6); if health_0 > 3 health_0 = health_0; end; d1 = fread(fid,1,'char'); d2 = fread(fid,1,'char');
168
d3 = fread(fid,1,'char'); sqrtA=(bitshift(d1, 16) + bitshift(d2, 8) + d3)*2^(-11);%корень квадратный из большей полуоси орбиты d1 = fread(fid,1,'bit8'); sign2 = sign(d1); d2 = fread(fid,1,'char'); d3 = fread(fid,1,'char'); omega0 =((bitshift(d1, 16) + bitshift(d2, 8) + d3)*2^(-23)) * pi; %долгота восходящего узла орбиты спут- ника d1 = fread(fid,1,'bit8'); sign1 = sign(d1); d2 = fread(fid,1,'char'); d3 = fread(fid,1,'char'); omega =((bitshift(d1, 16) + bitshift(d2, 8) + d3)*2^(-23)) * pi; %аргумент перигея d1 = fread(fid,1,'bit8'); d2 = fread(fid,1,'char'); d3 = fread(fid,1,'char');
M0 =((bitshift(d1, 16) + bitshift(d2, 8) + d3)*2^(-23)) * pi;%средняя аномалия на время привязки t0a af0_1 = fread(fid,1,'char'); af0_1_bin = dec2bin(af0_1); d2 = fread(fid,1,'char'); d2bin = dec2bin(d2); d3 = fread(fid,1,'char'); d3bin = dec2bin(d3); maska = bin2dec('00011100'); d4 = bitand(d3,maska); af0_2 = bitshift(d4, -2); af0_dop = bitshift(af0_1, 3) + af0_2; koef = 2^(-20); af0 = var_sign(af0_dop,11,koef);%коэффициент коррекции времени af0
% af1: maska = bin2dec('11100000'); d4 = bitand(d3,maska); af1_2 = bitshift(d2, 3); af1_1 = bitshift(d4, -5);
169
af1_dop = bitor(af1_2, af1_1); koef = 2^(-38); af1 = var_sign(af1_dop,11,koef);%коэффициент коррекции времени af1
%-------------------------------------------------- if SV_ID > 0 alm(SV_ID)= struct('SV_ID',SV_ID, 'week',week, 'health_0', health_0, 'health',health, 'ecc',ecc,...
't0a',t0a, 'incl_angle',incl_angle, 'omega_dot',omega_dot,...
'sqrtA',sqrtA, 'omega0',omega0, 'omega', omega, 'M0',M0,...
'af0',af0, 'af1',af1); end; end;
%Формирование вывода данных альманаха в формате YUMA for i=1:30 %63 if alm(i).SV_ID > 0 if alm(i).SV_ID < 10 fprintf(fw,'**** Week %i almaNAU for PRN-0%i **********\n',alm(i).week, alm(i).SV_ID); fprintf(fw,'ID: 0%i\n',alm(i).SV_ID); else fprintf(fw,'**** Week %i almaNAU for PRN-%i **********\n',alm(i).week, alm(i).SV_ID); fprintf(fw,'ID: %i\n',alm(i).SV_ID); end; if alm(i).health < 10 fprintf(fw,'Health: %i0%i\n', alm(i).health_0, alm(i).health); else fprintf(fw,'Health: %i%i\n', alm(i).health_0, alm(i).health); end; strdop = e_norm(alm(i).ecc, 10); fprintf(fw,'Eccentricity: %s\n', strdop); fprintf(fw,'Time of Applicability(s): %6.4f\n',alm(i).t0a); fprintf(fw,'Orbital Incluation(rad): %0.10f \n',alm(i).incl_angle); strdop = e_norm(alm(i).omega_dot, 10); fprintf(fw,'Rate of Right Ascen(r/s): %s\n', strdop);
170
fprintf(fw,'SQRT(A) (m^1/2): %4.7f \n',alm(i).sqrtA); strdop = e_norm(alm(i).omega0, 10); fprintf(fw,'Right Ascen at Week(rad): %s\n', strdop); if alm(i).omega < 0 fprintf(fw,'Argument of Perigee(rad): %1.10f \n',alm(i).omega); else fprintf(fw,'Argument of Perigee(rad): %1.10f \n',alm(i).omega); end; strdop = e_norm(alm(i).M0, 10); fprintf(fw,'Mean Anom(rad): %s\n', strdop); strdop = e_norm(alm(i).af0, 10); fprintf(fw,'Af0(s): %s\n', strdop); strdop = e_norm(alm(i).af1, 10); fprintf(fw,'Af1(s/s): %s\n', strdop); fprintf(fw,'week: %i \n',alm(i).week); fprintf(fw,' \n'); end; %if alm(i).SV_ID > 0 end; %i
%CR_C = fread(fid,1,'uint32');
%CR_C = dec2hex(CR_C); fclose(fid);%закрытия файла для чтения fclose(fw);%закрытия файла для записи
Скачать 4.49 Mb.