전체 로직 설명서 지구 시뮬레이션 프로젝트 — 과학 원리 · 수학 모델 · 알고리즘 · 데이터 흐름

1. 프로젝트 전체 아키텍처

1-1. 5계층 구조

1-2. 파일-계층 매핑

● = 핵심 사용, △ = 부분 참조

2. Layer 1: 궤도 역학 계산 엔진

밀란코비치 궤도 요소를 시간의 함수로 계산하는 핵심 엔진입니다. Berger (1978)의 해석적 풀이를 간략화한 푸리에 급수를 사용합니다.

2-1. 이심률 (Eccentricity) 모델

과학적 원리

지구 궤도의 이심률은 주로 목성과 토성의 중력 섭동에 의해 변합니다. 완전한 계산에는 N-body 시뮬레이션이 필요하지만, 장주기 변동은 준주기적이어서 푸리에 급수로 근사할 수 있습니다.

수학 모델

e(t) = e₀ + Σᵢ Aᵢ · cos(2π·t/Tᵢ + φᵢ)

사용된 항:
  e₀ = 0.028011          // 평균 이심률

  항 1: A=0.010738, T=413.0 kyr, φ=3.135   // 413kyr 장주기 (가장 큰 진폭)
  항 2: A=0.008263, T=100.0 kyr, φ=0.208   // 100kyr 빙하기 주기
  항 3: A=0.004095, T= 95.0 kyr, φ=1.950   // 95kyr (100kyr과 간섭 → 맥놀이)
  항 4: A=0.003633, T=123.8 kyr, φ=4.885   // 보조 항
  항 5: A=0.002425, T=131.0 kyr, φ=2.105   // 보조 항

왜 이 항들인가?

코드 구현

function calcEccentricity(tKyr) {
  return 0.028011
    + 0.010738 * Math.cos(2*Math.PI * tKyr/413.0 + 3.135)
    + 0.008263 * Math.cos(2*Math.PI * tKyr/100.0 + 0.208)
    + 0.004095 * Math.cos(2*Math.PI * tKyr/95.0  + 1.950)
    + 0.003633 * Math.cos(2*Math.PI * tKyr/123.8 + 4.885)
    + 0.002425 * Math.cos(2*Math.PI * tKyr/131.0 + 2.105);
}

2-2. 황도경사각 (Obliquity) 모델

과학적 원리

자전축 기울기는 달과 태양의 조석력, 그리고 다른 행성들의 중력 섭동에 의해 변합니다. 달이 없었다면 기울기 변동은 현재의 2.4° 범위가 아니라 0°~85°까지 혼돈적으로 변했을 것입니다.

수학 모델

ε(t) = ε₀ + Σᵢ Bᵢ · cos(2π·t/Tᵢ + ψᵢ)

  ε₀ = 23.320556°        // 평균 황도경사각

  항 1: B=1.000°, T=41.0 kyr, ψ=2.399   // 주 주기 (지배적)
  항 2: B=0.430°, T=39.73 kyr, ψ=3.810  // 보조 (41kyr과 간섭)
  항 3: B=0.280°, T=53.6 kyr, ψ=1.280   // 장주기 변조
  항 4: B=0.175°, T=40.3 kyr, ψ=5.200   // 미세 조정

2-3. 기후적 세차 지수 (Climatic Precession Index)

과학적 원리

세차운동 자체는 ~26,000년 주기로 자전축이 원뿔을 그리는 것이지만, 기후에 영향을 미치는 것은 근일점 경도(ω̃)와 이심률의 곱인 e·sin(ω̃)입니다. 이것이 어느 반구의 어느 계절에 태양이 더 가까운지를 결정합니다.

수학 모델

P(t) = e(t) · sin(ω̃(t))

ω̃(t) = 2π·t/21.7 + 1.2
      + 0.6·sin(2π·t/23.7 + 0.5)   // 23kyr 성분
      + 0.3·sin(2π·t/18.97 + 2.1)  // 19kyr 성분

2-4. 일사량 (Insolation) 모델

과학적 원리

65°N 하지 일사량은 밀란코비치 이론에서 빙하기 발생의 핵심 지표입니다. 이 값이 임계치 이하로 내려가면 여름에 전년도 눈이 녹지 않아 빙하가 축적됩니다.

수학 모델

Q(t) = Q₀ + (∂Q/∂ε)·(ε(t) - ε̄) + (∂Q/∂P)·P(t)

  Q₀ = 475.0 W/m²        // 기준 일사량
  ∂Q/∂ε = +8.0 W/m²/°     // 경사각 민감도
  ∂Q/∂P = -70.0 W/m²       // 세차 민감도
  ε̄ = 23.32°              // 기준 경사각

전개:
  Q(t) = 475 + 8·(ε(t) - 23.32) - 70·e(t)·sin(ω̃(t))

3. Layer 2: 대기순환 모델

지구의 대기순환을 위도·경도·계절의 함수로 계산하는 모델입니다. 핵심 함수는 getWind(lat, lon, dayOfYear, obliquity)이며, 임의의 지점에서의 바람 벡터(u, v)를 반환합니다.

3-1. 3셀 순환 구조 로직

알고리즘 흐름

function getWind(lat, lon, dayOfYear, obliquity) {

  // 1단계: 태양 적위 계산
  const decl = obliquity * sin(2π * (dayOfYear - 80) / 365);

  // 2단계: ITCZ 위치 결정 (적위의 60%를 따라감)
  const itcz = decl * 0.6;

  // 3단계: ITCZ 기준 상대 위도 계산
  const relLat = lat - itcz;
  const absRelLat = |relLat|;

  // 4단계: 순환 셀 판별 → 기본 바람 계산
  if (absRelLat < 30°)  → 해들리 셀: 무역풍 (u < 0, v → ITCZ)
  else if (absRelLat < 60°) → 페렐 셀:  편서풍 (u > 0)
  else                     → 극순환 셀: 극동풍 (u < 0)

  // 5단계: 코리올리 편향 추가
  u += sign(lat) * strength * 0.25;

  // 6단계: 육지 가열 효과 (몬순)
  if (isLand && summer && 10° < |lat| < 40°)
    v += sign(lat) * 2.5;  // 해양→대륙 유입풍

  // 7단계: 제트기류 추가 (30° 및 60° 부근)
  u += gaussian(|lat|, 31°, σ=2°) * 4;
  u += gaussian(|lat|, 60°, σ=3°) * 3;

  // 8단계: 경도별 미세 변동
  u += sin(lon * π/60) * 0.8;

  return { u, v };
}

3-2. 각 단계의 물리적 근거

3-3. ITCZ 동역학과 황도경사각 연동

// 황도경사각이 ITCZ에 미치는 영향 체인:

ε 증가 → δ_max 증가 → ITCZ 이동폭 증가 → 몬순 강화
          ↓
        고위도 여름 일사량 ↑
          ↓
        적도-극 온도차 감소 (고위도 가열)
          ↓
        편서풍 약화, 해들리 셀 확장

ε 감소 → δ_max 감소 → ITCZ 이동폭 감소 → 몬순 약화
          ↓
        고위도 여름 일사량 ↓
          ↓
        적도-극 온도차 증가 (고위도 냉각)
          ↓
        편서풍 강화, 해들리 셀 축소

4. Layer 3: 지구 물리 데이터

4-1. 대륙 폴리곤 시스템

17개 대륙/섬 폴리곤을 [경도, 위도] 좌표 배열로 정의합니다. 아프리카, 유럽, 아시아, 북미, 남미, 호주, 남극, 그린란드, 인도, 인도네시아, 일본, 영국, 뉴질랜드, 마다가스카르, 보르네오, 파푸아뉴기니, 중앙아메리카를 포함합니다.

// 예: 아프리카 폴리곤 (20개 꼭짓점)
[[-17,15], [-17,21], [-13,28], [-5,36], [10,37], [11,33],
 [25,32], [33,30], [35,32], [42,13], [51,12], [50,2],
 [42,-2], [40,-11], [35,-26], [28,-34], [18,-35], [15,-28],
 [30,-18], [33,-12], ...]

4-2. 육지 마스크 생성 알고리즘

// 해상도: 720×360 (0.5° 격자) = 259,200 픽셀
// 빌드 시 1회 계산 → Uint8Array에 저장 → O(1) 조회

for each 격자점 (row, col):
  lat = 90 - row × 0.5
  lon = -180 + col × 0.5
  for each 대륙 폴리곤:
    if pointInPolygon(lon, lat, polygon):
      landMask[row × 720 + col] = 1
      break

Point-in-Polygon 알고리즘 (Ray Casting)

// 레이 캐스팅: 점에서 오른쪽으로 반직선을 쏘아
// 폴리곤 변과의 교차 횟수가 홀수면 내부

function pointInPolygon(x, y, poly):
  inside = false
  for 각 변 (i, j):
    if (yi > y) ≠ (yj > y):        // 변이 y를 걸치는가?
      xIntersect = (xj-xi)·(y-yi)/(yj-yi) + xi
      if x < xIntersect:
        inside = !inside              // 토글
  return inside

시간복잡도: O(n) per point, n = 꼭짓점 수
전체 마스크 빌드: O(259,200 × 17 × avg_vertices) ≈ ~0.5초

4-3. 바이옴 색상 시스템

function getTerrainColor(lat, lon, isLand):

  if |lat| > 75°     → 빙설 [220,230,240]
  if !isLand:
    depth = f(lon, lat)  → 심해↔천해 그라데이션
  if |lat| > 65°     → 툰드라 [180,200,190]
  if |lat| > 50°     → 북방림 [50,100,40]
  if |lat| < 10°     → 열대우림 [20,90,20]
  if |lat| < 30°:
    dryness = f(lon)   → 사막 [160,140,80] or 사바나 [80,130,40]
  else               → 온대 [34,110,34]

5. Layer 4: 렌더링 엔진

5-1. 3D 구면 투영 (Orthographic Projection)

정투영 변환 수학

// 위도경도 → 3D 직교좌표
x = cos(lat) · sin(lon - rotLon)
y = sin(lat)
z = cos(lat) · cos(lon - rotLon)

// 자전축 기울기 적용 (X축 회전)
y' = y·cos(tilt) - z·sin(tilt)
z' = y·sin(tilt) + z·cos(tilt)

// 화면 좌표 (z' > 0인 점만 가시)
screenX = cx + x · R
screenY = cy - y' · R
visible = (z' > 0)

역투영 (화면→위경도, 지구본 렌더링용)

// 각 화면 픽셀에서 위경도를 역산
x  = (screenX - cx) / R
y' = -(screenY - cy) / R
z' = √(1 - x² - y'²)    // 구면 위의 점 (존재 조건: x²+y'² ≤ 1)

// 기울기 역변환
y = y'·cos(tilt) + z'·sin(tilt)
z = -y'·sin(tilt) + z'·cos(tilt)

lat = arcsin(y) · 180/π
lon = arctan2(x, z) · 180/π + rotLon

렌더링 파이프라인

5-2. 바람 입자 시스템

// 파티클 풀: 3,000개
// 각 파티클 속성:
{
  lat, lon,           // 현재 위치
  prevLat, prevLon,   // 이전 위치 (궤적 선분용)
  age, maxAge         // 수명 관리 (60~120 프레임)
}

// 업데이트 루프 (매 프레임):
for each particle:
  1. 이전 위치 저장
  2. getWind(lat, lon, day, obliquity) → {u, v}
  3. Δlon = u × dt × 0.15 / cos(lat)   // 메르카토르 보정
     Δlat = v × dt × 0.15
  4. 위치 업데이트 + 경도 래핑 (-180~180)
  5. age++ → maxAge 초과 시 랜덤 위치로 리셋

// 렌더링:
  - 지도: prevPos → curPos 선분 (풍속별 색상)
  - 지구본: 3D 투영 후 가시 반구만 그리기
  - 수명 기반 알파 페이드 (등장: 0→1, 퇴장: 1→0)
  - 풍속 색상 매핑: <4m/s 파랑 → 4~7 노랑 → 7~10 주황 → >10 빨강

5-3. 시계열 차트 렌더링

// 사전 계산: N=2400~3600 데이터 포인트
// 시간 범위: -500kyr~+100kyr (또는 -800kyr~+100kyr)

function drawChart(canvas, data, color, yMin, yMax):
  1. 그리드 라인 + Y축 라벨 (4등분)
  2. X축 라벨 (100kyr 간격)
  3. 간빙기 마커 (MIS 5e, 7, 9, 11, ...) 수직선
  4. 데이터 선 그리기 (Path2D)
  5. 면적 그라데이션 채우기 (color → transparent)
  6. 현재 시점 수직 점선 + 값 표시 점

6. Layer 5: 분석 및 비교 엔진

6-1. 분석 흐름

6-2. 고기후 데이터 포인트 (참조값)

7. 모듈 간 데이터 흐름

7-1. 실시간 시뮬레이션 프레임 루프 (earth-simulation.html)

requestAnimationFrame(animate)

function animate(timestamp):
  dt = timestamp - lastTime      // 프레임 간격 (ms)

  // ① 상태 업데이트
  orbitalAngle += dt × orbSpeed  // 공전 각도
  rotLon += dt × rotSpeed        // 자전 경도
  dayOfYear = f(orbitalAngle)    // 궤도각 → 날짜

  // ② 바람 입자 업데이트
  for each particle:
    wind = getWind(p.lat, p.lon, dayOfYear, obliquity)
    p.lon += wind.u × dt
    p.lat += wind.v × dt

  // ③ 렌더링 (3개 캔버스 동시)
  renderOrbital(ctx1, orbitalAngle, ecc, obliquity)
  renderGlobe(ctx2, rotLon, obliquity, dayOfYear)
  renderMap(ctx3, dayOfYear, obliquity)

  // ④ UI 업데이트
  displayValues(dayOfYear, obliquity, ecc, wind)

7-2. 전체 데이터 의존성 그래프

8. 핵심 알고리즘 상세 해설

8-1. 구면 위 바람 입자 이동의 메르카토르 보정

// 문제: 위도가 높아질수록 경선 간 실제 거리가 줄어듦
// cos(60°) = 0.5 → 경도 1°의 실거리가 적도의 절반
// 보정 없이 동일한 Δlon을 적용하면 고위도에서 비정상적 가속

Δlon = u × dt / max(cos(lat × π/180), 0.1)
                                         ↑ 극점 발산 방지
Δlat = v × dt
// 이 보정으로 고위도에서의 동서 이동이 실제 거리에 비례

8-2. 지구본의 조명 모델

// Lambert 반사 모델 (확산 반사)

// 태양 방향 벡터 (계절에 따라 변화)
sunDir = normalize(sin(seasonAngle)×0.3, -0.2, 1.0)

// 각 표면점의 법선 벡터 = (x, -y', z') (구면이므로 위치=법선)
light = dot(normal, sunDir) = nx·sx + ny·sy + nz·sz

// 최소 조도 보장 (암부도 완전히 검지 않게)
light = max(0.15, light)

// 대기 산란 효과 (가장자리에서 파란 빛)
edge = 1 - distance_from_center / R
atmosphere = (1 - edge)³ × 0.4
finalColor.b += atmosphere × 180

8-3. 기온 프록시 추정 (trend-analysis.html)

// 빙하코어 δ¹⁸O와 유사한 패턴을 만들기 위한 간략 모델
// 핵심: 일사량의 시간 적분 (빙하 볼륨은 일사량 적분에 비례)

function calcTemp(tKyr):
  Q = calcInsol(tKyr)

  // 빙하 볼륨 프록시: 과거 8kyr 일사량 부족분의 적분
  ice = Σ(t-8 to t) [475 - calcInsol(τ)] × 0.12

  // 기온 = 일사량 직접효과 + 빙하 피드백 + 경사각 효과
  T = (Q - 475) × 0.05      // 일사량 직접
    - ice × 0.3              // 빙하-알베도 양의 피드백
    + (ε - 23.3) × 0.4       // 경사각 고위도 효과

9. 과학적 근거 및 한계

9-1. 사용된 과학적 모델의 근거

9-2. 이 시뮬레이션의 한계

9-3. 더 정밀한 시뮬레이션을 위해 필요한 것