thinmpoleraddiffm

FINDTHINMPOLERADDIFFM

0001 function [B66, M, rout] = findthinmpoleraddiffm(rin, PolynomA, PolynomB, L, irho, E0, max_order)
0002 %FINDTHINMPOLERADDIFFM
0003 
0004 % Physical constants used in calculations
0005 persistent TWOPI CGAMMA M0C2 LAMBDABAR CER CU
0006 if isempty(TWOPI) %Initialize constansts on the first call
0007     TWOPI       = 2*pi;
0008     CGAMMA      = 8.846056192e-05;             % [m]/[GeV^3] Ref[1] (4.1)
0009     M0C2        = 5.10999060e5;              % Electron rest mass [eV]
0010     LAMBDABAR   = 3.86159323e-13;            % Compton wavelength/2pi [m]
0011     CER         = 2.81794092e-15;            % Classical electron radius [m]
0012     CU          = 1.323094366892892;         % 55/(24*sqrt(3))
0013 end
0014 
0015 
0016 % Calculate field from polynomial coefficients
0017 P1 = i*PolynomA(1:max_order+1)+PolynomB(1:max_order+1);
0018 Z1 = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order)]);
0019 S1 = sum(P1.*Z1);
0020 Bx = real(S1);
0021 By = imag(S1);
0022 
0023 B2P = B2perp([Bx; By+irho; 0], irho, rin);
0024 B3P = B2P^(3/2);
0025 p_norm = 1/(1+rin(5));
0026 p_norm2 = p_norm^2;
0027 
0028 CRAD = CGAMMA*E0^3/(TWOPI*1e27);
0029 BB = CU * CER * LAMBDABAR *  (E0/M0C2)^5 * L * B3P * (1+rin(5))^4*...
0030                 (1+rin(1)*irho + (rin(2)^2+rin(4)^2)*p_norm2/2);
0031 
0032 % Propagate particle
0033 rout = rin;
0034 
0035 % Loss of energy (dp/p) due to radiation
0036 rout(5) = rin(5) - CRAD*(1+rin(5))^2*B2P*...
0037     (1+rin(1)*irho + (rin(1)^2+rin(3)^2)*p_norm2/2)*L;
0038 
0039 % Change in transverse momentum due to radiation
0040 %   Angle does not change but dp/p changes due to radiation
0041 %   and therefore transverse canonical momentum changes
0042 %   px = x'*(1+dp/p)
0043 %   py = y'*(1+dp/p)
0044 rout([2 4]) = rin([2 4])*(1+rout(5))/(1+rin(5));
0045 
0046 % transverse kick due to magnetic field
0047 rout(2) = rout(2) - L*(Bx-(rin(5)-rin(1)*irho)*irho);
0048 rout(4) = rout(4) + L*By;
0049 
0050 % pathlength
0051 rout(6) = rout(6) + L*irho*rin(1); 
0052 
0053 
0054 % Calculate transfer matrix at rin
0055 P2 = i*PolynomA(2:max_order+1)+PolynomB(2:max_order+1);
0056 Z2 = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order-1)]);
0057 S2 = sum(P2.*(1:max_order).*Z2);
0058 
0059 M = eye(6);
0060 M(2,1)   = -L*real(S2);
0061 M(2,3)   =  L*imag(S2);
0062 M(4,1)   =  L*imag(S2);
0063 M(4,3)   =  L*real(S2);
0064 M(2,5)   =  L*irho;
0065 M(2,1)   =  M(2,1) - L*irho*irho;
0066 M(6,1)   =  L*irho;
0067 
0068 
0069 %    Calculate Ohmi's diffusion matrix of a thin multipole  element
0070 %    For elements with straight coordinate system irho = 0
0071 %    For curved elements the B polynomial (PolynomB in MATLAB)
0072 %    MUST NOT include the guide field  By0 = irho * E0 /(c*e)
0073 
0074 B66 = zeros(6);
0075 B66(2,2)    = BB*rin(2)^2*p_norm2;
0076 B66(2,4)    = BB*rin(2)*rin(4)*p_norm2;
0077 B66(4,2)    = B66(2,4);
0078 B66(4,4)    = BB*rin(4)^2*p_norm2;
0079 B66(5,2)    = BB*rin(2)*p_norm;
0080 B66(2,5)    = B66(5,2);
0081 B66(5,4)    = BB*rin(4)*p_norm;
0082 B66(4,5)    = B66(5,4);
0083 B66(5,5)    = BB;
0084 
0085 function b2 = B2perp(B, irho, rin)
0086 % Calculates sqr(|e x B|) , where 'e' is a unit vector in the direction of
0087 % velocity. Components of the  velocity vector:
0088 % ex = xpr;
0089 % ey = ypr;
0090 % ez = (1+x*irho);
0091 
0092 E = [rin(2)/(1+rin(5));rin(4)/(1+rin(5));1+rin(1)*irho];
0093 b2 = sum(cross(E/norm(E),B).^2);
0094