[4] | 1 | function orbit = findorbit4(RING,dP,varargin); |
---|
| 2 | %FINDORBIT4 finds closed orbit in the 4-d transverse phase |
---|
| 3 | % space by numerically solving for a fixed point of the one turn |
---|
| 4 | % map M calculated with LINEPASS |
---|
| 5 | % |
---|
| 6 | % (X, PX, Y, PY, dP2, CT2 ) = M (X, PX, Y, PY, dP1, CT1) |
---|
| 7 | % |
---|
| 8 | % under the CONSTANT MOMENTUM constraint, dP2 = dP1 = dP and |
---|
| 9 | % there is NO constraint on the 6-th coordinate CT |
---|
| 10 | % |
---|
| 11 | % IMPORTANT!!! FINDORBIT4 imposes a constraint on dP and relaxes |
---|
| 12 | % the constraint on the revolution frequency. A physical storage |
---|
| 13 | % ring does exactly the opposite: the momentum deviation of a |
---|
| 14 | % particle on the closed orbit settles at the value |
---|
| 15 | % such that the revolution is synchronous with the RF cavity |
---|
| 16 | % |
---|
| 17 | % HarmNumber*Frev = Frf |
---|
| 18 | % |
---|
| 19 | % To impose this artificial constraint in FINDORBIT4 |
---|
| 20 | % PassMethod used for any element SHOULD NOT |
---|
| 21 | % 1. change the longitudinal momentum dP (cavities , magnets with radiation) |
---|
| 22 | % 2. have any time dependence (localized impedance, fast kickers etc) |
---|
| 23 | % |
---|
| 24 | % FINDORBIT4(RING,dP) is 4x1 vector - fixed point at the |
---|
| 25 | % entrance of the 1-st element of the RING (x,px,y,py) |
---|
| 26 | % |
---|
| 27 | % FINDORBIT4(RING,dP,REFPTS) is 4-by-Length(REFPTS) |
---|
| 28 | % array of column vectors - fixed points (x,px,y,py) |
---|
| 29 | % at the entrance of each element indexed REFPTS array. |
---|
| 30 | % REFPTS is an array of increasing indexes that select elements |
---|
| 31 | % from the range 1 to length(RING)+1. |
---|
| 32 | % See further explanation of REFPTS in the 'help' for FINDSPOS |
---|
| 33 | % |
---|
| 34 | % FINDORBIT4(RING,dP,REFPTS,GUESS) - same as above but the search |
---|
| 35 | % for the fixed point starts at the initial condition GUESS |
---|
| 36 | % Otherwise the search starts from [0 0 0 0 0 0]'. |
---|
| 37 | % GUESS must be a 6-by-1 vector; |
---|
| 38 | % |
---|
| 39 | % [ORBIT, FIXEDPOINT] = FINDORBIT4( ... ) |
---|
| 40 | % The optional second return parameter is |
---|
| 41 | % a 6-by-1 vector of initial conditions |
---|
| 42 | % on the closed orbit at the entrance to the RING. |
---|
| 43 | % |
---|
| 44 | % See also FINDORBIT, FINDSYNCORBIT, FINDORBIT6. |
---|
| 45 | |
---|
| 46 | if ~iscell(RING) |
---|
| 47 | error('First argument must be a cell array'); |
---|
| 48 | end |
---|
| 49 | |
---|
| 50 | %d = sqrt(eps); % step size for numerical differentiation |
---|
| 51 | d = 1e-8; |
---|
| 52 | D5 = d*eye(4); |
---|
| 53 | |
---|
| 54 | max_iterations = 20; |
---|
| 55 | J = zeros(4); |
---|
| 56 | |
---|
| 57 | % Check if guess argument was supplied |
---|
| 58 | if nargin == 4 |
---|
| 59 | if isnumeric(varargin{2}) & all(eq(size(varargin{2}),[6,1])) |
---|
| 60 | Ri = varargin{2}; |
---|
| 61 | else |
---|
| 62 | error('The last argument GUESS must be a 6-by-1 vector'); |
---|
| 63 | end |
---|
| 64 | else |
---|
| 65 | Ri = zeros(6,1); |
---|
| 66 | end |
---|
| 67 | % Set the momentum component of Ri to the specified dP |
---|
| 68 | Ri(5) = dP; |
---|
| 69 | |
---|
| 70 | %First iteration |
---|
| 71 | RMATi= Ri*ones(1,5); |
---|
| 72 | for k = 1:4 |
---|
| 73 | RMATi(k,k) = RMATi(k,k)+d; |
---|
| 74 | end |
---|
| 75 | RMATf = linepass(RING,RMATi); |
---|
| 76 | Rf = RMATf(:,5); |
---|
| 77 | % compute the transverse part of the Jacobian |
---|
| 78 | J4 = [RMATf(1:4,1:4)-RMATf(1:4,5)*ones(1,4)]/d; |
---|
| 79 | % Replace matrix inversion with \ |
---|
| 80 | % B4 = inv(eye(4) - J4); |
---|
| 81 | % Ri_next = Ri + [B4*(Rf(1:4)-Ri(1:4)); 0; 0]; |
---|
| 82 | |
---|
| 83 | Ri_next = Ri + [(eye(4) - J4)\(Rf(1:4)-Ri(1:4)); 0; 0]; |
---|
| 84 | |
---|
| 85 | change = norm(Ri_next - Ri); |
---|
| 86 | Ri = Ri_next; |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | itercount = 1; |
---|
| 90 | |
---|
| 91 | while (change>eps) & (itercount < max_iterations) |
---|
| 92 | RMATi=[Ri Ri Ri Ri Ri]; |
---|
| 93 | for k = 1:4 |
---|
| 94 | RMATi(k,k) = RMATi(k,k)+d; |
---|
| 95 | end |
---|
| 96 | RMATf = linepass(RING,RMATi,'reuse'); |
---|
| 97 | Rf = RMATf(:,5); |
---|
| 98 | |
---|
| 99 | % compute the transverse part of the Jacobian |
---|
| 100 | J4 = [RMATf(1:4,1:4)-RMATf(1:4,5)*ones(1,4)]/d; |
---|
| 101 | % Replace matrix inversion with \ |
---|
| 102 | % B4 = inv(eye(4) - J4); |
---|
| 103 | % Ri_next = Ri + [B4*(Rf(1:4)-Ri(1:4)); 0; 0]; |
---|
| 104 | |
---|
| 105 | Ri_next = Ri + [(eye(4) - J4)\(Rf(1:4)-Ri(1:4)); 0; 0]; |
---|
| 106 | change = norm(Ri_next - Ri); |
---|
| 107 | Ri = Ri_next; |
---|
| 108 | itercount = itercount+1; |
---|
| 109 | end; |
---|
| 110 | |
---|
| 111 | if(nargin<3) | (varargin{1}==(length(RING)+1)) |
---|
| 112 | % return only the fixed point at the entrance of RING{1} |
---|
| 113 | orbit=Ri(1:4,1); |
---|
| 114 | else % 3-rd input argument - vector of reference points along the RING |
---|
| 115 | % is supplied - return orbit |
---|
| 116 | orb6 = linepass(RING,Ri,varargin{1},'reuse'); |
---|
| 117 | orbit = orb6(1:4,:); |
---|
| 118 | end |
---|
| 119 | |
---|
| 120 | if nargout==2 |
---|
| 121 | varargout{1}=Ri; |
---|
| 122 | end |
---|