Updated version of check_wronskian and added check_tmat scripts

utils/check_tmat.py

+#!/usr/bin/env python
+
+from numpy import *
+from matplotlib.pyplot import *
+
+t = loadtxt('tmat_atom_001_energ_001.dat')
+t = t[:,0]+1j*t[:,1]
+t = t.reshape(int(sqrt(len(t))), -1)
+t = mat(t)
+
+with open('tmat_atom_001_energ_001.dat') as f:
+    tmpline = f.readlines()[0]
+    etxt = tmpline.split('=')[-1].split()
+    e = float(etxt[0])+1j*float(etxt[1])
+
+figure(figsize=(6,8))
+subplot(4,2,1); imshow(real(t)); colorbar(); title('Re'); ylabel('t')
+subplot(4,2,2); imshow(imag(t)); colorbar(); title('Im')
+
+dt = (t.transpose().conjugate() - t)/2.
+subplot(4,2,3); imshow(real(dt)); colorbar(); ylabel('lhs=(t^dagger-t)/2')
+subplot(4,2,4); imshow(imag(dt)); colorbar()
+
+dt2 = 1j*sqrt(real(e))*(t.transpose().conjugate() * t)
+subplot(4,2,5); imshow(real(dt2)); colorbar(); ylabel('rhs=i*sqrt(e)*t^dagger*t')
+subplot(4,2,6); imshow(imag(dt2)); colorbar()
+
+dt = dt-dt2
+#dt[abs(dt2)<10**-2]=0
+subplot(4,2,7); imshow(log10(abs(real(dt)))); colorbar(); ylabel('log10(abs(lhs-rhs))')
+subplot(4,2,8); imshow(log10(abs(imag(dt)))); colorbar()
+
+suptitle('Checks for optical theorem')
+subplots_adjust(top=0.891, bottom=0.048, left=0.043, right=0.943, hspace=0.483, wspace=0.174)
+
+t = array(t)
+for il in range(len(t)):
+  print il, log10(abs( imag(t[il,il])-sqrt(real(e))*sum(abs(t[:,il])**2) ))
+
+show()

utils/check_wronskian.py

 #!/usr/bin/env python
 
-from numpy import loadtxt, sqrt, sum
-from matplotlib.pyplot import imshow, colorbar, show, figure, subplot, axis, subplots_adjust, title
+import sys
+from numpy import loadtxt, sqrt, sum, log, abs
+from matplotlib.pyplot import imshow, colorbar, show, figure, subplot, axis, subplots_adjust, title, plot, legend, gca, axvline
 
-figure(figsize=(12,6))
-subplots_adjust(left=0.04, bottom=0.25, right=0.99, top=0.75)
+if len(sys.argv)>1:
+  ir = int(sys.argv[1])
+else:
+  ir = 0
+if len(sys.argv)>2:
+  lin = False
+else:
+  lin = True
 
-# first kind
-subplot(1,2,1)
-title('wronskian of first kind')
-d = sum(loadtxt('test_wronskian')[:,2:], axis=1)
-d = d.reshape(int(sqrt(len(d))),-1)
-imshow(d, interpolation='nearest')
-axis('equal')
-colorbar()
+r = loadtxt('test_rmesh')
+rpan = loadtxt('test_rpan_intervall')
 
-# second kind
-subplot(1,2,2)
-title('wronskian of second kind')
-d = sum(loadtxt('test_wronskian2')[:,2:], axis=1)
-d = d.reshape(int(sqrt(len(d))),-1)
-imshow(d, interpolation='nearest')
-axis('equal')
-colorbar()
+def plot_wavefun(name, ir, lin, tadd, ipl):
+  d = loadtxt('test_'+name)
+
+  figure(4+icomp)
+  subplot(2,2,1+ipl)
+  title(name+tadd)
+  l0 = int(sqrt(len(d)))
+  d1 = d.reshape(l0,l0,-1)
+  if lin:
+    plot(r, d1[0, 0, 2+icomp::2], label='(lm1,lm2)=(1,1)')
+    plot(r, d1[1, 1, 2+icomp::2], label='(lm1,lm2)=(2,2)')
+    plot(r, d1[4, 4, 2+icomp::2], label='(lm1,lm2)=(5,5)')
+  else:
+    plot(r, abs(d1[0, 0, 2+icomp::2]), label='(lm1,lm2)=(1,1)')
+    plot(r, abs(d1[1, 1, 2+icomp::2]), label='(lm1,lm2)=(2,2)')
+    plot(r, abs(d1[4, 4, 2+icomp::2]), label='(lm1,lm2)=(5,5)')
+    ax = gca()
+    ax.set_yscale('log')
+  if ir!=0:
+    axvline(r[abs(ir)], color='r', lw=0.5)
+  legend()
+
+  figure(2+icomp)
+  subplot(2,2,1+ipl)
+  title(name+tadd)
+  if ir<=0:
+    d = sum(d[:,2+icomp+ir::2], axis=1)
+  else:
+    d = d[:,2+ir+icomp]
+  d = d.reshape(int(sqrt(len(d))),-1)
 
-figure()
+  if lin:
+    imshow(d, interpolation='nearest')
+  else:
+    imshow(log(abs(d)), interpolation='nearest')
+
+  axis('equal')
+  colorbar()
+
+
+#figsize=(12,6))
+#subplots_adjust(left=0.04, bottom=0.25, right=0.99, top=0.75)
+
+# first kind
+figure(100)
+d = loadtxt('test_wronskian')
+l0 = int(sqrt(len(d)))
+d1 = d.reshape(l0,l0,-1)
+icomp = 0
+if lin:
+  plot(r, d1[0, 0, 2+icomp::2], label='(lm1,lm2)=(1,1)')
+  plot(r, d1[1, 1, 2+icomp::2], label='(lm1,lm2)=(2,2)')
+  plot(r, d1[4, 4, 2+icomp::2], label='(lm1,lm2)=(5,5)')
+else:
+  #plot(r, abs(d1[0, 0, 2+icomp::2]-d1[0, 0, 2+icomp::2][-1]), label='(lm1,lm2)=(1,1)')
+  #plot(r, abs(d1[1, 1, 2+icomp::2]-d1[1, 1, 2+icomp::2][-1]), label='(lm1,lm2)=(2,2)')
+  #plot(r, abs(d1[4, 4, 2+icomp::2]-d1[4, 4, 2+icomp::2][-1]), label='(lm1,lm2)=(5,5)')
+  plot(r, abs(d1[0, 0, 2+icomp::2]), label='(lm1,lm2)=(1,1)')
+  plot(r, abs(d1[1, 1, 2+icomp::2]), label='(lm1,lm2)=(2,2)')
+  plot(r, abs(d1[4, 4, 2+icomp::2]), label='(lm1,lm2)=(5,5)')
+  ax = gca()
+  ax.set_xlim(r[1], r[-1])
+  ax.set_yscale('log')
+  ax.set_xscale('log')
+if ir!=0:
+  axvline(r[abs(ir)], color='r', lw=0.8)
+for ipan in rpan:
+  axvline(ipan, color='grey', lw=0.5)
+"""
+ax2 = gca().twiny()
+ax2.plot(range(len(r))) # Create a dummy plot
+ax2.cla()
+if not lin:
+  ax2.set_xlim(1, len(r)-1)
+  ax2.set_xscale('log')
+"""
+
+figure(1)
 subplot(2,2,1)
-title('rll')
-d = sum(loadtxt('test_rll')[:,2:], axis=1)
+title('wronskian of first kind')
+if ir<=0:
+  d = sum(d[:,2+ir::2], axis=1)
+else:
+  d = d[:,2+ir]
 d = d.reshape(int(sqrt(len(d))),-1)
-imshow(d, interpolation='nearest')
+if lin:
+  imshow(d, interpolation='nearest')
+else:
+  imshow(log(abs(d)), interpolation='nearest')
 axis('equal')
 colorbar()
+# imag
+figure(1)
 subplot(2,2,2)
-title('rllleft')
-d = sum(loadtxt('test_rllleft')[:,2:], axis=1)
+title('imaginary part')
+if ir<=0:
+  d = sum(loadtxt('test_wronskian')[:,3+ir::2], axis=1)
+else:
+  d = loadtxt('test_wronskian')[:,3+ir]
 d = d.reshape(int(sqrt(len(d))),-1)
-imshow(d, interpolation='nearest')
+if lin:
+  imshow(d, interpolation='nearest')
+else:
+  imshow(log(abs(d)), interpolation='nearest')
 axis('equal')
 colorbar()
+
+# second kind
+figure(1)
 subplot(2,2,3)
-title('sll')
-d = sum(loadtxt('test_sll')[:,2:], axis=1)
+title('wronskian of second kind')
+if ir<=0:
+  d = sum(loadtxt('test_wronskian2')[:,2+ir::2], axis=1)
+else:
+  d = loadtxt('test_wronskian2')[:,2+ir]
 d = d.reshape(int(sqrt(len(d))),-1)
-imshow(d, interpolation='nearest')
+if lin:
+  imshow(d, interpolation='nearest')
+else:
+  imshow(log(abs(d)), interpolation='nearest')
 axis('equal')
 colorbar()
+# imag
+figure(1)
 subplot(2,2,4)
-title('sllleft')
-d = sum(loadtxt('test_sllleft')[:,2:], axis=1)
+title('imaginary part')
+if ir<=0:
+  d = sum(loadtxt('test_wronskian2')[:,3+ir::2], axis=1)
+else:
+  d = loadtxt('test_wronskian2')[:,3+ir]
 d = d.reshape(int(sqrt(len(d))),-1)
-imshow(d, interpolation='nearest')
+if lin:
+  imshow(d, interpolation='nearest')
+else:
+  imshow(log(abs(d)), interpolation='nearest')
 axis('equal')
 colorbar()
 
+#"""
+for icomp in [0,1]: # for real/imag
+  if icomp==0:
+    tadd = ' real part'
+  else:
+    tadd = ' imag part'
+  plot_wavefun('rll', ir, lin, tadd, 0)
+  plot_wavefun('rllleft', ir, lin, tadd, 1)
+  plot_wavefun('sll', ir, lin, tadd, 2)
+  plot_wavefun('sllleft', ir, lin, tadd, 3)
+#"""
+ 
+
 show()