diff --git a/include/simulation/vcl_ray_propagator.hpp b/include/simulation/vcl_ray_propagator.hpp index 5f20a266..e096a6dc 100644 --- a/include/simulation/vcl_ray_propagator.hpp +++ b/include/simulation/vcl_ray_propagator.hpp @@ -35,6 +35,7 @@ #include #include #include +#include namespace calin { namespace simulation { namespace vcl_ray_propagator { @@ -201,4 +202,154 @@ template class alignas(VCLArchitecture::vec_bytes) Dav #endif }; +template class alignas(VCLArchitecture::vec_bytes) TrivialVCLFocalPlaneRayPropagator: + public VCLFocalPlaneRayPropagator +{ +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; + using Ray_vt = typename calin::math::ray::VCLRay; + + using RayTracer = calin::simulation::vcl_raytracer::VCLScopeRayTracer; + using TraceInfo = calin::simulation::vcl_raytracer::VCLScopeTraceInfo; + using RealRNG = calin::math::rng::VCLRealRNG; +#endif // not defined SWIG + + CALIN_TYPEALIAS(ArchRNG, calin::math::rng::VCLRNG); + + 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)), + ref_index_(ref_index) + { + // nothing to see here + } + + virtual ~TrivialVCLFocalPlaneRayPropagator() { + for(auto* scope : scopes_) { + delete scope; + } + } + + void add_telecope(const Eigen::Vector3d& r0, double radius, unsigned iobs, + double focal_length = 0, double field_of_view_radius = M_PI/2) { + TelescopeDetails* scope = new TelescopeDetails; + scope->r0 = r0; + scope->radius = radius; + scope->radius_squared = radius*radius; + scope->iobs = iobs; + scope->obs_dir = Eigen::Vector3d::UnitY(); + scope->field_of_view_radius = field_of_view_radius; + scope->field_of_view_uycut = std::cos(field_of_view_radius); + scope->global_to_fp = Eigen::Matrix3d::Identity(); + scope->focal_length = focal_length; + scopes_.push_back(scope); + } + + virtual std::vector detector_spheres() { + std::vector spheres; + spheres.reserve(scopes_.size()); + for(const auto* scope : scopes_) { + spheres.emplace_back(scope->r0, scope->radius, scope->obs_dir, scope->field_of_view_radius, + scope->iobs); + } + 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 >= scopes_.size()) { + throw std::out_of_range("Telescope number out of range"); + } + auto* scope = scopes_[iscope]; + scope->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()); + scope->obs_dir = scope->global_to_fp.transpose() * Eigen::Vector3d::UnitY(); + } + +#ifndef SWIG + double_bvt propagate_rays_to_focal_plane( + unsigned scope_id, Ray_vt& ray, double_bvt ray_mask, + VCLFocalPlaneParameters& fp_parameters) final { + if(scope_id >= scopes_.size()) { + throw std::out_of_range("Telescope number out of range"); + } + const auto* scope = scopes_[scope_id]; + ray.translate_origin(scope->r0.template cast()); + ray.rotate(scope->global_to_fp.template cast()); + double_vt uy = -ray.uy(); + ray_mask = ray.propagate_to_y_plane_with_mask(ray_mask, /* d= */ 0, /* time_reversal_ok = */ false, ref_index_); + ray_mask &= ray.x()*ray.x() + ray.z()*ray.z() < scope->radius_squared; + ray_mask &= uy >= scope->field_of_view_uycut; + double_vt F_over_uy = scope->focal_length/uy; + TraceInfo info; + fp_parameters.fplane_x = select(ray_mask, F_over_uy * ray.ux(), 0); + fp_parameters.fplane_y = 0.0; + fp_parameters.fplane_z = select(ray_mask, F_over_uy * ray.uz(), 0); + fp_parameters.fplane_ux = select(ray_mask, ray.ux(), 0); + fp_parameters.fplane_uy = select(ray_mask, uy, 0); + 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 { + using calin::util::string::double_to_string_with_commas; + double xmin = std::numeric_limits::infinity(); + double xmax = -std::numeric_limits::infinity(); + double ymin = std::numeric_limits::infinity(); + double ymax = -std::numeric_limits::infinity(); + for(const auto* scope : scopes_) { + xmin = std::min(xmin, scope->r0[0]); + xmax = std::max(xmax, scope->r0[0]); + ymin = std::min(ymin, scope->r0[1]); + ymax = std::max(ymax, scope->r0[1]); + } + double A = (xmax-xmin)*(ymax-ymin)*1e-10; + std::ostringstream stream; + stream << indent0 << "Array of " << this->scopes_.size() + << " pseudo-telescopes covering " + << double_to_string_with_commas(A,3) << " km^2.\n"; + return stream.str(); + } + +#ifndef SWIG +private: + + struct TelescopeDetails + { + Eigen::Vector3d r0; // Center of detector sphere [cm] + double radius; // Radius of sphere [cm] + double radius_squared; // Squared radius of sphere [cm^2] + unsigned iobs; // Observation layer associated with this detector + Eigen::Vector3d obs_dir; // Pointing direction of detector + double field_of_view_radius; // Field of view of detector [radians] + double 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 focal_length; // Focal length [cm] + }; + + RealRNG* rng_ = nullptr; + double ref_index_ = 1.0; + std::vector scopes_; +#endif +}; + } } } // namespace calin::simulations::vcl_ray_propagator diff --git a/swig/simulation/ray_propagator.i b/swig/simulation/ray_propagator.i index 1e1135a9..a468ce9f 100644 --- a/swig/simulation/ray_propagator.i +++ b/swig/simulation/ray_propagator.i @@ -50,3 +50,7 @@ %template (DaviesCottonVCLFocalPlaneRayPropagator128) calin::simulation::vcl_ray_propagator::DaviesCottonVCLFocalPlaneRayPropagator; %template (DaviesCottonVCLFocalPlaneRayPropagator256) calin::simulation::vcl_ray_propagator::DaviesCottonVCLFocalPlaneRayPropagator; %template (DaviesCottonVCLFocalPlaneRayPropagator512) calin::simulation::vcl_ray_propagator::DaviesCottonVCLFocalPlaneRayPropagator; + +%template (TrivialVCLFocalPlaneRayPropagator128) calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator; +%template (TrivialVCLFocalPlaneRayPropagator256) calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator; +%template (TrivialVCLFocalPlaneRayPropagator512) calin::simulation::vcl_ray_propagator::TrivialVCLFocalPlaneRayPropagator;