Hi Rune and Peter

Peter Popov wrote:

> On Sun, 13 Feb 2000 14:19:32 +0100, "Rune" <run### [at] inamecom>
> wrote:
>
> >Yes, I too have found that that's the only way.
> >Now I just need to know how to calculate all the
> >points of intersection. How do I do that?
> >
> >Greetings,
> >
> >Rune
>
> For each two circles, solve the system:
>
> (X-X1)^2+(Y-Y1)^2=R1^2
> (X-X2)^2+(Y-Y2)^2=R2^2
>
> where the two circles have centers (X1,Y1) and (X2,Y2) and radii R1
> and R2, respectively. The system will have no real roots if they do
> not intersect, one double real root if they touch, two distinct real
> roots if they intersect and infinitely many solutions if they
> coincide.
>
> Peter Popov
> pet### [at] tagpovrayorg
> ICQ: 15002700

I just had to see if I could solve the problem of finding these points
of
intersections without having to solve the above equation system.
(or maybe I did solve it without knowing it ? ;-)

Here is my graphical (and numerical) POV-Ray solution.

It can be used to solve the problem for intersection of spheres
(3D-problems)
or intersection of circles (2D-problems).

Best regards,
Tor Olav Kristensen

mailto:tor### [at] hotmailcom
http://www.crosswinds.net/~tok/tokrays.html

// ===================================================================
// Finding the points where two spheres (or circles) meet. By Tor Olav
Kristensen.
// ===================================================================

#version 3.1;
#include "colors.inc"

// ===================================================================

// Data for the spheres (3D) or circles (2D)
// The z co-ordinate in the sphere centres can be set to zero for a
// 2D-problem/solution in the XY-plane.
// If such an 2D solution is needed: Also select UpVector = z
//
// Note: If you choose the two spheres to have the same centre, then
// my solution will not work. If they also have the same radius, there
// will have infinitely many "intersections", else they will have none.
// It would however be simple to test for these cases before my code
// is used.

#declare Radius1 = 3.5;
//#declare Center1 = <-3, -2, 0>; // 2D-problem
#declare Center1 = <-3, -2, 2>; // 3D-problem

#declare Radius2 = 3;
//#declare Center2 = <1, -5, 0>; // 2D-problem
#declare Center2 = <1, -5, 1>; // 3D-problem

// ===================================================================

// My attempt to solve the problem

#declare Vector1to2 = Center2 - Center1;
#declare Distance = vlength(Vector1to2);
#declare C = (Radius1*Radius1 - Radius2*Radius2 + Distance*Distance)/
             (2*Distance);
#declare BetweenPoint = Center1 + C*vnormalize(Vector1to2);
#declare TR = Radius1*Radius1 - C*C;
//#declare UpVector = z; // for 2D-problem
#declare UpVector = <1, 3, -2>; // example for 3D-problem
// Note: Do not choose an UpVector that is parallel to Vector1to2
#declare SideVector = vcross(Vector1to2, UpVector);

#if (TR > 0)
  #declare SideVector = sqrt(TR)*vnormalize(SideVector) // Resizing
  #declare InterSection1 = BetweenPoint + SideVector;
  #declare InterSection2 = BetweenPoint - SideVector;
#else
  #if (TR = 0)
    #declare InterSection = BetweenPoint;
  #end
#end

#declare NewUpVector = vcross(SideVector, Vector1to2);
// I use this NewUpVector instead of UpVector just to make sure
// that I have an "Up"-vector that is perpendicular to the line
// between the two spherecentres.
// Should I have used vcross(Vector1to2, SideVector) instead?
// Maybe not: The numerical results seems to be the same.
// But the graphical view changes to "the other side"

// ===================================================================

// Graphical illustration of the problem

sphere { Center1, Radius1 pigment { color Green } }
sphere { Center2, Radius2 pigment { color Cyan  } }

plane {
  NewUpVector, 0
  translate BetweenPoint
  pigment { color Yellow }
}

// ===================================================================

// Graphical illustration of my proposed solution

#if (TR > 0)
  sphere { InterSection1, 0.8 pigment { color Red  } }
  sphere { InterSection2, 0.8 pigment { color Blue } }
#else
  #if (TR = 0)
    sphere { InterSection, 0.8 pigment { color Magenta } }
  #end
#end

// Note: A more correct graphical solution for a the 3D-problem
// would be to put a torus along the "circle of intersection" of
// the two spheres. Its major radius would then equal sqrt(TR)
// It's center should be at BetweenPoint. You can find an image here:
// http://www.crosswinds.net/~tok/Sphere2Sphere04.jpg

// ===================================================================

// Text output of the proposed solution to the debug window
// Press Alt M when parsing to view this text output

#macro PrintVector(Vector)

  #debug concat("<", str(Vector.x, 0, -1), ", ",
                     str(Vector.y, 0, -1), ", ",
                     str(Vector.z, 0, -1), " >")

#end // macro PrintVector


#if (TR > 0)
  #debug concat("\nFirst  intersection: ")
  PrintVector(InterSection1)
  #debug concat("\nSecond intersection: ")
  PrintVector(InterSection2)
#else
  #if (TR = 0)
    #debug concat("\nIntersection: ")
    PrintVector(InterSection)
  #else
    #debug "\nThere are no intersections"
  #end // if
#end // if

#debug "\n"

// ===================================================================

// Just my macros for axes
// X-axis is red
// Y-axis is green
// Z-axis is blue

#macro Arrow(StartPoint, EndPoint, Radius)

  #local a = Radius*3;
  #local b = Radius*2;
  #local c = Radius;
  #local Dirv = vnormalize(EndPoint-StartPoint);

  merge {
    difference {
      cone {
        EndPoint, 0,
        EndPoint-(a+b)*Dirv, Radius+c
      }
      cone {
        EndPoint-a*Dirv, Radius,
        EndPoint-(a+2*b)*Dirv, Radius+2*c
      }
    }
    cylinder {
      StartPoint, EndPoint-a*0.999*Dirv, Radius
    }
  }

#end // macro Arrow


#macro Axes(Size)

  #local Thickness = Size/40;

  union {
    object {
      Arrow(-Size*x, Size*x, Thickness)
      pigment { color Red }
    }
    object {
      Arrow(-Size*y, Size*y, Thickness)
      pigment { color Green }
    }
    object {
      Arrow(-Size*z, Size*z, Thickness)
      pigment { color Blue }
    }
  }

#end // macro Axes


Axes(5)

// ===================================================================

#declare CameraPosition = BetweenPoint + 3*(Radius1+Radius2)*
                          vnormalize(NewUpVector);

#if (TR > 0)
  light_source { CameraPosition + (Radius1+Radius2)*SideVector, White }
  light_source { CameraPosition - (Radius1+Radius2)*SideVector, White }
#else
  light_source { CameraPosition, White }
  light_source { CameraPosition, White }
#end

camera {
  location CameraPosition
  look_at BetweenPoint
}

// ===================================================================

Post a reply to this message