matlab 内d道计算问题

原因

所给链接的代码不完整，缺少ndd_fun函数。

顺便鄙视一下该链接的提供者，这么广为流传的东西下载居然还要积分，简直穷疯了。

代码

帮你好好找了一下，找到了完整的程序，供参考（全部代码保存到一个M文件运行即可，或直接下载附件）：

function ndd

%59nian130

A=0.87         %q(炮)膛横断面积A  dm^2

G=19%33.4           %d重  kg

W0=2.04        %药室容积  dm^3

l_g=25.0       %身管行程  dm

P_0 =30000       %起动压力  kpa

fai1=1.02        %次要功系数

K=1.03        %运动阻力系数φ1

theta =0.2       %火药热力系数

%=========================================

f=950000          %火药力  kg*dm/kg 

alpha=1         %余容  dm^3/kg

delta=1.6         %火药重度γ

%==================================

ome=2.2%12.9      %第一种装药量  kg

u1=5.0024*10^-5        %第一种装药烧速系数  dm^3/(s*kg)

n1=0.82         %第一种装药的压力指数n1

lambda=-0.0071     %第一种装药形状特征量λ1

lambda_s=0     %第一种装药分裂点形状特征量λ1s

chi=1.00716          %第一种装药形状特征量χ1

chi_s=0          %第一种装药分裂点形状特征量χ1s

mu=0            %第一种装药形状特征量μ1

et1=1.14*10^-2            %第一种装药药厚δ01

d1=2.5*10^-2             %第一种装药火药内径d1

Ro1=0             %药型系数α1

%=========================================

%常数与初值计算-----------------------------------------------------

l_0=W0/A

Delta=ome/W0

phi=K + ome/(3*G)

v_j=196*f*ome/(phi*theta*G)

v_j=sqrt(v_j)

B = 98*(et1*A)^2/( u1*u1*f*ome*phi*G )

B=B*(f*Delta)^(2-2*n1)

Z_s=1+Ro1*(d1/2+et1)/et1

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(1,12)

C(1)=chiC(2)=lambdaC(3)=lambda_sC(4)=chi_sC(5)=Z_s%

C(6)=thetaC(7)=BC(8)=n1C(9)=DeltaC(10)=deltaC(11)=alphaC(12)=mu

C

y0=[Z_000psi_0]

options = odeset('outputfcn','odeplot')

[tt,y] = ode45(@ndd_fun,0:100,[Z_000],options,C)

l = y(:,2)

l = l*l_0

fl = find(l>=l_g)

fl = min(fl)

[tt,y] = ode45(@ndd_fun,0:0.005:fl,[Z_000],options,C)

Z = y(:,1)lx = y(:,2) vx = y(:,3) 

psi = (Z>=0&Z<1).*( chi*Z.*(1 + lambda*Z + mu*Z) ) +...%%%%%%%%%

      (Z>=1&Z<Z_s).*( chi_s*Z.*(1 + lambda_s*Z) ) +...

      (Z>=Z_s)*1

l_psi = 1 - (Delta/delta)*(1-psi) - alpha*Delta*psi

px = ( psi - vx.*vx )./( lx + l_psi )

p = px*f*Delta/100

v = vx*v_j/10

l = lx*l_0

t = tt*l_0*1000/v_j

fl = find(l>=l_g)

fl = min(fl)+1

p(fl:end)=[]v(fl:end)=[]l(fl:end)=[]t(fl:end)=[]

pd=px*f*Delta/100/(1+ome/3/fai1/G)

pt=pd*(1+ome/2/fai1/G)

aa=max(px)

M=find(px==aa)

Pm=[tt(M)*l_0*1000/v_j lx(M)*l_0 vx(M)*v_j/10 px(M)*f*Delta/100 pt(M) pd(M) psi(M) Z(M)]

%ll=length(tt)

ran=find(Z>=1)

ran=min(ran)

Zf=[tt(ran)*l_0*1000/v_j lx(ran)*l_0 vx(ran)*v_j/10 px(ran)*f*Delta/100 pt(ran) pd(ran) psi(ran) Z(ran)]

jie=find(psi>=1)

jie=min(jie)

psij=[tt(jie)*l_0*1000/v_j lx(jie)*l_0 vx(jie)*v_j/10 px(jie)*f*Delta/100 pt(jie) pd(jie) psi(jie) Z(jie)]

pg=[tt(end)*l_0*1000/v_j lx(end)*l_0 vx(end)*v_j/10 px(end)*f*Delta/100 pt(end) pd(end) psi(end) Z(end)]

Ry1=[ZfpsijpgPm]

Ry2=[tt*l_0*1000/v_j lx*l_0 vx*v_j/10 px*f*Delta/100 pt pd psi Z]

subplot(2,2,1)

plot(t,p,'linewidth',2)

grid on

xlabel('\fontsize{8}\bft  (ms)')

ylabel('\fontsize{8}\bfp  (kg/cm^{2})')

title('\fontsize{8}\bft-p曲线')

subplot(2,2,2)

plot(t,v,'linewidth',2)

grid on

xlabel('\fontsize{8}\bft  (ms)')

ylabel('\fontsize{8}\bfv  (m/s)')

title('\fontsize{8}\bft-v曲线')

subplot(2,2,3)

plot(l,p,'linewidth',2)

grid on

xlabel('\fontsize{8}\bfl  (dm)')

ylabel('\fontsize{8}\bfp  (kg/cm^{2})')

title('\fontsize{8}\bfl-p曲线')

subplot(2,2,4)

plot(l,v,'linewidth',2)

grid on

xlabel('\fontsize{8}\bfl  (dm)')

ylabel('\fontsize{8}\bfv  (m/s)')

title('\fontsize{8}\bfl-v曲线')

tspan = length(t)/20

tspan = 1:ceil(tspan):length(t)

tspan(end) = length(t)

fprintf('        t(ms)     p(kg/cm^2)     v(m/s)       l(dm)')

format short g

Result = [t(tspan) p(tspan) v(tspan) l(tspan)]

format

%------------------------------------------------------------------

function dy = ndd_fun(t,y,C)

chi=C(1)lambda=C(2)lambda_s=C(3)chi_s=C(4)Z_s=C(5)mu=C(12)

theta=C(6)B=C(7)V=C(8)Delta=C(9)delta=C(10)alpha=C(11)

Z = y(1) l = y(2) v = y(3)

psi = (Z>=0&Z<1).*( chi*Z.*(1 + lambda*Z + mu*Z) ) +...

(Z>=1&Z<Z_s).*( chi_s*Z.*(1 + lambda_s*Z) ) +...

(Z>=Z_s)*1

l_psi = 1 - (Delta/delta)*(1-psi) - alpha*Delta*psi

p = ( psi - v*v )/( l + l_psi )

dy(1) = sqrt(theta/(2*B))*(p^V)*(Z>=0&Z<=Z_s)

dy(2) = v

dy(3) = theta*p/2

dy = [dy(1)dy(2)dy(3)]

结果

并未听说有独立销售的d道计算程序，也许是公安部不允许吧，毕竟这种软件流入不法分子手中后果是相当...的，但是5.11品牌的HRT系列手表集成了d道计算功能，输入风速、风向、d种等数据后手表会计算出相应的d道并显示出准星和照门的具体调整方式，不过手表目前不支持中文，而且价格也较高，会在1000元以上。

其实你可以手动计算，有几个公式（不全是公式）先介绍一下：

1.风的分类，你可以拿出指针式手表看一下，12点和6点这条线上的风叫零速风，9点和3点这一方向上的风叫全速风，其余时刻对应方向上的风为半速风。

2.风速的测量。当然，如果有专业的测量工具更好，无工具情况下，最常用的方法是目测旗帜等物品被风吹起的角度（无旗帜等物品时，可站定手持纸或其他类似物品，举到与头同高处，松手，然后测量物品掉落处与自己的距离，最后利用三角函数计算出角度）这个角度处以常数4即为风速，但须注意，此方法测出的风速单位为MPh而不是常用的KMh

3.将风速转化为分钟角度。1分钟角度=1/60度，也就是狙击镜刻度的1/20，差不多就是每100米1寸的偏差，如果你是狙击手，利用分钟角度去调节瞄准角度以进行高精度射击是必要的，公式如下（射程m/100*风速里每小时)/常数=分钟角度（怕你看不懂这种写法，口述一下：射程与100倍的风速的商除以常数等于分钟角度 射程单位为米、风速单位为里/小时），这里的常数是不固定的，100~500米常数为15 600米常数为14 700~800米常数为13 900米常数为12 1000米常数为11，例如射程700米 风速为10里 计算出的分钟角度为5.38 需要注意此算法是全速风的算法，如是半速风，只要将5.38除以2就行

以上只是d道计算的方法，实际应用中，影响子d飞行的因素还有很多，比如你的呼吸、你对扳机的控制LOS（视线）与BP（d道）的关系、所处的高度和湿度等，1000米以上的射击还要考虑地球自转对子d的影响，另外子d的品牌不同也会对d道造成影响，总之，d道学是很大的一个学科，要联系到很多知识，做到灵活运用的惟一方式就是多练习，所以如果你真的是部队的新Gunner，那么还是建议你多练习手动计算以积累经验，不要去想软件解决，毕竟软件那种线性的统计不能规划多学科融合的d道学的方方面面；如果你和我一样，只是个q械发烧友的话，那么我可以负责的讲，手动计算有着软件计算所难以企及的乐趣。

很抱歉你要的软件目前真的没有零售（也许会有内部使用的，那可没几人能弄来），你要是有1000左右的闲钱，可以考虑弄快5.11的HRT军表玩玩。

我还是觉得手动计算好些，不过不管怎样，只要自己喜欢就好不是？

呵呵，还是祝你在这个方面找到自己的乐趣所在，并乐于其中吧

通常用数值方法来求解导d的运动方程。d道计算在导d各设计阶段和导d试验以及射表编制工作中，是一项必不可少的工作。由于导d射程的增大和导d命中精度的提高，对d道计算的精度要求也愈来愈高。d道计算的精度，依赖于导d运动方程的描述精确度和数值计算方法。对有控制的导d来说，d道计算影响命中精度的主要因素，是制导特征量的计算精度和地球非球形摄动。

当导d射程达到1000公里左右时，就必须考虑地球扁率的摄动，而对洲际导d还必须顾及高阶项的摄动（如重力异常等），计算结果表明，高阶摄动对射程影响为1～2公里的量级。随着导d命中精度的提高，对这些因素的考虑必须愈来愈精确（如果导d采用末制导，情况则不同，这时命中精度只依赖于制导的精度）。d道计算与高速电子计算技术结合，使飞行仿真模拟成为现实。给定d道飞行模型，编制计算程序，通过计算机大量计算，可以得到各种随机过程的许多结果，例如导d最大射程的概率特征、落点散布等。用飞行仿真模拟可以部分代替昂贵的飞行试验。

---