Hamiltonian.cpp

Go to the documentation of this file.

//
// Class: Hamiltonian
//   Constructs thick lens Hamiltonian up to arbitrary order for beamline elements
//
// Copyright (c) 2018, Philippe Ganz, ETH Zürich
// All rights reserved
//
// Implemented as part of the Master thesis
// "s-based maps from TPS & Lie-Series applied to Proton-Therapy Gantries"
//
// 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/>.
//
 
#include "Hamiltonian.h"
 
Hamiltonian::Hamiltonian(int truncOrder) : truncOrder_m(truncOrder)
{
 
    series_t::setGlobalTruncOrder(truncOrder+1);
 
    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);
 
    return ( delta / beta0 )
            - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
                    - ( px*px )
                    - ( py*py )
                    - 1./( beta0 * beta0 * gamma0 * gamma0 ),truncOrder_m+1);
}
 
 
Hamiltonian::series_t Hamiltonian::rbend(double& beta0,
                            double& gamma0,
                            double& /*q*/,
                            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);
 
 
}
 
 
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,
                double& /*h*/,
                double& /*k0*/,
                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,
                                              const double& /*q*/,
                                              const double& k1)
{
    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
    )
    + 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
 
    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) );
}