Commit 3f56f1f9aabf59955e4227b2895129e1349a51e2
0 parents
initial commit
Showing
46 changed files
with
6531 additions
and
0 deletions
Show diff stats
| 1 | +++ a/CMakeLists.txt | |
| 1 | +#Specify the version being used aswell as the language | |
| 2 | +cmake_minimum_required(VERSION 2.8) | |
| 3 | +#Name your project here | |
| 4 | +project(bimsim) | |
| 5 | + | |
| 6 | +#set the module directory | |
| 7 | +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}") | |
| 8 | + | |
| 9 | +#find BOOST | |
| 10 | +set(Boost_USE_STATIC_LIBS ON) | |
| 11 | +set(Boost_USE_MULTITHREADED ON) | |
| 12 | +set(Boost_USE_STATIC_RUNTIME OFF) | |
| 13 | +find_package( Boost 1.49.0 COMPONENTS program_options ) | |
| 14 | + | |
| 15 | +#find the Qt4 | |
| 16 | +find_package(Qt4 REQUIRED) | |
| 17 | +include_directories(${QT_INCLUDE_DIRECTORY}) | |
| 18 | +include(${QT_USE_FILE}) | |
| 19 | + | |
| 20 | +#set up CUDA | |
| 21 | +find_package(CUDA) | |
| 22 | + | |
| 23 | +#find OpenGL | |
| 24 | +#find_package(OpenGL REQUIRED) | |
| 25 | + | |
| 26 | +#find GLUT | |
| 27 | +#set(GLUT_ROOT_PATH $ENV{GLUT_ROOT_PATH}) | |
| 28 | +#find_package(GLUT REQUIRED) | |
| 29 | + | |
| 30 | +#find GLEW | |
| 31 | +#find_package(GLEW REQUIRED) | |
| 32 | + | |
| 33 | +#add Qt OpenGL stuff | |
| 34 | +#set(QT_USE_QTOPENGL TRUE) | |
| 35 | + | |
| 36 | +#ask the user for the RTS location | |
| 37 | +set(RTS_ROOT_PATH $ENV{RTS_ROOT_PATH}) | |
| 38 | +find_package(RTS REQUIRED) | |
| 39 | + | |
| 40 | +#set the include directories | |
| 41 | +include_directories( | |
| 42 | + ${CMAKE_CURRENT_BINARY_DIR} | |
| 43 | + ${QT_INCLUDES} | |
| 44 | + ${QT_QTOPENGL_INCLUDE_DIR} | |
| 45 | +# ${OPENGL_INCLUDE_DIR} | |
| 46 | +# ${GLEW_INCLUDE_PATH} | |
| 47 | +# ${GLUT_INCLUDE_DIR} | |
| 48 | + ${RTS_INCLUDE_DIR} | |
| 49 | + ${Boost_INCLUDE_DIR} | |
| 50 | +) | |
| 51 | + | |
| 52 | +#enable warnings | |
| 53 | +if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX) | |
| 54 | + add_definitions(-Wall) | |
| 55 | +endif() | |
| 56 | + | |
| 57 | +#Assign source files to the appropriate variables | |
| 58 | +file(GLOB SRC_CPP "*.cpp") | |
| 59 | +file(GLOB SRC_H "*.h") | |
| 60 | +#file(GLOB SRC_UI "*.ui") | |
| 61 | +#file(GLOB SRC_QRC "*.qrc") | |
| 62 | +file(GLOB SRC_CU "*.cu") | |
| 63 | + | |
| 64 | +#set up copying data files | |
| 65 | +#configure_file(rtsSU8_k.txt ${CMAKE_CURRENT_BINARY_DIR}/rtsSU8_k.txt @ONLY) | |
| 66 | +#configure_file(rtsSU8_n.txt ${CMAKE_CURRENT_BINARY_DIR}/rtsSU8_n.txt @ONLY) | |
| 67 | +#configure_file(SurfaceMagnitude.glsl ${CMAKE_CURRENT_BINARY_DIR}/SurfaceMagnitude.glsl @ONLY) | |
| 68 | +#configure_file(SurfaceDeform.glsl ${CMAKE_CURRENT_BINARY_DIR}/SurfaceDeform.glsl @ONLY) | |
| 69 | + | |
| 70 | +#determine which source files have to be moc'd | |
| 71 | +#Qt4_wrap_cpp(UI_MOC ${SRC_H}) | |
| 72 | +#Qt4_wrap_ui(UI_H ${SRC_UI}) | |
| 73 | +#Qt4_add_resources(ALL_RCC ${ALL_QRC}) | |
| 74 | + | |
| 75 | +#moc the necessary files | |
| 76 | +#Qt4_automoc(${ALL_CPP}) | |
| 77 | + | |
| 78 | +#source_group(QtMoc FILES ${UI_MOC}) | |
| 79 | +#source_group(QtUI FILES ${SRC_UI}) | |
| 80 | + | |
| 81 | +#create an executable | |
| 82 | +cuda_add_executable(bimsim | |
| 83 | + ${SRC_CPP} | |
| 84 | + ${SRC_H} | |
| 85 | +# ${UI_H} | |
| 86 | +# ${UI_MOC} | |
| 87 | +# ${ALL_RCC} | |
| 88 | + ${SRC_CU} | |
| 89 | +) | |
| 90 | + | |
| 91 | +#set the link libraries | |
| 92 | +target_link_libraries(bimsim ${QT_LIBRARIES} ${QT_QTOPENGL_LIBRARY} ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ${GLEW_LIBRARY} ${CUDA_cufft_LIBRARY} ${CUDA_cublas_LIBRARY} ${GLUT_glut_LIBRARY} ${Boost_LIBRARIES}) | |
| 93 | + | |
| 94 | + | |
| 95 | + | ... | ... |
| 1 | +++ a/FindCUDASDK.cmake | |
| 1 | +# | |
| 2 | +# The script defines the following variables: | |
| 3 | +# | |
| 4 | +############################################################################## | |
| 5 | +# Note: Removed everything related to CUDA_SDK_ROOT_DIR and only left this as | |
| 6 | +# a possible environment variable to set the SDK directory. | |
| 7 | +# Include file will be: CUDA_CUT_INCLUDE_DIR | |
| 8 | +# Cutil library: CUDA_CUT_LIBRARY | |
| 9 | +############################################################################## | |
| 10 | +# | |
| 11 | +# | |
| 12 | +# CUDA_SDK_ROOT_DIR -- Path to the CUDA SDK. Use this to find files in the | |
| 13 | +# SDK. This script will not directly support finding | |
| 14 | +# specific libraries or headers, as that isn't | |
| 15 | +# supported by NVIDIA. If you want to change | |
| 16 | +# libraries when the path changes see the | |
| 17 | +# FindCUDA.cmake script for an example of how to clear | |
| 18 | +# these variables. There are also examples of how to | |
| 19 | +# use the CUDA_SDK_ROOT_DIR to locate headers or | |
| 20 | +# libraries, if you so choose (at your own risk). | |
| 21 | +# | |
| 22 | +# This code is licensed under the MIT License. See the FindCUDASDK.cmake script | |
| 23 | +# for the text of the license. | |
| 24 | + | |
| 25 | +# The MIT License | |
| 26 | +# | |
| 27 | +# License for the specific language governing rights and limitations under | |
| 28 | +# Permission is hereby granted, free of charge, to any person obtaining a | |
| 29 | +# copy of this software and associated documentation files (the "Software"), | |
| 30 | +# to deal in the Software without restriction, including without limitation | |
| 31 | +# the rights to use, copy, modify, merge, publish, distribute, sublicense, | |
| 32 | +# and/or sell copies of the Software, and to permit persons to whom the | |
| 33 | +# Software is furnished to do so, subject to the following conditions: | |
| 34 | +# | |
| 35 | +# The above copyright notice and this permission notice shall be included | |
| 36 | +# in all copies or substantial portions of the Software. | |
| 37 | +# | |
| 38 | +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS | |
| 39 | +# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | |
| 40 | +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL | |
| 41 | +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | |
| 42 | +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING | |
| 43 | +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER | |
| 44 | +# DEALINGS IN THE SOFTWARE. | |
| 45 | +# | |
| 46 | +############################################################################### | |
| 47 | + | |
| 48 | +# FindCUDASDK.cmake | |
| 49 | + | |
| 50 | +# # Check to see if the CUDA_SDK_ROOT_DIR has changed, | |
| 51 | +# # if it has then clear the cache variable, so that it will be detected again. | |
| 52 | +# if(NOT "${CUDA_SDK_ROOT_DIR}" STREQUAL "${CUDA_SDK_ROOT_DIR_INTERNAL}") | |
| 53 | +# # No specific variables to catch. Use this kind of code before calling | |
| 54 | +# # find_package(CUDA) to clean up any variables that may depend on this path. | |
| 55 | +# | |
| 56 | +# # unset(MY_SPECIAL_CUDA_SDK_INCLUDE_DIR CACHE) | |
| 57 | +# # unset(MY_SPECIAL_CUDA_SDK_LIBRARY CACHE) | |
| 58 | +# endif() | |
| 59 | +# | |
| 60 | +# ######################## | |
| 61 | +# # Look for the SDK stuff | |
| 62 | +# find_path(CUDA_SDK_ROOT_DIR cutil.h | |
| 63 | +# PATH_SUFFIXES "common/inc" "C/common/inc" | |
| 64 | +# "$ENV{NVSDKCUDA_ROOT}" | |
| 65 | +# "[HKEY_LOCAL_MACHINE\\SOFTWARE\\NVIDIA Corporation\\Installed Products\\NVIDIA SDK 10\\Compute;InstallDir]" | |
| 66 | +# "/Developer/GPU\ Computing/C" | |
| 67 | +# ) | |
| 68 | +# | |
| 69 | +# # fallback method for determining CUDA_SDK_ROOT_DIR in case the previous one failed! | |
| 70 | +# if (NOT CUDA_SDK_ROOT_DIR) | |
| 71 | +# find_path(CUDA_SDK_ROOT_DIR C/common/inc/cutil.h | |
| 72 | +# "$ENV{NVSDKCUDA_ROOT}" | |
| 73 | +# "[HKEY_LOCAL_MACHINE\\SOFTWARE\\NVIDIA Corporation\\Installed Products\\NVIDIA SDK 10\\Compute;InstallDir]" | |
| 74 | +# "/Developer/GPU\ Computing/C" | |
| 75 | +# ) | |
| 76 | +# endif() | |
| 77 | + | |
| 78 | +# Keep the CUDA_SDK_ROOT_DIR first in order to be able to override the | |
| 79 | +# environment variables. | |
| 80 | +set(CUDA_SDK_SEARCH_PATH | |
| 81 | + "${CUDA_SDK_ROOT_DIR}" | |
| 82 | + "${CUDA_TOOLKIT_ROOT_DIR}/local/NVSDK0.2" | |
| 83 | + "${CUDA_TOOLKIT_ROOT_DIR}/NVSDK0.2" | |
| 84 | + "${CUDA_TOOLKIT_ROOT_DIR}/NV_SDK" | |
| 85 | + "${CUDA_TOOLKIT_ROOT_DIR}/NV_CUDA_SDK" | |
| 86 | + "$ENV{HOME}/NVIDIA_CUDA_SDK" | |
| 87 | + "$ENV{HOME}/NVIDIA_CUDA_SDK_MACOSX" | |
| 88 | + "$ENV{HOME}/NVIDIA_GPU_Computing_SDK" | |
| 89 | + "/Developer/CUDA" | |
| 90 | + ) | |
| 91 | + | |
| 92 | +# Find include file from the CUDA_SDK_SEARCH_PATH | |
| 93 | + | |
| 94 | +find_path(CUDA_CUT_INCLUDE_DIR | |
| 95 | + cutil.h | |
| 96 | + PATHS ${CUDA_SDK_SEARCH_PATH} | |
| 97 | + PATH_SUFFIXES "common/inc" "C/common/inc" | |
| 98 | + DOC "Location of cutil.h" | |
| 99 | + NO_DEFAULT_PATH | |
| 100 | + ) | |
| 101 | +# Now search system paths | |
| 102 | +find_path(CUDA_CUT_INCLUDE_DIR cutil.h DOC "Location of cutil.h") | |
| 103 | + | |
| 104 | +# mark_as_advanced(CUDA_CUT_INCLUDE_DIR) | |
| 105 | + | |
| 106 | + | |
| 107 | +# Example of how to find a library in the CUDA_SDK_ROOT_DIR | |
| 108 | + | |
| 109 | +# cutil library is called cutil64 for 64 bit builds on windows. We don't want | |
| 110 | +# to get these confused, so we are setting the name based on the word size of | |
| 111 | +# the build. | |
| 112 | + | |
| 113 | +# New library might be called cutil_x86_64 ! | |
| 114 | + | |
| 115 | +if(CMAKE_SIZEOF_VOID_P EQUAL 8) | |
| 116 | + set(cuda_cutil_name cutil64) | |
| 117 | +else(CMAKE_SIZEOF_VOID_P EQUAL 8) | |
| 118 | + set(cuda_cutil_name cutil32) | |
| 119 | +endif(CMAKE_SIZEOF_VOID_P EQUAL 8) | |
| 120 | + | |
| 121 | +find_library(CUDA_CUT_LIBRARY | |
| 122 | + NAMES ${cuda_cutil_name} cutil cutil_x86_64 cutil_i386 | |
| 123 | + PATHS ${CUDA_SDK_SEARCH_PATH} | |
| 124 | + # The new version of the sdk shows up in common/lib, but the old one is in lib | |
| 125 | + # The very newest installation Path of the SDK is in subdirectory 'C'. Please add this Path to the possible suffixes. | |
| 126 | + PATH_SUFFIXES "C/lib" "common/lib" "lib" "C/common/lib" "common/lib" | |
| 127 | + DOC "Location of cutil library" | |
| 128 | + NO_DEFAULT_PATH | |
| 129 | + ) | |
| 130 | +# # Now search system paths | |
| 131 | +# find_library(CUDA_CUT_LIBRARY NAMES cutil ${cuda_cutil_name} DOC "Location of cutil library") | |
| 132 | +mark_as_advanced(CUDA_CUT_LIBRARY) | |
| 133 | +set(CUDA_CUT_LIBRARIES ${CUDA_CUT_LIBRARY}) | |
| 134 | + | |
| 135 | +############################# | |
| 136 | +# Check for required components | |
| 137 | +if(CUDA_CUT_INCLUDE_DIR) | |
| 138 | + set(CUDASDK_FOUND TRUE) | |
| 139 | +endif(CUDA_CUT_INCLUDE_DIR) | |
| 140 | + | |
| 141 | +# set(CUDA_SDK_ROOT_DIR_INTERNAL "${CUDA_SDK_ROOT_DIR}" CACHE INTERNAL | |
| 142 | +# "This is the value of the last time CUDA_SDK_ROOT_DIR was set successfully." FORCE) | ... | ... |
| 1 | +++ a/FindGLEW.cmake | |
| 1 | +# | |
| 2 | +# Try to find GLEW library and include path. | |
| 3 | +# Once done this will define | |
| 4 | +# | |
| 5 | +# GLEW_FOUND | |
| 6 | +# GLEW_INCLUDE_PATH | |
| 7 | +# GLEW_LIBRARY | |
| 8 | +# | |
| 9 | + | |
| 10 | +IF (WIN32) | |
| 11 | + FIND_PATH( GLEW_INCLUDE_PATH GL/glew.h | |
| 12 | + $ENV{PROGRAMFILES}/GLEW/include | |
| 13 | + ${PROJECT_SOURCE_DIR}/src/nvgl/glew/include | |
| 14 | + DOC "The directory where GL/glew.h resides") | |
| 15 | + FIND_LIBRARY( GLEW_LIBRARY | |
| 16 | + NAMES glew GLEW glew32 glew32s | |
| 17 | + PATHS | |
| 18 | + $ENV{PROGRAMFILES}/GLEW/lib | |
| 19 | + ${PROJECT_SOURCE_DIR}/src/nvgl/glew/bin | |
| 20 | + ${PROJECT_SOURCE_DIR}/src/nvgl/glew/lib | |
| 21 | + DOC "The GLEW library") | |
| 22 | +ELSE (WIN32) | |
| 23 | + FIND_PATH( GLEW_INCLUDE_PATH GL/glew.h | |
| 24 | + /usr/include | |
| 25 | + /usr/local/include | |
| 26 | + /sw/include | |
| 27 | + /opt/local/include | |
| 28 | + DOC "The directory where GL/glew.h resides") | |
| 29 | + FIND_LIBRARY( GLEW_LIBRARY | |
| 30 | + NAMES GLEW glew | |
| 31 | + PATHS | |
| 32 | + /usr/lib64 | |
| 33 | + /usr/lib | |
| 34 | + /usr/local/lib64 | |
| 35 | + /usr/local/lib | |
| 36 | + /sw/lib | |
| 37 | + /opt/local/lib | |
| 38 | + DOC "The GLEW library") | |
| 39 | +ENDIF (WIN32) | |
| 40 | + | |
| 41 | +IF (GLEW_INCLUDE_PATH) | |
| 42 | + SET( GLEW_FOUND 1 CACHE STRING "Set to 1 if GLEW is found, 0 otherwise") | |
| 43 | +ELSE (GLEW_INCLUDE_PATH) | |
| 44 | + SET( GLEW_FOUND 0 CACHE STRING "Set to 1 if GLEW is found, 0 otherwise") | |
| 45 | +ENDIF (GLEW_INCLUDE_PATH) | |
| 46 | + | |
| 47 | +MARK_AS_ADVANCED( | |
| 48 | + GLEW_FOUND | |
| 49 | + GLEW_INCLUDE_PATH | |
| 50 | + GLEW_LIBRARY | |
| 51 | +) | |
| 0 | 52 | \ No newline at end of file | ... | ... |
| 1 | +++ a/FindRTS.cmake | |
| 1 | +# Tries to find the RTS include directory | |
| 2 | + | |
| 3 | + FIND_PATH( RTS_INCLUDE_DIR NAMES rts_glShaderProgram.h | |
| 4 | + PATHS | |
| 5 | + ${CMAKE_CURRENT_SOURCE_DIR}/rts | |
| 6 | + ${RTS_ROOT_PATH} | |
| 7 | +) | |
| 8 | + | |
| 9 | +IF (RTS_FOUND) | |
| 10 | + #The following deprecated settings are for backwards compatibility with CMake1.4 | |
| 11 | + SET (RTS_INCLUDE_PATH ${RTS_INCLUDE_DIR}) | |
| 12 | +ENDIF(RTS_FOUND) | |
| 13 | + | |
| 14 | +FIND_PACKAGE_HANDLE_STANDARD_ARGS(RTS REQUIRED_VARS TRUE RTS_INCLUDE_DIR) | ... | ... |
| 1 | +++ a/README.md | |
| 1 | +BIM-Sim computes the interactions between an incident electromagnetic field and spherical scatterers using Mie theory. | |
| 2 | + | |
| 3 | +A Windows 32 executable is provided in the Win32 directory. All of the necessary runtime libraries are included. Start by running 'ricalc -h' for command-line options. | |
| 4 | + | |
| 5 | +This software is designed to be used with CMake, so we recommend using it to build the application. | |
| 6 | + | |
| 7 | +This software requires the following libraries for compiling: | |
| 8 | + | |
| 9 | +CUDA (including CUFFT and CUBLAS) - NVidia's GPU programming toolkit -- https://developer.nvidia.com/ | |
| 10 | + | |
| 11 | +RTS - Real-Time Scraps C++ library (my codebase) -- https://github.com/dmayerich/RTS | |
| 12 | + | |
| 13 | +Boost C++ Library (specifically the Program Options library) -- http://www.boost.org/ | |
| 14 | + | |
| 15 | +Qt (Nokia's GUI Library) -- https://qt-project.org/ | |
| 16 | + | |
| 17 | + | |
| 18 | + | |
| 19 | + | |
| 20 | + | |
| 21 | + | ... | ... |
| 1 | +++ a/bessel.h | |
| 1 | +#ifndef bessH | |
| 2 | +#define bessH | |
| 3 | +#define _USE_MATH_DEFINES | |
| 4 | +#include <math.h> | |
| 5 | +#include <complex> | |
| 6 | +using namespace std; | |
| 7 | +#define eps 1e-15 | |
| 8 | +#define el 0.5772156649015329 | |
| 9 | + | |
| 10 | +int msta1(double x,int mp); | |
| 11 | +int msta2(double x,int n,int mp); | |
| 12 | +int bessjy01a(double x,double &j0,double &j1,double &y0,double &y1, | |
| 13 | + double &j0p,double &j1p,double &y0p,double &y1p); | |
| 14 | +int bessjy01b(double x,double &j0,double &j1,double &y0,double &y1, | |
| 15 | + double &j0p,double &j1p,double &y0p,double &y1p); | |
| 16 | +int bessjyna(int n,double x,int &nm,double *jn,double *yn, | |
| 17 | + double *jnp,double *ynp); | |
| 18 | +int bessjynb(int n,double x,int &nm,double *jn,double *yn, | |
| 19 | + double *jnp,double *ynp); | |
| 20 | +int bessjyv(double v,double x,double &vm,double *jv,double *yv, | |
| 21 | + double *jvp,double *yvp); | |
| 22 | +int bessik01a(double x,double &i0,double &i1,double &k0,double &k1, | |
| 23 | + double &i0p,double &i1p,double &k0p,double &k1p); | |
| 24 | +int bessik01b(double x,double &i0,double &i1,double &k0,double &k1, | |
| 25 | + double &i0p,double &i1p,double &k0p,double &k1p); | |
| 26 | +int bessikna(int n,double x,int &nm,double *in,double *kn, | |
| 27 | + double *inp,double *knp); | |
| 28 | +int bessiknb(int n,double x,int &nm,double *in,double *kn, | |
| 29 | + double *inp,double *knp); | |
| 30 | +int bessikv(double v,double x,double &vm,double *iv,double *kv, | |
| 31 | + double *ivp,double *kvp); | |
| 32 | +int cbessjy01(complex<double> z,complex<double> &cj0,complex<double> &cj1, | |
| 33 | + complex<double> &cy0,complex<double> &cy1,complex<double> &cj0p, | |
| 34 | + complex<double> &cj1p,complex<double> &cy0p,complex<double> &cy1p); | |
| 35 | +int cbessjyna(int n,complex<double> z,int &nm,complex<double> *cj, | |
| 36 | + complex<double> *cy,complex<double> *cjp,complex<double> *cyp); | |
| 37 | +int cbessjynb(int n,complex<double> z,int &nm,complex<double> *cj, | |
| 38 | + complex<double> *cy,complex<double> *cjp,complex<double> *cyp); | |
| 39 | +int cbessik01(complex<double>z,complex<double>&ci0,complex<double>&ci1, | |
| 40 | + complex<double>&ck0,complex<double>&ck1,complex<double>&ci0p, | |
| 41 | + complex<double>&ci1p,complex<double>&ck0p,complex<double>&ck1p); | |
| 42 | +int cbessikna(int n,complex<double> z,int &nm,complex<double> *ci, | |
| 43 | + complex<double> *ck,complex<double> *cip,complex<double> *ckp); | |
| 44 | +int cbessiknb(int n,complex<double> z,int &nm,complex<double> *ci, | |
| 45 | + complex<double> *ck,complex<double> *cip,complex<double> *ckp); | |
| 46 | +int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv, | |
| 47 | + complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp); | |
| 48 | +int cbessikv(double v,complex<double>z,double &vm,complex<double> *civ, | |
| 49 | + complex<double> *ckv,complex<double> *civp,complex<double> *ckvp); | |
| 50 | + | |
| 51 | +#endif | ... | ... |
| 1 | +++ a/bessik.cpp | |
| 1 | +// bessik.cpp -- computation of modified Bessel functions In, Kn | |
| 2 | +// and their derivatives. Algorithms and coefficient values from | |
| 3 | +// "Computation of Special Functions", Zhang and Jin, John | |
| 4 | +// Wiley and Sons, 1996. | |
| 5 | +// | |
| 6 | +// (C) 2003, C. Bond. All rights reserved. | |
| 7 | +// | |
| 8 | +#define _USE_MATH_DEFINES | |
| 9 | +#include <math.h> | |
| 10 | +#include "bessel.h" | |
| 11 | + | |
| 12 | +double gamma(double x); | |
| 13 | + | |
| 14 | +int bessik01a(double x,double &i0,double &i1,double &k0,double &k1, | |
| 15 | + double &i0p,double &i1p,double &k0p,double &k1p) | |
| 16 | +{ | |
| 17 | + double r,x2,ca,cb,ct,ww,w0,xr,xr2; | |
| 18 | + int k,kz; | |
| 19 | + static double a[] = { | |
| 20 | + 0.125, | |
| 21 | + 7.03125e-2, | |
| 22 | + 7.32421875e-2, | |
| 23 | + 1.1215209960938e-1, | |
| 24 | + 2.2710800170898e-1, | |
| 25 | + 5.7250142097473e-1, | |
| 26 | + 1.7277275025845, | |
| 27 | + 6.0740420012735, | |
| 28 | + 2.4380529699556e1, | |
| 29 | + 1.1001714026925e2, | |
| 30 | + 5.5133589612202e2, | |
| 31 | + 3.0380905109224e3}; | |
| 32 | + static double b[] = { | |
| 33 | + -0.375, | |
| 34 | + -1.171875e-1, | |
| 35 | + -1.025390625e-1, | |
| 36 | + -1.4419555664063e-1, | |
| 37 | + -2.7757644653320e-1, | |
| 38 | + -6.7659258842468e-1, | |
| 39 | + -1.9935317337513, | |
| 40 | + -6.8839142681099, | |
| 41 | + -2.7248827311269e1, | |
| 42 | + -1.2159789187654e2, | |
| 43 | + -6.0384407670507e2, | |
| 44 | + -3.3022722944809e3}; | |
| 45 | + static double a1[] = { | |
| 46 | + 0.125, | |
| 47 | + 0.2109375, | |
| 48 | + 1.0986328125, | |
| 49 | + 1.1775970458984e1, | |
| 50 | + 2.1461706161499e2, | |
| 51 | + 5.9511522710323e3, | |
| 52 | + 2.3347645606175e5, | |
| 53 | + 1.2312234987631e7}; | |
| 54 | + | |
| 55 | + if (x < 0.0) return 1; | |
| 56 | + if (x == 0.0) { | |
| 57 | + i0 = 1.0; | |
| 58 | + i1 = 0.0; | |
| 59 | + k0 = 1e308; | |
| 60 | + k1 = 1e308; | |
| 61 | + i0p = 0.0; | |
| 62 | + i1p = 0.5; | |
| 63 | + k0p = -1e308; | |
| 64 | + k1p = -1e308; | |
| 65 | + return 0; | |
| 66 | + } | |
| 67 | + x2 = x*x; | |
| 68 | + if (x <= 18.0) { | |
| 69 | + i0 = 1.0; | |
| 70 | + r = 1.0; | |
| 71 | + for (k=1;k<=50;k++) { | |
| 72 | + r *= 0.25*x2/(k*k); | |
| 73 | + i0 += r; | |
| 74 | + if (fabs(r/i0) < eps) break; | |
| 75 | + } | |
| 76 | + i1 = 1.0; | |
| 77 | + r = 1.0; | |
| 78 | + for (k=1;k<=50;k++) { | |
| 79 | + r *= 0.25*x2/(k*(k+1)); | |
| 80 | + i1 += r; | |
| 81 | + if (fabs(r/i1) < eps) break; | |
| 82 | + } | |
| 83 | + i1 *= 0.5*x; | |
| 84 | + } | |
| 85 | + else { | |
| 86 | + if (x >= 50.0) kz = 7; | |
| 87 | + else if (x >= 35.0) kz = 9; | |
| 88 | + else kz = 12; | |
| 89 | + ca = exp(x)/sqrt(2.0*M_PI*x); | |
| 90 | + i0 = 1.0; | |
| 91 | + xr = 1.0/x; | |
| 92 | + for (k=0;k<kz;k++) { | |
| 93 | + i0 += a[k]*pow(xr,k+1); | |
| 94 | + } | |
| 95 | + i0 *= ca; | |
| 96 | + i1 = 1.0; | |
| 97 | + for (k=0;k<kz;k++) { | |
| 98 | + i1 += b[k]*pow(xr,k+1); | |
| 99 | + } | |
| 100 | + i1 *= ca; | |
| 101 | + } | |
| 102 | + if (x <= 9.0) { | |
| 103 | + ct = -(log(0.5*x)+el); | |
| 104 | + k0 = 0.0; | |
| 105 | + w0 = 0.0; | |
| 106 | + r = 1.0; | |
| 107 | + ww = 0.0; | |
| 108 | + for (k=1;k<=50;k++) { | |
| 109 | + w0 += 1.0/k; | |
| 110 | + r *= 0.25*x2/(k*k); | |
| 111 | + k0 += r*(w0+ct); | |
| 112 | + if (fabs((k0-ww)/k0) < eps) break; | |
| 113 | + ww = k0; | |
| 114 | + } | |
| 115 | + k0 += ct; | |
| 116 | + } | |
| 117 | + else { | |
| 118 | + cb = 0.5/x; | |
| 119 | + xr2 = 1.0/x2; | |
| 120 | + k0 = 1.0; | |
| 121 | + for (k=0;k<8;k++) { | |
| 122 | + k0 += a1[k]*pow(xr2,k+1); | |
| 123 | + } | |
| 124 | + k0 *= cb/i0; | |
| 125 | + } | |
| 126 | + k1 = (1.0/x - i1*k0)/i0; | |
| 127 | + i0p = i1; | |
| 128 | + i1p = i0-i1/x; | |
| 129 | + k0p = -k1; | |
| 130 | + k1p = -k0-k1/x; | |
| 131 | + return 0; | |
| 132 | +} | |
| 133 | + | |
| 134 | +int bessik01b(double x,double &i0,double &i1,double &k0,double &k1, | |
| 135 | + double &i0p,double &i1p,double &k0p,double &k1p) | |
| 136 | +{ | |
| 137 | + double t,t2,dtmp,dtmp1; | |
| 138 | + | |
| 139 | + if (x < 0.0) return 1; | |
| 140 | + if (x == 0.0) { | |
| 141 | + i0 = 1.0; | |
| 142 | + i1 = 0.0; | |
| 143 | + k0 = 1e308; | |
| 144 | + k1 = 1e308; | |
| 145 | + i0p = 0.0; | |
| 146 | + i1p = 0.5; | |
| 147 | + k0p = -1e308; | |
| 148 | + k1p = -1e308; | |
| 149 | + return 0; | |
| 150 | + } | |
| 151 | + if (x < 3.75) { | |
| 152 | + t = x/3.75; | |
| 153 | + t2 = t*t; | |
| 154 | + i0 = (((((0.0045813*t2+0.0360768)*t2+0.2659732)*t2+ | |
| 155 | + 1.2067492)*t2+3.0899424)*t2+3.5156229)*t2+1.0; | |
| 156 | + i1 = x*(((((0.00032411*t2+0.00301532)*t2+0.02658733*t2+ | |
| 157 | + 0.15084934)*t2+0.51498869)*t2+0.87890594)*t2+0.5); | |
| 158 | + } | |
| 159 | + else { | |
| 160 | + t = 3.75/x; | |
| 161 | + dtmp1 = exp(x)/sqrt(x); | |
| 162 | + dtmp = (((((((0.00392377*t-0.01647633)*t+0.026355537)*t-0.02057706)*t+ | |
| 163 | + 0.00916281)*t-0.00157565)*t+0.00225319)*t+0.01328592)*t+0.39894228; | |
| 164 | + i0 = dtmp*dtmp1; | |
| 165 | + dtmp = (((((((-0.00420059*t+0.01787654)*t-0.02895312)*t+0.02282967)*t- | |
| 166 | + 0.01031555)*t+0.00163801)*t-0.00362018)*t-0.03988024)*t+0.39894228; | |
| 167 | + i1 = dtmp*dtmp1; | |
| 168 | + } | |
| 169 | + if (x < 2.0) { | |
| 170 | + t = 0.5*x; | |
| 171 | + t2 = t*t; // already calculated above | |
| 172 | + dtmp = (((((0.0000074*t2+0.0001075)*t2+0.00262698)*t2+0.0348859)*t2+ | |
| 173 | + 0.23069756)*t2+0.4227842)*t2-0.57721566; | |
| 174 | + k0 = dtmp - i0*log(t); | |
| 175 | + dtmp = (((((-0.00004686*t2-0.00110404)*t2-0.01919402)*t2- | |
| 176 | + 0.18156897)*t2-0.67278578)*t2+0.15443144)*t2+1.0; | |
| 177 | + k1 = dtmp/x + i1*log(t); | |
| 178 | + } | |
| 179 | + else { | |
| 180 | + t = 2.0/x; | |
| 181 | + dtmp1 = exp(-x)/sqrt(x); | |
| 182 | + dtmp = (((((0.00053208*t-0.0025154)*t+0.00587872)*t-0.01062446)*t+ | |
| 183 | + 0.02189568)*t-0.07832358)*t+1.25331414; | |
| 184 | + k0 = dtmp*dtmp1; | |
| 185 | + dtmp = (((((-0.00068245*t+0.00325614)*t-0.00780353)*t+0.01504268)*t- | |
| 186 | + 0.0365562)*t+0.23498619)*t+1.25331414; | |
| 187 | + k1 = dtmp*dtmp1; | |
| 188 | + } | |
| 189 | + i0p = i1; | |
| 190 | + i1p = i0 - i1/x; | |
| 191 | + k0p = -k1; | |
| 192 | + k1p = -k0 - k1/x; | |
| 193 | + return 0; | |
| 194 | +} | |
| 195 | +int bessikna(int n,double x,int &nm,double *in,double *kn, | |
| 196 | + double *inp,double *knp) | |
| 197 | +{ | |
| 198 | + double bi0,bi1,bk0,bk1,g,g0,g1,f,f0,f1,h,h0,h1,s0; | |
| 199 | + int k,m,ecode; | |
| 200 | + | |
| 201 | + if ((x < 0.0) || (n < 0)) return 1; | |
| 202 | + if (x < eps) { | |
| 203 | + for (k=0;k<=n;k++) { | |
| 204 | + in[k] = 0.0; | |
| 205 | + kn[k] = 1e308; | |
| 206 | + inp[k] = 0.0; | |
| 207 | + knp[k] = -1e308; | |
| 208 | + } | |
| 209 | + in[0] = 1.0; | |
| 210 | + inp[1] = 0.5; | |
| 211 | + return 0; | |
| 212 | + } | |
| 213 | + nm = n; | |
| 214 | + ecode = bessik01a(x,in[0],in[1],kn[0],kn[1],inp[0],inp[1],knp[0],knp[1]); | |
| 215 | + if (n < 2) return 0; | |
| 216 | + bi0 = in[0]; | |
| 217 | + bi1 = in[1]; | |
| 218 | + bk0 = kn[0]; | |
| 219 | + bk1 = kn[1]; | |
| 220 | + if ((x > 40.0) && (n < (int)(0.25*x))) { | |
| 221 | + h0 = bi0; | |
| 222 | + h1 = bi1; | |
| 223 | + for (k=2;k<=n;k++) { | |
| 224 | + h = -2.0*(k-1.0)*h1/x+h0; | |
| 225 | + in[k] = h; | |
| 226 | + h0 = h1; | |
| 227 | + h1 = h; | |
| 228 | + } | |
| 229 | + } | |
| 230 | + else { | |
| 231 | + m = msta1(x,200); | |
| 232 | + if (m < n) nm = m; | |
| 233 | + else m = msta2(x,n,15); | |
| 234 | + f0 = 0.0; | |
| 235 | + f1 = 1.0e-100; | |
| 236 | + for (k=m;k>=0;k--) { | |
| 237 | + f = 2.0*(k+1.0)*f1/x+f0; | |
| 238 | + if (x <= nm) in[k] = f; | |
| 239 | + f0 = f1; | |
| 240 | + f1 = f; | |
| 241 | + } | |
| 242 | + s0 = bi0/f; | |
| 243 | + for (k=0;k<=m;k++) { | |
| 244 | + in[k] *= s0; | |
| 245 | + } | |
| 246 | + } | |
| 247 | + g0 = bk0; | |
| 248 | + g1 = bk1; | |
| 249 | + for (k=2;k<=nm;k++) { | |
| 250 | + g = 2.0*(k-1.0)*g1/x+g0; | |
| 251 | + kn[k] = g; | |
| 252 | + g0 = g1; | |
| 253 | + g1 = g; | |
| 254 | + } | |
| 255 | + for (k=2;k<=nm;k++) { | |
| 256 | + inp[k] = in[k-1]-k*in[k]/x; | |
| 257 | + knp[k] = -kn[k-1]-k*kn[k]/x; | |
| 258 | + } | |
| 259 | + return 0; | |
| 260 | +} | |
| 261 | +int bessiknb(int n,double x,int &nm,double *in,double *kn, | |
| 262 | + double *inp,double *knp) | |
| 263 | +{ | |
| 264 | + double s0,bs,f,f0,f1,sk0,a0,bkl,vt,r,g,g0,g1; | |
| 265 | + int k,kz,m,l; | |
| 266 | + | |
| 267 | + if ((x < 0.0) || (n < 0)) return 1; | |
| 268 | + if (x < eps) { | |
| 269 | + for (k=0;k<=n;k++) { | |
| 270 | + in[k] = 0.0; | |
| 271 | + kn[k] = 1e308; | |
| 272 | + inp[k] = 0.0; | |
| 273 | + knp[k] = -1e308; | |
| 274 | + } | |
| 275 | + in[0] = 1.0; | |
| 276 | + inp[1] = 0.5; | |
| 277 | + return 0; | |
| 278 | + } | |
| 279 | + nm = n; | |
| 280 | + if (n == 0) nm = 1; | |
| 281 | + m = msta1(x,200); | |
| 282 | + if (m < nm) nm = m; | |
| 283 | + else m = msta2(x,nm,15); | |
| 284 | + bs = 0.0; | |
| 285 | + sk0 = 0.0; | |
| 286 | + f0 = 0.0; | |
| 287 | + f1 = 1.0e-100; | |
| 288 | + for (k=m;k>=0;k--) { | |
| 289 | + f = 2.0*(k+1.0)*f1/x+f0; | |
| 290 | + if (k <= nm) in[k] = f; | |
| 291 | + if ((k != 0) && (k == 2*(int)(k/2))) { | |
| 292 | + sk0 += 4.0*f/k; | |
| 293 | + } | |
| 294 | + bs += 2.0*f; | |
| 295 | + f0 = f1; | |
| 296 | + f1 = f; | |
| 297 | + } | |
| 298 | + s0 = exp(x)/(bs-f); | |
| 299 | + for (k=0;k<=nm;k++) { | |
| 300 | + in[k] *= s0; | |
| 301 | + } | |
| 302 | + if (x <= 8.0) { | |
| 303 | + kn[0] = -(log(0.5*x)+el)*in[0]+s0*sk0; | |
| 304 | + kn[1] = (1.0/x-in[1]*kn[0])/in[0]; | |
| 305 | + } | |
| 306 | + else { | |
| 307 | + a0 = sqrt(M_PI_2/x)*exp(-x); | |
| 308 | + if (x >= 200.0) kz = 6; | |
| 309 | + else if (x >= 80.0) kz = 8; | |
| 310 | + else if (x >= 25.0) kz = 10; | |
| 311 | + else kz = 16; | |
| 312 | + for (l=0;l<2;l++) { | |
| 313 | + bkl = 1.0; | |
| 314 | + vt = 4.0*l; | |
| 315 | + r = 1.0; | |
| 316 | + for (k=1;k<=kz;k++) { | |
| 317 | + r *= 0.125*(vt-pow(2.0*k-1.0,2))/(k*x); | |
| 318 | + bkl += r; | |
| 319 | + } | |
| 320 | + kn[l] = a0*bkl; | |
| 321 | + } | |
| 322 | + } | |
| 323 | + g0 = kn[0]; | |
| 324 | + g1 = kn[1]; | |
| 325 | + for (k=2;k<=nm;k++) { | |
| 326 | + g = 2.0*(k-1.0)*g1/x+g0; | |
| 327 | + kn[k] = g; | |
| 328 | + g0 = g1; | |
| 329 | + g1 = g; | |
| 330 | + } | |
| 331 | + inp[0] = in[1]; | |
| 332 | + knp[0] = -kn[1]; | |
| 333 | + for (k=1;k<=nm;k++) { | |
| 334 | + inp[k] = in[k-1]-k*in[k]/x; | |
| 335 | + knp[k] = -kn[k-1]-k*kn[k]/x; | |
| 336 | + } | |
| 337 | + return 0; | |
| 338 | +} | |
| 339 | + | |
| 340 | +// The following program computes the modified Bessel functions | |
| 341 | +// Iv(x) and Kv(x) for arbitrary positive order. For negative | |
| 342 | +// order use: | |
| 343 | +// | |
| 344 | +// I-v(x) = Iv(x) + 2/pi sin(v pi) Kv(x) | |
| 345 | +// K-v(x) = Kv(x) | |
| 346 | +// | |
| 347 | +int bessikv(double v,double x,double &vm,double *iv,double *kv, | |
| 348 | + double *ivp,double *kvp) | |
| 349 | +{ | |
| 350 | + double x2,v0,piv,vt,a1,v0p,gap,r,bi0,ca,sum; | |
| 351 | + double f,f1,f2,ct,cs,wa,gan,ww,w0,v0n; | |
| 352 | + double r1,r2,bk0,bk1,bk2,a2,cb; | |
| 353 | + int n,k,kz,m; | |
| 354 | + | |
| 355 | + if ((v < 0.0) || (x < 0.0)) return 1; | |
| 356 | + x2 = x*x; | |
| 357 | + n = (int)v; | |
| 358 | + v0 = v-n; | |
| 359 | + if (n == 0) n = 1; | |
| 360 | + if (x == 0.0) { | |
| 361 | + for (k=0;k<=n;k++) { | |
| 362 | + iv[k] = 0.0; | |
| 363 | + kv[k] = -1e308; | |
| 364 | + ivp[k] = 0.0; | |
| 365 | + kvp[k] = 1e308; | |
| 366 | + } | |
| 367 | + if (v0 == 0.0) { | |
| 368 | + iv[0] = 1.0; | |
| 369 | + ivp[1] = 0.5; | |
| 370 | + } | |
| 371 | + vm = v; | |
| 372 | + return 0; | |
| 373 | + } | |
| 374 | + piv = M_PI*v0; | |
| 375 | + vt = 4.0*v0*v0; | |
| 376 | + if (v0 == 0.0) { | |
| 377 | + a1 = 1.0; | |
| 378 | + } | |
| 379 | + else { | |
| 380 | + v0p = 1.0+v0; | |
| 381 | + gap = gamma(v0p); | |
| 382 | + a1 = pow(0.5*x,v0)/gap; | |
| 383 | + } | |
| 384 | + if (x >= 50.0) kz = 8; | |
| 385 | + else if (x >= 35.0) kz = 10; | |
| 386 | + else kz = 14; | |
| 387 | + if (x <= 18.0) { | |
| 388 | + bi0 = 1.0; | |
| 389 | + r = 1.0; | |
| 390 | + for (k=1;k<=30;k++) { | |
| 391 | + r *= 0.25*x2/(k*(k+v0)); | |
| 392 | + bi0 += r; | |
| 393 | + if (fabs(r/bi0) < eps) break; | |
| 394 | + } | |
| 395 | + bi0 *= a1; | |
| 396 | + } | |
| 397 | + else { | |
| 398 | + ca = exp(x)/sqrt(2.0*M_PI*x); | |
| 399 | + sum = 1.0; | |
| 400 | + r = 1.0; | |
| 401 | + for (k=1;k<=kz;k++) { | |
| 402 | + r *= -0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x); | |
| 403 | + sum += r; | |
| 404 | + } | |
| 405 | + bi0 = ca*sum; | |
| 406 | + } | |
| 407 | + m = msta1(x,200); | |
| 408 | + if (m < n) n = m; | |
| 409 | + else m = msta2(x,n,15); | |
| 410 | + f2 = 0.0; | |
| 411 | + f1 = 1.0e-100; | |
| 412 | + for (k=m;k>=0;k--) { | |
| 413 | + f = 2.0*(v0+k+1.0)*f1/x+f2; | |
| 414 | + if (k <= n) iv[k] = f; | |
| 415 | + f2 = f1; | |
| 416 | + f1 = f; | |
| 417 | + } | |
| 418 | + cs = bi0/f; | |
| 419 | + for (k=0;k<=n;k++) { | |
| 420 | + iv[k] *= cs; | |
| 421 | + } | |
| 422 | + ivp[0] = v0*iv[0]/x+iv[1]; | |
| 423 | + for (k=1;k<=n;k++) { | |
| 424 | + ivp[k] = -(k+v0)*iv[k]/x+iv[k-1]; | |
| 425 | + } | |
| 426 | + ww = 0.0; | |
| 427 | + if (x <= 9.0) { | |
| 428 | + if (v0 == 0.0) { | |
| 429 | + ct = -log(0.5*x)-el; | |
| 430 | + cs = 0.0; | |
| 431 | + w0 = 0.0; | |
| 432 | + r = 1.0; | |
| 433 | + for (k=1;k<=50;k++) { | |
| 434 | + w0 += 1.0/k; | |
| 435 | + r *= 0.25*x2/(k*k); | |
| 436 | + cs += r*(w0+ct); | |
| 437 | + wa = fabs(cs); | |
| 438 | + if (fabs((wa-ww)/wa) < eps) break; | |
| 439 | + ww = wa; | |
| 440 | + } | |
| 441 | + bk0 = ct+cs; | |
| 442 | + } | |
| 443 | + else { | |
| 444 | + v0n = 1.0-v0; | |
| 445 | + gan = gamma(v0n); | |
| 446 | + a2 = 1.0/(gan*pow(0.5*x,v0)); | |
| 447 | + a1 = pow(0.5*x,v0)/gap; | |
| 448 | + sum = a2-a1; | |
| 449 | + r1 = 1.0; | |
| 450 | + r2 = 1.0; | |
| 451 | + for (k=1;k<=120;k++) { | |
| 452 | + r1 *= 0.25*x2/(k*(k-v0)); | |
| 453 | + r2 *= 0.25*x2/(k*(k+v0)); | |
| 454 | + sum += a2*r1-a1*r2; | |
| 455 | + wa = fabs(sum); | |
| 456 | + if (fabs((wa-ww)/wa) < eps) break; | |
| 457 | + ww = wa; | |
| 458 | + } | |
| 459 | + bk0 = M_PI_2*sum/sin(piv); | |
| 460 | + } | |
| 461 | + } | |
| 462 | + else { | |
| 463 | + cb = exp(-x)*sqrt(M_PI_2/x); | |
| 464 | + sum = 1.0; | |
| 465 | + r = 1.0; | |
| 466 | + for (k=1;k<=kz;k++) { | |
| 467 | + r *= 0.125*(vt-pow(2.0*k-1.0,2.0))/(k*x); | |
| 468 | + sum += r; | |
| 469 | + } | |
| 470 | + bk0 = cb*sum; | |
| 471 | + } | |
| 472 | + bk1 = (1.0/x-iv[1]*bk0)/iv[0]; | |
| 473 | + kv[0] = bk0; | |
| 474 | + kv[1] = bk1; | |
| 475 | + for (k=2;k<=n;k++) { | |
| 476 | + bk2 = 2.0*(v0+k-1.0)*bk1/x+bk0; | |
| 477 | + kv[k] = bk2; | |
| 478 | + bk0 = bk1; | |
| 479 | + bk1 = bk2; | |
| 480 | + } | |
| 481 | + kvp[0] = v0*kv[0]/x-kv[1]; | |
| 482 | + for (k=1;k<=n;k++) { | |
| 483 | + kvp[k] = -(k+v0)*kv[k]/x-kv[k-1]; | |
| 484 | + } | |
| 485 | + vm = n+v0; | |
| 486 | + return 0; | |
| 487 | +} | ... | ... |
| 1 | +++ a/bessjy.cpp | |
| 1 | +// bessjy.cpp -- computation of Bessel functions Jn, Yn and their | |
| 2 | +// derivatives. Algorithms and coefficient values from | |
| 3 | +// "Computation of Special Functions", Zhang and Jin, John | |
| 4 | +// Wiley and Sons, 1996. | |
| 5 | +// | |
| 6 | +// (C) 2003, C. Bond. All rights reserved. | |
| 7 | +// | |
| 8 | +// Note that 'math.h' provides (or should provide) values for: | |
| 9 | +// pi M_PI | |
| 10 | +// 2/pi M_2_PI | |
| 11 | +// pi/4 M_PI_4 | |
| 12 | +// pi/2 M_PI_2 | |
| 13 | +// | |
| 14 | +#define _USE_MATH_DEFINES | |
| 15 | +#include <math.h> | |
| 16 | +#include "bessel.h" | |
| 17 | + | |
| 18 | +double gamma(double x); | |
| 19 | +// | |
| 20 | +// INPUT: | |
| 21 | +// double x -- argument of Bessel function | |
| 22 | +// | |
| 23 | +// OUTPUT (via address rtsPointers): | |
| 24 | +// double j0 -- Bessel function of 1st kind, 0th order | |
| 25 | +// double j1 -- Bessel function of 1st kind, 1st order | |
| 26 | +// double y0 -- Bessel function of 2nd kind, 0th order | |
| 27 | +// double y1 -- Bessel function of 2nd kind, 1st order | |
| 28 | +// double j0p -- derivative of Bessel function of 1st kind, 0th order | |
| 29 | +// double j1p -- derivative of Bessel function of 1st kind, 1st order | |
| 30 | +// double y0p -- derivative of Bessel function of 2nd kind, 0th order | |
| 31 | +// double y1p -- derivative of Bessel function of 2nd kind, 1st order | |
| 32 | +// | |
| 33 | +// RETURN: | |
| 34 | +// int error code: 0 = OK, 1 = error | |
| 35 | +// | |
| 36 | +// This algorithm computes the above functions using series expansions. | |
| 37 | +// | |
| 38 | +int bessjy01a(double x,double &j0,double &j1,double &y0,double &y1, | |
| 39 | + double &j0p,double &j1p,double &y0p,double &y1p) | |
| 40 | +{ | |
| 41 | + double x2,r,ec,w0,w1,r0,r1,cs0,cs1; | |
| 42 | + double cu,p0,q0,p1,q1,t1,t2; | |
| 43 | + int k,kz; | |
| 44 | + static double a[] = { | |
| 45 | + -7.03125e-2, | |
| 46 | + 0.112152099609375, | |
| 47 | + -0.5725014209747314, | |
| 48 | + 6.074042001273483, | |
| 49 | + -1.100171402692467e2, | |
| 50 | + 3.038090510922384e3, | |
| 51 | + -1.188384262567832e5, | |
| 52 | + 6.252951493434797e6, | |
| 53 | + -4.259392165047669e8, | |
| 54 | + 3.646840080706556e10, | |
| 55 | + -3.833534661393944e12, | |
| 56 | + 4.854014686852901e14, | |
| 57 | + -7.286857349377656e16, | |
| 58 | + 1.279721941975975e19}; | |
| 59 | + static double b[] = { | |
| 60 | + 7.32421875e-2, | |
| 61 | + -0.2271080017089844, | |
| 62 | + 1.727727502584457, | |
| 63 | + -2.438052969955606e1, | |
| 64 | + 5.513358961220206e2, | |
| 65 | + -1.825775547429318e4, | |
| 66 | + 8.328593040162893e5, | |
| 67 | + -5.006958953198893e7, | |
| 68 | + 3.836255180230433e9, | |
| 69 | + -3.649010818849833e11, | |
| 70 | + 4.218971570284096e13, | |
| 71 | + -5.827244631566907e15, | |
| 72 | + 9.476288099260110e17, | |
| 73 | + -1.792162323051699e20}; | |
| 74 | + static double a1[] = { | |
| 75 | + 0.1171875, | |
| 76 | + -0.1441955566406250, | |
| 77 | + 0.6765925884246826, | |
| 78 | + -6.883914268109947, | |
| 79 | + 1.215978918765359e2, | |
| 80 | + -3.302272294480852e3, | |
| 81 | + 1.276412726461746e5, | |
| 82 | + -6.656367718817688e6, | |
| 83 | + 4.502786003050393e8, | |
| 84 | + -3.833857520742790e10, | |
| 85 | + 4.011838599133198e12, | |
| 86 | + -5.060568503314727e14, | |
| 87 | + 7.572616461117958e16, | |
| 88 | + -1.326257285320556e19}; | |
| 89 | + static double b1[] = { | |
| 90 | + -0.1025390625, | |
| 91 | + 0.2775764465332031, | |
| 92 | + -1.993531733751297, | |
| 93 | + 2.724882731126854e1, | |
| 94 | + -6.038440767050702e2, | |
| 95 | + 1.971837591223663e4, | |
| 96 | + -8.902978767070678e5, | |
| 97 | + 5.310411010968522e7, | |
| 98 | + -4.043620325107754e9, | |
| 99 | + 3.827011346598605e11, | |
| 100 | + -4.406481417852278e13, | |
| 101 | + 6.065091351222699e15, | |
| 102 | + -9.833883876590679e17, | |
| 103 | + 1.855045211579828e20}; | |
| 104 | + | |
| 105 | + if (x < 0.0) return 1; | |
| 106 | + if (x == 0.0) { | |
| 107 | + j0 = 1.0; | |
| 108 | + j1 = 0.0; | |
| 109 | + y0 = -1e308; | |
| 110 | + y1 = -1e308; | |
| 111 | + j0p = 0.0; | |
| 112 | + j1p = 0.5; | |
| 113 | + y0p = 1e308; | |
| 114 | + y1p = 1e308; | |
| 115 | + return 0; | |
| 116 | + } | |
| 117 | + x2 = x*x; | |
| 118 | + if (x <= 12.0) { | |
| 119 | + j0 = 1.0; | |
| 120 | + r = 1.0; | |
| 121 | + for (k=1;k<=30;k++) { | |
| 122 | + r *= -0.25*x2/(k*k); | |
| 123 | + j0 += r; | |
| 124 | + if (fabs(r) < fabs(j0)*1e-15) break; | |
| 125 | + } | |
| 126 | + j1 = 1.0; | |
| 127 | + r = 1.0; | |
| 128 | + for (k=1;k<=30;k++) { | |
| 129 | + r *= -0.25*x2/(k*(k+1)); | |
| 130 | + j1 += r; | |
| 131 | + if (fabs(r) < fabs(j1)*1e-15) break; | |
| 132 | + } | |
| 133 | + j1 *= 0.5*x; | |
| 134 | + ec = log(0.5*x)+el; | |
| 135 | + cs0 = 0.0; | |
| 136 | + w0 = 0.0; | |
| 137 | + r0 = 1.0; | |
| 138 | + for (k=1;k<=30;k++) { | |
| 139 | + w0 += 1.0/k; | |
| 140 | + r0 *= -0.25*x2/(k*k); | |
| 141 | + r = r0 * w0; | |
| 142 | + cs0 += r; | |
| 143 | + if (fabs(r) < fabs(cs0)*1e-15) break; | |
| 144 | + } | |
| 145 | + y0 = M_2_PI*(ec*j0-cs0); | |
| 146 | + cs1 = 1.0; | |
| 147 | + w1 = 0.0; | |
| 148 | + r1 = 1.0; | |
| 149 | + for (k=1;k<=30;k++) { | |
| 150 | + w1 += 1.0/k; | |
| 151 | + r1 *= -0.25*x2/(k*(k+1)); | |
| 152 | + r = r1*(2.0*w1+1.0/(k+1)); | |
| 153 | + cs1 += r; | |
| 154 | + if (fabs(r) < fabs(cs1)*1e-15) break; | |
| 155 | + } | |
| 156 | + y1 = M_2_PI * (ec*j1-1.0/x-0.25*x*cs1); | |
| 157 | + } | |
| 158 | + else { | |
| 159 | + if (x >= 50.0) kz = 8; // Can be changed to 10 | |
| 160 | + else if (x >= 35.0) kz = 10; // " " 12 | |
| 161 | + else kz = 12; // " " 14 | |
| 162 | + t1 = x-M_PI_4; | |
| 163 | + p0 = 1.0; | |
| 164 | + q0 = -0.125/x; | |
| 165 | + for (k=0;k<kz;k++) { | |
| 166 | + p0 += a[k]*pow(x,-2*k-2); | |
| 167 | + q0 += b[k]*pow(x,-2*k-3); | |
| 168 | + } | |
| 169 | + cu = sqrt(M_2_PI/x); | |
| 170 | + j0 = cu*(p0*cos(t1)-q0*sin(t1)); | |
| 171 | + y0 = cu*(p0*sin(t1)+q0*cos(t1)); | |
| 172 | + t2 = x-0.75*M_PI; | |
| 173 | + p1 = 1.0; | |
| 174 | + q1 = 0.375/x; | |
| 175 | + for (k=0;k<kz;k++) { | |
| 176 | + p1 += a1[k]*pow(x,-2*k-2); | |
| 177 | + q1 += b1[k]*pow(x,-2*k-3); | |
| 178 | + } | |
| 179 | + j1 = cu*(p1*cos(t2)-q1*sin(t2)); | |
| 180 | + y1 = cu*(p1*sin(t2)+q1*cos(t2)); | |
| 181 | + } | |
| 182 | + j0p = -j1; | |
| 183 | + j1p = j0-j1/x; | |
| 184 | + y0p = -y1; | |
| 185 | + y1p = y0-y1/x; | |
| 186 | + return 0; | |
| 187 | +} | |
| 188 | +// | |
| 189 | +// INPUT: | |
| 190 | +// double x -- argument of Bessel function | |
| 191 | +// | |
| 192 | +// OUTPUT: | |
| 193 | +// double j0 -- Bessel function of 1st kind, 0th order | |
| 194 | +// double j1 -- Bessel function of 1st kind, 1st order | |
| 195 | +// double y0 -- Bessel function of 2nd kind, 0th order | |
| 196 | +// double y1 -- Bessel function of 2nd kind, 1st order | |
| 197 | +// double j0p -- derivative of Bessel function of 1st kind, 0th order | |
| 198 | +// double j1p -- derivative of Bessel function of 1st kind, 1st order | |
| 199 | +// double y0p -- derivative of Bessel function of 2nd kind, 0th order | |
| 200 | +// double y1p -- derivative of Bessel function of 2nd kind, 1st order | |
| 201 | +// | |
| 202 | +// RETURN: | |
| 203 | +// int error code: 0 = OK, 1 = error | |
| 204 | +// | |
| 205 | +// This algorithm computes the functions using polynomial approximations. | |
| 206 | +// | |
| 207 | +int bessjy01b(double x,double &j0,double &j1,double &y0,double &y1, | |
| 208 | + double &j0p,double &j1p,double &y0p,double &y1p) | |
| 209 | +{ | |
| 210 | + double t,t2,dtmp,a0,p0,q0,p1,q1,ta0,ta1; | |
| 211 | + if (x < 0.0) return 1; | |
| 212 | + if (x == 0.0) { | |
| 213 | + j0 = 1.0; | |
| 214 | + j1 = 0.0; | |
| 215 | + y0 = -1e308; | |
| 216 | + y1 = -1e308; | |
| 217 | + j0p = 0.0; | |
| 218 | + j1p = 0.5; | |
| 219 | + y0p = 1e308; | |
| 220 | + y1p = 1e308; | |
| 221 | + return 0; | |
| 222 | + } | |
| 223 | + if(x <= 4.0) { | |
| 224 | + t = x/4.0; | |
| 225 | + t2 = t*t; | |
| 226 | + j0 = ((((((-0.5014415e-3*t2+0.76771853e-2)*t2-0.0709253492)*t2+ | |
| 227 | + 0.4443584263)*t2-1.7777560599)*t2+3.9999973021)*t2 | |
| 228 | + -3.9999998721)*t2+1.0; | |
| 229 | + j1 = t*(((((((-0.1289769e-3*t2+0.22069155e-2)*t2-0.0236616773)*t2+ | |
| 230 | + 0.1777582922)*t2-0.8888839649)*t2+2.6666660544)*t2- | |
| 231 | + 3.999999971)*t2+1.9999999998); | |
| 232 | + dtmp = (((((((-0.567433e-4*t2+0.859977e-3)*t2-0.94855882e-2)*t2+ | |
| 233 | + 0.0772975809)*t2-0.4261737419)*t2+1.4216421221)*t2- | |
| 234 | + 2.3498519931)*t2+1.0766115157)*t2+0.3674669052; | |
| 235 | + y0 = M_2_PI*log(0.5*x)*j0+dtmp; | |
| 236 | + dtmp = (((((((0.6535773e-3*t2-0.0108175626)*t2+0.107657607)*t2- | |
| 237 | + 0.7268945577)*t2+3.1261399273)*t2-7.3980241381)*t2+ | |
| 238 | + 6.8529236342)*t2+0.3932562018)*t2-0.6366197726; | |
| 239 | + y1 = M_2_PI*log(0.5*x)*j1+dtmp/x; | |
| 240 | + } | |
| 241 | + else { | |
| 242 | + t = 4.0/x; | |
| 243 | + t2 = t*t; | |
| 244 | + a0 = sqrt(M_2_PI/x); | |
| 245 | + p0 = ((((-0.9285e-5*t2+0.43506e-4)*t2-0.122226e-3)*t2+ | |
| 246 | + 0.434725e-3)*t2-0.4394275e-2)*t2+0.999999997; | |
| 247 | + q0 = t*(((((0.8099e-5*t2-0.35614e-4)*t2+0.85844e-4)*t2- | |
| 248 | + 0.218024e-3)*t2+0.1144106e-2)*t2-0.031249995); | |
| 249 | + ta0 = x-M_PI_4; | |
| 250 | + j0 = a0*(p0*cos(ta0)-q0*sin(ta0)); | |
| 251 | + y0 = a0*(p0*sin(ta0)+q0*cos(ta0)); | |
| 252 | + p1 = ((((0.10632e-4*t2-0.50363e-4)*t2+0.145575e-3)*t2 | |
| 253 | + -0.559487e-3)*t2+0.7323931e-2)*t2+1.000000004; | |
| 254 | + q1 = t*(((((-0.9173e-5*t2+0.40658e-4)*t2-0.99941e-4)*t2 | |
| 255 | + +0.266891e-3)*t2-0.1601836e-2)*t2+0.093749994); | |
| 256 | + ta1 = x-0.75*M_PI; | |
| 257 | + j1 = a0*(p1*cos(ta1)-q1*sin(ta1)); | |
| 258 | + y1 = a0*(p1*sin(ta1)+q1*cos(ta1)); | |
| 259 | + } | |
| 260 | + j0p = -j1; | |
| 261 | + j1p = j0-j1/x; | |
| 262 | + y0p = -y1; | |
| 263 | + y1p = y0-y1/x; | |
| 264 | + return 0; | |
| 265 | +} | |
| 266 | +int msta1(double x,int mp) | |
| 267 | +{ | |
| 268 | + double a0,f0,f1,f; | |
| 269 | + int i,n0,n1,nn; | |
| 270 | + | |
| 271 | + a0 = fabs(x); | |
| 272 | + n0 = (int)(1.1*a0)+1; | |
| 273 | + f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-mp; | |
| 274 | + n1 = n0+5; | |
| 275 | + f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-mp; | |
| 276 | + for (i=0;i<20;i++) { | |
| 277 | + nn = (int)(n1-(n1-n0)/(1.0-f0/f1)); | |
| 278 | + f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-mp; | |
| 279 | + if (abs(nn-n1) < 1) break; | |
| 280 | + n0 = n1; | |
| 281 | + f0 = f1; | |
| 282 | + n1 = nn; | |
| 283 | + f1 = f; | |
| 284 | + } | |
| 285 | + return nn; | |
| 286 | +} | |
| 287 | +int msta2(double x,int n,int mp) | |
| 288 | +{ | |
| 289 | + double a0,ejn,hmp,f0,f1,f,obj; | |
| 290 | + int i,n0,n1,nn; | |
| 291 | + | |
| 292 | + a0 = fabs(x); | |
| 293 | + hmp = 0.5*mp; | |
| 294 | + ejn = 0.5*log10(6.28*n)-n*log10(1.36*a0/n); | |
| 295 | + if (ejn <= hmp) { | |
| 296 | + obj = mp; | |
| 297 | + n0 = (int)(1.1*a0); | |
| 298 | + if (n0 < 1) n0 = 1; | |
| 299 | + } | |
| 300 | + else { | |
| 301 | + obj = hmp+ejn; | |
| 302 | + n0 = n; | |
| 303 | + } | |
| 304 | + f0 = 0.5*log10(6.28*n0)-n0*log10(1.36*a0/n0)-obj; | |
| 305 | + n1 = n0+5; | |
| 306 | + f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-obj; | |
| 307 | + for (i=0;i<20;i++) { | |
| 308 | + nn = (int)(n1-(n1-n0)/(1.0-f0/f1)); | |
| 309 | + f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-obj; | |
| 310 | + if (abs(nn-n1) < 1) break; | |
| 311 | + n0 = n1; | |
| 312 | + f0 = f1; | |
| 313 | + n1 = nn; | |
| 314 | + f1 = f; | |
| 315 | + } | |
| 316 | + return nn+10; | |
| 317 | +} | |
| 318 | +// | |
| 319 | +// INPUT: | |
| 320 | +// double x -- argument of Bessel function of 1st and 2nd kind. | |
| 321 | +// int n -- order | |
| 322 | +// | |
| 323 | +// OUPUT: | |
| 324 | +// | |
| 325 | +// int nm -- highest order actually computed (nm <= n) | |
| 326 | +// double jn[] -- Bessel function of 1st kind, orders from 0 to nm | |
| 327 | +// double yn[] -- Bessel function of 2nd kind, orders from 0 to nm | |
| 328 | +// double j'n[]-- derivative of Bessel function of 1st kind, | |
| 329 | +// orders from 0 to nm | |
| 330 | +// double y'n[]-- derivative of Bessel function of 2nd kind, | |
| 331 | +// orders from 0 to nm | |
| 332 | +// | |
| 333 | +// Computes Bessel functions of all order up to 'n' using recurrence | |
| 334 | +// relations. If 'nm' < 'n' only 'nm' orders are returned. | |
| 335 | +// | |
| 336 | +int bessjyna(int n,double x,int &nm,double *jn,double *yn, | |
| 337 | + double *jnp,double *ynp) | |
| 338 | +{ | |
| 339 | + double bj0,bj1,f,f0,f1,f2,cs; | |
| 340 | + int i,k,m,ecode; | |
| 341 | + | |
| 342 | + nm = n; | |
| 343 | + if ((x < 0.0) || (n < 0)) return 1; | |
| 344 | + if (x < 1e-15) { | |
| 345 | + for (i=0;i<=n;i++) { | |
| 346 | + jn[i] = 0.0; | |
| 347 | + yn[i] = -1e308; | |
| 348 | + jnp[i] = 0.0; | |
| 349 | + ynp[i] = 1e308; | |
| 350 | + } | |
| 351 | + jn[0] = 1.0; | |
| 352 | + jnp[1] = 0.5; | |
| 353 | + return 0; | |
| 354 | + } | |
| 355 | + ecode = bessjy01a(x,jn[0],jn[1],yn[0],yn[1],jnp[0],jnp[1],ynp[0],ynp[1]); | |
| 356 | + if (n < 2) return 0; | |
| 357 | + bj0 = jn[0]; | |
| 358 | + bj1 = jn[1]; | |
| 359 | + if (n < (int)0.9*x) { | |
| 360 | + for (k=2;k<=n;k++) { | |
| 361 | + jn[k] = 2.0*(k-1.0)*bj1/x-bj0; | |
| 362 | + bj0 = bj1; | |
| 363 | + bj1 = jn[k]; | |
| 364 | + } | |
| 365 | + } | |
| 366 | + else { | |
| 367 | + m = msta1(x,200); | |
| 368 | + if (m < n) nm = m; | |
| 369 | + else m = msta2(x,n,15); | |
| 370 | + f2 = 0.0; | |
| 371 | + f1 = 1.0e-100; | |
| 372 | + for (k=m;k>=0;k--) { | |
| 373 | + f = 2.0*(k+1.0)/x*f1-f2; | |
| 374 | + if (k <= nm) jn[k] = f; | |
| 375 | + f2 = f1; | |
| 376 | + f1 = f; | |
| 377 | + } | |
| 378 | + if (fabs(bj0) > fabs(bj1)) cs = bj0/f; | |
| 379 | + else cs = bj1/f2; | |
| 380 | + for (k=0;k<=nm;k++) { | |
| 381 | + jn[k] *= cs; | |
| 382 | + } | |
| 383 | + } | |
| 384 | + for (k=2;k<=nm;k++) { | |
| 385 | + jnp[k] = jn[k-1]-k*jn[k]/x; | |
| 386 | + } | |
| 387 | + f0 = yn[0]; | |
| 388 | + f1 = yn[1]; | |
| 389 | + for (k=2;k<=nm;k++) { | |
| 390 | + f = 2.0*(k-1.0)*f1/x-f0; | |
| 391 | + yn[k] = f; | |
| 392 | + f0 = f1; | |
| 393 | + f1 = f; | |
| 394 | + } | |
| 395 | + for (k=2;k<=nm;k++) { | |
| 396 | + ynp[k] = yn[k-1]-k*yn[k]/x; | |
| 397 | + } | |
| 398 | + return 0; | |
| 399 | +} | |
| 400 | +// | |
| 401 | +// Same input and output conventions as above. Different recurrence | |
| 402 | +// relations used for 'x' < 300. | |
| 403 | +// | |
| 404 | +int bessjynb(int n,double x,int &nm,double *jn,double *yn, | |
| 405 | + double *jnp,double *ynp) | |
| 406 | +{ | |
| 407 | + double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv; | |
| 408 | + double ec,bs,byk,p0,p1,q0,q1; | |
| 409 | + static double a[] = { | |
| 410 | + -0.7031250000000000e-1, | |
| 411 | + 0.1121520996093750, | |
| 412 | + -0.5725014209747314, | |
| 413 | + 6.074042001273483}; | |
| 414 | + static double b[] = { | |
| 415 | + 0.7324218750000000e-1, | |
| 416 | + -0.2271080017089844, | |
| 417 | + 1.727727502584457, | |
| 418 | + -2.438052969955606e1}; | |
| 419 | + static double a1[] = { | |
| 420 | + 0.1171875, | |
| 421 | + -0.1441955566406250, | |
| 422 | + 0.6765925884246826, | |
| 423 | + -6.883914268109947}; | |
| 424 | + static double b1[] = { | |
| 425 | + -0.1025390625, | |
| 426 | + 0.2775764465332031, | |
| 427 | + -1.993531733751297, | |
| 428 | + 2.724882731126854e1}; | |
| 429 | + | |
| 430 | + int i,k,m; | |
| 431 | + nm = n; | |
| 432 | + if ((x < 0.0) || (n < 0)) return 1; | |
| 433 | + if (x < 1e-15) { | |
| 434 | + for (i=0;i<=n;i++) { | |
| 435 | + jn[i] = 0.0; | |
| 436 | + yn[i] = -1e308; | |
| 437 | + jnp[i] = 0.0; | |
| 438 | + ynp[i] = 1e308; | |
| 439 | + } | |
| 440 | + jn[0] = 1.0; | |
| 441 | + jnp[1] = 0.5; | |
| 442 | + return 0; | |
| 443 | + } | |
| 444 | + if (x <= 300.0 || n > (int)(0.9*x)) { | |
| 445 | + if (n == 0) nm = 1; | |
| 446 | + m = msta1(x,200); | |
| 447 | + if (m < nm) nm = m; | |
| 448 | + else m = msta2(x,nm,15); | |
| 449 | + bs = 0.0; | |
| 450 | + su = 0.0; | |
| 451 | + sv = 0.0; | |
| 452 | + f2 = 0.0; | |
| 453 | + f1 = 1.0e-100; | |
| 454 | + for (k = m;k>=0;k--) { | |
| 455 | + f = 2.0*(k+1.0)/x*f1 - f2; | |
| 456 | + if (k <= nm) jn[k] = f; | |
| 457 | + if ((k == 2*(int)(k/2)) && (k != 0)) { | |
| 458 | + bs += 2.0*f; | |
| 459 | +// su += pow(-1,k>>1)*f/(double)k; | |
| 460 | + su += (-1)*((k & 2)-1)*f/(double)k; | |
| 461 | + } | |
| 462 | + else if (k > 1) { | |
| 463 | +// sv += pow(-1,k>>1)*k*f/(k*k-1.0); | |
| 464 | + sv += (-1)*((k & 2)-1)*(double)k*f/(k*k-1.0); | |
| 465 | + } | |
| 466 | + f2 = f1; | |
| 467 | + f1 = f; | |
| 468 | + } | |
| 469 | + s0 = bs+f; | |
| 470 | + for (k=0;k<=nm;k++) { | |
| 471 | + jn[k] /= s0; | |
| 472 | + } | |
| 473 | + ec = log(0.5*x) +0.5772156649015329; | |
| 474 | + by0 = M_2_PI*(ec*jn[0]-4.0*su/s0); | |
| 475 | + yn[0] = by0; | |
| 476 | + by1 = M_2_PI*((ec-1.0)*jn[1]-jn[0]/x-4.0*sv/s0); | |
| 477 | + yn[1] = by1; | |
| 478 | + } | |
| 479 | + else { | |
| 480 | + t1 = x-M_PI_4; | |
| 481 | + p0 = 1.0; | |
| 482 | + q0 = -0.125/x; | |
| 483 | + for (k=0;k<4;k++) { | |
| 484 | + p0 += a[k]*pow(x,-2*k-2); | |
| 485 | + q0 += b[k]*pow(x,-2*k-3); | |
| 486 | + } | |
| 487 | + cu = sqrt(M_2_PI/x); | |
| 488 | + bj0 = cu*(p0*cos(t1)-q0*sin(t1)); | |
| 489 | + by0 = cu*(p0*sin(t1)+q0*cos(t1)); | |
| 490 | + jn[0] = bj0; | |
| 491 | + yn[0] = by0; | |
| 492 | + t2 = x-0.75*M_PI; | |
| 493 | + p1 = 1.0; | |
| 494 | + q1 = 0.375/x; | |
| 495 | + for (k=0;k<4;k++) { | |
| 496 | + p1 += a1[k]*pow(x,-2*k-2); | |
| 497 | + q1 += b1[k]*pow(x,-2*k-3); | |
| 498 | + } | |
| 499 | + bj1 = cu*(p1*cos(t2)-q1*sin(t2)); | |
| 500 | + by1 = cu*(p1*sin(t2)+q1*cos(t2)); | |
| 501 | + jn[1] = bj1; | |
| 502 | + yn[1] = by1; | |
| 503 | + for (k=2;k<=nm;k++) { | |
| 504 | + bjk = 2.0*(k-1.0)*bj1/x-bj0; | |
| 505 | + jn[k] = bjk; | |
| 506 | + bj0 = bj1; | |
| 507 | + bj1 = bjk; | |
| 508 | + } | |
| 509 | + } | |
| 510 | + jnp[0] = -jn[1]; | |
| 511 | + for (k=1;k<=nm;k++) { | |
| 512 | + jnp[k] = jn[k-1]-k*jn[k]/x; | |
| 513 | + } | |
| 514 | + for (k=2;k<=nm;k++) { | |
| 515 | + byk = 2.0*(k-1.0)*by1/x-by0; | |
| 516 | + yn[k] = byk; | |
| 517 | + by0 = by1; | |
| 518 | + by1 = byk; | |
| 519 | + } | |
| 520 | + ynp[0] = -yn[1]; | |
| 521 | + for (k=1;k<=nm;k++) { | |
| 522 | + ynp[k] = yn[k-1]-k*yn[k]/x; | |
| 523 | + } | |
| 524 | + return 0; | |
| 525 | + | |
| 526 | +} | |
| 527 | + | |
| 528 | +// The following routine computes Bessel Jv(x) and Yv(x) for | |
| 529 | +// arbitrary positive order (v). For negative order, use: | |
| 530 | +// | |
| 531 | +// J-v(x) = Jv(x)cos(v pi) - Yv(x)sin(v pi) | |
| 532 | +// Y-v(x) = Jv(x)sin(v pi) + Yv(x)cos(v pi) | |
| 533 | +// | |
| 534 | +int bessjyv(double v,double x,double &vm,double *jv,double *yv, | |
| 535 | + double *djv,double *dyv) | |
| 536 | +{ | |
| 537 | + double v0,vl,vg,vv,a,a0,r,x2,bjv0,bjv1,bjvl,f,f0,f1,f2; | |
| 538 | + double r0,r1,ck,cs,cs0,cs1,sk,qx,px,byv0,byv1,rp,xk,rq; | |
| 539 | + double b,ec,w0,w1,bju0,bju1,pv0,pv1,byvk; | |
| 540 | + int j,k,l,m,n,kz; | |
| 541 | + | |
| 542 | + x2 = x*x; | |
| 543 | + n = (int)v; | |
| 544 | + v0 = v-n; | |
| 545 | + if ((x < 0.0) || (v < 0.0)) return 1; | |
| 546 | + if (x < 1e-15) { | |
| 547 | + for (k=0;k<=n;k++) { | |
| 548 | + jv[k] = 0.0; | |
| 549 | + yv[k] = -1e308; | |
| 550 | + djv[k] = 0.0; | |
| 551 | + dyv[k] = 1e308; | |
| 552 | + if (v0 == 0.0) { | |
| 553 | + jv[0] = 1.0; | |
| 554 | + djv[1] = 0.5; | |
| 555 | + } | |
| 556 | + else djv[0] = 1e308; | |
| 557 | + } | |
| 558 | + vm = v; | |
| 559 | + return 0; | |
| 560 | + } | |
| 561 | + if (x <= 12.0) { | |
| 562 | + for (l=0;l<2;l++) { | |
| 563 | + vl = v0 + l; | |
| 564 | + bjvl = 1.0; | |
| 565 | + r = 1.0; | |
| 566 | + for (k=1;k<=40;k++) { | |
| 567 | + r *= -0.25*x2/(k*(k+vl)); | |
| 568 | + bjvl += r; | |
| 569 | + if (fabs(r) < fabs(bjvl)*1e-15) break; | |
| 570 | + } | |
| 571 | + vg = 1.0 + vl; | |
| 572 | + a = pow(0.5*x,vl)/gamma(vg); | |
| 573 | + if (l == 0) bjv0 = bjvl*a; | |
| 574 | + else bjv1 = bjvl*a; | |
| 575 | + } | |
| 576 | + } | |
| 577 | + else { | |
| 578 | + if (x >= 50.0) kz = 8; | |
| 579 | + else if (x >= 35.0) kz = 10; | |
| 580 | + else kz = 11; | |
| 581 | + for (j=0;j<2;j++) { | |
| 582 | + vv = 4.0*(j+v0)*(j+v0); | |
| 583 | + px = 1.0; | |
| 584 | + rp = 1.0; | |
| 585 | + for (k=1;k<=kz;k++) { | |
| 586 | + rp *= (-0.78125e-2)*(vv-pow(4.0*k-3.0,2.0))* | |
| 587 | + (vv-pow(4.0*k-1.0,2.0))/(k*(2.0*k-1.0)*x2); | |
| 588 | + px += rp; | |
| 589 | + } | |
| 590 | + qx = 1.0; | |
| 591 | + rq = 1.0; | |
| 592 | + for (k=1;k<=kz;k++) { | |
| 593 | + rq *= (-0.78125e-2)*(vv-pow(4.0*k-1.0,2.0))* | |
| 594 | + (vv-pow(4.0*k+1.0,2.0))/(k*(2.0*k+1.0)*x2); | |
| 595 | + qx += rq; | |
| 596 | + } | |
| 597 | + qx *= 0.125*(vv-1.0)/x; | |
| 598 | + xk = x-(0.5*(j+v0)+0.25)*M_PI; | |
| 599 | + a0 = sqrt(M_2_PI/x); | |
| 600 | + ck = cos(xk); | |
| 601 | + sk = sin(xk); | |
| 602 | + | |
| 603 | + if (j == 0) { | |
| 604 | + bjv0 = a0*(px*ck-qx*sk); | |
| 605 | + byv0 = a0*(px*sk+qx*ck); | |
| 606 | + } | |
| 607 | + else if (j == 1) { | |
| 608 | + bjv1 = a0*(px*ck-qx*sk); | |
| 609 | + byv1 = a0*(px*sk+qx*ck); | |
| 610 | + } | |
| 611 | + } | |
| 612 | + } | |
| 613 | + jv[0] = bjv0; | |
| 614 | + jv[1] = bjv1; | |
| 615 | + djv[0] = v0*jv[0]/x-jv[1]; | |
| 616 | + djv[1] = -(1.0+v0)*jv[1]/x+jv[0]; | |
| 617 | + if ((n >= 2) && (n <= (int)(0.9*x))) { | |
| 618 | + f0 = bjv0; | |
| 619 | + f1 = bjv1; | |
| 620 | + for (k=2;k<=n;k++) { | |
| 621 | + f = 2.0*(k+v0-1.0)*f1/x-f0; | |
| 622 | + jv[k] = f; | |
| 623 | + f0 = f1; | |
| 624 | + f1 = f; | |
| 625 | + } | |
| 626 | + } | |
| 627 | + else if (n >= 2) { | |
| 628 | + m = msta1(x,200); | |
| 629 | + if (m < n) n = m; | |
| 630 | + else m = msta2(x,n,15); | |
| 631 | + f2 = 0.0; | |
| 632 | + f1 = 1.0e-100; | |
| 633 | + for (k=m;k>=0;k--) { | |
| 634 | + f = 2.0*(v0+k+1.0)*f1/x-f2; | |
| 635 | + if (k <= n) jv[k] = f; | |
| 636 | + f2 = f1; | |
| 637 | + f1 = f; | |
| 638 | + } | |
| 639 | + if (fabs(bjv0) > fabs(bjv1)) cs = bjv0/f; | |
| 640 | + else cs = bjv1/f2; | |
| 641 | + for (k=0;k<=n;k++) { | |
| 642 | + jv[k] *= cs; | |
| 643 | + } | |
| 644 | + } | |
| 645 | + for (k=2;k<=n;k++) { | |
| 646 | + djv[k] = -(k+v0)*jv[k]/x+jv[k-1]; | |
| 647 | + } | |
| 648 | + if (x <= 12.0) { | |
| 649 | + if (v0 != 0.0) { | |
| 650 | + for (l=0;l<2;l++) { | |
| 651 | + vl = v0 +l; | |
| 652 | + bjvl = 1.0; | |
| 653 | + r = 1.0; | |
| 654 | + for (k=1;k<=40;k++) { | |
| 655 | + r *= -0.25*x2/(k*(k-vl)); | |
| 656 | + bjvl += r; | |
| 657 | + if (fabs(r) < fabs(bjvl)*1e-15) break; | |
| 658 | + } | |
| 659 | + vg = 1.0-vl; | |
| 660 | + b = pow(2.0/x,vl)/gamma(vg); | |
| 661 | + if (l == 0) bju0 = bjvl*b; | |
| 662 | + else bju1 = bjvl*b; | |
| 663 | + } | |
| 664 | + pv0 = M_PI*v0; | |
| 665 | + pv1 = M_PI*(1.0+v0); | |
| 666 | + byv0 = (bjv0*cos(pv0)-bju0)/sin(pv0); | |
| 667 | + byv1 = (bjv1*cos(pv1)-bju1)/sin(pv1); | |
| 668 | + } | |
| 669 | + else { | |
| 670 | + ec = log(0.5*x)+el; | |
| 671 | + cs0 = 0.0; | |
| 672 | + w0 = 0.0; | |
| 673 | + r0 = 1.0; | |
| 674 | + for (k=1;k<=30;k++) { | |
| 675 | + w0 += 1.0/k; | |
| 676 | + r0 *= -0.25*x2/(k*k); | |
| 677 | + cs0 += r0*w0; | |
| 678 | + } | |
| 679 | + byv0 = M_2_PI*(ec*bjv0-cs0); | |
| 680 | + cs1 = 1.0; | |
| 681 | + w1 = 0.0; | |
| 682 | + r1 = 1.0; | |
| 683 | + for (k=1;k<=30;k++) { | |
| 684 | + w1 += 1.0/k; | |
| 685 | + r1 *= -0.25*x2/(k*(k+1)); | |
| 686 | + cs1 += r1*(2.0*w1+1.0/(k+1.0)); | |
| 687 | + } | |
| 688 | + byv1 = M_2_PI*(ec*bjv1-1.0/x-0.25*x*cs1); | |
| 689 | + } | |
| 690 | + } | |
| 691 | + yv[0] = byv0; | |
| 692 | + yv[1] = byv1; | |
| 693 | + for (k=2;k<=n;k++) { | |
| 694 | + byvk = 2.0*(v0+k-1.0)*byv1/x-byv0; | |
| 695 | + yv[k] = byvk; | |
| 696 | + byv0 = byv1; | |
| 697 | + byv1 = byvk; | |
| 698 | + } | |
| 699 | + dyv[0] = v0*yv[0]/x-yv[1]; | |
| 700 | + for (k=1;k<=n;k++) { | |
| 701 | + dyv[k] = -(k+v0)*yv[k]/x+yv[k-1]; | |
| 702 | + } | |
| 703 | + vm = n + v0; | |
| 704 | + return 0; | |
| 705 | +} | |
| 706 | + | ... | ... |
| 1 | +++ a/cbessik.cpp | |
| 1 | +// cbessik.cpp -- complex modified Bessel functions. | |
| 2 | +// Algorithms and coefficient values from "Computation of Special | |
| 3 | +// Functions", Zhang and Jin, John Wiley and Sons, 1996. | |
| 4 | +// | |
| 5 | +// (C) 2003, C. Bond. All rights reserved. | |
| 6 | +// | |
| 7 | +#include <complex> | |
| 8 | +using namespace std; | |
| 9 | +#include "bessel.h" | |
| 10 | + | |
| 11 | +static complex<double> cii(0.0,1.0); | |
| 12 | +static complex<double> czero(0.0,0.0); | |
| 13 | +static complex<double> cone(1.0,0.0); | |
| 14 | + | |
| 15 | +double gamma(double x); | |
| 16 | + | |
| 17 | +int cbessik01(complex<double>z,complex<double>&ci0,complex<double>&ci1, | |
| 18 | + complex<double>&ck0,complex<double>&ck1,complex<double>&ci0p, | |
| 19 | + complex<double>&ci1p,complex<double>&ck0p,complex<double>&ck1p) | |
| 20 | +{ | |
| 21 | + complex<double> z1,z2,zr,zr2,cr,ca,cb,cs,ct,cw; | |
| 22 | + double a0,w0; | |
| 23 | + int k,kz; | |
| 24 | + static double a[] = { | |
| 25 | + 0.125, | |
| 26 | + 7.03125e-2, | |
| 27 | + 7.32421875e-2, | |
| 28 | + 1.1215209960938e-1, | |
| 29 | + 2.2710800170898e-1, | |
| 30 | + 5.7250142097473e-1, | |
| 31 | + 1.7277275025845, | |
| 32 | + 6.0740420012735, | |
| 33 | + 2.4380529699556e1, | |
| 34 | + 1.1001714026925e2, | |
| 35 | + 5.5133589612202e2, | |
| 36 | + 3.0380905109224e3}; | |
| 37 | + static double b[] = { | |
| 38 | + -0.375, | |
| 39 | + -1.171875e-1, | |
| 40 | + -1.025390625e-1, | |
| 41 | + -1.4419555664063e-1, | |
| 42 | + -2.7757644653320e-1, | |
| 43 | + -6.7659258842468e-1, | |
| 44 | + -1.9935317337513, | |
| 45 | + -6.8839142681099, | |
| 46 | + -2.7248827311269e1, | |
| 47 | + -1.2159789187654e2, | |
| 48 | + -6.0384407670507e2, | |
| 49 | + -3.3022722944809e3}; | |
| 50 | + static double a1[] = { | |
| 51 | + 0.125, | |
| 52 | + 0.2109375, | |
| 53 | + 1.0986328125, | |
| 54 | + 1.1775970458984e1, | |
| 55 | + 2.1461706161499e2, | |
| 56 | + 5.9511522710323e3, | |
| 57 | + 2.3347645606175e5, | |
| 58 | + 1.2312234987631e7, | |
| 59 | + 8.401390346421e08, | |
| 60 | + 7.2031420482627e10}; | |
| 61 | + | |
| 62 | + a0 = abs(z); | |
| 63 | + z2 = z*z; | |
| 64 | + z1 = z; | |
| 65 | + if (a0 == 0.0) { | |
| 66 | + ci0 = cone; | |
| 67 | + ci1 = czero; | |
| 68 | + ck0 = complex<double> (1e308,0); | |
| 69 | + ck1 = complex<double> (1e308,0); | |
| 70 | + ci0p = czero; | |
| 71 | + ci1p = complex<double>(0.5,0.0); | |
| 72 | + ck0p = complex<double>(-1e308,0); | |
| 73 | + ck1p = complex<double>(-1e308,0); | |
| 74 | + return 0; | |
| 75 | + } | |
| 76 | + if (real(z) < 0.0) z1 = -z; | |
| 77 | + if (a0 <= 18.0) { | |
| 78 | + ci0 = cone; | |
| 79 | + cr = cone; | |
| 80 | + for (k=1;k<=50;k++) { | |
| 81 | + cr *= 0.25*z2/(double)(k*k); | |
| 82 | + ci0 += cr; | |
| 83 | + if (abs(cr/ci0) < eps) break; | |
| 84 | + } | |
| 85 | + ci1 = cone; | |
| 86 | + cr = cone; | |
| 87 | + for (k=1;k<=50;k++) { | |
| 88 | + cr *= 0.25*z2/(double)(k*(k+1.0)); | |
| 89 | + ci1 += cr; | |
| 90 | + if (abs(cr/ci1) < eps) break; | |
| 91 | + } | |
| 92 | + ci1 *= 0.5*z1; | |
| 93 | + } | |
| 94 | + else { | |
| 95 | + if (a0 >= 50.0) kz = 7; | |
| 96 | + else if (a0 >= 35.0) kz = 9; | |
| 97 | + else kz = 12; | |
| 98 | + ca = exp(z1)/sqrt(2.0*M_PI*z1); | |
| 99 | + ci0 = cone; | |
| 100 | + zr = 1.0/z1; | |
| 101 | + for (k=0;k<kz;k++) { | |
| 102 | + ci0 += a[k]*pow(zr,k+1.0); | |
| 103 | + } | |
| 104 | + ci0 *= ca; | |
| 105 | + ci1 = cone; | |
| 106 | + for (k=0;k<kz;k++) { | |
| 107 | + ci1 += b[k]*pow(zr,k+1.0); | |
| 108 | + } | |
| 109 | + ci1 *= ca; | |
| 110 | + } | |
| 111 | + if (a0 <= 9.0) { | |
| 112 | + cs = czero; | |
| 113 | + ct = -log(0.5*z1)-el; | |
| 114 | + w0 = 0.0; | |
| 115 | + cr = cone; | |
| 116 | + for (k=1;k<=50;k++) { | |
| 117 | + w0 += 1.0/k; | |
| 118 | + cr *= 0.25*z2/(double)(k*k); | |
| 119 | + cs += cr*(w0+ct); | |
| 120 | + if (abs((cs-cw)/cs) < eps) break; | |
| 121 | + cw = cs; | |
| 122 | + } | |
| 123 | + ck0 = ct+cs; | |
| 124 | + } | |
| 125 | + else { | |
| 126 | + cb = 0.5/z1; | |
| 127 | + zr2 = 1.0/z2; | |
| 128 | + ck0 = cone; | |
| 129 | + for (k=0;k<10;k++) { | |
| 130 | + ck0 += a1[k]*pow(zr2,k+1.0); | |
| 131 | + } | |
| 132 | + ck0 *= cb/ci0; | |
| 133 | + } | |
| 134 | + ck1 = (1.0/z1 - ci1*ck0)/ci0; | |
| 135 | + if (real(z) < 0.0) { | |
| 136 | + if (imag(z) < 0.0) { | |
| 137 | + ck0 += cii*M_PI*ci0; | |
| 138 | + ck1 = -ck1+cii*M_PI*ci1; | |
| 139 | + } | |
| 140 | + else if (imag(z) > 0.0) { | |
| 141 | + ck0 -= cii*M_PI*ci0; | |
| 142 | + ck1 = -ck1-cii*M_PI*ci1; | |
| 143 | + } | |
| 144 | + ci1 = -ci1; | |
| 145 | + } | |
| 146 | + ci0p = ci1; | |
| 147 | + ci1p = ci0-1.0*ci1/z; | |
| 148 | + ck0p = -ck1; | |
| 149 | + ck1p = -ck0-1.0*ck1/z; | |
| 150 | + return 0; | |
| 151 | +} | |
| 152 | +int cbessikna(int n,complex<double> z,int &nm,complex<double> *ci, | |
| 153 | + complex<double> *ck,complex<double> *cip,complex<double> *ckp) | |
| 154 | +{ | |
| 155 | + complex<double> ci0,ci1,ck0,ck1,ckk,cf,cf1,cf2,cs; | |
| 156 | + double a0; | |
| 157 | + int k,m,ecode; | |
| 158 | + a0 = abs(z); | |
| 159 | + nm = n; | |
| 160 | + if (a0 < 1.0e-100) { | |
| 161 | + for (k=0;k<=n;k++) { | |
| 162 | + ci[k] = czero; | |
| 163 | + ck[k] = complex<double>(-1e308,0); | |
| 164 | + cip[k] = czero; | |
| 165 | + ckp[k] = complex<double>(1e308,0); | |
| 166 | + } | |
| 167 | + ci[0] = cone; | |
| 168 | + cip[1] = complex<double>(0.5,0.0); | |
| 169 | + return 0; | |
| 170 | + } | |
| 171 | + ecode = cbessik01(z,ci[0],ci[1],ck[0],ck[1],cip[0],cip[1],ckp[0],ckp[1]); | |
| 172 | + if (n < 2) return 0; | |
| 173 | + ci0 = ci[0]; | |
| 174 | + ci1 = ci[1]; | |
| 175 | + ck0 = ck[0]; | |
| 176 | + ck1 = ck[1]; | |