hankel.cpp 7.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250
  1. /************************* HANKEL FUNCTIONS **********************/
  2. //#include "stdafx.h"
  3. #include <math.h>
  4. #include "complex.h"
  5. #include "hankel.h"
  6. #define abs(a) ((a)>=0 ? (a) : -(a))
  7. double recurrence( int, int, double, double* );
  8. /***** Function Hank20 - to calculate H^(2)_0(x), x - double ****/
  9. Complex Hank20( double x ) {
  10. double vj0, vy0;
  11. hank01( vj0, vy0, x, 1 );
  12. return( cmplx( vj0, -vy0 ) );
  13. }
  14. /***** Function Hank21 - to calculate H^(2)_1(x), x - double ****/
  15. Complex Hank21( double x ) {
  16. double vj1, vy1;
  17. hank11( vj1, vy1, x, 1 );
  18. return( cmplx( vj1, -vy1 ) );
  19. }
  20. /***** Function Hank10 - to calculate H^(1)_0(x), x - double ****/
  21. Complex Hank10( double x ) {
  22. double vj0, vy0;
  23. hank01( vj0, vy0, x, 1 );
  24. return( cmplx( vj0, vy0 ) );
  25. }
  26. /***** Function Hank11 - to calculate H^(1)_1(x), x - double ****/
  27. Complex Hank11( double x ) {
  28. double vj1, vy1;
  29. hank11( vj1, vy1, x, 1 );
  30. return( cmplx( vj1, vy1 ) );
  31. }
  32. int hank01( double& vj0, double& vy0, double xd, int n ) {
  33. int n1, n2;
  34. double x, y, z, fx, x1, x2, x3, x4;
  35. double xlg = 1.0e+70;
  36. double a[]={0.0,
  37. -.17e-18 , .1222e-16 ,
  38. -.75885e-15 , .4125321e-13 ,
  39. -.194383469e-11 , .7848696314e-10 ,
  40. -.267925353056e-08 , .7608163592419e-07 ,
  41. -.176194690776215e-05 , .3246032882100508e-04 ,
  42. -.46062616620627505e-03 , .48191800694676045e-02,
  43. -.34893769411408885e-01 , .15806710233209726 ,
  44. -.37009499387264978 , .26517861320333681 ,
  45. -.87234423528522213e-02 , .31545594294978024 ,
  46. -.1e-19 , .39e-18 ,
  47. -.2698e-16 , .164349e-14 ,
  48. -.8747341e-13 , .402633082e-11 ,
  49. -.15837552542e-09 , .524879478733e-08 ,
  50. -.14407233274019e-06 , .32065325376548e-05 ,
  51. -.5632079141056987e-04 , .75311359325777423e-03,
  52. -.72879624795520792e-02 , .47196689595763387e-01,
  53. -.17730201278114358 , .26156734625504664 ,
  54. .17903431407718266 ,-.27447430552974527 ,
  55. -.66292226406569883e-01 ,
  56. -.1e-19 , .2e-19 ,
  57. -.11e-18 , .55e-18 ,
  58. -.288e-17 , .1631e-16 ,
  59. -.10012e-15 , .67481e-15 ,
  60. -.506903e-14 , .4326596e-13 ,
  61. -.43045789e-12 , .516826239e-11 ,
  62. -.7864091377e-10 , .163064646352e-08 ,
  63. -.5170594537606e-07 , .307518478751947e-05 ,
  64. -.53652204681321174e-03 , .19989206986950373e01 ,
  65. .1e-19 ,-.3e-19 ,
  66. .13e-18 ,-.62e-18 ,
  67. .311e-17 ,-.1669e-16 ,
  68. .9662e-16 ,-.60999e-15 ,
  69. .425523e-14 ,-.3336328e-13 ,
  70. .30061451e-12 ,-.320674742e-11 ,
  71. .4220121905e-10 ,-.72719159369e-09 ,
  72. .1797245724797e-07 ,-.74144984110606e-06 ,
  73. .683851994261165e-04 ,-.31111709210674018e-01};
  74. x = xd;
  75. y = abs(x);
  76. z = y * .125;
  77. if(z == 0.0) {
  78. vj0 = 1.0;
  79. vj0 =-xlg;
  80. return 0;
  81. }
  82. if(z > 1.0) {
  83. z = 1.0 / z;
  84. x2= 4.0 * z*z - 2.0;
  85. n1= 38;
  86. n2= 55;
  87. x1 = recurrence(n1, n2, x2, a);
  88. n1 = 56;
  89. n2 = 73;
  90. fx = recurrence(n1, n2, x2, a);
  91. x2 = cos(y - 0.7853981633974483);
  92. x3 = sin(y - 0.7853981633974483);
  93. x4 = 0.7978845608028654 / sqrt(y);
  94. fx = fx * z;
  95. vj0= x4 * (x1*x2 - fx*x3);
  96. vy0= x4 * (fx*x2 + x1*x3);
  97. return 0;
  98. }
  99. else {
  100. x2 = 4.0 * z*z - 2.0;
  101. n1 = 1;
  102. n2 = 18;
  103. vj0= recurrence(n1, n2, x2, a);
  104. if(n <= 0) return 0;
  105. n1 = 19;
  106. n2 = 37;
  107. fx = recurrence(n1, n2, x2, a);
  108. vy0= .6366197723675813 * log(y) * vj0 + fx;
  109. return 0;
  110. }
  111. }
  112. int hank11( double& vj1, double& vy1, double xd, int n ) {
  113. int n1, n2;
  114. double x, y, z, fx, x1, x2, x3, x4;
  115. double xlg = 1.0e+70;
  116. double a[] = {0.0,
  117. -.4e-19 , .295e-17 ,
  118. -.19554e-15 , .1138572e-13 ,
  119. -.57774042e-12 , .2528123664e-10 ,
  120. -.94242129816e-9 , .2949707007278e-7 ,
  121. -.76175878054003e-6 , .1588701923993213e-4 ,
  122. -.26044438934858068e-3 , .32402701826838575e-2 ,
  123. -.29175524806154208e-1 , .17770911723972828e0 ,
  124. -.66144393413454325e0 , .12879940988576776e1 ,
  125. -.11918011605412169e1 , .12967175412105298e1 ,
  126. .9e-19 ,-.658e-17 ,
  127. .42773e-15 ,-.2440949e-13 ,
  128. .121143321e-11 ,-.5172121473e-10 ,
  129. .187547032473e-8 ,-.5688440039919e-7 ,
  130. .141662436449235e-5 ,-.283046401495148e-4 ,
  131. .44047862986709951e-3 ,-.51316411610610848e-2 ,
  132. .42319180353336904e-1 ,-.22662499155675492e0 ,
  133. .67561578077218767e0 ,-.76729636288664594e0 ,
  134. -.12869738438135000e0 , .40608211771868508e-1 ,
  135. .1e-19 ,-.2e-19 ,
  136. .12e-18 ,-.58e-18 ,
  137. .305e-17 ,-.1731e-16 ,
  138. .10668e-15 ,-.72212e-15 ,
  139. .545267e-14 ,-.4684224e-13 ,
  140. .46991955e-12 ,-.570486364e-11 ,
  141. .881689866e-10 ,-.187189074911e-8 ,
  142. .6177633960644e-7 ,-.398728430048891e-5 ,
  143. .89898983308594085e-3 , .20018060817200274e1 ,
  144. -.1e-19 , .3e-19 ,
  145. -.14e-18 , .65e-18 ,
  146. -.328e-17 , .1768e-16 ,
  147. -.10269e-15 , .65083e-15 ,
  148. -.456125e-14 , .3596777e-13 ,
  149. -.32643157e-12 , .351521879e-11 ,
  150. -.4686363688e-10 , .82291933277e-9 ,
  151. -.2095978138408e-7 , .91386152579555e-6 ,
  152. -.9627723549157079e-4 , .93555574139070650e-1};
  153. x = xd;
  154. y = abs(x);
  155. z = y * .125;
  156. if(z <= 0.0) {
  157. vj1 = 0.0;
  158. vj1 =-xlg;
  159. return 0;
  160. }
  161. if(z > 1) {
  162. z = 1.0 / z;
  163. x2= 4.0 * z*z - 2.0;
  164. n1= 37;
  165. n2= 54;
  166. x1= recurrence(n1, n2, x2, a);
  167. n1 = 55;
  168. n2 = 72;
  169. fx = recurrence(n1, n2, x2, a);
  170. x2 = cos(y - 2.356194490192345);
  171. x3 = sin(y - 2.356194490192345);
  172. x4 = 0.7978845608028654 / sqrt(y);
  173. fx = fx * z;
  174. vj1= x4 *(x1 * x2 - fx * x3);
  175. vy1= x4 *(fx * x2 + x1 * x3);
  176. return 0;
  177. }
  178. else {
  179. x2 = 4.0 * z*z - 2.0;
  180. n1 = 1;
  181. n2 = 18;
  182. fx = recurrence(n1, n2, x2, a);
  183. vj1= fx * z;
  184. if(n <= 0) return 0;
  185. n1 = 19;
  186. n2 = 36;
  187. fx = recurrence(n1, n2, x2, a);
  188. vy1= 0.6366197723675813 * (log(y) * vj1 - 1.0 / y) + fx * z;
  189. return 0;
  190. }
  191. }
  192. double recurrence( int n1, int n2, double x2, double* a ) {
  193. int i;
  194. double q1, q2 = 0.0, q3 = 0.0;
  195. for(i = n1; i <= n2; i++) {
  196. q1 = q2;
  197. q2 = q3;
  198. q3 = x2 * q2 - q1 + a[i];
  199. }
  200. return (q3 - q1) * 0.5;
  201. }