using TyBase#使用fopen等函数需要用到此库
using TyPlot#使用plot函数需要用到此库
using TyMath#使用math函数需要用到此库

定义梯形机翼参数

root_chord=2 # 根弦长度
tip_chord=0.8 # 梢弦长度
half_span=3 # 半展长
alpha=5 # 单位（度）
half_chord_sweep_angle=20 # 二分之一弦线后掠角（角度制）
half_chord_sweep_rad=deg2rad(half_chord_sweep_angle) # 转换为弧度

网格参数

num_chordwise=10 # 弦向网格数量
num_spanwise=10 # 半展向网格数量
total_grids = num_spanwise * num_chordwise * 2 #双侧网格总数

创建一个结构体，用于存储单个网格的信息

mutable struct Grids
x1
z1
x2
z2
cp_x
cp_z
end
grids = repelem(Grids(nothing, nothing, nothing,nothing,nothing,nothing),total_grids,1)

生成单侧展向坐标（右半部分）

z_right=LinRange(0, half_span, num_spanwise+1)

计算右侧前缘和后缘坐标

front_x_right=zeros(num_spanwise+1,1)
back_x_right=zeros(num_spanwise+1,1)
for i=1:num_spanwise + 1
# 计算根弦长度与梢弦长度的差值
chord_diff=root_chord-tip_chord
# 计算当前展向位置的弦长
current_chord=root_chord-chord_diff*(z_right[i]/half_span)
# 计算当前展向位置二分之一弦线的z坐标
half_chord_z=z_right[i]*tan(half_chord_sweep_rad)
# 计算前缘x坐标，反转方向
front_x_right[i]=-(half_chord_z+current_chord/2)
# 计算后缘x坐标，反转方向
back_x_right[i]=-(half_chord_z-current_chord/2)
end

生成右侧网格坐标

X_right=zeros(num_spanwise+1, num_chordwise+1)
Z_right=zeros(num_spanwise+1, num_chordwise+1)
for i=1:num_spanwise+1
X_right[i, :]=LinRange(front_x_right[i], back_x_right[i],num_chordwise+1)
Z_right[i, :]=repelem(z_right[i],1,num_chordwise+1)
end

生成左侧坐标（对称处理）

X_left=X_right
Z_left=-Z_right

合并左右两侧坐标

X=[X_left; X_right]
Z=[Z_left; Z_right]

绘制图形

figure()
hold("on")

绘制水平网格线（展向，与z轴平行）

for i=1:num_spanwise+1
plot(Z_left[i, :], X_left[i, :], "k-")
plot(Z_right[i, :], X_right[i, :], "k-")
end

绘制垂直网格线（弦向，与x轴平行）

for j=1:num_chordwise+1
plot([Z_left[:, j]; Z_right[:, j]], [X_left[:, j]; X_right[:, j]], "k-")
end

存储网格信息并计算影响系数

grid_num = 0 # 网格计数器
for side = 1:2 # 1=左侧，2=右侧
if side == 1
local X_side = X_left
local Z_side = Z_left
else
local X_side = X_right
local Z_side = Z_right
end

for i = 1:num_spanwise
    for j = 1:num_chordwise
    global grid_num
    grid_num +=1
    # 获取当前网格的四个顶点坐标
    x1 = X_side[i, j]
    z1 = Z_side[i, j]
    x2 = X_side[i, j+1]
    z2 = Z_side[i, j+1]
    x3 = X_side[i+1, j+1]
    z3 = Z_side[i+1, j+1]
    x4 = X_side[i+1, j]
    z4 = Z_side[i+1, j]
    
    # 计算 1/4 弦线和 3/4 弦线的坐标
    x_quarter_chord_1 = x1 .+ 0.75*(x2 - x1)
    z_quarter_chord_1 = z1 .+ 0.75*(z2 - z1)
    x_quarter_chord_2 = x4 .+ 0.75*(x3 - x4)
    z_quarter_chord_2 = z4 .+ 0.75*(z3 - z4)
    
    x_three_quarter_chord_1 = x1 .+ 0.25*(x2 - x1)
    z_three_quarter_chord_1 = z1 .+ 0.25*(z2 - z1)
    x_three_quarter_chord_2 = x4 .+ 0.25*(x3 - x4)
    z_three_quarter_chord_2 = z4 .+ 0.25*(z3 - z4)
    
    # 绘制附着涡线（1/4 弦线）
    plot([z_quarter_chord_1, z_quarter_chord_2], [x_quarter_chord_1, x_quarter_chord_2], "r--")
    
    # 计算控制点坐标（3/4 弦线中点）
    control_point_z = (z_three_quarter_chord_1 + z_three_quarter_chord_2) / 2
    control_point_x = (x_three_quarter_chord_1 + x_three_quarter_chord_2) / 2
    
    # 绘制控制点
    plot(control_point_z, control_point_x, "b.")
    # 存储网格信息
    grids[grid_num].x1 =x_quarter_chord_1 
    grids[grid_num].z1 =z_quarter_chord_1
    grids[grid_num].x2 =x_quarter_chord_2
    grids[grid_num].z2 = z_quarter_chord_2
    grids[grid_num].cp_x = control_point_x
    grids[grid_num].cp_z = control_point_z
    end
end

end

#标注坐标轴
xlabel("Z 坐标")
ylabel("X 坐标")
title("平板机翼二维网格划分")

计算所有影响系数 C_ij

num_grids = grid_num # 总网格数（双侧）
C = zeros(num_grids, num_grids)
for i = 1:num_grids
# 当前控制点坐标
xi = grids[i].cp_x
zi = grids[i].cp_z

for j = 1:num_grids
    # 马蹄涡j的四个顶点坐标
    x1j = grids[j].x1
    z1j = grids[j].z1
    x2j = grids[j].x2
    z2j = grids[j].z2

    #计算影响系数 Cij

    part0=-1 / ((xi - x1j) * (z2j - z1j) - (x2j - x1j) * (zi - z1j))

    part1_num = (x2j - x1j) * (xi - x1j) + (z2j - z1j) * (zi - z1j)
    part1_den = sqrt((xi - x1j)^2 + (zi - z1j)^2)
    part1 = part1_num / part1_den
    
    part2_num = (x2j - x1j) * (xi - x2j) + (z2j - z1j) * (zi - z2j)
    part2_den = sqrt((xi - x2j)^2 + (zi - z2j)^2)
    part2 = part2_num / part2_den

    term2 = part1 - part2

    term3_part1 = (1 + (xi - x1j) / sqrt((xi - x1j)^2 + (zi - z1j)^2)) / (z1j - zi)
    term3_part2 = (1 + (xi - x2j) / sqrt((xi - x2j)^2 + (zi - z2j)^2)) / (z2j - zi)

    term3 = term3_part1 - term3_part2

    C[i,j]=(1/4*pi*100)*(part0 * term2 + term3)
end

end
C
为什么运行出来Cij都长一样，哪里出了问题