TP1 : transmissivité d'une couche : fit du modèle de Malkmus

L'objectif du TP est de générer une courbe de la transmissivité d'une couche en fonction de son épaisseur, puis de l'analyser pour en déduire les paramètres du modèle de Malkmus fitté à cette courbe.

On suppose que l'on travaille dans un répertoire ${HOME}/mon_nom/tp_star-line/TP1 conformément aux instructions données dans la présentation générale du TP.

export LANG=en_US.UTF-8

Partie 1 : Générer la donnée du slab

Dans le cadre de ce TP, la donnée a été générée au préalable et est déjà disponible dans le répertoire data mais on redonne ici la commande qui a permis de générer le fichier mixture.txt qui décrit le mélange de gaz (attention si vous copiez collez la commande, il ne faut pas d'espace en début de ligne) :

cat > data/mixture.txt <<EOD
# molecule    concentration [mol/mol]     cutoff [cm^-1]
H2O           1.9E-2                      25
    # isotope-id    abundance
    161             9.97E-01
    181             2.00E-03
    171             1.00E-03

# molecule    concentration [mol/mol]     cutoff [cm^-1]
CO2           2.87E-4                     50
O3            6.15E-6                     25
EOD

Pour obtenir les donnés associées aux transitions moléculaires pour ce mélange, il faut aller sur le site de Hitran, créer un compte, commander la donnée à partir du formulaire et télécharger les fichiers résultant. Pour ce TP, cette étape a été faite à l'avance et les fichiers sont dans le répertoire data/. Il y a un fichier readme-1.txt qui documente la donnée, des fichiers bib contenant les références scientifiques associées à la base de donnée, et enfin les données commandées dans un fichier .par. A l'origine, le nom du fichier de donnée est le numéro de commande. Pour ce TP, on a renommé le fichier en mixture_lines.par.

Remarque :

Si vous changez la description du mélange de gaz ci dessus, par exemple en rajoutant des isotopes ou des molécules, il faudra regénérer les fichiers
.par fournis par Hitran. Ou par HITEMP si vous êtes en haute température.

https://hitran.org/lbl/ # to get the .par file
https://hitran.org/media/molparams.txt

Remarque :

Si on a besoin de compresser les données de transitions (fichier .par livré
par hitran), on peut utiliser l'utilitaire shtr pour générer un fichier
plus léger, dans un format compressé et qui sera lisible par les autres
utilitaires de star-line (le format shtr)

shtr -l data/mixture_lines.par -o data/mixture_lines.shtr -vvv -iH

Partie 2 : Créer l'arbre de transitions

Une fois la donnée brute générée, on va construire l'arbre structurant cette donnée. La donnée brute correspond à un grand nombre de transitions. L'algorithme de transfert radiatif par Monte Carlo, formulé en collisions nulles, a besoin d'un majorant des valeurs du coefficient d'absorption. Le coefficient d'absorption est lui-même une somme des contributions de toutes les transitions. On rappelle que dans le cadre théorique développé dans les thèses de M. Galtier et Y. Nyffenegger-Pere et qu'implémente star-line, on ne calcule jamais les valeurs du coefficient d'absorption, on échantillonne directement les transitions. Ainsi, l'arbre que l'on va construire a deux objectifs : (1) produire un spectre majorant, (2) structurer les données décrivant les transitions de façon à accélérer l'échantillonnage d'une transition parmis toutes celles de la donnée brute.

Pour la construction de l'arbre, on passera par l'utilitaire sln-build qui est distribué avec star-line, dont la syntaxe est donnée dans la page de manuel :

Syntaxe :

sln-build [-chsv] [-a arity] [-e polyline_opt[:polyline_opt ...]] \
[-L leaf_nlines] [-l line_profile] [-o accel_struct] \
[-t thread_count] -m molparams -P pressure -T temperature \
-x mixture [lines] \

Faire man sln-build et étudier la commande et ses différentes options.

On va construire ici un arbre pour le slab dont les concentrations gazeuses ont été décrites précédemment, et on y associe les conditions de pression et de température à l'aide des options -P et -T :

sln-build -P 1 -T 296 -o data/mixture_tree.bin \
           -m data/molparams.txt -x data/mixture.txt \
           data/mixture_lines.par

L'arbre ainsi construit a été stocké dans le fichier data/mixture_tree.bin

Partie 3 : Calcul de la transmissivité d'un slab

rm -f thickness.dat
rm -f transmissivities.dat

Pour calculer la transmissivité d'une couche, il faut commencer par définir une plage de nombre d'onde. Le programme sln-slab fait l'hypothèse d'une illumination uniforme sur la plage spectrale, ce qui est réaliste si celle-ci est mince.

nu_range=100,125 # [cm⁻¹]

Ensuite, on génère une liste d'épaisseurs de couche pour voir comment évolue la transmissivité en fonction de l'épaisseur. Si l'on modifie la plage spectrale, il faudra sûrement modifier également les épaisseurs car le coefficient d'atténuation a une forte dépendance spectrale.

thin_limit="$(seq -w 0.0005 0.0005 0.01)"
thickness="$(seq -w 0.05 0.05 1)"
thick_limit="$(seq -w 1 0.1 2)"
list_thickness="$thin_limit $thickness $thick_limit"

Enfin, on définit le nombre de réalisations pour sln-slab. Là encore, cette valeur sera à adapter à la plage spectrale choisie. En particulier, lorsque la couche est très optiquement mince, il va falloir beaucoup de réalisations pour faire converger la transmissivité.

samples_count=100000

On peut maintenant lancer sln-slab pour les paramètres choisis. On enregistre l'épaisseur de la couche, sa transmissivité et l'écart type de la transmissivité dans le fichier results.dat.

for thickness in $list_thickness ; do
  echo $thickness >> thickness.dat
  sln-slab -S $nu_range -T $thickness \
    -n $samples_count \
    -a data/mixture_tree.bin \
    -m data/molparams.txt \
    -sl data/mixture_lines.shtr \
    |grep transmiss|cut -d":" -f2|sed 's,+/-,,'|sed 's/;.*$//' \
    >> transmissivities.dat
done
paste thickness.dat transmissivities.dat > results.dat

Partie 4 : Fit du modèle de Malkmus sur les données générés par sln-slab

Le document ${TP_STAR_LINE}/Doc_spectral/fiche_malkmus.pdf (également disponible en ligne) décrit le modèle de Malkmus, qui donne l'expression de la transmissivité moyenne d'une couche d'épaisseur l sous hypothèses formulées sur la distribution des raies (cf. equation (22)).

Le modèle de Malkmus a deux parmètres, k_barre, le coefficient d'absorption moyen sur la bande, et phi, un paramètre de forme du spectre (phi est proche de zéro quand les raies sont très piquées et se recouvrent peu, et grand devant 1 lorsque les raies se recouvrent beaucoup).

Pour obtenir ces deux paramètres à partir de la courbe de transmissivité calculée par sln-slab, on va s'aider des limites faibles (lorsque l'épaisseur optique est petite, les centres de raies ne sont pas saturés) et fortes (lorsque l'épaisseur optique est grande, les centres de raies sont saturés et ce sont les ailes de raies qui vont déterminer l'évolution de la transmissivité).

On va utiliser les limites de l'équation (22) pour déterminer les coefficients k_barre et phi.

En limite faible (k_barre*l est petit devant 1), la transmissivité s'exprime comme :

T = 1-k_barre*l

On peut donc obtenir k_barre par régression linéaire de la transmissivité pour de petites épaisseurs de couche. On propose d'utiliser les 5 premières valeurs de transmissivité pour le calculer :

n=5
k=$(awk -v n=$n \
  '{ if(NR<n){sum += (1-$2)/$1; i++}} \
  END {print sum / i}' results.dat)

De même, lorsque k_barre*l est grand devant 1, l'équation (22) peut s'approximer par :

T = exp(-sqrt(2k_barrel*phi))

On propose donc d'utiliser les 10 dernières valeurs d'épaisseur optique données par sln-slab pour estimer phi :

n=40
phi=$(awk -v n=$n -v k=$k \
  '{ if(NR>n){sum += log($2)^2/$1; i++}} \
  END {print sum / (2*k*i)}' results.dat)

Enfin, on peut comparer la transmissivité donnée par sln-slab au modèle de Malkmus avec les paramètres k_barre et phi ainsi déterminés :

gnuplot -p -e "k=$k" -e "phi=$phi" \
  gnuplot_scripts/slab_transmissivity.gp