-
Notifications
You must be signed in to change notification settings - Fork 0
/
BeamOpticsTutorial.html
273 lines (261 loc) · 82.1 KB
/
BeamOpticsTutorial.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Beam optics tutorial (Section 3.3.1)</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left; }
.S1 { margin: 15px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left; }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S2 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
.S3 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
.S5 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left; }
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
.S7 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px; }
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement,.eoOutputWrapper .matrixElement { min-height: 18px; box-sizing: border-box;}
.embeddedOutputsMatrixElement .matrixElement,.eoOutputWrapper .matrixElement,.rtcDataTipElement .matrixElement { position: relative;}
.matrixElement .variableValue,.rtcDataTipElement .matrixElement .variableValue { white-space: pre; display: inline-block; vertical-align: top; overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper { display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body { padding-top: 0 !important; max-height: 100px;}
.inlineElement .matrixElement { max-height: 300px;}
.embeddedOutputsMatrixElement.rightPaneElement {}
.rightPaneElement .matrixElement,.rtcDataTipElement .matrixElement { overflow: hidden; padding-left: 9px;}
.rightPaneElement .matrixElement { margin-bottom: -1px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer,.eoOutputWrapper .matrixElement .valueContainer,.rtcDataTipElement .matrixElement .valueContainer { white-space: nowrap; margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis.hide,.eoOutputWrapper .matrixElement .verticalEllipsis.hide,.rtcDataTipElement .matrixElement .valueContainer .horizontalEllipsis.hide,.rtcDataTipElement .matrixElement .verticalEllipsis.hide { display: none;}
.embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .verticalEllipsis, .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .horizontalEllipsis { display:none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis { margin-bottom: -3px;}
.eoOutputWrapper .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer { cursor: default !important;}
.embeddedOutputsVariableElement { white-space: pre-wrap; word-wrap: break-word; min-height: 18px; max-height: 250px; overflow: auto;}
.variableElement {}
.embeddedOutputsVariableElement.inlineElement {}
.inlineElement .variableElement {}
.embeddedOutputsVariableElement.rightPaneElement { min-height: 16px;}
.rightPaneElement .variableElement { padding-top: 2px; padding-left: 9px;}
.variableNameElement { margin-bottom: 3px; display: inline-block;}
/* * Ellipses as base64 for HTML export. */.matrixElement .horizontalEllipsis,.rtcDataTipElement .matrixElement .horizontalEllipsis { display: inline-block; margin-top: 3px; /* base64 encoded version of images-liveeditor/HEllipsis.png */ width: 30px; height: 12px; background-repeat: no-repeat; background-image: url("data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAB0AAAAJCAYAAADO1CeCAAAAJUlEQVR42mP4//8/A70xw0i29BUDFPxnAEtTW37wWDqakIa4pQDvOOG89lHX2gAAAABJRU5ErkJggg==");}
.matrixElement .verticalEllipsis,.textElement .verticalEllipsis,.rtcDataTipElement .matrixElement .verticalEllipsis,.rtcDataTipElement .textElement .verticalEllipsis { margin-left: 35px; /* base64 encoded version of images-liveeditor/VEllipsis.png */ width: 12px; height: 30px; background-repeat: no-repeat; background-image: url("data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAoAAAAZCAYAAAAIcL+IAAAALklEQVR42mP4//8/AzGYgWyFMECMwv8QddRS+P//KyimlmcGUOFoOI6GI/UVAgDnd8Dd4+NCwgAAAABJRU5ErkJggg==");}
.S8 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px; }
.S9 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap; }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }</style></head><body><div class = rtcContent><div class = 'S0'><span>Companion software for "Volker Ziemann, </span><span style=' font-style: italic;'>Hands-on Accelerator physics using MATLAB, CRCPress, 2019</span><span>" (https://www.crcpress.com/9781138589940)</span></div><h1 class = 'S1'><span>Beam optics tutorial (Section 3.3.1)</span></h1><div class = 'S0'><span>Volker Ziemann, 211106, CC-BY-SA-4.0</span></div><div class = 'S0'><span>In this tutorial we'll develop a first beam optics code that allows us to track single particles, or "rays" through a sequence of drift spaces and thin lenses. Let's therefore first clear the MATLAB workspace and define anonymous functions (those with an '@') that return the transfer matrices. Here D(L) will return the 2x2 matrix for a drift space and Q(F) that of a thin lens with focal length F</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span >clear </span><span style="color: rgb(170, 4, 249);">all</span><span >; close </span><span style="color: rgb(170, 4, 249);">all</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >D=@(L)[1,L;0,1];</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre"><span >Q=@(F)[1,0;-1/F,1];</span></span></div></div></div><div class = 'S5'><span>Now it is easy to obtain the transfer matrix of a beam line consisting of a sequnece of such elements. For example, a beamline, where the first element is a drift space with a length of 2m, followed by a lens with a focal length of 1 m and a second drift space having a length of 3m.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div class = 'S6'><span style="white-space: pre"><span >A=D(3)*Q(1)*D(2)</span></span></div><div class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableMatrixElement" uid="AF17D69A" data-scroll-top="null" data-scroll-left="null" data-width="1334" data-testid="output_0" style="width: 1364px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="matrixElement veSpecifier saveLoad" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="veVariableName variableNameElement double" style="width: 1334px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="headerElementClickToInteract" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">A = </span><span class="veVariableValueSummary veMetaSummary" style="white-space: normal; font-style: normal; color: rgb(179, 179, 179); font-size: 12px;">2×2</span></div></div><div class="valueContainer" data-layout="{"columnWidth":44,"totalColumns":"2","totalRows":"2","charsPerColumn":6}" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="variableValue" style="width: 90px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> -2 -1
-1 -1
</div><div class="horizontalEllipsis hide" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div><div class="verticalEllipsis hide" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div></div></div></div></div></div></div><div class = 'S5'><span>I suggest to verify the numerical calculation by multiplying the individual matrices by hand on a peice of paper.</span></div><div class = 'S0'><span>Directly coding the matrix multiplications is inconvenient and we therefore introduce a very simple way to describe sequences of elements as an array, here called </span><span style=' font-family: monospace;'>fodo.</span><span> This array has four columns with the following interpretations</span></div><ul class = 'S8'><li class = 'S9'><span>First column: code, drift=1, quad=2</span></li><li class = 'S9'><span>Second column: repeat code for the element</span></li><li class = 'S9'><span>Third column: length of one segment</span></li><li class = 'S9'><span>Fourth column: strength of one segment, focal length for a thin quadrupole.</span></li></ul><div class = 'S0'><span>The second column with the repeat code is not strictly necessary but will make the graphs easier to interpret, because we can plot many intermediate points, rather than connecting values at the end of elements by a straight line.</span></div><div class = 'S0'><span>This particular sequence starts with an element having a 1 in the first column and is therefore a drift space. It consists of 5 segments, each 0.2 m long. The last column is not used here and contains a zero, The next element is a lens, beacuse the first code is 2 and it is defocusing, because the fourth column contains a negative number for the focal length. The third element is again a driftspace, this time with 10 segments of 0.2 m, followed by a focusing lens and another drift space. Such a sequence of focusing lens, drift space, defocusing lens and another drift space is usually referred to as a FODO (get it?) cell. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span >F=2.1;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >fodo=[1, 5, 0.2, 0;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > 2, 1, 0, -F;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > 1, 10, 0.2, 0;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > 2, 1, 0, F;</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S10'><span style="white-space: pre"><span > 1, 5, 0.2, 0] </span></span></div><div class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableMatrixElement" uid="C8B1D865" data-scroll-top="null" data-scroll-left="null" data-width="1334" data-testid="output_1" style="width: 1364px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="matrixElement veSpecifier saveLoad" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="veVariableName variableNameElement double" style="width: 1334px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="headerElementClickToInteract" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">fodo = </span><span class="veVariableValueSummary veMetaSummary" style="white-space: normal; font-style: normal; color: rgb(179, 179, 179); font-size: 12px;">5×4</span></div></div><div class="valueContainer" data-layout="{"columnWidth":72,"totalColumns":"4","totalRows":"5","charsPerColumn":10}" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="variableValue" style="width: 290px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> 1.0000 5.0000 0.2000 0
2.0000 1.0000 0 -2.1000
1.0000 10.0000 0.2000 0
2.0000 1.0000 0 2.1000
1.0000 5.0000 0.2000 0
</div><div class="horizontalEllipsis hide" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div><div class="verticalEllipsis hide" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div></div></div></div></div></div></div><div class = 'S5'><span>Now we build our beamline as a sequence of one or many FODO cells. The latter we achieve by simply concatenating</span><span style=' font-family: monospace;'> fodo</span><span> vertically (note the semicolon!). We will always use</span><span style=' font-family: monospace;'> beamline</span><span> as a reserved name of the beamline that we work on.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%beamline=fodo; % one cell</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%beamline=[fodo;fodo] % two cells</span></span></div></div><div class="inlineWrapper outputs"><div class = 'S10'><span style="white-space: pre"><span >beamline=repmat(fodo,10,1) </span><span style="color: rgb(2, 128, 9);">% ten cells</span></span></div><div class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableMatrixElement" uid="5CC5EFC4" data-scroll-top="null" data-scroll-left="null" data-width="1334" data-testid="output_2" style="width: 1364px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="matrixElement veSpecifier saveLoad" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="veVariableName variableNameElement double" style="width: 1334px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="headerElementClickToInteract" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">beamline = </span><span class="veVariableValueSummary veMetaSummary" style="white-space: normal; font-style: normal; color: rgb(179, 179, 179); font-size: 12px;">50×4</span></div></div><div class="valueContainer" data-layout="{"columnWidth":72,"totalColumns":"4","totalRows":"50","charsPerColumn":10}" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="variableValue" style="width: 290px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> 1.0000 5.0000 0.2000 0
2.0000 1.0000 0 -2.1000
1.0000 10.0000 0.2000 0
2.0000 1.0000 0 2.1000
1.0000 5.0000 0.2000 0
1.0000 5.0000 0.2000 0
2.0000 1.0000 0 -2.1000
1.0000 10.0000 0.2000 0
2.0000 1.0000 0 2.1000
1.0000 5.0000 0.2000 0
</div><div class="horizontalEllipsis hide" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div><div class="verticalEllipsis" style="white-space: nowrap; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"></div></div></div></div></div></div></div><div class = 'S5'><span>Before we can assemble the transfer matrices, we have to do a bit of accounting, like counting the number of elements </span><span style=' font-family: monospace;'>nlines</span><span> in the beamline description and the number </span><span style=' font-family: monospace;'>nmat</span><span> of matrices that we will need to store the matrices. This is just given by the sum of the segments in the second column plus one. Next we allocate space for the nmat </span><span texencoding="2\times 2" style="vertical-align:-5px"><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAEUAAAAkCAYAAADfJffWAAAEGklEQVRoQ+3YWei1UxQG8N9nyHRhHm7MhGSeFREyRYTMQyhzRIbIjTmRDKUolHkmIvOUKRdmmeKGC0oRkow9tY/e73zvOe8+///xfd/Fu29OnXftvdd+9rPWetaeox/zIDCnx2ReBHpQWljRg9KDUpcseqb0TJk9U8Kik7EvNsUy+BAv4Wr8VrfFArNaD+djc2yIb/ABrscb47waFT6r4kHsNGLyV9gfHy+wI4/f+MRy+KVbzP7BNThv1BKjQHkCO+A6vIclsRlOx/JlsXexHf5YyIDZtjDhGdyL77AODsSeDV8PwUNtvreBEuMbsE2hXHPeyngT65Y/D8d9CxkoH+FRXNzi19m4tvz/WQmrKkWbBW/FUyMOmxzzZPkWGp5bCcpquBRn4dfKOTHbBCcgB/q7Y17YfAu2R8JkeIQEb5cLz/dl8XObUfO/xfE0dh+zeRJuFsoGd+OoygNm3b3wSkneNcAkwb+AlXBSOfC47S7A++UMo+xSJAYXuT6+7AIl3xepuJEfsAKuwEWVoKyBl7E2XsU+HYxpAnIXjq3wKxfVxpCmi2eU9BC7pfB7DShdZ1yilOM4sEu5+a45g+9rFmDW6gAmYRCGrIh7cAz+qt2kw+5yXFj82LXNdiaKNhXnLXyKjWbgaIBJCOX3New9xJjoiucLIEniCc9pARJ3B2F8GO6fFig34bSSd3KbMxlhSkJpAExC6RdsUQBJaD6AI6YMSJJ9RNzrheWtoTYpU1bHJ7gTp84Ejcac5JYAk1wTxiQ3PVZyVfRDyv2fs9xjeHrU7HHYGl+MWntSUCLqIuRys9MQbU1gBj4+gkP/B0Ai6pLg98Nz48CeBJTI4iOxM36a4g0eUMRWlvyxVKf8TnOsUvJgdNLtXQvXgnIwLilx+H3XohN8b1aZMC86aRIdU7NVym6a2Ej+hE/nqAFlN9yIPfBt54r1BtEhL5Yqk5BJmXy25JgaHVOzU0COQk9ivbJmQmy6QNkRt5Wy+fWIRdMgpnJMkmOawuxhpDwmqTZzzGyBWbT0ZZ93CMy8CKRp/G+MA2WrQrkkpjRPbSPUTB+UZ4Qa2Z41moDkeSJlt1llhsv1sI6pufCc646So84cMyHN7wa4rAaUjUscXlV6ieF1s+lyOAW5ifzWjDR3CZn0MtEhSdxtZTfAJA/kN+V6oGNq9ojNzeVZI49MbU3kYuVy0klHgM6VFtqYkiYp1I3QqRlbIm8rXaMJSJRkABmnVJstwSTApHM/p8uZ8v1xpPrNNYZBSSudJ8eItJrxDlL/a0bCLM8OqQJHVyrVABPGJNccX1FOm+8lNT7Fn3meSLoSbc3CtTYBPPGdhmySXiaK96DyCli716zs5icos3J0fk7uQWlBuwelB6UuCHum9EzpmVKHQM+UOpz+BU0byyVjhDkSAAAAAElFTkSuQmCC" width="34.5" height="18" /></span><span> matrices in the array Racc and initialize the first matrix as the </span><span texencoding="2\times2" style="vertical-align:-5px"><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAEUAAAAkCAYAAADfJffWAAAEGklEQVRoQ+3YWei1UxQG8N9nyHRhHm7MhGSeFREyRYTMQyhzRIbIjTmRDKUolHkmIvOUKRdmmeKGC0oRkow9tY/e73zvOe8+///xfd/Fu29OnXftvdd+9rPWetaeox/zIDCnx2ReBHpQWljRg9KDUpcseqb0TJk9U8Kik7EvNsUy+BAv4Wr8VrfFArNaD+djc2yIb/ABrscb47waFT6r4kHsNGLyV9gfHy+wI4/f+MRy+KVbzP7BNThv1BKjQHkCO+A6vIclsRlOx/JlsXexHf5YyIDZtjDhGdyL77AODsSeDV8PwUNtvreBEuMbsE2hXHPeyngT65Y/D8d9CxkoH+FRXNzi19m4tvz/WQmrKkWbBW/FUyMOmxzzZPkWGp5bCcpquBRn4dfKOTHbBCcgB/q7Y17YfAu2R8JkeIQEb5cLz/dl8XObUfO/xfE0dh+zeRJuFsoGd+OoygNm3b3wSkneNcAkwb+AlXBSOfC47S7A++UMo+xSJAYXuT6+7AIl3xepuJEfsAKuwEWVoKyBl7E2XsU+HYxpAnIXjq3wKxfVxpCmi2eU9BC7pfB7DShdZ1yilOM4sEu5+a45g+9rFmDW6gAmYRCGrIh7cAz+qt2kw+5yXFj82LXNdiaKNhXnLXyKjWbgaIBJCOX3New9xJjoiucLIEniCc9pARJ3B2F8GO6fFig34bSSd3KbMxlhSkJpAExC6RdsUQBJaD6AI6YMSJJ9RNzrheWtoTYpU1bHJ7gTp84Ejcac5JYAk1wTxiQ3PVZyVfRDyv2fs9xjeHrU7HHYGl+MWntSUCLqIuRys9MQbU1gBj4+gkP/B0Ai6pLg98Nz48CeBJTI4iOxM36a4g0eUMRWlvyxVKf8TnOsUvJgdNLtXQvXgnIwLilx+H3XohN8b1aZMC86aRIdU7NVym6a2Ej+hE/nqAFlN9yIPfBt54r1BtEhL5Yqk5BJmXy25JgaHVOzU0COQk9ivbJmQmy6QNkRt5Wy+fWIRdMgpnJMkmOawuxhpDwmqTZzzGyBWbT0ZZ93CMy8CKRp/G+MA2WrQrkkpjRPbSPUTB+UZ4Qa2Z41moDkeSJlt1llhsv1sI6pufCc646So84cMyHN7wa4rAaUjUscXlV6ieF1s+lyOAW5ifzWjDR3CZn0MtEhSdxtZTfAJA/kN+V6oGNq9ojNzeVZI49MbU3kYuVy0klHgM6VFtqYkiYp1I3QqRlbIm8rXaMJSJRkABmnVJstwSTApHM/p8uZ8v1xpPrNNYZBSSudJ8eItJrxDlL/a0bCLM8OqQJHVyrVABPGJNccX1FOm+8lNT7Fn3meSLoSbc3CtTYBPPGdhmySXiaK96DyCli716zs5icos3J0fk7uQWlBuwelB6UuCHum9EzpmVKHQM+UOpz+BU0byyVjhDkSAAAAAElFTkSuQmCC" width="34.5" height="18" /></span><span>unit matrix eye(2). This additional matrix at the beginning accounts for the "plus one" in the definition of</span><span style=' font-family: monospace;'> nmat</span><span>. Finally we allocate the array </span><span style=' font-family: monospace;'>spos</span><span> to store the position of each segment along the beam line.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div class = 'S6'><span style="white-space: pre"><span >nlines=size(beamline,1)</span></span></div><div class = 'S7'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>nlines = 50</div></div></div><div class="inlineWrapper"><div class = 'S11'><span style="white-space: pre"><span >nmat=sum(beamline(:,2))+1;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >Racc=zeros(2,2,nmat);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >Racc(:,:,1)=eye(2);</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre"><span >spos=zeros(nmat,1);</span></span></div></div></div><div class = 'S5'><span>Now we are ready to loop through the beamline description by looping over the lines, labeled by</span><span style=' font-family: monospace;'> line</span><span>, and then over the segments, labeled by </span><span style=' font-family: monospace;'>seg</span><span>. As we do that we increment the counter</span><span style=' font-family: monospace;'> ic</span><span> that keeps track of which segment we're currently working on. Then we initialize the transfer matrix of the current segment </span><span style=' font-family: monospace;'>Rcurr</span><span> by the unit marix before looking at the first column to decide whether the element is a drift space or a lens and fill </span><span style=' font-family: monospace;'>Rcurr</span><span> with the appropriate transfer matrix. In</span><span style=' font-family: monospace;'> case</span><span> the first column (</span><span style=' font-family: monospace;'>beamline(line,1)</span><span>) is a 1, it is a drift and we use the value from the third column as the length of the drift space. In </span><span style=' font-family: monospace;'>case</span><span> the first column contains a 2, it is a lens and we use the focal length from the fourth column as the argument to the function </span><span style=' font-family: monospace;'>Q(F).</span><span> In case we encounter any other number in the first column, we display a warning. At this point </span><span style=' font-family: monospace;'>Rcurr</span><span> contains the transfer matrix for the current segment and we find </span><span style=' font-weight: bold;'>the transfer matrix from the start of the beamline to the end of the current element</span><span> </span><span style=' font-family: monospace;'>Racc(:,:,ic)</span><span> by left-multiplying Rcurr to </span><span style=' font-family: monospace;'>Racc(:,:,ic-1)</span><span> which already contains the transfer matrix from the start to the end of he previous segment. Finally we calculate the position </span><span style=' font-family: monospace;'>spos(ic)</span><span> of the current segment from that of the previous segment </span><span style=' font-family: monospace;'>spos(ic-1)</span><span> by adding the length of the current segment to it.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span >ic=1;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >line=1:nlines</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">for </span><span >seg=1:beamline(line,2)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > ic=ic+1;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Rcurr=eye(2);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">switch </span><span >beamline(line,1)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">case </span><span >1</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Rcurr=D(beamline(line,3));</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">case </span><span >2</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Rcurr=Q(beamline(line,4));</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">otherwise</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > disp(</span><span style="color: rgb(170, 4, 249);">'unknown element'</span><span >);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Racc(:,:,ic)=Rcurr*Racc(:,:,ic-1);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > spos(ic)=spos(ic-1)+beamline(line,3); </span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div class = 'S5'><span>At this point</span><span style=' font-family: monospace;'> Racc</span><span> contains all trasnfer matrices from the start of the beam line to the end of each and every segment while </span><span style=' font-family: monospace;'>spos</span><span> contains the positions of each segment.</span></div><div class = 'S0'><span>Now we can use the transfer matrices to plot the trajectory of a particle with starting coordinates </span><span texencoding="x_0 = 0" style="vertical-align:-6px"><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAFIAAAAoCAYAAABtla08AAAEeklEQVRoQ+3YV4gsRRTG8d81YsQMPqmYUUyIOWcFvWZFxBwQnxQTglkxPeiDiBgwIYgZc84ZTChmRH1QMIMJxcTn1qy9c7dnZnfby87S52WX6eqqU/8+dc53apbWGiEwq5FZ2km0IBsKghZkC7IhAg1N00ZkC7IhAg1N00bkWJDz4xAcjvmwCH7CNbgBf9Vxb0H+R2YZ3IXNsQ/uKY8OxE14ovz+y3gwW5AjVObFo9gOl+HELlhn4FzchgNakPX58Ehci9+wNH7uGro4vi3HfTc81D1VG5EjRD7EqngG29Twfgmb4MVy/McMa0GyDt4qVM7GOTUgL8Rp+Bsr4bPquF4g10PCOF/qVjxSeXEJnFCqWhZPZRtWOx5XFOf3xt01G0luDIfYYbixH8gFcD5+xK7YFO9irfJiqtur5avkp9m4d1gpohNp2cK2eLpmLztVgunyEkijQ/sd7S3wXBm9Oj7CzaWyRXNtiOvwa4MgXzFy3KZiyXnrDjhB9nNwGbs+3qx5byPEt9gd2K9fRFafz4MvsRyOLX+z0P0DOjmZYW8gaWUq9j7WHHCC6MPInlhy36c1762GD8qzx7HjREBmbETqXngd7+DQAR0clmHJ/Tm2sRXweY3jgfxJefZwSXujQ/sd7Qw8GZcgij4LfTMshAb0M0UjbWEsKeXtmvc2wGvlWdLZURONyE6e/BOL9ciHC+KsUpQiXpcvHyFRPJ3tYpxSHNyqUhO6fd4eOdKx83DmREEmfySPxAL1hXGoJJc+WORQnInW2r10C6mEqfrT1aryp5cC2Re3l02kXlw9EZBpjZ7FolgZp5Zj3g2ls0i6gnQHHQv0tF2dZD4IzLldtSPrOqcmgjsROp51+u08W6NSeP4d2y9HJn88gI1LI38f9hhnlYwJrAj1gOvYBTi9VNBU0kFsblft+BTfIu9SeHapcTL6cusuTT06tBfI/bEDjimiO9dK3yGCPEd3bUS8v1fEe2RRdGXV0incWVqrui89CNz/e0yO6lWloC6J37sWXBjfl/0egeu7HaqCXArRSjlaCffchiTKUq1zI/J1ieDkycC7CMdhRXyMxyoyorNO8uOTuBLJRdPVkuOfL11c/Iy/VTsJl5auJ0VnjgveKsi0eSkQ0VH5IjtXdFMmDahEaJ5HpCfavsCWJY/OofaRTiH6MwI+c09ni8pIikow7Vmp0PE7RSYnLmntq/E2UQWZI5weO9GV/7tlS6IrxzQQc3/3cplws1LJs1jSQdU6NytzCNhpSnQhHF1axsi9XPgmjaVWRDtW8/+YLfQrNoPsd5XSg+eGOVFctUihVPHklOSWGWtNgMxX/KGEfqp71aLLUqR63fPNCLhNgAyIyKIc/cifPypk0ukEYq9blRZkhcBBuAXJl7mS71i017KVu8wZAa1fsZnKJiMf0odGFqSyx3KjkkiNwH1qKpMPw7tNHe3sNaI1VT+tZGRRrp1y+1xtGYeByaR8bBLkpByYKS+1IBv6ki3IFmRDBBqapo3IFmRDBBqapo3IhkD+A9J8ziliFun/AAAAAElFTkSuQmCC" width="41" height="20" /></span><span> and </span><span texencoding="x_0'=10^{-3}" style="vertical-align:-10px"><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAHoAAAAyCAYAAACTUs/lAAAGiElEQVR4Xu2bdYhtVRSHv2cnBgYqFna32N2YGFjYIrZiF2JiIhao2IqJ3d2J3dgtKhhgN5+so3fOm3vfnjvnnLnzZq9/Zrh3n73XWr+z117rt/YdRZYR4YFRI8LKbCQZ6BHyEmSgM9AjxAMjxMy8o4c30NsCBwBzA58AZwDn92dSBnr4Ar0PsDRwKTAhcDCwErApcGPZrAz08AR6KuB0YKcW9acGvgYuK33+75AM9PAEegbge+CnFvXHjc+OA07OO3pgwI4P7AjsB8yf+Oh8wP7AUsAfwGTAY8CpwDuJc3QzbHngEmBl4IsMdJoLxwO2B44EZgN+ASZOeHQ74ALgiTgrvwOmBa6J83Rr4LaEeQY6ZFbgTmBz4I2mkjF3wQrAQwPVtkfGrwdsFODuGgCnAL06cC/wA7BgZMGFSVMC78Vc+uaFimx13mPiTDZyfAn4sqlHH6njjN4FOAGYviJjmp7GF/X3WPQOQODHBLQR4G1g9rDdSFAWz00z42eAZSoySvyMGDMCRosDga+AOYAfW9eoA2hDyPPAURUZM5TTeObtkAD0NsCVoeiqwMP9KL02cHd8vhZwX/zv3zUSjdwLOLfDWDfY4cCKwON1Aj1pJBwmJGaFw11SgbZu3QT4FTCcGgHKYmj9FnD3XwzsHAMEZbpER70UR0C74QsArwEbA7cMBOgpInRZmJu+G3pajVgnDLwp3lYXmBc4KVHxXh+WCrT16zTAK8AiHYx6M/zzYYT5qu03ZL8bFYJr/SedQve6gGHoN+DQAHoL4Pp4WurttPj/RWBx4HjgxFJ9V7UxTc6XAvREwM+h1CPAKh0UfBJYNr6X9DArr1JMxEwgjRJ9JPWMvh8wq7w6Dn2pNlN5QfWF8O15CmhNZAZjwAbAdYOZIJ7V8OLF7Ga6FKCLXeT8hkujWjsxf9Ffipn5690oBVgzGzU9Ai4H/gTmjP/NKUwMuwJ6T+CcKMQXBc4DNgP+6lLRMT1meXPzmAYlfL8lcG3CuHZDUoB29zwaE0g/6uh2chWwVXxpmWW93Y14PF4ILBTHheWULNlF7aJE6o5eGHg5NFK53QbxNnZj2FA9kwK0odiQrDi+lX8u6+2LYJRRLLEstRqRVKDHiTdlcuAU4JBGtBv6RVKAtnZ+P1Q1+7Z71E5uBTyWlJmBT5syMRVo9SnOac/lI5pScIjXSQFaarRoLjwYuUw7teW8Ddl/R2uxIGZqNzMVaMdJApiECfiatWvWGwukAK2mRXlVVB/ttLfGtda16SCb1ZikAm2TW5LfEkouV1LATK8uGU5Ztz4wszc51TeWTXatyuLlAAkTI4BNjiIpq8uHfeZNAdr2nBmezM/ngOf1EhUS8/0ZOpyybvU3u7f0VEzOnu7HKOvrotHTykf0BNATAA8EXWdtZmiyvNoXOCs0NIu0lhsbJTV0yx/YnTLBsqEh51wWLwT43QfAXDVHxNEW729Hm/bbIP8m6jLZHssC5UzAMF6c07uHgaO1xcYS1CVtJIZkBw29ncQul92uj4F5SlSxPLfUpJ0mmxv6r1EpA+0ZbOfJbPAj4IagPwulZGTsipg1Wju6249uVOPmFrOUfKslaZJ5ctd2kj2As+NynjdTinzGzpYdKzeGpEbjUgbaN0/wZgl6rb8yylrQ88ZwbQivMylr3CFBZNh2NCeZqUWBz4KWfTWIkXa6LRZ9YXMbOXA7epJNXuYrSKfG7UpJxgaqlLveJoiJm2SCTNqxEQUGOlceX5EHqgbanX4PsH401m1t2uyQ6tu7Ip3zNF14oGqgJQQkD2xvFiK5YrLm7cSC/O9C1fzIYDxQJdBeb302kjNDdSGTBE9uYtcoSTAYx4xtz1YJdHFfSVbr9pKjvHlh1iroWYbAA1UCfQXgj75Mxoq2XWGStbg8uXejDO1ZGvZAlUBblq3W5uaENy82BJaMOr1hM/NyVQJtwmXiZXem/GuB4pZko832DO//HqgSaNkfiYbloqRq9bN9WjNxfzoiRZilYQ9UCbSZtpf25XzvKtlhM0SmSFpR3jhLwx6oEmhbl88Bh5XuddsM8FqrREqnG5INmz6ylqsSaD1nK9M2nB2aQoo+bOM92JEFZWdrqwbaZMwa2r+yYM7vX8l9wbfrlWUIPFA10JpgiXVQ7Gwb8e5yG+7FrxmGwMy8ZB1AZ6/2oAcy0D0ISh0qZaDr8GoPzpmB7kFQ6lApA12HV3twzgx0D4JSh0oZ6Dq82oNz/gOrqDRCsivnugAAAABJRU5ErkJggg==" width="61" height="25" /></span><span> rad by mapping the initial (column) vector </span><span style=' font-family: monospace;'>x0 </span><span>with </span><span style=' font-family: monospace;'>Racc</span><span> to the ened of each segment and record both position</span><span style=' font-family: monospace;'> x(1)</span><span> and angle </span><span style=' font-family: monospace;'>x(2)</span><span> in the previously allocated array </span><span style=' font-family: monospace;'>data</span><span>. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span >data=zeros(nmat,2); </span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >x0=[0;1e-3];</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >k=1:nmat</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > x=Racc(:,:,k)*x0;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > data(k,1)=x(1);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > data(k,2)=x(2);</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div class = 'S5'><span>Now we can plot the trajectory and annotate the axes. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span >plot(spos,1000*data(:,1),</span><span style="color: rgb(170, 4, 249);">'k'</span><span >) </span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'s [m]'</span><span >)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'x [mm]'</span><span >) </span></span></div></div><div class="inlineWrapper outputs"><div class = 'S10'><span style="white-space: pre"><span >title(</span><span style="color: rgb(170, 4, 249);">'Trajectory'</span><span >)</span></span></div><div class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="D4409CDA" data-scroll-top="null" data-scroll-left="null" data-testid="output_4" style="width: 1364px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAYAAAAv7h+nAAAgAElEQVR4AezBfczddXn48Xe+lIoVCQ9zBo3rB6NcK8OtPkzdUv1eZz4MZjHSgWBFzvfWqRNbXfbgsix4OIv+ITKigcnQkPt7JFAVRhWFgGO9Px9GmOiwMtR5baS9CqglHfgAqKW0/HbimnX8aOnD/XTuc71e1RMhhBBCCCOmIoQQQghhxFSEEEIIIYyYihBCCCGEEVMRQgghhDBiKkIIIYQQRkxFCCGEEMKIqQghhBBCGDEVIYQQQggjpiKEEEIIYcRUhBBCCCGMmIoQQgghhBFTEUIIIYQwYipCCCGEEEZMRQhhLF188cWsWrWKVatWsWrVKlatWsWqVatYtWoVq1atotvtcrDe+ta3smrVKh5//HGm03333cd9991HCCFUhBDG0je+8Q3Wr1/P+vXrWb9+PevXr2f9+vWsX7+e9evXc+ONN3Kwrr/+etavX8+uXbuYLp/61Kd48YtfjJkRQggVIYSxdOWVV7J9+3a2b9/OZZddxtCZZ57J9u3b2b59O9///vc5WP/0T//ErbfeyqJFi5gu69evZ/v27YQQwlBFCGEsLVq0iMWLF7N48WIOP/xwhg477DAWL17M4sWLWbx4MStXruSMM87giiuu4Nhjj6XT6TD0ve99j7e85S0cddRRHHnkkSxfvpzJyUl2+/jHP87HPvYxdu3axdAjjzzCBz7wAZ773Ody9NFH87a3vY0tW7awp8cff5x+v8+JJ57IUUcdxfLly7niiisY6vf7fOtb32Lo/PPP5+KLL2a3T3/607ziFa/gqKOO4sQTT6Tf7/PYY4+x28qVKznjjDO44oorOPbYYznhhBNYuXIlN954I7v98Ic/ZOXKlbz1rW8lhDAaKkIIYS9uuOEGrr/+et773vfy2GOPkVLi8ccf53Wvex1f+tKX+L3f+z3+4A/+gO9+97u8853v5M4772Topptu4oYbbmDXrl0MnXbaaVxyySWklHj961/P5z73OV796lfz4IMPstu73vUuLrjgArZv384b3/hG7r33Xv7oj/6IK6+8ku9973s8/PDDDH3nO99h06ZNDPV6Pd773vfy7W9/G1XlJz/5CRdccAFvfvOb2e2GG27g+uuv573vfS+PPfYYz3ve87jhhhv41Kc+xW7XXHMNN9xwA89+9rMJIYyGihBC2IcdO3Zw8cUX88gjj3DZZZfx8MMPc9FFF3HZZZfxxS9+kS984QuceeaZDN1zzz082dTUFDlnXvrSl3LHHXdw7bXXcsEFF7B161YmJycZ2rJlC5/97Gd55jOfyTe/+U2uvfZaPvvZz/L7v//7PPDAA6xbt47XvOY1DF133XVceuml/PCHP+SjH/0ohx12GP/6r//K9ddfz/e+9z1e+MIXcvPNN/OVr3yF3Xbs2MHFF1/MI488wpe//GUOP/xwbrrpJn70ox8x9A//8A8MNU1DCGE0VIQQwtM466yzGDriiCM45phj+MM//EOOOeYY3vWud/GqV72Kq6++mr352te+xtAjjzzCu9/9bt797ndz2223MXTnnXcydOeddzL0+te/nuOOO46hlStXctNNN/Hnf/7nPJVbb72VnTt30ul0OPnkkxk65phjOO200xj60pe+xJ7OOussho499lje/va3s3PnTr74xS/ywAMPcOutt5JS4jWveQ0hhNFQEUIIT+O4445jtwcffBAR4eyzz+bee+/lzW9+M6rK3vzkJz9haMeOHTz44IM8+OCDPPvZz+b0009n+fLlDO3cuZOD9axnPYs9PetZz2Lo8ccfZ0/HHXccuzVNw9A111zD+vXrGTr33HMJIYyOihBCeBqLFi1itxtvvBF358wzz+Qf//Ef+eu//mue85znsDcvf/nLGXrRi17Eddddx3XXXccFF1zAueeey+rVqxl68YtfzNCGDRv4xS9+wdC3vvUtXvCCF/Dud7+bPe3atYuhZcuWMXTLLbfw4IMPstuGDRsYeu1rX8ueFi1axG51XfPCF76Qr371qwwGA4be/va3E0IYHRUhhHAAFi9ezNC//Mu/8MUvfpGPfexjXHPNNQw99thjPNlpp53G8573PG655Rb+5E/+hMFgwKmnnsrpp5/OXXfdxdDy5ct57Wtfy6OPPsqKFStYs2YNZ599Nvfffz/Pe97zGDr88MMZuvjii7n00kv5zd/8TU499VQeffRRXvOa13DeeedR1zVf+9rXEBHe/va3sy9N07Bz506+9rWv8epXv5oTTzyREMLoqAghhANw5plncvrpp3P//fdz+umnMxgM+Iu/+AuGNmzYwJMdccQR3HzzzbzkJS/hk5/8JE3T8JOf/IS//du/ZeXKlex27bXXcuqpp3LnnXfyd3/3d9xzzz28//3vp9frMfSWt7yFww47jJtvvpmbbrqJoS984Qu8//3v55577uGyyy7j1ltv5U1vehNTU1MsXryYfTn33HPZrdvtEkIYLRUhhLH3rne9iyeeeIJ169axpyeeeIInnniCPVVVxXXXXcfPf/5z/uu//ovvfve7XHjhhTzxxBNMTk6ya9cudu3axVBVVQydfPLJ/Nu//Rs///nP2bp1Kz/96U/50z/9U/b0nOc8hxtvvJHt27dz//3389hjj3HppZdSVRVD73nPe/jZz37G/fffz/XXX8/QkUceyaWXXsovfvELfvCDH7B9+3a+8pWvcPzxx7PbE088wRNPPMGTVVXFYYcdxuGHH85ZZ51FCGG0VIyphx56iA0bNnDHHXcQQjhwRxxxBMcddxx7uuWWW1i+fDk7duzgGc94BosWLWJPRxxxBM997nOpqoq9Wbx4Mc9//vOpqoonW7x4Mc9//vOpqoo9VVXF8ccfz+LFi3k63/jGNzj//PN505vexM6dOznrrLM45phjCCGMlooxVEph5cqV3HjjjVx00UWcc8457Nq1ixDCofn85z/P3XffzWGHHcZFF13EfPTAAw/wkY98hLvvvpuXvvSlXHjhhYQQRk/FmNm5cyd/9Vd/xSc+8QkuuugirrnmGn784x/z1a9+lRDCofnMZz7Dww8/zM9+9jPWrFnDfLRy5UoeffRRfv7zn/PNb36T448/nhDC6KkYM6UUnv/85/PKV76S3b7yla9wyimnEEI4dEceeSSLFy9mPluyZAlHHHEEIYTRVTFmfvSjH/GCF7yAD3/4w/zWb/0WL3vZy7jiiisIIYQQwuioGDP33HMPN998M7/xG7/BXXfdxbp16/j7v/97brvtNvbmHe94ByKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICCKCiCAiiAgigoggIogIIoKIICKICO94xztYiCrGzK/92q+xdOlSzjrrLIZEhDe84Q3ceOON7M3Xv/51zAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzBgyM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzY8jMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8yMITPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDOGzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzhswMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wYMjPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM2PIzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMjCEzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzvv71r7MQVYyZ4447jierqoqqqgghhBDCaKgYM51Oh4ceeoipqSmGHnroIf75n/+Z0047jRBCCCGMhooxc/jhh3PppZfyN3/zN5x99tm88Y1v5KyzzuJVr3oVC9WaNWsYZWvWrGGUrVmzhlG2Zs0aRt2aNWsYZWvWrGGUrVmzhlG2Zs0awvxTMYZe8YpXMDU1xeTkJHfccQfnnXceC9natWsZZWvXrmWUrV27llG2du1aRt3atWsZZWvXrmWUrV27llG2du1awvxTMcae+cxncthhhxFCCCGE0VIRQgghhDBiKkIIIYQQRkxFCCGEEMKIqQghhBBCGDEVIYQQQggjpiKEEEIIYcRUhBBCCCGMmIoQQgghhBFTEUIIIYQwYipCCCGMvE6nw8TEBCGMi4oQQggjz91p25YQxkVFCCGEkebuDKWUcHdCGAcVIewnd2diYgJ3J4Qwf+ScSSmhquScCWEcVISwn9ydtm3JORNCmD8GgwGqSl3XlFIIYRxUhLCf+v0+Q6UUQgjzh7tT1zVN09C2LSGMg4oQ9oO7k3OmaRpCCPOHuzOkqgypKm3bEsJCVxHCfnB3UkpMTk7Sti0hhPnB3UkpsZuqUkohhIWuIoT9MBgMaJqGoZQS7k4IYe71+31Uld2WLl1K27aEsNBVhLAf3J26rhlSVXLOhBDmh7qu2a1pGlJKuDshLGQVITyNnDPujqoyVNc1pRRCCHPL3XF3VJU9qSo5Z0JYyCpCeBqlFFSV3ZqmoW1bQghzy91JKfFkdV1TSiGEhawihKeRc6bb7bKnlBLuTghh7vT7fVSVJ2uahrZtCWEhqwhhH9wdd0dV2ZOqknMmhDC36rrmqaSUyDkTwkJVEcI+DAYDUko8WV3XlFIIIcwNd8fdUVWeiqoyGAwIYaGqCGEfcs70ej2eLKVE27aEEOaGu5NSYm/quiaEhawihL1wd9wdVeXJVJUhdyeEMPtKKagqe9M0DW3b4u6EsBBVhLAXOWdSSuxN0zTknAkhzL6cM3Vdsy+qSs6ZEBaiihD2YjAY0O122ZuUEqUUQgizy91xd1SVfel2u5RSCGEhqgjhKbg7OWeapmFvli5dirsTQphd7k5KiaejqrRtSwgLUUUIT8HdUVX2pWkacs6EEGZXKQVV5emklEgp4e6EsNBUhPAU+v0+qsrTSSmRcyaEMHtyztR1zf5QVXLOhLDQVITwJO5Ozpm6rnk6qkophRDC7HB33B1VZX/UdU0phRAWmooQnsTdSSmhqjyduq5xd0IIs8PdSSmxv5qmoW1bQlhoKkJ4klIKTdOwP5qmoW1bQgizo5SCqnIgUkq0bUsIC0lFCE+Sc6aua/ZXSgl3J4Qw83LO1HXNgVBVSimEsJBUhLAHd8fdUVX2l6qScyaEMLPcHXdHVTkQdV2TcyaEhaQihD0MBgNUlQNR1zWlFEIIM8vdSSlxoJqmwd1xd0JYKCpC2EPOmW63y4Fomoa2bQkhzKxSCqrKwWiahpwzISwUFSH8D3fH3VFVDlRKCXcnhDBzcs7Udc3BqOuaUgohLBQVIfyPnDMpJQ6GqpJzJoQwM9wdd0dVORhN09C2LSEsFBUh/I/BYECv1+Ng1HVNKYUQwsxwd1JKHIqUEu5OCAtBRQj/zd1xd1SVg5VzJoQwM0opqCqHQlXJORPCQlARwn9zd1JKHKymaXB33J0QwvTLOVPXNYeirmtKKYSwEFSE8N/6/T6qyqFQVXLOhBCml7vj7qgqh6JpGtq2JYSFoCKMPXcn50xd1xyKlBKlFEII08vdSSkxHVSVtm0JYdRVjLm77rqLbdu2Mc7cnZQSqsqhqOuaEML0K6WgqkyHbrdLKYUQRl3FGLvnnns455xzuOuuuxhnpRSapuFQNU1D27aEEKZXzpm6rpkubdsSwqirGFM7duzgz/7sz/iVX/kVxl3OmbqumQ4pJdydEML0cHfcHVVlOjRNQ0oJdyeEUVYxpi6++GJe97rXceKJJzLO3B13R1WZDqpKzpkQwvRwd1JKTCdVJedMCKOsYgx9/etf54477uADH/gA424wGKCqTJe6rimlEEKYHqUUVJXpVNc1pRRCGGUVY+anP/0pH/7wh7n44os5ECKCiCAiXHLJJSwUOWe63S7TpWka2rYlhDA9cs7Udc10apqGtm0JC9Mll1yCiCAiiAgLVcWYufDCCznppJPYsmULpRQeeughvvOd72Bm7IuZYWaYGWvXrmUm5Zxp25aZ5u64O6rKdEop4e6EEA6Nu5NzRlWZbiklcs6EhWft2rWYGWaGmbFQVYyZ5zznOTz66KNcffXVXH311Xz/+9+nlMLtt9/OfNHv9+n3++ScmUk5Z1JKTDdVJedMCOHQuDspJWaCqjIYDAhhVFWMmQ9+8INcfvnlXH755Vx++eW85CUv4bzzzmNiYoL5wN1xd6ampuh0OuScmSmDwYBer8d0q+uaUgohhENTSqFpGmZCXdeEMMoqwrzi7qSUSCkxNTXFxMQEM8HdcXdUlemWUqJtW0IIhybnTF3XzISmaWjbFncnhFFUMeYuv/xyXv/61zNflFJQVYZUlcnJSTqdDtPN3UkpMRNUlSF3J4RwcNydnDOqykxpmoacMyGMooowr+Scqeua3VQVVaXT6TCd+v0+3W6XmdI0DTlnQggHx91JKTGT6rqmlEIIo6gizCs5Z1JK7KnX6zHU7/eZDu5OzpmUEjMlpUQphRDCwSml0DQNM6lpGtq2JYRRVBHmDXcnpURKiSebmpoi50zOmUPl7qgqqspMWbp0Ke5OCOHg5Jyp65qZllLC3Qlh1FSEeWMwGJBSYm8mJyfpdDrknDkUg8EAVWUmNU1DzpkQwoFzd3LOpJSYaapKzpkQRk1FmDdyzqgqe5NSYmpqiomJCQ5Fzpm6rplpKSXatiWEcGDcnZQSKSVmWl3XlFIIYdRUhHmlrmv2RVWZnJyk0+lwMHLODKkqM01V2bJlCyGEA1NKQVWZDU3T0LYtIYyaijAvuDvujqrydFQVVaXT6XCgSik0TcNsqOsadyeEcGByznS7XWaLqtK2LSGMkoowL7g7KSX2V6/XY6jf73Mgcs7Udc1saJqGtm0JIew/dyfnTEqJ2aKqlFIIYZRUhHmhlIKqciCmpqbIOZNzZn+4O+6OqjJbUkq4OyGE/ePupJRIKTFbli5dStu2hDBKKsK8kHOmrmsO1OTkJJ1Oh5wzT2cwGKCqzCZVJedMCGH/lFJQVWZT0zSklHB3QhgVFWHOuTs5Z1SVA5VSYmpqiomJCZ5Ozplut8tsquuaUgohhP2Tc6bb7TLbVJWcM9PN3ck5E8J0qwhzzt1JKXGwVJXJyUk6nQ574+64O6rKbGqahrZtCSHsn5wzc6Gua0opTBd3p9/v0+l06HQ6uDshTKeKMOdKKTRNw6FQVVJKdDodnkrOmZQScyGlhLsTQtg3dyelhKoy25qmoW1bpoO7MzExQdu2TE1N0TQNOWdCmE4VYc7lnKnrmkM1OTnJUL/f58kGgwHdbpe5oKrknAkh7NtgMCClxFxJKZFz5mC5O/1+n06ng6qyefNmUkrUdU0phRCmU0WYczlnpsvU1BQ5Z3LO7ObuuDtN0zAX6rqmlEIIYd9yzvR6PeaKqlJK4WC0bcvExARDmzdvptfrsVvTNLRtSwjTqSLMKXcnpYSqMl0mJyfpdDrknBlyd1JKzKW2bQkh7FvOmblU1zXuzoFwdyYmJuj3+/R6PXq9Hk8lpYS7E8J0qQhzyt1JKTGdUkpMTU0xMTHBUL/fR1WZK03TMOTuhBCemruTUkJVmStN09C2Le7O03F3+v0+nU6HlBKbN29GVdkbVSXnTAjTpSLMqX6/j6oy3VSVyclJOp0OOWfqumYuqSo5Z0IIT20wGJBSYq6pKjln9sXd6XQ6uDtTU1P0ej2eTl3XlFIIYbpUhDlX1zUzQVVJKaGqqCpzSVUppRBCeGo5Z3q9HnOt2+1SSuGpuDv9fp9Op0Ov12NycpKUEvujaRratiWE6VIR5oy74+6oKjNlcnKSqakp5trSpUsJIeyduzMfqCpt2/Jk/X6fTqfD0ObNm2mahgOVUsLdCWE6VIQ54+6klBgHTdPQti0hhP+fuzOkqsy1lBIpJdydIXen0+nQti1TU1P0ej0OlqqScyaE6VAR5kwpBVVlXKSUyDkTQvi/cs6klJgvVJXBYEC/36fT6aCqbN68mZQSh6Kua0ophDAdKsKcyTlT1zXjQlVxd0II/9dgMEBVmS+63S4XXHABQ5s3b6bX6zFd2rYlhOlQEeaEu5NzRlUZF3VdU0ohhPB/uTt1XTNfqCqbN2+m1+sxnZqmYcjdCeFQVYQ54e6klBgnTdPQti0hhP/l7gypKvNJSomZ0DQNOWdCOFQVYU6UUkgpMW5SSrg7IYRfyjmTUmJcpJQopRDCoaoIcyLnTK/XY9yoKjlnQgi/NBgMUFXGxdKlS3F3QjhUFWFOuDvjqK5rSimEEH7J3anrmnHRNA05Z0I4VBVh1rk7Q6rKuEkp0bYtIQRwd4ZUlXGSUqJtW0I4FBVh1rk7KSXGkaqSUsLdCWHcuTspJcaNqrJlyxZCOBQVYdaVUlBVxpWqknMmhHHX7/dRVcZNXde4OyEcioow63LO1HXNuKrrmlIKIYwbd6dtWyYmJuh0Ogx1u13GTdM0tG1LCIeiIswqd8fdUVXGmbsTwkLn7rRty8TEBJ1Oh06nQymFuq7p9XpMTU2RUmIcpZRwd0I4WBVhVrk7KSXGWdM05Jxxd0JYSNydtm2ZmJig0+nQ6XQopVDXNb1ej82bNzM5OUnTNKgq40xVyTkTwsGqCLOqlIKqMu5SSuScCWGUuTtt2zIxMUGn06HT6TAYDEgp0ev12Lx5M5OTkzRNg6oS/ldd15RSCOFgVYRZlXOmrmvGnapSSiGEUdO2LRMTE3Q6HTqdDoPBgJQSvV6PzZs3MzU1Ra/XQ1UJe9c0DW3bEsLBqgizKudMSolxV9c1IYyanDP9fp+UEr1ej82bNzM1NUWv10NVCQcmpYS7E8LBqAizxt1JKZFSYtw1TUPbtoQwSkopqCq9Xg9VJRwaVSXnTAgHoyLMmsFgQEqJ8EspJXLOhDAqcs50u13C9KjrmlIKIRyMijBrcs6oKuGXpqammJiYIOdMCPOdu5NzRlUJ06NpGtq2JYSDURFmVV3XhF9KKTE1NUWn0yHnTAjzmbuTUiJMr5QS7k4IB6oizAp3x91RVcL/SimxefNm+v0+/X6fEOarUgpN0xCml6qScyaEA1URZoW7k1Ii/P9SSkxNTZFzptPpEMJ8lHOmrmvC9KrrmlIKIRyoijArSimoKmHvpqamUFU6nQ4hzCfujrujqoTp17YtIRyoijArcs7UdU3Yt16vR7fb5YQTTiDnTAjzgbuTUiJMv6ZpGHJ3QjgQFWHGuTs5Z1SV8PSapmFycpJOp0POmRDmWikFVSXMDFUl50wIB6JiTN1zzz3ccsstfPOb32SmuTspJcL+U1U2b95Mv9+n3+8TwlzKOVPXNWFmqCqlFEI4EBVj6CMf+Qjvec97uPnmm+n3+6xevZrt27czU0opNE1DODApJSYnJ8k5MzExQQhzwd1xd1SVMDOWLl2KuxPCgagYM//+7//O5z//ea677uFgr0wAACAASURBVDo+/vGP86UvfYmHH36YL3/5y8yUnDN1XRMOXEqJqakphjqdDiHMNncnpUSYOU3TkHMmhANRMWaOPvpoLr/8co4++mh2O+GEE/jBD37ATMk5Ew7N5OQkqsoJJ5yAuxPCbBkMBqgqYWallGjblhD2V8WYOf744/nd3/1ddtuyZQtTU1O84Q1vYCa4OyklVJVwaHq9HpOTk3Q6HXLOhDAb3J26rgkzS1XZsmULIeyvijH2wAMP0DQN5513HsuWLWNfRAQRQUS45JJL2F/uTkqJMD1UlampKTqdDjlnQphJ7o67o6qEmVXXNe5OOHSXXHIJIoKIICIsVBVj6u677+b000/n3HPP5X3vex9Px8wwM8yMtWvXsr/6/T6qSpg+KSU2b95Mv9+n3+8TwkzJOZNSIsy8pmlo25Zw6NauXYuZYWaYGQtVxRi6/fbbeec738kFF1zAxMQEM62ua8L0SikxNTVFzplOp0MIM2EwGKCqhNmRUsLdCWF/VIyZ++67jzVr1nDhhRfS6XTYsWMHO3bsYOfOnUw3d8fdUVXCzJiamkJVOeGEEwhhurk7dV0TZoeqknMmhP1RMWauvvpqHn30Uf74j/+Yk08+mZNPPpmTTz6Zj370o0w3dyelRJhZvV6Ppmno9/uEMF3cnSFVJcyOuq4ppRDC/qgYM3/5l3+JmWFmmBlmhpnx4Q9/mOlWSkFVCTOv1+txwQUXEMJ0GQwGpJQIs6dpGtq2JYT9URFmTM6Zuq4JsyOlRM6ZEKZDzpler0eYXSkl3J0Qnk5FmBHuTs4ZVSXMDlVlMBgQwqFyd3LOpJQIs0tVyTkTwtOpCDPC3UkpEWZPXde0bUsIh8rdSSmRUiLMrrquKaUQwtOpCDOilEJKiTB7mqYhpYS7E8KhKKWgqoTZp6q0bUsIT6cizIicM71ejzC7VJWcMyEcipwz3W6XMPtSSqSUcHdC2JeKMCPcnTD76rqmlEIIB8vdcXdUlTA3VJWcMyHsS0WYdu7OkKoSZlfTNLRtSwgHy91JKRHmTl3XlFIIYV8qwrRzd1JKhLmRUqJtW0I4GKUUVJUwt9q2JYR9qQjTrpSCqhLmRtM0lFII4WDknKnrmjB3mqZhyN0JYW8qwrTLOVPXNWFuLF26lLZtCeFAuTvujqoS5paqknMmhL2pCNPK3XF3VJUwN5qmIaWEuxPCgXB3UkqEuaeqlFIIYW8qwrRyd1JKhLmlquScCeFAlFJQVcLcW7p0Ke5OCHtTEaZVKQVVJcytuq4ppRDCgcg5U9c1Ye41TUPOmRD2piJMq5wzdV0T5lbTNLRtSwj7y91xd1SVMD+klGjblhCeSkWYVjlnUkqEuaeqtG1LCPvD3UkpEeYPVWXLli2E8FQqwrRxd1JKpJQIc6/b7VJKIYT90e/3UVXC/FHXNe5OCE+lIkybwWBASokwf7RtSwj7q65rwvzRNA1t2xLCU6kI0ybnjKoS5oemaUgp4e6EsC/ujrujqoT5JaWEuxPCk1WEaVXXNWH+UFVyzoSwLzlnUkqE+UdVyTkTwpNVhGnh7rg7qkqYP+q6ppRCCPsyGAzo9XqE+aeua0ophPBkFWFauDspJcL80jQNbdsSwr7knAnzU9M0tG1LCE9WEaZFKQVVJcw/qkrbtoTwVNydlBKqSpifUkq4OyHsqSJMi5wzdV0T5p9ut0sphRCeymAwQFUJ85eqknMmhD1VhEPm7rg7qkqYf1SVtm0J4anknOl2u4T5q65rSimEsKeKcMhyzqSUCPNTSomUEu5OCHtyd3LOpJQI81dKibZtCWFPFeGQDQYDer0eYf5SVXLOhLAndyelREqJMH+pKikl3J0QdqsIh8TdcXdUlTB/1XVNKYUQ9lRKoWkawvynquScCWG3inBIcs6klAjzW9M0tG1LCHvKOVPXNWH+q+uaUgoh7FYRDslgMKDX6xHmP1WlbVtCGHJ33B1VJYyGtm0JYbeKeWrbtm1s27aNbdu2sW3bNrZt28a2bdvYtm0b27ZtY9u2bWzbto1t27axbds2tm3bxrZt25hN7o67o6qE+a/b7VJKIYQhdyelRBgNTdMw5O6EMFQxT73tbW9jxYoVrFixghUrVrBixQpWrFjBihUrWLFiBStWrGDFihWsWLGCFStWsGLFCuq6ZjblnEkpEUZD0zS0bUsIQ6UUVJUwOlSVnDMhDFXMY5/73OfYuHEjGzduZOPGjWzcuJGNGzeyceNGNm7cyMaNG9m4cSMbN27ktttuY7YNBgN6vR5hdKSUcHdCyDlT1zVhdKgqpRRCGKqYp0444QSOPPJIlixZwpIlS1iyZAlLlixhyZIlLFmyhCVLlrBkyRKWLFnCkiVLeNaznkVKidni7rg7qkoYHapKzpkw3twdd0dVCaNj6dKlhLBbxTz1mc98hhe/+MXsdvfdd3POOeewevVqVq9ezerVq1m9ejXnnHMOQ0uWLOHGG29ktgwGA1SVMFrquqaUQhhv7k5KiTBamqahbVtCGKoYAR/84Ac544wzuO+++9i6dStbt25l69atbN26la1btzIXcs50u13CaGmahrZtcXfC+CqloKqE0ZNSom1bQqgYAaUUPv3pT1NKYcOGDWzYsIENGzawYcMGbrnlFmabu+PuqCph9KgqOWfC+Mo5U9c1YfSoKqUUQqgYAYsWLeKkk05ivhgMBqgqYTR1u11KKYTx5O64O6pKGD11XdO2LSFUjIDzzz+fD33oQ+zcuZP5IOdMt9sljKamaWjbljCecs6klAijqWkaUkq4O2G8VYyAU045hdtvv52TTjqJ5cuXs3z5cpYvX87y5ct52ctexmxyd9wdVSWMrpQS7k4YP4PBAFUljC5VJedMGG8VI+CUU07hmGOO4corr2TdunWsW7eOdevWsW7dOq666ipm02AwQFUJo01VGQwGhPHj7tR1TRhddV1TSiGMt4oR8NBDD7F+/Xpe+cpXsmzZMpYtW8ayZctYtmwZy5YtYzblnOl2u4TRVtc17k4YL+7OkKoSRlfTNLRtSxhvFSPguOOO495772WuuTvujqoSRlvTNLRti7sTxsdgMCClRBh9qkrbtoTxVTECrrrqKs4991wuu+wyNm3axKZNm9i0aRObNm1i06ZNzJbBYICqEhaGpmnIORPGR86ZXq9HGH3dbpdSCmF8VYyAs88+m6FPfOITnHrqqZx66qmceuqpnHrqqaxcuZLZknOm2+0SFoa6rimlEMaDu5NzJqVEGH1N09C2LWF8VYyAUgpmhplhZpgZZoaZ8d3vfpfZ4O64O6pKWBiapqFtW8J4cHdSSqSUCAtDSgl3J4ynirBfBoMBqkpYWFJKuDth4SuloKqEhUNVGQwGhPFUMSLuvvtuzjnnHFavXs3q1atZvXo1q1ev5pxzzmE25JzpdruEhUVVGQwGhIUv50y32yUsHHVd4+6E8VQxAj74wQ9yxhlncN9997F161a2bt3K1q1b2bp1K1u3buVg3Hfffdxyyy2YGU9nx44duDuqSlhY6rrG3QkLm7uTc0ZVCQtH0zS0bYu7E8ZPxQgopfDpT3+aUgobNmxgw4YNbNiwgQ0bNnDLLbdwoL785S9z9tlnc/PNN/O+972PT37yk+zLT3/6U1SVsPA0TUPbtrg7YeFyd1SVsPA0TUPOmTB+KkbAokWLOOmkk5gOO3fupNfrMRgM+PjHP861117L5OQk7s7e/OxnP6Pb7RIWpqZpyDkTFq5SCqpKWHjquqaUQhg/FSPg/PPP50Mf+hA7d+7kUN16660cffTRvOhFL2Lo2GOP5bWvfS233XYbT8Xdefzxx1FVwsJU1zWlFMLClXOmrmvCwtM0DW3bEsZPxQg45ZRTuP322znppJNYvnw5y5cvZ/ny5SxfvpyXvexlHIgf//jH/Pqv/zp7OvLII/mP//gP9ua4445DRBARRIRLLrmEsHA0TUPbtoSFyd1xd1SVsDCllMg5E37pkksuQUQQEUSEhapiBJxyyikcc8wxXHnllaxbt45169axbt061q1bx1VXXcWB2LlzJ1VVsaeqqti1axdPJaXEUUcdhZlhZpgZa9euJSwsKSVyzoSFx91JKREWLlVlMBgQfmnt2rWYGWaGmbFQVYyAhx56iPXr1/PKV76SZcuWsWzZMpYtW8ayZctYtmwZB+IZz3gGO3fuZE+7du1i0aJFhPGlqgwGA8LCU0pBVQkLV13XtG1LGC8VI+C4447j3nvvZTr86q/+Kt/+9rfZ049+9CNe/vKXE8ZXXde4O2HhyTlT1zVh4WqahpQS7k4YHxUj4KqrruLcc8/lsssuY9OmTWzatIlNmzaxadMmNm3axIH47d/+bYZKKQz953/+J7fffju/8zu/QxhfTdOQc8bdCQuHu+PuqCphYVNVcs6E8VH9v/bgBkar+s77/7sHvVVidomLUtl0OVbbr3Vx2eCqkVDPbxRc2EhYm10L1HXOWDMVqjWbPkUwjFOthgetyKzsaO1ckzYQu6yuYdMsLjrnCExWejmJItN8HCI/JG1K2h0JjV10crF3Jv/wz9zWB5Cn65r5vl40gAULFjDs0UcfZe7cucydO5e5c+cyd+5cbrzxRo5FkiSsXr2apUuX0tzczMKFC1mxYgUTJ07EjW15nlMUBW70iDGSpilu9MuyjLIscWNHQgMoyxJJSEISkpCEJPr7+zlWV199Ndu3b2fdunXs2LGDOXPm4FyWZZRliRs9yrIkhIAb/fI8p1Kp4MaOhDo1b948BgYGOFq///3vmT59Osdi/PjxJEmCc8PyPKdSqeBGj6IoyLIMNzaEEKhUKrixIaFO/c///A+SGBgYYGBggIGBAQYGBhgYGGBgYICBgQEGBgYYGBhgYGAASRw6dAjnjkeaphRFgWt8MUZijIQQcGNDc3MzZVnixoaEOpUkCd/5zneYP38+8+fPZ/78+cyfP5/58+czf/585s+fz/z585k/fz7z58/nK1/5Cp/61Kdw7niEEOju7sY1vhgjaZrixo48z6lUKrixIaFOPf/88/T399Pf309/fz/9/f309/fT399Pf38//f399Pf309/fT39/P/39/ezatQvnjkeWZRRFgWt87e3thBBwY0uapsQYcaNfgnPu/5fnOTFGYoy4xpdlGW5sCSHQ3t6OG/0SnHP/jzzPKYoC17hijMQYCSHgxpYsy4gx4ka/BOfc/yPLMsqyxDWuoihI0xQ39uR5TlEUxBhxo1tCg6vVajh3IuV5TqVSwTWu7u5uQgi4sSnPc4qiwI1uCQ2gqamJnTt38n49PT1cfvnlOHeipWlKjBHXmGKMZFmGG5uyLKMsS9zoltAArrnmGv7u7/6ONWvWcMTdd9/NHXfcwVe+8hWcO9FCCBRFgWs8MUaGhRBwY1Oe51QqFdzoltAAHnzwQZ566ikef/xxsixj2rRpbN++neeff55ly5bh3ImWZRllWeIaT3d3N2ma4sa2NE2pVCq40SuhQcycOZPvfOc7/PrXv+bQoUN873vfY8qUKTh3MuR5TqVSwTWeoihoa2vDjW1tbW2UZYkbvRIawODgIFmWsXLlStasWcN9993HP/7jP9Lc3IxzJ0uapsQYcY2lKAqcG1apVHCjV0IDuOmmm5gwYQKvvfYac+bMYeHChfz85z9n9+7dXHbZZTh3MoQQKIoC1zhijKRpSggBN7bleU6apsQYcaNTQgP41re+xXPPPcdZZ53FEX/0R3/E9u3baWlpwbmTIcsyyrLENY7u7m5CCDg3LIRAURS40SmhAcybN48P8+1vfxvnToYQApVKBdc4iqKgubkZ54ZlWUZZlrjRKcE594HSNCVNU2KMuPoXY6QoCtI0xblheZ5TqVSIMeJGnwTn3IcKIVAUBa7+xRhJ05Q0TXHuiDzPKYoCN/okOOc+VJZllGWJq39lWRJCwLmRsiyjLEvc6JPgnPtIlUoFV/+KoqC5uRnnRsrznEqlght9EpxzHyrPc4bFGHH1K8ZIURSEEHDu/dI0pSgK3OiS4Jz7SHmeUxQFrn7FGEnTFOc+SAiB7u5u3OiS4Jz7SGmaUpYlrn6VZUme5zj3QbIso1Kp4EaXBOfcR5oyZQoxRlz9KoqCLMtw7oPkeU6apsQYcaNHgnPuI+V5TlEUxBhx9SfGSIyREALOfZgQAkVR4EaPBOfcx0rTlKIocPUnxkiapjj3UbIsoyxL3OiR4Jz7WCEEyrLE1Z+yLAkh4NxHyfOcSqWCGz0SnHMfK8syXH0qioIsy3Du44QQqFQquNEhwTn3sfI8p1Kp4OpLjJEYIyEEnPs4zc3NlGWJGx0SnHNHJU1TiqLA1Y8YI2ma4tzRCCFQqVRwo0OCc+6ohBAoyxJXP8qyJISAc0cjTVPSNCXGiGt8Cc65o5JlGTFGXP0oioIsy3DuaIUQKIoC1/gSnHNHJc9zKpUKrj7EGIkxEkLAuaOVZRllWeIaX4Jz7qilaUqMEXf6xRhJ0xTnjkWe51QqFWKMuMaW4Jw7aiEEiqLAnX5lWRJCwLljlec5RVHgGluCc+6oZVlGWZa4068oCrIsw7ljlWUZZVniGluCc+6o5XlOpVLBnV4xRmKMhBBw7ljleU6lUsE1tgTn3DFJ05QYI+70iTGSpinOfVJpmlIUBa5xJTjnjkkIgaIocKdPe3s7IQSc+6RCCHR3d+MaV4Jz7phkWUZZlrjTK8synPuksiyjUqngGleCc+6Y5HlOpVLBnR4xRmKMhBBw7pPK85w0TYkx4hpTgnPumKVpSowRd+oVRUGapjh3vEIIFEWBa0wJzrljFkKgKArcqdfd3U0IAeeOV1tbG+3t7RRFgWs8Cc65Y5ZlGWVZ4k69GCNZluHc8UrTlK6uLlpaWiiKAtdYEpxzxyxNUyqVCu7UijEyLISAcydCCIGuri5aWlooigLXOBKcc8cshECapsQYcadOd3c3aZri3IkUQqCrq4uWlhaKosA1hoQxavfu3WzZsoW+vj6c+yRCCBRFgTt1iqKgra0N5060EAJdXV20tLRQFAWu/iWMQQ888ACtra1s3ryZ9vZ2Fi1axLvvvotzxyLLMsqyxJ06RVHg3MkSQqCrq4uWlhaKosDVt4Qx5he/+AVPP/00zzzzDKtWreK5557jd7/7HZs2bcK5Y1WpVHCnRoyRNE0JIeDcyRJCoKuri6amJoqiwNWvhDFmwoQJdHZ2MmHCBI646KKL+NWvfoVzxyLPc4bFGHEnX3d3NyEEnDvZQgjs2bOHpqYmiqLA1aeEMebCCy9kxowZHLF37156enqYPXs2H8XMMDPMjLVr1+LcsBACRVHgTr6iKGhubsa5UyFNU/bs2UNTUxNFUdBI1q5di5lhZpgZo1XCGLZ//37yPGfJkiV84Qtf4KNIQhKSuOuuu3BuWAiBsixxJ1eMkaIoSNMU506VNE3Zs2cPLS0tNJK77roLSUhCEqNVwij3wAMPMH36dKZPn84Xv/hFjti5cyc33XQTt956K4sXL8a5T2LKlCnEGHEnV4yRNE1J0xTnTqU0TdmzZw+u/iSMcosWLaKjo4OOjg4efvhhhvX29nLbbbdx33330dLSgnOfVJ7nFEWBO7nKsiSEgHPOHZEwyn32s59lxowZzJgxg6uuuop9+/Zx5513snLlSpqamhgaGmJoaIharYZzn0SaplQqFdzJUxQFzc3NOOfcEQljzPr163nnnXe44447mDp1KlOnTmXq1Kl8//vfx7lPIoRAWZa4kyPGSFEUhBBwzrkjEsaY7373u0hCEpKQhCSWL1+Oc59ElmW4kyfGSJqmOOfcSAnOueOS5zmVSgV3cpRlSZ7nOOfcSAnOueOWpilFUeBOvKIoyLIM55wbKcE5d9xCCMQYcSdWjJEYIyEEnHNupATn3HHLsoyyLHEnVoyRNE1xzrn3S3DOHbc8z6lUKrgTqyxLQgg459z7JTjnTog0TYkx4k6coijIsgznnHu/BOfcCRFCoCgK3IkRYyTGSAgB55x7vwTn3AmRZRllWeJOjBgjaZrinHMfJME5d0LkeU6lUsGdGGVZEkLAOec+SIJz7oRJ05QYI+74FUVBlmU459wHSXDOnTAhBIqiwB2fGCMxRkIIOOfcB0lwzp0wWZZRliXu+MQYSdMU55z7MAnOuRMmhEClUsEdmxgjMUaKoqBSqdDe3k4IAeec+zAJzrkTJk1T0jQlxoj7/8QYiTFSFAWVSoX29nZaWlpoamqiqamJT33qUzQ1NdHU1ER7eztlWRJCoK2tDeec+zAJzrkTKoRAURSMNUVRUKlUaGlpoaWlhaamJi666CIuuugimpqaaG9vpyxLYoxkWUZbWxttbW387//+L3v27GHPnj309PTQ1dVFW1sbzjn3URKccydUlmWUZclYEmOkpaWFsizJsowsy2hra6Onp4f//d//Zc+ePfT09NDV1UVXVxd5nhNCIISAc859EgnOuRMqTVMqlQpjSYyRNE3p6uoiz3PyPCeEQJqmOOfcyZDgnDuhQggMizEyVpRlSQgB55w7VRKccydcnucURcFYURQFWZbhnHOnSoJz7oRL05SyLBkLYozEGAkh4Jxzp0qCc+6EmzJlCjFGxoIYI2ma4pxzp1KCc+6Ey/OcoiiIMXIiVCoVmpqa+NSnPkVLSwv1pCxLQgg459yplOCcOylCCBRFwfGoVCo0NTXR3t5Oc3MzXV1dVCoV6klRFGRZhnPOnUoJzrmTIk1TyrLkWMUYqVQqNDU10d7eTnNzM3v27CHPc/I8J01TYozUgxgjMUZCCDjn3KmU4Jw7KbIs41jEGKlUKjQ1NdHe3k5bWxt79uwhz3NGCiFQFAX1IMZImqY459ypluCcOynyPKdSqfBxYoy0t7fT1NREd3c3XV1d7NmzhxACHyTLMsqypB6UZUkIAeecO9USnHMnTZqmFEXBB4kx0t7eTlNTEzFGurq66OnpIYTAR8nznEqlQj0oioIsy3DOuVMtwTl30oQQKMuSkWKMtLe309TURIyRrq4uurq6CCFwtNI0JcbI6RRjJMZICAHnnDvVEpxzJ02WZcQYGRZjpL29naamJmKM9PT00NXVRQiBYxVCoLu7m9Mpxkiapjjn3OmQ4Jw7afI8p1Kp0N7eTlNTE8N6enro6uoiTVM+qSzLiDFyOpVlSQgB55w7HRKccyfVfffdx7Cenh7a2tpI05Tjlec5lUqF06koCrIswznnTocE59xJ1dbWRltbG2maciKFEKhUKpwOMUZijIQQcM650yHBOdeQQgiUZcnpEGMkTVOcc+50SXDONaQpU6ZQqVQ4HcqyJISAc86dLgnOuYaU5zlpmhJj5FQrioIsy3DOudMlwTnXsEIIFEXBqRRjJMZICAHnnDtdEpxzDSvLMsqy5FSKMZKmKc45dzolOOcaVp7nVCoVTqWyLAkh4Jxzp1OCc66hpWlKjJFTpSgKsizDOedOpwTnXEMLIVAUBadCjJEYIyEEnHPudEpwzjW0LMsoy5JTIcZImqY459zpluCca2h5nlOpVDgVyrIkhIBzzp1uCc65hpemKZVKhZOtKAqyLMM55063BOdcw8vznLIsOZlijMQYCSHgnHOnW4JzruFNmTKFoig4mWKMpGmKc87VgwTnXMPL85wYIzFGTpayLAkh4Jxz9SBhjHv11Vf5zW9+g3ONLs9ziqLgZCmKgizLcM65epAwhu3evZtbbrmFV199FecaXZZllGXJyRBjJMZICAHnnKsHCWPU0NAQ3/zmN5k4cSLOjQZ5nlOpVDgZYoykaYpzztWLhDHqkUce4frrr+fzn/88zo0WaZoSY+REK8uSEALOOVcvEsagHTt28PLLL/ONb3yDo2VmmBlmxtq1a3GuHoUQKIqCE60oCrIswzlX/9auXYuZYWaYGaNVwhhz8OBBli9fziOPPMKxkIQkJHHXXXfhXD3KsoyyLDmRYozEGAkh4Jyrf3fddReSkIQkRquEUe6BBx5g+vTpTJ8+nS9+8YusXLmSyy67jL1791KWJYODg+zatQtJONfo8jynUqlwIsUYSdMU55yrJwmj3KJFi+jo6KCjo4OHH36Y888/n3feeYf169ezfv16fvnLX1KWJb29vTg3GqRpSlEUnChlWRJCwDnn6knCKPfZz36WGTNmMGPGDK666iruvvtuOjs76ezspLOzk8svv5wlS5bQ0tKCc6NBCIHu7m5OlKIoyLIM55yrJwnOuVElyzJijJwIMUZijIQQcM65epIwxnV2djJr1iycGy3yPKcoCmKMHK8YI2ma4pxz9SbBOTfq5HlOURQcr7IsCSHgnHP1JsE5N+pkWUZZlhyvoijIsgznnKs3Cc65USeEQKVS4XjEGIkxEkLAOefqTYJzbtRJ05Q0TYkx8knFGEnTFOecq0cJzrlRKYRAURR8UmVZEkLAOefqUYJzblTKsoyyLPmkiqIgyzKcc64eJTjnRqU8z6lUKnwSMUZijIQQcM65epTgnBu10jSlKAqOVYyRNE1xzrl6leCcG7VCCHR3d3OsyrIkhIBzztWrBOfcqJVlGZ9EURRkWYZzztWrBOfcqJXnOZVKhRgjRyvGSIyREALOOVevEpxzo1oIgaIoOFoxRtI0xTnn6lmCc25Ua25upixLjlZZloQQcM65epbgnBv1KpUKR6soCrIswznn6lmCc25Uy/OcNE2JMfJxYozEGAkh4Jxz9SzBOTfqhRAoioKPE2MkTVOcc67eJTjnRr0syyjLko9TliUhBJxzrt4lOOdGvTzPqVQqfJyiKMiyDOecq3cJzrkxIU1TYox8mBgjMUZCCDjnXL1LcM6NCSEEuru7+TAxRtI0xTnnGkGCc25MyLKMGCMfpixLQgg451wjSHDOjQl5nlOpVPgwRVGQZRnOOdcIEpxzY0YIgUqlwvvFGIkxEkLAOecaQYJzbswIIVCWJe8XYyRNU5xzrlEkOOfGjClTplCpVHi/siwJIeCcc40iCp7jQwAAGiZJREFUwTk3ZuR5TpqmxBgZqSgKsizDOecaRYJzbkwJIVAUBUfEGIkxEkLAOecaRYJzbkzJsoyyLDkixkiapjjnXCNJcM6NKXmeU6lUOKIsS0IIOOdcI0lwzo05aZoSY2RYURRkWYZzzjWSBOfcmBNCoCgKYozEGAkh4JxzjSTBOTfmZFlGWZbEGEnTFOecazQJzrkxJ89zKpUKZVkSQsA55xpNgnNuTErTlPvuu48sy3DOuUaT4Jwbk0IIpGlKCAHnnGs0Cc65Mam5uZm2tjacc64RJTjnxqQQAnme45xzjSjBOeecc67BJDjnnHPONZgE55xzzrkGk+Ccc84512ASnHPOOecaTIJzzjnnXINJcM4555xrMAnOOeeccw0mwTnnnHOuwSSMUYODg7z44ou8/PLLjHZr166lka1du5ZGtnbtWhrZ2rVraXRr166lka1du5ZGtnbtWhrZ2rVrcfUnYQwqy5Ibb7yRn/3sZ6xevZpbbrmFw4cPM1p1dHTQyDo6OmhkHR0dNLKOjg4aXUdHB42so6ODRtbR0UEj6+jowNWfhDGmVqtxzz338Oijj7J69Wr+5V/+hQMHDvD888/jnHPOucaQMMaUZcmf/umfctVVV3HEv//7vzNnzhycc8451xgSxpi3336bz3zmMyxfvpxp06Yxffp0nnrqKT7KVVddhZlhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBnDzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzBhmZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmxjAzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzhpkZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZgwzM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzY5iZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZw8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wwM8wMM8PMMDPMDDPDzDAzzAwzw8wYZmaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZoaZYWaYGWaGmWFmmBlmhplhZpgZZsZVV13FaJQwxuzevZvNmzfz53/+57z66qts2LCBf/7nf2bbtm18mB//+MdIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUlIQhKSkIQkJCEJSUhCEpKQhCQkIQlJSEISkpCEJCQhCUn8+Mc/ZjRKGOUeeOABpk+fzvTp0/niF7/In/3ZnzFlyhS+/OUvM8zMmD17Nj/72c9wzjnnXGNIGOUWLVpER0cHHR0dPPzww/zJn/wJ75ckCUmS4JxzzrnGkDDKffazn2XGjBnMmDGDq666iqamJgYHB+np6WHY4OAgW7duZd68eTjnnHOuMSSMMWeeeSYdHR1873vfY8GCBdxwww18+ctf5uqrr8Y555xzjSFhDPqrv/orenp66Orq4uWXX2bJkiU455xzrnEkjGHnnHMO48aNwznnnHONJcE555xzrsEkuA+1b98+tmzZgiQazeDgINVqlWq1SrVapVqtcvDgQRrF1q1beb99+/axZcsWJNEItm7dykiDg4NUq1Wq1SrVapVqtcrBgwepN7t372bLli309fXxfvv27WPLli1Iol7t3r2bLVu20NfXx0iDg4NUq1Wq1SrVapVqtcrBgwepN5LYsmULMUbeb9++fWzZsgVJ1DNJbNmyhRgjIw0ODlKtVqlWq1SrVarVKgcPHqRevfrqq/zmN79hpH379rFlyxYkUe9effVVfvOb33DE4OAg1WqVarVKtVqlWq1y8OBBGlWC+0CbNm1iwYIFbN68mcWLF7NmzRoaybPPPktzczOtra20trbS2trKa6+9RiN4/PHHWbp0KSNt2rSJBQsWsHnzZhYvXsyaNWuoZ48//jhLly5lpGeffZbm5mZaW1tpbW2ltbWV1157jXrywAMP0NrayubNm2lvb2fRokW8++67DNu0aRMLFixg8+bNLF68mDVr1lBvHnjgAVpbW9m8eTPt7e0sWrSId999l2HPPvsszc3NtLa20traSmtrK6+99hr15Ac/+AF33XUXL7zwArfffjudnZ0csWnTJhYsWMDmzZtZvHgxa9asoR794Ac/4K677uKFF17g9ttvp7OzkyOeffZZmpubaW1tpbW1ldbWVl577TXq0e7du7nlllt49dVXOWLTpk0sWLCAzZs3s3jxYtasWUO92r17N7fccguvvvoqRzz77LM0NzfT2tpKa2srra2tvPbaazSqBPcHarUabW1tdHd3s2rVKjZu3EhXVxcxRhrFrl27WLZsGX19ffT19dHX18fMmTOpZwcOHOCee+7hhz/8ISPVajXa2tro7u5m1apVbNy4ka6uLmKM1JsDBw5wzz338MMf/pD327VrF8uWLaOvr4++vj76+vqYOXMm9eIXv/gFTz/9NM888wyrVq3iueee43e/+x2bNm2iVqvR1tZGd3c3q1atYuPGjXR1dRFjpF784he/4Omnn+aZZ55h1apVPPfcc/zud79j06ZNDNu1axfLli2jr6+Pvr4++vr6mDlzJvViYGCAH/3oR/z0pz/loYceYsOGDaxZs4bBwUFqtRptbW10d3ezatUqNm7cSFdXFzFG6snAwAA/+tGP+OlPf8pDDz3Ehg0bWLNmDYODgwzbtWsXy5Yto6+vj76+Pvr6+pg5cyb1ZmhoiG9+85tMnDiRI2q1Gm1tbXR3d7Nq1So2btxIV1cXMUbqzdDQEN/85jeZOHEiI+3atYtly5bR19dHX18ffX19zJw5k0aV4P7ASy+9xIQJE7jkkksYdt5553Httdeybds2GkV/fz8XX3wxg4ODDA0N0QgeffRRzjvvPB588EFGeumll5gwYQKXXHIJw8477zyuvfZatm3bRr159NFHOe+883jwwQd5v/7+fi6++GIGBwcZGhqi3kyYMIHOzk4mTJjAERdddBG/+tWveOmll5gwYQKXXHIJw8477zyuvfZatm3bRr2YMGECnZ2dTJgwgSMuuugifvWrXzGsv7+fiy++mMHBQYaGhqg3F198Mc8++ywTJkxg2JlnnkmtVmNoaIiXXnqJCRMmcMkllzDsvPPO49prr2Xbtm3Uk4svvphnn32WCRMmMOzMM8+kVqsxNDTEsP7+fi6++GIGBwcZGhqiXj3yyCNcf/31fP7zn+eIl156iQkTJnDJJZcw7LzzzuPaa69l27Zt1JtHHnmE66+/ns9//vOM1N/fz8UXX8zg4CBDQ0M0ugT3Bw4cOMCll17KSOeeey5vvPEGjaBWq/HWW29x//33c+ONNzJt2jTuvfde6t3y5cv59re/zTnnnMNIBw4c4NJLL2Wkc889lzfeeIN6s3z5cr797W9zzjnnMFKtVuOtt97i/vvv58Ybb2TatGnce++91JMLL7yQGTNmcMTevXvp6elh9uzZHDhwgEsvvZSRzj33XN544w3qxYUXXsiMGTM4Yu/evfT09DB79mxqtRpvvfUW999/PzfeeCPTpk3j3nvvpZ4kScIll1xCrVbj6aefprm5ma9//etMmjSJAwcOcOmllzLSueeeyxtvvEE9SZKESy65hFqtxtNPP01zczNf//rXmTRpErVajbfeeov777+fG2+8kWnTpnHvvfdSb3bs2MHLL7/MN77xDUY6cOAAl156KSOde+65vPHGG9STHTt28PLLL/ONb3yDkWq1Gm+99Rb3338/N954I9OmTePee++lkSW4P1Cr1UiShJGSJOHw4cM0gv379zNr1iyeeOIJent76enpYevWrWzYsIF6liQJH6RWq5EkCSMlScLhw4epN0mS8EH279/PrFmzeOKJJ+jt7aWnp4etW7eyYcMG6tH+/fvJ85wlS5bwhS98gVqtRpIkjJQkCYcPH6Ye7d+/nzzPWbJkCV/4whfYv38/s2bN4oknnqC3t5eenh62bt3Khg0bqDeDg4O8++67XHDBBWzfvp0DBw5Qq9VIkoSRkiTh8OHD1KPBwUHeffddLrjgArZv386BAwfYv38/s2bN4oknnqC3t5eenh62bt3Khg0bqBcHDx5k+fLlPPLII7xfrVYjSRJGSpKEw4cPUy8OHjzI8uXLeeSRR3i//fv3M2vWLJ544gl6e3vp6elh69atbNiwgUaV4P7AWWedRa1WY6TDhw9zxhln0AgmT57MY489xuTJkxk2adIkZs+ezSuvvEIjOuuss6jVaox0+PBhzjjjDBrF5MmTeeyxx5g8eTLDJk2axOzZs3nllVeoNzt37uSmm27i1ltvZfHixQw766yzqNVqjHT48GHOOOMM6s3OnTu56aabuPXWW1m8eDHDJk+ezGOPPcbkyZMZNmnSJGbPns0rr7xCvTn//PO59dZbefLJJzn77LPp7u7mrLPOolarMdLhw4c544wzqEfnn38+t956K08++SRnn3023d3dTJ48mccee4zJkyczbNKkScyePZtXXnmFerFy5Uouu+wy9u7dS1mWDA4OsmvXLiRx1llnUavVGOnw4cOcccYZ1IuVK1dy2WWXsXfvXsqyZHBwkF27diGJyZMn89hjjzF58mSGTZo0idmzZ/PKK6/QqBLcH7jgggt4/fXXGentt9/miiuuoBHs3buXjRs3MtJ7773HuHHjaEQXXHABr7/+OiO9/fbbXHHFFTSKvXv3snHjRkZ67733GDduHPWkt7eX2267jfvuu4+WlhaOuOCCC3j99dcZ6e233+aKK66gnvT29nLbbbdx33330dLSwhF79+5l48aNjPTee+8xbtw46sWbb77JT37yE0b69Kc/za9//WsuuOACXn/9dUZ6++23ueKKK6gnb775Jj/5yU8Y6dOf/jS//vWv2bt3Lxs3bmSk9957j3HjxlEvzj//fN555x3Wr1/P+vXr+eUvf0lZlvT29nLBBRfw+uuvM9Lbb7/NFVdcQb04//zzeeedd1i/fj3r16/nl7/8JWVZ0tvby969e9m4cSMjvffee4wbN45GleD+wJVXXsmwsiwZNjAwQG9vL9dccw2N4NChQ7S1tbF7926G7d+/nxdeeIF58+bRiK688kqGlWXJsIGBAXp7e7nmmmtoFIcOHaKtrY3du3czbP/+/bzwwgvMmzePerFv3z7uvPNOVq5cSVNTE0NDQwwNDVGr1bjyyisZVpYlwwYGBujt7eWaa66hXuzbt48777yTlStX0tTUxNDQEENDQ9RqNQ4dOkRbWxu7d+9m2P79+3nhhReYN28e9aJWq/HQQw/x5ptvMuy3v/0t27ZtY/bs2Vx55ZUMK8uSYQMDA/T29nLNNddQT2q1Gg899BBvvvkmw37729+ybds2Zs+ezaFDh2hra2P37t0M279/Py+88ALz5s2jXtx99910dnbS2dlJZ2cnl19+OUuWLKGlpYUrr7ySYWVZMmxgYIDe3l6uueYa6sXdd99NZ2cnnZ2ddHZ2cvnll7NkyRJaWlo4dOgQbW1t7N69m2H79+/nhRdeYN68eTSqBPcHkiRh9erVLF26lObmZhYuXMiKFSuYOHEijcDMWLZsGTfffDPNzc3MnTuX22+/nZkzZ9KIkiRh9erVLF26lObmZhYuXMiKFSuYOHEijcLMWLZsGTfffDPNzc3MnTuX22+/nZkzZ1Iv1q9fzzvvvMMdd9zB1KlTmTp1KlOnTuX73/8+SZKwevVqli5dSnNzMwsXLmTFihVMnDiRerF+/Xreeecd7rjjDqZOncrUqVOZOnUq3//+9zEzli1bxs0330xzczNz587l9ttvZ+bMmdSLz33uc9x777186Utf4qtf/SqzZs3i1ltv5brrriNJElavXs3SpUtpbm5m4cKFrFixgokTJ1JPPve5z3HvvffypS99ia9+9avMmjWLW2+9leuuuw4zY9myZdx88800Nzczd+5cbr/9dmbOnEkjSJKE1atXs3TpUpqbm1m4cCErVqxg4sSJNAIzY9myZdx88800Nzczd+5cbr/9dmbOnEmjSnAf6Oqrr2b79u2sW7eOHTt2MGfOHBrJokWLqFarrFu3jmq1SktLC40iyzK2bt3KSFdffTXbt29n3bp17Nixgzlz5lDPsixj69atjLRo0SKq1Srr1q2jWq3S0tJCPfnud7+LJCQhCUlIYvny5Qy7+uqr2b59O+vWrWPHjh3MmTOHevLd734XSUhCEpKQxPLlyxm2aNEiqtUq69ato1qt0tLSQr1ZuHAhfX19rFixgldeeYWvfe1rHHH11Vezfft21q1bx44dO5gzZw71aOHChfT19bFixQpeeeUVvva1r3HEokWLqFarrFu3jmq1SktLC/Wss7OTWbNmccTVV1/N9u3bWbduHTt27GDOnDnUs87OTmbNmsURixYtolqtsm7dOqrVKi0tLTSyBPeRxo8fT5IkNKIkSRg/fjxJkjBajB8/niRJaFRJkjB+/HiSJKFRjR8/niRJaERJkjB+/HiSJKFeJUnCxIkTGTduHB9k/PjxJElCPUuShIkTJzJu3DjeL0kSxo8fT5IkNKrx48eTJAmNKEkSxo8fT5IkNLoE55xzzrkGk+Ccc84512ASnHPOOecaTIJzzjnnXINJcM4555xrMAnOOeeccw0mwTnnTrG7776bO++8k29961t8Ek8++SR33nknd955J7t378Y5N/YkOOfcKfbiiy/ymc98hhtuuIFP4oorruD666/nP//zPzlw4ADOubEnwTnnToMrrriCG264gU9i+vTp/PVf/zXOubErwTnnToAXX3yRv//7v2f69OnccMMNPP744xytr33ta2zZsoV/+Id/YPr06SxYsIC9e/fyr//6r8yaNYsrr7yShx56COecOyLBOeeO0969e1m8eDFf/vKXeemll7jnnnv44Q9/yMaNGzkaRVHQ1tbGzTffTEdHB++++y633HILzz//PN/73vdYunQpP/7xj3n++edxzrlhCc45d5xijIwbN44ZM2Zw7rnn0tTUxFNPPcXll1/O0brtttuYN28eM2bM4JZbbuG///u/efjhh5kxYwY33XQTf/EXf8HPf/5znHNuWIJzzh2nmTNncumllzJr1iwWLFjAY489xtlnn42ZcbSmTJnCEeeccw5nn3025557Lkf88R//MbVaDeecG5bgnHPHady4cWzcuJF/+qd/4rOf/Sz/9m//xt/+7d/y1FNP4ZxzJ0OCc84dpzfffJP/+I//oKmpiQcffJAXX3yRPM958skncc65kyHBOeeO029/+1u+9a1v8fLLLzPs8OHDvPXWW3zuc5/DOedOhgTnnDtOV111FUuWLOG2227jL//yL5k2bRr79+9n5cqVOOfcyZDgnHMnwJ133snOnTvZvHkz1WqVZ555hgsvvJCjIYlZs2ZxxN/8zd/Q19fHSJ2dnSxfvhznnBuW4JxzJ0iSJEyaNImzzjqLj3Po0CF+//vf80m8++67vPPOOzjnxq4E55w7xc4880zuuecerrvuOj6J73znO1x33XX8n//zf0iSBOfc2JPgnHOnWF9fHzt37uS//uu/+CTWrFnDzp072blzJ9OnT8c5N/YkOOecc841mATnnHPOuQaT4JxzzjnXYBKcc8455xpMgnPOOedcg/m/DRbCWcnMAIAAAAAASUVORK5CYII=" style="width: 560px;"></div></div></div></div></div><div class = 'S0'><span>We observe an oscillation of the particle's position along the beam line. This type of transverse oscillation is usually called </span><span style=' font-style: italic;'>betatron oscillation</span><span>.</span></div><div class = 'S0'><span>Now you can go back to the start of the script and change the focal length F and observe how the oscillation changes. You might also want to increase the length of the beam line by adding more copies of the FODO cells in the definition of beamline by increasing the number of copies in </span><span style=' font-family: monospace;'>repmat()</span><span>. Explore!</span></div><div class = 'S0'><span>Calculating transfer matrices is a task we repeatedly have to do and it is therefore convenient to encapsulate the code that does that in a dedicated function; we call it </span><span style=' font-family: monospace;'>calcmat().</span><span> It receives the </span><span style=' font-family: monospace;'>beamline</span><span> as input and returns </span><span style=' font-family: monospace;'>Racc, spos,nmat</span><span>, and</span><span style=' font-family: monospace;'> nlines</span><span>. Internally it calculates nmat and nlines and then allocates arrays for </span><span style=' font-family: monospace;'>Racc </span><span>and </span><span style=' font-family: monospace;'>spos</span><span> and then loops over all lines and segments, just as we did in the code above. We use such a function in the following way</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S12'><span style="white-space: pre"><span >[Racc,spos,nmat,nlines]=calcmat(beamline);</span></span></div></div></div><div class = 'S0'><span>and--viola--all transfer matrces and the positions of all segments are immediately available. Here is the definition of the function </span><span style=' font-family: monospace;'>calcmat()</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">function </span><span >[Racc,spos,nmat,nlines]=calcmat(beamline)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >nlines=size(beamline,1);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >nmat=sum(beamline(:,2))+1;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >Racc=zeros(2,2,nmat);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >Racc(:,:,1)=eye(2);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >spos=zeros(nmat,1);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >ic=1;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >line=1:nlines</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">for </span><span >seg=1:beamline(line,2)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > ic=ic+1;</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Rcurr=eye(2);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">switch </span><span >beamline(line,1)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">case </span><span >1</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Rcurr=D2(beamline(line,3));</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">case </span><span >2</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Rcurr=Q2(beamline(line,4));</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">otherwise</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > disp(</span><span style="color: rgb(170, 4, 249);">'unknown element'</span><span >);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > Racc(:,:,ic)=Rcurr*Racc(:,:,ic-1);</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > spos(ic)=spos(ic-1)+beamline(line,3); </span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span > </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div class = 'S5'><span>Inside this function the function D(L) an Q(F) are unknown. We could either copy the two lines with their definitions into calcmat, or define additional functions that provide the same functionality. First for the drift space</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">function </span><span >out=D2(L)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >out=[1,L;0,1];</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div class = 'S5'><span>and for the lens</span></div><div class="CodeBlock"><div class="inlineWrapper"><div class = 'S2'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">function </span><span >out=Q2(F)</span></span></div></div><div class="inlineWrapper"><div class = 'S3'><span style="white-space: pre"><span >out=[1,0;-1/F,1];</span></span></div></div><div class="inlineWrapper"><div class = 'S4'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div class = 'S5'><span>Here I have given new names (D2, Q2)to functions for drift and quad, because the names "D" and "Q" are already used for the matrices in the first few lines of this script.</span></div><div class = 'S0'><span>These service functions will make all further examples much more compact and readable, because we never have to copy and paste the enclosed code verbatim into our scripts. This simple tutorial only contains the bare essentials of beam optics calculations. Have a look at the web site where you picked this script up and check out the other more advanced examples. </span></div>
<br>
<!--
##### SOURCE BEGIN #####
%%
% Companion software for "Volker Ziemann, _Hands-on Accelerator physics using
% MATLAB, CRCPress, 2019_" (https://www.crcpress.com/9781138589940)
%% Beam optics tutorial (Section 3.3.1)
% Volker Ziemann, 211106, CC-BY-SA-4.0
%
% In this tutorial we'll develop a first beam optics code that allows us to
% track single particles, or "rays" through a sequence of drift spaces and thin
% lenses. Let's therefore first clear the MATLAB workspace and define anonymous
% functions (those with an '@') that return the transfer matrices. Here D(L) will
% return the 2x2 matrix for a drift space and Q(F) that of a thin lens with focal
% length F
clear all; close all
D=@(L)[1,L;0,1];
Q=@(F)[1,0;-1/F,1];
%%
% Now it is easy to obtain the transfer matrix of a beam line consisting of
% a sequnece of such elements. For example, a beamline, where the first element
% is a drift space with a length of 2m, followed by a lens with a focal length
% of 1 m and a second drift space having a length of 3m.
A=D(3)*Q(1)*D(2)
%%
% I suggest to verify the numerical calculation by multiplying the individual
% matrices by hand on a peice of paper.
%
% Directly coding the matrix multiplications is inconvenient and we therefore
% introduce a very simple way to describe sequences of elements as an array, here
% called |fodo.| This array has four columns with the following interpretations
%%
% * First column: code, drift=1, quad=2
% * Second column: repeat code for the element
% * Third column: length of one segment
% * Fourth column: strength of one segment, focal length for a thin quadrupole.
%%
% The second column with the repeat code is not strictly necessary but will
% make the graphs easier to interpret, because we can plot many intermediate points,
% rather than connecting values at the end of elements by a straight line.
%
% This particular sequence starts with an element having a 1 in the first column
% and is therefore a drift space. It consists of 5 segments, each 0.2 m long.
% The last column is not used here and contains a zero, The next element is a
% lens, beacuse the first code is 2 and it is defocusing, because the fourth column
% contains a negative number for the focal length. The third element is again
% a driftspace, this time with 10 segments of 0.2 m, followed by a focusing lens
% and another drift space. Such a sequence of focusing lens, drift space, defocusing
% lens and another drift space is usually referred to as a FODO (get it?) cell.
F=2.1;
fodo=[1, 5, 0.2, 0;
2, 1, 0, -F;
1, 10, 0.2, 0;
2, 1, 0, F;
1, 5, 0.2, 0]
%%
% Now we build our beamline as a sequence of one or many FODO cells. The latter
% we achieve by simply concatenating |fodo| vertically (note the semicolon!).
% We will always use |beamline| as a reserved name of the beamline that we work
% on.
%beamline=fodo; % one cell
%beamline=[fodo;fodo] % two cells
beamline=repmat(fodo,10,1) % ten cells
%%
% Before we can assemble the transfer matrices, we have to do a bit of accounting,
% like counting the number of elements |nlines| in the beamline description and
% the number |nmat| of matrices that we will need to store the matrices. This
% is just given by the sum of the segments in the second column plus one. Next
% we allocate space for the nmat $2\times 2$ matrices in the array Racc and initialize
% the first matrix as the $2\times2$unit matrix eye(2). This additional matrix
% at the beginning accounts for the "plus one" in the definition of |nmat|. Finally
% we allocate the array |spos| to store the position of each segment along the
% beam line.
nlines=size(beamline,1)
nmat=sum(beamline(:,2))+1;
Racc=zeros(2,2,nmat);
Racc(:,:,1)=eye(2);
spos=zeros(nmat,1);
%%
% Now we are ready to loop through the beamline description by looping over
% the lines, labeled by |line|, and then over the segments, labeled by |seg|.
% As we do that we increment the counter |ic| that keeps track of which segment
% we're currently working on. Then we initialize the transfer matrix of the current
% segment |Rcurr| by the unit marix before looking at the first column to decide
% whether the element is a drift space or a lens and fill |Rcurr| with the appropriate
% transfer matrix. In |case| the first column (|beamline(line,1)|) is a 1, it
% is a drift and we use the value from the third column as the length of the drift
% space. In |case| the first column contains a 2, it is a lens and we use the
% focal length from the fourth column as the argument to the function |Q(F).|
% In case we encounter any other number in the first column, we display a warning.
% At this point |Rcurr| contains the transfer matrix for the current segment and
% we find *the transfer matrix from the start of the beamline to the end of the
% current element* |Racc(:,:,ic)| by left-multiplying Rcurr to |Racc(:,:,ic-1)|
% which already contains the transfer matrix from the start to the end of he previous
% segment. Finally we calculate the position |spos(ic)| of the current segment
% from that of the previous segment |spos(ic-1)| by adding the length of the current
% segment to it.
ic=1;
for line=1:nlines
for seg=1:beamline(line,2)
ic=ic+1;
Rcurr=eye(2);
switch beamline(line,1)
case 1
Rcurr=D(beamline(line,3));
case 2
Rcurr=Q(beamline(line,4));
otherwise
disp('unknown element');
end
Racc(:,:,ic)=Rcurr*Racc(:,:,ic-1);
spos(ic)=spos(ic-1)+beamline(line,3);
end
end
%%
% At this point |Racc| contains all trasnfer matrices from the start of the
% beam line to the end of each and every segment while |spos| contains the positions
% of each segment.
%%
% Now we can use the transfer matrices to plot the trajectory of a particle
% with starting coordinates $x_0 = 0$ and $x_0'=10^{-3}$ rad by mapping the initial
% (column) vector |x0| with |Racc| to the ened of each segment and record both
% position |x(1)| and angle |x(2)| in the previously allocated array |data|.
data=zeros(nmat,2);
x0=[0;1e-3];
for k=1:nmat
x=Racc(:,:,k)*x0;
data(k,1)=x(1);
data(k,2)=x(2);
end
%%
% Now we can plot the trajectory and annotate the axes.
plot(spos,1000*data(:,1),'k')
xlabel('s [m]')
ylabel('x [mm]')
title('Trajectory')
%%
% We observe an oscillation of the particle's position along the beam line.
% This type of transverse oscillation is usually called _betatron oscillation_.
%
% Now you can go back to the start of the script and change the focal length
% F and observe how the oscillation changes. You might also want to increase the
% length of the beam line by adding more copies of the FODO cells in the definition
% of beamline by increasing the number of copies in |repmat()|. Explore!
%
% Calculating transfer matrices is a task we repeatedly have to do and it is
% therefore convenient to encapsulate the code that does that in a dedicated function;
% we call it |calcmat().| It receives the |beamline| as input and returns |Racc,
% spos,nmat|, and |nlines|. Internally it calculates nmat and nlines and then
% allocates arrays for |Racc| and |spos| and then loops over all lines and segments,
% just as we did in the code above. We use such a function in the following way
[Racc,spos,nmat,nlines]=calcmat(beamline);
%%
% andREPLACE_WITH_DASH_DASHviolaREPLACE_WITH_DASH_DASHall transfer matrces and the positions of all segments are immediately
% available. Here is the definition of the function |calcmat()|
function [Racc,spos,nmat,nlines]=calcmat(beamline)
nlines=size(beamline,1);
nmat=sum(beamline(:,2))+1;
Racc=zeros(2,2,nmat);
Racc(:,:,1)=eye(2);
spos=zeros(nmat,1);
ic=1;
for line=1:nlines
for seg=1:beamline(line,2)
ic=ic+1;
Rcurr=eye(2);
switch beamline(line,1)
case 1
Rcurr=D2(beamline(line,3));
case 2
Rcurr=Q2(beamline(line,4));
otherwise
disp('unknown element');
end
Racc(:,:,ic)=Rcurr*Racc(:,:,ic-1);
spos(ic)=spos(ic-1)+beamline(line,3);
end
end
end
%%
% Inside this function the function D(L) an Q(F) are unknown. We could either
% copy the two lines with their definitions into calcmat, or define additional
% functions that provide the same functionality. First for the drift space
function out=D2(L)
out=[1,L;0,1];
end
%%
% and for the lens
function out=Q2(F)
out=[1,0;-1/F,1];
end
%%
% Here I have given new names (D2, Q2)to functions for drift and quad, because
% the names "D" and "Q" are already used for the matrices in the first few lines
% of this script.
%
% These service functions will make all further examples much more compact and
% readable, because we never have to copy and paste the enclosed code verbatim
% into our scripts. This simple tutorial only contains the bare essentials of
% beam optics calculations. Have a look at the web site where you picked this
% script up and check out the other more advanced examples.
##### SOURCE END #####
-->
</div></body></html>