diff --git a/tracy/tracy/src/t2elem.cc b/tracy/tracy/src/t2elem.cc
index c38f5ca20b0e831d29470b65670dba51bedbee9d..0a31571745a5710a0758ac908c767c0fc0e2adb9 100644
--- a/tracy/tracy/src/t2elem.cc
+++ b/tracy/tracy/src/t2elem.cc
@@ -81,8 +81,14 @@ void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R,
 
  Purpose:
  Get the longitudinal momentum ps
-**********************************************************/
-
+ 
+    for approximation Hamitonian:
+    ps = 1+delta   
+    
+    For exact hamitonian:
+    
+    ps 
+ **********************************************************/
 template<typename T>
 inline T get_p_s(ss_vect<T> &x)
 {
@@ -588,27 +594,35 @@ void radiate(ss_vect<T> &x, const double L, const double h_ref, const T B[]) {
         is_tps<T>::emittance(B2_perp, ds, ps0, xp);
 }
 
+/********************************************************************
 static double get_psi(double irho, double phi, double gap) {
-    /* Correction for magnet gap (longitudinal fringe field)
+
+  Purpose:
+    Correction for magnet gap (longitudinal fringe field)
 
      irho h = 1/rho [1/m]
      phi  dipole edge angle
      gap  full gap between poles
 
-     2
-     K1*gap*h*(1 + sin phi)
-     psi = ----------------------- * (1 - K2*g*gap*tan phi)
-     cos phi
+                             2
+            K1*gap*h*(1 + sin phi)
+     psi = ----------------------- * (1 - K2*h*gap*tan phi)
+                  cos phi
 
      K1 is usually 1/2
-     K2 is zero here                                                  */
+     K2 is zero here                                                  
+
 
+
+*********************************************************************/
+static double get_psi(double irho, double phi, double gap) {
+   
     double psi;
 
     const double k1 = 0.5, k2 = 0.0;
 
-    if (phi == 0.0)
-        psi = 0.0;
+    if (phi == 0.0)  //hard edge
+        psi = 0.0;  
     else
         psi = k1 * gap * irho * (1.0 + sqr(sin(dtor(phi)))) / cos(dtor(phi))
                 * (1.0 - k2 * gap * irho * tan(dtor(phi)));
@@ -771,6 +785,21 @@ static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x) {
     }
 }
 
+
+/***************************************************************
+template<typename T>
+void p_rot(double phi, ss_vect<T> &x)
+
+  Purpose:
+     the effect of edge focusing of the dipole
+  Input:
+     phi:  entrance or exit angle of the dipole edge
+     x:    initial cooridinate of the particle
+  Output:
+  
+  Return:
+  
+****************************************************************/
 template<typename T>
 void p_rot(double phi, ss_vect<T> &x) {
     T c, s, ps, p;
@@ -787,6 +816,19 @@ void p_rot(double phi, ss_vect<T> &x) {
     x[ct_] += (1.0 + x1[delta_]) * x1[x_] * s / p;
 }
 
+/***************************************************************
+template<typename T>
+void bend_fringe(double hb, ss_vect<T> &x)
+
+  Purpose:
+     the effect of longitudinal fringe field using exact Hamitonian 
+  Input:
+  
+  Output:
+  
+  Return:
+  
+****************************************************************/
 template<typename T>
 void bend_fringe(double hb, ss_vect<T> &x) {
     T coeff, u, ps, ps2, ps3;
@@ -822,6 +864,18 @@ void bend_fringe(double hb, ss_vect<T> &x) {
  form. If keep up to the 4-th order nonlinear terms, the formula with goes to the
  following:
  (see E. Forest's book (Beam Dynamics: A New Attitude and Framework), p390):
+                b2
+ x = x0 +/- ---------- (x0^3 + 3*z0^2*x0)
+            12(1 + dP)
+
+                   b2
+ px = px0 +/-  ---------- (2*x0*y0*py0 - x0^2*px0 - y0^2*py0)
+               4(1 + dP)
+
+                b2
+ y = y0 -/+ ---------- (y0^3 + 3*x0^2*y0)
+            12(1 + dP)
+
                      b2
        x = x0 +/- ---------- (x0^3 + 3*z0^2*x0)
                   12(1 + dP)
@@ -845,6 +899,16 @@ void bend_fringe(double hb, ss_vect<T> &x) {
        tau = tau0 -/+ ----------- (y0^3*px - x0^3*py + 3*x0^2*y*py - 3*y0^2*x0*px)
                       12(1 + dP)^2
 
+                 b2
+ py = py0 -/+  ---------- (2*x0*y0*px0 - y0^2*py0 - x0^2*py0)
+               4(1 + dP)
+
+ dP = dP0;
+
+
+                  b2
+ tau = tau0 -/+ ----------- (y0^3*px - x0^3*py + 3*x0^2*y*py - 3*y0^2*x0*px)
+                12(1 + dP)^2
 
  Input:
  b2       = quadrupole strength
@@ -1920,6 +1984,32 @@ void Drift_Pass_M(CellType &Cell, Vector &xref, Matrix &X) {
     Drift(Cell.Elem.PL, xref);
 }
 
+/****************************************************************************/
+/* void thinkick_M(long Order, double *MB, double L, double irho, pthicktype pthick,
+                double *xref, vector *x)
+
+   Purpose:
+
+
+   Input:
+       none
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       thinkick
+
+   Comments:
+       none
+
+****************************************************************************/
 void thin_kick_M(int Order, double MB[], double L, double irho, Vector &xref,
         Matrix &x) {
     int i;
@@ -2374,6 +2464,31 @@ void Mpole_Setmatrix(int Fnum1, int Knum1, double K) {
     bendmat(M->AD55, elemp->PL / 2.0, M->Pirho, 0.0, M->PTx2, -M->Pgap, K);
 }
 
+/****************************************************************************/
+/* void Wiggler_Setmatrix(long Fnum1, long Knum1, double L, double lambda, double k0, double kx)
+
+   Purpose:
+
+
+   Input:
+       none
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+
+****************************************************************************/
 void Wiggler_Setmatrix(int Fnum1, int Knum1, double L, double kx, double kz,
         double k0) {
     double t, s, c, k, ky, LL;
@@ -2407,6 +2522,33 @@ void Wiggler_Setmatrix(int Fnum1, int Knum1, double L, double kx, double kz,
     mergeto4by5(W->W55, ah, av);
 }
 
+/****************************************************************************/
+/* void Mpole_Pass_M(CellType *Cell, double *xref, vector *x)
+
+   Purpose:
+
+
+   Input:
+       none
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       BUG does nothing if quad .... !!!! Meth_first not see
+       maybe just put a globval.MatMeth
+       30/03/03 : changed   Meth_First by Meth_Fourth
+
+****************************************************************************/
 void Mpole_Pass_M(CellType &Cell, Vector &xref, Matrix &x) {
     double k;
     elemtype *elemp;
@@ -2453,6 +2595,31 @@ void Mpole_Pass_M(CellType &Cell, Vector &xref, Matrix &x) {
     LtoG(xref, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1);
 }
 
+/****************************************************************************/
+/* void Wiggler_Pass_M(CellType *Cell, double *xref, vector *x)
+
+   Purpose:
+
+
+   Input:
+       none
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+
+****************************************************************************/
 void Wiggler_Pass_M(CellType &Cell, Vector &xref, Matrix &x) {
     elemtype *elemp;
     WigglerType *W;