1 %{
2 Finds closest point on a ellipse, for two cases,
3 to a given point.
4 Case 1: Ellipse is a circle. Point is outside or
5 inside of ellipse.
6 Case 2: Ellipse centered at origin with major
7 and minor axes the x and y axes, and given point
8 outside the ellipse.
9 Octave 5.1.0
10 John Bryan 2019
11 %}
12
13 echo off;
14 close all;
15 clear all;
16 clear;clf
17 graphics_toolkit fltk
18
19 function Case2Plot(x0,y0,a,b,c,d,Cx,Cy,m,filename)
20 figure (m);
21 x0=0; % horizontal center coordinate
22 y0=0; % vertical center coordinate
23 t=-pi:0.01:pi;
24 x=x0+a*cos(t);
25 y=y0+b*sin(t);
26 plot(x,y)
27 hold on
28 plot(Cx,Cy,"*","markersize",5,"markerfacecolor","auto");
29 plot(c,d,"*","markersize",5,"markerfacecolor","auto");
30 text(c,d,['(' num2str(c,"%3.2f") ',' num2str(d,"%3.2f") ')'])
31 text(Cx,Cy,['(' num2str(Cx,"%3.2f") ',' num2str(Cy,"%3.2f") ')'])
32 xmin=findmin(-a,c)-2;
33 xmax=findmax(a,c)+2;
34 ymin=findmin(-b,d)-2;
35 ymax=findmax(b,d)+2;
36 axis([xmin xmax ymin ymax]);
37 set(gcf,'InvertHardCopy','off')
38 print (gcf, "-S500,500", filename, '-dpng')
39 endfunction
40
41 function Case1Plot(Ax,Ay,Bx,By,Cx,Cy,r,N,i,filename)
42 figure (i);
43 theta=2*pi*linspace(0,1,200);
44 xc=Ax+r*cos(theta);
45 yc=Ay+r*sin(theta);
46 plot([Ax Bx],[Ay By])
47 hold on;
48 plot(Ax,Ay,"*","markersize",7,"markerfacecolor","auto");
49 plot(Bx,By,"*","markersize",7,"markerfacecolor","auto");
50 plot(Cx,Cy,"*","markersize",7,"markerfacecolor","auto");
51 h=plot(xc,yc,'linewidth',1);
52 plot([Ax Bx],[Ay By])
53 text(Ax,Ay,['(' num2str(Ax) ',' num2str(Ay) ')'])
54 text(Bx,By,['(' num2str(Bx) ',' num2str(By) ')'])
55 text(Cx,Cy,['(' num2str(Cx,"%3.2f") ',' num2str(Cy,"%3.2f") ')'])
56 xmin=findmin(Ax-r,Bx)-1;
57 xmax=findmax(Ax+r,Bx)+1;
58 ymin=findmin(Ay-r,By)-1;
59 ymax=findmax(Ay+r,By)+1;
60 axis([xmin xmax ymin ymax]);
61 axis equal;
62 set(gcf,'InvertHardCopy','off')
63 print (gcf, "-S500,500", filename, '-dpng')
64 endfunction
65
66 function xin = CalculateInitialGuess(c,d,a,b)
67 if (a>b) r=a; else r=b; endif;
68 [Cx,Cy]=CalculateCirclePoint(0,0,c,d,r);
69 xin=Cx;
70 return;
71 endfunction
72
73 function [Ex,Ey]=CalculateCirclePoint(Dx,Dy,Px,Py,r)
74 magPD=sqrt((Px-Dx)^2+(Py-Dy)^2)
75 if ((Px-Dx)==0) Ex=Dx;
76 else Ex=Dx+(r*(Px-Dx)/magPD); endif;
77 if ((Py-Dy)==0) Ey=Dy;
78 else Ey=Dy+(r*(Py-Dy)/magPD); endif;
79 endfunction
80
81 function max0 = findmax(a,b)
82 if(a > b)
83 max0=a;
84 return;
85 else
86 max0=b;
87 return;
88 endif
89 endfunction
90
91 function min0 = findmin(a,b)
92 if(a < b)
93 min0=a;
94 return;
95 else
96 min0=b;
97 return;
98 endif
99 endfunction
100
101 function Case1
102 N=4;
103 D_x=[0, 3, -2, 1];
104 D_y=[0, -1, 3, 2];
105 P_x=[5, 2, -3, 3];
106 P_y=[3, 1, -1, 2];
107 radius=[1, 2, 3, 4];
108 for i=[1:1:N]
109 Dx=D_x(i);
110 Dy=D_y(i);
111 Px=P_x(i);
112 Py=P_y(i);
113 r=radius(i);
114 [Ex,Ey]=CalculateCirclePoint(Dx,Dy,Px,Py,r)
115 filename=["plot" int2str(i) ".png"]
116 Case1Plot(Dx,Dy,Px,Py,Ex,Ey,r,N,i,filename);
117 endfor
118 endfunction
119
120 function Case2
121 C=[3, -4, -4, 5];
122 D=[3, 3, -3, -4];
123 A=[2, 1, 3, 2];
124 B=[1, 2, 1, 3];
125 N=4;
126 for m=[N+1:1:N+4]
127 c=C(m-N)
128 d=D(m-N)
129 a=A(m-N)
130 b=B(m-N)
131 x=zeros(1,20);
132 xin = CalculateInitialGuess(c,d,a,b)
133 x(1)=xin;
134 for n=[1:1:20]
135 x0=x(n);
136 fterm1num=b^2*d*x0;
137 fterm1den=(b^2-a^2)*x0+a^2*c;
138 fterm1=fterm1num/fterm1den;
139 if (d<0)
140 fterm2=b*sqrt(1-x0^2/a^2);
141 else
142 fterm2=-b*sqrt(1-x0^2/a^2);
143 endif
144 f=fterm1+fterm2;
145 fpterm1num=b^2*d*((b^2-a^2)*x0+a^2*c)-b^2*d*x0*(b^2-a^2);
146 fpterm1den=((b^2-a^2)*x0+a^2*c)^2;
147 fpterm1=fpterm1num/fpterm1den;
148 if (d<0)
149 fpterm2=-(2*b*x0)/(a^2*sqrt(1-x0^2/a^2));
150 else
151 fpterm2=(2*b*x0)/(a^2*sqrt(1-x0^2/a^2));
152 endif
153 fp=fpterm1+fpterm2;
154 x(n+1)=x(n)-f/fp;
155 f
156 endfor
157 x
158 x20=real(x(20))
159 if (d<0)
160 y20=-b*sqrt(1-x20^2/a^2)
161 else
162 y20=b*sqrt(1-x20^2/a^2)
163 endif
164 Ex=x20;
165 Ey=y20;
166 filename=["ellipseplot" int2str(m) ".png"]
167 Case2Plot(0,0,a,b,c,d,Ex,Ey,m,filename)
168 minus1_ProductOfSlopes_check=-(b^2*Ex*(Ey-d))/(a^2*Ey*(Ex-c))
169 endfor
170 endfunction
171
172 Case1
173 Case2