%% Visualisierung der Wignerfunktion 

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

%%
clear,close all, clc
aufl = 401;       % Legt die Bildauflösung fest 
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); % Berechnet die Hermite-Polynome
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);
%%  klassische Bahn
r = sqrt(2*n+1);
phi = 0:0.1:2*pi;
xcl = r*sin(phi);
pcl = r*cos(phi);
%% Berechnung der Wigner-Funktion
zwi1 = exp(2*i*posc(1)*Y);
Wqp = sum(Z.*zwi1);
for k=2:length(posc)
    zwi1 = exp(2*i*posc(k)*Y);
    Wqp = [Wqp;sum(Z.*zwi1)];
end
%% Visualisierung
aview = -44.5000;
bview = 80;
view(aview,bview);   % Orientierung der Abb.

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" enthaelt numerischen
% Summenanteil: --> exp(i*0*...) 
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;  % nicht plotten
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