그 동안 강의를 통해 배운 내용으로 화성과 지구의 날에 따른 거리를 구하는 프로그램을 짜보았다. 거리를 구하려면 먼저 태양을 초점으로 하는 화성과 지구 각각의 날짜에 따른 위치를 알아야한다. 이는 화성과 지구의 궤도를 구해서 날짜에 따른를 구한다는 말이 된다.궤도를 구하기 위해서 전통적인 6가지 요소를 사용하게 되는데, 그 첫 번째 요소로 장반경와 eccentricity, 그리고 perigee 통과시간가 있다. 이 세 요소가 주어지면 궤도면에서의 위성의 위치를 알 수 있다. 두 번째 요소 그룹은 궤도면을 표현하는데 사용하는데, 이는 ascending node의 적경(right ascension), 궤도면 경사각, 그리고 perigee의 위치를 나타내는이다. 이에 대한 자세한 내용은 모 기관의 홈페이지에서 찾았는데 그 값은 다음 표에 나타내었다.(km)(deg)(deg)(deg)(deg)Earth0.01671750100.464571102.9376810Mars0.09338651.8496914-4.55343-23.94362949.5595389table 1. 2000.1.1이 값들은 모 기관에서 측정한 자료들이기 때문에 내가 이용하기에는 적당치 않아 쓸 수 있는 수치로 바꾸어야한다. 위 자료를 통해 궤도를 구하는 방법은, mean anomaly를 통해 eccentric anomaly를 구해주고 거기서 나온 E로와를 구해주면 r을 구할 수 있고 이를 x, y 좌표 평면으로 구해주어 마지막으로 rotation을 해주면 각각 행성의 궤도를 구할 수 있다. 그 다음 각각의 시간에 따른 거리를 구해주면 된다.먼저 mean anomaly를 구하려면,;equ.1로와 mean anomaly를 구해 줄 수 있다. mean anomaly를 구하는 식에서가 없어지고가 나타난 이유는값을 찾기가 힘들므로의 값을 0으로 두면 mean anomaly의 값은 2000.1.1.의만큼 변위를 더해주어 그 위치를 시작점으로 하였다. 이를 통해 mean anomaly를 구할 수가 있었다.다음으로 eccentric anomaly를 구해보면,equ.2equ.2와 같은 비선형 방정식을 통해서 E를 구해야 한다. 이는 수치해석 방법을 통해서 구할수 있는데 나는 open-method의 일종인 fixed-point method를 이용해서 풀었다.그리고와는 책에 나오는 공식인,equ.3equ.4라는 두 공식으로 구할 수 있는데 여기서를 직접 구하지 않고 sin과 cos형태로 구하는 이유는 r을 구할때나 x, y좌표로 변환 시켜줄 때는가 필요하지 않고와가 필요하기 때문이다.x, y좌표로 변환 시키는 건 잘 알고 있으므로 설명하지 않고 r을 구하는 공식에 대해 써보면,equ.5로 쓸 수 있다. 마지막으로 궤도의 rotation은 찾은 자료에 나와 있는 것을 이용했다. 이는 첨부 자료의 8-34식을 보면 된다. 이러한 이론을 이용하여 알고리즘을 짰다.clear%%Initial valuesT=1:input('T =');t_e=100.46457166*pi/180; t_m=-4.55343205*pi/180 ; %Mean longitudea_e=1.496*10^8; a_m=2.2792e+008; %Semi-major axise_eccent=0.0167175; m_eccent=0.0933865; %% Orbit eccentricitymu=1.32715*10^11; %%Gravitational parameterW_e=102.93768193*pi/180; W_m=-23.94362959*pi/180;%Longitude of the ascending nodew_e=102.93768193*pi/180; w_m=-73.5031682*pi/180;o_e=0; o_m=49.55953891*pi/180; %Argument of perihelioni_e=0; i_m=1.84969142*pi/180; %Orbit inclinationt=T*3600*24;%Compute the mean anomalyM_e=t_e-W_e+sqrt(mu/a_e^3)*t;M_m=t_m-W_m+sqrt(mu/a_m^3)*t;%Compute the eccentric longituden_M=length(M_e);E_e=zeros(1,n_M);E_m=zeros(1,n_M);for a=1:n_M;x_e0=1;x_m0=1;dx=1;while 1if dx