%Wave packet plotter - this program takes in a user-specified wavelength
%and wavenumber fractional spectral width and uses these values to form a
%wave packet. The program then displays the wave packet as a function of
%wavenumber phi(k) and its Fourier transform, which is the position
%^wavefunction psi(x).
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 (Delta k/k) for the wavenumber spectrum

prompt={'Center Wavelength (\lambda_0) in meters)',...%Prompts user
    'Fractional spectral width \Delta k/k (0 to 1)'};
title1='Wavelength and fractional width'; %Title of dialog box
dims=[1,70]; %Dimensions of dialog box
definput={'1','0.1'}; %Default inputs
params=inputdlg(prompt,title1,dims,definput,options); %Gets user inputs
    lambda0=str2num(char(params(1)));
    kfrac=str2num(char(params(2)));
k0=2*pi/lambda0;%Determines center wavenumber based on center wavelength
deltak=k0*kfrac;%Finds width of wavenumber spectrum using fractional width

%The following code makes a vector (kvec) that runs from slightly above
%zero to twice the center wavenumber (k0).  Within that range, the
%wavenumber values will be zero between 0 and kleft, then the wavenumber
%values will be one between kleft and kright, and then the wavenumber
%values will return to zero between kright and kmax.  So the wavenumber
%spectrum is a rectangular pulse with a flat top.

kleft=k0-(deltak/2);%Finds lowest non-zero wavenumber
kright=k0+(deltak/2);%Finds highest non-zero wavenumber
kmax=2*k0;%Sets maximum wavenumber in spectrum (wavenumber spectrum
           %is the full range of wavenumbers, including wavenumbers with
           %value of zero as well as wavenumbers with value of one)
kinc=kmax/1024;%Sets wavenumber increment
kmin=kinc-(1e-6);%Sets minimum wavenumber in spectrum to slightly less
                    %than one increment
kvec=kmin:kinc:kmax;%Makes a vector of wavenumbers
phi=zeros(1,length(kvec));%Sets the value of all wavenumbers to zero
phi(floor(kleft/kinc):ceil(kright/kinc))=1;%Sets the values of the non-zero
                                            %wavenumbers to one

xinc=2*pi/kmax;%Sets the position increment
xwin=2*pi/kinc;%Sets the position window
xvec=-xwin/2+1e-12:xinc:(xwin/2);%Makes a position vector extending from
                    %the left side of the position window to the right side
                    %of the position window (the small offset of 1e-12
                    %avoids having a zero value in xvec, which appears in
                    %the denominator of the position wavefunction.

psi=(1/sqrt(2*pi))*deltak*(exp(i*k0*xvec)).*(sin(deltak*xvec/2)...
    ./(deltak*xvec/2));%Calculates the position wavefunction (I've already
          %taken the Fourier transform of phi(k), and this function
          %is the result - so the program doesn't have to calculate an FFT).

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

subplot(1,3,1);%Refers to the leftmost plot in horizontal group of 3 graphs
plot(kvec,phi,'k','LineWidth',1);%Plots the wavenumber spectrum (phi)
grid on;%Draws vertical and horizontal grid lines
plotlim_neg_x=k0-4*deltak;%Sets plot limit for negative x
plotlim_pos_x=k0+4*deltak;%Sets plot limit for positive x
plotlim_neg_y=0;%Sets plot limit for negative y
plotlim_pos_y=1.2*max(phi);%Sets plot limit for positive y
axis([plotlim_neg_x plotlim_pos_x plotlim_neg_y plotlim_pos_y]...
,'square');%Sets axis range
xlabel('k (rad/m)');%Labels x-axis4

ylabel('\phi(k)');%Labels y-axis
text(plotlim_neg_x,0.9*plotlim_pos_y,...
    ['\lambda_0= ',num2str(lambda0),' m','   k_0= ',num2str(k0,2),' rad/m']);...
    %Writes center wavelength and wavenumber on plot
text(plotlim_neg_x,0.8*plotlim_pos_y,...
    ['\Deltak/k= ',num2str(kfrac)]);%Writes fractional spectral width on plot
title('Wavenumber Spectrum \phi(k)');
pbaspect([1 1 1]);%Makes sure plots are square (1:1 aspect ratio)
hold on;%Allows multiple plots on same graph

subplot(1,3,2);%Refers to center plot in horizontal group of three graphs
plot(xvec,real(psi),'-k','LineWidth',1);%Plots real part of position wavefunction
grid on;%Draws vertical and horizontal grid lines
axis([-4*pi/deltak 4.001*pi/deltak -1.2*max(real(psi)) 1.2*max(real(psi))]);
xlabel('x (m)');
ylabel('Real \psi(x)');
title('Real Part of \psi(x)');
pbaspect([1 1 1]);%Makes sure plots are square (1:1 aspect ratio)

subplot(1,3,3);%Refers to rightmost plot in horizontal group of three graphs
plot(xvec,abs(psi),'-k','LineWidth',1);%Plots magnitude of position wavefunction
grid on;%Draws vertical and horizontal grid lines
axis([-4*pi/deltak 4.001*pi/deltak 0 1.2*max(abs(psi))]);
xlabel('x (m)');
ylabel('|\psi(x)|');
text(-4*pi/deltak,1.1*max(abs(psi)),...
    ['Full width at 1/e points \approx ',num2str(1.4*lambda0/kfrac),' m']);
    %Writes width of position wavefunction on graph
title('Magnitude of \psi(x)');
pbaspect([1 1 1]);%Makes sure plots are square (1:1 aspect ratio)