cornu.py 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
  1. from math import *
  2. # implementation adapted from cephes
  3. def polevl(x, coef):
  4. ans = coef[-1]
  5. for i in range(len(coef) - 2, -1, -1):
  6. ans = ans * x + coef[i]
  7. return ans
  8. sn = [
  9. -2.99181919401019853726E3,
  10. 7.08840045257738576863E5,
  11. -6.29741486205862506537E7,
  12. 2.54890880573376359104E9,
  13. -4.42979518059697779103E10,
  14. 3.18016297876567817986E11
  15. ]
  16. sn.reverse()
  17. sd = [
  18. 1.00000000000000000000E0,
  19. 2.81376268889994315696E2,
  20. 4.55847810806532581675E4,
  21. 5.17343888770096400730E6,
  22. 4.19320245898111231129E8,
  23. 2.24411795645340920940E10,
  24. 6.07366389490084639049E11
  25. ]
  26. sd.reverse()
  27. cn = [
  28. -4.98843114573573548651E-8,
  29. 9.50428062829859605134E-6,
  30. -6.45191435683965050962E-4,
  31. 1.88843319396703850064E-2,
  32. -2.05525900955013891793E-1,
  33. 9.99999999999999998822E-1
  34. ]
  35. cn.reverse()
  36. cd = [
  37. 3.99982968972495980367E-12,
  38. 9.15439215774657478799E-10,
  39. 1.25001862479598821474E-7,
  40. 1.22262789024179030997E-5,
  41. 8.68029542941784300606E-4,
  42. 4.12142090722199792936E-2,
  43. 1.00000000000000000118E0
  44. ]
  45. cd.reverse()
  46. fn = [
  47. 4.21543555043677546506E-1,
  48. 1.43407919780758885261E-1,
  49. 1.15220955073585758835E-2,
  50. 3.45017939782574027900E-4,
  51. 4.63613749287867322088E-6,
  52. 3.05568983790257605827E-8,
  53. 1.02304514164907233465E-10,
  54. 1.72010743268161828879E-13,
  55. 1.34283276233062758925E-16,
  56. 3.76329711269987889006E-20
  57. ]
  58. fn.reverse()
  59. fd = [
  60. 1.00000000000000000000E0,
  61. 7.51586398353378947175E-1,
  62. 1.16888925859191382142E-1,
  63. 6.44051526508858611005E-3,
  64. 1.55934409164153020873E-4,
  65. 1.84627567348930545870E-6,
  66. 1.12699224763999035261E-8,
  67. 3.60140029589371370404E-11,
  68. 5.88754533621578410010E-14,
  69. 4.52001434074129701496E-17,
  70. 1.25443237090011264384E-20
  71. ]
  72. fd.reverse()
  73. gn = [
  74. 5.04442073643383265887E-1,
  75. 1.97102833525523411709E-1,
  76. 1.87648584092575249293E-2,
  77. 6.84079380915393090172E-4,
  78. 1.15138826111884280931E-5,
  79. 9.82852443688422223854E-8,
  80. 4.45344415861750144738E-10,
  81. 1.08268041139020870318E-12,
  82. 1.37555460633261799868E-15,
  83. 8.36354435630677421531E-19,
  84. 1.86958710162783235106E-22
  85. ]
  86. gn.reverse()
  87. gd = [
  88. 1.00000000000000000000E0,
  89. 1.47495759925128324529E0,
  90. 3.37748989120019970451E-1,
  91. 2.53603741420338795122E-2,
  92. 8.14679107184306179049E-4,
  93. 1.27545075667729118702E-5,
  94. 1.04314589657571990585E-7,
  95. 4.60680728146520428211E-10,
  96. 1.10273215066240270757E-12,
  97. 1.38796531259578871258E-15,
  98. 8.39158816283118707363E-19,
  99. 1.86958710162783236342E-22
  100. ]
  101. gd.reverse()
  102. def fresnel(xxa):
  103. x = abs(xxa)
  104. x2 = x * x
  105. if x2 < 2.5625:
  106. t = x2 * x2
  107. ss = x * x2 * polevl(t, sn) / polevl(t, sd)
  108. cc = x * polevl(t, cn) / polevl(t, cd)
  109. elif x > 36974.0:
  110. ss = 0.5
  111. cc = 0.5
  112. else:
  113. t = pi * x2
  114. u = 1.0 / (t * t)
  115. t = 1.0 / t
  116. f = 1.0 - u * polevl(u, fn) / polevl(u, fd)
  117. g = t * polevl(u, gn) / polevl(u, gd)
  118. t = pi * .5 * x2
  119. c = cos(t)
  120. s = sin(t)
  121. t = pi * x
  122. cc = 0.5 + (f * s - g * c) / t
  123. ss = 0.5 - (f * c + g * s) / t
  124. if xxa < 0:
  125. cc = -cc
  126. ss = -ss
  127. return ss, cc
  128. def eval_cornu(t):
  129. spio2 = sqrt(pi * .5)
  130. s, c = fresnel(t / spio2)
  131. s *= spio2
  132. c *= spio2
  133. return s, c