clear all;clf

n=20;
theta=(-n:2:n)/n*pi;
phi=(-n:2:n)'/n*pi/2;
cosphi=cos(phi);cosphi(1)=0;cosphi(n+1)=0;
sintheta=sin(theta);sintheta(1)=0;sintheta(n+1)=0;

x=cosphi*cos(theta);
y=cosphi*sintheta;
z=sin(phi)*ones(1,n+1);

figure(1);clf;whitebg('w');hold on
view(3);
C1=[-0.5 0 1.5];R1=0.8;
Sat1=surf(0.1*x+C1(1),0.1*y+C1(2),0.1*z+C1(3),1*ones(size(z)));
S1=surf(R1*x+C1(1),R1*y+C1(2),R1*z+C1(3),10*ones(size(z)));

C2=[0.5 0 1.5];R2=0.9;
Sat2=surf(0.1*x+C2(1),0.1*y+C2(2),0.1*z+C2(3),1*ones(size(z)));
S2=surf(R2*x+C2(1),R2*y+C2(2),R2*z+C2(3),10*ones(size(z)));

%% Resultat avec 3eme satellite : 
% localisation
loc=1;
if loc==1
C3=[0. -1. 1.];R3=0.9;
Sat3=surf(0.1*x+C3(1),0.1*y+C3(2),0.1*z+C3(3),1*ones(size(z)));
S3=surf(R3*x+C3(1),R3*y+C3(2),R3*z+C3(3),10*ones(size(z)));

Cr=[0.5 0.5 0.5];Rr=0.1;
Sr=surf(Rr*x+Cr(1),Rr*y+Cr(2),Rr*z+Cr(3),1*ones(size(z)));
end


shading interp
alpha(0.4)
alpha(Sat1,1)
alpha(Sat2,1)
alpha(Sat3,1);
if loc==1;alpha(Sr,1);end
camlight;
lighting gouraud;
axis square equal
%axis off

%Intersec(C1,R1,C2,R2,10)

tht=0:0.1:2*pi;
xI=(R2*C1(1)+R1*C2(1))/(R1+R2);%xI=-0.4;
yI=(R2*C1(2)+R1*C2(2))/(R1+R2);%yI=0;
zI=(R2*C1(3)+R1*C2(3))/(R1+R2);%zI=1.5;

D1=(xI-C1(1))^2+(yI-C1(2))^2+(zI-C1(3))^2;
D1=sqrt(D1);
D2=(xI-C2(1))^2+(yI-C2(2))^2+(zI-C2(3))^2;
D2=sqrt(D2);
R=R1^2-D1^2;R=sqrt(R)
%R=R2^2-D2^2;R=sqrt(R)


%plot3(xI,yI,zI,'.r')
%plot3(xI+zeros(size(tht)),yI+R*cos(tht),zI+R*sin(tht),'-r')

pp=plot3(-0.5,-0.5,1.2,'xr');
set(pp,'markersize',10);
set(pp,'linewidth',2.5)
%% Earthmap : 
load('topo.mat','topo','topomap1');
contour(0:359,-89:90,topo,[0 0],'b')

[x,y,z] = sphere(50);

props.AmbientStrength = 0.1;
props.DiffuseStrength = 1;
props.SpecularColorReflectance = .5; 
props.SpecularExponent = 20;
props.SpecularStrength = 1;
props.FaceColor= 'texture';
props.EdgeColor = 'none';
props.FaceLighting = 'phong';
props.Cdata = topo;
surface(x,y,z,props);
light('position',[-1 0 1]);
light('position',[-1.5 0.5 -0.5], 'color', [.6 .2 .2]);
view(3)

axis([-1. 1. -1.5 1 -1 3])
axis off