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

#include "AbsBeamline/SBend.h"
kraus's avatar
kraus committed
22
#include "Algorithms/PartPusher.h"
23
#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
24
#include "AbsBeamline/BeamlineVisitor.h"
25
#include "Utilities/Options.h"
26
#include "Fields/Fieldmap.h"
gsell's avatar
gsell committed
27 28 29 30 31 32 33 34
#include <iostream>
#include <fstream>

// Class SBend
// ------------------------------------------------------------------------

SBend::SBend():
    Component(),
35
    pusher_m(),
36
    fileName_m(""),
gsell's avatar
gsell committed
37
    fieldmap_m(NULL),
kraus's avatar
kraus committed
38
    fast_m(false),
39 40 41 42 43
    angle_m(0.0),
    aperture_m(0.0),
    designEnergy_m(0.0),
    designRadius_m(0.0),
    fieldAmplitude_m(0.0),
44 45
    bX_m(0.0),
    bY_m(0.0),
46 47 48
    angleGreaterThanPi_m(false),
    entranceAngle_m(0.0),
    exitAngle_m(0.0),
49
    fieldIndex_m(0.0),
50
    elementEdge_m(0.0),
gsell's avatar
gsell committed
51
    startField_m(0.0),
kraus's avatar
kraus committed
52
    endField_m(0.0),
53 54
    reinitialize_m(false),
    recalcRefTraj_m(false),
gsell's avatar
gsell committed
55 56
    length_m(0.0),
    gap_m(0.0),
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
    refTrajMapSize_m(0),
    refTrajMapStepSize_m(0.0),
    entranceParameter1_m(0.0),
    entranceParameter2_m(0.0),
    entranceParameter3_m(0.0),
    exitParameter1_m(0.0),
    exitParameter2_m(0.0),
    exitParameter3_m(0.0),
    xOriginEngeEntry_m(0.0),
    zOriginEngeEntry_m(0.0),
    deltaBeginEntry_m(0.0),
    deltaEndEntry_m(0.0),
    polyOrderEntry_m(0),
    xExit_m(0.0),
    zExit_m(0.0),
    xOriginEngeExit_m(0.0),
    zOriginEngeExit_m(0.0),
    deltaBeginExit_m(0.0),
    deltaEndExit_m(0.0),
    polyOrderExit_m(0),
    cosEntranceAngle_m(1.0),
    sinEntranceAngle_m(0.0),
    exitEdgeAngle_m(0.0),
    cosExitAngle_m(1.0),
    sinExitAngle_m(0.0) {
gsell's avatar
gsell committed
82 83 84 85 86
    setElType(isDipole);
}

SBend::SBend(const SBend &right):
    Component(right),
87
    pusher_m(right.pusher_m),
88
    fileName_m(right.fileName_m),
gsell's avatar
gsell committed
89
    fieldmap_m(right.fieldmap_m),
kraus's avatar
kraus committed
90
    fast_m(right.fast_m),
91 92 93 94 95
    angle_m(right.angle_m),
    aperture_m(right.aperture_m),
    designEnergy_m(right.designEnergy_m),
    designRadius_m(right.designRadius_m),
    fieldAmplitude_m(right.fieldAmplitude_m),
96 97
    bX_m(right.bX_m),
    bY_m(right.bY_m),
98 99 100
    angleGreaterThanPi_m(right.angleGreaterThanPi_m),
    entranceAngle_m(right.entranceAngle_m),
    exitAngle_m(right.exitAngle_m),
101
    fieldIndex_m(right.fieldIndex_m),
102
    elementEdge_m(right.elementEdge_m),
gsell's avatar
gsell committed
103 104
    startField_m(right.startField_m),
    endField_m(right.endField_m),
105 106
    reinitialize_m(right.reinitialize_m),
    recalcRefTraj_m(right.recalcRefTraj_m),
gsell's avatar
gsell committed
107 108
    length_m(right.length_m),
    gap_m(right.gap_m),
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
    refTrajMapX_m(right.refTrajMapX_m),
    refTrajMapY_m(right.refTrajMapY_m),
    refTrajMapZ_m(right.refTrajMapZ_m),
    refTrajMapSize_m(right.refTrajMapSize_m),
    refTrajMapStepSize_m(right.refTrajMapStepSize_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),
    xOriginEngeEntry_m(right.xOriginEngeEntry_m),
    zOriginEngeEntry_m(right.zOriginEngeEntry_m),
    deltaBeginEntry_m(right.deltaBeginEntry_m),
    deltaEndEntry_m(right.deltaEndEntry_m),
    polyOrderEntry_m(right.polyOrderEntry_m),
    xExit_m(right.xExit_m),
    zExit_m(right.zExit_m),
    xOriginEngeExit_m(right.xOriginEngeExit_m),
    zOriginEngeExit_m(right.zOriginEngeExit_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),
    exitEdgeAngle_m(right.exitEdgeAngle_m),
    cosExitAngle_m(right.cosExitAngle_m),
    sinExitAngle_m(right.sinExitAngle_m) {

gsell's avatar
gsell committed
138 139
    setElType(isDipole);

140
}
gsell's avatar
gsell committed
141

Steve Russell's avatar
Steve Russell committed
142
SBend::SBend(const std::string &name):
143
    Component(name),
144
    pusher_m(),
145
    fileName_m(""),
gsell's avatar
gsell committed
146
    fieldmap_m(NULL),
kraus's avatar
kraus committed
147
    fast_m(false),
148 149 150 151 152
    angle_m(0.0),
    aperture_m(0.0),
    designEnergy_m(0.0),
    designRadius_m(0.0),
    fieldAmplitude_m(0.0),
153 154
    bX_m(0.0),
    bY_m(0.0),
155 156 157
    angleGreaterThanPi_m(false),
    entranceAngle_m(0.0),
    exitAngle_m(0.0),
158
    fieldIndex_m(0.0),
159
    elementEdge_m(0.0),
gsell's avatar
gsell committed
160 161
    startField_m(0.0),
    endField_m(0.0),
162 163
    reinitialize_m(false),
    recalcRefTraj_m(false),
gsell's avatar
gsell committed
164 165
    length_m(0.0),
    gap_m(0.0),
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
    refTrajMapSize_m(0),
    refTrajMapStepSize_m(0.0),
    entranceParameter1_m(0.0),
    entranceParameter2_m(0.0),
    entranceParameter3_m(0.0),
    exitParameter1_m(0.0),
    exitParameter2_m(0.0),
    exitParameter3_m(0.0),
    xOriginEngeEntry_m(0.0),
    zOriginEngeEntry_m(0.0),
    deltaBeginEntry_m(0.0),
    deltaEndEntry_m(0.0),
    polyOrderEntry_m(0),
    xExit_m(0.0),
    zExit_m(0.0),
    xOriginEngeExit_m(0.0),
    zOriginEngeExit_m(0.0),
    deltaBeginExit_m(0.0),
    deltaEndExit_m(0.0),
    polyOrderExit_m(0),
    cosEntranceAngle_m(1.0),
    sinEntranceAngle_m(0.0),
    exitEdgeAngle_m(0.0),
    cosExitAngle_m(1.0),
    sinExitAngle_m(0.0) {

gsell's avatar
gsell committed
192
    setElType(isDipole);
193

gsell's avatar
gsell committed
194 195 196 197 198 199 200 201 202
}

SBend::~SBend() {
}

void SBend::accept(BeamlineVisitor &visitor) const {
    visitor.visitSBend(*this);
}

203 204 205 206
/*
 * OPAL-MAP methods
 * ================
 */
gsell's avatar
gsell committed
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
double SBend::getNormalComponent(int n) const {
    return getField().getNormalComponent(n);
}

double SBend::getSkewComponent(int n) const {
    return getField().getSkewComponent(n);
}

void SBend::setNormalComponent(int n, double v) {
    getField().setNormalComponent(n, v);
}

void SBend::setSkewComponent(int n, double v) {
    getField().setSkewComponent(n, v);
}

223

224 225 226 227 228 229 230 231 232
/*
 * 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.
 */
233
bool SBend::apply(const size_t &i, const double &t, double E[], double B[]) {
234 235 236 237 238

    Vector_t Ev(0.0, 0.0, 0.0);
    Vector_t Bv(0.0, 0.0, 0.0);
    if(apply(RefPartBunch_m->R[i], RefPartBunch_m->get_rmean(), t, Ev, Bv))
        return true;
gsell's avatar
gsell committed
239 240 241 242 243 244 245 246 247 248 249

    E[0] = Ev(0);
    E[1] = Ev(1);
    E[2] = Ev(2);
    B[0] = Bv(0);
    B[1] = Bv(1);
    B[2] = Bv(2);

    return false;
}

250
bool SBend::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
gsell's avatar
gsell committed
251

252
    if(designRadius_m > 0.0) {
gsell's avatar
gsell committed
253

254 255 256
        // Shift position to magnet frame.
        Vector_t X = RefPartBunch_m->X[i];
        X(2) += startField_m - elementEdge_m;
gsell's avatar
gsell committed
257

258 259 260 261 262 263
        /*
         * Add in transverse bend displacements. (ds is already
         * accounted for.)
         */
        X(0) -= dx_m;
        X(1) -= dy_m;
gsell's avatar
gsell committed
264

265 266 267 268 269
        // Get field from field map.
        Vector_t eField(0.0, 0.0, 0.0);
        Vector_t bField(0.0, 0.0, 0.0);
        CalculateMapField(X, eField, bField);
        bField *= fieldAmplitude_m;
gsell's avatar
gsell committed
270

271 272 273
        B(0) += bField(0);
        B(1) += bField(1);
        B(2) += bField(2);
gsell's avatar
gsell committed
274

275
    }
gsell's avatar
gsell committed
276

277 278
    return false;
}
gsell's avatar
gsell committed
279

280 281 282 283 284
bool SBend::apply(const Vector_t &R,
                  const Vector_t &centroid,
                  const double &t,
                  Vector_t &E,
                  Vector_t &B) {
gsell's avatar
gsell committed
285

286
    if(designRadius_m > 0.0) {
gsell's avatar
gsell committed
287

288 289
        int index = static_cast<int>
                    (std::floor((R(2) - startField_m) / refTrajMapStepSize_m));
gsell's avatar
gsell committed
290

291
        if(index > 0 && index + 1 < refTrajMapSize_m) {
gsell's avatar
gsell committed
292

293 294 295 296 297 298 299 300
            // Find indices for position in pre-computed central trajectory map.
            double lever = (R(2) - startField_m) / refTrajMapStepSize_m - index;
            double x = (1.0 - lever) * refTrajMapX_m.at(index)
                       + lever * refTrajMapX_m.at(index + 1);
            double y = (1.0 - lever) * refTrajMapY_m.at(index)
                       + lever * refTrajMapY_m.at(index + 1);
            double z = (1.0 - lever) * refTrajMapZ_m.at(index)
                       + lever * refTrajMapZ_m.at(index + 1);
gsell's avatar
gsell committed
301

302 303 304 305 306
            // Adjust position relative to pre-computed central trajectory map.
            Vector_t X(0.0, 0.0, 0.0);
            X(0) = R(0) + x;
            X(1) = R(1) + y;
            X(2) = z;
gsell's avatar
gsell committed
307

308 309 310
            Vector_t tempE(0.0, 0.0, 0.0);
            Vector_t tempB(0.0, 0.0, 0.0);
            Vector_t XInBendFrame = RotateToBendFrame(X);
gsell's avatar
gsell committed
311

312 313 314 315 316 317
            /*
             * Add in transverse bend displacements. (ds is already
             * accounted for.)
             */
            XInBendFrame(0) -= dx_m;
            XInBendFrame(1) -= dy_m;
gsell's avatar
gsell committed
318

319 320
            CalculateMapField(XInBendFrame, tempE, tempB);
            tempB = fieldAmplitude_m * RotateOutOfBendFrame(tempB);
gsell's avatar
gsell committed
321

322 323 324
            B(0) += tempB(0);
            B(1) += tempB(1);
            B(2) += tempB(2);
gsell's avatar
gsell committed
325

326 327
        }
    }
gsell's avatar
gsell committed
328

329
    return false;
gsell's avatar
gsell committed
330

331
}
gsell's avatar
gsell committed
332

333 334 335
bool SBend::bends() const {
    return true;
}
gsell's avatar
gsell committed
336

337 338 339
void SBend::finalise() {
    online_m = false;
}
gsell's avatar
gsell committed
340

kraus's avatar
kraus committed
341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
void SBend::goOnline(const double &) {

    // Check if we need to reinitialize the bend field amplitude.
    if(reinitialize_m) {
        reinitialize_m = Reinitialize();
        recalcRefTraj_m = false;
    }

    /*
     * Always recalculate the reference trajectory on first call even
     * if we do not reinitialize the bend. The reference trajectory
     * has to be calculated at the same energy as the actual beam or
     * we do not get accurate values for the magnetic field in the output
     * file.
     */
    if(recalcRefTraj_m) {
        double angleX = 0.0;
        double angleY = 0.0;
        CalculateRefTrajectory(angleX, angleY);
        recalcRefTraj_m = false;
    }

    online_m = true;
}

366 367 368 369
void SBend::getDimensions(double &sBegin, double &sEnd) const {
    sBegin = startField_m;
    sEnd = endField_m;
}
gsell's avatar
gsell committed
370

371 372
ElementBase::ElementType SBend::getType() const {
    return SBEND;
373
}
gsell's avatar
gsell committed
374

375 376 377 378
void SBend::initialise(PartBunch *bunch,
                       double &startField,
                       double &endField,
                       const double &scaleFactor) {
gsell's avatar
gsell committed
379

380
    Inform msg("SBend ");
gsell's avatar
gsell committed
381

382
    if(InitializeFieldMap(msg)) {
gsell's avatar
gsell committed
383

384 385 386 387 388 389 390 391
        SetupPusher(bunch);
        ReadFieldMap(msg);
        SetupBendGeometry(msg, startField, endField);
        double bendAngleX = 0.0;
        double bendAngleY = 0.0;
        CalculateRefTrajectory(bendAngleX, bendAngleY);
        recalcRefTraj_m = true;
        Print(msg, bendAngleX, bendAngleY);
392 393

        // Pass start and end of field to calling function.
gsell's avatar
gsell committed
394 395 396 397
        startField = startField_m;
        endField = endField_m;

    } else {
398 399 400 401
        msg << " There is something wrong with your field map \""
            << fileName_m
            << "\"";
        endField = startField - 1.0e-3;
gsell's avatar
gsell committed
402 403 404
    }
}

405 406
double SBend::GetBendAngle() const {
    return angle_m;
gsell's avatar
gsell committed
407 408
}

409 410
double SBend::GetBendRadius() const {
    return designRadius_m;
gsell's avatar
gsell committed
411 412
}

413 414
double SBend::GetEffectiveCenter() const {
    return elementEdge_m + designRadius_m * angle_m / 2.0;
415 416
}

417 418
double SBend::GetEffectiveLength() const {
    return designRadius_m * angle_m;
gsell's avatar
gsell committed
419 420
}

421 422
std::string SBend::GetFieldMapFN() const {
    return fileName_m;
423 424
}

425 426
double SBend::GetStartElement() const {
    return elementEdge_m;
gsell's avatar
gsell committed
427 428
}

429 430
void SBend::SetAngleGreaterThanPiFlag(bool angleGreaterThanPi) {
    angleGreaterThanPi_m = angleGreaterThanPi;
gsell's avatar
gsell committed
431 432
}

433 434
void SBend::SetAperture(double aperture) {
    aperture_m = std::abs(aperture);
gsell's avatar
gsell committed
435 436
}

437 438
void SBend::SetBendAngle(double angle) {
    angle_m = angle;
gsell's avatar
gsell committed
439 440
}

441 442
void SBend::SetBeta(double beta) {
    Orientation_m(1) = beta;
gsell's avatar
gsell committed
443 444
}

445 446
void SBend::SetDesignEnergy(double energy) {
    designEnergy_m = std::abs(energy);
gsell's avatar
gsell committed
447 448
}

449 450 451
void SBend::SetEntranceAngle(double entranceAngle) {
    entranceAngle_m = entranceAngle;
}
gsell's avatar
gsell committed
452

453 454
void SBend::setExitAngle(double exitAngle) {
    exitAngle_m = exitAngle;
455
}
gsell's avatar
gsell committed
456

457 458 459
void SBend::SetFieldAmplitude(double k0, double k0s) {
    bY_m = k0;
    bX_m = k0s;
460
}
gsell's avatar
gsell committed
461

462 463
void SBend::SetFieldMapFN(std::string fileName) {
    fileName_m = fileName;
464
}
Steve Russell's avatar
Steve Russell committed
465

466 467
void SBend::SetFullGap(double gap) {
    gap_m = std::abs(gap);
468
}
gsell's avatar
gsell committed
469

470
void SBend::SetK1(double k1) {
471
    fieldIndex_m = k1;
472
}
gsell's avatar
gsell committed
473

474 475
void SBend::SetLength(double length) {
    length_m = std::abs(length);
476
}
Steve Russell's avatar
Steve Russell committed
477

478 479
void SBend::SetRotationAboutZ(double rotation) {
    Orientation_m(2) = rotation;
480
}
Steve Russell's avatar
Steve Russell committed
481

482
void SBend::AdjustFringeFields(double ratio) {
gsell's avatar
gsell committed
483

484 485 486 487 488 489 490 491 492 493 494
    double delta = std::abs(entranceParameter1_m - entranceParameter2_m);
    entranceParameter1_m = entranceParameter2_m - delta * ratio;

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

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

    delta = std::abs(exitParameter2_m - exitParameter3_m);
    exitParameter3_m = exitParameter2_m + delta * ratio;
Steve Russell's avatar
Steve Russell committed
495

gsell's avatar
gsell committed
496 497
}

498
double SBend::CalculateBendAngle() {
gsell's avatar
gsell committed
499 500

    const double mass = RefPartBunch_m->getM();
501
    const double gamma = designEnergy_m / mass + 1.0;
gsell's avatar
gsell committed
502
    const double betaGamma = sqrt(pow(gamma, 2.0) - 1.0);
503
    const double beta = betaGamma / gamma;
gsell's avatar
gsell committed
504 505 506
    const double deltaT = RefPartBunch_m->getdT();

    // Integrate through field for initial angle.
507 508 509 510
    Vector_t X(0.0, 0.0, startField_m - elementEdge_m);
    Vector_t P(0.0, 0.0, betaGamma);
    double deltaS = 0.0;
    double bendLength = endField_m - startField_m;
gsell's avatar
gsell committed
511

512
    while(deltaS < bendLength) {
gsell's avatar
gsell committed
513 514 515 516 517

        X /= Vector_t(Physics::c * deltaT);
        pusher_m.push(X, P, deltaT);
        X *= Vector_t(Physics::c * deltaT);

518 519 520 521
        Vector_t eField(0.0, 0.0, 0.0);
        Vector_t bField(0.0, 0.0, 0.0);
        CalculateMapField(X, eField, bField);
        bField = fieldAmplitude_m * bField;
gsell's avatar
gsell committed
522 523

        X /= Vector_t(Physics::c * deltaT);
524
        pusher_m.kick(X, P, eField, bField, deltaT);
gsell's avatar
gsell committed
525 526 527 528

        pusher_m.push(X, P, deltaT);
        X *= Vector_t(Physics::c * deltaT);

529
        deltaS += deltaT * beta * Physics::c;
gsell's avatar
gsell committed
530

531
    }
gsell's avatar
gsell committed
532

533
    double angle =  -atan2(P(0), P(2));
gsell's avatar
gsell committed
534

535
    return angle;
gsell's avatar
gsell committed
536

537
}
gsell's avatar
gsell committed
538

539 540 541 542
void SBend::CalcCentralField(Vector_t R,
                             double deltaX,
                             double angle,
                             Vector_t &B) {
gsell's avatar
gsell committed
543

544 545
    double nOverRho = fieldIndex_m / designRadius_m;
    double expFactor = exp(-nOverRho * deltaX);
546 547
    double bxBzFactor = expFactor * nOverRho * R(1);

548 549 550
    B(0) = -bxBzFactor * cos(angle);
    B(1) = expFactor * (1.0 - pow(nOverRho * R(1), 2.0) / 2.0);
    B(2) = -bxBzFactor * sin(angle);
gsell's avatar
gsell committed
551

552
}
gsell's avatar
gsell committed
553

554 555 556 557 558
void SBend::CalcEngeFunction(double zNormalized,
                             std::vector<double> engeCoeff,
                             int polyOrder,
                             double &engeFunc,
                             double &engeFuncDeriv,
559
                             double &engeFuncSecDerivNorm) {
560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576

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

    if(polyOrder >= 2) {

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

        for(int index = 2; index <= polyOrder; index++) {
            expSum += engeCoeff.at(index) * pow(zNormalized, index);
            expSumDeriv += index * engeCoeff.at(index)
                           * pow(zNormalized, index - 1);
            expSumSecDeriv += index * (index - 1) * engeCoeff.at(index)
                              * pow(zNormalized, index - 2);
gsell's avatar
gsell committed
577 578
        }

579
    } else if(polyOrder == 1) {
Steve Russell's avatar
Steve Russell committed
580

581 582 583
        expSum = engeCoeff.at(0)
                 + engeCoeff.at(1) * zNormalized;
        expSumDeriv = engeCoeff.at(1);
gsell's avatar
gsell committed
584

585 586 587 588 589
    } else
        expSum = engeCoeff.at(0);

    double engeExp = exp(expSum);
    engeFunc = 1.0 / (1.0 + engeExp);
590

Jianjun Yang's avatar
Jianjun Yang committed
591
    if(!std::isnan(engeFunc)) {
592 593 594 595 596 597 598 599 600

        expSumDeriv /= gap_m;
        expSumSecDeriv /= pow(gap_m, 2.0);
        double engeExpDeriv = expSumDeriv * engeExp;
        double engeExpSecDeriv = (expSumSecDeriv + pow(expSumDeriv, 2.0))
                                 * engeExp;
        double engeFuncSq = pow(engeFunc, 2.0);

        engeFuncDeriv = -engeExpDeriv * engeFuncSq;
Jianjun Yang's avatar
Jianjun Yang committed
601
        if (std::isnan(engeFuncDeriv) || std::isinf(engeFuncDeriv))
602 603 604 605 606
            engeFuncDeriv = 0.0;

        engeFuncSecDerivNorm = -engeExpSecDeriv * engeFunc
                               + 2.0 * pow(engeExpDeriv, 2.0)
                                 * engeFuncSq;
Jianjun Yang's avatar
Jianjun Yang committed
607
        if (std::isnan(engeFuncSecDerivNorm) || std::isinf(engeFuncSecDerivNorm))
608 609
            engeFuncSecDerivNorm = 0.0;

610
    } else {
611 612
        engeFunc = 0.0;
        engeFuncDeriv = 0.0;
613
        engeFuncSecDerivNorm = 0.0;
gsell's avatar
gsell committed
614

615
    }
616
}
gsell's avatar
gsell committed
617

618 619 620
void SBend::CalcEntranceFringeField(Vector_t REntrance,
                                    double deltaX,
                                    Vector_t &B) {
gsell's avatar
gsell committed
621

622 623 624
    double zNormalized = -REntrance(2) / gap_m;
    double engeFunc = 0.0;
    double engeFuncDeriv = 0.0;
625
    double engeFuncSecDerivNorm = 0.0;
gsell's avatar
gsell committed
626

627 628 629 630 631
    CalcEngeFunction(zNormalized,
                     engeCoeffsEntry_m,
                     polyOrderEntry_m,
                     engeFunc,
                     engeFuncDeriv,
632
                     engeFuncSecDerivNorm);
gsell's avatar
gsell committed
633

634 635
    double nOverRho = fieldIndex_m / designRadius_m;
    double expFactor = exp(-nOverRho * deltaX);
636 637
    double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;

638
    double bXEntrance = -engeFunc * nOverRho * expFactor* REntrance(1);
639 640
    double bYEntrance = expFactor * engeFunc
                        * (1.0  - trigFactor * pow(REntrance(1), 2.0) / 2.0);
641
    double bZEntrance = -expFactor * engeFuncDeriv * REntrance(1);
gsell's avatar
gsell committed
642

643 644 645
    B(0) = bXEntrance * cosEntranceAngle_m - bZEntrance * sinEntranceAngle_m;
    B(1) = bYEntrance;
    B(2) = bXEntrance * sinEntranceAngle_m + bZEntrance * cosEntranceAngle_m;
gsell's avatar
gsell committed
646 647
}

648 649 650 651 652
void SBend::CalcExitFringeField(Vector_t RExit, double deltaX, Vector_t &B) {

    double zNormalized = RExit(2) / gap_m;
    double engeFunc = 0.0;
    double engeFuncDeriv = 0.0;
653
    double engeFuncSecDerivNorm = 0.0;
654 655 656 657 658
    CalcEngeFunction(zNormalized,
                     engeCoeffsExit_m,
                     polyOrderExit_m,
                     engeFunc,
                     engeFuncDeriv,
659
                     engeFuncSecDerivNorm);
660

661 662
    double nOverRho = fieldIndex_m / designRadius_m;
    double expFactor = exp(-nOverRho * deltaX);
663 664
    double trigFactor = pow(nOverRho, 2.0) + engeFuncSecDerivNorm;

665
    double bXExit = -engeFunc * nOverRho * expFactor* RExit(1);
666 667
    double bYExit = expFactor * engeFunc
                    * (1.0 - trigFactor * pow(RExit(1), 2.0) / 2.0);
668
    double bZExit = expFactor * engeFuncDeriv * RExit(1);
669 670 671 672 673

    B(0) = bXExit * cosExitAngle_m - bZExit * sinExitAngle_m;
    B(1) = bYExit;
    B(2) = bXExit * sinExitAngle_m + bZExit * cosExitAngle_m;
}
gsell's avatar
gsell committed
674

675
void SBend::CalculateMapField(Vector_t R, Vector_t &E, Vector_t &B) {
gsell's avatar
gsell committed
676

677 678
    E = Vector_t(0.0);
    B = Vector_t(0.0);
gsell's avatar
gsell committed
679

680 681 682 683 684 685 686 687 688
    //    Vector_t REntrance(0.0, 0.0, 0.0);
    //    Vector_t RExit(0.0, 0.0, 0.0);
    //    if (IsPositionInEntranceField(R, REntrance)) {
    //        CalcEntranceFringeField(REntrance, 0.0, B);
    //    } else if (IsPositionInExitField(R, RExit)) {
    //        CalcExitFringeField(RExit, 0.0, B);
    //    } else {
    //        CalcCentralField(R, 0.0, 0.0, B);
    //    }
gsell's avatar
gsell committed
689

690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707
    double deltaXEntrance = 0.0;
    double deltaXExit = 0.0;
    bool inEntranceRegion = InMagnetEntranceRegion(R, deltaXEntrance);
    bool inExitRegion = InMagnetExitRegion(R, deltaXExit);

    if(!inEntranceRegion && !inExitRegion) {

        double deltaX = 0.0;
        double angle = 0.0;
        if(InMagnetCentralRegion(R, deltaX, angle)) {
            Vector_t REntrance(0.0, 0.0, 0.0);
            Vector_t RExit(0.0, 0.0, 0.0);
            if(IsPositionInEntranceField(R, REntrance))
                CalcEntranceFringeField(REntrance, deltaX, B);
            else if(IsPositionInExitField(R, RExit))
                CalcExitFringeField(RExit, deltaX, B);
            else
                CalcCentralField(R, deltaX, angle, B);
gsell's avatar
gsell committed
708

709
        }
gsell's avatar
gsell committed
710

711
    } else if(inEntranceRegion && !inExitRegion) {
gsell's avatar
gsell committed
712

713 714 715 716 717
        Vector_t REntrance(0.0, 0.0, 0.0);
        if(IsPositionInEntranceField(R, REntrance)) {
            CalcEntranceFringeField(REntrance, deltaXEntrance, B);
        } else if(REntrance(2) > 0.0)
            CalcCentralField(R, deltaXEntrance, 0.0, B);
gsell's avatar
gsell committed
718

719
    } else if(!inEntranceRegion && inExitRegion) {
gsell's avatar
gsell committed
720

721 722 723 724 725
        Vector_t RExit(0.0, 0.0, 0.0);
        if(IsPositionInExitField(R, RExit)) {
            CalcExitFringeField(RExit, deltaXExit, B);
        } else if(RExit(2) < 0.0)
            CalcCentralField(R, deltaXExit, angle_m, B);
gsell's avatar
gsell committed
726

727
    } else if(inEntranceRegion && inExitRegion) {
gsell's avatar
gsell committed
728

729 730 731 732 733 734 735 736 737 738
        /*
         * This is an unusual condition and should only happen with
         * a sector magnet that bends more than 180 degrees. Here, we
         * have the possibility that the particle sees both the
         * entrance and exit fringe fields.
         */
        Vector_t BEntrance(0.0, 0.0, 0.0);
        Vector_t REntrance(0.0, 0.0, 0.0);
        if(IsPositionInEntranceField(R, REntrance))
            CalcEntranceFringeField(REntrance, deltaXEntrance, BEntrance);
739

740 741 742 743
        Vector_t BExit(0.0, 0.0, 0.0);
        Vector_t RExit(0.0, 0.0, 0.0);
        if(IsPositionInExitField(R, RExit))
            CalcExitFringeField(RExit, deltaXExit, BExit);
744

745 746 747
        B(0) = BEntrance(0) + BExit(0);
        B(1) = BEntrance(1) + BExit(1);
        B(2) = BEntrance(2) + BExit(2);
748 749 750 751

    }
}

752
void SBend::CalculateRefTrajectory(double &angleX, double &angleY) {
753 754

    const double mass = RefPartBunch_m->getM();
755
    const double gamma = designEnergy_m / mass + 1.;
756 757 758
    const double betaGamma = sqrt(gamma * gamma - 1.);
    const double dt = RefPartBunch_m->getdT();

759 760
    Vector_t X(0.0, 0.0, startField_m - elementEdge_m);
    Vector_t P(0.0, 0.0, betaGamma);
761

762 763 764 765 766 767
    if(!refTrajMapX_m.empty())
        refTrajMapX_m.clear();
    if(!refTrajMapY_m.empty())
        refTrajMapY_m.clear();
    if(!refTrajMapZ_m.empty())
        refTrajMapZ_m.clear();
768

769 770 771
    refTrajMapX_m.push_back(X(0));
    refTrajMapY_m.push_back(X(1));
    refTrajMapZ_m.push_back(X(2));
772

773 774 775
    refTrajMapStepSize_m = betaGamma / gamma * Physics::c * dt;
    double deltaS = 0.0;
    double bendLength = endField_m - startField_m;
776

777
    while(deltaS < bendLength) {
778 779 780 781 782

        X /= Vector_t(Physics::c * dt);
        pusher_m.push(X, P, dt);
        X *= Vector_t(Physics::c * dt);

783 784 785 786 787 788 789 790 791 792 793 794 795 796
        Vector_t eField(0.0, 0.0, 0.0);
        Vector_t bField(0.0, 0.0, 0.0);
        Vector_t XInBendFrame = RotateToBendFrame(X);

        /*
         * Add in transverse bend displacements. (ds is already
         * accounted for.)
         */
        XInBendFrame(0) -= dx_m;
        XInBendFrame(1) -= dy_m;

        CalculateMapField(XInBendFrame, eField, bField);
        bField = fieldAmplitude_m * RotateOutOfBendFrame(bField);

797
        X /= Vector_t(Physics::c * dt);
798 799
        pusher_m.kick(X, P, eField, bField, dt);

800 801 802
        pusher_m.push(X, P, dt);
        X *= Vector_t(Physics::c * dt);

803 804 805 806 807 808
        refTrajMapX_m.push_back(X(0));
        refTrajMapY_m.push_back(X(1));
        refTrajMapZ_m.push_back(X(2));

        deltaS += refTrajMapStepSize_m;

809 810
    }

811
    refTrajMapSize_m = refTrajMapX_m.size();
812

813
    if(std::abs(Orientation_m(2)) == Physics::pi / 2.0
814 815 816 817
       || Orientation_m(2) == 3.0 * Physics::pi / 2.0)
        angleX = 0.0;
    else
        angleX = -atan2(P(0), P(2));
818

819 820 821 822 823
    if(Orientation_m(2) == 0.0
       || Orientation_m(2) == Physics::pi)
        angleY = 0.0;
    else
        angleY = atan2(P(1), P(2));
824 825 826

}

827 828 829
double SBend::EstimateFieldAdjustmentStep(double actualBendAngle,
        double mass,
        double betaGamma) {
830

831 832
    double amplitude1 = fieldAmplitude_m;
    double bendAngle1 = actualBendAngle;
833

834 835 836 837 838 839
    // Estimate field adjustment step.
    double effectiveLength = angle_m * designRadius_m;
    double fieldStep = (angle_m - bendAngle1) * betaGamma * mass / (2.0 * effectiveLength * Physics::c);
    if(pow(fieldAmplitude_m * effectiveLength * Physics::c / (betaGamma * mass), 2.0) < 1.0)
        fieldStep = (angle_m - bendAngle1) * betaGamma * mass / (2.0 * effectiveLength * Physics::c)
                    * std::sqrt(1.0 - pow(fieldAmplitude_m * effectiveLength * Physics::c / (betaGamma * mass), 2.0));
840

841
    fieldStep *= amplitude1 / std::abs(amplitude1);
842

843
    return fieldStep;
844 845 846

}

847
void SBend::FindBendEffectiveLength(double startField, double endField) {
848

849 850 851 852 853 854 855 856
    /*
     * Use an iterative procedure to set the width of the
     * default field map for the defined field amplitude
     * and bend angle.
     */
    SetEngeOriginDelta(0.0);
    SetFieldCalcParam(false);
    SetFieldBoundaries(startField, endField);
857

858 859
    double actualBendAngle = CalculateBendAngle();
    double error = std::abs(actualBendAngle - angle_m);
860 861

    if(error > 1.0e-6) {
gsell's avatar
gsell committed
862

863 864 865 866 867 868 869
        double deltaStep = 0.0;
        if(std::abs(actualBendAngle) < std::abs(angle_m))
            deltaStep = -gap_m / 2.0;
        else
            deltaStep = gap_m / 2.0;

        double delta1 = 0.0;
870 871
        double bendAngle1 = actualBendAngle;

872 873 874 875 876 877 878 879 880 881 882 883 884
        double delta2 = deltaStep;
        SetEngeOriginDelta(delta2);
        SetFieldCalcParam(false);
        SetFieldBoundaries(startField, endField);
        double bendAngle2 = CalculateBendAngle();

        if(std::abs(bendAngle1) > std::abs(angle_m)) {
            while(std::abs(bendAngle2) > std::abs(angle_m)) {
                delta2 += deltaStep;
                SetEngeOriginDelta(delta2);
                SetFieldCalcParam(false);
                SetFieldBoundaries(startField, endField);
                bendAngle2 = CalculateBendAngle();
885 886
            }
        } else {
887 888 889 890 891 892
            while(std::abs(bendAngle2) < std::abs(angle_m)) {
                delta2 += deltaStep;
                SetEngeOriginDelta(delta2);
                SetFieldCalcParam(false);
                SetFieldBoundaries(startField, endField);
                bendAngle2 = CalculateBendAngle();
893 894 895
            }
        }

896
        // Now we should have the proper field map width bracketed.
897
        unsigned int iterations = 1;
898 899
        double delta = 0.0;
        error = std::abs(actualBendAngle - angle_m);
900 901
        while(error > 1.0e-6 && iterations < 100) {

902 903 904 905 906
            delta = (delta1 + delta2) / 2.0;
            SetEngeOriginDelta(delta);
            SetFieldCalcParam(false);
            SetFieldBoundaries(startField, endField);
            double newBendAngle = CalculateBendAngle();
907

908
            error = std::abs(newBendAngle - angle_m);
909 910 911 912 913 914 915

            if(error > 1.0e-6) {

                if(bendAngle1 - angle_m < 0.0) {

                    if(newBendAngle - angle_m < 0.0) {
                        bendAngle1 = newBendAngle;
916
                        delta1 = delta;
917 918
                    } else {
                        bendAngle2 = newBendAngle;
919
                        delta2 = delta;
920 921 922 923 924 925
                    }

                } else {

                    if(newBendAngle - angle_m < 0.0) {
                        bendAngle2 = newBendAngle;
926
                        delta2 = delta;
927 928
                    } else {
                        bendAngle1 = newBendAngle;
929
                        delta1 = delta;
930 931 932 933 934 935 936
                    }
                }
            }
            iterations++;
        }
    }
}
937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 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 1044 1045 1046 1047 1048 1049 1050

void SBend::FindBendStrength(double mass,
                             double gamma,
                             double betaGamma,
                             double charge) {

    /*
     * Use an iterative procedure to set the magnet field amplitude
     * for the defined bend angle.
     */
    double actualBendAngle = CalculateBendAngle();
    double fieldStep = EstimateFieldAdjustmentStep(actualBendAngle,
                       mass,
                       betaGamma);
    double amplitude1 = fieldAmplitude_m;
    double bendAngle1 = actualBendAngle;

    double amplitude2 = fieldAmplitude_m + fieldStep;
    fieldAmplitude_m = amplitude2;
    double bendAngle2 = CalculateBendAngle();

    if(std::abs(bendAngle1) > std::abs(angle_m)) {
        while(std::abs(bendAngle2) > std::abs(angle_m)) {
            amplitude2 += fieldStep;
            fieldAmplitude_m = amplitude2;
            bendAngle2 = CalculateBendAngle();
        }
    } else {
        while(std::abs(bendAngle2) < std::abs(angle_m)) {
            amplitude2 += fieldStep;
            fieldAmplitude_m = amplitude2;
            bendAngle2 = CalculateBendAngle();
        }
    }

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

        fieldAmplitude_m = (amplitude1 + amplitude2) / 2.0;
        double newBendAngle = CalculateBendAngle();

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

        if(error > 1.0e-6) {

            if(bendAngle1 - angle_m < 0.0) {

                if(newBendAngle - angle_m < 0.0) {
                    bendAngle1 = newBendAngle;
                    amplitude1 = fieldAmplitude_m;
                } else {
                    bendAngle2 = newBendAngle;
                    amplitude2 = fieldAmplitude_m;
                }

            } else {

                if(newBendAngle - angle_m < 0.0) {
                    bendAngle2 = newBendAngle;
                    amplitude2 = fieldAmplitude_m;
                } else {
                    bendAngle1 = newBendAngle;
                    amplitude1 = fieldAmplitude_m;
                }
            }
        }
        iterations++;
    }
}

bool SBend::FindChordLength(Inform &msg,
                            double &chordLength,
                            bool &chordLengthFromMap) {

    /*
     * Find bend chord length. If this was not set by the user using the
     * L (length) attribute, infer it from the field map.
     */
    chordLength = length_m;
    if(chordLength > 0.0) {
        chordLengthFromMap = false;
        return true;
    } else {

        if(chordLength == 0.0)
            chordLength = exitParameter2_m - entranceParameter2_m;

        chordLengthFromMap = true;

        if(chordLength <= 0.0) {
            msg << "Magnet length inferred from field map is less than or equal"
                " to zero. Check your bend magnet input."
                << endl;
            return false;
        } else
            return true;

    }
}

bool SBend::FindIdealBendParameters(double chordLength) {

    double refMass = RefPartBunch_m->getM();
    double refGamma = designEnergy_m / refMass + 1.0;
    double refBetaGamma = sqrt(pow(refGamma, 2.0) - 1.0);
    double refCharge = RefPartBunch_m->getQ();

    if(angle_m != 0.0) {

        if(angle_m < 0.0) {
            // Negative angle is a positive bend rotated 180 degrees.
            angle_m = std::abs(angle_m);