Find circle radius & circle center in 3D from (3) XYZ pts 5

[DrMetal]

GIVEN: X-Y-Z coordinates of 3 points in 3D space. FIND: Radius of circle on those 3 points & Circle center co-ordinatesBookmark:
Question: Can someone provide a mathematical step by step procedure to calculate and determine a circle radius and the center point coordinates from three GIVEN 3D points with X-Y-Z coordinates?

- I have searched the web extensively, and found lots of proposed solutions, but my college math is quite rusty and after several days, I still can not come up with the correct detailed solution. I can do this easily from three points in 2D, in Excel using three simultaneous equations with TWO variables (i.e. find circle center @ h, j). But I can't get it to work in 3D with a similar spreadsheet approach with THREE variables (i.e. find circle center @ h-j-k). Seems I keep going in mathematical circles (no pun intended) with a modified Excel approach with the three quadratic equations.

- Is it possible to solve this mathematically with (3) simultaneous quadratic equations? Or do I need a 4th equation, determinates, etc?. Seems when I try the double elimination (or substitution) of variables approach, I keep coming up with unsolvable, or meaningless equations. With as much time as we have spent on this problem, we would be happy to compensate someone for their time to solve this for us.

- We have several drawings with bent pipes in 3D where we have only THREE 3D points given to us. We need to be able to determine the TRUE bend radius of the pipes. These are special alloy pipes @ $25k each, so we want to be sure we have the correct 2D bend radius before we bend them. I'd prefer an approach that could be put into a spreadsheet to make reiterative calcs go quickly. I'd prefer to have a detailed step-by-step solution vs. being directed to one of the web solutions. I have seen most of them, and still am stuck on this.
Thanks, MWP

 

[IDS]

I have posted a solution to this problem at:

http://newtonexcelbach.wordpress.com/2009/07/18/two-new-functions-for-ip-xls/

Including a spreadsheet with an on-spreadsheet solution, and two User Defined Functions that will find the 3D circle centre given either any 3 points on the circle, or two tangent points and the intersection of the tangents.

The method used for 3 points on the circle is:
1 Read the 3D coordinates for three points on the circle
2 Translate Point 2 and Point 3 for an origin at Point 1
3 Find polar coordinates of Point 3
4 Rotate Points 2 and 3 about the Z axis so that Point 3 is on the XZ plane
5 Rotate Points 2 and 3 about the Y axis so that Point 3 is on the X axis
6 Find angle of Point 2 from XY plane
7 Rotate Points 2 about the X axis so that Point 2 is on the XY plane
8 Find the XY coordinates of the mid-points of lines 1-2 and 3-2.
9 Find a second point on the perpendiculars through mid-points
10 Find the XY coordinates of the intersection of the perpendiculars. This is the centre of the circle.
11 Find the radius of the circle
12 Rotate and translate the centre point back to the original axes
13 Assign the 3D coordinates of the circle centre, and the circle radius, to the function return value as a 1×4 array


Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[IDS]

I've taken the liberty of modifying Kapitan's routine to work as a function, and added it to the IP spreadsheet, which can be downloaded from the link in the previous post, and I've posted the code below.

The functions give identical answers, and also take almost exactly the same time to run.

Function ArcCenP3_2(Point1 As Variant, Point2 As Variant, Point3 As Variant) As Variant
' Based on "CircumscribedRadius_04() posted by "Kapitan" on the Eng-Tips Forum, 5 Jul 09
' [URL unfurl="true"]http://eng-tips.com/viewthread.cfm?qid=248118&page=1[/URL]
Dim XA As Double, YA As Double, ZA As Double
Dim XB As Double, YB As Double, ZB As Double
Dim XC As Double, YC As Double, ZC As Double
Dim AB As Double, BC As Double, AC As Double
Dim ABi As Double, ABj As Double, ABk As Double
Dim ACi As Double, ACj As Double, ACk As Double
Dim cosBAC  As Double, AD As Double, CD As Double
Dim XD As Double, YD As Double, ZD As Double
Dim CDi As Double, CDj As Double, CDk As Double
Dim Ni As Double, Nj As Double, Nk As Double
Dim sinBAC As Double, R As Double, X2e As Double, Y2e As Double, Z2e As Double
Dim X_Centre As Double, Y_Centre As Double, Z_Centre As Double
Dim Cent(1 To 1, 1 To 4) As Double

    If TypeName(Point1) = "Range" Then Point1 = Point1.Value2
    If TypeName(Point2) = "Range" Then Point2 = Point2.Value2
    If TypeName(Point3) = "Range" Then Point3 = Point3.Value2
 '   ———————————————————————————————————————————————————————————————————
'   Given three points in space (Point1, Point2, Point3)
XA = Point1(1, 1)
YA = Point1(1, 2)
ZA = Point1(1, 3)
XB = Point2(1, 1)
YB = Point2(1, 2)
ZB = Point2(1, 3)
XC = Point3(1, 1)
YC = Point3(1, 2)
ZC = Point3(1, 3)
'   Lengths of AB, AC, AC
 AB = (((XA - XB) ^ 2) + ((YA - YB) ^ 2) + ((ZA - ZB) ^ 2)) ^ 0.5
BC = (((XB - XC) ^ 2) + ((YB - YC) ^ 2) + ((ZB - ZC) ^ 2)) ^ 0.5
AC = (((XA - XC) ^ 2) + ((YA - YC) ^ 2) + ((ZA - ZC) ^ 2)) ^ 0.5
'   Direction cosines of AB(ABi,ABj,ABk)
ABi = (XB - XA) / AB
ABj = (YB - YA) / AB
ABk = (ZB - ZA) / AB
'   Direction cosines of AC(ACi,ACj,ACk)
ACi = (XC - XA) / AC
ACj = (YC - YA) / AC
ACk = (ZC - ZA) / AC
'   Cosine of angle BAC
cosBAC = (AB ^ 2 + AC ^ 2 - BC ^ 2) / (2 * AB * AC)
AD = cosBAC * AC
CD = (AC ^ 2 - AD ^ 2) ^ 0.5
'   Position of point D, which is C projected normally onto AB
XD = XA + (AD * ABi)
YD = YA + (AD * ABj)
ZD = ZA + (AD * ABk)
'   Direction cosines of CD(Cdi,CDj,CDk)
CDi = (XC - XD) / CD
CDj = (YC - YD) / CD
CDk = (ZC - ZD) / CD
'   Direction cosines of normal to AB and CD
'   — to be used for rotations of circle centre
Ni = (ABk * CDj) - (ABj * CDk)
Nj = (ABi * CDk) - (ABk * CDi)
Nk = (ABj * CDi) - (ABi * CDj)
'   # Diameter of circumscribed circle of a triangle is equal to the
'   the length of any side divided by sine of the opposite angle.
'   This is done in a coordinate system where X is colinear with AB, Y is // to CD,
'   and Z is the normal (N) to X and Y, and the origin is point A
'   R = D / 2
sinBAC = (1 - cosBAC ^ 2) ^ 0.5
R = (BC / sinBAC) / 2
'   Centre of circumscribed circle is point E
X2e = AB / 2
Y2e = (R ^ 2 - X2e ^ 2) ^ 0.5
Z2e = 0
'   Transform matrix'                    Rotations                 Translations
'             ——————————————————————————————————————————————
'               ABi  ,   ABj  ,  ABk                 Xa
'               CDi  ,   CDj  ,  CDk                 Ya
'                Ni  ,    Nj  ,   Nk                 Za
'             ——————————————————————————————————————————————
'   Position of circle centre in absolute axis system
X_Centre = XA + (X2e * ABi) + (Y2e * CDi) + (Z2e * Ni)
Y_Centre = YA + (X2e * ABj) + (Y2e * CDj) + (Z2e * Nj)
Z_Centre = ZA + (X2e * ABk) + (Y2e * CDk) + (Z2e * Nk)

'   ———————————————————————————————————————————————————————————————————

Cent(1, 1) = X_Centre
Cent(1, 2) = Y_Centre
Cent(1, 3) = Z_Centre
Cent(1, 4) = R

ArcCenP3_2 = Cent
End Function

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[Kapitan]

An excellent solution, Doug — and many thanks for taking the time to update my old BASIC routine.

 

[rb1957]

IDS,

maybe i'm being unusually thick, and i haven't gone through your solution line-by-line, however, if i have three points on a plane i can make to different arcs through them depending on the tangency condition. i'm surprised that a math procedure can determine the right choice ?

 

[zekeman]

They are still doing the OP problem of circumscribing a circle, not the tangency problem, albeit too inefficiently.

 

[rb1957]

IDS's earlier post says "will find the 3D circle centre given either any 3 points on the circle, or two tangent points and the intersection of the tangents."

 

[IDS]

rb1957 - no, 3 points define a triangle on a plane (unless they are co-linear), and you can only draw one circle through the vertices of a triangle.

And yes, the spreadsheet has a function to find the circle defined by two tangent points and their intersection point (in fact it works with one tangent point and any point on the other tangent line past the tangent point). Thanks for reading

zekeman - as noted above my functions do address the two tangent question. If you read the posts you would also have seen that the calculation time using the approach of rotating the points onto a plane was almost identical to that using kapitan's routine, which used the direction cosine approach. To my mind the rotation approach makes it more intuitively obvious what is going on and why it works. I might have a go at putting your method into a function when I have some time, but not today.

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[williedawg]

After some fumbling around, I was able to get Kapitan's Code to work (This is the first time I've used VB)

Just as he suggested, it was pasted into a Macro in VB and then Run. In order to use some other 3D points in space, I invoked the Macro, then used "Step Into" to edit in some different values; Ran it, noted the new values, then went back and compared to a SW drawing. All looked OK. Is this the "normal" way to do this?

Now,... how do I get the 2nd Code written by Doug in this thread to work ?
It starts out as a Function, not a Sub. How is this implemented or made into a Macro in Excel Basic ?
Or is some other Excel method/procedure needed?

I looked at the IP spreadsheet and downloaded it, but am having a hard time separating out only what's needed just for the circle and how to initialize it with the variable 3 point inputs.

Thanks for all the time and effort you have put into this.

 

[GregLocock]

No that is a truly horrible way to use a user-function!

It should work exactly like a built in function. The 'help' on user functions describes what to do.

Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[IDS]

williedawg - There are two ways to use a VBA function, either as a function in the spreadsheet, just like the Excel built-in functions, or call it from a VBA subroutine.

You can either follow the instructions below using the IP.xls spreadsheet, or if you prefer you can copy and paste the code into a Visual Basic module in a new spreadsheet, but if you do that you will need to also copy the functions Rotate() and IP(), which are called from the other functions. It's probably easier just to use the IP spreadsheet for a start.

To use any of the three functions enter the X, Y, and Z coordinates of three points on three separate lines. You can use any sheet. Suppose you enter the data in the ranges: A1:C1, A2:C2, and A3:C3, with row 1 and 3 being points on the circle and line 2 either an intermediate point on the circle or the tangent intersection point. Then to use the functions enter in any convenient cells:

=ArcCenT2IP(A1:C1, A3:C3, A2:C2)
=ArcCenP3(A1:C1, A2:C2, A3:C3)
=ArcCenP3_2(A1:C1, A2:C2, A3:C3)

Note the different order for ArcCenT2IP

In each cell the X coordinate of the circle centre will appear. ArcCenT2IP will give a different result to the other two of course.

The functions return an array; 1x4 for ArcCenP3 and ArcCenP3_2 and 2x4 for ArcCenT2IP. To get the other results select 4 adjacent columns and 1 or 2 rows, with the function entered in the top left cell, press F2 to enter edit mode, then press ctrl-shift-enter. There are examples of this on the IP spreadsheet.

An example of using the function called from another subroutine is given below:

Sub ArcCen()
Dim TP1 As Variant, TP2 As Variant, TIP As Variant
Dim Cent As Variant, Res As Range

TP1 = Range("A1:C1").Value2
TP2 = Range("A2:C2").Value2
TIP = Range("A3:C3").Value2

Set Res = Range("A5:D6")

Cent = ArcCenT2IP(TP1, TP2, TIP)

Res = Cent

End Sub

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[DrMetal]

Hi Doug,
- See what I got you into! Heck-of-a thread on this subject. I had no idea that it would be this "intriguing" to all the posters. Thanks for all of your input on this, and your emails too, and everyone else for the postings. We are getting ready to do the original PIPE developments and BENDS as described in earlier posts. FYI ... these are 8" x XXS (0.875"THK) 2205 DUPLEX SSTL PIPES x 20'L to bend and put into these 3D positions. We'll be using this thread to finally figure out and mathematically verify the correct 2D bends. Will post the final results of our success (or degree of success) later. Thanks again to all. Welcome to keep commenting if there is a need. Mike P.

 

[williedawg]

(Phew)
Got the IP.xls to work.
Entered some values in a new Sheet and compared to SW drawing...accurate to 4 plcs.

"To get the other results select 4 adjacent columns and 1 or 2 rows, with the function entered in the top left cell, press F2 to enter edit mode, then press ctrl-shift-enter."

This part didnt seem to work for me. (maybe something I don't understand about the array function)

Instead of "select 4 adjacent columns", I selected the cell containing =ArcCenT2IP(A1:C1, A3:C3, A2:C2), then drag-selected 3 more columns to the right, but only 1 row, then did the F2 ctrl-shift-enter thing, and it filled the cells with the correct values for a 2pt circle with tangent lines intersecting at the 3rd pt.

For the next row down, I did the same thing, and it filled in the values for the 3pt circle situation

Doug, Kapitan and Greg, Thanks so much for the help on this topic and your patience with an Excel novice.

 

[IDS]

I selected the cell containing =ArcCenT2IP(A1:C1, A3:C3, A2:C2), then drag-selected 3 more columns to the right, but only 1 row, then did the F2 ctrl-shift-enter thing, and it filled the cells with the correct values for a 2pt circle with tangent lines intersecting at the 3rd pt.

That's what I meant. I usually hold down the shift key and use the arrow keys to select. So I start on the cell with the function in it, hold down shift and press -> 3 times then F2 and ctrl-shift-enter. The mouse works just as well of course. I need to find a way to describe the process more clearly.

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[DrMetal]

Hi Doug,
- Quick question ... I must be doing something wrong ... when I plug in new numbers into the points for XYZ in your spreadsheet, the results all return the "#NAME" error, as if it can't read my inputs. Do I need to format the numbers that I enter a certain way? Or something else simple that I am doing wrong? I was just playing around with it a little tonight and could not get the Excel sheet to show numbers. As soon as I change even one of your sample entries, I get nothing but cells of errors. Sorry for the basic question, but I am still an Excel novice, and remember very little about basic and almost nothing about VB for that matter. Thanks again for the work. Mike P.

 

[zekeman]

IDS,
Sorry, I didn't read your post but I don't like rotating and unrotating stuff which to me is unnecessary.

My take is:
The center of the circumscribing circle is simply at the intersection of 3 planes, namely
1) the plane determined by the 3 pts 1,2,3 in space
2) the plane perpendicular to line 1-2 containing the midpoint of 1-2
3) the plane perpendicular to line 2-3 containing the midpoint of 2-3.

First write the equation of the plane (pts1,2,3)
eq1 x+ay+ bz=c
where a,b and c are determined by entering the coordinates of the 3 points, giving 3 equations whose solution of a,b,c determines the plane of eq1.(from my previous post)

Now we get the direction cosines of the 2 lines 1-2 and 2-3
which are proportional to

line 1-2
(x1-x2),(y1-y2),z1-z2)

line 2-3
(x2-x3),y2-y3),z2-z3)

Next we write the equation of the plane perpendicular to line 1-2
eq2 (x1-x2)X+(y1-y2)Y+(z1-z2)Z=D1
where D1 is determined from the midpoint location (x1+x2)/2,y1+y2)/2,(z1+z2)/2 in this plane; so substituting for X,Y,Z I get D1 to be
[(x1^2-x2^2+y1^2-y2^2+z1^2-z2^2]/2=D1
Similarly the plane perpendicular to 2-3 containing the midpoint of 2-3 is
eq3)(x2-x3)X+(y2-y3)Y+(z2-z3)=D2
and
D2=[x2^-x3^2+y2^2-y3^2+z2^2-z3^2]/2
So the problem boils down to solving the 3x3 for for a,b,c, then evaluating D1 and D2 and finally solving for X.Y, Z from the 3x3 set of the 3 planes, which is formally more concise.
But I have no problem for those who need more clarity and don't feel comfortable working in space. and as far a time to do it, it's not that much greater, but why all the additional steps.

 

[IDS]

DrMetal - it sounds like you don't have macros enabled, which is easy to do in Excel 2007 because that is the default state. In 2007 when you start up you get a small message above the spreadsheet saying:

Security Warning Macros have been disabled [Options]

You need to click on the Options box and select "Trust this content"

In earlier versions the option to enable macros is much more obvious.

Pleae let me know what version you are using, and if enabling macros fixes the problem.

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[Kapitan]

I've modified my Basic subroutine to deal with the intersection point case — by a different method, and what I think is far simpler than before. It's all done in 3D using direction cosines: it uses no auxiliary axis system, no transform matrix, no rotations, and no plane equations.

The centre of the required bend radius lays somewhere along the bisector of the two tangent lines that create the point (B) of intersection of the tangents — and this distance, from point B, of this radius centre along this bisector, is the hypotenuse of a right-angled triangle (either ABD or CBD) — depending which is the shorter of the two lines. By finding the bisector of lines AB and BC between A and C, the direction cosines of the bisector are calculated, This distance from B to the radius centre is the length of the line from the closest given point (A or C) to the intersection point, divided by the cosine of half the angle ABC.

To avoid using the other of the two bisectors of the two lines, I've found the bisector line using the mid-point of two points spaced equal distances along the lines AB and BC (I've used an arbitary 1000 units, which is preferably longer than AB or BC.)

The position of the bend radius centre is found by multiplying the direction cosines of the bisector by the length of the hypotenuse - and adding them to the x,y,z values of the intersection point (B).

As before, it will run in an Excel VBA macro: I've tried to keep the thing 'readable' . I haven't (yet) done a check for the three point co-linear case, although it will flip through the straight-line situation. Once again, I've had someone check it against a 3D model — and it works — whether the next pipe bend is a LH or RH one. It also outputs identical results to Doug's method.

Sub CornerRadiusCentre_02()
        Range("A1:Z512").NumberFormat = "###0.000000"
'   ———————————————————————————————————————————————————————————————————
'   Given three points in space (A,B,C)
         Xa = -71: Ya = -37: Za = 47: '   Start Point
         Xb = 59: Yb = -40: Zb = 0 '      Tangent Intersection Point
         Xc = 100: Yc = 7: Zc = 88: '     Finish Point
'   Lengths of AB, AC, AC
         AB = (((Xa - Xb) ^ 2) + ((Ya - Yb) ^ 2) + ((Za - Zb) ^ 2)) ^ 0.5
         BC = (((Xb - Xc) ^ 2) + ((Yb - Yc) ^ 2) + ((Zb - Zc) ^ 2)) ^ 0.5
         AC = (((Xa - Xc) ^ 2) + ((Ya - Yc) ^ 2) + ((Za - Zc) ^ 2)) ^ 0.5
'   Direction cosines of AB
        ABi = (Xa - Xb) / AB
        ABj = (Ya - Yb) / AB
        ABk = (Za - Zb) / AB
'   Direction cosines of BC
        BCi = (Xc - Xb) / BC
        BCj = (Yc - Yb) / BC
        BCk = (Zc - Zb) / BC
'   From point B create two points (G and H) 1000 units along AB and BC
         Xg = Xb + (1000 * ABi)
         Yg = Yb + (1000 * ABj)
         Zg = Zb + (1000 * ABk)
'
         Xh = Xb + (1000 * BCi)
         Yh = Yb + (1000 * BCj)
         Zh = Zb + (1000 * BCk)
'   Find mid-point (M) of G and H
         Xm = (Xg + Xh) / 2
         Ym = (Yg + Yh) / 2
         Zm = (Zg + Zh) / 2
'   Length of line BM (Tangents bisector)
         BM = ((Xm - Xb) ^ 2 + (Ym - Yb) ^ 2 + (Zm - Zb) ^ 2) ^ 0.5
'   Direction cosines of BM
        BMi = (Xm - Xb) / BM
        BMj = (Ym - Yb) / BM
        BMk = (Zm - Zb) / BM
'   Cosine of angle ABM (Half angle ABC)
     cosABM = (ABi * BMi) + (ABj * BMj) + (ABk * BMk)
'   Find the shorter of AB and BC
        If AB < BC Then Base = AB Else Base = BC
'   Right-angled triange ABD or CBD
 Hypotenuse = Base / cosABM
'   Find coordinates of Bend Radius Centre
   X_centre = Xb + (Hypotenuse * BMi)
   Y_centre = Yb + (Hypotenuse * BMj)
   Z_centre = Zb + (Hypotenuse * BMk)
'   Find Bend Radius
   Bend_Rad = (Hypotenuse ^ 2 - Base ^ 2) ^ 0.5
'   Results:—
                    Cells(1, 1) = "Bend Radius": Cells(1, 2) = Bend_Rad
                    Cells(2, 1) = "X_centre   ": Cells(2, 2) = X_centre
                    Cells(3, 1) = "Y_centre   ": Cells(3, 2) = Y_centre
                    Cells(4, 1) = "Z_centre   ": Cells(4, 2) = Z_centre
                    Cells(6, 1) = "Given Points"
                    Cells(7, 1) = "Point A(x,y,z)": Cells(7, 2) = Xa
                    Cells(7, 3) = Ya: Cells(7, 4) = Za
                    Cells(8, 1) = "Point B(x,y,z)": Cells(8, 2) = Xb
                    Cells(8, 3) = Yb: Cells(8, 4) = Zb
                    Cells(9, 1) = "Point C(x,y,z)": Cells(9, 2) = Xc
                    Cells(9, 3) = Yc: Cells(9, 4) = Zc
             
End Sub

 

[rb1957]

from the direction cosines of AB and BC,
and the co-ords of A and C,
you know the line thru A (and C) normal to AB (and BC)
these two lines intersect at the centre (O),
and the radius is the distance OC.

 

[IDS]

I have also had a go at simplifying the whole thing, and done away with all the rotations. The code below makes use of the fact that there are three similar right angled triangles that allow the distance from the tangent intersection point to the arc centre to be found from the known coordinates by simple ratios.

More than half the code in the new function is to check that the two tangent points are the same distance from the IP, and if necessary moving the further one so that they are.

The new function runs about 7 times quicker than the earlier one.

Function ArcCenT2IP2(TP1 As Variant, TP2 As Variant, IntP As Variant) As Variant

' Find the centre and radius of a circular arc defined by two tangent points
' and the intersection point of the tangents

    Dim DL1(1 To 4) As Double, DL2(1 To 4) As Double, DL3(1 To 3) As Double, DL4(1 To 3) As Double
    Dim MP(1 To 3) As Double, a As Double, b As Double, c As Double, d As Double
    Dim NDims As Long, i As Long
    Dim Arccen(1 To 2, 1 To 4) As Double
    Const MaxDiff As Double = 0.0000000001

    If TypeName(TP1) = "Range" Then TP1 = TP1.Value2
    If TypeName(TP2) = "Range" Then TP2 = TP2.Value2
    If TypeName(IntP) = "Range" Then IntP = IntP.Value2

    NDims = UBound(TP1, 2)

    ' Check that tangent points are equidistant from IntP
    For i = 1 To NDims
        DL1(i) = IntP(1, i) - TP1(1, i)
        DL2(i) = IntP(1, i) - TP2(1, i)
    Next i
    DL1(4) = (DL1(1) ^ 2 + DL1(2) ^ 2 + DL1(3) ^ 2) ^ 0.5
    DL2(4) = (DL2(1) ^ 2 + DL2(2) ^ 2 + DL2(3) ^ 2) ^ 0.5

    ' If not, scale the longer line towards the IP
    If Abs(DL2(4) / DL1(4) - 1) > MaxDiff Then
        If DL2(4) > DL1(4) Then
            ScaleFact = DL1(4) / DL2(4)
            For i = 1 To NDims
                DL2(i) = DL2(i) * ScaleFact
                TP2(1, i) = IntP(1, i) - DL2(i)
                Arccen(2, i) = TP2(1, i)
            Next i
            DL2(i) = DL2(i) * ScaleFact
            Arccen(2, 4) = 2
        Else
            ScaleFact = DL2(4) / DL1(4)
            For i = 1 To NDims
                DL1(i) = DL1(i) * ScaleFact
                TP1(1, i) = IntP(1, i) - DL1(i)
                Arccen(2, i) = TP1(1, i)
            Next i
            DL1(i) = DL1(i) * ScaleFact
            Arccen(2, 4) = 1
        End If
    End If

' Point A = TP1, B = IP, C = mid-point of TP1-TP2
'Find side lengths of triangle ABC, a, b, c
'Find coordinates of Point C, midpoint of TP1-TP2

For i = 1 To NDims
MP(i) = (TP1(1, i) + TP2(1, i)) / 2
DL3(i) = (MP(i) - IntP(1, i))
DL4(i) = (MP(i) - TP1(1, i))
Next i

a = (DL3(1) ^ 2 + DL3(2) ^ 2 + DL3(3) ^ 2) ^ 0.5
b = (DL4(1) ^ 2 + DL4(2) ^ 2 + DL4(3) ^ 2) ^ 0.5
c = DL1(4)

' Radius = bc / a
' d = (distance arc centre to IP)/a = Rc / ba
R = b * c / a
d = R * c / (b * a)

For i = 1 To NDims
Arccen(1, i) = IntP(1, i) + DL3(i) * d
Next i
Arccen(1, 4) = R

    ArcCenT2IP2 = Arccen
End Function

Doug Jenkins
Interactive Design Services

http://newtonexcelbach.wordpress.com/

 

[rb1957]

"check that the two tangent points are the same distance from the IP, and if necessary moving the further one so that they are."

well that's pretty redundant, 'cause they have to be, congruent triangles.