%Infinite potential well probability density plotter. This program allows %the user to add two or more infinite-well position wavefunctions together %and to see how the combined wavefunction changes over time. clear all;%Clears any existing data options.WindowStyle='normal'; %Allows user to manipulate windows while ... %waiting for user input prompt={'Number of wavefunctions to add','Width of Well (meters)','Mass (kg)'}; %Prompts user for input title1='Infinite Potential Well';%Title of dialog box dims=[1,60];%Dimensions of dialog box definput={'2','1e-10','9.11e-31'};%Default inputs params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs for numfuncs=str2num(char(params(1))); %number of wavefunctions to add a=str2num(char(params(2))); %width of potential well in meters m=str2num(char(params(3))); %mass of particle in kg for count=1:numfuncs;%Loop to get quantum number and amount of each function prompt={['Quantum Number (n) of wavefunction #',num2str(count)],... ['Relative amount of wavefunction #',num2str(count),' (0 to 1)']}; %Prompts user for input title1='Wavefunction Numbers and Amounts';%Title of dialog box dims=[1,60];%Dimensions of dialog box definput={num2str(count),'0.5'};%Default inputs params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs for n(count)=str2num(char(params(1))); %quantum number of wavefunction amt(count)=str2num(char(params(2))); %amount of wavefunction end; total_amt=sum(amt);%Finds total of amounts for all wavefunctions frac=sqrt(amt/total_amt);%Finds fraction of total for each wavefunction hbar=1.05457e-34;%Planck's modified constant xinc=0.01*a;%Position (x) increment xvec=0:xinc:a+xinc;%Makes vector of position values xvec_norm=xvec/a;%Normalizes position vector by dividing by width of well kn=n*pi/a;%Finds wavenumber for each wavefunction En=hbar^2*kn.^2/(2*m);%Finds energy for each wavefunction %The following statements find the period of the composite wavefunction. %If only two wavefunctions are added, the period is related to the %energy difference between those two wavefunctions. If three or more %wavefunctions are added, the period depends on the least common multiplier %(lcm) of the inverse square of the quantum numbers (n^-2). if numfuncs==2;%Handles the case of only two wavefunctions being added omega=abs(En(2)-En(1))/hbar;%Finds angular frequency (omega) based on the %energy difference between the wavefunctions T=2*pi/omega;%Finds the period based on the angular frequency (omega) Sampletime=T/200;%Sets the time between samples as 1/200th of the period else;%Handles the case of more than two wavefunctions being added factor_lcm=double(lcm(sym(n.^-2)));%This statement finds a factor needed %for the period when 3 or more wavefunctions are added. %To find the least common multiple of the inverse %square of the quantum numbers, the MATLAB lcm function %requires n^-2 to be converted into a symbolic expression. %Once the lcm function operates has found the least %common multiple of n^-2, that values is converted into %a double-precision numerical value and stored as the %variable factor_lcm. omega=min(En)/hbar;%Finds the smallest angular frequency based on the %lowest energy T=2*pi*factor_lcm*min(n.^2)/omega;%Uses the least common multiple factor %and the smallest angular frequency %to find the period omega_max=max(En)/hbar;%Finds the biggest angular frequency of the %wavefunctions Min_period=2*pi/omega_max;%Find the shortest period of the wavefunctions Sampletime=Min_period/20;%Sets the time between samples as 1/20th of the %shortest period end; figure('units','normalized','outerposition',[0.25 0.42 0.5 0.5]);%Opens and %sizes figure window for t=(0:Sampletime:T);%Makes time variable (t) for count=1:numfuncs;%Counts through the number of wavefunctions %selected by the user psi_x(count,:)=frac(count)*sqrt(2/(a))*sin(kn(count)*xvec)... *exp(-i*En(count)*t/hbar);%Makes the position wavefunctions %as 2-dimensional matrices. The first %index specifies the wavefunction number %and the second index specifies the time, %from t=0 to t=T with step size equal to %Sampletime. end; psi_xtot=sum(psi_x);%Adds the wavefunctions together. This sum %operates on the first index, so it adds %wavefunction 1 to wavefunction 2 (and %wavefunctions 3, 4, etc. if the user has %specified more than 2 wavefunctions). x_probden_norm=a*(conj(psi_xtot).*psi_xtot);%The combined position %wavefunction (psi_xtot) is complex, so the %position probability density is formed by %multiplying the wavefunction by its complex %conjugate. It is normalized by multiplying %by the width of the potential well (a). plot(xvec_norm,x_probden_norm,'k');%Plots the normalized probabilty %density as a function of the normalized %position. axis([0 max(xvec_norm) 0 2*numfuncs]);%Sets the limits of the axes grid on;%Draws vertical and horizontal grid lines xlabel('Normalized Position (x/a)');%Labels the x-axis ylabel('Normalized Probability Density');%Labels the y-axis title('Infinite Potential Well Probability Density vs. Time'); %The following statements add text information to the plots. The %first "text" statement and the "for" loop write the combined %wavefunction as "Psi(x,t)=" followed by two or more wavefunctions %(fraction 1) psi_1 plus (fraction 2)psi_2 + ... (up to however %many wavefunctions the user specified). text(0.01,1.9*numfuncs,... ['\Psi(x,t) = ',num2str(frac(1),'%4.2f'),'\psi_{' int2str(n(1)) '}']); for count=2:numfuncs; text(0.1*count,1.9*numfuncs,... ['+',num2str(frac(count),'%4.2f'),'\psi_{' int2str(n(count)) '}']); end; %The following two "text" statements write the time and the period %on the plot as the time increments upward from 0 to T. text(0.01,1.75*numfuncs,['t= ',num2str(t,2),' sec']); text(0.6,1.9*numfuncs,['Period= ',num2str(T,2),' sec']); %The following three statements make three thick lines to show outline of %potential well on the position-density plot line([0 0],[0 2*numfuncs],'Linewidth',3,'Color', 'k'); line([1 1],[0 2*numfuncs],'Linewidth',3,'Color', 'k'); line([0 1],[0 0],'Linewidth',3,'Color', 'k'); pause(0.02); end;