star-gs

Literate program for a geometric sensitivity calculation
git clone git://git.meso-star.fr/star-gs.git
Log | Files | Refs | README | LICENSE

commit f62c7cca5ed19150bc9ecada6d88f68e4d0ffefc
parent 0a6987709d56bbf6f91afc0f81e5db769dd4c186
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed, 15 Feb 2023 21:51:55 +0100

Rédige le texte de l'annexe E

Cette annexe contient la structure de mise en oeuvre et n'avait
jusqu'alors que les seuls listings de code C. Dans cette validation est
ajouté le texte qui les accompagne tout en changeant l'ordre dans lequel
le code C est présenté.

Diffstat:
MMakefile | 1+
Msrc/sgs.c | 19+++++--------------
Msrc/sgs_c.h | 3+--
Msrc/sgs_compute_sensitivity_translation.nw | 149++++++++++++++++++++++++++++++++++++++++++++++++++-----------------------------
4 files changed, 101 insertions(+), 71 deletions(-)

diff --git a/Makefile b/Makefile @@ -62,6 +62,7 @@ sgs: $(OBJ) @$(SHELL) configure.sh program && echo "config done" > $@ || exit 1 $(TMPC): $(NOWEB) config.mk + @echo "TANGLE $@" @$(TANGLE) $(TANGLE_OPTS)\ -R$$(c="$@" && echo "$${c##*/}" | sed 's/_/\\_/g') $(NOWEB) | cpif $@ diff --git a/src/sgs.c b/src/sgs.c @@ -36,27 +36,18 @@ static res_T setup_scene(struct sgs* sgs, const struct sgs_args* args) { - double sz[3]; res_T res = RES_OK; res = sgs_geometry_box_create(sgs, &args->geom, &sgs->scene.geom); if(res != RES_OK) goto error; - d3_sub(sz, args->geom.upper, args->geom.lower); - - /* Setup the 2D size of the emitting surface (YZ frame) */ - sgs->scene.Se_sz[0] = sz[1]; - sgs->scene.Se_sz[1] = sz[2]; - - /* Setup the 2D size of the reflective surface (XY frame) */ - sgs->scene.Sr_sz[0] = sz[0]; - sgs->scene.Sr_sz[1] = sz[1]; + d3_sub(sgs->scene.D, args->geom.upper, args->geom.lower); /* Setup the 2D lower/upper bounds of the receiver */ - sgs->scene.recv_min[0] = sz[0]*1.0/8.0 + args->geom.lower[0]; - sgs->scene.recv_min[1] = sz[1]*3.0/8.0 + args->geom.lower[1]; - sgs->scene.recv_max[0] = sz[0]*3.0/8.0 + args->geom.lower[0]; - sgs->scene.recv_max[1] = sz[1]*5.0/8.0 + args->geom.lower[1]; + sgs->scene.recv_min[0] = sgs->scene.D[0]*1.0/8.0 + args->geom.lower[0]; + sgs->scene.recv_min[1] = sgs->scene.D[1]*3.0/8.0 + args->geom.lower[1]; + sgs->scene.recv_max[0] = sgs->scene.D[0]*3.0/8.0 + args->geom.lower[0]; + sgs->scene.recv_max[1] = sgs->scene.D[1]*5.0/8.0 + args->geom.lower[1]; exit: return res; diff --git a/src/sgs_c.h b/src/sgs_c.h @@ -29,8 +29,7 @@ struct sgs_scene { struct sgs_geometry* geom; - double Sr_sz[2]; /* 2D size of the reflective surface (in XY frame) */ - double Se_sz[2]; /* 2D size of the emitting surface (in YZ frame) */ + double D[3]; /* Dimension of the box */ double recv_min[2]; /* 2D lower bound of the receiver (in XY frame) */ double recv_max[2]; /* 2D upper bound of the receiver (in XY frame) */ }; diff --git a/src/sgs_compute_sensitivity_translation.nw b/src/sgs_compute_sensitivity_translation.nw @@ -211,6 +211,7 @@ paroi: \begin{equation} S_b = L^{eq}(T) \lbrack 1- \cos(2 \pi \frac{x}{Dx}) \rbrack \lbrack 1- \cos(2 \pi \frac{y}{Dy}) \rbrack +\label{eq:S_b} \end{equation} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -672,7 +673,7 @@ Sb_sensib = + rho * proj_chi_e.beta * d_u_cs_Sb; /* Poids de sensibilité */ -sensib = Sb_sensib * PI * get_Sr_area(scene); +sensib = Sb_sensib * PI * get_Sr(scene); @ <<renvoyer le poids>>= @@ -973,7 +974,7 @@ double weight_flux_part_spec; @ <<calcul du poids>>= -weight_flux_part_spec = Sb * rho * PI * get_Sr_area(scene); +weight_flux_part_spec = Sb * rho * PI * get_Sr(scene); @ <<renvoyer le poids>>= @@ -992,6 +993,7 @@ weight_flux_part_spec = 0; % Annexe fonction de calcul %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Fonction de calcul} +\label{sec:fonction_de_calul} \paragraph{} Le programme présenté jusqu'alors s'est concentré sur la mise en {\oe}uvre de @@ -1147,35 +1149,28 @@ run_realisation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Structure de mise oeuvre} -\paragraph{TODO} Décrire la structure du fichier C issu de noweb et son -articulation avec le reste des sources écrites en C notamment lors de -l'énumération des fichiers à inclure. +\paragraph{} +Cette partie décrit la structure du fichier C dans lequel l'algorithme de +calcul de sensibilité est mis en oeuvre: <<sgs\_compute\_sensitivity\_translation.c>>= <<liste des inclusions>> -<<constantes du système>> +<<constantes>> <<fonctions utilitaires>> <<fonction de réalisation>> <<calculer la sensibilité à la translation>> @ -<<constantes du système>>= -enum {X, Y, Z}; -@ - -<<liste des inclusions>>= -#include "sgs_c.h" -#include "sgs_geometry.h" -#include "sgs_log.h" - -#include <star/smc.h> -#include <star/ssp.h> - -#include <rsys/cstr.h> -#include <rsys/double2.h> -@ +\paragraph{} +En plus des fonctions de calcul écrites dans les sections précédents, notre +fichier contient des fonctions utilitaires notament utilisées pour interroger +les propriétés physiques du système (figure~\ref{fig:configuration}) telles que +$\rho$, la réflectivité de la paroi du haut (equation~\ref{eq:rho}), ou $S_b$, +l'émission thermique de la paroi de droite (équation~\ref{eq:S_b}). Le calcul +de leur gradient surfacique, utilisé pour évaluer leur dérivée surfacique +(section~\ref{subsec:poids}), est également ajouté comme fonction utilitaire: <<fonctions utilitaires>>= static double @@ -1186,8 +1181,8 @@ get_rho ASSERT(scene && pos); return 0.25 - * (1 - cos(2*PI*pos[X]/scene->Sr_sz[X])) - * (1 - cos(2*PI*pos[Y]/scene->Sr_sz[Y])); + * (1 - cos(2*PI*pos[X]/scene->D[X])) + * (1 - cos(2*PI*pos[Y]/scene->D[Y])); } static void @@ -1199,24 +1194,22 @@ get_grad_rho ASSERT(scene && pos && grad); grad[X] = 0.25 - * (((2*PI)/scene->Sr_sz[X])*sin(2*PI*pos[X]/scene->Sr_sz[X])) - * (1 - cos(2*PI*pos[Y]/scene->Sr_sz[Y])); + * (((2*PI)/scene->D[X])*sin(2*PI*pos[X]/scene->D[X])) + * (1 - cos(2*PI*pos[Y]/scene->D[Y])); grad[Y] = 0.25 - * (((2*PI)/scene->Sr_sz[Y])*sin(2*PI*pos[Y]/scene->Sr_sz[Y])) - * (1 - cos(2*PI*pos[X]/scene->Sr_sz[X])); + * (((2*PI)/scene->D[Y])*sin(2*PI*pos[Y]/scene->D[Y])) + * (1 - cos(2*PI*pos[X]/scene->D[X])); } -@ -<<fonctions utilitaires>>= static double get_Sb (const struct sgs_scene* scene, const double pos[2]) { - ASSERT(position && size); + ASSERT(scene && pos); return - (1 - cos(2*PI*pos[X]/scene->Se_sz[X])) - * (1 - cos(2*PI*pos[Y]/scene->Se_sz[Y])); + (1 - cos(2*PI*pos[X]/scene->D[Y])) + * (1 - cos(2*PI*pos[Y]/scene->D[Z])); } static void @@ -1228,23 +1221,49 @@ get_grad_Sb ASSERT(scene && pos && grad); grad[X] = - (((2*PI)/scene->Se_sz[X])*sin(2*PI*pos[X]/scene->Se_sz[X])) - * (1 - cos(2*PI*pos[Y]/scene->Se_sz[Y])); + (((2*PI)/scene->D[Y])*sin(2*PI*pos[X]/scene->D[Y])) + * (1 - cos(2*PI*pos[Y]/scene->D[Z])); grad[Y] = - (((2*PI)/scene->Se_sz[Y])*sin(2*PI*pos[Y]/scene->Se_sz[Y])) - * (1 - cos(2*PI*pos[X]/scene->Se_sz[X])); + (((2*PI)/scene->D[Z])*sin(2*PI*pos[Y]/scene->D[Z])) + * (1 - cos(2*PI*pos[X]/scene->D[Y])); } @ +\paragraph{} +Autres fonctions utilitaires, les fonctions utilisées lors du suivie de chemins +pour déterminer si un rayon a intersecté le récepteur (fonction +[[hit_receiver]]) ou la source (fonction [[hit_source]]). Pour la source il +suffit de s'assurer que le rayon intersecte le côté droit de la boîte, ici +identifié par la constante [[SGS_SURFACE_X_MAX]]. + <<fonctions utilitaires>>= -static double -get_Sr_area(const struct sgs_scene* scene) +static int +hit_source + (const struct sgs_scene* scene, + const double ray_org[3], + const double ray_dir[3], + const struct sgs_hit* hit) { - ASSERT(scene); - return scene->Sr_sz[X] * scene->Sr_sz[Y]; + ASSERT(scene && ray_org && ray_dir && hit); + /* Éviter l'avertissement de compilation "variable inutilisée" */ + (void)scene, (void)ray_org, (void)ray_dir; + + if(SGS_HIT_NONE(hit) || hit->surface != SGS_SURFACE_X_MAX) { + return 0; + } else { + return 1; /* L'intersection a lieu sur la source */ + } } @ +\paragraph{} +Pour le récepteur on teste d'abord si l'intersection a lieu sur la paroi de la +boîte sur laquelle celui-ci est positionné, une paroi identifiée dans le code +par la constante [[SGS_SURFACE_Z_MIN]]. Reste alors à déterminer si +l'intersection se situe sur le récepteur lui même. Pour cela il nous suffit de +tester si la position d'intersection appartient au sous domaine de la surface +sur lequel le récepteur est défini. + <<fonctions utilitaires>>= static int hit_receiver @@ -1261,6 +1280,7 @@ hit_receiver return 0; } + /* Calculer la position de l'intersection */ hit_pos[X] = ray_org[X] + hit->distance*ray_dir[X]; hit_pos[Y] = ray_org[Y] + hit->distance*ray_dir[Y]; hit_pos[Z] = ray_org[Z] + hit->distance*ray_dir[Z]; @@ -1277,26 +1297,45 @@ hit_receiver } @ +\paragraph{} +Dernière fonction utilitaire à écrire, la fonction qui calcule la surface $S_r$ +de la paroi spéculaire: + <<fonctions utilitaires>>= -static int -hit_source - (const struct sgs_scene* scene, - const double ray_org[3], - const double ray_dir[3], - const struct sgs_hit* hit) +static double +get_Sr(const struct sgs_scene* scene) { - ASSERT(scene && ray_org && ray_dir && hit); - /* Éviter l'avertissement de compilation "variable inutilisée" */ - (void)scene, (void)ray_org, (void)ray_dir; - - if(SGS_HIT_NONE(hit) || hit->surface != SGS_SURFACE_X_MAX) { - return 0; - } else { - return 1; - } + ASSERT(scene); + return scene->D[X] * scene->D[Y]; } @ +\paragraph{} +On complète notre fichier en définissant les constantes [[X]], [[Y]] et [[Z]] +utilisées tout du long pour simplifier la lecture du code lors des accès aux +éléments d'un vecteur. + +<<constantes>>= +enum {X, Y, Z}; +@ + +\paragraph{} +Enfin, on liste en début des sources les fichiers d'en-tête sur lequel notre +code s'appuie. + +<<liste des inclusions>>= +#include "sgs_c.h" +#include "sgs_geometry.h" +#include "sgs_log.h" + +#include <star/smc.h> /* Calcul Monte-Carlo */ +#include <star/ssp.h> /* Générateur de nombre aléatoire et distributions */ + +/* Manipuler des vecteurs de double à 2 et 3 dimensions */ +#include <rsys/double2.h> +#include <rsys/double3.h> +@ + \bibliographystyle{apalike} \bibliography{biblio}