| 1 |
gezelter |
4407 |
#!/opt/local/bin/python |
| 2 |
|
|
|
| 3 |
|
|
''' |
| 4 |
|
|
VectorFieldPlot - plots electric and magnetic fieldlines in svg |
| 5 |
|
|
http://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot |
| 6 |
|
|
|
| 7 |
|
|
Copyright (C) 2010 Geek3 |
| 8 |
|
|
|
| 9 |
|
|
This program is free software; you can redistribute it and/or modify |
| 10 |
|
|
it under the terms of the GNU General Public License as published by |
| 11 |
|
|
the Free Software Foundation; |
| 12 |
|
|
either version 3 of the License, or (at your option) any later version. |
| 13 |
|
|
|
| 14 |
|
|
This program is distributed in the hope that it will be useful, |
| 15 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 16 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. |
| 17 |
|
|
See the GNU General Public License for more details. |
| 18 |
|
|
|
| 19 |
|
|
You should have received a copy of the GNU General Public License |
| 20 |
|
|
along with this program; if not, see http://www.gnu.org/licenses/ |
| 21 |
|
|
''' |
| 22 |
|
|
|
| 23 |
|
|
version = '1.3' |
| 24 |
|
|
|
| 25 |
|
|
|
| 26 |
|
|
from math import * |
| 27 |
|
|
from lxml import etree |
| 28 |
|
|
import scipy as sc |
| 29 |
|
|
import scipy.optimize as op |
| 30 |
|
|
import scipy.integrate as ig |
| 31 |
|
|
import bisect |
| 32 |
|
|
|
| 33 |
|
|
|
| 34 |
|
|
|
| 35 |
|
|
# some helper functions |
| 36 |
|
|
def vabs(x): |
| 37 |
|
|
''' |
| 38 |
|
|
euclidian vector norm for any kind of vector |
| 39 |
|
|
''' |
| 40 |
|
|
return sqrt(sum([i**2 for i in x])) |
| 41 |
|
|
|
| 42 |
|
|
def vnorm(x): |
| 43 |
|
|
''' |
| 44 |
|
|
vector normalisation |
| 45 |
|
|
''' |
| 46 |
|
|
d = vabs(x) |
| 47 |
|
|
if d != 0.: return sc.array(x) / vabs(x) |
| 48 |
|
|
return sc.array(x) |
| 49 |
|
|
|
| 50 |
|
|
def rot(xy, phi): |
| 51 |
|
|
''' |
| 52 |
|
|
2D vector rotation |
| 53 |
|
|
''' |
| 54 |
|
|
s = sin(phi); c = cos(phi) |
| 55 |
|
|
return sc.array([c * xy[0] - s * xy[1], c * xy[1] + s * xy[0]]) |
| 56 |
|
|
|
| 57 |
|
|
def cosv(v1, v2): |
| 58 |
|
|
''' |
| 59 |
|
|
find the cosine of the angle between two vectors |
| 60 |
|
|
''' |
| 61 |
|
|
d1 = sum(v1**2); d2 = sum(v2**2) |
| 62 |
|
|
if d1 != 1.: d1 = sqrt(d1) |
| 63 |
|
|
if d2 != 1.: d2 = sqrt(d2) |
| 64 |
|
|
if d1 * d2 == 0.: return 1. |
| 65 |
|
|
return sc.dot(v1, v2) / (d1 * d2) |
| 66 |
|
|
|
| 67 |
|
|
def sinv(v1, v2): |
| 68 |
|
|
''' |
| 69 |
|
|
find the sine of the angle between two vectors |
| 70 |
|
|
''' |
| 71 |
|
|
d1 = sum(v1**2); d2 = sum(v2**2) |
| 72 |
|
|
if d1 != 1.: d1 = sqrt(d1) |
| 73 |
|
|
if d2 != 1.: d2 = sqrt(d2) |
| 74 |
|
|
if d1 * d2 == 0.: return 0. |
| 75 |
|
|
return (v1[0] * v2[1] - v1[1] * v2[0]) / (d1 * d2) |
| 76 |
|
|
|
| 77 |
|
|
def angle_dif(a1, a2): |
| 78 |
|
|
return ((a2 - a1 + pi) % (2. * pi)) - pi |
| 79 |
|
|
|
| 80 |
|
|
def list_interpolate(l, t): |
| 81 |
|
|
n = max(0, bisect.bisect_right(l, t) - 1) |
| 82 |
|
|
s = None |
| 83 |
|
|
if t < l[0]: |
| 84 |
|
|
if l[1] - l[0] == 0.: |
| 85 |
|
|
s = 0. |
| 86 |
|
|
else: |
| 87 |
|
|
s = (t - l[0]) / float(l[1] - l[0]) |
| 88 |
|
|
elif t >= l[-1]: |
| 89 |
|
|
n = len(l) - 2 |
| 90 |
|
|
if l[-1] - l[-2] == 0.: |
| 91 |
|
|
s = 1. |
| 92 |
|
|
else: |
| 93 |
|
|
s = (t - l[-2]) / float(l[-1] - l[-2]) |
| 94 |
|
|
else: |
| 95 |
|
|
s = (t - l[n]) / (l[n+1] - l[n]) |
| 96 |
|
|
return n, s |
| 97 |
|
|
|
| 98 |
|
|
def pretty_vec(p): |
| 99 |
|
|
return '{0:> 9.5f},{1:> 9.5f}'.format(p[0], p[1]) |
| 100 |
|
|
|
| 101 |
|
|
|
| 102 |
|
|
class FieldplotDocument: |
| 103 |
|
|
''' |
| 104 |
|
|
creates a svg document structure using lxml.etree |
| 105 |
|
|
''' |
| 106 |
|
|
def __init__ (self, name, width=800, height=600, digits=3.5, unit=100, |
| 107 |
|
|
center=None, licence='GFDL-cc', commons=False): |
| 108 |
|
|
self.name = name |
| 109 |
|
|
self.width = float(width) |
| 110 |
|
|
self.height = float(height) |
| 111 |
|
|
self.digits = float(digits) |
| 112 |
|
|
self.unit = float(unit) |
| 113 |
|
|
self.licence = licence |
| 114 |
|
|
self.commons = commons |
| 115 |
|
|
if center == None: self.center = [width / 2., height / 2.] |
| 116 |
|
|
else: self.center = [float(i) for i in center] |
| 117 |
|
|
|
| 118 |
|
|
# create document structure |
| 119 |
|
|
self.svg = etree.Element('svg', |
| 120 |
|
|
nsmap={None: 'http://www.w3.org/2000/svg', |
| 121 |
|
|
'xlink': 'http://www.w3.org/1999/xlink'}) |
| 122 |
|
|
self.svg.set('version', '1.1') |
| 123 |
|
|
self.svg.set('baseProfile', 'full') |
| 124 |
|
|
self.svg.set('width', str(int(width))) |
| 125 |
|
|
self.svg.set('height', str(int(height))) |
| 126 |
|
|
|
| 127 |
|
|
# title |
| 128 |
|
|
self.title = etree.SubElement(self.svg, 'title') |
| 129 |
|
|
self.title.text = self.name |
| 130 |
|
|
|
| 131 |
|
|
# description |
| 132 |
|
|
self.desc = etree.SubElement(self.svg, 'desc') |
| 133 |
|
|
self.desc.text = '' |
| 134 |
|
|
self.desc.text += self.name + '\n' |
| 135 |
|
|
self.desc.text += 'created with VectorFieldPlot ' + version + '\n' |
| 136 |
|
|
self.desc.text += 'http://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot\n' |
| 137 |
|
|
if commons: |
| 138 |
|
|
self.desc.text += """ |
| 139 |
|
|
about: http://commons.wikimedia.org/wiki/File:{0}.svg |
| 140 |
|
|
""".format(self.name) |
| 141 |
|
|
if self.licence == 'GFDL-cc': |
| 142 |
|
|
self.desc.text += """rights: GNU Free Documentation license, |
| 143 |
|
|
Creative Commons Attribution ShareAlike license\n""" |
| 144 |
|
|
self.desc.text += ' ' |
| 145 |
|
|
|
| 146 |
|
|
# background |
| 147 |
|
|
self.background = etree.SubElement(self.svg, 'rect') |
| 148 |
|
|
self.background.set('id', 'background') |
| 149 |
|
|
self.background.set('x', '0') |
| 150 |
|
|
self.background.set('y', '0') |
| 151 |
|
|
self.background.set('width', str(width)) |
| 152 |
|
|
self.background.set('height', str(height)) |
| 153 |
|
|
self.background.set('fill', '#ffffff') |
| 154 |
|
|
|
| 155 |
|
|
# image elements |
| 156 |
|
|
self.content = etree.SubElement(self.svg, 'g') |
| 157 |
|
|
self.content.set('id', 'image') |
| 158 |
|
|
self.content.set('transform', |
| 159 |
|
|
'translate({0},{1}) scale({2},-{2})'.format( |
| 160 |
|
|
self.center[0], self.center[1], self.unit)) |
| 161 |
|
|
|
| 162 |
|
|
self.arrow_geo = {'x_nock':0.3,'x_head':3.8,'x_tail':-2.2,'width':4.5} |
| 163 |
|
|
|
| 164 |
|
|
def __get_defs(self): |
| 165 |
|
|
if 'defs' not in dir(self): |
| 166 |
|
|
self.defs = etree.Element('defs') |
| 167 |
|
|
self.desc.addnext(self.defs) |
| 168 |
|
|
return self.defs |
| 169 |
|
|
|
| 170 |
|
|
def __check_fieldlines(self, linecolor='#000000', linewidth=1.): |
| 171 |
|
|
if 'fieldlines' not in dir(self): |
| 172 |
|
|
self.fieldlines = etree.SubElement(self.content, 'g') |
| 173 |
|
|
self.fieldlines.set('id', 'fieldlines') |
| 174 |
|
|
self.fieldlines.set('fill', 'none') |
| 175 |
|
|
self.fieldlines.set('stroke', linecolor) |
| 176 |
|
|
self.fieldlines.set('stroke-width', |
| 177 |
|
|
str(linewidth / self.unit)) |
| 178 |
|
|
self.fieldlines.set('stroke-linejoin', 'round') |
| 179 |
|
|
self.fieldlines.set('stroke-linecap', 'round') |
| 180 |
|
|
if 'count_fieldlines' not in dir(self): self.count_fieldlines = 0 |
| 181 |
|
|
|
| 182 |
|
|
def __check_symbols(self): |
| 183 |
|
|
if 'symbols' not in dir(self): |
| 184 |
|
|
self.symbols = etree.SubElement(self.content, 'g') |
| 185 |
|
|
self.symbols.set('id', 'symbols') |
| 186 |
|
|
if 'count_symbols' not in dir(self): self.count_symbols = 0 |
| 187 |
|
|
|
| 188 |
|
|
def __check_whitespot(self): |
| 189 |
|
|
if 'whitespot' not in dir(self): |
| 190 |
|
|
self.whitespot = etree.SubElement( |
| 191 |
|
|
self.__get_defs(), 'radialGradient') |
| 192 |
|
|
self.whitespot.set('id', 'white_spot') |
| 193 |
|
|
for attr, val in [['cx', '0.65'], ['cy', '0.7'], ['r', '0.75']]: |
| 194 |
|
|
self.whitespot.set(attr, val) |
| 195 |
|
|
for col, of, op in [['#ffffff', '0', '0.7'], |
| 196 |
|
|
['#ffffff', '0.1', '0.5'], ['#ffffff', '0.6', '0'], |
| 197 |
|
|
['#000000', '0.6', '0'], ['#000000', '0.75', '0.05'], |
| 198 |
|
|
['#000000', '0.85', '0.15'], ['#000000', '1', '0.5']]: |
| 199 |
|
|
stop = etree.SubElement(self.whitespot, 'stop') |
| 200 |
|
|
stop.set('stop-color', col) |
| 201 |
|
|
stop.set('offset', of) |
| 202 |
|
|
stop.set('stop-opacity', op) |
| 203 |
|
|
|
| 204 |
|
|
def __get_arrowname(self, fillcolor='#000000'): |
| 205 |
|
|
if 'arrows' not in dir(self): |
| 206 |
|
|
self.arrows = {} |
| 207 |
|
|
if fillcolor not in self.arrows.iterkeys(): |
| 208 |
|
|
arrow = etree.SubElement(self.__get_defs(), 'path') |
| 209 |
|
|
self.arrows[fillcolor] = arrow |
| 210 |
|
|
arrow.set('id', 'arrow' + str(len(self.arrows))) |
| 211 |
|
|
arrow.set('stroke', 'none') |
| 212 |
|
|
arrow.set('fill', fillcolor) |
| 213 |
|
|
arrow.set('transform', 'scale({0})'.format(1. / self.unit)) |
| 214 |
|
|
arrow.set('d', |
| 215 |
|
|
'M {0},0 L {1},{3} L {2},0 L {1},-{3} L {0},0 Z'.format( |
| 216 |
|
|
self.arrow_geo['x_nock'], self.arrow_geo['x_tail'], |
| 217 |
|
|
self.arrow_geo['x_head'], self.arrow_geo['width'] / 2.)) |
| 218 |
|
|
return self.arrows[fillcolor].get('id') |
| 219 |
|
|
|
| 220 |
|
|
def draw_charges(self, field, scale=1.): |
| 221 |
|
|
if 'monopoles' not in field.elements: return |
| 222 |
|
|
charges = field.elements['monopoles'] |
| 223 |
|
|
self.__check_symbols() |
| 224 |
|
|
self.__check_whitespot() |
| 225 |
|
|
|
| 226 |
|
|
for charge in charges: |
| 227 |
|
|
c_group = etree.SubElement(self.symbols, 'g') |
| 228 |
|
|
self.count_symbols += 1 |
| 229 |
|
|
c_group.set('id', 'charge{0}'.format(self.count_symbols)) |
| 230 |
|
|
c_group.set('transform', |
| 231 |
|
|
'translate({0},{1}) scale({2},{2})'.format( |
| 232 |
|
|
charge[0], charge[1], float(scale) / self.unit)) |
| 233 |
|
|
|
| 234 |
|
|
#### charge drawing #### |
| 235 |
|
|
c_bg = etree.SubElement(c_group, 'circle') |
| 236 |
|
|
c_shade = etree.SubElement(c_group, 'circle') |
| 237 |
|
|
c_symb = etree.SubElement(c_group, 'path') |
| 238 |
|
|
if charge[2] >= 0.: c_bg.set('style', 'fill:#ff0000; stroke:none') |
| 239 |
|
|
else: c_bg.set('style', 'fill:#0000ff; stroke:none') |
| 240 |
|
|
for attr, val in [['cx', '0'], ['cy', '0'], ['r', '14']]: |
| 241 |
|
|
c_bg.set(attr, val) |
| 242 |
|
|
c_shade.set(attr, val) |
| 243 |
|
|
c_shade.set('style', |
| 244 |
|
|
'fill:url(#white_spot); stroke:#000000; stroke-width:2') |
| 245 |
|
|
# plus sign |
| 246 |
|
|
if charge[2] >= 0.: |
| 247 |
|
|
c_symb.set('d', 'M 2,2 V 8 H -2 V 2 H -8 V -2' |
| 248 |
|
|
+ ' H -2 V -8 H 2 V -2 H 8 V 2 H 2 Z') |
| 249 |
|
|
# minus sign |
| 250 |
|
|
else: c_symb.set('d', 'M 8,2 H -8 V -2 H 8 V 2 Z') |
| 251 |
|
|
c_symb.set('style', 'fill:#000000; stroke:none') |
| 252 |
|
|
|
| 253 |
|
|
def draw_currents(self, field, scale=1.): |
| 254 |
|
|
if ('wires' not in field.elements |
| 255 |
|
|
and 'ringcurrents' not in field.elements): |
| 256 |
|
|
return |
| 257 |
|
|
self.__check_symbols() |
| 258 |
|
|
self.__check_whitespot() |
| 259 |
|
|
currents = [] |
| 260 |
|
|
if 'wires' in field.elements: |
| 261 |
|
|
for i in field.elements['wires']: |
| 262 |
|
|
currents.append(i) |
| 263 |
|
|
if 'ringcurrents' in field.elements: |
| 264 |
|
|
for i in field.elements['ringcurrents']: |
| 265 |
|
|
currents.append(list(i[:2] + rot([0., i[3]], i[2])) + [i[-1]]) |
| 266 |
|
|
currents.append(list(i[:2] - rot([0., i[3]], i[2])) + [-i[-1]]) |
| 267 |
|
|
|
| 268 |
|
|
for cur in currents: |
| 269 |
|
|
c_group = etree.SubElement(self.symbols, 'g') |
| 270 |
|
|
self.count_symbols += 1 |
| 271 |
|
|
if cur[-1] >= 0.: direction = 'out' |
| 272 |
|
|
else: direction = 'in' |
| 273 |
|
|
c_group.set('id', |
| 274 |
|
|
'current_{0}{1}'.format(direction, self.count_symbols)) |
| 275 |
|
|
c_group.set('transform', |
| 276 |
|
|
'translate({0},{1}) scale({2},{2})'.format( |
| 277 |
|
|
cur[0], cur[1], float(scale) / self.unit)) |
| 278 |
|
|
|
| 279 |
|
|
#### current drawing #### |
| 280 |
|
|
c_bg = etree.SubElement(c_group, 'circle') |
| 281 |
|
|
c_shade = etree.SubElement(c_group, 'circle') |
| 282 |
|
|
c_bg.set('style', 'fill:#b0b0b0; stroke:none') |
| 283 |
|
|
for attr, val in [['cx', '0'], ['cy', '0'], ['r', '14']]: |
| 284 |
|
|
c_bg.set(attr, val) |
| 285 |
|
|
c_shade.set(attr, val) |
| 286 |
|
|
c_shade.set('style', |
| 287 |
|
|
'fill:url(#white_spot); stroke:#000000; stroke-width:2') |
| 288 |
|
|
if cur[-1] >= 0.: # dot |
| 289 |
|
|
c_symb = etree.SubElement(c_group, 'circle') |
| 290 |
|
|
c_symb.set('cx', '0') |
| 291 |
|
|
c_symb.set('cy', '0') |
| 292 |
|
|
c_symb.set('r', '4') |
| 293 |
|
|
else: # cross |
| 294 |
|
|
c_symb = etree.SubElement(c_group, 'path') |
| 295 |
|
|
c_symb.set('d', 'M {1},-{0} L {0},-{1} L {2},{3} L {0},{1} \ |
| 296 |
|
|
L {1},{0} {3},{2} L -{1},{0} L -{0},{1} L -{2},{3} L -{0},-{1} L -{1},-{0} \ |
| 297 |
|
|
L {3},-{2} L {1},-{0} Z'.format(11.1, 8.5, 2.6, 0)) |
| 298 |
|
|
c_symb.set('style', 'fill:#000000; stroke:none') |
| 299 |
|
|
|
| 300 |
|
|
def draw_magnets(self, field): |
| 301 |
|
|
if 'coils' not in field.elements: return |
| 302 |
|
|
coils = field.elements['coils'] |
| 303 |
|
|
self.__check_symbols() |
| 304 |
|
|
|
| 305 |
|
|
for coil in coils: |
| 306 |
|
|
m_group = etree.SubElement(self.symbols, 'g') |
| 307 |
|
|
self.count_symbols += 1 |
| 308 |
|
|
m_group.set('id', 'magnet{0}'.format(self.count_symbols)) |
| 309 |
|
|
m_group.set('transform', |
| 310 |
|
|
'translate({0},{1}) rotate({2})'.format( |
| 311 |
|
|
coil[0], coil[1], degrees(coil[2]))) |
| 312 |
|
|
|
| 313 |
|
|
#### magnet drawing #### |
| 314 |
|
|
r = coil[3]; l = coil[4] |
| 315 |
|
|
colors = ['#00cc00', '#ff0000'] |
| 316 |
|
|
SN = ['S', 'N'] |
| 317 |
|
|
if coil[5] < 0.: |
| 318 |
|
|
colors.reverse() |
| 319 |
|
|
SN.reverse() |
| 320 |
|
|
m_defs = etree.SubElement(m_group, 'defs') |
| 321 |
|
|
m_gradient = etree.SubElement(m_defs, 'linearGradient') |
| 322 |
|
|
m_gradient.set('id', 'magnetGrad{0}'.format(self.count_symbols)) |
| 323 |
|
|
for attr, val in [['x1', '0'], ['x2', '0'], ['y1', str(coil[3])], |
| 324 |
|
|
['y2', str(-coil[3])], ['gradientUnits', 'userSpaceOnUse']]: |
| 325 |
|
|
m_gradient.set(attr, val) |
| 326 |
|
|
for col, of, opa in [['#000000', '0', '0.125'], |
| 327 |
|
|
['#ffffff', '0.07', '0.125'], ['#ffffff', '0.25', '0.5'], |
| 328 |
|
|
['#ffffff', '0.6', '0.2'], ['#000000', '1', '0.33']]: |
| 329 |
|
|
stop = etree.SubElement(m_gradient, 'stop') |
| 330 |
|
|
stop.set('stop-color', col) |
| 331 |
|
|
stop.set('offset', of) |
| 332 |
|
|
stop.set('stop-opacity', opa) |
| 333 |
|
|
for i in [0, 1]: |
| 334 |
|
|
rect = etree.SubElement(m_group, 'rect') |
| 335 |
|
|
for attr, val in [['x', [-l, 0][i]], ['y', -r], |
| 336 |
|
|
['width', [2*l, l][i]], ['height', 2 * r], |
| 337 |
|
|
['style', 'fill:{0}; stroke:none'.format(colors[i])]]: |
| 338 |
|
|
rect.set(attr, str(val)) |
| 339 |
|
|
rect = etree.SubElement(m_group, 'rect') |
| 340 |
|
|
for attr, val in [['x', -l], ['y', -r], |
| 341 |
|
|
['width', 2 * l], ['height', 2 * r], |
| 342 |
|
|
['style', 'fill:url(#magnetGrad{0}); stroke-width:{1}; stroke-linejoin:miter; stroke:#000000'.format(self.count_symbols, 4. / self.unit)]]: |
| 343 |
|
|
rect.set(attr, str(val)) |
| 344 |
|
|
for i in [0, 1]: |
| 345 |
|
|
text = etree.SubElement(m_group, 'text') |
| 346 |
|
|
for attr, val in [['text-anchor', 'middle'], ['y', -r], |
| 347 |
|
|
['transform', 'translate({0},{1}) scale({2},-{2})'.format( |
| 348 |
|
|
[-0.65, 0.65][i] * l, -0.44 * r, r / 100.)], |
| 349 |
|
|
['style', 'fill:#000000; stroke:none; ' + |
| 350 |
|
|
'font-size:120px; font-family:Bitstream Vera Sans']]: |
| 351 |
|
|
text.set(attr, str(val)) |
| 352 |
|
|
text.text = SN[i] |
| 353 |
|
|
|
| 354 |
|
|
def draw_line(self, fieldline, maxdist=10., linewidth=2., |
| 355 |
|
|
linecolor='#000000', attributes=[], arrows_style=None): |
| 356 |
|
|
''' |
| 357 |
|
|
draws a calculated fieldline from a FieldLine object |
| 358 |
|
|
to the FieldplotDocument svg image |
| 359 |
|
|
''' |
| 360 |
|
|
self.__check_fieldlines(linecolor, linewidth) |
| 361 |
|
|
self.count_fieldlines += 1 |
| 362 |
|
|
|
| 363 |
|
|
bounds = {} |
| 364 |
|
|
bounds['x0'] = -(self.center[0] + 0.5 * linewidth) / self.unit |
| 365 |
|
|
bounds['y0'] = -(self.height - self.center[1] + |
| 366 |
|
|
0.5 * linewidth) / self.unit |
| 367 |
|
|
bounds['x1'] = (self.width - self.center[0] + |
| 368 |
|
|
0.5 * linewidth) / self.unit |
| 369 |
|
|
bounds['y1'] = (self.center[1] + 0.5 * linewidth) / self.unit |
| 370 |
|
|
|
| 371 |
|
|
# fetch the polyline from the fieldline object |
| 372 |
|
|
polylines = fieldline.get_polylines(self.digits, maxdist, bounds) |
| 373 |
|
|
if len(polylines) == 0: return |
| 374 |
|
|
|
| 375 |
|
|
line = etree.Element('path') |
| 376 |
|
|
if self.fieldlines.get('stroke') != linecolor: |
| 377 |
|
|
line.set('stroke', linecolor) |
| 378 |
|
|
if self.fieldlines.get('stroke-width') != str(linewidth / self.unit): |
| 379 |
|
|
line.set('stroke-width', str(linewidth / self.unit)) |
| 380 |
|
|
for attr, val in attributes: |
| 381 |
|
|
line.set(attr, val) |
| 382 |
|
|
|
| 383 |
|
|
#### line drawing #### |
| 384 |
|
|
path_data = [] |
| 385 |
|
|
for polyline in polylines: |
| 386 |
|
|
line_points = polyline['path'] |
| 387 |
|
|
for i, p in enumerate(line_points): |
| 388 |
|
|
# go through all points, draw them if line segment is visible |
| 389 |
|
|
ptext = '{1:.{0}f},{2:.{0}f}'.format( |
| 390 |
|
|
int(ceil(self.digits)), p[0], p[1]) |
| 391 |
|
|
if i == 0: path_data.append('M ' + ptext) |
| 392 |
|
|
else: path_data.append('L ' + ptext) |
| 393 |
|
|
# close path if possible |
| 394 |
|
|
if (vabs(polylines[0]['path'][0] - polylines[-1]['path'][-1]) |
| 395 |
|
|
< .1**self.digits): |
| 396 |
|
|
closed = True |
| 397 |
|
|
if len(polylines) == 1: |
| 398 |
|
|
path_data.append('Z') |
| 399 |
|
|
elif len(polylines) > 1: |
| 400 |
|
|
# rearrange array cyclic |
| 401 |
|
|
path_data.pop(0) |
| 402 |
|
|
while path_data[0][0] != 'M': |
| 403 |
|
|
path_data.append(path_data.pop(0)) |
| 404 |
|
|
else: closed = False |
| 405 |
|
|
|
| 406 |
|
|
path = ' '.join(path_data) |
| 407 |
|
|
line.set('d', path) |
| 408 |
|
|
|
| 409 |
|
|
if arrows_style == None: |
| 410 |
|
|
# include path directly into document structure |
| 411 |
|
|
line.set('id', 'fieldline{0}'.format(self.count_fieldlines)) |
| 412 |
|
|
self.fieldlines.append(line) |
| 413 |
|
|
else: |
| 414 |
|
|
line_and_arrows = etree.SubElement(self.fieldlines, 'g') |
| 415 |
|
|
line_and_arrows.set('id', 'fieldline' + str(self.count_fieldlines)) |
| 416 |
|
|
line_and_arrows.append(line) |
| 417 |
|
|
line_and_arrows.append(self.__draw_arrows( |
| 418 |
|
|
arrows_style, linewidth, polylines, linecolor, closed)) |
| 419 |
|
|
|
| 420 |
|
|
def __draw_arrows(self, arrows_style, linewidth, polylines, |
| 421 |
|
|
linecolor='#000000', closed=False): |
| 422 |
|
|
''' |
| 423 |
|
|
draws arrows on polylines. |
| 424 |
|
|
options in "arrows_style": |
| 425 |
|
|
min_arrows: minimum number of arrows per segment |
| 426 |
|
|
max_arrows: maximum number of arrows per segment (None: no limit) |
| 427 |
|
|
dist: optimum distance between arrows |
| 428 |
|
|
scale: relative size of arrows to linewidth |
| 429 |
|
|
offsets [start_offset, mid_end, mid_start, end_offset] |
| 430 |
|
|
fixed_ends [True, False, False, True]: |
| 431 |
|
|
make first/last arrow distance invariable |
| 432 |
|
|
''' |
| 433 |
|
|
min_arrows = 1; max_arrows = None; arrows_dist = 1.; scale = linewidth |
| 434 |
|
|
offsets = 4 * [0.5]; fixed_ends = 4 * [False] |
| 435 |
|
|
if 'min_arrows' in arrows_style: |
| 436 |
|
|
min_arrows = arrows_style['min_arrows'] |
| 437 |
|
|
if 'max_arrows' in arrows_style: |
| 438 |
|
|
max_arrows = arrows_style['max_arrows'] |
| 439 |
|
|
if 'dist' in arrows_style: |
| 440 |
|
|
arrows_dist = arrows_style['dist'] |
| 441 |
|
|
if 'scale' in arrows_style: |
| 442 |
|
|
scale *= arrows_style['scale'] |
| 443 |
|
|
if 'offsets' in arrows_style: |
| 444 |
|
|
offsets = arrows_style['offsets'] |
| 445 |
|
|
if 'fixed_ends' in arrows_style: |
| 446 |
|
|
fixed_ends = arrows_style['fixed_ends'] |
| 447 |
|
|
if scale == 1.: scaletext = '' |
| 448 |
|
|
else: scaletext = ' scale({0})'.format(scale) |
| 449 |
|
|
|
| 450 |
|
|
arrows = etree.Element('g') |
| 451 |
|
|
arrows.set('id', 'arrows' + str(self.count_fieldlines)) |
| 452 |
|
|
for j, polyline in enumerate(polylines): |
| 453 |
|
|
line_points = polyline['path'] |
| 454 |
|
|
mina = min_arrows |
| 455 |
|
|
maxa = max_arrows |
| 456 |
|
|
# measure drawn path length |
| 457 |
|
|
lines_dist = [0.] |
| 458 |
|
|
for i in range(1, len(line_points)): |
| 459 |
|
|
lines_dist.append(lines_dist[-1] |
| 460 |
|
|
+ vabs(line_points[i] - line_points[i-1])) |
| 461 |
|
|
|
| 462 |
|
|
offs = [offsets[2], offsets[1]] |
| 463 |
|
|
fixed = [fixed_ends[2], fixed_ends[1]] |
| 464 |
|
|
if polyline['start']: |
| 465 |
|
|
offs[0] = offsets[0] |
| 466 |
|
|
fixed[0] = fixed_ends[0] |
| 467 |
|
|
if polyline['end']: |
| 468 |
|
|
offs[1] = offsets[3] |
| 469 |
|
|
fixed[1] = fixed_ends[3] |
| 470 |
|
|
|
| 471 |
|
|
d01 = [0., lines_dist[-1]] |
| 472 |
|
|
for i in [0, 1]: |
| 473 |
|
|
if fixed[i]: |
| 474 |
|
|
d01[i] += offs[i] * arrows_dist * [1., -1.][i] |
| 475 |
|
|
mina -= 1 |
| 476 |
|
|
if maxa != None: maxa -= 1 |
| 477 |
|
|
if d01[1] - d01[0] < 0.: break |
| 478 |
|
|
elif d01[1] - d01[0] == 0.: d_list = [d01[0]] |
| 479 |
|
|
else: |
| 480 |
|
|
d_list = [] |
| 481 |
|
|
if fixed[0]: d_list.append(d01[0]) |
| 482 |
|
|
if maxa > 0 or maxa == None: |
| 483 |
|
|
number_intervals = (d01[1] - d01[0]) / arrows_dist |
| 484 |
|
|
number_offsets = 0. |
| 485 |
|
|
for i in [0, 1]: |
| 486 |
|
|
if fixed[i]: number_offsets += .5 |
| 487 |
|
|
else: number_offsets += offs[i] - .5 |
| 488 |
|
|
n = int(number_intervals - number_offsets + 0.5) |
| 489 |
|
|
n = max(n, mina) |
| 490 |
|
|
if maxa != None: n = min(n, maxa) |
| 491 |
|
|
if n > 0: |
| 492 |
|
|
d = (d01[1] - d01[0]) / float(n + number_offsets) |
| 493 |
|
|
if fixed[0]: d_start = d01[0] + d |
| 494 |
|
|
else: d_start = offs[0] * d |
| 495 |
|
|
for i in range(n): |
| 496 |
|
|
d_list.append(d_start + i * d) |
| 497 |
|
|
if fixed[1]: d_list.append(d01[1]) |
| 498 |
|
|
|
| 499 |
|
|
geo = self.arrow_geo # shortcut |
| 500 |
|
|
#### arrow drawing #### |
| 501 |
|
|
for d1 in d_list: |
| 502 |
|
|
# calculate arrow position and direction |
| 503 |
|
|
if d1 < 0. or d1 > lines_dist[-1]: continue |
| 504 |
|
|
d0 = d1 + (geo['x_nock'] * scale + 2.5*linewidth * |
| 505 |
|
|
(geo['x_tail'] - geo['x_nock']) / geo['width']) / self.unit |
| 506 |
|
|
if closed and d0 < 0.: d0 = lines_dist[-1] + d0 |
| 507 |
|
|
d2 = d1 + (geo['x_head'] * scale + linewidth * |
| 508 |
|
|
(geo['x_tail'] - geo['x_head']) / geo['width']) / self.unit |
| 509 |
|
|
if closed and d2 > lines_dist[-1]: d1 -= lines_dist[-1] |
| 510 |
|
|
i0, s0 = list_interpolate(lines_dist, d0) |
| 511 |
|
|
i1, s1 = list_interpolate(lines_dist, d1) |
| 512 |
|
|
i2, s2 = list_interpolate(lines_dist, d2) |
| 513 |
|
|
p0 = line_points[i0] + s0 * (line_points[i0+1]-line_points[i0]) |
| 514 |
|
|
p1 = line_points[i1] + s1 * (line_points[i1+1]-line_points[i1]) |
| 515 |
|
|
p2 = line_points[i2] + s2 * (line_points[i2+1]-line_points[i2]) |
| 516 |
|
|
p = None; angle = None |
| 517 |
|
|
if vabs(p2-p1) <= .1**self.digits or (d2 <= d0 and not closed): |
| 518 |
|
|
v = line_points[i1+1] - line_points[i1] |
| 519 |
|
|
p = p1 |
| 520 |
|
|
angle = atan2(v[1], v[0]) |
| 521 |
|
|
else: |
| 522 |
|
|
v = p2 - p0 |
| 523 |
|
|
p = p0 + sc.dot(p1 - p0, v) * v / vabs(v)**2 |
| 524 |
|
|
angle = atan2(v[1], v[0]) |
| 525 |
|
|
|
| 526 |
|
|
arrow = etree.SubElement(arrows, 'use') |
| 527 |
|
|
arrow.set('{http://www.w3.org/1999/xlink}href', |
| 528 |
|
|
'#' + self.__get_arrowname(linecolor)) |
| 529 |
|
|
arrow.set('transform', ('translate({0:.' |
| 530 |
|
|
+ str(int(ceil(self.digits))) + 'f},{1:.' |
| 531 |
|
|
+ str(int(ceil(self.digits))) |
| 532 |
|
|
+ 'f}) rotate({2:.2f})').format(p[0], p[1], |
| 533 |
|
|
degrees(angle)) + scaletext) |
| 534 |
|
|
return arrows |
| 535 |
|
|
|
| 536 |
|
|
def draw_object(self, name, params, group=None): |
| 537 |
|
|
''' |
| 538 |
|
|
Draw arbitraty svg object. |
| 539 |
|
|
Params must be a dictionary of valid svg parameters. |
| 540 |
|
|
''' |
| 541 |
|
|
self.__check_symbols() |
| 542 |
|
|
if group == None: |
| 543 |
|
|
obj = etree.SubElement(self.symbols, name) |
| 544 |
|
|
else: |
| 545 |
|
|
obj = etree.SubElement(group, name) |
| 546 |
|
|
for i, j in params.iteritems(): |
| 547 |
|
|
obj.set(str(i), str(j)) |
| 548 |
|
|
return obj |
| 549 |
|
|
|
| 550 |
|
|
def write(self, filename=None): |
| 551 |
|
|
# put symbols on top |
| 552 |
|
|
if 'content' in dir(self): |
| 553 |
|
|
for element in self.content: |
| 554 |
|
|
if element.get('id').startswith('symbols'): |
| 555 |
|
|
self.content.append(element) |
| 556 |
|
|
|
| 557 |
|
|
# write content to file |
| 558 |
|
|
if filename == None: filename = self.name |
| 559 |
|
|
outfile = open(filename + '.svg', 'w') |
| 560 |
|
|
outfile.write(etree.tostring(self.svg, xml_declaration=True, |
| 561 |
|
|
pretty_print=True, encoding='utf-8')) |
| 562 |
|
|
outfile.close() |
| 563 |
|
|
print 'image written to', filename + '.svg' |
| 564 |
|
|
|
| 565 |
|
|
|
| 566 |
|
|
|
| 567 |
|
|
class FieldLine: |
| 568 |
|
|
''' |
| 569 |
|
|
calculates field lines |
| 570 |
|
|
''' |
| 571 |
|
|
def __init__(self, field, start_p, start_v=None, start_d=None, |
| 572 |
|
|
directions='forward', maxn=1000, maxr=300.0, hmax=1.0, |
| 573 |
|
|
pass_dipoles=0, path_close_tol=5e-3, bounds_func=None, |
| 574 |
|
|
stop_funcs=[None, None]): |
| 575 |
|
|
''' |
| 576 |
|
|
field: a field in which the line exists |
| 577 |
|
|
start_p: [x0, y0]: where the line starts |
| 578 |
|
|
start_v: [vx0, vy0]: optional start direction |
| 579 |
|
|
start_d: [dx0, dy0]: optional dipole start direction (slope to x=1) |
| 580 |
|
|
directions: forward, backward, both: bidirectional |
| 581 |
|
|
unit: estimation for the scale of the scene |
| 582 |
|
|
maxn: maximum number of steps |
| 583 |
|
|
maxr: maximum number of units to depart from start |
| 584 |
|
|
hmax: maximum number of units for stepsize |
| 585 |
|
|
pass_dipoles: number of dipoles to be passed through (-1 = infinite) |
| 586 |
|
|
''' |
| 587 |
|
|
self.field = field |
| 588 |
|
|
self.p_start = sc.array(start_p) |
| 589 |
|
|
self.first_point = self.p_start |
| 590 |
|
|
self.bounds_func = bounds_func |
| 591 |
|
|
self.stop_funcs = stop_funcs |
| 592 |
|
|
if start_v == None: self.v_start = None |
| 593 |
|
|
else: self.v_start = sc.array(start_v) |
| 594 |
|
|
if start_d == None: self.d_start = None |
| 595 |
|
|
else: self.d_start = sc.array(start_d) |
| 596 |
|
|
self.__create_nodes(directions, maxn, maxr, hmax, |
| 597 |
|
|
pass_dipoles, path_close_tol) |
| 598 |
|
|
|
| 599 |
|
|
def __get_nearest_pole(self, p, v=None): |
| 600 |
|
|
''' |
| 601 |
|
|
returns distance to nearest pole |
| 602 |
|
|
''' |
| 603 |
|
|
p_near = self.first_point |
| 604 |
|
|
d_near = vabs(self.first_point - p) |
| 605 |
|
|
if v != None: d_near *= 1.3 - cosv(v, self.first_point - p) |
| 606 |
|
|
type_near = 'start' |
| 607 |
|
|
mon = [] |
| 608 |
|
|
for ptype, poles in self.field.elements.iteritems(): |
| 609 |
|
|
if ptype not in ['monopoles', 'dipoles'] or len(poles) == 0: |
| 610 |
|
|
continue |
| 611 |
|
|
for pole in poles: |
| 612 |
|
|
d = vabs(pole[:2] - p) |
| 613 |
|
|
if v != None: d *= 1.3 - cosv(v, pole[:2] - p) |
| 614 |
|
|
if d < d_near: |
| 615 |
|
|
d_near = d |
| 616 |
|
|
p_near = pole |
| 617 |
|
|
type_near = ptype |
| 618 |
|
|
return sc.array(p_near), type_near |
| 619 |
|
|
|
| 620 |
|
|
def __rkstep(self, p, v, f, h): |
| 621 |
|
|
''' |
| 622 |
|
|
fourth order Runge Kutta step |
| 623 |
|
|
''' |
| 624 |
|
|
k1 = h * v |
| 625 |
|
|
k2 = h * f(p + k1 / 2.) |
| 626 |
|
|
k3 = h * f(p + k2 / 2.) |
| 627 |
|
|
k4 = h * f(p + k3) |
| 628 |
|
|
p1 = p + (k1 + 2. * (k2 + k3) + k4) / 6. |
| 629 |
|
|
return p1 |
| 630 |
|
|
|
| 631 |
|
|
def __create_nodes_part(self, sign, maxn, maxr, hmax, |
| 632 |
|
|
pass_dipoles, path_close_tol): |
| 633 |
|
|
''' |
| 634 |
|
|
executes integration from startpoint to one end |
| 635 |
|
|
''' |
| 636 |
|
|
# p is always the latest position |
| 637 |
|
|
# v is always the latest normalized velocity |
| 638 |
|
|
# h is always the latest step size |
| 639 |
|
|
# l is always the summarized length |
| 640 |
|
|
err = 5e-8 # error tolerance for integration |
| 641 |
|
|
f = None |
| 642 |
|
|
if sign >= 0.: f = self.field.Fn |
| 643 |
|
|
else: f = lambda r: -self.field.Fn(r) |
| 644 |
|
|
# first point |
| 645 |
|
|
p = self.p_start |
| 646 |
|
|
if self.v_start != None: |
| 647 |
|
|
v = vnorm(self.v_start) * sign |
| 648 |
|
|
else: |
| 649 |
|
|
v = f(p) |
| 650 |
|
|
nodes = [{'p':p.copy(), 'v_in':None}] |
| 651 |
|
|
xtol = 20. * err; ytol = path_close_tol |
| 652 |
|
|
# initialize loop |
| 653 |
|
|
h = (sqrt(5) - 1.) / 10.; h_old = h |
| 654 |
|
|
l = 0.; i = 0 |
| 655 |
|
|
while i < maxn and l < maxr: |
| 656 |
|
|
i += 1 |
| 657 |
|
|
if len(nodes) == 1 and self.d_start != None: |
| 658 |
|
|
# check for start from a dipole |
| 659 |
|
|
h = vabs(self.d_start) |
| 660 |
|
|
p = p + self.d_start |
| 661 |
|
|
v = f(p) |
| 662 |
|
|
nodes[-1]['v_out'] = h * vnorm(2.0 * vnorm(self.d_start) - v) |
| 663 |
|
|
nodes.append({'p':p.copy(), 'v_in':h * v}) |
| 664 |
|
|
elif len(nodes) > 1: |
| 665 |
|
|
# check for special cases |
| 666 |
|
|
nearest_pole, pole_type = self.__get_nearest_pole(p, v) |
| 667 |
|
|
vpole = nearest_pole[:2] - p |
| 668 |
|
|
dpole = vabs(vpole) |
| 669 |
|
|
vpole /= dpole |
| 670 |
|
|
|
| 671 |
|
|
cv = cosv(v, vpole); sv = sinv(v, vpole) |
| 672 |
|
|
if ((dpole < 0.1 or h >= dpole) |
| 673 |
|
|
and (cv > 0.9 or dpole < ytol)): |
| 674 |
|
|
# heading for some known special point |
| 675 |
|
|
if pole_type == 'start': |
| 676 |
|
|
# is the fieldline about to be closed? |
| 677 |
|
|
if ((dpole * abs(sv) < ytol) and |
| 678 |
|
|
(dpole * abs(cv) < xtol) and (l > 1e-3)): |
| 679 |
|
|
# path is closed |
| 680 |
|
|
nodes[-1]['v_out'] = None |
| 681 |
|
|
print 'closed at', pretty_vec(p) |
| 682 |
|
|
break |
| 683 |
|
|
elif (h > 0.99 * dpole and (cv > 0.9 or |
| 684 |
|
|
(cv > 0. and dpole * abs(sv) < ytol))): |
| 685 |
|
|
# slow down |
| 686 |
|
|
h = max(4.*err, dpole*cv * max(.9, 1-.1*dpole*cv)) |
| 687 |
|
|
|
| 688 |
|
|
if (pole_type == 'monopoles' and |
| 689 |
|
|
dpole < 0.01 and cv > .996): |
| 690 |
|
|
# approaching a monopole: end line with x**3 curve |
| 691 |
|
|
nodes[-1]['v_out'] = vnorm(v) * dpole |
| 692 |
|
|
v = vnorm(1.5 * vnorm(vpole) - |
| 693 |
|
|
.5 * vnorm(nodes[-1]['v_out'])) |
| 694 |
|
|
nodes.append({'p':nearest_pole[:2].copy(), |
| 695 |
|
|
'v_in':v * dpole, 'v_out':None}) |
| 696 |
|
|
l += h |
| 697 |
|
|
break |
| 698 |
|
|
|
| 699 |
|
|
if (pole_type == 'dipoles' and |
| 700 |
|
|
dpole < 0.01 and cv > .996): |
| 701 |
|
|
# approaching a dipole |
| 702 |
|
|
m = sign * vnorm(nearest_pole[2:4]) |
| 703 |
|
|
p = nodes[-1]['p'] + 2. * sc.dot(m, vpole) * m * dpole |
| 704 |
|
|
# approximation by a y=x**1.5 curve |
| 705 |
|
|
nodes[-1]['v_out'] = 2. * vnorm(v) * dpole |
| 706 |
|
|
nodes.append({'p':nearest_pole[:2].copy(), |
| 707 |
|
|
'v_in':sc.zeros(2), 'v_out':sc.zeros(2)}) |
| 708 |
|
|
l += h |
| 709 |
|
|
# check if the path is being closed |
| 710 |
|
|
v_end = self.first_point - p |
| 711 |
|
|
if ((dpole * abs(sinv(v, v_end)) < ytol) and |
| 712 |
|
|
(dpole * abs(cosv(v, v_end)) < xtol) and l > 1e-3): |
| 713 |
|
|
# path is closed |
| 714 |
|
|
nodes[-1]['v_out'] = None |
| 715 |
|
|
break |
| 716 |
|
|
if pass_dipoles == 0: |
| 717 |
|
|
nodes[-1]['v_out'] = None |
| 718 |
|
|
break |
| 719 |
|
|
if pass_dipoles > 0: |
| 720 |
|
|
pass_dipoles -= 1 |
| 721 |
|
|
v = f(p) |
| 722 |
|
|
nodes.append({'p':p.copy(), 'v_in':2.*vnorm(v)*dpole}) |
| 723 |
|
|
l += h |
| 724 |
|
|
continue |
| 725 |
|
|
|
| 726 |
|
|
# buckle detection at unknown places |
| 727 |
|
|
elif h < 0.01: |
| 728 |
|
|
# check change rate of curvature |
| 729 |
|
|
hh = h * 3. |
| 730 |
|
|
v0 = f(p + hh / 2. * v) |
| 731 |
|
|
v1 = f(p + hh * v) |
| 732 |
|
|
angle0 = atan2(v[1], v[0]) |
| 733 |
|
|
angle1 = atan2(v0[1], v0[0]) |
| 734 |
|
|
angle2 = atan2(v1[1], v1[0]) |
| 735 |
|
|
a0 = angle_dif(angle1, angle0) |
| 736 |
|
|
a1 = angle_dif(angle2, angle1) |
| 737 |
|
|
adif = angle_dif(a1, a0) |
| 738 |
|
|
corner_limit = 1e4 |
| 739 |
|
|
if abs(adif) / hh**2 > corner_limit: |
| 740 |
|
|
# assume a corner here |
| 741 |
|
|
if abs(a0) >= abs(a1): |
| 742 |
|
|
h0 = 0.; h1 = hh / 2. |
| 743 |
|
|
vm = vnorm(vnorm(v) + vnorm(v0)) |
| 744 |
|
|
else: |
| 745 |
|
|
h0 = hh / 2.; h1 = hh |
| 746 |
|
|
vm = vnorm(vnorm(v0) + vnorm(v1)) |
| 747 |
|
|
if vabs(vm)==0.: vm = vnorm(sc.array([v0[1], -v0[0]])) |
| 748 |
|
|
hc = op.brentq(lambda hc: sinv(f(p+hc*v), vm), h0, h1) |
| 749 |
|
|
v2 = f(p + hc / 2. * v) |
| 750 |
|
|
if sinv(f(p), vm) * sinv(f(p + 2.*hc*v2), vm) <= 0.: |
| 751 |
|
|
hc = op.brentq(lambda hc: sinv(f(p + hc * v2), |
| 752 |
|
|
vm), 0., 2. * hc) |
| 753 |
|
|
nodes[-1]['v_out'] = vnorm(nodes[-1]['v_in']) * hc |
| 754 |
|
|
# create a corner |
| 755 |
|
|
# use second-order formulas instead of runge-kutta |
| 756 |
|
|
p += hc * v2 |
| 757 |
|
|
print 'corner at', pretty_vec(p) |
| 758 |
|
|
v = vnorm(2. * v2 - v) |
| 759 |
|
|
nodes.append({'p':p.copy(),'v_in':v*hc,'corner':True}) |
| 760 |
|
|
l += h |
| 761 |
|
|
# check if the path is being closed |
| 762 |
|
|
v_end = self.first_point - p |
| 763 |
|
|
if ((dpole * abs(sinv(v, v_end)) < ytol) and |
| 764 |
|
|
(dpole * abs(cosv(v, v_end)) < xtol) and l > 1e-3): |
| 765 |
|
|
# path is closed |
| 766 |
|
|
nodes[-1]['v_out'] = None |
| 767 |
|
|
break |
| 768 |
|
|
# check area after the corner |
| 769 |
|
|
# lengths are chosen to ensure corner detection |
| 770 |
|
|
p0 = p + hh * .2 * f(p + hh * .2 * v1); va0 = f(p0) |
| 771 |
|
|
p1 = p0 + hh * .4 * va0; va1 = f(p1) |
| 772 |
|
|
p2 = p1 + hh * .4 * va1; va2 = f(p2) |
| 773 |
|
|
angle0 = atan2(va0[1], va0[0]) |
| 774 |
|
|
angle1 = atan2(va1[1], va1[0]) |
| 775 |
|
|
angle2 = atan2(va2[1], va2[0]) |
| 776 |
|
|
a0 = angle_dif(angle1, angle0) |
| 777 |
|
|
a1 = angle_dif(angle2, angle1) |
| 778 |
|
|
adif = angle_dif(a1, a0) |
| 779 |
|
|
if (abs(adif) / (.8*hh)**2 > corner_limit or |
| 780 |
|
|
abs(a0) + abs(a1) >= pi / 2.): |
| 781 |
|
|
print 'end edge at', pretty_vec(p) |
| 782 |
|
|
# direction after corner changes again -> end line |
| 783 |
|
|
nodes[-1]['v_out'] = None |
| 784 |
|
|
break |
| 785 |
|
|
vm = vnorm(1.25 * va1 - 0.25 * va2) |
| 786 |
|
|
v = f(p + hh * vm) |
| 787 |
|
|
nodes[-1]['v_out'] = vnorm(2. * vm - v) * hh |
| 788 |
|
|
p += vm * hh |
| 789 |
|
|
nodes.append({'p':p.copy(), 'v_in':v * hh}) |
| 790 |
|
|
l += h |
| 791 |
|
|
|
| 792 |
|
|
# make single and double runge-kutta step |
| 793 |
|
|
p11 = self.__rkstep(p, vnorm(v), f, h) |
| 794 |
|
|
p21 = self.__rkstep(p, vnorm(v), f, h / 2.) |
| 795 |
|
|
p22 = self.__rkstep(p21, f(p21), f, h / 2.) |
| 796 |
|
|
diff = vabs(p22 - p11) |
| 797 |
|
|
if diff < 2. * err: |
| 798 |
|
|
# accept step |
| 799 |
|
|
p = (16. * p22 - p11) / 15. |
| 800 |
|
|
nodes[-1]['v_out'] = vnorm(v) * h |
| 801 |
|
|
v = f(p) |
| 802 |
|
|
if vabs(v) == 0.: |
| 803 |
|
|
# field is zero, line is stuck -> end line |
| 804 |
|
|
nodes[-1]['v_out'] = None |
| 805 |
|
|
break |
| 806 |
|
|
if (len(nodes) >= 2 |
| 807 |
|
|
and vabs(nodes[-1]['p'] - nodes[-2]['p']) == 0.): |
| 808 |
|
|
if h > 2. * err: h /= 7. |
| 809 |
|
|
else: |
| 810 |
|
|
# point doesn_t move, line is stuck -> end line |
| 811 |
|
|
nodes = nodes[:-1] |
| 812 |
|
|
nodes[-1]['v_out'] = None |
| 813 |
|
|
break |
| 814 |
|
|
nodes.append({'p':p.copy(), 'v_in':v * h}) |
| 815 |
|
|
l += h |
| 816 |
|
|
|
| 817 |
|
|
# stop at the prohibited area |
| 818 |
|
|
if self.stop_funcs != None and self.stop_funcs != [None, None]: |
| 819 |
|
|
stop_fct = self.stop_funcs[{-1.0:0, 1.0:1}[sign]] |
| 820 |
|
|
if stop_fct(nodes[-1]['p']) > 0.0: |
| 821 |
|
|
while len(nodes) > 1 and stop_fct(nodes[-2]['p']) > 0.0: |
| 822 |
|
|
nodes = nodes[:-1] |
| 823 |
|
|
if len(nodes) > 1: |
| 824 |
|
|
p, p1 = nodes[-2]['p'], nodes[-1]['p'] |
| 825 |
|
|
t = op.brentq(lambda t: stop_fct(p + t * (p1 - p)), |
| 826 |
|
|
0.0, 1.0) |
| 827 |
|
|
nodes[-1]['p'] = p + t * (p1 - p) |
| 828 |
|
|
h = vabs(nodes[-1]['p'] - p) |
| 829 |
|
|
nodes[-2]['v_out'] = f(nodes[-2]['p']) * h |
| 830 |
|
|
nodes[-1]['v_in'] = f(nodes[-1]['p']) * h |
| 831 |
|
|
print 'stopped at', pretty_vec(nodes[-1]['p']) |
| 832 |
|
|
break |
| 833 |
|
|
|
| 834 |
|
|
# adapt step carefully |
| 835 |
|
|
if diff > 0.: |
| 836 |
|
|
factor = (err / diff) ** .25 |
| 837 |
|
|
if h < h_old: h_new = min((h + h_old) / 2., h * factor) |
| 838 |
|
|
else: h_new = h * max(0.5, factor) |
| 839 |
|
|
h_old = h |
| 840 |
|
|
h = h_new |
| 841 |
|
|
else: |
| 842 |
|
|
h_old = h |
| 843 |
|
|
h *= 2. |
| 844 |
|
|
h = min(hmax, max(err, h)) |
| 845 |
|
|
|
| 846 |
|
|
nodes[-1]['v_out'] = None |
| 847 |
|
|
if i == maxn: |
| 848 |
|
|
print maxn, 'integration steps exceeded at', pretty_vec(p) |
| 849 |
|
|
if l >= maxr: |
| 850 |
|
|
print 'integration boundary',str(maxr),'exceeded at',pretty_vec(p) |
| 851 |
|
|
return nodes |
| 852 |
|
|
|
| 853 |
|
|
def __is_loop(self, nodes, path_close_tol): |
| 854 |
|
|
if vabs(nodes[0]['p'] - nodes[-1]['p']) > max(5e-4, path_close_tol): |
| 855 |
|
|
return False |
| 856 |
|
|
length = 0. |
| 857 |
|
|
for i in range(1, len(nodes)): |
| 858 |
|
|
length += vabs(nodes[i]['p'] - nodes[i-1]['p']) |
| 859 |
|
|
if length > 5e-3: |
| 860 |
|
|
return True |
| 861 |
|
|
return False |
| 862 |
|
|
|
| 863 |
|
|
def __create_nodes(self, directions, |
| 864 |
|
|
maxn, maxr, hmax, pass_dipoles, path_close_tol): |
| 865 |
|
|
''' |
| 866 |
|
|
creates self.nodes from one or two parts |
| 867 |
|
|
wrapper for __self.create_nodes_part |
| 868 |
|
|
''' |
| 869 |
|
|
closed = False |
| 870 |
|
|
if (directions == 'forward'): |
| 871 |
|
|
self.nodes = self.__create_nodes_part( |
| 872 |
|
|
1., maxn, maxr, hmax, pass_dipoles, path_close_tol) |
| 873 |
|
|
else: |
| 874 |
|
|
nodes1 = self.__create_nodes_part( |
| 875 |
|
|
-1., maxn, maxr, hmax, pass_dipoles, path_close_tol) |
| 876 |
|
|
# reverse nodes1 |
| 877 |
|
|
nodes1.reverse() |
| 878 |
|
|
for node in nodes1: |
| 879 |
|
|
v_out = node['v_out'] |
| 880 |
|
|
if node['v_in'] == None: node['v_out'] = None |
| 881 |
|
|
else: node['v_out'] = -node['v_in'] |
| 882 |
|
|
if v_out == None: node['v_in'] = None |
| 883 |
|
|
else: node['v_in'] = -v_out |
| 884 |
|
|
self.nodes = nodes1 |
| 885 |
|
|
if len(self.nodes) > 0: self.first_point = self.nodes[0]['p'] |
| 886 |
|
|
if directions != 'backward': |
| 887 |
|
|
# is it already a closed loop? |
| 888 |
|
|
if not self.__is_loop(self.nodes, path_close_tol): |
| 889 |
|
|
nodes2 = self.__create_nodes_part( |
| 890 |
|
|
1., maxn, maxr, hmax, pass_dipoles, path_close_tol) |
| 891 |
|
|
self.nodes[-1]['v_out'] = nodes2[0]['v_out'] |
| 892 |
|
|
self.nodes += nodes2[1:] |
| 893 |
|
|
|
| 894 |
|
|
# append accumulated normalized sum |
| 895 |
|
|
self.nodes[0]['t'] = 0. |
| 896 |
|
|
for i in range(1, len(self.nodes)): |
| 897 |
|
|
self.nodes[i]['t'] = (self.nodes[i-1]['t'] |
| 898 |
|
|
+ vabs(self.nodes[i-1]['p'] - self.nodes[i]['p'])) |
| 899 |
|
|
length = self.nodes[-1]['t'] |
| 900 |
|
|
if length != 0.: |
| 901 |
|
|
for i in range(1, len(self.nodes)): |
| 902 |
|
|
self.nodes[i]['t'] /= length |
| 903 |
|
|
# add corner tag to all nodes |
| 904 |
|
|
for i, node in enumerate(self.nodes): |
| 905 |
|
|
if not node.has_key('corner'): |
| 906 |
|
|
self.nodes[i]['corner'] = False |
| 907 |
|
|
|
| 908 |
|
|
def get_position(self, t): |
| 909 |
|
|
''' |
| 910 |
|
|
dense output routine |
| 911 |
|
|
t: parameter, 0 <= t <= 1 |
| 912 |
|
|
''' |
| 913 |
|
|
nodes = self.nodes |
| 914 |
|
|
if len(nodes) == 1: |
| 915 |
|
|
return nodes[0]['p'] |
| 916 |
|
|
if len(nodes) <= 0: |
| 917 |
|
|
return sc.zeros(2) |
| 918 |
|
|
if t != 1.: t = t % 1. |
| 919 |
|
|
n, p = list_interpolate([i['t'] for i in nodes], t) |
| 920 |
|
|
p0, v0 = nodes[n]['p'], nodes[n]['v_out'] |
| 921 |
|
|
p1, v1 = nodes[n+1]['p'], nodes[n+1]['v_in'] |
| 922 |
|
|
# cubic bezier interpolation (hermite interpolation) |
| 923 |
|
|
q = 1. - p |
| 924 |
|
|
xy = q*p0 + p*p1 + p * q * ((p - q) * (p1 - p0) + (q*v0 - p*v1)) |
| 925 |
|
|
return xy |
| 926 |
|
|
|
| 927 |
|
|
def __bending(self, p0, p3, t0, t3): |
| 928 |
|
|
# calculate two extra points on intervall |
| 929 |
|
|
p1 = self.get_position((2.*t0 + t3) / 3.) |
| 930 |
|
|
p2 = self.get_position((t0 + 2.*t3) / 3.) |
| 931 |
|
|
# d1, d2: point distances from straight line |
| 932 |
|
|
d1 = (p1 - p0)[0] * (p3 - p0)[1] - (p1 - p0)[1] * (p3 - p0)[0] |
| 933 |
|
|
d1 /= vabs(p3 - p0) |
| 934 |
|
|
d2 = (p2 - p0)[0] * (p3 - p0)[1] - (p2 - p0)[1] * (p3 - p0)[0] |
| 935 |
|
|
d2 /= vabs(p3 - p0) |
| 936 |
|
|
dsum, ddif = d1 + d2, d1 - d2 |
| 937 |
|
|
d = 0. |
| 938 |
|
|
if abs(ddif) < 1e-5: |
| 939 |
|
|
d = 10. / 9. * (abs(d1) + abs(d2)) / 2. |
| 940 |
|
|
else: |
| 941 |
|
|
# calculate line bending as max distance of a deg-3 polynomial: |
| 942 |
|
|
y = lambda x: 13.5 * x * (1.-x) * (d1 * (2./3.-x) + d2 * (x-1./3.)) |
| 943 |
|
|
# all the factors come from the quadratic formula |
| 944 |
|
|
xm = .5 + dsum / (18. * ddif) |
| 945 |
|
|
xd = sqrt(27. * ddif**2 + dsum**2) / (18. * ddif) |
| 946 |
|
|
x1, x2 = min(xm + xd, xm - xd), max(xm + xd, xm - xd) |
| 947 |
|
|
if x1 > 0.: |
| 948 |
|
|
d = max(d, abs(y(x1))) |
| 949 |
|
|
if x2 < 1.: |
| 950 |
|
|
d = max(d, abs(y(x2))) |
| 951 |
|
|
return d |
| 952 |
|
|
|
| 953 |
|
|
def __get_polyline(self, t0, t1, digits=3.5, maxdist=10., mindist=4e-4): |
| 954 |
|
|
''' |
| 955 |
|
|
returns points of an adapted polyline, |
| 956 |
|
|
representing the fieldline to an accuracy of digits. |
| 957 |
|
|
no corner should be between t0 and t1. |
| 958 |
|
|
''' |
| 959 |
|
|
f = self.get_position |
| 960 |
|
|
t_list = sc.linspace(t0, t1, 10) |
| 961 |
|
|
value_list = [f(t) for t in t_list] |
| 962 |
|
|
|
| 963 |
|
|
# adapt t_list |
| 964 |
|
|
num = 0; num_success = 0 |
| 965 |
|
|
while len(t_list) > 2: |
| 966 |
|
|
ratios = []; delta_t = [] |
| 967 |
|
|
N_old = len(t_list) - 1 |
| 968 |
|
|
success = True |
| 969 |
|
|
# get bending |
| 970 |
|
|
for i in range(N_old): |
| 971 |
|
|
bend = self.__bending(value_list[i], value_list[i+1], |
| 972 |
|
|
t_list[i], t_list[i + 1]) |
| 973 |
|
|
d = vabs(value_list[i+1] - value_list[i]) |
| 974 |
|
|
# keep point distance smaller than maxdist |
| 975 |
|
|
ratio = d / maxdist |
| 976 |
|
|
if num > 10: exponent = 1. / (num - 8.) |
| 977 |
|
|
else: exponent = 0.5 |
| 978 |
|
|
# find best ratio, assuming bending is proportional to d**2 |
| 979 |
|
|
if bend != 0.: |
| 980 |
|
|
ratio = max(ratio, (bend / 0.1 ** digits)**exponent) |
| 981 |
|
|
ratio = min(ratio, d / mindist) |
| 982 |
|
|
if ratio > 1.1: # 1 + 0.1 for termination safety |
| 983 |
|
|
success = False |
| 984 |
|
|
ratio = min(max(.25, ratio), 4.) # prevent too big changes |
| 985 |
|
|
ratios.append(ratio) |
| 986 |
|
|
delta_t.append(t_list[i + 1] - t_list[i]) |
| 987 |
|
|
|
| 988 |
|
|
n = sum(ratios) |
| 989 |
|
|
N = max(1, ceil(n)) # new intervall number must be an integer |
| 990 |
|
|
num += 1 |
| 991 |
|
|
# check if we all intervalls are good enough and we are finished |
| 992 |
|
|
if success == True: num_success += 1 |
| 993 |
|
|
else: num_success = 0 |
| 994 |
|
|
if num_success > 2 and N < N_old: num_success = 2 |
| 995 |
|
|
if num_success >= 3: break |
| 996 |
|
|
if num >= 25: |
| 997 |
|
|
print 'polyline creation did not converge after', num, 'tries!' |
| 998 |
|
|
break |
| 999 |
|
|
ratios = [ratio * N / n for ratio in ratios] |
| 1000 |
|
|
|
| 1001 |
|
|
# rearrange t_list |
| 1002 |
|
|
t_list = [t0] # initialize again |
| 1003 |
|
|
N0 = 0; Nt = 0.; N1 = 0.; t = t0 |
| 1004 |
|
|
for i in range(N_old): |
| 1005 |
|
|
N1 += ratios[i] |
| 1006 |
|
|
while N1 - N0 >= 1.: |
| 1007 |
|
|
N0 += 1 |
| 1008 |
|
|
t += delta_t[i] * (N0 - Nt) / ratios[i] |
| 1009 |
|
|
Nt = N0 |
| 1010 |
|
|
if len(t_list) == N: |
| 1011 |
|
|
break |
| 1012 |
|
|
t_list.append(t) |
| 1013 |
|
|
t += delta_t[i] * (N1 - Nt) / ratios[i] |
| 1014 |
|
|
Nt = N1 |
| 1015 |
|
|
t_list.append(t1) |
| 1016 |
|
|
value_list = [f(t) for t in t_list] |
| 1017 |
|
|
return value_list, t_list |
| 1018 |
|
|
|
| 1019 |
|
|
def __out_of_bounds(self, p, bounds): |
| 1020 |
|
|
''' |
| 1021 |
|
|
returns a points distance to the drawing area |
| 1022 |
|
|
>0: outside; <=0: inside |
| 1023 |
|
|
''' |
| 1024 |
|
|
if self.bounds_func != None: |
| 1025 |
|
|
s = self.bounds_func(p) |
| 1026 |
|
|
if s > 0.: return s |
| 1027 |
|
|
if bounds == None: return -1. |
| 1028 |
|
|
if (p[0] < bounds['x0'] or p[1] < bounds['y0'] |
| 1029 |
|
|
or p[0] > bounds['x1'] or p[1] > bounds['y1']): |
| 1030 |
|
|
return sqrt((p[0] - bounds['x0'])**2 + (p[1] - bounds['y0'])**2 |
| 1031 |
|
|
+ (bounds['x1'] - p[0])**2 + (bounds['y1'] - p[1])**2) |
| 1032 |
|
|
else: |
| 1033 |
|
|
return max(bounds['x0'] - p[0], bounds['y0'] - p[1], |
| 1034 |
|
|
p[0] - bounds['x1'], p[1] - bounds['y1']) |
| 1035 |
|
|
|
| 1036 |
|
|
def get_polylines(self, digits=3.5, maxdist=10., bounds=None): |
| 1037 |
|
|
''' |
| 1038 |
|
|
returns polyline segments that are inside of bounds. |
| 1039 |
|
|
the path is represented as a set of adapted line segments |
| 1040 |
|
|
which are cut at the image bounds and at edges. |
| 1041 |
|
|
''' |
| 1042 |
|
|
if len(self.nodes) <= 1: return [] |
| 1043 |
|
|
|
| 1044 |
|
|
# search for all corners |
| 1045 |
|
|
corners = [] |
| 1046 |
|
|
for node in self.nodes: |
| 1047 |
|
|
if node['corner']: corners.append(node['t']) |
| 1048 |
|
|
if len(corners) == 0 or corners[0] != 0.: corners.insert(0, 0.) |
| 1049 |
|
|
if corners[-1] != 1.: corners.append(1.) |
| 1050 |
|
|
|
| 1051 |
|
|
# search for points where line intersects bounds |
| 1052 |
|
|
edges = []; parts_outside = False; inside1 = False; t1 = 0. |
| 1053 |
|
|
if self.__out_of_bounds(self.nodes[0]['p'], bounds) <= 0.: |
| 1054 |
|
|
inside1 = True |
| 1055 |
|
|
edges.append({'t0':0.}) |
| 1056 |
|
|
for i in range(1, len(self.nodes)): |
| 1057 |
|
|
t0 = t1; t1 = self.nodes[i]['t'] |
| 1058 |
|
|
p1 = self.nodes[i]['p'] |
| 1059 |
|
|
inside0 = inside1 |
| 1060 |
|
|
inside1 = (self.__out_of_bounds(p1, bounds) <= 0.) |
| 1061 |
|
|
if inside1: |
| 1062 |
|
|
if not inside0: |
| 1063 |
|
|
edges.append({'t0':op.brentq(lambda t: |
| 1064 |
|
|
self.__out_of_bounds(self.get_position(t), |
| 1065 |
|
|
bounds), t0, t1)}) |
| 1066 |
|
|
if i == len(self.nodes) - 1: |
| 1067 |
|
|
edges[-1]['t1'] = 1. |
| 1068 |
|
|
else: |
| 1069 |
|
|
parts_outside = True |
| 1070 |
|
|
if inside0: |
| 1071 |
|
|
edges[-1]['t1'] = (op.brentq(lambda t: |
| 1072 |
|
|
self.__out_of_bounds(self.get_position(t), |
| 1073 |
|
|
bounds), t0, t1)) |
| 1074 |
|
|
|
| 1075 |
|
|
# all points are outside the drawing area |
| 1076 |
|
|
if len(edges) == 0: return [] |
| 1077 |
|
|
|
| 1078 |
|
|
# join first and last segment |
| 1079 |
|
|
if (len(edges) > 1 and |
| 1080 |
|
|
edges[0]['t0'] == 0. and edges[-1]['t1'] == 1. and |
| 1081 |
|
|
vabs(self.get_position(1.) - self.get_position(0.)) <= 1e-5): |
| 1082 |
|
|
edges[0]['t0'] = edges[-1]['t0'] - 1. |
| 1083 |
|
|
edges = edges[:-1] |
| 1084 |
|
|
|
| 1085 |
|
|
# insert corners to all segments |
| 1086 |
|
|
for edge in edges: |
| 1087 |
|
|
edge['corners'] = [] |
| 1088 |
|
|
for c in corners: |
| 1089 |
|
|
if edge['t0'] < c and c < edge['t1']: |
| 1090 |
|
|
edge['corners'].append(c) |
| 1091 |
|
|
|
| 1092 |
|
|
# create final polylines |
| 1093 |
|
|
polyline = [] |
| 1094 |
|
|
for interval in edges: |
| 1095 |
|
|
line = [] |
| 1096 |
|
|
t_list = [interval['t0']] + interval['corners'] + [interval['t1']] |
| 1097 |
|
|
for i in range(1, len(t_list)): |
| 1098 |
|
|
pl = self.__get_polyline(t_list[i-1], t_list[i], |
| 1099 |
|
|
digits, maxdist)[0] |
| 1100 |
|
|
if i == 1: line += pl |
| 1101 |
|
|
else: line += pl[1:] |
| 1102 |
|
|
if len(line) >= 2: |
| 1103 |
|
|
polyline.append({'path':line, |
| 1104 |
|
|
'start':(interval['t0']==0.), 'end':(interval['t1']==1.)}) |
| 1105 |
|
|
return polyline |
| 1106 |
|
|
|
| 1107 |
|
|
|
| 1108 |
|
|
|
| 1109 |
|
|
class Field: |
| 1110 |
|
|
''' |
| 1111 |
|
|
represents an electromagnetic field together with |
| 1112 |
|
|
charges, potential, setup etc. |
| 1113 |
|
|
''' |
| 1114 |
|
|
def __init__ (self, elements={}): |
| 1115 |
|
|
self.elements = {} |
| 1116 |
|
|
for name, params in elements.iteritems(): |
| 1117 |
|
|
self.add_element(name, params) |
| 1118 |
|
|
|
| 1119 |
|
|
''' |
| 1120 |
|
|
possible types: |
| 1121 |
|
|
'homogeneous': [Fx, Fy] |
| 1122 |
|
|
'gradient': [g, Ax, Ay, Bx, By] |
| 1123 |
|
|
'monopoles': [x, y, charge] |
| 1124 |
|
|
'dipoles': [x, y, phi, q] |
| 1125 |
|
|
'quadrupoles': [x, y, phi, q] |
| 1126 |
|
|
'wires': [x, y, I] |
| 1127 |
|
|
'charged_planes': [x0, y0, x1, y1, charge] |
| 1128 |
|
|
'ringcurrents': [x0, y0, phi, R, I] |
| 1129 |
|
|
'coils': [x0, y0, phi, R, Lhalf, I] |
| 1130 |
|
|
'custom': user defined function |
| 1131 |
|
|
''' |
| 1132 |
|
|
|
| 1133 |
|
|
def add_element(self, name, params): |
| 1134 |
|
|
if len(params) >= 1 and str(type(params[0])) == "<type 'function'>": |
| 1135 |
|
|
if name in self.elements: self.elements[name] += params |
| 1136 |
|
|
else: self.elements[name] = params |
| 1137 |
|
|
else: |
| 1138 |
|
|
el = [sc.array(param, dtype='float') for param in params] |
| 1139 |
|
|
if name in self.elements: self.elements[name] += el |
| 1140 |
|
|
else: self.elements[name] = el |
| 1141 |
|
|
|
| 1142 |
|
|
def get_elements(self, name): |
| 1143 |
|
|
if name in self.elements: return self.elements[name] |
| 1144 |
|
|
else: return [] |
| 1145 |
|
|
|
| 1146 |
|
|
def F(self, xy): |
| 1147 |
|
|
''' |
| 1148 |
|
|
returns the field force as a vector |
| 1149 |
|
|
''' |
| 1150 |
|
|
Fxy = sc.zeros(2) |
| 1151 |
|
|
|
| 1152 |
|
|
# homogeneous: homogeneus field in a given direction |
| 1153 |
|
|
for hom in self.get_elements('homogeneous'): |
| 1154 |
|
|
Fxy += hom |
| 1155 |
|
|
|
| 1156 |
|
|
# gradient: gradient field in a given direction |
| 1157 |
|
|
for gra in self.get_elements('gradient'): |
| 1158 |
|
|
g = gra[0] |
| 1159 |
|
|
A = sc.array(gra[1:3]) |
| 1160 |
|
|
B = sc.array(gra[3:5]) |
| 1161 |
|
|
cpsi = sc.dot( A, B ) |
| 1162 |
|
|
r = xy |
| 1163 |
|
|
Fxy[0] += g*(2.0*(A[0]*B[0]-cpsi/2.0)*r[0] + (A[0]*B[1]+A[1]*B[0])*r[1])/2.0 |
| 1164 |
|
|
Fxy[1] += g*((A[0]*B[1]+A[1]*B[0])*r[0] + 2.0*(A[1]*B[1]-cpsi/2.0)*r[1])/2.0 |
| 1165 |
|
|
|
| 1166 |
|
|
# monopoles: electric charges and magnetic monopoles |
| 1167 |
|
|
for mon in self.get_elements('monopoles'): |
| 1168 |
|
|
r = xy - sc.array(mon[:2]) |
| 1169 |
|
|
d = vabs(r) |
| 1170 |
|
|
if d != 0.: |
| 1171 |
|
|
Fxy += mon[-1] * r / d**3 |
| 1172 |
|
|
|
| 1173 |
|
|
# dipoles: pointlike electric or magnetic dipole |
| 1174 |
|
|
for dip in self.get_elements('dipoles'): |
| 1175 |
|
|
r = xy - sc.array(dip[:2]) |
| 1176 |
|
|
d = vabs(r) |
| 1177 |
|
|
if d != 0.: |
| 1178 |
|
|
p = sc.array(dip[2:4]) |
| 1179 |
|
|
Fxy += (3. * sc.dot(p, r) * r - sc.dot(r, r) * p) / (4.*pi*d**5) |
| 1180 |
|
|
else: |
| 1181 |
|
|
# the sign of this is unphysical, but more useful |
| 1182 |
|
|
return dip[2:4] |
| 1183 |
|
|
|
| 1184 |
|
|
# quadrupoles: pointlike electric or magnetic quadrupoles |
| 1185 |
|
|
for quad in self.get_elements('quadrupoles'): |
| 1186 |
|
|
r = xy - sc.array(quad[:2]) |
| 1187 |
|
|
d = vabs(r) |
| 1188 |
|
|
r /= d |
| 1189 |
|
|
if d != 0.: |
| 1190 |
|
|
p = rot([0,1], quad[2]) |
| 1191 |
|
|
pr = sc.dot(p, r) |
| 1192 |
|
|
Fxy += (((5.*pr**2 - 1.) * r - 2.*pr * p) * |
| 1193 |
|
|
3.*quad[3] / (4.*pi * d**4)) |
| 1194 |
|
|
else: |
| 1195 |
|
|
return quad[2:4] |
| 1196 |
|
|
|
| 1197 |
|
|
# wires: infinite straight wire perpendicular to image plane |
| 1198 |
|
|
for wire in self.get_elements('wires'): |
| 1199 |
|
|
r = xy - sc.array(wire[:2]) |
| 1200 |
|
|
Fxy += wire[-1]/(2*pi) * sc.array([-r[1], r[0]]) / (r[0]**2 + r[1]**2) |
| 1201 |
|
|
|
| 1202 |
|
|
# charged_planes: rectangular plane with edges [x0,y0] and [x1,y1] |
| 1203 |
|
|
# perpendicular to image plane and infinite in z-direction |
| 1204 |
|
|
for plane in self.get_elements('charged_planes'): |
| 1205 |
|
|
r = xy - .5 * (plane[2:4] + plane[0:2]) |
| 1206 |
|
|
p = .5 * (plane[2:4] - plane[0:2]) |
| 1207 |
|
|
rsq = r[0]**2 + r[1]**2 |
| 1208 |
|
|
psq = p[0]**2 + p[1]**2 |
| 1209 |
|
|
X = r[0] * p[1] - p[0] * r[1] |
| 1210 |
|
|
Y = r[1] * p[1] + p[0] * r[0] |
| 1211 |
|
|
if X == 0.: Fa = 0. |
| 1212 |
|
|
else: Fa = atan((psq - Y) / X) + atan((psq + Y) / X) |
| 1213 |
|
|
Fb = atanh(2. * Y / (rsq + psq)) |
| 1214 |
|
|
Fxy += (plane[4] / (4*pi*psq) |
| 1215 |
|
|
* sc.array([Fa*p[1]+Fb*p[0], Fb*p[1]-Fa*p[0]])) |
| 1216 |
|
|
|
| 1217 |
|
|
# ringcurrents: round currentloop perpendicular to image plane |
| 1218 |
|
|
# caution: slow because of numerical integration |
| 1219 |
|
|
for ring in self.get_elements('ringcurrents'): |
| 1220 |
|
|
r = xy - sc.array(ring[:2]) |
| 1221 |
|
|
# change into a relative coordinate system with a and b |
| 1222 |
|
|
a, b = rot(r, -ring[2]) / ring[3] |
| 1223 |
|
|
c = 1. + a**2 + b**2 |
| 1224 |
|
|
def fa(t): h=cos(t); return (1.-b*h) / sqrt((c-2.*b*h)**3) |
| 1225 |
|
|
def fb(t): h=cos(t); return (a * h) / sqrt((c-2.*b*h)**3) |
| 1226 |
|
|
Fa = ig.quad(fa, 0., pi, epsrel=0., full_output=True)[0] |
| 1227 |
|
|
Fb = ig.quad(fb, 0., pi, epsrel=0., full_output=True)[0] |
| 1228 |
|
|
# backtransform |
| 1229 |
|
|
Fxy += rot([Fa, Fb], ring[2]) * (ring[-1] / ring[3]) |
| 1230 |
|
|
|
| 1231 |
|
|
# coil: dense cylinder coil or cylinder magnet |
| 1232 |
|
|
# caution: slow because of numerical integration |
| 1233 |
|
|
for coil in self.get_elements('coils'): |
| 1234 |
|
|
r = xy - sc.array(coil[:2]) |
| 1235 |
|
|
# change into a relative coordinate system with a and b |
| 1236 |
|
|
a, b = rot(r, -coil[2]) / coil[3] |
| 1237 |
|
|
c = 1. + b**2 |
| 1238 |
|
|
l = coil[4] / coil[3] |
| 1239 |
|
|
am = a - l; am2 = am ** 2 |
| 1240 |
|
|
ap = a + l; ap2 = ap ** 2 |
| 1241 |
|
|
def fa(t): |
| 1242 |
|
|
h = cos(t); |
| 1243 |
|
|
d = c - 2. * b * h |
| 1244 |
|
|
return (1. - b*h) / d * (ap / sqrt(d + ap2) |
| 1245 |
|
|
- am / sqrt(d + am2)) |
| 1246 |
|
|
def fb(t): |
| 1247 |
|
|
h = cos(t); |
| 1248 |
|
|
d = c - 2. * b * h |
| 1249 |
|
|
return h / sqrt(d + am2) - h / sqrt(d + ap2) |
| 1250 |
|
|
Fa = ig.quad(fa, 0., pi, full_output=True)[0] |
| 1251 |
|
|
Fb = ig.quad(fb, 0., pi, full_output=True)[0] |
| 1252 |
|
|
# backtransform |
| 1253 |
|
|
Fxy += rot([Fa, Fb], coil[2]) * (coil[-1] / (2. * coil[4])) |
| 1254 |
|
|
|
| 1255 |
|
|
# custom: user defined function |
| 1256 |
|
|
for cust in self.get_elements('custom'): |
| 1257 |
|
|
Fxy += cust(xy) |
| 1258 |
|
|
|
| 1259 |
|
|
return Fxy |
| 1260 |
|
|
|
| 1261 |
|
|
def Fn(self, xy): |
| 1262 |
|
|
''' |
| 1263 |
|
|
returns the normalized field force, i.e. direction of field lines |
| 1264 |
|
|
''' |
| 1265 |
|
|
force = self.F(xy) |
| 1266 |
|
|
d = vabs(force) |
| 1267 |
|
|
if (d != 0.): return force / d |
| 1268 |
|
|
return sc.array([0., 0.]) |
| 1269 |
|
|
|
| 1270 |
|
|
### append your specific field creation here ### |
| 1271 |
|
|
# see http://commons.wikimedia.org/wiki/File:VFPt_charges_plus_minus.svg for an example |
