内弹道计算程序

%内弹道计算程序
function ndd
clc;clear;
%构造与装填条件
S =1.681; %枪（炮）膛横断面积 dm^2
q = 10.24; %弹重 kg
W_0 = 10.35; %药室容积 dm^3
l_g = 62.25; %身管行程 dm
%进程的进行条件
P_0 = 45000; %起动压力 kpa
theta = 0.25; %火药热力系数
K = 1.222; %次要功系数
%装填条件
delta =1.6; %火药密度 kg/dm^3
f = 948000; %火药力 kg*dm/kg
alpha =1; %余容 dm^3/kg
omega = 10.35; %药量 kg
V = 0.9627; %烧速指数
u =0.0000088; %烧速系数 dm^3/（s*kg）
%火药特征（仅适用于多孔火药）
e_1 = 0.0088; %厚度（肉厚）/2 dm
d_0 = 0.022; %孔道直径 dm
D_0 = 0.1364; %药粒直径（梅花形不用管此项） dm
c = 1.3; %药粒长度/2 dm
n = 7; %孔数
flag = 1; %1:圆柱形多孔火药 2:梅花形七孔火药 3:梅花形十四孔火药 4:梅花形十九孔火药
%形状函数
Flag = [ 1 n 0 D_0 0 0.2956;
2 8 65.2968 d_0+4*e_1 d_0+2*e_1 0.1547;
8/3 47/3 141.4764 d_0+4*e_1 d_0+2*e_1 0.1547;
3 21 195.8903 d_0+4*e_1 d_0+2*e_1 0.1547];
F1 = Flag（flag1）;
F2 = Flag（flag2）;
F3 = Flag（flag3）;
F4 = Flag（flag4）;
F5 = Flag（flag5）;
F6 = Flag（flag6）;
Pi_1 = （ F1*F4 + F2*d_0 ）/（2*c）;
Q_1 = （ F3*F5*F5 + F1*F4*F4 - F2*d_0^2 ）/（ 4*c^2 ）;
beta = e_1/c;
chi = （ Q_1 + 2*Pi_1 ）*beta/Q_1;
lambda = （ n - 1 - 2*Pi_1 ）*beta/（ Q_1 + 2*Pi_1 ）;
mu = （ 1 - n ）*beta^2/（ Q_1 + 2*Pi_1 ）;
psi_s = chi*（ 1 + lambda + mu ）;
rho = F6*（ d_0/2 + e_1 ）;
Z_s = 1 + rho/e_1;
chi_s = （ psi_s*Z_s^2 - 1 ）/（ Z_s^2 - Z_s ）;
lambda_s = psi_s/chi_s - 1;
%常数与初值计算
l_0 = W_0/S;
Delta = omega/W_0;
phi = K + omega/（3*q）;
v_j = 196*f*omega/（phi*theta*q）;
v_j = sqrt（v_j）;
B = 98*（e_1*S）^2/（ u*u*f*omega*phi*q ）;
B = B*（f*Delta）^（2-2*V）;
p_0 = P_0/（f*Delta）;
psi_0 = （1/Delta - 1/delta）/（f/P_0 + alpha - 1/delta）;
Z_0 = （sqrt（1+4*psi_0*lambda/chi） - 1）/（2*lambda）;
C = zeros（111）;
C（1）=c