diff --git a/geopter/gui/Analysis/renderer_qcp.cpp b/geopter/gui/Analysis/renderer_qcp.cpp index d6f1d92..3b03247 100644 --- a/geopter/gui/Analysis/renderer_qcp.cpp +++ b/geopter/gui/Analysis/renderer_qcp.cpp @@ -50,7 +50,7 @@ void RendererQCP::set_current_cell(int row, int col) } -void RendererQCP::draw_line(Eigen::Vector2d p1, Eigen::Vector2d p2, const Rgb& color, int line_style) +void RendererQCP::draw_line(const Eigen::Vector2d& p1, const Eigen::Vector2d& p2, const Rgb& color, int line_style) { QCPAxisRect *axisRect = customPlot_->axisRect(current_cell_index_); @@ -84,7 +84,7 @@ void RendererQCP::draw_line(Eigen::Vector2d p1, Eigen::Vector2d p2, const Rgb& c line->setVisible(true); } -void RendererQCP::draw_polyline(std::vector &pts, const Rgb& color, int line_style) +void RendererQCP::draw_polyline(const std::vector &pts, const Rgb& color, int line_style) { QCPAxisRect *axisRect = customPlot_->axisRect(current_cell_index_); //QCPCurve* polyline = new QCPCurve(_customPlot->xAxis, _customPlot->yAxis); @@ -116,7 +116,7 @@ void RendererQCP::draw_polyline(std::vector &pts, const Rgb& co } -void RendererQCP::draw_polyline(std::vector &x, std::vector &y, const Rgb& color, int line_style) +void RendererQCP::draw_polyline(const std::vector &x, const std::vector &y, const Rgb& color, int line_style) { int pointCount = x.size(); QVector curveData(pointCount); @@ -148,7 +148,7 @@ void RendererQCP::draw_polyline(std::vector &x, std::vector &y, polyline->setVisible(true); } -void RendererQCP::draw_dots(std::vector &x, std::vector &y, const Rgb &color, double dot_size) +void RendererQCP::draw_dots(const std::vector &x, const std::vector &y, const Rgb &color, double dot_size) { int pointCount = x.size(); QVector dotsData(pointCount); @@ -197,9 +197,11 @@ void RendererQCP::set_y_axis_label(std::string label) void RendererQCP::set_aspect_ratio(double h_per_v) { - //_customPlot->yAxis->setScaleRatio(_customPlot->xAxis, h_per_v); + customPlot_->yAxis->setScaleRatio(customPlot_->xAxis, h_per_v); QCPAxis *xAxis = customPlot_->axisRect(current_cell_index_)->axis(QCPAxis::atBottom); - customPlot_->axisRect(current_cell_index_)->axis(QCPAxis::atLeft)->setScaleRatio(xAxis, h_per_v); + QCPAxis *yAxis = customPlot_->axisRect(current_cell_index_)->axis(QCPAxis::atLeft); + + xAxis->setScaleRatio(yAxis, h_per_v); } void RendererQCP::set_mouse_interaction(bool state) diff --git a/geopter/gui/Analysis/renderer_qcp.h b/geopter/gui/Analysis/renderer_qcp.h index 5ea1fbd..a53c53a 100644 --- a/geopter/gui/Analysis/renderer_qcp.h +++ b/geopter/gui/Analysis/renderer_qcp.h @@ -19,10 +19,10 @@ class RendererQCP : public Renderer void set_grid_layout(int rows, int cols) override; void set_current_cell(int row, int col) override; - void draw_line(Eigen::Vector2d p1, Eigen::Vector2d p2, const Rgb& color, int line_style=0) override; - void draw_polyline(std::vector& pts, const Rgb& color, int line_style=0) override; - void draw_polyline(std::vector& x, std::vector& y, const Rgb& color, int line_style=0) override; - void draw_dots(std::vector& x, std::vector& y, const Rgb& color, double dot_size) override; + void draw_line(const Eigen::Vector2d& p1, const Eigen::Vector2d& p2, const Rgb& color, int line_style=0) override; + void draw_polyline(const std::vector& pts, const Rgb& color, int line_style=0) override; + void draw_polyline(const std::vector& x, const std::vector& y, const Rgb& color, int line_style=0) override; + void draw_dots(const std::vector& x, const std::vector& y, const Rgb& color, double dot_size) override; void set_x_axis_range(double xmin, double xmax) override; void set_y_axis_range(double ymin, double ymax) override; diff --git a/geopter/optical/include/Assembly/surface.h b/geopter/optical/include/Assembly/surface.h index 31e3560..d536fa6 100644 --- a/geopter/optical/include/Assembly/surface.h +++ b/geopter/optical/include/Assembly/surface.h @@ -13,7 +13,6 @@ class Surface : public Interface Surface(double r); ~Surface(); - void update() override; /** * @brief Compute intersect point and distance from the previous interface diff --git a/geopter/optical/include/Assembly/surface_profile.h b/geopter/optical/include/Assembly/surface_profile.h index 0592a77..ecd8434 100644 --- a/geopter/optical/include/Assembly/surface_profile.h +++ b/geopter/optical/include/Assembly/surface_profile.h @@ -20,8 +20,6 @@ class SurfaceProfile SurfaceProfile(); virtual ~SurfaceProfile(); - virtual void update(); - virtual std::string name() const; virtual void set_cv(double c); diff --git a/geopter/optical/include/Paraxial/paraxial_ray.h b/geopter/optical/include/Paraxial/paraxial_ray.h index 5cdb048..4d42687 100644 --- a/geopter/optical/include/Paraxial/paraxial_ray.h +++ b/geopter/optical/include/Paraxial/paraxial_ray.h @@ -8,30 +8,75 @@ namespace geopter { /** Paraxial ray parameter at each surface */ -struct ParaxialRayAtSurface +class ParaxialRayAtSurface { +public: ParaxialRayAtSurface(){ - ht = 0.0; - slp = 0.0; - aoi = 0.0; - n = 1.0; + y_ = 0.0; + u_prime_ = 0.0; + i_ = 0.0; + n_prime_ = 0.0; + before_ = nullptr; + } + + ParaxialRayAtSurface(double y, double u_prime, double i, double n_prime){ + y_ = y; + u_prime_ = u_prime; + i_ = i; + n_prime_ = n_prime; + before_ = nullptr; + } + + void set_before(ParaxialRayAtSurface* before){ + before_ = before; + } + + double y(){ + return y_; + } + + /** Returns incoming slope */ + double u(){ + if(before_){ + return before_->u_prime(); + }else{ // assume at object + return u_prime_; + } + } + + /** Returns outgoing slope */ + double u_prime(){ + return u_prime_; + }; + + double i(){ + return i_; } - ParaxialRayAtSurface(double y, double u, double i, double ind){ - ht = y; - slp = u; - aoi = i; - n = ind; + double n(){ + if(before_){ + return before_->n_prime(); + }else{ + return n_prime_; + } } + double n_prime(){ + return n_prime_; + } + + /** Returns paraxial conjugate point distance */ double l_prime() { - return (- ht/slp); + return (- y_/u_prime_); } - double ht; - double slp; // outgoing - double aoi; // incident - double n; // incident +private: + double y_; + double u_prime_; // outgoing + double i_; // incident + double n_prime_; // outgoing + + ParaxialRayAtSurface *before_; }; @@ -43,10 +88,10 @@ class ParaxialRay ~ParaxialRay(); void prepend(std::shared_ptr par_ray_at_srf); - void prepend(double ht, double slp, double aoi= 0, double n = 1.0); + void prepend(double y, double u_prime, double i= 0, double n_prime = 1.0); void append(std::shared_ptr par_ray_at_srf); - void append(double ht, double slp, double aoi= 0, double n = 1.0); + void append(double y, double u_prime, double i= 0, double n_prime = 1.0); void clear(); diff --git a/geopter/optical/include/Renderer/renderer.h b/geopter/optical/include/Renderer/renderer.h index 19cc68c..0b5d30e 100644 --- a/geopter/optical/include/Renderer/renderer.h +++ b/geopter/optical/include/Renderer/renderer.h @@ -24,11 +24,11 @@ class Renderer virtual void set_grid_layout(int rows, int cols); virtual void set_current_cell(int row, int col); - virtual void draw_line(Eigen::Vector2d p1, Eigen::Vector2d p2, const Rgb& color, int line_style=0); - virtual void draw_polyline(std::vector& pts, const Rgb& color, int line_style=0); - virtual void draw_polyline(std::vector& x, std::vector& y, const Rgb& color, int line_style=0); + virtual void draw_line(const Eigen::Vector2d& p1, const Eigen::Vector2d& p2, const Rgb& color, int line_style=0); + virtual void draw_polyline(const std::vector& pts, const Rgb& color, int line_style=0); + virtual void draw_polyline(const std::vector& x, const std::vector& y, const Rgb& color, int line_style=0); virtual void draw_text(std::string str, const Rgb& color); - virtual void draw_dots(std::vector& x, std::vector& y, const Rgb& color, double dot_size); + virtual void draw_dots(const std::vector& x, const std::vector& y, const Rgb& color, double dot_size); virtual void draw_x_axis(bool state= true); virtual void draw_y_axis(bool state= true); diff --git a/geopter/optical/include/Sequential/ray.h b/geopter/optical/include/Sequential/ray.h index e19e576..65d55ff 100644 --- a/geopter/optical/include/Sequential/ray.h +++ b/geopter/optical/include/Sequential/ray.h @@ -47,6 +47,9 @@ struct RayAtSurface /** Ray height from the axis */ double height() const { return sqrt( pow(intersect_pt(0),2) + pow(intersect_pt(1),2) ); } + +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW }; diff --git a/geopter/optical/include/Sequential/sequential_trace.h b/geopter/optical/include/Sequential/sequential_trace.h index d9d9e22..c473eb4 100644 --- a/geopter/optical/include/Sequential/sequential_trace.h +++ b/geopter/optical/include/Sequential/sequential_trace.h @@ -18,17 +18,23 @@ class SequentialTrace SequentialTrace(OpticalSystem* sys); ~SequentialTrace(); - //std::shared_ptr trace_pupil_ray(Eigen::Vector2d pupil_crd, int fi, int wi); - std::shared_ptr trace_pupil_ray(const Eigen::Vector2d& pupil_crd, const Field* fld, double wvl); + /** + * @brief Trace patterned rays (grid, hexapolar, etc) + * @param rays Traced Ray vector + * @param pupils pupil coordinated list + * @param fld Field pointer + * @param wvl wavelength value + * @return number of valid rays + */ + int trace_pupil_pattern_rays(std::vector< std::shared_ptr >& rays, const std::vector< Eigen::Vector2d >& pupils, const Field* fld, double wvl); + /** Trace a ray throughout the given sequantial path */ 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(const Field* fld, double wvl); /** @@ -41,18 +47,15 @@ class SequentialTrace * @return Eigen::Vector2d aim point on paraxial entrance pupil plane */ 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 * @param normal normal of the surface * @return direction of refracted ray */ - Eigen::Vector3d bend(Eigen::Vector3d d_in, Eigen::Vector3d normal, double n_in, double n_out); + Eigen::Vector3d bend(const Eigen::Vector3d& d_in, const 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(const Field* fld); /** Get sequential path between start and end */ @@ -71,9 +74,9 @@ class SequentialTrace bool apply_vig_status() const; private: - double y_stop_coordinate(double y1, int ifcx, Eigen::Vector3d pt0, double dist, double wvl, double y_target); + double y_stop_coordinate(double y1, int ifcx, const Eigen::Vector3d& pt0, double dist, double wvl, double y_target); - double compute_vignetting_factor_for_pupil(Eigen::Vector2d full_pupil, const Field& fld); + double compute_vignetting_factor_for_pupil(const Eigen::Vector2d& full_pupil, const Field& fld); OpticalSystem *opt_sys_; diff --git a/geopter/optical/include/Spec/field.h b/geopter/optical/include/Spec/field.h index a388489..b2eeaaf 100644 --- a/geopter/optical/include/Spec/field.h +++ b/geopter/optical/include/Spec/field.h @@ -25,6 +25,8 @@ class Field double vuy() const; double weight() const; Rgb render_color() const; + Eigen::Vector2d aim_pt() const; + Eigen::Vector3d object_coord() const; void set_x(double x); void set_y(double y); @@ -34,26 +36,26 @@ class Field void set_vuy(double vuy); void set_weight(double wt); void set_render_color(const Rgb& color); - - Eigen::Vector2d aim_pt() const; void set_aim_pt(PupilCrd aim_pt); + void set_object_coord(Eigen::Vector3d obj_crd); PupilCrd apply_vignetting(PupilCrd pupil) const; - void update(); + EIGEN_MAKE_ALIGNED_OPERATOR_NEW protected: double x_; double y_; + double wt_; double vux_; double vuy_; double vlx_; double vly_; - double wt_; - Eigen::Vector2d aim_pt_; - + Eigen::Vector3d object_coord_; Rgb render_color_; + + }; } //namespace geopter diff --git a/geopter/optical/include/Spec/pupil_spec.h b/geopter/optical/include/Spec/pupil_spec.h index 6f44b6a..e4ff529 100644 --- a/geopter/optical/include/Spec/pupil_spec.h +++ b/geopter/optical/include/Spec/pupil_spec.h @@ -40,6 +40,8 @@ class PupilSpec void print(std::ostringstream& oss); + EIGEN_MAKE_ALIGNED_OPERATOR_NEW + private: int pupil_type_; double value_; diff --git a/geopter/optical/include/Spec/wvl.h b/geopter/optical/include/Spec/wvl.h index 313e312..c4b2d18 100644 --- a/geopter/optical/include/Spec/wvl.h +++ b/geopter/optical/include/Spec/wvl.h @@ -2,7 +2,6 @@ #define WVL_H #include "spectral_line.h" - #include "Renderer/rgb.h" namespace geopter{ @@ -30,6 +29,7 @@ class Wvl void set_render_color(Rgb color); + /** Returns render color */ Rgb render_color() const; private: diff --git a/geopter/optical/include/System/optical_system.h b/geopter/optical/include/System/optical_system.h index 3b91061..5d723b5 100644 --- a/geopter/optical/include/System/optical_system.h +++ b/geopter/optical/include/System/optical_system.h @@ -66,6 +66,7 @@ class OpticalSystem void update_reference_rays(); void update_semi_diameters(); void update_paraxial_data(); + void update_object_coords(); std::unique_ptr opt_assembly_; std::unique_ptr opt_spec_; diff --git a/geopter/optical/src/Analysis/aberration.cpp b/geopter/optical/src/Analysis/aberration.cpp index d60a8d2..974feb8 100644 --- a/geopter/optical/src/Analysis/aberration.cpp +++ b/geopter/optical/src/Analysis/aberration.cpp @@ -1,3 +1,6 @@ +#define _USE_MATH_DEFINES +#include + #include #include #include @@ -144,20 +147,16 @@ void Aberration::plot_longitudinal_spherical_aberration(double scale) int ref_wvl_idx = opt_sys_->optical_spec()->spectral_region()->reference_index(); int num_wvls = opt_sys_->optical_spec()->spectral_region()->wvl_count(); - double img_dst = opt_sys_->first_order_data().img_dist; // collect l_prime std::vector l_primes; { int num_srf = opt_sys_->optical_assembly()->surface_count(); - int img_srf = num_srf-1; 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++){ 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; + double l_prime = ax_ray->at(last_surf)->l_prime(); l_primes.push_back(l_prime); } } @@ -240,14 +239,12 @@ void Aberration::plot_astigmatism(double scale) std::vector l_primes; { int num_srf = opt_sys_->optical_assembly()->surface_count(); - int img_srf = num_srf-1; + //int img_srf = num_srf-1; 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++){ 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; + double l_prime = ax_ray->at(last_surf)->l_prime(); l_primes.push_back(l_prime); } } @@ -324,8 +321,8 @@ void Aberration::plot_chromatic_focus_shift(double lower_wvl, double higher_wvl) tracer->get_starting_coords(&y0, &u0, &ybar0, &ubar0); 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 y = ref_prx_ray->at(last_surf)->y(); + double u_prime = ref_prx_ray->at(img_srf)->u_prime(); double ref_l_prime = -y/u_prime; @@ -335,8 +332,8 @@ void Aberration::plot_chromatic_focus_shift(double lower_wvl, double higher_wvl) double wvl = lower_wvl + i*wvl_step; 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)->y(); + u_prime = prx_ray->at(img_srf)->u_prime(); double l_prime = -y/u_prime; ydata.push_back(l_prime - ref_l_prime); xdata.push_back(wvl); @@ -381,7 +378,7 @@ void Aberration::plot_spot_diagram(int pattern, int nrd, double scale, double do } std::vector ray_dx, ray_dy; - std::shared_ptr ray; + std::vector< std::shared_ptr > rays; const WvlSpec* wvl_spec = opt_sys_->optical_spec()->spectral_region(); @@ -399,12 +396,16 @@ void Aberration::plot_spot_diagram(int pattern, int nrd, double scale, double do for(int wi = 0; wi < num_wvl_; wi++){ double wvl = wvl_spec->wvl(wi)->value(); + + int valid_ray_count = tracer->trace_pupil_pattern_rays(rays, pupils, fld, wvl); + ray_dx.clear(); ray_dy.clear(); - for(auto& pupil : pupils){ - ray = tracer->trace_pupil_ray(pupil, fld, wvl); + ray_dx.reserve(valid_ray_count); + ray_dy.reserve(valid_ray_count); + for(auto& ray : rays){ 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); diff --git a/geopter/optical/src/Assembly/surface.cpp b/geopter/optical/src/Assembly/surface.cpp index 571887e..0ea39dc 100644 --- a/geopter/optical/src/Assembly/surface.cpp +++ b/geopter/optical/src/Assembly/surface.cpp @@ -23,11 +23,6 @@ Surface::~Surface() profile_.reset(); } -void Surface::update() -{ - profile_->update(); -} - void Surface::intersect(Eigen::Vector3d &pt, double &s, const Eigen::Vector3d& p0, const Eigen::Vector3d& d, double eps, double z_dir) { diff --git a/geopter/optical/src/Assembly/surface_profile.cpp b/geopter/optical/src/Assembly/surface_profile.cpp index 6ffc547..f79923d 100644 --- a/geopter/optical/src/Assembly/surface_profile.cpp +++ b/geopter/optical/src/Assembly/surface_profile.cpp @@ -13,11 +13,6 @@ SurfaceProfile::~SurfaceProfile() } -void SurfaceProfile::update() -{ - -} - std::string SurfaceProfile::name() const { return name_; @@ -100,7 +95,7 @@ void SurfaceProfile::intersect_spencer(Eigen::Vector3d& pt, double& s, const Eig int iter = 0; - while(delta > eps && iter < 50) + while(delta > eps && iter < 30) { p = p0 + s1*d; s2 = s1 - f(p)/d.dot(df(p)); diff --git a/geopter/optical/src/CMakeLists.txt b/geopter/optical/src/CMakeLists.txt index 8d21e3f..1dd192e 100644 --- a/geopter/optical/src/CMakeLists.txt +++ b/geopter/optical/src/CMakeLists.txt @@ -71,7 +71,6 @@ set(OPTICAL_SRCS Utility/utility.cpp ) - add_library(geopteroptical STATIC ${OPTICAL_SRCS}) target_include_directories(geopteroptical PUBLIC ${CMAKE_SOURCE_DIR}../include) diff --git a/geopter/optical/src/Paraxial/paraxial_ray.cpp b/geopter/optical/src/Paraxial/paraxial_ray.cpp index b9b86b7..dc7c8a9 100644 --- a/geopter/optical/src/Paraxial/paraxial_ray.cpp +++ b/geopter/optical/src/Paraxial/paraxial_ray.cpp @@ -33,21 +33,29 @@ void ParaxialRay::prepend(std::shared_ptr par_ray_at_srf) par_ray_at_srfs_.insert(itr, par_ray_at_srf); } -void ParaxialRay::prepend(double ht, double slp, double aoi, double n) +void ParaxialRay::prepend(double y, double u_prime, double i, double n_prime) { - auto prx_ray_at_srf = std::make_shared(ht,slp,aoi,n); + auto prx_ray_at_srf = std::make_shared(y, u_prime, i, n_prime); auto itr = par_ray_at_srfs_.begin(); par_ray_at_srfs_.insert( itr, prx_ray_at_srf); } void ParaxialRay::append(std::shared_ptr par_ray_at_srf) { + ParaxialRayAtSurface *before = par_ray_at_srfs_.back().get(); + par_ray_at_srf->set_before(before); par_ray_at_srfs_.push_back(par_ray_at_srf); } -void ParaxialRay::append(double ht, double slp, double aoi, double n) +void ParaxialRay::append(double y, double u_prime, double i, double n_prime) { - auto prx_ray_at_srf = std::make_shared(ht, slp, aoi, n); + auto prx_ray_at_srf = std::make_shared(y, u_prime, i, n_prime); + + if(par_ray_at_srfs_.size() > 0){ + ParaxialRayAtSurface *before = par_ray_at_srfs_.back().get(); + prx_ray_at_srf->set_before(before); + } + par_ray_at_srfs_.push_back(prx_ray_at_srf); } @@ -110,15 +118,11 @@ 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; - if( i == num_srf-1){ - 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]->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]->y(); + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->u_prime(); + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->n_prime()*par_ray_at_srfs_[i]->u_prime(); + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->i(); + oss << std::setw(val_w) << std::right << std::fixed << std::setprecision(prec) << par_ray_at_srfs_[i]->n() * par_ray_at_srfs_[i]->i(); 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 1669fa5..0e3bb68 100644 --- a/geopter/optical/src/Paraxial/paraxial_trace.cpp +++ b/geopter/optical/src/Paraxial/paraxial_trace.cpp @@ -41,7 +41,8 @@ std::shared_ptr ParaxialTrace::trace_paraxial_ray_from_object(doubl double u0_prime = u0; double i0 = y0*c0 + u0; double n0 = par_path.at(0).n; - par_ray->append(y0, u0_prime, i0, n0); + double n0_prime = n0; + par_ray->append(y0, u0_prime, i0, n0_prime); if(path_size == 1) { return par_ray; @@ -65,7 +66,7 @@ std::shared_ptr ParaxialTrace::trace_paraxial_ray_from_object(doubl 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_prime); // transfer to next d = par_path.at(k).t; @@ -204,18 +205,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()->y(); + double bk1 = q_ray->back()->y(); + double ck1 = n_k*p_ray->back()->u_prime(); + double dk1 = n_k*q_ray->back()->u_prime(); 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)->y(); + double bs1 = q_ray->at(stop)->y(); + double cs1 = n_s*p_ray->at(stop)->u_prime(); + double ds1 = n_s*q_ray->at(stop)->u_prime(); double ybar1 = -bs1; double ubar1 = as1; diff --git a/geopter/optical/src/Renderer/renderer.cpp b/geopter/optical/src/Renderer/renderer.cpp index 4fc2248..b4c2278 100644 --- a/geopter/optical/src/Renderer/renderer.cpp +++ b/geopter/optical/src/Renderer/renderer.cpp @@ -22,22 +22,22 @@ void Renderer::draw_text(std::string str, const Rgb& color) } -void Renderer::draw_line(Eigen::Vector2d p1, Eigen::Vector2d p2, const Rgb& color, int line_style) +void Renderer::draw_line(const Eigen::Vector2d& p1, const Eigen::Vector2d& p2, const Rgb& color, int line_style) { } -void Renderer::draw_polyline(std::vector& pts, const Rgb& color, int line_style) +void Renderer::draw_polyline(const std::vector& pts, const Rgb& color, int line_style) { } -void Renderer::draw_polyline(std::vector& x, std::vector& y, const Rgb& color, int line_style) +void Renderer::draw_polyline(const std::vector& x, const std::vector& y, const Rgb& color, int line_style) { } -void Renderer::draw_dots(std::vector &x, std::vector &y, const Rgb &color, double dot_size) +void Renderer::draw_dots(const std::vector &x, const std::vector &y, const Rgb &color, double dot_size) { } diff --git a/geopter/optical/src/Sequential/sequential_trace.cpp b/geopter/optical/src/Sequential/sequential_trace.cpp index 32e0109..dcaf251 100644 --- a/geopter/optical/src/Sequential/sequential_trace.cpp +++ b/geopter/optical/src/Sequential/sequential_trace.cpp @@ -36,15 +36,6 @@ SequentialTrace::~SequentialTrace() opt_sys_ = nullptr; } -/* -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); -} -*/ std::shared_ptr SequentialTrace::trace_pupil_ray(const Eigen::Vector2d& pupil_crd, const Field *fld, double wvl) { @@ -65,13 +56,9 @@ std::shared_ptr SequentialTrace::trace_pupil_ray(const Eigen::Vector2d& pup pt1(1) = eprad*vig_pupil(1) + aim_pt(1); pt1(2) = fod.obj_dist + fod.enp_dist; - Eigen::Vector3d dir0; - dir0 = pt1 - pt0; - //double length = dir0.norm(); - //dir0 = dir0/length; + Eigen::Vector3d dir0 = pt1 - pt0; dir0.normalize(); - //int img = opt_sys_->optical_assembly()->image_index(); SequentialPath path = sequential_path(0, image_index_, wvl); std::shared_ptr ray = trace_ray_throughout_path(path, pt0, dir0); @@ -80,11 +67,55 @@ std::shared_ptr SequentialTrace::trace_pupil_ray(const Eigen::Vector2d& pup return ray; } +int SequentialTrace::trace_pupil_pattern_rays(std::vector< std::shared_ptr >& rays, const std::vector& pupils, const Field* fld, double wvl) +{ + + + auto fod = opt_sys_->first_order_data(); + double eprad = fod.enp_radius; + PupilCrd aim_pt = fld->aim_pt(); + + //Eigen::Vector3d pt0 = this->object_coord(fld); + Eigen::Vector3d pt0 = fld->object_coord(); + + SequentialPath path = sequential_path(0, image_index_, wvl); + + rays.clear(); + rays.reserve(pupils.size()); + + int valid_count = 0; + + for(auto& pupil : pupils) + { + Eigen::Vector2d vig_pupil = pupil; + + if(do_apply_vig_){ + vig_pupil = fld->apply_vignetting(pupil); + } + + Eigen::Vector3d pt1; + pt1(0) = eprad*vig_pupil(0) + aim_pt(0); + pt1(1) = eprad*vig_pupil(1) + aim_pt(1); + pt1(2) = fod.obj_dist + fod.enp_dist; + + Eigen::Vector3d dir0 = pt1 - pt0; + dir0.normalize(); + + auto ray = trace_ray_throughout_path(path, pt0, dir0); + ray->set_wvl(wvl); + + if(ray->status() == PassThrough){ + valid_count++; + } + + rays.push_back(ray); + } + return valid_count; +} std::shared_ptr SequentialTrace::trace_ray_throughout_path(const SequentialPath& seq_path, const Eigen::Vector3d& pt0, const Eigen::Vector3d& dir0) { std::shared_ptr ray = std::make_shared(); - ray->set_status(RayStatus::PassThrough); if(seq_path.size() == 0) { return ray; // empty ray @@ -194,12 +225,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; + double y0 = opt_sys_->axial_ray(ref_wvl_idx_)->at(0)->y(); + double u0 = opt_sys_->axial_ray(ref_wvl_idx_)->at(0)->u(); 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)->y(); + double u_prime = ax_ray->at(img_srf)->u_prime(); double l_prime = -y/u_prime; delete parax_tracer; @@ -302,9 +333,6 @@ 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(); - assert(end <= num_srfs_); SequentialPath path; @@ -332,19 +360,9 @@ 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); - double wvl = opt_sys_->optical_spec()->spectral_region()->wavelength(wi); - - return aim_chief_ray(fld, wvl); -} -*/ Eigen::Vector2d SequentialTrace::aim_chief_ray(const Field *fld, double wvl) { - //double wvl = opt_sys_->optical_spec()->spectral_region()->wavelength(wi); int stop = opt_sys_->optical_assembly()->stop_index(); Eigen::Vector2d xy_target({0.0, 0.0}); Eigen::Vector2d aim_pt = search_aim_point(stop, xy_target, fld, wvl); @@ -352,15 +370,6 @@ 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); - double wvl = opt_sys_->optical_spec()->spectral_region()->wavelength(wi); - - return search_aim_point(srf_idx, xy_target, fld, wvl); -} -*/ Eigen::Vector2d SequentialTrace::search_aim_point(int srf_idx, const Eigen::Vector2d& xy_target, const Field *fld, double wvl) { @@ -416,7 +425,7 @@ Eigen::Vector2d SequentialTrace::search_aim_point(int srf_idx, const Eigen::Vect } -double SequentialTrace::y_stop_coordinate(double y1, int ifcx, Eigen::Vector3d pt0, double dist, double wvl, double y_target) +double SequentialTrace::y_stop_coordinate(double y1, int ifcx, const Eigen::Vector3d& pt0, double dist, double wvl, double y_target) { Eigen::Vector3d pt1({0.0, y1, dist}); Eigen::Vector3d dir0 = pt1 - pt0; @@ -453,7 +462,7 @@ double SequentialTrace::y_stop_coordinate(double y1, int ifcx, Eigen::Vector3d p } -Eigen::Vector3d SequentialTrace::bend(Eigen::Vector3d d_in, Eigen::Vector3d normal, double n_in, double n_out) +Eigen::Vector3d SequentialTrace::bend(const Eigen::Vector3d& d_in, const Eigen::Vector3d& normal, double n_in, double n_out) { try { double normal_len = normal.norm(); @@ -468,13 +477,6 @@ 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) { @@ -532,7 +534,7 @@ std::vector SequentialTrace::compute_vignetting_factors(const Field& fld } -double SequentialTrace::compute_vignetting_factor_for_pupil(Eigen::Vector2d full_pupil, const Field &fld) +double SequentialTrace::compute_vignetting_factor_for_pupil(const Eigen::Vector2d& full_pupil, const Field &fld) { bool orig_vig_state = do_apply_vig_; do_apply_vig_ = false; diff --git a/geopter/optical/src/Spec/field.cpp b/geopter/optical/src/Spec/field.cpp index ab9ff6c..7851255 100644 --- a/geopter/optical/src/Spec/field.cpp +++ b/geopter/optical/src/Spec/field.cpp @@ -13,6 +13,7 @@ Field::Field(double x, double y, double wt, Rgb color, double vuy, double vly, d render_color_(color) { aim_pt_ = Eigen::Vector2d::Zero(2); + object_coord_ = Eigen::Vector3d::Zero(3); } Field::~Field() @@ -110,6 +111,16 @@ void Field::set_aim_pt(PupilCrd aim_pt) aim_pt_ = aim_pt; } +Eigen::Vector3d Field::object_coord() const +{ + return object_coord_; +} + +void Field::set_object_coord(Eigen::Vector3d obj_crd) +{ + object_coord_ = obj_crd; +} + PupilCrd Field::apply_vignetting(PupilCrd pupil) const { PupilCrd vig_pupil = pupil; @@ -129,7 +140,4 @@ PupilCrd Field::apply_vignetting(PupilCrd pupil) const return vig_pupil; } -void Field::update() -{ - //_chief_ray = nullptr; -} + diff --git a/geopter/optical/src/System/optical_system.cpp b/geopter/optical/src/System/optical_system.cpp index 3fff5c3..12b676f 100644 --- a/geopter/optical/src/System/optical_system.cpp +++ b/geopter/optical/src/System/optical_system.cpp @@ -152,6 +152,50 @@ void OpticalSystem::update_aim_pt() } } +void OpticalSystem::update_object_coords() +{ + Eigen::Vector3d obj_pt; + + Eigen::Vector3d ang_dg; + Eigen::Vector3d img_pt; + Eigen::Vector3d dir_tan; + + const int field_type = opt_spec_->field_of_view()->field_type(); + + for(int fi = 0; fi < num_fld_; fi++){ + Field* fld = opt_spec_->field_of_view()->field(fi); + double fld_x = fld->x(); + double fld_y = fld->y(); + + switch (field_type) + { + case FieldType::OBJ_ANG: + ang_dg = Eigen::Vector3d({fld_x, fld_y, 0.0}); + dir_tan(0) = tan(ang_dg(0) * M_PI/180.0); + dir_tan(1) = tan(ang_dg(1) * M_PI/180.0); + dir_tan(2) = tan(ang_dg(2) * M_PI/180.0); + obj_pt = -dir_tan*(fod_.obj_dist + fod_.enp_dist); + break; + + case FieldType::OBJ_HT: + obj_pt(0) = fld_x; + obj_pt(1) = fld_y; + obj_pt(2) = 0.0; + break; + + case FieldType::IMG_HT: + img_pt = Eigen::Vector3d({fld_x, fld_y, 0.0}); + obj_pt = fod_.red*img_pt; + break; + + default: + obj_pt = Eigen::Vector3d::Zero(3); + } + + fld->set_object_coord(obj_pt); + } +} + void OpticalSystem::update_semi_diameters() { @@ -251,21 +295,21 @@ void OpticalSystem::update_paraxial_data() 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)->y(); + double bk1 = q_ray->at(img)->y(); + double ck1 = n_k*p_ray->at(img)->u_prime(); + double dk1 = n_k*q_ray->at(img)->u_prime(); // 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)->y()*pr_ray->at(0)->u_prime() - pr_ray->at(1)->y()*ax_ray->at(0)->u_prime() ); 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)->y() - 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)->u_prime()); fod_.m = ak1 + ck1*fod_.img_dist/n_k; fod_.red = dk1 + ck1*fod_.obj_dist; @@ -273,18 +317,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)->u_prime()); + fod_.obj_ang = atan(pr_ray->at(0)->u_prime()) * 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)->u_prime(); + fod_.enp_dist = -pr_ray->at(1)->y()/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)->y()/pr_ray->at(img)->u_prime() - fod_.img_dist); + fod_.exp_radius = abs( fod_.opt_inv/(n_k*pr_ray->at(img)->u_prime()) ); - 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)->u_prime()) ); + fod_.img_na = n_k*sin( atan(ax_ray->at(img)->u_prime()) ); } void OpticalSystem::update_reference_rays() @@ -357,6 +401,7 @@ void OpticalSystem::update_model() num_fld_ = opt_spec_->field_of_view()->field_count(); update_paraxial_data(); + update_object_coords(); update_aim_pt(); update_reference_rays(); update_semi_diameters();