""" J.A.C. Weideman: Computation of the Complex Error Function (SIAM - J. Num. Anal., 31, 1497-1518, 1994) """ from math import pi, sqrt sqrt2 = sqrt(2.0) recSqrtPi = 1./sqrt(pi) try: import numpy as np except ImportError, msg: raise SystemExit (str(msg) + '\nweideman.py: import numpy failed!') from humlicek import humlicek2 # select the functions to be imported (thus 'from weideman import *' will not import the constants!) __all__ = 'weideman16a weideman24a weideman32a hum1wei24a hum1wei32A hum1wei32a hum2wei32A hum2wei32a'.split() #################################################################################################################################### ##### ##### ##### Notes on polynomial evaluation: ##### ##### ##### ##### * poly1d(a,Z) and the explicit Horner scheme used here have approximately the same speed according to %timeit ##### ##### ##### ##### * poly1d(a,Z) does not work if Z (or z) is an array!!! ##### ##### When you first construct the polynomial with polyWeideman = poly1d(a) ##### ##### then you can evaluate the complex error function also with an array, e.g. ##### ##### w = (recSqrtPi + 2.0 * polyWeideman(Z) * recLmiZ) * recLmiZ ##### ##### ##### #################################################################################################################################### L16 = 4.0/sqrt(sqrt2) # sqrt(N)/2^(1/4) = 3.3635856610148585 for N=16 a16 = np.array([ 1.8976999933151775e+00, # a0 = L/sqrtPi 1.7483958860819617e+00, 1.3622408222719591e+00, 8.8644783020505524e-01, 4.6929090090360348e-01, 1.9124172674669498e-01, 5.1822402431611909e-02, 3.6825673170924143e-03, -3.8810151890227547e-03, -1.5276597401218833e-03, 8.7031584284902146e-05, 2.1071056396581287e-04, 2.1709867932466571e-05, -2.7346404624048783e-05, -5.5842334104101195e-06, 3.9812875760042709e-06, 9.9393225361232851e-07]) def weideman16a (z): """ Complex error function --- J.A.C. Weideman: SIAM J. Num. Anal. 31 1497-1518 (1994); equation (38.I) and table I. """ x, y = z.real, z.imag iz = 1j*x - y lpiz = L16 + iz # L - y + x*complex(0.,1.) lmiz = L16 - iz # L + y - x*complex(0.,1.) recLmiZ = 1.0 / lmiz Z = lpiz * recLmiZ polynom = a16[1]+(a16[2]+(a16[3]+(a16[4]+(a16[5]+(a16[6]+(a16[7]+(a16[8]+(a16[9]+(a16[10]+(a16[11]+(a16[12]+(a16[13]+(a16[14]+(a16[15]+a16[16]*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z w = (recSqrtPi + 2.0 * polynom * recLmiZ) * recLmiZ return w #################################################################################################################################### L24 = sqrt(24./sqrt2) # sqrt(N)/2^(1/4) = 4.1195342878142354 a24 = np.array([ 2.3241983342526162e+00, # a0 = L/sqrtPi 2.1978589365315417e+00, 1.8562864992055408e+00, 1.3948196733791203e+00, 9.2570871385886788e-01, 5.3611395357291292e-01, 2.6549639598807689e-01, 1.0838723484566792e-01, 3.3723366855316413e-02, 6.2150063629501763e-03, -4.9364269012806686e-04, -7.8166429956142650e-04, -2.0748431511424456e-04, 2.4331415462641969e-05, 3.0471066083243790e-05, 4.1394617248575527e-06, -3.0388931839840047e-06, -1.0856475790698251e-06, 2.5682641346701115e-07, 1.8738343486619108e-07, -1.9122258522976932e-08, -3.0082822811202271e-08, 1.3310461806370372e-09, 4.9048215867870488e-09, -1.5137461654527820e-10]) def weideman24a (z): """ Complex error function --- J.A.C. Weideman: SIAM J. Num. Anal. 31 1497-1518 (1994); equation (38.I) and table I. """ iz = 1j*z.real - z.imag lpiz = L24 + iz # wL - y + x*complex(0.,1.) lmiz = L24 - iz # wL + y - x*complex(0.,1.) recLmiZ = 1.0 / lmiz Z = lpiz * recLmiZ polynom = a24[1]+(a24[2]+(a24[3]+(a24[4]+(a24[5]+(a24[6]+(a24[7]+(a24[8]+(a24[9]+(a24[10]+(a24[11]+(a24[12]+(a24[13]+(a24[14]+(a24[15]+(a24[16]+(a24[17]+(a24[18]+(a24[19]+(a24[20]+(a24[21]+(a24[22]+(a24[23]+a24[24]*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z w = (recSqrtPi + 2.0 * polynom * recLmiZ) * recLmiZ return w #################################################################################################################################### L32 = sqrt(32./sqrt2) # sqrt(N)/2^(1/4) = 4.1195342878142354 a32 = np.array([ 2.6837530678616557e+00, # a0 = L/sqrtPi 2.5722534081245696e+00, 2.2635372999002676e+00, 1.8256696296324824e+00, 1.3455441692345453e+00, 9.0192548936480144e-01, 5.4601397206393498e-01, 2.9544451071508926e-01, 1.4060716226893769e-01, 5.7304403529837900e-02, 1.9006155784845689e-02, 4.5195411053501429e-03, 3.9259136070122748e-04, -2.4532980269928922e-04, -1.3075449254548613e-04, -2.1409619200870880e-05, 6.8210319440412389e-06, 4.4015317319048931e-06, 4.2558331390536872e-07, -4.1840763666294341e-07, -1.4813078891201116e-07, 2.2930439569075392e-08, 2.3797557105844622e-08, 8.1248960947953431e-10, -3.2080150458594088e-09, -5.2310170266050247e-10, 4.1537465934749353e-10, 1.1658312885903929e-10, -5.5441820344468828e-11, -2.1542618451370239e-11, 8.0314997274316680e-12, 3.7424975634801558e-12, -1.3031797863050087e-12]) def weideman32a (z): """ Complex error function --- J.A.C. Weideman: SIAM J. Num. Anal. 31 1497-1518 (1994); equation (38.I) and table I. """ iz = 1j*z.real - z.imag lpiz = L32 + iz # wL - y + x*complex(0.,1.) lmiz = L32 - iz # wL + y - x*complex(0.,1.) recLmiZ = 1.0 / lmiz Z = lpiz * recLmiZ polynom = a32[1]+(a32[2]+(a32[3]+(a32[4]+(a32[5]+(a32[6]+(a32[7]+(a32[8]+(a32[9]+(a32[10]+(a32[11]+(a32[12]+(a32[13]+(a32[14]+(a32[15]+(a32[16]+(a32[17]+(a32[18]+(a32[19]+(a32[20]+(a32[21]+(a32[22]+(a32[23]+(a32[24]+(a32[25]+(a32[26]+(a32[27]+(a32[28]+(a32[29]+(a32[30]+(a32[31]+a32[32]*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z w = (recSqrtPi + 2.0 * polynom * recLmiZ) * recLmiZ return w #################################################################################################################################### ################################ Combinations of Weideman and Humlicek w4 ################################ #################################################################################################################################### def hum1wei24a (z): """ Complex error function using Humlicek's and Weideman's (N=24) rational approximation: |x|+y>15: Humlicek (JQSRT, 1982) rational approximation for region I; else: J.A.C. Weideman (SIAM-NA 1994); equation (38.I) and table I. """ w = 1j*recSqrtPi *z / (z*z-0.5) # Humlicek (1982) approx 1 for s>15 # any points inside? mask = abs(z.real)+z.imag<15.0 if mask.any(): zSmall = z[mask] # only the subset of z inside np.place(w, mask, weideman24a(zSmall)) return w def hum1wei32a (z): """ Complex error function using Humlicek's and Weideman's (N=32) rational approximation: |x|+y>15: Humlicek (JQSRT, 1982) rational approximation for region I; else: J.A.C. Weideman (SIAM-NA 1994); equation (38.I) and table I. """ w = 1j*recSqrtPi *z / (z*z-0.5) # Humlicek (1982) approx 1 for s>15 # any points inside? mask = abs(z.real)+z.imag<15.0 if mask.any(): zSmall = z[mask] # only the subset of z inside np.place(w, mask, weideman32a(zSmall)) return w def hum2wei32a (z): """ Complex error function using Humlicek's and Weideman's (N=32) rational approximation: |x|+y>8.: Humlicek (JQSRT, 1982) rational approximation for region II; else: J.A.C. Weideman (SIAM-NA 1994); equation (38.I) and table I. """ w = humlicek2(z) # Humlicek (1982) approx 2 # any points inside? mask = abs(z.real)+z.imag<8.0 if mask.any(): zSmall = z[mask] # only the subset of z inside np.place(w, mask, weideman32a(zSmall)) return w #################################################################################################################################### ##### ##### ##### the versions above are slightly faster (the Weideman polynomial is only evaluated for the subset of z) ##### ##### BUT ##### ##### these versions do not work with the "matrix" of Z values (generated with meshgrid) for the contour plots in demo.py ##### ##### ##### ##### Numpy's vectorize can be used to make a "matrix-aware" function from a scalar version of hum1wei24a etc. ##### ##### BUT this is much slower than the versions above! ##### ##### ##### #################################################################################################################################### def hum1wei32A (z, sCut=15.0): """ Complex error function: Humlicek (JQSRT, 1982) region I for |x|+y>sCut and Weideman inside. """ return np.where(abs(z.real)+z.imag>sCut, 1j*recSqrtPi *z / (z*z-0.5), weideman32a(z)) def hum2wei32A (z, sCut=8): """ Complex error function: Humlicek (JQSRT, 1982) region II for |x|+y>sCut and Weideman inside. """ return np.where(abs(z.real)+z.imag>sCut, humlicek2(z), weideman32a(z))