Орбитальное движение навигационных спутников

Орбитальное движение навигационных спутников

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

Орбитальное движение спутников GPS

Цель работы: изучение основ вторичной обработки информации в навигационных приемниках, работающих со спутниками GPS.

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

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

Для изучения и исследования орбитального движения спутников GPS применяются следующие m-файлы и функции: ECEFLLH.m, LLHECEF.m - функции для преобразования координат; Tim.m - функция для расчета времени; Yuma_GPS_Alm1.m - функция считывания данных альманаха и m-файл Orbita_GPS.m, рассчитывающий орбиты спутников GPS по данным альманаха, представленного в формате YUMA. Более полные сведения о перечисленных программах даются в комментариях к ним.

Порядок выполнение работы.

Пакет программ для выполнения лабораторной работы расположен в папке ORBITA_GPSv1. Для выполнения работы в качестве входных данных потребуется альманах спутников GPS в формате YUMA, который можно получить с нескольких сайтов или с навигационного приемника, если приемник имеется в лаборатории.

1. Запишите альманах в формате YUMA на жесткий диск компьютера. По умолчанию в комплексе программ установлен альманах с именем almanac_yuma_week0371_589824.txt (альманах приведён в приложении). Имя альманаха внесите в отчет.

2. Создайте папку ORBITA_GPSv1_My и скопируйте в нее все программы из папки ORBITA_GPSv1 и альманах, с которым Вы будете работать.

3. Запустите MatLab.

4. Обратитесь к папке ORBITA_GPSv1_My и откройте ее.

5. Последовательно откройте функции и m-файлы из папки ORBITA_GPSv1: Orbita_GPS.m,Yuma_GPS_Alm1.m, Tim.m, LLHECEF.m, ECEFLLH.m,PR_Tim.m,PR_Yuma_GPS_Alm1.m, almanac_yuma_week0371_589824.txt. Изучите тексты программ и комментарии в каждой программе.

6. Выполните m-файл PR_Tim.m и изучите в командном окне полученный результат. Сформируйте на основе m-файла PR_Tim.m файл с именем PR_Tim_My.m. Введите в сформированный файл данные, соответствующие дате выполнения работы, сохраните файл и выполните файл. Результат из командного окна запишите в отчет.

Пример выполнения m-файла PR_Tim.m с данными: 10.03.2016, 15:20

ans =

Columns 1 through 2

1887 863

Columns 3 through 4

13213 4

Column 5

7. Выполните m-файлPR_Yuma_GPS_Alm1.m и изучите в командном окне полученный результат. Сопоставьте данные из командного окна с данными альманаха almanac_yuma_week0371_589824.txt. Сформируйте на основе m-файлаPR_Yuma_GPS_Alm1.m файл с именем PR_Yuma_GPS_Alm1_My.m. Введите в сформированный файл имя записанного альманаха, сохраните файл и выполните файл, сравните результат из командного окна с данными альманаха и результат запишите в отчет.

8. Очистите командное окно из папки ORBITA_GPSv1_My, откройте m-файл с именем Orbita_GPS.m, выполните m-файл, ознакомьтесь с графическим результатом.

Рис. 1 Результат выполнения Orbita_GPS.m

9. Присвойте переменной Dat в файле Orbita_GPS.m имя альманаха из п.1 и введите данные о начале отсчёта как указано в комментариях к файлу. Установитескорость вращения Земли равную нулю. Убедитесь, что все графики, кроме графика 1 закомментированы. Исполните файл. Результат выполнения график 1 запишите на диск, как приложение к отчёту. В отчет внесите описание графика 1.

Рис. 2 Результат выполнения Orbita_GPS.m с данными: 10.03.2016, 15:20

10. Установите скорость вращения Земли, равную 7.2921151467e-005. Выберите номер спутника, с которым будете работать.

Выберите номер спутника равный дню вашего рождения (число спутников GPS, как правило, близко к30). Внесите этот номер в отчет.

Рис. 3 Результат выполнения графика 1 с номером спутника 13

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

Рис.4 Результат выполнения графика 2 с номером спутника 13

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

Рис.5 Результат выполнения графика 3 с номером спутника 13

12. Разблокируй график 4, исполните файл. Объясните поведение проекции орбиты на плоскость XY и изменение дальности до спутника от времени. График запишите на сменный диск, а интерпретацию результатов в отчет.

Рис.6 Результат выполнения графика 4 с номером спутника 13

13. Выберите группу спутников, принадлежащих одной из 6 орбит, и повторите расчётные процедуры пунктов 9-12. График запишите на сменный диск, а интерпретацию результатов в отчет.

Рис.7 Пример шести орбит спутников GPS, при скорости вращения Земли равной нулю

Вопросы:

1. Объясните зависимость изменения доплеровской частоты при движении спутника по орбите.

2. Объясните зависимость изменения дальности до спутника от времени для неподвижного наблюдателя.

3. Для каких целей используются данные альманаха в спутниковых навигационных приёмниках.

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

 

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

Наблюдение спутников GPS

Цель работы: изучение основ вторичной обработки информации в навигационных приёмниках, работающих со спутниками GPS.

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

1. Запишите альманах в формате YUMA на жесткий диск компьютера. По умолчанию в комплексе программ установлен альманах с именем almanac_yuma_week0371_589824.txt (альманах приведён в приложении). Имя альманаха внесите в отчет.

2. Создайте папку ORBITA_GPSv2_My.

3. В папку ORBITA_GPSv2_My запишите: файл альманаха, с которым будете работать, скопируйте функции ECEFLLH, LLHECEF, Yuma_GPS_Alm1, файл Tim.m из папки ORBITA_GPSv1 и файл Orbita_GPS_1.m из папки ORBITA_GPSv2. Эти данные внесите в отчет в соответствующий раздел.

4. Пользуясь комментариями к файлу Orbita_GPS_1.m установите время, в которое должны наблюдаться спутники. Выполните файл и результаты, полученные на графике, занесите в отчет.

Задание: день, в который должны наблюдаться спутники, будет равен дате вашего рождения (пример 13 декабря 1994 будет – 12/13/1994), точное время задайте произвольно.

5. Во входных данных файла Orbita_GPS_1.m измените время на 6 часов, выполнить файл и полученную конфигурацию наблюдаемых спутников занести в отчет. Объяснить разницу между результатами п.4-5.

Рис.1 Полученная конфигурация видимости29 спутников 13.03.16 с разницей в 6 часов

6. Во входных данных файла Orbita_GPS_1.m, пользуясь комментариями установить вывод наблюдаемых спутников в таблицу, появляющуюся в командном окне после исполнения файла и сопоставить эту таблицу с данными, полученными с навигационного приемника на установленное в файле время.

Примечание: практически все навигационные приёмники имеют опцию для отображения на компьютере наблюдаемых спутников в полярной системе координат.

Вопросы:

1. В каких пределах изменяется угол видимости спутника?

2. В каких пределах изменяется азимут спутника?

3. Сколько спутников GPS могут одновременно находиться в зоне видимости?

 

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

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

Одним из отличий системы ГЛОНАСС от GPS является то, что размножение эфемерид в системе ГЛОНАСС производится через решение системы дифференциальных уравнений.

Приведённая программа иллюстрирует применение функции MatLab при решении системы дифференциальных уравнений методом Рунге-Кутта для размножения эфемерид спутников ГЛОНАСС.

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

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

1. Создайте папку ORBITA_GLONASS_My и скопируйте в неё из папки ORBITA_GLONASS m-файл ОRBITA_1.m и функцию orbit_GL.

2. Изучите программные процедуры и комментарии к m-файлу ОRBITA_1.m и функции orbit_GL.

3. С сайта Российского космического агентства или навигационного приемника, работающего с системой ГЛОНАСС, запишите данные координат и скорости работоспособного спутника и включите их в отчет по лабораторной работе.

 

Рис.1 Данные координат и скорости работоспособных спутников

4. Откройте файл ОRBITA_1.m из папки ORBITA_GLONASS_My и выполните его. Ознакомьтесь с полученным графическим изображением орбиты спутника ГЛОНАСС.

 

Рис.1 Примеры выполнения файла ОRBITA_1.m

 

5. Введите координаты и скорости спутника ГЛОНАСС из п. 3 в строки входных данных файла ОRBITA_1.m, исполните файл. Проанализируйте полученное графическое изображение орбиты спутника ГЛОНАСС. Результаты анализа внесите в отчет.

6. Дополните файл ОRBITA_1.m процедурой вывода координат и скорости спутника на одно из значений текущего времени. Исполните файл и запишите полученные значения координат и скоростей спутника в отчет.

7. Измените в п. 6 текущее время на 900 секунд (граница интервала размножения эфемерид), исполните файл, занесите результат выполнения файла в отчет. Проанализируйте и объясните разницу в результатах, полученных в п.6 и 7.

 

Вопросы:

1. Для каких целей требуется размножать координаты и скорости навигационных спутников?

2. Какие параметры являются начальными условиями при решении системы дифференциальных уравнении орбитального движения спутников ГЛОНАСС?

3. В какой системе координат передаются данные о координатах и скорости в спутниковой системе ГЛОНАСС?

4. В какой системе координат решаются дифференциальные уравнения орбитального движения спутников ГЛОНАСС?

 

Преобразование данных навигационных спутников

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

Преобразование данных альманаха приемника СН4701 в формат YUMA

Цель: работы преобразование(конвертирование) данных альманаха навигационных спутников в формат YUMA.

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

В спутниковых радионавигационных системах ГЛОНАСС и GPS эфемериды рассчитываются по принципиально разным формулам. В то же время весьма целесообразным является применение в задачах прогнозирования навигационных сеансов использовать одинаковый формат данных альманаха. Отсюда вытекает задача преобразования (конвертирования) альманаха системы ГЛОНАСС в широко распространенный формат YUMA без потери данных об альманахе спутников ГЛОНАСС.

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

1. Запишите альманах, с приемника CH4701 в папку In_dat. Если не имеется возможности получить данные альманаха с приемника, то можно воспользоваться альманахом с именем 002.txt из приложения.

2. Создайте папку ALM_CH4701_V3_My

3. Создайте файлы и функци convers_ALM_GLN_YUMA_GPS.m, read_Alm, Gln_data_from_NA, s0_Nut, gln_a, read_Alm, write_GPS_alm, GLN_satfind_YUMA, Gln_data_from_NA,JD_from_epohi, week_GLONAS_gps, semi_axis, JD_epohi, JD_data в папке ALM_CH4701_V3_My и присвойте переменной Dat имя альманаха из папки In_dat. Выполните файл. Результат выполнения запишется в файл AlmGG.yum. Это и будет преобразованный в формат YUMA альманах спутников GPS и ГЛОНАСС.

4. Загрузите альманах AlmGG.yum в программу прогнозирования PLANNING и убедитесь в его работоспособности.

 

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

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

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

Программа рассчитывает углы азимута и места навигационных спутников 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

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

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

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

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

%Полуоси земно