← Back to branch summary

~plane1/+junk/dissertation

« back to all changes in this revision

Viewing changes to src/Analysis/DLHPC/DLHPRF_rf_hammer.py

  • Committer: plane1
  • Date: 2016-06-23 20:36:15 UTC
  • Revision ID: plane1@hawk.iit.edu-20160623203615-dqlye6youiil2bbh
Changed directory names for HPRF cavities to the acronyms now used in the dissertation (i.e., HC and DLHC)

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Analysis/DLHPC/DLHPRF_rf_hammer.py
 
1
# -*- coding: utf-8 -*-
 
2
"""
 
3
Created on Sun Jan 25 20:32:38 2015
 
4
 
 
5
@author: plane
 
6
"""
 
7
import math
 
8
import numpy as np
 
9
from scipy import fftpack
 
10
import scipy.signal
 
11
import matplotlib.pyplot as plt  # Matplotlib's pyplot: MATLAB-like syntax
 
12
from bisect import bisect_left, bisect_right
 
13
from functools import partial
 
14
from operator import add
 
15
from pgl.superfish import parse_sf7
 
16
from pgl.curve import polynomials, fit_spliced_polynomials, gaussian, fit_gaussian
 
17
 
 
18
sf7_dir = 'C:\\Users\\plane\\Documents\\Superfish\\DL HPRF\\99.5% Donut\\'
 
19
sf7_file = 'DL_HPRF_100x100_field_map_20_MV_m.txt'
 
20
field_map = parse_sf7(''.join((sf7_dir, sf7_file)))
 
21
# print field_map.z_values()
 
22
# Calculate surface charge and current distributions from field map
 
23
# print field_map.data
 
24
# print field_map.Ez(-0.81, 0.228600)
 
25
 
 
26
"""
 
27
file = open('fields_values.txt', 'w')
 
28
for r in field_map.r_values():
 
29
  for z in field_map.z_values():
 
30
    file.write("".join((str(z), ' ', str(r), ' ', str(field_map.H(z, r)), '\n')))
 
31
file.close()
 
32
"""
 
33
inner_radius = 11.43  # cm (radius of inner wall of cavity)
 
34
donut_radius = 3.853  # cm
 
35
inner_height = 8.128  # cm (height of center flange)
 
36
z_offset = inner_height/2.0
 
37
 
 
38
z_values = field_map.z_values()-z_offset;
 
39
r_values = field_map.r_values();
 
40
 
 
41
print ''.join(('delta_z: ', str(z_values[1]-z_values[0])))
 
42
print ''.join(('delta_r: ', str(r_values[1]-r_values[0])))
 
43
"""
 
44
end_wall_r_values = r_values[bisect_right(r_values, button_radius):]
 
45
button_z_values = z_values[:bisect_right(z_values,
 
46
                                         -inner_height/2+button_base_height)]
 
47
"""
 
48
 
 
49
############## Calculate Normal Forces ##############
 
50
 
 
51
# Ez and Er are in MV/m. H is in A/m
 
52
f_E_w = np.zeros(len(r_values))  # electric force on wall
 
53
f_H_w = np.zeros(len(r_values))  # magnetic force on wall
 
54
f_Ez_d = np.zeros(len(r_values))  # z-axis electric force on donut
 
55
f_Er_d = np.zeros(len(r_values))  # r-axis electric force on donut
 
56
# c = 2.998e10  # cm/s
 
57
epsilon_0 = 8.85418782e-12  # m^-3 kg^-1 s^4 A^2
 
58
# epsilon_r = 9.84  # 98.5% purity alumina
 
59
epsilon_r = 9.57  # 99.5% purity alumina
 
60
mu_0 = 4*math.pi*1e-7  # kg m s-2 A^-2
 
61
for index,r in enumerate(r_values):
 
62
  # Electric normal force on end walls
 
63
  field_value = field_map.E(z_values[0]+z_offset, r)
 
64
  f_E_w[index] = epsilon_0/2 * (field_value*1e6)**2
 
65
 
 
66
  # Magnetic normal force on end walls
 
67
  field_value = field_map.H(z_values[0]+z_offset, r)
 
68
  f_H_w[index] = -mu_0/2. * field_value**2
 
69
 
 
70
  last_z = 0.0
 
71
  last_E = field_map.E(z_values[5]+z_offset, r)
 
72
  for z in z_values[5:]:
 
73
    current_E = field_map.E(z+z_offset, r)
 
74
    #print abs(current_E-last_E)/current_E
 
75
    if (abs(current_E-last_E)/current_E) > 0.25:
 
76
      # print ''.join(('z: ', str(last_z), '\tr: ', str(r)))
 
77
      # last_E is the surface field entering dielectric
 
78
      Ez_out = field_map.Ez(last_z+z_offset, r)*1e6
 
79
      Ez_in = field_map.Ez(z+z_offset, r)*1e6  # approximately
 
80
      Er_out = field_map.Er(last_z+z_offset, r)*1e6
 
81
      Er_in = field_map.Er(z+z_offset, r)*1e6  # approximately
 
82
      sigma = -epsilon_0 * math.sqrt((Ez_out-Ez_in)**2 + (Er_out-Er_in)**2)
 
83
      f_Ez_d[index] = 0.5 * sigma * Ez_out
 
84
      f_Er_d[index] = 0.5 * sigma * Er_out
 
85
      break
 
86
    last_z = z
 
87
    last_E = current_E
 
88
 
 
89
f_E_c = np.zeros(len(z_values))
 
90
f_H_c = np.zeros(len(z_values))
 
91
f_E_b = np.zeros(len(z_values))
 
92
f_H_b = np.zeros(len(z_values))
 
93
g_c = np.zeros(len(z_values))
 
94
g_b = np.zeros(len(z_values))
 
95
sigma_c = 1.6
 
96
sigma_b = 1.3
 
97
mu = 0
 
98
y_mirror = 0.0
 
99
for index,z in enumerate(z_values):
 
100
  # Electric normal force on circumference
 
101
  field_value = field_map.E(z+z_offset, r_values[::-1][0])
 
102
  f_E_c[index] = epsilon_0/2 * (field_value*1e6)**2
 
103
 
 
104
  # Magnetic normal force on circumference
 
105
  field_value = field_map.H(z+z_offset, r_values[::-1][0])
 
106
  f_H_c[index] = -mu_0/2. * field_value**2
 
107
 
 
108
  # Electric normal force on button
 
109
  """
 
110
  for r in r_values:
 
111
    field_value = field_map.E(z+z_offset, r)
 
112
    if (field_value > 0.):
 
113
      #       SI
 
114
      # f_E_r = -epsilon_0/2 * E^2
 
115
      f_E_b[index] = epsilon_0/2. * (field_value*1e6)**2
 
116
  """
 
117
  """
 
118
        The field in the gap between the buttons falls off a bit towards the
 
119
        center. This makes it difficult to get a good polynomial fit. If we
 
120
        flip the dip around, this provides a nice cap on the curve that
 
121
        produces a nice fit in the region we care about (on the buttons).
 
122
  """
 
123
  """
 
124
      print ''.join(("z: ", str(z)))
 
125
      if (abs(abs(z)-0.8128) < 1e-12):
 
126
        # save the y value right where the gap starts to use as an offset
 
127
        gap_y_offset_E_b = f_E_b[z_index]
 
128
        print ''.join(("gap_y_offset_E_b: ", str(gap_y_offset_E_b)))
 
129
      if (abs(z) <= 0.8128):
 
130
        # in the gap we flip the curve about the x axis and then slide
 
131
        # it back into place on the original side of the axis.
 
132
        f_E_b[z_index] = -f_E_b[z_index]+2*gap_y_offset_E_b
 
133
  """
 
134
 
 
135
############## Plynomial Curve Fitting ##############
 
136
# Electric normal pressure on end wall
 
137
orders = np.array((12, 12, 12, 12, 12, 12, 12, 12))
 
138
 
 
139
"""
 
140
domains = np.array((map(partial(add, r_values[22]), r_values),
 
141
                    z_values, z_values,
 
142
                    map(partial(add, r_values[22]), r_values),
 
143
                    z_values, z_values))
 
144
"""
 
145
domains = np.array((r_values, z_values, r_values,
 
146
                    r_values, z_values, r_values,
 
147
                    r_values, z_values))
 
148
# print r_values
 
149
"""
 
150
forces = np.array((np.roll(f_E_w, -23)[:-21], f_E_c, f_E_b,
 
151
                   np.roll(f_H_w, -23)[:-21], f_H_c, f_H_b))
 
152
"""
 
153
forces = np.array((f_E_w, f_E_c, f_Ez_d,
 
154
                   f_H_w, f_H_c, f_Er_d,
 
155
                   f_E_w+f_H_w, f_E_c+f_H_c))
 
156
segments_sets = [[[0,len(z_values)]],
 
157
                [[0,int(math.floor(len(z_values)/2.)-10)],
 
158
                 [int(math.floor(len(r_values)/2.)-10),int(math.floor(len(z_values)/2.)+11)],
 
159
                 [int(math.floor(len(r_values)/2.)+11),len(r_values)]],
 
160
                [[0,34], [34,42], [42,72], [72,79], [79,101]],
 
161
                [[0,len(z_values)]],
 
162
                [[0,len(r_values)]],
 
163
                [[0,34], [34,42], [42,72], [72,79], [79,101]],
 
164
                [[0,len(z_values)]],
 
165
                [[0,int(math.floor(len(z_values)/2.)-10)],
 
166
                 [int(math.floor(len(r_values)/2.)-10),int(math.floor(len(z_values)/2.)+11)],
 
167
                 [int(math.floor(len(r_values)/2.)+11),len(r_values)]]]
 
168
print segments_sets
 
169
for index,segments in enumerate(segments_sets):
 
170
  message = "["
 
171
  for segment in segments:
 
172
    message = ''.join((message, str([domains[index][segment[0]],
 
173
                                     domains[index][segment[1]-1]]),
 
174
                       ' '))
 
175
  print ''.join((message, ']'))
 
176
# print "".join(("Shape of forces: ", str(np.shape(forces))))
 
177
names = ("f_E_w", "f_E_c", "f_Ez_d",
 
178
         "f_H_w", "f_H_c", "f_Er_d",
 
179
         "f_E_w+f_H_w", "f_E_c+f_H_c")
 
180
coefficient_sets = map(fit_spliced_polynomials, domains, forces, orders,
 
181
                       segments_sets)
 
182
for i in range(np.shape(forces)[0]):
 
183
  coefficient_set = coefficient_sets[i]
 
184
  coeff_strings = []
 
185
  for coefficients in coefficient_set:
 
186
    coeff_string = ','.join(map(str, coefficients))
 
187
    coeff_strings += [''.join(('[', coeff_string, ']'))]
 
188
  message = "".join(("Polynomial coefficients for ", names[i], ":\n\t",
 
189
                     ",".join(coeff_strings)))
 
190
  print message
 
191
print len(segments_sets)
 
192
print len(coefficient_sets)
 
193
fits = map(lambda coefficient_set, segments, domain:
 
194
            polynomials(segments, coefficient_set, domain),
 
195
            coefficient_sets, segments_sets, domains)
 
196
##############   Plotting    ##############
 
197
fig = plt.figure(figsize=(16.0, 9.0))
 
198
 
 
199
i = int(0)
 
200
plt.subplot(3, 3, i+1)
 
201
#plt.plot(domains[i], forces[i]+forces[i+3])
 
202
plt.plot(domains[i], forces[i])
 
203
plt.plot(domains[i], fits[i])
 
204
plt.title('Normal Electric Pressure on End Wall')
 
205
plt.xlabel('Radius (cm)')
 
206
plt.ylabel('Pressure (N/m^2)')
 
207
#plt.plot(domains[i+3], np.abs(forces[i+3]))
 
208
#plt.plot(domains[i+3], -fits[i+3])
 
209
 
 
210
i += 1
 
211
plt.subplot(3, 3, i+1)
 
212
#plt.plot(domains[i], forces[i]+forces[i+3])
 
213
plt.plot(domains[i], forces[i])
 
214
plt.plot(domains[i], fits[i])
 
215
plt.title('Normal Electric Pressure on Circumference')
 
216
plt.xlabel('z (cm)')
 
217
plt.ylabel('Pressure (N/m^2)')
 
218
#plt.plot(domains[i+3], np.abs(forces[i+3]))
 
219
 
 
220
i += 1
 
221
plt.subplot(3, 3, i+1)
 
222
#plt.plot(domains[i], forces[i]+forces[i+3])
 
223
plt.plot(domains[i], forces[i])
 
224
plt.plot(domains[i], fits[i])
 
225
plt.title('Axial Electric Pressure on Donut Surface')
 
226
plt.xlabel('Radius (cm)')
 
227
plt.ylabel('Pressure (N/m^2)')
 
228
#plt.plot(domains[i+3], np.abs(forces[i+3]))
 
229
 
 
230
i += 1
 
231
plt.subplot(3, 3, i+1)
 
232
plt.plot(domains[i], forces[i])
 
233
plt.plot(domains[i], fits[i])
 
234
plt.title('Normal Magnetic Pressure on End Wall')
 
235
plt.xlabel('Radius (cm)')
 
236
plt.ylabel('Pressure (N/m^2)')
 
237
 
 
238
i += 1
 
239
plt.subplot(3, 3, i+1)
 
240
plt.plot(domains[i], forces[i])
 
241
plt.plot(domains[i], fits[i])
 
242
plt.title('Normal Magnetic Pressure on Circumference')
 
243
plt.xlabel('z (cm)')
 
244
plt.ylabel('Pressure (N/m^2)')
 
245
 
 
246
i += 1
 
247
plt.subplot(3, 3, i+1)
 
248
plt.plot(domains[i], forces[i])
 
249
plt.plot(domains[i], fits[i])
 
250
plt.title('Radial Electric Pressure on Donut')
 
251
plt.xlabel('Radius (cm)')
 
252
plt.ylabel('Pressure (N/m^2)')
 
253
 
 
254
i += 1
 
255
plt.subplot(3, 3, i+1)
 
256
plt.plot(domains[i], forces[i])
 
257
plt.plot(domains[i], fits[i])
 
258
plt.title('Normal E&M Pressure on End Wall')
 
259
plt.xlabel('Radius (cm)')
 
260
plt.ylabel('Pressure (N/m^2)')
 
261
 
 
262
i += 1
 
263
plt.subplot(3, 3, i+1)
 
264
plt.plot(domains[i], forces[i])
 
265
plt.plot(domains[i], fits[i])
 
266
plt.title('Normal E&M Pressure on Circumference')
 
267
plt.xlabel('z (cm)')
 
268
plt.ylabel('Pressure (N/m^2)')
 
269
"""
 
270
plt.plot(domains[7], forces[7])
 
271
plt.plot(domains[7], fits[7])
 
272
plt.title('Normal E&M Pressure on Circumference')
 
273
plt.xlabel('z (cm)')
 
274
plt.ylabel('Pressure (N/m^2)')
 
275
"""
 
276
plt.subplots_adjust(hspace=0.35, wspace=0.4)
 
277
 
 
278
plt.show()
 
 
b'\\ No newline at end of file'