Allow setting "dp" in 6D linear optics

[lfarv]

Up to now, in linear optics computations based on atlinopt6, "dp" could be set only for 4D lattices (so-called "radiation OFF"), since in 6D the off-momentum is fixed by the RF frequency.

With this PR, setting "dp" in atlinopt6 is also allowed in 6D, by first making a copy of the lattice with a modified the RF frequency.

This applies in turn to some other functions already based on atlinopt6:

• atplot
• tunechrom
• atfittune
• atfitchrom

and will apply in future conversions to atlinopt6.

This feature makes the behaviour for 6D lattices ("radiation ON") fully identical to to 4D.

Warning 1

This offers some new convenience, may cost extra computations: for a single call to atlinopt6 (for instance), the time spent by:

[rdata,edata]=atlinopt6(ring,'dp',dp,...);

is identical to the one spent by:

ring2=atsetcavity(ring,'Frequency','nominal','dp',dp);
[rdata,edata]=atlinopt6(ring2,...);

But if several functions are called in a row with an identical "dp" value, the frequency change will be executed several times, leading to additional useless computations. Better setting explicitly the desired frequency once at the beginning.

Warning 2

• For "new" functions where "dp" is a keyword parameter (atlinopt6 for instance), the new feature is only applied if "dp" is explicitly provided. Otherwise the lattice is left unchanged, as it is the case now.
• For "old" functions where "dp" is a positional parameter, whenever possible it becomes optional and the same rule applies. This is the case for atplot, tunechrom, atfittune, atfitchrom. But it may happen when converting other functions to atlinopt6 that this is not possible, if the "dp" argument cannot be distinguished from the next one. In that case, specifying 'dp',[] or 'dp',NaN will achieve the same result: "dp" will be set to 0 in 4D or left unchanged in 6D. These special values are accepted is all functions.

Warning 3

This behaviour is backward compatible with the previous one, with a minor difference:

in 6D, calling atplot(ring, 0.0,srange) used to ignore the 0.0 value. Now, it will effectively recompute the RF frequency to achieve 0.0, even it is was the case before. You must skip the "dp" argument and call atplot(ring,srange) to leave ring unchanged (see above). It's similar for the other functions.

in 6D, you must skip the "dp" argument or set it to [] or NaN to leave the lattice unchanged. Otherwise, the specified value will be honoured, with additional computation, which anyway correct!

Control

This feature is provided by a wrapper function wrapper6D which is called in all the previously mentioned functions. This allows to control the new behaviour with a global option, to disable it, enable it or enable it with a warning. This option and its default value is still to be defined.

[swhite2401]

Since you are modifying this @simoneliuzzo pointed out an issue with ringpara leading to wrong results in the present of synchrotron motion. It is due to the fact the it still calls atlinopt that leads to significant errors in the radiation integral computation.

We may want to integrate linopt6 in ringpara such that it returns consistent results in all configuration. Looking into this I also found that atlinopt return a column vector and atlinopt6 a row which is compatibility issue.

[lfarv]

@swhite2401: ringpara is of course in the list of functions to convert. Not on top of the list, atx and atfastring are the next ones (more critical).

You are right for the output of atlinopt and atlinopt6 (not exactly in fact: it's the opposite, atlinopt returns a row and atlinopt6 a column). I don't know what to do. A column makes more sense (the ring is a column), but I do not want to change atlinopt (compatibility !). It does not make any difference when extracting the different components: for instance beta=cat(1,lindata.beta); works for both. And anyway, the input and output are different. But turning the atlinopt6 output into a row would be possible, I think.

[swhite2401]

Well I ran into this problem earlier comparing DipoleRadiation(ring,lindata) -> I had to take the transpose for the atlinopt6 output, I am sure this issue will show up in other places when migrating form atlinopt to atlinopt6... anyway this is just a side comment

[lfarv]

• the elemdata output argument of atlinopt6 is a row vector, as in atlinopt,
• When specifying "dp" for a 6D ring, a warning may be emitted, depending on the 'WarningDp6D' option. The goal is to notify the user that it's better to handle the RF frequency explicitly to avoid multiple copies.
  At the moment, the default value of this option is "true", but this may be changed.
  In any case, the "dp" value will be correctly taken into account.

>> [rdata,edata]=atlinopt6(ringrad);             % No warning
>> [rdata,edata]=atlinopt6(ringrad,'dp',0.01);   % A warning is issued
Warning:
Specifying "dp" for a 6D lattice creates a copy with a modified RF frequency.
For a better efficiency, handle the RF frequency beforehand,
or to avoid this warning, use "setoption('WarningDp6D',false)" 
>>
>> setoption('WarningDp6D',false)
>> [rdata,edata]=atlinopt6(ringrad,'dp',0.01);   % No warning
>>

[lfarv]

The WarningDp6D option can also be applied for a single function call. The overwrites the default behaviour:

>> [rdata,edata]=atlinopt6(ringrad);                           % No warning
>> [rdata,edata]=atlinopt6(ringrad,dp=0.01);                   % A warning is emitted
Warning:
Specifying "dp" for a 6D lattice creates a copy with a modified RF frequency.
For a better efficiency, handle the RF frequency beforehand,
or to avoid this warning, use "setoption('WarningDp6D',false)" 
>>
>> [rdata,edata]=atlinopt6(ringrad,dp=0.01,WarningDp6D=false); % No warning for this statement
>> tune=tunechrom(ringrad,0.01);                               % But it is still emitted here
Warning:
Specifying "dp" for a 6D lattice creates a copy with a modified RF frequency.
For a better efficiency, handle the RF frequency beforehand,
or to avoid this warning, use "setoption('WarningDp6D',false)" 
>> 
>> setoption(WarningDp6D=false);                               % Turn off the warnings
>>
>> [rdata,edata]=atlinopt6(ringrad,dp=0.01);                   % No warning here
>> tune=tunechrom(ringrad,0.01);                               % Nor here
>> 

[lfarv]

@simoneliuzzo and @carmignani: this improvement is now ready and answers some of you questions. Do you have any comment?

[simoneliuzzo]

Dear Laurent,
Thank you for this developments.
I will need some time to test.
I try to do it tomorrow.
best regards
Simone

[simoneliuzzo]

Dear Laurent,

I prepared this small script to see if the usual commands I am used to use would work within this branch.
I noted on the side my observations.
I run the script both with master and atlinopt6 branches, forcing radiation on or off.
The list of observations becomes a bit difficult, so I will try to summarize better tomorrow.

close all
clear

setoption('WarningDp6D',false)

load('/Users/liuzzo/beamdyn/OPTICS/EBS/S28F_all_BM/LATTICE/AT/BMandCANTING/S28F_all_BM.mat')
load('/Users/liuzzo/beamdyn/OPTICS/EBS/S28F_all_BM/LATTICE/AT/BMandCANTING/S28F_PA.mat','ARCA');

figure('Name','rad off, dpp =0.0'); atplot(ARCA); % <-many slices lattice layout, missing dpp comment (ok in master)
figure('Name','rad off, dpp =0.01'); atplot(ARCA,0.01); % <- large BPMs in layout and missing dpp missing (ok in master)

ring=atSetRingProperties(r);

[~, ring]= check_radiation(ring,0,'force');  % make sure radiation is OFF 

figure('Name','rad off, dpp =0.0'); atplot(ring);  % <- ok

figure('Name','rad off, dpp =0.01'); atplot(ring,0.01); % <- ok

figure('Name','rad on, dpp =0.01'); atplot(atradon(ring),0.01); % <-  missing dpp comment  ( dpp ignored in master)

figure('Name','atradoff, dpp =0.01'); atplot(atradoff(ring),0.01); % <- missing layout (ok in master)

figure('Name','atsetcavity radoff, dpp =0.01'); atplot(atsetcavity(ring,6.5e6,0,992),0.01); % <- missing layout ( dpp ignored in master)


%% 
[ldp,tdp,cdp]=atlinopt(r,0.01,1);

rc = atfitchrom(r,0.01,cdp,'SF\w*','SD\w*'); % <- ask to do nothing

[ldpc,tdpc,cdpc]=atlinopt(rc,0.01,1);
disp("initial = requested at 0.01: " + cdp)
disp("final at 0.01: " + cdpc) % <- small change even if should have done nothing.

rc = atfitchrom(r,0.01,cdp+[1 1],'SF\w*','SD\w*');% <- ask to change by 1 unit in both planes at 1%
rc = atfitchrom(rc,0.01,cdp+[1 1],'SF\w*','SD\w*');
rc = atfitchrom(rc,0.01,cdp+[1 1],'SF\w*','SD\w*');
rc = atfitchrom(rc,0.01,cdp+[1 1],'SF\w*','SD\w*');
rc = atfitchrom(rc,0.01,cdp+[1 1],'SF\w*','SD\w*');

[ldpc,tdpc,cdpc]=atlinopt(rc,0.01,1);
disp("required +1: " + (cdp+[1 1]))
disp("final+1: " + cdpc) % <- after several iterations, still some visible discrepancy among request and final value


%% tune
[ldp,tdp,cdp]=atlinopt(ring,0.01,1);

rc = atfittune(r,0.01,tdp,'QF1\w*','QD2\w*'); % <- ask to do nothing

[ldpc,tdpc,cdpc]=atlinopt(rc,0.01,1);
disp("initial = requested at 0.01: " + tdp)
disp("final at 0.01: " + tdpc); % <- small change even if should have done nothing.

rc = atfittune(ring ,0.01,tdp+[1 1]/10,'QF1\w*','QD2\w*');% <- ask to change by 1 unit in both planes at 1%
rc = atfittune(rc,0.01,tdp+[1 1]/10,'QF1\w*','QD2\w*');
rc = atfittune(rc,0.01,tdp+[1 1]/10,'QF1\w*','QD2\w*');
rc = atfittune(rc,0.01,tdp+[1 1]/10,'QF1\w*','QD2\w*');
rc = atfittune(rc,0.01,tdp+[1 1]/10,'QF1\w*','QD2\w*');

[ldpc,tdpc,cdpc]=atlinopt(rc,0.01,1);
disp("required +1: " + (tdp+[1 1]/10))
disp("final+1: " + tdpc) % <- after several iterations, still some visible discrepancy among request and final value

[lfarv]

@simoneliuzzo: Thanks for the test Simone, I'll take your remarks one by one, but at some point, I'll need your lattices. For now:

load('/Users/liuzzo/beamdyn/OPTICS/EBS/S28F_all_BM/LATTICE/AT/BMandCANTING/S28F_all_BM.mat')
load('/Users/liuzzo/beamdyn/OPTICS/EBS/S28F_all_BM/LATTICE/AT/BMandCANTING/S28F_PA.mat','ARCA');

figure('Name','rad off, dpp =0.0'); atplot(ARCA); % <-many slices lattice layout, missing dpp comment (ok in master)
figure('Name','rad off, dpp =0.01'); atplot(ARCA,0.01); % <- large BPMs in layout and missing dpp missing (ok in master)

From this, I see:

• The "missing dpp comment" suggests that the 'ARCA' lattice you are loading has radiation ON, unlike what you say in the legend. Indeed, for a 6D lattice with a given frequency, you cannot easily derive dp/p. Since dp varies all along the lattice, the dp value of the closed orbit at any point is not an indication. To get it, you have to look at the closed orbit, record the dp value at the initial point (or any other point), then set the nominal frequency, record the dp at the same point and make the difference. I do not want atplot to make this sequence, so atplot does not display dp for 6D lattice. One could display the "desired" value, in case it is provided (like in your 2nd line), but the real value may differ slightly. And what about your 1st line where it's not specified?
• The slicing of synoptic elements looks like a real problem, I'll look at that.

ring=atSetRingProperties(r);

What is r ? Is it coming from S28F_all_BM.mat? You should always use the x=load(...) form so that you control the variable names and do not create unpredictable variables which may erase existing ones.

[ldp,tdp,cdp]=atlinopt(r,0.01,1);

If r has radiation ON, like ARCA, this whole sequence using atlinopt (and not atlinopt6) with radiation ON is doubtful…
More generally, for the following, you are alternating between r and ring and I do not see the rationale.

Thanks for your tests, they are really very useful!

[simoneliuzzo]

Dear Laurent,

I copyed different parts of scripts. r = ring.
I will send you via email the lattice files.
atlinopt is to be discarded forever? (then may be better to remove in AT3.0) I need to make sure it works, cause I call it many times basically everywhere in my code. So, if this branch fixes the backward compatibility of the atlinopt dp input, I will have less work to make the old code work.

best regards
Simone

[simoneliuzzo]

Dear Laurent,
concerning the missing dp/p comment, I understand your point. The frequency on the other hand is just one value. This value could define the dp/p. Or may be the frequency difference from nominal could be shown in the 6D case?
May be showing if the lattice is 6D or 4D or none of the two, could also be a useful content of the atplot comment.
best regards
Simone

[lfarv]

atlinopt is to be discarded forever?

No ! atlinopt will stay as it is. It's the only way to get the Sagan / Rubin analysis of coupled rings. And for that reason, it will be limited to 4D for ever. I'd just like to add the check_radiation test (when you are ready)

[swhite2401]

atlinopt is to be discarded forever?

No ! atlinopt will stay as it is. It's the only way to get the Sagan / Rubin analysis of coupled rings. And for that reason, it will be limited to 4D for ever. I'd just like to add the check_radiation test (when you are ready)

Isn't this atlinopt4?

[lfarv]

Isn't this atlinopt4?

You are right, but atlinop4 has the new interface. we need to keep atlinopt for backward compatibility.

[carmignani]

Hello @lfarv,
I just tried to do some tests on this new branch and I'm having a problem at the first atplot. Simone does not have the same problem and I don't understand why.
In line 69 of xplot there is ring0 that is not defined and so atplot fails.

[lfarv]

Hello @carmignani : I was too fast in trying to correct the bug reported by @simoneliuzzo… Pull the last revision and it should be OK. Sorry…

[lfarv]

@simoneliuzzo:

The frequency on the other hand is just one value. This value could define the dp/p. Or may be the frequency difference from nominal could be shown in the 6D case?

There is another way to derive dp for a 6D lattice: compute dct from df=f-f_nominal (trivial), and then use findsyncorbit (with rad. OFF) to get dp. And this suggest a new function, which, given one of dp, dct or df will use findorbit4 or findsyncorbit to output the 3 parameters. At the cost of a single orbit computation. This could be used for display, as you suggest, but also to replace the use of the momentum compaction factor to compute df for 6D lattices. More accurate (it does not assume a linear α) and faster (1 orbit computation instead of 2 in mcf).

I propose to leave this branch as it is (no dp display in atplot for 6D), and to introduce this new function (atoffmomentum ?) in another PR.

[carmignani]

I'm slowly checking the different modifications of the pull request.
Do we need to keep the warning in case the ringparam is not there? It looks like it is a quite frequent warning, but I'm not sure it is nice to have a ringparam element all the time in the lattices, especially if we work with only single cells, during lattice design phase. Can't we remove the warning?

[swhite2401]

I'm slowly checking the different modifications of the pull request. Do we need to keep the warning in case the ringparam is not there? It looks like it is a quite frequent warning, but I'm not sure it is nice to have a ringparam element all the time in the lattices, especially if we work with only single cells, during lattice design phase. Can't we remove the warning?

I agree this warning is very annoying and has to be removed, this is detrimental to user's experience: warnings should be avoided in general they are counter-productive

[lfarv]

Do we need to keep the warning in case the ringparam is not there?

We can remove it. The idea comes from time measurements I did: to get a ring parameter, for a test lattice (hmba I think), it takes 80 ms without RingParam, and 0.2 ms with. Huge gain. Another example: when requiring ring properties in each test, the standard test sequence of GitHub (matlab_tests) takes on my Mac 15s with RingParam, and 45s without (it was a modified test sequence tailored for that measurement, however "hmba" has the RingParam, so the whole difference come from the "dba" tests). There is nothing new there, it has always been like that, I just noticed it.

Even when working on single cells, and may be even more, it's useful to create the RIngParam element (which by the way provides the Periodicity property useful in that case).

So I agree, I'll remove the warning, but if you have ideas to recall users that it's beneficial, that would be nice.

Forgot to say: the advantage of the warning is that it tells exactly what to do !

[swhite2401]

Do we need to keep the warning in case the ringparam is not there?

We can remove it. The idea comes from time measurements I did: to get a ring parameter, for a test lattice (hmba I think), it takes 80 ms without RingParam, and 0.2 ms with. Huge gain. Another example: when requiring ring properties in each test, the standard test sequence of GitHub (matlab_tests) takes on my Mac 15s with RingParam, and 45s without (it was a modified test sequence tailored for that measurement, however "hmba" has the RingParam, so the whole difference come from the "dba" tests). There is nothing new there, it has always been like that, I just noticed it.

Even when working on single cells, and may be even more, it's useful to create the RIngParam element (which by the way provides the Periodicity property useful in that case).

So I agree, I'll remove the warning, but if you have ideas to recall users that it's beneficial, that would be nice.

Forgot to say: the advantage of the warning is that it tells exactly what to do !

Can't we built an ring=atloadlattice(filename, matkey) function, similar to the one we have in python that would create this element? We may use it in instead of the native load() function

[lfarv]

Can't we built an ring=atloadlattice(filename, matkey) function, similar to the one we have in python that would create this element? We may use it in instead of the native load() function

That's the best solution. And then, try to convince users to use it instead of the usual load.

[simoneliuzzo]

Ok for the load functions. One line instead of 3-4 is better, and it is a single change, likely not perturbing the following code. I imagine that the loaded lattices would be defaulted to 4D? This could give troubles when loading optics tuned for off energy operation. May be RingParam could remember this feature.
To convince users It will have to work also for structures that are not rings (arc, transfer lines, etc..), and eventually load all AT lattices in the file (how to recognize them among other variables?).

[simoneliuzzo]

Concerning atlinopt dp input, I have no objection in having a warning that something is wrong, (it could even be further detailed with: please use atlinopt6). I do not understand why the wrong atlinopt+6D+dp output changed compared to previous versions of AT. The correct 6D optics are given by atlinopt6, why changing also atlinopt?
If atlinopt may not be recovered, then I would perfer an error to the warning, as the values output of atlinopt+6D+dp are largely wrong.

[simoneliuzzo]

Sorry, I did not read well the outputs of my test functions (thanks @carmignani ). The functions that give the warning and the wrong results are only atfittune atplot and atfitchrom. atlinopt is ok.

[carmignani]

This morning I did a git pull in this branch and I had a lot of conflicts in several files. I tried to solve the problems for about half hour, then I recloned the repository from zero. I thought it was only a problem of my local repository, but I see that Simone had a very similar problem just now, so maybe there is something wrong in the last commits. By the way, the new clone works well, so I don't understand what was the problem.

[lfarv]

@simone:

The functions that give the warning and the wrong results are only atfittune atplot and atfitchrom. atlinopt is ok.

Where do you get a warning ?
What result may be wrong ?

[lfarv]

@carmignani:

This morning I did a git pull in this branch and I had a lot of conflicts in several files

Did you modify any file in your working directory? I don't understand what happened. But I also do not understand the 17 commits from me reported above. They look like old things, I did not work so much yesterday ! I hope everything is back to normal.

[lfarv]

@simoneliuzzo:

Concerning atlinopt dp input, I have no objection in having a warning that something is wrong, (it could even be further detailed with: please use atlinopt6)

Do I understand that rather than adding right now a "check_radiation" in atlinopt to prevent using it in 6D, you would put (temporarily) a warning mentioning the problem and pushing to use atlinopt6 instead?

That looks like a good way to diagnose and correct this problem step by step. In addition, the warning could be silenced, like what is done for "dp in 6D".

[carmignani]

@lfarv,
I didn't modify anything in my local copy of the repository. I did a git pull and I got an error saying that I had something like 15 commits above and 17 behind so I had to solve the conflicts before merging.
One conflict was about the function atfindparams I think. I can copy here the errors if you want tomorrow, now I'm on a different computer. I really could not fix the problem and finally I recloned a new copy from zero. Simone had the same problem and he also recloned.

[lfarv]

@carmignani: don't worry, I could not interpret the messages, I'm not a git expert. I probably did something wrong, I don' know what, but a far as I could see, the repository is correct…