From f123c55d87d0bf272a0cf29630b76d211369fd62 Mon Sep 17 00:00:00 2001
From: Philippe Ganz <gp@student.ethz.ch>
Date: Tue, 30 Jan 2018 13:41:36 +0100
Subject: [PATCH] Implement map generation in Ticktracker::execute
---
Doxyfile | 8 +-
doc/randomnotes.txt | 13 -
doc/refcard/opal-ref-card-bibliography.bib | 39 --
doc/refcard/opal-refcard.tex | 756 ---------------------
src/Algorithms/ThickTracker.cpp | 584 +++++++++++++---
src/Classic/Algorithms/PartData.h | 2 +-
src/Elements/OpalBeamline.h | 2 +-
tests/Maps/map.in | 27 +-
8 files changed, 513 insertions(+), 918 deletions(-)
delete mode 100644 doc/randomnotes.txt
delete mode 100644 doc/refcard/opal-ref-card-bibliography.bib
delete mode 100644 doc/refcard/opal-refcard.tex
diff --git a/Doxyfile b/Doxyfile
index 4770de8c4..73556de8e 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -51,7 +51,7 @@ PROJECT_BRIEF = "OPAL"
# and the maximum width should not exceed 200 pixels. Doxygen will copy the logo
# to the output directory.
-PROJECT_LOGO =
+PROJECT_LOGO = ./amas.png
# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) path
# into which the generated documentation will be written. If a relative path is
@@ -1572,7 +1572,7 @@ EXTRA_SEARCH_MAPPINGS =
# If the GENERATE_LATEX tag is set to YES doxygen will generate LaTeX output.
# The default value is: YES.
-GENERATE_LATEX = NO
+GENERATE_LATEX = YES
# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. If a
# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
@@ -1675,7 +1675,7 @@ PDF_HYPERLINKS = YES
# The default value is: YES.
# This tag requires that the tag GENERATE_LATEX is set to YES.
-USE_PDFLATEX = NO
+USE_PDFLATEX = YES
# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \batchmode
# command to the generated LaTeX files. This will instruct LaTeX to keep running
@@ -2328,4 +2328,4 @@ GENERATE_LEGEND = YES
# The default value is: YES.
# This tag requires that the tag HAVE_DOT is set to YES.
-DOT_CLEANUP = YES
\ No newline at end of file
+DOT_CLEANUP = YES
diff --git a/doc/randomnotes.txt b/doc/randomnotes.txt
deleted file mode 100644
index 302de8f68..000000000
--- a/doc/randomnotes.txt
+++ /dev/null
@@ -1,13 +0,0 @@
-
-A note on documentation etc. for collegues starting work on OPAL:
-
-1) The full documentation gor OPAL 2.x.x is on the Wiki: https://gitlab.psi.ch/OPAL/src/wikis/Manual/home
-
-2) Also the C++ source code must be documented please read: http://en.wikipedia.org/wiki/Doxygen and document
- your source code accordingly.
-
-3) We have to do some artwork for our paper(s) and documentation. Please use only TikZ and see the power of
- this package (and many examples to start) at http://www.texample.net/tikz/examples/ .
-
-
-
diff --git a/doc/refcard/opal-ref-card-bibliography.bib b/doc/refcard/opal-ref-card-bibliography.bib
deleted file mode 100644
index 07f7e80bf..000000000
--- a/doc/refcard/opal-ref-card-bibliography.bib
+++ /dev/null
@@ -1,39 +0,0 @@
-%% objective - bibliography file for femaxx reference card
-%% author - benedikt oswald
-%% modified - 2008 jun 11, creation, benedikt oswald
-
-@Article{arbenzetal2001,
- author = {Arbenz, P. and Geus, R. and Adam, S.},
- title = {Solving Maxwell eigenvalue problems for accelerating cavities},
- journal = {Physical Review Special Topics - Accelerators and Beams},
- year = {2001},
- volume = {4},
- pages = {022001-1:022001-10},
- annote = {Abstract: we investigate algorithms for computing steady state electromagnetic waves
- in cavities. The Maxwell equations for the strength of the electric field are solved by
- a mixed method with quadratic finite edge (N{\'e}d{\'e}lec) elements for the field
- values and corresponding node-based finite elements for the Lagrange multiplier, This
- approach avoids so-called spurious modes which are introduced if the divergence-free
- condition for the electric field is not treated properly. To compute a few of the smalles
- positive eigenvalues and corresponding eigenmodes of the resulting large sparse-matrix
- eigenvalue problems, two algorithms have been used: the implicitly restarted Lanczos
- algorithm and the Jacobi-Davidson algorithm, both with shit-and-invert spectral transformation.
- Two-level hierarchical basis preconditioners have been employed for the iterative solution
- of the resulting systems of equations.}
-}
-
-
-@PhdThesis{geus2002,
- author = {Geus, R.},
- title = {The Jacobi-Davidson Algorithm for solving large sparse symmetric eigenvalue problems with application to the design of accelerator cavities},
- school = {Swiss Federal Institute of Technology Zurich},
- year = {2002},
- OPTkey = {Diss. ETH NO. 14734},
- type = {},
- address = {},
- month = {},
- note = {},
- annote = {}
-}
-
-
diff --git a/doc/refcard/opal-refcard.tex b/doc/refcard/opal-refcard.tex
deleted file mode 100644
index 78facc009..000000000
--- a/doc/refcard/opal-refcard.tex
+++ /dev/null
@@ -1,756 +0,0 @@
-%% objective - reference card for the \opal eigenvalue computation code
-%% objective - list and explain \opal code usage
-%% modified - 2008 jun 10, creation from latex cheat sheet file, benedikt oswald
-%%
-%%
-
-\documentclass[10pt,landscape,a4paper]{article}
-\usepackage{color}
-\usepackage{multicol}
-\usepackage{calc}
-\usepackage{ifthen}
-
-\usepackage{amsmath} %% diverse Matheerweiterungen
-\usepackage{amssymb} %% diverse Matheerweiterungen, z.B. \mathbb{R}
-\usepackage{wasysym} %% diverse Matheerweiterungen, z.B. \oiint
-\usepackage{epsfig} %% um eps-Dateien einzubinden (\epsfig{file=...})
-\usepackage{longtable} %% fuer Tabellen ueber mehrere Seiten
-
-\usepackage{graphicx}
-\usepackage{makeidx}
-
-\usepackage{natbib}
-
-\usepackage[landscape]{geometry}
-% If you're reading this, be prepared for confusion. Making this was
-% a learning experience for me, and it shows. Much of the placement
-% was hacked in; if you make it better, let me know...
-
-
-% 2008-04
-% Changed page margin code to use the geometry package. Also added code for
-% conditional page margins, depending on paper size. Thanks to Uwe Ziegenhagen
-% for the suggestions.
-
-% 2006-08
-% Made changes based on suggestions from Gene Cooperman. <gene at ccs.neu.edu>
-
-
-% To Do:
-% \listoffigures \listoftables
-% \setcounter{secnumdepth}{0}
-
-
-% This sets page margins to .5 inch if using letter paper, and to 1cm
-% if using A4 paper. (This probably isn't strictly necessary.)
-% If using another size paper, use default 1cm margins.
-\ifthenelse{\lengthtest { \paperwidth = 11in}}
- { \geometry{top=.5in,left=.5in,right=.5in,bottom=.5in} }
- {\ifthenelse{ \lengthtest{ \paperwidth = 297mm}}
- {\geometry{top=1cm,left=1cm,right=1cm,bottom=1cm} }
- {\geometry{top=1cm,left=1cm,right=1cm,bottom=1cm} }
- }
-
-% Turn off header and footer
-\pagestyle{empty}
-\bibliographystyle{plain}
-
-% Redefine section commands to use less space
-\makeatletter
-\renewcommand{\section}{\@startsection{section}{1}{0mm}%
- {-1ex plus -.5ex minus -.2ex}%
- {0.5ex plus .2ex}%x
- {\normalfont\large\bfseries}}
-\renewcommand{\subsection}{\@startsection{subsection}{2}{0mm}%
- {-1explus -.5ex minus -.2ex}%
- {0.5ex plus .2ex}%
- {\normalfont\normalsize\bfseries}}
-\renewcommand{\subsubsection}{\@startsection{subsubsection}{3}{0mm}%
- {-1ex plus -.5ex minus -.2ex}%
- {1ex plus .2ex}%
- {\normalfont\small\bfseries}}
-
-
-\newcommand{\opal}{\textsc{OPAL }}
-\newcommand{\opalt}{\textsc{OPAL-t }}
-\newcommand{\opalcycl}{\textsc{OPAL-cycl }}
-\newcommand{\opalmap}{\textsc{OPAL-map }}
-
-\newcommand{\noopalt}{\textsc{\leftthumbsdown \opalt }}
-\newcommand{\noopalmap}{\textsc{\leftthumbsdown \opal }}
-\newcommand{\noopalcycl}{\textsc{\leftthumbsdown \opalcycl }}
-
-
-
-
-\makeatother
-
-% Define BibTeX command
-\def\BibTeX{{\rm B\kern-.05em{\sc i\kern-.025em b}\kern-.08em
- T\kern-.1667em\lower.7ex\hbox{E}\kern-.125emX}}
-
-% Don't print section numbers
-\setcounter{secnumdepth}{0}
-
-
-\setlength{\parindent}{0pt}
-\setlength{\parskip}{0pt plus 0.5ex}
-
-
-%% -----------------------------------------------------------------------
-%% begin \opal refcard
-%% -----------------------------------------------------------------------
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% New commands.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\newcommand{\unit}[1]{\mbox{$\mathrm{#1}$}} %% physical unit
-\renewcommand{\vec}[1]{\mathbf{#1}} %% vector
-\newcommand{\uvec}[1]{\mathbf{\hat{#1}}} %% unit vector
-\newcommand{\mat}[1]{\mathsf{#1}} %% matrix
-\newcommand{\func}[1]{\mbox{$\mathrm{#1}$}} %% math. function
-\newcommand{\cu}{i} %% complex unit
-%% replace this with an argument of resp. type.
-\newcommand{\replace}[1]{\texttt{\mbox{<\textsl{#1}>}}}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Colors.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\definecolor{CmdBoxBackground}{gray}{0.7}
-\definecolor{ExampleBoxBackground}{gray}{0.9}
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\begin{document}
-
-% \raggedright
-\footnotesize
-\begin{multicols}{3}
-
-
-% multicol parameters
-% These lengths are set only within the two main columns
-%\setlength{\columnseprule}{0.25pt}
-\setlength{\premulticols}{1pt}
-\setlength{\postmulticols}{1pt}
-\setlength{\multicolsep}{1pt}
-\setlength{\columnsep}{2pt}
-
-
-
-
-
-\begin{center}
- \LARGE{\large{\textbf{\textsf{\opal}}} reference card} \\
-\end{center}
-
-
-\begin{center}
- \includegraphics[width=70mm]{./pictures/torus200kmode0electricfield}
-% We have computed a few modes of a toroidal resonator structure, with
-% respective diameters of $1.0$ and $0.5$ meters. The plot shows the magnitude
-% of the electric field of the lowest mode with a resonance frequency
-% of $218.067+08$ Mhz and the associated quality factor $Q = 62'870e$ [-].
-% Mesh of the torus courtesy Markus Bopp, PSI.
-\end{center}
-
-
-
-
-
-
-
-
-
-\textsf{\textbf{Applicability:}} Commands in this reference card
-apply to the \textsf{\opal} eigenvalue code only. The electromagnetic
-frontend \textsf{heronion} is documented elsewhere.
-%%
-%%
-%%
-\textsf{\textbf{Iron Rule:}} parameters that are used for the solution of
-a specific electromagnetic problem are communicated to the code
-at one single location only, i.e. as command line parameters when
-the solver is run.
-%%
-%%
-\textsf{\textbf{Units:}} the \textsf{\opal} code uses the MKS system,
-aka. meter, kilogram, second, for all its computations.
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{\textsf{Usage Paradigm}}
-
-
-\texttt{\opal} is an advanced computational electrodynamics package consisting of three different,
-subsequently used programs.
-%%
-%%
-This document describes code usage. For detailed information on implementational questions
-studying the code is essential and the documentation generated by Doxygen is helpful.
-%%
-%%
-%%
-
-
-
-
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{\textsf{Physics}}
-
-The \textsf{\opal} code computes the eigenmodal fields of electromagnetic
-resonators by solving the electric field vector wave (\textit{curl-curl})
-equation in $3$-dimensional space, in frequency domain and on a tetrahedral
-mesh for scalar, non-dispersive dielectric and magnetic materials.
-%%
-%%
-\begin{equation}
-\label{eq_electric_curl_curl}
-%%
-\nabla \times \frac{1}{\mu_{r}} \nabla \times \vec{E}
-%%
-+ k_{0}^2 \epsilon_{r} \vec{E}
-%%
-= 0
-%%
-\end{equation}
-%%
-At present materials do not model electromagnetic losses. The computation
-of a resonator's quality figure uses an \textit{a posteriori}
-perturbation approach.
-%%
-%%
-The eigenvalue $\ell = k_{0}^2$ is defined in terms of electromagnetic constants
-and the resonance frequency $f_{\mathrm{res}}$.
-%%
-%%
-The resonance frequency $f_{\mathrm{res}}$ is computed from the eigenvalue $\ell$
-%%
-%%
-\begin{equation}
- \label{eq_fres_from_k0}
- f = \Biggl( \sqrt{\frac{\ell}{\epsilon_{0}\mu_{0}}} \: \: \Biggr) \frac{1}{2\pi}
-\end{equation}
-%%
-%%
-The eigenvalue shift $\sigma$ can also be understood as a frequency shift,
-using Eq. (\ref{eq_fres_from_k0}), cf. also the \texttt{--sigma} parameter.
-%%
-The theory and the implementation of \textsf{\opal} is explained
-in \cite[]{arbenzetal2001,geus2002}.
-
-
-
-
-
-
-
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\section{\textsf{Command Line Parameters}}
-
-
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{General}
-\label{section_general}
-
-The executable of the \textsf{\opal} code is \texttt{\opal\_driver}. Running the code
-depends on the operating system.
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{\opal\_driver}
-}
-
-\vspace{1mm}
-
-
-\noindent \textbf{\textsf{Nota bene:}} On a Linux system, it is common
-practice to use the \texttt{mpirun} command, e.g.
-
-\vspace{1mm}
-
-\noindent \texttt{mpirun -np 2 <path\_to\_\opal\_executable>/\opal\_driver ...}
-
-\vspace{1mm}
-
-The situation is different on Cray or IBM/bluegene systems.
-
-
-
-
-
-
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Help}
-\label{section_help}
-
-To request a description of all built-in parameters of \textsf{\opal}
-use
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--help}
-}
-
-\vspace{1mm}
-
-\noindent \textbf{\textsf{Example:}} \texttt{\opal\_driver --help}
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Tetrahedral mesh file}
-\label{section_tetrahedral_mesh_file}
-
-To specify the path to the file that contains the input tetrahedral volume
-and triangular surface mesh for the computation:
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--mesh=\replace{string}}
-}
-
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}} \texttt{--mesh=cavity.h5}
-
-\vspace{2mm}
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Eigenvalue shift}
-\label{section_eigenvalue_shift}
-
-To specify the eigenvalue shift $\sigma$ used in the
-iterative solution of the eigenvalue computation
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--sigma=\replace{double}}
-}
-
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}} \texttt{--sigma=2.0}
-
-\vspace{2mm}
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Number of computed eigenmodes}
-\label{section_numbero_of_computed_eigenmodes}
-
-To specify the desired number of eigenmodal solution (eigenmodes)
-to be computed, specify:
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--kmax=\replace{integer}}
-}
-
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}} \texttt{--kmax=10}
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Nota bene:}} the \textsf{\opal} code
-uses an iterative scheme to compute eigenmodes, i.e. a variable
-number of iterations is required to obtain an eigenmode with
-given accuracy (cf. parameter \texttt{--tol} parameter). By default
-$5$ eigenmodes are computed for which the default number of
-$200$iterations (cf. parameter \texttt{--jitmax}) is sufficient.
-%%
-%%
-If more eigensolutions are required, then \texttt{--jitmax} should
-be increased. The increase depends on mesh quality, complexity
-of the problem and the number of desired modes.
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Numerical precision}
-\label{section_numerical_precision}
-
-To specify the accuracy for the computation of
-the eigenmodes, write
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--tol=\replace{double}} ($1.0 \cdot 10^{-8}$)
-}
-
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}} \texttt{--tol=1.0e-6} which is often
-good enough or a good start to explore the eigenmodal space of a resonator
-structure.
-
-\vspace{2mm}
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Finite element approximation order}
-\label{section_finite_element_approximation_order}
-
-
-At present $1^{\mathrm{st}}$ and $2^{\mathrm{nd}}$ order
-base functions are available in \textsf{\opal}.
-%%
-%%
-To select either of them, write
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--order=1}
-}
-
-\vspace{2mm}
-
-\noindent or
-
-\vspace{2mm}
-
-\fbox{
-\texttt{--order=2}
-}
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Caveat:}} While $1^{\mathrm{st}}$
-order base functions are robust and reliable, computations
-using $2^{\mathrm{nd}}$ order base functions may depend
-on geometry and the calculation may not converge.
-
-\vspace{2mm}
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Preconditioner selection}
-\label{section_preconditioner_selection}
-
-\textsf{\opal} offers a choice of different preconditioners,
-used in the iterative solution of the eigenmodal system.
-%%
-%%
-Possible selections are, where \texttt{|} denotes the logical \texttt{OR}
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--aprec=ml | neumann | lu | diag | 2level | if } (ml)
-}
-
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Available preconditioner types:}}
-
-\noindent \texttt{--aprec=neumann} : requests a Neumann type preconditioner
-
-\noindent \texttt{--aprec=ml} : requests a multilevel type preconditioner
-
-\noindent \texttt{--aprec=lu} : requests a LU type preconditioner
-
-\noindent \texttt{--aprec=2level} : requests a 2level type preconditioner
-
-\noindent \texttt{--aprec=diag} : requests a diagonal type preconditioner
-
-\noindent \texttt{--aprec=if} : requests a if type preconditioner
-
-
-\vspace{2mm}
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Switch off log file}
-\label{section_switch_of_log_file}
-
-By default, a \textsf{\opal} run produces a log file with information
-on all aspects of the eigenmodal computation. To switch log file generation
-off, write
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--no-logfile}
-}
-
-
-\vspace{2mm}
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Symmetry plane configuration}
-\label{section_symmetry_plane_configuration}
-
-Often, the geometry of electromagnetic resonators exhibits symmetries which can
-be exploited to reduce the computational effort considerably.
-%%
-A computational mesh in \textsf{\opal} compliand HDF5 file also contains
-a list of surface triangles, each of which is attributed a surface tag whose
-values are either $1,2$ or $3$.
-%%
-%%
-The surface triangle tags, $1$, $2$ or $3$, denote the lowest three bits in a binary number;
-depending on the bits being either set ($1$) or not set ($0$),
-the resulting binary number can express decimal values in the range of
-$0$ through $7$.
-%%
-%%
-If the corresponding bit of a surface triangle tag is set, then triangles with this
-tag are assigned the perfect electric conductor (PEC) boundary condition.
-%%
-If the bit is not set, then triangles with this tag are assigned the perfect
-magnetic conductor (PMC) boundary condition, in short:
-
-\noindent Bit $=0$ : $\vec{E} \times \vec{n} = 0$ (PEC)
-
-\noindent Bit $=1$ : $\vec{E} \cdot \vec{n} = 0$ (PMC)
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--sym-plane-config=\replace{integer}}
-}
-
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}} \texttt{--sym-plane-config=5}
-%%
-means that on surface triangles with
-tag $1$ we have a PMC boundary condition, on surface triangles with tag $2$ we
-have a PEC boundary condition, and on surface triangles with surface tag $3$ we
-have, once again, a PMC boundary condition.
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Debug mode}
-\label{section_debug_mode}
-
-To switch on debug mode during code execution, use
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--debug}
-}
-
-
-\vspace{2mm}
-
-
-
-
-
-
-
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Result file}
-\label{section_result_file}
-
-To specify the name of the HDF5 file into which the eigenmodal
-calculations' results are written
-
-
-\vspace{1mm}
-
-\fbox{
-\texttt{--eigendatafile=\replace{string}}
-}
-
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}} \texttt{--eigendatafile=box\_output.h5}
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Cartesian sampling}
-\label{section_result_file}
-
-To specify an orthogonal cartesian region within which
-the eigenmodal solution is sampled, use
-
-\vspace{1mm}
-
-\begin{verbatim}
---cartesian-sampling-x-start = <double>
---cartesian-sampling-x-stop = <double>
---number-samples-x-axis = <integer>
---cartesian-sampling-y-start = <double>
---cartesian-sampling-y-stop = <double>
---number-samples-y-axis = <integer>
---cartesian-sampling-z-start = <double>
---cartesian-sampling-z-stop = <double>
---number-samples-z-axis = <integer>
-\end{verbatim}
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}}
-
-\begin{verbatim}
---cartesian-sampling-x-start = 0.0
---cartesian-sampling-x-stop = 1.0
---number-samples-x-axis = 11
---cartesian-sampling-y-start = 0.5
---cartesian-sampling-y-stop = 1.5
---number-samples-y-axis = 21
---cartesian-sampling-z-start = -1.5
---cartesian-sampling-z-stop = -0.5
---number-samples-z-axis = 41
-\end{verbatim}
-
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Cylinder sampling}
-\label{section_result_file}
-
-Specifying an orthogonal cyliner region within which
-the eigenmodal solution is sampled, is similar to the
-cartesian region selection. The cylindric region is
-defined through the axis of the sampling cylinder volume
-and the end position of its radius vector, w.r.t cylinder
-axis' start location.
-
-
-
-\begin{verbatim}
---cylinder-sampling-x-bottom = <double>
---cylinder-sampling-x-top = <double>
---cylinder-sampling-y-bottom = <double>
---cylinder-sampling-y-top = <double>
---cylinder-sampling-z-bottom = <double>
---cylinder-sampling-z-top = <double>
---number-samples-on-cylinder-axis = <integer>
---cylinder-radial-vector-end-position-x = <double>
---cylinder-radial-vector-end-position-y = <double>
---cylinder-radial-vector-end-position-z = <double>
---number-samples-on-cylinder-radial-vector = <integer>
---number-samples-in-cylinder-azimuthal-plane = <integer>
-\end{verbatim}
-
-\vspace{2mm}
-
-\noindent \textbf{\textsf{Example:}}
-
-\begin{verbatim}
---cylinder-sampling-x-bottom = 0.5
---cylinder-sampling-x-top = 0.5
---cylinder-sampling-y-bottom = 0.7
---cylinder-sampling-y-top = 0.7
---cylinder-sampling-z-bottom = 0.0
---cylinder-sampling-z-top = 1.6
---number-samples-on-cylinder-axis = 171
---cylinder-radial-vector-end-position-x = 1.0
---cylinder-radial-vector-end-position-y = 1.3
---cylinder-radial-vector-end-position-z = 0.0
---number-samples-on-cylinder-radial-vector = 11
---number-samples-in-cylinder-azimuthal-plane = 4
-\end{verbatim}
-
-\noindent which specifies a cylinder sampling volume where the cylinder's
-axis starts at $[x=0.5,y=0.7,z=0.0]$ and ends at $[x=0.5,y=0.7,z=1.6]$ and
-the end position of th radius vector is $[x=1.0,y=1.3,z=0.0]$ with
-as many as $[171,11,4]$ samples on the cylinder axis, the radius vector
-and on the cylinder azimuthal plane, respectively.
-
-
-{\scriptsize
-\bibliography{opal-ref-card-bibliography}
-}
-
-\vspace{2mm}
-
-\scriptsize
-%%
-%%
-For ontacts:
-%%
-\texttt{andreas.adelmann@psi.ch}.
-%%
-%%
-\noindent \textbf{DISCLAIMER} All data in this reference card have been collected and presented
-with great care. However, no responsibility is accepted for any consequences resulting
-from errors or omissions. Neither do we take any responsibility for using the codes.
-%%
-%%
-We encourage you to report bugs, problems and suggestions to us.
-%%
-%%
-\noindent On the internet type \texttt{www.amas.web.psi.ch} to find us.
-%%
-\end{multicols}
-\end{document}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/src/Algorithms/ThickTracker.cpp b/src/Algorithms/ThickTracker.cpp
index c32aa92e2..7a69420b5 100644
--- a/src/Algorithms/ThickTracker.cpp
+++ b/src/Algorithms/ThickTracker.cpp
@@ -18,6 +18,20 @@
//
// ------------------------------------------------------------------------
+
+/*
+#include <cmath>
+#include <exception>
+#include <functional>
+#include <iostream>
+#include <string>
+#include <vector>
+#include <fstream>
+*/
+
+
+#include <typeinfo>
+#include <fstream>
#include "Algorithms/ThickTracker.h"
#include "Algorithms/OrbitThreader.h"
#include "Algorithms/CavityAutophaser.h"
@@ -31,13 +45,24 @@
#include "Beamlines/FlaggedBeamline.h"
#include "Fields/BMultipoleField.h"
-#include "FixedAlgebra/FTps.h"
-#include "FixedAlgebra/FTpsMath.h"
-#include "FixedAlgebra/FVps.h"
+#include "Classic/FixedAlgebra/FTps.h"
+#include "Classic/FixedAlgebra/FTpsMath.h"
+#include "Classic/FixedAlgebra/FVps.h"
+
+#include "Classic/Fields/BSingleMultipoleField.h"
+
+
+#include "Classic/Algorithms/PartData.h" //for the beam reference
+
#include "Physics/Physics.h"
//#include "Utilities/NumToStr.h"
+
#include "Elements/OpalBeamline.h"
+#include <Classic/BeamlineCore/MultipoleRep.h>
+#include "Distribution/MapGenerator.h"
+
+#define DIM 3
class Beamline;
class PartData;
@@ -269,7 +294,6 @@ void ThickTracker::changeDT() {
itsBunch_m->dt[i] = itsBunch_m->getdT();
}
}
-
void ThickTracker::findStartPosition(const BorisPusher &pusher) {
double t = 0.0;
@@ -328,128 +352,496 @@ void ThickTracker::findStartPosition(const BorisPusher &pusher) {
}
}
+/**
+ * @brief Algorithm for Thick Map-Tracking
+ */
+
void ThickTracker::execute() {
- Inform msg("ThickTracker ", *gmsg);
- msg << "in execute " << __LINE__ << " " << __FILE__ << endl;
+ std::ofstream outfile;
+ std::ofstream tmap;
+ outfile.open ("/home/phil/Documents/ETH/MScProj/OPAL/src/tests/Maps/generatedMaps.txt");
+ tmap.open ("/home/phil/Documents/ETH/MScProj/OPAL/src/tests/Maps/TransferMap.txt");
+ outfile << std::setprecision(20);
- /*
+ Inform msg("ThickTracker", *gmsg);
+
+ msg << "in execute " << __LINE__ << " " << __FILE__ << endl;
+
+ /*
First some setup and general preparation. Mostly copied from ParalellTTracker.
Some of them we maybe do not need at all.
- */
+ */
- OpalData::getInstance()->setInPrepState(true);
+ OpalData::getInstance()->setInPrepState(true);
- BorisPusher pusher(itsReference);
- OpalData::getInstance()->setGlobalPhaseShift(0.0);
+ BorisPusher pusher(itsReference);
+ OpalData::getInstance()->setGlobalPhaseShift(0.0);
- dtCurrentTrack_m = itsBunch_m->getdT();
+ dtCurrentTrack_m = itsBunch_m->getdT();
- if (OpalData::getInstance()->hasPriorTrack() || OpalData::getInstance()->inRestartRun()) {
+ if (OpalData::getInstance()->hasPriorTrack() || OpalData::getInstance()->inRestartRun()) {
Options::openMode = Options::APPEND;
- }
-
- prepareSections();
+ }
- itsOpalBeamline_m.compute3DLattice();
- itsOpalBeamline_m.save3DLattice();
- itsOpalBeamline_m.save3DInput();
+ prepareSections();
- std::queue<double> timeStepSizes(dtAllTracks_m);
- std::queue<unsigned long long> numSteps(localTrackSteps_m);
- double minTimeStep = timeStepSizes.front();
- unsigned long long totalNumSteps = 0;
- while (timeStepSizes.size() > 0) {
- if (minTimeStep > timeStepSizes.front()) {
- totalNumSteps = std::ceil(totalNumSteps * minTimeStep / timeStepSizes.front());
- minTimeStep = timeStepSizes.front();
- }
- totalNumSteps += std::ceil(numSteps.front() * timeStepSizes.front() / minTimeStep);
-
- numSteps.pop();
- timeStepSizes.pop();
- }
-
- itsOpalBeamline_m.activateElements();
-
- if (OpalData::getInstance()->hasPriorTrack() ||
- OpalData::getInstance()->inRestartRun()) {
-
- referenceToLabCSTrafo_m = itsBunch_m->toLabTrafo_m;
- RefPartR_m = referenceToLabCSTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
- RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
-
- pathLength_m = itsBunch_m->get_sPos();
- zstart_m = pathLength_m;
-
- restoreCavityPhases();
- } else {
- RefPartR_m = Vector_t(0.0);
- RefPartP_m = euclidean_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
-
- if (itsBunch_m->getTotalNum() > 0) {
- if (!itsOpalBeamline_m.containsSource()) {
- RefPartP_m = OpalData::getInstance()->getP0() / itsBunch_m->getM() * Vector_t(0, 0, 1);
- }
-
- if (zstart_m > pathLength_m) {
- findStartPosition(pusher);
- }
-
- itsBunch_m->set_sPos(pathLength_m);
- }
- }
- Vector_t rmin, rmax;
- itsBunch_m->get_bounds(rmin, rmax);
-
- OrbitThreader oth(itsReference,
- referenceToLabCSTrafo_m.transformTo(RefPartR_m),
- referenceToLabCSTrafo_m.rotateTo(RefPartP_m),
- pathLength_m,
- -rmin(2),
- itsBunch_m->getT(),
- minTimeStep,
- totalNumSteps,
- zStop_m.back() + 2 * rmax(2),
- itsOpalBeamline_m);
-
- oth.execute();
-
+// itsOpalBeamline_m.compute3DLattice();
+// itsOpalBeamline_m.save3DLattice();
+// itsOpalBeamline_m.save3DInput();
+//
+// std::queue<double> timeStepSizes(dtAllTracks_m);
+// std::queue<unsigned long long> numSteps(localTrackSteps_m);
+// double minTimeStep = timeStepSizes.front();
+// unsigned long long totalNumSteps = 0;
+// while (timeStepSizes.size() > 0) {
+// if (minTimeStep > timeStepSizes.front()) {
+// totalNumSteps = std::ceil(totalNumSteps * minTimeStep / timeStepSizes.front());
+// minTimeStep = timeStepSizes.front();
+// }
+// totalNumSteps += std::ceil(numSteps.front() * timeStepSizes.front() / minTimeStep);
+//
+// numSteps.pop();
+// timeStepSizes.pop();
+// }
+//
+// itsOpalBeamline_m.activateElements();
+//
+// if (OpalData::getInstance()->hasPriorTrack() ||
+// OpalData::getInstance()->inRestartRun()) {
+//
+// referenceToLabCSTrafo_m = itsBunch_m->toLabTrafo_m;
+// RefPartR_m = referenceToLabCSTrafo_m.transformFrom(itsBunch_m->RefPartR_m);
+// RefPartP_m = referenceToLabCSTrafo_m.rotateFrom(itsBunch_m->RefPartP_m);
+//
+// pathLength_m = itsBunch_m->get_sPos();
+// zstart_m = pathLength_m;
+//
+// restoreCavityPhases();
+// } else {
+// RefPartR_m = Vector_t(0.0);
+// RefPartP_m = euclidean_norm(itsBunch_m->get_pmean_Distribution()) * Vector_t(0, 0, 1);
+//
+// if (itsBunch_m->getTotalNum() > 0) {
+// if (!itsOpalBeamline_m.containsSource()) {
+// RefPartP_m = OpalData::getInstance()->getP0() / itsBunch_m->getM() * Vector_t(0, 0, 1);
+// }
+//
+// if (zstart_m > pathLength_m) {
+// findStartPosition(pusher);
+// }
+//
+// itsBunch_m->set_sPos(pathLength_m);
+// }
+// }
+//
+// Vector_t rmin, rmax;
+// itsBunch_m->get_bounds(rmin, rmax);
+//
+// OrbitThreader oth(itsReference,
+// referenceToLabCSTrafo_m.transformTo(RefPartR_m),
+// referenceToLabCSTrafo_m.rotateTo(RefPartP_m),
+// pathLength_m,
+// -rmin(2),
+// itsBunch_m->getT(),
+// minTimeStep,
+// totalNumSteps,
+// zStop_m.back() + 2 * rmax(2),
+// itsOpalBeamline_m);
+//
+// oth.execute();
+//
/*
End of setup and general preparation.
*/
- msg << *itsBunch_m << endl;
+ msg << *itsBunch_m << endl;
- /*
+ /*
This is an example how one can loop
over all elements.
- */
-
- auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
-
- FieldList::iterator it = allElements.begin();
-
- const FieldList::iterator end = allElements.end();
-
- if (it == end)
- msg << "No element in lattice" << endl;
-
- for (; it != end; ++ it) {
- std::shared_ptr<Component> element = (*it).getElement();
- msg << "Element name " << element->getName() << endl;
+ */
+
+//=======================================================================
+
+
+ //generate the Hamiltonian
+ typedef FTps<double, 2 * DIM> Series;
+ typedef FVps<double, 2 * DIM> Map;
+
+ Series x = Series::makeVariable(0); //SIXVect::X);
+ Series px = Series::makeVariable(1); //SIXVect::PX);
+ Series y = Series::makeVariable(2); //SIXVect::Y);
+ Series py = Series::makeVariable(3); //SIXVect::PY);
+ //Series z = Series::makeVariable(4); //SIXVect::TT);
+ Series delta = Series::makeVariable(5); //SIXVect::PT);
+
+
+
+//Diese Parameter muessen per & uebernommen werden
+//======================================================
+
+ int order = 2; //actually already defined
+// define parameter for Hamiltonian
+//
+// double b = 1.; // field gradient [T/m]
+// double E = 200.; //PartData::getE(); //beam energy [MeV]
+// double ds = 0.4; //stepsize [m]
+//
+// double EProt = refPart.getM(); //938.2720813; //MeV/c
+
+ double gamma0 = (itsBunch_m->getInitialGamma()) ; //(EProt + E) / EProt;
+ double beta0 = (itsBunch_m->getInitialBeta()); //std::sqrt(gamma0 * gamma0 - 1.0) / gamma0;
+
+ double P0 = (itsBunch_m-> getP()); //beta0 * gamma0 * EProt * 1e6 / Physics::c;
+ double q = (itsBunch_m->getQ()); // particle change [e]
+
+//=====================================================
+ double E=(itsReference.getE() - itsReference.getM());
+ E = std::round(E);
+ double AMU = 1.66053873e-27;
+ double EZERO = 1.602176462e-19 ;
+ double CLIGHT= 2.99792458e8 ;
+ double AMUMEV = AMU*(CLIGHT*CLIGHT)/EZERO ;
+ double M0= 1.00727646688;
+
+
+ double ETA = E/(M0*AMUMEV) ;
+ double PHI = std::sqrt(ETA*(2+ETA)) ;
+
+ P0 = (AMUMEV*M0)*PHI;
+ gamma0=ETA+1;
+ beta0= std::sqrt(1-1/(gamma0*gamma0));
+
+
+//=====================================================
+
+
+ beta0 = (PHI / (1 + ETA));
+ gamma0 = PHI / beta0;
+
+ msg << std::setprecision(20);
+ msg << "BreamEnergy: "<< E << endl;
+ msg << "P0 COSY: " << P0 << endl;
+ msg << "E0 COSY: " << AMUMEV*M0 << endl;
+//=====================================================
+ Series::setGlobalTruncOrder(order);
+
+// double r0 = 1.; // get aperture, resp. the radius of magnet [m]
+
+
+ ///Creates the Hamiltonian for beam line element
+ /**\param element iterative pointer to the elements along the beam line
+ * \return Hamiltonian times elemenrLength
+ */
+ auto Hamiltonian = [&](std::shared_ptr<Component> element) {
+
+ Series H; //Hamiltonian
+ Series Hs;//Hamiltonian times ElementLength
+
+ switch(element->getType()) {
+
+ /**Driftspace
+ * \f[H_{Drift}= \frac{\delta}{\beta_0} -
+ * \sqrt{\left(\frac{1}{\beta_0} + \delta \right)^2 -p_x^2 -p_y^2 - \frac{1}{\left(\beta_0 \gamma_0\right)^2 } } \f]
+ */
+ case ElementBase::ElementType::DRIFT: {
+ Drift* pDrift= dynamic_cast<Drift*> (element.get());
+ //outfile << "0T, ";
+ //outfile << pDrift ->getElementLength() << "m" <<std::endl;
+ H=( delta / beta0 )
+ - sqrt((1./ beta0 + delta ) *(1./ beta0 + delta )
+ - ( px*px )
+ - ( py*py )
+ - 1./( beta0 * beta0 * gamma0 * gamma0 ),order
+ );
+ Hs = H * pDrift->getElementLength();
+ break;
+ }
+
+ /**Rectangular Bend
+ * \f[H_{Dipole}= \frac{\delta}{\beta_0} - \left( 1+ hx \right)
+ * \sqrt{\left(\frac{1}{\beta_0} + \delta \right)^2 -p_x^2 -p_y^2 - \frac{1}{\left(\beta_0 \gamma_0\right)^2 } } +
+ * \left( 1+ hx \right) k_0 \left(x - \frac{hx^2}{2 \left( 1+ hx \right)}\right) \f]
+ */
+ case ElementBase::ElementType::RBEND: {
+ RBend* pRBend= dynamic_cast<RBend*> (element.get());
+
+// double rho = P0 / (b * q); //bending radius [m]
+ double h = 1. / pRBend ->getBendRadius(); //inverse bending radius [1/m]
+ double K0= pRBend ->getB()*(Physics::c/itsReference.getP());
+
+// double phi = 45 * 0.0174532925; // tilt angle for dipole fringe field [rad]
+
+
+ //outfile << K0 << "T, ";
+ //outfile << pRBend ->getArcLength() << "m" <<std::endl;
+ H=( delta / beta0 )
+ - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+ - ( px*px )
+ - ( py*py )
+ - 1./( beta0*beta0 * gamma0*gamma0 ),order
+ ))
+ - (h * x)
+ * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+ - ( px*px )
+ - ( py*py )
+ - 1./( beta0*beta0 * gamma0*gamma0 ),order-1
+ ))
+ + K0 * x * (1. + 0.5 * h* x);
+
+ Hs = H * pRBend->getArcLength();
+
+ break;
+ }
+
+ /**Sector Bend
+ * \f[H_{Dipole}= \frac{\delta}{\beta_0} - \left( 1+ hx \right)
+ * \sqrt{\left(\frac{1}{\beta_0} + \delta \right)^2 -p_x^2 -p_y^2 - \frac{1}{\left(\beta_0 \gamma_0\right)^2 } } +
+ * \left( 1+ hx \right) k_0 \left(x - \frac{hx^2}{2 \left( 1+ hx \right)}\right) \f]
+ */
+ case ElementBase::ElementType::SBEND: {
+
+ SBend* pSBend= dynamic_cast<SBend*> (element.get());
+// double rho = P0 / (b * q); //bending radius [m]
+ double h = 1. / pSBend ->getBendRadius(); //inverse bending radius [1/m]
+
+ double K0= pSBend ->getB()*(Physics::c/itsReference.getP());
+ msg<< "Magnetic Field Sbend:" << pSBend ->getB() << "T" << endl;
+// double phi = 45 * 0.0174532925; // tilt angle for dipole fringe field [rad]
+ msg<< "Effecgtive Length:" << pSBend ->getEffectiveLength() << "m" << endl;
+ //outfile << pSBend ->getB() << "T, ";
+ //outfile << pSBend ->getEffectiveLength() << "m" <<std::endl;
+ H=( delta / beta0 )
+ - (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+ - ( px*px )
+ - ( py*py )
+ - 1./( beta0*beta0 * gamma0*gamma0 ),order
+ ))
+ - (h * x)
+ * (sqrt ((1./ beta0 + delta) *(1./ beta0 + delta)
+ - ( px*px )
+ - ( py*py )
+ - 1./( beta0*beta0 * gamma0*gamma0 ),order-1
+ ))
+ + K0 * x * (1. + 0.5 * h* x);
+
+ Hs= H * pSBend->getEffectiveLength();
+
+ break;
+ }
+
+ /**Quadrupole "in Multipolegroup"
+ * \f[H_{Quadrupole}= \frac{\delta}{\beta_0} -
+ * \sqrt{\left(\frac{1}{\beta_0} + \delta \right)^2 -p_x^2 -p_y^2 - \frac{1}{\left(\beta_0 \gamma_0\right)^2 } } +
+ * \frac{1}{2} k_1 \left( x^2 - y^2 \right) \f]
+ */
+ case ElementBase::ElementType::MULTIPOLE: {
+ Multipole* pMultipole= dynamic_cast<Multipole*> (element.get());
+
+ double K1= pMultipole->getField().getNormalComponent(2)*(Physics::c/P0);
+ K1= std::round(K1*1e6)/1e6;
+// outfile << K1 << "T/m, ";
+// outfile << pMultipole ->getElementLength() << "m, " ;
+// outfile << P0 << "MeV"<<std::endl;
+// K1 = (q * b) / (P0 * r0)
+
+ H= ( delta / beta0 )
+ - sqrt ((1./ beta0 + delta ) *(1./ beta0 + delta)
+ - ( px*px )
+ - ( py*py )
+ - 1./( beta0*beta0 * gamma0*gamma0 ),order
+ )
+ + 0.5 * (q * K1) *(Physics::c/P0) * (x*x - y*y);
+
+ Hs=H* pMultipole ->getElementLength();
+
+ break;
+ }
+
+ default:
+ throw LogicalError("ThickTracker::exercute,",
+ "Please use already defined beam line element: At this time just driftspace, multipoles and dipoles");
+ break;
+
+ }
+ return Hs;
+ };
+
+
+//=======================================================================
+
+
+
+ FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
+ struct sort_by_pos {
+ bool operator()(const ClassicField &a, const ClassicField &b)
+ {return a.getElement()->getElementPosition() < b.getElement()->getElementPosition();}
+ };
+
+ allElements.sort(sort_by_pos());
+
+
+
+ FieldList::iterator it = allElements.begin();
+ const FieldList::iterator end = allElements.end();
+ if (it == end) msg << "No element in lattice" << endl;
+
+ Map combinedMaps;
+
+
+
+
+
+// msg << "beam parameter " << itsReference.getE() - itsReference.getM() << endl;
+// msg << "beam line length: " << allElements.size() << endl;
+// //const int el = allElements.size();
+ msg << std::setprecision(20);
+ msg << "P0 OPAL: " << itsBunch_m-> getP() << endl;
+ msg << "E0 OPAL: " << itsBunch_m-> getM() << endl;
+
+
+
+// msg << "q: " << q << endl;
+
+
+ std::map<std::string, double> elementDic;
+ std::map<std::string, double>::iterator dicit;
+
+ std::vector<Map> mapVec;
+ double selement =0;
+ for (; it != end; ++ it) {
+ dicit = elementDic.begin();
+
+ std::shared_ptr<Component> element = (*it).getElement();
+
+
+ dicit= elementDic.find(element->getName());
+ if (dicit !=elementDic.end()){
+ throw LogicalError("ThickTracker::execute,",
+ "Same Element twice in beamline:"+ element->getName());
+
+ }else{
+ elementDic.insert(std::pair<std::string, double> (element->getName(), element->getElementPosition()));
+ }
+
+ if (selement > element->getElementPosition()){
+ msg << "There is an overlap! @ Element: " << element ->getName() <<
+ "-> starts at: " << element->getElementPosition()<<
+ " the previous ends at: " << selement << endl;
+// throw LogicalError("ThickTracker::exercute,",
+// "Overlap at element:"+ element->getName());
+
+ }else{
+ selement=element->getElementPosition()+ element ->getElementLength();
+ }
+
+ if (element ->getType()==ElementBase::MONITOR){
+ mapVec.insert(mapVec.end(), combinedMaps);
+ combinedMaps.identity();
+
+ }
+
+
+ // outfile << "Element name: " << element->getName() << std::endl;
+ // outfile << "Element type: " << element->getType()<< std::endl;
+ msg << "=============================="<< endl;
+ msg << "Element name: " << element->getName() << endl;
+ msg << "Element Length: " << element->getElementLength() << endl;
+// msg << "Element nummer: " << element->getType() << endl;
+ msg << "ElementPosition"<< element->getElementPosition()<< endl;
+ msg << "EntryPoint" << element->getEntrance()<< endl;
+
+
+// outfile << "Element Hamiltonian: \n" << Hamiltonian(element) << std::endl;
+// outfile << element->getName() << ", ";
+// outfile << E/1.0e6 << "MeV, ";
+
+ Map ElementMap=ExpMap(-Hamiltonian(element) ,order);
+// outfile << Hamiltonian(element);
+// outfile << "\n=========="<< std::endl;
+ //outfile << ElementMap << std::endl;
+ combinedMaps = ElementMap * combinedMaps;
+ //outfile << "CombinedMap: \n" << combinedMaps << std::endl;
+ //outfile << "------------------------------------------------------------" << std::endl;
+
+
}
-}
+ tmap << "A customized FODO" << std::endl;
+ tmap << combinedMaps << std::endl;
+// outfile << combinedMaps << std::endl;
+
+ FVector<double, 6> particle;
+ OpalParticle part;
+ outfile <<"P0: " << itsBunch_m->getP()/(M0*AMUMEV) << " "
+ << itsBunch_m->getInitialGamma()*itsBunch_m->getInitialBeta()<< std::endl;
+
+
+ for (unsigned int idx=0; idx< itsBunch_m->getTotalNum(); idx++){
+ part =itsBunch_m ->get_part(idx);
+
+ for(int dim=0; dim<6; dim++){
+ particle[dim]=itsBunch_m ->get_part(idx)[dim];
+ }
+ outfile << "=======================" << std::endl;
+ outfile << "Particle#: " << idx << " "<<itsBunch_m->getBeta(idx) << std::endl;
+ outfile << particle;
+// outfile << "particle[5]: " << particle[5] << std::endl;
+// outfile << "1st term= " << (particle[5]*M0*AMUMEV/ itsBunch_m->getP()) << std::endl;
+// outfile << "1st term= " << (particle[5]*itsBunch_m->getM()/ itsBunch_m->getP()) << std::endl;
+ particle[5]=
+ (particle[5]*itsBunch_m->getM()/itsBunch_m->getP())
+ * std::sqrt( 1./(particle[5]* particle[5]) +1)
+ -1./itsBunch_m->getInitialBeta();
+
+// outfile << "---> Unit Transformation <--- \n" << particle;
+ outfile << "---> The Map is applied <---"<< std::endl;
+ particle= combinedMaps * particle;
+// outfile << "particle[5]" << particle[5] << std::endl;
+//
+// outfile << "divident: " << particle[5]*M0*AMUMEV << std::endl;
+// outfile << "divisor: " << itsBunch_m->getP() << std::endl;
+// outfile << "M0: " << M0*AMUMEV << "OPAL -> "<< itsBunch_m->getM() << std::endl;
+// outfile << "P0: " << itsBunch_m->getP() << std::endl;
+// outfile << "gamma0: " << itsBunch_m->getInitialGamma() << std::endl;
+// outfile << "delta: " << particle[5] << std::endl;
+ particle[4]= selement;
+ particle[5]= (particle[5] + 1./itsBunch_m->getInitialBeta()) * itsBunch_m->getP()/itsBunch_m->getM()
+ /std::sqrt( 1./(itsBunch_m ->get_part(idx)[5]* itsBunch_m ->get_part(idx)[5]) +1) ;
+
+ outfile << particle << std::endl;
+ }
+
+
+
+
+ msg << part.y() << endl;
+
+ msg<< "number of Particles" << itsBunch_m->getTotalNum() << endl;
+
+
+ for(int i=0; i<6; i++){
+ msg << "OPALPart" << part[i]<< endl;
+ msg << "OPALBart" << particle[i]<< endl;
+ }
+
+ std::map <std::string, double>::iterator mi;
+ for (mi=elementDic.begin(); mi!=elementDic.end(); ++mi)
+ msg << mi->first << " = " << mi->second << endl;
+
+}
+
+
/*
namespace {
diff --git a/src/Classic/Algorithms/PartData.h b/src/Classic/Algorithms/PartData.h
index fbf626c80..735d54416 100644
--- a/src/Classic/Algorithms/PartData.h
+++ b/src/Classic/Algorithms/PartData.h
@@ -97,7 +97,7 @@ protected:
double charge; // Particle charge.
double mass; // Particle mass.
double beta; // particle velocity divided by c.
- double gamma; // particle energy divided by particle mas
+ double gamma; // particle energy divided by particle mass
};
diff --git a/src/Elements/OpalBeamline.h b/src/Elements/OpalBeamline.h
index c1181cae5..b90add286 100644
--- a/src/Elements/OpalBeamline.h
+++ b/src/Elements/OpalBeamline.h
@@ -303,4 +303,4 @@ inline
bool OpalBeamline::containsSource() {
return containsSource_m;
}
-#endif // OPAL_BEAMLINE_H
\ No newline at end of file
+#endif // OPAL_BEAMLINE_H
diff --git a/tests/Maps/map.in b/tests/Maps/map.in
index 78d05dede..fbbf3c8f4 100644
--- a/tests/Maps/map.in
+++ b/tests/Maps/map.in
@@ -4,7 +4,7 @@ OPTION, VERSION=10900;
TITLE, STRING="MAP Test";
-REAL Edes = 1.0;
+REAL Edes = 0.10;
REAL gamma = (Edes+PMASS)/PMASS;
REAL beta = sqrt(1-(1/gamma^2));
REAL gambet = gamma*beta;
@@ -14,13 +14,23 @@ REAL rf = 50.6328e6; //need to be confirmed
VALUE,{gamma,brho,Edes,beta,gambet};
-QP1: QUADRUPOLE, L=1.0, K1= 1.0, ELEMEDGE=0.000;
-QP2: QUADRUPOLE, L=1.0, K1=-1.0, ELEMEDGE=0.500;
-QP3: QUADRUPOLE, L=1.0, K1= 1.0, ELEMEDGE=2.000;
+QP1: QUADRUPOLE, L=1.0, K1= 2.00, ELEMEDGE=0.000;
+QP2: QUADRUPOLE, L=1.0, K1=-3.00, ELEMEDGE=1.000;
+QP3: QUADRUPOLE, L=1.0, K1= 4.00, ELEMEDGE=2.000;
-D: DRIFT, L=2.0, ELEMEDGE=3.000;
-QUADTEST: LINE=(QP1,QP2,QP3,D);
+D: DRIFT, L=2.0, ELEMEDGE=3.000 ;
+
+//BENDS: RBend, K0 = 1.000, FMAPFN = "1DPROFILE1-DEFAULT", ELEMEDGE = 5.0, DESIGNENERGY = Edes*1e3, L = 1.0, GAP = 0.02;
+
+BENDS: SBend, ANGLE = 30.0*Pi / 180.0, FMAPFN = "1DPROFILE1-DEFAULT", ELEMEDGE = 5.0, DESIGNENERGY = Edes*1e3, L = 1.0, GAP = 0.02;
+
+
+//As an example for a sextupole
+S11: SEXTUPOLE, L=0.4, K2=0.00134, ELEMEDGE=4.000;
+
+//QUADTEST: LINE=( QP1, QP2, D, QP3);
+QUADTEST: LINE=( QP1 );
D1: DISTRIBUTION, TYPE=GAUSS,
SIGMAX= 5*1.0e-03, SIGMAPX= 0.696*1.0e06, CORRX= 0.0,
@@ -38,7 +48,8 @@ BEAM1: BEAM, PARTICLE=PROTON, PC=P0, NPART=1000, BCURRENT=2.0e-03, BFREQ=rf, CHA
SELECT, LINE=QUADTEST;
TRACK, LINE=QUADTEST, BEAM=BEAM1, MAXSTEPS=1, DT=1.0e-11, ZSTOP=8.0;
- RUN, METHOD = "THICK", BEAM=BEAM1, FIELDSOLVER=FS1, DISTRIBUTION=D1;
+
+RUN, METHOD = "THICK", BEAM=BEAM1, FIELDSOLVER=FS1, DISTRIBUTION=D1;
ENDTRACK;
-STOP;
\ No newline at end of file
+STOP;
--
GitLab