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
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
|
import calculateField as magnet
import fieldManipulation as field
import numpy as np
import readData as read
import matplotlib.pyplot as plot
from matplotlib.backends.backend_pdf import PdfPages
def createTestField(outputFile):
"""
Create a field map to test other field manipulation routines on. Save it so that
it can be read back in later to speed up testing.
Inputs:
- outputFile: File name to save an approx. Spectrometer Solenoid map as
Returns:
- 2D meshed grid of (r, z) and (Br, Bz)
"""
layers = ([42, 28, 56, 20, 62])
turns = ([115, 114, 64, 768, 64])
currents = ([264.8, 285.6, 233.7, 275.7, 240.2])
radii = ([258.0e-3, 258.0e-3, 258.0e-3, 258.0e-3, 258.0e-3])
length = ([201.3e-3, 199.5e-3,110.6e-3, 1314.3e-3, 110.6e-3])
thickness = ([46.2e-3, 30.9e-3, 60.9e-3, 22.1e-3, 67.8e-3])
dLayer1 = thickness[0]/(1.0*layers[0] - 1.0)
dLayer2 = thickness[1]/(1.0*layers[1] - 1.0)
dLayer3 = thickness[2]/(1.0*layers[2] - 1.0)
dLayer4 = thickness[3]/(1.0*layers[3] - 1.0)
dLayer5 = thickness[4]/(1.0*layers[4] - 1.0)
dTurn1 = length[0]/(1.0*turns[0] - 1.0)
dTurn2 = length[1]/(1.0*turns[1] - 1.0)
dTurn3 = length[2]/(1.0*turns[2] - 1.0)
dTurn4 = length[3]/(1.0*turns[3] - 1.0)
dTurn5 = length[4]/(1.0*turns[4] - 1.0)
layerSeparation = ([dLayer1, dLayer2, dLayer3, dLayer4, dLayer5])
turnSeparation = ([dTurn1, dTurn2, dTurn3, dTurn4, dTurn5])
upstreamPositions = ([0.0, 440.9e-3, 885.3e-3, 1033.5e-3, 2385.3e-3])
r = ([0.0, 30.0e-3, 60.0e-3, 90.0e-3, 120.0e-3, 150.0e-3, 180.0e-3])
z = np.arange(-0.5, 3.0, 15.0e-3)
Z, R = np.meshgrid(z, r)
Br, Bz = magnet.makeMagnet(5, layers, turns, layerSeparation, turnSeparation,
upstreamPositions, currents, radii, R, Z)
description = '# Approximate field of one Spectrometer Solenoid using parameters from MICE Note 397:\n'
description += '#\t\tcurrent = '+str(currents)+' A\n'
description += '#\t\tinnerRadius = '+str(radii)+' m\n'
description += '#\t\t'+str(turns)+' turns of conductor per layer\n'
description += '#\t\t'+str(layers)+' layers of conductor per coil\n'
description += '#\t\tupstream locations = '+str(upstreamPositions)+' m'
magnet.printField(R, Z, Br, Bz, outputFile, description)
return R, Z, Br, Bz
def createInterpolationField(outputFile):
"""
Create a fine meshed field map to test interpolation routines on. Save it so that
it can be read back in later to speed up testing.
Inputs:
- outputFile: File name to save an approx. Spectrometer Solenoid map as
Returns:
- Nothing
"""
layers = ([42, 28, 56, 20, 62])
turns = ([115, 114, 64, 768, 64])
currents = ([264.8, 285.6, 233.7, 275.7, 240.2])
radii = ([258.0e-3, 258.0e-3, 258.0e-3, 258.0e-3, 258.0e-3])
length = ([201.3e-3, 199.5e-3,110.6e-3, 1314.3e-3, 110.6e-3])
thickness = ([46.2e-3, 30.9e-3, 60.9e-3, 22.1e-3, 67.8e-3])
dLayer1 = thickness[0]/(1.0*layers[0] - 1.0)
dLayer2 = thickness[1]/(1.0*layers[1] - 1.0)
dLayer3 = thickness[2]/(1.0*layers[2] - 1.0)
dLayer4 = thickness[3]/(1.0*layers[3] - 1.0)
dLayer5 = thickness[4]/(1.0*layers[4] - 1.0)
dTurn1 = length[0]/(1.0*turns[0] - 1.0)
dTurn2 = length[1]/(1.0*turns[1] - 1.0)
dTurn3 = length[2]/(1.0*turns[2] - 1.0)
dTurn4 = length[3]/(1.0*turns[3] - 1.0)
dTurn5 = length[4]/(1.0*turns[4] - 1.0)
layerSeparation = ([dLayer1, dLayer2, dLayer3, dLayer4, dLayer5])
turnSeparation = ([dTurn1, dTurn2, dTurn3, dTurn4, dTurn5])
upstreamPositions = ([0.0, 440.9e-3, 885.3e-3, 1033.5e-3, 2385.3e-3])
r = np.arange(0.0, 200.0e-3, 10.0e-3)
z = np.arange(-0.5, 3.0, 5.0e-3)
Z, R = np.meshgrid(z, r)
Br, Bz = magnet.makeMagnet(5, layers, turns, layerSeparation, turnSeparation,
upstreamPositions, currents, radii, R, Z)
description = '# Approximate field of one Spectrometer Solenoid using parameters from MICE Note 397:\n'
description += '#\t\tcurrent = '+str(currents)+' A\n'
description += '#\t\tinnerRadius = '+str(radii)+' m\n'
description += '#\t\t'+str(turns)+' turns of conductor per layer\n'
description += '#\t\t'+str(layers)+' layers of conductor per coil\n'
description += '#\t\tupstream locations = '+str(upstreamPositions)+' m'
magnet.printField(R, Z, Br, Bz, outputFile, description)
def test_convertToList(r, z, Br, Bz, inputFile):
"""
Convert 2d gridded data into a polarMeasurements list. This should compare exactly to
what we read in from Tests/test_fieldManipulation1.dat
Inputs:
- r: 2D gridded radial co-ordinates (m)
- z: 2D gridded longitudinal co-ordinates (m)
- Br: Radial field component calculated on the (r, z) grid
- Bz: Longitudinal field component calculated on the (r, z) grid
- inputFile: String denoting the file we'll compare against
Returns:
- Boolean value based on whether the converted gridded data is equivalent to the
data read from Tests/test_fieldManipulation1.dat
"""
listedData = field.convertToList(r, z, Br, Bz)
readData = read.readFile(inputFile)
# Both convertToList() and readData.readFile() should return data that's sorted according
# to increasing z.
if len(listedData) != len(readData):
print 'listedData length != readData length, halting test.'
return False
for i in range(0, len(listedData)):
if ((listedData[i].r - readData[i].r < 1.0e-6)
and (listedData[i].phi - readData[i].phi < 1.0e-6)
and (listedData[i].z - readData[i].z < 1.0e-6)
and (listedData[i].Br - readData[i].Br < 1.0e-6)
and (listedData[i].Bphi - readData[i].Bphi < 1.0e-6)
and (listedData[i].Bz - readData[i].Bz < 1.0e-6)):
return True
else:
if listedData[i].r != readData[i].r:
print 'listedData.r != readData.r ---> ', listedData[i].r, ' != ', readData[i].r
if listedData[i].phi != readData[i].phi:
print 'listedData.phi != readData.phi ---> ', listedData[i].phi, ' != ', readData[i].phi
if listedData[i].z != readData[i].z:
print 'listedData.z != readData.z ---> ', listedData[i].z, ' != ', readData[i].z
if listedData[i].Br != readData[i].Br:
print 'listedData.Br != readData.Br ---> ', listedData[i].Br, ' != ', readData[i].Br
if listedData[i].Bphi != readData[i].Bphi:
print 'listedData.Bphi != readData.Bphi ---> ', listedData[i].Bphi, ' != ', readData[i].Bphi
if listedData[i].Bz != readData[i].Bz:
print 'listedData.Bz != readData.Bz ---> ', listedData[i].Bz, ' != ', readData[i].Bz
return False
def test_convertToGrid(r, z, Br, Bz, inputFile):
"""
Test the conversion of a set of polarMeasurements to a 2D grid.
Inputs:
- r: 2D gridded radial co-ordinates to compare against (m)
- z: 2D gridded longitudinal co-ordinates to compare against (m)
- Br: Radial field components to compare against calculated on the (r, z) grid
- Bz: Longitudinal field components to compare against calculated on the (r, z) grid
- inputFile: Location of the file we want to test
Returns:
- Boolean statement of whether the gridded input file equals the original 2D gridded data
"""
grid_r, grid_z, grid_Br, grid_Bz = field.convertToGrid(inputFile)
rCheck = np.allclose(r, grid_r)
zCheck = np.allclose(z, grid_z)
BrCheck = np.allclose(Br, grid_Br)
BzCheck = np.allclose(Bz, grid_Bz)
if rCheck and zCheck and BrCheck and BzCheck:
return True
else:
print 'r vs. grid_r = ', rCheck
print 'z vs. grid_z = ', zCheck
print 'Br vs. grid_Br = ', BrCheck
print 'Bz vs. grid_Bz = ', BzCheck
return False
def test_interpolation(inputFile1, inputFile2):
"""
Check that interpolating the data on a broad grid gives us similar results to what
we calculate on a finer mesh.
Inputs:
- inputFile1: Coarse gridded data
- inputFile2: Finely gridded data, with points overlapping those of outputFile1
"""
pp = PdfPages('Tests/test_fieldManipulation.pdf')
coarseField = read.readFile(inputFile1)
fineField = read.readFile(inputFile2)
coarse_r, coarse_z, coarse_Br, coarse_Bz = field.convertToGrid(coarseField)
fine_r, fine_z, fine_Br, fine_Bz = field.convertToGrid(fineField)
# Now we check the interpolator
interpolated_r, interpolated_z, interpolated_Br, interpolated_Bz = field.interpolate(fine_r, fine_z, coarse_r, coarse_z, coarse_Br, coarse_Bz)
Br_diff = fine_Br - interpolated_Br
Bz_diff = fine_Bz - interpolated_Bz
plot.plot(fine_z, Bz_diff, 'k-')
plot.title('Interpolated field - calculated field')
plot.xlabel('z (m)')
plot.ylabel('dBz (T)')
plot.savefig(pp, format='pdf')
plot.show()
plot.plot(fine_z, Br_diff, 'k-')
plot.title('Interpolated field - calculated field')
plot.xlabel('z (m)')
plot.ylabel('dBr (T)')
plot.savefig(pp, format='pdf')
plot.show()
pp.close()
if np.all(Br_diff) < 1.0e-6 and np.all(Bz_diff) < 1.0e-6:
result = True
else:
result = False
return result
def test_scaleField(inputFile):
"""
Check that length and field scaling is working properly
"""
lengthScale = 2.0 # should double the length of the magnet
fieldScale = 2.0 # should double the field of the magnet
originalField = read.readFile(inputFile)
original_r, original_z, original_Br, original_Bz = field.convertToGrid(originalField)
scaledField = field.scaleField(originalField, lengthScale, fieldScale)
scaled_r, scaled_z, scaled_Br, scaled_Bz = field.convertToGrid(scaledField)
checkLength_z = np.multiply(2.0, original_z)
checkLength_r = np.multiply(2.0, original_r)
checkField_z = np.multiply(2.0, original_Bz)
checkField_r = np.multiply(2.0, original_Br)
if np.array_equal(scaled_r, checkLength_r) and np.array_equal(scaled_z, checkLength_z):
length = True
else:
length = False
if np.array_equal(scaled_Br, checkField_r) and np.array_equal(scaled_Bz, checkField_z):
scale = True
else:
scale = False
if length and scale:
return True
else:
return False
def test_addFields(inputFile):
"""
Check that we can add fields in proportion -- take the same field "twice" and
add them together in 50% proportion.
"""
originalField = read.readFile(inputFile)
newField = field.addFields(originalField, originalField, 0.5)
original_r, original_z, original_Br, original_Bz = field.convertToGrid(originalField)
new_r, new_z, new_Br, new_Bz = field.convertToGrid(newField)
if np.array_equal(original_r, new_r) and np.array_equal(original_z, new_z) and np.array_equal(original_Br, new_Br) and np.array_equal(original_Bz, new_Bz):
return True
else:
return False
outputFile1 = 'Tests/test_fieldManipulation1.dat'
outputFile2 = 'Tests/test_fieldManipulation2.dat'
r, z, Br, Bz = createTestField(outputFile1)
createInterpolationField(outputFile2)
result1 = test_convertToList(r, z, Br, Bz, outputFile1)
result2 = test_convertToGrid(r, z, Br, Bz, outputFile1)
result3 = test_interpolation(outputFile1, outputFile1)
result4 = test_scaleField(outputFile1)
result5 = test_addFields(outputFile1)
print 'test_convertToList() returns ', result1
print 'test_convertToGrid() returns ', result2
print 'test_interpolation() returus ', result3
print 'test_scaleField() returns ', result4
print 'test_addFields() returns ', result5
|