Skip to content

Commit

Permalink
Add TrivialVCLFocalPlaneRayPropagator
Browse files Browse the repository at this point in the history
  • Loading branch information
@sfegan
sfegan committed Oct 10, 2023
1 parent 4def29c commit be2c586
Show file tree
Hide file tree
Showing 2 changed files with 155 additions and 0 deletions.
151 changes: 151 additions & 0 deletions include/simulation/vcl_ray_propagator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
#include <simulation/ray_processor.hpp>
#include <simulation/vso_ray_processor.hpp>
#include <util/log.hpp>
#include <util/string.hpp>

namespace calin { namespace simulation { namespace vcl_ray_propagator {

Expand Down Expand Up @@ -201,4 +202,154 @@ template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) Dav
#endif
};

template<typename VCLArchitecture> class alignas(VCLArchitecture::vec_bytes) TrivialVCLFocalPlaneRayPropagator:
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 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)),
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<calin::simulation::ray_processor::RayProcessorDetectorSphere> detector_spheres() {
std::vector<calin::simulation::ray_processor::RayProcessorDetectorSphere> 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<VCLArchitecture>& 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<double_vt>());
ray.rotate(scope->global_to_fp.template cast<double_vt>());
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<double>::infinity();
double xmax = -std::numeric_limits<double>::infinity();
double ymin = std::numeric_limits<double>::infinity();
double ymax = -std::numeric_limits<double>::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<TelescopeDetails*> scopes_;
#endif
};

} } } // namespace calin::simulations::vcl_ray_propagator
4 changes: 4 additions & 0 deletions swig/simulation/ray_propagator.i
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,3 +50,7 @@
%template (DaviesCottonVCLFocalPlaneRayPropagator128) calin::simulation::vcl_ray_propagator::DaviesCottonVCLFocalPlaneRayPropagator<calin::util::vcl::VCL128Architecture>;
%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>;

0 comments on commit be2c586

Please sign in to comment.