1 | initial version |

This is more a comment and code review than an answer, but since it's a bit long I'm using the answer field.

Below I'll use `r`

for `sqrt(x^2+y^2)`

and `C_R`

for the circle of radius `R`

centered at the origin.

In the definition of

`C`

you don't need parametric plots. Use Sage's built-in`circle`

command.It looks like

`c1`

,`c3`

,`c4`

are meant to be discretized circles of radius 1, 3, and 4.In their definition, it seems that there should be no

`*u`

inside the`cos`

and`sin`

.Are you using discretized circles on purpose, or as a workaround for continuous circles?

If it's a workaround for continuous circles, note that

`dist((x,y),C_R)`

is just`abs(R-r)`

.To see that, draw the circles

`C_R`

and the circle of radius`r`

centered at`(x,y)`

.The quantity

`d = setdist(x,S)`

is always non-negative, so`min(1,1-d)`

is always`1-d`

.So your functions

`h0`

and`h1`

only differ if`1 <= r <= 3`

, and in this case`h0 = 0`

and`h1 = 1`

.Did you mean something else in the definition of

`h1`

?

Assuming discretized circles were a workaround for continuous circles, I would rewrite your code as follows.

```
# Use "real double field" floating-point values for efficiency
zero = RDF.zero()
one = RDF.one()
four = 4 * one
twopi = 2 * RDF.pi()
# Circles as graphic objects
C1 = circle((0,0), 1)
C3 = circle((0,0), 3)
C4 = circle((0,0), 4)
C = C1 + C3 + C4
# useful functions
def norm2((x,y)):
"""
Norm squared of `(x,y)`
"""
return x^2 + y^2
def norm((x,y)):
"""
Norm of `(x,y)`
"""
return norm2((x,y)).sqrt()
# h0 and h1
def h0(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
Let C_n be the circle of radius n, for n = 1, 3, 4.
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_1
* 1 between C_1 and C_3
* 4 - r between C_3 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < 1:
return r
if 1 <= r <= 3:
return one
if 3 <= r <= 4:
return four - r
return zero
def h1(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
Let C_n be the circle of radius n, for n = 1, 3, 4.
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_1
* 0 between C_1 and C_3
* 4 - r between C_3 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < 1:
return r
if 1 <= r <= 3:
return zero
if 3 <= r <= 4:
return four - r
return zero
def dx(x,y):
return y * h0(x,y)
def dy(x,y):
return -x * h0(x,y) + y * h1(x,y)
```

If discretized circles are relevant, varying the number of points may help illustrate your point, so use a parameter.

```
# discretized circles
n = 99
DC1 = [(cos(i*twopi/n), sin(i*twopi/n)) for i in xrange(n)]
DC3 = [(3*x, 3*y) for (x,y) in DC1]
DC4 = [(4*x, 4*y) for (x,y) in DC1]
```

We would also need to rewrite the functions `h0`

and `h1`

using these discretized circles.

For that, I would rewrite your auxiliary functions as follows:

```
# more useful functions
def dist2((x,y),(xx,yy)):
"""
Distance squared from `(x,y)` to `(xx,yy)`
"""
return norm2((x - xx, y - yy))
def dist((x,y),(xx,yy)):
"""
Distance squared from `(x,y)` to `(xx,yy)`
"""
return norm((x - xx, y - yy))
def setdist2(u,S):
"""
Distance squared from point to set
"""
return min(dist2(u,v) for v in S)
def setdist(u,S):
"""
Distance from point to set
"""
return setdist2(u,S).sqrt()
```

Once you answer some of the questions above, we can investigate the vector-field-plotting part.

```
x, y = var('x y')
VF = plot_vector_field((dx(x,y),dy(x,y)), (x,-4,4), (y,-4,4))
(C + VF).show(aspect_ratio=1)
```

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.