Commit ccebd9e7 authored by gsell's avatar gsell
Browse files

Merge branch '461-cleanup-beam-envelope-tracker' into 'master'

Resolve "cleanup: beam envelope tracker"

Closes #461

See merge request !273
parents 6e19d56c 8c42e1bf
......@@ -20,15 +20,10 @@ set (_SRCS
bet/EnvelopeSlice.cpp
bet/EnvelopeBunch.cpp
bet/profile.cpp
bet/math/integrate.cpp
bet/math/interpol.cpp
bet/math/rk.cpp
bet/math/functions.cpp
bet/math/linfit.cpp
bet/math/root.cpp
bet/math/savgol.cpp
bet/math/svdfit.cpp
bet/BetError.cpp
)
include_directories (
......@@ -56,18 +51,13 @@ set (HDRS
ThickMapper.h
ThickTracker.h
TransportMapper.h
bet/BetError.h
bet/EnvelopeBunch.h
bet/EnvelopeSlice.h
bet/profile.h
bet/math/functions.h
bet/math/integrate.h
bet/math/interpol.h
bet/math/linfit.h
bet/math/rk.h
bet/math/root.h
bet/math/savgol.h
bet/math/svdfit.h
)
install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Algorithms")
......
/* error.C
error handling
Project: Beam Envelope Tracker (BET)
Revision history
Date Description Programmer
------------ -------------------------------------------- --------------
07-03-06 Created Rene Bakker
Last Revision:
$Id: error.C 91 2007-05-06 17:47:45Z bakker $
*/
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <stdarg.h>
#include <errno.h>
#include "Algorithms/bet/BetError.h"
void initBetErrorMsg(int level, const char *fbase) {
printf("\n\n error.cpp: initErrorMesg \n\n");
} /* initErrorMsg() */
void initBetErrorFilename(const char *fbase) {
printf("\n\n error.cpp: initErrorFilename \n\n");
} /* initErrorFileName() */
void setBetReportLevel(int level) {
printf("\n\n error.cpp: setReportLevel \n\n");
}
void writeBetError(std::string err) {
std::cout << "\n\n error.cpp: writeError: " << err << " \n\n" << std::endl;
}
/* error.h
error handling definition
Project: Beam Envelope Tracker (BET)
Revision history
Date Description Programmer
------------ -------------------------------------------- --------------
07-03-06 Created Rene Bakker
Last Revision:
$Id: error.h 91 2007-05-06 17:47:45Z bakker $
*/
#ifndef _ERROR_DEF
#define _ERROR_DEF
#include <stdio.h>
#include <string>
#include <iostream>
void initBetErrorMsg(int level = 0, const char *fbase = NULL);
void initBetErrorFilename(const char *fbase);
void writeBetError(std::string err = "");
void setBetReportLevel(int level);
#endif
......@@ -53,11 +53,11 @@ static double rootValue = 0.0;
// used in setLShape for Gaussian
static void erfRoot(double x, double *fn, double *df) {
double v = erfc(fabs(x));
double v = std::erfc(std::abs(x));
double eps = 1.0e-05;
*fn = v - rootValue;
*df = (erfc(fabs(x) + eps) - v) / eps;
*df = (std::erfc(std::abs(x) + eps) - v) / eps;
}
......@@ -73,7 +73,7 @@ EnvelopeBunch::EnvelopeBunch(const PartData *ref):
}
EnvelopeBunch::EnvelopeBunch(const std::vector<OpalParticle> &rhs,
EnvelopeBunch::EnvelopeBunch(const std::vector<OpalParticle> &/*rhs*/,
const PartData *ref):
PartBunch(ref),
numSlices_m(0),
......@@ -370,7 +370,7 @@ void EnvelopeBunch::calcEmittance(double *emtnx, double *emtny, double *emtx, do
}
}
void EnvelopeBunch::calcEnergyChirp(double *g0, double *dgdt, double *gInc, int *nValid) {
void EnvelopeBunch::calcEnergyChirp(double *g0, double *dgdt, double *gInc, int */*nValid*/) {
std::vector<double> dtl(numMySlices_m, 0.0);
std::vector<double> b(numMySlices_m, 0.0);
std::vector<double> g(numMySlices_m, 0.0);
......@@ -448,12 +448,12 @@ void EnvelopeBunch::calcEnergyChirp(double *g0, double *dgdt, double *gInc, int
}
// chrip and uncorrelated energy sread
linfit(&dtG[0], &gG[0], nVTot, &gZero, &gt, &dum2, &dum3, &dum4);
linfit(&dtG[0], &gG[0], nVTot, gZero, gt, dum2, dum3, dum4);
*dgdt = gt;
rms = 0.0;
for(int i = 0; i < nVTot; i++) {
rms += pow(gG[i] - gZero - gt * dtG[i], 2);
rms += std::pow(gG[i] - gZero - gt * dtG[i], 2);
}
*gInc = sqrt(rms / nVTot);
......@@ -580,7 +580,7 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double
for(int i = 0; i < numMySlices_m; i++) {
slices_m[i]->p[SLI_z] = -(((numSlices_m - 1) - (mySliceStartOffset_m + i)) * bunch_width) / numSlices_m;
}
I0avg_m = Q_m * Physics::c / fabs(2.0 * emission_time_s);
I0avg_m = Q_m * Physics::c / std::abs(2.0 * emission_time_s);
break;
case bsGauss:
......@@ -589,7 +589,7 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double
for(int i = 1; i <= numSlices_m / 2; i++) {
rootValue = 1.0 - 2.0 * i * frac / (numSlices_m + 1);
v = fabs(emission_time_s) * sqr2 * findRoot(erfRoot, 0.0, 5.0, 1.0e-5) * (emission_time_s < 0.0 ? Physics::c : 1.0);
v = std::abs(emission_time_s) * sqr2 * findRoot(erfRoot, 0.0, 5.0, 1.0e-5) * (emission_time_s < 0.0 ? Physics::c : 1.0);
if((n2 + i) >= mySliceStartOffset_m && (n2 + i) <= mySliceEndOffset_m)
slices_m[n2+i-mySliceStartOffset_m]->p[SLI_z] = z0 + v * slices_m[n2+i-mySliceStartOffset_m]->p[SLI_beta];
......@@ -648,7 +648,7 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double
backup();
}
void EnvelopeBunch::setTShape(double enx, double eny, double rx, double ry, double b0) {
void EnvelopeBunch::setTShape(double enx, double eny, double rx, double ry, double /*b0*/) {
/*
Inform msg("setTshape");
msg << "set SLI_x to " << rx / 2.0 << endl;
......@@ -685,7 +685,7 @@ void EnvelopeBunch::setTOffset(double x0, double px0, double y0, double py0) {
void EnvelopeBunch::setEnergy(double E0, double dE) {
double g0 = 1.0 + (Physics::q_e * E0 / (Physics::EMASS * Physics::c * Physics::c));
double dg = fabs(dE) * Physics::q_e / (Physics::EMASS * Physics::c * Physics::c);
double dg = std::abs(dE) * Physics::q_e / (Physics::EMASS * Physics::c * Physics::c);
double z0 = zAvg();
for(int i = 0; i < numMySlices_m; i++) {
......@@ -794,7 +794,7 @@ void EnvelopeBunch::calcI() {
for(i = my_start; i <= vend; i++) {
j = 1;
do {
dz = fabs(z1[i+j] - z1[i-j]);
dz = std::abs(z1[i+j] - z1[i-j]);
j++;
} while((dz < dzMin * (j - 1)) && ((i + j) < n1) && ((i - j) >= 0));
......@@ -942,8 +942,8 @@ void EnvelopeBunch::cSpaceCharge() {
double R = 2 * sigma_av;
double A = R / L / getGamma(i);
double H1 = sqrt((1 - zeta / L) * (1 - zeta / L) + A * A) - sqrt((zeta / L) * (zeta / L) + A * A) - fabs(1 - zeta / L) + fabs(zeta / L);
double H2 = sqrt((1 - xi / L) * (1 - xi / L) + A * A) - sqrt((xi / L) * (xi / L) + A * A) - fabs(1 - xi / L) + fabs(xi / L);
double H1 = sqrt((1 - zeta / L) * (1 - zeta / L) + A * A) - sqrt((zeta / L) * (zeta / L) + A * A) - std::abs(1 - zeta / L) + std::abs(zeta / L);
double H2 = sqrt((1 - xi / L) * (1 - xi / L) + A * A) - sqrt((xi / L) * (xi / L) + A * A) - std::abs(1 - xi / L) + std::abs(xi / L);
//FIXME: Q_act or Q?
//double Q = activeSlices_m * Q_m / numSlices_m;
......@@ -1005,7 +1005,7 @@ void EnvelopeBunch::cSpaceCharge() {
for(int j = 0; j < numSlices_m; j++) {
double zj = z_m[j];
double dz = fabs(zj - z0);
double dz = std::abs(zj - z0);
if((dz > dzMin) && (zj > zCat)) {
double Aj = xi[j] / (dz * dz);
double v = 1.0 - (1.0 / sqrt(1.0 + Aj));
......@@ -1089,7 +1089,7 @@ ALT: SLI_z (commented by Rene)
// dYdt[SLI_z] = Y[SLI_beta]*c*cos(Y[SLI_px0]/Y[SLI_beta]/c)*cos(Y[SLI_py0]/Y[SLI_beta]/c);
*/
void EnvelopeBunch::derivs(double tc, double Y[], double dYdt[]) {
void EnvelopeBunch::derivs(double /*tc*/, double Y[], double dYdt[]) {
double g2 = 1.0 / (1.0 - Y[SLI_beta] * Y[SLI_beta]);
double g = sqrt(g2);
double g3 = g2 * g;
......@@ -1110,8 +1110,8 @@ void EnvelopeBunch::derivs(double tc, double Y[], double dYdt[]) {
if(solver_m & sv_radial) {
/// minimum spot-size due to emittance
/// \f[ \left(\frac{\epsilon_n c}{\gamma}\right)^2 \f]
double enxc2 = pow((emtnx0_m + emtbx0_m) * Physics::c / (Y[SLI_beta] * g), 2);
double enyc2 = pow((emtny0_m + emtby0_m) * Physics::c / (Y[SLI_beta] * g), 2);
double enxc2 = std::pow((emtnx0_m + emtbx0_m) * Physics::c / (Y[SLI_beta] * g), 2);
double enyc2 = std::pow((emtny0_m + emtby0_m) * Physics::c / (Y[SLI_beta] * g), 2);
#ifdef USE_HOMDYN_SC_MODEL
......@@ -1266,27 +1266,27 @@ void EnvelopeBunch::timeStep(double tStep, double _zCat) {
// set default allowed error for integration
double epsLocal = eps;
// mark that the integration was not successful yet
int ode_result = 1;
bool ode_result = false;
while(ode_result == 1) {
while(ode_result == false) {
if(solver_m & sv_fixedStep) {
rk4(&(sp->p[SLI_z]), SLNPAR, t_m, time_step_s, Gderivs);
ode_result = 0;
ode_result = true;
} else {
int nok, nbad;
activeSlice_m = i;
ode_result = odeint(&(sp->p[SLI_z]), SLNPAR, t_m, t_m + time_step_s, epsLocal, 0.1 * time_step_s, 0.0, &nok, &nbad, Gderivs);
ode_result = odeint(&(sp->p[SLI_z]), SLNPAR, t_m, t_m + time_step_s, epsLocal, 0.1 * time_step_s, 0.0, nok, nbad, Gderivs);
}
if(ode_result != 0) {
if(ode_result == false) {
// restore the backup
sp->restore();
epsLocal *= 10.0;
}
}
if(ode_result == 1) {
// :FIXME: the above loop cannot exit with ode_result == false
if(ode_result == false) {
// use fixed step integration if dynamic fails
rk4(&(sp->p[SLI_z]), SLNPAR, t_m, time_step_s, Gderivs);
......@@ -1351,7 +1351,7 @@ void EnvelopeBunch::timeStep(double tStep, double _zCat) {
} else {
msg << "EnvelopeBunch::run() No valid slices to subtract average" << endl;
}
beta = sqrt(1.0 - (1.0 / pow(ga, 2)));
beta = sqrt(1.0 - (1.0 / std::pow(ga, 2)));
fi_x = px0a / Physics::c / beta;
fi_y = py0a / Physics::c / beta;
......@@ -1371,7 +1371,10 @@ void EnvelopeBunch::timeStep(double tStep, double _zCat) {
}
}
void EnvelopeBunch::initialize(int num_slices, double Q, double energy, double width, double emission_time, double frac, double current, double bunch_center, double bX, double bY, double mX, double mY, double Bz0, int nbin) {
void EnvelopeBunch::initialize(int num_slices, double Q, double energy, double /*width*/,
double emission_time, double frac, double /*current*/,
double bunch_center, double bX, double bY, double mX,
double mY, double Bz0, int nbin) {
#ifdef USE_HOMDYN_SC_MODEL
*gmsg << "* Using HOMDYN space-charge model" << endl;
......@@ -1560,4 +1563,12 @@ Inform &EnvelopeBunch::slprint(Inform &os) {
//os.flags(ff);
}
return os;
}
\ No newline at end of file
}
// vi: set et ts=4 sw=4 sts=4:
// Local Variables:
// mode:c
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
// End:
\ No newline at end of file
/* fuctions.C
Collection of special functions
NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Revision history
Date Description Programmer
------------ -------------------------------------------- ------------
30/04/1992 Copied from "numerical recipies" Rene Bakker
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "Algorithms/bet/math/functions.h"
/* Bessel function
======================================================================= */
static double bessj0(double x) {
double ax, z;
double xx, y, ans, ans1, ans2;
if((ax = fabs(x)) < 8.0) {
y = x * x;
ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7
+ y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718
+ y * (59272.64853 + y * (267.8532712 + y * 1.0))));
ans = ans1 / ans2;
} else {
z = 8.0 / ax;
y = z * z;
xx = ax - 0.785398164;
ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
+ y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
+ y * (-0.6911147651e-5 + y * (0.7621095161e-6
- y * 0.934935152e-7)));
ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
}
return (double) ans;
} /* bessj0 */
static double bessj1(double x) {
double ax, z;
double xx, y, ans, ans1, ans2;
if((ax = fabs(x)) < 8.0) {
y = x * x;
ans1 = x * (72362614232.0 + y * (-7895059235.0 + y * (242396853.1
+ y * (-2972611.439 + y * (15704.48260 + y * (-30.16036606))))));
ans2 = 144725228442.0 + y * (2300535178.0 + y * (18583304.74
+ y * (99447.43394 + y * (376.9991397 + y * 1.0))));
ans = ans1 / ans2;
} else {
z = 8.0 / ax;
y = z * z;
xx = ax - 2.356194491;
ans1 = 1.0 + y * (0.183105e-2 + y * (-0.3516396496e-4
+ y * (0.2457520174e-5 + y * (-0.240337019e-6))));
ans2 = 0.04687499995 + y * (-0.2002690873e-3
+ y * (0.8449199096e-5 + y * (-0.88228987e-6
+ y * 0.105787412e-6)));
ans = sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
if(x < 0.0) ans = -ans;
}
return (double) ans;
} /* bessj1 */
/*************************\
* *
* External functions *
* *
\*************************/
#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10
double bessj(int n, double x) {
int j, jsum, m;
double ax, bj, bjm, bjp, sum, tox, ans;
if(n < 0) {
fprintf(stderr, "Illegal argument bessel function.\n");
exit(1);
}
if(n == 0) return bessj0(x);
if(n == 1) return bessj1(x);
ax = fabs(x);
if(ax == 0.0) return 0.0;
else if(ax > (double) n) {
tox = 2.0 / ax;
bjm = bessj0(ax);
bj = bessj1(ax);
for(j = 1; j < n; j++) {
bjp = j * tox * bj - bjm;
bjm = bj;
bj = bjp;
}
ans = bj;
} else {
tox = 2.0 / ax;
m = 2 * ((n + (int) sqrt(ACC * n)) / 2);
jsum = 0;
bjp = ans = sum = 0.0;
bj = 1.0;
for(j = m; j > 0; j--) {
bjm = j * tox * bj - bjp;
bjp = bj;
bj = bjm;
if(fabs(bj) > BIGNO) {
bj *= BIGNI;
bjp *= BIGNI;
ans *= BIGNI;
sum *= BIGNI;
}
if(jsum) sum += bj;
jsum = !jsum;
if(j == n) ans = bjp;
}
sum = 2.0 * sum - bj;
ans /= sum;
}
return x < 0.0 && n % 2 == 1 ? -ans : ans;
} /* bessj */
#undef ACC
#undef BIGNO
#undef BIGNI
/* Error function
======================================================================= */
/*
errfc()
Returns the complementary error function erfc(x) with fractional error
everywhere less than 1.2e-07. */
double errfc(double x) {
double t, z, ans;
z = fabs(x);
t = 1.0 / (1.0 + 0.5 * z);
ans = t * exp(-z * z - 1.26551223 +
t * (1.00002368 +
t * (0.37409196 +
t * (0.09678418 +
t * (-0.18628806 +
t * (0.27886807 +
t * (-1.13520398 +
t * (1.48851587 +
t * (-0.82215223 +
t * 0.17087277)))))))));
return x >= 0.0 ? ans : 2.0 - ans;
} /*erfc() */
/* functions.h
implementation of special functions
Project: Beam Envelope Tracker (BET)
Revision history
Date Description Programmer
------------ -------------------------------------------- --------------
09-03-06 Created Rene Bakker
Last Revision:
$Id: functions.h 29 2007-04-14 17:03:18Z l_bakker $
*/
#ifndef _FUNCTIONS_DEF
#define _FUNCTIONS_DEF
/* besselj()
Interger bessel function of the first kind
n = 0, 1, ......
*/
double bessj(int n, double x);
/* complementary error function */
double errfc(double x);
#endif
/* integration.C
integration routines using Rombergs method
NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Revision history
Date Description Programmer
------------ -------------------------------------------- ------------
10/03/2006 Copied from "numerical recipies" Rene Bakker
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "Algorithms/bet/BetError.h"
#include "Algorithms/bet/math/integrate.h"
/* internal functions
====================================================================== */
/* vector
allocate an array of doubles [0..n-1] and check memory */
static double *vector(int n) {
double *b;
b = (double *) malloc(sizeof(double) * n);
if(!b) {
writeBetError("Insufficient memory malloc bytes (rk.C)"); //,sizeof(double)*n);
}
return b;
} /* vector */
#define FUNC(x) ((*func)(x))
/* tranzd()
This routine computes the nth stage of refinement of an extended
trapezoidal rule. func is input as a pointer to the function to
be integrated between limits a and b, also input. When called with
n=1, the routine returns the crudest estimate of integral a->b f(x)dx.
Subsequent calls with n=2,3,... (in that sequential order) will improve
the accuracy by adding 2n-2 additional interior points.
*/
static double trapzd(double(*func)(double), double a, double b, int n) {
double x, tnm, sum, del;
static double s;
int it, j;
if(n == 1) {
return (s = 0.5 * (b - a) * (FUNC(a) + FUNC(b)));
} else {
for(it = 1, j = 1; j < n - 1; j++) it <<= 1;
tnm = it;
del = (b - a) / tnm;
x = a + 0.5 * del;
for(sum = 0.0, j = 1; j <= it; j++, x += del) sum += FUNC(x);
s = 0.5 * (s + (b - a) * sum / tnm);
return s;
}
} /* trapzd() */
#undef FUNC
/* polint
Given arrays xa[1..n] and ya[1..n], and given a value x,
this routine returns a value y, and an error estimate dy.
If P(x) is the polynomial of degree N - 1 such that
P(xa[i]) = ya[i], i = 1,..n, then the returned value y = P(x).
*/
static void polint
(
double xa[], double ya[], int n,
double x, double *y, double *dy) {
int i, m, ns = 1;
double den, dif, dift, ho, hp, w;
double *c, *d;
dif = fabs(x - xa[1]);
c = vector(n + 1);
d = vector(n + 1);
for(i = 1; i <= n; i++) {
if((dift = fabs(x - xa[i])) < dif) {
ns = i;
dif = dift;
}
c[i] = ya[i];
d[i] = ya[i];
}
*y = ya[ns--];
for(m = 1; m < n; m++) {
for(i = 1; i <= n - m; i++) {
ho = xa[i] - x;
hp = xa[i+m] - x;
w = c[i+1] - d[i];
if((den = ho - hp) == 0.0) {
int j;
fprintf(stderr, "Polint: n=%d, i=%d\n", n, i);