あなたの場合、「多角形の側面」がないため、に類似したものはなく、参照されるすべての場所と参照するすべての変数theta
を変更する必要があります( など)。theta
theta
phi
次のようなものが機能するはずです。リンクからコードをコピーして貼り付け、 と を削除し、theta
とphi
を再定義cdist
しpolyEq
、多角形の代わりに答えに単位円をプロットしました。これらの選択肢のいずれかが不明な場合は、質問してください。
M = 19; % number of circles to fit
toms r % radius of circles
x = tom('x', 2, M); % coordinates of circle centers
clear pi % 3.1415...
cdist = 1; % radius of unit circle
%%%%%% equations saying all circles are inside of unit circle
polyEq = (sqrt(x(1,:).^2 + x(2,:).^2) + r <= cdist);
% create a set of equations that say that no circles overlap
circEq = cell(M-1,1);
for i=1:M-1
circEq{i} = ( sqrt(sum((x(:,i+1:end)-repmat(x(:,i),1,M-i)).^2)) >= 2*r );
end
% starting guess
x0 = { r == 0.5*sqrt(1/M), x == 0.3*randn(size(x)) };
% solve the problem, maximizing r
options = struct;
% try multiple starting guesses and choose best result
options.solver = 'multimin';
options.xInit = 30; % number of different starting guesses
solution = ezsolve(-r,{polyEq,circEq},x0,options);
% plot result
x_vals = (0:100)./100;
plot(sqrt(1 - x_vals.^2),'-') % top of unit circle
plot(-1.*sqrt(1 - x_vals.^2),'-') % bottom of unit circle
axis image
hold on
alpha = linspace(0,2*pi,100);
cx = solution.r*cos(alpha);
cy = solution.r*sin(alpha);
for i=1:M
plot(solution.x(1,i)+cx,solution.x(2,i)+cy) % circle number i
end
hold off
title(['Maximum radius = ' num2str(solution.r)]);