diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp index 07a4c0eefd7452039779f084ff2910f2fb622ff6..ad28bc4871b6b03495a065c44787f8ce7543da5e 100755 --- a/corsika/modules/radio/CoREAS.hpp +++ b/corsika/modules/radio/CoREAS.hpp @@ -92,86 +92,206 @@ namespace corsika { // This is a SignalPathCollection auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)}; - // loop over both paths at once and directly compare 'start' and 'end' attributes - for (size_t i = (paths1.size() == paths2.size()) ? 0 : paths1.size(); // TODO: throw an exception if the sizes don't match - (i < paths1.size() && i < paths2.size()); i++) { - - // First start with the 'start' point - // calculate preDoppler factor - double preDoppler_{1. - paths1[i].refractive_index_source_ * // TODO: use the refractive index at source, not average! - beta_.dot(paths1[i].emit_)}; - std::cout << "preDoppler: " << preDoppler_<< std::endl; - - // calculate receive times for startpoint - auto startPointReceiveTime_{paths1[i].propagation_time_ + startTime_}; // TODO: total time -> propagation time - std::cout << "START RECEIVE TIME: " << startPointReceiveTime_ << std::endl; // TODO: time 0 is when the imaginary primary hits the ground - - // get receive unit vector at 'start' - auto ReceiveVectorStart_ {paths1[i].receive_}; - std::cout << "BETA DOT Path Emit PRE: " << beta_.dot(paths1[i].emit_) << std::endl; - - // calculate electric field vector for startpoint - ElectricFieldVector EV1_ = - paths1[i].receive_.cross(paths1[i].receive_.cross(beta_)) - .getComponents() / - (paths1[i].R_distance_ * preDoppler_) * - (1 / (4 * M_PI * track.getDuration())) * // TODO: divide by sample width not track.getDuration! - ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_; - - // Now continue with the 'end' point - // calculate postDoppler factor - double postDoppler_{ - 1. - paths2[i].refractive_index_source_ * - beta_.dot(paths2[i].emit_)}; // maybe this is path.receive_ (?) - std::cout << "postDoppler: " << postDoppler_<< std::endl; - - // calculate receive times for endpoint - auto endPointReceiveTime_{paths2[i].propagation_time_ + endTime_}; - std::cout << "END RECEIVE TIME: " << endPointReceiveTime_ << std::endl; - - // get receive unit vector at 'end' - auto ReceiveVectorEnd_ {paths2[i].receive_}; - std::cout << "BETA DOT Path Emit POST: " << beta_.dot(paths2[i].emit_) << std::endl; - - // calculate electric field vector for endpoint - ElectricFieldVector EV2_ = - paths2[i].receive_.cross(paths2[i].receive_.cross(beta_)) - .getComponents() / - (paths2[i].R_distance_ * postDoppler_) * - ((-1) / (4 * M_PI * track.getDuration())) * // TODO: divide by sample width not track.getDuration! - ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_; - - ////////////////////////////////////////////////////////////////////////////// - // start comparing stuff - if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) { - - std::cout << "Gets into if less than 1.e-9" << std::endl; - - auto gridResolution_ {antenna.duration_}; - auto deltaT_ { endPointReceiveTime_ - startPointReceiveTime_ }; - - if (std::fabs(deltaT_ / 1_s) < gridResolution_ / 1_s) { - - EV1_ = EV1_ * std::fabs(deltaT_ / gridResolution_); - EV2_ = EV2_ * std::fabs(deltaT_ / gridResolution_); - - const long startBin = static_cast<long>(std::floor(startPointReceiveTime_/gridResolution_+0.5l)); - const long endBin = static_cast<long>(std::floor(endPointReceiveTime_/gridResolution_+0.5l)); - const double startBinFraction = (startPointReceiveTime_/gridResolution_)-std::floor(startPointReceiveTime_/gridResolution_); - const double endBinFraction = (endPointReceiveTime_/gridResolution_)-std::floor(endPointReceiveTime_/gridResolution_); + // loop over both paths at once and directly compare 'start' and 'end' attributes + for (size_t i = (paths1.size() == paths2.size()) ? 0 : paths1.size(); // TODO: throw an exception if the sizes don't match + (i < paths1.size() && i < paths2.size()); i++) { + + // First start with the 'start' point + // calculate preDoppler factor + double preDoppler_{1. - paths1[i].refractive_index_source_ * // TODO: use the refractive index at source, not average! + beta_.dot(paths1[i].emit_)}; + std::cout << "preDoppler: " << preDoppler_<< std::endl; + + // calculate receive times for startpoint + auto startPointReceiveTime_{paths1[i].propagation_time_ + startTime_}; // TODO: total time -> propagation time + std::cout << "START RECEIVE TIME: " << startPointReceiveTime_ << std::endl; // TODO: time 0 is when the imaginary primary hits the ground + + // get receive unit vector at 'start' + auto ReceiveVectorStart_ {paths1[i].receive_}; + std::cout << "BETA DOT Path Emit PRE: " << beta_.dot(paths1[i].emit_) << std::endl; + + // calculate electric field vector for startpoint + ElectricFieldVector EV1_ = + paths1[i].receive_.cross(paths1[i].receive_.cross(beta_)) + .getComponents() / + (paths1[i].R_distance_ * preDoppler_) * + (1 / (4 * M_PI * track.getDuration())) * // TODO: divide by sample width not track.getDuration! + ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_; + + // Now continue with the 'end' point + // calculate postDoppler factor + double postDoppler_{ + 1. - paths2[i].refractive_index_source_ * + beta_.dot(paths2[i].emit_)}; // maybe this is path.receive_ (?) + std::cout << "postDoppler: " << postDoppler_<< std::endl; + + // calculate receive times for endpoint + auto endPointReceiveTime_{paths2[i].propagation_time_ + endTime_}; + std::cout << "END RECEIVE TIME: " << endPointReceiveTime_ << std::endl; + + // get receive unit vector at 'end' + auto ReceiveVectorEnd_ {paths2[i].receive_}; + std::cout << "BETA DOT Path Emit POST: " << beta_.dot(paths2[i].emit_) << std::endl; + + // calculate electric field vector for endpoint + ElectricFieldVector EV2_ = + paths2[i].receive_.cross(paths2[i].receive_.cross(beta_)) + .getComponents() / + (paths2[i].R_distance_ * postDoppler_) * + ((-1) / (4 * M_PI * track.getDuration())) * // TODO: divide by sample width not track.getDuration! + ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_; + + ////////////////////////////////////////////////////////////////////////////// + // start comparing stuff + if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) { + + std::cout << "Gets into if less than 1.e-9" << std::endl; + + auto gridResolution_ {antenna.duration_}; + auto deltaT_ { endPointReceiveTime_ - startPointReceiveTime_ }; + + if (std::fabs(deltaT_ / 1_s) < gridResolution_ / 1_s) { + + EV1_ = EV1_ * std::fabs(deltaT_ / gridResolution_); + EV2_ = EV2_ * std::fabs(deltaT_ / gridResolution_); + + const long startBin = static_cast<long>(std::floor(startPointReceiveTime_/gridResolution_+0.5l)); + const long endBin = static_cast<long>(std::floor(endPointReceiveTime_/gridResolution_+0.5l)); + const double startBinFraction = (startPointReceiveTime_/gridResolution_)-std::floor(startPointReceiveTime_/gridResolution_); + const double endBinFraction = (endPointReceiveTime_/gridResolution_)-std::floor(endPointReceiveTime_/gridResolution_); + + // only do timing modification if contributions would land in same bin + if (startBin == endBin) { + + // if startE arrives before endE + if (deltaT_ / 1_s >= 0) { + if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center + { + startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_; // shift EV1_ to previous gridpoint + } + else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center + { + endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_; // shift EV2_ to next gridpoint + } + else // points on both sides of bin center + { + const double leftDist = 1.0-startBinFraction; + const double rightDist = endBinFraction; + // check if asymmetry to right or left + if (rightDist >= leftDist) + { + endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_; // shift EV2_ to next gridpoint + } + else + { + startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_; // shift EV1_ to previous gridpoint + } + } + } + else // if endE arrives before startE + { + if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center + { + endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_; // shift EV2_ to previous gridpoint + } + else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center + { + startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_; // shift EV1_ to next gridpoint + } + else // points on both sides of bin center + { + const double leftDist = 1.0-endBinFraction; + const double rightDist = startBinFraction; + // check if asymmetry to right or left + if (rightDist >= leftDist) + { + startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_; // shift EV1_ to next gridpoint + } + else + { + endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_; // shift EV2_ to previous gridpoint + } + } + } // End of else statement + } // End of if for startbin == endbin + } // End of if deltaT < gridresolution + } // End of if that checks small doppler factors + + // TODO; fix this if with the one above, they should work together + // perform ZHS-like calculation close to Cherenkov angle + if (std::fabs(preDoppler_) <= approxThreshold_ || std::fabs(postDoppler_) <= approxThreshold_) { + + // clear the existing paths for this particle and track + paths1.clear(); + paths2.clear(); + + // get global simulation time for the middle point of that track. (This is my best guess for now) + auto midTime_{particle.getTime() - (track.getDuration() / 2)}; + + // get "mid" position of the track (that may not work properly) + auto midPoint_{track.getPosition(0.5)}; // TODO: get mid position geometrically + + // get the Path (path3) from the middle "endpoint" to the antenna. + // This is a SignalPathCollection + auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)}; + + // now loop over the paths for endpoint that we got above + for (auto const& path : paths3) { + + auto const midPointReceiveTime_{path.propagation_time_ + midTime_}; + auto midDoppler_{1. - path.average_refractive_index_ * beta_.dot(path.emit_)}; + + // change the values of the receive unit vectors of start and end + ReceiveVectorStart_ = path.receive_; + ReceiveVectorEnd_ = path.receive_; + + // CoREAS calculation -> get ElectricFieldVector3 for "midPoint" + ElectricFieldVector EVmid_ = + path.receive_.cross(path.receive_.cross(beta_)).getComponents() / + (path.R_distance_ * midDoppler_) * (1 / (4 * M_PI * track.getDuration())) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_; + + EV1_ = EVmid_; + EV2_= - EVmid_; + + auto deltaT_{(endPoint_ - startPoint_).getNorm() / (constants::c * beta_.getNorm() * std::fabs(midDoppler_))}; // TODO: Caution with this! + + if (startPointReceiveTime_ < endPointReceiveTime_) // EVstart_ arrives earlier + { + startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_; + endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_; + } + else // EVend_ arrives earlier + { + startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_; + endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_; + } + + const long double gridResolution_{antenna.duration_ / 1_s}; + deltaT_ = endPointReceiveTime_ - startPointReceiveTime_; + + // redistribute contributions over time scale defined by the observation time resolution + if (std::fabs(deltaT_ / 1_s) < gridResolution_) { + + EV1_ = EV1_ * std::fabs((deltaT_ / 1_s) / gridResolution_); + EV2_ = EV2_ * std::fabs((deltaT_ / 1_s) / gridResolution_); + + const long startBin = static_cast<long>(std::floor((startPointReceiveTime_ / 1_s)/gridResolution_+0.5l)); + const long endBin = static_cast<long>(std::floor((endPointReceiveTime_ / 1_s) /gridResolution_+0.5l)); + const double startBinFraction = ((startPointReceiveTime_ / 1_s)/gridResolution_)-std::floor((startPointReceiveTime_ / 1_s)/gridResolution_); + const double endBinFraction = ((endPointReceiveTime_ / 1_s)/gridResolution_)-std::floor((endPointReceiveTime_ / 1_s)/gridResolution_); // only do timing modification if contributions would land in same bin if (startBin == endBin) { // if startE arrives before endE - if (deltaT_ / 1_s >= 0) { + if ((deltaT_ / 1_s) >= 0) { if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center { - startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_; // shift EV1_ to previous gridpoint + startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_ * 1_s; // shift EV1_ to previous gridpoint } else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center { - endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_; // shift EV2_ to next gridpoint + endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_ * 1_s; // shift EV2_ to next gridpoint } else // points on both sides of bin center { @@ -180,11 +300,11 @@ namespace corsika { // check if asymmetry to right or left if (rightDist >= leftDist) { - endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_; // shift EV2_ to next gridpoint + endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_ * 1_s; // shift EV2_ to next gridpoint } else { - startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_; // shift EV1_ to previous gridpoint + startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_ * 1_s; // shift EV1_ to previous gridpoint } } } @@ -192,11 +312,11 @@ namespace corsika { { if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center { - endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_; // shift EV2_ to previous gridpoint + endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_ * 1_s; // shift EV2_ to previous gridpoint } else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center { - startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_; // shift EV1_ to next gridpoint + startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_ * 1_s; // shift EV1_ to next gridpoint } else // points on both sides of bin center { @@ -205,154 +325,34 @@ namespace corsika { // check if asymmetry to right or left if (rightDist >= leftDist) { - startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_; // shift EV1_ to next gridpoint + startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_ * 1_s; // shift EV1_ to next gridpoint } else { - endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_; // shift EV2_ to previous gridpoint + endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_ * 1_s; // shift EV2_ to previous gridpoint } } } // End of else statement } // End of if for startbin == endbin } // End of if deltaT < gridresolution - } // End of if that checks small doppler factors - - // TODO; fix this if with the one above, they should work together - // perform ZHS-like calculation close to Cherenkov angle - if (std::fabs(preDoppler_) <= approxThreshold_ || std::fabs(postDoppler_) <= approxThreshold_) { - - // clear the existing paths for this particle and track - paths1.clear(); - paths2.clear(); - - // get global simulation time for the middle point of that track. (This is my best guess for now) - auto midTime_{particle.getTime() - (track.getDuration() / 2)}; - - // get "mid" position of the track (that may not work properly) - auto midPoint_{track.getPosition(0.5)}; // TODO: get mid position geometrically - - // get the Path (path3) from the middle "endpoint" to the antenna. - // This is a SignalPathCollection - auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)}; - - // now loop over the paths for endpoint that we got above - for (auto const& path : paths3) { - - auto const midPointReceiveTime_{path.propagation_time_ + midTime_}; - auto midDoppler_{1. - path.average_refractive_index_ * beta_.dot(path.emit_)}; - - // change the values of the receive unit vectors of start and end - ReceiveVectorStart_ = path.receive_; - ReceiveVectorEnd_ = path.receive_; - - // CoREAS calculation -> get ElectricFieldVector3 for "midPoint" - ElectricFieldVector EVmid_ = - path.receive_.cross(path.receive_.cross(beta_)).getComponents() / - (path.R_distance_ * midDoppler_) * (1 / (4 * M_PI * track.getDuration())) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_; - EV1_ = EVmid_; - EV2_= - EVmid_; + } // End of looping over paths3 - auto deltaT_{(endPoint_ - startPoint_).getNorm() / (constants::c * beta_.getNorm() * std::fabs(midDoppler_))}; // TODO: Caution with this! - - if (startPointReceiveTime_ < endPointReceiveTime_) // EVstart_ arrives earlier - { - startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_; - endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_; - } - else // EVend_ arrives earlier - { - startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_; - endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_; - } - - const long double gridResolution_{antenna.duration_ / 1_s}; - deltaT_ = endPointReceiveTime_ - startPointReceiveTime_; - - // redistribute contributions over time scale defined by the observation time resolution - if (std::fabs(deltaT_ / 1_s) < gridResolution_) { - - EV1_ = EV1_ * std::fabs((deltaT_ / 1_s) / gridResolution_); - EV2_ = EV2_ * std::fabs((deltaT_ / 1_s) / gridResolution_); - - const long startBin = static_cast<long>(std::floor((startPointReceiveTime_ / 1_s)/gridResolution_+0.5l)); - const long endBin = static_cast<long>(std::floor((endPointReceiveTime_ / 1_s) /gridResolution_+0.5l)); - const double startBinFraction = ((startPointReceiveTime_ / 1_s)/gridResolution_)-std::floor((startPointReceiveTime_ / 1_s)/gridResolution_); - const double endBinFraction = ((endPointReceiveTime_ / 1_s)/gridResolution_)-std::floor((endPointReceiveTime_ / 1_s)/gridResolution_); - - // only do timing modification if contributions would land in same bin - if (startBin == endBin) { - - // if startE arrives before endE - if ((deltaT_ / 1_s) >= 0) { - if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center - { - startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_ * 1_s; // shift EV1_ to previous gridpoint - } - else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center - { - endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_ * 1_s; // shift EV2_ to next gridpoint - } - else // points on both sides of bin center - { - const double leftDist = 1.0-startBinFraction; - const double rightDist = endBinFraction; - // check if asymmetry to right or left - if (rightDist >= leftDist) - { - endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_ * 1_s; // shift EV2_ to next gridpoint - } - else - { - startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_ * 1_s; // shift EV1_ to previous gridpoint - } - } - } - else // if endE arrives before startE - { - if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center - { - endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_ * 1_s; // shift EV2_ to previous gridpoint - } - else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center - { - startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_ * 1_s; // shift EV1_ to next gridpoint - } - else // points on both sides of bin center - { - const double leftDist = 1.0-endBinFraction; - const double rightDist = startBinFraction; - // check if asymmetry to right or left - if (rightDist >= leftDist) - { - startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_ * 1_s; // shift EV1_ to next gridpoint - } - else - { - endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_ * 1_s; // shift EV2_ to previous gridpoint - } - } - } // End of else statement - } // End of if for startbin == endbin - } // End of if deltaT < gridresolution - - } // End of looping over paths3 - - std::cout << "RECEIVE using ZHS-like approximation" << std::endl; - antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_); - antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_); + std::cout << "RECEIVE using ZHS-like approximation" << std::endl; + antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_); + antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_); - } // end of ZHS-like approximation + } // end of ZHS-like approximation - std::cout << "RIGHT BEFORE RECEIVE INCIDENT :" << std::endl; - std::cout << "startTTimes_.at(index): " << startPointReceiveTime_ << std::endl; + std::cout << "RIGHT BEFORE RECEIVE INCIDENT :" << std::endl; + std::cout << "startTTimes_.at(index): " << startPointReceiveTime_ << std::endl; - antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_); - std::cout << "SECOND RECEIVE ! ! !" << std::endl; - std::cout << "endTTimes_.at(index): " << endPointReceiveTime_ << std::endl; - antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_); + antenna.receive(startPointReceiveTime_, ReceiveVectorStart_, EV1_); + std::cout << "SECOND RECEIVE ! ! !" << std::endl; + std::cout << "endTTimes_.at(index): " << endPointReceiveTime_ << std::endl; + antenna.receive(endPointReceiveTime_, ReceiveVectorEnd_, EV2_); - } // End of loop over both paths to get signal info + } // End of loop over both paths to get signal info } // End of looping over antennas } // End of simulate method diff --git a/corsika/modules/radio/RadioProcess.hpp b/corsika/modules/radio/RadioProcess.hpp index 1ad0e020a23095c4bfbdecefb0e0a52e13041846..ecc0943e25a5354edbcda5f4008941e31f5d9a81 100644 --- a/corsika/modules/radio/RadioProcess.hpp +++ b/corsika/modules/radio/RadioProcess.hpp @@ -29,7 +29,7 @@ namespace corsika { */ template <typename TRadioDetector, typename TRadioImpl, typename TPropagator> class RadioProcess : public ContinuousProcess< - RadioProcess<TRadioDetector, TRadioImpl, TPropagator>> { + RadioProcess<TRadioDetector, TRadioImpl, TPropagator>> { // using ParticleType = corsika::setup::Stack::particle_type; // using TrackType = corsika::LeapFrogTrajectory; @@ -112,26 +112,26 @@ namespace corsika { * TODO: This is placeholder so we can use text output while * we wait for the true output formatting to be ready. **/ - bool writeOutput() const { - // this for loop still has some issues - int i = 1; - for (auto& antenna : detector_.getAntennas()) { - - auto [t,E] = antenna.getWaveform(); - auto c = xt::hstack(xt::xtuple(t,E)); - std::ofstream out_file ("antenna" + to_string(i) + "_output.csv"); - xt::dump_csv(out_file, c); - out_file.close(); - ++i; - - } - // how this method should work: - // 1. Loop over the antennas in the collection - // 2. Get their waveforms - // 3. Create a text file for each antenna - // 4. and write out two columns, time and field. + bool writeOutput() const { + // this for loop still has some issues + int i = 1; + for (auto& antenna : detector_.getAntennas()) { + + auto [t,E] = antenna.getWaveform(); + auto c = xt::hstack(xt::xtuple(t,E)); + std::ofstream out_file ("antenna" + to_string(i) + "_output.csv"); + xt::dump_csv(out_file, c); + out_file.close(); + ++i; } + // how this method should work: + // 1. Loop over the antennas in the collection + // 2. Get their waveforms + // 3. Create a text file for each antenna + // 4. and write out two columns, time and field. + + } /** * Return the maximum step length for this particle and track. @@ -150,4 +150,4 @@ namespace corsika { }; // END: class RadioProcess -} // namespace corsika +} // namespace corsika \ No newline at end of file diff --git a/corsika/modules/radio/ZHS.hpp b/corsika/modules/radio/ZHS.hpp index 8b7662980bd85960ec75797269861cc093b24f45..1b49c826227363a0362eed151bc59a63fc8fe2ef 100644 --- a/corsika/modules/radio/ZHS.hpp +++ b/corsika/modules/radio/ZHS.hpp @@ -22,9 +22,9 @@ namespace corsika { template <typename TRadioDetector, typename TPropagator> class ZHS final : public RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator> { - using Base = RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator>; - using Base::detector_; - + using Base = RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator>; + using Base::detector_; + public: using ElectricFieldVector = QuantityVector<ElectricFieldType::dimension_type>; @@ -113,7 +113,7 @@ namespace corsika { */ template <typename Particle, typename Track> LengthType MaxStepLength(Particle const& particle, - Track const& track) const { + Track const& track) const { // TODO : This is where we control the maximum step size // of a particle track in order to maintain the accuracy @@ -127,4 +127,4 @@ namespace corsika { }; // END: class RadioProcess -} // namespace corsika +} // namespace corsika \ No newline at end of file diff --git a/corsika/modules/radio/antennas/Antenna.hpp b/corsika/modules/radio/antennas/Antenna.hpp index 9bb58f3f67d4aaabcd89c1122ecd81f491cbd48e..8fe910cc8dba7be74856e0ae67c88b38dcde93b2 100644 --- a/corsika/modules/radio/antennas/Antenna.hpp +++ b/corsika/modules/radio/antennas/Antenna.hpp @@ -31,7 +31,7 @@ namespace corsika { // this stores the polarization vector of an electric field using ElectricFieldVector = - QuantityVector<ElectricFieldType::dimension_type>; + QuantityVector<ElectricFieldType::dimension_type>; // using MagneticFieldVector = // QuantityVector<MagneticFieldType::dimension_type>; @@ -61,8 +61,8 @@ namespace corsika { * for the particular antenna implementation and usage. * */ - template <typename... TVArgs> - void receive(TVArgs&&... args); + template <typename... TVArgs> + void receive(TVArgs&&... args); /** * Get the location of this antenna.