Fix bug wiggler rounding

[oscarxblanco]

Dear all,
this PR solves issue #1012 where there was in inconsistent behaviour of a wiggler set to Bmax=0 when compared to a drift of the same length.

The number of iterations is now calculated using round. The user must take into account that the ratio between total wiggler length and the wiggle_period should be an integer.

In addition, wiggler elements cannot be tracked alone. I include a check on the definition of gamma that is only set by the lattice, according to pull request #894 .

[oscarxblanco]

@lfarv and @swhite2401 ,

atError is failing when used inside GWigSymplecticPass and GWigSymplecticRadPass methods. I have included it to avoid tracking a wiggler without energy because it is used latter in the integrator. For example, here

at/atintegrators/gwig.c

Line 157 in 780657c

double VCw = pWig->VCw_raw[i] * pWig->Aw / pWig->Po;

where pWig->Po has been previously set to:

at/atintegrators/gwig.c

Line 40 in 780657c

Wig->Po = sqrt(gamma*gamma - 1.0);

Do you prefer that I remove the energy check ? Or to solve this error ?

atintegrators/atelem.c:205:22: note: expanded from macro 'atError'
        205 | #define atError(...) return (struct elem *) PyErr_Format(PyExc_ValueError, __VA_ARGS__)

[oscarxblanco]

Dear @SebastienJoly ,
this should fix the issue #1012 that you reported. Could you check it ?

Wrt to the tracking of a single wiggler, I failed to understand how to modify atError to run macOS, it runs on linux but I can not test macOS locally, so, I will open a separate issue.

[lfarv]

@oscarxblanco: normally, to track a single element, you can use the energy keyword in element_track() or Element.track() to provide the necessary value.

However, it is useful to understand why atError fails. At compile time ? This macro must run in both python and Matlab, I can have a look if you open another issue for that.

[SebastienJoly]

@oscarxblanco I am afraid the problem remains. If I use the example in the issue, I get different tunes if the length of the last drift/undulator is equal to 0.425. But, I do obtain the exact same tunes if it is 0.4, 0.41, 0.42, 0.43, or even 0.435. Does it work with you?

[oscarxblanco]

@oscarxblanco: normally, to track a single element, you can use the energy keyword in element_track() or Element.track() to provide the necessary value.

However, it is useful to understand why atError fails. At compile time ? This macro must run in both python and Matlab, I can have a look if you open another issue for that.

Oh , great ! Thank you, it works.

>>> und.track(zin,energy=1e9)
array([[4.250e-06],
       [1.000e-05],
       [0.000e+00],
       [0.000e+00],
       [0.000e+00],
       [2.125e-11]])

[lfarv]

@oscarxblanco, @SebastienJoly, @swhite2401 : I think that we should be more strict. If the wiggler length is not a multiple of the period length, it is an error. With the rounding set here, there is still a problem: when tracking n_periods x period_length, you still do not reach the total length !

I would prefer raising an exception if total_length / period_length is not an integer, with a very small tolerance.

[oscarxblanco]

@oscarxblanco I am afraid the problem remains. If I use the example in the issue, I get different tunes if the length of the last drift/undulator is equal to 0.425. But, I do obtain the exact same tunes if it is 0.4, 0.41, 0.42, 0.43, or even 0.435. Does it work with you?

Yes, it works for me. Here is my script, I get the same tune with a Drift of or a Wiggler.
example2.py

[SebastienJoly]

I am confused, I confirm that the tune difference still appear in your example. Moreover, I get a difference of the order of 1e-8 and 1e-13 for the first and last element of the list obtained from:
print(zout2.T-np.ravel(zout).T)
The difference drops to 1e-20 when I use an undulator length which does not lead to a tune difference.

[oscarxblanco]

@SebastienJoly , could you try again ? As requested by @lfarv I have added a check of the wiggler total length to period length ratio. It works on pyat, not yet in AT matlab.

When I set the total length to 0.425 m, I get a difference in tune between the wiggler set Bmax=0 and a drift on the order of 10^-14, and a difference of 10^-22 or less on the tracking of a particle starting with an angle of 10e-6 rad.
I leave here below the exact values

tune difference : [ 3.38618023e-15 -1.83741911e-14             nan]
track difference : [[-8.47032947e-22  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00 -2.94036858e-25]]

[oscarxblanco]

pyat implements the Ltot/Lw check at the python level so the check written in c on the pass method is redundant.

For AT matlab I tested it and It works, and I think it is worth leaving the check at the c level because there are many cases when the lattice is saved as .mat and elements are used with any configuration they had a the time of creation of the file.

[oscarxblanco]

Hi @SebastienJoly ,
I tried again to compare the tune difference of a varying length undulator with Bmax=0 and a drift of the same length.
I see no issue now.
The tune difference changes a bit, but, it is on the order of 10^{-14}.
I am attaching a screenshot of the tune difference vs length calculated in matlab on 10 points and in pyat on 1000 points.
I also attach the script I use to produce the pyat plot (on the right).
And I leave here below the versions I am running.

Python 3.11.2 (main, Apr 28 2025, 14:11:48) [GCC 12.2.0]
at.__version__=0.7.2.dev56+g59a25cbd.d20251107
np.__version__=1.26.0

plot_result_example2.py

[oscarxblanco]

Dear @SebastienJoly and @lfarv ,
I have removed several changes that are not really needed. The only significant changes are:

• Pass methods now use the c function round to determine the number Nw.
• The tolerance in the atwiggler LENGTH/LW ratio is set to 10^{-9}.

Now I'll wait for @SebastienJoly to confirm what I get with this modification before trying anything else.

I think @SebastienJoly mentioned in #1012 to set LENGTH and NW directly in order to avoid problems with the lengths LENGTH/LW ratio. Unfortunately I think that would break compatibility or make the element a little too complicated (as is the case of the VariableThinMultipole, for example) because one would need to assure that redundant parameters are given the right order of priority and privilege one. But, in any case, if you consider this necessary it could be implemented.

[lfarv]

@oscarxblanco : If I understand, there is no more test at the C level, we rely on the tests at the creation of the element. That's enough for me, it just leaves the case of "wrong" wigglers store in .mat files…

@SebastienJoly : you are right, defining the wiggler by its number of periods rather than the period length would eliminate the problem. But I agree with @oscarxblanco that it would either break the compatibility or make the element creation much more complicated. I can add that the present definition is based on the fact that in real life, a wiggler is usually defined by its period length ( U35, W240…) rather than by its number of periods…

I can approve the PR, but I let @SebastienJoly decide !

[lfarv]

OK

[SebastienJoly]

Dear @oscarxblanco @lfarv,
Thanks for the additional checks and for convincing me this is a minor problem unlikely to impact my simulations when b_max != 0.
However, the tune difference remains, even in the example2.py you sent above. I asked a colleague to try running it using this branch, and he observed the same tune difference. I have exactly the same pyat and numpy versions as you.
Could it be a problem related to the machine precision then?

[oscarxblanco]

Hi @SebastienJoly ,
could you give more details on the difference you see ? Is it on the plot ? or is it a particular case you tried again ?
I wonder if it is the installation or still the pass method.

• If you have matlab, could you do atmexall and run an example that is failing ? That will point to the pass method.
• About the pyat installation, I usually set a virtual environment with python -m venv some_name, I source the activation file, move to my at code copy, do rm -rf build and do pip install --config-settings mpi=1 ".[mpi]". Do you do anything different ? I use Debian, which OS and architecture are you using ?

gcc (Debian 12.2.0-14+deb12u1) 12.2.0
Linux pc542 6.1.0-40-amd64 #1 SMP PREEMPT_DYNAMIC Debian 6.1.153-1 (2025-09-20) x86_64 GNU/Linux

In case it is the pass method,

• could you add a printf to check the number of iterations ? I would like to know if you still have the cases when the integration is missing some steps. This is the code now:
  

  at/atintegrators/GWigSymplecticPass.c

  Lines 51 to 59 in 3dc3686

  	static void GWigPass_2nd(struct gwig *pWig, double *X)
  	{
  	int Nstep = pWig->PN*(pWig->Nw);
  	double dl = pWig->Lw/(pWig->PN);
  	for (int i = 1; i <= Nstep; i++) {
  	GWigMap_2nd(pWig, X, dl);
  	}
  	}

  
  And I suggest you to add one line:

static void GWigPass_2nd(struct gwig *pWig, double *X)
{                                   
  int Nstep = pWig->PN*(pWig->Nw);
  double dl    = pWig->Lw/(pWig->PN);

  printf("Nstep=%d, dl=%.6g,PN=%d,Nw=%d,Lw=%.3g\n",Nstep,dl,pWig->PN,pWig->Nw,pWig->Lw);
                                              
  for (int i = 1; i <= Nstep; i++) {
    GWigMap_2nd(pWig, X, dl);
  }
}

You should see many lines, like this, confirming the parameters in the call to the integrator :

Nstep=1000, dl=0.00045,PN=10,Nw=100,Lw=0.0045

[SebastienJoly]

Hi @oscarxblanco
I followed your recommendation by creating a fresh python environment and reran the following script:

import numpy as np
import at

QF = at.Multipole("QF", 0.5, poly_b=np.array([0,1.2]), poly_a=np.zeros(2))

drift_length = 0.5
Dr = at.Drift("Dr", drift_length)
HalfDr = at.Drift("Dr2", drift_length/2)
Dr3 = at.Drift("Dr3", 0.075)
Dr4 = at.Drift("Dr4", 0.425)
und = at.Wiggler('und', length=Dr4.Length, wiggle_period=Dr4.Length/100, b_max=0,Nmeth=2,Nstep=10,energy=1e9)

QD = at.Quadrupole("QD", 0.5, -1.2)
Bend = at.Dipole("Bend", 1, 2 * np.pi / 20)

FODOcell = at.Lattice(
    [HalfDr, Bend, Dr, QF, Dr, Bend, Dr, QD, Dr3, Dr4],
    name="Simple FODO cell",
    energy=1e9,
)

print(FODOcell.radiation_parameters().tunes)

uFODOcell = FODOcell.deepcopy()
uFODOcell[-1] = und.deepcopy()

print(uFODOcell.radiation_parameters().tunes)

ring0 = at.Lattice([und],energy=1e9)
zin = np.zeros((6,1))
zin[1] = 1e-5
zout,*_ = ring0.track(zin)
zout2 = Dr4.track(zin)
print(zout2)
print(zout2.T-np.ravel(zout).T)

and here is the output I get

[0.32536863 0.2625984         nan]

[0.32498454 0.2623137         nan]

[[4.250e-06]
 [1.000e-05]
 [0.000e+00]
 [0.000e+00]
 [0.000e+00]
 [2.125e-11]]
[[4.250e-08 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.125e-13]]

To go one step further, this is what I get when I do print(und)

Wiggler:
       FamName: und
        Length: 0.425
    PassMethod: GWigSymplecticPass
          Bmax: 0.0
            Bx: []
            By: [[1.]
 [1.]
 [0.]
 [1.]
 [1.]
 [0.]]
            Lw: 0.00425
        NHharm: 1
        NVharm: 0
         Nmeth: 2
         Nstep: 10

I am on the git branch with your fix and use numpy 1.26.0
As you can see, I still observe the tune difference I reported in the issue. Everything is working, meaning I get no error from a PassMethod call but I do not get the expected result,

Unfortunately, I still do not have matlab. I will ask a colleague if he can run your tests.

Do you need more information?

[oscarxblanco]

Hi @SebastienJoly , there is still something we are doing differently. I get the correct tune, and a single particle tracking difference of 10^{-22}. See here below

tune with drift          [0.32536863 0.2625984         nan]
tune with wiggler Bmax=0 [0.32536863 0.2625984         nan]
[[4.250e-06]
 [1.000e-05]
 [0.000e+00]
 [0.000e+00]
 [0.000e+00]
 [2.125e-11]]
Particle tracking difference : [[-8.47032947e-22  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00 -2.94036858e-25]]

With respect to the output you show, the last line says there is an error of 10^{-8} in the single particle tracking. I copy it here below for clarity

[[4.250e-08 0.000e+00 0.000e+00 0.000e+00 0.000e+00 2.125e-13]]

I can reproduce your result if I use the master branch, where there is a mistake on the number of step during the integration. So, my suspicions are still about the installation or compilation of the pass method. Could you try to add one line to the pass method and check the output ? It is the line I wrote in a previous message.

printf("Nstep=%d, dl=%.6g,PN=%d,Nw=%d,Lw=%.3g\n",Nstep,dl,pWig->PN,pWig->Nw,pWig->Lw);

In that way we will know that the code you are editing is the one you execute, and we could see some of the integrator parameters. Also, you could also check that the round function is implemented in the gwig.c file.

Other than that, I have tested the installation with pip gpu flags, pip mpi flags, just pip inside my copy of the repo, and all of them give me the right tune.

[SebastienJoly]

Hi @oscarxblanco
I added the line you mentioned, but still cannot see an output. However, when opening the gwig.c file, I can see your change with the round() function.
I tried to remove the build/ folder and reinstall pyAT using pip, but the printf() is still not active.
You are right, there must be a problem with the c compilation,n but I do not understand what I am doing wrong.

For context, I have

gcc (Ubuntu 11.4.0-1ubuntu1~22.04.2) 11.4.0
Description:	Ubuntu 22.04.5 LTS

Can it be related to my gcc version?

[oscarxblanco]

gcc (Ubuntu 11.4.0-1ubuntu1~22.04.2) 11.4.0
Description:	Ubuntu 22.04.5 LTS

Can it be related to my gcc version?

Very unlikely.
I think you are installing and running pyat in different virtual environments, or, you are installing the default pip version of pyat so your local copy is never used.

are you using the dot ? do not use the -e version for the moment.

cd my_at_copy/at
pip install .

if you are working on a IDE, e.g. spyder, or jupyter, or conda, etc., you need to configure the IDE to use the virtual environment you prefer.

Have you tried on a terminal ?

(some_venv) $ deactivate # just needed if you have any venv activated by default
$ python -m venv new_venv
$ source new_venv/bin/activate
(new_venv) $ cd my_at_copy/at # move inside your copy
(new_venv) $ git pull
(new_venv) $ git checkout fix_bug_wiggler_rounding
(new_venv) $ pip install .
(new_venv) $ cd to_your_code_folder
(new_venv) $ python
(new_venv) >>> exec(open('the_code.py').read())

and that is all. I hope it is complete enough as I am doing it from memory.