[과제] COEs를 통해 Perifocal, ECI 좌표계에서 궤도 그리기

Kepler's Equation과 Newton Raphson방정식에 대해서 이해하면 M,E의 관계식을 해결할 수 있고 

그럼 궤도요소와 t를 받았을 때 M을 구하고, M과 뉴턴랩슨 방식을 이용해 E를 구한 뒤 v까지 구할 수 있다. 

 

이 개념을 코드에 적용시키면 시간에 따른 v, true anomaly를 구할 수 있다. 

그래서 지난 과제 마무리를 위한 포스팅을 작성해보려고 한다. 지난번엔 개념을 위주로 다루었다면 이번엔 이 개념을 코드에 적용시키는 방법을 다루어보겠다. 

 

[과제] Kepler's Equation과 Newton Raphson

 

과제 목표 : Molniya 궤도(이심률이 0.7로 매우큼)의 궤도 요소를 활용하여 2D, 3D좌표계 그리기 

Molniya Orbit

몰니야궤도란 고위도 지역(러시아)에서 통신 및 원격 탐사를 위한 커버리지를 제공하도록 설계된 궤도이다. 장축 반지름은 약 26,600km, 이심률은 0.737으로 원지점 고도가 40,000km에 달한다. 장축 반지름 통해 주기를 계산해보면 주기는 약 12시간이다. 기울기가 63.4(degree)인데, 이는 J2섭동을 긍정적으로 이용하여 argument of perigee를 고정시켜 특정 위치만 지나도록 설계된 아주 똑똑한 궤도임 ! 

 

출처 : Sellers understanding space

 

현재 자세히 설명은 못하나 이해한 바까지 작성해보자면 ,

 

J2섭동은 지구가 완전한 구가 아니라 적도 방향으로 불룩한 편평체이기 때문에 이 모양으로 인해 발생하는 외력이 발생하는 것이다. 그래서 대부분 인공위성이 그리는 궤도면이 흔들리는 세차 운동을 한다. 그런데 이 J2를 역이용하여 w(omega, argument of perigee)를 나타낼 수 있는데, w의 변화율이 0이 되는 inclination이 바로 63.4도이다 ! 그래서 기울기가 63.4일 때 w가 고정되므로 이를 이용하여 argument of perigee를 고정하여 항상 같은 곳을 바라보게 설계된 것임 

 

몰니야 궤도 시뮬레이션 

그래서 이 몰니야 궤도 요소를 가정하고 시뮬레이션을 만들어보았다. 

 

변수 정의

구분	변수	설명
INPUT	a	Semimajor axis = 26600km
	e	Eccentricity = 0.737
	i	Inclination = 63.4
	Omega	RAAN = 0 (임의 지정)
	omega	Argument of perigee(270)
	t	궤도가 출발한 뒤의 시간
OUTPUT		인공위성의 위치 벡터
		몰니야 궤도 애니메이션 – Perifocal, ECI좌표

 

알고리즘 

 

1. 초기 true anomaly와 지난 시간 hrs 입력 후 위상 구하기

* Newton-Raphson사용하여 M → E 구하기  

 

2.  t에 따라 변하는 궤도 표현을 위해 t 정의 및 t에 따른 M,E 계산

3. Perifocal 좌표계에서 궤도 표현

 

 

4. ECI 좌표계에서 궤도 표현

PQW좌표계와 IJK좌표계 변환하는 행렬 구하기

 

 

5. 2번에서 구한 R값에 A’를 곱해 IJK에서의 R,V로 표현

 

실행 코드

 

코드는 t(시간)에 따른 v를 구하고 v를 통해 2d(perifocal), 3d(eci)로 변환하는 과정만 나타내겠다. 

%PERIFOCAL FRAME ::  draw the orbit

T = 4*pi/n;
t = linspace(0,T,500);

% Time to Mean anomaly
M = Ini.M + n*(t-t0);
M = mod(M,2*pi);

% Mean anomaly to Eccentric anomaly
E = zeros(size(M));
for k = 1:length(M)
    E(k) = get_E(e,M(k));
end

% Eccenctic to True
%nu = true anomaly
nu = 2*atan2(sqrt(1+e).*sin(E/2), sqrt(1-e).*cos(E/2));
r = a*(1-e^2)./(1+e*cos(nu));

%position vector(2D) : [r*cosv r*sinv]

x_pqw = r.*cos(nu);
y_pqw = r.*sin(nu);

r_pqw = [x_pqw; y_pqw; zeros(size(x_pqw))];

figure;
axis equal;
grid on;
hold on;
set(gcf,'color',[1,1,1]);

xlabel('P axis(km)');
ylabel('Q axis(km)');
title('Molniya Orbit in Perifocal Frame');

plot(x_pqw,y_pqw,'b','LineWidth',1.2);
plot(0,0,'ko','MarkerFaceColor','k');

% h : specific angular Momentum vector
hSat = plot(x_pqw(1), y_pqw(1), 'ro', 'MarkerFaceColor','r');

for k = 1:length(nu)
    set(hSat, 'XData', x_pqw(k), 'YData',y_pqw(k));
    drawnow;
end



%ECI FRAME : draw the orbit
Omega = deg2rad(raan);
inc = deg2rad(i);
omega = deg2rad(aop);

ROT3 = @(ang)[ cos(ang) sin(ang) 0; -sin(ang) cos(ang) 0; 0 0 1];
ROT1 = @(ang)[ 1 0 0; 0 cos(ang) sin(ang);0 -sin(ang) cos(ang)];

A = ROT3(-Omega)*ROT1(-inc)*ROT3(-omega);
r_eci = A*r_pqw;



%Animation
figure;
axis equal;
grid on;
view(3); %3d axis
hold on;
set(gcf,'color','w');

xlabel('ECI X(km)');
ylabel('ECI Y(km)');
zlabel('ECI Z(km)');
title('Molniya Orbit in ECI Frame');

plot3(r_eci(1,:), r_eci(2,:), r_eci(3,:),'b','LineWidth',1);
plot3(0,0,0,'ko','MarkerFaceColor','r');
hSat_eci = plot3(r_eci(1,1), r_eci(2,1), r_eci(3,1),'ro','MarkerFaceColor','r');
for k=1:length(t)
    set(hSat_eci, 'XData', r_eci(1,k), 'YData', r_eci(2,k),'ZData',r_eci(3,k));
    drawnow;
end

 

아직 더러운 코드 같은데 ,, 

나중되면 더 나아지겠지 ? ? 

 

처음에는 true anomaly, v가 변함에 따라 달라지는 시뮬레이션을 그렸는데, 그렇게 그리니 원지점, 근지점에서의 속도가 반대가 되어버려서 시간에 따른 v를 구하고 R(위치)를 그려보았다. 

 

코드 실행 결과 

Perifocal 

 

eci