Модель движения навигационных спутников GPS и ГЛОНАСС

Цель: изучение орбитального движения навигационных спутников GPS и ГЛОНАСС с последующей визуализацией спутников на любой момент времени.

Теоретические сведения.

Развитие систем спутниковой радионавигации в настоящее время и в ближайшей перспективе ориентировано на применении нескольких систем. При этом в одном навигационном приемнике должны обрабатываться сигналы спутников систем, которые находятся в зоне видимости.

Программа рассчитывает углы азимута и места навигационных спутников GPS и ГЛОНАСС по данным совместного для обеих систем альманаха в формате YUMA. Файл альманаха мож- но получить на сайте Национального авиационного университета или запросить по элек- тронной почте cnsatm@nau. edu.ua.

 

Выполнение работы

1. Создайте папку Vsion_GLONASS_GPS_My.

2. Скопируйте в папку Vsion_GLONASS_GPS_My функции и файлы.

3. Запишите в папку Vsion_GLONASS_GPS_My альманах спутников GPS и ГЛОНАСС в формате YUMA

4. Разберите тексты функций и файлов в папке Vsion_GLONASS_GPS_My. Откройте файл Vision_GLONASS_GPS.m и внесите в него входные данные: имя альманаха и время, на которое выполняется моделирование (в комментариях к файлу это указано). Выполните файл, и результаты визуализации с полярной диаграммы внесите в отчет.

5. Измените время наблюдения на 6 часов, выполните файл Vision_GLONASS_GPS.m. Занесите результаты наблюдения в отчет и объясните при-чины изменения конфигурации спутников.

6. В файле Vision_GLONASS_GPS.m. измените дату наблюдения на два дня, сохранив время наблюдения, и выполните файл. Занесите результаты выполнения про- граммы в отчет и объясните причины, по которым спутники ГЛОНАСС изменили свое местоположения, а спутники GPS почти остались на том же месте.

Рис.1 Видимость спутников GPS и ГЛОНАСС

Вопросы:

1. Чему равен период обращения спутников GPS?.

2. Чему равен период обращения спутников ГЛОНАСС?

3. Применяя второй закон Кеплера и данные альманаха определите периоды обращения спутников GPS и спутников ГЛОНАСС.

Решение навигационной задачи

Лабораторная работа № 11

Решение навигационной задачи

Цель работы: изучение метода расчета координат потребителя по данным, полученным с навигационных спутников GPS, ГЛОНАСС и измеренным псевдодальностям

Теоретические сведения.

Главной задачей спутникового навигационного приемника является определения координат. Для решения этой задачи требуется обнаружить сигналы спутников, перейти в режим слежения за сигналами, принять и декодировать данные, поступаюшие в сообщениях спутников, измерить расстояния до спутников, обработать всю информацию и решить навигационную задачу.

Выполнение работы

1. Создайте папку Координаты приемника_My, запишите в нее файл positionV0.m, изучите программные процедуры и комментарии.

2. Последовательно выполните файл positionV0.m для расчета координат при участии в расчетах 15, 14, 13, 9, 8, 7, 6 спутников из разных созвездий. Проанализируйте полученные результаты. Расчеты позиции и результаты анализа запишите в отчет.

3. Добавьте к псевдодальностям до спутников, которые используются для решения навигационной задачи одинаковое приращение. Исполните файл. Проанализируйте и объясните полученный результат. Данные занесите в отчет.

Вопросы:

1. Какие входные данные требуются для решения навигационной задачи?

2. Какое минимальное количество навигационных спутников одного созвездия требуется для решения навигационной задачи?

3. Каким образом определяется псевдодальность до навигационного спутника?

4. Как определяется пространственный геометрический фактор (PDOP)?

5. Какие начальные условия могут быть при решении навигационной задачи?

Приложения

Функции и файлы из папки ORBITA_GPSv1

 

Функция ECEFLLH

function [Rx,Ry,Rz] = ECEFLLH(lon,lat,hr)

%Имя функции: ECEFLLH.m

%Функция выполняет преобразование координат

%Входные данные:lon-долгота,lat-широта,h-высота;a,b-большая

%и малая полуоси эллипсоида

%Выходные данные:Rx,Ry,Rz- координаты в ECEF

%Для WGS-84

a=6378137.0;

b=6356752.314;

n=a*a/sqrt(a*a*cos(lat)*cos(lat)+b*b*sin(lat)*sin(lat));

Rx=(n+hr)*cos(lat)*cos(lon);

Ry=(n+hr)*cos(lat)*sin(lon);

Rz=(b*b/(a*a)*n+hr)*sin(lat);

 

Функция LLHECEF

function [lons,lats,hrs] = LLHECEF(Xk,Yk,Zk)

%Имя функции: LLHECEF.m

%Функция выполняет преобразование координат.

%Входные данные:Rx,Ry,Rz- координаты в ECEF

%Выходные данные:lon-долгота,lat-широта,h-высота

%a,b-большая и малая полуоси эллипсоида

a=6378137.0;

b=6356752.314;

xy = sqrt(Xk*Xk + Yk*Yk);

thet = atan(Zk*a/(xy*b));

esq = 1.0-b*b/(a*a);

epsq = a*a/(b*b)-1.0;

lats = atan((Zk+epsq*b*(sin(thet)^3))/(xy-esq*a*(cos(thet)^3)));

lons = atan2(Yk,Xk);%!

if lons < 0

lons = 2*pi + lons;

end;

n = a*a/sqrt(a*a*cos(lats)*cos(lats) + b*b*sin(lats)*sin(lats));

hrs = xy/cos(lats)-n;

end

 

Функция Tim

function [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s)

%Имя функции:Tim.m

%Функция Tim.m работает совместно со стандартной функцией MatLabDAYSDIF.m и рассчитывает:

%week- текущую GPS-неделю, modeweek- модифицированную GPS неделю, d- количество дней,

%dweek- день недели,weeks- время GPS

%Входные данные:d2='10/23/2007' - 'номер месяца/номер дня месяца/номер года', h=23.0 -часы,

%;min=59.0- минуты, s=59.0- секунды на которые рассчитываются выходные данные

%d2='10/35/2003';h=23.0;min=59.0;s=59.0;

week = floor(DAYSDIF('1/6/1980',d2,3)/7);% текущая GPS-неделя

modeweek=week-1024;% модифицированная GOS-неделя

d = DAYSDIF('1/6/1980',d2,3);%количество дней

dweek=fix(d-week*7);%номер дня недели (нулевой день-воскресенье)

weeks=(dweek)*24*60*60+h*60*60+min*60+s;% время GPS в неделе (секунды)

 

 

Файл PR_Tim.m

%Пример применения функции Tim.m

%Имя m- файла: PR_Tim.m

d2= '10/13/2006';

h=22.0;

min=40.0;

s=11.0;

[week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s);

[week,modeweek,d,dweek,weeks] %=1396 372 9777 5 513611-результат расчета;

%1396-неделя GPS, отсчитываемая с ночи с 5 на 6 января 1980 года, 372=1396-1024- модифициро-ванная

% неделя GPS, 9777- количество дней прошедших с 6 января 1980 года, 5-пятый день недели (пятни-ца),

%считая с понедельника, 513611-количество секунд от начала текущей недели.

d = DAYSDIF('1/6/1980',d2,3);%функция MatLab

 

Функция Yuma_GPS_Alm1

 

function [alm,max_kol] = Yuma_GPS_Alm1(Dat)

%Имя функции:Yuma_GPS_Alm1.m

%Функция читает данные альманаха, записанные в формате YUMA

%Входные данные записываются в переменную Dat, например,

%Dat='Имя файла альманаха в формате YUMA'

%Выходные данные:1. Численные значения альманаха спутников GPS, представляемые

%в виде структуры в переменной alm=[%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)], где

% цифра в скобках обозначает порядковый номер параметра альманаха в формате YUMA

%Для чтения альманаха в m-файл фикция записывается в виде [alm,max_kol] = Yuma_GPS_Alm1(Dat).

%2. Количество спутников GPS записывается в переменную max_kol.

fori=1:32% цикл

alm(i).ID = 0;% обнуление массива

alm(i).Health=63;% обнуление массива

end;

fid =fopen(Dat,'rt');% открыть файл для чтения

%чтение данных из файла

max_kol = 0;

while not(feof(fid))

s1=fscanf(fid,'%s',6);

if not(feof(fid))

lenstr = length(s1);

while (fscanf(fid,'%s',1) == '*')

end

str1 = fscanf(fid,'%s',1);

lenstr = length(str1);

n_sv = sscanf(str1,'%d');

strID = str1(1:lenstr);

ID = sscanf(strID,'%d');

alm(ID).ID = ID;

t_2=fscanf(fid,'%s',1);

alm(ID).Health=fscanf(fid,'%d',1);

t_3=fscanf(fid,'%s',1);

alm(ID).e = fscanf(fid,'%g',1);

t_4=fscanf(fid,'%s',3);

alm(ID).TOA =fscanf(fid,'%g',1);

t_5=fscanf(fid,'%s',2);

alm(ID).deltai=fscanf(fid,'%g',1);%i0

t_6=fscanf(fid,'%s',4);

alm(ID).OMEGADOT=fscanf(fid,'%g',1);

while not(fscanf(fid,'%c',1) == ':')

end

alm(ID).A05=fscanf(fid,'%g',1);

t_8=fscanf(fid,'%s',4);

alm(ID).omega0 =fscanf(fid,'%g',1);

t_9=fscanf(fid,'%s',3);

alm(ID).omega=fscanf(fid,'%g',1);

t_10=fscanf(fid,'%s',2);

alm(ID).M0=fscanf(fid,'%g',1);

t_11=fscanf(fid,'%s',1);

alm(ID).Af0=fscanf(fid,'%g',1);

t_12=fscanf(fid,'%s',1);

alm(ID).Af1=fscanf(fid,'%g',1);

t_13=fscanf(fid,'%s',1);

alm(ID).Week=fscanf(fid,'%g',1);

max_kol = max_kol+1;

end

end

fclose(fid)

 

Файл Orbita_GPS.m

clear all

 

%Имя m-файла:Orbita_GPS.m

%Программа рассчитывает орбиты навигационных спутников GPS

%Входные данные:

%файл альманаха в формате Yuma,имя файла альманаха присваивается

%переменной "Dat",например,Dat = 'имя файла альманаха';

%данные о начале отсчета "d2",d2='месяц/день/год';h=час;min=минута;s=секунда;

%координаты позиции приемника –lat (широта в радианах),lon (долгота в радианах,

%hr (высота в метрах);

%шаг, с каким будут рассчитываться параметры орбит (step,секунды);

%количество точек (L), в которых будут рассчитываться параметры орбит

%L=12*3600/step,L читается так: количество часов (например,12)

%число секунд в часе (3600) деленное на шаг (step)

%В программе применяются функции: Yuma_GPS_Alm1.m- считывание данных альманаха,

%заданного в формате YUMA; ECEFLLH.m, LLHECEF.m - преобразование координат;Tim.m- расчет времени;

%Постоянные:

%скорость вращения Земли

OMEGAeDOT=7.2921151467e-005;

%или

%OMEGAeDOT=0;

mu=398600500000000;

F_CONST = 4.442807633E-10;

%Задание цветов для графики

j_color = 0;

color6(1:16) = [':' 'k' '.' 'r' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h'];

 

%Входные данные

Dat = 'almanac_yuma_week0371_589824.txt';

d2='10/06/2006';h=13.0;min=8.0;s=55.0;

lat = 0.88032730015257; %50 град; 26 мин.; 20.54 с

lon = 0.53109641675259;%30 град; 25 мин.; 46.4995 с

hr=184;%высота в метрах

step=300;

L=(10*3600)/step;

 

%Чтение альманаха

 

[alm,max_kol] = Yuma_GPS_Alm1(Dat);

kol = 0;

for i = 1: max_kol

id=alm(i).ID;

if id > 0

kol = kol + 1;

 

nom_ns(kol) = id;

end

nom_ns;%номер навигационного спутника

end

 

%Преобразование координат

 

[Rx,Ry,Rz] = ECEFLLH(lon, lat,hr);

%Rx=0;Ry=0;,Rz=0;%центр масс Земли

%Выбор спутников:

%для выбора спутников вводится параметры kol-количеество спутников для

%исследования и номера спутников, например, kol =4; nom_ns(1:kol) = [3 6 7 31],

%такая запись обозначает, что исследуются (рассчитываются орбиты 4 спутников

%с номерами 3,6,7,31; количество номеров спутников должно совпадать с kol

%Варианты (можно любые другие)

%kol =9

%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10];

%kol =5

%nom_ns(1:kol) = [1 13 14 26 29]; %1 спутники орбиты 1

%kol =5

%nom_ns(1:kol) = [2 5 22 28 30]; %2 спутники орбиты 2

%kol =4

%nom_ns(1:kol) = [3 6 7 31]; %3 спутники орбиты 3

%kol =5

%nom_ns(1:kol) = [4 11 15 17 24 ]; %4 спутники орбиты 4

%kol =4

%nom_ns(1:kol) = [8 9 25 27]; %5 спутники орбиты 5

 

%nom_ns(1:kol) = [10 18 20 21 23]; %6 спутники орбиты 6

%kol =29

%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31];

%kol =14;

%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10 11 13 14 15 16];% 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31];

kol =2;

nom_ns(1:kol)=[1 3 ];

str1 = num2str(nom_ns(1:kol));

for k = 1: kol

i = nom_ns(k);

%Начало отсчета текущего времени

[week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s);

%Расчет орбит спутников по формулам ()

for j = 1:L % 0:L

t(j)=weeks+step*j; %-step;

%t1(j) = t(j)/60; %изменение дискретности текущего времени

%d_wn = (week - alm(i).Week);%если в альманахе учтены 1024

%d_wn = 0;

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;

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));

n=n0;

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 %k

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;

nuk =atan2(sqrt(1-alm(i).e^2)*sin(Ek),(cos(Ek)-alm(i).e));

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) = xkk *cos(OMEGAk)-ykk*cos(ik)*sin(OMEGAk);

Yk(j) = xkk*sin(OMEGAk)+ykk*cos(ik)*cos(OMEGAk);

 

Zk(j) = ykk*sin(ik);

%Дальности до спутников

PR(j) = sqrt((Xk(j) - Rx)^2 + (Yk(j) - Ry)^2 + (Zk(j) - Rz)^2);

 

%Перевод в географическую систему если требуется

%[lons,lats,hrs] = LLHECEF(Xk,Yk,Zk);

%(Llon(j),Llat(j),Hhr(j)) = [lons,lats,hrs];

 

xls = Xk(j) - Rx;

yls = Yk(j) - Ry;

zls = Zk(j) - 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(j) =satang;

 

%Расчет угла азимута

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) = atan2(xaz,yaz);

end

if az(j) < 0

az(j) = az(j) + pi*2;

end

end % j

for j = 1:L

[Llon(j),Llat(j),Hhr(j)] = LLHECEF(Xk(j),Yk(j),Zk(j));%преобразование координат

if j > 1

if abs(Llon(j)-Llon(j-1)) > pi

Llon(j) = Llon(j) + 2*pi;

end

end

end

 

j_color = j_color + 1;

if (j_color > 14)

j_color = 1;

end

%F_ont=get(gcf,'CurrentAxes'),'FontSize',16,'FontName','TimesNewRoman';%формат текста на графиках

S = color6(j_color);

%Графика

 

%График 1 для вывода графика убрать символы %{ и %} относящиеся к данному графику

%{

h_F1=gca;

plot3(Xk(:),Yk(:),Zk(:),S),

hold on,

set(get(gcf,'CurrentAxes'),'FontSize',12,'FontName','TimesNewRoman');

set(h_F1,'Position',[0.1 0.1 0.85 0.9]);

xlabel('Координата X')

ylabel('Координата Y'),

zlabel ('Координата Z'),grid on

%}

%{

%График 2 для вывода графика убрать символы %{ и %} относящиеся к данному графику

subplot(2,1,1),plot(t,az(:),S),

set(get(gcf,'CurrentAxes'),'FontSize',12,'FontName','TimesNewRoman')

hold on,

xlabel('Время'),

ylabel('Угол азимута,радиан')

grid on

 

subplot(2,1,2),plot(t,TT(:),S),

set(get(gcf,'CurrentAxes'),'FontSize',12,'FontName','TimesNewRoman')

hold on,xlabel ('Время'),

ylabel('Угол видимости'),

grid on

%}

%{

 

%График 3 для вывода графика убрать символы %{ и %} относящиеся к данному графику

subplot(2,1,1),plot(t(1:(j-1)),doppler(:),S),

set(get(gcf,'CurrentAxes'),'FontSize',12,'FontName','TimesNewRoman')

hold on,

xlabel('Время '),

ylabel('Доплеровская частота'), grid on

subplot(2,1,2), plot(Llon(1:j),Llat(1:j),S),

set(get(gcf,'CurrentAxes'),'FontSize',12,'FontName','TimesNewRoman')

hold on, xlabel('Долгота'),ylabel('Широта')

grid on

%}

%График 4 для вывода графика убрать символы %{ и %} относящиеся к данному графику

%hF=figure('Color','w','MenuBar','none')

subplot(1,2,1),plot(Xk(:),Yk(:),S),

 

set(get(gcf,'CurrentAxes'),'FontSize',12,'FontName','Times New Roman')

holdon,

xlabel('Проекция орбит на плоскость XY')

grid on

subplot(1,2,2), plot(t,PR(:),S),

set(get(gcf,'CurrentAxes'),'FontSize',12,'FontName','Times New Roman')

hold on,

xlabel('Время'),

ylabel('Дальность,метр '),grid on

end % i

clear

 

Файл Orbita_GPS_1.m

clear all

%Имя m-файла:Orbita_GPS_1.m

%Программа рассчитывает углы видимости, азимута и положение видимых спутников на заданный момент времени навигационных спутников GPS

%Входные данные:

%файл альманаха в формате 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;

%или

%OMEGAeDOT=0;

mu=398600500000000;

F_CONST = 4.442807633E-10;

kt=3;%установка времени на титульной надписи графика, определяется параметрами d2'; h; min; s и j или L;

%Задание цветов для графики

j_color = 0;

color6(1:14) = ['k' 'b' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h'];

%Входные данные

Dat = 'almanac_yuma_week0371_589824.txt';

d2='03/13/2016'; h=16.0; min=20.0; s=30.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)];

step=3600;%шаг отсчета времени в секундах (300=5 минутам);шаг можно изменять

L=(24*3600) / step;% убрать % перед L для вывода таблицы улов видимости и азимута

%L=1;% установить % перед L для вывода таблицы улов видимости и азимута

%Чтение альманаха

[alm,max_kol] = Yuma_GPS_Alm1(Dat);

nom = 1;

i = 0;

k = 0;

for i = 1: max_kol

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)

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;

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:L % 0:L

t(j)=weeks+step*(j); %-step;

%t1(j) = t(j)/60; %изменение дискретности текущего времени

%d_wn = (week - alm(i).Week);

%d_wn = 0;

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;

for k = 1: kol

i = nom_ns(k);

%Расчет орбит спутников по формулам

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);

Mk = alm(i).M0 + n*tk;

e=alm(i).e;

%решение уравнения Кеплера

eps = 1.0E-15;

y = e * sin(Mk);

x1 = Mk;

x = y;

fork = 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);

%(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;

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;

%Построение полярной системы координат

ifEL(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)

j_color = j_color + 1;

if (j_color > 14)

j_color = 1;

end

S = color6(j_color);

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;

%Формирование таблицы вывода времени 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');

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');

gridon;

%Построение изображений видимых спутников на круговой диаграмме

i=1;

for k=1:kol

i = nom_ns(k);

plot(xt(kt,i),yt(kt,i), 'Marker','d','MarkerSize',20)

title(t_itle);

xlabel(X_label,'FontSize',12,'FontName','TimesNewRoman')

set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman')

hold on

str1 = num2str(i, 2);

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

 

Файл ОRBITA_1.m

%Имя m-файла:ORBITA_1.m %Программа иллюстрирует процедуру размножения эфемерид и орбиты спутника ГЛОНАСС

%(демонстрация упрощенного варианта решения системы дифференциальных уравнений

%движения спутника)

%Программа выполняется совместно с функцией orbit_GL, использующей функцию MatLab ode45

%для решения дифференциальных уравнений методом Рунге-Кутта

%Входные данные:

%вектор координат x, y, z спутника XYZ (размерность-метр);

%вектор скоростей спутника по осям x, y, z (размерность-м/с) VXYZ;

%текущее время t= "начальное время": "шаг": "конечное время= "время в часах"*3600"

%Выходные данные:

%координаты спутника X, Y,Z (x, y, z) в абсолютной (относительной) системах координат;

%скорости спутника Vx, Vy, Vz в абсолютной системе координат;

%вектор текущего времени T;

%вектор текущих координат и скоростей V

%Расчет вектора входных параметров y

omega = 0.7292115*10^(-4);%- скорость вращения Земли

t=0:360:23*3600;

S=-omega*3*3600;% угол

%XYZ=[21840.10466;-9006.95351;-9696.59786];%координнаты спутника

XYZ=[9795803.22265;-7174949.70703;22480344.23828 ];%координнаты спутника

mS=[cos(S) -sin(S) 0;sin(S) cos(S) 0;0 0 1]; %матрица преобразования координат

 

%VXYZ=[-1.19933288;0.58113958;-3.25131421];%скорости спутника

VXYZ=[2773.857116;1295.602798;-814.5313262];

ys1=[mS*XYZ]';%вектор преобразованных координат

ys2=[mS*VXYZ]'+omega*[-ys1(2) ys1(1) 0];% вектор преобразованных скоростей

y=[ys1 ys2];%вектор начальных условий

%Расчет орбиты спутника с помощью функции ode45

%[T,V] = ode45(@orbit_GL,[0:360:23*3600],[y],[]);

[T,V] = ode45(@orbit_GL,[t],[y],[]);

% Координаты и график орбиты спутника

X=V(:,1);

Y=V(:,2);

Z=V(:,3);

 

subplot(2,1,1), plot3(X,Y,Z),grid on

set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman');

%Vx=V(:,4);

%Vy=V(:,5);

%Vz=V(:,6);

%subplot(1,3,2), plot3(Vx,Vy,Vz)

% Координаты и график орбиты спутника в системе координат ПЗ90

S=omega*T;

x= X.*cos(S)+Y.*sin(S);

y =-X.*sin(S)+Y.*cos(S);

z =Z;

subplot(2,1,2), plot3(x,y,z),grid on

set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman');

function [dy1 y1]= orbit_GL(t,y1)

%Имя функции: orbit_GL

%Функция записи системы дифференциальных уравнений для решения с помощью стандартной

% программы MatLab

dy1 = zeros(6,1);

prom=398600.44*10^9/((y1(1)*y1(1)+y1(2)*y1(2)+y1(3)*y1(3))^1.5);

dy1=[y1(4) y1(5) y1(6) [-y1(1) -y1(2) -y1(3)]*prom]';

 

 

Функции и файлы из папки ALM_CH4701_V3

Файл convers_АLM_GLN_YUMA_GPS.m

%Имя файла: 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\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-84ellipsoid parameters

A_PZ90_M = 6378136.0; %6378 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;

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

ifstr1=='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;

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;% ifstr1== '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;

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);

 

Функция Gln_data_from_NA

function [time_UTC]= Gln_data_from_NA(leap_year, day_from_leap);

%Имя функции:Gln_data_from_NA

%function [time_UTC]= Gln_data_from_NA(leap_year, day_from_leap);

% преобразует NA - день привязки альманаха

% от ближайшего високосного года- leap_year

% в текущую дату: timeUTC (год, месяц, день)

DnMon= [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]; %количество дней в месяцах

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;

 

Функция write_GPS_alm

 

function []=write_GPS_alm(fw,alm)

%Имя функции:write_GPS_alm

%Функция записывает альманахGPS в совместный альманах спутников GPS и ГЛОНАСС в формате YUMA

i=0;

format long e;

for i=1:31

ifalm(i).ID > 0

%Заголовок альманаха

fprintf(fw,'**** Week%i almaNAU for PRN-0%i **********\n',alm(i).Week, alm(i).ID);

%Номер спутникаGPS

fprintf(fw,'ID: %i\n',alm(i).ID);

%Здоровье спутникаGPS

fprintf(fw,'Health: %i\n',alm(i).Health);

%Эксцентриситет орбиты спутникаGPS

strdop = e_norm(alm(i).ecc, 9);

fprintf(fw,'Eccentricity: %s\n',strdop);

%Время от начала неделиGPS, на которое заданы параметры альманаха

fprintf(fw,'Time of Applicability(s): %6.4f\n',alm(i).TOA);

%Наклонение орбиты спутника GPS

fprintf(fw,'Orbital Incluation(rad): %0.10f\n', alm(i).deltai);

%Скорость изменения восходящего узла орбиты спутника GPS

strdop= e_norm(alm(i).OMEGADOT, 9);

fprintf(fw,'Rate ofRight Ascen(r/s): %s\n',strdop);

%Корень квадратный из большой полуоси орбиты спутник аGPS

fprintf(fw,'SQRT(A) (m^1/2): %4.7f\n',alm(i).A05);

%Долгота восходящего узла орбиты спутникаGPS

strdop= e_norm(alm(i).omega0, 9);

fprintf(fw,'Right Ascen at Week(rad): %s\n',strdop);

%Аргумент перигея орбиты спутникаGPS

strdop = e_norm(alm(i).omega, 9);

fprintf(fw,'Argument of Perigee(rad): %s\n',strdop);

%Средняя аномалия спутникаGPS

strdop= e_norm(alm(i).M0, 9);

fprintf(fw,'Mean Anom(rad): %s\n', strdop);

%Коэффициенты полинома для учета поправок времени

strdop= e_norm(alm(i).Af0m, 9);

fprintf(fw,'Af0(s): %s\n',strdop);

strdop= e_norm(alm(i).Af01, 9);

fprintf(fw,'Af1(s/s): %s\n', strdop);

%Номер недели

fprintf(fw,'week: %i \n',alm(i).Week);

%Запись альманахаGPS в файл

fprintf(fw,'\n');

end; %if alm(i).ID > 0

end;%i

 

Функция GLN_satfind_YUMA

 

function []= GLN_satfind_YUMA(a,timeUTC_leap, alm_gln, fid);

%Имя функции:GLN_satfind_YUMA

%функция обрабатывае входные данные

KOL_GLN= 24;

A_PZ90_KM = a / 1000;

for (i = 1:KOL_GLN)

%alm_gln.Na -(сек) время привязки альманаха от начала предшествующего високосного года

day_from_leap = alm_gln(i).Na;

timeUTC= Gln_data_from_NA(timeUTC_leap.year,day_from_leap);

timeUTC.ti = 0.0;

nut=0;

S0 = s0_Nut(timeUTC, nut);

time_s0 = S0.s0_m_mod; %time_s0 - истинное звездное время в текущий момент обсервации

year = timeUTC.year;

% leap_year = fix(year / 4) * 4; % ближайший к текущему (предыдущий) високосный год

leap_year = timeUTC_leap.year; % ближайший к текущему(предыдущий) високосный год

%------------

% 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();

ti= alm_gln (i).TLambdaN; % текущее время обсервации от начала дня

n00= fix(ti / 86400);

% n0 - номер текущих суток внутри4-х летнего периода(от ближайшего високосного года)

n0= JD_from_epohi(leap_year, timeUTC) + n00+ 1;

prn =alm_gln (i).ID;

health = alm_gln (i).Health;

Hn=alm_gln (i).Hn;

if ((prn >0) &(Hn > 0))

gln_a(A_PZ90_KM, i, n0, ti, time_s0, alm_gln,timeUTC_leap, fid);

end; % f ((prn > 0) & (health == 1))

end; % for (i = 1: KOL_GLN)

fclose(fid);

 

 

Функцияs0_Nut

function [S0] = s0_Nut(timeUTC, nut)

%Имя функции: s0_Nut

%функция рассчитывает истинное и среднее звездное время по формулам ()

% среднее звездное время s0 на 0ч UTC

%year = 1993; mon = 1; day = 0;

%fprintf('function s0_m - start \n');

jd2000 = 2451545; % 12h UTC 1 января

% Применить функцию JD_data

[jd, day_year] = JD_data(timeUTC);

if (jd == NaN)

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;

 

Функция JD_from_epohi

function [jd] =JD_from_epohi(epoha, timeUTC);

%Имя функции:JD_from_epohi

%Функция вычисляет 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;

 

Функция 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

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;

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 -...

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 ~