Fortran Coder

查看: 121|回复: 1

[混编] 有没有大佬会MATLAB调用fortran生成的dll,下面是修改的代码

[复制链接]

1

帖子

1

主题

0

精华

新人

F 币
11 元
贡献
3 点
发表于 2024-9-28 12:13:18 | 显示全部楼层 |阅读模式
[Fortran] 纯文本查看 复制代码
        SUBROUTINE LINEEHL(X,P,H,N,X0,XE,W,E1,EDA0,R,US)
        !DEC$ ATTRIBUTES DLLEXPORT :: LINEEHL
        !DEC$ ATTRIBUTES REFERENCE :: X,P,H
        INTEGER :: N
        DOUBLE PRECISION :: X0,XE,W,E1,EDA0,R,US,ENDA,A1,A2,A3,HM0,DH,PH,B, PAI=3.14159265,Z=0.68,P0=1.96E8
        !COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH
        !COMMON/COM2/EDA0
        !COMMON/COM4/X0,XE
        !COMMON/COM3/E1,PH,B,R
       !real, parameter :: PAI=3.14159265,Z=0.68,P0=1.96E8
                !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/
                W1=W/(E1*R)
                PH=E1*SQRT(0.5*W1/PAI)
                A1=(LOG(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)
                CALL SUBAK(N)
                CALL EHL(N)
                END

                SUBROUTINE EHL(N)
        !DEC$ ATTRIBUTES DLLEXPORT :: EHL
        !IMPLICIT NONE
        INTEGER ::  N,MK,I,KK
        DOUBLE PRECISION :: DX,X0,XE,H2,P2,B,RR,H3,P3,A1,A2,A3,HM0,DH,E1,PH,ENDA,Z,ERP
        DOUBLE PRECISION :: X(N),P(N),H(N),RO(N),POLD(N),EPS(N),EDA(N),V(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)THEN
                 H2=1.E3
                P2=0.0
        ENDIF
        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)
        !DEC$ ATTRIBUTES DLLEXPORT :: OUTHP
        IMPLICIT NONE
        INTEGER ::  I,N
        DOUBLE PRECISION :: X(N),P(N),H(N)
        DO 10 I=1,N
            X(I)=X(I)
            P(I)=P(I)
            H(I)=H(I)
                !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)
        !DEC$ ATTRIBUTES DLLEXPORT :: HREE
     !   IMPLICIT NONE
        INTEGER :: I,N,KK
        DOUBLE PRECISION :: X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N),AK(0:1100)
        DOUBLE PRECISION :: H00,W1,C3,DX,G0,DW,HMIN,H0,DH,HM0,A1,A2,A3,Z,ENDA,EDA0,PAI1
        !COMMON /COM1/ENDA,A1,A2,A3,Z,HM0,DH/COM2/EDA0/
       ! COMMON /COMAK/AK
                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)
        !DEC$ ATTRIBUTES DLLEXPORT :: ITER
       ! IMPLICIT NONE
        INTEGER ::  N, I, K,KK
        DOUBLE PRECISION :: DX, PAI1
        DOUBLE PRECISION :: X(N),P(N),H(N),RO(N),EPS(N),EDA(N),V(N),AK(0:1100),D1,D2,D3,D8,D9,D10,D11,D12
     !   COMMON /COMAK/AK
                !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)
        !DEC$ ATTRIBUTES DLLEXPORT :: VI
       ! IMPLICIT NONE
        INTEGER ::  N, I, J, IJ
        DOUBLE PRECISION :: DX, PAI1, C
        DOUBLE PRECISION :: P(N), V(N)
        DOUBLE PRECISION AK(0:1100)   ! 显式声明 AK 为 DOUBLE PRECISION
      !  COMMON /COMAK/ AK
                PAI1=0.318309886
                C=LOG(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)
        !DEC$ ATTRIBUTES DLLEXPORT :: SUBAK
      !  IMPLICIT NONE
        INTEGER :: MM, I
        DOUBLE PRECISION:: AK(0:1100)   ! 显式声明 AK 为 DOUBLE PRECISION
    !    COMMON /COMAK/ AK
                DO 10 I=0,MM
10      AK(I)=(I+0.5)*(LOG(ABS(I+0.5))-1.)-(I-0.5)*(LOG(ABS(I-0.5))-1.)
                RETURN
                END

                SUBROUTINE FZ(N,P,POLD)
        !DEC$ ATTRIBUTES DLLEXPORT :: FZ
        INTEGER :: N, I
        DOUBLE PRECISION ::  P(N),POLD(N)
                DO 10 I=1,N
10      POLD(I)=P(I)
                RETURN
                END

                SUBROUTINE ERROP(N,P,POLD,ERP)
        !DEC$ ATTRIBUTES DLLEXPORT :: ERROP
        INTEGER :: N, I
        DOUBLE PRECISION ::  P(N),POLD(N)
        DOUBLE PRECISION ::  SD,SUM,ERP
                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


830

帖子

2

主题

0

精华

大宗师

F 币
3905 元
贡献
2323 点
发表于 6 天前 | 显示全部楼层
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-10-6 15:51

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表