Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Make track/twiss able to treat a thick solenoid exactly #1129

Merged
merged 7 commits into from Aug 15, 2022
Merged

Make track/twiss able to treat a thick solenoid exactly #1129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
fed4b43
Create a sinc function
Jul 15, 2022
88f29f7
Compute the exact solution for a uniform solenoid in TWISS
Jul 15, 2022
df05381
Remove debugging lines
Jul 15, 2022
0a26cbe
Add a test case for the exact solenoid solution
Jul 15, 2022
8ebb2e7
Make exact solenoid work with xtilt attribute
Jul 19, 2022
e2124e0
Tracking for solenoid is now exact
Jul 19, 2022
a460524
Merge branch 'master' into exact-solenoid
rdemaria Aug 15, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
36 changes: 19 additions & 17 deletions src/trrun.f90
View file
Open in desktop
Expand Up @@ -3387,6 +3387,11 @@ subroutine trsol(track,ktrack,dxt,dyt)
double precision :: curv, const, rfac
double precision :: beta_sqr, f_damp_t

real(kind(0d0)), external :: sinc
real(kind(0d0)) :: ptr2,si
real(kind(0d0)), dimension(2) :: pk
real(kind(0d0)), dimension(4) :: rps

!---- Get solenoid parameters
elrad = node_value('lrad ')
bvk = node_value('other_bv ')
Expand Down Expand Up @@ -3491,18 +3496,18 @@ subroutine trsol(track,ktrack,dxt,dyt)
pt_ = track(6,i)

! set up constants
onedp = sqrt(one + two*pt_*beti + pt_**2);
pk = [ px_+sk*y_, py_-sk*x_ ]
ptr2 = pk(1)**2 + pk(2)**2
onedp = sqrt(one + two*pt_*beti + pt_**2 - ptr2)

! set up constants
cosTh = cos(two*skl/onedp)
sinTh = sin(two*skl/onedp)
cosTh = cos(skl/onedp)
sinTh = sin(skl/onedp)
omega = sk/onedp;

! Store the kick for radiation calculations
pxf_ = (omega*((cosTh-one)*y_-sinTh*x_)+py_*sinTh+px_*(one+cosTh))/two;
pyf_ = (omega*((one-cosTh)*x_-sinTh*y_)-px_*sinTh+py_*(one+cosTh))/two;
dxt(i) = pxf_ - track(2,i);
dyt(i) = pyf_ - track(4,i);
dxt(i) = sinTh*(cosTh*pk(2)-sinTh*pk(1))
dyt(i) = -sinTh*(cosTh*pk(1)+sinTh*pk(2))

if ((step.eq.1).or.(step.eq.3)) then
if (radiate .and. elrad .gt. zero) then
Expand Down Expand Up @@ -3532,17 +3537,14 @@ subroutine trsol(track,ktrack,dxt,dyt)
endif
endif
else
! total path length traveled by the particle
bet = onedp / (beti + pt_);
length_ = length - half/(onedp**2)*(omega*(sinTh-two*length*omega)*(x_**2+y_**2)+&
two*(one-cosTh)*(px_*x_+py_*y_)-(sinTh/omega+two*length)*(px_**2+py_**2))/four;

! Thick transport
track(1,i) = ((one+cosTh)*x_+sinTh*y_+(px_*sinTh-py_*(cosTh-one))/omega)/two;
track(3,i) = ((one+cosTh)*y_-sinTh*x_+(py_*sinTh+px_*(cosTh-one))/omega)/two;
track(2,i) = pxf_;
track(4,i) = pyf_;
track(5,i) = z_ + length*beti - length_/bet;
si = sinc(skl/onedp)*length/onedp
rps = [ cosTh*x_ + sinTh*y_, cosTh*px_ + sinTh*py_, cosTh*y_ - sinTh*x_, cosTh*py_ - sinTh*px_ ]
track(1,i) = cosTh*rps(1) + si*rps(2)
track(2,i) = cosTh*rps(2) - sk*sinTh*rps(1)
track(3,i) = cosTh*rps(3) + si*rps(4)
track(4,i) = cosTh*rps(4) - sk*sinTh*rps(3)
track(5,i) = z_ + length*(pt_*(pt_+2d0*beti)/gammas**2-ptr2)/(betas*onedp*(1d0+onedp+betas*pt_))
endif
enddo ! step

Expand Down
173 changes: 125 additions & 48 deletions src/twiss.f90
View file
Open in desktop
Expand Up @@ -6742,6 +6742,7 @@ end SUBROUTINE tmsol

SUBROUTINE tmsol0(fsec,ftrk,orbit,fmap,el,ek,re,te)
use twissbeamfi, only : radiate, deltap, beta, gamma, dtbyds, arad
use twisslfi
use math_constfi, only : zero, one, two, three, six
use matrices, only : EYE
implicit none
Expand Down Expand Up @@ -6774,6 +6775,13 @@ SUBROUTINE tmsol0(fsec,ftrk,orbit,fmap,el,ek,re,te)
double precision :: rfac, kx, ky
double precision :: pt, bet0, bet_sqr, f_damp_t

integer :: i,j
real(kind(1d0)) :: beti,bg2iptr2,c2,elpsi1,elpsi3,elpsi5,ps,ps2,ptb,ptr2,s2
real(kind(1d0)), dimension(2) :: pk
real(kind(1d0)), dimension(4) :: vmu,vpi,vtau
real(kind(1d0)), dimension(4,4) :: ms1,mpi
real(kind(1d0)), external :: sinc ! defined in util.f90

bet0 = get_value('beam ','beta ')

beta0 = get_value('probe ','beta ')
Expand Down Expand Up @@ -6802,23 +6810,53 @@ SUBROUTINE tmsol0(fsec,ftrk,orbit,fmap,el,ek,re,te)
bvk = node_value('other_bv ')
sks = sks * bvk

pt = orbit(6)
!---- Set up C's and S's.
sk = sks / two / (one + deltap)
skl = sk * el
if (exact_expansion) then
beti = 1d0/beta
ptb = pt + beti
pk = [ orbit(2)+sk*orbit(3), orbit(4)-sk*orbit(1) ]
ptr2 = pk(1)**2 + pk(2)**2
ps2 = 1d0 + (2*beti + pt)*pt - ptr2
ps = sqrt(ps2)
bg2iptr2 = 1d0/(beta*gamma)**2 + ptr2
elpsi1 = el / ps
elpsi3 = elpsi1 / ps2
elpsi5 = elpsi3 / ps2
skl = sk * elpsi1
sibk = sinc(skl) * elpsi1
c2 = cos(2*skl)
s2 = sin(2*skl)
ms1 = elpsi3*reshape([ &
-sk*s2, c2, sk*c2, s2, &
-sk*sk*c2, -sk*s2, -sk*sk*s2, sk*c2, &
-sk*c2, -s2, -sk*s2, c2, &
sk*sk*s2, -sk*c2, -sk*sk*c2, -sk*s2 &
], [4,4], order=[2,1])
mpi = reshape([ &
sk*sk, 0d0, 0d0, -sk, &
0d0, 1d0, sk, 0d0, &
0d0, sk, sk*sk, 0d0, &
-sk, 0d0, 0d0, 1d0 &
], [4,4])
vmu(1) = elpsi3*(c2*pk(1)+s2*pk(2))
vmu(3) = elpsi3*(c2*pk(2)-s2*pk(1))
vmu(2) = sk*vmu(3)
vmu(4) = -sk*vmu(1)
vpi = [ -sk*pk(2), pk(1), sk*pk(1), pk(2) ]
else
skl = sk * el
sibk = sinc(skl) * el
end if
co = cos(skl)
si = sin(skl)
if (abs(skl) .lt. ten5m) then
sibk = (one - skl**2/six) * el
else
sibk = si/sk
endif

!---- Half radiation effect at entry.
if (radiate .and. ftrk) then
kx = ((sk**2)*orbit(1)-sk*orbit(4))*el;
ky = ((sk**2)*orbit(3)+sk*orbit(2))*el;
rfac = (arad * gamma**3 / three) * (kx**2 + ky**2) / el;
pt = orbit(6);
bet_sqr = (pt*pt + two*pt/bet0 + one) / (one/bet0 + pt)**2;
f_damp_t = sqrt(one + rfac*(rfac - two) / bet_sqr);
orbit(2) = orbit(2) * f_damp_t;
Expand All @@ -6844,41 +6882,68 @@ SUBROUTINE tmsol0(fsec,ftrk,orbit,fmap,el,ek,re,te)
re_s(3,2) = - re_s(1,4)
re_s(4,1) = sk * si**2
re_s(2,3) = - re_s(4,1)
re_s(5,6) = el/(beta*gamma)**2

ek_s(5) = el*dtbyds
if (exact_expansion) then
! re_s(1:4,1:4) will be updated later
ek_s(1:4) = [ sibk*(co*pk(1)+si*pk(2)), si*(co*pk(2)-si*pk(1)), sibk*(co*pk(2)-si*pk(1)), -si*(co*pk(1)+si*pk(2)) ]
ek_s(5) = el*dtbyds + elpsi1 * ( pt*(2*beti+pt)/gamma**2 - ptr2 ) / ( beta**2 * (beti*ps+ptb) )
re_s(1:4,6) = -ptb*vmu
re_s(5,1:4) = -elpsi3*ptb*vpi
re_s(5,6) = elpsi3 * bg2iptr2
else
re_s(5,6) = el/(beta*gamma)**2
ek_s(5) = el*dtbyds
end if

!---- Second-order terms.
if (fsec) then
temp = el * co * si / beta
te_s(1,4,6) = - temp
te_s(3,2,6) = temp
te_s(1,1,6) = temp * sk
te_s(2,2,6) = temp * sk
te_s(3,3,6) = temp * sk
te_s(4,4,6) = temp * sk
te_s(2,3,6) = temp * sk**2
te_s(4,1,6) = - temp * sk**2

temp = el * (co**2 - si**2) / (two * beta)
te_s(1,2,6) = - temp
te_s(3,4,6) = - temp
te_s(1,3,6) = - temp * sk
te_s(2,4,6) = - temp * sk
te_s(3,1,6) = temp * sk
te_s(4,2,6) = temp * sk
te_s(2,1,6) = temp * sk**2
te_s(4,3,6) = temp * sk**2

temp = el / (two * beta)
te_s(5,2,2) = - temp
te_s(5,4,4) = - temp
te_s(5,1,4) = temp * sk
te_s(5,2,3) = - temp * sk
te_s(5,1,1) = - temp * sk**2
te_s(5,3,3) = - temp * sk**2
te_s(5,6,6) = - three * re_s(5,6) / (two * beta)
call tmsymm(te_s)
if (exact_expansion) then
vtau(1) = (3d0*vmu(1))/ps2 - 2d0*sk*elpsi3*elpsi3*(s2*pk(1)-c2*pk(2))
vtau(3) = (3d0*vmu(3))/ps2 - 2d0*sk*elpsi3*elpsi3*(c2*pk(1)+s2*pk(2))
vtau(2) = sk*vtau(3)
vtau(4) = -sk*vtau(1)
te_s(1:4,6,6) = 0.5d0*((ptb*ptb)*vtau - vmu)
te_s(5,1:4,6) = (0.5d0*elpsi5*(3*ptb*ptb-ps2))*vpi
te_s(5,6,1:4) = te_s(5,1:4,6)
do i=1,4
te_s(1:4,i,6) = (-0.5d0*ptb)*(vpi(i)*vtau + ms1(:,i))
te_s(1:4,6,i) = te_s(1:4,i,6)
te_s(5,1:4,i) = (-1.5d0*elpsi5*ptb*vpi(i))*vpi - (0.5d0*elpsi3*ptb)*mpi(:,i)
do j=1,4
te_s(1:4,i,j) = 0.5d0*(ms1(:,i)*vpi(j) + ms1(:,j)*vpi(i) + mpi(i,j)*vmu + (vpi(i)*vpi(j))*vtau)
end do
end do
te_s(5,6,6) = -1.5d0*elpsi5*ptb*bg2iptr2
else
temp = el * co * si / beta
te_s(1,4,6) = - temp
te_s(3,2,6) = temp
te_s(1,1,6) = temp * sk
te_s(2,2,6) = temp * sk
te_s(3,3,6) = temp * sk
te_s(4,4,6) = temp * sk
te_s(2,3,6) = temp * sk**2
te_s(4,1,6) = - temp * sk**2

temp = el * (co**2 - si**2) / (two * beta)
te_s(1,2,6) = - temp
te_s(3,4,6) = - temp
te_s(1,3,6) = - temp * sk
te_s(2,4,6) = - temp * sk
te_s(3,1,6) = temp * sk
te_s(4,2,6) = temp * sk
te_s(2,1,6) = temp * sk**2
te_s(4,3,6) = temp * sk**2

temp = el / (two * beta)
te_s(5,2,2) = - temp
te_s(5,4,4) = - temp
te_s(5,1,4) = temp * sk
te_s(5,2,3) = - temp * sk
te_s(5,1,1) = - temp * sk**2
te_s(5,3,3) = - temp * sk**2
te_s(5,6,6) = - three * re_s(5,6) / (two * beta)
call tmsymm(te_s)
end if
endif


Expand All @@ -6898,9 +6963,19 @@ SUBROUTINE tmsol0(fsec,ftrk,orbit,fmap,el,ek,re,te)
re_t1(5,2) = -pxbeta
call tmtrak(ek_t1,re_t1,te_t1,orbit,orbit)
call tmcat(.true.,re_t1,te_t1,re,te,re,te)
end if


if (exact_expansion) then
orbit(1:4) = matmul(re_s(1:4,1:4),orbit(1:4))
orbit(5) = orbit(5) + ek_s(5)
do i=1,4
re_s(1:4,i) = re_s(1:4,i) + vpi(i)*vmu
end do
else
call tmtrak(ek_s,re_s,te_s,orbit,orbit) ! Calls the normal solenoid
endif

if(abs(xtilt_rad) > ten5m) then
call tmcat(.true.,re_s,te_s,re,te,re,te)

!To tilt it back
Expand All @@ -6914,17 +6989,20 @@ SUBROUTINE tmsol0(fsec,ftrk,orbit,fmap,el,ek,re,te)

call tmtrak(ek_t2,re_t2,te_t1 ,orbit,orbit)
call tmcat(.true.,re_t2,te_t1,re,te,re,te)

else
ek=ek_s
re=re_s
te=te_s
call tmtrak(ek,re,te,orbit,orbit)
endif
else
ek = ek_s
re = re_s
te = te_s
end if
else ! ftrk == .false.
ek=ek_s
re=re_s
te=te_s
if (exact_expansion) then
do i=1,4
re(1:4,i) = re(1:4,i) + vpi(i)*vmu
end do
end if
endif


Expand Down Expand Up @@ -9343,4 +9421,3 @@ SUBROUTINE twcptk_print(re,r0mat, e, f)


end SUBROUTINE twcptk_print

12 changes: 12 additions & 0 deletions src/util.f90
View file
Open in desktop
Expand Up @@ -2247,3 +2247,15 @@ double precision function bips(a, mmax1)
&a*(4576d0/dble(11+mmax1)+7d0*a*(-(624d0/dble(12+mmax1))+(598d0*&
&a)/dble(13+mmax1)-(575d0*a**2)/dble(14+mmax1)))))))))))))/8388608d0;
end function bips

! sin(x)/x, with prper handling of values near zero
function sinc(x)
real(kind(1d0)), intent(in) :: x
real(kind(1d0)) :: sinc
real(kind(1d0)), parameter :: e = 1.5d0*sqrt(epsilon(1d0))
if (abs(x) < e) then
sinc = 1d0 - (x*x)/6d0
else
sinc = sin(x) / x
end if
end function sinc
8 changes: 4 additions & 4 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0001.ref
View file
Open in desktop
@@ -1,10 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0001"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %16s "5.04.00 Linux 64"
@ DATE %08s "15/05/18"
@ TIME %08s "18.31.00"
@ ORIGIN %16s "5.08.01 Linux 64"
@ DATE %08s "19/07/22"
@ TIME %08s "14.39.30"
* NUMBER TURN X PX Y PY T PT S E
$ %d %d %le %le %le %le %le %le %le %le
1 0 0.001 0 0 0 0 -0.001 0 1
1 1 -0.0004179664169 -9.093718502e-05 0 0 -6.495195848e-08 -0.001 0 1
1 1 -0.000417966426 -9.08462474e-05 2.924416989e-20 -4.564552995e-21 -1.054305582e-07 -0.001 0 1
10 changes: 5 additions & 5 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0002.ref
View file
Open in desktop
@@ -1,10 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0002"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %16s "5.04.00 Linux 64"
@ DATE %08s "15/05/18"
@ TIME %08s "18.31.00"
@ ORIGIN %16s "5.08.01 Linux 64"
@ DATE %08s "19/07/22"
@ TIME %08s "14.39.30"
* NUMBER TURN X PX Y PY T PT S E
$ %d %d %le %le %le %le %le %le %le %le
2 0 0.001 0 0 0 0 2.220446049e-16 0 1
2 1 -0.0004161468365 -9.092974268e-05 0 3.388131789e-21 -5.946003867e-08 2.220446049e-16 0 1
2 0 0.001 0 0 0 0 0 0 1
2 1 -0.0004161468456 -9.092974227e-05 -2.773044558e-20 -8.222854172e-21 -1.000000138e-07 0 0 1
8 changes: 4 additions & 4 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0003.ref
View file
Open in desktop
@@ -1,10 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0003"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %16s "5.04.00 Linux 64"
@ DATE %08s "15/05/18"
@ TIME %08s "18.31.00"
@ ORIGIN %16s "5.08.01 Linux 64"
@ DATE %08s "19/07/22"
@ TIME %08s "14.39.30"
* NUMBER TURN X PX Y PY T PT S E
$ %d %d %le %le %le %le %le %le %le %le
3 0 0.001 0 0 0 0 0.001 0 1
3 1 -0.0004143292288 -9.092178557e-05 6.783039842e-20 0 -5.398501379e-08 0.001 0 1
3 1 -0.0004143292379 -9.101270695e-05 0 0 -9.458573659e-08 0.001 0 1
18 changes: 9 additions & 9 deletions tests/test-thick-solenoid/test-thick-solenoid.ref
View file
Open in desktop
@@ -1,9 +1,9 @@

++++++++++++++++++++++++++++++++++++++++++++
+ MAD-X 5.07.00 (64 bit, Linux) +
+ MAD-X 5.08.01 (64 bit, Linux) +
+ Support: mad@cern.ch, http://cern.ch/mad +
+ Release date: 2021.05.03 +
+ Execution date: 2021.12.10 13:04:18 +
+ Release date: 2022.02.25 +
+ Execution date: 2022.07.19 14:39:30 +
++++++++++++++++++++++++++++++++++++++++++++
option, -debug, -echo;

Expand All @@ -16,17 +16,17 @@ one pass is on
1 0 0.001 0
2 0 0.001 0
3 0 0.001 0
1 1 -0.0004179664169 -9.093718502e-05
2 1 -0.0004161468365 -9.092974268e-05
3 1 -0.0004143292288 -9.092178557e-05
1 1 -0.000417966426 -9.08462474e-05
2 1 -0.0004161468456 -9.092974227e-05
3 1 -0.0004143292379 -9.101270695e-05

y py t pt
0 0 0 -0.001
0 0 0 0
0 0 0 0.001
6.769487314e-20 3.388131789e-21 -6.495195493e-08 -0.001
0 0 -5.946003689e-08 0
6.783039842e-20 3.388131789e-21 -5.398501202e-08 0.001
2.924416989e-20 -4.564552995e-21 -1.054305582e-07 -0.001
-2.773044558e-20 -8.222854172e-21 -1.000000138e-07 0
0 0 -9.458573659e-08 0.001

s e
0 0
Expand Down