ClosedOrbitFinder.h 28.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
/**
 * @file ClosedOrbitFinder.h
 * 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.
 *
 * @author Matthias Frey
 * @version 1.0
 */

12 13 14 15 16 17
#ifndef CLOSEDORBITFINDER_H
#define CLOSEDORBITFINDER_H

#include <array>
#include <cmath>
#include <functional>
adelmann's avatar
adelmann committed
18
#include <limits>
19
#include <numeric>
adelmann's avatar
adelmann committed
20
#include <stdexcept>
adelmann's avatar
adelmann committed
21
#include <string>
22
#include <utility>
23 24 25 26
#include <vector>

#include "physics.h"

adelmann's avatar
adelmann committed
27
#include "MagneticField.h" // ONLY FOR STAND-ALONE PROGRAM
28 29 30 31 32 33 34


#include <fstream>

// include headers for integration
#include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>

35
/// @brief Finds a closed orbit of a cyclotron for a given energy
36 37 38
template<typename Value_type, typename Size_type, class Stepper>
class ClosedOrbitFinder
{
39 40 41 42 43 44 45 46 47 48 49 50 51 52
    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;

        /// Sets the initial values for the integration and calls findOrbit().
        /*!
         * @param E is the energy [MeV] to which the closed orbit should be found
         * @param wo is the nominal orbital frequency (see paper of Dr. C. Baumgarten: "Transverse-Longitudinal Coupling by Space Charge in Cyclotrons" (2012), formula (1))
adelmann's avatar
adelmann committed
53
         * @param N specifies the number of splits (2pi/N), i.e number of integration steps
54 55 56 57
         * @param accuracy specifies the accuracy of the closed orbit
         * @param maxit is the maximal number of iterations done. Program stops if closed orbit not found within this time.
         * @param Emin is the minimum energy [MeV] needed in cyclotron
         * @param Emax is the maximum energy [MeV] reached in cyclotron
adelmann's avatar
adelmann committed
58 59 60 61 62
         * @param nSector is the number of sectors (--> symmetry) of cyclotron
         * @param rmin is the minimal radius of the cyclotron, \f$ \left[r_{min}\right] = m \f$
         * @param ntheta is the number of angle splits (fieldmap variable)
         * @param nradial is the number of radial splits (fieldmap variable)
         * @param dr is the radial step size, \f$ \left[\Delta r\right] = m \f$
63
         * @param fieldmap is the location of the file that specifies the magnetic field
adelmann's avatar
adelmann committed
64
         * @param domain is a boolean (default: true). If "true" the closed orbit is computed over a single sector, otherwise over 2*pi.
65
         */
kraus's avatar
kraus committed
66
        ClosedOrbitFinder(value_type, value_type, size_type, value_type, size_type, value_type, value_type, size_type, value_type, size_type, size_type, value_type, const std::string&, bool = true);
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88

        /// Returns the inverse bending radius (size of container N+1)
        container_type& getInverseBendingRadius();

        /// Returns the step lengths of the path (size of container N+1)
        container_type& getPathLength();

        /// Returns the field index (size of container N+1)
        container_type& getFieldIndex();

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

        /// Returns the closed orbit (size of container N+1)
        container_type& getOrbit();

        /// Returns the relativistic factor gamma
        value_type getGamma();

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

adelmann's avatar
adelmann committed
89 90
        /// Returns the frequency error
        value_type getFrequencyError();
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113

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

    private:
        /// 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
         */
        bool findOrbit(value_type, size_type);

        /// Fills in the values of h_m, ds_m, fidx_m. It gets called by in by constructor.
        void computeOrbitProperties();

        /// 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);

adelmann's avatar
adelmann committed
114
        /// This function computes nzcross_ which is used to compute the tune in z-direction and the frequency error
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
        void computeVerticalOscillations();

        /// 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]
        /// Stores current position in longitudinal direction for the solutions of the ODE with different initial values
        std::array<value_type,2> z_m; // z_m = [z1, z2]
        /// Stores current momenta in longitudinal direction for the solutions of the ODE with different initial values
        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;
        /// Stores the radial momentum
        container_type pr_m;
        /// 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

        /// Is the energy for which the closed orbit should be found
        value_type E_m;
        /// Is the nominal orbital frequency
        value_type wo_m;
adelmann's avatar
adelmann committed
146
        /// Number of integration steps
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
        size_type N_m;
        /// Is the angle step size
        value_type dtheta_m;

        /// Is the relativistic factor
        value_type gamma_m;

        /// Is the average radius
        value_type ravg_m;

        /// Is the phase
        value_type phase_m;

        /// Boolean which tells if a closed orbit for this configuration could be found (get set by the function findOrbit)
        bool converged_m;

        /// Minimum energy needed in cyclotron
        value_type Emin_m;

        /// Maximum energy reached in cyclotron
        value_type Emax_m;
kraus's avatar
kraus committed
168

adelmann's avatar
adelmann committed
169 170
        /// Number of sectors (symmetry)
        size_type nSector_m;
kraus's avatar
kraus committed
171

adelmann's avatar
adelmann committed
172 173
        /// Minimal radius of cyclotron, \f$ \left[r_{min}\right] = m \f$
        value_type rmin_m;
kraus's avatar
kraus committed
174

adelmann's avatar
adelmann committed
175 176
        /// Number of angle splits (fieldmap)
        size_type ntheta_m;
kraus's avatar
kraus committed
177

adelmann's avatar
adelmann committed
178 179
        /// Number of radial steps (fieldmap)
        size_type nradial_m;
kraus's avatar
kraus committed
180

adelmann's avatar
adelmann committed
181 182
        /// Radial step size, \f$ \left[\Delta r\right] = m \f$
        value_type dr_m;
183 184 185 186 187 188

        /// 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)
        value_type lastOrbitVal_m;

        /// 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)
        value_type lastMomentumVal_m;
kraus's avatar
kraus committed
189

adelmann's avatar
adelmann committed
190
        /// Boolean which is true if computeVerticalOscillations() executed, otherwise false. This is used for checking in getTunes() and getFrequencyError().
191 192 193 194
        bool vertOscDone_m;

        /// Location of magnetic field
        std::string fieldmap_m;
kraus's avatar
kraus committed
195

adelmann's avatar
adelmann committed
196 197
        /// Boolean which is true by default. "true": orbit integration over one sector only, "false": integration over 2*pi
        bool domain_m;
198

199 200 201 202 203
        /// Defines the stepper for integration of the ODE's
        Stepper stepper_m;

        /// ONLY FOR STAND-ALONE PROGRAM
        float** bmag_m;
204 205 206 207 208 209
};

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

210
    template<typename Value_type, typename Size_type, class Stepper>
211
ClosedOrbitFinder<Value_type, Size_type, Stepper>::ClosedOrbitFinder(value_type E, value_type wo, size_type N, value_type accuracy,
kraus's avatar
kraus committed
212
        size_type maxit,value_type Emin, value_type Emax, size_type nSector, value_type rmin, size_type ntheta, size_type nradial, value_type dr, const std::string& fieldmap, bool domain)
adelmann's avatar
adelmann committed
213 214 215 216
    : nxcross_m(0), nzcross_m(0), E_m(E), wo_m(wo), N_m(N), dtheta_m(2.0*M_PI/value_type(N)),
    gamma_m(E/physics::E0+1.0), ravg_m(0), phase_m(0), converged_m(false), Emin_m(Emin), Emax_m(Emax), nSector_m(nSector),
    rmin_m(rmin), ntheta_m(ntheta), nradial_m(nradial), dr_m(dr), lastOrbitVal_m(0.0), lastMomentumVal_m(0.0),
    vertOscDone_m(false), fieldmap_m(fieldmap), domain_m(domain), stepper_m()
217
{
218 219

    if (Emin_m > Emax_m || E_m < Emin_m || E > Emax_m)
adelmann's avatar
adelmann committed
220
        throw std::domain_error("Error in ClosedOrbitFinder: Emin <= E <= Emax and Emin < Emax");
kraus's avatar
kraus committed
221

adelmann's avatar
adelmann committed
222 223 224
    // velocity: beta = v/c = sqrt(1-1/(gamma*gamma))
    if (gamma_m == 0)
        throw std::invalid_argument("Error in ClosedOrbitFinder: Relativistic factor equal zero.");
kraus's avatar
kraus committed
225

adelmann's avatar
adelmann committed
226 227 228 229
    // if domain_m = true --> integrate over a single sector
    if (domain_m) {
        N_m /=  nSector_m;
    }
kraus's avatar
kraus committed
230

231 232 233 234 235
    // 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
236 237
    r_m.reserve(N_m + 1);
    pr_m.reserve(N_m + 1);
kraus's avatar
kraus committed
238

239
    // reserve memory of N_m for the properties (--> size = 0, capacity = N_m)
adelmann's avatar
adelmann committed
240 241 242
    h_m.reserve(N_m);
    ds_m.reserve(N_m);
    fidx_m.reserve(N_m);
kraus's avatar
kraus committed
243

244
    // compute closed orbit
245
    converged_m = findOrbit(accuracy, maxit);
kraus's avatar
kraus committed
246

247 248 249 250 251 252
    // compute h, ds, fidx, rav (average radius)
    computeOrbitProperties();
}

template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type& ClosedOrbitFinder<Value_type, Size_type, Stepper>::getInverseBendingRadius() {
253
    return h_m;
254 255 256 257
}

template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type& ClosedOrbitFinder<Value_type, Size_type, Stepper>::getPathLength() {
258
    return ds_m;
259 260 261 262
}

template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type& ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFieldIndex() {
263
    return fidx_m;
264 265 266
}

template<typename Value_type, typename Size_type, class Stepper>
267 268 269
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);
kraus's avatar
kraus committed
270

271 272 273
    // compute nzcross_m
    if (!vertOscDone_m)
        computeVerticalOscillations();
kraus's avatar
kraus committed
274

275 276 277 278
    // compute vertical tune
    value_type nuz = computeTune(z_m,pz_m[1],nzcross_m);

    return std::make_pair(nur,nuz);
279 280 281 282
}

template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_type& ClosedOrbitFinder<Value_type, Size_type, Stepper>::getOrbit() {
283
    return r_m;
284 285 286 287
}

template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::getGamma() {
288
    return gamma_m;
289 290 291 292
}

template<typename Value_type, typename Size_type, class Stepper>
inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::getAverageRadius() {
293
    return ravg_m;
294 295 296
}

template<typename Value_type, typename Size_type, class Stepper>
adelmann's avatar
adelmann committed
297
typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::getFrequencyError() {
kraus's avatar
kraus committed
298

299 300 301
    // if the vertical oscillations aren't computed, we have to, since there we also compuote the frequency error.
    if(!vertOscDone_m)
        computeVerticalOscillations();
kraus's avatar
kraus committed
302

303
    return phase_m;
304 305 306 307
}

template<typename Value_type, typename Size_type, class Stepper>
inline bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::isConverged() {
308
    return converged_m;
kraus's avatar
kraus committed
309
}
310 311 312 313 314 315 316

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

template<typename Value_type, typename Size_type, class Stepper>
bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type accuracy, size_type maxit) {
317 318 319 320 321
    /* REMARK TO GORDON
     * q' = 1/b = 1/bcon
     * a' = a = acon
     */

adelmann's avatar
adelmann committed
322 323 324 325
    // READ IN MAGNETIC FIELD: ONLY FOR STAND-ALONE PROGRAM
    bmag_m = MagneticField::malloc2df(ntheta_m,nradial_m);
    MagneticField::ReadSectorMap(bmag_m,nradial_m,ntheta_m,1,fieldmap_m,0.0);
    MagneticField::MakeNFoldSymmetric(bmag_m,ntheta_m,nradial_m,ntheta_m/nSector_m,nSector_m);
326
    value_type bint, brint, btint;
kraus's avatar
kraus committed
327

328 329 330
    // 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);
kraus's avatar
kraus committed
331

332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
    // store acon and bcon locally
    value_type acon = physics::acon(wo_m);               // [acon] = m
    value_type invbcon = 1.0 / physics::bcon(wo_m);        // [bcon] = MeV*s/(C*m^2) = 10^6 T = 10^7 kG (kilo Gauss)

    // helper constants
    value_type p2;                                      // p^2 = p*p
    value_type pr2;                                     // squared radial momentum (pr^2 = pr*pr)
    value_type ptheta, invptheta;                       // Gordon, formula (5c)
    value_type invdenom;                                // denominator for computing dr,dpr
    value_type xold = 0.0;                              // for counting nxcross

    // 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
    auto store = [&](state_type& y, const value_type t)
kraus's avatar
kraus committed
347
    {
348 349 350 351 352 353
        r_m[idx] = y[0];
        pr_m[idx] = y[1];

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

355 356 357 358 359 360 361
        ++idx;
    };

    // define the six ODEs (using lambda function)
    std::function<void(const state_type&, state_type&, const double)> orbit_integration = [&](const state_type &y, state_type &dydt, const double theta){
        pr2 = y[1] * y[1];
        if (p2 < pr2)
adelmann's avatar
adelmann committed
362
            throw std::domain_error("Error in ClosedOrbitFinder::findOrbit: p_{r} > p^{2} (defined in Gordon paper)");
kraus's avatar
kraus committed
363

364 365 366 367 368
        // Gordon, formula (5c)
        ptheta = std::sqrt(p2 - pr2);
        invptheta = 1.0 / ptheta;

        // intepolate values of magnetic field
adelmann's avatar
adelmann committed
369
        MagneticField::interpolate(&bint,&brint,&btint,theta * 180 / M_PI,nradial_m,ntheta_m,y[0],rmin_m,dr_m,bmag_m);
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394
        bint *= invbcon;
        brint *= invbcon;

        // Gordon, formula (5a)
        dydt[0] = y[0] * y[1] * invptheta;
        // Gordon, formula (5b)
        dydt[1] = ptheta - y[0] * bint;
        // Gordon, formulas (9a) and (9b)
        for (size_type i = 2; i < 5; i += 2) {
            dydt[i] = (y[1] * y[i] + y[0] * p2 * y[i+1] * invptheta * invptheta) * invptheta;
            dydt[i+1] = - y[1] * y[i+1] * invptheta - (bint + y[0] * brint) * y[i];
        }
    };

    // define initial state container for integration: y = {r, pr, x1, px1, x2, px2}
    state_type y(6);

    container_type err(2);                                      // difference of last and first value of r (1. element) and pr (2. element)
    container_type delta = {0.0, 0.0};                          // correction term for initial values: r = r + dr, pr = pr + dpr; Gordon, formula (17)
    value_type error = std::numeric_limits<value_type>::max();  // amplitude of error; Gordon, formula (18) (a = a')
    size_type niterations = 0;                                  // if niterations > maxit --> stop iteration

    /*
     * Christian:
     * N = 1440 ---> N = 720 ---> dtheta = 2PI/720 --> nsteps = 721
kraus's avatar
kraus committed
395
     *
396
     * 0, 2, 4, ... ---> jeden zweiten berechnene: 1, 3, 5, ... interpolieren --> 1440 Werte
kraus's avatar
kraus committed
397
     *
398 399
     * Matthias:
     * N = 1440 --> dtheta = 2PI/1440 --> nsteps = 1441
kraus's avatar
kraus committed
400
     *
401
     * 0, 1, 2, 3, 4, 5, ... --> 1440 Werte
kraus's avatar
kraus committed
402
     *
403
     */
kraus's avatar
kraus committed
404

405 406
    // iterate until suggested energy (start with minimum energy)
    value_type E = Emin_m;
kraus's avatar
kraus committed
407 408

    // step size of energy
adelmann's avatar
adelmann committed
409
    value_type dE = (E_m - Emin_m) / (Emax_m - Emin_m);
kraus's avatar
kraus committed
410

adelmann's avatar
adelmann committed
411 412
    // energy increase not more than 0.25
    dE = (dE > 0.25) ? 0.25 : dE;
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428

    // energy dependent values
    value_type en = E / physics::E0;                      // en = E/E0 = E/(mc^2) (E0 is kinetic energy)
    value_type p = acon * std::sqrt(en * (2.0 + en));     // momentum [p] = m; Gordon, formula (3)
    value_type gamma2 = (1.0 + en) * (1.0 + en);          // = gamma^2
    value_type beta = std::sqrt(1.0 - 1.0 / gamma2);
    p2 = p * p;                                           // p^2 = p*p
    value_type invgamma4 = 1.0 / (gamma2 * gamma2);       // = 1/gamma^4

    // set initial values for radius and radial momentum for lowest energy Emin
    // orbit, [r] = m; Gordon, formula (20)
    // radial momentum; Gordon, formula (20)
    container_type init = {beta * acon, 0.0};

    // store initial values for updating values for higher energies
    container_type previous_init = {0.0, 0.0};
kraus's avatar
kraus committed
429

adelmann's avatar
adelmann committed
430
    do {
431 432

        // (re-)set inital values for r and pr
kraus's avatar
kraus committed
433
        r_m[0] = init[0];
434
        pr_m[0] = init[1];
kraus's avatar
kraus committed
435

436 437 438 439 440 441 442 443 444 445 446 447
        // integrate until error smaller than user-define accuracy
        do {
            // (re-)set inital values
            x_m[0]  = 1.0;               // x1; Gordon, formula (10)
            px_m[0] = 0.0;               // px1; Gordon, formula (10)
            x_m[1]  = 0.0;               // x2; Gordon, formula (10)
            px_m[1] = 1.0;               // px2; Gordon, formula (10)
            nxcross_m = 0;               // counts the number of crossings of x-axis (excluding first step)
            idx = 0;                     // index for looping over r and pr arrays

            // fill container with initial states
            y = {init[0],init[1], x_m[0], px_m[0], x_m[1], px_m[1]};
kraus's avatar
kraus committed
448

449 450
            // integrate from 0 to 2*pi (one has to get back to the "origin")
            boost::numeric::odeint::integrate_n_steps(stepper_m,orbit_integration,y,0.0,dtheta_m,N_m,store);
kraus's avatar
kraus committed
451

452 453 454 455 456
            // write new state
            x_m[0] = y[2];
            px_m[0] = y[3];
            x_m[1] = y[4];
            px_m[1] = y[5];
kraus's avatar
kraus committed
457

458 459 460 461
            // 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)
            err[0] = r_m[N_m] - r_m[0];      // Gordon, formula (14)
            err[1] = pr_m[N_m] - pr_m[0];    // Gordon, formula (14)
kraus's avatar
kraus committed
462

463 464 465 466
            // correct inital values of r and pr
            invdenom = 1.0 / (x_m[0] + px_m[1] - 2.0);
            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)
kraus's avatar
kraus committed
467

468 469 470
            // improved initial values; Gordon, formula (17) (here it's used for higher energies)
            init[0] += delta[0];
            init[1] += delta[1];
kraus's avatar
kraus committed
471

472 473
            // compute amplitude of the error
            error = std::sqrt(delta[0] * delta[0] + delta[1] * delta[1] * invgamma4) / r_m[0];
kraus's avatar
kraus committed
474

475
        } while (error > accuracy && niterations++ < maxit);
kraus's avatar
kraus committed
476

477 478
        // reset iteration counter
        niterations = 0;
kraus's avatar
kraus committed
479

480 481
        // reset correction term
        delta[0] = delta[1] = 0.0;
adelmann's avatar
adelmann committed
482 483 484 485 486 487

        // increase energy by dE
        if (E_m <= E + dE)
            E = E_m;
        else
            E += dE;
kraus's avatar
kraus committed
488

489 490 491 492 493 494
        // set constants for new energy E
        en = E / physics::E0;                     // en = E/E0 = E/(mc^2) (E0 is kinetic energy)
        p = acon * std::sqrt(en * (2.0 + en));    // momentum [p] = m; Gordon, formula (3)
        p2 = p * p;                               // p^2 = p*p
        gamma2 = (1.0 + en) * (1.0 + en);
        invgamma4 = 1.0 / (gamma2 * gamma2);
kraus's avatar
kraus committed
495 496


adelmann's avatar
adelmann committed
497
    } while (E != E_m);
kraus's avatar
kraus committed
498

499 500 501 502 503
    /* store last entry, since it is needed in computeVerticalOscillations(), because we have to do the same
     * number of integrations steps there.
     */
    lastOrbitVal_m = r_m[N_m];           // needed in computeVerticalOscillations()
    lastMomentumVal_m = pr_m[N_m];       // needed in computeVerticalOscillations()
kraus's avatar
kraus committed
504

505 506 507
    // 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();
kraus's avatar
kraus committed
508

509

510 511 512
    // returns true if converged, otherwise false
    return error < accuracy;
}
513 514 515

template<typename Value_type, typename Size_type, class Stepper>
Value_type ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeTune(const std::array<value_type,2>& y, value_type py2, size_type ncross) {
516
    // Y = [y1, y2; py1, py2]
kraus's avatar
kraus committed
517

518 519
    // cos(mu)
    value_type cos = 0.5 * (y[0] + py2);
kraus's avatar
kraus committed
520

521 522
    value_type twopi = 2.0 * M_PI;
    value_type mu;
kraus's avatar
kraus committed
523

524 525
    // sign of sin(mu) has to be equal to y2
    bool negative = std::signbit(y[1]);
kraus's avatar
kraus committed
526

527
    bool uneven = (ncross % 2);
kraus's avatar
kraus committed
528

529 530 531
    if (std::fabs(cos) > 1.0) {
        // store the number of crossings
        value_type n = ncross;
kraus's avatar
kraus committed
532

533 534
        if (uneven)
            n = ncross - 1;
kraus's avatar
kraus committed
535

536 537 538
        // Gordon, formula (36b)
        value_type muPrime = -std::acosh(std::fabs(cos));      // mu'
        mu = n * M_PI + muPrime;
kraus's avatar
kraus committed
539

540 541 542 543 544 545 546
    } else {
        value_type muPrime = (uneven) ? std::acos(-cos) : std::acos(cos);    // mu'
        /* It has to be fulfilled: 0<= mu' <= pi
        * But since |cos(mu)| <= 1, we have
        * -1 <= cos(mu) <= 1 --> 0 <= mu <= pi (using above programmed line), such
        * that condition is already fulfilled.
        */
kraus's avatar
kraus committed
547

548 549
        // Gordon, formula (36)
        mu = ncross * M_PI + muPrime;
kraus's avatar
kraus committed
550

551 552 553 554 555 556 557
        // if sign(y[1]) > 0 && sign(sin(mu)) < 0
        if (!negative && std::signbit(std::sin(mu))) {
            mu = ncross * M_PI - muPrime;
        } else if (negative && !std::signbit(std::sin(mu))) {
            mu = ncross * M_PI - muPrime + twopi;
        }
    }
kraus's avatar
kraus committed
558

559
    // nu = mu/theta, where theta = integration domain
kraus's avatar
kraus committed
560

adelmann's avatar
adelmann committed
561 562 563 564 565
    /* domain_m = true --> only integrated over a single sector --> multiply by nSector_m to
     * get correct tune.
     */
    if (domain_m)
        mu *= nSector_m;
kraus's avatar
kraus committed
566

567
    return mu / twopi;
568 569 570
}

template<typename Value_type, typename Size_type, class Stepper>
571
void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeOrbitProperties() {
kraus's avatar
kraus committed
572
    /*
573 574 575 576 577
     * 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
     */
kraus's avatar
kraus committed
578

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

582 583 584 585 586 587
    value_type invbcon = 1.0 / physics::bcon(wo_m);
    value_type en = E_m / physics::E0;                                  // en = E/E0 = E/(mc^2) (E0 is kinetic energy)
    value_type p = physics::acon(wo_m) * std::sqrt(en * (2.0 + en));    // momentum [p] = m; Gordon, formula (3)
    value_type p2 = p * p;
    value_type theta = 0.0;                                             // angle for interpolating
    value_type ptheta;
588

589 590 591 592 593 594 595
    // 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
adelmann's avatar
adelmann committed
596
        MagneticField::interpolate(&bint,&brint,&btint,theta * 180.0 / M_PI,nradial_m,ntheta_m,r_m[i],rmin_m,dr_m,bmag_m);
597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
        bint *= invbcon;
        brint *= invbcon;
        btint *= invbcon;

        // inverse bending radius
        h_m[i] = bint / p;

        // local field index
        ptheta = std::sqrt(p2 - pr_m[i] * pr_m[i]);
        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;
613
    }
614 615 616

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

template<typename Value_type, typename Size_type, class Stepper>
620
void ClosedOrbitFinder<Value_type, Size_type, Stepper>::computeVerticalOscillations() {
kraus's avatar
kraus committed
621

622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641
    vertOscDone_m = true;

    // READ IN MAGNETIC FIELD: ONLY FOR STAND-ALONE PROGRAM
    value_type bint, brint, btint; // B, dB/dr, dB/dtheta

    value_type en = E_m / physics::E0;                                  // en = E/E0 = E/(mc^2) with kinetic energy E0
    value_type p = physics::acon(wo_m) * std::sqrt(en *(en + 2.0));     // Gordon, formula (3)
    value_type p2 = p * p;                                              // p^2 = p*p
    size_type idx = 0;                                                  // index for going through container
    value_type pr2;                                                     // pr^2 = pr*pr
    value_type ptheta, invptheta;                                       // Gordon, formula (5c)
    value_type zold = 0.0;                                              // for counting nzcross

    // store bcon locally
    value_type invbcon = 1.0 / physics::bcon(wo_m);                     // [bcon] = MeV*s/(C*m^2) = 10^6 T = 10^7 kG (kilo Gauss)

    // define the ODEs (using lambda function)
    std::function<void(const state_type&, state_type&, const double)> vertical = [&](const state_type &y, state_type &dydt, const double theta){
        pr2 = y[1] * y[1];
        if (p2 < pr2)
adelmann's avatar
adelmann committed
642
            throw std::domain_error("Error in ClosedOrbitFinder::computeVerticalOscillations: p_{r} > p^{2} (defined in Gordon paper)");
kraus's avatar
kraus committed
643

644 645 646 647 648
        // Gordon, formula (5c)
        ptheta = std::sqrt(p2 - pr2);
        invptheta = 1.0 / ptheta;

        // intepolate values of magnetic field
adelmann's avatar
adelmann committed
649
        MagneticField::interpolate(&bint,&brint,&btint,theta * 180 / M_PI,nradial_m,ntheta_m,y[0],rmin_m,dr_m,bmag_m);
650 651 652
        bint *= invbcon;
        brint *= invbcon;
        btint *= invbcon;
kraus's avatar
kraus committed
653

654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
        // We have to integrate r and pr again, otherwise we don't have the Runge-Kutta of the B-field
        // Gordon, formula (5a)
        dydt[0] = y[0] * y[1] * invptheta;
        // Gordon, formula (5b)
        dydt[1] = ptheta - y[0] * bint;

        // Gordon, formulas (22a) and (22b)
        for (size_type i = 2; i < 5; i += 2) {
            dydt[i] = y[0] * y[i+1] * invptheta;
            dydt[i+1] = (y[0] * brint - y[1] * invptheta * btint) * y[i];
        }

        // integrate phase
        dydt[6] = y[0] * invptheta * gamma_m - 1;
    };

    // to get next index for r and pr (to iterate over container)
    auto next = [&](state_type& y, const value_type t) {
        // number of times z2 changes sign
        nzcross_m += (idx > 0) * (y[4] * zold < 0);
        zold = y[4];
        ++idx;
    };

    // set initial state container for integration: y = {r, pr, z1, pz1, z2, pz2, phase}
    state_type y = {r_m[0], pr_m[0], 1.0, 0.0, 0.0, 1.0, 0.0};

    // add last element for integration (since we have to return to the initial point (--> size = N_m+1, capacity = N_m+1)
    r_m.push_back(lastOrbitVal_m);
    pr_m.push_back(lastMomentumVal_m);
kraus's avatar
kraus committed
684

685 686
    // integrate: assume no imperfections --> only integrate over a single sector (dtheta_m = 2pi/N_m)
    boost::numeric::odeint::integrate_n_steps(stepper_m,vertical,y,0.0,dtheta_m,N_m,next);
kraus's avatar
kraus committed
687

688 689 690
    // remove last element again (--> size = N_m, capacity = N_m+1)
    r_m.pop_back();
    pr_m.pop_back();
kraus's avatar
kraus committed
691

692 693 694 695 696 697
    // write new state
    z_m[0] = y[2];
    pz_m[0] = y[3];
    z_m[1] = y[4];
    pz_m[1] = y[5];
    phase_m = y[6] / (2.0 * M_PI);
kraus's avatar
kraus committed
698

adelmann's avatar
adelmann committed
699 700 701 702 703
    /* domain_m = true --> only integrated over a single sector
     * --> multiply by nSector_m to get correct phase_m
     */
    if (domain_m)
        phase_m *= nSector_m;
704 705
}

kraus's avatar
kraus committed
706
#endif