'''

args = list(args)

fn = ([logging.debug,logging.info,logging.warning,logging.error])[args.pop(0)-1]

msg = ' '.join(map(str, args))

'''

'''

d = vabs(x)

if d != 0.: return sc.array(x) / vabs(x)

'''

s = sin(phi); c = cos(phi)

'''

d1 = sum(v1**2); d2 = sum(v2**2)

if d1 != 1.: d1 = sqrt(d1)

if d2 != 1.: d2 = sqrt(d2)

if d1 * d2 == 0.: return 1.

'''

d1 = sum(v1**2); d2 = sum(v2**2)

if d1 != 1.: d1 = sqrt(d1)

if d2 != 1.: d2 = sqrt(d2)

if d1 * d2 == 0.: return 0.

n = max(0, bisect.bisect_right(l, t) - 1)

s = None

if t < l[0]:

if l[1] - l[0] == 0.:

s = 0.

else:

s = (t - l[0]) / float(l[1] - l[0])

elif t >= l[-1]:

n = len(l) - 2

if l[-1] - l[-2] == 0.:

s = 1.

else:

s = (t - l[-2]) / float(l[-1] - l[-2])

else:

s = (t - l[n]) / (l[n+1] - l[n])

'''

def __init__ (self, name, width=800, height=600, digits=3.5, unit=100,

center=None, licence='GFDL-cc', commons=False):

self.name = name

self.width = float(width)

self.height = float(height)

self.digits = float(digits)

self.unit = float(unit)

self.licence = licence

self.commons = commons

if center == None: self.center = [width / 2., height / 2.]

else: self.center = [float(i) for i in center]

# create document structure

self.svg = etree.Element('svg',

nsmap={None: 'http://www.w3.org/2000/svg',

'xlink': 'http://www.w3.org/1999/xlink'})

self.svg.set('version', '1.1')

self.svg.set('baseProfile', 'full')

self.svg.set('width', str(int(width)))

self.svg.set('height', str(int(height)))

# title

self.title = etree.SubElement(self.svg, 'title')

self.title.text = self.name

# description

self.desc = etree.SubElement(self.svg, 'desc')

self.desc.text = ''

self.desc.text += self.name + '\n'

self.desc.text += 'created with VectorFieldPlot ' + version + '\n'

self.desc.text += 'http://commons.wikimedia.org/wiki/User:Geek3/VectorFieldPlot\n'

if commons:

self.desc.text += """

""".format(self.name)

if self.licence == 'GFDL-cc':

self.desc.text += """rights: GNU Free Documentation license,

self.desc.text += ' '

# background

self.background = etree.SubElement(self.svg, 'rect')

self.background.set('id', 'background')

self.background.set('x', '0')

self.background.set('y', '0')

self.background.set('width', str(width))

self.background.set('height', str(height))

self.background.set('fill', '#ffffff')

# image elements

self.content = etree.SubElement(self.svg, 'g')

self.content.set('id', 'image')

self.content.set('transform',

'translate({0},{1}) scale({2},-{2})'.format(

self.center[0], self.center[1], self.unit))

self.arrow_geo = {'x_nock':0.3,'x_head':3.8,'x_tail':-2.2,'width':4.5}

def __get_defs(self):

if 'defs' not in dir(self):

self.defs = etree.Element('defs')

self.desc.addnext(self.defs)

return self.defs

def __check_fieldlines(self, linecolor='#000000', linewidth=1.):

if 'fieldlines' not in dir(self):

self.fieldlines = etree.SubElement(self.content, 'g')

self.fieldlines.set('id', 'fieldlines')

self.fieldlines.set('fill', 'none')

self.fieldlines.set('stroke', linecolor)

self.fieldlines.set('stroke-width',

str(linewidth / self.unit))

self.fieldlines.set('stroke-linejoin', 'round')

self.fieldlines.set('stroke-linecap', 'round')

if 'count_fieldlines' not in dir(self): self.count_fieldlines = 0

def __check_symbols(self):

if 'symbols' not in dir(self):

self.symbols = etree.SubElement(self.content, 'g')

self.symbols.set('id', 'symbols')

if 'count_symbols' not in dir(self): self.count_symbols = 0

def __check_whitespot(self):

if 'whitespot' not in dir(self):

self.whitespot = etree.SubElement(

self.__get_defs(), 'radialGradient')

self.whitespot.set('id', 'white_spot')

for attr, val in [['cx', '0.65'], ['cy', '0.7'], ['r', '0.75']]:

self.whitespot.set(attr, val)

for col, of, op in [['#ffffff', '0', '0.7'],

['#ffffff', '0.1', '0.5'], ['#ffffff', '0.6', '0'],

['#000000', '0.6', '0'], ['#000000', '0.75', '0.05'],

['#000000', '0.85', '0.15'], ['#000000', '1', '0.5']]:

stop = etree.SubElement(self.whitespot, 'stop')

stop.set('stop-color', col)

stop.set('offset', of)

stop.set('stop-opacity', op)

def __get_arrowname(self, fillcolor='#000000'):

if 'arrows' not in dir(self):

self.arrows = {}

if fillcolor not in self.arrows.keys():

arrow = etree.SubElement(self.__get_defs(), 'path')

self.arrows[fillcolor] = arrow

arrow.set('id', 'arrow' + str(len(self.arrows)))

arrow.set('stroke', 'none')

arrow.set('fill', fillcolor)

arrow.set('transform', 'scale({0})'.format(1. / self.unit))

arrow.set('d',

'M {0},0 L {1},{3} L {2},0 L {1},-{3} L {0},0 Z'.format(

self.arrow_geo['x_nock'], self.arrow_geo['x_tail'],

self.arrow_geo['x_head'], self.arrow_geo['width'] / 2.))

return self.arrows[fillcolor].get('id')

def draw_charges(self, field, scale=1.):

if 'monopoles' not in field.elements: return

charges = field.elements['monopoles']

self.__check_symbols()

self.__check_whitespot()

for charge in charges:

c_group = etree.SubElement(self.symbols, 'g')

self.count_symbols += 1

c_group.set('id', 'charge{0}'.format(self.count_symbols))

c_group.set('transform',

'translate({0},{1}) scale({2},{2})'.format(

charge[0], charge[1], float(scale) / self.unit))

#### charge drawing ####

c_bg = etree.SubElement(c_group, 'circle')

c_shade = etree.SubElement(c_group, 'circle')

c_symb = etree.SubElement(c_group, 'path')

if charge[2] >= 0.: c_bg.set('style', 'fill:#ff0000; stroke:none')

else: c_bg.set('style', 'fill:#0000ff; stroke:none')

for attr, val in [['cx', '0'], ['cy', '0'], ['r', '14']]:

c_bg.set(attr, val)

c_shade.set(attr, val)

c_shade.set('style',

'fill:url(#white_spot); stroke:#000000; stroke-width:2')

# plus sign

if charge[2] >= 0.:

c_symb.set('d', 'M 2,2 V 8 H -2 V 2 H -8 V -2'

+ ' H -2 V -8 H 2 V -2 H 8 V 2 H 2 Z')

# minus sign

else: c_symb.set('d', 'M 8,2 H -8 V -2 H 8 V 2 Z')

c_symb.set('style', 'fill:#000000; stroke:none')

def draw_currents(self, field, scale=1.):

if ('wires' not in field.elements

and 'ringcurrents' not in field.elements):

return

self.__check_symbols()

self.__check_whitespot()

currents = []

if 'wires' in field.elements:

for i in field.elements['wires']:

currents.append(i)

if 'ringcurrents' in field.elements:

for i in field.elements['ringcurrents']:

currents.append(list(i[:2] + rot([0., i[3]], i[2])) + [i[-1]])

currents.append(list(i[:2] - rot([0., i[3]], i[2])) + [-i[-1]])

for cur in currents:

c_group = etree.SubElement(self.symbols, 'g')

self.count_symbols += 1

if cur[-1] >= 0.: direction = 'out'

else: direction = 'in'

c_group.set('id',

'current_{0}{1}'.format(direction, self.count_symbols))

c_group.set('transform',

'translate({0},{1}) scale({2},{2})'.format(

cur[0], cur[1], float(scale) / self.unit))

#### current drawing ####

c_bg = etree.SubElement(c_group, 'circle')

c_shade = etree.SubElement(c_group, 'circle')

c_bg.set('style', 'fill:#b0b0b0; stroke:none')

for attr, val in [['cx', '0'], ['cy', '0'], ['r', '14']]:

c_bg.set(attr, val)

c_shade.set(attr, val)

c_shade.set('style',

'fill:url(#white_spot); stroke:#000000; stroke-width:2')

if cur[-1] >= 0.: # dot

c_symb = etree.SubElement(c_group, 'circle')

c_symb.set('cx', '0')

c_symb.set('cy', '0')

c_symb.set('r', '4')

else: # cross

c_symb = etree.SubElement(c_group, 'path')

c_symb.set('d', 'M {1},-{0} L {0},-{1} L {2},{3} L {0},{1} \

L {3},-{2} L {1},-{0} Z'.format(11.1, 8.5, 2.6, 0))

c_symb.set('style', 'fill:#000000; stroke:none')

def draw_magnets(self, field):

if 'coils' not in field.elements: return

coils = field.elements['coils']

self.__check_symbols()

for coil in coils:

m_group = etree.SubElement(self.symbols, 'g')

self.count_symbols += 1

m_group.set('id', 'magnet{0}'.format(self.count_symbols))

m_group.set('transform',

'translate({0},{1}) rotate({2})'.format(

coil[0], coil[1], degrees(coil[2])))

#### magnet drawing ####

r = coil[3]; l = coil[4]

colors = ['#00cc00', '#ff0000']

SN = ['S', 'N']

if coil[5] < 0.:

colors.reverse()

SN.reverse()

m_defs = etree.SubElement(m_group, 'defs')

m_gradient = etree.SubElement(m_defs, 'linearGradient')

m_gradient.set('id', 'magnetGrad{0}'.format(self.count_symbols))

for attr, val in [['x1', '0'], ['x2', '0'], ['y1', str(coil[3])],

['y2', str(-coil[3])], ['gradientUnits', 'userSpaceOnUse']]:

m_gradient.set(attr, val)

for col, of, opa in [['#000000', '0', '0.125'],

['#ffffff', '0.07', '0.125'], ['#ffffff', '0.25', '0.5'],

['#ffffff', '0.6', '0.2'], ['#000000', '1', '0.33']]:

stop = etree.SubElement(m_gradient, 'stop')

stop.set('stop-color', col)

stop.set('offset', of)

stop.set('stop-opacity', opa)

for i in [0, 1]:

rect = etree.SubElement(m_group, 'rect')

for attr, val in [['x', [-l, 0][i]], ['y', -r],

['width', [2*l, l][i]], ['height', 2 * r],

['style', 'fill:{0}; stroke:none'.format(colors[i])]]:

rect.set(attr, str(val))

rect = etree.SubElement(m_group, 'rect')

for attr, val in [['x', -l], ['y', -r],

['width', 2 * l], ['height', 2 * r],

['style', 'fill:url(#magnetGrad{0}); stroke-width:{1}; stroke-linejoin:miter; stroke:#000000'.format(self.count_symbols, 4. / self.unit)]]:

rect.set(attr, str(val))

for i in [0, 1]:

text = etree.SubElement(m_group, 'text')

for attr, val in [['text-anchor', 'middle'], ['y', -r],

['transform', 'translate({0},{1}) scale({2},-{2})'.format(

[-0.65, 0.65][i] * l, -0.44 * r, r / 100.)],

['style', 'fill:#000000; stroke:none; ' +

'font-size:120px; font-family:Bitstream Vera Sans']]:

text.set(attr, str(val))

text.text = SN[i]

def draw_line(self, fieldline, maxdist=10., linewidth=2.,

linecolor='#000000', attributes=[], arrows_style=None):

'''

self.__check_fieldlines(linecolor, linewidth)

self.count_fieldlines += 1

bounds = {}

bounds['x0'] = -(self.center[0] + 0.5 * linewidth) / self.unit

bounds['y0'] = -(self.height - self.center[1] +

0.5 * linewidth) / self.unit

bounds['x1'] = (self.width - self.center[0] +

0.5 * linewidth) / self.unit

bounds['y1'] = (self.center[1] + 0.5 * linewidth) / self.unit

# fetch the polyline from the fieldline object

polylines = fieldline.get_polylines(self.digits, maxdist, bounds)

if len(polylines) == 0: return

line = etree.Element('path')

if self.fieldlines.get('stroke') != linecolor:

line.set('stroke', linecolor)

if self.fieldlines.get('stroke-width') != str(linewidth / self.unit):

line.set('stroke-width', str(linewidth / self.unit))

for attr, val in attributes:

line.set(attr, val)

#### line drawing ####

path_data = []

for polyline in polylines:

line_points = polyline['path']

for i, p in enumerate(line_points):

# go through all points, draw them if line segment is visible

ptext = '{1:.{0}f},{2:.{0}f}'.format(

int(ceil(self.digits)), p[0], p[1])

if i == 0: path_data.append('M ' + ptext)

else: path_data.append('L ' + ptext)

# close path if possible

if (vabs(polylines[0]['path'][0] - polylines[-1]['path'][-1])

< .1**self.digits):

closed = True

if len(polylines) == 1:

path_data.append('Z')

elif len(polylines) > 1:

# rearrange array cyclic

path_data.pop(0)

while path_data[0][0] != 'M':

path_data.append(path_data.pop(0))

else: closed = False

path = ' '.join(path_data)

line.set('d', path)

if arrows_style == None:

# include path directly into document structure

line.set('id', 'fieldline{0}'.format(self.count_fieldlines))

self.fieldlines.append(line)

else:

line_and_arrows = etree.SubElement(self.fieldlines, 'g')

line_and_arrows.set('id', 'fieldline' + str(self.count_fieldlines))

line_and_arrows.append(line)

line_and_arrows.append(self.__draw_arrows(

arrows_style, linewidth, polylines, linecolor, closed))

def draw_vectors(self, vectors, linewidth=2.,

linecolor='#000000', attributes=[], arrows_style=None):

arrowhead = etree.SubElement( self.__get_defs(), 'marker')

arrowhead.set('id', 'arrowhead')

arrowhead.set('orient', 'auto')

arrowhead.set('refX', '1')

arrowhead.set('refY', '5')

arrowhead.set('markerWidth', '3')

arrowhead.set('markerHeight', '3')

arrowhead.set('viewBox', '0 0 10 10')

path = etree.SubElement(arrowhead, 'path')

path.set('fill', linecolor)

path.set('d', 'M 0 0 L 10 5 L 0 10 z' )

self.fieldvectors = etree.SubElement(self.content, 'g')

self.fieldvectors.set('id', 'fieldvectors')

self.fieldvectors.set('fill', 'none')

self.fieldvectors.set('stroke', linecolor)

self.fieldvectors.set('stroke-width', str(linewidth / self.unit))

for x,y,F in vectors.scaled_vectors:

line = etree.Element('line')

line.set('x1', '%f'%x )

line.set('y1', '%f'%y )

line.set('x2', '%f'%(x+F[0]) )

line.set('y2', '%f'%(y+F[1]) )

line.set('marker-end', 'url(#arrowhead)' )

self.fieldvectors.append(line)

def __draw_arrows(self, arrows_style, linewidth, polylines,

linecolor='#000000', closed=False):

'''

min_arrows = 1; max_arrows = None; arrows_dist = 1.; scale = linewidth

offsets = 4 * [0.5]; fixed_ends = 4 * [False]

if 'min_arrows' in arrows_style:

min_arrows = arrows_style['min_arrows']

if 'max_arrows' in arrows_style:

max_arrows = arrows_style['max_arrows']

if 'dist' in arrows_style:

arrows_dist = arrows_style['dist']

if 'scale' in arrows_style:

scale *= arrows_style['scale']

if 'offsets' in arrows_style:

offsets = arrows_style['offsets']

if 'fixed_ends' in arrows_style:

fixed_ends = arrows_style['fixed_ends']

if scale == 1.: scaletext = ''

else: scaletext = ' scale({0})'.format(scale)

arrows = etree.Element('g')

arrows.set('id', 'arrows' + str(self.count_fieldlines))

for j, polyline in enumerate(polylines):

line_points = polyline['path']

mina = min_arrows

maxa = max_arrows

# measure drawn path length

lines_dist = [0.]

for i in range(1, len(line_points)):

lines_dist.append(lines_dist[-1]

+ vabs(line_points[i] - line_points[i-1]))

offs = [offsets[2], offsets[1]]

fixed = [fixed_ends[2], fixed_ends[1]]

if polyline['start']:

offs[0] = offsets[0]

fixed[0] = fixed_ends[0]

if polyline['end']:

offs[1] = offsets[3]

fixed[1] = fixed_ends[3]

d01 = [0., lines_dist[-1]]

for i in [0, 1]:

if fixed[i]:

d01[i] += offs[i] * arrows_dist * [1., -1.][i]

mina -= 1

if maxa is not None: maxa -= 1

if d01[1] - d01[0] < 0.: break

elif d01[1] - d01[0] == 0.: d_list = [d01[0]]

else:

d_list = []

if fixed[0]: d_list.append(d01[0])

if (maxa == None) or (maxa > 0):

number_intervals = (d01[1] - d01[0]) / arrows_dist

number_offsets = 0.

for i in [0, 1]:

if fixed[i]: number_offsets += .5

else: number_offsets += offs[i] - .5

n = int(number_intervals - number_offsets + 0.5)

n = max(n, mina)

if maxa is not None: n = min(n, maxa)

if n > 0:

d = (d01[1] - d01[0]) / float(n + number_offsets)

if fixed[0]: d_start = d01[0] + d

else: d_start = offs[0] * d

for i in range(n):

d_list.append(d_start + i * d)

if fixed[1]: d_list.append(d01[1])

geo = self.arrow_geo # shortcut

#### arrow drawing ####

for d1 in d_list:

# calculate arrow position and direction

if d1 < 0. or d1 > lines_dist[-1]: continue

d0 = d1 + (geo['x_nock'] * scale + 2.5*linewidth *

(geo['x_tail'] - geo['x_nock']) / geo['width']) / self.unit

if closed and d0 < 0.: d0 = lines_dist[-1] + d0

d2 = d1 + (geo['x_head'] * scale + linewidth *

(geo['x_tail'] - geo['x_head']) / geo['width']) / self.unit

if closed and d2 > lines_dist[-1]: d1 -= lines_dist[-1]

i0, s0 = list_interpolate(lines_dist, d0)

i1, s1 = list_interpolate(lines_dist, d1)

i2, s2 = list_interpolate(lines_dist, d2)

p0 = line_points[i0] + s0 * (line_points[i0+1]-line_points[i0])

p1 = line_points[i1] + s1 * (line_points[i1+1]-line_points[i1])

p2 = line_points[i2] + s2 * (line_points[i2+1]-line_points[i2])

p = None; angle = None

if vabs(p2-p1) <= .1**self.digits or (d2 <= d0 and not closed):

v = line_points[i1+1] - line_points[i1]

p = p1

angle = atan2(v[1], v[0])

else:

v = p2 - p0

p = p0 + sc.dot(p1 - p0, v) * v / vabs(v)**2

angle = atan2(v[1], v[0])

arrow = etree.SubElement(arrows, 'use')

arrow.set('{http://www.w3.org/1999/xlink}href',

'#' + self.__get_arrowname(linecolor))

arrow.set('transform', ('translate({0:.'

+ str(int(ceil(self.digits))) + 'f},{1:.'

+ str(int(ceil(self.digits)))

+ 'f}) rotate({2:.2f})').format(p[0], p[1],

degrees(angle)) + scaletext)

return arrows

def draw_object(self, name, params, group=None):

'''

self.__check_symbols()

if group == None:

obj = etree.SubElement(self.symbols, name)

else:

obj = etree.SubElement(group, name)

for i, j in params.items():

obj.set(str(i), str(j))

return obj

def write(self, filename=None):

# put symbols on top

if 'content' in dir(self):

for element in self.content:

if element.get('id').startswith('symbols'):

self.content.append(element)

# write content to file

if filename == None: filename = self.name

outfile = open(filename + '.svg', 'wb')

outfile.write(etree.tostring(self.svg, xml_declaration=True,

pretty_print=True, encoding='utf-8'))

outfile.close()

'''

def __init__(self, field, start_p, start_v=None, start_d=None,

directions='forward', maxn=1000, maxr=300.0, hmax=1.0,

pass_dipoles=0, path_close_tol=5e-3, bounds_func=None,

stop_funcs=[None, None]):

'''

self.field = field

self.p_start = sc.array(start_p)

self.first_point = self.p_start

self.bounds_func = bounds_func

self.stop_funcs = stop_funcs

if start_v is None: self.v_start = None

else: self.v_start = sc.array(start_v)

if start_d is None: self.d_start = None

else: self.d_start = sc.array(start_d)

self.__create_nodes(directions, maxn, maxr, hmax,

pass_dipoles, path_close_tol)

def __get_nearest_pole(self, p, v=None):

'''

p_near = self.first_point

d_near = vabs(self.first_point - p)

if not v is None: d_near *= 1.3 - cosv(v, self.first_point - p)

type_near = 'start'

mon = []

for ptype, poles in self.field.elements.items():

if ptype not in ['monopoles', 'dipoles'] or len(poles) == 0:

continue

for pole in poles:

d = vabs(pole[:2] - p)

if not v is None: d *= 1.3 - cosv(v, pole[:2] - p)

if d < d_near:

d_near = d

p_near = pole

type_near = ptype

return sc.array(p_near), type_near

def __rkstep(self, p, v, f, h):

'''

k1 = h * v

k2 = h * f(p + k1 / 2.)

k3 = h * f(p + k2 / 2.)

k4 = h * f(p + k3)

p1 = p + (k1 + 2. * (k2 + k3) + k4) / 6.

return p1

def __create_nodes_part(self, sign, maxn, maxr, hmax,

pass_dipoles, path_close_tol):

'''

# p is always the latest position

# v is always the latest normalized velocity

# h is always the latest step size

# l is always the summarized length

err = 5e-8 # error tolerance for integration

f = None

if sign >= 0.: f = self.field.Fn

else: f = lambda r: -self.field.Fn(r)

# first point

p = self.p_start

if not self.v_start is None:

v = vnorm(self.v_start) * sign

else:

v = f(p)

nodes = [{'p':p.copy(), 'v_in':None}]

xtol = 20. * err; ytol = path_close_tol

# initialize loop

h = (sqrt(5) - 1.) / 10.; h_old = h

l = 0.; i = 0

while i < maxn and l < maxr:

i += 1

if len(nodes) == 1 and not self.d_start is None:

# check for start from a dipole

h = vabs(self.d_start)

p = p + self.d_start

v = f(p)

nodes[-1]['v_out'] = h * vnorm(2.0 * vnorm(self.d_start) - v)

nodes.append({'p':p.copy(), 'v_in':h * v})

elif len(nodes) > 1:

# check for special cases

nearest_pole, pole_type = self.__get_nearest_pole(p, v)

vpole = nearest_pole[:2] - p

dpole = vabs(vpole)

vpole /= dpole

cv = cosv(v, vpole); sv = sinv(v, vpole)

if ((dpole < 0.1 or h >= dpole)

and (cv > 0.9 or dpole < ytol)):

# heading for some known special point

if pole_type == 'start':

# is the fieldline about to be closed?

if ((dpole * abs(sv) < ytol) and

(dpole * abs(cv) < xtol) and (l > 1e-3)):

# path is closed

nodes[-1]['v_out'] = None

message(1,'closed at', pretty_vec(p))

break

elif (h > 0.99 * dpole and (cv > 0.9 or

(cv > 0. and dpole * abs(sv) < ytol))):

# slow down

h = max(4.*err, dpole*cv * max(.9, 1-.1*dpole*cv))

if (pole_type == 'monopoles' and

dpole < 0.01 and cv > .996):

# approaching a monopole: end line with x**3 curve

nodes[-1]['v_out'] = vnorm(v) * dpole

v = vnorm(1.5 * vnorm(vpole) -

.5 * vnorm(nodes[-1]['v_out']))

nodes.append({'p':nearest_pole[:2].copy(),

'v_in':v * dpole, 'v_out':None})

l += h

break

if (pole_type == 'dipoles' and

dpole < 0.01 and cv > .996):

# approaching a dipole

m = sign * vnorm(nearest_pole[2:4])

p = nodes[-1]['p'] + 2. * sc.dot(m, vpole) * m * dpole

# approximation by a y=x**1.5 curve

nodes[-1]['v_out'] = 2. * vnorm(v) * dpole

nodes.append({'p':nearest_pole[:2].copy(),

'v_in':sc.zeros(2), 'v_out':sc.zeros(2)})

l += h

# check if the path is being closed

v_end = self.first_point - p

if ((dpole * abs(sinv(v, v_end)) < ytol) and

(dpole * abs(cosv(v, v_end)) < xtol) and l > 1e-3):

# path is closed

nodes[-1]['v_out'] = None

break

if pass_dipoles == 0:

nodes[-1]['v_out'] = None

break

if pass_dipoles > 0:

pass_dipoles -= 1

v = f(p)

nodes.append({'p':p.copy(), 'v_in':2.*vnorm(v)*dpole})

l += h

continue

# buckle detection at unknown places

elif h < 0.01:

# check change rate of curvature

hh = h * 3.

v0 = f(p + hh / 2. * v)

v1 = f(p + hh * v)

angle0 = atan2(v[1], v[0])

angle1 = atan2(v0[1], v0[0])

angle2 = atan2(v1[1], v1[0])

a0 = angle_dif(angle1, angle0)

a1 = angle_dif(angle2, angle1)

adif = angle_dif(a1, a0)

corner_limit = 1e4

if abs(adif) / hh**2 > corner_limit:

# assume a corner here

if abs(a0) >= abs(a1):

h0 = 0.; h1 = hh / 2.

vm = vnorm(vnorm(v) + vnorm(v0))

else:

h0 = hh / 2.; h1 = hh

vm = vnorm(vnorm(v0) + vnorm(v1))

if vabs(vm)==0.: vm = vnorm(sc.array([v0[1], -v0[0]]))

hc = op.brentq(lambda hc: sinv(f(p+hc*v), vm), h0, h1)

v2 = f(p + hc / 2. * v)

if sinv(f(p), vm) * sinv(f(p + 2.*hc*v2), vm) <= 0.:

hc = op.brentq(lambda hc: sinv(f(p + hc * v2),

vm), 0., 2. * hc)

nodes[-1]['v_out'] = vnorm(nodes[-1]['v_in']) * hc

# create a corner

# use second-order formulas instead of runge-kutta

p += hc * v2

message(1,'corner at', pretty_vec(p))

v = vnorm(2. * v2 - v)

nodes.append({'p':p.copy(),'v_in':v*hc,'corner':True})

l += h

# check if the path is being closed

v_end = self.first_point - p

if ((dpole * abs(sinv(v, v_end)) < ytol) and

(dpole * abs(cosv(v, v_end)) < xtol) and l > 1e-3):

# path is closed

nodes[-1]['v_out'] = None

break

# check area after the corner

# lengths are chosen to ensure corner detection

p0 = p + hh * .2 * f(p + hh * .2 * v1); va0 = f(p0)

p1 = p0 + hh * .4 * va0; va1 = f(p1)

p2 = p1 + hh * .4 * va1; va2 = f(p2)

angle0 = atan2(va0[1], va0[0])

angle1 = atan2(va1[1], va1[0])

angle2 = atan2(va2[1], va2[0])

a0 = angle_dif(angle1, angle0)

a1 = angle_dif(angle2, angle1)

adif = angle_dif(a1, a0)

if (abs(adif) / (.8*hh)**2 > corner_limit or

abs(a0) + abs(a1) >= pi / 2.):

message(1,'end edge at', pretty_vec(p))

# direction after corner changes again -> end line

nodes[-1]['v_out'] = None

break

vm = vnorm(1.25 * va1 - 0.25 * va2)

v = f(p + hh * vm)

nodes[-1]['v_out'] = vnorm(2. * vm - v) * hh

p += vm * hh

nodes.append({'p':p.copy(), 'v_in':v * hh})

l += h

# make single and double runge-kutta step

p11 = self.__rkstep(p, vnorm(v), f, h)

p21 = self.__rkstep(p, vnorm(v), f, h / 2.)