source: MML/trunk/at/doc_html/at/atphysics/thinmpoleraddiffm.html @ 4

Last change on this file since 4 was 4, checked in by zhangj, 10 years ago

Initial import--MML version from SOLEIL@2013
File size: 7.0 KB
Line 
1<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
2                "http://www.w3.org/TR/REC-html40/loose.dtd">
3<html>
4<head>
5  <title>Description of thinmpoleraddiffm</title>
6  <meta name="keywords" content="thinmpoleraddiffm">
7  <meta name="description" content="FINDTHINMPOLERADDIFFM">
8  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
9  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
10  <meta name="robots" content="index, follow">
11  <link type="text/css" rel="stylesheet" href="../../m2html.css">
12</head>
13<body>
14<a name="_top"></a>
15<div><a href="../../index.html">Home</a> &gt;  <a href="../index.html">at</a> &gt; <a href="index.html">atphysics</a> &gt; thinmpoleraddiffm.m</div>
16
17<!--<table width="100%"><tr><td align="left"><a href="../../index.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
18<td align="right"><a href="index.html">Index for at/atphysics&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->
19
20<h1>thinmpoleraddiffm
21</h1>
22
23<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
24<div class="box"><strong>FINDTHINMPOLERADDIFFM</strong></div>
25
26<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
27<div class="box"><strong>function [B66, M, rout] = findthinmpoleraddiffm(rin, PolynomA, PolynomB, L, irho, E0, max_order) </strong></div>
28
29<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
30<div class="fragment"><pre class="comment">FINDTHINMPOLERADDIFFM</pre></div>
31
32<!-- crossreference -->
33<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
34This function calls:
35<ul style="list-style-image:url(../../matlabicon.gif)">
36</ul>
37This function is called by:
38<ul style="list-style-image:url(../../matlabicon.gif)">
39</ul>
40<!-- crossreference -->
41
42<h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
43<ul style="list-style-image:url(../../matlabicon.gif)">
44<li><a href="#_sub1" class="code">function b2 = B2perp(B, irho, rin)</a></li></ul>
45<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
46<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [B66, M, rout] = findthinmpoleraddiffm(rin, PolynomA, PolynomB, L, irho, E0, max_order)</a>
470002 <span class="comment">%FINDTHINMPOLERADDIFFM</span>
480003
490004 <span class="comment">% Physical constants used in calculations</span>
500005 <span class="keyword">persistent</span> TWOPI CGAMMA M0C2 LAMBDABAR CER CU
510006 <span class="keyword">if</span> isempty(TWOPI) <span class="comment">%Initialize constansts on the first call</span>
520007     TWOPI       = 2*pi;
530008     CGAMMA      = 8.846056192e-05;             <span class="comment">% [m]/[GeV^3] Ref[1] (4.1)</span>
540009     M0C2        = 5.10999060e5;              <span class="comment">% Electron rest mass [eV]</span>
550010     LAMBDABAR   = 3.86159323e-13;            <span class="comment">% Compton wavelength/2pi [m]</span>
560011     CER         = 2.81794092e-15;            <span class="comment">% Classical electron radius [m]</span>
570012     CU          = 1.323094366892892;         <span class="comment">% 55/(24*sqrt(3))</span>
580013 <span class="keyword">end</span>
590014
600015
610016 <span class="comment">% Calculate field from polynomial coefficients</span>
620017 P1 = i*PolynomA(1:max_order+1)+PolynomB(1:max_order+1);
630018 Z1 = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order)]);
640019 S1 = sum(P1.*Z1);
650020 Bx = real(S1);
660021 By = imag(S1);
670022
680023 B2P = <a href="#_sub1" class="code" title="subfunction b2 = B2perp(B, irho, rin)">B2perp</a>([Bx; By+irho; 0], irho, rin);
690024 B3P = B2P^(3/2);
700025 p_norm = 1/(1+rin(5));
710026 p_norm2 = p_norm^2;
720027
730028 CRAD = CGAMMA*E0^3/(TWOPI*1e27);
740029 BB = CU * CER * LAMBDABAR *  (E0/M0C2)^5 * L * B3P * (1+rin(5))^4*<span class="keyword">...</span>
750030                 (1+rin(1)*irho + (rin(2)^2+rin(4)^2)*p_norm2/2);
760031
770032 <span class="comment">% Propagate particle</span>
780033 rout = rin;
790034
800035 <span class="comment">% Loss of energy (dp/p) due to radiation</span>
810036 rout(5) = rin(5) - CRAD*(1+rin(5))^2*B2P*<span class="keyword">...</span>
820037     (1+rin(1)*irho + (rin(1)^2+rin(3)^2)*p_norm2/2)*L;
830038
840039 <span class="comment">% Change in transverse momentum due to radiation</span>
850040 <span class="comment">%   Angle does not change but dp/p changes due to radiation</span>
860041 <span class="comment">%   and therefore transverse canonical momentum changes</span>
870042 <span class="comment">%   px = x'*(1+dp/p)</span>
880043 <span class="comment">%   py = y'*(1+dp/p)</span>
890044 rout([2 4]) = rin([2 4])*(1+rout(5))/(1+rin(5));
900045
910046 <span class="comment">% transverse kick due to magnetic field</span>
920047 rout(2) = rout(2) - L*(Bx-(rin(5)-rin(1)*irho)*irho);
930048 rout(4) = rout(4) + L*By;
940049
950050 <span class="comment">% pathlength</span>
960051 rout(6) = rout(6) + L*irho*rin(1);
970052
980053
990054 <span class="comment">% Calculate transfer matrix at rin</span>
1000055 P2 = i*PolynomA(2:max_order+1)+PolynomB(2:max_order+1);
1010056 Z2 = cumprod([1, (rin(1)+i*rin(3))*ones(1,max_order-1)]);
1020057 S2 = sum(P2.*(1:max_order).*Z2);
1030058
1040059 M = eye(6);
1050060 M(2,1)   = -L*real(S2);
1060061 M(2,3)   =  L*imag(S2);
1070062 M(4,1)   =  L*imag(S2);
1080063 M(4,3)   =  L*real(S2);
1090064 M(2,5)   =  L*irho;
1100065 M(2,1)   =  M(2,1) - L*irho*irho;
1110066 M(6,1)   =  L*irho;
1120067
1130068
1140069 <span class="comment">%    Calculate Ohmi's diffusion matrix of a thin multipole  element</span>
1150070 <span class="comment">%    For elements with straight coordinate system irho = 0</span>
1160071 <span class="comment">%    For curved elements the B polynomial (PolynomB in MATLAB)</span>
1170072 <span class="comment">%    MUST NOT include the guide field  By0 = irho * E0 /(c*e)</span>
1180073
1190074 B66 = zeros(6);
1200075 B66(2,2)    = BB*rin(2)^2*p_norm2;
1210076 B66(2,4)    = BB*rin(2)*rin(4)*p_norm2;
1220077 B66(4,2)    = B66(2,4);
1230078 B66(4,4)    = BB*rin(4)^2*p_norm2;
1240079 B66(5,2)    = BB*rin(2)*p_norm;
1250080 B66(2,5)    = B66(5,2);
1260081 B66(5,4)    = BB*rin(4)*p_norm;
1270082 B66(4,5)    = B66(5,4);
1280083 B66(5,5)    = BB;
1290084
1300085 <a name="_sub1" href="#_subfunctions" class="code">function b2 = B2perp(B, irho, rin)</a>
1310086 <span class="comment">% Calculates sqr(|e x B|) , where 'e' is a unit vector in the direction of</span>
1320087 <span class="comment">% velocity. Components of the  velocity vector:</span>
1330088 <span class="comment">% ex = xpr;</span>
1340089 <span class="comment">% ey = ypr;</span>
1350090 <span class="comment">% ez = (1+x*irho);</span>
1360091
1370092 E = [rin(2)/(1+rin(5));rin(4)/(1+rin(5));1+rin(1)*irho];
1380093 b2 = sum(cross(E/norm(E),B).^2);
1390094</pre></div>
140<hr><address>Generated on Mon 21-May-2007 15:26:45 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
141</body>
142</html>
Note: See TracBrowser for help on using the repository browser.