special-functions-test-generator.py 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130
  1. #!/usr/bin/env python3
  2. import mpmath as mp
  3. import numpy as np
  4. Du_test = [
  5. # // x, [Re(m), Im(m)], Qext, Qsca, test_name
  6. [0.099, [0.75,0], 7.417859e-06, 7.417859e-06, 'a'],
  7. [0.101, [0.75,0], 8.033538e-06, 8.033538e-06, 'b'],
  8. [10, [0.75,0], 2.232265, 2.232265, 'c'],
  9. [1000, [0.75,0], 1.997908, 1.997908, 'd'],
  10. [100, [1.33,-1e-5], 2.101321, 2.096594, 'e'],
  11. [10000, [1.33,-1e-5], 2.004089, 1.723857, 'f'],
  12. [0.055, [1.5, -1], 0.101491, 1.131687e-05, 'g'],
  13. [0.056, [1.5, -1], 0.1033467, 1.216311e-05, 'h'],
  14. [100, [1.5, -1], 2.097502, 1.283697, 'i'],
  15. [10000, [1.5, -1], 2.004368, 1.236574, 'j'],
  16. [1, [10, -10], 2.532993, 2.049405, 'k'],
  17. [100, [10, -10,], 2.071124, 1.836785, 'l'],
  18. [10000, [10, -10], 2.005914, 1.795393, 'm'],
  19. [80, [1.05, 1], 0, 0, 'Yang'],
  20. [1, [mp.pi, 1], 0, 0, 'pi'],
  21. [1, [mp.pi, -1], 0, 0, 'pi'],
  22. [1, [mp.pi, mp.pi], 0, 0, 'pi'],
  23. [1, [2*mp.pi, -1], 0, 0, 'pi'],
  24. [1, [2*mp.pi, mp.pi], 0, 0, 'pi'],
  25. [1, [2*mp.pi, 1], 0, 0, 'pi'],
  26. ]
  27. # // Dtest refractive index is m={1.05,1}, the size parameter is x = 80
  28. n_list = [0,1,30,50,60,70,75,80,85,90,99,116,130];
  29. def get_z_values(du_list):
  30. zlist = []
  31. for record in du_list:
  32. x = mp.mpf(str(record[0]))
  33. m = mp.mpc(str(record[1][0]), str(record[1][1]))
  34. z = x*m
  35. zlist.append(z)
  36. return zlist
  37. # APPLIED OPTICS / Vol. 53, No. 31 / 1 November 2014, eq(13)
  38. def LeRu_cutoff(z):
  39. x = mp.fabs(z)
  40. return int(x + 11 * x**(1/3) + 1)
  41. def get_n_list(z, max_number_of_elements = 10):
  42. nmax = LeRu_cutoff(z)
  43. factor = nmax**(1/(max_number_of_elements-2))
  44. n_list = [int(factor**i) for i in range(max_number_of_elements-1) ]
  45. n_list.append(0)
  46. n_set = set(n_list)
  47. return sorted(n_set)
  48. # Riccati-Bessel z*j_n(z)
  49. def psi(n,z):
  50. return mp.sqrt( (mp.pi * z)/2 ) * mp.autoprec(mp.besselj)(n+1/2,z)
  51. # to compare r(n,z) with Wolfram Alpha
  52. # n=49, z=1.3-2.1i, SphericalBesselJ[n-1,z]/SphericalBesselJ[n,z]
  53. def r(n,z):
  54. if n > 0:
  55. return psi(n-1,z)/psi(n,z)
  56. return mp.cos(z)/mp.sin(z)
  57. def D1(n,z):
  58. return r(n,z) - n/z
  59. def get_test_data_nlist(z_record, output_dps, n):
  60. x = str(z_record[0])
  61. mr = str(z_record[1][0])
  62. mi = str(z_record[1][1])
  63. z_str = ''
  64. try:
  65. z = mp.mpf(x)*mp.mpc(mr, mi)
  66. D1nz = D1(n,z)
  67. z_str=('{{'+
  68. mp.nstr(z.real, output_dps*2)+','+
  69. mp.nstr(z.imag, output_dps*2)+'},'+
  70. str(n)+',{'+
  71. mp.nstr(D1nz.real, output_dps)+','+
  72. mp.nstr(D1nz.imag, output_dps)+'},'+
  73. mp.nstr(mp.fabs(D1nz.real* 10**-output_dps),2)+',' +
  74. mp.nstr(mp.fabs(D1nz.imag* 10**-output_dps),2)+
  75. '},\n')
  76. except:
  77. pass
  78. return z_str
  79. def get_test_data(Du_test, output_dps, max_num_elements_of_n_list):
  80. output_str = ('// complex(z), n, complex(D1_n(z)), abs_err_real, abs_err_imag\n'+
  81. 'std::vector< std::tuple< std::complex<double>, int, std::complex<double>, double, double > >'+
  82. 'D1_test_'+str(output_dps)+'digits = {\n')
  83. for z_record in Du_test:
  84. x = str(z_record[0])
  85. mr = str(z_record[1][0])
  86. mi = str(z_record[1][1])
  87. mp.mp.dps = 20
  88. z = mp.mpf(x)*mp.mpc(mr, mi)
  89. n_list = get_n_list(z, max_num_elements_of_n_list)
  90. print(z, n_list)
  91. for n in n_list:
  92. mp.mp.dps = 20
  93. old_z_string = get_test_data_nlist(z_record, output_dps, n)
  94. mp.mp.dps = 37
  95. new_z_string = get_test_data_nlist(z_record, output_dps, n)
  96. while old_z_string != new_z_string:
  97. new_dps = int(mp.mp.dps * 1.41)
  98. if new_dps > 300: break
  99. mp.mp.dps = new_dps
  100. print("New dps = ", mp.mp.dps)
  101. old_z_string = new_z_string
  102. new_z_string = get_test_data_nlist(z_record, output_dps, n)
  103. output_str += new_z_string
  104. output_str += '};\n'
  105. return output_str
  106. def main ():
  107. output_dps = 10
  108. max_num_elements_of_nlist = 15
  109. out_filename = 'test_spec_functions_data.h'
  110. output_str = get_test_data(Du_test, output_dps, max_num_elements_of_nlist)
  111. with open(out_filename, 'w') as out_file:
  112. out_file.write(output_str)
  113. main()