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

Added support for thick RBEND in TRACK #697

Merged
merged 5 commits into from Dec 18, 2018
Merged

Added support for thick RBEND in TRACK #697

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
59d5e98
Added support for thick RBEND
alatina Nov 15, 2018
5b46232
Update trrun.f90
alatina Nov 16, 2018
048d83f
Merge branch 'master' into al-dev
alatina Nov 27, 2018
853ffb1
Addressing issue #704
alatina Nov 29, 2018
2749610
Added example for thick tracking in an RBEND
alatina Dec 17, 2018
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
2 changes: 1 addition & 1 deletion Makefile_test
View file
Open in desktop
Expand Up @@ -61,7 +61,7 @@ test-dynap \
test-c6t test-c6t-2 test-c6t-3 test-c6t-4 test-c6t-5 \
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-thick-rbend 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 test-match-10\
Expand Down
74 changes: 41 additions & 33 deletions src/trrun.f90
View file
Open in desktop
Expand Up @@ -471,10 +471,11 @@ subroutine trrun(switch, turns, orbit0, rt, part_id, last_turn, last_pos, &

if ( code.ne.code_drift .and. &
code.ne.code_quadrupole .and. &
code.ne.code_rbend .and. &
code.ne.code_sbend .and. &
code.ne.code_matrix .and. &
code.ne.code_solenoid .and. &
el.ne.zero ) then ! rbend missing ?
el.ne.zero ) then
!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 @@ -1063,16 +1064,14 @@ subroutine ttmult(track,ktrack,dxt,dyt,turn)
if (damp) then
do jtrk = 1,ktrack
curv = sqrt((dipr + dxt(jtrk))**2 + (dipi + dyt(jtrk))**2) / elrad

px = track(2,jtrk)
py = track(4,jtrk)
pt = track(6,jtrk)
if (quantum) then
call trphot(elrad,curv,rfac,deltas)
call trphot(elrad,curv,rfac,pt)
else
rfac = const * curv**2 * elrad
endif

px = track(2,jtrk)
py = track(4,jtrk)
pt = track(6,jtrk)
track(2,jtrk) = px - rfac * (one + pt) * px
track(4,jtrk) = py - rfac * (one + pt) * py
track(6,jtrk) = pt - rfac * (one + pt) ** 2
Expand Down Expand Up @@ -1115,16 +1114,14 @@ subroutine ttmult(track,ktrack,dxt,dyt,turn)
if (damp) then
do jtrk = 1,ktrack
curv = sqrt((dipr + dxt(jtrk))**2 + (dipi + dyt(jtrk))**2) / elrad

px = track(2,jtrk)
py = track(4,jtrk)
pt = track(6,jtrk)
if (quantum) then
call trphot(elrad,curv,rfac,deltas)
call trphot(elrad,curv,rfac,pt)
else
rfac = const * curv**2 * elrad
endif

px = track(2,jtrk)
py = track(4,jtrk)
pt = track(6,jtrk)
track(2,jtrk) = px - rfac * (one + pt) * px
track(4,jtrk) = py - rfac * (one + pt) * py
track(6,jtrk) = pt - rfac * (one + pt) ** 2
Expand Down Expand Up @@ -1786,10 +1783,10 @@ subroutine ttcorr(el,track,ktrack,turn)

if (damp) then !---- Full damping.
do i = 1, ktrack
if (quantum) call trphot(el,curv,rfac,deltas)
px = track(2,i)
py = track(4,i)
pt = track(6,i)
if (quantum) call trphot(el,curv,rfac,pt)
track(2,i) = px - rfac * (one + pt) * px
track(4,i) = py - rfac * (one + pt) * py
track(6,i) = pt - rfac * (one + pt) ** 2
Expand Down Expand Up @@ -1870,10 +1867,10 @@ subroutine ttcorr(el,track,ktrack,turn)
if (radiate .and. el .ne. 0) then
if (damp) then !---- Full damping.
do i = 1, ktrack
if (quantum) call trphot(el,curv,rfac,deltas)
px = track(2,i)
py = track(4,i)
pt = track(6,i)
if (quantum) call trphot(el,curv,rfac,pt)
track(2,i) = px - rfac * (one + pt) * px
track(4,i) = py - rfac * (one + pt) * py
track(6,i) = pt - rfac * (one + pt) ** 2
Expand All @@ -1882,7 +1879,6 @@ subroutine ttcorr(el,track,ktrack,turn)
rpx = rfac * (one + track(6,1)) * track(2,1)
rpy = rfac * (one + track(6,1)) * track(4,1)
rpt = rfac * (one + track(6,1)) ** 2

do i = 1, ktrack
track(2,i) = track(2,i) - rpx
track(4,i) = track(4,i) - rpy
Expand Down Expand Up @@ -3444,8 +3440,11 @@ subroutine trclor(switch,orbit0)
if (code .eq. code_placeholder) code = code_instrument
el = node_value('l ')
if (itra .eq. 1 .and. &
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 ?
code.ne.code_drift .and. &
code.ne.code_quadrupole .and. &
code.ne.code_rbend .and. &
code.ne.code_sbend .and. &
code.ne.code_matrix .and. el.ne.zero ) then
! .not.(is_drift() .or. is_thin() .or. is_quad() .or. is_dipole() .or. is_matrix()) ) then
print *,"\ncode: ",code," el: ",el," THICK ELEMENT FOUND\n"
print *,"Track dies nicely"
Expand Down Expand Up @@ -4159,7 +4158,7 @@ subroutine ttrfmult(track, ktrack, turn)
if (radiate .and. elrad .ne. zero) then
if (quantum) then
curv = sqrt(dpx**2+dpy**2) / elrad;
call trphot(elrad,curv,rfac,deltas)
call trphot(elrad,curv,rfac,pt)
else
rfac = arad * gammas**3 * (dpx**2+dpy**2) / (three*elrad)
endif
Expand Down Expand Up @@ -4188,7 +4187,7 @@ subroutine ttrfmult(track, ktrack, turn)
if (radiate .and. elrad .ne. zero) then
if (quantum) then
curv = sqrt(dpx**2+dpy**2) / elrad;
call trphot(elrad,curv,rfac,deltas)
call trphot(elrad,curv,rfac,pt)
else
rfac = arad * gammas**3 * (dpx**2+dpy**2) / (three*elrad)
endif
Expand Down Expand Up @@ -4415,7 +4414,7 @@ subroutine tttquad(track, ktrack)
hy = ( k1*y) / delta_plus_1;
if (quantum) then
curv = sqrt(hx**2+hy**2);
call trphot(length,curv,rfac,deltas)
call trphot(length,curv,rfac,pt)
else
rfac = (arad * gamma**3 * length / three) * (hx**2 + hy**2);
endif
Expand Down Expand Up @@ -4443,7 +4442,7 @@ subroutine tttquad(track, ktrack)
hy = ( k1*y) / delta_plus_1;
if (quantum) then
curv = sqrt(hx**2+hy**2);
call trphot(length,curv,rfac,deltas)
call trphot(length,curv,rfac,pt)
else
rfac = (arad * gamma**3 * length / three) * (hx**2 + hy**2);
endif
Expand Down Expand Up @@ -4489,6 +4488,8 @@ subroutine tttdipole(track, ktrack)
use twtrrfi
use trackfi
use math_constfi, only : zero, one, two, three, half
use code_constfi

implicit none
!-------------------------*
! Andrea Latina 2013-2014 *
Expand All @@ -4515,8 +4516,9 @@ subroutine tttdipole(track, ktrack)
double precision, external :: node_value, get_value
double precision :: f_errors(0:maxferr)
integer, external :: node_fd_errors
integer :: n_ferr
integer :: n_ferr, code

code = node_value('mad8_type ')
arad = get_value('probe ','arad ')
gamma = get_value('probe ','gamma ')
radiate = get_value('probe ','radiate ') .ne. zero
Expand All @@ -4536,6 +4538,11 @@ subroutine tttdipole(track, ktrack)
k0 = h;
k1 = node_value('k1 ');

if (code .eq. code_rbend) then
e1 = e1 + angle / two;
e2 = e2 + angle / two;
endif

!---- Apply errors
f_errors = zero
n_ferr = node_fd_errors(f_errors)
Expand Down Expand Up @@ -4567,7 +4574,7 @@ subroutine tttdipole(track, ktrack)
hy = ( k1*y) / delta_plus_1;
if (quantum) then
curv = sqrt(hx**2+hy**2);
call trphot(length * (one + h*x) - two * tan(e1)*x, curv, rfac, deltas);
call trphot(length * (one + h*x) - two * tan(e1)*x, curv, rfac, pt);
else
rfac = (arad * gamma**3 * two / three) * (hx**2 + hy**2) * (length / two * (one + h*x) - tan(e1)*x)
endif
Expand Down Expand Up @@ -4597,7 +4604,7 @@ subroutine tttdipole(track, ktrack)
hy = ( k1*y) / delta_plus_1;
if (quantum) then
curv = sqrt(hx**2+hy**2);
call trphot(length * (one + h*x) - two * tan(e2)*x, curv, rfac, deltas);
call trphot(length * (one + h*x) - two * tan(e2)*x, curv, rfac, pt);
else
rfac = (arad * gamma**3 * two / three) * (hx**2 + hy**2) * (length / two * (one + h*x) - tan(e2)*x)
endif
Expand Down Expand Up @@ -4635,9 +4642,10 @@ subroutine tttdipole(track, ktrack)

end subroutine tttdipole

subroutine trphot(el,curv,rfac,deltap)
subroutine trphot(el,curv,rfac,pt)
use math_constfi, only : zero, one, two, three, five, twelve
use phys_constfi, only : clight, hbar
use trackfi, only : arad, gammas, betas
implicit none
!----------------------------------------------------------------------*
! Purpose: *
Expand All @@ -4663,17 +4671,17 @@ subroutine trphot(el,curv,rfac,deltap)
! The random integers are generated in the range [0, MAXRAN). *
!---- Table definition: maxtab, taby(maxtab), tabxi(maxtab) *
!----------------------------------------------------------------------*
double precision :: el, curv, rfac, deltap
double precision :: el, curv, rfac, pt

integer :: i, ierror, j, nphot
double precision :: amean, dlogr, slope, ucrit, xi, sumxi, InvSynFracInt,rn_arg
integer, parameter :: maxtab=101
double precision :: tabxi(maxtab),taby(maxtab)
double precision :: arad, pc, gamma, amass
double precision :: pc, gamma, amass
character(len=20) text
integer, external :: get_option
integer :: synrad
double precision :: get_value, frndm
double precision :: get_value, frndm, delta_plus_1
double precision, parameter :: fac1=3.256223d0
! integer, parameter :: maxran=1000000000

Expand Down Expand Up @@ -4733,16 +4741,16 @@ subroutine trphot(el,curv,rfac,deltap)
2.64617491d0 /

!Get constants
arad = get_value('probe ','arad ')
pc = get_value('probe ','pc ')
amass = get_value('probe ','mass ')
gamma = get_value('probe ','gamma ')
deltap = get_value('probe ','deltap ')

gamma = (betas*pt + one)*gammas;
delta_plus_1 = sqrt(pt*pt + two*pt/betas + one);

!---- AMEAN is the average number of photons emitted.,
! NPHOT is the integer number generated from Poisson's law.
!-AL- AMEAN implicitly takes el / 2 (half the element length)
amean = five * sqrt(three) / (twelve * hbar * clight) * abs(arad * pc * (one+deltap) * el * curv)
amean = five * sqrt(three) / (twelve * hbar * clight) * abs(arad * pc * delta_plus_1 * el * curv)
ucrit = three/two * hbar * clight * gamma**3 * abs(curv)
sumxi = zero
if (amean.gt.3d-1) then
Expand Down
2 changes: 2 additions & 0 deletions tests/test-thick-rbend/test-thick-rbend.cfg
View file
Open in desktop
@@ -0,0 +1,2 @@
1-7 * skip # head
* * any abs=1e-10 rel=5e-7 # global
23 changes: 23 additions & 0 deletions tests/test-thick-rbend/test-thick-rbend.madx
View file
Open in desktop
@@ -0,0 +1,23 @@
option, -debug, -echo;

title, "thick-rbend tracking test";

lmb := 1.;
amb := 0.1;

mb1: rbend, l:=lmb, angle:=amb;

myseq: sequence, l=1.1, refer=centre;
mb1, at= 1.1/2;
endsequence;

beam,energy=1000;
use, sequence=myseq;

track, onepass, dump;
start, x=0, px=0, y=1e-6, py=0, t=0, pt=-1e-2;
start, x=0, px=0, y=1e-6, py=0, t=0, pt=0;
start, x=0, px=0, y=1e-6, py=0, t=0, pt=1e-2;
run, turns=1;
endtrack;
stop;
45 changes: 45 additions & 0 deletions tests/test-thick-rbend/test-thick-rbend.ref
View file
Open in desktop
@@ -0,0 +1,45 @@

++++++++++++++++++++++++++++++++++++++++++++
+ MAD-X 5.04.02 (64 bit, Darwin) +
+ Support: mad@cern.ch, http://cern.ch/mad +
+ Release date: 2018.10.03 +
+ Execution date: 2018.12.17 23:05:33 +
++++++++++++++++++++++++++++++++++++++++++++
option, -debug, -echo;

enter TRACK module
one pass is on

++++++ table: tracksumm

number turn x px
1 0 0 0
2 0 0 0
3 0 0 0
1 1 -0.0005550457922 -0.0009983173484
2 1 0 0
3 1 0.0005440646768 0.0009983506517

y py t pt
1e-06 0 0 -0.01
1e-06 0 0 0
1e-06 0 0 0.01
1e-06 0 1.663843534e-05 -0.01
1e-06 0 0 0
1e-06 0 -1.668780949e-05 0.01

s e
0 0
0 0
0 0
1.1 0
1.1 0
1.1 0
exit TRACK module


Number of warnings: 0

++++++++++++++++++++++++++++++++++++++++++++
+ MAD-X finished normally +
++++++++++++++++++++++++++++++++++++++++++++
4 changes: 4 additions & 0 deletions tests/test-thick-rbend/track.obs0001.p0001.cfg
View file
Open in desktop
@@ -0,0 +1,4 @@
4-6 * skip # head

* * any abs=1e-10 rel=5e-7
46 4 any abs=1e-10 rel=1.2e-6 # SLC6, GCC32, 2013.11.08
10 changes: 10 additions & 0 deletions tests/test-thick-rbend/track.obs0001.p0001.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0001"
@ TYPE %08s "TRACKOBS"
@ TITLE %25s "thick-rbend tracking test"
@ ORIGIN %17s "5.04.02 Darwin 64"
@ DATE %08s "17/12/18"
@ TIME %08s "23.05.33"
* NUMBER TURN X PX Y PY T PT S E
$ %d %d %le %le %le %le %le %le %le %le
1 0 0 0 1e-06 0 0 -0.01 0 1000
1 1 -0.0005550457922 -0.0009983173484 1e-06 0 1.663843534e-05 -0.01 0 1000
3 changes: 3 additions & 0 deletions tests/test-thick-rbend/track.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-rbend/track.obs0001.p0002.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0002"
@ TYPE %08s "TRACKOBS"
@ TITLE %25s "thick-rbend tracking test"
@ ORIGIN %17s "5.04.02 Darwin 64"
@ DATE %08s "17/12/18"
@ TIME %08s "23.05.33"
* NUMBER TURN X PX Y PY T PT S E
$ %d %d %le %le %le %le %le %le %le %le
2 0 0 0 1e-06 0 0 0 0 1000
2 1 0 0 1e-06 0 0 0 0 1000
3 changes: 3 additions & 0 deletions tests/test-thick-rbend/track.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-rbend/track.obs0001.p0003.ref
View file
Open in desktop
@@ -0,0 +1,10 @@
@ NAME %19s "TRACK.OBS0001.P0003"
@ TYPE %08s "TRACKOBS"
@ TITLE %25s "thick-rbend tracking test"
@ ORIGIN %17s "5.04.02 Darwin 64"
@ DATE %08s "17/12/18"
@ TIME %08s "23.05.33"
* NUMBER TURN X PX Y PY T PT S E
$ %d %d %le %le %le %le %le %le %le %le
3 0 0 0 1e-06 0 0 0.01 0 1000
3 1 0.0005440646768 0.0009983506517 1e-06 0 -1.668780949e-05 0.01 0 1000