exact Twiss Riccardo [WIP]

src/twiss.f90

@@ -3890,7 +3890,7 @@ SUBROUTINE tmbend(ftrk,fcentre,orbit,fmap,el,dl,ek,re,te,code)
 
   ! ugly kludge to mimic k0 integration in errors
   if (ktap .ne. 0) then
-     f_errors(0) = f_errors(0) + ktap*angle     
+     f_errors(0) = f_errors(0) + ktap*angle
   endif
 
   !h_k = h_k * (1+ktap) ! tapering applied to actual strength
@@ -4103,6 +4103,7 @@ subroutine tmwire(fsec,ftrk,orbit,fmap,el,ek,re,te)
 end subroutine
 
 SUBROUTINE tmsect(fsec,el,h,dh,sk1,sk2,ek,re,te)
+  use twisslfi, only: exact_expansion
   use twissbeamfi, only : beta, gamma, dtbyds
   use matrices, only: EYE
   use math_constfi, only : zero, one, two, three, four, six, nine, twelve, fifteen
@@ -4141,6 +4142,7 @@ SUBROUTINE tmsect(fsec,el,h,dh,sk1,sk2,ek,re,te)
   double precision, parameter :: s1=one, s2=one/six, s3=one/httwty, s4=one/5040d0
   double precision, parameter :: cg0=one/twty, cg1=5d0/840d0, cg2=21d0/60480d0
   double precision, parameter :: ch0=one/fvty6, ch1=14d0/4032d0, ch2=147d0/443520d0
+  double precision :: newdeltasplusone, ff, newbeta, newgamma, pt
 
   !---- Initialize.
   EK = zero
@@ -6151,7 +6153,6 @@ SUBROUTINE tmquad(fsec,ftrk,fcentre,plot_tilt,orbit,fmap,el,dl,ek,re,te)
   double precision, external :: node_value, get_value
   double precision :: bet0, bet_sqr, f_damp_t
 
-  !double precision :: newbet0, newdeltap, ff
   integer, external :: get_option
   character(len=name_len) el_name
 
@@ -6207,30 +6208,8 @@ 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
-
   call qdbody(fsec,ftrk,tilt,sk1,orbit,dl,ek,re,te)
 
-  ! restore pt
-  !sk1 =sk1/ff
-  !orbit(2)=orbit(2)/ff
-  !orbit(4)=orbit(4)/ff
-  !orbit(6)=pt
-  ! end restore pt
-
-
-
-
-
   if (fcentre) return
 
   !---- Half radiation effect at exit.
@@ -6257,6 +6236,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 +6262,32 @@ 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
 
   !---- 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 +6299,80 @@ 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
-
+  ek(5) = el*dtbyds ! to be checked
+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 (ftrk) call tmtrak(ek,re,te,orbit,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
+
   !---- Apply tilt.
   if (tilt .ne. zero) call tmtilt(fsec,tilt,ek,re,te)
 
@@ -6675,6 +6696,7 @@ SUBROUTINE tmsext(fsec,ftrk,fcentre,orbit,fmap,el,dl,ek,re,te)
 end SUBROUTINE tmsext
 
 SUBROUTINE sxbody(fsec,ftrk,tilt,sk2,orbit,el,ek,re,te)
+  use twisslfi, only: exact_expansion
   use twissbeamfi, only : beta, gamma, dtbyds
   use math_constfi, only : zero, two, three, four
   implicit none
@@ -6699,6 +6721,14 @@ SUBROUTINE sxbody(fsec,ftrk,tilt,sk2,orbit,el,ek,re,te)
   double precision :: orbit(6), ek(6), re(6,6), te(6,6,6)
 
   double precision :: skl, s1, s2, s3, s4
+  double precision :: x,px,y,py,t,pt,deltaplusone
+
+
+  if (exact_expansion) then
+     pt= orbit(6)
+     deltaplusone=sqrt(pt**2+2*pt/beta+1)
+  endif
+
 
   !---- First-order terms.
   re(1,2) = el
@@ -6745,6 +6775,15 @@ SUBROUTINE sxbody(fsec,ftrk,tilt,sk2,orbit,el,ek,re,te)
 
   !---- Track orbit.
   if (ftrk) call tmtrak(ek,re,te,orbit,orbit)
+
+  !if (exact_expansion) then
+  !    ! restore pt, px(deltas), py(deltas) using old deltap
+  !    orbit(6)=pt
+  !    orbit(2)=orbit(2);
+  !    orbit(4)=orbit(4);
+  !endif
+
+
   !---- Apply tilt.
   if (tilt .ne. zero) call tmtilt(fsec,tilt,ek,re,te)
 

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
+
+
+
+
+
+
+

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;