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