1 # this code is really a mess ...
5 from RandomArray
import *
6 from LinearAlgebra
import *
16 return math
.sqrt(reduce(lambda x
, y
: x
+y
*y
, v
, 0))
18 def blochvector(l
, k
):
19 assert k
>= 0 and k
<= l
/2
20 v
= cos(2*math
.pi
/l
*(k
)*arange(l
))
23 def localvector(l
, pos
):
28 def local2vector(l
, pos1
, pos2
):
36 def diag_anderson(l
, w
):
38 """diagonalize the anderson model"""
42 ham
= zeros((l
, l
), Float
)
47 evalues
, evectors
= eigenvectors(ham
)
48 sorter
= argsort(evalues
)
49 return take(evalues
, sorter
), take(evectors
, sorter
)
53 def diag_harper_full(l
, lam
):
55 """diagonalize the harper model by full diagonalization"""
59 while result
[-1] < tonumber
:
60 result
.append(result
[-2]+result
[-1])
63 l2
, ltest
= fibs(l
)[-2:]
64 if ltest
!= l
: raise ValueError("l isn't a fibonacci number")
66 pot
=lam
*cos((2*pi
*l2
*arange(l
))/l
)
67 ham
= zeros((l
, l
), Float
)
72 evalues
, evectors
= eigenvectors(ham
)
73 sorter
= argsort(evalues
)
74 return take(evalues
, sorter
), take(evectors
, sorter
)
78 def diag_harper(l
, lam
):
80 """diagonalize the harper model by diagonalizing the symmetric and antisymmetic parts"""
84 while result
[-1] < tonumber
:
85 result
.append(result
[-2]+result
[-1])
88 def createsymm(dest
, source
):
92 for i
in xrange(1, (l
+1)/2):
93 dest
[i
] = dest
[l
-i
] = source
[i
] / sqrt2
96 dest
[l
/2] = source
[l
/2]
97 for i
in xrange(1, l
/2):
98 dest
[i
] = dest
[l
-i
] = source
[i
] / sqrt2
100 def createasymm(dest
, source
):
104 for i
in xrange(1, (l
+1)/2):
105 dest
[i
] = source
[i
-1] / sqrt2
106 dest
[l
-i
] = -source
[i
-1] / sqrt2
108 dest
[0] = dest
[l
/2] = 0
109 for i
in xrange(1, l
/2):
110 dest
[i
] = -source
[i
-1] / sqrt2
111 dest
[l
-i
] = source
[i
-1] / sqrt2
113 l2
, ltest
= fibs(l
)[-2:]
114 if ltest
!= l
: raise ValueError("l isn't a fibonacci number")
116 # symmetric and antisymmetric diagonalization
119 symmsize
= (l
-1)/2 + 1
120 pot
= lam
*cos((2*pi
*l2
*arange(symmsize
))/l
)
121 pot
[symmsize
- 1] += 1
122 symmham
= zeros((symmsize
, symmsize
), Float
)
123 for i
in xrange(symmsize
):
124 symmham
[i
, i
] = pot
[i
]
126 symmham
[i
, i
-1] = symmham
[i
-1, i
] = 1
128 symmham
[i
, i
-1] = symmham
[i
-1, i
] = sqrt2
132 pot
= lam
*cos((2*pi
*l2
*arange(symmsize
))/l
)
133 symmham
= zeros((symmsize
, symmsize
), Float
)
134 for i
in xrange(symmsize
):
135 symmham
[i
, i
] = pot
[i
]
136 if i
> 1 and i
< symmsize
- 1:
137 symmham
[i
, i
-1] = symmham
[i
-1, i
] = 1
139 symmham
[i
, i
-1] = symmham
[i
-1, i
] = sqrt2
140 symmevalues
, symmevectors
= eigenvectors(symmham
)
145 pot
= lam
*cos((2*pi
*l2
*(arange(asymmsize
)+1))/l
)
146 pot
[asymmsize
-1] -= 1
147 asymmham
= zeros((asymmsize
, asymmsize
), Float
)
148 for i
in xrange(asymmsize
):
149 asymmham
[i
, i
] = pot
[i
]
151 asymmham
[i
, i
-1] = asymmham
[i
-1, i
] = 1
155 pot
= lam
*cos((2*pi
*l2
*(arange(asymmsize
)+1))/l
)
156 asymmham
= zeros((asymmsize
, asymmsize
), Float
)
157 for i
in xrange(asymmsize
):
158 asymmham
[i
, i
] = pot
[i
]
160 asymmham
[i
, i
-1] = asymmham
[i
-1, i
] = 1
161 asymmevalues
, asymmevectors
= eigenvectors(asymmham
)
163 # build complete solution
164 symmsort
= argsort(symmevalues
)
165 asymmsort
= argsort(asymmevalues
)
167 evalues
= zeros((l
,), Float
)
168 evectors
= zeros((l
, l
), Float
)
170 if j
!= symmsize
and (k
== asymmsize
or
171 symmevalues
[symmsort
[j
]] < asymmevalues
[asymmsort
[k
]]):
172 evalues
[i
] = symmevalues
[symmsort
[j
]]
173 createsymm(evectors
[i
], symmevectors
[symmsort
[j
]])
176 evalues
[i
] = asymmevalues
[asymmsort
[k
]]
177 createasymm(evectors
[i
], asymmevectors
[asymmsort
[k
]])
179 return evalues
, evectors
183 def test_harper(l
, lam
):
185 """compare diag_harper and diag_harper_full"""
187 res1
= diag_harper(l
, lam
)
188 res2
= diag_harper_full(l
, lam
)
189 for x1
, x2
in zip(res1
[0], res2
[0]):
190 assert math
.fabs(x1
-x2
) < 1e-8, "eigenvalues differ"
191 for v1
, v2
in zip(res1
[1], res2
[1]):
193 for x1
, x2
in zip(v1
, v2
):
196 for x1
, x2
in zip(v1
, v2
):
197 assert math
.fabs(x1
-x2
) < 1e-8, "eigenvectors differ\n%s\n%s" % (v1
, v2
)
199 for x1
, x2
in zip(v1
, v2
):
200 assert math
.fabs(x1
+x2
) < 1e-8, "eigenvectors differ\n%s\n%s" % (v1
, -v2
)
206 def wigner_slow(vector
):
208 """create wigner function of a state (direct way)"""
211 wig
= zeros((l
, l
), Float
)
213 # wigner function (direct way)
214 for x0loop
in xrange(l
):
216 for k0loop
in xrange(l
):
218 k0
= (k0loop
-0.5*l
+0.5)*2*pi
/l
220 k0
= (k0loop
-0.5*l
+1)*2*pi
/l
222 for yloop
in xrange(l
):
223 y
= yloop
- l
/2 + 1 - l
%2
224 v
= lambda x
: vector
[(x
+10*l
)%l
]
226 v1
= 0.5*(v(x0
-(y
-1)/2) + v(x0
-(y
+1)/2))
227 v2
= 0.5*(v(x0
+(y
-1)/2) + v(x0
+(y
+1)/2))
233 sum += v1
* v2
* exp(1j
*k0
*y
)
234 assert abs(sum.imag
) < 1e-10, "imaginary part should be zero"
235 wig
[x0loop
, k0loop
] = sum.real
243 """create wigner function of a state (fft)"""
246 wig
= zeros((l
, l
), Float
)
248 fftarray
= zeros(l
, Float
)
249 fftresult
= zeros(l
, Complex
)
250 wig
= zeros((l
, l
), Float
)
251 for x0loop
in xrange(l
):
253 for yloop
in xrange(l
):
254 y
= yloop
- l
/2 + 1 - l
%2
255 v
= lambda x
: vector
[(x
+10*l
)%l
]
257 v1
= 0.5*(v(x0
-(y
-1)/2) + v(x0
-(y
+1)/2))
258 v2
= 0.5*(v(x0
+(y
-1)/2) + v(x0
+(y
+1)/2))
264 fftarray
[(y
+10*l
)%l
] = v1
* v2
265 fftresult
= real_fft(fftarray
)
266 for k0loop
in xrange(l
):
268 index
= int(abs(k0loop
-0.5*l
+0.5))
270 index
= int(abs(k0loop
-0.5*l
+1))
271 wig
[x0loop
, k0loop
] = fftresult
[index
].real
277 def test_wigner(vector
):
279 """compare wigner_slow and wigner"""
281 res1
= wigner_slow(vector
)
282 res2
= wigner(vector
)
283 for v1
, v2
in zip(res1
, res2
):
284 for x1
, x2
in zip(v1
, v2
):
285 assert math
.fabs(x1
-x2
) < 1e-8, "wigner function differs\n%s\n%s" % (res1
, res2
)
287 # test_wigner(diag_anderson(10, 1)[1][5])
288 # test_wigner(diag_anderson(11, 1)[1][5])
291 def husimi_from_wigner(wig
):
293 """calculate the husimi function out of the wigner function"""
298 hus
= zeros((l
, l
), Float
)
301 sys
.stderr
.write("wigner -> husimi (very slow) ...%6.2f%%\r" % ((100.0*(x1
*l
+k1
))/l
/l
))
306 if dx
< -l
/2: dx
+= l
309 if dk
< -l
/2: dk
+= l
312 hus
[x1
, k1
] += 2.0/l
/l
* wig
[x2
, k2
] * exp(-0.5*dx
*dx
/sigma
/sigma
-2*sigma
*sigma
*dk
*dk
)
313 sys
.stderr
.write("wigner -> husimi (very slow) ... done. \n")
319 def husimi_slow(vector
):
322 hus
= zeros((l
, l
), Float
)
323 phases
= zeros((l
, l
), Float
)
325 factor
=1/sqrt(sqrt(2*pi
)*sigma
*l
)
327 for x0loop
in xrange(l
):
329 for k0loop
in xrange(l
):
331 k0
= (k0loop
-0.5*l
+0.5)*2*pi
/l
333 k0
= (k0loop
-0.5*l
+1)*2*pi
/l
337 sum += vector
[x
] * factor
* exp(-(x
-x0
-l
)*(x
-x0
-l
)/(4*sigma
*sigma
)) * phase
338 sum += vector
[x
] * factor
* exp(-(x
-x0
)*(x
-x0
)/(4*sigma
*sigma
)) * phase
339 sum += vector
[x
] * factor
* exp(-(x
-x0
+l
)*(x
-x0
+l
)/(4*sigma
*sigma
)) * phase
340 hus
[x0loop
, k0loop
] = abs(sum)*abs(sum)
341 phases
[x0loop
, k0loop
] = math
.atan2(sum.imag
, sum.real
)
350 hus
= zeros((l
, l
), Float
)
352 factor
=1/sqrt(sqrt(2*pi
)*sigma
*l
)
354 fftarray
= zeros(l
, Complex
)
355 fftresult
= zeros(l
, Complex
)
356 heights
= zeros((l
, l
), Float
)
357 phases
= zeros((l
, l
), Float
)
358 for x0loop
in xrange(l
):
360 for xloop
in xrange(l
):
362 while (x
-x0
< -0.5*l
): x
+= l
363 while (x
-x0
> 0.5*l
): x
-=l
364 fftarray
[xloop
] = vector
[xloop
] * factor
* exp(-(x
-x0
)*(x
-x0
)/(4*sigma
*sigma
))
365 #fftresult = real_fft(fftarray)
366 fftresult
= fft(fftarray
)
367 for k0loop
in xrange(l
):
369 index
= (int(k0loop
-0.5*l
+0.5) + 10*l
)%l
370 #index = int(abs(k0loop-0.5*l+0.5))
373 index
= int(abs(k0loop
-0.5*l
+1))
374 heights
[x0loop
, k0loop
] = abs(fftresult
[index
])*abs(fftresult
[index
])
375 phases
[x0loop
, k0loop
] = math
.atan2(fftresult
[index
].imag
, fftresult
[index
].real
)
376 # if 2*x0loop > l: phases[x0loop, k0loop] = -phases[x0loop, k0loop] + math.pi
377 # if 2*k0loop > l: phases[x0loop, k0loop] = -phases[x0loop, k0loop]
379 return heights
, phases
383 def test_husimi(vector
):
385 """compare husimi_slow and husimi"""
387 res1
= husimi_slow(vector
)
388 res2
= husimi(vector
)
389 for v1
, v2
in zip(res1
, res2
):
390 for x1
, x2
in zip(v1
, v2
):
391 assert math
.fabs(x1
-x2
) < 1e-8, "husimi function differs\n%s\n%s" % (res1
, res2
)
393 # test_husimi(diag_anderson(10, 1)[1][5])
394 # test_husimi(diag_anderson(11, 1)[1][5])
399 def plot_grid(grid
, filename
):
401 class griddata(datafile
._datafile
):
403 def __init__(self
, grid
):
405 for i
in range(len(grid
)):
406 for j
in range(len(grid
[i
])):
408 k
= (j
-0.5*l
+0.5)*2*pi
/l
410 k
= (j
-0.5*l
+1)*2*pi
/l
411 if l
%2 or j
< len(grid
[i
]) - 1:
412 data
.append([grid
[i
, j
], i
+0.5, i
+1.5, k
-math
.pi
/l
, k
+math
.pi
/l
])
414 data
.append([grid
[i
, j
], i
+0.5, i
+1.5, -math
.pi
, -math
.pi
+math
.pi
/l
])
415 data
.append([grid
[i
, j
], i
+0.5, i
+1.5, math
.pi
-math
.pi
/l
, math
.pi
])
416 datafile
._datafile
.__init
__(self
, ["color", "xmin", "xmax", "kmin", "kmax"], data
)
420 mypainter
= graph
.axispainter(baselineattrs
=None, zerolineattrs
=None, innerticklengths
=0, outerticklengths
="0.2")
423 t
= c
.insert(tex
.tex())
424 g
= c
.insert(graph
.graphxy(t
, width
=10, height
=10, x2
=None, y2
=None, backgroundattrs
=canvas
.stroked(),
425 x
=graph
.linaxis(min=0.5, max=len(grid
)+0.5, title
="$x$",
426 part
=graph
.manualpart(ticks
=("1", str(len(grid
))), texts
=("1", "$L$")),
428 y
=graph
.linaxis(min=-pi
, max=pi
, divisor
=pi
, suffix
=r
"\pi", title
="$k$",
429 part
=graph
.linpart("1"), painter
=mypainter
)))
430 g
.plot(graph
.data(hd
, color
="color", xmin
="xmin", xmax
="xmax", ymin
="kmin", ymax
="kmax"),
431 graph
.rect(color
.gradient
.ReverseRainbow
))
434 c
.writetofile(filename
)
439 theta
= 15 * math
.pi
/180
440 phi
= -25 * math
.pi
/180
443 distance
= 150*allscale
444 gridscale
= 1*allscale
450 def __init__(self
, x_3d
, y_3d
, z_3d
, x_2d
, y_2d
, depth
, phase
):
460 return "(%s, %s, %s), (%s, %s), %s, %s " % (self
.x_3d
, self
.y_3d
, self
.z_3d
, self
.x_2d
, self
.y_2d
, self
.depth
, self
.phase
)
463 f
= 1/math
.sqrt(reduce(lambda x
, y
: x
+ y
*y
, list, 0))
464 return map(lambda x
, f
=f
: x
*f
, list)
467 a
= (math
.sin(phi
), -math
.cos(phi
), 0)
468 b
= (-math
.cos(phi
)*math
.sin(theta
), -math
.sin(phi
)*math
.sin(theta
), math
.cos(theta
))
469 eye
= (distance
*math
.cos(phi
)*math
.cos(theta
), distance
*math
.sin(phi
)*math
.cos(theta
), distance
*math
.sin(theta
))
473 def point_3d(self
, x
, y
, z
, phase
=0):
474 d0
= float(a
[0]*b
[1]*(z
-eye
[2]) + a
[2]*b
[0]*(y
-eye
[1]) + a
[1]*b
[2]*(x
-eye
[0])
475 -a
[2]*b
[1]*(x
-eye
[0]) - a
[0]*b
[2]*(y
-eye
[1]) - a
[1]*b
[0]*(z
-eye
[2]))
476 da
= (eye
[0]*b
[1]*(z
-eye
[2]) + eye
[2]*b
[0]*(y
-eye
[1]) + eye
[1]*b
[2]*(x
-eye
[0])
477 -eye
[2]*b
[1]*(x
-eye
[0]) - eye
[0]*b
[2]*(y
-eye
[1]) - eye
[1]*b
[0]*(z
-eye
[2]))
478 db
= (a
[0]*eye
[1]*(z
-eye
[2]) + a
[2]*eye
[0]*(y
-eye
[1]) + a
[1]*eye
[2]*(x
-eye
[0])
479 -a
[2]*eye
[1]*(x
-eye
[0]) - a
[0]*eye
[2]*(y
-eye
[1]) - a
[1]*eye
[0]*(z
-eye
[2]))
480 dc
= (a
[0]*b
[1]*eye
[2] + a
[2]*b
[0]*eye
[1] + a
[1]*b
[2]*eye
[0]
481 -a
[2]*b
[1]*eye
[0] - a
[0]*b
[2]*eye
[1] - a
[1]*b
[0]*eye
[2])
482 return gridpoint(x
, y
, z
, da
/d0
, db
/d0
, -dc
/d0
, phase
)
484 def point_grid(self
, x
, y
, z
, phase
):
485 x
, y
, z
= x
- planecenter
[0], y
- planecenter
[1], z
- planecenter
[2]
486 return self
.point_3d(x
* gridscale
, y
* gridscale
, z
* zscale
, phase
)
490 values
, vectors
= diag_anderson(l
, 0.5)
492 heights
, phases
= husimi(vector
)
493 planecenter
= (0.5*(len(heights
) - 1), 0.5*(len(heights
[0]) - 1), 0)
494 datagrid
= [[None for j
in range(len(heights
[i
]))] for i
in range(len(heights
))]
498 for i
in range(len(heights
)):
499 for j
in range(len(heights
[i
])):
500 datagrid
[i
][j
] = mygrid
.point_grid(i
, j
, heights
[i
][j
]*1000, phases
[i
][j
])
502 def triangle(p1
, p2
, p3
, hue
):
503 ax
, ay
, az
= p2
.x_3d
- p1
.x_3d
, p2
.y_3d
- p1
.y_3d
, p2
.z_3d
- p1
.z_3d
504 bx
, by
, bz
= p3
.x_3d
- p1
.x_3d
, p3
.y_3d
- p1
.y_3d
, p3
.z_3d
- p1
.z_3d
505 axb
= norm(ay
*bz
- az
*by
, az
*bx
- ax
*bz
, ax
*by
-ay
*bx
)
506 e
= norm(eye
[0]-p1
.x_3d
, eye
[1]-p1
.y_3d
, eye
[2]-p1
.z_3d
)
507 axb_r
= axb
[0]*e
[0] + axb
[1]*e
[1] + axb
[2]*e
[2]
508 depthpoint
= mygrid
.point_3d((p1
.x_3d
+ p2
.x_3d
+ p3
.x_3d
)/3.0,
509 (p1
.y_3d
+ p2
.y_3d
+ p3
.y_3d
)/3.0, 0)
510 return [depthpoint
.depth
,
512 path
.path(path
.moveto(p1
.x_2d
, p1
.y_2d
),
513 path
.lineto(p2
.x_2d
, p2
.y_2d
),
514 path
.lineto(p3
.x_2d
, p3
.y_2d
),
518 def gethue(phase1
, phase2
):
519 x1
= 0.5*phase1
/math
.pi
520 x2
= 0.5*phase2
/math
.pi
521 while x1
-x2
> 0.5: x1
-= 1
522 while x1
-x2
< -0.5: x1
+= 1
524 while dx
> 1: dx
-= 1
525 while dx
< 0: dx
+= 1
529 for i
in range(len(heights
)-1):
530 for j
in range(len(heights
[i
])-1):
531 midpoint
= mygrid
.point_3d(0.25 * (datagrid
[i
][j
].x_3d
+ datagrid
[i
+1][j
].x_3d
+ datagrid
[i
][j
+1].x_3d
+ datagrid
[i
+1][j
+1].x_3d
),
532 0.25 * (datagrid
[i
][j
].y_3d
+ datagrid
[i
+1][j
].y_3d
+ datagrid
[i
][j
+1].y_3d
+ datagrid
[i
+1][j
+1].y_3d
),
533 0.25 * (datagrid
[i
][j
].z_3d
+ datagrid
[i
+1][j
].z_3d
+ datagrid
[i
][j
+1].z_3d
+ datagrid
[i
+1][j
+1].z_3d
))
534 triangles
.append(triangle(datagrid
[i
][j
], datagrid
[i
+1][j
], midpoint
, gethue(datagrid
[i
][j
].phase
, datagrid
[i
+1][j
].phase
)))
535 triangles
.append(triangle(datagrid
[i
+1][j
], datagrid
[i
+1][j
+1], midpoint
, gethue(datagrid
[i
+1][j
].phase
, datagrid
[i
+1][j
+1].phase
)))
536 triangles
.append(triangle(datagrid
[i
+1][j
+1], datagrid
[i
][j
+1], midpoint
, gethue(datagrid
[i
+1][j
+1].phase
, datagrid
[i
][j
+1].phase
)))
537 triangles
.append(triangle(datagrid
[i
][j
+1], datagrid
[i
][j
], midpoint
, gethue(datagrid
[i
][j
+1].phase
, datagrid
[i
][j
].phase
)))
544 colors
= [t
[1] for t
in triangles
if t
[1] > 0]
547 for triangle
in triangles
:
549 c
.fill(triangle
[2], [color
.gray
.black
])
551 usecolor
= color
.hsb(triangle
[3], 1, 0.25 + 0.65 * ((triangle
[1] - mincol
)/(maxcol
- mincol
)) ** 0.5 )
552 c
.fill(triangle
[2], [usecolor
])
554 c
.writeEPSfile("example")