ClosedOrbitFinder.h 35.2 KB
Newer Older
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
//
// Class ClosedOrbitFinder
//   This class finds a closed orbit of a cyclotron for a given energy.
//   The algorithm is based on the paper of M. M. Gordon: "Computation of
//   closed orbits and basic focusing properties for sector-focused cyclotrons
//   and the design of 'cyclops'" (1983)
//   As template arguments one chooses the type of the variables and the
//   integrator for the ODEs. The supported steppers can be found on
//   http://www.boost.org/ where it is part of the library Odeint.
//
// Copyright (c) 2014, Matthias Frey, ETH Zürich
// All rights reserved
//
// Implemented as part of the Semester thesis by Matthias Frey
// "Matched Distributions in Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
27 28 29
#ifndef CLOSEDORBITFINDER_H
#define CLOSEDORBITFINDER_H

30
#include <algorithm>
31 32 33
#include <array>
#include <cmath>
#include <functional>
adelmann's avatar
adelmann committed
34
#include <limits>
35
#include <numeric>
adelmann's avatar
adelmann committed
36
#include <string>
37
#include <utility>
38 39
#include <vector>

40
#include "Utilities/Options.h"
41 42
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
43
#include "Physics/Physics.h"
44

frey_m's avatar
frey_m committed
45 46
#include "AbstractObjects/OpalData.h"

47
#include "AbsBeamline/Cyclotron.h"
frey_m's avatar
frey_m committed
48

49 50
// include headers for integration
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
frey_m's avatar
frey_m committed
51 52 53
#include <boost/filesystem.hpp>

extern Inform *gmsg;
54 55 56 57

template<typename Value_type, typename Size_type, class Stepper>
class ClosedOrbitFinder
{
58 59 60 61 62 63 64 65 66
    public:
        /// Type of variables
        typedef Value_type value_type;
        /// Type for specifying sizes
        typedef Size_type size_type;
        /// Type of container for storing quantities (path length, orbit, etc.)
        typedef std::vector<value_type> container_type;
        /// Type for holding state of ODE values
        typedef std::vector<value_type> state_type;
snuverink_j's avatar
snuverink_j committed
67

frey_m's avatar
frey_m committed
68
        typedef std::function<void(const state_type&, state_type&, const double)> function_t;
69 70 71

        /// Sets the initial values for the integration and calls findOrbit().
        /*!
72
         * @param E0 is the potential energy (particle energy at rest) [MeV].
adelmann's avatar
adelmann committed
73
         * @param N specifies the number of splits (2pi/N), i.e number of integration steps
74
         * @param cycl is the cyclotron element
75 76
         * @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector,
         * otherwise over 2*pi.
77 78
         * @param Nsectors is an int (default: 1). Number of sectors that the field map is averaged over
         * in order to avoid first harmonic. Only valid if domain is false
79
         */
frey_m's avatar
frey_m committed
80
        ClosedOrbitFinder(value_type E0, size_type N,
81
                          Cyclotron* cycl, bool domain = true, int Nsectors = 1);
82 83

        /// Returns the inverse bending radius (size of container N+1)
84
        container_type getInverseBendingRadius(const value_type& angle = 0);
85 86

        /// Returns the step lengths of the path (size of container N+1)
87
        container_type getPathLength(const value_type& angle = 0);
88 89

        /// Returns the field index (size of container N+1)
90
        container_type getFieldIndex(const value_type& angle = 0);
91 92 93 94

        /// Returns the radial and vertical tunes (in that order)
        std::pair<value_type,value_type> getTunes();

95 96 97 98 99 100 101
        /// Returns the closed orbit (size of container N+1) starting at specific angle (only makes sense when computing
        /// the closed orbit for a whole turn) (default value: 0°).
        /// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
        /// part, i.e. it takes the next starting index below. There's no interpolation of the radius.
        /*!
         * @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
         */
102 103
        container_type getOrbit(value_type angle = 0);

104 105 106 107 108 109
        /// Returns the momentum of the orbit (size of container N+1)starting at specific angle (only makes sense when
        /// computing the closed orbit for a whole turn) (default value: 0°), \f$ \left[ p_{r} \right] = \si{m}\f$.
        /// Attention: It computes the starting index of the array. If it's not an integer it just cuts the floating point
        /// part, i.e. it takes the next starting index below. There's no interpolation of the momentum.
        /*!
         * @param angle is the start angle for the output. Has to be within [0°,360°[ (default: 0°).
110
         * @returns the momentum in \f$ \beta * \gamma \f$ units
111
         */
112
        container_type getMomentum(value_type angle = 0);
113 114 115 116

        /// Returns the average orbit radius
        value_type getAverageRadius();

adelmann's avatar
adelmann committed
117 118
        /// Returns the frequency error
        value_type getFrequencyError();
119 120 121 122 123 124 125 126

        /// Returns true if a closed orbit could be found
        bool isConverged();

        /// Computes the closed orbit
        /*!
         * @param accuracy specifies the accuracy of the closed orbit
         * @param maxit is the maximal number of iterations done for finding the closed orbit
frey_m's avatar
frey_m committed
127
         * @param ekin energy for which to find closed orbit (in tune mode: upper limit of range)
frey_m's avatar
frey_m committed
128
         * @param dE step increase [MeV]
129
         * @param rguess initial radius guess in [mm]
frey_m's avatar
frey_m committed
130
         * @param isTuneMode comptute tunes of all energies in one sweep
131
         */
132 133 134 135
        bool findOrbit(value_type accuracy, size_type maxit,
                       value_type ekin,
                       value_type dE = 1.0, value_type rguess = -1.0,
                       bool isTuneMode = false);
136

137
        /// Fills in the values of h_m, ds_m, fidx_m.
138
        void computeOrbitProperties(const value_type& E);
139

140
    private:
141 142 143 144 145 146 147 148
        /// This function is called by the function getTunes().
        /*! Transfer matrix Y = [y11, y12; y21, y22] (see Gordon paper for more details).
         * @param y are the positions (elements y11 and y12 of Y)
         * @param py2 is the momentum of the second solution (element y22 of Y)
         * @param ncross is the number of sign changes (\#crossings of zero-line)
         */
        value_type computeTune(const std::array<value_type,2>&, value_type, size_type);

frey_m's avatar
frey_m committed
149 150 151
        // Compute closed orbit for given energy
        bool findOrbitOfEnergy_m(const value_type&, container_type&, value_type&,
                                 const value_type&, size_type);
snuverink_j's avatar
snuverink_j committed
152

adelmann's avatar
adelmann committed
153
        /// This function computes nzcross_ which is used to compute the tune in z-direction and the frequency error
154
//         void computeVerticalOscillations();
snuverink_j's avatar
snuverink_j committed
155

156 157
        /// This function rotates the calculated closed orbit finder properties to the initial angle
        container_type rotate(value_type angle, container_type& orbitProperty);
158 159 160 161 162

        /// Stores current position in horizontal direction for the solutions of the ODE with different initial values
        std::array<value_type,2> x_m; // x_m = [x1, x2]
        /// Stores current momenta in horizontal direction for the solutions of the ODE with different initial values
        std::array<value_type,2> px_m; // px_m = [px1, px2]
frey_m's avatar
frey_m committed
163
        /// Stores current position in vertical direction for the solutions of the ODE with different initial values
164
        std::array<value_type,2> z_m; // z_m = [z1, z2]
frey_m's avatar
frey_m committed
165
        /// Stores current momenta in vertical direction for the solutions of the ODE with different initial values
166 167 168 169 170 171 172 173
        std::array<value_type,2> pz_m; // pz_m = [pz1, pz2]

        /// Stores the inverse bending radius
        container_type h_m;
        /// Stores the step length
        container_type ds_m;
        /// Stores the radial orbit (size: N_m+1)
        container_type r_m;
174 175
        /// Stores the vertical oribt (size: N_m+1)
        container_type vz_m;
176 177
        /// Stores the radial momentum
        container_type pr_m;
178 179
        /// Stores the vertical momentum
        container_type vpz_m;
180 181 182 183 184 185 186 187
        /// Stores the field index
        container_type fidx_m;

        /// Counts the number of zero-line crossings in horizontal direction (used for computing horizontal tune)
        size_type nxcross_m;
        /// Counts the number of zero-line crossings in vertical direction (used for computing vertical tune)
        size_type nzcross_m; //#crossings of zero-line in x- and z-direction

188
        /// Is the rest mass [MeV / c**2]
189
        value_type E0_m;
snuverink_j's avatar
snuverink_j committed
190

191
        /// Is the nominal orbital frequency
192 193 194
        /* (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal
         * Coupling by Space Charge in Cyclotrons" (2012), formula (1))
         */
195
        value_type wo_m;
adelmann's avatar
adelmann committed
196
        /// Number of integration steps
197 198 199 200 201 202 203 204 205 206
        size_type N_m;
        /// Is the angle step size
        value_type dtheta_m;

        /// Is the average radius
        value_type ravg_m;

        /// Is the phase
        value_type phase_m;

207
        /**
208 209 210 211
         * Stores the last orbit value (since we have to return to the beginning to check the convergence in the
         * findOrbit() function. This last value is then deleted from the array but is stored in lastOrbitVal_m to
         * compute the vertical oscillations)
         */
snuverink_j's avatar
snuverink_j committed
212
        /* value_type lastOrbitVal_m; */
213

214 215 216 217 218
        /**
         * Stores the last momentum value (since we have to return to the beginning to check the convergence in the
         * findOrbit() function. This last value is then deleted from the array but is stored in lastMomentumVal_m to
         * compute the vertical oscillations)
         */
snuverink_j's avatar
snuverink_j committed
219
        /* value_type lastMomentumVal_m; */
220 221

        /**
222 223 224
         * Boolean which is true by default. "true": orbit integration over one sector only, "false": integration
         * over 2*pi
         */
adelmann's avatar
adelmann committed
225
        bool domain_m;
226 227 228 229 230
        /**
         * Number of sectors to average the field map over
         * in order to avoid first harmonic. Only valid if domain is false
         */
        int  nSectors_m;
231

232 233
        /// Defines the stepper for integration of the ODE's
        Stepper stepper_m;
snuverink_j's avatar
snuverink_j committed
234

235
        /*!
snuverink_j's avatar
snuverink_j committed
236
         * This quantity is defined in the paper "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons"
237 238 239 240
         * of Dr. Christian Baumgarten (2012)
         * The lambda function takes the orbital frequency \f$ \omega_{o} \f$ (also defined in paper) as input argument.
         */
        std::function<double(double)> acon_m = [](double wo) { return Physics::c / wo; };
snuverink_j's avatar
snuverink_j committed
241

242 243 244 245 246 247 248
        /// Cyclotron unit \f$ \left[T\right] \f$ (Tesla)
        /*!
         * The lambda function takes the orbital frequency \f$ \omega_{o} \f$ as input argument.
         */
        std::function<double(double, double)> bcon_m = [](double e0, double wo) {
            return e0 * 1.0e7 / (/* physics::q0 */ 1.0 * Physics::c * Physics::c / wo);
        };
snuverink_j's avatar
snuverink_j committed
249

250
        Cyclotron* cycl_m;
251 252 253 254 255 256
};

// -----------------------------------------------------------------------------------------------------------------------
// PUBLIC MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------

257 258 259
template<typename Value_type, typename Size_type, class Stepper>
ClosedOrbitFinder<Value_type,
                  Size_type,
frey_m's avatar
frey_m committed
260
                  Stepper>::ClosedOrbitFinder(value_type E0,
261
                                              size_type N, Cyclotron* cycl,
262
                                              bool domain, int nSectors)
263 264 265
    : nxcross_m(0)
    , nzcross_m(0)
    , E0_m(E0)
266
    , wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
267 268 269 270
    , N_m(N)
    , dtheta_m(Physics::two_pi/value_type(N))
    , ravg_m(0)
    , phase_m(0)
snuverink_j's avatar
snuverink_j committed
271 272
    /* , lastOrbitVal_m(0.0) */
    /* , lastMomentumVal_m(0.0) */
273
    , domain_m(domain)
274
    , nSectors_m(nSectors)
275
    , stepper_m()
276
    , cycl_m(cycl)
277
{
snuverink_j's avatar
snuverink_j committed
278

279
    if ( cycl_m->getFMLowE() > cycl_m->getFMHighE() )
frey_m's avatar
frey_m committed
280 281
        throw OpalException("ClosedOrbitFinder::ClosedOrbitFinder()",
                            "Incorrect cyclotron energy (MeV) bounds: Maximum cyclotron energy smaller than minimum cyclotron energy.");
snuverink_j's avatar
snuverink_j committed
282

adelmann's avatar
adelmann committed
283 284
    // if domain_m = true --> integrate over a single sector
    if (domain_m) {
285
        N_m /=  cycl_m->getSymmetry();
adelmann's avatar
adelmann committed
286
    }
snuverink_j's avatar
snuverink_j committed
287

frey_m's avatar
frey_m committed
288 289
    cycl_m->read(cycl_m->getFieldFlag(cycl_m->getCyclotronType()),
                 cycl_m->getBScale());
290

291 292 293 294 295
    // reserve storage for the orbit and momentum (--> size = 0, capacity = N_m+1)
    /*
     * we need N+1 storage, since dtheta = 2pi/N (and not 2pi/(N-1)) that's why we need N+1 integration steps
     * to return to the origin (but the return size is N_m)
     */
adelmann's avatar
adelmann committed
296 297
    r_m.reserve(N_m + 1);
    pr_m.reserve(N_m + 1);
298 299
    vz_m.reserve(N_m + 1);
    vpz_m.reserve(N_m + 1);
300

301
    // reserve memory of N_m for the properties (--> size = 0, capacity = N_m)
adelmann's avatar
adelmann committed
302 303 304
    h_m.reserve(N_m);
    ds_m.reserve(N_m);
    fidx_m.reserve(N_m);
305 306 307
}

template<typename Value_type, typename Size_type, class Stepper>
308
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
309
    ClosedOrbitFinder<Value_type, Size_type, Stepper>::getInverseBendingRadius(const value_type& angle)
310
{
311 312 313 314
    if (angle != 0.0)
        return rotate(angle, h_m);
    else
        return h_m;
315 316 317
}

template<typename Value_type, typename Size_type, class Stepper>
318
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
319
    ClosedOrbitFinder<Value_type, Size_type, Stepper>::getPathLength(const value_type& angle)
320
{
321 322 323 324
    if (angle != 0.0)
        return rotate(angle, ds_m);
    else
        return ds_m;
325 326 327
}

template<typename Value_type, typename Size_type, class Stepper>
328
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
329
    ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFieldIndex(const value_type& angle)
330
{
331 332
    if (angle != 0.0)
        return rotate(angle, fidx_m);
frey_m's avatar
frey_m committed
333
    return fidx_m;
334 335 336
}

template<typename Value_type, typename Size_type, class Stepper>
337 338 339
std::pair<Value_type,Value_type> ClosedOrbitFinder<Value_type, Size_type, Stepper>::getTunes() {
    // compute radial tune
    value_type nur = computeTune(x_m,px_m[1],nxcross_m);
340

341 342 343 344
    // compute vertical tune
    value_type nuz = computeTune(z_m,pz_m[1],nzcross_m);

    return std::make_pair(nur,nuz);
345 346 347
}

template<typename Value_type, typename Size_type, class Stepper>
348
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
349
    ClosedOrbitFinder<Value_type, Size_type, Stepper>::getOrbit(value_type angle)
350
{
351 352 353 354
    if (angle != 0.0)
        return rotate(angle, r_m);
    else
        return r_m;
355 356 357 358
}

template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
359
    ClosedOrbitFinder<Value_type, Size_type, Stepper>::getMomentum(value_type angle)
360 361
{
    container_type pr = pr_m;
snuverink_j's avatar
snuverink_j committed
362

363 364
    if (angle != 0.0)
        pr = rotate(angle, pr);
snuverink_j's avatar
snuverink_j committed
365

366 367
    // change units from meters to \beta * \gamma
    /* in Gordon paper:
snuverink_j's avatar
snuverink_j committed
368
     *
369
     * p = \gamma * \beta * a
snuverink_j's avatar
snuverink_j committed
370
     *
371
     * where a = c / \omega_{0} with \omega_{0} = 2 * \pi * \nu_{0} = 2 * \pi * \nu_{rf} / h
snuverink_j's avatar
snuverink_j committed
372
     *
373 374 375
     * c: speed of light
     * h: harmonic number
     * v_{rf}: nomial rf frequency
snuverink_j's avatar
snuverink_j committed
376
     *
377
     * Units:
snuverink_j's avatar
snuverink_j committed
378
     *
379
     * [a] = m --> [p] = m
snuverink_j's avatar
snuverink_j committed
380
     *
381
     * The momentum in \beta * \gamma is obtained by dividing by "a"
382
     */
383
    value_type factor =  1.0 / acon_m(wo_m);
384
    std::for_each(pr.begin(), pr.end(), [factor](value_type& p) { p *= factor; });
snuverink_j's avatar
snuverink_j committed
385

386
    return pr;
387 388 389
}

template<typename Value_type, typename Size_type, class Stepper>
390 391 392
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
    ClosedOrbitFinder<Value_type, Size_type, Stepper>::getAverageRadius()
{
393
    return ravg_m;
394 395 396
}

template<typename Value_type, typename Size_type, class Stepper>
397 398
typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type
    ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFrequencyError()
399
{
400
    return phase_m;
401 402 403 404 405 406
}

// -----------------------------------------------------------------------------------------------------------------------
// PRIVATE MEMBER FUNCTIONS
// -----------------------------------------------------------------------------------------------------------------------

frey_m's avatar
frey_m committed
407 408


409
template<typename Value_type, typename Size_type, class Stepper>
410 411
bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type accuracy,
                                                                  size_type maxit,
frey_m's avatar
frey_m committed
412
                                                                  value_type ekin,
frey_m's avatar
frey_m committed
413
                                                                  value_type dE,
frey_m's avatar
frey_m committed
414 415
                                                                  value_type rguess,
                                                                  bool isTuneMode)
416
{
417 418 419 420
    /* REMARK TO GORDON
     * q' = 1/b = 1/bcon
     * a' = a = acon
     */
421

422 423 424
    // resize vectors (--> size = N_m+1, capacity = N_m+1), note: we do N_m+1 integration steps
    r_m.resize(N_m+1);
    pr_m.resize(N_m+1);
425 426
    vz_m.resize(N_m+1);
    vpz_m.resize(N_m+1);
427

428 429
    // store acon locally
    value_type acon = acon_m(wo_m);            // [acon] = m
430 431
    // amplitude of error; Gordon, formula (18) (a = a')
    value_type error = std::numeric_limits<value_type>::max();
432 433 434 435

    /*
     * Christian:
     * N = 1440 ---> N = 720 ---> dtheta = 2PI/720 --> nsteps = 721
436
     *
437
     * 0, 2, 4, ... ---> jeden zweiten berechnene: 1, 3, 5, ... interpolieren --> 1440 Werte
438
     *
439 440
     * Matthias:
     * N = 1440 --> dtheta = 2PI/1440 --> nsteps = 1441
441
     *
442
     * 0, 1, 2, 3, 4, 5, ... --> 1440 Werte
443
     *
444
     */
445

frey_m's avatar
frey_m committed
446 447 448
    value_type E         = cycl_m->getFMLowE(); // starting energy
    value_type E_fin     = ekin;                // final    energy
    const value_type eps = dE * 1.0e-1;         // articial constant for stopping criteria
449

frey_m's avatar
frey_m committed
450
    if (isTuneMode) {
frey_m's avatar
frey_m committed
451
        E_fin = cycl_m->getFMHighE();
frey_m's avatar
frey_m committed
452
    }
Andreas Adelmann's avatar
Andreas Adelmann committed
453

frey_m's avatar
frey_m committed
454 455
    namespace fs = boost::filesystem;
    fs::path dir = OpalData::getInstance()->getInputBasename();
456
    dir = dir.parent_path() / OpalData::getInstance()->getAuxiliaryOutputDirectory();
frey_m's avatar
frey_m committed
457
    std::string tunefile = (dir / "tunes.dat").string();
458

frey_m's avatar
frey_m committed
459 460
    if ( isTuneMode ) {
        std::ofstream out(tunefile, std::ios::out);
461

frey_m's avatar
frey_m committed
462
        out << std::left
463
            << std::setw(15) << "# energy[MeV]"
464
            << std::setw(15) << "radius_ini[m]"
frey_m's avatar
frey_m committed
465
            << std::setw(15) << "momentum_ini[Beta Gamma]"
466
            << std::setw(15) << "radius_avg[m]"
frey_m's avatar
frey_m committed
467 468 469 470
            << std::setw(15) << "nu_r"
            << std::setw(15) << "nu_z"
            << std::endl;
    }
471 472 473 474
    // initial guess
    container_type init;
    enum Guess {NONE, FIRST, SECOND};
    Guess guess = NONE;
frey_m's avatar
frey_m committed
475 476
    value_type rn1 = 0.0, pn1 = 0.0; // normalised r, pr value of previous closed orbit
    value_type rn2 = 0.0, pn2 = 0.0; // normalised r, pr value of second to previous closed orbit
477

478 479
    // iterate until suggested energy (start with minimum energy)
    // increase energy by dE
frey_m's avatar
frey_m committed
480 481 482 483 484
    *gmsg << level3 << "Start iteration to find closed orbit of energy " << E_fin << " MeV "
          << "in steps of " << dE << " MeV." << endl;

    for (; E <= E_fin + eps; E += dE) {

frey_m's avatar
frey_m committed
485
        error = std::numeric_limits<value_type>::max();
486 487 488 489 490 491

        // energy dependent values
        value_type en     = E / E0_m;        // en = E/E0 = E/(mc^2) (E0 is potential energy)
        value_type gamma  = en + 1.0;
        value_type gamma2 = gamma * gamma;   // = gamma^2
        value_type beta   = std::sqrt(1.0 - 1.0 / gamma2);
492
        value_type p = acon_m(wo_m) * std::sqrt(en * (2.0 + en));  // momentum [p] = m; Gordon, formula (3)
493 494 495 496 497 498 499 500 501 502 503 504 505 506

        if (guess == NONE) {
            // set initial values for radius and radial momentum for lowest energy Emin
            // orbit, [r] = m;  Gordon, formula (20)
            // radial momentum; Gordon, formula (20)
            //      r            pr   z    pz
            init = {beta * acon, 0.0, 0.0, 1.0};
            if (rguess >= 0.0) {
                init[0] = rguess * 0.001;
            }
            guess = FIRST;
        } else if (guess == FIRST) {
            // new initial values based on previous one, formula (21)
            init[0] = (beta*acon) * rn1;
507
            init[1] = p*pn1;
508 509 510 511
            guess = SECOND;
        } else if (guess == SECOND) {
            // second extrapolation, formula (21)
            init[0] = (beta*acon) * (rn1 + (rn1-rn2));
512
            init[1] = p*(pn1 + (pn1-pn2));
513 514 515 516 517 518 519
        }

        std::fill(  r_m.begin(),   r_m.end(), 0);
        std::fill( pr_m.begin(),  pr_m.end(), 0);
        std::fill( vz_m.begin(),  vz_m.end(), 0);
        std::fill(vpz_m.begin(), vpz_m.end(), 0);

520
        // (re-)set inital values for r and pr
521 522 523
        r_m[0]   = init[0];
        pr_m[0]  = init[1];
        vz_m[0]  = init[2];
524
        vpz_m[0] = init[3];
525

frey_m's avatar
frey_m committed
526 527
        *gmsg << level3 << "    Try to find orbit for " << E << " MeV ... ";

frey_m's avatar
frey_m committed
528
        if ( !this->findOrbitOfEnergy_m(E, init, error, accuracy, maxit) ) {
frey_m's avatar
frey_m committed
529
            *gmsg << endl << "ClosedOrbitFinder didn't converge for energy " + std::to_string(E) + " MeV." << endl;
530
            guess = NONE;
531
            continue;
frey_m's avatar
frey_m committed
532
        }
frey_m's avatar
frey_m committed
533 534 535

        *gmsg << level3 << "Successfully found." << endl;

536 537
        // store for next initial guess
        rn2 = rn1;
538 539 540
        pn2 = pn1;
        rn1 = r_m[0] / (acon * beta);
        pn1 = pr_m[0] / p;
541

frey_m's avatar
frey_m committed
542
        if ( isTuneMode ) {
543

544 545 546
            this->computeOrbitProperties(E);
            value_type reo = this->getOrbit(   cycl_m->getPHIinit())[0];
            value_type peo = this->getMomentum(cycl_m->getPHIinit())[0];
frey_m's avatar
frey_m committed
547 548
            std::pair<value_type , value_type > tunes = this->getTunes();

549 550
            *gmsg << std::left
                  << "* ----------------------------" << endl
frey_m's avatar
frey_m committed
551 552
                  << "* Closed orbit info (Gordon units):" << endl
                  << "*" << endl
553 554 555 556
                  << "* kinetic energy:   " << std::setw(12) << E      << " [MeV]" << endl
                  << "* average radius:   " << std::setw(12) << ravg_m << " [m]" << endl
                  << "* initial radius:   " << std::setw(12) << reo    << " [m]" << endl
                  << "* initial momentum: " << std::setw(12) << peo    << " [Beta Gamma]" << endl
frey_m's avatar
frey_m committed
557 558 559 560
                  << "* frequency error:  " << phase_m        << endl
                  << "* horizontal tune:  " << tunes.first    << endl
                  << "* vertical tune:    " << tunes.second   << endl
                  << "* ----------------------------" << endl << endl;
561

frey_m's avatar
frey_m committed
562 563 564 565
            std::ofstream out(tunefile, std::ios::app);
            out << std::left
                << std::setw(15) << E
                << std::setw(15) << reo
frey_m's avatar
frey_m committed
566
                << std::setw(15) << peo
567
                << std::setw(15) << ravg_m
frey_m's avatar
frey_m committed
568 569 570 571
                << std::setw(15) << tunes.first
                << std::setw(15) << tunes.second << std::endl;
            out.close();
        }
572
    }
573

frey_m's avatar
frey_m committed
574 575
    *gmsg << level3 << "Finished closed orbit finder successfully." << endl;

576 577 578
    /* store last entry, since it is needed in computeVerticalOscillations(), because we have to do the same
     * number of integrations steps there.
     */
snuverink_j's avatar
snuverink_j committed
579 580
    /* lastOrbitVal_m    = r_m[N_m];        // needed in computeVerticalOscillations() */
    /* lastMomentumVal_m = pr_m[N_m];       // needed in computeVerticalOscillations() */
581

582 583 584
    // remove last entry (since we don't have to store [0,2pi], but [0,2pi[)  --> size = N_m, capacity = N_m+1
    r_m.pop_back();
    pr_m.pop_back();
585

586 587 588 589 590
    /* domain_m = true --> only integrated over a single sector
     * --> multiply by cycl_m->getSymmetry() to get correct phase_m
     */
    if (domain_m)
        phase_m *= cycl_m->getSymmetry();
591

592 593 594
    // returns true if converged, otherwise false
    return error < accuracy;
}
595

frey_m's avatar
frey_m committed
596 597 598 599 600 601 602 603
template<typename Value_type, typename Size_type, class Stepper>
bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbitOfEnergy_m(
    const value_type& E,
    container_type& init,
    value_type& error,
    const value_type& accuracy,
    size_type maxit)
{
604 605 606
    /* *gmsg << "rguess : " << init[0] << endl; */
    /* *gmsg << "prguess: " << init[1] << endl; */

frey_m's avatar
frey_m committed
607 608
    value_type bint, brint, btint;
    value_type invbcon = 1.0 / bcon_m(E0_m, wo_m);      // [bcon] = MeV*s/(C*m^2) = 10^6 T = 10^7 kG (kilo Gauss)
snuverink_j's avatar
snuverink_j committed
609

frey_m's avatar
frey_m committed
610 611 612 613 614 615
    value_type xold = 0.0;                              // for counting nxcross
    value_type zold = 0.0;                              // for counting nzcross

    // index for reaching next element of the arrays r and pr (no nicer way found yet)
    size_type idx = 0;
    // observer for storing the current value after each ODE step (e.g. Runge-Kutta step) into the containers of r and pr
616
    auto store = [&](state_type& y, const value_type /*t*/)
frey_m's avatar
frey_m committed
617
    {
618 619 620
        r_m[idx]   = y[0];
        pr_m[idx]  = y[1];
        vz_m[idx]  = y[6];
frey_m's avatar
frey_m committed
621 622 623 624 625
        vpz_m[idx] = y[7];

        // count number of crossings (excluding starting point --> idx>0)
        nxcross_m += (idx > 0) * (y[4] * xold < 0);
        xold = y[4];
snuverink_j's avatar
snuverink_j committed
626

frey_m's avatar
frey_m committed
627 628 629
        // number of times z2 changes sign
        nzcross_m += (idx > 0) * (y[10] * zold < 0);
        zold = y[10];
snuverink_j's avatar
snuverink_j committed
630

frey_m's avatar
frey_m committed
631 632
        ++idx;
    };
snuverink_j's avatar
snuverink_j committed
633

frey_m's avatar
frey_m committed
634 635 636 637
    // define initial state container for integration: y = {r, pr, x1, px1, x2, px2,
    //                                                      z, pz, z1, pz1, z2, pz2,
    //                                                      phase}
    state_type y(11);
snuverink_j's avatar
snuverink_j committed
638

frey_m's avatar
frey_m committed
639 640 641 642 643 644
    // difference of last and first value of r (1. element) and pr (2. element)
    container_type err(2);
    // correction term for initial values: r = r + dr, pr = pr + dpr; Gordon, formula (17)
    container_type delta = {0.0, 0.0};
    // if niterations > maxit --> stop iteration
    size_type niterations = 0;
snuverink_j's avatar
snuverink_j committed
645

frey_m's avatar
frey_m committed
646
    // energy dependent values
647 648
    value_type en     = E / E0_m;                      // en = E/E0 = E/(mc^2) (E0 is potential energy)
    value_type gamma  = en + 1.0;
frey_m's avatar
frey_m committed
649
    value_type p = acon_m(wo_m) * std::sqrt(en * (2.0 + en));     // momentum [p] = m; Gordon, formula (3)
650 651
    value_type gamma2    = gamma * gamma;           // = gamma^2
    value_type invgamma4 = 1.0 / (gamma2 * gamma2); // = 1/gamma^4
652 653
    value_type p2 = p * p;                              // p^2 = p*p

frey_m's avatar
frey_m committed
654 655
    // helper constants
    value_type pr2;                                     // squared radial momentum (pr^2 = pr*pr)
656
    value_type ptheta, invptheta;                       // azimuthal momentum
frey_m's avatar
frey_m committed
657 658 659 660 661 662 663

    // define the six ODEs (using lambda function)
    function_t orbit_integration = [&](const state_type &y,
                                       state_type &dydt,
                                       const double theta)
    {
        pr2 = y[1] * y[1];
664 665
        if (p2 < pr2) {
            //*gmsg << theta << " " << p2 << " " << pr2 << endl;
666
            throw OpalException("ClosedOrbitFinder::findOrbitOfEnergy_m()",
frey_m's avatar
frey_m committed
667
                                "p_{r}^2 > p^{2} (defined in Gordon paper) --> Square root of negative number.");
668
        }
frey_m's avatar
frey_m committed
669
        // Gordon, formula (5c)
670
        ptheta    = std::sqrt(p2 - pr2);
frey_m's avatar
frey_m committed
671
        invptheta = 1.0 / ptheta;
672 673 674 675 676 677 678 679 680 681 682 683 684 685 686
        // average field over the number of sectors
        brint=0.0, btint=0.0, bint=0.0;
        for (int i = 0; i<nSectors_m; i++) {
            double angle = theta + i * Physics::two_pi / nSectors_m;
            double tmpbr, tmpbt, tmpb;
            // interpolate values of magnetic field
            cycl_m->apply(y[0], y[6], angle, tmpbr, tmpbt, tmpb);
            brint += tmpbr;
            btint += tmpbt;
            bint  += tmpb;
        }
        brint /= nSectors_m;
        btint /= nSectors_m;
        bint  /= nSectors_m;
        // multiply by 1 / b
snuverink_j's avatar
snuverink_j committed
687
        bint  *= invbcon;
frey_m's avatar
frey_m committed
688 689 690 691 692
        brint *= invbcon;
        btint *= invbcon;

        // Gordon, formula (5a)
        dydt[0] = y[0] * y[1] * invptheta;
693
        // Gordon, formula (5b) (typo in paper! second equal sign is a minus)
frey_m's avatar
frey_m committed
694 695 696
        dydt[1] = ptheta - y[0] * bint;
        // Gordon, formulas (9a) and (9b)
        for (size_type i = 2; i < 5; i += 2) {
697
            dydt[i]   = (y[1] * y[i] + y[0] * p2 * y[i+1] * invptheta * invptheta) * invptheta;
frey_m's avatar
frey_m committed
698 699
            dydt[i+1] = - y[1] * y[i+1] * invptheta - (bint + y[0] * brint) * y[i];
        }
700

frey_m's avatar
frey_m committed
701 702
        // Gordon, formulas (22a) and (22b)
        for (size_type i = 6; i < 12; i += 2) {
703
            dydt[i]   = y[0] * y[i+1] * invptheta;
frey_m's avatar
frey_m committed
704 705 706 707
            dydt[i+1] = (y[0] * brint - y[1] * invptheta * btint) * y[i];
        }

        // integrate phase
708
        dydt[12] = y[0] * invptheta * gamma - 1;
709

frey_m's avatar
frey_m committed
710 711
    };

712
    // integrate until error smaller than user-defined accuracy
frey_m's avatar
frey_m committed
713
    do {
714
        // (re-)set initial values
715
        x_m[0]  = 1.0;               // x1;  Gordon, formula (10)
frey_m's avatar
frey_m committed
716
        px_m[0] = 0.0;               // px1; Gordon, formula (10)
717
        x_m[1]  = 0.0;               // x2;  Gordon, formula (10)
frey_m's avatar
frey_m committed
718 719 720 721 722 723 724 725 726
        px_m[1] = 1.0;               // px2; Gordon, formula (10)
        z_m[0]  = 1.0;
        pz_m[0] = 0.0;
        z_m[1]  = 0.0;
        pz_m[1] = 1.0;
        phase_m = 0.0;
        nxcross_m = 0;               // counts the number of crossings of x-axis (excluding first step)
        nzcross_m = 0;
        idx = 0;                     // index for looping over r and pr arrays
727

frey_m's avatar
frey_m committed
728 729 730 731 732 733 734 735
        // fill container with initial states
        y = {init[0],init[1],
             x_m[0], px_m[0], x_m[1], px_m[1],
             init[2], init[3],
             z_m[0], pz_m[0], z_m[1], pz_m[1],
             phase_m
        };

736
        try {
737
            // integrate from 0 to 2*pi / nSectors (one has to get back to the "origin")
738
            boost::numeric::odeint::integrate_n_steps(stepper_m, orbit_integration,y,0.0,dtheta_m,N_m,store);
739 740
        } catch(OpalException & ex) {
            *gmsg << ex.where() << " " << ex.what() << endl;
741 742 743
            break;
        }

frey_m's avatar
frey_m committed
744
        // write new state
745
        x_m[0]  = y[2];
frey_m's avatar
frey_m committed
746
        px_m[0] = y[3];
747
        x_m[1]  = y[4];
frey_m's avatar
frey_m committed
748
        px_m[1] = y[5];
749 750

        z_m[0]  = y[8];
frey_m's avatar
frey_m committed
751
        pz_m[0] = y[9];
752
        z_m[1]  = y[10];
frey_m's avatar
frey_m committed
753 754
        pz_m[1] = y[11];
        phase_m = y[12] * Physics::u_two_pi; // / (2.0 * Physics::pi);
755

frey_m's avatar
frey_m committed
756 757
        // compute error (compare values of orbit and momentum for theta = 0 and theta = 2*pi)
        // (Note: size = N_m+1 --> last entry is N_m)
758
        err[0] =  r_m[N_m] -  r_m[0];    // Gordon, formula (14)
frey_m's avatar
frey_m committed
759 760
        err[1] = pr_m[N_m] - pr_m[0];    // Gordon, formula (14)

761
        // correct initial values of r and pr
frey_m's avatar
frey_m committed
762
        value_type invdenom = 1.0 / (x_m[0] + px_m[1] - 2.0);
763 764
        delta[0] = ((px_m[1] - 1.0) * err[0] -  x_m[1] * err[1]) * invdenom; // dr;  Gordon, formula (16a)
        delta[1] = (( x_m[0] - 1.0) * err[1] - px_m[0] * err[0]) * invdenom; // dpr; Gordon, formula (16b)
frey_m's avatar
frey_m committed
765 766 767 768

        // improved initial values; Gordon, formula (17) (here it's used for higher energies)
        init[0] += delta[0];
        init[1] += delta[1];
769 770

        // compute amplitude of the error (Gordon, formula (18)
frey_m's avatar
frey_m committed
771
        error = std::sqrt(delta[0] * delta[0] + delta[1] * delta[1] * invgamma4) / r_m[0];
772 773 774 775 776

        // *gmsg << "iteration " << niterations << " error: " << error << endl;

    } while ((error > accuracy) && (niterations++ < maxit));

777
    if (error > accuracy)
778
        *gmsg << "findOrbit not converged after " << niterations << " iterations with error: " << error << ". Needed accuracy " << accuracy << endl;
779

frey_m's avatar
frey_m committed
780 781 782
    return (error < accuracy);
}

783
template<typename Value_type, typename Size_type, class Stepper>
784 785 786
Value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeTune(const std::array<value_type,2>& y,
                                                                          value_type py2, size_type ncross)
{
787
    // Y = [y1, y2; py1, py2]
788

789 790
    // cos(mu)
    value_type cos = 0.5 * (y[0] + py2);
snuverink_j's avatar
snuverink_j committed
791

792
    value_type mu;
793

794 795
    // sign of sin(mu) has to be equal to y2
    bool negative = std::signbit(y[1]);
796

797
    bool uneven = (ncross % 2);
798

799
    value_type abscos = std::abs(cos);
800 801
    value_type muPrime;
    if (abscos > 1.0) {
802 803
        // store the number of crossings
        if (uneven)
804
            ncross = ncross + 1;
805

806
        // Gordon, formula (36b)
807
        muPrime = -std::acosh(abscos);    // mu'
808

809
    } else {
810
        muPrime = (uneven) ? std::acos(-cos) : std::acos(cos);    // mu'
811
        /* It has to be fulfilled: 0<= mu' <= pi
812 813 814 815 816
         * But since |cos(mu)| <= 1, we have
         * -1 <= cos(mu) <= 1 --> 0 <= mu <= pi (using above programmed line), such
         * that condition is already fulfilled.
         */
    }
817

818 819
    // Gordon, formula (36)
    mu = ncross * Physics::pi + muPrime;
820

821
    if (abscos < 1.0) {
822 823
        // if sign(y[1]) > 0 && sign(sin(mu)) < 0
        if (!negative && std::signbit(std::sin(mu))) {
824
            mu = ncross * Physics::pi - muPrime;
825
        } else if (negative && !std::signbit(std::sin(mu))) {
826
            mu = ncross * Physics::pi - muPrime + Physics::two_pi;
827 828
        }
    }
829

830
    // nu = mu/theta, where theta = integration domain
831

832
    /* domain_m = true --> only integrated over a single sector --> multiply by cycl_m->getSymmetry() to
adelmann's avatar
adelmann committed
833 834 835
     * get correct tune.
     */
    if (domain_m)
836
        mu *= cycl_m->getSymmetry();
837

838
    return mu * Physics::u_two_pi;
839 840 841
}

template<typename Value_type, typename Size_type, class Stepper>
842
void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties(const value_type& E) {
843
    /*
844 845 846 847 848
     * The formulas for h, fidx and ds are from the paper:
     * "Tranverse-Longitudinal Coupling by Space Charge in Cyclotrons"
     * written by Dr. Christian Baumgarten (2012, PSI)
     * p. 6
     */
849

adelmann's avatar
adelmann committed
850
    // READ IN MAGNETIC FIELD: ONLY FOR STAND-ALONE PROGRAM
851
    value_type bint, brint, btint; // B, dB/dr, dB/dtheta
852

853
    value_type invbcon = 1.0 / bcon_m(E0_m, wo_m);
854
    value_type en = E / E0_m;                                  // en = E/E0 = E/(mc^2) (E0 is potential energy)
855
    value_type p = acon_m(wo_m) * std::sqrt(en * (2.0 + en));  // momentum [p] = m; Gordon, formula (3)
856
    value_type p2 = p * p;
857
    value_type theta = 0.0;                                    // angle for interpolating
858
    value_type ptheta;
859

860 861 862 863 864 865 866
    // resize of container (--> size = N_m, capacity = N_m)
    h_m.resize(N_m);
    fidx_m.resize(N_m);
    ds_m.resize(N_m);

    for (size_type i = 0; i < N_m; ++i) {
        // interpolate magnetic field
867
        cycl_m->apply(r_m[i], vz_m[i], theta, brint, btint, bint);
snuverink_j's avatar
snuverink_j committed
868
        bint  *= invbcon;
869 870
        brint *= invbcon;
        btint *= invbcon;
snuverink_j's avatar
snuverink_j committed
871

872 873 874 875
        // inverse bending radius
        h_m[i] = bint / p;

        // local field index
876
        if (p < pr_m[i])
snuverink_j's avatar
snuverink_j committed
877
            throw OpalException("ClosedOrbitFinder::computeOrbitProperties()",
878 879
                                "p_{r}^2 > p^{2} " + std::to_string(p) + " " + std::to_string(pr_m[i]) + " (defined in Gordon paper) --> Square root of negative number.");

880
        ptheta = std::sqrt(p2 - pr_m[i] * pr_m[i]);
881

882 883 884 885 886 887 888
        fidx_m[i] = (brint * ptheta - btint * pr_m[i] / r_m[i]) / p2; //(bint*bint);

        // path length element
        ds_m[i] = std::hypot(r_m[i] * pr_m[i] / ptheta,r_m[i]) * dtheta_m; // C++11 function

        // increase angle
        theta += dtheta_m;
889
    }
890 891 892

    // compute average radius
    ravg_m = std::accumulate(r_m.begin(),r_m.end(),0.0) / value_type(r_m.size());
893 894
}

snuverink_j's avatar
snuverink_j committed
895
template<typename Value_type, typename Size_type, class Stepper>
896 897 898 899
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type
ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, container_type &orbitProperty) {

    container_type orbitPropertyCopy = orbitProperty;
snuverink_j's avatar
snuverink_j committed
900

901 902 903 904 905 906 907 908
    // compute the number of steps per degree
    value_type deg_step = N_m / 360.0;

    // compute starting point
    size_type start = deg_step * angle;

    // copy end to start
    std::copy(orbitProperty.begin() + start, orbitProperty.end(), orbitPropertyCopy.begin());
snuverink_j's avatar
snuverink_j committed
909

910 911 912 913 914 915 916
    // copy start to end
    std::copy_n(orbitProperty.begin(), start, orbitPropertyCopy.end() - start);

    return orbitPropertyCopy;

}

917
#endif