MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Sphere through 4 points

  • To: mathgroup at smc.vnet.net
  • Subject: [mg85542] Re: Sphere through 4 points
  • From: dh <dh at metrohm.ch>
  • Date: Wed, 13 Feb 2008 04:27:40 -0500 (EST)
  • References: <fojrcq$gme$1@smc.vnet.net>


Hi Steve,

by reduction of parameters and precomputing expressions you can save 

time. The following takes a list of quadruples of 3-vectors and returns 

a list of centers and radius:



getCircles[v : {{_, _, _, _} ..}] :=

  Module[{x1, x2, x3, y1, y2, y3, z1, z2, z3, v1, v2, v3},

   subs := {x1 -> v1[[1]], x2 -> v2[[1]], x3 -> v3[[1]], y1 -> v1[[2]],

      y2 -> v2[[2]], y3 -> v3[[2]], z1 -> v1[[3]], z2 -> v2[[3]],

     z3 -> v3[[3]]};

   ( {v1, v2,

        v3} = {#[[1]] - #[[4]], #[[2]] - #[[4]], #[[3]] - #[[4]]} ;

      x = (x3^2 (y2 z1 - y1 z2) + x2^2 (-y3 z1 + y1 z3) +

           y2^2 (-y3 z1 + y1 z3) +

           z2 (y3 (x1^2 + y1^2 - y1 y3 + z1 (z1 - z2)) + y1 z2 z3 -

              y1 z3^2) +

           y2 (y3^2 z1 - (x1^2 + y1^2 + z1^2) z3 +

              z1 z3^2))/(2 (x3 y2 z1 - x2 y3 z1 - x3 y1 z2 + x1 y3 z2 +

              x2 y1 z3 - x1 y2 z3)) /. subs;

      y = (x1 x3^2 z2 +

           x3 (y2^2 z1 - (x1^2 + y1^2 + z1^2) z2 + z1 z2^2) +

           x2^2 (x3 z1 - x1 z3) +

           x2 (-(x3^2 + y3^2) z1 + (x1^2 + y1^2 + z1^2) z3 - z1 z3^2) +

            x1 (y3^2 z2 - (y2^2 + z2^2) z3 + z2 z3^2))/(2 (x3 y2 z1 -

             x2 y3 z1 - x3 y1 z2 + x1 y3 z2 + x2 y1 z3 - x1 y2 z3)) /.

        subs;

      z = (x1^2 x3 y2 + x2^2 (-x3 y1 + x1 y3) +

           x3 (y1^2 y2 + y2 z1^2 - y1 (y2^2 + z2^2)) +

           x2 (x3^2 y1 - y3 (x1^2 + y1^2 - y1 y3 + z1^2) + y1 z3^2) -

           x1 (x3^2 y2 - y3 (y2^2 - y2 y3 + z2^2) +

              y2 z3^2))/(2 (x3 y2 z1 - x2 y3 z1 - x3 y1 z2 + x1 y3 z2 +

              x2 y1 z3 - x1 y2 z3)) /. subs;

      r = Sqrt[x^2 + y^2 + z^2];

      {{x, y, z} + #[[4]], r}

      ) & /@ v

   ]



here we compute 1000 circles:

dat=RandomReal[{-1,1},{1000,4,3}];

Timing[getCircles @ dat; ]

{1.188,Null}

hope this helps, Daniel



Steve Gray wrote:

> 	Given 4 points in 3-space not lying in a plane, I want the

> center and radius of the incident sphere. I know about the elegant

> determinant solution (4 or 5 determinants, each 4x4, plus a few

> arithmetic operations) and I'm using that, but maybe someone has coded

> a faster version. 

> 	Thanks for any information.

> 

> Steve Gray

> 




  • Prev by Date: Re: Fitting experimental data
  • Next by Date: Mathematica Book and Documentation
  • Previous by thread: Re: Sphere through 4 points
  • Next by thread: Partial Derivatives in 3D