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