test_bessel.cc 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321
  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_small = 1, n_big = 15;
  36. int n = n_big+2, n_print;
  37. std::complex<double> z_small(0.0012, 0.0034);
  38. std::complex<double> z_medium(1.0, 2.0);
  39. std::complex<double> z_big(18.0, 17.0);
  40. std::complex<double> z_big_real(81.0, 0.0);
  41. std::complex<double> z_big_imag(0.0, 13.0);
  42. std::complex<double> z;
  43. int nm;
  44. // csj and cdj - complex argument spherical bessel first kind and derivative
  45. // csy and cdy - complex argument spherical bessel second kind and derivative
  46. std::vector< std::complex<double> > csj, cdj, csy, cdy, csh;
  47. std::vector<std::complex<double> > Psi, Zeta, dPsi, dZeta;
  48. auto c_i = std::complex<double>(0.0,1.0);
  49. printf("===== Small ========\n");
  50. z = z_small;
  51. printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
  52. nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy );
  53. csh.resize(csj.size());
  54. for (unsigned int i = 0; i < csj.size(); ++i)
  55. csh[i] = csj[i]+c_i*csy[i];
  56. nmie::bessel::calcPsi(n, z, Psi, dPsi);
  57. nmie::bessel::calcZeta(n, z, Zeta, dZeta);
  58. // n_print = n_small;
  59. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print,
  60. // csj[n_print].real(), csj[n_print].imag(), n_print,
  61. // csy[n_print].real(), csy[n_print].imag());
  62. // printf("WA csj = e r e i; WA csy = e r e i\n");
  63. // n_print = n_big;
  64. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big,
  65. // csj[n_print].real(), csj[n_print].imag(), n_print,
  66. // csy[n_print].real(), csy[n_print].imag());
  67. // printf(" WA csj = e r e i; WA csy = e r e i\n");
  68. // n_print = n_small;
  69. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  70. // csh[n_print].real(), csh[n_print].imag());
  71. // printf("WA csj = e r e i;\n");
  72. // n_print = n_big;
  73. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  74. // csh[n_print].real(), csh[n_print].imag());
  75. // printf("WA csj = e r e i;\n");
  76. // n_print = n_small;
  77. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
  78. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  79. // Zeta[n_print].real(), Zeta[n_print].imag());
  80. // printf("WA Psi = e r e i; WA Zeta = e r e i\n");
  81. // n_print = n_big;
  82. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
  83. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  84. // Zeta[n_print].real(), Zeta[n_print].imag());
  85. // printf(" WA Psi = e r e i; WA Zeta = e r e i\n");
  86. n_print = n_small;
  87. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
  88. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  89. dZeta[n_print].real(), dZeta[n_print].imag());
  90. printf("WA dPsi = 0.800005e-03r+0.226667e-02i; WA dZeta = +4.8284 e+04r-5.98822e+04i\n");
  91. n_print = n_big;
  92. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
  93. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  94. dZeta[n_print].real(), dZeta[n_print].imag());
  95. printf(" WA dPsi = +1.75388e-53r-6.94429e-54i; WA dZeta = +8.58553e+55r+7.47397e+55i\n");
  96. printf("===== Medium ========\n");
  97. z = z_medium;
  98. printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
  99. nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy );
  100. csh.resize(csj.size());
  101. for (unsigned int i = 0; i < csj.size(); ++i)
  102. csh[i] = csj[i]+c_i*csy[i];
  103. nmie::bessel::calcPsi(n, z, Psi, dPsi);
  104. nmie::bessel::calcZeta(n, z, Zeta, dZeta);
  105. // n_print = n_small;
  106. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print,
  107. // csj[n_print].real(), csj[n_print].imag(), n_print,
  108. // csy[n_print].real(), csy[n_print].imag());
  109. // printf("WA csj = e r e i; WA csy = e r e i\n");
  110. // n_print = n_big;
  111. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big,
  112. // csj[n_print].real(), csj[n_print].imag(), n_print,
  113. // csy[n_print].real(), csy[n_print].imag());
  114. // printf(" WA csj = e r e i; WA csy = e r e i\n");
  115. // n_print = n_small;
  116. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  117. // csh[n_print].real(), csh[n_print].imag());
  118. // printf("WA csj = e r e i;\n");
  119. // n_print = n_big;
  120. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  121. // csh[n_print].real(), csh[n_print].imag());
  122. // printf("WA csj = e r e i;\n");
  123. // n_print = n_small;
  124. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
  125. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  126. // Zeta[n_print].real(), Zeta[n_print].imag());
  127. // printf("WA Psi = e r e i; WA Zeta = e r e i\n");
  128. // n_print = n_big;
  129. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
  130. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  131. // Zeta[n_print].real(), Zeta[n_print].imag());
  132. // printf(" WA Psi = e r e i; WA Zeta = e r e i\n");
  133. n_print = n_small;
  134. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
  135. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  136. dZeta[n_print].real(), dZeta[n_print].imag());
  137. printf("WA dPsi = 2.41792e r+1.27781 e i; WA dZeta = +1.99423e-01r-7.01483e-02i\n");
  138. n_print = n_big;
  139. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
  140. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  141. dZeta[n_print].real(), dZeta[n_print].imag());
  142. printf(" WA dPsi = -1.03313e-11r-1.13267e-11i; WA dZeta = -2.11402e+11r+8.34528e+10i\n");
  143. printf("===== Big ========\n");
  144. z = z_big;
  145. printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
  146. nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy );
  147. csh.resize(csj.size());
  148. for (unsigned int i = 0; i < csj.size(); ++i)
  149. csh[i] = csj[i]+c_i*csy[i];
  150. nmie::bessel::calcPsi(n, z, Psi, dPsi);
  151. nmie::bessel::calcZeta(n, z, Zeta, dZeta);
  152. // n_print = n_small;
  153. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print,
  154. // csj[n_print].real(), csj[n_print].imag(), n_print,
  155. // csy[n_print].real(), csy[n_print].imag());
  156. // printf("WA csj = e r e i; WA csy = e r e i\n");
  157. // n_print = n_big;
  158. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big,
  159. // csj[n_print].real(), csj[n_print].imag(), n_print,
  160. // csy[n_print].real(), csy[n_print].imag());
  161. // printf(" WA csj = e r e i; WA csy = e r e i\n");
  162. // n_print = n_small;
  163. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  164. // csh[n_print].real(), csh[n_print].imag());
  165. // printf("WA csj = e r e i;\n");
  166. // n_print = n_big;
  167. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  168. // csh[n_print].real(), csh[n_print].imag());
  169. // printf("WA csj = e r e i;\n");
  170. // n_print = n_small;
  171. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
  172. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  173. // Zeta[n_print].real(), Zeta[n_print].imag());
  174. // printf("WA Psi = e r e i; WA Zeta = e r e i\n");
  175. // n_print = n_big;
  176. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
  177. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  178. // Zeta[n_print].real(), Zeta[n_print].imag());
  179. // printf(" WA Psi = e r e i; WA Zeta = e r e i\n");
  180. n_print = n_small;
  181. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
  182. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  183. dZeta[n_print].real(), dZeta[n_print].imag());
  184. printf("WA dPsi = e r e i; WA dZeta = -3.11025e-08r-2.90558e-08i\n");
  185. n_print = n_big;
  186. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
  187. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  188. dZeta[n_print].real(), dZeta[n_print].imag());
  189. printf(" WA dPsi = -2.47740e+05r+3.05899e+05i; WA dZeta = -6.13497e-07r-1.14967e-06i\n");
  190. printf("===== Big real========\n");
  191. z = z_big_real;
  192. printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
  193. nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy );
  194. csh.resize(csj.size());
  195. for (unsigned int i = 0; i < csj.size(); ++i)
  196. csh[i] = csj[i]+c_i*csy[i];
  197. nmie::bessel::calcPsi(n, z, Psi, dPsi);
  198. nmie::bessel::calcZeta(n, z, Zeta, dZeta);
  199. // n_print = n_small;
  200. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print,
  201. // csj[n_print].real(), csj[n_print].imag(), n_print,
  202. // csy[n_print].real(), csy[n_print].imag());
  203. // printf("WA csj = e r e i; WA csy = e r e i\n");
  204. // n_print = n_big;
  205. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big,
  206. // csj[n_print].real(), csj[n_print].imag(), n_print,
  207. // csy[n_print].real(), csy[n_print].imag());
  208. // printf(" WA csj = e r e i; WA csy = e r e i\n");
  209. // n_print = n_small;
  210. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  211. // csh[n_print].real(), csh[n_print].imag());
  212. // printf("WA csj = e r e i;\n");
  213. // n_print = n_big;
  214. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  215. // csh[n_print].real(), csh[n_print].imag());
  216. // printf("WA csj = e r e i;\n");
  217. // n_print = n_small;
  218. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
  219. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  220. // Zeta[n_print].real(), Zeta[n_print].imag());
  221. // printf("WA Psi = e r e i; WA Zeta = e r e i\n");
  222. // n_print = n_big;
  223. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
  224. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  225. // Zeta[n_print].real(), Zeta[n_print].imag());
  226. // printf(" WA Psi = e r e i; WA Zeta = e r e i\n");
  227. n_print = n_small;
  228. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
  229. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  230. dZeta[n_print].real(), dZeta[n_print].imag());
  231. printf("WA dPsi = e r e i; WA dZeta = -6.20203e-01r-7.84344e-01i\n");
  232. n_print = n_big;
  233. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
  234. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  235. dZeta[n_print].real(), dZeta[n_print].imag());
  236. printf(" WA dPsi = e r e i; WA dZeta = -7.13982e-01r+6.86858e-01i\n");
  237. printf("===== Big imag========\n");
  238. z = z_big_imag;
  239. printf("z = %+10.5er%+10.5ei\n", z.real(), z.imag());
  240. nmie::bessel::csphjy (n, z, nm, csj, cdj, csy, cdy );
  241. csh.resize(csj.size());
  242. for (unsigned int i = 0; i < csj.size(); ++i)
  243. csh[i] = csj[i]+c_i*csy[i];
  244. nmie::bessel::calcPsi(n, z, Psi, dPsi);
  245. nmie::bessel::calcZeta(n, z, Zeta, dZeta);
  246. // n_print = n_small;
  247. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_print,
  248. // csj[n_print].real(), csj[n_print].imag(), n_print,
  249. // csy[n_print].real(), csy[n_print].imag());
  250. // printf("WA csj = e r e i; WA csy = e r e i\n");
  251. // n_print = n_big;
  252. // printf("csj[%i] = %+10.5er%+10.5ei; csy[%i] = %+10.5er%+10.5ei\n", n_big,
  253. // csj[n_print].real(), csj[n_print].imag(), n_print,
  254. // csy[n_print].real(), csy[n_print].imag());
  255. // printf(" WA csj = e r e i; WA csy = e r e i\n");
  256. // n_print = n_small;
  257. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  258. // csh[n_print].real(), csh[n_print].imag());
  259. // printf("WA csj = e r e i;\n");
  260. // n_print = n_big;
  261. // printf("csh[%i] = %+10.5er%+10.5ei\n", n_print,
  262. // csh[n_print].real(), csh[n_print].imag());
  263. // printf("WA csj = e r e i;\n");
  264. // n_print = n_small;
  265. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_print,
  266. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  267. // Zeta[n_print].real(), Zeta[n_print].imag());
  268. // printf("WA Psi = e r e i; WA Zeta = e r e i\n");
  269. // n_print = n_big;
  270. // printf("Psi[%i] = %+10.5er%+10.5ei; Zeta[%i] = %+10.5er%+10.5ei\n", n_big,
  271. // Psi[n_print].real(), Psi[n_print].imag(), n_print,
  272. // Zeta[n_print].real(), Zeta[n_print].imag());
  273. // printf(" WA Psi = e r e i; WA Zeta = e r e i\n");
  274. n_print = n_small;
  275. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_print,
  276. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  277. dZeta[n_print].real(), dZeta[n_print].imag());
  278. printf("WA dPsi = e r e i; WA dZeta = -2.47233e-22r-2.44758e-06i\n");
  279. n_print = n_big;
  280. printf("dPsi[%i] = %+10.5er%+10.5ei; dZeta[%i] = %+10.5er%+10.5ei\n", n_big,
  281. dPsi[n_print].real(), dPsi[n_print].imag(), n_print,
  282. dZeta[n_print].real(), dZeta[n_print].imag());
  283. printf(" WA dPsi = e r e i; WA dZeta = 8.39449e-19r+1.28767e-02i\n");
  284. }