EnvelopeBunch.cpp 47.8 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6
#include <iostream>
#include <cfloat>
#include <fstream>
#include <cmath>
#include <string>
#include <limits>
7 8 9
#include <algorithm>
#include <typeinfo>

gsell's avatar
gsell committed
10 11 12 13 14 15 16 17 18 19 20
//FIXME: replace with IPPL Vector_t
#include <vector>
//FIXME: remove
#include <mpi.h>

#include "Algorithms/bet/math/root.h"     // root finding routines
#include "Algorithms/bet/math/linfit.h"   // linear fitting routines
#include "Algorithms/bet/math/savgol.h"   // savgol smoothing routine
#include "Algorithms/bet/math/rk.h"       // Runge-Kutta Integration

#include "Algorithms/bet/EnvelopeBunch.h"
21
#include "Utilities/Util.h"
22
#include "Utility/IpplTimings.h"
gsell's avatar
gsell committed
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

#define USE_HOMDYN_SC_MODEL

extern Inform *gmsg;

/** for space charge
    ignore for too low energies
    beta = 0.05 -> 640 eV
    beta = 0.10 ->   3 keV
    beta = 0.20 ->  11 keV
    beta = 0.30 ->  25 keV
    beta = 0.40 ->  47 keV
    beta = 0.50 ->  90 keV
    beta = 0.60 -> 128 keV
    beta = 0.70 -> 205 keV
    beta = 0.75 -> 261 keV
    beta = 0.80 -> 341 keV
    beta = 0.85 -> 460 keV
*/
#define BETA_MIN1 0.30     // minimum beta-value for space-charge calculations: start
snuverink_j's avatar
snuverink_j committed
43 44

#ifndef USE_HOMDYN_SC_MODEL
gsell's avatar
gsell committed
45
#define BETA_MIN2 0.45     // minimum beta-value for space-charge calculations: full impact
snuverink_j's avatar
snuverink_j committed
46
#endif
gsell's avatar
gsell committed
47 48 49 50 51 52 53 54

// Hack allows odeint in rk.C to be called with a class member function
static EnvelopeBunch *thisBunch = NULL;  // pointer to access calling bunch
static void Gderivs(double t, double Y[], double dYdt[]) { thisBunch->derivs(t, Y, dYdt); }
static double rootValue = 0.0;

// used in setLShape for Gaussian
static void erfRoot(double x, double *fn, double *df) {
gsell's avatar
gsell committed
55
    double v = std::erfc(std::abs(x));
gsell's avatar
gsell committed
56 57 58
    double eps = 1.0e-05;

    *fn = v - rootValue;
gsell's avatar
gsell committed
59
    *df = (std::erfc(std::abs(x) + eps) - v) / eps;
gsell's avatar
gsell committed
60 61
}

62 63 64 65 66 67 68 69 70 71 72 73 74

EnvelopeBunch::EnvelopeBunch(const PartData *ref):
    PartBunch(ref),
    numSlices_m(0),
    numMySlices_m(0) {

    calcITimer_m = IpplTimings::getTimer("calcI");
    spaceChargeTimer_m = IpplTimings::getTimer("spaceCharge");
    isValid_m = true;

}


gsell's avatar
gsell committed
75
EnvelopeBunch::EnvelopeBunch(const std::vector<OpalParticle> &/*rhs*/,
frey_m's avatar
frey_m committed
76
                             const PartData *ref):
77 78 79 80 81 82 83 84 85
    PartBunch(ref),
    numSlices_m(0),
    numMySlices_m(0)
{}


EnvelopeBunch::~EnvelopeBunch() {
}

gsell's avatar
gsell committed
86 87 88 89
void EnvelopeBunch::calcBeamParameters() {
    Inform msg("calcBeamParameters");
    IpplTimings::startTimer(statParamTimer_m);
    double ex, ey, nex, ney, b0, bRms, bMax, bMin, g0, dgdt, gInc;
90
    double RxRms = 0.0, RyRms = 0.0, Px = 0.0, PxRms = 0.0, PxMax = 0.0, PxMin = 0.0, Py = 0.0, PyRms = 0.0, PyMax = 0.0, PyMin = 0.0;
gsell's avatar
gsell committed
91 92 93 94 95 96 97 98
    double Pz, PzMax, PzMin, PzRms;
    double x0Rms, y0Rms, zRms, zMax, zMin, I0, IRms, IMax, IMin;

    runStats(sp_beta, &b0, &bMax, &bMin, &bRms, &nValid_m);
    runStats(sp_I, &I0, &IMax, &IMin, &IRms, &nValid_m);
    runStats(sp_z, &z0_m, &zMax, &zMin, &zRms, &nValid_m);
    runStats(sp_Pz, &Pz, &PzMax, &PzMin, &PzRms, &nValid_m);

99
    if(solver_m & sv_radial) {
gsell's avatar
gsell committed
100 101 102 103 104 105 106
        runStats(sp_Rx, &Rx_m, &RxMax_m, &RxMin_m, &RxRms, &nValid_m);
        runStats(sp_Ry, &Ry_m, &RyMax_m, &RyMin_m, &RyRms, &nValid_m);
        runStats(sp_Px, &Px, &PxMax, &PxMin, &PxRms, &nValid_m);
        runStats(sp_Py, &Py, &PyMax, &PyMin, &PyRms, &nValid_m);
        calcEmittance(&nex, &ney, &ex, &ey, &nValid_m);
    }

107
    if(solver_m & sv_offaxis) {
gsell's avatar
gsell committed
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
        runStats(sp_x0, &x0_m, &x0Max_m, &x0Min_m, &x0Rms, &nValid_m);
        runStats(sp_y0, &y0_m, &y0Max_m, &y0Min_m, &y0Rms, &nValid_m);
    }

    calcEnergyChirp(&g0, &dgdt, &gInc, &nValid_m);
    double Bfz = AvBField();
    double Efz = AvEField();
    double mc2e = 1.0e-6 * Physics::EMASS * Physics::c * Physics::c / Physics::q_e;

    E_m = mc2e * (g0 - 1.0);
    dEdt_m = 1.0e-12 * mc2e * dgdt;
    Einc_m = mc2e * gInc;
    tau_m =  zRms / Physics::c;
    I_m = IMax;
    Irms_m = Q_m * nValid_m * Physics::c / (zRms * sqrt(Physics::two_pi) * numSlices_m);
    Px_m = Px / Physics::c;
    Py_m = Py / Physics::c;
    Ez_m = 1.0e-6 * Efz;
    Bz_m = Bfz;

    //in [mrad]
    emtn_m = Vector_t(ex, ey, 0.0);
    norm_emtn_m = Vector_t(nex, ney, 0.0);

    maxX_m = Vector_t(RxMax_m, RyMax_m, zMax);
    maxP_m = Vector_t(PxMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, PyMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, PzMax);

    //minX in the T-T sense is -RxMax_m
    minX_m = Vector_t(-RxMax_m, -RyMax_m, zMin);
    minP_m = Vector_t(-PxMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, -PyMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, PzMin);

    /* PxRms is the rms of the divergence in x direction in [rad]: x'
     * x' = Px/P0 -> Px = x' * P0 = x' * \beta/c*E (E is the particle total
     * energy)
     * Pz* in [\beta\gamma]
     * \pi from [rad]
     */
    sigmax_m = Vector_t(RxRms / 2.0, RyRms / 2.0, zRms);
    sigmap_m = Vector_t(PxRms * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi / 2.0, PyRms * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi / 2.0, PzRms);

    dEdt_m = PzRms * mc2e * b0;

    IpplTimings::stopTimer(statParamTimer_m);
}

void EnvelopeBunch::runStats(EnvelopeBunchParameter sp, double *xAvg, double *xMax, double *xMin, double *rms, int *nValid) {
    int i, nV = 0;
155
    std::vector<double> v(numMySlices_m, 0.0);
gsell's avatar
gsell committed
156 157 158 159 160

    //FIXME: why from 1 to n-1??
    switch(sp) {
        case sp_beta:      // normalized velocity (total) [-]
            for(i = 1; i < numMySlices_m - 1; i++) {
161
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_beta];
gsell's avatar
gsell committed
162 163 164 165
            }
            break;
        case sp_gamma:     // Lorenz factor
            for(i = 1; i < numMySlices_m - 1; i++) {
166
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->computeGamma();
gsell's avatar
gsell committed
167 168 169 170
            }
            break;
        case sp_z:         // slice position [m]
            for(i = 0; i < numMySlices_m - 0; i++) {
171
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_z];
gsell's avatar
gsell committed
172 173 174 175
            }
            break;
        case sp_I:         // slice position [m]
            for(i = 1; i < numMySlices_m - 1; i++) {
176 177
                if((slices_m[i]->p[SLI_beta] > BETA_MIN1) && ((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()))
                    v[nV++] = 0.0; //FIXME: currentProfile_m->get(slices_m[i]->p[SLI_z], itype_lin);
gsell's avatar
gsell committed
178 179 180 181
            }
            break;
        case sp_Rx:        // beam size x [m]
            for(i = 0; i < numMySlices_m - 0; i++) {
182
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = 2.* slices_m[i]->p[SLI_x];
gsell's avatar
gsell committed
183 184 185 186
            }
            break;
        case sp_Ry:        // beam size y [m]
            for(i = 0; i < numMySlices_m - 0; i++) {
187
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = 2. * slices_m[i]->p[SLI_y];
gsell's avatar
gsell committed
188 189 190 191
            }
            break;
        case sp_Px:        // beam divergence x
            for(i = 0; i < numMySlices_m - 0; i++) {
192
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_px];
gsell's avatar
gsell committed
193 194 195 196
            }
            break;
        case sp_Py:        // beam divergence y
            for(i = 0; i < numMySlices_m - 0; i++) {
197
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_py];
gsell's avatar
gsell committed
198 199 200 201
            }
            break;
        case sp_Pz:
            for(i = 1; i < numMySlices_m - 1; i++) {
202
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_beta] * slices_m[i]->computeGamma();
gsell's avatar
gsell committed
203 204 205 206
            }
            break;
        case sp_x0:        // position centroid x [m]
            for(i = 1; i < numMySlices_m - 1; i++) {
207
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_x0];
gsell's avatar
gsell committed
208 209 210 211
            }
            break;
        case sp_y0:        // position centroid y [m]
            for(i = 1; i < numMySlices_m - 1; i++) {
212
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_y0];
gsell's avatar
gsell committed
213 214 215 216
            }
            break;
        case sp_px0:       // angular deflection centroid x
            for(i = 1; i < numMySlices_m - 1; i++) {
217
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_px0];
gsell's avatar
gsell committed
218 219 220 221
            }
            break;
        case sp_py0:      // angular deflection centroid y
            for(i = 1; i < numMySlices_m - 1; i++) {
222
                if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) v[nV++] = slices_m[i]->p[SLI_py0];
gsell's avatar
gsell committed
223 224 225 226 227 228 229 230
            }
            break;
        default :
            throw OpalException("EnvelopeBunch", "EnvelopeBunch::runStats() Undefined label (programming error)");
            break;
    }

    int nVTot = nV;
231
    allreduce(&nVTot, 1, std::plus<int>());
gsell's avatar
gsell committed
232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
    if(nVTot <= 0) {
        *xAvg = 0.0;
        *xMax = 0.0;
        *xMin = 0.0;
        *rms = 0.0;
        *nValid = 0;
    } else {
        double M1 = 0.0;
        double M2 = 0.0;
        double maxv = std::numeric_limits<double>::min();
        double minv = std::numeric_limits<double>::max();
        if(nV > 0) {
            M1 = v[0];
            M2 = v[0] * v[0];
            maxv = v[0];
            minv = v[0];
            for(i = 1; i < nV; i++) {
                M1 += v[i];
                M2 += v[i] * v[i];
                maxv = std::max(maxv, v[i]);
                minv = std::min(minv, v[i]);
            }
        }

256 257 258 259
        allreduce(&M1, 1, std::plus<double>());
        allreduce(&M2, 1, std::plus<double>());
        allreduce(&maxv, 1, std::greater<double>());
        allreduce(&minv, 1, std::less<double>());
gsell's avatar
gsell committed
260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299

        *xAvg = M1 / nVTot;
        *xMax = maxv;
        *xMin = minv;

        //XXX: ok so this causes problems. E.g in the case of all transversal
        //components we cannot compare a rms_rx with rms_x (particles)
        //produced by the T-Tracker

        //in case of transversal stats we want to calculate the rms
        //*rms = sqrt(M2/nVTot);

        //else a sigma
        //*sigma = sqrt(M2/nVTot - M1*M1/(nVTot*nVTot));
        //this is sigma:
        //*rms = sqrt(M2/nVTot - M1*M1/(nVTot*nVTot));
        *nValid = nVTot;

        if(sp == sp_Rx ||
           sp == sp_Ry ||
           sp == sp_Px ||
           sp == sp_Py)
            *rms = sqrt(M2 / nVTot); //2.0: half of the particles
        else
            *rms = sqrt(M2 / nVTot - M1 * M1 / (nVTot * nVTot));
    }
}

void EnvelopeBunch::calcEmittance(double *emtnx, double *emtny, double *emtx, double *emty, int *nValid) {
    double sx = 0.0;
    double sxp = 0.0;
    double sxxp = 0.0;
    double sy = 0.0;
    double syp = 0.0;
    double syyp = 0.0;
    double betagamma = 0.0;

    // find the amount of active slices
    int nV = 0;
    for(int i = 0; i < numMySlices_m; i++) {
300
        if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid())
gsell's avatar
gsell committed
301 302 303 304 305 306 307 308
            nV++;
    }

    if(nV > 0) {
        int i1 = nV;
        nV = 0;
        double bg = 0.0;
        for(int i = 0; i < i1; i++) {
309
            if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid()) {
gsell's avatar
gsell committed
310 311
                nV++;

312
                if(solver_m & sv_radial) {
gsell's avatar
gsell committed
313
                    assert(i < numMySlices_m);
314
                    auto sp = slices_m[i];
315
                    bg = sp->p[SLI_beta] * sp->computeGamma();
gsell's avatar
gsell committed
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355

                    double pbc = bg * sp->p[SLI_px] / (sp->p[SLI_beta] * Physics::c);
                    sx   += sp->p[SLI_x] * sp->p[SLI_x];
                    sxp  += pbc * pbc;
                    sxxp += sp->p[SLI_x] * pbc;

                    pbc   = bg * sp->p[SLI_py] / (sp->p[SLI_beta] * Physics::c);
                    sy   += sp->p[SLI_y] * sp->p[SLI_y];
                    syp  += pbc * pbc;
                    syyp += sp->p[SLI_y] * pbc;

                    betagamma += sqrt(1 + sp->p[SLI_px] * sp->p[SLI_px] + sp->p[SLI_py] * sp->p[SLI_py]);
                }
            }
        }
    }

    int nVToT = nV;
    reduce(nVToT, nVToT, OpAddAssign());
    if(nVToT == 0) {
        *emtnx = 0.0;
        *emtny = 0.0;
        *emtx = 0.0;
        *emty = 0.0;
        *nValid = 0;
    } else {
        reduce(sx, sx, OpAddAssign());
        reduce(sy, sy, OpAddAssign());
        reduce(sxp, sxp, OpAddAssign());
        reduce(syp, syp, OpAddAssign());
        reduce(sxxp, sxxp, OpAddAssign());
        reduce(syyp, syyp, OpAddAssign());
        reduce(betagamma, betagamma, OpAddAssign());
        sx /= nVToT;
        sy /= nVToT;
        sxp /= nVToT;
        syp /= nVToT;
        sxxp /= nVToT;
        syyp /= nVToT;

356 357
        *emtnx = sqrt(sx * sxp - sxxp * sxxp + emtnx0_m * emtnx0_m + emtbx0_m * emtbx0_m);
        *emtny = sqrt(sy * syp - syyp * syyp + emtny0_m * emtny0_m + emtby0_m * emtby0_m);
gsell's avatar
gsell committed
358 359 360 361 362 363 364 365 366

        betagamma /= nVToT;
        betagamma *= sqrt(1.0 - (1 / betagamma) * (1 / betagamma));
        *emtx = *emtnx / betagamma;
        *emty = *emtny / betagamma;
        *nValid = nVToT;
    }
}

gsell's avatar
gsell committed
367
void EnvelopeBunch::calcEnergyChirp(double *g0, double *dgdt, double *gInc, int */*nValid*/) {
368 369 370
    std::vector<double> dtl(numMySlices_m, 0.0);
    std::vector<double> b(numMySlices_m, 0.0);
    std::vector<double> g(numMySlices_m, 0.0);
gsell's avatar
gsell committed
371 372 373 374 375 376 377 378 379 380

    // defaults
    *g0   = 1.0;
    *dgdt = 0.0;
    *gInc = 0.0;

    double zAvg = 0.0;
    double gAvg = 0.0;
    int j = 0;
    for(int i = 0; i < numMySlices_m; i++) {
381 382
        std::shared_ptr<EnvelopeSlice> cs = slices_m[i];
        if((cs->is_valid()) && (cs->p[SLI_z] > zCat_m)) {
gsell's avatar
gsell committed
383 384 385
            zAvg    += cs->p[SLI_z];
            dtl[i-j]  = cs->p[SLI_z];
            b[i-j]   = cs->p[SLI_beta];
386
            g[i-j]   = cs->computeGamma();
gsell's avatar
gsell committed
387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
            gAvg    += g[i-j];
        } else
            ++j;
    }
    int nV = numMySlices_m - j;

    int nVTot = nV;
    reduce(nVTot, nVTot, OpAddAssign());
    if(nVTot > 0) {
        reduce(gAvg, gAvg, OpAddAssign());
        reduce(zAvg, zAvg, OpAddAssign());
        gAvg = gAvg / nVTot;
        zAvg = zAvg / nVTot;
        *g0  = gAvg;
    }

    // FIXME: working with global arrays
    // make dtl, b and g global
    if(nVTot > 2) {
406 407 408
        std::vector<double> dtG(nVTot, 0.0);
        std::vector<double> bG(nVTot, 0.0);
        std::vector<double> gG(nVTot, 0.0);
gsell's avatar
gsell committed
409 410

        int numproc = Ippl::Comm->getNodes();
411 412 413 414
        std::vector<int> numsend(numproc, nV);
        std::vector<int> offsets(numproc);
        std::vector<int> offsetsG(numproc);
        std::vector<int> zeros(numproc, 0);
gsell's avatar
gsell committed
415

416
        MPI_Allgather(&nV, 1, MPI_INT, &offsets.front(), 1, MPI_INT, Ippl::getComm());
gsell's avatar
gsell committed
417 418 419 420 421 422 423 424
        offsetsG[0] = 0;
        for(int i = 1; i < numproc; i++) {
            if(offsets[i-1] == 0)
                offsetsG[i] = 0;
            else
                offsetsG[i] = offsetsG[i-1] + offsets[i-1];
        }

425 426 427 428 429 430 431
        MPI_Alltoallv(&dtl.front(), &numsend.front(), &zeros.front(), MPI_DOUBLE,
                      &dtG.front(), &offsets.front(), &offsetsG.front(), MPI_DOUBLE,
                      Ippl::getComm());

        MPI_Alltoallv(&b.front(), &numsend.front(), &zeros.front(), MPI_DOUBLE,
                      &bG.front(), &offsets.front(), &offsetsG.front(), MPI_DOUBLE,
                      Ippl::getComm());
gsell's avatar
gsell committed
432

433 434 435
        MPI_Alltoallv(&g.front(), &numsend.front(), &zeros.front(), MPI_DOUBLE,
                      &gG.front(), &offsets.front(), &offsetsG.front(), MPI_DOUBLE,
                      Ippl::getComm());
gsell's avatar
gsell committed
436 437 438 439 440 441 442 443 444

        double dum2, dum3, dum4, rms, gZero, gt;

        // convert z to t
        for(int i = 0; i < nVTot; i++) {
            dtG[i] = (dtG[i] - zAvg) / (bG[i] * Physics::c);
        }

        // chrip and uncorrelated energy sread
gsell's avatar
gsell committed
445
        linfit(&dtG[0], &gG[0], nVTot, gZero, gt, dum2, dum3, dum4);
gsell's avatar
gsell committed
446 447 448 449
        *dgdt = gt;

        rms = 0.0;
        for(int i = 0; i < nVTot; i++) {
450
            rms += std::pow(gG[i] - gZero - gt * dtG[i], 2);
gsell's avatar
gsell committed
451 452 453 454 455 456 457
        }

        *gInc = sqrt(rms / nVTot);
    }
}


458
void EnvelopeBunch::distributeSlices(int nSlice) {
gsell's avatar
gsell committed
459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484
    numSlices_m = nSlice;
    int rank = Ippl::Comm->myNode();
    int numproc = Ippl::Comm->getNodes();
    numMySlices_m = nSlice / numproc;
    if(numMySlices_m < 13) {
        if(rank == 0) {
            numMySlices_m = 14;
        } else {
            numMySlices_m = (nSlice - 14) / (numproc - 1);
            if(rank - 1 < (nSlice - 14) % (numproc - 1))
                numMySlices_m++;
        }
    } else {
        if(rank < nSlice % numproc)
            numMySlices_m++;
    }

    mySliceStartOffset_m = rank * ((int)numSlices_m / numproc);
    if(rank < numSlices_m % numproc)
        mySliceStartOffset_m += rank;
    else
        mySliceStartOffset_m += numSlices_m % numproc;
    mySliceEndOffset_m = mySliceStartOffset_m + numMySlices_m - 1;
}


485 486 487
void EnvelopeBunch::createBunch() {

    // Memory issue when going to larger number of processors (this is
gsell's avatar
gsell committed
488 489
    // also the reason for the 30 slices hard limit per core!)
    // using heuristic (until problem is located)
490 491 492 493 494 495 496 497 498 499 500
    //size_t nSlices = (3 / 100.0 + 1.0) * numMySlices_m;

    size_t nSlices = getLocalNum();

    KR = std::unique_ptr<Vector_t[]>(new Vector_t[nSlices]);
    KT = std::unique_ptr<Vector_t[]>(new Vector_t[nSlices]);
    EF = std::unique_ptr<Vector_t[]>(new Vector_t[nSlices]);
    BF = std::unique_ptr<Vector_t[]>(new Vector_t[nSlices]);

    z_m.resize(numSlices_m, 0.0);
    b_m.resize(numSlices_m, 0.0);
gsell's avatar
gsell committed
501 502 503 504 505 506 507 508 509

    for(unsigned int i = 0; i < nSlices; i++) {
        KR[i] = Vector_t(0);
        KT[i] = Vector_t(0);
        BF[i] = Vector_t(0);
        EF[i] = Vector_t(0);
    }

    // set default DE solver method
510
    solver_m = sv_radial | sv_offaxis | sv_lwakes | sv_twakes;
gsell's avatar
gsell committed
511

512
    //FIXME: WHY?
gsell's avatar
gsell committed
513 514 515 516
    if(numSlices_m < 14)
        throw OpalException("EnvelopeBunch::createSlices", "use more than 13 slices");

    // defaults:
517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
    Q_m        = 0.0;     // no charge
    t_m        = 0.0;     // t = 0 s
    t_offset_m = 0.0;     // offset time by tReset function
    emtnx0_m   = 0.0;
    emtny0_m   = 0.0;     // no intrinsic emittance
    emtbx0_m   = 0.0;
    emtby0_m   = 0.0;     // no intrinsic emittance Bush effect
    Bz0_m      = 0.0;     // no magnetic field on creation (cathode)
    dx0_m      = 0.0;
    dy0_m      = 0.0;     // no offset of coordinate system
    dfi_x_m    = 0.0;
    dfi_y_m    = 0.0;     // no rotation of coordinate system

    slices_m.resize(nSlices);
    for( auto & slice : slices_m ) {
532 533 534 535 536 537 538 539 540
        slice = std::shared_ptr<EnvelopeSlice>(new EnvelopeSlice());
    }
    //XXX: not supported by clang at the moment
    //std::generate(s.begin(), s.end(),
        //[]()
        //{
            //return std::shared_ptr<EnvelopeSlice>(new EnvelopeSlice());
        //});

541 542 543 544 545
    Esct_m.resize(nSlices, 0.0);
    G_m.resize(nSlices, 0.0);
    Exw_m.resize(nSlices, 0.0);
    Eyw_m.resize(nSlices, 0.0);
    Ezw_m.resize(nSlices, 0.0);
546

gsell's avatar
gsell committed
547
    for(unsigned int i = 0; i < nSlices; i++) {
548 549 550 551 552
        Esct_m[i] = 0.0;
        G_m[i]    = 0.0;
        Ezw_m[i]  = 0.0;
        Exw_m[i]  = 0.0;
        Eyw_m[i]  = 0.0;
gsell's avatar
gsell committed
553 554
    }

555
    currentProfile_m = NULL;
556 557
    I0avg_m = 0.0;
    dStat_m = ds_fieldsSynchronized | ds_slicesSynchronized;
gsell's avatar
gsell committed
558 559 560 561 562 563 564 565 566 567 568 569 570 571
}

//XXX: convention: s[n] is the slice closest to the cathode and all positions
//are negative (slices left of cathode)
//XXX: every processor fills its slices in bins (every proc holds all bins,
//some empty)
void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double emission_time_s, double frac) {
    unsigned int n2 = numSlices_m / 2;
    double sqr2 = sqrt(2.0), v;
    double bunch_width = 0.0;

    switch(shape) {
        case bsRect:

572
            bunch_width = Physics::c * emission_time_s * slices_m[0]->p[SLI_beta];
gsell's avatar
gsell committed
573
            for(int i = 0; i < numMySlices_m; i++) {
574
                slices_m[i]->p[SLI_z] = -(((numSlices_m - 1) - (mySliceStartOffset_m + i)) * bunch_width) / numSlices_m;
gsell's avatar
gsell committed
575
            }
576
            I0avg_m = Q_m * Physics::c / std::abs(2.0 * emission_time_s);
gsell's avatar
gsell committed
577 578 579 580
            break;

        case bsGauss:
            if(n2 >= mySliceStartOffset_m && n2 <= mySliceEndOffset_m)
581
                slices_m[n2-mySliceStartOffset_m]->p[SLI_z] = z0;
gsell's avatar
gsell committed
582 583 584

            for(int i = 1; i <= numSlices_m / 2; i++) {
                rootValue = 1.0 - 2.0 * i * frac / (numSlices_m + 1);
585
                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);
gsell's avatar
gsell committed
586 587

                if((n2 + i) >= mySliceStartOffset_m && (n2 + i) <= mySliceEndOffset_m)
588
                    slices_m[n2+i-mySliceStartOffset_m]->p[SLI_z] = z0 + v * slices_m[n2+i-mySliceStartOffset_m]->p[SLI_beta];
gsell's avatar
gsell committed
589 590

                if((n2 - i) >= mySliceStartOffset_m && (n2 - i) <= mySliceEndOffset_m)
591
                    slices_m[n2-i-mySliceStartOffset_m]->p[SLI_z] = z0 - v * slices_m[n2-i-mySliceStartOffset_m]->p[SLI_beta];
gsell's avatar
gsell committed
592 593
            }

594
            I0avg_m = 0.0;
gsell's avatar
gsell committed
595 596 597 598 599
            break;
    }

    double gz0 = 0.0, gzN = 0.0;
    if(Ippl::Comm->myNode() == 0)
600
        gz0 = slices_m[0]->p[SLI_z];
gsell's avatar
gsell committed
601
    if(Ippl::Comm->myNode() == Ippl::Comm->getNodes() - 1)
602
        gzN = slices_m[numMySlices_m-1]->p[SLI_z];
gsell's avatar
gsell committed
603 604 605
    reduce(gz0, gz0, OpAddAssign());
    reduce(gzN, gzN, OpAddAssign());

606
    hbin_m = (gzN - gz0) / nebin_m;
gsell's avatar
gsell committed
607 608 609 610 611 612 613 614 615 616

    // initialize all bins with an empty vector
    for(int i = 0; i < nebin_m; i++) {
        std::vector<int> tmp;
        this->bins_m.push_back(tmp);
    }

    // add all slices to appropriated bins
    int bin_i = 0, slice_i = 0;
    while(slice_i < numMySlices_m) {
617
        if((bin_i + 1) * hbin_m < slices_m[slice_i]->p[SLI_z] - gz0) {
gsell's avatar
gsell committed
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632
            bin_i++;
        } else {
            this->bins_m[bin_i].push_back(slice_i);
            slice_i++;
        }
    }

    //find first bin containing slices
    firstBinWithValue_m = 0;
    for(; firstBinWithValue_m < nebin_m; firstBinWithValue_m++)
        if(bins_m[firstBinWithValue_m].size() != 0) break;

    backup();
}

gsell's avatar
gsell committed
633
void EnvelopeBunch::setTShape(double enx, double eny, double rx, double ry, double /*b0*/) {
gsell's avatar
gsell committed
634
    // set emittances
635 636
    emtnx0_m = enx;
    emtny0_m = eny;
gsell's avatar
gsell committed
637
    //FIXME: is rx here correct?
638 639
    emtbx0_m = Physics::q_e * rx * rx * Bz0_m / (8.0 * Physics::EMASS * Physics::c);
    emtby0_m = Physics::q_e * ry * ry * Bz0_m / (8.0 * Physics::EMASS * Physics::c);
gsell's avatar
gsell committed
640 641 642

    //XXX: rx = 2*rms_x
    for(int i = 0; i < numMySlices_m; i++) {
643 644 645 646
        slices_m[i]->p[SLI_x] = rx / 2.0; //rms_x
        slices_m[i]->p[SLI_y] = ry / 2.0; //rms_y
        slices_m[i]->p[SLI_px] = 0.0;
        slices_m[i]->p[SLI_py] = 0.0;
gsell's avatar
gsell committed
647 648 649 650 651 652 653
    }

    backup();
}

void EnvelopeBunch::setTOffset(double x0, double px0, double y0, double py0) {
    for(int i = 0; i < numMySlices_m; i++) {
654 655 656 657
        slices_m[i]->p[SLI_x0] = x0;
        slices_m[i]->p[SLI_px0] = px0;
        slices_m[i]->p[SLI_y0] = y0;
        slices_m[i]->p[SLI_py0] = py0;
gsell's avatar
gsell committed
658 659 660 661 662
    }
}

void EnvelopeBunch::setEnergy(double E0, double dE) {
    double g0    = 1.0 + (Physics::q_e * E0 / (Physics::EMASS * Physics::c * Physics::c));
663
    double dg    = std::abs(dE) * Physics::q_e / (Physics::EMASS * Physics::c * Physics::c);
gsell's avatar
gsell committed
664 665 666
    double z0    = zAvg();

    for(int i = 0; i < numMySlices_m; i++) {
667 668
        double g = g0 + (slices_m[i]->p[SLI_z] - z0) * dg;
        slices_m[i]->p[SLI_beta] = sqrt(1.0 - (1.0 / (g * g)));
gsell's avatar
gsell committed
669 670 671 672 673 674 675 676 677 678 679 680
    }

    backup();
}

void EnvelopeBunch::synchronizeSlices() {

    for(int i = 0; i < numSlices_m; i++) {
        z_m[i] = 0.0;
        b_m[i] = 0.0;
    }
    for(int i = 0; i < numMySlices_m; i++) {
681 682
        b_m[mySliceStartOffset_m+i] = slices_m[i]->p[SLI_beta];
        z_m[mySliceStartOffset_m+i] = slices_m[i]->p[SLI_z];
gsell's avatar
gsell committed
683 684
    }

frey_m's avatar
frey_m committed
685 686
    allreduce(z_m.data(), numSlices_m, std::plus<double>());
    allreduce(b_m.data(), numSlices_m, std::plus<double>());
gsell's avatar
gsell committed
687 688 689 690 691 692
}

void EnvelopeBunch::calcI() {
    Inform msg("calcI ");

    static int already_called = 0;
693
    if((dStat_m & ds_currentCalculated) || (already_called && (Q_m <= 0.0)))
gsell's avatar
gsell committed
694 695 696
        return;
    already_called = 1;

697
    std::vector<double> z1(numSlices_m, 0.0);
698
    std::vector<double> b (numSlices_m, 0.0);
gsell's avatar
gsell committed
699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723
    double bSum = 0.0;
    double dz2Sum = 0.0;
    int n1 = 0;
    int my_start = 0, my_end = 0;

    for(int i = 0; i < numSlices_m; i++) {
        if(b_m[i] > 0.0) {
            if((unsigned int) i == mySliceStartOffset_m)
                my_start = n1;
            if((unsigned int) i == mySliceEndOffset_m)
                my_end = n1;

            b[n1] = b_m[i];
            z1[n1] = z_m[i];
            if(n1 > 0)
                dz2Sum += ((z1[n1] - z1[n1-1]) * (z1[n1] - z1[n1-1]));
            bSum += b_m[i];
            n1++;
        }
    }
    if(Ippl::Comm->myNode() == 0)
        my_start++;

    if(n1 < 2) {
        msg << "n1 (= " << n1 << ") < 2" << endl;
724
        currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
gsell's avatar
gsell committed
725 726 727 728 729 730
        return;
    }

    double sigma_dz = sqrt(dz2Sum / (n1 - 1));
    double beta = bSum / n1;

731 732 733 734 735 736 737 738 739 740 741 742 743 744
    //sort beta's according to z1 with temporary vector of pairs
    std::vector<std::pair<double,double>> z1_b(numSlices_m);
    for (int i=0; i<numSlices_m; i++) {
        z1_b[i] = std::make_pair(z1[i],b[i]);
    }

    std::sort(z1_b.begin(), z1_b.end(),
              [](const std::pair<double,double> &left, const std::pair<double,double> &right)
              {return left.first < right.first;});

    for (int i=0; i<numSlices_m; i++) {
        z1[i] = z1_b[i].first;
        b [i] = z1_b[i].second;
    }
gsell's avatar
gsell committed
745 746 747 748 749

    double q = Q_m > 0.0 ? Q_m / numSlices_m : Physics::q_e;
    double dz = 0.0;

    // 1. Determine current from the slice distance
750
    std::vector<double> I1(n1, 0.0);
gsell's avatar
gsell committed
751 752 753 754 755 756 757 758 759 760 761 762 763

    //limit the max current to 5x the sigma value to reduce noise problems
    double dzMin = 0.2 * sigma_dz;

    //FIXME: Instead of using such a complicated mechanism (store same value
    //again) to enforce only to use a subsample of all points, use a vector to
    //which all values get pushed, no need for elimination later.

    int vend = my_end;
    if(Ippl::Comm->myNode() == Ippl::Comm->getNodes() - 1)
        vend--;

    //XXX: THIS LOOP IS THE EXPENSIVE PART!!!
gsell's avatar
gsell committed
764
    int i, j;
gsell's avatar
gsell committed
765 766 767
#ifdef _OPENMP
#pragma omp parallel for private(i,j,dz)
#endif //_OPENMP
gsell's avatar
gsell committed
768 769
    for(i = my_start; i <= vend; i++) {
        j = 1;
gsell's avatar
gsell committed
770
        do {
771
            dz = std::abs(z1[i+j] - z1[i-j]);
gsell's avatar
gsell committed
772 773 774 775 776 777 778 779 780 781
            j++;
        } while((dz < dzMin * (j - 1)) && ((i + j) < n1) && ((i - j) >= 0));

        if((dz >= dzMin * (j - 1)) && ((i + j) < n1) && ((i - j) >= 0)) {
            I1[i] = 0.25 * q * Physics::c * (b[i+j] + b[i-j]) / (dz * (j - 1));   // WHY 0.25?
        } else {
            I1[i] = 0.0;
        }
    }

frey_m's avatar
frey_m committed
782
    allreduce(I1.data(), n1, std::plus<double>());
gsell's avatar
gsell committed
783 784 785 786 787 788 789 790 791 792 793 794
    for(int i = 1; i < n1 - 1; i++) {
        if(I1[i] == 0.0)
            I1[i] = I1[i-1];
    }
    I1[0]    = I1[1];
    I1[n1-1] = I1[n1-2];


    //2. Remove points with identical z-value and then smooth the current profile
    double zMin = zTail();
    double zMax = zHead();
    dz = (zMax - zMin) / numSlices_m; // create a window of the average slice distance
795 796
    std::vector<double> z2(n1, 0.0);
    std::vector<double> I2(n1, 0.0);
gsell's avatar
gsell committed
797 798 799
    double Mz1 = 0.0;
    double MI1 = 0.0;
    int np = 0;
gsell's avatar
gsell committed
800
    j = 0;
gsell's avatar
gsell committed
801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849

    // XXX: COMPUTE ON ALL NODES
    // first value
    while((j < n1) && ((z1[j] - z1[0]) <= dz)) {
        Mz1 += z1[j];
        MI1 += I1[j];
        ++j;
        ++np;
    }
    z2[0] = Mz1 / np;
    I2[0] = MI1 / np;

    // following values
    int k = 0;
    for(int i = 1; i < n1; i++) {
        // add new points
        int j = 0;
        while(((i + j) < n1) && ((z1[i+j] - z1[i]) <= dz)) {
            if((z1[i+j] - z1[i-1]) > dz) {
                Mz1 += z1[i+j];
                MI1 += I1[i+j];
                ++np;
            }
            ++j;
        }

        // remove obsolete points
        j = 1;
        while(((i - j) >= 0) && ((z1[i-1] - z1[i-j]) < dz)) {
            if((z1[i] - z1[i-j]) > dz) {
                Mz1 -= z1[i-j];
                MI1 -= I1[i-j];
                --np;
            }
            ++j;
        }
        z2[i-k] = Mz1 / np;
        I2[i-k] = MI1 / np;

        // make sure there are no duplicate z coordinates
        if(z2[i-k] <= z2[i-k-1]) {
            I2[i-k-1] = 0.5 * (I2[i-k] + I2[i-k-1]);
            k++;
        }
    }

    int n2 = n1 - k;
    if(n2 < 1) {
        msg << "Insufficient points to calculate the current (m = " << n2 << ")" << endl;
850
        currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
gsell's avatar
gsell committed
851 852 853
    } else {
        // 3. smooth further
        if(n2 > 40) {
854
            sgSmooth(&I2.front(), n2, n2 / 20, n2 / 20, 0, 1);
gsell's avatar
gsell committed
855 856 857
        }

        // 4. create current profile
858 859
        currentProfile_m = std::unique_ptr<Profile>(new Profile(&z2.front(),
                                                    &I2.front(), n2));
gsell's avatar
gsell committed
860 861 862 863 864 865 866 867 868 869

        /**5. Normalize profile to match bunch charge as a constant
         * However, only normalize for sufficient beam energy
         */
        thisBunch = this;

        double Qcalc = 0.0;
        double z = zMin;
        dz = (zMax - zMin) / 99;
        for(int i = 1; i < 100; i++) {
870
            Qcalc += currentProfile_m->get(z, itype_lin);
gsell's avatar
gsell committed
871 872 873
            z += dz;
        }
        Qcalc *= (dz / (beta * Physics::c));
874
        currentProfile_m->scale((Q_m > 0.0 ? Q_m : Physics::q_e) / Qcalc);
gsell's avatar
gsell committed
875 876
    }

877
    dStat_m |= ds_currentCalculated;
gsell's avatar
gsell committed
878 879 880 881 882 883 884 885 886 887 888 889 890
}

void EnvelopeBunch::cSpaceCharge() {
    Inform msg("cSpaceCharge");

#ifdef USE_HOMDYN_SC_MODEL
    int icON = 1;
    double zhead = zHead();
    double ztail = zTail();
    double L = zhead - ztail;

    for(int i = 0; i < numMySlices_m; i++) {

891 892 893 894
        if(slices_m[i]->p[SLI_z] > 0.0)  {
            double zeta = slices_m[i]->p[SLI_z] - ztail;
            double xi   = slices_m[i]->p[SLI_z] + zhead;
            double sigma_av = (slices_m[i]->p[SLI_x] + slices_m[i]->p[SLI_y]) / 2;
gsell's avatar
gsell committed
895 896 897
            double R = 2 * sigma_av;
            double A = R / L / getGamma(i);

898 899
            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);
gsell's avatar
gsell committed
900 901

            //FIXME: Q_act or Q?
902
            Esct_m[i] = (Q_m / 2 / Physics::pi / Physics::epsilon_0 / R / R) * (H1 - icON * H2);
gsell's avatar
gsell committed
903 904 905
            double G1 = (1 - zeta / L) / sqrt((1 - zeta / L) * (1 - zeta / L) + A * A) + (zeta / L) / sqrt((zeta / L) * (zeta / L) + A * A);
            double G2 = (1 - xi / L) / sqrt((1 - xi / L) * (1 - xi / L) + A * A) + (xi / L) / sqrt((xi / L) * (xi / L) + A * A);

906
            G_m[i] = (1 - getBeta(i) * getBeta(i)) * G1 - icON * (1 + getBeta(i) * getBeta(i)) * G2;
gsell's avatar
gsell committed
907 908 909 910 911 912 913 914
        }
    }
#else
    if(numSlices_m < 2) {
        msg << "called with insufficient slices (" << numSlices_m << ")" << endl;
        return;
    }

915
    if((Q_m <= 0.0) || (currentProfile_m->max() <= 0.0)) {
gsell's avatar
gsell committed
916 917 918 919 920 921 922 923 924
        msg << "Q or I_max <= 0.0, aborting calculation" << endl;
        return;
    }

    for(int i = 0; i < numMySlices_m; ++i) {
        Esct[i] = 0.0;
        G[i]    = 0.0;
    }

925
    double Imax = currentProfile_m->max();
gsell's avatar
gsell committed
926 927 928
    int nV = 0;
    double sm = 0.0;
    double A0 = 0.0;
929
    std::vector<double> xi(numSlices_m, 0.0);
gsell's avatar
gsell committed
930 931

    for(int i = 0; i < numMySlices_m; i++) {
932
        if(slices_m[i]->p[SLI_z] > 0.0) {
gsell's avatar
gsell committed
933
            nV++;
934
            A0 = 4.0 * slices_m[i]->p[SLI_x] * slices_m[i]->p[SLI_y];
gsell's avatar
gsell committed
935
            sm += A0;
936
            xi[i+mySliceStartOffset_m] = A0 * (1.0 - slices_m[i]->p[SLI_beta] * slices_m[i]->p[SLI_beta]); // g2
gsell's avatar
gsell committed
937 938 939 940
        }
    }

    int nVTot = nV;
941
    allreduce(&nVTot, 1, std::plus<int>());
gsell's avatar
gsell committed
942 943 944 945 946
    if(nVTot < 2) {
        msg << "Exiting, to few nV slices" << endl;
        return;
    }

frey_m's avatar
frey_m committed
947
    allreduce(xi.data(), numSlices_m, std::plus<double>());
948
    allreduce(&sm, 1, std::plus<double>());
gsell's avatar
gsell committed
949 950 951 952 953
    A0 = sm / nVTot;
    double dzMin = 5.0 * Physics::c * Q_m / (Imax * numSlices_m);

    for(int localIdx = 0; localIdx < numMySlices_m; localIdx++) {
        double z0 = z_m[localIdx + mySliceStartOffset_m];
954
        if(z0 > 0.0 && slices_m[localIdx]->p[SLI_beta] > BETA_MIN1) {
gsell's avatar
gsell committed
955 956 957 958 959

            sm = 0.0;

            for(int j = 0; j < numSlices_m; j++) {
                double zj = z_m[j];
960
                double dz = std::abs(zj - z0);
gsell's avatar
gsell committed
961 962 963 964 965 966 967 968 969 970 971 972
                if((dz > dzMin) && (zj > zCat)) {
                    double Aj = xi[j] / (dz * dz);
                    double v  = 1.0 - (1.0 / sqrt(1.0 + Aj));
                    if(zj > z0) {
                        sm -= v;
                    } else {
                        sm += v;
                    }
                }
            }

            // longitudinal effect
973
            double bt = slices_m[localIdx]->p[SLI_beta];
gsell's avatar
gsell committed
974 975 976
            double btsq = (bt - BETA_MIN1) / (BETA_MIN2 - BETA_MIN1) * (bt - BETA_MIN1) / (BETA_MIN2 - BETA_MIN1);
            Esct[localIdx] = (bt < BETA_MIN1 ? 0.0 : bt < BETA_MIN2 ? btsq : 1.0);
            Esct[localIdx] *= Q_m * sm / (Physics::two_pi * Physics::epsilon_0 * A0 * (nVTot - 1));
977
            G[localIdx] = currentProfile_m->get(z0, itype_lin) / Imax;
gsell's avatar
gsell committed
978 979 980

            // tweak to compensate for non-relativity
            if(bt < BETA_MIN2) {
981
                if(slices_m[localIdx]->p[SLI_beta] < BETA_MIN1)
gsell's avatar
gsell committed
982 983 984 985 986 987 988 989 990 991 992 993
                    G[localIdx] = 0.0;
                else
                    G[localIdx] *= btsq;
            }
        }
    }

    return;
#endif
}

double EnvelopeBunch::moveZ0(double zC) {
994
    zCat_m = zC;
gsell's avatar
gsell committed
995 996 997
    double dz = zC - zHead();
    if(dz > 0.0) {
        for(int i = 0; i < numMySlices_m; i++) {
998
            slices_m[i]->p[SLI_z] += dz;
gsell's avatar
gsell committed
999 1000
        }
        backup(); // save the new state
1001
        *gmsg << "EnvelopeBunch::moveZ0(): bunch moved with " << dz << " m to " << zCat_m << " m" << endl;
gsell's avatar
gsell committed
1002 1003 1004 1005 1006 1007 1008 1009 1010
    }

    return dz;
}

double EnvelopeBunch::tReset(double dt) {
    double new_dt = dt;

    if(dt == 0.0) {
1011 1012
        new_dt = t_m;
        *gmsg << "EnvelopeBunch time reset at z = " << zAvg() << " m with: " << t_m << " s, new offset: " << t_m + t_offset_m << " s";
gsell's avatar
gsell committed
1013 1014
    }

1015 1016
    t_offset_m += new_dt;
    t_m        -= new_dt;
gsell's avatar
gsell committed
1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043

    return new_dt;
}

/** derivs for RK routine
 * Definition of equation set:
 *  ==========================
 Y[SLI_z]    = z      dz/dt  = beta*c*cos(a)
                      cos(a) = 1 - (px0^2 + py0^2)/c2
 Y[SLI_beta] = beta   db/dt  = (e0/mc)*(E_acc + E_sc)/gamma^3
 Y[SLI_x]    = x      dx/dt  = px = Y[SLI_px]
 Y[SLI_px]   = px     dpx/dt = f(x,beta) - (beta*gamma^2(db/dt)*px + Kr*x)
 Y[SLI_y]    = y      dy/dt  = py = Y[SLI_py]
 Y[SLI_py]   = py     dpy/dt = f(y,beta) - (beta*gamma^2(db/dt)*py + Kr*y)
 Y[SLI_x0]   = x0     dx0/dt = px0 = Y[SLI_px0]
 Y[SLI_px0]  = px0    dpx0/dt= -(beta*gamma^2(db/dt)*px0) + Kt*x0
 Y[SLI_y0]   = y0     dy0/dt = py0 = Y[SLI_py0]
 Y[SLI_py0]  = py0    dpy0/dt= -(beta*gamma^2(db/dt)*py0) + Kt*y0

 Transversal space charge blowup:
 f(x,beta) = c^2*I/(2Ia)/(x*beta*gamma^3)

ALT: SLI_z (commented by Rene)
    // dYdt[SLI_z] = Y[SLI_beta]*c*sqrt(1.0 - (pow(Y[SLI_px0],2) + pow(Y[SLI_py0],2))/pow(c*Y[SLI_beta],2));
    // dYdt[SLI_z] = Y[SLI_beta]*c*cos(Y[SLI_px0]/Y[SLI_beta]/c)*cos(Y[SLI_py0]/Y[SLI_beta]/c);

 */
gsell's avatar
gsell committed
1044
void EnvelopeBunch::derivs(double /*tc*/, double Y[], double dYdt[]) {
gsell's avatar
gsell committed
1045 1046 1047 1048 1049 1050 1051 1052
    double g2 = 1.0 / (1.0 - Y[SLI_beta] * Y[SLI_beta]);
    double g  = sqrt(g2);
    double g3 = g2 * g;

    double alpha = sqrt(Y[SLI_px0] * Y[SLI_px0] + Y[SLI_py0] * Y[SLI_py0]) / Y[SLI_beta] / Physics::c;
    /// \f[ \dot{z} = \beta c cos(\alpha) \f]
    dYdt[SLI_z]  = Y[SLI_beta] * Physics::c * cos(alpha);

1053
    //NOTE: (reason for - when using Esl(2)) In OPAL we use e_0 with a sign.
gsell's avatar
gsell committed
1054 1055 1056
    // The same issue with K factors (accounted in K factor calculation)
    //
    /// \f[ \dot{\beta} = \frac{e_0}{mc\gamma^3} \left(E_{z,\text{ext}} + E_{z,\text{sc}} + E_{z, \text{wake}} \right)
1057
    dYdt[SLI_beta] = Physics::e0mc * (-Esl_m(2) + Esct_m[currentSlice_m] + Ezw_m[currentSlice_m]) / g3;
gsell's avatar
gsell committed
1058 1059 1060 1061

    /// \f[ \beta * \gamma^2 * \dot{\beta} \f]
    double bg2dbdt = Y[SLI_beta] * g2 * dYdt[SLI_beta];

1062
    if(solver_m & sv_radial) {
gsell's avatar
gsell committed
1063 1064
        /// minimum spot-size due to emittance
        /// \f[ \left(\frac{\epsilon_n c}{\gamma}\right)^2 \f]
1065 1066
        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);
gsell's avatar
gsell committed
1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079


#ifdef USE_HOMDYN_SC_MODEL
        //FIXME: Q_act or Q?
        double kpc  = 0.5 * Physics::c * Physics::c * (Y[SLI_beta] * Physics::c) * activeSlices_m * Q_m / numSlices_m / (curZHead_m - curZTail_m) / Physics::Ia;

        /// \f[ \dot{\sigma} = p \f]
        dYdt[SLI_x] = Y[SLI_px];
        dYdt[SLI_y] = Y[SLI_py];

        double sigma_av = (Y[SLI_x] + Y[SLI_y]) / 2;

        /// \f[ \ddot{\sigma} = -\gamma^2\beta\dot{\beta}\dot{\sigma} - K\sigma + 2c^2\left(\frac{I}{2I_0}\right)\frac{G}{\beta R^2}(1-\beta)^2 \sigma + \left(\frac{\epsilon_n c}{\gamma}\right)^2 \frac{1}{\sigma^3} \f]
1080 1081
        dYdt[SLI_px] = -bg2dbdt * Y[SLI_px] - KRsl_m(0) * Y[SLI_x] + (kpc / sigma_av / Y[SLI_beta] / g / 2) * G_m[currentSlice_m] + enxc2 / g3;
        dYdt[SLI_py] = -bg2dbdt * Y[SLI_py] - KRsl_m(1) * Y[SLI_y] + (kpc / sigma_av / Y[SLI_beta] / g / 2) * G_m[currentSlice_m] + enyc2 / g3;
gsell's avatar
gsell committed
1082 1083 1084 1085
#else
        // transverse space charge
        // somewhat strange: I expected: c*c*I/(2*Ia) (R. Bakker)
        //XXX: Renes version, I[cs] already in G
1086
        double kpc  = 0.5 * Physics::c * Physics::c * currentProfile_m->max() / Physics::Ia;
gsell's avatar
gsell committed
1087 1088 1089 1090

        /// \f[ \dot{\sigma} = p \f]
        dYdt[SLI_x]  = Y[SLI_px];
        dYdt[SLI_y]  = Y[SLI_py];
1091 1092
        dYdt[SLI_px] = (kpc * G_m[currentSlice_m] / (Y[SLI_x] * Y[SLI_beta] * g3)) + enxc2 / g3 - (KRsl_m(0) * Y[SLI_x]) - (bg2dbdt * Y[SLI_px]);
        dYdt[SLI_py] = (kpc * G_m[currentSlice_m] / (Y[SLI_y] * Y[SLI_beta] * g3)) + enyc2 / g3 - (KRsl_m(1) * Y[SLI_y]) - (bg2dbdt * Y[SLI_py]);
gsell's avatar
gsell committed
1093 1094 1095 1096 1097 1098 1099 1100 1101
#endif
    } else {
        /// \f[ \dot{\sigma} = p \f]
        dYdt[SLI_x]  = Y[SLI_px];
        dYdt[SLI_y]  = Y[SLI_py];
        dYdt[SLI_px] = 0.0;
        dYdt[SLI_py] = 0.0;
    }

1102
    if(solver_m & sv_offaxis) {
gsell's avatar
gsell committed
1103 1104
        dYdt[SLI_x0]  = Y[SLI_px0];
        dYdt[SLI_y0]  = Y[SLI_py0];
1105 1106
        dYdt[SLI_px0] = -KTsl_m(0) - (bg2dbdt * Y[SLI_px0]) + Physics::e0m * (g * Exw_m[currentSlice_m]);
        dYdt[SLI_py0] = -KTsl_m(1) - (bg2dbdt * Y[SLI_py0]) + Physics::e0m * (g * Eyw_m[currentSlice_m]);
gsell's avatar
gsell committed
1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136
    } else {
        dYdt[SLI_x0]  = Y[SLI_px0];
        dYdt[SLI_y0]  = Y[SLI_py0];
        dYdt[SLI_px0] = 0.0;
        dYdt[SLI_py0] = 0.0;
    }
}


void EnvelopeBunch::computeSpaceCharge() {
    IpplTimings::startTimer(selfFieldTimer_m);

    // Calculate the current profile
    if(Q_m > 0.0) {

        IpplTimings::startTimer(calcITimer_m);
#ifndef USE_HOMDYN_SC_MODEL
        //XXX: only used in BET-SC-MODEL
        // make sure all processors have all global betas and positions
        synchronizeSlices();
        calcI();
#endif
        IpplTimings::stopTimer(calcITimer_m);

        // the following assumes space-charges do not change significantly over nSc steps
        IpplTimings::startTimer(spaceChargeTimer_m);
        cSpaceCharge();
        IpplTimings::stopTimer(spaceChargeTimer_m);

    } else {
1137
        currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
gsell's avatar
gsell committed
1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152
    }

    IpplTimings::stopTimer(selfFieldTimer_m);
}


void EnvelopeBunch::timeStep(double tStep, double _zCat) {
    Inform msg("tStep");
    static int msgParsed = 0;

    // default accuracy of integration
    double eps = 1.0e-4;
    double time_step_s = tStep;

    thisBunch = this;
1153
    zCat_m     = _zCat;
gsell's avatar
gsell committed
1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165

    // backup last stage before new execution
    backup();

    //FIXME: HAVE A PROBLEM WITH EMISSION (EMISSION OFF GIVES GOOD RESULTS)
    activeSlices_m = numSlices_m;

    curZHead_m = zHead();
    curZTail_m = zTail();

    for(int i = 0; i < numMySlices_m; i++) {
        // make the current slice index global in this class
1166 1167
        currentSlice_m = i;
        std::shared_ptr<EnvelopeSlice> sp = slices_m[i];
gsell's avatar
gsell committed
1168 1169

        // Assign values of K for certain slices
1170 1171 1172 1173
        KRsl_m = KR[i];
        KTsl_m = KT[i];
        Esl_m = EF[i];
        Bsl_m = BF[i];
gsell's avatar
gsell committed
1174 1175 1176 1177 1178 1179

        // only for slices already emitted
        if(true /*sp->emitted*/) {
            // set default allowed error for integration
            double epsLocal = eps;
            // mark that the integration was not successful yet
gsell's avatar
gsell committed
1180
            bool ode_result = false;
gsell's avatar
gsell committed
1181

gsell's avatar
gsell committed
1182
            while(ode_result == false) {
gsell's avatar
gsell committed
1183

1184 1185
                if(solver_m & sv_fixedStep) {
                    rk4(&(sp->p[SLI_z]), SLNPAR, t_m, time_step_s, Gderivs);
gsell's avatar
gsell committed
1186
                    ode_result = true;
gsell's avatar
gsell committed
1187 1188
                } else {
                    int nok, nbad;
1189
                    activeSlice_m = i;
gsell's avatar
gsell committed
1190
                    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);
gsell's avatar
gsell committed
1191 1192
                }

gsell's avatar
gsell committed
1193
                if(ode_result == false) {
gsell's avatar
gsell committed
1194 1195 1196 1197 1198
                    // restore the backup
                    sp->restore();
                    epsLocal *= 10.0;
                }
            }
gsell's avatar
gsell committed
1199 1200
            // :FIXME: the above loop cannot exit with ode_result == false
            if(ode_result == false) {
gsell's avatar
gsell committed
1201
                // use fixed step integration if dynamic fails
1202
                rk4(&(sp->p[SLI_z]), SLNPAR, t_m, time_step_s, Gderivs);
gsell's avatar
gsell committed
1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213

                if(msgParsed < 2) {
                    msg << "EnvelopeBunch::run() Switched to fixed step RK rountine for solving of DE at slice " << i << endl;
                    msgParsed = 2;
                }
            } else if((epsLocal != eps) && (msgParsed == 0)) {
                msg << "EnvelopeBunch::run() integration accuracy relaxed to " << epsLocal << " for slice " << i << " (ONLY FIRST OCCURENCE MARKED!)" << endl;
                msgParsed = 1;
            }
        }

1214 1215 1216
        if(slices_m[i]->check()) {
            msg << "Slice " << i << " no longer valid at z = " <<  slices_m[i]->p_old[SLI_z] << " beta = " << slices_m[i]->p_old[SLI_beta] << endl;
            msg << "Slice " << i << " no longer valid at z = " <<  slices_m[i]->p[SLI_z] << " beta = " << slices_m[i]->p[SLI_beta] << endl;
gsell's avatar
gsell committed
1217 1218 1219 1220 1221
            isValid_m = false;
            return;
        }
    }
    // mark that slices might not be synchronized (and space charge accordingly)