Bug in Lattice.avlinopt

[lfarv]

The results of atavedata (Matlab) and Lattice.avlinopt (python) are different. Here are the <β_x> for the first 10 elements of the test file hmba.mat.

Matlab	python
6.9000	6.9000
6.9000	6.9000
6.9000	6.9000
7.2396	7.2132
7.9188	.9188
7.9353	7.8489
7.4451	7.4451
5.8628	7.3854
4.9281	4.9281
3.6244	4.5013
------	-------

The same appears for <β_z>. From a 1st analysis, Matlab results are correct and it seems that the problem is in the average value in drifts in python.

[ZeusMarti]

Thanks for pointing it out, sorry I did not compare numerically the results, just qualitatively. I'll look at it.

[ZeusMarti]

Hi @lfarv ,

I think the problem affects only the drifts (which actually I did not consider in my initial tests): given this definition:

https://github.com/atcollab/at/blob/cb61c7aa6f1696bb4a6e52399907db62f2e36aca/pyat/at/physics/linear.py#L1141-L1146 https://github.com/atcollab/at/blob/cb61c7aa6f1696bb4a6e52399907db62f2e36aca/pyat/at/physics/linear.py#L1235-L1235

the booleans that select the focusing long elements are including the drifts

https://github.com/atcollab/at/blob/cb61c7aa6f1696bb4a6e52399907db62f2e36aca/pyat/at/physics/linear.py#L1268-L1271

https://github.com/atcollab/at/blob/cb61c7aa6f1696bb4a6e52399907db62f2e36aca/pyat/at/physics/linear.py#L1278

https://github.com/atcollab/at/blob/cb61c7aa6f1696bb4a6e52399907db62f2e36aca/pyat/at/physics/linear.py#L1289

This could be solved by just using k=0.0 instead of numpy.finfo(numpy.float64).eps in line 1145. In this case the behavior is the same as in Matlab. Actually in Matlab, the equivalent of fff is foc: https://github.com/atcollab/at/blob/cb61c7aa6f1696bb4a6e52399907db62f2e36aca/atmat/atphysics/atavedata.m#L57-L58 which does not include any drift.

Then there is the question of why betafoc and betadrift give different results with K=eps. I checked that this is true in both Matlab and Python, so there may be some issue with the formula...

For the moment would you agree on just changing line 1145 as a quick fix or you want to wait until everything is more clear?

[lfarv]

Hello @ZeusMarti: I think the problem comes when k is very small. In the expression of betafoc, the 1st two terms become extremely large, but partially cancel, approaching to 2*β₀L when k tends to 0. It's a kind of Infinity - Infinity problem. This is numerically very unstable and should be avoided.

So I think the problem also appears also for quadrupoles with very small k (fortunately quite unusual).

It would probably be safer to return 0 from get_strenghif k is less that some ε: magnets with very small strength will be considered as drifts.

But I agree with returning 0 on line 1145: it should immediately solve the problem.

[ZeusMarti]

Yes, you are correct, betafoc tends to the right value and at some point for very small K values starts giving nonsense numbers:

I'm not sure what you do in this cases, should I create a new branch for this small fix?

[lfarv]

By the way, I discovered the problem by running pytest pyat/test_matlab. This test sequence runs by calling Matlab from python and comparing the results. It cannot be run automatically on GitHub because of license problems with the Matlab engine. So I periodically run it locally, which proves very useful.

When the problem is solved, I'll try to do the same test in the "Matlab tests" action: it works the other way by calling python from Matlab, but this is supported on GitHub and runs on each commit. This will improve the automatic test sequence.

[lfarv]

@ZeusMarti : yes, post a pull request: committing directly to the master branch is not allowed!

[lfarv]

@ZeusMarti : according to your plot, I would replace all K < 1.e-7 by 0 !