functions.f

Go to the documentation of this file.

 
 
      MODULE functions
 
      USE stel_kinds
 
      implicit none
 
      PUBLIC :: two_power, two_power_gs
      PRIVATE :: check
 
      CONTAINS
 
      REAL(rprec) function two_power(x, b)
 
      REAL(rprec), INTENT(in)                  :: x
      REAL(rprec), DIMENSION(0:20), INTENT(in) :: b
 
      two_power = b(0)*((1 - x**b(1))**b(2))
 
      END FUNCTION
 
      REAL(rprec) function two_power_gs(x, b)
 
      REAL(rprec), INTENT(in)                  :: x
      REAL(rprec), DIMENSION(0:20), INTENT(in) :: b
 
      INTEGER :: i
 
      two_power_gs = 1.0
      DO i = 3, 18, 3
         two_power_gs = two_power_gs +                                         &
     &                  b(i)*exp(-((x - b(i+1))/b(i+2))**2.0_dp)
      END DO
      two_power_gs = two_power_gs*two_power(x,b)
 
      END FUNCTION
 
      FUNCTION function_test()
 
      LOGICAL :: function_test
      REAL(rprec) :: result
      REAL(rprec), DIMENSION(0:20) :: b(0:20) = 0
 
      b(0:2) = (/ 1.0d+0, 10.0d+0, 2.0d+0 /)
      result = two_power(0.0d+0, b)
      function_test = check(1.0d+0,result,1,"two_power()")
      IF (.not.function_test) RETURN
 
      result = two_power(1.0d+0, b)
      function_test = check(0.0d+0,result,2,"two_power()")
      IF (.not.function_test) RETURN
 
      b(0:2) = (/ 1.0d+0, 1.0d+0, 1.0d+0 /)
      result = two_power(0.5d+0, b)
      function_test = check(0.5d+0,result,3,"two_power()")
      IF (.not.function_test) RETURN
 
      b(0:2) = (/ 1.0d+0, 1.0d+0, 2.0d+0 /)
      result = two_power(0.5d+0, b)
      function_test = check(0.25d+0,result,4,"two_power()")
      IF (.not.function_test) RETURN
 
      b(0:5) = (/ 1.0d+0, 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0, 1.0d+0 /)
      result = two_power_gs(0.4d+0, b)
      function_test = check(two_power(0.4d+0, b),                              &
     &                      result,1,"two_power_gs")
      IF (.not.function_test) RETURN
 
      b(0:5) = (/ 1.0d+0, 1.0d+0, 0.0d+0, 1.0d+0, 0.8d+0, 0.1d+0 /)
      result = two_power_gs(0.8d+0, b)
      function_test = check(2.0d+0,result,1,"two_power_gs")
      IF (.not.function_test) RETURN
 
      END FUNCTION
 
      FUNCTION check(expected, recieved, testNum, name)
 
      LOGICAL :: check
      REAL(rprec), INTENT(in) :: expected, recieved
      INTEGER, INTENT(in) :: testnum
      CHARACTER (LEN=*), INTENT(in) :: name
 
      check = expected .eq. recieved
      IF (.not.check) THEN
         write(*,*) "functions.f: ", name, " test", testnum,                   &
     &              "failed."
         write(*,*) "Expected", expected, "Recieved", recieved
      END IF
 
      END FUNCTION
 
      END MODULE