[特殊函数] 高斯-勒让德求积工具

!--------------------------------------------------------------------------------
!EXAMPLE INTEGERAL SIN(X), X FROM 0.0 TO PI
!--------------------------------------------------------------------------------
PROGRAM TEST_INTEGRAL
    IMPLICIT NONE

    REAL(KIND=8)::A,B
    REAL(KIND=8),EXTERNAL::YOUR_FUNC
    REAL(KIND=8),EXTERNAL::GUALEGINTG12
    INTEGER,PARAMETER::M=5

    ! A1   B1/A2 B2/A3  B3/A4...BM-1/AM  BM
    ! A|-----|-----|------|.......|-----|B

    REAL(KIND=8)::AX(M)
    REAL(KIND=8)::BX(M)
    REAL(KIND=8)::DX
    REAL(KIND=8)::OTHERS=0.D0
    REAL(KIND=8)::INTGVAL
    INTEGER::I

    WRITE(*,*)
    WRITE(*,*)
    WRITE(*,FMT='(4X,A)',ADVANCE='NO'),'A = '
    READ(*,FMT='(F15.8)')A
    WRITE(*,FMT='(4X,A)',ADVANCE='NO'),'B = '
    READ(*,FMT='(F15.8)')B

    DX=(B-A)/REAL(M,KIND=8)
    INTGVAL=0.D0
    DO I=1,M
        AX(I)=A+REAL(I-1,KIND=8)*DX
        BX(I)=AX(I)+DX
        INTGVAL=INTGVAL+GUALEGINTG12(YOUR_FUNC,OTHERS,AX(I),BX(I))
    ENDDO

    WRITE(*,FMT='(4X,A10,F12.9)'),'INTGVAL = ',INTGVAL
    READ(*,*)

END PROGRAM TEST_INTEGRAL
!--------------------------------------------------------------------------------

REAL(KIND=8) FUNCTION YOUR_FUNC(X,OTHERS)
    IMPLICIT NONE
    REAL(KIND=8),INTENT(IN)::X,OTHERS
    !................
    !FOR EXAMPLE
    YOUR_FUNC=SIN(X)+OTHERS !OR YOUR_FUNC=X*X+OTHERS
    RETURN
END FUNCTION YOUR_FUNC
!--------------------------------------------------------------------------------

REAL(KIND=8) FUNCTION GUALEGINTG12(YOUR_FUNC,OTHERS,A,B)
    IMPLICIT NONE
    REAL(KIND=8),INTENT(IN)::A,B
    REAL(KIND=8)::YOUR_FUNC,OTHERS
    EXTERNAL::YOUR_FUNC
    REAL(KIND=8)::RSLT
    REAL(KIND=8)::X,APB2,ASB2
    INTEGER::I
!--------------------------------------------------------------------------------
!12-POINT GAUSS-LEGENDRE WEIGHTS AND ABSCISSAS
!--------------------------------------------------------------------------------
    !GAUSS WEIGHTS
    REAL(KIND=8)::WEIGHTS(12)=(/&
     0.4717533638651182719461596148501741D-01,&
     0.1069393259953184309602547181939961D+00,&
     0.1600783285433462263346525295433591D+00,&
     0.2031674267230659217490644558097984D+00,&
     0.2334925365383548087608498989248781D+00,&
     0.2491470458134027850005624360429511D+00,&
     0.2491470458134027850005624360429511D+00,&
     0.2334925365383548087608498989248781D+00,&
     0.2031674267230659217490644558097984D+00,&
     0.1600783285433462263346525295433591D+00,&
     0.1069393259953184309602547181939961D+00,&
     0.4717533638651182719461596148501741D-01 &
    /)
    !GAUSS NODES IN INTERVAL[-1.0, 1.0]
    REAL(KIND=8)::NODES(12)=(/&
    -0.9815606342467192506905490901492809D+00,&
    -0.9041172563704748566784658661190962D+00,&
    -0.7699026741943046870368938332128180D+00,&
    -0.5873179542866174472967024189405343D+00,&
    -0.3678314989981801937526915366437176D+00,&
    -0.1252334085114689154724413694638531D+00,&
     0.1252334085114689154724413694638531D+00,&
     0.3678314989981801937526915366437176D+00,&
     0.5873179542866174472967024189405343D+00,&
     0.7699026741943046870368938332128180D+00,&
     0.9041172563704748566784658661190962D+00,&
     0.9815606342467192506905490901492809D+00 &
    /)
!--------------------------------------------------------------------------------
    RSLT=0.D0
    ASB2=0.5D0*(B-A)
    APB2=0.5D0*(B+A)
    DO I=1,12   
        X=ASB2*NODES(I)+APB2
        RSLT=RSLT+WEIGHTS(I)*YOUR_FUNC(X,OTHERS)
    ENDDO
    RSLT=RSLT*ASB2
    GUALEGINTG12=RSLT
    RETURN
END FUNCTION GUALEGINTG12   
!--------------------------------------------------------------------------------

PROGRAM test
  INTEGER, PARAMETER :: dp=KIND(1.0D0)
  REAL(KIND=dp) :: x=1.2, y=1.2D0     ! (KIND=8),  literal "8" is compiler dependent
  WRITE(*,*) x, y                     ! x, default precision only
END PROGRAM test

WRITE(*,'(4X,A)',ADVANCE='NO') 'A = ' ! [), 'A = '] ---> [) 'A = '],  no ","
READ(*,*) A                           ! (F15.8), formatted READ, dangerous

FUNCTION GUALEGINTG12(YOUR_FUNC, OTHERS, A, B)
  IMPLICIT NONE
  INTEGER, PARAMETER :: dp=KIND(1.0D0)
  REAL(KIND=dp)             :: GUALEGINTG12
  REAL(KIND=dp), EXTERNAL   :: YOUR_FUNC
  REAL(KIND=dp), INTENT(IN) :: OTHERS, A, B