Merge pull request #924 from tpersson/radiationIntegralsFinal

Radiation integrals final

doc/latexuguide/biblio-mad.bib

@@ -538,11 +538,31 @@ @TechReport{lebedev2010
   note =         {\url{http://cern.ch/madx/doc/1748-0221_5_10_P10010.pdf}}}
 
 @TechReport{kapin2013,
-	author        = "Valery Kapin and Frank Schmidt",
-	title         = "{MADX-SC Flag Description}",
-	month         = "Nov",
-	year          = "2013",
-	number        = {CERN-ACC-NOTE-2013-0036}}
+  author        = "Valery Kapin and Frank Schmidt",
+  title         = "{MADX-SC Flag Description}",
+  month         = "Nov",
+  year          = "2013",
+  number        = {CERN-ACC-NOTE-2013-0036}}
+
+@article{Helm1973,
+  author = {R. H. Helm and M. J. Lee and P.L. Morton and M. Sands},
+  title = "{Evaluation of synchrotron radiation integrals}",
+  number = {SLAC-PUB-1193},
+  journal = {IEEE Trans. Nucl. Sci.},
+  pages = {900-901},
+  year = {1973},
+  note =         {\href{https://www.slac.stanford.edu/pubs/slacpubs/1000/slac-pub-1193.pdf}{PDF}}}
+}
+
+@article{Jowett1986,
+  author        = {J. M. Jowett},
+  title         = {Introductory statistical mechanics for electron storage rings},
+  journal       = {AIP Conf. Proc.},
+  number        = {SLAC-PUB-4033},
+  pages         = {864-970},
+  year          = {1986},
+  note =         {\href{http://cds.cern.ch/record/175614}{PDF}}}
+}
 
 %% EOF
 

doc/latexuguide/conventions.tex

@@ -740,7 +740,9 @@ \subsection{Variables in the SUMM Table}
 	\ttitem{SYNCH\_2} Second synchrotron radiation integral  
 	\ttitem{SYNCH\_3} Third synchrotron radiation integral  
 	\ttitem{SYNCH\_4} Fourth synchrotron radiation integral  
-	\ttitem{SYNCH\_5} Fifth synchrotron radiation integral  
+	\ttitem{SYNCH\_5} Fifth synchrotron radiation integral    
+	\ttitem{SYNCH\_6} Sixth synchrotron radiation integral  
+	\ttitem{SYNCH\_8} Eighth synchrotron radiation integral  
 \end{madlist}
 
 Notice that in \madx, PT substitutes DELTAP as longitudinal
@@ -751,6 +753,24 @@ \subsection{Variables in the SUMM Table}
 number of time equal to the order of the derivative to find the
 functions given in the literature. 
 
+The first five radiation integrals, SYNCH\_1 through to SYNCH\_5, 
+are defined and computed using the expressions in~\cite{Helm1973}. 
+It should be noted that these integrals assume that the beam traverses 
+through the centre of each element and does not take into account effects 
+from non-zero closed orbits, including orbits resulting from chromatic effects. 
+These integrals are therefore useful for estimating various radiation 
+properties of an ideal machine but the EMIT module should be used for 
+a more general case. 
+
+The integrals six and eight,  SYNCH\_6 and SYNCH\_8, as defined 
+in~\cite{Jowett1986} have been implemented to capture phenomena 
+not considered by the first five integrals. Integral six is required 
+for computing the radiation in the quadrupoles and integral eight 
+can be used to understand how the partition function varies for off 
+momentum beams as detailed in~\cite{Jowett1986}. Again, these 
+estimates should be used with caution and the EMIT module should 
+be used for precise computations.
+
 \subsection{Variables in the TRACK Table}
 \label{subsec:tables-track}
 The command RUN writes tables with the following variables: 

src/mad_gcst.c

@@ -645,9 +645,11 @@ const int summ_table_types[] =
   2, 2, 2, 2, 2,
   2, 2, 2, 2, 2,
   2, 2, 2, 2, 2,
-  2, 2, 2, 2, 2,
+  2, 2, 2, 2, 
   //2, 2, 2, 2, 2,
-  2, 2, 2, 2, 2, //for nflips
+  2, 2, 2, 2, 2,
+  2, 2,
+  2, //for nflips
 };
 
 const char* const summ_table_cols[] =
@@ -657,7 +659,9 @@ const char* const summ_table_cols[] =
   "xcorms", "q2", "dq2", "betymax", "dymax",
   "dyrms", "ycomax", "ycorms", "deltap",
   //"synch_1","synch_2","synch_3","synch_4","synch_5",
-  "synch_1","synch_2","synch_3","synch_4","synch_5", "nflips", //for nflips
+  "synch_1","synch_2","synch_3","synch_4","synch_5",
+  "synch_6","synch_8",
+   "nflips", //for nflips
   " "  /* blank terminates */
 };
 

src/twiss.f90

@@ -67,6 +67,7 @@ SUBROUTINE twiss(rt,disp0,tab_name,sector_tab_name)
   mode_flip  =.false.
 
   synch_1=zero;  synch_2=zero;  synch_3=zero;  synch_4=zero;  synch_5=zero
+  synch_6=zero;  synch_8=zero
 
   suml=zero; circ=zero; eta=zero; alfa=zero; gamtr=zero; wgt=zero
 
@@ -2962,6 +2963,7 @@ SUBROUTINE twchgo
 
   !---- and synchrotron radiation integrals
   synch_1 = zero; synch_2 = zero; synch_3 = zero; synch_4 = zero; synch_5 = zero
+  synch_6=zero;  synch_8=zero
 
   !---- Loop over positions.
   i = restart_sequ()
@@ -3084,21 +3086,21 @@ SUBROUTINE tw_synch_int()
   use twisslfi
   use twisscfi
   use twissbeamfi, only : beta
-  use math_constfi, only : two
+  use math_constfi, only : zero, two
   implicit none
   !----------------------------------------------------------------------*
   !     Purpose:                                                         *
   !     Compute and add synchrotron radiation integral                   *
   !----------------------------------------------------------------------*
   double precision :: an, e1, e2, sk1, rhoinv, blen
-  double precision :: syncint(5)
+  double precision :: syncint(8)
 
   double precision, external :: node_value
 
   !---- Initialisation
   blen = node_value('blen ')
   rhoinv = node_value('rhoinv ')
-  sk1 = node_value('k1 ')
+  sk1 = node_value('k1 ') + node_value('k1tap ')
   e1 = node_value('e1 ')
   e2 = node_value('e2 ')
   an = node_value('angle ')
@@ -3108,17 +3110,20 @@ SUBROUTINE tw_synch_int()
   endif
 
   !---- Synchrotron radiation integrals through bending magnets.
-  if (rhoinv .ne. 0.d0) then
+
      ! Note that calcsyncint expects dx and dpx as derivatives wrt deltap.
      ! since MAD take disp(1) and disp(2) as derivatives wrt pt, they must be
      ! multiplied by beta before the call to calcsyncint.
+     syncint = zero
      call calcsyncint(rhoinv,blen,sk1,e1,e2,betx,alfx,disp(1)*beta,disp(2)*beta,syncint)
      synch_1 = synch_1 + syncint(1)
      synch_2 = synch_2 + syncint(2)
      synch_3 = synch_3 + syncint(3)
      synch_4 = synch_4 + syncint(4)
      synch_5 = synch_5 + syncint(5)
-  endif
+     synch_6 = synch_6 + syncint(6)
+     synch_8 = synch_8 + syncint(8)
+
 end subroutine tw_synch_int
 
 
@@ -3382,6 +3387,8 @@ SUBROUTINE tw_summ(rt,tt)
   call double_to_table_curr('summ ','synch_3 ' ,synch_3)
   call double_to_table_curr('summ ','synch_4 ' ,synch_4)
   call double_to_table_curr('summ ','synch_5 ' ,synch_5)
+  call double_to_table_curr('summ ','synch_6 ' ,synch_6)
+  call double_to_table_curr('summ ','synch_8 ' ,synch_8)
 
 end SUBROUTINE tw_summ
 
@@ -8976,7 +8983,7 @@ SUBROUTINE tmrfmult(fsec,ftrk,orbit,fmap,ek,re,te)
 end SUBROUTINE tmrfmult
 
 SUBROUTINE calcsyncint(rhoinv,blen,k1,e1,e2,betxi,alfxi,dxi,dpxi,I)
-  use name_lenfi
+  use math_constfi, only : zero, one, two
   implicit none
   !----------------------------------------------------------------------*
   !     Purpose:                                                         *
@@ -9017,15 +9024,16 @@ SUBROUTINE calcsyncint(rhoinv,blen,k1,e1,e2,betxi,alfxi,dxi,dpxi,I)
   !----------------------------------------------------------------------*
 
   double precision, intent(IN) :: rhoinv, blen, k1, e1, e2, betxi, alfxi, dxi, dpxi
-  double precision, intent(OUT) :: I(5)
+  double precision, intent(OUT) :: I(8)
 
   ! local variables
-  double precision :: dx2, gamx, dispaverage, curlyhaverage
-  double precision :: betx, alfx, dx, dpx
+  double precision :: dx2, gamx, dispaverage, curlyhaverage, lq
+  double precision :: betx, alfx, dx, dpx, u0x, u1x, u2x
+  double precision :: gammai, betxaverage, k1n
   double complex :: k2, k, kl
 
   integer, external :: get_option
-  character(len=name_len) :: name
+  double precision, external :: node_value
 
   betx = betxi
   dx = dxi
@@ -9056,26 +9064,54 @@ SUBROUTINE calcsyncint(rhoinv,blen,k1,e1,e2,betxi,alfxi,dxi,dpxi,I)
                        gamx*(3*kl - 4*sin(kl) + sin(kl)*cos(kl))/(2*k2*kl**3) &
                      - alfx*(1-cos(kl))**2/(k*kl**3) &
                      + betx*(kl-cos(kl)*sin(kl))/(2*kl**3)))
-
+  if (rhoinv .ne. 0.d0) then
   I(1) = dispaverage * rhoinv * blen
   I(2) = rhoinv*rhoinv * blen
   I(3) = abs(rhoinv)**3 * blen
   I(4) = dispaverage*rhoinv*(rhoinv**2 + 2*k1) * blen &
            - rhoinv*rhoinv*(dx*tan(e1) + dx2*tan(e2))
   I(5) = curlyhaverage * abs(rhoinv)**3 * blen
+  end if
+
+  if(k1 .ne. 0.d0) then
+    lq = node_value('l ')
+    gammai = (one+alfxi*alfxi)/betxi;
+    dx2 = real(dx*cos(kl) + dpx*sin(kl)/k)
+    dispaverage = real(dx * sin(kl)/kl &
+             + dpxi * (1 - cos(kl))/(k*kl))
+  if(k1 .ge. zero) then
+    k1n = k1
+    u0x = (one + sin(two*sqrt(k1n)*lq)/(two*sqrt(k1n)*lq))/two
+    u1x = sin(sqrt(k1n)*lq)**two/(k1n*lq)
+    u2x = (one - sin(two*sqrt(k1n)*lq)/(two*sqrt(k1n)*lq))/(two*k1n)
+    dx2 = cos(sqrt(k1n)*lq)*dxi + (one/sqrt(k1n))*sin(sqrt(k1n)*lq)*dpxi
+    dispaverage = (dxi+dx2)/two
+  else
+    k1n = -k1
+    u0x = (one + sinh(two*sqrt(k1n)*lq)/(two*sqrt(k1n)*lq))/two
+    u1x = sinh(sqrt(k1n)*lq)**two/(k1n*lq)
+    u2x = -(one - sinh(two*sqrt(k1n)*lq)/(two*sqrt(k1n)*lq))/(two*k1n)
+    dx2 = cosh(sqrt(k1n)*lq)*dxi + (one/sqrt(k1n))*sinh(sqrt(k1n)*lq)*dpxi
+    dispaverage = (dxi+dx2)/two
+  endif
+
+    betxaverage = betxi*u0x - alfxi*u1x  + gammai*u2x
+
+    I(6) = (k1n**2)*betxaverage*lq
+    I(8) = (k1n**2)*dispaverage**2*lq
+
+  endif
 
   if (get_option('debug ') .ne. 0) then
-    !LD: 22.12.2019
-    call element_name(name,len(name))
-    print *, 'Synchrotron integrals at exit of element ', name
-    print *, 'Input:  rhoinv = ', rhoinv, 'k1 = ', k1, 'e1 =', e1, 'e2 = ', e2, 'blen = ', blen
-    print *, '        betxi = ', betxi, 'alfxi = ', alfxi, 'dxi = ', dxi, 'dpxi = ', dpxi
-    print *, ' --> '
-    print *, '        k2 = ', k2, '  k = ', k, 'k*l = ', kl
-    print *, '        alfx = ', alfx, 'dpx = ', dpx, 'gamx = ', gamx, 'dx2 = ', dx2
-    print *, '        dispaverage = ', dispaverage, 'curlyhaverage = ', curlyhaverage
-    print *, 'Contributions to Radiation Integrals:', I(1), I(2), I(3), I(4), I(5)
-    print *, ' '
+     print *, ' '
+     print *, 'Input:  rhoinv = ', rhoinv, 'k1 = ', k1, 'e1 =', e1, 'e2 = ', e2, 'blen = ', blen
+     print *, '        betxi = ', betxi, 'alfxi = ', alfxi, 'dxi = ', dxi, 'dpxi = ', dpxi
+     print *, ' --> '
+     print *, '        k2 = ', k2, '  k = ', k, 'k*l = ', kl
+     print *, '        alfx = ', alfx, 'dpx = ', dpx, 'gamx = ', gamx, 'dx2 = ', dx2
+     print *, '        dispaverage = ', dispaverage, 'curlyhaverage = ', curlyhaverage
+     print *, 'Contributions to Radiation Integrals:', I(1), I(2), I(3), I(4), I(5)
+     print *, ' '
   endif
 
 END SUBROUTINE calcsyncint

src/util.f90

@@ -380,6 +380,7 @@ module twisscfi
   double precision :: wgt=0.d0, suml=0.d0, circ=0.d0, eta=0.d0, alfa=0.d0, gamtr=0.d0
   double precision :: wx=0.d0, phix=0.d0, dmux=0.d0, xix=0.d0, wy=0.d0, phiy=0.d0, dmuy=0.d0, xiy=0.d0
   double precision :: synch_1=0.d0, synch_2=0.d0, synch_3=0.d0, synch_4=0.d0, synch_5=0.d0
+  double precision :: synch_6=0.d0, synch_7=0.d0, synch_8=0.d0
   double precision :: gammacp=1.d0
   integer :: nmode_flip=0
 end module twisscfi

testing/updateFromList.py

@@ -31,7 +31,7 @@
 
 #copyfile(src, dst)
 mypath_orig = '../tests/'
-file1 = open('failedTests.txt', 'r') 
+file1 = open('../mynewfail.txt', 'r') 
 Lines = file1.readlines() 
 loc = os.getcwd()
 count = 0
@@ -57,4 +57,4 @@
 			cpname = f[:-4] + ".ref"
 			copyfile(f, cpname)
 
-	#print(onlyfiles)
+	#print(onlyfiles)

tests/test-aperture-2/ptc-twiss-table.ref

@@ -84,8 +84,8 @@
 @ PTCOMAX          %le                   0
 @ TITLE            %58s "LHC Version 6.503 -  July 2009 - version with RCOLLIMATORs"
 @ ORIGIN           %16s "5.05.02 Linux 64"
-@ DATE             %08s "20/07/20"
-@ TIME             %08s "12.12.52"
+@ DATE             %08s "29/07/20"
+@ TIME             %08s "11.42.25"
 * NAME                                S               BETX               BETY                  X                 PX                  Y                 PY APERTYPE                       APER_1            ON_ELEM 
 $ %s                                %le                %le                %le                %le                %le                %le                %le %s                                %le                %le 
  "TAS.1R1"                       20.915        50.76701663        50.76706051    -0.002000016125   -7.053266099e-10     0.003555550365    0.0001700000459 "CIRCLE"                             0                  0

tests/test-aperture-2/test-aperture-2.ref

@@ -3,7 +3,7 @@
   +     MAD-X 5.05.02  (64 bit, Linux)       +
   + Support: mad@cern.ch, http://cern.ch/mad +
   + Release   date: 2019.07.25               +
-  + Execution date: 2020.07.20 12:12:40      +
+  + Execution date: 2020.07.29 11:41:51      +
   ++++++++++++++++++++++++++++++++++++++++++++
 TITLE, "LHC Version 6.503 -  July 2009 - version with RCOLLIMATORs" ;
 
@@ -104,8 +104,8 @@ orbit:  -2.000002E-03 -6.829939E-10  3.774260E-11  1.700000E-04  0.000000E+00  0
            synch_2            synch_3            synch_4            synch_5 
                  0                  0                  0                  0 
 
-            nflips 
-                 0 
+           synch_6            synch_8             nflips 
+                 0                  0                  0 
 
 
 ! Initialize PTC

tests/test-aperture-2/twiss-table.ref

@@ -42,8 +42,8 @@
 @ SYNCH_5          %le                   0
 @ TITLE            %58s "LHC Version 6.503 -  July 2009 - version with RCOLLIMATORs"
 @ ORIGIN           %16s "5.05.02 Linux 64"
-@ DATE             %08s "20/07/20"
-@ TIME             %08s "12.12.41"
+@ DATE             %08s "29/07/20"
+@ TIME             %08s "11.41.55"
 * NAME                                S               BETX               BETY                  X                 PX                  Y                 PY APERTYPE                       APER_1            ON_ELEM 
 $ %s                                %le                %le                %le                %le                %le                %le                %le %s                                %le                %le 
  "TAS.1R1"                       20.915        50.76701356        50.76701003    -0.002000016307   -6.829939483e-10     0.003555549899    0.0001699999934 "CIRCLE"                             0                  0

tests/test-aperture-3/ap.tfs.ref

@@ -26,8 +26,8 @@
 @ AT_ELEMENT       %12s  "COL.PU_QF4:2"
 @ TITLE            %44s "AD HE 2010 optics. Anti-Protons - 3.57 GeV/c"
 @ ORIGIN           %16s "5.05.02 Linux 64"
-@ DATE             %08s "20/04/20"
-@ TIME             %08s "10.26.25"
+@ DATE             %08s "29/07/20"
+@ TIME             %08s "11.42.26"
 * NAME                                S               BETX               BETY                 DX                  X                  Y                 N1 APERTYPE                       APER_1             APER_2             APER_3             APER_4               XTOL               YTOL               RTOL 
 $ %s                                %le                %le                %le                %le                %le                %le                %le %s                                %le                %le                %le                %le                %le                %le                %le 
  "AD$START"                           0        3.850951181                  1       0.1151944018                  0                  0             999999 "CIRCLE"                             0                  0                  0                  0                  0                  0                  0

tests/test-aperture-3/test-aperture-3.ref

@@ -1,9 +1,9 @@
 
   ++++++++++++++++++++++++++++++++++++++++++++
-  +     MAD-X 5.03.07  (64 bit, Darwin)      +
+  +     MAD-X 5.05.02  (64 bit, Linux)       +
   + Support: mad@cern.ch, http://cern.ch/mad +
-  + Release   date: 2017.10.20               +
-  + Execution date: 2017.12.14 18:01:18      +
+  + Release   date: 2019.07.25               +
+  + Execution date: 2020.07.29 11:42:26      +
   ++++++++++++++++++++++++++++++++++++++++++++
 
 
@@ -113,8 +113,8 @@ final orbit vector:     0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
            synch_2            synch_3            synch_4            synch_5 
                  0                  0                  0                  0 
 
-            nflips 
-                 0 
+           synch_6            synch_8             nflips 
+                 0                  0                  0 
 
 
 

tests/test-aperture/ap.tfs.ref

@@ -26,8 +26,8 @@
 @ AT_ELEMENT       %05s  "MQ1:1"
 @ TITLE            %25s "LHC Cell APERTURE example"
 @ ORIGIN           %16s "5.05.02 Linux 64"
-@ DATE             %08s "20/04/20"
-@ TIME             %08s "10.05.02"
+@ DATE             %08s "29/07/20"
+@ TIME             %08s "11.41.51"
 * NAME                       S      BETX      BETY        DX         X         Y        N1 APERTYPE              APER_1    APER_2    APER_3    APER_4      XTOL      YTOL      RTOL 
 $ %s                       %le       %le       %le       %le       %le       %le       %le %s                       %le       %le       %le       %le       %le       %le       %le 
  "LHCCELL$START"       0.00000 177.45326  32.26597   2.18636   0.00000   0.00000 999999.00000 "CIRCLE"              0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000