Bend2D.cpp 64.1 KB
Newer Older
1
// ------------------------------------------------------------------------
2
// $RCSfile: Bend2D.cpp,v $
3 4 5 6 7 8
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
9
// Definitions for class: Bend2D
10 11 12 13 14 15 16 17 18 19 20
//   Defines the abstract interface for a sector bend magnet.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------

21
#include "AbsBeamline/Bend2D.h"
22
#include "Algorithms/PartBunchBase.h"
23 24
#include "AbsBeamline/BeamlineVisitor.h"
#include "Utilities/Options.h"
25
#include "Utilities/Util.h"
26 27
#include "Fields/Fieldmap.h"
#include "AbstractObjects/OpalData.h"
28 29 30 31
#include "Structure/MeshGenerator.h"

#include "gsl/gsl_poly.h"

32 33 34 35 36
#include <iostream>
#include <fstream>

extern Inform *gmsg;

37
// Class Bend2D
38 39
// ------------------------------------------------------------------------

40
Bend2D::Bend2D():
41 42
    Bend2D("")
{}
43

44
Bend2D::Bend2D(const Bend2D &right):
45
    BendBase(right),
46 47 48 49 50 51 52
    messageHeader_m(" * "),
    pusher_m(right.pusher_m),
    designRadius_m(right.designRadius_m),
    exitAngle_m(right.exitAngle_m),
    fieldIndex_m(right.fieldIndex_m),
    startField_m(right.startField_m),
    endField_m(right.endField_m),
53 54
    widthEntranceFringe_m(right.widthEntranceFringe_m),
    widthExitFringe_m(right.widthExitFringe_m),
55 56 57 58 59 60 61
    reinitialize_m(right.reinitialize_m),
    entranceParameter1_m(right.entranceParameter1_m),
    entranceParameter2_m(right.entranceParameter2_m),
    entranceParameter3_m(right.entranceParameter3_m),
    exitParameter1_m(right.exitParameter1_m),
    exitParameter2_m(right.exitParameter2_m),
    exitParameter3_m(right.exitParameter3_m),
62 63
    engeCoeffsEntry_m(right.engeCoeffsEntry_m),
    engeCoeffsExit_m(right.engeCoeffsExit_m),
64 65 66 67
    entryFieldValues_m(NULL),
    exitFieldValues_m(NULL),
    entryFieldAccel_m(NULL),
    exitFieldAccel_m(NULL),
68 69 70 71 72 73 74 75
    deltaBeginEntry_m(right.deltaBeginEntry_m),
    deltaEndEntry_m(right.deltaEndEntry_m),
    polyOrderEntry_m(right.polyOrderEntry_m),
    deltaBeginExit_m(right.deltaBeginExit_m),
    deltaEndExit_m(right.deltaEndExit_m),
    polyOrderExit_m(right.polyOrderExit_m),
    cosEntranceAngle_m(right.cosEntranceAngle_m),
    sinEntranceAngle_m(right.sinEntranceAngle_m),
76
    tanEntranceAngle_m(right.tanEntranceAngle_m),
77
    tanExitAngle_m(right.tanExitAngle_m),
78 79 80 81 82 83
    beginToEnd_m(right.beginToEnd_m),
    beginToEnd_lcs_m(right.beginToEnd_lcs_m),
    toEntranceRegion_m(right.toEntranceRegion_m),
    toExitRegion_m(right.toExitRegion_m),
    computeAngleTrafo_m(right.computeAngleTrafo_m),
    maxAngle_m(right.maxAngle_m),
kraus's avatar
kraus committed
84
    nSlices_m(right.nSlices_m){
85 86
}

87
Bend2D::Bend2D(const std::string &name):
88
    BendBase(name),
89 90 91 92 93 94 95
    messageHeader_m(" * "),
    pusher_m(),
    designRadius_m(0.0),
    exitAngle_m(0.0),
    fieldIndex_m(0.0),
    startField_m(0.0),
    endField_m(0.0),
96 97
    widthEntranceFringe_m(0.0),
    widthExitFringe_m(0.0),
98 99 100 101 102 103 104
    reinitialize_m(false),
    entranceParameter1_m(0.0),
    entranceParameter2_m(0.0),
    entranceParameter3_m(0.0),
    exitParameter1_m(0.0),
    exitParameter2_m(0.0),
    exitParameter3_m(0.0),
105 106 107 108
    entryFieldValues_m(NULL),
    exitFieldValues_m(NULL),
    entryFieldAccel_m(NULL),
    exitFieldAccel_m(NULL),
109 110 111 112 113 114 115 116
    deltaBeginEntry_m(0.0),
    deltaEndEntry_m(0.0),
    polyOrderEntry_m(0),
    deltaBeginExit_m(0.0),
    deltaEndExit_m(0.0),
    polyOrderExit_m(0),
    cosEntranceAngle_m(1.0),
    sinEntranceAngle_m(0.0),
117
    tanEntranceAngle_m(0.0),
118
    tanExitAngle_m(0.0),
119
    maxAngle_m(0.0),
kraus's avatar
kraus committed
120
    nSlices_m(1){
121 122
}

123
Bend2D::~Bend2D() {
124 125 126 127 128 129 130 131 132 133 134 135 136
    if (entryFieldAccel_m != NULL) {
        for (unsigned int i = 0; i < 3u; ++ i) {
            gsl_spline_free(entryFieldValues_m[i]);
            gsl_spline_free(exitFieldValues_m[i]);
            entryFieldValues_m[i] = NULL;
            exitFieldValues_m[i] = NULL;
        }
        delete[] entryFieldValues_m;
        delete[] exitFieldValues_m;

        gsl_interp_accel_free(entryFieldAccel_m);
        gsl_interp_accel_free(exitFieldAccel_m);
    }
137 138 139 140 141 142 143 144 145 146
}
/*
 * OPAL-T Methods.
 * ===============
 */

/*
 *  This function merely repackages the field arrays as type Vector_t and calls
 *  the equivalent method but with the Vector_t data types.
 */
147
bool Bend2D::apply(const size_t &i,
148 149 150
                 const double &t,
                 Vector_t &E,
                 Vector_t &B) {
151 152
    Vector_t P;
    return apply(RefPartBunch_m->R[i], P, t, E, B);
153 154
}

155
bool Bend2D::apply(const Vector_t &R,
156 157 158 159
                   const Vector_t &/*P*/,
                   const double &/*t*/,
                   Vector_t &/*E*/,
                   Vector_t &B) {
160

161
    if (designRadius_m > 0.0) {
162

163 164 165
        Vector_t X = R;
        Vector_t bField(0.0);
        if (!calculateMapField(X, bField)) {
166

167
            B += fieldAmplitude_m * bField;
168

169
            return false;
170
        }
171 172

        return true;
173 174 175 176 177
    }
    return false;

}

178
bool Bend2D::applyToReferenceParticle(const Vector_t &R,
179 180 181 182
                                    const Vector_t &P,
                                    const double &t,
                                    Vector_t &E,
                                    Vector_t &B) {
183
    return apply(R, P, t, E, B);
184 185
}

186
void Bend2D::goOnline(const double &) {
187
    online_m = true;
188 189
}

190
void Bend2D::initialise(PartBunchBase<double, 3> *bunch,
191 192
                      double &startField,
                      double &endField) {
193 194 195

    Inform msg(messageHeader_m.c_str(), *gmsg);

196
    if (initializeFieldMap()) {
197 198 199 200 201
        std::string name = Util::toUpper(getName()) + " ";
        msg << level2
            << "======================================================================\n"
            << "===== " << std::left << std::setw(64) << std::setfill('=') << name << "\n"
            << "======================================================================\n";
202

203 204
        setupPusher(bunch);
        readFieldMap(msg);
205
        setupBendGeometry(startField, endField);
206 207 208

        double bendAngleX = 0.0;
        double bendAngleY = 0.0;
209 210
        calculateRefTrajectory(bendAngleX, bendAngleY);
        print(msg, bendAngleX, bendAngleY);
211 212 213 214 215

        // Pass start and end of field to calling function.
        startField = startField_m;
        endField = endField_m;

216 217 218 219
        startField_m -= elementEdge_m;
        endField_m -= elementEdge_m;
        elementEdge_m = 0.0;

220
        setupFringeWidths();
221 222 223 224
        msg << level2
            << "======================================================================\n"
            << "======================================================================\n"
            << endl;
225 226 227 228 229 230 231 232
    } else {
        ERRORMSG("There is something wrong with your field map \""
                 << fileName_m
                 << "\"");
        endField = startField - 1.0e-3;
    }
}

233
void Bend2D::adjustFringeFields(double ratio) {
234
    findChordLength(chordLength_m);
235 236 237 238 239

    double delta = std::abs(entranceParameter1_m - entranceParameter2_m);
    entranceParameter1_m = entranceParameter2_m - delta * ratio;

    delta = std::abs(entranceParameter2_m - entranceParameter3_m);
240
    entranceParameter3_m = entranceParameter2_m + delta * ratio;
241 242

    delta = std::abs(exitParameter1_m - exitParameter2_m);
243
    exitParameter1_m = exitParameter2_m - delta * ratio;
244 245 246 247

    delta = std::abs(exitParameter2_m - exitParameter3_m);
    exitParameter3_m = exitParameter2_m + delta * ratio;

248
    setupFringeWidths();
249 250
}

251
double Bend2D::calculateBendAngle() {
252

253
    const double betaGamma = calcBetaGamma();
254
    // const double beta = betaGamma / gamma;
255
    const double deltaT = RefPartBunch_m->getdT();
256
    const double cdt = Physics::c * deltaT;
257
    // Integrate through field for initial angle.
258
    Vector_t oldX;
259 260
    Vector_t X = -deltaBeginEntry_m * Vector_t(tanEntranceAngle_m, 0.0, 1.0);
    Vector_t P = betaGamma * Vector_t(sinEntranceAngle_m, 0.0, cosEntranceAngle_m);
261 262
    double deltaS = 0.0;
    double bendLength = endField_m - startField_m;
263
    const Vector_t eField(0.0);
264

265
    while (deltaS < bendLength) {
266 267
        oldX = X;
        X /= cdt;
268
        pusher_m.push(X, P, deltaT);
269
        X *= cdt;
270 271

        Vector_t bField(0.0, 0.0, 0.0);
272
        calculateMapField(X, bField);
273 274
        bField = fieldAmplitude_m * bField;

275
        X /= cdt;
276 277 278
        pusher_m.kick(X, P, eField, bField, deltaT);

        pusher_m.push(X, P, deltaT);
279
        X *= cdt;
280

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
281
        deltaS += euclidean_norm(X - oldX);
282 283 284

    }

285
    double angle = - std::atan2(P(0), P(2)) + entranceAngle_m;
286 287 288 289

    return angle;
}

290
void Bend2D::calcEngeFunction(double zNormalized,
291 292 293 294 295 296 297 298 299 300
                             const std::vector<double> &engeCoeff,
                             int polyOrder,
                             double &engeFunc,
                             double &engeFuncDeriv,
                             double &engeFuncSecDerivNorm) {

    double expSum = 0.0;
    double expSumDeriv = 0.0;
    double expSumSecDeriv = 0.0;

301
    if (polyOrder >= 2) {
302 303 304 305 306

        expSum = engeCoeff.at(0)
            + engeCoeff.at(1) * zNormalized;
        expSumDeriv = engeCoeff.at(1);

307
        for (int index = 2; index <= polyOrder; index++) {
308
            expSum += engeCoeff.at(index) * std::pow(zNormalized, index);
309
            expSumDeriv += index * engeCoeff.at(index)
310
                * std::pow(zNormalized, index - 1);
311
            expSumSecDeriv += index * (index - 1) * engeCoeff.at(index)
312
                * std::pow(zNormalized, index - 2);
313 314
        }

315
    } else if (polyOrder == 1) {
316 317 318 319 320 321 322 323

        expSum = engeCoeff.at(0)
            + engeCoeff.at(1) * zNormalized;
        expSumDeriv = engeCoeff.at(1);

    } else
        expSum = engeCoeff.at(0);

324
    double engeExp = std::exp(expSum);
325 326
    engeFunc = 1.0 / (1.0 + engeExp);

327
    if (!std::isnan(engeFunc)) {
328 329

        expSumDeriv /= gap_m;
330
        expSumSecDeriv /= std::pow(gap_m, 2.0);
331
        double engeExpDeriv = expSumDeriv * engeExp;
332
        double engeExpSecDeriv = (expSumSecDeriv + std::pow(expSumDeriv, 2.0))
333
            * engeExp;
334
        double engeFuncSq = std::pow(engeFunc, 2.0);
335 336 337 338 339 340

        engeFuncDeriv = -engeExpDeriv * engeFuncSq;
        if (std::isnan(engeFuncDeriv) || std::isinf(engeFuncDeriv))
            engeFuncDeriv = 0.0;

        engeFuncSecDerivNorm = -engeExpSecDeriv * engeFunc
341
            + 2.0 * std::pow(engeExpDeriv, 2.0)
342 343 344 345 346 347 348 349 350 351 352 353
            * engeFuncSq;
        if (std::isnan(engeFuncSecDerivNorm) || std::isinf(engeFuncSecDerivNorm))
            engeFuncSecDerivNorm = 0.0;

    } else {
        engeFunc = 0.0;
        engeFuncDeriv = 0.0;
        engeFuncSecDerivNorm = 0.0;

    }
}

354 355
// :FIXME: is this correct?
Vector_t Bend2D::calcCentralField(const Vector_t &/*R*/, double /*deltaX*/) {
356

357
    Vector_t B(0, 0, 0);
ext-rogers_c's avatar
ext-rogers_c committed
358 359 360 361
    //double nOverRho = fieldIndex_m / designRadius_m;
    //double expFactor = exp(-nOverRho * deltaX);
    //double bxBzFactor = expFactor * nOverRho * R(1);
    //Vector_t rotationCenter(-designRadius_m, R(1), 0.0);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
362
    //double cosangle = dot(R - rotationCenter, Vector_t(1, 0, 0)) / euclidean_norm(R - rotationCenter);
363

ext-rogers_c's avatar
ext-rogers_c committed
364
    //B(0) = -bxBzFactor * cosangle;
365 366
    //B(1) = expFactor * (1.0 - std::pow(nOverRho * R(1), 2.0) / 2.0);
    //B(2) = -bxBzFactor * std::sqrt(1 - std::pow(cosangle, 2));
367

368
    B(1) = 1.0;
369

370 371
    return B;
}
372

373
Vector_t Bend2D::calcEntranceFringeField(const Vector_t &R,
374
                                         double /*deltaX*/) {
375

376 377 378
    const CoordinateSystemTrafo toEntranceRegion(Vector_t(0, 0, entranceParameter2_m),
                                                 Quaternion(0, 0, 1, 0));
    const Vector_t Rprime = toEntranceRegion.transformTo(R);
379

380
    Vector_t B(0.0);
381

382 383 384 385 386 387
    if (Rprime(2) <= -entranceParameter3_m) {
        B(1) = 1;
    } else if (Rprime(2) < -entranceParameter1_m) {
        double engeFunc = gsl_spline_eval(entryFieldValues_m[0], Rprime(2), entryFieldAccel_m);
        double engeFuncDeriv = gsl_spline_eval(entryFieldValues_m[1], Rprime(2), entryFieldAccel_m);
        double engeFuncSecDeriv = gsl_spline_eval(entryFieldValues_m[2], Rprime(2), entryFieldAccel_m);
388

ext-rogers_c's avatar
ext-rogers_c committed
389 390 391
        //double nOverRho = fieldIndex_m / designRadius_m;
        //double expFactor = exp(-nOverRho * deltaX);
        //double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;
392

ext-rogers_c's avatar
ext-rogers_c committed
393 394 395 396
        //double bXEntrance = -nOverRho * expFactor * Rprime(1) * engeFunc;
        //double bYEntrance = (expFactor * engeFunc *
        //                      (1.0  - 0.5 * trigFactor * pow(Rprime(1), 2.0)));
        //double bZEntrance = expFactor * Rprime(1) * engeFuncDeriv;
397 398


399
        // B(1) = (engeFunc *
kraus's avatar
kraus committed
400
        //  (1.0 - 0.5 * engeFuncSecDerivNorm * pow(Rprime(1), 2.0)));
ext-rogers_c's avatar
ext-rogers_c committed
401

402
        B(1) = (engeFunc - 0.5 * engeFuncSecDeriv * std::pow(Rprime(1), 2.0));
ext-rogers_c's avatar
ext-rogers_c committed
403

kraus's avatar
kraus committed
404
        B(2) = engeFuncDeriv * Rprime(1);
405
    }
406

407 408
    return toEntranceRegion.rotateFrom(B);
}
409

410
Vector_t Bend2D::calcExitFringeField(const Vector_t &R,
411
                                     double /*deltaX*/) {
412

413 414 415
    const CoordinateSystemTrafo fromEndToExitRegion(Vector_t(0, 0, exitParameter2_m),
                                                    Quaternion(1, 0, 0, 0));
    const CoordinateSystemTrafo toExitRegion = (fromEndToExitRegion *
416
                                                 getBeginToEnd_local());
417
    const Vector_t Rprime = toExitRegion.transformTo(R);
418

419
    Vector_t B(0.0);
420

421 422 423 424 425 426
    if (Rprime(2) <= exitParameter1_m) {
        B(1) = 1;
    } else if (Rprime(2) < exitParameter3_m) {
        double engeFunc = gsl_spline_eval(exitFieldValues_m[0], Rprime(2), exitFieldAccel_m);
        double engeFuncDeriv = gsl_spline_eval(exitFieldValues_m[1], Rprime(2), exitFieldAccel_m);
        double engeFuncSecDeriv = gsl_spline_eval(exitFieldValues_m[2], Rprime(2), exitFieldAccel_m);
427

ext-rogers_c's avatar
ext-rogers_c committed
428 429 430 431 432 433 434 435 436 437 438
        //double nOverRho = fieldIndex_m / designRadius_m;
        //double expFactor = exp(-nOverRho * deltaX);
        //double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;

        //double bXExit = -nOverRho * expFactor * Rprime(1) * engeFunc;
        //double bYExit = (expFactor * engeFunc *
        //                 (1.0 - 0.5 * trigFactor * pow(Rprime(1), 2.0)));
        //double bZExit = expFactor * Rprime(1) * engeFuncDeriv;

        //B(1) = (engeFunc *
        //        (1.0 - 0.5 * engeFuncSecDerivNorm * pow(Rprime(1), 2.0)));
439
        B(1) = (engeFunc - 0.5 * engeFuncSecDeriv * std::pow(Rprime(1), 2.0));
440

441
        B(2) = engeFuncDeriv * Rprime(1);
442
    }
443
    return toExitRegion.rotateFrom(B);
444 445
}

446
bool Bend2D::calculateMapField(const Vector_t &R, Vector_t &B) {
447

448 449 450 451
    B = Vector_t(0.0);
    bool verticallyInside = (std::abs(R(1)) < 0.5 * gap_m);
    bool horizontallyInside = false;
    Vector_t rotationCenter(-designRadius_m * cosEntranceAngle_m, R(1), designRadius_m * sinEntranceAngle_m);
452
    if (inMagnetCentralRegion(R)) {
453
        if (verticallyInside) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
454
            double deltaX = 0.0;//euclidean_norm(R - rotationCenter) - designRadius_m;
455 456 457 458 459 460 461 462 463 464 465 466 467 468
            bool inEntranceRegion = isPositionInEntranceField(R);
            bool inExitRegion = isPositionInExitField(R);

            if (inEntranceRegion && inExitRegion) {
                Vector_t Rp = transformToEntranceRegion(R);
                Vector_t Rpp = transformToExitRegion(R);

                if (std::abs(Rp[2]) < std::abs(Rpp[2])) {
                    inExitRegion = false;
                } else {
                    inEntranceRegion = false;
                }
            }
            if (inEntranceRegion) {
469
                B = calcEntranceFringeField(R, deltaX);
470
            } else if (inExitRegion) {
471 472 473 474
                B = calcExitFringeField(R, deltaX);
            } else {
                B = calcCentralField(R, deltaX);
            }
475

476 477 478 479
            return false;
        }
        return true;
    }
480

481 482
    Vector_t BEntrance(0.0), BExit(0.0);
    verticallyInside = (std::abs(R(1)) < gap_m);
483

484 485 486 487
    bool inEntranceRegion = inMagnetEntranceRegion(R);
    bool inExitRegion = inMagnetExitRegion(R);

    if (inEntranceRegion) {
488 489 490 491 492
        horizontallyInside = true;
        if (verticallyInside) {
            BEntrance = calcEntranceFringeField(R, R(0));
        }
    }
493

494
    if (inExitRegion) {
495 496 497
        horizontallyInside = true;
        if (verticallyInside) {
            Vector_t Rprime = getBeginToEnd_local().transformTo(R);
498

499 500 501
            BExit = calcExitFringeField(R, Rprime(0));
        }
    }
502

503
    B = BEntrance + BExit;
504

505 506 507
    bool hitMaterial = (horizontallyInside && (!verticallyInside));
    return hitMaterial;
}
508

509
void Bend2D::calculateRefTrajectory(double &angleX, double &/*angleY*/) {
510

511 512
    std::ofstream trajectoryOutput;
    if (Options::writeBendTrajectories && Ippl::myNode() == 0) {
513 514 515 516 517
        std::string fname = Util::combineFilePath({
            OpalData::getInstance()->getAuxiliaryOutputDirectory(),
            OpalData::getInstance()->getInputBasename() + "_" + getName() + "_traj.dat"
        });
        trajectoryOutput.open(fname);
518 519 520 521 522 523
        trajectoryOutput.precision(12);
        trajectoryOutput << "# " << std::setw(18) << "s"
                         << std::setw(20) << "x"
                         << std::setw(20) << "z"
                         << std::setw(20) << "By"
                         << std::endl;
524 525
    }

526 527 528 529
    const double gamma = calcGamma();
    const double betaGamma = calcBetaGamma();
    Vector_t X = -deltaBeginEntry_m * Vector_t(tanEntranceAngle_m, 0.0, 1.0);
    Vector_t P = betaGamma * Vector_t(sinEntranceAngle_m, 0.0, cosEntranceAngle_m);
kraus's avatar
kraus committed
530

531
    if (!refTrajMap_m.empty())
532
        refTrajMap_m.clear();
kraus's avatar
kraus committed
533

534
    refTrajMap_m.push_back(X);
kraus's avatar
kraus committed
535

536
    const Vector_t eField(0.0);
537 538
    const double dt = RefPartBunch_m->getdT();
    const Vector_t scaleFactor(Physics::c * dt);
539 540 541
    const double stepSize = betaGamma / gamma * Physics::c * dt;
    const double bendLength = endField_m - startField_m;
    double deltaS = 0.0;
542
    while (deltaS < bendLength) {
kraus's avatar
kraus committed
543

544
        X /= scaleFactor;
kraus's avatar
kraus committed
545
        pusher_m.push(X, P, dt);
546
        X *= scaleFactor;
kraus's avatar
kraus committed
547 548

        Vector_t bField(0.0, 0.0, 0.0);
549 550 551 552
        Vector_t XInBendFrame = X;

        calculateMapField(XInBendFrame, bField);
        bField = fieldAmplitude_m * bField;
kraus's avatar
kraus committed
553

554 555 556 557 558 559 560
        if (Options::writeBendTrajectories && Ippl::myNode() == 0) {
            trajectoryOutput << std::setw(20) << deltaS + 0.5 * stepSize
                             << std::setw(20) << X(0)
                             << std::setw(20) << X(2)
                             << std::setw(20) << bField(1)
                             << std::endl;
        }
kraus's avatar
kraus committed
561

562
        X /= scaleFactor;
kraus's avatar
kraus committed
563 564 565
        pusher_m.kick(X, P, eField, bField, dt);

        pusher_m.push(X, P, dt);
566
        X *= scaleFactor;
kraus's avatar
kraus committed
567

568
        refTrajMap_m.push_back(X);
kraus's avatar
kraus committed
569

570
        deltaS += stepSize;
kraus's avatar
kraus committed
571 572
    }

573
    angleX = -std::atan2(P(0), P(2)) + entranceAngle_m;
574 575
}

576
double Bend2D::estimateFieldAdjustmentStep(double actualBendAngle) {
577 578 579

    // Estimate field adjustment step.
    double effectiveLength = angle_m * designRadius_m;
580 581 582 583 584
    double effectiveFieldAmplitude = calcFieldAmplitude(2.0 * effectiveLength);
    double ratioSquared = std::pow(fieldAmplitude_m / effectiveFieldAmplitude, 2.0);
    double fieldStep = (angle_m - actualBendAngle) * effectiveFieldAmplitude;
    if (ratioSquared < 1.0) {
        fieldStep = (angle_m - actualBendAngle) * effectiveFieldAmplitude * std::sqrt(1.0 - ratioSquared);
585 586
    }
    fieldStep *= fieldAmplitude_m / std::abs(fieldAmplitude_m);
587 588 589 590 591

    return fieldStep;

}

592
void Bend2D::findBendEffectiveLength(double startField, double endField) {
593 594 595 596 597 598

    /*
     * Use an iterative procedure to set the width of the
     * default field map for the defined field amplitude
     * and bend angle.
     */
599 600 601
    setEngeOriginDelta(0.0);
    setFieldCalcParam();
    setFieldBoundaries(startField, endField);
602

603
    double actualBendAngle = calculateBendAngle();
604 605
    double error = std::abs(actualBendAngle - angle_m);

606
    if (error > 1.0e-6) {
607 608

        double deltaStep = 0.0;
609
        if (std::abs(actualBendAngle) < std::abs(angle_m))
610 611 612 613 614 615 616 617
            deltaStep = -gap_m / 2.0;
        else
            deltaStep = gap_m / 2.0;

        double delta1 = 0.0;
        double bendAngle1 = actualBendAngle;

        double delta2 = deltaStep;
618 619 620 621
        setEngeOriginDelta(delta2);
        setFieldCalcParam();
        setFieldBoundaries(startField, endField);
        double bendAngle2 = calculateBendAngle();
622

623 624
        if (std::abs(bendAngle1) > std::abs(angle_m)) {
            while (std::abs(bendAngle2) > std::abs(angle_m)) {
625
                delta2 += deltaStep;
626 627 628 629
                setEngeOriginDelta(delta2);
                setFieldCalcParam();
                setFieldBoundaries(startField, endField);
                bendAngle2 = calculateBendAngle();
630 631
            }
        } else {
632
            while (std::abs(bendAngle2) < std::abs(angle_m)) {
633
                delta2 += deltaStep;
634 635 636 637
                setEngeOriginDelta(delta2);
                setFieldCalcParam();
                setFieldBoundaries(startField, endField);
                bendAngle2 = calculateBendAngle();
638 639 640 641 642 643
            }
        }

        // Now we should have the proper field map width bracketed.
        unsigned int iterations = 1;
        error = std::abs(actualBendAngle - angle_m);
644
        while (error > 1.0e-6 && iterations < 100) {
645

snuverink_j's avatar
snuverink_j committed
646
            double delta = (delta1 + delta2) / 2.0;
647 648 649 650
            setEngeOriginDelta(delta);
            setFieldCalcParam();
            setFieldBoundaries(startField, endField);
            double newBendAngle = calculateBendAngle();
651 652 653

            error = std::abs(newBendAngle - angle_m);

654
            if (error > 1.0e-6) {
655

656
                if (bendAngle1 - angle_m < 0.0) {
657

658
                    if (newBendAngle - angle_m < 0.0) {
659 660 661
                        bendAngle1 = newBendAngle;
                        delta1 = delta;
                    } else {
snuverink_j's avatar
snuverink_j committed
662
                        // bendAngle2 = newBendAngle;
663 664 665 666 667
                        delta2 = delta;
                    }

                } else {

668
                    if (newBendAngle - angle_m < 0.0) {
snuverink_j's avatar
snuverink_j committed
669
                        // bendAngle2 = newBendAngle;
670 671 672 673 674 675 676 677 678 679 680 681
                        delta2 = delta;
                    } else {
                        bendAngle1 = newBendAngle;
                        delta1 = delta;
                    }
                }
            }
            iterations++;
        }
    }
}

682
void Bend2D::findBendStrength() {
683 684 685 686 687

    /*
     * Use an iterative procedure to set the magnet field amplitude
     * for the defined bend angle.
     */
688 689

    const double tolerance = 1e-7; // [deg]
690 691 692 693 694
    double actualBendAngle = calculateBendAngle();
    double error = std::abs(actualBendAngle - angle_m) * Physics::rad2deg;
    if (error < tolerance)
        return;

695
    double fieldStep = std::copysign(1.0, fieldAmplitude_m) * std::abs(estimateFieldAdjustmentStep(actualBendAngle));
696
    double amplitude1 = fieldAmplitude_m;
697
    double amplitude2 = amplitude1;
698
    double bendAngle1 = actualBendAngle;
699 700 701 702 703 704 705 706 707 708 709 710
    double bendAngle2 = bendAngle1;

    int stepSign = std::abs(bendAngle1) > std::abs(angle_m) ? -1: 1;
    while (true) {
        amplitude1 = amplitude2;
        bendAngle1 = bendAngle2;

        amplitude2 += stepSign * fieldStep;
        fieldAmplitude_m = amplitude2;
        bendAngle2 = calculateBendAngle();
        if (stepSign * (std::abs(bendAngle2) - std::abs(angle_m)) > 0) {
            break;
711 712 713 714 715
        }
    }

    // Now we should have the proper field amplitude bracketed.
    unsigned int iterations = 1;
716
    while (error > tolerance && iterations < 100) {
717 718

        fieldAmplitude_m = (amplitude1 + amplitude2) / 2.0;
719
        double newBendAngle = calculateBendAngle();
720

721
        error = std::abs(newBendAngle - angle_m) * Physics::rad2deg;
722

723
        if (error > tolerance) {
724

725
            if (bendAngle1 - angle_m < 0.0) {
726

727
                if (newBendAngle - angle_m < 0.0) {
728 729 730
                    bendAngle1 = newBendAngle;
                    amplitude1 = fieldAmplitude_m;
                } else {
snuverink_j's avatar
snuverink_j committed
731
                    // bendAngle2 = newBendAngle;
732 733 734 735 736
                    amplitude2 = fieldAmplitude_m;
                }

            } else {

737
                if (newBendAngle - angle_m < 0.0) {
snuverink_j's avatar
snuverink_j committed
738
                    // bendAngle2 = newBendAngle;
739 740 741 742 743 744 745 746 747 748 749
                    amplitude2 = fieldAmplitude_m;
                } else {
                    bendAngle1 = newBendAngle;
                    amplitude1 = fieldAmplitude_m;
                }
            }
        }
        iterations++;
    }
}

750
bool Bend2D::findIdealBendParameters(double chordLength) {
751

752
    bool reinitialize = false;
753

754
    if (angle_m != 0.0) {
755

756
        if (angle_m < 0.0) {
757
            // Negative angle is a positive bend rotated 180 degrees.
758
            setEntranceAngle(std::copysign(1, angle_m) * entranceAngle_m);
759
            setExitAngle(std::copysign(1, angle_m) * exitAngle_m);
760
            angle_m = std::abs(angle_m);
761
            rotationZAxis_m += Physics::pi;
762
        }
763 764
        designRadius_m = calcDesignRadius(chordLength, angle_m);
        fieldAmplitude_m = calcFieldAmplitude(designRadius_m);
765
        reinitialize = true;
766 767
    } else {

768 769 770
        rotationZAxis_m += std::atan2(fieldAmplitudeX_m, fieldAmplitudeY_m);
        double refCharge = RefPartBunch_m->getQ();
        if (refCharge < 0.0) {
771
            rotationZAxis_m -= Physics::pi;
772 773
        }

774 775 776
        fieldAmplitude_m = std::copysign(1.0, refCharge) * std::hypot(fieldAmplitudeX_m, fieldAmplitudeY_m);
        designRadius_m = calcDesignRadius(fieldAmplitude_m);
        angle_m = calcBendAngle(chordLength, designRadius_m);
777
        reinitialize = false;
778
    }
779
    return reinitialize;
780 781
}

782
bool Bend2D::initializeFieldMap() {
783 784 785

    fieldmap_m = Fieldmap::getFieldmap(fileName_m, fast_m);

786 787
    if (fieldmap_m != NULL) {
        if (fileName_m != "1DPROFILE1-DEFAULT")
788 789
            return true;
        else
790
            return setupDefaultFieldMap();
791 792 793 794 795 796

    } else
        return false;

}

797
bool Bend2D::inMagnetCentralRegion(const Vector_t &R) const {
798

799
    Vector_t rotationCenter(-designRadius_m * cosEntranceAngle_m, R(1), designRadius_m * sinEntranceAngle_m);
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
800
    double distFromRotCenter = euclidean_norm(R - rotationCenter);
801 802
    Vector_t Rprime = getBeginToEnd_local().transformTo(R);
    Vector_t Rpprime = computeAngleTrafo_m.transformTo(R);
803

804
    double effectiveAngle = std::fmod(Physics::two_pi - std::atan2(Rpprime(0), Rpprime(2)), Physics::two_pi);
805

806 807 808 809 810
    if (std::abs(distFromRotCenter - designRadius_m) < 0.5 * aperture_m.second[0] &&
        effectiveAngle >= 0.0 && effectiveAngle < maxAngle_m) {
        if (effectiveAngle < 0.5 * maxAngle_m) return R(2) >= 0.0;
        return Rprime(2) < 0.0;
    }
811

812
    return false;
813 814
}

815
bool Bend2D::inMagnetEntranceRegion(const Vector_t &R) const {
816

817 818 819
    return (R(2) >= entranceParameter1_m &&
            R(2) < 0.0 &&
            std::abs(R(0)) < aperture_m.second[0]);
820 821
}

822
bool Bend2D::inMagnetExitRegion(const Vector_t &R) const {
823

824
    Vector_t Rprime = getBeginToEnd_local().transformTo(R);
825

826 827 828
    return (Rprime(2) >= 0 &&
            Rprime(2) < exitParameter3_m &&
            std::abs(Rprime(0)) < aperture_m.second[0]);
829 830
}

831
bool Bend2D::isPositionInEntranceField(const Vector_t &R) const {
832

833 834 835
    return (polyOrderEntry_m >= 0 &&
            R(2) >= entranceParameter1_m &&
            R(2) < entranceParameter3_m);
836 837
}

838
bool Bend2D::isPositionInExitField(const Vector_t &R) const {
839
    Vector_t Rprime = getBeginToEnd_local().transformTo(R);
840