33 double det = -
pow(c_prim, 3) *
pow(s, 3) + 3 * c * s_prim *
pow(c_prim, 2)
34 *
pow(s, 2) - 3 *
pow(c, 2) * c_prim * s *
pow(s_prim, 2) +
pow(c, 3) *
pow(s_prim, 3);
36 M(0, 0) = (-c_prim * s *
pow(s_prim, 2) + c *
pow(s_prim, 3)) / det;
37 M(0, 1) = (-2 * c_prim *
pow(s, 2) * s_prim + 2 * c * s *
pow(s_prim, 2)) / det;
38 M(0, 2) = (-c_prim *
pow(s, 3) + c *
pow(s, 2) * s_prim) / det;
39 M(1, 0) = (-
pow(c_prim, 2) * s * s_prim + c * c_prim *
pow(s_prim, 2)) / det;
40 M(1, 1) = (-
pow(c_prim, 2) *
pow(s, 2) +
pow(c, 2) *
pow(s_prim, 2)) / det;
41 M(1, 2) = (-c * c_prim *
pow(s, 2) +
pow(c, 2) * s * s_prim) / det;
42 M(2, 0) = (-
pow(c_prim, 3) * s + c *
pow(c_prim, 2) * s_prim) / det;
43 M(2, 1) = (-2 * c *
pow(c_prim, 2) * s + 2 *
pow(c, 2) * c_prim * s_prim) / det;
44 M(2, 2) = (-
pow(c, 2) * c_prim * s +
pow(c, 3) * s_prim) / det;
58 int printWidth, printPrecision;
66 const ColDesc allColumns[] = {
79 const ColDesc defaultColumns[] = {
87 const ColDesc *findCol(
const Aperture &table,
const std::string &colName) {
88 for(
const ColDesc *col = allColumns; col->colName; ++col) {
89 if(colName == col->colName) {
96 "\" has no column named \"" + colName +
"\".");
111 Column(
const Aperture &tab,
const std::string &colName,
const ColDesc &desc);
113 Column(
const Column &);
120 virtual double evaluate()
const;
123 virtual void print(std::ostream &os,
int precedence = 99)
const;
129 const Column &operator=(
const Column &);
148 Column::Column(
const Aperture &tab,
const std::string &colName,
const ColDesc &desc):
149 itsTable(tab), colName(colName),
150 get(desc.get), ind_1(desc.ind_1), ind_2(desc.ind_2)
154 Column::Column(
const Column &rhs):
156 itsTable(rhs.itsTable), colName(rhs.colName),
157 get(rhs.get), ind_1(rhs.ind_1), ind_2(rhs.ind_2)
166 return new Column(*
this);
170 double Column::evaluate()
const {
171 return (itsTable.*
get)(itsTable.getCurrent(), ind_1, ind_2);
175 void Column::print(std::ostream &os,
int)
const {
183 Table(
SIZE,
"APERTURE",
"help"), itsTable() {
185 (
"TABLE",
"Name of the TWISS table to be used.");
187 (
"BEAM",
"The beam to be used",
"UNNAMED_BEAM");
189 (
"NSLICE",
"the number of slices of the interpolation inside the optic element", 10);
191 (
"STATIC",
"recalculation if static equal false",
true);
193 (
"DATA",
"The data needed for aperture calculation(co,deltap,BBeat,HALO)");
195 (
"DEFAULTAPERTURE",
"The default beam screen for markers and drift generated in sequences");
197 (
"FILE",
"Name of file to receive APERTURE output",
"APERTURE.dat");
203 Table(name, parent), itsTable(name)
207 Orb(a.Orb), Interpol(a.Interpol) {}
211 return Interpol[ind].Beta_x;
214 return Interpol[ind].Beta_y;
218 return Interpol[ind].Disp_x;
221 return Interpol[ind].Disp_x_prim;
224 return Interpol[ind].Disp_y;
227 return Interpol[ind].Disp_y_prim;
230 return Interpol[ind].apert;
258 for(
int i = 1;
i < nslice; ++
i) {
268 for(
int i = 1;
i < nslice; ++
i) {
279 for(
int i = 1;
i < nslice; ++
i) {
289 double K = f.
normal(2) * 0.299792458;
314 double K = f.
normal(2) * 0.299792458;
325 ERRORMSG(
"MSplit::visitCyclotron(const Cyclotron &cy) not implemented");
340 double K = f.
normal(2) * 0.299792458;
389 if(
data.
type.compare(
"rbend") == 0) {
437 }
else if(
data.
type.compare(
"sbend") == 0) {
488 for(
int j = 1; j <= nslice - 1; ++j) {
513 for(
int j = 1; j <= nslice - 1; ++j) {
534 for(
int j = 1; j <= nslice - 1; ++j) {
554 for(
int j = 1; j <= nslice - 1; ++j) {
576 for(
int j = 1; j <= nslice - 1; ++j) {
597 for(
int j = 1; j <= nslice - 1; ++j) {
618 for(
int j = 0; j < nslice; ++j) {
631 if(nam1 ==
"[DRIFT]") {
639 dat.erase(dat.begin(), dat.begin() + 4);
650 dat.erase(dat.begin(), dat.begin() + 4);
659 vec.second.erase(
vec.second.begin(),
vec.second.begin() + 4);
666 for(
int j = 0; j < nslice; ++j) {
675 for(
int j = 1; j < nslice; ++j) {
688 double *Pn =
new double[pt_halo], *Qn =
new double[pt_halo];
692 double *Apert =
new double[pt_halo*pt_delta];
694 double deltax, deltay;
696 for(
int k = 0; k < pt_delta; ++k) {
702 for(
int i = 0;
i < pt_halo; ++
i) {
705 double a_d = Qn[
i] / Pn[
i];
706 double b_d = deltay - a_d * deltax;
709 for(
int j = 0; j < size - 1; ++j) {
715 x = (b_d - b_arc) / (a_arc - a_d);
716 double y = a_d * x + b_d;
718 double ang =
atan2(y, x);
722 if((ang > ang1) && (ang <= ang2))
break;
724 Apert[
i+k*pt_halo] =
n1 * (x - deltax) / Pn[
i];
728 for(
int i = 1;
i < pt_halo * pt_delta; ++
i) {
746 if((r > -1.5) && (r < -0.5)) {
747 double ang =
atan2(v, h);
748 for(
int i = 0;
i < nb_pt; ++
i) {
759 }
else if((r > -10.5) && (r < -9.5)) {
762 while(iter != vec.end()) {
769 if((r != 0) && (h == 0) && (v == 0)) {
770 for(
int i = 0;
i < nb_pt; ++
i) {
775 }
else if((r != 0) && (h != 0) && (v == 0)) {
776 double ang =
acos(h / r);
777 for(
int i = 0;
i < nb_pt; ++
i) {
788 }
else if((r != 0) && (h == 0) && (v != 0)) {
789 double ang =
asin(v / r);
790 for(
int i = 0;
i < nb_pt; ++
i) {
801 }
else if((r != 0) && (h != 0) && (v != 0)) {
802 double ang1 =
acos(h / r);
803 double ang2 =
asin(v / r);
804 for(
int i = 0;
i < nb_pt; ++
i) {
820 for(
int i = 0;
i < nb_pt; ++
i) {
841 dat.erase(dat.begin(), dat.begin() + 4);
851 ofstream outFile(file_out.c_str());
853 cerr <<
"Aperture: unable to open output file:" << file_out <<
endl;
857 outFile.setf(std::ios::scientific, std::ios::floatfield);
861 outFile <<
"@ DATE %s " << timer.
date() <<
endl;
862 outFile <<
"@ TIME %s " << timer.
time() <<
endl;
863 outFile <<
"@ ORIGIN %s OPAL_9.5/7\n";
864 outFile <<
"@ COMMENT %s " <<
endl;
865 outFile <<
"@ NSLICE %e " << nslice <<
endl;
867 outFile <<
"*" << setw(15) <<
"NAME" << setw(15) <<
"Type" << setw(15) <<
"S" << setw(15) <<
"L" << setw(15) <<
"Apert" <<
endl;
868 outFile <<
"$" << setw(15) <<
"%23s" << setw(15) <<
"%12s" << setw(15) <<
"%e" << setw(15) <<
"%e" << setw(15) <<
"%e" <<
endl;
893 if(dat.size() != 0) {
896 const std::string nam =
i->getElement()->getName();
897 if(nam ==
"[DRIFT]") {
898 A_row row(*
i, static_cast<int>(nslice));
900 calcul(
i, row, static_cast<int>(nslice),
tp);
906 if((elem.
getApert().second.size() == 0) && (Typ !=
"MARKER"))
i++;
908 A_row row(*
i, static_cast<int>(nslice));
910 calcul(
i, row, static_cast<int>(nslice),
tp);
920 const std::string nam =
i->getElement()->getName();
921 if(nam ==
"[DRIFT]")
i++;
925 if(Typ ==
"MARKER")
i++;
926 else if(elem.
getApert().second.size() == 0)
i++;
928 A_row row(*
i, static_cast<int>(nslice));
930 calcul(
i, row, static_cast<int>(nslice),
tp);
950 const ColDesc *col = findCol(*
this, colName);
951 return (this->*(col->get))(row, col->ind_1, col->ind_2);
958 for(
const ColDesc *col = defaultColumns; col->colName; ++col) {
960 new Column(*
this, col->colName, *col);
961 columns.push_back(
Cell(expr, col->printWidth, col->printPrecision));
967 const ColDesc *col = findCol(*
this, colName);
970 std::vector<double> column;
972 for(A_Tline::const_iterator row =
begin(); row !=
end(); ++row) {
975 column.push_back((this->*(col->get))(*row, col->ind_1, col->ind_2));
983 std::vector<std::string> &cols) {
985 std::vector<double> result;
989 for(
const ColDesc *col = defaultColumns; col->colName != 0; ++col) {
990 result.push_back((this->*col->get)(row, col->ind_1, col->ind_2));
994 for(std::vector<std::string>::const_iterator iter = cols.begin();
995 iter != cols.end(); ++iter) {
996 const ColDesc *col = findCol(*
this, *iter);
997 result.push_back((this->*(col->get))(row, col->ind_1, col->ind_2));
1005 if(
itsLine == name)
return true;
1008 for(A_Tline::const_iterator row =
begin(); row !=
end(); ++row) {
1009 if(row->getElement()->getName() ==
name)
return true;
1017 const ColDesc *col = findCol(*
this, colname);
1018 return new Column(*
this, colname, *col);
1036 std::ostringstream os;
1037 os << row << std::ends;
1038 throw OpalException(
"Aperture::findRow()",
"A_row \"" + os.str() +
1039 "\" not found in twiss table \"" +
getOpalName() +
"\".");
virtual std::vector< double > getRow(const PlaceRep &pos, const std::vector< std::string > &cols)
Return a table row.
virtual double getExitFaceRotation() const =0
Get exit pole face rotation.
std::vector< Cell > CellArray
An array of cell descriptors.
void calcul_Apert(A_row &a, int slice, Twiss *tp)
virtual BMultipoleField & getField() override=0
Get multipole expansion of field.
virtual void fill()=0
Refill the buffer.
double getBeta_y(int ind)
double getDisp_y_prim(int ind)
virtual Expressions::PtrToScalar< double > makeColumnExpression(const std::string &colname) const
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
PETE_TUTree< FnArcCos, typename T::PETE_Expr_t > acos(const PETE_Expr< T > &l)
constexpr double e
The value of .
Interface for basic beam line object.
double getDisp_x_prim(int ind)
TLine::const_iterator end() const
Access to last row.
double normal(int) const
Get component.
virtual BMultipoleField & getField() override=0
Get multipole expansion of field.
A templated representation for matrices.
double getEY() const
Return emittance for mode 2.
virtual void visitSBend(const SBend &)
Apply the algorithm to a sector bend.
virtual BMultipoleField & getField() override=0
Get multipole field.
A_row & findRow(const PlaceRep &row)
T det(const AntiSymTenzor< T, 3 > &t)
Interface for a Cyclotron.
TLine::const_iterator begin() const
Access to first row.
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
void leave(const FlaggedElmPtr &) const
Leave an element or line.
A simple arc in the XZ plane.
The base class for all OPAL exceptions.
std::string getType_elm()
Tps< T > sin(const Tps< T > &x)
Sine.
const A_row & getCurrent() const
void execute()
Apply the algorithm to the top-level beamline.
virtual bool isDependent(const std::string &name) const
Find out if table depends on the object identified by [b]name[/b].
virtual double getElementLength() const
Get design length.
void enter(const FlaggedElmPtr &) const
Enter an element or line.
virtual PlanarArcGeometry & getGeometry() override=0
Get SBend geometry.
virtual const std::string & getName() const
Get element name.
Tps< T > tan(const Tps< T > &x)
Tangent.
virtual double getEntryFaceRotation() const =0
Get pole entry face rotation.
std::vector< Attribute > itsAttr
The object attributes (see Attribute.hh).
std::pair< ElementBase::ApertureType, std::vector< double > > getApert() const
bool isActive() const
Return status.
bool dynamic
Flag dynamic table.
Interface for general multipole.
std::string date() const
Return date.
bool getBool(const Attribute &attr)
Return logical value.
A_Tline::const_iterator begin() const
bool isActive() const
Test for active range.
Tps< T > cosh(const Tps< T > &x)
Hyperbolic cosine.
virtual double getExitFaceRotation() const =0
Get exit pole face rotation.
double getBendAngle() const
Get angle.
const Object * getBaseObject() const
Return the object's base type object.
Representation of a place within a beam line or sequence.
double getBETXMAX(const A_row &, int i1=0, int i2=0) const
virtual double getLength()
Return the length of the table.
const std::string & getOpalName() const
Return object name.
void initialize()
Initialise data for search.
FVector< double, 3 > Euler_x
A_Tline::const_iterator current
virtual std::vector< double > getColumn(const RangeRep &rng, const std::string &colName)
Return column [b]col[/b] of this table, limited by [b]range[/b].
constexpr double pi
The value of .
virtual void visitCyclotron(const Cyclotron &)
Apply the algorithm to an cyclotron.
virtual const Beamline * getLine() const
Return embedded CLASSIC beamline.
double skew(int) const
Get component.
struct Aperture::Data data
virtual double getElementLength() const
Get design length.
FMatrix< double, 6, 6 > Transf_mat
void enter(const FlaggedElmPtr &) const
Enter an element or line.
constexpr double c
The velocity of light in m/s.
std::vector< Aperture::coord > getShape(std::vector< double > vec)
static Beam * find(const std::string &name)
Find named BEAM.
A_Tline::const_iterator end() const
void registerOwnership(const AttributeHandler::OwnerType &itsClass) const
Base class for all beam line elements.
static Table * find(const std::string &name)
Find named Table.
std::vector< double > getRealArray(const Attribute &attr)
Get array value.
PETE_TBTree< FnArcTan2, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > atan2(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
double getDisp_x(int ind)
double getBETi(const Row &, int i1, int=0) const
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
std::vector< coord > ShapeHalo
virtual void fill()
Refill the buffer.
void initialize()
Initialise data for search.
Representation of a range within a beam line or sequence.
void leave(const FlaggedElmPtr &) const
Leave an element or line.
virtual double getElementLength() const
Get element length.
An abstract sequence of beam line components.
virtual double getElementLength() const
Get element length.
FVector< double, 2 > Dispx
double getDisp(const Row &, int i1, int=0) const
Dispersion.
Tps< T > sqrt(const Tps< T > &x)
Square root.
virtual double getEntryFaceRotation() const =0
Get pole entry face rotation.
The geometry for a RBend element.
virtual void visitMultipole(const Multipole &)
Apply the algorithm to a multipole.
double getP() const
The constant reference momentum per particle.
static Element * find(const std::string &name)
Find named Element.
double getDisp_y(int ind)
std::vector< pt_interpol > Interpol
virtual void applyDefault(const ElementBase &)
The magnetic field of a multipole.
std::string time() const
Return time.
double getEX() const
Return emittance for mode 1.
const PartData & getReference() const
Return the embedded CLASSIC PartData.
virtual Object * clone(const std::string &name)
Return a clone.
double getM() const
The constant mass per particle.
The base class for all OPAL objects.
Tps< T > cos(const Tps< T > &x)
Cosine.
Attribute makeBool(const std::string &name, const std::string &help)
Make logical attribute.
Descriptor for printing a table cell.
double getBETYMAX(const A_row &, int i1=0, int i2=0) const
virtual void append(const T &)
Append a T object.
virtual double getCell(const PlaceRep &place, const std::string &colName)
Return value in selected table cell.
Attribute makeRealArray(const std::string &name, const std::string &help)
Create real array attribute.
PETE_TUTree< FnArcSin, typename T::PETE_Expr_t > asin(const PETE_Expr< T > &l)
std::vector< coord > ShapeBeamScreen
virtual RBendGeometry & getGeometry() override=0
Get RBend geometry.
virtual double getBendAngle() const
Get angle.
double getBeta_x(int ind)
std::string::iterator iterator
double getAPERTMIN(const A_row &row, int i1=0, int i2=0) const
FVector< double, 2 > Dispy
virtual CellArray getDefault() const
Return the default print columns.
double getReal(const Attribute &attr)
Return real value.
virtual bool matches(Table *rhs) const
Check that [b]rhs[/b] is of same type as [b]this[/b].
The base class for all OPAL tables.
FVector< double, 3 > Euler_y
Attribute makeString(const std::string &name, const std::string &help)
Make string attribute.
ElementBase * getElement() const
Get the element pointer.
virtual void printTable(std::ostream &, const CellArray &) const
Print list for the table.
double getALFi(const Row &, int i1, int=0) const
Attribute makeReal(const std::string &name, const std::string &help)
Make real attribute.
void calcul(Twiss::TLine::iterator i, A_row &a, int order, Twiss *tp)
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
double getS(const Row &, int=0, int=0) const
Arc length for given row.
Tps< T > sinh(const Tps< T > &x)
Hyperbolic sine.
virtual void visitRBend(const RBend &)
Apply the algorithm to a rectangular bend.
Inform & endl(Inform &inf)
A section of a beam line.
std::string getString(const Attribute &attr)
Get string value.