commit 6901f6b770584b2175eae2e8e446c3f335693176
parent a8267b9da8eec0203d1840860a6ef9c9c51b2245
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 8 Sep 2017 11:38:35 +0200
Merge branch 'release_0.3'
Diffstat:
25 files changed, 699 insertions(+), 647 deletions(-)
diff --git a/README.md b/README.md
@@ -26,6 +26,17 @@ variable the install directories of its dependencies.
## Release notes
+### Version 0.3
+
+- Full rewrite of the estimated values. The global results report the cosine
+ factor, and the overall flux that is: absorbed by the receivers, atmosphere
+ or others entities; occluded before it reaches a primary entity; missed
+ because it does not reaches any surface. The per receiver results include the
+ incoming/absorbed flux in 3 situations: all phenomenons are taken into
+ account; the atmosphere is disabled; the material propagate the whole
+ incoming flux, i.e. they absorbed nothing.
+- Update the `ssol_solve` API. Streamed binary outputs are removed.
+
### Version 0.2.2
- Fix the estimation of the cosine factor for the sampled instances: it was
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -28,7 +28,7 @@ find_package(RSys 0.4 REQUIRED)
find_package(Star3D 0.4.1 REQUIRED)
find_package(Star3DUT 0.2 REQUIRED)
find_package(StarCPR 0.1 REQUIRED)
-find_package(StarSF 0.1 REQUIRED)
+find_package(StarSF 0.2 REQUIRED)
find_package(StarSP 0.4 REQUIRED)
find_package(OpenMP 1.2 REQUIRED)
@@ -50,8 +50,8 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D Star3DUT StarCPR StarSF Sta
# Configure and define targets
################################################################################
set(VERSION_MAJOR 0)
-set(VERSION_MINOR 2)
-set(VERSION_PATCH 2)
+set(VERSION_MINOR 3)
+set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
set(SSOL_FILES_SRC
@@ -133,7 +133,7 @@ if(NOT NO_TEST)
function(build_test _name)
add_executable(${_name} ${SSOL_SOURCE_DIR}/${_name}.c)
target_link_libraries(${_name}
- solstice-solver solstice-test RSys Star3D StarSP)
+ solstice-solver RSys Star3D StarSP)
endfunction()
function(register_test _name)
@@ -145,13 +145,7 @@ if(NOT NO_TEST)
build_test(${_name})
register_test(${_name} ${_name})
endfunction()
-
- add_library(solstice-test STATIC
- ${SSOL_SOURCE_DIR}/test_ssol_geometries.h
- ${SSOL_SOURCE_DIR}/test_ssol_materials.h
- ${SSOL_SOURCE_DIR}/test_ssol_utils.h
- ${SSOL_SOURCE_DIR}/test_ssol_utils.c)
-
+
new_test(test_ssol_atmosphere)
new_test(test_ssol_by_receiver_integration)
new_test(test_ssol_camera)
diff --git a/src/ssol.h b/src/ssol.h
@@ -28,8 +28,8 @@
#endif
/* Helper macro that asserts if the invocation of the Solstice function `Func'
- * returns an error. One should use this macro on Solstice function calls for which
- * no explicit error checking is performed */
+ * returns an error. One should use this macro on Solstice function calls for
+ * which no explicit error checking is performed */
#ifndef NDEBUG
#define SSOL(Func) ASSERT(ssol_ ## Func == RES_OK)
#else
@@ -178,7 +178,7 @@ static const struct ssol_quadric_parabol SSOL_QUADRIC_PARABOL_NULL =
struct ssol_quadric_hyperbol {
/* Define (x^2 + y^2) / a^2 - (z - 1/2)^2 / b^2 + 1 = 0; z > 0
- * with a^2 = f - f^2; b = f -1/2; f = real_focal / (img_focal + real_focal) */
+ * with a^2 = f - f^2; b = f -1/2; f = real_focal/(img_focal + real_focal) */
double img_focal, real_focal;
};
#define SSOL_QUADRIC_HYPERBOL_NULL__ { -1.0 , -1.0 }
@@ -322,26 +322,6 @@ struct ssol_thin_dielectric_shader {
static const struct ssol_thin_dielectric_shader
SSOL_THIN_DIELECTRIC_SHADER_NULL = SSOL_THIN_DIELECTRIC_SHADER_NULL__;
-/* The type of data produced on receiver hits as ssol_solve() write them on its
- * FILE* argument */
-struct ssol_receiver_data {
- uint64_t realization_id;
- uint32_t segment_id;
-
- /* Its absolute value is the identifier of an SSOL instance. A negative
- * value means for back faces primitive */
- int32_t receiver_id;
-
- float wavelength;
- float pos[3];
- float in_dir[3];
- float normal[3];
- double weight;
- float uv[2];
-
- /* TODO Add the geometry and primitive identifier */
-};
-
struct ssol_instantiated_shaded_shape {
struct ssol_shape* shape;
struct ssol_material* mtl_front;
@@ -388,11 +368,11 @@ static const struct ssol_mc_result SSOL_MC_RESULT_NULL = SSOL_MC_RESULT_NULL__;
struct ssol_mc_global {
struct ssol_mc_result cos_factor; /* [0 1] */
- struct ssol_mc_result absorbed; /* In W */
+ struct ssol_mc_result absorbed_by_receivers; /* In W */
struct ssol_mc_result shadowed; /* In W */
struct ssol_mc_result missing; /* In W */
- struct ssol_mc_result atmosphere; /* In W */
- struct ssol_mc_result reflectivity; /* In W */
+ struct ssol_mc_result absorbed_by_atmosphere; /* In W */
+ struct ssol_mc_result other_absorbed; /* In W */
};
#define SSOL_MC_GLOBAL_NULL__ { \
SSOL_MC_RESULT_NULL__, \
@@ -405,10 +385,16 @@ struct ssol_mc_global {
static const struct ssol_mc_global SSOL_MC_GLOBAL_NULL = SSOL_MC_GLOBAL_NULL__;
struct ssol_mc_receiver {
- struct ssol_mc_result integrated_irradiance; /* In W */
- struct ssol_mc_result integrated_absorbed_irradiance; /* In W */
- struct ssol_mc_result absorptivity_loss; /* In W */
- struct ssol_mc_result reflectivity_loss; /* In W */
+ struct ssol_mc_result incoming_flux; /* In W */
+ struct ssol_mc_result incoming_if_no_atm_loss; /* In W */
+ struct ssol_mc_result incoming_if_no_field_loss; /* In W */
+ struct ssol_mc_result incoming_lost_in_field; /* In W */
+ struct ssol_mc_result incoming_lost_in_atmosphere; /* In W */
+ struct ssol_mc_result absorbed_flux; /* In W */
+ struct ssol_mc_result absorbed_if_no_atm_loss; /* In W */
+ struct ssol_mc_result absorbed_if_no_field_loss; /* In W */
+ struct ssol_mc_result absorbed_lost_in_field; /* In W */
+ struct ssol_mc_result absorbed_lost_in_atmosphere; /* In W */
/* Internal data */
size_t N__;
@@ -420,6 +406,12 @@ struct ssol_mc_receiver {
SSOL_MC_RESULT_NULL__, \
SSOL_MC_RESULT_NULL__, \
SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
0, NULL, NULL \
}
static const struct ssol_mc_receiver SSOL_MC_RECEIVER_NULL =
@@ -430,6 +422,12 @@ static const struct ssol_mc_receiver SSOL_MC_RECEIVER_NULL =
{ -1, -1, -1 }, \
{ -1, -1, -1 }, \
{ -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
+ { -1, -1, -1 }, \
0, NULL, NULL \
}
@@ -449,15 +447,27 @@ struct ssol_mc_sampled {
};
struct ssol_mc_primitive {
- struct ssol_mc_result integrated_irradiance; /* In W */
- struct ssol_mc_result integrated_absorbed_irradiance; /* In W */
- struct ssol_mc_result absorptivity_loss; /* In W */
- struct ssol_mc_result reflectivity_loss; /* In W */
+ struct ssol_mc_result incoming_flux; /* In W */
+ struct ssol_mc_result incoming_if_no_atm_loss; /* In W */
+ struct ssol_mc_result incoming_if_no_field_loss; /* In W */
+ struct ssol_mc_result incoming_lost_in_field; /* In W */
+ struct ssol_mc_result incoming_lost_in_atmosphere; /* In W */
+ struct ssol_mc_result absorbed_flux; /* In W */
+ struct ssol_mc_result absorbed_if_no_atm_loss; /* In W */
+ struct ssol_mc_result absorbed_if_no_field_loss; /* In W */
+ struct ssol_mc_result absorbed_lost_in_field; /* In W */
+ struct ssol_mc_result absorbed_lost_in_atmosphere; /* In W */
};
#define SSOL_MC_PRIMITIVE_NULL__ { \
SSOL_MC_RESULT_NULL__, \
SSOL_MC_RESULT_NULL__, \
SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
+ SSOL_MC_RESULT_NULL__, \
SSOL_MC_RESULT_NULL__ \
}
static const struct ssol_mc_primitive SSOL_MC_PRIMITIVE_NULL =
@@ -712,7 +722,7 @@ ssol_shape_get_triangle_indices
const unsigned itri,
unsigned ids[3]);
-/* Define a punched surface in local space, i.e. no translation & no orientation */
+/* Define a punched surface in local space, i.e. no transformation */
SSOL_API res_T
ssol_punched_surface_setup
(struct ssol_shape* shape,
@@ -1061,8 +1071,8 @@ ssol_atmosphere_set_absorption
struct ssol_data* absorption);
/*******************************************************************************
-* Estimator API - Describe the state of a simulation.
-******************************************************************************/
+ * Estimator API - Describe the state of a simulation.
+ ******************************************************************************/
SSOL_API res_T
ssol_estimator_ref_get
(struct ssol_estimator* estimator);
@@ -1073,7 +1083,7 @@ ssol_estimator_ref_put
SSOL_API res_T
ssol_estimator_get_mc_global
- (const struct ssol_estimator* estimator,
+ (struct ssol_estimator* estimator,
struct ssol_mc_global* mc_global);
SSOL_API res_T
@@ -1172,7 +1182,6 @@ ssol_solve
struct ssp_rng* rng,
const size_t realisations_count,
const struct ssol_path_tracker* tracker, /* NULL<=>Do not record the paths */
- FILE* output, /* May be NULL <=> does not ouput ssol_receiver_data */
struct ssol_estimator** estimator);
SSOL_API res_T
diff --git a/src/ssol_c.h b/src/ssol_c.h
@@ -34,6 +34,7 @@ struct ray_data {
struct ssol_scene* scn; /* The scene into which the ray is traced */
struct s3d_primitive prim_from; /* Primitive from which the ray starts */
const struct ssol_instance* inst_from; /* Instance from which the ray starts */
+ const struct shaded_shape* sshape_from; /* Shape from which the ray starts */
enum ssol_side_flag side_from; /* Primitive side from which the ray starts */
short discard_virtual_materials; /* Define if virtual materials are not RT */
short reversed_ray; /* Define if the ray direction is reversed */
@@ -45,7 +46,7 @@ struct ray_data {
};
static const struct ray_data RAY_DATA_NULL = {
- NULL, S3D_PRIMITIVE_NULL__, NULL, SSOL_INVALID_SIDE, 0, 0, 0, {0,0,0}, FLT_MAX
+ NULL, S3D_PRIMITIVE_NULL__, NULL, NULL, SSOL_INVALID_SIDE, 0, 0, 0, {0,0,0}, FLT_MAX
};
diff --git a/src/ssol_estimator.c b/src/ssol_estimator.c
@@ -107,23 +107,26 @@ ssol_estimator_ref_put(struct ssol_estimator* estimator)
res_T
ssol_estimator_get_mc_global
- (const struct ssol_estimator* estimator,
+ (struct ssol_estimator* estimator,
struct ssol_mc_global* global)
{
if(!estimator || !global) return RES_BAD_ARG;
#define SETUP_MC_RESULT(Name) { \
const double N = (double)estimator->realisation_count; \
- const struct mc_data* data = &estimator->Name; \
- global->Name.E = data->weight / N; \
- global->Name.V = data->sqr_weight / N - global->Name.E*global->Name.E; \
- global->Name.SE = global->Name.V > 0 ? sqrt(global->Name.V / N) : 0; \
+ struct mc_data* data = &estimator->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ global->Name.E = weight / N; \
+ global->Name.V = sqr_weight / N - global->Name.E*global->Name.E; \
+ global->Name.V = global->Name.V > 0 ? global->Name.V : 0; \
+ global->Name.SE = sqrt(global->Name.V / N); \
} (void)0
SETUP_MC_RESULT(cos_factor);
- SETUP_MC_RESULT(absorbed);
+ SETUP_MC_RESULT(absorbed_by_receivers);
SETUP_MC_RESULT(shadowed);
SETUP_MC_RESULT(missing);
- SETUP_MC_RESULT(atmosphere);
- SETUP_MC_RESULT(reflectivity);
+ SETUP_MC_RESULT(absorbed_by_atmosphere);
+ SETUP_MC_RESULT(other_absorbed);
#undef SETUP_MC_RESULT
return RES_OK;
}
@@ -164,15 +167,24 @@ ssol_estimator_get_mc_sampled_x_receiver
mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back;
#define SETUP_MC_RESULT(Name) { \
const double N = (double)estimator->realisation_count; \
- const struct mc_data* data = &mc_rcv1->Name; \
- rcv->Name.E = data->weight / N; \
- rcv->Name.V = data->sqr_weight / N - rcv->Name.E*rcv->Name.E; \
- rcv->Name.SE = rcv->Name.V > 0 ? sqrt(rcv->Name.V / N) : 0; \
+ struct mc_data* data = &mc_rcv1->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ rcv->Name.E = weight / N; \
+ rcv->Name.V = sqr_weight / N - rcv->Name.E*rcv->Name.E; \
+ rcv->Name.V = rcv->Name.V > 0 ? rcv->Name.V : 0; \
+ rcv->Name.SE = sqrt(rcv->Name.V / N); \
} (void)0
- SETUP_MC_RESULT(integrated_irradiance);
- SETUP_MC_RESULT(integrated_absorbed_irradiance);
- SETUP_MC_RESULT(absorptivity_loss);
- SETUP_MC_RESULT(reflectivity_loss);
+ SETUP_MC_RESULT(incoming_flux);
+ SETUP_MC_RESULT(incoming_if_no_atm_loss);
+ SETUP_MC_RESULT(incoming_if_no_field_loss);
+ SETUP_MC_RESULT(incoming_lost_in_field);
+ SETUP_MC_RESULT(incoming_lost_in_atmosphere);
+ SETUP_MC_RESULT(absorbed_flux);
+ SETUP_MC_RESULT(absorbed_if_no_atm_loss);
+ SETUP_MC_RESULT(absorbed_if_no_field_loss);
+ SETUP_MC_RESULT(absorbed_lost_in_field);
+ SETUP_MC_RESULT(absorbed_lost_in_atmosphere);
#undef SETUP_MC_RESULT
rcv->mc__ = mc_rcv1;
rcv->N__ = mc_samp->nb_samples;
@@ -226,15 +238,18 @@ ssol_estimator_get_mc_sampled
mc = htable_sampled_find(&estimator->mc_sampled, &samp_instance);
if(!mc) return RES_BAD_ARG;
sampled->nb_samples = mc->nb_samples;
- #define SETUP_MC_RESULT(Name) { \
- const double N = (double)estimator->realisation_count; \
- const struct mc_data* data = &mc->Name; \
- sampled->Name.E = data->weight / N; \
- sampled->Name.V = data->sqr_weight/N - sampled->Name.E*sampled->Name.E; \
- sampled->Name.SE = sampled->Name.V > 0 ? sqrt(sampled->Name.V / N) : 0; \
+ #define SETUP_MC_RESULT(Name, Count) { \
+ const double N = (double)(Count); \
+ struct mc_data* data = &mc->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ sampled->Name.E = weight / N; \
+ sampled->Name.V = sqr_weight/N - sampled->Name.E*sampled->Name.E; \
+ sampled->Name.V = sampled->Name.V > 0 ? sampled->Name.V : 0; \
+ sampled->Name.SE = sqrt(sampled->Name.V / N); \
} (void)0
- SETUP_MC_RESULT(cos_factor);
- SETUP_MC_RESULT(shadowed);
+ SETUP_MC_RESULT(cos_factor, sampled->nb_samples);
+ SETUP_MC_RESULT(shadowed, estimator->realisation_count);
#undef SETUP_MC_RESULT
return RES_OK;
}
diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h
@@ -29,23 +29,77 @@ struct ssol_instance;
/* Monte carlo data */
struct mc_data {
- double weight;
- double sqr_weight;
+ size_t irealisation;
+ double tmp;
+
+ /* Internal data; use get() */
+ double weight__;
+ double sqr_weight__;
};
-#define MC_DATA_NULL__ { 0, 0 }
+#define MC_DATA_NULL__ { SIZE_MAX, 0, 0, 0 }
static const struct mc_data MC_DATA_NULL = MC_DATA_NULL__;
#define MC_RECEIVER_DATA \
- struct mc_data integrated_irradiance; /* In W */ \
- struct mc_data integrated_absorbed_irradiance; /* In W */ \
- struct mc_data absorptivity_loss; /* In W */ \
- struct mc_data reflectivity_loss; /* In W */
+ struct mc_data incoming_flux; /* In W */ \
+ struct mc_data incoming_if_no_atm_loss; /* In W */ \
+ struct mc_data incoming_if_no_field_loss; /* In W */ \
+ struct mc_data incoming_lost_in_field; /* In W */ \
+ struct mc_data incoming_lost_in_atmosphere; /* In W */ \
+ struct mc_data absorbed_flux; /* In W */ \
+ struct mc_data absorbed_if_no_atm_loss; /* In W */ \
+ struct mc_data absorbed_if_no_field_loss; /* In W */ \
+ struct mc_data absorbed_lost_in_field; /* In W */ \
+ struct mc_data absorbed_lost_in_atmosphere; /* In W */
-#define MC_RECEIVER_DATA_NULL__ \
- MC_DATA_NULL__, \
- MC_DATA_NULL__, \
- MC_DATA_NULL__, \
- MC_DATA_NULL__
+/*******************************************************************************
+ * Deferred MC data accumulators
+ ******************************************************************************/
+static INLINE void
+mc_data_init(struct mc_data* data)
+{
+ ASSERT(data);
+ *data = MC_DATA_NULL;
+}
+
+static INLINE void
+mc_data_flush(struct mc_data* data)
+{
+ ASSERT(data);
+ data->weight__ += data->tmp;
+ data->sqr_weight__ += data->tmp * data->tmp;
+ data->tmp = 0;
+}
+
+static INLINE void
+mc_data_add_weight(struct mc_data* data, size_t irealisation, double w)
+{
+ ASSERT(data);
+ ASSERT(irealisation != SIZE_MAX);
+ if(irealisation != data->irealisation) {
+ mc_data_flush(data);
+ data->irealisation = irealisation;
+ }
+ data->tmp += w;
+}
+
+static INLINE void
+mc_data_accum(struct mc_data* dst, struct mc_data* src)
+{
+ ASSERT(dst && src);
+ mc_data_flush(dst);
+ mc_data_flush(src);
+ dst->weight__ += src->weight__;
+ dst->sqr_weight__ += src->sqr_weight__;
+}
+
+static INLINE void
+mc_data_get(struct mc_data* data, double* weight, double* sqr_weight)
+{
+ ASSERT(data && weight && sqr_weight);
+ mc_data_flush(data);
+ *weight = data->weight__;
+ *sqr_weight = data->sqr_weight__;
+}
/*******************************************************************************
* One sided per shape MC data
@@ -54,7 +108,19 @@ struct mc_primitive_1side {
MC_RECEIVER_DATA
};
-#define MC_PRIMITIVE_1SIDE_NULL__ { MC_RECEIVER_DATA_NULL__ }
+#define MC_PRIMITIVE_1SIDE_NULL__ { \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__, \
+ MC_DATA_NULL__ \
+}
+
static const struct mc_primitive_1side MC_PRIMITIVE_1SIDE_NULL =
MC_PRIMITIVE_1SIDE_NULL__;
@@ -142,15 +208,38 @@ struct mc_receiver_1side {
struct htable_shape2mc shape2mc;
};
+static FINLINE void
+mc_receiver_1side_copy_mc_weights__
+ (struct mc_receiver_1side* dst, const struct mc_receiver_1side* src)
+{
+ ASSERT(dst && src);
+ dst->incoming_flux = src->incoming_flux;
+ dst->incoming_if_no_atm_loss = src->incoming_if_no_atm_loss;
+ dst->incoming_if_no_field_loss = src->incoming_if_no_field_loss;
+ dst->incoming_lost_in_atmosphere = src->incoming_lost_in_atmosphere;
+ dst->incoming_lost_in_field = src->incoming_lost_in_field;
+ dst->absorbed_flux = src->absorbed_flux;
+ dst->absorbed_if_no_atm_loss = src->absorbed_if_no_atm_loss;
+ dst->absorbed_if_no_field_loss = src->absorbed_if_no_field_loss;
+ dst->absorbed_lost_in_atmosphere = src->absorbed_lost_in_atmosphere;
+ dst->absorbed_lost_in_field = src->absorbed_lost_in_field;
+}
+
static INLINE void
mc_receiver_1side_init
(struct mem_allocator* allocator, struct mc_receiver_1side* mc)
{
ASSERT(mc);
- mc->integrated_irradiance = MC_DATA_NULL;
- mc->integrated_absorbed_irradiance = MC_DATA_NULL;
- mc->absorptivity_loss = MC_DATA_NULL;
- mc->reflectivity_loss = MC_DATA_NULL;
+ mc->incoming_flux = MC_DATA_NULL;
+ mc->incoming_if_no_atm_loss = MC_DATA_NULL;
+ mc->incoming_if_no_field_loss = MC_DATA_NULL;
+ mc->incoming_lost_in_atmosphere = MC_DATA_NULL;
+ mc->incoming_lost_in_field = MC_DATA_NULL;
+ mc->absorbed_flux = MC_DATA_NULL;
+ mc->absorbed_if_no_atm_loss = MC_DATA_NULL;
+ mc->absorbed_if_no_field_loss = MC_DATA_NULL;
+ mc->absorbed_lost_in_atmosphere = MC_DATA_NULL;
+ mc->absorbed_lost_in_field = MC_DATA_NULL;
htable_shape2mc_init(allocator, &mc->shape2mc);
}
@@ -166,10 +255,7 @@ mc_receiver_1side_copy
(struct mc_receiver_1side* dst, const struct mc_receiver_1side* src)
{
ASSERT(dst && src);
- dst->integrated_irradiance = src->integrated_irradiance;
- dst->integrated_absorbed_irradiance = src->integrated_absorbed_irradiance;
- dst->absorptivity_loss = src->absorptivity_loss;
- dst->reflectivity_loss = src->reflectivity_loss;
+ mc_receiver_1side_copy_mc_weights__(dst, src);
return htable_shape2mc_copy(&dst->shape2mc, &src->shape2mc);
}
@@ -178,10 +264,7 @@ mc_receiver_1side_copy_and_release
(struct mc_receiver_1side* dst, struct mc_receiver_1side* src)
{
ASSERT(dst && src);
- dst->integrated_irradiance = src->integrated_irradiance;
- dst->integrated_absorbed_irradiance = src->integrated_absorbed_irradiance;
- dst->absorptivity_loss = src->absorptivity_loss;
- dst->reflectivity_loss = src->reflectivity_loss;
+ mc_receiver_1side_copy_mc_weights__(dst, src);
return htable_shape2mc_copy_and_release(&dst->shape2mc, &src->shape2mc);
}
@@ -448,11 +531,11 @@ struct ssol_estimator {
/* Implicit MC computations */
struct mc_data cos_factor;
- struct mc_data absorbed;
+ struct mc_data absorbed_by_receivers;
struct mc_data shadowed;
struct mc_data missing;
- struct mc_data atmosphere;
- struct mc_data reflectivity;
+ struct mc_data absorbed_by_atmosphere;
+ struct mc_data other_absorbed;
struct htable_receiver mc_receivers; /* Per receiver MC */
struct htable_sampled mc_sampled; /* Per sampled instance MC */
diff --git a/src/ssol_material.c b/src/ssol_material.c
@@ -64,16 +64,6 @@ ssol_data_ceq(const struct ssol_data* a, const struct ssol_data* b)
return i;
}
-/* Define if the submitted media are *certainly* equals. Refer to the
- * check_ssol_data for more details. */
-static INLINE int
-media_ceq(const struct ssol_medium* a, const struct ssol_medium* b)
-{
- ASSERT(a && b);
- return ssol_data_ceq(&a->refractive_index, &b->refractive_index)
- && ssol_data_ceq(&a->absorption, &b->absorption);
-}
-
static void
shade_normal_default
(struct ssol_device* dev,
@@ -738,3 +728,12 @@ material_get_next_medium
return RES_OK;
}
+/* Define if the submitted media are *certainly* equals. Refer to the
+* check_ssol_data for more details. */
+int
+media_ceq(const struct ssol_medium* a, const struct ssol_medium* b)
+{
+ ASSERT(a && b);
+ return ssol_data_ceq(&a->refractive_index, &b->refractive_index)
+ && ssol_data_ceq(&a->absorption, &b->absorption);
+}
diff --git a/src/ssol_material_c.h b/src/ssol_material_c.h
@@ -92,5 +92,8 @@ material_get_next_medium
const struct ssol_medium* medium, /* Current mediu */
struct ssol_medium* next_medium);
+extern LOCAL_SYM int
+media_ceq(const struct ssol_medium* a, const struct ssol_medium* b);
+
#endif /* SSOL_MATERIAL_C_H */
diff --git a/src/ssol_mc_receiver.c b/src/ssol_mc_receiver.c
@@ -53,15 +53,27 @@ ssol_estimator_get_mc_receiver
mc_rcv1 = side == SSOL_FRONT ? &mc_rcv->front : &mc_rcv->back;
#define SETUP_MC_RESULT(Name) { \
const double N = (double)estimator->realisation_count; \
- const struct mc_data* data = &mc_rcv1->Name; \
- rcv->Name.E = data->weight / N; \
- rcv->Name.V = data->sqr_weight/N - rcv->Name.E*rcv->Name.E; \
- rcv->Name.SE = rcv->Name.V > 0 ? sqrt(rcv->Name.V / N) : 0; \
+ struct mc_data* data = &mc_rcv1->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ rcv->Name.E = weight / N; \
+ rcv->Name.V = sqr_weight/N - rcv->Name.E*rcv->Name.E; \
+ rcv->Name.V = rcv->Name.V > 0 ? rcv->Name.V : 0; \
+ rcv->Name.SE = sqrt(rcv->Name.V / N); \
} (void)0
- SETUP_MC_RESULT(integrated_irradiance);
- SETUP_MC_RESULT(integrated_absorbed_irradiance);
- SETUP_MC_RESULT(absorptivity_loss);
- SETUP_MC_RESULT(reflectivity_loss);
+ #define MC_SETUP_ALL { \
+ SETUP_MC_RESULT(incoming_flux); \
+ SETUP_MC_RESULT(incoming_if_no_atm_loss); \
+ SETUP_MC_RESULT(incoming_if_no_field_loss); \
+ SETUP_MC_RESULT(incoming_lost_in_atmosphere); \
+ SETUP_MC_RESULT(incoming_lost_in_field); \
+ SETUP_MC_RESULT(absorbed_flux); \
+ SETUP_MC_RESULT(absorbed_if_no_atm_loss); \
+ SETUP_MC_RESULT(absorbed_if_no_field_loss); \
+ SETUP_MC_RESULT(absorbed_lost_in_atmosphere); \
+ SETUP_MC_RESULT(absorbed_lost_in_field); \
+ } (void)0
+ MC_SETUP_ALL;
#undef SETUP_MC_RESULT
rcv->mc__ = mc_rcv1;
rcv->N__ = estimator->realisation_count;
@@ -108,10 +120,7 @@ ssol_mc_shape_get_mc_primitive
prim->Name.V = 0; \
prim->Name.SE = 0; \
} (void)0
- SETUP_MC_RESULT(integrated_irradiance);
- SETUP_MC_RESULT(integrated_absorbed_irradiance);
- SETUP_MC_RESULT(absorptivity_loss);
- SETUP_MC_RESULT(reflectivity_loss);
+ MC_SETUP_ALL;
#undef SETUP_MC_RESULT
} else {
struct s3d_attrib attr;
@@ -143,19 +152,20 @@ ssol_mc_shape_get_mc_primitive
#define SETUP_MC_RESULT(Name) { \
const double N = (double)shape->N__; \
- const struct mc_data* data = &mc_prim1->Name; \
- prim->Name.E = data->weight / N; \
- prim->Name.V = data->sqr_weight/N - prim->Name.E*prim->Name.E; \
- prim->Name.SE = prim->Name.V > 0 ? sqrt(prim->Name.V / N) : 0; \
+ struct mc_data* data = &mc_prim1->Name; \
+ double weight, sqr_weight; \
+ mc_data_get(data, &weight, &sqr_weight); \
+ prim->Name.E = weight / N; \
+ prim->Name.V = sqr_weight/N - prim->Name.E*prim->Name.E; \
+ prim->Name.V = prim->Name.V > 0 ? prim->Name.V : 0; \
+ prim->Name.SE = sqrt(prim->Name.V / N); \
prim->Name.E /= area; \
prim->Name.V /= area*area; \
prim->Name.SE /= area; \
} (void)0
- SETUP_MC_RESULT(integrated_irradiance);
- SETUP_MC_RESULT(integrated_absorbed_irradiance);
- SETUP_MC_RESULT(absorptivity_loss);
- SETUP_MC_RESULT(reflectivity_loss);
+ MC_SETUP_ALL;
#undef SETUP_MC_RESULT
+ #undef MC_SETUP_ALL
}
return RES_OK;
diff --git a/src/ssol_scene.c b/src/ssol_scene.c
@@ -85,6 +85,9 @@ ssol_scene_create
res = s3d_scene_create(dev->s3d, &scene->scn_samp);
if(res != RES_OK) goto error;
+ /* default air medium */
+ ssol_medium_copy(&scene->air, &SSOL_MEDIUM_VACUUM);
+
exit:
if(out_scene) *out_scene = scene;
return res;
@@ -231,6 +234,7 @@ ssol_scene_clear(struct ssol_scene* scene)
htable_instance_clear(&scene->instances_samp);
S3D(scene_clear(scene->scn_rt));
S3D(scene_clear(scene->scn_samp));
+ ssol_medium_clear(&scene->air);
if(scene->sun) SSOL(scene_detach_sun(scene, scene->sun));
if(scene->atmosphere) SSOL(scene_detach_atmosphere(scene, scene->atmosphere));
return RES_OK;
@@ -491,7 +495,9 @@ hit_filter_function
return 1;
}
hit_side = d3_dot(dir, N) < 0 ? SSOL_FRONT : SSOL_BACK;
- if(inst == rdata->inst_from && hit_side != rdata->side_from) {
+ if(inst == rdata->inst_from
+ && sshape == rdata->sshape_from
+ && hit_side != rdata->side_from) {
/* The intersected instance is the one from which the ray starts,
* ensure that the ray does not intersect the opposite side of the
* quadric
diff --git a/src/ssol_scene_c.h b/src/ssol_scene_c.h
@@ -45,6 +45,7 @@ struct ssol_scene {
struct ssol_sun* sun; /* Sun of the scene */
struct ssol_atmosphere* atmosphere; /* Atmosphere of the scene */
+ struct ssol_medium air; /* Defined according to atmosphere's properties */
struct ssol_device* dev;
ref_T ref;
diff --git a/src/ssol_solver.c b/src/ssol_solver.c
@@ -53,11 +53,11 @@ struct thread_context {
struct ssp_rng* rng;
struct ssf_bsdf* bsdf;
struct mc_data cos_factor;
- struct mc_data absorbed;
+ struct mc_data absorbed_by_receivers;
struct mc_data shadowed;
struct mc_data missing;
- struct mc_data atmosphere;
- struct mc_data reflectivity;
+ struct mc_data absorbed_by_atmosphere;
+ struct mc_data other_absorbed;
struct htable_receiver mc_rcvs;
struct htable_sampled mc_samps;
struct darray_path paths; /* paths */
@@ -107,11 +107,11 @@ thread_context_copy
dst->rng = src->rng;
dst->bsdf = src->bsdf;
dst->cos_factor = src->cos_factor;
- dst->absorbed = src->absorbed;
+ dst->absorbed_by_receivers = src->absorbed_by_receivers;
dst->shadowed = src->shadowed;
dst->missing = src->missing;
- dst->atmosphere = src->atmosphere;
- dst->reflectivity = src->reflectivity;
+ dst->absorbed_by_atmosphere = src->absorbed_by_atmosphere;
+ dst->other_absorbed = src->other_absorbed;
res = htable_receiver_copy(&dst->mc_rcvs, &src->mc_rcvs);
if(res != RES_OK) return res;
res = htable_sampled_copy(&dst->mc_samps, &src->mc_samps);
@@ -170,12 +170,26 @@ struct point {
double dir[3];
float uv[2];
double wl; /* Sampled wavelength */
- /* MC weights, before and after hit */
- double incoming_weight, weight;
+ const struct ssol_material* material;
+ /* tmp quantities to compute weights */
+ double kabs_at_pt;
+ double partial_atm, partial_recv, partial_other;
+ /* MC weights */
+ /* Set once */
+ double initial_flux; /* the initial flux*/
double cos_factor; /* local cos at the starting point */
- double absorbed_irradiance; /* current hit only */
- double absorptivity_loss_before, absorptivity_loss;
- double reflectivity_loss_before, reflectivity_loss;
+ /* outgoing weights at previous hit */
+ double prev_outgoing_flux;
+ double prev_outgoing_if_no_atm_loss;
+ double prev_outgoing_if_no_field_loss;
+ /* incoming weights at current hit */
+ double incoming_flux;
+ double incoming_if_no_atm_loss;
+ double incoming_if_no_field_loss;
+ /* outgoing weights at current hit */
+ double outgoing_flux;
+ double outgoing_if_no_atm_loss;
+ double outgoing_if_no_field_loss;
enum ssol_side_flag side;
};
@@ -189,7 +203,9 @@ struct point {
{0, 0, 0}, /* Direction */ \
{0, 0}, /* UV */ \
0, /* Wavelength */ \
- 0, 0, 0, 0, 0, 0, 0, 0, /* MC weights */ \
+ NULL, /* Material */ \
+ 0, 0, 0, 0, /* tmp values */ \
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, /* MC weights */ \
SSOL_FRONT /* Side */ \
}
static const struct point POINT_NULL = POINT_NULL__;
@@ -211,6 +227,7 @@ point_init
struct ranst_sun_dir* ran_sun_dir,
struct ranst_sun_wl* ran_sun_wl,
struct ssp_rng* rng,
+ struct ssol_medium* current_medium,
int* is_lit)
{
struct ssol_surface_fragment frag;
@@ -219,6 +236,7 @@ point_init
struct ray_data ray_data = RAY_DATA_NULL;
struct ssol_material* mtl;
double N[3], Np[3];
+ double w0;
float dir[3], pos[3], range[2] = { 0, FLT_MAX };
size_t id;
res_T res = RES_OK;
@@ -277,18 +295,23 @@ point_init
/* Initialise the Monte Carlo weight */
if(pt->sshape->shape->type != SHAPE_PUNCHED) {
double surface_sun_cos = fabs(d3_dot(Np, pt->dir));
- pt->weight = scn->sun->dni * sampled_area_proxy * surface_sun_cos;
+ w0 = scn->sun->dni * sampled_area_proxy * surface_sun_cos;
pt->cos_factor = surface_sun_cos;
} else {
double cos_ratio, surface_proxy_cos, surface_sun_cos;
surface_proxy_cos = fabs(d3_dot(pt->N, Np));
surface_sun_cos = fabs(d3_dot(Np, pt->dir));
cos_ratio = surface_sun_cos / surface_proxy_cos;
- pt->weight = scn->sun->dni * sampled_area_proxy * cos_ratio;
+ w0 = scn->sun->dni * sampled_area_proxy * cos_ratio;
pt->cos_factor = surface_sun_cos;
}
+
+ pt->partial_atm = pt->partial_recv = pt->partial_other = 0;
+ pt->initial_flux = w0;
+ pt->prev_outgoing_flux = w0;
+ pt->prev_outgoing_if_no_atm_loss = w0;
+ pt->prev_outgoing_if_no_field_loss = w0;
d3_set(pt->N, N);
- pt->absorptivity_loss = pt->reflectivity_loss = 0;
ASSERT(d3_dot(pt->N, pt->dir) <= 0);
/* Store sampled entity related weights */
@@ -296,10 +319,29 @@ point_init
if(res != RES_OK) goto error;
pt->mc_samp->nb_samples++;
+ /* Define the medium in which the sampled point lies */
+ pt->material = point_get_material(pt);
+ switch (pt->material->type) {
+ case SSOL_MATERIAL_DIELECTRIC:
+ case SSOL_MATERIAL_THIN_DIELECTRIC:
+ /* TODO: check sampled face role!!! */
+ ssol_medium_copy(current_medium,
+ (pt->side == SSOL_FRONT) ?
+ &pt->material->in_medium : &pt->material->out_medium);
+ break;
+ case SSOL_MATERIAL_MATTE:
+ case SSOL_MATERIAL_MIRROR:
+ case SSOL_MATERIAL_VIRTUAL:
+ ssol_medium_copy(current_medium, &scn->air);
+ break;
+ default: FATAL("Unreachable code\n"); break;
+ }
+
/* Initialise the ray data to avoid self intersection */
ray_data.scn = scn;
ray_data.prim_from = pt->prim;
ray_data.inst_from = pt->inst;
+ ray_data.sshape_from = pt->sshape;
ray_data.side_from = pt->side;
ray_data.discard_virtual_materials = 1; /* Do not intersect virtual mtl */
ray_data.reversed_ray = 1; /* The ray direction is reversed */
@@ -352,7 +394,9 @@ point_update_from_hit
case SHAPE_PUNCHED:
d3_normalize(pt->N, rdata->N);
d3_muld(tmp, pt->dir, rdata->dst);
- d3_add(pt->pos, pt->pos, tmp);
+ f3_set_d3(tmpf, tmp);
+ f3_add(tmpf, org, tmpf);
+ d3_set_f3(pt->pos, tmpf);
break;
default: FATAL("Unreachable code"); break;
}
@@ -366,23 +410,33 @@ point_update_from_hit
pt->side = SSOL_BACK;
d3_minus(pt->N, pt->N); /* Force the normal to look forward dir */
}
+
+ /* Update material */
+ pt->material = point_get_material(pt);
+}
+
+static FINLINE int
+point_is_receiver(const struct point* pt)
+{
+ return (pt->inst->receiver_mask & (int)pt->side) != 0;
}
static FINLINE res_T
point_shade
(struct point* pt,
struct ssf_bsdf* bsdf,
- struct ssol_medium* medium,
+ const struct ssol_medium* in_medium,
+ struct ssol_medium* out_medium,
struct ssp_rng* rng,
double dir[3])
{
struct ssol_material* mtl;
struct ssol_surface_fragment frag;
- double r = 1;
+ double propagated = 0;
double wi[3], N[3], pdf;
int type = 0;
res_T res;
- ASSERT(pt && bsdf && medium && rng && dir);
+ ASSERT(pt && bsdf && in_medium && out_medium && rng && dir);
/* TODO ensure that if `prim' was sampled, then the surface fragment setup
* remains valid in *all* situations, i.e. even though the point primitive
@@ -400,7 +454,7 @@ point_shade
/* Shade the surface fragment */
mtl = point_get_material(pt);
SSF(bsdf_clear(bsdf));
- res = material_setup_bsdf(mtl, &frag, pt->wl, medium, 0, bsdf);
+ res = material_setup_bsdf(mtl, &frag, pt->wl, in_medium, 0, bsdf);
if(res != RES_OK) return res;
/* Perturbe the normal */
@@ -411,44 +465,45 @@ point_shade
d3_minus(wi, pt->dir);
if(d3_dot(wi, N) <= 0) {
- r = 0;
+ propagated = 0;
} else {
double cos_dir_Ng;
- r = ssf_bsdf_sample(bsdf, rng, wi, N, dir, &type, &pdf);
- ASSERT(0 <= r && r <= 1);
+ propagated = ssf_bsdf_sample(bsdf, rng, wi, N, dir, &type, &pdf);
+ ASSERT(0 <= propagated && propagated <= 1);
/* Due to the perturbed normal, the sampled direction may point in the
* wrong direction wrt the sampled BSDF component. */
cos_dir_Ng = d3_dot(frag.Ng, dir);
if((cos_dir_Ng > 0 && (type & SSF_TRANSMISSION))
|| (cos_dir_Ng < 0 && (type & SSF_REFLECTION))) {
- r = 0;
+ propagated = 0;
}
}
- pt->incoming_weight = pt->weight;
- pt->absorptivity_loss_before = pt->absorptivity_loss;
- pt->reflectivity_loss_before = pt->reflectivity_loss;
- pt->absorbed_irradiance = (1 - r) * pt->weight;
- pt->reflectivity_loss += pt->absorbed_irradiance;
- pt->weight = pt->incoming_weight - pt->absorbed_irradiance;
-
- if(type & SSF_TRANSMISSION) material_get_next_medium(mtl, medium, medium);
+ pt->kabs_at_pt = (1 - propagated);
+ pt->outgoing_flux = pt->incoming_flux * propagated;
+ pt->outgoing_if_no_atm_loss = pt->incoming_if_no_atm_loss * propagated;
+ pt->outgoing_if_no_field_loss = point_is_receiver(pt)
+ ? pt->incoming_if_no_field_loss*propagated : pt->incoming_if_no_field_loss;
+
+ if(type & SSF_TRANSMISSION) {
+ material_get_next_medium(mtl, in_medium, out_medium);
+ } else {
+ ssol_medium_copy(out_medium, in_medium);
+ }
return RES_OK;
}
static FINLINE void
-point_hit_virtual(struct point* pt)
-{
- pt->absorbed_irradiance = 0;
- pt->incoming_weight = pt->weight;
- pt->absorptivity_loss_before = pt->absorptivity_loss;
- pt->reflectivity_loss_before = pt->reflectivity_loss;
-}
-
-static FINLINE int
-point_is_receiver(const struct point* pt)
+point_hit_virtual
+ (struct point* pt,
+ const struct ssol_medium* in_medium,
+ struct ssol_medium* out_medium)
{
- return (pt->inst->receiver_mask & (int)pt->side) != 0;
+ pt->kabs_at_pt = 0;
+ pt->outgoing_flux = pt->incoming_flux;
+ pt->outgoing_if_no_atm_loss = pt->incoming_if_no_atm_loss;
+ pt->outgoing_if_no_field_loss = pt->incoming_if_no_field_loss;
+ ssol_medium_copy(out_medium, in_medium);
}
static FINLINE int32_t
@@ -459,31 +514,6 @@ point_get_id(const struct point* pt)
return pt->side == SSOL_FRONT ? (int32_t)inst_id : -(int32_t)inst_id;
}
-static FINLINE res_T
-point_dump
- (const struct point* pt,
- const size_t irealisation,
- const size_t isegment,
- FILE* stream)
-{
- struct ssol_receiver_data out;
- size_t n;
-
- if(!stream) return RES_OK;
-
- out.realization_id = irealisation;
- out.segment_id = (uint32_t)isegment;
- out.receiver_id = point_get_id(pt);
- out.wavelength = (float)pt->wl;
- f3_set_d3(out.pos, pt->pos);
- f3_set_d3(out.in_dir, pt->dir);
- f3_set_d3(out.normal, pt->N);
- f2_set(out.uv, pt->uv);
- out.weight = pt->weight;
- n = fwrite(&out, sizeof(out), 1, stream);
- return n != 1 ? RES_IO_ERR : RES_OK;
-}
-
/*******************************************************************************
* Helper functions
******************************************************************************/
@@ -526,14 +556,21 @@ accum_mc_receivers_1side
res_T res = RES_OK;
ASSERT(dst && src);
- #define ACCUM_WEIGHT(Name) { \
- dst->Name.weight += src->Name.weight; \
- dst->Name.sqr_weight += src->Name.sqr_weight; \
+ #define ACCUM_ALL { \
+ ACCUM_WEIGHT(incoming_flux); \
+ ACCUM_WEIGHT(incoming_if_no_atm_loss); \
+ ACCUM_WEIGHT(incoming_lost_in_field); \
+ ACCUM_WEIGHT(incoming_lost_in_atmosphere); \
+ ACCUM_WEIGHT(incoming_if_no_field_loss); \
+ ACCUM_WEIGHT(absorbed_flux); \
+ ACCUM_WEIGHT(absorbed_if_no_atm_loss); \
+ ACCUM_WEIGHT(absorbed_if_no_field_loss); \
+ ACCUM_WEIGHT(absorbed_lost_in_field); \
+ ACCUM_WEIGHT(absorbed_lost_in_atmosphere); \
} (void)0
- ACCUM_WEIGHT(integrated_irradiance);
- ACCUM_WEIGHT(integrated_absorbed_irradiance);
- ACCUM_WEIGHT(absorptivity_loss);
- ACCUM_WEIGHT(reflectivity_loss);
+
+ #define ACCUM_WEIGHT(Name) mc_data_accum(&dst->Name, &src->Name)
+ ACCUM_ALL;
#undef ACCUM_WEIGHT
/* Merge the per shape MC */
@@ -563,20 +600,16 @@ accum_mc_receivers_1side
res = mc_shape_1side_get_mc_primitive(mc_shape1_dst, iprim, &mc_prim1_dst);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name) { \
- mc_prim1_dst->Name.weight += mc_prim1_src->Name.weight; \
- mc_prim1_dst->Name.sqr_weight += mc_prim1_src->Name.sqr_weight; \
- } (void)0
- ACCUM_WEIGHT(integrated_irradiance);
- ACCUM_WEIGHT(integrated_absorbed_irradiance);
- ACCUM_WEIGHT(absorptivity_loss);
- ACCUM_WEIGHT(reflectivity_loss);
+ #define ACCUM_WEIGHT(Name) \
+ mc_data_accum(&mc_prim1_dst->Name, &mc_prim1_src->Name)
+ ACCUM_ALL;
#undef ACCUM_WEIGHT
htable_prim2mc_iterator_next(&it_prim);
}
htable_shape2mc_iterator_next(&it_shape);
}
+ #undef ACCUM_ALL
exit:
return res;
@@ -594,10 +627,7 @@ accum_mc_sampled(struct mc_sampled* dst, struct mc_sampled* src)
mc_receiver_init(NULL, &mc_rcv_null);
- #define ACCUM_WEIGHT(Name) { \
- dst->Name.weight += src->Name.weight; \
- dst->Name.sqr_weight += src->Name.sqr_weight; \
- } (void)0
+ #define ACCUM_WEIGHT(Name) mc_data_accum(&dst->Name, &src->Name)
ACCUM_WEIGHT(cos_factor);
ACCUM_WEIGHT(shadowed);
#undef ACCUM_WEIGHT
@@ -635,42 +665,51 @@ error:
goto exit;
}
+/*
+ * FIXME Are the following accumulations OK when the radiative path bounces
+ * several times on the same receiver ? It seems weird to add the overall
+ * absorptivity and reflectivity losses to the corresponding per receiver
+ * accumulators since they already registered some losses.
+ */
static res_T
update_mc
- (const struct point* pt,
+ (struct point* pt,
const size_t irealisation,
- const size_t ibounce,
- struct thread_context* thread_ctx,
- FILE* output)
+ struct thread_context* thread_ctx)
{
struct mc_receiver_1side* mc_rcv1 = NULL;
struct mc_receiver_1side* mc_samp_x_rcv1 = NULL;
res_T res = RES_OK;
ASSERT(pt && thread_ctx && point_is_receiver(pt));
- res = point_dump(pt, irealisation, ibounce, output);
- if(res != RES_OK) goto error;
-
- /* Global MC accumulation */
- #define ACCUM_WEIGHT(Res, W) { \
- Res.weight += (W); \
- Res.sqr_weight += (W)*(W); \
- } (void)0
- ACCUM_WEIGHT(thread_ctx->absorbed, pt->absorbed_irradiance);
- #undef ACCUM_WEIGHT
+ pt->partial_recv += pt->incoming_flux - pt->outgoing_flux;
/* Per receiver MC accumulation */
res = get_mc_receiver_1side(&thread_ctx->mc_rcvs, pt->inst, pt->side, &mc_rcv1);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name, W) { \
- mc_rcv1->Name.weight += (W); \
- mc_rcv1->Name.sqr_weight += (W)*(W); \
+ #define ACCUM_ALL { \
+ ACCUM_WEIGHT(incoming_flux, pt->incoming_flux); \
+ ACCUM_WEIGHT(incoming_if_no_atm_loss, pt->incoming_if_no_atm_loss); \
+ ACCUM_WEIGHT(incoming_if_no_field_loss, pt->incoming_if_no_field_loss); \
+ ACCUM_WEIGHT(incoming_lost_in_field, \
+ pt->incoming_if_no_field_loss - pt->incoming_flux); \
+ ACCUM_WEIGHT(incoming_lost_in_atmosphere, \
+ pt->incoming_if_no_atm_loss - pt->incoming_flux); \
+ ACCUM_WEIGHT(absorbed_flux, pt->incoming_flux * pt->kabs_at_pt); \
+ ACCUM_WEIGHT(absorbed_if_no_atm_loss, \
+ pt->incoming_if_no_atm_loss * pt->kabs_at_pt); \
+ ACCUM_WEIGHT(absorbed_if_no_field_loss, \
+ pt->incoming_if_no_field_loss * pt->kabs_at_pt); \
+ ACCUM_WEIGHT(absorbed_lost_in_field, \
+ (pt->incoming_if_no_field_loss - pt->incoming_flux) * pt->kabs_at_pt); \
+ ACCUM_WEIGHT(absorbed_lost_in_atmosphere, \
+ (pt->incoming_if_no_atm_loss - pt->incoming_flux) * pt->kabs_at_pt); \
} (void)0
- ACCUM_WEIGHT(integrated_irradiance, pt->incoming_weight);
- ACCUM_WEIGHT(integrated_absorbed_irradiance, pt->absorbed_irradiance);
- ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss_before);
- ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss_before);
+
+ #define ACCUM_WEIGHT(Name, W) \
+ mc_data_add_weight(&mc_rcv1->Name, irealisation, W)
+ ACCUM_ALL;
#undef ACCUM_WEIGHT
/* Per-sampled/receiver MC accumulation */
@@ -678,14 +717,9 @@ update_mc
(pt->mc_samp, pt->inst, pt->side, &mc_samp_x_rcv1);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name, W) { \
- mc_samp_x_rcv1->Name.weight += (W); \
- mc_samp_x_rcv1->Name.sqr_weight += (W)*(W); \
- } (void)0
- ACCUM_WEIGHT(integrated_irradiance, pt->incoming_weight);
- ACCUM_WEIGHT(integrated_absorbed_irradiance, pt->absorbed_irradiance);
- ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss_before);
- ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss_before);
+ #define ACCUM_WEIGHT(Name, W) \
+ mc_data_add_weight(&mc_samp_x_rcv1->Name, irealisation, W)
+ ACCUM_ALL;
#undef ACCUM_WEIGHT
/* Per primitive receiver MC accumulation */
@@ -699,16 +733,12 @@ update_mc
res = mc_shape_1side_get_mc_primitive(mc_shape1, pt->prim.prim_id, &mc_prim1);
if(res != RES_OK) goto error;
- #define ACCUM_WEIGHT(Name, W) { \
- mc_prim1->Name.weight += (W); \
- mc_prim1->Name.sqr_weight += (W)*(W); \
- } (void)0
- ACCUM_WEIGHT(integrated_irradiance, pt->incoming_weight);
- ACCUM_WEIGHT(integrated_absorbed_irradiance, pt->absorbed_irradiance);
- ACCUM_WEIGHT(absorptivity_loss, pt->absorptivity_loss_before);
- ACCUM_WEIGHT(reflectivity_loss, pt->reflectivity_loss_before);
+ #define ACCUM_WEIGHT(Name, W) \
+ mc_data_add_weight(&mc_prim1->Name, irealisation, W)
+ ACCUM_ALL;
#undef ACCUM_WEIGHT
}
+ #undef ACCUM_ALL
exit:
return res;
@@ -718,7 +748,7 @@ error:
static res_T
trace_radiative_path
- (const size_t path_id, /* Unique id of the radiative path */
+ (const size_t irealisation, /* Unique id of the realisation */
const double sampled_area_proxy, /* Overall area of the sampled geometries */
struct thread_context* thread_ctx,
struct ssol_scene* scn,
@@ -726,11 +756,11 @@ trace_radiative_path
struct s3d_scene_view* view_rt,
struct ranst_sun_dir* ran_sun_dir,
struct ranst_sun_wl* ran_sun_wl,
- const struct ssol_path_tracker* tracker, /* May be NULL */
- FILE* output) /* May be NULL */
+ const struct ssol_path_tracker* tracker) /* May be NULL */
{
struct path path;
- struct ssol_medium medium = SSOL_MEDIUM_VACUUM;
+ struct ssol_medium in_medium = SSOL_MEDIUM_VACUUM;
+ struct ssol_medium out_medium = SSOL_MEDIUM_VACUUM;
struct s3d_hit hit = S3D_HIT_NULL;
struct point pt;
float org[3], dir[3], range[2] = { 0, FLT_MAX };
@@ -744,7 +774,8 @@ trace_radiative_path
/* Find a new starting point of the radiative random walk */
res = point_init(&pt, sampled_area_proxy, scn, &thread_ctx->mc_samps,
- view_samp, view_rt, ran_sun_dir, ran_sun_wl, thread_ctx->rng, &is_lit);
+ view_samp, view_rt, ran_sun_dir, ran_sun_wl, thread_ctx->rng,
+ &in_medium, &is_lit);
if(res != RES_OK) goto error;
if(tracker) {
@@ -759,68 +790,72 @@ trace_radiative_path
}
/* Register the init position onto the sampled geometry */
- res = path_add_vertex(&path, pt.pos, pt.weight);
+ res = path_add_vertex(&path, pt.pos, pt.initial_flux);
if(res != RES_OK) goto error;
}
- #define ACCUM_WEIGHT(Res, W) { \
- Res.weight += (W); \
- Res.sqr_weight += (W)*(W); \
- } (void)0
- ACCUM_WEIGHT(thread_ctx->cos_factor, pt.cos_factor);
- ACCUM_WEIGHT(pt.mc_samp->cos_factor, pt.cos_factor);
+ #define ACCUM_WEIGHT(Res, W) mc_data_add_weight(&Res, irealisation, W)
if(!is_lit) { /* The starting point is not lit */
- ACCUM_WEIGHT(pt.mc_samp->shadowed, pt.weight);
- ACCUM_WEIGHT(thread_ctx->shadowed, pt.weight);
+ ACCUM_WEIGHT(pt.mc_samp->shadowed, pt.initial_flux);
+ ACCUM_WEIGHT(thread_ctx->shadowed, pt.initial_flux);
if(tracker) path.type = SSOL_PATH_SHADOW;
} else {
- struct ssol_material* mtl = point_get_material(&pt);
-
- /* Define the medium in which the sampled point lies */
- switch(mtl->type) {
- case SSOL_MATERIAL_DIELECTRIC:
- case SSOL_MATERIAL_THIN_DIELECTRIC:
- ssol_medium_copy(&medium, &mtl->out_medium);
- break;
- default:
- if(scn->atmosphere)
- ssol_data_copy(&medium.absorption, &scn->atmosphere->absorption);
- break;
- }
-
/* Setup the ray as if it starts from the current point position in order
* to handle the points that start from a virtual material */
f3_set_d3(org, pt.pos);
f3_set_d3(dir, pt.dir);
- hit.distance = 0;
+ hit.distance = 0; /* first loop has no atmospheric absorption */
for(;;) { /* Here we go for the radiative random walk */
+ const int in_atm = media_ceq(&in_medium, &scn->air);
+ const int hit_receiver = point_is_receiver(&pt);
+ const int hit_virtual = pt.material->type == SSOL_MATERIAL_VIRTUAL;
+ int last_segment = 0;
struct ray_data ray_data = RAY_DATA_NULL;
- double absorption;
+ double trans = 1;
+
+ /* Compute medium absorption along the incoming segment. */
+ if(hit.distance > 0) {
+ const double kabs = ssol_data_get_value(&in_medium.absorption, pt.wl);
+ ASSERT(0 <= kabs && kabs <= 1);
+ if(kabs > 0) {
+ trans = exp(-kabs * hit.distance);
+ }
+ }
+ pt.incoming_flux = pt.prev_outgoing_flux * trans;
+ pt.incoming_if_no_atm_loss = in_atm ?
+ pt.prev_outgoing_if_no_atm_loss : pt.prev_outgoing_if_no_atm_loss * trans;
+ pt.incoming_if_no_field_loss = (!in_atm) ?
+ pt.prev_outgoing_if_no_field_loss : pt.prev_outgoing_if_no_field_loss * trans;
/* Compute interaction with material */
- mtl = point_get_material(&pt);
- if(mtl->type == SSOL_MATERIAL_VIRTUAL) {
- point_hit_virtual(&pt);
+ if(hit_virtual) {
+ point_hit_virtual(&pt, &in_medium, &out_medium);
} else {
- /* Modulate the point weight wrt its scattering functions and generate
- * an outgoing direction */
- res = point_shade(&pt, thread_ctx->bsdf, &medium, thread_ctx->rng, pt.dir);
+ /* Modulate the point weights wrt its scattering functions and generate
+ * an outgoing direction and set out_medium accordingly */
+ res = point_shade(&pt, thread_ctx->bsdf, &in_medium, &out_medium,
+ thread_ctx->rng, pt.dir);
if(res != RES_OK) goto error;
}
- if(point_is_receiver(&pt)) {
+ /* If receiver update MC results */
+ if(hit_receiver) {
hit_a_receiver = 1;
- res = update_mc(&pt, path_id, depth, thread_ctx, output);
+ res = update_mc(&pt, irealisation, thread_ctx);
if(res != RES_OK) goto error;
+ } else {
+ pt.partial_other += pt.incoming_flux * pt.kabs_at_pt;
}
- /* Stop the radiative random walk */
- if(pt.weight == 0) break;
+ /* Stop the radiative random walk if no more flux */
+ if(!pt.outgoing_flux) {
+ break;
+ }
/* Setup new ray parameters */
- if(mtl->type == SSOL_MATERIAL_VIRTUAL) {
+ if(hit_virtual) {
/* Note that for Virtual materials, the ray parameters 'org' & 'dir'
* are not updated to ensure that it pursues its traversal without any
* accuracy issue */
@@ -836,46 +871,70 @@ trace_radiative_path
ray_data.scn = scn;
ray_data.prim_from = pt.prim;
ray_data.inst_from = pt.inst;
+ ray_data.sshape_from = pt.sshape;
ray_data.side_from = pt.side;
ray_data.discard_virtual_materials = 0;
ray_data.reversed_ray = 0;
ray_data.range_min = range[0];
ray_data.dst = FLT_MAX;
S3D(scene_view_trace_ray(view_rt, org, dir, range, &ray_data, &hit));
- if(S3D_HIT_NONE(&hit)) {
+ if(S3D_HIT_NONE(&hit)) { /* The ray is lost! */
/* Add the point of the last path segment going to the infinite */
- if (tracker && tracker->infinite_ray_length > 0) {
+ if(tracker && tracker->infinite_ray_length > 0) {
double pos[3], wi[3];
d3_set_f3(wi, dir);
d3_add(pos, pt.pos, d3_muld(wi, wi, tracker->infinite_ray_length));
- res = path_add_vertex(&path, pos, pt.weight);
+ res = path_add_vertex(&path, pos, pt.outgoing_flux);
if (res != RES_OK) goto error;
}
- break;
+ last_segment = 1; /* Path reached its last segment */
+ ASSERT(in_atm);
}
- depth += mtl->type != SSOL_MATERIAL_VIRTUAL;
-
- /* Take into account the medium attenuation */
- absorption = ssol_data_get_value(&medium.absorption, pt.wl);
- if(absorption > 0 && hit.distance > 0) {
- const double transmissivity = exp(-absorption * hit.distance);
- ASSERT(0 < transmissivity && transmissivity <= 1);
- pt.absorptivity_loss += (1 - transmissivity) * pt.weight;
- pt.weight *= transmissivity;
+ /* Don't change prev_outgoing weigths nor record segment absorption until
+ * a non-virtual material is hit or this segment is the last one.
+ * This is because propagation is restarted from the same origin until
+ * a non-virtual material is hit or no further hit can be found. */
+ if(last_segment || !hit_virtual) {
+ if(in_atm) {
+ pt.partial_atm += pt.prev_outgoing_flux - pt.incoming_flux;
+ } else {
+ pt.partial_other += pt.prev_outgoing_flux - pt.incoming_flux;
+ }
+ if(last_segment) {
+ break;
+ }
+ pt.prev_outgoing_flux = pt.outgoing_flux;
+ pt.prev_outgoing_if_no_atm_loss = pt.outgoing_if_no_atm_loss;
+ pt.prev_outgoing_if_no_field_loss = pt.outgoing_if_no_field_loss;
}
+ depth += !hit_virtual;
+ ASSERT(depth < 100); /* FIXME: create a true cancel path for this MC sample */
+
/* Update the point */
point_update_from_hit(&pt, scn, org, dir, &hit, &ray_data);
if(tracker) {
- res = path_add_vertex(&path, pt.pos, pt.weight);
+ res = path_add_vertex(&path, pt.pos, pt.outgoing_flux);
if (res != RES_OK) goto error;
}
+
+ ssol_medium_copy(&in_medium, &out_medium);
}
- ACCUM_WEIGHT(thread_ctx->atmosphere, pt.absorptivity_loss);
- /* all the remaining weight is lost */
- ACCUM_WEIGHT(thread_ctx->missing, pt.weight);
+ /* Check conservation of energy */
+ ASSERT((double)depth * pt.initial_flux * DBL_EPSILON >=
+ fabs(pt.initial_flux -
+ (pt.outgoing_flux + pt.partial_recv + pt.partial_atm + pt.partial_other)));
+
+ /* Now that the sample ends successfully,
+ * record MC weights and register the remaining power as missing */
+ ACCUM_WEIGHT(pt.mc_samp->cos_factor, pt.cos_factor);
+ ACCUM_WEIGHT(thread_ctx->cos_factor, pt.cos_factor);
+ ACCUM_WEIGHT(thread_ctx->missing, pt.outgoing_flux);
+ ACCUM_WEIGHT(thread_ctx->absorbed_by_receivers, pt.partial_recv);
+ ACCUM_WEIGHT(thread_ctx->absorbed_by_atmosphere, pt.partial_atm);
+ ACCUM_WEIGHT(thread_ctx->other_absorbed, pt.partial_other);
#undef ACCUM_WEIGHT
if(tracker) {
@@ -888,7 +947,8 @@ trace_radiative_path
if(res != RES_OK) goto error;
}
exit:
- ssol_medium_clear(&medium);
+ ssol_medium_clear(&in_medium);
+ ssol_medium_clear(&out_medium);
if(tracker) path_release(&path);
return res;
error:
@@ -904,7 +964,6 @@ ssol_solve
struct ssp_rng* rng_state,
const size_t realisations_count,
const struct ssol_path_tracker* path_tracker,
- FILE* output,
struct ssol_estimator** out_estimator)
{
struct htable_receiver_iterator r_it, r_end;
@@ -946,6 +1005,12 @@ ssol_solve
res = scene_check(scn, FUNC_NAME);
if(res != RES_OK) goto error;
+
+ /* init air properties */
+ if(scn->atmosphere)
+ ssol_data_copy(&scn->air.absorption, &scn->atmosphere->absorption);
+ else
+ ssol_data_copy(&scn->air.absorption, &SSOL_MEDIUM_VACUUM.absorption);
/* Create data structures shared by all threads */
res = scene_create_s3d_views(scn, &view_rt, &view_samp, &sampled_area,
@@ -999,7 +1064,7 @@ ssol_solve
/* Execute a MC experiment */
res_local = trace_radiative_path((size_t)i, sampled_area_proxy, thread_ctx,
- scn, view_samp, view_rt, ran_sun_dir, ran_sun_wl, path_tracker, output);
+ scn, view_samp, view_rt, ran_sun_dir, ran_sun_wl, path_tracker);
if(res_local == RES_BAD_OP) {
if(ATOMIC_INCR(&nfailures) >= max_failures) {
@@ -1016,18 +1081,16 @@ ssol_solve
/* Merge per thread global MC estimations */
FOR_EACH(i, 0, nthreads) {
- const struct thread_context* thread_ctx;
- thread_ctx = darray_thread_ctx_cdata_get(&thread_ctxs)+i;
- #define ACCUM_WEIGHT(Name) { \
- estimator->Name.weight += thread_ctx->Name.weight; \
- estimator->Name.sqr_weight += thread_ctx->Name.sqr_weight; \
- } (void)0
+ struct thread_context* thread_ctx;
+ thread_ctx = darray_thread_ctx_data_get(&thread_ctxs)+i;
+ #define ACCUM_WEIGHT(Name) \
+ mc_data_accum(&estimator->Name, &thread_ctx->Name)
ACCUM_WEIGHT(cos_factor);
- ACCUM_WEIGHT(absorbed);
+ ACCUM_WEIGHT(absorbed_by_receivers);
ACCUM_WEIGHT(shadowed);
ACCUM_WEIGHT(missing);
- ACCUM_WEIGHT(atmosphere);
- ACCUM_WEIGHT(reflectivity);
+ ACCUM_WEIGHT(absorbed_by_atmosphere);
+ ACCUM_WEIGHT(other_absorbed);
estimator->realisation_count += thread_ctx->realisation_count;
#undef ACCUM_WEIGHT
}
diff --git a/src/test_ssol_by_receiver_integration.c b/src/test_ssol_by_receiver_integration.c
@@ -134,27 +134,27 @@ main(int argc, char** argv)
#define S_DNI_cos (4 * 1000 * cos(PI / 4))
#define GET_MC_RCV ssol_estimator_get_mc_receiver
#define GET_MC_SAMP_X_RCV ssol_estimator_get_mc_sampled_x_receiver
- CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator1), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator1), RES_OK);
CHECK(GET_MC_RCV(estimator1, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
CHECK(ssol_instance_set_receiver(heliostat, SSOL_FRONT, 0), RES_OK);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S_DNI_cos, S_DNI_cos * 2e-1), 1);
- CHECK(ssol_solve(scene, rng, 8 * N__, 0, NULL, &estimator2), RES_OK);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, S_DNI_cos, S_DNI_cos * 2e-1), 1);
+ CHECK(ssol_solve(scene, rng, 8 * N__, NULL, &estimator2), RES_OK);
CHECK(GET_MC_RCV(estimator2, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S_DNI_cos, S_DNI_cos * 5e-2), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, S_DNI_cos, S_DNI_cos * 5e-2), 1);
CHECK(ssol_estimator_ref_put(estimator1), RES_OK);
- CHECK(ssol_solve(scene, rng, 3 * N__, 0, NULL, &estimator1), RES_OK);
+ CHECK(ssol_solve(scene, rng, 3 * N__, NULL, &estimator1), RES_OK);
CHECK(GET_MC_RCV(estimator1, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S_DNI_cos, S_DNI_cos * 1e-1), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, S_DNI_cos, S_DNI_cos * 1e-1), 1);
CHECK(GET_MC_SAMP_X_RCV(estimator1, heliostat, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(heliostat=>target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S_DNI_cos, S_DNI_cos * 1e-1), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, S_DNI_cos, S_DNI_cos * 1e-1), 1);
/* Free data */
CHECK(ssol_instance_ref_put(heliostat), RES_OK);
diff --git a/src/test_ssol_solver1.c b/src/test_ssol_solver1.c
@@ -86,10 +86,8 @@ main(int argc, char** argv)
double transform2[12]; /* 3x4 column major matrix */
double dbl;
size_t i, count, fcount, scount;
- FILE* tmp = NULL;
double m, std;
double a_m, a_std;
- uint32_t r_id;
unsigned ntris;
(void) argc, (void) argv;
@@ -118,17 +116,18 @@ main(int argc, char** argv)
CHECK(ssol_sun_create_directional(dev, &sun), RES_OK);
CHECK(ssol_sun_set_direction(sun, d3(dir, 1, 0, -1)), RES_OK);
CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK);
- CHECK(ssol_sun_set_dni(sun, 1000), RES_OK);
+#define DNI 1000
+ CHECK(ssol_sun_set_dni(sun, DNI), RES_OK);
CHECK(ssol_scene_create(dev, &scene), RES_OK);
- CHECK(ssol_solve(NULL, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
- CHECK(ssol_solve(scene, NULL, 10, 0, NULL, &estimator), RES_BAD_ARG);
- CHECK(ssol_solve(scene, rng, 0, 0, NULL, &estimator), RES_BAD_ARG);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, NULL), RES_BAD_ARG);
+ CHECK(ssol_solve(NULL, rng, 10, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, NULL, 10, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 0, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 10, NULL, NULL), RES_BAD_ARG);
/* No geometry */
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG);
/* Create scene content */
CHECK(ssol_shape_create_mesh(dev, &dummy), RES_OK);
@@ -164,14 +163,14 @@ main(int argc, char** argv)
CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
/* No sun */
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG);
CHECK(ssol_scene_attach_sun(scene, sun), RES_OK);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_OK);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_ref_put(estimator), RES_OK);
CHECK(ssol_solve
- (scene, rng, 1, &SSOL_PATH_TRACKER_DEFAULT, NULL, &estimator), RES_OK);
+ (scene, rng, 1, &SSOL_PATH_TRACKER_DEFAULT, &estimator), RES_OK);
CHECK(ssol_estimator_get_tracked_paths_count(NULL, NULL), RES_BAD_ARG);
CHECK(ssol_estimator_get_tracked_paths_count(estimator, NULL), RES_BAD_ARG);
@@ -248,7 +247,7 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_instance_sample(secondary, 0), RES_OK);
CHECK(ssol_instance_sample(heliostat, 0), RES_OK);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG);
CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat, &sampled), RES_BAD_ARG);
CHECK(ssol_instance_sample(target, 1), RES_OK);
@@ -257,15 +256,15 @@ main(int argc, char** argv)
/* No attached sun */
CHECK(ssol_scene_detach_sun(scene, sun), RES_OK);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG);
CHECK(ssol_sun_ref_put(sun), RES_OK);
/* Sun with no spectrum */
CHECK(ssol_sun_create_directional(dev, &sun), RES_OK);
CHECK(ssol_sun_set_direction(sun, d3(dir, 1, 0, -1)), RES_OK);
- CHECK(ssol_sun_set_dni(sun, 1000), RES_OK);
+ CHECK(ssol_sun_set_dni(sun, DNI), RES_OK);
CHECK(ssol_scene_attach_sun(scene, sun), RES_OK);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG);
CHECK(ssol_scene_detach_sun(scene, sun), RES_OK);
CHECK(ssol_sun_ref_put(sun), RES_OK);
@@ -274,45 +273,37 @@ main(int argc, char** argv)
CHECK(ssol_sun_set_direction(sun, d3(dir, 1, 0, -1)), RES_OK);
CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK);
CHECK(ssol_scene_attach_sun(scene, sun), RES_OK);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG);
- CHECK(ssol_sun_set_dni(sun, 1000), RES_OK);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG);
+ CHECK(ssol_sun_set_dni(sun, DNI), RES_OK);
/* No receiver in scene */
CHECK(ssol_instance_set_receiver(heliostat, 0, 0), RES_OK);
CHECK(ssol_instance_set_receiver(secondary, 0, 0), RES_OK);
CHECK(ssol_instance_set_receiver(target, 0, 0), RES_OK);
- CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_OK);
+ CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_OK);
CHECK(ssol_instance_set_receiver(heliostat, SSOL_FRONT, 0), RES_OK);
CHECK(ssol_instance_set_receiver(secondary, SSOL_FRONT, 0), RES_OK);
CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 0), RES_OK);
CHECK(ssol_estimator_ref_put(estimator), RES_OK);
/* Can sample any geometry; variance is high */
- NCHECK(tmp = tmpfile(), 0);
#define N__ 10000
#define GET_MC_RCV ssol_estimator_get_mc_receiver
#define GET_MC_GLOBAL ssol_estimator_get_mc_global
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
CHECK(ssol_estimator_get_failed_count(estimator, &fcount), RES_OK);
CHECK(fcount, 0);
- printf("Ir = %g +/- %g; ", m, std);
#define COS cos(PI / 4)
-#define DNI 1000
#define DNI_cos (DNI * COS)
- CHECK(eq_eps(m, 4 * DNI_cos, MMAX(4 * DNI_cos * 1e-2, 2*std)), 1);
+ m = 4 * DNI_cos;
#define SQR(x) ((x)*(x))
dbl = sqrt((SQR(12 * DNI_cos) / 3 - SQR(4 * DNI_cos)) / (double)count);
- CHECK(eq_eps(std, dbl, dbl*1e-2), 1);
+ std = dbl;
/* Target was sampled but shadowed by secondary */
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE);
- printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(eq_eps(mc_global.shadowed.E, m, 2 * dbl), 1);
CHECK(eq_eps(mc_global.missing.E, 2*m, 2*mc_global.missing.SE), 1);
CHECK(GET_MC_RCV(NULL, NULL, SSOL_BACK, NULL), RES_BAD_ARG);
@@ -332,36 +323,30 @@ main(int argc, char** argv)
CHECK(GET_MC_RCV(NULL, target, SSOL_FRONT, &mc_rcv), RES_BAD_ARG);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 2 * std), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-1), 1);
CHECK(ssol_estimator_ref_put(estimator), RES_OK);
/* Sample primary mirror only; variance is low */
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_instance_sample(secondary, 0), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g; ", m, std);
- CHECK(eq_eps(m, 4 * DNI_cos, MMAX(4 * DNI_cos * 1e-2, std)), 1);
- CHECK(eq_eps(std, 0, 1e-4), 1);
+ m = 4 * DNI_cos;
+ std = 0;
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE);
- printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
CHECK(eq_eps(mc_global.missing.E, m, 1e-4), 1);
CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 1e-8), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-4), 1);
CHECK(ssol_estimator_ref_put(estimator), RES_OK);
/* Check atmosphere model; with no absorption result is unchanged */
@@ -371,29 +356,23 @@ main(int argc, char** argv)
CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK);
CHECK(ssol_scene_attach_atmosphere(scene, atm), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g; ", m, std);
- CHECK(eq_eps(m, 4 * DNI_cos, MMAX(4 * DNI_cos * 1e-2, std)), 1);
- CHECK(eq_eps(std, 0, 1e-4), 1);
+ m = 4 * DNI_cos;
+ std = 0;
CHECK(ssol_scene_detach_atmosphere(scene, atm), RES_OK);
CHECK(ssol_atmosphere_ref_put(atm), RES_OK);
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE);
- printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
CHECK(eq_eps(mc_global.missing.E, m, 1e-4), 1);
CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 1e-8), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-4), 1);
CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1);
CHECK(ssol_estimator_ref_put(estimator), RES_OK);
@@ -421,58 +400,30 @@ main(int argc, char** argv)
CHECK(ssol_scene_attach_atmosphere(scene, atm), RES_OK);
CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 1), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &a_m, &a_std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g; ", a_m, a_std);
#define K (exp(-KA * 4 * sqrt(2)))
- CHECK(eq_eps(a_m, REFLECTIVITY * 4 * K * DNI_cos,
- MMAX(4 * K * DNI_cos * 1e-1, a_std)), 1);
- CHECK(eq_eps(a_std, 0, 1e-4), 1);
+ a_m = REFLECTIVITY * 4 * K * DNI_cos;
+ a_std = 0;
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE);
- printf("Atmosphere = %g +/- %g; ", mc_global.atmosphere.E, mc_global.atmosphere.SE);
- printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
- CHECK(eq_eps(mc_global.missing.E + mc_global.atmosphere.E + mc_global.absorbed.E,
- m, 1e-4), 1);
+ CHECK(eq_eps(
+ mc_global.missing.E + mc_global.shadowed.E + mc_global.absorbed_by_receivers.E
+ + mc_global.absorbed_by_atmosphere.E + mc_global.other_absorbed.E,
+ m,
+ 1e-4), 1);
CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
- printf
- ("\tIr(target) = %g +/- %g (%.2g %%)\n",
- mc_rcv.integrated_irradiance.E,
- mc_rcv.integrated_irradiance.SE,
- 100 * mc_rcv.integrated_irradiance.E / m);
- printf
- ("\tAtmospheric Loss(target) = %g +/- %g (%.2g %%)\n",
- mc_rcv.absorptivity_loss.E,
- mc_rcv.absorptivity_loss.SE,
- 100 * mc_rcv.absorptivity_loss.E / m);
- printf
- ("\tReflectivity Loss(target) = %g +/- %g (%.2g %%)\n",
- mc_rcv.reflectivity_loss.E,
- mc_rcv.reflectivity_loss.SE,
- 100 * mc_rcv.reflectivity_loss.E / m);
-
+ PRINT_RCV(mc_rcv);
CHECK(ssol_estimator_get_sampled_count(estimator, &scount), RES_OK);
CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat, &sampled), RES_BAD_ARG);
CHECK(ssol_estimator_get_mc_sampled(estimator, heliostat2, &sampled), RES_OK);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, a_m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, a_std, 1e-4), 1);
- CHECK(eq_eps
- ( mc_rcv.integrated_irradiance.E
- + mc_rcv.absorptivity_loss.E
- + mc_rcv.reflectivity_loss.E, m, 1e-8), 1);
- CHECK(eq_eps
- ( mc_rcv.integrated_irradiance.E
- + mc_rcv.absorptivity_loss.E
- + mc_rcv.reflectivity_loss.E
- + (1 - mc_global.cos_factor.E) * 4 * DNI, 4 * DNI, 1e-8), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, a_m, 1e-4), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, a_std, 1e-4), 1);
+
CHECK(ssol_mc_receiver_get_mc_shape(NULL, NULL, NULL), RES_BAD_ARG);
CHECK(ssol_mc_receiver_get_mc_shape(&mc_rcv, NULL, NULL), RES_BAD_ARG);
CHECK(ssol_mc_receiver_get_mc_shape(NULL, square, NULL), RES_BAD_ARG);
@@ -505,10 +456,10 @@ main(int argc, char** argv)
area = d3_len(d3_cross(N, d3_sub(E1, v1, v0), d3_sub(E2, v2, v0))) * 0.5;
CHECK(ssol_mc_shape_get_mc_primitive(&mc_shape, (unsigned)i, &mc_prim), RES_OK);
- dbl += mc_prim.integrated_irradiance.E * area;
+ dbl += mc_prim.incoming_flux.E * area;
}
- CHECK(eq_eps(dbl, a_m, 1.e-6), 1);
+ CHECK(eq_eps(dbl, a_m, 1e-4), 1);
CHECK(ssol_estimator_ref_put(estimator), RES_OK);
CHECK(ssol_scene_detach_instance(scene, heliostat2), RES_OK);
@@ -523,7 +474,7 @@ main(int argc, char** argv)
CHECK(ssol_sun_create_directional(dev, &sun_mono), RES_OK);
CHECK(ssol_sun_set_direction(sun_mono, d3(dir, 1, 0, -1)), RES_OK);
CHECK(ssol_sun_set_spectrum(sun_mono, spectrum), RES_OK);
- CHECK(ssol_sun_set_dni(sun_mono, 1000), RES_OK);
+ CHECK(ssol_sun_set_dni(sun_mono, DNI), RES_OK);
CHECK(ssol_scene_detach_sun(scene, sun), RES_OK);
CHECK(ssol_scene_attach_sun(scene, sun_mono), RES_OK);
ka[1] = 0.2; ka[0] = ka[2] = 0.1;
@@ -536,28 +487,23 @@ main(int argc, char** argv)
abs_data.value.spectrum = abs_spectrum;
CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g; ", m, std);
#define K2 (exp(-0.121 * 4 * sqrt(2)))
- CHECK(eq_eps(m, 4 * K2 * DNI_cos, MMAX(4 * K2 * DNI_cos * 1e-4, std)), 1);
- CHECK(eq_eps(std, 0, 1e-4), 1);
+ m = 4 * K2 * DNI_cos;
+ std = 0;
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Shadows = %g +/- %g; ", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g; ", mc_global.missing.E, mc_global.missing.SE);
- printf("Cos = %g +/- %g; ", mc_global.cos_factor.E, mc_global.cos_factor.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
CHECK(eq_eps(mc_global.missing.E, m, 1e-4), 1);
CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ PRINT_RCV(mc_rcv);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 1e-4), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-4), 1);
/* Free data */
CHECK(ssol_instance_ref_put(heliostat2), RES_OK);
diff --git a/src/test_ssol_solver2.c b/src/test_ssol_solver2.c
@@ -82,9 +82,7 @@ main(int argc, char** argv)
double transform2[12]; /* 3x4 column major matrix */
double transform3[12]; /* 3x4 column major matrix */
size_t count;
- FILE* tmp;
double m, std;
- uint32_t r_id;
(void) argc, (void) argv;
@@ -177,21 +175,16 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 10000
#define GET_MC_RCV ssol_estimator_get_mc_receiver
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g\n", m, std);
#define COS cos(PI / 4)
#define DNI_cos (1000 * COS)
- CHECK(eq_eps(m, 4 * DNI_cos, 4 * DNI_cos * 1e-4), 1);
+ m = 4 * DNI_cos;
#define SQR(x) ((x)*(x))
- CHECK(eq_eps(std, 0, 1e-4), 1);
+ std = 0;
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
@@ -202,12 +195,12 @@ main(int argc, char** argv)
CHECK(GET_MC_RCV(estimator, heliostat1, SSOL_BACK, &mc_rcv), RES_BAD_ARG);
CHECK(GET_MC_RCV(estimator, secondary, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(secondary) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 1e-8), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-4), 1);
/* Free data */
CHECK(ssol_instance_ref_put(heliostat1), RES_OK);
diff --git a/src/test_ssol_solver2b.c b/src/test_ssol_solver2b.c
@@ -82,9 +82,7 @@ main(int argc, char** argv)
double transform2[12]; /* 3x4 column major matrix */
double transform3[12]; /* 3x4 column major matrix */
size_t count;
- FILE* tmp;
double m, std;
- uint32_t r_id;
(void) argc, (void) argv;
d33_splat(transform1, 0);
@@ -181,21 +179,15 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 50000
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g\n", m, std);
#define COS cos(PI / 4)
#define DNI_cos (1000 * COS)
- CHECK(eq_eps(m, 2 * DNI_cos, MMAX(2 * DNI_cos * 1e-2, std)), 1);
+ m = 2 * DNI_cos;
#define SQR(x) ((x)*(x))
- CHECK(eq_eps(std,
- sqrt((SQR(4 * DNI_cos) / 2 - SQR(2 * DNI_cos)) / (double)count), 1e-3), 1);
+ std = sqrt((SQR(4 * DNI_cos) / 2 - SQR(2 * DNI_cos)) / (double)count);
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
@@ -206,9 +198,9 @@ main(int argc, char** argv)
CHECK(ssol_estimator_get_mc_receiver
(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 2 * std), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-3), 1);
CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
CHECK(count, 0);
diff --git a/src/test_ssol_solver3.c b/src/test_ssol_solver3.c
@@ -72,9 +72,7 @@ main(int argc, char** argv)
double transform[12]; /* 3x4 column major matrix */
double area;
size_t count;
- FILE* tmp;
double m, std;
- uint32_t r_id;
(void) argc, (void) argv;
d3_splat(transform + 9, 0);
@@ -136,21 +134,15 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 20000
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g\n", m, std);
#define COS cos(PI / 4)
#define DNI_cos (1000 * COS)
- CHECK(eq_eps(m, 4 * DNI_cos, 4 * DNI_cos * 2e-1), 1);
+ m = 4 * DNI_cos;
#define SQR(x) ((x)*(x))
- CHECK(eq_eps(std,
- sqrt((SQR(400*DNI_cos) / 100 - SQR(4*DNI_cos)) / (double)count), 20), 1);
+ std = sqrt((SQR(400*DNI_cos) / 100 - SQR(4*DNI_cos)) / (double)count);
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
@@ -161,9 +153,9 @@ main(int argc, char** argv)
CHECK(ssol_estimator_get_mc_receiver
(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 2 * std), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 10), 1);
CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
CHECK(count, 0);
CHECK(ssol_instance_get_area(heliostat, &area), RES_OK);
diff --git a/src/test_ssol_solver4.c b/src/test_ssol_solver4.c
@@ -72,9 +72,7 @@ main(int argc, char** argv)
double dir[3];
double transform[12]; /* 3x4 column major matrix */
size_t count;
- FILE* tmp;
double m1, std1, m2, std2;
- uint32_t r_id1, r_id2;
(void) argc, (void) argv;
#define FOCAL 10
@@ -144,24 +142,17 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target2, 0), RES_OK);
CHECK(ssol_scene_attach_instance(scene, target2), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 100000
#define GET_MC_RCV ssol_estimator_get_mc_receiver
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target1, &r_id1), RES_OK);
- CHECK(ssol_instance_get_id(target2, &r_id2), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id1, count, &m1, &std1), RES_OK);
- CHECK(pp_sum(tmp, (int32_t)r_id2, count, &m2, &std2), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g\n", m1, std1);
#define COS cos(0)
#define DNI_cos (1000 * COS)
- CHECK(eq_eps(m1, 400 * DNI_cos, 400 * DNI_cos * 1e-4), 1);
- CHECK(eq_eps(std1, 0, 1), 1);
- CHECK(m1, m2);
- CHECK(std1, std2);
+ m1 = 400 * DNI_cos;
+ std1 = 0;
+ m2 = m1;
+ std2 = std1;
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
@@ -170,14 +161,14 @@ main(int argc, char** argv)
CHECK(eq_eps(mc_global.missing.E, 400 * DNI_cos, 1e-2), 1); /* nothing absorbed */
CHECK(GET_MC_RCV(estimator, target1, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target1) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m1, 1e-2), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std1, 1e-2), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m1, 1e-2), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std1, 1e-2), 1);
CHECK(GET_MC_RCV(estimator, target2, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target2) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m2, 1e-2), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std2, 1e-2), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m2, 1e-2), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std2, 1e-2), 1);
CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
CHECK(count, 0);
diff --git a/src/test_ssol_solver5.c b/src/test_ssol_solver5.c
@@ -71,9 +71,7 @@ main(int argc, char** argv)
double dir[3];
double transform[12]; /* 3x4 column major matrix */
size_t count;
- FILE* tmp;
double m, std;
- uint32_t r_id;
(void) argc, (void) argv;
#define FOCAL 10
@@ -136,19 +134,14 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 10000
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(pp_sum(tmp, (int32_t)r_id, count, &m, &std), RES_OK);
- CHECK(fclose(tmp), 0);
- printf("Ir = %g +/- %g\n", m, std);
#define COS cos(0)
#define DNI_cos (1000 * COS)
- CHECK(eq_eps(m, 400 * DNI_cos, 20), 1);
- CHECK(eq_eps(std, 0, 1), 1);
+ m = 400 * DNI_cos;
+ std = 0;
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
@@ -158,9 +151,9 @@ main(int argc, char** argv)
CHECK(ssol_estimator_get_mc_receiver
(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, m, 1e-8), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, std, 1e-4), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, m, 1e-8), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, std, 1e-4), 1);
CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
CHECK(count, 0);
diff --git a/src/test_ssol_solver6.c b/src/test_ssol_solver6.c
@@ -78,7 +78,6 @@ main(int argc, char** argv)
struct ssol_mc_receiver mc_rcv;
double dir[3];
double transform[12]; /* 3x4 column major matrix */
- FILE* tmp;
(void) argc, (void) argv;
@@ -175,30 +174,25 @@ main(int argc, char** argv)
transform[11] = 10;
CHECK(ssol_instance_set_transform(target2, transform), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 10000
#define GET_MC_RCV ssol_estimator_get_mc_receiver
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(fclose(tmp), 0);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
- printf("Absorbed = %g +/- %g\n", mc_global.absorbed.E, mc_global.absorbed.SE);
- printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
- printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
- printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(eq_eps(mc_global.shadowed.E, 100000, 2 * 100000/sqrt(N__)), 1);
CHECK(eq_eps(mc_global.missing.E, 0, 0), 1);
CHECK(GET_MC_RCV(estimator, target1, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target1) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, 100000, 2*100000/sqrt(N__)), 1);
- CHECK(mc_rcv.integrated_irradiance.E, mc_rcv.integrated_absorbed_irradiance.E);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, 100000, 2*100000/sqrt(N__)), 1);
+ CHECK(mc_rcv.incoming_flux.E, mc_rcv.absorbed_flux.E);
CHECK(GET_MC_RCV(estimator, target2, SSOL_BACK, &mc_rcv), RES_OK);
printf("Ir(target2) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, 0, 1), 1);
- CHECK(mc_rcv.integrated_irradiance.E, mc_rcv.integrated_absorbed_irradiance.E);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, 0, 1), 1);
+ CHECK(mc_rcv.incoming_flux.E, mc_rcv.absorbed_flux.E);
/* Free data */
CHECK(ssol_instance_ref_put(heliostat1), RES_OK);
diff --git a/src/test_ssol_solver7.c b/src/test_ssol_solver7.c
@@ -82,7 +82,6 @@ main(int argc, char** argv)
struct ssol_mc_receiver mc_rcv;
double dir[3], area;
double transform[12]; /* 3x4 column major matrix */
- FILE* tmp;
/* primary is a parabol */
struct ssol_quadric quadric1 = SSOL_QUADRIC_DEFAULT;
struct ssol_punched_surface punched1 = SSOL_PUNCHED_SURFACE_NULL;
@@ -192,30 +191,24 @@ main(int argc, char** argv)
#define TOTAL (HELIOSTAT_SZ * HELIOSTAT_SZ * DNI_cos)
#define GET_MC_RCV ssol_estimator_get_mc_receiver
- NCHECK(tmp = tmpfile(), 0);
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(fclose(tmp), 0);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK);
CHECK(ssol_estimator_get_sampled_area(estimator, &area), RES_OK);
printf("Total = %g\n", area * DNI_cos);
CHECK(eq_eps(area * DNI_cos, TOTAL, TOTAL * 1e-4), 1);
- printf("Absorbed = %g +/- %g\n", mc_global.absorbed.E, mc_global.absorbed.SE);
- printf("Cos = %g +/- %g\n", mc_global.cos_factor.E, mc_global.cos_factor.SE);
- printf("Shadows = %g +/- %g\n", mc_global.shadowed.E, mc_global.shadowed.SE);
- CHECK(eq_eps(mc_global.shadowed.E, 0, 1e-4), 1);
- printf("Missing = %g +/- %g\n", mc_global.missing.E, mc_global.missing.SE);
+ PRINT_GLOBAL(mc_global);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Abs(target1) = %g +/- %g\n",
- mc_rcv.integrated_absorbed_irradiance.E,
- mc_rcv.integrated_absorbed_irradiance.SE);
+ mc_rcv.absorbed_flux.E,
+ mc_rcv.absorbed_flux.SE);
printf("Ir(target1) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, TOTAL, TOTAL * 1e-4), 1);
- CHECK(eq_eps(mc_rcv.integrated_absorbed_irradiance.E,
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, TOTAL, TOTAL * 1e-4), 1);
+ CHECK(eq_eps(mc_rcv.absorbed_flux.E,
(1 - REFLECTIVITY) * TOTAL, (1 - REFLECTIVITY) *TOTAL * 1e-4), 1);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.SE, 0, 1e-2), 1);
+ CHECK(eq_eps(mc_rcv.incoming_flux.SE, 0, 1e-2), 1);
/* Free data */
CHECK(ssol_instance_ref_put(heliostat), RES_OK);
diff --git a/src/test_ssol_solver8.c b/src/test_ssol_solver8.c
@@ -75,8 +75,6 @@ main(int argc, char** argv)
double dir[3];
double transform[12]; /* 3x4 column major matrix */
size_t count;
- FILE* tmp;
- uint32_t r_id;
(void) argc, (void) argv;
d3_splat(transform + 9, 0);
@@ -139,14 +137,11 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 100000
#define GET_MC_RCV ssol_estimator_get_mc_receiver
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(fclose(tmp), 0);
CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
CHECK(count, 0);
#define S (sqrt(2) * X_SZ * Y_SZ)
@@ -159,9 +154,9 @@ main(int argc, char** argv)
3 * mc_global.missing.SE), 1); /* nothing absorbed */
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target1) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, S * DNI,
- 2 * mc_rcv.integrated_irradiance.SE), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, S * DNI,
+ 2 * mc_rcv.incoming_flux.SE), 1);
/* Free data */
CHECK(ssol_instance_ref_put(heliostat), RES_OK);
diff --git a/src/test_ssol_solver9.c b/src/test_ssol_solver9.c
@@ -77,8 +77,6 @@ main(int argc, char** argv)
double dir[3];
double transform[12]; /* 3x4 column major matrix */
size_t count;
- FILE* tmp;
- uint32_t r_id;
(void) argc, (void) argv;
d3_splat(transform + 9, 0);
@@ -133,14 +131,11 @@ main(int argc, char** argv)
CHECK(ssol_instance_sample(target, 0), RES_OK);
CHECK(ssol_scene_attach_instance(scene, target), RES_OK);
- NCHECK(tmp = tmpfile(), 0);
#define N__ 100000
#define GET_MC_RCV ssol_estimator_get_mc_receiver
- CHECK(ssol_solve(scene, rng, N__, 0, tmp, &estimator), RES_OK);
- CHECK(ssol_instance_get_id(target, &r_id), RES_OK);
+ CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK);
CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK);
CHECK(count, N__);
- CHECK(fclose(tmp), 0);
CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK);
CHECK(count, 0);
#define DNI_TGT_S (DNI * TGT_X * TGT_Y)
@@ -154,9 +149,9 @@ main(int argc, char** argv)
3 * mc_global.missing.SE), 1);
CHECK(GET_MC_RCV(estimator, target, SSOL_FRONT, &mc_rcv), RES_OK);
printf("Ir(target1) = %g +/- %g\n",
- mc_rcv.integrated_irradiance.E, mc_rcv.integrated_irradiance.SE);
- CHECK(eq_eps(mc_rcv.integrated_irradiance.E, MMIN(DNI_S, DNI_TGT_S),
- 3 * mc_rcv.integrated_irradiance.SE), 1);
+ mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE);
+ CHECK(eq_eps(mc_rcv.incoming_flux.E, MMIN(DNI_S, DNI_TGT_S),
+ 3 * mc_rcv.incoming_flux.SE), 1);
/* Free data */
CHECK(ssol_instance_ref_put(heliostat), RES_OK);
diff --git a/src/test_ssol_utils.c b/src/test_ssol_utils.c
@@ -1,62 +0,0 @@
-/* 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/>. */
-
-#define _POSIX_C_SOURCE 200809L /* snprintf support */
-
-#include "test_ssol_utils.h"
-#include "ssol_c.h"
-
-#include <rsys/math.h>
-
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-
-res_T
-pp_sum
- (FILE* f,
- const int32_t receiver_id,
- const size_t count,
- double* mean,
- double* std)
-{
- struct ssol_receiver_data hit;
- double sum = 0;
- double sum2 = 0;
- double E, V, SE;
-
- if(!f || !mean || !std || !count)
- return RES_BAD_ARG;
-
- rewind(f);
- while (1 == fread(&hit, sizeof(struct ssol_receiver_data), 1, f)) {
- if (ferror(f))
- return RES_BAD_ARG;
-
- if (receiver_id != hit.receiver_id)
- continue;
-
- sum += hit.weight;
- sum2 += hit.weight * hit.weight;
- }
-
- E = sum / (double)count;
- V = MMAX(sum2 / (double)count - E*E, 0);
- SE = sqrt(V / (double)count);
-
- *mean = E;
- *std = SE;
- return RES_OK;
-}
diff --git a/src/test_ssol_utils.h b/src/test_ssol_utils.h
@@ -30,12 +30,47 @@ check_memory_allocator(struct mem_allocator* allocator)
}
}
-extern LOCAL_SYM res_T
-pp_sum
- (FILE* f,
- const int32_t receiver_id,
- const size_t count,
- double* mean,
- double* std);
+#define PRINT_GLOBAL(Mc) { \
+ printf("Shadows = %g +/- %g; ", (Mc).shadowed.E, (Mc).shadowed.SE); \
+ printf("Missing = %g +/- %g; ", (Mc).missing.E, (Mc).missing.SE); \
+ printf("Receivers = %g +/- %g; ", \
+ (Mc).absorbed_by_receivers.E, (Mc).absorbed_by_receivers.SE); \
+ printf("Atmosphere = %g +/- %g; ", \
+ (Mc).absorbed_by_atmosphere.E, (Mc).absorbed_by_atmosphere.SE); \
+ printf("Other absorbed = %g +/- %g; ", \
+ (Mc).other_absorbed.E, (Mc).other_absorbed.SE); \
+ printf("Cos = %g +/- %g\n", (Mc).cos_factor.E, (Mc).cos_factor.SE); \
+} (void)0
+
+#define PRINT_RCV(Rcv) { \
+ printf("\tIncoming flux(target) = %g +/- %g \n", \
+ (Rcv).incoming_flux.E, (Rcv).incoming_flux.SE); \
+ printf("\tIncoming flux wo Atmosphere(target) = %g +/- %g (%.2g %%)\n", \
+ (Rcv).incoming_if_no_atm_loss.E, (Rcv).incoming_if_no_atm_loss.SE, \
+ 100 * (Rcv).incoming_if_no_atm_loss.E / (Rcv).incoming_flux.E); \
+ printf("\tIncoming flux wo Field Loss(target) = %g +/- %g (%.2g %%)\n", \
+ (Rcv).incoming_if_no_field_loss.E, (Rcv).incoming_if_no_field_loss.SE, \
+ 100 * (Rcv).incoming_if_no_field_loss.E / (Rcv).incoming_flux.E); \
+ printf("\tAtmospheric Loss on Incoming(target) = %g +/- %g (%.2g %%)\n", \
+ (Rcv).incoming_lost_in_atmosphere.E, (Rcv).incoming_lost_in_atmosphere.SE, \
+ 100 * (Rcv).incoming_lost_in_atmosphere.E / (Rcv).incoming_flux.E); \
+ printf("\tOptical Field Loss(target) on Incoming = %g +/- %g (%.2g %%)\n", \
+ (Rcv).incoming_lost_in_field.E, (Rcv).incoming_lost_in_field.SE, \
+ 100 * (Rcv).incoming_lost_in_field.E / (Rcv).incoming_flux.E); \
+ printf("\tAbsorbed flux(target) = %g +/- %g \n", \
+ (Rcv).absorbed_flux.E, (Rcv).absorbed_flux.SE); \
+ printf("\tAbsorbed flux wo Atmosphere(target) = %g +/- %g (%.2g %%)\n", \
+ (Rcv).absorbed_if_no_atm_loss.E, (Rcv).absorbed_if_no_atm_loss.SE, \
+ 100 * (Rcv).absorbed_if_no_atm_loss.E / (Rcv).absorbed_flux.E); \
+ printf("\tAbsorbed flux wo Field Loss(target) = %g +/- %g (%.2g %%)\n", \
+ (Rcv).absorbed_if_no_field_loss.E, (Rcv).absorbed_if_no_field_loss.SE, \
+ 100 * (Rcv).absorbed_if_no_field_loss.E / (Rcv).absorbed_flux.E); \
+ printf("\tAtmospheric Loss(target) on Absorbed = %g +/- %g (%.2g %%)\n", \
+ (Rcv).absorbed_lost_in_atmosphere.E, (Rcv).absorbed_lost_in_atmosphere.SE, \
+ 100 * (Rcv).absorbed_lost_in_atmosphere.E / (Rcv).absorbed_flux.E); \
+ printf("\tOptical Field Loss(target) on Absorbed = %g +/- %g (%.2g %%)\n", \
+ (Rcv).absorbed_lost_in_field.E, (Rcv).absorbed_lost_in_field.SE, \
+ 100 * (Rcv).absorbed_lost_in_field.E / (Rcv).absorbed_flux.E); \
+} (void)0
#endif /* TEST_SSOL_UTILS_H */
