Coverage for src / thunderfish / efield.py: 0%
181 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-15 17:50 +0000
1"""Simulations of spatial electric fields.
3## Electric monopoles
5For simulating the spatial geometry of electric fields generated by electric fishes
6and perturbed by objects, first generate monopoles and charges:
8- `efish_monopoles()`: monopoles for simulating the electric field of an electric fish.
9- `object_monopoles()`: monopoles for simulating a circular object.
11## Potential, electric field, and field lines
13- `epotential()`: simulation of electric field potentials.
14- `epotential_meshgrid()`: simulation of electric field potentials on a mesh grid.
15- `efield()`: simulation of electric field.
16- `efield_meshgrid()`: simulation of electric field on a mesh grid.
17- `projection()`: projection of electric field on surface normals.
18- `fieldline()`: compute an electric field line.
20## Visualization
22- `squareroot_transform()`: square-root transformation keeping the sign.
23- `plot_fieldlines()`: plot field lines with arrows.
24"""
26import numpy as np
27import matplotlib.pyplot as plt
29from matplotlib.patches import FancyArrowPatch
32def efish_monopoles(pos=(0, 0), direction=(1, 0), size=10.0, bend=0, nneg=1):
33 """Monopoles for simulating the electric field of an electric fish.
35 This implements the model published in
36 Chen, House, Krahe, Nelson (2005) "Modeling signal and background
37 components of electrosensory scenes", J Comp Physiol A 191: 331-345
39 Ten monopoles per unit size are uniformly distributed along the fish's body axis.
40 The first (tail) nneg monopoles get negative charges that equal in sum
41 the sum of the positive unit charges of the remaining (head) monopoles.
42 The strength of the dipole increases linearly with fish size.
44 Pass the returned monopole positions and charges on to the epotential() function
45 to simulate the resulting electric field potentials, to the efield() function
46 to simulate the electric field, or to object_monopoles() to add an object.
48 Parameters
49 ----------
50 pos: tuple of floats
51 Coordinates of the fish's position (its center).
52 The number of elements in pos set the number of dimensions to be used.
53 direction: tuple of floats
54 Coordinates of a vector defining the orientation of the fish.
55 Missing dimensions are filled in with zeros.
56 Note: currently only rotations in the x-y plane are implemented.
57 size: float
58 Size of the fish. Per size unit 10 monopols are distributed along
59 the fish's body axis.
60 bend: float
61 Bending angle of the fish's tail in degree.
62 nneg: int
63 Number of negative charges to be used. The remaining ones are positively charged.
65 Returns
66 -------
67 poles: 2D array of floats
68 Positions of the monopoles with n-dimensional coordinates
69 as specified by the number of elements in pos.
70 charges: array of floats
71 The charge of each monopole.
73 Example
74 -------
75 ```
76 fish1 = ((-8, -5), (1, 0.5), 18.0, -25)
77 poles1 = efish_monopoles(*fish1)
78 ```
79 """
80 n = int(10*size)
81 npos = n - nneg
82 ppx = 0.1
83 pos = np.asarray(pos)
84 dirv = np.zeros(len(pos))
85 dirv[:len(direction)] = direction
86 charges = np.ones(n)
87 charges[:nneg] = -float(npos)/nneg
88 poles = np.zeros((n, len(pos)))
89 poles[:,0] = np.arange(-n//2, -n//2+n)*ppx
90 if np.abs(bend) > 1.e-8:
91 xm = -np.min(poles[:,0]) # tip of fish tail
92 r = -180.0*xm/bend/np.pi # radius of circle on which to bend the tail
93 xp = poles[poles[:,0]<0.0,0] # all negative x coordinates of poles
94 beta = xp/r # angle on circle for each x coordinate
95 poles[poles[:,0]<0.0,0] = -np.abs(r*np.sin(beta)) # transformed x coordinates
96 poles[poles[:,0]<0.0,1] = r*(1.0-np.cos(beta)) # transformed y corrdinates
97 # rotation matrix:
98 theta = -np.arctan2(dirv[1], dirv[0])
99 c = np.cos(theta)
100 s = np.sin(theta)
101 rm = np.array(((c, -s), (s, c)))
102 # rotation:
103 poles[:,:2] = np.dot(poles[:,:2], rm)
104 # translation:
105 poles += pos
106 return poles, charges
109def object_monopoles(pos=(0, 0), radius=1.0, chi=1.0, *args):
110 """Monopoles for simulating a circular object.
112 The circular object is approximated by an induced dipole, as
113 sugested by Rasnow B (1996) "The effects of simple objects on the
114 electric field of Apteronotus", J Comp Physiol A 178:397-411 and
115 Chen, House, Krahe, Nelson (2005) "Modeling signal and background
116 components of electrosensory scenes", J Comp Physiol A 191: 331-345.
118 Pass the returned monopole positions and charges on to the
119 epotential() function to simulate the resulting electric field
120 potentials or to the efield() function to simulate the electric
121 field.
123 Two monopoles with charges q and -q separated by dx form a dipole
124 with dipole moment p = q dx. According to Chen et al (2005), this
125 dipole moment should equal chi*radius**3*E_obj, where E_obj is the
126 electric field generated by the fishes at the position of the
127 object. We normalize E_obj and multiply it with a small number eps
128 to get dx. Accordingly, we have to set q to chi*radius**3
129 |E_obj|/eps.
131 Parameters
132 ----------
133 pos: tuple of floats
134 Coordinates of the fish's position (its center).
135 The number of dimensions must be the same as the one of the poles
136 passed on in args.
137 radius: float
138 Radius of the small circular object.
139 chi: float
140 Electrical contrast. Unity for a perfect conductor, -0.5 for a
141 perfect insulator and zero if the electrical impedance of the sphere
142 matches that of the surrounding water.
143 args: list of tuples
144 Each tuple contains as the first argument the position of
145 monopoles (2D array of floats), and as the second argument the
146 corresponding charges (array of floats) of electric fish. Use
147 efish_monopoles() to generate monopoles and corresponding charges.
149 Returns
150 -------
151 poles: 2D array of floats
152 Positions of the monopoles.
153 charges: array of floats
154 The charge of each monopole.
156 Example
157 -------
158 ```
159 fish1 = ((-8, -5), (1, 0.5), 18.0, -25)
160 fish2 = ((12, 3), (0.8, 1), 20.0, 20)
161 poles1 = efish_monopoles(*fish1)
162 poles2 = efish_monopoles(*fish2)
163 poles3 = object_monopoles((-6, 0), 1.0, -0.5, poles1, poles2)
164 ```
165 """
166 eps = 0.1 # distance of the two monopoles
167 pos = np.asarray(pos)
168 # electric field at object position:
169 eobj = efield(pos, *args)
170 eobjnorm = np.linalg.norm(eobj)
171 # induced dipole:
172 charges = np.ones(2)*chi*radius**3*eobjnorm/eps
173 charges[0] = -charges[0]
174 poles = np.zeros((2, len(pos)))
175 poles[0,:] = -eobj*0.5*eps/eobjnorm # distance between monopoles
176 poles[1,:] = +eobj*0.5*eps/eobjnorm # distance between monopoles
177 poles += pos # translation to required position
178 return poles, charges
181def epotential(pos, *args):
182 """Simulation of electric field potentials.
184 Parameters
185 ----------
186 pos: 2D array of floats
187 Each row contains the coordinates (2D or 3D)
188 for which the potential should be computed.
189 args: list of tuples
190 Each tuple contains as the first argument the position of monopoles
191 (2D array of floats), and as the second argument the corresponding charges
192 (array of floats). Use efish_monopoles() to generate monopoles and
193 corresponding charges.
195 Returns
196 -------
197 pot: 1D array of float
198 The potential for each position in `pos`.
199 """
200 pos = np.asarray(pos)
201 pot = np.zeros(len(pos))
202 for poles, charges in args:
203 for p, c in zip(poles, charges):
204 r = pos - p
205 rnorm = np.linalg.norm(r, axis=1)
206 rnorm[np.abs(rnorm) < 1e-12] = 1.0e-12
207 pot += c/rnorm
208 return pot
211def epotential_meshgrid(xx, yy, zz, *args):
212 """Simulation of electric field potentials on a mesh grid.
214 This is a simple wrapper for epotential().
216 Parameters
217 ----------
218 xx: 2D array of floats
219 Range of x coordinates as returned by numpy.meshgrid().
220 yy: 2D array of floats
221 Range of y coordinates as returned by numpy.meshgrid().
222 zz: None or 2D array of floats
223 z coordinates on the meshgrid defined by xx and yy.
224 If provided, poles in args must be 3D.
225 If None then treat it as a 2D problem with poles in args providing 2D coordinate.
226 args: list of tuples
227 Each tuple contains as the first argument the position (2D or 3D) of monopoles
228 (2D array of floats), and as the second argument the corresponding charges
229 (array of floats). Use efish_monopoles() to generate monopoles and
230 corresponding charges.
232 Returns
233 -------
234 pot: 2D array of floats
235 The potential for the mesh grid defined by xx and yy and evaluated
236 at (xx, yy, zz).
238 Example
239 -------
240 ```
241 fig, ax = plt.subplots()
242 maxx = 30.0
243 maxy = 27.0
244 x = np.linspace(-maxx, maxx, 200)
245 y = np.linspace(-maxy, maxy, 200)
246 xx, yy = np.meshgrid(x, y)
247 fish1 = ((-8, -5), (1, 0.5), 18.0, -25)
248 fish2 = ((12, 3), (0.8, 1), 20.0, 20)
249 poles1 = efish_monopoles(*fish1)
250 poles2 = efish_monopoles(*fish2)
251 poles3 = object_monopoles((-6, 0), 1.0, -0.5, poles1, poles2)
252 allpoles = (poles1, poles2, poles3)
253 # potential:
254 pot = epotential_meshgrid(xx, yy, None, *allpoles)
255 thresh = 0.65
256 zz = squareroot_transform(pot/200, thresh)
257 levels = np.linspace(-thresh, thresh, 16)
258 ax.contourf(x, y, -zz, levels, cmap='RdYlBu')
259 ax.contour(x, y, -zz, levels, zorder=1, colors='#707070',
260 linewidths=0.1, linestyles='solid')
261 plt.show()
262 ```
263 """
264 pos = np.vstack((xx.ravel(), yy.ravel())).T
265 pot = epotential(pos, *args)
266 return pot.reshape(xx.shape)
269def efield(pos, *args):
270 """Simulation of electric field given a set of electric monopoles.
272 Parameters
273 ----------
274 pos: array of floats
275 Each row contains the coordinates (2D or 3D)
276 for which the potential should be computed.
277 A single (1D) position is also accepted.
278 args: list of tuples
279 Each tuple contains as the first argument the position of monopoles
280 (2D array of floats), and as the second argument the corresponding charges
281 (array of floats). Use efish_monopoles() to generate monopoles and
282 corresponding charges.
284 Returns
285 -------
286 field: array of floats
287 The electric field components for each position in `pos`.
288 """
289 pos = np.asarray(pos)
290 onedim = len(pos.shape) == 1
291 if onedim:
292 pos = pos.reshape(-1, len(pos))
293 field = np.zeros(pos.shape)
294 for poles, charges in args:
295 for p, c in zip(poles, charges):
296 r = pos - p
297 rnorm = np.linalg.norm(r, axis=1)
298 rnorm[np.abs(rnorm) < 1e-12] = 1.0e-12
299 fac = c/rnorm**3
300 field += r*fac[:,np.newaxis]
301 return field[0] if onedim else field
304def efield_meshgrid(xx, yy, zz, *args):
305 """Simulation of electric field on a mesh grid.
307 This is a simple wrapper for efield().
309 Parameters
310 ----------
311 xx: 2D array of floats
312 Range of x coordinates as returned by numpy.meshgrid().
313 yy: 2D array of floats
314 Range of y coordinates as returned by numpy.meshgrid().
315 zz: None or 2D array of floats
316 z coordinates on the meshgrid defined by xx and yy.
317 If provided, poles in args must be 3D.
318 If None then treat it as a 2D problem with poles in args providing 2D coordinate.
319 args: list of tuples
320 Each tuple contains as the first argument the position (2D or 3D) of monopoles
321 (2D array of floats), and as the second argument the corresponding charges
322 (array of floats). Use efish_monopoles() to generate monopoles and
323 corresponding charges.
325 Returns
326 -------
327 pot: 2D array of floats
328 The potential for the mesh grid defined by xx and yy and evaluated
329 at (xx, yy, zz).
331 Returns
332 -------
333 fieldx: 2D array of floats
334 The x-coordinate of the electric field for the mesh grid
335 defined by xx and yy and evaluated at (xx, yy, zz).
336 fieldy: 2D array of floats
337 The y-coordinate of the electric field for the mesh grid
338 defined by xx and yy and evaluated at (xx, yy, zz).
339 fieldz: 2D array of floats
340 The z-coordinate of the electric field for the mesh grid
341 defined by xx and yy and evaluated at (xx, yy, zz).
342 This is only returned if zz is not None.
344 Example
345 -------
346 ```
347 fig, ax = plt.subplots()
348 maxx = 30.0
349 maxy = 27.0
350 x = np.linspace(-maxx, maxx, 40)
351 y = np.linspace(-maxy, maxy, 40)
352 xx, yy = np.meshgrid(x, y)
353 fish1 = ((-8, -5), (1, 0.5), 18.0, -25)
354 fish2 = ((12, 3), (0.8, 1), 20.0, 20)
355 poles1 = efish_monopoles(*fish1)
356 poles2 = efish_monopoles(*fish2)
357 poles3 = object_monopoles((-6, 0), 1.0, -0.5, poles1, poles2)
358 allpoles = (poles1, poles2, poles3)
359 fieldx, fieldy = efield_meshgrid(xx, yy, None, *allpoles)
360 u = squareroot_transform(fieldx, 0)
361 v = squareroot_transform(fieldy, 0)
362 ax.quiver(qx, qy, u, v, units='xy', angles='uv', scale=2, scale_units='xy',
363 width=0.07, headwidth=5)
364 ```
365 """
366 if zz is None:
367 pos = np.vstack((xx.ravel(), yy.ravel())).T
368 ef = efield(pos, *args)
369 return ef[:,0].reshape(xx.shape), ef[:,1].reshape(xx.shape)
370 else:
371 pos = np.vstack((xx.ravel(), yy.ravel(), zz.ravel())).T
372 ef = efield(pos, *args)
373 return ef[:,0].reshape(xx.shape), ef[:,1].reshape(xx.shape), ef[:,2].reshape(xx.shape)
376def projection(ex, ey, ez, nx, ny, nz):
377 """Projection of electric field on surface normals.
379 Parameters
380 ----------
381 ex: array of floats
382 x-coordinates of the electric field.
383 ey: array of floats
384 y-coordinates of the electric field.
385 ez: array of floats
386 z-coordinates of the electric field.
387 nx: array of floats
388 x-coordinates of the surface normals.
389 ny: array of floats
390 y-coordinates of the surface normals.
391 nz: array of floats
392 z-coordinates of the surface normals.
393 """
394 ef = np.vstack((ex.ravel(), ey.ravel(), ez.ravel())).T
395 nf = np.vstack((nx.ravel(), ny.ravel(), nz.ravel())).T
396 proj = np.sum(ef*nf, axis=1)
397 return proj.reshape(ex.shape)
400def fieldline(pos0, bounds, *args, eps=0.1, maxiter=1000):
401 """Compute an electric field line.
403 From the initial position `pos0` the field line is computed in both directions
404 until it leaves the area defined by `bounds`.
406 Parameters
407 ----------
408 pos0: array of floats
409 Initial position for computing the field line.
410 bounds: None or 2D array of floats
411 If not None, stop integration of the field line if it exceeds bounds.
412 First row are the minimum coordinates and second row the maximum coordinates.
413 args: list of tuples
414 Each tuple contains as the first argument the position of monopoles
415 (2D array of floats), and as the second argument the corresponding charges
416 (array of floats). Use efish_monopoles() to generate monopoles and
417 corresponding charges.
418 eps: float
419 Stepsize in unit of the coordinates.
420 maxiter: int
421 Maximum number of iteration steps.
423 Returns
424 -------
425 fl: 2D array of floats
426 Coordinates of the computed field line.
428 Example
429 -------
430 ```
431 fig, ax = plt.subplots()
432 fish1 = ((-8, -5), (1, 0.5), 18.0, -25)
433 fish2 = ((12, 3), (0.8, 1), 20.0, 20)
434 poles1 = efish_monopoles(*fish1)
435 poles2 = efish_monopoles(*fish2)
436 fl = fieldline((0, -16), [[-maxx, -maxy], [maxx, maxy]], poles1, poles2)
437 plot_fieldlines(ax, [fl], 5, color='b', lw=2)
438 plt.show()
439 ```
440 """
441 bounds = np.asarray(bounds)
442 p = np.array(pos0)
443 n = maxiter//2
444 # forward integration:
445 flf = np.zeros((n, len(pos0)))
446 for i in range(len(flf)):
447 flf[i,:] = p
448 if np.any(p < bounds[0,:]) or np.any(p > bounds[1,:]) or (bounds is not None and
449 i >= 5 and np.all((flf[i,:] - flf[i-1,:])*(flf[i-1,:] - flf[i-2,:])<0)):
450 flf = flf[:i,:]
451 break
452 uv = efield(p, *args)
453 uv /= np.linalg.norm(uv)
454 p = p + eps*uv
455 # backward integration:
456 p = np.array(pos0)
457 flb = np.zeros((n, len(pos0)))
458 for i in range(len(flb)):
459 flb[i,:] = p
460 if np.any(p < bounds[0,:]) or np.any(p > bounds[1,:]) or (bounds is not None and
461 i >= 5 and np.all((flb[i,:] - flb[i-1,:])*(flb[i-1,:] - flb[i-2,:])<0)):
462 flb = flb[:i,:]
463 break
464 uv = efield(p, *args)
465 uv /= np.linalg.norm(uv)
466 p = p - eps*uv
467 fl = np.vstack((flb[::-2], flf[::2]))
468 return fl
471def squareroot_transform(values, thresh=0.0):
472 """Square-root transformation keeping the sign.
474 Takes the square root of positive values and takes the square root
475 of the absolute values of negative values and negates the results.
477 Then truncate symmetrically both positive and negative values to
478 a threshold.
480 The resulting transformed values give nice contour lines in a
481 contour plot.
483 Parameters
484 ----------
485 values: array of float
486 The values to be transformed, i.e. potentials or field strengths.
487 thresh: float or None
488 Maximum absolute value of the returned values.
489 Must be positive!
490 If thresh equals zero, then do not apply treshold.
491 If None, take the smaller of the maximum of the
492 positive values or of the absolute negative values.
494 Returns
495 -------
496 values: array of float
497 The transformed (square-rooted and thresholded) values.
498 """
499 values = np.array(values)
500 sel = values>=0.0
501 values[sel] = values[sel]**0.5
502 values[np.logical_not(sel)] = -((-values[np.logical_not(sel)])**0.5)
503 if thresh is None:
504 thresh = min(np.max(values), -np.min(values))
505 if thresh > 0:
506 values[values>thresh] = thresh
507 values[values<-thresh] = -thresh
508 return values
511def plot_fieldlines(ax, flines, pos=5, **kwargs):
512 """Plot field lines with arrows.
514 Parameters
515 ----------
516 ax: matplotlib axes
517 Axes in which to plot the field lines.
518 flines: list of 2D arrays
519 The field lines.
520 pos: float
521 The position of the arrow on the field line in units of the coordinates.
522 **kwargs: key word arguments
523 Passed on to plot().
524 Applies optional zorder argument also to arrow.
525 """
526 xmin, xmax = ax.get_xlim()
527 ymin, ymax = ax.get_ylim()
528 dx = 0.05*np.abs(xmax-xmin)
529 dy = 0.05*np.abs(ymax-ymin)
530 akwargs = dict()
531 if 'zorder' in kwargs:
532 akwargs['zorder'] = kwargs['zorder']
533 lw = 1
534 if 'lw' in kwargs:
535 lw = kwargs['lw']
536 if 'linewidth' in kwargs:
537 lw = kwargs['linewidth']
538 for fl in flines:
539 ax.plot(fl[:,0], fl[:,1], **kwargs)
540 # arrows:
541 d = np.diff(fl, axis=0)
542 dd = np.linalg.norm(d, axis=1)
543 dist = np.cumsum(dd)
544 if dist[-1] >= 6:
545 idx0 = np.argmin(np.abs(dist-pos))
546 if (np.abs(fl[0,0]-xmin)<dx or np.abs(fl[0,0]-xmax)<dx or
547 np.abs(fl[0,1]-ymin)<dy or np.abs(fl[0,1]-ymax)<dy):
548 idx0 = np.argmin(np.abs(dist[-1]-dist-pos))
549 idx1 = np.argmin(np.abs(dist-0.5*dist[-1]))
550 idx = min(idx0, idx1)
551 adx = fl[idx+1,:] - fl[idx,:]
552 ndx = np.linalg.norm(adx)
553 if ndx < 1e-10:
554 continue
555 adx /= ndx
556 posa = fl[idx+1,:] - 0.1*min(dx,dy)*adx
557 posb = fl[idx+1,:]
558 arrow = FancyArrowPatch(posA=posa, posB=posb, shrinkA=0, shrinkB=0,
559 arrowstyle='fancy', mutation_scale=8*lw,
560 connectionstyle='arc3', fill=True,
561 color=kwargs['color'], **akwargs)
562 ax.add_patch(arrow)
565def main():
566 fig, ax = plt.subplots()
567 maxx = 30.0
568 maxy = 27.0
569 x = np.linspace(-maxx, maxx, 200)
570 y = np.linspace(-maxy, maxy, 200)
571 xx, yy = np.meshgrid(x, y)
572 fish1 = ((-8, -5), (1, 0.5), 18.0, -25)
573 fish2 = ((12, 3), (0.8, 1), 20.0, 20)
574 poles1 = efish_monopoles(*fish1)
575 poles2 = efish_monopoles(*fish2)
576 poles3 = object_monopoles((-6, 0), 1.0, -0.5, poles1, poles2)
577 allpoles = (poles1, poles2, poles3)
578 # potential:
579 pot = epotential_meshgrid(xx, yy, None, *allpoles)
580 thresh = 0.65
581 zz = squareroot_transform(pot/200, thresh)
582 levels = np.linspace(-thresh, thresh, 16)
583 ax.contourf(x, y, -zz, levels, cmap='RdYlBu')
584 ax.contour(x, y, -zz, levels, zorder=1, colors='#707070',
585 linewidths=0.1, linestyles='solid')
586 # electric field vectors:
587 n = 5
588 qx, qy = np.meshgrid(x[n::2*n], y[n::2*n])
589 fieldx, fieldy = efield_meshgrid(qx, qy, None, *allpoles)
590 u = squareroot_transform(fieldx, 0)
591 v = squareroot_transform(fieldy, 0)
592 ax.quiver(qx, qy, u, v, units='xy', angles='uv', scale=2, scale_units='xy',
593 width=0.07, headwidth=5)
594 # field line:
595 bounds = [[-maxx, -maxy], [maxx, maxy]]
596 fl = fieldline((0, -16), bounds, *allpoles)
597 plot_fieldlines(ax, [fl], 5, color='b', lw=2)
598 plt.show()
601if __name__ == '__main__':
602 main()