Fix RadIntegrals in ringpara

[mashl]

dear all
I think the new RadIntegrals could solve the previous bugs especially in the dispersion issue

[lfarv]

Hello @mashl . RadIntegrals still crashes if a wiggler has more than 1 harmonic:
Operator '+' is not supported for operands of type 'function_handle'.
at line 39

[lfarv]

@TeresiaOlsson: could you share your Diamond-II lattice with wigglers? Ideally, we would need wigglers with more than 1 harmonic (it's a problem at the moment), and with horizontal dependence of the vertical field (Kx ~= 0). It would be helpful also to have the 5 integral contributions from a single wiggler, both from analytical formulae and from Tracy. That's a lot…

[mashl]

Hello @lfarv
Actually I wanted to ask you an example for a multi-harmonic wiggler then I saw your message to @TeresiaOlsson .
Meanwhile, I will start to work on the problem with the typical multi-harmonic wiggler but it would be better if we have a benchmark for validation of the results.

[TeresiaOlsson]

@lfarv I will check what we have. We have a list of integral contributions for most of our IDs that has been benchmarked between analytic, Elegant and Tracy, but I'm not sure we have an AT lattice with all of them in. I will check that and if not make one, but then it might take a couple of days.

[mashl]

Dear @lfarv .
I changed the way of definition of the magnetic field in RadIntegrals beside little change in ringpara and now it works for multi-harmonic wigglers. for the test, I checked the below wigglers on ILSF lattice:
% Horizontal wiggler
atwiggler('WIG', LWIG, LWIG/20,1, E,'GWigSymplecticPass','By',[1,1;0.8,0.2;0,0;1,2;1,2;0,0])

%Vertical wiggler
atwiggler('WIG', LWIG, LWIG/20,1, E,'GWigSymplecticPass','By',[],'Bx',[1,1;0.8,0.2;1,2;0,0;1,2;0,0])

% Elliptical Polarized Wiggler
atwiggler('WIG', LWIG, LWIG/20,1, E,'GWigSymplecticPass','By',[1,1;0.8,0.2;0,0;1,2;1,2;0,0],'Bx',[1,1;0.8,0.2;1,2;0,0;1,2;0,0])

the results make sense.
I change the functions on my fork, I'm not sure if they updated in this pull request or not.
https://github.com/mashl/at/tree/master/atmat/atphysics/ParameterSummaryFunctions

if the results pass the benchmarks, for the next step, I think the dispersion could be calculated through the ID with a more efficient way than FDM with periodic boundary condition

[lfarv]

@mashl : the pull request is based on your fork, so all modifications you do are taken into account !

I see a problem with the dispersion: you cannot carry the dispersion over a non-linear element as you are doing. The field expansion is valid next to the beam axis, while the dispersion can be tens of centimetres or even meters.
So to compute the dispersion along the wiggler, you need to track 2 particles. The first, let's take oplus starts from (d*ds0, d*dsp0, 0, 0, d ,0), the second, say ominus, starts at (-d*ds0, -d*dsp0, 0, 0, -d ,0) where ds0, dsp0 is the dispersion at the entrance of the wiggler and d is a small number (atlinopt uses 1.0e-6). Then you get the dispersion as (oplus - ominus)/ d / 2.
If Bx(3,i) is zero, the wiggler is infinitely wide and acts horizontally as a drift space, so you get the right result, but if Bx(3,i) is non-zero, your dispersion is wrong. An easy check is to look at your computed dispersion at the last point (end of the wiggler) and compare it to the lattice value at the beginning of the following drift: they should match!

[mashl]

Hello @lfarv ,
I made some changes on gwig.c and get output from a particle while it passes through the wiggler. By using the outputs from oplus and ominus particles, I defined the dispersion as you mentioned. I tried to define the dispersion in general form and add the initial vertical dispersion to it too.
Oplus =[ dDis0(1); dDis0(2); dDis0(3); dDis0(4); d;0];
Ominus=[-dDis0(1); -dDis0(2); -dDis0(3); -dDis0(4); -d;0];

the modified codes are work well for ILSF lattice.

[TeresiaOlsson]

Hi @lfarv and @mashl ,
I have put together a file with the synchrotron integrals for the Diamond-II lattice that I get from Elegant, but I have problems with putting the wigglers into the lattice without breaking it. I will try to pull you latest changes and see if how far that gets me with solving them.

[TeresiaOlsson]

No, that made the situation worse... I am currently trying to use the wiggler_integrals branch, but now my lattice doesn't even work without the wigglers. There seem to have been differences made to the global parameters in atsummary when merging the centralized calculation of the radiation integrals that doesn't work for our lattice. Do you perhaps have a lattice file that works with these changes so I can check how our lattice needs to be modified?

[TeresiaOlsson]

atsummary in the wiggler_integrals branch doesn't work for the lattices in machine_data either. It gives the error:

"Unrecognized function or variable 'E0'.
Error in atsummary (line 34)
elseif isfield(E0,'LatticeFile')"

[mashl]

Hello @TeresiaOlsson ,
Please download my fork: https://github.com/mashl/at
and give me your email so I can send you our lattice.
my lattice which includes Wiggler and Cavity works with this fork

[TeresiaOlsson]

Thank you. It looks like your fork has the same problem in atsummary that doesn't work for our lattice or the ones in machine_data, but at least by looking at your lattice I could hopefully figure out why. My email is [email protected].

I also currently have a new problem in ringpara that didn't exist before. That function doesn't give errors, but I get the wrong momentum compaction which previously was okay.

[lfarv]

@TeresiaOlsson : can you send me your lattice? I'll open a new pull request to correct that bug.

[TeresiaOlsson]

@lfarv I have attached it. This is the version without wigglers since the other one needs a bit of a clean up for it to be easier to understand for someone outside of Diamond. I'm working on fixing that so it's easier to use for the benchmark.

I think there might be a difference in the lattice syntax between @mashl s fork and the common repository that caused the problem. Especially how global variables are used. We would prefer to keep the syntax as it was if possible since we have made that compatible with our MML and is currently starting to test AT2 with MML in the control room.

In the meantime, I got a working lattice file from @mashl and will try to temporarily modify our lattice to his syntax and test the radiation integrals with wigglers on his fork.

M_H6BA_15_1_1.txt

This is the values we had before:

************* Summary for '' ************
Energy: 3.50000 [GeV]
Gamma: 6849.32968
Circumference: 560.57386 [m]
Revolution time: 1869.87311 [ns] (0.53480 [MHz])
Betatron tune H: 0.16602 (88.78730 [kHz])
V: 0.23976 (128.22392 [kHz])
Momentum Compaction Factor: 1.17430e-04
Chromaticity H: +1.26236
V: +1.15736

---

Synchrotron Integral 1: 6.58344e-02 [m]
2: 3.17270e-01 [m^-1]
3: 1.72713e-02 [m^-2]
4: -1.18616e-01 [m^-1]
5: 3.81224e-06 [m^-1]
6: 2.37422e-02 [m^-1]
Damping Partition H: 1.37386
V: 1.00000
E: 1.62614
Radiation Loss: 670.29756 [keV]
Natural Energy Spread: 7.75757e-04
Natural Emittance: 1.57224e-10 [mrad]
Radiation Damping H: 14.19503 [ms] or 7591.44 turns
V: 19.50205 [ms] or 10429.61 turns
E: 11.99288 [ms] or 6413.74 turns
Slip factor: -1.17409e-04
Momentum compaction factor: 1.17441e-04 (1.17430e-04)

---

Assuming cavities Voltage: 1660.00000 [kV]
Frequency: 499.49913 [MHz]
Harmonic Number: 934
Overvoltage factor: 2.47651
Synchronous Phase: 2.72593 [rad] (156.18444 [deg])
Linear Energy Acceptance: 3.514 %
Synchrotron Tune: 0.00275 (1.47169 kHz or 363.39 turns)
Bunch Length: 2.95270 [mm] (9.84233 ps)

---

Injection:
H: beta = 08.369 [m] alpha = +2.1e-04 eta = +0.000 [m] eta' = +5.6e-07
V: beta = 04.397 [m] alpha = +2.0e-05 eta = +0.000 [m] eta' = +0.0e+00
********** End of Summary for '' **********

[lfarv]

The bug causing the crash is corrected and will appear in the master branch soon.
@TeresiaOlsson : Results are very close to your output, the difference is due to analytical integration used now. the momentum compaction is identical between atsummary and ringpara.

[TeresiaOlsson]

Thank you so much @lfarv . I will pull the change and see if I can get it to work now.

[lfarv]

The bug is now corrected in the master branch.
@TeresiaOlsson, you can keep your files as they are, can you check the results?
@mashl, can you merge the master into your branch so that you get the correction? Thanks

[TeresiaOlsson]

My lattice seems to work now :) I will continue with putting in the wigglers and see if can solve the problems I had with that from before. It might just be because the input to atwiggler has changed.

[mashl]

@lfarv , of course, I will update my branch.
did you check the changes in RadIntegrals?

[TeresiaOlsson]

There seems to be a bug with the I1 integral depending on where in the lattice the cavity sits, which caused the wrong value I saw for the momentum compaction in ringpara.

Cavity at the end:

Synchrotron Integral 1: 6.58525e-02 [m]
2: 3.17270e-01 [m^-1]
3: 1.72713e-02 [m^-2]
4: -1.18617e-01 [m^-1]
5: 3.81404e-06 [m^-1]
6: 2.37382e-02 [m^-1]

This is the expected result.

Cavity somewhere in the middle:

Synchrotron Integral 1: 7.11824e-02 [m]
2: 3.17270e-01 [m^-1]
3: 1.72713e-02 [m^-2]
4: -1.18585e-01 [m^-1]
5: 3.81377e-06 [m^-1]
6: 2.37289e-02 [m^-1]

Do you also get this for your lattices?

[lfarv]

@TeresiaOlsson: atsummary and ringpara need cavities and radiation to be turned off! Otherwise dispersion and beta functions are not defined. Then there should be no problem.

[TeresiaOlsson]

Ah, yes. I had forgotten that AT doesn't have checks for that like pyAT does. Thank you.

[TeresiaOlsson]

However, there is something that is inconsistent. ringpara does calculation of coupled damping times and for that it seems the cavity has to be on otherwise you get the following warning:

Warning: Matrix is singular to working precision.

In findorbit6 (line 101)
In findm66 (line 43)
In ringpara (line 106)

Warning: failed coupled damping times computation

In ringpara (line 108)

So at the moment ringpara gives a warning if you have the cavity off and the wrong I1 integral if you have it on and your cavity isn't in the beginning or end of the lattice. For both cases I have the radiation off so it's just the cavity causing the problem.

[TeresiaOlsson]

I have now managed to put in a single ID in my lattice using @mashl s fork and the results looks as below for I1-I5:

Elegant: 6.585240E-02 | 3.308244E-01 | 1.861466E-02 | -1.186174E-01 | 3.814208E-06

analytic: 6.58524E-02 | 3.30820E-01 | 1.86140E-02 | -1.18617E-01 | 3.81421E-06

ringpara: 6.7365389E-2 | 3.308248E-1 | 1.861472E-2 | -1.186031E-01 | 3.82E-06

atsummary: 6.58525E-2 | 3.17270E-01 | 1.72713E-02 | -1.18617e-01 | 3.82460e-06

It looks fairly good for ringpara except I1 which I think has a too large change. The effect of that is that the momentum compaction changes with 2.69884E-06 from 1.17473E-4 to 1.2017e-04, which is much more than I expected.

Could perhaps the cause be the same as for why I1 change depending on where the cavity sits in the lattice? I noticed that for the radiation effects of the wiggler to be included I need to use the passmethod "GWigSymplecticRadPass" which means the radiation is turned on for the wiggler even if it's turned of for the rest of the lattice.

Should I expect ringpara and atsummary to give the same output for this fork or only after it has been merged?

I will continue with putting in some other IDs that are stronger or sit in sections of high dispersion. This was just a fairly weak ID in a low dispersion section.

[lfarv]

@TeresiaOlsson: Thanks for those 1st results. But I insist that to get correct results, you should not include radiation anywhere: the radiation integrals rely on the dispersion function which is wrong if you include radiation. The calculation done in RadIntegrals.m does not (or should not: what do you think, @mashl) require the use of GWigSymplecticRadPass.

It looks like the numbers you give are for the whole ring, dipole radiation included. Can you get the ID contribution alone, which would make the comparison more accurate? And do you have a reference for what you use for analytical checks ? Thanks you !

[mashl]

Dear @lfarv , As you correctly mentioned the radiation effect does not and should not be included in the radiation integrals. The radiation integrals are depend on the magnetic field of the IDs and the optical functions of the lattice. So, the pass method of the wiggler should be GWigSymplecticPass.

[TeresiaOlsson]

Seems like I might have been doing something wrong the first time because now I also get an effect for GWigSymplecticPass. Sorry about that. Maybe it would be a good idea to add a check to make sure the cavity and radiation is off everywhere? I'm guessing I'm not going to be the only one that will be confused or not know that.

Yes, the numbers are for the whole ring. Do you mean I can get only the ID contribution from the code? Or you just want me to calculate the difference between with and without ID? These are the differences I get now with cavity off, radiation off and GWigSymplecticPass:

Elegant: -1.100000E-07 | 1.355490E-02 | 1.343320E-03 | 0.000000E+00 | 1.650000E-10
Analytic: -1.06322E-07 | 1.35506E-02 | 1.34269E-03 | -3.62211E-10 | 1.68091E-10

ringpara: 6.43562E-06 | 1.355524E-02 | 1.343380E-03 | 9.954206E-08 | 8.943734E-12

So it looks like it was the radiation that was the problem.

This is a reference I have for the analytic formulas (chapter 5) http://cds.cern.ch/record/399409/files/p807.pdf. I have a more detailed reference that I use myself, but I have to ask my colleague for permission before circulating it because it's his and part of a paper that's not published yet. I however think the only difference is how far the expansion of delta_I5 is done.

[TeresiaOlsson]

Now I have finished the whole lattice. I1 looks a bit too large to me, but I guess that depends on which precision we expect since it's overall small. Otherwise I only saw problems with I5 in straights with high dispersion or a strong ID so there seems to still be a problem with the dispersion.

An example for the K02 for the ID contribution to I1-I5:
Elegant: -8.000000E-08 | 1.019620E-02 | 1.010460E-03 | 0.000000E+00 | 2.897680E-07
ringpara: 4.73745E-06 | 1.01964E-02 | 1.01051E-03 | 4.62846E-08 | 2.46286E-05

I have attached the lattice together with a spreadsheet over all results. The spreadsheet has four pages, one with just the ID parameters and three with the results for elegant, AT and the analytic formulas.

To run the lattice you need to first run the script create_ID_settings where you can turn on individual IDs by setting their state to 1 instead of zero.

create_ID_settings.txt
M_H6BA_15_1_1_IDs.txt
Diamond-II IDs sept 2020.xlsx

[lfarv]

@TeresiaOlsson : thanks for the work to provide this large amount of data! We have a lot now to crosscheck our results.

I hoped that you could get individual contributions from Elegant, rather that making differences: the is less accurate (for instance on I1, you can get only one significant digit). But anyway with have enough data for benchmarking.

I agree with you that a safety check must be added in atlinopt for cavities and radiation on. The errors that you got are surprisingly large, and many people just ignore this… I'll open another branch for that.

[TeresiaOlsson]

@lfarv I'm not sure if you can get individual contributions from Elegant, but I will investigate. Maybe it at least could be possible to run a smaller part of the lattice and get more digits.

Thank you so much. I would definitely had forgotten to turn the cavity off and not be aware that the results were wrong unless you had pointed that out and I knew from another code what to expect.

Let me know if you run into problems with running the lattice. It's a bit complicated with all the IDs. It worked for me, but I could easily have made a typo somewhere so it doesn't work for you. There are also 3 three-pole wigglers in the lattice, but those I should have turned off.

[lfarv]

@mashl : I see a few problems with the files created in gwig.c

1. it makes GWigSymplecticPass unacceptably slow for usual tracking,
2. it crashes if the user has no write permission in the current working directory,
3. it leaves files all around.

2 and 3 can be solved by putting the file in /tmp or equivalent (see tmpfile, tmpnam, mkstemp in <stdio.h>).
1 can only be solved by duplicating most of GWigSymplecticPass.c to create a new function (not a passmethod) which outputs positions within the wiggler instead of just at the end. But duplicated code is worrying for maintenance… Otherwise, one can go back to another integrating method, or neglect the self-induced dispersion and use only the lattice dispersion.

Anyway, to avoid mixing problems, I suggest to leave this as it is for the moment, since it works, and focus on the wrong results for I1 and I5. Did you check the dispersion curve that you get through the wiggler?

[lfarv]

@TeresiaOlsson and @mashl : We have now enough data to validate the results for simple planar wigglers. Then we'll have to check:

• wigglers with kx ~= 0,
• multi-harmonic wigglers,
• vertical wigglers,
• elliptical or helicoidal wigglers.

I guess that analytical results cannot be used anymore (except probably for vertical wigglers). @TeresiaOlsson, can Elegant handle those cases? I would then try modifying one of your IDs (sorry for DLS users…).

[TeresiaOlsson]

@lfarv I will ask my colleagues that have more experience with wigglers. There are several different wiggler models in Elegant so they might have tested them and know how reliable they are.

[mashl]

@lfarv , I have some ideas for solving the gwig.c problems but we get back to them at the appropriate time.
Now, I'm working on I1 and I5 results.
the first problem I faced is the way of calculation of U0 in ringpara. it doesn't use the I2, so the effect of wiggler not included in it and all parameters depend on U0. can you check this?

[lfarv]

@mashl : you are absolutely right concerning U0. Since I2 is available, the energy loss/turn could be computed from it, and that would be similar to what is done in atsummary. Feel free to correct. In addition, atenergy (which is presently used and which is used in many other places) must be modified accordingly to include the wigglers. I take note for the future.

[mashl]

@lfarv , there were two major bugs! on RadIntegral I tried to solve them.
The debugging procedure is attached here. I think these new files should be used in the validation process.
@TeresiaOlsson , I fixed the compatibility issue between different versions of AT. I put my new lattice here.
RadIntegral Debug.pdf
SR_02_01.txt

[lfarv]

@mashl: your scan for step and d is convincing. So it looks like you have now a good dispersion. nstep=25 would make tracking 5 times slower than necessary, so this is again pushing for a dedicated function rather than using the integrator. For the moment, you should avoid modifying all wiggler elements in the ring by replacing at line 83 of RadIntegrals.m

ring{Wigidx(1),1}.Nstep=100;  %increase integration precision
ring{Wigidx(1),1}.Nmeth=4;
WIGELEM=ring(Wigidx(1));

by

WIGELEM=ring(Wigidx(1));
WIGELEM.Nstep=25;

But you can also see that the self-induced dispersion is in the tens of microns range (80 um in your case), which is completely negligible compared to the lattice dispersion. I'll prepare a much simpler computation neglecting this self induced dispersion and use only the lattice dispersion. It avoids the problem of tracking the dispersion, so that almost all integrals can be analytically computed.

[lfarv]

Following the outcome of @mashl study of the self-induced dispersion in wigglers, I tried another way of computing the wiggler radiation integrals in #188. Neglecting the self-induced dispersion makes a lot of simplifications. I ran the set of IDs provided by @TeresiaOlsson, and I get:

• an almost perfect agreement with Elegant
• a faster computation by a factor 22

@TeresiaOlsson and @mashl : can you have a look at this branch? Thanks

[mashl]

@lfarv , Sorry, I was busy today and now I'm far behind the discussion.
I would like to point out two things:
first, for computing dispersion within the wiggler, the particle only passes through the wiggler not whole the ring:
linepass(WIGELEM,Oplus,1:2);
therefore, I think a significant increase in computation time is more related to the writing process (fprint) in gwig.c rather than the number of integration steps.
the next point is, the agreement between the results of AT and ELEGANT, in the case we get the dispersion constant in the wiggler, maybe is because of the similarity of the computation process. I mean maybe the variation of dispersion is not taken to the account in radiation Integral calculation in ELEGANT.
So, I think it's better to compare the results of these two branch #186 and #188 to each other especially for the cases in which wiggler is located at high dispersion sections

[lfarv]

@mashl:
For speed: I fully agree, but be careful to modify the wiggler Nstep locally (WIGELEM variable) but not in the ring.
For comparison: I also agree, since I do not know what's inside Elegant, though I would trust it…

[TeresiaOlsson]

@mashl and @lfarv: At Diamond we have been discussing whether if Elegant includes the dispersion correctly or not and our conclusion so far is that it does. That's mostly based on that we got good agreement between the analytic formula and Tracy where Johan Bengtsson did the implementation including the dispersion. He however used some feature of Tracy that as far as I understand doesn't currently exist in AT and is a lot of work to implement.

The elegant user manual for the wiggler element says: "The radiation integrals were computed analytically using Mathematica, including the variation of the horizontal beta function and dispersion. For an odd number of poles, half-strength end-poles are assumed in order to match the dispersion of the wiggler. For an even number of poles, half-length end poles are assumed (i.e., we start and end in the middle of a pole), for the same reason." That indicate that it is taken into account, but the manual is very vague on how it's done.

So I agree, we aren't completely sure what Elegant is doing and it could be that we are fooling ourselves by trusting it too much. I will try to look at the Elegant source code and see if I manage to find out how it's done so we can feel more confident about the results.

[mashl]

@lfarv , @TeresiaOlsson ,
If we could take a pause for some days maybe I could provide another reliable benchmark.
actually, I used the Ohmi envelope method to study the effect of IDs on lattice parameters (it was part of my Ph.D. thesis), and of course, write its diffusion matrix calculator code. But the problem is the codes are not compatible with AT 2. I need some time to update their syntax. I compared the results of Radiation integrals and Ohmi envelope for four different lattices (ILSF, ALS, ALBA, NSLS II) and the results are in good agreement with each other.

[TeresiaOlsson]

@mashl It sounds really good to me if you have time to look into that. I eventually want to be able to do collective effects simulations including wigglers using atfastring and if I understood correctly from @lfarv that function uses the Ohmi envelop method so perhaps that would also be a help towards that?

[mashl]

@TeresiaOlsson , I'm not familiar with atfastring but if it uses the Ohmi envelop then this modification could be helpful for you. I used the ohmienvelope function on atx.