%   this is a demo for spherical coordinates
% 
 
o = [0,0,0];   % origin will be shifted otherwise we won't see the arrows inside our ellipsoid

 arrow3(o, [1,0,0],'yellow');
 arrow3(o, [0,1,0],'yellow');
 arrow3(o, [0,0,1],'magenta');    % this plots the i,j,k, vectors
 
 zmin = 0;
 zmax = 4;       % plotting window
 
 
   %%  Example 1. Half a sphere with markings
   %
     
 dtheta = pi/48;  
 dhpi = dtheta;
 rho = 1;
 
 [theta, phi] = meshgrid(0:dtheta:2*pi,  0:dphi:(pi/2));       % get whole sphere if 0 < phi < pi
 
   X = rho*cos(theta).*sin(phi);
   Y = rho*sin(theta).*sin(phi);
   Z = rho*cos(phi);
          surf(X,Y,Z,'FaceColor',[0.8 1 0.8]);        % surface plotting
 
     %  a)  showing various traces phi = const
     % 
   
     t = 0:dtheta:(2*pi);
     theta = t;
     phi = pi/6;
     
   x = rho*cos(theta)*sin(phi);
   y = rho*sin(theta)*sin(phi);
   z = rho*cos(phi) + 0*x;     
   
     hold on
     plot3(x,y,z,'r',...
                'LineWidth',2);
    
     phi = pi/4;
   x = rho*cos(theta)*sin(phi);
   y = rho*sin(theta)*sin(phi);
   z = rho*cos(phi) + 0*x;     
   
     plot3(x,y,z,'m',...
                'LineWidth',2);
  
    %
    %   now adding slices where theta =const 
    
    
     t = 0:dphi:(pi/2);
     theta = 0;
     phi = t;
     
   x = rho*cos(theta)*sin(phi);
   y = rho*sin(theta)*sin(phi);
   z = rho*cos(phi) + 0*x;     
   
     hold on
     plot3(x,y,z,...
                'LineWidth',2);
    
    
     t = 0:dphi:(pi/2);
     theta = pi/2;
     phi = t;
     
   x = rho*cos(theta)*sin(phi);
   y = rho*sin(theta)*sin(phi);
   z = rho*cos(phi) + 0*x;     
   
     hold on
     plot3(x,y,z,...
                'LineWidth',2);