From 0d5b01185e08c8a22770e43de78cf6e493b88b1f Mon Sep 17 00:00:00 2001 From: Hiiragi Date: Sat, 25 Sep 2021 15:22:14 +0900 Subject: [PATCH] use pointer instead of object for Ray and ParaxialRay --- .../field_curvature_setting_dialog.cpp | 6 +- .../Analysis/longitudinal_setting_dialog.cpp | 6 +- .../gui/Analysis/paraxial_trace_dialog.cpp | 4 +- .../gui/Analysis/single_ray_trace_dialog.cpp | 12 +- .../Analysis/spot_diagram_setting_dialog.cpp | 9 +- .../Analysis/spot_diagram_setting_dialog.ui | 68 +++--- .../Analysis/transverse_ray_fan_dialog.cpp | 6 +- geopter/gui/main_window.cpp | 2 +- geopter/gui/plot_view_dock.cpp | 1 - geopter/optical/include/Analysis/aberration.h | 10 +- geopter/optical/include/Analysis/layout.h | 2 +- .../include/Assembly/even_polynomial.h | 4 +- geopter/optical/include/Assembly/interface.h | 2 +- geopter/optical/include/Assembly/spherical.h | 6 +- geopter/optical/include/Assembly/surface.h | 2 +- .../include/Assembly/surface_profile.h | 10 +- .../optical/include/Paraxial/paraxial_ray.h | 14 +- .../optical/include/Paraxial/paraxial_trace.h | 8 +- geopter/optical/include/Sequential/ray.h | 16 +- .../include/Sequential/sequential_trace.h | 21 +- .../optical/include/System/optical_system.h | 32 +-- geopter/optical/include/optical.h | 4 +- geopter/optical/src/Analysis/aberration.cpp | 210 +++++++++++------- geopter/optical/src/Analysis/layout.cpp | 18 +- geopter/optical/src/Analysis/longitudinal.cpp | 10 +- geopter/optical/src/Analysis/ray_fan.cpp | 23 +- .../optical/src/Assembly/even_polynomial.cpp | 6 +- geopter/optical/src/Assembly/interface.cpp | 2 +- geopter/optical/src/Assembly/spherical.cpp | 6 +- geopter/optical/src/Assembly/surface.cpp | 2 +- .../optical/src/Assembly/surface_profile.cpp | 10 +- geopter/optical/src/CMakeLists.txt | 4 +- geopter/optical/src/Paraxial/paraxial_ray.cpp | 31 +-- .../optical/src/Paraxial/paraxial_trace.cpp | 33 +-- geopter/optical/src/Sequential/ray.cpp | 52 +++-- .../src/Sequential/sequential_trace.cpp | 132 ++++++----- geopter/optical/src/Spec/field.cpp | 25 +-- geopter/optical/src/System/optical_system.cpp | 84 +++---- geopter/optical/test/triplet.cpp | 12 +- 39 files changed, 503 insertions(+), 402 deletions(-) diff --git a/geopter/gui/Analysis/field_curvature_setting_dialog.cpp b/geopter/gui/Analysis/field_curvature_setting_dialog.cpp index ffddbeb..ded1a14 100644 --- a/geopter/gui/Analysis/field_curvature_setting_dialog.cpp +++ b/geopter/gui/Analysis/field_curvature_setting_dialog.cpp @@ -34,9 +34,9 @@ void FieldCurvatureSettingDialog::updateParentDockContent() m_opticalSystem->update_model(); - Longitudinal *lon = new Longitudinal(m_opticalSystem, m_renderer); - lon->plot_ast(scale); - delete lon; + Aberration *abr = new Aberration(m_opticalSystem, m_renderer); + abr->plot_astigmatism(scale); + delete abr; m_renderer->update(); } diff --git a/geopter/gui/Analysis/longitudinal_setting_dialog.cpp b/geopter/gui/Analysis/longitudinal_setting_dialog.cpp index ed5d821..65e6409 100644 --- a/geopter/gui/Analysis/longitudinal_setting_dialog.cpp +++ b/geopter/gui/Analysis/longitudinal_setting_dialog.cpp @@ -34,9 +34,9 @@ void LongitudinalSettingDialog::updateParentDockContent() m_opticalSystem->update_model(); - Longitudinal *lon = new Longitudinal(m_opticalSystem, m_renderer); - lon->plot_lsa(scale); - delete lon; + Aberration *abr = new Aberration(m_opticalSystem, m_renderer); + abr->plot_longitudinal_spherical_aberration(scale); + delete abr; m_renderer->update(); } diff --git a/geopter/gui/Analysis/paraxial_trace_dialog.cpp b/geopter/gui/Analysis/paraxial_trace_dialog.cpp index 9b0d347..0651058 100644 --- a/geopter/gui/Analysis/paraxial_trace_dialog.cpp +++ b/geopter/gui/Analysis/paraxial_trace_dialog.cpp @@ -45,11 +45,11 @@ void ParaxialTraceDialog::updateParentDockContent() oss << std::endl; oss << "Axial Ray..." << std::endl; - m_opticalSystem->axial_ray(wi).print(oss); + m_opticalSystem->axial_ray(wi)->print(oss); oss << std::endl; oss << "Principle Ray..." << std::endl; - m_opticalSystem->principle_ray(wi).print(oss); + m_opticalSystem->principle_ray(wi)->print(oss); oss << std::endl; m_parentDock->setStringStreamToText(oss); diff --git a/geopter/gui/Analysis/single_ray_trace_dialog.cpp b/geopter/gui/Analysis/single_ray_trace_dialog.cpp index 2e3e1cf..2cf4f9e 100644 --- a/geopter/gui/Analysis/single_ray_trace_dialog.cpp +++ b/geopter/gui/Analysis/single_ray_trace_dialog.cpp @@ -91,15 +91,17 @@ void SingleRayTraceDialog::doPupilRayTrace() double px = ui->pupilXEdit->text().toDouble(); double py = ui->pupilYEdit->text().toDouble(); int fi = ui->fieldCombo->currentIndex(); + Field* fld = m_opticalSystem->optical_spec()->field_of_view()->field(fi); int wi = ui->wvlForPupilCombo->currentIndex(); + double wvl = m_opticalSystem->optical_spec()->spectral_region()->wvl(wi)->value(); //Field *fld = opt_sys_->optical_spec()->field_of_view()->field(fi); Eigen::Vector2d pupil_crd({px, py}); - double wvl = m_opticalSystem->optical_spec()->spectral_region()->wvl(wi)->value(); + // trace SequentialTrace *tracer = new SequentialTrace(m_opticalSystem); - Ray ray_trace_result = tracer->trace_pupil_ray(pupil_crd, fi, wi); + std::shared_ptr ray_trace_result = tracer->trace_pupil_ray(pupil_crd, fld, wvl); delete tracer; // construct output text @@ -109,7 +111,7 @@ void SingleRayTraceDialog::doPupilRayTrace() oss << "Field: " << fi << std::endl; oss << "Wavelength: " << wi << " " << wvl << std::endl; - ray_trace_result.print(oss); + ray_trace_result->print(oss); oss << std::endl; // write to textview dock @@ -136,7 +138,7 @@ void SingleRayTraceDialog::doObjectRayTrace() double wvl = m_opticalSystem->optical_spec()->spectral_region()->wvl(wi)->value(); SequentialTrace *tracer = new SequentialTrace(m_opticalSystem); - Ray ray_trace_result = tracer->trace_ray_throughout_path(tracer->overall_sequential_path(wi), p0, dir0); + auto ray_trace_result = tracer->trace_ray_throughout_path(tracer->overall_sequential_path(wi), p0, dir0); delete tracer; std::ostringstream oss; @@ -145,7 +147,7 @@ void SingleRayTraceDialog::doObjectRayTrace() oss << "Object Space Direction : " << "(" << L << ", " << M << ", " << N << ")" << std::endl; oss << "Wavelength: " << wi << " " << wvl << std::endl; - ray_trace_result.print(oss); + ray_trace_result->print(oss); oss << std::endl; m_parentDock->setStringStreamToText(oss); diff --git a/geopter/gui/Analysis/spot_diagram_setting_dialog.cpp b/geopter/gui/Analysis/spot_diagram_setting_dialog.cpp index 313a1af..4d484ca 100644 --- a/geopter/gui/Analysis/spot_diagram_setting_dialog.cpp +++ b/geopter/gui/Analysis/spot_diagram_setting_dialog.cpp @@ -13,8 +13,12 @@ SpotDiagramSettingDialog::SpotDiagramSettingDialog(OpticalSystem* sys, PlotViewD ui->setupUi(this); m_renderer = new RendererQCP(m_parentDock->customPlot()); + + ui->rayPatternCombo->clear(); + ui->rayPatternCombo->addItems(QStringList({"Grid", "Hexapolar"})); + ui->rayPatternCombo->setCurrentIndex(0); ui->nrdEdit->setValidator(new QIntValidator(1,1000,this)); - ui->nrdEdit->setText(QString::number(30)); + ui->nrdEdit->setText(QString::number(20)); ui->scaleEdit->setValidator(new QDoubleValidator(0.0,1000.0,4,this)); ui->scaleEdit->setText(QString::number(0.1)); ui->dotSizeEdit->setValidator(new QDoubleValidator(0.0,1000.0,4,this)); @@ -31,6 +35,7 @@ void SpotDiagramSettingDialog::updateParentDockContent() { m_renderer->clear(); + int pattern = ui->rayPatternCombo->currentIndex(); int nrd = ui->nrdEdit->text().toInt(); double scale = ui->scaleEdit->text().toDouble(); double dotSize = ui->dotSizeEdit->text().toDouble(); @@ -38,7 +43,7 @@ void SpotDiagramSettingDialog::updateParentDockContent() m_opticalSystem->update_model(); Aberration *abr = new Aberration(m_opticalSystem, m_renderer); - abr->plot_spot_diagram(nrd,scale,dotSize); + abr->plot_spot_diagram(pattern, nrd, scale, dotSize); delete abr; m_renderer->update(); diff --git a/geopter/gui/Analysis/spot_diagram_setting_dialog.ui b/geopter/gui/Analysis/spot_diagram_setting_dialog.ui index 1efa3f8..2f91638 100644 --- a/geopter/gui/Analysis/spot_diagram_setting_dialog.ui +++ b/geopter/gui/Analysis/spot_diagram_setting_dialog.ui @@ -14,24 +14,10 @@ Dialog - - - - Qt::Horizontal - - - QDialogButtonBox::Cancel|QDialogButtonBox::Ok - - - - - - - Ray Density: - - + + - + Qt::Horizontal @@ -44,7 +30,7 @@ - + @@ -54,14 +40,7 @@ - - - - Scale: - - - - + Qt::Horizontal @@ -74,7 +53,7 @@ - + @@ -84,16 +63,47 @@ - - + + + + Qt::Horizontal + + + QDialogButtonBox::Cancel|QDialogButtonBox::Ok + + + + + Scale: + + + + Dot Size: + + + + Ray Density: + + + + + + + + + + Pattern: + + + diff --git a/geopter/gui/Analysis/transverse_ray_fan_dialog.cpp b/geopter/gui/Analysis/transverse_ray_fan_dialog.cpp index 951a245..2bd7077 100644 --- a/geopter/gui/Analysis/transverse_ray_fan_dialog.cpp +++ b/geopter/gui/Analysis/transverse_ray_fan_dialog.cpp @@ -37,9 +37,9 @@ void TransverseRayFanDialog::updateParentDockContent() m_opticalSystem->update_model(); - RayFan *ray_fan = new RayFan(m_opticalSystem, m_renderer); - ray_fan->plot(scale, nrd); - delete ray_fan; + Aberration *abr = new Aberration(m_opticalSystem, m_renderer); + abr->plot_transverse_aberration(scale, nrd); + delete abr; m_renderer->update(); } diff --git a/geopter/gui/main_window.cpp b/geopter/gui/main_window.cpp index 39d38d0..92bae5f 100644 --- a/geopter/gui/main_window.cpp +++ b/geopter/gui/main_window.cpp @@ -249,7 +249,7 @@ void MainWindow::showSpotDiagram() PlotViewDock *dock = new PlotViewDock("Spot Diagram", opt_sys_.get()); dock->createSettingDialog(); m_dockManager->addDockWidgetFloating(dock); - dock->resize(300,200); + //dock->resize(300,200); dock->updatePlot(); } diff --git a/geopter/gui/plot_view_dock.cpp b/geopter/gui/plot_view_dock.cpp index b1991db..a6026af 100644 --- a/geopter/gui/plot_view_dock.cpp +++ b/geopter/gui/plot_view_dock.cpp @@ -13,7 +13,6 @@ PlotViewDock::PlotViewDock(QString label, OpticalSystem* sys, QWidget *parent): { this->setFeature(CDockWidget::DockWidgetDeleteOnClose, true); this->setMinimumSizeHintMode(CDockWidget::MinimumSizeHintFromDockWidget); - //this->resize(300,200); this->setMinimumSize(300,200); // QCustomPlot diff --git a/geopter/optical/include/Analysis/aberration.h b/geopter/optical/include/Analysis/aberration.h index 6a40edf..c56a673 100644 --- a/geopter/optical/include/Analysis/aberration.h +++ b/geopter/optical/include/Analysis/aberration.h @@ -25,9 +25,17 @@ class Aberration void plot_chromatic_focus_shift(double lower_wvl, double higher_wvl); - void plot_spot_diagram(int nrd, double scale, double dot_size); + void plot_spot_diagram(int pattern, int nrd, double scale, double dot_size); private: + std::vector create_grid_circle(int nrd); + std::vector create_hexapolar(int nrd); + + enum SpotRayPattern{ + Grid, + Hexapolar + }; + OpticalSystem* opt_sys_; Renderer* renderer_; int num_fld_; diff --git a/geopter/optical/include/Analysis/layout.h b/geopter/optical/include/Analysis/layout.h index 16952ac..8548225 100644 --- a/geopter/optical/include/Analysis/layout.h +++ b/geopter/optical/include/Analysis/layout.h @@ -27,7 +27,7 @@ class Layout void draw_fan_rays(int nrd= 10); /** Draw a single ray */ - void draw_single_ray(const Ray& ray, const Rgb& color); + void draw_single_ray(const std::shared_ptr& ray, const Rgb& color); void update(); diff --git a/geopter/optical/include/Assembly/even_polynomial.h b/geopter/optical/include/Assembly/even_polynomial.h index 9d9b2dc..c9f0a0e 100644 --- a/geopter/optical/include/Assembly/even_polynomial.h +++ b/geopter/optical/include/Assembly/even_polynomial.h @@ -30,8 +30,8 @@ class EvenPolynomial : public SurfaceProfile int coef_count() const; double sag(double x, double y) const override; - double f(Eigen::Vector3d p) const override; - Eigen::Vector3d df(Eigen::Vector3d p) const override; + double f(const Eigen::Vector3d& p) const override; + Eigen::Vector3d df(const Eigen::Vector3d& p) const override; double deriv_1st(double h) const override; double deriv_2nd(double h) const override; diff --git a/geopter/optical/include/Assembly/interface.h b/geopter/optical/include/Assembly/interface.h index 43d0d9c..bbc5e7d 100644 --- a/geopter/optical/include/Assembly/interface.h +++ b/geopter/optical/include/Assembly/interface.h @@ -57,7 +57,7 @@ class Interface /** Returns true if the given point(x,y) is inside of aperture */ bool point_inside(double x, double y) const; - bool point_inside(Eigen::Vector2d pt) const; + bool point_inside(const Eigen::Vector2d& pt) const; void set_local_transform(Transformation tfrm); void set_global_transform(Transformation tfrm); diff --git a/geopter/optical/include/Assembly/spherical.h b/geopter/optical/include/Assembly/spherical.h index 9f8fc11..e5fdb17 100644 --- a/geopter/optical/include/Assembly/spherical.h +++ b/geopter/optical/include/Assembly/spherical.h @@ -12,11 +12,11 @@ class Spherical : public SurfaceProfile Spherical(double c=0.0); ~Spherical(); - double f(Eigen::Vector3d p) const override; - Eigen::Vector3d df(Eigen::Vector3d p) const override; + double f(const Eigen::Vector3d& p) const override; + Eigen::Vector3d df(const Eigen::Vector3d& p) const override; double sag(double x, double y) const override; - void intersect(Eigen::Vector3d& pt, double& s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps=1.0e-12, int z_dir=1.0) override; + void intersect(Eigen::Vector3d& pt, double& s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps=1.0e-12, int z_dir=1.0) override; }; diff --git a/geopter/optical/include/Assembly/surface.h b/geopter/optical/include/Assembly/surface.h index 317149f..31e3560 100644 --- a/geopter/optical/include/Assembly/surface.h +++ b/geopter/optical/include/Assembly/surface.h @@ -24,7 +24,7 @@ class Surface : public Interface * @param eps precision * @param z_dir */ - void intersect(Eigen::Vector3d& pt, double& s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps=1.0e-12, double z_dir=1.0); + void intersect(Eigen::Vector3d& pt, double& s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps=1.0e-12, double z_dir=1.0); Eigen::Vector3d normal(Eigen::Vector3d p); diff --git a/geopter/optical/include/Assembly/surface_profile.h b/geopter/optical/include/Assembly/surface_profile.h index f4c6c6d..0592a77 100644 --- a/geopter/optical/include/Assembly/surface_profile.h +++ b/geopter/optical/include/Assembly/surface_profile.h @@ -35,13 +35,13 @@ class SurfaceProfile virtual double radius() const; /** Returns the value of the profile function at point @param{p} */ - virtual double f(Eigen::Vector3d p) const; + virtual double f(const Eigen::Vector3d& p) const; /** Returns the gradient of the profile function at point *p* */ - virtual Eigen::Vector3d df(Eigen::Vector3d p) const; + virtual Eigen::Vector3d df(const Eigen::Vector3d& p) const; /** Returns the unit normal of the profile at point *p* */ - virtual Eigen::Vector3d normal(Eigen::Vector3d p) const; + virtual Eigen::Vector3d normal(const Eigen::Vector3d& p) const; /** Returns the sagitta (z coordinate) of the surface at x, y */ virtual double sag(double x, double y) const; @@ -57,7 +57,7 @@ class SurfaceProfile * @param eps numeric tolerance for convergence of any iterative procedure * @param z_dir +1 if propagation positive direction, -1 if otherwise */ - virtual void intersect(Eigen::Vector3d& pt, double& s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps=1.0e-12, int z_dir=1); + virtual void intersect(Eigen::Vector3d& pt, double& s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps=1.0e-12, int z_dir=1); /** * @brief Intersect a profile, starting from an arbitrary point @@ -69,7 +69,7 @@ class SurfaceProfile * @param eps numeric tolerance for convergence of any iterative procedure * @param z_dir +1 if propagation positive direction, -1 if otherwise */ - virtual void intersect_spencer(Eigen::Vector3d& pt, double& s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps=1.0e-12, int z_dir=1.0); + virtual void intersect_spencer(Eigen::Vector3d& pt, double& s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps=1.0e-12, int z_dir=1.0); protected: diff --git a/geopter/optical/include/Paraxial/paraxial_ray.h b/geopter/optical/include/Paraxial/paraxial_ray.h index d8f432e..5cdb048 100644 --- a/geopter/optical/include/Paraxial/paraxial_ray.h +++ b/geopter/optical/include/Paraxial/paraxial_ray.h @@ -42,10 +42,10 @@ class ParaxialRay ParaxialRay(); ~ParaxialRay(); - void prepend(ParaxialRayAtSurface par_ray_at_srf); + void prepend(std::shared_ptr par_ray_at_srf); void prepend(double ht, double slp, double aoi= 0, double n = 1.0); - void append(ParaxialRayAtSurface par_ray_at_srf); + void append(std::shared_ptr par_ray_at_srf); void append(double ht, double slp, double aoi= 0, double n = 1.0); void clear(); @@ -57,16 +57,16 @@ class ParaxialRay int size() const; /** Access to data at the given index */ - ParaxialRayAtSurface at(int i) const; + std::shared_ptr at(int i) const; /** Returns the data at the beginning */ - ParaxialRayAtSurface front() const; + std::shared_ptr front() const; /** Return the data at the last surface */ - ParaxialRayAtSurface back() const; + std::shared_ptr back() const; /** Access to the data at the surface just before the image */ - ParaxialRayAtSurface at_before_image() const; + std::shared_ptr at_before_image() const; /** Write ray data to standard output */ void print() const; @@ -75,7 +75,7 @@ class ParaxialRay void print(std::ostringstream& oss) const; private: - std::vector par_ray_at_srfs_; + std::vector< std::shared_ptr > par_ray_at_srfs_; std::string name_; }; diff --git a/geopter/optical/include/Paraxial/paraxial_trace.h b/geopter/optical/include/Paraxial/paraxial_trace.h index 50b64bb..efa0598 100644 --- a/geopter/optical/include/Paraxial/paraxial_trace.h +++ b/geopter/optical/include/Paraxial/paraxial_trace.h @@ -21,15 +21,15 @@ class ParaxialTrace * @param y0 paraxial ray height at starting point * @param u0 slope of outgoing ray at starting point */ - ParaxialRay trace_paraxial_ray_from_object(double y0, double u0, double wvl); + std::shared_ptr trace_paraxial_ray_from_object(double y0, double u0, double wvl); - ParaxialRay trace_paraxial_ray_aiming_at_surface(int target_srf, double y_target, double u_prime_target ,double wvl); + std::shared_ptr trace_paraxial_ray_aiming_at_surface(int target_srf, double y_target, double u_prime_target ,double wvl); /** trace [y1,u1] = [1.0, 0.0] */ - ParaxialRay trace_paraxial_ray_parallel_to_axis_at_s1(); + std::shared_ptr trace_paraxial_ray_parallel_to_axis_at_s1(); /** trace [y1,u1] = [0.0, uq0] */ - ParaxialRay trace_paraxial_ray_with_slope_at_s1(); + std::shared_ptr trace_paraxial_ray_with_slope_at_s1(); /** * @brief Generate the paraxial path * diff --git a/geopter/optical/include/Sequential/ray.h b/geopter/optical/include/Sequential/ray.h index df2128f..e19e576 100644 --- a/geopter/optical/include/Sequential/ray.h +++ b/geopter/optical/include/Sequential/ray.h @@ -3,6 +3,7 @@ #include #include +#include #include "Eigen/Core" @@ -19,6 +20,7 @@ struct RayAtSurface after_dir = Eigen::Vector3d::Zero(3); distance_from_before = 0.0; optical_path_length = 0.0; + before = nullptr; } /** ray intersect point on the surface */ @@ -36,6 +38,8 @@ struct RayAtSurface /** optical path length from the previous to the current */ double optical_path_length; + RayAtSurface* before; + // aliases double x() const { return intersect_pt(0); } double y() const { return intersect_pt(1); } @@ -66,14 +70,14 @@ class Ray void clear(); /** Add data at the beginning */ - void prepend(RayAtSurface ray_at_srf); + void prepend(std::shared_ptr ray_at_srf); /** Add data at the last */ - void append(RayAtSurface ray_at_srf); + void append(std::shared_ptr ray_at_srf); - RayAtSurface at(int i) const; - RayAtSurface front() const; - RayAtSurface back() const; + RayAtSurface* at(int i) const; + RayAtSurface* front() const; + RayAtSurface* back() const; void set_wvl(double wvl); @@ -109,7 +113,7 @@ class Ray void print(); private: - std::vector ray_at_srfs_; + std::vector< std::shared_ptr > ray_at_srfs_; double wvl_; int status_; }; diff --git a/geopter/optical/include/Sequential/sequential_trace.h b/geopter/optical/include/Sequential/sequential_trace.h index e126479..d9d9e22 100644 --- a/geopter/optical/include/Sequential/sequential_trace.h +++ b/geopter/optical/include/Sequential/sequential_trace.h @@ -1,6 +1,8 @@ #ifndef SEQUENTIALTRACE_H #define SEQUENTIALTRACE_H +#include + #include "Eigen/Core" #include "Sequential/ray.h" @@ -16,16 +18,16 @@ class SequentialTrace SequentialTrace(OpticalSystem* sys); ~SequentialTrace(); - Ray trace_pupil_ray(Eigen::Vector2d pupil_crd, int fi, int wi); + //std::shared_ptr trace_pupil_ray(Eigen::Vector2d pupil_crd, int fi, int wi); - Ray trace_pupil_ray(Eigen::Vector2d pupil_crd, const Field* fld, double wvl); + std::shared_ptr trace_pupil_ray(const Eigen::Vector2d& pupil_crd, const Field* fld, double wvl); /** Trace a ray throughout the given sequantial path */ - Ray trace_ray_throughout_path(const SequentialPath& seq_path, Eigen::Vector3d pt0, Eigen::Vector3d dir0); + std::shared_ptr trace_ray_throughout_path(const SequentialPath& seq_path, const Eigen::Vector3d& pt0, const Eigen::Vector3d& dir0); Eigen::Vector2d trace_coddington(const Field* fld, double wvl); - Eigen::Vector2d aim_chief_ray(int fi, int wi); + //Eigen::Vector2d aim_chief_ray(int fi, int wi); Eigen::Vector2d aim_chief_ray(const Field* fld, double wvl); @@ -38,8 +40,8 @@ class SequentialTrace * @param wi wavelength index * @return Eigen::Vector2d aim point on paraxial entrance pupil plane */ - Eigen::Vector2d search_aim_point(int srf_idx, Eigen::Vector2d xy_target, const Field* fld, double wvl); - Eigen::Vector2d search_aim_point(int srf_idx, Eigen::Vector2d xy_target, int fi, int wi); + Eigen::Vector2d search_aim_point(int srf_idx, const Eigen::Vector2d& xy_target, const Field* fld, double wvl); + //Eigen::Vector2d search_aim_point(int srf_idx, Eigen::Vector2d xy_target, int fi, int wi); /** @brief Refract incoming direction, d_in, about normal * @param d_in incident direction @@ -49,7 +51,7 @@ class SequentialTrace Eigen::Vector3d bend(Eigen::Vector3d d_in, Eigen::Vector3d normal, double n_in, double n_out); /** Get object coordinate for the given field */ - Eigen::Vector3d object_coord(int fi); + //Eigen::Vector3d object_coord(int fi); Eigen::Vector3d object_coord(const Field* fld); @@ -80,6 +82,11 @@ class SequentialTrace double ref_wvl_val_; int ref_wvl_idx_; + int num_srfs_; + int num_gaps_; + int image_index_; + double object_distance_; + bool do_aperture_check_; bool do_apply_vig_; diff --git a/geopter/optical/include/System/optical_system.h b/geopter/optical/include/System/optical_system.h index 8bde28f..3b91061 100644 --- a/geopter/optical/include/System/optical_system.h +++ b/geopter/optical/include/System/optical_system.h @@ -47,11 +47,11 @@ class OpticalSystem void add_surface_and_gap(double r, double t, std::string mat_name); - ParaxialRay axial_ray(int wi) const; - ParaxialRay principle_ray(int wi) const; + std::shared_ptr axial_ray(int wi) const; + std::shared_ptr principle_ray(int wi) const; FirstOrderData first_order_data() const; - Ray reference_ray(int ri, int fi, int wi) const; + std::shared_ptr reference_ray(int ri, int fi, int wi) const; void update_vignetting_factors(); @@ -71,18 +71,22 @@ class OpticalSystem std::unique_ptr opt_spec_; std::unique_ptr material_lib_; + // fundamental data + int num_wvl_; + int num_fld_; + // -----> Paraxial Data /** parallel to axis at s1 */ - ParaxialRay p_ray_; + std::shared_ptr p_ray_; /** with slope at s1 */ - ParaxialRay q_ray_; + std::shared_ptr q_ray_; /** paraxial axial rays computed with all wavelengths */ - std::vector ax_rays_; + std::vector> ax_rays_; /** paraxial principle rays computed with all wavelengths */ - std::vector pr_rays_; + std::vector> pr_rays_; /** first order data computed with all wavelengths */ FirstOrderData fod_; @@ -92,19 +96,19 @@ class OpticalSystem // -----> Sequential Data /** reference rays */ - std::vector ref_rays1_; // chief ray - std::vector ref_rays2_; // upper meridional marginal - std::vector ref_rays3_; // lower meridional marginal - std::vector ref_rays4_; // upper sagittal marginal - std::vector ref_rays5_; // lower sagittal marginal + std::vector< std::shared_ptr > ref_rays1_; // chief ray + std::vector< std::shared_ptr > ref_rays2_; // upper meridional marginal + std::vector< std::shared_ptr > ref_rays3_; // lower meridional marginal + std::vector< std::shared_ptr > ref_rays4_; // upper sagittal marginal + std::vector< std::shared_ptr > ref_rays5_; // lower sagittal marginal /** focus shifts */ std::vector x_focus_shifts_; std::vector y_focus_shifts_; int to_ray_index(int fi, int wi) const { - int num_wvl = opt_spec_->spectral_region()->wvl_count(); - return num_wvl * (fi) + wi; + //int num_wvl = opt_spec_->spectral_region()->wvl_count(); + return num_wvl_ * (fi) + wi; } // <----- Sequential Data End diff --git a/geopter/optical/include/optical.h b/geopter/optical/include/optical.h index 673caf8..eb403c5 100644 --- a/geopter/optical/include/optical.h +++ b/geopter/optical/include/optical.h @@ -13,8 +13,8 @@ /** include all header files **/ #include "Analysis/aberration.h" #include "Analysis/layout.h" -#include "Analysis/ray_fan.h" -#include "Analysis/longitudinal.h" +//#include "Analysis/ray_fan.h" +//#include "Analysis/longitudinal.h" #include "Assembly/optical_assembly.h" #include "Assembly/aperture.h" diff --git a/geopter/optical/src/Analysis/aberration.cpp b/geopter/optical/src/Analysis/aberration.cpp index c6f963b..d60a8d2 100644 --- a/geopter/optical/src/Analysis/aberration.cpp +++ b/geopter/optical/src/Analysis/aberration.cpp @@ -38,43 +38,36 @@ Aberration::~Aberration() void Aberration::plot_transverse_aberration(double scale, double nrd) { - constexpr int num_interpolation = 50; - const int num_flds = opt_sys_->optical_spec()->field_of_view()->field_count(); - const int num_wvls = opt_sys_->optical_spec()->spectral_region()->wvl_count(); - - renderer_->set_grid_layout(num_flds, 2); + renderer_->set_grid_layout(num_fld_, 2); Eigen::Vector2d pupil; SequentialTrace *tracer = new SequentialTrace(opt_sys_); - tracer->set_aperture_check(true); - tracer->set_apply_vig(false); + tracer->set_aperture_check(false); + tracer->set_apply_vig(true); - for(int fi = 0; fi < num_flds; fi++) + for(int fi = 0; fi < num_fld_; fi++) { // trace chief ray - //Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); + Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); - Ray chief_ray = opt_sys_->reference_ray(1,fi,ref_wvl_idx_); - double x0 = chief_ray.back().x(); - double y0 = chief_ray.back().y(); + std::shared_ptr chief_ray = opt_sys_->reference_ray(1,fi,ref_wvl_idx_); + double x0 = chief_ray->back()->x(); + double y0 = chief_ray->back()->y(); // trace zonal rays for all wavelengths - for(int wi = 0; wi < num_wvls; wi++) + for(int wi = 0; wi < num_wvl_; wi++) { + double wvl = opt_sys_->optical_spec()->spectral_region()->wvl(wi)->value(); Rgb render_color = opt_sys_->optical_spec()->spectral_region()->wvl(wi)->render_color(); std::vector px, py; std::vector dx, dy; - /* + px.reserve(nrd); py.reserve(nrd); dx.reserve(nrd); dy.reserve(nrd); - */ - - Ray tangential_ray; - Ray sagittal_ray; for(int ri = 0; ri < nrd; ri++) { @@ -82,11 +75,11 @@ void Aberration::plot_transverse_aberration(double scale, double nrd) pupil(0) = 0.0; pupil(1) = -1.0 + (double)ri*2.0/(double)(nrd-1); - tangential_ray = tracer->trace_pupil_ray(pupil, fi, wi); - - if(tangential_ray.status() == RayStatus::PassThrough){ - double y = tangential_ray.back().y(); - py.push_back(pupil(1)); + auto tangential_ray = tracer->trace_pupil_ray(pupil, fld, wvl); + if(tangential_ray->status() == RayStatus::PassThrough){ + double y = tangential_ray->back()->y(); + Eigen::Vector2d vig_pupil = fld->apply_vignetting(pupil); + py.push_back(vig_pupil(1)); dy.push_back(y - y0); } @@ -94,48 +87,32 @@ void Aberration::plot_transverse_aberration(double scale, double nrd) pupil(0) = -1.0 + (double)ri*2.0/(double)(nrd-1); pupil(1) = 0.0; - sagittal_ray = tracer->trace_pupil_ray(pupil, fi, wi); - if(sagittal_ray.status() == RayStatus::PassThrough){ - double x = sagittal_ray.back().x(); - px.push_back(pupil(0)); + auto sagittal_ray = tracer->trace_pupil_ray(pupil, fld, wvl); + if(sagittal_ray->status() == RayStatus::PassThrough){ + double x = sagittal_ray->back()->x(); + Eigen::Vector2d vig_pupil = fld->apply_vignetting(pupil); + px.push_back(vig_pupil(0)); dx.push_back(x - x0); } } - // interpolation - tk::spline spline_tangential; - spline_tangential.set_points(py,dy); - - tk::spline spline_sagittal; - spline_sagittal.set_points(px,dx); - - std::vector spx(num_interpolation), spy(num_interpolation); - std::vector sdx(num_interpolation), sdy(num_interpolation); - double px_min = *std::min_element(px.begin(), px.end()); - double px_max = *std::max_element(px.begin(), px.end()); - double py_min = *std::min_element(py.begin(), py.end()); - double py_max = *std::max_element(py.begin(), py.end()); - for(int k = 0; k < num_interpolation; k++){ - spx[k] = px_min + (double)k*(px_max-px_min)/(double)(num_interpolation -1); - spy[k] = py_min + (double)k*(py_max-py_min)/(double)(num_interpolation -1); - sdx[k] = spline_sagittal(spx[k]); - sdy[k] = spline_tangential(spy[k]); - } // draw - renderer_->set_current_cell(num_flds - fi -1, 0); + renderer_->set_current_cell(num_fld_ - fi -1, 0); renderer_->set_y_axis_range(-scale, scale); renderer_->set_x_axis_range(-1.0, 1.0); renderer_->set_y_axis_label("dy"); - renderer_->draw_polyline(spy, sdy, render_color); + //renderer_->draw_polyline(spy, sdy, render_color); + renderer_->draw_polyline(py, dy, render_color); renderer_->draw_x_axis(); renderer_->draw_y_axis(); - renderer_->set_current_cell(num_flds - fi -1, 1); + renderer_->set_current_cell(num_fld_ - fi -1, 1); renderer_->set_y_axis_range(-scale, scale); renderer_->set_x_axis_range(-1.0, 1.0); renderer_->set_y_axis_label("dx"); - renderer_->draw_polyline(spx, sdx, render_color); + //renderer_->draw_polyline(spx, sdx, render_color); + renderer_->draw_polyline(px, dx, render_color); renderer_->draw_x_axis(); renderer_->draw_y_axis(); } @@ -162,7 +139,7 @@ void Aberration::plot_longitudinal_spherical_aberration(double scale) PupilCrd pupil; - Ray ray; + std::shared_ptr ray; Field* fld0 = opt_sys_->optical_spec()->field_of_view()->field(0); int ref_wvl_idx = opt_sys_->optical_spec()->spectral_region()->reference_index(); @@ -177,9 +154,9 @@ void Aberration::plot_longitudinal_spherical_aberration(double scale) int last_surf = num_srf - 1 -1; int num_wvl = opt_sys_->optical_spec()->spectral_region()->wvl_count(); for(int wi = 0; wi < num_wvl; wi++){ - ParaxialRay ax_ray = opt_sys_->axial_ray(wi); - double y = ax_ray.at(last_surf).ht; - double u_prime = ax_ray.at(img_srf).slp; + std::shared_ptr ax_ray = opt_sys_->axial_ray(wi); + double y = ax_ray->at(last_surf)->ht; + double u_prime = ax_ray->at(img_srf)->slp; double l_prime = -y/u_prime; l_primes.push_back(l_prime); } @@ -212,10 +189,10 @@ void Aberration::plot_longitudinal_spherical_aberration(double scale) //if(ray.status() == RayStatus::PassThrough){ py.push_back(pupil(1)); - int at_image = ray.size()-1; - double y = ray.y(at_image); - double M = ray.M(at_image); - double N = ray.N(at_image); + int at_image = ray->size()-1; + double y = ray->y(at_image); + double M = ray->M(at_image); + double N = ray->N(at_image); lsa.push_back(-y*N/M); y_list.push_back(y); @@ -267,9 +244,9 @@ void Aberration::plot_astigmatism(double scale) int last_surf = num_srf - 1 -1; int num_wvl = opt_sys_->optical_spec()->spectral_region()->wvl_count(); for(int wi = 0; wi < num_wvl; wi++){ - ParaxialRay ax_ray = opt_sys_->axial_ray(wi); - double y = ax_ray.at(last_surf).ht; - double u_prime = ax_ray.at(img_srf).slp; + std::shared_ptr ax_ray = opt_sys_->axial_ray(wi); + double y = ax_ray->at(last_surf)->ht; + double u_prime = ax_ray->at(img_srf)->slp; double l_prime = -y/u_prime; l_primes.push_back(l_prime); } @@ -346,9 +323,9 @@ void Aberration::plot_chromatic_focus_shift(double lower_wvl, double higher_wvl) ParaxialTrace *tracer = new ParaxialTrace(opt_sys_); tracer->get_starting_coords(&y0, &u0, &ybar0, &ubar0); - ParaxialRay ref_prx_ray = tracer->trace_paraxial_ray_from_object(y0, u0, ref_wvl); - double y = ref_prx_ray.at(last_surf).ht; - double u_prime = ref_prx_ray.at(img_srf).slp; + std::shared_ptr ref_prx_ray = tracer->trace_paraxial_ray_from_object(y0, u0, ref_wvl); + double y = ref_prx_ray->at(last_surf)->ht; + double u_prime = ref_prx_ray->at(img_srf)->slp; double ref_l_prime = -y/u_prime; @@ -356,10 +333,10 @@ void Aberration::plot_chromatic_focus_shift(double lower_wvl, double higher_wvl) for(int i = 0; i < num_data; i++) { double wvl = lower_wvl + i*wvl_step; - ParaxialRay prx_ray = tracer->trace_paraxial_ray_from_object(y0, u0, wvl); + std::shared_ptr prx_ray = tracer->trace_paraxial_ray_from_object(y0, u0, wvl); - y = prx_ray.at(last_surf).ht; - u_prime = prx_ray.at(img_srf).slp; + y = prx_ray->at(last_surf)->ht; + u_prime = prx_ray->at(img_srf)->slp; double l_prime = -y/u_prime; ydata.push_back(l_prime - ref_l_prime); xdata.push_back(wvl); @@ -382,7 +359,7 @@ void Aberration::plot_chromatic_focus_shift(double lower_wvl, double higher_wvl) } -void Aberration::plot_spot_diagram(int nrd, double scale, double dot_size) +void Aberration::plot_spot_diagram(int pattern, int nrd, double scale, double dot_size) { renderer_->set_grid_layout(num_fld_,1); @@ -390,40 +367,52 @@ void Aberration::plot_spot_diagram(int nrd, double scale, double dot_size) tracer->set_aperture_check(true); tracer->set_apply_vig(false); - Eigen::Vector2d pupil; + std::vector pupils; + + switch (pattern) { + case Aberration::SpotRayPattern::Grid : + pupils = create_grid_circle(nrd); + break; + case Aberration::SpotRayPattern::Hexapolar : + pupils = create_hexapolar(nrd); + break; + default: + throw "Undefined spot pattern"; + } + std::vector ray_dx, ray_dy; + std::shared_ptr ray; + + const WvlSpec* wvl_spec = opt_sys_->optical_spec()->spectral_region(); for(int fi = 0; fi < num_fld_; fi++){ renderer_->set_current_cell(num_fld_ - fi -1, 0); + Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); + // trace chief ray //pupil << 0.0, 0.0; //Ray chief_ray = tracer->trace_pupil_ray(pupil, fi, ref_wvl_idx_); - Ray chief_ray = opt_sys_->reference_ray(1, fi, ref_wvl_idx_); - double chief_ray_x = chief_ray.back().x(); - double chief_ray_y = chief_ray.back().y(); + std::shared_ptr chief_ray = opt_sys_->reference_ray(1, fi, ref_wvl_idx_); + double chief_ray_x = chief_ray->back()->x(); + double chief_ray_y = chief_ray->back()->y(); for(int wi = 0; wi < num_wvl_; wi++){ + double wvl = wvl_spec->wvl(wi)->value(); ray_dx.clear(); ray_dy.clear(); - for(int pi = 0; pi < nrd; pi++) { - for(int pj = 0; pj < nrd; pj++) { - pupil(0) = -1.0 + (double)pi*2.0/(double)(nrd-1); - pupil(1) = -1.0 + (double)pj*2.0/(double)(nrd-1); - - Ray ray = tracer->trace_pupil_ray(pupil,fi,wi); - - if(ray.status() == RayStatus::PassThrough){ - ray_dx.push_back(ray.back().x() - chief_ray_x); - ray_dy.push_back(ray.back().y() - chief_ray_y); + for(auto& pupil : pupils){ + ray = tracer->trace_pupil_ray(pupil, fld, wvl); - } + if(ray->status() == RayStatus::PassThrough){ + ray_dx.push_back(ray->back()->x() - chief_ray_x); + ray_dy.push_back(ray->back()->y() - chief_ray_y); } } - Rgb color = opt_sys_->optical_spec()->spectral_region()->wvl(wi)->render_color(); - renderer_->draw_dots(ray_dx, ray_dy, color, dot_size); + Rgb color = wvl_spec->wvl(wi)->render_color(); + renderer_->draw_dots(ray_dx, ray_dy, color, dot_size); } renderer_->set_y_axis_range(-scale, scale); @@ -440,3 +429,54 @@ void Aberration::plot_spot_diagram(int nrd, double scale, double dot_size) renderer_->update(); } + + + +std::vector Aberration::create_grid_circle(int nrd) +{ + std::vector pupils; + Eigen::Vector2d pupil; + + for(int pi = 0; pi < nrd; pi++) { + pupil(0) = -1.0 + (double)pi*2.0/(double)(nrd-1); + + for(int pj = 0; pj < nrd; pj++) { + pupil(1) = -1.0 + (double)pj*2.0/(double)(nrd-1); + + //double r = sqrt( pow(pupil(0),2) + pow(pupil(1),2) ); + double r2 = pow(pupil(0),2) + pow(pupil(1),2); + if(r2 <= 1.0){ + pupils.push_back(pupil); + } + } + } + + return pupils; +} + +std::vector Aberration::create_hexapolar(int nrd) +{ + std::vector pupils; + Eigen::Vector2d pupil; + + int half_num_rings = nrd/2; + for (int r = 0; r < nrd/2; r++){ + int num_rays_in_ring = 6*r; + if(num_rays_in_ring == 0){ + pupil(0) = 0.0; + pupil(1) = 0.0; + pupils.push_back(pupil); + continue; + } + + double ang_step = 2*M_PI/(double)num_rays_in_ring; + + for(int ai = 0; ai < num_rays_in_ring; ai++){ + pupil(0) = (double)r * 1.0/(half_num_rings) * cos((double)ai*ang_step); + pupil(1) = (double)r * 1.0/(half_num_rings) * sin((double)ai*ang_step); + pupils.push_back(pupil); + } + } + + return pupils; +} diff --git a/geopter/optical/src/Analysis/layout.cpp b/geopter/optical/src/Analysis/layout.cpp index 91d37e0..d98c506 100644 --- a/geopter/optical/src/Analysis/layout.cpp +++ b/geopter/optical/src/Analysis/layout.cpp @@ -64,7 +64,7 @@ void Layout::draw_reference_rays() int num_flds = opt_sys_->optical_spec()->field_of_view()->field_count(); Rgb color; - Ray ray; + std::shared_ptr ray; for(int fi = 0; fi < num_flds; fi++) { @@ -84,11 +84,11 @@ void Layout::draw_reference_rays() void Layout::draw_fan_rays(int nrd) { - int ref_wvl_idx = opt_sys_->optical_spec()->spectral_region()->reference_index(); + double ref_wvl_val = opt_sys_->optical_spec()->spectral_region()->reference_wvl(); int num_flds = opt_sys_->optical_spec()->field_of_view()->field_count(); Rgb color; - Ray ray; + std::shared_ptr ray; Eigen::Vector2d pupil; SequentialTrace *tracer = new SequentialTrace(opt_sys_); @@ -97,12 +97,14 @@ void Layout::draw_fan_rays(int nrd) for(int fi = 0; fi < num_flds; fi++) { + Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); + color = opt_sys_->optical_spec()->field_of_view()->field(fi)->render_color(); for(int ri = 0; ri < nrd; ri++) { pupil(0) = 0.0; pupil(1) = -1.0 + (double)ri*step; - ray = tracer->trace_pupil_ray(pupil, fi, ref_wvl_idx); + ray = tracer->trace_pupil_ray(pupil, fld, ref_wvl_val); draw_single_ray(ray, color); } @@ -248,13 +250,13 @@ void Layout::draw_surface(Surface* srf, double max_y, const Rgb& color) renderer_->draw_polyline(pts,color); } -void Layout::draw_single_ray(const Ray& ray, const Rgb& color) +void Layout::draw_single_ray(const std::shared_ptr& ray, const Rgb& color) { - for(int i = 1; i < ray.size(); i++) + for(int i = 1; i < ray->size(); i++) { - Eigen::Vector3d pt_to = ray.at(i).intersect_pt; + Eigen::Vector3d pt_to = ray->at(i)->intersect_pt; - Eigen::Vector3d pt_from = ray.at(i-1).intersect_pt; + Eigen::Vector3d pt_from = ray->at(i-1)->intersect_pt; Surface* cur_srf = opt_sys_->optical_assembly()->surface(i); Surface* prev_srf = opt_sys_->optical_assembly()->surface(i-1); diff --git a/geopter/optical/src/Analysis/longitudinal.cpp b/geopter/optical/src/Analysis/longitudinal.cpp index 5863352..82ff69d 100644 --- a/geopter/optical/src/Analysis/longitudinal.cpp +++ b/geopter/optical/src/Analysis/longitudinal.cpp @@ -60,7 +60,7 @@ void Longitudinal::plot_lsa(double scale) PupilCrd pupil; - Ray ray; + std::shared_ptr ray; Field* fld0 = opt_sys_->optical_spec()->field_of_view()->field(0); int ref_wvl_idx = opt_sys_->optical_spec()->spectral_region()->reference_index(); @@ -94,10 +94,10 @@ void Longitudinal::plot_lsa(double scale) //if(ray.status() == RayStatus::PassThrough){ py.push_back(pupil(1)); - int at_image = ray.size()-1; - double y = ray.y(at_image); - double M = ray.M(at_image); - double N = ray.N(at_image); + int at_image = ray->size()-1; + double y = ray->y(at_image); + double M = ray->M(at_image); + double N = ray->N(at_image); lsa.push_back(-y*N/M); y_list.push_back(y); diff --git a/geopter/optical/src/Analysis/ray_fan.cpp b/geopter/optical/src/Analysis/ray_fan.cpp index a298c89..536dde8 100644 --- a/geopter/optical/src/Analysis/ray_fan.cpp +++ b/geopter/optical/src/Analysis/ray_fan.cpp @@ -51,13 +51,14 @@ void RayFan::plot(double scale, int nrd) // trace chief ray Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); - Ray chief_ray = opt_sys_->reference_ray(1,fi,ref_wvl_idx); - double x0 = chief_ray.back().x(); - double y0 = chief_ray.back().y(); + std::shared_ptr chief_ray = opt_sys_->reference_ray(1,fi,ref_wvl_idx); + double x0 = chief_ray->back()->x(); + double y0 = chief_ray->back()->y(); // trace zonal rays for all wavelengths for(int wi = 0; wi < num_wvls; wi++) { + double wvl = opt_sys_->optical_spec()->spectral_region()->wvl(wi)->value(); Rgb render_color = opt_sys_->optical_spec()->spectral_region()->wvl(wi)->render_color(); std::vector px, py; @@ -67,8 +68,8 @@ void RayFan::plot(double scale, int nrd) dx.reserve(nrd); dy.reserve(nrd); - Ray tangential_ray; - Ray sagittal_ray; + std::shared_ptr tangential_ray; + std::shared_ptr sagittal_ray; for(int ri = 0; ri < nrd; ri++) { @@ -76,10 +77,10 @@ void RayFan::plot(double scale, int nrd) pupil(0) = 0.0; pupil(1) = -1.0 + (double)ri*step; - tangential_ray = tracer->trace_pupil_ray(pupil, fi, wi); + tangential_ray = tracer->trace_pupil_ray(pupil, fld, wvl); - if(tangential_ray.status() == RayStatus::PassThrough){ - double y = tangential_ray.back().y(); + if(tangential_ray->status() == RayStatus::PassThrough){ + double y = tangential_ray->back()->y(); vig_pupil = fld->apply_vignetting(pupil); py.push_back(vig_pupil(1)); dy.push_back(y - y0); @@ -89,9 +90,9 @@ void RayFan::plot(double scale, int nrd) pupil(0) = -1.0 + (double)ri*step; pupil(1) = 0.0; - sagittal_ray = tracer->trace_pupil_ray(pupil, fi, wi); - if(sagittal_ray.status() == RayStatus::PassThrough){ - double x = sagittal_ray.back().x(); + sagittal_ray = tracer->trace_pupil_ray(pupil, fld, wvl); + if(sagittal_ray->status() == RayStatus::PassThrough){ + double x = sagittal_ray->back()->x(); vig_pupil = fld->apply_vignetting(pupil); px.push_back(vig_pupil(0)); dx.push_back(x - x0); diff --git a/geopter/optical/src/Assembly/even_polynomial.cpp b/geopter/optical/src/Assembly/even_polynomial.cpp index faec75d..3438ba4 100644 --- a/geopter/optical/src/Assembly/even_polynomial.cpp +++ b/geopter/optical/src/Assembly/even_polynomial.cpp @@ -99,12 +99,12 @@ double EvenPolynomial::sag(double x, double y) const return (z + z_asp); } -double EvenPolynomial::f(Eigen::Vector3d p) const +double EvenPolynomial::f(const Eigen::Vector3d& p) const { return ( p(2) - this->sag(p(0), p(1)) ); } -Eigen::Vector3d EvenPolynomial::df(Eigen::Vector3d p) const +Eigen::Vector3d EvenPolynomial::df(const Eigen::Vector3d& p) const { //sphere + conic contribution double r2 = p(0)*p(0) + p(1)*p(1); @@ -185,7 +185,7 @@ double EvenPolynomial::deriv_2nd(double h) const double z5 = 0.0; for(int i = 0; i < num_coefs_; i++) { double numeric_coef = ( 2*(i+1) + 1 ) * ( 2*(i+1) + 2 ); - double z5_comp = numeric_coef * coefs_[i] * pow(h, 2*(i+1)); + //double z5_comp = numeric_coef * coefs_[i] * pow(h, 2*(i+1)); z5 += numeric_coef * coefs_[i] * pow(h, 2*(i+1)); } diff --git a/geopter/optical/src/Assembly/interface.cpp b/geopter/optical/src/Assembly/interface.cpp index 062c040..10b58c9 100644 --- a/geopter/optical/src/Assembly/interface.cpp +++ b/geopter/optical/src/Assembly/interface.cpp @@ -154,7 +154,7 @@ bool Interface::point_inside(double x, double y) const return true; } -bool Interface::point_inside(Eigen::Vector2d pt) const +bool Interface::point_inside(const Eigen::Vector2d& pt) const { if(clear_aperture_) { return clear_aperture_->point_inside(pt(0), pt(1)); diff --git a/geopter/optical/src/Assembly/spherical.cpp b/geopter/optical/src/Assembly/spherical.cpp index faf9ab3..7b7f074 100644 --- a/geopter/optical/src/Assembly/spherical.cpp +++ b/geopter/optical/src/Assembly/spherical.cpp @@ -15,12 +15,12 @@ Spherical::~Spherical() } -double Spherical::f(Eigen::Vector3d p) const +double Spherical::f(const Eigen::Vector3d& p) const { return p(2) - 0.5*cv_*p.dot(p); } -Eigen::Vector3d Spherical::df(Eigen::Vector3d p) const +Eigen::Vector3d Spherical::df(const Eigen::Vector3d& p) const { Eigen::Vector3d df; df(0) = -cv_*p(0); @@ -52,7 +52,7 @@ double Spherical::sag(double x, double y) const } -void Spherical::intersect(Eigen::Vector3d &pt, double &s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps, int z_dir) +void Spherical::intersect(Eigen::Vector3d &pt, double &s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps, int z_dir) { Eigen::Vector3d p = p0; double ax2 = cv_; diff --git a/geopter/optical/src/Assembly/surface.cpp b/geopter/optical/src/Assembly/surface.cpp index 1659272..571887e 100644 --- a/geopter/optical/src/Assembly/surface.cpp +++ b/geopter/optical/src/Assembly/surface.cpp @@ -29,7 +29,7 @@ void Surface::update() } -void Surface::intersect(Eigen::Vector3d &pt, double &s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps, double z_dir) +void Surface::intersect(Eigen::Vector3d &pt, double &s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps, double z_dir) { profile_->intersect(pt,s,p0,d,eps,z_dir); } diff --git a/geopter/optical/src/Assembly/surface_profile.cpp b/geopter/optical/src/Assembly/surface_profile.cpp index b5e76e3..6ffc547 100644 --- a/geopter/optical/src/Assembly/surface_profile.cpp +++ b/geopter/optical/src/Assembly/surface_profile.cpp @@ -55,17 +55,17 @@ double SurfaceProfile::radius() const return r; } -double SurfaceProfile::f(Eigen::Vector3d p) const +double SurfaceProfile::f(const Eigen::Vector3d& p) const { return 0.0; } -Eigen::Vector3d SurfaceProfile::df(Eigen::Vector3d p) const +Eigen::Vector3d SurfaceProfile::df(const Eigen::Vector3d& p) const { return Eigen::Vector3d::Zero(3); } -Eigen::Vector3d SurfaceProfile::normal(Eigen::Vector3d p) const +Eigen::Vector3d SurfaceProfile::normal(const Eigen::Vector3d& p) const { return df(p).normalized(); } @@ -85,13 +85,13 @@ double SurfaceProfile::deriv_2nd(double h) const return 0.0; } -void SurfaceProfile::intersect(Eigen::Vector3d& pt, double& s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps, int z_dir) +void SurfaceProfile::intersect(Eigen::Vector3d& pt, double& s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps, int z_dir) { intersect_spencer(pt, s, p0, d, eps, z_dir); } -void SurfaceProfile::intersect_spencer(Eigen::Vector3d& pt, double& s, Eigen::Vector3d p0, Eigen::Vector3d d, double eps, int z_dir) +void SurfaceProfile::intersect_spencer(Eigen::Vector3d& pt, double& s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps, int z_dir) { Eigen::Vector3d p = p0; double s1 = -f(p)/d.dot(df(p)); diff --git a/geopter/optical/src/CMakeLists.txt b/geopter/optical/src/CMakeLists.txt index 2ed25c7..8d21e3f 100644 --- a/geopter/optical/src/CMakeLists.txt +++ b/geopter/optical/src/CMakeLists.txt @@ -6,8 +6,8 @@ SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) set(OPTICAL_SRCS Analysis/aberration.cpp - Analysis/ray_fan.cpp - Analysis/longitudinal.cpp + #Analysis/ray_fan.cpp + #Analysis/longitudinal.cpp Analysis/layout.cpp Assembly/optical_assembly.cpp diff --git a/geopter/optical/src/Paraxial/paraxial_ray.cpp b/geopter/optical/src/Paraxial/paraxial_ray.cpp index 318a525..b9b86b7 100644 --- a/geopter/optical/src/Paraxial/paraxial_ray.cpp +++ b/geopter/optical/src/Paraxial/paraxial_ray.cpp @@ -27,7 +27,7 @@ int ParaxialRay::size() const return (int)par_ray_at_srfs_.size(); } -void ParaxialRay::prepend(ParaxialRayAtSurface par_ray_at_srf) +void ParaxialRay::prepend(std::shared_ptr par_ray_at_srf) { auto itr = par_ray_at_srfs_.begin(); par_ray_at_srfs_.insert(itr, par_ray_at_srf); @@ -35,19 +35,20 @@ void ParaxialRay::prepend(ParaxialRayAtSurface par_ray_at_srf) void ParaxialRay::prepend(double ht, double slp, double aoi, double n) { + auto prx_ray_at_srf = std::make_shared(ht,slp,aoi,n); auto itr = par_ray_at_srfs_.begin(); - par_ray_at_srfs_.insert( itr, ParaxialRayAtSurface(ht,slp,aoi,n) ); + par_ray_at_srfs_.insert( itr, prx_ray_at_srf); } -void ParaxialRay::append(ParaxialRayAtSurface par_ray_at_srf) +void ParaxialRay::append(std::shared_ptr par_ray_at_srf) { par_ray_at_srfs_.push_back(par_ray_at_srf); } void ParaxialRay::append(double ht, double slp, double aoi, double n) { - ParaxialRayAtSurface par_ray_at_srf(ht, slp, aoi, n); - par_ray_at_srfs_.push_back(par_ray_at_srf); + auto prx_ray_at_srf = std::make_shared(ht, slp, aoi, n); + par_ray_at_srfs_.push_back(prx_ray_at_srf); } std::string ParaxialRay::name() const @@ -60,17 +61,17 @@ void ParaxialRay::set_name(std::string name) name_ = name; } -ParaxialRayAtSurface ParaxialRay::at(int i) const +std::shared_ptr ParaxialRay::at(int i) const { return par_ray_at_srfs_[i]; } -ParaxialRayAtSurface ParaxialRay::front() const +std::shared_ptr ParaxialRay::front() const { return par_ray_at_srfs_.front(); } -ParaxialRayAtSurface ParaxialRay::back() const +std::shared_ptr ParaxialRay::back() const { return par_ray_at_srfs_.back(); } @@ -82,7 +83,7 @@ void ParaxialRay::print() const std::cout << oss.str() << std::endl; } -ParaxialRayAtSurface ParaxialRay::at_before_image() const +std::shared_ptr ParaxialRay::at_before_image() const { int img_idx = par_ray_at_srfs_.size() - 1; return par_ray_at_srfs_[img_idx - 1]; @@ -109,15 +110,15 @@ void ParaxialRay::print(std::ostringstream& oss) const for(int i = 0; i < num_srf; i++) { oss << std::setw(idx_w) << std::right << i; - oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i].ht; - oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i].slp; + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->ht; + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->slp; if( i == num_srf-1){ - oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i].slp; + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->slp; }else{ - oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i+1].n*par_ray_at_srfs_[i].slp; + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i+1]->n*par_ray_at_srfs_[i]->slp; } - oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i].aoi; - oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i].n * par_ray_at_srfs_[i].aoi; + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->aoi; + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->n * par_ray_at_srfs_[i]->aoi; oss << std::endl; } oss << std::endl; diff --git a/geopter/optical/src/Paraxial/paraxial_trace.cpp b/geopter/optical/src/Paraxial/paraxial_trace.cpp index 86aceb1..1669fa5 100644 --- a/geopter/optical/src/Paraxial/paraxial_trace.cpp +++ b/geopter/optical/src/Paraxial/paraxial_trace.cpp @@ -23,9 +23,9 @@ ParaxialTrace::~ParaxialTrace() opt_sys_ = nullptr; } -ParaxialRay ParaxialTrace::trace_paraxial_ray_from_object(double y0, double u0, double wvl) +std::shared_ptr ParaxialTrace::trace_paraxial_ray_from_object(double y0, double u0, double wvl) { - ParaxialRay par_ray; + auto par_ray = std::make_shared(); int img = opt_sys_->optical_assembly()->image_index(); ParaxialPath par_path = paraxial_path(0, img, wvl); @@ -41,7 +41,7 @@ ParaxialRay ParaxialTrace::trace_paraxial_ray_from_object(double y0, double u0, double u0_prime = u0; double i0 = y0*c0 + u0; double n0 = par_path.at(0).n; - par_ray.append(y0, u0_prime, i0, n0); + par_ray->append(y0, u0_prime, i0, n0); if(path_size == 1) { return par_ray; @@ -65,7 +65,7 @@ ParaxialRay ParaxialTrace::trace_paraxial_ray_from_object(double y0, double u0, i_prime = i*(n/n_prime); u_prime = i_prime - y*c; - par_ray.append(y, u_prime, i, n); + par_ray->append(y, u_prime, i, n); // transfer to next d = par_path.at(k).t; @@ -77,7 +77,7 @@ ParaxialRay ParaxialTrace::trace_paraxial_ray_from_object(double y0, double u0, return par_ray; } -ParaxialRay ParaxialTrace::trace_paraxial_ray_parallel_to_axis_at_s1() +std::shared_ptr ParaxialTrace::trace_paraxial_ray_parallel_to_axis_at_s1() { const double ref_wvl = opt_sys_->optical_spec()->spectral_region()->reference_wvl(); const int img = opt_sys_->optical_assembly()->image_index(); @@ -86,7 +86,8 @@ ParaxialRay ParaxialTrace::trace_paraxial_ray_parallel_to_axis_at_s1() constexpr double y1 = 1.0; constexpr double u1 = 0.0; - ParaxialRay par_ray; + + auto par_ray = std::make_shared(); // at object double t0 = par_path.at(0).t; @@ -103,7 +104,7 @@ ParaxialRay ParaxialTrace::trace_paraxial_ray_parallel_to_axis_at_s1() } -ParaxialRay ParaxialTrace::trace_paraxial_ray_with_slope_at_s1() +std::shared_ptr ParaxialTrace::trace_paraxial_ray_with_slope_at_s1() { const int img = opt_sys_->optical_assembly()->image_index(); @@ -114,7 +115,7 @@ ParaxialRay ParaxialTrace::trace_paraxial_ray_with_slope_at_s1() ParaxialPath par_path = paraxial_path(0, img, ref_wvl); - ParaxialRay p_ray_bar; + auto p_ray_bar = std::make_shared(); double ybar1 = 0.0; double ubar1 = uq0; @@ -203,18 +204,18 @@ void ParaxialTrace::compute_starting_data() auto p_ray = trace_paraxial_ray_parallel_to_axis_at_s1(); auto q_ray = trace_paraxial_ray_with_slope_at_s1(); - double ak1 = p_ray.back().ht; - double bk1 = q_ray.back().ht; - double ck1 = n_k*p_ray.back().slp; - double dk1 = n_k*q_ray.back().slp; + double ak1 = p_ray->back()->ht; + double bk1 = q_ray->back()->ht; + double ck1 = n_k*p_ray->back()->slp; + double dk1 = n_k*q_ray->back()->slp; int stop = opt_sys_->optical_assembly()->stop_index(); double n_s = opt_sys_->optical_assembly()->gap(stop)->material()->rindex(ref_wvl); - double as1 = p_ray.at(stop).ht; - double bs1 = q_ray.at(stop).ht; - double cs1 = n_s*p_ray.at(stop).slp; - double ds1 = n_s*q_ray.at(stop).slp; + double as1 = p_ray->at(stop)->ht; + double bs1 = q_ray->at(stop)->ht; + double cs1 = n_s*p_ray->at(stop)->slp; + double ds1 = n_s*q_ray->at(stop)->slp; double ybar1 = -bs1; double ubar1 = as1; diff --git a/geopter/optical/src/Sequential/ray.cpp b/geopter/optical/src/Sequential/ray.cpp index b97045e..f8e1a1e 100644 --- a/geopter/optical/src/Sequential/ray.cpp +++ b/geopter/optical/src/Sequential/ray.cpp @@ -14,6 +14,9 @@ Ray::Ray() Ray::~Ray() { + for(auto &r : ray_at_srfs_){ + r.reset(); + } ray_at_srfs_.clear(); } @@ -22,14 +25,19 @@ int Ray::size() const return (int)ray_at_srfs_.size(); } -void Ray::prepend(RayAtSurface ray_at_srf) +void Ray::prepend(std::shared_ptr ray_at_srf) { ray_at_srfs_.insert(ray_at_srfs_.begin(), ray_at_srf); + ray_at_srfs_.front()->before = ray_at_srfs_[1].get(); } -void Ray::append(RayAtSurface ray_at_srf) +void Ray::append(std::shared_ptr ray_at_srf) { ray_at_srfs_.push_back(ray_at_srf); + int len = ray_at_srfs_.size(); + if(len > 1){ + ray_at_srfs_.back()->before = ray_at_srfs_[len-1-1].get(); + } } @@ -38,25 +46,25 @@ void Ray::set_wvl(double wvl) wvl_ = wvl; } -RayAtSurface Ray::at(int i) const +RayAtSurface* Ray::at(int i) const { assert(i >= 0); if( i < (int)ray_at_srfs_.size() ) { - return ray_at_srfs_[i]; + return ray_at_srfs_[i].get(); } else { throw "Out of index"; } } -RayAtSurface Ray::front() const +RayAtSurface* Ray::front() const { - return ray_at_srfs_.front(); + return ray_at_srfs_.front().get(); } -RayAtSurface Ray::back() const +RayAtSurface* Ray::back() const { - return ray_at_srfs_.back(); + return ray_at_srfs_.back().get(); } void Ray::set_status(int s) @@ -71,32 +79,32 @@ int Ray::status() const double Ray::x(int i) const { - return ray_at_srfs_[i].intersect_pt(0); + return ray_at_srfs_[i]->intersect_pt(0); } double Ray::y(int i) const { - return ray_at_srfs_[i].intersect_pt(1); + return ray_at_srfs_[i]->intersect_pt(1); } double Ray::z(int i) const { - return ray_at_srfs_[i].intersect_pt(2); + return ray_at_srfs_[i]->intersect_pt(2); } double Ray::L(int i) const { - return ray_at_srfs_[i].after_dir(0); + return ray_at_srfs_[i]->after_dir(0); } double Ray::M(int i) const { - return ray_at_srfs_[i].after_dir(1); + return ray_at_srfs_[i]->after_dir(1); } double Ray::N(int i) const { - return ray_at_srfs_[i].after_dir(2); + return ray_at_srfs_[i]->after_dir(2); } double Ray::aoi(int i) const @@ -105,11 +113,11 @@ double Ray::aoi(int i) const Eigen::Vector3d normal; if(i > 0) { - inc_dir = ray_at_srfs_[i-1].after_dir; - normal = ray_at_srfs_[i].normal; + inc_dir = ray_at_srfs_[i-1]->after_dir; + normal = ray_at_srfs_[i]->normal; }else{ - inc_dir = ray_at_srfs_[0].after_dir; - normal = ray_at_srfs_[0].normal; + inc_dir = ray_at_srfs_[0]->after_dir; + normal = ray_at_srfs_[0]->normal; } // We need signed value @@ -121,8 +129,8 @@ double Ray::aoi(int i) const double Ray::aor(int i) { - Eigen::Vector3d after_dir = ray_at_srfs_[i].after_dir; - Eigen::Vector3d normal = ray_at_srfs_[i].normal; + Eigen::Vector3d after_dir = ray_at_srfs_[i]->after_dir; + Eigen::Vector3d normal = ray_at_srfs_[i]->normal; double tanU1 = after_dir(1)/after_dir(2); double tanU2 = normal(1)/normal(2); double tanI_prime = (tanU1 - tanU2)/(1.0 + tanU1*tanU2); @@ -150,13 +158,13 @@ void Ray::print(std::ostringstream& oss) for(int si = 0; si < num_srfs; si++) { - Eigen::Vector3d intercept = ray_at_srfs_[si].intersect_pt; + Eigen::Vector3d intercept = ray_at_srfs_[si]->intersect_pt; oss << std::setw(idx_w) << std::right << si; oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << intercept(0); oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << intercept(1); oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << intercept(2); - Eigen::Vector3d after_dir = ray_at_srfs_[si].after_dir; + Eigen::Vector3d after_dir = ray_at_srfs_[si]->after_dir; oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << after_dir(0); oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << after_dir(1); oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << after_dir(2); diff --git a/geopter/optical/src/Sequential/sequential_trace.cpp b/geopter/optical/src/Sequential/sequential_trace.cpp index a9cd705..32e0109 100644 --- a/geopter/optical/src/Sequential/sequential_trace.cpp +++ b/geopter/optical/src/Sequential/sequential_trace.cpp @@ -22,6 +22,11 @@ SequentialTrace::SequentialTrace(OpticalSystem* sys): ref_wvl_idx_ = opt_sys_->optical_spec()->spectral_region()->reference_index(); ref_wvl_val_ = opt_sys_->optical_spec()->spectral_region()->reference_wvl(); + num_srfs_ = opt_sys_->optical_assembly()->surface_count(); + num_gaps_ = opt_sys_->optical_assembly()->gap_count(); + image_index_ = opt_sys_->optical_assembly()->image_index(); + object_distance_ = opt_sys_->optical_assembly()->gap(0)->thi(); + do_aperture_check_ = false; do_apply_vig_ = true; } @@ -31,15 +36,17 @@ SequentialTrace::~SequentialTrace() opt_sys_ = nullptr; } -Ray SequentialTrace::trace_pupil_ray(Eigen::Vector2d pupil_crd, int fi, int wi) +/* +std::shared_ptr SequentialTrace::trace_pupil_ray(Eigen::Vector2d pupil_crd, int fi, int wi) { Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); double wvl = opt_sys_->optical_spec()->spectral_region()->wvl(wi)->value(); return trace_pupil_ray(pupil_crd, fld, wvl); } +*/ -Ray SequentialTrace::trace_pupil_ray(Eigen::Vector2d pupil_crd, const Field *fld, double wvl) +std::shared_ptr SequentialTrace::trace_pupil_ray(const Eigen::Vector2d& pupil_crd, const Field *fld, double wvl) { Eigen::Vector2d vig_pupil = pupil_crd; @@ -49,10 +56,7 @@ Ray SequentialTrace::trace_pupil_ray(Eigen::Vector2d pupil_crd, const Field *fld auto fod = opt_sys_->first_order_data(); double eprad = fod.enp_radius; - PupilCrd aim_pt; - //aim_pt(0) = 0.0; - //aim_pt(1) = 0.0; - aim_pt = fld->aim_pt(); + PupilCrd aim_pt = fld->aim_pt(); Eigen::Vector3d pt0 = this->object_coord(fld); @@ -67,20 +71,20 @@ Ray SequentialTrace::trace_pupil_ray(Eigen::Vector2d pupil_crd, const Field *fld //dir0 = dir0/length; dir0.normalize(); - int img = opt_sys_->optical_assembly()->image_index(); - SequentialPath path = sequential_path(0, img, wvl); + //int img = opt_sys_->optical_assembly()->image_index(); + SequentialPath path = sequential_path(0, image_index_, wvl); - Ray ray = trace_ray_throughout_path(path, pt0, dir0); - ray.set_wvl(wvl); + std::shared_ptr ray = trace_ray_throughout_path(path, pt0, dir0); + ray->set_wvl(wvl); return ray; } -Ray SequentialTrace::trace_ray_throughout_path(const SequentialPath& seq_path, Eigen::Vector3d pt0, Eigen::Vector3d dir0) +std::shared_ptr SequentialTrace::trace_ray_throughout_path(const SequentialPath& seq_path, const Eigen::Vector3d& pt0, const Eigen::Vector3d& dir0) { - Ray ray; - ray.set_status(RayStatus::PassThrough); + std::shared_ptr ray = std::make_shared(); + ray->set_status(RayStatus::PassThrough); if(seq_path.size() == 0) { return ray; // empty ray @@ -100,14 +104,14 @@ Ray SequentialTrace::trace_ray_throughout_path(const SequentialPath& seq_path, E // first surface - RayAtSurface ray_at_srf; - ray_at_srf.intersect_pt = pt0; - ray_at_srf.distance_from_before = 0.0; - ray_at_srf.after_dir = dir0; - ray_at_srf.normal = seq_path.at(0).srf->normal(pt0); - ray_at_srf.optical_path_length = 0.0; + auto ray_at_srf = std::make_shared(); + ray_at_srf->intersect_pt = pt0; + ray_at_srf->distance_from_before = 0.0; + ray_at_srf->after_dir = dir0; + ray_at_srf->normal = seq_path.at(0).srf->normal(pt0); + ray_at_srf->optical_path_length = 0.0; - ray.append(ray_at_srf); + ray->append(ray_at_srf); // trace ray till the image int path_size = seq_path.size(); @@ -138,17 +142,18 @@ Ray SequentialTrace::trace_ray_throughout_path(const SequentialPath& seq_path, E n_out = seq_path.at(i).rind; after_dir = bend(before_dir, srf_normal, n_in, n_out); - ray_at_srf.intersect_pt = intersect_pt; - ray_at_srf.after_dir = after_dir.normalized(); - ray_at_srf.distance_from_before = distance_from_before; - ray_at_srf.optical_path_length = n_in*distance_from_before; - ray_at_srf.normal = srf_normal; + ray_at_srf = std::make_shared(); + ray_at_srf->intersect_pt = intersect_pt; + ray_at_srf->after_dir = after_dir.normalized(); + ray_at_srf->distance_from_before = distance_from_before; + ray_at_srf->optical_path_length = n_in*distance_from_before; + ray_at_srf->normal = srf_normal; - ray.append(ray_at_srf); + ray->append(ray_at_srf); if(do_aperture_check_) { if( !cur_srf->point_inside(intersect_pt(0),intersect_pt(1)) ){ - ray.set_status(RayStatus::Blocked); + ray->set_status(RayStatus::Blocked); return ray; break; } @@ -161,10 +166,10 @@ Ray SequentialTrace::trace_ray_throughout_path(const SequentialPath& seq_path, E } } catch (TraceMissedSurfaceError& ray_miss) { - ray.set_status(RayStatus::MissedSurface); + ray->set_status(RayStatus::MissedSurface); return ray; } catch (TraceTIRError& ray_tir){ - ray.set_status(RayStatus::TotalReflection); + ray->set_status(RayStatus::TotalReflection); return ray; } @@ -189,12 +194,12 @@ Eigen::Vector2d SequentialTrace::trace_coddington(const Field *fld, double wvl) ParaxialTrace *parax_tracer = new ParaxialTrace(opt_sys_); - double y0 = opt_sys_->axial_ray(ref_wvl_idx_).at(0).ht; - double u0 = opt_sys_->axial_ray(ref_wvl_idx_).at(0).slp; - ParaxialRay ax_ray = parax_tracer->trace_paraxial_ray_from_object(y0, u0, wvl); + double y0 = opt_sys_->axial_ray(ref_wvl_idx_)->at(0)->ht; + double u0 = opt_sys_->axial_ray(ref_wvl_idx_)->at(0)->slp; + auto ax_ray = parax_tracer->trace_paraxial_ray_from_object(y0, u0, wvl); - double y = ax_ray.at(last_surf).ht; - double u_prime = ax_ray.at(img_srf).slp; + double y = ax_ray->at(last_surf)->ht; + double u_prime = ax_ray->at(img_srf)->slp; double l_prime = -y/u_prime; delete parax_tracer; @@ -203,7 +208,7 @@ Eigen::Vector2d SequentialTrace::trace_coddington(const Field *fld, double wvl) // Here, we have l_prime // off axis - Ray ray = trace_pupil_ray(PupilCrd({0.0, 0.0}), fld, wvl); + std::shared_ptr ray = trace_pupil_ray(PupilCrd({0.0, 0.0}), fld, wvl); double s_before, t_before, s_after, t_after; double n_before, n_after; @@ -216,31 +221,32 @@ Eigen::Vector2d SequentialTrace::trace_coddington(const Field *fld, double wvl) s_before = std::numeric_limits::infinity(); t_before = std::numeric_limits::infinity(); }else{ - auto dir0 = ray.at(0).after_dir; + auto dir0 = ray->at(0)->after_dir; double cosUpr = dir0(2); - double B = opt_sys_->optical_assembly()->gap(0)->thi(); - double Zpr = ray.z(1); + //double B = opt_sys_->optical_assembly()->gap(0)->thi(); + double B = object_distance_; + double Zpr = ray->z(1); s_before = (B - Zpr)/cosUpr; t_before = s_before; } - for(int i = 1; i < ray.size()-1; i++) { + for(int i = 1; i < ray->size()-1; i++) { Surface* surf = path.at(i).srf; n_before = path.at(i-1).rind; n_after = path.at(i).rind; - cosI = cos(ray.aoi(i)); - cosI_prime = cos(ray.aor(i)); + cosI = cos(ray->aoi(i)); + cosI_prime = cos(ray->aor(i)); //sinI = sqrt(1.0 - cosI*cosI); - sinI = sin(ray.aoi(i)); + sinI = sin(ray->aoi(i)); sinI_prime = sinI * n_before/n_after; //cosI_prime = sqrt(1.0 - sinI_prime*sinI_prime); - cosU = ray.at(i-1).after_dir(2);//ray.at(i-1).after_dir.norm(); + cosU = ray->at(i-1)->after_dir(2);//ray.at(i-1).after_dir.norm(); sinU = sqrt(1.0 - cosU*cosU); - cosU_prime = ray.at(i).after_dir(2);//ray.at(i).after_dir.norm(); + cosU_prime = ray->at(i)->after_dir(2);//ray.at(i).after_dir.norm(); if(surf->profile()->name() == "SPH") { double c = surf->profile()->cv(); @@ -248,7 +254,7 @@ Eigen::Vector2d SequentialTrace::trace_coddington(const Field *fld, double wvl) obl_pwr_t = obl_pwr_s; }else{ // aspherical double sinI_U = sinI*cosU - cosI*sinU; - double y = ray.y(i); + double y = ray->y(i); double cs = sinI_U/y; obl_pwr_s = cs*(n_after*cosI_prime - n_before*cosI); @@ -259,8 +265,8 @@ Eigen::Vector2d SequentialTrace::trace_coddington(const Field *fld, double wvl) } - double z1 = ray.z(i); - double z2 = ray.z(i+1); + double z1 = ray->z(i); + double z2 = ray->z(i+1); double d = path.at(i).d; double D = (d + z2 - z1)/cosU_prime; @@ -275,7 +281,7 @@ Eigen::Vector2d SequentialTrace::trace_coddington(const Field *fld, double wvl) } // Closing Equation - double z = ray.z(ray.size()-1-1); + double z = ray->z(ray->size()-1-1); double zs = s_after*cosU_prime + z - l_prime; double zt = t_after*cosU_prime + z - l_prime; @@ -296,23 +302,23 @@ SequentialPath SequentialTrace::overall_sequential_path(double wvl) SequentialPath SequentialTrace::sequential_path(int start, int end, double wvl) { - int num_srfs = opt_sys_->optical_assembly()->surface_count(); - int num_gaps = opt_sys_->optical_assembly()->gap_count(); + //int num_srfs = opt_sys_->optical_assembly()->surface_count(); + //int num_gaps = opt_sys_->optical_assembly()->gap_count(); - assert(end <= num_srfs); + assert(end <= num_srfs_); SequentialPath path; SequentialPathComponent path_comp; for(int i = start; i <= end; i++) { - if ( i < num_srfs ) { + if ( i < num_srfs_ ) { path_comp.srf = opt_sys_->optical_assembly()->surface(i); }else{ path_comp.srf = nullptr; } - if( i < num_gaps ) { + if( i < num_gaps_ ) { path_comp.d = opt_sys_->optical_assembly()->gap(i)->thi(); path_comp.rind = opt_sys_->optical_assembly()->gap(i)->material()->rindex(wvl); }else { @@ -326,6 +332,7 @@ SequentialPath SequentialTrace::sequential_path(int start, int end, double wvl) return path; } +/* Eigen::Vector2d SequentialTrace::aim_chief_ray(int fi, int wi) { Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); @@ -333,6 +340,7 @@ Eigen::Vector2d SequentialTrace::aim_chief_ray(int fi, int wi) return aim_chief_ray(fld, wvl); } +*/ Eigen::Vector2d SequentialTrace::aim_chief_ray(const Field *fld, double wvl) { @@ -344,6 +352,7 @@ Eigen::Vector2d SequentialTrace::aim_chief_ray(const Field *fld, double wvl) return aim_pt; } +/* Eigen::Vector2d SequentialTrace::search_aim_point(int srf_idx, Eigen::Vector2d xy_target, int fi, int wi) { Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); @@ -351,8 +360,9 @@ Eigen::Vector2d SequentialTrace::search_aim_point(int srf_idx, Eigen::Vector2d x return search_aim_point(srf_idx, xy_target, fld, wvl); } +*/ -Eigen::Vector2d SequentialTrace::search_aim_point(int srf_idx, Eigen::Vector2d xy_target, const Field *fld, double wvl) +Eigen::Vector2d SequentialTrace::search_aim_point(int srf_idx, const Eigen::Vector2d& xy_target, const Field *fld, double wvl) { assert(srf_idx > 0); @@ -415,7 +425,7 @@ double SequentialTrace::y_stop_coordinate(double y1, int ifcx, Eigen::Vector3d p //std::vector ray; std::vector ray_at_srfs; - Ray ray_trace_result; + std::shared_ptr ray_trace_result; try { //RayTraceResult ray_tr_rslt = RayTrace::trace(seq_model, pt0, dir0, wvl); ray_trace_result = trace_ray_throughout_path(overall_sequential_path(wvl), pt0, dir0); @@ -435,7 +445,7 @@ double SequentialTrace::y_stop_coordinate(double y1, int ifcx, Eigen::Vector3d p } - double y_ray = ray_trace_result.at(ifcx).intersect_pt(1); + double y_ray = ray_trace_result->at(ifcx)->intersect_pt(1); //double y_ray = ray_at_srfs[ifcx].intersect_pt(1); return (y_ray - y_target); @@ -458,11 +468,13 @@ Eigen::Vector3d SequentialTrace::bend(Eigen::Vector3d d_in, Eigen::Vector3d norm } } +/* Eigen::Vector3d SequentialTrace::object_coord(int fi) { Field* fld = opt_sys_->optical_spec()->field_of_view()->field(fi); return object_coord(fld); } +*/ Eigen::Vector3d SequentialTrace::object_coord(const Field* fld) { @@ -537,9 +549,9 @@ double SequentialTrace::compute_vignetting_factor_for_pupil(Eigen::Vector2d full vig_pupil(0) = full_pupil(0)*(1.0 - a); vig_pupil(1) = full_pupil(1)*(1.0 - a); - Ray ray_full_marginal = trace_pupil_ray(vig_pupil, &fld, ref_wvl_val_); + std::shared_ptr ray_full_marginal = trace_pupil_ray(vig_pupil, &fld, ref_wvl_val_); - if(ray_full_marginal.status() == RayStatus::PassThrough){ // no vignetting + if(ray_full_marginal->status() == RayStatus::PassThrough){ // no vignetting return a; // 0.0 } @@ -555,9 +567,9 @@ double SequentialTrace::compute_vignetting_factor_for_pupil(Eigen::Vector2d full m = (a+b)/2.0; vig_pupil(0) = full_pupil(0)*(1.0 - m); vig_pupil(1) = full_pupil(1)*(1.0 - m); - Ray ray_m = trace_pupil_ray(vig_pupil, &fld, ref_wvl_val_); + std::shared_ptr ray_m = trace_pupil_ray(vig_pupil, &fld, ref_wvl_val_); - if(ray_m.status() == RayStatus::PassThrough){ + if(ray_m->status() == RayStatus::PassThrough){ b = m; }else{ a = m; diff --git a/geopter/optical/src/Spec/field.cpp b/geopter/optical/src/Spec/field.cpp index dde41b7..ab9ff6c 100644 --- a/geopter/optical/src/Spec/field.cpp +++ b/geopter/optical/src/Spec/field.cpp @@ -114,25 +114,16 @@ PupilCrd Field::apply_vignetting(PupilCrd pupil) const { PupilCrd vig_pupil = pupil; - if(pupil(0) < 0.0){ - if(vlx_ != 0.0){ - vig_pupil(0) *= (1.0 - vlx_); - } - }else{ - if(vux_ != 0.0){ - vig_pupil(0) *= (1.0 - vux_); - } + if(pupil(0) < 0.0){ + vig_pupil(0) *= (1.0 - vlx_); + }else{ + vig_pupil(0) *= (1.0 - vux_); } - if(pupil(1) < 0.0){ - if(vly_ != 0.0){ - vig_pupil(1) *= (1.0 - vly_); - } - } - else{ - if(vuy_ != 0.0){ - vig_pupil(1) *= (1.0 - vuy_); - } + if(pupil(1) < 0.0){ + vig_pupil(1) *= (1.0 - vly_); + }else{ + vig_pupil(1) *= (1.0 - vuy_); } return vig_pupil; diff --git a/geopter/optical/src/System/optical_system.cpp b/geopter/optical/src/System/optical_system.cpp index ebeef59..3fff5c3 100644 --- a/geopter/optical/src/System/optical_system.cpp +++ b/geopter/optical/src/System/optical_system.cpp @@ -93,12 +93,12 @@ void OpticalSystem::add_surface_and_gap(double r, double t, std::string mat_name opt_assembly_->add_surface_and_gap(s, g); } -ParaxialRay OpticalSystem::axial_ray(int wi) const +std::shared_ptr OpticalSystem::axial_ray(int wi) const { return ax_rays_[wi]; } -ParaxialRay OpticalSystem::principle_ray(int wi) const +std::shared_ptr OpticalSystem::principle_ray(int wi) const { return pr_rays_[wi]; } @@ -108,7 +108,7 @@ FirstOrderData OpticalSystem::first_order_data() const return fod_; } -Ray OpticalSystem::reference_ray(int ri, int fi, int wi) const +std::shared_ptr OpticalSystem::reference_ray(int ri, int fi, int wi) const { switch (ri) { @@ -138,12 +138,13 @@ void OpticalSystem::update_aim_pt() if(opt_assembly_->surface_count() > 2) { int field_count = opt_spec_->field_of_view()->field_count(); - int ref_wi = opt_spec_->spectral_region()->reference_index(); + double ref_wvl = opt_spec_->spectral_region()->reference_wvl(); SequentialTrace *tracer = new SequentialTrace(this); for(int fi = 0; fi < field_count; fi++){ - auto aim_pt = tracer->aim_chief_ray(fi, ref_wi); + Field* fld = opt_spec_->field_of_view()->field(fi); + auto aim_pt = tracer->aim_chief_ray(fld, ref_wvl); opt_spec_->field_of_view()->field(fi)->set_aim_pt(aim_pt); } @@ -164,7 +165,7 @@ void OpticalSystem::update_semi_diameters() } // update semi diameter - Ray chief_ray, mer_upper_ray, mer_lower_ray, sag_upper_ray, sag_lower_ray; + std::shared_ptr chief_ray, mer_upper_ray, mer_lower_ray, sag_upper_ray, sag_lower_ray; for(int fi = 0; fi < num_fld; fi++) { @@ -176,11 +177,11 @@ void OpticalSystem::update_semi_diameters() for(int si = 0; si < num_srf; si++) { - double chief_ray_ht = chief_ray.at(si).height(); - double mer_upper_ray_ht = mer_upper_ray.at(si).height(); - double mer_lower_ray_ht = mer_lower_ray.at(si).height(); - double sag_upper_ray_ht = sag_upper_ray.at(si).height(); - double sag_lower_ray_ht = sag_lower_ray.at(si).height(); + double chief_ray_ht = chief_ray->at(si)->height(); + double mer_upper_ray_ht = mer_upper_ray->at(si)->height(); + double mer_lower_ray_ht = mer_lower_ray->at(si)->height(); + double sag_upper_ray_ht = sag_upper_ray->at(si)->height(); + double sag_lower_ray_ht = sag_lower_ray->at(si)->height(); double ray_ht_for_cur_fld = std::max({chief_ray_ht, mer_upper_ray_ht, mer_lower_ray_ht, sag_upper_ray_ht, sag_lower_ray_ht}); @@ -226,11 +227,11 @@ void OpticalSystem::update_paraxial_data() for(int wi = 0; wi < num_wvls; wi++) { double wvl = opt_spec_->spectral_region()->wvl(wi)->value(); - ParaxialRay ax_ray_for_wvl = tracer->trace_paraxial_ray_from_object(y0, u0, wvl); - ax_ray_for_wvl.set_name("axial ray traced with W" + std::to_string(wi)); + auto ax_ray_for_wvl = tracer->trace_paraxial_ray_from_object(y0, u0, wvl); + ax_ray_for_wvl->set_name("axial ray traced with W" + std::to_string(wi)); - ParaxialRay pr_ray_for_wvl = tracer->trace_paraxial_ray_from_object(ybar0, ubar0, wvl); - pr_ray_for_wvl.set_name("principle ray traced with W" + std::to_string(wi)); + auto pr_ray_for_wvl = tracer->trace_paraxial_ray_from_object(ybar0, ubar0, wvl); + pr_ray_for_wvl->set_name("principle ray traced with W" + std::to_string(wi)); ax_rays_.push_back(ax_ray_for_wvl); pr_rays_.push_back(pr_ray_for_wvl); @@ -244,27 +245,27 @@ void OpticalSystem::update_paraxial_data() double n_0 = opt_assembly_->gap(0)->material()->rindex(ref_wvl); double n_k = opt_assembly_->image_space_gap()->material()->rindex(ref_wvl); - ParaxialRay ax_ray = ax_rays_[ref_wi]; - ParaxialRay pr_ray = pr_rays_[ref_wi]; + auto ax_ray = ax_rays_[ref_wi]; + auto pr_ray = pr_rays_[ref_wi]; - ParaxialRay p_ray = p_ray_; - ParaxialRay q_ray = q_ray_; + std::shared_ptr p_ray = p_ray_; + std::shared_ptr q_ray = q_ray_; - double ak1 = p_ray.at(img).ht; - double bk1 = q_ray.at(img).ht; - double ck1 = n_k*p_ray.at(img).slp; - double dk1 = n_k*q_ray.at(img).slp; + double ak1 = p_ray->at(img)->ht; + double bk1 = q_ray->at(img)->ht; + double ck1 = n_k*p_ray->at(img)->slp; + double dk1 = n_k*q_ray->at(img)->slp; // fill in the content of first order data - fod_.opt_inv = n_0 * ( ax_ray.at(1).ht*pr_ray.at(0).slp - pr_ray.at(1).ht*ax_ray.at(0).slp ); + fod_.opt_inv = n_0 * ( ax_ray->at(1)->ht*pr_ray->at(0)->slp - pr_ray->at(1)->ht*ax_ray->at(0)->slp ); fod_.obj_dist = opt_assembly_->gap(0)->thi(); fod_.img_dist = opt_assembly_->image_space_gap()->thi(); fod_.efl = -1.0/ck1; fod_.pp1 = (dk1 - 1.0)*(n_0/ck1); - fod_.ppk = (p_ray.at(img-1).ht - 1.0)*(n_k/ck1); + fod_.ppk = (p_ray->at(img-1)->ht - 1.0)*(n_k/ck1); fod_.ffl = fod_.pp1 - fod_.efl; fod_.bfl = fod_.efl - fod_.ppk; - fod_.fno = -1.0/(2.0*n_k*ax_ray.at(img).slp); + fod_.fno = -1.0/(2.0*n_k*ax_ray->at(img)->slp); fod_.m = ak1 + ck1*fod_.img_dist/n_k; fod_.red = dk1 + ck1*fod_.obj_dist; @@ -272,18 +273,18 @@ void OpticalSystem::update_paraxial_data() fod_.n_obj = n_0; fod_.n_img = n_k; - fod_.img_ht = -fod_.opt_inv/(n_k*ax_ray.at(img).slp); - fod_.obj_ang = atan(pr_ray.at(0).slp) * 180.0/M_PI; + fod_.img_ht = -fod_.opt_inv/(n_k*ax_ray->at(img)->slp); + fod_.obj_ang = atan(pr_ray->at(0)->slp) * 180.0/M_PI; - double nu_pr0 = n_0*pr_ray.at(0).slp; - fod_.enp_dist = -pr_ray.at(1).ht/nu_pr0; + double nu_pr0 = n_0*pr_ray->at(0)->slp; + fod_.enp_dist = -pr_ray->at(1)->ht/nu_pr0; fod_.enp_radius = abs(fod_.opt_inv/nu_pr0); - fod_.exp_dist = -(pr_ray.at(img).ht/pr_ray.at(img).slp - fod_.img_dist); - fod_.exp_radius = abs( fod_.opt_inv/(n_k*pr_ray.at(img).slp) ); + fod_.exp_dist = -(pr_ray->at(img)->ht/pr_ray->at(img)->slp - fod_.img_dist); + fod_.exp_radius = abs( fod_.opt_inv/(n_k*pr_ray->at(img)->slp) ); - fod_.obj_na = n_0*sin( atan(ax_ray.at(0).slp) ); - fod_.img_na = n_k*sin( atan(ax_ray.at(img).slp) ); + fod_.obj_na = n_0*sin( atan(ax_ray->at(0)->slp) ); + fod_.img_na = n_k*sin( atan(ax_ray->at(img)->slp) ); } void OpticalSystem::update_reference_rays() @@ -309,12 +310,14 @@ void OpticalSystem::update_reference_rays() tracer->set_apply_vig(true); for(int fi = 0; fi < num_fld; fi++) { + Field* fld = opt_spec_->field_of_view()->field(fi); for(int wi = 0; wi < num_wvl; wi++) { - ref_rays1_.push_back( tracer->trace_pupil_ray(pupil_chief, fi, wi) ); - ref_rays2_.push_back( tracer->trace_pupil_ray(pupil_upper_mer, fi, wi) ); - ref_rays3_.push_back( tracer->trace_pupil_ray(pupil_lower_mer, fi, wi) ); - ref_rays4_.push_back( tracer->trace_pupil_ray(pupil_upper_sag, fi, wi) ); - ref_rays5_.push_back( tracer->trace_pupil_ray(pupil_lower_sag, fi, wi) ); + double wvl = opt_spec_->spectral_region()->wvl(wi)->value(); + ref_rays1_.push_back( tracer->trace_pupil_ray(pupil_chief, fld, wvl) ); + ref_rays2_.push_back( tracer->trace_pupil_ray(pupil_upper_mer, fld, wvl) ); + ref_rays3_.push_back( tracer->trace_pupil_ray(pupil_lower_mer, fld, wvl) ); + ref_rays4_.push_back( tracer->trace_pupil_ray(pupil_upper_sag, fld, wvl) ); + ref_rays5_.push_back( tracer->trace_pupil_ray(pupil_lower_sag, fld, wvl) ); // set name @@ -350,6 +353,9 @@ void OpticalSystem::update_model() { opt_assembly_->update_model(); + num_wvl_ = opt_spec_->spectral_region()->wvl_count(); + num_fld_ = opt_spec_->field_of_view()->field_count(); + update_paraxial_data(); update_aim_pt(); update_reference_rays(); diff --git a/geopter/optical/test/triplet.cpp b/geopter/optical/test/triplet.cpp index 9cc00d7..542c00f 100644 --- a/geopter/optical/test/triplet.cpp +++ b/geopter/optical/test/triplet.cpp @@ -84,8 +84,8 @@ int main() auto pr_ray = sys->principle_ray(ref_wi); auto fod = sys->first_order_data(); - ax_ray.print(oss); - pr_ray.print(oss); + ax_ray->print(oss); + pr_ray->print(oss); fod.print(oss); @@ -94,16 +94,16 @@ int main() */ oss << "Upper marginal ray at F0" << std::endl; - Ray ray = sys->reference_ray(2, 0, ref_wi); - ray.print(oss); + auto ray = sys->reference_ray(2, 0, ref_wi); + ray->print(oss); oss << "Upper marginal ray at F2" << std::endl; ray = sys->reference_ray(2, 2, ref_wi); - ray.print(oss); + ray->print(oss); oss << "Lower marginal ray at F2" << std::endl; ray = sys->reference_ray(3, 2, ref_wi); - ray.print(oss); + ray->print(oss); // output to text std::ofstream fout;