18544749821

Committed 15 Oct 2025 10:47PM UTC coverage: 89.946% (+45.8%) from 44.166%

Build # 18544749821

Build Type

push

github

Committed by

web-flow

Commit Message

Restructured backend for SolTrace (#63)

* Initial API for refactored SolTrace

* Updates to various APIs; first pass at element/ray source container implementation

* Use existing matrix-vector code as backend for matrix and vector classes; start implementing APIs

* Created new test on power tower example file with 50000 rays. Changed CMake to use cpp 17.

* Added ground results csv to new folder in google-tests directory

* Removed register prefix from variables in mtrand.h. More errors to come.

* Removed lines in CMakeLists causing linux failure, fixed bug in power tower test

* Created parabola.stinput test and a test using powertower optimization on the powertower sample

* Commented out parabola test for now, added nonexistent branch to CI to see how it runs

* Fixed syntax error

* Changed tests based on PR feedback

* Fixed error in CMakeLists that caused failure on windows

* Removed output messages from st_sim_run_test

* Fixed silly mistake in previous commit

* Added comments and removed unused libraries

* Changed naming conventions of tests

* Added high flux solar furnace test changes, changed parabola test file name

* Fixed typo

* Move refactored data to its own directory; add to refactored classes to build

* Add missed files

* Add unit tests for basic linear algebra and container template

* More tests; Start implementation of composite element

* Add more unit tests on element

* Add missed headers

* Remove target from cmake build command in CI

* Initial virtual element implementation; smoke test for virtual element

* Correct aperture spelling; add tests for virtual elements

* Add some unit tests for simulation data

* Move to static library for simulation data (temporary); update cmake files to get correct mac architecture

* Attempt to fix coverage failure

* Attempt to capture inline functions in code coverage; add a few more explicit tests

* Fix build

* Separate storage of CompositeElement and SingleElement; update te... (continued)

Coverage Stats

4357 of 4844 new or added lines in 62 files covered. (89.95%)

4357 of 4844 relevant lines covered (89.95%)

8565975.82 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.97

/coretrace/simulation_data/cst_templates/linear_fresnel.cpp

#include "linear_fresnel.hpp"

#include <stdexcept>
#include <sstream>

#include "aperture.hpp"
#include "constants.hpp"
#include "element.hpp"
#include "surface.hpp"

#include "cst_templates/utilities.hpp"

namespace SolTrace::Data {

LinearFresnel::LinearFresnel()
    : CompositeElement(),
      initialized(false),
      receiver_height(-1.0),
      azimuth(-1.0),
      tilt(-1.0),
      focused_panels(false),
      aperture_size_x(-1.0),
      aperture_size_y(-1.0),
      num_panels_x(-1),
      num_panels_y(-1),
      gap_x(-1.0),
      gap_y(-1.0),
      gap_center(-1.0),
      abs_diameter(-1.0),
      env_diameter(-1.0),
      env_thickness(-1.0),
      // tracking_angle(0.0),
      tracking_limit_lower(-180.0),
      tracking_limit_upper(180.0)
{
    this->optics_absorber.set_ideal_absorption();
    this->optics_mirror.set_ideal_reflection();
    this->optics_env_out.set_ideal_transmission();
    this->optics_env_in.set_ideal_transmission();

    this->tracking_origin.set_values(1.0, 0.0, 0.0);
    this->rotation_axis.set_values(0.0, 1.0, 0.0);
    this->neutral_normal.set_values(0.0, 0.0, 1.0);
}

LinearFresnel::~LinearFresnel()
{
    // Clear vectors
    this->mirrors.clear();
    this->absorbers.clear();
    this->envelope.clear();

    return;
}

// double LinearFresnel::get_tracking_angle_degrees() const
// {
//     return this->tracking_angle;
// }

// double LinearFresnel::get_tracking_angle_radians() const
// {
//     return D2R * this->get_tracking_angle_degrees();
// }

void LinearFresnel::set_angles(double azimuth, double tilt)
{
    if (azimuth < -180.0 || azimuth > 180.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_angles: Invalid azimuth angle ("
           << azimuth << "). Must be between -180 and 180 degrees.";
        throw std::invalid_argument(ss.str());
    }

    if (tilt < 0.0 || tilt > 90.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_angles: Invalid tilt angle ("
           << tilt << "). Must be between 0 and 90 degrees.";
        throw std::invalid_argument(ss.str());
    }

    this->coordinates_initialized = false;
    this->azimuth = azimuth;
    this->tilt = tilt;

    // Convert input angles to spherical coordinate angles
    double az = azimuth * D2R;
    double el = tilt * D2R;
    double pol = 0.5 * PI - az;
    double inc = 0.5 * PI - el;

    // Convert spherical coordinates to cartesian coordinates
    // y-axis
    this->rotation_axis.set_values(sin(inc) * cos(pol),
                                   sin(inc) * sin(pol),
                                   cos(inc));

    // z-axis
    this->neutral_normal.set_values(sin(-el) * cos(pol),
                                    sin(-el) * sin(pol),
                                    cos(-el));

    cross_product(this->rotation_axis,
                  this->neutral_normal,
                  this->tracking_origin);

    return;
}

void LinearFresnel::set_aperture_size(double len_x, double len_y)
{
    if (len_x <= 0.0 || len_y <= 0.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_aperture_size: Invalid aperture dimensions. "
           << "len_x (" << len_x << ") and len_y (" << len_y
           << ") must be positive.";
        throw std::invalid_argument(ss.str());
    }

    this->initialized = false;
    this->aperture_size_x = len_x;
    this->aperture_size_y = len_y;

    return;
}

void LinearFresnel::set_focused_panels(bool focused)
{
    this->initialized = false;
    this->focused_panels = focused;
    return;
}

void LinearFresnel::set_gaps(double gap_x,
                             double gap_y,
                             double gap_center)
{
    if (gap_x < 0.0 || gap_y < 0.0 || gap_center < 0.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_gaps: Invalid gap dimensions. "
           << "gap_x (" << gap_x << "), gap_y (" << gap_y
           << "), and gap_center (" << gap_center
           << ") must be non-negative.";
        throw std::invalid_argument(ss.str());
    }

    this->initialized = false;
    this->gap_x = gap_x;
    this->gap_y = gap_y;
    this->gap_center = gap_center;

    return;
}

void LinearFresnel::set_number_panels(int_fast64_t num_x,
                                      int_fast64_t num_y)
{
    if (num_x <= 0 || num_y <= 0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_number_panels: Invalid panel count. "
           << "num_x (" << num_x << ") and num_y (" << num_y
           << ") must be positive.";
        throw std::invalid_argument(ss.str());
    }

    this->initialized = false;
    this->num_panels_x = num_x;
    this->num_panels_y = num_y;

    return;
}

void LinearFresnel::set_optics(const OpticalProperties &mirror,
                               const OpticalProperties &absorber,
                               const OpticalProperties &envelop_outer,
                               const OpticalProperties &envelop_inner)
{
    this->optics_mirror = mirror;
    this->optics_absorber = absorber;
    this->optics_env_out = envelop_outer;
    this->optics_env_in = envelop_inner;
    return;
}

void LinearFresnel::set_receiver_height(double height)
{
    if (height <= 0.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_receiver_height: Invalid receiver height ("
           << height << "). Height must be positive.";
        throw std::invalid_argument(ss.str());
    }

    this->initialized = false;
    this->receiver_height = height;

    return;
}

void LinearFresnel::set_receiver_dimensions(double absorber_diameter,
                                            double envelop_diameter,
                                            double envelop_thickness)
{
    if (absorber_diameter <= 0.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_receiver_dimensions: "
           << "Invalid absorber diameter (" << absorber_diameter
           << "). Must be positive.";
        throw std::invalid_argument(ss.str());
    }

    if (envelop_diameter <= 0.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_receiver_dimensions: "
           << "Invalid envelope diameter (" << envelop_diameter
           << "). Must be positive.";
        throw std::invalid_argument(ss.str());
    }

    if (envelop_thickness < 0.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_receiver_dimensions: "
           << "Invalid envelope thickness (" << envelop_thickness
           << "). Must be non-negative.";
        throw std::invalid_argument(ss.str());
    }

    if (absorber_diameter >= envelop_diameter)
    {
        std::stringstream ss;
        ss << "LinearFresnel::set_receiver_dimensions: "
           << "Absorber diameter (" << absorber_diameter
           << ") must be smaller than envelope diameter ("
           << envelop_diameter << ").";
        throw std::invalid_argument(ss.str());
    }

    this->initialized = false;
    this->abs_diameter = absorber_diameter;
    this->env_diameter = envelop_diameter;
    this->env_thickness = envelop_thickness;

    return;
}

void LinearFresnel::set_tracking_limits(double lower, double upper)
{
    if (lower > upper)
    {
        std::stringstream ss;
        ss << "ParabolicTrough::set_tracking_limits: Invalid tracking "
           << "limits. Lower limit (" << lower
           << ") must be less than or equal to upper limit ("
           << upper << ").";
        throw std::invalid_argument(ss.str());
    }

    this->tracking_limit_lower = lower;
    this->tracking_limit_upper = upper;

    return;
}

void LinearFresnel::create_geometry()
{
    if (this->initialized)
    {
        return;
    }

    this->enforce_user_fields_set();

    this->mirrors.clear();
    this->absorbers.clear();
    this->envelope.clear();
    this->clear();

    // Create mirrors
    if (this->num_panels_x == 1 && this->gap_center != 0.0)
    {
        // TODO: This should generate a warning.
        this->gap_center = 0.0;
        this->gap_x = 0.0;
    }

    double panel_len_x = this->aperture_size_x -
                         this->gap_x * (this->num_panels_x - 1);
    panel_len_x -= this->gap_center - this->gap_x;
    panel_len_x /= this->num_panels_x;
    double panel_x = -0.5 * this->aperture_size_x + 0.5 * panel_len_x;

    double panel_len_y = this->aperture_size_y -
                         this->gap_y * (this->num_panels_y - 1);
    panel_len_y /= this->num_panels_y;

    element_id sts;
    Vector3d origin;
    Vector3d aim;
    // TODO: Should mirrors aim at the receiver center or origin?
    Vector3d receiver_pos(0.0,
                          0.0,
                          this->receiver_height - 0.5 * this->abs_diameter);
    double panel_y, flen;
    single_element_ptr mirror;
    aperture_ptr ap;
    surface_ptr surf;
    Vector3d khat(0.0, 0.0, 1.0);

    for (int_fast64_t i = 0; i < this->num_panels_x; ++i)
    {
        panel_y = -0.5 * this->aperture_size_y + 0.5 * panel_len_y;
        for (int_fast64_t j = 0; j < this->num_panels_y; ++j)
        {
            mirror = make_element<SingleElement>();

            origin.set_values(panel_x, panel_y, 0.0);
            vector_add(1.0, receiver_pos, -1.0, origin, aim);
            make_unit_vector(aim);
            vector_add(0.5, khat, 0.5, aim);
            vector_add(1.0, origin, 1.0, aim);
            make_unit_vector(aim);
            aim.scalar_mult(1000.0);
            mirror->set_reference_frame_geometry(origin, aim, 0.0);

            ap = make_aperture<Rectangle>(panel_len_x, panel_len_y);
            mirror->set_aperture(ap);

            if (this->focused_panels)
            {
                flen = sqrt(panel_x * panel_x +
                            this->receiver_height * this->receiver_height);
                surf = make_surface<Parabola>(flen, 0.0);
            }
            else
            {
                surf = make_surface<Flat>();
            }
            mirror->set_surface(surf);

            mirror->set_front_optical_properties(this->optics_mirror);
            mirror->set_back_optical_properties(this->optics_mirror);
            mirror->enable();

            std::stringstream name;
            name << "Mirror_" << this->num_panels_y * i + j;
            mirror->set_name(name.str());

            this->mirrors.push_back(mirror);
            sts = this->add_element(mirror);
            if (!Element::is_success(sts))
            {
                // TODO: Make this more helpful
                throw std::runtime_error("Failed to add mirror element");
            }

            panel_y += panel_len_y + this->gap_y;
        }
        panel_x += panel_len_x;
        if (i == this->num_panels_x / 2 - 1)
        {
            panel_x += this->gap_center;
        }
        else
        {
            panel_x += this->gap_x;
        }
    }

    // Absorber
    // TODO: Single element at the moment. Break up into multiple tubes.
    auto abs = make_element<SingleElement>();
    abs->set_name("Absorber");
    origin.set_values(0.0, 0.0,
                      this->receiver_height - 0.5 * this->abs_diameter);
    aim.set_values(0.0, 0.0, 1.0);
    vector_add(1.0, origin, 1.0, aim);
    abs->set_reference_frame_geometry(origin, aim, 0.0);
    abs->set_aperture(make_aperture<Rectangle>(this->abs_diameter,
                                               this->aperture_size_y));
    abs->set_surface(make_surface<Cylinder>(0.5 * this->abs_diameter));
    abs->set_front_optical_properties(this->optics_absorber);
    abs->set_back_optical_properties(this->optics_absorber);
    abs->enable();

    this->absorbers.push_back(abs);
    sts = this->add_element(abs);
    if (!Element::is_success(sts))
    {
        // TODO: Make this more helpful
        throw std::runtime_error("Failed to add absorber element");
    }

    // Envelope -- Outer
    auto envout = make_element<SingleElement>();
    envout->set_name("EnvelopeOuter");
    origin.set_values(0.0, 0.0,
                      this->receiver_height - 0.5 * this->env_diameter);
    aim.set_values(0.0, 0.0, 1.0);
    vector_add(1.0, origin, 1.0, aim);
    envout->set_reference_frame_geometry(origin, aim, 0.0);
    envout->set_aperture(make_aperture<Rectangle>(this->env_diameter,
                                                  this->aperture_size_y));
    envout->set_surface(make_surface<Cylinder>(0.5 * this->env_diameter));
    envout->set_front_optical_properties(this->optics_env_out);
    envout->set_back_optical_properties(this->optics_env_out);
    envout->enable();

    this->envelope.push_back(envout);
    sts = this->add_element(envout);
    if (!Element::is_success(sts))
    {
        throw std::runtime_error("Failed to add outer envelope element");
    }

    // Envelope -- Inner
    auto envin = make_element<SingleElement>();
    envin->set_name("EnvelopeInner");
    origin.set_values(0.0, 0.0,
                      this->receiver_height -
                          0.5 * this->env_diameter + this->env_thickness);
    aim.set_values(0.0, 0.0, 1.0);
    vector_add(1.0, origin, 1.0, aim);
    envin->set_reference_frame_geometry(origin, aim, 0.0);
    double ap_x = this->env_diameter - 2 * this->env_thickness;
    double ap_y = this->aperture_size_y;
    envin->set_aperture(make_aperture<Rectangle>(ap_x, ap_y));
    envin->set_surface(make_surface<Cylinder>(0.5 * ap_x));
    envin->set_front_optical_properties(this->optics_env_in);
    envin->set_back_optical_properties(this->optics_env_in);
    envin->enable();
    this->envelope.push_back(envin);
    sts = this->add_element(envin);
    if (!Element::is_success(sts))
    {
        throw std::runtime_error("Failed to add inner envelope element");
    }

    this->initialized = true;

    return;
}

void LinearFresnel::update_geometry(double azimuth, double elevation)
{
    if (elevation < 0.0 || elevation > 90.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::update_geometry: Invalid elevation ("
           << elevation << "). Elevation must lie between 0 and 90 degrees.";
        throw std::invalid_argument(ss.str());
    }

    if (azimuth < -180.0 || azimuth > 180.0)
    {
        std::stringstream ss;
        ss << "LinearFresnel::update_geometry: Invalid azimuth ("
           << elevation << "). Azimuth must lie between -180 and 180 degrees.";
        throw std::invalid_argument(ss.str());
    }

    if (!this->initialized)
    {
        std::stringstream ss;
        ss << "LinearFresnel::update_geometry: Uninitialized. "
           << "Call create_geometry() first.";
        throw std::invalid_argument(ss.str());
    }

    // NOTE: Setting the aim point and z-rotation only need to be done
    // once. But they need to be done after the LinearFresnel object
    // has been added to its reference element (if any) so we do it
    // here at the cost of repeating ourselves.

    // Set aim point
    Vector3d z_axis_ref, y_axis_ref;
    this->convert_global_to_reference(z_axis_ref, this->neutral_normal);
    this->convert_global_to_reference(y_axis_ref, this->rotation_axis);
    vector_add(1000.0, z_axis_ref, 1.0, this->origin, this->aim);

    // Set z-rotation
    double beta = asin(z_axis_ref[1]);
    double gamma = acos(y_axis_ref[1] / cos(beta));
    this->set_zrot_radians(gamma);

    // Update coordinate conversions so we can use them below
    this->coordinates_initialized = false;
    this->compute_coordinate_rotations();

    // std::cout << "Rotation axis: " << this->rotation_axis
    //           << "\nNeutral normal: " << this->neutral_normal
    //           << std::endl;

    // std::cout << "Global to Local: " << this->get_global_to_local()
    //           << std::endl;

    // std::cout << "z axis ref: " << z_axis_ref
    //           << "\ny axis ref: " << y_axis_ref
    //           << "\nBeta: " << beta
    //           << "\nGamma: " << gamma
    //           << std::endl;

    // Sun position projected into rotation plane and converted
    // to LinearFresnel object coordinates
    Vector3d sun_pos, sun_proj_local;
    sun_position_vector_degrees(sun_pos, azimuth, elevation);
    this->convert_global_to_local(sun_proj_local, sun_pos);
    // Project to rotation plane
    sun_proj_local[1] = 0.0;
    sun_proj_local.make_unit();

    // std::cout << "Sun Position: " << sun_pos
    //           //   << "\nSun Proj Global: " << sun_proj
    //           //   << "\nSun Proj Local: " << temp
    //           << "\nSun Proj Local: " << sun_proj_local
    //           << std::endl;

    // Set aimpoint for mirrors
    Vector3d aim_mirror_ref;
    // Absorber position projected into the plane of rotation (the xz-plane)
    Vector3d origin_abs_proj(0.0,
                             0.0,
                             this->receiver_height - 0.5 * this->abs_diameter);
    // origin_abs_proj.make_unit();
    for (auto iter : this->mirrors)
    {
        // Get mirror to to receiver vector
        vector_add(1.0, origin_abs_proj,
                   -1.0, iter->get_origin_ref(),
                   aim_mirror_ref);

        // Project onto the rotation plane
        aim_mirror_ref[1] = 0.0;
        aim_mirror_ref.make_unit();
        // std::cout << "Origin: " << iter->get_origin_ref()
        //           << "\nMirror to Receiver: " << aim_mirror_ref
        //           << std::endl;

        // Take bisector vector with the sun
        vector_add(1.0, sun_proj_local, 1.0, aim_mirror_ref);
        aim_mirror_ref.make_unit();
        // std::cout << "Sun Proj Local: " << sun_proj_local
        //           << "\nAim Mirror Ref: " << aim_mirror_ref
        //           << std::endl;

        // Dot product with [0, 0, 1]
        double theta = acos(aim_mirror_ref[2]) * R2D;
        // Dot product with [1, 0, 0]
        if (aim_mirror_ref[0] < 0)
            theta = -theta;
        // std::cout << "Theta: " << theta << std::endl;

        if (theta < this->tracking_limit_lower)
        {
            theta = this->tracking_limit_lower * D2R;
            aim_mirror_ref.set_values(sin(theta), 0.0, cos(theta));
        }
        else if (theta > this->tracking_limit_upper)
        {
            theta = this->tracking_limit_upper * D2R;
            aim_mirror_ref.set_values(sin(theta), 0.0, cos(theta));
        }

        // std::cout << "Theta 2: " << theta * R2D
        //           << "\nAim Mirror Ref: " << aim_mirror_ref
        //           << std::endl;

        // Add origin of mirror
        vector_add(1.0, iter->get_origin_ref(), 1000.0, aim_mirror_ref);
        // aim_mirror_ref.scalar_mult(1000.0);

        // Set aim point
        iter->set_aim_vector(aim_mirror_ref);
        iter->compute_coordinate_rotations();
    }

    return;
}

void LinearFresnel::enforce_user_fields_set() const
{
    if (this->initialized)
    {
        CompositeElement::enforce_user_fields_set();
    }

    // Validate that all required parameters have been set
    if (this->aperture_size_x <= 0.0 || this->aperture_size_y <= 0.0)
    {
        throw std::invalid_argument(
            "LinearFresnel: Aperture size must be set "
            "before creating geometry.");
    }

    if (this->num_panels_x <= 0 || this->num_panels_y <= 0)
    {
        throw std::invalid_argument(
            "LinearFresnel: Number of panels must be set "
            "before creating geometry.");
    }

    if (this->receiver_height <= 0.0)
    {
        throw std::invalid_argument(
            "LinearFresnel: Receiver height must be set "
            "before creating geometry.");
    }

    if (this->abs_diameter <= 0.0 ||
        this->env_diameter <= 0.0)
    {
        throw std::invalid_argument(
            "LinearFresnel: Receiver dimensions must be set "
            "before creating geometry.");
    }

    return;
}

} // namespace SolTrace::Data

1	#include "linear_fresnel.hpp"
2
3	#include <stdexcept>
4	#include <sstream>
5
6	#include "aperture.hpp"
7	#include "constants.hpp"
8	#include "element.hpp"
9	#include "surface.hpp"
10
11	#include "cst_templates/utilities.hpp"
12
13	namespace SolTrace::Data {
14
15	LinearFresnel::LinearFresnel()	15✔
16	: CompositeElement(),
17	initialized(false),	15✔
18	receiver_height(-1.0),	15✔
19	azimuth(-1.0),	15✔
20	tilt(-1.0),	15✔
21	focused_panels(false),	15✔
22	aperture_size_x(-1.0),	15✔
23	aperture_size_y(-1.0),	15✔
24	num_panels_x(-1),	15✔
25	num_panels_y(-1),	15✔
26	gap_x(-1.0),	15✔
27	gap_y(-1.0),	15✔
28	gap_center(-1.0),	15✔
29	abs_diameter(-1.0),	15✔
30	env_diameter(-1.0),	15✔
31	env_thickness(-1.0),	15✔
32	// tracking_angle(0.0),
33	tracking_limit_lower(-180.0),	15✔
34	tracking_limit_upper(180.0)	15✔
35	{
36	this->optics_absorber.set_ideal_absorption();	15✔
37	this->optics_mirror.set_ideal_reflection();	15✔
38	this->optics_env_out.set_ideal_transmission();	15✔
39	this->optics_env_in.set_ideal_transmission();	15✔
40
41	this->tracking_origin.set_values(1.0, 0.0, 0.0);	15✔
42	this->rotation_axis.set_values(0.0, 1.0, 0.0);	15✔
43	this->neutral_normal.set_values(0.0, 0.0, 1.0);	15✔
44	}	15✔
45
46	LinearFresnel::~LinearFresnel()	15✔
47	{
48	// Clear vectors
49	this->mirrors.clear();	15✔
50	this->absorbers.clear();	15✔
51	this->envelope.clear();	15✔
52
53	return;	15✔
54	}	15✔
55
56	// double LinearFresnel::get_tracking_angle_degrees() const
57	// {
58	// return this->tracking_angle;
59	// }
60
61	// double LinearFresnel::get_tracking_angle_radians() const
62	// {
63	// return D2R * this->get_tracking_angle_degrees();
64	// }
65
66	void LinearFresnel::set_angles(double azimuth, double tilt)	12✔
67	{
68	if (azimuth < -180.0 \|\| azimuth > 180.0)	12✔
69	{
70	std::stringstream ss;	2✔
71	ss << "LinearFresnel::set_angles: Invalid azimuth angle ("	2✔
72	<< azimuth << "). Must be between -180 and 180 degrees.";	2✔
73	throw std::invalid_argument(ss.str());	2✔
74	}	2✔
75
76	if (tilt < 0.0 \|\| tilt > 90.0)	10✔
77	{
78	std::stringstream ss;	2✔
79	ss << "LinearFresnel::set_angles: Invalid tilt angle ("	2✔
80	<< tilt << "). Must be between 0 and 90 degrees.";	2✔
81	throw std::invalid_argument(ss.str());	2✔
82	}	2✔
83
84	this->coordinates_initialized = false;	8✔
85	this->azimuth = azimuth;	8✔
86	this->tilt = tilt;	8✔
87
88	// Convert input angles to spherical coordinate angles
89	double az = azimuth * D2R;	8✔
90	double el = tilt * D2R;	8✔
91	double pol = 0.5 * PI - az;	8✔
92	double inc = 0.5 * PI - el;	8✔
93
94	// Convert spherical coordinates to cartesian coordinates
95	// y-axis
96	this->rotation_axis.set_values(sin(inc) * cos(pol),	8✔
97	sin(inc) * sin(pol),	8✔
98	cos(inc));
99
100	// z-axis
101	this->neutral_normal.set_values(sin(-el) * cos(pol),	8✔
102	sin(-el) * sin(pol),	8✔
103	cos(-el));
104
105	cross_product(this->rotation_axis,	8✔
106	this->neutral_normal,	8✔
107	this->tracking_origin);	8✔
108
109	return;	8✔
110	}
111
112	void LinearFresnel::set_aperture_size(double len_x, double len_y)	15✔
113	{
114	if (len_x <= 0.0 \|\| len_y <= 0.0)	15✔
115	{
116	std::stringstream ss;	5✔
117	ss << "LinearFresnel::set_aperture_size: Invalid aperture dimensions. "
118	<< "len_x (" << len_x << ") and len_y (" << len_y	5✔
119	<< ") must be positive.";	5✔
120	throw std::invalid_argument(ss.str());	5✔
121	}	5✔
122
123	this->initialized = false;	10✔
124	this->aperture_size_x = len_x;	10✔
125	this->aperture_size_y = len_y;	10✔
126
127	return;	10✔
128	}
129
130	void LinearFresnel::set_focused_panels(bool focused)	5✔
131	{
132	this->initialized = false;	5✔
133	this->focused_panels = focused;	5✔
134	return;	5✔
135	}
136
137	void LinearFresnel::set_gaps(double gap_x,	10✔
138	double gap_y,
139	double gap_center)
140	{
141	if (gap_x < 0.0 \|\| gap_y < 0.0 \|\| gap_center < 0.0)	10✔
142	{
143	std::stringstream ss;	3✔
144	ss << "LinearFresnel::set_gaps: Invalid gap dimensions. "
145	<< "gap_x (" << gap_x << "), gap_y (" << gap_y	3✔
146	<< "), and gap_center (" << gap_center	3✔
147	<< ") must be non-negative.";	3✔
148	throw std::invalid_argument(ss.str());	3✔
149	}	3✔
150
151	this->initialized = false;	7✔
152	this->gap_x = gap_x;	7✔
153	this->gap_y = gap_y;	7✔
154	this->gap_center = gap_center;	7✔
155
156	return;	7✔
157	}
158
159	void LinearFresnel::set_number_panels(int_fast64_t num_x,	14✔
160	int_fast64_t num_y)
161	{
162	if (num_x <= 0 \|\| num_y <= 0)	14✔
163	{
164	std::stringstream ss;	4✔
165	ss << "LinearFresnel::set_number_panels: Invalid panel count. "
166	<< "num_x (" << num_x << ") and num_y (" << num_y	4✔
167	<< ") must be positive.";	4✔
168	throw std::invalid_argument(ss.str());	4✔
169	}	4✔
170
171	this->initialized = false;	10✔
172	this->num_panels_x = num_x;	10✔
173	this->num_panels_y = num_y;	10✔
174
175	return;	10✔
176	}
177
178	void LinearFresnel::set_optics(const OpticalProperties &mirror,	4✔
179	const OpticalProperties &absorber,
180	const OpticalProperties &envelop_outer,
181	const OpticalProperties &envelop_inner)
182	{
183	this->optics_mirror = mirror;	4✔
184	this->optics_absorber = absorber;	4✔
185	this->optics_env_out = envelop_outer;	4✔
186	this->optics_env_in = envelop_inner;	4✔
187	return;	4✔
188	}
189
190	void LinearFresnel::set_receiver_height(double height)	12✔
191	{
192	if (height <= 0.0)	12✔
193	{
194	std::stringstream ss;	2✔
195	ss << "LinearFresnel::set_receiver_height: Invalid receiver height ("	2✔
196	<< height << "). Height must be positive.";	2✔
197	throw std::invalid_argument(ss.str());	2✔
198	}	2✔
199
200	this->initialized = false;	10✔
201	this->receiver_height = height;	10✔
202
203	return;	10✔
204	}
205
206	void LinearFresnel::set_receiver_dimensions(double absorber_diameter,	17✔
207	double envelop_diameter,
208	double envelop_thickness)
209	{
210	if (absorber_diameter <= 0.0)	17✔
211	{
212	std::stringstream ss;	2✔
213	ss << "LinearFresnel::set_receiver_dimensions: "
214	<< "Invalid absorber diameter (" << absorber_diameter	2✔
215	<< "). Must be positive.";	2✔
216	throw std::invalid_argument(ss.str());	2✔
217	}	2✔
218
219	if (envelop_diameter <= 0.0)	15✔
220	{
221	std::stringstream ss;	2✔
222	ss << "LinearFresnel::set_receiver_dimensions: "
223	<< "Invalid envelope diameter (" << envelop_diameter	2✔
224	<< "). Must be positive.";	2✔
225	throw std::invalid_argument(ss.str());	2✔
226	}	2✔
227
228	if (envelop_thickness < 0.0)	13✔
229	{
230	std::stringstream ss;	1✔
231	ss << "LinearFresnel::set_receiver_dimensions: "
232	<< "Invalid envelope thickness (" << envelop_thickness	1✔
233	<< "). Must be non-negative.";	1✔
234	throw std::invalid_argument(ss.str());	1✔
235	}	1✔
236
237	if (absorber_diameter >= envelop_diameter)	12✔
238	{
239	std::stringstream ss;	2✔
240	ss << "LinearFresnel::set_receiver_dimensions: "
241	<< "Absorber diameter (" << absorber_diameter	2✔
242	<< ") must be smaller than envelope diameter ("	2✔
243	<< envelop_diameter << ").";	2✔
244	throw std::invalid_argument(ss.str());	2✔
245	}	2✔
246
247	this->initialized = false;	10✔
248	this->abs_diameter = absorber_diameter;	10✔
249	this->env_diameter = envelop_diameter;	10✔
250	this->env_thickness = envelop_thickness;	10✔
251
252	return;	10✔
253	}
254
255	void LinearFresnel::set_tracking_limits(double lower, double upper)	4✔
256	{
257	if (lower > upper)	4✔
258	{
259	std::stringstream ss;	1✔
260	ss << "ParabolicTrough::set_tracking_limits: Invalid tracking "
261	<< "limits. Lower limit (" << lower	1✔
262	<< ") must be less than or equal to upper limit ("	1✔
263	<< upper << ").";	1✔
264	throw std::invalid_argument(ss.str());	1✔
265	}	1✔
266
267	this->tracking_limit_lower = lower;	3✔
268	this->tracking_limit_upper = upper;	3✔
269
270	return;	3✔
271	}
272
273	void LinearFresnel::create_geometry()	11✔
274	{
275	if (this->initialized)	11✔
276	{
NEW 277	return;	×
278	}
279
280	this->enforce_user_fields_set();	11✔
281
282	this->mirrors.clear();	7✔
283	this->absorbers.clear();	7✔
284	this->envelope.clear();	7✔
285	this->clear();	7✔
286
287	// Create mirrors
288	if (this->num_panels_x == 1 && this->gap_center != 0.0)	7✔
289	{
290	// TODO: This should generate a warning.
NEW 291	this->gap_center = 0.0;	×
NEW 292	this->gap_x = 0.0;	×
293	}
294
295	double panel_len_x = this->aperture_size_x -	7✔
296	this->gap_x * (this->num_panels_x - 1);	7✔
297	panel_len_x -= this->gap_center - this->gap_x;	7✔
298	panel_len_x /= this->num_panels_x;	7✔
299	double panel_x = -0.5 * this->aperture_size_x + 0.5 * panel_len_x;	7✔
300
301	double panel_len_y = this->aperture_size_y -	7✔
302	this->gap_y * (this->num_panels_y - 1);	7✔
303	panel_len_y /= this->num_panels_y;	7✔
304
305	element_id sts;
306	Vector3d origin;	7✔
307	Vector3d aim;	7✔
308	// TODO: Should mirrors aim at the receiver center or origin?
309	Vector3d receiver_pos(0.0,
310	0.0,
311	this->receiver_height - 0.5 * this->abs_diameter);	7✔
312	double panel_y, flen;
313	single_element_ptr mirror;	7✔
314	aperture_ptr ap;	7✔
315	surface_ptr surf;	7✔
316	Vector3d khat(0.0, 0.0, 1.0);	7✔
317
318	for (int_fast64_t i = 0; i < this->num_panels_x; ++i)	31✔
319	{
320	panel_y = -0.5 * this->aperture_size_y + 0.5 * panel_len_y;	24✔
321	for (int_fast64_t j = 0; j < this->num_panels_y; ++j)	89✔
322	{
323	mirror = make_element<SingleElement>();	65✔
324
325	origin.set_values(panel_x, panel_y, 0.0);	65✔
326	vector_add(1.0, receiver_pos, -1.0, origin, aim);	65✔
327	make_unit_vector(aim);	65✔
328	vector_add(0.5, khat, 0.5, aim);	65✔
329	vector_add(1.0, origin, 1.0, aim);	65✔
330	make_unit_vector(aim);	65✔
331	aim.scalar_mult(1000.0);	65✔
332	mirror->set_reference_frame_geometry(origin, aim, 0.0);	65✔
333
334	ap = make_aperture<Rectangle>(panel_len_x, panel_len_y);	65✔
335	mirror->set_aperture(ap);	65✔
336
337	if (this->focused_panels)	65✔
338	{
339	flen = sqrt(panel_x * panel_x +	48✔
340	this->receiver_height * this->receiver_height);	48✔
341	surf = make_surface<Parabola>(flen, 0.0);	48✔
342	}
343	else
344	{
345	surf = make_surface<Flat>();	17✔
346	}
347	mirror->set_surface(surf);	65✔
348
349	mirror->set_front_optical_properties(this->optics_mirror);	65✔
350	mirror->set_back_optical_properties(this->optics_mirror);	65✔
351	mirror->enable();	65✔
352
353	std::stringstream name;	65✔
354	name << "Mirror_" << this->num_panels_y * i + j;	65✔
355	mirror->set_name(name.str());	65✔
356
357	this->mirrors.push_back(mirror);	65✔
358	sts = this->add_element(mirror);	65✔
359	if (!Element::is_success(sts))	65✔
360	{
361	// TODO: Make this more helpful
NEW 362	throw std::runtime_error("Failed to add mirror element");	×
363	}
364
365	panel_y += panel_len_y + this->gap_y;	65✔
366	}	65✔
367	panel_x += panel_len_x;	24✔
368	if (i == this->num_panels_x / 2 - 1)	24✔
369	{
370	panel_x += this->gap_center;	6✔
371	}
372	else
373	{
374	panel_x += this->gap_x;	18✔
375	}
376	}
377
378	// Absorber
379	// TODO: Single element at the moment. Break up into multiple tubes.
380	auto abs = make_element<SingleElement>();	7✔
381	abs->set_name("Absorber");	14✔
382	origin.set_values(0.0, 0.0,	7✔
383	this->receiver_height - 0.5 * this->abs_diameter);	7✔
384	aim.set_values(0.0, 0.0, 1.0);	7✔
385	vector_add(1.0, origin, 1.0, aim);	7✔
386	abs->set_reference_frame_geometry(origin, aim, 0.0);	7✔
387	abs->set_aperture(make_aperture<Rectangle>(this->abs_diameter,	14✔
388	this->aperture_size_y));	7✔
389	abs->set_surface(make_surface<Cylinder>(0.5 * this->abs_diameter));	7✔
390	abs->set_front_optical_properties(this->optics_absorber);	7✔
391	abs->set_back_optical_properties(this->optics_absorber);	7✔
392	abs->enable();	7✔
393
394	this->absorbers.push_back(abs);	7✔
395	sts = this->add_element(abs);	7✔
396	if (!Element::is_success(sts))	7✔
397	{
398	// TODO: Make this more helpful
NEW 399	throw std::runtime_error("Failed to add absorber element");	×
400	}
401
402	// Envelope -- Outer
403	auto envout = make_element<SingleElement>();	7✔
404	envout->set_name("EnvelopeOuter");	14✔
405	origin.set_values(0.0, 0.0,	7✔
406	this->receiver_height - 0.5 * this->env_diameter);	7✔
407	aim.set_values(0.0, 0.0, 1.0);	7✔
408	vector_add(1.0, origin, 1.0, aim);	7✔
409	envout->set_reference_frame_geometry(origin, aim, 0.0);	7✔
410	envout->set_aperture(make_aperture<Rectangle>(this->env_diameter,	14✔
411	this->aperture_size_y));	7✔
412	envout->set_surface(make_surface<Cylinder>(0.5 * this->env_diameter));	7✔
413	envout->set_front_optical_properties(this->optics_env_out);	7✔
414	envout->set_back_optical_properties(this->optics_env_out);	7✔
415	envout->enable();	7✔
416
417	this->envelope.push_back(envout);	7✔
418	sts = this->add_element(envout);	7✔
419	if (!Element::is_success(sts))	7✔
420	{
NEW 421	throw std::runtime_error("Failed to add outer envelope element");	×
422	}
423
424	// Envelope -- Inner
425	auto envin = make_element<SingleElement>();	7✔
426	envin->set_name("EnvelopeInner");	14✔
427	origin.set_values(0.0, 0.0,	7✔
428	this->receiver_height -	7✔
429	0.5 * this->env_diameter + this->env_thickness);	7✔
430	aim.set_values(0.0, 0.0, 1.0);	7✔
431	vector_add(1.0, origin, 1.0, aim);	7✔
432	envin->set_reference_frame_geometry(origin, aim, 0.0);	7✔
433	double ap_x = this->env_diameter - 2 * this->env_thickness;	7✔
434	double ap_y = this->aperture_size_y;	7✔
435	envin->set_aperture(make_aperture<Rectangle>(ap_x, ap_y));	7✔
436	envin->set_surface(make_surface<Cylinder>(0.5 * ap_x));	7✔
437	envin->set_front_optical_properties(this->optics_env_in);	7✔
438	envin->set_back_optical_properties(this->optics_env_in);	7✔
439	envin->enable();	7✔
440	this->envelope.push_back(envin);	7✔
441	sts = this->add_element(envin);	7✔
442	if (!Element::is_success(sts))	7✔
443	{
NEW 444	throw std::runtime_error("Failed to add inner envelope element");	×
445	}
446
447	this->initialized = true;	7✔
448
449	return;	7✔
450	}	7✔
451
452	void LinearFresnel::update_geometry(double azimuth, double elevation)	8✔
453	{
454	if (elevation < 0.0 \|\| elevation > 90.0)	8✔
455	{
456	std::stringstream ss;	2✔
457	ss << "LinearFresnel::update_geometry: Invalid elevation ("	2✔
458	<< elevation << "). Elevation must lie between 0 and 90 degrees.";	2✔
459	throw std::invalid_argument(ss.str());	2✔
460	}	2✔
461
462	if (azimuth < -180.0 \|\| azimuth > 180.0)	6✔
463	{
464	std::stringstream ss;	2✔
465	ss << "LinearFresnel::update_geometry: Invalid azimuth ("	2✔
466	<< elevation << "). Azimuth must lie between -180 and 180 degrees.";	2✔
467	throw std::invalid_argument(ss.str());	2✔
468	}	2✔
469
470	if (!this->initialized)	4✔
471	{
472	std::stringstream ss;	1✔
473	ss << "LinearFresnel::update_geometry: Uninitialized. "
474	<< "Call create_geometry() first.";	1✔
475	throw std::invalid_argument(ss.str());	1✔
476	}	1✔
477
478	// NOTE: Setting the aim point and z-rotation only need to be done
479	// once. But they need to be done after the LinearFresnel object
480	// has been added to its reference element (if any) so we do it
481	// here at the cost of repeating ourselves.
482
483	// Set aim point
484	Vector3d z_axis_ref, y_axis_ref;	3✔
485	this->convert_global_to_reference(z_axis_ref, this->neutral_normal);	3✔
486	this->convert_global_to_reference(y_axis_ref, this->rotation_axis);	3✔
487	vector_add(1000.0, z_axis_ref, 1.0, this->origin, this->aim);	3✔
488
489	// Set z-rotation
490	double beta = asin(z_axis_ref[1]);	3✔
491	double gamma = acos(y_axis_ref[1] / cos(beta));	3✔
492	this->set_zrot_radians(gamma);	3✔
493
494	// Update coordinate conversions so we can use them below
495	this->coordinates_initialized = false;	3✔
496	this->compute_coordinate_rotations();	3✔
497
498	// std::cout << "Rotation axis: " << this->rotation_axis
499	// << "\nNeutral normal: " << this->neutral_normal
500	// << std::endl;
501
502	// std::cout << "Global to Local: " << this->get_global_to_local()
503	// << std::endl;
504
505	// std::cout << "z axis ref: " << z_axis_ref
506	// << "\ny axis ref: " << y_axis_ref
507	// << "\nBeta: " << beta
508	// << "\nGamma: " << gamma
509	// << std::endl;
510
511	// Sun position projected into rotation plane and converted
512	// to LinearFresnel object coordinates
513	Vector3d sun_pos, sun_proj_local;	3✔
514	sun_position_vector_degrees(sun_pos, azimuth, elevation);	3✔
515	this->convert_global_to_local(sun_proj_local, sun_pos);	3✔
516	// Project to rotation plane
517	sun_proj_local[1] = 0.0;	3✔
518	sun_proj_local.make_unit();	3✔
519
520	// std::cout << "Sun Position: " << sun_pos
521	// // << "\nSun Proj Global: " << sun_proj
522	// // << "\nSun Proj Local: " << temp
523	// << "\nSun Proj Local: " << sun_proj_local
524	// << std::endl;
525
526	// Set aimpoint for mirrors
527	Vector3d aim_mirror_ref;	3✔
528	// Absorber position projected into the plane of rotation (the xz-plane)
529	Vector3d origin_abs_proj(0.0,
530	0.0,
531	this->receiver_height - 0.5 * this->abs_diameter);	3✔
532	// origin_abs_proj.make_unit();
533	for (auto iter : this->mirrors)	11✔
534	{
535	// Get mirror to to receiver vector
536	vector_add(1.0, origin_abs_proj,	8✔
537	-1.0, iter->get_origin_ref(),	16✔
538	aim_mirror_ref);
539
540	// Project onto the rotation plane
541	aim_mirror_ref[1] = 0.0;	8✔
542	aim_mirror_ref.make_unit();	8✔
543	// std::cout << "Origin: " << iter->get_origin_ref()
544	// << "\nMirror to Receiver: " << aim_mirror_ref
545	// << std::endl;
546
547	// Take bisector vector with the sun
548	vector_add(1.0, sun_proj_local, 1.0, aim_mirror_ref);	8✔
549	aim_mirror_ref.make_unit();	8✔
550	// std::cout << "Sun Proj Local: " << sun_proj_local
551	// << "\nAim Mirror Ref: " << aim_mirror_ref
552	// << std::endl;
553
554	// Dot product with [0, 0, 1]
555	double theta = acos(aim_mirror_ref[2]) * R2D;	8✔
556	// Dot product with [1, 0, 0]
557	if (aim_mirror_ref[0] < 0)	8✔
558	theta = -theta;	6✔
559	// std::cout << "Theta: " << theta << std::endl;
560
561	if (theta < this->tracking_limit_lower)	8✔
562	{
563	theta = this->tracking_limit_lower * D2R;	6✔
564	aim_mirror_ref.set_values(sin(theta), 0.0, cos(theta));	6✔
565	}
566	else if (theta > this->tracking_limit_upper)	2✔
567	{
568	theta = this->tracking_limit_upper * D2R;	2✔
569	aim_mirror_ref.set_values(sin(theta), 0.0, cos(theta));	2✔
570	}
571
572	// std::cout << "Theta 2: " << theta * R2D
573	// << "\nAim Mirror Ref: " << aim_mirror_ref
574	// << std::endl;
575
576	// Add origin of mirror
577	vector_add(1.0, iter->get_origin_ref(), 1000.0, aim_mirror_ref);	8✔
578	// aim_mirror_ref.scalar_mult(1000.0);
579
580	// Set aim point
581	iter->set_aim_vector(aim_mirror_ref);	8✔
582	iter->compute_coordinate_rotations();	8✔
583	}	8✔
584
585	return;	6✔
586	}	3✔
587
588	void LinearFresnel::enforce_user_fields_set() const	13✔
589	{
590	if (this->initialized)	13✔
591	{
592	CompositeElement::enforce_user_fields_set();	2✔
593	}
594
595	// Validate that all required parameters have been set
596	if (this->aperture_size_x <= 0.0 \|\| this->aperture_size_y <= 0.0)	13✔
597	{
598	throw std::invalid_argument(
599	"LinearFresnel: Aperture size must be set "
600	"before creating geometry.");	1✔
601	}
602
603	if (this->num_panels_x <= 0 \|\| this->num_panels_y <= 0)	12✔
604	{
605	throw std::invalid_argument(
606	"LinearFresnel: Number of panels must be set "
607	"before creating geometry.");	1✔
608	}
609
610	if (this->receiver_height <= 0.0)	11✔
611	{
612	throw std::invalid_argument(
613	"LinearFresnel: Receiver height must be set "
614	"before creating geometry.");	1✔
615	}
616
617	if (this->abs_diameter <= 0.0 \|\|	10✔
618	this->env_diameter <= 0.0)	9✔
619	{
620	throw std::invalid_argument(
621	"LinearFresnel: Receiver dimensions must be set "
622	"before creating geometry.");	1✔
623	}
624
625	return;	9✔
626	}
627
628	} // namespace SolTrace::Data

NREL / SolTrace / 18544749821

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous