Skip to content

Commit

Permalink
Merge pull request #1129 from jsberg-bnl/exact-solenoid
Browse files Browse the repository at this point in the history 
Make track/twiss able to treat a thick solenoid exactly
  • Loading branch information
@rdemaria
rdemaria committed Aug 15, 2022
2 parents e19ca7b + a460524 commit 602c8f7
Show file tree
Hide file tree
Showing 20 changed files with 1,072 additions and 87 deletions.
36 changes: 19 additions & 17 deletions src/trrun.f90
View file
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
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
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
@@ -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
@@ -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
@@ -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
@@ -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

0 comments on commit 602c8f7

Please sign in to comment.