ni = 200
for i=10 to ni do { x = i/ni*3
findroot( f1, 0, 0.2, 1e-8, y )
P[i-10]: x,y }
np = ni-10
for deg=0 to 359 by 45 do {
for i=0 to np do { Z[i]: rot_(P[i].x,P[i].y,deg*dtor_) }
fitcurve(Z,np,thick 2 outlined "blue")
for i=0 to np do { Z[i]: rot_(P[i].x,-P[i].y,deg*dtor_) }
fitcurve(Z,np,thick 2 outlined "blue")
}