-
Notifications
You must be signed in to change notification settings - Fork 31
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Update momentum aperture matlab (#786)
* updates matlab momentum acceptance to use the closed orbit * updates the Momentum Aperture arguments * improves verbose
- Loading branch information
1 parent
cb61c7a
commit 6cee58a
Showing
2 changed files
with
232 additions
and
98 deletions.
There are no files selected for viewing
134 changes: 116 additions & 18 deletions
134
atmat/atphysics/TouschekPiwinski/MomAperture_allRing.m
100755 → 100644
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 | ||
| map_l = maps(:,2); | ||
| map_h = maps(:,1); | ||
| return |
196 changes: 116 additions & 80 deletions
196
atmat/atphysics/TouschekPiwinski/momentum_aperture_at.m
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 && delta<deltalimit) || ( delta<=0 && delta>deltalimit) | ||
|
|
||
| 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 && delta<deltalimit) || ( delta<=0 && delta>deltalimit ) | ||
| 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; |