Issue202 spectro x2 d crystal

Notes_Upgrades/SpectroX2D/SpectroX2D_ConeCircleMagAxis.tex

@@ -0,0 +1,133 @@
+\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}
+\newcommand{\lra}{\Leftrightarrow}
+\DeclareMathOperator{\n}{\underline{n}}
+\DeclareMathOperator{\ex}{\underline{e}_x}
+\DeclareMathOperator{\ey}{\underline{e}_y}
+\DeclareMathOperator{\ez}{\underline{e}_z}
+\DeclareMathOperator{\bragg}{\theta_{bragg}}
+\DeclareMathOperator{\Z}{Z_C}
+\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 circle (magnetic axis)}
+\author{Didier VEZINET}
+\date{04.12.2019}
+\maketitle
+
+\tableofcontents
+
+\chapter{Geometry}
+\label{chap:Geometry}
+
+\section{Generic cone and plane}
+
+Let's consider a cartesian frame $(O, \ex, \ey, \ez)$.
+Let's consider a half-cone defined by its axis $(S, \n)$ and half-opening $\alpha = \pi/2 - \bragg$.
+The coordinates of $S$ are $(x_S, y_S, z_S)$.
+The coordinates of $\n$ are $(n_x, n_y, n_z)$.
+Let's consider a circle of axis $(O, \ez)$, of radius $R$, centered on $C$ or coordinates $(0, 0, \Z)$.
+
+Let's consider point $M$ of coordinates $(x, y, z)$ and $(R, \theta, z)$ belonging to both the cone and the circle.
+
+
+$$
+\lt\{
+	\begin{array}{ll}
+		M \in circle & \lra \ud{OM} = \Z\ez + R(\cos(\theta)\ex + \sin(\theta)\ey)\\
+		M \in cone & \lra \ud{SM}.\n = \cos(\alpha)\|\ud{SM}\|
+	\end{array}
+\rt.
+$$
+
+
+
+\section{Intersection}
+
+If $M$ belongs to both the circle and the cone, then:
+
+$$
+\begin{array}{ll}
+	& \lt[\lt(\ud{SO}+\ud{OM}\rt).\n\rt]^2 = \cos^2(\alpha)\|\ud{SO} + \ud{OM}\|^2\\
+	\lra &
+	\lt(\ud{SO}.\n\rt)^2 + \lt(\ud{OM}.\n\rt)^2 + 2\lt(\ud{SO}.\n\rt)\lt(\ud{OM}.\n\rt)
+	= \cos^2(\alpha)\lt[\|\ud{SO}\|^2 + \|\ud{OM}\|^2 + 2\ud{SO}.\ud{OM}\rt]\\
+	\lra &
+	\lt(\ud{OM}.\n\rt)^2 + 2\lt(\ud{SO}.\n\rt)\lt(\ud{OM}.\n\rt)
+	- \cos^2(\alpha)\|\ud{OM}\|^2 - 2\cos^2(\alpha)\ud{SO}.\ud{OM} + A = 0
+\end{array}
+$$
+
+Where we have introduced $A = \lt(\ud{SO}.\n\rt)^2 - \cos^2(\alpha)\|\ud{SO}\|^2$
+
+Now, we can write:
+$$
+\lt\{
+	\begin{array}{lll}
+		\|\ud{OM}\|^2 & = & \Z^2 + R^2\\
+		\ud{OM}.\n & = & \Z n_z + R\cos(\theta)n_x + R\sin(\theta)n_y\\
+		\lt(\ud{OM}.\n\rt)^2 & = & \lt(\Z n_z\rt)^2 + \lt(R\cos(\theta)n_x\rt)^2 + \lt(R\sin(\theta)n_y\rt)^2\\
+				     &&  + 2\Z R\cos(\theta)n_xn_Z + 2\Z R\sin(\theta)n_yn_z + 2R^2\cos(\theta)\sin(\theta)n_xn_y\\
+		\ud{SO}.\ud{OM} & = & -\Z z_S - Rx_S\cos(\theta) - Ry_S\sin(\theta)
+	\end{array}
+\rt.
+$$
+
+Hence:
+$$
+\begin{array}{ll}
+	& \lt(\ud{OM}.\n\rt)^2 + 2\lt(\ud{SO}.\n\rt)\lt(\ud{OM}.\n\rt)
+	- \cos^2(\alpha)\|\ud{OM}\|^2 - 2\cos^2(\alpha)\ud{SO}.\ud{OM} + A\\
+	= &
+	\lt(\Z n_z\rt)^2 + \lt(R\cos(\theta)n_x\rt)^2 + \lt(R\sin(\theta)n_y\rt)^2\\
+	& + 2\Z R\cos(\theta)n_xn_Z + 2\Z R\sin(\theta)n_yn_z + 2R^2\cos(\theta)\sin(\theta)n_xn_y\\
+	& +2\lt(\ud{SO}.\n\rt)\lt(\Z n_z + R\cos(\theta)n_x + R\sin(\theta)n_y\rt)\\
+	& - \cos^2(\alpha)\Z^2 - \cos^2(\alpha)R^2\\
+	& + 2\cos^2(\alpha)\lt(\Z z_S + Rx_S\cos(\theta) + Ry_S\sin(\theta)\rt) + A
+\end{array}
+$$
+
+
+
+
+
+
+
+\section{Parametric equation}
+
+\subsection{From bragg angle and parameter to local cartesian coordinates}
+
+
+
+
+
+
+
+\appendix
+\chapter{Appendices}
+
+\section{Section}
+\subsection{Subsection}
+
+\end{document}

tofu/geom/_comp_optics.py

@@ -106,25 +106,70 @@ def get_lamb_from_bragg(bragg, d, n=None):
 #           Approximate solution
 # ###############################################
 
-def get_approx_detector_rel(rcurve, bragg):
+def get_approx_detector_rel(rcurve, bragg, tangent_to_rowland=None):
+
+    if tangent_to_rowland is None:
+        tangent_to_rowland = True
 
     # 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
+    n_crystdet_rel = np.r_[-np.sin(bragg), np.cos(bragg), 0.]
+    if tangent_to_rowland:
+        bragg2 = 2.*bragg
+        det_nout_rel = np.r_[-np.cos(bragg2), -np.sin(bragg2), 0.]
+        det_ei_rel = np.r_[np.sin(bragg2), -np.cos(bragg2), 0.]
+    else:
+        det_nout_rel = -n_crystdet_rel
+        det_ei_rel = np.r_[np.cos(bragg), np.sin(bragg), 0]
+    return det_dist, n_crystdet_rel, det_nout_rel, det_ei_rel
 
-def get_det_abs_from_rel(det_dist, det_nout_rel, det_ei_rel,
-                         summit, nout, e1, e2):
+
+def get_det_abs_from_rel(det_dist, n_crystdet_rel, det_nout_rel, det_ei_rel,
+                         summit, nout, e1, e2,
+                         ddist=None, di=None, dj=None,
+                         dtheta=None, dpsi=None, tilt=None):
+    # Reference
     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
+
+    # Apply translation of center (ddist, di, dj)
+    if ddist is None:
+        ddist = 0.
+    if di is None:
+        di = 0.
+    if dj is None:
+        dj = 0.
+    det_dist += ddist
+
+    n_crystdet = (n_crystdet_rel[0]*nout
+                  + n_crystdet_rel[1]*e1 + n_crystdet_rel[2]*e2)
+    det_cent = summit + det_dist*n_crystdet + di*det_ei + dj*det_ej
+
+    # Apply angles on unit vectors with respect to themselves
+    if dtheta is None:
+        dtheta = 0.
+    if dpsi is None:
+        dpsi = 0.
+    if tilt is None:
+        tilt = 0.
+
+    # dtheta and dpsi
+    det_nout2 = ((np.cos(dpsi)*det_nout
+                 + np.sin(dpsi)*det_ei)*np.cos(dtheta)
+                 + np.sin(dtheta)*det_ej)
+    det_ei2 = (np.cos(dpsi)*det_ei - np.sin(dpsi)*det_nout)
+    det_ej2 = np.cross(det_nout2, det_ei2)
+
+    # tilt
+    det_ei3 = np.cos(tilt)*det_ei2 + np.sin(tilt)*det_ej2
+    det_ej3 = np.cross(det_nout2, det_ei3)
+
+    return det_cent, det_nout2, det_ei3, det_ej3
 
 
 # ###############################################
@@ -144,6 +189,7 @@ def checkformat_vectang(Z, nn, frame_cent, 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)
@@ -155,6 +201,7 @@ def get_e1e2_detectorplane(nn, nIn):
     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)
@@ -168,35 +215,115 @@ def calc_xixj_from_braggphi(summit, det_cent, det_nout, det_ei, det_ej,
     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]
+def calc_braggphi_from_xixjpts(det_cent, det_ei, det_ej,
+                               summit, nin, e1, e2,
+                               xi=None, xj=None, pts=None):
+
+    if pts is None:
+        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
     else:
+        assert pts.ndim == 2
+        pts = pts[:, :, None]
         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
+
+
+# ###############################################
+#           From plasma pts
+# ###############################################
+
+def calc_psidthetaphi_from_pts_lamb(pts, center, rcurve,
+                                    bragg, nlamb, npts,
+                                    nout, e1, e2, extenthalf, ntheta=None):
+
+    if ntheta is None:
+        ntheta = 100
+
+    scaPCem = np.full((nlamb, npts), np.nan)
+    dtheta = np.full((nlamb, npts, ntheta), np.nan)
+    psi = np.full((nlamb, npts, ntheta), np.nan)
+
+    # Get to scalar product
+    PC = center[:, None] - pts
+    PCnorm2 = np.sum(PC**2, axis=0)
+    cos2 = np.cos(bragg)**2
+    deltaon4 = (rcurve**2*cos2[:, None]**2
+                - (rcurve**2*cos2[:, None]
+                   - PCnorm2[None, :]*np.sin(bragg)[:, None]**2))
+
+    # Get two relevant solutions
+    ind = deltaon4 >= 0.
+    cos2 = np.repeat(cos2[:, None], npts, axis=1)[ind]
+    PCnorm = np.tile(np.sqrt(PCnorm2), (nlamb, 1))[ind]
+    sol1 = -rcurve*cos2 - np.sqrt(deltaon4[ind])
+    sol2 = -rcurve*cos2 + np.sqrt(deltaon4[ind])
+    # em is a unit vector and ...
+    ind1 = (np.abs(sol1) <= PCnorm) & (sol1 >= -rcurve)
+    ind2 = (np.abs(sol2) <= PCnorm) & (sol2 >= -rcurve)
+    assert not np.any(ind1 & ind2)
+    sol1 = sol1[ind1]
+    sol2 = sol2[ind2]
+    indn = ind.nonzero()
+    ind1 = [indn[0][ind1], indn[1][ind1]]
+    ind2 = [indn[0][ind2], indn[1][ind2]]
+    scaPCem[ind1[0], ind1[1]] = sol1
+    scaPCem[ind2[0], ind2[1]] = sol2
+    ind = ~np.isnan(scaPCem)
+
+    # Get equation on PCem
+    X = np.sum(PC*nout[:, None], axis=0)
+    Y = np.sum(PC*e1[:, None], axis=0)
+    Z = np.sum(PC*e2[:, None], axis=0)
+
+    scaPCem = np.repeat(scaPCem[..., None], ntheta, axis=-1)
+    ind = ~np.isnan(scaPCem)
+    XYnorm = np.repeat(np.repeat(np.sqrt(X**2 + Y**2)[None, :],
+                                 nlamb, axis=0)[..., None],
+                       ntheta, axis=-1)[ind]
+    Z = np.repeat(np.repeat(Z[None, :], nlamb, axis=0)[..., None],
+                  ntheta, axis=-1)[ind]
+    angextra = np.repeat(
+        np.repeat(np.arctan2(Y, X)[None, :], nlamb, axis=0)[..., None],
+        ntheta, axis=-1)[ind]
+    dtheta[ind] = np.repeat(
+        np.repeat(extenthalf[1]*np.linspace(-1, 1, ntheta)[None, :],
+                  npts, axis=0)[None, ...],
+        nlamb, axis=0)[ind]
+
+    psi[ind] = (np.arccos(
+        (scaPCem[ind] - Z*np.sin(dtheta[ind]))/(XYnorm*np.cos(dtheta[ind])))
+                + angextra)
+    psi[ind] = np.arctan2(np.sin(psi[ind]), np.cos(psi[ind]))
+    indnan = (~ind) | (np.abs(psi) > extenthalf[0])
+    psi[indnan] = np.nan
+    dtheta[indnan] = np.nan
+    return dtheta, psi, indnan

tofu/geom/_core.py

@@ -4141,6 +4141,11 @@ def _complete_dX12(self, dgeom):
                         pinhole = dgeom['D'][:, 0] + k[0]*u[:, 0]
                         dgeom['pinhole'] = pinhole
 
+            if np.any(np.isnan(dgeom['D'])):
+                msg = ("Some LOS have nan as starting point !\n"
+                       + "The geometry may not be provided !")
+                raise Exception(msg)
+
             # Test if all D are on a common plane or line
             va = dgeom["D"] - dgeom["D"][:, 0:1]
 
@@ -4209,9 +4214,6 @@ def _complete_dX12(self, dgeom):
                         dgeom["dX12"] = {}
                     dgeom["dX12"].update({"nIn": nIn, "e1": e1, "e2": e2})
 
-                    if not self._is2D():
-                        return dgeom
-
                     # Test binning
                     if dgeom["pinhole"] is not None:
                         k1ref = -np.sum(