UniqueUserID General User
Joined: 15 Sep 2009 Posts: 34

Posted: Mon Nov 02, 2009 12:52 am Post subject: Kelvin Functions Ber(x), Bei(x), Ber'(x), Bei'(x) 


If you ever need these functions, this should save you a bit of time:
Code: 
REM ****** Kelvin Functions *********
Function Ber(ByVal x as double) as double
' Calculates Ber(x) function
' Uses recurrence relation based on series expansion formula 820.3
' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"
' 4th Edition, 1961, MacMillan
' Basic code Copyright 2009, Robert Weaver
' Revised: 20091113
' Permission is hereby granted to copy and distribute verbatim copies of this software.
' This software is without warranty. The user accepts all risks incurred from its use.
Dim i As integer,sign As integer
Dim sum As Double,termi As Double,term As Double
if x=0 then
Ber()= 1
else
term=1
sign=1
sum=term*sign
'300 iterations is enough to calculate any value within the range of double precision
for i = 1 to 300
sign=sign
termi=((x/(4*i)) * (x/(4*i2)))^2
term=term*termi
sum=sum+sign*term
'Skip out of loop if current term is less than 1e12 of value
if abs(term/sum)<1e12 then exit for
next
Ber()= sum
end if
end Function
Function Ber_(ByVal x as double) as double
' Calculates Ber'(x) function  d/dx Ber(x)
' Uses recurrence relation based on series expansion formula 820.5
' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"
' 4th Edition, 1961, MacMillan
' Basic code Copyright 2009, Robert Weaver
' Revised: 20091113
' Permission is hereby granted to copy and distribute verbatim copies of this software.
' This software is without warranty. The user accepts all risks incurred from its use.
Dim i As integer,sign As integer
Dim sum As Double,termi As Double,term As Double
if x=0 then
Ber_()= 0
else
term=(x/2)^3/2
sign=1
sum=term*sign
'300 iterations is enough to calculate any value within the range of double precision
for i = 1 to 300
sign=sign
termi=(x/(4*i)) * (x/(4*i+2))^2 * (x/(4*i+4))
term=term*termi
sum=sum+sign*term
'Skip out of loop if current term is less than 1e12 of value
if abs(term/sum)<1e12 then exit for
next
Ber_()= sum
end if
end Function
Function Bei(ByVal x as double) as double
Dim i As integer,sign As integer
Dim sum As Double,termi As Double,term As Double
' Calculates Bei(x) function
' Uses recurrence relation based on series expansion formula 820.4
' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"
' 4th Edition, 1961, MacMillan
' Basic code Copyright 2009, Robert Weaver
' Revised: 20091113
' Permission is hereby granted to copy and distribute verbatim copies of this software.
' This software is without warranty. The user accepts all risks incurred from its use.
if x=0 then
Bei()= 0
else
term=x*x/4
sign=1
sum=term*sign
'300 iterations is enough to calculate any value within the range of double precision
for i = 1 to 300
sign=sign
termi=((x/(4*i)) * (x/(4*i+2)))^2
term=term*termi
sum=sum+sign*term
'Skip out of loop if current term is less than 1e12 of value
if abs(term/sum)<1e12 then exit for
next
Bei()= sum
end if
end Function
Function Bei_(ByVal x as double) as double
Dim i As integer,sign As integer
Dim sum As Double,term As Doublei,term As Double
' Calculates Bei'(x) function  d/dx Bei(x)
' Uses recurrence relation based on series expansion formula 820.6
' from H. B. Dwight, "Tables of Integrals and Other Mathematical Data"
' 4th Edition, 1961, MacMillan
' Basic code Copyright 2009, Robert Weaver
' Revised: 20091113
' Permission is hereby granted to copy and distribute verbatim copies of this software.
' This software is without warranty. The user accepts all risks incurred from its use.
if x=0 then
Bei_()= 0
else
term=x/2
sign=1
sum=term*sign
'300 iterations is enough to calculate any value within the range of double precision
for i = 1 to 300
sign=sign
termi=(x/(4*i2)) * (x/(4*i))^2 * (x/(4*i+2))
term=term*termi
sum=sum+sign*term
'Skip out of loop if current term is less than 1e12 of value
if abs(term/sum)<1e12 then exit for
next
Bei_()= sum
end if
end Function

20091113, Edit Note:
I recently discovered that my original dimension statements
Dim i,sign As integer
Dim sum,termi,term As Double
do not declare ALL variables on the respective line to be all integers or all doubles, just the variable immediately preceding the 'As'.
Therefore I've edited the Dim statements so that all variables are declared correctly.
The failure to declare the floating point variables as doubles didn't seem to have hurt accuracy however. I compared the results of original versions of the functions to published tables and they agree to 10 decimal places. However, it's safer to have things declared properly. 
