一倍是多少：外弹道计算程序

2020年11月20日发(作者：万毅)
Program dM6
C—————6自由度弹丸外弹道飞行计算程序—————
C= = = 积分变量注释
C Y（1）————————————————弹丸飞行时间 T
C Y（2）————————————————弹丸飞行速度 V
C Y（3）————————————————弹道倾角 Q
C Y（4）————————————————质心的X轴坐标 X
C Y（5）————————————————质心的Y轴坐标 Y
C Y（6）————————————————质心的Z轴坐标 Z
C Y（7）————————————————弹丸自转角 r
C Y（8）————————————————弹丸自转角速度 drdt
C Y（9）————————————————偏角 Pasi1
C Y（10）————————————————偏角 Pasi2
C Y（11）————————————————摆动角 fai1
C Y（12）————————————————摆动角 fai2
C Y（13）————————————————摆动角速度 d（fai1）dt
C Y（14）————————————————摆动角速度 d（fai2）dt
C = = = = = = = = = = = = = =
C Pasi=Pasi1+i*Pasi2
C Fai=fai1+i*fai2
C = = = = =* * * * * * * = = = = = = * * * * * * * = = = = = =

REAL M0，Mzz0,My0,Kz0,Ia,Ic
Real Mzzb,Myb,Mxzb
DIMENSION Y(14),DY(14),YA(14),YB(14),Y1(14),Y2(14)
DIMENSION Y3(14)
COMMONSHUMMM
COMMONWBBI1,I3,I0
COMMONWBCN1,SXm(12),Cx0b(12),Cnb(12),Xpb(12)
COMMONWBDN2,DXm(12),Czb(12),Mzzb(12),Myb(12),Mxzb( 12)
COMMONWBEM0,Ic,Ia,Cd0k,S,D,NH,Xcg
COMMONWBFSg,Sd,Cx,Cy0,Mzz0,My0,Kz0
COMMONWBPFIP0,PSI0
C = = = = = = = = = = = = = = FORMAT-01 BEGAIN = = = = = = = = = = = =
1000 FORMAT(2(),13X,8( ‘*’),2X, ‘TYPE OF PROGECTILE IS:’ ,
&2X,A6，2X, 8( ‘*’))
1007 FORMAT(3(),5X,‘ G= ’,F6.2，‘ CD0K= ’,F6.3，‘ BD1= ’,F7.5，
&‘ H01= ’,F6.5，‘ Ic= ’,F7.4)
1008 FORMAT(14X,49(‘=’))
1009 FORMAT(2(),8X,‘T’,6X,‘V’，7X,‘Q’,7X ,‘X’,8X,‘Y’,9X,‘Z’,
&5X,‘drdt’,5X,‘ Cx ’)
1010 FORMAT(3X,F6.2,2X,F6.1,2X,F5.1,3X,F7.1,3X ,F7.1,1X,F7.1,
&3X,F6.1,1X,F8.4)
1019 FORMAT(2(),8X,‘T’,4X,‘fi1p’，4X,‘fi2p’,4X ,‘psi1’,5X,‘ SD ’,
&5X,‘dlta1’, 3X,‘dlta2’，3X,‘ SGg ’,4X , ‘ Sgd ’)
1020 FORMA T(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)
1500 FORMAT(3X,F7.3,2X,F6.1,2X,F7.3,2X,F7.3,2X,F7.3,2X,
&F7.1, 3X,F7.1)

1022 FORMAT(78X)
1090 FORMAT( ‘END’ )
910 FORMAT(A7)
920 FORMAT(2A7)
C = = = = = = = = = = = = = = FORMAT-01 END = = = = = = = = = = = =

OPEN(12,FILE=‘OUT1’) ! 输出文件
OPEN(14,FILE=‘OUT3’)
OPEN(15,FILE=‘’)
OPEN(13,FILE=‘’) ! 打开数据文件
C = = = = 读取数据文件及各数据意义注释= = = =

READ(13，*)M0，V0，Q0
C 弹丸质量（kg）,初速（ms）,射角（度）
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.1415926180.0
S=3.1415926*D*D4.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*V0DTwest
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)
CALL ATM(Y(5))
CALL AERO(Y)
Sg=(Y(8)(2*IaIc*Y(2)))**2Kz0
Sd=2.0*(Cy0-D*D*My0*M0Ic2)(Cy0-Cx+D*D*M0*Mzz0Ia2)
Sgd=1Sg-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
CALL RK(n,h,0,y,dy,ya,yb)
30 CALL RK(n,H,1,Y,DY,YA,YB)
CALL OUTPUT(Y)
nn=nn+1
do 1 i=1, N
y1(i)=y2(i)
y2(i)=y3(i)
1 y3(i)=y(i)
if(.3)then
Ey=5.*nx
Call laq1(n,y1,y2,y3,4,Ey)
d1=57.3*(y3(11)-y3(9))
d2=57.3*(y3(12)-y3(10))
2 write( 14,1071)y3(1),y3(2),y3(3),y3(4),y3(5),y3(6),d1,d2
nx=nx+1
endif
endif
1071 format(3x,f7.3,f7.2,f8.1,f8.5,f8.5,f9.3,f9.3)
C = = = = = = = = = 计算终止条件设定= = = = = = = =
IF(Y(1).GT...Y(5).LE.0.)THEN
WRITE(*,*)Y(4),Y(5),Y(6)
PAUSE
STOP
ENDIF
C = = = = = = = = = = = = = = = = = == = = = = = = =
IF(Y(3).GT.0.) GO TO 30
DO 201 I=1，n
201 Y1(I)=Y(I)
H=-H2
CALL RK(n,H,1,Y,DY,YA,YB)
DO 202 I=1,n
202 Y2(I)=Y(I)
CALL RK(n,H,1,Y,DY,YA,YB)
CALL LAQ1(n,Y1,Y2,Y,3,0.)
MMM=1
CALL OUTPUT(Y)
MMM=0
DO 203 I=1,n
203 Y(I)=Y1(I)
H=-2*H
CALL RK(n,H,0,Y,DY,YA,YB)
20 CALL RK(n,H,1,Y,DY,YA,YB)
IF(Y(5).GT.0.)THEN
CALL OUTPUT(Y)
GO TO 20
END IF
MMM=1
DO 204 I=1，n
204 Y1(I)=Y(I)
H=-2H
CALL RK(n,H,1,Y,DY,YA,YB)
DO 205 I=1,n
205 Y2(I)=Y(I)
CALL RK(n,H,1,Y,DY,YA,YB)
CALL LAQ1(n,Y1,Y2,Y,5,0.)
CALL OUTPUT(Y)
WRITE(15,1090)
STOP
END

C 拉格朗日插值子程序

SUBROUTINE LAQ1(N,U1,U2,U3,K,E)
DIMENSION U1(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)-U 3(K)))
P3=((E-U1(K))*(E-U2(K))((U3(K)-U1(K))*( U3(K)-U2(K)))
DO 10 I=1，n
10 U3(I)=P1*U1(I)+P2*U2(I)+P3*U3(I)
END

C 输出子程序，可以根据需要修改

SUBROUTINE OUTPUT(Y)
REAL M0,Ia,Ic,Mzz0,My0,Kz0
DIMENSION Y(14)
COMMONSHUMMM
COMMONWBBI1,I3,i0
COMMONWBEM0,Ic,Ia,Cd0k,S,D,NH,Xcg
COMMONWBFSg,Sd,Cx,Cy0,Mzz0,My0,Kz0
C = = = = = = = = = = = = = = FORMAT-02 BEGAIN = = = = = = = = = = = =
1011 FORMAT(3X,F6.2,2X,F6.1,2X,F5.1,3 X,F9.3,3X,F9.3,1X,F9.3,
&3X,F6.1,1X,F8.4)
1020 FORMAT(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)
1022 FORMAT(78X)
C = = = = = = = = = = = = = = FORMAT-02 END = = = = = = = = = = = =
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.03.14159265
Sg=(Y(8)(2*IaIc*Y(2)))**2Kz0
Sd=2.*(Cy0-D*D*My0*M0Ic2)(Cy0-Cx+D*D*M0*Mzz0Ia2)
Sgd=1Sg-Sd*(2-Sd)
C = = = = = = = = = = = BEGIN IF = = = = = = = = = = = = = = = = = = =
IF(.0.)THEN
C if(.1) go to 9
9 IF(.0)RETURN
IF(.5)THEN
WRITE( * ,1022)
WRITE(15,1022)
WRITE(12,1022)
I1=0
END IF
I3=I3+1
IF ()THEN
WRITE(15,101 0)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
END IF
ELSE
IF(.5)THEN
WRITE( * ,1022)
WRITE(15,1022)
WRITE(12,1022)
I1=0
END IF
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
END IF
Return
end
C 龙格—库塔法积分子程序
SUBROUTINE RK(N,H,L,Y,DY,B,C)
DIMENSION A(4),Y(N),DY(N),B(N),C(N)
IF(.0) THEN
CALL ATM(Y(5))
CALL AERO(Y)
END IF
A(1)=H2
A(2)=A(1)
A(3)=H
A(4)=H
IF(L)2，2，4
2 DO 3 I=1，N
3 B(I)=Y(I)
CALL DERY(N,Y,DY)
RETURN
4 DO 5 K=1，3
DO 6 I=1，N
C(I)=B(I)+A(K)*DY(I)
6 Y(I)=Y(I)+A(K+1)*DY(I)3.0
5 CALL DERY(N,C,DY)
DO 7 I=1,N
7 Y(I)=Y(I)+A(1)*DY(I)3.0
CALL ATM(Y(5))
CALL AERO(Y)
GO TO 2
END


SUBROUTINE ATM9(Z)
COMMONWBZRO,AC,PA,XM
IF(.9300.)THEN
TM=288.9-0.006328*Z
PP=(1.0-0.000021905*Z)**5.4
ELSE IF(.12000.)THEN
TM=230.-0.006328*(Z-9300.)+0.0 00001172(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-9 300.)-6328.)&32221.057)+0.1939252
)
ELSE
TM=221.5
PP=0.1937254*EXP(-(Z-12000.)6483.305)
END IF
AC=SQRT(1.404*9.806*29.27*TM)
RO=1.206*PP*288.9TM
PA=1.e5*PP
END



SUBROUTINE AERO(Y)
MO,Ia,Ic,My0,Mzz0,Kz0
REAL Mzzb,Myb,Mxzb
REAL Mz,Kz,Kzz,Kxz,Ky,Mzz,My,Mxz
DIMENSION Y(16)
COMMONWBZRO,AC,PA,XM
COMMONWBGbx,by,Kz,Kzz,Kxz,Ky,bz
COMMONWBCN1,SXm(12),Cx0b(12),Cnb(12),Xpb(12)
COMMONWBDN2,DXm(12),Czb(12), Mzzb(12),Myb(12),Mxzb(12)
COMMONWBEM0,Ic,Ia,Cd0k,S,D,NH,Xcg
COMMONWBFSg,SD,Cx,Cy0,Mzz0,My0,Kz0
Xm=y(2)AC
CALL LAGR(SXm,CX0b,N1,Xm,Cx0)
CALL LAGR(SXm,CNb,N1,Xm,CN)
CALL LAGR(SXm,Xpb,N1,Xm,Xcp)
CALL LAGR(DXm,Czb,N2,Xm,Cz)
CALL LAGR(DXm,Mzzb,N2,Xm,Mzz)
CALL LAGR(DXm,Myb,N2,Xm,My)
CALL LAGR(DXm,Mxzb,N2,Xm,Mxz)
DELTA22=(Y(11)-Y(9))* *2+(Y(12)-Y(10))* *2
EF=0.5*RO*SMO
Cx2=2.*Cn
Cx=Cx0+cx2*DELTA22
DELTA0=SQRT(DELTA22)
Cy=Cx0+cx
Mz=Cn*(XcgD-Xcp)
bx=FF*Cx
by=FF*Cy
bz=EF*Cz*D2.0
Kz=FF*M0*D*MzIa
Kzz=FF*D*D*Mzz*M0Ia2.0
Kxz=FF*D*D*Mxz*M0Ia2.
Ky=FF*D*D*My*M0Ia2.0
Cy0=Cy
Mzz0=Mzz
My0=My
Kz0=Kz
END

SUBROUTINE DERY(N,Y,DY)
DIMENSION Y(N),DY(N)
REAL Kz,Kzz,Kxz,Ky,Ia,Ic
COMMONWBEMO,Ic,Ia,Cd0K,S,D,NH,Xcg
COMMONWBGbx,by,Kz,Kzz,Kxz,Ky,bz

CALL ATM(Y(5))
CALL AERO(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
&+IcIa*Y(8)*DY(3)
END


SUBROUTINE LAGR(X0,Y0,N,X,Y)
DIMENSION X0(N),Y0(N)
IF(.X0(N),Y0(N))goto10
Y=Y0(N)
RETURN
IF(.X0(1))goto20
Y=Y0(1)
RETURN
J=1
DO 30 I=1,N-1
IF(.X0(I+1))GOTO30
J=I
GOTO 40
CONTINUE
CONTINUE
Y=Y0(J)+(Y0(J+1)-Y0(J))(X0(J+1)-X0(J))*+(X-X0(J))
RETURN
END