% définition du cercle à l'équateur
t=0:pi/100:2*pi;
r=3;
C=[r*cos(t)
   r*sin(t)
   zeros(size(t))];

% rotation autour de Ox de 55°
a=deg2rad(55);
Rx=[1   0       0
    0 cos(a) -sin(a)
	0 sin(a)  cos(a)];
C=Rx*C;

% rotations autour de Oz de 60°
b=pi/3;
Rz=[cos(b) -sin(b) 0
    sin(b)  cos(b) 0
	  0       0    1];
	  
% représentation graphique
% terre
load topo
bgcolor=[.8 .8 .8];
f.fig = figure('Tag','Earthmap','colormap',topomap1,'Color',bgcolor);

[x,y,z] = sphere(50);
hg.axes3.DataAspectRatio = [1 1 1];
hg.axes3.PlotBoxAspectRatioMode = 'auto';
   
f.ax(3) = axes('Tag','3D','Visible','off', hg.axes3);
   
mat.dull.AmbientStrength = 0.1;
mat.dull.DiffuseStrength = 1;
mat.dull.SpecularColorReflectance = .5; 
mat.dull.SpecularExponent = 20;
mat.dull.SpecularStrength = 1;
f.surf = surface(x,y,z, ...
                  mat.dull, ...
                  'FaceColor','texture',...
                  'EdgeColor','none',...
                  'FaceLighting','phong',...
                  'Cdata',topo,...
                  'Parent',f.ax(3));
f.lite(1) = light('position',[-1 0 1]);
f.lite(2) = light('position',[-1.5 0.5 -0.5], 'color', [.6 .2 .2]);
view(3) 

% satellites
ind=[  1  51 101 151
      41  91 141 191
	  81 131 181  30
	 121 171  20  70
	 161 201  50 100
	 181  30  80 130];
hold on
for i=0:5
	C=Rz*C;
	plot3(C(1,:),C(2,:),C(3,:))
	j=ind(i+1,:);
	plot3(C(1,j),C(2,j),C(3,j),'or')
end

axis equal
axis off
hold off


