The entire codebase in Matlab - Codebase for calculating the area and perimeter of a variable non-convex polygons and animates burn Working Progress - CEA Implimentation
clf
clear
Radius = input('Radius of Fueslage:');
t = linspace(0,2*pi,1000);
plot(Radius*cos(t),Radius*sin(t),'b')
axis equal
% point input
loop = 0;
while loop < 2
loop = input('How many points :');
end
disp('Pick points inside the circle to make a non - convex polgon')
temp = zeros(2,loop);
hold on
cord = ones(2,(19*loop));
for i = 1:loop
temp(:,i) =[ 1000000000; 1000000000];
while (temp(1,i)^2 + temp(2,i)^2 > Radius^2)
[temp(1,i),temp(2,i)] = ginput(1);
if (temp(1,i)^2 + temp(2,i)^2 > Radius^2)
disp('Point shold be inside circle')
end
end
plot(temp(1,i),temp(2,i),'ro')
text(temp(1,i),temp(2,i),int2str(i),'FontSize',14)
end
%ordering the points
if loop > 3
x_m = mean(temp(1,:));
y_m = mean(temp(2,:));
d = zeros(loop,loop);
for i = 1:loop
for j = 1:loop
if j ~= i
d(i,j) = mid_dist(temp(1,i),temp(2,i),temp(1,j),temp(2,j),x_m,y_m);
end
end
end
%creating point map based on farthest distance of line from centroid
pointmap = 1:1:loop;
i = 1;
[m ,pointmap(i)] = max(d(:,i));
d(i,:) = zeros(1,loop);
i = pointmap(i);
while i ~= 1
[m ,pointmap(i)] = max(d(:,i));
d(i,:) = zeros(1,loop);
i = pointmap(i);
end
temp1 = temp;
% creating lines from the points
i = pointmap(1);
for k = 2:1:loop
temp(1,k) = temp1(1,i);
temp(2,k) = temp1(2,i);
i = pointmap(i);
end
end
% creating lines from the points
for i = 1:loop
if i > 1
segment_x = linspace(temp(1,i-1),temp(1,i),20);
segment_y = linspace(temp(2,i-1),temp(2,i),20);
else
segment_x = linspace(temp(1,loop),temp(1,1),20);
segment_y = linspace(temp(2,loop),temp(2,1),20);
end
segment_x = segment_x(2:end);
segment_y = segment_y(2:end);
cord(1,((19*(i-1)+1):19*i)) = segment_x;
cord(2,((19*(i-1)+1):19*i)) = segment_y;
end
temp = [temp temp(:,1)];
hold on
% plot(temp(1,:),temp(2,:),'r')
plot(cord(1,:),cord(2,:),'go')
%Burn code
[t_b,O] = burn_time(cord,0.1,Radius);
plot(O(1),O(2),'k*');
r = 0.1 * ones(1,t_b);
[x,y] = burn(cord,r,Radius,t_b);
%plot cord
h = plot(x(1,:),y(1,:),'go');
axis equal
for j = 1:t_b
set(h,'XData',x(j,:));
set(h,'YData',y(j,:));
drawnow;
pause(0.075)
end
function[d] = mid_dist(x_1,y_1,x_2,y_2,x_0,y_0)
x_m = (x_1 + x_2)/2;
y_m = (y_1 + y_2)/2;
d = sqrt((x_0 - x_m)^2 + (y_0 - y_m)^2);
end
function[d] = perp_dist(x_1,y_1,x_2,y_2,x_0,y_0)
d = abs(((x_2 - x_1)*(y_1 - y_0) - (x_1 - x_0)*(y_2 - y_1))/(sqrt((x_2 - x_1)^2 + (y_2 - y_1)^2))) ;
end
function [P] = intersection(a_1,b_1,a_d,b_d,r)
syms y_1 x_1 x_d y_d m y x
if a_1 == a_d
P(1) = a_1;
temp_y = sqrt(r^2 - a_1^2 );
if temp_y > b_d
temp_y = -temp_y;
end
P(2) = temp_y;
else
m = (y_d - y_1)/(x_d - x_1);
eqn = (y_d - y)/(x_d - x) == m;
S = solve(eqn,y);
eq2 = x^2 + S^2 == r^2;
K = solve(eq2,x);
x_1 = a_1;
y_1 = b_1;
x_d = a_d;
y_d = b_d;
L = subs(K);
M = subs(S);
if L(1) > a_1 && L(1) < a_d
P(1) = L(1);
else
P(1) = L(2);
end
P(2) = subs(M,x,P(1));
end
end
function[ t,c ] = burn_time(cord,r_b,R)
x_cord = [cord(1,:)];
y_cord = [cord(2,:)];
x_m = mean(x_cord);
y_m = mean(y_cord);
x_cord_new = x_cord;
y_cord_new = y_cord;
k = size(x_cord);
n = (x_cord.^2 + y_cord.^2).^(0.5);
[M,index] = min(n);
i = index;
if (i > 1) && (i < k(2)-1)
tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(i+1) - x_cord(i)))/2....
, -((y_cord(i) - y_cord(i-1)) + (y_cord(i+1) - y_cord(i)))/2 ] ;
elseif i == 1
tangent_dir = [((x_cord(i) - x_cord(k(2))) + (x_cord(i+1) - x_cord(i)))/2....
, -((y_cord(i) - y_cord(k(2))) + (y_cord(i+1) - y_cord(i)))/2 ] ;
elseif i == k(2)
tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(1) - x_cord(i)))/2....
, -((y_cord(i) - y_cord(i-1)) + (y_cord(1) - y_cord(i)))/2 ] ;
end
temp_burn_dir = flip(tangent_dir) / norm(tangent_dir);
point = [x_cord(i) - x_m , y_cord(i) - y_m] ;
new_point = point + temp_burn_dir ;
if norm(new_point) < norm(point)
burn_dir = -temp_burn_dir;
else
burn_dir = temp_burn_dir;
end
x_cord_new(i) = r_b*burn_dir(1) + x_cord(i);
y_cord_new(i) = r_b*burn_dir(2) + y_cord(i);
t = 1;
while x_cord_new(i)^2 + y_cord_new(i)^2 < R^2
x_cord_new(i) = r_b*burn_dir(1) + x_cord_new(i);
y_cord_new(i) = r_b*burn_dir(2) + y_cord_new(i);
t = t+1;
end
t = t+10;
c = [x_cord(i) ; y_cord(i)];
end
function[x,y] = burn(cord,r_b,R,t_b)
x_cord = [cord(1,:)];
y_cord = [cord(2,:)];
plot(x_cord,y_cord,'ro')
x_cord_new = x_cord;
y_cord_new = y_cord;
k = size(x_cord);
x = ones(t_b,k(2));
y = ones(t_b,k(2));
x(1,:) = x_cord;
y(1,:) = y_cord;
end_flag = zeros(1,k(2));
for j = 1:t_b
for i = 1:k(2)
if end_flag(i) == 0
if (i > 1) && (i < k(2)-1)
tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(i+1) - x_cord(i)))/2....
, -((y_cord(i) - y_cord(i-1)) + (y_cord(i+1) - y_cord(i)))/2 ] ;
elseif i == 1
tangent_dir = [((x_cord(i) - x_cord(k(2))) + (x_cord(i+1) - x_cord(i)))/2....
, -((y_cord(i) - y_cord(k(2))) + (y_cord(i+1) - y_cord(i)))/2 ] ;
elseif i == k(2)
tangent_dir = [((x_cord(i) - x_cord(i-1)) + (x_cord(1) - x_cord(i)))/2....
, -((y_cord(i) - y_cord(i-1)) + (y_cord(1) - y_cord(i)))/2 ] ;
end
temp_burn_dir = r_b(j)*flip(tangent_dir) / norm(tangent_dir);
point = [x_cord(i) - mean(x_cord), y_cord(i)-mean(y_cord)] ;
new_point = point + temp_burn_dir ;
if norm(new_point) < norm(point)
burn_dir = -temp_burn_dir;
else
burn_dir = temp_burn_dir;
end
x_cord_new(i) = burn_dir(1)+x_cord(i);
y_cord_new(i) = burn_dir(2)+y_cord(i);
if sqrt(x_cord_new(i)^2 + y_cord_new(i)^2) > R
K = intersection(x_cord(i),y_cord(i),x_cord_new(i),y_cord_new(i),R);
x_cord_new(i) = K(1);
y_cord_new(i) = K(2);
end_flag(i) = 1;
end
else
x_cord_new(i) = x_cord(i);
y_cord_new(i) = y_cord(i);
end
end
x_cord = x_cord_new;
y_cord = y_cord_new;
x(j,:) = x_cord_new;
y(j,:) = y_cord_new;
end
end