Функция 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
Скачать 4.49 Mb.