joeydumont / complex_bessel

A C++ library to evaluate Bessel functions of all kinds.
GNU Lesser General Public License v3.0
52 stars 25 forks source link

Caclulating Bessel functions of a first kind with a real and imaginary component yields a complex vice real number #24

Open karlinterestedincomplexbessel opened 3 months ago

karlinterestedincomplexbessel commented 3 months ago

First Kudos to D.E. Amos, and Dumont and McNogginson for making the original Fortran code and supplying a wrapper for c+. Since c / c++ is the only code I know at a moderately advanced level, your work has been of great benefit to me and I hope many others. I did need to translate the Amos code given in the github link into c / c++ as I spent a lot of time trying to get the wrapper to work using Visual Studio Community 2022 (64-bit) Version 19.9.6 but was not able to do so as my understanding of doing so is rather limited.

After completing all subroutines (now functions completely called by reference) required for Bessel Functions of the 1st kind, as well as the partially redundant Hanke1 and 2 functions, I tested the c code (both the translated Fortran as well as the original Fortran. They seem to work almost alright, but I I am having a problem when the argument to these functions is a pure real number. In this case, the return value for bessj() always contains an imaginary component (for example argument 1.328, with fnu = 1 and N = 100). This is a disappointment for me because I have been looking for a robust manner to calculate complex Bessel Functions for some time in order to design laser light scattering program that will work with absorbing particles in order to design and interpret phase Doppler Anamometry experiments for particle sizing. Can anyone offer me an explanation for the behavior I am observing as well as a way around the problem. Thanks again.

joeydumont commented 3 months ago

Hi there!

Could you show me the code you are running? I couldn't reproduce your issue by calling besselJ(100, std::complex<double>(1.328,0)), which results in (1.75755e-176,0), which exactly agrees with WolframAlpha.

karlinterestedincomplexbessel commented 3 months ago

Thank you for responding! I have attached a zipped folder. If you have visual studio, you can decompress and immediately run FortranLearn1\x64\Debug\FortranLearn1.exe which will print out Bessel functions of the first kind from order 100 through 110 (in increments of 1) for z = 1.328. The source code is in the folder C:\Users\karlk\Documents\Visual Studio 2022\Projects\FortranLearn1\Source.

The subroutines are compiled from the file “zbesh.for” which was downloaded from your github files. Very strange that my results are different from what you are getting. I am a novice at fortran.

Results from this program give values for the real part identical to Wolf FortranLearn1.zip ram Mathematica, but also yield an imaginary part. For order 100, the result is 1.75755e-176 + 3.35846e-190i.

joeydumont commented 3 months ago

Ah, you're not using the complex_bessel itself, just your own wrappers around the Fortran code. Unfortunately, I cannot see the source code for the binary that calls the zbesh.for subroutines to compute these values.

Second, while 3.35846e-190i is indeed nonzero, it might as well be. The underlying library does not set the result to zero when the argument is known to return a real value. I suggest you either set a manual cutoff for what you consider "zero", or handle the case there z is real yourself.

I cannot run your compiled binary because I am on Linux, not Windows, and I don't have access to MSVC to compile your code.