IAP GITLAB

Skip to content
Snippets Groups Projects
Commit b911ca6b authored by Juan Ammerman's avatar Juan Ammerman Committed by Nikos Karastathis
Browse files

Draft implementation of ZHS

parent cfca107c
No related branches found
No related tags found
1 merge request!329Radio interface
Hide whitespace changes
Inline Side-by-side
corsika/modules/radio/ZHS.hpp
+ 125
− 39
View file @ b911ca6b
......@@ -12,7 +12,7 @@
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/modules/radio/propagators/SignalPath.hpp>
#include <bits/stdc++.h>
#include <cmath>
namespace corsika {
......@@ -26,8 +26,8 @@ namespace corsika {
using Base::detector_;
public:
using ElectricFieldVector =
QuantityVector<ElectricFieldType::dimension_type>;
// using PotentialVector = QuantityVector<PotentialVectorType::dimension_type>;
using ElectricFieldVector = QuantityVector<ElectricFieldType::dimension_type>;
/**
* Construct a new ZHS instance.
*
......@@ -51,55 +51,141 @@ namespace corsika {
template <typename Particle, typename Track>
ProcessReturn simulate(Particle& particle, Track const& track) const {
//get global simulation time for that track. (This is my best guess for now)
auto startTime_ {particle.getTime()
// TODO: think if we reuse these variables for the case of not being in the Fraunhoffer approx.
auto const startTime_ {particle.getTime()
- track.getDuration()}; // time at start point of track.
auto endTime_ {particle.getTime()}; // time at end point of track.
auto midTime_ {(startTime_ + endTime_) / 2};
auto const endTime_ {particle.getTime()}; // time at end point of track.
auto startPoint_ = track.getPosition(0);
auto const startPoint_{track.getPosition(0)};
auto const endPoint_{track.getPosition(1)};
LengthType trackLength_ {(startPoint_ - endPoint_).getNorm()};
// track velocity
auto trackVelocity_ {(track.getVelocity(0) + track.getVelocity(1)) / 2};
auto const trackVelocity_ {(track.getVelocity(0) + track.getVelocity(1)) / 2};
// beta is defined as velocity / speed of light
auto beta_ { trackVelocity_ / constants::c};
auto const beta_ { trackVelocity_ / constants::c};
// get particle charge
auto const charge_ {get_charge(particle.getPID())};
// get "mid" position of the track geometrically
auto const midVector_ { (startPoint_ - endPoint_) / 2 };
auto const midPoint_ {Point(midVector_.getCoordinateSystem(), midVector_.getComponents().getX(),
midVector_.getComponents().getY(), midVector_.getComponents().getZ())};
//changed if deltaT1_ > deltaT2_, so not const
auto constants {charge_/ (4 * M_PI)/ (constants::epsilonZero) / constants::c};
// we loop over each antenna in the collection
for (auto& antenna : detector_.getAntennas()) {
// Probably will be reused in the case of not being in Fraunhoffer.
auto midPaths{this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)};
//Loop over midPaths to first check Fraunhoffer limit
for (size_t i{0}; i < midPaths.size(); i++){
double const betaTimesK {beta_.dot(midPaths[i].emit_)/beta_.getNorm()};
double const slitSize_ {1. - betaTimesK * betaTimesK * trackLength_ * trackLength_ /1_m /1_m};
//Parameter that determines the limit for the Fraunhoffer limit (probably related to antenna sampling rate
double const lambda {0.1};
double const fraunhLimit {slitSize_/midPaths[i].R_distance_/lambda * 1_m};
if (fraunhLimit > 1.) // Checks if we are in fraunhoffer domain (maybe it should be less?)
{
///code for dividing track and calculating field.
std::cout << "This code hasnt been implemented!" << std::endl;
}
else // Calculate vector potential of whole track
{
// paths from start of track to antenna
auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_m)};
// paths from end of track to antenna
auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
// First implementation, need to check if there are the same number of paths, but once its working.
//modified later if deltaT1_ > deltaT2_
auto deltaT1_ {startTime_ + paths1[i].propagation_time_};
auto deltaT2_ {endTime_ + paths2[i].propagation_time_};
long double const gridResolution_{antenna.duration_ / 1_s};
//make deltaT1_ be the smallest time => changes step function order so
// constants is changed to account for it
if (deltaT1_ > deltaT2_)
{
deltaT1_ = {endTime_ + paths2[i].propagation_time_};
deltaT2_ = {startTime_ + paths1[i].propagation_time_};
constants = - constants;
}
long const startBin {std::floor((deltaT1_ / 1_s)/gridResolution_)};
long const endBin {std::floor((deltaT2_ / 1_s) /gridResolution_)};
auto const betaPerp_ {midPaths[i].emit_.cross(beta_.cross(midPaths[i].emit_))};
double const denominator {1 - midPaths[i].refractive_index_source_ *
beta_.dot(midPaths[i].emit_)};
long const numberOfBins{endBin - startBin};
//IMPORTANT!!!!! Using ElectricFieldVector instead of PotentialVector until antenna is adapted
for (long int it{0}; it <= numberOfBins; ++it)
{
if (startBin == endBin)//track contained in bin
{
// std::cout << "Track in one bin !" << std::endl;
if (std::fabs(denominator)>1.e-15L)//if not in Cerenkov angle then
{
double const f {std::fabs((deltaT2_ - deltaT1_)/1_s)/gridResolution_};
//should be PotentialVector const Vp_ = betaPerp_.getComponents()/denominator/
// midPaths[i].R_distance_ * constants * f;
// but to make it compile until antenna is adapted it stays like that to test it.
ElectricFieldVector const Vp_ = betaPerp_.getComponents()/denominator/
midPaths[i].R_distance_ * constants * f / 1_s;
// std::cout << "Vector Potential NOT in Cerenkov Angle: " << Vp_ << std::endl;
// std::cout << "Time " << deltaT2_ << std::endl;
antenna.receive(deltaT2_, betaPerp_, Vp_);
}
else//If emission in Cerenkov angle => approximation
{
double const f {(endTime_ - startTime_)/gridResolution_/1_s};
ElectricFieldVector const Vp_ = betaPerp_.getComponents()/
midPaths[i].R_distance_ * constants * f / 1_s;
// std::cout << "Vector Potential in Cerenkov Angle: " << Vp_ << std::endl;
// std::cout << "Time " << deltaT1_ << std::endl;
antenna.receive(deltaT2_, betaPerp_, Vp_);
}//end if Cerenkov angle approx
}
//TODO: should we check for Cerenkov angle?
else //Track is contained in more than one bin
{
// std::cout << "Track in multiple bins !" << std::endl;
if (it == 0)//start of track
{
double const f {std::fabs(startTime_/1_s / gridResolution_ - startBin)};
ElectricFieldVector const Vp_ = betaPerp_.getComponents() * f *
constants/denominator/midPaths[i].R_distance_ / 1_s;
antenna.receive(deltaT1_, betaPerp_, Vp_);
}
else if (it == numberOfBins)//end of track
{
double const f {std::fabs(endTime_/1_s / gridResolution_ - endBin)};
ElectricFieldVector const Vp_ = betaPerp_.getComponents() * f *
constants / denominator / midPaths[i].R_distance_ / 1_s;
antenna.receive(deltaT2_, betaPerp_, Vp_);
}
else
{
ElectricFieldVector const Vp_ = betaPerp_.getComponents() * constants /
denominator/midPaths[i].R_distance_ / 1_s;
antenna.receive(deltaT1_ + gridResolution_ * it * 1_s, betaPerp_, Vp_);
}
}
}//end loop over bins in which potential vector is not zero
}//finish if statement of track in fraunhoffer or not
}//end loop over mid paths
} // END: loop over antennas
return ProcessReturn::Ok;
} //end simulate
// get the Path from the track to the antenna
// This is a SignalPathCollection
auto paths{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_nm)};
// now loop over the paths that we got above
// Note: for the StraightPropagator, there will only be a single
// path but other propagators may return more than one.
for (auto const& path : paths) {
auto t1 = path.total_time_;
QuantityVector<ElectricFieldType::dimension_type> v11{10_V / 1_m, 10_V / 1_m, 10_V / 1_m};
antenna.receive(t1, path.emit_, v11);
// calculate the ZHS formalism for this particle-antenna
// ElectricFieldVector EV_ = ((- constants) * trackVelocity_.dot(path.emit)) *
// ((midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
// acos(track.getDirection(0).dot(path.emit_))) * startTime_)
// - (midTime_ + path.total_time_ - (1 - path.average_refractivity_ * beta_ *
// acos(track.getDirection(0).dot(path.emit_))) * endTime_))
// / (1 - path.average_refractivity_ * beta_ * acos(track.getDirection(0).dot(path.emit_)));
//
// // pass it to the antenna
// antenna.receive(midTime_ + path.total_time_, path.receive_, EV_);
} // END: loop over paths
} // END: loop over antennas
}
/**
* Return the maximum step length for this particle and track.
......
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

data privacy - impressum