% ===================== 局部函数：面积-马赫数关系反解 =====================
function M = solveAreaMach(areaRatio, gamma, branch)

% 方程：
% (A/A*)^2 = (1/M^2)*[ 2/(gamma+1)*(1+(gamma-1)/2*M^2) ]^((gamma+1)/(gamma-1))

% 由面积比 A/A* 反解马赫数 M

% branch    : 'subsonic' 或 'supersonic'
%


    f = @(M) (1./M.^2) .* ...
        ((2/(gamma+1)) * (1 + (gamma-1)/2 * M.^2)).^((gamma+1)/(gamma-1)) ...
        - areaRatio^2;

    switch lower(branch)
        case 'subsonic'
            % 亚声支：0 < M < 1
            M = fzero(f, [1e-6, 0.999]);
        case 'supersonic'
            % 超声支：M > 1
            M = fzero(f, [1.001, 20]);
        otherwise
            error('branch 必须为 subsonic 或 supersonic');
    end
end