EnvelopeBunch.cpp 54.6 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

// 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); }
51
//static double Gcur(double z) { return thisBunch->currentProfile_m->get(z, itype_lin); }
gsell's avatar
gsell committed
52 53 54 55
static double rootValue = 0.0;

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

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

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

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


EnvelopeBunch::~EnvelopeBunch() {
}

gsell's avatar
gsell committed
87 88 89 90
void EnvelopeBunch::calcBeamParameters() {
    Inform msg("calcBeamParameters");
    IpplTimings::startTimer(statParamTimer_m);
    double ex, ey, nex, ney, b0, bRms, bMax, bMin, g0, dgdt, gInc;
91
    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
92 93 94 95 96 97 98 99
    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);

100
    if(solver_m & sv_radial) {
gsell's avatar
gsell committed
101 102 103 104 105 106 107
        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);
    }

108
    if(solver_m & sv_offaxis) {
gsell's avatar
gsell committed
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
        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;
126 127 128 129
    // dx02_m = dx0_m;
    // dy02_m = dy0_m;
    // dfi_x_m = 180.0 * dfi_x / Physics::pi;
    // dfi_y_m = 180.0 * dfi_y / Physics::pi;
gsell's avatar
gsell committed
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 155 156 157 158 159 160
    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);
    //minP_m = Vector_t(-maxP_m[0], -maxP_m[1], 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;
161
    std::vector<double> v(numMySlices_m, 0.0);
gsell's avatar
gsell committed
162 163 164 165 166

    //FIXME: why from 1 to n-1??
    switch(sp) {
        case sp_beta:      // normalized velocity (total) [-]
            for(i = 1; i < numMySlices_m - 1; i++) {
167
                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
168 169 170 171
            }
            break;
        case sp_gamma:     // Lorenz factor
            for(i = 1; i < numMySlices_m - 1; i++) {
172
                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
173 174 175 176
            }
            break;
        case sp_z:         // slice position [m]
            for(i = 0; i < numMySlices_m - 0; i++) {
177
                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
178 179 180 181
            }
            break;
        case sp_I:         // slice position [m]
            for(i = 1; i < numMySlices_m - 1; i++) {
182 183
                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
184 185 186 187
            }
            break;
        case sp_Rx:        // beam size x [m]
            for(i = 0; i < numMySlices_m - 0; i++) {
188
                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
189 190 191 192
            }
            break;
        case sp_Ry:        // beam size y [m]
            for(i = 0; i < numMySlices_m - 0; i++) {
193
                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
194 195 196 197
            }
            break;
        case sp_Px:        // beam divergence x
            for(i = 0; i < numMySlices_m - 0; i++) {
198
                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
199 200 201 202
            }
            break;
        case sp_Py:        // beam divergence y
            for(i = 0; i < numMySlices_m - 0; i++) {
203
                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
204 205 206 207
            }
            break;
        case sp_Pz:
            for(i = 1; i < numMySlices_m - 1; i++) {
208
                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
209 210 211 212
            }
            break;
        case sp_x0:        // position centroid x [m]
            for(i = 1; i < numMySlices_m - 1; i++) {
213
                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
214 215 216 217
            }
            break;
        case sp_y0:        // position centroid y [m]
            for(i = 1; i < numMySlices_m - 1; i++) {
218
                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
219 220 221 222
            }
            break;
        case sp_px0:       // angular deflection centroid x
            for(i = 1; i < numMySlices_m - 1; i++) {
223
                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
224 225 226 227
            }
            break;
        case sp_py0:      // angular deflection centroid y
            for(i = 1; i < numMySlices_m - 1; i++) {
228
                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
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 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 300 301 302 303 304 305
            }
            break;
        default :
            throw OpalException("EnvelopeBunch", "EnvelopeBunch::runStats() Undefined label (programming error)");
            break;
    }

    int nVTot = nV;
    MPI_Allreduce(MPI_IN_PLACE, &nVTot, 1, MPI_INT, MPI_SUM, Ippl::getComm());
    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]);
            }
        }

        MPI_Allreduce(MPI_IN_PLACE, &M1, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
        MPI_Allreduce(MPI_IN_PLACE, &M2, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
        MPI_Allreduce(MPI_IN_PLACE, &maxv, 1, MPI_DOUBLE, MPI_MAX, Ippl::getComm());
        MPI_Allreduce(MPI_IN_PLACE, &minv, 1, MPI_DOUBLE, MPI_MIN, Ippl::getComm());

        *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++) {
306
        if((slices_m[i]->p[SLI_z] > zCat_m) && slices_m[i]->is_valid())
gsell's avatar
gsell committed
307 308 309 310 311 312 313 314
            nV++;
    }

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

318
                if(solver_m & sv_radial) {
gsell's avatar
gsell committed
319
                    assert(i < numMySlices_m);
320
                    auto sp = slices_m[i];
321
                    bg = sp->p[SLI_beta] * sp->computeGamma();
gsell's avatar
gsell committed
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 356 357 358 359 360 361

                    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;

362 363
        *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
364 365 366 367 368 369 370 371 372

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

gsell's avatar
gsell committed
373
void EnvelopeBunch::calcEnergyChirp(double *g0, double *dgdt, double *gInc, int */*nValid*/) {
374 375 376
    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
377 378 379 380 381 382 383 384 385 386

    // 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++) {
387 388
        std::shared_ptr<EnvelopeSlice> cs = slices_m[i];
        if((cs->is_valid()) && (cs->p[SLI_z] > zCat_m)) {
gsell's avatar
gsell committed
389 390 391
            zAvg    += cs->p[SLI_z];
            dtl[i-j]  = cs->p[SLI_z];
            b[i-j]   = cs->p[SLI_beta];
392
            g[i-j]   = cs->computeGamma();
gsell's avatar
gsell committed
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
            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) {
412 413 414
        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
415 416

        int numproc = Ippl::Comm->getNodes();
417 418 419 420
        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
421

422
        MPI_Allgather(&nV, 1, MPI_INT, &offsets.front(), 1, MPI_INT, Ippl::getComm());
gsell's avatar
gsell committed
423 424 425 426 427 428 429 430
        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];
        }

431 432 433 434 435 436 437
        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
438

439 440 441
        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
442 443 444 445 446 447 448 449 450

        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
451
        linfit(&dtG[0], &gG[0], nVTot, gZero, gt, dum2, dum3, dum4);
gsell's avatar
gsell committed
452 453 454 455
        *dgdt = gt;

        rms = 0.0;
        for(int i = 0; i < nVTot; i++) {
456
            rms += std::pow(gG[i] - gZero - gt * dtG[i], 2);
gsell's avatar
gsell committed
457 458 459 460 461 462 463
        }

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


464
void EnvelopeBunch::distributeSlices(int nSlice) {
gsell's avatar
gsell committed
465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490
    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;
}


491 492 493
void EnvelopeBunch::createBunch() {

    // Memory issue when going to larger number of processors (this is
gsell's avatar
gsell committed
494 495
    // also the reason for the 30 slices hard limit per core!)
    // using heuristic (until problem is located)
496 497 498 499 500 501 502 503 504 505 506
    //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
507 508 509 510 511 512 513 514 515

    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
516
    solver_m = sv_radial | sv_offaxis | sv_lwakes | sv_twakes;
gsell's avatar
gsell committed
517

518
    //FIXME: WHY?
gsell's avatar
gsell committed
519 520 521 522
    if(numSlices_m < 14)
        throw OpalException("EnvelopeBunch::createSlices", "use more than 13 slices");

    // defaults:
523 524 525 526 527 528 529 530 531 532 533 534 535 536 537
    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 ) {
538 539 540 541 542 543 544 545 546
        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());
        //});

547 548 549 550 551
    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);
552

gsell's avatar
gsell committed
553
    for(unsigned int i = 0; i < nSlices; i++) {
554 555 556 557 558
        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
559 560
    }

561
    currentProfile_m = NULL;
562 563
    I0avg_m = 0.0;
    dStat_m = ds_fieldsSynchronized | ds_slicesSynchronized;
gsell's avatar
gsell committed
564 565 566 567 568 569 570 571 572 573 574 575 576 577
}

//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:

578
            bunch_width = Physics::c * emission_time_s * slices_m[0]->p[SLI_beta];
kraus's avatar
kraus committed
579
//            std::cout << "bunch_width = " << bunch_width << " SLI_beta= " << slices_m[0]->p[SLI_beta] << std::endl;
gsell's avatar
gsell committed
580
            for(int i = 0; i < numMySlices_m; i++) {
581
                slices_m[i]->p[SLI_z] = -(((numSlices_m - 1) - (mySliceStartOffset_m + i)) * bunch_width) / numSlices_m;
gsell's avatar
gsell committed
582
            }
583
            I0avg_m = Q_m * Physics::c / std::abs(2.0 * emission_time_s);
gsell's avatar
gsell committed
584 585 586 587
            break;

        case bsGauss:
            if(n2 >= mySliceStartOffset_m && n2 <= mySliceEndOffset_m)
588
                slices_m[n2-mySliceStartOffset_m]->p[SLI_z] = z0;
gsell's avatar
gsell committed
589 590 591

            for(int i = 1; i <= numSlices_m / 2; i++) {
                rootValue = 1.0 - 2.0 * i * frac / (numSlices_m + 1);
592
                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
593 594

                if((n2 + i) >= mySliceStartOffset_m && (n2 + i) <= mySliceEndOffset_m)
595
                    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
596 597

                if((n2 - i) >= mySliceStartOffset_m && (n2 - i) <= mySliceEndOffset_m)
598
                    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
599 600
            }

601
            I0avg_m = 0.0;
gsell's avatar
gsell committed
602 603 604 605 606
            break;
    }

    double gz0 = 0.0, gzN = 0.0;
    if(Ippl::Comm->myNode() == 0)
607
        gz0 = slices_m[0]->p[SLI_z];
gsell's avatar
gsell committed
608
    if(Ippl::Comm->myNode() == Ippl::Comm->getNodes() - 1)
609
        gzN = slices_m[numMySlices_m-1]->p[SLI_z];
gsell's avatar
gsell committed
610 611 612
    reduce(gz0, gz0, OpAddAssign());
    reduce(gzN, gzN, OpAddAssign());

613
    hbin_m = (gzN - gz0) / nebin_m;
gsell's avatar
gsell committed
614 615 616 617 618 619 620 621 622 623

    // 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) {
624
        if((bin_i + 1) * hbin_m < slices_m[slice_i]->p[SLI_z] - gz0) {
gsell's avatar
gsell committed
625 626 627 628 629 630 631 632
            bin_i++;
        } else {
            this->bins_m[bin_i].push_back(slice_i);
            slice_i++;
        }
    }

    ////XXX: debug output
adelmann's avatar
adelmann committed
633 634 635
    /*
    for(unsigned int i = 0; i < bins_m.size(); i++) {
      if(bins_m[i].size() > 0) {
kraus's avatar
kraus committed
636 637 638 639
    std::cout << Ippl::Comm->myNode() << ": Bin " << i << ": ";
    for(unsigned int j = 0; j < bins_m[i].size(); j++)
      std::cout << " " << mySliceStartOffset_m + bins_m[i][j] << "(" << slices_m[j]->p[SLI_z] << ")";
    std::cout << std::endl;
adelmann's avatar
adelmann committed
640 641 642
      }
    }
    */
gsell's avatar
gsell committed
643 644 645 646 647 648 649 650
    //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
651
void EnvelopeBunch::setTShape(double enx, double eny, double rx, double ry, double /*b0*/) {
adelmann's avatar
adelmann committed
652
  /*
gsell's avatar
gsell committed
653 654 655
    Inform msg("setTshape");
    msg << "set SLI_x to " << rx / 2.0 << endl;
    msg << "set SLI_y to " << ry / 2.0 << endl;
adelmann's avatar
adelmann committed
656
  */
gsell's avatar
gsell committed
657
    // set emittances
658 659
    emtnx0_m = enx;
    emtny0_m = eny;
gsell's avatar
gsell committed
660
    //FIXME: is rx here correct?
661 662
    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);
adelmann's avatar
adelmann committed
663 664
    //msg << "set emtbx0 to " << emtbx0 << endl;
    //msg << "set emtby0 to " << emtby0 << endl;
gsell's avatar
gsell committed
665 666 667

    //XXX: rx = 2*rms_x
    for(int i = 0; i < numMySlices_m; i++) {
668 669 670 671
        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
672 673 674 675 676 677 678
    }

    backup();
}

void EnvelopeBunch::setTOffset(double x0, double px0, double y0, double py0) {
    for(int i = 0; i < numMySlices_m; i++) {
679 680 681 682
        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
683 684 685 686 687
    }
}

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

    for(int i = 0; i < numMySlices_m; i++) {
692 693
        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
694 695 696 697 698 699 700 701 702 703 704 705
    }

    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++) {
706 707
        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
708 709 710 711 712 713 714 715 716 717
    }

    MPI_Allreduce(MPI_IN_PLACE, &(z_m[0]), numSlices_m, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
    MPI_Allreduce(MPI_IN_PLACE, &(b_m[0]), numSlices_m, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
}

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

    static int already_called = 0;
718
    if((dStat_m & ds_currentCalculated) || (already_called && (Q_m <= 0.0)))
gsell's avatar
gsell committed
719 720 721
        return;
    already_called = 1;

722
    std::vector<double> z1(numSlices_m, 0.0);
723
    std::vector<double> b (numSlices_m, 0.0);
gsell's avatar
gsell committed
724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749
    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;
        //throw OpalException("EnvelopeBunch", "Insufficient points to calculate the current (n1)");
750
        currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
gsell's avatar
gsell committed
751 752 753 754 755 756
        return;
    }

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

757 758 759 760 761 762 763 764 765 766 767 768 769 770
    //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
771 772 773 774 775

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

    // 1. Determine current from the slice distance
776
    std::vector<double> I1(n1, 0.0);
gsell's avatar
gsell committed
777 778 779 780 781 782 783 784 785 786 787 788 789

    //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
790
    int i, j;
gsell's avatar
gsell committed
791 792 793
#ifdef _OPENMP
#pragma omp parallel for private(i,j,dz)
#endif //_OPENMP
gsell's avatar
gsell committed
794 795
    for(i = my_start; i <= vend; i++) {
        j = 1;
gsell's avatar
gsell committed
796
        do {
797
            dz = std::abs(z1[i+j] - z1[i-j]);
gsell's avatar
gsell committed
798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820
            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;
        }
    }

    MPI_Allreduce(MPI_IN_PLACE, &(I1[0]), n1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
    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
821 822
    std::vector<double> z2(n1, 0.0);
    std::vector<double> I2(n1, 0.0);
gsell's avatar
gsell committed
823 824 825
    double Mz1 = 0.0;
    double MI1 = 0.0;
    int np = 0;
gsell's avatar
gsell committed
826
    j = 0;
gsell's avatar
gsell committed
827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896

    // 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++) {
        //for(int i = my_start; i <= my_end; 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_temp[i] = Mz1 / np;
        //I2_temp[i] = MI1 / np;
        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++;
        }
    }

    //MPI_Allreduce(MPI_IN_PLACE, &(z2_temp[0]), n1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
    //MPI_Allreduce(MPI_IN_PLACE, &(I2_temp[0]), n1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());

    ////FIXME: we dont need copy of z2 and I2: z2[i-k] = z2[i];
    //int k = 0;
    //for(int i = 1; i < n1; i++) {
    //z2[i-k] = z2_temp[i];
    //I2[i-k] = I2_temp[i];
    //// 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++;
    //}
    //}

    //free(z2_temp);
    //free(I2_temp);

    int n2 = n1 - k;
    if(n2 < 1) {
        msg << "Insufficient points to calculate the current (m = " << n2 << ")" << endl;
897
        currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
gsell's avatar
gsell committed
898 899 900
    } else {
        // 3. smooth further
        if(n2 > 40) {
901
            sgSmooth(&I2.front(), n2, n2 / 20, n2 / 20, 0, 1);
gsell's avatar
gsell committed
902 903 904
        }

        // 4. create current profile
905 906
        currentProfile_m = std::unique_ptr<Profile>(new Profile(&z2.front(),
                                                    &I2.front(), n2));
gsell's avatar
gsell committed
907 908 909 910 911 912 913 914 915 916

        /**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++) {
917
            Qcalc += currentProfile_m->get(z, itype_lin);
gsell's avatar
gsell committed
918 919 920
            z += dz;
        }
        Qcalc *= (dz / (beta * Physics::c));
921
        currentProfile_m->scale((Q_m > 0.0 ? Q_m : Physics::q_e) / Qcalc);
gsell's avatar
gsell committed
922 923
    }

924
    dStat_m |= ds_currentCalculated;
gsell's avatar
gsell committed
925 926 927 928 929 930 931 932 933 934 935 936 937
}

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++) {

938 939 940 941
        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
942 943 944
            double R = 2 * sigma_av;
            double A = R / L / getGamma(i);

945 946
            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
947 948 949

            //FIXME: Q_act or Q?
            //double Q = activeSlices_m * Q_m / numSlices_m;
950
            Esct_m[i] = (Q_m / 2 / Physics::pi / Physics::epsilon_0 / R / R) * (H1 - icON * H2);
gsell's avatar
gsell committed
951 952 953
            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);

954
            G_m[i] = (1 - getBeta(i) * getBeta(i)) * G1 - icON * (1 + getBeta(i) * getBeta(i)) * G2;
gsell's avatar
gsell committed
955 956 957 958 959 960 961 962
        }
    }
#else
    if(numSlices_m < 2) {
        msg << "called with insufficient slices (" << numSlices_m << ")" << endl;
        return;
    }

963
    if((Q_m <= 0.0) || (currentProfile_m->max() <= 0.0)) {
gsell's avatar
gsell committed
964 965 966 967 968 969 970 971 972
        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;
    }

973
    double Imax = currentProfile_m->max();
gsell's avatar
gsell committed
974 975 976
    int nV = 0;
    double sm = 0.0;
    double A0 = 0.0;
977
    std::vector<double> xi(numSlices_m, 0.0);
gsell's avatar
gsell committed
978 979

    for(int i = 0; i < numMySlices_m; i++) {
980
        if(slices_m[i]->p[SLI_z] > 0.0) {
gsell's avatar
gsell committed
981
            nV++;
982
            A0 = 4.0 * slices_m[i]->p[SLI_x] * slices_m[i]->p[SLI_y];
gsell's avatar
gsell committed
983
            sm += A0;
984
            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
985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001
        }
    }

    int nVTot = nV;
    MPI_Allreduce(MPI_IN_PLACE, &nVTot, 1, MPI_INT, MPI_SUM, Ippl::getComm());
    if(nVTot < 2) {
        msg << "Exiting, to few nV slices" << endl;
        return;
    }

    MPI_Allreduce(MPI_IN_PLACE, &xi[0], numSlices_m, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
    MPI_Allreduce(MPI_IN_PLACE, &sm, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
    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];
1002
        if(z0 > 0.0 && slices_m[localIdx]->p[SLI_beta] > BETA_MIN1) {
gsell's avatar
gsell committed
1003 1004 1005 1006 1007

            sm = 0.0;

            for(int j = 0; j < numSlices_m; j++) {
                double zj = z_m[j];
1008
                double dz = std::abs(zj - z0);
gsell's avatar
gsell committed
1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020
                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
1021
            double bt = slices_m[localIdx]->p[SLI_beta];
gsell's avatar
gsell committed
1022 1023 1024
            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));
1025
            G[localIdx] = currentProfile_m->get(z0, itype_lin) / Imax;
gsell's avatar
gsell committed
1026 1027 1028

            // tweak to compensate for non-relativity
            if(bt < BETA_MIN2) {
1029
                if(slices_m[localIdx]->p[SLI_beta] < BETA_MIN1)
gsell's avatar
gsell committed
1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041
                    G[localIdx] = 0.0;
                else
                    G[localIdx] *= btsq;
            }
        }
    }

    return;
#endif
}

double EnvelopeBunch::moveZ0(double zC) {
1042
    zCat_m = zC;
gsell's avatar
gsell committed
1043 1044 1045
    double dz = zC - zHead();
    if(dz > 0.0) {
        for(int i = 0; i < numMySlices_m; i++) {
1046
            slices_m[i]->p[SLI_z] += dz;
gsell's avatar
gsell committed
1047 1048
        }
        backup(); // save the new state
1049
        *gmsg << "EnvelopeBunch::moveZ0(): bunch moved with " << dz << " m to " << zCat_m << " m" << endl;
gsell's avatar
gsell committed
1050 1051 1052 1053 1054 1055 1056 1057 1058
    }

    return dz;
}

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

    if(dt == 0.0) {
1059 1060
        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
1061 1062
    }

1063 1064
    t_offset_m += new_dt;
    t_m        -= new_dt;
gsell's avatar
gsell committed
1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091

    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
1092
void EnvelopeBunch::derivs(double /*tc*/, double Y[], double dYdt[]) {
gsell's avatar
gsell committed
1093 1094 1095 1096 1097 1098 1099 1100
    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);

1101
    //NOTE: (reason for - when using Esl(2)) In OPAL we use e_0 with a sign.
gsell's avatar
gsell committed
1102 1103 1104
    // 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)
1105
    dYdt[SLI_beta] = Physics::e0mc * (-Esl_m(2) + Esct_m[currentSlice_m] + Ezw_m[currentSlice_m]) / g3;
gsell's avatar
gsell committed
1106 1107 1108 1109

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

1110
    if(solver_m & sv_radial) {
gsell's avatar
gsell committed
1111 1112
        /// minimum spot-size due to emittance
        /// \f[ \left(\frac{\epsilon_n c}{\gamma}\right)^2 \f]
1113 1114
        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
1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128


#ifdef USE_HOMDYN_SC_MODEL
        //FIXME: Q_act or Q?
        //double Q = activeSlices_m * Q_m / numSlices_m;
        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]
1129 1130
        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
1131 1132 1133 1134
#else
        // transverse space charge
        // somewhat strange: I expected: c*c*I/(2*Ia) (R. Bakker)
        //XXX: Renes version, I[cs] already in G
1135
        double kpc  = 0.5 * Physics::c * Physics::c * currentProfile_m->max() / Physics::Ia;
gsell's avatar
gsell committed
1136 1137 1138 1139

        /// \f[ \dot{\sigma} = p \f]
        dYdt[SLI_x]  = Y[SLI_px];
        dYdt[SLI_y]  = Y[SLI_py];
1140 1141
        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
1142 1143 1144 1145 1146 1147 1148 1149 1150
#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;
    }

1151
    if(solver_m & sv_offaxis) {
gsell's avatar
gsell committed
1152 1153
        dYdt[SLI_x0]  = Y[SLI_px0];
        dYdt[SLI_y0]  = Y[SLI_py0];
1154 1155
        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
1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190
    } 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);

        //if(numSlices_m > 40) {
        // smooth Esct to get rid of numerical noise
        //double Esc = 0;
        //sgSmooth(Esct,n,n/20,n/20,0,4);
        //}
    } else {
1191
        currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
gsell's avatar
gsell committed
1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206
    }

    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;
1207
    zCat_m     = _zCat;
gsell's avatar
gsell committed
1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222

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

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

    //if(lastEmittedBin_m < nebin_m) {

    //time_step_s = emission_time_step_;
    //size_t nextBin = nebin_m - 1 - lastEmittedBin_m;
    //if(Ippl::Comm->myNode() == 0) cout << "Emitting bin " << nextBin << endl;

    //double gz0 = 0.0;
    //if(Ippl::Comm->myNode() == 0)
1223
    //gz0 = slices_m[0]->p[SLI_z];
gsell's avatar
gsell committed
1224 1225 1226 1227
    //MPI_Allreduce(MPI_IN_PLACE, &gz0, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());

    ////XXX: bin holds local slice number
    //for(int j = 0; j < bins_m[nextBin].size(); j++) {
1228 1229 1230 1231
    //slices_m[bins_m[nextBin][j]]->p[SLI_z] -= gz0;
    //slices_m[bins_m[nextBin][j]]->p[SLI_z] -= hbin_m * nextBin;
    //slices_m[bins_m[nextBin][j]]->emitted = true;
    //cout << "\tEmitting slice " << mySliceStartOffset_m + bins_m[nextBin][j] << " at " << slices_m[bins_m[nextBin][j]]->p[SLI_z] << endl;
gsell's avatar
gsell committed
1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248
    //}

    //lastEmittedBin_m++;
    //activeSlices_m += bins_m[nextBin].size();
    //}

    //activeSlices_m = min(activeSlices_m+1, numSlices_m);
    //if(activeSlices_m != numSlices_m)
    //time_step_s = emission_time_step_;
    //else
    //time_step_s = tStep;

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

    for(int i = 0; i < numMySlices_m; i++) {
        // make the current slice index global in this class
1249 1250
        currentSlice_m = i;
        std::shared_ptr<EnvelopeSlice> sp = slices_m[i];
gsell's avatar
gsell committed
1251 1252 1253 1254 1255 1256 1257 1258

        //if(i + mySliceStartOffset_m == numSlices_m - activeSlices_m) {
        //sp->emitted = true;
        //sp->p[SLI_z] = 0.0;
        //std::cout << "emitting " << i + mySliceStartOffset_m << std::endl;
        //}

        // Assign values of K for certain slices
1259 1260 1261 1262
        KRsl_m = KR[i];
        KTsl_m = KT[i];
        Esl_m = EF[i];
        Bsl_m = BF[i];
gsell's avatar
gsell committed
1263 1264 1265 1266 1267 1268

        // 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
1269