旋翼翼型参数化拟合

function [upper_coeff, lower_coeff] = fit_oa309_airfoil()
    % 读取翼型坐标数据
    data = readmatrix('oa309.dat.txt', 'NumHeaderLines', 2);
    x = data(:, 1);
    y = data(:, 2);
    
    % 分离上表面和下表面坐标
    [x_upper, y_upper, x_lower, y_lower] = separate_airfoil_surfaces(x, y);
    
    % 拟合上表面系数
    upper_coeff = fit_cst_surface(x_upper, y_upper);
    
    % 拟合下表面系数
    lower_coeff = fit_cst_surface(x_lower, y_lower);
    
    % 验证拟合结果
    visualize_fit(x_upper, y_upper, x_lower, y_lower, upper_coeff, lower_coeff);
end

function [x_upper, y_upper, x_lower, y_lower] = separate_airfoil_surfaces(x, y)
    % 分离上表面和下表面坐标
    idx = find(y < 0, 1);  % 找到第一个负y值点
    x_upper = x(1:idx-1);
    y_upper = y(1:idx-1);
    x_lower = x(idx:end);
    y_lower = y(idx:end);
    
    % 添加前缘点(0,0)到两个表面
    x_upper = [0; x_upper];
    y_upper = [0; y_upper];
    x_lower = [0; x_lower];
    y_lower = [0; y_lower];
end

function coeff = fit_cst_surface(x, y)
    % 使用CST参数化方法拟合翼型表面
    n = 7;  % CST阶次
    binoms = arrayfun(@(i) nchoosek(n, i), 0:n);  % 二项式系数
    
    % 构建设计矩阵
    A = zeros(length(x), n+2);  % 包括A0-A7和y_tail
    
    for k = 1:length(x)
        xk = x(k);
        term = xk^0.5 * (1-xk)^1.0;  % 前缘形状函数
        
        % 计算CST基函数
        for i = 0:n
            A(k, i+1) = term * binoms(i+1) * (1-xk)^(n-i) * xk^i;
        end
        
        % 添加y_tail项
        A(k, n+2) = xk;  % x*y_tail项
    end
    
    % 求解最小二乘问题
    coeff = A \ y;
end

function visualize_fit(x_upper, y_upper, x_lower, y_lower, upper_coeff, lower_coeff)
    % 计算拟合曲线
    y_upper_fit = evaluate_cst(x_upper, upper_coeff);
    y_lower_fit = evaluate_cst(x_lower, lower_coeff);
    
    % 计算拟合误差
    upper_error = max(abs(y_upper - y_upper_fit));
    lower_error = max(abs(y_lower - y_lower_fit));
    
    fprintf('上表面最大拟合误差: %.6f\n', upper_error);
    fprintf('下表面最大拟合误差: %.6f\n', lower_error);
    
    % 绘制结果
    figure;
    hold on;
    grid on;
    axis equal;
    
    % 绘制原始数据
    plot(x_upper, y_upper, 'ro', 'MarkerSize', 4, 'DisplayName', '上表面原始数据');
    plot(x_lower, y_lower, 'bo', 'MarkerSize', 4, 'DisplayName', '下表面原始数据');
    
    % 绘制拟合曲线
    plot(x_upper, y_upper_fit, 'r-', 'LineWidth', 2, 'DisplayName', '上表面拟合');
    plot(x_lower, y_lower_fit, 'b-', 'LineWidth', 2, 'DisplayName', '下表面拟合');
    
    % 标注关键参数
    text(0.05, 0.1, sprintf('前缘半径: %.5f', 0.00277), 'FontSize', 10);
    text(0.05, 0.08, sprintf('最大厚度: 12%% @ 41%%弦长'), 'FontSize', 10);
    text(0.05, 0.06, sprintf('最大弯度: 1.22%% @ 18.6%%弦长'), 'FontSize', 10);
    
    title('OA309翼型CST参数化拟合');
    xlabel('弦向位置 (x/c)');
    ylabel('法向位置 (y/c)');
    legend('Location', 'best');
    hold off;
end

function y_fit = evaluate_cst(x, coeff)
    % 使用CST系数计算翼型坐标
    n = 7;
    binoms = arrayfun(@(i) nchoosek(n, i), 0:n);
    
    y_fit = zeros(size(x));
    
    for k = 1:length(x)
        xk = x(k);
        term = xk^0.5 * (1-xk)^1.0;
        sum_cst = 0;
        
        % 计算CST求和项
        for i = 0:n
            sum_cst = sum_cst + coeff(i+1) * binoms(i+1) * (1-xk)^(n-i) * xk^i;
        end
        
        % 添加y_tail项
        y_fit(k) = term * sum_cst + xk * coeff(n+2);
    end
end