# Fourier transform via Hamilton from Hermite functions
rmax = 128;
nH= rmax;
# Gruenbaum's tridiagonal Hamilton operator
# a Hamiltonian is just an operator commuting with Fourier transform
# having distinct eigenvalues
sm = sin(pi*(0:nH)/nH) ;
am = 2 * (shift(sm,-1) .* fliplr(sm)) (1:nH);
bm = (sm .* fliplr(sm))(2:(nH)) ;
T = diag(bm,-1) + diag(am) + diag(bm,1) ;
# we are interested in the eigenvectors only
[Vx,~] = eig(T) ;
#plot(Vx*diag(exp((-nM:nM)*0.5i*pi))*Vx') ; # -1i*S*F
#plot((Vx*(eye(nH)+diag(0.01i*pi*(0:nH-1)/nH))*Vx')^(50*nH)) ;
# norm(Vx*diag(exp(((-nM):nM)*0.5i*pi))*Vx'-1i*F')
# the real Hamilton operator, we force the eigenvalues we need
Ht = (Vx*diag((nH-1:-1:0))*Vx') ;
Ht(127,1) = 0; Ht(1,127) = 0; Ht(128,1) = 0; Ht(1,128) = 0;
Hts = sparse(Ht.*(abs(Ht)>0.1));
# this Hamiltonian is not minimal sparse, more exploration is needed
# make the Hamiltonian for 2-dimensional domain
Hts2= kron(Hts,speye(rmax)) + kron(speye(rmax),Hts);

# now we verify, that the solution of Schroedinger equation is the Ft
N =5000; % reicht auch eine kleinere Zahl?
H0x = 0.5i/N*pi*Hts2;
# T2: Korrektur für bessere Orthogonalität
T2 = (speye(rmax^2)-H0x)*(speye(rmax^2)-H0x')*(speye(rmax^2)-H0x) - (speye(rmax^2)-H0x);
H0x += 0.5*T2 .* spones(H0x);
T2 = diag(sumsq(speye(rmax^2)-H0x))-speye(rmax^2);
H0x += 0.5*T2;
if exist("Ifovf")
  D=conj(vec(Ifovf(257+(-64:63),257+(-64:63)))/4);
  for i=1:N; D -= H0x*D;endfor
  Df = D;
  for i=1:N; D -= H0x*D;endfor
else # just the backtransform
  D=vec(I1f/4);
  for i=1:N; D -= H0x*D;endfor
endif
nshow(reshape(D,size(I1f)))

#[~,~,Q] = laplacian([128,128]);
#qq = spones(Q*Q);

#Q0 = speye(rmax^2) - H0x;
#Q1 = Q0;
#Q1 = 0.5 *Q1 .* (abs(Q1)>1/100000); Q1 += (Q0 * Q1');
#Q0 = Q1 * Q1.';
#Q0 = Q0 .* (abs(Q0)>1/100000);
#Q0 = speye(rmax^2) - Q0

if(0)
# trying to make the Hamilton sparser
H1 = 0.5*H0x;
H1.*=(abs(H1)>0.0002);
H1 = 0.5*(H1+(Hx+H0x)*(H1+Hx)'-Hx).*spones(H1);
D=conj(vec(Ifovf(257+(-64:63),257+(-64:63)))/4);
for i=1:N; D -= H1*D;endfor
for i=1:N; D -= H1*D;endfor
nshow(reshape(D,size(I1f)))
endif
