Drift expanded around the closed orbit instead of the ideal orbit

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 \

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

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()

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], "

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
 

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.
 
 

tests/test-twiss-exact/madx_exact.tfs.cfg

@@ -0,0 +1,3 @@
+44-48   *     skip         # origin, date and time
+*       *     abs=1e-298
+

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

tests/test-twiss-exact/ptc.tfs.cfg

@@ -0,0 +1,2 @@
+85-88   *   skip        # date, version
+* * abs=1e-12