Matlab frequency control

[lfarv]

This applies to Matlab the modifications of #540: in 6D, the functions atlinopt6, findm66, findorbit6 and tunechrom now accept a non-zero dp, dct or df argument. The RF frequency will be temporarily modified to provide the desired off-momentum value.

[lfarv]

@simoneliuzzo : Any comment on this ?

[simoneliuzzo]

Dear @lfarv,

I did not look at this PR yet. I will do it soon.

[simoneliuzzo]

Dear @lfarv,

I think this will require quite some tests, and I do not yet know exactly which ones...

Could we add (at least in pyAT) a test (those run at "git push" to the repository) for each of the functions mentioned above making sure off-energy tracking is correct?

For example showing that the same results (parameters, optics, etc..) are produced for the:

• on-energy lattice when equivalent df, dp, dct are applied
• off-energy lattice (with a given df already set in the lattice) when equivalent df, dp, dct are applied

[lfarv]

@simoneliuzzo : I added a test checking that tunes, chromaticities and beta functions are identical in 4D and 6D when varying dp. Is it ok for you ?

[lfarv]

Rebased on master

[lfarv]

@simoneliuzzo : can you check that the implemented tests are good enough ? I intend to merge this week.

[simoneliuzzo]

Dear Laurent,

I will likely not have time to test this week. May be an other reviewer could do this in time?

best regards
Simone

[lfarv]

I need someone to approve this PR pending for more than 6 weeks now…

[swhite2401]

@lfarv, for me the principle is ok since we implemented it in python. However, I am not a matlab user so I am not able to validate this part... is this blocking you for other developments?

[lfarv]

@swhite2401 : the frequency_control function introduced here is necessary to go on on #585… I understand that you cannot spend too much time on testing and validating things… But I'd like to settle this before leaving for a few weeks!

[simoneliuzzo]

Dear @lfarv,

I am postponing testing this PR as It appears like a lot of work to do and I never feel like I have enough time in a row to do the job.

There are 4 major functions, with 3 inputs to be tested each for at least 2 different optics and at least for the on- energy and off- energy case. It makes 4 * 3 * 2 * 2 = 48 cases to test.

If all these cases could be included in the test then it would be easy to just verify that the tests are correct and approve the PR.

[lfarv]

@simoneliuzzo : I understand ! But the integrated test covers roughly what you mention, including 2 different off-momentum values (positive and negative), apart from the 2 different optics (I have only one with 6D capability: hmba).

[simoneliuzzo]

The test could be extended to also include off-energy lattices then?
The input lattices would have frequency not set to the nominal value. As for the EBS booster and SR off-energy operation.

Also I see that the test looks only at dp. What about ct and df?

[lfarv]

@simoneliuzzo : The test compares the 4d and 6d lattices for 3 dp values: -0.01, 0 and 0.01. But you are right, it does not take dct and df into account. The problem is choice of the range: the 3 dp values work for ~any lattices. For dct and df. the ranges vary a lot because of large differences of α between lattices. It would require indexing the df and dct values with the lattice name… I'll look at that !

[lfarv]

@simoneliuzzo : the test sequence now includes 6D orbit and optics tests varying dp, dct and df. 6D orbits are checked comparing the Matlab and python results on 2 different lattices, optics results are checked by comparing 4D and 6D results on single lattice. This represents 27 cases among your 42 theoretical cases, but the others (findm66) are implicitly included.

Is this good enough for you ?

[simoneliuzzo]

Dear @lfarv,

the updated methods would work also for off-energy optics?

We are recently frequently working with optics that are matched off-energy, with dp=1% set in the lattice with a frequency shift.

Would the results be the same in pyAT and matlab AT also in this case, where the input lattice has already a shift in RF frequency?

thank you
best regards
Simone

[lfarv]

@simoneliuzzo: This PR concerns the introduction of dp, dct and df as keyword parameters of several optics functions for 6D lattices: in this case, the RF frequency will be tuned to satisfy the requirements. This was recently provided in python in #540, and of course the results are identical in both versions. Up to know, you would get a warning "In 6D, dp is not taken into account...". But if you don't provide any of these arguments, the behaviour is unchanged: the lattice is taken as it is, without modifying the frequency. This is like this since the introduction of atlinopt6. With again identical results in Matlab and python. Note that this is different from specifying dp=0 !

[simoneliuzzo]

Dear @lfarv,

I cloned and checkout this branch. I will try to use the new inputs and give you feedback.

[simoneliuzzo]

Dear @lfarv,

I am running some tests. How to input dp, dct, df for findorbit6? I looked at the help but could not find out.
may be the help is not updated?

I have just cloned a new AT and switched to matlab_frequency_control branch.
I have also set my path to use this AT and run atmexall.

[simoneliuzzo]

Dear @lfarv,

I figured out the input for findoribt6. I think the help should be updated as for atlinopt6.

The warning issued when using dp, dct or df input always mentions dp.

[simoneliuzzo]

Dear @lfarv
I have tested findorbit6 with these inputs:

indrf = find(atgetcells(ring,'Frequency'))';
alpha = mcf(ring);
f0 = atgetfieldvalues(ring,indrf,'Frequency');
dp = 0.01;
df = -f0(1)*dp*alpha;
C = findspos(ring, length(ring)+1);
dct=-C*df/(f0+df);

I wanted to check that for equivalent dp, df and dct the output of findorbit is the same

    disp("dp = "+dp)
    orb6 = findorbit6(ring,1:10,'dp',dp);
    disp(orb6)

    disp("dct = "+dct)
    orb6 = findorbit6(ring,1:10,'dct',dct);
    disp(orb6)

    disp("df = "+ df +"Hz")
    orb6 = findorbit6(ring,1:10,'df',df);
    disp(orb6)

When I look at the printed data this is the output:

dp = 0.01
  Columns 1 through 8

   0.001477818094902   0.001479723275425   0.001479723275425   0.001480139561676   0.001480139561676   0.001479837883947   0.001479564318121   0.001230100305299
   0.000000630614112   0.000000630614112   0.000000630614112   0.000000630614112   0.000000630614112  -0.000002154515213  -0.000002154515213  -0.000516781169136
                   0                   0                   0                   0                   0                   0                   0                   0
                   0                   0                   0                   0                   0                   0                   0                   0
   0.010409572282057   0.010409572282057   0.010409572282057   0.010409572282057   0.010409572282057   0.010409572282057   0.010409572282057   0.010409571390395
  -0.092564275562642  -0.092564275562048  -0.092564275562048  -0.092564275561918  -0.092564275561918  -0.092564275561677  -0.092564275561386  -0.092564232258413

  Columns 9 through 10

   0.001181360999074   0.001181360999074
  -0.000516781169136  -0.000516781169136
                   0                   0
                   0                   0
   0.010409571390395   0.010409571390395
  -0.092564219794381  -0.092564219794381

dct = 0.0015008
  Columns 1 through 8

   0.001463697889157   0.001465578137679   0.001465578137679   0.001465988976231   0.001465988976231   0.001465694475254   0.001465426560187   0.001218346599528
   0.000000622307415   0.000000622307415   0.000000622307415   0.000000622307415   0.000000622307415  -0.000002109827884  -0.000002109827884  -0.000511820432195
                   0                   0                   0                   0                   0                   0                   0                   0
                   0                   0                   0                   0                   0                   0                   0                   0
   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321557358685
  -0.092530461362482  -0.092530461361903  -0.092530461361903  -0.092530461361777  -0.092530461361777  -0.092530461361546  -0.092530461361267  -0.092530418880150

  Columns 9 through 10

   0.001170070951314   0.001170070951314
  -0.000511820432195  -0.000511820432195
                   0                   0
                   0                   0
   0.010321557358685   0.010321557358685
  -0.092530406652131  -0.092530406652131

df = -626.0032Hz
  Columns 1 through 8

   0.001463697889157   0.001465578137679   0.001465578137679   0.001465988976231   0.001465988976231   0.001465694475254   0.001465426560187   0.001218346599528
   0.000000622307415   0.000000622307415   0.000000622307415   0.000000622307415   0.000000622307415  -0.000002109827884  -0.000002109827884  -0.000511820432195
                   0                   0                   0                   0                   0                   0                   0                   0
                   0                   0                   0                   0                   0                   0                   0                   0
   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321558233241   0.010321557358685
  -0.092530461362482  -0.092530461361903  -0.092530461361903  -0.092530461361777  -0.092530461361777  -0.092530461361546  -0.092530461361267  -0.092530418880150

  Columns 9 through 10

   0.001170070951314   0.001170070951314
  -0.000511820432195  -0.000511820432195
                   0                   0
                   0                   0
   0.010321557358685   0.010321557358685
  -0.092530406652131  -0.092530406652131

for df and dct the values are identical, while for dp, they are not (but almost).

May be I do some error when computing the equivalent dct and df?

best regards
Simone

[simoneliuzzo]

Dear @lfarv
I tried also

orb6 = findorbit6(ring,1:10,'dp',dp,'dct',dct,'df',df);

and

orb6 = findorbit6(ring,1:10,'dp',dp,'dct',dct*3,'df',df*0.1);

in the first case the result is identical to the one when giving as input only dct or only df

in the second case the result is a different orbit, not seen in any of the other cases. How are the 3 different inputs handled when requested together? May be such requests should not be possible for the user?

[simoneliuzzo]

Dear @lfarv,

in the master branch the dp, dct and df inputs are ignored. I thought they would error for wrong input given as other matlab functions.

>> reshape(magic(5),1,25)
ans =
    17    23     4    10    11    24     5     6    12    18     1     7    13    19    25     8    14    20    21     2    15    16    22     3     9
>> reshape(magic(5),1,25,'badinput',98)
Error using reshape
Size arguments must be integer scalars.

This PR will solves this issue for the three inputs dp, dct and df, giving them a meaning.

However there are an infinity of possible wrong input keywords. Using the standard matlab input parser this is easily achieved.

function testinputs(r,varargin)

p = inputParser;

addRequired(p,'r',@iscell);
addOptional(p,'mux',0.0,@isnumeric);

parse(p,r,varargin{:});

Here the ouptut of a request for a non existing input

>> testinputs({1,2},'mux',23)
>> testinputs({1,2},'mux',23,'pippo',nan)
Error using testinputs (line 35)
'pippo' is not a recognized parameter.

[simoneliuzzo]

Dear @lfarv

not only for findorbit6 but also for atlinopt6, tunechrom and findm66 the results of input dp are different from those of input dct and df. The dp input appears to be the one that is different from the result obtained from a manual frequency change.

[simoneliuzzo]

Dear @lfarv
I attach here the two functions I have set up to test this branch.
test_matlab_frequency_control.zip

[lfarv]

@simoneliuzzo: tanks for the feedback ! I'll answer one by one:

I figured out the input for findoribt6. I think the help should be updated as for atlinopt6.
The warning issued when using dp, dct or df input always mentions dp.

Right. I'll update the help !

in the master branch the dp, dct and df inputs are ignored. I thought they would error for wrong input given as other matlab functions.

The choice between warning and errors was discussed at length: see #373. The question is now solved !

[lfarv]

@simoneliuzzo:

for df and dct the values are identical, while for dp, they are not (but almost).
May be I do some error when computing the equivalent dct and df?

Your computation is approximate because alpha is non linear. The correct computation is:

        function [dct,df]=dctdf(r4, dpp)
            [frf,l]=atGetRingProperties(r4,'rf_frequency','cell_length');
            [~,o0]=findorbit4(r4,dp=dpp);
            o1=ringpass(r4, o0);
            dct=o1(6);
            df=-frf*dct/(l+dct);
        end

Here you get the exact path lengthening from tracking. The results are then much closer (though they will never be exactly identical).

[lfarv]

Update for my previous post: the r4 input in the function is r4=atdisable_6d(ring): that's how you get the path lengthening for a given dp

[simoneliuzzo]

function [dct,df]=dctdf(r4, dpp)
            [frf,l]=atGetRingProperties(r4,'rf_frequency','cell_length');
            [~,o0]=findorbit4(r4,dp=dpp);
            o1=ringpass(r4, o0);
            dct=o1(6);
            df=-frf*dct/(l+dct);
        end

Dear @lfarv,

this function could be usefull for AT. As a common user, I always convert dp in df using the momentum compaction factor.

[lfarv]

@simoneliuzzo:

orb6 = findorbit6(ring,1:10,'dp',dp,'dct',dct,'df',df);

As you mentioned, such an input is stupid, so you get a stupid answer… Seriously, the results depends on the following priority order: df, dct and finally dp.

It would be better to forbid such duplicate contraints. But that's not limited to the functions modified in this PR. It can be a topic for another pull request

[simoneliuzzo]

@simoneliuzzo:

orb6 = findorbit6(ring,1:10,'dp',dp,'dct',dct,'df',df);

As you mentioned, such an input is stupid, so you get a stupid answer… Seriously, the results depends on the following priority order: df, dct and finally dp.

It would be better to forbid such duplicate contraints. But that's not limited to the functions modified in this PR. It can be a topic for another pull request

I agree. Also as I mentioned in a previous message, matlab offers already the solution to this problem, but it would need a lot of work.

[simoneliuzzo]

@simoneliuzzo:

for df and dct the values are identical, while for dp, they are not (but almost).
May be I do some error when computing the equivalent dct and df?

Your computation is approximate because alpha is non linear. The correct computation is:

        function [dct,df]=dctdf(r4, dpp)
            [frf,l]=atGetRingProperties(r4,'rf_frequency','cell_length');
            [~,o0]=findorbit4(r4,dp=dpp);
            o1=ringpass(r4, o0);
            dct=o1(6);
            df=-frf*dct/(l+dct);
        end

Here you get the exact path lengthening from tracking. The results are then much closer (though they will never be exactly identical).

Dear @lfarv

using this function for the conversion all comparisons are correct.

Only the help of findorbit6 is missing the description of the dp, df, dct arguments.

[lfarv]

this function could be usefull for AT.

The idea in this pull request is that this conversion should not be necessary any more: you can use the way you find the most convenient.

There is a better implementation in atsetcavity:

            hh=props_harmnumber(harmcell,cell_h);
            if isfinite(df)
                frev = frev + df/hh;
            elseif isfinite(dct)
                frev=frev * lcell/(lcell+dct);
            elseif isfinite(dp)
                % Find the path lengthening for dp
                [~,rnorad]=check_radiation(ring,false,'force');
                [~,orbitin]=findorbit4(rnorad,dp);
                orbitout=ringpass(rnorad,orbitin);
                dct=orbitout(6);
                frev=frev * lcell/(lcell+dct);
            end
            frequency = hh * frev;

And that's where the priority is defined if you give several keywords.

But if you find this function useful, it can be cleaned up and added.

[lfarv]

I agree. Also as I mentioned in a previous message, matlab offers already the solution to this problem, but it would need a lot of work.

Agreed. But the best way is not to use the Matlab parser: practically, in AT, a function "consumes" the keyword it needs and transmits the remaining ones to the sub functions it calls (example: atlinopt6 -> findm66 -> findorbit6). So the upper level function does not have the full list of authorised keywords. It would be the task of the lowest level functions to check in the end, after all keywords have been consumed, that the remaining list is empty.

As you said: a lot of work…

[lfarv]

@simoneliuzzo : the help of findorbit6 is corrected, and there is now an error if more than one keyword in dp, dct or df is specified.

Thank for your feedback: what you pointed out is out of the capabilities of automatic test sequences.

Because of the new approval rules, you'll have to approve again after my commit…