궤도적분 응용 : 3차원 궤도 적분과 검증

˚어제 완료했던 과제를 3차원 궤도 적분으로 이어서 진행해본다. 

 

추가로 받은 과제 

+ addpath를 사용하여 constants가 정의된 객체에서 gm(중력상수, mu = 3.986*1e5 km**3/s**2)을 정의하기 

addpath('my_path')

constants; %객체 가져옴 
constants().mm %객체 접근 후 파라미터 가져옴

 

2D(Perifocal 좌표계)에서 3D(ECI 좌표계)로 변환을 이어나가기 위해서 기존 4개를 받던 행렬에서 6개로 늘리고 

i, Omega, omega를 통해서 좌표 변환을 마무리해야한다. 

 

전체 Flow 

 

main : tbp_3d.m 

subroutine 1 : eulode_3d.m

subroutine 2 : tbp_3d_dydt.m 

 

1. 2D → 3D 행렬 변환 

하나씩 차근차근해보는 이유가 그냥 냅다 고쳐버리면 1000%문제가 생긴다. (필자는 초보라 그럼. 이 글을 보시는 고수 분들은 아닐수도 있음) 2D를 3D로 변환했을 때 z=0이고 xy평면에 붙어있어야하는데 이렇게 일자 모양이 나왔다. 

 

dydt = @(t,Y)[
    Y(3);
    Y(4);
    0
    -mu*Y(1)/((sqrt(Y(1)^2+Y(2)^2))^3) ; %a of x axis
    -mu*Y(2)/((sqrt(Y(1)^2+Y(2)^2))^3) ; %a of y axis
    0 %a of z axis
];


dydt = @(t,Y)[
    Y(4);
    Y(5);
    Y(6); %v_Z
    -mu*Y(1)/((sqrt(Y(1)^2+Y(2)^2+Y(3)^2))^3) ; %a of x axis
    -mu*Y(2)/((sqrt(Y(1)^2+Y(2)^2+Y(3)^2))^3) ; %a of y axis
    -mu*Y(3)/((sqrt(Y(1)^2+Y(2)^2+Y(3)^2))^3) ; %a of z axis
];

바보같은 실수를 했다.

인덱스도 밀리고 나중에 함수를 이용해서 계산을 해야하는데 초기값인 0을 그냥 저기 넣어주다니 .. 

 

이것만 잘 고쳐주면 기초 공사 완성 !

이제 이 3D틀을 가지고 기울기, RAAN, Argument of perigee를 통해

2차원 Perifocal 궤도를 3차원 ECI 궤도로 변환해주면 된다. 

 

 

 

 

 

2. Perifocal → ECI 좌표 변환 

Euler적용해서 만드는 것도 헷갈리는데 subroutine이 2개 더 있으니 헷갈려서 FlowChart를 만들어보았다. 

flowchart

 

좌표 변환을 위해 각 Epoch마다 x,y,z좌표를 바꾸어야하니 euler법을 적용한 솔버인 eulode 함수를 변형시킨 eulode_3d 함수를 만들었다. 

%ECI Frame
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);

 

이렇게 Frame 적용할 수 있게 매개변수를 받아와서 변환 행렬도 만들어주고 .. 

 

개같이 멸망 

이게 뭣이란 말임 !! 

 

그런데 의외로 쉽게 해결되었다 .. 변수 하나 더 설정하고 변수에 넣어주면 됨 .. ^^ (시간만 날려따 ㅎ.ㅎ)

 

%Before
for i = 1:n-1

    slope = dydt(t(i),y(i,:));
    dt = t(i+1)-t(i);
    y(i+1, : ) = y(i, : ) + slope' * dt;

    % each epoch - how can I get 3D coordinate ?
    y(i, 1:3) = y(i,1:3) * A'
    
end

%After
for i = 1:n-1

    slope = dydt(t(i),y(i,:));
    dt = t(i+1)-t(i);
    y(i+1, : ) = y(i, : ) + slope' * dt;

    % each epoch - how can I get 3D coordinate ?
    y_eci(i,:) = y(i,1:3)*A';

end

 

변수 업데이트할 때는 꼬옥... 새로운 변수 만들어서 넣어주기 

(오염도 방지하고 데이터가 꼬이지 않아 정확한 결과가 나온다 !) 

각각 h= 10, h = 100일 때

이제서야 잘 나오는 ECI좌표계에서의 궤도 

 

 

 

3. 검증 

그럼 내가 예측한 궤도가 얼마나 잘 맞는가에 대해 확인하는 것까지 RK4를 사용하여 검증해보자 !  

Runge-Kutta방식을 간략히 설명해보자면, 오일러 방식에서 평균을 내어 조금 더 정확하게 예측하는 방법이다. 

(이 부분에 대해서는 수치해석을 배우며 좀 더 자세히 다루어보겠다.) 

 

매트랩 내장함수에는 RK를  사용할 수 있는 ode45솔버 함수가 정의되어 있다. 

[t,rk4_result] = ode45(@(t,Y) tbp_3d_dydt(mu,t,Y), tspan, y0);

그래서 이렇게 정의해주면 바로 사용 가능 ! 

 

 

이제 3차원에서 63.4˚ 기울어져있는 원을 예측한 값을 볼 수 있다. h값(step값)을 줄이면 빨간 원(해석해)에 가까워지는 것을 확인할 수 있다. 당연히 반대로 h를 크게 설정하면 거의 경로를 이탈해버린다. 

 

추가적으로 Specific Mechanical Energy 식을 활용하여 에너지 보존이 잘 되고 있나...를 그래프로 그려보았다. 

역시 Euler법으로 예측한 에너지는 발산하지만, RK4에서는 잘 보존되고 있는 것을 확인할 수 있다.