Телескопы покупают здесь


A A A A Автор Тема: Программа для перерасчета координат спутника ГЛОНАСС  (Прочитано 423 раз)

0 Пользователей и 1 Гость просматривают эту тему.

Оффлайн EVolАвтор темы

  • Новичок
  • *
  • Сообщений: 1
  • Благодарностей: 0
  • Мне нравится этот форум!
Доброго времени суток. Нуждаюсь в вашей помощи. Выдали задание - реализовать программу перерасчета координат спутника ГЛОНАСС в Matlab. Реализовывал программу по интерфейсному контрольному документу ГЛОНАСС 2016 г. страницы 57-62 ссылка - http://russianspacesystems.ru/wp-content/uploads/2016/08/IKD.-Obshh.-opis.-Red.-1.0-2016.pdf . Дошел до последнего пункта, до самого интегрирования, прикладываю скрин формул. Вопрос заключается в том, что я не могу понять как тут происходит интегрирование, если в правой части получается численное значение, я как понял нужно интегрировать по времени, в которое нужно рассчитать координаты. Можете подсказать/помочь как это реализуется. Может я что то неправильно сделал или неправильно понял. Также прикладываю код программы в Matlab.

%начальные условия
tb=11700;
ti=12300;
S0=7500;
JD0=2456177.625;
xtb=7003.008789;
ytb=-12206.626953;
ztb=21280.765625;
Vxtb=0.7835417;
Vytb=-2.8042530;
Vztb=1.3525150;


%расчет начальных улсовия
wz=7.2921151467*10^(-5);
Stb=S0+wz*(tb-10800);
x0tb=xtb*cos(Stb)-ytb*sin(Stb);
y0tb=xtb*sin(Stb)+ytb*cos(Stb);
z0tb=ztb;
xt0tb=Vxtb*cos(Stb)-Vytb*sin(Stb)-wz*y0tb;
yt0tb=Vxtb*sin(Stb)+Vytb*cos(Stb)+wz*x0tb;
zt0tb=Vztb;

%направляющие косинусы и радиус-вектор на луну
T=(JD0+(tb-10800)/86400-2451545)/36525;
e=0.4090926006-0.0002270711*T;
G=1.4547885346 + 71.0176852437*T-0.0001801481*(T^2);
sigmal=2.1824391966-33.7570459536*T + 0.0000362262*(T^2);
ql=2.3555557435 + 8328.6914257190*T + 0.0001545547*(T^2);
al=3.84385243e+05;% большая полуось орбиты Луны
el=0.054900489;% эксцентриситет лунной орбиты
il=0.0898041080;%среднее наклонение орбиты Луны к плоскости эклиптики
Sezv=1-(cos(sigmal))^2 *(1-cos(il));
Etazv=sin(sigmal)*sin(il);
Jzv=cos(sigmal)*sin(il);
Se11=sin(sigmal)*cos(sigmal)*(1-cos(il));
Se12=1-(sin(sigmal))^2 * (1-cos(il));
Eta11=Sezv*cos(e)-Jzv*sin(e);
Eta12=Se11*cos(e)+Etazv*sin(e);
J11=Sezv*sin(e)+Jzv*cos(e);
J12=Se11*sin(e)-Etazv*cos(e);

%Решение уравнения кеплера решается методом итераций, пока Е Е (пред) л л ? не будет меньше 10–8,
eps = 10^(-8); % точность
x0 = 0.00; % начальное приближение
N = 100; % количество итераций, чтобы не было зацикливаний
g = inline('ql+el*sin(x0)','x0','ql','el');
x1 = g(x0,ql,el); % первое значение
for i = 1 : N % делаем максимум 100 итераций
    if abs(x1 - x0) <= eps
        break
    end
    x0 = x1;
    x1 = g(x0,ql,el);
end
El=x1;

sinVl=(sqrt(1-el^2)*sin(El))/(1-el*cos(El));
cosVl=(cos(El)-el)/(1-el*cos(El));

Sel=(sinVl*cos(G)+cosVl*sin(G))*Se11+(cosVl*cos(G)-sinVl*sin(G))*Se12;
Etal=(sinVl*cos(G)+cosVl*sin(G))*Eta11+(cosVl*cos(G)-sinVl*sin(G))*Eta12;
Jl=(sinVl*cos(G)+cosVl*sin(G))*J11+(cosVl*cos(G)-sinVl*sin(G))*J12;
rl=al*(1-el*cos(El));

%направляющие косинусы и радиус-вектор на солнце
ws = -7.6281824375 + 0.0300101976*T + 0.0000079741*(T^2);
qs=6.2400601269+628.3019551714*T- 0.0000026820*(T^2);
es=0.016719;
as=1.49598*10^8;

%Решение уравнения кеплера решается методом итераций, пока Е Е (пред) л л ? не будет меньше 10–8,
eps = 10^(-8); % точность
y0 = 0.00; % начальное приближение
N = 100; % количество итераций, чтобы не было зацикливаний
g = inline('qs+es*sin(y0)','y0','qs','es');
y1 = g(y0,qs,es); % первое значение
for i = 1 : N % делаем максимум 100 итераций
    if abs(y1 - y0) <= eps
        break
    end
    y0 = y1;
    y1 = g(y0,qs,es);
end
Es=y1;

sinVs=(sqrt(1-es^2)*sin(Es))/(1-es*cos(Es));
cosVs=(cos(Es)-es)/(1-es*cos(Es));

Ses=cosVs*cos(ws)-sinVs*sin(ws);
Etas=(sinVs*cos(ws)+cosVs*sin(ws))*cos(e);
Js=(sinVs*cos(ws)+cosVs*sin(ws))*sin(e);;
rs=as*(1-es*cos(Es));

%Ускорения от лунных и солнечных гравитационных возмущений
Gl=4902.799*10^(9);
Gs=13271244*10^(13);

xkl=x0tb/rl;
ykl=y0tb/rl;
zkl=z0tb/rl;
Gkl=Gl/(rl^2);
deltal=((Sel-xkl)^2+(Etal-ykl)^2+(Jl-zkl)^2)^(3/2);
jx0l=Gkl*(((Sel-xkl)/deltal)-Sel);
jy0l=Gkl*(((Etal-ykl)/deltal)-Etal);
jz0l=Gkl*(((Jl-zkl)/deltal)-Jl);

xks=x0tb/rs;
yks=y0tb/rs;
zks=z0tb/rs;
Gks=Gs/(rs^2);
deltas=((Ses-xks)^2+(Etas-yks)^2+(Js-zks)^2)^(3/2);
jx0s=Gks*(((Ses-xks)/deltas)-Ses);
jy0s=Gks*(((Etas-yks)/deltas)-Etas);
jz0s=Gks*(((Js-zks)/deltas)-Js);

J02=1082625.75*10^(-9);
GM=398600441.8*10^6;
ae=6378136;

r0=sqrt(x0tb^2+y0tb^2+z0tb^2);
GMk=GM/(r0)^2;
xk0=x0tb/r0;
yk0=y0tb/r0;
zk0=z0tb/r0;
ro=ae/r0;

fun1=-GMk*xk0 - (3/2)*J02*GMk*xk0*ro^2 *(1-5*zk0^2)+jx0s+jx0l;
fun2=-GMk*yk0 - (3/2)*J02*GMk*yk0*ro^2 *(1-5*zk0^2)+jy0s+jy0l;
fun3=-GMk*zk0 - (3/2)*J02*GMk*zk0*ro^2 *(3-5*zk0^2)+jz0s+jz0l;