📄 crse111.f90
字号:
PROGRAM CRSE
USE DFLIB
TYPE (rccoord) curpos
PARAMETER (PI=3.1415926)
COMMON w
INTEGER(2) i,i1,i2,i3,i4,i5,i6,i7,i8,qi,jn,jnmax
REAL(8) n0,nk,dn,eps,E,g1,g2,h9,d5,fi,n,n1,n2,mb1,mb2,cf,cs
REAL(8) dx(0:61),dd(0:60),p(0:60),b(0:60),m(0:60),j(0:60),c(0:60) !fortran 的数组是从1开始的C的数组是从0开始的所以要加1
! to initialize */
OPEN(100,FILE="crse.dat")
OPEN(200,FILE='crse.res')
READ(100,*)w,jnmax,i1,i2,i3,i4,i5,i6,i7,i8
READ(100,*)n0,nk,dn,eps,E,h9,mb1,mb2,cf,cs,p(i5),p(i6),p(i7),p(i8),j(i3),j(i4)
DO i=0,w
READ(100,*) dx(i),dd(i)
END DO
CLOSE(100)
DO i=0,w
p(i)=0.0
j(i)=0.0
b(i)=0.0
m(i)=0.0
c(i)=0.0
END DO
!to calculate the dimension of b,and m
DO i=1,w-1
b(i)=dx(i)/(E*PI/64.0*dd(i)*dd(i)*dd(i)*dd(i))
END DO
DO i=0,w-1
m(i)=(PI/8.0*h9*(dd(i)*dd(i)*dx(i)+dd(i+1)*dd(i+1)*dx(i+1))+p(i))/981.0
END DO
OPEN(UNIT=jn,FILE='USER')
! to search for the c.r.s
n=n0
qi=1
jn=0
d5=sss(jn,n,i1,i2,qi,dx,dd,b,m,c,j,cf,cs,mb1,mb2)
fi=d5
DO WHILE((n+dn.LE.nk).AND.(jn.LT.jnmax))
!OPEN(UNIT=jn+13,FILE='USER')
OPEN(UNIT=jn,FILE='USER')
2000 n1=n
n=n+dn
n2=n
d5=sss(jn,n,i1,i2,qi,dx,dd,b,m,c,j,cf,cs,mb1,mb2)
IF(d5*fi.GT.0.0)THEN
GO TO 2000
END IF
DO WHILE(((n1+n2)/2.0-n1).GT.eps)
n=(n1+n2)/2.0
d5=sss(jn,n,i1,i2,qi,dx,dd,b,m,c,j,cf,cs,mb1,mb2)
IF(d5*fi.LT.0.0)THEN
n2=n
ELSE
n1=n
END IF
END DO
IF(abs(d5).GT.1e04)THEN
fi=-fi
CYCLE
END IF
qi=-1
jn=jn+1
d5=sss(jn,n,i1,i2,qi,dx,dd,b,m,c,j,cf,cs,mb1,mb2)
fi=-fi
qi=1
END DO
IF((n+dn).GT.nk)THEN
WRITE(jn,*)"*** The c.r.s. doesn't exit whthin nk ***",d5
ELSE
WRITE(jn,*)"********** The normal end **********"
END IF
CLOSE(200)
END PROGRAM CRSE
!***** The function sss to calculate the rest c.q. ******/
REAL(8) FUNCTION sss(jn,n,i1,i2,qi,dx,d,b,m,c,j,cf,cs,mb1,mb2)
PARAMETER(PI=3.1415926)
COMMON w
INTEGER(2) qi,i1,i2,jn
REAL(8) n,mb1,mb2,cf,cs
REAL(8) dx(0:60),d(0:60),b(0:60),m(0:60),c(0:60),j(0:60)
INTEGER(2) i
REAL(8) d5,s,mm1,q1,yy1,mp1,th1,mm2,q2,yy2,mp2,th2
REAL(8) a1(0:60),a2(0:60),y(0:60)
! to calculate the rotary speed and stiffness
!W=11
s=PI/30.0*n
s=s*s
c(i1)=cf*(cs-mb1*s/981.0)/(cf+cs-mb1*s/981.0)
c(i2)=cf*(cs-mb2*s/981.0)/(cf+cs-mb2*s/981.0)
! to calculate the first group of coefficients
mm1=0.0
q1=0.0
yy1=0.0
th1=1.0
DO i=1,w
q1=q1+m(i-1)*yy1*s
mp1=mm1+th1*j(i-1)*s
mm1=mp1+q1*dx(i)
yy1=yy1+th1*dx(i)+b(i)*dx(i)*(mm1/6.0+mp1/3.0)
th1=th1+b(i)*(mm1+mp1)/2.0
q1=q1-c(i)*yy1
IF(qi.LT.0) THEN
a1(i)=yy1
END IF
END DO
! to calculate the second group of coefficients
mm2=0.0
q2=0.0
th2=0.0
yy2=1.0
DO i=1,w
q2=q2+m(i-1)*yy2*s
mp2=mm2+th2*j(i-1)*s
mm2=mp2+q2*dx(i)
yy2=yy2+th2*dx(i)+b(i)*dx(i)*(mm2/6.0+mp2/3.0)
th2=th2+b(i)*(mm2+mp2)/2.0
q2=q2-c(i)*yy2
IF(qi.LT.0) THEN
a2(i)=yy2
END IF
END DO
! to calculate the rest c.q. and output the c.r.s
d5=mm1-mm2*q1/q2
WRITE(jn,*)'D5=',d5,'N=',n
IF(qi.LT.0)THEN
WRITE(jn-1,*)"***** Nkp(",jn,")=",n
y(0)=-q1/q2
DO i=1,w
y(i)=a1(i)-a2(i)*q1/q2
END DO
IF(jn.EQ.1)THEN
WRITE(200,*)'********************************************************'
WRITE(200,*)' THE RESULTS OF CALCULATING THE CRITICAL ROTARY SPEED'
WRITE(200,*)'*******************************************************'
END IF
WRITE(200,*)' ***** Nkp(',jn,')=',n,'*****'
DO i=0,w,3
WRITE(200,*)'y( ',i,')=',y(i)
IF(i+1.GT.w)THEN
CYCLE
ELSE
WRITE(200,*)'y(',i+1,')=',y(i+1)
END IF
IF(i+2.GT.w)THEN
CYCLE
ELSE
WRITE(200,*)'y(',i+2,')=',y(i+2)
END IF
END DO
OPEN(UNIT=jn+3,FILE='USER')
CALL drawing(d,dx,y) ! output the results by drawing
IF(jn.EQ.3)THEN
WRITE(200,*)'----------------------- THE END -----------------------'
END IF
END IF
!CALL drawing(d,dx,y) ! output the results by drawing
sss=d5
RETURN
END FUNCTION SSS
!The function to draw
REAL(8) FUNCTION drawing(d,dx,y)
USE DFLIB
!TYPE (rccoord) curpos !to output the text
TYPE (rccoord) rc
PARAMETER(PI=3.1415926)
COMMON w
REAL(8) d(0:60),dx(0:60),y(0:60)
INTEGER(4) i,l1,l2,h1,h2,h3,h4
REAL(8) s1,s2,s3,gdx,dmax,ymax
INTEGER(2) status !to draw line
TYPE (xycoord) xy !to draw line
INTEGER(2) oldcolor
!!!!!!!!!!!!创建一个子窗口
! To calculate the relative coefficients
gdx=0.0
dmax=0.0
ymax=0.0
l1=50
h3=250
DO i=0,w
gdx=gdx+dx(i)
IF(ymax.LT.abs(y(i)))THEN
ymax=abs(y(i))
END IF
IF(dmax.LT.d(i))THEN
dmax=d(i)
END IF
END DO
s1=540.0/gdx
s2=50.0/dmax
s3=abs(70.0/ymax)
! to draw the central line and write letters */
CALL MOVETO(INT2(30), INT2(80), xy)
status = LINETO(INT2(165), INT2(80))
CALL MOVETO(INT2(169), INT2(80), xy)
status = LINETO(INT2(171), INT2(80))
CALL MOVETO(INT2(175), INT2(80), xy)
status = LINETO(INT2(310), INT2(80))
CALL MOVETO(INT2(314), INT2(80), xy)
status = LINETO(INT2(316), INT2(80))
CALL MOVETO(INT2(320), INT2(80), xy)
status = LINETO(INT2(455), INT2(80))
CALL MOVETO(INT2(459), INT2(80), xy)
status = LINETO(INT2(461), INT2(80))
CALL MOVETO(INT2(465), INT2(80), xy)
status = LINETO(INT2(600), INT2(80))
oldcolor = SETTEXTCOLOR(INT2(3))
CALL SETTEXTPOSITION (INT2(11), INT2(6), rc)
CALL OUTTEXT( '1')
CALL SETTEXTPOSITION (INT2(20), INT2(6), rc)
CALL OUTTEXT( '0')
CALL SETTEXTPOSITION (INT2(15), INT2(6), rc)
CALL OUTTEXT( '-1')
CALL SETTEXTPOSITION (INT2(9.5), INT2(6), rc)
CALL OUTTEXT( 'Y/Ymax')
!to draw the shaft and the curve
h3=240-int(y(0)*s3)
DO i=1,w-1
l2=l1+INT2(dx(i)*s1)
h1=80-INT2(d(i)*s2)
h2=80+INT2(d(i)*s2)
IF(d(i-1).NE.d(i))THEN
CALL MOVETO(INT2(l1), INT2(h1), xy)
status=LINETO(INT2(l1), INT2(h2))
END IF
CALL MOVETO(INT2(l1), INT2(h2), xy)
status=LINETO(INT2(l2), INT2(h2))
CALL MOVETO(INT2(l1), INT2(h1), xy)
status=LINETO(INT2(l2), INT2(h1))
h4=240-INT2(y(i)*s3)
CALL MOVETO(INT2(l1), INT2(h3), xy)
status=LINETO(INT2(l2), INT2(h4))
l1=l2
h3=h4
IF(d(i).EQ.d(i+1)) CYCLE
CALL MOVETO(INT2(l2), INT2(h1), xy)
status=LINETO(INT2(l2), INT2(h2))
END DO
CALL MOVETO(INT2(50), INT2(240), xy)
status=LINETO(INT2(l2), INT2(240))
status=RECTANGLE($GBORDER,INT2(50),INT2(170),INT2(l2),INT2(310))
RETURN
END FUNCTION drawing
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -