Merge pull request #495 from alatina/master

Implemented thick solenoid tracking in module TRACK
MethodicalAcceleratorDesign · Nov 7, 2017 · 3b48a75 · 3b48a75
2 parents eeb5e1c + 814ae67
commit 3b48a75
Show file tree

Hide file tree

Showing 17 changed files with 294 additions and 54 deletions.
diff --git 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
@@ -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"
@@ -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:                                                             *
@@ -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)
@@ -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)
 

diff --git a/tests/test-thick-solenoid/madx-thick.obs0001.p0001.cfg 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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;
+