C     ============================================
C     Fortran 参考测试程序 - 生成 Rust 对比数据
C     ============================================

      PROGRAM TEST_REF
      IMPLICIT REAL*8 (A-H,O-Z), LOGICAL*1 (L)

C     所有声明必须在可执行语句之前
      DIMENSION XL(3), YL(3)
      DIMENSION A(3), B(3), C(3), R(3), U(3)
      DIMENSION XX(4), YY(4)
      INTEGER I

C     初始化 yint 测试数据
      XL(1) = 0.0D0
      XL(2) = 1.0D0
      XL(3) = 2.0D0
      YL(1) = 0.0D0
      YL(2) = 1.0D0
      YL(3) = 4.0D0

C     初始化 tridag 测试数据
      A(1) = 0.0D0
      A(2) = 1.0D0
      A(3) = 1.0D0
      B(1) = 2.0D0
      B(2) = 2.0D0
      B(3) = 2.0D0
      C(1) = 1.0D0
      C(2) = 1.0D0
      C(3) = 0.0D0
      R(1) = 5.0D0
      R(2) = 6.0D0
      R(3) = 5.0D0

C     初始化 ylintp 测试数据
      XX(1) = 0.0D0
      XX(2) = 1.0D0
      XX(3) = 2.0D0
      XX(4) = 3.0D0
      YY(1) = 1.0D0
      YY(2) = 3.0D0
      YY(3) = 5.0D0
      YY(4) = 7.0D0

C     测试 expo
      WRITE(*,*) '=== EXPO ==='
      DO 100 I = -5, 5
        XVAL = DBLE(I) * 20.0D0
        Y = EXPO(XVAL)
        WRITE(*,'(F10.2,1X,ES23.16)') XVAL, Y
  100 CONTINUE

C     测试 yint
      WRITE(*,*) '=== YINT ==='
      DO 200 I = -1, 6
        X0 = DBLE(I) * 0.5D0
        Y = YINT(XL, YL, X0)
        WRITE(*,'(F10.2,1X,ES23.16)') X0, Y
  200 CONTINUE

C     测试 lagran
      WRITE(*,*) '=== LAGRAN ==='
      DO 300 I = -1, 6
        X0 = DBLE(I) * 0.5D0
        CALL LAGRAN(0.0D0, 1.0D0, 2.0D0, 0.0D0, 1.0D0, 4.0D0, X0, Y)
        WRITE(*,'(F10.2,1X,ES23.16)') X0, Y
  300 CONTINUE

C     测试 eint
      WRITE(*,*) '=== EINT ==='
      DO 400 I = 1, 10
        TVAL = DBLE(I) * 0.5D0
        CALL EINT(TVAL, E1, E2, E3)
        WRITE(*,'(F10.2,3(1X,ES23.16))') TVAL, E1, E2, E3
  400 CONTINUE

C     测试 tridag
      WRITE(*,*) '=== TRIDAG ==='
      CALL TRIDAG(A, B, C, R, U, 3)
      WRITE(*,'(3(1X,ES23.16))') U(1), U(2), U(3)

C     测试 gntk
      WRITE(*,*) '=== GNTK ==='
      DO 500 I = 0, 4
        FR = 1.0D0 + DBLE(I) * 0.5D0
        G1 = GNTK(1, FR)
        G2 = GNTK(2, FR)
        G3 = GNTK(3, FR)
        G4 = GNTK(4, FR)
        WRITE(*,'(F10.2,4(1X,ES23.16))') FR, G1, G2, G3, G4
  500 CONTINUE

C     测试 raph
      WRITE(*,*) '=== RAPH ==='
      GAM = 1.0D0
      Z1 = 0.5D0
      Z2 = 0.3D0
      A1VAL = 1.0D0
      A2VAL = 2.0D0
      DGAM = RAPH(GAM, Z1, Z2, A1VAL, A2VAL)
      WRITE(*,'(6(1X,ES23.16))') GAM, Z1, Z2, A1VAL, A2VAL, DGAM

C     测试 ffcros
      WRITE(*,*) '=== FFCROS ==='
      FF1 = FFCROS(1, 1, 5000.0D0, 1.0D15)
      FF2 = FFCROS(0, 1, 5000.0D0, 1.0D15)
      FF3 = FFCROS(1, 0, 5000.0D0, 1.0D15)
      WRITE(*,'(3(1X,ES23.16))') FF1, FF2, FF3

C     测试 erfcx
      WRITE(*,*) '=== ERFCX ==='
      DO 600 I = 0, 14
        XVAL = DBLE(I)
        Y = ERFCX(XVAL)
        WRITE(*,'(F10.2,1X,ES23.16)') XVAL, Y
  600 CONTINUE

C     测试 erfcin
      WRITE(*,*) '=== ERFCIN ==='
      DO 700 I = 1, 9
        XVAL = DBLE(I) * 0.1D0
        Y = ERFCIN(XVAL)
        WRITE(*,'(F10.2,1X,ES23.16)') XVAL, Y
  700 CONTINUE

C     测试 sghe12
      WRITE(*,*) '=== SGHE12 ==='
      DO 800 I = 1, 5
        FR = 1.0D15 * DBLE(I)
        Y = SGHE12(FR)
        WRITE(*,'(ES23.16,1X,ES23.16)') FR, Y
  800 CONTINUE

C     测试 ylintp
      WRITE(*,*) '=== YLINTP ==='
      DO 900 I = -1, 8
        X0 = DBLE(I) * 0.5D0
        Y = YLINTP(X0, XX, YY, 4, 4)
        WRITE(*,'(F10.2,1X,ES23.16)') X0, Y
  900 CONTINUE

      END