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