Skip to content

GitLab

Commit a9709d69 authored by Graham, Aaron's avatar Graham, Aaron
Browse files

Add FunctionTableGenerator

parent 2ec09ea4
Hide whitespace changes
Inline Side-by-side
src/CMakeLists.txt
View file @ a9709d69
......@@ -103,6 +103,7 @@ TRIBITS_ADD_LIBRARY(Utils
SchemaParser.f90
SpeciesElements.f90
FutilityComputingEnvironment.f90
FunctionTableGenerator.f90
NonLinearSolver.f90
SolutionAcceleration.f90
AndersonAcceleration.f90
......
src/ExtendedMath.f90
View file @ a9709d69
......@@ -22,6 +22,7 @@ PUBLIC :: GreatestCommonDivisor
PUBLIC :: LeastCommonMultiple
PUBLIC :: ATAN2PI
PUBLIC :: RotateQtrClockwise
PUBLIC :: ThirdOrderBickley
INTERFACE LeastCommonMultiple
MODULE PROCEDURE LeastCommonMultiple_scalar
......@@ -245,7 +246,7 @@ ELEMENTAL FUNCTION ATAN2PI(x,y) RESULT(theta)
ENDFUNCTION ATAN2PI
!
!-------------------------------------------------------------------------------
!> @brief This routine will rotate the input x and y coordinates @nrot number of
!> @brief This routine will rotate the input x and y coordinates @nrot number of
!> clockwise quarter rotations.
!> @param x the x-coordinate
!> @param y the y-coordinate
......@@ -253,7 +254,7 @@ ENDFUNCTION ATAN2PI
!>
!> @note: Values from [-3,-1] will be handled as well, even though they
!> represent counter clockwise rotations. If nrot == 0, no rotations are
!> performed. The results of these three rotations were taken from the
!> performed. The results of these three rotations were taken from the
!> rotation matrix:
!> [cos(theta) -sin(theta) * [x
!> sin(theta) cos(theta)] y] where theta is degrees, CW.
......@@ -279,4 +280,83 @@ ELEMENTAL SUBROUTINE RotateQtrClockwise(x,y,nrot)
ENDSELECT
ENDSUBROUTINE RotateQtrClockwise
!
!-------------------------------------------------------------------------------
!> @brief It calculates third order Bickley function
!> @param x The argument of the third order Bickley function
!> @returns y The result of the third order Bickley function
!>
PURE FUNCTION ThirdOrderBickley(x) RESULT(y)
REAL(SRK),INTENT(IN) :: x
REAL(SRK) :: y
!local variables
REAL(SRK) :: xx
xx=ABS(x)
IF(xx < 0.05_SRK) THEN !Most common case
y=(0.7266088_SRK*xx-0.9990226_SRK)*xx+0.7853961_SRK
ELSEIF(xx > 9.0_SRK) THEN !Second most common case
y=0.0_SRK
ELSEIF(xx < 0.10_SRK) THEN
y=(0.6466375_SRK*xx-0.9912340_SRK)*xx+0.7852024_SRK
ELSEIF(xx < 0.15_SRK) THEN
y=(0.5856605_SRK*xx-0.9791293_SRK)*xx+0.7845986_SRK
ELSEIF(xx < 0.20_SRK) THEN
y=(0.5346648_SRK*xx-0.9638914_SRK)*xx+0.7834577_SRK
ELSEIF(xx < 0.25_SRK) THEN
y=(0.4907827_SRK*xx-0.9463843_SRK)*xx+0.7817094_SRK
ELSEIF(xx < 0.30_SRK) THEN
y=(0.4521752_SRK*xx-0.9271152_SRK)*xx+0.7793031_SRK
ELSEIF(xx < 0.35_SRK) THEN
y=(0.4177388_SRK*xx-0.9064822_SRK)*xx+0.7762107_SRK
ELSEIF(xx < 0.40_SRK) THEN
y=(0.3869945_SRK*xx-0.8849865_SRK)*xx+0.7724519_SRK
ELSEIF(xx < 0.45_SRK) THEN
y=(0.3590753_SRK*xx-0.8626685_SRK)*xx+0.7679903_SRK
ELSEIF(xx < 0.50_SRK) THEN
y=(0.3338676_SRK*xx-0.8400133_SRK)*xx+0.7628988_SRK
ELSEIF(xx < 0.60_SRK) THEN
y=(0.2998569_SRK*xx-0.8054172_SRK)*xx+0.7540982_SRK
ELSEIF(xx < 0.70_SRK) THEN
y=(0.2609154_SRK*xx-0.7587821_SRK)*xx+0.7401279_SRK
ELSEIF(xx < 0.80_SRK) THEN
y=(0.2278226_SRK*xx-0.7125290_SRK)*xx+0.7239594_SRK
ELSEIF(xx < 0.90_SRK) THEN
y=(0.1994999_SRK*xx-0.6672761_SRK)*xx+0.7058777_SRK
ELSEIF(xx < 1.00_SRK) THEN
y=(0.1751248_SRK*xx-0.6234536_SRK)*xx+0.6861762_SRK
ELSEIF(xx < 1.40_SRK) THEN
y=((-0.05337485_SRK*xx+0.3203223_SRK)*xx-0.7538355_SRK)*xx+0.7247294_SRK
ELSEIF(xx < 1.80_SRK) THEN
y=((-0.03146833_SRK*xx+0.2295280_SRK)*xx-0.6279752_SRK)*xx+0.6663720_SRK
ELSEIF(xx < 2.20_SRK) THEN
y=((-0.01906198_SRK*xx+0.1631667_SRK)*xx-0.5094124_SRK)*xx+0.5956163_SRK
ELSEIF(xx < 2.60_SRK) THEN
y=((-0.01174752_SRK*xx+0.1152418_SRK)*xx-0.4046007_SRK)*xx+0.5191031_SRK
ELSEIF(xx < 3.0_SRK) THEN
y=((-0.007328415_SRK*xx+0.08097913_SRK)*xx-0.3159648_SRK)*xx+0.4425954_SRK
ELSEIF(xx < 3.40_SRK) THEN
y=((-0.004617254_SRK*xx+0.05669960_SRK)*xx-0.2434341_SRK)*xx+0.3703178_SRK
ELSEIF(xx < 3.80_SRK) THEN
y=(0.007923547_SRK*xx-0.07158569_SRK)*xx+0.1684022_SRK
ELSEIF(xx < 4.20_SRK) THEN
y=(0.005095111_SRK*xx-0.05016344_SRK)*xx+0.1278307_SRK
ELSEIF(xx < 4.60_SRK) THEN
y=(0.003286040_SRK*xx-0.03501524_SRK)*xx+0.09611422_SRK
ELSEIF(xx < 5.00_SRK) THEN
y=(0.002126242_SRK*xx-0.02437465_SRK)*xx+0.07170491_SRK
ELSEIF(xx < 5.80_SRK) THEN
y=(0.001123687_SRK*xx-0.01425519_SRK)*xx+0.04616317_SRK
ELSEIF(xx < 6.60_SRK) THEN
y=(4.762937e-4_SRK*xx-6.810124e-3_SRK)*xx+0.02475115_SRK
ELSEIF(xx < 7.40_SRK) THEN
y=(2.031843e-4_SRK*xx-3.232035e-3_SRK)*xx+0.01302864_SRK
ELSEIF(xx < 8.20_SRK) THEN
y=(8.701440e-5_SRK*xx-1.524126e-3_SRK)*xx+6.749972e-3_SRK
ELSE !Range 8.2 <= xx <= 9.0, but the > 9.0 portion was already handled
y=(3.742673e-5_SRK*xx-7.157367e-4_SRK)*xx+3.454768e-3_SRK
ENDIF
ENDFUNCTION ThirdOrderBickley
!
ENDMODULE ExtendedMath
src/FunctionTableGenerator.f90 0 → 100644
View file @ a9709d69
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!   Futility Development Group    !
!              All rights reserved.           !
!                         !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license !
! can be found in LICENSE.txt in the head directory of this repository.   !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!> @brief This module defines the @c FunctionTable object and a factory subroutine
!> to generate the table
!>
!> The @c generateFunctionTable function should be called to retrieve a pointer
!> to the new @c FunctionTable. The function to be tabluated is provided as an
!> argument. A parameter list provides such parameters as the minimum and maximum
!> values for which the function should be tabulated, behavior for values outside
!> that range, and generation of the table. The table can be generated using
!> fixed spacing, or it can be dynamically refined until client-specified criteria
!> such as maximum error are achieved.
!>
!> Once the @c FunctionTable pointer is returned to the client, that client is
!> now the owner of the pointer and the memory. This module provides no
!> mechanisms for memory management of the table. The client should call
!> @c evaluate on the table to obtain the value of the function for an input.
!> Once the client is done with the table, it should call @c clear, then
!> deallocate the pointer.
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE FunctionTableGeneratorModule
#include "Futility_DBC.h"
USE Futility_DBC
USE IntrType
USE Constants_Conversion
USE IO_Strings
USE ParameterLists
USE FutilityComputingEnvironmentModule
IMPLICIT NONE
PRIVATE
PUBLIC :: FunctionTable
PUBLIC :: generateFunctionTable
#ifdef UNIT_TEST
PUBLIC :: generateTableData
PUBLIC :: runTestFunctionTableGenerator
INTERFACE
MODULE SUBROUTINE runTestFunctionTableGenerator()
ENDSUBROUTINE runTestFunctionTableGenerator
ENDINTERFACE
#endif
INTERFACE generateFunctionTable
MODULE PROCEDURE generateFunctionTable1D
ENDINTERFACE generateFunctionTable
!> @brief abstract interface for a 1D tabulated function
!> @param x the input value
!> @returns y the output value
ABSTRACT INTERFACE
FUNCTION tabulatedFunction1D(x) RESULT(y)
IMPORT :: SRK
REAL(SRK),INTENT(IN) :: x
REAL(SRK) :: y
ENDFUNCTION tabulatedFunction1D
ENDINTERFACE
!> @brief Stores tabulated data for a function and returns the approximate function value
TYPE :: FunctionTable
PRIVATE
!> Logical indicating initialization status
LOGICAL(SBK) :: isInit=.FALSE.
!> Logical indicating if the table has been populated with data
LOGICAL(SBK) :: hasData=.FALSE.
!> Logical indicating if an error should be thrown for input values greater than @c maxRange
LOGICAL(SBK) :: boundsErrorHigh=.TRUE.
!> Logical indicating if an error should be thrown for input values less than @c minRange
LOGICAL(SBK) :: boundsErrorLow=.TRUE.
!> Number of data points contained in the table
INTEGER(SIK) :: nPoints=0_SIK
!> The minimum input value stored on the table
REAL(SRK) :: minRange=HUGE(ZERO)
!> The maximum input value stored on the table
REAL(SRK) :: maxRange=-HUGE(ZERO)
!> The return value for inputs greater than @c maxRange if @c boundsErrorHigh is false
REAL(SRK) :: highDefaultValue=ZERO
!> The return value for inputs greater than @c minRange if @c boundsErrorLow is false
REAL(SRK) :: lowDefaultValue=ZERO
!> The spacing between input values stored by the table
REAL(SRK) :: spacing=-HUGE(ZERO)
!> The inverse of @c spacing; stored for optimization purposes
REAL(SRK) :: inverseSpacing=-HUGE(ZERO)
!> The list of input values for which the function has been explicitly evaluated
REAL(SRK),ALLOCATABLE :: funcCoordinates(:)
!> The list of function values calculated from @c funcCoordinates and used for interpolation
REAL(SRK),ALLOCATABLE :: funcValues(:)
!> Pointer to the function used for tabulation
PROCEDURE(tabulatedFunction1D),NOPASS,POINTER :: func => NULL()
!> Pointer to the computing environment
CLASS(FutilityComputingEnvironment),POINTER :: ce => NULL()
!
CONTAINS
!> @copybrief FunctionTableGeneratorModule::initFunctionTable1D
!> @copydetails FunctionTableGeneratorModule::initFunctionTable1D
PROCEDURE,PASS :: init => initFunctionTable1D
!> @copybrief FunctionTableGeneratorModule::clearFunctionTable1D
!> @copydetails FunctionTableGeneratorModule::clearFunctionTable1D
PROCEDURE,PASS :: clear => clearFunctionTable1D
!> @copybrief FunctionTableGeneratorModule::evaluateFunctionTable1D
!> @copydetails FunctionTableGeneratorModule::evaluateFunctionTable1D
PROCEDURE,PASS :: evaluate => evaluateFunctionTable1D
ENDTYPE FunctionTable
CHARACTER(LEN=*),PARAMETER :: modName='FunctionTableGenerator'
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
!> @brief generates a @c FunctionTable object and returns a pointer to it
!> @param ce the computing environment
!> @param params the input parameter list controlling the generation of the table
!> @param func the procedure to be tabulated; must have an interface matching @c tabulatedFunction1D
!> @returns table pointer to the initialized function table
!>
FUNCTION generateFunctionTable1D(ce,params,func) RESULT(table)
CLASS(FutilityComputingEnvironment),TARGET,INTENT(INOUT) :: ce
CLASS(ParamType),INTENT(IN) :: params
PROCEDURE(tabulatedFunction1D) :: func
CLASS(FunctionTable),POINTER :: table
!
REAL(SRK) :: maxRelError,maxRelErrorData,maxAbsError,maxAbsErrorData
REQUIRE(params%has('spacing') /= (params%has('maxRelError') .OR. params%has('maxAbsError')))
ALLOCATE(table)
CALL table%init(ce,params,func)
!Now setup the table. If spacing is specified, just call setup once for that spacing
IF(params%has('spacing')) THEN
CALL params%get('spacing',table%spacing)
REQUIRE(table%spacing > ZERO)
CALL generateTableData(table)
!If error is specified, start with 1000 points and refine until the error is small enough
ELSE
IF(params%has('maxRelError')) THEN
CALL params%get('maxRelError',maxRelError)
REQUIRE(maxRelError > ZERO)
ELSE
maxRelError=HUGE(maxRelError)
ENDIF
IF(params%has('maxAbsError')) THEN
CALL params%get('maxAbsError',maxAbsError)
REQUIRE(maxAbsError > ZERO)
ELSE
maxAbsError=HUGE(maxAbsError)
ENDIF
table%spacing = (table%maxRange - table%minRange)*0.001_SRK
maxRelErrorData = HUGE(maxRelError)
maxAbsErrorData = HUGE(maxAbsError)
DO WHILE(maxRelErrorData > maxRelError .OR. maxAbsErrorData > maxAbsError)
CALL generateTableData(table)
CALL evaluateTableQuality(table,maxRelErrorData,maxAbsErrorData)
ENDDO
ENDIF
ENDFUNCTION generateFunctionTable1D
!
!-------------------------------------------------------------------------------
!> @brief Initializes a function table without populating its data
!> @param this the @c FunctionTable to initialize
!> @param ce the computing environment
!> @param params the input parameter list controlling the initialization of the table
!> @param func the procedure to be tabulated; must have an interface matching @c tabulatedFunction1D
!>
SUBROUTINE initFunctionTable1D(this,ce,params,func)
CLASS(FunctionTable),INTENT(INOUT) :: this
CLASS(FutilityComputingEnvironment),TARGET,INTENT(INOUT) :: ce
CLASS(ParamType),INTENT(IN) :: params
PROCEDURE(tabulatedFunction1D) :: func
REQUIRE(params%has('minRange'))
REQUIRE(params%has('maxRange'))
REQUIRE(params%has('boundsErrorHigh'))
REQUIRE(params%has('boundsErrorLow'))
!Set the boundaries and default values (or error behavior) outside the boundaries
CALL params%get('boundsErrorHigh',this%boundsErrorHigh)
CALL params%get('boundsErrorLow',this%boundsErrorLow)
IF(.NOT.this%boundsErrorHigh) THEN
IF(params%has('highDefaultValue')) THEN
CALL params%get('highDefaultValue',this%highDefaultValue)
ENDIF
ENDIF
IF(.NOT.this%boundsErrorLow) THEN
IF(params%has('lowDefaultValue')) THEN
CALL params%get('lowDefaultValue',this%lowDefaultValue)
ENDIF
ENDIF
CALL params%get('minRange',this%minRange)
CALL params%get('maxRange',this%maxRange)
REQUIRE(this%minRange < this%maxRange)
!Associate pointers
this%func => func
this%ce => ce
this%isInit=.TRUE.
this%hasData=.FALSE.
ENDSUBROUTINE initFunctionTable1D
!
!-------------------------------------------------------------------------------
!> @brief Clears a @c FunctionTable object
!> @param this the table to clear
!>
SUBROUTINE clearFunctionTable1D(this)
CLASS(FunctionTable),INTENT(INOUT) :: this
this%isInit = .FALSE.
this%hasData = .FALSE.
this%boundsErrorHigh = .TRUE.
this%boundsErrorLow = .TRUE.
this%nPoints = 0_SIK
this%minRange = HUGE(ZERO)
this%maxRange = -HUGE(ZERO)
this%highDefaultValue = ZERO
this%lowDefaultValue = ZERO
this%spacing = -HUGE(ZERO)
this%inverseSpacing = -HUGE(ZERO)
IF(ALLOCATED(this%funcCoordinates)) DEALLOCATE(this%funcCoordinates)
IF(ALLOCATED(this%funcValues)) DEALLOCATE(this%funcValues)
this%func => NULL()
this%ce => NULL()
ENDSUBROUTINE clearFunctionTable1D
!
!-------------------------------------------------------------------------------
!> @brief Approximately evaluates the tabulated function at @c x
!> @param this the @c FunctionTable object
!> @param x the input value for the function
!> @returns y the approximate function values at @c x
!>
FUNCTION evaluateFunctionTable1D(this,x) RESULT(y)
CLASS(FunctionTable),INTENT(IN) :: this
REAL(SRK),INTENT(IN) :: x
REAL(SRK) :: y
!
CHARACTER(LEN=*),PARAMETER :: myName='evaluateFunctionTable1D'
INTEGER :: tableIndex
REAL(SRK) :: interval, fraction
REQUIRE(this%isInit)
REQUIRE(this%hasData)
REQUIRE(this%inverseSpacing > ZERO)
REQUIRE(ALLOCATED(this%funcCoordinates))
REQUIRE(ALLOCATED(this%funcValues))
interval = x - this%minRange
tableIndex = FLOOR(interval * this%inverseSpacing) + 1
IF(tableIndex < 1) THEN
IF(this%boundsErrorLow) THEN
CALL this%ce%exceptHandler%raiseError(modName//'::'//myName// &
' - Value '//str(x)//' is below the minimum table value '//str(this%minRange)//'!')
ELSE
y = this%lowDefaultValue
ENDIF
ELSEIF(tableIndex >= this%nPoints) THEN !Use >= to avoid index out of bounds in ELSE statement
IF(this%boundsErrorHigh) THEN
CALL this%ce%exceptHandler%raiseError(modName//'::'//myName// &
' - Value '//str(x)//' is above the maximum table value '//str(this%maxRange)//'!')
ELSE
y = this%highDefaultValue
ENDIF
ELSE
fraction = (x - this%funcCoordinates(tableIndex)) * this%inverseSpacing
y = (ONE - fraction) * this%funcValues(tableIndex) + fraction * this%funcValues(tableIndex+1)
ENDIF
ENDFUNCTION evaluateFunctionTable1D
!
!-------------------------------------------------------------------------------
!> @brief Generates the data for a @c FunctionTable object
!> @param table the object to populate data for
!>
!> If @c hasData is true, then it is assumed that the function table data is being
!> refined. The spacing will be cut in half and the data regenerated. If false, then
!> new data will be populated using whatever the current value of @c spacing is.
!>
SUBROUTINE generateTableData(table)
CLASS(FunctionTable),INTENT(INOUT) :: table
!
INTEGER(SIK) :: i,stride
REAL(SRK) :: x
REAL(SRK),ALLOCATABLE :: oldValues(:),oldCoordinates(:)
REQUIRE(ASSOCIATED(table%func))
REQUIRE(table%maxRange > -HUGE(table%maxRange))
REQUIRE(table%minRange < HUGE(table%minRange))
REQUIRE(table%spacing > ZERO)
IF(table%hasData) THEN
!Assume we've already generated the table once and are refining by 2x
IF((SIZE(table%funcValues) > 2) .AND. (MOD(SIZE(table%funcValues),2) == 1)) THEN
CALL MOVE_ALLOC(table%funcValues,oldValues)
CALL MOVE_ALLOC(table%funcCoordinates,oldCoordinates)
stride = 2
table%spacing = table%spacing * HALF
!Generated some table that doesn't conform to how we're going to generate it this time,
! so just throw out the old values
ELSE
DEALLOCATE(table%funcValues)
DEALLOCATE(table%funcCoordinates)
stride = 1
ENDIF
ELSE
stride = 1
ENDIF
table%inverseSpacing = ONE / table%spacing
table%nPoints = NINT((table%maxRange - table%minRange) * table%inverseSpacing) + 1
ALLOCATE(table%funcValues(table%nPoints))
ALLOCATE(table%funcCoordinates(table%nPoints))
IF(stride == 2) THEN
table%funcValues(1:table%nPoints:2)=oldValues(:)
table%funcCoordinates(1:table%nPoints:2)=oldCoordinates(:)
DEALLOCATE(oldValues)
DEALLOCATE(oldCoordinates)
ENDIF
DO i=stride,table%nPoints,stride
x = table%minRange + REAL(i-1,SRK)*table%spacing
table%funcCoordinates(i) = x
table%funcValues(i) = table%func(x)
ENDDO !i
table%hasData=.TRUE.
ENDSUBROUTINE generateTableData
!
!-------------------------------------------------------------------------------
!> @brief Evaluates the quality of a @c FunctionTable compared to the real function
!> @param table the @c FunctionTable to evaluate
!> @returns maxRelError the maximum relative error of the function table
!> @returns maxAbsError the maximum absolute error of the function table
!>
!> If the function returns 0.0 within double precision tolerance, then the value
!> of the table at that input is used to calculate the relative error. This may
!> be of importance when generating the table dynamically.
!>
SUBROUTINE evaluateTableQuality(table,maxRelError,maxAbsError)
CLASS(FunctionTable),INTENT(IN) :: table
REAL(SRK),INTENT(OUT) :: maxRelError
REAL(SRK),INTENT(OUT) :: maxAbsError
!
INTEGER(SIK) :: i
REAL(SRK) :: x,yTab,yRef
REQUIRE(ALLOCATED(table%funcValues))
REQUIRE(ASSOCIATED(table%func))
REQUIRE(table%maxRange > -HUGE(table%maxRange))
REQUIRE(table%minRange < HUGE(table%minRange))
REQUIRE(table%spacing > ZERO)
maxRelError = ZERO
maxAbsError = ZERO
DO i=1,table%nPoints-1
x = table%minRange + (REAL(i-1,SRK)+HALF) * table%spacing
yTab = table%evaluate(x)
yRef = table%func(x)
IF((yRef .APPROXEQA. ZERO) .AND. (.NOT.(yTab .APPROXEQA. ZERO))) THEN
maxRelError = MAX(maxRelError,yTab)
ELSE
maxRelError = MAX(maxRelError,ABS((yTab - yRef)/yRef))
ENDIF
maxAbsError = MAX(maxAbsError,ABS(yTab - yRef))
ENDDO
ENDSUBROUTINE evaluateTableQuality
!
ENDMODULE FunctionTableGeneratorModule
\ No newline at end of file
unit_tests/CMakeLists.txt
View file @ a9709d69
......@@ -76,6 +76,7 @@ SET(UNIT_TEST_NAMES
testSchemaParser
testSpeciesElements
testRSORprecon
testFunctionTableGenerator
)
IF("${PROJECT_NAME}" STREQUAL "Futility" OR "${PROJECT_NAME}" STREQUAL "MPACT")
......
unit_tests/testExtendedMath/testExtendedMath.f90
View file @ a9709d69
......@@ -23,6 +23,7 @@ REGISTER_SUBTEST('GCD',testGCD)
REGISTER_SUBTEST('GCD',testLCM)
REGISTER_SUBTEST('ATAN2PI',testATAN2PI)
REGISTER_SUBTEST('RotateQtrClockwise',testRotateQtrClockwise)
REGISTER_SUBTEST("testThirdOrderBickley", testThirdOrderBickley)
FINALIZE_TEST()
!
......@@ -247,4 +248,72 @@ FUNCTION oldATAN2PI(x,y) RESULT(theta)
IF(theta < 0) theta=2.0_SRK*PI+theta
ENDFUNCTION oldATAN2PI
!
!-------------------------------------------------------------------------------
SUBROUTINE testThirdOrderBickley()
ASSERT(SOFTEQ(ThirdOrderBickley(0.025_SRK),0.760874665500000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.025)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.075_SRK),0.714497185937500_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.075)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.125_SRK),0.671358382812500_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.125)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.175_SRK),0.631150814500000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.175)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.225_SRK),0.593618806687500_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.225)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.275_SRK),0.558542169500000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.275)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.325_SRK),0.525727645750000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.325)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.375_SRK),0.495003064062500_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.375)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.425_SRK),0.466214163562500_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.425)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.475_SRK),0.439221359750000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.475)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.550_SRK),0.401825452250000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.550)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.650_SRK),0.357156291500000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.650)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.750_SRK),0.317712862500000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.750)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.850_SRK),0.282831692750000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.850)")
ASSERT(SOFTEQ(ThirdOrderBickley(0.950_SRK),0.251945412000000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(0.950)")
ASSERT(SOFTEQ(ThirdOrderBickley(1.200_SRK),0.189159171200000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(1.200)")
ASSERT(SOFTEQ(ThirdOrderBickley(1.600_SRK),0.120309080320000_SRK ,1.0e-12_SRK),"ThirdOrderBickley(1.600)")
ASSERT(SOFTEQ(ThirdOrderBickley(2.000_SRK),7.696245999999995E-002_SRK,1.0e-12_SRK),"ThirdOrderBickley(2.000)")
ASSERT(SOFTEQ(ThirdOrderBickley(2.400_SRK),4.945647152000016E-002_SRK,1.0e-12_SRK),"ThirdOrderBickley(2.400)")
ASSERT(SOFTEQ(ThirdOrderBickley(2.800_SRK),3.189697312000006E-002_SRK,1.0e-12_SRK),"ThirdOrderBickley(2.800)")
ASSERT(SOFTEQ(ThirdOrderBickley(3.200_SRK),2.063440492800001E-002_SRK,1.0e-12_SRK),"ThirdOrderBickley(3.200)")
ASSERT(SOFTEQ(ThirdOrderBickley(3.600_SRK),1.338288511999999E-002_SRK,1.0e-12_SRK),"ThirdOrderBickley(3.600)")
ASSERT(SOFTEQ(ThirdOrderBickley(4.000_SRK),8.698716000000009E-003_SRK,1.0e-12_SRK),"ThirdOrderBickley(4.000)")
ASSERT(SOFTEQ(ThirdOrderBickley(4.400_SRK),5.664898399999979E-003_SRK,1.0e-12_SRK),"ThirdOrderBickley(4.400)")
ASSERT(SOFTEQ(ThirdOrderBickley(4.800_SRK),3.695205679999994E-003_SRK,1.0e-12_SRK),"ThirdOrderBickley(4.800)")
ASSERT(SOFTEQ(ThirdOrderBickley(5.400_SRK),1.951856920000009E-003_SRK,1.0e-12_SRK),"ThirdOrderBickley(2.4)")
ASSERT(SOFTEQ(ThirdOrderBickley(6.200_SRK),8.371110279999996E-004_SRK,1.0e-12_SRK),"ThirdOrderBickley(2.8)")
ASSERT(SOFTEQ(ThirdOrderBickley(7.000_SRK),3.604257000000000E-004_SRK,1.0e-12_SRK),"ThirdOrderBickley(3.2)")
ASSERT(SOFTEQ(ThirdOrderBickley(7.800_SRK),1.557452960000002E-004_SRK,1.0e-12_SRK),"ThirdOrderBickley(3.6)")
ASSERT(SOFTEQ(ThirdOrderBickley(8.600_SRK),6.751333079999975E-005_SRK,1.0e-12_SRK),"ThirdOrderBickley(4.0)")
!data from Wolfram Alpha 76 significant figures
!arg x ThirdOrderBickley(x) result of the Bickley function of order three
!0.025 0.76087435435473089481415940877597033571853054821901991516400325474643284773336
!0.075 0.71449679782598526588821547470452465456072719612047630320980195234993342864256
!0.125 0.67135807577109342817946318906174187092787964671445454738446543799021539000978
!0.175 0.63115047862169504452063910448507825276853540965153164809546500702205218310449
!0.225 0.59361841911938945552648920262348309470941047039495859963279516413794004545305
!0.275 0.55854178728401589214049334376515267256572163215350210072437632033519341209916
!0.325 0.52572724070349487181924281894203869507925667656046621747294445654651987447062
!0.375 0.49500276632364780191816272471601602357514824260519750266189586060847625690883
!0.425 0.46621393979863575392216528721577559310239218091302181577479117780150980744750
!0.475 0.43922118686469524800396457044600487137488630880560313204930013463782482792567
!0.550 0.40182535282391698297458635029025432128965486542492333850466147236063725415729
!0.650 0.35715627270140456636611676768402689893215759424244353827223333051537222578432
!0.750 0.31771277656952221958158362255061189690906158712914606465171107896701959742579
!0.850 0.28283162553269040885783128737810736981905570575929112100383119692506037271062
!0.950 0.25194540869425338958124027964104324074752288880999752520768477667293821782424
!1.200 0.18916287763245597978706605859244691850573157613366372399064946356597412542183
!1.600 0.12031089184783858071717638522312464645541380493830736933748211196063082301064
!2.000 0.076963590311658422876818821208590636542634695527976466961577389725761406511519
!2.400 0.049457028031707628064213758059820561481339021775308114578972671841384222775386
!2.800 0.031897511889064199646903466297702116064408073380184568612045684723723587333865
!3.200 0.020634782845562449157993465589641699827465512521288053099219320252709113760499
!3.600 0.013382991329065580861327812619087558102372841626198681240934154676282966267723
!4.000 0.0086987886861997011003542329497674945951430192610029923790805463110000398579947
!4.400 0.0056648909039604982498722037905876778734163941370103855884253919179207346896954
!4.800 0.0036953063069731639358811565454225299312839796120822391929484878605617712940590
!5.400 0.0019521925233853647800369181229233195067620714546474160712559007214062890525572
!6.200 0.00083732548774562206899852981353429892451414336783756347601246685001654830130831
!7.000 0.00036061979247139406228437660036841834297077580646658610363966045026193115652369
!7.800 0.00015584424160148445631661363976906910135930089110821666602687874703788143473585