OPAL (Object Oriented Parallel Accelerator Library)  2021.1.99
OPAL
Hamiltonian.cpp
Go to the documentation of this file.
1 //
2 // Class: Hamiltonian
3 // Constructs thick lens Hamiltonian up to arbitrary order for beamline elements
4 //
5 // Copyright (c) 2018, Philippe Ganz, ETH Zürich
6 // All rights reserved
7 //
8 // Implemented as part of the Master thesis
9 // "s-based maps from TPS & Lie-Series applied to Proton-Therapy Gantries"
10 //
11 // This file is part of OPAL.
12 //
13 // OPAL is free software: you can redistribute it and/or modify
14 // it under the terms of the GNU General Public License as published by
15 // the Free Software Foundation, either version 3 of the License, or
16 // (at your option) any later version.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
20 //
21 
22 #include "Hamiltonian.h"
23 
24 Hamiltonian::Hamiltonian(int truncOrder) : truncOrder_m(truncOrder)
25 {
26 
27  series_t::setGlobalTruncOrder(truncOrder+1);
28 
35 }
36 
38 {
39  double beta0 = this->getBeta_m(gamma0);
40 
41  return ( delta / beta0 )
42  - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
43  - ( px*px )
44  - ( py*py )
45  - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1);
46 }
47 
48 
50  double& gamma0,
51  double& /*q*/,
52  double& h,
53  double& k0)
54 {
55 
56 
57 
58  return ( delta / beta0 )
59  - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
60  - ( px*px )
61  - ( py*py )
62  - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m+1
63  ))
64  - (h * x)
65  * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
66  - ( px*px )
67  - ( py*py )
68  - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m
69  ))
70  + k0 * x * (1. + 0.5 * h* x);
71 
72 
73 }
74 
75 
77  const double& h,
78  const double& K0)
79 {
80  double beta0 = this->getBeta_m(gamma0);
81 
82  return ( delta / beta0 )
83  - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
84  - ( px*px )
85  - ( py*py )
86  - 1./( beta0*beta0 * gamma0*gamma0 ),(truncOrder_m+1)
87  ))
88  - (h * x)
89  * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
90  - ( px*px )
91  - ( py*py )
92  - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m
93  ))
94  + K0 * x * (1. + 0.5 * h* x);
95 }
96 
97 
99  double& beta0,
100  double& gamma0,
101  double& /*h*/,
102  double& /*k0*/,
103  series_t& ax,
104  series_t& az)
105 {
106  if (truncOrder_m == 2){
107  return ( delta / beta0 )
108  - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
109  - ( px*px)
110  - ( py*py )
111  - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1
112  ) - az;
113  }else{
114  return ( delta / beta0 )
115  - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
116  - ( px*px - 2.0*px*ax - ax*ax)
117  - ( py*py )
118  - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1
119  ) - az;
120  }
121  //std::cout << H << std::endl;
122  //std::cout << H.getMaxOrder() << std::endl;
123  /*H=( delta / beta0 )
124  - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
125  - ( px - ax )*( px - ax )
126  - ( py*py )
127  - 1./( beta0*beta0 * gamma0*gamma0 ),(truncOrder_m+1)))
128  - (h * x)
129  * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
130  - ( px - ax )*( px - ax )
131  - ( py*py )
132  - 1./( beta0*beta0 * gamma0*gamma0 ), truncOrder_m))
133  + h * x * az;*/
134 }
135 
136 
138  const double& /*q*/,
139  const double& k1)
140 {
141  double beta0 = this->getBeta_m(gamma0);
142 
143  return ( delta / beta0 )
144  - sqrt ((1./ beta0 + delta ) *(1./ beta0 + delta)
145  - ( px*px )
146  - ( py*py )
147  - 1./( beta0*beta0 * gamma0*gamma0 ),truncOrder_m+1
148  )
149  + 0.5 * k1 * (x*x - y*y);
150 }
151 
152 
154  const double& k0)
155 {
156  if ( truncOrder_m > 1 ) {
157  //TODO higher order thin lens approximation
158  }
159  // else
160 
161  double angle = phi;
162  if ( k0 < 0.0 )
163  angle = -phi;
164  return -0.5 * (x * x - y * y) * k0 * std::tan(angle);
165 }
166 
167 
168 double Hamiltonian::getBeta_m(const double& gamma) {
169  return std::sqrt(1.0 - 1.0 / (gamma * gamma) );
170 }
Tps< T > tan(const Tps< T > &x)
Tangent.
Definition: TpsMath.h:147
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
series_t delta
Definition: Hamiltonian.h:45
Hamiltonian::series_t sbend(const double &gamma0, const double &h, const double &k0)
Definition: Hamiltonian.cpp:76
Hamiltonian::series_t drift(const double &gamma0)
Definition: Hamiltonian.cpp:37
Hamiltonian::series_t rbend(double &beta0, double &gamma0, double &q, double &h, double &K0)
Definition: Hamiltonian.cpp:49
series_t z
Definition: Hamiltonian.h:44
Hamiltonian::series_t quadrupole(const double &gamma0, const double &q, const double &k1)
Hamiltonian(int truncOrder)
Definition: Hamiltonian.cpp:24
series_t y
Definition: Hamiltonian.h:42
series_t py
Definition: Hamiltonian.h:43
series_t x
Definition: Hamiltonian.h:40
Hamiltonian::series_t bendFringe(double &beta0, double &gamma0, double &h, double &k0, series_t &ax, series_t &az)
Definition: Hamiltonian.cpp:98
Hamiltonian::series_t fringeField(const double &phi, const double &k0)
series_t px
Definition: Hamiltonian.h:41
static FTps makeVariable(int var)
Make variable.
Definition: FTps.hpp:481
static void setGlobalTruncOrder(int order)
Set the global truncation order.
Definition: FTps.hpp:419