Skip to content

Commit

Permalink
Add AllSkyVCLFocalPlaneRayPropagator128
Browse files Browse the repository at this point in the history
  • Loading branch information
@sfegan
sfegan committed Oct 12, 2023
1 parent c122794 commit dff89d6
Show file tree
Hide file tree
Showing 3 changed files with 152 additions and 17 deletions.
42 changes: 36 additions & 6 deletions include/simulation/vcl_iact_array.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,8 @@ template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) VCL

CALIN_TYPEALIAS(FocalPlaneRayPropagator, calin::simulation::vcl_ray_propagator::VCLFocalPlaneRayPropagator<VCLArchitecture>);
CALIN_TYPEALIAS(DaviesCottonVCLFocalPlaneRayPropagator, calin::simulation::vcl_ray_propagator::DaviesCottonVCLFocalPlaneRayPropagator<VCLArchitecture>);
CALIN_TYPEALIAS(TrivialVCLFocalPlaneRayPropagator, calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator<VCLArchitecture>);
CALIN_TYPEALIAS(PerfectOpticsVCLFocalPlaneRayPropagator, calin::simulation::vcl_ray_propagator::PerfectOpticsVCLFocalPlaneRayPropagator<VCLArchitecture>);
CALIN_TYPEALIAS(AllSkyVCLFocalPlaneRayPropagator, calin::simulation::vcl_ray_propagator::AllSkyVCLFocalPlaneRayPropagator<VCLArchitecture>);

VCLIACTArray(calin::simulation::atmosphere::LayeredRefractiveAtmosphere* atm,
const calin::simulation::detector_efficiency::AtmosphericAbsorption& atm_abs,
Expand Down Expand Up @@ -289,13 +290,19 @@ template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) VCL
const std::string& propagator_name, SplinePEAmplitudeGenerator* pe_generator = nullptr,
bool adopt_pe_processor = false, bool adopt_pe_generator = false);

TrivialVCLFocalPlaneRayPropagator* add_trivial_propagator(
PerfectOpticsVCLFocalPlaneRayPropagator* add_perfect_optics_propagator(
const Eigen::VectorXd& x, const Eigen::VectorXd& y, const Eigen::VectorXd& z,
double radius, double focal_length, double field_of_view_radius, PEProcessor* pe_processor,
const DetectionEfficiency& detector_efficiency, const AngularEfficiency& fp_angular_efficiency,
const std::string& propagator_name, SplinePEAmplitudeGenerator* pe_generator = nullptr,
bool adopt_pe_processor = false, bool adopt_pe_generator = false);

AllSkyVCLFocalPlaneRayPropagator* add_all_sky_propagator(
Eigen::VectorXd& r0, double radius, PEProcessor* pe_processor,
const DetectionEfficiency& detector_efficiency, const AngularEfficiency& fp_angular_efficiency,
const std::string& propagator_name, SplinePEAmplitudeGenerator* pe_generator = nullptr,
bool adopt_pe_processor = false, bool adopt_pe_generator = false);

void point_telescope_az_el_phi_deg(unsigned iscope, double az_deg, double el_deg, double phi_deg);
void point_telescope_az_el_deg(unsigned iscope, double az_deg, double el_deg);

Expand Down Expand Up @@ -574,8 +581,8 @@ VCLIACTArray<VCLArchitecture>::add_davies_cotton_propagator(
}

template<typename VCLArchitecture>
calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator<VCLArchitecture>*
VCLIACTArray<VCLArchitecture>::add_trivial_propagator(
calin::simulation::vcl_ray_propagator::PerfectOpticsVCLFocalPlaneRayPropagator<VCLArchitecture>*
VCLIACTArray<VCLArchitecture>::add_perfect_optics_propagator(
const Eigen::VectorXd& x, const Eigen::VectorXd& y, const Eigen::VectorXd& z,
double radius, double focal_length, double field_of_view_radius, PEProcessor* pe_processor,
const DetectionEfficiency& detector_efficiency, const AngularEfficiency& fp_angular_efficiency,
Expand All @@ -585,8 +592,8 @@ VCLIACTArray<VCLArchitecture>::add_trivial_propagator(
if(x.size() != y.size() or x.size() != z.size()) {
throw std::runtime_error("Number of telescope x, y and z coordinates must be equal");
}
auto* propagator = new calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator<VCLArchitecture>(
this->rng_, ref_index_, /* adopt_rng= */ false);
auto* propagator = new calin::simulation::vcl_ray_propagator::PerfectOpticsVCLFocalPlaneRayPropagator<VCLArchitecture>(
ref_index_);
for(unsigned iscope=0; iscope<x.size(); ++iscope) {
Eigen::Vector3d r0;
r0 << x[iscope], y[iscope], z[iscope];
Expand All @@ -604,6 +611,29 @@ VCLIACTArray<VCLArchitecture>::add_trivial_propagator(
return propagator;
}

template<typename VCLArchitecture>
calin::simulation::vcl_ray_propagator::AllSkyVCLFocalPlaneRayPropagator<VCLArchitecture>*
VCLIACTArray<VCLArchitecture>::add_all_sky_propagator(
Eigen::VectorXd& r0, double radius, PEProcessor* pe_processor,
const DetectionEfficiency& detector_efficiency, const AngularEfficiency& fp_angular_efficiency,
const std::string& propagator_name, SplinePEAmplitudeGenerator* pe_generator,
bool adopt_pe_processor, bool adopt_pe_generator)
{
auto* propagator = new calin::simulation::vcl_ray_propagator::AllSkyVCLFocalPlaneRayPropagator<VCLArchitecture>(
config_.observation_level(), r0, radius, /* field_of_view_radius = */ M_PI/2, ref_index_);

auto* bandwidth_manager = new VCLDCBandwidthManager<VCLArchitecture>(
&atm_abs_, detector_efficiency, fp_angular_efficiency, zobs_,
config_.detector_energy_lo(), config_.detector_energy_hi(),
config_.detector_energy_bin_width(), propagator_name);

add_propagator(propagator, pe_processor, bandwidth_manager, pe_generator, propagator_name,
/* adopt_propagator = */ true, adopt_pe_processor, adopt_pe_generator);

return propagator;

}

template<typename VCLArchitecture> void VCLIACTArray<VCLArchitecture>::
point_telescope_az_el_phi_deg(unsigned iscope,
double az_deg, double el_deg, double phi_deg)
Expand Down
117 changes: 109 additions & 8 deletions include/simulation/vcl_ray_propagator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,7 @@ template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) Dav
#endif
};

template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) TrivialVCLFocalPlaneRayPropagator:
template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) PerfectOpticsVCLFocalPlaneRayPropagator:
public VCLFocalPlaneRayPropagator<VCLArchitecture>
{
public:
Expand All @@ -214,22 +214,18 @@ template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) Tri
using Real_vt = calin::util::vcl::VCLDoubleReal<VCLArchitecture>;
using Ray_vt = typename calin::math::ray::VCLRay<Real_vt>;

using RayTracer = calin::simulation::vcl_raytracer::VCLScopeRayTracer<Real_vt>;
using TraceInfo = calin::simulation::vcl_raytracer::VCLScopeTraceInfo<Real_vt>;
using RealRNG = calin::math::rng::VCLRealRNG<Real_vt>;
#endif // not defined SWIG

CALIN_TYPEALIAS(ArchRNG, calin::math::rng::VCLRNG<VCLArchitecture>);

TrivialVCLFocalPlaneRayPropagator(ArchRNG* rng = nullptr, double ref_index = 1.0, bool adopt_rng = false):
rng_(new RealRNG(rng==nullptr ? new ArchRNG(__PRETTY_FUNCTION__) : rng,
rng==nullptr ? true : adopt_rng)),
PerfectOpticsVCLFocalPlaneRayPropagator(double ref_index = 1.0):
ref_index_(ref_index)
{
// nothing to see here
}

virtual ~TrivialVCLFocalPlaneRayPropagator() {
virtual ~PerfectOpticsVCLFocalPlaneRayPropagator() {
for(auto* scope : scopes_) {
delete scope;
}
Expand Down Expand Up @@ -346,10 +342,115 @@ template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) Tri
double focal_length; // Focal length [cm]
};

RealRNG* rng_ = nullptr;
double ref_index_ = 1.0;
std::vector<TelescopeDetails*> scopes_;
#endif
};

template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) AllSkyVCLFocalPlaneRayPropagator:
public VCLFocalPlaneRayPropagator<VCLArchitecture>
{
public:
#ifndef SWIG
using int64_vt = typename VCLArchitecture::int64_vt;
using double_bvt = typename VCLArchitecture::double_bvt;
using double_vt = typename VCLArchitecture::double_vt;
using Vector3d_vt = typename VCLArchitecture::Vector3d_vt;
using Real_vt = calin::util::vcl::VCLDoubleReal<VCLArchitecture>;
using Ray_vt = typename calin::math::ray::VCLRay<Real_vt>;

using TraceInfo = calin::simulation::vcl_raytracer::VCLScopeTraceInfo<Real_vt>;
#endif // not defined SWIG

CALIN_TYPEALIAS(ArchRNG, calin::math::rng::VCLRNG<VCLArchitecture>);

AllSkyVCLFocalPlaneRayPropagator(unsigned observation_level,
const Eigen::VectorXd& sphere_center, double sphere_radius, double field_of_view_radius = M_PI/2,
double ref_index = 1.0):
observation_level_(observation_level), sphere_center_(sphere_center),
sphere_radius_(sphere_radius), sphere_radius_squared_(sphere_radius * sphere_radius),
sphere_field_of_view_radius_(field_of_view_radius),
sphere_field_of_view_uycut_(-std::cos(field_of_view_radius)),
global_to_fp_(Eigen::Matrix3d::Identity()), ref_index_(ref_index)
{
// nothing to see here
}

virtual ~AllSkyVCLFocalPlaneRayPropagator() {
// nothing to see here
}

virtual std::vector<calin::simulation::ray_processor::RayProcessorDetectorSphere> detector_spheres() {
std::vector<calin::simulation::ray_processor::RayProcessorDetectorSphere> spheres;
spheres.emplace_back(sphere_center_, sphere_radius_,
global_to_fp_.transpose() * Eigen::Vector3d::UnitY(), M_PI/2, observation_level_);
return spheres;
}

void start_propagating() final {
// nothing to see here
}

void point_telescope_az_el_phi_deg(unsigned iscope,
double az_deg, double el_deg, double phi_deg) final {
if(iscope >= 1) {
throw std::out_of_range("Telescope number out of range");
}
global_to_fp_ =
Eigen::AngleAxisd(phi_deg*M_PI/180.0, Eigen::Vector3d::UnitZ()) *
Eigen::AngleAxisd(-el_deg*M_PI/180.0, Eigen::Vector3d::UnitX()) *
Eigen::AngleAxisd(az_deg*M_PI/180.0, Eigen::Vector3d::UnitZ());
}

#ifndef SWIG
double_bvt propagate_rays_to_focal_plane(
unsigned scope_id, Ray_vt& ray, double_bvt ray_mask,
VCLFocalPlaneParameters<VCLArchitecture>& fp_parameters) final {
if(scope_id >= 1) {
throw std::out_of_range("Telescope number out of range");
}
ray.translate_origin(sphere_center_.template cast<double_vt>());
ray_mask = ray.propagate_to_z_plane_with_mask(ray_mask, /* d= */ sphere_center_.z(),
/* time_reversal_ok = */ false, ref_index_);
ray.rotate(global_to_fp_.template cast<double_vt>());
ray_mask &= ray.x()*ray.x() + ray.z()*ray.z() < sphere_radius_squared_;
ray_mask &= ray.uy() <= sphere_field_of_view_uycut_;
TraceInfo info;
fp_parameters.fplane_x = select(ray_mask, x, 0);
fp_parameters.fplane_y = select(ray_mask, y, 0); // ** UNUSED **
fp_parameters.fplane_z = select(ray_mask, z, 0);
fp_parameters.fplane_ux = select(ray_mask, ray.ux(), 0);
fp_parameters.fplane_uy = select(ray_mask, ray.uy(), 0); // ** UNUSED **
fp_parameters.fplane_uz = select(ray_mask, ray.uz(), 0);
fp_parameters.fplane_t = select(ray_mask, ray.time(), 0);
fp_parameters.pixel_id = 0;
fp_parameters.detection_prob = 1.0;
return ray_mask;
}
#endif

void finish_propagating() final {
// nothing to see here
}

std::string banner(const std::string& indent0 = "", const std::string& indentN = "") const final {
std::ostringstream stream;
stream << indent0 << "All-sky detector.\n";
return stream.str();
}

#ifndef SWIG
private:
unsigned observation_level_; // Observation layer associated with this detector
Eigen::Vector3d sphere_center_; // Center of detector sphere [cm]
double sphere_radius_; // Radius of sphere [cm]
double sphere_radius_squared_; // Squared radius of sphere [cm^2]
double sphere_field_of_view_radius_; // Field of view of detector [radians]
double sphere_field_of_view_uycut_; // Value of cos(uy) that coresponds to FoV radius
Eigen::Matrix3d global_to_fp_; // Rotation matrix from global to focal plane

double ref_index_ = 1.0;
#endif
};

} } } // namespace calin::simulations::vcl_ray_propagator
10 changes: 7 additions & 3 deletions swig/simulation/ray_propagator.i
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,10 @@
%template (DaviesCottonVCLFocalPlaneRayPropagator256) calin::simulation::vcl_ray_propagator::DaviesCottonVCLFocalPlaneRayPropagator<calin::util::vcl::VCL256Architecture>;
%template (DaviesCottonVCLFocalPlaneRayPropagator512) calin::simulation::vcl_ray_propagator::DaviesCottonVCLFocalPlaneRayPropagator<calin::util::vcl::VCL512Architecture>;

%template (TrivialVCLFocalPlaneRayPropagator128) calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator<calin::util::vcl::VCL128Architecture>;
%template (TrivialVCLFocalPlaneRayPropagator256) calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator<calin::util::vcl::VCL256Architecture>;
%template (TrivialVCLFocalPlaneRayPropagator512) calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator<calin::util::vcl::VCL512Architecture>;
%template (PerfectOpticsVCLFocalPlaneRayPropagator128) calin::simulation::vcl_ray_propagator::PerfectOpticsVCLFocalPlaneRayPropagator<calin::util::vcl::VCL128Architecture>;
%template (PerfectOpticsVCLFocalPlaneRayPropagator256) calin::simulation::vcl_ray_propagator::PerfectOpticsVCLFocalPlaneRayPropagator<calin::util::vcl::VCL256Architecture>;
%template (PerfectOpticsVCLFocalPlaneRayPropagator512) calin::simulation::vcl_ray_propagator::PerfectOpticsVCLFocalPlaneRayPropagator<calin::util::vcl::VCL512Architecture>;

%template (AllSkyVCLFocalPlaneRayPropagator128) calin::simulation::vcl_ray_propagator::AllSkyVCLFocalPlaneRayPropagator<calin::util::vcl::VCL128Architecture>;
%template (AllSkyVCLFocalPlaneRayPropagator256) calin::simulation::vcl_ray_propagator::AllSkyVCLFocalPlaneRayPropagator<calin::util::vcl::VCL256Architecture>;
%template (AllSkyVCLFocalPlaneRayPropagator512) calin::simulation::vcl_ray_propagator::AllSkyVCLFocalPlaneRayPropagator<calin::util::vcl::VCL512Architecture>;

0 comments on commit dff89d6

Please sign in to comment.