Monte Carlo & MEX 가속

N_MC = 1000;  % Monte Carlo 횟수

for i = 1:N_MC
    % 1. 파라미터 샘플링 (확률적)
    thrust_scale = 1.0 + (unc.thrust_3sig/3) * randn();
    Cd0 = nom.Cd0 * (1.0 + (unc.Cd0_3sig/3) * randn());
    
    % 2. 시뮬레이션 실행
    [apogee, max_speed, ...] = rocket_sim_core(thrust_scale, Cd0, ...);
    
    % 3. 결과 저장
    results.apogee(i) = apogee;
    results.max_speed(i) = max_speed;
end

% 4. 통계 분석
apogee_mean = mean(results.apogee);
apogee_std = std(results.apogee);

%% 불확실성 정의 (3-sigma 값)
unc.thrust_3sig = 0.15;       % 추력 ±15%
unc.Cd0_3sig = 0.20;          % Cd0 ±20%
unc.CN_alpha_3sig = 0.15;     % CN_alpha ±15%
unc.launch_angle_3sig = 2.0;  % 발사각 ±2°
unc.mass_3sig = 0.05;         % 질량 ±5%
unc.wind_max = 3.0;           % 최대 풍속 [m/s]

%% 공칭값
nom.Cd0 = 0.224;
nom.CN_alpha = 12.0;
nom.launch_angle = 86.0;

% randn()은 표준정규분포 (μ=0, σ=1)에서 샘플링
x = randn();  % -∞ ~ +∞ 범위의 값

% 원하는 분포로 변환: X = μ + σ * randn()
% 3σ 불확실성을 σ로 변환: σ = 3σ값 / 3

% 예: 추력 스케일 (공칭=1.0, 3σ=±15%)
sigma_thrust = unc.thrust_3sig / 3;  % 0.15/3 = 0.05
thrust_scale = 1.0 + sigma_thrust * randn();

% 예: Cd0 (공칭=0.224, 3σ=±20%)
Cd0 = nom.Cd0 * (1.0 + (unc.Cd0_3sig/3) * randn());

% rand()는 [0, 1) 범위의 균등분포

% 풍속: 0 ~ 3 m/s
wind_speed = unc.wind_max * rand();  % [0, 3)

% 풍향: 0 ~ 2π rad (모든 방향 동등)
wind_dir = 2*pi * rand();  % [0, 2π)

% 매번 다른 결과
rng('shuffle');

% 재현 가능한 결과
rng(42);  % 시드 = 42

% 바람 벡터 (관성계, NED)
wind_I = [wind_speed * cos(wind_dir); ...
          wind_speed * sin(wind_dir); ...
          0];

% Body frame으로 변환
wind_B = R_I2B * wind_I;

% 상대 속도 계산
V_rel = vel - wind_B;
V = norm(V_rel);

% 받음각 계산 (상대 속도 기준!)
alpha = atan2(V_rel(3), V_rel(1));
beta = atan2(V_rel(2), V_rel(1));

% 기본 통계
apogee_mean = mean(results.apogee);
apogee_std = std(results.apogee);

fprintf('평균: %.1f m\n', apogee_mean);
fprintf('표준편차: %.1f m\n', apogee_std);
fprintf('범위: %.1f ~ %.1f m\n', min(results.apogee), max(results.apogee));

% 3σ 범위 (99.7% 신뢰구간)
fprintf('3σ 범위: %.1f ~ %.1f m\n', ...
        apogee_mean - 3*apogee_std, ...
        apogee_mean + 3*apogee_std);

% 특정 값 이상 확률
percentile_340 = sum(results.apogee >= 339.57) / N_MC * 100;
fprintf('339.57m 이상 확률: %.1f%%\n', percentile_340);

% 경험적 CDF 계산
[f, x] = ecdf(results.apogee);

% 플롯
plot(x, f*100, 'b-', 'LineWidth', 2);
xlabel('최대 고도 [m]');
ylabel('누적 확률 [%]');
title('최대 고도 CDF');

histogram(results.apogee, 30, ...
          'FaceColor', [0.3 0.6 0.9], ...
          'EdgeColor', 'w');
hold on;

% 평균선
xline(apogee_mean, 'r-', 'LineWidth', 2);

% ±3σ 범위
xline(apogee_mean - 3*apogee_std, 'r--');
xline(apogee_mean + 3*apogee_std, 'r--');

% 비교값
xline(249, 'g-', 'OpenRocket');
xline(339.57, 'm-', '실측');

param_names = {'추력', 'Cd0', 'CN_alpha', '발사각', '질량', '풍속'};
param_data = [params.thrust_scale, params.Cd0, params.CN_alpha, ...
              params.launch_angle, params.mass_scale, params.wind_speed];

correlations = zeros(6, 1);
for i = 1:6
    r = corrcoef(param_data(:,i), results.apogee);
    correlations(i) = r(1,2);
    fprintf('  %s: %.3f\n', param_names{i}, correlations(i));
end

% 막대 그래프로 시각화
barh(correlations);
set(gca, 'YTickLabel', param_names);
xlabel('상관계수');
title('민감도 분석');

function [apogee, max_speed, ...] = rocket_sim_core(thrust_scale, ...)
%#codegen  % ← 이 주석이 Coder 호환성 검사 활성화

% Coder 호환: 배열 크기 명시적 초기화
pos = zeros(3, 1);
vel = zeros(3, 1);
quat = zeros(4, 1);

% int32로 타입 명시
N_steps = int32(floor(t_end / dt)) + 1;

%% 입력 타입 정의 (샘플 값)
thrust_scale = 1.0;
Cd0 = 0.15;
CN_alpha = 12.0;
launch_angle = 86.0;
mass_scale = 1.0;
wind_speed = 0.0;
wind_dir = 0.0;
Q_force = 1.0;
Q_moment = 0.01;
dt = 0.002;
t_end = 20.0;

%% codegen 설정
cfg = coder.config('mex');
cfg.GenerateReport = true;      % HTML 리포트 생성
cfg.ReportPotentialDifferences = false;

% 최적화 설정
cfg.SaturateOnIntegerOverflow = false;  % 오버플로 체크 해제
cfg.IntegrityChecks = false;             % 범위 체크 해제
cfg.ResponsivenessChecks = false;        % 응답성 체크 해제

%% 컴파일
input_args = {thrust_scale, Cd0, CN_alpha, launch_angle, mass_scale, ...
              wind_speed, wind_dir, Q_force, Q_moment, dt, t_end};

codegen('rocket_sim_core', '-config', cfg, '-args', input_args);

% MEX 버전 호출 (함수명에 _mex 접미사)
[apogee, time, speed, mach, impact, dx, dy] = ...
    rocket_sim_core_mex(1.0, 0.15, 12.0, 86.0, 1.0, ...
                        0.0, 0.0, 1.0, 0.01, 0.002, 20.0);

% MATLAB 버전과 동일한 인터페이스!

% ❌ 나쁜 예: 루프에서 배열 크기 증가
for i = 1:N
    results(i) = compute(i);  % 매번 메모리 재할당!
end

% ✅ 좋은 예: 미리 할당
results = zeros(N, 1);  % 한 번만 할당
for i = 1:N
    results(i) = compute(i);
end

% ❌ 나쁜 예: for 루프
for i = 1:N
    y(i) = x(i)^2 + 2*x(i) + 1;
end

% ✅ 좋은 예: 벡터화 (훨씬 빠름)
y = x.^2 + 2*x + 1;

function yi = interp1_coder(x, y, xi)
    % 선형 보간 (Coder 호환)
    n = length(x);
    
    % 경계 처리
    if xi <= x(1)
        yi = y(1);
        return;
    end
    if xi >= x(n)
        yi = y(n);
        return;
    end
    
    % 구간 찾기
    idx = 1;
    for k = 1:n-1
        if xi >= x(k) && xi < x(k+1)
            idx = k;
            break;
        end
    end
    
    % 선형 보간
    t = (xi - x(idx)) / (x(idx+1) - x(idx));
    yi = y(idx) + t * (y(idx+1) - y(idx));
end

function [rho, a] = atmos_isa(alt)
    % ISA 대기 모델 (Coder 호환)
    T0 = 288.15;      % 해수면 온도 [K]
    P0 = 101325;      % 해수면 압력 [Pa]
    L = 0.0065;       % 기온 감률 [K/m]
    R = 287.058;      % 기체 상수 [J/(kg·K)]
    g = 9.81;         % 중력 [m/s²]
    gamma = 1.4;      % 비열비
    
    T = T0 - L * alt;
    P = P0 * (T / T0)^(g / (R * L));
    rho = P / (R * T);
    a = sqrt(gamma * R * T);
end

target_apogee = 339.57;

% 가장 가까운 케이스 찾기
[~, best_idx] = min(abs(results.apogee - target_apogee));

% 시그마 편차 계산
thrust_sigma = (params.thrust_scale(best_idx) - 1.0) / (unc.thrust_3sig/3);
Cd0_sigma = (params.Cd0(best_idx) - nom.Cd0) / (nom.Cd0 * unc.Cd0_3sig/3);

fprintf('추력 스케일: %.3f (%+.2fσ)\n', params.thrust_scale(best_idx), thrust_sigma);
fprintf('Cd0: %.4f (%+.2fσ)\n', params.Cd0(best_idx), Cd0_sigma);