Release fixes for 5.09

Changes.md

@@ -9,14 +9,15 @@ MAD-X master
     *   [PR 1119](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1119) Fix tapering issues (R. De Maria **WIP**)
 *   Track:
     *   [PR 1148](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1148) Fix to get Dynap to work with the ny space Charge module (T. Persson)
-    *   [PR 1109](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1109) Space charge (F. Schimdt, H. Renshal, A. Latina)   
+    *   [PR 1109](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1109) Space charge (F. Schimdt, H. Renshal, A. Latina)
+    *   [PR_1159](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1159) Fix calls to pro_input from trupdate, adding null termination (S. Berg)
 *   PTC:
     *   [PR 1153](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1153), [PR 1154](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1154), [PR 1156](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1156) Fix mdump for map debugging and sync with MAD-NG (L. Deniau)
     *   [PR 1142](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1142) Fix issue-1140 generating function wrong sign from ptc_normal (P. Skowonronski) 
     *   [PR 1111](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1111) Ensure that exact_model is set before zero_key is called in PTC (J. S. Berg)
     *   [PR 1095](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1095) Implement more robust, optional, PTC DA map output (L. Deniau)
     *   [PR 1131](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1131) Add spin table (T. Persson and P. Skowonronski)
-    *   [PR 1114](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1114) Exact option in PTC translation (J. S. Berg)             
+    *   [PR 1114](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1114) Exact option in PTC translation (J. S. Berg)
 *   Misc:
     *   [PR 1139](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1139) Remove inform forcing for seq generation (R. De Maria)
     *   [PR 1144](https://github.com/MethodicalAcceleratorDesign/MAD-X/pull/1144) Flag to throw fatal error when missing file in interactive mode (T. Persson)
@@ -33,7 +34,7 @@ MAD-X master
 
 MAD-X release 5.08.01 (2022.02.25)
 
-*    Fix radiation of mulitpoles in TRACK. #1079. (Riccardo)
+*    Fix radiation of multipoles in TRACK. #1079. (Riccardo)
 *    Fix in exact drift transfer map in TWISS #1077. (Tobias)
 *    Allows disabling the scaling of TWISS in PTC by default but possible to activate it. #1073. (Tobias)
 *    Fixes a bug in the dqmin calculation. #1075. (Tobias) 

VERSION

@@ -1,3 +1,3 @@
-VERSION      = 5.08.01
-VERSION_NUM  = 50801
-VERSION_DATE = 2022.02.25
+VERSION      = 5.09.00
+VERSION_NUM  = 50900
+VERSION_DATE = 2023.03.08

doc/latexuguide/pitfalls.tex

@@ -139,9 +139,14 @@ \chapter{\madx recipes and pitfalls}
   parameters via a transformation from the Teng-Edwards Twiss
   parameters. Obviously this fails when the Teng-Edwards Twiss
   parameters are ill defined. In this case one has to rely on
-  \hyperref[chap:ptc-twiss]{\texttt{PTC\_TWISS}}.    
-  \item[MAKETHIN might invalidate a sequence]
-   Several \madx commands such as {\texttt{ELIGN}}, {\texttt{EFCOMP}} won't work on sequences produced by {\texttt{MAKETHIN}}. In order to use such commands on thin sequences, it is advisable to save the sequence on a file, and then re-load it.
+  \hyperref[chap:ptc-twiss]{\texttt{PTC\_TWISS}}.
+
+\item[MAKETHIN might invalidate a sequence]
+   Several \madx commands such as {\texttt{EALIGN}}, {\texttt{EFCOMP}}
+   won't work on sequences produced by {\texttt{MAKETHIN}}. In order to
+   use such commands on thin sequences, it is advisable to save the
+   sequence on a file, and then re-load it.
+
 
 
 \end{description}

doc/latexuguide/taper.tex

@@ -20,8 +20,8 @@ \chapter{TAPER: strength tapering for energy variation}
   can cycle \texttt{ITERATE} times through the dipole map to converge to the 
   appropriate tapering value of the element corresponding to the energy
   loss through the adjusted bending magnet.  \\
-  For quadrupoles and sextupoles, a tracking through the map is
-  done only once for any positive non-zero value of
+  For quadrupoles, sextupoles, octupoles and multipoles, a tracking
+  through the map is done only once for any positive non-zero value of
   \texttt{ITERATE}. \\  
   If \texttt{ITERATE=0} is specified, the tapering is only calculated
   from the momentum deviation at the entrance of the element and no
@@ -73,7 +73,8 @@ \chapter{TAPER: strength tapering for energy variation}
 
 The \texttt{KTAP} attribute can be calculated and stored, or applied
 from input file, for each node of type \texttt{BEND},
-\texttt{QUADRUPOLE} or \texttt{SEXTUPOLE} in the active sequence. 
+\texttt{QUADRUPOLE}, \texttt{SEXTUPOLE}, \texttt{OCTUPOLE} or
+\texttt{MULTIPOLE} in the active sequence.  
 
 If \texttt{STEPSIZE} is not specified or specfied as zero, the
 calculated value is stored and each node has a unique value 

madx/__init__.py

@@ -0,0 +1,3 @@
+from .madx import Madx
+
+print("Work in progress...")

madx/madx.py

@@ -0,0 +1,90 @@
+import cpymad.madx
+import re
+import numpy as np
+import scipy.optimize
+from collections import ChainMap
+
+from .tablemixin import TableMixIn
+
+
+class Table(cpymad.madx.Table, TableMixIn):
+    pass
+
+
+class TableMap(cpymad.madx.TableMap):
+    def __getitem__(self, name):
+        try:
+            t = Table(name, self._libmadx)
+            t.set_index("name")
+            return t
+        except ValueError:
+            raise KeyError("Table not found {!r}".format(name)) from None
+
+
+class Madx(cpymad.madx.Madx):
+    def __init__(
+        self,
+        libmadx=None,
+        command_log=None,
+        stdout=None,
+        history=None,
+        prompt=None,
+        **Popen_args,
+    ):
+        super().__init__(libmadx, command_log, stdout, history, prompt, **Popen_args)
+
+        self.table = TableMap(self._libmadx)
+
+    def track_single_4d(mad, coord=np.zeros(4)):
+        mad.track(onepass=True)
+        mad.start(
+            x=coord[0], px=coord[1], y=coord[2], py=coord[3])
+        mad.run()
+        mad.endtrack()
+        t = mad.table["track.obs0001.p0001"]
+        out = [t[n][-1] for n in ["x", "px", "y", "py"]]
+        return np.array(out)
+
+    def track_single_6d(mad, coord=np.zeros(6)):
+        mad.track(onepass=True)
+        mad.start(
+            x=coord[0], px=coord[1], y=coord[2], py=coord[3], t=coord[4], pt=coord[5])
+        mad.run()
+        mad.endtrack()
+        t = mad.table["track.obs0001.p0001"]
+        out = [t[n][-1] for n in ["x", "px", "y", "py", "t", "pt"]]
+        return np.array(out)
+
+
+    def find_closed_orbit_track_4d(mad, x0=None, xtol=1e-10,mode='4d'):
+        def ftosolve(x):
+            if abs(x).max()>1:
+                raise ValueError(f"{x} too large")
+            return mad.track_single(x) - x
+
+        res = scipy.optimize.fsolve(ftosolve, x0)
+        return res
+
+    def track_single_6d(mad, coord):
+        mad.track(onepass=True)
+        mad.start(
+            x=coord[0], px=coord[1], y=coord[2], py=coord[3], t=coord[4], pt=coord[5]
+        )
+        mad.run()
+        mad.endtrack()
+        t = mad.table["track.obs0001.p0001"]
+        out = [t[n][-1] for n in ["x", "px", "y", "py", "t", "pt"]]
+        return np.array(out)
+
+    def find_closed_orbit_track_4d(mad, x0=np.zeros(6), xtol=1e-10):
+        def ftosolve(x):
+            if abs(x).max()>1:
+                raise ValueError(f"{x} too large")
+            return mad.track_single(x) - x
+
+        res = scipy.optimize.fsolve(ftosolve, x0)
+        return res
+
+    def find_closed_orbit_twiss(mad):
+        t = mad.twiss()
+        return np.array([t[n][0] for n in ["x", "px", "y", "py", "t", "pt"]])

release_python.sh

@@ -0,0 +1,20 @@
+#!/bin/bash
+set -euo pipefail; IFS=$'\n\t'
+
+NAME=$( python setup.py --name )
+VER=$( python setup.py --version )
+
+echo "========================================================================"
+echo "Tagging $NAME v$VER"
+echo "========================================================================"
+
+#git tag v$VER
+#git push origin v$VER
+
+echo "========================================================================"
+echo "Releasing $NAME v$VER on PyPI"
+echo "========================================================================"
+
+python setup.py sdist
+twine upload dist/*
+rm -r dist/ *.egg-info

setup.py

@@ -0,0 +1,12 @@
+from setuptools import setup, find_packages
+
+setup(
+        name='madx',
+        version='0.0.0',
+        description='Methodical Accelerator Design MAD-X',
+        author='CERN',
+        author_email='riccardo.de.maria@cern.ch',
+        url='https://github.com/MethodicalAcceleratorDesign/mad-x',
+        packages=find_packages(),
+        install_requires=['numpy','matplotlib','cpymad'],
+)

src/emit.f90

@@ -305,9 +305,9 @@ subroutine emdamp(code, deltap, em1, em2, orb1, orb2, re)
         ktap = node_value('ktap ')
         sks = zero
         if (sk0.ne.0)  then
-           h = sk0 * (1 + ktap) ! tapering
+           h = sk0 * (one + ktap) ! tapering
         else
-           h = an / el * (1 + ktap) ! tapering
+           h = an / el * (one + ktap) ! tapering
         endif
 
         !---- Refer orbit and eigenvectors to magnet midplane.
@@ -489,19 +489,20 @@ subroutine emdamp(code, deltap, em1, em2, orb1, orb2, re)
         sk2 = zero
         sk3 = zero
         sksol = zero
+        ktap = node_value('ktap ')
         select case (code)
         case (code_quadrupole)  !---- Quadrupole
-           sk1 = bvk * node_value('k1 ') * (1 + node_value('ktap '))
+           sk1 = bvk * node_value('k1 ') * (one + ktap) ! tapering
            str  = sk1
            n    = 1
            twon = two
         case (code_sextupole)   !---- Sextupole
-           sk2 = bvk * node_value('k2 ') * (1 + node_value('ktap '))
+           sk2 = bvk * node_value('k2 ') * (one + ktap) ! tapering
            str  = sk2 / two
            n    = 2
            twon = four
         case (code_octupole)   !---- Octupole
-           sk3 = bvk * node_value('k3 ')
+           sk3 = bvk * node_value('k3 ') * (one + ktap) ! tapering
            str  = sk3 / six
            n    = 3
            twon = six
@@ -617,11 +618,13 @@ subroutine emdamp(code, deltap, em1, em2, orb1, orb2, re)
         an = node_value('angle ')
         if (an .ne. 0) f_errors(0) = f_errors(0) + normal(0) - an
 
+        ktap = node_value('ktap ')
+
         !---- Other components and errors.
         nord = 0
         do i = 0, max(nn, ns, n_ferr/2-1)
-           f_errors(2*i)   = bvk * (normal(i) + f_errors(2*i))   / (one + deltap)
-           f_errors(2*i+1) = bvk * (skew(i)   + f_errors(2*i+1)) / (one + deltap)
+           f_errors(2*i)   = bvk * (normal(i) * (one + ktap) + f_errors(2*i))   / (one + deltap) ! tapering
+           f_errors(2*i+1) = bvk * (skew(i)   * (one + ktap) + f_errors(2*i+1)) / (one + deltap) ! tapering
            ! get the maximum effective order; loop runs over maximum of user given values
            if (f_errors(2*i) .ne. zero .or. f_errors(2*i+1) .ne. zero)  nord = i
         enddo

src/mad_def.h

@@ -28,7 +28,7 @@
 #define MAX_TAG 50          /* for SXF output */
 #define CHAR_BUFF_SIZE 100000 /* size of dynamic char_buff members */
 #define IN_BUFF_SIZE 500000 /* initial size of buffer for command groups */
-#define LINE_FILL 240        /* max. line length -2 for "save" output */
+#define LINE_FILL 2400       /* max. line length -2 for "save" output */
 #define LINE_F_MAD8 70      /* the same, for mad-8 format */
 #define MADX_LINE_MAX 78         /* for SXF output */
 #define MATCH_WORK 10       /* no. of work spaces in matching */

src/mad_dict.c

@@ -1834,6 +1834,7 @@ const char *const_element_def =
 "tilt     = [r, 0],  "
 "k3       = [r, 0],  "
 "k3s      = [r, 0],  "
+"ktap     = [r, 0],  "  
 "mech_sep = [r, 0], "
 "v_pos = [r, 0], "
 "knl      = [r, {0}], "
@@ -1883,6 +1884,7 @@ const char *const_element_def =
 "lrad     = [r, 0], "
 "knl      = [r, {0}], "
 "ksl      = [r, {0}], "
+"ktap     = [r, 0], "  
 "angle    = [r, 0], "
 "apertype = [s, circle, circle], "
 "aperture = [r, {0}],  "

src/mad_node.c

@@ -406,17 +406,16 @@ void set_closed_orb_node(int *index, double *pos){
 }
 
 
-void get_tt_multipoles(int *nn, double *knl, int *ns, double *ksl){
+void get_tt_multipoles(int *nn, double *knl, int *ns, double *ksl, double ktap){
     nn[0]=current_node->p_elem->multip->nn;
     ns[0]=current_node->p_elem->multip->ns;
     for(int i=0;i<*nn;i++){
-      knl[i] = current_node->p_elem->multip->knl[i];
+      knl[i] = current_node->p_elem->multip->knl[i]* (1 + ktap);
     }
     for(int i=0;i<*ns;i++){
-      ksl[i] = current_node->p_elem->multip->ksl[i];
+      ksl[i] = current_node->p_elem->multip->ksl[i] * (1 + ktap);
     }
 
-
 }
 void alloc_tt_attrib(int *length){
   current_node->p_elem->tt_attrib = mycalloc("tmp_array_tt", (*length+1), sizeof (*current_node->p_elem->tt_attrib));
@@ -538,13 +537,10 @@ store_node_value(const char* par, double* value)
   else if (strcmp(lpar, "k2") == 0) store_comm_par_value("k2",*value,el->def);
   // The inform is to make sure they are written out to a new sequence. 
   else if (strcmp(lpar, "ktap") == 0) {
-    store_comm_par_value("ktap",*value,el->def);
-//    el->def->par_names->inform[9] = 1;  // for quads and sextupoles
-//    el->def->par_names->inform[23] = 1; // for bends
+    set_command_par_value("ktap",el->def, *value); //mark it set
   }
   else if (strcmp(lpar, "lagtap") == 0) {
-    store_comm_par_value("lagtap",*value,el->def);
-//    el->def->par_names->inform[9] = 1;
+      set_command_par_value("lagtap",el->def, *value); //mark it set
   }
   else if (strcmp(lpar, "lag") == 0) store_comm_par_value("lag",*value,el->def);
 

src/mad_node.h

@@ -134,7 +134,7 @@ void    alloc_tt_attrib(int *length);
 void    set_tt_attrib(int *index, double *value);
 double  get_tt_attrib(int *index);
 void    set_tt_multipoles(int *maxmul);
-void    get_tt_multipoles(int *nn, double *knl, int *ns, double *ksl);
+void    get_tt_multipoles(int *nn, double *knl, int *ns, double *ksl, double ktap);
 double  node_obs_point(void);
 void    store_orbit_correctors(void);
 double  get_closed_orb_node(int *index);

src/taper.f90

@@ -100,7 +100,7 @@ subroutine taperreset(error)
 
   select case (code)
 
-  case (code_rbend, code_sbend, code_quadrupole, code_sextupole)
+  case (code_rbend, code_sbend, code_quadrupole, code_sextupole, code_octupole, code_multipole)
      call store_node_value('ktap ', 0.d0)
 
   end select