% Spherical Coordinates
% 2021/11/22 by T. Hirano
clear; clc;
%
fig = figure;
hold on;
xmin=-2; xmax=2;
ymin=-2; ymax=2;
zmin=-2; zmax=2;
axis([xmin xmax ymin ymax zmin zmax]);
pbaspect([1 1 1]);
axis off;
view(60,20);
camzoom(1.7);
% Draw axes
quiver3(-2,0,0,4,0,0,0,'k','filled'); % x-axis
quiver3(0,-2,0,0,4,0,0,'k','filled'); % y-axis
quiver3(0,0,-2,0,0,4,0,'k','filled'); % z-axis
text(0,0,0,'O');
text(xmax,0,0,'x');
text(0,ymax,0,'y');
text(0,0,zmax,'z');
%
x = 1;
y = 0.5;
z = 1.5;
r=sqrt(x^2+y^2+z^2);
theta=atan2(sqrt(x^2+y^2),z);
phi=atan2(y,x);
quiver3(0,0,0,x,y,z,0,'r','filled'); % Position Vector
% Unit vectors
quiver3(x,y,z,(r/3)*sin(theta)*cos(phi),(r/3)*sin(theta)*sin(phi),(r/3)*cos(theta),0,'b','filled'); % r
quiver3(x,y,z,(r/3)*cos(theta)*cos(phi),(r/3)*cos(theta)*sin(phi),-(r/3)*sin(theta),0,'b','filled'); % theta
quiver3(x,y,z,-(r/3)*sin(phi),(r/3)*cos(phi),0,0,'b','filled'); % phi
text(x+(r/3)*sin(theta)*cos(phi),y+(r/3)*sin(theta)*sin(phi),z+(r/3)*cos(theta),'r\^');
text(x+(r/3)*cos(theta)*cos(phi),y+(r/3)*cos(theta)*sin(phi),z-(r/3)*sin(theta),'theta\^');
text(x-(r/3)*sin(phi),y+(r/3)*cos(phi),z,'phi\^');
% Aux. lines
l1_1=line([0 0 x x x 0 0 x x],[y 0 0 0 y y y y 0],[z z z 0 0 0 z z z]);
l1_2=line([0 x x],[0 y y],[0 0 z]);
l1_1.LineStyle = '--';
l1_2.LineStyle = '--';
l1_1.Color = [0.5 0.5 0.5];
l1_2.Color = [0.5 0.5 0.5];
% Draw sphere
r1=0.95*r;
t1=atan2(sqrt(r^2-r1^2),r1);
np=12; nt=24;
for i=1:np
    p=(i-1)*(2*pi)/np;
    t=t1:(pi-2*t1)/nt:pi-t1;
    l2(i)=line(r.*sin(t).*cos(p),r.*sin(t).*sin(p),r.*cos(t));
    l2(i).LineStyle = '--';
    l2(i).Color = [0.5 0.5 0.5];
end
nz=8; np=32;
for i=1:nz+1
    zc=-r1+(i-1)*(2*r1)/nz;
    rc=sqrt(r^2-zc^2);
    p=0:2*pi/np:2*pi;
    l3(i)=line(rc.*cos(p),rc.*sin(p),zc*ones([1 length(p)]));
    l3(i).LineStyle = '--';
    l3(i).Color = [0.5 0.5 0.5];
end
% Arc for angle
% theta
na=24;
for i=1:na
    t2=pi/2:-theta/na:pi/2-theta;
    p_tmp=[(sqrt(x^2+y^2)/2).*cos(t2); zeros([1 length(t2)]); (sqrt(x^2+y^2)/2).*sin(t2)]; % Draw arc in x-z plane
    mrot_z = [cos(phi), -sin(phi), 0; sin(phi), cos(phi), 0; 0, 0, 1]; % Rotation matrix along z
    p_tmp=mrot_z*p_tmp; % Rotate points along z
    l4(i)=line(p_tmp(1,:),p_tmp(2,:),p_tmp(3,:));
    l4(i).Color = 'm';
end
% phi
na=24;
for i=1:na
    p=0:phi/na:phi;
    l5(i)=line((sqrt(x^2+y^2)/2).*cos(p),(sqrt(x^2+y^2)/2).*sin(p),zeros([1 length(p)]));
    l5(i).Color = 'm';
end