18984697163


#include "simdata_io.hpp"

#include <array>
#include <cstring>
#include <exception>
#include <sstream>
#include <string>

#include "constants.hpp"
#include "ray_source.hpp"
#include "simulation_data.hpp"
#include "single_element.hpp"
#include "stage_element.hpp"
#include "sun.hpp"
#include "surface.hpp"

namespace SolTrace::Data {

int st_sun_position(double lat, double day, double hour,
                    double *x, double *y, double *z)
{
    /*
    computes the sun vector xyz given arguments
    lat : [deg] latitude
    day : [] day of the year
    hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time

    xyz coordinate system:
        x: +west
        y: +zenith
        z: +north
    */

    double Declination, HourAngle, Elevation, Azimuth;
    // Use D2R and R2D from constants.hpp
    constexpr double deg_to_rad = D2R;
    constexpr double rad_to_deg = R2D;

    Declination = rad_to_deg * asin(0.39795 * cos(0.98563 * deg_to_rad * (day - 173)));
    HourAngle = 15 * (hour - 12);
    Elevation = rad_to_deg * asin(sin(Declination * deg_to_rad) * sin(lat * deg_to_rad) + cos(Declination * deg_to_rad) * cos(HourAngle * deg_to_rad) * cos(lat * deg_to_rad));
    Azimuth = rad_to_deg * acos((sin(deg_to_rad * Declination) * cos(deg_to_rad * lat) - cos(deg_to_rad * Declination) * sin(deg_to_rad * lat) * cos(deg_to_rad * HourAngle)) / cos(deg_to_rad * Elevation) + 0.0000000001);
    if (sin(HourAngle * deg_to_rad) > 0.0)
        Azimuth = 360 - Azimuth;
    *x = -sin(Azimuth * deg_to_rad) * cos(Elevation * deg_to_rad);
    *y = sin(Elevation * deg_to_rad);
    *z = cos(Azimuth * deg_to_rad) * cos(Elevation * deg_to_rad);

    return 1;
}

DistributionType char_to_distribution(const char dist_char)
{
    switch (dist_char)
    {
    case ('g'):
    {
        return DistributionType::GAUSSIAN;
    }
    case ('p'):
    {
        return DistributionType::PILLBOX;
    }
    case ('f'):
    {
        return DistributionType::DIFFUSE;
    }
    case ('d'):
    {
        return DistributionType::USER_DEFINED;
    }
    default:
    {
        return DistributionType::GAUSSIAN;
    }
    }
}

SunShape char_to_sunshape(const char dist_char)
{
    switch (dist_char)
    {
    case ('g'):
    {
        return SunShape::GAUSSIAN;
    }
    case ('p'):
    {
        return SunShape::PILLBOX;
    }
    case ('d'):
    {
        return SunShape::USER_DEFINED;
    }
    default:
    {
        return SunShape::GAUSSIAN;
    }
    }
}

InteractionType int_to_interaction(const int interaction_int)
{
    switch (interaction_int)
    {
    case (1):
    {
        return InteractionType::REFRACTION;
    }
    case (2):
    {
        return InteractionType::REFLECTION;
    }
    default:
    {
        return InteractionType::REFLECTION;
    }
    }
}

ApertureType char_to_aperture(const char aperture_char)
{
    switch (aperture_char)
    {
    case ('c'):
    {
        return ApertureType::CIRCLE;
    }
    case ('h'):
    {
        return ApertureType::HEXAGON;
    }
    case ('t'):
    {
        return ApertureType::EQUILATERAL_TRIANGLE;
    }
    case ('r'):
    {
        return ApertureType::RECTANGLE;
    }
    case ('a'):
    {
        return ApertureType::ANNULUS;
    }
    case ('l'):
    {
        return ApertureType::SINGLE_AXIS_CURVATURE_SECTION;
    }
    case ('i'):
    {
        return ApertureType::IRREGULAR_TRIANGLE;
    }
    case ('q'):
    {
        return ApertureType::IRREGULAR_QUADRILATERAL;
    }
    default:
    {
        return ApertureType::APERTURE_UNKNOWN;
    }
    }
}

SurfaceType char_to_surface(const char surface_char)
{
    switch (surface_char)
    {
    case ('s'):
        return SurfaceType::SPHERE;
    case ('p'):
        return SurfaceType::PARABOLA;
    case ('o'):
        return SurfaceType::HYPER;
    case ('g'):
        return SurfaceType::GENERAL_SPENCER_MURTY;
    case ('f'):
        return SurfaceType::FLAT;
    case ('c'):
        return SurfaceType::CONE;
    case ('t'):
        return SurfaceType::CYLINDER;
    case ('d'):
        return SurfaceType::TORUS;
    default:
        return SurfaceType::SURFACE_UNKNOWN;
    }
}

static void read_line(char *buf, int len, FILE *fp)
{
    fgets(buf, len, fp);
    int nch = strlen(buf);
    if (nch > 0 && buf[nch - 1] == '\n')
        buf[nch - 1] = 0;
    if (nch - 1 > 0 && buf[nch - 2] == '\r')
        buf[nch - 2] = 0;
}

std::vector<std::string> split(const std::string &str,
                               const std::string &delim,
                               bool ret_empty,
                               bool ret_delim)
{
    std::vector<std::string> list;

    char cur_delim[2] = {0, 0};
    std::string::size_type m_pos = 0;
    std::string token;

    while (m_pos < str.length())
    {
        std::string::size_type pos = str.find_first_of(delim, m_pos);
        if (pos == std::string::npos)
        {
            cur_delim[0] = 0;
            token.assign(str, m_pos, std::string::npos);
            m_pos = str.length();
        }
        else
        {
            cur_delim[0] = str[pos];
            std::string::size_type len = pos - m_pos;
            token.assign(str, m_pos, len);
            m_pos = pos + 1;
        }

        if (token.empty() && !ret_empty)
            continue;

        list.push_back(token);

        if (ret_delim && cur_delim[0] != 0 && m_pos < str.length())
            list.push_back(std::string(cur_delim));
    }

    return list;
}

bool process_sun(FILE *fp, SimulationData &sd)
{
    char buf[1024];

    // Read Sun info
    int bi = 0, count = 0;
    char cshape = 'g';
    double Sigma, HalfWidth;
    bool PointSource;
    double X, Y, Z, Latitude, Day, Hour;
    bool UseLDHSpec;

    read_line(buf, 1023, fp);
    sscanf(buf, "SUN\tPTSRC\t%d\tSHAPE\t%c\tSIGMA\t%lg\tHALFWIDTH\t%lg",
           &bi, &cshape, &Sigma, &HalfWidth);
    PointSource = (bi != 0);
    cshape = tolower(cshape);

    read_line(buf, 1023, fp);

    sscanf(buf, "XYZ\t%lg\t%lg\t%lg\tUSELDH\t%d\tLDH\t%lg\t%lg\t%lg",
           &X, &Y, &Z, &bi, &Latitude, &Day, &Hour);
    UseLDHSpec = (bi != 0);

    read_line(buf, 1023, fp);
    sscanf(buf, "USER SHAPE DATA\t%d", &count);
    std::vector<double> angle_vec;
    std::vector<double> intensity_vec;
    if (count > 0)
    {
        for (int i = 0; i < count; i++)
        {
            double x, y;
            read_line(buf, 1023, fp);
            sscanf(buf, "%lg\t%lg", &x, &y);
            angle_vec.push_back(x);
            intensity_vec.push_back(y);
        }
    }

    // Make sun
    auto sun = make_ray_source<Sun>();

    // Define sun position
    if (UseLDHSpec)
    {
        // sun->set_position(Latitude, Day, Hour);
        st_sun_position(Latitude, Day, Hour, &X, &Y, &Z);
    }
    sun->set_position(X, Y, Z);
    // else
    // {
    //  sun->set_position(X, Y, Z);
    // }

    // Define sun shape
    SunShape sun_shape = char_to_sunshape(cshape);
    sun->set_shape(sun_shape, Sigma, HalfWidth, 0.0, angle_vec, intensity_vec);
    // TOD: Buie sun shape not implemented here

    // TODO set point source

    // Attach sun to simulation data
    sd.add_ray_source(sun);
    return true;
}

bool read_optic_surface(FILE *fp,
                        OpticalProperties &optics,
                        int &OpticalSurfaceNumber)
{
    if (!fp)
        return false;
    char buf[1024];
    read_line(buf, 1023, fp);
    std::vector<std::string> parts = split(std::string(buf), "\t", true, false);
    if (parts.size() < 15)
    {
        printf("too few tokens for optical surface: %zu\n", parts.size());
        printf("\t>> %s\n", buf);
        return false;
    }

    char ErrorDistribution = 'g';
    if (parts[1].length() > 0)
        ErrorDistribution = parts[1][0];

    int ApertureStopOrGratingType = atoi(parts[2].c_str());
    OpticalSurfaceNumber = atoi(parts[3].c_str());
    int DiffractionOrder = atoi(parts[4].c_str());
    double Reflectivity = atof(parts[5].c_str());
    double Transmissivity = atof(parts[6].c_str());
    double RMSSlope = atof(parts[7].c_str());
    double RMSSpecularity = atof(parts[8].c_str());
    double RefractionIndexReal = atof(parts[9].c_str());
    double RefractionIndexImag = atof(parts[10].c_str());
    double GratingCoeffs[4];
    GratingCoeffs[0] = atof(parts[11].c_str());
    GratingCoeffs[1] = atof(parts[12].c_str());
    GratingCoeffs[2] = atof(parts[13].c_str());
    GratingCoeffs[3] = atof(parts[14].c_str());

    bool UseReflectivityTable = false;
    int refl_npoints = 0;
    double *refl_angles = 0;
    double *refls = 0;

    bool UseTransmissivityTable = false;
    int trans_npoints = 0;
    double *trans_angles = 0;
    double *transs = 0;

    if (parts.size() >= 17)
    {
        UseReflectivityTable = (atoi(parts[15].c_str()) > 0);
        refl_npoints = atoi(parts[16].c_str());
        if (parts.size() >= 19)
        {
            UseTransmissivityTable = (atoi(parts[17].c_str()) > 0);
            trans_npoints = atoi(parts[18].c_str());
        }
    }

    if (UseReflectivityTable)
    {
        refl_angles = new double[refl_npoints];
        refls = new double[refl_npoints];

        for (int i = 0; i < refl_npoints; i++)
        {
            read_line(buf, 1023, fp);
            sscanf(buf, "%lg %lg", &refl_angles[i], &refls[i]);
        }
    }
    if (UseTransmissivityTable)
    {
        trans_angles = new double[trans_npoints];
        transs = new double[trans_npoints];

        for (int i = 0; i < trans_npoints; i++)
        {
            read_line(buf, 1023, fp);
            sscanf(buf, "%lg %lg", &trans_angles[i], &transs[i]);
        }
    }

    // Define optical properties
    InteractionType interaction = InteractionType::REFLECTION;
    DistributionType dist = char_to_distribution(ErrorDistribution);
    optics = OpticalProperties(interaction, dist, Transmissivity,
                               Reflectivity, RMSSlope, RMSSpecularity,
                               RefractionIndexReal, RefractionIndexImag);

    if (refl_angles != 0)
        delete[] refl_angles;
    if (refls != 0)
        delete[] refls;
    if (trans_angles != 0)
        delete[] trans_angles;
    if (transs != 0)
        delete[] transs;
    return true;
}

bool process_optics(
    FILE *fp,
    std::map<std::string, std::array<OpticalProperties, 2>> &optics_map)
{
    char buf[1024];

    // Read number of optics
    int count = 0;
    read_line(buf, 1023, fp);
    sscanf(buf, "OPTICS LIST COUNT\t%d", &count);

    // Define each optics
    for (int i = 0; i < count; i++)
    {
        // Read optical pair info line
        read_line(buf, 1023, fp);

        if (strncmp(buf, "OPTICAL PAIR", 12) == 0)
        {
            // int iopt = st_add_optic(cxt, (const char*)(buf + 13));
            std::string optics_name = std::string(buf + 13);
            OpticalProperties optics_front, optics_back;
            int OpticalSurfaceNumber = 0;
            read_optic_surface(fp, optics_front, OpticalSurfaceNumber);
            read_optic_surface(fp, optics_back, OpticalSurfaceNumber);

            optics_map[optics_name][0] = optics_front;
            optics_map[optics_name][1] = optics_back;
        }
        else
            return false;
    }

    return true;
}

bool read_element(
    FILE *fp,
    std::map<std::string, std::array<OpticalProperties, 2>> &optics_map,
    element_ptr &el)
{
    char buf[1024];
    read_line(buf, 1023, fp);

    std::vector<std::string> tok = split(buf, "\t", true, false);
    if (tok.size() < 29)
    {
        printf("too few tokens for element: %zu\n", tok.size());
        printf("\t>> %s\n", buf);
        return false;
    }

    bool enabled = atoi(tok[0].c_str()) ? 1 : 0;
    double xyz[3] = {
        atof(tok[1].c_str()),
        atof(tok[2].c_str()),
        atof(tok[3].c_str())};
    double aim[3] = {
        atof(tok[4].c_str()),
        atof(tok[5].c_str()),
        atof(tok[6].c_str())};
    double zrot = atof(tok[7].c_str());

    char ShapeIndex = ' ';
    if (tok[8].length() > 0)
    {
        ShapeIndex = tok[8][0];
    }
    else
    {
        printf("no aperture index specified for element\n");
        return false;
    }

    std::vector<double> aperture_params;
    for (int i = 0; i < 8; i++)
    {
        aperture_params.push_back(atof(tok[i + 9].c_str()));
    }

    char SurfaceIndex = ' ';
    if (tok[17].length() > 0)
    {
        SurfaceIndex = tok[17][0];
    }
    else
    {
        printf("no surface index specified for element\n");
        return false;
    }

    std::vector<double> surface_params;
    for (int i = 0; i < 8; i++)
    {
        surface_params.push_back(atof(tok[i + 18].c_str()));
    }

    // Skipping surface file for now
    std::string SurfaceFile = tok[26];
    std::string optics_name = tok[27].c_str();
    InteractionType interaction = int_to_interaction(atoi(tok[28].c_str()));

    // Create element
    el = make_element<SingleElement>();

    // Make aperture
    ApertureType aperture_type = char_to_aperture(ShapeIndex);
    if (aperture_type == ApertureType::APERTURE_UNKNOWN)
    {
        std::stringstream ss;
        ss << "Aperture character " << ShapeIndex
           << " returned unknown aperture type " << aperture_type;
        throw std::invalid_argument(ss.str());
    }
    
    aperture_ptr ap_ptr = Aperture::make_aperture_from_type(
        aperture_type, aperture_params);
    if (ap_ptr == nullptr)
    {
        std::stringstream ss;
        ss << "Unable to make aperture pointer -- "
           << "\nChar: " << ShapeIndex
           << "\nType: " << aperture_type
           << "\nParams: [";
        for (auto cit = aperture_params.cbegin();
             cit != aperture_params.cend();
             ++cit)
        {
            ss << *cit << ", ";
        }
        ss << "]" << std::endl;
        throw std::runtime_error(ss.str());
    }
    el->set_aperture(ap_ptr);

    // Make surface
    SurfaceType surface_type = char_to_surface(SurfaceIndex);
    if (surface_type == SurfaceType::SURFACE_UNKNOWN)
    {
        std::stringstream ss;
        ss << "Unknown surface type " << surface_type;
        throw std::invalid_argument(ss.str());
    }
    surface_ptr surf_ptr = make_surface_from_type(surface_type, surface_params);
    el->set_surface(surf_ptr);
    if (surface_type == SurfaceType::CYLINDER)
    {
        ap_ptr = el->get_aperture();
        auto rect = std::dynamic_pointer_cast<Rectangle>(ap_ptr);
        auto cyl = std::dynamic_pointer_cast<Cylinder>(surf_ptr);
        if (rect == nullptr || cyl == nullptr)
        {
            throw std::invalid_argument("This should not happen!");
        }
        rect->x_length = 2.0 * cyl->radius;
    }

    // Set element position and orientation
    el->set_reference_frame_geometry(Vector3d(xyz),
                                     Vector3d(aim),
                                     zrot);

    // Set optical properties
    OpticalProperties optics_front = optics_map[optics_name][0];
    OpticalProperties optics_back = optics_map[optics_name][1];
    el->set_front_optical_properties(optics_front);
    el->set_back_optical_properties(optics_back);

    // Set optical interaction type
    el->get_front_optical_properties()->my_type = interaction;
    el->get_back_optical_properties()->my_type = interaction;

    return true;
}

bool process_stages(
    FILE *fp,
    SimulationData &sd,
    std::map<std::string, std::array<OpticalProperties, 2>> &optics_map)
{
    char buf[1024];

    // Loop through stages
    int count_stage = 0;
    read_line(buf, 1023, fp);
    sscanf(buf, "STAGE LIST COUNT\t%d", &count_stage);

    for (int i_stage = 0; i_stage < count_stage; i_stage++)
    {
        int virt = 0, multi = 1, count_element = 0, tr = 0;
        double X, Y, Z, AX, AY, AZ, ZRot;

        read_line(buf, 1023, fp);
        sscanf(buf, "STAGE\tXYZ\t%lg\t%lg\t%lg\tAIM\t%lg\t%lg\t%lg\tZROT\t%lg\tVIRTUAL\t%d\tMULTIHIT\t%d\tELEMENTS\t%d\tTRACETHROUGH\t%d",
               &X, &Y, &Z,
               &AX, &AY, &AZ,
               &ZRot,
               &virt,
               &multi,
               &count_element,
               &tr);

        read_line(buf, 1023, fp); // read name

        // Make stage
        stage_ptr stage = make_stage(i_stage);
        stage->set_origin(X, Y, Z);
        stage->set_aim_vector(AX, AY, AZ);
        stage->set_zrot(ZRot);
        stage->compute_coordinate_rotations();

        // Loop through elements
        for (int i_element = 0; i_element < count_element; i_element++)
        {
            element_ptr el;
            read_element(fp, optics_map, el);
            el->set_name(std::to_string(i_element));
            // TODO make virtual if stage is virtual?

            if (!Element::is_success(stage->add_element(el)))
            {
                std::cout << "Failed to add element to stage" << std::endl;
            }
        }

        if (!Element::is_success(sd.add_stage(stage)))
        {
            std::cout << "Failed to add stage to SimulationData" << std::endl;
        }
    }

    return true;
}

bool process_sim_par(FILE *fp, SimulationData &sd)
{
    char buf[1024];

    // Check if end of file
    if (feof(fp))
        return false;

    // Check for simulation parameters
    read_line(buf, 1023, fp);
    if (strncmp(buf, "TRACE", 5) != 0)
        return false;

    // Get simulation parameters
    int n_rays, n_rays_sun, n_cpu, seed, ss, err, pf;
    int n = sscanf(buf, "TRACE\tNRAY\t%d\tNSUN\t%d\tCPU\t%d\tSEED\t%d\tSUNSHAPE\t%d\tERRORS\t%d\tPTFOCUS\t%d",
                   &n_rays, &n_rays_sun, &n_cpu, &seed, &ss, &err, &pf);

    // Assign simulation parameters
    SimulationParameters &par = sd.get_simulation_parameters();
    par.number_of_rays = n_rays;
    par.max_number_of_rays = n_rays_sun;
    par.seed = seed;
    par.include_sun_shape_errors = ss;
    par.include_optical_errors = err;

    // TODO Assign number CPUs, point focus?
    return true;
}

bool load_stinput_file(SimulationData &sd, std::string filename)
{
    // TODO: Reset simulation data?

    // Read in file
    FILE *fp = fopen(filename.data(), "r");
    if (!fp)
    {
        printf("failed to open system input file: %s\n",
               filename.data());
        return false;
    }

    // Buffer to store read line
    char buf[1024];

    // Get version info (if first line starts with '#')
    int vmaj = 0, vmin = 0, vmic = 0;
    char c = fgetc(fp);
    if (c == '#')
    {
        read_line(buf, 1023, fp);
        sscanf(buf, " SOLTRACE VERSION %d.%d.%d INPUT FILE",
               &vmaj, &vmin, &vmic);

        // unsigned int file_version = vmaj*10000 + vmin*100 + vmic;

        printf("loading input file version %d.%d.%d\n",
               vmaj, vmin, vmic);
    }
    else
    {
        ungetc(c, fp);
    }

    // Read in Sun
    process_sun(fp, sd);

    // Read in Optics
    std::map<std::string, std::array<OpticalProperties, 2>> optics_map;
    process_optics(fp, optics_map);

    // Read in Stages
    process_stages(fp, sd, optics_map);

    // Read in simulation parameters (if any)
    process_sim_par(fp, sd);

    return true;
}

} // namespace SolTrace::Data

1
2	#include "simdata_io.hpp"
3
4	#include <array>
5	#include <cstring>
6	#include <exception>
7	#include <sstream>
8	#include <string>
9
10	#include "constants.hpp"
11	#include "ray_source.hpp"
12	#include "simulation_data.hpp"
13	#include "single_element.hpp"
14	#include "stage_element.hpp"
15	#include "sun.hpp"
16	#include "surface.hpp"
17
18	namespace SolTrace::Data {
19
20	int st_sun_position(double lat, double day, double hour,	1✔
21	double x, double y, double *z)
22	{
23	/*
24	computes the sun vector xyz given arguments
25	lat : [deg] latitude
26	day : [] day of the year
27	hour : [hour] solar time. 12.00 corresponds to sun at maximum elevation and does not necessarily match local time
28
29	xyz coordinate system:
30	x: +west
31	y: +zenith
32	z: +north
33	*/
34
35	double Declination, HourAngle, Elevation, Azimuth;
36	// Use D2R and R2D from constants.hpp
37	constexpr double deg_to_rad = D2R;	1✔
38	constexpr double rad_to_deg = R2D;	1✔
39
40	Declination = rad_to_deg * asin(0.39795 * cos(0.98563 * deg_to_rad * (day - 173)));	1✔
41	HourAngle = 15 * (hour - 12);	1✔
42	Elevation = rad_to_deg * asin(sin(Declination * deg_to_rad) * sin(lat * deg_to_rad) + cos(Declination * deg_to_rad) * cos(HourAngle * deg_to_rad) * cos(lat * deg_to_rad));	1✔
43	Azimuth = rad_to_deg * acos((sin(deg_to_rad * Declination) * cos(deg_to_rad * lat) - cos(deg_to_rad * Declination) * sin(deg_to_rad * lat) * cos(deg_to_rad * HourAngle)) / cos(deg_to_rad * Elevation) + 0.0000000001);	1✔
44	if (sin(HourAngle * deg_to_rad) > 0.0)	1✔
45	Azimuth = 360 - Azimuth;	×
46	x = -sin(Azimuth deg_to_rad) * cos(Elevation * deg_to_rad);	1✔
47	y = sin(Elevation deg_to_rad);	1✔
48	z = cos(Azimuth deg_to_rad) * cos(Elevation * deg_to_rad);	1✔
49
50	return 1;	1✔
51	}
52
53	DistributionType char_to_distribution(const char dist_char)	10✔
54	{
55	switch (dist_char)	10✔
56	{
57	case ('g'):	10✔
58	{
59	return DistributionType::GAUSSIAN;	10✔
60	}
UNCOV 61	case ('p'):	×
62	{
UNCOV 63	return DistributionType::PILLBOX;	×
64	}
NEW 65	case ('f'):	×
66	{
NEW 67	return DistributionType::DIFFUSE;	×
68	}
UNCOV 69	case ('d'):	×
70	{
UNCOV 71	return DistributionType::USER_DEFINED;	×
72	}
73	default:	×
74	{
75	return DistributionType::GAUSSIAN;	×
76	}
77	}
78	}
79
80	SunShape char_to_sunshape(const char dist_char)	3✔
81	{
82	switch (dist_char)	3✔
83	{
NEW 84	case ('g'):	×
85	{
NEW 86	return SunShape::GAUSSIAN;	×
87	}
88	case ('p'):	2✔
89	{
90	return SunShape::PILLBOX;	2✔
91	}
92	case ('d'):	1✔
93	{
94	return SunShape::USER_DEFINED;	1✔
95	}
NEW 96	default:	×
97	{
NEW 98	return SunShape::GAUSSIAN;	×
99	}
100	}
101	}
102
103	InteractionType int_to_interaction(const int interaction_int)	6,312✔
104	{
105	switch (interaction_int)	6,312✔
106	{
107	case (1):	×
108	{
109	return InteractionType::REFRACTION;	×
110	}
111	case (2):	6,312✔
112	{
113	return InteractionType::REFLECTION;	6,312✔
114	}
115	default:	×
116	{
117	return InteractionType::REFLECTION;	×
118	}
119	}
120	}
121
122	ApertureType char_to_aperture(const char aperture_char)	6,312✔
123	{
124	switch (aperture_char)	6,312✔
125	{
126	case ('c'):	×
127	{
128	return ApertureType::CIRCLE;	×
129	}
130	case ('h'):	26✔
131	{
132	return ApertureType::HEXAGON;	26✔
133	}
134	case ('t'):	×
135	{
136	return ApertureType::EQUILATERAL_TRIANGLE;	×
137	}
138	case ('r'):	6,285✔
139	{
140	return ApertureType::RECTANGLE;	6,285✔
141	}
142	case ('a'):	×
143	{
144	return ApertureType::ANNULUS;	×
145	}
146	case ('l'):	1✔
147	{
148	return ApertureType::SINGLE_AXIS_CURVATURE_SECTION;	1✔
149	}
150	case ('i'):	×
151	{
152	return ApertureType::IRREGULAR_TRIANGLE;	×
153	}
154	case ('q'):	×
155	{
156	return ApertureType::IRREGULAR_QUADRILATERAL;	×
157	}
158	default:	×
159	{
160	return ApertureType::APERTURE_UNKNOWN;	×
161	}
162	}
163	}
164
165	SurfaceType char_to_surface(const char surface_char)	6,312✔
166	{
167	switch (surface_char)	6,312✔
168	{
169	case ('s'):	26✔
170	return SurfaceType::SPHERE;	26✔
171	case ('p'):	6,283✔
172	return SurfaceType::PARABOLA;	6,283✔
173	case ('o'):	×
174	return SurfaceType::HYPER;	×
175	case ('g'):	×
176	return SurfaceType::GENERAL_SPENCER_MURTY;	×
177	case ('f'):	2✔
178	return SurfaceType::FLAT;	2✔
179	case ('c'):	×
180	return SurfaceType::CONE;	×
181	case ('t'):	1✔
182	return SurfaceType::CYLINDER;	1✔
183	case ('d'):	×
184	return SurfaceType::TORUS;	×
185	default:	×
186	return SurfaceType::SURFACE_UNKNOWN;	×
187	}
188	}
189
190	static void read_line(char buf, int len, FILE fp)	6,388✔
191	{
192	fgets(buf, len, fp);	6,388✔
193	int nch = strlen(buf);	6,388✔
194	if (nch > 0 && buf[nch - 1] == '\n')	6,388✔
195	buf[nch - 1] = 0;	6,386✔
196	if (nch - 1 > 0 && buf[nch - 2] == '\r')	6,388✔
197	buf[nch - 2] = 0;	×
198	}	6,388✔
199
200	std::vector<std::string> split(const std::string &str,	6,322✔
201	const std::string &delim,
202	bool ret_empty,
203	bool ret_delim)
204	{
205	std::vector<std::string> list;	6,322✔
206
207	char cur_delim[2] = {0, 0};	6,322✔
208	std::string::size_type m_pos = 0;	6,322✔
209	std::string token;	6,322✔
210
211	while (m_pos < str.length())	189,528✔
212	{
213	std::string::size_type pos = str.find_first_of(delim, m_pos);	183,206✔
214	if (pos == std::string::npos)	183,206✔
215	{
216	cur_delim[0] = 0;	10✔
217	token.assign(str, m_pos, std::string::npos);	10✔
218	m_pos = str.length();	10✔
219	}
220	else
221	{
222	cur_delim[0] = str[pos];	183,196✔
223	std::string::size_type len = pos - m_pos;	183,196✔
224	token.assign(str, m_pos, len);	183,196✔
225	m_pos = pos + 1;	183,196✔
226	}
227
228	if (token.empty() && !ret_empty)	183,206✔
229	continue;	×
230
231	list.push_back(token);	183,206✔
232
233	if (ret_delim && cur_delim[0] != 0 && m_pos < str.length())	183,206✔
234	list.push_back(std::string(cur_delim));	×
235	}
236
237	return list;	12,644✔
238	}	6,322✔
239
240	bool process_sun(FILE *fp, SimulationData &sd)	3✔
241	{
242	char buf[1024];
243
244	// Read Sun info
245	int bi = 0, count = 0;	3✔
246	char cshape = 'g';	3✔
247	double Sigma, HalfWidth;
248	bool PointSource;
249	double X, Y, Z, Latitude, Day, Hour;
250	bool UseLDHSpec;
251
252	read_line(buf, 1023, fp);	3✔
253	sscanf(buf, "SUN\tPTSRC\t%d\tSHAPE\t%c\tSIGMA\t%lg\tHALFWIDTH\t%lg",	3✔
254	&bi, &cshape, &Sigma, &HalfWidth);
255	PointSource = (bi != 0);	3✔
256	cshape = tolower(cshape);	3✔
257
258	read_line(buf, 1023, fp);	3✔
259
260	sscanf(buf, "XYZ\t%lg\t%lg\t%lg\tUSELDH\t%d\tLDH\t%lg\t%lg\t%lg",	3✔
261	&X, &Y, &Z, &bi, &Latitude, &Day, &Hour);
262	UseLDHSpec = (bi != 0);	3✔
263
264	read_line(buf, 1023, fp);	3✔
265	sscanf(buf, "USER SHAPE DATA\t%d", &count);	3✔
266	std::vector<double> angle_vec;	3✔
267	std::vector<double> intensity_vec;	3✔
268	if (count > 0)	3✔
269	{
270	for (int i = 0; i < count; i++)	30✔
271	{
272	double x, y;
273	read_line(buf, 1023, fp);	28✔
274	sscanf(buf, "%lg\t%lg", &x, &y);	28✔
275	angle_vec.push_back(x);	28✔
276	intensity_vec.push_back(y);	28✔
277	}
278	}
279
280	// Make sun
281	auto sun = make_ray_source<Sun>();	3✔
282
283	// Define sun position
284	if (UseLDHSpec)	3✔
285	{
286	// sun->set_position(Latitude, Day, Hour);
287	st_sun_position(Latitude, Day, Hour, &X, &Y, &Z);	1✔
288	}
289	sun->set_position(X, Y, Z);	3✔
290	// else
291	// {
292	// sun->set_position(X, Y, Z);
293	// }
294
295	// Define sun shape
296	SunShape sun_shape = char_to_sunshape(cshape);	3✔
297	sun->set_shape(sun_shape, Sigma, HalfWidth, 0.0, angle_vec, intensity_vec);	3✔
298	// TOD: Buie sun shape not implemented here
299
300	// TODO set point source
301
302	// Attach sun to simulation data
303	sd.add_ray_source(sun);	3✔
304	return true;	3✔
305	}	3✔
306
307	bool read_optic_surface(FILE *fp,	10✔
308	OpticalProperties &optics,
309	int &OpticalSurfaceNumber)
310	{
311	if (!fp)	10✔
312	return false;	×
313	char buf[1024];
314	read_line(buf, 1023, fp);	10✔
315	std::vector<std::string> parts = split(std::string(buf), "\t", true, false);	30✔
316	if (parts.size() < 15)	10✔
317	{
318	printf("too few tokens for optical surface: %zu\n", parts.size());	×
319	printf("\t>> %s\n", buf);	×
320	return false;	×
321	}
322
323	char ErrorDistribution = 'g';	10✔
324	if (parts[1].length() > 0)	10✔
325	ErrorDistribution = parts[1][0];	10✔
326
327	int ApertureStopOrGratingType = atoi(parts[2].c_str());	10✔
328	OpticalSurfaceNumber = atoi(parts[3].c_str());	10✔
329	int DiffractionOrder = atoi(parts[4].c_str());	10✔
330	double Reflectivity = atof(parts[5].c_str());	10✔
331	double Transmissivity = atof(parts[6].c_str());	10✔
332	double RMSSlope = atof(parts[7].c_str());	10✔
333	double RMSSpecularity = atof(parts[8].c_str());	10✔
334	double RefractionIndexReal = atof(parts[9].c_str());	10✔
335	double RefractionIndexImag = atof(parts[10].c_str());	10✔
336	double GratingCoeffs[4];
337	GratingCoeffs[0] = atof(parts[11].c_str());	10✔
338	GratingCoeffs[1] = atof(parts[12].c_str());	10✔
339	GratingCoeffs[2] = atof(parts[13].c_str());	10✔
340	GratingCoeffs[3] = atof(parts[14].c_str());	10✔
341
342	bool UseReflectivityTable = false;	10✔
343	int refl_npoints = 0;	10✔
344	double *refl_angles = 0;	10✔
345	double *refls = 0;	10✔
346
347	bool UseTransmissivityTable = false;	10✔
348	int trans_npoints = 0;	10✔
349	double *trans_angles = 0;	10✔
350	double *transs = 0;	10✔
351
352	if (parts.size() >= 17)	10✔
353	{
354	UseReflectivityTable = (atoi(parts[15].c_str()) > 0);	2✔
355	refl_npoints = atoi(parts[16].c_str());	2✔
356	if (parts.size() >= 19)	2✔
357	{
358	UseTransmissivityTable = (atoi(parts[17].c_str()) > 0);	2✔
359	trans_npoints = atoi(parts[18].c_str());	2✔
360	}
361	}
362
363	if (UseReflectivityTable)	10✔
364	{
365	refl_angles = new double[refl_npoints];	×
366	refls = new double[refl_npoints];	×
367
368	for (int i = 0; i < refl_npoints; i++)	×
369	{
370	read_line(buf, 1023, fp);	×
371	sscanf(buf, "%lg %lg", &refl_angles[i], &refls[i]);	×
372	}
373	}
374	if (UseTransmissivityTable)	10✔
375	{
376	trans_angles = new double[trans_npoints];	×
377	transs = new double[trans_npoints];	×
378
379	for (int i = 0; i < trans_npoints; i++)	×
380	{
381	read_line(buf, 1023, fp);	×
382	sscanf(buf, "%lg %lg", &trans_angles[i], &transs[i]);	×
383	}
384	}
385
386	// Define optical properties
387	InteractionType interaction = InteractionType::REFLECTION;	10✔
388	DistributionType dist = char_to_distribution(ErrorDistribution);	10✔
389	optics = OpticalProperties(interaction, dist, Transmissivity,	10✔
390	Reflectivity, RMSSlope, RMSSpecularity,
391	RefractionIndexReal, RefractionIndexImag);
392
393	if (refl_angles != 0)	10✔
394	delete[] refl_angles;	×
395	if (refls != 0)	10✔
396	delete[] refls;	×
397	if (trans_angles != 0)	10✔
398	delete[] trans_angles;	×
399	if (transs != 0)	10✔
400	delete[] transs;	×
401	return true;	10✔
402	}	10✔
403
404	bool process_optics(	3✔
405	FILE *fp,
406	std::map<std::string, std::array<OpticalProperties, 2>> &optics_map)
407	{
408	char buf[1024];
409
410	// Read number of optics
411	int count = 0;	3✔
412	read_line(buf, 1023, fp);	3✔
413	sscanf(buf, "OPTICS LIST COUNT\t%d", &count);	3✔
414
415	// Define each optics
416	for (int i = 0; i < count; i++)	8✔
417	{
418	// Read optical pair info line
419	read_line(buf, 1023, fp);	5✔
420
421	if (strncmp(buf, "OPTICAL PAIR", 12) == 0)	5✔
422	{
423	// int iopt = st_add_optic(cxt, (const char*)(buf + 13));
424	std::string optics_name = std::string(buf + 13);	5✔
425	OpticalProperties optics_front, optics_back;	5✔
426	int OpticalSurfaceNumber = 0;	5✔
427	read_optic_surface(fp, optics_front, OpticalSurfaceNumber);	5✔
428	read_optic_surface(fp, optics_back, OpticalSurfaceNumber);	5✔
429
430	optics_map[optics_name][0] = optics_front;	5✔
431	optics_map[optics_name][1] = optics_back;	5✔
432	}	5✔
433	else
434	return false;	×
435	}
436
437	return true;	3✔
438	}
439
440	bool read_element(	6,312✔
441	FILE *fp,
442	std::map<std::string, std::array<OpticalProperties, 2>> &optics_map,
443	element_ptr &el)
444	{
445	char buf[1024];
446	read_line(buf, 1023, fp);	6,312✔
447
448	std::vector<std::string> tok = split(buf, "\t", true, false);	18,936✔
449	if (tok.size() < 29)	6,312✔
450	{
451	printf("too few tokens for element: %zu\n", tok.size());	×
452	printf("\t>> %s\n", buf);	×
453	return false;	×
454	}
455
456	bool enabled = atoi(tok[0].c_str()) ? 1 : 0;	6,312✔
457	double xyz[3] = {
458	atof(tok[1].c_str()),	6,312✔
459	atof(tok[2].c_str()),	6,312✔
460	atof(tok[3].c_str())};	12,624✔
461	double aim[3] = {
462	atof(tok[4].c_str()),	6,312✔
463	atof(tok[5].c_str()),	6,312✔
464	atof(tok[6].c_str())};	12,624✔
465	double zrot = atof(tok[7].c_str());	6,312✔
466
467	char ShapeIndex = ' ';	6,312✔
468	if (tok[8].length() > 0)	6,312✔
469	{
470	ShapeIndex = tok[8][0];	6,312✔
471	}
472	else
473	{
474	printf("no aperture index specified for element\n");	×
475	return false;	×
476	}
477
478	std::vector<double> aperture_params;	6,312✔
479	for (int i = 0; i < 8; i++)	56,808✔
480	{
481	aperture_params.push_back(atof(tok[i + 9].c_str()));	50,496✔
482	}
483
484	char SurfaceIndex = ' ';	6,312✔
485	if (tok[17].length() > 0)	6,312✔
486	{
487	SurfaceIndex = tok[17][0];	6,312✔
488	}
489	else
490	{
491	printf("no surface index specified for element\n");	×
492	return false;	×
493	}
494
495	std::vector<double> surface_params;	6,312✔
496	for (int i = 0; i < 8; i++)	56,808✔
497	{
498	surface_params.push_back(atof(tok[i + 18].c_str()));	50,496✔
499	}
500
501	// Skipping surface file for now
502	std::string SurfaceFile = tok[26];	6,312✔
503	std::string optics_name = tok[27].c_str();	6,312✔
504	InteractionType interaction = int_to_interaction(atoi(tok[28].c_str()));	6,312✔
505
506	// Create element
507	el = make_element<SingleElement>();	6,312✔
508
509	// Make aperture
510	ApertureType aperture_type = char_to_aperture(ShapeIndex);	6,312✔
511	if (aperture_type == ApertureType::APERTURE_UNKNOWN)	6,312✔
512	{
513	std::stringstream ss;	×
514	ss << "Aperture character " << ShapeIndex
515	<< " returned unknown aperture type " << aperture_type;	×
516	throw std::invalid_argument(ss.str());	×
517	}	×
518
519	aperture_ptr ap_ptr = Aperture::make_aperture_from_type(
520	aperture_type, aperture_params);	6,312✔
521	if (ap_ptr == nullptr)	6,312✔
522	{
523	std::stringstream ss;	×
524	ss << "Unable to make aperture pointer -- "
525	<< "\nChar: " << ShapeIndex
526	<< "\nType: " << aperture_type	×
527	<< "\nParams: [";	×
528	for (auto cit = aperture_params.cbegin();	×
529	cit != aperture_params.cend();	×
530	++cit)	×
531	{
532	ss << *cit << ", ";	×
533	}
534	ss << "]" << std::endl;	×
535	throw std::runtime_error(ss.str());	×
536	}	×
537	el->set_aperture(ap_ptr);	6,312✔
538
539	// Make surface
540	SurfaceType surface_type = char_to_surface(SurfaceIndex);	6,312✔
541	if (surface_type == SurfaceType::SURFACE_UNKNOWN)	6,312✔
542	{
543	std::stringstream ss;	×
544	ss << "Unknown surface type " << surface_type;	×
545	throw std::invalid_argument(ss.str());	×
546	}	×
547	surface_ptr surf_ptr = make_surface_from_type(surface_type, surface_params);	6,312✔
548	el->set_surface(surf_ptr);	6,312✔
549	if (surface_type == SurfaceType::CYLINDER)	6,312✔
550	{
551	ap_ptr = el->get_aperture();	1✔
552	auto rect = std::dynamic_pointer_cast<Rectangle>(ap_ptr);	1✔
553	auto cyl = std::dynamic_pointer_cast<Cylinder>(surf_ptr);	1✔
554	if (rect == nullptr \|\| cyl == nullptr)	1✔
555	{
556	throw std::invalid_argument("This should not happen!");	×
557	}
558	rect->x_length = 2.0 * cyl->radius;	1✔
559	}	1✔
560
561	// Set element position and orientation
562	el->set_reference_frame_geometry(Vector3d(xyz),	12,624✔
563	Vector3d(aim),	12,624✔
564	zrot);
565
566	// Set optical properties
567	OpticalProperties optics_front = optics_map[optics_name][0];	6,312✔
568	OpticalProperties optics_back = optics_map[optics_name][1];	6,312✔
569	el->set_front_optical_properties(optics_front);	6,312✔
570	el->set_back_optical_properties(optics_back);	6,312✔
571
572	// Set optical interaction type
573	el->get_front_optical_properties()->my_type = interaction;	6,312✔
574	el->get_back_optical_properties()->my_type = interaction;	6,312✔
575
576	return true;	6,312✔
577	}	6,312✔
578
579	bool process_stages(	3✔
580	FILE *fp,
581	SimulationData &sd,
582	std::map<std::string, std::array<OpticalProperties, 2>> &optics_map)
583	{
584	char buf[1024];
585
586	// Loop through stages
587	int count_stage = 0;	3✔
588	read_line(buf, 1023, fp);	3✔
589	sscanf(buf, "STAGE LIST COUNT\t%d", &count_stage);	3✔
590
591	for (int i_stage = 0; i_stage < count_stage; i_stage++)	9✔
592	{
593	int virt = 0, multi = 1, count_element = 0, tr = 0;	6✔
594	double X, Y, Z, AX, AY, AZ, ZRot;
595
596	read_line(buf, 1023, fp);	6✔
597	sscanf(buf, "STAGE\tXYZ\t%lg\t%lg\t%lg\tAIM\t%lg\t%lg\t%lg\tZROT\t%lg\tVIRTUAL\t%d\tMULTIHIT\t%d\tELEMENTS\t%d\tTRACETHROUGH\t%d",	6✔
598	&X, &Y, &Z,
599	&AX, &AY, &AZ,
600	&ZRot,
601	&virt,
602	&multi,
603	&count_element,
604	&tr);
605
606	read_line(buf, 1023, fp); // read name	6✔
607
608	// Make stage
609	stage_ptr stage = make_stage(i_stage);	6✔
610	stage->set_origin(X, Y, Z);	6✔
611	stage->set_aim_vector(AX, AY, AZ);	6✔
612	stage->set_zrot(ZRot);	6✔
613	stage->compute_coordinate_rotations();	6✔
614
615	// Loop through elements
616	for (int i_element = 0; i_element < count_element; i_element++)	6,318✔
617	{
618	element_ptr el;	6,312✔
619	read_element(fp, optics_map, el);	6,312✔
620	el->set_name(std::to_string(i_element));	6,312✔
621	// TODO make virtual if stage is virtual?
622
623	if (!Element::is_success(stage->add_element(el)))	6,312✔
624	{
625	std::cout << "Failed to add element to stage" << std::endl;	×
626	}
627	}	6,312✔
628
629	if (!Element::is_success(sd.add_stage(stage)))	6✔
630	{
631	std::cout << "Failed to add stage to SimulationData" << std::endl;	×
632	}
633	}	6✔
634
635	return true;	3✔
636	}
637
638	bool process_sim_par(FILE *fp, SimulationData &sd)	3✔
639	{
640	char buf[1024];
641
642	// Check if end of file
643	if (feof(fp))	3✔
644	return false;	×
645
646	// Check for simulation parameters
647	read_line(buf, 1023, fp);	3✔
648	if (strncmp(buf, "TRACE", 5) != 0)	3✔
649	return false;	2✔
650
651	// Get simulation parameters
652	int n_rays, n_rays_sun, n_cpu, seed, ss, err, pf;
653	int n = sscanf(buf, "TRACE\tNRAY\t%d\tNSUN\t%d\tCPU\t%d\tSEED\t%d\tSUNSHAPE\t%d\tERRORS\t%d\tPTFOCUS\t%d",	1✔
654	&n_rays, &n_rays_sun, &n_cpu, &seed, &ss, &err, &pf);
655
656	// Assign simulation parameters
657	SimulationParameters &par = sd.get_simulation_parameters();	1✔
658	par.number_of_rays = n_rays;	1✔
659	par.max_number_of_rays = n_rays_sun;	1✔
660	par.seed = seed;	1✔
661	par.include_sun_shape_errors = ss;	1✔
662	par.include_optical_errors = err;	1✔
663
664	// TODO Assign number CPUs, point focus?
665	return true;	1✔
666	}
667
668	bool load_stinput_file(SimulationData &sd, std::string filename)	3✔
669	{
670	// TODO: Reset simulation data?
671
672	// Read in file
673	FILE *fp = fopen(filename.data(), "r");	3✔
674	if (!fp)	3✔
675	{
676	printf("failed to open system input file: %s\n",	×
677	filename.data());
678	return false;	×
679	}
680
681	// Buffer to store read line
682	char buf[1024];
683
684	// Get version info (if first line starts with '#')
685	int vmaj = 0, vmin = 0, vmic = 0;	3✔
686	char c = fgetc(fp);	3✔
687	if (c == '#')	3✔
688	{
689	read_line(buf, 1023, fp);	3✔
690	sscanf(buf, " SOLTRACE VERSION %d.%d.%d INPUT FILE",	3✔
691	&vmaj, &vmin, &vmic);
692
693	// unsigned int file_version = vmaj10000 + vmin100 + vmic;
694
695	printf("loading input file version %d.%d.%d\n",	3✔
696	vmaj, vmin, vmic);
697	}
698	else
699	{
700	ungetc(c, fp);	×
701	}
702
703	// Read in Sun
704	process_sun(fp, sd);	3✔
705
706	// Read in Optics
707	std::map<std::string, std::array<OpticalProperties, 2>> optics_map;	3✔
708	process_optics(fp, optics_map);	3✔
709
710	// Read in Stages
711	process_stages(fp, sd, optics_map);	3✔
712
713	// Read in simulation parameters (if any)
714	process_sim_par(fp, sd);	3✔
715
716	return true;	3✔
717	}	3✔
718
719	} // namespace SolTrace::Data

NREL / SolTrace / 18984697163

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