Commit bb20f319 authored by kraus's avatar kraus
Browse files

first steps towards collimator with flexible shape

parent a24ff73a
This diff is collapsed.
......@@ -13,7 +13,6 @@
#include <streambuf>
#include <cstdlib>
#include <cmath>
#include <tuple>
std::string UDouble("([0-9]+\\.?[0-9]*([Ee][+-]?[0-9]+)?)");
std::string Double("(-?[0-9]+\\.?[0-9]*([Ee][+-]?[0-9]+)?)");
......@@ -22,11 +21,7 @@ std::string FCall("([a-z]*)\\((.*)");
typedef std::string::iterator iterator;
struct BoundingVolume {
virtual bool isInside(const Vector_t &X) const = 0;
};
struct BoundingBox: public BoundingVolume {
struct BoundingBox {
Vector_t center_m;
double width_m;
double height_m;
......@@ -44,7 +39,7 @@ struct BoundingBox: public BoundingVolume {
height_m(urc[1] - llc[1])
{ }
virtual bool isInside(const Vector_t &X) const {
bool isInside(const Vector_t &X) const {
if (2 * std::abs(X[0] - center_m[0]) <= width_m &&
2 * std::abs(X[1] - center_m[1]) <= height_m)
return true;
......@@ -69,26 +64,6 @@ struct BoundingBox: public BoundingVolume {
}
};
struct BoundingSphere: public BoundingVolume {
Vector_t center_m;
double radiusSqr_m;
BoundingSphere(const Vector_t &c,
double r):
center_m(c),
radiusSqr_m(r*r)
{ }
virtual bool isInside(const Vector_t &X) const {
if (std::pow(X[0] - center_m[0], 2) +
std::pow(X[1] - center_m[1], 2) < radiusSqr_m)
return true;
return false;
}
};
struct AffineTransformation: public Tenzor<double, 3> {
AffineTransformation(const Vector_t& row0,
const Vector_t& row1):
......@@ -161,12 +136,9 @@ bool parse(iterator &it, const iterator &end, Function* &fun);
struct Base: public Function {
AffineTransformation trafo_m;
BoundingBox bb_m;
// std::tuple<unsigned int,
// double,
// double> repeat_m;
Base():
trafo_m()
// , repeat_m(std::make_tuple(0u, 0.0, 0.0))
{ }
virtual Base* clone() = 0;
......@@ -433,7 +405,6 @@ struct Repeat: public Function {
virtual void apply(std::vector<Base*> &bfuncs) {
AffineTransformation shift(Vector_t(1.0, 0.0, -shiftx_m),
Vector_t(0.0, 1.0, -shifty_m));
std::cout << shiftx_m << "\t" << shift(0, 2) << std::endl;
func_m->apply(bfuncs);
const unsigned int size = bfuncs.size();
......@@ -822,13 +793,11 @@ int main(int argc, char *argv[])
std::ofstream out("test.gpl");
for (Base* bfun: baseBlocks) {
bfun->print(0);
bfun->computeBoundingBox();
std::cout << std::endl;
bfun->computeBoundingBox();
bfun->writeGnuplot(out);
}
out.close();
out.open("particles.gpl");
if (baseBlocks.size() > 0) {
Vector_t llc, urc;
Base* first = baseBlocks.front();
......@@ -848,6 +817,55 @@ int main(int argc, char *argv[])
urc[1] = std::max(urc[1], bb.center_m[1] + 0.5 * bb.height_m);
}
struct CoordinateComp {
double delta_m;
CoordinateComp(double delta):
delta_m(delta)
{ }
bool operator()(double a, double b) const {
return a < b && std::abs(a - b) > delta_m;
}
};
std::set<double, CoordinateComp> xCoordinates(CoordinateComp(1e-6 * (urc[0] - llc[0])));
std::set<double, CoordinateComp> yCoordinates(CoordinateComp(1e-6 * (urc[1] - llc[1])));
for (Base *func: baseBlocks) {
Vector_t center = func->trafo_m.transformFrom(Vector_t(0.0));
xCoordinates.insert(center[0]);
yCoordinates.insert(center[1]);
}
std::set<double> dualXCoordinates;
auto it = xCoordinates.begin();
auto next = std::next(it);
dualXCoordinates.insert(std::min(llc[0], 1.5 * (*it) - 0.5 * (*next)));
for (; next != xCoordinates.end(); ++ it, ++ next) {
dualXCoordinates.insert(0.5 * ((*next) + (*it)));
}
dualXCoordinates.insert(std::max(urc[0], 1.5 * (*it) - 0.5 * (*std::prev(it))));
std::set<double> dualYCoordinates;
it = yCoordinates.begin();
next = std::next(it);
dualYCoordinates.insert(std::min(llc[1], 1.5 * (*it) - 0.5 * (*next)));
for (; next != yCoordinates.end(); ++ it, ++ next) {
dualYCoordinates.insert(0.5 * ((*next) + (*it)));
}
dualXCoordinates.insert(std::max(urc[0], 1.5 * (*it) - 0.5 * (*std::prev(it))));
for (auto it = dualXCoordinates.begin(); it != dualXCoordinates.end(); ++ it) {
out << std::setw(14) << *it << std::setw(14) << *(dualYCoordinates.begin()) << "\n"
<< std::setw(14) << *it << std::setw(14) << *(dualYCoordinates.rbegin()) << "\n"
<< std::endl;
}
for (auto it = dualYCoordinates.begin(); it != dualYCoordinates.end(); ++ it) {
out << std::setw(14) << *(dualXCoordinates.begin()) << std::setw(14) << *it << "\n"
<< std::setw(14) << *(dualXCoordinates.rbegin()) << std::setw(14) << *it << "\n"
<< std::endl;
}
out.close();
out.open("particles.gpl");
gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
for (unsigned int i = 0; i < 100000; ++ i) {
......@@ -867,11 +885,13 @@ int main(int argc, char *argv[])
gsl_rng_free(rng);
}
for (Base* func: baseBlocks) {
delete func;
for (Base* func: baseBlocks) {
delete func;
}
}
out.close();
}
delete fun;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment