commit 43e47fdf51a03eb559a82c921c67886502795f0c
parent 214403f413875775f5b346d568ba96a34b3ea9fd
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Fri, 8 Dec 2017 16:00:25 +0100
Add gaussian shunshape ans associated test.
Diffstat:
8 files changed, 355 insertions(+), 0 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -183,6 +183,7 @@ if(NOT NO_TEST)
new_test(test_ssol_solver9)
new_test(test_ssol_solver10)
new_test(test_ssol_solver11)
+ new_test(test_ssol_solver12)
new_test(test_ssol_sun)
build_test(test_ssol_draw)
diff --git a/src/ssol.h b/src/ssol.h
@@ -996,6 +996,12 @@ ssol_sun_create_pillbox
(struct ssol_device* dev,
struct ssol_sun** sun);
+/* The sun disk intensity has a gaussian shape */
+SSOL_API res_T
+ssol_sun_create_gaussian
+ (struct ssol_device* dev,
+ struct ssol_sun** sun);
+
/* The sun disk intensity is controlled by a circumsolar ratio. From the paper
* "Sunshape distributions for terrestrial solar simulations". D. Buie, A.G.
* Monger, C.J. Dey */
@@ -1045,6 +1051,11 @@ ssol_sun_pillbox_set_half_angle
const double half_angle); /* In ]0, PI/2], in radian */
SSOL_API res_T
+ssol_sun_gaussian_set_std_dev
+ (struct ssol_sun* sun,
+ const double std_dev); /* In ]0, +Inf[, in radian */
+
+SSOL_API res_T
ssol_sun_set_buie_param
(struct ssol_sun* sun,
const double param); /* In ]0, 1[ */
diff --git a/src/ssol_ranst_sun_dir.c b/src/ssol_ranst_sun_dir.c
@@ -44,6 +44,11 @@ struct ran_pillbox_state {
double basis[9];
};
+struct ran_gaussian_state {
+ double std_dev;
+ double basis[9];
+};
+
struct ran_dirac_state {
double dir[3];
};
@@ -56,6 +61,7 @@ struct ranst_sun_dir {
union {
struct ran_buie_state buie;
struct ran_pillbox_state pillbox;
+ struct ran_gaussian_state gaussian;
struct ran_dirac_state dirac;
} state;
@@ -260,6 +266,30 @@ ran_pillbox_get
}
/*******************************************************************************
+* Gaussian random variate
+******************************************************************************/
+static double*
+ran_gaussian_get
+ (const struct ranst_sun_dir* ran,
+ struct ssp_rng* rng,
+ double dir[3])
+{
+ double pt[3];
+ double phi, theta, sin_theta;
+
+ ASSERT(ran && 0 <= ran->state.gaussian.std_dev);
+ theta = ssp_ran_gaussian(rng, 0, ran->state.gaussian.std_dev);
+ sin_theta = sin(theta);
+ phi = ssp_rng_uniform_double(rng, 0, 2 * PI);
+ pt[0] = cos(phi) * sin_theta;
+ pt[1] = sin(phi) * sin_theta;
+ pt[2] = cos(theta);
+ d33_muld3(dir, ran->state.pillbox.basis, pt);
+ d3_normalize(dir, dir);
+ return dir;
+}
+
+/*******************************************************************************
* Dirac distribution
******************************************************************************/
static double*
@@ -358,6 +388,19 @@ ranst_sun_dir_pillbox_setup
}
res_T
+ranst_sun_dir_gaussian_setup
+ (struct ranst_sun_dir* ran,
+ const double std_dev,
+ const double dir[3])
+{
+ if(!ran || !dir || std_dev <= 0) return RES_BAD_ARG;
+ ran->get = ran_gaussian_get;
+ ran->state.gaussian.std_dev = std_dev;
+ d33_basis(ran->state.pillbox.basis, dir);
+ return RES_OK;
+}
+
+res_T
ranst_sun_dir_dirac_setup
(struct ranst_sun_dir* ran,
const double dir[3])
diff --git a/src/ssol_ranst_sun_dir.h b/src/ssol_ranst_sun_dir.h
@@ -55,6 +55,12 @@ ranst_sun_dir_pillbox_setup
const double dir[3]);
extern LOCAL_SYM res_T
+ranst_sun_dir_gaussian_setup
+ (struct ranst_sun_dir* ran,
+ const double std_dev,
+ const double dir[3]);
+
+extern LOCAL_SYM res_T
ranst_sun_dir_dirac_setup
(struct ranst_sun_dir* ran,
const double dir[3]);
diff --git a/src/ssol_sun.c b/src/ssol_sun.c
@@ -93,6 +93,12 @@ ssol_sun_create_pillbox(struct ssol_device* dev, struct ssol_sun** out_sun)
}
res_T
+ssol_sun_create_gaussian(struct ssol_device* dev, struct ssol_sun** out_sun)
+{
+ return sun_create(dev, out_sun, SUN_GAUSSIAN);
+}
+
+res_T
ssol_sun_create_buie
(struct ssol_device* dev, struct ssol_sun** out_sun)
{
@@ -179,6 +185,15 @@ ssol_sun_pillbox_set_half_angle(struct ssol_sun* sun, const double half_angle)
}
res_T
+ssol_sun_gaussian_set_std_dev(struct ssol_sun* sun, const double std_dev)
+{
+ if(!sun || std_dev <= 0 || sun->type != SUN_GAUSSIAN)
+ return RES_BAD_ARG;
+ sun->data.gaussian.std_dev = std_dev;
+ return RES_OK;
+}
+
+res_T
ssol_sun_set_buie_param
(struct ssol_sun* sun,
const double ratio)
@@ -213,6 +228,10 @@ sun_create_direction_distribution
res = ranst_sun_dir_pillbox_setup
(ran_dir, sun->data.pillbox.half_angle, sun->direction);
break;
+ case SUN_GAUSSIAN:
+ res = ranst_sun_dir_gaussian_setup
+ (ran_dir, sun->data.gaussian.std_dev, sun->direction);
+ break;
case SUN_BUIE:
res = ranst_sun_dir_buie_setup
(ran_dir, sun->data.csr.ratio, sun->direction);
diff --git a/src/ssol_sun_c.h b/src/ssol_sun_c.h
@@ -26,6 +26,7 @@ struct ranst_sun_wl;
enum sun_type {
SUN_DIRECTIONAL,
SUN_PILLBOX,
+ SUN_GAUSSIAN,
SUN_BUIE,
SUN_TYPES_COUNT__
};
@@ -34,6 +35,10 @@ struct pillbox {
double half_angle;
};
+struct gaussian {
+ double std_dev;
+};
+
struct buie {
double ratio;
};
@@ -46,6 +51,7 @@ struct ssol_sun {
enum sun_type type;
union {
struct pillbox pillbox;
+ struct gaussian gaussian;
struct buie csr;
} data;
diff --git a/src/test_ssol_solver12.c b/src/test_ssol_solver12.c
@@ -0,0 +1,202 @@
+/* Copyright (C) CNRS 2016-2017
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "ssol.h"
+#include "test_ssol_utils.h"
+#define REFLECTIVITY 0
+#include "test_ssol_materials.h"
+
+#include <rsys/math.h>
+
+#define X_SZ 10
+#define Y_SZ 10
+#define PLANE_NAME SQUARE
+#define HALF_X (X_SZ / 2)
+#define HALF_Y (Y_SZ / 2)
+STATIC_ASSERT((HALF_X * 2 == X_SZ), ONLY_ENVEN_VALUES_FOR_SQUARE_X);
+STATIC_ASSERT((HALF_Y * 2 == Y_SZ), ONLY_ENVEN_VALUES_FOR_SQUARE_Y);
+#include "test_ssol_rect_geometry.h"
+
+#define POLYGON_NAME POLY
+#define HALF_X (X_SZ / 2)
+#define HALF_Y (Y_SZ / 2)
+STATIC_ASSERT((HALF_X * 2 == X_SZ), ONLY_ENVEN_VALUES_FOR_X_SZ);
+STATIC_ASSERT((HALF_Y * 2 == Y_SZ), ONLY_ENVEN_VALUES_FOR_Y_SZ);
+#include "test_ssol_rect2D_geometry.h"
+
+#include <rsys/double33.h>
+
+#include <star/s3d.h>
+#include <star/ssp.h>
+
+static void
+get_wlen(const size_t i, double* wlen, double* data, void* ctx)
+{
+ double wavelengths[3] = { 1, 2, 3 };
+ double intensities[3] = { 1, 0.8, 1 };
+ CHECK(i < 3, 1);
+ (void) ctx;
+ *wlen = wavelengths[i];
+ *data = intensities[i];
+}
+
+int
+main(int argc, char** argv)
+{
+ struct mem_allocator allocator;
+ struct ssol_device* dev;
+ struct ssp_rng* rng;
+ struct ssol_scene* scene;
+ struct ssol_shape* square;
+ struct ssol_shape* parabol;
+ struct ssol_vertex_data attribs[1] = { SSOL_VERTEX_DATA_NULL__ };
+ struct ssol_material* m_mtl;
+ struct ssol_matte_shader shader = SSOL_MATTE_SHADER_NULL;
+ struct ssol_object* m_object;
+ struct ssol_object* q_object;
+ struct ssol_carving carving = SSOL_CARVING_NULL;
+ struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT;
+ struct ssol_punched_surface punched = SSOL_PUNCHED_SURFACE_NULL;
+ struct ssol_instance* geom1;
+ struct ssol_instance* geom2;
+ struct ssol_sun* sun;
+ struct ssol_spectrum* spectrum;
+ struct ssol_estimator* estimator;
+ struct ssol_mc_global mc_global1;
+ struct ssol_mc_global mc_global2;
+ struct ssol_mc_receiver mc_rcv1;
+ struct ssol_mc_receiver mc_rcv2;
+ double dir[3];
+ double transform[12]; /* 3x4 column major matrix */
+ size_t count;
+
+ (void) argc, (void) argv;
+ mem_init_proxy_allocator(&allocator, &mem_default_allocator);
+
+ d33_set_identity(transform);
+
+ CHECK(ssol_device_create
+ (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev), RES_OK);
+
+#define DNI 1000
+ CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK);
+ CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK);
+ CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK);
+ CHECK(ssol_sun_create_gaussian(dev, &sun), RES_OK);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, 0.01), RES_OK);
+ CHECK(ssol_sun_set_direction(sun, d3(dir, 0, 0, -1)), RES_OK);
+ CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK);
+ CHECK(ssol_sun_set_dni(sun, DNI), RES_OK);
+ CHECK(ssol_scene_create(dev, &scene), RES_OK);
+ CHECK(ssol_scene_attach_sun(scene, sun), RES_OK);
+
+ /* Create scene content */
+
+ CHECK(ssol_shape_create_mesh(dev, &square), RES_OK);
+ attribs[0].usage = SSOL_POSITION;
+ attribs[0].get = get_position;
+ CHECK(ssol_mesh_setup(square, SQUARE_NTRIS__, get_ids,
+ SQUARE_NVERTS__, attribs, 1, (void*) &SQUARE_DESC__), RES_OK);
+
+ CHECK(ssol_material_create_matte(dev, &m_mtl), RES_OK);
+ shader.normal = get_shader_normal;
+ shader.reflectivity = get_shader_reflectivity_2;
+ CHECK(ssol_matte_setup(m_mtl, &shader), RES_OK);
+
+ CHECK(ssol_object_create(dev, &m_object), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(m_object, square, m_mtl, m_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(m_object, &geom1), RES_OK);
+ CHECK(ssol_instance_set_receiver(geom1, SSOL_FRONT, 0), RES_OK);
+ d3_splat(transform + 9, 0);
+ transform[9] = -10;
+ CHECK(ssol_instance_set_transform(geom1, transform), RES_OK);
+ CHECK(ssol_scene_attach_instance(scene, geom1), RES_OK);
+
+#define N1__ 10000
+#define GET_MC_RCV ssol_estimator_get_mc_receiver
+ CHECK(ssol_solve(scene, rng, N1__, 0, NULL, &estimator), RES_OK);
+ CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
+ CHECK(count, N1__);
+ CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
+ CHECK(count, 0);
+#define DNI_S (DNI * X_SZ * Y_SZ)
+ CHECK(ssol_estimator_get_mc_global(estimator, &mc_global1), RES_OK);
+ CHECK(GET_MC_RCV(estimator, geom1, SSOL_FRONT, &mc_rcv1), RES_OK);
+
+ PRINT_GLOBAL(mc_global1);
+ PRINT_RCV(mc_rcv1);
+ CHECK(mc_global1.cos_factor.E, 1);
+ CHECK(mc_global1.cos_factor.SE, 0);
+ CHECK(mc_global1.absorbed_by_receivers.E, DNI_S);
+ CHECK(mc_global1.absorbed_by_receivers.SE, 0);
+
+ CHECK(ssol_shape_create_punched_surface(dev, ¶bol), RES_OK);
+ carving.get = get_polygon_vertices;
+ carving.operation = SSOL_AND;
+ carving.nb_vertices = POLY_NVERTS__;
+ carving.context = &POLY_EDGES__;
+ quadric.type = SSOL_QUADRIC_PARABOL;
+ quadric.data.parabol.focal = 10;
+ punched.nb_carvings = 1;
+ punched.quadric = &quadric;
+ punched.carvings = &carving;
+ CHECK(ssol_punched_surface_setup(parabol, &punched), RES_OK);
+
+ CHECK(ssol_object_create(dev, &q_object), RES_OK);
+ CHECK(ssol_object_add_shaded_shape(q_object, parabol, m_mtl, m_mtl), RES_OK);
+ CHECK(ssol_object_instantiate(q_object, &geom2), RES_OK);
+ CHECK(ssol_instance_set_receiver(geom2, SSOL_FRONT, 0), RES_OK);
+ d3_splat(transform + 9, 0);
+ transform[9] = +10;
+ CHECK(ssol_instance_set_transform(geom2, transform), RES_OK);
+ CHECK(ssol_scene_attach_instance(scene, geom2), RES_OK);
+
+ CHECK(ssol_scene_detach_instance(scene, geom1), RES_OK);
+ CHECK(ssol_estimator_ref_put(estimator), RES_OK);
+
+#define N2__ 100000
+ CHECK(ssol_solve(scene, rng, N2__, 0, NULL, &estimator), RES_OK);
+ CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
+ CHECK(count, N2__);
+ CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
+ CHECK(count, 0);
+ CHECK(ssol_estimator_get_mc_global(estimator, &mc_global2), RES_OK);
+ CHECK(GET_MC_RCV(estimator, geom2, SSOL_FRONT, &mc_rcv2), RES_OK);
+
+ PRINT_GLOBAL(mc_global2);
+ PRINT_RCV(mc_rcv2);
+ CHECK(eq_eps(mc_global2.absorbed_by_receivers.E, DNI_S, 3 * mc_global2.absorbed_by_receivers.SE), 1);
+
+ /* Free data */
+ CHECK(ssol_instance_ref_put(geom1), RES_OK);
+ CHECK(ssol_instance_ref_put(geom2), RES_OK);
+ CHECK(ssol_object_ref_put(m_object), RES_OK);
+ CHECK(ssol_object_ref_put(q_object), RES_OK);
+ CHECK(ssol_shape_ref_put(square), RES_OK);
+ CHECK(ssol_shape_ref_put(parabol), RES_OK);
+ CHECK(ssol_material_ref_put(m_mtl), RES_OK);
+ CHECK(ssol_estimator_ref_put(estimator), RES_OK);
+ CHECK(ssol_device_ref_put(dev), RES_OK);
+ CHECK(ssol_scene_ref_put(scene), RES_OK);
+ CHECK(ssp_rng_ref_put(rng), RES_OK);
+ CHECK(ssol_spectrum_ref_put(spectrum), RES_OK);
+ CHECK(ssol_sun_ref_put(sun), RES_OK);
+
+ check_memory_allocator(&allocator);
+ mem_shutdown_proxy_allocator(&allocator);
+ CHECK(mem_allocated_size(), 0);
+
+ return 0;
+}
diff --git a/src/test_ssol_sun.c b/src/test_ssol_sun.c
@@ -81,6 +81,10 @@ main(int argc, char** argv)
CHECK(ssol_sun_pillbox_set_half_angle(sun, 999), RES_BAD_ARG);
CHECK(ssol_sun_pillbox_set_half_angle(sun, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(NULL, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, -0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, 0.1), RES_BAD_ARG);
+
CHECK(ssol_sun_set_buie_param(NULL, 0.1), RES_BAD_ARG);
CHECK(ssol_sun_set_buie_param(sun, -0.1), RES_BAD_ARG);
CHECK(ssol_sun_set_buie_param(sun, 999), RES_BAD_ARG);
@@ -131,6 +135,65 @@ main(int argc, char** argv)
CHECK(ssol_sun_pillbox_set_half_angle(sun, 0.1), RES_OK);
CHECK(ssol_sun_pillbox_set_half_angle(sun, 0.1), RES_OK);
+ CHECK(ssol_sun_gaussian_set_std_dev(NULL, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, -0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, 0.1), RES_BAD_ARG);
+
+ CHECK(ssol_sun_set_buie_param(NULL, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_set_buie_param(sun, -0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_set_buie_param(sun, 999), RES_BAD_ARG);
+ CHECK(ssol_sun_set_buie_param(sun, 0.1), RES_BAD_ARG);
+
+ CHECK(ssol_sun_ref_put(sun), RES_OK);
+
+ CHECK(ssol_sun_create_gaussian(NULL, &sun), RES_BAD_ARG);
+ CHECK(ssol_sun_create_gaussian(dev, NULL), RES_BAD_ARG);
+ CHECK(ssol_sun_create_gaussian(dev, &sun), RES_OK);
+
+ CHECK(ssol_sun_ref_get(NULL), RES_BAD_ARG);
+ CHECK(ssol_sun_ref_get(sun), RES_OK);
+
+ CHECK(ssol_sun_ref_put(NULL), RES_BAD_ARG);
+ CHECK(ssol_sun_ref_put(sun), RES_OK);
+
+ CHECK(ssol_sun_set_spectrum(NULL, spectrum), RES_BAD_ARG);
+ CHECK(ssol_sun_set_spectrum(sun, NULL), RES_BAD_ARG);
+ CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK);
+ CHECK(ssol_sun_set_spectrum(sun, spectrum2), RES_OK);
+ CHECK(ssol_sun_set_spectrum(sun, spectrum2), RES_OK);
+
+ CHECK(ssol_sun_set_direction(NULL, dir), RES_BAD_ARG);
+ CHECK(ssol_sun_set_direction(sun, NULL), RES_BAD_ARG);
+ CHECK(ssol_sun_set_direction(sun, dir0), RES_BAD_ARG);
+ CHECK(ssol_sun_set_direction(sun, dir), RES_OK);
+ CHECK(ssol_sun_set_direction(sun, dir), RES_OK);
+
+ CHECK(ssol_sun_get_direction(NULL, tmp), RES_BAD_ARG);
+ CHECK(ssol_sun_get_direction(sun, NULL), RES_BAD_ARG);
+ CHECK(ssol_sun_get_direction(sun, tmp), RES_OK);
+ CHECK(d3_eq(dir, tmp), 1);
+
+ CHECK(ssol_sun_set_dni(NULL, 1000), RES_BAD_ARG);
+ CHECK(ssol_sun_set_dni(sun, 0), RES_BAD_ARG);
+ CHECK(ssol_sun_set_dni(sun, 1000), RES_OK);
+ CHECK(ssol_sun_set_dni(sun, 1000), RES_OK);
+
+ CHECK(ssol_sun_get_dni(NULL, &dni), RES_BAD_ARG);
+ CHECK(ssol_sun_get_dni(sun, NULL), RES_BAD_ARG);
+ CHECK(ssol_sun_get_dni(sun, &dni), RES_OK);
+ CHECK(dni, 1000);
+
+ CHECK(ssol_sun_pillbox_set_half_angle(NULL, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_pillbox_set_half_angle(sun, -0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_pillbox_set_half_angle(sun, 999), RES_BAD_ARG);
+ CHECK(ssol_sun_pillbox_set_half_angle(sun, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_pillbox_set_half_angle(sun, 0.1), RES_BAD_ARG);
+
+ CHECK(ssol_sun_gaussian_set_std_dev(NULL, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, -0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, 0.1), RES_OK);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, 0.1), RES_OK);
+
CHECK(ssol_sun_set_buie_param(NULL, 0.1), RES_BAD_ARG);
CHECK(ssol_sun_set_buie_param(sun, -0.1), RES_BAD_ARG);
CHECK(ssol_sun_set_buie_param(sun, 999), RES_BAD_ARG);
@@ -179,6 +242,10 @@ main(int argc, char** argv)
CHECK(ssol_sun_pillbox_set_half_angle(sun, 999), RES_BAD_ARG);
CHECK(ssol_sun_pillbox_set_half_angle(sun, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(NULL, 0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, -0.1), RES_BAD_ARG);
+ CHECK(ssol_sun_gaussian_set_std_dev(sun, 0.1), RES_BAD_ARG);
+
CHECK(ssol_sun_set_buie_param(NULL, 0.1), RES_BAD_ARG);
CHECK(ssol_sun_set_buie_param(sun, -0.1), RES_BAD_ARG);
CHECK(ssol_sun_set_buie_param(sun, 999), RES_BAD_ARG);
