Add yaml definitions of MARKER, MATRIX and MCDESINT

zgoubi_metadata/data/elements_yaml/MARKER.yaml

@@ -1,3 +1,9 @@
 zgoubi_name: MARKER
 params:
 doc: |
+  Marker
+
+  .. rubric:: Zgoubi manual description
+
+  Just a marker.
+  No data ’.plt’ as a second LABEL will cause storage of current coordinates into zgoubi.plt

zgoubi_metadata/data/elements_yaml/MATRIX.yaml

@@ -1,8 +1,135 @@
 zgoubi_name: MATRIX
 params:
-- {name: IORD, type: I, default: 0, doc: ""}
-- {name: IFOC, type: I, default: 0, doc: ""}
+- IORD: {type: I, default: 0, doc: ""}
+- IFOC: {type: I, default: 0, doc: ""}
+
 template:
 - - IORD
   - IFOC
 doc: |
+
+  Calculation of transfer coefficients, periodic parameters
+
+  .. rubric:: Zgoubi manual description
+
+  ``MATRIX`` causes the calculation of the transfer coefficients through the optical structure, from the ``OBJET`` down to the location
+  where ``MATRIX`` is introduced in the structure, or, upon option, down to the horizontal focus closest to that location. In this last
+  case the position of the focus is calculated automatically in the same way as the position of the waist in ``IMAGE``. Depending
+  on ``OBJET`` and on option ``IFOC``, ``MATRIX`` also delivers the beam matrix and betatron phase advances or (case of a periodic
+  structure) periodic beam matrix and tunes, chromaticities and other global parameters.
+  Depending on the value of option ``IORD``, different procedures follow
+
+  * If ``IORD`` = 0, ``MATRIX`` is inhibited (equivalent to ``FAISCEAU``, whatever ``IFOC``).
+
+  * If ``IORD`` = 1, using ``OBJET``, ``KOBJ`` = 5[.N], the first order transfer matrix [Rij] is calculated, from a third order approximation
+    of the coordinates. For instance
+
+    .. math::
+
+      \begin{align}
+        Y^+ &= (\frac{Y}{T_0}) T_0 + (\frac{Y}{T_0^2}) T_0^2 +  (\frac{Y}{T_0^3}) T_0^3 \\
+        Y^- &= -(\frac{Y}{T_0}) T_0 + (\frac{Y}{T_0^2}) T_0^2 -  (\frac{Y}{T_0^3}) T_0^3 \\
+      \end{align}
+
+    will yield, neglecting third order terms,
+
+    .. math::
+
+      R_{11} = (\frac{Y}{T_0}) = \frac{Y^+ - Y^-}{2T_0}
+
+    This is repeated N times if ``OBJET``, ``KOBJ`` = 5.N (2 ≤ N ≤ 99) is used, so delivering first order data for each 11-particle
+    set, with for each set the reference trajectory being trajectory number 1+(N-1)x11 in the Nx11-particle list (keyword
+    ``FAISCEAU`` will print out that list in zgoubi.res, if desired).
+    If ``OBJET``, ``KOBJ`` = 5.1 is used instead (hence introducing initial optical function values, :math:`\alpha_{Y, Z}`,
+    :math:`\alpha_{Y, Z}`, :math:`\D_{Y, Z}`, :math:`\D'_{Y, Z}`), then, using the Rij above,
+     ``MATRIX`` will transport the optical functions and phase advances :math:`\phi_{Y}` , :math:`\phi_{Z}`, following
+
+     .. math::
+
+       \begin{pmatrix}
+       \beta \\
+       \alpha \\
+       \gamma
+       \end{pmatrix}_{at\,MATRIX}
+       =
+       \begin{pmatrix}
+       R_{11}^2 & -2R_{11}R_{12} & R_{12}^2 \\
+       -R_{11}R_{21} & R_{12}R_{21} & R_{11}R_{12}  \\
+       R_{21}^2 & -2R_{21}R_{22} & R_{22}^2
+       \end{pmatrix}_{at\,MATRIX}
+       \begin{pmatrix}
+       \beta \\
+       \alpha \\
+       \gamma
+       \end{pmatrix}_{at\,OBJET}
+
+    .. math::
+
+      \begin{align}
+      \Delta \phi_Y &= Atan \frac{R_{12}}{(R_{11} \beta_{Y, objet} - R_{12} \alpha_{Y, objet})} \\
+      \Delta \phi_Z &= Atan \frac{R_{34}}{(R_{33} \beta_{Z, objet} - R_{34} \alpha_{Z, objet})} \\
+      \end{align}
+
+    .. math::
+
+      \phi_{Y, Z} \rightarrow \phi_{Y, Z} + 2\pi \text{ if } \phi_{Y, Z} < 0, \text{ given [0, $\pi$] Atan determination}
+
+    and print these out.
+
+  * If ``IORD`` = 2, using ``OBJET``, ``KOBJ`` = 6, fifth order Taylor expansions are used for the calculation of the first order transfer
+    matrix [Rij ] and of the second order matrix [Tijk]. Other higher order coefficients are also calculated.
+
+  * If ``IORD`` = 3, using ``OBJET``, ``KOBJ`` = 6.1, transport coefficients up to 3rd order are computed using 102 rays (after routines
+    developed for ``RAYTRACE`` [65, 66]).
+
+  An automatic generation of an appropriate object for the use ofMATRIX can be obtained using the procedureOBJET (pages 59, 282),
+  as follows
+
+  * if ``IORD`` = 1, use ``OBJET``(``KOBJ`` = 5[.N, 2 ≤ N ≤ 99]), that generates up to 99*11 initial coordinates. In this case, up to 99
+    matrices may be calculated, each one wrt. to the reference trajectory of concern (trajectory number 1, 12, 23, ... 1+(N-1)x11
+    respectively).
+
+  * if ``IORD`` = 2, use ``OBJET``(``KOBJ`` = 6) that generates 61 initial coordinates.
+
+  * if ``IORD`` = 3, use ``OBJET``(``KOBJ`` = 6.1) that generates 102 initial coordinates
+
+  The next option, IFOC, acts as follows
+
+  * If ``IFOC`` = 0, the transfer coefficients are calculated at the location ofMATRIX, and with respect to the reference trajectory.
+    For instance, :math:`Y^+` and :math:`T^+` above are defined for particle number i as :math:`Y^+` = :math:`Y^+(i) - Y(Ref)`,
+    and :math:`T^+` = :math:`T^+(i) - T(Ref)`.
+
+  * If ``IFOC`` = 1, the transfer coefficients are calculated at the horizontal focus closest to ``MATRIX`` (determined automatically),
+    while the reference direction is that of the reference particle. For instance, :math:`Y^+` is defined for particle number i as
+    :math:`Y^+ = Y^+(i) − Y_{focus}`, while :math:`T^+` is defined as :math:`T^+` = :math:`T^+(i) - T(Ref)`).
+
+    * If ``IFOC`` = 2, no change of reference frame is performed : the coordinates refer to the current frame. Namely, :math:`Y^+ = Y^+(i)`,
+      :math:`T^+ = T^+(i)`, etc.
+
+  ** Periodic Structures **
+
+  * If ``IFOC`` = 10 +NPeriod, then, from the 1-turn transport matrix as obtained in the way described above, ``MATRIX`` calculates
+  periodic parameters characteristic of the structure such as optical functions and tune numbers, assuming that it is NPeriod -
+  periodic. In non-coupled hypothesis, this is performed by a simple identification
+
+  .. math::
+
+    [T_{ij}] = I cos \mu + J sin \mu
+
+  In the coupled hypothesis, the computation uses the Edwards-Teng method [37].
+  Matrix computation is repeated N times if ``OBJET``, ``KOBJ`` = 5.N (2 ≤ N ≤ 99) is used, so delivering first order data for
+  each of the N 11-particle sets.
+
+  If ``IORD`` = 2 (using ``OBJET``, ``KOBJ`` = 6) additional periodic parameters are computed such as chromaticities, beta-function
+  momentum dependence, etc.
+
+  ``PRINT`` : Addition of ``PRINT`` following ``IORD``, ``IFOC`` [, coupled] will cause stacking of ``MATRIX`` output data into
+  zgoubi.MATRIX.out file (convenient for use with e.g. gnuplot type of data treatment software).
+  Addition of “coupled ” next to ``IORD``, ``IFOC`` [, ``PRINT``], in the case of periodic beam matrix
+  request (i.e., ``IFOC`` = 10 + NPeriod ) will cause use of coupled formalism.
+
+  :underline:`About the source code`:
+
+  The program matric computes the transport matrix coefficients, it is called by zgoubi when the keyword ``MATRIX`` is met
+  along the zgoubi.dat sequence. matimp, called by matric, ensures the print out of the matrix (and possibly the beam matrix)
+  in zgoubi.res, and upon ``PRINT`` option in zgoubi.MATRIX.out as well.

zgoubi_metadata/data/elements_yaml/MCDESINT.yaml

@@ -1,14 +1,104 @@
 zgoubi_name: MCDESINT
+
 params:
-- {name: M2, type: E, unit: "", default: 0, doc: ""}
-- {name: M3, type: E, unit: "", default: 0, doc: ""}
-- {name: I1, type: I, default: 0, doc: ""}
-- {name: I2, type: I, default: 0, doc: ""}
-- {name: I3, type: I, default: 0, doc: ""}
+- INFO: {type: A, unit: "", default: 0, doc: "Switch"}
+- M2: {type: E, unit: "MeV_c2", default: 0, doc: "Masses of the two decay products"}
+- M3: {type: E, unit: "MeV_c2", default: 0, doc: "Masses of the two decay products"}
+- TAU2: {type: E, unit:"s", default: 0, doc: "COM lifetime of particle 2"}
+- I1: {type: I, default: 0, doc: "Seeds for random number generators"}
+- I2: {type: I, default: 0, doc: "Seeds for random number generators"}
+- I3: {type: I, default: 0, doc: "Seeds for random number generators"}
+
 template:
-- - M2
-  - M3
+- - [INFO, M2, M3, TAU2]
 - - I1
   - I2
   - I3
+
 doc: |
+
+  Monte-Carlo simulation of in-flight decay
+
+  .. rubric:: Zgoubi manual description
+
+  When ``MCDESINT``[38] is met in a structure (normally, after ``OBJET`` or after ``CIBLE``), in-flight decay simulation
+  starts. It must be preceded by ``PARTICUL`` for the definition of mass ``M1`` and ``COM`` lifetime :math:`\tau_1` of the
+  parent particle.
+  The two-body decay simulated is
+
+  .. math::
+
+    1 \rightarrow  2 + 3
+
+  The decay is isotropic in the center of mass. 1 is the incoming particle, with mass M1, momentum :math:`p1 = \gamma_1 M_1 \beta_1 c`
+  (relative momentum :math:`D1 = \frac{p_1}{q}\frac{1}{BORO} with BORO = reference rigidity, defined in ``[MC]OBJET``).
+  2 and 3 are decay products with respective masses and momenta M2,M3 and :math:`p2 = \gamma_2 M_2 \beta_2 c`, :math:`p3 = \gamma_3 M_3 \beta_3 c`.
+  The average distance s1 at which a particle will decay is related to its center of mass lifetime :math:`\tau_1` by
+
+  .. math::
+
+    s_1 = c \tau_1 \sqrt{\gamma_1^2 -1}
+
+  The actual path length s up to the decay point is calculated, for each one of the particles defined by
+  ``[MC]OBJET`` and prior t ray-tracing, from a random number :math:`0 < R_1 \leq 1` by using the exponential decay
+  formula
+
+  .. math::
+
+    s = s_1 \ln(R_1)
+
+  After decay of the parent particle 1, particle 2 will be ray-traced with assumed positive charge, while
+  particle 3 is discarded. Its scattering angles in the center of mass :math:`\theta^*` and :math:`\phi` are generated from two other
+  random numbers :math:`0 < R_2 \leq 1` and :math:`0 < R_3 \leq 1` by
+
+  .. math::
+
+    \begin{align}
+      \theta^* &= acos(1-2R_2) \text{ ( $0 < \theta^* \leq \pi$)}\\
+      \phi &= 2 \pi R_3 \text{ ( $0 < \phi \leq 2\pi$)}\\
+    \end{align}
+
+  :math:`\phi` is a relativistic invariant, and θ in the laboratory frame (Fig. 43) is given by
+
+  .. math::
+
+    \tan(\theta) = \frac{1}{\gamma_1} \frac{\sin\theta^*}{\frac{\beta_1}{\beta_2^*} + \cos\theta^*}
+
+  :math:`\beta_2^*` and momentum :math:`\p_2` are given by
+
+  .. math::
+
+    \begin{align}
+      \gamma_2^* &= \frac{M_1^2 + M_2^2 - M_3^2}{2M_1M_2}\\
+      \gamma_2 &= \gamm_1\gamma_2^* (1+\beta_1\beta_2\cos\theta^*)\\
+      \beta_2 &= (1 - \frac{1}{\gamma_2^2})^{1/2}\\
+      p_2 = M_2 \sqrt{\gamma_2^2 -1}
+    \end{align}
+
+  Finally, :math:`\theta` and :math:`\phi` are transformed into the angles T2 and P2 in the zgoubi frame, and the relative momentum
+  takes the value :math:`D2 = \frac{p_2}{q}\frac{1}{BORO} with BORO = reference rigidity, see ``OBJET``), while the starting
+  position of M2 is the very location of the parent particle decay, (Y1,Z1, s1).
+  The decay simulation by zgoubi satisfies the following procedures. In optical elements and field maps,
+  after each integration step ``XPAS``, the actual path length of the particle, F(6, I), is compared to its limit path
+  length s. If s is passed, then the particle is considered as having decayed at :math:`F(6, I) − \frac{XPAS}{2}`, at a position
+  obtained by a linear translation from the position at F(6, I). Presumably, the smaller ``XPAS``, the smaller the
+  error on position and angles at the decay point. In ``ESL`` and ``CHANGREF``, F(6, I) is compared to s at the end of the element.
+
+  If the decay occurs inside the element, the particle is considered as having decayed at its actual limit path length s, thus its coordinates at
+  s are recalculated by translation.
+  The limit path length of all particles (I = 1, ``IMAX``) is stored in the array ``FDES`` (6, I). For further statistical
+  purposes (e.g., use of ``HISTO``) the daughter particle 2 is tagged with an ′S′ standing for “secondary”. When
+  a particle decays, its coordinates D, Y , T, Z, P, s, time at the decay point are stored in ``FDES`` (J, I), J = 1, 7.
+
+  * A note on negative drifts* :
+
+  The use of negative drifts with ``MCDESINT`` is allowed and correct. For instance, negative drifts may occur in a structure
+  for some of the particles when using ``CHANGREF`` (due to the Z-axis rotation or a negative ``XCE``),
+  or when using ``DRIFT`` with XL < 0. Provision has been made to take it into account during the
+  ``MCDESINT`` procedure, as follows.
+  If, due to a negative drift, a secondary particle reaches back the decay location of its parent particle, then
+  the parent particle is “resurrected “ with its original coordinates at that location, the secondary particle is
+  discarded, and ray-tracing resumes in a regular way for the parent particle which is again allowed to decay,
+  after the same path length. This procedure is made possible by prior storage of the coordinates of the parent
+  particles (in array ``FDES`` (J, I)) each time a decay occurs.
+  Negative steps (``XPAS``< 0) in optical elements are not compatible with ``MCDESINT``.