From 1c7780d03cd72745c7646dd07daaed98e15a8769 Mon Sep 17 00:00:00 2001 From: Andrea Latina Date: Mon, 6 Nov 2017 15:37:53 +0100 Subject: [PATCH 1/6] Implemented thick solenoid tracking --- src/trrun.f90 | 123 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 86 insertions(+), 37 deletions(-) diff --git a/src/trrun.f90 b/src/trrun.f90 index 2663942ec..ace570a9d 100644 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -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" @@ -3081,7 +3085,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 implicit none !----------------------------------------------------------------------* ! Purpose: * @@ -3096,66 +3101,105 @@ subroutine trsol(track,ktrack) integer :: ktrack integer :: i - double precision :: beta + double precision :: bet0, bet double precision :: sk, skl, sks, sksl, 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 ') + bet0 = get_value('probe ','beta ') ! deltap = get_value('probe ','deltap ') ! !---- Get solenoid parameters ! elrad = node_value('lrad ') sksl = node_value('ksi ') sks = node_value('ks ') + length = node_value('l '); !---- 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 + + if (length.eq.zero) then + + !---- 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 + 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(skl/onedp) + sinTh = sin(skl/onedp) + omega = sk/onedp; + + ! total path length traveled by the particle + length_ = length -half/(onedp**2)*(& + omega**2*(cosTh*sinTh-length*omega)*(x_**2+y_**2)+& + omega*(one-(cosTh**2-sinTh**2))*(px_*x_+py_*y_)+& + (-cosTh*sinTh-length*omega)*(px_**2+py_**2))/(two*omega); + + track(1,i) = cosTh*sinTh*y_+cosTh**2*x_+(py_*sinTh**2)/omega+(px_*cosTh*sinTh)/omega; + track(2,i) = -omega*sinTh**2*y_-omega*cosTh*sinTh*x_+cosTh*sinTh*py_+cosTh**2*px_; + track(3,i) = cosTh**2*y_-cosTh*sinTh*x_-(px_*sinTh**2)/omega+(py_*cosTh*sinTh)/omega; + track(4,i) = -omega*cosTh*sinTh*y_+omega*sinTh**2*x_+cosTh**2*py_-cosTh*sinTh*px_; + track(5,i) = z_ + length/bet0 - length_/bet; - !---- 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 - enddo + endif end subroutine trsol subroutine tttrans(track,ktrack) @@ -4625,6 +4669,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) From 72bf11111f24234a963fa976ddcdbdaa17f04500 Mon Sep 17 00:00:00 2001 From: Andrea Latina Date: Mon, 6 Nov 2017 15:41:42 +0100 Subject: [PATCH 2/6] Cosmetics... --- src/trrun.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/trrun.f90 b/src/trrun.f90 index ace570a9d..19e30ec7e 100644 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -3134,7 +3134,6 @@ subroutine trsol(track,ktrack) skl = sksl / two if (length.eq.zero) then - !---- Loop over particles do i = 1, ktrack ! Ripken formulae p.28 (3.35 and 3.36) @@ -3166,7 +3165,8 @@ subroutine trsol(track,ktrack) track(5,i) = (sigf + (xf*pyf - yf*pxf)*Z) / bet0 ! track(6,i) = psigf*bet0 enddo - else + else + !---- Loop over particles do i = 1, ktrack ! initial phase space coordinates x_ = track(1,i); From f9d9349d5c937fc24c161ece17286cf0b167477b Mon Sep 17 00:00:00 2001 From: Andrea Latina Date: Mon, 6 Nov 2017 18:56:37 +0100 Subject: [PATCH 3/6] Improved thick solenoid map and added test --- Makefile_test | 1 + src/trrun.f90 | 58 +++++++------- .../madx-thick.obs0001.p0001.cfg | 3 + .../madx-thick.obs0001.p0001.ref | 10 +++ .../madx-thick.obs0001.p0002.cfg | 3 + .../madx-thick.obs0001.p0002.ref | 10 +++ .../madx-thick.obs0001.p0003.cfg | 3 + .../madx-thick.obs0001.p0003.ref | 10 +++ .../madx-thin.obs0001.p0001.cfg | 3 + .../madx-thin.obs0001.p0001.ref | 10 +++ .../madx-thin.obs0001.p0002.cfg | 3 + .../madx-thin.obs0001.p0002.ref | 10 +++ .../madx-thin.obs0001.p0003.cfg | 3 + .../madx-thin.obs0001.p0003.ref | 10 +++ .../test-thick-solenoid.cfg | 2 + .../test-thick-solenoid.madx | 38 +++++++++ .../test-thick-solenoid.ref | 77 +++++++++++++++++++ 17 files changed, 222 insertions(+), 32 deletions(-) create mode 100644 tests/test-thick-solenoid/madx-thick.obs0001.p0001.cfg create mode 100644 tests/test-thick-solenoid/madx-thick.obs0001.p0001.ref create mode 100644 tests/test-thick-solenoid/madx-thick.obs0001.p0002.cfg create mode 100644 tests/test-thick-solenoid/madx-thick.obs0001.p0002.ref create mode 100644 tests/test-thick-solenoid/madx-thick.obs0001.p0003.cfg create mode 100644 tests/test-thick-solenoid/madx-thick.obs0001.p0003.ref create mode 100644 tests/test-thick-solenoid/madx-thin.obs0001.p0001.cfg create mode 100644 tests/test-thick-solenoid/madx-thin.obs0001.p0001.ref create mode 100644 tests/test-thick-solenoid/madx-thin.obs0001.p0002.cfg create mode 100644 tests/test-thick-solenoid/madx-thin.obs0001.p0002.ref create mode 100644 tests/test-thick-solenoid/madx-thin.obs0001.p0003.cfg create mode 100644 tests/test-thick-solenoid/madx-thin.obs0001.p0003.ref create mode 100644 tests/test-thick-solenoid/test-thick-solenoid.cfg create mode 100644 tests/test-thick-solenoid/test-thick-solenoid.madx create mode 100644 tests/test-thick-solenoid/test-thick-solenoid.ref diff --git a/Makefile_test b/Makefile_test index f53c0d775..857877b70 100644 --- a/Makefile_test +++ b/Makefile_test @@ -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\ diff --git a/src/trrun.f90 b/src/trrun.f90 index 1e5a651dc..6d5a94384 100644 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -3150,7 +3150,7 @@ end subroutine ttdpdg subroutine trsol(track,ktrack) - use math_constfi, only : zero, one, two, half + use math_constfi, only : zero, one, two, half, four implicit none !----------------------------------------------------------------------* ! Purpose: * @@ -3166,7 +3166,7 @@ subroutine trsol(track,ktrack) integer :: i double precision :: bet0, bet - double precision :: sk, skl, sks, sksl, cosTh, sinTh, Q, R, Z + double precision :: sk, skl, cosTh, sinTh, Q, R, Z double precision :: xf, yf, pxf, pyf, sigf, psigf, bvk double precision :: onedp, fpsig, fppsig @@ -3183,21 +3183,14 @@ subroutine trsol(track,ktrack) ! !---- Get solenoid parameters ! elrad = node_value('lrad ') - sksl = node_value('ksi ') - sks = node_value('ks ') - length = node_value('l '); - - !---- 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 + 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) @@ -3229,40 +3222,41 @@ subroutine trsol(track,ktrack) track(5,i) = (sigf + (xf*pyf - yf*pxf)*Z) / bet0 ! track(6,i) = psigf*bet0 enddo + else + + 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); + 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(skl/onedp) - sinTh = sin(skl/onedp) + 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**2*(cosTh*sinTh-length*omega)*(x_**2+y_**2)+& - omega*(one-(cosTh**2-sinTh**2))*(px_*x_+py_*y_)+& - (-cosTh*sinTh-length*omega)*(px_**2+py_**2))/(two*omega); - - track(1,i) = cosTh*sinTh*y_+cosTh**2*x_+(py_*sinTh**2)/omega+(px_*cosTh*sinTh)/omega; - track(2,i) = -omega*sinTh**2*y_-omega*cosTh*sinTh*x_+cosTh*sinTh*py_+cosTh**2*px_; - track(3,i) = cosTh**2*y_-cosTh*sinTh*x_-(px_*sinTh**2)/omega+(py_*cosTh*sinTh)/omega; - track(4,i) = -omega*cosTh*sinTh*y_+omega*sinTh**2*x_+cosTh**2*py_-cosTh*sinTh*px_; + 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 - endif end subroutine trsol diff --git a/tests/test-thick-solenoid/madx-thick.obs0001.p0001.cfg b/tests/test-thick-solenoid/madx-thick.obs0001.p0001.cfg new file mode 100644 index 000000000..a34c9f0de --- /dev/null +++ b/tests/test-thick-solenoid/madx-thick.obs0001.p0001.cfg @@ -0,0 +1,3 @@ +4-6 * skip # head + +* * any abs=1e-10 rel=5e-7 diff --git a/tests/test-thick-solenoid/madx-thick.obs0001.p0001.ref b/tests/test-thick-solenoid/madx-thick.obs0001.p0001.ref new file mode 100644 index 000000000..3ef71d4f3 --- /dev/null +++ b/tests/test-thick-solenoid/madx-thick.obs0001.p0001.ref @@ -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 diff --git a/tests/test-thick-solenoid/madx-thick.obs0001.p0002.cfg b/tests/test-thick-solenoid/madx-thick.obs0001.p0002.cfg new file mode 100644 index 000000000..73364783d --- /dev/null +++ b/tests/test-thick-solenoid/madx-thick.obs0001.p0002.cfg @@ -0,0 +1,3 @@ +4-6 * skip # head + +* * any abs=1e-10 rel=5e-7 # SLC6, GCC32 diff --git a/tests/test-thick-solenoid/madx-thick.obs0001.p0002.ref b/tests/test-thick-solenoid/madx-thick.obs0001.p0002.ref new file mode 100644 index 000000000..bb9f235aa --- /dev/null +++ b/tests/test-thick-solenoid/madx-thick.obs0001.p0002.ref @@ -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 diff --git a/tests/test-thick-solenoid/madx-thick.obs0001.p0003.cfg b/tests/test-thick-solenoid/madx-thick.obs0001.p0003.cfg new file mode 100644 index 000000000..73364783d --- /dev/null +++ b/tests/test-thick-solenoid/madx-thick.obs0001.p0003.cfg @@ -0,0 +1,3 @@ +4-6 * skip # head + +* * any abs=1e-10 rel=5e-7 # SLC6, GCC32 diff --git a/tests/test-thick-solenoid/madx-thick.obs0001.p0003.ref b/tests/test-thick-solenoid/madx-thick.obs0001.p0003.ref new file mode 100644 index 000000000..02219aa83 --- /dev/null +++ b/tests/test-thick-solenoid/madx-thick.obs0001.p0003.ref @@ -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 diff --git a/tests/test-thick-solenoid/madx-thin.obs0001.p0001.cfg b/tests/test-thick-solenoid/madx-thin.obs0001.p0001.cfg new file mode 100644 index 000000000..a34c9f0de --- /dev/null +++ b/tests/test-thick-solenoid/madx-thin.obs0001.p0001.cfg @@ -0,0 +1,3 @@ +4-6 * skip # head + +* * any abs=1e-10 rel=5e-7 diff --git a/tests/test-thick-solenoid/madx-thin.obs0001.p0001.ref b/tests/test-thick-solenoid/madx-thin.obs0001.p0001.ref new file mode 100644 index 000000000..db9bac581 --- /dev/null +++ b/tests/test-thick-solenoid/madx-thin.obs0001.p0001.ref @@ -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 diff --git a/tests/test-thick-solenoid/madx-thin.obs0001.p0002.cfg b/tests/test-thick-solenoid/madx-thin.obs0001.p0002.cfg new file mode 100644 index 000000000..73364783d --- /dev/null +++ b/tests/test-thick-solenoid/madx-thin.obs0001.p0002.cfg @@ -0,0 +1,3 @@ +4-6 * skip # head + +* * any abs=1e-10 rel=5e-7 # SLC6, GCC32 diff --git a/tests/test-thick-solenoid/madx-thin.obs0001.p0002.ref b/tests/test-thick-solenoid/madx-thin.obs0001.p0002.ref new file mode 100644 index 000000000..005978748 --- /dev/null +++ b/tests/test-thick-solenoid/madx-thin.obs0001.p0002.ref @@ -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 diff --git a/tests/test-thick-solenoid/madx-thin.obs0001.p0003.cfg b/tests/test-thick-solenoid/madx-thin.obs0001.p0003.cfg new file mode 100644 index 000000000..73364783d --- /dev/null +++ b/tests/test-thick-solenoid/madx-thin.obs0001.p0003.cfg @@ -0,0 +1,3 @@ +4-6 * skip # head + +* * any abs=1e-10 rel=5e-7 # SLC6, GCC32 diff --git a/tests/test-thick-solenoid/madx-thin.obs0001.p0003.ref b/tests/test-thick-solenoid/madx-thin.obs0001.p0003.ref new file mode 100644 index 000000000..bc4f106cc --- /dev/null +++ b/tests/test-thick-solenoid/madx-thin.obs0001.p0003.ref @@ -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 diff --git a/tests/test-thick-solenoid/test-thick-solenoid.cfg b/tests/test-thick-solenoid/test-thick-solenoid.cfg new file mode 100644 index 000000000..2ae61b9be --- /dev/null +++ b/tests/test-thick-solenoid/test-thick-solenoid.cfg @@ -0,0 +1,2 @@ +1-7 * skip # head +* * any abs=1e-10 rel=5e-7 # global diff --git a/tests/test-thick-solenoid/test-thick-solenoid.madx b/tests/test-thick-solenoid/test-thick-solenoid.madx new file mode 100644 index 000000000..2d7909e1d --- /dev/null +++ b/tests/test-thick-solenoid/test-thick-solenoid.madx @@ -0,0 +1,38 @@ +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; + +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; + diff --git a/tests/test-thick-solenoid/test-thick-solenoid.ref b/tests/test-thick-solenoid/test-thick-solenoid.ref new file mode 100644 index 000000000..7c71aa9db --- /dev/null +++ b/tests/test-thick-solenoid/test-thick-solenoid.ref @@ -0,0 +1,77 @@ + + ++++++++++++++++++++++++++++++++++++++++++++ + + MAD-X 5.03.07 (64 bit, Darwin) + + + Support: mad@cern.ch, http://cern.ch/mad + + + Release date: 2017.10.20 + + + Execution date: 2017.11.06 18:38:10 + + ++++++++++++++++++++++++++++++++++++++++++++ +option, -debug, -echo; + +enter TRACK module +one pass is on + +++++++ table: tracksumm + + number turn x px + 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 + + y py t pt + 0 0 0 -0.001 + 0 0 0 2.220446049e-16 + 0 0 0 0.001 + 0 0 -6.495195848e-08 -0.001 + 0 3.388131789e-21 -5.946003867e-08 2.220446049e-16 + 6.783039842e-20 0 -5.398501379e-08 0.001 + + s e + 0 0 + 0 0 + 0 0 + 20 0 + 20 0 + 20 0 +exit TRACK module + +makethin: style chosen : teapot +makethin: slicing sequence : myseq +enter TRACK module +one pass is on + +++++++ table: tracksumm + + number turn x px + 1 0 0.001 0 + 2 0 0.001 0 + 3 0 0.001 0 + 1 1 -0.0004179968082 -9.084288303e-05 + 2 1 -0.0004161771647 -9.092639127e-05 + 3 1 -0.0004143594937 -9.100936927e-05 + + y py t pt + 0 0 0 -0.001 + 0 0 0 2.220446049e-16 + 0 0 0 0.001 + -2.711952481e-19 -3.599890026e-20 -1.054344426e-07 -0.001 + -2.050193423e-19 5.145725155e-20 -1.000038811e-07 2.220446049e-16 + -4.805886845e-19 -2.922263668e-20 -9.458958447e-08 0.001 + + s e + 0 0 + 0 0 + 0 0 + 20 0 + 20 0 + 20 0 +exit TRACK module + + + Number of warnings: 0 + + ++++++++++++++++++++++++++++++++++++++++++++ + + MAD-X finished normally + + ++++++++++++++++++++++++++++++++++++++++++++ From 270fd79fee5baf4cec5dfdcc7028d00c0b3c139a Mon Sep 17 00:00:00 2001 From: Andrea Date: Mon, 6 Nov 2017 20:43:37 +0100 Subject: [PATCH 4/6] Implemented drift transfer map when thick solenoid has strength zero --- src/trrun.f90 | 35 +++++++++++++++++++++++------------ 1 file changed, 23 insertions(+), 12 deletions(-) diff --git a/src/trrun.f90 b/src/trrun.f90 index 6d5a94384..d514b8820 100644 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -3241,21 +3241,32 @@ subroutine trsol(track,ktrack) 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; + if (skl.ne.zero) then + + ! 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; - ! 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; - 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; + else + + length_ = length/sqrt(onedp**2-px_**2-py_**2); ! length/pz + track(1,i) = x_ + px_*length_; + track(3,i) = y_ + py_*length_; + track(5,i) = z_ + length/bet0 - (one/bet0+pt_)*length_; + endif + enddo endif end subroutine trsol From 18605e4820de74de504083efb99f646dd54fa635 Mon Sep 17 00:00:00 2001 From: Andrea Date: Mon, 6 Nov 2017 23:42:28 +0100 Subject: [PATCH 5/6] Thick solenoid: if strength == 0 then it just calls ttdft() --- src/trrun.f90 | 61 ++++++++++++++++++++------------------------------- 1 file changed, 24 insertions(+), 37 deletions(-) diff --git a/src/trrun.f90 b/src/trrun.f90 index d514b8820..7aa153ddc 100644 --- a/src/trrun.f90 +++ b/src/trrun.f90 @@ -3176,11 +3176,8 @@ subroutine trsol(track,ktrack) double precision :: x_, y_, z_, px_, py_, pt_ !---- Initialize. - ! dtbyds = get_value('probe ','dtbyds ') - ! gamma = get_value('probe ','gamma ') bet0 = get_value('probe ','beta ') - ! deltap = get_value('probe ','deltap ') - ! + !---- Get solenoid parameters ! elrad = node_value('lrad ') bvk = node_value('other_bv ') @@ -3222,33 +3219,30 @@ subroutine trsol(track,ktrack) track(5,i) = (sigf + (xf*pyf - yf*pxf)*Z) / bet0 ! track(6,i) = psigf*bet0 enddo - else - - 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_); - - if (skl.ne.zero) then - - ! set up constants + 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 + ! 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; @@ -3257,17 +3251,10 @@ subroutine trsol(track,ktrack) 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; - - else - - length_ = length/sqrt(onedp**2-px_**2-py_**2); ! length/pz - track(1,i) = x_ + px_*length_; - track(3,i) = y_ + py_*length_; - track(5,i) = z_ + length/bet0 - (one/bet0+pt_)*length_; - - endif - - enddo + enddo + else + call ttdrf(length,track,ktrack); + endif endif end subroutine trsol From 814ae67c2476d15a03f35f4c4c944f902193d0a4 Mon Sep 17 00:00:00 2001 From: Andrea Date: Tue, 7 Nov 2017 00:14:01 +0100 Subject: [PATCH 6/6] Preparing for thorough test with mkthin --- tests/test-thick-solenoid/test-thick-solenoid.madx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test-thick-solenoid/test-thick-solenoid.madx b/tests/test-thick-solenoid/test-thick-solenoid.madx index 2d7909e1d..b77925c8b 100644 --- a/tests/test-thick-solenoid/test-thick-solenoid.madx +++ b/tests/test-thick-solenoid/test-thick-solenoid.madx @@ -27,6 +27,9 @@ 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;