xppaqwe
xppaqwe
2020-06-11 23:41

用MATLAB调用Fortran程序

  • matlab
  • fortran
  • 仿真

我要崩溃了,我不管怎么修改都无法成功运行!求助各位大神!救救孩子!老师给的程序如下

异步自励磁仿真过程计算程序
异步自励磁仿真计算程序ZLC是按上述步骤编写的,用于计算隐极发电机异步自励磁。如需计算同步自励磁,就应对该程序作适当的修改,主要修改计及饱和特性的部分,因为凸极机受饱和的影响。以下是程序ZLC的几点说明。
功能。计算隐极发电机异步自励磁的动态过程。
输入数据。所有的原始数据均存放在数据文件ZLC.DAT中,按自由格式存放,其顺序为
、、、、、、
、、、、、、R、
t、、UN、NDYJG、NBHTX
BU(1), BI(1)
BU(2), BI(2)

         BU(NBHTX), BI(NBHTX)

(3)输出数据。所有的输出数据均存放在文件ZLC.OUT中,共有六列数据,从左至右分别为:时间t(s),(pu), (pu),(pu),(pu),
I(pu),其中I为电流综合相量的模值。
(4)有关变量及数据的意义。与程序SXDL和DKFC兼容,此外还有如下补充:NDYJG为一正数,表示每隔NDYJG个步长输出一次计算结果。因自励磁过程可能很长,如每次计算都输出结果,则输出数据文件可能回十分庞大。
NBHTX为饱和曲线的点数。
BU(50)、BI(50)为饱和特性曲线上各点的纵坐标和横坐标最多可存放50个点。注意,和均需归算为标腰值。
RE为发电机的外接电阻R
为电容C的倒数
PROGRAM main
!integer::KWD
COMMON/BKG/XD,XAD,RXF,RXD,RF,RD,RXFD,XQ,XAQ,RXH,RXQ,RH,RQ,RR
DIMENSION C(6,6),D1(6),D2(6),E(6),X(6),G(3,3),F(3),BU(50),BI(50)
DATA PAI/3.1415927/OMG/1.0/X/6*0./T/0./TY/0./UD/0./UQ/0./UF/0./
DATA X(2)/0.0012/NDJS/0/
OPEN(5,FILE='zlc.txt',STATUS='OLD')
READ(5,*) XD,XAD,RXF,RXD,RF,RD,RXFD
READ(5,*) XQ,XAQ,RXH,RXQ,RH,RQ,RR,RE,XC
READ(5,*) DLTTY,TMAX,UN,NDYJG,NBHTX
OPEN(7,FILE='zlc1.txt',STATUS='OLD')
READ(7,*) (BU(I),I=1,NBHTX)
OPEN(8,FILE='zlc2.txt',STATUS='OLD')
READ(8,*) (BI(I),I=1,NBHTX)
OPEN(6,FILE='zlcout.txt')
WRITE(6,'(A)')' T Id Iq If Ih AA'
DLTT=DLTTY*314.15927
DLTT2=DLTT/2.0
A11=DLTT2*XC+RE
A12=-DLTT2*RE
A21=DLTT2*RE
A22=DLTT2*XC+RE
AAA=A11*A22-A12*A21
A11N=A22/AAA
A12N=-A12/AAA
A21N=-A21/AAA
A22N=A11/AAA
A11P=A11N+DLTT2*A12N
A12P=-DLTT2*A11N+A12N
A21P=A21N+DLTT2*A22N
A22P=-DLTT2*A21N+A22N
100 T=T+DLTT
TY=TY+DLTTY
C1=-UD+RE*X(1)-DLTT2*UQ+DLTT2*RE*X(2)-DLTT2*XC*X(1)
C2=-UQ+RE*X(2)+DLTT2*UD-DLTT2*RE*X(1)-DLTT2*XC*X(2)
C1P=A11N*C1+A12N*C2
C2P=A21N*C1+A22N*C2
AID=X(3)+X(4)-X(1)
AIQ=X(5)+X(6)-X(2)
AII=SQRT(AID**2+AIQ**2)
DO 20 I=2,NBHTX-1
IF(AII.LE.BI(I+1).and.AII.GE.BI(I)) THEN
UU=BU(I)+(BU(I+1)-BU(I))*(AII-BI(I))/(BI(I+1)-BI(I))
XADP=XAD
XAQP=XAQ
XAD=UU/AII
XAQ=UU/AII
XD=XD-XADP+XAD
RXF=RXF-XADP+XAD
RXFD=RXFD-XADP+XAD
RXD=RXD-XADP+XAD
XQ=XQ-XAQP+XAQ
RXH=RXH-XAQP+XAQ
RXQ=RXQ-XAQP+XAQ
END IF
20 CONTINUE
CALL DKFC(0,C,D1,D2,X,E,DLTT2,UD,UQ,UF,OMG)
CALL DKFC(1,C,D1,D2,X,E,DLTT2,UD,UQ,UF,OMG)
AG11=D1(1)
AG12=D2(1)
AG21=D1(2)
AG22=D2(2)
BG1=E(1)
BG2=E(2)
G11=AG11-A11P
G12=AG12-A12P
G21=AG21-A21P
G22=AG22-A22P
F1=C1P-BG1
F2=C2P-BG2
GGG=G11*G22-G12*G21
G11N=G22/GGG
G12N=-G12/GGG
G21N=-G21/GGG
G22N=G11/GGG
UD=G11N*F1+G12N*F2
UQ=G21N*F1+G22N*F2
X(1)=D1(1)*UD+D2(1)*UQ+E(1)
X(2)=D1(2)*UD+D2(2)*UQ+E(2)
X(3)=D1(3)*UD+D2(3)*UQ+E(3)
X(4)=D1(4)*UD+D2(4)*UQ+E(4)
X(5)=D1(5)*UD+D2(5)*UQ+E(5)
X(6)=D1(6)*UD+D2(6)*UQ+E(6)
AA=SQRT(X(1)**2+X(2)**2)
NDJS=NDJS+1
IF(NDJS.EQ.NDYJG) THEN
WRITE(*,'(F7.2,2F9.4)') TY,XAD,AA
WRITE(6,'(F7.2,7F8.4)') TY,X(1),X(2),X(3),X(5),AA
NDJS=0
END IF
IF(TY.LT.TMAX) GOTO 100
STOP
END
!$INCLUDE:'DFKC.F90'
SUBROUTINE DKFC(KWD,C,D1,D2,X,E,DLTT2,UD,UQ,UF,OMG)
!integer::KWD

COMMON/BKG/XD,XAD,RXF,RXD,RF,RD,RXFD,XQ,XAQ,RXH,RXQ,RH,RQ,RR
DIMENSION C(6,6),D1(6),D2(6),E(6),X(6)
IF(KWD.NE.O) GOTO 200
DO 10 I=1,6
DO 10 J=1,6
10 C(I,J)=0.0
C(1,1)=-(XD+DLTT2*RR)
C(1,2)=DLTT2*OMG*XQ
C(1,3)=XAD
C(1,4)=XAD
C(1,5)=-DLTT2*OMG*XAQ
C(1,6)=-DLTT2*OMG*XAQ
C(2,1)=-DLTT2*OMG*XD
C(2,2)=-(XQ+DLTT2*RR)
C(2,3)=DLTT2*OMG*XAD
C(2,4)=DLTT2*OMG*XAD
C(2,5)=XAQ
C(2,6)=XAQ
C(3,1)=-XAD
C(3,3)=RXF+DLTT2*RF
C(3,4)=RXFD
C(4,1)=-XAD
C(4,3)=RXFD
C(4,4)=RXD+DLTT2*RD
C(5,2)=-XAQ
C(5,5)=RXH+DLTT2*RH
C(5,6)=XAQ
C(6,2)=-XAQ
C(6,5)=XAQ
C(6,6)=RXQ+DLTT2*RQ
DO 20 I=2,6
DO 20 J=1,I-1
DO 20 K=J+1,6
20 C(I,K)=C(I,K)-C(I,J)*C(J,K)/C(J,J)
RETURN
200 DO 25 I=1,6
D1(I)=0.0
25 D2(I)=0.0
D1(1)=DLTT2
D2(2)=DLTT2
E(1)=-(XD-DLTT2*RR)*X(1)-DLTT2*OMG*XQ*X(2)+XAD*(X(3)+X(4))+DLTT2*OMG*XAQ*(X(5)+X(6))+DLTT2*UD
E(2)=DLTT2*OMG*XD*X(1)-(XQ-DLTT2*RR)*X(2)-DLTT2*OMG*XAD*(X(3)+X(4))+XAQ*(X(5)+X(6))+DLTT2*UQ
E(3)=-XAD*X(1)+(RXF-DLTT2*RF)*X(3)+RXFD*X(4)+DLTT2*2.0*UF
E(4)=-XAD*X(1)+RXFD*X(3)+(RXD-DLTT2*RD)*X(4)
E(5)=-XAQ*X(2)+(RXH-DLTT2*RH)*X(5)+XAQ*X(6)
E(6)=-XAQ*X(2)+XAQ*X(5)+(RXQ-DLTT2*RQ)*X(6)
DO 30 I=2,6
DO 30 J=1,I-1
D1(I)=D1(I)-C(I,J)*D1(J)/C(J,J)
D2(I)=D2(I)-C(I,J)*D2(J)/C(J,J)
30 E(I)=E(I)-C(I,J)*E(J)/C(J,J)
DO 40 I=6,1,-1
DO 50 J=I+1,6
D1(I)=D1(I)-C(I,J)*D1(J)
D2(I)=D2(I)-C(I,J)*D2(J)
50 E(I)=E(I)-C(I,J)*E(J)
D1(I)=D1(I)/C(I,I)
D2(I)=D2(I)/C(I,I)
40 E(I)=E(I)/C(I,I)
RETURN
END

  • 点赞
  • 回答
  • 收藏
  • 复制链接分享

2条回答