RBend.cpp 21.6 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 22 23 24 25 26 27 28 29 30 31 32 33 34 35
// ------------------------------------------------------------------------
// $RCSfile: RBend.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: RBend
//   Defines the abstract interface for a rectangular bend magnet.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------

#include "AbsBeamline/RBend.h"
#include "AbsBeamline/BeamlineVisitor.h"
#include "Fields/Fieldmap.hh"
#include <iostream>
#include <fstream>

extern Inform *gmsg;

// Class RBend
// ------------------------------------------------------------------------

int RBend::RBend_counter_m = 0;

RBend::RBend():
    Component(),
36 37 38
    length_m(0.0),
    gap_m(0.0),
    reinitialize_m(false),
gsell's avatar
gsell committed
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
    filename_m(""),
    fieldmap_m(NULL),
    amplitude_m(0.0),
    field_orientation_m(1.0, 0.0),
    startField_m(0.0),
    endField_m(0.0),
    fast_m(false),
    sin_face_alpha_m(0.0),
    cos_face_alpha_m(1.0),
    tan_face_alpha_m(0.0),
    sin_face_beta_m(0.0),
    cos_face_beta_m(1.0),
    tan_face_beta_m(0.0),
    design_energy_m(0.0),
    angle_m(0.0),
    map_m(NULL),
    map_size_m(0),
    map_step_size_m(0.0),
    pusher_m(),
    startElement_m(0.0),
    R_m(0.0),
    effectiveLength_m(0.0),
    effectiveCenter_m(0.0) {
    setElType(isDipole);
}


RBend::RBend(const RBend &right):
    Component(right),
68 69 70
    length_m(right.length_m),
    gap_m(right.gap_m),
    reinitialize_m(right.reinitialize_m),
gsell's avatar
gsell committed
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
    filename_m(right.filename_m),
    fieldmap_m(right.fieldmap_m),
    amplitude_m(right.amplitude_m),
    field_orientation_m(right.field_orientation_m),
    startField_m(right.startField_m),
    endField_m(right.endField_m),
    fast_m(right.fast_m),
    sin_face_alpha_m(right.sin_face_alpha_m),
    cos_face_alpha_m(right.cos_face_alpha_m),
    tan_face_alpha_m(right.tan_face_alpha_m),
    sin_face_beta_m(right.sin_face_beta_m),
    cos_face_beta_m(right.cos_face_beta_m),
    tan_face_beta_m(right.tan_face_beta_m),
    design_energy_m(right.design_energy_m),
    angle_m(right.angle_m),
    map_size_m(right.map_size_m),
    map_step_size_m(right.map_step_size_m),
    pusher_m(right.pusher_m),
    startElement_m(right.startElement_m),
    R_m(right.R_m),
    effectiveLength_m(right.effectiveLength_m),
    effectiveCenter_m(right.effectiveCenter_m) {
    setElType(isDipole);
    if(map_size_m > 0) {
        map_m = new double[map_size_m + 1];
        for(int i = 0; i < map_size_m + 1; ++i)
            map_m[i] = right.map_m[i];
    } else {
        map_m = NULL;
    }
}


Steve Russell's avatar
Steve Russell committed
104
RBend::RBend(const std::string &name):
gsell's avatar
gsell committed
105
    Component(name),
106 107 108
    length_m(0.0),
    gap_m(0.0),
    reinitialize_m(false),
gsell's avatar
gsell committed
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 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 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
    filename_m(""),
    fieldmap_m(NULL),
    amplitude_m(0.0),
    field_orientation_m(1.0, 0.0),
    startField_m(0.0),
    endField_m(0.0),
    fast_m(false),
    sin_face_alpha_m(0.0),
    cos_face_alpha_m(1.0),
    tan_face_alpha_m(0.0),
    sin_face_beta_m(0.0),
    cos_face_beta_m(1.0),
    tan_face_beta_m(0.0),
    design_energy_m(0.0),
    angle_m(0.0),
    map_m(NULL),
    map_size_m(0),
    map_step_size_m(0.0),
    pusher_m(),
    startElement_m(0.0),
    R_m(0.0),
    effectiveLength_m(0.0),
    effectiveCenter_m(0.0) {
    setElType(isDipole);
}


RBend::~RBend() {
    if(map_m)
        delete[] map_m;
}


void RBend::accept(BeamlineVisitor &visitor) const {
    visitor.visitRBend(*this);
}


double RBend::getNormalComponent(int n) const {
    return getField().getNormalComponent(n);
}


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


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


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

//Implemented BET functions for dipole
//ff
//transverse kick
void RBend::addKR(int i, double t, Vector_t &K) {
    Inform msg("RBend::addK()");

    Vector_t tmpE(0.0, 0.0, 0.0);
    Vector_t tmpB(0.0, 0.0, 0.0);
    Vector_t tmpE_diff(0.0, 0.0, 0.0);
    Vector_t tmpB_diff(0.0, 0.0, 0.0);
    double pz = RefPartBunch_m->getZ(i) - startField_m - ds_m;
    const Vector_t tmpA(RefPartBunch_m->getX(i) - dx_m, RefPartBunch_m->getY(i) - dy_m, pz);

    DiffDirection zdir(DZ);
    myFieldmap->getFieldstrength(tmpA, tmpE, tmpB);
    myFieldmap->getFieldstrength_fdiff(tmpA, tmpE_diff, tmpB_diff, zdir);

    double g = RefPartBunch_m->getGamma(i);

    if(fabs(scale_m * tmpB_diff(2)) > 0.1) {
        double cf = Physics::q_e * tmpB(2) / (g * Physics::EMASS);
        K += Vector_t(-pow(cf * scale_m * tmpB(0), 2) / 3.0, -pow(cf * scale_m * tmpB(1), 2) / 3.0, 0.0);
    }
}

void RBend::addKT(int i, double t, Vector_t &K) {
    Inform msg("RBend::addK()");

    Vector_t tmpE(0.0, 0.0, 0.0);
    Vector_t tmpB(0.0, 0.0, 0.0);
    double pz = RefPartBunch_m->getZ(i) - startField_m - ds_m;
    const Vector_t tmpA(RefPartBunch_m->getX(i) - dx_m, RefPartBunch_m->getY(i) - dy_m, pz);
    myFieldmap->getFieldstrength(tmpA, tmpE, tmpB);

    double b = RefPartBunch_m->getBeta(i);
    double g = 1 / sqrt(1 - b * b);

    double cf = -Physics::q_e * Physics::c * b * tmpB(2) * scale_m / (g * Physics::EMASS);
    Vector_t temp(cf * tmpB(1), cf * tmpB(0), 0.0);

    //FIXME: K += ??
}

bool RBend::apply(const int &i, const double &t, double E[], double B[]) {
    Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
    if(apply(i, t, Ev, Bv)) return true;

    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;
}

bool RBend::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {

Steve Russell's avatar
Steve Russell committed
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
    // If this is the first call, the bend angle is specified in the input
    // file and the design energy of the bend is different from the average
    // energy of the beam, we reinitialize the bend.
    if(reinitialize_m) {
        if(design_energy_m != RefPartBunch_m->get_meanEnergy() * 1.0e6) {
            design_energy_m = RefPartBunch_m->get_meanEnergy() * 1.0e6;

            setBendStrength();

            double zBegin = 0.0;
            double zEnd = 0.0;
            double rBegin = 0.0;
            double rEnd = 0.0;
            fieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);
            calculateRefTrajectory(zBegin);

            Inform msg("RBend ");
            msg << "Bend design energy changed to: " << design_energy_m * 1.0e-6 << " MeV" << endl;
            msg << "Field amplitude:               " << amplitude_m << " T" << endl;
        }

        reinitialize_m = false;
    }

gsell's avatar
gsell committed
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 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
    const Vector_t &X = RefPartBunch_m->X[i];
    Vector_t strength(0.0), info(0.0);

    fieldmap_m->getFieldstrength(X, strength, info);

    if(info(0) > 0.99) {

        B(1) += amplitude_m * (strength(0) - strength(2) / 2. * X(1) * X(1));
        if(info(1) > 0.99) {
            B(2) += amplitude_m * strength(1) * X(1);
        } else {
            B(2) -= amplitude_m * strength(1) * X(1);
        }
    } else if(fabs(info(0)) < 0.01)
        B(1) += amplitude_m;

    return false;
}

bool RBend::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {

    int index = (int)floor((R(2) - startField_m) / map_step_size_m);
    if(index > 0 && index + 1 < map_size_m) {

        // Find indices for z position in pre-computed central trajectory map.
        double lever = (R(2) - startField_m) / map_step_size_m - index;
        double z = (1. - lever) * map_m[index] + lever * map_m[index + 1];

        // Rotate x and y to the the bend's local coordinate system.
        //
        // 1) Rotate about the z axis by angle negative Orientation_m(2).
        // 2) Rotate about the y axis by angle negative Orientation_m(0).
        // 3) Rotate about the x axis by angle Orientation_m(1).

        const double sina = sin(Orientation_m(0));
        const double cosa = cos(Orientation_m(0));
        const double sinb = sin(Orientation_m(1));
        const double cosb = cos(Orientation_m(1));
        const double sinc = sin(Orientation_m(2));
        const double cosc = cos(Orientation_m(2));

        Vector_t X(0.0);
        X(0) = (cosa * cosc) * R(0) + (cosa * sinc) * R(1) - sina *        R(2);
        X(1) = (-cosb * sinc - sina * sinb * cosc) * R(0) + (cosb * cosc - sina * sinb * sinc) * R(1) - cosa * sinb * R(2);
        X(2) = z;

        Vector_t strength(0.0), info(0.0);

        fieldmap_m->getFieldstrength(X, strength, info);
        Vector_t tempB(0.0);

        if(info(0) > 0.99) {

            tempB(1) = amplitude_m * (strength(0) - strength(2) / 2. * X(1) * X(1));
            if(info(1) > 0.99) {
                tempB(2) = amplitude_m * strength(1) * X(1);
            } else {
                tempB(2) = -amplitude_m * strength(1) * X(1);
            }
        } else if(fabs(info(0)) < 0.01)
            tempB(1) = amplitude_m;

        // Rotate field out of the bend's local coordinate system and back to lab frame.
        //
        // 1) Rotate about the x axis by angle Orientation_m(1).
        // 2) Rotate about the y axis by angle Orientation_m(0).
        // 3) Rotate about the z axis by angle negative Orientation_(2).

        B(0) +=  cosa * cosc * tempB(0) + (-sina * sinb * cosc - cosb * sinc) * tempB(1) + (sina * cosb * cosc - sinb * sinc) * tempB(2);
        B(1) +=  cosa * sinc * tempB(0) + (-sina * sinb * sinc + cosb * cosc) * tempB(1) + (sina * cosb * sinc + sinb * cosc) * tempB(2);
        B(2) += -sina *        tempB(0) + (-cosa * sinb) * tempB(1) + (cosa * cosb) * tempB(2);
    }
    return false;
}

void RBend::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {

    Inform msg("RBend ");

    double zBegin = 0.0;
    double zEnd = 0.0;
    double rBegin = 0.0;
    double rEnd = 0.0;
Steve Russell's avatar
Steve Russell committed
332

gsell's avatar
gsell committed
333 334 335 336 337 338 339 340 341 342
    startElement_m = startField;

    RefPartBunch_m = bunch;
    pusher_m.initialise(bunch->getReference());

    fieldmap_m = Fieldmap::getFieldmap(filename_m, fast_m);
    fieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);

    if((fieldmap_m != NULL) && (zEnd > zBegin)) {

Steve Russell's avatar
Steve Russell committed
343
        // Read in field map.
gsell's avatar
gsell committed
344 345 346 347
        msg << getName() << " using file ";
        fieldmap_m->getInfo(&msg);
        Fieldmap::readMap(filename_m);

Steve Russell's avatar
Steve Russell committed
348 349 350 351 352
        // Check that the design energy is greater than zero.
        if(design_energy_m <= 0.0) {
            msg << "The bend must have a design enregy greater than zero set in the input file." << endl;
            return;
        }
gsell's avatar
gsell committed
353

Steve Russell's avatar
Steve Russell committed
354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369
        // If using default field map, set length and gap.
        if(filename_m == "1DPROFILE1-DEFAULT") {
            if(gap_m <= 0.0 || length_m <= 0.0) {
                msg << "If using \"1DPROFILE1-DEFAULT\" field map you must set GAP (full magnet gap) and L (length) in the OPAL input file." << endl;
                return;
            } else {
                fieldmap_m->setFieldGap(gap_m);
                fieldmap_m->setFieldLength(length_m);
                fieldmap_m->adjustFringeFields();
                msg << "Adjusted fringe field parameters." << endl;
                fieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);
                fieldmap_m->getInfo(&msg);
                startElement_m = startField;
            }
        }
        fieldmap_m->setEdgeConstants(0.0, 0.0, 0.0);
gsell's avatar
gsell committed
370 371

        length_m = zEnd - zBegin;
Steve Russell's avatar
Steve Russell committed
372 373 374 375
        if(length_m < 0.0) {
            // There is probably something wrong with the fieldmap.
            return;
        }
gsell's avatar
gsell committed
376

Steve Russell's avatar
Steve Russell committed
377 378 379 380 381 382 383
        // If the bend angle is specified, find proper field strength. If only
        // the field strength is given, just calculate the bend angle. This also
        // sets the bend exit angle appropriately.
        if(angle_m != 0.0) {
            if(angle_m < 0.0) {
                angle_m *= -1.0;
                Orientation_m(2) += Physics::pi;
gsell's avatar
gsell committed
384
            }
Steve Russell's avatar
Steve Russell committed
385 386 387 388 389
            setBendStrength();
            reinitialize_m = true;
        } else {
            angle_m = calculateBendAngle(length_m);
            reinitialize_m = false;
gsell's avatar
gsell committed
390 391
        }

Steve Russell's avatar
Steve Russell committed
392 393
        // Calculate the reference particle trajectory map.
        double bendAngle = calculateRefTrajectory(zBegin);
gsell's avatar
gsell committed
394

Steve Russell's avatar
Steve Russell committed
395 396
        startField = startField_m;
        endField = endField_m;
gsell's avatar
gsell committed
397

Steve Russell's avatar
Steve Russell committed
398 399 400 401 402 403 404
        msg << "Start:            " << startField_m << " m (in floor coordinates)" << endl;
        msg << "End:              " << endField_m << " m (in floor coordinates)" << endl;
        msg << "Bend angle:       " << bendAngle * 180.0 / Physics::pi << " degrees" << endl;
        msg << "Field amplitude:  " << amplitude_m << " T" << endl;
        msg << "Bend radius:      " << R_m << " m" << endl;
        msg << "Effective length: " << effectiveLength_m << " m (in s coordinates)" << endl;
        msg << "Effective center: " << effectiveCenter_m << " m (in s coordinates with respect to bend field map start position)" << endl;
gsell's avatar
gsell committed
405 406 407 408 409 410 411 412 413 414 415 416 417 418

    } else {
        endField = startField - 1e-3;
    }
}

void RBend::finalise() {
    online_m = false;
}

bool RBend::bends() const {
    return true;
}

Steve Russell's avatar
Steve Russell committed
419 420 421 422
void RBend::setBendAngle(const double &angle) {
    angle_m = angle * Physics::pi / 180.0;
}

gsell's avatar
gsell committed
423 424 425 426
void RBend::setAmplitudem(double vPeak) {
    amplitude_m = vPeak;
}

Steve Russell's avatar
Steve Russell committed
427 428 429 430 431 432 433 434 435
void RBend::setFullGap(const double &gap) {
    gap_m = gap;
}

void RBend::setLength(const double &length) {
    length_m = length;
}

void RBend::setFieldMapFN(std::string fmapfn) {
gsell's avatar
gsell committed
436 437 438
    filename_m = fmapfn;
}

Steve Russell's avatar
Steve Russell committed
439
std::string RBend::getFieldMapFN() const {
gsell's avatar
gsell committed
440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
    return filename_m;
}

void RBend::getDimensions(double &zBegin, double &zEnd) const {
    zBegin = startField_m;
    zEnd = endField_m;
}

double RBend::getEffectiveLength() const {
    return effectiveLength_m;
}

double RBend::getEffectiveCenter() const {
    return effectiveCenter_m;
}

double RBend::getBendAngle() const {
    return angle_m;
}

double RBend::getStartElement() const {
    return startElement_m;
}

double RBend::getR() const {
    return R_m;
}


Steve Russell's avatar
Steve Russell committed
469 470
const std::string &RBend::getType() const {
    static const std::string type("RBend");
gsell's avatar
gsell committed
471 472 473
    return type;
}

Steve Russell's avatar
Steve Russell committed
474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623
void RBend::setBendStrength() {
    // This routine uses an iterative procedure to set the bend strength
    // so that the bend angle is the one we want.
    //
    // This is a primitive approach in that it is not very efficient. But it
    // is stable and is only called a few times during a simulation.

    // Estimate bend field magnitude.
    const double mass = RefPartBunch_m->getM();
    const double gamma = design_energy_m / mass + 1.0;
    const double betaGamma = sqrt(pow(gamma, 2.0) - 1.0);
    const double charge = RefPartBunch_m->getQ();

    calculateEffectiveLength();
    amplitude_m = (charge / fabs(charge)) * betaGamma * mass / (Physics::c * (effectiveLength_m * cos_face_alpha_m) / sin(angle_m));

    // Find initial angle.
    double zBegin = 0.0;
    double zEnd = 0.0;
    double rBegin = 0.0;
    double rEnd = 0.0;
    fieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);
    double actualBendAngle = calculateBendAngle(zEnd - zBegin);

    // Adjust field amplitude to get desired bend angle.
    int iterations = 1;
    double fieldAdjustment = amplitude_m / 10.0;

    if(fabs(actualBendAngle) > fabs(angle_m))
        fieldAdjustment *= -1.0;

    bool lastGreater = true;
    if(fabs(actualBendAngle) < fabs(angle_m))
        lastGreater = false;

    while(fabs(actualBendAngle - angle_m) > 1.0e-8 && iterations <= 1000) {

        actualBendAngle = calculateBendAngle(zEnd - zBegin);
        iterations++;

        if((!lastGreater && fabs(actualBendAngle) > fabs(angle_m)) || (lastGreater && fabs(actualBendAngle) < fabs(angle_m)))
            fieldAdjustment /= -10.0;

        if(fabs(actualBendAngle) > fabs(angle_m)) lastGreater = true;
        else lastGreater = false;

        amplitude_m += fieldAdjustment;
    }
}

double RBend::calculateBendAngle(double bendLength) {
    // This routine calculates the bend angle using an iterative process.

    // Calculate angle.
    const double mass = RefPartBunch_m->getM();
    const double gamma = design_energy_m / mass + 1.0;
    const double betaGamma = sqrt(pow(gamma, 2.0) - 1.0);
    const double deltaT = RefPartBunch_m->getdT();

    // Integrate through field for initial angle.
    Vector_t X(0.0, 0.0, 0.0);
    Vector_t P(-betaGamma * sin_face_alpha_m, 0.0, betaGamma * cos_face_alpha_m);
    Vector_t strength(0.0, 0.0, 0.0);
    Vector_t bField(0.0, 0.0, 0.0);
    Vector_t temp(0.0, 0.0, 0.0);

    while(P(2) > 0.0 && X(2) < bendLength) {

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

        fieldmap_m->getFieldstrength(X, strength, temp);
        bField(1) = amplitude_m * strength(0);
        temp = Vector_t(0.0);

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

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

    }

    double angle =  -atan2(P(0), P(2)) - Orientation_m(0);

    return angle;
}

double RBend::calculateRefTrajectory(const double zBegin) {

    // Calculate the reference trajectory map.
    const double mass = RefPartBunch_m->getM();
    const double gamma = design_energy_m / mass + 1.;
    const double betagamma = sqrt(gamma * gamma - 1.);
    const double dt = RefPartBunch_m->getdT();
    int j = 0;

    Vector_t tmp(0.0);
    Vector_t Bfield(0.0);
    Vector_t strength(0.0);
    Vector_t X(0.0);
    Vector_t P(-betagamma * sin_face_alpha_m, 0.0, betagamma * cos_face_alpha_m); // TODO: make it 3D

    bool EntryFringe_passed = false;
    double PathLengthEntryFringe = 0.0;  // in S coordinates. This value is different from zBegin due to the curvature!

    if(map_m != NULL) delete map_m;

    map_step_size_m = betagamma / gamma * Physics::c * dt;
    map_size_m = static_cast<int>(floor(length_m / 2. * Physics::pi / map_step_size_m));
    map_m = new double[map_size_m + 1];
    map_m[0] = 0.0;

    while(map_m[j] < length_m && j < map_size_m) {
        strength = Vector_t(0.0);
        X /= Vector_t(Physics::c * dt);
        pusher_m.push(X, P, dt);
        X *= Vector_t(Physics::c * dt);

        fieldmap_m->getFieldstrength(X, strength, tmp);
        if(X(2) >= fabs(zBegin) && !EntryFringe_passed) {

            // This is the point where we pass ELEMEDGE
            // not the end of the entry fringe field as the
            // name suggests.

            EntryFringe_passed = true;
            PathLengthEntryFringe = j * map_step_size_m;
        }
        Bfield(1) = amplitude_m * strength(0);
        tmp = Vector_t(0.0);
        X /= Vector_t(Physics::c * dt);
        pusher_m.kick(X, P, tmp, Bfield, dt);
        pusher_m.push(X, P, dt);
        X *= Vector_t(Physics::c * dt);

        map_m[++j] = X(2);

    }

    map_size_m = j;
    double angle = -atan2(P(0), P(2)) - Orientation_m(0);

    startField_m = startElement_m - PathLengthEntryFringe;
    endField_m = startField_m + map_step_size_m * j;

    // Set "ideal" bend radius and effective length.
    R_m = fabs(betagamma * mass / (Physics::c * amplitude_m));
624
    effectiveLength_m = R_m * std::abs(angle);
Steve Russell's avatar
Steve Russell committed
625 626 627 628
    calculateEffectiveCenter();

    return angle;
}
gsell's avatar
gsell committed
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674

void RBend::calculateEffectiveLength() {

    double zBegin = 0.0;
    double zEnd = 0.0;
    double rBegin = 0.0;
    double rEnd = 0.0;
    fieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);

    // Uses Simpson's rule to integrate field. Make step size about 1 mm.

    // This must be odd.
    unsigned int numberOfIntSteps = 2 * static_cast<unsigned int>(floor((zEnd - zBegin) * 1000.0 / 2.0)) + 1;

    double deltaZ = (zEnd - zBegin) / numberOfIntSteps;
    effectiveLength_m = 0.0;

    for(unsigned int integralIndex = 1; integralIndex <= (numberOfIntSteps - 1) / 2; integralIndex++) {

        Vector_t strength(0.0);
        Vector_t info(0.0);
        Vector_t X(0.0);
        X(2) = (2 * integralIndex - 1) * deltaZ;
        fieldmap_m->getFieldstrength(X, strength, info);
        double field1 = strength(0);

        X(2) = 2 * integralIndex * deltaZ;
        fieldmap_m->getFieldstrength(X, strength, info);
        double field2 = strength(0);

        X(2) = (2 * integralIndex + 1) * deltaZ;
        fieldmap_m->getFieldstrength(X, strength, info);
        double field3 = strength(0);

        effectiveLength_m += deltaZ * (field1 + 4.0 * field2 + field3) / 3.0;
    }
}

void RBend::calculateEffectiveCenter() {

    double zBegin = 0.0;
    double zEnd = 0.0;
    double rBegin = 0.0;
    double rEnd = 0.0;
    fieldmap_m->getFieldDimensions(zBegin, zEnd, rBegin, rEnd);

Steve Russell's avatar
Steve Russell committed
675 676
    // Initial guess for effective center.
    double effectiveCenter = fabs(R_m * angle_m / 2.0) - zBegin;
gsell's avatar
gsell committed
677

Steve Russell's avatar
Steve Russell committed
678 679
    // Find initial angle.
    double actualBendAngle = calculateBendAngle(effectiveCenter);
gsell's avatar
gsell committed
680

Steve Russell's avatar
Steve Russell committed
681 682 683
    // Adjust effective center to get a bend angle 0.5 times the full bend angle.
    int iterations = 1;
    double lengthAdjustment = effectiveCenter / 10.0;
gsell's avatar
gsell committed
684

Steve Russell's avatar
Steve Russell committed
685 686
    if(fabs(actualBendAngle) > fabs(angle_m / 2.0))
        lengthAdjustment *= -1.0;
gsell's avatar
gsell committed
687

Steve Russell's avatar
Steve Russell committed
688 689 690
    bool lastGreater = true;
    if(fabs(actualBendAngle) < fabs(angle_m / 2.0))
        lastGreater = false;
gsell's avatar
gsell committed
691

Steve Russell's avatar
Steve Russell committed
692
    while(fabs(actualBendAngle - angle_m / 2.0) > 1.0e-8 && iterations <= 100) {
gsell's avatar
gsell committed
693

Steve Russell's avatar
Steve Russell committed
694 695
        actualBendAngle = calculateBendAngle(effectiveCenter);
        iterations++;
gsell's avatar
gsell committed
696

Steve Russell's avatar
Steve Russell committed
697 698
        if((!lastGreater && fabs(actualBendAngle) > fabs(angle_m / 2.0)) || (lastGreater && fabs(actualBendAngle) < fabs(angle_m / 2.0)))
            lengthAdjustment /= -10.0;
gsell's avatar
gsell committed
699

Steve Russell's avatar
Steve Russell committed
700 701
        if(fabs(actualBendAngle) > fabs(angle_m / 2.0)) lastGreater = true;
        else lastGreater = false;
gsell's avatar
gsell committed
702

Steve Russell's avatar
Steve Russell committed
703
        effectiveCenter += lengthAdjustment;
gsell's avatar
gsell committed
704 705

    }
Steve Russell's avatar
Steve Russell committed
706
    effectiveCenter_m = effectiveCenter - R_m * sin(angle_m / 2.0) + R_m * angle_m / 2.0;
gsell's avatar
gsell committed
707
}