<?php
function supprDecimales($dec) {
    if ($dec>0 )
        return floor($dec);
        else
            return ceil($dec);
}

function dec2dms($dec) {
    $d = supprDecimales($dec);
    $m = supprDecimales(($dec - $d) * 60);
    $s = abs(round(((($dec - $d) * 60) - $m) * 60));
    $m = abs($m);
    return $d.'°'.$m."'".$s.'"';
}

/**/
function ALG0001($phi,$e) {
    $temp = ( 1 - ( $e * sin( $phi ) ) ) / ( 1 + ( $e * sin( $phi ) ) );
    $L = log ( tan ( (pi()/4) + ($phi/2) ) * pow ($temp, ($e/2) ));
    return $L;
}

/** Calcul de la latitude à partir de la latitude isométrique
 **/
function ALG0002($L,$e,$epsilon) {
    $phi[0] = 2 * atan(exp($L)) - (pi()/2);

    $i=0;
    do {
        $i++;
        $temp = ( 1 + ( $e * sin( $phi[$i-1] ) ) ) / ( 1 - ( $e * sin( $phi[$i-1] ) ) );
        $phi[$i] = 2 * atan ( pow ($temp, ($e/2)) * exp ($L) ) - pi()/2;
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);

    return $phi[$i];
}

/** Transformation de coordonnées en projection conique conforme de Lambert, en coordonnées géographiques
 ** @param double $X		Coordonnée X en projection conique conforme de Lambert du point
 ** @param double $Y		Coordonnée Y en projection conique conforme de Lambert du point
 ** @param double $n		Exposant de la projection
 ** @param double $c		Constante de la projection
 ** @param double $Xs		Coordonnée Xs en projection du pôle
 ** @param double $Ys		Coordonnée Ys en projection du pôle
 ** @param double $lambdac	Longitude de l'origine par rapport au méridien origine
 ** @param double $e		Première excentricité de l'ellipsoïde
 ** @param double $epsilon	Tolérance de convergence
 ** @return array	lambda	Longitude par rapport au méridien origine
 **					phi		Latitude
 **/
function ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon) {
    $R = sqrt( pow(($X - $Xs),2) + pow(($Y - $Ys),2) );
    $gamma = atan(($X - $Xs)/($Ys - $Y));

    $lambda = $lambdac + ($gamma / $n);

    $L = (-1 / $n) * log(abs($R/$c));

    $phi = ALG0002($L,$e,$epsilon);

    $coords['lambda']=$lambda;
    $coords['phi']=$phi;

    return $coords;
}

/** Transformation des coordonnées géographiques ellipsoïdales en coordonnées cartésiennes
 **
 **/
function ALG0009($lambda,$phi,$he,$a,$e) {
    $N = ALG0021($phi,$a,$e);

    $X = ($N + $he) * cos($phi) * cos($lambda);

    $Y = ($N + $he) * cos($phi) * sin($lambda);

    $Z = ($N * (1 - $e*$e)  + $he) * sin ($phi);

    $coords['X']=$X;
    $coords['Y']=$Y;
    $coords['Z']=$Z;

    return $coords;
}

/**/
function ALG0012($X,$Y,$Z,$a,$e,$epsilon) {
    $lambda = atan ($Y/$X);

    $P = sqrt($X*$X + $Y*$Y);
    $phi[0] = atan ($Z/ ($P * (1 - ( ($a*$e*$e)/sqrt($X*$X + $Y*$Y + $Z*$Z) ) ) ) );

    $i = 0;
    do
    {
        $i++;
        $temp =  pow((1 - ( $a * $e*$e * cos($phi[$i - 1] )/( $P * sqrt(1 - $e*$e * sin($phi[$i - 1])*sin($phi[$i - 1]))))),-1);
        $phi[$i] = atan( $temp * $Z / $P );
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);

    $phix = $phi[$i];

    $he = ($P/cos($phix)) - ($a/sqrt(1 - $e*$e * sin($phix)*sin($phix)));

    $coords['lambda']=$lambda;
    $coords['phi']=$phix;
    $coords['he']=$he;

    return $coords;
}

/** Transformation d'un jeu de 7 paramètres entre 2 systèmes géodésiques
 ** @param	double $Tx	Translation suivant l'axe des X (de 1 vers 2)
 ** @param	double $Ty	Translation suivant l'axe des Y (de 1 vers 2)
 ** @param	double $Tz	Translation suivant l'axe des Z (de 1 vers 2)
 ** @param	double $D	Facteur d'échelle de 1 vers 2
 ** @param	double $Rx	Angle de rotation autour de l'axe des X en radian (de 1 vers 2)
 ** @param	double $Ry	Angle de rotation autour de l'axe des Y en radian (de 1 vers 2)
 ** @param	double $Rz	Angle de rotation autour de l'axe des Z en radian (de 1 vers 2)
 ** @param	double $U	Vecteur de coordonnées cartésiennes tridimensionnelles dans le système 1 : U=(Ux,Uy,Uz)
 ** @return array Vecteur de coordonnées cartésiennes tridimensionnelles dans le système 2, V=(Vx,Vy,Vz)
 */
function ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$U) {
    $V=array();
    $V['X'] = $Tx + $U['X'] * (1 + $D) + $U['Z'] * $Ry - $U['Y'] * $Rz;
    $V['Y'] = $Ty + $U['Y'] * (1 + $D) + $U['X'] * $Rz - $U['Z'] * $Rx;
    $V['Z'] = $Tz + $U['Z'] * (1 + $D) + $U['Y'] * $Rx - $U['X'] * $Ry;
    return $V;
}

/** Détermination des paramètres de calcul d'une projection Lambert conique
 ** dans le cas tangent, avec ou sans facteur d'échelle en fonction des paramètres de définition usuels
 ** @return array n, C, lambdac, Xs, Ys
 **/
function ALG0019($lambda0,$phi0,$k0,$X0,$Y0,$a,$e) {
    $lambdac = $lambda0;
    $n = sin($phi0);
    $C = $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) * exp ( $n * ALG0001($phi0,$e) );
    $Xs = $X0;
    $Ys = $Y0 + $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) ;

    $tab ['e'] = $e;
    $tab ['n'] = $n;
    $tab ['C'] = $C;
    $tab ['lambdac'] = $lambdac;
    $tab ['Xs'] = $Xs;
    $tab ['Ys'] = $Ys;

    return $tab;

}

/** Calcul de la grande normale de l'ellipsoïde
 ** @param  double $phi		Latitude
 ** @param  double $a		Demi-grand axe de l'ellipsoïde
 ** @param  double $e		Première excentricité de l'ellipsoïde
 ** @return double
 **/
function ALG0021($phi,$a,$e) {
    $N = $a/sqrt( 1 - $e * $e * sin($phi) * sin($phi) );
    return $N;
}

/** Calcul des constantes d'une projection Lambert conique conforme dans le cas sécant
 ** @param	double $lambda0	Longitude origine en radian par rapport au méridien origine
 ** @param	double $phi0	Latitude origine
 ** @param	double $X0		Coordonnée X en projection du point origine
 ** @param	double $Y0		Coordonnée Y en projection du point origine
 ** @param	double $phi1	Latitude en radian du 1er parallèle automécoïque
 ** @param	double $phi2	Latitude en radian du 2ème parallèle automécoïque
 ** @param  double $a		Demi-grand axe de l'ellipsoïde
 ** @param	double $e		Première excentricité de l'ellipsoïde
 ** @return array n, C, lambdac, Xs, Ys
 **/
function ALG0054($lambda0,$phi0,$X0,$Y0,$phi1,$phi2,$a,$e) {
    $lambdac = $lambda0;
    $n = ( (log( (ALG0021($phi2,$a,$e)*cos($phi2))/(ALG0021($phi1,$a,$e)*cos($phi1)) ))/(ALG0001($phi1,$e) - ALG0001($phi2,$e) ));
    $C = ((ALG0021($phi1,$a,$e)* cos($phi1))/$n) * exp($n * ALG0001($phi1,$e));

    if ($phi0 == (pi()/2)) {
        $Xs = $X0;
        $Ys = $Y0;
    } else {
        $Xs = $X0;
        $Ys = $Y0 + $C * exp(-1 * $n * ALG0001($phi0,$e));
    }
    $tab=array(	'e' => $e,
        'n' => $n,
        'C' => $C,
        'lambdac' => $lambdac,
        'Xs'=> $Xs,
        'Ys'=> $Ys);
    return $tab;
}

function Lambert2WGS84($X,$Y,$orig='L93') {
    $epsilon = 0.00000000001;

    $lambdac = 0.04079234433;	// pour greenwich
    $e = 0.08248325676; 		// première excentricité de l ellipsoïde Clarke 1880 français
    $he = 100;
    $a = 6378249.2;				// demi-grand axe de l'ellipsoide

    $Tx = -168;
    $Ty = -60;
    $Tz = +320;
    $D  = 0;
    $Rx = $Ry = $Rz = 0;

    $orig=strtoupper($orig);
    switch ($orig) {
        case 'LI':
        case 'L1':		$n = 0.7604059656;
        $c = 11603796.98;
        $Xs = 600000;
        $Ys = 5657616.674;
        break;
        case 'LII':
        case 'LIIE':
        case 'L2E':
        case 'L2':		$n = 0.7289686274;
        $c = 11745793.39;
        $Xs = 600000;
        if ($orig=='L2E' || $orig=='LIIE')
            $Ys = 8199695.768;
            else	$Ys = 6199695.768;
            break;
        case 'LIII':
        case 'L3':		$n = 0.6959127966;
        $c = 11947992.52;
        $Xs = 600000;
        $Ys = 6791905.085;
        break;
        case 'LIV':
        case 'L4':		$n = 0.6712679322;
        $c = 12136281.99;
        $Xs = 234.358;
        $Ys = 7239161.542;
        break;
        case 'L93':
        default:		$n = 0.7256077650;
        $c = 11745255.426;
        $Xs = 700000;
        $Ys = 12655612.050;
        break;
    }

    $coords = ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon);

    $coords = ALG0009($coords['lambda'],$coords['phi'],$he,$a,$e);

    $coords = ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$coords);

    $a = 6378137.0; // ellipsoïdes WGS84
    $f = 1/298.257223563;
    $b = $a*(1-$f);
    $e = sqrt(($a*$a - $b*$b)/($a*$a));

    $X = $coords['X'];
    $Y = $coords['Y'];
    $Z = $coords['Z'];
    $coords = ALG0012($X,$Y,$Z,$a,$e,$epsilon);

    $xy['long']  = rad2deg($coords['lambda']);
    $xy['lat']   = rad2deg($coords['phi']);
    return $xy;
}

/**
 * Function to convert geographic coordinates to "lambert 2 etendue" coordinates
 * from: http://www.forumsig.org/showthread.php?p=64050#post64050
 **/
function geos2lambert($latitude,$longitude) {
    //0)degres-minutes-secondes + orientation (d,m,s,o) en radian
    $lambda_w = $longitude * pi()/180 ;
    $phi_w =  $latitude * pi()/180 ;

    //1) coordonnées géographiques WGS84 (phi_w,lambda_w) en coordonnées cartésiennes WGS84 (X_w,Y_w,Z_w)
    $a_w = 6378137.0;
    $b_w = 6356752.314;

    $e2_w = ($a_w*$a_w-$b_w*$b_w)/($a_w*$a_w);

    //et on a la grande normale de l'ellipsoide WGS84
    $N = $a_w/sqrt(1-$e2_w*pow(sin($phi_w),2));

    //ainsi on obtient
    $X_w = $N * cos($phi_w) * cos($lambda_w);
    $Y_w = $N * cos($phi_w) * sin($lambda_w);
    $Z_w = $N * (1-$e2_w) * sin($phi_w);
    //Ref.: http://www.ign.fr/telechargement/MPr...RCE/NTG_80.pdf and http://de.wikipedia.org/wiki/WGS84

    //2) coordonnées cartésiennes WGS84 (X_w,Y_w,Z_w) en coordonnées cartésiennes NTF (X_n,Y_n,Z_n)

    $dX = 168.0;
    $dY = 60.0;
    $dZ = -320.0;

    $X_n = $X_w + $dX;
    $Y_n = $Y_w + $dY;
    $Z_n = $Z_w + $dZ;
    //ref.: http://support.esrifrance.fr/Documents/Generalites/Projections/Generalites/Generalites.htm#2

    //3) coordonnées cartésiennes NTF (X_n,Y_n,Z_n) en coordonnées géographiques NTF (phi_n,lambda_n)
    $a_n = 6378249.2;
    $b_n = 6356515.0;

    $e2_n = ($a_n*$a_n-$b_n*$b_n)/($a_n*$a_n);
    //on définit une tolérance de convergence
    $epsilon = pow(10,-10);

    //puis on amorce une boucle de calcul
    $p0=atan($Z_n/sqrt($X_n*$X_n+$Y_n*$Y_n)*(1-($a_n*$e2_n)/(sqrt($X_n*$X_n+$Y_n*$Y_n+$Z_n*$Z_n))));
    $p1=atan(($Z_n/sqrt($X_n*$X_n+$Y_n*$Y_n))/(1-($a_n*$e2_n*cos($p0))/(sqrt(($X_n*$X_n+$Y_n*$Y_n)*(1-$e2_n*pow(sin($p0),2))))));
    while(!(abs($p1-$p0)<$epsilon)){
        $p0 = $p1;
        $p1 = atan(($Z_n/sqrt($X_n*$X_n+$Y_n*$Y_n))/(1-($a_n*$e2_n*cos($p0))/(sqrt(($X_n*$X_n+$Y_n*$Y_n)*(1-$e2_n*pow(sin($p0),2))))));
    }
    $phi_n = $p1;
    $lambda_n = atan($Y_n/$X_n);

    //4) coordonnées géographiques NTF (phi_n,lambda_n) en coordonnées projetées en Lambert II étendu (X_l2e, Y_l2e)
    $n = 0.7289686274;
    $c = 11745793.39;
    $Xs = 600000.0;
    $Ys = 8199695.768;

    $e_n = sqrt($e2_n);
    $lambda0 = 0.04079234433198;   //correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich
    //puis on calcule la latitude isométrique
    $L = log(tan(pi()/4 + $phi_n/2) * pow(((1-$e_n*sin($phi_n))/(1+$e_n*sin($phi_n))),($e_n/2)));

    //enfin on projette
    $X_l2e = $Xs + $c*exp((-$n*$L))*sin($n*($lambda_n-$lambda0));
    $Y_l2e = $Ys - $c*exp((-$n*$L))*cos($n*($lambda_n-$lambda0));

    return array("x_l2e"=>$X_l2e,"y_l2e"=>$Y_l2e);
}


/** Conversion WGS84 en Lambert93
 ** NTG_71.pdf
 ** http://www.ign.fr/affiche_rubrique.asp?rbr_id=1700&lng_id=FR
 **/
function geos2lambert93($latitude,$longitude) {

    // Système WGS84
    $a=6378137; 		// demi grand axe de l'ellipsoide (m)
    $e=0.08181919106;	// première excentricitè de l'ellipsoide

    // Paramètres de projections
    $l0=$lc=deg2rad(3);		// longitude de référence
    $phi0=deg2rad(46.5);	// latitude d'origine en radian
    $phi1=deg2rad(44);		// 1er parallele automÃ©coÃ¯que
    $phi2=deg2rad(49);		// 2eme parallele automÃ©coÃ¯que

    $x0=700000;		//coordonnées à l'origine
    $y0=6600000;	//coordonnées à l'origine

    // coordonnées du point à traduire
    $phi=deg2rad($latitude);
    $l=deg2rad($longitude);

    // calcul des grandes normales
    $gN1=$a/sqrt(1-$e*$e*sin($phi1)*sin($phi1));
    $gN2=$a/sqrt(1-$e*$e*sin($phi2)*sin($phi2));

    //calculs de slatitudes isométriques
    $gl1=log(tan(pi()/4+$phi1/2)*pow((1-$e*sin($phi1))/(1+$e*sin($phi1)),$e/2));
    $gl2=log(tan(pi()/4+$phi2/2)*pow((1-$e*sin($phi2))/(1+$e*sin($phi2)),$e/2));
    $gl0=log(tan(pi()/4+$phi0/2)*pow((1-$e*sin($phi0))/(1+$e*sin($phi0)),$e/2));
    $gl=log(tan(pi()/4+$phi/2)*pow((1-$e*sin($phi))/(1+$e*sin($phi)),$e/2));

    //calcul de l'exposant de la projection
    $n=(log(($gN2*cos($phi2))/($gN1*cos($phi1))))/($gl1-$gl2);//ok

    //calcul de la constante de projection
    $c=(($gN1*cos($phi1))/$n)*exp($n*$gl1);//ok

    //calcul des coordonnées
    $ys=$y0+$c*exp(-1*$n*$gl0);

    //calcul des coordonnées lambert
    $x93=$x0+$c*exp(-1*$n*$gl)*sin($n*($l-$lc));
    $y93=$ys-$c*exp(-1*$n*$gl)*cos($n*($l-$lc));

    return array("x_93"=>$x93,"y_93"=>$y93);
}