test_bessel.cc 3.5 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677
  1. //**********************************************************************************//
  2. // Copyright (C) 2009-2015 Ovidio Pena <ovidio@bytesfall.com> //
  3. // Copyright (C) 2013-2015 Konstantin Ladutenko <kostyfisik@gmail.com> //
  4. // //
  5. // This file is part of scattnlay //
  6. // //
  7. // This program is free software: you can redistribute it and/or modify //
  8. // it under the terms of the GNU General Public License as published by //
  9. // the Free Software Foundation, either version 3 of the License, or //
  10. // (at your option) any later version. //
  11. // //
  12. // This program is distributed in the hope that it will be useful, //
  13. // but WITHOUT ANY WARRANTY; without even the implied warranty of //
  14. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
  15. // GNU General Public License for more details. //
  16. // //
  17. // The only additional remark is that we expect that all publications //
  18. // describing work using this software, or all commercial products //
  19. // using it, cite the following reference: //
  20. // [1] O. Pena and U. Pal, "Scattering of electromagnetic radiation by //
  21. // a multilayered sphere," Computer Physics Communications, //
  22. // vol. 180, Nov. 2009, pp. 2348-2354. //
  23. // //
  24. // You should have received a copy of the GNU General Public License //
  25. // along with this program. If not, see <http://www.gnu.org/licenses/>. //
  26. //**********************************************************************************//
  27. #include <complex>
  28. #include <cmath>
  29. #include <cstdio>
  30. #include <stdexcept>
  31. #include <vector>
  32. #include "../../bessel.h"
  33. #include "../../nmie.h"
  34. int main(int argc, char *argv[]) {
  35. int n = 10;
  36. std::complex<double> z(1.0, 2.0);
  37. int nm;
  38. std::vector< std::complex<double> > csj, cdj, csy, cdy;
  39. nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy );
  40. auto result = csj;
  41. printf("===========Calculate and compare against Wolfram Alpha\n");
  42. printf("j(0,1+i2) = re(%16.18f)\n im(%16.18f)\n",
  43. real(result[0]),
  44. imag(result[0]));
  45. printf("WA j() = re(1.4169961192118759)\n im(-0.874391197002)\n");
  46. printf("j(1,1+i2) = re(%16.18f)\n im(%16.18f)\n",
  47. real(result[1]),
  48. imag(result[1]));
  49. printf("WA j() = re(0.74785726329830368)\n im(0.68179207555304)\n");
  50. printf("j(4,1+i2) = re(%16.18f)\n im(%16.18f)\n",
  51. real(result[4]),
  52. imag(result[4]));
  53. printf("WA j() = re(-0.01352410550046)\n im(-0.027169663050653)\n");
  54. n = 20;
  55. std::complex<double> z1(0.001, 0.0);
  56. nmie::bessel::csphjy (n, z1, nm, csj, cdj, csy, cdy );
  57. result = csj;
  58. printf("$$$$ REAL $$$$$$ Calculate and compare against Wolfram Alpha\n");
  59. printf("y(0,1) = %16.18f\n", real(result[0]));
  60. printf("y(1,1) = %16.18f\n", real(result[1]));
  61. printf("y(2,1) = %.14g\n", real(result[2]));
  62. }