solstice-solver

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

commit 1d88319669c99645f1df3031478cd838e7fa307a
parent 8009404285eecf4930527cc627e4e2bc141783e3
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 12 Jul 2016 15:04:13 +0200

Rename the ssol_randst_sun_dir structure and its functions

Diffstat:
Mcmake/CMakeLists.txt | 4++--
Dsrc/ssol_distributions.c | 364-------------------------------------------------------------------------------
Dsrc/ssol_distributions_c.h | 66------------------------------------------------------------------
Asrc/ssol_ranst_sun_dir.c | 365+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Asrc/ssol_ranst_sun_dir.h | 63+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Msrc/ssol_solver.c | 6+++---
Msrc/ssol_solver_c.h | 4++--
7 files changed, 435 insertions(+), 437 deletions(-)

diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt @@ -50,11 +50,11 @@ set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) set(SSOL_FILES_SRC ssol_device.c - ssol_distributions.c ssol_image.c ssol_material.c ssol_object.c ssol_object_instance.c + ssol_ranst_sun_dir.c ssol_scene.c ssol_shape.c ssol_spectrum.c @@ -66,11 +66,11 @@ set(SSOL_FILES_INC_API set(SSOL_FILES_INC ssol_device_c.h - ssol_distributions_c.h ssol_image_c.h ssol_material_c.h ssol_object_c.h ssol_object_instance_c.h + ssol_ranst_sun_dir.h ssol_scene_c.h ssol_shape_c.h ssol_spectrum_c.h diff --git a/src/ssol_distributions.c b/src/ssol_distributions.c @@ -1,364 +0,0 @@ -/* Copyright (C) CNRS 2016 - * - * 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 "ssol_distributions_c.h" - -#include <star/ssp.h> - -#include <rsys/double33.h> -#include <rsys/math.h> -#include <rsys/mem_allocator.h> -#include <rsys/rsys.h> -#include <rsys/ref_count.h> - -/******************************************************************************* - * Distributions types - ******************************************************************************/ -struct ran_buie_state { - double thetaSD; - double deltaThetaCSSD; - double gamma; - double etokTimes1000toGamma; - double alpha; - double hRect1; - double hRect2; - double proba_rect1; - double basis[9]; -}; - -struct ran_pillbox_state { - double radius; - double basis[9]; -}; - -struct ran_dirac_state { - double dir[3]; -}; - -/* One single type for all distributions. Only the state type depends on the - * distribution type */ -struct ssol_ranst_sun_dir { - double*(*get) - (const struct ssol_ranst_sun_dir* ran, struct ssp_rng* rng, double dir[3]); - union { - struct ran_buie_state buie; - struct ran_pillbox_state pillbox; - struct ran_dirac_state dirac; - } state; - - ref_T ref; - struct mem_allocator* allocator; -}; - -/******************************************************************************* - * Helper functions - ******************************************************************************/ -static void -distrib_sun_release(ref_T* ref) -{ - struct ssol_ranst_sun_dir* ran; - ASSERT(ref); - ran = CONTAINER_OF(ref, struct ssol_ranst_sun_dir, ref); - MEM_RM(ran->allocator, ran); -} - -/******************************************************************************* - * Buie distribution - ******************************************************************************/ -static double -chi_value(const double csr) -{ - if (csr > 0.145) { - return -0.04419909985804843 + - csr * (1.401323894233574 + - csr * (-0.3639746714505299 + - csr * (-0.9579768560161194 + 1.1550475450828657 * csr))); - } - if (csr > 0.035) { - return 0.022652077593662934 + - csr * (0.5252380349996234 + - csr * (2.5484334534423887 - 0.8763755326550412 * csr)); - } - return 0.004733749294807862 + - csr * (4.716738065192151 + - csr * (-463.506669149804 + - csr * (24745.88727411664 + - csr * (-606122.7511711778 + 5521693.445014727 * csr)))); -} - -static double -phi_solar_disk(const double theta) -{ - /* The parameter theta is the zenith angle in radians */ - return cos(326 * theta) / cos(308 * theta); -} - -static double -phi_circum_solor_region(const double theta, const struct ran_buie_state* state) -{ - /* The parameter theta is the zenith angle in radians */ - return state->etokTimes1000toGamma * pow(theta, state->gamma); -} - -static double -phi(const double theta, const struct ran_buie_state* state) -{ - /* The parameter theta is the zenith angle in radians */ - if (theta < state->thetaSD) return phi_solar_disk(theta); - else return phi_circum_solor_region(theta, state); -} - -static double -pdf_theta(const double theta, const struct ran_buie_state* state) -{ - /* The parameter theta is the zenith angle in radians */ - return state->alpha * phi(theta, state) * sin(theta); -} - -static double -gamma_value(const double chi) -{ - /* gamma is the gradient of the 2nd part of the curve in log/log space*/ - return 2.2 * log(0.52 * chi) * pow(chi, 0.43) - 0.1; -} - -static double -k_value(const double chi) -{ - /* k is the intercept of the 2nd part of the curve at an angular displacement - * of zero */ - return 0.9 * log(13.5 * chi) * pow(chi, -0.3); -} - -static double -integral_B - (const double k, - const double gamma, - const double thetaCS, - const double thetaSD) -{ - const double g2 = gamma + 2.0; - return exp(k) * pow(1000, gamma) / g2 * (pow(thetaCS, g2) - pow(thetaSD, g2)); -} - -static double -proba_rect1 - (const double widthR1, - const double heightR1, - const double widthR2, - const double heightR2) -{ - const double areaR1 = widthR1 * heightR1; - const double areaR2 = widthR2 * heightR2; - return areaR1 / (areaR1 + areaR2); -} - -static double -zenith_angle(struct ssp_rng* rng, const struct ran_buie_state* state) -{ - double theta; - double value; - ASSERT(state); - - do { - if (ssp_rng_canonical(rng) < state->proba_rect1) { - theta = ssp_rng_canonical(rng) * state->thetaSD; - value = ssp_rng_canonical(rng) * state->hRect1; - } - else { - theta = state->thetaSD + ssp_rng_canonical(rng) * state->deltaThetaCSSD; - value = ssp_rng_canonical(rng) * state->hRect2; - } - } while (value > pdf_theta(theta, state)); - - return theta; -} - -static void -fill_buie_state(struct ran_buie_state* s, const double p) -{ - double thetaCS, integralA, chi, k, integralB; - ASSERT(s); - ASSERT(0 <= p && p < 1); - - s->thetaSD = 0.00465; - thetaCS = 0.0436; - s->deltaThetaCSSD = thetaCS - s->thetaSD; - integralA = 9.224724736098827E-6; - chi = chi_value(p); - k = k_value(chi); - s->gamma = gamma_value(chi); - s->etokTimes1000toGamma = exp(k) * pow(1000, s->gamma); - integralB = integral_B(k, s->gamma, thetaCS, s->thetaSD); - s->alpha = 1 / (integralA + integralB); - /* The value 0.0038915695846209047 radians is the value for which the - * probability density function is a maximum. The value of the maximum varies - * with the circumsolar ratio, but the value of theta where the maximum is - * obtained remains fixed. */ - s->hRect1 = 1.001 * pdf_theta(0.0038915695846209047, s); - s->hRect2 = pdf_theta(s->thetaSD, s); - s->proba_rect1 = proba_rect1(s->thetaSD, s->hRect1, s->deltaThetaCSSD, s->hRect2); -} - -static double* -ssol_ran_buie_get - (const struct ssol_ranst_sun_dir* ran, struct ssp_rng* rng, double dir[3]) -{ - double phi, theta, sinTheta, cosTheta, cosPhi, sinPhi; - ASSERT(ran->state.buie.thetaSD > 0); - phi = ssp_rng_uniform_double(rng, 0, 2 * PI); - theta = zenith_angle(rng, &ran->state.buie); - sinTheta = sin(theta); - cosTheta = cos(theta); - cosPhi = cos(phi); - sinPhi = sin(phi); - dir[0] = sinTheta * sinPhi; - dir[1] = -cosTheta; - dir[2] = sinTheta * cosPhi; - d33_muld3(dir, ran->state.buie.basis, dir); - return dir; -} - -/******************************************************************************* - * Pillbox distribution - ******************************************************************************/ -static double* -ssol_ran_pillbox_get - (const struct ssol_ranst_sun_dir* ran, - struct ssp_rng* rng, - double dir[3]) -{ - double pt[3]; - - ASSERT(ran && ran->state.pillbox.radius > 0); - ssp_ran_uniform_disk(rng, ran->state.pillbox.radius, pt); - pt[2] = 1; - d33_muld3(dir, ran->state.pillbox.basis, pt); - d3_normalize(dir, pt); - return dir; -} - -/******************************************************************************* -* Dirac distribution -******************************************************************************/ -static double* -ssol_ran_dirac_get -(const struct ssol_ranst_sun_dir* ran, - struct ssp_rng* rng, - double dir[3]) -{ - (void) rng; - ASSERT(d3_is_normalized(ran->state.dirac.dir)); - d3_set(dir, ran->state.dirac.dir); - return dir; -} - -/******************************************************************************* - * Local functions - ******************************************************************************/ -res_T -ssol_ranst_sun_dir_create - (struct mem_allocator* allocator, - struct ssol_ranst_sun_dir** out_ran) -{ - struct ssol_ranst_sun_dir* ran = NULL; - - if (!out_ran) return RES_BAD_ARG; - - allocator = allocator ? allocator : &mem_default_allocator; - - ran = MEM_CALLOC(allocator, 1, sizeof(struct ssol_ranst_sun_dir)); - if (!ran) return RES_MEM_ERR; - - ref_init(&ran->ref); - ran->allocator = allocator; - *out_ran = ran; - - return RES_OK; -} - -res_T -ssol_ranst_sun_dir_ref_get(struct ssol_ranst_sun_dir* ran) -{ - if (!ran) return RES_BAD_ARG; - ref_get(&ran->ref); - return RES_OK; -} - -res_T -ssol_ranst_sun_dir_ref_put(struct ssol_ranst_sun_dir* ran) -{ - if (!ran) return RES_BAD_ARG; - ref_put(&ran->ref, distrib_sun_release); - return RES_OK; -} - -double* -ssol_ranst_sun_dir_get - (const struct ssol_ranst_sun_dir* ran, - struct ssp_rng* rng, - double dir[3]) -{ - ASSERT(ran); - return ran->get(ran, rng, dir); -} - -res_T -ssol_ranst_sun_dir_buie_setup - (struct ssol_ranst_sun_dir* ran, - double param, - const double dir[3]) -{ - const double minCRSValue = 1E-6; - const double maxCRSValue = 0.849; - if (!ran || !dir || param < minCRSValue || param > maxCRSValue) - return RES_BAD_ARG; - - ran->get = ssol_ran_buie_get; - d33_basis(ran->state.buie.basis, dir); - fill_buie_state(&ran->state.buie, param); - return RES_OK; -} - -res_T -ssol_ranst_sun_dir_pillbox_setup - (struct ssol_ranst_sun_dir* ran, - double aperture, - const double dir[3]) -{ - double radius; - if (!ran || !dir || aperture <= 0 || aperture >= PI ) - return RES_BAD_ARG; - radius = tan(0.5 * aperture); - ran->get = ssol_ran_pillbox_get; - ran->state.pillbox.radius = radius; - d33_basis(ran->state.pillbox.basis, dir); - return RES_OK; -} - -res_T -ssol_ranst_sun_dir_dirac_setup - (struct ssol_ranst_sun_dir* ran, - const double dir[3]) -{ - if (!ran || !dir) return RES_BAD_ARG; - if (0 == d3_normalize(ran->state.dirac.dir, dir)) - /* zero vector */ - return RES_BAD_ARG; - ran->get = ssol_ran_dirac_get; - return RES_OK; -} - diff --git a/src/ssol_distributions_c.h b/src/ssol_distributions_c.h @@ -1,66 +0,0 @@ -/* Copyright (C) CNRS 2016 - * - * 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/>. */ - -#ifndef SSOL_DISTRIBUTIONS_C_H -#define SSOL_DISTRIBUTIONS_C_H - -/******************************************************************************* - * Sun distributions - ******************************************************************************/ - -/* Forward declaration of opaque types */ -struct ssol_ranst_sun_dir; - -/* Forward declaration of external types */ -struct ssp_rng; -struct mem_allocator; - -res_T -ssol_ranst_sun_dir_create - (struct mem_allocator* allocator, - struct ssol_ranst_sun_dir** ran); - -res_T -ssol_ranst_sun_dir_buie_setup - (struct ssol_ranst_sun_dir* ran, - double param, - const double dir[3]); - -res_T -ssol_ranst_sun_dir_pillbox_setup - (struct ssol_ranst_sun_dir* ran, - double aperture, /* apparent size in radians */ - const double dir[3]); - -res_T -ssol_ranst_sun_dir_dirac_setup - (struct ssol_ranst_sun_dir* ran, - const double dir[3]); - -res_T -ssol_ranst_sun_dir_ref_get - (struct ssol_ranst_sun_dir* ran); - -res_T -ssol_ranst_sun_dir_ref_put - (struct ssol_ranst_sun_dir* ran); - -double* -ssol_ranst_sun_dir_get - (const struct ssol_ranst_sun_dir* ran, - struct ssp_rng* rng, - double dir[3]); - -#endif /* SSOL_DISTRIBUTIONS_C_H */ diff --git a/src/ssol_ranst_sun_dir.c b/src/ssol_ranst_sun_dir.c @@ -0,0 +1,365 @@ +/* Copyright (C) CNRS 2016 + * + * 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 "ssol_ranst_sun_dir.h" + +#include <star/ssp.h> + +#include <rsys/double33.h> +#include <rsys/math.h> +#include <rsys/mem_allocator.h> +#include <rsys/rsys.h> +#include <rsys/ref_count.h> + +/******************************************************************************* + * Distributions types + ******************************************************************************/ +struct ran_buie_state { + double thetaSD; + double deltaThetaCSSD; + double gamma; + double etokTimes1000toGamma; + double alpha; + double hRect1; + double hRect2; + double proba_rect1; + double basis[9]; +}; + +struct ran_pillbox_state { + double radius; + double basis[9]; +}; + +struct ran_dirac_state { + double dir[3]; +}; + +/* One single type for all distributions. Only the state type depends on the + * distribution type */ +struct ranst_sun_dir { + double*(*get) + (const struct ranst_sun_dir* ran, struct ssp_rng* rng, double dir[3]); + union { + struct ran_buie_state buie; + struct ran_pillbox_state pillbox; + struct ran_dirac_state dirac; + } state; + + ref_T ref; + struct mem_allocator* allocator; +}; + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static void +distrib_sun_release(ref_T* ref) +{ + struct ranst_sun_dir* ran; + ASSERT(ref); + ran = CONTAINER_OF(ref, struct ranst_sun_dir, ref); + MEM_RM(ran->allocator, ran); +} + +/******************************************************************************* + * Buie distribution + ******************************************************************************/ +static FINLINE double +chi_value(const double csr) +{ + if (csr > 0.145) { + return -0.04419909985804843 + + csr * (1.401323894233574 + + csr * (-0.3639746714505299 + + csr * (-0.9579768560161194 + 1.1550475450828657 * csr))); + } + if (csr > 0.035) { + return 0.022652077593662934 + + csr * (0.5252380349996234 + + csr * (2.5484334534423887 - 0.8763755326550412 * csr)); + } + return 0.004733749294807862 + + csr * (4.716738065192151 + + csr * (-463.506669149804 + + csr * (24745.88727411664 + + csr * (-606122.7511711778 + 5521693.445014727 * csr)))); +} + +static FINLINE double +phi_solar_disk(const double theta) +{ + /* The parameter theta is the zenith angle in radians */ + return cos(326 * theta) / cos(308 * theta); +} + +static FINLINE double +phi_circum_solor_region(const double theta, const struct ran_buie_state* state) +{ + /* The parameter theta is the zenith angle in radians */ + return state->etokTimes1000toGamma * pow(theta, state->gamma); +} + +static FINLINE double +phi(const double theta, const struct ran_buie_state* state) +{ + /* The parameter theta is the zenith angle in radians */ + if (theta < state->thetaSD) return phi_solar_disk(theta); + else return phi_circum_solor_region(theta, state); +} + +static FINLINE double +pdf_theta(const double theta, const struct ran_buie_state* state) +{ + /* The parameter theta is the zenith angle in radians */ + return state->alpha * phi(theta, state) * sin(theta); +} + +static FINLINE double +gamma_value(const double chi) +{ + /* gamma is the gradient of the 2nd part of the curve in log/log space*/ + return 2.2 * log(0.52 * chi) * pow(chi, 0.43) - 0.1; +} + +static FINLINE double +k_value(const double chi) +{ + /* k is the intercept of the 2nd part of the curve at an angular displacement + * of zero */ + return 0.9 * log(13.5 * chi) * pow(chi, -0.3); +} + +static FINLINE double +integral_B + (const double k, + const double gamma, + const double thetaCS, + const double thetaSD) +{ + const double g2 = gamma + 2.0; + return exp(k) * pow(1000, gamma) / g2 * (pow(thetaCS, g2) - pow(thetaSD, g2)); +} + +static FINLINE double +proba_rect1 + (const double widthR1, + const double heightR1, + const double widthR2, + const double heightR2) +{ + const double areaR1 = widthR1 * heightR1; + const double areaR2 = widthR2 * heightR2; + return areaR1 / (areaR1 + areaR2); +} + +static FINLINE double +zenith_angle(struct ssp_rng* rng, const struct ran_buie_state* state) +{ + double theta; + double value; + ASSERT(state); + + do { + if (ssp_rng_canonical(rng) < state->proba_rect1) { + theta = ssp_rng_canonical(rng) * state->thetaSD; + value = ssp_rng_canonical(rng) * state->hRect1; + } + else { + theta = state->thetaSD + ssp_rng_canonical(rng) * state->deltaThetaCSSD; + value = ssp_rng_canonical(rng) * state->hRect2; + } + } while (value > pdf_theta(theta, state)); + + return theta; +} + +static void +fill_buie_state(struct ran_buie_state* s, const double p) +{ + double thetaCS, integralA, chi, k, integralB; + ASSERT(s); + ASSERT(0 <= p && p < 1); + + s->thetaSD = 0.00465; + thetaCS = 0.0436; + s->deltaThetaCSSD = thetaCS - s->thetaSD; + integralA = 9.224724736098827E-6; + chi = chi_value(p); + k = k_value(chi); + s->gamma = gamma_value(chi); + s->etokTimes1000toGamma = exp(k) * pow(1000, s->gamma); + integralB = integral_B(k, s->gamma, thetaCS, s->thetaSD); + s->alpha = 1 / (integralA + integralB); + /* The value 0.0038915695846209047 radians is the value for which the + * probability density function is a maximum. The value of the maximum varies + * with the circumsolar ratio, but the value of theta where the maximum is + * obtained remains fixed. */ + s->hRect1 = 1.001 * pdf_theta(0.0038915695846209047, s); + s->hRect2 = pdf_theta(s->thetaSD, s); + s->proba_rect1 = proba_rect1 + (s->thetaSD, s->hRect1, s->deltaThetaCSSD, s->hRect2); +} + +static double* +ran_buie_get + (const struct ranst_sun_dir* ran, struct ssp_rng* rng, double dir[3]) +{ + double phi, theta, sinTheta, cosTheta, cosPhi, sinPhi; + ASSERT(ran->state.buie.thetaSD > 0); + phi = ssp_rng_uniform_double(rng, 0, 2 * PI); + theta = zenith_angle(rng, &ran->state.buie); + sinTheta = sin(theta); + cosTheta = cos(theta); + cosPhi = cos(phi); + sinPhi = sin(phi); + dir[0] = sinTheta * sinPhi; + dir[1] = -cosTheta; + dir[2] = sinTheta * cosPhi; + d33_muld3(dir, ran->state.buie.basis, dir); + return dir; +} + +/******************************************************************************* + * Pillbox random variate + ******************************************************************************/ +static double* +ran_pillbox_get + (const struct ranst_sun_dir* ran, + struct ssp_rng* rng, + double dir[3]) +{ + double pt[3]; + + ASSERT(ran && ran->state.pillbox.radius > 0); + ssp_ran_uniform_disk(rng, ran->state.pillbox.radius, pt); + pt[2] = 1; + d33_muld3(dir, ran->state.pillbox.basis, pt); + d3_normalize(dir, pt); + return dir; +} + +/******************************************************************************* + * Dirac distribution + ******************************************************************************/ +static double* +ran_dirac_get +(const struct ranst_sun_dir* ran, + struct ssp_rng* rng, + double dir[3]) +{ + (void) rng; + ASSERT(d3_is_normalized(ran->state.dirac.dir)); + d3_set(dir, ran->state.dirac.dir); + return dir; +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +ranst_sun_dir_create + (struct mem_allocator* allocator, + struct ranst_sun_dir** out_ran) +{ + struct ranst_sun_dir* ran = NULL; + + if (!out_ran) return RES_BAD_ARG; + + allocator = allocator ? allocator : &mem_default_allocator; + + ran = MEM_CALLOC(allocator, 1, sizeof(struct ranst_sun_dir)); + if (!ran) return RES_MEM_ERR; + + ref_init(&ran->ref); + ran->allocator = allocator; + *out_ran = ran; + + return RES_OK; +} + +res_T +ranst_sun_dir_ref_get(struct ranst_sun_dir* ran) +{ + if (!ran) return RES_BAD_ARG; + ref_get(&ran->ref); + return RES_OK; +} + +res_T +ranst_sun_dir_ref_put(struct ranst_sun_dir* ran) +{ + if (!ran) return RES_BAD_ARG; + ref_put(&ran->ref, distrib_sun_release); + return RES_OK; +} + +double* +ranst_sun_dir_get + (const struct ranst_sun_dir* ran, + struct ssp_rng* rng, + double dir[3]) +{ + ASSERT(ran); + return ran->get(ran, rng, dir); +} + +res_T +ranst_sun_dir_buie_setup + (struct ranst_sun_dir* ran, + double param, + const double dir[3]) +{ + const double minCRSValue = 1E-6; + const double maxCRSValue = 0.849; + if (!ran || !dir || param < minCRSValue || param > maxCRSValue) + return RES_BAD_ARG; + + ran->get = ran_buie_get; + d33_basis(ran->state.buie.basis, dir); + fill_buie_state(&ran->state.buie, param); + return RES_OK; +} + +res_T +ranst_sun_dir_pillbox_setup + (struct ranst_sun_dir* ran, + double aperture, + const double dir[3]) +{ + double radius; + if (!ran || !dir || aperture <= 0 || aperture >= PI ) + return RES_BAD_ARG; + radius = tan(0.5 * aperture); + ran->get = ran_pillbox_get; + ran->state.pillbox.radius = radius; + d33_basis(ran->state.pillbox.basis, dir); + return RES_OK; +} + +res_T +ranst_sun_dir_dirac_setup + (struct ranst_sun_dir* ran, + const double dir[3]) +{ + if (!ran || !dir) return RES_BAD_ARG; + if (0 == d3_normalize(ran->state.dirac.dir, dir)) + /* zero vector */ + return RES_BAD_ARG; + ran->get = ran_dirac_get; + return RES_OK; +} + diff --git a/src/ssol_ranst_sun_dir.h b/src/ssol_ranst_sun_dir.h @@ -0,0 +1,63 @@ +/* Copyright (C) CNRS 2016 + * + * 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/>. */ + +#ifndef SSOL_RANST_SUN_DIR_H +#define SSOL_RANST_SUN_DIR_H + +/* Random variate state of a sun direction */ +struct ranst_sun_dir; + +/* External types */ +struct ssp_rng; +struct mem_allocator; + +extern LOCAL_SYM res_T +ranst_sun_dir_create + (struct mem_allocator* allocator, + struct ranst_sun_dir** ran); + +extern LOCAL_SYM res_T +ranst_sun_dir_buie_setup + (struct ranst_sun_dir* ran, + const double param, + const double dir[3]); + +extern LOCAL_SYM res_T +ranst_sun_dir_pillbox_setup + (struct ranst_sun_dir* ran, + const double aperture, /* Apparent size in radians */ + const double dir[3]); + +extern LOCAL_SYM res_T +ranst_sun_dir_dirac_setup + (struct ranst_sun_dir* ran, + const double dir[3]); + +extern LOCAL_SYM res_T +ranst_sun_dir_ref_get + (struct ranst_sun_dir* ran); + +extern LOCAL_SYM res_T +ranst_sun_dir_ref_put + (struct ranst_sun_dir* ran); + +extern LOCAL_SYM double* +ranst_sun_dir_get + (const struct ranst_sun_dir* ran, + struct ssp_rng* rng, + double dir[3]); + +#endif /* SSOL_RANST_SUN_DIR_H */ + diff --git a/src/ssol_solver.c b/src/ssol_solver.c @@ -188,14 +188,14 @@ set_sun_distributions /* then the direction distribution */ switch (sun->type) { case SUN_DIRECTIONAL: - res = ssol_ranst_sun_dir_dirac_setup(data->sun_dir_ran, sun->direction); + res = ranst_sun_dir_dirac_setup(data->sun_dir_ran, sun->direction); break; case SUN_PILLBOX: - res = ssol_ranst_sun_dir_pillbox_setup + res = ranst_sun_dir_pillbox_setup (data->sun_dir_ran, sun->data.pillbox.aperture, sun->direction); break; case SUN_CSR: - res = ssol_ranst_sun_dir_buie_setup + res = ranst_sun_dir_buie_setup (data->sun_dir_ran, sun->data.csr.ratio, sun->direction); break; default: diff --git a/src/ssol_solver_c.h b/src/ssol_solver_c.h @@ -16,7 +16,7 @@ #ifndef SSOL_SOLVER_C_H #define SSOL_SOLVER_C_H -#include "ssol_distributions_c.h" +#include "ssol_ranst_sun_dir.h" #include <rsys/ref_count.h> #include <rsys/list.h> @@ -46,7 +46,7 @@ struct solver_data { /* the 3D scene used for raytracing */ struct s3d_scene *scene; /* the random distributions for sun sampling */ - struct ssol_ranst_sun_dir* sun_dir_ran; + struct ranst_sun_dir* sun_dir_ran; struct ssp_ranst_piecewise_linear* sun_spectrum_ran; };