Momentum aperture on a closed orbit in AT matlab

[oscarxblanco]

Dear all,
I have had a look on the momentum aperture functions in AT matlab and they seem to be unmodified over the last 11 years.
https://github.com/atcollab/at/blob/master/atmat/atphysics/TouschekPiwinski/momentum_aperture_at.m
https://github.com/atcollab/at/blob/master/atmat/atphysics/TouschekPiwinski/MomAperture_allRing.m

They might still be used for nominal lattices without errors and a cavity timelag that sets the closed orbit to zero, but, when errors are present I think the calculation should be done using small (but not negligible, e.g. 1$\sigma$) offsets around the closed orbit.

I would like to share a new version of these two .m scripts that I wrote in order to be able to pass the closed orbit, and arbitrary offset on every point.

In the doc info you will find also the boundary width, a new explanation of the algorithm, and many new parameters that previously were fixed (or adapted to the needs in local AT copies).

Let me know if this of any use, so I would open a pull request.

Best regards,
o
momentum_aperture_at.m ( download available on the next comment)
MomAperture_allRing.m.txt

[oscarxblanco]

Dear all,
I realized today that one of the files I uploaded yesterday has an error that I just fixed.
It is inside momentum_aperture_at.m where I missed a call to reforbit, and therefore misses to use the given reference orbit recursively.

Now it is added in line 96.

# file momentum_aperture_at.m
                 ...
                'reforbit',par.reforbit ... # line 96
                 ...

I am deleting the file uploaded yesterday, and uploading the new file here. Please use this file
momentum_aperture_at.m.txt

I have compared the result in matlab wrt the momentum acceptance output in pyat for one point in a ring, with small rectangular physical aperture of 3 mm (to avoid large energy offsets), and the numerical result is the same. They differ by a small amount within the resolution I tested ($O(10^{-5})$).

% matlab result
1000 turns, point 1 of 1, 100%	dpP: 0.028	dpN:-0.030
dpP =   0.027820000000000
dpN =  -0.030260000000000
# python result
boundary= [-0.030267,  0.027827]

In case anybody would like to try it out here is matlab and python test codes:

%% AT matlab
Nturns=1000
corbit = findorbit6(THERING,positions);
offsets = corbit;
offsets(1,:) = offsets(1,:) + 1e-5;
offsets(3,:) = offsets(3,:) + 1e-5;
%%
[dpN,dpP] = MomAperture_allRing(THERING, ...
                    positions, ...
                    Nturns, ...
                    'verbose',1, ...
                    'reforbit',offsets, ...
                    'xyinitoffsets',0*zeros(length(positions),2), ...
                    'ntimessplit',3);

# pyat grid 2
resolution = 1e-5
amplitude = 0.01
divider = 10

nturns = 1000
orbit0,orbit=at.find_orbit6(ring,refpts)
offset = orbit
offset[:,0] = offset[:,0] + 1e-5
offset[:,2] = offset[:,2] + 1e-5
boundary, tracked, survive = at.get_momentum_acceptance(ring,
                                resolution,
                                amplitude,
                                nturns=nturns,
                                refpts=refpts,
                                offset=offset,
                                grid_mode=at.GridMode(2), # 2 is recursive, 1 is cartesian
                                divider=divider,
                                use_mp=True,
                                verbose=True,
                                shift_zero=0)

[oscarxblanco]

Included in #786