% rotator.m - scatter plot of rotator states
clear all; close all;

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.7*fullscreen(3) 0.7 *fullscreen(4)])  % determines the position and size on the screen

subplot(1,2,1);                % s distribution
Phi=@(phi) 1/sqrt(2*pi);       % azimuthal part of wave function
Th=@(theta) 1/sqrt(2);         % polar part of wave function
%Th=@(costheta) costheta*sqrt(1.5); % polar wave function as f. of cos(theta)
npr=1;                         % principal quantum number
R=@(r) 2*exp(-r);      % radial wave function (a=1)
%R=@(r) r*exp(-r/2)/sqrt(24);
  
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)
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)



NR=500;                 % 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);

r=4; %radius
% Cartesian coordinates:
z=r*costhrand;
sinthrand=sqrt(1-costhrand.^2);
x=r*sinthrand.*cos(phirand);
y=r*sinthrand.*sin(phirand);

hold on;
%plot3(x,y,z,'.');
rotate3d on;
axis equal;
view(-230,10);
axis off;
na=2; xx=1.5*r; t=linspace(-xx,xx,na);
zzz=zeros(1,na);
line(t,zzz,zzz); text(xx+0.2,0,0,'x'); 
line(zzz,t,zzz); text(0,xx+0.2,0,'y');
line(zzz,zzz,t); text(0,0,xx+0.2,'z');
for it=1:length(x)
	 plot3(x(it),y(it),z(it),'.');
pause(0.01);
end
xe=r*cos(vphi); ye=r*sin(vphi);
line(xe,ye,vphi-vphi);
hold off;




subplot(1,2,2);                % pz distribution
Phi=@(phi) 1/sqrt(2*pi);       % azimuthal part of wave function
%Th=@(theta) 1/sqrt(2);         % polar part of wave function
Th=@(costheta) costheta*sqrt(1.5); % polar wave function as f. of cos(theta)
npr=1;                         % principal quantum number
R=@(r) 2*exp(-r);      % radial wave function (a=1)
%R=@(r) r*exp(-r/2)/sqrt(24);
  
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)
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)



%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);

r=4; %radius
% Cartesian coordinates:
z=r*costhrand;
sinthrand=sqrt(1-costhrand.^2);
x=r*sinthrand.*cos(phirand);
y=r*sinthrand.*sin(phirand);



hold on;
%plot3(x,y,z,'.');
rotate3d on;
axis equal;
view(-230,10);
axis off;
na=2; xx=1.5*r; t=linspace(-xx,xx,na);
zzz=zeros(1,na);
line(t,zzz,zzz); text(xx+0.2,0,0,'x'); 
line(zzz,t,zzz); text(0,xx+0.2,0,'y');
line(zzz,zzz,t); text(0,0,xx+0.2,'z');
for it=1:length(x)
	 plot3(x(it),y(it),z(it),'.');
pause(0.01);
end
xe=r*cos(vphi); ye=r*sin(vphi);
line(xe,ye,vphi-vphi);
hold off;
% print -dpng figrotator.png

break;
[xs,ys,zs]=sphere(40);
xsr=r*xs; ysr=r*ys; zsr=r*zs;
surf(xsr,ysr,zsr);
