Skip to content

Commit

Permalink
Merge pull request #609 from MethodicalAcceleratorDesign/ld-dev
Browse files Browse the repository at this point in the history
Ld dev: integration of synrad from H. Burkhardt
  • Loading branch information
madcern committed Mar 27, 2018
2 parents 561d4e7 + af3ee07 commit b7f4133
Show file tree
Hide file tree
Showing 6 changed files with 140 additions and 29 deletions.
7 changes: 7 additions & 0 deletions doc/latexuguide/biblio-mad.bib
Original file line number Diff line number Diff line change
Expand Up @@ -491,6 +491,13 @@ @Article{roy1990
volume = {A298},
pages = {128-133}}

@TechReport{hbu2007,
author = "Helmut Burkhardt",
title = "{Monte Carlo Generation of the Energy Spectrum of Synchrotron Radiation}",
month = "Jun",
year = "2007",
number = {CERN-OPEN-2007-018}}

@TechReport{kapin2006,
author = {V. Kapin and F. Schmidt},
title = {{PTC modules for MAD-X code}},
Expand Down
2 changes: 1 addition & 1 deletion doc/latexuguide/format.tex
Original file line number Diff line number Diff line change
Expand Up @@ -352,7 +352,7 @@ \subsubsection{Random Number Generators}

The command \hyperref[sec:option]{\texttt{OPTION}} has two attributes to control the new random generator of MAD-X. The attribute \texttt{RAND} can take the value \texttt{default} (for backward compatibility) or \texttt{best} for the new generator. The latter supports 10 multiple independent generators numbered from 1 to 10, and can be selected with the attribute \texttt{RANDID=id} in combination to \texttt{RAND}. If \texttt{RANDID} is not specified and \texttt{RAND=best}, then the last \texttt{id} set is used (init: 1). The user can switch between the 10 independent generators at will to use different streams for the next generated random numbers. When a new seed is set through to command \hyperref[sec:eoption]{\texttt{EOPTION}}, \hyperref[sec:coption]{\texttt{COPTION}} or \hyperref[sec:track]{\texttt{TRACK}}, it will only affect the currently selected stream. The special initial seed for the \texttt{best} generator can be restored by using the value zero instead of the commands default 123456789.

The following example displays the first 10 random numbers of the 10 streams of the \texttt{best} generator, using the intial seed.
The following example displays the first 10 random numbers of the 10 streams of the \texttt{best} generator, using the default initial seed.

\begin{verbatim}
option, -info;
Expand Down
5 changes: 4 additions & 1 deletion doc/latexuguide/thintrack.tex
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,10 @@ \section{TRACK}
and flag \texttt{RADIATE} in the \texttt{BEAM} command). \\ (Default:~false)

\ttitem{QUANTUM} flag to introduce quantum excitation via random
number generator and tables for photon emission. \\ (Default:~false)
number generator and tables look-up (\texttt{SYNRAD} $=1$, see ref. \cite{roy1990}) or polynomial
interpolation (\texttt{SYNRAD} $=2$, see ref. \cite{hbu2007}) for photon emission.
The choice of the generator can be selected via the command
\hyperref[sec:option]{\texttt{OPTION}} attribute \texttt{SYNRAD}. \\ (Default:~2)

\ttitem{SEED} If \texttt{QUANTUM} is true, it selects a particular sequence of random values.
A \texttt{SEED} value is an integer in the range [0...999999999] (default:
Expand Down
1 change: 1 addition & 0 deletions src/mad_dict.c
Original file line number Diff line number Diff line change
Expand Up @@ -276,6 +276,7 @@ const char *const_command_def =
"rbarc = [l, true, true], "
"rand = [s, default, best], "
"randid = [i, 0, 1], "
"synrad = [i, 2, 2], " /* ld */
"cfsbend = [l, true, true], " /* combined function sbend */
"thin_foc = [l, true, true], "
"sympl = [l, true, true], "
Expand Down
85 changes: 85 additions & 0 deletions src/mad_synrad.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
// mad_synrad.cpp SR photon spectrum generator H.Burkhardt CERN-OPEN-2007-018
// also used in Geant4 G4SynchrotronRadiation
// modified here to make routine name and argument FORTRAN call compatible

#include <cmath>

inline double Chebyshev(double a,double b,const double c[],int m,double x)
{
double y;
double y2=2.0*(y=(2.0*x-a-b)/(b-a)); // Change of variable.
double d=0,dd=0,sv;
for (int j=m-1;j>=1;j--) // Clenshaw's recurrence.
{
sv=d;
d=y2*d-dd+c[j];
dd=sv;
}
return y*d-dd+0.5*c[0];
}

extern "C" double invsynfracint_(double *fortx)
{
double x=fortx[0];
// from 0 to 0.7
static const double aa1=0 ,aa2=0.7;
static const int ncheb1=27;
static const double cheb1[] =
{
1.22371665676046468821,0.108956475422163837267,0.0383328524358594396134,0.00759138369340257753721,
0.00205712048644963340914,0.000497810783280019308661,0.000130743691810302187818,0.0000338168760220395409734,
8.97049680900520817728e-6,2.38685472794452241466e-6,6.41923109149104165049e-7,1.73549898982749277843e-7,
4.72145949240790029153e-8,1.29039866111999149636e-8,3.5422080787089834182e-9,9.7594757336403784905e-10,
2.6979510184976065731e-10,7.480422622550977077e-11,2.079598176402699913e-11,5.79533622220841193e-12,
1.61856011449276096e-12,4.529450993473807e-13,1.2698603951096606e-13,3.566117394511206e-14,1.00301587494091e-14,
2.82515346447219e-15,7.9680747949792e-16};
// from 0.7 to 0.9132260271183847
static const double aa3=0.9132260271183847;
static const int ncheb2=27;
static const double cheb2[] =
{
1.1139496701107756,0.3523967429328067,0.0713849171926623,0.01475818043595387,0.003381255637322462,
0.0008228057599452224,0.00020785506681254216,0.00005390169253706556,0.000014250571923902464,3.823880733161044e-6,
1.0381966089136036e-6,2.8457557457837253e-7,7.86223332179956e-8,2.1866609342508474e-8,6.116186259857143e-9,
1.7191233618437565e-9,4.852755117740807e-10,1.3749966961763457e-10,3.908961987062447e-11,1.1146253766895824e-11,
3.1868887323415814e-12,9.134319791300977e-13,2.6211077371181566e-13,7.588643377757906e-14,2.1528376972619e-14,
6.030906040404772e-15,1.9549163926819867e-15};
// Chebyshev with exp/log scale
// a = -Log[1 - SynFracInt[1]]; b = -Log[1 - SynFracInt[7]];
static const double aa4=2.4444485538746025480,aa5=9.3830728608909477079;
static const int ncheb3=28;
static const double cheb3[] =
{
1.2292683840435586977,0.160353449247864455879,-0.0353559911947559448721,0.00776901561223573936985,
-0.00165886451971685133259,0.000335719118906954279467,-0.0000617184951079161143187,9.23534039743246708256e-6,
-6.06747198795168022842e-7,-3.07934045961999778094e-7,1.98818772614682367781e-7,-8.13909971567720135413e-8,
2.84298174969641838618e-8,-9.12829766621316063548e-9,2.77713868004820551077e-9,-8.13032767247834023165e-10,
2.31128525568385247392e-10,-6.41796873254200220876e-11,1.74815310473323361543e-11,-4.68653536933392363045e-12,
1.24016595805520752748e-12,-3.24839432979935522159e-13,8.44601465226513952994e-14,-2.18647276044246803998e-14,
5.65407548745690689978e-15,-1.46553625917463067508e-15,3.82059606377570462276e-16,-1.00457896653436912508e-16};
static const double aa6=33.122936966163038145;
static const int ncheb4=27;
static const double cheb4[] =
{
1.69342658227676741765,0.0742766400841232319225,-0.019337880608635717358,0.00516065527473364110491,
-0.00139342012990307729473,0.000378549864052022522193,-0.000103167085583785340215,0.0000281543441271412178337,
-7.68409742018258198651e-6,2.09543221890204537392e-6,-5.70493140367526282946e-7,1.54961164548564906446e-7,
-4.19665599629607704794e-8,1.13239680054166507038e-8,-3.04223563379021441863e-9,8.13073745977562957997e-10,
-2.15969415476814981374e-10,5.69472105972525594811e-11,-1.48844799572430829499e-11,3.84901514438304484973e-12,
-9.82222575944247161834e-13,2.46468329208292208183e-13,-6.04953826265982691612e-14,1.44055805710671611984e-14,
-3.28200813577388740722e-15,6.96566359173765367675e-16,-1.294122794852896275e-16};

if(x<aa2) return x*x*x*Chebyshev(aa1,aa2,cheb1,ncheb1,x);
else if(x<aa3) return Chebyshev(aa2,aa3,cheb2,ncheb2,x);
else if(x<1-0.0000841363)
{
double y=-log(1-x);
return y*Chebyshev(aa4,aa5,cheb3,ncheb3,y);
}
else
{
double y=-log(1-x);
return y*Chebyshev(aa5,aa6,cheb4,ncheb4,y);
}
}

69 changes: 42 additions & 27 deletions src/trrun.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1254,7 +1254,7 @@ subroutine ttyrot(track,ktrack)
py = TRACK(4,i)
t = TRACK(5,i)
pt = TRACK(6,i)

pz = sqrt(one + two*pt/bet0i + pt**2 - px**2 - py**2)
ptt = 1 - ta*px/pz
track(1,i) = x/(ca*ptt)
Expand Down Expand Up @@ -3239,12 +3239,12 @@ subroutine trsol(track,ktrack)
cosTh = cos(two*skl/onedp)
sinTh = sin(two*skl/onedp)
omega = sk/onedp;

! total path length traveled by the particle
bet = onedp / (one/bet0 + pt_);
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;
Expand Down Expand Up @@ -4378,12 +4378,12 @@ subroutine tttquad(track, ktrack)
else
tilt = zero
endif

if (k1.eq.zero) then
call ttdrf(length,track,ktrack);
return
endif

!---- Prepare to calculate the kick and the matrix elements
do jtrk = 1, ktrack
!---- The particle position
Expand Down Expand Up @@ -4541,8 +4541,8 @@ subroutine tttdipole(track, ktrack)
n_ferr = node_fd_errors(f_errors)
if (k0.ne.0) f_errors(0) = f_errors(0) + k0*length - angle;
k0 = k0 + f_errors(0) / length ! dipole term
k1 = k1 + f_errors(2) / length ! quad term
k1 = k1 + f_errors(2) / length ! quad term

!---- Apply entrance dipole edge effect
if (node_value('kill_ent_fringe ') .eq. zero) &
call ttdpdg_map(track, ktrack, e1, h1, hgap, fint, zero)
Expand Down Expand Up @@ -4641,12 +4641,18 @@ subroutine trphot(el,curv,rfac,deltap)
implicit none
!----------------------------------------------------------------------*
! Purpose: *
! option synrad=1 *
! Generate random energy loss for photons, using a look-up table to *
! invert the function Y. Ultra-basic interpolation computed; *
! leads to an extrapolation outside the table using the two outmost *
! point on each side (low and high). *
! Assumes ultra-relativistic particles (beta = 1). *
! Author: Ghislain Roy *
! Author: Ghislain Roy *
! option synrad=2 *
! Generate random energy loss for photons, using the *
! Synchrotron radiation spectrum generator *
! described in CERN-OPEN-2007-018 by Helmut Burkhardt * *

! Input: *
! EL (double) Element length. *
! CURV (double) Local curvature of orbit. *
Expand All @@ -4660,12 +4666,13 @@ subroutine trphot(el,curv,rfac,deltap)
double precision :: el, curv, rfac, deltap

integer :: i, ierror, j, nphot
double precision :: amean, dlogr, slope, ucrit, xi, sumxi
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
character(len=20) text

integer, external :: get_option
integer :: synrad
double precision :: get_value, frndm
double precision, parameter :: fac1=3.256223d0
! integer, parameter :: maxran=1000000000
Expand Down Expand Up @@ -4753,23 +4760,31 @@ subroutine trphot(el,curv,rfac,deltap)

!---- For all photons, sum the radiated photon energy in units of UCRIT
if (nphot .ne. 0) then
do i = 1, nphot

!---- Find a uniform random number in the range [ 0,3.256223 ].
! Note that the upper limit is not exactly 15*sqrt(3)/8
! because of imprecision in the integration of F.
dlogr = log(fac1 * frndm())
!---- Now look for the energy of the photon in the table TABY/TABXI
do j = 2, maxtab
if (dlogr .le. taby(j) ) go to 20
enddo

!---- Perform linear interpolation and sum up energy lost.
20 slope = (dlogr - taby(j-1)) / (taby(j) - taby(j-1))
xi = dexp(tabxi(j-1) + slope * (tabxi(j) - tabxi(j-1)))
sumxi = sumxi + xi
! write(60,*) xi ! 2016-Mar-16 18:22:30 ghislain: dump individual photons to file
enddo
synrad = get_option('synrad ')
if (synrad .eq. 1) then
do i = 1, nphot
!---- Find a uniform random number in the range [ 0,3.256223 ].
! Note that the upper limit is not exactly 15*sqrt(3)/8
! because of imprecision in the integration of F.
dlogr = log(fac1 * frndm())
!---- Now look for the energy of the photon in the table TABY/TABXI
do j = 2, maxtab
if (dlogr .le. taby(j) ) go to 20
enddo
!---- Perform linear interpolation and sum up energy lost.
20 slope = (dlogr - taby(j-1)) / (taby(j) - taby(j-1))
xi = dexp(tabxi(j-1) + slope * (tabxi(j) - tabxi(j-1)))
sumxi = sumxi + xi
! write(60,*) xi ! 2016-Mar-16 18:22:30 ghislain: dump individual photons to file
enddo
else !-- synrad .eq. 2
do i = 1, nphot
rn_arg = frndm()
xi = InvSynFracInt(rn_arg)
!hbutest print *,xi,rn_arg," InvSynFracIntdata"
sumxi = sumxi + xi
enddo
endif
endif
endif

Expand Down

0 comments on commit b7f4133

Please sign in to comment.