
function demo_plot

phi=[0 pi/4 pi/2]; phi=phi';
theta=[pi/2 0 pi/4]; theta=theta';

P=[cos(phi).*sin(theta) sin(phi).*sin(theta) cos(theta)];
P1a=(P(1,:))'; P2a=(P(2,:))'; P3a=(P(3,:))';


% ... barycenter ...
RR=norm(P1a);
vert0=[P1a P2a P3a]/RR; 
CC=1/3*(P1a+P2a+P3a); CC=CC/norm(CC);

% .................... rotation matrix to the north pole ..................

[az,el,r] = cart2sph(CC(1),CC(2),CC(3));
phi=az; theta=pi/2-el;
cp=cos(phi); sp=sin(phi); ct=cos(theta); st=sin(theta);
R1=[ct 0 -st; 0 1 0; st 0 ct]; R2=[cp sp 0; -sp cp 0; 0 0 1]; 
R=R1*R2; invR=R';

% ............... vertices of the triangle at the North Pole ..............

vert1=R*vert0; 

P1=vert1(:,1);
P2=vert1(:,2);
P3=vert1(:,3);

% cap analysis
d(1)=norm(P1-[0; 0; 1]);
d(2)=norm(P2-[0; 0; 1]);
d(3)=norm(P3-[0; 0; 1]);

dmax=max(d);
thCAP=acos(dmax/2);
hCAP=dmax*cos(thCAP);
rCAP=dmax*sin(thCAP);
zCAP=1-hCAP;

XT=cub_sphpgon_QMC(10000,vert1');

sphere(20);
hold on;
plot3(XT(:,1),XT(:,2),XT(:,3),'ro','MarkerSize',3,'MarkerEdgeColor','r');

plot_cap(zCAP,rCAP);
hold off


function plot_cap(z,r)

th=linspace(0,2*pi,1000); th=th';
plot3(r*cos(th),r*sin(th),z*ones(size(th)),'m-','LineWidth',3);



function plot_side(P1,P2)

theta=linspace(0,1,1000);
for k=1:size(P1,1)
Px=theta*P1(1)+(1-theta)*P2(1);
Py=theta*P1(2)+(1-theta)*P2(2);
Pz=theta*P1(3)+(1-theta)*P2(3);
end

Pnorm=sqrt(Px.^2+Py.^2+Pz.^2);
% P=[Px./Pnorm Py./Pnorm Pz./Pnorm];
plot3(Px./Pnorm,Py./Pnorm,Pz./Pnorm,'k-','LineWidth',3);


