Topic: Draw ellipse when two tangent points are known

[hmspe]

Does anyone have a routine to create an ellipse based on two points, with the tangent of the ellipse at each point known?  The intent is to do something similar to a fillet, but using an elliptical arc instead of an arc.  The center of the ellipse will lie on the normal from the point that is farthest from the intersection of the tangents.

Thanks.

[ribarm]

This is an old one of mine... See if it can help you - ellipse always have center at point (0 0 0)...

Code: [Select]

(defun c:tan-pt-el ( / o ox oy osm ptt1 ptt2 pt m z xo yo a b pta ptb )
  (setq osm (getvar 'osmode))
  (setvar "osmode" 0)
  (setq o (list 0 0 0))
  (setq ox (list 1 0 0))
  (setq oy (list 0 1 0))
  (command "xline" o ox "")
  (command "xline" o oy "")
  (command "zoom" "c" o "" "")
  (command "regen")
  (prompt "\nEnter point of tangent at x axis (first tangent definition) : ")
  (setq ptt1 (getpoint o))
  (prompt "\nEnter point of tangent at y axis (second tangent definition) : ")
  (setq ptt2 (getpoint o))
  (command "xline" ptt1 ptt2 "")
  (setvar "osmode" 512)
  (setq pt (getpoint "\nPick point on tangent through witch I draw ellipse : "))
  (command "point" pt "")
  (setvar "pdmode" 3)
  (setvar "pdsize" 0)
  (setvar "osmode" 0)
  (setq m (car ptt1))
  (setq z (cadr ptt2))
  (if (minusp m) (setq m (- m)))
  (if (minusp z) (setq z (- z)))
  (if (minusp (car pt)) (setq xo (- (car pt))) (setq xo (car pt)))
  (if (minusp (cadr pt)) (setq yo (- (cadr pt))) (setq yo (cadr pt)))
  (setq b (sqrt (* z yo)))
  (setq a (sqrt (/ (* (expt b 2) m xo) (* z yo))))
  (setq pta (list a 0 0))
  (setq ptb (list 0 b 0))
  (command "ellipse" "c" o pta ptb "")
  (setvar 'osmode osm)
  (princ)
)

M.R.

[hmspe]

Thanks. 

[Stefan]

Here is my (unfinished) version. Works only in 2D and draw a full ellipse.

Code - Auto/Visual Lisp: [Select]

1. (defun C:TEST ( / p1 p11 p2 p22)
2.   (if
3.     (and
4.       (setq p1 (getpoint "\nFirst point: "))
5.       (setq p11 (getpoint p1 "\nDirection: "))
6.       (setq p2 (getpoint "\nSecond point: "))
7.       (setq p22 (getpoint p2 "\nDirection: "))
8.       )
9.     (2pEllipse p1 (mapcar '- p11 p1)
10.                p2 (mapcar '- p22 p2)
11.       )
12.     )
13.   (princ)
14.   )
15.  
16. (defun 2pEllipse (p1 d1 p2 d2 / u1 u2 int p u d n m a b o x y)
17.   (if
18.     (setq u1  (atan (cadr d1) (car d1))
19.           u2  (atan (cadr d2) (car d2))
20.           int (inters
21.                 p1 (polar p1 u1 1.0)
22.                 p2 (polar p2 u2 1.0)
23.                 nil
24.               )
25.     )
26.     (progn
27.       (if
28.         (> (distance int p1) (distance int p2))
29.         (setq p p1 u u1 d d2)
30.         (setq p p2 u u2 d d1 p2 p1)
31.       )
32.       (setq n  (list (list (cos u) (sin u)) (list (- (sin u)) (cos u)))
33.             p2 (mxv n (mapcar '- p2 p))
34.             x  (car p2)
35.             y  (cadr p2)
36.             m  ((lambda (l) (/ (cadr l) (car l))) (mxv n d))
37.             b  (* y (/ (- y (* m x)) (- (* 2 y) (* m x))))
38.             a  (/ (* b x) (sqrt (- (* 2 b y) (* y y))))
39.             o  (mapcar '+ p (mxv (apply 'mapcar (cons 'list n)) (list 0.0 b)))
40.       )
41.       (command "ELLIPSE" "c" "non" o "non" (polar o u a) "non" p)
42.     )
43.   )
44. )
45.  
46. (defun mxv (m v)
47.   (mapcar '(lambda (u) (apply '+ (mapcar '* u v))) m)
48.   )
49.  

[Lee Mac]

Here is my (unfinished) version.

Very good Stefan!

I was initially going to use the property that a line extending from the intersection of the tangents and passing through the midpoint of the line between the two given points will also pass through the ellipse center in order to determine the center location and minor radius:



And then use a rearrangement of the cartesian equation of an ellipse with the minor radius and second point in order to find the major radius.

However, your code is far more succinct.

Though, I'm rather embarassed to confess that, when studying your code, I could not derive the equations that you are using to calculate the minor & major radii after you rotate the ellipse to be orthogonal to the WCS axes:

Code - Auto/Visual Lisp: [Select]

1.            b  (* y (/ (- y (* m x)) (- (* 2 y) (* m x))))
2.            a  (/ (* b x) (sqrt (- (* 2 b y) (* y y))))

I understand that variable m is the tangent of the angle formed by the second tangent direction, but could not derive the above equations geometrically from the diagram.

Could you possibly shed some light?

PS: I should think you could alternatively achieve the rotation of the ellipse by exploiting an appropriate change of coordinates using the trans function, but I do like your use of the rotation matrix.

[hmspe]

Stephan:  Thank you.

Lee:  Thanks for the reminder on finding the center of the ellipse.  Lots of years since my last math class.


Martin

[ribarm]

I've tried differently, but I keep getting undetermined result for x...

x^2*0+x*0+0=0
(this is how equation looks like after ACAD calculate it through lisp)

Code: [Select]

tgA=z/x=l/y=s/p=(x+m)/h (1)
s=p*tgA
z=x*tgA
l=y*tgA
z*h=x*(x+m)
m=(z*h-x^2)/x
y=l/tgA

h^2=x*y                 (2)
h=sqrt(x*y)

(x+m)^2=z*l             (3)

p^2=x^2+h^2             (4)

p*r=h*(x+y)             (5)
r=h*(x+y)/p

r^2=h^2+y^2             (6)

q^2=m^2+(l-h)^2         (7)

s^2=z^2+(x+m)^2         (8)

(p+q)^2=l^2+(x+m)^2     (9)
(p+q)^2=l^2+z*l         (...3)
(p+q)^2=l*(l+z)

s^2=(l+z)^2-(p+q)^2     (10)
s^2=(l+z)^2-l*(l+z)     (...9)
s^2=l^2+2*l*z+z^2-l^2-l*z
s^2=z^2+l*z

(x+y)^2=k^2+k^2-2*k*k*cos(180-2A)     (11)
(x+y)^2=2*k^2-2*k^2*cos(180-2A)
(x+y)^2=2*k^2*(1-cos(180-2A))
(x+y)^2=2*k^2*(1+cos2A)
2*k^2=(x+y)^2/(1+cos2A)
4*k^2=(2*x^2+4*x*y+2*y^2)/(1+cos2A)
4*k^2=(2*x^2+4*x*(l/tgA)+2*(l/tgA)^2)/(1+cos2A)

---------------------------

(s+r)^2+q^2=(2*k)^2         (12)
s^2+2*s*r+r^2+q^2=4*k^2         (...10)
z^2+l*z+2*p*tgA*h*(x+y)/p+h^2+y^2+q^2=4*k^2         (...1...5...6)
x^2*tgA^2+l*x*tgA+2*tgA*h*(x+y)+h^2+y^2+q^2=4*k^2         (...1)
x^2*tgA^2+l*x*tgA+2*tgA*(sqrt(x*y))*(x+y)+x*y+y^2+q^2=4*k^2         (...2)
                  ----------------------- ---
                  + tgA^2*(x+y)^2 - tgA^2*(x+y)^2 = 0
                  ---------------
                 +(tgA*(x+y)+(sqrt(x*y)))^2 - tgA^2*(x+y)^2
x^2*tgA^2+l*x*tgA+(tgA*(x+y)+(sqrt(x*y)))^2-tgA^2*(x+y)^2+y^2+q^2=4*k^2
(tgA*(x+y)+(sqrt(x*y)))^2-x^2*tgA^2-2*x*y*tgA^2-y^2*tgA^2+y^2+q^2+x^2*tgA^2+l*x*tgA=4*k^2
((x*tgA+y*tgA)+(sqrt(x*y)))^2-x^2*tgA^2-2*x*y*tgA^2-y^2*tgA^2+y^2+q^2+x^2*tgA^2+l*x*tgA=4*k^2         (...1...2)
((z+l)+h)^2-z^2-2*h^2*tgA^2-l^2+(l/tgA)^2+q^2+z^2+l*z=4*k^2
           ----                              ----
(z+l+h)^2-2*h^2*tgA^2-l^2+(l/tgA)^2+q^2+l*z=4*k^2
z^2+l^2+h^2+2*z*l+2*z*h+2*l*h-2*h^2*tgA^2-l^2+(l/tgA)^2+q^2+l*z=4*k^2
   ----    ------                        ----              ----
z^2+h^2+3*z*l+2*z*h+2*l*h-2*h^2*tgA^2+(l/tgA)^2+q^2=4*k^2         (...7)
z^2+h^2+3*z*l+2*z*h+2*l*h-2*h^2*tgA^2+(l/tgA)^2+m^2+(l-h)^2=4*k^2
z^2+h^2+3*z*l+2*z*h+2*l*h-2*h^2*tgA^2+(l/tgA)^2+m^2+l^2-2*l*h+h^2=4*k^2
                   ------                              ------
z^2+h^2+3*z*l+2*z*h-2*h^2*tgA^2+(l/tgA)^2+m^2+l^2+h^2=4*k^2
z^2+h^2+3*z*l+2*z*h-2*h^2*tgA^2+(l/tgA)^2+((z*h-x^2)/x)^2+l^2+h^2=4*k^2         (...1)
z^2+h^2+3*z*l+2*z*h-2*h^2*tgA^2+(l/tgA)^2+z^2*h^2/x^2-2*z*h*x^2/x^2+x^4/x^2+l^2+h^2=4*k^2
z^2+h^2+3*z*l+2*z*h-2*h^2*tgA^2+(l/tgA)^2+z^2*h^2/x^2-2*z*h+x^2+l^2+h^2=4*k^2
x^2*tgA^2+x*y+3*x*tgA*l+2*x*tgA*(sqrt(x*y))-2*x*y*tgA^2+(l/tgA)^2+x^2*tgA^2*x*y/x^2-2*x*tgA*(sqrt(x*y))+x^2+l^2+x*y=4*k^2         (...1...2)
x^2*tgA^2+x*(l/tgA)+3*x*tgA*l+2*x*tgA*(sqrt(x*y))-2*x*(l/tgA)*tgA^2+(l/tgA)^2+x^2*tgA^2*x*(l/tgA)/x^2-2*x*tgA*(sqrt(x*y))+x^2+l^2+x*(l/tgA)=4*k^2         (...1)
                             --------------------                                                    --------------------
x^2*tgA^2+x*(l/tgA)+3*x*tgA*l-2*x*(l/tgA)*tgA^2+(l/tgA)^2+x*l*tgA+x^2+l^2+x*(l/tgA)=4*k^2
                   ----------                            --------
x^2*tgA^2+x*(l/tgA)+4*x*tgA*l-2*x*l*tgA+(l/tgA)^2+x^2+l^2+x*(l/tgA)=4*k^2
                   --------------------
x^2*tgA^2+x*(l/tgA)+2*x*tgA*l+(l/tgA)^2+x^2+l^2+x*(l/tgA)=4*k^2
         ----------                            ----------
x^2*tgA^2+2*x*(l/tgA)+2*x*tgA*l+(l/tgA)^2+x^2+l^2=(2*x^2+4*x*(l/tgA)+2*(l/tgA)^2)/(1+cos2A)         (...11)
x^2*tgA^2+x^2*tgA^2*cos2A+2*x*(l/tgA)+2*x*(l/tgA)*cos2A+2*x*tgA*l+2*x*tgA*l*cos2A+(l/tgA)^2+(l/tgA)^2*cos2A+x^2+x^2*cos2A+l^2+l^2*cos2A-2*x^2-4*x*(l/tgA)-2*(l/tgA)^2=0
-------------------------                                                                                  --------------              ------       
x^2*(tgA^2+tgA^2*cos2A+1+cos2A-2)+2*x*(l/tgA)+2*x*(l/tgA)*cos2A+2*x*tgA*l+2*x*tgA*l*cos2A+(l/tgA)^2+(l/tgA)^2*cos2A+l^2+l^2*cos2A-4*x*(l/tgA)-2*(l/tgA)^2=0
                                 --------------------------------------------------------                                        ------------
x^2*(tgA^2+tgA^2*cos2A+1+cos2A-2)+x*(2*(l/tgA)+2*(l/tgA)*cos2A+2*tgA*l+2*tgA*l*cos2A-4*(l/tgA))+(l/tgA)^2+(l/tgA)^2*cos2A+l^2+l^2*cos2A-2*(l/tgA)^2=0
x^2*((1+cos2A)*(tgA^2+1)-2)+x*((1+cos2A)(2*(l/tgA)+2*l*tgA)-4*(l/tgA))+((1+cos2A)*((l/tgA)^2+l^2)-2*(l/tgA)^2)=0
----------------------------------------------------------------------------------------------------------------

    -((1+cos2A)*(2*(l/tgA)+2*l*tgA)-4*(l/tgA))[+-](sqrt(((1+cos2A)(2*(l/tgA)+2*l*tgA)-4*(l/tgA))^2-4*((1+cos2A)*(tgA^2+1)-2)*((1+cos2A)*((l/tgA)^2+l^2)-2*(l/tgA)^2)))
x = ------------------------------------------------------------------------------------------------------------------------------------------------------------------
    2*((1+cos2A)*(tgA^2+1)-2)
   
    -b[+-](sqrt(b^2-4*a*c))
x = -----------------------
    2*a


Code - Auto/Visual Lisp: [Select]

1. (defun tan ( a )
2.   (/ (sin a) (cos a))
3. )
4.  
5. (defun makeline ( p1 p2 )
6.   (entmake
7.     (list
8.       '(0 . "LINE")
9.       (cons 10 p1)
10.       (cons 11 p2)
11.     )
12.   )
13. )
14.  
15. (defun c:test ( / p1 p11 p2 p22 p1v p1n p111 p222 p3 l an a b c x ce p&q p4 )
16.   (if
17.     (and
18.       (setq p1 (getpoint "\nPick first point : "))
19.       (setq p11 (getpoint p1 "\nDirection : "))
20.       (setq p2 (getpoint "\nPick second point : "))
21.       (setq p22 (getpoint p2 "\nDirection : "))
22.       (setq p1v (mapcar '- p11 p1))
23.       (setq p1n (polar '(0.0 0.0 0.0) (+ (/ pi 2.0) (angle '(0.0 0.0 0.0) p1v)) 1.0))
24.       (setq p111 (mapcar '+ p1 p1n))
25.       (setq p222 (mapcar '+ p2 p1v))
26.       (setq p3 (inters p2 p22 p1 p111 nil))
27.       (setq l (distance p2 (inters p1 p111 p2 p222 nil)))
28.       (if (< (angle p1 p111) (angle p2 p3)) (setq an (- (angle p2 p3) (angle p1 p111))) (setq an (- (angle p1 p111) (angle p2 p3))))
29.       (setq b (- (* (+ 1.0 (cos (* 2.0 an))) (+ (* 2.0 (/ l (tan an))) (* 2.0 l (tan an)))) (* 4.0 (/ l (tan an)))))
30.       (setq a (- (* (+ 1.0 (cos (* 2.0 an))) (+ (expt (tan an) 2.0) 1.0)) 2.0))
31.       (setq c (- (* (+ 1.0 (cos (* 2.0 an))) (+ (expt (/ l (tan an)) 2.0) (expt l 2.0))) (* 2.0 (expt (/ l (tan an)) 2.0))))
32.       (setq x (/ (+ (- b) (sqrt (- (expt b 2.0) (* 4.0 a c)))) (* 2.0 a)))
33.       (setq ce (polar (polar p2 (angle p1 p111) x) (angle p2 p222) l))
34.       (setq p&q (sqrt (* l (+ l (* x (tan an))))))
35.       (setq p4 (polar ce (angle p222 p2) p&q))
36.     )
37.     (progn
38.       (makeline p1 p11)
39.       (makeline p2 p22)
40.       (command "_.ellipse" "_C" "_non" ce "_non" p1 "_non" p4)
41.     )
42.   )
43.   (princ)
44. )
45.  

No use, keep getting divide by zero...

At least I tried differently, regards, M.R.

[Stefan]

Stephan:  Thank you.

You're welcome Martin.

Code - Auto/Visual Lisp: [Select]

1.            b  (* y (/ (- y (* m x)) (- (* 2 y) (* m x))))
2.            a  (/ (* b x) (sqrt (- (* 2 b y) (* y y))))

Could you possibly shed some light?

I will try Lee, as much as my English allows me.

[Lee Mac]

Fantastic explanation Stefan, thanks 
- it was the equation for the tangent to the ellipse that was eluding me.

[Lee Mac]

Here is my version, using the observations noted in my earlier post:

Code: [Select]

(defun _ellipse ( p1 d1 p2 d2 / i1 i2 r1 tmp )
    (if (setq i1 (inters p1 (mapcar '+ p1 d1) p2 (mapcar '+ p2 d2) nil))
        (progn
            (if (< (distance p1 i1) (distance p2 i1))
                (setq tmp p1 p1 p2 p2 tmp)
            )
            (setq i2 (inters i1 (mapcar '(lambda ( a b ) (/ (+ a b) 2.0)) p1 p2) p1 (polar p1 (+ (angle p1 i1) (/ pi 2.0)) 1.0) nil)
                  r1 (distance p1 i2)
                  p2 (trans (mapcar '- p2 i2) 0 (mapcar '- p1 i2))
                  p2 (mapcar '* p2 p2)
                  p2 (trans (list (sqrt (/ (car p2) (- 1.0 (/ (caddr p2) (* r1 r1))))) 0.0) (mapcar '- p1 i2) 0)
            )
            (command "_.ellipse" "_C" "_non" i2 "_non" p1 "_non" (mapcar '+ i2 p2))
        )
    )
)
Call with the same arguments as Stefan's function:

Code: [Select]

(defun c:test ( / d1 d2 p1 p2 )
    (if
        (and
            (setq p1 (getpoint "\n1st point: "))
            (setq d1 (getpoint "\n1st direction: " p1))
            (setq p2 (getpoint "\n2nd point: "))
            (setq d2 (getpoint "\n2nd direction: " p2))
        )
        (_ellipse p1 (mapcar '- d1 p1) p2 (mapcar '- d2 p2))
    )
    (princ)
)

[hmspe]

Thanks for all the help.  I'm taking another run at outlining LWPOLYLINEs with arc segments with varying width.  I work with architects who use lwpolylines for things like handrails for stairs.  For my electrical drawings I want to convert those lwpolylines to outlines so that the electrical components are not obscured by background entities.  Autocad and Bricscad apparently use a constant rate of increase in width per unit along the arc segment, similar to a spiral, for the edges of arc segments with varying width.  This causes anomalies in the entity.  My thought is that this is because a spiral centered at (0,0) is never perpendicular where it cross either axis. 

Consider the lwpolyline created by

Code: [Select]

(command "pline" (list 0. 0. 0.) "w" 2.0 2.0 (list 2. 0. 0.) "a" "w" 2.0 4.0 (list 4. 3. 0.) "l" (list 2.3247915 7. 0.) "")As shown in the attachment the resulting lwpolyline is inelegant, the alignment of the arcs in the outline with the adjacent line segments is poor, and the width at the wide end of the arc is incorrect.  The code in PPL.LSP (available on this site) provides a reasonable approximation of the actual outline, anomalies included.  The image on the right is using either Stefan's or Lee's routine to create the ellipses.