diff --git a/Makefile_test b/Makefile_test
index abd565662..0a39d1bf3 100644
--- a/Makefile_test
+++ b/Makefile_test
@@ -51,7 +51,7 @@ test-track-7 test-track-8 test-track-9 test-track-10 test-track-11 test-track-12
 test-track-acd test-track-rotations \
 test-twiss test-twiss-2 test-twiss-3 test-twiss-4 test-twiss-5 \
 test-twiss-6 test-twiss-8 test-twiss-9 test-twiss-10 test-twiss-11 test-tapering \
-test-twiss-12 test-twiss-13 test-twiss-14 test-twiss-15 test-translation test-crabcavity \
+test-twiss-12 test-twiss-13 test-twiss-14 test-twiss-15 test-twiss-exact test-translation test-crabcavity \
 test-xrotation test-yrotation test-rotations test-interpolate test-rf-fringe \
 test-cororbit test-cororbit-2 test-cororbit-3 test-cororbit-4 \
 test-emit test-emit-2 \
diff --git a/doc/latexuguide/twiss.tex b/doc/latexuguide/twiss.tex
index b1d4e3ea2..c02d87825 100644
--- a/doc/latexuguide/twiss.tex
+++ b/doc/latexuguide/twiss.tex
@@ -30,7 +30,7 @@ \chapter{Twiss Module}
        \>SECTORPURE=logical,\\
        \>EIGENVECTOR=tabname, EIGENFILE=filename,\\
        \>KEEPORBIT=name, USEORBIT=name,\\
-       \>COUPLE=logical,\\
+       \>COUPLE=logical, exact=logical, \\
        \>RIPKEN=logical, TAPERING=logical \\;
 }
 
@@ -172,7 +172,10 @@ \chapter{Twiss Module}
   longer be set since TWISS in \madx is always calculated in
   coupled mode. \madx computes the coupled functions in the
   sense of Edwards and Teng \cite{edwards1973}. 
-  For the uncoupled cases they reduce to the C and S functions.    
+  For the uncoupled cases they reduce to the C and S functions. 
+
+  \ttitem{EXACT} If this is used the dirft is expanded around the actual
+  closed orbit instead of the ideal orbit. (Default: false)   
   
   \textit{ Twiss calculation is 4D only!} : The \texttt{TWISS}
   command calculates an approximate 6D closed orbit when the
diff --git a/scripts/expand_map.py b/scripts/expand_map.py
new file mode 100644
index 000000000..8000fafb9
--- /dev/null
+++ b/scripts/expand_map.py
@@ -0,0 +1,97 @@
+from sympy import *
+import numpy as np
+
+def make_fortran_str(in_exp):
+    in_str = str(in_exp)
+    for i in range(0,10):
+        print(in_str)
+        print(str(i))
+        in_str = in_str.replace(str(i),str(i)+'d0')
+    return in_str
+
+
+x,px,y,py,t,pt = symbols('x px y py t pt') # The canonical variables 
+xf,pxf,yf,pyf,tf,ptf = symbols('xf pxf yf pyf tf ptf') # The canonical variables 
+dl, beta = symbols('dl beta', positive = True, real = True)
+l_pz, myzero, csq = symbols('l_pz, myzero csq', real = True)
+
+
+dmapi = [x,py,y,py,t,pt]
+msq = 1 + 2*pt*beta + pt**2 - px**2 - py**2
+l_pz = dl / sqrt(msq)
+xf = px*l_pz
+pxf = px
+yf = py*l_pz
+pyf = py
+tf = (dl*beta - (beta + pt) * l_pz)
+ptf = pt
+myzero = 0
+
+dmapi = [x,px,y,py,t,pt]
+dmapf = [xf,pxf,yf,pyf,tf,ptf]
+te = np.empty((6,6,6),dtype=object)
+re = eye(6)
+
+re_out = open("re_out.txt", "w")
+te_out = open("te_out.txt", "w")
+
+for i in range(0,6):
+    for j in range(0,6):
+        re[i,j] = diff(dmapf[i],dmapi[j])
+        if(re[i,j]!=0):
+            m_i = i + 1
+            m_j = j + 1
+            re[i,j].simplify()
+            tmp = re[i,j]
+            tmp = tmp.subs(msq,csq)
+            print(f're({m_i},{m_j}) = '+ str(tmp), file=re_out)
+            for k in range(0,6):
+                te[i,j,k] = diff(re[i,j],dmapi[k])
+                m_k = k + 1
+                if(te[i,j,k]!=myzero and te[i,j,k] is not None):
+                    tmp = te[i,j,k]
+                    #tmp.simplify()
+                    tmp = tmp.subs(msq,csq)/2
+                    print(f'te({m_i},{m_j},{m_k}) = '+ make_fortran_str(tmp) , file=te_out)
+
+re_out.close()                    
+
+J = Matrix([[0, 1, 0, 0, 0, 0 ],
+           [-1, 0, 0, 0, 0, 0 ],
+           [0 , 0, 0, 1, 0, 0 ],
+           [0,  0,-1, 0, 0, 0 ],
+           [0,  0, 0, 0, 0, 1 ],
+           [0,  0, 0, 0,-1, 0 ]])
+#t1 = te[0,:,:]
+#f1 = t1*J*re + re.transpose()*J*t1
+#f1.simplify()
+
+res = re.transpose()*J*re
+res.simplify()
+exit()
+
+
+with open('toMADX', 'w') as f:
+    print('re(2,1) = '+ str(dpxdx), file=f)
+    print('re(2,3) = '+ str(dpxdy), file=f)
+    print('re(4,1) = '+ str(dpydx), file=f)
+    print('re(4,3) = '+ str(dpydy), file=f)
+
+
+fin = open("toMADX", "rt")
+#output file to write the result to
+fout = open("outToMADX.txt", "w")
+#for each line in the input file
+for line in fin:
+    #read replace the string and write to output file
+    a = line.replace('(L + Lint - Abs(L - Lint))','Leff')
+    a = a.replace('(x**2 + y**2)','R')
+    a = a.replace('(L + Lint)**2','Lsum')
+    a = a.replace('(-L + Lint)**2','Ldiff')
+    a = a.replace('4*x**2 + 4*y**2','4d0*R')
+
+    fout.write(a)
+
+#close input and output files
+fin.close()
+fout.close()
\ No newline at end of file
diff --git a/src/mad_dict.c b/src/mad_dict.c
index 566cb7ab7..d91f9a1f2 100644
--- a/src/mad_dict.c
+++ b/src/mad_dict.c
@@ -29,7 +29,7 @@ const char *const_constant_def =
 "const emass = 0.51099895000e-3; "
 "const mumass = 0.1056583755; "
 "const nmass = 0.93956542052; "
-"const pmass = 0.93827208816; "
+"const pmass = 0.93827208816; " 
 "const clight = 299792458; "
 "const qelect = 1.602176634e-19; "
 "const hbar = 6.582119569e-25; " /* GeV*s */
@@ -1076,6 +1076,7 @@ const char *const_command_def =
 "keeporbit= [s, default, default], "
 "tolerance= [r, 1.e-6], "
 "deltap   = [s, none], "
+"exact = [l, false, true], "
 "eigenvector = [l, false, true], "
 "eigenfile = [s, eigenvectors, eigenvectors], "
 "tapering = [l, false, true], "
diff --git a/src/twiss.f90 b/src/twiss.f90
index c7c6f3a58..e4a6fa895 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -40,6 +40,7 @@ SUBROUTINE twiss(rt,disp0,tab_name,sector_tab_name)
   fsecarb=.false.
 
   ORBIT0 = zero
+  exact_expansion = get_value('twiss ','exact ') .ne. zero
   call get_node_vector('orbit0 ', 6, orbit0)
   RT = EYE
   RW = EYE
@@ -6903,6 +6904,7 @@ end SUBROUTINE tmyrot
 
 SUBROUTINE tmdrf(fsec,ftrk,orbit,fmap,dl,ek,re,te)
   use twissbeamfi, only : beta, gamma, dtbyds
+  use twisslfi
   use matrices, only : EYE
   use math_constfi, only : zero, two, three
   implicit none
@@ -6924,34 +6926,90 @@ SUBROUTINE tmdrf(fsec,ftrk,orbit,fmap,dl,ek,re,te)
   logical :: fsec, ftrk, fmap
   double precision :: dl
   double precision :: orbit(6), ek(6), re(6,6), te(6,6,6)
+  double precision :: px, py, pt, csq, l_pz, c3sq, c52sq
 
   !---- Initialize.
   EK = zero
   RE = EYE
   if (fsec) TE = zero
   fmap = dl .ne. zero
-
   !---- First-order terms.
+  if(exact_expansion) then
+      px = orbit(2)
+      py = orbit(4)
+      pt = orbit(6)
+
+      csq = 1 + 2*pt*beta + pt**2 - px**2 - py**2
+      l_pz = dl / sqrt(csq)
+      c3sq = csq**(3d0/2d0)
+      c52sq =csq**(5d0/2d0)
+      re(1,2) = dl/sqrt(csq) + dl*px**2/c3sq;
+      re(1,4) = dl*px*py/c3sq;
+      re(1,6) = dl*px*(-beta - pt)/c3sq;
+      re(3,2) = dl*px*py/c3sq;
+      re(3,4) = dl/sqrt(csq) + dl*py**2/c3sq;
+      re(3,6) = dl*py*(-beta - pt)/c3sq;
+      re(5,2) = -dl*px*(beta + pt)/c3sq;
+      re(5,4) = -dl*py*(beta + pt)/c3sq;
+      re(5,6) = -dl/sqrt(csq) - dl*(-beta - pt)*(beta + pt)/c3sq;
+
+      if (fsec) then
+         te(1,2,2) = 3d0*dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**3d0/(2d0*c52sq)
+         te(1,2,4) = dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*py/(2d0*c52sq)
+         te(1,2,6) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + dl*px**2d0*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
+         te(1,4,2) = dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*py/(2d0*c52sq)
+         te(1,4,4) = dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px*py**2d0/(2d0*c52sq)
+         te(1,4,6) = dl*px*py*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
+         te(1,6,2) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*(-beta - pt)/(2d0*c52sq)
+         te(1,6,4) = 3d0*dl*px*py*(-beta - pt)/(2d0*c52sq)
+         te(1,6,6) = -dl*px/(2d0*csq**(3d0/2d0)) + dl*px*(-3d0*beta - 3d0*pt)*(-beta - pt)/(2d0*c52sq)
+         te(3,2,2) = dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*py/(2d0*c52sq)
+         te(3,2,4) = dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px*py**2d0/(2d0*c52sq)
+         te(3,2,6) = dl*px*py*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
+         te(3,4,2) = dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px*py**2d0/(2d0*c52sq)
+         te(3,4,4) = 3d0*dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*py**3d0/(2d0*c52sq)
+         te(3,4,6) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + dl*py**2d0*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
+         te(3,6,2) = 3d0*dl*px*py*(-beta - pt)/(2d0*c52sq)
+         te(3,6,4) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + 3d0*dl*py**2d0*(-beta - pt)/(2d0*c52sq)
+         te(3,6,6) = -dl*py/(2d0*csq**(3d0/2d0)) + dl*py*(-3d0*beta - 3d0*pt)*(-beta - pt)/(2d0*c52sq)
+         te(5,2,2) = -dl*(beta + pt)/(2d0*csq**(3d0/2d0)) - 3d0*dl*px**2d0*(beta + pt)/(2d0*c52sq)
+         te(5,2,4) = -3d0*dl*px*py*(beta + pt)/(2d0*c52sq)
+         te(5,2,6) = -dl*px/(2d0*csq**(3d0/2d0)) - dl*px*(-3d0*beta - 3d0*pt)*(beta + pt)/(2d0*c52sq)
+         te(5,4,2) = -3d0*dl*px*py*(beta + pt)/(2d0*c52sq)
+         te(5,4,4) = -dl*(beta + pt)/(2d0*csq**(3d0/2d0)) - 3d0*dl*py**2d0*(beta + pt)/(2d0*c52sq)
+         te(5,4,6) = -dl*py/(2d0*csq**(3d0/2d0)) - dl*py*(-3d0*beta - 3d0*pt)*(beta + pt)/(2d0*c52sq)
+         te(5,6,2) = -dl*px/(2d0*csq**(3d0/2d0)) - 3d0*dl*px*(-beta - pt)*(beta + pt)/(2d0*c52sq)
+         te(5,6,4) = -dl*py/(2d0*csq**(3d0/2d0)) - 3d0*dl*py*(-beta - pt)*(beta + pt)/(2d0*c52sq)
+         te(5,6,6) = -dl*(-beta - pt)/csq**(3d0/2d0) + dl*(beta + pt)/(2d0*csq**(3d0/2d0)) &
+         - dl*(-3d0*beta - 3d0*pt)*(-beta - pt)*(beta + pt)/(2d0*c52sq)
+      endif
 
-  re(1,2) = dl
-  re(3,4) = dl
-  re(5,6) = dl/(beta*gamma)**2
+      orbit(1) = orbit(1) + px*l_pz
+      orbit(3) = orbit(3) + py*l_pz
+      orbit(5) = orbit(5) + (dl*beta - (beta + pt) * l_pz)
+  else
 
-  ek(5) = dl*dtbyds
+     re(1,2) = dl
+     re(3,4) = dl
+     re(5,6) = dl/(beta*gamma)**2
+
+     ek(5) = dl*dtbyds
+
+     !---- Second-order terms.
+     if (fsec) then
+        te(1,2,6) = - dl / (two * beta)
+        te(1,6,2) = te(1,2,6)
+        te(3,4,6) = te(1,2,6)
+        te(3,6,4) = te(3,4,6)
+        te(5,2,2) = te(1,2,6)
+        te(5,4,4) = te(5,2,2)
+        te(5,6,6) = te(1,2,6) * three / (beta * gamma) ** 2
+     endif
+     
+     !---- Track orbit.
+     if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
 
-  !---- Second-order terms.
-  if (fsec) then
-     te(1,2,6) = - dl / (two * beta)
-     te(1,6,2) = te(1,2,6)
-     te(3,4,6) = te(1,2,6)
-     te(3,6,4) = te(3,4,6)
-     te(5,2,2) = te(1,2,6)
-     te(5,4,4) = te(5,2,2)
-     te(5,6,6) = te(1,2,6) * three / (beta * gamma) ** 2
   endif
-  
-  !---- Track orbit.
-  if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
 
 end SUBROUTINE tmdrf
 
diff --git a/src/util.f90 b/src/util.f90
index 7b066e0a8..b26a2f823 100644
--- a/src/util.f90
+++ b/src/util.f90
@@ -352,7 +352,7 @@ module twisslfi
   public
   logical :: centre=.false., first
   logical :: rmatrix=.false., sectormap=.false., ripken=.false.
-  logical :: mode_flip=.false.
+  logical :: mode_flip=.false., exact_expansion=.false.
   logical :: ele_body=.false.
 
 
diff --git a/tests/test-twiss-exact/madx_exact.tfs.cfg b/tests/test-twiss-exact/madx_exact.tfs.cfg
new file mode 100644
index 000000000..787d44248
--- /dev/null
+++ b/tests/test-twiss-exact/madx_exact.tfs.cfg
@@ -0,0 +1,3 @@
+44-48   *     skip         # origin, date and time
+*       *     abs=1e-298
+
diff --git a/tests/test-twiss-exact/madx_exact.tfs.ref b/tests/test-twiss-exact/madx_exact.tfs.ref
new file mode 100644
index 000000000..312b031fe
--- /dev/null
+++ b/tests/test-twiss-exact/madx_exact.tfs.ref
@@ -0,0 +1,61 @@
+@ NAME             %05s "TWISS"
+@ TYPE             %05s "TWISS"
+@ SEQUENCE         %03s "SEQ"
+@ PARTICLE         %06s "PROTON"
+@ MASS             %le  0.93827208816
+@ CHARGE           %le  1.00000000000
+@ ENERGY           %le  2.00000000000
+@ PC               %le  1.76625181913
+@ GAMMA            %le  2.13157784958
+@ KBUNCH           %le  1.00000000000
+@ BCURRENT         %le  0.00000000000
+@ SIGE             %le  0.00100000000
+@ SIGT             %le  1.00000000000
+@ NPART            %le  1.00000000000
+@ EX               %le  1.00000000000
+@ EY               %le  1.00000000000
+@ ET               %le  0.00100000000
+@ BV_FLAG          %le  1.00000000000
+@ LENGTH           %le  10.00000000000
+@ ALFA             %le  0.00000000000
+@ ORBIT5           %le  -0.00000000000
+@ GAMMATR          %le  0.00000000000
+@ Q1               %le  0.23414041683
+@ Q2               %le  0.21860189884
+@ DQ1              %le  0.00000000000
+@ DQ2              %le  0.00000000000
+@ DXMAX            %le  0.00883657376
+@ DYMAX            %le  0.17673147524
+@ XCOMAX           %le  0.01000200560
+@ YCOMAX           %le  0.20004011206
+@ BETXMAX          %le  101.04031624715
+@ BETYMAX          %le  52.06009815221
+@ XCORMS           %le  0.00638635935
+@ YCORMS           %le  0.12772718707
+@ DXRMS            %le  0.00564222194
+@ DYRMS            %le  0.11284443887
+@ DELTAP           %le  0.00000000000
+@ SYNCH_1          %le  0.00000000000
+@ SYNCH_2          %le  0.00000000000
+@ SYNCH_3          %le  0.00000000000
+@ SYNCH_4          %le  0.00000000000
+@ SYNCH_5          %le  0.00000000000
+@ TITLE            %08s "no-title"
+@ ORIGIN           %16s "5.07.00 Linux 64"
+@ DATE             %08s "25/08/21"
+@ TIME             %08s "00.30.01"
+* NAME                           S             X             Y            PT          BETX          BETY 
+$ %s                           %le           %le           %le           %le           %le           %le 
+ "SEQ$START"         0.00000000000 0.00000000000 0.00000000000 0.00000000000 1.00000000000 2.00000000000
+ "DRIFT_0"           1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
+ "M1"                1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
+ "DRIFT_1"           3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
+ "M2"                3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
+ "DRIFT_2"           5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
+ "M3"                5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
+ "DRIFT_3"           7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
+ "M4"                7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
+ "DRIFT_4"           9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
+ "M5"                9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
+ "DRIFT_5"           10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
+ "SEQ$END"           10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
diff --git a/tests/test-twiss-exact/ptc.tfs.cfg b/tests/test-twiss-exact/ptc.tfs.cfg
new file mode 100644
index 000000000..f3647ba00
--- /dev/null
+++ b/tests/test-twiss-exact/ptc.tfs.cfg
@@ -0,0 +1,2 @@
+85-88   *   skip        # date, version
+* * abs=1e-12
diff --git a/tests/test-twiss-exact/ptc.tfs.ref b/tests/test-twiss-exact/ptc.tfs.ref
new file mode 100644
index 000000000..ea4f9c134
--- /dev/null
+++ b/tests/test-twiss-exact/ptc.tfs.ref
@@ -0,0 +1,103 @@
+@ NAME             %05s "TWISS"
+@ TYPE             %05s "TWISS"
+@ SEQUENCE         %03s "SEQ"
+@ PARTICLE         %06s "PROTON"
+@ MASS             %le  0.93827208816
+@ CHARGE           %le  1.00000000000
+@ ENERGY           %le  2.00000000000
+@ PC               %le  1.76625181913
+@ GAMMA            %le  2.13157784958
+@ KBUNCH           %le  1.00000000000
+@ BCURRENT         %le  0.00000000000
+@ SIGE             %le  0.00100000000
+@ SIGT             %le  1.00000000000
+@ NPART            %le  1.00000000000
+@ EX               %le  1.00000000000
+@ EY               %le  1.00000000000
+@ ET               %le  0.00100000000
+@ DELTAP           %le  0.00000000000
+@ LENGTH           %le  10.00000000000
+@ ALPHA_C          %le  -1000000.00000000000
+@ ALPHA_C_P        %le  -1000000.00000000000
+@ ALPHA_C_P2       %le  -1000000.00000000000
+@ ALPHA_C_P3       %le  -1000000.00000000000
+@ ETA_C            %le  -1000000.00000000000
+@ GAMMA_TR         %le  -1000000.00000000000
+@ Q1               %le  0.23414041683
+@ Q2               %le  0.21860189884
+@ DQ1              %le  -1000000.00000000000
+@ DQ2              %le  -1000000.00000000000
+@ QS               %le  -1000000.00000000000
+@ BETA_X_MIN       %le  1.00000000000
+@ BETA_X_MAX       %le  101.04031624715
+@ BETA_Y_MIN       %le  2.00000000000
+@ BETA_Y_MAX       %le  52.06009815221
+@ BETA11MIN        %le  1.00000000000
+@ BETA11MAX        %le  101.04031624715
+@ BETA12MIN        %le  0.00000000000
+@ BETA12MAX        %le  0.00000002002
+@ BETA13MIN        %le  0.00000000000
+@ BETA13MAX        %le  0.00002567481
+@ BETA21MIN        %le  0.00000000000
+@ BETA21MAX        %le  0.00000004005
+@ BETA22MIN        %le  2.00000000000
+@ BETA22MAX        %le  52.06009815221
+@ BETA23MIN        %le  0.00000000000
+@ BETA23MAX        %le  0.01026992554
+@ BETA31MIN        %le  0.00000000000
+@ BETA31MAX        %le  0.00012837407
+@ BETA32MIN        %le  0.00000000000
+@ BETA32MAX        %le  0.02567481386
+@ BETA33MIN        %le  5.00000000000
+@ BETA33MAX        %le  6.59915529106
+@ DISP1MIN         %le  -0.00506703206
+@ DISP1MAX         %le  0.00000000000
+@ DISP2MIN         %le  0.00000000000
+@ DISP2MAX         %le  0.00000000000
+@ DISP3MIN         %le  -0.10134064113
+@ DISP3MAX         %le  0.00000000000
+@ DISP4MIN         %le  0.00000000000
+@ DISP4MAX         %le  0.00000000000
+@ ORBIT_X          %le  0.00000000000
+@ ORBIT_PX         %le  0.00100000000
+@ ORBIT_Y          %le  0.00000000000
+@ ORBIT_PY         %le  0.02000000000
+@ ORBIT_PT         %le  0.00000000000
+@ ORBIT_T          %le  0.00000000000
+@ XCORMS           %le  0.00638635935
+@ PXCORMS          %le  0.00100000000
+@ YCORMS           %le  0.12772718707
+@ PYCORMS          %le  0.02000000000
+@ TCORMS           %le  0.00145006893
+@ PTCORMS          %le  0.00000000000
+@ XCOMIN           %le  0.00000000000
+@ XCOMAX           %le  0.01000200560
+@ PXCOMIN          %le  0.00100000000
+@ PXCOMAX          %le  0.00100000000
+@ YCOMIN           %le  0.00000000000
+@ YCOMAX           %le  0.20004011206
+@ PYCOMIN          %le  0.02000000000
+@ PYCOMAX          %le  0.02000000000
+@ TCOMIN           %le  -0.00227102748
+@ TCOMAX           %le  -0.00000000000
+@ PTCOMIN          %le  0.00000000000
+@ PTCOMAX          %le  0.00000000000
+@ TITLE            %08s "no-title"
+@ ORIGIN           %16s "5.07.00 Linux 64"
+@ DATE             %08s "25/08/21"
+@ TIME             %08s "00.30.01"
+* NAME                           S             X             Y            PT          BETX          BETY 
+$ %s                           %le           %le           %le           %le           %le           %le 
+ "SEQ$START"         0.00000000000 0.00000000000 0.00000000000 0.00000000000 1.00000000000 2.00000000000
+ "DRIFT_0"           1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
+ "M1"                1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
+ "DRIFT_1"           3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
+ "M2"                3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
+ "DRIFT_2"           5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
+ "M3"                5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
+ "DRIFT_3"           7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
+ "M4"                7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
+ "DRIFT_4"           9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
+ "M5"                9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
+ "DRIFT_5"           10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
+ "SEQ$END"           10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
diff --git a/tests/test-twiss-exact/test-twiss-exact.cfg b/tests/test-twiss-exact/test-twiss-exact.cfg
new file mode 100644
index 000000000..f455f4806
--- /dev/null
+++ b/tests/test-twiss-exact/test-twiss-exact.cfg
@@ -0,0 +1,4 @@
+1-8     *     skip         # head
+213-215 *     skip         # tail
+*       *     abs=1e-298
+
diff --git a/tests/test-twiss-exact/test-twiss-exact.madx b/tests/test-twiss-exact/test-twiss-exact.madx
new file mode 100644
index 000000000..cc12db1e8
--- /dev/null
+++ b/tests/test-twiss-exact/test-twiss-exact.madx
@@ -0,0 +1,30 @@
+
+seq: sequence, l=10;
+m1:marker, at=1;
+m2:marker, at=3;
+m3:marker, at=5;
+m4:marker, at=7;
+m5:marker, at=9;
+endsequence;
+set,  format="13.11f";
+beam,particle=proton,energy=2;
+use, sequence=seq;
+option, sympl=false;
+px0 = 0.001;
+py0 = 0.02;
+betx0 = 1;
+bety0 = 2;
+select, flag = twiss, column=name,s,x,y,pt,betx,bety;
+
+twiss, file="madx.tfs", betx=betx0, bety=bety0, px=px0, py=py0, pt=pt0;
+twiss, file="madx_exact.tfs", betx=betx0, bety=bety0, px=px0, py=py0, exact;
+
+
+ptc_create_universe;
+ptc_create_layout,model=3, method=2, nst=1, exact=true;
+
+ptc_setswitch, debuglevel=0, fringe=false, totalpath=false	;
+PTC_TWISS, table=twiss, icase=6, no=3, betx=betx0, bety=bety0, betz=5, px=px0,py=py0;
+write, table=twiss, file="ptc.tfs";
+PTC_END;
+stop;
diff --git a/tests/test-twiss-exact/test-twiss-exact.ref b/tests/test-twiss-exact/test-twiss-exact.ref
new file mode 100644
index 000000000..15ab76b9b
--- /dev/null
+++ b/tests/test-twiss-exact/test-twiss-exact.ref
@@ -0,0 +1,205 @@
+
+  ++++++++++++++++++++++++++++++++++++++++++++
+  +     MAD-X 5.07.00  (64 bit, Linux)       +
+  + Support: mad@cern.ch, http://cern.ch/mad +
+  + Release   date: 2021.05.03               +
+  + Execution date: 2021.08.25 00:34:10      +
+  ++++++++++++++++++++++++++++++++++++++++++++
+
+
+seq: sequence, l=10;
+
+m1:marker, at=1;
+
+m2:marker, at=3;
+
+m3:marker, at=5;
+
+m4:marker, at=7;
+
+m5:marker, at=9;
+
+endsequence;
+
+set,  format="13.11f";
+
+beam,particle=proton,energy=2;
+
+use, sequence=seq;
+
+option, sympl=false;
+
+px0 = 0.001;
+
+py0 = 0.02;
+
+betx0 = 1;
+
+bety0 = 2;
+
+select, flag = twiss, column=name,s,x,y,pt,betx,bety;
+
+
+
+twiss, file="madx.tfs", betx=betx0, bety=bety0, px=px0, py=py0, pt=pt0;
+
+enter Twiss module
+ Found 3 initial orbit vector values from twiss command. 
+  
+open line - error with deltap:   0.000000E+00
+initial orbit vector:   0.000000E+00  1.000000E-03  0.000000E+00  2.000000E-02  0.000000E+00  0.000000E+00
+final orbit vector:     1.000000E-02  1.000000E-03  2.000000E-01  2.000000E-02 -2.270344E-03  0.000000E+00
+
+++++++ table: summ
+
+       length        orbit5          alfa       gammatr 
+10.00000000000 -0.00000000000 0.00000000000 0.00000000000 
+
+           q1           dq1       betxmax         dxmax 
+0.23413724128 0.00000000000 101.00000000000 0.01132341367 
+
+        dxrms        xcomax        xcorms            q2 
+0.00723008881 0.01000000000 0.00638507876 0.21858352091 
+
+          dq2       betymax         dymax         dyrms 
+0.00000000000 52.00000000000 0.22646827347 0.14460177627 
+
+       ycomax        ycorms        deltap       synch_1 
+0.20000000000 0.12770157520 0.00000000000 0.00000000000 
+
+      synch_2       synch_3       synch_4       synch_5 
+0.00000000000 0.00000000000 0.00000000000 0.00000000000 
+
+      synch_6       synch_8        nflips 
+0.00000000000 0.00000000000 0.00000000000 
+twiss, file="madx_exact.tfs", betx=betx0, bety=bety0, px=px0, py=py0, exact;
+
+enter Twiss module
+ Found 2 initial orbit vector values from twiss command. 
+  
+open line - error with deltap:   0.000000E+00
+initial orbit vector:   0.000000E+00  1.000000E-03  0.000000E+00  2.000000E-02  0.000000E+00  0.000000E+00
+final orbit vector:     1.000201E-02  1.000000E-03  2.000401E-01  2.000000E-02 -1.771200E-03  0.000000E+00
+
+++++++ table: summ
+
+       length        orbit5          alfa       gammatr 
+10.00000000000 -0.00000000000 0.00000000000 0.00000000000 
+
+           q1           dq1       betxmax         dxmax 
+0.23414041683 0.00000000000 101.04031624715 0.00883657376 
+
+        dxrms        xcomax        xcorms            q2 
+0.00564222194 0.01000200560 0.00638635935 0.21860189884 
+
+          dq2       betymax         dymax         dyrms 
+0.00000000000 52.06009815221 0.17673147524 0.11284443887 
+
+       ycomax        ycorms        deltap       synch_1 
+0.20004011206 0.12772718707 0.00000000000 0.00000000000 
+
+      synch_2       synch_3       synch_4       synch_5 
+0.00000000000 0.00000000000 0.00000000000 0.00000000000 
+
+      synch_6       synch_8        nflips 
+0.00000000000 0.00000000000 0.00000000000 
+
+
+
+
+ptc_create_universe;
+
+Determined SECTOR NMUL MAX :  0
+ptc_create_layout,model=3, method=2, nst=1, exact=true;
+
+ MAD-X Beam Parameters
+                  Energy : 0.200000E+01
+          Kinetic Energy : 0.106173E+01
+      Particle Rest Mass : 0.938272E+00
+                Momentum : 0.176625E+01
+ Setting MADx with 
+     energy           2.0000000000000000     
+     method                   2
+     Num. of steps            1
+     charge           1.0000000000000000     
+  Length of machine:    10.000000000000000     
+ The machine is a RING
+ ------------------------------------ PTC Survey ------------------------------------
+ Before start:    0.0000000000000000        0.0000000000000000        0.0000000000000000     
+ Before   end:    0.0000000000000000        0.0000000000000000        0.0000000000000000     
+ After  start:    0.0000000000000000        0.0000000000000000        0.0000000000000000     
+ After    end:    0.0000000000000000        0.0000000000000000        10.000000000000000     
+
+
+ptc_setswitch, debuglevel=0, fringe=false, totalpath=false	;
+
+PTC_TWISS, table=twiss, icase=6, no=3, betx=betx0, bety=bety0, betz=5, px=px0,py=py0;
+
+++++++ warning: my_state:   no cavity - dimensionality reduced 6 -> 5 and 1/2
+ Reduced SUMM Table (Inital parameters specified)
+
+++++++ table: ptc_twiss_summary
+
+       length       alpha_c     alpha_c_p    alpha_c_p2 
+10.00000000000 -1000000.00000000000 -1000000.00000000000 -1000000.00000000000 
+
+   alpha_c_p3         eta_c      gamma_tr            q1 
+-1000000.00000000000 -1000000.00000000000 -1000000.00000000000 0.23414041683 
+
+           q2           dq1           dq2            qs 
+0.21860189884 -1000000.00000000000 -1000000.00000000000 -1000000.00000000000 
+
+   beta_x_min    beta_x_max    beta_y_min    beta_y_max 
+1.00000000000 101.04031624715 2.00000000000 52.06009815221 
+
+    beta11min     beta11max     beta12min     beta12max 
+1.00000000000 101.04031624715 0.00000000000 0.00000002002 
+
+    beta13min     beta13max     beta21min     beta21max 
+0.00000000000 0.00002567481 0.00000000000 0.00000004005 
+
+    beta22min     beta22max     beta23min     beta23max 
+2.00000000000 52.06009815221 0.00000000000 0.01026992554 
+
+    beta31min     beta31max     beta32min     beta32max 
+0.00000000000 0.00012837407 0.00000000000 0.02567481386 
+
+    beta33min     beta33max      disp1min      disp1max 
+5.00000000000 6.59915529106 -0.00506703206 0.00000000000 
+
+     disp2min      disp2max      disp3min      disp3max 
+0.00000000000 0.00000000000 -0.10134064113 0.00000000000 
+
+     disp4min      disp4max        deltap       orbit_x 
+0.00000000000 0.00000000000 0.00000000000 0.00000000000 
+
+     orbit_px       orbit_y      orbit_py      orbit_pt 
+0.00100000000 0.00000000000 0.02000000000 0.00000000000 
+
+      orbit_t        xcorms        ycorms       pxcorms 
+0.00000000000 0.00638635935 0.12772718707 0.00100000000 
+
+      pycorms        tcorms       ptcorms        xcomax 
+0.02000000000 0.00145006893 0.00000000000 0.01000200560 
+
+       ycomax       pxcomax       pycomax        tcomax 
+0.20004011206 0.00100000000 0.02000000000 -0.00000000000 
+
+      ptcomax        xcomin        ycomin       pxcomin 
+0.00000000000 0.00000000000 0.00000000000 0.00100000000 
+
+      pycomin        tcomin       ptcomin 
+0.02000000000 -0.00227102748 0.00000000000 
+write, table=twiss, file="ptc.tfs";
+
+PTC_END;
+
+stop;
+
+
+  Number of warnings: 1
+0 in C and 1 in Fortran
+
+  ++++++++++++++++++++++++++++++++++++++++++++
+  +          MAD-X finished normally         +
+  ++++++++++++++++++++++++++++++++++++++++++++