function [drabf,drbcf,dracf,thetaf,vrabf,vrbcf,vracf,tf,arab,arbc,arac,ktot] = lepnorm(drab,drbc,theta,Frab,Frbc,Frac,vrabi,vrbci,vraci,hr1r1,hr1r2,hr1r3,hr2r2,hr2r3,hr3r3,ma,mb,mc,ti,dt,MEP)

      hesm = [hr1r1 hr1r2 hr1r3; hr1r2 hr2r2 hr2r3; hr1r3 hr2r3 hr3r3];
      
      % Reduced masses.
      mab = (ma*mb)/(ma+mb);
      mbc = (mb*mc)/(mb+mc);
      mac = (ma*mc)/(ma+mc);

      prab = vrabi*mab;
      prbc = vrbci*mbc;
      prac = vraci*mac;
      
      % G-Matrix
      y1 = 1/ma;
      y2 = 1/mb;
      y3 = 1/mc;

      GM(1,1)=y1+y2;
      GM(1,2)=y2*cos(theta);
      GM(2,1)=GM(1,2);
      GM(2,2)=y3+y2;
      GM(1,3)=-y2*sin(theta)/drbc;
      GM(3,1)=GM(1,3);
      GM(2,3)=-y2*sin(theta)/drab;
      GM(3,2)=GM(2,3);
      gmt=1/(drab.^2)+1/(drbc.^2)-(2*cos(theta)/(drab*drbc));
      GM(3,3)=y1/(drab.^2)+y3/(drbc.^2)+y2*gmt;
      
      [GMVec,GMVal] = eig(GM);
      GMVal1 = GMVal.^0.5;
      GMVal2 = diag(diag(GMVal).^-0.5);
      GRR = GMVec*GMVal1*GMVec';
      GROOT = GMVec*GMVal2*GMVec';

      % G-Matrix Weighted Hessian;
      MWH = GRR*hesm*GRR;
      [ALT,W2] = eig(MWH); %ALT is antisymmetric version in Fort code but that does not give the right G-Matrix!!!!
      W2 = diag(W2);
      
      % Gradient Vector in mass-weighted coordinates
      GRAD = [-Frab ; -Frbc ; -Frac];
      GRADN = ALT'*(GRR*GRAD);
      
      % Momentum Vector in Normal Coordinates
      MOM = [prab;prbc;prac];
      PCMO = ALT'*(GRR*MOM);
      
      ktot = 0.5*(PCMO(1).^2+PCMO(2).^2+PCMO(3).^2);
      
      for i = 1:3;
          if W2(i) < -1.0E-15;
              wmod=sqrt(abs(W2(i)));
              wmt=wmod*dt;
              q(i)=PCMO(i)*sinh(wmt)/wmod+GRADN(i)*(1-cosh(wmt))/(wmod.^2);
              PCMO(i) = PCMO(i)*cosh(wmt)-GRADN(i)*sinh(wmt)/wmod;
          elseif abs(W2(i)) < 1.0E-15;
              q(i)=PCMO(i)*dt-(0.5*GRADN(i)*(dt.^2));
              PCMO(i) = PCMO(i)-GRADN(i)*dt;
          else
              wroot = sqrt(W2(i));
              tfn1 = GRADN(i)*(1-cos(wroot*dt))/(wroot.^2);
              q(i) = PCMO(i)*sin(wroot*dt)/wroot-tfn1;
              PCMO(i) = PCMO(i)*cos(wroot*dt)-GRADN(i)*sin(wroot*dt)/wroot;
              
          end
      end
      
      XX = GRR*(ALT*q');
      
      if MEP == true;
        XX = XX.*5;
      end
      
      drabf = drab + XX(1);
      drbcf = drbc + XX(2);
      thetaf = theta + XX(3);
      dracf = ((drab.^2)+(drbc.^2)-2*drab*drbc*cos(thetaf)).^0.5;
      
      MOM = GROOT*ALT*PCMO;
      
      tf = ti + dt;

      vrabf = MOM(1)/mab;
      vrbcf = MOM(2)/mbc;
      vracf = 0;
%       if (vrabf+vrbcf < 0)
%           vracf = -((vrabf.^2)+(vrbcf.^2)-2*vrabf*vrbcf*cos(thetaf)).^0.5;
%       else
%           vracf = ((vrabf.^2)+(vrbcf.^2)-2*vrabf*vrbcf*cos(thetaf)).^0.5;
%       end
      
      arab = Frab/mab;
      arbc = Frbc/mbc;
      if (arab+arbc < 0)
          arac = -((arab.^2)+(arbc.^2)-2*arab*arbc*cos(thetaf)).^0.5;
      else
          arac = ((arab.^2)+(arbc.^2)-2*arab*arbc*cos(thetaf)).^0.5;
      end