kelvin.cpp 8.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395
  1. /****Calculate the value of Kelvin Function****/
  2. #include "stdafx.h"
  3. extern double Pi;
  4. extern Complex Ic;
  5. Complex seta(double);
  6. Complex dphi(double);
  7. //-------------------------------------------
  8. // Pi = 4.0 * atan(1.0);
  9. // s_pi_4 = Pi / 4.0;
  10. //-------------------------------------------
  11. static double s_pi_4 = atan(1.0);
  12. static double s_sq_2 = sqrt(2.0);
  13. // return the modified Bessel function I0(x) for real x
  14. double ber(double x) {
  15. double ber;
  16. if(fabs(x) <= 8.0) {
  17. double r2, r3, r4, r5, r6, r7;
  18. double x8; //,p1,p2,p3,p4,p5,p6,p7,p8;
  19. /****
  20. p1 = 1.0;
  21. p2 = 64.0;
  22. p3 = 113.77777774;
  23. p4 = 32.36345652;
  24. p5 = 2.64191397;
  25. p6 = 0.08349609;
  26. p7 = 0.00122552;
  27. p8 = 0.00000901;
  28. ********/
  29. x8 = x/8.0;
  30. r2 = x8 * x8 * x8 * x8;
  31. r3 = r2 * r2;
  32. r4 = r3 * r2;
  33. r5 = r4 * r2;
  34. r6 = r5 * r2;
  35. r7 = r6 * r2;
  36. ber = 1.0 - 64.0*r2 + 113.77777774*r3 -
  37. 32.36345652*r4 + 2.64191397*r5 -
  38. 0.08349609*r6 + 0.00122552*r7 - 0.00000901*r7*r2;
  39. } else {
  40. cerr << "Ber function error" << endl;
  41. }
  42. return (ber);
  43. }
  44. double bei(double x) {
  45. double bei;
  46. if(fabs(x) <= 8.0) {
  47. double tp, r1, r2, r3, r4, r5, r6;
  48. double x8; //,p1,p2,p3,p4,p5,p6,p7;
  49. /****
  50. p1 = 16.0;
  51. p2 = 113.77777774;
  52. p3 = 72.81777742;
  53. p4 = 10.56765779;
  54. p5 = 0.52185615;
  55. p6 = 0.01103667;
  56. p7 = 0.00011346;
  57. *****/
  58. x8 = x/8.0;
  59. r1 = x8 * x8;
  60. tp = x8 * x8 * x8 * x8;
  61. r2 = r1 * tp;
  62. r3 = r2 * tp;
  63. r4 = r3 * tp;
  64. r5 = r4 * tp;
  65. r6 = r5 * tp;
  66. bei = 16.0*r1 - 113.77777774*r2 + 72.81777742*r3 -
  67. 10.56765779*r4 + 0.52185615*r5 - 0.01103667*r6 +
  68. 0.00011346*r6*tp;
  69. } else {
  70. cerr << "Bei function error" << endl;
  71. }
  72. return(bei);
  73. }
  74. double ker(double x) {
  75. double ker,c1;
  76. if(x <= 0.0) {
  77. cerr << "Ker function error" << endl;
  78. } else if(x <= 8.0) {
  79. double x8,c1; // ,p1,p2,p3,p4,p5,p6,p7,p8;
  80. double r2, r3, r4, r5, r6, r7;
  81. /******
  82. p1 = 0.57721566;
  83. p2 = 59.05819744;
  84. p3 = 171.36272133;
  85. p4 = 60.60977451;
  86. p5 = 5.65539121;
  87. p6 = 0.19636347;
  88. p7 = 0.00309699;
  89. p8 = 0.00002458;
  90. ******/
  91. x8 = x/8.0;
  92. r2 = x8 * x8 * x8 * x8;
  93. r3 = r2 * r2;
  94. r4 = r3 * r2;
  95. r5 = r4 * r2;
  96. r6 = r5 * r2;
  97. r7 = r6 * r2;
  98. c1 = log(x/2.0);
  99. ker = -c1*ber(x) + s_pi_4*bei(x) - 0.57721566 -
  100. 59.05819744*r2 + 171.36272133*r3 - 60.60977451*r4 +
  101. 5.65539121*r5 - 0.19636347*r6 + 0.00309699*r7 -
  102. 0.00002458*r7*r2;
  103. } else {
  104. Complex fx,d1;
  105. c1 = sqrt(Pi/(2.0*x));
  106. d1 = (1.0+Ic)/s_sq_2;
  107. fx = c1*exp(-d1*x+seta(-x));
  108. ker = real(fx);
  109. }
  110. return(ker);
  111. }
  112. double kei( double x ) {
  113. double kei,c1;
  114. Complex fx,d1;
  115. if(x <= 0.0) {
  116. cerr << "Kei function error" << endl;
  117. } else if (x <= 8.0) {
  118. double x8; // ,p1,p2,p3,p4,p5,p6,p7;
  119. double tp, r1, r2, r3, r4, r5, r6;
  120. /********
  121. p1 = 6.76454936;
  122. p2 = 142.91827687;
  123. p3 = 124.23569650;
  124. p4 = 21.30060904;
  125. p5 = 1.17509064;
  126. p6 = 0.02695875;
  127. p7 = 0.00029532;
  128. *****/
  129. x8 = x/8.0;
  130. r1 = x8 * x8;
  131. tp = x8 * x8 * x8 * x8;
  132. r2 = r1 * tp;
  133. r3 = r2 * tp;
  134. r4 = r3 * tp;
  135. r5 = r4 * tp;
  136. r6 = r5 * tp;
  137. c1 = log(x/2.0);
  138. kei = -c1*bei(x) - s_pi_4*ber(x) + 6.76454936*r1 -
  139. 142.91827687*r2 + 124.23569650*r3 - 21.30060904*r4 +
  140. 1.17509064*r5 - 0.02695875*r6 + 0.00029532*r6*tp;
  141. } else {
  142. c1 = sqrt(Pi/(2.0*x));
  143. d1 = (1.0+Ic)/s_sq_2;
  144. fx = c1*exp(-d1*x+seta(-x));
  145. kei = imag(fx);
  146. }
  147. return(kei);
  148. }
  149. double berp(double x) {
  150. double berp;
  151. if(fabs(x) <= 8.0) {
  152. double x8; // ,p1,p2,p3,p4,p5,p6,p7;
  153. double r1, tp, r2, r3, r4, r5, r6;
  154. /*******
  155. p1 = 4.0;
  156. p2 = 14.22222222;
  157. p3 = 6.06814810;
  158. p4 = 0.66047849;
  159. p5 = 0.02609253;
  160. p6 = 0.00045957;
  161. p7 = 0.00000394;
  162. *****/
  163. x8 = x/8.0;
  164. r1 = x8 * x8;
  165. tp = x8 * x8 * x8 * x8;
  166. r2 = r1 * tp;
  167. r3 = r2 * tp;
  168. r4 = r3 * tp;
  169. r5 = r4 * tp;
  170. r6 = r5 * tp;
  171. berp = x*(r1*(-4.0) + 14.22222222*r2 - 6.06814810*r3 +
  172. 0.66047849*r4 - 0.02609253*r5 +
  173. 0.00045957*r6 - 0.00000394*r6*tp);
  174. } else {
  175. cerr << "Berp function error." << endl;
  176. }
  177. return(berp);
  178. }
  179. double beip(double x) {
  180. double beip;
  181. if(fabs(x) <= 8.0) {
  182. double x8; // ,p1,p2,p3,p4,p5,p6,p7;
  183. double r2, r3, r4, r5, r6;
  184. /******
  185. p1 = 0.5;
  186. p2 = 10.66666666;
  187. p3 = 11.37777772;
  188. p4 = 2.31167514;
  189. p5 = 0.14677204;
  190. p6 = 0.00379386;
  191. p7 = 0.00004609;
  192. ******/
  193. x8 = x/8.0;
  194. r2 = x8 * x8 * x8 * x8;
  195. r3 = r2 * r2;
  196. r4 = r3 * r2;
  197. r5 = r4 * r2;
  198. r6 = r5 * r2;
  199. beip = x*(0.5 - 10.66666666*r2 + 11.37777772*r3 -
  200. 2.31167514*r4 + 0.14677204*r5 - 0.00379386*r6 +
  201. 0.00004609*r6*r2);
  202. } else {
  203. cerr << "Beip function error" << endl;
  204. }
  205. return(beip);
  206. }
  207. double kerp(double x) {
  208. double kerp,c1;
  209. if(x <= 0.0) {
  210. cerr << "Kerp function error" << endl;
  211. } else if(x <= 8.0) {
  212. double x8,c1,c2; // ,p1,p2,p3,p4,p5,p6,p7;
  213. double r1, tp, r2, r3, r4, r5, r6;
  214. /******
  215. p1 = 3.69113734;
  216. p2 = 21.42034017;
  217. p3 = 11.36433272;
  218. p4 = 1.41384780;
  219. p5 = 0.06136358;
  220. p6 = 0.00116137;
  221. p7 = 0.00001075;
  222. *****/
  223. x8 = x/8.0;
  224. r1 = x8 * x8;
  225. tp = r1 * r1;
  226. r2 = r1 * tp;
  227. r3 = r2 * tp;
  228. r4 = r3 * tp;
  229. r5 = r4 * tp;
  230. r6 = r5 * tp;
  231. c1 = log(x/2.0);
  232. c2 = 1.0/x;
  233. kerp = -c1*berp(x) - c2*ber(x) + s_pi_4*beip(x) + x*(r1*(-3.69113734) +
  234. 21.42034017*r2 - 11.36433272*r3 + 1.41384780*r4 -
  235. 0.06136358*r5 + 0.00116137*r6 - 0.00001075*r6*tp);
  236. } else {
  237. Complex fx,d1;
  238. c1 = sqrt(Pi/(2.0*x));
  239. d1 = (1.0+Ic)/s_sq_2;
  240. fx = -c1*exp(-d1*x+seta(-x))*dphi(-x);
  241. kerp = real(fx);
  242. }
  243. return(kerp);
  244. }
  245. double keip(double x) {
  246. double keip,c1;
  247. if(x <= 0.0) {
  248. cerr << "Keip function error" << endl;
  249. } else if(x <= 8.0) {
  250. double x8,c2; // ,p1,p2,p3,p4,p5,p6,p7;
  251. double r2, r3, r4, r5, r6;
  252. /*******
  253. p1 = 0.21139217;
  254. p2 = 13.39858846;
  255. p3 = 19.41182758;
  256. p4 = 4.65950823;
  257. p5 = 0.33049424;
  258. p6 = 0.00926707;
  259. p7 = 0.00011997;
  260. ******/
  261. x8 = x/8.0;
  262. c1 = log(x/2.0);
  263. c2 = 1.0/x;
  264. r2 = x8 * x8 * x8 * x8;
  265. r3 = r2 * r2;
  266. r4 = r3 * r2;
  267. r5 = r4 * r2;
  268. r6 = r5 * r2;
  269. keip = -c1*beip(x) - c2*bei(x) - s_pi_4*berp(x) + x*(0.21139217 -
  270. 13.39858846*r2 + 19.41182758*r3 -
  271. 4.65950823*r4 + 0.33049424*r5 - 0.00926707*r6 +
  272. 0.00011997*r6*r2);
  273. } else {
  274. Complex fx,d1;
  275. c1 = sqrt(Pi/(2.0*x));
  276. d1 = (1.0+Ic)/s_sq_2;
  277. fx = -c1*exp(-d1*x+seta(-x))*dphi(-x);
  278. keip = imag(fx);
  279. }
  280. return(keip);
  281. }
  282. static Complex seta_p1 = cmplx(0.0,-0.3926991);
  283. static Complex seta_p2 = cmplx(0.0110486,0.0);
  284. static Complex seta_p3 = cmplx(0.0,-0.0009765);
  285. static Complex seta_p4 = cmplx(-0.0000906,-0.0000901);
  286. static Complex seta_p5 = cmplx(-0.0000252,0.0);
  287. static Complex seta_p6 = cmplx(-0.0000034,0.0000051);
  288. static Complex seta_p7 = cmplx(0.0000006,.0000019);
  289. static Complex dphi_p1 = cmplx(0.7071068,0.7071068);
  290. static Complex dphi_p2 = cmplx(-0.0625001,-0.0000001);
  291. static Complex dphi_p3 = cmplx(-0.0013813,0.0013811);
  292. static Complex dphi_p4 = cmplx(0.0000005,0.0002452);
  293. static Complex dphi_p5 = cmplx(0.0000346,0.0000338);
  294. static Complex dphi_p6 = cmplx(0.0000117,-0.0000024);
  295. static Complex dphi_p7 = cmplx(0.0000016,-0.0000032);
  296. Complex seta(double x) {
  297. Complex seta;
  298. if(x == 0.0) {
  299. cerr << "Seta function error" << endl;
  300. } else {
  301. double x8, r3, r4, r5, r6;
  302. x8 = 8.0/x;
  303. r3 = x8 * x8;
  304. r4 = r3 * x8;
  305. r5 = r4 * x8;
  306. r6 = r5 * x8;
  307. seta = seta_p1 + seta_p2*x8 + seta_p3*r3 + seta_p4*r4 + seta_p5*r5
  308. + seta_p6*r6 + seta_p7*r6*x8;
  309. }
  310. return(seta);
  311. }
  312. Complex dphi(double x) {
  313. Complex dphi;
  314. if(x == 0.0) {
  315. cerr << "dphi function error" << endl;
  316. } else {
  317. double x8, r3, r4, r5, r6;
  318. x8 = 8.0/x;
  319. r3 = x8 * x8;
  320. r4 = r3 * x8;
  321. r5 = r4 * x8;
  322. r6 = r5 * x8;
  323. dphi = dphi_p1 + dphi_p2*x8 + dphi_p3*r3 + dphi_p4*r4 + dphi_p5*r5 +
  324. dphi_p6*r6 + dphi_p7*r6*x8;
  325. }
  326. return(dphi);
  327. }