function [x,y]=Projection_Lambert93(lambda,phi)

// Fonctions de conversion des angles ///////////////////
deff('angle=rad2deg(angle0)','angle=(angle0./%pi).*180');
deff('angle=deg2rad(angle0)','angle=(angle0./180).*%pi');

lambda=(lambda.*%pi)./(180);
phi=(phi.*%pi)./(180);

Phi0=46.5*%pi/180; //rad
Phi1=49*%pi/180;   //rad
Phi2=44*%pi/180;   //rad
lambda0=3*%pi/180; //rad
f=1/298.257222101;//aplatissment
a=6378137.0;        //demi-grand axe de l'ellipsoide
e=sqrt(2*f-f^2);
x0=700000;
y0=6600000;

// Calcul de n //////////////////////////
S0=e.*sin(Phi0);
S1=e.*sin(Phi1);
S2=e.*sin(Phi2);
e2=e/2;
m0=cos(Phi0)./sqrt(1-S0.^2);
m1=cos(Phi1)./sqrt(1-S1.^2);
m2=cos(Phi2)./sqrt(1-S2.^2);
t=tan((%pi/4)-(phi./2))./((1-S0)./(1+S0)).^e2;
t0=tan((%pi/4)-(Phi0/2))./((1-S0)./(1+S0)).^e2;
t1=tan((%pi/4)-(Phi1/2))./((1-S1)./(1+S1)).^e2;
t2=tan((%pi/4)-(Phi2/2))./((1-S2)./(1+S2)).^e2;
n=(log(m1)-log(m2))./(log(t1)-log(t2));
// Calcul de rho0 //////////////////////
F=m1./(n.*t1.^n);
rho0=a.*F.*t0.^n;
rho=a.*F.*t.^n;
/////////////////////////////////////////


gama=n.*(lambda-lambda0);
y=y0+rho0-rho.*cos(gama);
x=x0+rho.*sin(gama);


endfunction