commit ad5ab69e0e416f9f6fd52b2e94a9c62452c27a6d
parent 35309912f04233443f9dbe897286cb22fc9c4df2
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Wed, 17 Apr 2024 11:39:09 +0200
Add a diffuse radiance to the external source
It simulates the radiance of the source that is scattered at least once
in the environment. This is a new functor on the spherical source that
depends on the time and direction along which this diffuse radiance is
interrogated. Once defined, the external diffuse flux takes into
account not only the radiance emitted by the external source and
reflected on surfaces but also, when the path goes to infinity, the
radiance of the source scattered in the environment which is returned by
this functor.
In fact, it's difficult to make it intelligible. And you have to take
care with the comments, which should help the reader to understand
what's going on. To sum up, external flux is now made up of 3 elements:
- a direct component
- a diffuse component coming from the source and reflected on surfaces
- a diffuse component originating from the source and scattered by the
environment, which is the novelty of this commit.
Green's function has been updated accordingly, since the external flux
is no longer linked solely to the source's power: one of its components
is linked to the source's diffuse radiance. So there is no longer one
constant term per Green's path, but a list of 2 terms depending on time
and direction. One term relates to the power of the source, the other to
the diffuse radiance of the source. The API for querying the Green
function has been updated to take account of these changes. In fact,
querying external flux now resembles querying other additive terms such
as power density or imposed flux, i.e. it relies on iterative functions.
Diffstat:
8 files changed, 416 insertions(+), 113 deletions(-)
diff --git a/src/sdis.h b/src/sdis.h
@@ -171,24 +171,39 @@ static const struct sdis_info SDIS_INFO_NULL = SDIS_INFO_NULL__;
/* Type of functor used to retrieve the source's position relative to time */
typedef void
(*sdis_get_position_T)
- (const double time,
+ (const double time, /* [s] */
double pos[3],
struct sdis_data* data);
/* Type of functor used to retrieve the source's power relative to time */
typedef double
(*sdis_get_power_T)
- (const double time,
+ (const double time, /* [s] */
+ struct sdis_data* data);
+
+/* Type of functor used to retrieve the diffuse part of the external radiance */
+typedef double /* [W/perpendicular m^2/sr] */
+(*sdis_get_diffuse_radiance_T)
+ (const double time, /* [s] */
+ const double dir[3],
struct sdis_data* data);
/* Input arguments of the sdis_spherical_source_create function */
struct sdis_spherical_source_create_args {
- sdis_get_position_T position; /* [m] */
+ sdis_get_position_T position; /* [m/fp_to_meter] */
sdis_get_power_T power; /* Total power [W] */
+
+ /* Describes the diffuse part of the source's radiance, i.e. the radiance
+ * emitted by the source and scattered at least once in the environment. This
+ * parameter is actually used to approximate a semi-transparent medium. Its
+ * value can be NULL, meaning that the source has not been scattered by the
+ * environment, or, to put it another way, that the source is in a vacuum. */
+ sdis_get_diffuse_radiance_T diffuse_radiance; /* [W/m^2/sr] */
+
struct sdis_data* data; /* Data sent to the position functor */
double radius; /* [m] */
};
-#define SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__ {NULL, NULL, 0, 0}
+#define SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__ {NULL, NULL, NULL, 0, 0}
static const struct sdis_spherical_source_create_args
SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL =
SDIS_SPHERICAL_SOURCE_CREATE_ARGS_NULL__;
@@ -441,27 +456,49 @@ struct sdis_green_path_end {
static const struct sdis_green_path_end SDIS_GREEN_PATH_END_NULL =
SDIS_GREEN_PATH_END_NULL__;
-/* Functor used to process the paths registered against the green function */
+struct sdis_green_external_flux_terms {
+ /* Term relative to source power [K/W] */
+ double term_wrt_power;
+
+ /* Term relative to diffuse source radiance [K/W/m^2/sr] */
+ double term_wrt_diffuse_radiance;
+
+ double time; /* [s] */
+ double dir[3]; /* Direction on which term_wrt_diffuse_radiance depends */
+};
+#define SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL__ {0,0,0,{0,0,0}}
+static const struct sdis_green_external_flux_terms
+SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL = SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL__;
+
+/* Function profile used to process the paths stored in the green function */
typedef res_T
(*sdis_process_green_path_T)
(struct sdis_green_path* path,
void* context);
-/* Functor used to process the power factor registered along a green path for a
- * given medium */
+/* Function profile used to process power factors registered along a green path
+ * for a given medium */
typedef res_T
(*sdis_process_medium_power_term_T)
(struct sdis_medium* medium,
- const double power_term,
+ const double power_term, /* [K/W] */
void* context);
-/* Functor used to process the flux factor registered along a green path for a
- * given interface side */
+/* Function profile used to process flux factors recorded along a green path for
+ * a given interface side */
typedef res_T
(*sdis_process_interface_flux_term_T)
(struct sdis_interface* interf,
const enum sdis_side side,
- const double flux_term,
+ const double flux_term, /* [K/W/m^2] */
+ void* context);
+
+/* Function profile used to process external flux factors recorded along a green
+ * path */
+typedef res_T
+(*sdis_process_external_flux_terms_T)
+ (struct sdis_source* source,
+ const struct sdis_green_external_flux_terms* terms,
void* context);
/*******************************************************************************
@@ -1118,11 +1155,18 @@ SDIS_API res_T
sdis_source_ref_put
(struct sdis_source* source);
-SDIS_API double
+SDIS_API double /* [W] */
sdis_source_get_power
(struct sdis_source* source,
const double time); /* [s] */
+/* Return the source radiance that is diffused in the environment */
+SDIS_API double /* [W/m^2/sr*] */
+sdis_source_get_diffuse_radiance
+ (struct sdis_source* source,
+ const double time, /* [s] */
+ const double dir[3]);
+
SDIS_API unsigned
sdis_source_get_id
(const struct sdis_source* source);
@@ -1486,6 +1530,11 @@ sdis_green_function_get_flux_terms_count
(const struct sdis_green_path* path,
size_t* nterms);
+SDIS_API res_T
+sdis_green_function_get_external_flux_terms_count
+ (const struct sdis_green_path* path,
+ size_t* nterms);
+
/* Iterate over all "power terms" associated to the path. Multiply each term
* by the power of their associated medium, that is assumed to be constant in
* time and space, gives the medium power registered along the path. */
@@ -1504,13 +1553,12 @@ sdis_green_path_for_each_flux_term
sdis_process_interface_flux_term_T func,
void* context);
-/* Return the external flux term, i.e. the relative net flux along the path from
- * the external source. Multiply it by the power of the source to obtain its
- * contribution to the path. */
+/* Iterate over all external flux terms associated to the path */
SDIS_API res_T
-sdis_green_path_get_external_flux_term
+sdis_green_path_for_each_external_flux_terms
(struct sdis_green_path* path,
- double* external_flux_term); /* [W/m^2] */
+ sdis_process_external_flux_terms_T func,
+ void* context);
/*******************************************************************************
* Heat path API
diff --git a/src/sdis_green.c b/src/sdis_green.c
@@ -26,6 +26,7 @@
#include <star/ssp.h>
+#include <rsys/cstr.h>
#include <rsys/dynamic_array.h>
#include <rsys/hash_table.h>
#include <rsys/mem_allocator.h>
@@ -83,9 +84,24 @@ flux_term_init(struct mem_allocator* allocator, struct flux_term* term)
#define DARRAY_FUNCTOR_INIT flux_term_init
#include <rsys/dynamic_array.h>
+static INLINE void
+extflux_terms_init
+ (struct mem_allocator* allocator,
+ struct sdis_green_external_flux_terms* terms)
+{
+ ASSERT(terms); (void)allocator;
+ *terms = SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL;
+}
+
+/* Generate the dynamic array of the external flux terms */
+#define DARRAY_NAME extflux_terms
+#define DARRAY_DATA struct sdis_green_external_flux_terms
+#define DARRAY_FUNCTOR_INIT extflux_terms_init
+#include <rsys/dynamic_array.h>
+
struct green_path {
double elapsed_time;
- double external_flux_term; /* [W/m^2] */
+ struct darray_extflux_terms extflux_terms; /* List of external flux terms */
struct darray_flux_term flux_terms; /* List of flux terms */
struct darray_power_term power_terms; /* List of volumic power terms */
union {
@@ -106,10 +122,10 @@ static INLINE void
green_path_init(struct mem_allocator* allocator, struct green_path* path)
{
ASSERT(path);
+ darray_extflux_terms_init(allocator, &path->extflux_terms);
darray_flux_term_init(allocator, &path->flux_terms);
darray_power_term_init(allocator, &path->power_terms);
path->elapsed_time = -INF;
- path->external_flux_term = 0;
path->limit.vertex = SDIS_RWALK_VERTEX_NULL;
path->limit.fragment = SDIS_INTERFACE_FRAGMENT_NULL;
path->limit.ray = SDIS_RADIATIVE_RAY_NULL;
@@ -133,12 +149,13 @@ green_path_copy(struct green_path* dst, const struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
- dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
dst->ilast_medium = src->ilast_medium;
dst->ilast_interf = src->ilast_interf;
+ res = darray_extflux_terms_copy(&dst->extflux_terms, &src->extflux_terms);
+ if(res != RES_OK) return res;
res = darray_flux_term_copy(&dst->flux_terms, &src->flux_terms);
if(res != RES_OK) return res;
res = darray_power_term_copy(&dst->power_terms, &src->power_terms);
@@ -152,12 +169,14 @@ green_path_copy_and_clear(struct green_path* dst, struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
- dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
dst->ilast_medium = src->ilast_medium;
dst->ilast_interf = src->ilast_interf;
+ res = darray_extflux_terms_copy_and_clear
+ (&dst->extflux_terms, &src->extflux_terms);
+ if(res != RES_OK) return res;
res = darray_flux_term_copy_and_clear(&dst->flux_terms, &src->flux_terms);
if(res != RES_OK) return res;
res = darray_power_term_copy_and_clear(&dst->power_terms, &src->power_terms);
@@ -172,12 +191,14 @@ green_path_copy_and_release(struct green_path* dst, struct green_path* src)
res_T res = RES_OK;
ASSERT(dst && src);
dst->elapsed_time = src->elapsed_time;
- dst->external_flux_term = src->external_flux_term;
dst->limit = src->limit;
dst->limit_id = src->limit_id;
dst->end_type = src->end_type;
dst->ilast_medium = src->ilast_medium;
dst->ilast_interf = src->ilast_interf;
+ res = darray_extflux_terms_copy_and_release
+ (&dst->extflux_terms, &src->extflux_terms);
+ if(res != RES_OK) return res;
res = darray_flux_term_copy_and_release(&dst->flux_terms, &src->flux_terms);
if(res != RES_OK) return res;
res = darray_power_term_copy_and_release(&dst->power_terms, &src->power_terms);
@@ -201,7 +222,11 @@ green_path_write(const struct green_path* path, FILE* stream)
/* Write elapsed time */
WRITE(&path->elapsed_time, 1);
- WRITE(&path->external_flux_term, 1);
+
+ /* Write the list of external flux terms */
+ sz = darray_extflux_terms_size_get(&path->extflux_terms);
+ WRITE(&sz, 1);
+ WRITE(darray_extflux_terms_cdata_get(&path->extflux_terms), sz);
/* Write the list of flux terms */
sz = darray_flux_term_size_get(&path->flux_terms);
@@ -252,7 +277,12 @@ green_path_read(struct green_path* path, FILE* stream)
/* Read elapsed time */
READ(&path->elapsed_time, 1);
- READ(&path->external_flux_term, 1);
+
+ /* Read the list of external flux terms */
+ READ(&sz, 1);
+ res = darray_extflux_terms_resize(&path->extflux_terms, sz);
+ if(res != RES_OK) goto error;
+ READ(darray_extflux_terms_data_get(&path->extflux_terms), sz);
/* Read the list of flux terms */
READ(&sz, 1);
@@ -402,14 +432,105 @@ green_function_fetch_interf
return *pinterf;
}
+static double /* [K] */
+green_path_power_contribution
+ (struct sdis_green_function* green,
+ const struct green_path* path)
+{
+ double temperature = 0; /* [K] */
+ size_t i=0, n=0;
+
+ ASSERT(green && path);
+
+ n = darray_power_term_size_get(&path->power_terms);
+ FOR_EACH(i, 0, n) {
+ struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
+ const struct power_term* power_term = NULL;
+ const struct sdis_medium* medium = NULL;
+ double power = 0; /* [W] */
+
+ power_term = darray_power_term_cdata_get(&path->power_terms) + i;
+ medium = green_function_fetch_medium(green, power_term->id);
+
+ /* Dummy argument used only to satisfy the function profile used to recover
+ * power. Its position is unused, since power is assumed to be constant in
+ * space, and its time is set to infinity, since the green function is
+ * assumed to be evaluated at steady state */
+ vtx.time = INF;
+ power = solid_get_volumic_power(medium, &vtx);
+
+ temperature += power_term->term * power; /* [K] */
+ }
+ return temperature; /* [K] */
+}
+
+static double /* [K] */
+green_path_flux_contribution
+ (struct sdis_green_function* green,
+ const struct green_path* path)
+{
+ double temperature = 0;
+ size_t i=0, n=0;
+
+ ASSERT(green && path);
+
+ n = darray_flux_term_size_get(&path->flux_terms);
+ FOR_EACH(i, 0, n) {
+ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
+ const struct flux_term* flux_term = NULL;
+ const struct sdis_interface* interf = NULL;
+ double flux = 0; /* [W/m^2] */
+
+ flux_term = darray_flux_term_cdata_get(&path->flux_terms) + i;
+ interf = green_function_fetch_interf(green, flux_term->id);
+
+ /* Interface fragment. Its position is unused, since flux is assumed to be
+ * constant in space, and its time is set to infinity, since the green
+ * function is assumed to be evaluated at steady state */
+ frag.time = INF;
+ frag.side = flux_term->side;
+ flux = interface_side_get_flux(interf, &frag);
+
+ temperature += flux_term->term * flux; /* [K] */
+ }
+ return temperature; /* [K] */
+}
+
+static double /* [K] */
+green_path_external_flux_contribution
+ (struct sdis_green_function* green,
+ const struct green_path* path)
+{
+ const struct sdis_source* extsrc = NULL;
+ double value = 0;
+ size_t i=0, n=0;
+
+ ASSERT(green && path);
+
+ if((extsrc = green->scn->source) == NULL) return 0;
+
+ n = darray_extflux_terms_size_get(&path->extflux_terms);
+ FOR_EACH(i, 0, n) {
+ const struct sdis_green_external_flux_terms* extflux = NULL;
+ double power = 0; /* [W] */
+ double diffrad = 0; /* [W/m^2/sr] */
+
+ extflux = darray_extflux_terms_cdata_get(&path->extflux_terms)+i;
+ power = source_get_power(extsrc, extflux->time);
+ diffrad = source_get_diffuse_radiance(extsrc, extflux->time, extflux->dir);
+
+ value += extflux->term_wrt_power * power; /* [K] */
+ value += extflux->term_wrt_diffuse_radiance * diffrad; /* [K] */
+ }
+ return value; /* [K] */
+}
+
static res_T
green_function_solve_path
(struct sdis_green_function* green,
const size_t ipath,
double* weight)
{
- const struct power_term* power_terms = NULL;
- const struct flux_term* flux_terms = NULL;
const struct green_path* path = NULL;
const struct sdis_medium* medium = NULL;
const struct sdis_interface* interf = NULL;
@@ -421,7 +542,6 @@ green_function_solve_path
double flux;
double external_flux;
double end_temperature;
- size_t i, n;
res_T res = RES_OK;
ASSERT(green && ipath < darray_green_path_size_get(&green->paths) && weight);
@@ -431,35 +551,9 @@ green_function_solve_path
goto error;
}
- /* Compute medium powers */
- power = 0;
- n = darray_power_term_size_get(&path->power_terms);
- power_terms = darray_power_term_cdata_get(&path->power_terms);
- FOR_EACH(i, 0, n) {
- vtx.time = INF;
- medium = green_function_fetch_medium(green, power_terms[i].id);
- power += power_terms[i].term * solid_get_volumic_power(medium, &vtx);
- }
-
- /* Compute interface fluxes */
- flux = 0;
- n = darray_flux_term_size_get(&path->flux_terms);
- flux_terms = darray_flux_term_cdata_get(&path->flux_terms);
- FOR_EACH(i, 0, n) {
- frag.time = INF;
- frag.side = flux_terms[i].side;
- interf = green_function_fetch_interf(green, flux_terms[i].id);
- flux += flux_terms[i].term * interface_side_get_flux(interf, &frag);
- }
-
- /* Compute external flux */
- external_flux = 0;
- if(green->scn->source) {
- /* NOTE: The power of the source is assumed to be constant over time and is
- * therefore recovered at steady state */
- external_flux =
- path->external_flux_term * source_get_power(green->scn->source, INF);
- }
+ power = green_path_power_contribution(green, path);
+ flux = green_path_flux_contribution(green, path);
+ external_flux = green_path_external_flux_contribution(green, path);
/* Compute path's end temperature */
switch(path->end_type) {
@@ -1279,6 +1373,33 @@ error:
}
res_T
+sdis_green_function_get_external_flux_terms_count
+ (const struct sdis_green_path* path_handle,
+ size_t* nterms)
+{
+ const struct green_path* path = NULL;
+ struct sdis_green_function* green = NULL;
+ res_T res = RES_OK;
+
+ if(!path_handle || !nterms) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ green = path_handle->green__; (void)green;
+ ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths));
+
+ path = darray_green_path_cdata_get(&green->paths) + path_handle->id__;
+
+ *nterms = darray_extflux_terms_size_get(&path->extflux_terms);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
sdis_green_path_for_each_power_term
(struct sdis_green_path* path_handle,
sdis_process_medium_power_term_T func,
@@ -1354,25 +1475,41 @@ error:
}
res_T
-sdis_green_path_get_external_flux_term
+sdis_green_path_for_each_external_flux_terms
(struct sdis_green_path* path_handle,
- double* external_flux_term)
+ sdis_process_external_flux_terms_T func,
+ void* context)
{
const struct green_path* path = NULL;
struct sdis_green_function* green = NULL;
+ size_t i, n;
res_T res = RES_OK;
- if(!path_handle || !external_flux_term) {
+ if(!path_handle || !func) {
res = RES_BAD_ARG;
goto error;
}
- green = path_handle->green__; (void)green;
+ green = path_handle->green__;
ASSERT(path_handle->id__ < darray_green_path_size_get(&green->paths));
path = darray_green_path_cdata_get(&green->paths) + path_handle->id__;
- *external_flux_term = path->external_flux_term;
+ n = darray_extflux_terms_size_get(&path->extflux_terms);
+ if(n && !green->scn->source) {
+ /* In can't have external flux terms without an external source */
+ log_err(green->scn->dev, "%s: the external source is missing\n", FUNC_NAME);
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ FOR_EACH(i, 0, n) {
+ const struct sdis_green_external_flux_terms* terms = NULL;
+ terms = darray_extflux_terms_cdata_get(&path->extflux_terms) + i;
+
+ res = func(green->scn->source, terms, context);
+ if(res != RES_OK) goto error;
+ }
exit:
return res;
@@ -1743,11 +1880,23 @@ error:
}
res_T
-green_path_add_external_flux_term
+green_path_add_external_flux_terms
(struct green_path_handle* handle,
- const double val) /* [W/m^2/sr] */
+ const struct sdis_green_external_flux_terms* terms)
{
- ASSERT(handle);
- handle->path->external_flux_term += val;
- return RES_OK;
+ res_T res = RES_OK;
+ ASSERT(handle && terms);
+
+ res = darray_extflux_terms_push_back(&handle->path->extflux_terms, terms);
+ if(res != RES_OK) {
+ log_err(handle->green->scn->dev,
+ "%s: cannot store external flux terms -- %s\n",
+ FUNC_NAME, res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
}
diff --git a/src/sdis_green.h b/src/sdis_green.h
@@ -108,9 +108,9 @@ green_path_add_flux_term
const double term);
extern LOCAL_SYM res_T
-green_path_add_external_flux_term
+green_path_add_external_flux_terms
(struct green_path_handle* handle,
- const double term); /* [W/m^2/sr] */
+ const struct sdis_green_external_flux_terms* terms);
#endif /* SDIS_GREEN_H */
diff --git a/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h b/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h
@@ -50,6 +50,19 @@ struct brdf {
#define BRDF_NULL__ {0, 0}
static const struct brdf BRDF_NULL = BRDF_NULL__;
+/* Incident diffuse flux is made up of two components. One corresponds to the
+ * diffuse flux due to the reflection of the source on surfaces. The other is
+ * the diffuse flux due to the source's radiation scattering at least once in
+ * the environment. */
+struct incident_diffuse_flux {
+ double reflected; /* [W/m^2] */
+ double scattered; /* [W/m^2] */
+ double dir[3]; /* Direction along wich the scattered part was retrieved */
+};
+#define INCIDENT_DIFFUSE_FLUX_NULL__ {0, 0, {0,0,0}}
+static const struct incident_diffuse_flux INCIDENT_DIFFUSE_FLUX_NULL =
+ INCIDENT_DIFFUSE_FLUX_NULL__;
+
/* Reflect the V wrt the normal N. By convention V points outward the surface.
* In fact, this function is a double-precision version of the reflect_3d
* function. TODO Clean this "repeat" */
@@ -314,15 +327,14 @@ XD(compute_incident_diffuse_flux)
const double in_N[DIM], /* Surface normal. (Away from the surface) */
const double time,
const struct sXd(hit)* in_hit, /* Current intersection */
- double* out_flux) /* [W/m^2] */
+ struct incident_diffuse_flux* diffuse_flux) /* [W/m^2] */
{
struct sXd(hit) hit = SXD_HIT_NULL;
double pos[3] = {0}; /* In 3D for ray tracing ray to the source */
double dir[3] = {0}; /* Incident direction (toward the surface). Always 3D.*/
double N[3] = {0}; /* Surface normal. Always 3D */
- double incident_diffuse_flux = 0; /* [W/m^2] */
res_T res = RES_OK;
- ASSERT(in_pos && in_N && in_hit);
+ ASSERT(in_pos && in_N && in_hit && diffuse_flux);
/* Local copy of input argument */
dX(set)(pos, in_pos);
@@ -332,6 +344,8 @@ XD(compute_incident_diffuse_flux)
/* Sample a diffusive direction in 3D */
ssp_ran_hemisphere_cos(rng, N, dir, NULL);
+ *diffuse_flux = INCIDENT_DIFFUSE_FLUX_NULL;
+
for(;;) {
/* External sources */
struct source_sample src_sample = SOURCE_SAMPLE_NULL;
@@ -345,7 +359,7 @@ XD(compute_incident_diffuse_flux)
struct brdf_sample brdf_sample = BRDF_SAMPLE_NULL;
/* Miscellaneous */
- double L = 0; /* incident direct flux to bounce position */
+ double L = 0; /* incident flux to bounce position */
double wi[3] = {0}; /* Incident direction (outward the surface). Always 3D */
double vec[DIM] = {0}; /* Temporary variable */
@@ -353,7 +367,21 @@ XD(compute_incident_diffuse_flux)
/* Find the following surface along the direction of propagation */
XD(trace_ray)(scn, pos, dir, INF, &hit, &hit);
- if(SXD_HIT_NONE(&hit)) break; /* No surface */
+ if(SXD_HIT_NONE(&hit)) {
+ /* No surface. Handle the radiance emitted by the source and scattered at
+ * least once in the environment. Note that the value returned is not the
+ * actual scattered component of the incident diffuse flux: it relates
+ * to the radiance of the source scattered along the input dir at the
+ * given instant. It must therefore be multiplied by this radiance to
+ * obtain its real contribution. This trick makes it possible to manage
+ * the external flux in the green function. */
+ const double Ld = source_get_diffuse_radiance(scn->source, time, dir);
+ diffuse_flux->scattered = Ld * PI; /* [W/m^2] */
+ diffuse_flux->dir[0] = dir[0];
+ diffuse_flux->dir[1] = dir[1];
+ diffuse_flux->dir[2] = dir[2];
+ break;
+ }
/* Retrieve the current position and normal */
dX(add)(pos, pos, dX(muld)(vec, dir, hit.distance));
@@ -385,7 +413,7 @@ XD(compute_incident_diffuse_flux)
if(!SOURCE_SAMPLE_NONE(&src_sample)) {
const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit);
- L = Ld; /* [W/m^2] */
+ L = Ld; /* [W/m^2/sr] */
}
/* Calculate the direct contribution of the rebound is diffuse */
@@ -401,21 +429,19 @@ XD(compute_incident_diffuse_flux)
/* The source is behind the surface */
if(cos_theta <= 0) {
- L = 0; /* [W/m^2] */
+ L = 0; /* [W/m^2/sr] */
/* The source is above the surface */
} else {
const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit);
- L = Ld * cos_theta / (PI * src_sample.pdf); /* [W/m^2] */
+ L = Ld * cos_theta / (PI * src_sample.pdf); /* [W/m^2/sr] */
}
}
- incident_diffuse_flux += L;
+ diffuse_flux->reflected += L; /* [W/m^2/sr] */
}
-
- incident_diffuse_flux *= PI;
+ diffuse_flux->reflected *= PI; /* [W/m^2] */
exit:
- *out_flux = incident_diffuse_flux;
return res;
error:
goto exit;
@@ -431,15 +457,19 @@ XD(handle_external_net_flux)
const struct XD(handle_external_net_flux_args)* args,
struct XD(temperature)* T)
{
+ /* Terms to be registered in the green function */
+ struct sdis_green_external_flux_terms green =
+ SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL;
+
/* Sampling external sources */
struct source_sample src_sample = SOURCE_SAMPLE_NULL;
/* External flux */
+ struct incident_diffuse_flux incident_flux_diffuse = INCIDENT_DIFFUSE_FLUX_NULL;
double incident_flux = 0; /* [W/m^2] */
- double incident_flux_diffuse = 0; /* [W/m^2] */
double incident_flux_direct = 0; /* [W/m^2] */
double net_flux = 0; /* [W/m^2] */
- double external_flux_term = 0; /* [W/m^2] */
+ double net_flux_sc = 0; /* [W/m^2] */
/* Sampled path */
double N[3] = {0}; /* Normal. Always in 3D */
@@ -448,6 +478,7 @@ XD(handle_external_net_flux)
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
/* Miscellaneous */
+ double sum_h = 0;
double emissivity = 0; /* Emissivity */
double Ld = 0; /* Incident radiance [W/m^2/sr] */
double cos_theta = 0;
@@ -495,27 +526,50 @@ XD(handle_external_net_flux)
(scn, rng, frag.P, N, frag.time, args->hit, &incident_flux_diffuse);
if(res != RES_OK) goto error;
- /* Calculate the global incident flux */
- incident_flux = incident_flux_direct + incident_flux_diffuse; /* [W/m^2] */
-
- /* Calculate the net flux */
+ /* Calculate the incident flux without the part scattered by the environment.
+ * The latter depends on the source's diffuse radiance, not on its power. On
+ * the other hand, both the direct incident flux and the diffuse incident flux
+ * reflected by surfaces are linear with respect to the source power. This
+ * term can therefore be recorded in the green function in relation to this
+ * power, whereas the incident diffused flux coming from the scattered source
+ * radiance depends on the diffuse radiance of the source */
+ incident_flux = /* [W/m^2] */
+ incident_flux_direct + incident_flux_diffuse.reflected;
+
+ /* Calculate the net flux [W/m^2] */
src_id = sdis_source_get_id(scn->source);
emissivity = interface_side_get_emissivity(args->interf, src_id, &frag);
res = interface_side_check_emissivity(scn->dev, emissivity, frag.P, frag.time);
if(res != RES_OK) goto error;
net_flux = incident_flux * emissivity; /* [W/m^2] */
- /* Until now, the net flux was calculated in relation to the source power.
- * What is calculated is the external flux term of the green function. This
- * must be multiplied by the source power to obtain the actual external flux*/
- external_flux_term = net_flux / (args->h_radi + args->h_conv + args->h_cond);
-
- /* Update the Monte Carlo weight */
- T->value += external_flux_term * source_get_power(scn->source, frag.time);
+ /* Calculate the net flux from the radiance source scattered at least once by
+ * the medium */
+ net_flux_sc = incident_flux_diffuse.scattered * emissivity; /* [W/m^2] */
+
+ /* Until now, net flux has been calculated on the basis of source power and
+ * source diffuse radiance. What is actually calculated are the external flux
+ * terms of the green function. These must be multiplied by the source power
+ * and the source diffuse radiance, then added together to give the actual
+ * external flux */
+ sum_h = (args->h_radi + args->h_conv + args->h_cond);
+ green.term_wrt_power = net_flux / sum_h; /* [K/W] */
+ green.term_wrt_diffuse_radiance = net_flux_sc / sum_h; /* [K/W/m^2/sr] */
+ green.time = frag.time; /* [s] */
+ green.dir[0] = incident_flux_diffuse.dir[0];
+ green.dir[1] = incident_flux_diffuse.dir[1];
+ green.dir[2] = incident_flux_diffuse.dir[2];
+
+ T->value += green.term_wrt_power * source_get_power(scn->source, green.time);
+ if(green.term_wrt_diffuse_radiance) {
+ T->value +=
+ green.term_wrt_diffuse_radiance
+ * source_get_diffuse_radiance(scn->source, green.time, green.dir);
+ }
- /* Register the external net flux term */
+ /* Register the external net flux terms */
if(args->green_path) {
- res = green_path_add_external_flux_term(args->green_path, external_flux_term);
+ res = green_path_add_external_flux_terms(args->green_path, &green);
if(res != RES_OK) goto error;
}
diff --git a/src/sdis_source.c b/src/sdis_source.c
@@ -133,10 +133,18 @@ sdis_source_ref_put(struct sdis_source* src)
double
sdis_source_get_power(struct sdis_source* src, const double time /* [s] */)
{
- ASSERT(src);
return source_get_power(src, time);
}
+double
+sdis_source_get_diffuse_radiance
+ (struct sdis_source* src,
+ const double time, /* [s] */
+ const double dir[3])
+{
+ return source_get_diffuse_radiance(src, time, dir);
+}
+
unsigned
sdis_source_get_id(const struct sdis_source* source)
{
@@ -321,6 +329,20 @@ source_get_power(const struct sdis_source* src, const double time /* [s] */)
return src->spherical.power(time, src->spherical.data);
}
+double /* [W/perpendicular m^2/sr] */
+source_get_diffuse_radiance
+ (const struct sdis_source* src,
+ const double time /* [s] */,
+ const double dir[3])
+{
+ ASSERT(src);
+ if(src->spherical.diffuse_radiance == NULL) {
+ return 0;
+ } else {
+ return src->spherical.diffuse_radiance(time, dir, src->spherical.data);
+ }
+}
+
void
source_compute_signature(const struct sdis_source* src, hash256_T hash)
{
diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h
@@ -66,6 +66,12 @@ source_get_power
(const struct sdis_source* source,
const double time); /* [s] */
+extern LOCAL_SYM double /* [W/m^2/sr] */
+source_get_diffuse_radiance
+ (const struct sdis_source* source,
+ const double time, /* [s] */
+ const double dir[3]);
+
extern LOCAL_SYM void
source_compute_signature
(const struct sdis_source* source,
diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c
@@ -191,6 +191,17 @@ source_get_power(const double time, struct sdis_data* data)
return SOURCE_POWER; /* [W] */
}
+static double
+source_get_diffuse_radiance
+ (const double time,
+ const double dir[3],
+ struct sdis_data* data)
+{
+ (void)time, (void)data; /* Avoid the "unusued variable" warning */
+ CHK(d3_is_normalized(dir));
+ return 50;
+}
+
static struct sdis_source*
create_source(struct sdis_device* sdis)
{
@@ -199,6 +210,7 @@ create_source(struct sdis_device* sdis)
args.position = source_get_position;
args.power = source_get_power;
+ args.diffuse_radiance = source_get_diffuse_radiance;
args.data = NULL;
args.radius = 3e-1; /* [m] */
OK(sdis_spherical_source_create(sdis, &args, &source));
diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c
@@ -42,7 +42,7 @@ accum_power_terms(struct sdis_medium* mdm, const double power_term, void* ctx)
struct sdis_solid_shader shader = SDIS_SOLID_SHADER_NULL;
struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL;
struct sdis_data* data = NULL;
- double* power = ctx;
+ double* power = ctx; /* Power contribution [K] */
CHK(mdm && ctx);
CHK(sdis_medium_get_type(mdm) == SDIS_SOLID);
@@ -65,7 +65,7 @@ accum_flux_terms
struct sdis_interface_shader shader = SDIS_INTERFACE_SHADER_NULL;
struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL;
struct sdis_data* data = NULL;
- double* flux = ctx;
+ double* flux = ctx; /* Flux contribution [K] */
double phi;
CHK(interf && ctx);
@@ -84,6 +84,27 @@ accum_flux_terms
}
static res_T
+accum_extflux
+ (struct sdis_source* source,
+ const struct sdis_green_external_flux_terms* terms,
+ void* ctx)
+{
+ double* extflux = ctx; /* External flux contribution [K] */
+ double power = 0; /* [W] */
+ double diffuse_radiance = 0; /* [W/m^2/sr] */
+
+ CHK(source && terms && ctx);
+
+ power = sdis_source_get_power(source, terms->time);
+ diffuse_radiance = sdis_source_get_diffuse_radiance
+ (source, terms->time, terms->dir);
+
+ *extflux += terms->term_wrt_power * power;
+ *extflux += terms->term_wrt_diffuse_radiance * diffuse_radiance;
+ return RES_OK;
+}
+
+static res_T
solve_green_path(struct sdis_green_path* path, void* ctx)
{
struct sdis_green_path_end end = SDIS_GREEN_PATH_END_NULL;
@@ -99,13 +120,12 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
struct sdis_green_function* green = NULL;
struct sdis_scene* scn = NULL;
- struct sdis_source* source = NULL;
struct green_accum* acc = NULL;
struct sdis_data* data = NULL;
enum sdis_medium_type type;
double power = 0;
double flux = 0;
- double external_flux = 0; /* [W/m^2] */
+ double extflux = 0;
double time, temp = 0;
double weight = 0;
CHK(path && ctx);
@@ -130,23 +150,15 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
BA(sdis_green_path_for_each_flux_term(path, NULL, &acc));
OK(sdis_green_path_for_each_flux_term(path, accum_flux_terms, &flux));
+ BA(sdis_green_path_for_each_external_flux_terms(NULL, &accum_extflux, &extflux));
+ BA(sdis_green_path_for_each_external_flux_terms(path, NULL, &extflux));
+ OK(sdis_green_path_for_each_external_flux_terms(path, &accum_extflux, &extflux));
+
BA(sdis_green_path_get_elapsed_time(NULL, NULL));
BA(sdis_green_path_get_elapsed_time(path, NULL));
BA(sdis_green_path_get_elapsed_time(NULL, &time));
OK(sdis_green_path_get_elapsed_time(path, &time));
- BA(sdis_green_path_get_external_flux_term(NULL, &external_flux));
- BA(sdis_green_path_get_external_flux_term(path, NULL));
- OK(sdis_green_path_get_external_flux_term(path, &external_flux));
- OK(sdis_scene_get_source(scn, &source));
- if(source == NULL) {
- CHK(external_flux == 0);
- } else {
- /* NOTE: source power is assumed constant in time and is therefore retrieved
- * at steady state*/
- external_flux *= sdis_source_get_power(source, INF); /* [W] */
- }
-
BA(sdis_green_path_get_end(NULL, NULL));
BA(sdis_green_path_get_end(NULL, &end));
BA(sdis_green_path_get_end(path, NULL));
@@ -184,7 +196,7 @@ solve_green_path(struct sdis_green_path* path, void* ctx)
default: FATAL("Unreachable code.\n"); break;
}
- weight = temp + power + external_flux + flux;
+ weight = temp + power + extflux + flux;
acc->sum += weight;
acc->sum2 += weight*weight;
