www/png/example_mkdata.py

   1 # this code is really a mess ...
   2
   3 import sys
   4 from Numeric import *
   5 from RandomArray import *
   6 from LinearAlgebra import *
   7 from FFT import *
   8 from functools import reduce
   9
  10
  11 sqrt2 = math.sqrt(2)
  12
  13
  14 def norm(v):
  15     return math.sqrt(reduce(lambda x, y: x+y*y, v, 0))
  16
  17 def blochvector(l, k):
  18     assert k >= 0 and k <= l/2
  19     v = cos(2*math.pi/l*(k)*arange(l))
  20     return v/norm(v)
  21
  22 def localvector(l, pos):
  23     v = zeros(l)
  24     v[pos] = 1
  25     return v/norm(v)
  26
  27 def local2vector(l, pos1, pos2):
  28     v = zeros(l)
  29     v[pos1] = 1
  30     v[pos2] = 1
  31     return v/norm(v)
  32
  33
  34
  35 def diag_anderson(l, w):
  36
  37     """diagonalize the anderson model"""
  38
  39     seed(x=1705, y=1111)
  40     pot=w*random(l)-0.5*w
  41     ham = zeros((l, l), Float)
  42     for i in range(l):
  43         ham[i, i] = pot[i]
  44         ham[i, (i+1)%l] = -1
  45         ham[i, (i-1)%l] = -1
  46     evalues, evectors = eigenvectors(ham)
  47     sorter = argsort(evalues)
  48     return take(evalues, sorter), take(evectors, sorter)
  49
  50
  51
  52 def diag_harper_full(l, lam):
  53
  54     """diagonalize the harper model by full diagonalization"""
  55
  56     def fibs(tonumber):
  57         result = [1, 1]
  58         while result[-1] < tonumber:
  59             result.append(result[-2]+result[-1])
  60         return result
  61
  62     l2, ltest = fibs(l)[-2:]
  63     if ltest != l: raise ValueError("l isn't a fibonacci number")
  64
  65     pot=lam*cos((2*pi*l2*arange(l))/l)
  66     ham = zeros((l, l), Float)
  67     for i in range(l):
  68         ham[i, i] = pot[i]
  69         ham[i, (i+1)%l] = 1
  70         ham[i, (i-1)%l] = 1
  71     evalues, evectors = eigenvectors(ham)
  72     sorter = argsort(evalues)
  73     return take(evalues, sorter), take(evectors, sorter)
  74
  75
  76
  77 def diag_harper(l, lam):
  78
  79     """diagonalize the harper model by diagonalizing the symmetric and antisymmetic parts"""
  80
  81     def fibs(tonumber):
  82         result = [1, 1]
  83         while result[-1] < tonumber:
  84             result.append(result[-2]+result[-1])
  85         return result
  86
  87     def createsymm(dest, source):
  88         l = len(dest)
  89         if l % 2:
  90             dest[0] = source[0]
  91             for i in range(1, (l+1)/2):
  92                 dest[i] = dest[l-i] = source[i] / sqrt2
  93         else:
  94             dest[0] = source[0]
  95             dest[l/2] = source[l/2]
  96             for i in range(1, l/2):
  97                 dest[i] = dest[l-i] = source[i] / sqrt2
  98
  99     def createasymm(dest, source):
 100         l = len(dest)
 101         if l % 2:
 102             dest[0] = 0
 103             for i in range(1, (l+1)/2):
 104                 dest[i] = source[i-1] / sqrt2
 105                 dest[l-i] = -source[i-1] / sqrt2
 106         else:
 107             dest[0] = dest[l/2] = 0
 108             for i in range(1, l/2):
 109                 dest[i] = -source[i-1] / sqrt2
 110                 dest[l-i] = source[i-1] / sqrt2
 111
 112     l2, ltest = fibs(l)[-2:]
 113     if ltest != l: raise ValueError("l isn't a fibonacci number")
 114
 115     # symmetric and antisymmetric diagonalization
 116     if l%2:
 117         # odd
 118         symmsize = (l-1)/2 + 1
 119         pot = lam*cos((2*pi*l2*arange(symmsize))/l)
 120         pot[symmsize - 1] += 1
 121         symmham = zeros((symmsize, symmsize), Float)
 122         for i in range(symmsize):
 123             symmham[i, i] = pot[i]
 124             if i > 1:
 125                 symmham[i, i-1] = symmham[i-1, i] = 1
 126             elif i:
 127                 symmham[i, i-1] = symmham[i-1, i] = sqrt2
 128     else:
 129         # even
 130         symmsize = l/2 + 1
 131         pot = lam*cos((2*pi*l2*arange(symmsize))/l)
 132         symmham = zeros((symmsize, symmsize), Float)
 133         for i in range(symmsize):
 134             symmham[i, i] = pot[i]
 135             if i > 1 and i < symmsize - 1:
 136                 symmham[i, i-1] = symmham[i-1, i] = 1
 137             elif i:
 138                 symmham[i, i-1] = symmham[i-1, i] = sqrt2
 139     symmevalues, symmevectors = eigenvectors(symmham)
 140
 141     if l%2:
 142         # odd
 143         asymmsize = (l-1)/2
 144         pot = lam*cos((2*pi*l2*(arange(asymmsize)+1))/l)
 145         pot[asymmsize-1] -= 1
 146         asymmham = zeros((asymmsize, asymmsize), Float)
 147         for i in range(asymmsize):
 148             asymmham[i, i] = pot[i]
 149             if i:
 150                 asymmham[i, i-1] = asymmham[i-1, i] = 1
 151     else:
 152         # even
 153         asymmsize = l/2 - 1
 154         pot = lam*cos((2*pi*l2*(arange(asymmsize)+1))/l)
 155         asymmham = zeros((asymmsize, asymmsize), Float)
 156         for i in range(asymmsize):
 157             asymmham[i, i] = pot[i]
 158             if i:
 159                 asymmham[i, i-1] = asymmham[i-1, i] = 1
 160     asymmevalues, asymmevectors = eigenvectors(asymmham)
 161
 162     # build complete solution
 163     symmsort = argsort(symmevalues)
 164     asymmsort = argsort(asymmevalues)
 165     j = k = 0
 166     evalues = zeros((l,), Float)
 167     evectors = zeros((l, l), Float)
 168     for i in range(l):
 169         if j != symmsize and (k == asymmsize or
 170                               symmevalues[symmsort[j]] < asymmevalues[asymmsort[k]]):
 171             evalues[i] = symmevalues[symmsort[j]]
 172             createsymm(evectors[i], symmevectors[symmsort[j]])
 173             j += 1
 174         else:
 175             evalues[i] = asymmevalues[asymmsort[k]]
 176             createasymm(evectors[i], asymmevectors[asymmsort[k]])
 177             k += 1
 178     return evalues, evectors
 179
 180
 181
 182 def test_harper(l, lam):
 183
 184     """compare diag_harper and diag_harper_full"""
 185
 186     res1 = diag_harper(l, lam)
 187     res2 = diag_harper_full(l, lam)
 188     for x1, x2 in zip(res1[0], res2[0]):
 189         assert math.fabs(x1-x2) < 1e-8, "eigenvalues differ"
 190     for v1, v2 in zip(res1[1], res2[1]):
 191         sum = 0
 192         for x1, x2 in zip(v1, v2):
 193             sum += x1 * x2
 194         if sum > 0:
 195             for x1, x2 in zip(v1, v2):
 196                 assert math.fabs(x1-x2) < 1e-8, "eigenvectors differ\n%s\n%s" % (v1, v2)
 197         else:
 198             for x1, x2 in zip(v1, v2):
 199                 assert math.fabs(x1+x2) < 1e-8, "eigenvectors differ\n%s\n%s" % (v1, -v2)
 200
 201 # test_harper(13, 5)
 202
 203
 204
 205 def wigner_slow(vector):
 206
 207     """create wigner function of a state (direct way)"""
 208
 209     l = len(vector)
 210     wig = zeros((l, l), Float)
 211
 212     # wigner function (direct way)
 213     for x0loop in range(l):
 214         x0 = x0loop
 215         for k0loop in range(l):
 216             if l%2:
 217                 k0 = (k0loop-0.5*l+0.5)*2*pi/l
 218             else:
 219                 k0 = (k0loop-0.5*l+1)*2*pi/l
 220             sum = 0j
 221             for yloop in range(l):
 222                 y = yloop - l/2 + 1 - l%2
 223                 v = lambda x: vector[(x+10*l)%l]
 224                 if y%2:
 225                     v1 = 0.5*(v(x0-(y-1)/2) + v(x0-(y+1)/2))
 226                     v2 = 0.5*(v(x0+(y-1)/2) + v(x0+(y+1)/2))
 227                 #    v1 = 0
 228                 #    v2 = 0
 229                 else:
 230                     v1 = v(x0-y/2)
 231                     v2 = v(x0+y/2)
 232                 sum += v1 * v2 * exp(1j*k0*y)
 233             assert abs(sum.imag) < 1e-10, "imaginary part should be zero"
 234             wig[x0loop, k0loop] = sum.real
 235
 236     return wig
 237
 238
 239
 240 def wigner(vector):
 241
 242     """create wigner function of a state (fft)"""
 243
 244     l = len(vector)
 245     wig = zeros((l, l), Float)
 246
 247     fftarray = zeros(l, Float)
 248     fftresult = zeros(l, Complex)
 249     wig = zeros((l, l), Float)
 250     for x0loop in range(l):
 251         x0 = x0loop
 252         for yloop in range(l):
 253             y = yloop - l/2 + 1 - l%2
 254             v = lambda x: vector[(x+10*l)%l]
 255             if y%2:
 256                 v1 = 0.5*(v(x0-(y-1)/2) + v(x0-(y+1)/2))
 257                 v2 = 0.5*(v(x0+(y-1)/2) + v(x0+(y+1)/2))
 258             #    v1 = 0
 259             #    v2 = 0
 260             else:
 261                 v1 = v(x0-y/2)
 262                 v2 = v(x0+y/2)
 263             fftarray[(y+10*l)%l] = v1 * v2
 264         fftresult = real_fft(fftarray)
 265         for k0loop in range(l):
 266             if l%2:
 267                 index = int(abs(k0loop-0.5*l+0.5))
 268             else:
 269                 index = int(abs(k0loop-0.5*l+1))
 270             wig[x0loop, k0loop] = fftresult[index].real
 271
 272     return wig
 273
 274
 275
 276 def test_wigner(vector):
 277
 278     """compare wigner_slow and wigner"""
 279
 280     res1 = wigner_slow(vector)
 281     res2 = wigner(vector)
 282     for v1, v2 in zip(res1, res2):
 283         for x1, x2 in zip(v1, v2):
 284             assert math.fabs(x1-x2) < 1e-8, "wigner function differs\n%s\n%s" % (res1, res2)
 285
 286 # test_wigner(diag_anderson(10, 1)[1][5])
 287 # test_wigner(diag_anderson(11, 1)[1][5])
 288
 289
 290 def husimi_from_wigner(wig):
 291
 292     """calculate the husimi function out of the wigner function"""
 293
 294     l = len(wig)
 295     sigma = sqrt(l/pi/4)
 296
 297     hus = zeros((l, l), Float)
 298     for x1 in range(l):
 299         for k1 in range(l):
 300             sys.stderr.write("wigner -> husimi (very slow) ...%6.2f%%\r" % ((100.0*(x1*l+k1))/l/l))
 301             hus[x1, k1] = 0
 302             for x2 in range(l):
 303                 for k2 in range(l):
 304                     dx = x1-x2
 305                     if dx < -l/2: dx += l
 306                     if dx > l/2: dx -= l
 307                     dk = k1-k2
 308                     if dk < -l/2: dk += l
 309                     if dk > l/2: dk -= l
 310                     dk *= 2*pi/l
 311                     hus[x1, k1] += 2.0/l/l * wig[x2, k2] * exp(-0.5*dx*dx/sigma/sigma-2*sigma*sigma*dk*dk)
 312     sys.stderr.write("wigner -> husimi (very slow) ... done.   \n")
 313
 314     return hus
 315
 316
 317
 318 def husimi_slow(vector):
 319
 320     l = len(vector)
 321     hus = zeros((l, l), Float)
 322     phases = zeros((l, l), Float)
 323     sigma=sqrt(l/pi/4)
 324     factor=1/sqrt(sqrt(2*pi)*sigma*l)
 325
 326     for x0loop in range(l):
 327         x0 = x0loop
 328         for k0loop in range(l):
 329             if l%2:
 330                 k0 = (k0loop-0.5*l+0.5)*2*pi/l
 331             else:
 332                 k0 = (k0loop-0.5*l+1)*2*pi/l
 333             sum = 0j
 334             for x in range(l):
 335                 phase = exp(1j*k0*x)
 336                 sum += vector[x] * factor * exp(-(x-x0-l)*(x-x0-l)/(4*sigma*sigma)) * phase
 337                 sum += vector[x] * factor * exp(-(x-x0  )*(x-x0  )/(4*sigma*sigma)) * phase
 338                 sum += vector[x] * factor * exp(-(x-x0+l)*(x-x0+l)/(4*sigma*sigma)) * phase
 339             hus[x0loop, k0loop] = abs(sum)*abs(sum)
 340             phases[x0loop, k0loop] = math.atan2(sum.imag, sum.real)
 341
 342     return hus, phases
 343
 344
 345
 346 def husimi(vector):
 347
 348     l = len(vector)
 349     hus = zeros((l, l), Float)
 350     sigma=sqrt(l/pi/4)
 351     factor=1/sqrt(sqrt(2*pi)*sigma*l)
 352
 353     fftarray = zeros(l, Complex)
 354     fftresult = zeros(l, Complex)
 355     heights = zeros((l, l), Float)
 356     phases = zeros((l, l), Float)
 357     for x0loop in range(l):
 358         x0 = x0loop
 359         for xloop in range(l):
 360             x = xloop
 361             while (x-x0 < -0.5*l): x += l
 362             while (x-x0 > 0.5*l): x -=l
 363             fftarray[xloop] = vector[xloop] * factor * exp(-(x-x0)*(x-x0)/(4*sigma*sigma))
 364         #fftresult = real_fft(fftarray)
 365         fftresult = fft(fftarray)
 366         for k0loop in range(l):
 367             if l%2:
 368                 index = (int(k0loop-0.5*l+0.5) + 10*l)%l
 369                 #index = int(abs(k0loop-0.5*l+0.5))
 370             else:
 371                 raise
 372                 index = int(abs(k0loop-0.5*l+1))
 373             heights[x0loop, k0loop] = abs(fftresult[index])*abs(fftresult[index])
 374             phases[x0loop, k0loop] = math.atan2(fftresult[index].imag, fftresult[index].real)
 375             # if 2*x0loop > l: phases[x0loop, k0loop] = -phases[x0loop, k0loop] + math.pi
 376             # if 2*k0loop > l: phases[x0loop, k0loop] = -phases[x0loop, k0loop]
 377
 378     return heights, phases
 379
 380
 381
 382 def test_husimi(vector):
 383
 384     """compare husimi_slow and husimi"""
 385
 386     res1 = husimi_slow(vector)
 387     res2 = husimi(vector)
 388     for v1, v2 in zip(res1, res2):
 389         for x1, x2 in zip(v1, v2):
 390             assert math.fabs(x1-x2) < 1e-8, "husimi function differs\n%s\n%s" % (res1, res2)
 391
 392 # test_husimi(diag_anderson(10, 1)[1][5])
 393 # test_husimi(diag_anderson(11, 1)[1][5])
 394
 395 l=101
 396 values, vectors = diag_anderson(l, 0.5)
 397 vector = vectors[24]
 398 heights, phases = husimi(vector)
 399
 400 f = open("example.dat", "w")
 401 for x in range(l):
 402     for y in range(l):
 403         f.write("%i %i %+20.15e %+20.15e\n" % (x, y, heights[x, y], phases[x, y]))