fn f_circumcenter p1 p2 p3 =
(
   dp1=[0,0,0]
   dp1.x=(p2.x-p1.x)
   dp1.y=(p2.y-p1.y)
   dp1.z=(p2.z-p1.z)

   dp2=[0,0,0]
   dp2.x=(p3.x-p1.x)
   dp2.y=(p3.y-p1.y)
   dp2.z=(p3.z-p1.z)

   tnx = dp1.y * dp2.z - dp2.y * dp1.z
   tny = dp1.z * dp2.x - dp2.z * dp1.x
   tnz = dp1.x * dp2.y - dp2.x * dp1.y

   sp1=[0,0,0]
   sp1.x=(p1.x+p3.x)/2.0
   sp1.y=(p1.y+p3.y)/2.0
   sp1.z=(p1.z+p3.z)/2.0

   dp1.x=(sp1.x-p1.x)
   dp1.y=(sp1.y-p1.y)
   dp1.z=(sp1.z-p1.z)

   dp2.x=tnx
   dp2.y=tny
   dp2.z=tnz

   nx = dp1.y * dp2.z - dp2.y * dp1.z
   ny = dp1.z * dp2.x - dp2.z * dp1.x
   nz = dp1.x * dp2.y - dp2.x * dp1.y

   sp2=[0,0,0]
   sp2.x=sp1.x+nx
   sp2.y=sp1.y+ny
   sp2.z=sp1.z+nz

   sp3=[0,0,0]
   sp3.x=(p1.x+p2.x)/2.0
   sp3.y=(p1.y+p2.y)/2.0
   sp3.z=(p1.z+p2.z)/2.0

   dp1.x=(sp3.x-p1.x)
   dp1.y=(sp3.y-p1.y)
   dp1.z=(sp3.z-p1.z)

   dp2.x=tnx
   dp2.y=tny
   dp2.z=tnz

   nx = dp1.y * dp2.z - dp2.y * dp1.z
   ny = dp1.z * dp2.x - dp2.z * dp1.x
   nz = dp1.x * dp2.y - dp2.x * dp1.y

   sx4=sp3.x+nx
   sy4=sp3.y+ny
   sz4=sp3.z+nz


   ax=sp2.x-sp1.x
   ay=sp2.y-sp1.y
   az=sp2.z-sp1.z

   bx=sx4-sp3.x
   byy=sy4-sp3.y
   bz=sz4-sp3.z

   cx=sp3.x-sp1.x
   cy=sp3.y-sp1.y
   cz=sp3.z-sp1.z

   qp1=[0,0,0]
   qp1.x = cy * bz - byy * cz
   qp1.y = cz * bx - bz * cx
   qp1.z = cx * byy - bx * cy

   qp2=[0,0,0]
   qp2.x = ay * bz - byy * az
   qp2.y = az * bx - bz * ax
   qp2.z = ax * byy - bx * ay

   dotp=qp1.x*qp2.x+qp1.y*qp2.y+qp1.z*qp2.z

   lee=qp2.x*qp2.x+qp2.y*qp2.y+qp2.z*qp2.z
   lee=sqrt lee

   si=dotp/(lee * lee)

   circumcenter=[0,0,0]
   circumcenter.x=sp1.x+ax*si
   circumcenter.y=sp1.y+ay*si
   circumcenter.z=sp1.z+az*si
   
   return circumcenter
)