궤도역학 실험실
케플러 궤도 전파부터 좌표변환, 실측 지상국 가시성, 호만 전이, J2/3체/태양복사압 섭동, 란베르트 문제, 다중 위성 성좌 커버리지, 위성간 링크 가시선, 행성간 궤적(Patched Conic), 도킹/랑데부(Clohessy-Wiltshire), 자세동역학/PID 제어, CR3BP까지, 궤도역학 자체를 1급 시뮬레이션 대상으로 처음부터 직접 구현한 27개의 점진적 시뮬레이션.
한눈에 보기
핵심 물리 검증 한눈에 보기
이 프로젝트의 핵심 차별점은 근사 대신 실제 궤도역학 공식을 직접 구현하고, 그 결과가 잘 알려진 물리 법칙 및 표준 문헌값과 일치하는지 매 스크립트에서 검증한다는 것이다.
▸ 기초 전파 — 궤도의 가장 기본이 되는 수학/좌표계 4개 · 01–04
원궤도는 초기값 자체가 정답, 이심률이 커질수록 반복이 늘어남
2체 문제의 핵심 물리 법칙이 수치정밀도 수준에서 보존됨
RK4가 이론대로 4차 정확도(O(h^4))로 수렴함을 확인
ECI→SEZ 회전 부호를 실제로 겪은 버그를 통해 고정
▸ 임무 설계 — 궤도 전이와 접촉, 운용 시나리오 9개 · 05,06,09,12–15,22,27
사각파 근사가 아니라 실제 궤도 전파로 계산된 접촉 창
01번 실제 궤도 전파로 전이 도착 시각을 재검증함
유니버설 변수 반복법과 닫힌 공식이 서로 다른 방법으로 같은 답에 수렴
06번 호만 전이 공식을 태양 GM으로 그대로 재사용해 실제 임무 값 재현
Clohessy-Wiltshire 선형화는 근접 운용에서만 정확함을 실측으로 확인
노이즈 낀 레이더 관측만으로 02번 궤도요소를 거꾸로 복원, 관측이 늘수록 안정화
델타-V는 161.00m/s로 호만(160.76m/s)과 비슷해도 시간은 압도적으로 오래 걸림
06번 호만 전이와 21번 PID 자세제어를 새 물리 없이 하나의 임무 타임라인으로 잇는다 — 델타-V를 다 쓴 궤도전이도 몇 초~몇십 초의 자세 재정렬 없이는 쓸모없는 자산으로 남는다
14번이 "가우스법 등은 범위가 커진다"며 미뤘던 일반 타원궤도 결정 — 02번 범용 궤도요소 변환을 재사용해 14번이 다루지 못한 이심률·근점편각까지 복원
▸ 섭동 — 완전한 2체 문제에서 벗어나는 요인들 8개 · 07×2,17,18,23,24,26,28
몰니야 궤도가 이 경사각을 쓰는 이유를 수치로 재현
지구 편평도(J2)가 궤도면 자체를 서서히 회전시킴
07번 J2 세차의 잔여 오차를 평면 변경 기동으로 주기 보정 — 1년에 6회, 기동당 13.17m/s
15번 저추력 상승 나선의 거울상 — 대기항력이 에너지를 빼앗아 재진입 근처에서 하강 속도가 2.20배 가속
17번이 "모델링하지 않는 요인"으로만 언급했던 태양/달 중력을 실제로 계산한다 — 달의 섭동이 태양보다 2.24배 크고, LEO에서는 무시 가능하지만 GEO 근처에서는 J2를 능가한다
9번(대기항력)과 거울상 — 고도가 높아지면 대기는 사라져도 태양빛은 남는다. 저고도에서는 항력이 압도적이지만, 이 고도를 넘으면 SRP가 유일하게 남는 비중력 섭동이 된다
24번의 캐논볼 모델(자세 무관 단면적)을 넘어, 21번의 쿼터니언 자세를 결합해 평판 패널의 코사인 법칙 유효 단면적을 계산 — 텀블링하면 SRP가 진동하지만 태양지향 제어는 캐논볼 최댓값 근처를 유지
23번이 "가우스 섭동방정식 기반 평균화가 필요한 별도 범위"로 남겨둔 가속도→RAAN 변화율 변환을 구현 — 23번의 가속도 크기 비율(144%)과 이 스크립트의 RAAN 변화율 비율은 서로 다른 물리량(검증오차 0.0225%로 07번과 교차검증)
▸ 다중 위성 — 성좌 커버리지와 위성간 링크 2개 · 10,11
"평면을 나누면 항상 좋다"는 직관이 지점 하나 기준에서는 뒤집힘(감추지 않음)
두 위성을 잇는 선분과 지구 중심 사이 최근접 거리로 기하학적 판정
▸ 3체 문제 — 라그랑주점과 그 주위 주기궤도 2개 · 19,20
12번이 피했던 진짜 3체 문제를 원형 제한 3체 문제(CR3BP)로 정면 돌파 — 회전좌표계 선형화 고유값과 RK4 직접 적분이 같은 결론을 준다
19번의 라그랑주점이 "점"이었다면, 이건 그 주위를 도는 "주기궤도" — 선형 고유벡터 추정치를 상태천이행렬 기반 미분수정으로 완성한다
▸ 자세동역학 — 궤도가 아닌 위성 자체의 회전 3개 · 16,21,25
토크 없는 자유회전에서 중간축만 작은 섭동에도 크게 발산함을 오일러 방정식으로 재현(중간축 정리)
16번이 "순수 동역학만, 제어는 범위 밖"이라 정했던 경계를 의도적으로 뒤집는다 — 쿼터니언 오차 기반 PD 피드백이 중간축 불안정 텀블링을 거의 완전히 정지시킨다
22번이 "반작용휠 디새추레이션 같은 추가 델타-V는 모델링하지 않음"이라 명시한 갭을 채운다 — 휠 포화 후 디새추레이션은 작지만 0이 아닌 델타-V가 든다
스물여덟 단계의 확장
각 단계는 처음부터 새로 만들지 않고 이전 단계의 인프라를 그대로 재사용한다 — 예를 들어 05, 06, 07번은 모두 01번의 solve_kepler_equation과 propagate_orbit을 import해서 쓴다.
▸ 기초 전파 — 궤도의 가장 기본이 되는 수학/좌표계 5개 · 01–04,08
케플러 방정식과 궤도 전파
M = E - e*sin(E)를 뉴턴-랍슨으로 풀어 시간과 궤도 위치를 잇는다. 이후 모든 스크립트가 이 모듈의 핵심 함수를 재사용한다.
궤도요소 ↔ 상태벡터, 보존량
궤도면 2차원 위치를 3번의 회전으로 3차원 ECI 상태벡터로 확장하고, vis-viva 방정식과 비에너지/비각운동량 보존을 검증한다.
2체 문제 RK4 수치적분
해석해 대신 운동방정식을 RK4로 직접 적분해 01/02번의 해석해와 대조하고, O(h^4) 수렴 차수와 에너지 드리프트를 확인한다 — 07번 섭동 모델의 신뢰성 기준선이 된다.
좌표계 변환 (ECI/ECEF/SEZ)
지구 자전(ECI→ECEF)과 관측자 위치(위경도)를 반영해 궤도 위치를 방위각/고도각으로 바꾼다. 05번이 이 변환을 재사용한다.
궤도, 접촉 창, 전이궤도 그리기
시뮬레이션 없이 앞선 스크립트들의 CSV만 읽어 궤도 형상, 고도각 시계열, 접촉 창 간트차트, 전이궤도, J2 세차를 그래프로 그린다.
▸ 임무 설계 — 궤도 전이와 접촉, 운용 시나리오 9개 · 05,06,09,12–15,22,27
실측 지상국 가시성
사각파 근사 대신, 01번의 실제 궤도 전파 + 04번의 좌표변환으로 매 순간 고도각을 계산해 실측 접촉 창(AOS/LOS)을 찾는다.
호만 전이
최소 에너지로 궤도를 변경하는 두 번의 임펄스를 계산하고, LEO→GEO 델타-V가 교과서 값과 일치하는지, 01번 전파로 전이 도착 시각이 정확한지 검증한다.
란베르트 문제
궤도요소가 주어졌을 때 위치를 구하던 순문제를 뒤집어, 두 위치·비행시간으로 속도벡터를 구하는 역문제를 유니버설 변수·스텀프 함수로 푼다. 06번 호만 전이와 교차검증한다.
행성간 궤적 (Patched Conic)
지구 하나의 중력장을 벗어나, 지구 탈출-태양 중심 전이-화성 포획을 세 개의 독립된 2체 문제로 이어붙인다. 06번의 호만 전이 공식을 태양 GM으로 그대로 재사용하고, 중력권(SOI) 반지름이 전이거리에 비해 실제로 작은지 스스로 검증한다.
도킹/랑데부 (Clohessy-Wiltshire)
09번의 큰 스케일 전이 다음 단계 — 목표 위성을 원점으로 한 상대좌표계에서 근접 운용 궤적을 선형화된 닫힌 해(상태천이행렬)로 구한다. 제로 드리프트 조건 없이는 진행 방향으로 선형 발산한다는 것과, 이 근사가 상대거리가 커지면 실제 비선형 궤도에서 벗어난다는 것을 01/02번 재전파로 검증한다.
궤도 결정 (Orbit Determination)
02번(궤도요소→상태벡터)과 정반대 방향 — 지상 레이더의 노이즈 낀 거리/방위각/고도각 관측값"만"으로 원궤도의 반장축·경사각·RAAN을 역산한다. 04번 관측 기하 변환을 정확히 뒤집어 순간 위치를 복원하고, 여러 시점의 위치에 SVD로 궤도면을 최소자승 피팅한다. 관측 노이즈가 커질수록, 관측 개수가 적을수록 추정이 불안정해지는 것을 직접 확인한다.
저추력 전기추진 궤도 전이 (Low-Thrust Transfer)
06번 호만 전이가 "순간 임펄스 두 번"으로 궤도를 바꿨다면, 이 스크립트는 정반대 축 — 이온엔진 같은 전기추진으로 훨씬 작은 가속도를 몇 시간〜며칠에 걸쳐 연속으로 가하는 전이를 다룬다. 03번 RK4 적분기의 중심력 가속도에 접선방향(진행방향) 추력 가속도를 더해 원궤도가 나선형으로 상승하는 과정을 직접 적분하고, 닫힌 형태 근사(Edelbaum 근사, Δv=|v1-v2|)와 대조해 검증한다. 06번과 같은 반지름 쌍에 대해 호만 전이와 델타-V/시간을 나란히 비교하면, 델타-V는 비슷한 수준이어도 걸리는 시간은 자릿수가 다르다는 실제 임무 설계의 트레이드오프가 드러난다.
궤도+자세 통합 임무 시나리오
01~21번은 궤도역학(병진운동)과 자세동역학(회전운동)을 서로 독립된 두 영역으로 다뤘다. 이 스크립트는 새 물리를 전혀 추가하지 않고, 이미 완성된 두 스크립트(06번 호만 전이, 21번 PID 자세제어)를 하나의 임무 타임라인으로 이어붙인다 — 위성이 주차궤도에서 임의의 텀블링 자세로 시작해, 호만 전이로 목표 임무궤도에 오른 뒤("on station"), PD 제어기로 지구지향 자세로 재정렬한다. 궤도 전이(수십 분)와 자세 재정렬(수십 초)의 시간축 스케일이 1000배 가까이 달라, 자세 재정렬 단계는 궤도 전이 완료 시각을 오프셋으로 기록해 별도 시간축을 하나의 타임라인에 잇는다. 자세 재정렬 단계의 첫 쿼터니언이 궤도 전이 완료 시점에 선언한 바로 그 임의 자세와 정확히 같다는 것으로 상태가 단계 경계를 실제로 통과함을 확인하고, 이 스크립트의 순수 토크 제어 모델에서는 자세 재정렬이 델타-V를 전혀 소모하지 않음을 보여 두 단계가 물리적으로 독립된 자원(추진제 vs 제어 시간)을 소모한다는 것을 증명한다. 지상국 접촉창(05번) 단계는 처음엔 포함하려 했지만, 궤도면 전체 기하를 다시 복제해야 해서 이 스크립트의 핵심 주장과 무관한 범위 확장이라 판단해 뺐다.
가우스법 궤도결정 (Gauss's Method)
14번(orbit_determination.py)은 "이심률이 0인 원궤도만 다룬다 — 완전히 일반적인 케플러 궤도 피팅(가우스법 등)은 라그랑주 계수와 고차 다항식 풀이가 추가로 필요해 범위가 커진다"고 명시적으로 범위를 그었다. 이 스크립트는 바로 그 가우스법을 구현해, 방위각/고도각만(거리 없이) 측정한 3회 관측값으로 일반 타원궤도를 복원한다. 세 시선벡터와 관측자 위치로 세운 라그랑주 계수(f,g) 3차 테일러 전개가 미지의 속도를 소거해, 중간 관측 시점의 거리 r2 하나짜리 8차 다항식으로 귀결된다는 것이 핵심 구조다. 공식 자체는 Curtis 교재 Example 5.11의 정확한 수치로 먼저 재현해 검증했다(r2≈8402km, 기대값 ~8000km 근방과 합치). 그런데 합성 LEO 관측 시나리오로 넓게 스윕하자, "반복 없는" 1차 패스는 세 시선벡터의 스칼라 삼중곱 D0(공면성 척도)가 조금만 작아져도 계수가 수십만 단위로 터지며 어떤 근도 참값과 맞지 않는 극심한 불안정성을 보였다 — 특정 관측 간격 하나에서만 좁고 날카롭게 성공하고 그 근처에서도 수십~수백% 오차로 무너졌다. 계획 단계에서는 이 반복 자체를 "범위 밖"으로 미루려 했으나, 실측이 그 계획을 뒤집었다 — 최초 r2 추정값으로 라그랑주 계수를 반복 재계산하는 단순 고정점 반복(스텀프 함수 기반 완전 보편변수 공식 없이) 하나만 추가해도, 넓은 관측 간격 구간 전체에서 반장축 오차가 1~3%대로 안정되는 극적인 개선을 확인했다. 수렴한 위치/속도벡터를 02번의 범용 궤도요소 변환에 그대로 넘기면 14번이 전혀 다루지 못했던 이심률·근점편각까지 포함한 6개 요소 전부가 나온다. D0가 작은 관측 기하에서는 가우스법이 근본적으로 불안정하다는 것을 감추지 않고, 전용 데모로 그 한계를 직접 보여준다(9번 란베르트의 180도 특이점과 같은 이 프로젝트의 태도). 이것으로 14번의 원궤도 한계를 넘어 일반 타원궤도 궤도결정이 완성됐다.
▸ 섭동 — 완전한 2체 문제에서 벗어나는 요인들 7개 · 07,17,18,23,24,26,28
J2 섭동 (RAAN/근점편각 세차)
지구 편평도가 일으키는 궤도면 회전을 계산한다. 임계경사각에서 근점편각 세차가 0이 되는 것, 극궤도에서 RAAN 세차가 0이 되는 것을 직접 확인한다.
궤도 유지 비용 (Station-Keeping)
07번(J2 섭동)은 궤도면(RAAN)이 서서히 흐트러진다는 것을 보여줬지만, 그 드리프트를 방치한 결과만 계산했다. 실제 위성은 이걸 방치하지 않고 주기적으로 작은 기동(임펄스)으로 되돌리는데, 07번의 예측 세차율 자체는(태양동기궤도라면) 오히려 의도된 값이라 보정 대상이 아니다 — 실제로 보정하는 것은 그 예측과 실측 사이의 작은 잔여 오차뿐이다. 이 스크립트는 RAAN을 되돌리는 데 필요한 평면 변경(plane-change) 델타-V 공식(Δv=2*v*sin(각도/2), 06번 호만 전이의 반지름 변경 공식과는 다른 종류)을 새로 정의하고, 허용 오차에 도달할 때마다 기동해 RAAN 오차가 톱니파 형태로 유계에 머무는 것을 확인한다. 허용 오차를 10배까지 바꿔도 임무 기간 총 델타-V 예산은 거의 그대로다 — 자주 보정하면 한 번에 쓰는 델타-V는 작지만 횟수가 늘고, 드물게 보정하면 그 반대라 서로 상쇄되기 때문이다.
궤도 재진입/대기항력 감쇠 (Orbital Decay)
15번(저추력 전기추진)이 접선방향 연속 추력으로 원궤도를 서서히 밀어 올렸다면, 이 스크립트는 정반대다 — 저궤도 위성은 완전한 진공이 아니라 극도로 희박한 대기 속을 지나가므로, 속도 반대방향으로 작용하는 미세한 항력이 계속 에너지를 빼앗아 고도가 서서히 떨어지다 결국 재진입한다. 지수함수 대기밀도 근사(ρ(h)=ρ0·exp(-(h-h0)/H))와 탄도계수 기반 항력 감속(a=-0.5·ρ(h)·v²/BC)을 15번의 "접선 추력" 항을 부호만 반전한 구조로 재사용한다. 대기밀도가 지수함수이므로 재진입 직전으로 갈수록 하강이 급격히 가속화되는 비선형 패턴이 나타나고, 탄도계수가 클수록(무겁고 단면적이 작은 위성) 수명이 길어지며, 초기 고도가 조금만 높아져도 수명이 기하급수적으로 늘어난다 — "폐기위성을 낮은 고도로만 유도해도 자연스럽게 몇 년 안에 재진입한다"는 실제 우주 쓰레기 정책의 물리적 근거를 수치로 보여준다.
3체(태양/달) 섭동
17번(station_keeping.py)의 docstring은 "세차율 계산 자체의 불확실성(더 높은 차수의 섭동, 대기항력, 태양/달의 중력 등 이 스크립트가 모델링하지 않는 요인들)"이라며 태양/달의 중력을 명시적으로 범위 밖으로 남겨뒀다. 이 스크립트는 그 요인을 실제로 계산 가능한 값으로 만든다 — 지구 중심 좌표계에서 위성이 받는 태양/달의 중력은 지구 자신도 같은 천체에 끌리므로 그 공통항을 뺀 상대(차분) 섭동이어야 한다(a_3rd = mu_body*[(r_body-r_sat)/|r_body-r_sat|^3 - r_body/|r_body|^3]) — 위성이 지구 중심에 있으면 두 항이 정확히 상쇄돼 0이 되는 항등식으로 공식 자체를 검증한다. 태양의 GM은 달의 2700만 배지만, 거리의 세제곱에 반비례하는 상대 섭동 공식 때문에 지구에 훨씬 가까운 달의 섭동이 실제로는 태양의 약 2.2배 크다는 교과서적 결과를 직접 재현한다. 07번(J2 섭동)과 같은 단위(가속도 크기)로 비교하면, LEO에서는 J2가 압도적이지만(비율 0.01% 수준) 고도가 높아질수록 격차가 줄어 GEO 부근에서는 3체 섭동이 J2를 능가한다 — 실제 GEO 위성 운용에서 태양/달 섭동이 J2보다 더 중요한 보정 대상인 이유를 수치로 보여준다. 태양과 달 모두 이심률 0인 원궤도, 지구 적도면과 같은 평면으로 단순화한다(황도 기울기와 달의 궤도경사는 무시) — 방향 정확도는 희생하지만 이 스크립트가 보여주려는 핵심(크기와 시간에 따른 변화)은 그대로 보존된다. 17번의 RAAN 세차율(rad/s)과 이 스크립트의 가속도(km/s²)는 단위가 달라 직접 변환하지 않는다 — 가우스 섭동방정식 기반 평균화가 필요한 별도 범위이기 때문이다.
태양복사압(SRP) 섭동
9번(orbital_decay.py)은 저궤도에서 지배적인 대기항력 감쇠를 다뤘다. 이 스크립트는 그 거울상이다 — 대기가 사실상 없는 고궤도에서는 태양복사압(SRP)이 남는 유일한 비중력 섭동이 된다. 단순화된 구형(cannonball) 모델에서 가속도는 a_srp = P_sr·Cr·(A/m)·unit(r_sat-r_sun)으로 주어진다 — 방향은 태양→위성 방향(태양 반대쪽)이며, 부호가 뒤집히면 태양 쪽으로 끌리는 물리적 오류가 되므로 이 부호 자체를 항등식 테스트로 검증했다. 처음엔 23번처럼 J2나 3체 섭동과 비교하려 했지만, 실제로 계산해보니 SRP는 태양까지의 거리(~1AU)에만 의존해 위성의 지구 중심 거리와 거의 무관하다는 것을 발견했다 — GEO에서도 J2 대비 비율이 1% 미만에 불과해 23번 같은 흥미로운 교차점이 나오지 않았다. 대신 9번의 지수함수 대기밀도 모델을 그대로 재사용해 항력과 비교하니, 200km에서는 항력이 압도적이지만 3000km에서는 SRP가 압도적이며 약 550~800km 사이에서 실제 교차점이 나온다는 것을 확인했다 — 실제 위성 궤도수명 연구가 약 600~800km 이하에서만 항력을 신경 쓰는 이유와 정확히 일치한다. 또한 태양-지구 축 기준 원통형 그림자 근사(두 조건 모두 필요: 태양 반대쪽 AND 축에서 수직거리가 지구 반지름 이내)를 도입해, 궤도면을 태양 방향과 일치시킨 특수 기하(베타각 0도)에서 궤도의 약 36%가 그림자에 들어가 SRP가 정확히 0이 되는 것을 직접 적분으로 확인했다 — 이 식 비율은 일반적인 궤도의 주장이 아니라 이 특수 기하에만 해당한다는 것을 명시한다.
자세 의존 태양복사압(SRP)
24번(solar_radiation_pressure.py)은 명시적으로 캐논볼 모델(구형 위성, 자세 무관 단면적)만 다뤘다 — 그 docstring이 "실제 위성의 비구형 형상에 따른 자세 의존 단면적 변화는 다루지 않는다"고 범위를 그었다. 이 스크립트는 그 갭을 연다 — 21번(pid_attitude_control.py)의 쿼터니언 자세 상태를 24번의 SRP 가속도 계산에 결합해, 평판 태양전지판처럼 자세에 따라 태양을 향한 유효 단면적이 달라지는(코사인 법칙, A_eff=A_max·max(0,cosθ)) 경우를 다룬다. 24번의 srp_acceleration/srp_acceleration_with_shadow는 area_to_mass_ratio_m2_kg를 단순 스칼라 곱으로만 쓴다는 것을 확인했으므로, 24번 자체는 전혀 수정하지 않고 자세로부터 계산한 "유효 A/m"만 그 함수에 그대로 넘긴다 — 22번이 06번/21번을 수정 없이 얇은 오케스트레이션 계층으로 결합한 것과 같은 설계 원칙이다. 쿼터니언으로 몸체좌표계 벡터를 관성좌표계로 회전시키는 함수(쿼터니언 샌드위치 v_inertial=q⊗[0,v]⊗conj(q))는 21번에도 없던 이 프로젝트 최초의 인프라라, 항등 회전과 알려진 90도 회전 두 경계 케이스로 먼저 검증했다. 토크 없는 자유회전으로 계속 텀블링하는 위성은 패널이 태양을 향했다 등졌다 반복해 유효 SRP 가속도가 뚜렷이 진동하지만(캐논볼의 "항상 일정"이라는 가정이 깨짐), 21번의 PD 제어로 패널을 태양 쪽으로 계속 향하게 유지하면 유효 SRP가 캐논볼 최댓값의 95% 이상에서 거의 일정하게 유지된다 — 태양지향 자세제어가 전력 생산뿐 아니라 SRP 예측 가능성에도 영향을 준다는 현실적 함의다. 궤도 주기(약 5829초)와 텀블링 회전 주기(약 13~63초)의 격차가 22번이 겪은 궤도전이-자세정착의 약 1000배 격차와 달리 두 자릿수(약 50~450배) 안에 들어온다는 것을 사전에 계산으로 확인해, 22번의 phase 컬럼 long-format CSV 없이 하나의 공유 시간축·단일 RK4로 궤도와 자세를 함께 적분했다. 5개 궤도 주기를 적분하면 텀블링 위성의 경로가 캐논볼 모델 경로에서 서서히 벌어지지만, SRP 자체의 절대 크기가 원래 작아(~1e-10km/s²) 그 편차도 반장축(7000km) 대비 0.0002% 수준(약 11.6m)에 머문다 — 과장하지 않는다.
GEO 통합 궤도유지 (J2+3체섭동)
17번(station_keeping.py)은 J2 세차의 "잔여 오차"만 보정하는 LEO/SSO 전용 궤도유지 비용을 계산했고, 23번(third_body_perturbation.py)은 GEO 고도에서 3체 섭동/J2 "가속도 크기" 비율이 약 144%(J2보다 큼)라는 것을 보였지만, 그 가속도를 17번이 쓰는 RAAN 변화율(rad/s)로 변환하는 작업은 "가우스 섭동방정식 기반 평균화가 필요한 별도 범위"라며 명시적으로 남겨뒀다. 이 스크립트는 바로 그 변환을 구현한다 — 임의의 섭동가속도를 RSW(반경/진행방향/궤도면법선) 좌표계로 분해하면 RAAN 변화율은 법선 성분(a_h)에만 의존한다는 가우스 섭동방정식(dΩ/dt=(r·sinu)/(n·a²·√(1-e²)·sini)·a_h)을 구현하고, J2의 RSW 성분(Vallado 표준식)을 새로 유도해 07번의 기존 해석적 세차율 공식과 대조해 검증했다(0.0225% 오차로 일치 — 밑바닥부터 재구현한 공식이 이미 신뢰받는 공식을 재현한 가장 강력한 증거). 같은 검증된 공식을 23번의 3체 섭동 가속도에 적용하자, 144%(가속도 크기 비율)와는 전혀 다른 24~95%(RAAN 변화율 비율, RAAN 각도에 따라 변함)가 나왔다 — 둘 중 하나가 틀린 게 아니라, "가속도 크기"와 "그 가속도가 실제로 RAAN을 바꾸는 효율"이 원칙적으로 다른 양이기 때문이다. J2는 (a,e,i)에만 의존하는 순수 영년(secular) 효과라 세차 전체를 원하거나(SSO) 정밀 예측으로 사전 보상할 수 있지만, 3체 섭동 RAAN 변화율은 RAAN·시각에 의존하는 비영년 효과라(RAAN=0/180도에서 0, 90/270도에서 최댓값인 주기 패턴을 실측으로 확인) "의도적으로 받아들일" 설계 선택지 자체가 없다 — 그래서 이 스크립트는 J2는 17번과 같은 관례로 잔여만, 3체 섭동은 전체를 보정 대상으로 삼는 비대칭 설계를 채택했다. GEO에서 3체 섭동이 J2 잔여 대비 약 1000배 큰 것으로 나타나 GEO 궤도유지는 J2가 아니라 3체 섭동이 지배적인 보정 대상임을 보이고, RAAN 평균 시나리오로 계산한 연간 델타-V(약 127m/s)가 문헌에서 흔히 인용되는 GEO 남북 궤도유지 예산(약 45~50m/s/year)과 자릿수가 맞는 것을 sanity check로 확인한다 — 이심률·경도 드리프트 보정은 다루지 않는 단순화 모델이라 정밀도 주장은 아니다.
▸ 다중 위성 — 성좌 커버리지와 위성간 링크 2개 · 10,11
성좌 커버리지
워커 델타 패턴으로 여러 위성을 배치하고, 05번의 접촉 창 계산을 "하나 이상 보임" 조건으로 일반화해 커버리지 공백을 계산한다. "평면을 나누면 항상 좋다"는 직관이 지점 하나 기준에서는 성립하지 않는다는 것을 실측으로 확인한다.
위성간 링크(ISL) 가시선
지상국이 아니라 위성-위성 사이의 시선을, 선분과 구(지구) 사이의 최근접 거리라는 순수 기하 문제로 판정한다. 정반대편(반드시 가려짐), 같은 평면(대부분 보임) 경계 케이스로 검증한다.
▸ 3체 문제 — 라그랑주점과 그 주위 주기궤도 2개 · 19,20
3체 문제와 라그랑주점 (CR3BP)
01~18번은 모두 "하나의 중심천체 + 작은 섭동/추력/항력"이라는 공통 골격이었고, 12번(Patched Conic)은 태양-지구-화성이라는 진짜 3체 문제를 세 개의 독립된 2체 문제로 이어붙여 피했다. 이 스크립트는 12번이 피했던 문제를 원형 제한 3체 문제(CR3BP)로 정면으로 연다 — 지구-달처럼 두 개의 큰 천체가 만드는 회전 중력장 안에서, 질량이 무시할 만한 세 번째 물체가 정지해 있을 수 있는 5개의 평형점(라그랑주점)을 찾는다. 회전좌표계 운동방정식에 이 프로젝트 최초로 원심력과 코리올리 힘이 등장하고, 무차원화(질량비 mu 하나로 계 전체가 결정됨)를 쓴다는 점에서 이전 스크립트들과 근본적으로 다르다. L1/L2/L3(콜린스)는 이분법으로 푼 5차 방정식의 근, L4/L5(삼각)는 두 주천체와 항상 정삼각형을 이루는 닫힌 공식으로 구하고, 각 점에서 선형화한 야코비 행렬의 고유값으로 안정성을 판정한다 — 지구-달 질량비에서 L4/L5는 안정(고유값이 순허수), L1/L2/L3는 불안정(고유값에 양의 실수부)임을 확인하고, 같은 크기의 섭동을 준 RK4 직접 적분 궤적으로도 L4는 유계, L1은 훨씬 크게 발산하는 것을 재확인한다.
평면 리아푸노프 궤도 (Lyapunov Orbits)
19번은 라그랑주점의 위치와 선형 안정성까지만 다뤘다 — 그 스크립트 docstring이 명시적으로 "헤일로 궤도 같은 CR3BP의 실전 응용은 다루지 않는다"고 범위를 그었다. 이 스크립트는 그다음 단계로, 라그랑주점 자체가 아니라 그 주위를 도는 주기궤도 가족이 존재한다는 것을 보여준다 — 실제 L1/L2 임무(SOHO, JWST 등)가 정확히 라그랑주점에 위성을 두지 않고 그 주위의 작은 주기궤도에 두는 이유의 수학적 근거다. 3차원 헤일로 궤도는 Richardson 3차 근사 계수 전사 실수 위험이 커서, 범위를 평면(2차원) 리아푸노프 궤도로 좁혔다. 19번이 만든 4x4 평면 야코비안의 순허수 고유벡터(중심 방향) 섭동을 초기 추정치로 삼아, 상태천이행렬(STM)을 상태와 함께 20차원으로 결합 적분하고, x축 대칭(y=0에서 vx=0)을 이용한 단일 슈팅 미분수정으로 진짜 비선형 주기궤도를 완성한다. 진폭이 커질수록 실제 주기가 선형 예측(2π/ω)에서 점점 벗어나는 것을 직접 확인해, "선형 근사는 시작점일 뿐, 진짜 주기궤도는 미분수정으로 완성해야 한다"는 핵심 메시지를 수치로 보여준다.
▸ 자세동역학 — 궤도가 아닌 위성 자체의 회전 3개 · 16,21,25
자세동역학 (Torque-Free Rigid Body, 중간축 정리)
01〜15번은 모두 위성을 점 질량으로 다뤘다 — 궤도상의 위치/속도(병진운동)만 계산했지, 위성 자체가 어느 방향을 향하는지(자세)는 다루지 않았다. 이 스크립트는 처음으로 궤도역학과 독립적인 새 물리 영역, 즉 강체의 회전운동을 연다. 외부 토크가 전혀 없는 강체의 오일러 회전방정식을 03/15번과 같은 RK4 구조로 직접 수치적분해, 관성모멘트가 가장 크거나 가장 작은 축 주위 회전은 작은 섭동에도 안정적이지만 중간 크기 축 주위 회전은 섭동이 기하급수적으로 커져 축 자체가 뒤집힌다는 "중간축 정리"(테니스 라켓 정리)를 재현한다. 자세제어(PID 피드백)는 이 프로젝트에서 계속 범위 밖으로 남겨두고, 순수한 동역학(외부 입력 없는 자유회전)만 다룬다.
PID 자세제어 (Quaternion Attitude Control)
16번은 토크 없는 강체 자유회전(오일러 방정식)만 다뤘고, 그 스크립트 docstring이 명시적으로 "쿼터니언 자세 전파나 피드백 제어는 범위 밖으로 남겨둔다"고 선을 그었다. 이 스크립트는 그 경계를 의도적으로 뒤집어, 목표 자세로 수렴시키는 PD 피드백 제어 루프를 연다. 16번은 각속도만 적분했지 자세(방향) 자체는 추적하지 않았으므로, scalar-first 쿼터니언과 그 운동학(dq/dt)이 이 스크립트에서 처음 등장하는 새 인프라다. 현재 자세와 목표 자세의 오차는 q_err=conj(q_target)⊗q로 정의하는데, q와 -q가 같은 자세를 나타내는 이중성 때문에 q_err의 스칼라부가 음수면 전체 부호를 뒤집어야 한다 — 안 그러면 제어기가 먼 길로 돌아간다. PD 게인은 임의의 값이 아니라 원하는 정착시간과 감쇠비에서 2차 시스템 공식(ωn=4/(ζ·ts), Kp=2·I·ωn², Kd=2·ζ·ωn·I)으로 직접 역산한다. 가장 직접적인 증명은 16번과 정확히 같은 중간축 섭동 설정에 PD 제어기를 적용하는 것 — 제어 없이는 180도까지 반복적으로 텀블링하던 회전이, 제어가 있으면 몇 도 이내로 거의 완전히 정지한다.
반작용휠 모멘텀 저장과 디새추레이션
22번(orbit_raise_and_reorient.py)은 "자세 재정렬은 이 스크립트의 PD 토크 제어 모델에서 델타-V를 전혀 소모하지 않는다(반작용휠 디새추레이션 같은 추가 델타-V는 모델링하지 않음)"이라고 명시적으로 인정한 갭이 있다. 이 스크립트는 그 갭을 실제로 채운다 — PD 법칙이 계산한 "몸체가 원하는 토크"는 추력기가 아니라 휠 모터가 만들고, 뉴턴 제3법칙에 의해 몸체는 정반대 부호의 반작용 토크를 받는다(torque_cmd = -tau_body_desired). 처음 구현에서 이 부호를 반대로 둬서 각속도가 10초 안에 발산하는 버그를 실제로 겪었다 — 올바른 부호로 고치자 21번과 동일한 수렴이 재현됐다. 외부 토크가 없으면 몸체+휠 총 각운동량의 "크기"가 보존되는지가 가장 강력한 검증이다(벡터 자체는 몸체가 회전하며 방향이 바뀌지만 크기는 보존). 상수 외란(예: 보정 안 된 중력경도 토크)이 지속되면 휠 모멘텀이 거의 정확히 선형으로 쌓이는 동안 포인팅 오차는 거의 0으로 유지돼 "조용히 쌓이는 문제"가 된다 — 포화 시점은 용량/외란 비율과 정확히 일치하고(실측 100.05초), 포화 후에는 포인팅이 완만하게가 아니라 180도까지 급격히 무너진다(하드 실패). 포화 축 클램프도 처음엔 "모멘텀을 줄이는 방향이면 통과"시키는 방향 인식 방식을 시도했는데, 자세 오차가 180도 근처(특이점)에서는 PD 법칙이 아주 큰 토크를 명령해 한 적분 스텝 만에 반대 극성으로 오버슈트하는 버그를 실측으로 발견해, 포화된 축은 방향 무관하게 전면 차단하는 방식으로 바꿨다. 포화 직후 신속하게 디새추레이션하면(몸체 각속도도 함께 리셋) 포인팅이 거의 완벽하게 유지되는 반면, 10초 넘게 방치하면 회복에 실패한다는 것도 확인했다. 디새추레이션 비용은 1.667mm/s로, 22번 기본 호만 전이(207.69m/s)의 약 0.0008%다 — 작지만 "정확히 0"이라는 22번의 주장에 현실적인 예외를 만든다. 외란이 멈추지 않는 한 휠은 같은 주기(용량/외란)마다 다시 포화되므로, 디새추레이션도 17번의 궤도 유지처럼 주기적으로 반복돼야 한다. 이것으로 자세동역학 3종(자유회전/PID 제어/휠 디새추레이션)이 모두 갖춰져 완성됐다.
결과 그래프
물리 법칙/표준값 대조 요약
각 스크립트가 자체 pytest로 검증하는 핵심 주장을 한 곳에 모았다 — 근사가 아니라 실제 공식/수치적분/알려진 문헌값과의 일치를 확인한다.
| 검증 항목 | 측정 기준 | 계산값 | 기준값 | 판정 |
|---|---|---|---|---|
| 케플러 방정식: 원궤도(e=0) 즉시 수렴 | 뉴턴-랍슨 반복 횟수 | 1 | 1 | 일치 |
| 2체 문제: 비에너지/비각운동량 보존 | 궤도 전체 최대 오차 | 0.00e+00 | < 1e-6 | 보존됨 |
| RK4 수치적분: O(h^4) 수렴 차수 | 스텝 절반 시 오차 감소율 | 약 16.5배 | 약 16배 | 근접 |
| 좌표변환: 천정 고도각 | 고도각(도) | 90.0000 | 90.0 | 일치 |
| 호만 전이: LEO→GEO 델타-V | 총 델타-V (km/s) | 3.893 | 약 3.9 | 교과서 값과 일치 |
| J2 섭동: 임계경사각(63.4도) 근점편각 세차 | 세차율(도/일) | 0.000000 | 0.0 | 이론과 일치 |
| 란베르트 문제: 06번 호만 전이와 교차검증 | 델타-V 오차(km/s) | 0.0000 | < 0.05 | 닫힌 공식과 일치 |
| 위성간 링크(ISL): 지구 정반대편 경계 케이스 | 시선-지구중심 최근접 거리(km) | 0.0000 | 0.0 | 기하 판정 함수 검증됨 |
| 행성간 궤적: 지구→화성 전이시간 | 전이시간(일) | 258.9 | 약 259 | 문헌값과 일치 |
| 행성간 궤적: SOI 반지름이 전이거리에서 차지하는 비율 | 지구 SOI 비율(%) | 1.18 | < 2.0 | Patched Conic 근사 타당 |
| 도킹/랑데부: 랑데부 델타-V 역산 후 목표 도달 | 최종 위치 오차(km) | 1.77m/s 소요 | < 1e-6km | 상태천이행렬 역산 검증됨 |
| 도킹/랑데부: CW 근사 오차(근접 vs 원거리) | 상대오차(%) | 0.0004% | 2.19%(500km에서) | 근접 근사임을 확인 |
| 궤도 결정: 관측값↔ECI 위치 왕복 변환 | 최대 왕복 오차(km) | 9.32e-12 | < 1e-6 | 04번 정변환의 정확한 역함수 |
| 궤도 결정: 작은 노이즈에서 원궤도 요소 복원 | 반장축/경사각 오차(%/도) | 0.0041% / 0.0043도 | < 1.0 각각 | 노이즈 낀 관측만으로 복원됨 |
| 궤도 결정: 관측 노이즈 대비 정확도 | 경사각 오차(작은 노이즈 vs 큰 노이즈, 도) | 0.0001도 vs 0.2117도 | 큰 노이즈 > 작은 노이즈 | 단조 증가 확인 |
| 저추력 전이: 나선 전이 후 반장축/이심률 | 반장축 오차(%) / 이심률 | 0.0063% / 0.002454 | < 1.0% / < 0.01 | 목표 반지름 도달, 원궤도 유지 |
| 저추력 전이: Edelbaum 근사 vs 수치적분 델타-V | 상대오차(%) | 0.14% | < 10.0% | 닫힌 형태 근사와 일치 |
| 자세동역학: 회전 에너지/각운동량 보존 | 최대 편차 | 4.29e-12 / 4.58e-12 | < 1e-6 | 토크 없는 자유회전, 수치오차 수준 보존 |
| 자세동역학: 중간축 정리(최소/최대축 vs 중간축) | 섭동 성장 배율 | 2.0배 / 1.4배 vs 100.0배 | 중간축 > 나머지 x5 | 중간축만 불안정함을 확인 |
| 궤도 유지: 주기적 보정 후 5년 방치 대비 오차 | 보정 없는 5년 오차(도) | 1.6927도 | 방치 시 계속 증가 | 보정 없이는 선형으로 계속 커짐 |
| 궤도 유지: 허용 오차에 따른 총 델타-V 편차 | 편차(%) | 10.0% | < 20% | 허용 오차 선택과 거의 무관 |
| 궤도 재진입: 가속 하강(전반부 대비 후반부 감소율) | 배율 | 2.20배 | 후반부 > 전반부 x2 | 지수함수 밀도 모델대로 재진입 근처에서 가속 |
| 궤도 재진입: 탄도계수/초기고도별 수명 범위 | 최소~최대(일) | 1.57〜15.74 / 1.19〜8.51 | 둘 다 단조 증가 | 탄도계수는 선형, 고도는 기하급수적 증가 |
| 라그랑주점: L1 위치(문헌값과 비교) | 지구 중심 기준 | 321,710km | 약 320,000~330,000km | 이분법으로 푼 5차 방정식 근이 문헌값과 일치 |
| 라그랑주점: L4/L1 선형 안정성(최대 고유값 실수부) | L4 / L1 | 0.000000 / 2.9321 | L4≈0, L1>0 | L4는 순허수(안정), L1은 양의 실수부(불안정) |
| 리아푸노프 궤도: L1 주기(선형 예측 대비 최대 편차) | 진폭 스윕 최대 편차 | 0.244% | 진폭 클수록 편차 증가 | 비선형 효과로 선형 예측에서 점점 벗어남 |
| 리아푸노프 궤도: L1/L2 주기 비율(같은 진폭) | L2/L1 | 1.2530 | 둘 다 수렴, 같은 자릿수 | 같은 방법이 L1/L2 양쪽에서 일반화됨 |
| PID 자세제어: point-and-hold 최종 오차(35도 초기 오차) | 3배 정착시간 후 오차 | 0.000도 | < 2도 | 정착시간/감쇠비 공식대로 매끄럽게 수렴 |
| PID 자세제어: 16번 중간축 텀블링 억제(최대 오차) | 제어없음 / 제어있음 | 180.0도 / 6.8도 | 5배 이상 억제 | PD 제어가 불안정한 텀블링을 거의 완전히 정지시킴 |
| 통합 임무: 궤도 전이 완료 후 자세 재정렬 최종 오차 | 정착시간 10.0초의 3배 적분 후 | 0.000도 | < 2도 | 핸드오프 상태에서 시작해도 동일한 수렴 기준을 만족 |
| 통합 임무: 총 임무시간(궤도전이+자세정렬) | 궤도 전이 0.84시간 + 자세 재정렬 | 0.85시간 | 두 단계 시간의 합과 일치 | 총 임무시간 = 궤도 전이시간 + 자세 재정렬 적분시간 |
| 3체 섭동: 달이 태양보다 큰 가속도 | 비율(달/태양) | 2.24배 | 약 2배(교과서값) | 질량은 작지만 거리의 세제곱 효과로 태양보다 섭동이 큼 |
| 3체 섭동: GEO에서 J2 대비 비율 | LEO 0.0129% / GEO 144.31% | GEO가 LEO보다 훨씬 큼 | 단조증가 | LEO에서는 무시 가능, GEO 근처에서는 J2를 능가 |
| SRP: 항력을 넘어서는 교차 고도 | 비율(SRP/항력) ≥ 1 | 800.0km | 약 550~800km | 9번 항력 모델과 직접 교차점이 존재 |
| SRP: 베타각 0도 원형 궤도의 식(eclipse) 비율 | 그림자 안 시간 비중 | 36.4% | 약 30~40%(실제 LEO 전형값) | 원통형 그림자 모델이 실측 LEO 식 비율과 일치 |
| 반작용휠: 외부 토크 없을 때 총 각운동량 크기 보존 | |L_total(t)| 최대 드리프트 | 5.1×10⁻¹⁶ | < 1×10⁻⁶ | 몸체-휠 반작용 부호가 맞아야만 성립하는 가장 엄격한 항등식 |
| 반작용휠: 포화 직후 디새추레이션으로 포인팅 회복 | 디새추레이션 후 오차 / 비용 비교 | 100.0초에 포화 → 0.007도로 회복 | 0 < 비용 < 호만 전이 델타-V | 대응이 10~20초만 늦어도 회복 실패 — 타이밍에 민감 |
| 자세 의존 SRP: 코사인 법칙 경계 케이스(0/60/90/180도) | 유효 A/m 대 max(0,cosθ)·max_A/m | 0도=100%, 60도=50%, 90도=0%, 180도=0% | 코사인 법칙과 정확히 일치 | 패널 뒷면(cosθ<0)은 항상 0으로 클램프 |
| 자세 의존 SRP: 텀블링 vs 태양지향 제어 | 유효 SRP 변동폭 / 태양지향 유지 비율 | 텀블링 31.8% 변동 / 태양지향 최소 99.91% 유지 | 텀블링은 뚜렷이 진동, 태양지향은 95% 이상 유지 | 캐논볼의 "항상 일정" 가정이 텀블링에서는 깨짐 |
| 가우스법: 이심률·근점편각 복원(14번이 못 다룬 요소) | 이심률 절대오차 / 근점편각 오차 | 0.0163 / 2.29도 | 이심률 < 0.05, 근점편각 < 5도 | 02번 범용 궤도요소 변환 재사용, 14번은 추정 자체가 불가능했던 요소 |
| 가우스법: 거의 공면인 관측 기하에서 ValueError | |D0|가 임계값 미만일 때 동작 | D0=-0.1547 (기본 시나리오) | 임계값 미만이면 조용히 틀린 값 대신 명시적 예외 | 가우스법의 근본적 기하 민감성을 감추지 않음 |
| GEO 궤도유지: 가우스 방정식이 07번 해석적 세차율을 재현 | 수치 평균 vs 해석적 공식 상대오차 | 0.0225% | 0.1% 미만 | 밑바닥부터 재구현한 공식이 이미 신뢰받는 07번 공식을 재현 |
| GEO 궤도유지: 3체 섭동이 J2 잔여를 지배 | 3체 섭동/J2 잔여 RAAN 변화율 비 | 970.9배 | 50배 이상(실측 약 970배) | GEO 궤도유지는 J2가 아니라 3체 섭동이 지배적 보정 대상 |