solstice-solver

Solver library of the solstice app
git clone git://git.meso-star.com/solstice-solver.git
Log | Files | Refs | README | LICENSE

commit c5e9c450a0c4e72ea4477801476165b2133b265d
parent 2a6b0445ced27282e89e8133299826c1ed47231f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 27 Sep 2017 13:43:54 +0200

Merge branch 'release_0.6'

Diffstat:
MREADME.md | 30++++++++++++++++++++++--------
Mcmake/CMakeLists.txt | 3++-
Msrc/ssol.h | 19++++++++++---------
Msrc/ssol_atmosphere.c | 24++++++++++++------------
Msrc/ssol_atmosphere_c.h | 2+-
Msrc/ssol_draw_pt.c | 10+++++-----
Msrc/ssol_estimator.c | 2+-
Msrc/ssol_estimator_c.h | 2+-
Msrc/ssol_material.c | 15++++++++-------
Msrc/ssol_ranst_sun_dir.c | 25++++++++++++++++---------
Msrc/ssol_ranst_sun_dir.h | 2+-
Msrc/ssol_scene.c | 9+++------
Msrc/ssol_scene_c.h | 7++++---
Msrc/ssol_solver.c | 133+++++++++++++++++++++++++++++++++++++++++++++++++------------------------------
Msrc/ssol_sun.c | 10+++++-----
Msrc/ssol_sun_c.h | 2+-
Msrc/test_ssol_atmosphere.c | 32++++++++++++++++----------------
Msrc/test_ssol_by_receiver_integration.c | 6+++---
Msrc/test_ssol_material.c | 14+++++++-------
Msrc/test_ssol_scene.c | 10+++++-----
Msrc/test_ssol_solver1.c | 60++++++++++++++++++++++++++++++------------------------------
Asrc/test_ssol_solver10.c | 202+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/test_ssol_solver2.c | 2+-
Msrc/test_ssol_solver2b.c | 2+-
Msrc/test_ssol_solver3.c | 2+-
Msrc/test_ssol_solver4.c | 2+-
Msrc/test_ssol_solver5.c | 2+-
Msrc/test_ssol_solver6.c | 2+-
Msrc/test_ssol_solver7.c | 2+-
Msrc/test_ssol_solver8.c | 4++--
Msrc/test_ssol_solver9.c | 2+-
Msrc/test_ssol_sun.c | 26+++++++++++++-------------
Msrc/test_ssol_utils.h | 2+-
33 files changed, 461 insertions(+), 206 deletions(-)

diff --git a/README.md b/README.md @@ -26,20 +26,34 @@ variable the install directories of its dependencies. ## Release notes +### Version 0.6 + +- Fix the integration for non parallel sun: the angle between the principal sun + direction and the sampled direction was not correctly taken into account + leading to a wrong initial weight for the optical paths. +- Fix the integration with shapes having perturbed normals: perturbed normals + must be taken into account in the bounces of the optical paths only, not in + the energy computations. +- Fix the distribution of the pillbox sun: the pdf was wrong. +- Fix the `ssol_sun_pillbox_aperture` function and rename it in + `ssol_sun_pillbox_set_theta_max`. The submitted parameter, i.e. `theta_max`, + is the angular radius but was treated as the angular diameter. +- Update the `ssol_solve` API: add a parameter that controls the number of + realisations than can fail before an error occurs. + ### Version 0.5 - Improve the performances up to 50% by optimising the allocation of the BSDF - during the radiative random walk. Performance gains are mainly observed in - situations where the radiative paths are deep, i.e. when they bounce on many - surfaces. + during the optical paths. Performance gains are mainly observed in situations + where the optical paths are deep, i.e. when they bounce on many surfaces. ### Version 0.4.2 -- Energy conservation property might not be ensured when the radiative paths - were fully absorbed. -- Handle infinite radiative paths, i.e. paths that bounces infinitely due to - the material properties and/or numerical inaccuracies. Use a Russian roulette - to stop the radiative random walk without bias. +- Energy conservation property might not be ensured when the optical paths were + fully absorbed. +- Handle infinite optical paths, i.e. paths that bounces infinitely due to the + material properties and/or numerical inaccuracies. Use a Russian roulette to + stop the optical random walk without bias. ### Version 0.4.1 diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -50,7 +50,7 @@ rcmake_append_runtime_dirs(_runtime_dirs RSys Star3D Star3DUT StarCPR StarSF Sta # Configure and define targets ################################################################################ set(VERSION_MAJOR 0) -set(VERSION_MINOR 5) +set(VERSION_MINOR 6) set(VERSION_PATCH 0) set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) @@ -181,6 +181,7 @@ if(NOT NO_TEST) new_test(test_ssol_solver7) new_test(test_ssol_solver8) new_test(test_ssol_solver9) + new_test(test_ssol_solver10) new_test(test_ssol_sun) build_test(test_ssol_draw) diff --git a/src/ssol.h b/src/ssol.h @@ -259,7 +259,7 @@ struct ssol_data { static const struct ssol_data SSOL_DATA_NULL = SSOL_DATA_NULL__; struct ssol_medium { - struct ssol_data absorption; + struct ssol_data extinction; struct ssol_data refractive_index; }; #define SSOL_MEDIUM_VACUUM__ {{SSOL_DATA_REAL, {0}}, {SSOL_DATA_REAL, {1}}} @@ -372,7 +372,7 @@ struct ssol_mc_global { 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 absorbed_by_atmosphere; /* In W */ + struct ssol_mc_result extinguished_by_atmosphere; /* In W */ struct ssol_mc_result other_absorbed; /* In W */ }; #define SSOL_MC_GLOBAL_NULL__ { \ @@ -1040,9 +1040,9 @@ ssol_sun_set_spectrum struct ssol_spectrum* spectrum); SSOL_API res_T -ssol_sun_set_pillbox_aperture +ssol_sun_pillbox_set_theta_max (struct ssol_sun* sun, - const double angle); /* In radian */ + const double theta_max); /* In radian */ SSOL_API res_T ssol_sun_set_buie_param @@ -1052,7 +1052,7 @@ ssol_sun_set_buie_param /******************************************************************************* * Atmosphere API - Describe an atmosphere model. ******************************************************************************/ -/* The atmosphere describes absorption along the light paths */ +/* The atmosphere describes extinction along the light paths */ SSOL_API res_T ssol_atmosphere_create (struct ssol_device* dev, @@ -1067,9 +1067,9 @@ ssol_atmosphere_ref_put (struct ssol_atmosphere* atmosphere); SSOL_API res_T -ssol_atmosphere_set_absorption +ssol_atmosphere_set_extinction (struct ssol_atmosphere* atmosphere, - struct ssol_data* absorption); + struct ssol_data* extinction); /******************************************************************************* * Estimator API - Describe the state of a simulation. @@ -1182,6 +1182,7 @@ ssol_solve (struct ssol_scene* scn, struct ssp_rng* rng, const size_t realisations_count, + const size_t max_failed_count, const struct ssol_path_tracker* tracker, /* NULL<=>Do not record the paths */ struct ssol_estimator** estimator); @@ -1242,7 +1243,7 @@ static FINLINE struct ssol_medium* ssol_medium_clear(struct ssol_medium* medium) { ASSERT(medium); - ssol_data_clear(&medium->absorption); + ssol_data_clear(&medium->extinction); ssol_data_clear(&medium->refractive_index); return medium; } @@ -1251,7 +1252,7 @@ static FINLINE struct ssol_medium* ssol_medium_copy(struct ssol_medium* dst, const struct ssol_medium* src) { ASSERT(dst && src); - ssol_data_copy(&dst->absorption, &src->absorption); + ssol_data_copy(&dst->extinction, &src->extinction); ssol_data_copy(&dst->refractive_index, &src->refractive_index); return dst; } diff --git a/src/ssol_atmosphere.c b/src/ssol_atmosphere.c @@ -33,25 +33,25 @@ atmosphere_release(ref_T* ref) ASSERT(ref); dev = atmosphere->dev; ASSERT(dev && dev->allocator); - ssol_data_clear(&atmosphere->absorption); + ssol_data_clear(&atmosphere->extinction); MEM_RM(dev->allocator, atmosphere); SSOL(device_ref_put(dev)); } static INLINE int -check_absorption(const struct ssol_data* absorption) +check_extinction(const struct ssol_data* extinction) { - if(!absorption) return 0; + if(!extinction) return 0; - /* Check absorptivity in [0, INF) */ - switch(absorption->type) { + /* Check extinction in [0, INF) */ + switch(extinction->type) { case SSOL_DATA_REAL: - if(absorption->value.real < 0 || absorption->value.real > 1) + if(extinction->value.real < 0 || extinction->value.real > 1) return 0; break; case SSOL_DATA_SPECTRUM: - if(!absorption->value.spectrum - || !spectrum_check_data(absorption->value.spectrum, 0, 1)) + if(!extinction->value.spectrum + || !spectrum_check_data(extinction->value.spectrum, 0, 1)) return 0; break; default: FATAL("Unreachable code\n"); break; @@ -115,13 +115,13 @@ ssol_atmosphere_ref_put } res_T -ssol_atmosphere_set_absorption +ssol_atmosphere_set_extinction (struct ssol_atmosphere* atmosphere, - struct ssol_data* absorption) + struct ssol_data* extinction) { - if(!atmosphere || !absorption || !check_absorption(absorption)) + if(!atmosphere || !extinction || !check_extinction(extinction)) return RES_BAD_ARG; - ssol_data_copy(&atmosphere->absorption, absorption); + ssol_data_copy(&atmosphere->extinction, extinction); return RES_OK; } diff --git a/src/ssol_atmosphere_c.h b/src/ssol_atmosphere_c.h @@ -23,7 +23,7 @@ struct ssol_scene; struct ssol_atmosphere { struct ssol_scene* scene_attachment; - struct ssol_data absorption; + struct ssol_data extinction; struct ssol_device* dev; ref_T ref; diff --git a/src/ssol_draw_pt.c b/src/ssol_draw_pt.c @@ -166,11 +166,11 @@ Li(struct ssol_scene* scn, wl = ranst_sun_wl_get(ctx->ran_wl, ctx->rng); if(scn->atmosphere) { - ssol_data_copy(&medium.absorption, &scn->atmosphere->absorption); + ssol_data_copy(&medium.extinction, &scn->atmosphere->extinction); } for(;;) { - double absorption; + double extinction; S3D(scene_view_trace_ray (view, ray_org, ray_dir, ray_range, &ray_data, &hit)); @@ -179,9 +179,9 @@ Li(struct ssol_scene* scn, break; } - absorption = ssol_data_get_value(&medium.absorption, wl); - if(absorption > 0) { - throughput *= exp(-absorption * hit.distance); + extinction = ssol_data_get_value(&medium.extinction, wl); + if(extinction > 0) { + throughput *= exp(-extinction * hit.distance); if(throughput <= 0) break; } diff --git a/src/ssol_estimator.c b/src/ssol_estimator.c @@ -125,7 +125,7 @@ ssol_estimator_get_mc_global SETUP_MC_RESULT(absorbed_by_receivers); SETUP_MC_RESULT(shadowed); SETUP_MC_RESULT(missing); - SETUP_MC_RESULT(absorbed_by_atmosphere); + SETUP_MC_RESULT(extinguished_by_atmosphere); SETUP_MC_RESULT(other_absorbed); #undef SETUP_MC_RESULT return RES_OK; diff --git a/src/ssol_estimator_c.h b/src/ssol_estimator_c.h @@ -556,7 +556,7 @@ struct ssol_estimator { struct mc_data absorbed_by_receivers; struct mc_data shadowed; struct mc_data missing; - struct mc_data absorbed_by_atmosphere; + struct mc_data extinguished_by_atmosphere; struct mc_data other_absorbed; struct htable_receiver mc_receivers; /* Per receiver MC */ diff --git a/src/ssol_material.c b/src/ssol_material.c @@ -233,8 +233,9 @@ create_thin_dielectric_bsdf eta_i = ssol_data_get_value(&mtl->out_medium.refractive_index, wavelength); eta_t = ssol_data_get_value (&mtl->data.thin_dielectric.slab_medium.refractive_index, wavelength); + /* Here extinction is absorption only */ absorption = ssol_data_get_value - (&mtl->data.thin_dielectric.slab_medium.absorption, wavelength); + (&mtl->data.thin_dielectric.slab_medium.extinction, wavelength); thickness = mtl->data.thin_dielectric.thickness; /* Setup the BxDF */ @@ -286,15 +287,15 @@ check_medium(const struct ssol_medium* medium) { if(!medium) return 0; - /* Check absorption in [0, INF) */ - switch(medium->absorption.type) { + /* Check extinction in [0, INF) */ + switch(medium->extinction.type) { case SSOL_DATA_REAL: - if(medium->absorption.value.real < 0) + if(medium->extinction.value.real < 0) return 0; break; case SSOL_DATA_SPECTRUM: - if(!medium->absorption.value.spectrum - || !spectrum_check_data(medium->absorption.value.spectrum, 0, DBL_MAX)) + if(!medium->extinction.value.spectrum + || !spectrum_check_data(medium->extinction.value.spectrum, 0, DBL_MAX)) return 0; break; default: FATAL("Unreachable code\n"); break; @@ -713,5 +714,5 @@ 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); + && ssol_data_ceq(&a->extinction, &b->extinction); } diff --git a/src/ssol_ranst_sun_dir.c b/src/ssol_ranst_sun_dir.c @@ -40,7 +40,7 @@ struct ran_buie_state { }; struct ran_pillbox_state { - double radius; + double sin2_theta_max; double basis[9]; }; @@ -243,10 +243,17 @@ ran_pillbox_get double dir[3]) { double pt[3]; + double phi, sin2_theta, sin_theta, cos_theta; - ASSERT(ran && ran->state.pillbox.radius > 0); - ssp_ran_uniform_disk(rng, ran->state.pillbox.radius, pt); - pt[2] = 1; + ASSERT(ran && 0 <= ran->state.pillbox.sin2_theta_max + && ran->state.pillbox.sin2_theta_max <= 1); + sin2_theta = ssp_rng_uniform_double(rng, 0, ran->state.pillbox.sin2_theta_max); + sin_theta = sqrt(sin2_theta); + cos_theta = sqrt(1 - sin2_theta); + phi = ssp_rng_uniform_double(rng, 0, 2 * PI); + pt[0] = cos(phi) * sin_theta; + pt[1] = sin(phi) * sin_theta; + pt[2] = cos_theta; d33_muld3(dir, ran->state.pillbox.basis, pt); d3_normalize(dir, dir); return dir; @@ -337,15 +344,15 @@ ranst_sun_dir_buie_setup res_T ranst_sun_dir_pillbox_setup (struct ranst_sun_dir* ran, - const double aperture, + const double theta_max, const double dir[3]) { - double radius; - if (!ran || !dir || aperture <= 0 || aperture >= PI ) + double s; + if (!ran || !dir || theta_max <= 0 || theta_max > PI ) return RES_BAD_ARG; - radius = tan(0.5 * aperture); ran->get = ran_pillbox_get; - ran->state.pillbox.radius = radius; + s = sin(theta_max); + ran->state.pillbox.sin2_theta_max = s * s; d33_basis(ran->state.pillbox.basis, dir); return RES_OK; } diff --git a/src/ssol_ranst_sun_dir.h b/src/ssol_ranst_sun_dir.h @@ -51,7 +51,7 @@ ranst_sun_dir_buie_setup extern LOCAL_SYM res_T ranst_sun_dir_pillbox_setup (struct ranst_sun_dir* ran, - const double aperture, /* Apparent size in radians */ + const double theta_max, /* In radians */ const double dir[3]); extern LOCAL_SYM res_T diff --git a/src/ssol_scene.c b/src/ssol_scene.c @@ -346,9 +346,7 @@ res_T scene_create_s3d_views (struct ssol_scene* scn, struct s3d_scene_view** out_view_rt, - struct s3d_scene_view** out_view_samp, - double* out_sampled_area, - double* out_sampled_area_proxy) + struct s3d_scene_view** out_view_samp) { struct htable_instance_iterator it, end; struct s3d_scene_view* view_rt = NULL; @@ -359,7 +357,6 @@ scene_create_s3d_views int has_receiver = 0; res_T res = RES_OK; ASSERT(scn && out_view_rt && out_view_samp); - ASSERT(out_sampled_area && out_sampled_area_proxy); S3D(scene_clear(scn->scn_samp)); htable_instance_clear(&scn->instances_samp); @@ -417,8 +414,8 @@ scene_create_s3d_views exit: *out_view_rt = view_rt; *out_view_samp = view_samp; - *out_sampled_area = sampled_area; - *out_sampled_area_proxy = sampled_area_proxy; + scn->sampled_area = sampled_area; + scn->sampled_area_proxy = sampled_area_proxy; return res; error: S3D(scene_clear(scn->scn_samp)); diff --git a/src/ssol_scene_c.h b/src/ssol_scene_c.h @@ -43,6 +43,9 @@ struct ssol_scene { struct s3d_scene* scn_rt; /* S3D scene to ray trace */ struct s3d_scene* scn_samp; /* S3D scene to sample */ + double sampled_area; /* area of the geometry in scn_rt */ + double sampled_area_proxy; /* area of the geometry in scn_samp */ + 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 */ @@ -57,9 +60,7 @@ extern LOCAL_SYM res_T scene_create_s3d_views (struct ssol_scene* scn, struct s3d_scene_view** view_rt, - struct s3d_scene_view** view_samp, - double* sampled_area, /* Area of the instance set as "samplable" */ - double* sampled_area_proxy); /* Area of the sampled geometries */ + struct s3d_scene_view** view_samp); extern LOCAL_SYM res_T scene_check diff --git a/src/ssol_solver.c b/src/ssol_solver.c @@ -43,9 +43,6 @@ #include <limits.h> #include <omp.h> -/* How many percent of random walk realisations may fail before an error occurs */ -#define MAX_PERCENT_FAILURES 0.01 - /******************************************************************************* * Thread context ******************************************************************************/ @@ -55,7 +52,7 @@ struct thread_context { struct mc_data absorbed_by_receivers; struct mc_data shadowed; struct mc_data missing; - struct mc_data absorbed_by_atmosphere; + struct mc_data extinguished_by_atmosphere; struct mc_data other_absorbed; struct htable_receiver mc_rcvs; struct htable_sampled mc_samps; @@ -97,7 +94,7 @@ thread_context_copy dst->absorbed_by_receivers = src->absorbed_by_receivers; dst->shadowed = src->shadowed; dst->missing = src->missing; - dst->absorbed_by_atmosphere = src->absorbed_by_atmosphere; + dst->extinguished_by_atmosphere = src->extinguished_by_atmosphere; dst->other_absorbed = src->other_absorbed; res = htable_receiver_copy(&dst->mc_rcvs, &src->mc_rcvs); if(res != RES_OK) return res; @@ -209,7 +206,6 @@ point_get_material(const struct point* pt) static res_T point_init (struct point* pt, - const double sampled_area_proxy, struct ssol_scene* scn, struct htable_sampled* sampled, struct s3d_scene_view* view_samp, @@ -220,12 +216,15 @@ point_init struct ssol_medium* current_medium, int* is_lit) { - struct ssol_surface_fragment frag; struct s3d_attrib attr; struct s3d_hit hit; struct ray_data ray_data = RAY_DATA_NULL; - struct ssol_material* mtl; - double N[3], Np[3]; + double N[3]; + double surface_sun_cos; + double surface_sun0_cos; + double sun0_sun_cos; + double surface_proxy_cos; + double cos_ratio; double w0; float dir[3], pos[3], range[2] = { 0, FLT_MAX }; size_t id; @@ -259,7 +258,7 @@ point_init /* Sample a sun direction */ ranst_sun_dir_get(ran_sun_dir, rng, pt->dir); - + /* Sample a wavelength */ pt->wl = ranst_sun_wl_get(ran_sun_wl, rng); @@ -280,25 +279,16 @@ point_init d3_minus(N, N); /* Force the normal to look forward dir */ } - /* Perturb the normal */ - surface_fragment_setup(&frag, pt->pos, pt->dir, N, &pt->prim, pt->uv); - mtl = point_get_material(pt); - material_shade_normal(mtl, &frag, pt->wl, Np); - /* Initialise the Monte Carlo weight */ - if(pt->sshape->shape->type != SHAPE_PUNCHED) { - double surface_sun_cos = fabs(d3_dot(Np, pt->dir)); - 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; - w0 = scn->sun->dni * sampled_area_proxy * cos_ratio; - pt->cos_factor = surface_sun_cos; - } - + surface_sun_cos = d3_dot(N, pt->dir); + surface_sun0_cos = fabs(d3_dot(scn->sun->direction, N)); + sun0_sun_cos = d3_dot(scn->sun->direction, pt->dir); + surface_proxy_cos = + (pt->sshape->shape->type == SHAPE_MESH) ? 1 : fabs(d3_dot(pt->N, N)); + cos_ratio = fabs(surface_sun_cos / (surface_proxy_cos * sun0_sun_cos)); + w0 = scn->sun->dni * scn->sampled_area_proxy * cos_ratio; + pt->cos_factor = scn->sampled_area_proxy / scn->sampled_area + * surface_sun0_cos / surface_proxy_cos; pt->energy_loss = w0; pt->initial_flux = w0; pt->prev_outgoing_flux = w0; @@ -513,6 +503,46 @@ point_get_id(const struct point* pt) /******************************************************************************* * Helper functions ******************************************************************************/ +static INLINE void +check_energy_conservation + (struct ssol_scene* scn, + struct ssol_estimator* estimator, + const int64_t nrealisations) +{ + struct ssol_mc_global global; + double dni; + double dni_s, pot; + double cos, rcv, atm, other, shadow, miss; + double cos_err, rcv_err, atm_err, other_err, shadow_err, miss_err; + double err, max_loss; + ASSERT(scn && estimator); + + if(RES_OK != ssol_estimator_get_mc_global(estimator, &global)) return; + if(RES_OK != ssol_sun_get_dni(scn->sun, &dni)) return; + + /* Fetch data */ + cos = global.cos_factor.E; + rcv = global.absorbed_by_receivers.E; + atm = global.extinguished_by_atmosphere.E; + other = global.other_absorbed.E; + shadow = global.shadowed.E; + miss = global.missing.E; + cos_err = global.cos_factor.SE; + rcv_err = global.absorbed_by_receivers.SE; + atm_err = global.extinguished_by_atmosphere.SE; + other_err = global.other_absorbed.SE; + shadow_err = global.shadowed.SE; + miss_err = global.missing.SE; + + /* Check energy conservation */ + dni_s = dni * scn->sampled_area; + pot = cos * dni_s; + err = dni_s * cos_err + rcv_err + atm_err + other_err + shadow_err + miss_err; + max_loss = 3 * err + (double)nrealisations * pot * DBL_EPSILON; + if(fabs(pot - (rcv + atm + other + shadow + miss)) > max_loss) + FATAL("error: the energy conservation property is not verified\n"); +} + /* Compute an empirical length of the path segment coming from/going to the * infinite, wrt the scene bounding box */ static INLINE double @@ -769,7 +799,7 @@ apply_factor_mc /* Cancel global MC estimations */ mc_data_apply_factor(&thread_ctx->cos_factor, irealisation, factor); - mc_data_apply_factor(&thread_ctx->absorbed_by_atmosphere, irealisation, factor); + mc_data_apply_factor(&thread_ctx->extinguished_by_atmosphere, irealisation, factor); mc_data_apply_factor(&thread_ctx->absorbed_by_receivers, irealisation, factor); mc_data_apply_factor(&thread_ctx->other_absorbed, irealisation, factor); mc_data_apply_factor(&thread_ctx->missing, irealisation, factor); @@ -829,7 +859,6 @@ cancel_mc static res_T trace_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, struct s3d_scene_view* view_samp, @@ -858,7 +887,7 @@ trace_radiative_path roulette_interval = 4 * typical_max_depth; /* First roulette */ /* Find a new starting point of the radiative random walk */ - res = point_init(&pt, sampled_area_proxy, scn, &thread_ctx->mc_samps, + res = point_init(&pt, scn, &thread_ctx->mc_samps, view_samp, view_rt, ran_sun_dir, ran_sun_wl, thread_ctx->rng, &in_medium, &is_lit); if(res != RES_OK) goto error; @@ -891,7 +920,7 @@ trace_radiative_path * 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; /* first loop has no atmospheric absorption */ + hit.distance = 0; /* first loop has no atmospheric extinction */ for(;;) { /* Here we go for the radiative random walk */ const int in_atm = media_ceq(&in_medium, &scn->air); @@ -902,12 +931,12 @@ trace_radiative_path struct ray_data ray_data = RAY_DATA_NULL; double trans = 1; - /* Compute medium absorption along the incoming segment. */ + /* Compute medium extinction 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); + const double k_ext = ssol_data_get_value(&in_medium.extinction, pt.wl); + ASSERT(0 <= k_ext && k_ext <= 1); + if(k_ext > 0) { + trans = exp(-k_ext * hit.distance); } } pt.incoming_flux = pt.prev_outgoing_flux * trans; @@ -990,14 +1019,14 @@ trace_radiative_path } } - /* Don't change prev_outgoing weigths nor record segment absorption until + /* Don't change prev_outgoing weigths nor record segment extinction 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(weight_is_zero || last_segment || !hit_virtual) { const double absorbed = pt.prev_outgoing_flux - pt.incoming_flux; if(in_atm) { - ACCUM_WEIGHT(thread_ctx->absorbed_by_atmosphere, absorbed); + ACCUM_WEIGHT(thread_ctx->extinguished_by_atmosphere, absorbed); } else { ACCUM_WEIGHT(thread_ctx->other_absorbed, absorbed); } @@ -1042,11 +1071,11 @@ trace_radiative_path ssol_medium_copy(&in_medium, &out_medium); } - /* Register the remaining flux as missing */ ACCUM_WEIGHT(thread_ctx->missing, pt.outgoing_flux); pt.energy_loss -= pt.outgoing_flux; + if(tracker) { path.type = hit_a_receiver ? SSOL_PATH_SUCCESS : SSOL_PATH_MISSING; } @@ -1092,6 +1121,7 @@ ssol_solve (struct ssol_scene* scn, struct ssp_rng* rng_state, const size_t realisations_count, + const size_t max_failed_count, const struct ssol_path_tracker* path_tracker, struct ssol_estimator** out_estimator) { @@ -1105,8 +1135,6 @@ ssol_solve struct ssol_estimator* estimator = NULL; struct ssol_path_tracker tracker; struct ssp_rng_proxy* rng_proxy = NULL; - double sampled_area; - double sampled_area_proxy; int64_t nrealisations = 0; int64_t max_failures = 0; int nthreads = 0; @@ -1124,12 +1152,12 @@ ssol_solve /* CL compiler supports OpenMP parallel loop whose indices are signed. The * following line ensures that the unsigned number of realisations does not * overflow the realisation index. */ - if(realisations_count > INT64_MAX) { + if(realisations_count > INT64_MAX || max_failed_count > INT64_MAX) { res = RES_BAD_ARG; goto error; } nrealisations = (int64_t)realisations_count; - max_failures = (int64_t)((double)nrealisations * MAX_PERCENT_FAILURES); + max_failures = (int64_t)max_failed_count; nthreads = (int)scn->dev->nthreads; res = scene_check(scn, FUNC_NAME); @@ -1137,13 +1165,12 @@ ssol_solve /* init air properties */ if(scn->atmosphere) - ssol_data_copy(&scn->air.absorption, &scn->atmosphere->absorption); + ssol_data_copy(&scn->air.extinction, &scn->atmosphere->extinction); else - ssol_data_copy(&scn->air.absorption, &SSOL_MEDIUM_VACUUM.absorption); + ssol_data_copy(&scn->air.extinction, &SSOL_MEDIUM_VACUUM.extinction); /* Create data structures shared by all threads */ - res = scene_create_s3d_views(scn, &view_rt, &view_samp, &sampled_area, - &sampled_area_proxy); + res = scene_create_s3d_views(scn, &view_rt, &view_samp); if(res != RES_OK) goto error; res = sun_create_direction_distribution(scn->sun, &ran_sun_dir); if(res != RES_OK) goto error; @@ -1192,7 +1219,7 @@ ssol_solve thread_ctx = darray_thread_ctx_data_get(&thread_ctxs) + ithread; /* Execute a MC experiment */ - res_local = trace_radiative_path((size_t)i, sampled_area_proxy, thread_ctx, + res_local = trace_radiative_path((size_t)i, thread_ctx, scn, view_samp, view_rt, ran_sun_dir, ran_sun_wl, path_tracker); if(res_local != RES_OK) { /* Cancel partial MC results */ @@ -1221,7 +1248,7 @@ ssol_solve ACCUM_WEIGHT(absorbed_by_receivers); ACCUM_WEIGHT(shadowed); ACCUM_WEIGHT(missing); - ACCUM_WEIGHT(absorbed_by_atmosphere); + ACCUM_WEIGHT(extinguished_by_atmosphere); ACCUM_WEIGHT(other_absorbed); estimator->realisation_count += thread_ctx->realisation_count; #undef ACCUM_WEIGHT @@ -1292,9 +1319,13 @@ ssol_solve } } - estimator->sampled_area = sampled_area; + estimator->sampled_area = scn->sampled_area; if(mt_res != RES_OK) res = (res_T)mt_res; + #ifndef NDEBUG + check_energy_conservation(scn, estimator, nrealisations); + #endif + exit: darray_thread_ctx_release(&thread_ctxs); if(view_rt) S3D(scene_view_ref_put(view_rt)); diff --git a/src/ssol_sun.c b/src/ssol_sun.c @@ -170,14 +170,14 @@ ssol_sun_set_spectrum(struct ssol_sun* sun, struct ssol_spectrum* spectrum) } res_T -ssol_sun_set_pillbox_aperture(struct ssol_sun* sun, const double angle) +ssol_sun_pillbox_set_theta_max(struct ssol_sun* sun, const double theta_max) { if(!sun - || angle <= 0 - || angle > PI * 0.5 + || theta_max <= 0 + || theta_max > PI || sun->type != SUN_PILLBOX) return RES_BAD_ARG; - sun->data.pillbox.aperture = angle; + sun->data.pillbox.theta_max = theta_max; return RES_OK; } @@ -214,7 +214,7 @@ sun_create_direction_distribution break; case SUN_PILLBOX: res = ranst_sun_dir_pillbox_setup - (ran_dir, sun->data.pillbox.aperture, sun->direction); + (ran_dir, sun->data.pillbox.theta_max, sun->direction); break; case SUN_BUIE: res = ranst_sun_dir_buie_setup diff --git a/src/ssol_sun_c.h b/src/ssol_sun_c.h @@ -31,7 +31,7 @@ enum sun_type { }; struct pillbox { - double aperture; + double theta_max; }; struct buie { diff --git a/src/test_ssol_atmosphere.c b/src/test_ssol_atmosphere.c @@ -32,7 +32,7 @@ main(int argc, char** argv) { struct mem_allocator allocator; struct ssol_device* dev; - struct ssol_data absorption, absorption2; + struct ssol_data extinction, extinction2; struct ssol_atmosphere* atm; struct ssol_spectrum* spectrum; (void) argc, (void) argv; @@ -42,11 +42,11 @@ main(int argc, char** argv) CHECK(ssol_device_create (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev), RES_OK); - absorption2.type = SSOL_DATA_SPECTRUM; + extinction2.type = SSOL_DATA_SPECTRUM; CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK); CHECK(ssol_spectrum_setup(spectrum, get_wlen, 2, NULL), RES_OK); - absorption2.type = SSOL_DATA_SPECTRUM; - absorption2.value.spectrum = spectrum; + extinction2.type = SSOL_DATA_SPECTRUM; + extinction2.value.spectrum = spectrum; CHECK(ssol_atmosphere_create(NULL, &atm), RES_BAD_ARG); CHECK(ssol_atmosphere_create(dev, NULL), RES_BAD_ARG); @@ -58,21 +58,21 @@ main(int argc, char** argv) CHECK(ssol_atmosphere_ref_put(NULL), RES_BAD_ARG); CHECK(ssol_atmosphere_ref_put(atm), RES_OK); - absorption.type = SSOL_DATA_REAL; - absorption.value.real = 0.1; - CHECK(ssol_atmosphere_set_absorption(NULL, &absorption), RES_BAD_ARG); - CHECK(ssol_atmosphere_set_absorption(atm, NULL), RES_BAD_ARG); - CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_OK); - CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_OK); - CHECK(ssol_atmosphere_set_absorption(atm, &absorption2), RES_OK); + extinction.type = SSOL_DATA_REAL; + extinction.value.real = 0.1; + CHECK(ssol_atmosphere_set_extinction(NULL, &extinction), RES_BAD_ARG); + CHECK(ssol_atmosphere_set_extinction(atm, NULL), RES_BAD_ARG); + CHECK(ssol_atmosphere_set_extinction(atm, &extinction), RES_OK); + CHECK(ssol_atmosphere_set_extinction(atm, &extinction), RES_OK); + CHECK(ssol_atmosphere_set_extinction(atm, &extinction2), RES_OK); - /* absorption values out of range */ - absorption.value.real = 2; - CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_BAD_ARG); + /* extinction values out of range */ + extinction.value.real = 2; + CHECK(ssol_atmosphere_set_extinction(atm, &extinction), RES_BAD_ARG); CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK); - CHECK(ssol_atmosphere_set_absorption(atm, &absorption2), RES_BAD_ARG); + CHECK(ssol_atmosphere_set_extinction(atm, &extinction2), RES_BAD_ARG); - CHECK(ssol_spectrum_ref_put(absorption2.value.spectrum), RES_OK); + CHECK(ssol_spectrum_ref_put(extinction2.value.spectrum), RES_OK); CHECK(ssol_device_ref_put(dev), RES_OK); CHECK(ssol_atmosphere_ref_put(atm), RES_OK); diff --git a/src/test_ssol_by_receiver_integration.c b/src/test_ssol_by_receiver_integration.c @@ -134,7 +134,7 @@ 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__, NULL, &estimator1), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator1), RES_OK); CHECK(GET_MC_RCV(estimator1, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE); @@ -142,13 +142,13 @@ main(int argc, char** argv) CHECK(eq_eps (mc_rcv.incoming_flux.E, S_DNI_cos, mc_rcv.incoming_flux.SE*3), 1); - CHECK(ssol_solve(scene, rng, 8 * N__, NULL, &estimator2), RES_OK); + CHECK(ssol_solve(scene, rng, 8 * N__, 0, NULL, &estimator2), RES_OK); CHECK(GET_MC_RCV(estimator2, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE); CHECK(eq_eps(mc_rcv.incoming_flux.E, S_DNI_cos, mc_rcv.incoming_flux.SE*3), 1); CHECK(ssol_estimator_ref_put(estimator1), RES_OK); - CHECK(ssol_solve(scene, rng, 3 * N__, NULL, &estimator1), RES_OK); + CHECK(ssol_solve(scene, rng, 3 * N__, 0, NULL, &estimator1), RES_OK); CHECK(GET_MC_RCV(estimator1, target, SSOL_FRONT, &mc_rcv), RES_OK); printf("Ir(target) = %g +/- %g\n", mc_rcv.incoming_flux.E, mc_rcv.incoming_flux.SE); diff --git a/src/test_ssol_material.c b/src/test_ssol_material.c @@ -162,17 +162,17 @@ test_thin_dielectric(struct ssol_device* dev) CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); shader.normal = get_shader_normal; - ssol_data_set_real(&mdm0.absorption, -1); + ssol_data_set_real(&mdm0.extinction, -1); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); - ssol_data_copy(&mdm0.absorption, &SSOL_MEDIUM_VACUUM.absorption); + ssol_data_copy(&mdm0.extinction, &SSOL_MEDIUM_VACUUM.extinction); ssol_data_set_real(&mdm0.refractive_index, 0); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); ssol_data_copy(&mdm0.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); - ssol_data_set_real(&mdm1.absorption, -1); + ssol_data_set_real(&mdm1.extinction, -1); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); - ssol_data_copy(&mdm1.absorption, &SSOL_MEDIUM_VACUUM.absorption); + ssol_data_copy(&mdm1.extinction, &SSOL_MEDIUM_VACUUM.extinction); ssol_data_set_real(&mdm1.refractive_index, 0); CHECK(ssol_thin_dielectric_setup(mtl, &shader, &mdm0, &mdm1, 1), RES_BAD_ARG); @@ -229,11 +229,11 @@ test_dielectric(struct ssol_device* dev) CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); ssol_data_copy(&mdm1.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); - ssol_data_set_real(&mdm0.absorption, -1); + ssol_data_set_real(&mdm0.extinction, -1); CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); - ssol_data_copy(&mdm0.absorption, &SSOL_MEDIUM_VACUUM.refractive_index); + ssol_data_copy(&mdm0.extinction, &SSOL_MEDIUM_VACUUM.refractive_index); - ssol_data_set_real(&mdm1.absorption, -1); + ssol_data_set_real(&mdm1.extinction, -1); CHECK(ssol_dielectric_setup(NULL, &dielectric, &mdm0, &mdm1), RES_BAD_ARG); ssol_data_copy(&mdm1.refractive_index, &SSOL_MEDIUM_VACUUM.refractive_index); diff --git a/src/test_ssol_scene.c b/src/test_ssol_scene.c @@ -78,7 +78,7 @@ main(int argc, char** argv) struct ssol_scene* scene2; struct ssol_atmosphere* atm; struct ssol_atmosphere* atm2; - struct ssol_data absorption; + struct ssol_data extinction; struct ssol_vertex_data vdata; struct scene_ctx ctx; struct desc desc; @@ -169,11 +169,11 @@ main(int argc, char** argv) CHECK(ssol_scene_detach_sun(scene2, sun), RES_OK); CHECK(ssol_atmosphere_create(dev, &atm), RES_OK); - absorption.type = SSOL_DATA_REAL; - absorption.value.real = 0.1; - CHECK(ssol_atmosphere_set_absorption(atm, &absorption), RES_OK); + extinction.type = SSOL_DATA_REAL; + extinction.value.real = 0.1; + CHECK(ssol_atmosphere_set_extinction(atm, &extinction), RES_OK); CHECK(ssol_atmosphere_create(dev, &atm2), RES_OK); - CHECK(ssol_atmosphere_set_absorption(atm2, &absorption), RES_OK); + CHECK(ssol_atmosphere_set_extinction(atm2, &extinction), RES_OK); CHECK(ssol_scene_attach_atmosphere(NULL, atm), RES_BAD_ARG); CHECK(ssol_scene_attach_atmosphere(scene, NULL), RES_BAD_ARG); diff --git a/src/test_ssol_solver1.c b/src/test_ssol_solver1.c @@ -67,7 +67,7 @@ main(int argc, char** argv) struct ssol_sun* sun_mono; struct ssol_spectrum* spectrum; struct ssol_spectrum* abs_spectrum; - struct ssol_data abs_data; + struct ssol_data extinction; struct ssol_atmosphere* atm; struct ssol_estimator* estimator; struct ssol_mc_sampled sampled; @@ -120,14 +120,14 @@ main(int argc, char** argv) CHECK(ssol_sun_set_dni(sun, DNI), RES_OK); CHECK(ssol_scene_create(dev, &scene), RES_OK); - 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); + 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); /* No geometry */ - CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG); + CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); /* Create scene content */ CHECK(ssol_shape_create_mesh(dev, &dummy), RES_OK); @@ -163,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, NULL, &estimator), RES_BAD_ARG); + CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); CHECK(ssol_scene_attach_sun(scene, sun), RES_OK); - CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_ref_put(estimator), RES_OK); CHECK(ssol_solve - (scene, rng, 1, &SSOL_PATH_TRACKER_DEFAULT, &estimator), RES_OK); + (scene, rng, 1, 0, &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); @@ -247,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, NULL, &estimator), RES_BAD_ARG); + CHECK(ssol_solve(scene, rng, 10, 0, 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); @@ -256,7 +256,7 @@ main(int argc, char** argv) /* No attached sun */ CHECK(ssol_scene_detach_sun(scene, sun), RES_OK); - CHECK(ssol_solve(scene, rng, 10, NULL, &estimator), RES_BAD_ARG); + CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); CHECK(ssol_sun_ref_put(sun), RES_OK); /* Sun with no spectrum */ @@ -264,7 +264,7 @@ main(int argc, char** argv) CHECK(ssol_sun_set_direction(sun, d3(dir, 1, 0, -1)), 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, NULL, &estimator), RES_BAD_ARG); + CHECK(ssol_solve(scene, rng, 10, 0, NULL, &estimator), RES_BAD_ARG); CHECK(ssol_scene_detach_sun(scene, sun), RES_OK); CHECK(ssol_sun_ref_put(sun), RES_OK); @@ -273,14 +273,14 @@ 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, NULL, &estimator), RES_BAD_ARG); + CHECK(ssol_solve(scene, rng, 10, 0, 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, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, 10, 0, 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); @@ -290,7 +290,7 @@ main(int argc, char** argv) #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__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(ssol_estimator_get_failed_count(estimator, &fcount), RES_OK); @@ -332,7 +332,7 @@ main(int argc, char** argv) CHECK(ssol_instance_sample(target, 0), RES_OK); CHECK(ssol_instance_sample(secondary, 0), RES_OK); - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); m = 4 * DNI_cos; @@ -349,14 +349,14 @@ main(int argc, char** argv) 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 */ + /* Check atmosphere model; with no extinction result is unchanged */ CHECK(ssol_atmosphere_create(dev, &atm), RES_OK); - abs_data.type = SSOL_DATA_REAL; - abs_data.value.real = 0; - CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK); + extinction.type = SSOL_DATA_REAL; + extinction.value.real = 0; + CHECK(ssol_atmosphere_set_extinction(atm, &extinction), RES_OK); CHECK(ssol_scene_attach_atmosphere(scene, atm), RES_OK); - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); m = 4 * DNI_cos; @@ -392,15 +392,15 @@ main(int argc, char** argv) CHECK(ssol_scene_attach_instance(scene, heliostat2), RES_OK); #define KA 0.03 - abs_data.value.real = KA; + extinction.value.real = KA; CHECK(ssol_spectrum_create(dev, &abs_spectrum), RES_OK); CHECK(ssol_spectrum_setup(abs_spectrum, get_wlen, 3, &desc), RES_OK); CHECK(ssol_atmosphere_create(dev, &atm), RES_OK); - CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK); + CHECK(ssol_atmosphere_set_extinction(atm, &extinction), RES_OK); CHECK(ssol_scene_attach_atmosphere(scene, atm), RES_OK); CHECK(ssol_instance_set_receiver(target, SSOL_FRONT, 1), RES_OK); - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); #define K (exp(-KA * 4 * sqrt(2))) @@ -411,7 +411,7 @@ main(int argc, char** argv) CHECK(eq_eps(mc_global.shadowed.E, 0, 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, + + mc_global.extinguished_by_atmosphere.E + mc_global.other_absorbed.E, m, 1e-4), 1); CHECK(eq_eps(mc_global.cos_factor.E, COS, 1e-4), 1); @@ -483,11 +483,11 @@ main(int argc, char** argv) desc.count = 3; CHECK(ssol_spectrum_setup(abs_spectrum, get_wlen, 3, &desc), RES_OK); CHECK(ssol_spectrum_setup(abs_spectrum, get_wlen, 2, &desc), RES_OK); - abs_data.type = SSOL_DATA_SPECTRUM; - abs_data.value.spectrum = abs_spectrum; - CHECK(ssol_atmosphere_set_absorption(atm, &abs_data), RES_OK); + extinction.type = SSOL_DATA_SPECTRUM; + extinction.value.spectrum = abs_spectrum; + CHECK(ssol_atmosphere_set_extinction(atm, &extinction), RES_OK); - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); #define K2 (exp(-0.121 * 4 * sqrt(2))) diff --git a/src/test_ssol_solver10.c b/src/test_ssol_solver10.c @@ -0,0 +1,202 @@ +/* Copyright (C) CNRS 2016-2017 + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. */ + +#include "ssol.h" +#include "test_ssol_utils.h" +#define REFLECTIVITY 0 +#include "test_ssol_materials.h" + +#include <rsys/math.h> + +#define X_SZ 10 +#define Y_SZ 10 +#define PLANE_NAME SQUARE +#define HALF_X (X_SZ / 2) +#define HALF_Y (Y_SZ / 2) +STATIC_ASSERT((HALF_X * 2 == X_SZ), ONLY_ENVEN_VALUES_FOR_SQUARE_X); +STATIC_ASSERT((HALF_Y * 2 == Y_SZ), ONLY_ENVEN_VALUES_FOR_SQUARE_Y); +#include "test_ssol_rect_geometry.h" + +#define POLYGON_NAME POLY +#define HALF_X (X_SZ / 2) +#define HALF_Y (Y_SZ / 2) +STATIC_ASSERT((HALF_X * 2 == X_SZ), ONLY_ENVEN_VALUES_FOR_X_SZ); +STATIC_ASSERT((HALF_Y * 2 == Y_SZ), ONLY_ENVEN_VALUES_FOR_Y_SZ); +#include "test_ssol_rect2D_geometry.h" + +#include <rsys/double33.h> + +#include <star/s3d.h> +#include <star/ssp.h> + +static void +get_wlen(const size_t i, double* wlen, double* data, void* ctx) +{ + double wavelengths[3] = { 1, 2, 3 }; + double intensities[3] = { 1, 0.8, 1 }; + CHECK(i < 3, 1); + (void) ctx; + *wlen = wavelengths[i]; + *data = intensities[i]; +} + +int +main(int argc, char** argv) +{ + struct mem_allocator allocator; + struct ssol_device* dev; + struct ssp_rng* rng; + struct ssol_scene* scene; + struct ssol_shape* square; + struct ssol_shape* parabol; + struct ssol_vertex_data attribs[1] = { SSOL_VERTEX_DATA_NULL__ }; + struct ssol_material* m_mtl; + struct ssol_matte_shader shader = SSOL_MATTE_SHADER_NULL; + struct ssol_object* m_object; + struct ssol_object* q_object; + struct ssol_carving carving = SSOL_CARVING_NULL; + struct ssol_quadric quadric = SSOL_QUADRIC_DEFAULT; + struct ssol_punched_surface punched = SSOL_PUNCHED_SURFACE_NULL; + struct ssol_instance* geom1; + struct ssol_instance* geom2; + struct ssol_sun* sun; + struct ssol_spectrum* spectrum; + struct ssol_estimator* estimator; + struct ssol_mc_global mc_global1; + struct ssol_mc_global mc_global2; + struct ssol_mc_receiver mc_rcv1; + struct ssol_mc_receiver mc_rcv2; + double dir[3]; + double transform[12]; /* 3x4 column major matrix */ + size_t count; + + (void) argc, (void) argv; + mem_init_proxy_allocator(&allocator, &mem_default_allocator); + + d33_set_identity(transform); + + CHECK(ssol_device_create + (NULL, &allocator, SSOL_NTHREADS_DEFAULT, 0, &dev), RES_OK); + +#define DNI 1000 + CHECK(ssp_rng_create(&allocator, &ssp_rng_threefry, &rng), RES_OK); + CHECK(ssol_spectrum_create(dev, &spectrum), RES_OK); + CHECK(ssol_spectrum_setup(spectrum, get_wlen, 3, NULL), RES_OK); + CHECK(ssol_sun_create_pillbox(dev, &sun), RES_OK); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 1), RES_OK); + CHECK(ssol_sun_set_direction(sun, d3(dir, 0, 0, -1)), RES_OK); + CHECK(ssol_sun_set_spectrum(sun, spectrum), RES_OK); + CHECK(ssol_sun_set_dni(sun, DNI), RES_OK); + CHECK(ssol_scene_create(dev, &scene), RES_OK); + CHECK(ssol_scene_attach_sun(scene, sun), RES_OK); + + /* Create scene content */ + + CHECK(ssol_shape_create_mesh(dev, &square), RES_OK); + attribs[0].usage = SSOL_POSITION; + attribs[0].get = get_position; + CHECK(ssol_mesh_setup(square, SQUARE_NTRIS__, get_ids, + SQUARE_NVERTS__, attribs, 1, (void*) &SQUARE_DESC__), RES_OK); + + CHECK(ssol_material_create_matte(dev, &m_mtl), RES_OK); + shader.normal = get_shader_normal; + shader.reflectivity = get_shader_reflectivity_2; + CHECK(ssol_matte_setup(m_mtl, &shader), RES_OK); + + CHECK(ssol_object_create(dev, &m_object), RES_OK); + CHECK(ssol_object_add_shaded_shape(m_object, square, m_mtl, m_mtl), RES_OK); + CHECK(ssol_object_instantiate(m_object, &geom1), RES_OK); + CHECK(ssol_instance_set_receiver(geom1, SSOL_FRONT, 0), RES_OK); + d3_splat(transform + 9, 0); + transform[9] = -10; + CHECK(ssol_instance_set_transform(geom1, transform), RES_OK); + CHECK(ssol_scene_attach_instance(scene, geom1), RES_OK); + +#define N1__ 10000 +#define GET_MC_RCV ssol_estimator_get_mc_receiver + CHECK(ssol_solve(scene, rng, N1__, 0, NULL, &estimator), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); + CHECK(count, N1__); + CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); + CHECK(count, 0); +#define DNI_S (DNI * X_SZ * Y_SZ) + CHECK(ssol_estimator_get_mc_global(estimator, &mc_global1), RES_OK); + CHECK(GET_MC_RCV(estimator, geom1, SSOL_FRONT, &mc_rcv1), RES_OK); + + PRINT_GLOBAL(mc_global1); + PRINT_RCV(mc_rcv1); + CHECK(mc_global1.cos_factor.E, 1); + CHECK(mc_global1.cos_factor.SE, 0); + CHECK(mc_global1.absorbed_by_receivers.E, DNI_S); + CHECK(mc_global1.absorbed_by_receivers.SE, 0); + + CHECK(ssol_shape_create_punched_surface(dev, &parabol), RES_OK); + carving.get = get_polygon_vertices; + carving.operation = SSOL_AND; + carving.nb_vertices = POLY_NVERTS__; + carving.context = &POLY_EDGES__; + quadric.type = SSOL_QUADRIC_PARABOL; + quadric.data.parabol.focal = 10; + punched.nb_carvings = 1; + punched.quadric = &quadric; + punched.carvings = &carving; + CHECK(ssol_punched_surface_setup(parabol, &punched), RES_OK); + + CHECK(ssol_object_create(dev, &q_object), RES_OK); + CHECK(ssol_object_add_shaded_shape(q_object, parabol, m_mtl, m_mtl), RES_OK); + CHECK(ssol_object_instantiate(q_object, &geom2), RES_OK); + CHECK(ssol_instance_set_receiver(geom2, SSOL_FRONT, 0), RES_OK); + d3_splat(transform + 9, 0); + transform[9] = +10; + CHECK(ssol_instance_set_transform(geom2, transform), RES_OK); + CHECK(ssol_scene_attach_instance(scene, geom2), RES_OK); + + CHECK(ssol_scene_detach_instance(scene, geom1), RES_OK); + CHECK(ssol_estimator_ref_put(estimator), RES_OK); + +#define N2__ 100000 + CHECK(ssol_solve(scene, rng, N2__, 0, NULL, &estimator), RES_OK); + CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); + CHECK(count, N2__); + CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); + CHECK(count, 0); + CHECK(ssol_estimator_get_mc_global(estimator, &mc_global2), RES_OK); + CHECK(GET_MC_RCV(estimator, geom2, SSOL_FRONT, &mc_rcv2), RES_OK); + + PRINT_GLOBAL(mc_global2); + PRINT_RCV(mc_rcv2); + CHECK(eq_eps(mc_global2.absorbed_by_receivers.E, DNI_S, 3 * mc_global2.absorbed_by_receivers.SE), 1); + + /* Free data */ + CHECK(ssol_instance_ref_put(geom1), RES_OK); + CHECK(ssol_instance_ref_put(geom2), RES_OK); + CHECK(ssol_object_ref_put(m_object), RES_OK); + CHECK(ssol_object_ref_put(q_object), RES_OK); + CHECK(ssol_shape_ref_put(square), RES_OK); + CHECK(ssol_shape_ref_put(parabol), RES_OK); + CHECK(ssol_material_ref_put(m_mtl), RES_OK); + CHECK(ssol_estimator_ref_put(estimator), RES_OK); + CHECK(ssol_device_ref_put(dev), RES_OK); + CHECK(ssol_scene_ref_put(scene), RES_OK); + CHECK(ssp_rng_ref_put(rng), RES_OK); + CHECK(ssol_spectrum_ref_put(spectrum), RES_OK); + CHECK(ssol_sun_ref_put(sun), RES_OK); + + check_memory_allocator(&allocator); + mem_shutdown_proxy_allocator(&allocator); + CHECK(mem_allocated_size(), 0); + + return 0; +} diff --git a/src/test_ssol_solver2.c b/src/test_ssol_solver2.c @@ -177,7 +177,7 @@ main(int argc, char** argv) #define N__ 10000 #define GET_MC_RCV ssol_estimator_get_mc_receiver - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); #define COS cos(PI / 4) diff --git a/src/test_ssol_solver2b.c b/src/test_ssol_solver2b.c @@ -180,7 +180,7 @@ main(int argc, char** argv) CHECK(ssol_scene_attach_instance(scene, target), RES_OK); #define N__ 50000 - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); #define COS cos(PI / 4) diff --git a/src/test_ssol_solver3.c b/src/test_ssol_solver3.c @@ -135,7 +135,7 @@ main(int argc, char** argv) CHECK(ssol_scene_attach_instance(scene, target), RES_OK); #define N__ 20000 - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); #define COS cos(PI / 4) diff --git a/src/test_ssol_solver4.c b/src/test_ssol_solver4.c @@ -144,7 +144,7 @@ main(int argc, char** argv) #define N__ 100000 #define GET_MC_RCV ssol_estimator_get_mc_receiver - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); #define COS cos(0) diff --git a/src/test_ssol_solver5.c b/src/test_ssol_solver5.c @@ -135,7 +135,7 @@ main(int argc, char** argv) CHECK(ssol_scene_attach_instance(scene, target), RES_OK); #define N__ 10000 - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); #define COS cos(0) diff --git a/src/test_ssol_solver6.c b/src/test_ssol_solver6.c @@ -176,7 +176,7 @@ main(int argc, char** argv) #define N__ 10000 #define GET_MC_RCV ssol_estimator_get_mc_receiver - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_mc_global(estimator, &mc_global), RES_OK); PRINT_GLOBAL(mc_global); CHECK(eq_eps(mc_global.shadowed.E, 100000, 2 * 100000/sqrt(N__)), 1); diff --git a/src/test_ssol_solver7.c b/src/test_ssol_solver7.c @@ -191,7 +191,7 @@ main(int argc, char** argv) #define TOTAL (HELIOSTAT_SZ * HELIOSTAT_SZ * DNI_cos) #define GET_MC_RCV ssol_estimator_get_mc_receiver - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, 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); diff --git a/src/test_ssol_solver8.c b/src/test_ssol_solver8.c @@ -111,7 +111,7 @@ main(int argc, char** argv) carving.nb_vertices = POLY_NVERTS__; carving.context = &POLY_EDGES__; quadric.type = SSOL_QUADRIC_PARABOLIC_CYLINDER; - quadric.data.parabol.focal = 1; + quadric.data.parabolic_cylinder.focal = 1; punched.nb_carvings = 1; punched.quadric = &quadric; punched.carvings = &carving; @@ -139,7 +139,7 @@ main(int argc, char** argv) #define N__ 100000 #define GET_MC_RCV ssol_estimator_get_mc_receiver - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); diff --git a/src/test_ssol_solver9.c b/src/test_ssol_solver9.c @@ -133,7 +133,7 @@ main(int argc, char** argv) #define N__ 100000 #define GET_MC_RCV ssol_estimator_get_mc_receiver - CHECK(ssol_solve(scene, rng, N__, NULL, &estimator), RES_OK); + CHECK(ssol_solve(scene, rng, N__, 0, NULL, &estimator), RES_OK); CHECK(ssol_estimator_get_realisation_count(estimator, &count), RES_OK); CHECK(count, N__); CHECK(ssol_estimator_get_failed_count(estimator, &count), RES_OK); diff --git a/src/test_ssol_sun.c b/src/test_ssol_sun.c @@ -76,10 +76,10 @@ main(int argc, char** argv) CHECK(ssol_sun_get_dni(sun, &dni), RES_OK); CHECK(dni, 1000); - CHECK(ssol_sun_set_pillbox_aperture(NULL, 0.1), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, -0.1), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, 999), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, 0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(NULL, 0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, -0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 999), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 0.1), RES_BAD_ARG); CHECK(ssol_sun_set_buie_param(NULL, 0.1), RES_BAD_ARG); CHECK(ssol_sun_set_buie_param(sun, -0.1), RES_BAD_ARG); @@ -125,11 +125,11 @@ main(int argc, char** argv) CHECK(ssol_sun_get_dni(sun, &dni), RES_OK); CHECK(dni, 1000); - CHECK(ssol_sun_set_pillbox_aperture(NULL, 0.1), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, -0.1), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, 999), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, 0.1), RES_OK); - CHECK(ssol_sun_set_pillbox_aperture(sun, 0.1), RES_OK); + CHECK(ssol_sun_pillbox_set_theta_max(NULL, 0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, -0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 999), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 0.1), RES_OK); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 0.1), RES_OK); CHECK(ssol_sun_set_buie_param(NULL, 0.1), RES_BAD_ARG); CHECK(ssol_sun_set_buie_param(sun, -0.1), RES_BAD_ARG); @@ -174,10 +174,10 @@ main(int argc, char** argv) CHECK(ssol_sun_get_dni(sun, &dni), RES_OK); CHECK(dni, 1000); - CHECK(ssol_sun_set_pillbox_aperture(NULL, 0.1), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, -0.1), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, 999), RES_BAD_ARG); - CHECK(ssol_sun_set_pillbox_aperture(sun, 0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(NULL, 0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, -0.1), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 999), RES_BAD_ARG); + CHECK(ssol_sun_pillbox_set_theta_max(sun, 0.1), RES_BAD_ARG); CHECK(ssol_sun_set_buie_param(NULL, 0.1), RES_BAD_ARG); CHECK(ssol_sun_set_buie_param(sun, -0.1), RES_BAD_ARG); diff --git a/src/test_ssol_utils.h b/src/test_ssol_utils.h @@ -36,7 +36,7 @@ check_memory_allocator(struct mem_allocator* allocator) 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); \ + (Mc).extinguished_by_atmosphere.E, (Mc).extinguished_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); \