diff --git a/Debug_JINTRACMesh.svg b/Debug_JINTRACMesh.svg new file mode 100644 index 000000000..890803549 --- /dev/null +++ b/Debug_JINTRACMesh.svg @@ -0,0 +1,21 @@ + + + + + + + + + + + + + diff --git a/Notes_Upgrades/Communications/PyconFR2019/2019_LM_PyconFR.tex b/Notes_Upgrades/Communications/PyconFR2019/2019_LM_PyconFR.tex index df05ec538..8b9e983bf 100644 --- a/Notes_Upgrades/Communications/PyconFR2019/2019_LM_PyconFR.tex +++ b/Notes_Upgrades/Communications/PyconFR2019/2019_LM_PyconFR.tex @@ -50,7 +50,7 @@ %{ % \leavevmode% % \hbox{% -% \begin{beamercolorbox}[wd=.6\paperwidth,ht=2.25ex,dp=1ex,center]{institute in head/foot}% +% \begin{beamercolorbox}[wd=.6\paperwidth,ht=2.25ex,dp=1ex,center]{institute in head/foot}% % \usebeamerfont{institute in head/foot}\insertshortinstitute % \end{beamercolorbox}% % \begin{beamercolorbox}[wd=.4\paperwidth,ht=2.25ex,dp=1ex,center]{date in head/foot}% @@ -94,7 +94,7 @@ \inst{1}% INRIA Grand-Est, TONUS Team, Strasbourg, France\\ - + \inst{2}% CEA, Cadarache, France @@ -119,7 +119,7 @@ %\usebackgroundtemplate{\includegraphics[height=\paperheight,width=\paperwidth]{figures/background.png}}% \begin{frame} \titlepage -\end{frame} +\end{frame} % ================================= % ================================= @@ -162,27 +162,27 @@ \section{Context} % ================================= \begin{frame} \frametitle{Context: Controlled fusion and magnetic confinement} - + \begin{columns} - \begin{column}{0.65\textwidth} - \begin{center} + \begin{column}{0.65\textwidth} + \begin{center} \textbf{D-T Fusion reaction\;\;\;}\vspace{0.3cm} %\includegraphics[width=3.5cm]{figures/DT_reaction3.png} - \resizebox{0.5\textwidth}{!}{\input{figures/reaction.tex}} - \end{center} - \vspace{-0.5cm} + \resizebox{0.5\textwidth}{!}{\input{figures/reaction.tex}} + \end{center} + \vspace{-0.5cm} \begin{itemize}%[leftmargin=1em] - \item Gas $>$ 100 Million\degree K composed of positive ions and negative electrons: plasma - \item Confinement using electromagnetic fields + \item Gas $>$ 100 Million\degree K composed of positive ions and negative electrons: plasma + \item Confinement using electromagnetic fields \item break-even not obtained yet \end{itemize} - \end{column} - - \begin{column}{0.35\textwidth} - \includegraphics[width=1.\textwidth]{figures/ITER.pdf} - \end{column} - + \end{column} + + \begin{column}{0.35\textwidth} + \includegraphics[width=1.\textwidth]{figures/ITER.pdf} + \end{column} + \end{columns} %\begin{itemize}%[leftmargin=1em] @@ -204,8 +204,8 @@ \section{Tomography diagnostics} % ================================= \begin{frame} \frametitle{Tokamak diagnostics to measure plasma quantities} - - \metroset{block=fill} + +\metroset{block=fill} \begin{alertblock}{Diagnostics} Set of instruments to measure plasma quantities, for understanding, control, optimization. \end{alertblock} @@ -214,7 +214,7 @@ \section{Tomography diagnostics} $\quad \Rightarrow$ cameras (1D or 2D) for measuring light in various wavelengths - + % %\begin{itemize} @@ -239,20 +239,20 @@ \section{Tomography diagnostics} %% ================================= %\begin{frame} %\frametitle{A tokamak as a poloidal + horizontal projections} -% +% %\begin{columns} -% +% % \begin{column}{0.25\textwidth} % \includegraphics[width=1.\textwidth]{figures/ITER.pdf} % \end{column} -% +% % \begin{column}{0.75\textwidth} % \begin{center} % \textbf{tofu.geom.Config class} % \end{center} % \includegraphics[width=1.\textwidth]{figures/Config01.pdf} -% \end{column} -% +% \end{column} +% %\end{columns} % % @@ -268,41 +268,41 @@ \section{Tomography diagnostics} $$M_i(t) = \iiint\limits_{V_i} \vv{\varepsilon(x,t)}\cdot\vv{n} \,\Omega_i \;{d}V$$ \vspace{-0.75cm} \begin{columns} - \begin{column}{0.325\textwidth} - - %\def\svgwidth{\linewidth} - %\import{figures/}{detectors.pdf_tex} - \includegraphics[width=\linewidth,trim={0 0 1cm 0}]{figures/Tomography.pdf} - - \end{column} - - \begin{column}{0.675\textwidth} - \begin{center} - - %DIRECT PROBLEM - \begin{block}{} + \begin{column}{0.325\textwidth} + + %\def\svgwidth{\linewidth} + %\import{figures/}{detectors.pdf_tex} + \includegraphics[width=\linewidth,trim={0 0 1cm 0}]{figures/Tomography.pdf} + + \end{column} + + \begin{column}{0.675\textwidth} + \begin{center} + + %DIRECT PROBLEM + \begin{block}{} \begin{itemize} \item \textcolor{myblue}{\textbf{Direct problem} (synthetic diagnostic):\\ } \end{itemize} \end{block} - %\textbf{Direct problem} (synthetic diagnostic):\\ + %\textbf{Direct problem} (synthetic diagnostic):\\ Simulated emissivity $\longrightarrow$ integrated measurements\\ %INVERSE PROBLEM - \begin{block}{} + \begin{block}{} \begin{itemize} \item \textcolor{myblue}{\textbf{Inverse problem} (tomography):\\ } \end{itemize} \end{block} - %\textbf{Direct problem} (synthetic diagnostic):\\ + %\textbf{Direct problem} (synthetic diagnostic):\\ Integrated measurements $\longrightarrow$ Reconstructed emissivity \\ - \end{center} - - \end{column} -% + \end{center} + + \end{column} +% \end{columns} @@ -319,44 +319,44 @@ \section{Tomography diagnostics} $$M_i(t) = \iiint\limits_{V_i} \vv{\varepsilon(x,t)}\cdot\vv{n} \,\Omega_i \;{d}V$$ \vspace{-1cm} \begin{columns} - \begin{column}{0.325\textwidth} - - %\def\svgwidth{\linewidth} - %\import{figures/}{detectors.pdf_tex} - \includegraphics[width=\linewidth,trim={0 0 1cm 0}]{figures/Tomography.pdf} - - \end{column} - - \begin{column}{0.675\textwidth} - \begin{center} - - %DIRECT PROBLEM - \begin{block}{} + \begin{column}{0.325\textwidth} + + %\def\svgwidth{\linewidth} + %\import{figures/}{detectors.pdf_tex} + \includegraphics[width=\linewidth,trim={0 0 1cm 0}]{figures/Tomography.pdf} + + \end{column} + + \begin{column}{0.675\textwidth} + \begin{center} + + %DIRECT PROBLEM + \begin{block}{} \begin{itemize} \item \textcolor{myblue}{\textbf{Direct problem} (synthetic diagnostic):\\ } \end{itemize} \end{block} - %\textbf{Direct problem} (synthetic diagnostic):\\ + %\textbf{Direct problem} (synthetic diagnostic):\\ Simulated emissivity $\longrightarrow$ measurements\\ \alert{Spatial integration} \vspace{-0.5cm} %INVERSE PROBLEM - \begin{block}{} + \begin{block}{} \begin{itemize} \item \textcolor{myblue}{\textbf{Inverse problem} (tomography):\\ } \end{itemize} \end{block} - %\textbf{Direct problem} (synthetic diagnostic):\\ + %\textbf{Direct problem} (synthetic diagnostic):\\ Integrated measurements $\longrightarrow$ Reconstructed emissivity \\ \alert{Mesh and basis functions construction, spatial integration, data filtering, inversion routines, etc.} - - \end{center} - - \end{column} -% + + \end{center} + + \end{column} +% \end{columns} \pause @@ -393,7 +393,7 @@ \section{The ToFu package} \item low traceability, reproducibility \item low standardization, unclear assumptions / methods \end{itemize} - + \end{frame} % ================================= @@ -415,7 +415,7 @@ \section{The ToFu package} \begin{column}{0.35\textwidth} \vspace{-1cm} \begin{center} - \includegraphics[width=0.8\linewidth]{figures/ci.pdf} + \includegraphics[width=0.8\linewidth]{figures/ci.pdf} \end{center} \end{column} \end{columns} @@ -425,9 +425,9 @@ \section{The ToFu package} \begin{center} \alert{\textbf{ToFu}\footnote{repository: \url{https://github.com/ToFuProject/tofu}}\footnote{documentation: \url{https://tofuproject.github.io/tofu/index.html} }\footcite{didier2016}} = \alert{\textbf{To}mography for \textbf{Fu}sion} \end{center} -\vspace{-5mm} +\vspace{-5mm} \end{block} - + \end{frame} % ================================= @@ -459,12 +459,12 @@ \section{The ToFu package} \end{column} \begin{column}{0.35\textwidth} \begin{center} - \includegraphics[width=\textwidth]{figures/qr-code.png} + \includegraphics[width=\textwidth]{figures/qr-code.png} \end{center} \end{column} \end{columns} \begin{center} - \includegraphics[width=\textwidth]{figures/badges.png} + \includegraphics[width=\textwidth]{figures/badges.png} \end{center} \end{frame} @@ -475,9 +475,9 @@ \section{The ToFu package} \frametitle{Tofu's structure} \begin{center} - \includegraphics[width=0.8\linewidth]{figures/tofu.pdf} + \includegraphics[width=0.8\linewidth]{figures/tofu.pdf} \end{center} - + \end{frame} % ================================= @@ -487,9 +487,9 @@ \section{The ToFu package} \frametitle{Tofu's geometry representation} \begin{center} - \includegraphics[width=0.6\linewidth]{figures/ITER.pdf} + \includegraphics[width=0.6\linewidth]{figures/ITER.pdf} \end{center} - + \end{frame} % ================================= @@ -498,9 +498,9 @@ \section{The ToFu package} \frametitle{Tofu's geometry representation} \begin{center} - \includegraphics[width=\linewidth]{figures/newiter.png} + \includegraphics[width=\linewidth]{figures/newiter.png} \end{center} - + \end{frame} % ================================= @@ -528,7 +528,7 @@ \section{The ToFu package} %\frametitle{What ToFu can do: 3D modeling of a 2D camera} % \begin{center} % \includegraphics[width=\textwidth]{figures/cam2d.png} -% \end{center} +% \end{center} %\end{frame} %% ================================= @@ -555,7 +555,7 @@ \section{The ToFu package} \begin{center} \textbf{2D Camera\\(tofu.geom.CamLOS2D)} \includegraphics[width=\textwidth]{figures/cam2d_bis.pdf} - \end{center} + \end{center} \end{column} \end{columns} \end{frame} @@ -582,7 +582,7 @@ \section{The ToFu package} \begin{center} \textbf{2D Camera with reflexions\\(tofu.geom.CamLOS2D)} \includegraphics[width=\textwidth]{figures/cam2d_bis_ref.pdf} - \end{center} + \end{center} \end{column} \end{columns} \end{frame} @@ -591,7 +591,7 @@ \section{The ToFu package} %% ================================= %\begin{frame} %\frametitle{What ToFu can do: computing synthetic signals} -% +% %\end{frame} % ================================= @@ -636,9 +636,9 @@ \section{Code Optimization} \item Geometry defined with minimal data polygon $(R,Z)$\\ extruded along $\varphi$ \item Symmetry of vessel along $\varphi$ - + \end{itemize} - + \vspace{0.5cm} \begin{center} \includegraphics[height=2.3cm]{figures/confB3_wp_view2.png}% @@ -647,7 +647,7 @@ \section{Code Optimization} ~ \includegraphics[height=2.3cm]{figures/confB3_wp_view3.png}% \end{center} - + \end{frame} % ================================= @@ -668,12 +668,12 @@ \section{Code Optimization} \item[{$\Rightarrow$}] Light memory-wise \end{itemize} \item[{$\Rightarrow$}] Equivalent to: set of truncated cones (frustums) of generatrix $A_iB_i$ - + \end{itemize} \end{column} \begin{column}{0.35\textwidth} \begin{center} - \includegraphics[width=\linewidth]{figures/tore_cones1.pdf} + \includegraphics[width=\linewidth]{figures/tore_cones1.pdf} \end{center} \end{column} \end{columns} @@ -691,7 +691,7 @@ \section{Code Optimization} DM = k u \end{array}\right. $$ - + \end{frame} % ================================= @@ -736,9 +736,9 @@ \section{Code Optimization} \frametitle{Tofu's structure} \begin{center} - \includegraphics[width=0.8\linewidth]{figures/tofu.pdf} + \includegraphics[width=0.8\linewidth]{figures/tofu.pdf} \end{center} - + \end{frame} % ================================= @@ -764,7 +764,7 @@ \section{Code Optimization} \item hybrid: compromise \end{itemize} \end{itemize} - + \end{frame} % ================================= @@ -800,7 +800,7 @@ \section{Code Optimization} % \item Integration method: \textbf{sum} (Cython or numpy) on midpoint % \end{itemize} % \end{frame} -% +% % %% ================================= @@ -860,10 +860,10 @@ \section{Code Optimization} \item use cython decorators: \textbf{@cython.boundscheck(False)}, \textbf{@cython.wraparound(False)} \item release the gil (when possible)! \item see more: \href{https://github.com/ToFuProject/tofu/blob/master/Notebooks/Cython_speedup_notes.ipynb}{\textcolor{blue}{my notes}}, \href{https://github.com/jeremiedbb/tutorial-euroscipy-2019}{\textcolor{blue}{Jeremie du Boisberranger's tutorial}} - + \end{itemize} \end{itemize} - + \end{frame} % ================================= @@ -876,9 +876,9 @@ \section{What's next} \frametitle{Tofu's structure} \begin{center} - \includegraphics[width=0.8\linewidth]{figures/tofu.pdf} + \includegraphics[width=0.8\linewidth]{figures/tofu.pdf} \end{center} - + \end{frame} % ================================= @@ -931,7 +931,7 @@ \section{What's next} % \end{tabular} %\end{table} % -% +% %\end{frame} %% ================================= @@ -957,7 +957,7 @@ \section{What's next} \begin{column}{0.4\textwidth} \begin{center} - \hspace{-0.5cm}\includegraphics[width=\linewidth]{figures/cones.png} + \hspace{-0.5cm}\includegraphics[width=\linewidth]{figures/cones.png} \end{center} \end{column} \end{columns} @@ -981,7 +981,7 @@ \section{What's next} %\begin{center} % \includegraphics[width=0.6\linewidth]{figures/meshes.png} -%\end{center} +%\end{center} \end{frame} % ================================= @@ -1020,7 +1020,7 @@ \section{What's next} \begin{center} - \includegraphics[width=0.85\linewidth]{figures/meshes.png} + \includegraphics[width=0.85\linewidth]{figures/meshes.png} \end{center} \vspace{-0.5cm} \end{frame} @@ -1046,7 +1046,7 @@ \subsection*{Splines interpolation} \begin{frame} \frametitle{B(asis)-Splines basis*} - B-Splines of degree $d$ are defined by the \textbf{recursion} formula: + B-Splines of degree $d$ are defined by the \textbf{recursion} formula: \begin{equation} B_j^{d+1}(x)= \dfrac{x - x_j}{x_{j+d}-x_j} B_j^d(x)+ \dfrac{x_{j+1} - x}{x_{j+d+1} - x_{j+1}} B_{j+1}^d (x) @@ -1055,7 +1055,7 @@ \subsection*{Splines interpolation} Some important properties about B-splines: \begin{itemize} - \item Piece-wise polynomials of degree $d \quad \Rightarrow$ \textbf{smoothness} + \item Piece-wise polynomials of degree $d \quad \Rightarrow$ \textbf{smoothness} \item Compact support $\Rightarrow$ \textbf{sparse matrix system} \item Partition of unity $\sum_j Bj (x) = 1$, $\forall x \quad \Rightarrow$ \textbf{conservation laws} \end{itemize} @@ -1071,4 +1071,3 @@ \subsection*{Splines interpolation} \end{document} - diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.out b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.out new file mode 100644 index 000000000..59e3f2f75 --- /dev/null +++ b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.out @@ -0,0 +1,11 @@ +\BOOKMARK [0][-]{chapter.1}{Geometry}{} +\BOOKMARK [1][-]{section.1.1}{Generic cone and plane}{chapter.1} +\BOOKMARK [1][-]{section.1.2}{Intersection}{chapter.1} +\BOOKMARK [1][-]{section.1.3}{Parametric equation}{chapter.1} +\BOOKMARK [2][-]{subsection.1.3.1}{From bragg angle and parameter to local cartesian coordinates}{section.1.3} +\BOOKMARK [2][-]{subsection.1.3.2}{From local cartesian coordinates to bragg angle}{section.1.3} +\BOOKMARK [1][-]{section.1.4}{Generalization}{chapter.1} +\BOOKMARK [2][-]{subsection.1.4.1}{Direct problem}{section.1.4} +\BOOKMARK [0][-]{appendix.A}{Appendices}{} +\BOOKMARK [1][-]{section.A.1}{Section}{appendix.A} +\BOOKMARK [2][-]{subsection.A.1.1}{Subsection}{section.A.1} diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.pdf b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.pdf new file mode 100644 index 000000000..8db47b61d Binary files /dev/null and b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.pdf differ diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.tex b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.tex new file mode 100644 index 000000000..7c61e2d67 --- /dev/null +++ b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.tex @@ -0,0 +1,534 @@ +\documentclass[a4paper,11pt,twoside,titlepage,openright]{book} + +\usepackage[english]{babel} +\usepackage{color} +\usepackage{graphicx} +\usepackage{amsmath} +\numberwithin{equation}{section} +\usepackage[margin=3cm]{geometry} +\usepackage{hyperref} +\usepackage{epsfig,amsfonts} +\usepackage{xcolor,import} + + +\pagestyle{plain} + +\newcommand{\ud}[1]{\underline{#1}} +\newcommand{\e}[1]{\underline{e}_{#1}} +\newcommand{\lt}{\left} +\newcommand{\rt}{\right} +\DeclareMathOperator{\n}{\underline{n}} +\DeclareMathOperator{\ei}{\underline{e}_1} +\DeclareMathOperator{\et}{\underline{e}_2} +\DeclareMathOperator{\ex}{\underline{e}_x} +\DeclareMathOperator{\ey}{\underline{e}_y} +\DeclareMathOperator{\ez}{\underline{e}_z} +\DeclareMathOperator{\nin}{\underline{n}_{in}} +\DeclareMathOperator{\nout}{\underline{n}_{out}} +\DeclareMathOperator{\np}{\underline{n}_{P}} +\DeclareMathOperator{\bragg}{\theta_{bragg}} +\DeclareMathOperator{\DD}{\cos(\theta)^2 - \sin(\psi)^2} +\newcommand{\wdg}{\wedge} +\newcommand{\hypot}[1]{\textbf{\textcolor{green}{#1}}} + + +\begin{document} + +\title{ToFu geometric tools\\ Intersection of a cone with a plane} +\author{Didier VEZINET} +\date{15.10.2019} +\maketitle + +\tableofcontents + +\chapter{Geometry} +\label{chap:Geometry} + +\section{Generic cone and plane} + +Let's consider a half-cone $C_1$ (defined only for $z > 0$), with summit on the cartesian frame's origin $(O, \ex, \ey, \ez)$. +The cone's axis is the $(O,\ez)$ axis. +It's angular opening is $\theta$. + +Let's consider plane $P_1$, of normal $\n$, intersection axis $(O,\ez)$ at point $P$ of coordinates $(0,0,Z_P)$. +Vector $\n$ is oriented by angles $\phi$ and $\psi$ such that one can define the local frame $(P, \ei, \et, \n)$: +$$ +\lt\{ + \begin{array}{ll} + \ei & = \cos(\phi)\ex + \sin(\phi)\ey\\ + \et & = \lt(-\sin(\phi)\ex + \cos(\phi)\ey\rt)\cos(\psi) + \sin(\psi)\ez\\ + \n & = \ei \wdg \et\\ + & = \lt( \sin(\phi)\ex - \cos(\phi)\ey \rt)\sin(\psi) + \cos(\psi)\ez + \end{array} +\rt. +$$ + +We want to find all points $M$ of coordinates $(x, y, z)$ and $(x_1, x_2)$ belonging both to the cone $C_1$ and the plane $P_1$. + +$$ +M \in C_1 \Leftrightarrow \underline{OM}.\ez = \cos(\theta) \|\underline{OM}\| +$$ + +$$ +M \in P_1 \Leftrightarrow \underline{PM}.\n = 0 +$$ + + +\section{Intersection} + +If $M$ belongs to both $P_1$ and $C_1$, then: +$$ +(\underline{OM}.\ez)^2 = \cos(\theta)^2 \|\underline{OM}\|^2 +$$ + +Given that: +$$ +\begin{array}{lll} + \underline{OM} & = \underline{OP} + \underline{PM}\\ + & = Z_P\ez + x_1\ei + x_2\et\\ + & = Z_P\ez + x_1\lt(\cos(\phi)\ex + \sin(\phi)\ey\rt) + x_2\lt(\lt(-\sin(\phi)\ex + \cos(\phi)\ey\rt)\cos(\psi) + \sin(\psi)\ez\rt)\\ + & = Z_P\ez + x_1\cos(\phi)\ex + x_1\sin(\phi)\ey - x_2\sin(\phi)\cos(\psi)\ex + x_2\cos(\phi)\cos(\psi)\ey + x_2\sin(\psi)\ez\\ + & = \lt(x_1\cos(\phi) - x_2\sin(\phi)\cos(\psi)\rt)\ex + \lt(x_1\sin(\phi) + x_2\cos(\phi)\cos(\psi)\rt)\ey + \lt(Z_P + x_2\sin(\psi)\rt)\ez\\ +\end{array} +$$ + +We have: +$$ +\begin{array}{lll} + (\underline{OM}.\ez)^2 & = \lt(Z_P + x_2\sin(\psi)\rt)^2\\ + & = Z_P^2 + 2Z_Px_2\sin(\psi) + x_2^2\sin(\psi)^2 +\end{array} +$$ + +And: +$$ +\begin{array}{lll} + \|\underline{OM}\|^2 & = & \| \lt(x_1\cos(\phi) - x_2\sin(\phi)\cos(\psi)\rt)\ex + \lt(x_1\sin(\phi) + x_2\cos(\phi)\cos(\psi)\rt)\ey + \lt(Z_P + x_2\sin(\psi)\rt)\ez\|^2\\ + & = & \lt( x_1\cos(\phi) - x_2\sin(\phi)\cos(\psi) \rt)^2\\ + & & + \lt(x_1\sin(\phi) + x_2\cos(\phi)\cos(\psi)\rt)^2\\ + & & + \lt(Z_P + x_2\sin(\psi)\rt)^2\\ + & = & x_1^2\cos(\phi)^2 - 2x_1x_2\cos(\phi)\sin(\phi)\cos(\psi) + x_2^2\sin(\phi)^2\cos(\psi)^2\\ + & & + x_1^2\sin(\phi)^2 + 2x_1x_2\sin(\phi)\cos(\phi)\cos(\psi) + x_2^2\cos(\phi)^2\cos(\psi)^2\\ + & & + Z_P^2 + 2Z_Px_2\sin(\psi) + x_2^2\sin(\psi)^2\\ + & = & x_1^2 + x_2^2\cos(\psi)^2\\ + & & + Z_P^2 + 2Z_Px_2\sin(\psi) + x_2^2\sin(\psi)^2\\ + & = & x_1^2 + x_2^2 + 2Z_Px_2\sin(\psi) + Z_P^2\\ +\end{array} +$$ + +Thus: +$$ +\begin{array}{lll} + & (\underline{OM}.\ez)^2 = \cos(\theta)^2 \|\underline{OM}\|^2\\ + \Leftrightarrow & Z_P^2 + 2Z_Px_2\sin(\psi) + x_2^2\sin(\psi)^2 = \cos(\theta)^2 \lt(x_1^2 + x_2^2 + 2Z_Px_2\sin(\psi) + Z_P^2\rt)\\ + \Leftrightarrow & Z_P^2\lt(1-\cos(\theta)^2\rt) + 2Z_Px_2\sin(\psi)\lt(1-\cos(\theta)^2\rt) = x_1^2\cos(\theta)^2 + x_2^2\lt(\cos(\theta)^2 - \sin(\psi)^2\rt)\\ + \Leftrightarrow & Z_P^2\sin(\theta)^2 + 2Z_Px_2\sin(\psi)\sin(\theta)^2 = x_1^2\cos(\theta)^2 + x_2^2\lt(\cos(\theta)^2 - \sin(\psi)^2\rt)\\ +\end{array} +$$ + +Considering that by hypothesis $\theta>0$: +$$ +\begin{array}{lll} + & (\underline{OM}.\ez)^2 = \cos(\theta)^2 \|\underline{OM}\|^2\\ + \Leftrightarrow & x_1^2\cos(\theta)^2 + x_2^2\lt(\DD\rt) - 2Z_Px_2\sin(\psi)\sin(\theta)^2 - Z_P^2\sin(\theta)^2 = 0\\ + \Leftrightarrow & x_1^2\frac{\cos(\theta)^2}{\DD} + x_2^2 - 2x_2Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD} - Z_P^2\frac{\sin(\theta)^2}{\DD} = 0\\ + \Leftrightarrow & x_1^2\frac{\cos(\theta)^2}{\DD} + \lt(x_2 - Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD}\rt)^2 - Z_P^2\frac{\sin(\psi)^2\sin(\theta)^4}{\lt(\DD\rt)^2} - Z_P^2\frac{\sin(\theta)^2}{\DD} = 0\\ + \Leftrightarrow & x_1^2\frac{\cos(\theta)^2}{\DD} + \lt(x_2 - Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD}\rt)^2 = Z_P^2 \frac{\sin(\theta)^2}{\lt(\DD\rt)^2} \lt(\sin(\psi)^2\sin(\theta)^2 + \DD \rt)\\ + \Leftrightarrow & x_1^2\frac{\cos(\theta)^2}{\DD} + \lt(x_2 - Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD}\rt)^2 = Z_P^2 \frac{\sin(\theta)^2}{\lt(\DD\rt)^2} \lt(-\sin(\psi)^2\cos(\theta)^2 + \cos(\theta)^2 \rt)\\ + \Leftrightarrow & x_1^2\frac{\cos(\theta)^2}{\DD} + \lt(x_2 - Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD}\rt)^2 = Z_P^2 \frac{\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2}{\lt(\DD\rt)^2} \\ + \Leftrightarrow & \frac{x_1^2}{ Z_P^2\frac{\sin(\theta)^2\cos(\psi)^2}{\DD} } + \frac{\lt(x_2-Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD}\rt)^2}{ Z_P^2\frac{\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2}{\lt(\DD\rt)^2} } = 1 +\end{array} +$$ + +Or, in a reduced conic form: + +$$ +\frac{x_1^2}{ a^2 } + \frac{\lt(x_2-x_2(C)\rt)^2}{b^2} = 1 +$$ + + +With: +$$ +\lt\{ + \begin{array}{lll} + x_2(C) & = Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD} & \text{ $x_2$ coordinate of the center} \\ + a^2 & = Z_P^2\frac{\sin(\theta)^2\cos(\psi)^2}{\DD} & \text{ squared minor radius} \\ + b^2 & = Z_P^2\frac{\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2}{\lt(\DD\rt)^2} & \text{ squared major radius} \\ + b^2 & = a^2\frac{\cos(\theta)^2}{\DD} \Leftrightarrow a^2 = b^2\lt(1-\frac{\sin(\psi)^2}{\cos(\theta)^2}\rt) + \end{array} +\rt. +$$ + +The distance $d_{CF}$ between the center $C$ and the focal point $F$ can be deduced from: +$$ +\begin{array}{lll} + d_{CF}^2 & = b^2-a^2\\ + & = b^2\frac{\sin(\psi)^2}{\cos(\theta)^2}\\ + & = Z_P^2\frac{\sin(\theta)^2\cos(\psi)^2\sin(\psi)^2}{\lt(\DD\rt)^2} +\end{array} +$$ + +Hence, the $x_2$ coordinate of $F$ is: +$$ +\begin{array}{lll} + x_2(F) & = x_2(C) \pm d_{CF}\\ + & = Z_P\frac{\sin(\psi)\sin(\theta)^2}{\DD} \pm Z_P\frac{\sin(\theta)\cos(\psi)\sin(\psi)}{\DD}\\ + & = Z_P\frac{\sin(\psi)\sin(\theta)^2 \pm \sin(\theta)\cos(\psi)\sin(\psi)}{\DD}\\ + & = Z_P\frac{\sin(\psi)\sin(\theta)}{\DD}\lt(\sin(\theta) \pm \cos(\psi)\rt)\\ +\end{array} +$$ + +It is worth noticing that the neither the focal point nor the center correspond to the intersection between the axes and the plane $P$. + + +\section{Parametric equation} + +In our case, only the axes $(O, \ez)$, fixed by the crystal's summit and normal, is independent from the cone's angular opening $\theta$. +It makes sense to design an ad-hoc coordinate system centered on the ellipse's center $C$ to use its parameterized equation. + +Knowing all geometrical parameters, it is possible to compute all points on the ellipse parameterizing them with $t$: +$$ +\lt\{ + \begin{array}{lll} + x_1 = a\cos(t)\\ + x_2 = x_2(C) + b\sin(t) + \end{array} +\rt. +$$ + +\emph{BEWARE} : parameter $t$ is not the angle of the point with respect to the ellipse's center. + +Now, we would like to parameterize the ellipse not with $t$ but with the angle $\beta$ with respect to the point $P$, because it is the physically relevant angle since it is taken with respect to the axis $(O, \ez)$ and relates to the impact point of the photon beam on the crystal's center. Also, it is the only common element to all ellipses. +The angle $\epsilon$ taken with respect to the center is not relevant because each ellipse has a different center. + +In this perspective: +$$ +\lt\{ + \begin{array}{lll} + x_1 = l(\beta)\cos(\beta)\\ + x_2 = l(\beta)\sin(\beta) + \end{array} +\rt. +$$ + +Keeping in mind that the ellipse is defined as: +$$ +\frac{x_1^2}{a^2} + \frac{(x_2-x_2(C))^2}{b^2} = 1 +$$ + +We can write: +$$ +\begin{array}{lll} + & l^2b^2\cos(\beta)^2 + a^2\lt(l^2\sin(\beta)^2 - 2lx_2(C)\sin(\beta) + x_2(C)^2\rt) = a^2b^2\\ + \Leftrightarrow + & l^2\lt(b^2\cos(\beta)^2 + a^2\sin(\beta)^2\rt) - 2la^2x_2(C)\sin(\beta) + a^2x_2(C)^2 - a^2b^2 = 0 +\end{array} +$$ + +Has solutions if: +$$ +\begin{array}{lll} + & \Delta = 4a^4x_2(C)^2\sin(\beta)^2 - 4\lt(b^2\cos(\beta)^2 + a^2\sin(\beta)^2\rt)\lt(a^2x_2(C)^2 - a^2b^2\rt) \geq 0\\ + \Leftrightarrow + & \Delta = 4a^2\lt[a^2x_2(C)^2\sin(\beta)^2 - \lt(b^2\cos(\beta)^2 + a^2\sin(\beta)^2\rt)\lt(x_2(C)^2 - b^2\rt)\rt] \geq 0\\ + \Leftrightarrow + & \Delta = 4a^2\lt[a^2x_2(C)^2\sin(\beta)^2 - b^2x_2(C)^2\cos(\beta)^2 - a^2x_2(C)^2\sin(\beta)^2 + b^4\cos(\beta)^2 + a^2b^2\sin(\beta)^2\rt] \geq 0\\ + \Leftrightarrow + & \Delta = 4a^2\lt[-b^2x_2(C)^2\cos(\beta)^2 + b^4\cos(\beta)^2 + a^2b^2\sin(\beta)^2\rt] \geq 0\\ + \Leftrightarrow + & \Delta = 4a^2b^2\lt[-x_2(C)^2\cos(\beta)^2 + b^2\cos(\beta)^2 + a^2\sin(\beta)^2\rt] \geq 0\\ +\end{array} +$$ + +Which is equivalent to, keeping in mind that $b^2-a^2 = d_{CF}^2$: +$$ +\begin{array}{lll} + & \Delta = 4a^2b^2\lt[a^2 + \lt(b^2-a^2-x_2(C)^2\rt)\cos(\beta)^2\rt] \geq 0\\ + \Leftrightarrow + & \lt(b^2-a^2-x_2(C)^2\rt)\cos(\beta)^2 \geq -a^2 + \Leftrightarrow + & \lt(d_{CF}^2-x_2(C)^2\rt)\cos(\beta)^2 \geq -a^2 +\end{array} +$$ + +If $d_{CF}^2-x_2(C)^2 \ geq 0$, this is true for all $\beta$ values, and this condition is met if: +$$ +\begin{array}{lll} + & b^2-a^2-x_2(C)^2 \geq 0\\ + \Leftrightarrow + & \frac{Z_P^2}{\lt(\DD\rt)^2}\lt(\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2 - \sin(\theta)^2\cos(\psi)^2\cos(\theta)^2 + \sin(\theta)^2\cos(\psi)^2\sin(\psi)^2 - \sin(\psi)^2\sin(\theta)^4\rt) \geq 0\\ + \Leftrightarrow + & \sin(\theta)^2\sin(\psi)^2\lt(\cos(\psi)^2 -\sin(\theta)^2\rt) \geq 0 +\end{array} +$$ + +Which is true if we have an ellipse, which is the only case of interest. +Hence $\Delta=4a^2b^2\lt[a^2 + \lt(b^2-a^2-x_2(C)^2\rt)\cos(\beta)^2\rt] \geq 0$, so: +$$ +\begin{array}{lll} + & l_{1,2} = \frac{2a^2x_2(C)\sin(\beta) \pm \sqrt{\Delta}}{2\lt(b^2\cos(\beta)^2 + a^2\sin(\beta)^2\rt)}\\ + \Leftrightarrow + & l_{1,2} = \frac{a^2x_2(C)\sin(\beta) \pm ab\sqrt{a^2 + \lt(b^2-a^2-x_2(C)^2\rt)\cos(\beta)^2}}{b^2\cos(\beta)^2 + a^2\sin(\beta)^2}\\ +\end{array} +$$ + +And we only want the positive solution: +$$ +l = \frac{a^2x_2(C)\sin(\beta) + ab\sqrt{a^2 + \lt(b^2-a^2-x_2(C)^2\rt)\cos(\beta)^2}}{b^2\cos(\beta)^2 + a^2\sin(\beta)^2} +$$ + + + + +\subsection{From bragg angle and parameter to local cartesian coordinates} + +Keep in mind that the frame $(P, \ei, \et)$ is, by definition aligned on the minor and major axes of the ellipse. +Hence, for an arbitrary frame $(R, \ud{e}_i, \ud{e}_j)$ on plane $P_1$, translated and rotated by $\alpha$ with respect to $(P, \ei, \et)$: + +$$ +\lt\{ + \begin{array}{lll} + \ud{e}_i = \cos(\alpha)\ei + \sin(\alpha)\et\\ + \ud{e}_j = -\sin(\alpha)\ei + \cos(\alpha)\et\\ + \ei = \cos(\alpha)\ud{e_i} - \sin(\alpha)\ud{e}_j + \et = \sin(\alpha)\ud{e_i} + \cos(\alpha)\ud{e}_j + \end{array} +\rt. +$$ + +Or, in coordinate tranforms: +$$ +\lt\{ + \begin{array}{lll} + x_1 = x_1(R) + x_i\cos(\alpha) - x_j\sin(\alpha)\\ + x_2 = x_2(R) + x_i\sin(\alpha) + x_j\cos(\alpha)\\ + x_i = (x_1-x_1(R))\cos(\alpha) + (x_2-x_2(R))\sin(\alpha)\\ + x_j = -(x_1-x_1(R))\sin(\alpha) + (x_2-x_2(R))\cos(\alpha) + \end{array} +\rt. +$$ + +Hence: +$$ +\lt\{ + \begin{array}{lll} + x_i = \lt(a\cos(\epsilon)-x_1(R)\rt)\cos(\alpha) + \lt(x_2(C)-x_2(R) + b\sin(\epsilon)\rt)\sin(\alpha)\\ + x_j = -\lt(a\cos(\epsilon)-x_1(R)\rt)\sin(\alpha) + \lt(x_2(C)-x_2(R) + b\sin(\epsilon)\rt)\cos(\alpha) + \end{array} +\rt. +$$ + +But +$$ +\lt\{ + \begin{array}{lll} + \|\ud{PM}\|^2 = x_1^2 + x_2^2 = \\ + x_1 = \|\ud{PM}\|\cos(\beta)\\ + x_2 = \|\ud{PM}\|\sin(\beta) + \end{array} +\rt. +$$ + + + +\subsection{From local cartesian coordinates to bragg angle} + +Knowing $(x_i, x_j)$ and all geometric parameters, we now want to derive $(\theta, \epsilon)$. + +From the previous equation, we can write: +$$ +\lt\{ + \begin{array}{lll} + x_i\cos(\alpha) - x_j\sin(\alpha) = a\cos(\epsilon)-x_1(R) & (1)\\ + x_i\sin(\alpha) + x_j\cos(\alpha) = x_2(C)-x_2(R) + b\sin(\epsilon) & + (2) + \end{array} +\rt. +$$ + +The dependency in $\theta$ is hidden in the expressions of $a$, $b$ and +$x_2(C)$. + +By squaring and summing, it is possible to get rid of the $\epsilon$ +dependency: +$$ +\lt\{ + \begin{array}{lll} + a^2\cos(\epsilon)^2 = \lt(x_i\cos(\alpha) - x_j\sin(\alpha) + + x_1(R)\rt)^2\\ + b^2\sin(\epsilon)^2 = \lt(x_i\sin(\alpha) + x_j\cos(\alpha) - x_2(C) + + x_2(R)\rt)^2 + \end{array} +\rt. +$$ + +Hence, keeping in mind that $a^2 = b^2\frac{\DD}{\cos(\theta)^2}$ and re-using the definitions of $x_1$ and $x_2$ which do not depend on the unknowns $(\theta, \epsilon)$: +$$ +\begin{array}{lllll} + & b^2 & = & \frac{\cos(\theta)^2}{\DD}x_1^2 + \lt(x_2 - x_2(C)\rt)^2\\ + \Leftrightarrow + & b^2\lt(\DD\rt) & = & \cos(\theta)^2x_1^2 + \lt(\DD\rt)\lt(x_2^2 - 2x_2x_2(C) + x_2(C)^2\rt)\\ + \Leftrightarrow + & Z_P^2\frac{\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2}{\DD} & = & \cos(\theta)^2x_1^2 + \lt(\DD\rt)\lt(x_2^2 - 2x_2x_2(C) + x_2(C)^2\rt)\\ + \Leftrightarrow + & Z_P^2\frac{\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2}{\DD} & = & \cos(\theta)^2\lt(x_1^2+x_2^2\rt) - \sin(\psi)^2x_2^2 - \lt(\DD\rt)\lt(2x_2x_2(C) - x_2(C)^2\rt)\\ + \Leftrightarrow + & Z_P^2\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2 & = & \cos(\theta)^2\lt(\DD\rt)\lt(x_1^2+x_2^2\rt) - \sin(\psi)^2\lt(\DD\rt)x_2^2\\ + \Leftrightarrow + & & & - \lt(\DD\rt)^2\lt(2x_2x_2(C) - x_2(C)^2\rt)\\ + \Leftrightarrow + & Z_P^2\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2 & = &\cos(\theta)^4\lt(x_1^2+x_2^2\rt) - \cos(\theta)^2\sin(\psi)^2\lt(x_1^2+2x_2^2\rt) + \sin(\psi)^4x_2^2\\ + & & & - \lt(\DD\rt)^2\lt(2x_2x_2(C) - x_2(C)^2\rt)\\ + \Leftrightarrow + & Z_P^2\sin(\theta)^2\cos(\psi)^2\cos(\theta)^2 & = & \cos(\theta)^4\lt(x_1^2+x_2^2\rt) - \cos(\theta)^2\sin(\psi)^2\lt(x_1^2+2x_2^2\rt) + \sin(\psi)^4x_2^2\\ + & & & - \lt[2x_2Z_P\sin(\psi)\sin(\theta)^2\lt(\DD\rt) - Z_P^2\sin(\psi)^2\sin(\theta)^4\rt]\\ + \Leftrightarrow + & Z_P^2\cos(\psi)^2\lt(\cos(\theta)^2-\cos(\theta)^4\rt) & = &\cos(\theta)^4\lt(x_1^2+x_2^2\rt) - \cos(\theta)^2\sin(\psi)^2\lt(x_1^2+2x_2^2\rt) + \sin(\psi)^4x_2^2\\ + & & & - 2x_2Z_P\sin(\psi)\sin(\theta)^2\lt(\DD\rt)\\ + & & & + Z_P^2\sin(\psi)^2\lt(1-2\cos(\theta)^2+\cos(\theta)^4\rt)\\ + \Leftrightarrow + & Z_P^2\cos(\psi)^2\lt(\cos(\theta)^2-\cos(\theta)^4\rt) & = &\cos(\theta)^4\lt(x_1^2+x_2^2\rt) - \cos(\theta)^2\sin(\psi)^2\lt(x_1^2+2x_2^2\rt) + \sin(\psi)^4x_2^2\\ + & & & - 2x_2Z_P\sin(\psi)\lt(\cos(\theta)^2 - \cos(\theta)^4 - \sin(\psi)^2\lt(1-\cos(\theta)^2\rt)\rt)\\ + & & & + Z_P^2\sin(\psi)^2\lt(1-2\cos(\theta)^2+\cos(\theta)^4\rt) +\end{array} +$$ + +Then introducing $X = \cos(\theta)^2$: +$$ +\begin{array}{lllll} + & b^2 & = & \frac{\cos(\theta)^2}{\DD}x_1^2 + \lt(x_2 - x_2(C)\rt)^2\\ + \Leftrightarrow + & Z_P^2\cos(\psi)^2\lt(X-X^2\rt) & = & X^2\lt(x_1^2+x_2^2\rt) - X\sin(\psi)^2\lt(x_1^2+2x_2^2\rt) + \sin(\psi)^4x_2^2\\ + & & & - 2x_2Z_P\sin(\psi)\lt(X - X^2 - \sin(\psi)^2 + X\sin(\psi)^2\rt)\\ + & & & + Z_P^2\sin(\psi)^2\lt(1 - 2X + X^2\rt) +\end{array} +$$ + +Which boils down to: +$$ +\begin{array}{lllll} + & 0 & = & X^2\lt[ x_1^2+x_2^2 + 2x_2Z_P\sin(\psi) + Z_P^2\sin(\psi)^2 + Z_P^2\cos(\psi)^2 \rt]\\ + & & & + X\lt[ -\sin(\psi)^2\lt(x_1^2+2x_2^2\rt) - \lt(1+\sin(\psi)^2\rt)2x_2Z_P\sin(\psi) - 2Z_P^2\sin(\psi)^2 - Z_P^2\cos(\psi)^2\rt]\\ + & & & + \sin(\psi)^4x_2^2 + 2x_2Z_P\sin(\psi)^3 + Z_P^2\sin(\psi)^2\\ + \Leftrightarrow + & 0 & = & X^2\lt[ x_1^2 + \lt(x_2 + Z_P\sin(\psi)\rt)^2 + Z_P^2\cos(\psi)^2 \rt]\\ + & & & + X\lt[ -\sin(\psi)^2\lt(x_1^2+2x_2^2\rt) - \lt(1+\sin(\psi)^2\rt)2x_2Z_P\sin(\psi) - 2Z_P^2\sin(\psi)^2 - Z_P^2\cos(\psi)^2\rt]\\ + & & & + \sin(\psi)^2\lt(\sin(\psi)^2x_2^2 + 2x_2Z_P\sin(\psi) + Z_P^2\rt)\\ + \Leftrightarrow + & 0 & = & X^2\lt[ x_1^2 + \lt(x_2 + Z_P\sin(\psi)\rt)^2 + Z_P^2\cos(\psi)^2 \rt]\\ + & & & - X\lt[ \sin(\psi)^2\lt(x_1^2+2x_2^2\rt) + \lt(1+\sin(\psi)^2\rt)2x_2Z_P\sin(\psi) + 2Z_P^2\sin(\psi)^2 + Z_P^2\cos(\psi)^2\rt]\\ + & & & + \sin(\psi)^2\lt(x_2\sin(\psi) + Z_P\rt)\\ + \Leftrightarrow + & 0 & = & X^2\lt[ x_1^2 + \lt(x_2 + Z_P\sin(\psi)\rt)^2 + Z_P^2\cos(\psi)^2 \rt]\\ + & & & - X\lt[ \sin(\psi)^2\lt(x_1^2+x_2^2\rt) + x_2^2\sin(\psi)^2 + 2x_2Z_P\sin(\psi) + 2x_2Z_P\sin(\psi)^3 + Z_P^2 + Z_P^2\sin(\psi)^2\rt]\\ + & & & + \sin(\psi)^2\lt(x_2\sin(\psi) + Z_P\rt)\\ + \Leftrightarrow + & 0 & = & X^2\lt[ x_1^2 + \lt(x_2 + Z_P\sin(\psi)\rt)^2 + Z_P^2 - Z_P^2\sin(\psi)^2 \rt]\\ + & & & - X\lt[ \sin(\psi)^2\lt(x_1^2+x_2^2\rt) + \lt(x_2\sin(\psi) + Z_P\rt)^2 + Z_P\sin(\psi)^2\lt(2x_2\sin(\psi) + Z_P\rt)\rt]\\ + & & & + \sin(\psi)^2\lt(x_2\sin(\psi) + Z_P\rt)\\ + \Leftrightarrow + & 0 & = & X^2\lt[ x_1^2 + x_2^2 + \lt(x_2\sin(\psi) + Z_P\rt)^2 - x_2^2\sin(\psi)^2 \rt]\\ + & & & - X\lt[ \sin(\psi)^2\lt(x_1^2+x_2^2\rt) + \lt(x_2\sin(\psi) + Z_P\rt)^2 + Z_P\sin(\psi)^2\lt(2x_2\sin(\psi) + Z_P\rt)\rt]\\ + & & & + \sin(\psi)^2\lt(x_2\sin(\psi) + Z_P\rt)\\ + \Leftrightarrow + & 0 & = & AX^2 - BX + C +\end{array} +$$ + +With: + +$$ +\lt\{ + \begin{array}{lll} + A = x_1^2 + x_2^2 + \lt(x_2\sin(\psi) + Z_P\rt)^2 - x_2^2\sin(\psi)^2\\ + B = \sin(\psi)^2\lt(x_1^2+x_2^2\rt) + \lt(x_2\sin(\psi) + Z_P\rt)^2 + Z_P\sin(\psi)^2\lt(2x_2\sin(\psi) + Z_P\rt)\\ + C = \sin(\psi)^2\lt(x_2\sin(\psi) + Z_P\rt) + \end{array} +\rt. +$$ + +Solutions exist if: +$$ +\begin{array}{lllll} + & \Delta & \geq & 0\\ + \Leftrightarrow + & B^2 - 4AC & \geq & 0\\ +\end{array} +$$ + +In which case, only solutions in $\lt[0,1\rt]$ are acceptable: +$$ +\cos(\theta)^2 = \frac{B \pm \sqrt(\Delta)}{2A} \in \lt[0,1\rt] +$$ + +And by definition, $\theta \in \lt[0,\frac{\pi}{2}\rt]$, hence $\cos(\theta)\geq0$ and: +$$ +\theta = \arccos\lt(\sqrt{\frac{B \pm \sqrt(\Delta)}{2A}}\rt) +$$ + +\subsubsection{Alternative method for $\theta$} + +By definition: +$$ +\underline{OM}.\ez = \cos(\theta)\|OM\| +$$ + +And: +$$ +\begin{array}{lllll} + \underline{OM} & = & \underline{OP} + \underline{PR} + \underline{RM}\\ + & = & \\ +\end{array} +$$ + +\newpage +\section{Generalization} + +This time, the crystal of curvature radius $R$ has center $C$ of coordinates $(x_C, y_C, z_C)$ in the tokamak's frame $(O, \ex, \ey, \ez)$.\\ +The direct orthonormal systems are: + +\begin{figure}[hbt] + \includegraphics[width=\linewidth]{SpectroX2D_Crystal.pdf} + \caption{Definition of the generilazed geometry} + \label{fig:sketch} +\end{figure} + + +$$ +\lt\{ + \begin{array}{lll} + (O, \ex, \ey, \ez)\\ + (C, \nout, \e{1}, \e{2})\\ + (P, \np, \e{a}, \e{b})\\ + (R, \np, \e{i}, \e{j}) + \end{array} +\rt. +$$ + +\subsection{Direct problem} + +We know all geometrical parameters, in particular, we know: +$$ +\lt\{ + \begin{array}{lll} + \ud{OC} = x(C)\ex + y(C)\ey + z(C)\ez\\ + \ud{CS} = R\nout\\ + \ud{SP} = -L\nout\\ + \ud{PR} = x_a(R)\e{a} + x_b(R)\e{b}\\ + \ud{RM} = x_i\e{i} + x_j\e{j} + \end{array} +\rt. +$$ + + + + + + +\appendix +\chapter{Appendices} + +\section{Section} +\subsection{Subsection} + +\end{document} diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.toc b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.toc new file mode 100644 index 000000000..bc3fa5b3d --- /dev/null +++ b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConePlaneEllipse.toc @@ -0,0 +1,13 @@ +\select@language {english} +\contentsline {chapter}{\numberline {1}Geometry}{5}{chapter.1} +\contentsline {section}{\numberline {1.1}Generic cone and plane}{5}{section.1.1} +\contentsline {section}{\numberline {1.2}Intersection}{5}{section.1.2} +\contentsline {section}{\numberline {1.3}Parametric equation}{7}{section.1.3} +\contentsline {subsection}{\numberline {1.3.1}From bragg angle and parameter to local cartesian coordinates}{8}{subsection.1.3.1} +\contentsline {subsection}{\numberline {1.3.2}From local cartesian coordinates to bragg angle}{8}{subsection.1.3.2} +\contentsline {subsubsection}{Alternative method for $\theta $}{10}{section*.2} +\contentsline {section}{\numberline {1.4}Generalization}{11}{section.1.4} +\contentsline {subsection}{\numberline {1.4.1}Direct problem}{11}{subsection.1.4.1} +\contentsline {chapter}{\numberline {A}Appendices}{13}{appendix.A} +\contentsline {section}{\numberline {A.1}Section}{13}{section.A.1} +\contentsline {subsection}{\numberline {A.1.1}Subsection}{13}{subsection.A.1.1} diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.pdf b/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.pdf new file mode 100644 index 000000000..6e6da4d54 Binary files /dev/null and b/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.pdf differ diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.png b/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.png new file mode 100644 index 000000000..0cf50e450 Binary files /dev/null and b/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.png differ diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.svg b/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.svg new file mode 100644 index 000000000..914e52b03 --- /dev/null +++ b/Notes_Upgrades/SpectroX2D/SpectroX2D_Crystal.svg @@ -0,0 +1,1860 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + nout + nin + e1 + e2 + C + P + nP + ei + ej + + + M (xi,xj) + ex + ey + ez + O + v + + nout + psi + + R + + S + + + + + + + + + + e1 + e2 + phi + S + + + + + + + + + nin + v + theta_bragg + S + + + + eb + ea + + + + + + + + ea + ksi + + P + ei + eb + + ej + + + + nout + + P + + psi + nP + + + + + R + L + + + + nP + + ez + + + + + + + + e2 + C + + nout + + + + extent[1] + + S + R + + + + + e1 + + + + diff --git a/XICS_mask.npz b/XICS_mask.npz new file mode 100644 index 000000000..6f349523f Binary files /dev/null and b/XICS_mask.npz differ diff --git a/doc/source/contributing.rst b/doc/source/contributing.rst index d49cafd61..5fbfff79a 100644 --- a/doc/source/contributing.rst +++ b/doc/source/contributing.rst @@ -114,7 +114,6 @@ page `__, from your branch, accept it. - Continuous integration """""""""""""""""""""" diff --git a/doc/source/tofu.data.rst b/doc/source/tofu.data.rst index 423d59c38..420dfcc32 100644 --- a/doc/source/tofu.data.rst +++ b/doc/source/tofu.data.rst @@ -56,4 +56,3 @@ tofu.data.\_plot module :members: :undoc-members: :show-inheritance: - diff --git a/doc/source/tofu.dumpro.rst b/doc/source/tofu.dumpro.rst index 18f2bfe1f..dc987c474 100644 --- a/doc/source/tofu.dumpro.rst +++ b/doc/source/tofu.dumpro.rst @@ -56,4 +56,3 @@ tofu.dumpro.\_plot module :members: :undoc-members: :show-inheritance: - diff --git a/doc/source/tofu.dust.rst b/doc/source/tofu.dust.rst index 5fdc6c022..c1349f951 100644 --- a/doc/source/tofu.dust.rst +++ b/doc/source/tofu.dust.rst @@ -32,4 +32,3 @@ tofu.dust.\_plot module :members: :undoc-members: :show-inheritance: - diff --git a/doc/source/tofu.geom.rst b/doc/source/tofu.geom.rst index e640d87a7..b8a94d4a8 100644 --- a/doc/source/tofu.geom.rst +++ b/doc/source/tofu.geom.rst @@ -96,4 +96,3 @@ tofu.geom.utils module :members: :undoc-members: :show-inheritance: - diff --git a/doc/source/tofu.rst b/doc/source/tofu.rst index cf238008d..570f29bc7 100644 --- a/doc/source/tofu.rst +++ b/doc/source/tofu.rst @@ -58,4 +58,3 @@ tofu.version module :members: :undoc-members: :show-inheritance: - diff --git a/tofu/data/_core.py b/tofu/data/_core.py index 8570e2b05..adbfd68e7 100644 --- a/tofu/data/_core.py +++ b/tofu/data/_core.py @@ -25,11 +25,13 @@ import tofu.data._plot as _plot import tofu.data._def as _def import tofu.data._physics as _physics + import tofu.data._spectrafit2d as _spectrafit2d except Exception: from . import _comp as _comp from . import _plot as _plot from . import _def as _def from . import _physics as _physics + from . import _spectrafit2d as _spectrafit2d __all__ = ['DataCam1D','DataCam2D', 'DataCam1DSpectral','DataCam2DSpectral', diff --git a/tofu/data/_spectrafit2d.py b/tofu/data/_spectrafit2d.py new file mode 100644 index 000000000..7979fb7d7 --- /dev/null +++ b/tofu/data/_spectrafit2d.py @@ -0,0 +1,1484 @@ + +# Built-in +import os +import warnings + +# Common +import numpy as np +import scipy.optimize as scpopt +import scipy.constants as scpct +import scipy.sparse as sparse +from scipy.interpolate import BSpline +import matplotlib.pyplot as plt + + +# -------------- +# TO BE MOVED TO tofu.data WHEN FINISHED !!!! +# -------------- + + +_NPEAKMAX = 12 + +########################################################### +########################################################### +# +# Preliminary +# utility tools for 1d spectral fitting +# +########################################################### +########################################################### + + +def remove_bck(x, y): + # opt = np.polyfit(x, y, deg=0) + opt = [np.nanmin(y)] + return y-opt[0], opt[0] + + +def get_peaks(x, y, nmax=None): + + if nmax is None: + nmax = _NPEAKMAX + + # Prepare + ybis = np.copy(y) + A = np.empty((nmax,), dtype=y.dtype) + x0 = np.empty((nmax,), dtype=x.dtype) + sigma = np.empty((nmax,), dtype=y.dtype) + def gauss(xx, A, x0, sigma): return A*np.exp(-(xx-x0)**2/sigma**2) + def gauss_jac(xx, A, x0, sigma): + jac = np.empty((xx.size, 3), dtype=float) + jac[:, 0] = np.exp(-(xx-x0)**2/sigma**2) + jac[:, 1] = A*2*(xx-x0)/sigma**2 * np.exp(-(xx-x0)**2/sigma**2) + jac[:, 2] = A*2*(xx-x0)**2/sigma**3 * np.exp(-(xx-x0)**2/sigma**2) + return jac + + dx = np.nanmin(np.diff(x)) + + # Loop + nn = 0 + while nn < nmax: + ind = np.nanargmax(ybis) + x00 = x[ind] + if np.any(np.diff(ybis[ind:], n=2) >= 0.): + wp = min(x.size-1, + ind + np.nonzero(np.diff(ybis[ind:],n=2)>=0.)[0][0] + 1) + else: + wp = ybis.size-1 + if np.any(np.diff(ybis[:ind+1], n=2) >= 0.): + wn = max(0, np.nonzero(np.diff(ybis[:ind+1],n=2)>=0.)[0][-1] - 1) + else: + wn = 0 + width = x[wp]-x[wn] + assert width>0. + indl = np.arange(wn, wp+1) + sig = np.ones((indl.size,)) + if (np.abs(np.mean(np.diff(ybis[ind:wp+1]))) + > np.abs(np.mean(np.diff(ybis[wn:ind+1])))): + sig[indl < ind] = 1.5 + sig[indl > ind] = 0.5 + else: + sig[indl < ind] = 0.5 + sig[indl > ind] = 1.5 + p0 = (ybis[ind], x00, width)#,0.) + bounds = (np.r_[0., x[wn], dx/2.], + np.r_[5.*ybis[ind], x[wp], 5.*width]) + try: + (Ai, x0i, sigi) = scpopt.curve_fit(gauss, x[indl], ybis[indl], + p0=p0, bounds=bounds, jac=gauss_jac, + sigma=sig, x_scale='jac')[0] + except Exception as err: + print(str(err)) + import ipdb + ipdb.set_trace() + pass + + ybis = ybis - gauss(x, Ai, x0i, sigi) + A[nn] = Ai + x0[nn] = x0i + sigma[nn] = sigi + + + nn += 1 + return A, x0, sigma + +def get_p0bounds_all(x, y, nmax=None, lamb0=None): + + yflat, bck = remove_bck(x, y) + amp, x0, sigma = get_peaks(x, yflat, nmax=nmax) + lamb0 = x0 + nmax = lamb0.size + + p0 = amp.tolist() + [0 for ii in range(nmax)] + sigma.tolist() + [bck] + + lx = [np.nanmin(x), np.nanmax(x)] + Dx = np.diff(lx) + dx = np.nanmin(np.diff(x)) + + bamp = (np.zeros(nmax,), np.full((nmax,),3.*np.nanmax(y))) + bdlamb = (np.full((nmax,), -Dx/2.), np.full((nmax,), Dx/2.)) + bsigma = (np.full((nmax,), dx/2.), np.full((nmax,), Dx/2.)) + bbck0 = (0., np.nanmax(y)) + + bounds = (np.r_[bamp[0], bdlamb[0], bsigma[0], bbck0[0]], + np.r_[bamp[1], bdlamb[1], bsigma[1], bbck0[1]]) + if not np.all(bounds[0] spectrum {0} / {1}".format(ii+1, nspect)) + + try: + popt, pcov = scpopt.curve_fit(func_sca, x, spectra[ii,:], + jac=func_sca_jac, + p0=p0, bounds=bounds, + max_nfev=max_nfev, xtol=xtol, + x_scale='jac', + verbose=verbose) + except Exception as err: + msg = " Convergence issue for {0} / {1}\n".format(ii+1, nspect) + msg += " => %s\n"%str(err) + msg += " => Resetting initial guess and bounds..." + print(msg) + try: + p0, bounds, _ = get_p0bounds(x, spectra[ii,:], + nmax=nmax, lamb0=lamb0) + popt, pcov = scpopt.curve_fit(func_sca, x, spectra[ii,:], + jac=func_sca_jac, + p0=p0, bounds=bounds, + max_nfev=max_nfev, xtol=xtol, + x_scale='jac', + verbose=verbose) + p0 = popt + popt, pcov = scpopt.curve_fit(func_sca, x, spectra[ii,:], + jac=func_sca_jac, + p0=p0, bounds=bounds, + max_nfev=max_nfev, xtol=xtol, + x_scale='jac', + verbose=verbose) + lch.append(ii) + except Exception as err: + print(str(err)) + import ipdb + ipdb.set_trace() + raise err + + out = np.split(popt, indsplit) + outstd = np.split(np.sqrt(np.diag(pcov)), indsplit) + if forcelamb is True: + amp[ii, :], sigma[ii, :], bck[ii] = out + ampstd[ii, :], sigmastd[ii, :], bckstd[ii] = outstd + else: + amp[ii, :], dlamb[ii, :], sigma[ii, :], bck[ii] = out + ampstd[ii,:], dlambstd[ii,:], sigmastd[ii,:], bckstd[ii] = outstd + fit[ii, :] = func_sca(x, *popt) + p0[:] = popt[:] + + if plot_debug and ii in [0,1]: + fit = func_vect(x, amp[ii,:], x0[ii,:], sigma[ii,:], bck0[ii]) + + plt.figure() + ax0 = plt.subplot(2,1,1) + ax1 = plt.subplot(2,1,2, sharex=ax0, sharey=ax0) + ax0.plot(x,spectra[ii,:], '.k', + x, np.sum(fit, axis=0), '-r') + ax1.plot(x, fit.T) + + std = np.sqrt(np.sum((spectra-fit)**2, axis=1)) + + dout = {'fit': fit, 'lamb0': lamb0, 'std': std, 'lch': lch, + 'amp': amp, 'ampstd': ampstd, + 'sigma': sigma, 'sigmastd': sigmastd, + 'bck': bck, 'bckstd': bckstd, + 'dlamb': dlamb, 'dlambstd': dlambstd} + + return dout + + + +########################################################### +########################################################### +# +# 1d spectral fitting with physics parameters +# +########################################################### +########################################################### + + +def get_lamb0_from_dlines(dlines): + lamb0, ions = zip(*[(vv['lamb0'], + np.full((len(vv['lamb0']),), kk)) + for kk, vv in dlines.items()]) + lamb0 = np.r_[lamb0] + ind = np.argsort(lamb0) + return lamb0[ind], np.concatenate(ions)[ind] + + +def get_dindx(bckdeg=None, dlines=None, nbs=None): + + nbck = bckdeg + 1 + if nbs is None: + # 1d spectral fit + nbs = 1 + + i0 = nbck + lk = ['sigma', 'dlamb', 'amp', 'ntot', 'nlamb'] + dindx= {'bck': np.r_[:nbck], + 'ions':dict.fromkeys(sorted(dlines.keys())), + 'nbs': nbs} + for kk in dindx['ions'].keys(): + dindx['ions'][kk] = dict.fromkeys(lk) + dindx['ions'][kk]['sigma'] = i0 + np.r_[:nbs] + dindx['ions'][kk]['dlamb'] = i0+nbs + np.r_[:nbs] + nlamb = len(dlines[kk]['lamb0']) + dindx['ions'][kk]['amp'] = i0+2*nbs + np.r_[:nlamb*nbs] + dindx['ions'][kk]['nlamb'] = nlamb + dindx['ions'][kk]['ntot'] = (2 + nlamb)*nbs + i0 += dindx['ions'][kk]['ntot'] + dindx['nall'] = i0 + return dindx + + +def get_x0_bounds(x01d=None, dlines=None, dindx=None, + lamb=None, data=None): + + x0 = np.zeros((dindx['nall'],), dtype=float) + if x01d is None: + # Get average spectral width and separation + lamb0_Delta = lamb0.max() - lamb0.min() + nlamb0 = lamb0.size + lamb0_delta = lamb0_Delta / nlamb0 + + nbs = dindx['nbs'] + + x0[dindx['bck']] = np.zeros((dindx['bck'].size,)) + for kk in dindx['ions'].keys(): + # sigma + x0[dindx[kk]['sigma']] = lamb0_delta + # dlamb + x0[dindx[kk]['dlamb']] = 0. + # amp + x0[dindx[kk]['amp']] = ampmean + + else: + x0[dindx['bck']] = x01d[dindx['bck']] + i0 = dindx['bck'].size + for kk in dindx['ions'].keys(): + # TBF + # x0[dindx[kk]['sigma']] = x01d[] + pass + + # Get bounds + lamb_delta = np.mean(np.abs(np.diff(np.unique(lamb)))) + datamax = np.nanmax(data) + bampup = min(datamax, np.nanmean(data) + np.nanstd(data)) + + bounds0 = np.zeros((dindx['nall'],), dtype=float) + bounds1 = np.zeros((dindx['nall'],), dtype=float) + if dindx['bck'].size == 1: + bounds0[dindx['bck']] = 0. + bounds1[dindx['bck']] = bampup + elif dindx['bck'].size == 2: + bounds0[dindx['bck'][0]] = 0. + bounds1[dindx['bck'][0]] = bampup + bounds0[dindx['bck'][0]] = 0. # TBC + bounds1[dindx['bck'][0]] = bampup # TBC + for kk in dindx['ions'].keys(): + bounds0[dindx[kk]['sigma']] = 2.*lamb_delta + bounds1[dindx[kk]['sigma']] = lamb0_delta*5. + bounds0[dindx[kk]['dlamb']] = -3.*lamb0_delta + bounds1[dindx[kk]['dlamb']] = 3.*lamb0_delta + bounds0[dindx[kk]['amp']] = 0. + bounds1[dindx[kk]['amp']] = datamax + + return x0, bounds + + + + +def get_funccostjac(): + + def func(): + pass + + def cost(): + pass + + def jac(): + pass + + return func, cost, jac + + + + + + + + + + +########################################################### +########################################################### +# +# 2d spectral fitting +# +########################################################### +########################################################### + +def get_knots_nbs_for_bsplines(knots_unique, deg): + if deg > 0: + knots = np.r_[[knots_unique[0]]*deg, knots_unique, + [knots_unique[-1]]*deg] + else: + knots = knots_unique + nbknotsperbs = 2 + deg + nbs = knots_unique.size - 1 + deg + assert nbs == knots.size - 1 - deg + return knots, nbknotsperbs, nbs + + +def get_2dspectralfit_func(lamb0, forcelamb=False, + deg=None, knots=None): + + lamb0 = np.atleast_1d(lamb0).ravel() + nlamb = lamb0.size + knots = np.atleast_1d(knots).ravel() + nknots = knots.size + nbsplines = np.unique(knots).size - 1 + deg + + # Define function + def func(lamb, phi, + camp=None, cwidth=None, cshift=None, + lamb0=lamb0, nlamb=nlamb, + knots=knots, deg=deg, forcelamb=forcelamb, + nbsplines=nbsplines, mesh=True): + assert phi.ndim in [1, 2] + if camp is not None: + assert camp.shape[0] == nbsplines + bsamp = BSpline(knots, camp, deg, + extrapolate=False, axis=0) + if csigma is not None: + assert csigma.shape[0] == nbsplines + bssigma = BSpline(knots, csigma, deg, + extrapolate=False, axis=0) + if mesh or phi.ndim == 2: + lamb0 = lamb0[None, None, :] + else: + lamb0 = lamb0[None, :] + if forcelamb: + if mesh: + assert angle.ndim == lamb.ndim == 1 + # shape (lamb, angle, lines) + return np.sum(bsamp(phi)[None,:,:] + * np.exp(-(lamb[:,None,None] + - lamb0)**2 + /(bssigma(phi)[None,:,:]**2)), axis=-1) + else: + assert phi.shape == lamb.shape + lamb = lamb[..., None] + # shape (lamb/angle, lines) + return np.sum(bsamp(phi) + * np.exp(-(lamb + - lamb0)**2 + /(bssigma(phi)**2)), axis=-1) + else: + if cdlamb is not None: + assert cdlamb.shape[0] == nbsplines + bsdlamb = BSpline(knots, cdlamb, deg, + extrapolate=False, axis=0) + + return func + + +def get_multigaussianfit2d_costfunc(lamb=None, phi=None, data=None, std=None, + lamb0=None, forcelamb=None, + deg=None, knots=None, + nlamb0=None, nkperbs=None, nbs=None, + nc=None, debug=None): + assert lamb.shape == phi.shape == data.shape + assert lamb.ndim == 1 + assert nc == nbs*nlamb0 + + if forcelamb is None: + forcelamb = False + if debug is None: + debug = False + + # Define func assuming all inpus properly formatted + if forcelamb: + # x = [camp[1-nbs,...,nbs*(nlamb0-1)-nc}, csigma[1-nc]] + def func(x, + lamb=lamb, phi=phi, data=data, std=std, + lamb0=lamb0, knots=knots, deg=deg, nc=nc): + amp = BSpline(knots, x[:nc], deg, + extrapolate=False, axis=0)(phi) + sigma = BSpline(knots, x[nc:], deg, + extrapolate=False, axis=0)(phi) + val = np.sum(amp[:, None] + * np.exp(-(lamb[:, None] - lamb0[None, :])**2 + /(sigma[:, None]**2)), axis=-1) + return (val-data)/(std*data.size) + + def jac(x, + lamb=lamb, phi=phi, std=std, + lamb0=lamb0, knots=knots, deg=deg, + nlamb0=nlamb0, nkperbs=nkperbs, nbs=nbs, nc=nc): + amp = BSpline(knots, x[:nc], deg, + extrapolate=False, axis=0)(phi) + sigma = BSpline(knots, x[nc:], deg, + extrapolate=False, axis=0)(phi) + jacx = sparse.csr_matrix((phi.size, 2*nc), dtype=float) + #jacx = np.zeros((phi.size, 2*nc), dtype=float) + for ii in range(nlamb0): + expi = np.exp(-(lamb-lamb0[ii])**2/sigma**2) + for jj in range(nbs): + ind = ii*nbs + jj + indk = np.r_[jj*nkperbs:(jj+1)*nkperbs] + # all bsplines are the same, only coefs (x) are changing + bj = BSpline.basis_element(knots[indk], + extrapolate=False)(phi) + #bj[np.isnan(bj)] = 0. + indok = ~np.isnan(bj) + # Differentiate wrt camp + jacx[indok, ind] = (bj * expi)[indok] + # Differentiate wrt csigma + jacx[indok, nc+ind] = ( + amp * (2*(lamb-lamb0[ii])**2*bj/sigma**3) * expi + )[indok] + return jacx/(std*phi.size) + else: + # x = [camp1-nbs*nlamb, csigma1-nbs*nlamb, cdlamb1-nbs*nlamb] + def func(x, + lamb=lamb, phi=phi, data=data, std=std, + lamb0=lamb0, knots=knots, deg=deg, + nbs=nbs, nlamb0=nlamb0, nc=nc, debug=debug): + amp = BSpline(knots, x[:nc].reshape((nbs, nlamb0), order='F'), + deg, extrapolate=False, axis=0)(phi) + sigma = BSpline(knots, x[nc:2*nc].reshape((nbs, nlamb0), order='F'), + deg, extrapolate=False, axis=0)(phi) + dlamb = BSpline(knots, x[2*nc:-1].reshape((nbs, nlamb0), order='F'), + deg, extrapolate=False, axis=0)(phi) + val = np.nansum(amp + * np.exp(-(lamb[:, None] - (lamb0[None, :]+dlamb))**2 + / sigma**2), + axis=-1) + x[-1] + if debug: + vmin, vmax = 0, np.nanmax(data) + fig = plt.figure(figsize=(14, 10)); + ax0 = fig.add_axes([0.05,0.55,0.25,0.4]) + ax1 = fig.add_axes([0.35,0.55,0.25,0.4], sharex=ax0, sharey=ax0) + ax2 = fig.add_axes([0.65,0.55,0.25,0.4], sharex=ax0, sharey=ax0) + ax3 = fig.add_axes([0.05,0.05,0.25,0.4], sharex=ax0, sharey=ax0) + ax4 = fig.add_axes([0.35,0.05,0.25,0.4], sharex=ax0, sharey=ax0) + ax5 = fig.add_axes([0.65,0.05,0.25,0.4], sharex=ax0, sharey=ax0) + ax0.scatter(lamb, phi, c=data, s=2, marker='s', edgecolors='None', + vmin=vmin, vmax=vmax) # DB + ax1.scatter(lamb, phi, c=val, s=2, marker='s', edgecolors='None', # DB + vmin=vmin, vmax=vmax) # DB + errmax = np.nanmax(np.abs((val-data) / (std*data.size))) + ax2.scatter(lamb, phi, c=(val-data) / (std*data.size), + s=2, marker='s', edgecolors='None', # DB + vmin=-errmax, vmax=errmax, cmap=plt.cm.seismic) # DB + dlamb0_amp = np.max(np.diff(lamb0))/np.nanmax(np.abs(amp)) + dlamb0_sigma = np.max(np.diff(lamb0))/np.nanmax(np.abs(sigma)) + dlamb0_dlamb = np.max(np.diff(lamb0))/np.nanmax(np.abs(dlamb)) + for ii in range(nlamb0): + ax3.axvline(lamb0[ii], ls='--', c='k') + ax4.axvline(lamb0[ii], ls='--', c='k') + ax5.axvline(lamb0[ii], ls='--', c='k') + ax3.plot(lamb0[ii] + dlamb0_amp*amp[:, ii], phi, '.', ms=4) + ax4.plot(lamb0[ii] + dlamb0_sigma*sigma[:, ii], phi, '.', ms=4) + ax5.plot(lamb0[ii] + dlamb0_dlamb*dlamb[:, ii], phi, '.', ms=4) + import ipdb # DB + ipdb.set_trace() # DB + return (val-data) / (std*data.size) + + def jac(x, + lamb=lamb, phi=phi, std=std, + lamb0=lamb0, knots=knots, deg=deg, + nlamb0=nlamb0, nkperbs=nkperbs, nbs=nbs, nc=nc): + amp = BSpline(knots, x[:nc], deg, + extrapolate=False, axis=0)(phi) + sigma = BSpline(knots, x[nc:2*nc], deg, + extrapolate=False, axis=0)(phi) + dlamb = BSpline(knots, x[2*nc:], deg, + extrapolate=False, axis=0)(phi) + #jacx = sparse.csr_matrix((phi.size, 2*nc), dtype=float) + jacx = np.zeros((phi.size, 3*nc+1), dtype=float) + for ii in range(nlamb0): + expi = np.exp(-(lamb-(lamb0[ii]+dlamb))**2/sigma**2) + indlamb = expi > 0.001 + for jj in range(nbs): + kk = ii*nbs + jj + indk = jj + np.r_[:nkperbs] + # all bsplines are the same, only coefs (x) are changing + bj = BSpline.basis_element(knots[indk], + extrapolate=False)(phi) + # bj[np.isnan(bj)] = 0. + indok = (~np.isnan(bj)) & indlamb + # Differentiate wrt camp + jacx[indok, kk] = (bj[indok] * expi[indok]) + # Differentiate wrt csigma + jacx[indok, nc+kk] = ( + amp * 2*(lamb-(lamb0[ii]+dlamb))**2*bj/sigma**3 * expi + )[indok] + # Differentiate wrt dlamb + jacx[indok, 2*nc+kk] = ( + amp * 2*(lamb-(lamb0[ii]+dlamb))*bj/sigma**2 * expi + )[indok] + jacx[:, -1] = 1. + return jacx/(std*phi.size) + return func, jac + +def multigaussianfit2d(lamb, phi, data, std=None, + lamb0=None, forcelamb=None, + knots=None, deg=None, nbsplines=None, + x0=None, bounds=None, + method=None, max_nfev=None, + xtol=None, ftol=None, gtol=None, + loss=None, verbose=0, debug=None): + + # Check inputs + if deg is None: + deg = 3 + if nbsplines is None: + nbsplines = 5 + if method is None: + method = 'trf' + # Only 2 method for pb. with bounds + assert method in ['trf', 'dogbox'], method + if xtol is None: + xtol = 1.e-6 + if ftol is None: + ftol = 1.e-6 + if gtol is None: + gtol = 1.e-8 + if loss is None: + loss = 'linear' + if max_nfev is None: + max_nfev = None + if std is None: + std = 0.1*np.nanmean(data) + assert lamb0 is not None + + # Get knots + if knots is None: + phimin, phimax = np.nanmin(phi), np.nanmax(phi) + knots = np.linspace(phimin, phimax, nbsplines+1-deg) + knots, nkperbs, nbs = get_knots_nbs_for_bsplines(np.unique(knots), deg) + + # Scaling + lambmin = np.nanmin(lamb) + lamb0Delta = np.max(lamb0) - np.min(lamb0) + nlamb0 = np.size(lamb0) + nc = nbs*nlamb0 + + dlambscale = lamb0Delta / nlamb0 + ampscale = np.nanmean(data) + np.nanstd(data) + + datascale = data / ampscale + lambscale = lamb / dlambscale + lamb0scale = lamb0 / dlambscale + stdscale = std / ampscale + + # Get cost function and jacobian + func, jac = get_multigaussianfit2d_costfunc(lamb=lambscale, + phi=phi, + data=datascale, + std=stdscale, + lamb0=lamb0scale, + forcelamb=forcelamb, + deg=deg, knots=knots, + nlamb0=nlamb0, nbs=nbs, + nkperbs=nkperbs, nc=nc, + debug=debug) + + # Get initial guess + if x0 is None: + x0 = np.r_[np.ones((nc,)), np.ones((nc,))] + if not forcelamb: + x0 = np.r_[x0, np.zeros((nc,))] + x0 = np.r_[x0, 0.] + + # Get bounds + if bounds is None: + bounds = (np.r_[np.zeros((nc,)), + np.full((nc,), nlamb0/100)], + np.r_[np.full((nc,), np.nanmax(data)/ampscale), + np.full((nc,), 3.)]) + if not forcelamb: + bounds = (np.r_[bounds[0], -np.full((nc,), 2.)], + np.r_[bounds[1], np.full((nc,), 2.)]) + bounds = (np.r_[bounds[0], 0.], + np.r_[bounds[1], 0.1*np.nanmax(data)/ampscale]) + + # Minimize + res = scpopt.least_squares(func, x0, jac=jac, bounds=bounds, + method=method, ftol=ftol, xtol=xtol, + gtol=gtol, x_scale=1.0, f_scale=1.0, loss=loss, + diff_step=None, tr_solver=None, + tr_options={}, jac_sparsity=None, + max_nfev=max_nfev, verbose=verbose, + args=(), kwargs={}) + + # Separate and reshape output + camp = res.x[:nc].reshape((nlamb0, nbs)) * ampscale + csigma = res.x[nc:2*nc].reshape((nlamb0, nbs)) * dlambscale + if forcelamb: + cdlamb = None + else: + cdlamb = res.x[2*nc:3*nc].reshape((nlamb0, nbs)) * dlambscale + + # Create output dict + dout = {'camp': camp, 'csigma': csigma, 'cdlamb': cdlamb, 'bck':res.x[-1], + 'fit':(func(res.x)*stdscale*data.size + datascale) * ampscale, + 'lamb0':lamb0, 'knots': knots, 'deg':deg, 'nbsplines': nbsplines, + 'cost': res.cost, 'fun': res.fun, 'active_mask': res.active_mask, + 'nfev': res.nfev, 'njev': res.njev, 'status': res.status} + return dout + + + + +########################################################### +# +# From DataCam2D +# +########################################################### +########################################################### + + +# DEPRECATED +def fit_spectra2d_x0_per_row(): + + # Loop from centre to edges + for jj in range(ny): + out = multiplegaussianfit(x, datat[jj,:], nmax=nmax, p0=p0u, bounds=None, + max_nfev=None, xtol=1.e-8, verbose=0, + percent=20, plot_debug=False) + x0[jj,:], x0_std[jj,:] = out[:2] + + for jj in range(nybis): + # Upper half + ju = indy1[jj] + out = multiplegaussianfit(x, spect, nmax=nmax, p0=p0u, bounds=None, + max_nfev=None, xtol=1.e-8, verbose=0, + percent=20, plot_debug=False) + x0[ju,:], x0_std[ju,:] = out[:2] + p0u[:nmax], p0u[nmax:2*nmax], = amp[ii,ju,:], x0[ii,ju,:] + p0u[2*nmax:3*nmax], p0u[3*nmax:] = sigma[ii,ju,:], bck0 + + # Lower half + jl = indy2[jj] + return x0 + + + +def get_func2d(y0, y1, x0_y, bspl_n, bspl_deg): + knots = np.linspace(y0,y1, 6) + bspliney = scpinterp.LSQUnivariateSpline() + def func(x, y, ampy_coefs, sigy_coefs, bcky_coefs): + amp_bs = BSpline(knots, ampy_coefs, k=bspl_deg, + extrapolate=False, axis=0) + amp = amp_bs(y) + x0y = x0_y(y) + return np.sum(amp*np.exp((x-xoy)**2/sig**2) + bck0, axis=-1) + return func + + + +def fit_spectra_2d(data2d, indt=None, nbin_init=None, + nmax=None, bck=None, nbsplines=None): + """ Return fitted spectra + + Can handle unique or multiple time + Takes already formatted 2d data: + - (nx, ny) + - (nt, nx, ny) + x being the horizontal / spectral direction (lambda) + + """ + + ##################### + # Check / format input + ##################### + + # Check data + assert isinstance(data, np.ndarray) + assert data.ndim in [2,3] + if data.ndim == 2: + data = np.reshape((1,data.shape[0],data.shape[1])) + if indt is not None: + data = data[indt,...] + + # Set bck type + if bck is None: + bck = 0 + assert type(bck) in [int, str] + if type(bck) is int: + nbck = bck + 1 + elif bck == 'exp': + nbck = 3 + + # Extract shape + nt = data.shape[0] + nlamb, ny = data.shape[1:] + x = np.arange(0,nlamb) + + # max number of spectral lines (gaussians) + if nmax is None: + nmax = 10 + + # Check nbin_init vs ny + if nbin_init is None: + nbin_init = 100 + if ny % 2 != nbin_init % 2: + nbin_init += 1 + + # get indybin + indybin = np.arange(0, nbin_init) + if ny % 2 == 0: + indybin += int(ny/2 - nbin_init/2) + else: + indybin += int((ny-1)/2 - (nbin_init-1)/2) + + # get indybis + if ny % 2 == 0: + indy1 = np.arange(ny/2-1, -1, -1) + indy2 = np.arange(ny/2, ny, 1) + else: + indy1 = np.arange((ny-1)/2-1, -1, -1) + indy2 = np.arange((ny-1)/2+1, ny, 1) + nybis = indy1.size + assert nybis == indy2.size + + ##################### + # initialize output + ##################### + + # results + amp = np.full((nt, ny, nmax), np.nan) + x0 = np.full((ny, nmax), np.nan) + sigma = np.full((nt, ny, nmax), np.nan) + bck0 = np.full((nt, ny, nbck), np.nan) + p0u = np.full((nmax*3 + nbck,), np.nan) + p0l = np.full((nmax*3 + nbck,), np.nan) + + # results std + amp_std = np.full((nt, ny, nmax), np.nan) + x0_std = np.full((nt, ny, nmax), np.nan) + sigma_std = np.full((nt, ny, nmax), np.nan) + bck0_std = np.full((nt, ny, nbck), np.nan) + std = np.full((nt,), np.nan) + + + ##################### + # Determine x0 for each row + ##################### + + datat = np.sum(data, axis=0) + + # Start with central binned (initial guess) + + # if ny is odd => use binned result for central line + + + x0 = None + + # Deduce ellipses + + # Derive x0_func(y, x0) + + + ##################### + # Define 2d multi-gaussian + ##################### + + # Define x0(y) law + + + + + ##################### + # compute fits, with fixed x0 + ##################### + + # Get binned spectra to deduce initial guess + databin = None + + # loop over y to fit spectra, start from middle, extent to edges + for ii in range(nt): + + out = fit2dspectra(x, y, data[ii,:,:]) + + p0u[:nmax] = amp[ii-1,indy1[0],:] + p0u[nmax:2*nmax] = x0[ii-1,indy1[0],:] + p0u[2*nmax:3*nmax] = sigma[ii,indy1[0],:] + p0u[3*nmax:] = bck0[ii,indy1[0]] + + p0l[:nmax] = amp[ii-1,indy2[0],:] + p0l[nmax:2*nmax] = x0[ii-1,indy2[0],:] + p0l[2*nmax:3*nmax] = sigma[ii,indy2[0],:] + p0l[3*nmax:] = bck0[ii,indy2[0]] + + + return (x0, x0_std), ( ) + + +########################################################### +########################################################### +# +# For inspiration +# +########################################################### +########################################################### + + + +# _indymod = {0: np.arange(0,195), + # 1: np.arange(212,407), + # 2: np.arange(424,619), + # 3: np.arange(636,831), + # 4: np.arange(848,1043), + # 5: np.arange(1060,1255), + # 6: np.arange(1272,1467)} + +# def get_binnedspectra(imp='FeXXV', mod=3, path='./', + # dshots=_dshots, indymod=_indymod, save=True): + + # # Prepare + # indy = indymod[mod] + # dimp = dict(dshots[imp]) + + # # Loop on shots + # lextra = ['Te0','ne0','dg16'] + # ls = sorted(dimp.keys()) + # for shot in ls: + + # # Load data + # if 'tlim' in dimp[shot].keys(): + # tlim = dimp[shot]['tlim'] + # else: + # tlim = None + # try: + # xics, kh = tfw.SpectroX2D.load_data(shot, + # tlim=tlim, plot=False) + # except Exception as err: + # dimp[shot]['err'] = str(err) + # continue + + # # Bin data + # data, t = xics._Ref['data'], xics._Ref['t'] + # dimp[shot]['nt'] = t.size + # spectra = np.empty((dimp[shot]['nt'], 487), dtype=float) + # for ii in range(0,dimp[shot]['nt']): + # spectra[ii,:] = np.nanmean(data[ii,:].reshape(1467,487)[indy,:], + # axis=0) + # dimp[shot]['spectra'] = spectra + # dimp[shot]['t'] = t + # #dimp[shot]['date'] = IRFMtb. + + # # Dextra + # for ss in lextra: + # try: + # indt = np.digitize(xics.dextra[ss]['t'], 0.5*(t[1:]+t[:-1])) + # val = np.empty((dimp[shot]['nt'],),dtype=float) + # std = np.empty((dimp[shot]['nt'],),dtype=float) + # ssum = np.empty((dimp[shot]['nt'],),dtype=float) + # for ii in range(0,dimp[shot]['nt']): + # val[ii] = np.nanmean(xics.dextra[ss]['data'][indt==ii]) + # std[ii] = np.nanstd(xics.dextra[ss]['data'][indt==ii]) + # ssum[ii] = np.nansum(xics.dextra[ss]['data'][indt==ii]) + # dimp[shot][ss] = {'mean':val, 'std':std, 'sum':ssum} + # except Exception as err: + # dimp[shot][ss] = {'mean':np.nan, 'std':np.nan, 'sum':np.nan} + + # # Reshape dict for np.savez and pandas DataFrame + # nt = np.array([dimp[shot]['nt'] for shot in ls]) + # nttot = np.sum(nt) + # ntcum = np.cumsum(nt) + # lk = ['shot','angle','spectra','t', + # 'Te0-mean','Te0-std','ne0-mean','ne0-std', + # 'dg16-sum'] + # dk = {} + # for k in lk: + # shape = (nttot,487) if k=='spectra' else (nttot,) + # dk[k] = np.full(shape, np.nan) + + # i0 = 0 + # for ii in range(0,len(ls)): + # ind = np.arange(i0,i0+nt[ii]) + # dk['shot'][ind] = ls[ii] + # dk['angle'][ind] = dimp[ls[ii]]['ang'] + # dk['spectra'][ind,:] = dimp[ls[ii]]['spectra'] + # dk['t'][ind] = dimp[ls[ii]]['t'] + # dk['Te0-mean'][ind] = dimp[ls[ii]]['Te0']['mean'] + # dk['Te0-std'][ind] = dimp[ls[ii]]['Te0']['std'] + # dk['ne0-mean'][ind] = dimp[ls[ii]]['ne0']['mean'] + # dk['ne0-std'][ind] = dimp[ls[ii]]['ne0']['std'] + # dk['dg16-sum'][ind] = dimp[ls[ii]]['dg16']['sum'] + # i0 = ntcum[ii] + + # # Saving + # if save: + # name = '%s_spectra'%imp + # path = os.path.abspath(path) + # pfe = os.path.join(path,name+'.npz') + # try: + # np.savez(pfe, **dk) + # print("Saved in :", pfe) + # except: + # import ipdb + # ipdb.set_trace() + # pass + # return dk + + + +# #################################################################### +# # spectral fit +# #################################################################### + + +# def remove_bck(x, y): + # #opt = np.polyfit(x, y, deg=0) + # opt = [np.nanmin(y)] + # ybis = y - opt[0] + # return ybis, opt[0] + + +# def get_peaks(x, y, nmax=10): + + # # Prepare + # ybis = np.copy(y) + # A = np.empty((nmax,),dtype=y.dtype) + # x0 = np.empty((nmax,),dtype=x.dtype) + # sigma = np.empty((nmax,),dtype=y.dtype) + # gauss = lambda xx, A, x0, sigma: A*np.exp(-(xx-x0)**2/sigma**2) + # def gauss_jac(xx, A, x0, sigma): + # jac = np.empty((xx.size,3),dtype=float) + # jac[:,0] = np.exp(-(xx-x0)**2/sigma**2) + # jac[:,1] = A*2*(xx-x0)/sigma**2 * np.exp(-(xx-x0)**2/sigma**2) + # jac[:,2] = A*2*(xx-x0)**2/sigma**3 * np.exp(-(xx-x0)**2/sigma**2) + # return jac + + # dx = np.nanmin(np.diff(x)) + + # # Loop + # nn = 0 + # while nn=0.): + # wp = min(x.size-1, + # ind + np.nonzero(np.diff(ybis[ind:],n=2)>=0.)[0][0] + 1) + # else: + # wp = ybis.size-1 + # if np.any(np.diff(ybis[:ind+1],n=2)>=0.): + # wn = max(0, np.nonzero(np.diff(ybis[:ind+1],n=2)>=0.)[0][-1] - 1) + # else: + # wn = 0 + # width = x[wp]-x[wn] + # assert width>0. + # indl = np.arange(wn,wp+1) + # sig = np.ones((indl.size,)) + # if (np.abs(np.mean(np.diff(ybis[ind:wp+1]))) + # > np.abs(np.mean(np.diff(ybis[wn:ind+1])))): + # sig[indlind] = 0.5 + # else: + # sig[indlind] = 1.5 + # p0 = (ybis[ind],x00,width)#,0.) + # bounds = (np.r_[0.,x[wn],dx/2.], + # np.r_[5.*ybis[ind],x[wp],5.*width]) + # try: + # (Ai, x0i, sigi) = scpopt.curve_fit(gauss, x[indl], ybis[indl], + # p0=p0, bounds=bounds, jac=gauss_jac, + # sigma=sig, x_scale='jac')[0] + # except Exception as err: + # print(str(err)) + # import ipdb + # ipdb.set_trace() + # pass + + # ybis = ybis - gauss(x, Ai, x0i, sigi) + # A[nn] = Ai + # x0[nn] = x0i + # sigma[nn] = sigi + + + # nn += 1 + # return A, x0, sigma + +# def get_p0bounds(x, y, nmax=10): + + # yflat, bck = remove_bck(x,y) + # A, x0, sigma = get_peaks(x, yflat, nmax=nmax) + + # p0 = A.tolist() + x0.tolist() + sigma.tolist() + [bck] + + # lx = [np.nanmin(x), np.nanmax(x)] + # Dx = np.diff(lx) + # dx = np.nanmin(np.diff(x)) + + # bA = (np.zeros(nmax,), np.full((nmax,),3.*np.nanmax(y))) + # bx0 = (np.full((nmax,),lx[0]-Dx/2.), np.full((nmax,),lx[1]+Dx/2.)) + # bsigma = (np.full((nmax,),dx/2.), np.full((nmax,),Dx/2.)) + # bbck0 = (0., np.nanmax(y)) + + # bounds = (np.r_[bA[0],bx0[0],bsigma[0], bbck0[0]], + # np.r_[bA[1],bx0[1],bsigma[1], bbck0[1]]) + # if not np.all(bounds[0] spectrum {0} / {1}".format(ii,nspect)) + + # try: + # popt, pcov = scpopt.curve_fit(func_sca, x, spectra[ii,:], + # jac=func_sca_jac, + # p0=p0, bounds=bounds, + # max_nfev=max_nfev, xtol=xtol, + # x_scale='jac', + # verbose=verbose) + # except Exception as err: + # msg = " Convergence issue for {0} / {1}\n".format(ii,nspect) + # msg += " => %s\n"%str(err) + # msg += " => Resetting initial guess and bounds..." + # print(msg) + # try: + # p0, bounds = get_p0bounds(x, spectra[ii,:], nmax=nmax) + # popt, pcov = scpopt.curve_fit(func_sca, x, spectra[ii,:], + # jac=func_sca_jac, + # p0=p0, bounds=bounds, + # max_nfev=max_nfev, xtol=xtol, + # x_scale='jac', + # verbose=verbose) + # p0 = popt + # popt, pcov = scpopt.curve_fit(func_sca, x, spectra[ii,:], + # jac=func_sca_jac, + # p0=p0, bounds=bounds, + # max_nfev=max_nfev, xtol=xtol, + # x_scale='jac', + # verbose=verbose) + # lch.append(ii) + # except Exception as err: + # print(str(err)) + # import ipdb + # ipdb.set_trace() + # continue + + + # A[ii,:] = popt[:nmax] + # x0[ii,:] = popt[nmax:2*nmax] + # sigma[ii,:] = popt[2*nmax:3*nmax] + # bck0[ii] = popt[3*nmax] + + # stdi = np.sqrt(np.diag(pcov)) + # Astd[ii,:] = stdi[:nmax] + # x0std[ii,:] = stdi[nmax:2*nmax] + # sigmastd[ii,:] = stdi[2*nmax:3*nmax] + # bck0std[ii] = stdi[3*nmax] + # std[ii] = np.sqrt(np.sum((spectra[ii,:]-func_sca(x,*popt))**2)) + + # p0[:] = popt[:] + + # if plot_debug and ii in [0,1]: + # fit = func_vect(x, A[ii,:], x0[ii,:], sigma[ii,:], bck0[ii]) + + # plt.figure() + # ax0 = plt.subplot(2,1,1) + # ax1 = plt.subplot(2,1,2, sharex=ax0, sharey=ax0) + # ax0.plot(x,spectra[ii,:], '.k', + # x, np.sum(fit,axis=0), '-r') + # ax1.plot(x, fit.T) + + # return (A,Astd), (x0,x0std), (sigma,sigmastd), (bck0,bck0std), std, lch + + + +# def add_gaussianfits(dimp, nmax=10, verbose=0, percent=20, + # path='./', save=False): + # assert type(dimp) in [str,dict] + + # # Prepare + # if type(dimp) is str: + # imp = str(dimp) + # if '_' in imp: + # imp = imp.split('_')[0] + # dimp = dict(np.load(dimp+'_spectra.npz')) + # inds = np.argsort(dimp['angle']) + # for k in dimp.keys(): + # if inds.size in dimp[k].shape: + # dimp[k] = dimp[k][inds] if dimp[k].ndim==1 else dimp[k][inds,:] + # else: + # imp = None + # save = False + + + # # Compute + # ind = np.arange(3,437) + # x = np.arange(0,dimp['spectra'].shape[1]) + # spectra = dimp['spectra'][:,ind] + # A, x0, sigma, bck0, std, lch = multiplegaussianfit(x[ind], spectra, + # nmax=nmax, percent=percent, + # verbose=verbose) + # # Store + # dimp['nmax'] = nmax + # dimp['indch'] = np.r_[lch] + # dimp['x'] = x + # dimp['ind'] = ind + # dimp['A'] = A[0] + # dimp['A-std'] = A[1] + # dimp['x0'] = x0[0] + # dimp['x0-std'] = x0[1] + # dimp['sigma'] = sigma[0] + # dimp['sigma-std'] = sigma[1] + # dimp['bck0'] = bck0[0] + # dimp['bck0-std'] = bck0[1] + # dimp['std'] = std + # if imp is not None: + # dimp['imp'] = imp + + # if save: + # name = '{0}_fitted{1}'.format(imp,nmax) + # path = os.path.abspath(path) + # pfe = os.path.join(path,name+'.npz') + # np.savez(pfe, **dimp) + # print("Saved in :", pfe) + + # return dimp + + + +# def plot_gaussianfits(dimp, ind): + + # # Prepare + # x = dimp['x'] + # spectra = dimp['spectra'][ind,:] + # func_vect, func_sca, func_sca_jac = get_func(dimp['nmax']) + # fit = func_vect(x, dimp['A'][ind,:], dimp['x0'][ind,:], + # dimp['sigma'][ind,:], dimp['bck0'][ind]) + + # # Plot + # plt.figure(); + # ax0 = plt.subplot(2,1,1) + # ax1 = plt.subplot(2,1,2, sharex=ax0, sharey=ax0) + + # ax0.plot(x, spectra, '.k', label='spectrum') + # ax0.plot(x, np.sum(fit,axis=0), '-r', label='fit') + # ax1.plot(x, fit.T) + + # ax1.set_xlabel(r'x') + # ax0.set_ylabel(r'data') + # ax1.set_ylabel(r'data') + + +# def plot_allraw(dimp): + + # # Prepare + # x = dimp['x'] + # spectra = dimp['spectra'] + # nspect = spectra.shape[0] + # spectranorm = spectra / np.nanmax(spectra,axis=1)[:,np.newaxis] + + # # Plot + # plt.figure(); + # ax0 = plt.subplot(2,1,1) + # ax1 = plt.subplot(2,1,2, sharex=ax0) + + # ax0.imshow(spectranorm, cmap=plt.cm.viridis, + # extent=(x.min(),x.max(),0,nspect), + # origin='lower', + # interpolation='bilinear') + + +# def extract_lines(dimp): + + # nspect, nx = dimp['spectra'].shape + # if dimp['imp'] == 'ArXVII': + # dlines = {'w': {'range':np.arange(280,nx)}, + # 'z': {'range':np.arange(0,200)}} + # for k in dlines.keys(): + # spect = np.copy(dimp['spectra']) + # spect[:,dlines[k]['range']] = 0. + # dlines[k].update({'x':np.full((nspect,),np.nan), + # 'x0':np.full((nspect,),np.nan), + # 'A':np.full((nspect,),np.nan), + # 'sigma':np.full((nspect,),np.nan)}) + # for ii in range(0,dimp['spectra'].shape[0]): + # xl = dimp['x'][np.argmax(spect[ii,:])] + # ind = np.argmin(np.abs(dimp['x0'][ii,:]-xl)) + # dlines[k]['x'][ii] = xl + # dlines[k]['x0'][ii] = dimp['x0'][ii,ind] + # dlines[k]['A'][ii] = dimp['A'][ii,ind] + # dlines[k]['sigma'][ii] = dimp['sigma'][ii,ind] + + # dimp.update(**dlines) + # return dimp diff --git a/tofu/geom/__init__.py b/tofu/geom/__init__.py index 40ed5f2e4..b2c91ec1c 100644 --- a/tofu/geom/__init__.py +++ b/tofu/geom/__init__.py @@ -11,6 +11,7 @@ from tofu.geom._core import PlasmaDomain, Ves, PFC, CoilPF, CoilCS, Config from tofu.geom._core import Rays, CamLOS1D, CamLOS2D +from tofu.geom._core_optics import * from . import utils __all__ = ['_GG', '_comp', '_plot', '_def', 'utils'] diff --git a/tofu/geom/_comp_optics.py b/tofu/geom/_comp_optics.py new file mode 100644 index 000000000..4917835b5 --- /dev/null +++ b/tofu/geom/_comp_optics.py @@ -0,0 +1,202 @@ + +import numpy as np +import scipy.interpolate as scpinterp + +# ############################################### +# ############################################### +# CrystalBragg +# ############################################### +# ############################################### + +# ############################################### +# sampling +# ############################################### + +def CrystBragg_sample_outline_sphrect(dpsi, dtheta, npsi=None, ntheta=None): + psi = dpsi*np.linspace(-1, 1., npsi) + theta = np.pi/2. + dtheta*np.linspace(-1, 1., ntheta) + psimin = np.full((ntheta,), psi[0]) + psimax = np.full((ntheta,), psi[-1]) + thetamin = np.full((npsi,), theta[0]) + thetamax = np.full((npsi,), theta[-1]) + psi = np.concatenate((psi, psimax, + psi[::-1], psimin)) + theta = np.concatenate((thetamin, theta, + thetamax, theta[::-1])) + return psi, theta + +def CrystBragg_get_noute1e2_from_psitheta(nout, e1, e2, psi, theta, + e1e2=True): + vout = ((np.cos(psi)[None, :]*nout[:, None] + + np.sin(psi)[None, :]*e1[:, None])*np.sin(theta)[None, :] + + np.cos(theta)[None, :]*e2[:, None]) + if e1e2: + ve1 = (-np.sin(psi)[None, :]*nout[:, None] + np.cos(psi)[None, :]*e1[:, None]) + ve2 = np.array([vout[1, :]*ve1[2, :] - vout[2, :]*ve1[1, :], + vout[2, :]*ve1[0, :] - vout[0, :]*ve1[2, :], + vout[0, :]*ve1[1, :] - vout[1, :]*ve1[0, :]]) + return vout, ve1, ve2 + else: + return vout + +def CrystBragg_sample_outline_plot_sphrect(center, nout, e1, e2, + rcurve, extenthalf, res=None): + dpsi, dtheta = extenthalf + if res is None: + res = np.min(extenthalf)/5. + npsi = 2*int(np.ceil(dpsi / res)) + 1 + ntheta = 2*int(np.ceil(dtheta / res)) + 1 + psi, theta = CrystBragg_sample_outline_sphrect(dpsi, dtheta, + npsi=npsi, ntheta=ntheta) + vout = CrystBragg_get_noute1e2_from_psitheta(nout, e1, e2, psi, theta, + e1e2=False) + return center[:, None] + rcurve*vout + +def CrystBragg_sample_outline_Rays(center, nout, e1, e2, + rcurve, extenthalf, + bragg, phi): + dpsi, dtheta = extenthalf + psi, theta = CrystBragg_sample_outline_sphrect(dpsi, dtheta, + npsi=3, ntheta=3) + psi = np.append(psi, [0]) + theta = np.append(theta, [np.pi/2.]) + npts = psi.size + + # add repetitions for rays + nrays = phi.size + psi = np.repeat(psi, nrays) + theta = np.repeat(theta, nrays) + + # add tiling for pts + bragg = np.tile(bragg, npts) + phi = np.tile(phi, npts) + + # Compute local vectors + vout, ve1, ve2 = CrystBragg_get_noute1e2_from_psitheta(nout, e1, e2, + psi, theta) + # Deduce D, u + D = center[:, None] + rcurve*vout + u = (-np.sin(bragg)*vect + + np.cos(bragg)*(np.cos(phi)*ve1 + np.sin(phi)*ve2)) + return D, u + + +# ############################################### +# lamb <=> bragg +# ############################################### + +def get_bragg_from_lamb(lamb, d, n=None): + """ n*lamb = 2d*sin(bragg) """ + if n is None: + n = 1 + bragg= np.full(lamb.shape, np.nan) + sin = n*lamb/(2.*d) + indok = np.abs(sin) <= 1. + bragg[indok] = np.arcsin(sin[indok]) + return bragg + +def get_lamb_from_bragg(bragg, d, n=None): + """ n*lamb = 2d*sin(bragg) """ + if n is None: + n = 1 + return 2*d*np.sin(bragg) / n + + +# ############################################### +# Approximate solution +# ############################################### + +def get_approx_detector_rel(rcurve, bragg): + + # distance crystal - det_center + det_dist = rcurve*np.sin(bragg) + + # det_nout and det_e1 in (nout, e1, e2) (det_e2 = e2) + det_nout_rel = np.r_[np.sin(bragg), -np.cos(bragg), 0.] + det_ei_rel = np.r_[np.cos(bragg), np.sin(bragg), 0] + return det_dist, det_nout_rel, det_ei_rel + +def get_det_abs_from_rel(det_dist, det_nout_rel, det_ei_rel, + summit, nout, e1, e2): + det_nout = (det_nout_rel[0]*nout + + det_nout_rel[1]*e1 + det_nout_rel[2]*e2) + det_ei = (det_ei_rel[0]*nout + + det_ei_rel[1]*e1 + det_ei_rel[2]*e2) + det_ej = np.cross(det_nout, det_ei) + det_cent = summit - det_dist*det_nout + return det_cent, det_nout, det_ei, det_ej + + +# ############################################### +# Coordinates transforms +# ############################################### + +def checkformat_vectang(Z, nn, frame_cent, frame_ang): + # Check / format inputs + nn = np.atleast_1d(nn).ravel() + assert nn.size == 3 + nn = nn / np.linalg.norm(nn) + Z = float(Z) + + frame_cent = np.atleast_1d(frame_cent).ravel() + assert frame_cent.size == 2 + frame_ang = float(frame_ang) + + return Z, nn, frame_cent, frame_ang + +def get_e1e2_detectorplane(nn, nIn): + e1 = np.cross(nn, nIn) + e1n = np.linalg.norm(e1) + if e1n < 1.e-10: + e1 = np.array([nIn[2], -nIn[1], 0.]) + e1n = np.linalg.norm(e1) + e1 = e1 / e1n + e2 = np.cross(nn, e1) + e2 = e2 / np.linalg.norm(e2) + return e1, e2 + +def calc_xixj_from_braggphi(summit, det_cent, det_nout, det_ei, det_ej, + nout, e1, e2, bragg, phi): + sp = (det_cent - summit) + vect = (-np.sin(bragg)[None, :]*nout[:, None] + + np.cos(bragg)[None, :]*(np.cos(phi)[None, :]*e1[:, None] + + np.sin(phi)[None, :]*e2[:, None])) + k = np.sum(sp*det_nout) / np.sum(vect*det_nout[:, None], axis=0) + pts = summit[:, None] + k[None, :]*vect + + xi = np.sum((pts - det_cent[:, None])*det_ei[:, None], axis=0) + xj = np.sum((pts - det_cent[:, None])*det_ej[:, None], axis=0) + return xi, xj + +def calc_braggphi_from_xixj(xi, xj, det_cent, det_ei, det_ej, + summit, nin, e1, e2): + + xi = xi[None, ...] + xj = xj[None, ...] + if xi.ndim == 1: + summit = summit[:, None] + det_cent = det_cent[:, None] + det_ei, det_ej = det_ei[:, None], det_ej[:, None] + nin, e1, e2 = nin[:, None], e1[:, None], e2[:, None] + else: + summit = summit[:, None, None] + det_cent = det_cent[:, None, None] + det_ei, det_ej = det_ei[:, None, None], det_ej[:, None, None] + nin, e1, e2 = nin[:, None, None], e1[:, None, None], e2[:, None, None] + + + pts = det_cent + xi*det_ei + xj*det_ej + + vect = pts - summit + vect = vect / np.sqrt(np.sum(vect**2, axis=0))[None, ...] + bragg = np.arcsin(np.sum(vect*nin, axis=0)) + + phi = np.arctan2(np.sum(vect*e2, axis=0), np.sum(vect*e1, axis=0)) + return bragg, phi + +def get_lambphifit(lamb, phi, nxi, nxj): + lambD = lamb.max()-lamb.min() + lambfit = lamb.min() +lambD*np.linspace(0, 1, nxi) + phiD = phi.max() - phi.min() + phifit = phi.min() + phiD*np.linspace(0, 1, nxj) + return lambfit, phifit diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 1f5fd37e0..60ce6d4bd 100644 --- a/tofu/geom/_core.py +++ b/tofu/geom/_core.py @@ -52,7 +52,6 @@ _PHITHETAPROJ_NPHI = 2000 _PHITHETAPROJ_NTHETA = 1000 _RES = 0.005 -_NTHREADS = 16 _DREFLECT = {"specular": 0, "diffusive": 1, "ccube": 2} """ @@ -4083,20 +4082,19 @@ def _complete_dX12(self, dgeom): # Test if unique starting point if dgeom["case"] in ["A", "B", "C"]: # Test if pinhole - if dgeom["case"] in ["A", "B"]: - u = dgeom["u"][:, 0:1] - sca2 = np.sum(dgeom["u"][:, 1:] * u, axis=0) ** 2 - if np.all(sca2 < 1.0 - 1.0e-9): - DDb = dgeom["D"][:, 1:] - dgeom["D"][:, 0:1] - k = np.sum( - DDb * (u - np.sqrt(sca2) * dgeom["u"][:, 1:]), axis=0 - ) - k = k / (1.0 - sca2) - if k[0] > 0 and np.allclose( - k, k[0], atol=1.0e-3, rtol=1.0e-6 - ): - pinhole = dgeom["D"][:, 0] + k[0] * u[:, 0] - dgeom["pinhole"] = pinhole + if dgeom['D'].shape[1] == 1 and dgeom['nRays'] > 1: + dgeom['pinhole'] = dgeom['D'].ravel() + elif dgeom['case'] in ['A', 'B']: + u = dgeom['u'][:, 0:1] + sca2 = np.sum(dgeom['u'][:, 1:]*u, axis=0)**2 + if np.all(sca2 < 1.0 - 1.e-9): + DDb = dgeom['D'][:, 1:]-dgeom['D'][:, 0:1] + k = np.sum(DDb*(u - np.sqrt(sca2)*dgeom['u'][:, 1:]), + axis=0) + k = k / (1.0-sca2) + if k[0] > 0 and np.allclose(k, k[0], atol=1.e-3, rtol=1.e-6): + pinhole = dgeom['D'][:, 0] + k[0]*u[:, 0] + dgeom['pinhole'] = pinhole # Test if all D are on a common plane or line va = dgeom["D"] - dgeom["D"][:, 0:1] @@ -4984,13 +4982,14 @@ def D(self): @property def u(self): - if self.isPinhole: - u = self._dgeom["pinhole"][:, np.newaxis] - self._dgeom["D"] - u = u / np.sqrt(np.sum(u ** 2, axis=0))[np.newaxis, :] - elif self._dgeom["u"].shape[1] < self._dgeom["nRays"]: - u = np.tile(self._dgeom["u"], self._dgeom["nRays"]) - else: - u = self._dgeom["u"] + if (self._dgeom['u'] is not None + and self._dgeom['u'].shape[1] == self._dgeom['nRays']): + u = self._dgeom['u'] + elif self.isPinhole: + u = self._dgeom['pinhole'][:, None] - self._dgeom['D'] + u = u / np.sqrt(np.sum(u**2, axis=0))[None, :] + elif self._dgeom['u'].shape[1] < self._dgeom['nRays']: + u = np.tile(self._dgeom['u'], self._dgeom['nRays']) return u @property @@ -6248,7 +6247,7 @@ def funcbis(*args, **kwdargs): DL = None ani = quant is None if num_threads is None: - num_threads = _NTHREADS + num_threads = _NUM_THREADS if np.all(indok): D, u = self.D, self.u diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py new file mode 100644 index 000000000..e58138373 --- /dev/null +++ b/tofu/geom/_core_optics.py @@ -0,0 +1,1119 @@ + +""" +This module is the geometrical part of the ToFu general package +It includes all functions and object classes necessary for tomography on Tokamaks +""" + +# Built-in +import sys +import os +import warnings +import copy + + +# Common +import numpy as np +import datetime as dtm +import matplotlib.pyplot as plt +import matplotlib as mpl + +# ToFu-specific +from tofu import __version__ as __version__ +import tofu.pathfile as tfpf +import tofu.utils as utils +try: + import tofu.geom._def as _def + import tofu.geom._GG as _GG + import tofu.geom._comp_optics as _comp_optics + import tofu.geom._plot_optics as _plot_optics +except Exception: + from . import _def as _def + from . import _GG as _GG + from . import _comp_optics as _comp_optics + from . import _plot_optics as _plot_optics + + + +__all__ = ['CrystalBragg'] + + + +_Type = 'Tor' +_NTHREADS = 16 + +""" +############################################################################### +############################################################################### + Ves class and functions +############################################################################### +""" + + + +class CrystalBragg(utils.ToFuObject): + """ A class defining crystals for Bragg diffraction + + A crystal can be of Type flat, cylindrical or spherical + It is characterized by its: + - geometry (Type, dimensions, curvature radii and position/orientation) + - Material and lattice + - Bragg parameters (angle vs lambda) + + + Parameters + ---------- + Id : str / tfpf.ID + A name string or a pre-built tfpf.ID class to be used to identify this + particular instance, if a string is provided, it is fed to tfpf.ID() + dgeom : dict + An array (2,N) or (N,2) defining the contour of the vacuum vessel in a + cross-section, if not closed, will be closed automatically + dspectral: str + Flag indicating whether the vessel will be a torus ('Tor') or a linear + device ('Lin') + SavePath : None / str + If provided, forces the default saving path of the object to the + provided value + + """ + + # Fixed (class-wise) dictionary of default properties + _ddef = {'Id':{'shot': 0, 'Exp': 'dummy', 'Diag': 'dummy', + 'include':['Mod', 'Cls', 'Exp', 'Diag', + 'Name', 'shot', 'version']}, + 'dgeom':{'Type': 'sph', 'Typeoutline': 'rect'}, + 'dmat':{}, + 'dbragg':{}, + 'dmisc':{'color':'k'}} + _dplot = {'cross':{'Elt':'P', + 'dP':{'color':'k','lw':2}, + 'dI':{'color':'k','ls':'--','marker':'x','ms':8,'mew':2}, + 'dBs':{'color':'b','ls':'--','marker':'x','ms':8,'mew':2}, + 'dBv':{'color':'g','ls':'--','marker':'x','ms':8,'mew':2}, + 'dVect':{'color':'r','scale':10}}, + 'hor':{'Elt':'P', + 'dP':{'color':'k','lw':2}, + 'dI':{'color':'k','ls':'--'}, + 'dBs':{'color':'b','ls':'--'}, + 'dBv':{'color':'g','ls':'--'}, + 'Nstep':50}, + '3d':{}} + _DEFLMOVEOK = ['rotate'] + _DEFLAMB = 3.971561e-10 + _DEFNPEAKS = 12 + # _DREFLECT_DTYPES = {'specular':0, 'diffusive':1, 'ccube':2} + + + # Does not exist beofre Python 3.6 !!! + def __init_subclass__(cls, color='k', **kwdargs): + # Python 2 + super(CrystalBragg,cls).__init_subclass__(**kwdargs) + # Python 3 + #super().__init_subclass__(**kwdargs) + cls._ddef = copy.deepcopy(CrystalBragg._ddef) + cls._dplot = copy.deepcopy(CrystalBragg._dplot) + cls._set_color_ddef(cls._color) + + @classmethod + def _set_color_ddef(cls, color): + cls._ddef['dmisc']['color'] = mpl.colors.to_rgba(color) + + def __init__(self, dgeom=None, dmat=None, dbragg=None, + Id=None, Name=None, Exp=None, Diag=None, shot=None, + fromdict=None, sep=None, + SavePath=os.path.abspath('./'), + SavePath_Include=tfpf.defInclude, color=None): + + # To replace __init_subclass__ for Python 2 + if sys.version[0]=='2': + self._dstrip = utils.ToFuObjectBase._dstrip.copy() + self.__class__._strip_init() + + # Create a dplot at instance level + self._dplot = copy.deepcopy(self.__class__._dplot) + + kwdargs = locals() + del kwdargs['self'] + # super() + super(CrystalBragg,self).__init__(**kwdargs) + + def _reset(self): + # super() + super(CrystalBragg,self)._reset() + self._dgeom = dict.fromkeys(self._get_keys_dgeom()) + self._dmat = dict.fromkeys(self._get_keys_dmat()) + self._dbragg = dict.fromkeys(self._get_keys_dbragg()) + self._dmisc = dict.fromkeys(self._get_keys_dmisc()) + #self._dplot = copy.deepcopy(self.__class__._ddef['dplot']) + + @classmethod + def _checkformat_inputs_Id(cls, Id=None, Name=None, + Exp=None, Diag=None, shot=None, Type=None, + include=None, + **kwdargs): + if Id is not None: + assert isinstance(Id,utils.ID) + Name, Exp, Type = Id.Name, Id.Exp, Id.Type + if Type is None: + Type = cls._ddef['dgeom']['Type'] + if Exp is None: + Exp = cls._ddef['Id']['Exp'] + if Diag is None: + Diag = cls._ddef['Id']['Diag'] + if shot is None: + shot = cls._ddef['Id']['shot'] + if include is None: + include = cls._ddef['Id']['include'] + + dins = {'Name':{'var':Name, 'cls':str}, + 'Exp': {'var':Exp, 'cls':str}, + 'Diag': {'var':Diag, 'cls':str}, + 'shot': {'var':shot, 'cls':int}, + 'Type': {'var':Type, 'in':['sph']}, + 'include':{'var':include, 'listof':str}} + dins, err, msg = cls._check_InputsGeneric(dins) + if err: + raise Exception(msg) + + kwdargs.update({'Name':Name, 'shot':shot, + 'Exp':Exp, 'Diag':Diag, 'Type':Type, + 'include':include}) + return kwdargs + + ########### + # Get largs + ########### + + @staticmethod + def _get_largs_dgeom(sino=True): + largs = ['dgeom'] + return largs + + @staticmethod + def _get_largs_dmat(): + largs = ['dmat'] + return largs + + @staticmethod + def _get_largs_dbragg(): + largs = ['dbragg'] + return largs + + @staticmethod + def _get_largs_dmisc(): + largs = ['color'] + return largs + + ########### + # Get check and format inputs + ########### + + @classmethod + def _checkformat_dgeom(cls, dgeom=None): + if dgeom is None: + return + assert isinstance(dgeom, dict) + lkok = cls._get_keys_dgeom() + assert all([isinstance(ss, str) and ss in lkok for ss in dgeom.keys()]) + for kk in cls._ddef['dgeom'].keys(): + dgeom[kk] = dgeom.get(kk, cls._ddef['dgeom'][kk]) + for kk in lkok: + dgeom[kk] = dgeom.get(kk, None) + if dgeom['center'] is not None: + dgeom['center'] = np.atleast_1d(dgeom['center']).ravel() + assert dgeom['center'].size == 3 + else: + assert dgeom['summit'] is not None + assert dgeom['rcurve'] is not None + if dgeom['summit'] is not None: + dgeom['summit'] = np.atleast_1d(dgeom['summit']).ravel() + assert dgeom['summit'].size == 3 + else: + assert dgeom['center'] is not None + assert dgeom['rcurve'] is not None + if dgeom['extenthalf'] is not None: + dgeom['extenthalf'] = np.atleast_1d(dgeom['extenthalf']) + assert dgeom['extenthalf'].size == 2 + if dgeom['rcurve'] is not None: + dgeom['rcurve'] = float(dgeom['rcurve']) + if dgeom['nout'] is not None: + dgeom['nout'] = np.atleast_1d(dgeom['nout']) + dgeom['nout'] = dgeom['nout'] / np.linalg.norm(dgeom['nout']) + assert dgeom['nout'].size == 3 + if dgeom['nin'] is not None: + dgeom['nin'] = np.atleast_1d(dgeom['nin']) + dgeom['nin'] = dgeom['nin'] / np.linalg.norm(dgeom['nin']) + assert dgeom['nin'].size == 3 + if dgeom['e1'] is not None: + dgeom['e1'] = np.atleast_1d(dgeom['e1']) + dgeom['e1'] = dgeom['e1'] / np.linalg.norm(dgeom['e1']) + assert dgeom['e1'].size == 3 + assert dgeom['e2'] is not None + if dgeom['e2'] is not None: + dgeom['e2'] = np.atleast_1d(dgeom['e2']) + dgeom['e2'] = dgeom['e2'] / np.linalg.norm(dgeom['e2']) + assert dgeom['e2'].size == 3 + if dgeom['e1'] is not None: + assert np.abs(np.sum(dgeom['e1']*dgeom['e2'])) < 1.e-12 + if dgeom['nout'] is not None: + assert np.abs(np.sum(dgeom['e1']*dgeom['nout'])) < 1.e-12 + assert np.abs(np.sum(dgeom['e2']*dgeom['nout'])) < 1.e-12 + assert np.linalg.norm(np.cross(dgeom['e1'], dgeom['e2']) + - dgeom['nout']) < 1.e-12 + if dgeom['nin'] is not None: + assert np.abs(np.sum(dgeom['e1']*dgeom['nin'])) < 1.e-12 + assert np.abs(np.sum(dgeom['e2']*dgeom['nin'])) < 1.e-12 + assert np.linalg.norm(np.cross(dgeom['e1'], dgeom['e2']) + + dgeom['nin']) < 1.e-12 + if dgeom['rotateaxis'] is not None: + rotax = np.asarray(dgeom['rotateaxis'], dtype=float) + assert rotax.shape == (2, 3) + rotax[1,:] = rotax[1,:] / np.linalg.norm(rotax[1,:]) + dgeom['rotateaxis'] = rotax + if dgeom['rotateangle'] is not None: + dgeom['rotateangle'] = float(dgeom['rotateangle']) + if dgeom['mobile'] is None: + if dgeom['rotateaxis'] is not None: + dgeom['mobile'] = 'rotate' + else: + dgeom['mobile'] = False + else: + assert isinstance(dgeom['mobile'], str) + assert dgeom['mobile'] in cls._DEFLMOVEOK + return dgeom + + @classmethod + def _checkformat_dmat(cls, dmat=None): + if dmat is None: + return + assert isinstance(dmat, dict) + lkok = cls._get_keys_dmat() + assert all([isinstance(ss, str) for ss in dmat.keys()]) + assert all([ss in lkok for ss in dmat.keys()]) + for kk in cls._ddef['dmat'].keys(): + dmat[kk] = dmat.get(kk, cls._ddef['dmat'][kk]) + if dmat['d'] is not None: + dmat['d'] = float(dmat['d']) + if dmat['formula'] is not None: + assert isinstance(dmat['formula'], str) + if dmat['density'] is not None: + dmat['density'] = float(dmat['density']) + if dmat['lengths'] is not None: + dmat['lengths'] = np.atleast_1d(dmat['lengths']).ravel() + if dmat['angles'] is not None: + dmat['angles'] = np.atleast_1d(dmat['angles']).ravel() + if dmat['cut'] is not None: + dmat['cut'] = np.atleast_1d(dmat['cut']).ravel().astype(int) + assert dmat['cut'].size <= 4 + return dmat + + @classmethod + def _checkformat_dbragg(cls, dbragg=None): + if dbragg is None: + return + assert isinstance(dbragg, dict) + lkok = ['angle'] + assert all([isinstance(ss, str) for ss in dbragg.keys()]) + assert all([ss in lkok for ss in dbragg.keys()]) + + for kk in cls._ddef['dbragg'].keys(): + dbragg[kk] = dbragg.get(kk, cls._ddef['dbragg'][kk]) + return dbragg + + @classmethod + def _checkformat_inputs_dmisc(cls, color=None): + if color is None: + color = mpl.colors.to_rgba(cls._ddef['dmisc']['color']) + assert mpl.colors.is_color_like(color) + return tuple(mpl.colors.to_rgba(color)) + + ########### + # Get keys of dictionnaries + ########### + + @staticmethod + def _get_keys_dgeom(): + lk = ['Type', 'Typeoutline', + 'summit', 'center', 'extenthalf', 'surface', + 'nin', 'nout', 'e1', 'e2', 'rcurve', + 'mobile', 'rotateaxis', 'rotateangle'] + return lk + + @staticmethod + def _get_keys_dmat(): + lk = ['formula', 'density', 'symmetry', + 'lengths', 'angles', 'cut', 'd'] + return lk + + @staticmethod + def _get_keys_dbragg(): + lk = ['rockingcurve'] + return lk + + @staticmethod + def _get_keys_dmisc(): + lk = ['color'] + return lk + + ########### + # _init + ########### + + def _init(self, dgeom=None, dmat=None, dbragg=None, + color=None, **kwdargs): + allkwds = dict(locals(), **kwdargs) + largs = self._get_largs_dgeom() + kwds = self._extract_kwdargs(allkwds, largs) + self.set_dgeom(**kwds) + largs = self._get_largs_dmat() + kwds = self._extract_kwdargs(allkwds, largs) + self.set_dmat(**kwds) + largs = self._get_largs_dbragg() + kwds = self._extract_kwdargs(allkwds, largs) + self._set_dbragg(**kwds) + largs = self._get_largs_dmisc() + kwds = self._extract_kwdargs(allkwds, largs) + self._set_dmisc(**kwds) + self._dstrip['strip'] = 0 + + ########### + # set dictionaries + ########### + + def set_dgeom(self, dgeom=None): + dgeom = self._checkformat_dgeom(dgeom) + if dgeom['e1'] is not None: + if dgeom['nout'] is None: + dgeom['nout'] = np.cross(dgeom['e1'], dgeom['e2']) + if dgeom['nin'] is None: + dgeom['nin'] = -dgeom['nout'] + if dgeom['center'] is None: + dgeom['center'] = dgeom['summit'] + dgeom['nin']*dgeom['rcurve'] + if dgeom['summit'] is None: + dgeom['summit'] = dgeom['center'] + dgeom['nout']*dgeom['rcurve'] + + if dgeom['extenthalf'] is not None: + if dgeom['Type'] == 'sph' and dgeom['Typeoutline'] == 'rect': + ind = np.argmax(dgeom['extenthalf']) + dphi = dgeom['extenthalf'][ind] + sindtheta = np.sin(dgeom['extenthalf'][ind-1]) + dgeom['surface'] = 4.*dgeom['rcurve']**2*dphi*sindtheta + self._dgeom = dgeom + + def set_dmat(self, dmat=None): + dmat = self._checkformat_dmat(dmat) + self._dmat = dmat + + def _set_dbragg(self, dbragg=None): + dbragg = self._checkformat_dbragg(dbragg) + self._dbragg = dbragg + + def _set_color(self, color=None): + color = self._checkformat_inputs_dmisc(color=color) + self._dmisc['color'] = color + self._dplot['cross']['dP']['color'] = color + self._dplot['hor']['dP']['color'] = color + # self._dplot['3d']['dP']['color'] = color + + def _set_dmisc(self, color=None): + self._set_color(color) + + ########### + # strip dictionaries + ########### + + def _strip_dgeom(self, lkeep=None): + lkeep = self._get_keys_dgeom() + utils.ToFuObject._strip_dict(self._dgeom, lkeep=lkeep) + + def _strip_dmat(self, lkeep=None): + lkeep = self._get_keys_dmat() + utils.ToFuObject._strip_dict(self._dmat, lkeep=lkeep) + + def _strip_dbragg(self, lkeep=None): + lkeep = self._get_keys_dbragg() + utils.ToFuObject._strip_dict(self._dbragg, lkeep=lkeep) + + def _strip_dmisc(self, lkeep=['color']): + utils.ToFuObject._strip_dict(self._dmisc, lkeep=lkeep) + + ########### + # rebuild dictionaries + ########### + + def _rebuild_dgeom(self, lkeep=None): + lkeep = self._get_keys_dgeom() + reset = utils.ToFuObject._test_Rebuild(self._dgeom, lkeep=lkeep) + if reset: + utils.ToFuObject._check_Fields4Rebuild(self._dgeom, + lkeep=lkeep, dname='dgeom') + self._set_dgeom(dgeom=self._dgeom) + + def _rebuild_dmat(self, lkeep=None): + lkeep = self._get_keys_dmat() + reset = utils.ToFuObject._test_Rebuild(self._dmat, lkeep=lkeep) + if reset: + utils.ToFuObject._check_Fields4Rebuild(self._dmat, + lkeep=lkeep, dname='dmat') + self.set_dmat(self._dmat) + + def _rebuild_dbragg(self, lkeep=None): + lkeep = self._get_keys_dbragg() + reset = utils.ToFuObject._test_Rebuild(self._dbragg, lkeep=lkeep) + if reset: + utils.ToFuObject._check_Fields4Rebuild(self._dbragg, + lkeep=lkeep, dname='dbragg') + self.set_dbragg(self._dbragg) + + def _rebuild_dmisc(self, lkeep=['color']): + reset = utils.ToFuObject._test_Rebuild(self._dmisc, lkeep=lkeep) + if reset: + utils.ToFuObject._check_Fields4Rebuild(self._dmisc, + lkeep=lkeep, dname='dmisc') + self._set_dmisc(color=self.dmisc['color']) + + ########### + # _strip and get/from dict + ########### + + @classmethod + def _strip_init(cls): + cls._dstrip['allowed'] = [0,1] + nMax = max(cls._dstrip['allowed']) + doc = """ + 1: Remove nothing""" + doc = utils.ToFuObjectBase.strip.__doc__.format(doc,nMax) + if sys.version[0]=='2': + cls.strip.__func__.__doc__ = doc + else: + cls.strip.__doc__ = doc + + def strip(self, strip=0): + # super() + super(CrystalBragg, self).strip(strip=strip) + + def _strip(self, strip=0): + if strip==0: + self._rebuild_dgeom() + self._rebuild_dmat() + self._rebuild_dbragg() + self._rebuild_dmisc() + else: + self._strip_dgeom() + self._strip_dmat() + self._strip_dbragg() + self._strip_dmisc() + + def _to_dict(self): + dout = {'dgeom':{'dict':self._dgeom, 'lexcept':None}, + 'dmat':{'dict':self._dmat, 'lexcept':None}, + 'dbragg':{'dict':self._dbragg, 'lexcept':None}, + 'dmisc':{'dict':self._dmisc, 'lexcept':None}, + 'dplot':{'dict':self._dplot, 'lexcept':None}} + return dout + + def _from_dict(self, fd): + self._dgeom.update(**fd.get('dgeom', {})) + self._dmat.update(**fd.get('dmat', {})) + self._dbragg.update(**fd.get('dbragg', {})) + self._dmisc.update(**fd.get('dmisc', {})) + self._dplot.update(**fd.get('dplot', {})) + + # ----------- + # Properties + # ----------- + + @property + def Type(self): + """Return the type of structure """ + return self._Id.Type + + @property + def dgeom(self): + return self._dgeom + + @property + def dmat(self): + """Return the polygon defining the structure cross-section""" + return self._dmat + + @property + def dbragg(self): + """Return the polygon defining the structure cross-section""" + return self._dbragg + + @property + def dmisc(self): + return self._dmisc + + @property + def nin(self): + return self._dgeom['nin'] + + @property + def nout(self): + return self._dgeom['nout'] + + @property + def e1(self): + return self._dgeom['e1'] + + @property + def e2(self): + return self._dgeom['e2'] + + @property + def summit(self): + return self._dgeom['summit'] + + @property + def center(self): + return self._dgeom['center'] + + # ----------------- + # methods for color + # ----------------- + + def set_color(self, col): + self._set_color(col) + + def get_color(self): + return self._dmisc['color'] + + # ----------------- + # methods for printing + # ----------------- + + def get_summary(self, sep=' ', line='-', just='l', + table_sep=None, verb=True, return_=False): + """ Summary description of the object content """ + + # ----------------------- + # Build material + col0 = ['formula', 'symmetry', 'cut', 'density', + 'd (A)', 'bragg(%s A) (deg)'%str(self._DEFLAMB)] + ar0 = [self._dmat['formula'], self._dmat['symmetry'], + str(self._dmat['cut']), str(self._dmat['density']), + '{0:5.3f}'.format(self._dmat['d']*1.e10), + str(self.get_bragg_from_lamb(self._DEFLAMB)[0]*180./np.pi)] + + # ----------------------- + # Build geometry + col1 = ['Type', 'Type outline', 'surface (cm^2)', 'rcurve', + 'half-extent', 'summit', 'center', 'nin', 'e1', 'e2'] + ar1 = [self._dgeom['Type'], self._dgeom['Typeoutline'], + '{0:5.1f}'.format(self._dgeom['surface']*1.e4), + '{0:5.2f}'.format(self._dgeom['rcurve']), + str(np.round(self._dgeom['extenthalf'], decimals=3)), + str(np.round(self._dgeom['summit'], decimals=2)), + str(np.round(self._dgeom['center'], decimals=2)), + str(np.round(self._dgeom['nin'], decimals=2)), + str(np.round(self._dgeom['e1'], decimals=2)), + str(np.round(self._dgeom['e2'], decimals=2))] + + if self._dmisc.get('color') is not None: + col1.append('color') + ar1.append(str(self._dmisc['color'])) + + lcol = [col0, col1] + lar = [ar0, ar1] + + # ----------------------- + # Build mobile + if self._dgeom['mobile'] != False: + if self._dgeom['mobile'] == 'rotate': + col2 = ['Mov. type', 'axis pt.', 'axis vector', 'pos. (deg)'] + ar2 = [self._dgeom['mobile'], + str(np.round(self._dgeom['rotateaxis'][0], decimals=2)), + str(np.round(self._dgeom['rotateaxis'][1], decimals=2)), + '{0:7.3f}'.format(self._dgeom['rotateangle']*180./np.pi)] + lcol.append(col2) + lar.append(ar2) + return self._get_summary(lar, lcol, + sep=sep, line=line, table_sep=table_sep, + verb=verb, return_=return_) + # ----------------- + # methods for moving + # ----------------- + + @staticmethod + def _rotate_vector(vect, dangle, u, u1, u2): + c1 = np.sum(vect*u1) + c2 = np.sum(vect*u2) + return (np.sum(vect*u)*u + + (c1*np.cos(dangle) - c2*np.sin(dangle))*u1 + + (c2*np.cos(dangle) + c1*np.sin(dangle))*u2) + + + def _rotate(self, angle=None, dangle=None): + assert 'rotateaxis' in self._dgeom.keys() + lc = [angle is not None, dangle is not None] + assert np.sum(lc) == 1 + + if lc[0]: + assert 'rotateangle' in self._dgeom.keys() + dangle = (angle - self._dgeom['rotateangle'])%(2*np.pi) + dangle = np.arctan2(np.sin(dangle), np.cos(dangle)) + angle = (self._dgeom['rotateangle'] + dangle)%(2.*np.pi) + angle = np.arctan2(np.sin(angle), np.cos(angle)) + + # Define local frame (u, u1, u2) + OS = self._dgeom['summit'] - self._dgeom['rotateaxis'][0] + u = self._dgeom['rotateaxis'][1] + u = u / np.linalg.norm(u) + Z = np.sum(OS*u) + u1 = OS - Z*u + u1 = u1 / np.linalg.norm(u1) + u2 = np.cross(u, u1) + assert np.abs(np.linalg.norm(u2) - 1.) < 1.e-9 + + # Deduce constant distance from axis + dist = np.sum(OS*u1) + + summit = (self._dgeom['rotateaxis'][0] + + dist*(np.cos(dangle)*u1 + np.sin(dangle)*u2) + Z*u) + nout = self._rotate_vector(self._dgeom['nout'], dangle, u, u1, u2) + e1 = self._rotate_vector(self._dgeom['e1'], dangle, u, u1, u2) + e2 = self._rotate_vector(self._dgeom['e2'], dangle, u, u1, u2) + center = summit - self._dgeom['rcurve']*nout + assert np.abs(np.sum(nout*e1)) < 1.e-12 + assert np.abs(np.sum(nout*e2)) < 1.e-12 + assert np.abs(np.sum(e1*e2)) < 1.e-12 + self._dgeom.update({'summit': summit, 'center':center, + 'nin': -nout, 'nout':nout, 'e1': e1, 'e2': e2, + 'rotateangle': angle}) + + def move(self, kind=None, **kwdargs): + if self._dgeom['mobile'] == False: + return + if kind is None: + kind = self._dgeom['mobile'] + assert kind in self._DEFLMOVEOK + if kind == 'rotate': + self._rotate(**kwdargs) + + # ----------------- + # methods for surface and contour sampling + # ----------------- + + def sample_outline_plot(self, res=None): + if self._dgeom['Type'] == 'sph': + if self._dgeom['Typeoutline'] == 'rect': + func = _comp_optics.CrystBragg_sample_outline_plot_sphrect + outline = func(self._dgeom['center'], self.dgeom['nout'], + self._dgeom['e1'], self._dgeom['e2'], + self._dgeom['rcurve'], self._dgeom['extenthalf'], + res) + else: + raise NotImplementedError + else: + raise NotImplementedError + return outline + + def sample_outline_Rays(self, res=None): + if self._dgeom['Type'] == 'sph': + if self._dgeom['Typeoutline'] == 'rect': + func = _comp_optics.CrystBragg_sample_outline_Rays + pts, phi, theta = func() + else: + raise NotImplementedError + else: + raise NotImplementedError + return outline + + + + # ----------------- + # methods for surface and contour sampling + # ----------------- + + def _checkformat_bragglamb(self, bragg=None, lamb=None, n=None): + lc = [lamb is not None, bragg is not None] + assert np.sum(lc) == 1, "Provide lamb xor bragg!" + if lc[0]: + bragg = self.get_bragg_from_lamb(np.atleast_1d(lamb), + n=n) + else: + bragg = np.atleast_1d(bragg) + return bragg + + def _checkformat_get_Rays_from(self, phi=None, bragg=None): + assert phi is not None + assert bragg is not None + bragg = np.atleast_1d(bragg) + phi = np.atleast_1d(phi) + nrays = max(phi.size, bragg.size) + if not phi.shape == bragg.shape: + if phi.size == 1: + phi = np.full(bragg.shape, phi[0]) + elif bragg.size == 1: + bragg = np.full(phi.shape, bragg[0]) + else: + msg = "phi and bragg/lamb must have the same shape!\n" + msg += " phi.shape: %s\n"%str(phi.shape) + msg += " bragg/lamb.shape: %s\n"%str(bragg.shape) + raise Exception(msg) + return phi, bragg + + def get_Rays_from_summit(self, phi=None, bragg=None, + lamb=None, n=None, + returnas=object, config=None, name=None): + + # Check inputs + bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb) + phi, bragg = self._checkformat_get_Rays_from(phi=phi, bragg=bragg) + # assert phi.ndim == 1 + phi = phi[None, ...] + bragg = bragg[None, ...] + + # Prepare + shape = tuple([3] + [1 for ii in range(phi.ndim)]) + nin = self._dgeom['nin'].reshape(shape) + e1 = self._dgeom['e1'].reshape(shape) + e2 = self._dgeom['e2'].reshape(shape) + + # Compute start point (D) and unit vectors (us) + D = self._dgeom['summit'] + us = (np.sin(bragg)*nin + + np.cos(bragg)*(np.cos(phi)*e1 + np.sin(phi)*e2)) + + # Format output + if returnas == tuple: + return (D, us) + elif returnas == object: + from ._core import CamLOS1D + if name is None: + name = self.Id.Name + 'ExtractCam' + if us.ndim > 2: + us = us.reshape((3, phi.size)) + return CamLOS1D(dgeom=(D, us), Name=name, Diag=self.Id.Diag, + Exp=self.Id.Exp, shot=self.Id.shot, config=config) + + def get_Rays_envelop(self, + phi=None, bragg=None, lamb=None, n=None, + returnas=object, config=None, name=None): + # Check inputs + phi, bragg = self._checkformat_get_Rays_from(phi=phi, bragg=bragg, + lamb=lamb, n=n) + assert phi.ndim == 1 + + # Compute + func = _comp_optics.CrystBragg_sample_outline_Rays + D, us = func(self._dgeom['center'], self._dgeom['nout'], + self._dgeom['e1'], self._dgeom['e2'], + self._dgeom['rcurve'], self._dgeom['extenthalf'], + bragg, phi) + + # Format output + if returnas == tuple: + return (D, us) + elif returnas == object: + from ._core import CamLOS1D + if name is None: + name = self.Id.Name + 'ExtractCam' + return CamLOS1D(dgeom=(D, us), Name=name, Diag=self.Id.Diag, + Exp=self.Id.Exp, shot=self.Id.shot, config=config) + + + + # ----------------- + # methods for general plotting + # ----------------- + + def plot(self, lax=None, proj=None, res=None, element=None, + color=None, + dP=None, dI=None, dBs=None, dBv=None, + dVect=None, dIHor=None, dBsHor=None, + dBvHor=None, dleg=None, + draw=True, fs=None, wintit=None, Test=True): + kwdargs = locals() + lout = ['self'] + for k in lout: + del kwdargs[k] + return _plot_optics.CrystalBragg_plot(self, **kwdargs) + + # ----------------- + # methods for generic first-approx + # ----------------- + + def get_bragg_from_lamb(self, lamb, n=None): + """ Braggs' law: n*lamb = 2dsin(bragg) """ + if self._dmat['d'] is None: + msg = "Interplane distance d no set !\n" + msg += " => self.set_dmat({'d':...})" + raise Exception(msg) + return _comp_optics.get_bragg_from_lamb(np.atleast_1d(lamb), + self._dmat['d'], n=n) + + def get_lamb_from_bragg(self, bragg, n=None): + """ Braggs' law: n*lamb = 2dsin(bragg) """ + if self._dmat['d'] is None: + msg = "Interplane distance d no set !\n" + msg += " => self.set_dmat({'d':...})" + raise Exception(msg) + return _comp_optics.get_lamb_from_bragg(np.atleast_1d(bragg), + self._dmat['d'], n=n) + + def get_approx_detector_frame(self, bragg=None, lamb=None, + rcurve=None, n=None, plot=False): + """ See notes for details on notations """ + + # Check / format inputs + if rcurve is None: + rcurve = self._dgeom['rcurve'] + bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) + + # Compute crystal-centered parameters in (nout, e1, e2) + func = _comp_optics.get_approx_detector_rel + det_dist, det_nout_rel, det_ei_rel = func(rcurve, bragg) + + # Deduce absolute position in (x, y, z) + func = _comp_optics.get_det_abs_from_rel + det_cent, det_nout, det_ei, det_ej = func(det_dist, + det_nout_rel, det_ei_rel, + self._dgeom['summit'], + self._dgeom['nout'], + self._dgeom['e1'], + self._dgeom['e2']) + + if plot: + dax = self.plot() + p0 = np.repeat(det_cent[:,None], 3, axis=1) + vv = np.vstack((det_nout, det_ei, det_ej)).T + dax['cross'].plot(np.hypot(det_cent[0], det_cent[1]), + det_cent[2], 'xb') + dax['hor'].plot(det_cent[0], det_cent[1], 'xb') + dax['cross'].quiver(np.hypot(p0[0, :], p0[1, :]), p0[2, :], + np.hypot(vv[0, :], vv[1, :]), vv[2, :], + units='xy', color='b') + dax['hor'].quiver(p0[0, :], p0[1, :], vv[0, :], vv[1, :], + units='xy', color='b') + return det_cent, det_nout, det_ei, det_ej + + def get_local_noute1e2(self, theta, psi): + if np.allclose([theta, psi], [np.pi/2., 0.]): + summit = self._dgeom['summit'] + nout = self._dgeom['nout'] + e1, e2 = self._dgeom['e1'], self._dgeom['e2'] + else: + func = _comp_optics.CrystBragg_get_noute1e2_from_psitheta + nout, e1, e2 = func(self._dgeom['nout'], + self._dgeom['e1'], self._dgeom['e2'], + [psi], [theta]) + nout, e1, e2 = nout[:, 0], e1[:, 0], e2[:, 0] + summit = self._dgeom['center'] + self._dgeom['rcurve']*nout + return summit, nout, e1, e2 + + + def calc_xixj_from_phibragg(self, phi=None, + bragg=None, lamb=None, n=None, + theta=None, psi=None, + det_cent=None, det_nout=None, + det_ei=None, det_ej=None, + data=None, plot=True, ax=None): + """ Assuming crystal's summit as frame origin + + According to [1], this assumes a local frame centered on the crystal + + These calculations are independent from the tokamak's frame: + The origin of the local frame is the crystal's summit + The (O, ez) axis is the crystal's normal + The crystal is tangent to (O, ex, ey) + + [1] tofu/Notes_Upgrades/SpectroX2D/SpectroX2D_EllipsesOnPlane.pdf + + Parameters: + ----------- + Z: float + Detector's plane intersection with (O, ez) axis + n: np.ndarray + (3,) array containing local (x,y,z) coordinates of the plane's + normal vector + """ + # Check / format inputs + bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) + phi = np.atleast_1d(phi) + assert bragg.ndim == phi.ndim == 1 + nphi, nbragg = phi.size, bragg.size + npts = max(nphi, nbragg) + assert nbragg in [1, npts] and nphi in [1, npts], (nbragg, nphi) + if nbragg < npts: + bragg = np.full((npts,), bragg[0]) + elif nphi < npts: + phi = np.full((npts,), phi[0]) + + lc = [det_cent is None, det_nout is None, + det_ei is None, det_ej is None] + assert all(lc) or not any(lc) + if all(lc): + func = self.get_approx_detector_frame + det_cent, det_nout, det_ei, det_ej = func(lamb=self._DEFLAMB) + + # Get local summit nout, e1, e2 if non-centered + if theta is None: + theta = np.pi/2. + if psi is None: + psi = 0. + summit, nout, e1, e2 = self.get_local_noute1e2(theta, psi) + + # Compute + xi, xj = _comp_optics.calc_xixj_from_braggphi(summit, + det_cent, det_nout, + det_ei, det_ej, + nout, e1, e2, + bragg, phi) + if plot: + func = _plot_optics.CrystalBragg_plot_approx_detector_params + ax = func(bragg, xi, xj, data, ax) + return xi, xj + + @staticmethod + def _checkformat_xixj(xi, xj): + xi = np.atleast_1d(xi) + xj = np.atleast_1d(xj) + if xi.shape == xj.shape: + return xi, xj, (xi, xj) + else: + return xi, xj, np.meshgrid(xi, xj) + + def calc_phibragg_from_xixj(self, xi, xj, n=None, + det_cent=None, det_ei=None, det_ej=None, + theta=None, psi=None, + plot=True, ax=None, **kwdargs): + + # Check / format inputs + xi, xj, (xii, xjj) = self._checkformat_xixj(xi, xj) + + lc = [det_cent is None, det_ei is None, det_ej is None] + assert all(lc) or not any(lc) + if all(lc): + func = self.get_approx_detector_frame + det_cent, _, det_ei, det_ej = func(lamb=self._DEFLAMB) + + # Get local summit nout, e1, e2 if non-centered + if theta is None: + theta = np.pi/2. + if psi is None: + psi = 0. + summit, nout, e1, e2 = self.get_local_noute1e2(theta, psi) + + # Compute + func = _comp_optics.calc_braggphi_from_xixj + bragg, phi = func(xii, xjj, det_cent, det_ei, det_ej, + summit, -nout, e1, e2) + + if plot != False: + func = _plot_optics.CrystalBragg_plot_braggangle_from_xixj + lax = func(xi=xii, xj=xjj, + ax=ax, plot=plot, + bragg=bragg * 180./np.pi, + angle=phi * 180./np.pi, + braggunits='deg', angunits='deg', **kwdargs) + return bragg, phi + + def plot_johannerror(self, lamb=None): + raise NotImplementedError + + def plot_data_vs_lambphi(self, xi=None, xj=None, data=None, + det_cent=None, det_ei=None, det_ej=None, + theta=None, psi=None, n=None, + plot=True, fs=None, + cmap=None, vmin=None, vmax=None): + # Check / format inputs + assert data is not None + xi, xj, (xii, xjj) = self._checkformat_xixj(xi, xj) + nxi = xi.size if xi is not None else np.unique(xii).size + nxj = xj.size if xj is not None else np.unique(xjj).size + + # Compute lamb / phi + func = self.calc_phibragg_from_xixj + bragg, phi = func(xii, xjj, n=n, + det_cent=det_cent, det_ei=det_ei, det_ej=det_ej, + theta=theta, psi=psi, plot=False) + assert bragg.shape == phi.shape == data.shape + lamb = self.get_lamb_from_bragg(bragg, n=n) + + # Compute lambfit / phifit and spectrum1d + lambfit, phifit = _comp_optics.get_lambphifit(lamb, phi, nxi, nxj) + lambfitbins = 0.5*(lambfit[1:] + lambfit[:-1]) + ind = np.digitize(lamb, lambfitbins) + spect1d = np.array([np.nanmean(data[ind==ii]) for ii in np.unique(ind)]) + + # plot + func = _plot_optics.CrystalBragg_plot_data_vs_lambphi + ax = func(xi, xj, bragg, lamb, phi, data, + lambfit=lambfit, phifit=phifit, spect1d=spect1d, + cmap=cmap, vmin=vmin, vmax=vmax, fs=fs) + return ax + + def plot_data_fit2d(self, xi=None, xj=None, data=None, mask=None, + det_cent=None, det_ei=None, det_ej=None, + theta=None, psi=None, n=None, + nlamb=None, lamb0=None, forcelamb=False, + deg=None, knots=None, nbsplines=None, + method=None, max_nfev=None, + xtol=None, ftol=None, gtol=None, + loss=None, verbose=0, debug=None, + plot=True, fs=None, dlines=None, dmoments=None, + cmap=None, vmin=None, vmax=None): + # Check / format inputs + assert data is not None + xi, xj, (xii, xjj) = self._checkformat_xixj(xi, xj) + nxi = xi.size if xi is not None else np.unique(xii).size + nxj = xj.size if xj is not None else np.unique(xjj).size + + # Compute lamb / phi + func = self.calc_phibragg_from_xixj + bragg, phi = func(xii, xjj, n=n, + det_cent=det_cent, det_ei=det_ei, det_ej=det_ej, + theta=theta, psi=psi, plot=False) + assert bragg.shape == phi.shape == data.shape + lamb = self.get_lamb_from_bragg(bragg, n=n) + + # Compute lambfit / phifit and spectrum1d + lambfit, phifit = _comp_optics.get_lambphifit(lamb, phi, nxi, nxj) + lambfitbins = 0.5*(lambfit[1:] + lambfit[:-1]) + ind = np.digitize(lamb, lambfitbins) + spect1d = np.array([np.nanmean(data[ind==ii]) for ii in np.unique(ind)]) + + # Compute fit for spect1d to get lamb0 if not provided + import tofu.data._spectrafit2d as _spectrafit2d + + func = _spectrafit2d.multiplegaussianfit1d + dfit1d = func(lambfit, spect1d, + nmax=nlamb, lamb0=lamb0, forcelamb=forcelamb, + p0=None, bounds=None, + max_nfev=None, xtol=xtol, verbose=0, + percent=20, plot_debug=False) + + # Reorder wrt lamb0 + ind = np.argsort(dfit1d['lamb0']) + for kk in ['lamb0', 'amp', 'ampstd', + 'sigma', 'sigmastd', 'dlamb', 'dlambstd']: + if dfit1d[kk].ndim == 1: + dfit1d[kk] = dfit1d[kk][ind] + else: + dfit1d[kk] = dfit1d[kk][0,ind] + + + # Compute dfit2d + if mask is None: + mask = np.ones(data.shape, dtype=bool) + func = _spectrafit2d.multigaussianfit2d + dfit2d = func(lamb[mask].ravel(), phi[mask].ravel(), data[mask].ravel(), + std=dfit1d['std'], + lamb0=dfit1d['lamb0'], forcelamb=forcelamb, + deg=deg, nbsplines=nbsplines, + method=method, max_nfev=max_nfev, + xtol=xtol, ftol=ftol, gtol=gtol, + loss=loss, verbose=verbose, debug=debug) + + + # plot + func = _plot_optics.CrystalBragg_plot_data_vs_fit + ax = func(xi, xj, bragg, lamb, phi, data, mask=mask, + lambfit=lambfit, phifit=phifit, spect1d=spect1d, + dfit1d=dfit1d, dfit2d=dfit2d, lambfitbins=lambfitbins, + cmap=cmap, vmin=vmin, vmax=vmax, + fs=fs, dmoments=dmoments) + return ax, dfit1d, None diff --git a/tofu/geom/_plot_optics.py b/tofu/geom/_plot_optics.py new file mode 100644 index 000000000..07e669218 --- /dev/null +++ b/tofu/geom/_plot_optics.py @@ -0,0 +1,527 @@ + + +import numpy as np +from scipy.interpolate import BSpline +import matplotlib.pyplot as plt +from matplotlib.colors import ListedColormap +import matplotlib.gridspec as gridspec +from matplotlib.axes._axes import Axes +from mpl_toolkits.mplot3d import Axes3D + +# tofu +from tofu.version import __version__ +from . import _def as _def + +_GITHUB = 'https://github.com/ToFuProject/tofu/issues' +_WINTIT = 'tofu-%s report issues / requests at %s'%(__version__, _GITHUB) + +_QUIVERCOLOR = plt.cm.viridis(np.linspace(0, 1, 3)) +_QUIVERCOLOR = np.array([[1., 0., 0., 1.], + [0., 1., 0., 1.], + [0., 0., 1., 1.]]) +_QUIVERCOLOR = ListedColormap(_QUIVERCOLOR) + + +# Generic +def _check_projdax_mpl(dax=None, proj=None, fs=None, wintit=None): + + # ---------------------- + # Check proj + if proj is None: + proj = 'all' + assert isinstance(proj, str) + proj = proj.lower() + lproj = ['cross', 'hor', '3d'] + assert proj in lproj + ['all'] + if proj == 'all': + proj = ['cross', 'hor'] + else: + proj = [proj] + + # ---------------------- + # Check dax + lc = [dax is None, issubclass(dax.__class__, Axes), + isinstance(dax, dict), isinstance(dax, list)] + assert any(lc) + if lc[0]: + dax = dict.fromkeys(proj) + elif lc[1]: + assert len(proj) == 1 + dax = {proj[0]: dax} + elif lc[2]: + lcax = [(kk in proj + and (ax is None or issubclass(ax.__class__, Axes))) + for kk, ax in dax.items()] + if not all(lcax): + msg = "Wrong key or axes in dax:\n" + msg += " - proj = %s"%str(proj) + msg += " - dax = %s"%str(dax) + raise Exception(msg) + else: + assert len(dax) == 2 + assert all([ax is None or issubclass(ax.__class__, Axes) + for ax in dax]) + dax = {'cross': dax[0], 'hor': dax[1]} + + # Populate with default axes if necessary + if 'cross' in proj and 'hor' in proj: + if 'cross' in proj and 'hor' in proj: + if dax['cross'] is None: + assert dax['hor'] is None + lax = _def.Plot_LOSProj_DefAxes('all', fs=fs, wintit=wintit) + dax['cross'], dax['hor'] = lax + elif 'cross' in proj and dax['cross'] is None: + dax['cross'] = _def.Plot_LOSProj_DefAxes('cross', fs=fs, + wintit=wintit) + elif 'hor' in proj and dax['hor'] is None: + dax['hor'] = _def.Plot_LOSProj_DefAxes('hor', fs=fs, + wintit=wintit) + elif '3d' in proj and dax['3d'] is None: + dax['3d'] = _def.Plot_3D_plt_Tor_DefAxes(fs=fs, + wintit=wintit) + for kk in lproj: + dax[kk] = dax.get(kk, None) + return dax + + + +# ################################################################# +# ################################################################# +# Generic geometry plot +# ################################################################# +# ################################################################# + +def CrystalBragg_plot(cryst, lax=None, proj=None, res=None, element=None, + color=None, dP=None, + dI=None, dBs=None, dBv=None, + dVect=None, dIHor=None, dBsHor=None, dBvHor=None, + dleg=None, indices=False, + draw=True, fs=None, wintit=None, tit=None, Test=True): + + # --------------------- + # Check / format inputs + + if Test: + msg = "Arg proj must be in ['cross','hor','all','3d'] !" + assert type(draw) is bool, "Arg draw must be a bool !" + assert cryst.__class__.__name__ == 'CrystalBragg' + if wintit is None: + wintit = _WINTIT + if dleg is None: + dleg = _def.TorLegd + + # --------------------- + # call plotting functions + + kwa = dict(fs=fs, wintit=wintit, Test=Test) + if proj == '3d': + # Temporary matplotlib issue + dleg = None + else: + dax = _CrystalBragg_plot_crosshor(cryst, proj=proj, res=res, dax=lax, element=element, + color=color) + + # recompute the ax.dataLim + ax0 = None + for kk, vv in dax.items(): + if vv is None: + continue + dax[kk].relim() + dax[kk].autoscale_view() + if dleg is not None: + dax[kk].legend(**dleg) + ax0 = vv + + # set title + if tit != False: + ax0.figure.suptitle(tit) + if draw: + ax0.figure.canvas.draw() + return dax + +def _CrystalBragg_plot_crosshor(cryst, proj=None, dax=None, element=None, res=None, + Pdict=_def.TorPd, Idict=_def.TorId, Bsdict=_def.TorBsd, + Bvdict=_def.TorBvd, Vdict=_def.TorVind, + color=None, ms=None, quiver_cmap=None, + LegDict=_def.TorLegd, indices=False, + draw=True, fs=None, wintit=None, Test=True): + if Test: + assert type(Pdict) is dict, 'Arg Pdict should be a dictionary !' + assert type(Idict) is dict, "Arg Idict should be a dictionary !" + assert type(Bsdict) is dict, "Arg Bsdict should be a dictionary !" + assert type(Bvdict) is dict, "Arg Bvdict should be a dictionary !" + assert type(Vdict) is dict, "Arg Vdict should be a dictionary !" + assert type(LegDict) is dict or LegDict is None, 'Arg LegDict should be a dictionary !' + + # --------------------- + # Check / format inputs + + if element is None: + element = 'oscvr' + element = element.lower() + if 'v' in element and quiver_cmap is None: + quiver_cmap = _QUIVERCOLOR + if color is None: + if cryst._dmisc.get('color') is not None: + color = cryst._dmisc['color'] + else: + color = 'k' + if ms is None: + ms = 6 + + # --------------------- + # Prepare axe and data + + dax = _check_projdax_mpl(dax=dax, proj=proj, fs=fs, wintit=wintit) + + if 's' in element or 'v' in element: + summ = cryst._dgeom['summit'] + if 'c' in element: + cent = cryst._dgeom['center'] + if 'r' in element: + ang = np.linspace(0, 2.*np.pi, 200) + rr = 0.5*cryst._dgeom['rcurve'] + row = cryst._dgeom['center'] + rr*cryst._dgeom['nout'] + row = (row[:, None] + + rr*(np.cos(ang)[None, :]*cryst._dgeom['nout'][:, None] + + np.sin(ang)[None, :]*cryst._dgeom['e1'][:, None])) + + # --------------------- + # plot + + if 'o' in element: + cont = cryst.sample_outline_plot(res=res) + if dax['cross'] is not None: + dax['cross'].plot(np.hypot(cont[0,:], cont[1,:]), cont[2,:], + ls='-', c=color, marker='None', + label=cryst.Id.NameLTX+' contour') + if dax['hor'] is not None: + dax['hor'].plot(cont[0,:], cont[1,:], + ls='-', c=color, marker='None', + label=cryst.Id.NameLTX+' contour') + if 's' in element: + if dax['cross'] is not None: + dax['cross'].plot(np.hypot(summ[0], summ[1]), summ[2], + marker='^', ms=ms, c=color, + label=cryst.Id.NameLTX+" summit") + if dax['hor'] is not None: + dax['hor'].plot(summ[0], summ[1], + marker='^', ms=ms, c=color, + label=cryst.Id.NameLTX+" summit") + if 'c' in element: + if dax['cross'] is not None: + dax['cross'].plot(np.hypot(cent[0], cent[1]), cent[2], + marker='o', ms=ms, c=color, + label=cryst.Id.NameLTX+" center") + if dax['hor'] is not None: + dax['hor'].plot(cent[0], cent[1], + marker='o', ms=ms, c=color, + label=cryst.Id.NameLTX+" center") + if 'r' in element: + if dax['cross'] is not None: + dax['cross'].plot(np.hypot(row[0,:], row[1,:]), row[2,:], + ls='--', color=color, marker='None', + label=cryst.Id.NameLTX+' rowland') + if dax['hor'] is not None: + dax['hor'].plot(row[0,:], row[1,:], + ls='--', color=color, marker='None', + label=cryst.Id.NameLTX+' rowland') + if 'v' in element: + nin = cryst._dgeom['nin'] + e1, e2 = cryst._dgeom['e1'], cryst._dgeom['e2'] + p0 = np.repeat(summ[:,None], 3, axis=1) + v = np.concatenate((nin[:, None], e1[:, None], e2[:, None]), axis=1) + if dax['cross'] is not None: + dax['cross'].quiver(np.hypot(p0[0,:], p0[1,:]), p0[2,:], + np.hypot(v[0,:], v[1,:]), v[2,:], + np.r_[0., 0.5, 1.], cmap=quiver_cmap, + angles='xy', scale_units='xy', + label=cryst.Id.NameLTX+" unit vect", **Vdict) + if dax['hor'] is not None: + dax['hor'].quiver(p0[0,:], p0[1,:], + v[0,:], v[1,:], + np.r_[0., 0.5, 1.], cmap=quiver_cmap, + angles='xy', scale_units='xy', + label=cryst.Id.NameLTX+" unit vect", **Vdict) + + return dax + + +# ################################################################# +# ################################################################# +# Bragg diffraction plot +# ################################################################# +# ################################################################# + +# Deprecated ? re-use ? +def CrystalBragg_plot_approx_detector_params(Rrow, bragg, d, Z, + frame_cent, nn): + + R = 2.*Rrow + L = 2.*R + ang = np.linspace(0., 2.*np.pi, 100) + + fig = plt.figure() + ax = fig.add_axes([0.1,0.1,0.8,0.8], aspect='equal') + + ax.axvline(0, ls='--', c='k') + ax.plot(Rrow*np.cos(ang), Rrow + Rrow*np.sin(ang), c='r') + ax.plot(R*np.cos(ang), R + R*np.sin(ang), c='b') + ax.plot(L*np.cos(bragg)*np.r_[-1,0,1], + L*np.sin(bragg)*np.r_[1,0,1], c='k') + ax.plot([0, d*np.cos(bragg)], [Rrow, d*np.sin(bragg)], c='r') + ax.plot([0, d*np.cos(bragg)], [Z, d*np.sin(bragg)], 'g') + ax.plot([0, L/10*nn[1]], [Z, Z+L/10*nn[2]], c='g') + ax.plot(frame_cent[1]*np.cos(2*bragg-np.pi), + Z + frame_cent[1]*np.sin(2*bragg-np.pi), c='k', marker='o', ms=10) + + ax.set_xlabel(r'y') + ax.set_ylabel(r'z') + ax.legend(loc='upper left', bbox_to_anchor=(1.05, 1.), frameon=False) + return ax + + +def CrystalBragg_plot_xixj_from_braggangle(bragg=None, xi=None, xj=None, + data=None, ax=None): + if ax is None: + fig = plt.figure() + ax = fig.add_axes([0.1,0.1,0.8,0.8], aspect='equal') + + for ii in range(len(bragg)): + deg ='{0:07.3f}'.format(bragg[ii]*180/np.pi) + ax.plot(xi[:,ii], xj[:,ii], '.', label='bragg %s'%deg) + + ax.set_xlabel(r'xi') + ax.set_ylabel(r'yi') + ax.legend(loc='upper left', bbox_to_anchor=(1.05, 1.), frameon=False) + return ax + + + + +def CrystalBragg_plot_braggangle_from_xixj(xi=None, xj=None, + bragg=None, angle=None, + ax=None, plot=None, + braggunits='rad', angunits='rad', + **kwdargs): + + if isinstance(plot, bool): + plot = 'contour' + + if ax is None: + fig = plt.figure() + ax0 = fig.add_axes([0.1, 0.1, 0.35, 0.8], aspect='equal') + ax1 = fig.add_axes([0.55, 0.1, 0.35, 0.8], aspect='equal') + ax = [ax0, ax1] + if plot == 'contour': + if 'levels' in kwdargs.keys(): + lvls = kwdargs['levels'] + del kwdargs['levels'] + obj0 = ax[0].contour(xi, xj, bragg, lvls, **kwdargs) + obj1 = ax[1].contour(xi, xj, angle, lvls, **kwdargs) + else: + obj0 = ax[0].contour(xi, xj, bragg, **kwdargs) + obj1 = ax[1].contour(xi, xj, angle, **kwdargs) + elif plot == 'imshow': + extent = (xi.min(), xi.max(), xj.min(), xj.max()) + obj0 = ax[0].imshow(bragg, extent=extent, aspect='equal', + adjustable='datalim', **kwdargs) + obj1 = ax[1].imshow(angle, extent=extent, aspect='equal', + adjustable='datalim', **kwdargs) + elif plot == 'pcolor': + obj0 = ax[0].pcolor(xi, xj, bragg, **kwdargs) + obj1 = ax[1].pcolor(xi, xj, angle, **kwdargs) + ax[0].set_xlabel(r'xi') + ax[1].set_xlabel(r'xi') + ax[0].set_ylabel(r'yi') + ax[1].set_ylabel(r'yi') + cax0 = plt.colorbar(obj0, ax=ax[0]) + cax1 = plt.colorbar(obj1, ax=ax[1]) + cax0.ax.set_ylabel(r'$\theta_{bragg}$ (%s)'%braggunits) + cax1.ax.set_ylabel(r'$ang$ (%s)'%angunits) + return ax + + +def CrystalBragg_plot_data_vs_lambphi(xi, xj, bragg, lamb, phi, data, + lambfit=None, phifit=None, spect1d=None, + cmap=None, vmin=None, vmax=None, + fs=None, dmargin=None, + angunits='deg'): + + # Check inputs + # ------------ + + if fs is None: + fs = (14,8) + if cmap is None: + cmap = plt.cm.viridis + if dmargin is None: + dmargin = {'left':0.03, 'right':0.99, + 'bottom':0.05, 'top':0.92, + 'wspace':None, 'hspace':0.4} + assert angunits in ['deg', 'rad'] + if angunits == 'deg': + bragg = bragg*180./np.pi + phi = phi*180./np.pi + phifit = phifit*180./np.pi + + + # pre-compute + # ------------ + + # extent + extent = (xi.min(), xi.max(), xj.min(), xj.max()) + extent2 = (lambfit.min(), lambfit.max(), phifit.min(), phifit.max()) + + # Plot + # ------------ + + fig = fig = plt.figure(figsize=fs) + gs = gridspec.GridSpec(4, 3, **dmargin) + ax0 = fig.add_subplot(gs[:3, 0], aspect='equal', adjustable='datalim') + ax1 = fig.add_subplot(gs[:3, 1], aspect='equal', adjustable='datalim', + sharex=ax0, sharey=ax0) + axs1 = fig.add_subplot(gs[3, 1], sharex=ax0) + ax2 = fig.add_subplot(gs[:3, 2]) + axs2 = fig.add_subplot(gs[3, 2], sharex=ax2, sharey=axs1) + + ax0.set_title('Coordinates transform') + ax1.set_title('Camera image') + ax2.set_title('Camera image transformed') + + ax0.set_ylabel(r'incidence angle ($deg$)') + + ax0.contour(xi, xj, bragg, 10, cmap=cmap) + ax0.contour(xi, xj, phi, 10, cmap=cmap, ls='--') + ax1.imshow(data, extent=extent, aspect='equal', + origin='lower', vmin=vmin, vmax=vmax) + axs1.plot(xi, np.sum(data, axis=0), c='k', ls='-') + ax2.scatter(lamb.ravel(), phi.ravel(), c=data.ravel(), s=2, + marker='s', edgecolors='None', + cmap=cmap, vmin=vmin, vmax=vmax) + axs2.plot(lambfit, spect1d, c='k', ls='-') + + ax2.set_xlim(extent2[0], extent2[1]) + ax2.set_ylim(extent2[2], extent2[3]) + + return [ax0, ax1] + + +def CrystalBragg_plot_data_vs_fit(xi, xj, bragg, lamb, phi, data, mask=None, + lambfit=None, phifit=None, spect1d=None, + dfit1d=None, dfit2d=None, lambfitbins=None, + cmap=None, vmin=None, vmax=None, + fs=None, dmargin=None, + angunits='deg', dmoments=None): + + # Check inputs + # ------------ + + if fs is None: + fs = (16, 9) + if cmap is None: + cmap = plt.cm.viridis + if dmargin is None: + dmargin = {'left':0.03, 'right':0.99, + 'bottom':0.05, 'top':0.92, + 'wspace':None, 'hspace':0.4} + assert angunits in ['deg', 'rad'] + if angunits == 'deg': + bragg = bragg*180./np.pi + phi = phi*180./np.pi + phifit = phifit*180./np.pi + + + # pre-compute + # ------------ + + # extent + extent = (xi.min(), xi.max(), xj.min(), xj.max()) + extent2 = (lambfit.min(), lambfit.max(), phifit.min(), phifit.max()) + + ind = np.digitize(lamb[mask].ravel(), lambfitbins) + spect2dmean = np.zeros((lambfitbins.size+1,)) + for ii in range(lambfitbins.size+1): + indi = ind==ii + if np.any(indi): + spect2dmean[ii] = np.nanmean(dfit2d['fit'][indi]) + + # Plot + # ------------ + + fig = plt.figure(figsize=fs) + gs = gridspec.GridSpec(4, 6, **dmargin) + ax0 = fig.add_subplot(gs[:3, 0], aspect='equal', adjustable='datalim') + ax1 = fig.add_subplot(gs[:3, 1], aspect='equal', adjustable='datalim', + sharex=ax0, sharey=ax0) + axs1 = fig.add_subplot(gs[3, 1], sharex=ax0) + ax2 = fig.add_subplot(gs[:3, 2]) + axs2 = fig.add_subplot(gs[3, 2], sharex=ax2, sharey=axs1) + ax3 = fig.add_subplot(gs[:3, 3], sharex=ax2, sharey=ax2) + axs3 = fig.add_subplot(gs[3, 3], sharex=ax2)#, sharey=axs1) + ax4 = fig.add_subplot(gs[:3, 4], sharex=ax2, sharey=ax2) + axs4 = fig.add_subplot(gs[3, 3], sharex=ax2)#, sharey=axs1) + ax5 = fig.add_subplot(gs[:3, 5], sharey=ax2) + + ax0.set_title('Coordinates transform') + ax1.set_title('Camera image') + ax2.set_title('Camera image transformed') + ax3.set_title('2d spectral fit') + ax4.set_title('2d error') + ax5.set_title('Moments') + + ax4.set_xlabel('%s'%angunits) + ax0.set_ylabel(r'incidence angle ($deg$)') + + ax0.contour(xi, xj, bragg, 10, cmap=cmap) + ax0.contour(xi, xj, phi, 10, cmap=cmap, ls='--') + ax1.imshow(data, extent=extent, aspect='equal', + origin='lower', vmin=vmin, vmax=vmax) + axs1.plot(xi, np.nanmean(data, axis=0), c='k', ls='-') + ax2.scatter(lamb.ravel(), phi.ravel(), c=data.ravel(), s=1, + marker='s', edgecolors='None', + cmap=cmap, vmin=vmin, vmax=vmax) + axs2.plot(lambfit, spect1d, c='k', ls='None', marker='.', ms=4) + axs2.plot(lambfit, dfit1d['fit'].ravel(), c='r', ls='-', label='fit') + for ll in dfit1d['lamb0']: + axs2.axvline(ll, c='k', ls='--') + + # dfit2d + ax3.scatter(lamb[mask].ravel(), phi[mask].ravel(), c=dfit2d['fit'], s=1, + marker='s', edgecolors='None', + cmap=cmap, vmin=vmin, vmax=vmax) + axs3.plot(lambfit, spect1d, c='k', ls='None', marker='.') + axs3.plot(lambfit, spect2dmean, c='b', ls='-') + err = dfit2d['fit'] - data[mask].ravel() + errmax = np.max(np.abs(err)) + ax4.scatter(lamb[mask].ravel(), phi[mask].ravel(), c=err, s=1, + marker='s', edgecolors='None', + cmap=plt.cm.seismic, vmin=-errmax, vmax=errmax) + + # Moments + if dmoments is not None: + if dmoments.get('ratio') is not None: + ind = dmoments['ratio'].get('ind') + if ind is None: + ind = [np.argmin(np.abs(dfit2d['lamb0']-ll)) + for ll in dmoments['ratio']['lamb']] + for indi in ind: + axs3.axvline(dfit2d['lamb0'][indi], c='k', ls='--') + amp0 = BSpline(dfit2d['knots'], dfit2d['camp'][ind[0],:], dfit2d['deg'])(phifit) + amp1 = BSpline(dfit2d['knots'], dfit2d['camp'][ind[1],:], dfit2d['deg'])(phifit) + lab = dmoments['ratio']['name'] + '{} / {}' + ratio = (amp0 / amp1) / np.nanmax(amp0 / amp1) + ax5.plot(amp0 / amp1, phifit, ls='-', c='k', label=lab) + if dmoments.get('sigma') is not None: + ind = dmoments['sigma'].get('ind') + if ind is None: + ind = np.argmin(np.abs(dfit2d['lamb0']-dmoments['sigma']['lamb'])) + axs3.axvline(dfit2d['lamb0'][ind], c='b', ls='--') + sigma = BSpline(dfit2d['knots'], dfit2d['csigma'][ind,:], dfit2d['deg'])(phifit) + lab = r'$\sigma({} A)$'.format(np.round(dfit2d['lamb0'][ind]*1.e10), + 4) + ax5.plot(sigma/np.nanmax(sigma), phifit, ls='-', c='b', label=lab) + + ax2.set_xlim(extent2[0], extent2[1]) + ax2.set_ylim(extent2[2], extent2[3]) + return [ax0, ax1] diff --git a/tofu/geom/inputs/test_data/SpectroX2D_WEST_55562_t39.145s.npz b/tofu/geom/inputs/test_data/SpectroX2D_WEST_55562_t39.145s.npz new file mode 100644 index 000000000..c2b8d33bd Binary files /dev/null and b/tofu/geom/inputs/test_data/SpectroX2D_WEST_55562_t39.145s.npz differ diff --git a/tofu/version.py b/tofu/version.py index fad71f20a..99db56810 100644 --- a/tofu/version.py +++ b/tofu/version.py @@ -1,2 +1,2 @@ # Do not edit, pipeline versioning governed by git tags! -__version__ = '1.4.1-345-ga146ee8' +__version__ = '1.4.1-396-g806dbc9'