哪位高手有matlabd道仿真程序

哪位高手有matlabd道仿真程序,第1张

第一个例子

导d仿真Matlab源代码(missile simulation)

http://www.matlabsky.com/thread-284-1-1.html

第二例子

%by dynamic

%see also http://www.matlabsky.com

%2009.2.18

%比例导引理想d道胡改仿真凳做备Matlab代码

%给出了比例导引法的差分方程, 采用Matlab语言, 对比例导引法的导d理想d枣毁道进行了三维数据仿真, 绘制出三维理想d道。随导d与目标参数变化情况, 对理想d道的特性进行了分析。

clear all

clc

clear

tt=0.1

sm=0.6*tt

st=0.42*tt

x(1)=0y(1)=0z(1)=0

pmr(:,1)=[x(1)y(1)z(1)]

ptr(:,1)=[25510]

m=3

q(1)=0

o(1)=0

a(1)=0

for(k=2:600)

ptr(:,k)=[25-0.42*cos(pi/6)*tt*k510+0.42*sin(pi/6)*k*tt]

r(k-1)=sqrt((ptr(1,k-1)-pmr(1,k-1))^2+(ptr(2,k-1)-pmr(2,k-1))^2+(ptr(3,k-1)-pmr(3,k-1))^2)

c=sqrt((ptr(1,k)-pmr(1,k-1))^2+(ptr(2,k)-pmr(2,k-1))^2+(ptr(3,k)-pmr(3,k-1))^2)

b=acos((r(k-1)^2+st^2-c^2)/(2*r(k-1)*st))

dq=acos((r(k-1)^2-st^2+c^2)/(2*r(k-1)*c))

if abs(imag(b))>0

b=0.0000001

end

if abs(imag(dq))>0

dq=0.0000001

end

q(k)=q(k-1)+dq

o(k)=o(k-1)+m*dq

a(k)=o(k)-q(k)

c1=r(k-1)*sin(b)/sin(a(k)+b)

c2=r(k-1)*sin(a(k))/sin(a(k)+b)

c3=sqrt((c1-sm)^2+(c2-st)^2+2*(c1-sm)*(c2-st)*cos(a(k)+b))

dq=a(k)-acos(((c1-sm)^2+c3^2-(c2-st)^2)/(2*(c1-sm)*c3))

if abs(imag(dq))>0

dq=0.0000001

end

q(k)=q(k-1)+dq

o(k)=o(k-1)+m*dq

a(k)=o(k)-q(k)

c1=r(k-1)*sin(b)/sin(a(k)+b)

c2=r(k-1)*sin(a(k))/sin(a(k)+b)

c3=sqrt((c1-sm)^2+(c2-st)^2+2*(c1-sm)*(c2-st)*cos(a(k)+b))

dq=a(k)-acos(((c1-sm)^2+c3^2-(c2-st)^2)/(2*(c1-sm)*c3))

if abs(imag(dq))>0

dq=0.0000001

end

q(k)=q(k-1)+dq

o(k)=o(k-1)+m*dq

a(k)=o(k)-q(k)

c1=r(k-1)*sin(b)/sin(a(k)+b)

c2=r(k-1)*sin(a(k))/sin(a(k)+b)

c3=sqrt((c1-sm)^2+(c2-st)^2+2*(c1-sm)*(c2-st)*cos(a(k)+b))

x1(k)=ptr(1,k-1)+c2/st*(ptr(1,k)-ptr(1,k-1))

y1(k)=ptr(2,k-1)+c2/st*(ptr(2,k)-ptr(2,k-1))

z1(k)=ptr(3,k-1)+c2/st*(ptr(3,k)-ptr(3,k-1))

x(k)=pmr(1,k-1)+sm/c1*(x1(k)-pmr(1,k-1))

y(k)=pmr(2,k-1)+sm/c1*(y1(k)-pmr(2,k-1))

z(k)=pmr(3,k-1)+sm/c1*(z1(k)-pmr(3,k-1))

pmr(:,k)=[x(k)y(k)z(k)]

r(k)=sqrt((ptr(1,k)-pmr(1,k))^2+(ptr(2,k)-pmr(2,k))^2+(ptr(3,k)-pmr(3,k))^2)

if r(k)<0.06

break

end

end

sprintf('遭遇时间:%3.1f',0.1*k)

figure(1)

plot3(pmr(1,1:k),pmr(2,1:k),pmr(3,1:k),'k',ptr(1,:),ptr(2,:),ptr(3,:))

axis([0 25 0 5 0 25])

text(x(80),y(80),z(80),'\leftarrow 比例导引')

grid on

原因

所给链接的代码不完整,缺少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)]

结果


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/12343751.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-24
下一篇 2023-05-24

发表评论

登录后才能评论

评论列表(0条)

保存