求黄平弹性流体动压润滑一书中fortran程序转成matlab程序，关于弹流润滑的程序，急

PROGRAM LINEEHL
        COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COM4/X0,XE/COM3/E1,PH,B,R
        DATA PAI,Z,P0/3.14159265,0.68,1.96E8/
        DATA N,X0,XE,W,E1,EDA0,R,Us/130,-4.0,1.5,1.0E5,2.2E11,0.05,0.05,1.5/
        OPEN(8,FILE='OUT.DAT',STATUS='UNKNOWN')
        W1=W/(E1*R)
        PH=E1*SQRT(0.5*W1/PAI)
        A1=(ALOG(EDA0)+9.67)
        A2=PH/P0
        A3=0.59/(PH*1.E-9)
        B=4.*R*PH/E1
        ALFA=Z*A1/P0
        G=ALFA*E1
        U=EDA0*US/(2.*E1*R)
        CC1=SQRT(2.*U)
        AM=2.*PAI*(PH/E1)**2/CC1
        ENDA=3.*(PAI/AM)**2/8.
        HM0=1.6*(R/B)**2*G**0.6*U**0.7*W1**(-0.13) 
        WRITE(*,*)N,X0,XE,W,E1,EDA0,R,US
        CALL SUBAK(N)
        CALL EHL(N)
        STOP
        END
        SUBROUTINE EHL(N)
        DIMENSION X(1100),P(1100),H(1100),RO(1100),POLD(1100),EPS(1100),EDA(1100),V(1100)
        COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM4/X0,XE
        COMMON /COM3/E1,PH,B,RR
        MK=1
        DX=(XE-X0)/(N-1.0)
        DO 10 I=1,N
        X(I)=X0+(I-1)*DX
        IF(ABS(X(I)).GE.1.0)P(I)=0.0
        IF(ABS(X(I)).LT.1.0)P(I)=SQRT(1.-X(I)*X(I))
10      CONTINUE
        CALL HREE(N,DX,X,P,H,RO,EPS,EDA,V)
        CALL FZ(N,P,POLD)
14        KK=19
        CALL ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
        MK=MK+1
        CALL ERROP(N,P,POLD,ERP)
        WRITE(*,*)'ERP=',ERP
        IF(ERP.GT.1.E-5.AND.DH.GT.1.E-7)THEN
        IF(MK.GE.50)THEN
        MK=1
        DH=0.5*DH
        ENDIF
        GOTO 14
        ENDIF
        IF(DH.LE.1.E-7)WRITE(*,*)'Pressures are not convergent!!!'
         H2=1.E3
        P2=0.0
        DO 106 I=1,N
        IF(H(I).LT.H2)H2=H(I)
        IF(P(I).GT.P2)P2=P(I)
106     CONTINUE
        H3=H2*B*B/RR
        P3=P2*PH
110        FORMAT(6(1X,E12.6))
120        CONTINUE
        WRITE(*,*)'P2,H2,P3,H3=',P2,H2,P3,H3
        CALL OUTHP(N,X,P,H)
        RETURN
        END
        SUBROUTINE OUTHP(N,X,P,H)
        DIMENSION X(N),P(N),H(N)
        DO 10 I=1,N
        WRITE(8,20)X(I),P(I),H(I)
10      CONTINUE
20      FORMAT(1X,6(E12.6,1X))
        RETURN
        END
        SUBROUTINE HREE(N,DX,X,P,H,RO,EPS,EDA,V)
        DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
        COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/COMAK/AK(0:1100)
        DATA KK,PAI1,G0/0,0.318309886,1.570796325/
        IF(KK.NE.0)GOTO 3
        H00=0.0
3        W1=0.0
        DO 4 I=1,N
4       W1=W1+P(I)
        C3=(DX*W1)/G0
        DW=1.-C3
        CALL VI(N,DX,P,V)
        HMIN=1.E3
        DO 30 I=1,N
        H0=0.5*X(I)*X(I)+V(I)
        IF(H0.LT.HMIN)HMIN=H0
        H(I)=H0
30        CONTINUE
        IF(KK.NE.0)GOTO 32
        KK=1
        DH=0.005*HM0
        H00=-HMIN+HM0
32        IF(DW.LT.0.0)H00=H00+DH
        IF(DW.GT.0.0)H00=H00-DH
        DO 60 I=1,N
        H(I)=H00+H(I)
        EDA(I)=EXP(A1*(-1.+(1.+A2*P(I))**Z))
        RO(I)=(A3+1.35*P(I))/(A3+P(I))
        EPS(I)=RO(I)*H(I)**3/(ENDA*EDA(I))
60        CONTINUE
        RETURN
        END
        SUBROUTINE ITER(N,KK,DX,X,P,H,RO,EPS,EDA,V)
        DIMENSION X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N)
        COMMON /COMAK/AK(0:1100)
        DATA PAI1/0.318309886/
        DO 100 K=1,KK
        D2=0.5*(EPS(1)+EPS(2))
        D3=0.5*(EPS(2)+EPS(3))
        DO 70 I=2,N-1
        D1=D2
        D2=D3
        IF(I.NE.N-1)D3=0.5*(EPS(I+1)+EPS(I+2))
        D8=RO(I)*AK(0)*PAI1
        D9=RO(I-1)*AK(1)*PAI1
        D10=1.0/(D1+D2+(D9-D8)*DX)
        D11=D1*P(I-1)+D2*P(I+1)
        D12=(RO(I)*H(I)-RO(I-1)*H(I-1)+(D8-D9)*P(I))*DX
        P(I)=(D11-D12)*D10
        IF(P(I).LT.0.0)P(I)=0.0
70      CONTINUE
        CALL HREE(N,DX,X,P,H,RO,EPS,EDA,V)
100     CONTINUE
        RETURN
        END
        SUBROUTINE VI(N,DX,P,V)
        DIMENSION P(N),V(N)
        COMMON /COMAK/AK(0:1100)
        PAI1=0.318309886
        C=ALOG(DX)
        DO 10 I=1,N
        V(I)=0.0
        DO 10 J=1,N
        IJ=IABS(I-J)
10        V(I)=V(I)+(AK(IJ)+C)*DX*P(J)
        DO I=1,N
        V(I)=-PAI1*V(I)
        ENDDO
        RETURN
        END

        SUBROUTINE SUBAK(MM)
        COMMON /COMAK/AK(0:1100)
        DO 10 I=0,MM
10      AK(I)=(I+0.5)*(ALOG(ABS(I+0.5))-1.)-(I-0.5)*(ALOG(ABS(I-0.5))-1.)
        RETURN
        END

        SUBROUTINE FZ(N,P,POLD)
        DIMENSION P(N),POLD(N)
        DO 10 I=1,N
10      POLD(I)=P(I)
        RETURN
        END

        SUBROUTINE ERROP(N,P,POLD,ERP)
        DIMENSION P(N),POLD(N)
        SD=0.0
        SUM=0.0
        DO 10 I=1,N
        SD=SD+ABS(P(I)-POLD(I))
        POLD(I)=P(I)
10      SUM=SUM+P(I)
        ERP=SD/SUM
        RETURN
        END

你好，我是有问必答小助手。为了技术专家团更好地为您解答问题，烦请您补充下（1）问题背景详情，（2）您想解决的具体问题，（3）问题相关代码图片或者报错信息。便于技术专家团更好地理解问题，并给出解决方案。

您可以点击问题下方的【编辑】，进行补充修改问题。

你好，可以给我一份么？717387814@qq.com,谢谢

您好，想一起探讨关于黄平书籍弹流润滑方面的知识。加我q314339464哈