外弹道计算程序.docx
- 文档编号:6738924
- 上传时间:2023-01-09
- 格式:DOCX
- 页数:14
- 大小:18.54KB
外弹道计算程序.docx
《外弹道计算程序.docx》由会员分享,可在线阅读,更多相关《外弹道计算程序.docx(14页珍藏版)》请在冰豆网上搜索。
外弹道计算程序
ProgramdM6
C-----6自由度弹丸外弹道飞行计算程序-----
C===积分变量注释
CY
(1)----------------弹丸飞行时间T
CY
(2)----------------弹丸飞行速度V
CY(3)----------------弹道倾角Q
CY(4)----------------质心的X轴坐标X
CY(5)----------------质心的Y轴坐标Y
CY(6)----------------质心的Z轴坐标Z
CY(7)----------------弹丸自转角r
CY(8)----------------弹丸自转角速度dr/dt
CY(9)----------------偏角Pasi1
CY(10)----------------偏角Pasi2
CY(11)----------------摆动角fai1
CY(12)----------------摆动角fai2
CY(13)----------------摆动角速度d(fai1)/dt
CY(14)----------------摆动角速度d(fai2)/dt
C==============
CPasi=Pasi1+i*Pasi2
CFai=fai1+i*fai2
C=====*******======*******======
REALM0,Mzz0,My0,Kz0,Ia,Ic
RealMzzb,Myb,Mxzb
DIMENSIONY(14),DY(14),YA(14),YB(14),Y1(14),Y2(14)
DIMENSIONY3(14)
COMMON/SHU/MMM
COMMON/WBB/I1,I3,I0
COMMON/WBC/N1,SXm(12),Cx0b(12),Cnb(12),Xpb(12)
COMMON/WBD/N2,DXm(12),Czb(12),Mzzb(12),Myb(12),Mxzb(12)
COMMON/WBE/M0,Ic,Ia,Cd0k,S,D,NH,Xcg
COMMON/WBF/Sg,Sd,Cx,Cy0,Mzz0,My0,Kz0
COMMON/WBP/FIP0,PSI0
C==============FORMAT-01BEGAIN============
1000FORMAT(2(/),13X,8('*'),2X,'TYPEOFPROGECTILEIS:
',
&2X,A6,2X,8('*'))
1007FORMAT(3(/),5X,'G=',F6.2,'CD0K=',F6.3,'BD1=',F7.5,
&'H01=',F6.5,'Ic=',F7.4)
1008FORMAT(14X,49('='))
1009FORMAT(2(/),8X,'T',6X,'V',7X,'Q',7X,'X',8X,'Y',9X,'Z',
&5X,'dr/dt',5X,'Cx'//)
1010FORMAT(3X,F6.2,2X,F6.1,2X,F5.1,3X,F7.1,3X,F7.1,1X,F7.1,
&3X,F6.1,1X,F8.4)
1019FORMAT(2(/),8X,'T',4X,'fi1p',4X,'fi2p',4X,'psi1',5X,'SD',
&5X,'dlta1',3X,'dlta2',3X,'SGg',4X,'Sgd'/)
1020FORMAT(3X,F7.3,2X,F7.4,1X,F7.4,2X,F7.4,2X,F7.4,1X,F7.3,&2X,F7.3,2X,F6.3,1X,F6.3)
1500FORMAT(3X,F7.3,2X,F6.1,2X,F7.3,2X,F7.3,2X,F7.3,2X,
&F7.1,3X,F7.1)
1022FORMAT(78X)
1090FORMAT('END')
910FORMAT(A7)
920FORMAT(2A7)
C==============FORMAT-01END============
OPEN(12,FILE='OUT1')!
输出文件
OPEN(14,FILE='OUT3')
OPEN(15,FILE='OUT.TXT')
OPEN(13,FILE='d6.dat')!
打开数据文件
C====读取数据文件及各数据意义注释====
READ(13,*)M0,V0,Q0
C弹丸质量(kg),初速(m/s),射角(度)
READ(13,*)D,Ia,Ic,Twst,Xcg
C弹丸直径(m),赤道和极转动惯量(kg.m2),缠度,质心位置(至弹顶m)
READ(13,*)H,NH,IPROP,CD0K
C步长(s),(其余填0,1,1)
READ(13,*)BD1,H01,BA1,BR1
C全部填0,0,0,0
READ(13,*)FIP01,fip02,psi01,psi02
C起始扰动角速度(弧度/s),起始偏角(弧度)
READ(13,*)N1,(SXM(I),I=1,N1)
C静态系数M数的个数与M(马赫数)
READ(13,*)(Cx0b(I),I=1,N1)
C阻力系数
READ(13,*)(Cnb(I),I=1,N1)
C法向力系数
READ(13,*)(Xpb(I),I=1,N1)
C压心位置
READ(13,*)N2,(DXM(I),I=1,N2)
C动态系数M数的个数与M(马赫数)
READ(13,*)(Czb(I),I=1,N2)
C马氏力系数
READ(13,*)(Mzzb(I),I=1,N2)
C赤道阻压力矩系数
READ(13,*)(Myb(I),I=1,N2)
C马氏力矩系数的导数
READ(13,*)(Mxzb(I),I=1,N2)
C极阻压力矩系数的导数
C====================================
Close(13)
N=14
Q1=Q0*3.1415926/180.0
S=3.1415926*D*D/4.0
WRITE(*,1022)
T0=0.
Y
(1)=T0
Y
(2)=V0
Y(3)=Q1
Y(4)=0.0
Y(5)=0.0
Y(6)=0.0
Y(7)=0.0
Y(8)=2.0*3.1415926*V0/D/Twest
Y(9)=psi01
Y(10)=PSI02
Y(11)=0.0
Y(12)=0.0
Y(13)=FIP01
Y(14)=fip02
I0=0
I1=0
I2=0
I3=0
WRITE(15,1007)m0,cd0k,bd1,h01,ic
WRITE(15,1008)
WRITE(*,1009)
DLTA1=Y(11)-Y(9)
DLTA2=Y(12)-Y(10)
DELTA=SQRT(DLTA1**2+DLTA2**2)
CALLATM(Y(5))
CALLAERO(Y)
Sg=(Y(8)/(2*Ia/Ic*Y
(2)))**2/Kz0
Sd=2.0*(Cy0-D*D*My0*M0/Ic/2)/(Cy0-Cx+D*D*M0*Mzz0/Ia/2)
Sgd=1/Sg-Sd*(2-Sd)
WRITE(*,1010)T0,V0,Y(4),Y(5),Y(6),Y(8),Cx
WRITE(15,1010)T0,V0,Y(4),Y(5),Y(6),Y(8),Cx
WRITE(12,1008)
WRITE(12.1019)
Y9=Y(9)*57.3
Y10=Y(10)*57.3
DLTA1=DLTA1*57.3
DLTA2=DLTA2*57.3
WRITE(12,1020)T0,y(13),y(14),Y9,Sd,DLTA1,DLTA2,Sg,Sgd
I1=i1=1
MMM=0
nn=0
nx=1
CALLRK(n,h,0,y,dy,ya,yb)
30CALLRK(n,H,1,Y,DY,YA,YB)
CALLOUTPUT(Y)
nn=nn+1
do1i=1,N
y1(i)=y2(i)
y2(i)=y3(i)
1y3(i)=y(i)
if(nn.ge.3)then
Ey=5.*nx
Calllaq1(n,y1,y2,y3,4,Ey)
d1=57.3*(y3(11)-y3(9))
d2=57.3*(y3(12)-y3(10))
2write(14,1071)y3
(1),y3
(2),y3(3),y3(4),y3(5),y3(6),d1,d2
nx=nx+1
endif
endif
1071format(3x,f7.3,f7.2,f8.1,f8.5,f8.5,f9.3,f9.3)
C=========计算终止条件设定========
IF(Y
(1).GT..002.AND.Y(5).LE.0.)THEN
WRITE(*,*)Y(4),Y(5),Y(6)
PAUSE
STOP
ENDIF
C==========================
IF(Y(3).GT.0.)GOTO30
DO201I=1,n
201Y1(I)=Y(I)
H=-H/2
CALLRK(n,H,1,Y,DY,YA,YB)
DO202I=1,n
202Y2(I)=Y(I)
CALLRK(n,H,1,Y,DY,YA,YB)
CALLLAQ1(n,Y1,Y2,Y,3,0.)
MMM=1
CALLOUTPUT(Y)
MMM=0
DO203I=1,n
203Y(I)=Y1(I)
H=-2*H
CALLRK(n,H,0,Y,DY,YA,YB)
20CALLRK(n,H,1,Y,DY,YA,YB)
IF(Y(5).GT.0.)THEN
CALLOUTPUT(Y)
GOTO20
ENDIF
MMM=1
DO204I=1,n
204Y1(I)=Y(I)
H=-2/H
CALLRK(n,H,1,Y,DY,YA,YB)
DO205I=1,n
205Y2(I)=Y(I)
CALLRK(n,H,1,Y,DY,YA,YB)
CALLLAQ1(n,Y1,Y2,Y,5,0.)
CALLOUTPUT(Y)
WRITE(15,1090)
STOP
END
C拉格朗日插值子程序
SUBROUTINELAQ1(N,U1,U2,U3,K,E)
DIMENSIONU1(14),U2(14),U3(14)
P1=((E-U2(K))*(E-U3(K))/((U1(K)-U2(K))*(U1(K)-U3(K)))
P2=((E-U1(K))*(E-U3(K))/((U2(K)-U1(K))*(U2(K)-U3(K)))
P3=((E-U1(K))*(E-U2(K))/((U3(K)-U1(K))*(U3(K)-U2(K)))
DO10I=1,n
10U3(I)=P1*U1(I)+P2*U2(I)+P3*U3(I)
END
C输出子程序,可以根据需要修改
SUBROUTINEOUTPUT(Y)
REALM0,Ia,Ic,Mzz0,My0,Kz0
DIMENSIONY(14)
COMMON/SHU/MMM
COMMON/WBB/I1,I3,i0
COMMON/WBE/M0,Ic,Ia,Cd0k,S,D,NH,Xcg
COMMON/WBF/Sg,Sd,Cx,Cy0,Mzz0,My0,Kz0
C==============FORMAT-02BEGAIN============
1011FORMAT(3X,F6.2,2X,F6.1,2X,F5.1,3X,F9.3,3X,F9.3,1X,F9.3,&3X,F6.1,1X,F8.4)
1020FORMAT(3X,F7.3,2X,F7.4,1X,F7.4,2X,F7.4,2X,F7.4,1X,F7.3
&2X,F7.1,2X,F6.3,1X,F6.3)
1022FORMAT(78X)
C==============FORMAT-02END============
DLTA1=(Y(11)-Y(9))*57.3
DLTA2=(Y(12)-Y(10))*57.3
DLTA=SQRT(DLTA1**2+DELT2**2)
Y(9)=Y(9)*57.3
Y(10)=Y(10)*57.3
QQ=Y(3)*180.0/3.14159265
Sg=(Y(8)/(2*Ia/Ic*Y
(2)))**2/Kz0
Sd=2.*(Cy0-D*D*My0*M0/Ic/2)/(Cy0-Cx+D*D*M0*Mzz0/Ia/2)
Sgd=1/Sg-Sd*(2-Sd)
C===========BEGINIF===================
IF(MMM.EQ.0.)THEN
Cif(i0.eq.1)goto9
9IF(NH.EQ.0)RETURN
IF(I1.EQ.5)THEN
WRITE(*,1022)
WRITE(15,1022)
WRITE(12,1022)
I1=0
ENDIF
I3=I3+1
IF(I3.EQ.NH)THEN
WRITE(15,1010)Y
(1),Y
(2),QQ,Y(4),Y(5),Y(6),Y(8),Cx
WRITE(12,1020)y
(1),y(13),y(14),Y9,Sd,DLTA1,DLTA2,Sg,Sgd
I1=I1+1
I3=0
ENDIF
ELSE
IF(I1.EQ.5)THEN
WRITE(*,1022)
WRITE(15,1022)
WRITE(12,1022)
I1=0
ENDIF
WRITE(*,1010)Y
(1),Y
(2),QQ,Y(4),Y(5),Y(6),Y(8),Cx
WRITE(15,1010)Y
(1),Y
(2),QQ,Y(4),Y(5),Y(6),Y(8),Cx
WRITE(12,1020)y
(1),y(13),y(14),Y9,Sd,DLTA1,DLTA2,Sg,Sgd
I1=I1+1
I3=0
ENDIF
Return
end
C龙格-库塔法积分子程序
SUBROUTINERK(N,H,L,Y,DY,B,C)
DIMENSIONA(4),Y(N),DY(N),B(N),C(N)
IF(L.EQ.0)THEN
CALLATM(Y(5))
CALLAERO(Y)
ENDIF
A
(1)=H/2
A
(2)=A
(1)
A(3)=H
A(4)=H
IF(L)2,2,4
2DO3I=1,N
3B(I)=Y(I)
CALLDERY(N,Y,DY)
RETURN
4DO5K=1,3
DO6I=1,N
C(I)=B(I)+A(K)*DY(I)
6Y(I)=Y(I)+A(K+1)*DY(I)/3.0
5CALLDERY(N,C,DY)
DO7I=1,N
7Y(I)=Y(I)+A
(1)*DY(I)/3.0
CALLATM(Y(5))
CALLAERO(Y)
GOTO2
END
SUBROUTINEATM9(Z)
COMMON/WBZ/RO,AC,PA,XM
IF(Z.LE.9300.)THEN
TM=288.9-0.006328*Z
PP=(1.0-0.000021905*Z)**5.4
ELSEIF(Z.LE.12000.)THEN
TM=230.-0.006328*(Z-9300.)+0.000001172(Z-9300.)*(Z-9300.)
PP=0.2922575*EXP(-2.1206426*(ATAN((2.344*(Z-9300.)-6328.)/
&32221.057)+0.1939252))
ELSE
TM=221.5
PP=0.1937254*EXP(-(Z-12000.)THEN
TM=230.-0.006328*(Z-9300.)+0.000001172*(Z-9300.)*(Z-9300)
PP=0.2922575*EXP(-2.1206426)*(ATAN((2.344*(Z-9300.)-6328.)/&32221.057)+0.1939252)
ELSE
TM=221.5
PP=0.1937254*EXP(-(Z-12000.)/6483.305)
ENDIF
AC=SQRT(1.404*9.806*29.27*TM)
RO=1.206*PP*288.9/TM
PA=1.e5*PP
END
SUBROUTINEAERO(Y)
MO,Ia,Ic,My0,Mzz0,Kz0
REALMzzb,Myb,Mxzb
REALMz,Kz,Kzz,Kxz,Ky,Mzz,My,Mxz
DIMENSIONY(16)
COMMON/WBZ/RO,AC,PA,XM
COMMON/WBG/bx,by,Kz,Kzz,Kxz,Ky,bz
COMMON/WBC/N1,SXm(12),Cx0b(12),Cnb(12),Xpb(12)
COMMON/WBD/N2,DXm(12),Czb(12),Mzzb(12),Myb(12),Mxzb(12)
COMMON/WBE/M0,Ic,Ia,Cd0k,S,D,NH,Xcg
COMMON/WBF/Sg,SD,Cx,Cy0,Mzz0,My0,Kz0
Xm=y
(2)/AC
CALLLAGR(SXm,CX0b,N1,Xm,Cx0)
CALLLAGR(SXm,CNb,N1,Xm,CN)
CALLLAGR(SXm,Xpb,N1,Xm,Xcp)
CALLLAGR(DXm,Czb,N2,Xm,Cz)
CALLLAGR(DXm,Mzzb,N2,Xm,Mzz)
CALLLAGR(DXm,Myb,N2,Xm,My)
CALLLAGR(DXm,Mxzb,N2,Xm,Mxz)
DELTA22=(Y(11)-Y(9))**2+(Y(12)-Y(10))**2
EF=0.5*RO*S/MO
Cx2=2.*Cn
Cx=Cx0+cx2*DELTA22
DELTA0=SQRT(DELTA22)
Cy=Cx0+cx
Mz=Cn*(Xcg/D-Xcp)
bx=FF*Cx
by=FF*Cy
bz=EF*Cz*D/2.0
Kz=FF*M0*D*Mz/Ia
Kzz=FF*D*D*Mzz*M0/Ia/2.0
Kxz=FF*D*D*Mxz*M0/Ia/2.
Ky=FF*D*D*My*M0/Ia/2.0
Cy0=Cy
Mzz0=Mzz
My0=My
Kz0=Kz
END
SUBROUTINEDERY(N,Y,DY)
DIMENSIONY(N),DY(N)
REALKz,Kzz,Kxz,Ky,Ia,Ic
COMMON/WBE/MO,Ic,Ia,Cd0K,S,D,NH,Xcg
COMMON/WBG/bx,by,Kz,Kzz,Kxz,Ky,bz
CALLATM(Y(5))
CALLAERO(Y)
DY
(1)=1
DY
(2)=-bx*Y
(2)*Y
(2)-9.806*SIN(Y(3)+y(9))
DY(3)=-9.806*COS(Y(3))/Y
(2)
DY(4)=Y
(2)*COS(Y(10))*COS(Y(3)+Y(9))
DY(5)=Y
(2)*COS(Y(10))*SIN(Y(3)+Y(9))
DY(6)=Y
(2)*SIN(Y(10))
DY(7)=Y(8)
DY(8)=-Kxz*Y
(2)*Y(8)
DLT1=Y(11)-Y(9)
DLT2=Y(12)-Y(10)
DY(9)=by*Y
(2)*DLT1+bz*Y(8)*DLT2+9.806*SIN(Y(3))*Y(9)/Y
(2)
DY(10)=by*Y
(2)*DLT2-bz*Y(8)*DLT1+9.806*SIN(Y(3))*Y(10)/Y
(2)
DY(11)=Y(13)
DY(12)=Y(14)
DY(13)=-Ic*Y(8)*Y(14)/Ia+Kz*Y
(2)*Y
(2)*DLT1-Kzz*Y
(2)*Y(13)
&+Ky*Y
(2)*Y(8)*DLT2
&-Kzz*Y
(2)*DY(3)-9.086*(Y
(2)*DY(3)*SIN(Y(3))
&+DY
(2)*COS(Y(3)))/(Y
(2)*Y
(2))
DY(14)=Ic*Y(8)*Y(13)/Ia+Kz*Y
(2)*Y
(2)*DLT2-Kzz*Y
(2)*Y(14)
&-Ky*Y
(2)*Y(8)*DLT1
&+Ic/Ia*Y(8)*DY(3)
END
SUBROUTINELAGR(X0,Y0,N,X,Y)
DIMENSIONX0(N),Y0(N)
IF(X.LT.X0(N),Y0(N))goto10
Y=Y0(N)
RETURN
IF(X.GT.X0
(1))goto20
Y=Y0
(1)
RETURN
J=1
DO30I=1,N-1
IF(X.GT.X0(I+1))GOTO30
J=I
GOTO40
CONTINUE
CONTINUE
Y=Y0(J)+(Y0(J+1)-Y0(J))/(X0(J+1)-X0(J))*+(X-X0(J))
RETURN
END
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 弹道 计算 程序