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