Back to November Calendar            ←Previous Entry                             Next Entry→

November 22, 2005

 

An Example  a Bernstein Glob

 

The goal of this paper is to develop software for drawing Bernstein Globs.  Remember them?

 

The first step for drawing an n-glob is to divide the perimeter of the unit circle into n equal segments by the points (say for 16 points)

 

>> n=16;

>> t = 0:2*pi/n:2*pi

>> circleNodes = [cos(t);sin(t)];

 

   Then we’d like to find the radius and the coordinates of the center of the circle perpendicular to the unit circle at  and .  The radius is relatively simple: rotate the points so that θ1 = 0.  It is then evident (see the diagram at right) that the radius we seek is the length of

AC = tan(θ /2) = tan((θ2  θ1)/2).  The half-angle formulas further yield . 

 

For the coordinates of the center of the perpendicular circle, consider the lines tangent to the circle at  and .  They each have the form  ,  so equating the y coordinates, we have:

                            

We can then substitute for x and solve for y:

                                  

This seems justified by symmetry.  Note also that  …

 

If you feel compelled to verify that the distance from the center of the perpendicular circle  to either of the points of intersection with the unit circle agrees with the aforementioned radius formula for the radius, you could express the square of this radius length with the somewhat ungainly expression:

                   

whose simplification I have, in assignation, left to the Voyage 200 (as shown at right.)

That’s relatively simple!  Indeed:

               

Since a circle centered at (h, k) with radius R can be parameterized by

                                                                 

In this case, we have expanded form for the parameterization of the perpendicular circle

                       

or, if you prefer, the more compact form using addition identities:

                                            

 

Suppose you just want to draw the portion of this circle either the inside or the outside the unit circle, you’d want to know what value of t corresponds to the point on the circle.  This means solving equations like x(t) = cos(θ1):

                      

or

                                            

If you want to use the Voyage 200 to solve this equation, be sure it’s in expanded form.  As the two screenshots below attempt to demonstrate, instead of the simple sin-1(sin(θ1)), you get

                                         

Go figure!

for attempts this, it returns the following (not the “and” conditions are all met):

                            

which suggests .  Note that  while .

 

Similarly, ,  and  

Well, I’ve got some doubts about this now, so I turn to a simple example.  Let θ1 = 0 and   so that

                                                           

For what t does this circle intersect the unit circle?  The coordinates of intersection are (1, 0) and .  Evidently,  will take us to (1,0) whereas  brings us to .  

 

Here’s a little table of values:

t

0

π/4

π/2

3π/4

π

5π/4

3π/2

7π/4

x

 

 

1

 

 

 

1

 

y

 

 

 

 

 

 

0

 

 

…and the Voyage 200 graph in Zoom Square mode:

 

 

Let’s put it to the test in MatLab.

 

Here’re the m-files for the general x and y parametric functions

 

function y = ycoord(a1, a2, t)

% the y-coordinate of a parameterized section of the circle

% which intersects the unit circle orthogonally at polar

y=(cos(a1)-cos(a2)+(1 - cos(a2-a1)).*sin(t))/sin(a2-a1)

 

function x = xcoord(a1, a2, t)

% the xcoordinate of a parameterized section of the circle

% which intersects the unit circle orthogonally at polar

x=(sin(a2)-sin(a1)+(1 - cos(a2-a1)).*cos(t))/sin(a2-a1)

 

First, I plot the unit circle for reference:

 

>> figure

>> hold on

>> theta = 0:pi/36:2*pi;

>> plot(cos(theta),sin(theta))

 

Then I fix the domain for the outnub of the blob and plot it:

 

>> t1=-pi/2:pi/100:3*pi/4;

>> plot(xcoord(0,pi/4,t1),
        ycoord(0,pi/4,t1),'g','LineWidth',2)

 

…and the beginnings of a blob takes shape:

 

Now suppose we want to draw the inside the unit-circle arc (an innub, if you like) of a circle orthogonal to the unit circle at .  Then I want to draw the reflection of the first nub through the y-axis.  Here are the MatLab commands and their result:

 

>> t2=5*pi/4:pi/100:7*pi/4;

>> plot(xcoord(pi/4,3*pi/4,t2), 
        ycoord(pi/4,3*pi/4,t2),'g','LineWidth',2)

>> t3=pi/4:pi/100:3*pi/2;

>> plot(xcoord(3*pi/4,pi,t3),ycoord(3*pi/4,pi,t3),'g','LineWidth',2)

 

 

We continue in this style to draw the closed curve which forms a 3-blob:

>> t4=7*pi/4:pi/100:5*pi/2;

>> plot(xcoord(pi,5*pi/4,t4),
        ycoord(pi,5*pi/4,t4),'g','LineWidth',2)

>> t5=3*pi/4:pi/100:2*pi;

>> plot(xcoord(5*pi/4,3*pi/2,t5),
        ycoord(5*pi/4,3*pi/2,t5),'g','LineWidth',2)

>> t6=pi/2:pi/100:pi;

>> plot(xcoord(3*pi/2,2*pi,t6),
        ycoord(3*pi/2,2*pi,t6),'g','LineWidth',2)

 

 

It’s then pretty straightforward to fit in the unit blobs to complete the 6-glob:

 

>> t7=theta;

>> plot(xcoord(pi/4,pi/2,t7),
        ycoord(pi/4,pi/2,t7),'r','LineWidth',2)

>> plot(xcoord(pi/2,3*pi/4,t7),

        ycoord(pi/2,3*pi/4,t7),'r','LineWidth',2)

>> plot(xcoord(5*pi/4,3*pi/2,t7),

        ycoord(5*pi/4,3*pi/2,t7),'r','LineWidth',2)

>> plot(xcoord(pi/4,pi/2,t7),

        ycoord(pi/4,pi/2,t7),'r','LineWidth',2)

>> plot(xcoord(pi/2,3*pi/4,t7),

        ycoord(pi/2,3*pi/4,t7),'r','LineWidth',2)

>> plot(xcoord(5*pi/4,3*pi/2,t7),

        ycoord(5*pi/4,3*pi/2,t7),'r','LineWidth',2)

>> plot(xcoord(pi/4,pi/2,t7),

        ycoord(pi/4,pi/2,t7),'r','LineWidth',2)

>> plot(xcoord(pi/2,3*pi/4,t7),

        ycoord(pi/2,3*pi/4,t7),'r','LineWidth',2)

>> plot(xcoord(5*pi/4,3*pi/2,t7),

        ycoord(5*pi/4,3*pi/2,t7),'r','LineWidth',2)

>> plot(xcoord(pi/4,pi/2,t7),

        ycoord(pi/4,pi/2,t7),'r','LineWidth',2)

>> plot(xcoord(pi/2,3*pi/4,t7),

        ycoord(pi/2,3*pi/4,t7),'r','LineWidth',2)

>> plot(xcoord(pi,5*pi/4,t7),

        ycoord(pi,5*pi/4,t7),'r','LineWidth',2)

>> plot(xcoord(3*pi/2,7*pi/4,t7),

        ycoord(3*pi/2,7*pi/4,t7),'r','LineWidth',2)

>> plot(xcoord(7*pi/4,2*pi,t7),

        ycoord(7*pi/4,2*pi,t7),'r','LineWidth',2)

 

I took the liberty of filling the blobs to make a more pleasing glob show:

 

 

Miscellaneous notes:

 

To be sure, if cos(t) = sin(θ1) then, substituting, expanding and collecting, shows that the equation is true for almost all unequal angles:

       

 

 

But I’m thinking there’s a better approach here.  Let’s look at the simple case of the circle which intersects the unit circle orthogonally at A(1,0) and at B(cos(θ), sin(θ)). In the diagram below, you can see that the center of this circle is at the intersection of the lines tangent to the unit circle at A and B.  The equations of these tangent lines are x = 1 and , and examining the y-coordinate of C gives a nice proof of the identity, .

The circle is centered at (h, k) =  and has a radius R=  and so the parametric equations are either

                                                             

or

                                                           

Rotating this circle through some multiple of θ can be obtained by multiplying the (x,y) vector by the rotation matrix , so it’s parameterization is

                          

 

Well that doesn’t seem so great…