star-gs

sgs_compute_sensitivity_translation.nw (69105B)

---

 % Copyright (C) 2021-2023 Centre National de la Recherche Scientifique
 % Copyright (C) 2021-2023 INSA Lyon
 % Copyright (C) 2021-2023 Institut Mines Télécom Albi-Carmaux
 % Copyright (C) 2021-2023 |Méso|Star> (contact@meso-star.com)
 % Copyright (C) 2021-2023 Institut Pascal
 % Copyright (C) 2021-2023 PhotonLyX (info@photonlyx.com)
 % Copyright (C) 2021-2023 Université de Lorraine
 % Copyright (C) 2021-2023 Université Paul Sabatier
 % Copyright (C) 2021-2023 Université Toulouse - Jean Jaurès
 %
 % 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/>.
 
 \documentclass[]{article}
 \usepackage[utf8]{inputenc}
 \usepackage[french]{babel}
 \usepackage[T1]{fontenc} % Pour la césure des mots avec accents
 \usepackage{url}
 \usepackage{amsfonts} % \mathbb
 \usepackage{amsmath}
 \usepackage{graphicx}
 \usepackage{noweb}
 \usepackage{tikz}
 \usepackage{pgfplots} % \begin{axis}
 \usetikzlibrary{patterns}
 \noweboptions{smallcode,longchunks}
 
 % Remove clearpage after nwfilename
 \makeatletter
 \def\nw@laterfilename#1{\endgroup \markboth{#1}{#1}}
 \makeatother
 
 % Avoid a page break after a code chunk
 \def\nwendcode{\endtrivlist \endgroup}
 \let\nwdocspar=\par
 %\def\nwendcode{\endtrivlist \endgroup \vfil\penalty10\vfilneg}
 %\let\nwdocspar=\smallbreak
 
 %Math et symboles
 \newcommand{\PI}{\ddot \pi}
 \newcommand{\etc}{\textit{etc.}}
 \newcommand{\ie}{\textit{i.e.}}
 \newcommand{\eg}{\textit{e.g.}}
 
 \begin{document}
 \pagestyle{noweb}
 
 \title{Sensibilité à la translation}
 \maketitle
 
 Le but du présent document est d'illustrer la mise en {\oe}uvre algorithmique
 d'un calcul de sensibilité géométrique sur l'exemple simple d'un
 parallélépipède (figure~\ref{fig:configuration}). Nous nous intéressons ici à
 la déformation géométrique liée à la translation de sa paroi supérieure et
 étudions l'impact de cette translation sur le flux reçu par un récepteur situé
 sur sa paroi inférieure.  Dans cet exercice, la sensibilité du flux est estimée
 par un algorithme de Monte Carlo analogue à la physique du transport de la
 sensibilité géométrique.  Cette pratique des algorithmes de Monte Carlo est
 largement utilisée en transferts radiatifs et consiste à imiter numériquement
 le transport de photons, en échantillonnant les sources du problème, puis en
 les propageant dans le milieu selon ses propriétés et la statistique donnée par
 les lois d'extinctions (Beer-Lambert) et de diffusions (fonction de phase).
 L'extension de cette pratique à la sensibilité géométrique est possible en
 s'appuyant sur le modèle de sensibilité, qui décrit les sources de sensibilité
 géométrique et la phénoménologie de leur transport dans le milieu.  La démarche
 suivie dans ce document est la suivante:
 \begin{itemize}
   \item le problème de sensibilité géométrique est posé;
   \item les sources de sensibilité sont identifiées, elles dépendent de
   quantité décrites par d'autres physiques ({\ie} la luminance);
   \item les couplages par les sources sont exposés;
   \item un algorithme de Monte Carlo analogue est utilisé pour résoudre le
   problème couplé (via la \textit{double randomization});
   \item la mise en {\oe}uvre de l'algorithme est explicite tout au long du
   document.
 \end{itemize}
 
 L'exemple présenté dans ce document est illustratif et sans aucune ambition de
 généralité. Il montre une application du modèle de sensibilité et une
 utilisation de la méthode de Monte Carlo qui lui sont spécifiques. Cela se
 traduit dans les choix de notations des variables (indice $h$ pour paroi du
 haut, $d$ pour paroi de droite {\etc}), dans la simplification des équations du
 modèle, et enfin dans l'écriture de l'algorithme de Monte Carlo. En effet,
 contrairement à une pratique plus conventionnelle, l'ensemble des données qui
 décrivent le système (la configuration géométrique et ses propriétés physiques)
 ne sont pas séparées de la procédure d'échantillonnage des chemins. Autrement
 dit, l'ensemble des chemins seront suivis relativement aux parois de la scène
 et à leurs propriétés physiques.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Description du problème
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Description du problème}
 \label{sec:probleme}
 
 L'observable radiative de notre problème est le flux $\varphi$ perçu par le
 récepteur et s'exprime comme ci-dessous:
 \begin{equation}
 \varphi = \int_{A_r} dS \int_{\mathcal{H}^-} d\vec{\omega}
 \ (\vec{\omega} \cdot \vec{n}) L(\vec{x},\vec{\omega},\PI)
 \label{eq:flux}
 \end{equation}
 avec $\PI$ le paramètre géométrique, $A_r$ la surface du récepteur et
 $\mathcal{H}^-$ l'hémisphère orientée par $\vec{n} = - \vec{e}_z$
 (figure~\ref{fig:configuration}).
 
 \begin{figure}[h!]
   \centering
   \begin{tikzpicture}
 
     % Paroi avant
     \draw (0,0) -- (4,0) -- (4,3) -- (0,3) -- cycle;
 
     % Paroi de droite
     \draw[pattern=dots] (4,0) -- (5,1) -- (5,4) --
       node[below=35pt, preaction={fill, white}]{$A_d$} (4,3);
 
     % Paroi arrière
     \draw[dashed] (1,1) -- (5,1);
     \draw[dashed] (1,4) -- (4.5,4);
     \draw (4.5,4) -- (5,4);
 
     % Paroi de gauche
     \draw[dashed] (0,0) -- (1,1) -- (1,3);
     \draw[fill]  (1,3) -- (1,3.5);
     \draw[dashed] (1,3.5) -- (1,4) -- (0.5,3.5);
     \draw (0.5,3.5) -- (0,3);
 
     % Paroi du haut, i.e. source de sensibilité
     \draw[pattern=north west lines] (0,3.5) -- (4,3.5) --
       node[left=50pt, preaction={fill, white}]{$A_h$} (5,4.5) --
       (1,4.5) -- cycle;
 
     % Récepteur
     \draw[pattern=crosshatch]
       (1.25, 0.25) -- (2.75, 0.25) -- (3.25, 0.75) -- (1.75, 0.75) -- cycle;
 
     % Dimension du récepteur
     \draw [latex-latex] (1.75, 0.9) --
       node[above, preaction={fill, white}] {$R_x$} (3.25, 0.9);
     \draw [latex-latex] (2.95, 0.25) -- node[right=2pt] {$R_y$} (3.45, 0.75);
 
     % Dimensions de la boîte
     \draw [thin, dotted] (1, 4.5) -- (1, 5);
     \draw [thin, dotted] (5, 4.5) -- (5, 5);
     \draw [thin, dotted] (0, 3.5) -- (-0.5, 3.5);
     \draw [thin, dotted] (1, 4.5) -- (0.5, 4.5);
     \draw [thin, dotted] (5, 1) -- (6.25, 1);
     \draw [thin, dotted] (5,4) -- (5.5,4);
     \draw [thin, dotted] (5,4.5) -- (6.25,4.5);
     \draw [latex-latex] (1, 5) -- node[above] {$D_x$} (5,5);
     \draw [latex-latex] (-0.5, 3.5) -- node[left=2pt] {$D_y$} (0.5, 4.5);
     \draw [latex-latex] (5.5, 1) -- node[right] {$h$} (5.5, 4);
     \draw [latex-latex] (5.5, 4) -- node[right] {$\PI$} (5.5, 4.5);
     \draw [latex-latex] (6.25, 1) -- node[right] {$D_z(\PI)$} (6.25, 4.5);
 
     % Repère
     \draw [-latex, thick] (0,0) -- node[below] {$\vec{e}_x$} (1.25,0);
     \draw [-latex, thick] (0,0) -- node[right=2pt] {$\vec{e}_y$} (0.75,0.75);
     \draw [-latex, thick] (0,0) -- node[left]  {$\vec{e}_z$} (0,1.25);
   \end{tikzpicture}
 
   \flushleft
   \caption{\textbf{La configuration géométrique} est un parallélépipède de
     dimension $D_x \times D_y \times D_z(\PI)$, avec $D_z(\PI) = h+\PI$, repéré
     dans la base $(\vec{e}_x, \vec{e}_y, \vec{e}_z)$. Le paramètre $\PI$ est le
     paramètre géométrique défini sur $\mathbb{R}^{+}$. En le modifiant, la
     position de la paroi supérieure est translatée vers le haut, ou dit
     autrement ``la boîte s'ouvre''. Les parois latérales ne dépendent pas de
     $\PI$.  Enfin, le récepteur de surface $A_r$ est positionné sur la paroi du
     bas et a pour dimension $R_x \times R_y$.}
   \label{fig:configuration}
 
 \end{figure}
 
 \paragraph{La configuration radiative}
 Les parois du parallélépipède sont toutes noires à l'exception de la paroi du
 haut (de surface $A_h = D_x \times D_y$) qui est froide, spéculaire, et dont le
 coefficient de réflexion $\rho$ est donné par:
 \begin{equation}
 \rho(\vec{x},-\vec{\omega}) = 0.25 \Bigl[ 1- \cos \left(2 \pi \frac{x}{Dx}
 \right) \Bigr] \Bigl[ 1- \cos \left(2 \pi \frac{y}{Dy} \right) \Bigr]
 \label{eq:rho}
 \end{equation}
 avec les constantes $D_x$ et $D_y$ décrites en figure~\ref{fig:configuration}.
 Nous ne discutons pas ici du profil spécifique de $\rho$ mais soulignons
 néanmoins qu'il simplifie le problème et permet ainsi de resserrer le présent
 document sur le seul développement d'un algorithme analogue de calcul de
 sensibilité.
 
 En ce qui concerne les sources, seule la paroi noire de droite est émettrice.
 Sa surface vaut $A_d = D_y \times h$ et la condition à la limite en luminance
 correspondante est décrite par:
 \begin{equation}
 L(\vec{x},\vec{\omega},\PI) = S_b(\vec{x}) \quad \quad \quad \vec{x} \in A_d \ ; \
 \vec{\omega} \cdot \vec{n}_d > 0
 \label{eq:cl_rad}
 \end{equation}
 avec $\vec{n}_d$ la normale de la paroi de droite et $S_b$ la source radiative
 de surface qui correspond ici à son émission thermique:
 \begin{equation}
 S_b(\vec{x}) = L^{eq}(T) \Bigl[ 1- \cos \left(2 \pi \frac{x}{Dx} \right) \Bigr]
 \Bigl[ 1- \cos \left(2 \pi \frac{y}{Dy} \right) \Bigr]
 \label{eq:S_b}
 \end{equation}
 Enfin, le milieu englobant notre ``boîte'' est considéré comme étant froid et
 transparent.
 
 \paragraph{La sensibilité géométrique du flux}
 On s'intéresse à la sensibilité du flux (équation \ref{eq:flux}) par rapport à
 la variation de la position de la paroi spéculaire:
 \begin{equation}
 \partial_{\PI}\varphi = \int_{A_r} dS \int_{\mathcal{H}^-} d\vec{\omega}
 \ (\vec{\omega} \cdot \vec{n}) \partial_{\PI} L(\vec{x},\vec{\omega},\PI)
 \label{eq:sensib_flux}
 \end{equation}
 Pour évaluer le flux $\varphi$ nous avons besoin de connaître la luminance
 $L(\vec{x},\vec{\omega},\PI)$, dans toutes les directions entrantes et en tout
 point du récepteur. De façon similaire, pour évaluer $\partial_{\PI} \varphi$
 nous avons besoin de connaître la sensibilité géométrique $\partial_{\PI}
 L(\vec{x},\vec{\omega},\PI)$, dans toutes les directions entrantes et en tout
 point du récepteur.
 
 La suite du document se concentre sur l'évaluation de cette sensibilité en
 commençant par énnoncer le modèle physique qui décrit les sources et le
 transport de la sensibilité géométrique (section~\ref{sec:modele_sensib}).
 Nous développons alors un algorithme Monte Carlo pour résoudre le problème que
 nous venons de poser en suivant la propagation des sources de sensibilité via
 l'échantillonnage de chemins qui partent directement de ces sources, en
 l'occurrence ici la seule paroi du haut (section~\ref{sec:monte_carlo}).
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Le modèle de sensibilité
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Modèle de sensibilité géométrique}
 \label{sec:modele_sensib}
 
 La sensibilité géométrique de la luminance est définie telle que:
 \begin{equation}
 s(\vec{x},\vec{\omega},\PI) = \partial_{3} L(\vec{x},\vec{\omega},\PI) =
 \partial_{\PI} L(\vec{x},\vec{\omega},\PI)
 \end{equation}
 Elle est considérée comme une quantité physique à part entière dont la
 phénoménologie est décrite par une équation de transfert radiatif dans le
 domaine et par des contraintes aux frontières (conditions aux limites) sur la
 sensibilité entrante dans le domaine.
 
 \paragraph{Équation de transport} Dans \cite{papier_sensib} l'équation de la
 sensibilité dans le milieu est donnée en toute généralité. Dans notre exemple
 le milieu est transparent ($k_a = k_s = 0$) et peut donc se résumer à
 l'équation \ref{eq:ETR-S}, avec $s = s(\vec{x},\vec{\omega},\PI)$:
 \begin{equation}
  \vec{w} \cdot \vec{\nabla} s = 0
 \label{eq:ETR-S}
 \end{equation}
 
 \paragraph{Les conditions aux limites} Toujours dans \cite{papier_sensib} la
 condition à la limite de la sensibilité géométrique est donnée en toute
 généralité et dans les cas spécifiques des parois noires, parois spéculaires et
 parois diffuses. Dans notre configuration, seule la paroi supérieure du
 parallélépipède est paramétrée par $\PI$. Elle est donc la seule source de
 sensibilité géométrique:
 \begin{equation}
 \begin{aligned}
 s(\vec{x},\vec{\omega},\PI) & = 0 \quad \quad \quad \vec{x} \notin A_h \\
 s(\vec{x},\vec{\omega},\PI) & = S_{b,\PI} + \rho(\vec{x},-\vec{\omega})
 s(\vec{x},\vec{\omega}_{spec},\PI) \quad \quad \quad \vec{x} \in A_h \\
 \end{aligned}
 \end{equation}
 avec $S_{b,\PI}$ la source de sensibilité et $\rho(\vec{x},-\vec{\omega})
 s(\vec{x},\vec{\omega}_{spec},\PI)$ la réflection de la sensibilité incidente à
 la paroi dans la direction spéculaire. Dans notre exemple, le milieu est
 transparent et toutes les autres conditions aux limites de sensibilité sont
 nulles. Il n'y a donc pas de sensibilité géométrique incidente à la paroi du
 haut spéculaire.
 %On se contente donc d'ignorer la sensibilité incidente à la paroi. On
 %rappelle par ailleurs que la paroi du haut est spéculaire.
 La source de sensibilité $S_{b,\PI}$ est donc définie comme ci-dessous en
 renvoyant le lecteur à l'annexe \ref{ann:cl_sensib} pour les développements qui
 mènent à cette expression:
 \begin{equation}
 \begin{aligned}
 S_{b,\PI} = & - \beta_{\vec{\chi},h} [\partial_{1,\vec{u}_h} \
 \rho(\vec{x},-\vec{\omega})] L(\vec{x},\vec{\omega}_{spec},\PI) \quad \quad
 \quad \quad \quad \vec{x} \in A_h \ ; \ \vec{\omega} \cdot \vec{n}_h > 0\\
 & + \rho(\vec{x},-\vec{\omega})
 \partial_{1,\vec{\chi}}L(\vec{x},\vec{\omega}_{spec},\PI)\\
 & - \rho(\vec{x},-\vec{\omega}) \beta_{\vec{\chi},h} \partial_{1,\vec{u}_h}
 L(\vec{x},\vec{\omega}_{spec},\PI)
 \end{aligned}
 \label{eq:clsensib}
 \end{equation}
 avec $\vec{\omega}_{spec} = \vec{\omega} - 2 (\vec{\omega} \cdot \vec{n}_h)
 \vec{n}_h$ et $\beta_{\vec{\chi},h}$ le coefficient issu de la décomposition de
 $\vec{\chi}$ en deux vecteurs, l'un orienté par $\vec{\omega}$ et l'autre
 orienté par un vecteur $\vec{u}_h$ tangent à la paroi du haut (voir annexe
 \ref{ann:proj}). Enfin $\beta_{\vec{\chi},h}$ est la norme du vecteur
 $\vec{\chi}$ projeté sur $\vec{u}_h$. La dérivée spatiale
 $\partial_{1,\vec{\gamma}} f(\vec{x},\vec{\omega}) = \vec{\gamma} \cdot
 \vec{\nabla}_{\vec{x}} f(\vec{x},\vec{\omega})$ est la dérivée directionnelle
 dans la direction $\vec{\gamma}$.
 
 \paragraph{La source de sensibilité} est dans notre cas une source de surface,
 émise par la paroi du haut et donnée par la condition à la limite décrite par
 l'expression~\ref{eq:clsensib}. Dans cette équation on note que la condition à
 la limite de sensibilité dépend de:
 \begin{itemize}
   \item la luminance incidente à la paroi dans la direction de transport
   spéculaire;
   \item la dérivée spatiale de la luminance dans la direction de dérivation
   $\vec{u}$, incidente à la paroi dans la direction de transport spéculaire;
   \item et la dérivée spatiale de la luminance dans la direction de dérivation
   $\vec{\chi}$, incidente à la paroi dans la direction de transport
   spéculaire.
 \end{itemize}
 En résumé, la source de sensibilité émise par la paroi spéculaire est le
 résultat du couplage entre le modèle de sensibilité, le modèle de transfert
 radiatif et le modèle de dérivée spatiale. La dérivée spatiale de la luminance
 est simplement considérée comme une quantité physique, au même titre que la
 sensibilité géométrique (voir~\cite{papier_sensib} pour la description de son
 modèle). Résoudre notre problème de sensibilité géométrique revient donc à
 résoudre un problème de transport couplé qui dépend à la fois des source
 radiatives (à travers $L(\vec{x},\vec{\omega}_{spec},\PI)$), des sources de
 dérivées spatiales dans la direction $\vec{u}$ (à travers $\partial_{1,\vec{u}}
 L(\vec{x},\vec{\omega}_{spec},\PI)$) et des sources de dérivées spatiales dans
 la direction $\vec{\chi}$ (à travers $\partial_{1,\vec{\chi}}
 L(\vec{x},\vec{\omega}_{spec},\PI)$).
 
 \paragraph{Les sources du problème couplé}
 La seule source radiative de notre configuration est donnée en
 section~\ref{sec:probleme} par l'équation \ref{eq:cl_rad}. Elle correspond à
 l'émission thermique $S_b$ de la paroi de droite de surface $A_d$. Pour les
 sources de dérivée spatiale, le modèle de dérivée spatiale élaboré dans
 \cite{papier_sensib} autorise des sources volumiques, des sources de surfaces
 et des sources locales situées sur les arrêtes d'une géométrie triangulées.
 Dans notre exemple ce modèle se simplifie de sorte que ces sources se résument
 aux seules sources émises par la surface du haut $A_h$ et la surface de droite
 $A_d$ (figure~\ref{fig:configuration}).
 
 Pour la dérivée spatiale dans la direction $\vec{u}_h$, la source de la paroi
 du haut est donnée par la condition à la limite:
 \begin{equation}
 \partial_{1,\vec{u}_h} L(\vec{x},\vec{\omega},\PI) = \beta_{\vec{u}_h,h} [
 \partial_{1,\vec{u}_s} \rho(\vec{x},-\vec{\omega}) ]
 L(\vec{x},\vec{\omega}_{spec},\PI) \quad \quad \quad \vec{x} \in A_h \ ; \
 \vec{\omega} \cdot \vec{n}_h > 0
 \label{eq:cl_duL_haut}
 \end{equation}
 et la source de la paroi de droite est donnée par la condition à la limite:
 \begin{equation}
 \partial_{1,\vec{u}_h} L(\vec{x},\vec{\omega},\PI) = \beta_{\vec{u}_h,d}
 \partial_{1,\vec{u}_{hd}} S_b(\vec{x}) \quad \quad \quad \vec{x} \in A_d \ ;
 \ \vec{\omega} \cdot \vec{n}_d > 0
 \label{eq:cl_duL_droite}
 \end{equation}
 
 Pour la dérivée spatiale dans la direction $\vec{\chi}$, la source de la paroi
 du haut est donnée par la condition à la limite:
 \begin{equation}
 \partial_{1,\vec{\chi}} L(\vec{x},\vec{\omega},\PI) = \beta_{\vec{\chi},h} [
 \partial_{1,\vec{u}_h} \rho(\vec{x},-\vec{\omega}) ]
 L(\vec{x},\vec{\omega}_{spec},\PI) \quad \quad \quad \vec{x} \in A_h \ ; \
 \vec{\omega} \cdot \vec{n}_h > 0
 \label{eq:cl_dchiL_haut}
 \end{equation}
 et la source de la paroi de droite est donnée par la condition à la limite:
 \begin{equation}
 \partial_{1,\vec{\chi}} L(\vec{x},\vec{\omega},\PI) = \beta_{\vec{\chi},d}
 \partial_{1,\vec{u}_{\chi,d}} S_b(\vec{x}) \quad \quad \quad \vec{x} \in A_d \
 ; \ \vec{\omega} \cdot \vec{n}_d > 0
 \label{eq:cl_dchiL_droite}
 \end{equation}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Algorithme Direct
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Résolution par Monte Carlo}
 \label{sec:monte_carlo}
 
 Dans cette section nous écrivons un algorithme de Monte Carlo pour résoudre le
 problème décrit en section~\ref{sec:probleme} à partir du modèle de sensibilité
 géométrique décrit en section~\ref{sec:modele_sensib}. Notre algorithme est
 construit de manière analogue avec pour support la mise en {\oe}uvre pratique
 de sa [[<<Fonction de réalisation>>]] en langage C~\cite{C_language}. Pour
 cela, nous échantillonnerons des chemins qui démarrent directement de la source
 de sensibilité, à savoir la paroi du haut, et qui la propage jusqu'au récepteur
 positionné sur la paroi du bas. La mise en {\oe}uvre ici proposée est néanmoins
 particulière dans le sens où chaque chemin du problème couplé est d'abord
 échantillonné et conservé en totalité (section~\ref{subsec:chemin}). Son poids
 n'est calculé qu'a posteriori à partir du chemin ainsi construit
 (section~\ref{subsec:poids}). Cette proposition diffère d'un algorithme Monte
 Carlo plus conventionnel où la position et la direction courante du chemin ne
 sont que des données locales à chaque étape de sa construction; son poids étant
 mis à jour si nécessaire à chacune de ces étapes. Ce choix singulier est dicté
 par la volonté de garder à tout instant une vue d'ensemble du problème pour
 s'éviter un travail d'abstraction et ainsi faciliter l'écriture analogue de
 notre algorithme. Une démarche adaptée en raison de la configuration simplifiée
 à laquelle nous nous intéressons.
 
 <<Fonction de réalisation>>=
 static res_T
 realisation
   (struct ssp_rng* rng,
    const struct sgs_scene* scene,
    double* w)
 {
   <<Données locales à la fonction de réalisation>>
   res_T res = RES_OK;
 
   <<Initialiser le poids>>
 
   <<Échantillonner un chemin du problème couplé>>
   <<Calcul du poids>>
 
 exit:
   <<Nettoyer les données locales de la fonction>>
   <<Renvoyer le poids>>
   return res;
 }
 @
 
 Notre fonction de réalisation prend en entrée un générateur de nombres
 aléatoires ([[rng]]) et un pointeur vers les données du système ([[scene]]).
 Dans la variable [[w]] sera renvoyé le poids de la sensibilité a $\PI$.
 
 \subsection{Le chemin}
 \label{subsec:chemin}
 
 Pour construire un chemin du problème couplé on échantillonne d'abord un chemin
 de sensibilité partant d'une position quelconque sur la surface de la source de
 sensibilité $A_h$. Ce chemin est alors complété par l'échantillonnage d'un
 chemin de dérivée spatiale dont la couplage avec le chemin de sensibilité forme
 un chemin du problème couplé.
 
 <<Échantillonner un chemin du problème couplé>>=
 <<Échantillonner une position sur la source de sensibilité>>
 <<Échantillonner un chemin de sensibilité>>
 <<Échantillonner un chemin de dérivée spatiale>>
 @
 
 Comme point de départ du chemin du problème couplé, on commence donc par
 échantillonner uniformément un point sur la surface émettrice de sensibilité à
 l'aide de la fonction [[sgs_geometry_sample_sensitivity_source]], et on stocke
 dans les variables [[pos_h]] et [[normal_h]] sa position et la normale
 correspondante. On récupère également dans [[surf_A_h]] l'identifiant de la
 surface que l'on vient d'échantillonner, dans notre cas la surface supérieure
 $A_h$ identifié dans le code par la constante [[SGS_SURFACE_Z_MAX]] (voir
 figure~\ref{fig:configuration}).
 
 <<Échantillonner une position sur la source de sensibilité>>=
 /* Échantillonner uniformément une position sur la source de sensibilité */
 sgs_geometry_sample_sensitivity_source(scene->geom, rng, &frag);
 d3_set(pos_h, frag.position);
 d3_set(normal_h, frag.normal);
 surf_A_h = frag.surface; /* surf_A_h == SGS_SURFACE_Z_MAX */
 @
 
 On rappelle que les sources de sensibilité proviennent des parois perturbées
 par une modification du paramètre $\PI$. En toute hypothèse, toute source de
 sensibilité incidente à $A_h$ serait réfléchie de façon spéculaire. Or, dans
 notre cas, nous n'avons qu'une seule source de sensibilité, la surface $A_h$
 elle-même, et par conséquent nous n'avons pas à tenir compte de ces
 sensibilités réfléchies. Nous n'échantillons donc que le seul chemin qui
 propage l'émission de sensibilité par $A_h$. Ce chemin sera notre chemin de
 sensibilité. Pour cela, il suffit d'échantillonner une direction d'émission
 lambertienne [[dir_emit_h]] autour de la normale [[normal_h]] de la surface
 $A_h$, et de lancer un rayon dans cette direction. Nous stockons alors dans
 [[hit0]] l'intersection de ce rayon avec la géométrie de la scène.
 
 <<Échantillonner un chemin de sensibilité>>=
 /* Échantillonner une direction d'émission de sensibilité  */
 ssp_ran_hemisphere_cos(rng, normal_h, dir_emit_h, NULL);
 /* Lancer le rayon qui propage l'émission de sensibilité */
 TRACE_RAY(pos_h, dir_emit_h, surf_A_h, &hit0);
 @
 
 La source de sensibilité est donnée dans la condition à la limite décrite par
 l'équation~\ref{eq:clsensib}. Elle dépend de la dérivée spatiale selon
 $\vec{\chi}$ incidente dans la direction spéculaire $\vec{\omega}_{spec}$ et de
 la dérivée spatiale selon $\vec{u}$ incidente à la même direction spéculaire.
 Ces contributions à l'émission de sensibilité ne sont pas connues et sont ici
 échantillonnées par \textit{double randomization}. Comme ces deux dérivées
 spatiales sont incidentes à la même direction $\vec{\omega}_{spec}$ nous
 pouvons nous contenter de ne suivre qu'un seul chemin dans cette direction. Ce
 chemin sera notre chemin de dérivée spatiale. Pour échantillonner ce chemin
 nous calculons d'abord $\vec{\omega}_{spec}$ ([[dir_spec_h]]) par réflexion
 spéculaire de la direction d'émission $\vec{\omega}$ ([[dir_emit_h]]) avant de
 lancer un rayon dans cette direction jusqu'à l'intersection [[hit1]] avec une
 surface .
 
 <<Échantillonner un chemin de dérivée spatiale>>=
 /* Calculer la direction spéculaire à dir_emit_h */
 reflect(dir_spec_h, dir_emit_h, normal_h);
 /* Tracer le rayon spéculaire partant de surf_A_h */
 TRACE_RAY(pos_h, dir_spec_h, surf_A_h, &hit1);
 @
 
 \subsection{Le poids}
 \label{subsec:poids}
 
 Dans notre problème, un chemin du problème couplé est composé d'un chemin de
 sensibilité et d'un chemin de dérivée spatiale, chacun d'entre eux se résumant
 à un segment dont l'origine commune se situe sur la paroi du haut, source de
 sensibilité. Or, on peut dès à présent déterminer qu'un chemin couplé aura une
 contribution nulle si le chemin de sensibilité n'atteint pas le récepteur ou si
 le chemin de dérivé spatiale n'atteint pas la source radiative, à savoir la
 paroi de droite.
 
 <<Initialiser le poids>>=
 sensib = 0;
 @
 <<Échantillonner un chemin de sensibilité>>=
 if(!hit_receiver(scene, pos_h, dir_emit_h, &hit0)) {
   goto exit;
 }
 @
 <<Échantillonner un chemin de dérivée spatiale>>=
 if(!hit_source(scene, pos_h, dir_spec_h, &hit1)) {
   goto exit;
 }
 @
 
 En conséquence nous pouvons assumer dans la suite de la fonction que nous
 n'aurons à calculer le poids que des seuls chemin couplés dont la contribution
 est non nulle. Dès lors, [[hit1]] représente une intersection sur la source
 radiative. On stocke dans [[normal_d]] la normale de la paroi correspondante
 dont on aura besoin pour le calcul du poids. De même, on initialise la variable
 [[dir_emit_d]] à $-\vec{\omega}_s$, cette direction nous sera également utile
 pour évaluer la source de dérivée spatiale.
 
 <<Échantillonner un chemin de dérivée spatiale>>=
 d3_normalize(normal_d, hit1.normal);
 d3_minus(dir_emit_d, dir_spec_h);
 @
 
 Dans notre problème couplé, la contribution du chemin, {\ie} le poids Monte
 Carlo, va s'exprimer à travers la condition à la limite de sensibilité
 (équation \ref{eq:clsensib}) et des sources de chacun de ses couplages.
 
 <<Calcul du poids>>=
 <<Décomposition du vecteur de déformation $\vec{\chi}$>>
 
 <<Calcul de la dérivée surfacique de $\rho$>>
 <<Calcul des sources de dérivées spatiales>>
 
 <<Calculer le poids de sensibilité>>
 @
 
 La décomposition du vecteur de déformation $\vec{\chi}$ permet d'obtenir le
 vecteur tangent $\vec{u}$ nécessaire dans l'expression de la source de
 sensibilité et de ses dérivées surfaciques (équation \ref{eq:clsensib}).
 
 <<Décomposition du vecteur de déformation $\vec{\chi}$>>=
 decomposition(chi, normal_h, dir_emit_h, &proj_chi_h);
 @
 
 Étant donné que le coefficient de réflection n'est défini qu'en frontière, à
 savoir sur un plan en deux dimensions, calculer la dérivée surfacique de $\rho$
 revient à travailler dans le plan. Et cette dérivée surfacique est le
 produit scalaire entre le gradient surfacique de $\rho$ et la direction de
 dérivation $\vec{u}$ transformée dans ce plan ([[u_2d]]).
 
 <<Calcul de la dérivée surfacique de $\rho$>>=
 <<Récupérer le gradient surfacique de $\rho$>>
 
 /* Transformer u dans le plan XY */
 u_2d[0] = proj_chi_h.u[X];
 u_2d[1] = proj_chi_h.u[Y];
 
 /* Calculer la dérivée surfacique de rho */
 d_rho = d2_dot(grad_rho_2d, u_2d);
 @
 
 Pour récupérer $\rho$ et son gradient, il nous suffit de transformer dans le
 plan la position d'émission [[pos_h]], et d'interroger les données associées à
 la position ainsi transformée ([[pos_h_2d]]).
 
 <<Récupérer le gradient surfacique de $\rho$>>=
 /* Transformer pos_h dans le plan XY */
 pos_h_2d[0] = pos_h[X];
 pos_h_2d[1] = pos_h[Y];
 
 rho = get_rho(scene, pos_h_2d);
 get_grad_rho(scene, pos_h_2d, grad_rho_2d);
 @
 
 On rappelle que la sensibilité est couplée à deux dérivées spatiales (selon
 $\vec{\chi}$ et $\vec{u}$) dont les sources sont données par les équations
 \ref{eq:cl_duL_haut}, \ref{eq:cl_duL_droite}, \ref{eq:cl_dchiL_haut} et
 \ref{eq:cl_dchiL_droite}.
 
 <<Calcul des sources de dérivées spatiales>>=
 <<Décomposition du vecteur $\vec{\chi}$>>
 <<Décomposition du vecteur $\vec{u}$>>
 <<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_e$>>
 <<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_{cs}$>>
 @
 
 Sur la paroi de droite, la décomposition du vecteur de déformation $\vec{\chi}$
 permet d'obtenir le vecteur tangent $\vec{u}$ nécessaire dans l'expression de
 la source de de la dérivée spatiale dans la direction $\vec{\chi}$ et de sa
 dérivée surfacique (équation \ref{eq:cl_dchiL_droite}).
 
 <<Décomposition du vecteur $\vec{\chi}$>>=
 decomposition(chi, normal_d, dir_emit_d, &proj_chi_d);
 @
 
 De la même façon, la décomposition du vecteur $\vec{u}$ permet d'obtenir le
 vecteur tangent $\vec{u}_e$ nécessaire dans l'expression de la source de la
 dérivée spatiale dans la direction $\vec{u}$ et de sa dérivée surfacique
 (equation \ref{eq:cl_duL_droite}).
 
 <<Décomposition du vecteur $\vec{u}$>>=
 decomposition(proj_chi_h.u, normal_d, dir_emit_d, &proj_uh_d);
 @
 
 La dérivée surfacique de $S_b$ dans la direction $\vec{u}_{e}$ ([[dSb_uhd]])
 est le produit scalaire entre le gradient surfacique de $S_b$ et la direction
 de dérivation $\vec{u}_{e}$ transformée dans le plan de la paroi de droite
 ([[u_hd_2d]]). Mais pour effectuer ce calcul nous devons au préalable récupérer
 le gradient de $S_b$ ([[grad_Sb_2d]]). Pour cela il nous suffit là encore de
 transformer dans le plan de la paroi de droite la position d'émission
 [[pos_d]], et d'interroger les données associées à la position ainsi
 transformée ([[pos_d_2d]]). Et nous en profitons au passage pour récupérer
 $S_b$ qui nous sera utile pour le calcul du poids. Pour rappel, dans notre
 configuration les deux dérivées spatiales $\partial_{\vec{\chi}} I$ et
 $\partial_{\vec{u}} I$ partagent une même position d'émission ([[pos_d]]).
 
 <<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_e$>>=
 /* Calculer la position sur l'émetteur */
 pos_d[X] = pos_h[X] + dir_spec_h[X]*hit1.distance;
 pos_d[Y] = pos_h[Y] + dir_spec_h[Y]*hit1.distance;
 pos_d[Z] = pos_h[Z] + dir_spec_h[Z]*hit1.distance;
 
 /* Transformer pos_d dans le plan YZ */
 pos_d_2d[0] = pos_d[Y];
 pos_d_2d[1] = pos_d[Z];
 
 /* Récupérer le gradient surfacique de Sb */
 Sb = get_Sb(scene, pos_d_2d);
 get_grad_Sb(scene, pos_d_2d, grad_Sb_2d);
 
 /* Transformer u_e dans le plan YZ */
 u_hd_2d[0] = proj_uh_d.u[Y];
 u_hd_2d[1] = proj_uh_d.u[Z];
 
 /* Dérivée surfacique de Sb dans la direction u_e */
 dSb_uhd = d2_dot(grad_Sb_2d, u_hd_2d);
 @
 
 Ne reste plus qu'à calculer la dérivée surfacique de $S_b$ dans la direction
 $\vec{u}_{cs}$ ([[dSb_uchid]]) comme étant le produit scalaire entre le
 gradient surfacique de $S_b$ et la direction $\vec{u}_{cs}$ transformée dans le
 plan de la paroi de droite ([[u_chid_2d]]).
 
 <<Calcul de la dérivée surfacique de $Sb$ dans la direction $\vec{u}_{cs}$>>=
 /* Transformer u_cs dans le plan YZ */
 u_chid_2d[0] = proj_chi_d.u[Y];
 u_chid_2d[1] = proj_chi_d.u[Z];
 
 /* Dérivée surfacique de Sb dans la direction u_cs */
 dSb_uchid = d2_dot(grad_Sb_2d, u_chid_2d);
 @
 
 On est alors en mesure d'évaluer le poids de sensibilité comme la contribution
 portée par le chemin du problème couplé, soit les sources de dérivées spatiales
 propagées jusqu'à la paroi du haut puis intégrées dans les sources de
 sensibilité propagées ensuite vers le récepteur.  Cette contribution est
 ensuite multipliée par la surface $A_h$ et l'angle $\pi$ ([[PI]]) étant donné
 l'échantillonnage des densités de probabilités de la position sur la surface et
 de la direction dans l'hémisphère sortante de la paroi.
 
 <<Calculer le poids de sensibilité>>=
 /* Calcul de la contribution du chemin couplé */
 Sb_sensib =
   - proj_chi_h.beta * d_rho * Sb
   - rho * proj_chi_h.beta * proj_uh_d.beta * dSb_uhd
   + rho * proj_chi_d.beta * dSb_uchid;
 
 /* Poids de sensibilité */
 sensib = Sb_sensib * PI * get_Sr(scene);
 @
 
 <<Renvoyer le poids>>=
 w[0] = sensib;
 @
 
 \subsection{Variables locales et macro}
 
 Lors de l'échantillonnage des chemins (section~\ref{subsec:chemin}) nous nous
 sommes appuyé sur la fonction [[TRACE_RAY]] nous permettant de lancer un rayon
 dans la scène. Cette fonction est en réalité un macro, locale à la fonction,
 qui encapsule l'appel à la fonction qui lance un rayon. En plus de l'origine et
 de la direction du rayon, cette fonction nécessite en entrée une plage des
 distances d'intersection possible ([[range]]), dans notre cas toujours définie
 à $[0, \infty]$. À noter également le paramètre d'entrée [[StartFrom]] qui
 stocke le triangle sur lequel se trouve l'origine du rayon, une donnée d'entrée
 utilisée pour éviter une auto-intersection, c'est à dire l'intersection du
 rayon avec le triangle dont il est issu.
 
 <<Données locales à la fonction de réalisation>>=
 /* Macro utilisée comme sucre syntaxique */
 #define TRACE_RAY(Org, Dir, StartFrom, Hit) {                          \
   double range[2];                                                     \
   range[0] = 0;                                                        \
   range[1] = INF;                                                      \
   sgs_geometry_trace_ray                                               \
     (scene->geom, (Org), (Dir), range, (StartFrom), (Hit));\
 } (void)0
 @
 
 <<Nettoyer les données locales de la fonction>>=
 #undef TRACE_RAY
 @
 
 \paragraph{}
 Enfin, nous déclarons ci-après l'ensemble des variables locales nécessaires à
 notre [[<<Fonction de réalisation>>]]:
 
 <<Données locales à la fonction de réalisation>>=
 /* Variables de la source de sensibilité */
 double dir_emit_h[3];
 double pos_h[3];
 double normal_h[3];
 double dir_spec_h[3];
 double pos_h_2d[2];
 enum sgs_surface_type surf_A_h;
 struct sgs_fragment frag; /* Position échantillonnée */
 
 /* Variables de la source radiative */
 double pos_d[3];
 double normal_d[3];
 double dir_emit_d[3];
 double pos_d_2d[2];
 
 /* Propriétés radiatives des surfaces et leur dérivée */
 double rho;
 double d_rho;
 double grad_rho_2d[2];
 double Sb;
 double dSb_uhd;
 double dSb_uchid;
 double grad_Sb_2d[2];
 
 /* Vecteurs des déformations et variables de projection */
 const double chi[3] = {0, 0, 1};
 double u_2d[2];
 double u_hd_2d[2];
 double u_chid_2d[2];
 struct projection proj_chi_h;
 struct projection proj_chi_d;
 struct projection proj_uh_d;
 
 /* Pour le calcul du poids */
 double Sb_sensib;
 double sensib;
 
 /* Intersections avec la géométrie */
 struct sgs_hit hit0;
 struct sgs_hit hit1;
 @
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Résultats
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Résultats}
 \label{sec:results}
 
 Au delà du simple calcul de sensibilité nous proposons ici une vérification des
 résultats par le calcul de différences finies. La différentiation est effectuée
 à partir des estimations du flux $\varphi(\PI)$ reçu par le capteur pour
 différentes valeurs du paramètre géométrique $\PI$, soit pour différentes
 hauteurs de la paroi spéculaire. Le calcul du poids associé au calcul du flux
 est décrit en annexe \ref{annex:flux}. La figure \ref{fig:resultats} présente les
 estimations de la sensibilité du flux et des différences finies correspondantes.
 Les paramètres des simulations Monte Carlo et des calculs en différences finies,
 et plus généralement le code source des scripts à l'origine de ces résultats
 sont données en annexe~\ref{annex:scripts}.
 
 \begin{figure}[h!]
   \centering
   \begin{tikzpicture}
     \begin{axis}[
       xlabel=$\frac{\PI}{h}$,
       ylabel=$\partial_{\PI} \hat{\varphi}$,
       width=0.7\linewidth,
       legend style={at={(0.95,0.3)}}
     ]
       \addplot[
         color=black,
         mark=square*,
         only marks,
         mark size=2pt,
         every mark/.append style={solid, fill=black},
         error bars/.cd,
         y dir=both,
         y explicit,
         error mark=|,
         error mark options={mark size=1pt},
         error bar style={
            line width=0.6pt,
            color=black}
       ]
       table[x=pi_over_h, y=sen_mc, y error=err_mc, col sep=space]{results.fd};
       \addlegendentry{Monte Carlo}
 
       \addplot[
         color=black,
         mark=*,
         only marks,
         mark size=1.5pt,
         every mark/.append style={solid, fill=white},
         error bars/.cd,
         y dir=both,
         y explicit,
         error mark=|,
         error mark options={mark size=1pt},
         error bar style={
            line width=0.6pt,
            color=black}
       ]
       table[x=pi_over_h, y=sen_fd, y error=err_fd, col sep=space]{results.fd};
       \addlegendentry{Différences Finies}
     \end{axis}
 
   \end{tikzpicture}
   \flushleft
   \caption{Estimations des sensibilités du flux par Monte Carlo et comparaison
     avec les différences finies du flux. Les résultats sont adimensionnalisés et
     affichés en fonction de l'ouverture du parallélépipède, qui est paramétrée
     par $\PI$. La sensibilité du flux $\partial_{\PI} \hat{\varphi} =
     \frac{\partial_{\PI} \varphi}{\varphi_{max}} h$ avec $\varphi_{max} =
     \varphi_{spec}(\PI=0)$. Dans cette expression, $\varphi_{spec}$ correspond
     uniquement à la partie du flux qui arrive au récepteur après avoir été
     réfléchie sur la paroi spéculaire (voir annexe~\ref{annex:flux}). Le nombre
     d'échantillonnage du poids de Monte Carlo nécessaire à la reproduction de ce
     résultat est $10^{8}$.}
   \label{fig:resultats}
 \end{figure}
 
 \appendix
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Annexe CL de sensibilité
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Détails de la condition à la limite de sensibilité}
 \label{ann:cl_sensib}
 
 Nous récupérons ici la condition à la limite pour une paroi spéculaire donnée
 dans \cite{papier_sensib} $s = s(\vec{x},\vec{\omega},\PI)$:
 \begin{equation}
 s = \mathcal{C}_b \lbrack s \rbrack + S_{b,\PI} \lbrack L, \partial_{1,\vec{u}}L,
 \partial_{1,\vec{\chi}} L, \partial_{2,\vec{\gamma}_t}L \rbrack \quad \quad \quad
 \vec{x} \in \partial G(\PI) ; \vec{\omega} \cdot \vec{n} > 0
 \label{eq:cl_sensib_gen}
 \end{equation}
 avec $S_{b,\PI}$ la source surfacique de sensibilité:
 \begin{equation}
 \begin{aligned}
 S_{b,\PI} = & - \alpha (\mathcal{C}[L] + S) \\
 & - \beta \partial_{1,\vec{u}} S_b - \partial_{2,\vec{\gamma}} S_b +
 \partial_{\PI} S_b \\
 & - \beta \partial_{1,\vec{u}} \mathcal{C}_b[L] + \partial_{\PI} \mathcal{C}_b
 [L] \\
 & - \partial_{2,\vec{\gamma}} \rho(\vec{x},-\vec{\omega}) \int_{H'}
 p_{\Omega'}(-\vec{\omega}'|\vec{x},-\vec{\omega})d\vec{\omega}' L \\
 & - \beta \mathcal{C}_b[\partial_{1,\vec{u}}L] +
 \mathcal{C}_b[\partial_{1,\vec{\chi}} L] \\
 & + 2 \mu \partial_{2,\vec{\gamma}_t} L(\vec{x},\vec{\omega}_{spec},\PI)
 \end{aligned}
 \end{equation}
 
 Dans cette équation $\mathcal{C}$ est l'opérateur collisionnel du milieu:
 \begin{equation}
 \mathcal{C}[L] = -k_a L(\vec{x},\vec{\omega},\PI) - k_s
 L(\vec{x},\vec{\omega},\PI) + k_s \int_{\mathcal{S}}
 p_{\Omega'}(\vec{\omega}'|\vec{x},\vec{\omega}) d\vec{\omega}'
 L(\vec{x},\vec{\omega}',\PI)
 \end{equation}
 La source $S$ est la source radiative du milieu (émission thermique $k_a
 L^{eq}(T)$). On trouve également $\mathcal{C}_b$ l'opérateur collisionnel de la
 surface. Appliqué à la luminance et dans le cas d'une paroi spéculaire il
 s'écrit:
 \begin{equation}
 \mathcal{C}_b[L] = \rho(\vec{x},-\vec{\omega}) \int_{H'} \delta(\vec{\omega}' -
 \vec{\omega}_{spec}) L(\vec{x},\vec{\omega}',\PI) d\vec{\omega}'
 \end{equation}
 avec $\vec{\omega}_{spec} = \vec{\omega} - 2(\vec{\omega} \cdot
 \vec{n})\vec{n}$
 
 \paragraph{Condition à la limite de notre exemple}
 Pour commencer, seule la paroi spéculaire est source de sensibilité géométrique.
 Nous voyons dans l'équation \ref{eq:cl_sensib_gen} que le terme collisionnel
 $\mathcal{C}_b[s]$ traduit la réflexion de la sensibilité incidente à la paroi
 spéculaire. Ce terme est obligatoirement nul puisque il n'y a aucune autre
 source qui pourrait émettre une sensibilité géométrique et aucune autre paroi
 réfléchissante qui pourrait réfléchir la sensibilité émise par la paroi
 spéculaire.
 
 Dans notre exemple le milieu est transparent, les termes
 $\mathcal{C}[L]$ et $S$ sont donc nuls. La paroi spéculaire est
 froide, la source surfacique $S_b$ qui dans cet exemple
 correspondrait à l'émission thermique de la paroi est donc aussi
 nulle.  L'opérateur collisionnel de la surface $\mathcal{C}_b$ est
 indépendant de $\PI$, la dérivée $\partial_{\PI} \mathcal{C}_b$ est
 donc nulle. Pour finir la déformation géométrique de la paroi
 spéculaire est une translation, l'axe de rotation $\vec{\gamma}$ est
 donc nul et toutes les dérivées angulaires n'ont plus lieux d'être
 dans la condition à la limite.
 
 La condition à la limite de sensibilité de la paroi spéculaire de la boite
 devient donc:
 \begin{equation}
 \begin{aligned}
 s = & - \beta \partial_{1,\vec{u}} \mathcal{C}_b[L] \\
 & + \mathcal{C}_b[\partial_{1,\vec{\chi}}L - \beta \partial_{1,\vec{u}} L]
 \end{aligned}
 \end{equation}
 avec
 \begin{equation}
 \beta \partial_{1,\vec{u}} \mathcal{C}_b [L] = \beta \left(\partial_{1,\vec{u}}
 \rho(\vec{x},-\vec{\omega}) \right) \int_{H'}
 \delta(\vec{\omega}'-\vec{\omega}_{spec}) L(\vec{x},\vec{\omega}',\PI)
 d\vec{\omega}'
 \end{equation}
 En prenant en compte le fait que:
 \begin{equation}
 \int_{H'} \delta(\vec{\omega}' - \vec{\omega}_{spec}) f(\vec{\omega}')
 d\vec{\omega}' = f(\vec{\omega}_{spec})
 \end{equation}
 on trouve finalement:
 \begin{equation}
 \begin{aligned}
 s(\vec{x},\vec{\omega},\PI) = & - \beta
 \left(\partial_{1,\vec{u}}\rho(\vec{x},-\vec{\omega}) \right)
 L(\vec{x},\vec{\omega}_{spec},\PI)  \\
 & + \partial_{1,\vec{\chi}}L(\vec{x},\vec{\omega}_{spec},\PI) - \beta
 \partial_{1,\vec{u}} L(\vec{x},\vec{\omega}_{spec},\PI)
 \end{aligned}
 \end{equation}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Annexe décomposition
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Décomposition du vecteur de déformation}
 \label{ann:proj}
 
 Dans le modèle de sensibilité la déformation est caractérisée par le
 vecteur de déformation $\vec{\chi}$. La condition à la limite de
 sensibilité dépend alors de la dérivée spatiale
 $\partial_{1,\vec{\chi}}L$ (équation \ref{eq:clsensib}).  Le champs
 de luminance n'étant pas connu il n'existe pas de solution
 analytique à cette dérivée. À la frontière nous avons donc choisi de
 décomposer la direction $\vec{\chi}$ ([[chi]]) en deux directions
 distinctes, une direction tangente à la frontière $\vec{u}$ ([[u]])
 et la direction de transport $\vec{\omega}$ ([[omega]]). Ainsi la
 dérivée de $L$ selon $\vec{\chi}$ devient une composition de
 dérivées de $L$ le long de la paroi (sur laquelle la luminance est
 connue) et le long de la direction de transport (retrouvant ainsi le
 terme de transport de l'ETR).
 
 La décomposition de $\vec{\chi}$ s'écrit:
 \begin{equation}
 \vec{\chi} = \alpha \vec{\omega} + \beta \vec{u}
 \label{eq:chi_decomp}
 \end{equation}
 avec $\alpha$ ([[alpha]]) et $\beta$ ([[beta]]) les coefficients issus de la
 projection de $\vec{\chi}$ sur $\vec{n}$ ([[normal]]) et $\vec{u}$.  Pour plus
 de précisions sur la base non-orthogonale utilisée pour cette décomposition,
 nous renvoyons le lecteur vers \cite{papier_sensib}.  Nous nous contentons de
 donner ici les résultats:
 \begin{equation}
 \alpha = \frac{\vec{\chi}.\vec{n}}{\vec{\omega}.\vec{n}} \ \; \quad \beta =
 \|\vec{\chi} - \alpha \vec{\omega} \| \ \; \quad \vec{u} = \frac{\vec{\chi} -
 \alpha \vec{\omega}}{\beta}
 \end{equation}
 Ce qui permet d'obtenir:
 \begin{equation}
 \partial_{1,\vec{\chi}}L = \alpha \partial_{1,\vec{\omega}} L + \beta
 \partial_{1,\vec{u}} L
 \end{equation}
 et de résoudre la dérivée en estimant $\partial_{1,\vec{u}}L$ le long de la
 surface et en utilisant l'ETR pour résoudre $\partial_{1,\vec{\omega}}L$:
 \begin{equation}
 \partial_{1,\vec{\omega}}L = \mathcal{C} \lbrack L \rbrack
 \end{equation}
 
 La fonction [[decomposition]] réalise cette décomposition et les résultats
 ([[alpha]], [[beta]] et [[u]]) sont retournés via les variables membres de la
 structure [[projection]]:
 
 <<Fonctions utilitaires>>=
 struct projection {
   double alpha;
   double beta;
   double u[3];
 };
 
 static void
 decomposition
   (const double chi[3],
    const double normal[3],
    const double omega[3],
    struct projection* projection)
 {
   ASSERT(chi && normal && omega && projection);
   ASSERT(d3_is_normalized(normal));
   ASSERT(d3_is_normalized(omega));
 
   projection->alpha = d3_dot(chi, normal) / d3_dot(omega, normal);
   projection->u[X] = chi[X] - projection->alpha*omega[X];
   projection->u[Y] = chi[Y] - projection->alpha*omega[Y];
   projection->u[Z] = chi[Z] - projection->alpha*omega[Z];
   projection->beta = d3_normalize(projection->u, projection->u);
 }
 @
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Annexe "luminance"
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Calcul de la contribution de la luminance qui dépend de $\PI$}
 \label{annex:flux}
 
 Le flux reçu par le capteur est décrit par l'équation \ref{eq:flux}. Dans notre
 problème, la luminance $L(\vec{x},\vec{\omega},\PI)$ incidente au récepteur
 dépend de deux sortes de chemins radiatifs qui traduisent le transport de la
 source émise par la paroi de droite :
 \begin{itemize}
   \item soit la source est transportée directement depuis la paroi de droite
   jusqu'au récepteur;
   \item soit la source est transportée en direction de la paroi spéculaire puis
   réfléchie jusqu'au récepteur.
 \end{itemize}
 Dans notre configuration le paramètre géométrique $\PI$ n'a d'influence que sur
 la hauteur de la paroi spéculaire. Ainsi, la contribution radiative de la source
 émise par la paroi de droite, transportée directement vers le récepteur, reste
 identique pour toutes valeurs de $\PI$ (c'est à dire quelle que soit la hauteur
 de la paroi spéculaire). Du point de vue de la sensibilité
 cette contribution n'aura donc aucune influence.
 
 En pratique cela signifie que, en vue d'une validation via un calcul des
 différences finies, le calcul par Monte Carlo du flux au récepteur n'est pas
 entièrement nécessaire. Nous pouvons nous contenter d'estimer l'unique partie du
 flux qui sera perturbée par une variation de $\PI$, soit les contributions
 portées par les chemins réfléchis par la paroi spéculaire. En terme d'algorithme
 nous pouvons donc réutiliser les chemins déjà échantillonnés pour le problème
 couplé puisque leur statistique correspond exactement à celle d'un chemin émis
 par la paroi de droite et réfléchi par la paroi du haut.
 
 Le poids ([[weight_flux_part_spec]]) correspondant à la partie du flux venant
 de la réflection sur la paroi spéculaire est donc calculé au même moment que
 celui de la sensibilité.
 
 <<Données locales à la fonction de réalisation>>=
 double weight_flux_part_spec;
 @
 
 Le poids de cette contribution correspond à la source radiative $S_b$
 (equation~\ref{eq:cl_rad}) représentée par la variable [[Sb]] multipliées par
 l'angle $\pi$ ([[PI]]), la surface $A_h$ et le coefficient de réflection $\rho$
 ([[rho]]).
 
 <<Calcul du poids>>=
 weight_flux_part_spec = Sb * rho * PI * get_Sr(scene);
 @
 
 <<Renvoyer le poids>>=
 w[1] = weight_flux_part_spec;
 @
 
 Ce poids est également initialisé en même temps que celui de la sensibilité de
 sorte que les chemins réfléchis n'atteignant pas le récepteur aient une
 contribution nulle.
 
 <<Initialiser le poids>>=
 weight_flux_part_spec = 0;
 @
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Annexe fonction de calcul
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Fonction de calcul}
 \label{sec:fonction_de_calul}
 
 Le programme présenté jusqu'alors s'est concentré sur la mise en {\oe}uvre de
 la seule fonction de réalisation de notre algorithme de Monte Carlo; la boucle
 d'intégration, l'accumulation des poids ou encore l'affichage des résultats y
 sont absents. Dans cette section nous détaillons ces étapes manquantes afin de
 compléter le programme écrit jusqu'ici et ainsi proposer une mise en {\oe}uvre
 complète de l'algorithme de Monte Carlo objet du présent document. Ces étapes
 sont regroupées dans la fonction qui suit:
 
 <<Calculer la sensibilité à la translation>>=
 res_T
 compute_sensitivity_translation(struct sgs* sgs)
 {
   <<Variables locales au calcul de sensibilité>>
   res_T res = RES_OK;
 
   <<Exécuter l'intégration Monte Carlo>>
   <<Afficher les résultats de l'estimation>>
 
 exit:
   <<Libérer les variables locales au calcul de sensibilité>>
   return res;
 error:
   goto exit;
 }
 @
 
 Dans cette fonction on commence par exécuter l'intégration Monte Carlo. Cette
 étape consiste à invoquer notre [[<<Fonction de réalisation>>]] autant de fois
 que de nombre de réalisations demandées, et à accumuler les poids qu'elle
 retourne. D'apparence triviale, cette simple boucle s'avère plus compliquée
 pour qui souhaite paralléliser son calcul. Au delà des questions propres à une
 exécution parallèle, d'aucun doit s'assurer que chaque processus dispose d'une
 séquence de nombre aléatoires qui lui est propre. C'est pourquoi nous utilisons
 ici la bibliothèque \texttt{Star-MonteCarlo} en charge de cette intégration
 parallèle. Pour cela nous créons d'abord un système \texttt{Star-MonteCarlo},
 c'est à dire une variable qui matérialise la bibliothèque à l'échelle de notre
 programme ([[smc]]). Et nous le configurons pour qu'il utilise le journal
 d'évènement ([[logger]]), l'allocateur mémoire ([[allocator]]) et le nombre de
 processus légers ([[nthreads]]) soumis en entrée de la fonction via la variable
 [[sgs]].  Nous configurons enfin le type de générateur aléatoire à utiliser, en
 l'occurrence le générateur pseudo alétoire
 Mersenne-Twister~\cite{matsumoto1998mersenne} ([[SSP_RNG_MT19937_64]]). Nous
 pouvons alors lancer l'intégration Monte Carlo à proprement parler. Pour cela
 nous définissons un intégrateur qui détermine la fonction à appeler à chaque
 réalisation ([[run_realisation]]), le type de poids calculés ([[smc_doubleN]],
 {\ie} un vecteur de réels) et le nombre total de réalisations
 ([[nrealisations]]) défini comme paramètre d'entrée via la variable [[sgs]].
 
 <<Exécuter l'intégration Monte Carlo>>=
 /* Configurer et créer le système Star-MonteCarlo */
 smc_args.logger = &sgs->logger;
 smc_args.allocator = sgs->allocator;
 smc_args.nthreads_hint = sgs->nthreads;
 smc_args.rng_type = SSP_RNG_MT19937_64;
 res = smc_device_create(&smc_args, &smc);
 if(res != RES_OK) goto error;
 
 /* Configurer l'intégrateur */
 integrator.integrand = run_realisation;
 integrator.type = &smc_doubleN;
 integrator.max_realisations = sgs->nrealisations;
 ctx.count = 2; /* Nombre de poids calculés (Sensibilité & "luminance") */
 ctx.integrand_data = sgs; /* Données d'entrée de la fonction integrand */
 
 /* Intégration Monte Carlo */
 res = smc_solve(smc, &integrator, &ctx, &estimator);
 if(res != RES_OK) goto error;
 @
 
 En sortie de l'intégration Monte Carlo (fonction [[smc_solve]]) nous disposons
 d'un estimateur de nos variables aléatoires, en l'occurrence la sensibilité à
 la translation et la fraction de la luminance qui dépend du paramètre {$\PI$}.
 Nous pouvons dès lors récupérer l'état de notre estimateur pour afficher
 l'espérance et l'écart type de ces variables.
 
 <<Afficher les résultats de l'estimation>>=
 res = smc_estimator_get_status(estimator, &status);
 if(res != RES_OK) goto error;
 
 sgs_log(sgs, "Sensibilité ~ %g +/- %g\n",
   SMC_DOUBLEN(status.E)[0], /* Espérance */
   SMC_DOUBLEN(status.SE)[0]); /* Écart type */
 
 sgs_log(sgs, "Luminance ~ %g +/- %g\n",
   SMC_DOUBLEN(status.E)[1], /* Espérance */
   SMC_DOUBLEN(status.SE)[1]); /* Écart type */
 @
 
 Ne reste plus qu'à déclarer les variables locales utilisées par notre fonction
 de calcul et de libérer en sortie l'espace mémoire allouée dynamiquement pour
 ces variables.
 
 <<Variables locales au calcul de sensibilité>>=
 /* Système */
 struct smc_device_create_args smc_args = SMC_DEVICE_CREATE_ARGS_DEFAULT;
 struct smc_device* smc = NULL;
 
 /* Intégrateur */
 struct smc_integrator integrator = SMC_INTEGRATOR_NULL;
 struct smc_doubleN_context ctx = SMC_DOUBLEN_CONTEXT_NULL;
 
 /* Résultat de l'estimatation */
 struct smc_estimator* estimator = NULL;
 struct smc_estimator_status status = SMC_ESTIMATOR_STATUS_NULL;
 @
 
 <<Libérer les variables locales au calcul de sensibilité>>=
 if(estimator) smc_estimator_ref_put(estimator);
 if(smc) smc_device_ref_put(smc);
 @
 
 Le lecteur attentif aura remarqué que l'intégrateur utilise la fonction
 [[run_realisation]] et non directement la [[<<Fonction de réalisation>>]]
 développée dans ce document (voir [[<<Exécuter l'intégration Monte Carlo>>]]).
 [[run_realisation]] est une fonction intermédiaire qui ne fait qu'appeler la
 fonction de réalisation. [[run_realisation]] est donc parfaitement dispensable
 sauf à la bibliothèque \texttt{Star-MonteCarlo} qui nous impose la signature de
 la fonction à utiliser, c'est à dire le type des paramètres d'entrées et de
 sorties. En d'autres termes, l'utilisation de cette fonction intermédiaire nous
 permet de faciliter l'écriture et la lecture de la fonction de réalisation en
 la libérant de contraintes fonctionnelles imposées par
 \texttt{Star-MonteCarlo}.
 
 <<Fonctions utilitaires>>=
 static res_T
 realisation
   (struct ssp_rng* rng,
    const struct sgs_scene* scene,
    double* weight);
 
 static res_T
 run_realisation
   (void* output,
    struct ssp_rng* rng,
    const unsigned ithread,
    void* ctx)
 {
   struct smc_doubleN_context* context = NULL;
   struct sgs* sgs = NULL;
   ASSERT(ctx && output);
   (void)ithread; /* Éviter l'avertissement "variable inutilisée" */
   context = ctx;
   sgs = context->integrand_data;
   return realisation(rng, &sgs->scene, SMC_DOUBLEN(output));
 }
 @
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Annexe structure du C
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Structure de mise {\oe}uvre}
 
 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>>
 
 <<Fonctions utilitaires>>
 <<Fonction de réalisation>>
 <<Calculer la sensibilité à la translation>>
 @
 
 En plus des fonctions de calcul écrites dans les sections précédents, notre
 fichier contient des fonctions utilitaires notamment 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
 get_rho
   (const struct sgs_scene* scene,
    const double pos[2])
 {
   ASSERT(scene && pos);
 
   return 0.25
     * (1 - cos(2*PI*pos[X]/scene->D[X]))
     * (1 - cos(2*PI*pos[Y]/scene->D[Y]));
 }
 
 static void
 get_grad_rho
   (const struct sgs_scene* scene,
    const double pos[2],
    double grad[2])
 {
   ASSERT(scene && pos && grad);
 
   grad[X] = 0.25
     * (((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->D[Y])*sin(2*PI*pos[Y]/scene->D[Y]))
     * (1 - cos(2*PI*pos[X]/scene->D[X]));
 }
 
 static double
 get_Sb
   (const struct sgs_scene* scene,
    const double pos[2])
 {
   ASSERT(scene && pos);
   return
     (1 - cos(2*PI*pos[X]/scene->D[Y]))
   * (1 - cos(2*PI*pos[Y]/scene->D[Z]));
 }
 
 static void
 get_grad_Sb
   (const struct sgs_scene* scene,
    const double pos[2],
    double grad[2])
 {
   ASSERT(scene && pos && grad);
 
   grad[X] =
       (((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->D[Z])*sin(2*PI*pos[Y]/scene->D[Z]))
     * (1 - cos(2*PI*pos[X]/scene->D[Y]));
 }
 @
 
 Autres fonctions utilitaires, les fonctions utilisées lors du suivi des 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 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 && 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 */
   }
 }
 @
 
 Pour le récepteur on teste d'abord si l'intersection a lieu sur la paroi du
 parallélépipède 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
   (const struct sgs_scene* scene,
    const double ray_org[3],
    const double ray_dir[3],
    const struct sgs_hit* hit)
 {
   double hit_pos[3];
   ASSERT(scene && ray_org && ray_dir && hit);
 
   /* Le rayon n'intersecte pas la surface où le récepteur se trouve */
   if(SGS_HIT_NONE(hit) || hit->surface != SGS_SURFACE_Z_MIN) {
     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];
 
   /* Le rayon n'intersecte pas le récepteur*/
   if(hit_pos[X] < scene->recv_min[X] || hit_pos[X] > scene->recv_max[X]
   || hit_pos[Y] < scene->recv_min[Y] || hit_pos[Y] > scene->recv_max[Y]) {
     return 0;
 
   /* Le rayon intersecte le récepteur*/
   } else {
     return 1;
   }
 }
 @
 
 Dernière fonction utilitaire à écrire, la fonction qui calcule la surface $A_h$
 de la paroi spéculaire:
 
 <<Fonctions utilitaires>>=
 static double
 get_Sr(const struct sgs_scene* scene)
 {
   ASSERT(scene);
   return scene->D[X] * scene->D[Y];
 }
 @
 
 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};
 @
 
 Enfin, on liste en début des sources les fichiers d'en-tête utilisés par notre
 code.
 
 <<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>
 @
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Annexe script de résultat
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{Script d'exécution du calcul et résultats}
 \label{annex:scripts}
 
 Nous écrivons ici les scripts \textit{shell} qui calculent la sensibilité du
 flux à $\PI$ d'abord par Monte Carlo puis par différences finies en vue d'une
 validation croisée des résulats.  Le premier script [[<<mc.sh>>]] lance
 plusieurs fois le programme dont le présent document décrit la fonction de
 réalisation; l'objet étant de calculer à différentes valeurs de $\PI$ la
 sensibilité du flux $\varphi$ reçu par le capteur (équation~\ref{eq:flux})
 ainsi que la fraction de ce même flux qui dépend de $\PI$
 (annexe~\ref{annex:flux}).
 
 <<mc.sh>>=
 #!/bin/sh -e
 
 # Estimation par Monte Carlo de la sensibilité et de la fraction du
 # flux qui dépend de pi
 
 <<Configuration géométrique>>
 <<Paramètres des calculs Monte Carlo>>
 
 <<Pour différentes valeurs de $\PI$>>
 do
   <<Lancer un calcul Monte Carlo>>
   <<Post traiter le résultat>>
   <<Passer à la valeur de $\PI$ suivante>>
 done
 @
 
 Chaque calcul Monte Carlo se résume à exécuter ledit programme nommé [[sgs]]
 (acronyme de \textit{Star Geometric Sensitivity}) pour une valeur de $\PI$
 considérée. On notera que la configuration géométrique décrite en
 figure~\ref{fig:configuration} est indépendante de ses dimensions réelles. Dans
 notre script ces dimensions sont fixées via deux variables qui définissent le
 coin inférieur ([[lower]]) et le coin supérieur ([[upper]]) d'un parallélépipède
 aligné aux axes, deux variables passées en arguments du programme via l'option
 [[-b]].
 
 <<Lancer un calcul Monte Carlo>>=
 out=$(./sgs \
   -n "${nrealisations}" \
   -b low="${lower}":upp="${upper}":pi="${pi}")
 @
 
 avec [[nrealisations]] le nombre de réalisations utilisées par le calcul Monte
 Carlo. Si ce document décrit en détail les sources C de sa
 [[<<Fonction de réalisation>>]], nous renvoyons le lecteur vers les autres
 fichiers C et l'aide du programme [[sgs]] (affichée via l'option [[-h]]) pour
 plus d'informations quant aux fonctionnement et options du programme.
 
 Une fois le calcul terminé, en post-traiter le résultat se résume à afficher
 sur la sortie standard la valeur de $\PI$ courante suivie de l'estimation et de
 l'écart type de la sensibilité et du flux que nous venons d'estimer. Ci-après
 nous utilisons la commande [[sed]] pour extraire et afficher ces valeurs
 stockées dans la variable [[out]] à l'issu de notre calcul Monte Carlo.
 
 <<Post traiter le résultat>>=
 prompt="[^~]\{1,\}~ "
 estimation="[^[:blank:]]\{1,\}"
 error="[^\n$]\{1,\}"
 line0="${prompt}\(${estimation}\) +/- \(${error}\)\n" # Sensibilité
 line1="${prompt}\(${estimation}\) +/- \(${error}\)$" # Flux
 printf "%s " "${pi}"
 echo "${out}" | sed -n "1{N;s#^${line0}${line1}#\1 \2 \3 \4#p}"
 @
 
 Ces deux étapes, à savoir le calcul Monte Carlo et son post-traitement, sont
 lancées [[nsteps]] fois pour différentes valeurs de $\PI$, la valeur de $\PI$ à
 chaque itération étant simplement sa valeur précédente incrémentée d'un pas
 constant égal à [[pi_step]].
 
 <<Pour différentes valeurs de $\PI$>>=
 i=0
 pi=0
 while [ "${i}" -lt "${nsteps}" ]
 @
 <<Passer à la valeur de $\PI$ suivante>>=
 i=$((i + 1))
 pi=$(printf "%s + %s\n" "${pi}" "${pi_step}" | bc)
 @
 
 Pour compléter le script, ne reste plus qu'à définir les variables qui
 caractérisent notre configuration géométrique ainsi que les paramètres qui
 pilotent nos différents calculs, tels que le nombre de calculs à lancer ou
 encore le nombre de réalisations par calcul.
 
 <<Configuration géométrique>>=
 h=1 # Hauteur du parallélépipède
 lower="0,0,0"
 upper="1,1,${h}"
 @
 <<Paramètres des calculs Monte Carlo>>=
 nrealisations=100000000
 nsteps=14
 pi_step=0.1
 @
 
 Dans un nouveau script [[fd.sh]] nous pouvons alors calculer par différences
 finies la sensibilité du flux à $\PI$ à partir des estimations Monte Carlo en
 sortie de [[<<mc.sh>>]], des résultats lus via l'entrée standard. Pour pouvoir
 comparer les résultats, ce second script recopie en sortie les sensibilités
 estimées par Monte Carlo ([[sen_mc]]) en plus d'écrire ces mêmes sensibilités
 calculées cette fois par différences finies ([[sen_fd]]).  Leur écarts types
 respectif ([[err_mc]] et [[err_fd]]) sont également des données de sortie. On
 notera enfin que l'ensemble des résultats sont adimentionnalisés et sont par
 conséquent donnés pour une valeur de $\PI$ indépendante de la hauteur du
 parallélépipède ([[pi_over_h]]).
 
 <<fd.sh>>=
 #!/bin/sh -e
 
 # Calcule la sensibilité par différences finies à partir des
 # estimations par Monte Carlo de la fraction du flux qui dépend de pi
 
 # Liste des données pour chaque ligne écrite en sortie
 printf "pi_over_h sen_mc err_mc sen_fd err_fd\n"
 
 <<Fonctions utilitaires à [[fd.sh]]>>
 
 <<Lire les paramètres d'entrée>>
 
 <<Pour chaque valeur de $\PI$ considérée>>
 do
   <<Lire la valeur du flux autour de $\PI$>>
   <<Calculer par différences finies la sensibilité du flux à $\PI$>>
   <<Lire l'estimation Monte Carlo de la sensibilité du flux à $\PI$>>
   <<Écrire les résultats adimensionnaliser>>
   <<Passer au calcul suivant>>
 done
 @
 
 Pour les calculs en différences finies nous avons besoin de connaître les
 paramètres utilisés par les estimations Monte Carlo, comme le $\delta$ entre
 chaque valeur de $\PI$ ([[pi_step]]) ou encore la hauteur du parallélépipède
 nécessaire pour adimensionnaliser les résultats ([[h]]). Pour passer ces
 données d'un script à l'autre, on ajoute aux sorties de [[mc.sh]] un en-tête
 contenant ces paramètres, en-tête qui peut alors être lu par le script
 [[fd.sh]].
 
 <<Paramètres des calculs Monte Carlo>>=
 printf "%s %s\n" "${h}" "${pi_step}"
 @
 <<Lire les paramètres d'entrée>>=
 read -r header
 h=$(echo "${header}" | cut -d' ' -f1)
 pi_step=$(echo "${header}" | cut -d' ' -f2)
 @
 
 Autre donnée nécessaire par l'adimensionnement, la valeur maximale du flux.
 Celle-ci correspond à la valeur de $\varphi$ lorsque la boîte est fermée, c'est
 à dire lorsque $\PI$ vaut $0$, soit la valeur de $\PI$ du premier calcul Monte
 Carlo lancé par [[mc.sh]]. Après avoir lu l'en-tête de ses sorties nous lisons
 donc ce premier résultat que l'on stocke dans la variable [[p]] avant d'en
 extraire la valeur $\varphi$ ([[phi_max]]).
 
 <<Lire les paramètres d'entrée>>=
 read -r p
 phi_max=$(echo "${p}" | cut -d' ' -f4 | float_to_bc)
 @
 
 La boucle principale du script consiste à lire l'entrée standard jusqu'à ce
 qu'il n'y ait plus de résultat Monte Carlo à traiter. À chaque itération, le
 résultat Monte Carlo pour la valeur de $\PI$ courante est stocké dans la
 variable [[c]] tandis que le résultat précédent et suivant sont respectivement
 stockés dans les variables [[p]] et [[n]].
 
 <<Pour chaque valeur de $\PI$ considérée>>=
 read -r c
 while read -r n
 @
 <<Passer au calcul suivant>>=
 p="${c}"
 c="${n}"
 @
 
 Il suffit alors d'extraire le flux ([[phi]]) et son erreur ([[err]]) aux
 valeurs de $\PI$ précédente ([[p]]) et suivante ([[n]]) pour calculer par
 différences finies la sensibilité et l'erreur associée pour la valeur de $\PI$
 courante:
 
 <<Lire la valeur du flux autour de $\PI$>>=
 phi_p=$(echo "${p}" | cut -d' ' -f4 | float_to_bc)
 err_p=$(echo "${p}" | cut -d' ' -f5 | float_to_bc)
 phi_n=$(echo "${n}" | cut -d' ' -f4 | float_to_bc)
 err_n=$(echo "${n}" | cut -d' ' -f5 | float_to_bc)
 @
 <<Calculer par différences finies la sensibilité du flux à $\PI$>>=
 sen_fd=$(echo "(${phi_n}-${phi_p})/(2*${pi_step})" | bc_cmd)
 err_fd=$(echo "(${err_n}+${err_p})/(2*${pi_step})" | bc_cmd)
 @
 
 On utilise alors le résultat Monte Carlo au $\PI$ considéré pour retrouver non
 seulement la valeur de $\PI$ mais aussi pour extraire l'estimation Monte Carlo
 de la sensibilité de $\varphi$ à $\PI$:
 
 <<Lire l'estimation Monte Carlo de la sensibilité du flux à $\PI$>>=
 pi=$(echo "${c}" | cut -d' ' -f1 | float_to_bc)
 sen_mc=$(echo "${c}" | cut -d' ' -f2 | float_to_bc)
 err_mc=$(echo "${c}" | cut -d' ' -f3 | float_to_bc)
 @
 
 Ne reste plus qu'à écrire l'ensemble des résultats attendus:
 
 <<Écrire les résultats adimensionnaliser>>=
 pi_over_h=$(echo "${pi}/${h}" | bc_cmd)
 sen_mc=$(echo "${sen_mc}/${phi_max}*${h}" | bc_cmd)
 err_mc=$(echo "${err_mc}/${phi_max}*${h}" | bc_cmd)
 sen_fd=$(echo "${sen_fd}/${phi_max}*${h}" | bc_cmd)
 err_fd=$(echo "${err_fd}/${phi_max}*${h}" | bc_cmd)
 printf "%s %s %s %s %s\n" \
   "${pi_over_h}" "${sen_mc}" "${err_mc}" "${sen_fd}" "${err_fd}"
 @
 
 On complète finalement notre script par les fonctions utilitaires utilisées
 tout du long, à savoir la fonction [[float_to_bc]] qui convertie un nombre à
 virgule flottante dans le format attendu par le programme \texttt{bc}, et la
 fonction [[bc_cmd]] qui exécute le programme \texttt{bc} et en post-traite le
 résultat pour supprimer les zéros qui suivent le dernier chiffre significatif.
 
 <<Fonctions utilitaires à [[fd.sh]]>>=
 float_to_bc()
 {
   sed 's/\([+-]\{0,1\}[0-9]\{0,\}\.\{0,1\}[0-9]\{1,\}\)'\
 '[eE]+\{0,1\}\(-\{0,1\}\)\([0-9]\{1,\}\)/(\1*10^\2\3)/g'
 }
 
 bc_cmd()
 {
   bc -l | sed '/\./s/\.\{0,\}0\{1,\}$//'
 }
 @
 
 \bibliographystyle{apalike}
 \bibliography{biblio}
 
 \end{document}
 \nwfilename{test}