Hamiltonian.cpp 5.06 KB
Newer Older
1 2 3 4 5 6 7 8 9
//
// Source file of the Hamiltonian class,
//   Constructs thick lens Hamiltonians up to arbitrary order.
//
// Copyright (c) 2018, Philippe Ganz, ETH Zürich
// All rights reserved
//
// OPAL is licensed under GNU GPL version 3.

snuverink_j's avatar
snuverink_j committed
10 11 12 13
#include "Hamiltonian.h"

Hamiltonian::Hamiltonian(int truncOrder) : truncOrder_m(truncOrder)
{
14

snuverink_j's avatar
snuverink_j committed
15
    series_t::setGlobalTruncOrder(truncOrder+1);
16

snuverink_j's avatar
snuverink_j committed
17 18 19 20 21 22 23 24 25 26 27
    x     = series_t::makeVariable(0);
    px    = series_t::makeVariable(1);
    y     = series_t::makeVariable(2);
    py    = series_t::makeVariable(3);
    z     = series_t::makeVariable(4);
    delta = series_t::makeVariable(5);
}

Hamiltonian::series_t Hamiltonian::drift(const double& gamma0)
{
    double beta0 = this->getBeta_m(gamma0);
28

snuverink_j's avatar
snuverink_j committed
29 30 31 32 33 34 35
    return ( delta / beta0 )
            - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
                    - ( px*px )
                    - ( py*py )
                    - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1);
}

36

snuverink_j's avatar
snuverink_j committed
37 38
Hamiltonian::series_t Hamiltonian::rbend(double& beta0,
                            double& gamma0,
39
                            double& /*q*/,
snuverink_j's avatar
snuverink_j committed
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61
                            double& h,
                            double& k0)
{



        return ( delta / beta0 )
        - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
                    - ( px*px )
                    - ( py*py )
                    - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m+1
            ))
    - (h * x)
    * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
                    - ( px*px )
                    - ( py*py )
                    - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m
            ))
    + k0 * x * (1. + 0.5 * h* x);


}
62 63


snuverink_j's avatar
snuverink_j committed
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
Hamiltonian::series_t Hamiltonian::sbend(const double& gamma0,
                                         const double& h,
                                         const double& K0)
{
    double beta0 = this->getBeta_m(gamma0);

    return ( delta / beta0 )
                    - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
                                    - ( px*px )
                                    - ( py*py )
                                    - 1./( beta0*beta0 * gamma0*gamma0 ),(truncOrder_m+1)
                            ))
                    - (h * x)
                    * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
                                    - ( px*px )
                                    - ( py*py )
                                    - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m
                            ))
                    + K0 * x * (1. + 0.5 * h* x);
}


Hamiltonian::series_t Hamiltonian::bendFringe(
                double& beta0,
                double& gamma0,
89 90
                double& /*h*/,
                double& /*k0*/,
snuverink_j's avatar
snuverink_j committed
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
                series_t& ax,
                series_t& az)
{
    if (truncOrder_m == 2){
        return ( delta / beta0 )
                        - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
                                - ( px*px)
                                - ( py*py )
                                - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1
                        ) - az;
    }else{
        return ( delta / beta0 )
            - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
                    - ( px*px - 2.0*px*ax - ax*ax)
                    - ( py*py )
                    - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1
            ) - az;
    }
    //std::cout << H << std::endl;
    //std::cout << H.getMaxOrder() << std::endl;
    /*H=( delta / beta0 )
                        - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
                                        - ( px - ax )*( px - ax )
                                        - ( py*py )
                                        - 1./( beta0*beta0 * gamma0*gamma0 ),(truncOrder_m+1)))
                        - (h * x)
                        * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
                                        - ( px - ax )*( px - ax )
                                        - ( py*py )
                                        - 1./( beta0*beta0 * gamma0*gamma0 ), truncOrder_m))
                         + h * x * az;*/
}


Hamiltonian::series_t Hamiltonian::quadrupole(const double& gamma0,
126
                                              const double& /*q*/,
snuverink_j's avatar
snuverink_j committed
127 128 129
                                              const double& k1)
{
    double beta0 = this->getBeta_m(gamma0);
130

snuverink_j's avatar
snuverink_j committed
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
    return ( delta / beta0 )
    - sqrt ((1./ beta0 + delta ) *(1./ beta0 + delta)
            - ( px*px )
            - ( py*py )
            - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m+1
    )
    + 0.5 * k1 * (x*x - y*y);
}


Hamiltonian::series_t Hamiltonian::fringeField(const double& phi,
                                               const double& k0)
{
    if ( truncOrder_m > 1 ) {
        //TODO higher order thin lens approximation
    }
    // else
148

snuverink_j's avatar
snuverink_j committed
149 150 151 152 153 154 155 156 157 158
    double angle = phi;
    if ( k0 < 0.0 )
        angle = -phi;
    return -0.5 * (x * x - y * y) * k0 * std::tan(angle);
}


double Hamiltonian::getBeta_m(const double& gamma) {
    return std::sqrt(1.0 - 1.0 / (gamma * gamma) );
}