diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_ConeCircleMagAxis.pdf b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConeCircleMagAxis.pdf new file mode 100644 index 000000000..136140d9c Binary files /dev/null and b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConeCircleMagAxis.pdf differ diff --git a/Notes_Upgrades/SpectroX2D/SpectroX2D_ConeCircleMagAxis.tex b/Notes_Upgrades/SpectroX2D/SpectroX2D_ConeCircleMagAxis.tex new file mode 100644 index 000000000..65fad2cb5 --- /dev/null +++ b/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} diff --git a/TFG_CrystalBragg_ExpWEST_DgXICS_ArXVII_sh00000_Vers1.4.1-174-g453d6a3.npz b/TFG_CrystalBragg_ExpWEST_DgXICS_ArXVII_sh00000_Vers1.4.1-174-g453d6a3.npz new file mode 100644 index 000000000..1620a6b6c Binary files /dev/null and b/TFG_CrystalBragg_ExpWEST_DgXICS_ArXVII_sh00000_Vers1.4.1-174-g453d6a3.npz differ diff --git a/tofu/geom/_comp_optics.py b/tofu/geom/_comp_optics.py index 4917835b5..3f31e95d1 100644 --- a/tofu/geom/_comp_optics.py +++ b/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,25 +215,32 @@ 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)) @@ -194,9 +248,82 @@ def calc_braggphi_from_xixj(xi, xj, det_cent, det_ei, det_ej, 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 diff --git a/tofu/geom/_core.py b/tofu/geom/_core.py index 4e4021a11..f1384dc2e 100644 --- a/tofu/geom/_core.py +++ b/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( diff --git a/tofu/geom/_core_optics.py b/tofu/geom/_core_optics.py index e58138373..8762194d9 100644 --- a/tofu/geom/_core_optics.py +++ b/tofu/geom/_core_optics.py @@ -13,6 +13,7 @@ # Common import numpy as np +import scipy.interpolate as scpinterp import datetime as dtm import matplotlib.pyplot as plt import matplotlib as mpl @@ -292,17 +293,17 @@ def _checkformat_dmat(cls, dmat=None): 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: + if dmat.get('d') is not None: dmat['d'] = float(dmat['d']) - if dmat['formula'] is not None: + if dmat.get('formula') is not None: assert isinstance(dmat['formula'], str) - if dmat['density'] is not None: + if dmat.get('density') is not None: dmat['density'] = float(dmat['density']) - if dmat['lengths'] is not None: + if dmat.get('lengths') is not None: dmat['lengths'] = np.atleast_1d(dmat['lengths']).ravel() - if dmat['angles'] is not None: + if dmat.get('angles') is not None: dmat['angles'] = np.atleast_1d(dmat['angles']).ravel() - if dmat['cut'] is not None: + if dmat.get('cut') is not None: dmat['cut'] = np.atleast_1d(dmat['cut']).ravel().astype(int) assert dmat['cut'].size <= 4 return dmat @@ -312,12 +313,69 @@ def _checkformat_dbragg(cls, dbragg=None): if dbragg is None: return assert isinstance(dbragg, dict) - lkok = ['angle'] + lkok = cls._get_keys_dbragg() 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]) + if dbragg.get('rockingcurve') is not None: + assert isinstance(dbragg['rockingcurve'], dict) + drock = dbragg['rockingcurve'] + lkeyok = ['sigam', 'deltad', 'Rmax', 'dangle', 'lamb', 'value', + 'type', 'source'] + lkout = [kk for kk in drock.keys() if kk not in lkeyok] + if len(lkout) > 0: + msg = ("Unauthorized keys in dbrag['rockingcurve']:\n" + + "\t-" + "\n\t-".join(lkout)) + raise Exception(msg) + try: + if drock.get('sigma') is not None: + dbragg['rockingcurve']['sigma'] = float(drock['sigma']) + dbragg['rockingcurve']['deltad'] = float( + drock.get('deltad', 0.)) + dbragg['rockingcurve']['Rmax'] = float( + drock.get('Rmax', 1.)) + dbragg['rockingcurve']['type'] = 'lorentz-log' + elif drock.get('dangle') is not None: + c2d = (drock.get('lamb') is not None + and drock.get('value').ndim == 2) + if c2d: + if drock['value'].shape != (drock['dangle'].size, + drock['lamb'].size): + msg = ("Tabulated 2d rocking curve should be:\n" + + "\tshape = (dangle.size, lamb.size)") + raise Exception(msg) + dbragg['rockingcurve']['dangle'] = np.r_[ + drock['dangle']] + dbragg['rockingcurve']['lamb'] = np.r_[drock['lamb']] + dbragg['rockingcurve']['value'] = drock['value'] + dbragg['rockingcurve']['type'] = 'tabulated-2d' + else: + if drock.get('lamb') is None: + msg = ("Please also specify the lamb for which " + + "the rocking curve was tabulated!") + raise Exception(msg) + dbragg['rockingcurve']['lamb'] = float(drock['lamb']) + dbragg['rockingcurve']['dangle'] = np.r_[ + drock['dangle']] + dbragg['rockingcurve']['value'] = np.r_[drock['value']] + dbragg['rockingcurve']['type'] = 'tabulated-1d' + if drock.get('source') is None: + msg = "Unknonw source for the tabulated rocking curve!" + warnings.warn(msg) + dbragg['rockingcurve']['source'] = drock.get('source') + except Exception as err: + msg = ("Provide the rocking curve as a dict with either:\n" + + "\t- parameters of a lorentzian in log10:\n" + + "\t\t{'sigma': float,\n" + + "\t\t 'deltad': float,\n" + + "\t\t 'Rmax': float}\n" + + "\t- tabulated (dangle, value), with source (url...)" + + "\t\t{'dangle': np.ndarray,\n" + + "\t\t 'value': np.ndarray,\n" + + "\t\t 'source': str}") + raise Exception(msg) return dbragg @classmethod @@ -404,7 +462,7 @@ def set_dmat(self, dmat=None): dmat = self._checkformat_dmat(dmat) self._dmat = dmat - def _set_dbragg(self, dbragg=None): + def set_dbragg(self, dbragg=None): dbragg = self._checkformat_dbragg(dbragg) self._dbragg = dbragg @@ -570,6 +628,13 @@ def summit(self): def center(self): return self._dgeom['center'] + @property + def rockingcurve(self): + if self._dbragg.get('rockingcurve') is not None: + if self._dbragg['rockingcurve'].get('type') is not None: + return self._dbragg['rockingcurve'] + raise Exception("rockingcurve was not set!") + # ----------------- # methods for color # ----------------- @@ -591,11 +656,17 @@ def get_summary(self, sep=' ', line='-', just='l', # ----------------------- # Build material col0 = ['formula', 'symmetry', 'cut', 'density', - 'd (A)', 'bragg(%s A) (deg)'%str(self._DEFLAMB)] + 'd (A)', 'bragg({:7.4} A) (deg)'.format(self._DEFLAMB*1e10), + 'rocking curve'] 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)] + try: + ar0.append(self.rockingcurve['type']) + except Exception as err: + ar0.append('None') + # ----------------------- # Build geometry @@ -692,6 +763,54 @@ def move(self, kind=None, **kwdargs): if kind == 'rotate': self._rotate(**kwdargs) + # ----------------- + # methods for rocking curve + # ----------------- + + def get_rockingcurve_func(self, lamb=None, n=None): + drock = self.rockingcurve + if drock['type'] == 'tabulated-1d': + if lamb is not None and lamb != drock['lamb']: + msg = ("rocking curve was tabulated only for:\n" + + "\tlamb = {} m\n".format(lamb) + + " => Please let lamb=None") + raise Exception(msg) + bragg = self._checkformat_bragglamb(lamb=drock['lamb'], n=n) + func = scpinterp.interp1d(drock['dangle']+bragg, drock['value'], + kind='linear', bounds_error=False, + fill_value=0, assume_sorted=True) + + elif drock['type'] == 'tabulated-2d': + lmin, lmax = drock['lamb'].min(), drock['lamb'].max() + if lamb is None or lamb < lmin or lamb > lmax: + msg = ("rocking curve was tabulated only in interval:\n" + + "\tlamb in [{}; {}] m\n".format(lmin, lmax) + + " => Please set lamb accordingly") + raise Exception(msg) + bragg = self._checkformat_bragglamb(lamb=lamb, n=n) + + def func(angle, lamb=lamb, bragg=bragg, drock=drock): + return scpinterp.interp2d(drock['dangle']+bragg, drock['lamb'], + drock['value'], kind='linear', + bounds_error=False, fill_value=0, + assume_sorted=True)(angle, lamb) + + else: + def func(angle, d=d, delta_bragg=delta_bragg, + Rmax=drock['Rmax'], sigma=drock['sigma']): + core = sigma**2/((angle - (bragg+delta_bragg))**2 + sigma**2) + if Rmax is None: + return core/(sigma*np.pi) + else: + return Rmax*core + return func + + def plot_rockingcurve(self, lamb=None, + fs=None, ax=None): + drock = self.rockingcurve + func = self.get_rockingcurve_func(bragg=bragg, n=n) + return _plot.CrystalBragg_plot_rockingcurve(func, fs=fs, ax=ax) + # ----------------- # methods for surface and contour sampling # ----------------- @@ -729,6 +848,9 @@ def sample_outline_Rays(self, res=None): def _checkformat_bragglamb(self, bragg=None, lamb=None, n=None): lc = [lamb is not None, bragg is not None] + if not any(lc): + lamb = self._DEFLAMB + lc[0] = True assert np.sum(lc) == 1, "Provide lamb xor bragg!" if lc[0]: bragg = self.get_bragg_from_lamb(np.atleast_1d(lamb), @@ -821,7 +943,8 @@ def get_Rays_envelop(self, # ----------------- def plot(self, lax=None, proj=None, res=None, element=None, - color=None, + color=None, det_cent=None, + det_nout=None, det_ei=None, det_ej=None, dP=None, dI=None, dBs=None, dBv=None, dVect=None, dIHor=None, dBsHor=None, dBvHor=None, dleg=None, @@ -836,6 +959,19 @@ def plot(self, lax=None, proj=None, res=None, element=None, # methods for generic first-approx # ----------------- + def get_phi_from_magaxis_summit(self, r, z, lamb=None, bragg=None, n=None): + # Check / format input + r = np.atleast_1d(r) + z = np.atleast_1d(z) + assert r.shape == z.shape + bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) + + # Compute phi + + return phi + + + def get_bragg_from_lamb(self, lamb, n=None): """ Braggs' law: n*lamb = 2dsin(bragg) """ if self._dmat['d'] is None: @@ -854,27 +990,66 @@ def get_lamb_from_bragg(self, bragg, n=None): 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 """ + def get_detector_approx(self, bragg=None, lamb=None, + rcurve=None, n=None, + ddist=None, di=None, dj=None, + dtheta=None, dpsi=None, tilt=None, + tangent_to_rowland=None, plot=False): + """ Return approximate ideal detector geometry + + Assumes infinitesimal and ideal crystal + Assumes detector center tangential to Rowland circle + Assumes detector center matching lamb (m) / bragg (rad) + + Detector described by center position, and (nout, ei, ej) unit vectors + By convention, nout = np.cross(ei, ej) + Vectors (ei, ej) define an orthogonal frame in the detector's plane + + Return: + ------- + det_cent: np.ndarray + (3,) array of (x, y, z) coordinates of detector center + det_nout: np.ndarray + (3,) array of (x, y, z) coordinates of unit vector + perpendicular to detector' surface + oriented towards crystal + det_ei: np.ndarray + (3,) array of (x, y, z) coordinates of unit vector + defining first coordinate in detector's plane + det_ej: np.ndarray + (3,) array of (x, y, z) coordinates of unit vector + defining second coordinate in detector's plane + """ # Check / format inputs if rcurve is None: rcurve = self._dgeom['rcurve'] bragg = self._checkformat_bragglamb(bragg=bragg, lamb=lamb, n=n) + if np.all(np.isnan(bragg)): + msg = ("There is no available bragg angle!\n" + + " => Check the vlue of self.dmat['d'] vs lamb") + raise Exception(msg) + + lf = ['summit', 'nout', 'e1', 'e2'] + lc = [rcurve is None] + [self._dgeom[kk] is None for kk in lf] + if any(lc): + msg = ("Some missing fields in dgeom for computation:" + + "\n\t-" + "\n\t-".join(['rcurve'] + lf)) + raise Exception(msg) # 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) + (det_dist, n_crystdet_rel, + det_nout_rel, det_ei_rel) = _comp_optics.get_approx_detector_rel( + rcurve, bragg, tangent_to_rowland=tangent_to_rowland) # 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']) + det_cent, det_nout, det_ei, det_ej = _comp_optics.get_det_abs_from_rel( + det_dist, n_crystdet_rel, det_nout_rel, det_ei_rel, + self._dgeom['summit'], + self._dgeom['nout'], self._dgeom['e1'], self._dgeom['e2'], + ddist=ddist, di=di, dj=dj, + dtheta=dtheta, dpsi=dpsi, tilt=tilt) if plot: dax = self.plot() @@ -891,6 +1066,31 @@ def get_approx_detector_frame(self, bragg=None, lamb=None, return det_cent, det_nout, det_ei, det_ej def get_local_noute1e2(self, theta, psi): + """ Return (nout, e1, e2) associated to pts on the crystal's surface + + All points on the spherical crystal's surface are identified + by (theta, psi) coordinates, where: + - theta = np.pi/2 for the center + - psi = 0 for the center + They are the spherical coordinates from a sphere centered on the + crystal's center of curvature. + + Return the pts themselves and the 3 perpendicular unti vectors + (nout, e1, e2), where nout is towards the outside of the sphere and + nout = np.cross(e1, e2) + + Return: + ------- + summit: np.ndarray + (3,) array of (x, y, z) coordinates of the points on the surface + nout: np.ndarray + (3,) array of (x, y, z) coordinates of outward unit vector + e1: np.ndarray + (3,) array of (x, y, z) coordinates of first unit vector + e2: np.ndarray + (3,) array of (x, y, z) coordinates of second unit vector + + """ if np.allclose([theta, psi], [np.pi/2., 0.]): summit = self._dgeom['summit'] nout = self._dgeom['nout'] @@ -946,7 +1146,7 @@ def calc_xixj_from_phibragg(self, phi=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 + func = self.get_detector_approx det_cent, det_nout, det_ei, det_ej = func(lamb=self._DEFLAMB) # Get local summit nout, e1, e2 if non-centered @@ -967,18 +1167,45 @@ def calc_xixj_from_phibragg(self, phi=None, ax = func(bragg, xi, xj, data, ax) return xi, xj + @staticmethod + def _checkformat_pts(pts): + pts = np.atleast_1d(pts) + if pts.ndim == 1: + pts = pts.reshape((3, 1)) + if 3 not in pts.shape or pts.ndim != 2: + msg = "pts must be a (3, npts) array of (X, Y, Z) coordinates!" + raise Exception(msg) + if pts.shape[0] != 3: + pts = pts.T + return pts + @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) + @staticmethod + def _checkformat_psidtheta(psi=None, dtheta=None, + psi_def=0., dtheta_def=0.): + if psi is None: + psi = psi_def + if dtheta is None: + dtheta = dtheta_def + psi = np.r_[psi] + dtheta = np.r_[dtheta] + if psi.size != dtheta.size: + msg = "psi and dtheta must be 1d arrays fo same size!" + raise Exception(msg) + + def calc_phibragg_from_xixj(self, xi, xj, n=None, det_cent=None, det_ei=None, det_ej=None, - theta=None, psi=None, + dtheta=None, psi=None, plot=True, ax=None, **kwdargs): # Check / format inputs @@ -987,38 +1214,309 @@ def calc_phibragg_from_xixj(self, xi, xj, n=None, 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) + det_cent, _, det_ei, det_ej = self.get_detector_approx( + 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) + bragg, phi = self.calc_phibragg_from_pts(pts) 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) + lax = _plot_optics.CrystalBragg_plot_braggangle_from_xixj( + 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 + # DEPRECATED ??? + def calc_phibragg_from_pts_on_summit(self, pts, n=None): + """ Return the bragg angle and phi of pts from crystal summit - 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): + The pts are provided as a (x, y, z) coordinates array + The bragg angle and phi are computed from the crystal's summit + + """ + # Check / format inputs + pts = self._checkformat_pts(pts) + + # Compute + vect = pts - self._dgeom['summit'][:, None] + vect = vect / np.sqrt(np.sum(vect**2, axis=0)) + bragg = np.pi/2 - np.arccos(np.sum(vect*self._dgeom['nin'][:, None], + axis=0)) + v1 = np.sum(vect*self._dgeom['e1'][:, None], axis=0) + v2 = np.sum(vect*self._dgeom['e2'][:, None], axis=0) + phi = np.arctan2(v2, v1) + return bragg, phi + + def plot_line_tracing_on_det(self, lamb=None, n=None, + xi_bounds=None, xj_bounds=None, nphi=None, + det_cent=None, det_nout=None, + det_ei=None, det_ej=None, + johann=False, lpsi=None, ltheta=None, + rocking=False, fs=None, dmargin=None, + wintit=None, tit=None): + """ Visualize the de-focusing by ray-tracing of chosen lamb + """ + # Check / format inputs + if lamb is None: + lamb = self._DEFLAMB + lamb = np.atleast_1d(lamb).ravel() + nlamb = lamb.size + + det = np.array([[xi_bounds[0], xi_bounds[1], xi_bounds[1], + xi_bounds[0], xi_bounds[0]], + [xj_bounds[0], xj_bounds[0], xj_bounds[1], + xj_bounds[1], xj_bounds[0]]]) + + # Compute lamb / phi + _, phi = self.calc_phibragg_from_xixj( + det[0, :], det[1, :], n=n, + det_cent=det_cent, det_ei=det_ei, det_ej=det_ej, + theta=None, psi=None, plot=False) + phimin, phimax = np.nanmin(phi), np.nanmax(phi) + phimin, phimax = phimin-(phimax-phimin)/10, phimax+(phimax-phimin)/10 + del phi + + # Get reference ray-tracing + if nphi is None: + nphi = 300 + phi = np.linspace(phimin, phimax, nphi) + bragg = self._checkformat_bragglamb(lamb=lamb, n=n) + + xi = np.full((nlamb, nphi), np.nan) + xj = np.full((nlamb, nphi), np.nan) + for ll in range(nlamb): + xi[ll, :], xj[ll, :] = self.calc_xixj_from_phibragg( + bragg=bragg[ll], phi=phi, n=n, + det_cent=det_cent, det_nout=det_nout, + det_ei=det_ei, det_ej=det_ej, plot=False) + + # Get johann-error raytracing (multiple positions on crystal) + xi_er, xj_er = None, None + if johann and not rocking: + if lpsi is None or ltheta is None: + lpsi = np.linspace(-1., 1., 15) + ltheta = np.linspace(-1., 1., 15) + lpsi, ltheta = np.meshgrid(lpsi, ltheta) + lpsi = lpsi.ravel() + ltheta = ltheta.ravel() + lpsi = self._dgeom['extenthalf'][0]*np.r_[lpsi] + ltheta = np.pi/2 + self._dgeom['extenthalf'][1]*np.r_[ltheta] + npsi = lpsi.size + assert npsi == ltheta.size + + xi_er = np.full((nlamb, npsi*nphi), np.nan) + xj_er = np.full((nlamb, npsi*nphi), np.nan) + for l in range(nlamb): + for ii in range(npsi): + i0 = np.arange(ii*nphi, (ii+1)*nphi) + xi_er[l, i0], xj_er[l, i0] = self.calc_xixj_from_phibragg( + phi=phi, bragg=bragg[l], lamb=None, n=n, + theta=ltheta[ii], psi=lpsi[ii], + det_cent=det_cent, det_nout=det_nout, + det_ei=det_ei, det_ej=det_ej, plot=False) + + # Get rocking curve error + if rocking: + pass + + # Plot + ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( + lamb, xi, xj, xi_er, xj_er, det=det, + johann=johann, rocking=rocking, + fs=fs, dmargin=dmargin, wintit=wintit, tit=tit) + + def calc_johannerror(self, xi=None, xj=None, err=None, + det_cent=None, det_ei=None, det_ej=None, n=None, + lpsi=None, ltheta=None, + plot=True, fs=None, cmap=None, + vmin=None, vmax=None, tit=None, wintit=None): + """ Plot the johann error + + The johann error is the error (scattering) induced by defocalization + due to finite crystal dimensions + There is a johann error on wavelength (lamb => loss of spectral + resolution) and on directionality (phi) + If provided, lpsi and ltheta are taken as normalized variations with + respect to the crystal summit and to its extenthalf. + Typical values are: + - lpsi = [-1, 1, 1, -1] + - ltheta = [-1, -1, 1, 1] + They must have the same len() + + """ + + # Check / format inputs + 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 + bragg, phi = self.calc_phibragg_from_xixj( + xii, xjj, n=n, + det_cent=det_cent, det_ei=det_ei, det_ej=det_ej, + theta=None, psi=None, plot=False) + assert bragg.shape == phi.shape + lamb = self.get_lamb_from_bragg(bragg, n=n) + + if lpsi is None: + lpsi = np.r_[-1., 0., 1., 1., 1., 0., -1, -1] + lpsi = self._dgeom['extenthalf'][0]*np.r_[lpsi] + if ltheta is None: + ltheta = np.r_[-1., -1., -1., 0., 1., 1., 1., 0.] + ltheta = np.pi/2 + self._dgeom['extenthalf'][1]*np.r_[ltheta] + npsi = lpsi.size + assert npsi == ltheta.size + lamberr = np.full(tuple(np.r_[npsi, lamb.shape]), np.nan) + phierr = np.full(lamberr.shape, np.nan) + for ii in range(npsi): + bragg, phierr[ii, ...] = self.calc_phibragg_from_xixj( + xii, xjj, n=n, + det_cent=det_cent, det_ei=det_ei, det_ej=det_ej, + theta=ltheta[ii], psi=lpsi[ii], plot=False) + lamberr[ii, ...] = self.get_lamb_from_bragg(bragg, n=n) + err_lamb = np.nanmax(np.abs(lamb[None, ...] - lamberr), axis=0) + err_phi = np.nanmax(np.abs(phi[None, ...] - phierr), axis=0) + if plot is True: + ax = _plot_optics.CrystalBragg_plot_johannerror( + xi, xj, lamb, phi, err_lamb, err_phi, err=err, + cmap=cmap, vmin=vmin, vmax=vmax, fs=fs, tit=tit, wintit=wintit) + return err_lamb, err_phi + + def calc_phibragg_from_pts(self, pts, dtheta=None, dpsi=None): + + # Check / format pts + pts = self._checkformat_pts(pts) + + # Get local summit nout, e1, e2 if non-centered + psi, dtheta = self._checkformat_psidtheta(psi=psi, dtheta=dtheta) + summit, nout, e1, e2 = self.get_local_noute1e2(dtheta, psi) + + # Compute + bragg, phi = _comp_optics.calc_braggphi_from_xixjpts( + det_cent, det_ei, det_ej, + summit, -nout, e1, e2, pts=pts) + return phi, bragg + + def get_lamb_avail_from_pts(self, pts): + pass + + def calc_thetapsi_from_lambpts(self, pts=None, lamb=None, n=None, + ntheta=None): + + # Check / Format inputs + pts = self._checkformat_pts(pts) + npts = pts.shape[1] + + if lamb is None: + lamb = self._DEFLAMB + lamb = np.r_[lamb] + nlamb = lamb.size + bragg = self.get_bragg_from_lamb(lamb, n=n) + + dtheta, psi, indnan = _comp_optics.calc_psidthetaphi_from_pts_lamb( + pts, self._dgeom['center'], self._dgeom['rcurve'], + bragg, nlamb, npts, + self._dgeom['nout'], self._dgeom['e1'], self._dgeom['e2'], + self._dgeom['extenthalf'], ntheta=ntheta) + + # import ipdb; ipdb.set_trace() # DB + bragg = np.repeat(np.repeat(bragg[:, None], npts, axis=-1)[..., None], + ntheta, axis=-1) + bragg[indnan] = np.nan + phi[ind] = self.calc_braggphi_from_pts(pts, + dtheta=dtheta[ind], + dpsi=dpsi[ind]) + + return dtheta, psi, phi, bragg + + def plot_line_from_pts_on_det(self, lamb=None, pts=None, + xi_bounds=None, xj_bounds=None, nphi=None, + det_cent=None, det_nout=None, + det_ei=None, det_ej=None, + johann=False, lpsi=None, ltheta=None, + rocking=False, fs=None, dmargin=None, + wintit=None, tit=None): + """ Visualize the de-focusing by ray-tracing of chosen lamb + """ + # Check / format inputs + if lamb is None: + lamb = self._DEFLAMB + lamb = np.atleast_1d(lamb).ravel() + nlamb = lamb.size + + det = np.array([[xi_bounds[0], xi_bounds[1], xi_bounds[1], + xi_bounds[0], xi_bounds[0]], + [xj_bounds[0], xj_bounds[0], xj_bounds[1], + xj_bounds[1], xj_bounds[0]]]) + + # Compute lamb / phi + _, phi = self.calc_phibragg_from_xixj( + det[0, :], det[1, :], n=n, + det_cent=det_cent, det_ei=det_ei, det_ej=det_ej, + theta=None, psi=None, plot=False) + phimin, phimax = np.nanmin(phi), np.nanmax(phi) + phimin, phimax = phimin-(phimax-phimin)/10, phimax+(phimax-phimin)/10 + del phi + + # Get reference ray-tracing + if nphi is None: + nphi = 300 + phi = np.linspace(phimin, phimax, nphi) + bragg = self._checkformat_bragglamb(lamb=lamb, n=n) + + xi = np.full((nlamb, nphi), np.nan) + xj = np.full((nlamb, nphi), np.nan) + for ll in range(nlamb): + xi[ll, :], xj[ll, :] = self.calc_xixj_from_phibragg( + bragg=bragg[ll], phi=phi, n=n, + det_cent=det_cent, det_nout=det_nout, + det_ei=det_ei, det_ej=det_ej, plot=False) + + # Get johann-error raytracing (multiple positions on crystal) + xi_er, xj_er = None, None + if johann and not rocking: + if lpsi is None or ltheta is None: + lpsi = np.linspace(-1., 1., 15) + ltheta = np.linspace(-1., 1., 15) + lpsi, ltheta = np.meshgrid(lpsi, ltheta) + lpsi = lpsi.ravel() + ltheta = ltheta.ravel() + lpsi = self._dgeom['extenthalf'][0]*np.r_[lpsi] + ltheta = np.pi/2 + self._dgeom['extenthalf'][1]*np.r_[ltheta] + npsi = lpsi.size + assert npsi == ltheta.size + + xi_er = np.full((nlamb, npsi*nphi), np.nan) + xj_er = np.full((nlamb, npsi*nphi), np.nan) + for l in range(nlamb): + for ii in range(npsi): + i0 = np.arange(ii*nphi, (ii+1)*nphi) + xi_er[l, i0], xj_er[l, i0] = self.calc_xixj_from_phibragg( + phi=phi, bragg=bragg[l], lamb=None, n=n, + theta=ltheta[ii], psi=lpsi[ii], + det_cent=det_cent, det_nout=det_nout, + det_ei=det_ei, det_ej=det_ej, plot=False) + + # Get rocking curve error + if rocking: + pass + + # Plot + ax = _plot_optics.CrystalBragg_plot_line_tracing_on_det( + lamb, xi, xj, xi_er, xj_er, det=det, + johann=johann, rocking=rocking, + fs=fs, dmargin=dmargin, wintit=wintit, tit=tit) + + def plot_data_vs_lambphi(self, xi=None, xj=None, data=None, mask=None, + det_cent=None, det_ei=None, det_ej=None, + theta=None, psi=None, n=None, + nlambfit=None, nphifit=None, + magaxis=None, npaxis=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) @@ -1026,24 +1524,58 @@ def plot_data_vs_lambphi(self, xi=None, xj=None, data=None, 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) + bragg, phi = self.calc_phibragg_from_xixj( + 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) + if mask is not None: + data[~mask] = np.nan + if nlambfit is None: + nlambfit = nxi + if nphifit is None: + nphifit = nxj + lambfit, phifit = _comp_optics.get_lambphifit(lamb, phi, + nlambfit, nphifit) 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)]) + spect1d = np.array([np.nanmean(data[ind == ii]) + for ii in np.unique(ind)]) + phifitbins = 0.5*(phifit[1:] + phifit[:-1]) + ind = np.digitize(phi, phifitbins) + vertsum1d = np.array([np.nanmean(data[ind == ii]) + for ii in np.unique(ind)]) + + # Get phiref from mag axis + lambax, phiax = None, None + if magaxis is not None: + if npaxis is None: + npaxis = 1000 + thetacryst = np.arctan2(self._dgeom['summit'][1], + self._dgeom['summit'][0]) + thetaax = thetacryst + np.pi/2*np.linspace(-1, 1, npaxis) + pts = np.array([magaxis[0]*np.cos(thetaax), + magaxis[0]*np.sin(thetaax), + np.full((npaxis,), magaxis[1])]) + braggax, phiax = self.calc_phibragg_from_pts(pts) + lambax = self.get_lamb_from_bragg(braggax) + phiax = np.arctan2(np.sin(phiax-np.pi), np.cos(phiax-np.pi)) + ind = ((lambax >= lambfit[0]) & (lambax <= lambfit[-1]) + & (phiax >= phifit[0]) & (phiax <= phifit[-1])) + lambax, phiax = lambax[ind], phiax[ind] + ind = np.argsort(lambax) + lambax, phiax = lambax[ind], phiax[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) + if plot: + ax = _plot_optics.CrystalBragg_plot_data_vs_lambphi( + xi, xj, bragg, lamb, phi, data, + lambfit=lambfit, phifit=phifit, spect1d=spect1d, + vertsum1d=vertsum1d, lambax=lambax, phiax=phiax, + cmap=cmap, vmin=vmin, vmax=vmax, fs=fs) return ax def plot_data_fit2d(self, xi=None, xj=None, data=None, mask=None, diff --git a/tofu/geom/_plot_optics.py b/tofu/geom/_plot_optics.py index 07e669218..9c35ad098 100644 --- a/tofu/geom/_plot_optics.py +++ b/tofu/geom/_plot_optics.py @@ -93,6 +93,7 @@ def _check_projdax_mpl(dax=None, proj=None, fs=None, wintit=None): def CrystalBragg_plot(cryst, lax=None, proj=None, res=None, element=None, color=None, dP=None, + det_cent=None, det_nout=None, det_ei=None, det_ej=None, dI=None, dBs=None, dBv=None, dVect=None, dIHor=None, dBsHor=None, dBvHor=None, dleg=None, indices=False, @@ -118,8 +119,10 @@ def CrystalBragg_plot(cryst, lax=None, proj=None, res=None, element=None, # Temporary matplotlib issue dleg = None else: - dax = _CrystalBragg_plot_crosshor(cryst, proj=proj, res=res, dax=lax, element=element, - color=color) + dax = _CrystalBragg_plot_crosshor(cryst, proj=proj, res=res, dax=lax, + element=element, color=color, + det_cent=det_cent, det_nout=det_nout, + det_ei=det_ei, det_ej=det_ej) # recompute the ax.dataLim ax0 = None @@ -139,7 +142,11 @@ def CrystalBragg_plot(cryst, lax=None, proj=None, res=None, element=None, ax0.figure.canvas.draw() return dax -def _CrystalBragg_plot_crosshor(cryst, proj=None, dax=None, element=None, res=None, + +def _CrystalBragg_plot_crosshor(cryst, proj=None, dax=None, + element=None, res=None, + det_cent=None, det_nout=None, + det_ei=None, det_ej=None, Pdict=_def.TorPd, Idict=_def.TorId, Bsdict=_def.TorBsd, Bvdict=_def.TorBvd, Vdict=_def.TorVind, color=None, ms=None, quiver_cmap=None, @@ -232,21 +239,77 @@ def _CrystalBragg_plot_crosshor(cryst, proj=None, dax=None, element=None, res=No 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,:], + pr = np.hypot(p0[0, :], p0[1, :]) + vr = np.hypot(p0[0, :]+v[0, :], p0[1, :]+v[1, :]) - pr + dax['cross'].quiver(pr, p0[2, :], + vr, 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,:], + 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) + # Detector + sc = None + if det_cent is not None: + if dax['cross'] is not None: + dax['cross'].plot(np.hypot(det_cent[0], det_cent[1]), det_cent[2], + marker='x', ms=ms, c=color, label="det_cent") + if dax['hor'] is not None: + dax['hor'].plot(det_cent[0], det_cent[1], + marker='x', ms=ms, c=color, label="det_cent") + + if det_nout is not None: + assert det_ei is not None and det_ej is not None + p0 = np.repeat(det_cent[:, None], 3, axis=1) + v = np.concatenate((det_nout[:, None], det_ei[:, None], + det_ej[:, None]), axis=1) + if dax['cross'] is not None: + pr = np.hypot(p0[0, :], p0[1, :]) + vr = np.hypot(p0[0, :]+v[0, :], p0[1, :]+v[1, :]) - pr + dax['cross'].quiver(pr, p0[2, :], + vr, v[2, :], + np.r_[0., 0.5, 1.], cmap=quiver_cmap, + angles='xy', scale_units='xy', + label="det 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="det unit vect", **Vdict) + return dax +# ################################################################# +# ################################################################# +# Rocking curve plot +# ################################################################# +# ################################################################# + +def CrystalBragg_plot_rockingcurve(Rmax=None, sigma=None, + bragg=None, delta_bragg=None, npts=None): + + # Prepare + if npts is None: + npts = 1000 + angle = bragg + delta_bragg + 3.*sigma*np.linspace(-1, 1, npts) + curve = func(angle) + + # Plot + if ax is None: + fig = plt.figure() + ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) + ax.plot(angle, curve, ls='-', lw=1., c='k') + ax.axvline(bragg, ls='--', lw=1, c='k') + return ax + + # ################################################################# # ################################################################# # Bragg diffraction plot @@ -342,8 +405,140 @@ def CrystalBragg_plot_braggangle_from_xixj(xi=None, xj=None, return ax +def CrystalBragg_plot_line_tracing_on_det(lamb, xi, xj, xi_err, xj_err, + det=None, + johann=None, rocking=None, + fs=None, dmargin=None, + wintit=None, tit=None): + + # Check inputs + # ------------ + + if fs is None: + fs = (6, 8) + if dmargin is None: + dmargin = {'left': 0.05, 'right': 0.99, + 'bottom': 0.06, 'top': 0.92, + 'wspace': None, 'hspace': 0.4} + + if wintit is None: + wintit = _WINTIT + if tit is None: + tit = "line tracing" + if johann is True: + tit += " - johann error" + if rocking is True: + tit += " - rocking curve" + + plot_err = johann is True or rocking is True + + # Plot + # ------------ + + fig = fig = plt.figure(figsize=fs) + gs = gridspec.GridSpec(1, 1, **dmargin) + ax0 = fig.add_subplot(gs[0, 0], aspect='equal', adjustable='datalim') + + ax0.plot(det[0, :], det[1, :], ls='-', lw=1., c='k') + for l in range(lamb.size): + lab = r'$\lambda$'+' = {:6.3f} A'.format(lamb[l]*1.e10) + l0, = ax0.plot(xi[l, :], xj[l, :], ls='-', lw=1., label=lab) + if plot_err: + ax0.plot(xi_err[l, ...], xj_err[l, ...], + ls='None', lw=1., c=l0.get_color(), + marker='.', ms=1) + + ax0.legend() + + if wintit is not False: + fig.canvas.set_window_title(wintit) + if tit is not False: + fig.suptitle(tit, size=14, weight='bold') + return [ax0] + + +def CrystalBragg_plot_johannerror(xi, xj, lamb, phi, err_lamb, err_phi, + cmap=None, vmin=None, vmax=None, + fs=None, dmargin=None, wintit=None, tit=None, + angunits='deg', err=None): + + # Check inputs + # ------------ + + if fs is None: + fs = (14, 8) + if cmap is None: + cmap = plt.cm.viridis + if dmargin is None: + dmargin = {'left': 0.05, 'right': 0.99, + 'bottom': 0.06, '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 + err_phi = err_phi*180./np.pi + + if err is None: + err = 'abs' + if err == 'rel': + err_lamb = 100.*err_lamb / (np.nanmax(lamb) - np.nanmin(lamb)) + err_phi = 100.*err_phi / (np.nanmax(phi) - np.nanmin(phi)) + err_lamb_units = '%' + err_phi_units = '%' + else: + err_lamb_units = 'm' + err_phi_units = angunits + + if wintit is None: + wintit = _WINTIT + if tit is None: + tit = False + + # pre-compute + # ------------ + + # extent + extent = (xi.min(), xi.max(), xj.min(), xj.max()) + + # Plot + # ------------ + + fig = fig = plt.figure(figsize=fs) + gs = gridspec.GridSpec(1, 3, **dmargin) + ax0 = fig.add_subplot(gs[0, 0], aspect='equal', adjustable='datalim') + ax1 = fig.add_subplot(gs[0, 1], aspect='equal', adjustable='datalim', + sharex=ax0, sharey=ax0) + ax2 = fig.add_subplot(gs[0, 2], aspect='equal', adjustable='datalim', + sharex=ax0, sharey=ax0) + + ax0.set_title('Iso-lamb and iso-phi at crystal summit') + ax1.set_title('Focalization error on lamb ({})'.format(err_lamb_units)) + ax2.set_title('Focalization error on phi ({})'.format(err_phi_units)) + + ax0.contour(xi, xj, lamb, 10, cmap=cmap) + ax0.contour(xi, xj, phi, 10, cmap=cmap, ls='--') + imlamb = ax1.imshow(err_lamb, extent=extent, aspect='equal', + origin='lower', interpolation='nearest', + vmin=vmin, vmax=vmax) + imphi = ax2.imshow(err_phi, extent=extent, aspect='equal', + origin='lower', interpolation='nearest', + vmin=vmin, vmax=vmax) + + plt.colorbar(imlamb, ax=ax1) + plt.colorbar(imphi, ax=ax2) + if wintit is not False: + fig.canvas.set_window_title(wintit) + if tit is not False: + fig.suptitle(tit, size=14, weight='bold') + + return [ax0, ax1, ax2] + + def CrystalBragg_plot_data_vs_lambphi(xi, xj, bragg, lamb, phi, data, - lambfit=None, phifit=None, spect1d=None, + lambfit=None, phifit=None, + spect1d=None, vertsum1d=None, + lambax=None, phiax=None, cmap=None, vmin=None, vmax=None, fs=None, dmargin=None, angunits='deg'): @@ -364,6 +559,9 @@ def CrystalBragg_plot_data_vs_lambphi(xi, xj, bragg, lamb, phi, data, bragg = bragg*180./np.pi phi = phi*180./np.pi phifit = phifit*180./np.pi + if phiax is not None: + phiax = 180*phiax/np.pi + # pre-compute @@ -377,29 +575,36 @@ def CrystalBragg_plot_data_vs_lambphi(xi, xj, bragg, lamb, phi, data, # ------------ fig = fig = plt.figure(figsize=fs) - gs = gridspec.GridSpec(4, 3, **dmargin) + gs = gridspec.GridSpec(4, 4, **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], sharey=ax2) ax0.set_title('Coordinates transform') ax1.set_title('Camera image') ax2.set_title('Camera image transformed') - ax0.set_ylabel(r'incidence angle ($deg$)') + ax2.set_ylabel(r'incidence angle ($deg$)') + ax2.set_xlabel(r'$\lambda$ ($m$)') + axs2.set_xlabel(r'$\lambda$ ($m$)') + ax3.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='-') + axs1.plot(xi, np.nanmean(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='-') + ax3.plot(vertsum1d, phifit, c='k', ls='-') + if phiax is not None: + ax2.plot(lambax, phiax, c='r', ls='-', lw=1.) ax2.set_xlim(extent2[0], extent2[1]) ax2.set_ylim(extent2[2], extent2[3]) diff --git a/tofu/imas2tofu/_core.py b/tofu/imas2tofu/_core.py index 2a73d5be6..612a17b67 100644 --- a/tofu/imas2tofu/_core.py +++ b/tofu/imas2tofu/_core.py @@ -1508,10 +1508,10 @@ def add_ids_base(self, occ=None, idd=None, user=user, tokamak=tokamak, version=version, ref=ref, isget=isget, get=get) - def add_ids_for_synthdiag(self, ids=None, occ=None, idd=None, - shot=None, run=None, refshot=None, refrun=None, - user=None, tokamak=None, version=None, - ref=None, isget=None, get=None): + def add_ids_synthdiag(self, ids=None, occ=None, idd=None, + shot=None, run=None, refshot=None, refrun=None, + user=None, tokamak=None, version=None, + ref=None, isget=None, get=None): """ Add pre-tabulated input ids necessary for calculating synth. signal The necessary input ids are given by self.get_inputs_for_synthsignal() @@ -3978,11 +3978,11 @@ def _save_to_imas_Config( obj, idd=None, shotfile=None, raise Exception(msg) if description_2d is None: - if nS == 1 and lcls[0] in ['Ves','PlasmaDomain']: + if nS == 1 and lcls[0] in ['Ves', 'PlasmaDomain']: description_2d = 0 else: - descrption_2d = 2 - assert description_2d in [0,2] + description_2d = 2 + assert description_2d in [0, 2] # Make sure StructIn is last (IMAS requirement) ind = lcls.index(lclsIn[0]) diff --git a/tofu/utils.py b/tofu/utils.py index 06784657a..4b6289428 100644 --- a/tofu/utils.py +++ b/tofu/utils.py @@ -142,7 +142,7 @@ def flatten_dict(d, parent_key='', sep=None, deep='ref', if k not in lexcept_key: if issubclass(v.__class__, ToFuObjectBase): if deep=='dict': - v = v.to_dict(deep='dict') + v = v.to_dict(sep=sep, deep='dict') elif deep=='copy': v = v.copy(deep='copy') new_key = parent_key + sep + k if parent_key else k @@ -320,7 +320,15 @@ def save(obj, path=None, name=None, sep=None, deep=False, mode='npz', # Get stripped dictionnary deep = 'dict' if deep else 'ref' if sep is None: - sep = _SEP + if mode == 'mat': + sep = '_' + else: + sep = _SEP + if mode == 'mat' and sep == '.': + msg = ("sep='.' cannot be used when mode='mat' (incompatible)\n" + + "Matlab would interpret variables as structures") + raise Exception(msg) + dd = obj.to_dict(strip=strip, sep=sep, deep=deep) pathfileext = os.path.join(path,name+'.'+mode) @@ -639,7 +647,7 @@ def _get_exception(q, ids, qtype='quantity'): def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, - ids=None, Name=None, out=None, tlim=None, config=None, + ids=None, Name=None, returnas=None, tlim=None, config=None, occ=None, indch=None, indDescription=None, equilibrium=None, dsig=None, data=None, X=None, t0=None, dextra=None, plot=True, plot_sig=None, plot_X=None, @@ -657,14 +665,16 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, raise Exception(msg) lok = ['Config', 'Plasma2D', 'Cam', 'Data'] - c0 = out is None or out in lok + c0 = returnas is None or returnas in lok if not c0: - msg = "Arg out must be in %s"%str(lok) + msg = "Arg returnas must be in {}".format(str(lok)) raise Exception(msg) # ------------------- # Prepare ids - assert ids is None or type(ids) in [list,str] + if type(ids) not in [list, str]: + msg = "Please specify an ids to load data from!" + raise Exception(msg) if type(ids) is str: ids = [ids] if type(ids) is list: @@ -826,17 +836,17 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, # ------------------- - # Prepare out - loutok = ['Config','Plasma2D','Cam','Data'] - c0 = out is None - c1 = out in loutok - c2 = type(out) is list and all([oo is None or oo in loutok - for oo in out]) + # Prepare returnas + loutok = ['Config', 'Plasma2D', 'Cam', 'Data'] + c0 = returnas is None + c1 = returnas in loutok + c2 = type(returnas) is list and all([oo is None or oo in loutok + for oo in returnas]) assert c0 or c1 or c2 if c0: - out = [None for _ in ids] + returnas = [None for _ in ids] elif c1: - out = [str(out) for _ in ids] + returnas = [str(returnas) for _ in ids] # Temporary caveat if nids > 1: @@ -852,35 +862,36 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, # Config if ids[ii] == 'wall': - assert out[ii] in [None,'Config'] - out[ii] = 'Config' - if out[ii] == 'Config': - assert ids[ii] in [None,'wall'] + assert returnas[ii] in [None, 'Config'] + returnas[ii] = 'Config' + if returnas[ii] == 'Config': + assert ids[ii] in [None, 'wall'] # Plasma2D lids = imas2tofu.MultiIDSLoader._lidsplasma if ids[ii] in lids: - assert out[ii] in [None,'Plasma2D'] - out[ii] = 'Plasma2D' - if out[ii] == 'Plasma2D': + assert returnas[ii] in [None, 'Plasma2D'] + returnas[ii] = 'Plasma2D' + if returnas[ii] == 'Plasma2D': assert ids[ii] in lids # Cam or Data lids = imas2tofu.MultiIDSLoader._lidsdiag if ids[ii] in lids: - assert out[ii] in [None,'Cam','Data'] - if out[ii] is None: - out[ii] = 'Data' - if out[ii] in ['Cam','Data']: + assert returnas[ii] in [None, 'Cam', 'Data'] + if returnas[ii] is None: + returnas[ii] = 'Data' + if returnas[ii] in ['Cam', 'Data']: assert ids[ii] in lids - dout = {shot[jj]: {oo:[] for oo in set(out)} for jj in range(0,nshot)} + dout = {shot[jj]: {oo: [] for oo in set(returnas)} + for jj in range(0, nshot)} # ------------------- # Prepare plot_ and complement ids - lPla = [ii for ii in range(0,nids) if out[ii] == 'Plasma2D'] - lCam = [ii for ii in range(0,nids) if out[ii] == 'Cam'] - lDat = [ii for ii in range(0,nids) if out[ii] == 'Data'] + lPla = [ii for ii in range(0, nids) if returnas[ii] == 'Plasma2D'] + lCam = [ii for ii in range(0, nids) if returnas[ii] == 'Cam'] + lDat = [ii for ii in range(0, nids) if returnas[ii] == 'Data'] nPla, nCam, nDat = len(lPla), len(lCam), len(lDat) if nDat > 1: plot_ = False @@ -960,23 +971,23 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, # export to instances for ii in range(0,nids): - if out[ii] == 'Config': + if returnas[ii] == 'Config': dout[ss]['Config'].append(multi.to_Config(Name=Name, occ=occ, indDescription=indDescription, plot=False)) - elif out[ii] == 'Plasma2D': + elif returnas[ii] == 'Plasma2D': dout[ss]['Plasma2D'].append(multi.to_Plasma2D(Name=Name, occ=occ, tlim=tlim, dsig=dsig, t0=t0, plot=False, plot_sig=plot_sig, dextra=dextra, plot_X=plot_X, config=config, bck=bck)) - elif out[ii] == 'Cam': + elif returnas[ii] == 'Cam': dout[ss]['Cam'].append(multi.to_Cam(Name=Name, occ=occ, ids=lids[ii], indch=indch, config=config, plot=False)) - elif out[ii] == "Data": + elif returnas[ii] == "Data": dout[ss]['Data'].append(multi.to_Data(Name=Name, occ=occ, ids=lids[ii], tlim=tlim, dsig=dsig, config=config, data=data, X=X, indch=indch, @@ -990,7 +1001,7 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, # Config & Cam for ss in shot: - for k0 in set(['Config','Cam']).intersection(out): + for k0 in set(['Config', 'Cam']).intersection(returnas): for ii in range(0, len(dout[ss][k0])): dout[ss][k0][ii].plot() @@ -1020,7 +1031,7 @@ def load_from_imas(shot=None, run=None, user=None, tokamak=None, version=None, dout = dout[shot[0]]['Data'][0] elif nshot == 1 and nPla == 1: dout = dout[shot[0]]['Plasma2D'][0] - return out + return dout @@ -1625,7 +1636,7 @@ def to_dict(self, strip=None, sep=None, deep='ref'): # Call class-specific dd = self._to_dict() # --------------------- - dd['dId'] = self._get_dId() + dd['dId'] = self._get_dId(sep=sep) dd['dstrip'] = {'dict':self._dstrip, 'lexcept':None} dout = {} @@ -1646,7 +1657,7 @@ def to_dict(self, strip=None, sep=None, deep='ref'): dout = flatten_dict(dout, parent_key='', sep=sep, deep=deep) return dout - def _get_dId(self): + def _get_dId(self, sep=None): """ To be overloaded """ return {'dict':{}} @@ -1875,8 +1886,8 @@ def Id(self): """ return self._Id - def _get_dId(self): - return {'dict':self.Id.to_dict()} + def _get_dId(self, sep=None): + return {'dict': self.Id.to_dict(sep=sep)} def _reset(self): if hasattr(self,'_Id'): diff --git a/tofu/version.py b/tofu/version.py index 52ca3f491..4fc937bc2 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.2-a5-15-g231dacd' +__version__ = '1.4.2-a5-48-gac0016b'