a = -1; b=1; c=a; d=b; n=40; 

h=(b-a)/n;   % Domain [a b] [c d]

for i=1:n+1,
  x(i) = a + (i-1)*h;  
end
y=x; z=zeros(n+1,n+1);

for i=1:n+1
  for j=1:n+1
    phi1 = sqrt( (x(i)+0.3)^2/0.25^2 + (y(j)+0.4)^2/0.1^2) -1;
    phi2 = sqrt( (x(i)-0.31)^2/0.15^2 + (y(j)-0.4)^2/0.2^2) -1;
    phi(i,j) = min(phi1,phi2);
  end
end 

figure(1); contour(x,y,phi,z,'b'); axis('square')
figure(2); mesh(x,y,phi)
figure(3); contour(x,y,phi,30)


% k-stars:  

r0 = 0.5; k=5; e = 0.15;

for i=1:n+1
  for j=1:n+1
    if abs(x(i)) < 1.0e-13
	if y(j) >= 0.0
	   theta = pi/2.0;
        else
           theta = -pi/2.0;
        end
    else
        theta = atan(y(j)/x(i));
        if x(i) <= 0.0
           theta = pi + theta;
	end
    end

    r = r0 + e*sin(k*theta);
    r1 = sqrt(x(i)*x(i) + y(j)*y(j) );
    phi(i,j) = r1 - r;
  end
end

figure(4); contour(x,y,phi,z,'b'); axis('square')