//Programmes de conversion de coordonnées géographiques
function utm2geo(xutm,yutm) {
//
//Ellipsoide WGS84
//
a2 = 6378137.0;
f2 = 1/298.257223563;
b2 = a2*(1-f2);
e2 = Math.sqrt(Math.pow(a2,2)-Math.pow(b2,2))/a2;
//
//Calcul des parametres de projection
//
fuseau = 32;
//Uniquement pour la region des Alpes
longitude0 = (6*fuseau-183)*Math.PI/180;
Xs = 500000;
Ys = 0;
//Hémisphre Nord seulement
k0 = 0.9996;
n = k0*a2;
longitude_c = longitude0;
//
xe1 = Math.pow(e2,0);
xe2 = Math.pow(e2,2);
xe3 = Math.pow(e2,4);
xe4 = Math.pow(e2,6);
xe5 = Math.pow(e2,8);
//
//Coefficients de projection inverse
//
C1 = xe1 - xe2/4 - 3*xe3/64 - 5*xe4/256 - 175*xe5/16384;
C2 = xe2/8 + xe3/48 + 7*xe4/2048 + xe5/61440;
C3 = xe3/768 + 3*xe4/1280 + 559*xe5/368640;
C4 = 17*xe4/30720 + 283*xe5/430080;
 C5 = 4397*xe5/41287680;
//
zpy = (yutm - Ys)/(n*C1);
zpx = (xutm - Xs)/(n*C1);

//
latiso = zpy - C2*Math.sin(2*zpy)*cosh(2*zpx) - C3*Math.sin(2*2*zpy)*cosh(2*2*zpx) - C4*Math.sin(2*3*zpy)*cosh(2*3*zpx) - C5*Math.sin(2*4*zpy)*cosh(2*4*zpx);
latiso_s = zpx - C2*Math.cos(2*zpy)*sinh(2*zpx) - C3*Math.cos(2*2*zpy)*sinh(2*2*zpx) - C4*Math.cos(2*3*zpy)*sinh(2*3*zpx) - C5*Math.cos(2*4*zpy)*sinh(2*4*zpx);
//
//Latitude isometrique inverse
//
longitude = longitude_c + Math.atan(sinh(latiso_s)/Math.cos(latiso));
phi = Math.asin(Math.sin(latiso)/cosh(latiso_s));
liso = Math.log(Math.tan(Math.PI/4 + phi/2));
//
eps = 0.0000000000001;
testconv = 1;
phi0 = eps;
while (testconv > eps) {
phi1 = 2*Math.atan(Math.pow((1 + e2*Math.sin(phi0))/(1 - e2*Math.sin(phi0)),e2/2)*Math.exp(liso)) - Math.PI/2;
testconv = Math.abs((phi1 - phi0)/phi0);
phi0 = phi1;
}
latitude = phi0; 
//
return latitude;
return longitude;}

function geo2lambert(longitude,latitude) {
//
//Parametres lambert IIe
//
xs = 600000.000;
ys = 8199695.768;
n = 0.7289686274;
C = 11745793.39;
lambda0 = 0.04079234433;
//
//Ellipsoide Clarke 1880
//
a1 = 6378249.2;
f1 = 1/293.466021;
b1 = a1*(1-f1);
e1 = Math.sqrt(Math.pow(a1,2)-Math.pow(b1,2))/a1;
//
//Ellipsoide WGS84
//
a2 = 6378137.0;
f2 = 1/298.257223563;
b2 = a2*(1-f2);
e2 = Math.sqrt(Math.pow(a2,2)-Math.pow(b2,2))/a2;
//
phi = latitude;
lambda = longitude;
h = 0;
//
W = Math.sqrt(1 - Math.pow(e2,2)*Math.pow(Math.sin(phi),2));
N = a2/W;
//
X = (N + h)*Math.cos(phi)*Math.cos(lambda);
Y = (N + h)*Math.cos(phi)*Math.sin(lambda);
Z = (N*(1 - Math.pow(e2,2)) + h)*Math.sin(phi);

//
X1 = X + 168;
Y1 = Y + 60;
Z1 = Z - 320;

//
lambda = 2*Math.atan(Y1/(X1 + Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2))));
phi0 = Math.atan(Z1/(Math.pow(X1,2) + Math.pow(Y1,2)*(1 - a1*Math.pow(e2,2)/Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2) + Math.pow(Z1,2)))));
testconv = 1;
eps = 0.00000000001;

//
while (testconv > eps) {
phi1 = Math.atan((Z1/Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2)))/(1 - a1*Math.pow(e1,2)*Math.cos(phi0)/Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2))/Math.sqrt(1 - Math.pow(e1,2)*Math.pow(Math.sin(phi0),2))));
testconv = Math.abs((phi1 - phi0)/phi0);
phi0 = phi1;
}
phi = phi0;
//
h = Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2))/Math.cos(phi) - a1/Math.sqrt(1 - Math.pow(e1,2)*Math.pow(Math.sin(phi),2));
//
//Projection Lambert 
//
isolat = 0.5*Math.log((1 + Math.sin(phi))/(1 - Math.sin(phi))) - e1*0.5*Math.log((1 + e1*Math.sin(phi))/(1 - e1*Math.sin(phi)));
R = C*Math.exp(-n*isolat);
       gamma = n*(lambda - lambda0);
x = xs + R*Math.sin(gamma);
y = ys - R*Math.cos(gamma);
//
xlambert2e = Math.round(x);
//UnitŽ en m
ylambert2e = Math.round(y);
//UnitŽ en m
//
return xlambert2e;
return ylambert2e;}function lambert2geo(x,y) {
//
//Parametres Lambert IIe
//
xs = 600000.000;
ys = 8199695.768;
n = 0.7289686274;
C = 11745793.39;
lambda0 = 0.04079234433;
//
//Parametres ellipsoide Clarke 1880
//
a1 = 6378249.2;
f1 = 1/293.466021;
b1 = a1*(1-f1);
e1 = Math.sqrt(Math.pow(a1,2) - Math.pow(b1,2))/a1;
//
//Parametres ellipsoide WGS84
//
a2 = 6378137.0;
f2 = 1/298.257223563;
b2 = a2*(1-f2);	e2 = Math.sqrt(Math.pow(a2,2) - Math.pow(b2,2))/a2;
//
//Parametre NTF
//
R = Math.sqrt(Math.pow((x - xs),2) + Math.pow((y - ys),2));
gamma = Math.atan((x - xs)/(ys - y));
lambda = lambda0 + gamma/n;
liso = -Math.log(R/C)/n;
//
phi0 = 2*Math.atan(Math.exp(liso)) - Math.PI/2;	test = 1;
eps = 0.00000000001;
//
while (test > eps) {
	phi1 = 2*Math.atan(Math.pow((1 + e1*Math.sin(phi0))/(1 - e1*Math.sin(phi0)),e1/2)*Math.exp(liso)) - Math.PI/2;
	test = Math.abs((phi1 - phi0)/phi0);
	phi0 = phi1;
}
phi = phi0;	h = 0;
W = Math.sqrt(1 - Math.pow(e1,2)*Math.pow(Math.sin(phi),2));
N = a1/W;
X = (N + h)*Math.cos(phi)*Math.cos(lambda);	Y = (N + h)*Math.cos(phi)*Math.sin(lambda);	Z = (N*(1 - Math.pow(e1,2)) + h)*Math.sin(phi);
X1 = X - 168;
Y1 = Y - 60;
Z1 = Z + 320;
lambda = 2*Math.atan(Y1/(X1 + Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2))));	phi0 = Math.atan(Z1/(Math.pow(X1,2) + Math.pow(Y1,2)*(1 - a2*Math.pow(e2,2)/Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2) + Math.pow(Z1,2)))));
test1 = 1;
while (test1 > eps) {
	phi1 = Math.atan((Z1/Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2)))/(1 - a2*Math.pow(e2,2)*Math.cos(phi0)/Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2))/Math.sqrt(1 - Math.pow(e2,2)*Math.pow(Math.sin(phi0),2))));
	test1 = Math.abs((phi1 - phi0)/phi0);
	phi0 = phi1;
}
phi = phi0;
h = Math.sqrt(Math.pow(X1,2) + Math.pow(Y1,2))/Math.cos(phi) - a2/Math.sqrt(1 - Math.pow(e2,2)*Math.pow(Math.sin(phi),2));
latitude = phi;
longitude = lambda;

return latitude;
return longitude;}

function geo2utm(longitude,latitude) {
//
//Parametres ellipsoide WGS84
a1 = 6378137.0;
f1 = 1/298.257223563;
b1 = a1*(1-f1);
e1 = Math.sqrt(Math.pow(a1,2) - Math.pow(b1,2))/a1;
//
//Calcul du fuseau UTM
fuseau = Math.floor((longitude*180/Math.PI + 180)/6)+1;
longitude0 = (6*fuseau - 183)*Math.PI/180;
//Longitude origine
latitude0 = 0;
//Latitude origine
k0 = 0.9996;
//Facteur d'Žchelle
X0 = 500000;
//False easting
Y0 = 0;
//False northing, hŽmisphre nord
//fuseau_out = [num2str(fuseau),'x'];
n = k0*a1;
//Rayon de la sphre intermŽdiaire
longitude_c = longitude0;
//Longitude origine par rapport au mŽridien origine
//
//Coefficients pour l'arc de mŽridien
xe1 = Math.pow(e1,0);
xe2 = Math.pow(e1,2);
xe3 = Math.pow(e1,4);
xe4 = Math.pow(e1,6);
xe5 = Math.pow(e1,8);
C1 = xe1 - xe2/4 - 3*xe3/64 - 5*xe4/256 - 175*xe5/16384;
C2 = -3*xe2/8 - 3*xe3/32 - 45*xe4/1024 - 105*xe5/4096;
C3 = 15*xe3/256 + 45*xe4/1024 + 525*xe5/16384;
C4 = -35*xe4/3072 - 175*xe5/12288;
C5 = 315*xe5/131072;
betastar = C1*latitude + C1*Math.sin(2*latitude) + C3*Math.sin(2*2*latitude) + C4*Math.sin(2*3*latitude) + C5*Math.sin(2*4*latitude);
//
//Constantes sur X,Y
Xs = X0;
Ys = Y0;
//
//Coefficients de projection
C1 = xe1 - xe2/4 - 3*xe3/64 - 5*xe4/256 - 175*xe5/16384;
C2 = xe2/8 - 1*xe3/96 - 9*xe4/1024 - 901*xe5/184320;
C3 = 13*xe3/768 + 17*xe4/5120 - 311*xe5/737280;
C4 = 61*xe4/15360 + 899*xe5/430080;
C5 = 49561*xe5/41287680;
//
//Latitude isomŽtrique
latiso = Math.log(Math.tan(Math.PI/4 + latitude/2)*Math.pow((1 - e1*Math.sin(latitude))/(1 + e1*Math.sin(latitude)),e1/2));	phi = Math.asin(Math.sin(longitude - longitude_c)/cosh(latiso));
latiso_s = Math.log(Math.tan(Math.PI/4 + phi/2));
eta = Math.atan(sinh(latiso)/Math.cos(longitude - longitude_c));
xutm = Math.round(n*C1*latiso_s + n*(C2*Math.cos(2*eta)*sinh(2*latiso_s) + C3*Math.cos(2*2*eta)*sinh(2*2*latiso_s) + C4*Math.cos(2*3*eta)*sinh(2*3*latiso_s) + C5*Math.cos(2*4*eta)*sinh(2*4*latiso_s)) + Xs);
yutm = Math.round(n*C1*eta + n*(C2*Math.sin(2*eta)*cosh(2*latiso_s) + C3*Math.sin(2*2*eta)*cosh(2*2*latiso_s) + C4*Math.sin(2*3*eta)*cosh(2*3*latiso_s) + C5*Math.sin(2*4*eta)*cosh(2*4*latiso_s)) + Ys);
return xutm;
return yutm;}

function cosh(xx) {return (Math.exp(xx) + Math.exp(-xx))/2;}function sinh(xx) {  return (Math.exp(xx) - Math.exp(-xx))/2;}
