diff --git a/atmat/atphysics/TouschekPiwinski/MomAperture_allRing.m b/atmat/atphysics/TouschekPiwinski/MomAperture_allRing.m old mode 100755 new mode 100644 index 071894982..59eab2fd3 --- a/atmat/atphysics/TouschekPiwinski/MomAperture_allRing.m +++ b/atmat/atphysics/TouschekPiwinski/MomAperture_allRing.m @@ -1,33 +1,131 @@ function [map_l,map_h]=MomAperture_allRing(THERING,points,varargin) -% all Ring momentum aperture -%points=1:10:length(THERING); +%MomAperture_allRing returns positive and negative momentum aperture +% boundaries where the particle is still alive. +% +% The boundary width (i.e. the uncertainty) is equal to +% energystep / (splitdivisor^ntimessplit), +% meaning that one more step of this size makes the +% particle unstable. +% +% [map_l,map_h] ... +% = MomAperture_allRing( +% THERING, ... +% POINTS ... +% ) +% [...] = MomAperture_allRing(..., NTURNS) +% Tracks over NTURNS to get the momentum aperture. Default 100 +% e.g [dppM,dppP]=MomAperture_allRing(THERING,positions,nturns) +% +% [...] = MomAperture_allRing(..., 'reforbit',ORBITIN) +% The initial particle coordinates are taken from ORBITIN. +% Default zeros(6,length(POINTS)) +% +% [...] = MomAperture_allRing(..., 'xyinitoffsets',[x y]) +% The transverse offsets to add to the reference orbit. +% Default 1e-5*ones(length(POINTS),2) +% +% [...] = MomAperture_allRing(..., 'deltalimits',[deltapos deltaneg]) +% The energy offset limits. Default [0.1 -0.1] +% +% [...] = MomAperture_allRing(..., 'initialguess',[posguess negguess]) +% The starting point of the recursive energy offsets. +% Default [0.0 0.0] +% +% [...] = MomAperture_allRing(..., 'energysteps',[posstep negstep]) +% The positive and negative initial energy steps. +% Default [0.01 -0.01] +% +% [...] = MomAperture_allRing(..., 'ntimessplit',NSPLIT) +% The number of recursive calls reducing the step size. Default 2 +% +% [...] = MomAperture_allRing(..., 'splitdivisor',SPLITTDIVISOR) +% The step divisor every time we split the step. Default 10. +% +% [...] = MomAperture_allRing(..., 'verbose',VERBOSE) +% Print info the current position. Default 0. +% If set to 1 it will print info at every reference point. +% If set to 2 it will print info at each energy step. +% +% ex: [map_l,map_h] = MomAperture_allRing(THERING,points,nturns); +% ex: [map_l,map_h] = ... +% MomAperture_allRing(THERING,points,nturns,'reforbit',findorbit6(THERING,points)); + +% 2024apr09 oblancog at ALBA : adds parameters, and doc info map_h=zeros(length(points),1); map_l=zeros(length(points),1); -if numel(varargin)>0 +% nturns is historically a positional argument, we work around to find it +if mod(numel(varargin),2) == 1 nturns=varargin{1}; + varargin(1) = []; else nturns=100; end +p = inputParser; +addOptional(p,'reforbit',zeros(6,length(points))); +addOptional(p,'xyinitoffsets',10e-6*ones(length(points),2)); +addOptional(p,'deltalimits', [0.1 -0.1]); +addOptional(p,'initialguess', [0.0 0.0]); +addOptional(p,'energysteps', [0.01 -0.01]); +addOptional(p,'ntimessplit',2); +addOptional(p,'splitdivisor',10); +addOptional(p,'verbose',0); +parse(p,varargin{:}); +par = p.Results; + +verbose = par.verbose; +deepverbose = false; +if verbose == 2 + deepverbose = true; +end + +if ~isequal(size(par.reforbit),[6,length(points)]) + fprintf('ERROR: The size of the reference orbit needs to be (6,number of points)\n'); + return; +end + +lengthpoints = length(points); +maps = zeros(lengthpoints,2); +deltalimits = par.deltalimits; +esteps = par.energysteps; for i=1:length(points) - % disp([i length(points) i/length(points)*100]) %cycle ring - THERING_cycl=[THERING(points(i):end); THERING(1:points(i)-1)]'; - try - map_h(i)=momentum_aperture_at(THERING_cycl,+0.1,[10^-6 10^-6],0,0,+0.01,2,10,nturns); - catch - map_h(i)=0; - end - + THERING_cycl=[THERING(points(i):end); THERING(1:points(i)-1)]'; + + for j = 1:length(deltalimits) + % try catch here is used to return a valid value even if the momap + % calculation has problems. At the end it is very expensive to + % throw an error and exit, better return zero if everything fails try - map_l(i)=momentum_aperture_at(THERING_cycl,-0.1,[10^-6 10^-6],0,0,-0.01,2,10,nturns); + maps(i,j)=momentum_aperture_at(... + THERING_cycl, ... + deltalimits(j), ... + par.xyinitoffsets(i,:), ... + par.initialguess(j), ... + par.initialguess(j), ... + esteps(j), ... + par.ntimessplit, ... + par.splitdivisor, ... + nturns, ... + 'reforbit',par.reforbit(:,i), ... + 'verbose', deepverbose ... + ); catch - map_l(i)=0; + maps(i,j)=0; end - - + end + if verbose + fprintf('%d turns, point %d of %d, %.0f%%\tdppP:%6.3f\tdppN:%6.3f\n', ... + nturns, ... + i, ... + lengthpoints, ... + i/lengthpoints*100, ... + maps(i,1), ... + maps(i,2) ... + ); + end end - - -return \ No newline at end of file +map_l = maps(:,2); +map_h = maps(:,1); +return diff --git a/atmat/atphysics/TouschekPiwinski/momentum_aperture_at.m b/atmat/atphysics/TouschekPiwinski/momentum_aperture_at.m index 9b5b640ca..2599bffab 100755 --- a/atmat/atphysics/TouschekPiwinski/momentum_aperture_at.m +++ b/atmat/atphysics/TouschekPiwinski/momentum_aperture_at.m @@ -1,80 +1,116 @@ -function [deltamax]=momentum_aperture_at(THERING,deltalimit,initcoord,delta,precdelta,deltastepsize,splits,split_step_divisor,nturns) -% function [deltamin, deltamax... -% ]=momentum_aperture_at(THERING,... -% deltalimit,... [min max] -% initcoord,... [x y] initial coordinate -% delta,... -% precdelta,... -% deltastepsize,... -% splits,... % number of splitting -% split_step_divisor) % divide the step size at every split -% -% following the ELEGANT routine: -% Start with ? = 0, i.e., zero momentum offset. -% 2. Track a particle to see if it gets lost. If so, proceed to step 4. -% 3. Increase ? by step size ?? and return to step 2. -% 4. If no splitting steps remain, proceed to the next step. Otherwise: -% (a) Change ? to deltas ? sb??., where ?s is the largest ? for which the particle survived, and sb is the steps_back parameter. -% (b) Divide the step size by split_step_divisor to get a new step size ??. -% (c) Set?=?+??. -% (d) Decrement the ?splits remaining? counter by 1. -% (e) Continue from step 2. -% 5. Stop. The momentum aperture is ?s -% -% ex: [deltamax]=momentum_aperture_at(THERING,0.1,[10^-6 10^-6],0,0,0.01,3,10,100) -% ex: [deltamin]=momentum_aperture_at(THERING,-0.1,[10^-6 10^-6],0,0,-0.01,3,10,100) - -%disp([delta splits]) - -if ( delta>=0 && deltadeltalimit) - - if splits>-1 - - % track for this delta - - [~, LOSS] =ringpass(THERING,[initcoord(1) 0 initcoord(2) 0 delta 0]',nturns); - - if LOSS~=1 % if NOT LOST go to next step - - [deltamax... - ]=momentum_aperture_at(THERING,... - deltalimit,... [min max] - initcoord,... [x y] - delta+deltastepsize,... % delta center - delta,... - deltastepsize,... - splits,... % number of splitting - split_step_divisor,... - nturns); - - else % if LOST reduce stepsize - - [deltamax... - ]=momentum_aperture_at(THERING,... - deltalimit,... [min max] - initcoord,... [x y] - precdelta+deltastepsize/split_step_divisor,... % go back to previous delta center and increase of smaller step - precdelta,... - deltastepsize/split_step_divisor,... - splits-1,... % number of splitting - split_step_divisor,... - nturns); - - end - else - - % no splitting steps remain - deltamax=delta-deltastepsize; - - end - -else - % limit reached - deltamax=delta; -end - - -return; - - - +function deltamax = momentum_aperture_at( ... + THERING, ... + deltalimit, ... + initcoord, ... + delta, ... + precdelta, ... + deltastepsize, ... + splits, ... + split_step_divisor, ... + nturns, ... + varargin ... + ) +%momentum_aperture_at recursively offsets the particle energy and checks +% for survival over n turns of tracking. +% Returns the stable energy boundary. +% +% deltamax ... +% = momentum_aperture_at( ... +% THERING,... +% deltalimit,... [min max] +% initcoord,... [x y] initial coordinate +% delta,... current energy offset +% precdelta,... previous energy offset +% deltastepsize,... +% splits,... number of times splitting the deltastepsize +% split_step_divisor, divides the step size at every split +% nturns +% ) +% +% ... = momentum_aperture_at(..., 'reforbit',ORBITIN) +% Use ORBITIN to define the reference. Useful when the closed orbit +% is not zero. +% +% Adapted routine based on ELEGANT +% 1. Start with delta = 0, i.e., zero momentum offset. +% 2. If the limit has been reached stop, otherwise +% If the number of step divisions is done, stop. Otherwise ... +% Track the particle +% If it survives, increase the energy by one step, and start 2) again. +% Otherwise, go back one step in energy, and divide the step. +% Count the number of times the step has been divided. +% Start 2) with the new step. +% +% Debugging info prints are commented to avoid speed reduction, +% +% The ELEGANT routine: +% https://ops.aps.anl.gov/manuals/elegant_latest/elegantsu53.html +% +% ex: [deltamax]=momentum_aperture_at(THERING,0.1,[10^-6 10^-6],0,0,0.01,3,10,100) +% ex: [deltamin]=momentum_aperture_at(THERING,-0.1,[10^-6 10^-6],0,0,-0.01,3,10,100) + +% 2024apr09 oblanco at ALBA : adds reforbit option. + +p = inputParser; +addOptional(p,'reforbit',zeros(6,1)); +addOptional(p,'verbose',false); +parse(p,varargin{:}); +par = p.Results; + +verbose = par.verbose; + +offset6d = par.reforbit + [initcoord(1); 0; initcoord(2); 0; delta; 0]; +if ( delta>=0 && deltadeltalimit ) + if splits>-1 + % track for this delta + [~, LOSS] =ringpass(THERING,offset6d,nturns); + if LOSS~=1 % if NOT LOST go to next step + thedelta = delta+deltastepsize; + thepreviousdelta = delta; + thedeltastepsize = deltastepsize; + thesplits = splits; + else % if LOST reduce stepsize + % go back to previous delta center and increase of smaller step + thedelta = precdelta+deltastepsize/split_step_divisor; + thepreviousdelta = precdelta; + thedeltastepsize = deltastepsize/split_step_divisor; + thesplits = splits - 1; + end + if verbose + if LOSS~=1 + deadalivemsg = 'alive'; + else + deadalivemsg = 'dead'; + end + fprintf( ... + 'split%d\tdelta%.5f\tprecdelta%.5f\tdeltastepsize%.5f\t%s\n', ... + splits, ... + delta, ... + precdelta, ... + deltastepsize, ... + deadalivemsg ... + ); + end + deltamax ... + = momentum_aperture_at( ... + THERING, ... + deltalimit, ... [min max] + initcoord, ... [x y] + thedelta, ... % delta center + thepreviousdelta, ... + thedeltastepsize, ... + thesplits, ... % number of splitting + split_step_divisor, ... + nturns, ... + 'reforbit',par.reforbit, ... + 'verbose', verbose ... + ); + else + % no splitting steps remain + deltamax=delta-deltastepsize; + end +else + % limit reached + deltamax=delta; +end +return;