OPAL (Object Oriented Parallel Accelerator Library) 2024.2
OPAL
Ellipse.cpp
Go to the documentation of this file.
4#include "Physics/Physics.h"
5#include "Physics/Units.h"
6
7#include <iostream>
8#include <fstream>
9#include <iomanip>
10#include <vector>
11#include <memory>
12#include <cmath>
13#include <string>
14
15namespace mslang {
16 void Ellipse::print(int indentwidth) {
17 std::string indent(indentwidth, ' ');
18 std::string indent2(indentwidth + 8, ' ');
19 Vector_t origin = trafo_m.getOrigin();
20 double angle = trafo_m.getAngle() * Units::rad2deg;
21 std::cout << indent << "ellipse, \n"
22 << indent2 << "w: " << width_m << ", \n"
23 << indent2 << "h: " << height_m << ", \n"
24 << indent2 << "origin: " << origin[0] << ", " << origin[1] << ",\n"
25 << indent2 << "angle: " << angle << "\n"
26 << indent2 << std::setw(14) << trafo_m(0, 0) << std::setw(14) << trafo_m(0, 1) << std::setw(14) << trafo_m(0, 2) << "\n"
27 << indent2 << std::setw(14) << trafo_m(1, 0) << std::setw(14) << trafo_m(1, 1) << std::setw(14) << trafo_m(1, 2) << "\n"
28 << indent2 << std::setw(14) << trafo_m(2, 0) << std::setw(14) << trafo_m(2, 1) << std::setw(14) << trafo_m(2, 2)
29 << std::endl;
30 }
31
32 void Ellipse::writeGnuplot(std::ofstream &out) const {
33 const unsigned int N = 101;
34 const double dp = Physics::two_pi / (N - 1);
35 const unsigned int colwidth = out.precision() + 8;
36
37 double phi = 0;
38 for (unsigned int i = 0; i < N; ++ i, phi += dp) {
39 Vector_t pt(0.0);
40 pt[0] = std::copysign(sqrt(std::pow(height_m * width_m * 0.25, 2) /
41 (std::pow(height_m * 0.5, 2) +
42 std::pow(width_m * 0.5 * tan(phi), 2))),
43 cos(phi));
44 pt[1] = pt[0] * tan(phi);
45 pt = trafo_m.transformFrom(pt);
46
47 out << std::setw(colwidth) << pt[0]
48 << std::setw(colwidth) << pt[1]
49 << std::endl;
50 }
51 out << std::endl;
52
53 for (auto item: divisor_m) {
54 item->writeGnuplot(out);
55 }
56
57 // bb_m.writeGnuplot(out);
58 }
59
60 void Ellipse::apply(std::vector<std::shared_ptr<Base> > &bfuncs) {
61 bfuncs.emplace_back(this->clone());
62 }
63
64 std::shared_ptr<Base> Ellipse::clone() const{
65 std::shared_ptr<Ellipse> elps(new Ellipse);
66 elps->width_m = width_m;
67 elps->height_m = height_m;
68 elps->trafo_m = trafo_m;
69 elps->bb_m = bb_m;
70
71 for (auto item: divisor_m) {
72 elps->divisor_m.emplace_back(item->clone());
73 }
74
75 return std::static_pointer_cast<Base>(elps);
76 }
77
79 Vector_t llc(0.0), urc(0.0);
80 const Vector_t e_x({1.0, 0.0, 0.0}), e_y({0.0, 1.0, 0.0});
81 const Vector_t center = trafo_m.transformFrom(Vector_t(0.0));
82 const Vector_t e_xp = trafo_m.transformFrom(e_x) - center;
83 const Vector_t e_yp = trafo_m.transformFrom(e_y) - center;
84 const double &M11 = e_xp[0];
85 const double &M12 = e_yp[0];
86 const double &M21 = e_xp[1];
87 const double &M22 = e_yp[1];
88
89 double t = atan2(height_m * M12, width_m * M11);
90 double halfwidth = 0.5 * (M11 * width_m * cos(t) +
91 M12 * height_m * sin(t));
92 llc[0] = center[0] - std::abs(halfwidth);
93 urc[0] = center[0] + std::abs(halfwidth);
94
95 t = atan2(height_m * M22, width_m * M21);
96
97 double halfheight = 0.5 * (M21 * width_m * cos(t) +
98 M22 * height_m * sin(t));
99
100 llc[1] = center[1] - std::abs(halfheight);
101 urc[1] = center[1] + std::abs(halfheight);
102
103 bb_m = BoundingBox2D(llc, urc);
104
105 for (auto item: divisor_m) {
106 item->computeBoundingBox();
107 }
108 }
109
110 bool Ellipse::isInside(const Vector_t &R) const {
111 if (!bb_m.isInside(R)) return false;
112
114 if (4 * (std::pow(X[0] / width_m, 2) + std::pow(X[1] / height_m, 2)) <= 1) {
115
116 for (auto item: divisor_m) {
117 if (item->isInside(R)) return false;
118 }
119
120 return true;
121 }
122
123 return false;
124 }
125
127
128 ArgumentExtractor arguments(std::string(it, end));
129 double width, height;
130 try {
131 width = parseMathExpression(arguments.get(0));
132 height = parseMathExpression(arguments.get(1));
133 } catch (std::runtime_error &e) {
134 std::cout << e.what() << std::endl;
135 return false;
136 }
137
138 Ellipse *elps = static_cast<Ellipse*>(fun);
139 elps->width_m = width;
140 elps->height_m = height;
141
142 if (elps->width_m < 0.0) {
143 std::cout << "Ellipse: a negative width provided '"
144 << arguments.get(0) << " = " << elps->width_m << "'"
145 << std::endl;
146 return false;
147 }
148 if (elps->height_m < 0.0) {
149 std::cout << "Ellipse: a negative height provided '"
150 << arguments.get(1) << " = " << elps->height_m << "'"
151 << std::endl;
152 return false;
153 }
154
155 it += (arguments.getLengthConsumed() + 1);
156
157 return true;
158 }
159}
PartBunchBase< T, Dim >::ConstIterator end(PartBunchBase< T, Dim > const &bunch)
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition TpsMath.h:129
Tps< T > tan(const Tps< T > &x)
Tangent.
Definition TpsMath.h:147
Tps< T > sin(const Tps< T > &x)
Sine.
Definition TpsMath.h:111
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition TpsMath.h:91
PETE_TBTree< FnArcTan2, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > atan2(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
#define X(arg)
Definition fftpack.cpp:106
constexpr double two_pi
The value of.
Definition Physics.h:33
constexpr double rad2deg
Definition Units.h:146
double parseMathExpression(const std::string &str)
Definition matheval.cpp:4
std::string::iterator iterator
Definition MSLang.h:13
std::vector< std::shared_ptr< Base > > divisor_m
Definition MSLang.h:40
AffineTransformation trafo_m
Definition MSLang.h:38
BoundingBox2D bb_m
Definition MSLang.h:39
Vector_t transformTo(const Vector_t &v) const
Vector_t transformFrom(const Vector_t &v) const
std::string get(unsigned int i) const
unsigned int getLengthConsumed() const
bool isInside(const Vector_t &X) const
virtual void computeBoundingBox()
Definition Ellipse.cpp:78
virtual void print(int indentwidth)
Definition Ellipse.cpp:16
virtual bool isInside(const Vector_t &R) const
Definition Ellipse.cpp:110
double width_m
Definition Ellipse.h:8
virtual std::shared_ptr< Base > clone() const
Definition Ellipse.cpp:64
virtual void writeGnuplot(std::ofstream &out) const
Definition Ellipse.cpp:32
virtual void apply(std::vector< std::shared_ptr< Base > > &bfuncs)
Definition Ellipse.cpp:60
double height_m
Definition Ellipse.h:9
static bool parse_detail(iterator &it, const iterator &end, Function *fun)
Definition Ellipse.cpp:126
Vektor< double, 3 > Vector_t
Definition Vektor.h:6