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

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

首先，开火方只把开火的信息传送给服务器，服务器把某架飞机开火的消息通知其它所有客户端，其它客户机根据开火飞机开火时的位置和朝向（其它飞机的姿态由服务器通知各客户端）画出子d轨迹（很可能此时画的轨迹并不计算物理特性，只是简单地按直线画图，或者，只进行简单的重力计算，并不会做碰撞计算），所以，所有计算机都能看到某架飞机开火，并看到d道，只是各计算机看到的d道并不完全一致，存在微小的误差。但这并不影响游戏效果。至此，解决了你所提出的“其它人不会观测到你的发射的d药”的问题。

另外，关于我说“没有必要去计算d道和通知客户机”，是对服务器而言的。因为，考虑到服务器的性能和负载能力，计算128架飞机（IL2好像可能128人连线）的d道对服务器要求太高，不太现实。这就是服务器没有必要计算d道的原因。关于“服务器没有必要把子d飞行轨迹通知各客户机”，是因为：1.这会造成大量的网络通信。2.其它客户机不会计算开火方是否命中（这会由开火方自行计算），它们只关心这次射击的结果（为了正确呈现图像，它们必需知道这个信息），并不关心子d飞行的详细过程。”

看来我是这里唯一认为服务器要计算d道轨迹的人，写的有点长，都看看吧，讨论一下

我尽量说的逻辑性和完整性以及前后呼应

我不同意你的观点，应为漏洞太多，要让程序运行是要符合逻辑的，你作了太多的假设，论证一个理论是否成立可以用假设和推断，但是如果有多个论证那么取信那个所用假设小少者，明白？

我总结一下我的理论

1 我没有否认过客户机要负责计算自己的d道，也就是客户机要计算自己的d道，并显示其路径。

2 服务器同样要做计算，但是有可能和客户机的不一样（延时，或其他原因造成），服务器计算的目的是为其他客户机显示，有人问为什么？因为如果只由客户机计算的话那d道数据完全依赖于客户机的专递，如果客户机网络不好，d道将无法表现，我问就会看到空中有直线飞行的d药，如果有人能观测到这个现象并录像我就改变观点。当前的sky可以说网络状况不完美，但为什么我没有发现这个现象？不能完整解释我的问题的，我将始终保持我的观点，如果给我解释也一定要有足够的逻辑性，不要用可能、大概、也许这样的字眼，我本人是程序员，要有完整的逻辑性，我知道程序中不认这样的语句，凡事都要有确定结果，除非是做抽奖程序才会用到随机算法。

3 服务器除了计算d道以外，还要负责数据传输的工作，如所有战机的位置数据，游戏规则，积分。。。，这点大家意见基本统一。

4 有人说服务器负责计算所有d道数据会占用大量带宽和计算量，这个不完全对，要看程序是怎么设计的，算法如何，如果做的糟糕就会出现这种现象。

首先计算量和传输量不成比例关系，我大个比方，每天的天气预报是个计算量繁重的任务吧需要亿次计算机才能完成，但是如果你是客户机查询天气，那么服务器会返回什么数据呢，是不是传输量并不很惊人？当然这是个比喻，实际情况要比这个复杂的。

再者有人说负责所有人的d道计算，就会造成负载过大，这个就是程序优化的问题，从当今cpu频率来说，这不是什么负担。有负担也不是什么不可完成的任务，记住算法，一个软件特别是游戏软件的命根子就是算法，无论什么引擎，图像质量，3d效果，都以算法为基础，如果不优化那我们每人有能力购置一套能够勉强大大游戏的机器，去用IBM服务器吧，我简单说一下优化方法，当然不会是IL2用的，只是给大家提个醒。

以下有点长请耐心的看看

d道是什么？d道是个矢量线段，对吧，它不会改变方向，它的抛物线很容易用2次曲线来描述，也就说数学模型比较容易，但是如果实时的不停计算这个函数也会给机器带来负担，所以设计游戏是要讲究效率和质量平衡的。在3维游戏里，为体现游戏场景要建立坐标系（分为左手和右手两种3D3用的是右手（好像）），又分为，世界坐标系，视觉坐标系，坐标分为X,Y,Z，在视觉坐标系中X,Y是你面到方向的平面坐标，Z就是景深，还有一个W值也用来表示景物的远近关系，W在矩阵转换中用到，关键的是Z轴，2D都不具备理解Z轴的能力，有个Z-BUFFER听说过吧3D卡的标志？既然是坐标系就有单位可以按比例表现真实距离，IL2中肯定是这样做的。还有一个最小单位用来表现粒度的，也就是刻画的真实性，我们不可能将现实的东西完全模拟出来。好说了一大堆该回到d道优化的问题上了，无论是机q还机炮的轨迹如果等分成小线段，都可以看成是直线，所以直到分段合理就可以乱真了，而且也没必要没个时刻都算，只要在线段的分界处算就行了，进一步优化干脆都没必要算，只要预先把分界点出的数值计算成表格，使用时从表里查询就行了。在一步优化，我们不用每个d药到计算，射击会打出一串d药，那我就一起算因为他们飞行的路径大致相同，我们只要预先确定一个散布值做每颗d药的分量即可，这样看来就算有上百人同时开火，我们的计算量大吗？

这还是我的一点小想法，IL2的不知道要精明多少倍，我认为服务器要计算的理由就是，除此以外，都有各种漏洞没有解决，服务器计算并通知但不做命中判断，只为了给除射击方以外的电脑显示用，这样做到最公平，子d一旦射出就不再受到双方网络影响，作弊？你只能看到你到的电脑中敌机不断中d，但是它照飞不误全无反应！

最后无论是那方计算d道，都存在延时，这是当今无法避免的，但是我劝某些人别想利用这个来个所谓提前躲避，一句话没用！因为延时情况是多变的，每一时刻都不同，在者自己和对方延时不一样时采取策略也不同，又时因该是滞后 *** 作才对，所以如果潜心研究这个倒不如把心思用在提高技术上，认同否？

---