%Gaussian Wave Packet Plotter - this program allows the user to form 
%Gaussian wave packets by taking in a user-specified center wavelength
%and fractional Gaussian spectral width (sigma_p/p0) in momentum space.
%These values are used to form a wave packet phi_tilde(p), which is then
%displayed in the momentum domain, and the corresponding position
%wavefunction psi(x) is computed and displayed.

clear all;%Clears any existing data

options.WindowStyle='normal'; %Allows user to manipulate windows while ...
                                %waiting for user input
options.Interpreter='tex'; %Allows use of LaTeX notation in titles

%The following code gets user input for the center wavelength (lambda0)
%and the fractional spectral width (the ratio of the standard deviation
%sigma_p over the center momentum p0) for the momentum spectrum

prompt={'Center wavelength (\lambda_0) in meters',...%Prompts user
    'Momentum fractional spectral width (\sigma_p over p_0)'};
title1='Wavelength and Fractional Momentum Gaussian width';%Dialog box title
dims=[1,80]; %Dimensions of dialog box
definput={'1e-9','0.015'}; %Default inputs
params=inputdlg(prompt,title1,dims,definput,options); %Gets user inputs
    lambda0=str2num(char(params(1)));%for center wavelength and
    frac_p=str2num(char(params(2)));%fractional width of momentum Gaussian 
hbar=1.0546e-34;%Planck's modified constant
p0=hbar*(2*pi/lambda0);%Determines center momentum based on center wavelength
sigma_p=frac_p*p0;%Finds momentum standard deviation
sigma_x=hbar/sigma_p;%Finds width of position wavefunction

%The following code makes a vector (pvec) that runs from slightly above
%zero to twice the center momentum (p0).

pmax=2*p0;%Maximum momentum set to twice center momentum
pinc=pmax/1024;%Increment of momentum
pvec=pinc:pinc:pmax;%Makes momentum vector starting at 0+pinc and extending
                    %to pmax in steps of pinc

figure('units','normalized','outerposition',[0.1 0.1 0.8 0.8]);%Opens and
                            %sets the size & position of a figure window

A_p=(1/(sigma_p^2*pi))^0.25;%Sets amplitude of momentum Gaussian
phi_tilde=A_p*exp(-((pvec-p0).^2)/(2*(sigma_p^2)));%Makes momentum Gaussian
                                                    %wavefunction

subplot(1,3,1);%Refers to leftmost plot in group of 3 horizontal plots
plot(pvec,phi_tilde,'k');%Plots momentum wavefunction
axis([p0-3*sigma_p p0+3*sigma_p 0 1.3*max(phi_tilde)],'square');%Sets axis range
xlabel('p (Js/m)');%Labels x-axis
ylabel('$\tilde{\phi}$(p)','Interpreter','latex');%Labels y-axis
text(p0-3*sigma_p,1.2*max(phi_tilde),...
    ['\lambda_0= ',num2str(lambda0),' m','     p_0= ',num2str(p0,2),' Js/m']);
    %Writes center wavelength and center momentum on plot
text(p0-3*sigma_p,1.1*max(phi_tilde),...
    ['\sigma_p= ',num2str(sigma_p,2),' Js/m']);%Writes momentum standard
                                                %deviation on plot
title('   Momentum Spectrum');%Title of subplot
pbaspect([1 1 1]);%Makes sure plots are square (1:1 aspect ratio)
grid on;%Draws horizontal and vertical gridlines


delta_x=hbar*(2*pi/max(pvec));%Sets step size in position (x)
xwin=hbar*(2*pi/pinc);%Sets size of position window
xvec=-xwin/2:delta_x:(xwin/2);%Makes position vector
A_x=(1/(sigma_x*sqrt(pi))^0.5);%Sets amplitude of position wavefunction
psi=A_x*(exp(-xvec.^2./(2*sigma_x ^2))).*(exp(i*p0*xvec/hbar));%Makes
                                                %position wavefunction

subplot(1,3,2);%Refers to middle plot of 3 horizontal plots
plot(xvec,real(psi),'-k');%Plots real part of position wavefunction
axis([-3*sigma_x 3*sigma_x -1.2*max(real(psi)) 1.2*max(real(psi))]);
                        %Sets axis limits
xlabel('x (m)');%Labels x-axis
ylabel('Real \psi(x)');%Labels y-axis
title('Real Part of \psi(x)');%Adds title to middle subplot
grid on;%Draws horizontal and vertical gridlines
pbaspect([1 1 1]);%Makes sure plots are square (1:1 aspect ratio)

subplot(1,3,3);%Refers to rightmost plot of 3 horizontal plots
plot(xvec,abs(psi),'-k');%Plots magnitude of position wavefunction
axis([-3*sigma_x 3*sigma_x 0 1.3*max(abs(psi))]);%Sets axis limits
xlabel('x (m)');%Labels x-axis
ylabel('|\psi(x)|');%Labels y-axis
text(-3*sigma_x,1.2*max(abs(psi)),...
    ['$\sigma_x=$ ',num2str(sigma_x,2),' m'],'Interpreter','latex');
        %Writes standard deviation of position wavefunction on plot
text(-3*sigma_x,1.1*max(abs(psi)),...
 ['$\sigma_x \sigma_p=$ ',num2str(sigma_x*sigma_p,3),' $ \textrm{Js}=\hbar$'],...
 'Interpreter','latex');%Shows that product (sigma_x)(sigma_p)equals
                            %modified Planck constant

title('Magnitude of \psi(x)');%Adds title to rightmostsubplot
grid on;%Draws horizontal and vertical gridlines
pbaspect([1 1 1]);%Makes sure plots are square (1:1 aspect ratio)