function [LambertX,LambertY]=WGS84toLambert(opt,lg,lt)

%% Conversion données GPS ellipsoïde:
%Transformation standard: WGS84
%Transformation grille: RGF93
%opt (indication zone Lambert) : 'I', 'II', 'III', 'IV', 'IIE', '93'
%Coordonnées d'entrée en dégrés décimaux
%[lg]: vecteur colonne des longitudes
%[lt]: vecteur colonne des latitudes
%Remarque: il n'y a pas de prise en compte de l'altitude

%% Mise en forme des données vecteur colonne
[mlg,nlg]=size(lg);
[mlt,nlt]=size(lt);
if mlg<nlg
    lg=lg';
end
if mlt<nlt
    lt=lt';
end
[mlg,~]=size(lg);
[mlt,~]=size(lt);
if mlg~=mlt
    disp('lg et lt: Pas même nombre de données')
    return
end

%% altitude non prise en compte
Z=zeros(mlg,1);

%% Parallèles automécoïques
%ligne1: Lambert I   : 48.59852278	50.39591167
%ligne2: Lambert II  : 45.89891889	47.69601444
%ligne3: Lambert III : 43.19929139	44.99609389
%ligne4: Lambert IV  : 41.56038778	42.76766333
%ligne5: Lambert IIE : 45.89891889  47.69601444
%ligne6: Lambert 93  : 44.00000000  49.00000000
mat=[48.59852278	50.39591167
    45.89891889	47.69601444
    43.19929139	44.99609389
    41.56038778	42.76766333
    45.89891889 47.69601444
    44.00000000 49.00000000];

%% Constantes
pm=pi/180;

if strcmp(opt,'93')==1
    a=6378137;
    e=0.08181919112;
    f=1/298.257222101;
    lambda0dec=3;
else
    a=6378249.2; %demi grand axe de l'ellipsoide
    b=6356515.0; %petit demi axe
    e=((a*a-b*b)/(a*a))^0.5; %Première excentricité
    f=(a-b)/a;
    lambda0dec=2.33722916666667;%Origine par rapport au méridien de Greenwich: 2°20'14.025
end

e2=e^2;

%% Limites Lambert et Constantes des zones
switch opt
    case 'I'
        ind=1;
        %     n=0.7604059656;%exposant de la projection
        %     c=11603796.9767;
        Xs=600000;
        Ys=5657616.674;
    case 'II'
        ind=2;
        %     n=0.7289686274;%exposant de la projection
        %     c=11745793.39;
        Xs=600000;
        Ys=6199695.768;
    case 'III'
        ind=3;
        %     n=0.6959127966;%exposant de la projection
        %     c=11947992.52;
        Xs=600000;
        Ys=6791905.085;
    case 'IV'
        ind=4;
        %     n=0.6712679322;%exposant de la projection
        %     c=12136281.99;
        Xs=234.358;
        Ys=7239161.542;
    case 'IIE'
        ind=5;
        %     n=0.7289686274;%exposant de la projection
        %     c=11745793.39;
        Xs=600000;
        Ys=8199695.768;
    case '93'
        ind=6;
        %         n=0.725607765;%exposant de la projection
        %         c=11754255.426;
        Xs=700000;
        Ys=12655612.050;
end
mat=mat(ind,:);
if min(lt)<min(mat) || max(lt)>max(mat)
    msgbox(['Hors zone_' num2str(opt)])
end

%% En radian
mat=mat*pi/180;
esur2=e/2;
lambda0=pi*lambda0dec/180;%longitude origine(en radian)

%% Calcul de l'exposant "n" et de la constante "c" de projection
LatISO=zeros(1,2);
N=zeros(1,2);
for i=1:2
    esur2=e/2;
    esinphi=e*sin(mat(i));
    v=(1-esinphi)/(1+esinphi);
    v=exp((esur2)*log(v));
    LatISO(i)=log(tan(pi/4+mat(i)/2)*v);
    N(i)=(a/(1-e^2*(sin(mat(i))^2))^0.5);
end
LatISO=reshape(LatISO,1,2);
N=reshape(N,1,2);
n=log((N(1,2)*cos(mat(1,2)))/(N(1,1)*cos(mat(1,1))))/(LatISO(1,1)-LatISO(1,2));
c=((N(1,1)*cos(mat(1,1)))/n(1))*exp(n(1)*LatISO(1,1));

%% preallocation
phi=zeros(mlg,1);lambda=zeros(mlg,1);LatISOdata=zeros(mlg,1);
X=zeros(mlg,1);Y=zeros(mlg,1);
slat=zeros(mlg,1);clat=zeros(mlg,1);slon=zeros(mlg,1);clon=zeros(mlg,1);sslat=zeros(mlg,1);
rn=zeros(mlg,1);rm=zeros(mlg,1);
dlat=zeros(mlg,1);dlon=zeros(mlg,1);
T1=zeros(mlg,2);T2=zeros(mlg,2);T3=zeros(mlg,2);T4=zeros(mlg,2);
ind_lg_T1=zeros(mlg,111);ind_lg_T2=zeros(mlg,111);ind_lg_T3=zeros(mlg,111);ind_lg_T4=zeros(mlg,111);
n_line_T1=zeros(mlg,1);n_line_T2=zeros(mlg,1);n_line_T3=zeros(mlg,1);n_line_T4=zeros(mlg,1);
dX_vec=zeros(mlg,4);dY_vec=zeros(mlg,4);dZ_vec=zeros(mlg,4);
xgr3d3=zeros(mlg,1);ygr3d3=zeros(mlg,1);
dX=zeros(mlg,1);dY=zeros(mlg,1);dZ=zeros(mlg,1);

%Grille IGN chargée pas de 0.1°
pas_grille_RG3D3=0.1;

if strcmp(opt,'93')==1
    for k=1:mlg
        T1(k,:)=[(fix(lg(k)*10))/10 (fix(lt(k)*10))/10];
        if T1(k,1)==lg(k) || T1(k,2)==lt(k)
            T2(k,:)=T1(k,:);
            T3(k,:)=T1(k,:);
            T4(k,:)=T1(k,:);
        else
            T2(k,:)=[T1(k,1) T1(k,2)+pas_grille_RG3D3];
            T3(k,:)=[T1(k,1)+pas_grille_RG3D3 T1(k,2)];
            T4(k,:)=[T1(k,1)+pas_grille_RG3D3 T1(k,2)+pas_grille_RG3D3];
        end
        load GR3D3
        ind_lg_T1(k,:)=find(abs(GR3D3(:,1)-T1(k,1))<10^-4)'; %#ok<NODEF>
        n_line_T1(k)=find((abs((GR3D3(ind_lg_T1(k,:),2)-T1(k,2))))<10^-4)+(min(ind_lg_T1(k))-1);
       
        ind_lg_T2(k,:)=find(abs(GR3D3(:,1)-T2(k,1))<10^-4)';
        n_line_T2(k)=find((abs((GR3D3(ind_lg_T2(k,:),2)-T2(k,2))))<10^-4)+(min(ind_lg_T2(k))-1);

        ind_lg_T3(k,:)=find(abs(GR3D3(:,1)-T3(k,1))<10^-4)';
        n_line_T3(k)=find((abs((GR3D3(ind_lg_T3(k,:),2)-T3(k,2))))<10^-4)+(min(ind_lg_T3(k))-1);
        
        ind_lg_T4(k,:)=find(abs(GR3D3(:,1)-T4(k,1))<10^-4)';
        n_line_T4(k)=find((abs((GR3D3(ind_lg_T4(k,:),2)-T4(k,2))))<10^-4)+(min(ind_lg_T4(k))-1);
        
        dX_vec(k,:)=[GR3D3(n_line_T1(k),3) GR3D3(n_line_T2(k),3) GR3D3(n_line_T3(k),3) GR3D3(n_line_T4(k),3)];
        dY_vec(k,:)=[GR3D3(n_line_T1(k),4) GR3D3(n_line_T2(k),4) GR3D3(n_line_T3(k),4) GR3D3(n_line_T4(k),4)];
        dZ_vec(k,:)=[GR3D3(n_line_T1(k),5) GR3D3(n_line_T2(k),5) GR3D3(n_line_T3(k),5) GR3D3(n_line_T4(k),5)];
        
        %Interpolation bilénaire
        xgr3d3(k)=(lg(k)-T1(k,1))/pas_grille_RG3D3;
        ygr3d3(k)=(lt(k)-T1(k,2))/pas_grille_RG3D3;
        dX(k)=(1-xgr3d3(k))*(1-ygr3d3(k))*dX_vec(1)+(1-xgr3d3(k))*ygr3d3(k)*dX_vec(2)+xgr3d3(k)*(1-ygr3d3(k))*dX_vec(3)+xgr3d3(k)*ygr3d3(k)*dX_vec(4);
        dY(k)=(1-xgr3d3(k))*(1-ygr3d3(k))*dY_vec(1)+(1-xgr3d3(k))*ygr3d3(k)*dY_vec(2)+xgr3d3(k)*(1-ygr3d3(k))*dY_vec(3)+xgr3d3(k)*ygr3d3(k)*dY_vec(4);
        dZ(k)=(1-xgr3d3(k))*(1-ygr3d3(k))*dZ_vec(1)+(1-xgr3d3(k))*ygr3d3(k)*dZ_vec(2)+xgr3d3(k)*(1-ygr3d3(k))*dZ_vec(3)+xgr3d3(k)*ygr3d3(k)*dZ_vec(4);
    end
end
%% Changement de système geodesique WGS84 à NTF et calcul en Lambert
for i=1:mlg
    if strcmp(opt,'93')==0
        %Changement de système geodesique WGS84 NTF
        %Paramètres:
        da=6378137-a;
        df=((6378137-6356752.314245)/6378137)-f;
        if strcmp(opt,'93')==0
            dX=zeros(mlg,1)+(-168);
            dY=zeros(mlg,1)+(-60);
            dZ=zeros(mlg,1)+320;
        end
        pm=pi/180;
        slat(i)=sin(lt(i)*pm);
        clat(i)=cos(lt(i)*pm);
        slon(i)=sin(lg(i)*pm);
        clon(i)=cos(lg(i)*pm);
        sslat(i)=slat(i)^2;
        adb=1/(1-f);
        rn(i)=a/(1-e2*sslat(i))^0.5;
        rm(i)=a*(1-e2)/(1-e2*sslat(i))^1.5;
        Z(i)
        dlat(i)=(((((-dX(i)*slat(i)*clon(i)-dY(i)*slat(i)*slon(i))+dZ(i)*clat(i))+(da*((rn(i)*e2*slat(i)*clat(i))/a)))+(df*(rm(i)*adb+rn(i)/adb)*slat(i)*clat(i))))/(rm(i)+Z(i));
        dlon(i)=(-dX(i)*slon(i)+dY(i)*clon(i))/((rn(i)+Z(i))*clat(i));
        %Calcul des nouvelles coordonnées corrigées degrés 
        lg(i)=(lg(i)*pi/180-dlon(i))*180/pi;
        lt(i)=(lt(i)*pi/180-dlat(i))*180/pi;
    end
    
    phi(i)=lt(i)*pm;lambda(i)=lg(i)*pm; %Conversion en radian
    esinphi(i)=e*sin(phi(i));
    v(i)=(1-esinphi(i))/(1+esinphi(i));
    v(i)=exp((esur2)*log(v(i)));
    LatISOdata(i)=log(tan(pi/4+phi(i)/2)*v(i));
    X(i)=Xs+c*exp(-n*LatISOdata(i))*sin(n*(lambda(i)-lambda0));
    Y(i)=Ys-c*exp(-n*LatISOdata(i))*cos(n*(lambda(i)-lambda0));
end
LambertX=X;
LambertY=Y;
disp(['Coordonnées convertis en Lambert_' opt])
disp(num2str(LambertX,'X = %11.2f m'))
disp(num2str(LambertY,'Y = %11.2f m'))