From 7d1943022432caeb61e0a63c13dea3f708587eb7 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Thu, 29 Jun 2023 15:26:13 +0200
Subject: [PATCH 01/14] test of the variables changes for the sextupole for the
 tapering

---
 src/twiss.f90 | 88 +++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 86 insertions(+), 2 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 1748f195d..ada23117c 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -6592,6 +6592,7 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
   integer, external :: el_par_vector, node_fd_errors
   double precision, external :: node_value, get_value
   double precision :: bet0, bet_sqr, f_damp_t
+  double precision :: newbet0, deltaplusone
 
   integer, external :: get_option
   character(len=name_len) el_name
@@ -6646,8 +6647,91 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
      orbit(6) = orbit(6) * (one - rfac) - rfac / bet0;
   endif
 
-  call sxbody(fsec,ftrk,tilt,sk2,orbit,dl,ek,re,te)
-  if (fcentre) return
+  if (exact_expansion) then
+     
+     pt = orbit(6)
+     deltaplusone = sqrt(pt**2+2*pt/bet0+1)
+     newbet0 = (deltaplusone)/ (1/bet0+pt)
+     orbit(2) = orbit(2)/deltaplusone
+     orbit(4) = orbit(4)/deltaplusone
+     orbit(5) = orbit(5)*(bet0/newbet0)
+     sk2 = sk2/deltaplusone
+     orbit(6) = 0
+
+     call sxbody(fsec,ftrk,tilt,sk2,orbit,dl,ek,re,te)
+
+     orbit(2) = orbit(2)*deltaplusone
+     orbit(4) = orbit(4)*deltaplusone
+     orbit(5) = orbit(5)*(newbet0/bet0)
+     sk2 = sk2*deltaplusone
+     orbit(6) = pt
+
+     !---- First order terms
+     
+     re(1,2) = re(1,2)*deltaplusone
+     re(3,4) = re(3,4)*deltaplusone
+     re(5,6) = re(5,6)*(bet0/newbet0)
+
+     ek(2) = ek(2)/deltaplusone
+     ek(4) = ek(4)/deltaplusone
+     ek(5) = ek(5)*(bet0/newbet0)
+
+     !---- Second order terms
+
+     te(1,1,2) = te(1,1,2)*deltaplusone
+     te(1,2,2) = te(1,2,2)*deltaplusone**2
+     te(1,3,4) = te(1,3,4)*deltaplusone
+     te(1,4,4) = te(1,4,4)*deltaplusone**2
+
+     te(2,2,2) = te(2,2,2)*deltaplusone
+     te(2,3,3) = te(2,3,3)/deltaplusone
+     te(2,4,4) = te(2,4,4)*deltaplusone
+
+     te(3,1,4) = te(3,1,4)*deltaplusone
+     te(3,2,3) = te(3,2,3)*deltaplusone
+     te(3,2,4) = te(3,2,4)*deltaplusone**2
+
+     te(4,1,3) = te(4,1,3)/deltaplusone
+     te(4,2,4) = te(4,2,4)*deltaplusone
+
+     if (fcentre) return
+
+     !---- First order terms
+
+     re(1,2) = re(1,2)/deltaplusone
+     re(3,4) = re(3,4)/deltaplusone
+     re(5,6) = re(5,6)/(bet0/newbet0)
+
+     ek(2) = ek(2)*deltaplusone
+     ek(4) = ek(4)*deltaplusone
+     ek(5) = ek(5)/(bet0/newbet0)
+
+     !---- Second order terms
+
+     te(1,1,2) = te(1,1,2)/deltaplusone
+     te(1,2,2) = te(1,2,2)/deltaplusone**2
+     te(1,3,4) = te(1,3,4)/deltaplusone
+     te(1,4,4) = te(1,4,4)/deltaplusone**2
+
+     te(2,2,2) = te(2,2,2)/deltaplusone
+     te(2,3,3) = te(2,3,3)*deltaplusone
+     te(2,4,4) = te(2,4,4)/deltaplusone
+
+     te(3,1,4) = te(3,1,4)/deltaplusone
+     te(3,2,3) = te(3,2,3)/deltaplusone
+     te(3,2,4) = te(3,2,4)/deltaplusone**2
+
+     te(4,1,3) = te(4,1,3)*deltaplusone
+     te(4,2,4) = te(4,2,4)/deltaplusone
+
+  else
+
+     call sxbody(fsec,ftrk,tilt,sk2,orbit,dl,ek,re,te)
+     if (fcentre) return
+
+  endif
+
+
 
   !---- Half radiation effects at exit.
   if (ftrk) then

From 1d9fd5a1bd610f82f1cb2121cce9dccafa7b1e99 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Thu, 29 Jun 2023 18:00:15 +0200
Subject: [PATCH 02/14] test variables changes for sextupoles (corrected)

---
 src/twiss.f90 | 42 +++++++-----------------------------------
 1 file changed, 7 insertions(+), 35 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index ada23117c..54de7b44e 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -6651,10 +6651,10 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
      
      pt = orbit(6)
      deltaplusone = sqrt(pt**2+2*pt/bet0+1)
-     newbet0 = (deltaplusone)/ (1/bet0+pt)
+     !newbet0 = (deltaplusone)/ (1/bet0+pt)
      orbit(2) = orbit(2)/deltaplusone
      orbit(4) = orbit(4)/deltaplusone
-     orbit(5) = orbit(5)*(bet0/newbet0)
+     orbit(5) = orbit(5)!*(bet0/newbet0)
      sk2 = sk2/deltaplusone
      orbit(6) = 0
 
@@ -6662,49 +6662,19 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
 
      orbit(2) = orbit(2)*deltaplusone
      orbit(4) = orbit(4)*deltaplusone
-     orbit(5) = orbit(5)*(newbet0/bet0)
+     orbit(5) = orbit(5)!*(newbet0/bet0)
      sk2 = sk2*deltaplusone
      orbit(6) = pt
 
-     !---- First order terms
-     
-     re(1,2) = re(1,2)*deltaplusone
-     re(3,4) = re(3,4)*deltaplusone
-     re(5,6) = re(5,6)*(bet0/newbet0)
-
-     ek(2) = ek(2)/deltaplusone
-     ek(4) = ek(4)/deltaplusone
-     ek(5) = ek(5)*(bet0/newbet0)
-
-     !---- Second order terms
-
-     te(1,1,2) = te(1,1,2)*deltaplusone
-     te(1,2,2) = te(1,2,2)*deltaplusone**2
-     te(1,3,4) = te(1,3,4)*deltaplusone
-     te(1,4,4) = te(1,4,4)*deltaplusone**2
-
-     te(2,2,2) = te(2,2,2)*deltaplusone
-     te(2,3,3) = te(2,3,3)/deltaplusone
-     te(2,4,4) = te(2,4,4)*deltaplusone
-
-     te(3,1,4) = te(3,1,4)*deltaplusone
-     te(3,2,3) = te(3,2,3)*deltaplusone
-     te(3,2,4) = te(3,2,4)*deltaplusone**2
-
-     te(4,1,3) = te(4,1,3)/deltaplusone
-     te(4,2,4) = te(4,2,4)*deltaplusone
-
-     if (fcentre) return
-
      !---- First order terms
 
      re(1,2) = re(1,2)/deltaplusone
      re(3,4) = re(3,4)/deltaplusone
-     re(5,6) = re(5,6)/(bet0/newbet0)
+     re(5,6) = re(5,6)!/(bet0/newbet0)
 
      ek(2) = ek(2)*deltaplusone
      ek(4) = ek(4)*deltaplusone
-     ek(5) = ek(5)/(bet0/newbet0)
+     ek(5) = ek(5)!/(bet0/newbet0)
 
      !---- Second order terms
 
@@ -6724,6 +6694,8 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
      te(4,1,3) = te(4,1,3)*deltaplusone
      te(4,2,4) = te(4,2,4)/deltaplusone
 
+     if (fcentre) return
+
   else
 
      call sxbody(fsec,ftrk,tilt,sk2,orbit,dl,ek,re,te)

From a2d8a5278be2b6659d83df04097b2b3d64fc53ea Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Thu, 27 Jul 2023 14:08:59 +0200
Subject: [PATCH 03/14] new method applied for both quadrupoles and sextupoles

---
 src/twiss.f90 | 191 +++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 143 insertions(+), 48 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 54de7b44e..7fd291cc2 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -6112,8 +6112,8 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
   use twtrrfi
   use twisslfi
   use twiss_elpfi
-  use twissbeamfi, only : radiate, deltap, gamma, arad
-  use math_constfi, only : zero, one, two, three
+  use twissbeamfi, only : radiate, deltap, beta, gamma, arad
+  use math_constfi, only : zero, one, two, three, four
   use name_lenfi
   implicit none
   !----------------------------------------------------------------------*
@@ -6142,14 +6142,14 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
 
   logical :: cplxy
   integer :: i, j, n_ferr, elpar_vl
-  double precision :: ct, st, tmp
+  double precision :: ct, st, tmp, biby4
   double precision :: f_errors(0:maxferr)
   double precision :: tilt, sk1, pt, sk1s, bvk, rfac
 
   integer, external :: node_fd_errors
   integer, external :: el_par_vector
   double precision, external :: node_value, get_value
-  double precision :: bet0, bet_sqr, f_damp_t
+  double precision :: bet0, bet_sqr, f_damp_t, oneplusdelta
 
   !double precision :: newbet0, newdeltap, ff
   integer, external :: get_option
@@ -6207,31 +6207,65 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
      orbit(6) = orbit(6) * (one - rfac) - rfac / bet0;
   endif
 
-  ! absorb pt in deltas
-  !pt= orbit(6)
-  !newdeltap=sqrt(pt**2+2*pt/bet0+1)-1
-  !newbet0= (1+newdeltap)/ (1/bet0+pt)
-  !ff= (one + deltap) / ( deltap+newdeltap) ! ratio
-  !orbit(2)=orbit(2)*ff
-  !orbit(4)=orbit(4)*ff
-  !orbit(6)=0
-  !sk1 =sk1*ff
-  ! start absorb pt in deltas
+  if (exact_expansion) then
 
-  call qdbody(fsec,ftrk,tilt,sk1,orbit,dl,ek,re,te)
+     pt = orbit(6)
+     oneplusdelta = sqrt(pt*pt + 2*pt/bet0 + 1)
+     orbit(2) = orbit(2)/oneplusdelta
+     orbit(4) = orbit(4)/oneplusdelta
+     sk1 = sk1/oneplusdelta
 
-  ! restore pt
-  !sk1 =sk1/ff
-  !orbit(2)=orbit(2)/ff
-  !orbit(4)=orbit(4)/ff
-  !orbit(6)=pt
-  ! end restore pt
+     call qdbody(fsec,ftrk,tilt,sk1,orbit,dl,ek,re,te)
 
+     orbit(2) = orbit(2)*oneplusdelta
+     orbit(4) = orbit(4)*oneplusdelta
+     sk1 = sk1*oneplusdelta
+     orbit(6) = pt
 
+     !---- First order terms
 
+     re(1,2) = re(1,2)/oneplusdelta
+     re(2,1) = re(2,1)*oneplusdelta
+     re(3,4) = re(3,4)/oneplusdelta
+     re(4,3) = re(4,3)*oneplusdelta
 
+     !---- Second order terms
 
-  if (fcentre) return
+     if (fsec) then
+
+        biby4 = one / (four * beta)
+
+        te(1,2,6) = te(1,2,6)/oneplusdelta
+        te(1,6,2) = te(1,6,2)/oneplusdelta
+
+        te(2,1,6) = te(2,1,6)*oneplusdelta
+        te(2,6,1) = te(2,6,1)*oneplusdelta
+
+        te(3,4,6) = te(3,4,6)/oneplusdelta
+        te(3,6,4) = te(3,6,4)/oneplusdelta
+
+        te(4,3,6) = te(4,3,6)*oneplusdelta
+        te(4,6,3) = te(4,6,3)*oneplusdelta
+
+        te(5,1,2) = te(5,1,2)/oneplusdelta
+        te(5,2,1) = te(5,2,1)/oneplusdelta
+        te(5,2,2) = te(5,2,2)/oneplusdelta**2
+
+        te(5,3,4) = te(5,3,4)/oneplusdelta
+        te(5,4,3) = te(5,4,3)/oneplusdelta
+        te(5,4,4) = te(5,4,4)/oneplusdelta**2
+
+     endif
+
+     if (fcentre) return
+
+  else
+
+     call qdbody(fsec,ftrk,tilt,sk1,orbit,dl,ek,re,te)
+
+     if (fcentre) return
+
+  endif
 
   !---- Half radiation effect at exit.
   if (radiate .and. ftrk) then
@@ -6257,6 +6291,7 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
 end SUBROUTINE tmquad
 
 SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
+  use twisslfi, only: exact_expansion
   use twissbeamfi, only : beta, gamma, dtbyds
   use math_constfi, only : zero, one, two, four, six, ten3m
   implicit none
@@ -6282,17 +6317,34 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
 
   double precision :: qk, qkl, qkl2
   double precision :: cx, sx, cy, sy, biby4
+  double precision :: x,px,y,py,t,pt,deltaplusone,nk1
+
+
+  !if (exact_expansion) then
+  !   x= orbit(1)
+  !   px= orbit(2)
+  !   y= orbit(3)
+  !   py= orbit(4)
+  !   t= orbit(5)
+  !   pt= orbit(6)
+  !   deltaplusone=sqrt(pt**2+2*pt/beta+1)
+  !   nk1=sk1/deltaplusone
+  !else
+  !   nk1= sk1
+  !endif
+
+  nk1 = sk1
 
   !---- Set up c's and s's.
-  qk = sqrt(abs(sk1))
+  qk = sqrt(abs(nk1))
   qkl = qk * el
   if (abs(qkl) .lt. ten3m) then
-     qkl2 = sk1 * el**2
+     qkl2 = nk1 * el**2
      cx = (one - qkl2 / two)
      sx = (one - qkl2 / six) * el
      cy = (one + qkl2 / two)
      sy = (one + qkl2 / six) * el
-  else if (sk1 .gt. zero) then
+  else if (nk1 .gt. zero) then
      cx = cos(qkl)
      sx = sin(qkl) / qk
      cy = cosh(qkl)
@@ -6304,54 +6356,85 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
      sy = sin(qkl) / qk
   endif
 
+  !if (exact_expansion) then
+  !---- First-order terms.
+  !re(1,1) = cx
+  !re(1,2) = sx/deltaplusone
+  !re(2,1) = - nk1 * sx * deltaplusone
+  !re(2,2) = cx
+  !re(3,3) = cy
+  !re(3,4) = sy /deltaplusone
+  !re(4,3) = + nk1 * sy *deltaplusone
+  !re(4,4) = cy
+  !re(5,6) = el/(beta*gamma)**2
+
+  !ek(5) = el*dtbyds ! to be checked
+
+  !else
+  
   !---- First-order terms.
   re(1,1) = cx
   re(1,2) = sx
-  re(2,1) = - sk1 * sx
+  re(2,1) = - nk1 * sx
   re(2,2) = cx
   re(3,3) = cy
   re(3,4) = sy
-  re(4,3) = + sk1 * sy
+  re(4,3) = + nk1 * sy
   re(4,4) = cy
   re(5,6) = el/(beta*gamma)**2
 
   ek(5) = el*dtbyds
 
+  !endif
+
   !---- Second-order terms.
   if (fsec) then
      biby4 = one / (four * beta)
 
-     te(1,1,6) = + sk1 * el * sx * biby4
+     te(1,1,6) = + nk1 * el * sx * biby4
      te(1,6,1) = te(1,1,6)
      te(2,2,6) = te(1,1,6)
      te(2,6,2) = te(1,1,6)
      te(1,2,6) = - (sx + el*cx) * biby4
      te(1,6,2) = te(1,2,6)
-     te(2,1,6) = - sk1 * (sx - el*cx) * biby4
+     te(2,1,6) = - nk1 * (sx - el*cx) * biby4
      te(2,6,1) = te(2,1,6)
 
-     te(3,3,6) = - sk1 * el * sy * biby4
+     te(3,3,6) = - nk1 * el * sy * biby4
      te(3,6,3) = te(3,3,6)
      te(4,4,6) = te(3,3,6)
      te(4,6,4) = te(3,3,6)
      te(3,4,6) = - (sy + el*cy) * biby4
      te(3,6,4) = te(3,4,6)
-     te(4,3,6) = + sk1 * (sy - el*cy) * biby4
+     te(4,3,6) = + nk1 * (sy - el*cy) * biby4
      te(4,6,3) = te(4,3,6)
 
-     te(5,1,1) = - sk1 * (el - sx*cx) * biby4
-     te(5,1,2) = + sk1 * sx**2 * biby4
+     te(5,1,1) = - nk1 * (el - sx*cx) * biby4
+     te(5,1,2) = + nk1 * sx**2 * biby4
      te(5,2,1) = te(5,1,2)
      te(5,2,2) = - (el + sx*cx) * biby4
-     te(5,3,3) = + sk1 * (el - sy*cy) * biby4
-     te(5,3,4) = - sk1 * sy**2 * biby4
+     te(5,3,3) = + nk1 * (el - sy*cy) * biby4
+     te(5,3,4) = - nk1 * sy**2 * biby4
      te(5,4,3) = te(5,3,4)
      te(5,4,4) = - (el + sy*cy) * biby4
      te(5,6,6) = (- six * re(5,6)) * biby4
   endif
 
   !---- Track orbit.
+  !if (exact_expansion) then
+
+  !   orbit(1)=cx*x + sx*px/deltaplusone
+  !   orbit(2)=-nk1 * sx * x*deltaplusone + cx*px
+  !   orbit(3)=cy*y + sy*py/deltaplusone
+  !   orbit(4)=nk1 * sy * y*deltaplusone + cy*py
+  !   orbit(5)=el/(beta*gamma)**2*pt
+  !   re(5,1)=re(5,1) + te(5,1,1)*x + te(5,1,2)*px + te(5,1,3)*y + te(5,1,4)*py + te(5,1,5)*t + te(5,1,6)*pt
+  !else
+  !   if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
+  !endif
+  
   if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
+  
   !---- Apply tilt.
   if (tilt .ne. zero) call tmtilt(fsec,tilt,ek,re,te)
 
@@ -6587,7 +6670,7 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
   integer :: i, j, n_ferr, elpar_vl
   double precision :: ct, st, tmp
   double precision :: f_errors(0:maxferr)
-  double precision :: tilt, sk2, pt, sk2s, bvk, rfac
+  double precision :: tilt, sk2, pt, sk2s, bvk, rfac, skl
 
   integer, external :: el_par_vector, node_fd_errors
   double precision, external :: node_value, get_value
@@ -6678,21 +6761,33 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
 
      !---- Second order terms
 
-     te(1,1,2) = te(1,1,2)/deltaplusone
-     te(1,2,2) = te(1,2,2)/deltaplusone**2
-     te(1,3,4) = te(1,3,4)/deltaplusone
-     te(1,4,4) = te(1,4,4)/deltaplusone**2
+     if (fsec) then 
+        skl = sk2 * el
+        if (skl .ne. zero) then
+
+           te(1,1,2) = te(1,1,2)/deltaplusone
+           te(1,2,2) = te(1,2,2)/deltaplusone**2
+           te(1,3,4) = te(1,3,4)/deltaplusone
+           te(1,4,4) = te(1,4,4)/deltaplusone**2
 
-     te(2,2,2) = te(2,2,2)/deltaplusone
-     te(2,3,3) = te(2,3,3)*deltaplusone
-     te(2,4,4) = te(2,4,4)/deltaplusone
+           te(2,2,2) = te(2,2,2)/deltaplusone
+           te(2,3,3) = te(2,3,3)*deltaplusone
+           te(2,4,4) = te(2,4,4)/deltaplusone
 
-     te(3,1,4) = te(3,1,4)/deltaplusone
-     te(3,2,3) = te(3,2,3)/deltaplusone
-     te(3,2,4) = te(3,2,4)/deltaplusone**2
+           te(3,1,4) = te(3,1,4)/deltaplusone
+           te(3,2,3) = te(3,2,3)/deltaplusone
+           te(3,2,4) = te(3,2,4)/deltaplusone**2
 
-     te(4,1,3) = te(4,1,3)*deltaplusone
-     te(4,2,4) = te(4,2,4)/deltaplusone
+           te(4,1,3) = te(4,1,3)*deltaplusone
+           te(4,2,4) = te(4,2,4)/deltaplusone
+           
+        endif
+
+        te(1,2,6) = te(1,2,6)/deltaplusone
+        te(3,4,6) = te(3,4,6)/deltaplusone
+        te(5,2,2) = te(5,2,2)/deltaplusone**2
+        te(5,4,4) = te(5,4,4)/deltaplusone**2
+     endif
 
      if (fcentre) return
 

From cad3ec08c4336b9af037aa28f5c14911926df32a Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Wed, 9 Aug 2023 15:44:37 +0200
Subject: [PATCH 04/14] first try bending magnets

---
 src/twiss.f90 | 195 ++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 181 insertions(+), 14 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 7fd291cc2..2f97221f9 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -3814,7 +3814,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
   !     re(6,6)   (double)  transfer matrix.                             *
   !     te(6,6,6) (double)  second-order terms.                          *
   !----------------------------------------------------------------------*
-  logical :: ftrk, fmap, fcentre
+  logical :: fsec, ftrk, fmap, fcentre
   double precision :: orbit(6), ek(6), re(6,6), te(6,6,6), el, dl
 
   logical :: cplxy
@@ -3824,13 +3824,13 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
   double precision :: f_errors(0:maxferr)
   double precision :: rw(6,6), tw(6,6,6), ek0(6)
   double precision :: x, y
-  double precision :: an, sk0, sk1, sk2, sks, tilt, e1, e2, h, h1, h2, hgap, fint, fintx, rhoinv, blen, bvk, ktap, angle
+  double precision :: an, sk0, sk1, sk2, sks, tilt, e1, e2, h, h1, h2, hgap, fint, fintx, rhoinv, blen, bvk, ktap, angle, sig
   double precision :: dh, corr, ct, st, hx, hy, rfac, pt, h_k
 
   integer, external :: el_par_vector, node_fd_errors
   double precision, external :: node_value, get_value
   character(len=name_len) :: name
-  double precision :: bet0, bet_sqr, f_damp_t
+  double precision :: bet0, bet_sqr, f_damp_t, oneplusdelta
 
   integer, external :: get_option
   
@@ -3927,21 +3927,187 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
   
   !---- Body of the dipole.
   !---- Get map for body section
-  call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
+
+!  if (exact_expansion) then
+     
+!     pt = orbit(6)
+!     oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+!     orbit(2) = orbit(2)/oneplusdelta
+!     orbit(4) = orbit(4)/oneplusdelta
+!     sk1 = sk1/oneplusdelta
+!     sk2 = sk2/oneplusdelta
+!     orbit(6) = 0
+
+!     call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
+
+!     orbit(2) = orbit(2)*oneplusdelta
+!     orbit(4) = orbit(4)*oneplusdelta
+!     sk1 = sk1*oneplusdelta
+!     sk2 = sk2*oneplusdelta
+!     orbit(6) = pt
+
+     !----First order terms
+
+!     re(1,2) = re(1,2)/oneplusdelta
+!     re(2,1) = re(2,1)*oneplusdelta
+!     re(2,6) = re(2,6)/oneplusdelta
+!     re(3,4) = re(3,4)/oneplusdelta
+!     re(5,2) = re(5,2)*oneplusdelta
+
+!     ek(2) = ek(2)*oneplusdelta
+
+     !----Second order terms
+
+!     if (fsec) then
+
+!        te(1,2,2) = te(1,2,2)/oneplusdelta**2
+!        te(1,2,6) = te(1,2,6)/oneplusdelta**2
+!        te(2,1,1) = te(2,1,1)*oneplusdelta
+!        te(2,2,2) = te(2,2,2)/oneplusdelta
+!        te(2,1,6) = te(2,1,6)*oneplusdelta
+!        te(2,6,6) = te(2,6,6)*oneplusdelta
+!        te(5,1,2) = te(5,1,2)/oneplusdelta
+!        te(5,2,2) = te(5,2,2)/oneplusdelta**2
+!        te(5,2,6) = te(5,2,6)/oneplusdelta
+
+        !----Mixed terms
+
+!        te(1,3,4) = te(1,3,4)/oneplusdelta
+!        te(1,4,4) = te(1,4,4)/oneplusdelta**2
+!        te(2,3,3) = te(2,3,3)*oneplusdelta
+!        te(2,4,4) = te(2,4,4)/oneplusdelta
+!        te(3,1,4) = te(3,1,4)/oneplusdelta
+!        te(3,2,3) = te(3,2,3)/oneplusdelta
+!        te(3,2,4) = te(3,2,4)/oneplusdelta**2
+ !       te(3,4,6) = te(3,4,6)/oneplusdelta
+!        te(4,1,3) = te(4,1,3)*oneplusdelta
+!        te(4,2,4) = te(4,2,4)/oneplusdelta
+!        te(4,3,6) = te(4,3,6)*oneplusdelta
+!        te(5,3,4) = te(5,3,4)/oneplusdelta
+!        te(5,4,4) = te(5,4,4)/oneplusdelta**2
+
+!     endif  
+
+!  else
+
+     call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
+
+!  endif
   
   !---- Get map for entrance fringe field and concatenate
-  if (.not.kill_ent_fringe) then
-     corr = (h_k + h_k) * hgap * fint
-     call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
-     call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
-  endif
+!  if (exact_expansion) then
+     
+!     if(.not.kill_ent_fringe) then
+!        corr = (h_k + h_k) * hgap * fint
+!        pt = orbit(6)
+!        oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+!        orbit(2) = orbit(2)/oneplusdelta
+!        orbit(4) = orbit(4)/oneplusdelta
+!        sk1 = sk1/oneplusdelta
+!        orbit(6) = 0
+        
+!        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+
+!        orbit(2) = orbit(2)*oneplusdelta
+!        orbit(4) = orbit(4)*oneplusdelta
+!        sk1 = sk1*oneplusdelta
+!        orbit(6) = pt
+
+        !----First order terms
+
+!        re(2,1) = re(2,1)*oneplusdelta
+!        re(4,3) = re(4,3)*oneplusdelta
+
+        !----Second order terms
+
+!        if (fsec) then
+
+!           te(2,1,1) = te(2,1,1)*oneplusdelta
+!           te(2,3,3) = te(2,3,3)*oneplusdelta
+
+!           if (sig .gt. zero) then
+
+!              te(2,3,3) = te(2,3,3)*oneplusdelta
+
+!           else
+
+!              te(2,1,1) = te(2,1,1)*oneplusdelta
+!              te(4,1,3) = te(4,1,3)*oneplusdelta
+
+!           endif
+!        endif
+
+!        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+
+!     endif
+
+!  else
+       
+     if (.not.kill_ent_fringe) then
+        corr = (h_k + h_k) * hgap * fint
+        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+     endif
+
+!  endif
   
   !---- Get map for exit fringe fields and concatenate
-  if (.not.kill_exi_fringe) then
-     if (fintx .lt. 0) fintx = fint
-     corr = (h_k + h_k) * hgap * fintx
-     call tmfrng(.true.,h_k,sk1,e2,h2,-one,corr,rw,tw)
-     call tmcat1(.true.,ek0,rw,tw,ek,re,te,ek,re,te)
+
+  if (exact_expansion) then
+     
+     if(.not.kill_exi_fringe) then
+        if (fintx .lt. 0) fintx = fint
+        corr = (h_k + h_k) * hgap * fint
+        pt = orbit(6)
+        oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+        orbit(2) = orbit(2)/oneplusdelta
+        orbit(4) = orbit(4)/oneplusdelta
+        sk1 = sk1/oneplusdelta
+        
+        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+
+        orbit(2) = orbit(2)*oneplusdelta
+        orbit(4) = orbit(4)*oneplusdelta
+        sk1 = sk1*oneplusdelta
+        orbit(6) = pt
+
+        !----First order terms
+
+        re(2,1) = re(2,1)*oneplusdelta
+        re(4,3) = re(4,3)*oneplusdelta
+
+        !----Second order terms
+
+        if (fsec) then
+
+           te(2,1,1) = te(2,1,1)*oneplusdelta
+           te(2,3,3) = te(2,3,3)*oneplusdelta
+
+           if (sig .gt. zero) then
+
+              te(2,3,3) = te(2,3,3)*oneplusdelta
+
+           else
+
+              te(2,1,1) = te(2,1,1)*oneplusdelta
+              te(4,1,3) = te(4,1,3)*oneplusdelta
+
+           endif
+        endif
+
+        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+
+     endif
+
+  else
+       
+     if (.not.kill_exi_fringe) then
+        if (fintx .lt. 0) fintx = fint
+        corr = (h_k + h_k) * hgap * fint
+        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+     endif
+
   endif
   
   !---- Apply tilt.
@@ -6214,6 +6380,7 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
      orbit(2) = orbit(2)/oneplusdelta
      orbit(4) = orbit(4)/oneplusdelta
      sk1 = sk1/oneplusdelta
+     orbit(6) = 0
 
      call qdbody(fsec,ftrk,tilt,sk1,orbit,dl,ek,re,te)
 

From 70665396d165e467d673277584343d8ddfaffbe1 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Mon, 14 Aug 2023 16:01:51 +0200
Subject: [PATCH 05/14] first failed test for dipoles

---
 src/twiss.f90 | 157 +++++++++++++++++++++++++-------------------------
 1 file changed, 79 insertions(+), 78 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 2f97221f9..9fc4f4964 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -3928,120 +3928,120 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
   !---- Body of the dipole.
   !---- Get map for body section
 
-!  if (exact_expansion) then
+  if (exact_expansion) then
      
-!     pt = orbit(6)
-!     oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-!     orbit(2) = orbit(2)/oneplusdelta
-!     orbit(4) = orbit(4)/oneplusdelta
-!     sk1 = sk1/oneplusdelta
-!     sk2 = sk2/oneplusdelta
-!     orbit(6) = 0
-
-!     call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
-
-!     orbit(2) = orbit(2)*oneplusdelta
-!     orbit(4) = orbit(4)*oneplusdelta
-!     sk1 = sk1*oneplusdelta
-!     sk2 = sk2*oneplusdelta
-!     orbit(6) = pt
+     pt = orbit(6)
+     oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+     orbit(2) = orbit(2)/oneplusdelta
+     orbit(4) = orbit(4)/oneplusdelta
+     sk1 = sk1/oneplusdelta
+     sk2 = sk2/oneplusdelta
+     orbit(6) = 0
+
+     call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
+
+     orbit(2) = orbit(2)*oneplusdelta
+     orbit(4) = orbit(4)*oneplusdelta
+     sk1 = sk1*oneplusdelta
+     sk2 = sk2*oneplusdelta
+     orbit(6) = pt
 
      !----First order terms
 
-!     re(1,2) = re(1,2)/oneplusdelta
-!     re(2,1) = re(2,1)*oneplusdelta
-!     re(2,6) = re(2,6)/oneplusdelta
-!     re(3,4) = re(3,4)/oneplusdelta
-!     re(5,2) = re(5,2)*oneplusdelta
+     re(1,2) = re(1,2)/oneplusdelta
+     re(2,1) = re(2,1)*oneplusdelta
+     re(2,6) = re(2,6)/oneplusdelta
+     re(3,4) = re(3,4)/oneplusdelta
+     re(5,2) = re(5,2)*oneplusdelta
 
-!     ek(2) = ek(2)*oneplusdelta
+     ek(2) = ek(2)*oneplusdelta
 
      !----Second order terms
 
-!     if (fsec) then
+     if (fsec) then
 
-!        te(1,2,2) = te(1,2,2)/oneplusdelta**2
-!        te(1,2,6) = te(1,2,6)/oneplusdelta**2
-!        te(2,1,1) = te(2,1,1)*oneplusdelta
-!        te(2,2,2) = te(2,2,2)/oneplusdelta
-!        te(2,1,6) = te(2,1,6)*oneplusdelta
-!        te(2,6,6) = te(2,6,6)*oneplusdelta
-!        te(5,1,2) = te(5,1,2)/oneplusdelta
-!        te(5,2,2) = te(5,2,2)/oneplusdelta**2
-!        te(5,2,6) = te(5,2,6)/oneplusdelta
+        te(1,2,2) = te(1,2,2)/oneplusdelta**2
+        te(1,2,6) = te(1,2,6)/oneplusdelta**2
+        te(2,1,1) = te(2,1,1)*oneplusdelta
+        te(2,2,2) = te(2,2,2)/oneplusdelta
+        te(2,1,6) = te(2,1,6)*oneplusdelta
+        te(2,6,6) = te(2,6,6)*oneplusdelta
+        te(5,1,2) = te(5,1,2)/oneplusdelta
+        te(5,2,2) = te(5,2,2)/oneplusdelta**2
+        te(5,2,6) = te(5,2,6)/oneplusdelta
 
         !----Mixed terms
 
-!        te(1,3,4) = te(1,3,4)/oneplusdelta
-!        te(1,4,4) = te(1,4,4)/oneplusdelta**2
-!        te(2,3,3) = te(2,3,3)*oneplusdelta
-!        te(2,4,4) = te(2,4,4)/oneplusdelta
-!        te(3,1,4) = te(3,1,4)/oneplusdelta
-!        te(3,2,3) = te(3,2,3)/oneplusdelta
-!        te(3,2,4) = te(3,2,4)/oneplusdelta**2
- !       te(3,4,6) = te(3,4,6)/oneplusdelta
-!        te(4,1,3) = te(4,1,3)*oneplusdelta
-!        te(4,2,4) = te(4,2,4)/oneplusdelta
-!        te(4,3,6) = te(4,3,6)*oneplusdelta
-!        te(5,3,4) = te(5,3,4)/oneplusdelta
-!        te(5,4,4) = te(5,4,4)/oneplusdelta**2
+        te(1,3,4) = te(1,3,4)/oneplusdelta
+        te(1,4,4) = te(1,4,4)/oneplusdelta**2
+        te(2,3,3) = te(2,3,3)*oneplusdelta
+        te(2,4,4) = te(2,4,4)/oneplusdelta
+        te(3,1,4) = te(3,1,4)/oneplusdelta
+        te(3,2,3) = te(3,2,3)/oneplusdelta
+        te(3,2,4) = te(3,2,4)/oneplusdelta**2
+        te(3,4,6) = te(3,4,6)/oneplusdelta
+        te(4,1,3) = te(4,1,3)*oneplusdelta
+        te(4,2,4) = te(4,2,4)/oneplusdelta
+        te(4,3,6) = te(4,3,6)*oneplusdelta
+        te(5,3,4) = te(5,3,4)/oneplusdelta
+        te(5,4,4) = te(5,4,4)/oneplusdelta**2
 
-!     endif  
+     endif  
 
-!  else
+  else
 
      call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
 
-!  endif
+  endif
   
   !---- Get map for entrance fringe field and concatenate
-!  if (exact_expansion) then
+  if (exact_expansion) then
      
-!     if(.not.kill_ent_fringe) then
-!        corr = (h_k + h_k) * hgap * fint
-!        pt = orbit(6)
-!        oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-!        orbit(2) = orbit(2)/oneplusdelta
-!        orbit(4) = orbit(4)/oneplusdelta
-!        sk1 = sk1/oneplusdelta
-!        orbit(6) = 0
+     if(.not.kill_ent_fringe) then
+        corr = (h_k + h_k) * hgap * fint
+        pt = orbit(6)
+        oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+        orbit(2) = orbit(2)/oneplusdelta
+        orbit(4) = orbit(4)/oneplusdelta
+        sk1 = sk1/oneplusdelta
+        orbit(6) = 0
         
-!        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
 
-!        orbit(2) = orbit(2)*oneplusdelta
-!        orbit(4) = orbit(4)*oneplusdelta
-!        sk1 = sk1*oneplusdelta
-!        orbit(6) = pt
+        orbit(2) = orbit(2)*oneplusdelta
+        orbit(4) = orbit(4)*oneplusdelta
+        sk1 = sk1*oneplusdelta
+        orbit(6) = pt
 
         !----First order terms
 
-!        re(2,1) = re(2,1)*oneplusdelta
-!        re(4,3) = re(4,3)*oneplusdelta
+        re(2,1) = re(2,1)*oneplusdelta
+        re(4,3) = re(4,3)*oneplusdelta
 
         !----Second order terms
 
-!        if (fsec) then
+        if (fsec) then
 
-!           te(2,1,1) = te(2,1,1)*oneplusdelta
-!           te(2,3,3) = te(2,3,3)*oneplusdelta
+           te(2,1,1) = te(2,1,1)*oneplusdelta
+           te(2,3,3) = te(2,3,3)*oneplusdelta
 
-!           if (sig .gt. zero) then
+           if (sig .gt. zero) then
 
-!              te(2,3,3) = te(2,3,3)*oneplusdelta
+              te(2,3,3) = te(2,3,3)*oneplusdelta
 
-!           else
+           else
 
-!              te(2,1,1) = te(2,1,1)*oneplusdelta
-!              te(4,1,3) = te(4,1,3)*oneplusdelta
+              te(2,1,1) = te(2,1,1)*oneplusdelta
+              te(4,1,3) = te(4,1,3)*oneplusdelta
 
-!           endif
-!        endif
+           endif
+        endif
 
-!        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
 
-!     endif
+     endif
 
-!  else
+  else
        
      if (.not.kill_ent_fringe) then
         corr = (h_k + h_k) * hgap * fint
@@ -4049,7 +4049,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
         call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
      endif
 
-!  endif
+  endif
   
   !---- Get map for exit fringe fields and concatenate
 
@@ -4063,6 +4063,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
         orbit(2) = orbit(2)/oneplusdelta
         orbit(4) = orbit(4)/oneplusdelta
         sk1 = sk1/oneplusdelta
+        orbit(6) = 0
         
         call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
 

From 9f8957473af8740a1f4206dede8c5e15079a92db Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Thu, 17 Aug 2023 10:43:24 +0200
Subject: [PATCH 06/14] test + correction dip/quad/sext

---
 src/twiss.f90 | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 9fc4f4964..4ba74d490 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -3938,7 +3938,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
      sk2 = sk2/oneplusdelta
      orbit(6) = 0
 
-     call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
+     call tmsect(.true.,dl,h,dh,sk0,sk1,sk2,ek,re,te)
 
      orbit(2) = orbit(2)*oneplusdelta
      orbit(4) = orbit(4)*oneplusdelta
@@ -3990,7 +3990,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
 
   else
 
-     call tmsect(.true.,dl,h,dh,sk1,sk2,ek,re,te)
+     call tmsect(.true.,dl,h,dh,sk0,sk1,sk2,ek,re,te)
 
   endif
   
@@ -4269,7 +4269,7 @@ subroutine tmwire(fsec,ftrk,orbit,fmap,el,ek,re,te)
 endif
 end subroutine
 
-SUBROUTINE tmsect(fsec,el,h,dh,sk1,sk2,ek,re,te)
+SUBROUTINE tmsect(fsec,el,h,dh,sk0,sk1,sk2,ek,re,te)
   use twissbeamfi, only : beta, gamma, dtbyds
   use matrices, only: EYE
   use math_constfi, only : zero, one, two, three, four, six, nine, twelve, fifteen
@@ -4290,7 +4290,7 @@ SUBROUTINE tmsect(fsec,el,h,dh,sk1,sk2,ek,re,te)
   !     te(6,6,6) (double)  second order terms.                          *
   !----------------------------------------------------------------------*
   logical, intent(IN) :: fsec
-  double precision :: el, h, dh, sk1, sk2
+  double precision :: el, h, dh, sk0, sk1, sk2
   double precision :: ek(6), re(6,6), te(6,6,6)
 
   double precision :: bi, bi2, bi2gi2
@@ -4319,7 +4319,12 @@ SUBROUTINE tmsect(fsec,el,h,dh,sk1,sk2,ek,re,te)
   bi2gi2 = one / (beta * gamma) ** 2
 
   !---- Horizontal.
-  xksq = h**2 + sk1
+  if (sk0 .eq. h) then
+     xksq = h**2 + sk1
+  else
+     xksq = h*sk0 + sk1
+  endif
+  !xksq = h**2 + sk1
   xk = sqrt(abs(xksq))
   xkl = xk * el
   xklsq = xksq * el**2

From 55442e1e4afc539b906c0fa237408f05e8ebf954 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Tue, 12 Sep 2023 15:27:04 +0200
Subject: [PATCH 07/14] set back without dipole fix but still with quad and
 sext fix

---
 src/twiss.f90 | 237 ++++++++++++++++++++++++++------------------------
 1 file changed, 121 insertions(+), 116 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 4ba74d490..d8d18fbdc 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -3896,6 +3896,11 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
   !h_k = h_k * (1+ktap) ! tapering applied to actual strength
   
   !---- Change coefficients using DELTAP.
+  !if (sk0 .eq. h) then
+  !   dh = (- h * deltap + bvk * f_errors(0) / el) / (one + deltap)
+  !else
+  !   dh = (- sk0 * deltap + bvk * f_errors(0) / el) / (one + deltap)
+  !endif
   dh = (- h * deltap + bvk * f_errors(0) / el) / (one + deltap) ! dipole term
   sk1 = bvk * (sk1 + f_errors(2) / el) / (one + deltap) ! quad term
   sk2 = bvk * (sk2 + f_errors(4) / el) / (one + deltap) ! sext term
@@ -3928,120 +3933,120 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
   !---- Body of the dipole.
   !---- Get map for body section
 
-  if (exact_expansion) then
+  !if (exact_expansion) then
      
-     pt = orbit(6)
-     oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-     orbit(2) = orbit(2)/oneplusdelta
-     orbit(4) = orbit(4)/oneplusdelta
-     sk1 = sk1/oneplusdelta
-     sk2 = sk2/oneplusdelta
-     orbit(6) = 0
-
-     call tmsect(.true.,dl,h,dh,sk0,sk1,sk2,ek,re,te)
-
-     orbit(2) = orbit(2)*oneplusdelta
-     orbit(4) = orbit(4)*oneplusdelta
-     sk1 = sk1*oneplusdelta
-     sk2 = sk2*oneplusdelta
-     orbit(6) = pt
+  !   pt = orbit(6)
+  !   oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+  !   orbit(2) = orbit(2)/oneplusdelta
+  !   orbit(4) = orbit(4)/oneplusdelta
+  !   sk1 = sk1/oneplusdelta
+  !   sk2 = sk2/oneplusdelta
+  !   orbit(6) = 0
+
+  !   call tmsect(.true.,dl,h,dh,sk0,sk1,sk2,ek,re,te)
+
+  !   orbit(2) = orbit(2)*oneplusdelta
+  !   orbit(4) = orbit(4)*oneplusdelta
+  !   sk1 = sk1*oneplusdelta
+  !   sk2 = sk2*oneplusdelta
+  !   orbit(6) = pt
 
      !----First order terms
 
-     re(1,2) = re(1,2)/oneplusdelta
-     re(2,1) = re(2,1)*oneplusdelta
-     re(2,6) = re(2,6)/oneplusdelta
-     re(3,4) = re(3,4)/oneplusdelta
-     re(5,2) = re(5,2)*oneplusdelta
+  !   re(1,2) = re(1,2)/oneplusdelta
+  !   re(2,1) = re(2,1)*oneplusdelta
+  !   re(2,6) = re(2,6)/oneplusdelta
+  !   re(3,4) = re(3,4)/oneplusdelta
+  !   re(5,2) = re(5,2)*oneplusdelta
 
-     ek(2) = ek(2)*oneplusdelta
+  !   ek(2) = ek(2)*oneplusdelta
 
      !----Second order terms
 
-     if (fsec) then
+  !   if (fsec) then
 
-        te(1,2,2) = te(1,2,2)/oneplusdelta**2
-        te(1,2,6) = te(1,2,6)/oneplusdelta**2
-        te(2,1,1) = te(2,1,1)*oneplusdelta
-        te(2,2,2) = te(2,2,2)/oneplusdelta
-        te(2,1,6) = te(2,1,6)*oneplusdelta
-        te(2,6,6) = te(2,6,6)*oneplusdelta
-        te(5,1,2) = te(5,1,2)/oneplusdelta
-        te(5,2,2) = te(5,2,2)/oneplusdelta**2
-        te(5,2,6) = te(5,2,6)/oneplusdelta
+  !      te(1,2,2) = te(1,2,2)/oneplusdelta**2
+  !      te(1,2,6) = te(1,2,6)/oneplusdelta**2
+  !      te(2,1,1) = te(2,1,1)*oneplusdelta
+  !      te(2,2,2) = te(2,2,2)/oneplusdelta
+  !      te(2,1,6) = te(2,1,6)*oneplusdelta
+  !      te(2,6,6) = te(2,6,6)*oneplusdelta
+  !      te(5,1,2) = te(5,1,2)/oneplusdelta
+  !      te(5,2,2) = te(5,2,2)/oneplusdelta**2
+  !      te(5,2,6) = te(5,2,6)/oneplusdelta
 
         !----Mixed terms
 
-        te(1,3,4) = te(1,3,4)/oneplusdelta
-        te(1,4,4) = te(1,4,4)/oneplusdelta**2
-        te(2,3,3) = te(2,3,3)*oneplusdelta
-        te(2,4,4) = te(2,4,4)/oneplusdelta
-        te(3,1,4) = te(3,1,4)/oneplusdelta
-        te(3,2,3) = te(3,2,3)/oneplusdelta
-        te(3,2,4) = te(3,2,4)/oneplusdelta**2
-        te(3,4,6) = te(3,4,6)/oneplusdelta
-        te(4,1,3) = te(4,1,3)*oneplusdelta
-        te(4,2,4) = te(4,2,4)/oneplusdelta
-        te(4,3,6) = te(4,3,6)*oneplusdelta
-        te(5,3,4) = te(5,3,4)/oneplusdelta
-        te(5,4,4) = te(5,4,4)/oneplusdelta**2
-
-     endif  
+  !      te(1,3,4) = te(1,3,4)/oneplusdelta
+  !      te(1,4,4) = te(1,4,4)/oneplusdelta**2
+  !      te(2,3,3) = te(2,3,3)*oneplusdelta
+  !      te(2,4,4) = te(2,4,4)/oneplusdelta
+  !      te(3,1,4) = te(3,1,4)/oneplusdelta
+  !      te(3,2,3) = te(3,2,3)/oneplusdelta
+  !      te(3,2,4) = te(3,2,4)/oneplusdelta**2
+  !      te(3,4,6) = te(3,4,6)/oneplusdelta
+  !      te(4,1,3) = te(4,1,3)*oneplusdelta
+  !      te(4,2,4) = te(4,2,4)/oneplusdelta
+  !      te(4,3,6) = te(4,3,6)*oneplusdelta
+  !      te(5,3,4) = te(5,3,4)/oneplusdelta
+  !      te(5,4,4) = te(5,4,4)/oneplusdelta**2
+
+  !   endif  
 
-  else
+  !else
 
      call tmsect(.true.,dl,h,dh,sk0,sk1,sk2,ek,re,te)
 
-  endif
+  !endif
   
   !---- Get map for entrance fringe field and concatenate
-  if (exact_expansion) then
+  !if (exact_expansion) then
      
-     if(.not.kill_ent_fringe) then
-        corr = (h_k + h_k) * hgap * fint
-        pt = orbit(6)
-        oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-        orbit(2) = orbit(2)/oneplusdelta
-        orbit(4) = orbit(4)/oneplusdelta
-        sk1 = sk1/oneplusdelta
-        orbit(6) = 0
+  !   if(.not.kill_ent_fringe) then
+  !      corr = (h_k + h_k) * hgap * fint
+  !      pt = orbit(6)
+  !      oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+  !      orbit(2) = orbit(2)/oneplusdelta
+  !      orbit(4) = orbit(4)/oneplusdelta
+  !      sk1 = sk1/oneplusdelta
+  !      orbit(6) = 0
         
-        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+  !      call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
 
-        orbit(2) = orbit(2)*oneplusdelta
-        orbit(4) = orbit(4)*oneplusdelta
-        sk1 = sk1*oneplusdelta
-        orbit(6) = pt
+  !      orbit(2) = orbit(2)*oneplusdelta
+  !      orbit(4) = orbit(4)*oneplusdelta
+  !      sk1 = sk1*oneplusdelta
+  !      orbit(6) = pt
 
         !----First order terms
 
-        re(2,1) = re(2,1)*oneplusdelta
-        re(4,3) = re(4,3)*oneplusdelta
+  !      re(2,1) = re(2,1)*oneplusdelta
+  !      re(4,3) = re(4,3)*oneplusdelta
 
         !----Second order terms
 
-        if (fsec) then
+  !      if (fsec) then
 
-           te(2,1,1) = te(2,1,1)*oneplusdelta
-           te(2,3,3) = te(2,3,3)*oneplusdelta
+  !         te(2,1,1) = te(2,1,1)*oneplusdelta
+  !         te(2,3,3) = te(2,3,3)*oneplusdelta
 
-           if (sig .gt. zero) then
+  !         if (sig .gt. zero) then
 
-              te(2,3,3) = te(2,3,3)*oneplusdelta
+  !            te(2,3,3) = te(2,3,3)*oneplusdelta
 
-           else
+  !         else
 
-              te(2,1,1) = te(2,1,1)*oneplusdelta
-              te(4,1,3) = te(4,1,3)*oneplusdelta
+  !            te(2,1,1) = te(2,1,1)*oneplusdelta
+  !            te(4,1,3) = te(4,1,3)*oneplusdelta
 
-           endif
-        endif
+  !         endif
+  !      endif
 
-        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+  !      call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
 
-     endif
+  !   endif
 
-  else
+  !else
        
      if (.not.kill_ent_fringe) then
         corr = (h_k + h_k) * hgap * fint
@@ -4049,58 +4054,58 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
         call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
      endif
 
-  endif
+  !endif
   
   !---- Get map for exit fringe fields and concatenate
 
-  if (exact_expansion) then
+  !if (exact_expansion) then
      
-     if(.not.kill_exi_fringe) then
-        if (fintx .lt. 0) fintx = fint
-        corr = (h_k + h_k) * hgap * fint
-        pt = orbit(6)
-        oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-        orbit(2) = orbit(2)/oneplusdelta
-        orbit(4) = orbit(4)/oneplusdelta
-        sk1 = sk1/oneplusdelta
-        orbit(6) = 0
+  !   if(.not.kill_exi_fringe) then
+  !      if (fintx .lt. 0) fintx = fint
+  !      corr = (h_k + h_k) * hgap * fint
+  !      pt = orbit(6)
+  !      oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
+  !      orbit(2) = orbit(2)/oneplusdelta
+  !      orbit(4) = orbit(4)/oneplusdelta
+  !      sk1 = sk1/oneplusdelta
+  !      orbit(6) = 0
         
-        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+  !      call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
 
-        orbit(2) = orbit(2)*oneplusdelta
-        orbit(4) = orbit(4)*oneplusdelta
-        sk1 = sk1*oneplusdelta
-        orbit(6) = pt
+  !      orbit(2) = orbit(2)*oneplusdelta
+  !      orbit(4) = orbit(4)*oneplusdelta
+  !      sk1 = sk1*oneplusdelta
+  !      orbit(6) = pt
 
         !----First order terms
 
-        re(2,1) = re(2,1)*oneplusdelta
-        re(4,3) = re(4,3)*oneplusdelta
+  !      re(2,1) = re(2,1)*oneplusdelta
+  !      re(4,3) = re(4,3)*oneplusdelta
 
         !----Second order terms
 
-        if (fsec) then
+  !      if (fsec) then
 
-           te(2,1,1) = te(2,1,1)*oneplusdelta
-           te(2,3,3) = te(2,3,3)*oneplusdelta
+  !         te(2,1,1) = te(2,1,1)*oneplusdelta
+  !         te(2,3,3) = te(2,3,3)*oneplusdelta
 
-           if (sig .gt. zero) then
+  !         if (sig .gt. zero) then
 
-              te(2,3,3) = te(2,3,3)*oneplusdelta
+  !            te(2,3,3) = te(2,3,3)*oneplusdelta
 
-           else
+  !         else
 
-              te(2,1,1) = te(2,1,1)*oneplusdelta
-              te(4,1,3) = te(4,1,3)*oneplusdelta
+  !            te(2,1,1) = te(2,1,1)*oneplusdelta
+  !            te(4,1,3) = te(4,1,3)*oneplusdelta
 
-           endif
-        endif
+  !         endif
+  !      endif
 
-        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+  !      call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
 
-     endif
+  !   endif
 
-  else
+  !else
        
      if (.not.kill_exi_fringe) then
         if (fintx .lt. 0) fintx = fint
@@ -4109,7 +4114,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
         call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
      endif
 
-  endif
+  !endif
   
   !---- Apply tilt.
   if (tilt .ne. zero) then
@@ -4319,12 +4324,12 @@ SUBROUTINE tmsect(fsec,el,h,dh,sk0,sk1,sk2,ek,re,te)
   bi2gi2 = one / (beta * gamma) ** 2
 
   !---- Horizontal.
-  if (sk0 .eq. h) then
-     xksq = h**2 + sk1
-  else
-     xksq = h*sk0 + sk1
-  endif
-  !xksq = h**2 + sk1
+  !if (sk0 .eq. h) then
+  !   xksq = h**2 + sk1
+  !else
+  !   xksq = h*sk0 + sk1
+  !endif
+  xksq = h**2 + sk1
   xk = sqrt(abs(xksq))
   xkl = xk * el
   xklsq = xksq * el**2

From 3c27cf07a52f8237df7c4c7ba733c32c3d48f581 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Tue, 23 Jan 2024 16:41:21 +0100
Subject: [PATCH 08/14] add tests

---
 tests/test-twiss-exactquad/quad.mad           | 73 +++++++++++++++++++
 .../test-twiss-exactktap.madx                 | 60 +++++++++++++++
 .../test-twiss-exactquad.madx                 | 50 +++++++++++++
 3 files changed, 183 insertions(+)
 create mode 100644 tests/test-twiss-exactquad/quad.mad
 create mode 100644 tests/test-twiss-exactquad/test-twiss-exactktap.madx
 create mode 100644 tests/test-twiss-exactquad/test-twiss-exactquad.madx

diff --git a/tests/test-twiss-exactquad/quad.mad b/tests/test-twiss-exactquad/quad.mad
new file mode 100644
index 000000000..a1bb87897
--- /dev/null
+++ b/tests/test-twiss-exactquad/quad.mad
@@ -0,0 +1,73 @@
+local beam, sequence, twiss, track, beta0, damap in MAD
+local sqrt in MAD.gmath
+local quadrupole in MAD.element
+
+
+local function drift_adj (elm, m, l) -- [KICKPATH] drift adjustment          -- checked
+  m.atdebug(elm, m, l, 'drift_adj:0')
+
+  local T in m
+  local _beta = 1/m.beam.beta
+
+  for i=1,m.npar do
+    local x, px, y, py, t, pt, beam in m[i]
+    local _beta = beam and 1/beam.beta or _beta
+    local  l_pz = l/sqrt(1 + (2*_beta)*pt + pt^2 - px^2 - py^2)
+
+    m[i].x = x + px*(l_pz-l)
+    m[i].y = y + py*(l_pz-l)
+    m[i].t = t - l_pz*(_beta+pt) + (1-T)*(l*_beta)
+  end
+
+  m.atdebug(elm, m, l, 'drift_adj:1')
+end
+
+
+
+local q1 = quadrupole 'mq1' {l=1, k1=1.2}
+local q2 = quadrupole 'mq2' {l=1, k1=-1.3}
+
+local bb = beam {particle="positron", pc=0.0001}
+
+local ss = sequence 'ss' {
+    beam =  bb,
+    q1 {at=0.5},
+    q2 {at=1.5}
+}
+local x0,y0, px0,py0,pt0 =1e-3, 1e-3,1e-3,1e-3, 0.1
+
+local mtbl, mflw
+
+mtbl, mflw = twiss {
+       sequence=ss,
+       X0=beta0 {x=x0,y=y0,px=px0,py=py0,pt=pt0, beta11=1, beta22=1},
+       fringe=0,
+       model="DKD",
+       method=8,
+       nslice=80
+   }
+
+print(mtbl.mu1[3],mtbl.mu2[3])
+
+--[[
+tw,m,dr=loadfile("quad.mad")()
+tw.x:print()
+tw.beta11:print()
+tw.mu1:print()
+]]
+
+
+local m={
+    --(damap {xy=2}):set0{x0, y0, px0, py0, 0, pt0},
+    (damap {xy=2}):set0{0, 0, 0, 0, 0, 0},
+    beam=bb, T=0, npar=1, atdebug=\->()
+}
+
+return mtbl,m,drift_adj
+
+
+
+
+
+
+
diff --git a/tests/test-twiss-exactquad/test-twiss-exactktap.madx b/tests/test-twiss-exactquad/test-twiss-exactktap.madx
new file mode 100644
index 000000000..2ad41c0ae
--- /dev/null
+++ b/tests/test-twiss-exactquad/test-twiss-exactktap.madx
@@ -0,0 +1,60 @@
+q1: quadrupole, l=1, k1=1.2;
+q2: quadrupole, l=1, k1=-1.3;
+
+ss: sequence,l=1;
+q1,at=0.5;
+q2,at=1.5;
+endsequence;
+
+!beam,pc=0.000001;
+beam,pc=100;
+use,sequence=ss;
+
+pt0=0.1;
+deltas=sqrt(pt0^2 + 2*pt0/beam->beta+1)-1;
+x0=1e-3; y0=1e-3;
+px0=1e-3; py0=1e-3;
+p0c=beam->pc;
+psc=p0c*(1+deltas);
+es=sqrt(beam->mass**2 + psc**2);
+
+
+twiss,betx=1,bety=1,x=x0, y=y0, px=px0, py=py0;
+qx1=table(twiss,ss$end,mux);qy1=table(twiss,ss$end,muy);
+betx1=table(twiss,ss$end,betx); bety1=table(twiss,ss$end,bety);
+x1=table(twiss,ss$end,x);
+y1=table(twiss,ss$end,y);
+xp1=table(twiss,ss$end,px);
+yp1=table(twiss,ss$end,py);
+print,text="trackend";
+
+q1,ktap=deltas;
+q2,ktap=deltas;
+
+twiss,betx=1/(1+deltas),bety=1/(1+deltas),pt=pt0,x=x0, y=y0, px=px0*(1+deltas), py=py0*(1+deltas);
+qx2=table(twiss,ss$end,mux);qy2=table(twiss,ss$end,muy);
+betx2=table(twiss,ss$end,betx)/(1+deltas); bety2=table(twiss,ss$end,bety)/(1+deltas);
+x2=table(twiss,ss$end,x);
+y2=table(twiss,ss$end,y);
+xp2=table(twiss,ss$end,px)/(1+deltas);
+yp2=table(twiss,ss$end,py)/(1+deltas);
+print,text="trackend";
+
+twiss,betx=1/(1+deltas),bety=1/(1+deltas),pt=pt0,x=x0, y=y0, px=px0*(1+deltas), py=py0*(1+deltas),exact;
+qx3=table(twiss,ss$end,mux);qy3=table(twiss,ss$end,muy);
+betx3=table(twiss,ss$end,betx)/(1+deltas); bety3=table(twiss,ss$end,bety)/(1+deltas);
+x3=table(twiss,ss$end,x);
+y3=table(twiss,ss$end,y);
+xp3=table(twiss,ss$end,px)/(1+deltas);
+yp3=table(twiss,ss$end,py)/(1+deltas);
+
+value,deltas;
+value,qx1,qx2,qx3;
+value,qy1,qy2,qy3;
+value,betx1,betx2,betx3;
+value,bety1,bety2,bety3;
+value,x1,x2,x3;
+value,y1,y2,y3;
+value,xp1,xp2,xp3;
+value,yp1,yp2,yp3;
+value,beam->beta;
diff --git a/tests/test-twiss-exactquad/test-twiss-exactquad.madx b/tests/test-twiss-exactquad/test-twiss-exactquad.madx
new file mode 100644
index 000000000..2119ff4d4
--- /dev/null
+++ b/tests/test-twiss-exactquad/test-twiss-exactquad.madx
@@ -0,0 +1,50 @@
+q1: quadrupole, l=1, k1=1.2;
+q2: quadrupole, l=1, k1=-1.3;
+
+ss: sequence,l=1;
+q1,at=0.5;
+q2,at=1.5;
+endsequence;
+
+beam,pc=0.0001;
+use,sequence=ss;
+
+x0=1e-3; y0=1e-3;
+px0=1e-3; py0=1e-3;
+pt0=0.1;
+deltas=sqrt(pt0^2 + 2*pt0/beam->beta+1)-1;
+
+p0c=beam->pc;
+psc=p0c*(1+deltas);
+es=sqrt(beam->mass**2 + psc**2);
+
+
+twiss,betx=1,bety=1,deltap=deltas,x=x0, y=y0, px=px0/(1+deltas), py=py0/(1+deltas);
+qx1=table(twiss,ss$end,mux);qy1=table(twiss,ss$end,muy);
+x1=table(twiss,ss$end,x);
+y1=table(twiss,ss$end,y);
+px1=table(twiss,ss$end,px)*psc;
+py1=table(twiss,ss$end,py)*psc;
+
+
+twiss,betx=1,bety=1,pt=pt0,x=x0, y=y0, px=px0, py=py0;
+qx2=table(twiss,ss$end,mux);qy2=table(twiss,ss$end,muy);
+x2=table(twiss,ss$end,x);
+y2=table(twiss,ss$end,y);
+px2=table(twiss,ss$end,px)*p0c;
+py2=table(twiss,ss$end,py)*p0c;
+
+twiss,betx=1,bety=1,pt=pt0,x=x0, y=y0, px=px0, py=py0,exact;
+qx3=table(twiss,ss$end,mux);qy3=table(twiss,ss$end,muy);
+x3=table(twiss,ss$end,x);
+y3=table(twiss,ss$end,y);
+px3=table(twiss,ss$end,px)*p0c;
+py3=table(twiss,ss$end,py)*p0c;
+
+value,beam->beta;
+value,qx1,qx2,qx3;
+value,qy1,qy2,qy3;
+value,x1,x2,x3;
+value,y1,y2,y3;
+value,px1,px2,px3;
+value,py1,py2,py3;

From a83065e6ea78891dc6ba43f4027201497b3d8320 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Tue, 23 Jan 2024 17:06:28 +0100
Subject: [PATCH 09/14] small corrections

---
 src/twiss.f90 | 82 +++++++++++----------------------------------------
 1 file changed, 17 insertions(+), 65 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index d8d18fbdc..1d01eefaa 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -4109,8 +4109,8 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
        
      if (.not.kill_exi_fringe) then
         if (fintx .lt. 0) fintx = fint
-        corr = (h_k + h_k) * hgap * fint
-        call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
+        corr = (h_k + h_k) * hgap * fintx
+        call tmfrng(.true.,h_k,sk1,e2,h2,-one,corr,rw,tw)
         call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
      endif
 
@@ -6435,16 +6435,14 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
 
      endif
 
-     if (fcentre) return
-
   else
 
      call qdbody(fsec,ftrk,tilt,sk1,orbit,dl,ek,re,te)
 
-     if (fcentre) return
-
   endif
 
+  if (fcentre) return
+
   !---- Half radiation effect at exit.
   if (radiate .and. ftrk) then
      rfac = (arad * gamma**3 * sk1**2 * el / three) * (orbit(1)**2 + orbit(3)**2)
@@ -6497,32 +6495,16 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
   double precision :: cx, sx, cy, sy, biby4
   double precision :: x,px,y,py,t,pt,deltaplusone,nk1
 
-
-  !if (exact_expansion) then
-  !   x= orbit(1)
-  !   px= orbit(2)
-  !   y= orbit(3)
-  !   py= orbit(4)
-  !   t= orbit(5)
-  !   pt= orbit(6)
-  !   deltaplusone=sqrt(pt**2+2*pt/beta+1)
-  !   nk1=sk1/deltaplusone
-  !else
-  !   nk1= sk1
-  !endif
-
-  nk1 = sk1
-
   !---- Set up c's and s's.
-  qk = sqrt(abs(nk1))
+  qk = sqrt(abs(sk1))
   qkl = qk * el
   if (abs(qkl) .lt. ten3m) then
-     qkl2 = nk1 * el**2
+     qkl2 = sk1 * el**2
      cx = (one - qkl2 / two)
      sx = (one - qkl2 / six) * el
      cy = (one + qkl2 / two)
      sy = (one + qkl2 / six) * el
-  else if (nk1 .gt. zero) then
+  else if (sk1 .gt. zero) then
      cx = cos(qkl)
      sx = sin(qkl) / qk
      cy = cosh(qkl)
@@ -6533,83 +6515,53 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
      cy = cos(qkl)
      sy = sin(qkl) / qk
   endif
-
-  !if (exact_expansion) then
-  !---- First-order terms.
-  !re(1,1) = cx
-  !re(1,2) = sx/deltaplusone
-  !re(2,1) = - nk1 * sx * deltaplusone
-  !re(2,2) = cx
-  !re(3,3) = cy
-  !re(3,4) = sy /deltaplusone
-  !re(4,3) = + nk1 * sy *deltaplusone
-  !re(4,4) = cy
-  !re(5,6) = el/(beta*gamma)**2
-
-  !ek(5) = el*dtbyds ! to be checked
-
-  !else
   
   !---- First-order terms.
   re(1,1) = cx
   re(1,2) = sx
-  re(2,1) = - nk1 * sx
+  re(2,1) = - sk1 * sx
   re(2,2) = cx
   re(3,3) = cy
   re(3,4) = sy
-  re(4,3) = + nk1 * sy
+  re(4,3) = + sk1 * sy
   re(4,4) = cy
   re(5,6) = el/(beta*gamma)**2
 
   ek(5) = el*dtbyds
 
-  !endif
 
   !---- Second-order terms.
   if (fsec) then
      biby4 = one / (four * beta)
 
-     te(1,1,6) = + nk1 * el * sx * biby4
+     te(1,1,6) = + sk1 * el * sx * biby4
      te(1,6,1) = te(1,1,6)
      te(2,2,6) = te(1,1,6)
      te(2,6,2) = te(1,1,6)
      te(1,2,6) = - (sx + el*cx) * biby4
      te(1,6,2) = te(1,2,6)
-     te(2,1,6) = - nk1 * (sx - el*cx) * biby4
+     te(2,1,6) = - sk1 * (sx - el*cx) * biby4
      te(2,6,1) = te(2,1,6)
 
-     te(3,3,6) = - nk1 * el * sy * biby4
+     te(3,3,6) = - sk1 * el * sy * biby4
      te(3,6,3) = te(3,3,6)
      te(4,4,6) = te(3,3,6)
      te(4,6,4) = te(3,3,6)
      te(3,4,6) = - (sy + el*cy) * biby4
      te(3,6,4) = te(3,4,6)
-     te(4,3,6) = + nk1 * (sy - el*cy) * biby4
+     te(4,3,6) = + sk1 * (sy - el*cy) * biby4
      te(4,6,3) = te(4,3,6)
 
-     te(5,1,1) = - nk1 * (el - sx*cx) * biby4
-     te(5,1,2) = + nk1 * sx**2 * biby4
+     te(5,1,1) = - sk1 * (el - sx*cx) * biby4
+     te(5,1,2) = + sk1 * sx**2 * biby4
      te(5,2,1) = te(5,1,2)
      te(5,2,2) = - (el + sx*cx) * biby4
-     te(5,3,3) = + nk1 * (el - sy*cy) * biby4
-     te(5,3,4) = - nk1 * sy**2 * biby4
+     te(5,3,3) = + sk1 * (el - sy*cy) * biby4
+     te(5,3,4) = - sk1 * sy**2 * biby4
      te(5,4,3) = te(5,3,4)
      te(5,4,4) = - (el + sy*cy) * biby4
      te(5,6,6) = (- six * re(5,6)) * biby4
   endif
-
-  !---- Track orbit.
-  !if (exact_expansion) then
-
-  !   orbit(1)=cx*x + sx*px/deltaplusone
-  !   orbit(2)=-nk1 * sx * x*deltaplusone + cx*px
-  !   orbit(3)=cy*y + sy*py/deltaplusone
-  !   orbit(4)=nk1 * sy * y*deltaplusone + cy*py
-  !   orbit(5)=el/(beta*gamma)**2*pt
-  !   re(5,1)=re(5,1) + te(5,1,1)*x + te(5,1,2)*px + te(5,1,3)*y + te(5,1,4)*py + te(5,1,5)*t + te(5,1,6)*pt
-  !else
-  !   if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
-  !endif
   
   if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
   

From 30c69ced9709a15476b1d269eee3e29a670112a6 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Tue, 23 Jan 2024 17:12:08 +0100
Subject: [PATCH 10/14] small corrections

---
 src/twiss.f90 | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 1d01eefaa..4461dbb16 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -6319,7 +6319,7 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
 
   logical :: cplxy
   integer :: i, j, n_ferr, elpar_vl
-  double precision :: ct, st, tmp, biby4
+  double precision :: ct, st, tmp
   double precision :: f_errors(0:maxferr)
   double precision :: tilt, sk1, pt, sk1s, bvk, rfac
 
@@ -6411,8 +6411,6 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
 
      if (fsec) then
 
-        biby4 = one / (four * beta)
-
         te(1,2,6) = te(1,2,6)/oneplusdelta
         te(1,6,2) = te(1,6,2)/oneplusdelta
 

From a38a309f5f55fbfb69023a6275ffb1d3149de9bb Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Tue, 23 Jan 2024 17:16:12 +0100
Subject: [PATCH 11/14] small corrections

---
 src/twiss.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 4461dbb16..450aa1609 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -6491,7 +6491,7 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
 
   double precision :: qk, qkl, qkl2
   double precision :: cx, sx, cy, sy, biby4
-  double precision :: x,px,y,py,t,pt,deltaplusone,nk1
+  double precision :: x,px,y,py,t,pt,deltaplusone
 
   !---- Set up c's and s's.
   qk = sqrt(abs(sk1))

From 2dba4825662298f5d5f0515199fbedadd0f8b483 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Tue, 23 Jan 2024 17:20:36 +0100
Subject: [PATCH 12/14] small corrections

---
 src/twiss.f90 | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 450aa1609..04617fa2f 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -6465,7 +6465,6 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
 end SUBROUTINE tmquad
 
 SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
-  use twisslfi, only: exact_expansion
   use twissbeamfi, only : beta, gamma, dtbyds
   use math_constfi, only : zero, one, two, four, six, ten3m
   implicit none
@@ -6491,7 +6490,6 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
 
   double precision :: qk, qkl, qkl2
   double precision :: cx, sx, cy, sy, biby4
-  double precision :: x,px,y,py,t,pt,deltaplusone
 
   !---- Set up c's and s's.
   qk = sqrt(abs(sk1))
@@ -6527,7 +6525,6 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
 
   ek(5) = el*dtbyds
 
-
   !---- Second-order terms.
   if (fsec) then
      biby4 = one / (four * beta)
@@ -6561,8 +6558,8 @@ SUBROUTINE qdbody(fsec,ftrk,tilt,sk1,orbit,el,ek,re,te)
      te(5,6,6) = (- six * re(5,6)) * biby4
   endif
   
+  !---- Track orbit.
   if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
-  
   !---- Apply tilt.
   if (tilt .ne. zero) call tmtilt(fsec,tilt,ek,re,te)
 
@@ -6917,15 +6914,13 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
         te(5,4,4) = te(5,4,4)/deltaplusone**2
      endif
 
-     if (fcentre) return
-
   else
 
      call sxbody(fsec,ftrk,tilt,sk2,orbit,dl,ek,re,te)
-     if (fcentre) return
 
   endif
 
+  if (fcentre) return
 
 
   !---- Half radiation effects at exit.

From 618b17c1005b247d4ff427d88a45e3a4b7570d1f Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Thu, 25 Jan 2024 16:20:42 +0100
Subject: [PATCH 13/14] tests corrections

---
 src/twiss.f90 | 165 --------------------------------------------------
 1 file changed, 165 deletions(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 04617fa2f..7b68f0cab 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -3933,179 +3933,15 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
   !---- Body of the dipole.
   !---- Get map for body section
 
-  !if (exact_expansion) then
-     
-  !   pt = orbit(6)
-  !   oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-  !   orbit(2) = orbit(2)/oneplusdelta
-  !   orbit(4) = orbit(4)/oneplusdelta
-  !   sk1 = sk1/oneplusdelta
-  !   sk2 = sk2/oneplusdelta
-  !   orbit(6) = 0
-
-  !   call tmsect(.true.,dl,h,dh,sk0,sk1,sk2,ek,re,te)
-
-  !   orbit(2) = orbit(2)*oneplusdelta
-  !   orbit(4) = orbit(4)*oneplusdelta
-  !   sk1 = sk1*oneplusdelta
-  !   sk2 = sk2*oneplusdelta
-  !   orbit(6) = pt
-
-     !----First order terms
-
-  !   re(1,2) = re(1,2)/oneplusdelta
-  !   re(2,1) = re(2,1)*oneplusdelta
-  !   re(2,6) = re(2,6)/oneplusdelta
-  !   re(3,4) = re(3,4)/oneplusdelta
-  !   re(5,2) = re(5,2)*oneplusdelta
-
-  !   ek(2) = ek(2)*oneplusdelta
-
-     !----Second order terms
-
-  !   if (fsec) then
-
-  !      te(1,2,2) = te(1,2,2)/oneplusdelta**2
-  !      te(1,2,6) = te(1,2,6)/oneplusdelta**2
-  !      te(2,1,1) = te(2,1,1)*oneplusdelta
-  !      te(2,2,2) = te(2,2,2)/oneplusdelta
-  !      te(2,1,6) = te(2,1,6)*oneplusdelta
-  !      te(2,6,6) = te(2,6,6)*oneplusdelta
-  !      te(5,1,2) = te(5,1,2)/oneplusdelta
-  !      te(5,2,2) = te(5,2,2)/oneplusdelta**2
-  !      te(5,2,6) = te(5,2,6)/oneplusdelta
-
-        !----Mixed terms
-
-  !      te(1,3,4) = te(1,3,4)/oneplusdelta
-  !      te(1,4,4) = te(1,4,4)/oneplusdelta**2
-  !      te(2,3,3) = te(2,3,3)*oneplusdelta
-  !      te(2,4,4) = te(2,4,4)/oneplusdelta
-  !      te(3,1,4) = te(3,1,4)/oneplusdelta
-  !      te(3,2,3) = te(3,2,3)/oneplusdelta
-  !      te(3,2,4) = te(3,2,4)/oneplusdelta**2
-  !      te(3,4,6) = te(3,4,6)/oneplusdelta
-  !      te(4,1,3) = te(4,1,3)*oneplusdelta
-  !      te(4,2,4) = te(4,2,4)/oneplusdelta
-  !      te(4,3,6) = te(4,3,6)*oneplusdelta
-  !      te(5,3,4) = te(5,3,4)/oneplusdelta
-  !      te(5,4,4) = te(5,4,4)/oneplusdelta**2
-
-  !   endif  
-
-  !else
-
      call tmsect(.true.,dl,h,dh,sk0,sk1,sk2,ek,re,te)
-
-  !endif
-  
-  !---- Get map for entrance fringe field and concatenate
-  !if (exact_expansion) then
-     
-  !   if(.not.kill_ent_fringe) then
-  !      corr = (h_k + h_k) * hgap * fint
-  !      pt = orbit(6)
-  !      oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-  !      orbit(2) = orbit(2)/oneplusdelta
-  !      orbit(4) = orbit(4)/oneplusdelta
-  !      sk1 = sk1/oneplusdelta
-  !      orbit(6) = 0
-        
-  !      call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
-
-  !      orbit(2) = orbit(2)*oneplusdelta
-  !      orbit(4) = orbit(4)*oneplusdelta
-  !      sk1 = sk1*oneplusdelta
-  !      orbit(6) = pt
-
-        !----First order terms
-
-  !      re(2,1) = re(2,1)*oneplusdelta
-  !      re(4,3) = re(4,3)*oneplusdelta
-
-        !----Second order terms
-
-  !      if (fsec) then
-
-  !         te(2,1,1) = te(2,1,1)*oneplusdelta
-  !         te(2,3,3) = te(2,3,3)*oneplusdelta
-
-  !         if (sig .gt. zero) then
-
-  !            te(2,3,3) = te(2,3,3)*oneplusdelta
-
-  !         else
-
-  !            te(2,1,1) = te(2,1,1)*oneplusdelta
-  !            te(4,1,3) = te(4,1,3)*oneplusdelta
-
-  !         endif
-  !      endif
-
-  !      call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
-
-  !   endif
-
-  !else
        
      if (.not.kill_ent_fringe) then
         corr = (h_k + h_k) * hgap * fint
         call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
         call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
      endif
-
-  !endif
   
   !---- Get map for exit fringe fields and concatenate
-
-  !if (exact_expansion) then
-     
-  !   if(.not.kill_exi_fringe) then
-  !      if (fintx .lt. 0) fintx = fint
-  !      corr = (h_k + h_k) * hgap * fint
-  !      pt = orbit(6)
-  !      oneplusdelta = sqrt(pt**2 + 2*pt/bet0 + 1)
-  !      orbit(2) = orbit(2)/oneplusdelta
-  !      orbit(4) = orbit(4)/oneplusdelta
-  !      sk1 = sk1/oneplusdelta
-  !      orbit(6) = 0
-        
-  !      call tmfrng(.true.,h_k,sk1,e1,h1,one,corr,rw,tw)
-
-  !      orbit(2) = orbit(2)*oneplusdelta
-  !      orbit(4) = orbit(4)*oneplusdelta
-  !      sk1 = sk1*oneplusdelta
-  !      orbit(6) = pt
-
-        !----First order terms
-
-  !      re(2,1) = re(2,1)*oneplusdelta
-  !      re(4,3) = re(4,3)*oneplusdelta
-
-        !----Second order terms
-
-  !      if (fsec) then
-
-  !         te(2,1,1) = te(2,1,1)*oneplusdelta
-  !         te(2,3,3) = te(2,3,3)*oneplusdelta
-
-  !         if (sig .gt. zero) then
-
-  !            te(2,3,3) = te(2,3,3)*oneplusdelta
-
-  !         else
-
-  !            te(2,1,1) = te(2,1,1)*oneplusdelta
-  !            te(4,1,3) = te(4,1,3)*oneplusdelta
-
-  !         endif
-  !      endif
-
-  !      call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
-
-  !   endif
-
-  !else
        
      if (.not.kill_exi_fringe) then
         if (fintx .lt. 0) fintx = fint
@@ -4114,7 +3950,6 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
         call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
      endif
 
-  !endif
   
   !---- Apply tilt.
   if (tilt .ne. zero) then

From 6ec7ef64a83e2ef3139c3fe7cf3743afcee354a0 Mon Sep 17 00:00:00 2001
From: Guillaume Simon <gusimon@cern.ch>
Date: Thu, 25 Jan 2024 16:25:29 +0100
Subject: [PATCH 14/14] tests corrections

---
 src/twiss.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/twiss.f90 b/src/twiss.f90
index 7b68f0cab..59fd94d9b 100644
--- a/src/twiss.f90
+++ b/src/twiss.f90
@@ -3947,7 +3947,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
         if (fintx .lt. 0) fintx = fint
         corr = (h_k + h_k) * hgap * fintx
         call tmfrng(.true.,h_k,sk1,e2,h2,-one,corr,rw,tw)
-        call tmcat1(.true.,ek,re,te,ek0,rw,tw,ek,re,te)
+        call tmcat1(.true.,ek0,rw,tw,ek,re,te,ek,re,te)
      endif