悬链线FORTRAN代码

      PROGRAM main
      IMPLICIT none
      REAL Z,T,A,E,TEMP,ET,XLO,WO
      INTEGER MEMNO,NPTS,IMPRE
            MEMNO=1
            NPTS=100
            IMPRE=0
            
            Z=60.
            T=60.
            A=1
            E=30000000
            TEMP=100
            ET=0.0000065
            XLO=100
            WO=1

            CALL PCAFX2(MEMNO,Z,T,A,E,TEMP,ET,XLO,WO,IMPRE,NPTS)

            stop
      end program main

      SUBROUTINE PCAFX2(MEMNO,Z,T,A,E,TEMP,ET,XLO,WO,IMPRE,NPTS)
            IMPLICIT NONE
            
            REAL Z,T,A,E,TEMP,ET,XLO,WO
            INTEGER MEMNO,NPTS,IMPRE
            REAL EPS,PI,H,V,XL,XLIN,W,CORD,AMBDA
            REAL F01,F02,F03,F04,COT,DF1,DF2,DET
            REAL TI,TJ,F,FF,G
            REAL AAH,AH,BV,A1,A2,B2,CA,CB,ACA,ACB
            REAL XLAFST,ELONG,ELONGF,SUBXL,MN
            INTEGER ICODE,LL,KK,MM
            REAL XCOOR(3,NPTS),FOC(4)
            INTEGER I,J

            EPS=0.00001
            DO 20 I=1,3
            DO 20 J=1,NPTS
20          XCOOR(I,J)=0.
            PI=3.14159265
            H=Z
            V=T
            ICODE=3
            
            XL=XLO*(1.+ET*TEMP)
            XLIN=XL
            W=WO*XLO/XL

            IF(IMPRE.EQ.1) GO TO 1
            WRITE(6,500)MEMNO
500         FORMAT(1H1,10X,'ANALYSIS OF CABLE NO ',I5)
            WRITE(6,501)
501         FORMAT(10X,'-------------------------')
            WRITE(6,350)
350         FORMAT(/1HO,'CABLE DATA IN THE PLANE OF THE CABLE'/)
            WRITE(6,502)H,V,A,E
502         FORMAT(3X,'H=',F8.2,2X,'V=',F8.2,2X,'A=',F8.6,2X,'E=',E8.2)
            WRITE(6,503)TEMP,ET,XLO,WO
503   FORMAT(3X,'TEMP=',F6.2,2X,'ET=',F8.2,2X,'XLO=',F8.2,2X,'WO=',F8.5)
            WRITE(6,504)EPS,ICODE,NPTS
504         FORMAT(3X,'EPS=',F8.2,2X,'ICODE=',I1,2X,'NPTS=',I5,/)
            WRITE(6,505)XL,W
505   FORMAT(2X,'AFTER TEMPERATURE EFFECT,XL=',F8.2,2X,'W=',F8.5,/)

1           LL=0

            KK=0
            IF(V.LE.0.) GO TO 2
            KK=1
            V=-V
            H=-H

2           CORD=SQRT(H*H+V*V)
            AMBDA=1000000
            IF(H.EQ.0.) GO TO 4
            IF(XL.LE.CORD) GO TO 3
            AMBDA=SQRT(((XL*XL-V*V)/(H*H)-1.)*3.)
            GO TO 4
3           AMBDA=0.20

4           F01=-(W*H)/(2*AMBDA)
            COT=1./TANH(AMBDA)
            F02=(W/2.)*(-V*COT+XL)
            DF1=0.
            DF2=0.

5           F01=F01+DF1
            F02=F02+DF2
6           F04=W*XL-F02
            F03=-F01
            TI=SQRT(F01*F01+F02*F02)
            TJ=SQRT(F03*F03+F04*F04)
            F=F04+TJ
            FF=TI-F02
            IF(FF.LT.0.0001) FF=0.0001
            G=F/FF
            IF(G.LT.0.0001) G=0.0001

            AAH=(1./W)*ALOG(G)+XL/(A*E)
            AH=-F01*AAH
            BV=(TJ*TJ+TI*TI)/(2.*E*A*W)+(TJ-TI)/W
            IF(LL.EQ.1) GO TO 16

            CA=H-AH
            CB=V-BV
            ACA=ABS(CA)
            ACB=ABS(CB)
            IF(ACA.LE.EPS.AND.ACB.LE.EPS) GO TO 10

            IF(TJ.LT.0.0001) TJ=0.0001
            B2=-(F02/TI+F04/TJ)/W-XL/(E*A)
            A1=-AAH-B2-XL/(E*A)
            A2=(F01/W)*(1./TJ-1./TI)
            DET=A1*B2-A2*A2
            DF1=(CA*B2-CB*A2)/DET
            DF2=(A1*CB-A2*CA)/DET
            GO TO 5

10          XLAFST=XL+(F04*TJ+F02*TI*F01*F01*ALOG(G))/(1.*E*A*W)
            FOC(1)=F01*(1.-2.*FLOAT(KK))
            FOC(2)=F02+FLOAT(KK)*(F04-F02)
            FOC(3)=F03*(1.-2.*FLOAT(KK))
            FOC(4)=F04+FLOAT(KK)*(F02-F04)
            IF(IMPRE.EQ.1) GO TO 18
            WRITE(6,508) XLAFST
508   FORMAT(/2X,'LENGTH AFTER TEMPERATURE EFFECT AND STRETCHING=',F8.2)
            ELONG=((XLAFST-XLIN)/XLIN)*100.
            WRITE(6,509) ELONG
509         FORMAT(2X,'ELONGATION =',F8.4,'PER CENT')
            ELONGF=XLAFST/XLIN
            WRITE(6,510) ELONGF
510         FORMAT(2X,'LENGTH FACTOR =',F8.4,/)
            LL=1

            SUBXL=XLIN/(FLOAT(NPTS)-1.)
            XL=-SUBXL
            F01=FOC(1)-(FOC(1)-FOC(3))*FLOAT(KK)
            F02=FOC(2)-(FOC(2)-FOC(4))*FLOAT(KK)
            DO 17 MM=1,NPTS
            GO TO 6
16          MN=MM+(NPTS-2*MM+1)*KK
            XCOOR(1,MN)=AH+Z*FLOAT(KK)
            XCOOR(2,MN)=BV+T*FLOAT(KK)
17          CONTINUE
18          RETURN
      END