Le portail francophone de la géomatique

        const double C = 11754255.426096; //constante de la projection
        const double E = 0.0818191910428158; //première excentricité de l'ellipsoïde
        const double N = 0.725607765053267; //exposant de la projection
        const double Xs = 700000.0; //coordonnées en projection du pole
        const double Ys = 12655612.049876; //coordonnées en projection du pole
        const double L0 = 3; // longitude de référence

        public void ConvertLatLongToLambert93(double lat, double lon, ref double px, ref double py)
        {
            double lat_rad = lat / 180 * Math.PI;
            double lat_iso = Math.Atanh(Math.Sin(lat_rad)) - E * Math.Atanh(E * Math.Sin(lat_rad)); //latitude isométrique
            px = ((C * Math.Exp(-N * (lat_iso))) * Math.Sin(N * (lon - L0) / 180 * Math.PI) + Xs);
            py = (Ys - (C * Math.Exp(-N * (lat_iso))) * Math.Cos(N * (lon - L0) / 180 * Math.PI));
        }

        public void ConvertLambert93ToLatLong(double px, double py, ref double lat, ref double lon)
        {
            double x = px - Xs;
            double y = py - Ys;
            double latiso = Math.Log(C / Math.Sqrt((x * x) + (y * y))) / N;
            double sinPhi = Math.Tanh(latiso + E * Math.Atanh(E * Math.Sin(1)));
            for (int i = 0; i < 6; i++) sinPhi = Math.Tanh(latiso + E * Math.Atanh(E * sinPhi));
            lat = Math.Asin(sinPhi) * 180.0 / Math.PI;
            lon = ((Math.Atan(-x / y)) / N + L0 / 180.0 * Math.PI) / Math.PI * 180.0;
        }