clc;?clear;
%%%%%%%%%火炮構造及彈丸諸元
S=2.66e-3;???????????????%槍（炮）膛橫斷面積S??m^2
V0=1.51e-3;???????????????%藥室容積V0??m^3
l_g=3.624;?????????????%身管行程lg??m
m=2.8;?????????????????%彈重m??kg
%%%%%%%%裝藥條件
f=950e3;????????????????%火藥力f??J/kg?
omega=1.16;?????????????%裝藥量ω??kg
alpha=0.001;????????????%余容α??m^3/kg
rou_p=1600;?????????????%ρp???kg/m^3
theta?=0.25;????????????%火藥熱力系數
u1=5.127e-8;????????????%燃速系數??dm^3/（s*kg）
n=0.83;?????????????????%燃速指數?n
e1=5.5e-4;??????????????%裝藥藥厚e1?dm
d=5.5e-4;???????????????%?
chi=0.75;???????????????%裝藥形狀特征量?χ
lambda=0.12;????????????%裝藥形狀特征量?λ
mu=0.0;?????????????????%裝藥形狀特征量?μ
chi_s=1.696;????????????%裝藥分裂點形狀特征量?χs
lambda_s?=?-0.4104;?????%裝藥分裂點形狀特征量?λs
%%%%%%%%%%起始條件
p_0?=3e7;????????????????????%起動壓力??Pa
%%%%%%%%%%計算常數???????用于計算φ，本程序中已提供該參數，所以以下兩個參數未用到
phi1=1.02;???????????????????
lambda2=0;
%K=1.03;??????????????????????%運動阻力系數K
%%%%%%%%%%%%常量計算
%phi=phi1+lambda2*omega/m;
%phi=K+omega/3/m;??????????????????
phi=1.168;????????????????????????%次要功系數φ
Delta=omega/V0;???????????????????%裝填密度
l_0?=?V0/S;???????????????????????%
vj=sqrt（2*f*omega/（theta*phi*m））;?%
B?=?S^2*e1^2*（f*Delta）^（2-2*n）/（f*omega*phi*m*u1^2）;
lg1=l_g/l_0;
psi_s=chi*（1+lambda+mu）;
%es=（1+chi_s-sqrt（（1+chi_s）^2-4*chi_s*psi_s））/（2*chi_s）;
%rou=e1/es-e1;
%rou=0.2956*（d/2+e1）;
%Zk=（e1+rou）/e1;
Zk=（chi_s+sqrt（chi_s^2+4*psi_s*（1-chi_s）））/（2*psi_s）;??%%%%需確認的部分
rou=Zk*e1-e1;??????????????????????????????????????????%%%%%%%%%%應先求ρ，進而確定Zk
%%%%%%%%%%%%%%%%初值計算
psi_0=（1/Delta-1/rou_p）/（f/p_0