% function [k,W]=Wigner(y,x)   %------------------------------------------
 function W = Wigner(y,x)      %------------------------------------------

%calculates the Wigner distribution of a general 1D function.
%Written by: Nadav Katz, nadav.katz@weizmann.ac.il +972-8-9343543
%Last update: 16.10.02 
%please contact with any corrections or comments.
%syntax: [W,p]=Wigner(y,x)
%input
%y - 1D wavefunction
%x - 1D x dimension for y(x)
%IMPORTANT: assumes y to be of even length Ni, and the zero point is taken at the Ni/2+1 point
%W - Wigner distribution of y(x), i.e. Int(y(x-s/2)*conj(y(x+s/2))exp(i*p*s), -Infinity, Infinity)
%k - vector of resulting fourier "k-space" dimension, corresponds to the y-axis (row number) of W

%measure length of wavefunction
Ni=length(y);
%keyboard;
%create zero padding vector
u=zeros(1,Ni);

%pad with zeros
y=[u y u];

%conjugate variable
yc=conj(y);

%allocate Ni*Ni matrices, which will contain the indices we need from y
A=zeros(Ni);
B=zeros(Ni);

%calculate indices, note shift in i, due to padding. j is not shifted, since it is related to s
%from the above definition of the wigner distrib.
for i=Ni+1:2*Ni,
    for j=-Ni/2:Ni/2-1,
        A(j+1+Ni/2,i-Ni)=i-j;
        B(j+1+Ni/2,i-Ni)=i+j;
    end;
end;

%use matrices to calculate values of matrices for wigner
Ay=y(A);
By=yc(B);
W1=Ay.*By;

%create fourier vector
dx=x(2)-x(1);
dk=pi/dx/Ni;

k=Ni/2*dk:-dk:-(Ni/2-1)*dk;  %----------------------------------at the original pro. this was returned.

%calculate fourier transform along columns (with required fftshifts to place zero at center
%of the matrix, as far as the fft is concerned. The real value is taken since the Wigner is
%known to be real, and divided by pi/dx in order properly normalize
W=real(fftshift(ifft(fftshift(W1,1)),1))/pi*dx;
return;