#!/usr/bin/env python

"""pert_cat.py

   Module to compute the eigenvalues for
   the quantized perturbed cat map.

   Version 0.3, 13.03.2002

   Acknowledgment: Many thanks to Fernando Perez
   who pointed out some improvements on the code.

   
   Copyright (C) 2002 Arnd Baecker <arnd.baecker@physik.uni-ulm.de>


   Usage:
   ------

    a)  # call with N and kappa given
        python pert_cat.py 51 0.3

    b)  # call without N and kappa, 
        # fallback to defaults N=101 and kappa=0.3
        python pert_cat.py 

    c)  # Use as a module (e.g. interactive session)
          
        start python and type


from AnalyseData import histogram,store_histogram
import Numeric
from math import pi
import pert_cat

N=53
N=501
kappa=0.3 
(evals,phases)=pert_cat.compute_evals_pcat(N,kappa);
# sort and unfold phases
s_phases=Numeric.sort(phases)*N/(2.0*pi)

# determine Level spacing
# (by computing the difference of the shifted eigenphases)
spacings=s_phases[1:]-s_phases[0:N]

(x_histogram,y_histogram)=histogram(spacings,0.0,10.0,100)
store_histogram(x_histogram,y_histogram,"histogram.dat")


Then use your favourite plotting programm to plot the spacing
histogram. For gnuplot you could do

goe_approx(x)=pi/2.0*x*exp(-pi/4*x*x)
gue_approx(x)=32/pi/pi*x*x*exp(-4/pi*x*x)
plot "histogram.dat" w l,goe_approx(x),gue_approx(x),exp(-x)


"""


import LinearAlgebra
from Numeric import *



def quantum_cat(N,kappa):
    """For a given N and kappa this functions returns
       the corresponding unitary matrix U 
       of the quantized perturbed cat map.
    """

    print "# Fill matrix ... N=%d  kappa=%g" % (N,kappa)
    # One could use the following lines:
    # I=1j                        # predefine sqrt(-1)
    # # First initialize complex matrix with NxN elements
    # mat=zeros((N,N), Complex)
    # for k in range(0,N):
    #     for l in range(0,N):
    #         mat[k,l]=cmath.exp(2.0*I*pi/N*(k*k-k*l+l*l)+ \
    #                        I*kappa*N/2.0/pi*sin(2.0*pi/N*l))/sqrt(N)
    #
    # However, the following is much faster (thanks to Fernando)
    alpha = 2.0*pi/N
    mat_fn = lambda k,l: alpha*(k*k-k*l+l*l)
    phi = fromfunction(mat_fn,(N,N)) + (kappa/alpha)*sin(alpha*arange(N))
    mat = exp(1j*phi)/sqrt(N)    # note that 1j=sqrt(-1)
    print "# Fill done"
    return(mat)

       
def compute_evals_pcat(N,kappa):
    """For a given N and kappa this functions 
       returns the eigenvalues and eigenphase
       of the unitary matrix U filled via quantum_cat(N,kappa).
    """


    # fill matrix U_N
    matU=quantum_cat(N,kappa)
    # determine the eigenvalues of the matrix U_N
    evals=LinearAlgebra.eigenvalues(matU)
        
    # determine phase \in [0,2\pi] from the eigenvalues
    phases_N = arctan2(evals.imag,evals.real) + pi
    # useful to determine level-spacing
    phases = concatenate( [phases_N,[phases_N[0]+2.0*pi] ] )
    
    return(evals,phases)    
           



####################################################################
# If called as main, then perform a test:
####################################################################

if __name__ == '__main__': 
    from string import atoi,atof
    import sys

    # check, if enough arguments are supplied
    if(len(sys.argv)<3):        # use default values
        N=101
        kappa=0.3
    else:                       # get values from commandline 
        N=atoi(sys.argv[1])       
        kappa=atof(sys.argv[2])   

    # Determine eigenvalues
    (evals,phases)=compute_evals_pcat(N,kappa)

    # print eigenvalues and eigenphases
    print "#   eval.real     eval.imag           phase        abs(eval)"
    for k in range(0,N):
        print("% e % e   % e   % e ") %  \
             (evals[k].real,evals[k].imag,phases[k],abs(evals[k]))