IAP GITLAB

Commit 0b3ca6f3 authored by André Schmidt's avatar André Schmidt

Merge branch 'andre-master' into 'master'

Added Intersection

See merge request !7
parents 12437888 974b0364
Pipeline #2015 failed with stages
......@@ -97,7 +97,7 @@ int main(int argc, char** argv) {
node1->SetModelProperties<UniformMagneticField<HomogeneousMedium<IEnvModel>>>(
B, 1_g / (1_cm * 1_cm * 1e9_cm), composition);
auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
/*auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
std::make_unique<geometry::Sphere>(center, seaLevel + 100.0_km));
node2->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
B, center, 540.1778_g / (1_cm * 1_cm) / 772170.16_cm, -772170.16_cm, composition,
......@@ -124,7 +124,7 @@ int main(int argc, char** argv) {
node4->AddChild(std::move(node5));
node3->AddChild(std::move(node4));
node2->AddChild(std::move(node3));
node1->AddChild(std::move(node2));
node1->AddChild(std::move(node2));*/
env.GetUniverse()->AddChild(std::move(node1));
// setup particle stack, and add primary particle
......
......@@ -229,26 +229,24 @@ namespace corsika::cascade {
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
LengthType const steplength = min_distance /
(directionBefore + directionBefore.cross(magneticfield) * k / 2).GetNorm();
// First Movement
//assuming magnetic field does not change during movement
auto position = vParticle.GetPosition() + directionBefore * Steplength / 2;
auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
min_distance * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - vParticle.GetPosition()) /
(position - vParticle.GetPosition()).GetNorm();
vParticle.SetMomentum(direction * vParticle.GetMomentum().GetNorm());
position = position + directionAfter * min_distance / 2;
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
vParticle.SetPosition(position);
} else {
vParticle.SetPosition(step.PositionFromArclength(min_distance));
}
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ...
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
step.LimitEndTo(min_distance);
......
......@@ -61,76 +61,36 @@ namespace corsika::process {
// << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
// determine direction of the particle after adding magnetic field
auto const* currentLogicalVolumeNode = p.GetNode();
int chargeNumber;
if(corsika::particles::IsNucleus(p.GetPID())) {
chargeNumber = p.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
}
if(chargeNumber != 0) {
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
std::cout << "TrackingLine B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
//determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - .GetCenter()) * k + 1) /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 * 1_m * 1_m / 4);
double b = directionBefore.dot(currentPosition - .GetCenter()) * 2 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 * 1_m / 4);
double c = ((currentPosition - .GetCenter()).GetSquaredNorm() - .GetRadius^2) /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k^2 / 4);
std::complex<double>* solutions = solve_quartic(0, a, b, c);
std::vector<double> tmp;
for(int i = 0; i < 4; i++) {
if(solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType const Steplength;
if(tmp.size() > 0) {
Steplength = *std::min_element(tmp.begin(),tmp.end()) * 1_m;
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection with anything!" << std::endl;
//what to do when this happens?
}
// First Movement
//assuming magnetic field does not change during movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity = direction * velocity.GetNorm();
std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
<< " GeV " << std::endl;
} else {
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
}
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
geometry::Line line(currentPosition, velocity);
//auto const* currentLogicalVolumeNode = p.GetNode();
auto const* currentLogicalVolumeNode = p.GetNode();
//~ auto const* currentNumericalVolumeNode =
//~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
auto const numericallyInside =
currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false") << std::endl;
auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
//charge of the particle
int chargeNumber;
if(corsika::particles::IsNucleus(p.GetPID())) {
chargeNumber = p.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
}
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
std::cout << " Magnetic Field: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
geometry::Vector<SpeedType::dimension_type> velocity1 = velocity;
geometry::Vector<SpeedType::dimension_type> velocity2 = velocity;
// for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
......@@ -138,6 +98,47 @@ namespace corsika::process {
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
//creating Line with magnetic field
if(chargeNumber != 0) {
//determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m);
double c = ((currentPosition - sphere.GetCenter()).GetSquaredNorm() -
(sphere.GetRadius() * sphere.GetRadius())) * 4 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
std::complex<double>* solutions = solve_quartic(0, a, b, c);
std::vector<double> tmp;
for(int i = 0; i < 4; i++) {
if(solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType Steplength;
if(tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection (1)!" << std::endl;
//what to do when this happens?
}
// First Movement
//assuming magnetic field does not change during movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity1 = direction * velocity.GetNorm();
} //instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
//line has some errors
geometry::Line line(currentPosition, velocity1);
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
......@@ -160,6 +161,47 @@ namespace corsika::process {
currentLogicalVolumeNode->GetVolume());
// for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
//creating Line with magnetic field
if(chargeNumber != 0) {
//determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
double b = directionBefore.dot(currentPosition - sphere.GetCenter()) * 8 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m);
double c = ((currentPosition - sphere.GetCenter()).GetSquaredNorm() -
(sphere.GetRadius() * sphere.GetRadius())) * 4 /
((directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k * 1_m * 1_m * 1_m * 1_m);
std::complex<double>* solutions = solve_quartic(0, a, b, c);
std::vector<double> tmp;
for(int i = 0; i < 4; i++) {
if(solutions[i].imag() == 0 && solutions[i].real() > 0) {
tmp.push_back(solutions[i].real());
}
}
LengthType Steplength;
if(tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
std::cout << "s = " << Steplength << std::endl;
} else {
std::cout << "no intersection (2)!" << std::endl;
//what to do when this happens?
}
// First Movement
//assuming magnetic field does not change during movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity2 = direction * velocity.GetNorm();
}//instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
geometry::Line line(currentPosition, velocity2);
[[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
[[maybe_unused]] auto dummy_t1 = t1;
intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
......@@ -183,6 +225,20 @@ namespace corsika::process {
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
// determine direction of the particle after adding magnetic field
//assuming magnetic field does not change during movement
// First Movement
auto position = currentPosition + velocity * min / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + velocity.cross(magneticfield) *
min * k;
// Second Movement
position = position + directionAfter * velocity.norm() * min / 2;
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity = direction * velocity.GetNorm();
geometry::Line line(currentPosition, velocity);
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
velocity.norm() * min, minIter->second);
......
......@@ -115,6 +115,11 @@ namespace corsika::process {
std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
<< " GeV " << std::endl;
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
double test =((directionBefore.cross(magneticfield)).dot(position-currentPosition) * k + 1) / (1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
std::cout << "Test: " << test << k << std::endl;
} else {
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
......@@ -143,6 +148,8 @@ namespace corsika::process {
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
std::cout << "Mittelpunkt: " << sphere.GetCenter().GetCoordinates() << std::endl;
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
......
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