File:Ironfilings ringmagnet.svg

#!/usr/bin/python3
# -*- coding: utf8 -*-

try:
    import svgwrite
except ImportError:
    print('requires svgwrite library: https://pypi.org/project/svgwrite/')
    # documentation at https://svgwrite.readthedocs.io/
    exit(1)


import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import minimize
from math import *


name = 'Ironfilings_ringmagnet'
size = 600, 600
R1 = 100
R2 = 200
L = 100
n_filings = 2500
l_filings = 24
w_filings = 3.5


def cel(kc, p, a, b):
    """
    Bulirsch complete elliptic integral
    """
    
    if kc == 0.:
        return nan
    
    tol = 1e-9 # actual relative error will be tol**2
    k = kc = fabs(kc)
    m = 1.
    
    if p > 0.:
        p = sqrt(p)
        b /= p
    else:
        f = kc * kc
        g = 1. - p
        q = (1. - f) * (b - a * p)
        f -= p
        p = sqrt(f / g)
        a = (a - b) / g
        b = a * p - q / (g * g * p)
    
    for i in range(11):
        f = a
        a += b / p
        g = k / p
        b = 2. * (b + f * g)
        p += g
        g = m
        m += kc
        
        if fabs(g - kc) <= g * tol:
            break
        
        kc = 2. * sqrt(k)
        k = kc * m
    
    return pi * .5 * (a * m + b) / (m * (m + p))


def Bfield_barmagnet(xy, R, L, M):
    # www.doi.org/10.1119/1.3256157
    rho, z = xy
    Rm, Rp = R - rho, R + rho
    zm, zp = z - L / 2, z + L / 2
    Rmzm = hypot(Rm, zm)
    Rmzp = hypot(Rm, zp)
    Rpzm = hypot(Rp, zm)
    Rpzp = hypot(Rp, zp)
    g = Rm / Rp
    km = Rmzm / Rpzm
    kp = Rmzp / Rpzp
    Frhom = cel(km, 1., 1., -1.) / Rpzm
    Frhop = cel(kp, 1., 1., -1.) / Rpzp
    Fzm = cel(km, g * g, 1., g) * zm / Rpzm
    Fzp = cel(kp, g * g, 1., g) * zp / Rpzp
    return M * R / pi * np.array((Frhop - Frhom, (Fzp - Fzm) / Rp))


def Bfield(xy):
    return Bfield_barmagnet(xy, R2, L, 1.) - Bfield_barmagnet(xy, R1, L, 1.)


def inside(p, w1, w2, h2):
    return fabs(p[0]) >= w1 and fabs(p[0]) <= w2 and fabs(p[1]) <= h2


def vnorm(x):
    return x / hypot(*x)


def bezier(field, p, l):
    # returns control points of a Bezier curve approximating the Bfield at p
    f = lambda t, xy: vnorm(field(xy))
    p2, p0 = solve_ivp(f, (0, -l / 2), p, t_eval=[-l / 4, -l / 2]).y.T
    p3, p1 = solve_ivp(f, (0, l / 2), p, t_eval=[l / 4, l / 2]).y.T
    v0 = vnorm(field(p0))
    v1 = vnorm(field(p1))
    
    if hypot(*(v1 - v0)) < 0.01:
        l0, l1 = l / 3, l / 3
    else:
        def err(x):
            l0, l1, t2, t4, t3 = x
            dist = 0.
            for t, pref in ((t2, p2), (t4, p), (t3, p3)):
                pc = (1-t)**3 * p0 + 3*(1-t)**2*t * (p0+l0*v0) + 3*(1-t)*t**2 * (p1-l1*v1) + t**3 * p1
                dist += hypot(*(pc - pref))**2
            return dist
        
        l0, l1 = minimize(err, [l / 3, l / 3, 0.25, 0.5, 0.75],
            bounds=((0, l), (0, l), (0, 1), (0, 1), (0, 1))).x[:2]
    
    p0c = p0 + v0 * l0
    p1c = p1 - v1 * l1
    
    return [p0, p0c, p1c, p1]


doc = svgwrite.Drawing(name + '.svg', profile='full', size=size)
doc.set_desc(name, 'https://commons.wikimedia.org/wiki/File:' + name +
    '.svg\nrights: Creative Commons Attribution ShareAlike license')
clip = doc.defs.add(doc.clipPath(id='image_clip'))
clip.add(doc.rect(insert=(-size[0]/2., -size[1]/2.), size=size))
doc.add(doc.rect(id='background', insert=(0, 0), size=size, fill='#ffffff', stroke='none'))
g = doc.add(doc.g(id='image', clip_path='url(#image_clip)',
    transform='translate({:.0f}, {:.0f}) scale(1,-1)'.format(size[0]/2., size[1]/2.)))
magnet_back = g.add(doc.g(id='magnet_back'))
filings = g.add(doc.g(id='iron-filings', fill='none', stroke='black',
    stroke_width=2, stroke_linecap='round'))
magnet_front = g.add(doc.g(id='magnet_front'))
mgrad = doc.defs.add(doc.linearGradient(id="magnetGrad",
    start=(0,0), end=(1,0), gradientUnits="objectBoundingBox"))
for c, of, op in [['#000000', '0', '0.33'], ['#ffffff', '0.4', '0.2'],
        ['#ffffff', '0.75', '0.5'], ['#ffffff', '0.93', '0.125'],
        ['#000000', '1', '0.125']]:
    mgrad.add_stop_color(of, c, op)
for x0, x1, isbg in [[-R2, -R1, False], [-R1, R1, True], [R1, R2, False]]:
    if isbg:
        magnet = magnet_back
        fill = 'url(#magnetGrad)'
        colors = ['#00cc00', '#ff0000']
    else:
        magnet = magnet_front
        fill = 'none'
        colors = ['#49da49', '#ff4949']
    for i in [0, 1]:
        magnet.add(doc.rect(insert=(x0, [-L/2, 0][i]), size=(x1-x0, [L, L/2][i]),
            fill=colors[i], stroke='none'))
    magnet.add(doc.rect(insert=(x0, -L/2), size=(x1 - x0, L), fill=fill,
            stroke='#000000', stroke_width=4, stroke_linejoin='miter'))


Bmax = 1.2 * hypot(*Bfield([R1 - 1e-6, 0]))

i_filings = 0
while i_filings < n_filings:
    x = np.random.uniform(-size[0]/2. - l_filings/2., size[0]/2. + l_filings/2.)
    y = np.random.uniform(-size[1]/2. - l_filings/2., size[1]/2. + l_filings/2.)
    l = np.random.uniform(l_filings*0.5, l_filings)
    
    if inside([x, y], R1, R2, L/2. - l/2):
        continue
    
    B = Bfield([x, y])
    Brel = min(hypot(*B) / Bmax, 1)
    line_density = Brel**(2/3)
    line_width = w_filings * Brel**(1/3)
    
    # use rejection sampling to reproduce field line density
    if np.random.random() >= line_density:
        continue
    
    points = bezier(Bfield, [x, y], l)
    
    if all([inside(p, R1, R2, L/2.) for p in [[x, y]] + points]):
        continue
    if all([not inside(p, -1, size[0]/2., size[1]/2.) for p in [[x, y]] + points]):
        continue
    
    filings.add(doc.path(
        d='M {:.1f},{:.1f} C {:.1f},{:.1f} {:.1f},{:.1f} {:.1f},{:.1f}'.format(
        *points[0], *points[1], *points[2], *points[3]),
        stroke_width='{:.1f}'.format(line_width)))
    i_filings += 1
    print(i_filings, end=' ', flush=True)


doc.save(pretty=True)