% Definiera området i r-phi-theta-koor
phi1 = pi/4;
theta1 = pi/4;
r1= 1;
phi2 = 1.2*phi1;
theta2 = 1.2*theta1;
r2 = 1.2*r1;
phi=phi1:pi/80:phi2;
theta=theta1:pi/80:theta2;
r = r1:.025:r2;

% Beräkna hörnpunkter i x1-x2-x3-koor
x1=r1.*sin(phi1).*cos(theta1);
y1=r1.*sin(phi1).*sin(theta1);
z1=r1.*cos(phi1);
x2=r1.*sin(phi2).*cos(theta1);
y2=r1.*sin(phi2).*sin(theta1);
z2=r1.*cos(phi2);
x3=r1.*sin(phi1).*cos(theta2);
y3=r1.*sin(phi1).*sin(theta2);
z3=r1.*cos(phi1);

% Rita de sex randytorna till området
hold off
colormap('gray')
% r=r1 konst
[PHI,THETA]=meshgrid(phi,theta);
X=r1.*sin(PHI).*cos(THETA);
Y=r1.*sin(PHI).*sin(THETA);
Z=r1.*cos(PHI);
S1=surfl(X,Y,Z);
hold on
% r=r2 konst
X=r2.*sin(PHI).*cos(THETA);
Y=r2.*sin(PHI).*sin(THETA);
Z=r2.*cos(PHI);
S2=surfl(X,Y,Z);
% phi=phi1 konst
[R,THETA]=meshgrid(r,theta);
X=R.*sin(phi1).*cos(THETA);
Y=R.*sin(phi1).*sin(THETA);
Z=R.*cos(phi1);
S3=surfl(X,Y,Z);
% phi=phi2 konst
X=R.*sin(phi2).*cos(THETA);
Y=R.*sin(phi2).*sin(THETA);
Z=R.*cos(phi2);
S4=surfl(X,Y,Z);
% theta=theta1 konst
[R,PHI]=meshgrid(r,phi);
X=R.*sin(PHI).*cos(theta1);
Y=R.*sin(PHI).*sin(theta1);
Z=R.*cos(PHI);
S5=surfl(X,Y,Z);
% theta=theta2 konst
X=R.*sin(PHI).*cos(theta2);
Y=R.*sin(PHI).*sin(theta2);
Z=R.*cos(PHI);
S6=surfl(X,Y,Z);

% Linjen till punkten (r1,phi1,theta1)
l=plot3([0 x1],[0 y1],[0 z1],'g*-');
set(l,'linewidth',2)
t=text(.7*r1*sin(.98*phi1)*cos(theta1),.7*r1*sin(.98*phi1)*sin(theta1),.7*r1*cos(.98*phi1),'R');
set(t,'color','g')

% Lite stödlinjer
plot3([0 x2],[0 y2],[0 z2],'k--')
plot3([0 x3],[0 y3],[0 z3],'k--')
plot3([0 x1],[0 y1],[0 0],'r--')
t=text(.6*r1*sin(phi1)*cos(.98*theta1),.6*r1*sin(phi1)*sin(.98*theta1),0,'R sin(\phi)');
set(t,'color','r')
plot3([0 x3],[0 y3],[0 0],'b--')
plot3([x1 x1],[y1 y1],[0 z1],'b-.')
plot3([x3 x3],[y3 y3],[0 z3],'b-.')
plot3([0 0],[0 0],[0 1],'k')

% Vinkelmarkeringarna
plot(.15*cos(0:.1*theta1:theta1),.15*sin(0:.1*theta1:theta1),'b')
t=text(.17*cos(.5*theta1),.17*sin(.5*theta1),0,'\theta');
set(t,'color','b')
plot(.15*cos(theta1:.1*(theta2-theta1):theta2),.15*sin(theta1:.1*(theta2-theta1):theta2),'b')
t=text(.2*cos(.5*(theta2-theta1)+theta1),.2*sin(.5*(theta2-theta1)+theta1),0,'\Delta\theta');
set(t,'color','b')
plot3(.15*sin(0:.1*phi1:phi1)*cos(theta1),.15*sin(0:.1*phi1:phi1)*cos(theta1),.15*cos(0:.1*phi1:phi1),'k')
text(.2*sin(.5*phi1)*cos(theta1),.17*sin(.5*phi1)*sin(theta1),.17*cos(.5*phi1),'\phi')
plot3(.15*sin(phi1:.1*(phi2-phi1):phi2)*cos(theta1),.15*sin(phi1:.1*(phi2-phi1):phi2)*cos(theta1),...
    .15*cos(phi1:.1*(phi2-phi1):phi2),'b')
text(.2*sin(.5*(phi2-phi1)+phi1)*cos(theta1),.2*sin(.5*(phi2-phi1)+phi1)*sin(theta1), ...
    .2*cos(.5*(phi2-phi1)+phi1),'\Delta\phi')

plot(r1*sin(phi1)*cos(theta1:.1*(theta2-theta1):theta2),r1*sin(phi1)*sin(theta1:.1*(theta2-theta1):theta2),'r--')
t=text(1.1*r1*sin(phi1)*cos(.5*(theta2-theta1)+theta1),1.1*r1*sin(phi1)*sin(.5*(theta2-theta1)+theta1),0,...
    'R sin(\phi) \Delta\theta');
set(t,'color','r')

% Kanterna i volymen och deras längder
l=plot3(r1*sin(phi1)*cos(theta1:.1*(theta2-theta1):theta2), ...
    r1*sin(phi1)*sin(theta1:.1*(theta2-theta1):theta2), ...
    r1*cos(phi1)*ones(size(theta1:.1*(theta2-theta1):theta2)),'r');
set(l,'linewidth',3)
t=text(r1*sin(.98*phi1)*cos(.5*(theta2-theta1)+theta1),r1*sin(.98*phi1)*sin(.5*(theta2-theta1)+theta1),...
    r1*cos(.98*phi1),'R sin(\phi) \Delta\theta');
set(t,'color','r')
l=plot3(r1*sin(phi1:.1*(phi2-phi1):phi2)*cos(theta1),r1*sin(phi1:.1*(phi2-phi1):phi2)*cos(theta1),...
    r1*cos(phi1:.1*(phi2-phi1):phi2),'r');
set(l,'linewidth',3)
t=text(r1*sin(.5*(phi2-phi1)+phi1)*cos(.98*theta1),r1*sin(.5*(phi2-phi1)+phi1)*sin(.98*theta1),...
    r1*cos(.5*(phi2-phi1)+phi1),'R \Delta\phi');
set(t,'color','r')
l=plot3([r1,r2]*sin(phi1)*cos(theta1),[r1,r2]*sin(phi1)*cos(theta1),...
    [r1,r2]*cos(phi1),'r');
set(l,'linewidth',3)
t=text((.5*(r2-r1)+r1)*sin(.98*phi1)*cos(.98*theta1),(.5*(r2-r1)+r1)*sin(.98*phi1)*sin(.98*theta1),...
    (.5*(r2-r1)+r1)*cos(.98*phi1),'\DeltaR');
set(t,'color','r')

% Axlar och synvinkel
set(gca,'xlim',[0 .8],'ylim',[0 .8], 'zlim',[0 .9])
xlabel('x_1')
ylabel('x_2')
zlabel('x_3')
view(-6,26)
grid on

rotate3d on
