Файл: Учебное пособие ( Лабораторный практикум на компьютере ) Київ 2008 1.pdf

Скачать файл
Заказать решение

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 11.01.2024

Просмотров: 295

Скачиваний: 5

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.

СОДЕРЖАНИЕ

Функция Read_GL_Alm function [alm,max_kol] = Read_GL_Alm(Dat); %Имя функции: Read_GL_Alm %Функция читает данные альманаха спутников ГЛОНАСС. %Входные данные: %Dat - имя файла, содержащего альманах спутников ГЛОНАСС , например, Dat = 'GLN_all.alm'; %Выходные данные: %alm - данные альманаха спутников ГЛОНАСС (структура) , %max_kol- максимальное количество спутников в альманахе fid =fopen(Dat,'rt');%открытие файла max_kol = 0; %чтение данных из файла while not(feof(fid)) s1=fscanf(fid,'%s',6); if not(feof(fid)) lenstr = length(s1); max_kol = max_kol + 1; % while (fscanf(fid,'%s',1) == '-') end str1 = fscanf(fid,'%s',1); % lenstr = length(str1); str2 = fscanf(fid,'%s',1); lenstr = length(str2); n_sv = sscanf(str2,'%d'); strID = str2(1:lenstr); ID = sscanf(strID,'%d'); alm(ID).ID = ID; %номер спутника (1...24) %%%%%%%%%%% t_1=fscanf(fid,'%s',1); alm(ID).Hn=fscanf(fid,'%d'); % номер частоты спутника t_2=fscanf(fid,'%s',1); alm(ID).Health=fscanf(fid,'%d',1);%здоровье (0- спутник здоров) t_3=fscanf(fid,'%s',1); alm(ID).ecc = fscanf(fid,'%g',1);% эксцентриситет орбиты спутника t_4=fscanf(fid,'%s',1); 108 alm(ID).Na =fscanf(fid,'%g',1 );% номер суток к которым относятся данные альманаха, отсчиты- ваемых %от ближайшего високосного года t_5=fscanf(fid,'%s',3); alm(ID).deltai = fscanf(fid,'%g',1); % наклонение орбиты спутника (радианы) t_6=fscanf(fid,'%s',2); alm(ID).LambdaN = fscanf(fid,'%g',1); % долгота восходящего узла орбиты спутника (радианы) while not(fscanf(fid,'%c',1) == ':') end % t_7=fscanf(fid,'%s',1) alm(ID).TLambdaN = fscanf(fid,'%g',1) - 10800;% время прохождения восходящего узла орбиты спутника, %к которому относятся данные альманаха, приведенное к времени UTC (секунда) t_8=fscanf(fid,'%s',2); alm(ID).omegan = fscanf(fid,'%g',1); %аргумент перигея (радианы) t_11=fscanf(fid,'%s',1); alm(ID).Tdr=fscanf(fid,'%g',1); % драконический период (секунда/виток) t_12=fscanf(fid,'%s',1); alm(ID).dTdr=fscanf(fid,'%g',1); % скорость изменения драконического периода (секунда/виток в квадрате) t_13=fscanf(fid,'%s',2); alm(ID).tau_n=fscanf(fid,'%g',1) / 1000.0;%коэффициет коррекции шкалы времени (сдвиг времени спутника %относительно системного времени ГЛОНАСС в секундах) end; end fclose(fid); Функция Gln_data_from_NA function [time_UTC] = Gln_data_from_NA(leap_year, day_from_leap); %Имя функции:Gln_data_from_NA %Функция предназначена для преобразования номера дня NA (день привязки альманаха % от ближайшего високосного года) в текущую дату %function [time_UTC] = Gln_data_from_NA(leap_year, day_from_leap); %Входные данные: % leap_year - ближайший високосный год %day_from_leap (NA) - день привязки альманаха %Выходные данные: %структура timeUTC (год, месяц, день) - текущая дата DnMon = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]; %количество дней в месяцах 109 n4 = mod(leap_year, 4); n100 = mod(leap_year, 100); n400 = mod(leap_year, 400); if (n4 == 0) n_leap = 1; else n_leap = 0; end if ((n100 == 0) & (n400 > 0)) n_leap = 0; end if (day_from_leap > (365 + n_leap)) day = day_from_leap - (365 + n_leap); k =fix (day / 365); day = mod(day, 365); else day = day_from_leap; k = -1; end; god = leap_year + k + 1; if (god > leap_year) n_leap = 0; end; mon = 1; mon_day = 31; while (day > mon_day) day = day - mon_day; mon = mon + 1; % LK mon_day = DnMon(mon); % mon = mon+1; if (mon == 2) mon_day = mon_day + n_leap; end end time_UTC.year = god; time_UTC.mon = mon; time_UTC.day = day; Функция JD_data function [JD, day_year] = JD_data(timeUTC) %Имя:JD_data % Функция JD_data(timeUTC) вычисляет : 110 %JD - номер юлианского дня, day_year - номер дня года. %Входные данные: %Структура timeUTC %timeUTC.year - год, % timeUTC.mon - месяц, % timeUTC.day - день. %Выходные данные: %JD - юлианский день; %day_year- день от начала года. %количество дней в месяцах 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; Функция LLH_to_ECEF function [XYZ] = LLH_to_ECEF( a, b, llh) %Имя функции:LLH_to_ECEF.m %Функция преобразования географических координат в прямоугольную геоцентрическую систему координат (ECEF) %Входные данные: %Структура llh %llh.lon-долгота (радиан), %llh.lat-широта (радиан), %llh.h-высота (метр); %a, b-большая и малая полуоси земного эллипсоида в WGS-84 (метр); %Выходные данные: 111 %Структура XYZ %XYZ.x - координата x в ECEF; %XYZ.y - координата y в ECEF; %XYZ.z - координата z в ECEF; a2=a * a; b2=b * b; r = a2 / sqrt(a2 * cos(llh.lat) * cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat)); XYZ.x = (r + llh.h) * cos(llh.lat) * cos(llh.lon); XYZ.y = (r + llh.h) * cos(llh.lat) * sin(llh.lon); XYZ.z = (b2 / a2 * r + llh.h) * sin(llh.lat); Функция WGS84_to_PZ90 function [pos_pz90] = WGS84_to_PZ90( pos_wgs84) %Имя функции: WGS84_to_PZ90.m %Функция преобразует координаты из системы WGS 84 в систему ПЗ 90 p = (1.0 - 0.12e-6); x = pos_wgs84.x; y = pos_wgs84.y; z = pos_wgs84.z; pos_pz90.x = (p * x + p * 0.82e-6 * y + 1.4 ); pos_pz90.y = (-p * 0.82e-6 * x + p * y + 1.4); pos_pz90.z = (p * z + 0.9); Функция ECEF_to_LLH function [llh] = ECEF_to_LLH(a, b, XYZ) %Имя функции:ECEF_to_LLH %Функция преобразует координаты из системы ECEF в географическую систему %Входные данные: %a, b-большая и малая полуоси земного эллипсоида (метр); ( a=6378137.0; %b=6356752.314; метры) %Структура XYZ %XYZ.x - координата x в ECEF; %XYZ.y - координата y в ECEF; %XYZ.z- координата z в ECEF; %Выходные данные: %Структура llh %llh.lon - долгота (радиан); %llh.lat - широта (радиан); %llh.h-высота (метр); a2 = a * a; b2 = b * b; xy = sqrt(XYZ.x * XYZ.x + XYZ.y * XYZ.y); thet = atan2(XYZ.z * a , (xy * b)); 112 %thet = atan(XYZ.z * a / (xy * b)); esq = 1.0 - b2 / a2; epsq = a2 / b2 - 1.0; llh.lat = atan((XYZ.z + epsq * b * (sin(thet) ^3)) / (xy - esq * a * (cos(thet) ^3))); llh.lon = atan2(XYZ.y, XYZ.x);%! while (llh.lon < 0) llh.lon = 2*pi + llh.lon; end; %if llh.lon < 0 % llh.lon = 2*pi + llh.lon; %end ; r = a2 / sqrt(a2 * cos(llh.lat) * cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat)); llh.h = xy / cos(llh.lat) - r; end Функция GLN_satfind function [kol_gln_a, satvis_gln, satpos_gln_ns, satpos_gln_ecef] = GLN_satfind(a, b, timeUTC, cur- rent_loc_pz90, alm_gln); %Имя функции:GLN_satfind %Функция вычисляет координаты спутников ГЛОНАСС, углы видимости и количество видимых спутников %Входные данные: %a, b- большая и малая полуоси земного эллипсоида; %timeUTC- параметры времени; %current_loc_pz90- широта, долгота и высота точки, из которой наблюдаются спутники (структура); % alm_gln-альманах спутников ГЛОНАСС (структура) %Выходные данные: %kol_gln_a- количество видимых спутников ГЛОНАСС, %satvis_gln- углы видимости спутников (структура), %satpos_gln_ns-координаты спутников (структура) в абсолютной системе координат (ECI) %satpos_gln_ecef- координаты спутников (структура) в подвижной системе координат (ECEF) RAD_TO_DEG = 180.00 / pi ; KOL_GLN = 24; nom_ns = 1; Min_Elev = 0; % минимальный угол видимости kol_gln_a = 0; [nom_sat_gln_a, trac_sat_gln_a, sat_pos, satvis_gln] = init_data(KOL_GLN); % число дней от фундаментальной эпохи 2000г от 12h, 0, January: % jd_up_epoh = JD_from_epohi( 2000, timeUTC) - 0.5; % from 12h UTC 2000 nut = 0; 113 S0 = s0_Nut( timeUTC, nut); time_s0 = S0.s0_m_mod ; %time_s0 - истинное звездное время в день обсервации, 0h UTC year = timeUTC.year; leap_year = fix( year / 4) * 4; % ближайший к текущему (предыдущий) високосный год ti = timeUTC.ti; % текущее время обсервации от начала дня n00 = fix(ti / 86400); % n0 - номер текущих суток внутри 4-х летнего периода (от ближайшего високосного года) n0 = JD_from_epohi(leap_year, timeUTC) + n00 + 1; [eci_current_loc, eci_rec_pos_xyz] = llh_to_eci(a, b, ti, time_s0, current_loc_pz90); satpos_eci = init_satpos_gln(); for ( i = 1 : KOL_GLN) prn = alm_gln (i).ID; health = alm_gln (i).Health; if ( (prn > 0) & (health == 0)) satpos_eci = gln_a_1(i, n0, ti, time_s0, alm_gln); [satpos_eci, satpos_gln_a] = satpos_eci_in_metr(satpos_eci); %Transformation coordinates from "Km" in "m" satpos_gln_ns(prn) = satpos_gln_a; %----- [satpos_ecef] =eci_to_ecef(time_s0, ti, satpos_eci); satpos_gln_ecef(prn) = satpos_ecef; % calculates angles of visibility by almanac: [top] = top_coord(prn, eci_current_loc, eci_rec_pos_xyz,satpos_gln_a); satvis_gln(prn) = top; if (satvis_gln(prn).r > 0.0) el = RAD_TO_DEG * satvis_gln(prn).el; if ( el >= Min_Elev) % fprintf('prn=%2i el=%f az=%f \n ', prn, el, RAD_TO_DEG * satvis_gln(prn).az); satpos_gln( prn) = satpos_eci; [satpos_ecef] =eci_to_ecef(time_s0, ti, satpos_eci); [trac_sat_gln_a] = rewrite_satpos(nom_ns, satpos_ecef); trac_sat_gln_a( nom_ns).r = satvis_gln( prn).r; trac_sat_gln_a (nom_ns).prn = prn; nom_sat_gln_a(nom_ns) = prn; nom_ns = nom_ns + 1; % number of visible satellites end; % if ( el >= Min_Elev) end; % if ( satvis_gln[prn-1].r > 0.0) 114 end; % f ( (prn > 0) & (health == 1)) end; % for ( i = 1 : KOL_GLN) kol_gln_a = nom_ns - 1; Функция ris_vis_sat function ris_vis_sat(max_n, nj, ti_start, step_t, ris_vid) %Имя функции:ris_vis_sat %Функция строит графики спутников, видимых из точки наблюдения mm = step_t/3600.0; j_color = 0; color7(1:7) = ['b' 'g' 'r' 'c' 'm' 'y' 'k']; for i = 1 : max_n str1 = ris_vid(i, 1: nj ); k = find(str1); if k > 0 len = length(k); if len == 1 hPlot = plot(k(1):k(1),ris_vid(i,k(1):k(1)), s); set(hPlot, 'LineWidth', 3); hold on; else n0 = k(1); j_color = j_color + 1; if (j_color > 7 ) j_color = 1; end s = color7(j_color); for j = 2 : len if (k(j) > (k(j-1) + 1)) | (j == len) k0 = k(j - 1); hPlot = plot(n0*mm : mm : k0*mm,ris_vid(i,n0:k0), s); set(hPlot, 'LineWidth', 3); hold on; n0 = k(j); end; % if end; % for j = 2 : len end; % else end; % if k > 0 end; % for i = 1 : max_n 115 grid on hAxes=gca; My = [0:2:25]; kol_x = ti_start + nj * step_t; Mx = [0 : 2 : (kol_x - ti_start) / 3600]; set(hAxes, 'ytick', My,'xtick',Mx) ; hText = gca; set(hText,'FontSize',12,'FontName','TimesNewRoman') xlabel('Время наблюдения спутников ГЛОНАСС, час') ylabel('Номер навигационного спутника') Функция 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); Функция init_data function[nom_sat_gln_a, trac_sat, sat_pos, satvis_gln] = init_data(kol) %Имя функции: init_data %Функция предназначена для инициализации массивов sat_pos.prn = 0; sat_pos.x = -1.0; sat_pos.y = -1.0; sat_pos.z = -1.0; sat_pos.vx = -1.0; sat_pos.vy = -1.0; sat_pos.vz = -1.0; 116 sat_pos.dvx = -1.0; sat_pos.dvy = -1.0; sat_pos.dvz = -1.0; sat_pos.range_m = -1.0; for( i = 1 : kol) nom_sat_gln_a(i) = 0; satvis_gln(i).el = -1; satvis_gln(i).az = -1; satvis_gln(i).r = 0; trac_sat(i).prn = 0; trac_sat(i).x = -1.0; trac_sat(i).y = -1.0; trac_sat(i).z = -1.0; trac_sat(i).vx = -1.0; trac_sat(i).vy = -1.0; trac_sat(i).vz = -1.0; trac_sat(i).range_m = -1.0; end; Функция JD_from_epohi function [jd] =JD_from_epohi( epoha, timeUTC); %Имя функции:JD_from_epohi.m %Функция вычисляет jd - количество дней от указанного года (epoha) % до текущей даты, указанной в структуре timeUTC, представленной в виде % (timeUTC.year, timeUTC.mon, timeUTC.day) jd0 = JD_epohi(epoha) + 1; % 12h, 1 den January [day, day_year] = JD_data(timeUTC); jd = day - jd0; Функция s0_Nut function [S0] = s0_Nut( timeUTC, nut) %Имя функции: s0_Nut %Функция рассчитывает истинное или среднее звездное время на 0ч UTC %Входные данные: %timeUTC - дата, на которую требуется рассчитать истинное или среднее звездное время %nut- признак (если nut= 0, то вычисляется среднее звездное время без учета нутации, иначе вычис- ляется %истинное звездное время) %Выходные данные: %S0 - истинное или среднее звездное время на 0ч UTC jd2000 = 2451545; % 12h UTC 1 января % Применить функцию JD_data [jd, day_year] = JD_data( timeUTC); if (jd == NaN) 117 s0_mod = NaN; h = NaN; min = NaN; sec = NaN; fprintf('function s0_m - end0 \n'); return; end; jd = jd - 0.5; d = jd - jd2000; t = d / 36525.0; % 36525 - юлианский период 100 лет t2 = t * t; h1 = 24110.54841; %h1=6.0*3600.0+41.0*60.0+50.54841; % h2 = 236.555367908 * d; h2 = 8640184.812866 * t ; h3 = 0.093104 * t2; h4 = t2 * t * 6.2E-6; if ( nut == 0) na = 0; else na = utc_nut(t); end; s0_m = h1 + h2 + h3 - h4; S0.s0_nut = s0_m + na; S0.s0_m_mod = mod(s0_m, 86400); s0_day = floor(s0_m / 86400); S0.s0_m_hour = S0.s0_m_mod / 3600.0; S0.s0_m_hour = floor(S0.s0_m_mod / 3600); sec_min = S0.s0_m_mod - S0.s0_m_hour * 3600; S0.s0_m_min = floor(sec_min / 60); S0.s0_m_sec = sec_min - S0.s0_m_min * 60; S0.s0_nut_mod = mod(S0.s0_nut, 86400); s0_day = floor(S0.s0_nut / 86400); S0.s0_nut_hour = S0.s0_nut_mod / 3600.0; S0.s0_nut_hour = floor(S0.s0_nut_mod / 3600); sec_min = S0.s0_nut_mod - S0.s0_nut_hour * 3600; S0.s0_nut_min = floor(sec_min / 60); S0.s0_nut_sec = sec_min - S0.s0_nut_min * 60; Функция llh_to_eci function [eci_llh, eci_xyz] = llh_to_eci(a, b, ti, time_s0, llh_loc) ; %Имя функции:llh_to_eci %Функция вычисляет позицию приемника в абсолютной геоцентрической системе координат (ECI) %Входные данные: %a, b-большая и малая полуоси земного эллипсоида (метр); 118 %ti- текущее время (секунды) , %time_s0- истиное звездное время, %Структура llh_loc - координаты приемника; { %llh_loc.lon-долгота (радиан); %llh_loc.lat-широта (радиан); %llh_loc.h-высота (метр); %Выходные данные: %Структура eci_llh - географические координаты приемника в абсолютной геоцентрической системе координат (ECI) %eci_llh.lon - долгота (радиан); %eci_llh.lat - широта (радиан); %eci_llh.h - высота (метр); %Структура eci_xyz- координаты приемника в абсолютной прямоугольной геоцентрической системе координат (ECI) % eci_xyz.x - координата x; % eci_xyz.y - координата y; % eci_xyz.z - координата z; OMEGA_Z = 0.7292115e-4; % угловая скорость вращения Земли, (рад/cek) SEC_IN_RAD = 7.2722052166430e-005; % (PI / 43200.0 ) // Number radian eci_llh.lon = llh_loc.lon + ti * OMEGA_Z + SEC_IN_RAD * time_s0; eci_llh.lat = llh_loc.lat; eci_llh.h = llh_loc.h; eci_xyz = llh_to_ecef( a, b, eci_llh); Функция init_satpos_gln function [sat_pos] = init_satpos_gln() ; %Имя функции:init_satpos_gln %Функция предназначена для инициализации структуры sat_pos sat_pos.prn = 0; sat_pos.x = -1.0; sat_pos.y = -1.0; sat_pos.z = -1.0; sat_pos.vx = -1.0; sat_pos.vy = -1.0; sat_pos.vz = -1.0; sat_pos.dvx = -1.0; sat_pos.dvy = -1.0; sat_pos.dvz = -1.0; sat_pos.range_m = -1.0; Функция gln_a_1 function [satpos_xyz_gln] = gln_a_1(ns, n0, ti_current, time_s0, alm_gln); 119 %Имя функции:gln_a_1 %Функция рассчитывает координаты и скорости спутников ГЛОНАСС в соответствии с интерфейс- ным %контрольным документом ГЛОНАСС %Входные данные: %ns- номер спутника , %n0 - номер текущих суток внутри 4-х летнего периода (от ближайшего високосного года), %ti_current - текущее время обсервации от начала дня, %time_s0 - истинное звездное время в текущий момент обсервации, %alm_gln - альманах спутников ГЛОНАСС %Выходные данные: %Структура satpos_xyz_gln %satpos_xyz_gln.x - координата по оси x; %satpos_xyz_gln.y- координата по оси y; %satpos_xyz_gln.z- координата по оси z; %satpos_xyz_gln.vx- скорость по оси x; %satpos_xyz_gln.vy - скорость по оси y; %satpos_xyz_gln.vz- скорость по оси z ; % I_MID - Mean value of an inclination of a plane of orbit of a satellite I_MID = 1.0995574287564; %(double)(PI * 63.0 / 180.0) T_DR_MID = 43200.0; %Mean value a dragon of cycle time of a satellite A_PZ90_KM = 6378.136; %(Km) Equatorial radius of the Earth MU = 398600.44; % (Km^3/cek^2) constant of a gravitational OMEGA_Z = 0.7292115e-4; %скорость вращения Земли (Angular speed of rotation of the Earth, рад/cek) 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; i_incl = alm_gln(ns).deltai; ecc2_1 = 1.0 - alm_gln(ns).ecc * alm_gln(ns).ecc; % t_dr = T_DR_MID + alm_gln(ns).Tdr; t_dr = alm_gln(ns).Tdr; 120 n_dr = pi * 2 / t_dr; %1. Вычисление полуоси a_n: a_n = semi_axis_1( 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_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; 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); 121 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 -... 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 *... 122 (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; omega_i = omega_0 + dd_omega; ii_incl = i_incl + ddi; ll_z = ma + alm_gln(ns).omegan + n_dr * tau + dd_ll; 123 % ll_z = ma + alm_gln(ns).omegan + n_dr * (ti - t_lambda_k) + dd_ll; mm_i = ll_z - w_i; % 6. Coordinates and components of satellite velocity vector % in absolute coordinate system OXaYaZa at the moment of ti time : % Ei(n) = Mi + ei * sin Ei(n-1), Ei(0) = Mi, | Ei(n) - Ei(n-1) | <= 10-8, eps_kepler = 1E-8; ee_i = kepler( mm_i, ecc_i, eps_kepler);% ee_i - эксцентрическая аномалия % alm_mca(ns).ek = ee_i; v_i = 2.0 * atan2(sqrt(1.0 + ecc_i) * tan(ee_i / 2.0) , sqrt(1.0 - ecc_i) ); u_i = v_i + w_i; ri = ai * (1.0 - ecc_i * cos(ee_i)) ; ecc_i2 = 1.0 - ecc_i * ecc_i; sqrt_mu_ai = sqrt( MU / ai); vr_i = sqrt_mu_ai * (ecc_i - sin(v_i) ) * ecc_i2; vu_i = sqrt_mu_ai * (1.0 + ecc_i * cos(v_i) ) * ecc_i2; cos_ui = cos(u_i); sin_ui = sin(u_i); cos_ii = cos(ii_incl); sin_ii = sin(ii_incl); cos_omega_i = cos(omega_i); sin_omega_i = sin(omega_i); x_dop = (cos_ui * cos_omega_i - sin_ui * sin_omega_i * cos_ii); y_dop = (cos_ui * sin_omega_i + sin_ui * cos_omega_i * cos_ii); satpos_xyz_gln.x = ri * x_dop; satpos_xyz_gln.y = ri * y_dop; satpos_xyz_gln.z = ri * sin_ui * sin_ii; satpos_xyz_gln.vx = vr_i * x_dop - vu_i * (sin_ui * cos_omega_i +... cos_ui * sin_omega_i * cos_ii); satpos_xyz_gln.vy = vr_i * y_dop - vu_i * (sin_ui * sin_omega_i -... cos_ui * cos_omega_i * cos_ii); satpos_xyz_gln.vz = vr_i * sin_ui * sin_ii + vu_i * cos_ui * sin_ii; Функция satpos_eci_in_metr function [satpos_eci, satpos_gln] = satpos_eci_in_metr(satpos_eci); %Имя функции: satpos_eci_in_metr %Функция преобразует координаты satpos_eci (структура), заданные в километрах в координаты %satpos_eci, satpos_gln в метрах satpos_eci.x = satpos_eci.x * 1000.0; satpos_eci.y = satpos_eci.y * 1000.0; 124 satpos_eci.z = satpos_eci.z * 1000.0; satpos_eci.vx = satpos_eci.vx * 1000.0; satpos_eci.vy = satpos_eci.vy * 1000.0; satpos_eci.vz = satpos_eci.vz * 1000.0; satpos_gln.x = satpos_eci.x; satpos_gln.y = satpos_eci.y; satpos_gln.z = satpos_eci.z; Функция init_satpos_gln function [sat_pos] = init_satpos_gln() ; %Имя функции:init_satpos_gln %Функция предназначена для инициализации структуры sat_pos sat_pos.prn = 0; sat_pos.x = -1.0; sat_pos.y = -1.0; sat_pos.z = -1.0; sat_pos.vx = -1.0; sat_pos.vy = -1.0; sat_pos.vz = -1.0; sat_pos.dvx = -1.0; sat_pos.dvy = -1.0; sat_pos.dvz = -1.0; sat_pos.range_m = -1.0; Функция eci_to_ecef function [satpos_ecef] =eci_to_ecef(s0, ti, satpos_eci) %Имя функции:eci_to_ecef %Функция преобразования координат %Входные данные: s0 - истинное звездное время в текущий момент обсервации , %ti - текущее время; satpos_eci %Структура satpos_eci %satpos_eci.x - координата x в абсолютной неподвижной системе координат (ECI); %satpos_eci.y - координата y в абсолютной неподвижной системе координат (ECI); %satpos_eci.z - координата z в абсолютной неподвижной системе координат (ECI); %satpos_eci.vx - скорость vx в абсолютной неподвижной системе координат (ECI); %satpos_eci.vy - скорость vy в абсолютной неподвижной системе координат (ECI); % satpos_eci.vz - скорость vz в абсолютной неподвижной системе координат (ECI); %Выходные данные: % Структура satpos_ecef %satpos_ecef.x - координата x в подвижной системе координат (ECEF); %satpos_ecef.y - координата y в подвижной системе координат (ECEF); %satpos_ecef.z - координата z в подвижной системе координат (ECEF); 125 %satpos_ecef.vx - скорость по оси x в подвижной системе координат (ECEF); %satpos_ecef.vy - скорость по оси z в подвижной системе координат (ECEF); %satpos_ecef.vz- скорость по оси z в подвижной системе координат (ECEF); %Коэффициенты % SEC_IN_RAD - коэффициент преобразования секунд в радианы % s0(radian) = s0 (sek) * SEC_IN_RAD, where % SEC_IN_RAD = 2 * pi / (24 * 3600) = pi / 43200 SEC_IN_RAD = pi / 43200; OMEGA_Z = 0.7292115e-4; %( скорость вращения Земли (angular speed of rotation of the Earth, рад/cek) s_zv = s0 * SEC_IN_RAD + OMEGA_Z * ti; cos_s = cos(s_zv); sin_s = sin(s_zv); satpos_ecef.x = satpos_eci.x * cos_s + satpos_eci.y * sin_s; satpos_ecef.y = -satpos_eci.x * sin_s + satpos_eci.y * cos_s; satpos_ecef.z = satpos_eci.z; satpos_ecef.vx = satpos_eci.vx * cos_s + satpos_eci.vy * sin_s + OMEGA_Z * satpos_ecef.y; satpos_ecef.vy = -satpos_eci.vx * sin_s + satpos_eci.vy * cos_s - OMEGA_Z * satpos_ecef.x; satpos_ecef.vz = satpos_eci.vz; Функция top_coord function [satvis] = top_coord(prn, rec_llh, rec_xyz, nlo_xyz) % Имя функции: top_coord % Функция выполняет расчет топоцентрических координат объекта по заданным %географическим (долгота, широта, высота) и геоцентрическим (x, y, z) %координатам приемника, а также геоцентрическим координатам объекта(x, y, z) %Входные данные: %prn - номер спутника %структура % rec_llh.lat - широта (рад) приемника %rec_llh.lon - долгота (рад) приемника %rec_llh.h - высота (м) приемника %структура %прямоугольные геоцентрические координаты приемника (м) % rec_xyz.x %rec_xyz.y %rec_xyz.z %структура 126 % прямоугольные геоцентрические координаты объекта (м) %nlo_xyz.x %nlo_xyz.y %nlo_xyz.z %Выходные данные: %структура satvis % top.s - проекция вектора дальности на ось, направленную на Юг (South) % top.e - проекция вектора дальности на ось, направленную на Восток (East) % top.z - проекция вектора дальности на ось, направленную в Зенит % top.daln - дальность до объекта % top.az - угол азимута объекта % top.el - угол видимости объекта rx = nlo_xyz.x - rec_xyz.x; ry = nlo_xyz.y - rec_xyz.y; rz = nlo_xyz.z - rec_xyz.z; r_sat = sqrt( rx * rx + ry * ry + rz * rz); r_rec = sqrt((rec_xyz.x)^2 + (rec_xyz.y)^2+ (rec_xyz.z)^2); top.r = r_sat; rx1 = rx; ry1 = ry; rz1 = rz; sin_lat = sin(rec_llh.lat); cos_lat = cos(rec_llh.lat); sin_lon = sin(rec_llh.lon); cos_lon = cos(rec_llh.lon); % Projections of vector of range in topocentric coordinate system: top.e = -sin_lon * rx1 + cos_lon * ry1; top.s = cos_lon * sin_lat * rx1 + sin_lon * sin_lat * ry1 - cos_lat * rz1; top.z = cos_lat * cos_lon * rx1 + cos_lat * sin_lon * ry1 + sin_lat * rz1; % azimut: отсчет по часовой стрелке от оси направленной на Север (N or -S) (-top.s) eps = 10e-10; if ( (abs(top.e) < eps) || (abs(top.s) < eps)) top.az = 0.0; else top.az = atan2(top.e,-top.s); end; if (top.az < 0.0) top.az = top.az + pi * 2; end; % elevation: cos_el_top = (rec_xyz.x * rx + rec_xyz.y * ry + rec_xyz.z * rz) / (r_sat * r_rec); if ( cos_el_top >= 1.00 ) 127 el = 0.0; else if ( cos_el_top <= -1.00 ) el = pi; else el = acos(cos_el_top); end; end; top.el = pi / 2.0 - el; satvis.el = top.el; satvis.az = top.az; satvis.r = r_sat; Функция rewrite_satpos function [result] = rewrite_satpos(nom, satpos) %Имя функции:rewrite_satpos %Функция предназначена для перезаписи структуры satpos в массив структур result result(nom).x = satpos.x; result(nom).y = satpos.y; result(nom).z = satpos.z; result(nom).vx = satpos.vx; result(nom).vy = satpos.vy; result(nom).vz = satpos.vz; Функция utc_nut function nut = utc_nut(t) %Имя: utc_nut %Функция предназначена для расчета нутации %Входные данные: %t =6.023472005475702e+002; %Выходные данные: %nut - нутация R = 1296000; % ( 1r=360grad=1 296 000 cek) RAD_SEK_ANGL = pi/(3600*180); t2 = t * t; t3 = t2 * t; l = 485866.733 + (1325.0 * R + 715922.633) * t + 31.310 * t2 + 0.064 * t3;%1.034807679476340e+012 l1 = 1287099.804 + (99 * R + 1292581.224) * t - 0.577 * t2 - 0.012 * t3; f = 335778.877 + (1342 * R + 295263.137) * t - 13.257 * t2 + 0.011 * t3; dd = 1072261.307 + (1236 * R + 1105601.328) * t - 6.891 * t2 + 0.019 * t3; omega = 450160.280 - (5 * R + 482890.539) * t + 7.455 * t2 + 0.008 * t3; eps0 = 84381.448 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3; % eps0 = 84381.447996 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3; 128 eps_d = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'd','e'); eps_k = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'k','e'); eps = eps0 + eps_d + eps_k; cos_eps = cos(RAD_SEK_ANGL * eps ) / 15.0; d_fi = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'd', 'f'); k_fi = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'k', 'f'); nut1 = d_fi * cos_eps; nut2 = k_fi * cos_eps; % nut3 = 0.00264 * sin(omega) + 0.000063 * sin(2.0 * omega) nut3 = 0; nut1_dop = nut1; nut2_dop = nut2; nut3_dop = nut3; nut = nut1 + nut2 + nut3; Функция utc_nut_fi_ep function nut_fi_eps = utc_nut_fi_eps(t, l, l1, f, dd, omega, typ_nut, fi_eps) % Имя функции:utc_nut_fi_eps %Функция предназначена для расчета нутации по епсилон и фи %Входные данные: %t - значение параметра приведено в функции utc_nut; %l - значение параметра приведено в функции utc_nut; %l1 - значение параметра приведено в функции utc_nut; %f -значение параметра приведено в функции utc_nut; %dd - значение параметра приведено в функции utc_nut; % оmega - значение параметра приведено в функции utc_nut; %typ_nut - признак параметра приведен в функции utc_nut; %fi_eps - признак параметра приведен в функции utc_nut; %Выходные данные: %nut_fi_eps - значение нутации по по епсилон и фи % применить функцию koef [koef_id, koef_abd, koef_ik, koef_abk] = koef; RAD_SEK_ANGL = pi/(3600*180); if (typ_nut == 'd') n = 30; else n = 76; end; sum_a = 0; sum_b = 0; for i = 1 : n if (typ_nut == 'd') 129 s1 = koef_id(i,1) * l + koef_id(i,2) * l1 + koef_id(i,3) * f + koef_id(i,4) * dd + koef_id(i,5) * omega; if (fi_eps == 'f') a = koef_abd(i,1) * 1e-4; bt = koef_abd(i,2) * 1e-4; else a = koef_abd(i,3) * 1e-4; bt = koef_abd(i,4) * 1e-4; end; else s1 = koef_ik(i,1) * l + koef_ik(i,2) * l1 + koef_ik(i,3) * f + koef_ik(i,4) * dd + koef_ik(i,5) * omega; if (fi_eps == 'f') a = koef_abk(i,1) * 1e-4; bt = koef_abk(i,2) * 1e-4; else a = koef_abk(i,3) * 1e-4; bt = koef_abk(i,4) * 1e-4; end; end; if (fi_eps == 'f') sin_s1 = sin(RAD_SEK_ANGL * s1); sa = a * sin_s1; sb = bt * sin_s1; else cos_s1 = cos(RAD_SEK_ANGL * s1); sa = a * cos_s1; sb = bt * cos_s1; end; arg = RAD_SEK_ANGL * s1; sum_a = sum_a + sa; sum_b = sum_b + sb; end; nut_fi_eps = sum_a + sum_b * t; Функция koef function[koef_id, koef_abd, koef_ik, koef_abk] = koef % Имя функции:koef % Функция предназначена для ввода коэффициентов при расчете нутации koef_id = [ 0, 0, 0, 0, 1 ; % 1 0, 0, 0, 0, 2 ; % 2 -2, 0, 2, 0, 1 ; % 3 2, 0,-2, 0, 0 ; % 4 -2, 0, 2, 0, 2 ; % 5 130 1,-1, 0,-1, 0 ; % 6 0,-2, 2,-2, 1 ; % 7 2, 0,-2, 0, 1 ; % 8 0, 0, 2,-2, 2 ; % 9 0, -1, 0, 0, 0 ; % 10 % исправлено 0, 1, 2,-2, 2 ; % 11 0,-1, 2,-2, 2 ; % 12 0, 0, 2,-2, 1 ; % 13 -2, 0, 0,2, 0 ; % 14 % исправлено 0, 0, 2,-2, 0 ; % 15 0, 2, 0, 0, 0 ; % 16 0, 1, 0, 0, 1 ; % 17 0, 2, 2,-2, 2 ; % 18 0,-1, 0, 0, 1 ; % 19 -2, 0, 0, 2, 1 ; % 20 0,-1, 2,-2, 1 ; % 21 2, 0, 0,-2, 1 ; % 22 0, 1, 2,-2, 1 ; % 23 1, 0, 0,-1, 0 ; % 24 2, 1, 0,-2, 0 ; % 25 0, 0,-2, 2, 1 ; % 26 0, 1,-2, 2, 0 ; % 27 0, 1, 0, 0, 2 ; % 28 -1, 0, 0, 1, 1 ; % 29 0, 1, 2,-2, 0 ];% 30 koef_abd = [ -171996.0,-174.2, 92025.0, 8.9; % 1 2062.0, 0.2, -895.0, 0.5; % 2 46.0, 0.0, -24.0, 0.0; % 3 11.0, 0.0, 0.0, 0.0; % 4 -3.0, 0.0, 1.0, 0.0; % 5 -3.0, 0.0, 0.0, 0.0; % 6 -2.0, 0.0, 1.0, 0.0; % 7 1.0, 0.0, 0.0, 0.0; % 8 -13187.0, -1.6, 5736.0,-3.1; % 9 -1426.0, 3.4, 54.0,-0.1; % 10 % исправлено -517.0, 1.2, 224.0,-0.6; % 11 217.0, -0.5, -95.0, 0.3; % 12 129.0, 0.1, -70.0, 0.0; % 13 -48.0, 0.0, 1.0, 0.0; % 14 % исправлено -22.0, 0.0, 0.0, 0.0; % 15 17.0, -0.1, 0.0, 0.0; % 16 131 -15.0, 0.0, 9.0, 0.0; % 17 -16.0, 0.1, 7.0, 0.0; % 18 -12.0, 0.0, 6.0, 0.0; % 19 -6.0, 0.0, 3.0, 0.0; % 20 -5.0, 0.0, 3.0, 0.0; % 21 4.0, 0.0, -2.0, 0.0; % 22 4.0, 0.0, -2.0, 0.0; % 23 -4.0, 0.0, 0.0, 0.0; % 24 1.0, 0.0, 0.0, 0.0; % 25 1.0, 0.0, 0.0, 0.0; % 26 -1.0, 0.0, 0.0, 0.0; % 27 1.0, 0.0, 0.0, 0.0; % 28 1.0, 0.0, 0.0, 0.0; % 29 -1.0, 0.0, 0.0, 0.0]; % 30 koef_ik = [ 0, 0, 2, 0, 2; % 31 1, 0, 0, 0, 0; % 32 0, 0, 2, 0, 1; % 33 1, 0, 2, 0, 2; % 34 1, 0, 0,-2, 0; % 35 -1, 0, 2, 0, 2; % 36 0, 0, 0, 2, 0; % 37 1, 0, 0, 0, 1; % 38 -1, 0, 0, 0, 1; % 39 -1, 0, 2, 2, 2; % 40 1, 0, 2, 0, 1; % 41 0, 0, 2, 2, 2; % 42 2, 0, 0, 0, 0; % 43 1, 0, 2,-2, 2; % 44 2, 0, 2, 0, 2; % 45 0, 0, 2, 0, 0; % 46 -1, 0, 2, 0, 1; % 47 -1, 0, 0, 2, 1; % 48 1, 0, 0,-2, 1; % 49 -1, 0, 2, 2, 1; % 50 1, 1, 0,-2, 0; % 51 0, 1, 2, 0, 2; % 52 0,-1, 2, 0, 2; % 53 1, 0, 2, 2, 2; % 54 1, 0, 0, 2, 0; % 55 2, 0, 2,-2, 2; % 56 0, 0, 0, 2, 1; % 57 0, 0, 2, 2, 1; % 58 132 1, 0, 2,-2, 1; % 59 0, 0, 0,-2, 1; % 60 1,-1, 0, 0, 0; % 61 2, 0, 2, 0, 1; % 62 0, 1, 0,-2, 0; % 63 1, 0,-2, 0, 0; % 64 0, 0, 0, 1, 0; % 65 1, 1, 0, 0, 0; % 66 1, 0, 2, 0, 0; % 67 1,-1, 2, 0, 2; % 68 -1,-1, 2, 2, 2; % 69 -2, 0, 0, 0, 1; % 70 3, 0, 2, 0, 2; % 71 0,-1, 2, 2, 2; % 72 1, 1, 2, 0, 2; % 73 -1, 0, 2,-2, 1; % 74 2, 0, 0, 0, 1; % 75 1, 0, 0, 0, 2; % 76 3, 0, 0, 0, 0; % 77 0, 0, 2, 1, 2; % 78 -1, 0, 0, 0, 2; % 79 1, 0, 0,-4, 0; % 80 -2, 0, 2, 2, 2; % 81 -1, 0, 2, 4, 2; % 82 2, 0, 0,-4, 0; % 83 1, 1, 2,-2, 2; % 84 1, 0, 2, 2, 1; % 85 -2, 0, 2, 4, 2; % 86 -1, 0, 4, 0, 2; % 87 1,-1, 0,-2, 0; % 88 2, 0, 2,-2, 1; % 89 2, 0, 2, 2, 2; % 90 1, 0, 0, 2, 1; % 91 0, 0, 4,-2, 2; % 92 3, 0, 2,-2, 2; % 93 1, 0, 2,-2, 0; % 94 0, 1, 2, 0, 1; % 95 -1,-1, 0, 2, 1; % 96 0, 0,-2, 0, 1; % 97 0, 0, 2,-1, 2; % 98 0, 1, 0, 2, 0; % 99 1, 0,-2,-2, 0; % 100 133 0,-1, 2, 0, 1; % 101 1, 1, 0,-2, 1; % 102 1, 0,-2, 2, 0; % 103 2, 0, 0, 2, 0; % 104 0, 0, 2, 4, 2; % 105 0, 1, 0, 1, 0]; % 106 koef_abk = [-2274.0, -0.2, 977.0, -0.5; % 31 712.0, 0.1, -7.0, 0.0; % 32 -386.0, -0.4, 200.0, 0.0; % 33 -301.0, 0.0, 129.0, -0.1; % 34 -158.0, 0.0, -1.0, 0.0; % 35 123.0, 0.0, -53.0, 0.0; % 36 63.0, 0.0, -2.0, 0.0; % 37 63.0, 0.1, -33.0, 0.0; % 38 -58.0, -0.1, 32.0, 0.0; % 39 -59.0, 0.0, 26.0, 0.0; % 40 -51.0, 0.0, 27.0, 0.0; % 41 -38.0, 0.0, 16.0, 0.0; % 42 29.0, 0.0, -1.0, 0.0; % 43 29.0, 0.0, -12.0, 0.0; % 44 -31.0, 0.0, 13.0, 0.0; % 45 26.0, 0.0, -1.0, 0.0; % 46 21.0, 0.0, -10.0, 0.0; % 47 16.0, 0.0, -8.0, 0.0; % 48 -13.0, 0.0, 7.0, 0.0; % 49 -10.0, 0.0, 5.0, 0.0; % 50 -7.0, 0.0, 0.0, 0.0; % 51 7.0, 0.0, -3.0, 0.0; % 52 -7.0, 0.0, 3.0, 0.0; % 53 -8.0, 0.0, 3.0, 0.0; % 54 6.0, 0.0, 0.0, 0.0; % 55 6.0, 0.0, -3.0, 0.0; % 56 -6.0, 0.0, 3.0, 0.0; % 57 -7.0, 0.0, 3.0, 0.0; % 58 6.0, 0.0, -3.0, 0.0; % 59 -5.0, 0.0, 3.0, 0.0; % 60 5.0, 0.0, 0.0, 0.0; % 61 -5.0, 0.0, 3.0, 0.0; % 62 -4.0, 0.0, 0.0, 0.0; % 63 4.0, 0.0, 0.0, 0.0; % 64 -4.0, 0.0, 0.0, 0.0; % 65 -3.0, 0.0, 0.0, 0.0; % 66 134 3.0, 0.0, 0.0, 0.0; % 67 -3.0, 0.0, 1.0, 0.0; % 68 -3.0, 0.0, 1.0, 0.0; % 69 -2.0, 0.0, 1.0, 0.0; % 70 -3.0, 0.0, 1.0, 0.0; % 71 -3.0, 0.0, 1.0, 0.0; % 72 2.0, 0.0, -1.0, 0.0; % 73 -2.0, 0.0, 1.0, 0.0; % 74 2.0, 0.0, -1.0, 0.0; % 75 -2.0, 0.0, 1.0, 0.0; % 76 2.0, 0.0, 0.0, 0.0; % 77 2.0, 0.0, -1.0, 0.0; % 78 1.0, 0.0, -1.0, 0.0; % 79 -1.0, 0.0, 0.0, 0.0; % 80 1.0, 0.0, -1.0, 0.0; % 81 -2.0, 0.0, 1.0, 0.0; % 82 -1.0, 0.0, 0.0, 0.0; % 83 1.0, 0.0, -1.0, 0.0; % 84 -1.0, 0.0, 1.0, 0.0; % 85 -1.0, 0.0, 1.0, 0.0; % 86 1.0, 0.0, 0.0, 0.0; % 87 1.0, 0.0, 0.0, 0.0; % 88 1.0, 0.0, -1.0, 0.0; % 89 -1.0, 0.0, 0.0, 0.0; % 90 -1.0, 0.0, 0.0, 0.0; % 91 1.0, 0.0, 0.0, 0.0; % 92 1.0, 0.0, 0.0, 0.0; % 93 -1.0, 0.0, 0.0, 0.0; % 94 1.0, 0.0, 0.0, 0.0; % 95 1.0, 0.0, 0.0, 0.0; % 96 -1.0, 0.0, 0.0, 0.0; % 97 -1.0, 0.0, 0.0, 0.0; % 98 -1.0, 0.0, 0.0, 0.0; % 99 -1.0, 0.0, 0.0, 0.0; % 100 -1.0, 0.0, 0.0, 0.0; % 101 -1.0, 0.0, 0.0, 0.0; % 102 -1.0, 0.0, 0.0, 0.0; % 103 1.0, 0.0, 0.0, 0.0; % 104 -1.0, 0.0, 0.0, 0.0; % 105 1.0, 0.0, 0.0, 0.0]; % 106 Функция semi_axis_1 135 % double semi_axis(double t_dr, double i_incl, double ecc, double omega_n) function [a_npp] = semi_axis_1(t_dr, i_incl, ecc, omega_n); %Имя функции:semi_axis_1 %Функция вычисляет радиус орбиты спутника ГЛОНАСС в соответствии с интерфейсным контроль- ным %документом ГЛОНАСС %Входные данные: %t_dr - драконический период обращения спутника ГЛОНАСС (секунды) %i_incl - наклонение орбиты спутника ГЛОНАСС (радиан) %ecc - эксцентриситет %omega_n - аргумент перигея орбиты спутника ГЛОНАСС (радиан) %Выходные данные: %a_npp - большая полуось орбиты спутника ГЛОНАСС (километр) 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; 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 136 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; Функция kepler function [Ek] = kepler(Mk, ecc, eps); %Имя функции:kepler %Функция предназначена для решения уравнения Кеплера % eps = 1.0E-15; y = ecc * sin(Mk); x1 = Mk; x = y; for k = 0 : 15 x2 = x1; x1 = x; y1 = y; y = Mk - (x - ecc * sin(x)); if (abs(y - y1) < eps) break end x = (x2 * y - x * y1) / (y - y1); end %k Ek = x; Функция map function map(N) %Имя функции:map %Применена функция MatLab для внесения в графики орбитального движения изображения Земли load('topo.mat','topo','topomap1'); [x,y,z] = sphere(50); cla reset %axis square off props.AmbientStrength = 0.1; 137 props.DiffuseStrength = 1; props.SpecularColorReflectance = .5; props.SpecularExponent = 20; props.SpecularStrength = 1; props.FaceColor= 'texture'; props.EdgeColor = 'none'; props.FaceLighting = 'phong'; props.Cdata = topo; surface(x*N,y*N,z*N,props); light('position',[-1 0 1]); light('position',[-1.5 0.5 -0.5], 'color', [.6 .2 .2]); view(3) ; grid on 4.3.5 Примеры выполнения комплекса программ ORBITA_GL_NAVIOR Некоторые результаты, иллюстрирующие работу программ, изображены на рис. 4.8 - рис. 4.11. Рис. 4.8. Орбиты спутников ГЛОНАСС в системе координат ECEF 138 Рис. 4.9. Орбиты спутников ГЛОНАСС в системе координат ECI Рис. 4.10. Орбита спутника ГЛОНАСС № 24 за 7 суток 139 0 24 68 10 12 14 16 18 20 22 24 02 46 810 12 14 16 18 20 22 24Время наблюдения спутников ГЛОНАСС, часНо ме р нав игац ионн ого сп ут ник аРис. 4.11. Видимость спутников ГЛОНАСС Программы имеют и другие возможности в зависимости от задач, которые может поставить специалист, исследующий орбитальное движение спутников. 140 1411   ...   4   5   6   7   8   9   10   11   ...   14

Функция 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);%закрытия файла для записи 1   ...   6   7   8   9   10   11   12   13   14

Цель лабораторной работы изучение орбитального движения навигационных спутников GPS, ГЛОНАСС, GALILEO с последующей визуализацией спутников на любой момент времени. 5.4.2 Лабораторная работа 5. 4 «Модель движения и визуализация спутников GPS, ГЛОНАСС, GALILEO » Комплекс программ для выполнения настоящей работы основан на пакете Vsion_GLONASS_GPS, вкотором файл Vision_GLONASS_GPS.m заменяется файлом Vision_GL_GPS_GAL.m, а входными данными является альманах спутников трех сис- тем- GPS, ГЛОНАСС, GALILEO приведенный в приложении. Рекомендуется следующий порядок выполнения лабораторной работы. 1. Создайте папку Vision_GLONASS_GPS_GALILEO _My. 2. Скопируйте в папку Vision_GLONASS_GPS_GALILEO _My изVsion_GLONASS_GPS функции и файлы, из папки Vi-sion_GLONASS_GPS_GALILEO файл Vision_GL_GPS_GAL.m, из приложения файл альманаха AlmGGG.yum. Комплекс готов к работе. Выполните задания 1, 2. 3. Задание 1 .Откройте файл Vision_GL_GPS_GAL.m, убедитесь в правильности введен- ных входных данных ( альманах , дата и время) исполните файл и результаты визуали- зации по каждой спутниковой системе внесите в отчет. 4. Задание 2. Измените в Vision_GL_GPS_GAL.m дату на 1 сутки, выполните файл, вне- сите результаты в отчет. Объясните причины, изменения местоположения спутников по каждой из систем. 5.4.3 Задание для самоподготовки Подготовьте сравнительную характеристику спутниковых радионавигационных систем GPS, ГЛОНАСС, GALILEO 186 5.4.4 Файл из папки Vision_GLONASS_GPS_GALILEO Файла Vision_GL_GPS_GAL.m clear all %Имя m-файла:Vision_GL_GPS_GAL.m %Программа рассчитывает углы видимости навигационных спутников GPS, ГЛОНАСС, GALLILEO %(азимута и места) спутников на заданный момент времени %Входные данные: %файл альманаха в формате Yuma,имя файла альманаха присваивается %переменной "Dat",например,Dat = 'имя файла альманаха'; %данные о начале отсчета "d2",d2='месяц/день/год';h=час;min=минута;s=секунда; %координаты позиции приемника -lat(широта в радианах),lon (долгота в радианах, %hr (высота в метрах); %шаг с каким будут рассчитываться параметры (step,секунды); %количество точек (L), в которых будут рассчитываться параметры орбит %L=12*3600/step,L читается так: количество часов (например,12) %число секунд в часе (3600) деленное на шаг (step) %Постоянные: %скорость вращения Земли OMEGAeDOT=7.2921151467e-005; %радиусы земного эллипсоида A_WGS84=6378137.0; B_WGS84=6356752.314; %константы mu=398600500000000; F_CONST = 4.442807633E-10; kt=1;%установка времени на титульной надписи графика, определяется параметрами d2'; h; min; s и j или L; %Задание цветов для графики color6(1:14) = ['k' 'b' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h']; %Входные данные Dat = 'AlmGGG.yum';%альманах спутников GPS, ГЛОНАСС, GALLILEO d2='12/21/2006'; h=14; min=52.0; s=58.0; %координаты точки, из которой наблюдаются спутники lat = 0.88032730015257;%50 градусов 26минут 20.54 секунд lon = 0.53109641675259;%30 градусов 25 минут 46.4995секунд hr=187.488;% метров X_label=['Широта' ':' num2str(lat) ';' 'долгота' ':' num2str(lon) ';' 'высота' ':' num2str(hr)]; 187 step=300;%шаг отсчета времени в секундах (300=5 минутам);шаг можно изменять step=0; %L=(24*3600) / step;% убрать % перед L для вывода таблицы улов видимости и азимута L=1;% установить % перед L для вывода таблицы улов видимости и азимута %Чтение альманаха [alm,max_kol] = Yuma_GPS_GLONAS_Alm(Dat); %max_kol = 32; nom = 1; i = 0; k = 0; i = 0; % for i = 1 : max_kol while ( i <91) id = alm(nom).ID; Health = alm(nom).Health; %fprintf('1: i=%i k=%i nom=%i id=%i Health=%i \n', i, k, nom, id, Health); if ( id > 0) Health = alm(nom).Health; if ( Health == 0) k = k + 1; nom_ns(k) = id; % fprintf('2: i=%i k=%i nom=%i id=%i Health=%i \n', i, k, nom, id, Health); nom = nom + 1; else nom = nom + 1; end; else nom = nom + 1; end; i = i + 1; end; % i kol = k; fprintf('kol=%i \n', kol); nom_ns % - номер навигационного спутника [Rx,Ry,Rz] = ECEFLLH(lon, lat,hr); %Rx=0;Ry=0;,Rz=0;%центр масс Земли %Начало отсчета текущего времени [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s); for j = 1% 0:L 188 % for j = 1:L % 0:L t( j )=weeks+step*(j); %-step; %t1(j) = t(j)/60; %изменение дискретности текущего времени %d_wn = (week - alm(i).Week); %d_wn = 0; for k = 1 : kol i = nom_ns(k) ; % input "i" !!! % fprintf('i=%i alm(i).Week= %i \n',i, alm(i).Week) d_wn=(modeweek - alm(i).Week);%если в альманахе не учтено 1024 tk = t(j) + d_wn * 604800 - alm(i).TOA; d_wn = abs(modeweek - alm(i).Week); dd = 302400.0 + d_wn * 604800; if ( ( alm(i).A05 > 0 ) & ( alm(i).Health == 0 ) ) while (abs(tk) > dd) if tk > dd tk = tk - 604800; else if tk < -dd tk = tk + 604800; end end % if end % while %Справочник по альманаху- цифра в скобках обозначает порядковый номер %параметра альманаха в формате YUMA %alm(ID).ID(1); alm(ID).Health(2); alm(ID).e(3); alm(ID).TOA(4); alm(ID).deltai(5); %alm(ID).OMEGADOT(6); alm(ID).A05(7); alm(ID).omega0(8); alm(ID).omega(9); %alm(ID).M0(10); alm(ID).Af0(11); alm(ID).Af1(12); alm(ID).Week(13); n0=sqrt((mu) / (alm(i).A05^6)); j2 = 1082.68E-6; re = (A_WGS84 + B_WGS84) / 2.; sin55 = sin(55.0 * pi / 180.0); dn = 1.5 * j2 * re * re / (alm(i).A05^4 ) * (1. - 1.5 * sin55 * sin55); %dn = 0; n=n0 * (1 + dn); 189 Mk = alm(i).M0 + n*tk; e=alm(i).e; %решение уравнения Кеплера eps = 1.0E-15; y = e * sin(Mk); x1 = Mk; x = y; for k = 0 : 15 % количество итераций x2 = x1; x1 = x; y1 = y; y = Mk - (x - e * sin(x)); if (abs(y - y1) < eps) break end x = (x2 * y - x * y1) / (y - y1); end % kepler Ek = x; deltr = F_CONST * alm(i).e * alm(i).A05 * sin(Ek); dt1 = alm(i).Af0 + alm(i).Af1 * tk + deltr; tk = tk - dt1; vd = 1. - alm(i).e * cos(Ek); nuk =atan2(sqrt(1-alm(i).e^2)*sin(Ek) / vd,(cos(Ek)-alm(i).e) / vd); Ek = acos((alm(i).e+cos(nuk))/(1+alm(i).e*cos(nuk))); Fk =nuk + alm(i).omega; uk =Fk; ik=alm(i).deltai; % alm(i) rk =(alm(i).A05^2)*(1.0-alm(i).e*cos(Ek)); xkk =rk*cos(uk); ykk =rk*sin(uk); OMEGAk =alm(i).omega0+(alm(i).OMEGADOT-OMEGAeDOT)*tk-OMEGAeDOT*alm(i).TOA; % fprintf('i=%i j=%i \n ',i, j); Xk(j,i) = xkk*cos(OMEGAk)-ykk*cos(ik)*sin(OMEGAk); Yk(j,i) = xkk*sin(OMEGAk)+ykk*cos(ik)*cos(OMEGAk); Zk(j,i) = ykk*sin(ik); %Дальности до спутников PR(j,i) = sqrt((Xk(j,i) - Rx)^2 + (Yk(j,i) - Ry)^2 + (Zk(j,i) - Rz)^2); 190 %Перевод в географическую систему если требуется %[lons,lats,hrs] = LLHECEF(Xk,Yk,Zk); %(Llon(j),Llat(j),Hhr(j)) = [lons,lats,hrs]; %расчет угла видимости спутника xls = Xk(j,i) - Rx; yls = Yk(j,i) - Ry; zls = Zk(j,i) - Rz; range1 = sqrt(xls*xls+yls*yls+zls*zls); if j>1 doppler(j-1) = (range1 - range2) * 5.2514 / step; end range2 = range1; P = sqrt(Rx * Rx + Ry * Ry + Rz * Rz); tdot = ( Rx*xls+Ry*yls+Rz*zls)/range1/P; xll = xls/range1; yll = yls/range1; zll = zls/range1; if tdot >= 1.00 b = 0.0; elseif tdot <= -1.00 b = pi; else b = acos( tdot); end satang = pi/2.0 - b; TT =satang; TT(j,i) =TT;%угол видимости спутников %расчет угла азимута спутников xn =-cos(lon)*sin(lat); yn =-sin(lon)*sin(lat); zn = cos(lat); xe =-sin(lon); ye = cos(lon); xaz = xe*xll + ye*yll; yaz = xn*xll + yn*yll + zn*zll; if (xaz == 0) or (yaz == 0) az(j)= 0; else az(j,i) = atan2(xaz,yaz); end 191 if az(j,i) < 0 az(j,i) = az(j,i) + pi*2; end AZ(j,i) =az(j,i) *180/pi;%угла азимута спутников в градусах EL(j,i) = TT(j,i) *180/pi;%угла видимости спутников в градусах % ПЕРЕСЧЕТ ВРЕМЕНИ A(j)=mod(t(j),86400); her(j)=floor(A(j)/3600); m(j)=floor((A(j)-her(j)*3600)/60); sek(j)=A(j)-her(j)*3600-m(j)*60; %Построение полярной системы координат if EL(j,i) < 0 elp = 180; else elp = (EL(j,i) - 90); end; azp = (AZ(j,i) + 90.0); rad = pi / 180; x0 = 0; y0 = 0; xt(j,i) = (elp * cos(azp * rad)); yt(j,i) = -(elp * sin(azp * rad)); % fprintf('i=%i j=%i \n' , i, j); end % i = ns end; % if ( alm(i).A05 > 0 ) end % j = time %ВНИМАНИЕ. Для вывода времени визуализации спутников на график установите kt t_itle=[d2 ' ' num2str(her(kt)) ':' num2str(m(kt)) ':' num2str(sek(kt))]; %X_label=['Широта' ':' num2str(lat) ';' 'долгота' ':' num2str(lon) ';' 'высота' ':' num2str(hr)]; %num2ctr(lat) %num2str(her(kt)) %X_label=['66' ':']; n = 6; max_n = max(nom_ns); n_end = mod(max(nom_ns),n); n_end = mod(kol, n); n2 = fix(kol / n) * n - n +1; 192 %Формирование таблицы вывода времени UTC (Time), GPS (Tgps в секундах), номера спутника (nns), % углов видимости и азимута от времени и номера спутника for i=1:n:kol fprintf(' Time Tgps '); for k=1: n nns = nom_ns(i+k-1); fprintf(' %2i ', nns); end; fprintf(' \n'); for j=1:L fprintf('%2i:%2i:%2i %i ',her(j),m(j),sek(j), t(j)); for k=1: n nns = nom_ns(i+k-1); fprintf('%6.1f *%6.1f ', EL(j,nns), AZ(j,nns) ); end; fprintf(' \n'); end ; % J=1:L if (i) == (n2) n = n_end; end; end% i hold on %Окружности уровней на круговой диаграмме видимости спутников k1 = 10; k2 = 30; k3 = 50; k4 = 70; k5 = 85; k6=90; n=0; for k=1:5:365 n=n+1; m1 = pi / 180; x(n)=cos(k*m1); y(n)=sin(k*m1); end; %График круговой диаграммы plot(k1*x(:),k1*y(:),'k:', k2*x(:),k2*y(:),'k:', k3*x(:),k3*y(:),'k:',k4*x(:),k4*y(:),'k:', k5*x(:),k5*y(:),'r', k6*x(:),k6*y(:),'r:'); text(5, 10,'80','FontSize',12,'FontName','TimesNewRoman'); 193 text(18, 23,'60','FontSize',12,'FontName','TimesNewRoman'); text(32, 37,'40','FontSize',12,'FontName','TimesNewRoman'); text(45, 50,'20','FontSize',12,'FontName','TimesNewRoman'); text(55, 60,'5','FontSize',12,'FontName','TimesNewRoman'); text(62, 67,'0','FontSize',12,'FontName','TimesNewRoman'); grid on; %Построение изображений видимых спутников на круговой диаграмме i=1; for k=1:kol i = nom_ns(k) ; if (i < 31) plot(xt(kt,i),yt(kt,i), 'Marker' ,'o','MarkerSize',20,'MarkerFaceColor','g' ) ; end; if ((i > 31) & ( i < 62)) plot(xt(kt,i),yt(kt,i), 'Marker' ,'o','MarkerSize',20,'MarkerFaceColor','r' ) ; end; if (i > 61) plot(xt(kt,i),yt(kt,i), 'Marker' ,'o','MarkerSize',20,'MarkerFaceColor','y' ) ; end; title(t_itle); xlabel(X_label,'FontSize',12,'FontName','TimesNewRoman') set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman') hold on str1 = num2str( i, 2); %{ if i<32 str1 = num2str( i, 2); else str1 = num2str( (i), 2); end %} text(xt(kt,i), yt(kt,i),str1,'FontSize',14,'FontName','TimesNewRoman','HorizontalAlignment','center' ); axis( [-100 100 -100 100]); %axis( [-90 90 -90 90]); end clear Результат выполнения программы изображен на рис. 5.5. 194 -100-50 050 100-100-80-60-40-20 020 40 60 80 100 80 60 40 20 50 12/21/2006 14:52:58Широта:0.88033;долгота:0.5311;высота:187.488 37 14 16 18 19 21 22 26 29 38 44 45 59 61 64 65 66 72 73 74 75 80 81 82 83 91Рис. 5.5. Видимость спутников GPS, ГЛОНАСС, GALILEO 5.4.5 Лабораторная работа 5. 5 «Орбиты спутников GPS, ГЛОНАСС, GALILEO » Цель лабораторной работы моделирование орбит спутников GPS, ГЛОНАСС, GALILEOДля одновременного расчета орбит спутников трех систем GPS, ГЛОНАСС, GALI-LEO пакет программ, сформированный в лабораторной работе 5. 4 может быть дополнен файлами ORBITA_GGG.m и map. m из папки ORBIT_GGG. Рекомендуется следующий порядок выполнения лабораторной работы. 1. Создайте папку ORBIT_GGG _My. 2. Из папки Vision_GLONASS_GPS_GALILEO скопируйте в папку ORBIT_GGG _My функции ECEFLLH, LLHECEF, Tim; файлы Yuma_GPS_GLONAS_Alm.m, AlmGGG.yum; из папки ORBIT_GGG файл ORBITA_GGG.m. 3. Запустите MatLab [7, 8]. 4. Обратитесь к папке ORBIT_GGG _My , откройте ее, изучите функции, файлы и про- граммные процедуры по комментариям и выполните задания 1, 2. 5. Задание 1. Откройте файл ORBITA_GGG.m. Сформируйте желаемую конфигурацию спутников, орбиты которых требуется рассчитать (см. комментарий к файлу). Устано- вите скорость вращения Земли, равную 0 (OMEGAeDOT=0). Выполните файл. Дайте описание, полученного графического изображения орбит. Результаты внесите в отчет. 6. Задание 2. В файле ORBITA_GGG.m сформируйте конфигурацию спутников: 1 спут- ник GPS, 1 спутник ГЛОНАСС, 1 спутник GALILEO. Установите 195 OMEGAeDOT=7.2921151467e-005. Выполните файл. Дайте описание, полученного графического изображения орбит. Результаты внесите в отчет. 5.4.6 Задания и вопросы для самоподготовки 1. Чем различаются орбиты спутников GPS, ГЛОНАСС, GALILEO? 2. По каким параметрам можно определить принадлежность спутника той или иной ор- бите? 3. Сформулируйте алгоритм определения дальности до спутника из точки наблюдения. 4. В какой системе координат рассчитываются углы видимости спутника? 5.4.7 Листинг файла ORBITA_GGG.m Файл ORBITA_GGG.m clear all %Имя m-файла:ORBITA_GGG.m %Программа рассчитывает орбиты спутников GPS, ГЛОНАСС, GALLILEO %Входные данные: %файл альманаха в формате Yuma,имя файла альманаха присваивается %переменной "Dat",например,Dat = 'имя файла альманаха'; %данные о начале отсчета "d2",d2='месяц/день/год';h=час;min=минута;s=секунда; %координаты позиции приемника -lat(широта в радианах),lon (долгота в радианах, %hr (высота в метрах); %шаг с каким будут рассчитываться параметры (step,секунды); %количество точек (L), в которых будут рассчитываться параметры орбит %L=12*3600/step,L читается так: количество часов (например,12) %число секунд в часе (3600) деленное на шаг (step) %Выходные данные: %Xk(j,i) - координата X спутника с номером i на момент времени j; %Yk(j,i) - координата Y спутника с номером i на момент времени j; %Zk(j,i) - координата Z спутника с номером i на момент времени j; %Примечание: входные и выходные данные могут быть дополнены / изменены закомментированны- ми %параметрами, помеченными строкой %%%%%%%%%%%%%%%%%%%%%%%; %Постоянные: %скорость вращения Земли %OMEGAeDOT=7.2921151467e-005; OMEGAeDOT=0;%%%%%%%%%%%%%%для орбитальных плоскостей %радиусы земного эллипсоида A_WGS84=6378137.0; B_WGS84=6356752.314; %константы mu=398600500000000; 196 F_CONST = 4.442807633E-10; %Задание цветов для графики j_color = 0; color6(1:11) = ['B' 'B' 'B' 'B' 'B' 'B' 'o' 'o' '+' '+' '+'];%%%%%%%%%%%%%%%%%%%%%%% %Входные данные Dat = 'AlmGGG.yum';%альманах спутников GPS, ГЛОНАСС, GALLILEO N=6378136;% радиус Земли (используется, как нормирующий коэффициент map(N);%функция выводит на графики Землю d2='12/21/2006'; %%%%%%%%%%%%%%%%%%%%%%% h=14; min=52.0; s=58.0;%%%%%%%%%%%%%%%%%%%%%%% %координаты точки, из которой наблюдаются спутники lat = 0.88032730015257;%50 градусов 26минут 20.54 секунд %%%%%%%%%%%%%%%%%%%%%% lon = 0.53109641675259;%30 градусов 25 минут 46.4995секунд %%%%%%%%%%%%%%%%%%%%% hr=187.488;% метров %%%%%%%%%%%%%%%%%%%%%%% step=600;%шаг отсчета времени в секундах (300=5 минутам);шаг можно изменять %step=0; L=(12*3600) / step;% убрать перед L для вывода таблицы улов видимости и азимута %%%%%%%%%% %L=1;% установить % перед L для вывода таблицы улов видимости и азимута %Чтение альманаха [alm,max_kol] = Yuma_GPS_GLONAS_Alm(Dat);%%%%%%%%%%%%%%%%%%%%%%% nom = 1; i = 0; k = 0; i = 0; while ( i <91) id = alm(nom).ID; Health = alm(nom).Health; if ( id > 0) Health = alm(nom).Health; if ( Health == 0) k = k + 1; nom_ns(k) = id; nom = nom + 1; else nom = nom + 1; end; else 197 nom = nom + 1; end; i = i + 1; end; % i kol = k; nom_ns ; % - номер навигационного спутника [Rx,Ry,Rz] = ECEFLLH(lon, lat,hr);%%%%%%%%%%%%%%%%%%%%%%% [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s); for j = 1:L % for j = 1:L % 0:L t( j )=weeks+step*(j); %-step;%%%%%%%%%%%%%%%%%%%%%%% %t1(j) = t(j)/60; %изменение дискретности текущего времени %d_wn = (week - alm(i).Week); %d_wn = 0; kol=11;%%%%%%%%%%%%%%%%%%%%%%% nom_ns(1:kol) =[1 3 4 5 9 10 38 58 62 72 82 ];%%%%%%%%%%%%%%%%%%%%%%% for k = 1 : kol i = nom_ns(k) ; % input "i" !!! d_wn=(modeweek - alm(i).Week);%если в альманахе не учтено 1024 tk = t(j) + d_wn * 604800 - alm(i).TOA; d_wn = abs(modeweek - alm(i).Week); dd = 302400.0 + d_wn * 604800; if ( ( alm(i).A05 > 0 ) & ( alm(i).Health == 0 ) ) while (abs(tk) > dd) if tk > dd tk = tk - 604800; else if tk < -dd tk = tk + 604800; end end % if end % while %Справочник по альманаху- цифра в скобках обозначает порядковый номер %параметра альманаха в формате YUMA %alm(ID).ID(1); alm(ID).Health(2); alm(ID).e(3); alm(ID).TOA(4); alm(ID).deltai(5); %alm(ID).OMEGADOT(6); alm(ID).A05(7); alm(ID).omega0(8); alm(ID).omega(9); %alm(ID).M0(10); alm(ID).Af0(11); alm(ID).Af1(12); alm(ID).Week(13); n0=sqrt((mu) / (alm(i).A05^6)); j2 = 1082.68E-6; re = (A_WGS84 + B_WGS84) / 2.; sin55 = sin(55.0 * pi / 180.0); dn = 1.5 * j2 * re * re / (alm(i).A05^4 ) * (1. - 1.5 * sin55 * sin55); 198 %dn = 0; n=n0 * (1 + dn); Mk = alm(i).M0 + n*tk; e=alm(i).e; %решение уравнения Кеплера eps = 1.0E-15; y = e * sin(Mk); x1 = Mk; x = y; for k = 0 : 15 % количество итераций x2 = x1; x1 = x; y1 = y; y = Mk - (x - e * sin(x)); if (abs(y - y1) < eps) break end x = (x2 * y - x * y1) / (y - y1); end % kepler Ek = x; deltr = F_CONST * alm(i).e * alm(i).A05 * sin(Ek); dt1 = alm(i).Af0 + alm(i).Af1 * tk + deltr; tk = tk - dt1; vd = 1. - alm(i).e * cos(Ek); nuk =atan2(sqrt(1-alm(i).e^2)*sin(Ek) / vd,(cos(Ek)-alm(i).e) / vd); Ek = acos((alm(i).e+cos(nuk))/(1+alm(i).e*cos(nuk))); Fk =nuk + alm(i).omega; uk =Fk; ik=alm(i).deltai; rk =(alm(i).A05^2)*(1.0-alm(i).e*cos(Ek)); xkk =rk*cos(uk); ykk =rk*sin(uk); OMEGAk =alm(i).omega0+(alm(i).OMEGADOT-OMEGAeDOT)*tk-OMEGAeDOT*alm(i).TOA; Xk(j,i) = xkk*cos(OMEGAk)-ykk*cos(ik)*sin(OMEGAk); Yk(j,i) = xkk*sin(OMEGAk)+ykk*cos(ik)*cos(OMEGAk); Zk(j,i) = ykk*sin(ik); %%%%%%%%%%%%%%%%%%%%%%% %{ %Дальности до спутников PR(j,i) = sqrt((Xk(j,i) - Rx)^2 + (Yk(j,i) - Ry)^2 + (Zk(j,i) - Rz)^2); %Перевод в географическую систему если требуется %[lons,lats,hrs] = LLHECEF(Xk,Yk,Zk); 199 %(Llon(j),Llat(j),Hhr(j)) = [lons,lats,hrs]; %расчет угла видимости спутника xls = Xk(j,i) - Rx; yls = Yk(j,i) - Ry; zls = Zk(j,i) - Rz; range1 = sqrt(xls*xls+yls*yls+zls*zls); if j>1 doppler(j-1) = (range1 - range2) * 5.2514 / step; end range2 = range1; P = sqrt(Rx * Rx + Ry * Ry + Rz * Rz); tdot = ( Rx*xls+Ry*yls+Rz*zls)/range1/P; xll = xls/range1; yll = yls/range1; zll = zls/range1; if tdot >= 1.00 b = 0.0; elseif tdot <= -1.00 b = pi; else b = acos( tdot); end satang = pi/2.0 - b; TT =satang; TT(j,i) =TT;%угол видимости спутников %расчет угла азимута спутников xn =-cos(lon)*sin(lat); yn =-sin(lon)*sin(lat); zn = cos(lat); xe =-sin(lon); ye = cos(lon); xaz = xe*xll + ye*yll; yaz = xn*xll + yn*yll + zn*zll; if (xaz == 0) or (yaz == 0) az(j)= 0; else az(j,i) = atan2(xaz,yaz); end if az(j,i) < 0 az(j,i) = az(j,i) + pi*2; 200 end AZ(j,i) =az(j,i) *180/pi;%угла азимута спутников в градусах EL(j,i) = TT(j,i) *180/pi;%угла видимости спутников в градусах % ПЕРЕСЧЕТ ВРЕМЕНИ A(j)=mod(t(j),86400); her(j)=floor(A(j)/3600); m(j)=floor((A(j)-her(j)*3600)/60); sek(j)=A(j)-her(j)*3600-m(j)*60; %Построение полярной системы координат if EL(j,i) < 0 elp = 180; else elp = (EL(j,i) - 90); end; azp = (AZ(j,i) + 90.0); rad = pi / 180; x0 = 0; y0 = 0; xt(j,i) = (elp * cos(azp * rad)); yt(j,i) = -(elp * sin(azp * rad)); %} %%%%%%%%%%%%%%%%%%%%%%%% end % i = ns end; % if ( alm(i).A05 > 0 ) end % j = time for (i=1:kol) j_color = j_color + 1; if (j_color > 11 ) j_color = 1; end S = color6(j_color); prn = nom_ns(i); %%%%%%%%%%%%%%%%%%%%%%% hold on h_F1 = gca; plot3(Xk(:,prn),Yk(:,prn),Zk(:,prn),S,'LineWidth',0.5); axis([ -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7)]); set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); set(h_F1,'Position',[0.1 0.1 0.85 0.9]) ; xlabel('Координата X') 201 ylabel('Координата Y'), zlabel('Координата Z'),grid on str1 = num2str( prn); text(Xk(j,prn), Yk(j,prn),Zk(j,prn),str1,'FontSize',14,'FontName','TimesNewRoman','HorizontalAlignment','center' ); hold on %%%%%%%%%%%%%%%%%%%%%%% end clear 5.4.8 Пример выполнения файла:ORBITA_GGG.m Результат выполнения программы изображен на рис. 5.6. Рис. 5.6. Орбиты спутников GPS, ГЛОНАСС, GALILEO На рис. 5.6 сплошными линиями показаны 6 орбит спутников GPS, маркером «о» 2 орбиты ГЛОНАСС, маркером «+» 3 орбиты GALILEO; цифрами обозначены номера спутников соответствующих систем. 5.5 Декодирование и расшифровка данных спутников ГЛОНАСС 5.5.1 Краткие сведения из теории В данном подразделе приводится программа декодирования данных спутников ГЛОНАСС. Данные, которые декодируются получены экспериментально с реального на-202 вигационного приемника после корреляционной обработки в виде последовательности символов «0» и «1». Поскольку объем данных относительной большой около 60000 сим- волов, то приводятся два вспомогательных m-файла. Эти вспомогательные программы по- зволяют сформировать данные на бумажном носителе в виде последовательности шестна- дцатеричных чисел, из которых можно восстановить m-файл для декодирования. Экспе- риментальные данные рассматриваемые далее представляют собой символы длительно- стью 10 миллисекунд и являются составляющими информационных символов сложенных по модулю 2 с меандром имеющим частоту 100 Гц. Примененный алгоритм декодирования основан на строгом соответствии данных, передаваемых с навигационного спутника ГЛОНАСС, которые соответствуют формату сообщений. Более подробно формат сообщений и данные, передаваемые в кадрах сообще- ний приводятся в интерфейсном контрольном документе ГЛОНАССи в книге [1] под- раздел 2. 2, стр. 77- 87). Комплекс программ декодирования и расшифровки находится в папке Decod_GLONASS. Цель лабораторной работы освоение и применение методики декодирования и расшифровки данных спутниковой системы ГЛОНАСС. 5.5.2 Лабораторная работа 5. 6 «Декодирование данных спутников ГЛОНАСС в на- вигационном приемнике» Рекомендуется следующий порядок выполнения лабораторной работы. 1. Создайте папку Decod_GLONASS _My и скопируете в ее файлы и функции из папки Decod_GLONASS. 2. Запустите MatLab. 3. Обратитесь к папке Decod_GLONASS _My , откройте ее, изучите функции, файлы и программные процедуры по комментариям и выполните задания 1, 2. 4. Задание 1. Откройте файл Hex_GLONASS.m. Изучите программные процедуры фай- ла, используя его текст и расширенный комментарий.Файл преобразования экспериментальных данных в шестнадцатеричные симво-лы Hex_GLONASS.m %Имя m- файла: Hex_GLONASS.m %Файл Hex_GLONASS.m предназначен для преобразования экспериментальных данных, записанных %в фйле с именем raw_bits_SV18J.dat в файл SV1.dat, записанный в шестнадцатеричных символах name_in = 'raw_bits_SV18J.dat'; % имя экспериментального файла fw = fopen('SV1.dat','Wt');%имя файла в шестнадцатеричных символах M=dlmread(name_in);%чтение экспериментального файла в матрицу M %чтение двоичных символов в 1875 строк по 16 символов в строке [ ](203 for i=1:1875 str_danD(i,:)=sprintf('%d',M(i+15*(i-1):(15+i+15*(i-1)))); end bin_dec=bin2dec(str_danD);%преобразование строк символов в десятичные числа dec_hex=dec2base(bin_dec,16)';%преобразование десятичных чисел в шестнадцатеричные fprintf(fw,'%c',dec_hex);%запись файла в шестнадцатеричном представлении fclose(fw);%закрытие записанного файла Расширенный комментарий к файлу Hex_GLONASS.m.При условии, что имеется файл данных спутника с именем raw_bits_SV18J.dat ре- зультат выполнения m- файла: Hex_GLONASS.m, записанный под именем SV1.dat и про- читанный клавишей F3 файла в Total Commander в шестнадцатеричном представлении имеет вид AAD34B52CAAB3534B52CAB5332ACD4AD3F1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CB52D3555334ACAAD4B55552B32CACAAB4AB2B34D3F1BA84B2CCB4B54D52AAAB2B2CAACB535554AB4ACD54B52CB3F1BA84B354AD3554ACCB2D2AAAACAAAD4CB2D35554ACAD4D2BF1BA84B352CCCCD4CD52D534D534CB535554D53532AD2AB34BF1BA84B34B52CD52B4CD52AD54B52AAD4CB34CAAD4B352D2B3F1BA84B34D32AAACCD4B5532CB5334ACAAAAD2D532AD34AACBF1BA84B32B5354AAB4CD2B2CAB552AAD534B35555532CD2D4BF1BA84B32CAAB2B2AAB32B55554CCB535554D535352ACAB4ABF1BA84B33555555555 5555555555555555555555555554D2ABF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534B54B4B34CB5334CAB54D2CB2AD32D2D2B33F1BA84B2B554ACB554AB3 534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B32AAAB55552CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAAAD55555555555555555555555555555554D4B3F1BA84B2CD555555555555555555555555555555555554D2B3F1BA84B355554AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB534BF1BA84B352AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB3F1BA84B3 4B554CACAAB4D2D4B4AD3552D4D32D54AD4B54B4ABF1BA84B34CCB34CD33552B2AAACD34AD3333 552D4CACAACCCBF1BA84B32AAAD555555555555555555555555555555554D4B3F1BA84B32D5555555555 555555555555555555555555552D53F1BA84B334AAD34B2AAB52D4CCB32AAD2B2D55554D2D2D54D3F1BA84B332CACB534AAB2CCACAAD34AD3333552D4D2AD4ACD3F1BA84B2AD2B534B534B4CB34ACCB35 4AB2D34D52CD2D332ABF1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B32AAAB55552CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD555555555555 555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAACAAAD32AD3554D2CCAAD2B2ACAAAD4D2B33333F1BA84B2CD4CCAD32D4AACCAB4B2CB52CCCD334D4D2D4AB4CBF1BA84B355553354B4D552CCB54ACAACAAD333554B4CB4CAD3F1BA84B352D4B4B32D4D34B532D4CB52CCCD4B52B2D4CAB54BF1BA84B34AAB554B55354ACAD3553553554CACAACAACACD2ABF1BA84B34D4B2B32B52B4D55354D34AD3332AAD2B2D2AACD4BF1BA84B32AAB52AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB554BF1BA84B32D5555555555555555555555555555555555552D53F1BA84B33554B5555555555555555555555555555555552D4BF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534B4B4B4CB34ACCB354AB2D34D52CD2D354CBF1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B32AAAB55552CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAB4D55555555555555555555555555555555 2B53F1BA84B2CD555555555555555555555555555555555554D2B3F1BA84B35554D555555555555555555555 5555555555552D4BF1BA84B352AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB3F1BA84B34AAB2D55555555555555555555555555555554D4B3F1BA84B34D555555555555555555555555555555555555 2D4BF1BA84B32AAB35555555555555555555555555555555552B53F1BA84B32D555555555555555555555555 5555555555552D53F1BA84B33554CD55555555555555555555555555555554D2ABF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534B4CB4B34CB5334CAB54D2CB2AD32D2D4D53F1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B32AAAB555 52CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B204 2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAD55555555555555 555555555555555555552B33F1BA84B2CD555555555555555555555555555555555554D2B3F1BA84B354AD52AAAD4CAB54CD54AAAD534AB4AAACB2B2B533F1BA84B352ACD2B52D2D5554AAACCB535554B2D53 52B4B4B53F1BA84B34B52B4AB4D332D54CAD34AAD52CD4D52AAD34D4CB3F1BA84B34CCAAD554ACACCACCD3334AD3334ACCAAB2D2AAD4BF1BA84B32AAD4D555555555555555555555555555555552B2BF1BA84B32D5555555555555555555555555555555555552D53F1BA84B33552D55555555555555555555555555 5555554D2CBF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534B2B4B4CB34ACCB354AB2D34D52CD2D354D3F1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B2D5554AAAAD34B52CAAB3534B52CAB5332ACD4AD3F1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CB52D3555334ACAAD4B55552B32CACAAB4AB2B34D3F1BA84B2CCB4B54D52AAAB2B2CAACB535554AB4ACD54B52CB3F1BA84B354AD3554ACCB2D2AAAACAAAD4CB2D35554ACAD4D2BF1BA84B352CCCCD4CD52D534D534CB535554D53532AD2AB34BF1BA84B34B52CD52B4CD52AD54B52AAD4CB34CAAD4B352D2B3F1BA84B34D32AAACCD4B5532CB5334ACAAAAD2D532AD34AACBF1BA8 4B32B5354AAB4CD2B2CAB552AAD534B35555532CD2D4BF1BA84B32CAAB2B2AAB32B55554CCB53555 4D535352ACAB4ABF1BA84B335555555555555555555555555555555555554D2ABF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534B2CB4B34CB5334CAB54D2CB2AD32D2D4D4BF1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B32AAAB555 52CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAAAD55555555555 555555555555555555554D4B3F1BA84B2CD555555555555555555555555555555555554D2B3F1BA84B355554AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB534BF1BA84B352AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB3F1BA84B34B554CACAAB4D2D4B4AD3552D4D32D54AD4B54B4ABF1BA84B34CCB34CD33552B2AAACD34AD3333552D4CACAACCCBF1BA84B32AAAD555555555555555555555555 555555554D4B3F1BA84B32D5555555555555555555555555555555555552D53F1BA84B334AAD34B2AAB52D4CCB32AAD2B2D55554D2D2D54D3F1BA84B332CACB534AAB2CCACAAD34AD3333552D4D2AD4ACD3F1BA84B2AD2B534B34B4B34CB5334CAB54D2CB2AD32D2D2B2BF1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B32AAAB55552CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAACAAAD32AD3554D2CCAAD2B2ACAAAD4D2B33333F1BA84B2CD4CCAD32D4AACCAB4B2CB52CCCD334D4D2D4AB4CBF1BA84B355553354B4D552CCB54ACAACAAD333 554B4CB4CAD3F1BA84B352D4B4B32D4D34B532D4CB52CCCD4B52B2D4CAB54BF1BA84B34AAB554B5 5354ACAD3553553554CACAACAACACD2ABF1BA84B34D4B2B32B52B4D55354D34AD3332AAD2B2D2AACD4BF1BA84B32AAB52AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB554BF1BA84B32D55555555 55555555555555555555555555552D53F1BA84B33554B5555555555555555555555555555555552D4BF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534B334B4CB34ACCB354AB2D34D52CD2D332B3F1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA8 4B2B32AAAB55552CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD5555555 55555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAB4D555555555555555555555555555555552B53F1BA84B2CD555555555555555555555555555555555554D2B3F1BA84B35554D5555555555555555555555555555555552D4BF1BA84B352AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAB3F1BA84B34AAB2D55555555555555555555555555555554D4B3F1BA84B34D5555555555555555555555555555555555552D4BF1BA84B32AAB35555555555555555555555555555555552B53F1BA84B32D5555555555555555555555555555555555552D53F1BA84B33554CD55555555555555555555555 555555554D2ABF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534D54B4B34CB5334CAB54D2CB2AD32D2CD4ABF1BA84B2B554ACB554AB3534B4D4D4AAACAB4ACAD4B54D32D3F1BA84B2B32AAAB55552CB4AD3554CACB4AD354ACCD532CB4BF1BA84B2D4AB4D4B54ACD54AAD555555555555555555552B2BF1BA84B2D34ACB4D5555555554B54AAAAAAAAAAAAAAAAACB33F1BA84B2CAAD55555555555555555555555555555555552B33F1BA84B2CD5555555555555555555 55555555555555554D2B3F1BA84B354AD52AAAD4CAB54CD54AAAD534AB4AAACB2B2B533F1BA84B35 2ACD2B52D2D5554AAACCB535554B2D5352B4B4B53F1BA84B34B52B4AB4D332D54CAD34AAD52CD4D52AAD34D4CB3F1BA84B34CCAAD554ACACCACCD3334AD3334ACCAAB2D2AAD4BF1BA84B32AAD4D555555555555555555555555555555552B2BF1BA84B32D5555555555555555555555555555555555552D53F1BA84B33552D555555555555555555555555555555554D2CBF1BA84B332AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABF1BA84B2AD2B534D534B4CB34ACCB354AB2D34D52CD2D2CD33F1BA84B2B554ACB4000000000000000000000000000000000000000000000000000 Этот результат можно было получить при условии, что имеется эксперименталь- ный файл raw_bits_SV18J.dat. Если экспериментального файла нет, то из его шестнадцате- ричного представления, приведенного выше можно сформировать файл для последующе-205 го декодирования. Для этого скопируем приведенное шестнадцатеричное представление файла и запишем его в виде файла с именем FilGL1.txt с помощью текстового редактора. В форму для декодирования (эта форма в точности соответствует файлу raw_bits_SV18J.dat файл FilGL1.txt преобразуется с помощью программы HexBin_GLONASS.m. Файл HexBin_GLONASS.m %Имя m- файла: HexBin_GLONASS.m %Файл HexBin_GLONASS предназначен для преобразования файла FilGL1.txt %в шестнадцатеричных символах в файл FilGL1_bin.txt в символах "0" и "1" name_in1='FilGL1.txt';%имя файла в шестнадцатеричном представлении fid = fopen(name_in1,'rt'); %открытие файла fw_bin=fopen('FilGL1_bin.txt','wb');%запись преобразованного файла %чтение и переформатирование файла MAS = fscanf(fid,'%s'); for i=1:1875 str_MAS(1:4,i)=sprintf('%s',MAS(i+3*(i-1):(3+i+3*(i-1)))); end MATR=reshape(MAS,4,1875)'; hex_dec= base2dec(MATR,16);%преобразование шестнадцатеричных чисел в десятичные dec_bin= dec2bin(hex_dec)';%преобразование десятичных чисел в двоичные MATR_Char = reshape(dec_bin,1,1875*16);%преобразование в строку %fprintf(fw_bin,'%c ',MATR_Char);%запись файла в символах "0" и "1"(вариант, пробел после с имеет значение) fprintf(fw_bin,'%c\n',MATR_Char);%запись файла в символах "0" и "1" (вариант) fclose(fw_bin);% закрытие файла Файлы FilGL1.txt (восстановленный из шестнадцетиричного) и raw_bits_SV18J.dat (экспериментальный) одинаковы и применяя любой из них можно проводить декодирова- ние данных спутников ГЛОНАСС. Прежде чем приступить к обработке данных файла FilGL1.txt рассмотрим процесс формирования и декодирования в навигационном приемнике информационного навигаци- онного сигнала. Этот процесс может быть представлен в виде следующей блок – схемы (рис. 5.7). Информация, передаваемая с каждого навигационного спутника представляет собой непрерывный поток информационных символов. Начала отсчетов передачи данных нуль часов Московского времени (начало передачи первого суперкадра каждые сутки). Структурно поток цифровой информации (ЦИ) формируется в виде непрерывно по- вторяющихся суперкадров (рис. 5.8). 206 207Рис. 5.7. Формирование и декодирование радионавигационного сигнала 1. Формирование навигационной информации 2. Форматирование навигационной информации 3. Преобразование навигационной информации в двоичный код 4. Формирование информационных сигналов, соот- ветствующих двоичному коду методом двухпози- ционной фазовой манипуляции с частотой 50 Гц 5. Преобразование информационных сигналов в сигналы с относительной фазовой манипуляцией 6. Сложение по модулю 2 информационных сигналов с меандровым колебанием частоты 100 Гц и выделе- ние 170 символов длительностью 10 миллисекунд 7. Дополнение 170 символов тридцатью символами кода метки времени (символ кода метки времени име- ет длитальность 10 миллисеконд). Окончание формирования одной строки навигацион- ногокадра Итог: 200 символов по 10 миллисекунд (2 секунды- время передачи 1 строки) А 208Рис. 5.7. Формирование и декодирование радионавигационного сигнала (продолжение) А 9. Модуляция несущих частот спутника полученным сигналом 10. Демодуляция радионавигационного сигнала 11.Корреляционная обработка демодулированного сигнала. Итог: поток символов длительностью 10 миллисекунд 12. Определение границ строк Выделяются строки длиной по 170 бит 13. Выделение информационных символов каждой строки Выделяются 77 символов информации и 8 провероч- ных символов (код Хемминга) 14. Проверка правильности считанной информации и исправление ошибок 15. Дешифровка двоичных символом 15. Декодирование двоичных символом 16. Перевод декодированных сигналов в десятичный эквивалент 17. Конец декодирования 8. Сложение по модулю 2 информационных сигналов с псевдослучайным кодом стандартной точности с частотой 0.511 МГц 209Рис. 5.8. Изображение потока данных навигационного спутника В приведенной далее программе декодирования навигационных данных спутников ГЛОНАСС выполняются процедуры п. п. 12… 17 (рис. 5.7) по данным п. 11 в соответст- вии с ИКД ГЛОНАСС. 5. Задание 2. Откройте файл Decod_GL3.m. Изучите программные процедуры файла, ис- пользуя его текст и расширенный комментарий. Выполните файл и в командном окне MatLab прочитайте результат декодирования и занесите его в отчет.5.5.3 Задания и вопросы для самоподготовки 1. Изучите форматы сообщений спутников ГЛОНАСС? 2. Что понимается под суперкадром, кадром, строкой и словом в сообщениях спутников ГЛОНАСС? 3. Что понимается под оперативной и неоперативной информацией, передаваемой со спутников ГЛОНАСС? 4. Как контролируется достоверность информации, передаваемой со спутников ГЛО-НАСС? 5. Как выполняется преобразование между системами счисления: переход из двоичной системы счисления в десятичную и наоборот; переход из шестнадцатеричной системы счисления в десятичную и наоборот; переход из двоичной системы счисления в шест- надцатеричную и наоборот? 6. Что такой бит, байт, младший разряд, старший разряд? …Метеки времени Информационные символы строк (строка передается в течение 1.7 секунды) Начало суток Окончание суток 5.5.4 Листинг файла:Decod_GL3.m %Имя m- файла:Decod_GL3.m %Программа декодирования данных спутников ГЛОНАСС %Символы метки времени Met=[1 1 1 1 1 0 0 0 1 1 0 1 1 1 0 1 0 1 0 0 0 0 1 0 0 1 0 1 1 0]; %Открытие файла для считывания данных %fid = fopen('raw_bits_SV18J.dat','rb') %fid = fopen( 'SV18fw_bin.txt','rb') fid = fopen('FilGL1_bin.txt','rb') if fid=-1 kol = 0;%Порядковый номер считываемых данных while (kol < 24) %Поиск и обнаружение меток времени в массиве данных for i = 1:30 MAS(i) = fscanf(fid,'%d',1); end C = xor(MAS, Met); D = any(C); kk = 0; while ( (feof(fid)==0) & (D == 1) ) MAS(1:29) = MAS(2:30); ch = fscanf(fid,'%d',1); MAS(30) = ch; C = xor(MAS, Met); D = any(C); %0 - all 0 kk = kk + 1; end % while %Считывание данных строк kol = kol + 1 for i = 1:170 danDDD(i)= (fscanf(fid,'%d',1)); end; str_danDDD=sprintf('%d',danDDD(1:170)); %Сложение данных по модулю 2 с меандром for i = 1:170 meandr = rem((i-1), 2); %остаток от деления на 2 ((i-1)/2) danDD(i) = xor(danDDD(i), meandr); end; % for i = 1:170 210 % DD = danDD %Выделение информационных символов dd1 = danDD(1:2:169); dd2 = danDD(2:2:170); ddxor = xor(dd1, dd2); Dany = any(ddxor); str1(1) = dd1(1); %Перевод данных из относительного кода for i = 2:85 str1(i) = xor(dd1(i-1), dd1(i)); end; % Контроля четности (правильности) принятых данных Skod(1,1:85)=[0 0 0 0 0 0 0 0 1 1 0 1 1 0 1 0 1 0 1 1 0 1 0 1 0 1 ... 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 ... 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 1 0 1 0 1 0 1 ... 0 1 0 1 0 1 0 1 0 1 0 1 0]; Skod(2,1:85)=[0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 0 0 1 1 0 1 1 0 0 1 1 ... 0 0 1 1 0 0 1 1 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 ... 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 1 1 0 0 1 1 ... 0 0 1 1 0 0 1 1 0 0 1 1 0]; Skod(3,1:85)=[0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 1 1 1 0 0 0 1 1 1 1 ... 0 0 0 0 1 1 1 1 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 ... 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 1 1 1 1 ... 0 0 0 0 1 1 1 1 0 0 0 0 1]; Skod(4,1:85)=[0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 ... 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 ... 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 ... 1 1 1 1 1 1 1 1 0 0 0 0 0]; Skod(5,1:85)=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 1 1 1 1 0]; Skod(6,1:85)=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0 0 0]; 211 Skod(7,1:85)=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1 1 1]; Skod(8,1:85)=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ... 1 1 1 1 1 1 1 1 1 1 1 1 1]; E = rot90(eye(85)); for k=1:8 Vi1 = Skod(k,1:85) * E; STR1Vi = and(str1,Vi1); a1 = sum(STR1Vi); n = floor(a1 / 2); c1 = a1 - n * 2; if (k <8) cc_all(kol,k)=xor(c1,str1(86 - k)); else cc_all(kol,k)=c1; end; end; %Декодирование номеров строк num4 = str1(2:5); NN1 = sprintf('%d',num4); num(kol)=bin2dec(NN1); %Вывод данных строк str_danDDD1=sprintf('%d',num(kol),danDDD(1:85)); str_danDDD2=sprintf('%d',num(kol),danDDD(86:170)); s_danDDD=[str_danDDD1;str_danDDD2]; %++++++++++++++++++++++++++++++++++++++++++ %Вывод номера считываемых данных (kol), номера строки %[num(kol)] и коэффициентов четности cc_all(kol,1:8)= %C8 C7 C6 C5 C4 C3 C2 C1 out1 = [kol num(kol) cc_all(kol,8:-1:1)] %Вывод считываемых данных строки SS1=sprintf('%d',str1); %Декодирование данных строк согласно их номеров 212 % в кадре [например, num(kol) == 1- cтрока номер 1] if num(kol) == 1 zx=1; zxdot=1; zx2dot=1; Nstr=sprintf('%d',str1(2:5)); R1= sprintf('%d',str1(6:7)); P1=sprintf('%d',str1(8:9)); tk5=sprintf('%d',str1(10:14)); tk6=sprintf('%d',str1(15:20)); tk1=sprintf('%d',str1(21)); if sprintf('%d',str1(22))=='1' zxdot = -1; end; Xdot=sprintf('%d',str1(23:45)); if sprintf('%d',str1(46))=='1' zx2dot = -1; end; X2dot=sprintf('%d',str1(47:50)); if sprintf('%d',str1(51))=='1' zx = -1; end; X=sprintf('%d',str1(52:77)); St1=bin2dec({Nstr tk5 tk6 tk1 X Xdot X2dot}); Mng=[1 1 1 30 zx*2^(-11) zxdot*2^(-20) zx2dot*2^(-30)]; Str1=(St1'.*Mng)'%вывод декодированных данных end; if num(kol) == 2 zy=1; zydot=1; zy2dot=1; Nstr=sprintf('%d',str1(2:5)); Bn= sprintf('%d',str1(6:8)); P2=sprintf('%d',str1(9)); tb=sprintf('%d',str1(10:16)); R2=sprintf('%d',str1(17:21)); if sprintf('%d',str1(22))=='1' zydot = -1; end; Ydot=sprintf('%d',str1(23:45)); if sprintf('%d',str1(46))=='1' zy2dot = -1; end; 213 Y2dot=sprintf('%d',str1(47:50)); if sprintf('%d',str1(51))=='1' zy = -1; end; Y=sprintf('%d',str1(52:77)); St2=bin2dec({Nstr Bn P2 tb Y Ydot Y2dot}); Mng=[1 1 1 15*60 zy*2^(-11) zydot*2^(-20) zy2dot*2^(-30)]; Str2=(St2'.*Mng)'%вывод декодированных данных end; if num(kol) == 3 zz=1; zzdot=1; zz2dot=1; zgamma=1; Nstr=sprintf('%d',str1(2:5)); P3= sprintf('%d',str1(6)); if sprintf('%d',str1(7))=='1' zgamma = -1; end; Gamma=sprintf('%d',str1(8:17)); R3=sprintf('%d',str1(18:19)); P=sprintf('%d',str1(20)); ln=sprintf('%d',str1(21)); if sprintf('%d',str1(22))=='1' zzdot = -1; end; Zdot=sprintf('%d',str1(23:45)); if sprintf('%d',str1(46))=='1' zz2dot = -1; end; Z2dot=sprintf('%d',str1(47:50)); if sprintf('%d',str1(51))=='1' zz = -1; end; Z=sprintf('%d',str1(52:77)); St3=bin2dec({Nstr P3 ln Gamma Z Zdot Z2dot}); Mng=[1 1 1 zgamma*2^(-40) zz*2^(-11) zzdot*2^(-20) zz2dot*2^(-30)]; Str3=(St3'.*Mng)'%вывод декодированных данных end; 214 if num(kol) == 4 ztaun=1; zdtaun=1; Nstr=sprintf('%d',str1(2:5)); if sprintf('%d',str1(6))=='1' ztaun = -1; end; taun= sprintf('%d',str1(7:27)); if sprintf('%d',str1(28))=='1' zdtaun = -1; end; deltatau=sprintf('%d',str1(29:32)); En=sprintf('%d',str1(33:37)); R4=sprintf('%d',str1(38:51)); P4=sprintf('%d',str1(52)); FT=sprintf('%d',str1(53:56)); R4i4=sprintf('%d',str1(57:59)); NT=sprintf('%d',str1(60:70)); nSV=sprintf('%d',str1(71:75)); modifSV=sprintf('%d',str1(76:77)); St4=bin2dec({Nstr taun deltatau En FT NT nSV modifSV}); Mng=[1 ztaun*2^(-30) zdtaun*2^(-30) 1 1 1 1 1]; Str4=(St4'.*Mng)'%вывод декодированных данных end; if num(kol) == 5 ztauc=1; ztauGPS=1; Nstr=sprintf('%d',str1(2:5)); NA= sprintf('%d',str1(6:16)); if sprintf('%d',str1(17))=='1' ztauc = -1; end; tauc=sprintf('%d',str1(18:48)); R5=sprintf('%d',str1(49)); N4=sprintf('%d',str1(50:54)); if sprintf('%d',str1(55))=='1' ztauGPS = -1; end; tauGPS=sprintf('%d',str1(56:76)); ln=sprintf('%d',str1(77)); St5=bin2dec({Nstr NA tauc N4 tauGPS ln Nstr Nstr}); Mng=[1 1 ztauc*2^(-27) 1 ztauGPS*2^(-30) 1 0 0]; 215 Str5=(St5'.*Mng)'%вывод декодированных данных end; if num(kol) == 6 ztaun=1; zlambda=1;zdeltai=1; Nstr=sprintf('%d',str1(2:5)); CNA= sprintf('%d',str1(6)); Mn=sprintf('%d',str1(7:8)); nSV=sprintf('%d',str1(9:13)); if sprintf('%d',str1(14))=='1' ztaun = -1; end; taun=sprintf('%d',str1(15:23)); if sprintf('%d',str1(24))=='1' zlambda = -1; end; lambda=sprintf('%d',str1(25:44)); if sprintf('%d',str1(45))=='1' zdeltai= -1; end; deltai=sprintf('%d',str1(46:62)); ecs=sprintf('%d',str1(63:77)); St6=bin2dec({Nstr Mn nSV taun lambda deltai ecs}); Mng=[1 1 1 ztaun*2^(-18) zlambda*2^(-20) zdeltai*2^(-20) 2^(-20)]; Str6=(St6'.*Mng)'%вывод декодированных данных end; if num(kol) == 7 zomegan=1;zdeltaT=1; zdeltaTdot=1; Nstr=sprintf('%d',str1(2:5)); if sprintf('%d',str1(6))=='1' zomegan = -1; end; omegan= sprintf('%d',str1(7:21)); tlambda=sprintf('%d',str1(22:44)); if sprintf('%d',str1(43))=='1' zdeltaT = -1; end; deltaT=sprintf('%d',str1(44:64)); if sprintf('%d',str1(65))=='1' zdeltaTdot= -1; end; 216 deltaTdot=sprintf('%d',str1(66:71)); Hn=sprintf('%d',str1(72:76)); ln=sprintf('%d',str1(77)); St7=bin2dec({Nstr omegan tlambda deltaT deltaTdot Hn ln}); Mng=[1 zomegan*2^(-15) 2^(-5) zdeltaT*2^(-9) zdeltaTdot*2^(-14) 1 1]; Str7=(St7'.*Mng)'%вывод декодированных данных end; if num(kol) == 8 ztaun=1; zlambda=1;zdeltai=1; Nstr=sprintf('%d',str1(2:5)); CNA= sprintf('%d',str1(6)); Mn=sprintf('%d',str1(7:8)); nSV=sprintf('%d',str1(9:13)); if sprintf('%d',str1(14))=='1' ztaun = -1; end; taun=sprintf('%d',str1(15:23)); if sprintf('%d',str1(24))=='1' zlambda = -1; end; lambda=sprintf('%d',str1(25:44)); if sprintf('%d',str1(45))=='1' zdeltai= -1; end; deltai=sprintf('%d',str1(46:62)); ecs=sprintf('%d',str1(63:77)); St8=bin2dec({Nstr Mn nSV taun lambda deltai ecs}); Mng=[1 1 1 ztaun*2^(-18) zlambda*2^(-20) zdeltai*2^(-20) 2^(-20)]; Str8=(St8'.*Mng)' %вывод декодированных данных end; if num(kol) == 9 zomegan=1;zdeltaT=1; zdeltaTdot=1; Nstr=sprintf('%d',str1(2:5)); if sprintf('%d',str1(6))=='1' zomegan = -1; end; omegan= sprintf('%d',str1(7:21)); tlambda=sprintf('%d',str1(22:44)); if sprintf('%d',str1(43))=='1' zdeltaT = -1; 217 end; deltaT=sprintf('%d',str1(44:64)); if sprintf('%d',str1(65))=='1' zdeltaTdot= -1; end; deltaTdot=sprintf('%d',str1(66:71)); Hn=sprintf('%d',str1(72:76)); ln=sprintf('%d',str1(77)); St9=bin2dec({Nstr omegan tlambda deltaT deltaTdot Hn ln}); Mng=[1 zomegan*2^(-15) 2^(-5) zdeltaT*2^(-9) zdeltaTdot*2^(-14) 1 1]; Str9=(St9'.*Mng)'%вывод декодированных данных end; if num(kol) == 10 ztaun=1; zlambda=1;zdeltai=1; Nstr=sprintf('%d',str1(2:5)); CNA= sprintf('%d',str1(6)); Mn=sprintf('%d',str1(7:8)); nSV=sprintf('%d',str1(9:13)); if sprintf('%d',str1(14))=='1' ztaun = -1; end; taun=sprintf('%d',str1(15:23)); if sprintf('%d',str1(24))=='1' zlambda = -1; end; lambda=sprintf('%d',str1(25:44)); if sprintf('%d',str1(45))=='1' zdeltai= -1; end; deltai=sprintf('%d',str1(46:62)); ecs=sprintf('%d',str1(63:77)); St10=bin2dec({Nstr Mn nSV taun lambda deltai ecs}); Mng=[1 1 1 ztaun*2^(-18) zlambda*2^(-20) zdeltai*2^(-20) 2^(-20)]; Str10=(St10'.*Mng)'%вывод декодированных данных end; if num(kol) == 11 zomegan=1;zdeltaT=1; zdeltaTdot=1; Nstr=sprintf('%d',str1(2:5)); if sprintf('%d',str1(6))=='1' zomegan = -1; 218 end; omegan= sprintf('%d',str1(7:21)); tlambda=sprintf('%d',str1(22:44)); if sprintf('%d',str1(43))=='1' zdeltaT = -1; end; deltaT=sprintf('%d',str1(44:64)); if sprintf('%d',str1(65))=='1' zdeltaTdot= -1; end; deltaTdot=sprintf('%d',str1(66:71)); Hn=sprintf('%d',str1(72:76)); ln=sprintf('%d',str1(77)); St11=bin2dec({Nstr omegan tlambda deltaT deltaTdot Hn ln}); Mng=[1 zomegan*2^(-15) 2^(-5) zdeltaT*2^(-9) zdeltaTdot*2^(-14) 1 1]; Str11=(St11'.*Mng)'%вывод декодированных данных end; if num(kol) == 12 ztaun=1; zlambda=1;zdeltai=1; Nstr=sprintf('%d',str1(2:5)); CNA= sprintf('%d',str1(6)); Mn=sprintf('%d',str1(7:8)); nSV=sprintf('%d',str1(9:13)); if sprintf('%d',str1(14))=='1' ztaun = -1; end; taun=sprintf('%d',str1(15:23)); if sprintf('%d',str1(24))=='1' zlambda = -1; end; lambda=sprintf('%d',str1(25:44)); if sprintf('%d',str1(45))=='1' zdeltai= -1; end; deltai=sprintf('%d',str1(46:62)); ecs=sprintf('%d',str1(63:77)); St12=bin2dec({Nstr Mn nSV taun lambda deltai ecs}); Mng=[1 1 1 ztaun*2^(-18) zlambda*2^(-20) zdeltai*2^(-20) 2^(-20)]; Str12=(St12'.*Mng)'%вывод декодированных данных end; 219 if num(kol) == 13 zomegan=1;zdeltaT=1; zdeltaTdot=1; Nstr=sprintf('%d',str1(2:5)); if sprintf('%d',str1(6))=='1' zomegan = -1; end; omegan= sprintf('%d',str1(7:21)); tlambda=sprintf('%d',str1(22:44)); if sprintf('%d',str1(43))=='1' zdeltaT = -1; end; deltaT=sprintf('%d',str1(44:64)); if sprintf('%d',str1(65))=='1' zdeltaTdot= -1; end; deltaTdot=sprintf('%d',str1(66:71)); Hn=sprintf('%d',str1(72:76)); ln=sprintf('%d',str1(77)); St13=bin2dec({Nstr omegan tlambda deltaT deltaTdot Hn ln}); Mng=[1 zomegan*2^(-15) 2^(-5) zdeltaT*2^(-9) zdeltaTdot*2^(-14) 1 1]; Str13=(St13'.*Mng)'%вывод декодированных данных end; if num(kol) == 14 ztaun=1; zlambda=1;zdeltai=1; Nstr=sprintf('%d',str1(2:5)); CNA= sprintf('%d',str1(6)); Mn=sprintf('%d',str1(7:8)); nSV=sprintf('%d',str1(9:13)); if sprintf('%d',str1(14))=='1' ztaun = -1; end; taun=sprintf('%d',str1(15:23)); if sprintf('%d',str1(24))=='1' zlambda = -1; end; lambda=sprintf('%d',str1(25:44)); if sprintf('%d',str1(45))=='1' zdeltai= -1; end; deltai=sprintf('%d',str1(46:62)); ecs=sprintf('%d',str1(63:77)); 220 St14=bin2dec({Nstr Mn nSV taun lambda deltai ecs}); Mng=[1 1 1 ztaun*2^(-18) zlambda*2^(-20) zdeltai*2^(-20) 2^(-20)]; Str14=(St14'.*Mng)'%вывод декодированных данных end; %Строка 14 для 5 кадра if num(kol) == 14 zB1=1; zB2=1; Nstr=sprintf('%d',str1(2:5)); if sprintf('%d',str1(6))=='1' zB1 = -1; end; B1=sprintf('%d',str1(7:16)); if sprintf('%d',str1(17))=='1' zB2 = -1; end; B2=sprintf('%d',str1(18:26)); Kp=sprintf('%d',str1(27:28)); St14_5=bin2dec({Nstr B1 B2 Kp}); Mng=[1 zB1*2^(-10) zB2*2^(-16) 1]; Str14_5=(St14_5'.*Mng)'%вывод декодированных данных end; if num(kol) == 15 zomegan=1;zdeltaT=1; zdeltaTdot=1; Nstr=sprintf('%d',str1(2:5)); if sprintf('%d',str1(6))=='1' zomegan = -1; end; omegan= sprintf('%d',str1(7:21)); tlambda=sprintf('%d',str1(22:44)); if sprintf('%d',str1(43))=='1' zdeltaT = -1; end; deltaT=sprintf('%d',str1(44:64)); if sprintf('%d',str1(65))=='1' zdeltaTdot= -1; end; deltaTdot=sprintf('%d',str1(66:71)); Hn=sprintf('%d',str1(72:76)); ln=sprintf('%d',str1(77)); 221 St15=bin2dec({Nstr omegan tlambda deltaT deltaTdot Hn ln}); Mng=[1 zomegan*2^(-15) 2^(-5) zdeltaT*2^(-9) zdeltaTdot*2^(-14) 1 1]; Str15=(St15'.*Mng)'%вывод декодированных данных end; str1; end; %while (feof(fid)==0) | % while (kol < 14) end % if fid=-1 fclose(fid); %Технологические данные Qout=[num' cc_all] %Номер строки и проверка на четность %Qoutnum1=Qout(1,2:9); %SS1=sprintf('%d',str1); %Qout123=[Str1 Str2 Str3] %Qout45=[Str4 Str5] %Qout6789=[Str6 Str7 Str8 Str9] %Qout10_11_12=[Str10 Str11 Str12] %Qout13_14_15=[Str13 Str14 Str15] clear; 5.5.5 Пример выполнения файла:Decod_GL3.m В результате выполнения программы в командном окне Command Window MatLab отображаются декодированные данные по каждой из строк в последовательности и соответствии ИКД ГЛОНАСС. Так, например, структура строки 1согласно ИКД ГЛОНАСС имеет вид 1 2 - 5 6 - 7 8 - 9 10 - 21 22 - 45 85 84 -81 80 - 79 78 - 77 76 - 65 64 - 410 m P1 t k)t(nXb.Номер строки в навигацион- ном кадре Резерв Признак смены оператив- ной инфо- рмации Время начала кадра внутри текущих суток, исчисляемое в шкале бортового времени. Начало суток по бортовому времени спут- ника совпадает с началом очередного су- перкадра. Составляющие вектора скорости n-го спутника в системе коорди- нат ПЗ-90 на момент времени tb46 - 50 51 - 77 1 - 8 1 - 30 40 - 36 35 - 9 8 - 1 1 - 30 )t(nXb..)t(nXbСоставляющие ускорения n-го спутника на момент времени tb, Координаты n-го спутника в системе координат ПЗ-90 на момент времени t b; Код Хемминга Метка времени Рис. 5.9. - Структура 1 строки для (1- 4) кадров суперкадра 222 Декодированная информация из командного окна MatLab приведена на рис. 5.10. Рис. 5.10 Декодированные данные спутника ГЛОНАСС Данные, изображенные на рис. 5.10 соответствуют рис. 5.9, но записаны в столбец. 223 32>1   ...   6   7   8   9   10   11   12   13   14

id = alm(nom).ID; while ( i < 24) id = alm(nom).ID;
Health = alm(nom).Health;
%fprintf('1: i=%i k=%i nom=%i id=%i Health=%i \n', i, k, nom, id, Health); if ( id > 0)
Health = alm(nom).Health; if ( Health == 0) k = k + 1; nom_ns(k) = id;
% fprintf('2: i=%i k=%i nom=%i id=%i Health=%i \n', i, k, nom, id, Health); nom = nom + 1; else nom = nom + 1; end; else nom = nom + 1; end; i = i + 1; end; % i kol = k;
% fprintf('kol=%i \n', kol); % для вывода номеров здоровых спутников перед fprintf убрать "%" nom_ns % - номера здоровых навигационных спутников,
%количество спутников, для которых строятся графики орбитального движения kol = 1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nom_ns(1:kol)=[24];%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%nom_ns(1:kol)=[1 2 3 4 7 8 19 21 22 23 24 ];%номера спутников (количество номеров должно
%равняться количеству спутников
KOL_GLN = 24;
DEG_TO_RAD = 0.017453292519943; % (PI / 180.00 )
RAD_TO_DEG = 57.295779513082; % (180.00 / PI )
A_WGS84_M = 6378137.0 ; % WGS-84 ellipsoid parameters
B_WGS84_M = 6356752.314; % WGS-84 ellipsoid parameters
A_PZ90_M = 6378136.0; %6 378 136 м - Equatorial radius of the Earth - ,большая полуось эллипсоида
B_PZ90_M = 6356751.36174571344; %AP_LAND (m) Polar radius of the Earth
%A_PZ90_KM = 6378.136; %(Km) Equatorial radius of the Earth
B_PZ90_KM = 6356.75136174571344; % AP_LAND (Km) Polar radius of the Earth
FACTOR_PZ90 = 1.0/298.257839303;% Коэффициент сжатия эллипсоида
C_LIGHT_M = 2.99792458E8; % m/sec Speed of light
RAD_IN_SEC = 7.2722052166430e-5 ;%(PI/43200.0) Number radian in second of time
2*PI/(24*3600)=PI/43200;
104
dt_lsf = 14;
% Вводим координаты приемника:
% current_loc_wgs - Pos. reciver (WGS-84): current_loc_wgs.lat = DEG_TO_RAD * (50. + 29.0 / 60.0 + 36.78 /
3600.0);%%%%%%%%%%%%%%%% current_loc_wgs.lon = DEG_TO_RAD * (30. + 27.0 / 60.0 + 50.5 /
3600.0);%%%%%%%%%%%%%%%% current_loc_wgs.h = 120.9; % метр, WGS 84 %%%%%%%%%%%%%%%%
%current_loc_wgs.hae = 0.1229; % Km PZ90
%Преобразование координат rec_pos_xyz_wgs = LLH_to_ECEF( A_WGS84_M, B_WGS84_M, current_loc_wgs); rec_pos_xyz_pz90_m = WGS84_to_PZ90(rec_pos_xyz_wgs); current_loc_pz90 = ECEF_to_LLH(A_PZ90_M, B_PZ90_M, rec_pos_xyz_pz90_m); step_t = 600;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%L = 12*24;
L =24*1 * 3600;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ti_start = 0; ti_start = Data_observ.ti + day_observ_alm * 86400; ti_end = ti_start + L; j = 0; for ( ti = ti_start : step_t : ti_end ) j = j + 1;
%ПЕРЕСЧЕТ ВРЕМЕНИ из секунд в часы, минуты и секунды: t(j) = ti;
A(j) = mod(t(j), 86400); hour(j) = floor(A(j) / 3600); m(j) = floor((A(j) - hour(j)*3600) / 60); sek(j) = A(j) - hour(j) * 3600 - m(j) * 60; for (n=1:24) % обнуление массивов в текущий момент времени для всех спутников x(j, n) = 0; y(j, n) = 0; z(j, n)= 0; ris_vid(n, j) = 0; end;
%timeUTC.ti = mod( ti, 86400); % текущее время обсервации от начала суток timeUTC.ti = ti;
105

[kol_gln_a, satvis_gln, satpos_gln_a, satpos_gln_ecef] = GLN_satfind( A_PZ90_M, B_PZ90_M, timeUTC, current_loc_pz90, alm); for (i=1: KOL_GLN) if ( satvis_gln(i).el >= 0) el = RAD_TO_DEG * satvis_gln(i).el; az = RAD_TO_DEG * satvis_gln(i).az;
AZ(j,i) = az;%уголы азимута спутников в градусах
EL(j,i) = el;%уголы видимости спутников в градусах ris_vid(i, j) = i;
% fprintf('i=%2i el=%6.2f az=%7.2f x=%12.2f y=%12.2f z=%12.2f \n', i, el, az, satpos_gln_a(i).x, sat- pos_gln_a(i).y, satpos_gln_a(i).z); else
AZ(j,i) = 0.0;
EL(j,i) = 0.0; end; end; for (i=1:KOL_GLN) prn = alm(i).ID; health = alm(i).Health; if ( (prn > 0) & (health == 0))
%Координаты спутников в системе ECEF. Для визуализации добавить в вывод графиков и
% снять блокировку. x(j,i)= satpos_gln_ecef(i).x; y(j,i)= satpos_gln_ecef(i).y; z(j,i)= satpos_gln_ecef(i).z;
%Координаты спутников в системе ECI. Для визуализации добавить в вывод графиков. x1(j,i)= satpos_gln_a(i).x; y1(j,i)= satpos_gln_a(i).y; z1(j,i)= satpos_gln_a(i).z;
%Скорости спутников. Для визуализации добавить в вывод графиков. x2(j,i)= satpos_gln_ecef(i).vx; y2(j,i)= satpos_gln_ecef(i).vy; z2(j,i)= satpos_gln_ecef(i).vz;
%Углы видимости и дальности до спутников. Для визуализации добавить в вывод графиков. x3(j,i)= satvis_gln(i).el*180/pi; y3(j,i)= satvis_gln(i).az; z3(j,i)= satvis_gln(i).r; end; end; end; % ti ( j )
106
max_n = 24;
%для блокировки вывода графика "Время наблюдения спутников ГЛОНАСС" перед функцией
%ris_vis_sat в следующей строке поставить % ris_vis_sat(max_n, j, ti_start, step_t, ris_vid); for (i=1:kol) j_color = j_color + 1; if (j_color > 14 ) j_color = 1; end
S = color6(j_color); prn = nom_ns(i); hold on h_F1 = gca;
%{
%ГРАФИКИ ПАРАМЕТРОВ ОРБИТАЛЬНОГО ДВИЖЕНИЯ СПУТНИКОВ ГЛОНАСС
%График орбит спутников ГЛОНАСС системе ECEF plot3(x(:,prn),y(:,prn),z(:,prn),S,'LineWidth',0.5); axis([ -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7)]); set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); set(h_F1,'Position',[0.1 0.1 0.85 0.9]) ; xlabel('Координата X') ylabel('Координата Y'), zlabel ('Координата Z'),grid on str1 = num2str( prn); text(x(j,prn), y(j,prn),z(j,prn),str1,'FontSize',14,'FontName','TimesNewRoman','HorizontalAlignment','center' );
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
%График орбит спутников ГЛОНАСС системе ECI plot3(x1(:,prn),y1(:,prn),z1(:,prn),S, 'LineWidth',1); axis([ -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7)]) set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); set(h_F1,'Position',[0.1 0.1 0.85 0.9]) ; xlabel('Координата X') ylabel('Координата Y'), zlabel('Координата Z'),grid on str1 = num2str( prn);
107
text(x1(1,prn), y1(1,prn),z1(1,prn),str1,'FontSize',14,'FontName','TimesNewRoman','VerticalAlignment',
'cap' );
%} end
1   2   3   4   5   6   7   8   9   10   ...   14

Функция Read_GL_Alm
function [alm,max_kol] = Read_GL_Alm(Dat);
%Имя функции: Read_GL_Alm
%Функция читает данные альманаха спутников ГЛОНАСС.
%Входные данные:
%Dat - имя файла, содержащего альманах спутников ГЛОНАСС , например, Dat = 'GLN_all.alm';
%Выходные данные:
%alm - данные альманаха спутников ГЛОНАСС (структура) ,
%max_kol- максимальное количество спутников в альманахе fid =fopen(Dat,'rt');%открытие файла max_kol = 0;
%чтение данных из файла while not(feof(fid)) s1=fscanf(fid,'%s',6); if not(feof(fid)) lenstr = length(s1); max_kol = max_kol + 1;
% while (fscanf(fid,'%s',1) == '-') end str1 = fscanf(fid,'%s',1);
% lenstr = length(str1); str2 = fscanf(fid,'%s',1); lenstr = length(str2); n_sv = sscanf(str2,'%d'); strID = str2(1:lenstr);
ID = sscanf(strID,'%d'); alm(ID).ID = ID; %номер спутника (1...24)
%%%%%%%%%%% t_1=fscanf(fid,'%s',1); alm(ID).Hn=fscanf(fid,'%d'); % номер частоты спутника t_2=fscanf(fid,'%s',1); alm(ID).Health=fscanf(fid,'%d',1);%здоровье (0- спутник здоров) t_3=fscanf(fid,'%s',1); alm(ID).ecc = fscanf(fid,'%g',1);% эксцентриситет орбиты спутника t_4=fscanf(fid,'%s',1);
108
alm(ID).Na =fscanf(fid,'%g',1 );% номер суток к которым относятся данные альманаха, отсчиты- ваемых
%от ближайшего високосного года t_5=fscanf(fid,'%s',3); alm(ID).deltai = fscanf(fid,'%g',1); % наклонение орбиты спутника (радианы) t_6=fscanf(fid,'%s',2); alm(ID).LambdaN = fscanf(fid,'%g',1); % долгота восходящего узла орбиты спутника (радианы) while not(fscanf(fid,'%c',1) == ':') end
% t_7=fscanf(fid,'%s',1) alm(ID).TLambdaN = fscanf(fid,'%g',1) - 10800;% время прохождения восходящего узла орбиты спутника,
%к которому относятся данные альманаха, приведенное к времени UTC (секунда) t_8=fscanf(fid,'%s',2); alm(ID).omegan = fscanf(fid,'%g',1); %аргумент перигея (радианы) t_11=fscanf(fid,'%s',1); alm(ID).Tdr=fscanf(fid,'%g',1); % драконический период (секунда/виток) t_12=fscanf(fid,'%s',1); alm(ID).dTdr=fscanf(fid,'%g',1); % скорость изменения драконического периода (секунда/виток в квадрате) t_13=fscanf(fid,'%s',2); alm(ID).tau_n=fscanf(fid,'%g',1) / 1000.0;%коэффициет коррекции шкалы времени (сдвиг времени спутника
%относительно системного времени ГЛОНАСС в секундах) end; end fclose(fid);
Функция Gln_data_from_NA
function [time_UTC] = Gln_data_from_NA(leap_year, day_from_leap);
%Имя функции:Gln_data_from_NA
%Функция предназначена для преобразования номера дня NA (день привязки альманаха
% от ближайшего високосного года) в текущую дату
%function [time_UTC] = Gln_data_from_NA(leap_year, day_from_leap);
%Входные данные:
% leap_year - ближайший високосный год
%day_from_leap (NA) - день привязки альманаха
%Выходные данные:
%структура timeUTC (год, месяц, день) - текущая дата
DnMon = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]; %количество дней в месяцах
109
n4 = mod(leap_year, 4); n100 = mod(leap_year, 100); n400 = mod(leap_year, 400); if (n4 == 0) n_leap = 1; else n_leap = 0; end if ((n100 == 0) & (n400 > 0)) n_leap = 0; end if (day_from_leap > (365 + n_leap)) day = day_from_leap - (365 + n_leap); k =fix (day / 365); day = mod(day, 365); else day = day_from_leap; k = -1; end; god = leap_year + k + 1; if (god > leap_year) n_leap = 0; end; mon = 1; mon_day = 31; while (day > mon_day) day = day - mon_day; mon = mon + 1; % LK mon_day = DnMon(mon);
% mon = mon+1; if (mon == 2) mon_day = mon_day + n_leap; end end time_UTC.year = god; time_UTC.mon = mon; time_UTC.day = day;
Функция JD_data
function [JD, day_year] = JD_data(timeUTC)
%Имя:JD_data
% Функция JD_data(timeUTC) вычисляет :
110

%JD - номер юлианского дня, day_year - номер дня года.
%Входные данные:
%Структура timeUTC
%timeUTC.year - год,
% timeUTC.mon - месяц,
% timeUTC.day - день.
%Выходные данные:
%JD - юлианский день;
%day_year- день от начала года.
%количество дней в месяцах
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;
Функция LLH_to_ECEF
function [XYZ] = LLH_to_ECEF( a, b, llh)
%Имя функции:LLH_to_ECEF.m
%Функция преобразования географических координат в прямоугольную геоцентрическую систему координат (ECEF)
%Входные данные:
%Структура llh
%llh.lon-долгота (радиан),
%llh.lat-широта (радиан),
%llh.h-высота (метр);
%a, b-большая и малая полуоси земного эллипсоида в WGS-84 (метр);
%Выходные данные:
111

%Структура XYZ
%XYZ.x - координата x в ECEF;
%XYZ.y - координата y в ECEF;
%XYZ.z - координата z в ECEF; a2=a * a; b2=b * b; r = a2 / sqrt(a2 * cos(llh.lat) * cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat));
XYZ.x = (r + llh.h) * cos(llh.lat) * cos(llh.lon);
XYZ.y = (r + llh.h) * cos(llh.lat) * sin(llh.lon);
XYZ.z = (b2 / a2 * r + llh.h) * sin(llh.lat);
Функция WGS84_to_PZ90
function [pos_pz90] = WGS84_to_PZ90( pos_wgs84)
%Имя функции: WGS84_to_PZ90.m
%Функция преобразует координаты из системы WGS 84 в систему ПЗ 90 p = (1.0 - 0.12e-6); x = pos_wgs84.x; y = pos_wgs84.y; z = pos_wgs84.z; pos_pz90.x = (p * x + p * 0.82e-6 * y + 1.4 ); pos_pz90.y = (-p * 0.82e-6 * x + p * y + 1.4); pos_pz90.z = (p * z + 0.9);
Функция ECEF_to_LLH
function [llh] = ECEF_to_LLH(a, b, XYZ)
%Имя функции:ECEF_to_LLH
%Функция преобразует координаты из системы ECEF в географическую систему
%Входные данные:
%a, b-большая и малая полуоси земного эллипсоида (метр); ( a=6378137.0; %b=6356752.314; метры)
%Структура XYZ
%XYZ.x - координата x в ECEF;
%XYZ.y - координата y в ECEF;
%XYZ.z- координата z в ECEF;
%Выходные данные:
%Структура llh
%llh.lon - долгота (радиан);
%llh.lat - широта (радиан);
%llh.h-высота (метр); a2 = a * a; b2 = b * b; xy = sqrt(XYZ.x * XYZ.x + XYZ.y * XYZ.y); thet = atan2(XYZ.z * a , (xy * b));
112

%thet = atan(XYZ.z * a / (xy * b)); esq = 1.0 - b2 / a2; epsq = a2 / b2 - 1.0; llh.lat = atan((XYZ.z + epsq * b * (sin(thet) ^3)) / (xy - esq * a * (cos(thet) ^3))); llh.lon = atan2(XYZ.y, XYZ.x);%! while (llh.lon < 0) llh.lon = 2*pi + llh.lon; end;
%if llh.lon < 0
% llh.lon = 2*pi + llh.lon;
%end ; r = a2 / sqrt(a2 * cos(llh.lat) * cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat)); llh.h = xy / cos(llh.lat) - r; end
Функция GLN_satfind
function [kol_gln_a, satvis_gln, satpos_gln_ns, satpos_gln_ecef] = GLN_satfind(a, b, timeUTC, cur- rent_loc_pz90, alm_gln);
%Имя функции:GLN_satfind
%Функция вычисляет координаты спутников ГЛОНАСС, углы видимости и количество видимых спутников
%Входные данные:
%a, b- большая и малая полуоси земного эллипсоида;
%timeUTC- параметры времени;
%current_loc_pz90- широта, долгота и высота точки, из которой наблюдаются спутники (структура);
% alm_gln-альманах спутников ГЛОНАСС (структура)
%Выходные данные:
%kol_gln_a- количество видимых спутников ГЛОНАСС,
%satvis_gln- углы видимости спутников (структура),
%satpos_gln_ns-координаты спутников (структура) в абсолютной системе координат (ECI)
%satpos_gln_ecef- координаты спутников (структура) в подвижной системе координат (ECEF)
RAD_TO_DEG = 180.00 / pi ;
KOL_GLN = 24; nom_ns = 1;
Min_Elev = 0; % минимальный угол видимости kol_gln_a = 0;
[nom_sat_gln_a, trac_sat_gln_a, sat_pos, satvis_gln] = init_data(KOL_GLN);
% число дней от фундаментальной эпохи 2000г от 12h, 0, January:
% jd_up_epoh = JD_from_epohi( 2000, timeUTC) - 0.5; % from 12h UTC 2000 nut = 0;
113

S0 = s0_Nut( timeUTC, nut); time_s0 = S0.s0_m_mod ; %time_s0 - истинное звездное время в день обсервации, 0h UTC year = timeUTC.year; leap_year = fix( year / 4) * 4; % ближайший к текущему (предыдущий) високосный год ti = timeUTC.ti; % текущее время обсервации от начала дня n00 = fix(ti / 86400);
% n0 - номер текущих суток внутри 4-х летнего периода (от ближайшего високосного года) n0 = JD_from_epohi(leap_year, timeUTC) + n00 + 1;
[eci_current_loc, eci_rec_pos_xyz] = llh_to_eci(a, b, ti, time_s0, current_loc_pz90); satpos_eci = init_satpos_gln(); for ( i = 1 : KOL_GLN) prn = alm_gln (i).ID; health = alm_gln (i).Health; if ( (prn > 0) & (health == 0)) satpos_eci = gln_a_1(i, n0, ti, time_s0, alm_gln);
[satpos_eci, satpos_gln_a] = satpos_eci_in_metr(satpos_eci); %Transformation coordinates from
"Km" in "m" satpos_gln_ns(prn) = satpos_gln_a;
%-----
[satpos_ecef] =eci_to_ecef(time_s0, ti, satpos_eci); satpos_gln_ecef(prn) = satpos_ecef;
% calculates angles of visibility by almanac:
[top] = top_coord(prn, eci_current_loc, eci_rec_pos_xyz,satpos_gln_a); satvis_gln(prn) = top; if (satvis_gln(prn).r > 0.0) el = RAD_TO_DEG * satvis_gln(prn).el; if ( el >= Min_Elev)
% fprintf('prn=%2i el=%f az=%f \n ', prn, el, RAD_TO_DEG * satvis_gln(prn).az); satpos_gln( prn) = satpos_eci;
[satpos_ecef] =eci_to_ecef(time_s0, ti, satpos_eci);
[trac_sat_gln_a] = rewrite_satpos(nom_ns, satpos_ecef); trac_sat_gln_a( nom_ns).r = satvis_gln( prn).r; trac_sat_gln_a (nom_ns).prn = prn; nom_sat_gln_a(nom_ns) = prn; nom_ns = nom_ns + 1; % number of visible satellites end; % if ( el >= Min_Elev) end; % if ( satvis_gln[prn-1].r > 0.0)
114
end; % f ( (prn > 0) & (health == 1)) end; % for ( i = 1 : KOL_GLN) kol_gln_a = nom_ns - 1;
Функция ris_vis_sat
function ris_vis_sat(max_n, nj, ti_start, step_t, ris_vid)
%Имя функции:ris_vis_sat
%Функция строит графики спутников, видимых из точки наблюдения mm = step_t/3600.0; j_color = 0; color7(1:7) = ['b' 'g' 'r' 'c' 'm' 'y' 'k']; for i = 1 : max_n str1 = ris_vid(i, 1: nj ); k = find(str1); if k > 0 len = length(k); if len == 1 hPlot = plot(k(1):k(1),ris_vid(i,k(1):k(1)), s); set(hPlot, 'LineWidth', 3); hold on; else n0 = k(1); j_color = j_color + 1; if (j_color > 7 ) j_color = 1; end s = color7(j_color); for j = 2 : len if (k(j) > (k(j-1) + 1)) | (j == len) k0 = k(j - 1); hPlot = plot(n0*mm : mm : k0*mm,ris_vid(i,n0:k0), s); set(hPlot, 'LineWidth', 3); hold on; n0 = k(j); end; % if end; % for j = 2 : len end; % else end; % if k > 0 end; % for i = 1 : max_n
115
grid on hAxes=gca;
My = [0:2:25]; kol_x = ti_start + nj * step_t;
Mx = [0 : 2 : (kol_x - ti_start) / 3600]; set(hAxes, 'ytick', My,'xtick',Mx) ; hText = gca; set(hText,'FontSize',12,'FontName','TimesNewRoman') xlabel('Время наблюдения спутников ГЛОНАСС, час') ylabel('Номер навигационного спутника')
Функция 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);
Функция init_data
function[nom_sat_gln_a, trac_sat, sat_pos, satvis_gln] = init_data(kol)
%Имя функции: init_data
%Функция предназначена для инициализации массивов sat_pos.prn = 0; sat_pos.x = -1.0; sat_pos.y = -1.0; sat_pos.z = -1.0; sat_pos.vx = -1.0; sat_pos.vy = -1.0; sat_pos.vz = -1.0;
116
sat_pos.dvx = -1.0; sat_pos.dvy = -1.0; sat_pos.dvz = -1.0; sat_pos.range_m = -1.0; for( i = 1 : kol) nom_sat_gln_a(i) = 0; satvis_gln(i).el = -1; satvis_gln(i).az = -1; satvis_gln(i).r = 0; trac_sat(i).prn = 0; trac_sat(i).x = -1.0; trac_sat(i).y = -1.0; trac_sat(i).z = -1.0; trac_sat(i).vx = -1.0; trac_sat(i).vy = -1.0; trac_sat(i).vz = -1.0; trac_sat(i).range_m = -1.0; end;
Функция JD_from_epohi
function [jd] =JD_from_epohi( epoha, timeUTC);
%Имя функции:JD_from_epohi.m
%Функция вычисляет jd - количество дней от указанного года (epoha)
% до текущей даты, указанной в структуре timeUTC, представленной в виде
% (timeUTC.year, timeUTC.mon, timeUTC.day) jd0 = JD_epohi(epoha) + 1; % 12h, 1 den January
[day, day_year] = JD_data(timeUTC); jd = day - jd0;
Функция s0_Nut
function [S0] = s0_Nut( timeUTC, nut)
%Имя функции: s0_Nut
%Функция рассчитывает истинное или среднее звездное время на 0ч UTC
%Входные данные:
%timeUTC - дата, на которую требуется рассчитать истинное или среднее звездное время
%nut- признак (если nut= 0, то вычисляется среднее звездное время без учета нутации, иначе вычис- ляется
%истинное звездное время)
%Выходные данные:
%S0 - истинное или среднее звездное время на 0ч UTC jd2000 = 2451545; % 12h UTC 1 января
% Применить функцию JD_data
[jd, day_year] = JD_data( timeUTC); if (jd == NaN)
117
s0_mod = NaN; h = NaN; min = NaN; sec = NaN; fprintf('function s0_m - end0 \n'); return; end; jd = jd - 0.5; d = jd - jd2000; t = d / 36525.0; % 36525 - юлианский период 100 лет t2 = t * t; h1 = 24110.54841;
%h1=6.0*3600.0+41.0*60.0+50.54841;
% h2 = 236.555367908 * d; h2 = 8640184.812866 * t ; h3 = 0.093104 * t2; h4 = t2 * t * 6.2E-6; if ( nut == 0) na = 0; else na = utc_nut(t); end; s0_m = h1 + h2 + h3 - h4;
S0.s0_nut = s0_m + na;
S0.s0_m_mod = mod(s0_m, 86400); s0_day = floor(s0_m / 86400);
S0.s0_m_hour = S0.s0_m_mod / 3600.0;
S0.s0_m_hour = floor(S0.s0_m_mod / 3600); sec_min = S0.s0_m_mod - S0.s0_m_hour * 3600;
S0.s0_m_min = floor(sec_min / 60);
S0.s0_m_sec = sec_min - S0.s0_m_min * 60;
S0.s0_nut_mod = mod(S0.s0_nut, 86400); s0_day = floor(S0.s0_nut / 86400);
S0.s0_nut_hour = S0.s0_nut_mod / 3600.0;
S0.s0_nut_hour = floor(S0.s0_nut_mod / 3600); sec_min = S0.s0_nut_mod - S0.s0_nut_hour * 3600;
S0.s0_nut_min = floor(sec_min / 60);
S0.s0_nut_sec = sec_min - S0.s0_nut_min * 60;
Функция llh_to_eci
function [eci_llh, eci_xyz] = llh_to_eci(a, b, ti, time_s0, llh_loc) ;
%Имя функции:llh_to_eci
%Функция вычисляет позицию приемника в абсолютной геоцентрической системе координат (ECI)
%Входные данные:
%a, b-большая и малая полуоси земного эллипсоида (метр);
118

%ti- текущее время (секунды) ,
%time_s0- истиное звездное время,
%Структура llh_loc - координаты приемника; {
%llh_loc.lon-долгота (радиан);
%llh_loc.lat-широта (радиан);
%llh_loc.h-высота (метр);
%Выходные данные:
%Структура eci_llh - географические координаты приемника в абсолютной геоцентрической системе координат (ECI)
%eci_llh.lon - долгота (радиан);
%eci_llh.lat - широта (радиан);
%eci_llh.h - высота (метр);
%Структура eci_xyz- координаты приемника в абсолютной прямоугольной геоцентрической системе координат (ECI)
% eci_xyz.x - координата x;
% eci_xyz.y - координата y;
% eci_xyz.z - координата z;
OMEGA_Z = 0.7292115e-4; % угловая скорость вращения Земли, (рад/cek)
SEC_IN_RAD = 7.2722052166430e-005; % (PI / 43200.0 ) // Number radian eci_llh.lon = llh_loc.lon + ti * OMEGA_Z + SEC_IN_RAD * time_s0; eci_llh.lat = llh_loc.lat; eci_llh.h = llh_loc.h; eci_xyz = llh_to_ecef( a, b, eci_llh);
Функция init_satpos_gln
function [sat_pos] = init_satpos_gln() ;
%Имя функции:init_satpos_gln
%Функция предназначена для инициализации структуры sat_pos sat_pos.prn = 0; sat_pos.x = -1.0; sat_pos.y = -1.0; sat_pos.z = -1.0; sat_pos.vx = -1.0; sat_pos.vy = -1.0; sat_pos.vz = -1.0; sat_pos.dvx = -1.0; sat_pos.dvy = -1.0; sat_pos.dvz = -1.0; sat_pos.range_m = -1.0;
Функция gln_a_1
function [satpos_xyz_gln] = gln_a_1(ns, n0, ti_current, time_s0, alm_gln);
119

%Имя функции:gln_a_1
%Функция рассчитывает координаты и скорости спутников ГЛОНАСС в соответствии с интерфейс- ным
%контрольным документом ГЛОНАСС
%Входные данные:
%ns- номер спутника ,
%n0 - номер текущих суток внутри 4-х летнего периода (от ближайшего високосного года),
%ti_current - текущее время обсервации от начала дня,
%time_s0 - истинное звездное время в текущий момент обсервации,
%alm_gln - альманах спутников ГЛОНАСС
%Выходные данные:
%Структура satpos_xyz_gln
%satpos_xyz_gln.x - координата по оси x;
%satpos_xyz_gln.y- координата по оси y;
%satpos_xyz_gln.z- координата по оси z;
%satpos_xyz_gln.vx- скорость по оси x;
%satpos_xyz_gln.vy - скорость по оси y;
%satpos_xyz_gln.vz- скорость по оси z ;
% I_MID - Mean value of an inclination of a plane of orbit of a satellite
I_MID = 1.0995574287564; %(double)(PI * 63.0 / 180.0)
T_DR_MID = 43200.0; %Mean value a dragon of cycle time of a satellite
A_PZ90_KM = 6378.136; %(Km) Equatorial radius of the Earth
MU = 398600.44; % (Km^3/cek^2) constant of a gravitational
OMEGA_Z = 0.7292115e-4; %скорость вращения Земли (Angular speed of rotation of the Earth, рад/cek)
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; i_incl = alm_gln(ns).deltai; ecc2_1 = 1.0 - alm_gln(ns).ecc * alm_gln(ns).ecc;
% t_dr = T_DR_MID + alm_gln(ns).Tdr; t_dr = alm_gln(ns).Tdr;
120
n_dr = pi * 2 / t_dr;
%1. Вычисление полуоси a_n: a_n = semi_axis_1( 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_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; 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);
121
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 -...
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 *...
122

(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; omega_i = omega_0 + dd_omega; ii_incl = i_incl + ddi; ll_z = ma + alm_gln(ns).omegan + n_dr * tau + dd_ll;
123

% ll_z = ma + alm_gln(ns).omegan + n_dr * (ti - t_lambda_k) + dd_ll; mm_i = ll_z - w_i;
% 6. Coordinates and components of satellite velocity vector
% in absolute coordinate system OXaYaZa at the moment of ti time :
% Ei(n) = Mi + ei * sin Ei(n-1), Ei(0) = Mi, | Ei(n) - Ei(n-1) | <= 10-8, eps_kepler = 1E-8; ee_i = kepler( mm_i, ecc_i, eps_kepler);% ee_i - эксцентрическая аномалия
% alm_mca(ns).ek = ee_i; v_i = 2.0 * atan2(sqrt(1.0 + ecc_i) * tan(ee_i / 2.0) , sqrt(1.0 - ecc_i) ); u_i = v_i + w_i; ri = ai * (1.0 - ecc_i * cos(ee_i)) ; ecc_i2 = 1.0 - ecc_i * ecc_i; sqrt_mu_ai = sqrt( MU / ai); vr_i = sqrt_mu_ai * (ecc_i - sin(v_i) ) * ecc_i2; vu_i = sqrt_mu_ai * (1.0 + ecc_i * cos(v_i) ) * ecc_i2; cos_ui = cos(u_i); sin_ui = sin(u_i); cos_ii = cos(ii_incl); sin_ii = sin(ii_incl); cos_omega_i = cos(omega_i); sin_omega_i = sin(omega_i); x_dop = (cos_ui * cos_omega_i - sin_ui * sin_omega_i * cos_ii); y_dop = (cos_ui * sin_omega_i + sin_ui * cos_omega_i * cos_ii); satpos_xyz_gln.x = ri * x_dop; satpos_xyz_gln.y = ri * y_dop; satpos_xyz_gln.z = ri * sin_ui * sin_ii; satpos_xyz_gln.vx = vr_i * x_dop - vu_i * (sin_ui * cos_omega_i +... cos_ui * sin_omega_i * cos_ii); satpos_xyz_gln.vy = vr_i * y_dop - vu_i * (sin_ui * sin_omega_i -... cos_ui * cos_omega_i * cos_ii); satpos_xyz_gln.vz = vr_i * sin_ui * sin_ii + vu_i * cos_ui * sin_ii;
Функция satpos_eci_in_metr
function [satpos_eci, satpos_gln] = satpos_eci_in_metr(satpos_eci);
%Имя функции: satpos_eci_in_metr
%Функция преобразует координаты satpos_eci (структура), заданные в километрах в координаты
%satpos_eci, satpos_gln в метрах satpos_eci.x = satpos_eci.x * 1000.0; satpos_eci.y = satpos_eci.y * 1000.0;
124
satpos_eci.z = satpos_eci.z * 1000.0; satpos_eci.vx = satpos_eci.vx * 1000.0; satpos_eci.vy = satpos_eci.vy * 1000.0; satpos_eci.vz = satpos_eci.vz * 1000.0; satpos_gln.x = satpos_eci.x; satpos_gln.y = satpos_eci.y; satpos_gln.z = satpos_eci.z;
Функция init_satpos_gln
function [sat_pos] = init_satpos_gln() ;
%Имя функции:init_satpos_gln
%Функция предназначена для инициализации структуры sat_pos sat_pos.prn = 0; sat_pos.x = -1.0; sat_pos.y = -1.0; sat_pos.z = -1.0; sat_pos.vx = -1.0; sat_pos.vy = -1.0; sat_pos.vz = -1.0; sat_pos.dvx = -1.0; sat_pos.dvy = -1.0; sat_pos.dvz = -1.0; sat_pos.range_m = -1.0;
Функция eci_to_ecef
function [satpos_ecef] =eci_to_ecef(s0, ti, satpos_eci)
%Имя функции:eci_to_ecef
%Функция преобразования координат
%Входные данные: s0 - истинное звездное время в текущий момент обсервации ,
%ti - текущее время; satpos_eci
%Структура satpos_eci
%satpos_eci.x - координата x в абсолютной неподвижной системе координат (ECI);
%satpos_eci.y - координата y в абсолютной неподвижной системе координат (ECI);
%satpos_eci.z - координата z в абсолютной неподвижной системе координат (ECI);
%satpos_eci.vx - скорость vx в абсолютной неподвижной системе координат (ECI);
%satpos_eci.vy - скорость vy в абсолютной неподвижной системе координат (ECI);
% satpos_eci.vz - скорость vz в абсолютной неподвижной системе координат (ECI);
%Выходные данные:
% Структура satpos_ecef
%satpos_ecef.x - координата x в подвижной системе координат (ECEF);
%satpos_ecef.y - координата y в подвижной системе координат (ECEF);
%satpos_ecef.z - координата z в подвижной системе координат (ECEF);
125

%satpos_ecef.vx - скорость по оси x в подвижной системе координат (ECEF);
%satpos_ecef.vy - скорость по оси z в подвижной системе координат (ECEF);
%satpos_ecef.vz- скорость по оси z в подвижной системе координат (ECEF);
%Коэффициенты
% SEC_IN_RAD - коэффициент преобразования секунд в радианы
% s0(radian) = s0 (sek) * SEC_IN_RAD, where
% SEC_IN_RAD = 2 * pi / (24 * 3600) = pi / 43200
SEC_IN_RAD = pi / 43200;
OMEGA_Z = 0.7292115e-4; %( скорость вращения Земли (angular speed of rotation of the Earth, рад/cek) s_zv = s0 * SEC_IN_RAD + OMEGA_Z * ti; cos_s = cos(s_zv); sin_s = sin(s_zv); satpos_ecef.x = satpos_eci.x * cos_s + satpos_eci.y * sin_s; satpos_ecef.y = -satpos_eci.x * sin_s + satpos_eci.y * cos_s; satpos_ecef.z = satpos_eci.z; satpos_ecef.vx = satpos_eci.vx * cos_s + satpos_eci.vy * sin_s + OMEGA_Z * satpos_ecef.y; satpos_ecef.vy = -satpos_eci.vx * sin_s + satpos_eci.vy * cos_s - OMEGA_Z * satpos_ecef.x; satpos_ecef.vz = satpos_eci.vz;
Функция top_coord
function [satvis] = top_coord(prn, rec_llh, rec_xyz, nlo_xyz)
% Имя функции: top_coord
% Функция выполняет расчет топоцентрических координат объекта по заданным
%географическим (долгота, широта, высота) и геоцентрическим (x, y, z)
%координатам приемника, а также геоцентрическим координатам объекта(x, y, z)
%Входные данные:
%prn - номер спутника
%структура
% rec_llh.lat - широта (рад) приемника
%rec_llh.lon - долгота (рад) приемника
%rec_llh.h - высота (м) приемника
%структура
%прямоугольные геоцентрические координаты приемника (м)
% rec_xyz.x
%rec_xyz.y
%rec_xyz.z
%структура
126

% прямоугольные геоцентрические координаты объекта (м)
%nlo_xyz.x
%nlo_xyz.y
%nlo_xyz.z
%Выходные данные:
%структура satvis
% top.s - проекция вектора дальности на ось, направленную на Юг (South)
% top.e - проекция вектора дальности на ось, направленную на Восток (East)
% top.z - проекция вектора дальности на ось, направленную в Зенит
% top.daln - дальность до объекта
% top.az - угол азимута объекта
% top.el - угол видимости объекта rx = nlo_xyz.x - rec_xyz.x; ry = nlo_xyz.y - rec_xyz.y; rz = nlo_xyz.z - rec_xyz.z; r_sat = sqrt( rx * rx + ry * ry + rz * rz); r_rec = sqrt((rec_xyz.x)^2 + (rec_xyz.y)^2+ (rec_xyz.z)^2); top.r = r_sat; rx1 = rx; ry1 = ry; rz1 = rz; sin_lat = sin(rec_llh.lat); cos_lat = cos(rec_llh.lat); sin_lon = sin(rec_llh.lon); cos_lon = cos(rec_llh.lon);
% Projections of vector of range in topocentric coordinate system: top.e = -sin_lon * rx1 + cos_lon * ry1; top.s = cos_lon * sin_lat * rx1 + sin_lon * sin_lat * ry1 - cos_lat * rz1; top.z = cos_lat * cos_lon * rx1 + cos_lat * sin_lon * ry1 + sin_lat * rz1;
% azimut: отсчет по часовой стрелке от оси направленной на Север (N or -S) (-top.s) eps = 10e-10; if ( (abs(top.e) < eps) || (abs(top.s) < eps)) top.az = 0.0; else top.az = atan2(top.e,-top.s); end; if (top.az < 0.0) top.az = top.az + pi * 2; end;
% elevation: cos_el_top = (rec_xyz.x * rx + rec_xyz.y * ry + rec_xyz.z * rz) / (r_sat * r_rec); if ( cos_el_top >= 1.00 )
127
el = 0.0; else if ( cos_el_top <= -1.00 ) el = pi; else el = acos(cos_el_top); end; end; top.el = pi / 2.0 - el; satvis.el = top.el; satvis.az = top.az; satvis.r = r_sat;
Функция rewrite_satpos
function [result] = rewrite_satpos(nom, satpos)
%Имя функции:rewrite_satpos
%Функция предназначена для перезаписи структуры satpos в массив структур result result(nom).x = satpos.x; result(nom).y = satpos.y; result(nom).z = satpos.z; result(nom).vx = satpos.vx; result(nom).vy = satpos.vy; result(nom).vz = satpos.vz;
Функция utc_nut
function nut = utc_nut(t)
%Имя: utc_nut
%Функция предназначена для расчета нутации
%Входные данные:
%t =6.023472005475702e+002;
%Выходные данные:
%nut - нутация
R = 1296000; % ( 1r=360grad=1 296 000 cek)
RAD_SEK_ANGL = pi/(3600*180); t2 = t * t; t3 = t2 * t; l = 485866.733 + (1325.0 * R + 715922.633) * t + 31.310 * t2 + 0.064 * t3;%1.034807679476340e+012 l1 = 1287099.804 + (99 * R + 1292581.224) * t - 0.577 * t2 - 0.012 * t3; f = 335778.877 + (1342 * R + 295263.137) * t - 13.257 * t2 + 0.011 * t3; dd = 1072261.307 + (1236 * R + 1105601.328) * t - 6.891 * t2 + 0.019 * t3; omega = 450160.280 - (5 * R + 482890.539) * t + 7.455 * t2 + 0.008 * t3; eps0 = 84381.448 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3;
% eps0 = 84381.447996 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3;
128
eps_d = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'd','e'); eps_k = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'k','e'); eps = eps0 + eps_d + eps_k; cos_eps = cos(RAD_SEK_ANGL * eps ) / 15.0; d_fi = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'd', 'f'); k_fi = utc_nut_fi_eps(t, l, l1, f, dd, omega, 'k', 'f'); nut1 = d_fi * cos_eps; nut2 = k_fi * cos_eps;
% nut3 = 0.00264 * sin(omega) + 0.000063 * sin(2.0 * omega) nut3 = 0; nut1_dop = nut1; nut2_dop = nut2; nut3_dop = nut3; nut = nut1 + nut2 + nut3;
Функция utc_nut_fi_ep
function nut_fi_eps = utc_nut_fi_eps(t, l, l1, f, dd, omega, typ_nut, fi_eps)
% Имя функции:utc_nut_fi_eps
%Функция предназначена для расчета нутации по епсилон и фи
%Входные данные:
%t - значение параметра приведено в функции utc_nut;
%l - значение параметра приведено в функции utc_nut;
%l1 - значение параметра приведено в функции utc_nut;
%f -значение параметра приведено в функции utc_nut;
%dd - значение параметра приведено в функции utc_nut;
% оmega - значение параметра приведено в функции utc_nut;
%typ_nut - признак параметра приведен в функции utc_nut;
%fi_eps - признак параметра приведен в функции utc_nut;
%Выходные данные:
%nut_fi_eps - значение нутации по по епсилон и фи
% применить функцию koef
[koef_id, koef_abd, koef_ik, koef_abk] = koef;
RAD_SEK_ANGL = pi/(3600*180); if (typ_nut == 'd') n = 30; else n = 76; end; sum_a = 0; sum_b = 0; for i = 1 : n if (typ_nut == 'd')
129
s1 = koef_id(i,1) * l + koef_id(i,2) * l1 + koef_id(i,3) * f + koef_id(i,4) * dd + koef_id(i,5) * omega; if (fi_eps == 'f') a = koef_abd(i,1) * 1e-4; bt = koef_abd(i,2) * 1e-4; else a = koef_abd(i,3) * 1e-4; bt = koef_abd(i,4) * 1e-4; end; else s1 = koef_ik(i,1) * l + koef_ik(i,2) * l1 + koef_ik(i,3) * f + koef_ik(i,4) * dd + koef_ik(i,5) * omega; if (fi_eps == 'f') a = koef_abk(i,1) * 1e-4; bt = koef_abk(i,2) * 1e-4; else a = koef_abk(i,3) * 1e-4; bt = koef_abk(i,4) * 1e-4; end; end; if (fi_eps == 'f') sin_s1 = sin(RAD_SEK_ANGL * s1); sa = a * sin_s1; sb = bt * sin_s1; else cos_s1 = cos(RAD_SEK_ANGL * s1); sa = a * cos_s1; sb = bt * cos_s1; end; arg = RAD_SEK_ANGL * s1; sum_a = sum_a + sa; sum_b = sum_b + sb; end; nut_fi_eps = sum_a + sum_b * t;
Функция koef
function[koef_id, koef_abd, koef_ik, koef_abk] = koef
% Имя функции:koef
% Функция предназначена для ввода коэффициентов при расчете нутации koef_id = [ 0, 0, 0, 0, 1 ; % 1 0, 0, 0, 0, 2 ; % 2
-2, 0, 2, 0, 1 ; % 3 2, 0,-2, 0, 0 ; % 4
-2, 0, 2, 0, 2 ; % 5 130

1,-1, 0,-1, 0 ; % 6 0,-2, 2,-2, 1 ; % 7 2, 0,-2, 0, 1 ; % 8 0, 0, 2,-2, 2 ; % 9 0, -1, 0, 0, 0 ; % 10 % исправлено
0, 1, 2,-2, 2 ; % 11 0,-1, 2,-2, 2 ; % 12 0, 0, 2,-2, 1 ; % 13
-2, 0, 0,2, 0 ; % 14 % исправлено
0, 0, 2,-2, 0 ; % 15 0, 2, 0, 0, 0 ; % 16 0, 1, 0, 0, 1 ; % 17 0, 2, 2,-2, 2 ; % 18 0,-1, 0, 0, 1 ; % 19
-2, 0, 0, 2, 1 ; % 20 0,-1, 2,-2, 1 ; % 21 2, 0, 0,-2, 1 ; % 22 0, 1, 2,-2, 1 ; % 23 1, 0, 0,-1, 0 ; % 24 2, 1, 0,-2, 0 ; % 25 0, 0,-2, 2, 1 ; % 26 0, 1,-2, 2, 0 ; % 27 0, 1, 0, 0, 2 ; % 28
-1, 0, 0, 1, 1 ; % 29 0, 1, 2,-2, 0 ];% 30 koef_abd = [ -171996.0,-174.2, 92025.0, 8.9; % 1 2062.0, 0.2, -895.0, 0.5; % 2 46.0, 0.0, -24.0, 0.0; % 3 11.0, 0.0, 0.0, 0.0; % 4
-3.0, 0.0, 1.0, 0.0; % 5
-3.0, 0.0, 0.0, 0.0; % 6
-2.0, 0.0, 1.0, 0.0; % 7 1.0, 0.0, 0.0, 0.0; % 8
-13187.0, -1.6, 5736.0,-3.1; % 9
-1426.0, 3.4, 54.0,-0.1; % 10 % исправлено
-517.0, 1.2, 224.0,-0.6; % 11 217.0, -0.5, -95.0, 0.3; % 12 129.0, 0.1, -70.0, 0.0; % 13
-48.0, 0.0, 1.0, 0.0; % 14 % исправлено
-22.0, 0.0, 0.0, 0.0; % 15 17.0, -0.1, 0.0, 0.0; % 16 131

-15.0, 0.0, 9.0, 0.0; % 17
-16.0, 0.1, 7.0, 0.0; % 18
-12.0, 0.0, 6.0, 0.0; % 19
-6.0, 0.0, 3.0, 0.0; % 20
-5.0, 0.0, 3.0, 0.0; % 21 4.0, 0.0, -2.0, 0.0; % 22 4.0, 0.0, -2.0, 0.0; % 23
-4.0, 0.0, 0.0, 0.0; % 24 1.0, 0.0, 0.0, 0.0; % 25 1.0, 0.0, 0.0, 0.0; % 26
-1.0, 0.0, 0.0, 0.0; % 27 1.0, 0.0, 0.0, 0.0; % 28 1.0, 0.0, 0.0, 0.0; % 29
-1.0, 0.0, 0.0, 0.0]; % 30 koef_ik = [ 0, 0, 2, 0, 2; % 31 1, 0, 0, 0, 0; % 32 0, 0, 2, 0, 1; % 33 1, 0, 2, 0, 2; % 34 1, 0, 0,-2, 0; % 35
-1, 0, 2, 0, 2; % 36 0, 0, 0, 2, 0; % 37 1, 0, 0, 0, 1; % 38
-1, 0, 0, 0, 1; % 39
-1, 0, 2, 2, 2; % 40 1, 0, 2, 0, 1; % 41 0, 0, 2, 2, 2; % 42 2, 0, 0, 0, 0; % 43 1, 0, 2,-2, 2; % 44 2, 0, 2, 0, 2; % 45 0, 0, 2, 0, 0; % 46
-1, 0, 2, 0, 1; % 47
-1, 0, 0, 2, 1; % 48 1, 0, 0,-2, 1; % 49
-1, 0, 2, 2, 1; % 50 1, 1, 0,-2, 0; % 51 0, 1, 2, 0, 2; % 52 0,-1, 2, 0, 2; % 53 1, 0, 2, 2, 2; % 54 1, 0, 0, 2, 0; % 55 2, 0, 2,-2, 2; % 56 0, 0, 0, 2, 1; % 57 0, 0, 2, 2, 1; % 58 132

1, 0, 2,-2, 1; % 59 0, 0, 0,-2, 1; % 60 1,-1, 0, 0, 0; % 61 2, 0, 2, 0, 1; % 62 0, 1, 0,-2, 0; % 63 1, 0,-2, 0, 0; % 64 0, 0, 0, 1, 0; % 65 1, 1, 0, 0, 0; % 66 1, 0, 2, 0, 0; % 67 1,-1, 2, 0, 2; % 68
-1,-1, 2, 2, 2; % 69
-2, 0, 0, 0, 1; % 70 3, 0, 2, 0, 2; % 71 0,-1, 2, 2, 2; % 72 1, 1, 2, 0, 2; % 73
-1, 0, 2,-2, 1; % 74 2, 0, 0, 0, 1; % 75 1, 0, 0, 0, 2; % 76 3, 0, 0, 0, 0; % 77 0, 0, 2, 1, 2; % 78
-1, 0, 0, 0, 2; % 79 1, 0, 0,-4, 0; % 80
-2, 0, 2, 2, 2; % 81
-1, 0, 2, 4, 2; % 82 2, 0, 0,-4, 0; % 83 1, 1, 2,-2, 2; % 84 1, 0, 2, 2, 1; % 85
-2, 0, 2, 4, 2; % 86
-1, 0, 4, 0, 2; % 87 1,-1, 0,-2, 0; % 88 2, 0, 2,-2, 1; % 89 2, 0, 2, 2, 2; % 90 1, 0, 0, 2, 1; % 91 0, 0, 4,-2, 2; % 92 3, 0, 2,-2, 2; % 93 1, 0, 2,-2, 0; % 94 0, 1, 2, 0, 1; % 95
-1,-1, 0, 2, 1; % 96 0, 0,-2, 0, 1; % 97 0, 0, 2,-1, 2; % 98 0, 1, 0, 2, 0; % 99 1, 0,-2,-2, 0; % 100 133

0,-1, 2, 0, 1; % 101 1, 1, 0,-2, 1; % 102 1, 0,-2, 2, 0; % 103 2, 0, 0, 2, 0; % 104 0, 0, 2, 4, 2; % 105 0, 1, 0, 1, 0]; % 106 koef_abk = [-2274.0, -0.2, 977.0, -0.5; % 31 712.0, 0.1, -7.0, 0.0; % 32
-386.0, -0.4, 200.0, 0.0; % 33
-301.0, 0.0, 129.0, -0.1; % 34
-158.0, 0.0, -1.0, 0.0; % 35 123.0, 0.0, -53.0, 0.0; % 36 63.0, 0.0, -2.0, 0.0; % 37 63.0, 0.1, -33.0, 0.0; % 38
-58.0, -0.1, 32.0, 0.0; % 39
-59.0, 0.0, 26.0, 0.0; % 40
-51.0, 0.0, 27.0, 0.0; % 41
-38.0, 0.0, 16.0, 0.0; % 42 29.0, 0.0, -1.0, 0.0; % 43 29.0, 0.0, -12.0, 0.0; % 44
-31.0, 0.0, 13.0, 0.0; % 45 26.0, 0.0, -1.0, 0.0; % 46 21.0, 0.0, -10.0, 0.0; % 47 16.0, 0.0, -8.0, 0.0; % 48
-13.0, 0.0, 7.0, 0.0; % 49
-10.0, 0.0, 5.0, 0.0; % 50
-7.0, 0.0, 0.0, 0.0; % 51 7.0, 0.0, -3.0, 0.0; % 52
-7.0, 0.0, 3.0, 0.0; % 53
-8.0, 0.0, 3.0, 0.0; % 54 6.0, 0.0, 0.0, 0.0; % 55 6.0, 0.0, -3.0, 0.0; % 56
-6.0, 0.0, 3.0, 0.0; % 57
-7.0, 0.0, 3.0, 0.0; % 58 6.0, 0.0, -3.0, 0.0; % 59
-5.0, 0.0, 3.0, 0.0; % 60 5.0, 0.0, 0.0, 0.0; % 61
-5.0, 0.0, 3.0, 0.0; % 62
-4.0, 0.0, 0.0, 0.0; % 63 4.0, 0.0, 0.0, 0.0; % 64
-4.0, 0.0, 0.0, 0.0; % 65
-3.0, 0.0, 0.0, 0.0; % 66 134

3.0, 0.0, 0.0, 0.0; % 67
-3.0, 0.0, 1.0, 0.0; % 68
-3.0, 0.0, 1.0, 0.0; % 69
-2.0, 0.0, 1.0, 0.0; % 70
-3.0, 0.0, 1.0, 0.0; % 71
-3.0, 0.0, 1.0, 0.0; % 72 2.0, 0.0, -1.0, 0.0; % 73
-2.0, 0.0, 1.0, 0.0; % 74 2.0, 0.0, -1.0, 0.0; % 75
-2.0, 0.0, 1.0, 0.0; % 76 2.0, 0.0, 0.0, 0.0; % 77 2.0, 0.0, -1.0, 0.0; % 78 1.0, 0.0, -1.0, 0.0; % 79
-1.0, 0.0, 0.0, 0.0; % 80 1.0, 0.0, -1.0, 0.0; % 81
-2.0, 0.0, 1.0, 0.0; % 82
-1.0, 0.0, 0.0, 0.0; % 83 1.0, 0.0, -1.0, 0.0; % 84
-1.0, 0.0, 1.0, 0.0; % 85
-1.0, 0.0, 1.0, 0.0; % 86 1.0, 0.0, 0.0, 0.0; % 87 1.0, 0.0, 0.0, 0.0; % 88 1.0, 0.0, -1.0, 0.0; % 89
-1.0, 0.0, 0.0, 0.0; % 90
-1.0, 0.0, 0.0, 0.0; % 91 1.0, 0.0, 0.0, 0.0; % 92 1.0, 0.0, 0.0, 0.0; % 93
-1.0, 0.0, 0.0, 0.0; % 94 1.0, 0.0, 0.0, 0.0; % 95 1.0, 0.0, 0.0, 0.0; % 96
-1.0, 0.0, 0.0, 0.0; % 97
-1.0, 0.0, 0.0, 0.0; % 98
-1.0, 0.0, 0.0, 0.0; % 99
-1.0, 0.0, 0.0, 0.0; % 100
-1.0, 0.0, 0.0, 0.0; % 101
-1.0, 0.0, 0.0, 0.0; % 102
-1.0, 0.0, 0.0, 0.0; % 103 1.0, 0.0, 0.0, 0.0; % 104
-1.0, 0.0, 0.0, 0.0; % 105 1.0, 0.0, 0.0, 0.0]; % 106
Функция semi_axis_1
135

% double semi_axis(double t_dr, double i_incl, double ecc, double omega_n) function [a_npp] = semi_axis_1(t_dr, i_incl, ecc, omega_n);
%Имя функции:semi_axis_1
%Функция вычисляет радиус орбиты спутника ГЛОНАСС в соответствии с интерфейсным контроль- ным
%документом ГЛОНАСС
%Входные данные:
%t_dr - драконический период обращения спутника ГЛОНАСС (секунды)
%i_incl - наклонение орбиты спутника ГЛОНАСС (радиан)
%ecc - эксцентриситет
%omega_n - аргумент перигея орбиты спутника ГЛОНАСС (радиан)
%Выходные данные:
%a_npp - большая полуось орбиты спутника ГЛОНАСС (километр)
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; 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
136
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;
Функция kepler
function [Ek] = kepler(Mk, ecc, eps);
%Имя функции:kepler
%Функция предназначена для решения уравнения Кеплера
% eps = 1.0E-15; y = ecc * sin(Mk); x1 = Mk; x = y; for k = 0 : 15 x2 = x1; x1 = x; y1 = y; y = Mk - (x - ecc * sin(x)); if (abs(y - y1) < eps) break end x = (x2 * y - x * y1) / (y - y1); end %k
Ek = x;
Функция map
function map(N)
%Имя функции:map
%Применена функция MatLab для внесения в графики орбитального движения изображения Земли load('topo.mat','topo','topomap1');
[x,y,z] = sphere(50); cla reset
%axis square off props.AmbientStrength = 0.1;
137
props.DiffuseStrength = 1; props.SpecularColorReflectance = .5; props.SpecularExponent = 20; props.SpecularStrength = 1; props.FaceColor= 'texture'; props.EdgeColor = 'none'; props.FaceLighting = 'phong'; props.Cdata = topo; surface(x*N,y*N,z*N,props); light('position',[-1 0 1]); light('position',[-1.5 0.5 -0.5], 'color', [.6 .2 .2]); view(3) ; grid on
4.3.5 Примеры выполнения комплекса программ ORBITA_GL_NAVIOR
Некоторые результаты, иллюстрирующие работу программ, изображены на рис. 4.8 - рис. 4.11.
Рис. 4.8. Орбиты спутников ГЛОНАСС в системе координат ECEF
138

Рис. 4.9. Орбиты спутников ГЛОНАСС в системе координат ECI
Рис. 4.10. Орбита спутника ГЛОНАСС № 24 за 7 суток
139

0 2
4 6
8 10 12 14 16 18 20 22 24 0
2 4
6 8
10 12 14 16 18 20 22 24
Время наблюдения спутников ГЛОНАСС, час
Но м
е р
нав и
гац и
онн о
го сп ут ник а
Рис. 4.11. Видимость спутников ГЛОНАСС
Программы имеют и другие возможности в зависимости от задач, которые может поставить специалист, исследующий орбитальное движение спутников.
140

141
1   ...   4   5   6   7   8   9   10   11   ...   14

РАЗДЕЛ 5 Преобразование данных навигационных спутников
5.1 Преобразование данных альманаха приемника СН 4701 в формат YUMA
5.1.1 Краткие сведения из теории
В спутниковых радионавигационных системах ГЛОНАСС и GPS эфемериды рассчи- тываются по принципиально разным формулам. В то же время весьма целесообразным является применение в задачах прогнозирования навигационных сеансов использовать одинаковый формат данных альманаха. Отсюда вытекает задача преобразования (конвер- тирования) альманаха системы ГЛОНАСС в широко распространенный формат YUMA без потери данных об альманахе спутников ГЛОНАСС.
Программы представленные в данном подразделе выполняют преобразование аль- манаха спутников GPS и ГЛОНАСС в формат YUMA, который воспринимается програм- мой PLANNING и обеспечивает возможность планирования сеансов навигационных оп- ределения при применении спутниковых систем GPS и ГЛОНАСС. Программный ком- плекс подготовлен в точном соответствии с интерфейсными контрольными документами [
].
Рис. 5.1.
Блок- схема программы преобразования альманаха ГЛОНАСС/GPS в формат YUMA convers_АLM_GLN_YUMA_GPS.m
Входные данные 'In_dat\002.txt'
Функция «Чтение альманаха»
Read_GL_Alm
Функция преобразования дня при- вязки данных альманаха
Gln_data_from_NA
Функция read_Alm(Dat)
Функция write_GPS_alm
Функция
GLN_satfind_YUMA
Функция
JD_epohi
Функция
Gln_data_from_NA
Функция
JD_from_epohi
Функция s0_Nut
Функция
JD_data
Функция week_GLONAS_gps
Функция gln_a
Функция semi_axis
Выходные данные:
AlmGG.yum

Для системы ГЛОНАСС Программы, записанные на CD-диск, находятся в папке
ALM_CH4701_V3. Схема построения программного комплекса изображена на рис. 5.1.
Основные комментарии и пояснения к программам приведены в листингах.
Цель лабораторной работы преобразование (конвертирование) данных альманаха навигационных спутников в формат YUMA.
5.1.2 Лабораторная работа 5. 1 «Конвертирование данных альманаха GPS и ГЛО-
НАСС в формат » YUMA
Пакет программ для выполнения лабораторной работы расположен в папке
ALM_CH4701_V3. Для выполнения работы в качестве входных данных потребуется аль- манах спутников GPS и ГЛОНАСС, который можно получить с нескольких сайтов или с навигационного приемника, если приемник имеется в лаборатории.
Рекомендуется следующий порядок выполнения лабораторной работы.
1. Запишите альманах, с приемника CH4701 в папку In_dat. Если не имеется возможно- сти получить данные альманаха с приемника, то можно воспользоваться альманахом с именем 002.txt из приложения к данной книге.
2. Создайте папку ALM_CH4701_V3_My и скопируйтев ее все программы из папки
ALM_CH4701_V3 и папку In_dat.
3. Запустите MatLab.
4. Обратитесь к папке ALM_CH4701_V3_My , откройте ее, изучите функции, файлы и программные процедуры по комментариям и выполните задания 1 и 2.
5. Задание 1. Откройте файл convers_ALM_GLN_YUMA_GPS.m в папке
ALM_CH4701_V3_My и присвойте переменной Dat имя альманаха из папки In_dat.
Выполните файл. Результат выполнения запишется в файл AlmGG.yum. Это и будет преобразованный в формат YUMA альманах спутников GPS и ГЛОНАСС.
6. Задание 2. Загрузите альманах AlmGG.yum в программу прогнозирования PLANNING и убедитесь в его работоспособности.
5.1.3 Задание для самоподготовки
Ознакомтесь и освойте работу с программой PLANNING
5.1.4 Функции и файлы из папки ALM_CH4701_V3
Файл convers_АLM_GLN_YUMA_GPS.m
142

%Имя файла:convers_АLM_GLN_YUMA_GPS.m
%от conversion- преобразование; программа преобразует данные альманаха GPS и ГЛОНАСС, полученные
%с приемника "СН 4701" (производитель ГП "ОРИЗОН-НАВИГАЦИЯ" в формат YUMA
%Входные данные: альманах- записан в папке с именем In_dat, например,
%Dat = 'In_dat\GLN_all_23_10_06_17_45.alm';
%Выходные данные:записываются в файл с именем 'AlmGG.yum'
%Dat = 'In_dat\GLN_all8_11_06.alm';
%Dat = 'In_dat\GLN06_11_06.alm';
%Dat = 'In_dat\GG_24_03_03.alm';% GG24
Dat = 'In_dat\002.txt';% CN
%Имя файла для записи
Name = 'AlmGG.yum';
%Открытие файла для записи fid =fopen(Name,'wt');
% Чтение альманаха
[alm_GPS,max_kol_GPS,alm_gln,max_kol_GLN] = read_Alm(Dat);
%Запись в файл альманаха GPS write_GPS_alm(fid, alm_GPS);
%Полуоси земного эллипсоида
A_WGS84_M = 6378137.0 ; % WGS-84 ellipsoid parameters
A_PZ90_M = 6378136.0; %6 378 136 м - Equatorial radius of the Earth - ,большая полуось эллипсоида
A_PZ90_M = A_WGS84_M;
%Смещение времени GPS от UTC в секундах
%dt_lsf = 13; leap_year = 2004;% Високосный год timeUTC_leap.year = leap_year; timeUTC_leap.mon = 0; timeUTC_leap.day = 0;
GLN_satfind_YUMA( A_PZ90_M, timeUTC_leap, alm_gln, fid);
Функция read_Alm
function [alm,max_kol_GPS,alm_gln,max_kol_GLN] = read_Alm(Dat)
%Имя функции:read_Alm
%Функция читает данные альманаха GPS и ГЛОНАСС, полученные
%с приемника "СН 4701" (производитель ГП "ОРИЗОН-НАВИГАЦИЯ"
%Dat = '001.txt'; i=0; for i=1:61 alm(i).ID = 0; alm(i).Health=63;
143
alm_gln(i).ID=0; alm_gln(i).Health=255; end; fid =fopen(Dat,'rt'); max_kol_GLN = 0; max_kol_GPS = 0; while not(feof(fid)) str1=fscanf(fid,'%s',1); % GPS | GLN len1 = length(str1); if len1 == 4 % GPS if str1 == 'GPS:' max_kol_GPS = max_kol_GPS + 1; str1=fscanf(fid,'%c',4);
ID=fscanf(fid,'%d',1); alm(ID).ID = ID; str2=fscanf(fid,'%c',13);
Health=fscanf(fid,'%d',1); alm(ID).Health=Health; str3=fscanf(fid,'%c',19); eccentricity=fscanf(fid,'%g',1); alm(ID).ecc =eccentricity; str4=fscanf(fid,'%c',27); deltai=fscanf(fid,'%g',1) ; %i0 alm(ID).deltai= deltai ; %i0 str5=fscanf(fid,'%c',30);
OMEGADOT=fscanf(fid,'%g',1)*1000; alm(ID).OMEGADOT=OMEGADOT; str6=fscanf(fid,'%c',37);
A05=fscanf(fid,'%g',1); alm(ID).A05=A05; str7=fscanf(fid,'%c',32); omega0 =fscanf(fid,'%g',1); alm(ID).omega0 =omega0; str8=fscanf(fid,'%c',26); omega =fscanf(fid,'%g',1); alm(ID).omega=omega;
144
str9=fscanf(fid,'%c',19);
M0=fscanf(fid,'%g',1); alm(ID).M0=M0; str10=fscanf(fid,'%c',11);
Af0m=fscanf(fid,'%g',1)/1000; alm(ID).Af0m=Af0m; str11=fscanf(fid,'%c',10);
Af01=fscanf(fid,'%g',1); alm(ID).Af01=Af01; str12=fscanf(fid,'%c',11);
Af0=fscanf(fid,'%g',1); alm(ID).Af0=Af0; str13=fscanf(fid,'%c',19);
Week=fscanf(fid,'%g',1)+1024; alm(ID).Week=Week; str14=fscanf(fid,'%c',27);
TOA=fscanf(fid,'%d')/1000; alm(ID).TOA =TOA; end;% if str1 == 'GPS' else % GLONASS len1 = length(str1); if len1 == 8 max_kol_GLN = max_kol_GLN+1; str2=fscanf(fid,'%c',4);
ID=fscanf(fid,'%d',1); alm_gln(ID).ID = ID; str3=fscanf(fid,'%c',15);
Health=fscanf(fid,'%d',1); alm_gln(ID).Health=Health; str4=fscanf(fid,'%c',19);
Hn=fscanf(fid,'%d',1); alm_gln(ID).Hn=Hn;
145
str5=fscanf(fid,'%c',34);
TaUGL=fscanf(fid,'%g',1); alm_gln(ID).tau_n=TaUGL/1000; str6=fscanf(fid,'%c',52);
LambdaN=fscanf(fid,'%g',1); alm_gln(ID). LambdaN=LambdaN; str7=fscanf(fid,'%c',24); deltai=fscanf(fid,'%g',1); alm_gln(ID).deltai=deltai; str8=fscanf(fid,'%c',21); eccentricity=fscanf(fid,'%g',1); alm_gln(ID).ecc =eccentricity; str9=fscanf(fid,'%c',39); omegan=fscanf(fid,'%g',1); alm_gln(ID).omegan =omegan; str10=fscanf(fid,'%c',48);
TLambdaN=fscanf(fid,'%g',1); alm_gln(ID).TLambdaN =TLambdaN/1000-10800; str11=fscanf(fid,'%c',27);
Tdr=fscanf(fid,'%g',1); alm_gln(ID).Tdr =Tdr/1000; str12=fscanf(fid,'%c',48); dTdr=fscanf(fid,'%g',1); alm_gln(ID).dTdr =dTdr;%?единицы измерения str13=fscanf(fid,'%c',16);
Na=fscanf(fid,'%g',1); alm_gln(ID).Na =Na; end % if len1 == 4 % GPS end % while end fclose(fid);
146

Смотрите также файлы