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

// S.Ayrinhac (02/2012)
// lambda : longitude en dégrés décimaux
// phi : latitude en degrés décimaux

format("e",16) //format long e

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

// [lambda,phi]=Projection_Lambert93(1029705.083,272723.849)
// [lambda,phi]=Projection_Lambert93(600043,7087060)

//// système géodésique      : RGF93
//// ellipsoide de référence : GRS80

//Ellipsoid 	Semi-major axis ( a ) 	Inverse flattening ( 1 / f )
//GRS80 	6 378 137 m 	298.257 222 101
//International 1924 	6 378 388 m 	297
//WGS84 	6 378 137 m 	298.257 223 563

//// CoordSys Earth Projection 3, 999, 0, 0, 0  , 0,  "m", 3, 46.5, 44, 49.0, 700000, 6600000
//// "Lambert 93",             3, 999, 0, 0 , 0 , 0 , 7 , 3 , 46.5 , 44.0, 49.0, 700000, 6600000
//// BD carthage : CoordSys E  3, 999, 0, 0.0000, 0.0000, 0.0000, "m", 3.000000000000, 46.5, 44.0, 49.0, 700000.000, 6600000.000

Phi1=deg2rad(49);   //rad
Phi2=deg2rad(44);   //rad
lambda0=deg2rad(3); //rad
Phi0=deg2rad(46.5); //rad
a=6378137.0;        //demi-grand axe de l'ellipsoide
f=1/298.257222101;//aplatissment
b=a.*(1-f);       //demi-petit axe
e=sqrt(2*f-f^2);//1.55e-4;// //0.08181919106; //e=sqrt(a.^2-b.^2)./a
     // http://www.umr-lisah.fr/CrSig/Documents/FichesTech/Coordonn//C3//A9es/projfr.pdf
     // http://sgcaf.free.fr/pages/techniques/ign_coordonnees.htm
     // par défaut j'ai e=0
x0=700000;
y0=6600000;
epsilon=1e-15;
//n=0.7256077650;


//////////////////////////////////////
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);
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));

F=m1./(n.*t1.^n);
rho0=a.*F.*t0.^n;

////calcul de phi //////////////////////////////////////////////
E=x;
E0=x0;
Ep=E-E0;

N0=y0;
N=y;
Np=N-N0;

rhop=sign(n).*sqrt(Ep.^2+(rho0-Np).^2);
tp=(rhop./(a.*F)).^(1./n);
gammap=atan(Ep./(rho0-Np));

lambda=gammap./n+lambda0;


phix(1)=%pi./2-2.*atan(tp);
i=2;
while 1~=0
    phip=phix(i-1);
    
    G=((1-e.*sin(phip))./(1+e.*sin(phip))).^e2;
    phix(i)=%pi/2-2.*atan(tp.*G);
    if abs(phix(i)-phix(i-1))<epsilon
        break
    else
        i=i+1;
    end
end
phi=phix(i);

//conversion en degrés décimaux ////////////////
phi=rad2deg(phi);
lambda=rad2deg(lambda);

endfunction