% three_D_oscillator.m - generates scatter plot of p_z orbital for three-dimensional oscillator
clear all; close all;

Phi=@(phi) 1/sqrt(2*pi);       % azimuthal part of wave function
Th=@(costheta) sqrt(3/2)*costheta;         % polar part of wave function, as function of cos(theta)
R=@(r) sqrt(8/3)*pi^(-1/4)*r*exp(-r^2/2);              % radial wave function (m\omega/hbar=1)
  
nphi=99;
vphi=linspace(0,2*pi,nphi+1);  % vector of phi-values
for n=1:nphi+1
	vPhi(n)=Phi(vphi(n));  % vector of Phi-values
end

nth=99;
vcosth=linspace(-1,1,nth+1);   % vector of cos(theta)-values
for n=1:nth+1
	vTh(n)=Th(vcosth(n));  % vector of Theta-values
end

phi_interval=2*pi/nphi;        % phi interval
costh_interval=2/nth;          % interval in cos(theta)
% cumulative sum of azimuthal probability:
CPphi=(cumsum(vPhi.^2)-vPhi(1).^2)*phi_interval; 
%CPphi(1)                      % checking that the cumulative sum  starts with zero
CPphi(nphi+1)=1;               % make sure that it ends with 1.
CPth=(cumsum(vTh.^2)-vTh(1).^2)*costh_interval;
%CPth(1)
CPth(nth+1)=1;
 CPth % checking monotonous increase (required for interpolation below)
% break;

vr(1)=0; CPr(1)=0; rstep1=6/100; nmax=100;

for n=2:nmax
vr(n)=vr(n-1)+rstep1;          % vector of r-values
CPr(n)=CPr(n-1)+rstep1*(vr(n)*R(vr(n)))^2; % cumulative radial probability
end
fprintf(' n= %4i rstep= %8.3f r= %8.3f CPr= %8.7f\n ', n,rstep1,vr(n),CPr(n));
CPr(nmax)=1;                    % make sure that the cumulative sum ends with 1.

NR=1000;                       % number of points
CPphirand=rand(1,NR);          % generates NR random values of integrated azimuthal probability
% Find the corresponding phi-values by interpolation:  
phirand=interp1(CPphi,vphi,CPphirand);
CPthrand=rand(1,NR);           % generates NR random values of integrated polar probability 
% Find the corresponding values of cos(theta) by interpolation:
costhrand=interp1(CPth,vcosth,CPthrand);
CPrrand=rand(1,NR);            % generates NR random values of integrated radial probability
% Find the corresponding values of r by interpolation:
rrand=interp1(CPr,vr,CPrrand);

% Cartesian coordinates of the NR positions:
z=rrand.*costhrand;
sinthrand=sqrt(1-costhrand.^2);
x=rrand.*sinthrand.*cos(phirand);
y=rrand.*sinthrand.*sin(phirand);

fullscreen = get(0,'ScreenSize');
% fullscreen is a matrix. fullscreen(1,2) are the pixel coordinates of the lower left corner, 
% fullscreen(3) is the width, fullscreen(4) is the height, all in pixels 

figure('Position',[0 0 0.5*fullscreen(3) 0.9 *fullscreen(4)])  % determines the position and size on the screen


rarray=[];
zarray=[];
sintharray=[];
phiarray=[];
deltaphi=2*pi/NR;
xx=max([abs(x) abs(y) abs(z)]);
rotate3d on;

% display positions one after one, while slowly rotating the system, to get a 3D impression:
for n=1:NR
Dphi=(n)*deltaphi+pi+0.1;
zarray=[zarray z(n)];                          % z-coordinates 1 to n (adding one at a time)
sintharray=[sintharray sinthrand(n)];
phiarray=[phiarray phirand(n)];
rarray=[rarray rrand(n)];
xarray=rarray.*sintharray.*cos(phiarray+Dphi); % x-coordinates 1 to n 
yarray=rarray.*sintharray.*sin(phiarray+Dphi); % y-coordinates 1 to n
plot3(xarray,yarray,zarray,'.');               % plots positions 1 to n
axis([-xx xx -xx xx -xx xx]);
%axis equal; 
axis off;
ax=xx*cos(Dphi); ay=xx*sin(Dphi); zzz=zeros(1:2);
tx=[-ax ax]; ty=[-ay ay];
line(tx,ty,zzz ); text(1.1*ax,1.1*ay,0,'x');   % x-axis
line(-ty,tx,zzz); text(-1.1*ay,1.1*ax,0,'y');  % y-axis
line(zzz,zzz,[-xx xx]); text(0,0,1.1*xx,'z');  % z-axis
pause(0.01);
%% view(-50,10);
M(n)=getframe;
end
%%% print -dpng fig4215til9-11.png % used together with the view statement three lines above 
                                   % to make the scatter plot at the beginning of Tillegg 5 
% nrep=1; movie(M,nrep,10); % repeat movie nrep times, with 10 frames per second
	   t1=text(0.3*xx,0,xx,'1000 points'); set(t1,'FontSize',16)
	   axis equal;






