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

Implemented thick solenoid tracking in module TRACK #495

Merged
merged 7 commits into from Nov 7, 2017
Merged

Implemented thick solenoid tracking in module TRACK #495

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
1c7780d
Implemented thick solenoid tracking
Nov 6, 2017
72bf111
Cosmetics...
Nov 6, 2017
8712227
Merge branch 'master' of github.com:MethodicalAcceleratorDesign/MAD-X
Nov 6, 2017
f9d9349
Improved thick solenoid map and added test
Nov 6, 2017
270fd79
Implemented drift transfer map when thick solenoid has strength zero
Nov 6, 2017
18605e4
Thick solenoid: if strength == 0 then it just calls ttdft()
Nov 6, 2017
814ae67
Preparing for thorough test with mkthin
Nov 6, 2017
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
1 change: 1 addition & 0 deletions Makefile_test
View file
Open in desktop
Expand Up @@ -60,6 +60,7 @@ test-dynap \
test-c6t test-c6t-2 test-c6t-3 test-c6t-4 \
test-thick-quad test-thick-quad-2 test-thick-quad-3 \
test-thick-dipole test-thick-dipole-2 test-thick-dipole-3 \
test-thick-solenoid \
test-jacobian test-jacobian-2 test-jacobian-knobs \
test-match test-match-2 test-match-3 test-match-4 \
test-match-5 test-match-6 test-match-7 test-match-8 test-match-9\
Expand Down
149 changes: 95 additions & 54 deletions src/trrun.f90
View file
Open in desktop
Expand Up @@ -469,8 +469,12 @@ subroutine trrun(switch, turns, orbit0, rt, part_id, last_turn, last_pos, &
l_buf(nlm+1) = el
call element_name(el_name,len(el_name))

if ( code.ne.code_drift .and. code.ne.code_quadrupole .and. code.ne.code_sbend &
.and. code.ne.code_matrix .and. el.ne.zero ) then ! rbend missing ?
if ( code.ne.code_drift .and. &
code.ne.code_quadrupole .and. &
code.ne.code_sbend .and. &
code.ne.code_matrix .and. &
code.ne.code_solenoid .and. &
el.ne.zero ) then ! rbend missing ?
!if (.not. (is_drift() .or. is_thin() .or. is_quad() .or. is_dipole() .or. is_matrix()) ) then
print *," "
print *,"code: ",code," el: ",el," THICK ELEMENT FOUND"
Expand Down Expand Up @@ -3145,7 +3149,8 @@ subroutine ttdpdg(track, ktrack)
end subroutine ttdpdg

subroutine trsol(track,ktrack)
use math_constfi, only : one, two, half

use math_constfi, only : zero, one, two, half, four
implicit none
!----------------------------------------------------------------------*
! Purpose: *
Expand All @@ -3160,66 +3165,97 @@ subroutine trsol(track,ktrack)
integer :: ktrack

integer :: i
double precision :: beta
double precision :: sk, skl, sks, sksl, cosTh, sinTh, Q, R, Z
double precision :: bet0, bet
double precision :: sk, skl, cosTh, sinTh, Q, R, Z
double precision :: xf, yf, pxf, pyf, sigf, psigf, bvk
double precision :: onedp, fpsig, fppsig

double precision :: get_value, node_value

double precision :: omega, length, length_
double precision :: x_, y_, z_, px_, py_, pt_

!---- Initialize.
! dtbyds = get_value('probe ','dtbyds ')
! gamma = get_value('probe ','gamma ')
beta = get_value('probe ','beta ')
! deltap = get_value('probe ','deltap ')
!
bet0 = get_value('probe ','beta ')

!---- Get solenoid parameters
! elrad = node_value('lrad ')
sksl = node_value('ksi ')
sks = node_value('ks ')

!---- BV flag
bvk = node_value('other_bv ')
sks = sks * bvk
sksl = sksl * bvk

!---- Set up strengths
! sk = sks / two / (one + deltap)
sk = sks / two
skl = sksl / two

!---- Loop over particles
do i = 1, ktrack
! Ripken formulae p.28 (3.35 and 3.36)
xf = track(1,i)
yf = track(3,i)
psigf = track(6,i) / beta

! We do not use a constant deltap!!!!! WE use full 6D formulae!
onedp = sqrt( one + 2*psigf + (beta**2)*(psigf**2) )
fpsig = onedp - one
fppsig = ( one + (beta**2)*psigf ) / onedp

! Set up C,S, Q,R,Z
cosTh = cos(skl/onedp)
sinTh = sin(skl/onedp)
Q = -skl * sk / onedp
R = fppsig / (onedp**2) * skl * sk
Z = fppsig / (onedp**2) * skl

pxf = track(2,i) + xf*Q
pyf = track(4,i) + yf*Q
sigf = track(5,i)*beta - half*(xf**2 + yf**2)*R

! Ripken formulae p.29 (3.37)
track(1,i) = xf * cosTh + yf * sinTh
track(2,i) = pxf * cosTh + pyf * sinTh
track(3,i) = -xf * sinTh + yf * cosTh
track(4,i) = -pxf * sinTh + pyf * cosTh
track(5,i) = (sigf + (xf*pyf - yf*pxf)*Z) / beta
! track(6,i) = psigf*beta

enddo
sk = bvk * node_value('ks ') / two
length = node_value('l ')

if (length.eq.zero) then

skl = bvk * node_value('ksi ') / two

!---- Loop over particles
do i = 1, ktrack
! Ripken formulae p.28 (3.35 and 3.36)
xf = track(1,i)
yf = track(3,i)
psigf = track(6,i) / bet0

! We do not use a constant deltap!!!!! WE use full 6D formulae!
onedp = sqrt( one + two*psigf + (bet0**2)*(psigf**2) )
fpsig = onedp - one
fppsig = ( one + (bet0**2)*psigf ) / onedp

! Set up C,S, Q,R,Z
cosTh = cos(skl/onedp)
sinTh = sin(skl/onedp)
Q = -skl * sk / onedp
R = fppsig / (onedp**2) * skl * sk
Z = fppsig / (onedp**2) * skl

pxf = track(2,i) + xf*Q
pyf = track(4,i) + yf*Q
sigf = track(5,i)*bet0 - half*(xf**2 + yf**2)*R

! Ripken formulae p.29 (3.37)
track(1,i) = xf * cosTh + yf * sinTh
track(2,i) = pxf * cosTh + pyf * sinTh
track(3,i) = -xf * sinTh + yf * cosTh
track(4,i) = -pxf * sinTh + pyf * cosTh
track(5,i) = (sigf + (xf*pyf - yf*pxf)*Z) / bet0
! track(6,i) = psigf*bet0
enddo
else
if (sk.ne.zero) then
skl = sk*length

!---- Loop over particles
do i = 1, ktrack
! initial phase space coordinates
x_ = track(1,i)
y_ = track(3,i)
px_ = track(2,i)
py_ = track(4,i)
z_ = track(5,i)
pt_ = track(6,i)

! set up constants
onedp = sqrt(one + two*pt_/bet0 + pt_**2);
bet = onedp / (one/bet0 + pt_);

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

! total path length traveled by the particle
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;

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) = (omega*((cosTh-one)*y_-sinTh*x_)+py_*sinTh+px_*(one+cosTh))/two;
track(4,i) = (omega*((one-cosTh)*x_-sinTh*y_)-px_*sinTh+py_*(one+cosTh))/two;
track(5,i) = z_ + length/bet0 - length_/bet;
enddo
else
call ttdrf(length,track,ktrack);
endif
endif
end subroutine trsol

subroutine tttrans(track,ktrack)
Expand Down Expand Up @@ -4689,6 +4725,11 @@ subroutine trphot(el,curv,rfac,deltap)
amean = five * sqrt(three) / (twelve * hbar * clight) * abs(arad * pc * (one+deltap) * el * curv)
ucrit = three/two * hbar * clight * gamma**3 * abs(curv)
sumxi = zero
if (amean.gt.3d-1) then
print *,"More than 0.3 photons emitted in element."
print *,"You might want to consider increasing the number of slices to reduce this number."
endif

if (amean .gt. zero) then
call dpoissn(amean, nphot, ierror)

Expand Down
3 changes: 3 additions & 0 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0001.cfg
View file
Open in desktop
@@ -0,0 +1,3 @@
4-6 * skip # head

* * any abs=1e-10 rel=5e-7
10 changes: 10 additions & 0 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0001.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0001"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %17s "5.03.07 Darwin 64"
@ DATE %08s "06/11/17"
@ TIME %08s "18.38.10"
* 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
3 changes: 3 additions & 0 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0002.cfg
View file
Open in desktop
@@ -0,0 +1,3 @@
4-6 * skip # head

* * any abs=1e-10 rel=5e-7 # SLC6, GCC32
10 changes: 10 additions & 0 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0002.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0002"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %17s "5.03.07 Darwin 64"
@ DATE %08s "06/11/17"
@ TIME %08s "18.38.10"
* 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
3 changes: 3 additions & 0 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0003.cfg
View file
Open in desktop
@@ -0,0 +1,3 @@
4-6 * skip # head

* * any abs=1e-10 rel=5e-7 # SLC6, GCC32
10 changes: 10 additions & 0 deletions tests/test-thick-solenoid/madx-thick.obs0001.p0003.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0003"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %17s "5.03.07 Darwin 64"
@ DATE %08s "06/11/17"
@ TIME %08s "18.38.10"
* 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 changes: 3 additions & 0 deletions tests/test-thick-solenoid/madx-thin.obs0001.p0001.cfg
View file
Open in desktop
@@ -0,0 +1,3 @@
4-6 * skip # head

* * any abs=1e-10 rel=5e-7
10 changes: 10 additions & 0 deletions tests/test-thick-solenoid/madx-thin.obs0001.p0001.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0001"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %17s "5.03.07 Darwin 64"
@ DATE %08s "06/11/17"
@ TIME %08s "18.38.10"
* 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.0004179968082 -9.084288303e-05 -2.711952481e-19 -3.599890026e-20 -1.054344426e-07 -0.001 0 1
3 changes: 3 additions & 0 deletions tests/test-thick-solenoid/madx-thin.obs0001.p0002.cfg
View file
Open in desktop
@@ -0,0 +1,3 @@
4-6 * skip # head

* * any abs=1e-10 rel=5e-7 # SLC6, GCC32
10 changes: 10 additions & 0 deletions tests/test-thick-solenoid/madx-thin.obs0001.p0002.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0002"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %17s "5.03.07 Darwin 64"
@ DATE %08s "06/11/17"
@ TIME %08s "18.38.10"
* 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.0004161771647 -9.092639127e-05 -2.050193423e-19 5.145725155e-20 -1.000038811e-07 2.220446049e-16 0 1
3 changes: 3 additions & 0 deletions tests/test-thick-solenoid/madx-thin.obs0001.p0003.cfg
View file
Open in desktop
@@ -0,0 +1,3 @@
4-6 * skip # head

* * any abs=1e-10 rel=5e-7 # SLC6, GCC32
10 changes: 10 additions & 0 deletions tests/test-thick-solenoid/madx-thin.obs0001.p0003.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0003"
@ TYPE %08s "TRACKOBS"
@ TITLE %28s "thick-solenoid tracking test"
@ ORIGIN %17s "5.03.07 Darwin 64"
@ DATE %08s "06/11/17"
@ TIME %08s "18.38.10"
* 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.0004143594937 -9.100936927e-05 -4.805886845e-19 -2.922263668e-20 -9.458958447e-08 0.001 0 1
2 changes: 2 additions & 0 deletions tests/test-thick-solenoid/test-thick-solenoid.cfg
View file
Open in desktop
@@ -0,0 +1,2 @@
1-7 * skip # head
* * any abs=1e-10 rel=5e-7 # global
41 changes: 41 additions & 0 deletions tests/test-thick-solenoid/test-thick-solenoid.madx
View file
Open in desktop
@@ -0,0 +1,41 @@
option, -debug, -echo;

title, "thick-solenoid tracking test";

L_sol := 10.;
Ks_sol := 0.2;

sol1: solenoid, l=L_sol, KS=Ks_sol, KSI=Ks_sol*L_sol;
sol2: solenoid, l=L_sol, KS=-Ks_sol, KSI=Ks_sol*L_sol;

myseq: sequence, l=2*L_sol;
sol1, at:= L_sol/2;
sol2, at:= L_sol+L_sol/2;
endsequence;

beam,energy=1; !
use, sequence=myseq;

track, onepass, dump, file= "madx-thick";
start, x=1e-3, y=0, px=0, py=0, t=0, pt=-1e-3;
start, x=1e-3, y=0, px=0, py=0, t=0, pt=0.0;
start, x=1e-3, y=0, px=0, py=0, t=0, pt=1e-3;
run, turns=1;
endtrack;

use, sequence=myseq;
select, flag=makethin, class=solenoid, thick=false, slice=50;
makethin, sequence=myseq, style=teapot;

!save, sequence= myseq,file= "myseq_thin.seq";
!call, file= "myseq_thin.seq";

use, sequence=myseq;
track, onepass, dump, file= "madx-thin";
start, x=1e-3, y=0, px=0, py=0, t=0, pt=-1e-3;
start, x=1e-3, y=0, px=0, py=0, t=0, pt=0.0;
start, x=1e-3, y=0, px=0, py=0, t=0, pt=1e-3;
run, turns=1;
endtrack;
stop;