%% Visualisierung der Wignerfunktion - Tests

% Vers. 2.8; 07.2015 
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

%%

clear,close all, clc
%clear
aufl = 401;% 201 genuegt 401 fuer Abb. Buch
x = linspace(-10.5,10.5,aufl);
y = linspace(-10.5/2,10.5/2,aufl);
[X,Y] = meshgrid(x,y);
%%
n=3;%0;%4;%6
p = hermitpoly(n);
fpoly = @(x) polyval(p,x);
%%
tic
Z = 1/(sqrt(pi)*2^n*factorial(n)).*arrayfun(fpoly,X-Y).*exp(-1/2*(X-Y).^2).*exp(-1/2*(X+Y).^2).*arrayfun(fpoly,X+Y);
toc
%%
posc = linspace(-10.5,10.5,aufl);
[X,P] = meshgrid(x,posc);
%%
r = sqrt(2*n+1);
phi = 0:0.1:2*pi;
xcl = r*sin(phi);
pcl = r*cos(phi);
%%
tic
zwi1 = exp(2*i*posc(1)*Y);
Wqp = sum(Z.*zwi1);
toc
tic
for k=2:length(posc)
    zwi1 = exp(2*i*posc(k)*Y);
    Wqp = [Wqp;sum(Z.*zwi1)];
end
fh = figure;
surf(X,P,abs(Wqp)), shading interp, shg    % Wignerfunktion - Absolutwert
hold on
zzz = zeros(size(xcl));
plot3(xcl,pcl,zzz)
%view(37.5,90)
hold off
figure,contour(X,P,abs(Wqp))               % Wignerfunktion Konturplot
hold on

plot(xcl,pcl,'k');
hold off
toc
%%
% Abb fuer Buch:
aview = -44.5000;
bview = 80;
view(aview,bview);


WqpPlot = abs(Wqp);
mitte=round(aufl/2);
maxa=mean(mean(WqpPlot(mitte-3:mitte+3,mitte-3:mitte+3)))
WqpPlot(WqpPlot > maxa) = maxa;
% Exakte Spitze und Umgebung aus Visualisierungsgruenden gedaempft
% Achtung Fouriertransformationen: "0-Punkt" enth?lt numerischen
% Summenanteil: --> exp(i*0*...) 
%maxa=mean(mean(WqpPlot(mitte-10:mitte+10,mitte-10:mitte+10)))
%WqpPlot(WqpPlot > maxa) = maxa;
ml=mitte-90;
mr=mitte+90;
WqpPlot=WqpPlot(ml:mr,ml:mr);
Xp = X(ml:mr,ml:mr);
Pp = P(ml:mr,ml:mr);
WqpPlot(1:floor(length(WqpPlot)/2),:)=NaN;
figure, surf(Xp,Pp,WqpPlot,'FaceColor','interp','EdgeColor','none','FaceLighting','phong')
hold on
%colormap gray
%load CM
%colormap(CM)
view(aview,bview);
plot(xcl,pcl,'k');
axis tight
xlabel('q')
ylabel('p')
set(gca,'ZTickLabel',[]);
%title(['n = ',num2str(n)])
hold off
%%
% figure, surf(Z), shading interpxlabel('q')
% figure, contour(Z),shg
% Zfft = fft(Z,[],1);
% %figure,surf(abs(Zfft))
% Zpl = real(fftshift(Zfft,1));
% figure,surf(Zpl),shading interp%,colormap(gray)
% %axis equal
% figure, contour(sqrt(Zpl)),shg
% figure(1)