Jupyter QtConsole 4.3.1

Python 3.6.9 (default, Oct 29 2019, 10:39:36) [GCC]

Type 'copyright', 'credits' or 'license' for more information

IPython 6.2.1 -- An enhanced Interactive Python. Type '?' for help.


In [1]: %cd /run/media/franz/DAC4-8310/

/run/media/franz/DAC4-8310


In [2]: from scipy.special import wofz

   ...: from scipy import pi

   ...: recPi = 1/pi

   ...:

   ...: from cpfX import cpf12a, cpf12b, cpf12w, cpfX, zpf16p, hum1zpf16w, hum2zpf16w, hum1zpf16m, hum2zpf16m

   ...: from racef import weideman24a, hum1wei24, hum2wei32, hum1wei24w, hum2wei32w

   ...:

   ...:

   ...: from matplotlib.pyplot import *

   ...: from matplotlib import ticker # possibly its already imported automtically

   ...: contourLevels = np.logspace(-12,4,17)

   ...: set_cmap('RdBu_r') # if you have already closed the figure this will give a RuntimeError, ignore it

   ...: rc('figure.subplot', top=0.95, bottom=0.07, left=0.08, right=0.975, wspace=0.15)


In [3]: def relErrorContour (X, Y, wTest, wRef, iml=False, errorY=False, cBar=True, yLabel=True):

   ...: """ Color contour plot of relative function difference (default real part of complex functions). """

   ...: if iml: # plot error of imaginary part

   ...: relError = abs(wTest.imag-wRef.imag)/wRef.imag

   ...: maxRelErrorOther = np.max(abs(wTest.real-wRef.real)/wRef.real)

   ...: cBarLabel = r'$\Delta L/L_{\rm ref}$'

   ...: else:

   ...: relError = abs(wTest.real-wRef.real)/wRef.real

   ...: ijMax = np.argmax(relError); iMax=ijMax%relError.shape[1]; jMax=int(ijMax/relError.shape[1])

   ...: print ('%s %g @ %i %i %i x = %.3f y = %.1g' %

   ...: ('max rel error:', np.max(relError), ijMax, iMax, jMax, X[0,iMax], Y[jMax,0]))

   ...: cBarLabel = r'$\Delta K/K_{\rm ref}$'

   ...: # max rel. error of imaginary part

   ...: relErrorIm = np.where(abs(wRef.imag)>0, abs(wTest.imag-wRef.imag)/wRef.imag, 0.0)

   ...: maxRelErrorOther = np.max(relErrorIm)

   ...:

   ...: if errorY:

   ...: axes([0.05,0.1,0.99,0.85]) # axes([0.05,0.1,0.75,0.85])

   ...: cPad=0.25

   ...: else:

   ...: cPad=0.02

   ...:

   ...: contourf(X,np.log10(Y), relError, contourLevels, locator=ticker.LogLocator())

   ...: xlabel('$x$')

   ...: if yLabel: ylabel(r'$\log_{10}(y)$')

   ...: if cBar: colorbar(label=cBarLabel, pad=cPad, extend='both') # pad=0.01

   ...: title ('max error %.3e (%.2e)' % (np.max(relError), maxRelErrorOther))

   ...:

   ...: if errorY:

   ...: maxRelErrorY = np.max(relError,1)

   ...: logY = np.log10(Y[:,0])

   ...: axes([0.67,0.1,0.19,0.85]) # axes([0.78,0.1,0.2,0.85])

   ...: semilogx (maxRelErrorY,logY); ylim(logY.min(),logY.max()); yticks(visible=0); xlabel(cBarLabel)

   ...: return maxRelErrorY


In [4]: rcParams['figure.figsize']=[12.,8.]


In [5]: xx=np.linspace(0.,25.,251)

   ...: yy=np.logspace(-8,2,101)

   ...:

   ...: # kind of 'matrices' required for the contour plots (note Python is case sensitive!)

   ...: XX, YY = np.meshgrid(xx, yy)

   ...: ZZ = XX +1j*YY

   ...:

   ...: # the reference code

   ...: wRef = wofz(ZZ)

   ...: kRef = wRef.real

   ...:

   ...: w_zpf16 = zpf16p(ZZ)

   ...: w_hum1zpf16 = hum1zpf16w(ZZ)

   ...: w_hum2zpf16 = hum2zpf16w(ZZ)

   ...: w_hum1wei24 = hum1wei24w(XX,YY) # note the two arguments!

   ...: w_hum2wei32 = hum2wei32w(XX,YY)


In [6]: # plot the real part of the complex error function for various y

   ...: for j,y in enumerate(yy[::10]): semilogy (xx, w_zpf16[10*j,:].real, label='%i %g' % (j, y))

   ...: xlabel('$x$'); ylabel(r'$K(x,y)$'); legend()

Out[6]: <matplotlib.legend.Legend at 0x7f95f249f630>


In [7]: # plot relative errors as contour plot and for selected y

   ...: wTest = w_hum1wei24

   ...: kTest = wTest.real

   ...: subplot(121)

   ...: for j,y in enumerate(yy[::10]):

   ...: kDelta = abs(kTest[10*j,:]-kRef[10*j,:])/kRef[10*j,:]

   ...: semilogy (xx, kDelta, label = '%.1e %.3g' % (y, max(kDelta)))

   ...: xlabel('$x$'); ylabel(r'$\Delta K/K_{\rm ref}$')

   ...: legend()

   ...: subplot(122)

   ...: relErrorContour(XX,YY, wTest, wRef)

   ...: suptitle('humlicek1 weideman24', fontsize='x-large')

max rel error: 0.258756 @ 126 126 0 x = 12.600 y = 1e-08


Out[7]: <matplotlib.text.Text at 0x7f95f028a630>


In [8]: # compare relative errors for several approximations

   ...: for j,y in enumerate(yy[::-20]):

   ...: subplot(6,1,j+1)

   ...: kRef = wofz(xx+1j*y).real

   ...: semilogy (xx, abs(zpf16p(xx+1j*y).real-kRef)/kRef, label='zpf16')

   ...: # milogy (xx, abs(hum1zpf16w(xx+1j*y).real-kRef)/kRef, label='hum1zpf16w')

   ...: semilogy (xx, abs(hum1zpf16m(xx, y).real-kRef)/kRef, label='hum1zpf16m')

   ...: # milogy (xx, abs(hum1wei24w(xx, y).real-kRef)/kRef, label='hum1wei24w')

   ...: semilogy (xx, abs(hum1wei24 (xx, y).real-kRef)/kRef, label='hum1wei24m')

   ...: semilogy (xx, abs(hum2wei32w(xx, y).real-kRef)/kRef, label='hum2wei32')

   ...: title('$y$ = %.1e' % y)

   ...: xlabel('$x$'); ylabel(r'$\Delta K/K_{\rm ref}$'); legend()

Out[8]: <matplotlib.legend.Legend at 0x7f95ee5800f0>


In [9]: