FR 3033063 A1 20160826 FR 1551589 A 20150224 1 PROCEDE DE CALCUL NUMERIQUE DE LA DIFFRACTION D'UNE STRUCTURE L'invention concerne la diffusion/diffraction d'une onde électromagnétique par les structures matérielles complexes, et en particulier la modélisation de propriétés de diffusion/diffraction et le calcul numérique de diffusion/diffraction pour des structures hétérogènes et d'épaisseur non négligeable par rapport à la longueur d'onde. Dans de nombreux domaines techniques, tels que la photolithographie ou les éléments optiques diffractifs, il s'avère primordial de pouvoir modéliser et calculer précisément la réponse d'un composant à un faisceau lumineux incident. La modélisation vise à calculer la diffraction d'une structure caractérisée par la répartition spatiale de la permittivité diélectrique définie suivant différents axes en chaque point de la structure. La modélisation électromagnétique d'une structure est effectuée dans l'état de la technique pour permettre soit un procédé de résolution approchée, soit un procédé de résolution statistique, soit un procédé de résolution exact. Les procédés de résolution approchés permettent d'obtenir un résultat de calcul de diffraction relativement rapide. Cependant, de tels procédés présentent 20 une précision insuffisante pour obtenir un résultat de calcul satisfaisant. Les procédés de résolution exacts de l'état de la technique sont très lents et requièrent une très grande capacité de mémoire numérique, ce qui limite leur application à des éléments optiques de faible complexité. Les procédés de résolution exacts sont donc le plus souvent utilisés uniquement comme moyen de 25 vérification ponctuelle d'une modélisation effectuée par un procédé de résolution approchée. Parmi les procédés de résolution exacts connus de l'état de la technique, on peut notamment citer le procédé de calcul au moyen des matrices de diffusion S. Pour cela, la structure à modéliser est décomposé en un ensemble de N 30 couches jointives et M ordres de diffraction. La matrice de diffusion Si de chaque couche d'indice i est calculée et permet d'exprimer les amplitudes de champ sortantes en fonction des amplitudes incidentes sur les deux faces de cette couche. On peut ainsi noter qu'une amplitude sortante d'une couche constitue une amplitude incidente pour une couche adjacente. 35 'il =5' 'il' ) Une matrice S pour le composant est formée en combinant l'ensemble des matrices Si des différentes couches. 3033063 2 S So Si 0 ... 0 SN_I La modélisation permet alors de déterminer les amplitudes de champ sortantes de la structure formée de l'ensemble des couches en fonction des 5 amplitudes de champ incidentes sur cette structure : (JV S bo bN Cependant, cette combinaison des matrices Si n'est pas un simple produit 10 mais un calcul complexe. Le temps de calcul est proportionnel à M3 et en général linéaire en N. Une méthode alternative à la méthode de la matrice S est connue par la publication originale de Bremmer et par le développement suivant de Sluijter pour le calcul de la transmission et la réflexion d'un système de couches uniformes 15 superposées par le calcul de la transmission et de la réflexion de couche à couche d'une onde plane d'un bord à l'autre du multicouche suivi du même calcul en sens inverse avec mémorisation des amplitudes intermédiaires, cet aller et retour étant répété de façon itérative et les résultats des itérations successives étant sommés jusqu'à ce que la somme converge. Cette méthode de propagation aller et retour à travers la structure a été étendue à des structures diffractantes composées de couches microstructurées latéralement, chaque pas dans un aller et retour calculant ici l'amplitude des modes de chaque couche obtenus par la méthode décrite dans le document de D.M. Pai and K.A. Awada, "Analysis of dielectric gratings of arbitrary profiles and thicknesses", J. Opt. Soc. Am. A 8, 755-762. La validité de cette méthode a été évaluée dans la publication M. Nevière and F. Montiel, "Deep gratings: a combination of the differential theory and the multiple reflection series", Opt. Commun. 108, 1-7 qui conclut que cette sommation explicite des résultats itératifs ne converge que pour des structures où la modulation d'indice de réfraction ou des interfaces est faible et représente une perturbation d'une structure de base. Le document EP2302360 décrit un procédé de modélisation des propriétés de diffraction de structures microscopiques périodiques et un procédé de calcul de diffraction associé. Le procédé exprime le problème par une forme intégrale volumique d'un champ vectoriel en remplacement du champ électrique. Le champ vectoriel est obtenu à partir du champ électrique par un changement de base, de façon à présenter une continuité aux limites du matériau. Des convolutions sont réalisées sur le champ vectoriel en utilisant des opérateurs de convolution, selon les séries finies de Laurent. On peut ainsi réaliser des produits de matrices au moyen de transformations de Fourier rapides. Un opérateur de convolution et de 3033063 3 changement de base est configuré pour transformer le champ vectoriel vers le champ électrique recherché, par l'intermédiaire d'un changement de base fonction des propriétés de matériau et de géométrie de la structure périodique. Le document 'New fast and memory-sparing method for rigorous 5 electromagnetic analysis of 2D periodic dielectric structures', par Ms Shcherbakov et Tishchenko, publié dans Journal of Quantitative Spectroscopy & Radiative Transfer aux pages 158-171, décrit un autre procédé de résolution exacte d'une modélisation de propriétés électromagnétiques d'une structure. Ce procédé est en particulier adapté aux structures diélectriques présentant une périodicité dans 10 le plan, tels que des réseaux de diffraction. Le procédé de modélisation découpe la structure en différentes couches planes parallèles à un plan XY dans un système de coordonnées cartésien XYZ. Chaque couche a une épaisseur respective selon la direction Z. Un réseau de diffraction bidimensionnel est modélisé par une variation périodique de la 15 permittivité diélectrique d'une couche selon deux directions différentes définies par des vecteurs du réseau dans le plan XY. Le procédé de résolution transforme l'équation d'onde sous forme d'une équation intégrale implicite. Le procédé de résolution convertit ensuite l'équation d'ondes dans l'espace réciproque selon les axes X et Y. Du fait du découpage de 20 la structure sous forme de couches de mêmes épaisseurs, la formulation intégrale peut être exprimée sous forme d'équations matricielles, correspondant à des sommes d'harmoniques et exprimant l'empilement des couches sous forme de sommes. Les composantes du champ électrique normale et tangente sont alors séparées. Une forme bloc-Tôplitz des matrices peut être établie sans nécessiter 25 d'inversions matricielles. La résolution consiste alors à réaliser une inversion matricielle d'une matrice qui peut être exprimée comme des produits de matrices bloc-diagonales et de matrices bloc-Tôplitz. L'inversion matricielle est en pratique calculée par des multiplications matricielles, avec une méthode de résolution d'équations linéaires itérative de type GMRES, plutôt que par un calcul direct 30 d'inversion matricielle. La forme block-Tôplitz permet de réaliser des calculs par transformation rapide de Fourier, avec un temps de calcul alors sensiblement proportionnel à M, et une utilisation de mémoire numérique sensiblement proportionnelle à M également. Le temps de calcul obtenu avec ce procédé de résolution est proportionnel 35 à N pour des structures simples et peu épaisses. Cependant, le temps de calcul et l'occupation mémoire augmentent rapidement avec l'épaisseur des couches de la structure à modéliser. Un grand nombre d'itérations est alors en effet nécessaire pour obtenir une convergence de la méthode itérative, ce qui se traduit par une augmentation sensible du temps de calcul et de la mémoire numérique 3033063 4 nécessaire. En outre, la méthode de résolution s'avère mal adaptée pour des structures comportant des hétérogénéités de couches très différentes dans la structure, ce qui remet en cause certaines hypothèses de la méthode de résolution. Par ailleurs, le calcul reste exigeant en quantité de mémoire numérique 5 utilisée, les données de résolution de l'ensemble de la structure devant rester mémorisées durant tout le calcul. Cette quantité de mémoire de calcul nécessaire limite fortement l'épaisseur des structures pouvant être modélisées. L'invention vise à résoudre un ou plusieurs de ces inconvénients. L'invention porte ainsi sur un procédé de calcul numérique de la diffraction d'une 10 structure dont un modèle numérique définit la permittivité diélectrique en chaque point, comprenant les étapes de : -définition d'une modélisation numérique d'un modèle de diffraction d'une structure, incluant : -la scission du modèle numérique de permittivité diélectrique de la structure 15 en un nombre N de modèles numériques de permittivité diélectrique de couches planes respectives superposées selon une direction et ordonnées selon un indice i compris entre 1 et N selon ladite direction ; -la détermination du modèle de diffraction de chacune desdites couches à partir du modèle numérique de permittivité diélectrique de cette couche, le 20 nombre N de couches étant déterminé pour que l'occupation de mémoire numérique pour le calcul de diffraction par le modèle de diffraction de chaque couche soit inférieur à K* M, avec M un nombre d'ordres de diffraction de cette couche et K un facteur indépendant de M; -calcul numérique de la diffraction de la structure, incluant : 25 -une itération initiale durant laquelle on réalise : -une propagation d'onde dans le sens des couches 1 vers N, incluant : -le calcul d'une onde réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice 1 ; -pour une couche d'indice i compris entre 2 et N, calcul d'une onde 30 réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice i, cette onde incidente étant l'onde transmise calculée pour la couche d'indice i-1 ; -mémorisation des ondes réfléchies calculées pour chacune des couches d'indice i et de l'onde transmise calculée pour la couche d'indice N ; 35 -et une propagation d'onde dans le sens des couches N vers 1, incluant : -une onde réfléchie et une onde transmise sont calculées par application d'une onde incidente au modèle de diffraction de la couche d'indice N ; -pour une couche d'indice i compris entre N-1 et 1, une onde réfléchie et une onde transmise sont calculées par application d'une onde incidente au 3033063 5 modèle de diffraction de la couche d'indice i, cette onde incidente incluant l'onde transmise calculée pour la couche d'indice i+1 ; -mémorisation des ondes réfléchies calculées pour chacune des couches d'indice i et de l'onde transmise calculée pour la couche d'indice N ; 5 -des itérations ultérieures, chaque itération ultérieure d'indice k incluant : -une propagation d'onde dans le sens des couches 2 vers N, incluant : -pour une couche d'indice i compris entre 2 et N, calcul d'une onde réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice i, cette onde incidente incluant l'onde 10 transmise calculée pour la couche d'indice i-1 durant cette itération et incluant l'onde réfléchie calculée pour la couche d'indice i-1 durant la dernière propagation dans le sens opposé; -mémorisation des ondes réfléchies calculées pour chacune des couches d'indice i et de l'onde transmise calculée pour la couche d'indice N ; 15 -une propagation d'onde dans le sens des couches N-1 vers 1, incluant : -pour une couche d'indice i compris entre N-1 et 1, calcul d'une onde réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice i, cette onde incidente incluant l'onde transmise calculée pour la couche d'indice i+1 durant cette itération et incluant 20 l'onde réfléchie calculée pour la couche d'indice i+1 durant la dernière propagation dans le sens opposé; -mémorisation des ondes réfléchies calculées pour chacune des couches d'indice i et de l'onde transmise calculée pour la couche d'indice 1 ; -un test de convergence de la solution Vk formée des ondes réfléchies et 25 transmises calculées lors de l'itération d'indice k, par application d'une méthode itérative de résolution de systèmes d'équations linéaires à la résolution de l'équation implicite de convergence Vk=P.Vk, avec P un opérateur correspondant à la transformation de la solution de toute solution Vk-1 en une solution Vk à l'itération k. 30 Selon une variante, ladite méthode itérative de résolution de systèmes d'équations linéaires est une méthode de type GMRES. Selon une autre variante, ladite méthode itérative de résolution de systèmes d'équations linéaires est une méthode du gradient biconjugué stabilisée. 35 Selon encore une variante, ladite détermination du modèle de diffraction de chacune desdites couche est mise en oeuvre par un procédé GSM. Selon encore une autre variante, N est au moins égal à 5. 3033063 6 Selon une variante, le calcul de propagation dans le sens des couches 2 à N et le calcul de propagation dans le sens des couches N-1 à 1 sont réalisés en parallèle. Selon encore une variante, le calcul de propagation dans le sens des 5 couches 2 à N et le calcul de propagation dans le sens des couches N-1 à 1 sont réalisés séquentiellement. Selon une autre variante, le calcul de propagation comprend des calculs en parallèle pour des couches différentes. 10 D'autres caractéristiques et avantages de l'invention ressortiront clairement de la description qui en est faite ci-après, à titre indicatif et nullement limitatif, en référence aux dessins annexés, dans lesquels : -la figure 1 est une vue en coupe d'un exemple de structure à modéliser scindé en différentes couches ; 15 -la figure 2 est une vue en coupe schématique d'un exemple de structure simplifiée de type OLED, à modéliser, scindée en différentes couches ; -la figure 3 représente schématiquement un système configuré pour modéliser et calculer numériquement la diffraction d'une structure ; -les figures 4 et 6 représentent schématiquement différents paramètres 20 calculés durant une première phase d'un procédé de calcul ; -les figures 5 et 7 représentent schématiquement différents paramètres calculés durant une phase ultérieure du procédé de calcul; -les figures 8 et 9 représentent schématiquement différents paramètres calculés durant une première itération d'un procédé de calcul selon une autre 25 variante ; -les figures 10 et 11 représentent schématiquement différents paramètres calculés durant des itérations ultérieures d'un procédé de calcul selon cette autre variante ; -les figures 12 et 13 représentent schématiquement différents paramètres 30 calculés durant une première itération d'un procédé de calcul selon encore une autre variante ; -la figure 14 illustre un logigramme d'un exemple de procédé de calcul numérique de la diffraction d'une structure. 35 La figure 1 est une vue en coupe d'un exemple d'une structure 1 à modéliser. La structure 1 présente ici une structure parallélépipédique avec une face supérieure et une face inférieure planes sur lesquelles des ondes électromagnétiques peuvent être incidentes. La structure 1 comporte ici 3033063 7 différentes zones par exemple réalisées en des matériaux différents présentant des répartitions de permittivité diélectrique distinctes. Un modèle numérique de la permittivité diélectrique de la structure 1 en chaque point (pour une résolution donnée du modèle numérique) est connu. La 5 permittivité diélectrique de la structure 1 est donc déterminée en chacun de ses points, ou déterminable par une loi (de répartition par exemple) en chacun de ses points. Un premier aspect de l'invention vise à définir un modèle numérique de la diffraction de la structure 1 pour l'application d'ondes électromagnétiques 10 incidentes sur la face supérieure et/ou la face inférieure. L'invention vise notamment à définir un tel modèle numérique pour permettre un calcul de diffraction dont la quantité de mémoire numérique utilisée est proportionnelle à un nombre M d'ordres de diffraction considérés. 15 L'invention peut être mise en oeuvre par un système de traitement numérique 2 illustré à la figure 3. Un tel système de traitement numérique 2 peut par exemple inclure un dispositif de calcul 21 (par exemple un serveur muni d'un système d'exploitation et d'applications de calcul appropriées), un dispositif de mémorisation 22 du modèle numérique de permittivité diélectrique de la structure 20 1, et un dispositif de mémorisation 23 de résultats de calcul. Le dispositif de mémorisation 23 peut par exemple mémoriser des modèles numériques (détaillés par la suite) de diffraction de différentes couches de la structure 1 ou le modèle numérique de diffraction de l'ensemble de la structure 1. Un exemple de procédé de modélisation numérique selon le premier 25 aspect de l'invention peut être le suivant. Le modèle numérique de permittivité diélectrique de la structure 1 est tout d'abord scindé en un nombre N de modèles numériques de permittivité diélectrique. Chacun de ces modèles numériques de permittivité diélectrique correspond à une couche plane respective de la structure 1, les différentes couches étant superposées selon une direction perpendiculaire 30 aux faces supérieure et inférieure de la structure 1. Chacune de ces couches est identifiée par la suite par son indice i, l'indice i croissant entre la face inférieure et la face supérieure de la structure 1 avec des valeurs comprises entre 1 et N, comme illustré à la figure 1. Dans un souci de simplification, les différentes couches présentent ici une 35 même épaisseur. L'invention pourrait cependant également être appliquée à une scission de la structure 1 en des couches de différentes épaisseurs. Par exemple, la figure 2 est une vue en coupe d'un autre exemple d'une structure 1 à modéliser. La figure 2 correspond à un schéma simplifié d'une structure de type OLED. Les différentes couches du modèle numérique sont ici choisies selon les fonctions des 3033063 8 différentes couches de la structure 1, ces différentes couches présentant typiquement des épaisseurs différentes. Les modèles numériques des différentes couches de la structure 1 pourront donc par exemple inclure un modèle numérique d'une couche d'air 11, un modèle numérique d'une couche de verre 5 12, un modèle numérique d'une couche de diffusion 13, un modèle numérique d'une électrode transparente 14, un modèle numérique d'une couche de polymère 15, un modèle numérique d'une couche émettrice 16, un modèle numérique d'une couche de polymère 17 et un modèle numérique d'une électrode métallique 18. Le nombre de couches N est choisi de sorte que chaque couche est 10 suffisamment fine pour qu'un modèle de diffraction de chaque couche puisse être déterminé à partir de son modèle numérique de permittivité diélectrique, et pour que le modèle de diffraction de chaque couche puisse être calculé avec une occupation de mémoire numérique inférieure ou égale à K*M. Avantageusement, ce nombre de couches N est choisi de sorte que chaque couche est suffisamment 15 fine pour qu'un modèle de diffraction de chaque couche puisse être calculé en un temps inférieur ou égal à K* M * log(M). K est un facteur indépendant de M, typiquement une constante. La détermination du modèle de diffraction de chaque couche peut par exemple être mis en oeuvre par le procédé appelé GSM (pour Generalized Source 20 Method en langue anglaise) par la suite, et décrit dans le document 'New fast and memory-sparing method for rigorous electromagnetic analysis of 2D periodic dielectric structures', Shcherbakov et Tishchenko, publié dans Journal of Quantitative spectroscopy & Radiative Transfer aux pages 158-171. Le modèle de diffraction de chaque couche i peut par exemple être noté 25 comme un opérateur linéaire Li tel que : = L, 30 avec fi et bi+i les amplitudes des M ordres de diffraction des ondes incidentes respectives sur les faces de la couche i, et fi et b, les amplitudes des M ordres diffractés par la couche i. L'opérateur Li peut être obtenu par la méthode des sources généralisées (GSM) décrite dans A. A. Shcherbakov, A. V. Tishchenko, 35 "New fast and memory-sparing method for rigorous electromagnetic analysis of 2D periodic dielectric structures", J. Quant. Spectrosc. Rad. Transfer 113, 158171 (2012) ou par une formulation analytique pour des couches très minces par rapport à la longueur d'onde décrite dans A. V. Tishchenko, "Analytical solutions of 2D grating diffraction: GSM versus Rayleigh hypothesis," Proc. SPIE 5249 p. 3033063 9 683-694 (2004) ou peut même être la matrice S dans le cas de couches simples où S est de type Tiiplitz ou bloc-diagonale. Selon un deuxième aspect de l'invention, on réalise un calcul de diffraction 5 d'au moins une onde incidente, à partir des modèles de diffraction des N couches. Ainsi, après la détermination des modèles de diffraction des N couches, des calculs de diffraction peuvent être réalisés pour différentes ondes incidentes à partir de ces modèles de diffraction. 10 Selon une première variante de calcul de diffraction, on réalise en parallèle un calcul de propagation par application d'une onde incidente sur la couche 1 et un calcul de propagation par application d'une onde incidente sur la couche N. Si l'utilisation de la mémoire numérique peut être encore davantage optimisée qu'avec cette première variante (ici une occupation de mémoire proportionnelle à 15 2M*(N+1)), celle-ci s'avère cependant particulièrement appropriée pour être mise en oeuvre par des moyens de calcul parallèles, par exemple des systèmes à processeurs ou cartes graphiques multiples. Pour le calcul de propagation basé sur l'application d'une onde incidente 20 sur la couche 1, une itération initiale (ou itération d'indice 0) du processus est illustrée en référence à la figure 4. M ordres de diffraction incidents dont les amplitudes sont contenues dans le vecteur fo° sont appliqués au modèle de diffraction de la couche 1, sur sa face 25 externe. Les M ordres de diffraction transmis dont les amplitudes sont contenues dans le vecteur fl° sont calculés avec ce modèle de diffraction, et les M ordres réfléchis bi° sont calculés et mémorisés. Ensuite, pour chaque couche d'indice i compris entre 2 et N, on applique les M ordres transmis calculés f°1 dans le modèle de diffraction de la couche 30 d'indice i. On calcule les M ordres transmis fi°, et on calcule et on mémorise les M ordres réfléchis b°. Pour la couche d'indice N, on mémorise en outre les M ordres transmis fi,°, . Les éléments mémorisés sont illustrés à l'intérieur de cercles en pointillés à la figure 4. Durant cette itération initiale, on applique en pratique l'opérateur Li de la 35 façon suivante, en l'absence d'incidence b11 : = L. 3033063 10 Pour le calcul de propagation basé sur l'application des M ordres incidents sur la couche N, une itération initiale du processus est illustrée en référence à la figure 6. Les ordres incidents cN°,1 sont appliqués au modèle de diffraction de la 5 couche N, sur sa face externe. Les ordres transmis civ° sont calculés avec ce modèle de diffraction, et les ordres réfléchis giv° sont calculés et mémorisés. Ensuite, pour chaque couche d'indice i compris entre N-1 et 1, on applique les ordres transmises calculés c°1 dans le modèle de diffraction de la couche d'indice i. On calcule les ordres transmis c°, et on calcule et on mémorise les 10 ordres réfléchis e. Pour la couche d'indice 1, on mémorise en outre les ordres transmis ci°. Les éléments mémorisés sont illustrés à l'intérieur de cercles en pointillés à la figure 6. Durant cette itération initiale, on applique en pratique l'opérateur Li de la façon suivante, en l'absence d'incidence ei : 15 = L. 20 On calcule les amplitudes de sortie de l'itération initiale, par addition des amplitudes des ordres incidents sur une face externe et des amplitudes des ordres réfléchis par cette même face : f0=fNO±gNO b° = bi° + ci° On note y° un vecteur d'amplitudes de base, formant une solution de base contenant les amplitudes à l'intérieur de la structure : V° = } I/Oi+1 avec tous les i compris entre 1 et N-1. Cette itération initiale nécessite un temps de calcul proportionnel au nombre N de couches et nécessite une occupation mémoire proportionnelle à 2M 35 * (N+1). Cette solution de base y° ne correspond pas à la solution retenue ultérieurement. La solution retenue est déterminée suite à des itérations supplémentaires du calcul de diffraction détaillées par la suite. 25 30 3033063 11 Le vecteur y des amplitudes de la solution finale contenant les amplitudes de tous les ordres de diffraction de toutes les couches est déterminé à partir de la solution de base y° par des itérations successives. Chaque itération est ici suivie d'un test de convergence, pour déterminer si le vecteur de solution de 5 l'itération est une solution finale. Chaque itération peut s'exprimer au moyen d'un opérateur P sous la forme vk = p sk-1 avec k l'indice de l'itération. En principe, la solution recherchée additionne l'ensemble des vecteurs des itérations. La solution vectorielle y recherchée est alors la solution de l'équation implicite y = y° + P.v 10 Cette équation implicite est vraie pour tout système électromagnétique linéaire et résulte de l'application du théorème d'Ewald-Oseen décrit dans le document Born, Max; Wolf, Emil (1999), Principles of Optics (7th ed.), Cambridge: Cambridge University Press, § 2.4. 15 Lors de chaque itération d'indice k>0, on réalise une propagation montante entre les couches 2 et N de chaque vecteur g ,k-11 additionné à chaque vecteur fki , et on réalise une propagation descendante de chaque vecteur bik+11 additionné à chaque vecteur c,k+1, entre les couches N-1 et 1. Ainsi, pour la propagation montante, on applique en pratique l'opérateur Li 20 de la façon suivante, comme illustré à la figure 5: 14.k\ k-1 k gi-1 bk =Li 0 25 avec fik = 0 . On mémorise alors bik et fN . Pour la propagation descendante, on applique en pratique l'opérateur Li de la façon suivante, comme illustré à la figure 7 : 30 = L ( o k-1 k i+1 C i+1 35 avec cNk = 0 . On mémorise alors g ,k et et . Chaque itération nécessite un temps de calcul proportionnel au nombre N de couches et nécessite une occupation mémoire proportionnelle à 2M * (N+1). 3033063 12 La solution vectorielle y finale est la solution de l'équation implicite y = y° + P.v . On réalise un test de convergence à l'issue de chaque itération. Le test de convergence fait partie d'une méthode de résolution d'un système 5 d'équations linéaires défini par l'équation implicite ci-dessus. La résolution de l'équation implicite est par exemple basée sur la méthode GMRES, méthode usuellement utilisée pour obtenir de façon itérative une solution numérique d'un système d'équations linéaires sans inversion de matrice. La méthode GMRES est notamment décrite dans le document 'Iterative Methods for 10 Sparse Linear Systems' de Monsieur Saad, deuxième édition de 'Society for Industrial and Applied Mathematics', publié en 2003 (ISBN 978-0-89871-534-7). La méthode GMRES cherche usuellement à résoudre le système d'équations linéaires de type: 15 A.x = b La matrice A est supposée inversible et de taille (m x m). De plus, on suppose que b est normé, i.e., 'pli= 1, 11.11 représentant ici la norme euclidienne. Le n-ième espace de Krylov pour ce problème est défini ainsi : Kn = {Vect}{ b, A.b, A2.b, ..., Aln-11.b } 20 Où Vect correspond au sous-espace vectoriel engendré. La méthode GMRES donne une approximation de la solution exacte de A.x = b par le vecteur xn E Kn qui minimise la norme du résidu : IIA.xn - b11. Pour garantir le caractère linéairement indépendant aux vecteurs b, A.b, ..., 25 An-1.b, on utilise la méthode d'Arnoldi pour trouver des vecteurs orthonormaux qi, q2,... , qn qui constituent une base de Kn. Ainsi, le vecteur xn E Kn peut s'écrire xn = Qnyn avec yn E Rn, et Qn une matrice de taille (m x n) formée des qi, q2,... , qn. La méthode d'Arnoldi produit aussi une matrice de Hessenberg supérieure 30 fl de taille (n+1).xn avec A.Qn = Qn+1. fin Comme Qn est orthogonale, on a 11A.xn - bll=11fin .yn -I3.e-i I I 35 Où ei = (1,0,0,... ,0) est le premier vecteur de la base canonique de Rn+-1, et p = IIb-A.xoII, avec xo un vecteur d'initialisation (par exemple nul). Ainsi, xn peut être trouvé en minimisant la norme du résidu rn = I. yn - 13.e1 3033063 13 Chaque itération de l'algorithme de la méthode GMRES inclut usuellement effectuer une étape de l'algorithme d'Arnoldi ; trouver yn qui minimise lirnil ; 5 calculer xn = Qn..yn ; recommencer tant que le résidu est plus grand qu'une quantité (dite tolérance) choisie arbitrairement au début de l'algorithme. Pour l'application de la méthode GMRES au test de convergence de 10 l'invention, on substitue l'opération (x-P.x) à l'opération de multiplication A.x mentionnée précédemment. L'espace vectoriel Wk = (vo v1 vk ), V =1 de Krylov est formé en utilisant 15 V 0 =-v-v Ilv°011 ° L'astérisque * après un vecteur désigne sa transposition avec la conjugaison. L'application d'un pas itératif Vk - P.vk permet de créer un nouveau 20 vecteurv k+1 et d'agrandir l'espace de Krylov jusqu'à Wk_ri L'opération de normalisation des vecteurs vk définit la matrice de Hessenberg Hl k 25 =1 P .V k), i 0, on additionne les amplitudes fN et Cik avec le vecteur Vk pour trouver un vecteur agrandi. La solution finale trouvée 20 ( g, c, 3033063 16 par une telle méthode itérative comprend les amplitudes fN et el . Les amplitudes f et b à la sortie de la structure sont trouvées par addition avec les amplitudes réfléchies de l'itération initiale : 5 g Ifb} Ib:No }±{fci Cette modification résulte en une occupation de mémoire plus grande, proportionnelle à (N+1) pour chaque itération, désormais on n'a pas besoin de l'itération finale, le nombre totale des itérations est donc diminué. 10 Selon une deuxième variante, on réalise séquentiellement l'application d'une propagation dans un sens, puis on applique les ordres réfléchis dans le sens opposé. Cette seconde variante permet d'optimiser l'utilisation de la mémoire numérique avec une occupation proportionnelle à M*(N+1). 15 Lors d'une itération initiale (ou itération d'indice 0) du processus illustrée en référence à la figure 8, on réalise par exemple une propagation des ordres de diffraction incidents sur la couche 1. Le vecteur des ordres incidents fo° est appliqué au modèle de diffraction de la couche 1, sur sa face externe. Les ordres transmis fl° sont calculés avec 20 ce modèle de diffraction, et les ordres réfléchis bl° sont calculés et mémorisés. Ensuite, pour chaque couche d'indice i compris entre 2 et N, on applique les ordres transmis calculés f°1 dans le modèle de diffraction de la couche d'indice i. On calcule les ordres transmis fi°, et on calcule et on mémorise les amplitudes des ordres réfléchis bi°. Pour la couche d'indice N, on mémorise en 25 outre les amplitudes des ordres transmis f'°, . Les éléments mémorisés sont illustrés à l'intérieur de cercles en pointillés à la figure 8. Durant cette itération initiale, on applique en pratique l'opérateur Li de la façon suivante, en l'absence d'incidence bi°,1 : 30 = L'itération initiale se poursuit par l'application des ordres incidents cN°,1 au modèle de diffraction de la couche N, sur sa face externe. Les ordres transmis c°N 35 sont calculés avec ce modèle diffraction, et les ordres réfléchi giv° sont calculés et mémorisés. Ensuite, pour chaque couche d'indice i compris entre N-1 et 1, on additionne les ordres transmis calculés c,°,1 et les ordres réfléchis calculés précédemment b,°,1 dans le modèle de diffraction de la couche d'indice i. 3033063 17 On calcule les amplitude des ordres transmis c,°, et on calcule et on mémorise les amplitudes des ordres réfléchis e . Pour la couche d'indice 1, on mémorise en outre les amplitudes des ordres transmis ci°. Les éléments mémorisés sont illustrés à l'intérieur de cercles en pointillés à la figure 9. Afin 5 d'optimiser la quantité de mémoire numérique utilisée, après le calcul d'un vecteur des ordres transmis c° et d'un vecteur e , on peut libérer la mémoire numérique occupée par les amplitudes des ordres bi°,1 calculées précédemment. Durant cette itération initiale, on applique en pratique l'opérateur Li de la façon suivante, en l'absence d'incidence e 1 : 10 ( 0 \Cio+1 bi°+1) 15 On calcule les amplitudes de sortie de l'itération initiale, par addition entre les amplitudes des ordres incidents sur une face externe et les amplitudes des ordres réfléchis par cette face : f 0 = fNO gNO b° +c° On note y° un vecteur de base, formant une solution de base avec les amplitudes à l'intérieur de la structure : v° = avec tous les i compris entre 1 et N-1. Cette itération initiale nécessite un temps de calcul proportionnel au nombre N de couches. Du fait de la libération de la mémoire numérique occupée 30 par les ordres calculés bi°,1, cette itération initiale nécessite une occupation mémoire proportionnelle à M * (N+1). Cette solution de base y° ne correspond pas à la solution retenue ultérieurement. La solution retenue est déterminée suite à des itérations 35 supplémentaires détaillées par la suite. Comme pour la première variante, le vecteur y de solution finale est déterminé à partir de la solution de base y° par des itérations successives. Chaque itération est suivie d'un test de convergence tel que détaillé 20 25 3033063 18 précédemment, pour déterminer si le vecteur de solution de l'itération est une solution finale. Chaque itération peut s'exprimer au moyen d'un opérateur P sous la forme vk = p sk-1 avec k l'indice de l'itération. La solution vectorielle y est la solution 5 de l'équation implicite y = y° + P.v Lors de chaque itération d'indice k>0, on réalise une propagation montante entre les couches 2 et N de chaque vecteur fiki , en prenant fik =0, additionné à chaque vecteur e-11 et on mémorise chaque vecteur bik résultant. En suite, on 10 réalise une propagation descendante de chaque vecteur ck,1 additionné à chaque vecteur bik+1, entre les couches N-1 et 1, en prenant cNk =0, et on mémorise chaque vecteur e résultant. Ainsi, pour la propagation montante, on applique en pratique l'opérateur Li de la façon suivante, comme illustré à la figure 10: 15 ( k-1k g i-1 + fi-1 0 avec fik =0 . 20 On mémorise alors bik . Après chaque calcul de bik , on libère la mémoire numérique occupée par k-I g i-I - Pour la propagation descendante, on applique en pratique l'opérateur Li de 25 la façon suivante, comme illustré à la figure 11 : ( ( o = L. g k VCi ) i+I C +1i 30 avec C Nk = 0 . On mémorise alors e . Après chaque calcul de e on libère la mémoire numérique occupée par Chaque itération nécessite un temps de calcul proportionnel au nombre N 35 de couches et nécessite une occupation mémoire proportionnelle à M * (N+1). Comme pour la première variante, on réalise un test de convergence à l'issue de chaque itération pour résoudre l'équation implicite y = y° + P.v . Le test de convergence utilise également une méthode de résolution d'un système = L. 3033063 19 d'équations linéaires pour rechercher la solution de cette équation implicite sans inversion de matrice. Lors de l'itération finale, on réalise une propagation montante entre les 5 couches 2 et N de chaque vecteur g, 1 additionné à chaque vecteur fi, et on réalise une propagation descendante de chaque vecteur bi+1 additionné à chaque vecteur ci+1, entre les couches N-1 et 1. Pour la propagation montante, on applique en pratique l'opérateur Li de la façon suivante, comme illustré à la figure 10: 10 ( ( g + f =Li 0 ) avec f1 =0. 15 On mémorise alors bi. Après chaque calcul de bi , on libère la mémoire numérique occupée par - 20 Pour la propagation descendante, on applique en pratique l'opérateur Li de la façon suivante, comme illustré à la figure 11 : ( ( 0 =L, bi+ +ci+1) 25 avec cN =0. On mémorise alors gi. Après chaque calcul de gi, on libère la mémoire numérique occupée par 30 Cette itération finale nécessite un temps de calcul proportionnel au nombre N de couches et nécessite une occupation mémoire proportionnelle à M * (N+1). On trouve en suite la solution vectorielle à la sortie de la structure IfRe±e}±{fi 35 b bi° + Selon une modification de l'algorithme, le vecteur de base y° est agrandi en comprenant les amplitudes f° et c°. Lors de chaque itération d'indice k>0, on additionne les amplitudes fNk et et avec le vecteur vk pour trouver un vecteur 3033063 20 élargi. Les amplitudes f et b à la sortie de la structure sont trouvées par addition avec les amplitudes réfléchies de l'itération initiale : g Ifb} Ib:No }±{ fN Cl Cette modification résulte en une occupation de mémoire plus grande, proportionnelle à (N+1) pour chaque itération, désormais on n'a pas besoin de l'itération finale, le nombre totale des itérations est donc diminué. 10 Selon une troisième variante, on réalise en parallèle un calcul de propagation par application d'une onde incidente sur la couche 1 et un calcul de propagation par application d'une onde incidente sur la couche N. L'utilisation de la mémoire numérique reste la même que dans la première variante (une occupation de mémoire proportionnelle à 2M*(N+1)), le calcul de propagation 15 dans les couches est davantage parallélisable, ce qui s'avère particulièrement approprié pour être mis en oeuvre pour un grand nombre de moyens de calcul parallèles (par exemple au moins 4, ou avec un nombre de moyens de calcul parallèles au moins égal au quart du nombre de couches N), par exemple des systèmes à processeurs ou cartes graphiques multiples. Pour simplifier, on va 20 maintenant détailler un exemple dans lequel on réalise autant de calculs en parallèle que de couches N. Pour le calcul de propagation basé sur l'application d'une onde incidente sur la couche 1, une itération initiale (ou itération d'indice 0) du processus est 25 illustrée en référence à la figure 12. M ordres de diffraction incidents dont les amplitudes sont contenues dans le vecteur fo° sont appliqués au modèle de diffraction de la couche 1, sur sa face externe. Les M ordres de diffraction transmis dont les amplitudes sont contenues dans le vecteur fi° sont calculés avec ce modèle de diffraction et mémorisés, et 30 les M ordres réfléchis bi° sont calculés et mémorisés : (fo b° = \ 1) Pour cette itération initiale on ne calcule que les ordres transmis bi°, et les 35 ordres réfléchis gZ. On suppose pour toutes les autres amplitudes des ordres transmis : bi° = 0 pour tout i compris entre 2 et N. 5 3033063 21 Pour le calcul de propagation basé sur l'application des M ordres incidents sur la couche N, une opération initiale du processus est illustrée en référence à la figure 13. Les ordres incidents cN°,1 sont appliqués au modèle de diffraction de la 5 couche N, sur sa face externe. Les ordres transmis civ° sont calculés avec ce modèle de diffraction et mémorisés, et les ordres réfléchis giv° sont calculés et mémorisés. ( 0 ( 0 gN 0 - 0 \CN \CN+1) On suppose pour toutes les autres amplitudes des ordres réfléchis : ) g (i =0 pour tout i compris entre 1 et N-1 15 On note y° un vecteur de base, formant une solution de base avec les amplitudes f °et b° 20 pour tout i compris entre 1 et N-1. Cette itération initiale nécessite un temps de calcul d'une couche puisque le calcul est parallélisé et nécessite une occupation mémoire proportionnelle à 4M. Pour trouver la solution exacte du problème de diffraction, des itérations 25 sont nécessaires et sont détaillées par la suite. Comme pour la première et la seconde variante, le vecteur v des amplitudes de la solution finale est déterminé à partir de la solution de base y° par des itérations successives. Chaque itération est suivie d'un test de convergence tel que détaillé précédemment pour déterminer si le vecteur de 30 solution de l'itération peut être considéré comme une solution finale. Chaque itération peut s'exprimer au moyen d'un opérateur P sous la forme vk = p sk-1 avec k l'indice de l'itération. La solution vectorielle v est la solution de l'équation implicite v = v° +P.v 35 Lors de chaque itération d'indice k>0, on réalise une propagation de chaque vecteur additionné à chaque vecteur bik+i vers une couche i. Un tel calcul est réalisé parallèlement pour chaque couche d'indice i entre 1 et N : 10 Vo =le} bi°,1 3033063 22 ( ( gi bk ) bik+-11) avec go" =0 et bNk-,11 = 0 . Chaque itération nécessite un temps de calcul d'une couche et nécessite 5 une occupation mémoire proportionnelle à 2M * (N-1). Comme pour les deux premières variantes, on réalise un test de convergence à l'issue de chaque itération pour résoudre l'équation implicite = y° +P.v . Le test de convergence utilise également une méthode de résolution 10 d'un système d'équations linéaires pour rechercher la solution de cette équation implicite sans inversion de matrice. Lors de l'itération finale, on réalise une propagation de vecteur g N 1 vers la couche N : g N-1 ° et une propagation de vecteur c1 vers la couche 1 : 20 = LN 1) Un tel calcul est réalisé parallèlement pour les deux couches 1 et N avec k1 go =0 et bNk-,11 = . Cette itération nécessite un temps de calcul d'une couche et nécessite une 25 occupation mémoire proportionnelle à 2M. La solution vectorielle à la sortie de la structure est trouvée comme la somme 30 b b1° +b1 Selon une modification de l'algorithme, le vecteur de base y° est agrandi en comprenant aussi les amplitudes f° et c°. Lors de chaque itération d'indice k>0, on additionne les amplitudes fNk et et avec le vecteur vk pour trouver un 15 ( = LN N ) 3033063 23 vecteur agrandi. Les amplitudes f et b à la sortie de la structure sont trouvées par addition avec les amplitudes réfléchies de l'itération initiale : g Ifb}-{b:No }±{ fN Cl Cette modification résulte en une occupation de mémoire plus grande, proportionnelle à (N+1) pour chaque itération, désormais on n'a pas besoin de l'itération finale, le nombre totale des itérations est donc diminué. 10 La figure 14 représente schématiquement un exemple de suite d'étapes mises en oeuvre dans un procédé de calcul numérique de la diffraction d'une structure. Lors de l'étape 301, on définit un modèle numérique de la permittivité diélectrique d'une structure en chacun de ses points. 15 Lors de l'étape 302, le modèle numérique de permittivité diélectrique de la structure est scindé en N modèles numériques de permittivité diélectrique, chacun de ces modèles numériques correspondant chacun à une couche plane de la structure, ces couches planes étant superposées selon une direction perpendiculaire aux faces supérieure et inférieure de la structure. Les couches, 20 correspondant à la scission du modèle numérique de permittivité diélectrique de la structure, sont déterminées de sorte que le modèle de diffraction de chacune de ces couches puisse être calculé à partir de son modèle numérique de permittivité diélectrique en un temps inférieur ou égal à K* M * log(M). Lors de l'étape 303, on détermine un modèle de diffraction pour chacune 25 des N couches à partir de son modèle numérique de permittivité diélectrique. Le modèle de diffraction de chaque couche i peut par exemple être noté comme un opérateur Li tel que : 30 ( ft = L, avec fi et bi+i les amplitudes incidentes sur les faces de la couche i, et f et b, les amplitudes diffractées par la couche i. Lors de l'étape 304, on réalise une itération initiale de propagation des ordres incidents à travers les modèles de diffraction des couches de la structure. 35 Lors de l'étape 305, on réalise une itération de propagation des ordres diffractés et calculés lors de l'itération précédente à travers les modèles de diffraction des couches de la structure. Lors de l'étape 306, on réalise un test de convergence de la dernière itération de propagation d'ordres diffractés en utilisant une méthode de résolution 5 3033063 24 d'un système d'équations linéaires sans inversion de matrice et en l'appliquant à la résolution de l'équation implicite recherchée pour la convergence des itérations. Si la condition du test de convergence n'est pas remplie, on exécute à nouveau l'étape 305 pour une nouvelle itération de propagation des ordres diffractés. 5 Lors de l'étape 307, on réalise une ultime itération de propagation des ordres diffractés et calculés lors des itérations précédentes, à travers les modèles de diffraction des couches de la structure. On en déduit le résultat de la diffraction de la structure. 10 Le procédé de calcul numérique de la diffraction d'une structure détaillé dans les exemples qui précèdent peut être appliqué dans différents domaines techniques, dont une liste non-exhaustive est présentée par la suite à titre d'exemples. 15 En optique diffractante, pour la conception de structures comportant des micro- et nano-structures (micro et nano-structures dont la dimension caractéristique est plus petite que 3 à 5 fois la longueur d'onde), les propriétés de diffraction dépendent de la polarisation et les méthodes scalaires usuelles les modélisent de façon erronée. Les procédés de calcul selon l'invention permettent 20 de simuler de larges sections de telles structures périodiques ou non-périodiques, voire la structure complète. Un tel procédé de calcul peut par exemple être appliqué à une micro-lentille réfracto-diffractante de dimension minimale et de poids minimal utilisée dans un système de lecture ou d'écriture optique à balayage ultra-rapide. 25 Un procédé de calcul selon l'invention peut également être appliqué en optique diffusante. Un tel procédé de calcul peut en particulier être appliqué pour l'optimisation de la répartition souhaitée de lumière mono- ou poly-chromatique par une couche diffusante. Une telle couche diffusante comprend typiquement un 30 milieu hôte contenant des microsphères ou micro-polyèdres d'indice de réfraction différent de celui du milieu hôte. La couche diffusante illumine généralement un écran de façon uniforme, éclairant de façon désirée un volume à partir d'une source de lumière localisée ou faisant partie d'une paroi lumineuse. En permettant de modéliser exactement de larges sections de telles couches diffusantes, le 35 procédé de calcul décrit prend en compte les effets de diffraction multiple dans une couche et les effets de cohérence. Un tel procédé de calcul peut également être utilisé en optique diffusante pour l'optimisation d'une couche diffusante destinée à extraire efficacement la 3033063 25 lumière produite, par exemple, par une source de lumière du type Diode ElectroLuminescente (LED en anglais) ou d'une diode électroluminescente organique. Un procédé de calcul selon l'invention peut également être appliqué en 5 microélectronique pour l'optimisation de différents processus de photolithographie, en particulier lorsque la dimension caractéristique des motifs au niveau du réticule (le masque) et/ou au niveau de la plaque de silicium est de l'ordre de grandeur ou plus petite que la longueur d'onde de projection. C'est en particulier le cas pour les noeuds technologiques de 45, de 30 nm et en dessous 10 aux longueurs d'onde des lasers à excimère KrF (248 nm) et ArF (193 nm). Un des processus à effectuer au niveau du réticule est généralement la correction de proximité optique (OPC ou Optical Proximity Correction en langue anglaise). Une telle correction OPC suppose le calcul exact de la transmission du réticule sous incidence optique variable, le réticule incluant des motifs de différentes 15 profondeurs de matériaux tels que SiO2, chrome, MoSiO, TaO. Les méthodes de calcul de transmission usuelles pour de telles structures sont les méthodes scalaires avec correctifs euristiques inspirés de calculs exacts locaux. L'usage de ces correctifs permet de conserver la rapidité de calcul des méthodes scalaires mais leur domaine de validité est limité. Un procédé de calcul 20 selon l'invention permet de repousser très loin la frontière entre le domaine des structures modélisables avec exactitude et celui des structures calculables approximativement. Un autre processus en photolithographie en microélectronique est celui de la modélisation de l'image latente dans une couche de résine photosensible 25 déposée sur un substrat. La topographie et la composition de la surface résulte d'un certain nombre de processus technologiques précédents, ce qui affecte notablement la distribution de puissance lumineuse incidente par diffraction en réflexion. Les méthodes scalaires de l'état de la technique sont confrontées à de grandes difficultés car la structure en dessous de la couche de photorésine où se 30 forme l'image latente projetée peut être très épaisse et très complexe. Un procédé de calcul selon l'invention permet de prendre en compte la structure réelle mise en oeuvre. REVENDICATIONS1. Procédé de calcul numérique de la diffraction d'une structure (1) dont un modèle numérique définit la permittivité diélectrique en chaque point, comprenant les étapes de : -définition d'une modélisation numérique d'un modèle de diffraction d'une structure, incluant : -la scission du modèle numérique de permittivité diélectrique de la structure en un nombre N de modèles numériques de permittivité diélectrique de couches planes respectives superposées selon une direction et ordonnées selon un indice i compris entre 1 et N selon ladite direction ; -la détermination du modèle de diffraction de chacune desdites couches à partir du modèle numérique de permittivité diélectrique de cette couche, le nombre N de couches étant déterminé pour que l'occupation de mémoire numérique pour le calcul de diffraction par le modèle de diffraction de chaque couche soit inférieur à K* M, avec M un nombre d'ordres de diffraction de cette couche et K un facteur indépendant de M; -calcul numérique de la diffraction de la structure, incluant : -une itération initiale durant laquelle on réalise : -une propagation d'onde dans le sens des couches 1 vers N, incluant : -le calcul d'une onde réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice 1 ; -pour une couche d'indice i compris entre 2 et N, calcul d'une onde réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice i, cette onde incidente étant l'onde transmise calculée pour la couche d'indice i-1 ; -mémorisation des ondes réfléchies calculées pour chacune des couches d'indice i et de l'onde transmise calculée pour la couche d'indice N ; -et une propagation d'onde dans le sens des couches N vers 1, incluant : -une onde réfléchie et une onde transmise sont calculées par application d'une onde incidente au modèle de diffraction de la couche d'indice N ; -pour une couche d'indice i compris entre N-1 et 1, une onde réfléchie et une onde transmise sont calculées par application d'une onde incidente au modèle de diffraction de la couche d'indice i, cette onde incidente incluant l'onde transmise calculée pour la couche d'indice i+1 ; 3033063 27 -mémorisation des ondes réfléchies calculées pour chacune des couches d'indice i et de l'onde transmise calculée pour la couche d'indice N ; -des itérations ultérieures, chaque itération ultérieure d'indice k incluant : 5 -une propagation d'onde dans le sens des couches 2 vers N, incluant : -pour une couche d'indice i compris entre 2 et N, calcul d'une onde réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice 10 i, cette onde incidente incluant l'onde transmise calculée pour la couche d'indice i-1 durant cette itération et incluant l'onde réfléchie calculée pour la couche d'indice i-1 durant la dernière propagation dans le sens opposé; -mémorisation des ondes réfléchies calculées pour chacune 15 des couches d'indice i et de l'onde transmise calculée pour la couche d'indice N ; -une propagation d'onde dans le sens des couches N-1 vers 1, incluant : -pour une couche d'indice i compris entre N-1 et 1, calcul 20 d'une onde réfléchie et d'une onde transmise par application d'une onde incidente au modèle de diffraction de la couche d'indice i, cette onde incidente incluant l'onde transmise calculée pour la couche d'indice i+1 durant cette itération et incluant l'onde réfléchie calculée pour la couche d'indice i+1 25 durant la dernière propagation dans le sens opposé; -mémorisation des ondes réfléchies calculées pour chacune des couches d'indice i et de l'onde transmise calculée pour la couche d'indice 1 ; -un test de convergence de la solution Vk formée des ondes 30 réfléchies et transmises calculées lors de l'itération d'indice k, par application d'une méthode itérative de résolution de systèmes d'équations linéaires à la résolution de l'équation implicite de convergence Vk=P.Vk, avec P un opérateur correspondant à la transformation de la solution de toute solution Vk-1 en une solution Vk 35 à l'itération k. 2. Procédé de calcul numérique selon la revendication 1, dans lequel ladite méthode itérative de résolution de systèmes d'équations linéaires est une méthode de type GMRES. 40 3033063 28 3. Procédé de calcul numérique selon la revendication 1, dans lequel ladite méthode itérative de résolution de systèmes d'équations linéaires est une méthode du gradient biconjugué stabilisée. 5 4. Procédé de calcul numérique selon l'une quelconque des revendications précédentes, dans lequel ladite détermination du modèle de diffraction de chacune desdites couche est mis en oeuvre par un procédé GSM. 5. Procédé de calcul numérique selon l'une quelconque des revendications 10 précédentes, dans lequel N est au moins égal à 5. 6. Procédé de calcul numérique selon l'une quelconque des revendications précédentes, dans lequel le calcul de propagation dans le sens des couches 2 à N et le calcul de propagation dans le sens des couches N-1 à 1 sont réalisés 15 en parallèle. 7. Procédé de calcul numérique selon l'une quelconque des revendications 1 à 5, dans lequel le calcul de propagation dans le sens des couches 2 à N et le calcul de propagation dans le sens des couches N-1 à 1 sont réalisés 20 séquentiellement. 8. Procédé de calcul numérique selon l'une quelconque des revendications 1 à 5, dans lequel le calcul de propagation comprend des calculs en parallèle pour des couches différentes. 25