Cylindrical crystal + tuto on CrystalBragg

[Didou09]

Motivations:

• The CrystalBragg class should handle cylindrical crystals
• Most of the routines should be easy to adapt, but a thorough review of all of them is necessary to make sure non is left aside

Main routines to update:

• The CrystalBragg class is defined in tofu/geom/_core_optics.py

• Don't be frightened away by the look of the class, it was coded quite some time away and there is some legacy code in there, but I'll help and all this will be refactored / cleaned-up / harmonized on the mid-term (some of my initial ideas turned out not so useful and others need implementing).

• The main methods to start with are:

  • set_dgeom(), which calls _check_optics._checkformat_dgeom()
  • set_dmat(), which calls _check_optics._checkformat_dmat()
  • set_dbragg(), which calls _check_optics._checkformat_dbragg()

These are called at instanciation.

• Then you'll have to check that some utility methods behave properly (and update them if necessary):
  • sample_outline_plot(): samples the outline of the crystal, for plotting
  • get_local_noute1e2(): compute, for each ray impacting the crystal, the local unit vectors defining a direct orthonormal frame (nout, e1, e2) locally at each impact point, taking into acccount the local curvature and optionally a miscut.
  • _get_rays_from_cryst() and get_rays_from_cryst(): used for ray-tracing
  • split(): to split the crystal in halves (useful only for WEST, low-priority)
  • plot(): plotting
  • calc_meridional_sagital_focus(): get meridional / sagital distances from summit
  • calc_xixj_from_braggphi(): compute the local coordinates on a detector (xi, xj) of rays coming from the crystal to the detector.
  • plot_line_on_det_tracing():
  • calc_johannerror(): rayt-tracing on the whole surface of the crystal, to a detector, and computation of the maximum wavelength discrepancy for each pixel.
  • get_lamb_avail_from_pts(): for each poits in the plasma, compute the wavelength interval that is accessible as seen from the cryst + detector pair.
  • calc_raytracing_from_lambpts(): ray-tracing from plasma points at given wavelengths

There are other methods, but these are the main ones (we'll see later if we add some more). There are quite a few, but I think that if the lower-level ones are correctly handled, most of the upper-level ones will work automatically, so actually I don't think there is so much work.

• Then, some methods might require some more adaptation:
  • get_detector_ideal(): compute the position of an ideally-positionned detector for a given wavelenght (and a given source position, the real detector position is usually a small variation in the neighbourhood if that ideal detector)

You will find methods about rocking curve and fitting, no need to look at these for the moment (the ones about rocking curves are within the scope of issue #654.

Suggestions:

• Create a dedicated branch Issue656_CylindricalCrystal to address this issue
• @cjperks7 proposed to take care of this, with support from me @Didou09

Current use of the CrystalBragg class:

A CrystalBragg instance is defined based on 3 dictionaries: - dgeom: holds the purely geometrical parameters (type of geometry, position, orientation with unit vectors, dimensions as half-extent, ...) - dmat: holds everything specific to the material of the crystal (formula, density, mesh symmetry, mesh lengths and angles, Miller indices of the cut, optional miscut angles, orientation of the lattice - if there is a miscut and the lattice is not parallel to the optical surface - with unit vectors...) - dbragg: holds everything specific to bragg diffraction (tabulated rocking curve, wavelength of reference...)

All positions / unit vectors in tofu are always defined in (X, Y, Z) cartesian coordinates taken from the absolute frame of the tokamak (where, at a toroidal angle of phi = 0, the unit vector of X equals the unit vector of R).

Some of the quantities in these dictionaries do not need to be provided by the user and are instead automatically computed from others (ex: unit vectors of the lattice orientation are computed for the values of miscut angles alpha and beta, which are 0 by default).

In addition to the 3 disctionnary, the current version of tofu asks for an Exp (name of the tokamak, in upper case), Diag (name fo the diagnostic) and a Name (name of that particular crystal, in case your diagnostic has several crystals). This may evolve in the future.

Direct copy-pasting into an ipython console should work for all the snippets below, provided you don't include the In [1]:.

# import tofu
In [1]: import tofu as tf

The crystal is oriented in 3d using a direct orthonormal base (nout, e1, e2) at its summit

# Defining the direct orthonormal base
In [2]: nout = np.r_[-0.39925082, -0.91684071, -0.0013777 ]
In [3]: e1 = np.r_[ -nout[1], nout[0], 0]
In [4]: e2 = np.cross(nout, e1)

Typically one of the unit vectors is horizontal (here e1).

# Defining dgeom
In [5]: dgeom = {
    ...:     'Type': 'sph',             # vs 'cyl', 'flat'
    ...:     'Typeoutline': 'rect',     # vs 'circ'
    ...:     # summit = center of the crystal (the keyword 'center' is reserved for another use)
    ...:     'summit': np.r_[4.6497750e-01, -8.8277925e+00,  3.5125000e-03],   # X, Y, Z
    ...:     # center = center of curvature, the center of the sphere for a spherical crystal
    ...:     'center': np.r_[1.560921  , -6.31106476,  0.00729429],    # X, Y, Z
    ...:     # extenthalf: half-width of the crystal, in both directions
    ...:     'extenthalf': [0.01457195, 0.01821494],       # rad in a curved direction, m in a linear direction
    ...:     # 'surface' will be computed
    ...:     # 'nout': normalized unit vector normal to the crystal's surface at its summit, pointing outwards (away from the sphere / cylinder center)
    ...:     'nout': nout,
    ...:     # 'nin' = -nout will be computed ('nin' is pointing inwards, towards the center of the sphere / cylinder)
    ...:     # 'e1': normalized unit vector, tangent to the crystal's surface at its summit, perpendicular to 'nin'
    ...:     'e1': e1,
    ...:     # 'e2': normalized unit vector, forming a direct orthonormal basis (nout, e1, e2)
    ...:     'e2': e2,
    ...:     'rcurve': 2.745,     # radius of curvature
    ...:     'move': None, # optional, in case your crystal can move and that movement can be parameterized by a single variable, name of the function to call
    ...: }

# Defining dmat
In [6]: dmat = {
    ...:     'formula': 'Quartz',
    ...:     'density': 2.6576,
    ...:     'symmetry': 'hexagonal',
    ...:     # 3 basis lengths of the hexagonal mesh
    ...:     'lengths': [4.9079e-10, 4.9079e-10, 5.3991e-10],
    ...:     # 3 basis angles of the hexagonal mesh
    ...:     'angles': [1.57079633, 1.57079633, 2.0943951],
    ...:     # Miller indices of the cut
    ...:     'cut': [1,1,0],
    ...:     'd': 2.4539499999999996e-10, # interplanar distance, is computed
    ...:     'alpha': 0., # miscut angle
    ...:     'beta': 0., # miscut orientation
    ...:     # 'nin', 'e1', 'e2' are computed from alpha and beta
    ...: }

# Defining dbragg
In [7]: dbragg = {
    ...:     'rockingcurve': {
    ...:         'dangle': None,   # 1d array of angles delta (rad) around the bragg angle
    ...:         'value': None,    # 1d array of reflected intensity, same shape as dangle
    ...:         'source': None,   # str, indicating where the rtabulated rocking curve was taken from
    ...:         'lamb': None,     # float (m) the wavelength for which the rocking curve was computed
    ...:         'type': None,     # str (ex: 'tabulated-1d', indicating the type of rocking curve
    ...:         },
    ...:     'lambref': 3.96e-10,  # wavelength of reference (i.e.: the default wavelength used for ray-tracing when the wavelength is not specified by the user
    ...:     # 'braggref'  bragg angle of reference, computed from 'lambref'
    ...: }

Instanciation and use:

In [8]: cryst = tf.geom.CrystalBragg(dgeom=dgeom, dmat=dmat, dbragg=dbragg, Exp='WEST', Diag='XICS', Name='ArXVII')

# just typing the variable pretty-prints a summary of the crystal's characterisctics
In [9]: cryst
Out[9]: 
formula  symmetry   cut      density  d (A)  bragg(     3.96 A) (deg)  Type  outline  surface (cm²)  rcurve  rocking curve
-------  ---------  -------  -------  -----  ------------------------  ----  -------  -------------  ------  -------------
Quartz   hexagonal  [1 1 0]  2.6576   2.454  53.79050296650186         sph   rect      80.0           2.745  None         

half-extent    summit               center               nout                    e1                      alpha  beta  color               
-------------  -------------------  -------------------  ----------------------  ----------------------  -----  ----  --------------------
[0.015 0.018]  [ 0.46 -8.83  0.  ]  [ 1.56 -6.31  0.01]  [-0.399 -0.917 -0.001]  [ 0.917 -0.399  0.   ]  0.0    0.0   (0.0, 0.0, 0.0, 1.0)

Plotting: the plotting method (plot()) provides a static figure with:

• A 2d poloidal projection
• A 2d horizontal (seen from above) projection
• A 3d view
• A view on an optional detector (for ray-tracing)

This method will have to be adapted for cylindrical crystal because by default it plots a rowland circle. It returns a dict of axes, and may raise some harmless warnings.

In [10]: dax = cryst.plot()

Optionally, it can include the plot of a detector (defined using a dict with center position cent, unit vectors (nout, ei, ej) and 2d outline as (2, npts) array), I recommend you save it as .npz and load it when needed (in a near-future a specific class will be created).

Optionnaly you can also add ray-tracing from pts in the plasma at chosen wavelength.

# build pts in the plasma along a line taken from the crystal summit at 3 different distances form it:
In [11]: pts0, vect = cryst.get_rays_from_cryst(phi=np.pi, returnas='(pts, vect)')

In [12]: pts = pts0 + np.r_[1.5, 3, 7]*vect[:, 0:1, 0]

# load a predefined detector
In [13]: det = dict(np.load('tofu_west/SpectroX2D/Inputs/det37_CTVD_incC4_New.npz', allow_pickle=True))

# plot for 2 wavelengths
In [14]: dax = cryst.plot(pts=pts, det=det, lamb=np.r_[3.95e-10, 3.96e-10, 3.97e-10])

The ray-tracing is computed for the grid {all pts} x {all lamb}. Coloring is done: - per pts (rays_color='pts', default) - per wavelength (rays_color='lamb')

[cjperks7]

Hey Didier, wanted to flag this error I get when I run essentially this same code you provided, but edit dgeom to dgeom{'Type' : 'flat'}. Seems to be an issue with sample_outline_plot() having an if statement with instructions only for the spherical crystal. In general, I get the sense from reading the source code that all the plotting tools are optimized for the spherical crystal. Isn't surprising.

I wanted to start with a flat crystal because I expected it to be so simple everything already worked. If you're aware of a workaround then feel to let me know. I'd have to deal with the same plotting upgrades for the cylindrical crystal and would like to validate my math in the limit of the crystal becoming flat, so not an inconvenience if you'd just recommend getting the plotting for both to work

[Didou09]

No problem, will try to handle this tomorrow morning

[cjperks7]

@Didou09 I'm going to make some changes in the ray-tracing calculations to handle flat and cylindrical crystal reflections. This will be focusing on functions such as get_local_noute1e2(), get_rays_from_cryst(), _get_rays_from_cryst(), and calc_raytracing_from_lambpts(). I won't change much that's already there for the spherical crystal. Just encompass what's there with an if statement based on the geometry.

I will need to change how the code handles points on the crystal surface as it's hardcoded to be in terms of two rotational distances (dtheta, psi) from the summit, so native to spherical. The cylindrical crystal will do something similar, but in terms of a rotational and translational distance (will define the cylinder axis to be along e1) while the flat will be in terms of two translational variables (along e1, e2). Since these variables are often passed within many functions, I think it might be easiest to group these variables into a single dictionary to be passed around with the variables contained returning to their original usage within the relevant geometry.