Порядок выполнения лабораторной работы
Создайте отдельную папку и скопируйте в нее папку «Орбита_GPS».
Запустите систему «Matlab» и настройте ее на доступ к Вашей папке.
Откройте в системе «Matlab» файл «Orbita_GPS1.m».
Выполните файл «Orbita_GPS1.m». Программа при первом запуске настроена таким образом, что результатом ее выполнения будет графическое окно (рис. 1).
 |
Рисунок 1.- Графическое окно файла «Orbita_GPS1.m» после первого
включения программы
Откройте текст (листинг m- файла «Orbita_GPS1.m»), внимательно изучите его и все комментарии, начинающиеся с символа «%».
Обратите внимание на то, что в файле «Orbita_GPS1.m», имеются следующие программные функции:
Yuma_GPS_Alm1(Dat)- функция чтения данных альманаха,
ECEFLLH(lon, lat,hr) – функция преобразования координат,
Tim(d2,h,min,s) – функция расчета начала отсчета времени.
Перед тем, как продолжить работать с программой «Orbita_GPS1.m», необходимо ознакомиться с данными программными функциями. Для этой цели служит m- файл «Example_Orbita.m».
Откройте файл «Example_Orbita.m» и выполните указания, приведенные в строке 2. В окне MatLab «Command Window» Вы увидите результат выполнения программы в виде цифр, являющихся данными альманаха, но записанными в другом порядке. Сопоставьте полученные данные с данными альманаха и убедитесь в их соответствии.
Выполните указания, приведенные в строке 10. В окне MatLab «Command Window» Вы увидите результат выполнения программы в виде преобразованных координат.
Выполните указания, приведенные в строке 13. В окне MatLab «Command Window» Вы увидите результат выполнения программы в виде данных начала отсчета изучаемого орбитального движения навигационных спутников.
Подставте в файл «Example_Orbita.m» исходные данные с которыми Вы предполагаете работать и повторите все процедуры.
Вернитесь к файлу «Orbita_GPS1.m»), подставьте в него Ваши исходные данные и приступайте к исследованию орбитального движения спутников, последовательно выполняя указания в комментариях к программе.
Примеры исследования
Пример 1. Орбиты 27 спутников за 12 часов
Пример 2. Орбита спутника 19 за 23 часа
Пример 3. Доплеровская частота спутника 6 за 23 часа
Пример 4. Азимут спутника 23 за 2 недели
 |
Пример 5. Азимут спутника 30 за 2 недели
 |
Пример 6. Видимость 27 спутников за 24 часа
 |
Пример 7. Следы 6 спутников за 12 часов
 |
Пример 8 Проекции орбит 27 спутников на плоскость XY за 12 часов
 |
Пример 9. Дальность до спутника 10.
 |
Задания для самостоятельной подготовки и вопросы для самопроверки
5. 1. Ознакомтесь с разделом 20 ICD-GPS-200C.
5. 2. Изучите алгоритм расчета орбиты спутника и размерность составляющих, входящих в формулы (g. 1- g. 16).
5. 3. Подготовтесь к работе с программным пакетом MatLab.
5. 4. Согласуйте данные таблиц 2 и 3.
5. 5. Рассчитайте количество секунд в GPS неделе.
5. 5. Рассчитайте на дату дня Вашего рождения week – номер GPS недели, modeweek- модифицированную GPS неделю ,d – номер дня в году, dweek – номно дня недели,weeks – количество секунд в неделе в день Вашего рождения.
5. 6. Преобразуйте координаты любого города из системы «широта, долгота, высота» в систему X, Y, Z.
5. 7. Найдите в справке системы MatLab функции: floor, fix, DAYSDIF, fopen, feof, fscanf, length, clear all, color, sqrt, sin, cos, acos, subplot, plot, hold on,xlabel, ylabel, grid on; операторы: for, if, end, while, else, break, elseif . Дайте им объяснение.
5. 8. Объясните характер изменения доплеровской частота от времени.
5. 9. Дайте качественное описание связи между доплеровской частотой и дальностью до навигационного спутника.
5. 10. Ответьте на вопрос: на скольких орбитах расположены спутники GPS.
5. 11. Сколько спутников располагаются на каждой из орбит для альманаха, который Вы используете при исследовании орбитального движения спутников.
5. 12. По какому признаку альманаха можно определить принадлежность спутника к одной или другой орбите.
Отчет
Отчет по лабораторной работе должен содержать:
3. 1. Пояснения к программам MatLab.
3. 2. Графическое и цифровое отображение результатов выполнения программы.
3. 3. Выводы по результатам моделирования орбиты навигационного спутника.
Результаты работы и отчет могут предоставляться в электронном виде.
Приложения
7.1 Приложение 1. Листинг программы «Orbita_GPS1.m».
clear all
%Входные данные
%Для работы с данной программой необходимо иметь текстовый
%файл альманаха в формате Yuma и входным данным Dat присвоить
% имя этого файла
Dat = '02_11_15_30.alm';
%Dat = 'имя файла альманаха, с которым Вы быдете работать';
%Данные о начале отсчета
d2='11/02/2004';h=10.0;min=0.0;s=0.0;
%В строку 11 или 12 введите дату и время,с которым Вы быдете работать
%d2='12/30/2004';h=12.0;min=0.0;s=0.0;
%d2='03/21/2003';h=9.0;min=49.0;s=40.0;
%Координаты позиции приемника:lat,lon -(радианы),hr-метры
lat=0.881278698506528;
%lat= pi/4;
lon=0.53169758803674;
hr=122.899802776054;
%В строки 18, 19, 20 введите координаты приемника,с которым Вы быдете работать
%lat=0.881278698506528;
%lon=0.53169758803674;
%hr=122.899802776054;
%Постоянные
%Скорость вращения Земли
OMEGAeDOT=7.2921151467e-005;
%OMEGAeDOT=0;
mu=398600500000000;
%Шаг с каким будут рассчитываться параметры орбит (step,секунды)
%step=60;
%Количество точек (L), в которых будут рассчитываться параметры орбит
%L=12*3600/step,L читается так: количество часов (например,12)
%*число секунд в часе (3600) деленное на шаг (step)
%L=12*3600/step;
step=60;
L=24*1*3600/step;
%В строки 35, 36 введите step и L, с которыми Вы будете работать
%step = 3600;
%L=24*3600/step;
%Задание цветов для графики
j_color = 0;
color6(1:14) = ['k' 'b' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h'];
%Чтение альманаха
[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=3504451.023;Ry=2061316.876;Rz=4897990.975;
%Rx=0;Ry=0;,Rz=0;
%Выбор спутников
kol =1;
%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10];
%В строку 63 введите количество спутников, которые Вы хотите исследовать
%kol =9
%В одну из строк 65-71 введите номнра спутников,которые Вы хотите исследовать
%Количество номеров спутников должно совпадать с kol из строки 63
%nom_ns(1:kol) = [1 13 14 26 29]; %1
%nom_ns(1:kol) = [2 5 22 28 30]; %2
%nom_ns(1:kol) = [3 6 7 31]; %3
%nom_ns(1:kol) = [4 11 15 17 24 ]; %4
%nom_ns(1:kol) = [8 9 25 27]; %5
%nom_ns(1:kol) = [10 18 20 21 23]; %6
%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10 11 13 14 15 17 18 19 20 21 23 24 25 26 27 28 29 30 31] ;
nom_ns(1:kol)=[10];
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 = abs(week - alm(i).Week);
tk = t(j) + d_wn * 604800 - alm(i).TOA;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
end % while
%%%%%%%%%%%%%%%%
%Справочник
%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 = C(10,i)+n*tk;
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;
%--------------------
F_CONST = 4.442807633E-10;
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];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DETERMINE IF A CLEAR LINE OF SIGHT EXISTS
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(i) = ralt
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;
% El(i,j)=satang;
TT(j) =satang;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Вариант 2 (азимут)
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);
%az(j) = atan2(yaz,xaz);
%az(j) = atan(yaz/xaz);
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
%axis([349258,352798, 2,15]);
%axis([-2.6e7 2.6e7 -2.6e7 2.6e7 -2.6e7 2.6e7]);
% color6(1:14) = ['k' 'b' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h'];
j_color = j_color + 1;
if (j_color > 14 )
j_color = 1;
end
S = color6(j_color);
%Вывод всех графиков на одну панель командой subplot
%График 1
subplot(4,2,1),plot(t,az(:),S),hold on,xlabel('Время'),ylabel('Угол азимута,радиан')
grid on
%График 2
subplot(4,2,2),plot(t(1:(j-1)),doppler(:),S),hold on,xlabel('Время'),
ylabel('Доплеровская частота'),grid on
%График 3
subplot(4,2,3),polar(az(1:j),t,S),hold on, grid on
%График 4
subplot(4,2,4),plot(t,TT(:),S),hold on,xlabel('Время'), ylabel('Угол видимости'),grid on
%График 5
subplot(4,2,5),plot3( Xk(:),Yk(:),Zk(:),S),hold on, xlabel('Координата X')
ylabel('Координата Y'),zlabel('Координата Z'),grid on
%График 6
subplot(4,2,6),plot(Llon(1:j),Llat(1:j),S),hold on, xlabel('Долгота'),ylabel('Широта')
grid on
%График 7
subplot(4,2,7),plot(Xk(:),Yk(:),S),hold on,grid on
%График 8
subplot(4,2,8),plot(t,PR(:),S),hold on, xlabel('Время'),ylabel('Дальность,метр '),grid on
hold on
% subplot(3,2,5)
%hold on
%subplot(3,2,6),grid on
% hold on
%1.График зависимости азимута спутников от времени
%для вывода графика убрать символ %
% figure
%plot(t,az(:),t,TT(:),S),hold on,xlabel('Время'),ylabel('Угол азимута,радиан'),grid on
% figure
%2.График зависимости доплеровской частоты от времени
%для вывода графика убрать символ %
%plot(t(1:(j-1)),doppler(:),S),hold on,xlabel('Время'),ylabel('Доплеровская частота')
%grid on
%figure
%3.График зависимости азимута спутников от времени в полярной
%системе координат для вывода графика убрать символ %
%figure
%polar(az(1:j),t,S),hold on, grid on
%figure
%polar(TT(1:j),t,S),hold on, grid on
%4. График зависимости угла видимости спутника от времени
%для вывода графика убрать символ %
%plot(t,TT(:),S),hold on,xlabel('Время'), ylabel('Угол видимости'),grid on
%figure
%5. График орбит спутников для вывода графика убрать символ %
%plot3( Xk(:),Yk(:),Zk(:),S),hold on, xlabel('Координата X')
%ylabel('Координата Y'),zlabel('Координата Z'),grid on
%figure
%6. График следа спутника в координатах широта долгота
%для вывода графика убрать символ %
%plot(Llon(1:j),Llat(1:j),S),hold on, xlabel('Долгота'),ylabel('Широта')
%grid on
%figure;
%7.График проекции орбиты спутника на плоскость XY
%для вывода графика убрать символ %
%plot(Xk(:),Yk(:),S),,hold on,grid on
%figure;
%8.График зависимости дальности до спутника от времени
%для вывода графика убрать символ %
%plot(t,PR(:),S),hold on, xlabel('Время'),ylabel('Дальность,метр '),grid on
%figure
%hold on
%set(hplot, 'LineWidth',1.5);
%axis off;
%end
%get(hplot,'LineWidth')
%end
[xn,yn,zn,xe,ye]';
%[xn1,yn1,zn1,xe1,ye1]'
end % i
%Данные для котроля
Out= [week,alm(31).Week,max_kol];
clear
7.2 Приложение 2. Листинг программы «Yuma_GPS_Alm1.m».
function [alm,max_kol] = Yuma_GPS_Alm1(Dat)
fid =fopen(Dat,'rt');
max_kol=28;
max_kol = 1;
while not(feof(fid))
max_kol = max_kol + 1;
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);
end
%C=A';
%Week_Alm(ID) = alm(ID).Week;
end
7.3 Приложение 3. Листинг программы «Tim.m».
function [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s)
%d2='10/35/2003';h=23.0;min=59.0;s=59.0;
%week = ceil(DAYSDIF('1/6/1980',d2,3)/7);% текущая неделя
week = floor(DAYSDIF('1/6/1980',d2,3)/7);% предыдущая неделя
modeweek=week-1024;
d = DAYSDIF('1/6/1980',d2,3);
%dweek=6+fix(d-week*7);
dweek=fix(d-week*7);
weeks=(dweek)*24*60*60+h*60*60+min*60+s;
7.4 Приложение 4. Листинг программы «ECEFLLH.m».
function [Rx,Ry,Rz] = ECEFLLH(lon, lat,hr)
%Преобразование координат
%Входные данные:lon-долгота,lat-широта,h-высота;a,b-большая
%и малая оси эллипсоида
%Выходные данные:Rx,Ry,Rz- координаты в ECEF
%Для WGS-84;
%Rx=3504451.023;Ry=2061316.876;Rz=4897990.975;
%lat=0.881278698506528;lon=0.53169758803674;
%hr=122.899802776054;
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);
7.5 Приложение 5. Листинг программы «LLHECEF.m».
function [lons,lats,hrs] = LLHECEF(Xk,Yk,Zk)
a=6378137.0;
b=6356752.314;
%X=3504451.023;Y=2061316.876;Z=4897990.975;
%lat=0.881278698506528;lon=0.53169758803674;
%hr=122.899802776054;
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;
%[lat,lon,hr]'
end
7.6 Приложение 6. Листинг программы «Example_Orbita.m».
clear all
%Закомментируйте все строки, кроме 3 и 9 и исполните программу
%Dat = '02_11_15_30.alm';
%[alm,max_kol] = Yuma_GPS_Alm1(Dat);
%[alm(1:31).ID, alm(1:31).Health,alm(1:31).e, alm(1:31).TOA,...
%alm(1:31).deltai, alm(1:31).OMEGADOT,alm(1:31).A05, alm(1:31).omega0,...
%alm(1:31).omega, alm(1:31).OMEGADOT,alm(1:31).M0, alm(1:31).Af0,...
%alm(1:31).Af1, alm(1:31).Week]'
%max_kol
%Закомментируйте все строки, кроме 11 и 12 и исполните программу
%lat=0.881278698506528;lon=0.53169758803674;hr=122.899802776054;
%[Rx,Ry,Rz] = ECEFLLH(lon, lat,hr)
%Закомментируйте все строки, кроме 14 и 15 и исполните программу
d2='11/02/2004';h=10.0;min=0.0;s=0.0;
[week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s)
Рекомендуемая литература
1. ICD-GPS-200C (есть электронная копия на компьютере), 160 с.
2. Дьяконов В. П., Абраменкова И. В., Круглов В. В. Matlab 5 с пакетами расширений.- М.:Нолидж. – 2001.- 878 с.
3. Бабак В. П., Конін В. В., Харченко В. П. Супутникова радіонавигація. К. : Техніка.-2004.- 327 с.