OpalBeamline.cpp 24.8 KB
Newer Older
gsell's avatar
gsell committed
1
#include "Elements/OpalBeamline.h"
2

3 4 5
#include "Utilities/Options.h"
#include "Utilities/Util.h"
#include "AbstractObjects/OpalData.h"
6
#include "AbsBeamline/Bend2D.h"
7 8 9 10 11
#include "Structure/MeshGenerator.h"

#include <boost/filesystem.hpp>
#include <boost/regex.hpp>
#include <fstream>
gsell's avatar
gsell committed
12 13 14

OpalBeamline::OpalBeamline():
    elements_m(),
15
    prepared_m(false),
16 17 18 19 20
    containsSource_m(false)
{
}

OpalBeamline::OpalBeamline(const Vector_t& origin,
21
                           const Quaternion& rotation):
22 23 24
    elements_m(),
    prepared_m(false),
    containsSource_m(false),
25
    coordTransformationTo_m(origin, rotation)
26
{
gsell's avatar
gsell committed
27 28 29 30 31 32
}

OpalBeamline::~OpalBeamline() {
    elements_m.clear();
}

33 34 35 36 37 38 39
std::set<std::shared_ptr<Component>> OpalBeamline::getElements(const Vector_t &x) {
    std::set<std::shared_ptr<Component> > elementSet;
    FieldList::iterator it = elements_m.begin();
    const FieldList::iterator end = elements_m.end();
    for (; it != end; ++ it) {
        std::shared_ptr<Component> element = (*it).getElement();
        Vector_t r = (*it).getCoordTransformationTo().transformTo(x);
40

41 42
        if (element->isInside(r)) {
            elementSet.insert(element);
gsell's avatar
gsell committed
43 44 45
        }
    }

46
    return elementSet;
gsell's avatar
gsell committed
47 48 49 50 51 52
}

unsigned long OpalBeamline::getFieldAt(const unsigned int &index, const Vector_t &pos, const long &sindex, const double &t, Vector_t &E, Vector_t &B) {

    unsigned long rtv = 0x00;

53 54
    return rtv;
}
gsell's avatar
gsell committed
55

56 57 58 59 60 61
unsigned long OpalBeamline::getFieldAt(const Vector_t &position,
                                       const Vector_t &momentum,
                                       const double &t,
                                       Vector_t &Ef,
                                       Vector_t &Bf) {
    unsigned long rtv = 0x00;
gsell's avatar
gsell committed
62

63
    std::set<std::shared_ptr<Component>> elements = getElements(position);
64

65 66
    std::set<std::shared_ptr<Component>>::const_iterator it = elements.begin();
    const std::set<std::shared_ptr<Component>>::const_iterator end = elements.end();
gsell's avatar
gsell committed
67

68 69 70 71
    for (; it != end; ++ it) {
        ElementBase::ElementType type = (*it)->getType();
        if (type == ElementBase::MONITOR ||
            type == ElementBase::MARKER ||
72
            type == ElementBase::CCOLLIMATOR ||
73
            type == ElementBase::DIAGNOSTIC) continue;
gsell's avatar
gsell committed
74

75 76 77
        Vector_t localR = transformToLocalCS(*it, position);
        Vector_t localP = rotateToLocalCS(*it, momentum);
        Vector_t localE(0.0), localB(0.0);
gsell's avatar
gsell committed
78

79
        (*it)->applyToReferenceParticle(localR, localP, t, localE, localB);
gsell's avatar
gsell committed
80

81 82
        Ef += rotateFromLocalCS(*it, localE);
        Bf += rotateFromLocalCS(*it, localB);
gsell's avatar
gsell committed
83 84
    }

85 86 87
    //         if(section.hasWake()) {
    //             rtv |= BEAMLINE_WAKE;
    //         }
88
    //         if(section.hasParticleMatterInteraction()) {
89
    //             rtv |= BEAMLINE_PARTICLEMATTERINTERACTION;
90
    //         }
gsell's avatar
gsell committed
91 92 93 94

    return rtv;
}

kraus's avatar
kraus committed
95
void OpalBeamline::switchElements(const double &min, const double &max, const double &kineticEnergy, const bool &nomonitors) {
96 97

    FieldList::iterator fprev;
gsell's avatar
gsell committed
98 99 100
    for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) {
        // don't set online monitors if the centroid of the bunch is allready inside monitor
        // or if explicitly not desired (eg during auto phasing)
101
        if(flit->getElement()->getType() == ElementBase::MONITOR) {
gsell's avatar
gsell committed
102 103 104
            double spos = (max + min) / 2.;
            if(!nomonitors && spos < (*flit).getStart()) {
                if(!(*flit).isOn() && max > (*flit).getStart()) {
kraus's avatar
kraus committed
105
                    (*flit).setOn(kineticEnergy);
gsell's avatar
gsell committed
106 107 108 109
                }
            }

        } else {
110
            if(!(*flit).isOn() && max > (*flit).getStart() && min < (*flit).getEnd()) {
kraus's avatar
kraus committed
111
                (*flit).setOn(kineticEnergy);
gsell's avatar
gsell committed
112
            }
113

114 115 116 117 118 119 120 121 122
            //check if multiple degraders follow one another with no other elements in between
            //if element is off and it is a degrader
            if (!(*flit).isOn() && flit->getElement()->getType() == ElementBase::DEGRADER) {
                //check if previous element: is on, is a degrader, ends where new element starts
                if ((*fprev).isOn() && fprev->getElement()->getType() == ElementBase::DEGRADER &&
                    ((*fprev).getEnd() + 0.01 > (*flit).getStart()) ) {
                    (*flit).setOn(kineticEnergy);
                }
            }
gsell's avatar
gsell committed
123
        }
124

kraus's avatar
kraus committed
125
        fprev = flit;
gsell's avatar
gsell committed
126 127 128
    }
}

129 130
void OpalBeamline::switchElementsOff(const double &min, ElementBase::ElementType eltype) {
    if(eltype == ElementBase::ANY) {
gsell's avatar
gsell committed
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151
        for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) {
            if((*flit).isOn() && min >= (*flit).getEnd()) {
                (*flit).setOff();
            }

        }
    } else {
        for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit) {
            if((*flit).isOn() && min >= (*flit).getEnd() && (*flit).getElement()->getType() == eltype) {
                (*flit).setOff();
            }
        }
    }
}

void OpalBeamline::switchElementsOff() {
    for(FieldList::iterator flit = elements_m.begin(); flit != elements_m.end(); ++ flit)
        (*flit).setOff();
}

void OpalBeamline::prepareSections() {
152 153 154 155
    if (elements_m.size() == 0) {
        prepared_m = true;
        return;
    }
156 157 158 159 160
    prepared_m = true;
}

void OpalBeamline::print(Inform &msg) const {
}
161

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
void OpalBeamline::swap(OpalBeamline & rhs) {
    std::swap(elements_m, rhs.elements_m);
    std::swap(prepared_m, rhs.prepared_m);
    std::swap(coordTransformationTo_m, rhs.coordTransformationTo_m);
}

void OpalBeamline::merge(OpalBeamline &rhs) {
    elements_m.insert(elements_m.end(),
                      rhs.elements_m.begin(),
                      rhs.elements_m.end());
    prepared_m = false;
    containsSource_m = containsSource_m || rhs.containsSource_m;
}


FieldList OpalBeamline::getElementByType(ElementBase::ElementType type) {
    if (type == ElementBase::ANY) {
        return elements_m;
gsell's avatar
gsell committed
180 181
    }

182 183 184 185
    FieldList elements_of_requested_type;
    for(FieldList::iterator fit = elements_m.begin(); fit != elements_m.end(); ++ fit) {
        if((*fit).getElement()->getType() == type) {
            elements_of_requested_type.push_back((*fit));
gsell's avatar
gsell committed
186 187
        }
    }
188 189
    return elements_of_requested_type;
}
gsell's avatar
gsell committed
190

191 192 193 194 195 196 197 198 199 200
void OpalBeamline::compute3DLattice() {
    static unsigned int order = 0;
    const FieldList::iterator end = elements_m.end();

    unsigned int minOrder = order;
    {
        double endPriorPathLength = 0.0;
        CoordinateSystemTrafo currentCoordTrafo = coordTransformationTo_m;

        elements_m.sort([](const ClassicField& a, const ClassicField& b) {
201 202 203 204 205 206
                double edgeA = 0.0, edgeB = 0.0;
                if (a.getElement()->isElementPositionSet())
                    edgeA = a.getElement()->getElementPosition();

                if (b.getElement()->isElementPositionSet())
                    edgeB = b.getElement()->getElementPosition();
207 208 209

                return edgeA < edgeB;
            });
210
        FieldList::iterator it = elements_m.begin();
211 212 213 214 215 216 217
        for (; it != end; ++ it) {
            if ((*it).isPositioned()) {
                continue;
            }
            (*it).order_m = minOrder;
            std::shared_ptr<Component> element = (*it).getElement();

218 219 220 221 222 223 224 225 226 227 228 229 230
            if (element->getType() != ElementBase::SBEND &&
                element->getType() != ElementBase::RBEND &&
                element->getType() != ElementBase::RBEND3D) {
                continue;
            }

            double beginThisPathLength = element->getElementPosition();
            Vector_t beginThis3D(0, 0, beginThisPathLength - endPriorPathLength);
            BendBase * bendElement = static_cast<BendBase*>(element.get());
            double thisLength = bendElement->getChordLength();
            double bendAngle = bendElement->getBendAngle();
            double entranceAngle = bendElement->getEntranceAngle();
            double arcLength = (thisLength * std::abs(bendAngle) / (2 * sin(std::abs(bendAngle) / 2)));
231

232 233 234
            double rotationAngleAboutZ = bendElement->getRotationAboutZ();
            Quaternion_t rotationAboutZ(cos(0.5 * rotationAngleAboutZ),
                                        sin(-0.5 * rotationAngleAboutZ) * Vector_t(0, 0, 1));
235

236 237
            Vector_t effectiveRotationAxis = rotationAboutZ.rotate(Vector_t(0, -1, 0));
            effectiveRotationAxis /= euclidean_norm(effectiveRotationAxis);
238

239 240 241 242 243 244
            Quaternion_t rotationAboutAxis(cos(0.5 * bendAngle),
                                           sin(0.5 * bendAngle) * effectiveRotationAxis);
            Quaternion_t halfRotationAboutAxis(cos(0.25 * bendAngle),
                                               sin(0.25 * bendAngle) * effectiveRotationAxis);
            Quaternion_t entryFaceRotation(cos(0.5 * entranceAngle),
                                           sin(0.5 * entranceAngle) * effectiveRotationAxis);
245

246 247 248 249 250 251 252 253 254
            if (!Options::idealized) {
                std::vector<Vector_t> truePath = bendElement->getDesignPath();
                Quaternion_t directionExitHardEdge(cos(0.5 * (0.5 * bendAngle - entranceAngle)),
                                                   sin(0.5 * (0.5 * bendAngle - entranceAngle)) * effectiveRotationAxis);
                Vector_t exitHardEdge = thisLength * directionExitHardEdge.rotate(Vector_t(0, 0, 1));
                double distanceEntryHETruePath = euclidean_norm(truePath.front());
                double distanceExitHETruePath = euclidean_norm(rotationAboutZ.rotate(truePath.back()) - exitHardEdge);
                double pathLengthTruePath = (*it).getEnd() - (*it).getStart();
                arcLength = pathLengthTruePath - distanceEntryHETruePath - distanceExitHETruePath;
gsell's avatar
gsell committed
255
            }
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270

            Vector_t chord = thisLength * halfRotationAboutAxis.rotate(Vector_t(0, 0, 1));
            Vector_t endThis3D = beginThis3D + chord;
            double endThisPathLength = beginThisPathLength + arcLength;

            CoordinateSystemTrafo fromEndLastToBeginThis(beginThis3D,
                                                         (entryFaceRotation * rotationAboutZ).conjugate());
            CoordinateSystemTrafo fromEndLastToEndThis(endThis3D,
                                                       rotationAboutAxis.conjugate());

            (*it).setCoordTransformationTo(fromEndLastToBeginThis * currentCoordTrafo);

            currentCoordTrafo = (fromEndLastToEndThis * currentCoordTrafo);

            endPriorPathLength = endThisPathLength;
gsell's avatar
gsell committed
271 272 273
        }
    }

274 275 276
    double endPriorPathLength = 0.0;
    CoordinateSystemTrafo currentCoordTrafo = coordTransformationTo_m;

277
    FieldList::iterator it = elements_m.begin();
278 279 280 281 282 283
    for (; it != end; ++ it) {
        if ((*it).isPositioned()) continue;

        (*it).order_m = order ++;

        std::shared_ptr<Component> element = (*it).getElement();
284
        double beginThisPathLength = element->getElementPosition();
285 286 287
        double thisLength = element->getElementLength();
        Vector_t beginThis3D(0, 0, beginThisPathLength - endPriorPathLength);

288 289 290
        if (element->getType() == ElementBase::SOURCE) {
            beginThis3D(2) -= thisLength;
        }
291 292 293 294 295 296 297 298 299 300 301 302

        Vector_t endThis3D;
        if (element->getType() == ElementBase::SBEND ||
            element->getType() == ElementBase::RBEND ||
            element->getType() == ElementBase::RBEND3D) {

            BendBase * bendElement = static_cast<BendBase*>(element.get());
            thisLength = bendElement->getChordLength();
            double bendAngle = bendElement->getBendAngle();

            double rotationAngleAboutZ = bendElement->getRotationAboutZ();
            Quaternion_t rotationAboutZ(cos(0.5 * rotationAngleAboutZ),
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
303 304 305
                                        sin(-0.5 * rotationAngleAboutZ) * Vector_t(0, 0, 1));

            Vector_t effectiveRotationAxis = rotationAboutZ.rotate(Vector_t(0, -1, 0));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
306
            effectiveRotationAxis /= euclidean_norm(effectiveRotationAxis);
307 308

            Quaternion_t rotationAboutAxis(cos(0.5 * bendAngle),
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
309
                                           sin(0.5 * bendAngle) * effectiveRotationAxis);
310
            Quaternion halfRotationAboutAxis(cos(0.25 * bendAngle),
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
311
                                             sin(0.25 * bendAngle) * effectiveRotationAxis);
312 313 314 315 316 317 318

            double arcLength = (thisLength * std::abs(bendAngle) /
                                (2 * sin(bendAngle / 2)));
            if (!Options::idealized) {
                std::vector<Vector_t> truePath = bendElement->getDesignPath();
                double entranceAngle = bendElement->getEntranceAngle();
                Quaternion_t directionExitHardEdge(cos(0.5 * (0.5 * bendAngle - entranceAngle)),
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
319
                                                   sin(0.5 * (0.5 * bendAngle - entranceAngle)) * effectiveRotationAxis);
320
                Vector_t exitHardEdge = thisLength * directionExitHardEdge.rotate(Vector_t(0, 0, 1));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
321
                double distanceEntryHETruePath = euclidean_norm(truePath.front());
322
                double distanceExitHETruePath = euclidean_norm(rotationAboutZ.rotate(truePath.back()) - exitHardEdge);
323 324
                double pathLengthTruePath = (*it).getEnd() - (*it).getStart();
                arcLength = pathLengthTruePath - distanceEntryHETruePath - distanceExitHETruePath;
gsell's avatar
gsell committed
325
            }
326 327 328 329 330 331 332 333 334 335 336 337 338

            endThis3D = (beginThis3D +
                         halfRotationAboutAxis.rotate(Vector_t(0, 0, thisLength)));
            CoordinateSystemTrafo fromEndLastToEndThis(endThis3D,
                                                       rotationAboutAxis.conjugate());
            currentCoordTrafo = fromEndLastToEndThis * currentCoordTrafo;

            endPriorPathLength = beginThisPathLength + arcLength;
        } else {
            double rotationAngleAboutZ = (*it).getElement()->getRotationAboutZ();
            Quaternion_t rotationAboutZ(cos(0.5 * rotationAngleAboutZ),
                                        sin(-0.5 * rotationAngleAboutZ) * Vector_t(0, 0, 1));

339
            CoordinateSystemTrafo fromLastToThis(beginThis3D, rotationAboutZ);
340 341

            (*it).setCoordTransformationTo(fromLastToThis * currentCoordTrafo);
gsell's avatar
gsell committed
342
        }
343 344 345 346 347 348 349 350

        (*it).fixPosition();
    }

    elements_m.sort(ClassicField::SortAsc);
}

void OpalBeamline::save3DLattice() {
351
    if (Ippl::myNode() != 0 || OpalData::getInstance()->isOptimizerRun()) return;
352 353 354 355 356 357 358 359 360

    elements_m.sort([](const ClassicField& a, const ClassicField& b) {
            return a.order_m < b.order_m;
        });

    FieldList::iterator it = elements_m.begin();
    FieldList::iterator end = elements_m.end();

    std::ofstream pos;
361
    std::string fileName = "data/" + OpalData::getInstance()->getInputBasename() + "_ElementPositions.txt";
362 363
    if (OpalData::getInstance()->getOpenMode() == OpalData::OPENMODE::APPEND &&
        boost::filesystem::exists(fileName)) {
364 365 366 367 368 369 370 371
        pos.open(fileName, std::ios_base::app);
    } else {
        pos.open(fileName);
    }

    MeshGenerator mesh;
    for (; it != end; ++ it) {
        std::shared_ptr<Component> element = (*it).getElement();
372 373 374
        CoordinateSystemTrafo toBegin = element->getEdgeToBegin() * (*it).getCoordTransformationTo();
        CoordinateSystemTrafo toEnd = element->getEdgeToEnd() * (*it).getCoordTransformationTo();
        Vector_t entry3D = toBegin.getOrigin();
375 376 377 378 379 380 381
        Vector_t exit3D = toEnd.getOrigin();

        mesh.add(*(element.get()));

        if (element->getType() == ElementBase::SBEND ||
            element->getType() == ElementBase::RBEND) {

382
            Bend2D * bendElement = static_cast<Bend2D*>(element.get());
383
            std::vector<Vector_t> designPath = bendElement->getDesignPath();
384 385
            toEnd = bendElement->getBeginToEnd_local() * (*it).getCoordTransformationTo();
            exit3D = toEnd.getOrigin();
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406

            unsigned int size = designPath.size();
            unsigned int minNumSteps = std::max(20.0,
                                                std::abs(bendElement->getBendAngle() / Physics::pi * 180));
            unsigned int frequency = std::floor((double)size / minNumSteps);

            pos << std::setw(30) << std::left << std::string("\"ENTRY EDGE: ") + element->getName() + std::string("\"")
                << std::setw(18) << std::setprecision(10) << entry3D(2)
                << std::setw(18) << std::setprecision(10) << entry3D(0)
                << std::setw(18) << std::setprecision(10) << entry3D(1)
                << "\n";

            Vector_t position = (*it).getCoordTransformationTo().transformFrom(designPath.front());
            pos << std::setw(30) << std::left << std::string("\"BEGIN: ") + element->getName() + std::string("\"")
                << std::setw(18) << std::setprecision(10) << position(2)
                << std::setw(18) << std::setprecision(10) << position(0)
                << std::setw(18) << std::setprecision(10) << position(1)
                << std::endl;

            for (unsigned int i = frequency; i + 1 < size; i += frequency) {

407
                position = (*it).getCoordTransformationTo().transformFrom(designPath[i]);
408 409 410 411
                pos << std::setw(30) << std::left << std::string("\"MID: ") + element->getName() + std::string("\"")
                    << std::setw(18) << std::setprecision(10) << position(2)
                    << std::setw(18) << std::setprecision(10) << position(0)
                    << std::setw(18) << std::setprecision(10) << position(1)
412
                    << std::endl;
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
            }

            position = (*it).getCoordTransformationTo().transformFrom(designPath.back());
            pos << std::setw(30) << std::left << std::string("\"END: ") + element->getName() + std::string("\"")
                << std::setw(18) << std::setprecision(10) << position(2)
                << std::setw(18) << std::setprecision(10) << position(0)
                << std::setw(18) << std::setprecision(10) << position(1)
                << std::endl;

            pos << std::setw(30) << std::left << std::string("\"EXIT EDGE: ") + element->getName() + std::string("\"")
                << std::setw(18) << std::setprecision(10) << exit3D(2)
                << std::setw(18) << std::setprecision(10) << exit3D(0)
                << std::setw(18) << std::setprecision(10) << exit3D(1)
                << std::endl;
        } else {
            pos << std::setw(30) << std::left << std::string("\"BEGIN: ") + element->getName() + std::string("\"")
                << std::setw(18) << std::setprecision(10) << entry3D(2)
                << std::setw(18) << std::setprecision(10) << entry3D(0)
                << std::setw(18) << std::setprecision(10) << entry3D(1)
                << "\n";

            pos << std::setw(30) << std::left << std::string("\"END: ") + element->getName() + std::string("\"")
                << std::setw(18) << std::setprecision(10) << exit3D(2)
                << std::setw(18) << std::setprecision(10) << exit3D(0)
                << std::setw(18) << std::setprecision(10) << exit3D(1)
                << std::endl;
        }
gsell's avatar
gsell committed
440
    }
441 442
    elements_m.sort(ClassicField::SortAsc);
    mesh.write(OpalData::getInstance()->getInputBasename());
gsell's avatar
gsell committed
443 444
}

445 446 447 448 449 450 451 452 453 454 455 456 457 458
namespace {
    std::string parseInput() {

        std::ifstream in(OpalData::getInstance()->getInputFn());
        std::string source("");
        std::string str;
        char testBit;
        const std::string commentFormat("");
        const boost::regex empty("^[ \t]*$");
        const boost::regex lineEnd(";");
        const std::string lineEndFormat(";\n");
        const boost::regex cppCommentExpr("//.*");
        const boost::regex cCommentExpr("/\\*.*?\\*/"); // "/\\*(?>[^*/]+|\\*[^/]|/[^*])*(?>(?R)(?>[^*/]+|\\*[^/]|/[^*])*)*\\*/"
        bool priorEmpty = true;
459 460

        in.get(testBit);
461 462 463 464 465 466 467 468 469 470 471 472 473
        while (!in.eof()) {
            in.putback(testBit);

            std::getline(in, str);
            str = boost::regex_replace(str, cppCommentExpr, commentFormat);
            str = boost::regex_replace(str, empty, commentFormat);
            if (str.size() > 0) {
                source += str;// + '\n';
                priorEmpty = false;
            } else if (!priorEmpty) {
                source += "##EMPTY_LINE##";
                priorEmpty = true;
            }
gsell's avatar
gsell committed
474

475 476
            in.get(testBit);
        }
477

478 479
        source = boost::regex_replace(source, cCommentExpr, commentFormat);
        source = boost::regex_replace(source, lineEnd, lineEndFormat, boost::match_default | boost::format_all);
480

481
        return source;
gsell's avatar
gsell committed
482
    }
483

484 485 486
    unsigned int getMinimalSignificantDigits(double num, const unsigned int maxDigits) {
        char buf[32];
        snprintf(buf, 32, "%.*f", maxDigits + 1, num);
487
        std::string numStr(buf);
488 489 490 491 492 493 494 495
        unsigned int length = numStr.length();

        unsigned int numDigits = maxDigits;
        unsigned int i = 2;
        while (i < maxDigits + 1 && numStr[length - i] == '0') {
            --numDigits;
            ++ i;
        }
496

497 498
        return numDigits;
    }
499

500 501
    std::string round2string(double num, const unsigned int maxDigits) {
        char buf[64];
502

503 504 505 506
        snprintf(buf, 64, "%.*f", getMinimalSignificantDigits(num, maxDigits), num);

        return std::string(buf);
    }
gsell's avatar
gsell committed
507 508
}

509
void OpalBeamline::save3DInput() {
510
    if (Ippl::myNode() != 0 || OpalData::getInstance()->isOptimizerRun()) return;
511 512 513 514 515

    FieldList::iterator it = elements_m.begin();
    FieldList::iterator end = elements_m.end();

    std::string input = parseInput();
516
    std::ofstream pos("data/" + OpalData::getInstance()->getInputBasename() + "_3D.opal");
517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533

    for (; it != end; ++ it) {
        std::string element = (*it).getElement()->getName();
        const boost::regex replacePSI("(" + element + "\\s*:[^\\n]*)PSI\\s*=[^,;]*,?", boost::regex::icase);
        input = boost::regex_replace(input, replacePSI, "\\1\\2");

        const boost::regex replaceELEMEDGE("(" + element + "\\s*:[^\\n]*)ELEMEDGE\\s*=[^,;]*(.)", boost::regex::icase);

        CoordinateSystemTrafo cst = (*it).getCoordTransformationTo();
        Vector_t origin = cst.getOrigin();
        Vector_t orient = Util::getTaitBryantAngles(cst.getRotation().conjugate(), element);
        for (unsigned int d = 0; d < 3; ++ d)
            orient(d) *= Physics::rad2deg;

        std::string x = (std::abs(origin(0)) > 1e-10? "X = " + round2string(origin(0), 10) + ", ": "");
        std::string y = (std::abs(origin(1)) > 1e-10? "Y = " + round2string(origin(1), 10) + ", ": "");
        std::string z = (std::abs(origin(2)) > 1e-10? "Z = " + round2string(origin(2), 10) + ", ": "");
gsell's avatar
gsell committed
534

535 536 537 538 539 540 541 542 543 544 545 546 547 548
        std::string theta = (orient(0) > 1e-10? "THETA = " + round2string(orient(0), 6) + " * PI / 180, ": "");
        std::string phi = (orient(1) > 1e-10? "PHI = " + round2string(orient(1), 6) + " * PI / 180, ": "");
        std::string psi = (orient(2) > 1e-10? "PSI = " + round2string(orient(2), 6) + " * PI / 180, ": "");
        std::string coordTrafo = x + y + z + theta + phi + psi;
        if (coordTrafo.length() > 2) {
            coordTrafo = coordTrafo.substr(0, coordTrafo.length() - 2); // remove last ', '
        }

        std::string position = ("\\1" + coordTrafo + "\\2");

        input = boost::regex_replace(input, replaceELEMEDGE, position);

        if ((*it).getElement()->getType() == ElementBase::RBEND ||
            (*it).getElement()->getType() == ElementBase::SBEND) {
549
            const Bend2D* dipole = static_cast<const Bend2D*>((*it).getElement().get());
550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567
            double angle = dipole->getBendAngle();
            double E1 = dipole->getEntranceAngle();
            double E2 = dipole->getExitAngle();

            const boost::regex angleR("(" + element + "\\s*:[^\\n]*ANGLE\\s*=)[^,;]*(.)");
            const std::string angleF("\\1 " + round2string(angle * 180 / Physics::pi, 6) + " / 180 * PI\\2");
            const boost::regex E1R("(" + element + "\\s*:[^\\n]*E1\\s*=)[^,;]*(.)");
            const std::string E1F("\\1 " + round2string(E1 * 180 / Physics::pi, 6) + " / 180 * PI\\2");
            const boost::regex E2R("(" + element + "\\s*:[^\\n]*E2\\s*=)[^,;]*(.)");
            const std::string E2F("\\1 " + round2string(E2 * 180 / Physics::pi, 6) + " / 180 * PI\\2");
            const boost::regex noRotation("(" + element + "\\s*:[^\\n]*),\\s*ROTATION\\s*=[^,;]*(.)");
            const std::string noRotationFormat("\\1\\2  ");

            input = boost::regex_replace(input, angleR, angleF);
            input = boost::regex_replace(input, E1R, E1F);
            input = boost::regex_replace(input, E2R, E2F);
            input = boost::regex_replace(input, noRotation, noRotationFormat);
        }
gsell's avatar
gsell committed
568
    }
569 570 571 572 573 574

    const boost::regex empty("##EMPTY_LINE##");
    const std::string emptyFormat("\n");
    input = boost::regex_replace(input, empty, emptyFormat);

    pos << input << std::endl;
gsell's avatar
gsell committed
575 576
}

577 578 579 580 581 582 583 584 585
void OpalBeamline::activateElements() {
    auto it = elements_m.begin();
    const auto end = elements_m.end();

    double designEnergy = 0.0;
    for (; it != end; ++ it) {
        std::shared_ptr<Component> element = (*it).getElement();
        if (element->getType() == ElementBase::SBEND ||
            element->getType() == ElementBase::RBEND) {
586
            Bend2D * bendElement = static_cast<Bend2D*>(element.get());
587 588 589 590 591
            designEnergy = bendElement->getDesignEnergy() * 1e-6;
        }
        (*it).setOn(designEnergy);
        // element->goOnline(designEnergy);
    }
592
}