Ellipse Intersect Line

[STxAxTIC]

In response to

https://www.qb64.org/forum/index.php?topic=2132

... all I could think is "why stop at circles when you can do ellipses?"

So this program calculates and plots the intersection points between any ellipse and any line. It's not polished whatsoever - you may break something if you try and tweak it. I didn't spend the time to make this more software-like. That can be on you guys if needed. Here's the hard part though:

Code: QB64: [Select]

1. SCREEN 12
2. PAINT (1, 1), 15
3.  
4. xorig = 0
5. yorig = 50
6.  
7. CALL cline(xorig, yorig, xorig + _WIDTH, yorig, 8)
8. CALL cline(xorig, yorig, xorig + -_WIDTH, yorig, 8)
9. CALL cline(xorig, yorig, xorig, yorig + _HEIGHT, 8)
10. CALL cline(xorig, yorig, xorig, yorig - _HEIGHT, 8)
11.  
12. xzoom = 20
13. yzoom = 20
14.  
15. ' Line
16. b = -2
17. vy = 1
18. vx = -.25
19. v = SQR(vx ^ 2 + vy ^ 2)
20. vx = vx / v
21. vy = vy / v
22. m = vy / vx
23. FOR alpha = -20 TO 20 STEP .001
24.     x = alpha * vx
25.     y = alpha * vy + b
26.     CALL ccircle(xorig + x * xzoom, yorig + y * yzoom, 1, 1)
27. NEXT
28.  
29. ' Ellipse
30. x0 = 2
31. y0 = -2
32. ax = 10
33. ay = 10
34. bx = -5
35. by = 5
36. FOR t = 0 TO 6.284 STEP .001
37.     x = x0 + ax * COS(t) + bx * SIN(t)
38.     y = y0 + ay * COS(t) + by * SIN(t)
39.     CALL ccircle(xorig + x * xzoom, yorig + y * yzoom, 1, 4)
40. NEXT
41.  
42. ' Intersections
43. a2 = ax ^ 2 + ay ^ 2
44. b2 = bx ^ 2 + by ^ 2
45. av = ax * vx + ay * vy
46. bv = bx * vx + by * vy
47. rbx = 0 - x0
48. rby = b - y0
49. adbr = ax * rbx + ay * rby
50. bdbr = bx * rbx + by * rby
51. aa = av ^ 2 / a2 ^ 2 + bv ^ 2 / b2 ^ 2
52. bb = 2 * (av * adbr / a2 ^ 2 + bv * bdbr / b2 ^ 2)
53. cc = adbr ^ 2 / a2 ^ 2 + bdbr ^ 2 / b2 ^ 2 - 1
54. alpha1 = (-bb + SQR(bb ^ 2 - 4 * aa * cc)) / (2 * aa)
55. alpha2 = (-bb - SQR(bb ^ 2 - 4 * aa * cc)) / (2 * aa)
56. x1 = alpha1 * vx
57. x2 = alpha2 * vx
58. y1 = alpha1 * vy + b
59. y2 = alpha2 * vy + b
60. CALL ccircle(xorig + x1 * xzoom, yorig + y1 * yzoom, 10, 1)
61. CALL ccircle(xorig + x2 * xzoom, yorig + y2 * yzoom, 10, 1)
62.  
63. END
64.  
65. SUB cline (x1, y1, x2, y2, col)
66.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
67. END SUB
68.  
69. SUB ccircle (x1, y1, r, col)
70.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), r, col
71. END SUB

[STxAxTIC]

Alright, I made it into "software"... Now you can change everything on the fly. Let me know if any bugs show up. (Note the line is assumed to have infinite length, it's not a ray.)

Code: QB64: [Select]

1. SCREEN 12
2.  
3. xorig = 0
4. yorig = 0
5.  
6. CALL cline(xorig, yorig, xorig + _WIDTH, yorig, 8)
7. CALL cline(xorig, yorig, xorig + -_WIDTH, yorig, 8)
8. CALL cline(xorig, yorig, xorig, yorig + _HEIGHT, 8)
9. CALL cline(xorig, yorig, xorig, yorig - _HEIGHT, 8)
10.  
11. xzoom = 20
12. yzoom = 20
13.  
14. ' Initialize line
15. b = -2
16. d = 2
17. lineang = .1
18. vx = COS(lineang)
19. vy = SIN(lineang)
20. m = vy / vx
21.  
22. ' Initialize ellipse
23. x0 = 2
24. y0 = -2
25. ellipsearg = .2
26. amag = 10
27. ax = amag * COS(ellipsearg)
28. ay = amag * SIN(ellipsearg)
29. bmag = 5
30. bx = bmag * COS(ellipsearg + 3.14 / 2)
31. by = bmag * SIN(ellipsearg + 3.14 / 2)
32.  
33. DO
34.  
35.     DO WHILE _MOUSEINPUT
36.         x = _MOUSEX
37.         y = _MOUSEY
38.         IF ((x > 0) AND (x < _WIDTH) AND (y > 0) AND (y < _HEIGHT)) THEN
39.             IF _MOUSEBUTTON(1) THEN
40.                 x = _MOUSEX
41.                 y = _MOUSEY
42.                 x0 = (x - _WIDTH / 2) / xzoom
43.                 y0 = (-y + _HEIGHT / 2) / yzoom
44.             END IF
45.             IF _MOUSEBUTTON(2) THEN
46.                 x = _MOUSEX
47.                 y = _MOUSEY
48.                 d = (x - _WIDTH / 2) / xzoom
49.                 b = (-y + _HEIGHT / 2) / yzoom
50.             END IF
51.             IF _MOUSEWHEEL > 0 THEN
52.                 lineang = lineang + .01
53.                 vx = COS(lineang)
54.                 vy = SIN(lineang)
55.                 m = vy / vx
56.             END IF
57.             IF _MOUSEWHEEL < 0 THEN
58.                 lineang = lineang - .01
59.                 vx = COS(lineang)
60.                 vy = SIN(lineang)
61.                 m = vy / vx
62.             END IF
63.         END IF
64.     LOOP
65.  
66.     SELECT CASE _KEYHIT
67.         CASE 18432
68.             bmag = bmag + .1
69.             bx = bmag * COS(ellipsearg + 3.14 / 2)
70.             by = bmag * SIN(ellipsearg + 3.14 / 2)
71.         CASE 20480
72.             bmag = bmag - .1
73.             bx = bmag * COS(ellipsearg + 3.14 / 2)
74.             by = bmag * SIN(ellipsearg + 3.14 / 2)
75.         CASE 19200
76.             ellipsearg = ellipsearg - .1
77.             ax = amag * COS(ellipsearg)
78.             ay = amag * SIN(ellipsearg)
79.             bx = bmag * COS(ellipsearg + 3.14 / 2)
80.             by = bmag * SIN(ellipsearg + 3.14 / 2)
81.         CASE 19712
82.             ellipsearg = ellipsearg + .1
83.             ax = amag * COS(ellipsearg)
84.             ay = amag * SIN(ellipsearg)
85.             bx = bmag * COS(ellipsearg + 3.14 / 2)
86.             by = bmag * SIN(ellipsearg + 3.14 / 2)
87.     END SELECT
88.  
89.     ' Intersections
90.     a2 = ax ^ 2 + ay ^ 2
91.     b2 = bx ^ 2 + by ^ 2
92.     av = ax * vx + ay * vy
93.     bv = bx * vx + by * vy
94.     rbx = d - x0
95.     rby = b - y0
96.     adbr = ax * rbx + ay * rby
97.     bdbr = bx * rbx + by * rby
98.     aa = av ^ 2 / a2 ^ 2 + bv ^ 2 / b2 ^ 2
99.     bb = 2 * (av * adbr / a2 ^ 2 + bv * bdbr / b2 ^ 2)
100.     cc = adbr ^ 2 / a2 ^ 2 + bdbr ^ 2 / b2 ^ 2 - 1
101.     arg = bb ^ 2 - 4 * aa * cc
102.     IF (arg > 0) THEN
103.         alpha1 = (-bb + SQR(arg)) / (2 * aa)
104.         alpha2 = (-bb - SQR(arg)) / (2 * aa)
105.         x1 = alpha1 * vx + d
106.         x2 = alpha2 * vx + d
107.         y1 = alpha1 * vy + b
108.         y2 = alpha2 * vy + b
109.     ELSE
110.         x1 = -999
111.         y1 = -999
112.         x2 = -999
113.         y2 = -999
114.     END IF
115.  
116.     GOSUB draweverything
117.  
118.     _LIMIT 60
119.     _DISPLAY
120. LOOP
121.  
122. END
123.  
124. draweverything:
125. CLS
126. PAINT (1, 1), 15
127. COLOR 0, 15
128. LOCATE 1, 1: PRINT "LClick=Move ellipse, RClick=Move line, Scroll=Tilt line, Arrows=Shift ellipse"
129. FOR alpha = -20 TO 20 STEP .001
130.     x = alpha * vx + d
131.     y = alpha * vy + b
132.     CALL ccircle(xorig + x * xzoom, yorig + y * yzoom, 1, 1)
133. NEXT
134. FOR t = 0 TO 6.284 STEP .001
135.     x = x0 + ax * COS(t) + bx * SIN(t)
136.     y = y0 + ay * COS(t) + by * SIN(t)
137.     CALL ccircle(xorig + x * xzoom, yorig + y * yzoom, 1, 4)
138. NEXT
139. CALL ccircle(xorig + x1 * xzoom, yorig + y1 * yzoom, 10, 1)
140. CALL ccircle(xorig + x2 * xzoom, yorig + y2 * yzoom, 10, 1)
141. RETURN
142.  
143. SUB cline (x1, y1, x2, y2, col)
144.     LINE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2)-(_WIDTH / 2 + x2, -y2 + _HEIGHT / 2), col
145. END SUB
146.  
147. SUB ccircle (x1, y1, r, col)
148.     CIRCLE (_WIDTH / 2 + x1, -y1 + _HEIGHT / 2), r, col
149. END SUB

[bplus]

Ah, having done similar, something I can understand and thus appreciate,

A+ man!

[STxAxTIC]

For completeness, attached is the derivation of the math enclosed above.

To re-derive what's going on in bplus's code next door, let a=b=r and simplify.