[Matlab] FFT Gaussian x-space -> p-space scalings.

[apoptosa]

Hi There,
Im having a little trouble with scalings required for a fft performed in Matlab. So here is the story:
I construct a position space gaussian, plot this guy out to 6 stdDev. My frontfactor is such that my probability density is normalized.
I know what the actual Fourier transform should be in momentum space, thus when I use fft in MATLAB I should get something that lies ontop of this solution.

My problem is that I can't, also I cannot really see where things are going wrong. So here is some code to peruse if you're curious about helping,
Thanks in advance :D


%A gaussian wavepacket
asimParam %physical constants [SI units]

%experimental parameters
k0=0; %planewave with this wavenumber is modulated by a Gaussian
w=100e-6; %Gaussian width parameter

%simulation parameters
stdDevNum=6; %number of standard deviations before truncation
disp(['Proportion of total Gaussian within (' num2str(stdDevNum) ') std dev. :'...
num2str(erf(stdDevNum/sqrt(2)))]);

%calculate the maximum momentum to use before truncating the Gaussian, this
%is done in terms of the maximum number of std. dev. to consider. This can
%be done because an explicit form of the std dev. is known (See Sakurai)
stdDevX=w/sqrt(2);
xMax=stdDevNum*stdDevX;
grX=1024; %make sure this is base 2pwr as this optimises fft

x=linspace(-xMax,xMax,grX)';
dx=x(2,1)-x(1,1);

%calculate the x-space wavefunction |psi(x,t=0)>
xPsi=1/(pi^(1/4)*sqrt(w))*exp(i*k0*x-x.^2/(2*w^2));

rhoPsiX=conj(xPsi).*xPsi; %rho as a func of x
rhoPsiXint=sum(dx*rhoPsiX); %sum of density (step size uniform);


%Momentum space wavefunction: |phi(p,t=0)>
%First use the explicit form given in Sakurai, then compare the fft
%function in Matlab to ensure things are working correctly.
pPhiAn=@(pIn)sqrt(w/(hb*sqrt(pi)))*exp(-(pIn-hb*k0).^2*w.^2/(2*hb^2));
stdDevP=hb/(sqrt(2)*w);
N=grX; %same number of points as for x-space
dk=1/(N*dx); %such that Fourier transform has correct units
dp=hb*dk; %now rescale to momentum

%pMax=stdDevP*stdDevNum;
pMax=dp*N/2;
pIn=linspace(-(pMax+hb*k0), (pMax+hb*k0), N)';
pPhiS=pPhiAn(pIn);
rhoPhiPs=conj(pPhiS).*pPhiS;
rhoPhiPsInt=sum(dp*rhoPhiPs); %sum of density (step size uniform);

%Now using Matlab's fft
pPhi=fftshift(fft(xPsi));
rhoPhiP=conj(pPhi).*pPhi;

%test if it works back to x-space
xPsiF=ifft(fft(xPsi));
rhoPsiXF=conj(xPsiF).*xPsiF;

%The rest of this m-file is comprised of plotting

 

[NateD]

the wavepackets and%obtaining various statistical properties of the wavepackets.%Plotting x-spacefigure(1)clf;plot(x,real(xPsi),'r',x,imag(xPsi),'b')title('x-space: Wavefunction real(r) and Imaginary(b)');xlabel('x [m]');ylabel('\psi(x) [a.u]');legend('Re','Im');grid on%plotting p-spacefigure(2)clf;plot(pIn,real(pPhiAn(pIn)),'r',pIn,imag(pPhiAn(pIn)),'b')title('p-space: Wavefunction real(r) and Imaginary(b)');xlabel('p [kg m/s]');ylabel('\phi(p) [a.u]');legend('Re','Im');grid on%plotting x-space densityfigure(3)clf;plot(x,rhoPsiX)title('x-space: \rho(x)');xlabel('x [m]');ylabel('\rho(x) [a.u]');grid on%plotting p-space densityfigure(4)clf;plot(pIn,rhoPhiP)title('p-space: \rho(p)');xlabel('p [kg m/s]');ylabel('\rho(p) [a.u]');grid on</code>A:There seem to be two issues with the code you provided. One is the scaling error, and the other is the incorrect calculation of $\phi$.The scaling error can be fixed simply by multiplying $p\phi$ by a normalization factor:<blockquote><p><span class="math-container">$$p\phi \rightarrow \frac{p}{