1 | function [B66, M, rout] = findthinmpoleraddiffm(rin, PolynomA, PolynomB, L, irho, E0, max_order) |
---|
2 | %FINDTHINMPOLERADDIFFM |
---|
3 | |
---|
4 | % Physical constants used in calculations |
---|
5 | persistent TWOPI CGAMMA M0C2 LAMBDABAR CER CU |
---|
6 | if isempty(TWOPI) %Initialize constansts on the first call |
---|
7 | TWOPI = 2*pi; |
---|
8 | CGAMMA = 8.846056192e-05; % [m]/[GeV^3] Ref[1] (4.1) |
---|
9 | M0C2 = 5.10999060e5; % Electron rest mass [eV] |
---|
10 | LAMBDABAR = 3.86159323e-13; % Compton wavelength/2pi [m] |
---|
11 | CER = 2.81794092e-15; % Classical electron radius [m] |
---|
12 | CU = 1.323094366892892; % 55/(24*sqrt(3)) |
---|
13 | end |
---|
14 | |
---|
15 | |
---|
16 | % Calculate field from polynomial coefficients |
---|
17 | P1 = i*PolynomA(1:max_order+1)+PolynomB(1:max_order+1); |
---|
18 | Z1 = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order)]); |
---|
19 | S1 = sum(P1.*Z1); |
---|
20 | Bx = real(S1); |
---|
21 | By = imag(S1); |
---|
22 | |
---|
23 | B2P = B2perp([Bx; By+irho; 0], irho, rin); |
---|
24 | B3P = B2P^(3/2); |
---|
25 | p_norm = 1/(1+rin(5)); |
---|
26 | p_norm2 = p_norm^2; |
---|
27 | |
---|
28 | CRAD = CGAMMA*E0^3/(TWOPI*1e27); |
---|
29 | BB = CU * CER * LAMBDABAR * (E0/M0C2)^5 * L * B3P * (1+rin(5))^4*... |
---|
30 | (1+rin(1)*irho + (rin(2)^2+rin(4)^2)*p_norm2/2); |
---|
31 | |
---|
32 | % Propagate particle |
---|
33 | rout = rin; |
---|
34 | |
---|
35 | % Loss of energy (dp/p) due to radiation |
---|
36 | rout(5) = rin(5) - CRAD*(1+rin(5))^2*B2P*... |
---|
37 | (1+rin(1)*irho + (rin(1)^2+rin(3)^2)*p_norm2/2)*L; |
---|
38 | |
---|
39 | % Change in transverse momentum due to radiation |
---|
40 | % Angle does not change but dp/p changes due to radiation |
---|
41 | % and therefore transverse canonical momentum changes |
---|
42 | % px = x'*(1+dp/p) |
---|
43 | % py = y'*(1+dp/p) |
---|
44 | rout([2 4]) = rin([2 4])*(1+rout(5))/(1+rin(5)); |
---|
45 | |
---|
46 | % transverse kick due to magnetic field |
---|
47 | rout(2) = rout(2) - L*(Bx-(rin(5)-rin(1)*irho)*irho); |
---|
48 | rout(4) = rout(4) + L*By; |
---|
49 | |
---|
50 | % pathlength |
---|
51 | rout(6) = rout(6) + L*irho*rin(1); |
---|
52 | |
---|
53 | |
---|
54 | % Calculate transfer matrix at rin |
---|
55 | P2 = i*PolynomA(2:max_order+1)+PolynomB(2:max_order+1); |
---|
56 | Z2 = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order-1)]); |
---|
57 | S2 = sum(P2.*(1:max_order).*Z2); |
---|
58 | |
---|
59 | M = eye(6); |
---|
60 | M(2,1) = -L*real(S2); |
---|
61 | M(2,3) = L*imag(S2); |
---|
62 | M(4,1) = L*imag(S2); |
---|
63 | M(4,3) = L*real(S2); |
---|
64 | M(2,5) = L*irho; |
---|
65 | M(2,1) = M(2,1) - L*irho*irho; |
---|
66 | M(6,1) = L*irho; |
---|
67 | |
---|
68 | |
---|
69 | % Calculate Ohmi's diffusion matrix of a thin multipole element |
---|
70 | % For elements with straight coordinate system irho = 0 |
---|
71 | % For curved elements the B polynomial (PolynomB in MATLAB) |
---|
72 | % MUST NOT include the guide field By0 = irho * E0 /(c*e) |
---|
73 | |
---|
74 | B66 = zeros(6); |
---|
75 | B66(2,2) = BB*rin(2)^2*p_norm2; |
---|
76 | B66(2,4) = BB*rin(2)*rin(4)*p_norm2; |
---|
77 | B66(4,2) = B66(2,4); |
---|
78 | B66(4,4) = BB*rin(4)^2*p_norm2; |
---|
79 | B66(5,2) = BB*rin(2)*p_norm; |
---|
80 | B66(2,5) = B66(5,2); |
---|
81 | B66(5,4) = BB*rin(4)*p_norm; |
---|
82 | B66(4,5) = B66(5,4); |
---|
83 | B66(5,5) = BB; |
---|
84 | |
---|
85 | function b2 = B2perp(B, irho, rin) |
---|
86 | % Calculates sqr(|e x B|) , where 'e' is a unit vector in the direction of |
---|
87 | % velocity. Components of the velocity vector: |
---|
88 | % ex = xpr; |
---|
89 | % ey = ypr; |
---|
90 | % ez = (1+x*irho); |
---|
91 | |
---|
92 | E = [rin(2)/(1+rin(5));rin(4)/(1+rin(5));1+rin(1)*irho]; |
---|
93 | b2 = sum(cross(E/norm(E),B).^2); |
---|
94 | |
---|