% lecture des données
load coast
t1=-pi/2:pi/12:pi/2;
t2=-pi:2*pi/12:pi;
t3=-pi:pi/50:pi;
t4=-pi/2:pi/100:pi/2;
n=max(size(t1));
m=max(size(t4));
parax=(ones(n,1)*t3)';
paray=ones(m,1)*t1;
merix=ones(m,1)*t2;
meriy=(ones(n,1)*t4)';
paraz=pi/2-paray;	% colatitude pour simplifier
meriz=pi/2-meriy;

% conversion degrés -> radians
latr=deg2rad(lat);	% theta
longr=deg2rad(long); % phi
z=pi/2-latr;

% projection cylindrique perspective

% parallèles et méridiens
plot(parax,tan(paray),':m')
hold on
plot(merix,tan(meriy),':m')
% contours
x=longr;			% la longitude est la même
y=tan(latr);		% sur une sphère la tangente donne la verticale du cylindre
plot(x,y)
title('Projection cylindrique perspective')


axis equal
axis([-4 4 -3 3])
axis off
set(gcf,'color','w')
hold off
figure(gcf)
pause

% projection quadratique (cylindrique équidistante)

% parallèles et méridiens
plot(parax,paray,':m')
hold on
plot(merix,meriy,':m')
% contours
x=longr;			
y=latr;		
plot(x,y)
title('Projection quadratique')


axis equal
axis off
hold off
pause

% projection isocylindrique de Lambert (cylindrique équivalente)

% parallèles et méridiens
plot(parax,sin(paray),':m')
hold on
plot(merix,sin(meriy),':m')
% contours
x=longr;			
y=sin(latr);	
plot(x,y)
title('Projection cylindrique de Lambert')


axis equal
axis off
hold off
pause

% projection cylindrique conforme (Mercator)

% parallèles et méridiens
plot(parax,log(tan(paray/2+pi/4)+eps),':m')
hold on
plot(merix,log(tan(meriy/2+pi/4)+eps),':m')
% contours
x=longr;			
y=log(tan(latr/2+pi/4)+eps);
plot(x,y)
title('Projection de Mercator')

axis equal
axis([-4 4 -3 3])
axis off
hold off
pause

% projection méricylindrique équivalente (sinusoïdale de Sanson)

% parallèles et méridiens
plot(parax.*cos(paray),paray,':m')
hold on
plot(merix.*cos(meriy),meriy,':m')
% contours
x=longr.*cos(latr);			
y=latr;		
plot(x,y)
title('Projection de sinusoïdale (Sanson)')

axis equal
axis off
hold off
pause

% projection azimutale perspective

% parallèles et méridiens
[x,y]=pol2cart(parax,tan(paraz));
plot(x,y,':m')
hold on
[x,y]=pol2cart(merix,tan(meriz));
plot(x,y,':m')
% contours
[x,y]=pol2cart(longr,tan(z));
i=find(latr<0);
y(i)=NaN.*y(i);
plot(x,y)
title('Projection gnomonique')

axis equal
axis([- 3 3 -3 3])
axis off
hold off
pause

% projection azimutale perspective

% parallèles et méridiens
[x,y]=pol2cart(parax,sin(paraz));
plot(x,y,':m')
hold on
[x,y]=pol2cart(merix,sin(meriz));
plot(x,y,':m')
% contours
[x,y]=pol2cart(longr,sin(z));
i=find(latr<0);
y(i)=NaN.*y(i);
plot(x,y)
title('Projection orthographique')

axis equal
axis off
hold off
pause

% projection azimutale équidistante (Guillaume Postel)

% parallèles et méridiens
[x,y]=pol2cart(parax,paraz);
plot(x,y,':m')
hold on
[x,y]=pol2cart(merix,meriz);
plot(x,y,':m')
% contours
[x,y]=pol2cart(longr,z);
plot(x,y)
title('Projection de Guillaume Postel')

axis equal
axis([-3 3 -3 3])
axis off
hold off
pause

% projection azimutale équivalente (azimutale de Lambert)

% parallèles et méridiens
[x,y]=pol2cart(parax,sqrt(2-2*cos(paraz)));
plot(x,y,':m')
hold on
[x,y]=pol2cart(merix,sqrt(2-2*cos(meriz)));
plot(x,y,':m')
% contours
[x,y]=pol2cart(longr,sqrt(2-2*cos(z)));
plot(x,y)
title('Projection azimutale de Lambert')

axis equal
axis off
hold off
pause

% projection azimutale conforme (stéréographique)

% parallèles et méridiens
[x,y]=pol2cart(parax,2*tan(paraz./2));
plot(x,y,':m')
hold on
[x,y]=pol2cart(merix,2*tan(meriz./2));
plot(x,y,':m')
% contours
[x,y]=pol2cart(longr,2*tan(z./2));
plot(x,y)
title('Projection stéréographique')

axis equal
axis([-3 3 -3 3])
axis off
hold off
pause

% projection conique équidistante de Delisle

% parallèles et méridiens
t1=pi/4;
t2=pi/3;
n=(cos(t1)-cos(t2))/(t2-t1);
r=(t2*cos(t1)-t1*cos(t2))/(cos(t1)-cos(t2));
[x,y]=pol2cart(n*parax,r-paray);
plot(x,y,':m')
hold on
[x,y]=pol2cart(n*merix,r-meriy);
plot(x,y,':m')
% contours
[x,y]=pol2cart(n*longr,r-latr);
plot(x,y)
title('Projection conique de Delisle')

axis equal
axis off
hold off
pause

% projection conique équivalente d'Albers

% parallèles et méridiens
t1=pi/4;
t2=pi/3;
n=(cos(t1)^2-cos(t2)^2)/(2*(sin(t2)-sin(t1)));
C=(cos(t1)^2+2*n*sin(t1))/(n^2);
[x,y]=pol2cart(n*parax,sqrt(C-(2/n)*sin(paray)));
plot(x,y,':m')
hold on
[x,y]=pol2cart(n*merix,sqrt(C-(2/n)*sin(meriy)));
plot(x,y,':m')
% contours
[x,y]=pol2cart(n*longr,sqrt(C-(2/n)*sin(latr)));
plot(x,y)
title('Projection conique d''Albers')

axis equal
axis off
hold off
pause

% projection conique conforme de Lambert

% parallèles et méridiens
t1=pi/4;
t2=pi/3;
n=(log(cos(t1)-log(cos(t2))))/(log(tan(pi/4-t1/2))-log(tan(pi/4-t2/2)));
k=(cos(t1)/n)*(tan(pi/4+t1/2))^n;
[x,y]=pol2cart(n*parax,k*tan(pi/4-paray/2).^n);
plot(x,y,':m')
hold on
[x,y]=pol2cart(n*merix,k*tan(pi/4-meriy/2).^n);
plot(x,y,':m')
% contours
[x,y]=pol2cart(n*longr,k*tan(pi/4-latr/2).^n);
plot(x,y)
title('Projection conique de Lambert')

axis equal
axis([-3 3 -3 3])
axis off
hold off
pause


